一種基于高分辨率雷達影像以及外部DEM輔助的復雜地形制圖方法

盧麗君，張繼賢，王 騰

1.中國測繪科學研究院，北京100830；2.武漢大學測繪遙感信息工程國家重點實驗室，湖北武漢430079

1 引 言

干涉合成孔徑雷達（interferometric synthetic aperture radar，InSAR）作為一種主動的微波遙感技術在地形測繪已經顯現(xiàn)出越來越大的優(yōu)勢［1－2］。針對國家一些重大的測繪生產任務，如西部測圖任務，在坡度陡、高差大的西部地區(qū)進行高精度的DEM制圖已經成為一項極具挑戰(zhàn)性的任務。同時，由于TerraSAR－X、COSMO－SkyMed等高分率雷達遙感衛(wèi)星的相繼發(fā)射，這些具有強干涉能力的雷達傳感器使高精度的DEM地形制圖已經成為可能。

在應用InSAR技術生成DEM過程中，主要的誤差源來自于相位誤差［3－5］。相位誤差主要由三部分組成：① 軌道誤差引起的相位趨勢誤差；② 由于時間、空間去相關和熱噪聲造成的誤差；③由于主從影像成像時間產生的大氣差的擾動所引起的誤差。為了更加精確地去除相位誤差，外部DEM已經被應用到雷達干涉測量中，文獻［6］提出了一種用粗的DEM方法減少干涉圖的局部相位帶寬，較低的相位帶寬減少了相位殘余，從而使相位解纏更加容易；文獻［7］則有用外部的DEM減少高程的搜索范圍；文獻［8］使用SRTM通過構建線性模型最大程度地去除了相位誤差。然而，針對坡度陡、高差大的復雜地形地區(qū)，單一的線性模型已無法擬合整個地形趨勢，同時，在高分辨率雷達數(shù)據處理步驟中，如建立外部DEM模擬相位和干涉解纏相位的對應關系，在軌道參數(shù)精度低的條件限制下也嚴重影響最終的DEM精度。

本文針對坡度陡、高差大的復雜地形地區(qū)高精度DEM制圖，首先根據雷達干涉影像的誤差相位構建線性模型，在此基礎上，研究一種外部DEM輔助下的DEM精化方法，最后應用于高分辨雷達影像（COSMO數(shù)據）在云南德欽地區(qū)的1∶50 000DEM制圖。

2 相位誤差分析

雷達影像干涉相位的表達式為

式中，φflt代表干涉相位中由于平地效應引起的相位；φz代表外部DEM產生的干涉相位；φn代表系統(tǒng)噪聲引起的相位；φaps代表主從影像成像時間產生的大氣差的殘余相位；φtrd能夠擬合由于軌道誤差產生的相位誤差和地形所產生的大氣相位差的線性趨勢，可以表示為［9］

式中，c、l1、l2和l3代表線性趨勢φtrd的線性系數(shù)。在把φflt和φz從干涉相位減去后，得到的相位差表示為

式中，φerr代表外部DEM引起的相位誤差。這部分相位差將通過線性回歸分析進行去除。在回歸分析中，假設外部DEM引起的相位誤差服從均值為零，標準方差為常量的標準正態(tài)分布。由于在回歸分析有大量的樣本點參與計算，φn和φaps的影響也基本可以忽略。因此，相位差的主要部分是φtrd。

除了以上建模的相位誤差以外，在影像的某些像素上還存在解纏相位殘余、雷達陰影等引起的相位粗差，這些粗差也直接影響最終DEM的產品精度。這部分粗差根據生成的DEM與外部DEM的距離范圍進行濾波處理。

根據上述分析，外部DEM輔助下的DEM精化方法的核心是線性回歸分析和系統(tǒng)粗差的濾除。

3 外部DEM輔助下的DEM精化方法

通過第二節(jié)對干涉相位的誤差分析，外部DEM輔助下的DEM精化方法主要由建立外部DEM模擬相位和干涉解纏相位的對應關系，多模型的線性回歸分析和高程點的濾波，如圖1所示。

3.1 干涉圖和外部DEM的精確對應關系建立

首先，通過應用軌道參數(shù)和高程值計算多普勒、距離和橢球方程來建立干涉圖和外部DEM的初始的對應關系。然后，將真實的幅度和模擬幅度影像進行配準從而得到更加精確的對應關系。具體步驟為：

（1）初始地理編碼表建立。應用軌道參數(shù)和高程值聯(lián)合求解多普勒、距離和橢球方程將外部DEM從地理坐標的轉化到雷達距離－多普勒坐標系下［10］，外部DEM的地理坐標和雷達坐標的初始對應關系被存儲在地理編碼表內。

（2）將雷達距離－多普勒坐標系下的外部DEM進行幅度影像模擬。

圖1 外部DEM輔助下的DEM精化處理流程Fig.1 Working flow of DEM refinement facilitated by external DEM

Muhleman模型被用于描述雷達的后向散射系數(shù)σ［11］，對于不同的局部入射角η，有

式中，K是常量，在雷達影像的某些區(qū)域，如疊掩和陰影區(qū)，Ps僅能反映后向散射的部分能量。

（3）真實的幅度影像和模擬幅度影像配準。多級配準技術將被應用于真實的幅度影像和模擬幅度影像配準。多級配準分為影像的粗配準和精配準。粗配準主要根據影像互相關函數(shù)的統(tǒng)計特性，通過在空間域內尋找兩幅影像的互相關函數(shù)最大值進行配準，其精度可達到像素級。精配準對每一對粗配準的數(shù)據塊在頻率域內運用互相關函數(shù)進行配準，其精度可達到子像素級。一旦得到每個數(shù)據塊的準確偏移量，雙線性的配準多項式便可以確定。

（4）地理編碼表的精化。真實的幅度影像和模擬幅度影像的精確對應關系確定下來后，便可用確定的二次多項式精化初始的地理編碼表。利用精化后的地理編碼表將外部DEM從地理坐標重新轉化到雷達距離－多普勒坐標系下，外部DEM和干涉影像便建立了準確的一一對應關系。

3.2 多模型的線性回歸分析

考慮到在高山地區(qū)地形的復雜性，在這種情況下多個回歸模型將被采用。通過區(qū)域增長算法［12］把解纏后的干涉相位分割成幾個子影像塊，子塊數(shù)目一般不宜過多（不超過10個），主要依據地形的變化進行分割。在每一個影像塊上，按照下面的步驟進行回歸分析。首先，把雷達距離－多普勒坐標系下的外部DEM按照式（6）轉化成干涉相位，即

式中，φz是從外部DEM轉化而來的干涉解纏相位；B是基線長度；R是斜距方向上目標到傳感器的距離；θ是入射角；α是基線與水平方向的夾角；Z是外部DEM的高程值。

然后，選取相干性大于一定閾值（一般為0.5）的像素點作為可靠的點去擬合線性模型。對于分割的不同影像塊，平均的相干性作為參考的閾值。一旦得到影像上每個像素的相位差Δφ，多個線性回歸模型可寫為

n和m分別代表每個模型像素點的個數(shù)和模型的個數(shù)，模型擬合系數(shù)m）能夠由最小二乘擬合得到。為了正確地估計相位趨勢，相位解纏中的殘余相位將通過迭代過程逐步地被去除。即殘余相位大于兩倍標準差的像素點將被剔除，并且相位趨勢迭代的被估計值到所有的像素點的殘余相位都小于一定的閾值。應用多模型的線性回歸，相位趨勢能被逐像素地、準確地從干涉圖中估計并去除。

3.3 不可靠的高程點的濾除

在去除了相位趨勢以后，解纏干涉圖能夠轉化成高程圖，高程通過式（8）計算得到

對于解纏干涉圖的系統(tǒng)粗差的去除，主要依賴于InSAR得到的高程值與外部DEM的距離范圍?；贑hebyshev理論，無論誤差服從什么分布，任何誤差在區(qū)間μ±4δ的概率至少是94%（μ和δ分別是誤差的均值和標準偏差）。在該方法中，設置4倍的外部DEM誤差作為去除不可靠高程點的標準。計算由InSAR和對應的外部DEM得到的高程點的絕對高程差，選擇絕對高程差小于4倍的外部DEM誤差的高程點并用于重構最終的DEM。

不可靠的高程點濾除后，由于InSAR影像在某些區(qū)域相干性較低，容易在影像中產生大的“空洞”（一系列連續(xù)的不可靠的高程點）。這里定義超過20個像素的連續(xù)區(qū)域定義為影像空洞，不可靠的高程點將被保留。而在較小的連續(xù)區(qū)域，其高程點將通過雙線性插值方法得到。

4 試驗與分析

4.1 試驗區(qū)和數(shù)據介紹

以云南的德欽地區(qū)作為試驗區(qū)，該地區(qū)的地形非常復雜：高差達到5 000m以上，平均坡度達到40°以上，同時該地區(qū)常年被冰雪覆蓋。雷達干涉測量作為該地區(qū)的重要DEM測量手段和方法，選取相隔一天的3m分辨率COSMOSkyMed雷達干涉數(shù)據對，覆蓋云南德欽地區(qū)約為1 300km2，具體數(shù)據參數(shù)見表1。

表1 數(shù)據集的基本信息Tab.1 Basic information of data sets

4.2 試驗結果

在該試驗中引入1∶100 000的外部DEM，應用軌道參數(shù)和高程值建立的初始地理編碼表將外部DEM從地理坐標轉化到雷達距離－多普勒坐標系下，同時進行SAR幅度影像模擬，如圖2所示。從圖2可以看出，經過影像模擬以后，模擬的SAR幅度影像和真實的SAR在細節(jié)上具有相似的紋理。

圖2 模擬的SAR幅度影像和真實的SAR幅度影像比較Fig.2 Comparision of SAR simulated amplitude image and real amplitude image

利用模擬的SAR幅度影像和真實的SAR幅度影像的對應關系精化地理編碼表，得到在幾何關系上和干涉圖一一對應的外部DEM，進而將外部DEM轉化成解纏的干涉相位，如圖3所示。最小費用流（MCF）算法用于干涉圖的解纏［13］，得到的解纏后的干涉相位圖，如圖4所示。

圖3 外部DEM的干涉相位圖Fig.3 Phase with external DEM

圖4 解纏后的干涉相位圖Fig.4 Unwrapping phase

引入的1∶100 000外部DEM在高山地區(qū)的高程精度優(yōu)于28m［14］，InSAR高程點的濾波閾值設置為4倍的外部DEM精度，約100m。將解纏后的干涉圖根據區(qū)域增長算法分割成3個子影像塊，如圖5所示。在每個子數(shù)據塊上選取相干性大于0.5的像素點進行線性回歸分析，逐步地把相位趨勢去除。在相位向高程轉化過程中，InSAR生成的DEM和對應的外部DEM相減得到的絕對高程差大于100m的高程點將被濾除。理論上，有4 801 458個高程點（總共有16 444 755個高程點）將被濾除，但由于影像上存在“空洞”，實際上由4 695 825個點被濾除，其余的點保留InSAR生成的高程結果。最終的1∶50 000DEM的結果如圖6所示，從局部的影像信息可以看出，精化后的DEM比未經過精化處理的DEM在細節(jié)信息上有了較大改進。

圖5 解纏干涉相位分割結果Fig.5 Region segmentation result of phase

圖6 1∶50 000 DEM結果 Fig.6 1∶50 000DEM result

4.3 結果驗證

為了精確地驗證所生成的1∶50 000DEM，23個驗證點從GPS測量得到，分布如圖6所示（三角形標記驗證點），GPS點的測量精度達到毫米級。為了比較DEM精化方法和傳統(tǒng)InSAR方法得到的DEM精度，在23個GPS檢查點上進行高程的誤差比較。如圖7所示，由DEM精化方法得到高程精度比傳統(tǒng)的InSAR方法有了一定程度的改善，同時高程誤差的分布也更加均勻。通過在這些檢查點上計算得到，傳統(tǒng)InSAR方法得到的DEM的標準偏差為40.9m，DEM精化方法得到的DEM的標準偏差為19.7m。

圖7 在GPS檢查點上的高程誤差分布Fig.7 Height errors distribution on check points

4.4 計算效率分析

為了分析該方法的計算效率，將SAR影像分成不同面積大小的數(shù)據塊，使用配置為雙核2.66GHz的CPU，內存2G的計算機進行處理。如圖8所示，DEM精化方法的處理時間是隨著SAR影像的處理面積而增大的，處理一整景SAR影像的時間約為8 000s。在DEM精化處理中，多模型的線性回歸分析所占的時間最長，約為整個處理時間的三分之二，它主要取決于選取參與回歸計算的像素點數(shù)目。

圖8 DEM精化方法的處理時間Fig.8 Processing time of DEM refinement

5 結 論

針對在地形復雜區(qū)域的高精度雷達制圖的需要，發(fā)展了一種基于高分率雷達影像的外部DEM輔助下的復雜地形制圖方法。本方法的特點之一是實現(xiàn)了外部DEM和干涉圖精確的一一對應關系，另一特點則是運用多模型的線性回歸策略擬合了地形趨勢并去除了相位誤差，也利用了引入的外部DEM很好地去除了系統(tǒng)粗差。通過精度的比較和分析，本方法得到的DEM精度比傳統(tǒng)的InSAR方法有了一定程度的提高，其高程精度（19.7m）能滿足國家1∶50 000DEM地形制圖要求［14］，但計算效率有待提高，今后可采用多任務并行的計算方式改善其計算效率。云南德欽地區(qū)作為復雜地形的代表試驗區(qū)，其試驗結果證實了該方法在復雜地形區(qū)域應用高分辨雷達影像的制圖能力。同時，該方法也已經推廣至國家西部測圖項目的生產任務中，適用于平地、丘陵和高山等各種地形區(qū)域，在應用星載TerraSAR、Radarsat－2和COSMO－SkyMed等高分辨干涉SAR影像進行DEM制圖上已經取得了較好的效果。

［1］ GUO Huadong.Theories and Application of Radar for Earth Observation［M］.Beijing：Science Press，2000：131－135.（郭華東.雷達對地觀測理論與應用［M］.北京：科學出版社，2000：131－135.）

［2］ ZHOU Jianmin，LI Zhen，LI Xinwu.Research on Rules of the Valley Glacier Motion in Western China Based on ALOS／PALSAR Interferometry［J］.Acta Geodaetica et Cartographica Sinica，2009，38（4）：341－347.（周建民，李震，李新武.基于ALOS／PALSAR雷達干涉數(shù)據的中國西部山谷冰川冰流運動規(guī)律研究［J］.測繪學報，2009，38（4）：341－347.）

［3］ ZEBKER H A，VILLASENOR J.Decorrelation in Interferometric Radar Echoes［J］.IEEE Transactions on Geoscience and Remote Sensing，1992，30（5）：950－959.

［4］ ZEBKER H A，WERNER C L，ROSEN P A，et al.Accuracy of Topographic Maps Derived from ERS－1Interferometric Radar［J］.IEEE Transactions on Geoscience and Remote Sensing，1994，32（4）：823－836.

［5］ ROSEN P A，HENSLEY S，JOUGHIN I R，et al.Synthetic Aperture Radar Interferometry［EB／OL］.［2010－01－12］.http：∥www.gps.caltech.edu／classes／ge167／file／Rosen＿IEEE2000＿Insar.pdf.

［6］ SEYMOUR M S，CUMMING I G.InSAR Terrain Height Estimation Using Low－quality Sparse DEM’s［C／OL］∥Proceedings of the 3rd ESA Scientific Symposium on ERS Satellites，F(xiàn)lorence：［s.n.］，1997［2010－01－12］.http：∥florence97.ers－symposium.org／.

［7］ EINEDER M，ADAM N.A Maximum－likelihood Estimator to Simultaneously Unwrap，Geocode，and Fuse SAR Interferograms from Different Viewing Geometries into One Digital Elevation Model［J］.IEEE Transactions on Geoscience and Remote Sensing，2005，43（1）：24－36.

［8］ LIAO Mingsheng，WANG Teng，LU Lijun，et al.Reconstruction of DEMs from ERS－1／2Tandem Data in Mountainous Area Facilitated by SRTM Data［J］.IEEE Transactions on Geoscience and Remote Sensing，2007，45（7）：2325－2335.

［9］ HANSSEN R F.Radar Interferometry：Data Interpretation and Error Analysis［M］.Dordrecht：Kluwer Academic Publishers，2001.

［10］ GEUDTNER D.The Interferometric Processing of ERS－1 SAR Data［R］.［S.l.］：European Space Agency.2006.

［11］ MUHLEMAN D O.Radar Scattering from Venus and the Moon［J］.The Astronomical Journal，1964，69：34－41.

［12］ ADAMS R，BISCHOF L.Seeded Region Growing［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1994，16（6）：641－647.

［13］ COSTANTIN M.A Novel Phase Unwrapping Method Based on Network Programming［J］.IEEE Transactions on Geoscience and Remote Sensing，1998，36（3）：813－821.

［14］ State Bureau of Surveying and Mapping.CH＼T1008－2001 Digital Products of Fundamental Geographic Information 1∶10 000，1∶50 000Digital Elevation Models［S］.Beijing：Surveying and Mapping Press，2001.（國家測繪局.CH＼T1008－2001基礎地理信息數(shù)字產品1∶10 000、1∶50 000數(shù)字高程模型［S］.北京：測繪出版社，2001.）