一种集成对流层大气延迟改正的时间序列InSAR分析方法与流程

本发明涉及雷达干涉测量技术领域，提出实现了一种面向地形起伏条件下滑坡形变监测的时间序列insar对流层大气延迟改正方法。

背景技术：

近三十年来，随着雷达卫星系统硬件的快速发展、观测数据的迅速积累和软件技术的日益进步，雷达遥感技术在对地观测领域取得了广泛的应用，其中，雷达干涉测量(insar)技术在大范围高精度地形测图中具有较大的优势；差分干涉测量(dinsar)技术能够准确提取地表微小形变信息，可用于监测地震、火山、滑坡等地质灾害。

然而，常规dinsar技术容易受到时空去相干和对流层大气延迟扰动的影响。为了解决这个难题，二十世纪初有学者提出了时间序列insar技术，通过分析相干点目标来避免时空去相干的影响。时间序列insar技术假设对流层大气延迟是空间相关和时间不相关的，采用时空组合滤波的方法来改正对流层大气延迟。

一般对于地形起伏的山区来说，重复轨道insar观测中的对流层大气延迟可分为湍流混合延迟分量和垂直分层延迟分量，其中，湍流混合延迟分量具有空间相关和时间不相关的统计分布特性；而垂直分层延迟分量与高程相关，并呈现与时间相关的季节性波动趋势。所以，传统时间序列insar技术中的时空组合滤波方法仅对湍流混合延迟分量有效，无法消除垂直分层延迟分量。后来，依据对流层延迟的时空分布特点和所用外部数据种类，有学者提出了滤波法、经验法和预测法三类方法估计和改正对流层延迟，这些方法各有优劣，但本领域尚未有效果足够理想的解决方案出现。

技术实现要素：

本发明目的在于克服现有技术不足，提出了一种集成对流层大气延迟改正的时间序列insar分析方法，在已有的时间序列insar数据处理技术流程中，加入对流层大气延迟改正的过程，分别改正垂直分层大气延迟和湍流混合延迟。其中，迭代线性相位-高程回归模型能够有效地避免发生形变的相干点的干扰；采用数值气象模型和迭代线性相位-高程回归模型联合估计垂直分层大气延迟，比仅用单一手段更加准确有效。

本发明提供的技术方案提供一种集成对流层大气延迟改正的时间序列insar分析方法，所述insar表示雷达干涉测量，其特征在于：在时间序列insar数据处理流程中，加入对流层大气延迟改正的过程，包括改正垂直分层大气延迟和湍流混合延迟，其中，采用数值气象模型和迭代线性相位-高程回归模型联合估计垂直分层大气延迟；所述迭代线性相位-高程回归模型，是采用迭代执行线性相位-高程回归模型，剔除发生形变的相干点，仅保留未发生形变的相干点来解算线性模型。

而且，实现过程包括以下步骤：

步骤1、对于包含n幅sar单视复数影像的数据集，选取其中一景为主影像，将其它n-1幅辅影像都与这幅主影像进行亚像元级精确配准；

步骤2、对这n-1幅辅影像进行重采样，与主影像进行共轭相乘生成干涉图，采用外部数字高程模型模拟干涉图的地形相位，并从干涉图中减去地形相位生成差分干涉图；

步骤3、进行相干点目标筛选；

步骤4、将选出的相干点目标连接成delauney三角网，采用三维相位解缠方法进行相位解缠，得到每个相干点的解缠后差分干涉相位；

步骤5、联合数值气象模型和迭代线性相位-高程模型估计垂直分层大气延迟；

步骤6、湍流混合延迟分量由时空组合滤波进行改正，包括在减去垂直分层延迟分量后，首先，对每个相干点的时间维差分相位序列，进行时间维的高通滤波，获得高通分量；然后，对每幅干涉图，进行空间维的低通滤波，获得低通分量，作为湍流混合延迟分量；

步骤7、改正垂直分层大气延迟和湍流混合延迟后，实现计算每个相干点的年平均线性形变速率和形变时间序列。

而且，步骤3中，相干点目标筛选的实现方式为，首先，使用幅度离差指数准则来选取初始候选点；然后，基于形变信号的空间相关性假设来分析每个候选点的相位稳定性，使用带通滤波器对相邻候选点的相位值进行滤波，得到残余相位噪声，计算相位稳定性指标gamma，进一步从候选点中筛选出最终的相干点目标。

而且，步骤5中，联合数值气象模型和迭代线性相位-高程模型估计垂直分层大气延迟，实现方式包括以下步骤，

步骤5.1、使用数值大气模型估计垂直分层大气延迟，得到垂直分层大气延迟l1；

步骤5.2、使用迭代线性相位-高程回归模型估计垂直分层大气延迟，得到垂直分层大气延迟l2；

步骤5.3、将数值气象模型和迭代线性相位-高程模型分别估计的垂直分层大气延迟进行融合，获取最终的垂直分层大气延迟，实现方式为，

针对某对差分干涉图，分别减去使用数值气象模型和迭代线性相位-高程模型估计获取的垂直分层大气延迟l1和l2，计算残余相位的均方根，保留均方根较小的相应垂直分层大气延迟作为最终的垂直分层大气延迟。

而且，步骤5.2中，使用迭代线性相位-高程回归模型估计垂直分层大气延迟，实现方式包括以下步骤，

步骤5.2.1、如果为第一次迭代，所有相干点上的垂直分层延迟相位设为零；否则使用上一次迭代得到的线性模型估计出的垂直分层延迟相位，从所有相干点的解缠相位中减去垂直分层延迟相位，计算临时线性形变速率；

步骤5.2.2、根据预设的形变速率阈值δrms筛选未发生形变的稳定相干点；

步骤5.2.3、基于步骤5.2.2所得的未发生形变的稳定相干点，利用最小二乘方法求解线性模型；

步骤5.2.4、计算所有相干点的垂直分层延迟相位；

步骤5.2.5，当未满足收敛条件时重复步骤5.2.1-5.2.4，直至迭代收敛，收敛条件为相邻两次迭代计算的延迟相位之差的均方根均值小于预设的迭代终止条件参数ε。

而且，步骤2中，预设的迭代终止条件参数ε设置为0.1弧度。

而且，用于改正地形陡峭山区地表形变监测中的对流层大气延迟相位误差。

本发明实现了时间序列insar分析方法中的对流层大气延迟改正，利用全球分布的数值气象模型数据容易获取的优点，迭代线性相位-高程回归模型不受发生形变的相干点目标影响，有效地削弱了大气延迟对山区地表形变信号提取的干扰。

附图说明

图1为本发明实施例的集成对流层大气延迟改正的时间序列insar技术整体流程图。

图2为本发明实施例的迭代线性相位-高程回归模型流程图。

图3为本发明实施例的迭代线性相位高程回归模型与数值气象模型联合估计垂直分层大气延迟。

具体实施方式

以下结合附图和实施例详细说明本发明的技术方案。然而，本领域的技术人员容易理解，实施例所描述的内容仅用于说明本发明，而不应当也不会限制权利要求书中所详细描述的本发明。

本发明提出一种集成对流层大气延迟改正的时间序列insar分析方法，用于改正地形陡峭山区地表形变监测中的对流层大气延迟相位误差。将对流层大气延迟分为垂直分层延迟分量和湍流混合延迟分量两项，首先，在已有线性相位-高程经验关系模型的基础上，发展一种迭代线性相位-高程回归模型，用于估计垂直分层延迟分量；然后，融合数值气象模型和迭代线性相位-高程模型，估计并改正垂直分层延迟分量；最后，采用时间序列insar分析方法中时间和空间滤波处理，改正湍流混合延迟分量。

参见图1，本发明实施例提出的一种集成对流层大气延迟改正的时间序列insar分析方法，包含以下步骤：

步骤1、对于包含n幅sar单视复数影像的数据集，选取其中一幅为主影像，选取标准为保证干涉对的时间基线和空间基线都尽量小，将其它n-1幅辅影像都与这幅主影像进行精确配准；其中，sar为合成孔径雷达，n为数据集中影像数量。

步骤2、对这n-1幅辅影像进行重采样，与主影像进行共轭相乘生成干涉图，采用外部数字高程模型(如srtm)模拟干涉图的地形相位，并从干涉图中减去地形相位生成差分干涉图。

步骤3、相干点目标筛选。

实施例中，相干点目标筛选实现如下：

首先使用幅度离差指数来选取初始候选点，某个像素的幅度离差指数da可表示为：

式中，μa和σa分别为幅度序列的平均值和标准差。当像素da小于某一个阈值时(优选为0.4)，该像素被选为初始候选点。

然后，基于形变信号的空间相关性假设来分析每个候选点的相位稳定性，使用带通滤波器对相邻候选点的相位值进行滤波，得到残余相位噪声。相位稳定性指标可以通过时间相干性gamma定义，其表达式如下：

式中，γ为时间相干性gamma，n为时序sar影像数量，为第i幅影像的残余相位噪声，e为数学常量，j为虚数单位符号。当某候选点的相干性gamma大于相应预设阈值时(优选为0.8)，该候选点被确定为最终的相干点目标。

步骤4、将选出的相干点目标连接成delauney三角网，采用三维相位解缠方法进行相位解缠，得到每个相干点的解缠差分干涉相位。三维相位解缠方法为现有技术，本发明不予赘述。

步骤5、联合数值气象模型和迭代线性相位-高程回归模型来估计垂直分层大气延迟。

步骤5.1、使用数值大气模型估计垂直分层大气延迟，得到垂直分层大气延迟：

基于数值大气模型的垂直分层大气延迟改正方法是利用全球气象再分析资料，包括温度、气压、湿度、位势等信息，根据大气物理模型公式计算sar影像采集时刻的天顶方向大气延迟量。insar差分干涉相位在空间和时间上都是相对量，需要将估计的垂直分层大气延迟在时间维和空间维做差。此外，insar信号传播路径为雷达视线方向(los)，不同于大气延迟的天顶方向，需要通过投影转换把天顶方向单路径相对延迟量转换为los方向双路径相对延迟相位。由于当前的全球气象再分析资料的空间分辨率较低，一般为几十公里，对于小范围的地表形变监测，数值气象模型方法无法估计小尺度的湍流混合延迟，所估计的对流层大气延迟仅为垂直分层大气延迟。

实施例中，基于数值气象模型的垂直分层大气延迟估计方法是利用全球气象再分析资料，包括温度、气压、湿度、位势等信息，根据公式(3)计算sar影像采集时刻的天顶方向大气延迟量。

式中，l为天顶方向大气延迟量，z为高程，z0为地表高程，p(z0)为表面气压，e为水汽分压，t为绝对温度，gm为重力加速度g在大气中的平均值，rd为干空气常数，rd＝287.05j/kg/k,rv为水汽常数，rv＝461.495j/kg/k。几个常数变量分别为：k1＝0.776k/pa,k2＝0.716k/pa,k3＝3.75×10³k²/pa。

本发明实施例使用的气象再分析资料是由欧洲中期天气预报中心(ecmwf)提供的全球均匀分布的格网化数值产品era5，空间分辨率约为30km，时间分辨率为1h，根据气压分层，每个格网点在垂直方向上有37层气象数据。由于era5的空间分辨率较低，对于小范围的地表形变监测，数值气象模型方法无法估计小尺度的湍流混合延迟，所估计的对流层大气延迟仅为垂直分层大气延迟。

insar差分干涉相位在空间和时间上都是相对量。假设形成干涉对的两幅影像a、b的采集时刻分别为t1和t2，干涉图上某两个像素a、b的高程分别为z1和z2，则这两个像素间的相对对流层延迟量δl(z1,z2,t1,t2)为：

δl(z1,z2,t1,t2)＝[l(z1,t1)-l(z1,t2)]-[l(z2,t1)-l(z2,t2)](4)

其中，l(z1,t1)为影像a上像素a的天顶方向大气延迟量，l(z1,t2)为影像a上像素a的天顶方向大气延迟量，l(z2,t1)为影像b上像素a的天顶方向大气延迟量，l(z2,t2)为影像b上像素b的天顶方向大气延迟量。

此外，insar信号传播路径为雷达视线方向(los)，不同于大气延迟的天顶方向，需要通过以下投影转换，把天顶方向单路径相对延迟量转换为los方向双路径相对延迟相位δφstrat(z1,z2,t1,t2)：

式中，θ为雷达传感器视角，λ为雷达系统采用的微波波长。

然后，使用迭代线性相位-高程回归模型估计垂直分层大气延迟。常规线性相位-高程回归模型利用大气延迟相位的时空分布特性，建立大气延迟相位与高程的线性关系模型，利用最小二乘求解与高程相关的垂直分层大气延迟分量δφstrat，

δφstrat＝kδφh+δφ0(6)

式中，h为高程，kδφ为转换系数，δφ0为全局常数偏移项。但是常规线性相位-高程回归模型的解算容易受到发生形变的相干点的影响。

步骤5.2、使用迭代线性相位-高程回归模型估计垂直分层大气延迟，得到垂直分层大气延迟：

常规线性相位-高程回归模型利用大气延迟相位的时空分布特性，建立大气延迟相位与高程的线性关系模型来估计求解与高程相关的垂直分层大气延迟分量，但是容易受到发生形变的相干点的影响。

本发明提出的迭代线性相位-高程回归模型(简称：迭代线性模型)在常规线性相位-高程回归模型(简称：线性模型)的基础上，采用迭代计算来排除发生形变的相干点对常规线性模型解算的影响，其核心思想是：剔除发生形变的相干点，仅保留未发生形变的相干点来解算线性模型。流程图如图2所示，具体迭代步骤如下：

步骤5.2.1、如果为第一次迭代，所有相干点上的垂直分层延迟相位δφstrat设为零；否则使用上一次迭代执行步骤5.2.4得到的线性模型估计出的垂直分层延迟相位，从所有相干点的解缠相位中减去垂直分层延迟相位，计算临时线性形变速率v；

步骤5.2.2、根据预设的形变速率阈值vthr来筛选未发生形变的相干点；

步骤5.2.3、基于本次迭代执行步骤5.2.2中筛选出的未发生形变的相干点，利用最小二乘方法求解公式(6)中的线性模型；

步骤5.2.4、计算所有相干点的垂直分层延迟相位δφstrat；

步骤5.2.5、重复步骤5.2.1-步骤5.2.4直至迭代收敛，收敛条件为相邻两次迭代计算的延迟相位之差的均方根均值δrms小于预设的迭代终止条件参数ε，即δrms<ε。

迭代线性回归模型中有两个参数的设置非常重要，第一个是步骤2中的线性形变速率阈值vthr将直接决定选出的非形变相干点数量，具体实施时该参数的选择可以综合考虑雷达波长、研究区域环境等因素进行设置，对于c波段雷达数据来说，优选设定δrms＝5mm/yr；第二个参数是步骤5中的迭代终止条件参数ε，迭代终止条件参数ε应取为一个非常小的正数，本实施例中延迟相位单位为弧度，参数ε取值优选设置为0.1弧度，更大的值会减少迭代次数，但是精度降低；更小的值会增加迭代次数，精度提升，但当ε减少到一定值时，精度改善不明显。

步骤5.3、将数值气象模型和迭代线性相位-高程模型分别估计的垂直分层大气延迟进行融合，获取最终的垂直分层大气延迟，包括针对某对差分干涉图，分别减使用数值气象模型和迭代线性相位-高程模型估计获取的垂直分层大气延迟，计算残余相位的均方根，保留均方根较小的相应垂直分层大气延迟作为最终的垂直分层大气延迟。实施例中的实现流程如图3所示。针对某一差分干涉图inf，减去数值气象模型估计的垂直分层大气延迟l1，计算残余相位的均方根δrms1；减去迭代线性模型估计的垂直分层大气延迟l2，计算残余相位的均方根δrms2。如果δrms1<δrms2，则保留数值气象模型估计的垂直分层大气延迟l1作为差分干涉图inf的垂直分层大气延迟，否则，保留迭代线性模型估计的垂直分层大气延迟l2作为差分干涉图inf的垂直分层大气延迟。将垂直分层大气延迟从差分干涉相位中减去。

步骤6、湍流混合延迟分量由时间和空间组合滤波进行改正。首先，对每个相干点的时间维差分相位序列，进行时间维的低通高斯滤波，从差分相位序列中减去低通滤波分量，获得高通滤波分量；然后，对每幅干涉图的所有相干点，进行空间维的低通滤波，获得低通分量，即为湍流混合延迟分量。将湍流混合延迟分量从差分干涉相位减去。

步骤7、在改正垂直分层大气延迟和湍流混合延迟后，从改正后的差分干涉序列计算每个相干点的年平均线性形变速率和形变量时间序列。基于改正后的差分干涉序列计算所需结果，可根据用户需求采用现有计算实现，本发明不予赘述。

具体实施时，以上流程可采用计算机软件方式实现自动运行。

本文中所描述的具体实施例仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。