一种稳健的基于雷达分布式目标的地面沉降监测方法与流程

本发明属于城市地质灾害预测技术领域，具体涉及一种稳健的基于雷达分布式目标的地面沉降监测方法。

背景技术：

地面沉降是在自然因素和人为因素的共同作用下，所产生的一种地表高程缓慢降低的地质现象，是一种难以补偿的永久性环境和资源损失，并能够诱发一系列其它的环境灾害，其形成过程缓慢，导致短时间内不易被察觉，且微量的沉降需要靠精密仪器才能被发现，易被人们忽视，地面沉降在很多情况下是发生在人口密集、工业发达的城市，其危害面广，破坏严重，且沉降过程是不可逆的。

近年来，时序insar方法在沿海城市沉降监测中取得较好的应用效果。insar技术即以合成孔径雷达复数据提取的相位信息为信息源获取地表的三维信息和变化信息，通过两副天线同时观测(单轨模式)，或两次近平行的观测(重复轨道模式)，获取地表同一景观的复图像对。由于目标与两天线位置的几何关系，在复图像上产生了相位差，形成干涉条纹图。干涉条纹图中包含了斜距向上的点与两天线位置之差的精确信息。利用传感器高度、雷达波长、波束视向及天线基线距之间的几何关系，可以测量出图像上每一点的三维位置和变化信息。

常规时序insar方法如永久散射体合成孔径雷达干涉测量(ps-insar)对监测点目标的选取十分严格，要求在所涉及的干涉纹图中都保持较高相干特性，造成在围垦填海区提取到的高相干点目标稀疏，难以准确的反映地面沉降的时间或空间变化，限制了这些方法在沿海区域地面沉降监测中的应用。ferretti等在2011年提出了新一代分布式散射体干涉测量(ds-insar)时序分析算法—squeesar，可以对ps及ds目标进行联合处理，增加监测目标点密度，提高监测精度。

然而，squeesar求解非线性最优相位需要迭代，严重影响运算效率。为此，我们提出一种稳健的基于雷达分布式目标的地面沉降监测方法，采用基于ks检验和特征分解的ds-insar时序分析方法，构建合理的地面沉降预测模型，对地面沉降进行预测，在很大程度上提前发现异常地质现象，防止灾害发生，预防地面沉降造成的地质灾害进而减小经济损失。

技术实现要素：

本发明的目的在于提供一种稳健的基于雷达分布式目标的地面沉降监测方法，通过ks检验方法计算并选取合适的ds点目标，建立协方差矩阵统计干涉相位和相干系数信息，提取同分布目标和估计相关参数；通过建立的协方差矩阵中提取与最大特征值对应的特征向量的相位来估计最优相位；根据拟合优度检测回归方程的整体拟合度，用于确定最终的ds点目标；最终联合ps点目标进行时序insar处理，通过误差改正、基线精化和相位解缠得到最后形变速率和时间序列，进而更好的对地面沉降进行精确的预测，算法稳健可靠，以解决上述背景技术中提出现有技术中的问题。

为实现上述目的，本发明采用了如下技术方案：一种稳健的基于雷达分布式目标的地面沉降监测方法，具体包括以下步骤：

s1：sar时间序列影像的获得：利用合成孔径雷达全天时、全天候对地实施观测，按照时间序列获得多组单视复数影像(slc)，选取一组单视复数影像(slc)为主影像，剩余的各组单视复数影像(slc)为从影像，将各组从影像按照同名像素点的方式进行匹配，实现配准过程，获得sar时间序列影像；

s2：ds候选点的识别：

s201：定义目标函数：假设ap和aq分别为sar时间序列影像的p和q两个样本的幅度值，两者的表达式为：

①和②式中，f(a)是累积分布函数，当有n组sar时间序列影像数据，ap的经验累积密度分布函数(cdf)表达式为：

③式中，是一个指标函数，当ap,i≤a时为1，反之则为0；

s202：ks检验：定义p和q两个样本的经验分布函数之间的最大距离，两个样本ap和aq及其对应的经验累积密度分布函数(cdf)和ks检验统计量定义式为：

④式中，是距离集合的最大值，当n>8时检验统计量tks分布收敛于柯尔莫哥洛夫分布，给定显著性水平α，

当满足⑤式时，ap假设被拒，⑤式中kα是柯尔莫哥洛夫分布的α分位数，kα通过柯尔莫哥洛夫cdf计算，计算式为：

s203：ds候选点：当不满足⑤式时，ap和aq的假设满足要求，sar时间序列影像的p和q两个样本符合相同的概率分布函数，即p和q两个样本均为后向散射特性相同的sar时间序列影像像元点，识别为ds候选点；

s3：协方差矩阵的构建：步骤s203中，sar时间序列影像的某一像元点协方差矩阵的数学表达式为：

t＝e[y⑦

⑦式中，e[]表示为期望，y为像元时间序列上的标准化后的复数观测向量，y＝[y1，y2，y3….yn]，且e[|yj|²]＝1，在⑦的基础上，对步骤s203获得的ds候选点的相位信息进行自适应滤波，构建协方差矩阵，协方差矩阵数学表达式为：

⑧式中，非对角线上对应的元素的绝对值为空间自适应滤波后的相干性值，相位为空间自适应滤波后的干涉相位，获得更加精确的干涉相位和相干性；

s4：估算协方差矩阵最优相位值：对步骤s3中获得的协方差矩阵⑧进行特征分解，特征分解公式为：

⑨式中，λi为特征值，按降序排序λ1≥λ2≥,...,≥λn,ui为对应的特征向量，假设u1是最大特征分量，约束条件为则可表述为最大值约束问题，u1表述公式为：

特征分解的相位估算值表达式为：

式中，θm,n＝θm-θn，rm,n和ζm,n分别是协方差矩阵t和矩阵|u1||u1|^t中第m行n列的元素，符号代表哈达玛积，

按照上述分解方式，从协方差矩阵中提取与最大特征值对应的特征向量的相位来估计最优相位；

s5：确定最终的ds点目标：假设估算后的最优相位为θ＝[θ1，θ2，...，θn]，为第n景和第k景slc的对应的原始干涉相位，则

式中，r为给定的阈值，大于给定阈值的ds候选点被保留并视为最终的ds点目标；

s6：获得ps点目标：对永久散射体合成孔径雷达干涉测量，获得ps点目标；

s7：最终形变速率和时序形变序列的确定：对步骤s5中获得的ds点目标和ps点目标联合进行时序insar处理，经误差改正、基线精化和相位解缠得到最终形变速率和时序形变序列。

优选的，所述步骤s7中，误差改正包括地形改正和大气改正，地形改正的方式为：时间序列干涉相位与影像空间基线进行线性回归模型分析，其斜率即是地形改正，大气改正的方式：为通过对干涉相位在时间域和空间域进行不同的带通滤波消除大气相位，进行大气改正。

优选的，所述步骤s7中，基线精化的方式为：结合地面控制点，通过最小二乘等算法对空间基线进行精化估计，得到更为准确的空间基线。

优选的，所述步骤s7中，相位解缠的方式为：获得的干涉相位位于[-π，+π]之间，与地面真实相位之间相差nπ个整周模糊度，恢复干涉相位的整周模糊度，并生成一副连续的解缠相位图，实现相位解缠。

本发明的技术效果和优点：本发明提出的一种稳健的基于雷达分布式目标的地面沉降监测方法，与现有技术相比，具有以下优点：

本发明通过ks检验方法计算并选取合适的ds点目标，建立协方差矩阵统计干涉相位和相干系数信息，提取同分布目标和估计相关参数；通过建立的协方差矩阵中提取与最大特征值对应的特征向量的相位来估计最优相位；根据拟合优度检测回归方程的整体拟合度，用于确定最终的ds点目标；最终联合ps点目标进行时序insar处理，通过误差改正、基线精化和相位解缠得到最后形变速率和时间序列，进而更好的对地面沉降进行精确的预测，算法稳健可靠。

附图说明

图1为本发明的一种稳健的基于雷达分布式目标的地面沉降监测方法的流程图；

图2为本发明的一种稳健的基于雷达分布式目标的地面沉降监测方法s202中两样本ks检验示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

请参阅图1-2，本发明提供了一种稳健的基于雷达分布式目标的地面沉降监测方法，具体包括以下步骤：

s1：sar时间序列影像的获得：利用合成孔径雷达全天时、全天候对地实施观测，按照时间序列获得多组单视复数影像(slc)，选取一组单视复数影像(slc)为主影像，剩余的各组单视复数影像(slc)为从影像，将各组从影像按照同名像素点的方式进行匹配，实现配准过程，获得sar时间序列影像；

s2：ds候选点的识别：

s201：定义目标函数：假设ap和aq分别为sar时间序列影像的p和q两个样本的幅度值，两者的表达式为：

①和②式中，f(a)是累积分布函数，当有n组sar时间序列影像数据，ap的经验累积密度分布函数(cdf)表达式为：

③式中，是一个指标函数，当ap,i≤a时为1，反之则为0；

s202：ks检验：定义p和q两个样本的经验分布函数之间的最大距离，两个样本ap和aq及其对应的经验累积密度分布函数(cdf)和ks检验统计量定义式为：

④式中，是距离集合的最大值，当n>8时检验统计量tks分布收敛于柯尔莫哥洛夫分布，给定显著性水平α，

当满足⑤式时，ap假设被拒，⑤式中kα是柯尔莫哥洛夫分布的α分位数，kα通过柯尔莫哥洛夫cdf计算，计算式为：

s203：ds候选点：当不满足⑤式时，ap和aq的假设满足要求，sar时间序列影像的p和q两个样本符合相同的概率分布函数，即p和q两个样本均为后向散射特性相同的sar时间序列影像像元点，识别为ds候选点；

s3：协方差矩阵的构建：步骤s203中，sar时间序列影像的某一像元点协方差矩阵的数学表达式为：

t＝e[y⑦

⑦式中，e[]表示为期望，y为像元时间序列上的标准化后的复数观测向量，y＝[y1，y2，y3….yn]，且e[|yj|²]＝1，在⑦的基础上，对步骤s203获得的ds候选点的相位信息进行自适应滤波，构建协方差矩阵，协方差矩阵数学表达式为：

⑧式中，非对角线上对应的元素的绝对值为空间自适应滤波后的相干性值，相位为空间自适应滤波后的干涉相位，获得更加精确的干涉相位和相干性；

s4：估算协方差矩阵最优相位值：对步骤s3中获得的协方差矩阵⑧进行特征分解，特征分解公式为：

⑨式中，λi为特征值，按降序排序λ1≥λ2≥,...,≥λn,ui为对应的特征向量，假设u1是最大特征分量，约束条件为则可表述为最大值约束问题，u1表述公式为：

特征分解的相位估算值表达式为：

式中，θm,n＝θm-θn，rm,n和ζm,n分别是协方差矩阵t和矩阵|u1||u1|^t中第m行n列的元素，符号o代表哈达玛积，

按照上述分解方式，从协方差矩阵中提取与最大特征值对应的特征向量的相位来估计最优相位；

s5：确定最终的ds点目标：假设估算后的最优相位为θ＝[θ1，θ2，...，θn]，为第n景和第k景slc的对应的原始干涉相位，则

式中，r为给定的阈值，大于给定阈值的ds候选点被保留并视为最终的ds点目标；

s6：获得ps点目标：对永久散射体合成孔径雷达干涉测量，获得ps点目标；

s7：最终形变速率和时序形变序列的确定：对步骤s5中获得的ds点目标和ps点目标联合进行时序insar处理，经误差改正、基线精化和相位解缠得到最终形变速率和时序形变序列。

较佳地，步骤s7中，误差改正包括地形改正和大气改正，地形改正的方式为：时间序列干涉相位与影像空间基线进行线性回归模型分析，其斜率即是地形改正，大气改正的方式：为通过对干涉相位在时间域和空间域进行不同的带通滤波消除大气相位，进行大气改正。

较佳地，步骤s7中，基线精化的方式为：结合地面控制点，通过最小二乘等算法对空间基线进行精化估计，得到更为准确的空间基线。

较佳地，所述步骤s7中，相位解缠的方式为：获得的干涉相位位于[-π，+π]之间，与地面真实相位之间相差nπ个整周模糊度，恢复干涉相位的整周模糊度，并生成一副连续的解缠相位图，实现相位解缠。

综上所述：本发明通过ks检验方法计算并选取合适的ds点目标，建立协方差矩阵统计干涉相位和相干系数信息，提取同分布目标和估计相关参数；通过建立的协方差矩阵中提取与最大特征值对应的特征向量的相位来估计最优相位；根据拟合优度检测回归方程的整体拟合度，用于确定最终的ds点目标；最终联合ps点目标进行时序insar处理，通过误差改正、基线精化和相位解缠得到最后形变速率和时间序列，进而更好的对地面沉降进行精确的预测，算法稳健可靠。

最后应说明的是：以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。