地基SAR在时间去相关严重区域的形变监测方法

摘要

本发明一种地基SAR在时间去相关严重区域的形变监测方法，具体步骤包括：S1永久散射体点选取；S2三角网格生成；S3空时三维相位解缠；S4大气相位补偿；S5形变反演；S6孤立PS点对应形变量剔除；其中步骤S1的具体过程为：分别进行了幅度离差筛选、幅度门限筛选和相位信息筛选。本发明在传统PS点选取幅度离差法的基础上增加了幅度门限筛选和相位信息筛选以提炼出更高质量的PS点，从而使得形变监测具有更高的准确性。

说明书

技术领域

本发明属于合成孔径雷达技术领域，具体涉及一种地基SAR在时间去相关 严重区域的形变监测方法。

背景技术

形变检测在结构工程和大跨径工程比如桥梁、高塔、大坝等建筑领域中有 重要的应用，而且在山体滑坡和雪崩预警等场景中也有广泛的应用。

而地基SAR由于其配置灵活，便于运输，重访周期短等优点，可以在不同 的场景下对目标进行长时间的形变监测。在实际过程中，首先，将雷达传感器 放置在精密滑轨上，以固定时间为周期，进行匀速且完全重复的移动，以构成 合成孔径；然后，通过分析长时间序列中目标点的相位数据，获取其形变量信 息，通常采用的数据分析的方法为永久散射体（PS）算法，其主要步骤包括： 1.PS点选取；2.Delaunay三角网格生成；3.空时三维相位解缠；4.大气相 位补偿；5.形变反演和6.孤立PS点对应的型变量剔除。

在传统的PS算法中，首先要选取那些能够在长时间序列中保持相关性的点 进行分析，即选取散射性能稳定的点（PS点）进行分析，步骤1.PS点的选取 一般是采用幅度离差准则，即通过判断目标幅度在长时间序列中的稳定性来选 取PS点。幅度离差的计算公式如下：

$D_A=\frac{{\sigma }_A}{m_A}---\left(1\right)$

式（1）中，σ_A表示某点A在整个观测时间内的幅度标准差，m_A表示点A在整 个观测时间内的幅度均值。

幅度离差准则在长时间内相关性较好的区域是一个效果很好的判断准则，但 是在时间去相关严重的区域，由于数据信噪比较低且受到严重去相关的影响， 幅度离差法的使用会导致更多非PS点被选为PS点。

同时，传统PS算法中的步骤4.大气相位补偿，其采用的非线性形变相位模 型如下：

式（2）中n表示第n幅干涉图像，i表示第i个PS点，等号右边的第一项表示 两次观测之间目标形变引起的相位，第二、三项代表轨道、大气相位等误差因 素导致的相位，第四项代表散射等其他因素导致的相位。大气相位补偿这一步 主要是选用无形变的PS点，通过最小二乘法估算每幅干涉图像中的[A_n,B_n]两个 参数，但是由于某些PS点事实上存在形变，所以直接一次估计这两个参数可能 带来较大的误差，特别是在时间去相关严重的区域，PS点的分布较不均匀的情 况下。

因此，在时间去相关严重的区域，传统的PS算法的性能有所降低，其已无 法准确对形变进行监测。

发明内容

有鉴于此，本发明针对PS点选取技术的不足和大气相位补偿处理方法的缺 陷提出一种地基SAR在时间去相关严重区域的形变监测方法。

该方案主要是这样实现的：

一种地基SAR在时间去相关严重区域的形变监测方法，具体步骤包括：S1 永久散射体点选取；S2三角网格生成；S3空时三维相位解缠；S4大气相位补偿； S5形变反演；S6孤立PS点对应形变量剔除；

步骤S1的具体过程为：

步骤101、幅度离差筛选，即

对N幅地基SAR干涉图像中所有的PS点采用幅度离差准则进行筛选，得到 点集P₁；

步骤102、幅度门限筛选，即

对点集P₁中的每一PS点，计算N幅干涉图像中该PS点处信号的平均幅度， 选取平均幅度大于设定幅度门限的PS点构成点集P₂；

步骤103、相位信息筛选，即

（1）对干涉图像进行网格划分；

（2）在每一幅干涉图像中，对点集P₂中每一PS点处信号进行低通滤波处 理；

（3）针对每幅干涉图像中的每个网格，计算网格内属于点集P2的PS点处 信号干涉相位的平均值L_j,n，其中n=1,2…N，j=1,2…J，J为干涉图像所划分的 网格总数，L_j,n表示第n幅图像第j个网格对应的干涉相位的平均值；

（4）针对点集P₂中的每个PS点，用PS点处信号的干涉相位减去PS点 所在网格的干涉相位平均值L_j,n，得到干涉相位差值i=1,2…I，I表示 点集P₂中PS点的总数，表示第n幅干涉图像第i个PS点的干涉相位，表示第n幅干涉图像第i个PS点的干涉相位差值；

（5）、根据所述计算点集P₂中每个PS点的相关系数α_i，

其中j为虚数单位；

（6）选择相关系数α_i大于设定系数门限A的PS点构成点集P₃。

较佳的，本发明还包括一种通过迭代计算方法获得几乎无形变的PS点，从 而实现对大气相位进行精确估计的大气相位补偿方法，具体方案包括：

步骤201、设点集P₃中所有PS点的形变量都为零，根据退化后的无形变相位 模型：

采用最小二乘法估计大气参数[A_n,B_n]，并将估计的结果记为为第 n幅干涉图像中第i个PS点的相位值，r_i为第i个PS点到雷达的距离，ω_i,n为 相位测量随机误差量，在大量观测数据中可认为其平均值为零；

步骤202、针对点集P₃中的每一PS点，计算其对应的相关系数γ_i；

其中，real表示取实部；

步骤203、选取相关系数γ_i高于设定系数门限B的PS点构成点集P₄；

步骤204、针对步骤203得到的点集P₄，按照步骤201至步骤203的方式重 复执行，直至得到的点集所包含的PS点对应的相关系数γ_i都大于设定门限值C 或得到的点集包含的PS点的数量小于设定阈值，根据当前得到的点集计算的 记为最终的大气相位估计结果[A_n,B_n]。

有益效果：

第一、本发明通过对现有PS点选取方法进行改进，在传统PS点选取幅度 离差法的基础上增加了幅度门限筛选和相位信息筛选以提炼出更高质量的PS点， 从而使得改进后的PS点选取方法可以在时间去相关严重的区域选择出高质量的 PS点，使后续的PS点数据处理做好准备。

第二、本发明通过迭代的大气相位补偿方法估计大气相位参数，从而保证 本发明可以更加准确的估计大气相位参数[A_n,B_n]，为后续的形变反演等步骤提 供更准确的相位信息。

附图说明

图1为本发明形变监测方法的流程图；

图2为(a)场景照片(b)SAR图像(P1高塔、P2桥、P3房子群)；

图3为未处理PS点时的干涉相位图（时间基线为15min）；

图4为以幅度离差算法选出PS点的干涉相位（时间基线15min）；

图5为幅度离差法、幅度法和相位分析法共同选出PS点的干涉相位（时间 基线15min）；

图6为Delaunay三角网格图；

图7为大气相位去除前的干涉相位图（时间基线3h）；

图8为大气相位去除后的干涉相位图（时间基线3h）；

图9为PS点形变反演图（时间基线11h）；

图10为高塔、桥梁和房子群区域时间序列形变反演结果图。

具体实施方式

下面结合附图和实施例，对本发明方法进行详细说明。

本发明提供了一种地基SAR在时间去相关严重区域的形变监测方法，为了 在时间去相关严重的区域取得良好的监测效果，该方法对传统形变监测永久散 射体（PS）算法主要做了如下改进：在传统PS点选取幅度离差法的基础上增加 了幅度门限筛选和相位信息筛选两个步骤以提炼出更高质量的PS点，从而最大 限度的减少时间去相关严重的干扰；除此之外该方法对大气相位补偿采用迭代 的算法，从而挑选出几乎无形变的PS点用于估计大气相位参数，从而提高大气 相位估计的精度。

如图1所示，本发明形变监测方法的具体步骤包括：S1永久散射体点选取； S2三角网格生成；S3空时三维相位解缠；S4大气相位补偿；S5形变反演；S6 孤立PS点剔除；

S1永久散射体点选取；

步骤101、幅度离差筛选。

对地基SAR干涉图像中的所有PS点进行幅度离差计算，每一个PS点的幅 度离差计算公式如公式（1）所示，设定合适的门限值，选出那些幅度离差值小 于门限值的点，作为PS点的点集P₁。

图2（a）为场景的示意图，图2（b）为SAR图像，图3为未处理PS点时 的时间基线为15min的干涉相位图，可以看出，数据受到严重的去相关影响， 几乎什么也辨别不出来。在进行幅度离差筛选得到的点集P₁中的PS点以15min 为时间基线的干涉相位图如图4所示，可以看到图4有了更好的可辨别度，但 是栅瓣和旁瓣的影响还是清晰可见。这是由于幅度离差定义为幅值的标准差与 幅度均值的比值，主要衡量的是幅值数据的稳定度。但是，在实际处理数据过 程中，强目标的栅瓣和旁瓣也满足稳定度的条件，只是它们的幅度均值较小。 仅仅使用幅度离差准则，可能会将更多的栅瓣和旁瓣位置处的点选入PS点。所 需要的PS点的本来真实相位数据就被强目标的栅瓣或旁瓣的相位信息所掩盖了， 进而会影响后续的空时三维相位解缠和大气相位估计。所以要再对初始选取的 PS点进行幅度门限筛选。

步骤102、幅度门限筛选。

对所述点集P₁中的每一个PS点求取N幅干涉图像中该PS点处信号的平均 幅度，选出平均幅度大于设定幅度门限值的PS点构成点集P₂，在幅度门限值的 选取的过程中，一般先将幅值图像转化为dB图像，再以最亮点的-40dB为门限 进行筛选。

步骤103、相位信息筛选。

由于对信噪比较低的目标仅仅使用幅度离差法选取PS点还可能会将很多噪 声点选入，在一些时间去相关严重的区域比如植被丰富的区域这种现象尤为严 重。这时在幅度门限筛选PS点完成之后，要再结合目标的相位信息进行进一步 的筛选。根据干涉相位理论，干涉相位可以写成以下的形式：

分析式（6）中四项相位的相位特性可以得知，等号右边的前三项相位和满足一定的空间相关性，且空间临近的像素点相关性较强。而第 四项相位则是随机噪声，无空间相关性。这时可采用低通滤波的方法， 得到第四项随机噪声的信息。具体步骤如下：

（1）、对干涉图像进行网格划分，就相当于对干涉图像对应的整个照射场景 进行网格划分。

（2）对点集P₂中每一PS点处信号进行低通滤波处理，即在每一幅干涉图像 中，对点集P₂的每一PS点处信号进行低通滤波处理。

（3）针对每幅干涉图像中的每个网格,计算该网格内属于点集P₂的PS点处 信号干涉相位的平均值L_j,n，其中n=1,2…N，N为干涉图像的总数，j=1,2…J， J为干涉图像所划分的网格总数，L_j,n表示第n幅图像第j个网格对应的干涉相 位的平均值。

（4）对每幅干涉图像，针对点集P₂中的每个PS点，用该PS点处信号的干 涉相位值减去该干涉图像中该PS点所在网格的干涉相位平均值L_j,n，得到 即每个PS点由于随机噪声等空间无相关性的因素造成的相位差，和 的下标n和i分别表示第n幅干涉图样和第i个PS点；

（5）、根据得到的计算点集P₂中每个PS点的相关系数α_i，其定义如 式（3），

其中j为虚数单位。

（6）选择相关系数α_i大于设定系数门限值A的PS点构成点集P₃。在对PS 点进行选择时，其宗旨为点越多越好，精度越高越好；当门限值A选择较低时， 该步骤中选择的PS点较多，但此时所选PS点的精度就较低；经过多次实验， 在所需精度允许的范围内，将门限值A选择为0.97较为合适。

点集P₃中的PS点以15min为时间基线的干涉相位图如图5所示。可以看到 虽然仍有一些小的噪声点，但是干涉相位图的整体可辨别性已经大大提高。

S2建立Delaunay三角网络。

该步骤为现有技术，在此对其进行简单说明：

在选择完PS点后，处理的对象由规则网络变为一些空间离散分布的点。为 了后续的相位解缠处理，需要在这些离散点之间建立联系。这时使用Delaunay 三角算法以所有PS点为顶点，建立Delaunay三角网络，将所有的PS点进行连 接。由步骤S1选出的PS点所建立的三角网络如图6所示。

S3空时三维相位解缠。

该步骤为现有技术，在此对其进行简单说明：

建立Delaunay三角网络后，需要对所有PS点进行空时三维相位解缠。空时 三维相位解缠分为两个步骤：

I、时间一维相位解缠：由于整个观测时间很长，长时间间隔中，“边”上 的相位可能会超过模糊周期，而在实际过程中，只能获取其在[-π,π]的区间内 的值。因此，需要首先恢复“边”上的相位，即每幅图像中任意相邻两PS点的 相位差。时间维相位解缠的基本元素不是PS点，而是由Delaunay三角建立的“边”， 即由“边”所连接的两个端点的相位差。

II、空间二维相位解缠。由于误差等其他因素的影响，某些“边”的相位 在时间一维相位解缠后可能会出现整数周期的相位误差。这时，在空间维度上， 采用基于最小费用流（MCF）的相位解缠。最大限度地纠正时间解缠后的相位误 差。同时，“边”的长度作为输入参数，参与误差估计。一般假设较长的“边” 上的相位可能出现误差的概率较高，有更大的可能性需要被修正。

S4大气相位补偿。

三维相位解缠后，恢复了目标的真实相位。但是，由于还存在大气相位等 误差的影响，仍无法直接提取形变相位，因此需要估计然后消除大气相位误差。 传统非线性形变相位模型如式（2），一般选择在整个照射过程中无形变的点， 此时式（2）退化为式（4），然后采用最小二乘法估计大气参数[A_n,B_n]。但是传 统的算法无法保证选出的点在整个照射过程中无形变，所以会引入较大的误差， 这种误差在时间去相关严重的区域尤其严重。本步骤结合迭代算法进行处理。

具体步骤包括：

步骤201、设点集P₃中所有PS点的形变量都为零，根据退化后的无形变相位 模型,如式（4），

采用最小二乘法估计大气参数[A_n,B_n]，估计结果记为此时i=1,2,…I， I为点集P₃中包含PS点的个数。

步骤202、针对点集P₃中的每一PS点，计算其对应的相关系数γ_i,

其中，real表示取实部。

步骤203、判断第i个点的相关系数γ_i是否高于设定的系数门限B，若是， 则这个PS点被保留，其余的点被舍弃，所有被保留的PS点构成点集P₄；在实 际过程中，门限B一般选择为0.9，其中门限B的选取规则与门限A的选取规定 相似。

步骤204、对步骤203中选出的PS点重复步骤201至步骤203的操作（在 执行完步骤201至203后，相关系数会发生变化），直到所有的PS点的相关系 数γ_i都大于一个设定好的门限值C（0.9）或者得到的点集包含的PS点的数量小 于设定阈值（1000）后结束迭代计算，得到的点集记为点集P，该点集P计算得 到的记为最终的大气相位估计结果[A_n,B_n]。时间基线为3h，大气相位 去除前后的干涉相位图如图7图8所示。其中门限值C的选取规则与门限A的 选取规定相似。

S5形变反演。

将包含大气相位的式（4）估计出来，从式（2）中去除后，剩余的相位可以简单转化为视线方向的形变Δd，转变关系如下：

S6剔除孤立PS点对应的形变。

该步骤为现有技术，先对其进行简单说明：

在实际处理过程中，如果某些PS点（或PS点集）离其他所有PS点的距离 都较远，那么它们在相位解缠和大气相位补偿的过程中有更大的概率可能会出 现错误。为了尽可能剔除误差的影响，需要剔除这些孤立PS点所对应的形变。 为了寻找孤立点PS，先将所有PS点建立Delaunay三角网格，并通过计算机图 论中聚类的算法，分布集中的PS点划归为一类，判断为是某个类的条件是某个 PS点至少与该类中任意一PS点的距离小于某个门限值，进而按照距离将所有的 PS点划分为多个类。如果某一类中的PS点的数量低于门限值，这一类中所有 PS点将会被定义为孤立PS点，并最终剔除。

最终时间基线为11h的PS点形变反演结果如图9所示，对桥梁、高塔和房 子群区域的时间序列分析如图10所示。桥梁和房子区域形变反演量与理论值基 本吻合，整个观测时间内无形变，反演精度在0.3mm以内。对于高塔区域，由 于可能存在热胀冷缩的影响，在观测时间内存在一定的形变。

通过本实施例的结果显示改良之后的地基SAR永久散射体技术在时间去相 关严重的区域也可以进行高精度目标形变反演。可见利用本发明方法，在时间 去相关严重的区域能够精确处理地基SAR模式下的形变监测数据，通过数据处 理，可以反演出场景微小形变信息。为地基SAR形变监测提供技术支持。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步 详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于 限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等 同替换、改进等，均应包含在本发明的保护范围之内。