基于網絡法時序InSAR 大氣誤差校正方法研究

李永生 張景發 姜文亮 羅 毅 王曉醉

1 中國地震局地殼應力研究所（地殼動力學重點實驗室），北京市安寧莊路1號，100085

2 中國地震局工程力學研究所，哈爾濱市學府路29號，150080

3 武漢大學科學技術發展研究院，武漢市珞喻路129號，430079

單幅干涉圖的測量精度受很多噪聲源的影響，其中之一是雷達信號在穿過大氣和電離層時所引起的相位延遲［1］。大氣相位延遲往往與地表形變信號纏繞在一起，湮沒真正的地面形變信號。因此，對大氣延遲信號的精確估計是高精度形變反演的重要內容［2］。目前去除InSAR 中大氣延遲主要依據地面氣象觀測資料［3－6］、GPS 數據［7－9］、MODIS或者MERIS數據等［10－13］，而更多的是發掘InSAR 數據集的內部聯系以及大氣信號的物理特性進行大氣誤差的校正。通常假設大氣延遲信號在時間維度上不相關，而在空間維度上相關，通過濾波方式對單個干涉圖進行處理［14－21］。本文采用短基線網絡解算方法，對每一幅影像獲取時刻的大氣延遲誤差進行估計，再模擬每一幅干涉圖中的大氣延遲相位，最后從解纏干涉圖上去除。根據延遲成因，將大氣誤差延遲分量分成3類，針對不同的分量類型，在短基線方法基礎上采用網絡法分別進行校正，并用一維協方差函數對大氣延遲誤差校正效果進行估計。

1 大氣延遲誤差

大氣相位延遲對干涉圖的影響包含兩部分［22］：1）與大氣三維異質分布有關，與地形不相關。該影響對山區和平原地區相似，可以利用二階平穩過程來描述，可細分為長波長誤差延遲和短波長誤差延遲（也稱湍流大氣延遲誤差）。2）與大氣的垂直特性相關。可以假設相位和地形呈線性關系［23］，并利用該假設對垂直大氣相位延遲進行校正。

如果大氣延遲在SAR 覆蓋范圍內分布均勻，可以認為其相位屬于長波長信號。該相位延遲在干涉圖上常常表現為一階或二階相位斜面。在實際處理中，長波長大氣相位誤差延遲在表現形式上和軌道誤差相似，可以當成軌道誤差信號，在軌道誤差處理時大部分可以同時消減。

湍流大氣主要發生在大氣中水汽變化劇烈的地區。水汽在垂直和橫向上短距離內發生劇烈變化，導致干涉圖產生隨機的干涉紋理。湍流大氣延遲誤差可以看成由各種不同尺度的大氣渦旋疊合而成的流動，這些漩渦的大小及旋轉軸方向呈隨機分布，也使得大氣延遲誤差呈隨機分布。該類延遲誤差是干涉處理中最為常見的形式。

大氣垂直分層延遲相位（topography correcting atmospheric delay，TCAD）由主從影像獲取時刻大氣層最低部的平均水汽含量變化所引起［22，24］。一般認為，空氣的溫度和壓力在空間上垂直分層。在平坦地區，垂直分層中大氣水汽含量在空間上相對均勻。如果在影像上出現高程變化，濕延遲隨高程以一定的比率變化。因此，雷達影像上的高程變化會導致大的相位延遲變化［23，25－28］。

2 網絡法大氣誤差校正

傳統的大氣延遲相位估計是基于單個干涉圖進行的［14－16］。本文將利用網絡法先對每一幅影像獲取時刻的大氣延遲誤差進行估計，再利用這些估計得到的單個時刻的延遲誤差重新構建每一干涉圖的大氣延遲誤差。假設δφj（x，y）為第j幅干 涉圖上（x，y）處 的 相 位 值，φ（tA，x，y）和φ（tB，x，y）分別表示在tA和tB成像時刻（x，y）處的相位值，每一幅干涉圖可以表示為：

基于短基線集網絡可以構建如下方程組：

式中，A表示M×N矩陣。矩陣A的元素Akl按照如下規則定義：如果l＝tB，則Akl＝1；如果l＝tA，則Akl＝－1；否則Akl＝0。tA和tB表示干涉圖k中主影像和從影像的成像時間。δφ是M維已知矢量，表示M幅干涉圖；φ是N維未知矢量，表示N個成像時刻的相位值。由于上述方程組中矩陣A秩虧，所以無法得到唯一解，可采用奇異值分解方法解算。解算式（2），可以得到每一時刻的大氣延遲，進而利用式（1）模擬出每一幅干涉圖的相位值。

對上述3種大氣延遲分量分別進行估計，其中每幅干涉圖中主從影像i和j上的長波長大氣相位延遲平面可以分別表示為φflati＝aix＋biy＋ci和φflatj＝ajx＋bjy＋cj，則干涉圖上大氣延遲相位平面可以表示為：

大氣垂直分層延遲相位（TCAD）可以表示為兩個獲取時刻的相位／高程的線性關系坡度差，影像i和j可以分別表示為φtropoi＝kiZ和φtropoj＝kjZ。其中k為未知常數，Z表示高程。構建下面的公式：

主從影像獲取時刻的對流層中湍流大氣延遲相位分別描述為φturbui和φturbuj，則對應的干涉圖湍流大氣延遲相位可以表示為φturbuij＝φturbui－φturbuj。

3 實驗結果

本文使用20 景ASAR 降軌數據，軌道為405，時間范圍為2007－02－20～2010－09－07。基于短基線集網絡構建原則，共生成60個干涉對，見圖1。該地區屬于多山地形，天氣變化劇烈，所以在一幅干涉圖上有可能由一個或者多個大氣因素引起大氣延遲誤差。

圖1 短基線集干涉圖組合方式Fig.1 Temporal and perpendicular baselines for the interferograms used in this study

一般的，長波長相位誤差包含軌道誤差、長波長大氣相位延遲誤差以及地殼長波長活動信號等。這里暫不考慮地殼長波長活動信號，采用網絡校正法同時對前兩種誤差進行去除。圖2顯示了長波長相位誤差去除后的估計結果。圖2（a）顯示了原始解纏干涉圖，從干涉圖上可以看出，在NW－SE向存在一個明顯的相位坡度，該坡度主要由軌道誤差和長波長大氣相位誤差造成。圖2（b）為經過長波長誤差相位校正后的解纏干涉圖。在圖2（a）、（b）相同位置分別選取一條剖面線進行對比可以發現，校正長波長相位延遲誤差后，NWSE向的相位坡度得到顯著改善，見圖2（c）。

與地形相關的大氣延遲誤差校正結果如圖3。圖3（a）顯示了校正前的解纏干涉圖，在河谷、湖泊邊緣與高山區域出現了明顯的與地形起伏相關的大氣延遲相位信號；圖3（b）為校正后的解纏干涉圖；圖3（c1）、（c2）分別表示地形相關大氣延遲誤差校正前后，LOS向相位和高程變化的線性關系。從圖3（c1）可以看出，該干涉圖中LOS向相位變化存在與高程線性相關的趨勢。在相位誤差校正后，線性坡度關系已經得到明顯修正，見圖3（c2）。

采用網絡法對湍流大氣延遲誤差進行模擬和校正，見圖4。圖4（a）中，經過長波長相位誤差校正之后，大氣延遲誤差主要由湍流大氣引起。對構成干涉的主從影像進行湍流大氣延遲模擬，并重構干涉圖的湍流大氣延遲誤差，見圖4（b）。與圖4（a）對比可以發現，網絡法可以精確模擬出干涉圖中主要的大氣湍流延遲誤差。圖4（c）為大氣校正后的解纏干涉圖，可以看出干涉圖上的大部分湍流延遲引起的相位誤差都得到了校正。

圖2 干涉圖（2007－07－10～2007－08－14）中長波長相位延遲校正Fig.2 The example of long wavelengthartifacts correction for the interferogram

圖3 干涉圖（2007－05－01～2007－08－14）中地形相關大氣延遲誤差校正Fig.3 The elimination example of topography correlated atmospheric delay（TCAD）for interferogram

圖4 干涉圖（2007－07－10～2007－08－14）中湍流大氣延遲誤差校正Fig.4 The example of atmospheric turbulence correction for the interferogram

為了估計干涉圖大氣相位校正前后的誤差水平，本文利用了一維協方差函數估計［29－30］：

其中，Cij表示像元i和j的協方差，σ2表示方差，dij表示兩個像元之間的距離，α為一維協方差函數的e－folding波長。

圖5顯示了所有解纏干涉圖在大氣延遲誤差校正前和校正后的一維協方差函數。圖5（a）中的紅色線表示校正前，藍色表示校正之后。在大氣誤差校正前，σ2和α的平均值分別為3.1mm2和1.5 km，見圖5（b）、（c）。在大氣誤差校正之后，兩者的平均值分別降低到了0.6mm2和0.21km，見圖5（d）、（e）。平均方差和e－folding波長分別降低了80%和86%，說明網絡法能有效地消減干涉圖中的大氣延遲誤差貢獻。

圖5 （a）大氣延遲誤差校正前后的指數協方差函數對比圖；（b）、（d）表示指數函數中方差在大氣誤差校正前后的直方圖顯示；（c）、（e）表示指數函數中e－folding的波長在大氣誤差校正前后的直方圖顯示Fig.5 （a）The comparison diagram of exponential covariance functions of the interferograms before and after APS correction.（b）and（d）Histogram of variance of the exponential function before and after the APS correction；（c）and（e）Histogram of e－folding wavelength of the exponential function before and after the APS correction

4 結 語

本文將大氣延遲相位誤差分量分成3類：長波長大氣相位信號、地形相關的大氣延遲誤差信號以及湍流大氣延遲誤差信號，并采用網絡法分別對3類大氣延遲誤差進行估計。從實驗結果分析，采用網絡法能精確地模擬出每種類型的大氣延遲誤差。采用一維協方差函數對大氣延遲誤差進行估計，平均方差經過大氣延遲校正從原來的3.1mm2降低到0.6mm2，降低了80%；而e－folding波長從原來的1.5km 減低到0.21km，減低了86%。

目前采用該方法還需要注意以下幾點：

1）在很多干涉圖上，3種大氣延遲誤差并不同時存在，常常只存在一種或者兩種，所以過度校正會引入新的誤差。

2）在對單幅干涉圖進行孤立的長波長信號處理時，會將由地殼活動導致的長波長信號當成軌道誤差或長波長大氣延遲誤差并被去除。一般的，時間基線較短的干涉對上，地殼活動導致的長波長信號可以忽略，但在網絡法計算中可以作為一種約束條件，改善長波長誤差相位的估計精度［21］。

3）針對地形相關大氣延遲誤差，采用相位和高程的全局線性回歸策略進行校正。該方法適用于局部地區（小范圍區域內）的TCAD 校正，但在全局區域（整景數據）該關系并不一定滿足。下一步將采用小波變換算法，對斜距相位變化和DEM進行小波多分辨率分解，進而進行地形相關大氣延遲誤差校正［22］。

致謝：感謝國家科技部、歐洲空間局龍計劃三期提供的ENVISAT ASAR原始數據。

［1］Williams S，Bock Y，Fang P.Integrated Satellite Interferometry：Tropospheric Noise，GPS Estimates and Implications for Interferometric Synthetic Aperture Radar Products［J］.Geophys Res，1998，103（B11）：51－67

［2］Barnhart W D，Lohman R B.Characterizing and Estimating Noise in InSAR and InSAR Time Series with MODIS［J］.Geochem Geophys Geosyst，2013，14（10）：4 121－4 132

［3］Zebker H A，Rosen P A，Hensley S.Atmospheric Effects in Interferometric Synthetic Aperture Radar Surface Deformation and Topographic Maps［J］.Geophys Res，1997，102：7 547－7 563

［4］Delacourt C，Achache B P.Tropospheric Corrections of SAR Interferograms with Strong Topography：Application to Etna［J］.Geophys Res Lett，1998，25：2 849－2 852

［5］Bonforte A，Ferretti A，Prati C，et al.Calibration of Atmospheric Effects on SAR Interferograms by GPS and Local Atmosphere Models：First Results［J］.Atmos J Terr Phys，2001，63：1 343－1 357

［6］Li Z W.Comparative Study of Empirical Tropospheric Models for the Hong Kong Region［J］.Survey Review，2008，40（310）：328－341

［7］Onn F，Zebker H A.Correction for Interferometric Synthetic Aperture Radar Atmospheric Phase Artifacts Using Time Series of Zenith Wet Delay Observations from a GOS Network［J］.Geophys Res，2006，111：B09102

［8］Onn F.Modeling Water Vapor Using GPS with Application to Mitigating InSAR Atmospheric Distortions［D］.Stanford University，2006

［9］Li Z H，Fielding E J，Cross P，et al.Interferometric Synthetic Aperture Radar Atmospheric Correction：GPS Topography－Dependent Turbulence Model［J］.Geophys Res，2006，111：B02404

［10］Gao B C，Kaufman Y J.Water Vapor Retrievals Using Moderate Resolution Imaging Spectroradiometer（MODIS）Near Infrared Channels［J］.Geophys Res，2003，108：4 389－4 398

［11］Li Z H，Muller J P，Cross P.Tropospheric Correction Techniques in Repeat－pass SAR Interferometry［C］.FRINGE 2003 Workshop，ESA ESRIN，Frascati，Italy，2003

［12］Li Z H，Muller J P，Cross P，et al.Interferometric Synthetic Aperture Radar（InSAR）Atmospheric Correction：GPS，Moderate Resolution Imaging Spectroradiometer（MODIS），and In－SAR Integration［J］.Geophys Res，2005，110：B03410

［13］Li Z H，Muller J P，Cross P，et al.Assessment of the Potential of MERIS Near－infrared Water Vapour Products to Correct ASAR Interferometric Measurements［J］.Int J Remote Sens，2006，27：349－365

［14］Ferretti A，Prati C，Rocca F.Nonlinear Subsidence Rate Estimation Using Permanent Scatters in Differential SAR Interferometry［J］.IEEE T Geosci Remote Sens，2000，38：2 202－2 212

［15］Ferretti A，Prati C，Rocca F.Permanent Scatters InSAR Interferometry［J］.IEEE T Geosci Remote Sens，2001，39：8－20

［16］Hooper A，Zebker H，Segall P，et al.A New Method for Measuring Deformation on Volcanoes and Other Natural Terrains Using InSAR Persistent Scatterer［J］.Geophys.Res Lett，2004，31：L23611

［17］Berardino P，Fornaro G，Lanari R，et al.A New Algorithm for Surface Deformation Monitoring Based on Small Baseline Differential SAR Interferograms［J］.IEEE Transactions on Geoscience and Remote Sensing，2002，40：2 375－2 383

［18］Schmidt D A，Bürgmann R.Time Dependent Land Uplift and Subsidence in the Santa Clara Valley，California，from a Large Interferometric Synthetic Aperture Radar Data Set［J］.Journal of Geophysical Research，2003，108（B9）：2 416

［19］Mora O，Mallorqui J J，Broquetas A.Linear and Nonlinear Terrain Deformation Maps from a Reduced Set of Interferometric SAR Images［J］.IEEE Transactions on Geoscience and Remote Sensing，2003，41：2 243－2 253

［20］Prati C，Ferretti A，Perissin D.Recent Advances on Surface Ground Deformation Measurement by Means of Repeated Space Borne SAR Observations［J］.Journal of Geodynamics，2010，49：161－170

［21］Juliet B.Multi－interferogram Method for Measuring Interseismic Deformation：Denali Fault，Alaska［J］.Geophysical Journal International，2007，170（3）：1 165－1 179

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

［23］CavaliéO.Ground Motion Measurement in the Lake Mead Area，Nevada，by Differential Synthetic Aperture Radar Interferometry Time Series Analysis：Probing the Lithosphere Rheological Structure［J］.Journal of Geophysical Research，2007，112（B3）：B03403

［24］Lin Y N.A Multiscale Approach to Estimating Topographically Correlated Propagation Delays in Radar Interferograms［J］.Geochemistry Geophysics Geosystems，2010，11（9）：Q09002

［25］Shirzaei M，Bürgmann R.Topography Correlated Atmospheric Delay Correction in Radar Interferometry Using Wavelet Transforms［J］.Geophys Res Lett，2012，39：L01305

［26］Romain J.Systematic InSAR Tropospheric Phase Delay Corrections from Global Meteorological Reanalysis Data［J］.Geophysical Research Letters，2011，38（17）：L17311

［27］Doin M P.Corrections of Stratified Tropospheric Delays in SAR Interferometry：Validation with Global Atmospheric Models［J］.Journal of Applied Geophysics，2009，69（1）：35－50

［28］CavaliéO.Measurement of Interseismic Strain Across the Haiyuan Fault（Gansu，China），by InSAR［J］.Earth and Planetary Science Letters，2008，275（3）：246－257

［29］Parsons B.The 1994Sefidabeh（Eastern Iran）Earthquakes Revisited：New Evidence from Satellite Radar Interferometry and Carbonate Dating about the Growth of an Active Fold above a Blind Thrust Fault［J］.Geophysical Journal International，2006，164（1）：202－217

［30］Wang Hua.InSAR Reveals Coastal Subsidence in the Pearl River Delta，China［J］.Geophysical Journal International，2012，191（3）：1 119－1 128