章 鑫, 劉世杰, 李彬彬, 童小華

(1.同濟(jì)大學(xué) 測繪與地理信息學(xué)院，上海 200092; 2.同濟(jì)大學(xué) 土木工程防災(zāi)國家重點(diǎn)實(shí)驗(yàn)室，上海 200092)

合成孔徑雷達(dá)干涉測量(Interferometric Synthetic Aperture Radar，InSAR)是獲取數(shù)字高程模型(Digital Elevation Model，DEM)數(shù)據(jù)的重要手段，能夠描述地面的起伏形態(tài)特征，廣泛地應(yīng)用于地形分析、地質(zhì)、水文監(jiān)測、生態(tài)、農(nóng)業(yè)以及軍事應(yīng)用等領(lǐng)域[1-2]，但其精度受到軌道定位、SAR影像的配準(zhǔn)、相位解纏等精度的影響[3-4]。為了提高DEM的利用價值，需要對DEM進(jìn)行精度提升。

大量學(xué)者對于InSAR 得到的DEM產(chǎn)品進(jìn)行了質(zhì)量提升研究。Crosetto[5]利用控制點(diǎn)及升降軌數(shù)據(jù)重疊區(qū)的公共點(diǎn)對ERS-1(Earth Remote Sensing Satellite-1)數(shù)據(jù)進(jìn)行優(yōu)化，分析并提出了高程誤差模型，DEM精度由18.1 m提升至11.4 m，但是沒提及DEM平面位置精度。Wessel等[6]總結(jié)了InSAR生成的DEM的校準(zhǔn)方法，給出了TanDEM-X任務(wù)的DEM誤差校正方案，在模擬與實(shí)際DEM中進(jìn)行了DEM的區(qū)域網(wǎng)平差處理，也只實(shí)現(xiàn)了高程精度的提升。Gruber等[7-8]以激光測高(Ice、Cloud and Land Elevation Satellite，ICESat)數(shù)據(jù)的高程作為地面控制點(diǎn)，使用區(qū)域網(wǎng)平差的方式對TanDEM-X DEM進(jìn)行系統(tǒng)誤差校正，結(jié)果證明了通過DEM區(qū)域網(wǎng)平差可顯著提高相鄰DEM間的相對高程精度，但是沒考慮DEM初始定位誤差對ICESat約束的影響。和會等[9]利用ICESat 對南極ERS數(shù)據(jù)的 InSAR DEM做了糾正，但采用的是簡單的線性高程糾正模型，需要密集的控制點(diǎn)分布，且沒有利用DEM重疊區(qū)的一致性約束。Hueso等[3]結(jié)合E-SAR系統(tǒng)對TanDEM-X的校正方法進(jìn)行實(shí)驗(yàn)驗(yàn)證，并給出TanDEM-X的校正方法過程，分析了DEM系統(tǒng)誤差的成分特性。Hossein等[10]用ICESat數(shù)據(jù)作為控制點(diǎn)，利用移動平均插值算法得到DEM間的差值面來提高先進(jìn)星載熱發(fā)射和反射輻射儀傳感器生成的全球DEM的高程精度，只利用了控制數(shù)據(jù)來進(jìn)行高程上的誤差擬合。張永俊[11]利用了多源控制點(diǎn)基于定標(biāo)的方法逐景對星載分布式InSAR DEM做了矯正。已有研究采用了簡單線性高程糾正模型，依賴密集的控制點(diǎn)，且沒有考慮重疊區(qū)的約束。而在引入?yún)^(qū)域網(wǎng)平差方法時，直接將控制點(diǎn)的平面位置綁定到DEM對應(yīng)位置的柵格進(jìn)行高程校正，沒有考慮DEM初始平面誤差造成的影響，這會使控制點(diǎn)產(chǎn)生錯誤的高程約束。

ICESat-2[12-13]是2018年發(fā)射的搭載了光子計數(shù)激光測高儀的傳感器，相比2003年發(fā)射的ICESat所搭載的全波形激光雷達(dá)系統(tǒng),具有更高的精度和更小的激光足印。筆者同時考慮InSAR DEM的平面和高程改正，基于星載InSAR DEM誤差特性分析，構(gòu)建InSAR DEM區(qū)域網(wǎng)平差方法，利用ICESat-2激光數(shù)據(jù)產(chǎn)品進(jìn)行高程約束，對DEM的平面、高程精度進(jìn)行校正，輸出大范圍的高精度的DEM產(chǎn)品。并對比了有無平面控制情況下，區(qū)域網(wǎng)平差優(yōu)化的DEM結(jié)果。

1 算法基礎(chǔ)

DEM平差流程如圖1所示。首先收集Sentinel A IW影像數(shù)據(jù)對，進(jìn)行干涉處理生成對應(yīng)的DEM數(shù)據(jù)；再對收集的ICESat-2激光點(diǎn)篩選生成區(qū)域網(wǎng)平差所需的高程控制點(diǎn)和檢核點(diǎn)；最后針對生成的DEM，先后進(jìn)行連接點(diǎn)匹配、區(qū)域網(wǎng)平差、DEM糾正處理過程。本文將區(qū)域網(wǎng)平差分解為平面定位六參數(shù)的優(yōu)化和高程誤差模型的平差求解兩個獨(dú)立的平差過程，聯(lián)同連接點(diǎn)的地理坐標(biāo)約束和高程約束來提高DEM的平面定位和高程精度。

圖1 DEM平差流程

1.1 激光點(diǎn)篩選

針對目標(biāo)DEM范圍內(nèi)的ICESat-2點(diǎn)，首先根據(jù)ATL08產(chǎn)品的屬性標(biāo)簽進(jìn)行篩選，基于Li等[14]提出的篩選流程來提取初始的ATL08高精度數(shù)據(jù)，提取的ATL08點(diǎn)分成兩部分，分別用于控制和檢核。此外，本文引入兩個新的約束：① 由于DEM可能存在初始平面定位誤差，保留的高程控制數(shù)據(jù)應(yīng)該位于平坦的地形上，且激光點(diǎn)對應(yīng)的DEM柵格應(yīng)該處于平坦區(qū)域[15]，如式(1)所示，通過設(shè)定閾值對DEM柵格周圍像素高程的方差進(jìn)行限制。② 過于密集的點(diǎn)對平差結(jié)果沒有明顯改善，因此本文通過格網(wǎng)均勻地選取高精度的激光高程控制數(shù)據(jù)。對于部分格網(wǎng)內(nèi)沒有控制點(diǎn)分布的情況，本文引入虛擬控制點(diǎn)作為輔助約束。

(1)

式中：R和C為目標(biāo)柵格周邊鄰域的行列號；μ為目標(biāo)柵格的高程值。由此來統(tǒng)計其鄰域柵格的高程差異變化。

1.2 DEM 連接點(diǎn)匹配

DEM在重疊區(qū)域存在高程誤差和定位誤差，提取連接點(diǎn)需綜合考慮兩者誤差的影響。DEM特征包括高程、坡度和紋理等，其中高程信息最為直接地顯示DEM地形信息。利用滑動匹配過程對DEM重疊區(qū)公共連接點(diǎn)進(jìn)行初始提取，對于窗口內(nèi)獲得的整像素，進(jìn)一步利用雙線性插值方法，使定位誤差校正于亞像素級，其中匹配測度函數(shù)采用比較穩(wěn)健且效率較高的高程差標(biāo)準(zhǔn)差[3]，實(shí)現(xiàn)連接點(diǎn)高精度提取。

1.3 區(qū)域網(wǎng)平差

DEM平面平差是利用連接點(diǎn)的平面約束進(jìn)行平面區(qū)域網(wǎng)平差。對于DEM，每個柵格的地理坐標(biāo)由角點(diǎn)的大地坐標(biāo)以及像元分辨率確定。平差前后的DEM的分辨率不會改變，因此對DEM整體做平面剛體變換糾正，如式(2)所示。

(2)

式中：X和Y分別為DEM某個柵格的地理坐標(biāo)；x和y為對應(yīng)的行列號；dx和dy為DEM的角點(diǎn)經(jīng)緯度，即平面區(qū)域網(wǎng)平差過程將優(yōu)化DEM的定位六參數(shù)中角點(diǎn)的平面坐標(biāo);矩陣A為旋轉(zhuǎn)參數(shù)和比例參數(shù)(分辨率)形成的矩陣。此時，通過引入DEM連接點(diǎn)的平面坐標(biāo)一致性約束，建立約束方程為

(3)

式中:X和Y分別為地理坐標(biāo)，包含同名點(diǎn)的DEMJ和K對應(yīng)的地理坐標(biāo)在優(yōu)化后應(yīng)該達(dá)成一致。通過至少兩對連接點(diǎn)建立的約束方程，即可通過最小二乘求解對應(yīng)的DEM平面定位參數(shù)糾正模型的系數(shù)。

DEM高程區(qū)域網(wǎng)平差是為了校正不同數(shù)據(jù)集之間的相對誤差，各DEM的高程誤差改正模型參數(shù)通過高程控制點(diǎn)進(jìn)行約束優(yōu)化，同時利用重疊區(qū)的連接點(diǎn)高程值一致條件建立約束。通過Gruber等[7]的研究，針對InSAR生成 DEM，其誤差分布模型符合：

(4)

式中：x為DEM的列號；y為行號。即由于SAR的成像模式，導(dǎo)致其方位向和距離向的誤差分布規(guī)律不一致，在任意方位時刻干涉的各項參數(shù)誤差恒定，因此距離方向高程糾正采用一階模型，而方位向的階數(shù)由于衛(wèi)星的位置變化通常采用多項式模型[6]。對于DEMJ和K上某個確定的連接點(diǎn)對(xJ,yJ)和(xK,yK)，平差約束則為它們的原始高程加上對應(yīng)的高程改正項達(dá)成一致，如式(5)所示。

hDEM_J(xJ,yJ)-gJ(xJ,yJ)=hDEM_K(xK,yK)-gK(xK,yK)

(5)

式中：gJ(xJ,yJ)和gK(xK,yK)分別為DEMJ和K的高程改正模型；hDEM_J(xJ,yJ)為DEMJ在(xJ,yJ)處的初始高程；hDEM_K(xK,yK)為DEMK在(xK,yK)處的初始高程。結(jié)合對控制點(diǎn)的約束，觀測方程可寫為

V=AX-L

(6)

式中：X為高程校正模型的系數(shù)向量。矩陣A為

(7)

式中：A的前兩行表示 ICESat-2 測高數(shù)據(jù)，分別用于改進(jìn)DEMJ和K。第3行表示相鄰 DEM 連接點(diǎn)之間的高程限制。然后利用最小二乘原理消除各DEM的高程系統(tǒng)誤差，得到誤差校正模型系數(shù)。針對高程控制點(diǎn)缺失的部分，為了解決缺少控制點(diǎn)部分導(dǎo)致的局部精度低的問題且使平差穩(wěn)定，引入高程虛擬控制點(diǎn)，降低對應(yīng)的權(quán)重，來達(dá)到虛擬約束。

2 實(shí)驗(yàn)結(jié)果

2.1 實(shí)驗(yàn)數(shù)據(jù)

本文收集了9對間隔12天的Sentinel A IW影像數(shù)據(jù)，通過干涉處理生成了對應(yīng)的DEM數(shù)據(jù)，如圖2所示。初始的DEM由于初始系統(tǒng)和干涉過程引入的相位誤差，導(dǎo)致DEM重疊部分存在明顯的紋理不一致現(xiàn)象。

2.2 實(shí)驗(yàn)結(jié)果

首先對收集的ICESat-2激光點(diǎn)進(jìn)行篩選，篩選前后DEM范圍內(nèi)的激光點(diǎn)分布對比如圖3所示。篩選前6866個點(diǎn)，篩選后設(shè)置200柵格大小的格網(wǎng)約束?？紤]到平坦地區(qū)的坡度要求為小于2°，且DEM分辨率為20 m，本文設(shè)置DEM柵格的鄰域方差閾值為(20×tan2°)2。由此保留了1230個激光點(diǎn)，其中按照三分之一的比例從每景DEM上分布的激光點(diǎn)中隨機(jī)選擇檢核點(diǎn)，剩余的作為控制點(diǎn)。隨后匹配出508對同名連接點(diǎn)。平面和高程區(qū)域網(wǎng)平差依次進(jìn)行后，得到了平差糾正后的DEM數(shù)據(jù)，平差前后的DEM對比如圖4所示。DEM重疊區(qū)的高程不一致現(xiàn)象被明顯糾正，通過連接點(diǎn)的高程一致性約束，平差的DEM在重疊區(qū)高程達(dá)到了很好的吻合。同時，為了對比本文方法和傳統(tǒng)不考慮平面的平差方法，利用同一數(shù)據(jù)只進(jìn)行了高程區(qū)域網(wǎng)平差，平差后DEM不一致性也能被良好地糾正。但在不考慮平面誤差時，由于初始的平面誤差，激光檢核點(diǎn)分布的位置也會存在一定的錯位。圖5為DEM平差后局部一致性對比，選取兩個DEM重疊區(qū)紋理進(jìn)行了放大。圖5(a)中初始的DEM平面和高程不一致性十分明顯。不考慮平面偏移時，經(jīng)過高程平差后的圖5(b)結(jié)果表明高程一致性得到一定提升，但是平面上仍存在著初始偏移。而由圖5(c)可知，考慮平面偏移的平差后DEM重疊區(qū)的高程和平面紋理的一致性都十分吻合。區(qū)域二的初始高程誤差更為明顯，整體對比也更為直觀。圖6為DEM平差前后檢核點(diǎn)高程誤差直方圖對比，不考慮平面誤差時，高程檢核點(diǎn)精度也有提升，誤差趨于0均值的正態(tài)分布，但是殘差接近于0的部分相較本文方法相對較少。表1為DEM平差結(jié)果的精度評估，統(tǒng)計了檢核點(diǎn)誤差和連接點(diǎn)定位坐標(biāo)的誤差RMSE。經(jīng)過平面平差，連接點(diǎn)的平面定位相對精度也有明顯提升，從初始的182.462 m提升至14.887 m，優(yōu)于一個像素。而只進(jìn)行高程平差時，DEM的高程誤差RMSE從平差前的19.372 m提升到了7.865 m，這是由于檢核點(diǎn)點(diǎn)位分布仍存在錯誤的定位，經(jīng)過本文提出的考慮平面平差的一體化平差方法，最后的高程精度提升至3.459 m。

圖2 初始9景Sentinel DEM數(shù)據(jù)

圖3 區(qū)域內(nèi)所有ICESat-2數(shù)據(jù)

圖4 平差前后DEM局部提升效果對比

表1 DEM平差結(jié)果的精度評估

圖5 DEM平差后局部一致性對比

圖6 DEM平差前后檢核點(diǎn)高程誤差直方圖對比

4 結(jié)束語

提出了一種針對InSAR DEM數(shù)據(jù)的區(qū)域網(wǎng)平差糾正方法，同時考慮了DEM的平面位置偏移和高程誤差，并通過ICESat-2激光測高數(shù)據(jù)進(jìn)行高程約束和檢核。對哨兵1A數(shù)據(jù)干涉生成的DEM數(shù)據(jù)進(jìn)行了驗(yàn)證，通過和傳統(tǒng)的不考慮平面誤差的DEM糾正方法進(jìn)行對比，證明本方法顯著地提升了DEM的高程精度(提升至82%以上)以及平面的紋理一致性(提升至92%)。后續(xù)考慮將試驗(yàn)推廣到更多傳感器生成的DEM數(shù)據(jù)中，研究多傳感器數(shù)據(jù)間的區(qū)域網(wǎng)平差方法。