一种基于迭代分解的大气信号改正方法

技术领域

本发明属于测绘科学与技术领域，涉及一种基于迭代分解的大气信号改正方法，主要应用于地基雷达干涉测量。

背景技术

地基合成孔径雷达(GBSAR)是一项相对较新的地面遥感技术，具备强大的地表形变监测能力。根据硬件设备的不同，GBSAR的观测范围有所不同。通常情况下，其观测距离可达到4-10公里，观测范围可达数平方公里到数十平方公里。在理想操作环境下，GBSAR可以达到亚毫米的测量精度。GBSAR可以探测微变信号，并可以广泛适用于自然目标以及人工建筑的变形监测任务。与其他遥感监测技术相比，GBSAR监测技术有其独特的优势。GBSAR形变监测依赖于InSAR(基雷达干涉测量)技术。众所周知，InSAR技术的形变监测精度受到大气对流层气象条件变化(尤其是大气湿度)的影响很大。即使观测发生的时间间隔很短，大气的影响也很显著。为了获得精确的测量位移，必须谨慎地改正大气变化，去除大气信号。

本节按时间顺序介绍GBSAR变形监测中与大气校正有关的、国内外的、具有代表性的研究工作以及技术，以及其中存在的不足。

早在2004年，Luzi等人使用毫米波传播模型，结合气象数据(湿度，温度和压力)，补偿大气变化(参见：Luzi,G.；Pieraccini,M.；Mecatti,D.；Noferini,L.；Guidi,G.；Moia,F.；Atzeni,C.Ground-based radar interferometry for landslides monitoring:Atmospheric and instrumental decorrelation sources on experimental data.IEEEtransactions on geoscience and remote sensing 2004,42,2454-2466)。他们指出，对于存在明显地形变化的监测区域，大气参数在空间上是不均衡的，从单个位置观测得到的气象数据无法应用于改正不同范围或坡度位置的大气变化。

Pipia等人于2006年提出了三种地基雷达InSAR大气改正方法(参见：Pipia,L.；Fabregas Canovas,F.J.；Aguasca Solé,A.；MallorquíFranquet,J.J.In A comparisonof different techniques for atmospheric artefact compensation in gbsardifferential acquisitions,IGARSS 2006:IEEE International Geoscience andRemote Sensing Symposium:31July-4 August 2006:Denver,Colorado,2006；IEEE:pp3722-3725)。其中前两种方法有一个相似的过程：即通过选择了多个稳定且高度相干的测量点，并使用它们的解缠相位来恢复大气相位和观测距离之间的线性回归模型。这两种方法的整体策略是相似的，唯一的不同之处在于相位解缠，其中一种使用一维解缠，另一种使用二维解缠。在这两种方法中，大气沿观测距离的线性变化关系，是在大气空间均匀性的假设下建立的，模型由斜率系数和偏移量组成。第三种方法是基于干涉图的光谱变化，利用二维Chirp-Z变换去除低频大气信号。这种方法的主要缺点是它需要一个足够宽，以及高度一致和稳定的区域。因此，这种方法取决于GBSAR监测范围的大小和数据采集之间的时间延迟。

Pipia等人在2008年使用连续零基线GBSAR极化测量法研究大气变化对位移测量的影响(参见：Pipia,L.；Fabregas,X.；Aguasca,A.；Lopez-Martinez,C.Atmosphericartifact compensation in ground-based dinsar applications.IEEE Geoscience andRemote Sensing Letters 2008,5,88-92)。该研究假定雷达信号通过均匀介质传播，设定不同观测距离的四个高度相干点确定为参考点，将这些参考点的实际相位的时间演化与模拟相位进行了比较。该研究发现，实际相位与模拟相位曲线，在变化形状、相对误差(即误差与距离比)、以及角度波动幅度方面具有强烈的相似性。基于这些发现，他们提出了一种基于相干性的方法来去除时间序列大气相位。在大气空间同质性的假设下，利用稳定的控制点来恢复大气与观测距离的线性回归模型。该方法并未考虑大气在高层方向的变化。

Iglesias等人在2014年指出在地形变化剧烈的山区，由于温度，压力和湿度等大气参数的波动，大气空间均匀性的假设不再满足，大气的折射率被认为是垂直分层变化的(参见：Iglesias,R.；Fabregas,X.；Aguasca,A.；Mallorqui,J.J.；Lopez-Martinez,C.；Gili,J.；Corominas,J.Atmospheric phase screen compensation in ground-based sarwith a multiple-regression model over mountainous regions.IEEE Transactionson Geoscience and Remote Sensing 2014,52,2436-2449)。因此，在地形陡峭的情况下，该研究将大气变化扩展为与观测距离和测量目标高度均相关的函数，通过选取一组高相干点进行校正。该方法被应用于改正零基准GBSAR数据在短期和长期InSAR处理中的大气变化。实际应用中，该方法无法改正由于大气湍流引起的空间上的长波大气信号。

Butt等人于2016年声称，可以将大气变化的时间序列作为无二阶平稳性的自相关变量(参见：Butt,J.；Wieser,A.；Conzett,S.In Intrinsic random functions formitigation of atmospheric effects in ground based radar interferometry,Proc.of JISDM,2016)。大气相位由两部分组成：大气相位、折射率和气象观测量之间的确定性关系；并从气象测量中推导出未模拟的残差相位与大气相位的随机关系。在模拟大气相位的确定性成分时，首先探测出观测区域内的永久性散射点(PS)，并将其作为稳定的控制点。PS点测量中通常的二阶平稳性假设被一个相对较弱的假设所取代：即一个k阶的固有随机函数。该研究基于该假设来模拟大气随机成分，并结合传统最小二乘估计方法来模拟大气中的确定性成分，从而实现大气改正。该方法认为大气信号在一个空间上是静态的，以及在时间维度上是自相关的，没有考虑由于大气湍流引起的空间上的动态长波大气信号。

综上所述，在GBSAR监测中，大气变化通常被视为一个静态信号，通过一个与观测距离和(或)高度相关的模型进行改正。该类模型回归方法在GBSAR滑坡监测中已被广泛使用，通常情况下，该类方法能够有效改正大气变化。但是，当强烈的动态大气湍流的发生情况下，该类方法会失效。实际中，当观测区域位于河流之上，或者观测距离相对较远(如接近或者达到10km时)时大气的动态变化经常很强烈，大气在空间上的变化难以模拟。

此外，在卫星InSAR领域，Yu等人在2017年提出使用一种迭代分解大气信号的方法，将卫星雷达干涉中的大气信号分解为高层相关部分和空间低频长波部分。但是该方法需要使用基于全球卫星定位系统数据(GPS)计算得出的大气总延迟或外部气象观测数据来估计估计高层相关的大气信号(参见：Yu C,Li Z,Penna NT.Interferometric syntheticaperture radar atmospheric correction using a GPS-based iterativetropospheric decomposition model.Remote Sensing of Environment.2018 Jan 1；204:109-21.)。对外部数据的依赖是该方法的主要缺陷，导致该方法在GBSAR干涉测量中并不能实现。因此，地基雷达干涉测量技术中特别需要一种能够兼顾静态大气信号与动态大气信号、且无需外部依赖数据的大气改正方法。

发明内容

本发明的目的是提供一种无需使用任何外部观测数据(如GPS或气象数据)、基于迭代分解的方法，同时改正地基雷达干涉测量中的静态大气信号和动态大气信号。为达到上述目的，本发明提供的一种基于迭代分解的大气信号改正方法，包括步骤：

步骤1，干涉图数据处理。

步骤2，雷达干涉相位中的大气信号总量分解为静态大气成分和动态大气成分。

步骤3，进一步将大气信号总量分为与观测距离和高程相关的大气成分和空间低频长波成分。

步骤4，最小二乘法估计与观测距离和高程相关的大气成分，实现第一步大气改正。

步骤5，二维高斯滤波模拟空间低频长波信号，实现第二部大气改正。

步骤6，重复步骤4和5直到迭代结束。

步骤7，改正最终估计得到的与观测距离和高程相关的大气成分和空间低频长波成分。

步骤1中所述干涉图数据处理，其特征在于，其进一步包括子步骤：对缠绕的干涉图进行滤波降噪；对滤波后干涉图进行相位解缠；以及干涉图相干性估计。

步骤2中所述雷达干涉相位中的大气信号总量分解为静态大气成分和动态大气成分，其特征在于：

2.1其中信号总量记为

2.2静态大气成分在空间中随观测距离R和相对高程H线性变化，即

2.3动态大气成分

步骤3中所述进一步将大气信号总量分为与观测距离和高程相关的大气成分和空间低频长波成分，其特征在于：

3.1静态大气成分

3.2本发明将上述大气信号总量进一步分解为与观测距离R和相对高程H的成分

步骤4中所述最小二乘法估计与观测距离和高程相关的大气成分，实现第一步大气改正，其特征在于，利用一定数量的相干点，通过相干点对应的解缠相位、观测距离、相对高程，使用最小二乘法估计与观测距离和高程相关大气成分

步骤5中所述二维高斯滤波模拟空间低频长波信号，实现第二部大气改正，其特征在于，从原有干涉图中减去每个点的与观测距离和高程相关大气成分

步骤6中所述重复步骤4和5直到迭代结束，其特征在于：

6.1空间低频长波成分

6.2迭代次数达到一定次数或者前后两次最小二乘估计与观测距离R和相对高程H的成分

步骤7中所述改正最终估计得到的与观测距离和高程相关的大气成分和空间低频长波成分，其特征在于：

7.1从干涉图中减去最终估计得到的成分

7.2重复该过程，可实现多幅干涉图的大气改正。

本发明将地基雷达干涉中的大气信号分解为与观测距离和相对高程相关的静态成分和由于局部气流变化引起的动态分量。本发明无需使用任何外部观测数据，仅使用干涉图上携带的相干点，首先估计出静态成分中与高程相关的分量c

和现有技术相比，本发明具有如下特点：

1、具有更广泛的适用性。通常地基雷达在近距离高频率观测，观测区域水汽变化相对稳定时，大气信号通常呈现为一个静态信号；在远距离观测或观测区域存在强烈的动态大气湍流的情况下，大气信号还包括一个额外的动态信号。本发明方法可以同时改正两种信号。

2、具备易实现性。本发明无需使用任何外部观测数据，如GPS、气象观测数据等，因而在实际测量中更灵活，更容易实现。

3、如果观测区域存在大面积形变(例如形变区域大小接近去除动态大气成分时候的截断波长)，与其他常规方法一样，本发明方法需要人工制作形变区域的掩膜，或者通过常规的InSAR层叠(stacking)方法确定的移动区域自动制作形变区域掩膜。

附图说明

图1本发明提出的基于迭代分解的大气信号改正方法流程图；

图2一幅滤波后缠绕干涉相位图；

图3缠绕的高程相关大气相位图；

图4缠绕的低频信号大气相位图；

图5缠绕的大气改正后干涉图；

图6金沙江白格滑坡2018年12月4日形变速率图；

图7金沙江白格滑坡2018年12月4日形变时间序列图。

具体实施方式

下面将结合附图进一步说明本发明的技术方案。图1显示了本发明提出的基于迭代分解的大气信号改正方法流程图，为实现本发明方法，具体实施主要包括以下步骤：

(1)本发明方法的处理的单个对象是一副缠绕的干涉相位图，如果有多幅干涉图，那对于所有干涉图的处理是各自独立的。本发明首先采用常规的滤波降噪方法对缠绕的干涉相位图进行滤波以降低噪声，图2显示了一幅滤波后缠绕干涉相位图。滤波后的缠绕干涉图，其相位值处于[-π,π]之间，本发明采用常规的二维或者三维相位解缠方法进行相位解缠，解缠后的干涉图相位值范围通常在缠绕范围之外。

(2)在生成缠绕相位干涉图时，本发明同时采用常规方法估计干涉图对应的相干性。该步骤与前一步骤是本发明方法的数据准备阶段。

(3)高程相关大气信号成分估计。本发明首先根据预定义的相干性阈值(通常为0.4–0.8之间)，结合预先定义的或者通过常规的InSAR层叠方法确定的移动区域掩膜，选择相干点。利用这些相干点的解缠相位、观测距离、相对高程，根据最小二乘方法，即可估计出与观测距离R和相对高程H的成分

(4)低频长波信号模拟。该步骤的处理对象是高程相关大气改正后干涉图

(5)判断迭代是否结束，即步骤(3)和(4)的重复次数是否达到设定的最大次数或者前后两次最小二乘估计

通过步骤(1)至步骤(5)，即可实现本发明提供的基于迭代分解的大气信号改正方法。此处给出在2018年12月4日白格滑坡中，使用GAMMA GPRI-II系统采集的61幅GBSAR影像。数据覆盖的时长为8.1小时，观测距离范围为6.5–9km。通过使用本发明提供的大气信号改正方法和常规的地基InSAR时序估计方法，得到观测区域的形变速率图(见图6：金沙江白格滑坡2018年12月4日形变速率图)和图6中标记的8个点的形变时间序列图(见图7：金沙江白格滑坡2018年12月4日形变时间序列图)。在该案例中，形变速率图和形变时间序列图中，负数值代表对应的点在靠近地基雷达，即该点有向下滑落的趋势。图7中，点1-4有表现为下滑，点5有泥土累积而表现微小的上升，点6-8表现为稳定。