面向复杂大气扰动的GB-InSAR相位误差补偿方法

金重阳 刘 毓 邓云开 田卫明，3 胡政权

（1.重庆三峡学院电子与信息工程学院，重庆 404130；2.北京理工大学信息与电子学院，北京 100081；3.北京理工大学重庆创新中心，重庆 401120）

1 引言

GB-InSAR 差分干涉处理分析形变时，需要以高质量像素点作为研究对象，即PS 点。PS 点干涉相位主要包含形变相位、大气相位及噪声相位［1］。雷达发射的电磁波在大气传输过程中，均会产生一定程度的折射，不同时刻所对应的大气条件发生变化，改变原有的传播方向与路径，产生大气相位延迟［2］。理论数据表明在温度为25 ℃，斜距1 km 处，仅1%的湿度变化可引起约1.3 mm形变测量误差。大气相位为GB-InSAR形变测量的主要误差源［3-4］。

针对大气相位的补偿，目前主要有GCP（Ground Control Points，地面控制点）补偿［5］、气象数据校正［6］，以及参数模型法［7］三种方法。第1 种方法主要通过人工在场景内布设角反射器作为GCP，通过估计GCP的大气相位实现监测目标点的相位补偿。在大气匀质性假设条件下，一般采用距离线性模型拟合各PS点的大气相位。对于较为复杂的监测场景，大气不满足匀质性假设时，采用高阶斜距模型或空间维插值实现区域内PS点的大气相位估计。2015年，汪学琴等学者［8］以不同分布的GCP剔除了目标点的大部分大气相位，与正垂线所得数据进行了比对分析，验证了该方法的有效性。2017 年，黄其欢等学者［9］利用GCP 推算出的大气扰动因子，根据距离加权获取所有待测点大气扰动相位，实现了大气相位补偿。考虑到传统的GCP补偿需要人工布设角反射点，依赖于布设点的数量与质量，常常面临场景中存在危险区域难以布设反射点的情况。第2种方法需要在观测场景内建立气象站，获得雷达数据对应时刻的气象参数（温度、湿度、大气压），基于大气折射率模型，定量求解大气相位，实现相位补偿。2008 年，Pipia 等［10］分析了大气对电磁波传输的影响，证实了气象数据的变化与干涉相位有直接关系。2014年，董杰等学者［6］使用角反射器作为监测目标，利用观测区域内的气象数据构建了大气折射模型，定量求解了大气扰动相位，修正了监测目标的干扰相位。气象数据校正依赖于监测区域中气象站的分布，数量较少的气象站难以准确获得各区域的气象数据，补偿效果较差。第3 种方法则是对PS 点的大气干涉相位数学建模，利用最小二乘回归估计出模型参数，实现相位误差的补偿。2005年，Noferini等［7］首次提出了基于PS技术，估计大气相位模型参数的方法。2014年，Iglesias等［11］分析了高程对于气象条件的影响，将高程引入大气相位参数模型，提高了测量精度。2016 年，徐亚明等［12］将该方法有效地应用于施工期高危边坡的监测中。当雷达监测区域过大或天气恶劣时，通常表现为较强的时空变性，部分干涉图无法建立多参数模型来有效估计大气相位，参数模型法误差较大，不再适用［13］。

因此，本文提出了基于PS 技术的改进方法，首先分析各邻近PS 对的差分相位序列标准差，实现PS 集合中噪声点以及形变点滤除，对剩余PS 采用MIC 结合空间维相干系数的双阈值实现稳定PS 的选取。利用稳定PS空间非均匀分布的特点，通过空间维插值拟合大气相位，从而实现大气相位的有效补偿。采用直线式扫描边坡监测雷达作为实验平台，对30 幅雷达图像进行分析，对比验证了本文方法在复杂大气条件下相位补偿的有效性。

2 补偿方法

2.1 参数模型补偿法

该方法基于PS 技术，利用雷达图像中PS 点估计模型参数。目前对于PS 点的选取广泛采用幅度离差法，通过评估图像中像素点的幅度稳定性来代替其相位稳定性，实现PS点的选择。

任意PS点的大气干涉相位Δφatm建模为［14］：

其中，λ为发射电磁波的波长，Δn(r，t)为与时间t和空间r相关的大气折射指数，L为电磁波传输路径。在大气空间同质性良好的前提下，Δn在空间r不发生变化，因此公式（1）可改写为：

式中，β0、β1分别表示线性方程的常量与线性相关系数，R表示发射天线与PS点的斜距。

基于所有PS点建立线性方程组

ΔΦ为n个PS 点解缠相位构成的n维向量，β为待估计线性系数所构成的2 维向量，ε为随机误差向量，表示一次模型的误差相位。利用最小二乘回归对β0、β1进行估计，得到

其中，T表示矩阵的转置。因此，线性模型估计的大气相位为：

由于部分PS 点中解缠相位存在噪声或变形的影响，导致模型与观测相位之间存在偏差。采用观测相位和估计相位差的无偏估计来去除有偏点，即不可靠PS点［15］。

通过式（6）准则剔除不可靠PS 点，再次迭代估计，直至不可靠PS点全部剔除，得到β0、β1准确估计。

ΔTatm取值一般在0.1～0.2 rad内。

除上述线性模型外，常用模型还包括方位角-斜距模型、二阶斜距模型等，如式（7）、（8）。

其中，β0、β1、β2为待估计的模型参数，R、θ分别为PS点的斜距与方位角。当大气变化较为复杂时，上述模型均无法做到有效补偿。

2.2 气象数据校正法

由式（1）可知，大气折射指数n仅与时间t和空间r相关。监测区域较小或天气良好时，大气随空间均匀变化，气象参数可视为保持一致。t时刻回波中的大气相位可以写为：

式中，φatm大气扰动相位只与距离、波长、大气折射指数有关。

根据大气折射经验模型，可得t时刻大气折射指数n(t)表达式为［16］：

式中，T为温度，单位为K。P为大气压强，单位为mbar。h为相对湿度，单位为%。

2.3 改进方法

2.3.1 噪点剔除

通过求解空间中相邻PS 点差分干涉相位序列标准差，设定合理阈值，剔除受噪声影响较大的PS点。

任意PS点差分干涉相位Δφ建模为：

其中，φdefo为形变相位分量，φatm为大气相位分量，φnoi为噪声相位分量。考虑到三个分量彼此独立，PS点的干涉相位序列标准差σInP可以表示为

其中，σdefo、σatm和σnoi分别对应各相位分量标准差。假定空间中相邻PS 点p1、p2，其差分相位Δφ(p1，p2)可表示为

在天气良好条件下，大气相位呈现大面积区域的高度空间相干性，即使复杂大气条件下存在空变性区域，局部空间中仍然呈现高度相干性，即趋近于0。而噪声相位在空间中则呈现无规则性，则临近PS 点对p1，p2 差分相位序列标准差可以表示为［17］

基于Delaunay 算法对PS 点构建三角形网络。任意一个PS 点均能由一个或多个三角形找到临近PS 点，定义为临近PS 点对。去除三角形中较长的边，且保证任意PS 点至少存在一个临近PS 点，通过计算临近PS点对差分相位序列标准差，设定合理阈值，去除与其相邻PS 点相位差异较大的点，即受噪声影响较大或部分发生形变的PS点。

2.3.2 基于双阈值的稳定点选择

由上述方法剔除受噪声影响较大的PS点后，剩余PS点干涉相位可表示为：

其中，φnoi为噪声相位分量，数值一般在0.1 rad 左右。若PS 点不发生形变，即为稳定PS 点，则φdefo趋近于0，干涉相位Δφ与大气相位φatm呈高度相关。通过计算PS 点大气相位序列与差分干涉相位序列的MIC，设定合理阈值，实现部分稳定PS点选择。

MIC，即最大互信息系数，用于衡量两个变量X和Y之间的相关程度，包括线性、非线性的强度。相较于Pearson 相关系数只能衡量变量线性相关程度的局限，MIC能够捕获到变量之间各种线性、非线性关系，对噪声不敏感，且具有更高的准确度，符合复杂大气条件下，各区域大气折射率呈现各种线性或非线性关系的实际情况，能实现空间大部分区域稳定PS点选择。其表达式为［18］：

其中，X、Y为离散随机变量，I[X；Y]为X、Y的互信息，表示为

式中，ρ(X)、ρ(Y)分别为X和Y的边缘概率，ρ(X，Y)是X和Y的联合概率。其计算方式通过将X、Y以散点图的形式离散在二维空间中，将当前空间沿X、Y方向上划分为|X|，|Y|个格子，以网格中当前散点个数的频率作为边缘概率与联合概率。B是变量，通常为n0.6，其中n为数据量的个数。

式中，X、Y分别为各PS 点的大气相位序列与差分干涉相位序列，通过计算MIC 并设定合理的阈值，实现部分稳定PS点的选择。由于大气常常表现为复杂时空变性，部分区域部分时间段内出现空变性，导致单一的MIC 阈值过低，出现部分区域缺少稳定PS点，对空间大气相位插值拟合过程中将会出现较大误差。

为此，文章引入相干系数，评估PS 点的稳定程度。PS 点主要有两种失相干因素，即时间失相干（大气改变）、空间失相干（位置改变）。部分区域在监测过程中的部分时间段呈现空变性，为避免时间失相干，通过求解多组连续获取的两张雷达图像的相干系数，设定较高阈值，实现部分空变区域稳定PS点选择。

相干系数法主要思路为在干涉相位对移动窗口并计算窗口内全体像元的相干信息作为当前窗口中心像元的相干系数值。其数学表达式为［19］：

其中S1，S2分别为干涉相位对的两个窗口，*指的是共轭相乘，m，n分别对应窗口的行、列。

通过设定MIC 与相干系数的合理阈值，求取并集，实现空间各区域范围内稳定PS 点选择，避免存在部分空变区域缺失稳定PS 点而引起较大的插值拟合误差。

2.3.3 Kriging插值拟合

选出的稳定PS 点在图像各区域中呈现非均匀分布，通过空间维插值拟合，实现所有PS 点大气相位估计。本文采用Kriging插值拟合。

Kriging 插值方法基于变异函数和结构分析理论。本文根据待求PS 点一定邻域范围内的稳定PS点大气相位，在考虑了其自身大小、形状以及与待求样点之间的空间位置信息，通过变异函数提供的结构信息对待求PS 点进行的一种最优的线性无偏估计。其数学表达式为［20］：

其中，φatm(si)为第i个稳定PS点处的大气相位，N表示稳定PS点个数。λi为对应样本点的未知权重，其满足估计值与真实值φatm(s0)的差值最小，即

其中，γ(h)为半变异函数，表示待求PS 点与稳定PS点距离的不同，影响程度不同。一般来说，在对气象要素场插值时半变异函数选用球状模型，其兼顾了储层参数的随机性与相关性。数学表达式为：

其中，C为拱高，C0为块金值，a为变程，表示空间中具有相关性的范围。

通过稳定PS 点空间位置以及差分干涉相位信息进行Kriging 插值拟合，实现所有PS 点大气相位的线性无偏最优估计。

3 实验信息

本次实验区域为重庆市万州区九道拐，所监测的山体边坡纵向高程约170 m、横向宽度约300 m，且地处长江北岸，时常受江风袭扰，该区域大气条件较为复杂［21］。

如图1（a）所示，采用实验平台为北京理工雷科公司研发的直线式扫描边坡监测雷达。该雷达工作波段为Ku，测量周期3～10 min，合成孔径约1.8 m，1 km 处空间分辨率为0.3 m×4 m，测量精度较高，达±0.1 mm。图1（b）、（c）所示分别为雷达监测区域与气象站。

图1 实验设备与场地Fig.1 Experimental equipments and venue

本次采用数据为边坡雷达连续获取的30 张雷达图像，时间从2020年12月18日10点40分～2020年12月18日15点13分，每幅图像测量周期约为9分钟。通过气象站获得该时间段内的气象数据（温度、湿度、大气压），其变化曲线如图2（a）、（b）、（c）所示。

图2 气象数据变化曲线图Fig.2 Meteorological data change curve

图3（a）所示为雷达图像，采用最大的幅值作为参考值进行了dB 处理。采用幅度离差法对30 幅图像进行PS 点选择，图3（b）为PS 点分布情况。通过设置幅度离差门限0.2，幅度门限−35 dB，共筛选出14842个PS点。

图3 雷达图像与PS点选择结果Fig.3 Radar image and PSs selection result

分析雷达图像时，采用连续获取的两张雷达图像作为主、副影像进行差分干涉，本次采用30 幅雷达图像前后差分干涉，得到29张干涉相位图。

短时间内，PS点一般不发生形变或形变量较小，即φdefo在0 rad 左右。如图4（a）、（b）所示，分别对应连续两幅雷达图像的干涉相位图与干涉相位散点图。

图4 相邻时间A组干涉相位图与散点图Fig.4 Interferometric phase diagram and scatter plot of group A at adjacent times

可以看出，PS点的干涉相位主要分布在-0.5～0.5 rad 范围内，随斜距发生线性变化。这说明此段时间范围内，大气变化较小，气象条件较为稳定，可以将大气延迟带来的相位误差建模为随斜距线性变化的分量，红色实线代表的线性大气相位估计可以起到很好的相位补偿作用。

图5（a）、（c）所示分别对应B、C 两组相邻时间干涉相位图。可以看出，干涉相位图中呈明显空变性且两组干涉相位散点图（b）、（d）中PS点的干涉相位较为分散，均随斜距发生了非线性改变，且变化趋势各不相同。说明在两组图像的获取时间范围内，大气具有较强的时变与空变性，无法通过多项式模型估计大气相位，常规方法不再适用。

图5 相邻时间B、C组干涉相位图与散点图Fig.5 Interference phase diagram and scatter diagram of two groups B and C at adjacent times

4 实验结果

4.1 噪点剔除

本次通过Delaunay三角网建立临近PS点对，以所有三角形边长的平均值长度作为参考，剔除超过该参考值（4.35 m）的长边，且保证任意PS点至少存在一个临近PS点。

根据式（14）计算出的所有PS点的临近标准差，若一个PS 点存在多个临近PS 点，则取其平均。如图6（a）所示，为所有PS 点临近标准差图。图6（b）所示为临近标准差散点图，由于随着斜距不断增加，PS 点信噪比逐渐降低，临近标准差也会相应增大。因此临近标准差采用随距离呈线性相关的阈值门限。本次在最短观测距离400 m采用0.35 rad，最远距离1100 m 处采用0.4 rad，如图中红色虚线所示。共筛选出597 个点，图6（c）中黑色像素点即为PS点中筛选出的噪点或部分发生形变的点，无规则分布在各个区域中，符合噪声特性。

图6 PS点临近标准差图与噪点选择结果Fig.6 PSs near standard deviation map and noise selection results

4.2 基于双阈值的稳定PS点选择

经过噪点剔除后，将剩余PS点差分干涉相位序列分别与对应斜距下的大气干涉相位序列进行MIC计算。如图7（a）所示为大气折射指数差分曲线，由式（10）计算得到，根据式（9）得到对应PS 点大气干涉相位序列。图7（b）所示为剩余PS 点MIC 图，设定阈值为0.92，得到与大气相位变化高度相关的稳定PS 点8436 个，其分布图如图7（c）所示。图中可以看出，稳定PS 点非均匀分布在图像中，但图中红色椭圆部分受大气空变性影响无法选出稳定PS点，在对其空间维插值补偿时，将会产生较大误差。

图7 大气折射指数差分曲线与高相关稳定PS点选择结果Fig.7 Atmospheric refractive index difference curve and high correlation stable PSs selection results

本文通过对29 张连续时间的干涉相位对的相干系数取平均，实现部分空变区域高相干稳定PS点的选择。图8（a）所示，为图像中所有像素点的相干系数图。设定阈值为0.99，在经噪点剔除后的剩余PS 点中共选取到高相干稳定PS 点4288 个，如图8（b）所示，图中可以看出在红色椭圆区域均选到稳定PS 点，弥补MIC 单阈值选点的不足。稳定PS点总数由双阈值取并集得到，共计10185 个，占比68.62%。图8（c）所示，即为本次选取的稳定PS 点分布图，稳定点非均匀分布在图像各个区域。

图8 相干系数图与稳定PS点选择结果Fig.8 Coherence coefficient map and stable PSs selection results

4.3 插值拟合结果

以本文所提方法筛选出稳定PS点后，利用其位置信息（X、Y坐标值）以及29 张差分干涉相位序列作为已知样本点，半变异函数选用球状模型，参考式（22），进行Kriging 插值拟合实现所有PS 点大气相位的线性无偏最优估计，从而进行相位补偿。

如图9（a）所示为本次对29 张干涉相位图补偿后的干涉相位累积图，相较于图9（b）常规线性补偿与图9（c）气象数据校正方法下的干涉相位累积图，所有PS 点补偿后的相位均在0 rad 左右，大气的空变性得到了有效改善。而两种常规方法补偿效果较差，噪点较多，且图中左侧PS 点干涉相位在-0.5 rad 左右，右侧点在0.5～1 rad 之间，此为29张干涉相位序列中大气相位补偿精度较低，残余大气相位不断叠加而导致。

图9 不同方法下的干涉图补偿结果Fig.9 Interferogram compensation results under different methods

为了更加直观地看出残余大气相位的叠加影响，在不同斜距、方位角处选择出幅度离差最小的PS点作为参考点，分析随时间序列的干涉相位变化图。本次共选取6 个参考点，其空间分布图如图10（a）所示。图10（b）为常规线性补偿后参考点的相位变化曲线，可以看出参考点相位明显偏离0 rad，且变化趋势各不相同。这说明常规线性补偿误差较大、精度较低，存在严重的过、欠补偿现象，由大气带来的相位误差随时间和空间剧烈变化，波动幅度在-1.5～2 rad 之间，约2.9 mm 形变误差。图10（c）为气象数据校正后参考点相位变化曲线，其波动较常规线性补偿更为剧烈，补偿效果最差，波动幅度在-1～3 rad 之间，约4.3 mm 形变误差。如图10（d）所示为本次改进方法下的参考点相位变化曲线，图中看出参考点相位波动均在0 rad 左右，且相位序列变化随机，不受大气相位叠加影响，补偿效果最优。其中波动原因主要来自于噪声，噪声相位一般在±0.1 rad之间，影响较小。

图10 不同方法下的参考点相位累积变化图Fig.10 The cumulative change diagram of the reference point phase under different methods

5 结论

GB-InSAR 形变监测中，大气变化所带来的相位延迟是主要的误差源之一。在复杂大气条件下，常规相位补偿方法不再适用，文章所提的一种基于PS技术的改进方法，解决了常规方法无法对大气波动剧烈时间段内的干涉相位图进行有效补偿的问题。通过对九道拐实测数据分析，比较了改进方法与常规方法下的参考点的相位变化曲线和累积相位图，验证了本文方法在复杂大气条件下相位补偿的有效性。

本文方法仍存在一些欠缺之处。首先在稳定PS点选择过程中，需要人为凭经验设定MIC 以及相干系数阈值；其次当雷达图像较少时，MIC无法有效捕捉到PS 点的大气相位序列与差分干涉相位序列之间的线性、非线性关系，导致稳定PS点的漏选；最后在插值拟合过程中，Kriging 插值算法复杂度较高，插值拟合效率较慢，对于应用于时序GB-InSAR图像实时处理中还需要进一步研究。