基于網(wǎng)絡(luò)法時(shí)序InSAR 大氣誤差校正方法研究

李永生 張景發(fā) 姜文亮 羅 毅 王曉醉

1 中國(guó)地震局地殼應(yīng)力研究所（地殼動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室），北京市安寧莊路1號(hào)，100085

2 中國(guó)地震局工程力學(xué)研究所，哈爾濱市學(xué)府路29號(hào)，150080

3 武漢大學(xué)科學(xué)技術(shù)發(fā)展研究院，武漢市珞喻路129號(hào)，430079

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

1 大氣延遲誤差

大氣相位延遲對(duì)干涉圖的影響包含兩部分［22］：1）與大氣三維異質(zhì)分布有關(guān)，與地形不相關(guān)。該影響對(duì)山區(qū)和平原地區(qū)相似，可以利用二階平穩(wěn)過(guò)程來(lái)描述，可細(xì)分為長(zhǎng)波長(zhǎng)誤差延遲和短波長(zhǎng)誤差延遲（也稱(chēng)湍流大氣延遲誤差）。2）與大氣的垂直特性相關(guān)。可以假設(shè)相位和地形呈線(xiàn)性關(guān)系［23］，并利用該假設(shè)對(duì)垂直大氣相位延遲進(jìn)行校正。

如果大氣延遲在SAR 覆蓋范圍內(nèi)分布均勻，可以認(rèn)為其相位屬于長(zhǎng)波長(zhǎng)信號(hào)。該相位延遲在干涉圖上常常表現(xiàn)為一階或二階相位斜面。在實(shí)際處理中，長(zhǎng)波長(zhǎng)大氣相位誤差延遲在表現(xiàn)形式上和軌道誤差相似，可以當(dāng)成軌道誤差信號(hào)，在軌道誤差處理時(shí)大部分可以同時(shí)消減。

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

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

2 網(wǎng)絡(luò)法大氣誤差校正

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

基于短基線(xiàn)集網(wǎng)絡(luò)可以構(gòu)建如下方程組：

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

對(duì)上述3種大氣延遲分量分別進(jìn)行估計(jì)，其中每幅干涉圖中主從影像i和j上的長(zhǎng)波長(zhǎng)大氣相位延遲平面可以分別表示為φf(shuō)lati＝aix＋biy＋ci和φf(shuō)latj＝ajx＋bjy＋cj，則干涉圖上大氣延遲相位平面可以表示為：

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

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

3 實(shí)驗(yàn)結(jié)果

本文使用20 景ASAR 降軌數(shù)據(jù)，軌道為405，時(shí)間范圍為2007－02－20～2010－09－07。基于短基線(xiàn)集網(wǎng)絡(luò)構(gòu)建原則，共生成60個(gè)干涉對(duì)，見(jiàn)圖1。該地區(qū)屬于多山地形，天氣變化劇烈，所以在一幅干涉圖上有可能由一個(gè)或者多個(gè)大氣因素引起大氣延遲誤差。

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

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

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

采用網(wǎng)絡(luò)法對(duì)湍流大氣延遲誤差進(jìn)行模擬和校正，見(jiàn)圖4。圖4（a）中，經(jīng)過(guò)長(zhǎng)波長(zhǎng)相位誤差校正之后，大氣延遲誤差主要由湍流大氣引起。對(duì)構(gòu)成干涉的主從影像進(jìn)行湍流大氣延遲模擬，并重構(gòu)干涉圖的湍流大氣延遲誤差，見(jiàn)圖4（b）。與圖4（a）對(duì)比可以發(fā)現(xiàn)，網(wǎng)絡(luò)法可以精確模擬出干涉圖中主要的大氣湍流延遲誤差。圖4（c）為大氣校正后的解纏干涉圖，可以看出干涉圖上的大部分湍流延遲引起的相位誤差都得到了校正。

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

圖3 干涉圖（2007－05－01～2007－08－14）中地形相關(guān)大氣延遲誤差校正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

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

其中，Cij表示像元i和j的協(xié)方差，σ2表示方差，dij表示兩個(gè)像元之間的距離，α為一維協(xié)方差函數(shù)的e－folding波長(zhǎng)。

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

圖5 （a）大氣延遲誤差校正前后的指數(shù)協(xié)方差函數(shù)對(duì)比圖；（b）、（d）表示指數(shù)函數(shù)中方差在大氣誤差校正前后的直方圖顯示；（c）、（e）表示指數(shù)函數(shù)中e－folding的波長(zhǎng)在大氣誤差校正前后的直方圖顯示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 結(jié) 語(yǔ)

本文將大氣延遲相位誤差分量分成3類(lèi)：長(zhǎng)波長(zhǎng)大氣相位信號(hào)、地形相關(guān)的大氣延遲誤差信號(hào)以及湍流大氣延遲誤差信號(hào)，并采用網(wǎng)絡(luò)法分別對(duì)3類(lèi)大氣延遲誤差進(jìn)行估計(jì)。從實(shí)驗(yàn)結(jié)果分析，采用網(wǎng)絡(luò)法能精確地模擬出每種類(lèi)型的大氣延遲誤差。采用一維協(xié)方差函數(shù)對(duì)大氣延遲誤差進(jìn)行估計(jì)，平均方差經(jīng)過(guò)大氣延遲校正從原來(lái)的3.1mm2降低到0.6mm2，降低了80%；而e－folding波長(zhǎng)從原來(lái)的1.5km 減低到0.21km，減低了86%。

目前采用該方法還需要注意以下幾點(diǎn)：

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

2）在對(duì)單幅干涉圖進(jìn)行孤立的長(zhǎng)波長(zhǎng)信號(hào)處理時(shí)，會(huì)將由地殼活動(dòng)導(dǎo)致的長(zhǎng)波長(zhǎng)信號(hào)當(dāng)成軌道誤差或長(zhǎng)波長(zhǎng)大氣延遲誤差并被去除。一般的，時(shí)間基線(xiàn)較短的干涉對(duì)上，地殼活動(dòng)導(dǎo)致的長(zhǎng)波長(zhǎng)信號(hào)可以忽略，但在網(wǎng)絡(luò)法計(jì)算中可以作為一種約束條件，改善長(zhǎng)波長(zhǎng)誤差相位的估計(jì)精度［21］。

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

致謝：感謝國(guó)家科技部、歐洲空間局龍計(jì)劃三期提供的ENVISAT ASAR原始數(shù)據(jù)。

［1］Williams S，Bock Y，F(xiàn)ang 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，F(xiàn)erretti 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，F(xiàn)ielding 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，F(xiàn)rascati，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，F(xiàn)ornaro 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，F(xiàn)erretti 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