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

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

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

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

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

1 算法基礎(chǔ)

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

圖1 DEM平差流程

1.1 激光點(diǎn)篩選

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

(1)

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

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

DEM在重疊區(qū)域存在高程誤差和定位誤差，提取連接點(diǎn)需綜合考慮兩者誤差的影響。DEM特征包括高程、坡度和紋理等，其中高程信息最為直接地顯示DEM地形信息。利用滑動(dòng)匹配過(guò)程對(duì)DEM重疊區(qū)公共連接點(diǎn)進(jìn)行初始提取，對(duì)于窗口內(nèi)獲得的整像素，進(jìn)一步利用雙線性插值方法，使定位誤差校正于亞像素級(jí)，其中匹配測(cè)度函數(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)平差。對(duì)于DEM，每個(gè)柵格的地理坐標(biāo)由角點(diǎn)的大地坐標(biāo)以及像元分辨率確定。平差前后的DEM的分辨率不會(huì)改變，因此對(duì)DEM整體做平面剛體變換糾正，如式(2)所示。

(2)

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

(3)

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

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

(4)

式中：x為DEM的列號(hào)；y為行號(hào)。即由于SAR的成像模式，導(dǎo)致其方位向和距離向的誤差分布規(guī)律不一致，在任意方位時(shí)刻干涉的各項(xiàng)參數(shù)誤差恒定，因此距離方向高程糾正采用一階模型，而方位向的階數(shù)由于衛(wèi)星的位置變化通常采用多項(xiàng)式模型[6]。對(duì)于DEMJ和K上某個(gè)確定的連接點(diǎn)對(duì)(xJ,yJ)和(xK,yK)，平差約束則為它們的原始高程加上對(duì)應(yīng)的高程改正項(xià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é)合對(duì)控制點(diǎn)的約束，觀測(cè)方程可寫(xiě)為

V=AX-L

(6)

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

(7)

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

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

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

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

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

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

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

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

圖4 平差前后DEM局部提升效果對(duì)比

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

圖5 DEM平差后局部一致性對(duì)比

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

4 結(jié)束語(yǔ)

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