时序合成孔径雷达干涉测量（interferometric synthetic aperture radar，InSAR）技术是测量地表形变的重要手段，可以用来获取大面积缓慢地表形变时变信息，在众多地球科学研究领域中得到了广泛应用^[1-2]。然而，时序InSAR技术的精度受多种误差因素影响，其中，大气延迟误差校正是时序InSAR技术面临的最大挑战^[3]。卫星发射的微波信号在穿过大气介质时，受大气水汽含量、温度和压力等条件的影响会发生不同程度的相位延迟，即大气延迟效应，通常包括电离层延迟和对流层延迟^[4]。由于C波段波长较短且研究区域位于中纬度，电离层延迟影响较小^[5]，本文仅研究对流层延迟。

对流层大气参数（如水汽含量、温度和压力）的时空变化可引起高达几十厘米量级的InSAR延迟信号^[6]，微小尺度的构造运动引起的地表位移可能会被掩盖。已有研究表明，阿尔金断裂滑动速率自西向东沿断层走向由15 mm/a逐渐减小至1 mm/a^[7]，绝大部分区域左旋走滑速率小于10 mm/a，表示10 a间隔的地表位移不足10 cm，投影到InSAR视线向（line of sight，LOS）则更小，比大气延迟低一个量级^[8]。

对流层延迟包含垂直分层延迟和湍流混合延迟两部分，大气折射率沿垂直方向的变化引起垂直分层延迟（本文将大气等效成无数个水平薄层，假设同一层折射率相同），而大气对流过程会引起湍流混合延迟（大气对流导致大气折射率在不同时刻的空间分布不同）。目前，国内外很多专家学者对大气对流层延迟改正方法展开研究，所得成果可以分为3类：（1）主要通过叠加、平均或时空滤波等方法去除具有时空随机分布性质的湍流混合延迟。文献[9]发现相位叠加技术能够有效降低InSAR数据的数字高程模型（digital elevation model，DEM）残差、大气相位及轨道误差的影响。根据大气信号在时间上高频、空间上低频的特性，时序InSAR通过时空滤波即可减少大气延迟的影响^[10]。（2）利用气象模型、多光谱观测或与其他外部数据相结合的方法估计对流层延迟。文献[11]利用动态的气象模型估计和去除大气水汽引起的信号延迟，并与GPS的反演结果进行比较，结果表明两者测量的长波长空间差分结果吻合。搭载在卫星平台的多光谱仪器也可以用于评估大气对干涉相位的影响，文献[12]验证了中分辨率成像光谱仪（medium resolution imaging spectrometer，MERIS）数据集校正大气水汽引起相位延迟的能力。此外，文献[13]利用MERIS多光谱数据与中尺度气象模型MM5联合估计对流层延迟相位，有效地减少了大气效应的影响。近年来，欧洲中期天气预报中心（European Centre for Medium-Range Weather Forecasts，ECMWF）提供的全球中分辨率数值大气再分析产品也被广泛应用，如文献[14]利用ECMWF再分析资料（ECMWF reanalysis，ERA）-Interim（ERA-I）数据集来估计对流层延迟^[14]，在后续研究中又利用多种全球大气模型进行验证^[15]。文献[16]结合GPS观测与高分辨率ECMWF资料开发了一种全球覆盖、近实时的通用InSAR大气改正在线服务（generic atmospheric correction online service for InSAR，GACOS）^[16]，目前得到了广泛应用。（3）根据地形与相位间的相关关系进行改正的经验模型改正法，如线性模型或幂律模型改正。文献[17]通过相位/高程相关分析，对每幅干涉图的对流层延迟相位进行估计和校正，利用ERA40全球大气再分析数据验证了校正结果。文献[18]提出了一种幂律校正方法，该方法考虑了大气性质的空间变化特性，通过幂律关系描述了相位延迟随高度的变化。

由于每种改正方法都存在各自的局限性，不同改正方法对垂直分层和湍流混合两个延迟部分的敏感程度不同（例如，基于气象模型的改正方法通常无法很好地估计湍流效应引起的延迟分量^[19]），难以做到在任何区域和时间均能有效削弱对流层延迟的影响。本文以青藏高原西北缘为研究区域，基于小基线集（small baseline subset，SBAS）原理，利用哨兵1号（Sentinel-1）宽幅影像数据进行时序InSAR处理，考虑到研究区域范围广、高程差异大、气象站点分布稀疏，选择了常用的经验模型线性改正、空间分辨率最高的GACOS改正和时间分辨率最高的ERA5数据集改正3种对流层延迟改正方法，讨论它们在该研究区域形变速率反演中的效果和适用性。

1 研究区域与数据源

青藏高原位于亚洲中部，是世界上最高最大的造山带，其北部存在一条规模巨大的左旋走滑断裂——阿尔金断裂带，该断裂带长达2 000 km^[20]，在印度-欧亚板块碰撞中占据独特的地位，它现今的运动特征对识别青藏高原变形模式具有独特的意义。本文选取了3个位于青藏高原西北缘的研究区域（D136、D165、D19），其中D136位于76.5°E~79.7°E，D165位于80.5°E~83.7°E，D19位于84.9°E~88.1°E。

如图 1所示，北东东向的阿尔金断裂带将青藏高原与塔里木盆地分隔开来，地形起伏强烈，高程差异可达5 000 m。自北向南，随着海拔增高，塔里木暖温带荒漠逐渐过渡为高山荒漠，在海拔5 km以上出现永久冰川。印度-欧亚板块之间的相对运动被青藏高原内部的地壳缩短和走滑剪切调整吸收，北部的柴达木盆地、祁连山脉和阿尔金山脉吸收15%~17%^[21]。尽管阿尔金断裂带构造活动明显，但历史记录到的强震（Mw > 7）的数量很少，2017—2019年研究区域内未发生过强震，研究区域及周边的地震分布情况见图 1。

图  1  研究区域与研究数据覆盖范围

Figure  1.  Research Area and Data Coverage

下载: 全尺寸图片 幻灯片

本文研究数据采用Sentinel-1宽幅（interferometric wide swath，IW）降轨数据，基于GAMMA平台对D136、D165、D19轨道近3 a（2017—2019年）的单视复数（single look complex，SLC）影像数据进行差分干涉处理，采用90 m分辨率的DEM去除差分干涉图中的地形相位，利用MintPy（Miami InSAR time-series software in Python）程序包进行时序分析并去除对流层延迟相位，所使用Sentinel-1数据的时空基线如图 2所示。

图  2  时空基线网络图

Figure  2.  Spatial-Temporal Baseline Networks

下载: 全尺寸图片 幻灯片

气象模型数据为ECMWF最新发布的全球天气再分析数据集ERA5和GACOS。ERA5数据集融合了地面观测、探空气球、卫星等气象观测数据，将地表到80 km高度的大气分为137层，网格空间分辨率为30 km（约0.25°×0.25°），时间分辨率精确到1 h，可以提供1979年至今每小时的大气、陆地和海洋气候变化数据。相比于ECMWF的前一代产品ERA-I数据集，ERA5具有更高的时间分辨率和空间分辨率^[22]，使得气象模型估计更接近于真实的大气情况。GACOS提供了SAR影像获取时间对应的对流层延迟数据，其数据源包含最高分辨率的ECMWF数值天气模型（水平分辨率0.1°×0.1°，时间分辨率6 h）、GNSS对流层延迟产品、空间分辨率90 m的航天飞机雷达地形测绘使命（shuttle radar topography mission，SRTM）与先进星载热发射和反射辐射仪全球数字高程模型（advanced spaceborne thermal emission and reflection radiometer global digital elevation model，ASTER GDEM）地形数据。

2 研究方法

2.1   时序InSAR技术

时序InSAR技术的出现大大提高了InSAR地表观测的精度，目前在时序InSAR技术上发展了很多软件与算法，如永久散射体差分干涉测量（permanent scatterer InSAR，PS-InSAR）^[23]、SBAS^[24]、斯坦福永久散射体干涉（Stanford method for persistent scatterers，StaMPS）^[25]等，不同方法具有各自的优缺点，应根据实际条件（研究区域变形特点、误差分布情况等）选出最优方法进行应用。由于本文研究区域位于青藏高原西北缘，地形差异明显，部分区域存在冰雪覆盖，相干性受到影响，因此采用基于SBAS原理的时序处理算法进行InSAR时序分析。

基于SBAS的时序InSAR技术，每个像素点的相位φ可以表示为：

式中，${\varphi }_{\mathrm{d}\mathrm{i}\mathrm{s}}$、${\varphi }_{\mathrm{u}\mathrm{n}\mathrm{w}}$、${\varphi }_{\mathrm{a}\mathrm{t}\mathrm{m}}$、${\varphi }_{\mathrm{o}\mathrm{r}\mathrm{b}}$、${\varphi }_{\mathrm{d}\mathrm{e}\mathrm{m}}$和${\varphi }_{\mathrm{n}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{e}}$分别为地表位移相位、解缠误差相位、大气延迟相位、轨道误差相位、DEM相位和其他噪声相位。数据处理时采用统计费用流方法（statistical⁃cost network⁃flow algorithm for phase unwrapping，SNAPHU）^[26]进行解缠，在解缠过程中进行了解缠误差的校正；在卫星数据处理时利用精密轨道数据去除轨道误差；由于DEM数据自身带有误差，在数据处理中不可避免地会将这些误差带入差分干涉图中，假定研究区域的变形为线性变化，常规的SBAS改正方法可获取可靠的高精度InSAR时序结果；大气误差主要考虑了对流层延迟的影响，本文采用地形线性相关的经验方法、GACOS改正和ERA5数据集改正3种方法校正对流层延迟相位。

对Sentinel-1数据提供的SAR影像差分干涉成图并进行时间序列分析，评估3种对流层延迟校正方法的性能，技术流程如图 3所示。时序分析基于MintPy实现，该软件是InSAR时序分析的开源程序包，通过读取解缠后的干涉图信息来生成LOS的地表形变信息^[27]。

图  3  时序InSAR对流层延迟改正方法流程

Figure  3.  Flowchart of Time-Series InSAR Tropospheric Delay Correction Method

下载: 全尺寸图片 幻灯片

2.2   大气改正方法

对流层延迟主要受大气中水汽成分的影响，其中垂直分层效应引入了高程相关的延迟分量，湍流效应引入了局部延迟分量。垂直分层效应被认为是特定区域和时间范围内的静态分层，而湍流效应则是动态的过程。如果在一段时间内干涉图主要受静态对流层延迟影响，其延迟信号可以认为是主从影像对应延迟信号的差值，由此可得到简单的线性关系^[28]：

式中，${\phi }_C$为单幅干涉图的整体相位偏移量；K为高程h与对流层延迟$\mathrm{\Delta }{\phi }_{\mathrm{t}\mathrm{r}\mathrm{o}\mathrm{p}}$的传递函数。

在进行对流层延迟校正时，分别采用了线性改正、GACOS改正和ERA5数据集改正3种方法。由于InSAR相位延迟与高程存在一定的线性关系，通过拟合延迟相位与地面高程的线性关系可以得到对流层延迟相位图；GACOS改正根据SAR影像的获取时间，利用对流层分解的迭代插值模型将对流层延迟分解为垂直分层延迟和湍流混合延迟，生成高分辨率的天顶对流层延迟图^[29]，再引入卫星的入射角信息即可获得LOS延迟图；ERA5数据集改正根据雷达影像的获取时间，提取温度、湿度、压力等气象参数，结合研究区域的DEM和入射角等信息，利用双线性插值计算不同高度的天顶大气延迟，然后投影到LOS，从而获取每个雷达影像历元的对流层延迟^[30]。

2.3   评估标准

标准差是衡量大气改正效果的一个重要指标，地表未发生形变时可以用来估计整体的大气改正效果^[31]。由于断层附近存在形变，直接计算标准差会引入误差，因此，本文将断层附近的形变区域掩模后计算非形变区的标准差。一般来说，地表形变主要发生在距离断层3D_G的范围以内，D_G为断层的闭锁深度^[32]，根据现有研究，本文研究区域内断层的闭锁深度约为6~10 km，因此将断层两侧30 km作为掩模区域，掩模后标准差计算公式为：

式中，$\sigma$表示标准差；M为像素点个数；V_i为第i个像素点的形变速率。

然而，标准差虽然能够衡量干涉图整体的噪声水平，但不能反映信号的空间结构，文献[2]研究了随机窗口下不同空间尺度的平均标准差，结果表明信号在空间尺度存在一些特征，需要更多指标来描述InSAR大气延迟的空间特性。本文增加了形变速率与地形的相关系数R作为评价对流层延迟改正效果的指标，其计算公式为：

式中，V_i和h分别表示第i个像素点的形变速率和高程；Cov为协方差；Var为方差。

3 青藏高原西北缘对流层延迟校正性能分析

3.1   D136区域对流层延迟校正

D136区域地形变化剧烈且有积雪覆盖，导致该地区相干性较差，以沿LOS从地面到卫星方向为正，图 4展示了D136原始形变速率图与经验模型线性改正、GACOS改正以及气象模型ERA5数据集改正的结果。

图  4  D136区域对流层延迟校正前后形变速率图

Figure  4.  Deformation Velocity Map Before and After Tropospheric Delay Correction for D136

下载: 全尺寸图片 幻灯片

从图 4中可以看出，对流层延迟改正后，D136东北部原本显示的形变速率大幅削减，但ERA5改正后仍残留了一部分大气误差。线性改正、GACOS改正和ERA5改正后形变速率标准差分别减少了41.05%、38.76%和10.05%。由此可见，3种大气校正方法具有一定的改正效果，其中线性改正和GACOS改正效果更加明显，尤其是线性改正效果最佳。

在时序处理过程中，湍流延迟影响被削弱，垂直分层延迟占主导地位，原始形变速率与地形呈现较大的相关性。在D136区域选取A₁$A_{\mathrm{1}}^{\text{'}}$剖面进行分析，分析结果如图 5所示，灰色部分为地形剖面，黑色、红色、绿色、蓝色虚线分别代表原始形变速率剖面、线性改正、GACOS改正、ERA5改正后的形变速率剖面。由图 5可以看出，剖面A₁$A_{\mathrm{1}}^{\text{'}}$原始形变速率与地形存在明显的相关性，且海拔较低的区域对流层延迟改正效果明显，随着高程增加，改正前后形变速率变化逐渐减小，尤其是海拔高于4 000 m的地区改正效果不明显。

图  5  地形与形变速率剖面（A₁$A_{\mathrm{1}}^{\text{'}}$）图 Profile of Topography and Deformation Velocity A₁$A_{\mathrm{1}}^{\text{'}}$

下载: 全尺寸图片 幻灯片

为了更加准确地评估3种方法的改正效果，本文计算了剖面上形变速率与地形的相关系数R、形变速率标准差σ以及它们的减少量（图 6）。从相关系数来看，GACOS改正效果更好，改正后相关系数降低了76.71%；从标准差来看，线性改正效果最优，改正后标准差降低了50.48%。

图  6  剖面A₁$A_{\mathrm{1}}^{\text{'}}$上对流层延迟改正方法统计结果

Figure  6.  Statistics Results of Different Tropospheric Correction Methods on Profile A₁$A_{\mathrm{1}}^{\text{'}}$

下载: 全尺寸图片 幻灯片

3.2   D165区域对流层延迟校正

D165位于阿尔金断裂带西段，断层活动明显，2008年和2012年在该区域于田地区发生了Mw 7.1和Mw 7.3两次强震。图 7展示了D165原始形变速率图与经验模型线性改正、GACOS模型改正以及气象模型ERA5数据集改正的结果，同样计算了非形变区的标准差。

图  7  D165区域对流层延迟校正前后形变速率图

Figure  7.  Deformation Velocity Map Before and After Tropospheric Delay Correction for D165

下载: 全尺寸图片 幻灯片

从图 7中可以看出，经过对流层延迟校正后，塔里木盆地原本的高形变速率明显减小，阿尔金断裂带附近的形变特征更加明显。非形变区标准差显示，线性改正、GACOS改正和ERA5改正分别使标准差减少了59.21%、55.97%和30.11%，同样是线性改正效果最好。

结合剖面分析，在D165区域选取A₂$A_{\mathrm{2}}^{\text{'}}$剖面，低海拔区域的对流层延迟改正效果十分显著，随高程增加改正效果逐渐减弱（见图 8）。计算剖面上形变速率与地形的相关系数R以及形变速率标准差σ，线性改正后R与σ分别减少了82.67%和64.16%，GACOS改正使两者减少了60%左右，ERA5改正效果不突出（见图 9）。

图  8  地形与形变速率剖面图

Figure  8.  Profile of Topography and Deformation Velocity

下载: 全尺寸图片 幻灯片

图  9  剖面A₂$A_{\mathrm{2}}^{\text{'}}$上对流层延迟改正方法统计结果

Figure  9.  Statistics Results of Different Tropospheric Correction Methods on Profile A₂$A_{\mathrm{2}}^{\text{'}}$

下载: 全尺寸图片 幻灯片

3.3   D19区域对流层延迟校正

D19位于阿尔金断裂带中段，该区域地震活动性并不高。图 10展示了D19原始形变速率图与经验模型线性改正、GACOS模型改正以及气象模型ERA5数据集改正的结果，并计算了非形变区的标准差。

图  10  D19区域对流层延迟校正前后形变速率图

Figure  10.  Deformation Velocity Map Before and After Tropospheric Delay Correction for D19

下载: 全尺寸图片 幻灯片

对流层延迟校正后，西北部塔里木盆地以及断层附近的LOS向形变速率更加平缓，与D136和D165区域不同的是，D19区域线性改正、GACOS改正和ERA5改正后标准差分别减少了25.13%、30.73%和20.15%，GACOS改正在D19区域的适用性更好。

选取剖面A₃$A_{\mathrm{3}}^{\text{'}}$，如图 11所示，原始形变速率与地形相关性呈现一定的相关性，但与D136、D165相比较小，对流层延迟改正后形变速率变化出现了不同程度的减缓，GACOS改正效果更加明显。通过计算，GACOS改正后相关系数与标准差分别减少了59.27%和53.42%，线性改正和ERA5改正效果略差。

图  11  地形与形变速率剖面图

Figure  11.  Profile of Topography and Deformation Velocity

下载: 全尺寸图片 幻灯片

图  12  剖面A₃$A_{\mathrm{3}}^{\text{'}}$上对流层延迟改正方法统计结果

Figure  12.  Statistics Results of Different Tropospheric Correction Methods on Profile A₃$A_{\mathrm{3}}^{\text{'}}$

下载: 全尺寸图片 幻灯片

3.4   GPS数据对比

从3个区域的对流层延迟校正结果来看，线性改正与GACOS改正在青藏高原西北缘的适用性更好。为了更加准确地评估3种改正方法，本文计算了研究区域内GPS站点^[33]投影到LOS的形变速率并与对流层延迟改正前后的结果进行比较。

研究区域内共选取了20个GPS站点，分布如图 13所示，图 13中展示的形变速率为GACOS改正的结果。将每个站点LOS向形变速率与InSAR形变速率进行对比，图 14展示了InSAR形变速率与GPS形变速率的对比结果，图 14（a）~14（c）分别为线性改正、GACOS改正、ERA5改正前后InSAR反演的形变速率与GPS形变速率的对比，图 14（d）表示各个站点InSAR与GPS形变速率差值的绝对值。

图  13  GPS站点分布图

Figure  13.  GPS Stations Distribution

下载: 全尺寸图片 幻灯片

图  14  对流层延迟校正前后InSAR形变速率与GPS形变速率对比

Figure  14.  Comparison of InSAR Deformation Rate and GPS Deformation Rate Before and AfterTropospheric Delay Correction

下载: 全尺寸图片 幻灯片

图14（a）~14（c）显示原始形变速率的均方根误差（root mean square error，RMSE）为5.72 mm/a，线性改正、GACOS改正和ERA5改正后均方根误差分别为为3.08 mm/a、2.79 mm/a和3.69 mm/a，分别降低了46.07%、51.28%和35.51%，GACOS改正效果最佳。从图 14（d）中可以看出，未进行对流层延迟改正的形变速率与GPS观测值相差很大，最多可以造成约10 mm/a的误差，对流层延迟改正后各点残差减小且变化趋于平缓，改正前各点残差的平均值为5.1 mm/a，改正后残差平均值分别为2.26、2.31、3.27 mm/a。综合均方根误差与残差的结果，线性改正与GACOS改正效果更好。

4 讨论与分析

在本文的时间序列内，研究区域未发生过Mw > 5地震等会引起较大形变信号的事件，可以认为远离断层区域的信号是大气延迟、剩余轨道误差、时空失相干等噪声相位的叠加，而断层附近还包含由断层运动引起的信号。塔里木盆地地质构造稳定，该区域在原始形变速率图中显示的形变信号主要由对流层延迟引起，因此，对流层延迟改正在形变速率反演中十分重要。

D136、D165、D19这3个研究区域的非形变区标准差表明线性改正和GACOS改正效果更适用于当前研究区域。从剖面分析来看，A₂$A_{\mathrm{2}}^{\text{'}}$剖面的原始相关系数最大，线性改正效果突出，而A₃$A_{\mathrm{3}}^{\text{'}}$的原始相关系数最小，GACOS改正效果则优于线性改正，ERA5改正在3个研究区域剖面中的适用性均较差。这主要是由于时序InSAR分析中湍流效应的影响被削弱，垂直分层延迟起主导作用，且A₂$A_{\mathrm{2}}^{\text{'}}$剖面形变速率与高程的相关系数高达0.93，这使得线性模型改正具有更大的优势，可以更好地拟合形变速率与地形的相关关系，从而获得相对较好的改正效果；A₃$A_{\mathrm{3}}^{\text{'}}$剖面的相关系数为0.60，GACOS以其高空间分辨率和多源数据集融合的特点发挥了更大的优势；而A₁$A_{\mathrm{1}}^{\text{'}}$剖面的相关系数为0.68，线性改正与GACOS改正效果相当。ERA5数据集在众多气象模型数据中时间分辨率最优，但根据ECMWF在网站发布的全球监测站分布可知，在本文的研究区域内，地面天气报告测站分布稀疏，缺少充足的气象观测数据，并且极端天气条件也会对观测数据产生影响，导致气象再分析资料ERA5数据集的质量较低，从而限制了ERA5数据集改正的效果。

从3个剖面中还可以看出，随海拔高度增加，对流层延迟校正前后形变速率变化逐渐减小。由于对流层延迟主要取决于大气水汽含量，随高度增加水汽含量降低，对流层延迟效应也会随之减小，因此剖面呈现出随海拔高度增加对流层延迟改正效果减弱的情况。

此外，InSAR与GPS数据比较时发现存在一些站点对流层延迟改正后残差变大。对于线性改正来说，式（2）中系数K由全局回归分析所得，本文研究区域范围较广且地形变化剧烈，其改正结果对局部区域不完全适用。对于GACOS改正和ERA5改正来说，两者都会受到时空分辨率的影响，无法准确反映SAR影像获取时间的真实大气条件，从而导致部分站点误差增大。

5 结语

本文利用Sentinel-1 IW影像数据研究了青藏高原西北缘2017—2019年的震间形变速率，评估了经验模型线性改正、GACOS数据集改正和ERA5改正3类方法在对流层延迟改正中的效果。位于76.5°E~79.7°E的D136区域对流层改正后非形变区标准差分别降低了41.05%、38.76%、10.05%，剖面原始形变速率与地形相关系数为0.68，剖面上GACOS与线性改正效果相当；位于80.5°E~83.7°E的D165区域对流层延迟改正后非形变区标准差分别降低了59.21%、55.97%、30.11%，剖面原始形变速率与地形相关系数为0.93，剖面上线性改正效果最优；位于84.9°E~88.1°E的D19区域对流层延迟改正后非形变区标准差分别降低了25.13%、30.73%、20.15%，剖面原始形变速率与地形相关系数为0.60，剖面上GACOS改正效果最优。与GPS数据比较发现，3种方法改正后RMSE分别降低了46.07%、51.28%、35.51%，残差平均值分别降低了55.69%、54.71%和35.88%。

研究结果表明，3种方法均可以削弱青藏高原西北缘的对流层延迟影响，其中线性改正和GACOS改正效果好，ERA5受限于测站的分布密度，改正效果稍差，因此线性改正和GACOS改正方法在青藏高原西北缘的适用性更强。由于地形、海拔和地理位置等因素都会对大气延迟造成影响，因此不同区域所适用的大气改正方法和改正效果都需要结合具体情况进行综合研究分析。