利用相干系數改進的數字高程模型反演方法

樂 穎，夏元平，錢文龍

(1.東華理工大學測繪工程學院，江西 南昌 330013；2.福州市勘測院，福建 福州 350108)

0 引言

數字高程模型(digital elevation model，DEM)蘊含了豐富的地表高程信息，在測繪、地質地貌、工程建設、環境治理等諸多領域具有極其重要的作用[1-4]。由于我國南北方氣候差異明顯，尤其是在我國西南部區域，經常伴有多云多雨天氣，通過攝影測量、光學遙感等常規手段很難獲取這些區域的DEM[5]。此外，相較于光學遙感，合成孔徑雷達干涉測量(interferometric synthetic aperture radar，InSAR)反演DEM時，不受天氣狀況的限制，可全天候全天時進行數據采集，并且具有高精度、高分辨率的優勢，從而成為了獲取 DEM 的常用手段之一[6]。針對不同地區、不同應用領域，分析其影響因素，并嘗試各種改進方法獲取DEM，從而改善DEM精度仍是目前研究的熱點[7]。

InSAR技術是目前獲取DEM的常用手段，通過獲取同一區域兩幅SAR影像，并利用影像上的相位信息進行差分干涉，從而反演出DEM。然而，InSAR技術由于軌道誤差、大氣影響、空間與時間失相關等固有缺陷，在實際應用中受到一定的影響。一些學者基于InSAR技術進行了改進研究，從而提高 InSAR 數據處理的效率和精度，并將其有效地用來反演DEM[8]。文獻[9]以DEM數據為基礎，結合強局部加權回歸算法，利用反距離加權方法進行空間插值，從而獲得了可靠性更好的河床重構數據。文獻[10]提出了一種基于DEM輔助后向投影模型的InSAR高程反演方法，顯著降低了干涉條紋的密度，有效地提高了DEM反演精度。文獻[11]采用InSAR技術獲取礦區DEM，并對定南縣嶺北離子型稀土礦區數據進行實驗，得到了該礦區的地形地貌特征。但是，上述文獻在InSAR反演DEM過程中，主要依靠人工憑借其經驗選取軌道精煉控制點，這對操作者的專業性要求比較嚴格，否則難以選取誤差較小的控制點。鑒于此，本文提出利用相干系數改進的數字高程模型反演方法。

1 傳統InSAR技術反演DEM原理

SAR影像經由信號采集與數據處理之后，影像中每一像素都將蘊含與斜距相關的相位信息以及地面分辨元的雷達后向散射強度信息。InSAR技術主要從干涉相位與干涉相關數據中獲取感興趣的可用信息。而干涉相位一般與傳感器到目標的距離有直接關系，是InSAR數據處理與信號提取的關鍵部分。在實際應用中，InSAR干涉相位ψ主要包括參考橢球面相位φref、地形相位φtop、形變相位φdef、大氣相位φatm和噪聲相位φnoi等多個分量。其中φref、φtop和φdef為主要組成部分，通常采用較短的時間基線和較長的空間基線干涉對；φnoi可采用低通相位濾波的方法有效抑制，φatm可根據相位的高通濾波進行一定的削弱。因此，干涉相位ψ可表示為[12]：

ψ=φref+φtop+φdef+φatm+φnoi。

(1)

φref是由參考橢球面本身引起的干涉相位，計算公式為：

(2)

式(2)中，λ表示波長，r1和r2分別為影像中任一像元等斜距投影到參考橢球面后對應于主、副影像的斜距，R1為任一像元對應的主影像斜距。

由于B∥=δR=r1-r2=Bsin(θ0-α)，將其代入式(2)，φref則可表示為：

(3)

式(3)中，θ0為主影像衛星雷達至地面點目標的側視角，B和α分別是經由衛星軌道數據和配準結果計算所得到的對應基線長度和基線傾角。

φtop是由目標大地高h所引起的相位分量，計算公式為：

(4)

由于斜距r和衛星高度H可遠大于觀測目標大地高h，即由目標高程引起的側視角變化量δθ將很小，因此，δθ可表示為：

(5)

將式(5)代入式(4)中，得到φtop分量的表達式為：

(6)

式(6)中，B⊥為干涉垂直基線長度。

φdef計算的是沿雷達視線方向的一維地表形變量，計算公式為：

(7)

式(7)中，Δr是由于觀測目標位移造成的沿雷達視線方向的斜距變化量。

最后，通過獲取同一區域的兩景SAR影像形成一個干涉對，而對于地面上的每一分辨元來說，ψ直接對應著地面分辨元與兩傳感器間的距離之差的精確信息。因此，可利用傳感器位置數據、SAR系統參數和基線參數等，借助于ψ計算每一地面分辨元的三維位置，即為DEM的反演。

2 利用相干系數改進的DEM反演方法

2.1 利用相干系數的DEM反演方法原理

將同一區域的主、副兩幅SAR影像所對應的像素相位值進行差分，便可得獲得一個一次差分相位圖，即干涉相位圖。在傳統InSAR技術反演DEM的數據處理過程中，主要包括以下幾個步驟：影像配準、干涉圖生成、自適應濾波和相干性計算、最小費用流相位解纏、控制點選取、軌道精煉和重去平、相位轉形變及地理編碼等，處理流程圖如圖1所示。

圖1 實驗流程圖Fig.1 Experimental flow chart

上述每個數據處理步驟過程中產生的誤差都將對干涉結果的精度和可靠性產生巨大的影響，而利用相干系數改進的數字高程模型反演方法的具體操作步驟為：

1) 影像配準。利用軌道數據和參考DEM獲取局部非參數偏移估計，并計算每個窗口的互相關函數，取最大值作為所選位置的參考位移；然后依據方位角和距離像元位置信息，經由多項式計算殘差參數偏移，將其與原始局部非參數偏移估計相加，從而獲得細化的殘差參數偏移；最后，將局部非參數偏移和改進的基于互相關擬合進行配準，得到亞像素級精度。

2) 干涉圖生成。根據兩景配準后的SAR影像之間的相位差(Φ)，可計算得到地球上某地面目標與傳感器位置之間的距離差。干涉相位的表達式

(8)

式(8)中，Image(I)和Real(I)分別表示干涉圖的虛部和實部。

3) 自適應濾波和相干性計算。每個穩定區域中所包含的像素選擇，主要來源于該穩定區域中所有像素的均值(Mall)與新像素值(Mnew)之間的差值，該差值對Mall進行了歸一化處理，是可供選擇的候選值。采用區域增長方法識別新的候選像素，相似均值因子是由一個數字在線性尺度上的表示，公式為：

(9)

一般情況下，需要多次過程迭代，才可獲取最佳相似平均因子。

4) 最小費用流相位解纏。SAR影像中所有像素所對應的相位皆有整周模糊度問題，干涉圖的相位只能是2π的倍數。在干涉處理中，本文采用最小費用流(minimum cost flow，MCF)方法為所有像素確定干涉相位的整周未知數。此類方法采用正方形的格網，綜合了圖像上所有的像元，并將相干性小于閾值的像元進行掩膜處理。將解纏最小相干性閾值設置為 0.15，若像元的相干系數大于0.15，則將該像元進行解纏。相較于區域增長法，當存在大范圍的低相干或其他限制增長的因素導致解纏困難時，最小費用流算法可以取得更優的結果。

5) 控制點選取。在缺少地面控制點信息并對研究區形變沒有先驗認識的情況，盲目選擇地面控制點將大概率引入誤差，從而導致最終結果產生偏差。因此，為了減少軌道誤差對DEM反演的影響，本文通過借助相干系數來剔除相干性低的像元點，保留相干性高的像元點作為控制點進行軌道精煉[13]。相干圖可描繪兩幅影像的相關程度，通過計算干涉條紋圖的相干系數，能夠定量評價SAR影像精確配準后所得到干涉條紋圖的質量。相干系數于1993年最初由普拉蒂提出，相干系數的取值范圍為[0，1]，值越接近1說明相干性越高，其理論模型公式為[14]：

(10)

式(10)中，E[·]表示數學期望，u*表示共軛復數，u1與u2為主副圖像的信號。在實際的InSAR處理過程中，在同一像素中無法獲取計算所需數量的采樣值。因此，根據式(10)的理論模型，基于SAR影像復數數據計算獲得相干系數的標準表達式為[15]：

(11)

式(11)中，u1(n，m)、u2(n，m)分別為主、副影像數據塊內某個坐標(n，m)處的復數值；|·|2為對應的二階范數；M與N分別為計算相干性的數據塊尺寸大?。籱與n為數據塊內對應的行列號。

6) 軌道精煉和重去平。該步驟主要是通過輸入控制點對軌道形態進行估算，而對于未解纏的相位信息，是否能夠正確轉換為高度(或位移)值也取決于該步驟是否成功。它可用于改進軌道(即糾正可能的誤差)，也可用于計算相位偏移(即獲得絕對相位值)，以及消除可能存在的相位偏移。

7)控制點選擇精度評定。為了更系統地驗證本文方法選取控制點的可靠性，首先統計干涉相干性系數圖中符合要求的像元數量，然后根據相干系數閾值選取軌道精煉控制點，最后計算各點絕對殘余誤差，計算公式為：

(12)

式(12)中，γ是干涉相干性測度，H是高程模糊度。

8) 相位轉形變及地理編碼。絕對校準和解纏相位與合成相位重新組合，并將其轉換為高度并編碼成地圖投影。通過對比距離-多普勒方法、相關的大地測量與制圖轉換方法，發現距離-多普勒方程可同時應用于兩個天線，不僅能夠提供各像素的高度，而且可得到其在給定的地圖和大地參考系統中的具體位置。

2.2 評價指標

本文以標準差和均方根偏差作為精度評價指標，對DEM的精度進行定量評定[16]。標準差(standard deviation，SD)是方差的算術平方根，用σ表示，在概率統計中最常作為統計分布程度上的測量依據。標準差主要反映數據的離散程度，值越小，則說明數據越穩定。其計算公式如下：

(13)

均方根偏差(root mean square error，RMSE)是一種常用的測量數值之間差異的量度，主要體現觀測值與真值之間的誤差，值越小，則說明方法所得結果越好。其計算公式如下：

(14)

3 實例驗證

3.1 研究區域概況

贛州市地處江西省南部，占據著省內最大的地域面積，總面積大約為整個江西省總面積的1/4，并且在江西省內擁有最多的區縣數量。定南縣作為贛州的區縣之一，總面積達到了1 321平方公里，坐落于江西省最南部。定南縣離子型稀土資源儲量非常豐富，涵蓋了20多種礦藏，并且縣區內河流分布廣泛，主要分屬于定南水和桃江兩水系[17]。

基于此，本文將定南縣嶺北鎮作為研究區域，地理坐標位置為東經114°47′49″E～115°22′48″E，北緯 24°33′37″N～25°03′21″N。圖2是贛州市定南縣的土壤類型分布圖，數據來源于世界土壤數據庫[18]。盡管定南縣的稀土種類豐富、品位高，但受限于落后的開采技術，單一的監測技術，亂采亂挖、釆富棄貧的現象仍然大規模的存在。長期以往導致大量礦產資源浪費，礦山地質環境愈發嚴重。因此，利用InSAR技術對DEM進行反演，從而獲取最新的數字高程模型，對于地表形變監測具有十分重要的意義。

圖2 贛州市定南縣土壤類型分布Fig.2 Distribution of soil types in Dingnan County, Ganzhou City

3.2 數據源

3.2.1雷達影像數據

本文使用的雷達影像數據為C波段的Sentinel-1A衛星的SLC數據，選取了2020年兩景影像數據，雷達影像數據的基本參數如表1所示。哨兵1號(Sentinel-1)衛星是歐洲航天局哥白尼計劃中的地球觀測衛星，由兩顆衛星組成，載有C波段合成孔徑雷達，可提供連續圖像。軌道信息是InSAR數據處理中非常重要的信息，從最初的圖像配準到最后的形變圖像生成都有著重要的作用。含有誤差的軌道信息造成基線誤差以殘差條紋的形式存在于干涉圖中。因此，本文使用了哨兵1號相應的衛星精密軌道數據對軌道信息進行修正，可去除因軌道誤差引起的系統性誤差。

表1 雷達影像數據的基本參數Tab.1 Basic parameters of radar image data

3.2.2外部高程數據

在雷達影像數據處理過程中需要借助外部DEM，用來提供參考地形或者參考地理坐標系。本文選取的是SRTM 30 m數據，SRTM(shuttle radar topography mission)即航天飛機雷達地形測繪使命，是美國太空總署(NASA)和國防部國家測繪局(NIMA)以及德國與意大利航天機構共同合作完成聯合測量，由美國發射的“奮進”號航天飛機上搭載SRTM系統完成[19]。SRTM系統共有約9.8萬億字節雷達影像的數據量，經過數據處理，制成了數字地形高程模型；于2003年開始公開發布，目前最新的版本為V4.1版本；分為30 m分辨率的SRTM1和90 m分辨率的SRTM3[20]。

3.3 實驗結果分析

3.3.1DEM反演結果

本文根據相干系數閾值選取了35個軌道精煉控制點，然后分別計算傳統InSAR技術反演DEM方法與本文方法選取軌道精煉控制點的絕對殘余誤差，結果對比如圖3所示。從圖3中可以明顯看出，本文方法用于選取軌道精煉控制點整體的絕對殘余誤差較小，只有個別點的絕對殘余誤差偏大，誤差符合要求；而傳統InSAR技術反演DEM方法選取控制點的絕對殘余誤差整體偏大，后續需要篩選絕對殘余誤差較大的點。此外，文中為了對比分析InSAR反演DEM方法與本文方法用于選取軌道精煉控制點的準確性，引入了標準差指標，并計算得到兩種方法絕對殘余誤差的標準差分別為0.71和0.34，足以定量說明本文方法準確性更高。因此，在缺少地面控制點信息或研究區形變先驗認識的情況下，若需要選擇可靠的地面控制點，可借助相干性系數進行剔除低相干性點。最終，通過軌道精煉和重去平步驟得到利用相干性系數的InSAR技術反演DEM與InSAR技術反演DEM這兩種方法的軌道精煉多項式分別為-4.67+0.000 017R+0.002 5A和-10.37-0.000 041R+0.000 88A，其中，R代表距離向，A代表方位向。

圖3 軌道精煉控制點精度對比Fig.3 Precision comparison of track refining control points

利用InSAR技術獲取的DEM分辨率為15 m，數據整體表現優良，對地形的細節刻畫明顯優于SRTM DEM[5]。圖4(a)、(b)和(c)分別為SRTM DEM、InSAR反演的DEM和本文方法反演的DEM；(d)、(e)和(f)分別是SRTM DEM山體陰影圖、InSAR反演DEM山體陰影圖和本文方法反演DEM山體陰影圖。由于時間跨度不同具有部分差異，但是這三種DEM結果在空間分布上整體一致，說明利用InSAR技術反演得到最新的DEM是可行的。在相同分辨率情況下，本文方法比InSAR反演的DEM方法結果在紋理信息和高程變化分布上與SRTM DEM更相符。

3.3.2結果分析

為了對比分析InSAR技術反演DEM結果和本文方法反演DEM結果與SRTM DEM之間的差異，分別從定性和定量兩方面對這兩種方法反演的DEM結果進行對比分析。

1) 定性分析

相干系數圖對DEM反演結果的質量有很大的影響，倘若整體的相干系數都處于一個很低的情況下，最終很難反演出準確的DEM結果。為了使得本文研究區域的相干性系數盡量高，本文選取了秋冬季節的兩景影像，避免了夏季受植被大范圍覆蓋的干擾，圖5為本文方法反演DEM的相干性系數圖。相干系數的取值范圍為在0和1之間，且值越接近1，代表著兩幅影像的相干性越高。從圖5中可以看出，整個研究區域整體的相干性都較高，從而說明，相干系數圖的質量較優，所得DEM反演結果準確性較高。與此同時，將本文方法反演DEM的等高線疊加在SRTM DEM上，結果見圖6。圖6表明這兩者的DEM結果高程起伏形態相同，且與SRTM DEM相比，本文方法反演DEM的分辨率更高，紋理信息更豐富，利用InSAR技術進行DEM反演結果是可靠的。利用相干系數的InSAR技術反演DEM可在減少軌道誤差的情況下，更準確且精細地反演出贛州市定南縣嶺北鎮最新的數字高程模型。

圖5 本文方法反演的DEM相干性系數圖Fig.5 DEM coherence coefficient diagram retrieved by this method

圖6 本文方法反演的DEM等高線疊加SRTM DEMFig.6 DEM contour line inversion by this method is superimposed with SRTM DEM

2) 擬合分析

為了方便分析兩組數據之間的內在聯系，常常通過數學模型進行曲線擬合。首先將InSAR技術反演DEM結果、本文方法反演DEM結果和SRTM DEM裁剪至同一區域范圍，并重采樣使這三種結果的行列數相同，最后利用線性模型、二次模型、指數模型、幾何模型、雙曲模型和對數平方模型這六種常見數學模型對三種DEM結果分別進行兩兩擬合，并比較選取最優擬合模型，其結果見圖7。本文方法反演DEM結果和SRTM DEM最優的擬合模型為對數平方模型，其表達式為Y=7 982.60-6 711.10×alog(x)+1 452.20×alog2(x)，擬合結果見圖7(a)；InSAR技術反演DEM結果和SRTM DEM最優的擬合模型為二次模型，其表達式為Y=238.92+0.35x+0.000 4x2，擬合結果見圖7(b)；本文方法反演DEM結果和InSAR技術反演DEM結果最優的擬合模型為二次模型，其表達式為Y=-159.20+1.12x-0.000 1x2，擬合結果見圖7(c)。結果表明，本文方法反演DEM結果比InSAR技術反演DEM結果與SRTM DEM更相似，相關性更高。

圖7 三種DEM間的擬合分析結果Fig.7 Fitting analysis results of three kinds of DEM

3.4 結果評估

利用InSAR技術對DEM進行反演可以獲取最新的數字高程模型，對于地表形變監測具有十分重要的意義。由于贛州市定南縣嶺北鎮的礦產資源豐富，因此本文結合已有礦權數據選取33個特征點代表InSAR技術反演稀土礦區DEM的精度對結果進一步定量評估，特征點的分布情況見圖8。以SRTM DEM作為參考DEM，分別對比InSAR技術反演DEM結果、本文方法反演DEM結果和參考DEM之間的標準差和均方根差，其結果對比見圖9。從圖中可以看出，這些特征點的高程趨勢高度一致，本文方法反演DEM結果比InSAR技術反演DEM結果更接近參考DEM的高程。SRTM DEM、InSAR技術反演DEM方法和本文方法的標準差分別為86.92、83.28和76.17，贛州市定南縣多為山區和居民地，因此地勢起伏差異會很大。InSAR技術反演DEM方法相較于SRTM DEM的均方根偏差和標準差分別是152.42 m和24.27 m；本文方法相較于SRTM DEM的均方根偏差和標準差分別是36.92 m和24.07 m。結果表明，相較于InSAR技術反演DEM方法，本文方法的均方根偏差和標準差分別降低了75.78 %和8.53%。由于InSAR反演DEM方法和本文方法所用數據和處理方法一致，僅僅是選取的控制點不同，因此，可以將這部分歸因于軌道精煉誤差導致的。利用相干系數改進的數字高程模型反演方法可以減少軌道精煉引起的誤差，使得最終反演的DEM結果更準確。

圖8 特征點選取分布圖Fig.8 Distribution diagram of feature points selection

圖9 SRTM、InSAR反演DEM和本文方法反演DEM結果對比Fig.9 Comparison of DEM inversion results from SRTM and InSAR and DEM inversion results from this method

4 結論

本文提出利用相干系數改進的數字高程模型反演方法。該方法通過設置一定的相干系數閾值，篩選出適宜的高相干性像元作為控制點進行軌道精煉處理，最后反演出DEM結果。以贛州市定南縣嶺北鎮作為研究區域，借助 2020 年兩景的Sentinel-1A衛星數據，以SRTM DEM作為參考DEM進行對比驗證，得到以下結論：1) 利用InSAR技術進行DEM反演是可行的，其結果在空間分布和紋理信息上與SRTM DEM結果一致；2) 利用相干系數改進的數字高程模型反演方法可以減少軌道精煉引起的誤差，其絕對殘余誤差的標準差為0.34，與InSAR技術反演DEM方法相比降低了52.11%，引入本文方法選取的控制點進行數據處理可以使最終反演出的DEM結果更準確；3) 以SRTM DEM為參考DEM，通過選取33個特征點對InSAR技術反演DEM方法和本文方法分別進行精度驗證，結果表明，相較于InSAR技術反演DEM方法，本文方法的均方根偏差和標準差分別降低了75.78%和8.53%。本實驗僅用免費的雷達數據進行了InSAR技術的反演DEM研究，對于其他雷達數據反演DEM的適用性還有待驗證。后續可收集其他雷達衛星的影像數據進行DEM反演研究，或者將利用相干系數反演的DEM作為參考DEM，用于DInSAR技術中進行地表形變監測研究。