基于网络法时序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