薛東劍，鄭 潔，李成繞，李婉秋，蘇璐璐

(1.成都理工大學(xué) 地球科學(xué)學(xué)院，四川 成都 610059；2.中國林業(yè)科學(xué)研究院 資源信息研究所,北京 100091；3.國土資源部地學(xué)空間信息技術(shù)重點實驗室,四川 成都 610059)

利用L波段星載重復(fù)軌道干涉SAR提取DEM及大氣效應(yīng)分析

薛東劍1,2,3，鄭 潔1,3，李成繞1,3，李婉秋1,3，蘇璐璐1,3

(1.成都理工大學(xué) 地球科學(xué)學(xué)院，四川 成都 610059；2.中國林業(yè)科學(xué)研究院 資源信息研究所,北京 100091；3.國土資源部地學(xué)空間信息技術(shù)重點實驗室,四川 成都 610059)

以四川省茂縣、安縣地區(qū)的ALOS PALSAR L波段雷達影像為數(shù)據(jù)源，在對地形復(fù)雜區(qū)的配準(zhǔn)、解纏等算法研究的基礎(chǔ)上，結(jié)合SAR特有的成像幾何結(jié)構(gòu)，對兩幅SAR圖像進行快速自動配準(zhǔn)，配準(zhǔn)誤差小于0.2個像元；在去除平地引起的干涉相位變化后，運用MCF對纏繞相位進行解纏，在此基礎(chǔ)上提取研究區(qū)的數(shù)字高程模型(DEM)，分析發(fā)現(xiàn)精度受多種因素影響，其中波長較長的L波段數(shù)據(jù)比波長較短的時間相關(guān)性好,此外為消除大氣效應(yīng)，采用改進的相位累積法去除大氣的影響，在一定程度上提高DEM精度。其誤差范圍在±10 m左右，對快速提取地形信息具有一定的借鑒意義。

PALSAR；InSAR；配準(zhǔn)；相位解纏；大氣效應(yīng)

以往運用遙感技術(shù)提取地形信息，多采用光學(xué)圖像，通過立體相對進行提??；圖像畸變及天氣條件嚴(yán)重限制了其應(yīng)用。到20世紀(jì)60年代末，合成孔徑雷達干涉測量技術(shù)(InSAR)應(yīng)運而生，其能夠全天候、全天時工作，精確度高，具有光學(xué)圖像無法比擬的優(yōu)勢，有著廣闊的應(yīng)用前景。Graham(1974)等[1]利用機載SAR數(shù)據(jù)獲取滿足1∶25 萬地形圖要求的高程數(shù)據(jù)，開創(chuàng)InSAR技術(shù)在對地觀測中獲取三維信息的先河。隨著星載合成孔徑雷達的不斷升空，產(chǎn)生大量豐富的SAR數(shù)據(jù)，雷達干涉被廣泛地應(yīng)用于提取地形信息[2-5]。研究表明, 波長較長的L波段SAR 數(shù)據(jù)比波長較短的X、C波段數(shù)據(jù), 在干涉測量時, 時間相關(guān)性好, 更適合于地形起伏大、植被覆蓋密、形變梯度大的區(qū)域[3]，尤其有利于西南地區(qū)高植被覆蓋區(qū)地形提取。并且ALOS PALSAR FBS數(shù)據(jù)的干涉臨界基線距較大，能夠使用更多的影像參與干涉處理。由于SAR特有的成像幾何結(jié)構(gòu)，其處理中的很多關(guān)鍵技術(shù)還一直處于不斷研究與發(fā)展中，尤其針對不同的研究區(qū)，其影響提取精度的處理是目前研究的焦點與熱點(如解纏、大氣效應(yīng)等)。采用茂縣地區(qū)的L波段ALOS PALSAR數(shù)據(jù)，在對配準(zhǔn)、解纏算法等研究的基礎(chǔ)上，提取了數(shù)字高程模型(DEM)，并重點分析了大氣效應(yīng)對提取精度的影響。采用該技術(shù)可快速提取地形信息，這對于地形變化監(jiān)測及由地形引起的災(zāi)害監(jiān)測預(yù)報、環(huán)境治理等問題具有重要的實用價值。

1 研究區(qū)概況及數(shù)據(jù)預(yù)處理

研究區(qū)位于四川省東北部的茂縣、安縣地區(qū)，地處龍門山中南段與四川盆地相結(jié)合區(qū)域，地理坐標(biāo)：東經(jīng)103°35′47″～104°22′40″，北緯 31°23′48″～32°00′13″，是“5.12”地震中受災(zāi)較為嚴(yán)重地帶，對其震后地形的提取具有一定的應(yīng)用價值。研究區(qū)地形主要為平壩、丘陵(臺地)、中低山三種類型，既以大光包斜沖斷層和北川沖斷層為界，西北部屬四川西部地槽區(qū)的前龍門山褶斷帶，系龍門山脈，地勢較高，最高峰位于九頂山主峰高4 989 m，為區(qū)內(nèi)最高點，最低點位于安縣秀水鎮(zhèn)，海拔為548 m。

研究中采用2010-01-24與2010-03-11 1.1級L波段(波長λ：0.236 m)ALOS PALSAR數(shù)據(jù)，每景數(shù)據(jù)覆蓋范圍約為2 565.2 km2。主輔單視復(fù)數(shù)(SLC)影像距離向像元大小4.68 m ，方位向3.18 m，9 344列，18 432行(見表1)，其中以2010-01-24數(shù)據(jù)為基準(zhǔn)(升軌)，進行干涉處理(見圖1，為研究區(qū)ALOS PALSAR HH數(shù)據(jù)，為顯示方便，在此進行了多視處理及地理編碼)；參考DEM數(shù)據(jù)選擇美國NASA和NIMA聯(lián)合測量的SRTM數(shù)據(jù)。

圖1 研究區(qū)ALOS PALSAR數(shù)據(jù)示意

成像時間中心經(jīng)緯度/(°)分辨率/m緯度經(jīng)度距離向方位向入射角/(°)近距點斜距/km極化方式時間基線/d空間基線/m2010012431．6943648103．97570492010031131．6944510103．9831671468318387402387474846717HH46479

2 InSAR關(guān)鍵處理技術(shù)

利用星載SAR重復(fù)軌道方式進行干涉測量，主要是通過重復(fù)兩次觀測，獲得同一區(qū)域的單視復(fù)數(shù)影像對，根據(jù)兩幅天線和觀測目標(biāo)之間的幾何關(guān)系(見圖2)，結(jié)合觀測平臺的軌道參數(shù)，在對相位差解算的基礎(chǔ)上，提取地面高程。

圖2 雷達干涉測量原理示意圖

由圖中幾何關(guān)系可推算高程：

h=H-ρcosθ.

(1)

(ρ+Δρ)2=ρ2+B2-2ρBcos(90-θ+α).

(2)

(3)

(4)

圖2中，A1，A2為天線在兩個成像時刻的位置；H是A1的高度；B是基線，α為基線與水平方向的夾角；θ為天線在A1處的視角，λ為波長；ρ，ρ+Δρ為A1,A2到觀測點Z的斜距；h為目標(biāo)點高程。利用重復(fù)軌道獲取的相位差僅是相位差值中的小數(shù)部分，需聯(lián)合相位解纏求出φt。由于SAR特有的成像方式，使其在數(shù)據(jù)獲取上具有一定要求，且處理上也有一定的難度，尤其是對兩幅圖像的配準(zhǔn)及相位解纏的好壞直接關(guān)系到提取地形的精度。

2.1 SAR圖像自動配準(zhǔn)

對不同天線獲得數(shù)據(jù)進行精確配準(zhǔn)，估測兩幅SLC影像之間距離和方位向的偏差，是利用SAR進行干涉測量提取DEM的關(guān)鍵步驟。為了準(zhǔn)確計算對應(yīng)的地面目標(biāo)的干涉相位或形變相位，需要實現(xiàn)雷達圖像之間亞像元級精度的配準(zhǔn)[12]。在光學(xué)圖像上，自動配準(zhǔn)的基本思想是相同的地物具有相似的光譜特征，通過對兩個圖像做相對移動，找出其相似性量度值最大或差別最小的位置作為圖像配準(zhǔn)的位置。SAR圖像的配準(zhǔn)也是類似的原理，目前國內(nèi)外大部分配準(zhǔn)方法基本上是在確定控制的基礎(chǔ)上，以控制點坐標(biāo)擬合多項式：

(5)

其中，xr,yr為參考影像中控制點坐標(biāo)；Xt,Yt為待配準(zhǔn)數(shù)據(jù)相應(yīng)控制點坐標(biāo)。如何精準(zhǔn)確定控制點，是配準(zhǔn)處理的前提，也是關(guān)鍵，文中針對ALOS PALSAR成像特征，綜合分析相位差平均強度梯度法、相干系數(shù)法、絕對差值法、最大頻譜法的基礎(chǔ)上，采用相干系數(shù)法進行了控制點選取，此方法充分利用了振幅與相位信息，在一定的搜索范圍內(nèi)，逐一計算每一點的相關(guān)系數(shù)，以最大值為配準(zhǔn)點的位置。為了避免產(chǎn)生模糊問題，取得精確的偏差估計，在多視估計后，進行一次單視處理，并采用改進的最小二乘法生成配準(zhǔn)多項式(表2)，從而使配準(zhǔn)誤差小于0.2個像元，提高干涉相關(guān)性。

表2 干涉相對配準(zhǔn)擬合參數(shù)

2.2 干涉圖生成及相位解纏

對干涉像對進行配準(zhǔn)后，SLC1與SLC2影像的數(shù)據(jù)可表示為

μ1=|μ1|ejφ1，μ2=|μ2|ejφ2.

(6)

通過共軛相乘，可以計算出同名點上的相位差：

|μ2||μ2|e-j(φ2-φ1).

(7)

其中，μ*表示復(fù)數(shù)的共軛，μint的相位是每一同名點上相應(yīng)于兩個觀測點的相位差(φ)，所得相位的主值數(shù)據(jù)就是被纏繞相位數(shù)據(jù)，以干涉圖形式表示。運用InSAR技術(shù)提取DEM的精度是否成功，關(guān)鍵是干涉像對的相位相關(guān)性的高低。經(jīng)處理，在本次實驗中，相干系數(shù)最大值是0.989，最小值是0.001 4，平均值是0.691 6， 0.4以上相關(guān)性占78%，該試驗區(qū)干涉像對的整體相干性較好(圖3)。此時生成的干涉圖中是由多個成分組成，需去平處理，把干涉圖中的地球彎曲距離和方位向的相位趨勢去除，去除了平地相位之后，干涉條紋由密集變?yōu)橄∈?，避免了將平地相位信息作為誤差引入高程和形變信息。從圖4可以看出，去平后干涉圖，在地形比較平坦的地區(qū)(圖4右上角地區(qū))干涉條紋已經(jīng)去除，這與實際地形比較相符。

圖3 相干圖與強度合成

圖4 去平后干涉

相位解纏在干涉處理中具有重要的地位，而方法的選擇會很大程度影響高程的精度。常用的算法有區(qū)域增長法，最小費用流量法等，其中，區(qū)域增長算法，是根據(jù)相位的質(zhì)量來規(guī)劃積分路徑，先選擇高質(zhì)量的區(qū)域進行積分，然后再向低質(zhì)量的區(qū)域，從而得到解纏圖像。這種方法在相位圖質(zhì)量較高時，解纏結(jié)果好，速度快，但無法識別殘差點，會引入相位誤差影響精度。最小費用流法是一種基于網(wǎng)絡(luò)規(guī)劃的相位解纏方法，其基本思想是最小化解纏后相位的導(dǎo)數(shù)與纏繞相位的導(dǎo)數(shù)之間的差異，該方法可兼顧精度與處理效率，尤其對范圍較大地區(qū)的解纏具有一定的優(yōu)勢。在此以最小費用流量(MCF)算法設(shè)置的閥值作為標(biāo)準(zhǔn)(在此相干閾值設(shè)為0.3，圖5)，將所有相干性滿足條件的相位添加到一個集合，在集合中建立三角網(wǎng)，構(gòu)造對偶圖，用最小費用流法連接各殘差點，求出最小費用流集合，再對相位圖進行積分，根據(jù)相位質(zhì)量從高到低進行解纏。

圖5 相位解纏結(jié)果

3 地形信息提取及大氣效應(yīng)分析

經(jīng)過相位解纏圖像(圖5)的相位主值已經(jīng)轉(zhuǎn)化成了真實的相位，采用式3計算出每個像元的高程，在此基礎(chǔ)上將獲取的高程信息由雷達坐標(biāo)轉(zhuǎn)化為地理坐標(biāo)(圖6)。為了分析提取出的DEM精度，采用SRTM DEM數(shù)據(jù)做參考，進行比較分析。SRTM數(shù)據(jù)在平原區(qū)精度優(yōu)于1∶250 000 DME,中誤差為4.921。

圖6 InSAR提取的DEM與SRTM DEM疊合分析

運用InSAR技術(shù)提取的DEM與SRTM DEM數(shù)據(jù)比較發(fā)現(xiàn)，整體趨勢比較吻合(圖6)，平坦地區(qū)誤差范圍較小，而坡度較大的地區(qū)，誤差相對較大，存在部分異常值，但差值大的像元比較少。①研究區(qū)大部位于山區(qū)，有植被覆蓋的山區(qū)受時間失相關(guān)的影響，對提取的精度造成一定的影響，而城市等平坦地區(qū)受時間失相關(guān)的影響較小, 其相關(guān)性較好；②雷達干涉測量中，除了時間去相干外，大氣效應(yīng)也是重要的影響因素，不同時間獲取的圖像受大氣影響產(chǎn)生不同的相位延遲，干涉相對中兩個SAR圖像回波信號的相位可表示為

(8)

(9)

其中，ρ1，ρ2為斜距，Δρ1，Δρ2為雷達電磁波穿過大氣層(電離層與對流層)，受電子濃度、氣壓、溫度、水汽等影響造成的大氣延遲，λ為波長。Δφ為兩次觀測大氣引起的相位差：

(10)

大氣造成的延遲主要發(fā)生在電離層與對流層。以往研究表明，大氣影響主要是由于水汽的變化引起的[17-18]，表征大氣水汽的物理量主要為可降水汽含量PWV(Precipitable Water Vapor)，可以轉(zhuǎn)換為天頂濕度延遲ZWD(Zenith Wet Delay)：

(11)

其中，K為折射常數(shù)，ρw為液態(tài)水密度，TM為對流層的加權(quán)平均溫度，Rv為大氣常數(shù)。針對重復(fù)軌道干涉測量中高程與相位之間的參量關(guān)系，可推導(dǎo)出：

(12)

σh為高程誤差，σδ為天頂濕度延遲誤差，θ為入射角，B⊥為垂直基線。

綜合目前研究現(xiàn)狀[15-20]，InSAR測量中常用的減少大氣影響的方法，除了數(shù)據(jù)本身因素外，主要有相位累積法，PS技術(shù)，基于GPS數(shù)據(jù)的改正方法以及利用外部水汽產(chǎn)品(如MERIS、FY、MODIS等)校正法等，其中PS技術(shù)、GPS校正及外部水汽校正法，雖然普適性較高，但主要適用于地表形變監(jiān)測中[21]。結(jié)合ALOS PALSAR數(shù)據(jù)特點，且考慮到外部數(shù)據(jù)獲取的同步性及相位解纏的難易等因素，文中借鑒相位累積法，采用2010年ALOS PALSAR數(shù)據(jù)生成三對獨立干涉圖進行相位梯度累加分析，相位梯度可表示為

(13)

(14)

為減少噪聲的影響，對干涉圖的相位梯度減去平均值和垂直基線的乘積,得到相位梯度變化

(15)

則大氣引起的相位延遲異常▽φatm為

(16)

▽φvari為相位梯度異常,γi為空間相關(guān)系數(shù)。通過該相位梯度累積平均一定程度上降低了大氣的影響。大氣改正后的高程差異值明顯低于未經(jīng)過大氣改成的高程差異值，且高程差異分布相對比較集中，高程差異的標(biāo)準(zhǔn)差由15.65 m提高到了6.02 m，尤其是平坦地區(qū)誤差范圍在±10 m以內(nèi)的置信度由68.78%提高到了93.21%，DEM精度得到提高。

4 結(jié)束語

以ALOS PALSAR數(shù)據(jù)為例，在對InSAR處理流程及算法詳細(xì)分析的基礎(chǔ)上，研究了地形復(fù)雜區(qū)域SAR圖像配準(zhǔn)、解纏等核心算法；以最小費用流(MCF)算法為基礎(chǔ)，將相干性滿足閾值條件的相位添加到一個集合，構(gòu)造對偶圖，用最小費用流法連接各殘差點，對相位進行積分，根據(jù)相位質(zhì)量從高到低進行解纏，并提取了DEM。研究表明利用雷達干涉測量提取地形信息具有自動化程度高、效率較高等特點，但失相關(guān)是限制INSAR技術(shù)應(yīng)用的一個較為嚴(yán)重的問題，尤其植被覆蓋度大、地形復(fù)雜地區(qū)，時間去相關(guān)明顯；除此外大氣效應(yīng)也是重要的影響因素，在地形提取中，運用相位梯度累加在一定程度上可以消除大氣的影響，提高DEM的精度； ALOS數(shù)據(jù)由于波長較長，對地表穿透深度大，干涉性相對較好，有利于西南地區(qū)高植被覆蓋區(qū)地形提取。

[1] GRAHAM L C.Synthetic Interferometer Radar for Topographic Mapping [A].Proceeding of the IEEE,1974,62(6)[C]. [ s.l.]:[ s.n.],1974.763-770.

[2] GOLDSTEIN R M, ZERBER H A,WERNER C L.Satellite Radar Interferometry :Twodimensional PhaseUnwrapping[J]. Radio Science,1988,23(4):713-720.

[3] RAUCOULES D, COLESANTI C, CARNEC C. Use of SAR interfer ometr y fo r detecting and assessing g ro und subsidence[ J] . C.R. Geo science, 2007( 339): 289-302.

[4] ZEBKER H A, ROSEN P A, Hensley S. Atmospheric Effects in Interferometric Synthetic Aperture Radar Surface Deformation and Topographic Maps [J]. Journal of Geophysical Research,1997,103( B4): 7547-7563.

[5] HSIEH Chia-Sheng,SHIH Tian-Yuan,CHING Jyr,et al. Using differential SAR interferometry to map land subsidence: a case study in the Pingtung Plain of SW Taiwan[J], Nat Hazards (2011) 58:1311-1332.

[6] KOCH A,HEIPKE C.Quality assessment of digital surface models derived from the Shuttle Radar Topography Mission(SRTM)[A].the IEEE 2001 International Geoscience and Remote Sensing Symposium[C].Sydney, Australia:2001.

[7] DEO R,MANICKAM S,RAO Y S,et al.Evaluation of interferometric SAR DEMS generated using TanDEM-X data.2013 IEEE International Geoscience and Remote Sensing Symposium(IGARSS)[C].Melbourne,VIC:IEEE,2013:2079-2082.

[8] 陳爾學(xué),李增元,車學(xué)儉.SAR 干涉測量數(shù)字高程模型提取與高程誤差校正[J].高技術(shù)通訊,2000(7):57-63.

[9] 劉國祥,丁曉利,李志林,等.InSAR建立DEM的試驗研究[J].測繪學(xué)報,2001,30(4):337-343.

[10] 吳宏安，張永紅.ALOS/PALSAR InSAR 數(shù)學(xué)模型研究[J].遙感信息,2011(2):106-109.

[11] 沈強,喬學(xué)軍,金銀龍,等.ALOS PALSAR雷達影像InSAR數(shù)據(jù)處理中的基線和地形誤差分析[J].大地測量與地球動力學(xué),2012,32(2):1-6.

[12] 劉廣,郭華東,范景輝.基于外部DEM的InSAR圖像配準(zhǔn)方法研究[J].遙感技術(shù)與應(yīng)用,2008,23(1):72-76.

[13] LEE J S.Digital Image Smoothing and the Sigma Filter.ComputerVision[J]. Graphies and Image Proeessing,1983, 24:255-269.

[14] WANG Qi-jie. Generalized functional model of maximum and minimum detectable deformation gradient for PALSAR interferometry[J].Trans. Nonferrous Met. Soc. China 24(2014) 824-832.

[15] HANSSEN R F．Radar Interferometry：Data Interpretation and Error Analysis[M].Dordreeht：Kluwer Academic Press,2001.

[16] 徐佳,關(guān)澤群,何秀鳳.重復(fù)軌道雷達干涉測量中的大氣影響及其研究展望[J].國土資源遙感,2007,72(2):1-7.

[17] ZEBKER H A,ROSEN P A,HENSLEY S. Atmospheric effects in interferometric synthetic aperture radar surface deformation and topographic maps[J]. Journal of Geophysical Research,1997,102(4):7547-7563.

[18] 孫廣通,張永紅,吳宏安.合成孔徑雷達干涉測量大氣改正研究綜述[J]. 遙感信息,2011(4):111-116.

[19] SONG Xiaogang,LI Deren,SHAN Xinjia,et al. Correction of atmospheric effect in ASAR interferogram using GPS and MODIS data[J]. Chinese Journal of Geophysics, 2009,52(6):1457-1464.

[20] PUYSSEGUR B,MICHEL R,AVOUAC J P,et al.Tropospheric Phase Delay in Interferometric Synthetic Aperture Radar Estimated from Meterological Model and Multispectral Imagery[J].Joural of Geophysical Research,2007,112:120-135.

[21] 崔喜愛,曾琪明,童慶禧.重軌星載InSAR測量中的大氣校正方法綜述[J].遙感技術(shù)與應(yīng)用,2014,29(1):9-17.

[22] GOLDSTEIN R. Atmospheric limitations to repeat-track radar interferometry,G eophys R. es.L ett., 1995,22(18):2517-2520.

[23] SANDWELL D T, PRICE E J. Phase Gradient Approach to Stacking Interferograms [J].Journal of Geophysical Research, 1998, 103 ( B12 ): 30183-30204.

[24] FERRETTI A, PRATI C, ROCCA F.Multibaseline InSAR DEM Reconstruction:the Wavelet Appro ach[J] .IEEE Transactions on Geoscience and Remote Sensing, 1999, 37 (2):705-715.

[25] Ramon Brcic, Alessandro Parizzi, Michael Eineder,et al. Estimation and compensation of onospheric delay for SAR interferometry[C].IEEE International Geoscience & Remote Sensing Symposium,2010,38(1):2908-2911.

UsingrepeatorbitL-bandinterferometricSARtoextractDEMandanalyzingtheatmosphericeffects

XUE Dongjian1,2,3, ZHENG Jie1,3,LI Chengrao1,3, LI Wanqiu1,3, SU Lulu1,3

(1. College of Earth Sciences of CDUT, Chengdu 610059,China；2. Chinese Academy of Forestry, Beijing 100091, China;3.Key Laboratory of Geoscience Spatial Information Technology of Ministry of Land and Resources,Chengdu 610059,China)

This paper, taking ALOS PALSAR L-band radar images as the data source in the counties of Maoxian and Anxian,Sichuan Province，adopts a coherent coefficient method aiming at two SAR images automatic registration and MCF algorithm at phase unwrapping based on the research of registration, unwrapping algorithms, extracted study area’s Digital Elevation Model(DEM). The key technology of interferometric synthetic aperture radar is focused on the basis of the analysis of the current research status. At the time of interferometry, the SAR data with longer wavelength L-band has better temporal coherence than the shorter wavelength.Then the accuracy and characteristics of DEM data are analyzed, whose most error range is about ± 10 m, and the influence of atmospheric effects in interferometriv aperture radar topographic map is discussed, which can provide the basic data for the research work of subsequent surface deformation.

PALSAR;InSAR; image registration;phase unwrapping; atmospheric effect

2016-11-23

四川省教育廳重點資助項目(16ZA0100)；國土資源部地學(xué)空間信息技術(shù)重點實驗室開放基金(KLGSIT2013-06)

薛東劍(1977-)，男，博士.

著錄：薛東劍，鄭潔，李成繞，等.利用L波段星載重復(fù)軌道干涉SAR提取DEM及大氣效應(yīng)分析[J].測繪工程，2018，27(1)：5-9，14.

10.19349/j.cnki.issn1006-7949.2018.01.002

TP75

A

1006-7949(2018)01-0005-05

李銘娜]