基于精細DEM的InSAR大氣相位改正試驗研究

李 曼，夏 耶，2，葛大慶，張 玲，范景輝，王 艷

(1．中國國土資源航空物探遙感中心，北京 100083;2．德國地球科學研究中心，波茨坦 14473)

0 引言

我國是世界上滑坡災害最嚴重的國家之一。特別是庫區滑坡，滑坡變形破壞直接危害坡體上人員的生命財產安全;滑坡體入江減少水庫的有效庫容，形成堰塞湖并堵塞河道，同時滑坡體高速失穩入江后所激發的次生涌浪對航運、沿岸居民或建筑物也造成威脅［1］。因此，庫區滑坡的變形破壞及穩定性問題已成為制約水電開發的主要工程地質問題之一。開展庫區滑坡長期監測逐漸成為防止和減輕庫岸滑坡災害的重要手段之一。傳統上對滑坡的監測和預警主要是采用位移計、應力計、水位計等現場儀器觀測，或利用GPS和TDR等遠程電子設備監測，但這些技術方法很難應用到植被茂盛、人員難以進入地區的大面積滑坡的監測上［2］，且因監測成本高、監測網絡密度低，對變形緩慢的滑坡監測靈敏度不是很理想。

自20世紀90年代中期開始，合成孔徑雷達干涉測量(InSAR)技術以其監測成本低、監測網絡密度高且可以監測到整個滑坡體的形變全貌等獨特優勢，而被應用于對緩變形滑坡體形變的監測中。但由于早期的星載雷達數據分辨率多在20～30 m之間，重復周期一般為20～50 d，對于面積不大的滑坡體，其時間去相關現象非常嚴重;另一方面，大氣波動、尤其是對流層濕度和溫度的變化造成大氣延遲在差分相位圖上表現出非均一性，過大的大氣延遲相位甚至會“淹沒”研究區域的形變相位，使得雷達干涉技術在滑坡監測領域遠遠沒有取得實際應用的價值。隨著TerraSAR－X和Cosmo－SkyMed等星載高分辨率雷達數據的出現，以及相應高分辨率數字高程模型數據的成功獲取，短時間跨距雷達數據的時間去相關現象已迎刃而解［3］，而大氣波動延遲成為雷達干涉技術監測緩變形滑坡短時間跨距形變場的主要制約因素。如何從干涉相位或差分相位中抑制或者分離出大氣延遲相位以提高滑坡體待解算形變量的準確性，成為監測緩變形滑坡亟待解決的關鍵問題。

目前，在差分干涉合成孔徑雷達(D－InSAR)技術的應用過程中，消除或抑制大氣延遲效應的技術大致可以分為2類:①在數據處理過程中，引入外部大氣數據，如實測地表氣象參數、大氣數值模型、地基GPS水汽反演結果以及MODIS水汽反演結果等，這些數據都能提供大氣層中水汽含量信息，將水汽含量轉換為大氣相位，進而從差分干涉相位中消除［4－5］;②充分考慮大氣在時間域呈高頻、在空間域呈相對低頻的特點，對多時相SAR數據進行處理，采用SAR數據集進行自身大氣校正，進而達到抑制或分離、消除大氣干擾相位的目的。線性組合法、干涉圖疊加法、隨機濾波法以及長時序D－InSAR等方法在實際數據處理過程中經常被采用，其中永久散射體干涉測量技術(permanent scatterer interferometry，PSInSAR)［6－7］和短基線集技術(small baseline subset，SBAS)是長時序 D－InSAR技術的典型代表，也是目前最為先進的SAR數據集自身校正法。可見，現有消除SAR圖像中大氣延遲的技術方法都能在一定程度上削弱大氣擾動引起的相位貢獻，但它們又都有各自的適用條件和自身無法克服的缺點，在實際應用過程中都受到這樣或那樣不同程度的限制，如MODIS和MERIS的水汽反演數據不能在夜晚以及白天云量較多的情況下發揮作用;線性組合法和干涉圖疊加法是假設地表形變符合線性模型;隨機濾波法需要先確定相鄰的穩定區域等［4，7］。

但在地形起伏較大、降水豐沛的高山峽谷區用InSAR技術進行地表形變監測時，如果無法獲得外部大氣數據，而雷達數據量又較少的情況下，采用大氣濕延遲相位與地面高程之間的相關模型削弱大氣擾動相位的方法逐漸得到了國內外學者的關注。Delacourt等［8］在應用InSAR技術監測意大利Etna火山時指出，大氣對流層中濕延遲與大氣層的溫度和水蒸汽的分壓之間滿足一定的函數關系。Xia等［9］分析大氣效應引起的延遲相位依賴于研究區的氣候和高程，如果2次獲取雷達圖像時的氣溫基本相同，則由這2幅圖像得到的差分干涉圖中大氣延遲效應產生的相位很小;反之，如果2次獲取雷達圖像時的研究區溫度差別較大，則干涉圖中大氣相位與研究區的高程相關。蔣延臣［10］基于濕延遲與高程的線性回歸模型校正于田地震ScanSAR差分干涉紋圖的大氣效應，在一定程度上改善了干涉紋圖的質量，但在DEM梯度變化較大的地區，大氣擾動的去除效果仍不是很好。基于此，本文以三峽庫區秭歸縣樹坪滑坡區為研究對象，根據大氣濕延遲相位與地面高程的相關性，通過提取樹坪滑坡體以外區域高相干點上的高程值及解纏后的差分干涉相位，從數學角度尋求差分干涉相位(大氣濕延遲相位)與高程之間的最優函數模型。根據這一模型，由樹坪滑坡體的DEM得出該區域的大氣效應相位，將其從干涉圖中去除，最后保留的相位主要是由地表形變引起。這種改正大氣效應的方法，可有效地改善干涉圖的質量，提高雷達干涉測量技術的形變監測精度。

1 數據源

1．1 高分SAR數據

本文采用德國DLR的Terra SAR－X聚束模式的高空間分辨率的雷達數據進行差分干涉處理。Terra-SAR衛星重訪周期為11 d，載波使用X波段，波長3．1 cm。該衛星定位精度高，空間基線控制能力強，為利用傳統InSAR方法測定滑坡體短期內的小尺度形變量和形變場的非均勻分布提供了重要保障［11］。

1．2 高程數據

高程數據通常用于補償InSAR數據處理過程中的高程相位。本文所用DEM數據是美國萊卡公司生產的機載激光雷達系統ALS50－II獲取的，經過了數據解算、系統誤差檢校、航帶系統誤差檢查修正及激光點云分類、質量控制及DEM精度評價等處理，其高程絕對精度一般為10～20 cm，平面精度通常約0．5～1 m。由于客觀條件限制，三峽庫區大部分地形起伏很大，其平面精度會相應略有變化，但高程精度通常變化較小［12］。

圖1 時間跨距11 d的差分干涉相位圖Fig．1 Differential interferometric phase map with 11 d interval

2 數據處理及相位特性分析

2．1 InSAR數據處理

采用兩軌法提取雷達目標的相位信息，獲取地表形變場。將重軌獲取的2景雷達圖像精確配準并生成干涉圖后，利用外部高精度DEM數據和雷達衛星成像的幾何信息模擬生成高程相位，從干涉相位中減去平地相位、高程相位，此時得到的差分干涉紋圖(圖1)中任一像元x的相位φdiff(x)可視為由形變相位Δφd(x)、大氣延遲相位Δφatm(x)、因DEM數據不準確引起的殘余地形相位Δφε－h(x)、因空間基線不準確引起的殘差相位Δφb(x)以及系統熱噪聲等引起的誤差相位 Δw(x)組成［7，13］，即

式中W{}表示纏繞相位。

2．2 差分干涉相位特性分析

差分干涉紋圖中各類型相位的量級大小以及在干涉圖中的分布特征與研究區的氣候條件、選取數據的性能指標以及InSAR技術的處理方法等密切相關。

本研究選取的TerraSAR雷達衛星分辨率為1 m，其基線精度高(通常小于10 cm)，影像覆蓋范圍為10 km×5 km。因此，由空間基線不準確產生的殘差相位Δφb(x)會很小。

因DEM數據不準確引起的殘余地形相位Δφε－h(x)也很小。這是因為，研究中采用了1 m分辨率的DEM數據對高程相位進行了補償。假定DEM數據誤差為1 m，根據雷達微波的基本參數以及SAR干涉測量成像的幾何關系，通過

分別計算出圖1兩幅相位圖中心點處由高程誤差引起的干涉相位值分別為－0．03 rad(圖1(a))和0．07 rad(圖 1(b))。式(2)中:Δφε－h為殘余地形相位;B⊥為2次獲取數據時雷達衛星所處位置的垂直基線長度;RL為衛星與圖像中心點之間的距離;H為衛星與星下點之間的距離;Δε為高程誤差;θ0為衛星入射角;λ為雷達波的波長。

通過計算發現，高程誤差相位在整幅干涉圖中的值基本不變。說明差分干涉圖中的相位主要并不是由高程誤差引起的。這一點也可以從另一角度進一步印證:由式(2)可知，高程誤差引起的相位與垂直基線B⊥長度呈正比，如果差分干涉圖中的相位主要由高程誤差引起，那么垂直基線B⊥=－84．96 m的差分干涉圖(圖1(b))中的相位值約為垂直基線B⊥=35．20 m的差分干涉圖(圖1(a))中相位值的2．41倍，但圖1(a)(b)兩幅差分干涉紋圖中的相位值并不存在這種倍數關系。于是，可以肯定地說，在差分干涉相位圖中，高程誤差及空間基線不準確引起的相位值很小;同時，系統熱噪聲和數據處理誤差產生的相位雖然無法避免，但也相對較小，基本可以忽略不計。這樣，圖1差分干涉圖中的相位主要由地表形變和大氣擾動引起，即

對樹坪滑坡區地表形變的長期監測結果表明，只在樹坪滑坡體的局部區域存在緩慢形變跡象，其他區域都相對比較穩定。因此，圖1(a)中大范圍明顯的相位基本都是由大氣延遲引起的，樹坪滑坡區的形變相位完全被大氣延遲相位所“淹沒”。但由于雷達影像沿距離向覆蓋范圍僅10 km，雷達衛星入射角的變化也較小(約0．6°)，入射角變化引起的大氣擾動相位也不明顯，干涉圖中這一明顯的大氣相位基本是由大氣層本身的局部相關波動引起。在距離雷達衛星較遠區域(圖1(a)右側)，相位值基本不變，這可能與相應區域的大氣層厚度變化較小或者基本沒有變化有關;在距雷達衛星較近區域(圖1(a)左側)，大氣延遲相位梯度變化明顯，且與入射角無關，隨著高程的增加，大氣層厚度變化明顯。

3 大氣延遲相位的改正方法

以樹坪滑坡體附近穩定點CR08(見第4節圖2)為相位解纏基準點，采用最小費用流算法對差分干涉圖進行相位解纏。解纏后的相位仍是由大氣相位和形變相位組成，只有將大氣相位去除后，才能得到研究區沿雷達視線向的地表形變信息。

由于大氣厚度大、成分復雜且垂直成層分布，根據大氣對重復軌道干涉測量的影響，可以將地表和遙感器之間的大氣分為電離層、中間大氣(平流層和中間層)和對流層3部分，其中，對流層中的空氣對流很明顯，云、霧、雨等現象都發生在這一層，水蒸氣也幾乎都在這一層內。對流層又可以細分為干延遲和濕延遲2部分，干延遲由非水蒸氣部分的大氣延遲產生，在時域內比較穩定，在空域內具有大尺度變化的特性;濕延遲是由云、雨等水蒸氣引起，水蒸汽分布極不均勻，且與溫度密切相關，其變化在時間域呈高頻，在空間域呈相對低頻［11］。

根據大氣各組分在時域、空域中的變化特性不同，雷達影像經干涉、差分及相位解纏后，差分干涉相位中的大氣延遲相位主要由對流層中的濕延遲引起，其他層對形變相位的影響都已經基本消除。對流層中的濕延遲取決于研究區的溫度和水蒸氣壓力，而溫度、水蒸氣壓力與所處區域的高程密切相關［8］。因此，濕延遲相位與相應區域高程呈間接相關，依此可以有效消除解纏相位中大氣濕延遲引起的干擾相位。另外，研究區山體的坡度、坡向相對高程來說，對濕延遲的影響相對很小，且擾動程度和影響范圍也不盡相同，因此本文不考慮坡度及坡向對濕延遲相位的影響。

大氣濕延遲在空間上的相關距離相對較小，在整幅差分干涉圖中(圖1(a))的分布特征完全不同，左側大氣相位梯度大，且隨高程升高而呈規律性變化;右側相對較小，甚至與位置也存在一定的相關性。于是根據不同區域大氣相位的變化特征不同，分別提取相應區域內的濕延遲相位、高程信息和沿距離向的坐標，基于最小二乘法，對濕延遲相位與高程、沿距離向的坐標值進行曲線或曲面擬合，在保證擬合曲線或曲面在高程(0，800)m范圍內沒有明顯的轉折，殘差都處于［－2，2］rad之間，且在誤差平方和最小的前提下，建立相應區域的大氣效應校正模型。最后，將模擬的大氣相位從解纏相位中去除，即可得出樹坪滑坡沿雷達視線向的形變信息。

4 樹坪滑坡形變信息提取

樹坪滑坡的變形歷史較長，1996年以前，坡體局部就出現過變形，從2003年6月三峽水庫蓄水后，滑坡變形加劇，在滑坡體的上、中、下部都出現了明顯的地面裂縫和房屋開裂現象。目前，樹坪滑坡體局部區域仍有比較明顯的變形，但變形速率相對較小，每年僅數10 cm，尤其在三峽水庫蓄水及水位消落期間，再加上大氣降水的影響，滑坡體局部區域變形速率較大［14］。本文采用高分辨率TerraSAR－X數據，通過傳統的雷達干涉技術提取樹坪滑坡形變信息。

圖2是樹坪滑坡在2012年1月份、時間跨距為11 d的差分干涉相位圖。由于在2012年1月2日至1月13日期間，樹坪滑坡區以小雨或陰天為主，且氣溫變化明顯，因而圖2顯示大氣濕延遲引起的擾動相位較大，甚至掩蓋了樹坪滑坡體的形變信息。在圖2中Ⅰ區域，大氣濕延遲相位僅與該區的高程密切相關，隨著高程的升高，大氣濕延遲相位變化明顯;而在Ⅱ區域及樹坪滑坡體區域，濕延遲相位不僅與研究區的高程有關，還與雷達與目標間的距離(即沿距離向的坐標值)有關，這就需要分別建立大氣濕延遲效應校正模型。

圖2 樹坪滑坡時間跨距為11 d的差分干涉相位圖(20120102—20120113)Fig．2 Differential interferometric phase map of Shuping landslide with 11 d interval

在圖2Ⅰ區域內，由相干系數≥0．3點上的濕延遲相位和高程值得出該區域的最優函數校正模型，即根據Ⅰ區域的DEM數據，由式(4)即可模擬出該區域所有點的大氣濕延遲相位，且大氣殘差相位都介于［－1．5，1．0］之間。

但在圖2的II區域及樹坪滑坡區，濕延遲相位不僅是高程的二次函數，同時也是圖像距離向坐標X的線性函數，且濕延遲在這2個區域的分布特征基本一致，由Ⅱ區域高相干點(相干系數≥0．3)上的濕延遲相位、高程值HⅡ及距離向坐標X，通過曲面擬合，模擬出Ⅱ區域的最優函數校正模型為式(5)，即

樹坪滑坡區的濕延遲相位分布特征也遵循式(5)的表達規律。因此，基于Ⅱ區域及樹坪滑坡區內所有點上的DEM值及沿距離向的坐標值X，由式(5)就可模擬出Ⅱ區域及樹坪滑坡體上所有點上的大氣濕延遲相位。濕延遲相位的殘差相位Δφresi，Ⅱwet介于［－2．0，1．0］之間。在對圖2差分干涉相位解纏時，解纏相位參考點CR08位于Ⅰ和Ⅱ區域的共同區域A中，解纏后A區域內濕延遲相位值接近于0。通過濕延遲相位校正模型(4)(5)，分別得出A區域的模擬大氣相位值差很小，取2次模擬相位的平均值作為A區域的濕延遲相位，即可得到圖2整個區域平滑的濕延遲相位分布狀態，將其從解纏相位中去除，進而恢復出了樹坪滑坡體在2012年1月2日至1月13日之間11 d內的形變信息(圖3)。

圖3 樹坪滑坡11 d內的形變信息(mm)Fig．3 Deformation of Shuping landslide with 11 d interval(mm)

從圖3可以看出，盡管樹坪滑坡周圍某些局部區域大氣濕延遲相位并沒有完全去除，但圖中紅、黃區域仍十分清楚地框出了樹坪滑坡體形變場在11 d中的位置、大小及其分布，且形變場是非均勻的。紅色區域為最大下沉區，集中分布在樹坪滑坡上部，沿雷達視線向最大沉降值達－4．4mm;黃色區域的形變值相對較小，在滑坡體中部及坡體前緣都有分布，樹坪滑坡這種局部變形特性與所處地區的氣候及三峽水庫水位消落密切相關。2012年1月初，秭歸小雨連綿不斷，降雨對滑坡的誘發作用從滑坡體上部開始，由于坡體后緣存在較寬裂縫，雨水入滲使滑帶巖土體發生軟化、泥化作用，造成滑帶的抗剪強度降低。另一方面，雨水在坡體內積聚形成較大的靜水壓力，削弱滑動面上的有效正應力，導致坡體抗滑力減小，這都會造成滑坡體上部局部區域出現較大變形。三峽水庫水位消落直接誘發滑坡體前緣出現明顯形變，三峽水庫于2011年12月31日開始啟動新一輪生態補水，2012年1—2月，三峽水庫下泄量按6 000 m3/s左右控制。可見，在這一時期三峽水庫水位一直在持續消落，而坡體組成物質的透水性弱，坡體內孔隙水壓力消散速度滯后于庫水位的降低速度，使得坡體前緣動水壓力增大，進而導致滑坡體的下滑力增加。樹坪滑坡體局部區域雖有宏觀變形現象，但其內并沒有形成貫通的滑動面，目前滑坡整體仍處于穩定狀態。

5 結論

基于樹坪滑坡區特殊的地理位置、復雜的地形地貌特征以及多雨潮濕的氣候條件，本文通過分析大氣濕延遲在空間的分布特征，分別建立空間局部區域大氣效應的最優函數校正模型，將大氣濕延遲相位從解纏相位中去除，最終恢復出樹坪滑坡體的形變場。

需要說明的是，這種大氣改正方法并沒有考慮研究區山體的坡向和坡度角變化對濕延遲產生的影響;同時又假設了對流層中濕延遲相位在空間局部區域內具有相關性。任何函數模型都不能完全確切地描述濕延遲相位的空間分布規律，無法準確無誤地將差分干涉圖中的濕延遲相位完全去除，在樹坪滑坡體的形變場中仍有濕延遲效應的影響，但由于樹坪滑坡體局部區域變形較大，通過這種去除大氣效應的方法，清楚地框出了樹坪滑坡體形變場在11 d中的位置、大小及其分布，較為真實地恢復出了樹坪滑坡體的形變信息，這有利于推動D－InSAR技術在三峽庫區大面積緩變形滑坡形變監測及穩定性調查的廣泛應用。

志謝:感謝中國國土資源航空物探遙感中心鄭雄偉提供了高精度DEM數據精度評價資料。

［1］ 李遠耀．三峽庫區漸進式庫岸滑坡的預測預報研究［D］．武漢:中國地質大學，2010:1－30．Li Y Y．Research on prediction and forecast of progressive bank landslide in the Three Gorges Reservoir［D］．Wuhan:China University of Geosciences，2010:1－30．

［2］ 王桂杰，謝謨文，邱 聘，等．D－InSAR技術在大范圍滑坡監測中的應用［J］．巖土力學，2010，31(4):1337－1344．Wang G J，Xie M W，Qiu P，et al．Application of D－InSAR technique to landslide monitoring［J］．Rock and Soil Mechanics，2010，31(4):1337－1344．

［3］ 夏 耶，范景輝，李 曼，等．高分辨率雷達數據三峽庫區滑坡監測技術［C］//第十八屆中國遙感大會論文集，2012:161－169．Xia Y，Fan JH，LiM，etal．Monitoring technique of Shuping landslide with high－resolution Radar images in Three Gorges Reservoir area，in China［C］//18th China Symposium on Remote Sensing，2012:161－169．

［4］ 吳云孫，李振洪，劉經南，等．InSAR觀測值大氣改正方法的研究進展［J］．武漢大學學報:信息科學版，2006，31(10):862－867．Wu Y S，Li ZH，Liu JN，et al．Atmospheric correction models for InSAR measurements［J］．Geomatics and Information Science of Wuhan University，2006，31(10):862－867．

［5］ 范青松，湯翠蓮，陳 于，等．GPS與InSAR技術在滑坡監測中的應用研究［J］．測繪科學，2006，31(5):60－62．Fan Q S，Tang C L，Chen Y，et al．Applications of GPS and InSAR in monitoring of landslide studies［J］．Science of Surveying and Mapping，2006，31(5):60－62．

［6］ Ferretti A，Prati C，Rocca F．Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry［J］．IEEE Transactions on Geoscience and Remote Sensing，2000，38(5):2202－2212．

［7］ 范景輝．基于相干目標的D－InSAR技術地表形變監測研究與應用［D］．北京:中國科學院研究生院，2008，136－139．Fan JH．Research on coherent targets based on D－InSAR technique and its application to surface deformation monitoring［D］．Beijing:Chinese Academy of Sciences，2008，136－139．

［8］ Delacourt C，Briole P，Achache J．Tropospheric corrections of SAR interferograms with strong topography，application to Etna［J］．Geophysical Research Letters，1998，25(15):2849－2852．

［9］ Xia Y，Kaufmann K，Guo X F．Landslide monitoring in the Three Gorges area using D－InSAR and corner reflectors［J］．Photogrammetric Engineering and Remote Sensing，2004，70(10):1167－1172．

［10］ 蔣延臣．星載寬幅合成孔徑雷達干涉測量形變監測理論與應用［D］．武漢:武漢大學，2010:82－90．Jiang Y C．Study on scanning synthetic aperture radar interferometry theory and application for deformation monitoring［D］．Wuhan:Wuhan University，2010:82－90．

［11］ 李小凡，李 穎，曾琪明，等．應用與ASAR同步的MERIS對重復軌道InSAR進行大氣校正［J］．北京大學學報:自然科學版，2009，45(6):1012－1018．Li X F，Li Y，Zeng Q M，et al．Correction of atmospheric effects on repeat－pass interferometric SAR using MERIS and ASAR synchronous data［J］．Acta Scientiarum Naturalium Universitatis Pekinensis，2009，45(6):1012－1018．

［12］ 劉圣偉，郭大海，陳偉濤，等．機載激光雷達技術在長江三峽工程庫區滑坡災害調查和監測中的應用研究［J］．中國地質，2012，39(2):507－517．Liu SW，Guo D H，Chen H T，et al．The application of airborne lidar technology in landslide investigation and monitoring of Three Gorges Reservoir area［J］．Geology in China，2012，39(2):507－517．

［13］ 葛大慶，王 艷，范景輝，等．地表形變D－InSAR監測方法及關鍵問題分析［J］．國土資源遙感，2007，19(4):14－22．Ge D Q，Wang Y，Fan JH，et al．A study of surface deformation monitoring using differential SAR interferometry technique and an analysis of its key problems［J］．Remote Sensing for Land and Resources，2007，19(4):14－22．

［14］ 汪發武，張業明，王功輝，等．三峽庫區樹坪滑坡受庫水位變化產生的變形特征［J］．巖石力學與工程學報，2007，26(3):509－517．Wang FW，Zhang Y M，Wang G H，et al．Deformation features of Shuping landslide caused by water level changes in Three Gorges Reservoir area，China［J］．Chinese Journal of Rock Mechanics and Engineering，2007，26(3):509－517．