盧麗君，張繼賢，王 騰

1.中國(guó)測(cè)繪科學(xué)研究院，北京100830；2.武漢大學(xué)測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北武漢430079

1 引 言

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

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

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

2 相位誤差分析

雷達(dá)影像干涉相位的表達(dá)式為

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

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

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

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

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

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

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

3.1 干涉圖和外部DEM的精確對(duì)應(yīng)關(guān)系建立

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

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

（2）將雷達(dá)距離－多普勒坐標(biāo)系下的外部DEM進(jìn)行幅度影像模擬。

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

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

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

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

（4）地理編碼表的精化。真實(shí)的幅度影像和模擬幅度影像的精確對(duì)應(yīng)關(guān)系確定下來后，便可用確定的二次多項(xiàng)式精化初始的地理編碼表。利用精化后的地理編碼表將外部DEM從地理坐標(biāo)重新轉(zhuǎn)化到雷達(dá)距離－多普勒坐標(biāo)系下，外部DEM和干涉影像便建立了準(zhǔn)確的一一對(duì)應(yīng)關(guān)系。

3.2 多模型的線性回歸分析

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

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

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

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

3.3 不可靠的高程點(diǎn)的濾除

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

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

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

4 試驗(yàn)與分析

4.1 試驗(yàn)區(qū)和數(shù)據(jù)介紹

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

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

4.2 試驗(yàn)結(jié)果

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

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

利用模擬的SAR幅度影像和真實(shí)的SAR幅度影像的對(duì)應(yīng)關(guān)系精化地理編碼表，得到在幾何關(guān)系上和干涉圖一一對(duì)應(yīng)的外部DEM，進(jìn)而將外部DEM轉(zhuǎn)化成解纏的干涉相位，如圖3所示。最小費(fèi)用流（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高程點(diǎn)的濾波閾值設(shè)置為4倍的外部DEM精度，約100m。將解纏后的干涉圖根據(jù)區(qū)域增長(zhǎng)算法分割成3個(gè)子影像塊，如圖5所示。在每個(gè)子數(shù)據(jù)塊上選取相干性大于0.5的像素點(diǎn)進(jìn)行線性回歸分析，逐步地把相位趨勢(shì)去除。在相位向高程轉(zhuǎn)化過程中，InSAR生成的DEM和對(duì)應(yīng)的外部DEM相減得到的絕對(duì)高程差大于100m的高程點(diǎn)將被濾除。理論上，有4 801 458個(gè)高程點(diǎn)（總共有16 444 755個(gè)高程點(diǎn)）將被濾除，但由于影像上存在“空洞”，實(shí)際上由4 695 825個(gè)點(diǎn)被濾除，其余的點(diǎn)保留InSAR生成的高程結(jié)果。最終的1∶50 000DEM的結(jié)果如圖6所示，從局部的影像信息可以看出，精化后的DEM比未經(jīng)過精化處理的DEM在細(xì)節(jié)信息上有了較大改進(jìn)。

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

圖6 1∶50 000 DEM結(jié)果 Fig.6 1∶50 000DEM result

4.3 結(jié)果驗(yàn)證

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

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

4.4 計(jì)算效率分析

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

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

5 結(jié) 論

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

［1］ GUO Huadong.Theories and Application of Radar for Earth Observation［M］.Beijing：Science Press，2000：131－135.（郭華東.雷達(dá)對(duì)地觀測(cè)理論與應(yīng)用［M］.北京：科學(xué)出版社，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雷達(dá)干涉數(shù)據(jù)的中國(guó)西部山谷冰川冰流運(yùn)動(dòng)規(guī)律研究［J］.測(cè)繪學(xué)報(bào)，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.（國(guó)家測(cè)繪局.CH＼T1008－2001基礎(chǔ)地理信息數(shù)字產(chǎn)品1∶10 000、1∶50 000數(shù)字高程模型［S］.北京：測(cè)繪出版社，2001.）