地基大视场SAR形变监测的非均匀大气相位校正方法

白泽朝 王彦平 王振海 胡 俊 李 洋 林 赟

①(北方工业大学电气与控制工程学院 北京 100144)

②(北方工业大学信息学院雷达监测技术实验室 北京 100144)

③(中南大学地球科学与信息物理学院 长沙 410083)

1 引言

滑坡是我国主要地质灾害之一，其中特大型滑坡，其滑坡体量巨大，往往会造成严重的财产损失和人员伤亡，滑坡地表变形是实现地质灾害监测和预警的重要辅助信息。地基干涉合成孔径雷达(Ground-Based Interferometric Synthetic Aperture Radar,GB-InSAR)是一种高精度的地表形变监测设备，雷达天线在线性轨道上移动，形成方位合成孔径，有效提升了影像方位向分辨率，通过对同一场景的重复扫描和成像，采用干涉测量方法提取目标区域的形变信息[1,2]。与星载干涉合成孔径雷达(Interferometric Synthetic Aperture Radar,InSAR)[3,4]相比，GB-InSAR具有观测角度灵活、测量精度高(可达亚毫米)的优点，是局部区域高精度形变监测的重要监测手段。对于特大型滑坡监测，不仅要求GB-InSAR平台实现远距离覆盖，而且在方位向上具有大视场角度成像能力，来满足特大滑坡广域范围监测需求。

此外，观测范围在方位向上的扩大，使得场景内环境更为复杂，特别是观测中大气变化的影响，大气会引起电磁波折射，从而造成电磁波传播路径和方向的变化。在两次成像时刻大气条件(温度、气压和湿度)的变化会导致不同的传播延迟，并且在相位中将存在大气相位分量(通常称为“大气相位屏”)[5-7]。目前的GB-InSAR大气相位校正方法主要分为3类，一类基于气象数据的方法[5]，依据电磁波传播的折射率与场景内温度、气压和湿度数据之间的关系来校正大气相位。一类基于函数模型方法[8-12]，通常大气影响在空间上具有强相关性，当满足大气在空间均匀性假设时，稳定的永久散射体点相位差在距离方向上的投影呈现线性函数关系，通过求取函数关系即可校正大气相位。此外，考虑了陡峭地形条件下大气折射率会随高程发生改变，有学者进一步考虑了高程对大气延迟系数的影响，给出了距离-高程的多元回归模型[13,14]或者兼顾水平位置和高程等信息的多参数函数模型[15,16]。另一类空间插值的方法[17-21]，当大气在空间非均匀变化时，依据稳定的特征点，通过插值的方式进行大气相位校正。

基于稳定永久散射体(Persistent Scatterers,PS)线性函数模型和距离-高程的多元回归模型的方法易于实现，并且已广泛应用在GB-InSAR形变监测领域。然而，在特大滑坡变形监测场景下，需要在方位向上具有较大的视场角度(>100°)，GB-In-SAR系统方位向范围的扩大，使得观测场景内环境更为复杂，可能存在明显的湍流或者局部水汽变化[13]，导致大气介质在方位方向呈现非均匀的特性[21]。因此，传统的基于稳定PS点函数模型的方法无法在复杂大气条件下准确校正大气相位。

在方位向大视场角度山区复杂大气条件下，大气参数的随机变化，可能在空间上存在线性和非线性两个特定特征的异常相位波动。针对上述问题，本文提出一种两阶段半经验模型校正线性和非线性大气相位。选择具有方位向较大视场角度和跨越长江观测的三峡库区新铺与藕塘滑坡实验区，详细介绍了大气校正的处理流程，并采用本文提出的两阶段半经验模型和常规方法进行大气相位校正实验，验证了算法的有效性。

2 两阶段半经验模型大气相位校正方法

2.1 PS点选择

PS点选择是进行大气相位校正的必要前提，需要根据实际观测场景，避免植被波动和水域等造成误选点情况，选择合适的选点策略。一般根据点的幅度离散和相位标准差的统计特性，在高信噪比的情况下，目标点的稳定性直接由幅度离差指数DA表示为[22]

其中，δA是幅度标准偏差，mA是幅度平均值。幅度离差指数越小，幅度信息越稳定，然而场景中存在水域时，往往会出现误选等情况，可以通过估算相干系数来估计每个像素点的干涉相位质量，滤除水域的误选点，公式为[23]

其中，S1和S2是与干涉图相同像素所对应的复值，|·|为绝对值运算符，E{}为期望值。将符合幅度离差阈值和相干性阈值的点选作PS点。

2.2 距离-高程模型初补偿

GB-InSAR目标点的干涉相位可以表示为[1]

其中，k为整数表示模糊度，φdef表示形变相位分量，φatm表示大气相位分量，φnoise表示误差相位分量。为了获取高精度的形变相位φdef，需要有效的大气相位φatm校正方法。

通过沿路径L积分得到雷达传感器与目标点之间大气相位为[10]

其中，n表示折射率的变化，它随时间t和距离γ而变化，L表示信号的传输路径，λ是波长。一般情况下，假设大气在空间上是匀质的，n在距离γ上不发生变化，在时间t上随机变化，因此可以将干涉图中大气相位分量表示为随距离变化的公式[8]

其中，β0为常数项，β1表示与距离相关的线性系数，γ表示传感器与目标点之间的距离。另外，考虑到地形对大气相位的影响，兼顾距离和地形的大气相位可以表示为[16]

其中，β0为常数项，β1表示与距离相关的线性系数，β2表示与距离、地形相关的系数，γ表示传感器与目标点之间的距离，h表示目标点的高程值。可以建立大气相位与距离和地形的方程组

其中，φ为n个PS点干涉相位构成的n×1维向量，X为常数1和n个PS点距离与高程值构成的n×3维向量，β为待估计的3×1维 向量，ε为相位误差构成的n×1维 向量。采用最小二乘对β进行求解

则大气相位的估计值为

将干涉相位φ减掉估计的大气相位φatm即为大气校正后相位。在线性大气相位校正阶段，基于相干性和幅度离差进行PS点选择，设置较高的阈值选作高质量PS点，采用距离-高程函数模型估算大气相位。此外，高质量PS点中可能存在一些不可靠点，比如形变区的点和相位跳变的点，为了确保估算结果的准确性，将大于二倍相位标准差σ的点剔除，迭代求解距离-高程函数模型参数，重新校正大气相位。

2.3 空间维插值补偿

在校正随高程变化的大气相位后，计算观测场景内高质量PS点形变，设定形变阈值掩模选择出稳定PS点。稳定PS点的相位主要包含大气相位和噪声相位，大气相位在空间上表现出非匀质特征，但在较小的距离范围内可以视为匀质特征，噪声相位在空间上表现为随机特征，不具备空间相关性。在一定的空间范围内进行均值滤波，则可以有效滤除噪声相位影响，同时基于稳定的PS点，采用反距离权重插值算法，插值出采用较低阈值选取的PS点大气相位，从而有效校正非线性大气相位，本文所研究的高质量和低阈值设置在4.2节详细介绍。其中反距离权重插值算法是利用已知点与待插值点的距离来定义权重，根据周边的已知点加权计算待插值点的相位值，距离越近权重越大，其公式为

其中，Z(x,y)为插值点结果，(x,y)为待插值点的空间坐标，Zi为第i个已知点相位值，|di|为第i个已知点与待插值点的空间距离，µ为权重的幂指数，一般取2[20]。此外，考虑到大气相位的空间相关性和计算效率，只选取待插值点周边最近的3个已知点参加计算。

本文算法的大气相位校正流程如图1所示，首先，利用幅度离差阈值和相干性阈值选取PS点，设置较高的阈值选作高质量PS点，参与后续的线性大气相位估计。然后，根据经验判断干涉图中大气误差影响范围和特点，基于高质量PS点采用距离-高程的函数模型对干涉相位校正高程相关的大气影响，在计算中剔除大于二倍相位标准差PS点，迭代求解距离-高程函数模型参数。最后计算序列图像形变，设定阈值选择出稳定PS点，采用反距离权重插值算法估计低阈值选取的PS点大气相位，从而有效校正非线性大气相位。

图1 本文算法大气相位校正流程Fig.1 The workflow for atmospheric phase screen processing steps

理论上，在地形复杂区域距离-高程的函数模型可以提高大气校正的准确性。然而，在实践中干涉图可能同时存在不同特征的大气相位，简单的函数模型可能校正效果无法满足监测需求。为了在复杂大气条件下保持形变测量的可靠性，拓展GB-In-SAR在大视场观测的应用场景，本文方法同时兼顾了线性和非线性大气相位，实现不同特征大气相位的有效校正。

3 研究区描述和系统参数

3.1 研究区描述

研究区为新铺和藕塘滑坡实验区，位于长江三峡库区重庆市奉节县安坪镇。图2为研究区范围，其中雷达布设在长江北岸，跨江观测长江南岸的橘色半椭圆区，观测区长度约为3500 m，宽度约为2500 m，雷达在方位向上的观测视角为120°。新铺滑坡位于长江南岸(右岸)斜坡地带，北抵长江河床，南至泰山庙山脊，总体斜长约2 km，相对高差约300 m，斜坡坡度15°～20°。藕塘滑坡位于新铺滑坡的西侧，属浅中切割单斜低山河谷地貌，岩层倾向与坡向近于一致，长江从藕塘北边由西流向北东，与岩层走向夹角为10°～15°。

图2 研究区范围Fig.2 The study area

3.2 地基大视场SAR系统

实验采用自主研制的地基大视场SAR系统对新铺和藕塘滑坡进行了观测，连续获取了29幅雷达图像，时间从2021年7月27日17:44～22:56。该GBInSAR系统采用角度关联合成大视场成像技术，实现方位向120°以上大视场角度范围的形变监测。如图3所示，其工作在Ku波段，波长为1.74 cm，其方位分辨率在1 km处为5.4 m，距离分辨率为0.37 m。

图3 GB-InSAR系统Fig.3 The GB-InSAR system

4 实验分析

4.1 大气相位分析

本文提出的大气相位校正模型应用于三峡库区新铺和藕塘滑坡形变观测中，并对大气校正效果进行了评估。首先，对相邻图像进行干涉处理，每个干涉图的时间基线约为10 min。根据观测区情况，选择同时满足幅度离差指数和相干系数的点作为PS点，其中设置幅度离差指数0.15，相干系数0.9，筛选出25837个高质量PS点。图4(a)和图4(b)分别为平均幅度图和平均相干系数图，图4(c)为高质量PS点的幅度离差图，可以看出坡体上大部分像素点的幅度稳定性和相干系数都很高，不存在水域和植被区的误选点。

图5显示了受大气影响小的干涉图，图6和图7分别显示了受到随高程变化的大气和非匀质大气影响的干涉图，使用不同的色标表示-π～π范围内观察到的相位变化。干涉对A的影像采集时间分别为19:10与19:20，干涉对B的影像采集时间分别为20:04与20:14，干涉对C的影像采集时间分别为19:53与20:04。在数据采集的间隔时间内，图5(a)中的干涉相位没有超过色标范围，也没有显示出重复的色标变化，图5(b)的高质量PS相位散点也都在0 rad附近波动，表明在2500 m的最大范围内没有观察到相位缠绕现象。另外，由于GB-InSAR的成像几何，可能在雷达近处存在有强目标的干扰，会在部分干涉图零方位角上出现一条相位变化十分均匀的条纹，后续通过PS点选择可以避免对大气相位校正的影响。

图6 干涉对B相位Fig.6 Interferogram B

大气引起的相位变化如图6(a)和图7(a)所示，可以将大气根据空间分布分为两种不同情况。图6(a)为受随高程变化大气影响的干涉相位，从其相位散点图6(b)中可以看出大气在1000～1500 m间具有恒定的相位递减，由于大气条件变化存在明显大气相位分量，该距离范围内干涉图中高质量PS点相位值位于正相位区域的0～1 rad。相位变化沿距离和地形具有恒定的递减，可以通过应用距离-高程模型对大气进行校正。

如图7(a)所示，干涉对C存在更为复杂的非匀质大气影响。由于地基大视场SAR观测场景的扩大，在受到风或者湍流大气影响时，场景内匀质大气的假设不再成立，从图7(b)相位散点可以看出，大气在1000～1500 m处呈现0～1 rad的相位递减，而在1500～2500 m处存在明显的非线性相位波动，因此采用传统的距离-高程的大气相位校正方法建模将难以有效去除该区域大气影响。

传统的基于函数模型的大气校正方法仅适用于小范围的匀质大气，即当气压、温度和湿度沿距离或者地形均匀变化时。因此，传统的大气校正方法不能作为所有干涉图大气校正的通用方法，必须设计一套能够处理复杂大气条件下的算法，使得可以适应地基大视场SAR观测场景，尤其是跨江等水汽复杂的观测中。

4.2 大气相位校正效果分析

由于观测范围广且地形复杂，大气条件变化明显，相应的大气相位误差也较大(如图6和图7)。在第1阶段采用距离-高程模型校正大气相位，图8和图9为在1000～1500 m范围校正随高程变化大气后的干涉图和对应高质量PS点相位散点图，PS点相位在0 rad波动。与图6和图7对比，干涉对B和干涉对C都有效校正了距离-高程相关的线性大气相位，最大的校正相位分量约为1 rad。但是，由于环境复杂，大气在空间上呈现不均匀变化情况，即使校正线性大气相位后，两个干涉图仍然存在非线性大气相位变化。如图9(b)所示，随着距离的变化，高质量PS点的相位呈现非线性变化趋势。

图7 干涉对C相位Fig.7 Interferogram C

图8 距离-高程校正后干涉对B相位Fig.8 Interferogram B after distance-elevation correction

图9 距离-高程校正后干涉对C相位Fig.9 Interferogram C after distance-elevation correction

在第2阶段采用非线性大气相位校正方法，对选取的高质量PS点计算累计形变量，设定形变阈值为5 mm，选择出稳定PS点，然后选取同时满足幅度离差指数0.25和相干系数0.8的75104个较低阈值的PS点，基于稳定的PS点通过插值的方式估算所有PS点的大气相位。如图10和图11为校正后的干涉对B和干涉对C，在校正前干涉对B在沿江边距离向约1000 m位置存在约1 rad相位波动，在距离向约1500 m位置存在约-0.7 rad相位波动(图10(a))。经过非线性大气校正后，补偿了明显的大气相位，干涉对B整体相位在0 rad附近波动。另外干涉对C存在更为复杂的大气相位，在校正前，不仅沿江边约1000 m,1500 m位置存在明显大气相位，在距离向约2500 m最远距离位置也存在约-0.5 rad相位波动(图11(a))，非线性大气校正后，干涉对C得到了有效相位补偿。

图10 非匀质大气校正后干涉对B相位Fig.10 Interferogram B after correcting nonlinear atmospheric

图11 非匀质大气校正后干涉对C相位Fig.11 Interferogram C after correcting nonlinear atmospheric

与小范围观测不同，GB-InSAR在远距离、大视场的观测中显示出更为复杂的大气相位变化，但可以通过选取稳定的有效测量点，提高形变信息的准确性。图12为获取的观测区累计形变量，如图12(a)所示，未校正大气相位时在江沿岸约1000 m和远距离2000～2500 m位置处存在明显大气相位引起的约5 mm形变误差，如图12(b)所示，仅使用距离-高程函数方法在沿江和远距离处大气相位有所削弱，但沿江边区域仍然存在大气相位，造成了约-3 mm的形变误差，如图12(c)所示，使用本文方法有效补偿了沿江和远距离区域大气相位，表明本文算法的大气相位校正结果相位误差显著降低。理论上，在复杂的场景根据大气相位的特征本文方法可以提高大气估计的准确性和形变提取的鲁棒性。从累计形变量图分析，监测区较稳定，在距离向约2000 m位置处，如图12(c)中红色矩形框所示，存在约4.2 mm的形变区，结合已有研究[24]，推测在监测时间段内，滑坡的后缘区存在明显的形变，即在藕塘滑坡存在向下滑动的情况。此外，在形变区外存在零星分布的少量形变点，推测可能是植被扰动造成的。

图12 累计形变量Fig.12 Cumulative deformation

为验证大气相位校正的准确性，挑选了在不同距离向和不同滑坡位置具有代表性的4个像元点变形曲线进行分析，图12(c)中标记了像元点位置，图13为像元点变形曲线图。P1位于新铺滑坡的高程较大位置，经过距离-高程函数校正后反而误差更大，最大误差约为2.4 mm。P2位于藕塘滑坡的高程较大位置，在未进行大气相位校正前，形变的波动幅度较大，最大误差约为-1.6 mm，经过距离-高程函数校正后，形变波动有所变小，但仍然存在约为-1.1 mm误差，本文方法校正后，变形序列变化更为随机且在0.5 mm范围内波动。P3位于监测到的藕塘形变区，总体呈现向下滑动的趋势，由于距离-高程函数方法主要校正沿江和远距离区域随高程变化的大气相位，该区域未进行大气校正，而本文方法大气校正后没有将真实的形变补偿掉，仍然保留了原有的形变趋势。P4位于新铺滑坡靠近长江的沿岸，经过距离-高程函数校正后最大误差由1.5 mm减小为1.1 mm，本文方法校正后，变形序列变化更为随机且在0.5 mm范围内波动。因此，本文方法不仅保留了真实的形变区，而且在非形变区实现了误差的有效降低。

图13 典型点变形曲线对比Fig.13 Cumulative deformation map

通过三峡库区新铺和藕塘特大滑坡的地基大视场雷达图像实验，验证和分析了本文方法的有效性与准确性。在校正江沿岸和远距离位置随高程变化的大气相位，本文方法可以达到与距离-高程函数方法相同甚至更优的精度。在校正非均匀大气相位，函数模型方法往往会失效，采用大量的稳定PS点作为参考点，空间维插值的方式，可以有效保证非均匀大气相位的校正精度。

5 结语

本文采用地基大视场SAR系统识别新铺和藕塘滑坡的形变区，观测到了由大气条件变化引起的干涉相位变化。由于该地区地形和大气条件复杂，传统的大气相位校正方法通常会提供不准确的形变结果。复杂大气环境下的干涉相位不仅存在随距离-高程变化的线性大气相位成分，还存在非线性大气相位成分。因此，提出了一种两阶段半经验模型估计各层大气相位，以消除空间域中的大气相位变化。该方法基于空间模型和统计补偿方法，在不使用温度、气压或湿度等大气参数的情况下，对空间非均匀性大气相位进行补偿，对于场景环境信息较少的GB-InSAR观测，可以有效地消除大气效应和准确地获取形变信息。本文方法不需要在站点和实验阶段收集与处理气象数据，拓展了GB-InSAR的适用性且能更好地对大范围下的大气误差进行校正。

本文方法还存在着一些不足，由于地基SAR系统工作在Ku波段，容易受到植被扰动影响，在形变区外仍然存在零星分布的由于植被扰动造成的误差点，需要在后续工作中进一步研究如何避免植被扰动带来的误差。