合成孔径雷达干涉测量（InSAR）技术是一种基于卫星遥感的空间大地测量新技术，具有全天时、全天候、高精度、覆盖范围广等特点，在地表形变监测中得到广泛应用^[1]. 当前，针对城市建构筑物，如桥梁^[2]、地铁^[3]、高铁^[4]等的高精度精细化时序InSAR形变监测技术是研究的热点. 然而，InSAR技术不可避免地受到大气相位延迟的影响，极大地降低了InSAR技术提取地表形变的可靠性^[5]. 大气延迟主要表现为三类：第一类是长波长大气延迟，由SAR影像覆盖范围内均匀分布的大气水汽引起，因其表现形式与轨道误差相似，可在校正轨道误差时一并去除^[6]；第二类是与地形相关的大气延迟，由一定范围内近地表高程变化导致的平均水汽含量变化引起^[7]，针对此类延迟，一般通过建立相位与高程的函数模型^[8]等方法予以去除；第三类是湍流大气延迟，由一定范围内大气水汽剧烈变化引起，在干涉图上表现为数公里范围的随机干涉纹理图^[9]，其在时间域上表现为高频信号，在空间域上表现为低频信号，因其在时空域上的随机特性，湍流大气延迟很难使用确定性模型描述^[10]，是处理InSAR大气延迟的难点.

目前，针对湍流大气延迟的改正方法主要有外部数据法^[11-12]和时空滤波法^[13]. 外部数据法可提供大范围、较精确的大气改正结果^[14]，但其改正精度受外部数据源的站点距离、采集时间、时空分辨率、观测条件等多种因素限制^[15]，难以达到预期效果. 此外，基于外部数据构建的大气模型精度严重依赖模型的初始条件与边界条件，难以确保校正结果的稳定性^[16]. 时空滤波法一般能有效缓解平坦区域的湍流延迟^[17]，该方法基于时间域的相位随机性假设，通过滤波的方法对大气延迟相位进行平滑^[18]，然而，当湍流延迟影响较大时，仅使用时空滤波法无法有效将其滤除，而过度的平滑可能导致形变信息的错误估计^[19].

本文考虑湍流大气延迟在时空域上的随机特性，提出一种引入Baarda粗差探测的InSAR大气校正方法，采用Baarda粗差探测理论，识别并去除对时序形变影响大的观测相位. 论文采用模拟数据和实测SAR影像开展时序InSAR形变监测实验研究，对本文方法与仅时空滤波法进行对比分析，最后，将本文方法应用于巴基斯坦拉合尔橙线轨道交通沿线的变形监测.

1.   引入Baarda粗差探测的InSAR大气校正方法

Baarda粗差探测^[20]是一种将粗差纳入函数模型的有效粗差探测与识别方法. 考虑InSAR时序形变相位中湍流大气延迟的随机特性和对相位的剧烈影响，本文提出引入Baarda粗差探测法的大气校正方法，其流程如图1所示. 该方法主要分两步：

图  1  引入Baarda粗差探测的InSAR大气校正方法

Figure  1.  InSAR tropospheric correction method incorporating Baarda data snooping

下载: 全尺寸图片 幻灯片

步骤1　基于Baarda粗差探测理论的湍流大气延迟相位识别与剔除.

1）在时间域上，归一化像元$i$的时序原始相位$ {\phi _i} $和成像时间$ t $，获取对应的归一化时序相位$ {\phi _{i,{\mathrm{n}}}} $和时间$t_{\mathrm{n}}$. 初始化/更新像元$ i $在时间序列上的粗差点集$ {{\boldsymbol{P}}_i} $. 基于先验变形模型$ {\phi _{i,{\mathrm{n}}}} = {\boldsymbol{G}}(t_n) $，${\boldsymbol{G}}(t_{\mathrm{n}}) $为变形模型矩阵，构建像元$ i $的间接平差模型，获取模型残差$ {{\boldsymbol{v}}_i} $，如式（1）所示.

式中：y为观测值向量，G为与$ t' $相关的变形模型矩阵，x为待估参数向量，e为观测误差向量；$ \widehat {{{\boldsymbol{v}}_i}} $、$ \widehat {\boldsymbol{x}} $分别为$ {{\boldsymbol{v}}_i} $和x的估计值向量.

2）构建像元$i$的时序正态分布向量$ {{\boldsymbol{u}}_i} $为

式中：$ \sigma {_i} $为像元$ i $的标准差，$ {{\boldsymbol{Q}}_{{{\boldsymbol{v}}_i}{{\boldsymbol{v}}_i}}} $为协因数矩阵，$I$为单位矩阵，$ n $为像元$ i $在时间序列上不含粗差的影像数，$ s $为自由度.

3）对于$ {{\boldsymbol{u}}_i} $，假设不存在粗差，选定显著水平α进行假设检验，$ {k_0} = {\mu _{1 - \alpha /2}} $为该显著水平下的正态分布值. 若$ \forall {{\boldsymbol{u}}_i} $≤$ {k_0} $，假设为真，$ {\phi _{i,{\mathrm{n}}}} $中不存在粗差，否则假设为假，$ {\phi _{i,{\mathrm{n}}}} $中存在粗差.

4）若存在粗差，将$ \max (\left| {{{\boldsymbol{u}}_i}} \right|) $对应的点记为粗差，将其剔除后重复上述步骤. 若不存在粗差，保存像元$ i $在时间序列上的粗差点集$ {{\boldsymbol{P}}_i} $.

5）重复上述步骤，直到获取全部像元在时间序列上的粗差点集$ {\boldsymbol{P}} $.

步骤2　基于时空滤波的大气延迟校正.

1）通过时间域一维高斯滤波，将相位$ \phi $分解为时间域低频相位$ {\phi _{{\text{lpt}}}} $和时间域高频相位$ {\phi _{{\text{hpt}}}} $.

式中：$ \Delta {t_j} $为$ t $与时刻$j$的差值，$ {M_{ij}} $为像元$ i $在$j$时刻的序数，$ {w_i}(\Delta {t_j}) $为像元$i$在$j$时刻的一维高斯核函数，$ \beta $为时间域滤波阈值，$k$为时间序列影像数.

2）通过空间域二维高斯滤波，将$ {\phi _{{\text{hpt}}}} $分解为空间域低频相位$ {\phi _{l{\text{p}}s}} $和空间域高频相位$ {\phi _{hps}} $.

3）将$ \phi $分解为形变相位$ {\phi _{def}} $和大气延迟相位$ {\phi _{atm}} $，完成局部湍流大气的识别与校正.

2.   基于模拟数据的大气延迟实验

2.1   SAR干涉图模拟

数据模拟理论和模型依据基于文献[8，21]. 模拟的SAR干涉图覆盖范围为10 km × 10 km，空间分辨率约80 m，以12 d为重访周期，共生成36景SAR时序干涉相位图. 以90 m空间分辨率的宁夏、银川周边航天飞机雷达地形测绘任务（SRTM）提供的数字高程模型（DEM）为参考地形，该区域海拔为1821～1965 m，模拟与地形相关的大气延迟相位$ {\phi _{ver}} $，范围为−π～π rad；以Peaks函数模拟线性的形变相位$ {\phi _{def}} $，模拟观测周期内最大形变速率约100 mm/年；以Kolmogorov模型模拟湍流大气相位$ {\phi _{tub}} $，范围为−4π～4π rad；模拟随机噪声相位$ {\phi _{noi}} $，范围为−0.5～0.5 rad. 图2显示了第17期模拟干涉图各分量及总相位$ {\phi _{{\text{tol}}}} $，$\phi _{{\text{tol}}} $如式（6）所示，图中，$ (x,y) $构成空间位置坐标.

图  2  第17期模拟SAR干涉相位组成

Figure  2.  Components of 17th simulated SAR interferometric phase

下载: 全尺寸图片 幻灯片

2.2   大气校正与评价

基于引入Baarda粗差探测的InSAR大气校正方法，以线性形变模型估计模拟时序变形序列的残差，构建统计量，采用99%的置信度进行Baarda粗差探测，对剔除粗差后的时序像元进行时空滤波，依照图1流程滤除大气延迟相位，得到校正后的变形时间序列. 为验证本文方法的大气校正效果，设置相同的时空滤波参数，采用无粗差探测机制的时空滤波法进行大气校正，得到变形时间序列.

采用线性模型对获取的变形时间序列进行拟合，得到原始形变时序、本文方法和仅用时空滤波法校正后的时序形变速率如图3（a）～（c） 所示，图3 （d）～（f）为对应的速率残差. 以形变速率±1 mm/年为阈值，将研究区域划分为稳定区域和形变区域.

图  3  模拟变形速率及残差

Figure  3.  Simulated deformation rates and residuals

下载: 全尺寸图片 幻灯片

图4展示了原始相位、本文方法和仅时空滤波法校正后在稳定区域和形变区域的速率残差统计直方图及其正态分布期望μ与标准差σ. 可以看出，在稳定区域和形变区域，本文方法的速率残差期望均处于0 mm/年附近，标准差相较于仅时空滤波法和原始相位分别降低了约25.8%、28.6%和16.0%、16.4%. 本文方法在形变区域的校正精度相较于稳定区域略有降低，但结果仍显著优于仅时空滤波法与原始相位结果.

图  4  模拟数据速率残差直方图

Figure  4.  Histogram of simulated data rate residuals

下载: 全尺寸图片 幻灯片

为验证本文方法对时间序列形变提取的效果，在稳定区域随机选取3000个像元作为原始序列并取均值，绘制本文方法与仅时空滤波法获取的时序变形相位均值，如图5所示. 可以看出，本文方法有效剔除了时间序列中受大气延迟影响严重的极端值，累计变形量处于0 rad附近，而采用仅时空滤波法的累计变形明显受到极端值的影响；相比仅时空滤波法，本文方法将相位标准差由0.86 rad下降至0.33 rad，降低了约62.1%.

图  5  稳定区域时间序列形变结果对比

Figure  5.  Comparison of time-series deformation in stable region

下载: 全尺寸图片 幻灯片

3.   Sentinel-1 SAR数据的大气延迟实验

3.1   研究区域与数据处理

拉合尔位于巴基斯坦旁遮普省，地处印度河流域冲积平原，地表沉积物覆盖厚度达300 m，地势平坦，平均海拔约200 m. 拉合尔属于亚热带大陆气候，气温和降水量随季节显著变化^[22-23]. 拉合尔地理位置如图6所示，黄色矩形为Sentinel-1 A影像覆盖范围，红色矩形为研究区域.

图  6  研究区及Sentinel-1A影像覆盖范围

Figure  6.  Study area and coverage of Sentinel-1A SAR imagery

下载: 全尺寸图片 幻灯片

本文使用2021年4月2日至2022年9月12日共计44景Sentinel-1A SAR 升轨影像数据，开展拉合尔时序InSAR变形监测研究，监测周期528 d，数据入射角约为39.3°. 时序干涉变形监测采用SBAS-InSAR方法，选取2021年4月2日SAR影像作为主影像，使用精密轨道数据辅助影像配准，为尽量避免相位解缠错误，使监测结果更有效反映地面沉降和空间分布变化趋势^[24-25]，将所有影像两两进行全分辨率干涉处理，共获取干涉对946对，进而采用30 m空间分辨率SRTM DEM数据进行DEM相位去除，通过设置振幅离差指数小于0.25，获取35385个相干像元，得到研究区域原始时序变形序列.

3.2   大气校正与评价

基于引入Baarda粗差探测的InSAR大气校正方法，以线性形变模型估计时序变形残差，构建统计量，采用97%的置信度进行Baarda粗差探测，对剔除粗差后的时序像元进行时空间滤波，依照图1流程滤除大气延迟相位，得到校正后的变形时间序列. 为验证本文方法的大气校正效果，设置相同的时空滤波参数，以仅时空滤波法校正结果为对照，图7展示了约3个月间隔随机选取的时序处理结果.

图  7  仅时空滤波的时序形变相位及统计直方图

Figure  7.  Time-series deformation phase and statistical histogram using only spatiotemporal filtering

下载: 全尺寸图片 幻灯片

用相同间隔绘制本文方法时序处理流程与统计结果，如图8所示. 实测Sentinel-1 SAR数据的大气校正. 可以看出，本文方法在真实场景下有效、准确识别出时序影像中受大气延迟影响大的像元并剔除，且探测粗差点的迭代次数与该像元受大气延迟影响程度密切相关，迭代次数越少，大气延迟污染越严重，2021年4月14日影像没有探测到粗差，该影像几乎未受到大气延迟污染. 对比图8（a）、（d）可以看出，原始时序相位统计直方图无明显规律，完成本文方法的大气校正后，随着形变量的增大，统计直方图出现两个明显的峰值，其中一个峰值位于0 rad附近，代表稳定区域，另一个峰值随时间累积逐渐增大，代表变形区域. 结合图7、8（d）可以看出，仅时空滤波法能够有效削弱湍流大气影响，但由于部分像元受湍流大气延迟过于严重，其结果在拉合尔城市周边沉降稳定区域出现了错误的形变信息. 本文方法将受湍流延迟影响较大的像元判定为时间序列上的粗差并剔除，从而获得了更准确的形变结果.

图  8  实测Sentinel-1 SAR数据的大气校正

Figure  8.  Atmospheric correction of measured Sentinel-1 SAR data

下载: 全尺寸图片 幻灯片

基于稳定区域像元相位是评价时序InSAR大气延迟相位校正效果的主要方法之一^[26]. 本文以雷达视线向线性形变速率±10 mm/年为稳定像元，构建评价数据集，采用相位标准差和半变异函数2种指标^[27]对本文方法进行评价.

相位标准差表示为

式中：$ m $为评价数据集像元数，$ {\phi ^*} $、$ \overline {{\phi ^*}} $分别为评价数据集像元相位及其均值.

半变异函数$ S(\cdot) $表示为两个空间点之间差异的平方平均值，如式（6）所示.

式中：d为滞后距离，设置为500 m； $ {\phi ^*}(p) $、$ {\phi ^*}(p + d) $为位置p处的评价数据集像元相位值.

为提高计算效率，在计算半变异函数时对评价数据集内稳定像元进行了1∶5的降采样.

图9展示了原始变形时序、仅时空滤波法和本文方法的相位标准差和半变异函数. 从图9（a）可以看出，原始变形时序具有较大的相位标准差，而仅时空滤波法和本文方法获取的相位标准差明显减小. 本文方法相较于原始时序相位的相位标准差平均降低约41.2%，且在整个时序中保持一定的稳定性. 本文方法获取的时序稳定点标准差略小于仅时空滤波的结果. 从图9 （b）可以看出，本文方法与仅时空滤波法结果明显优于原始形变时间序列，在5 km以内的空间尺度上，两种变形时序结果精度相当. 随着空间尺度的增大，仅时空滤波法与本文方法校正结果的半变异函数趋于稳定，本文方法校正结果的半变异函数相较于同空间尺度下原始相位半变异函数降低约74%，优于仅时空滤波法的65%.

图  9  统计指标评价结果

Figure  9.  Evaluation results of statistical indicators

下载: 全尺寸图片 幻灯片

4.   拉合尔市精细化变形监测应用

采用本文方法获取的拉合尔地区视线向年均形变速率结果如图10所示. 从图中可以看出，拉合尔地区存在明显的地面沉降，东西向沉降范围约10 km，南北向沉降范围约14 km，沉降区域与非沉降区边界明显，沉降区域与城市建筑区具有一致性，最大年均沉降速率位于城区南部，视线向年均最大沉降速率约为−72.8 mm/年.

图  10  拉合尔市年视线向变形速率

Figure  10.  Annual line of sight （LOS） deformation rate in Lahore

下载: 全尺寸图片 幻灯片

拉合尔市橙线轨道交通总长约26 km，全线共设车站26座^[28]. 通过对图10结果进行克里金插值，并将视线向变形转换为垂直向，设置200 m的缓冲区，本文获得橙线沿线沉降速率，如图11所示. 可以看出，橙线最大沉降发生在Pakistan Mint站附近，沉降速率约为−38.7 mm/年. 据中华人民共和国国家标准《地质灾害危险性评估规范》（GB/T 40112-2021），橙线地面沉降发育程度中处于强发育区、中等发育区和弱发育区分别占线路总长的17.6%、23.1%和59.3%.

图  11  橙线沿线沉降速率及沉降梯度

Figure  11.  Settlement rate and gradient along orange line

下载: 全尺寸图片 幻灯片

线路不均匀沉降直接影响乘客舒适性并对列车运营构成威胁，图11展示了橙线沉降梯度. 可以看出，橙线最大差异沉降出现在GPO站附近，最大沉降梯度约为0.03 mm/（m·年）. 据林辉等^[29]有关桥梁损坏严重等级判定，橙线全线处于基本未受影响等级，不均匀沉降对线路安全影响有限，但随着时间增加，线路不均匀沉降风险也会增大，需持续监测.

5.   结　论

针对时序InSAR大气校正中，传统时空滤波法无法准确估计湍流延迟影响较大像元的形变相位问题，提出一种引入Baarda粗差探测的InSAR大气校正方法，将湍流延迟影响较大的像元判定为时间序列上的粗差并剔除. 本文主要结论如下：

1） Baarda粗差探测法能准确识别时序影像中受湍流大气影响较大的像元. 粗差点与湍流大气延迟相位空间分布一致，且探测粗差点的迭代次数与该像元受大气噪声影响程度密切相关，迭代次数越少，大气相位污染越严重.

2） 大气延迟校正结果表明，本文方法和仅时空滤波法校正结果明显优于原始形变时序，且本文方法优于仅时空滤波法. 模拟数据结果显示，相较于仅时空滤波法，变形速率在稳定区域和形变区域的标准差分别下降了约25.8%和16.0%，稳定区相位标准差下降了约62.1%. 拉合尔Sentinel-1 SAR数据结果显示，本文方法获取的半变异函数相较于同空间尺度下原始相位半变异函数降低约74%，优于仅时空滤波法的65%.

3） 本文方法在巴基斯坦拉合尔地区得到成功应用，监测结果表明，拉合尔地区存在大面积地面沉降，沉降区域与非沉降区边界明显，视线向年均沉降速率达−72.8 mm/年，该市橙线沿线约17.6%的区域处于地面沉降发育程度强发育区间.

致谢：欧洲航天局（ESA）提供Sentinel-1A SAR数据支持；河海大学优硕培育项目（422003519）为本研究提供部分资助.