一种结合路径跟踪的多基线InSAR高程重建方法与流程

本发明涉及多基线insar高程重建技术领域，具体是一种结合路径跟踪的多基线insar高程重建方法。

背景技术：

干涉合成孔径雷达(syntheticapertureradarinterferometry，insar)可以高精度，高可靠性地获取地表三维信息和高程变化信息，被广泛应用于海洋监控、火山监测、地震检测和数字高程重建等领域。相位解缠又是insar技术应用的关键步骤，传统的单基线相位解缠技术受限于itoh相位连续性假设，通常难以有效获取不连续地形干涉图的真实干涉相位。有鉴于此，各种不受itoh相位连续性条件限制的多通道insar相位解缠技术相继被提出。

1994年，xu提出了基于中国余数定理(chineseremaindertheorem，crt)的多通道相位解缠技术，用多幅基线满足两两互质的干涉图构建同余方程组，利用crt求得唯一解，能有效地克服了相位模糊问题，可从理想状态的无噪声干涉图中获取精确的解缠相位。袁志辉等提出了闭环鲁棒的crt多通道相位解缠技术，改善了crt多通道相位解缠技术噪声鲁棒性，可有效地从噪声较小的多幅干涉图获取目标高程信息。ghiglia等提出了多基线最大似然函数(maximumlikelihood，ml)法，该方法利用真实相位与缠绕相位符合圆高斯分布的性质，使用贝叶斯条件概率对真实相位进行估计，在低噪声情况下能得到较好的解缠结果，但在噪声较大时，解缠结果出现大量“毛刺”点，算法鲁棒性较差。为了消除ml估计中的“毛刺”点，提高ml算法的鲁棒性，ferraiuolo，pascazio等人提出基于高斯马尔科夫随机场(markovrandomfields，mrfs)的最大后验概率(maximumaposteriori，map)估计方法，利用马尔科夫随机场统计分布模型来描述高程的先验分布，将相位解缠转化成为能量函数最小化的优化问题，再使用迭代条件模型(icm)进行求解。基于马尔科夫随机场的map算法，需要估计复杂的超参数，耗时较长，giampaolo提出一种全变分(totalvariation，tv)模型的map算法，避免了复杂的超参数估计，大大缩短了运行时间，但由于缺少局部信息，算法精度有所下降。

技术实现要素：

本发明的目的在于克服现有技术的不足，而提供一种结合路径跟踪的多基线insar高程重建方法，在路径跟踪策略的指导下通过边缘检测区分出连续与不连续区域，在连续区域则引入邻域约束，提高算法的抗噪性；在不连续区域则不引入邻域约束，避免邻域约束引起的误差传递现象。本发明将传统单基线算法的框架引入多基线算法中，提高了多基线算法的鲁棒性。本发明能够有效消除ml估计中的“毛刺”现象，重建高程精度高于map算法，且运行速度相比map算法得到大大提升。

实现本发明目的的技术方案是：

一种结合路径跟踪的多基线insar高程重建方法，包括如下步骤；

1)获取干涉相位质量图

获取干涉合成孔径雷达sar图像中最长基线对应的二维干涉图，即缠绕相位图，对二维干涉图进行预处理，获得干涉图像中反映干涉像元质量的像元的相位质量值，即干涉相位质量图，指导相位展开路径，干涉相位质量图由干涉图直接获得，相位质量值ρm,n的计算公式为：

公式(1)中，是以像素点(m,n)为中心的局部窗口内缠绕相位值，k×k为局部窗口的尺寸；

2)从步骤1)得到的干涉相位质量图中，选择相位质量值最大的非边界点，使用最大似然算法对该非边界点进行最大似然估计，得到估计高程值，作为路径跟踪策略的起始像素点，最大似然率和联合最大似然率fml(φ|h)的计算公式如下：

公式(2)和公式(3)中，κ为高度灵敏度，h为对应的地表高程值，代表第i幅干涉图的缠绕相位γi为第i幅干涉图的相干系数；

3)在步骤2)中获得已得到估计高程值的起始像素点，将起始像素点4邻域中未解缠的点，按相位质量从大到小排列作为邻接列；

4)取出邻接列中质量最好的点，在以该点为中心3*3的窗口内，计算该窗口内已解缠像素点间的被估计点的高程值若最大高程差hdiff小于阈值hth，则认为该点保持连续；否则认为该点为突变点，并将该点4邻域中未解缠的点加入邻接列；

5)判断邻接列是否为空，若非空，则继续执行步骤4)；若为空，则算法结束，得到重建的高程图。

步骤1)中，利用生成相位质量图，将单基线相位解缠算法中的路径跟踪策略引入到多基线相位解缠算法中。

步骤3)中，所述的4邻域中未解缠的点，是指该点上下左右四个相邻点。

步骤4)中，所述的被估计点的高程值计算公式为：

公式(4)和公式(5)中，β是归一化参数；y∈nx表示y为像素点x邻域的像素点集合，一般取8连通域；表示以x点为中心的3*3矩阵中已解缠像素点的平均高程值，是一个标量；表示待解缠像素点的可能高程值，是一个向量；阈值hth设置为最长基线对应的高度模糊数，公式(4)为一个迭代递推的估计运算。

与现有技术相比，如经典的ml算法，该算法利用公式(3)对地形进行高程重建，但在噪声较大时会产生大量的“毛刺”点。map算法则在ml算法的基础上引入了高斯马尔科夫随机场，利用高斯马尔科夫随机场统计分布模型作为先验函数将相位解缠问题转化为能量函数最小化的优化问题问题，由于要进行复杂的超参数估计，因此map算法运行十分缓慢。tv+map算法采用了全变分能量函数模型，避免了复杂超参数的估计，由于缺少局部信息算法精度不如map算法。为了提高多基线insar高程重建的精度并且保持良好的运算效率，本申请将传统单基线相位解缠技术框架引入到多基线相位解缠技术中。传统多基线相位解缠技术大多只考虑单个像素点的解缠，没有引入邻域约束从而导致算法抗噪性差，而传统单基线算法通过邻域约束可以获得很好的鲁棒性，但在不连续区域由于邻域约束导致严重的误差传递。

本发明的优点为：1、利用单基线相位解缠框架中的路径跟踪策略指导多基线高程重建，提高解缠质量的可靠性；2、利用边缘检测技术检测出不连续边界点，在连续区域引入邻域约束提高算法的鲁棒性，在不连续区域则不引入邻域约束，避免邻域约束引起的误差传递现象。3、有效的消除了ml算法估计中的毛刺现象，同时保证了运算效率。

附图说明

图1为本发明提出的一种结合路径跟踪的多基线insar高程重建方法实施流程图；

图2为路径跟踪策略生成邻接列示意图，黑色像素点表示已解缠点，灰色像素点表示已解缠点4邻域中未解缠像素点，由其构成邻接列；

图3(a)-图3(f)分别为实施例1高斯模拟的地形高程图、加入噪声的长基线干涉图、ml算法重建高程图、map算法重建高程图、tv+map重建高程图、本发明重建高程图；

图4(a)-图4(f)分别为实施例2美国isola公园真实地形高程图、加入噪声的长基线干涉图、ml算法重建高程图、map算法重建高程图、tv+map重建高程图、本发明重建高程图；

图5(a)-图5(c)分别为实施例1对高斯模拟地形中map算法重建高程误差、tv+map算法重建误差，以及本发明重建高程误差；

图6(a)-图6(c)分别为实施例2对美国isola公园真实地形中map算法重建高程误差、tv+map算法重建误差，以及本发明重建高程误差。

具体实施方式

下面结合附图和实施例对本发明内容做进一步阐述，但不是对本发明的限定。

实施例1：

针对多山地形图进行相位展开，图3(a)为高斯模拟地形高程，其横坐标为距离向像素，纵坐标为方位向像素，右侧彩条代表高程值，单位为米；图3(b)为多山地形场景的含噪声长基线缠绕相位图，纵横坐标和图3(a)相同，其信噪比为7.4db。

一种结合路径跟踪的多基线insar高程重建方法，如图1所示，包括如下步骤：

1)获取干涉相位质量图

对图3(b)二维干涉图进行处理获得干涉相位质量图，本实施例1的干涉相位质量图为相位的微分偏差图。

获取干涉合成孔径雷达sar图像中最长基线对应的二维干涉图，即缠绕相位图，对二维干涉图进行预处理，获得干涉图像中反映干涉像元质量的像元的相位质量值，即干涉相位质量图，指导相位展开路径，干涉相位质量图由干涉图直接获得，相位质量值ρm,n的计算公式为：

公式(1)中，是以像素点(m,n)为中心的局部窗口内缠绕相位值，k×k为局部窗口的尺寸；

2)从步骤1)得到的干涉相位质量图中，选择相位质量值最大的非边界点，使用最大似然算法对该非边界点进行最大似然估计，得到估计高程值，作为路径跟踪策略的起始像素点，最大似然率和联合最大似然率fml(φ|h)的计算公式如下：

公式(2)和公式(3)中，κ为高度灵敏度，h为对应的地表高程值，代表第i幅干涉图的缠绕相位γi为第i幅干涉图的相干系数；

3)在步骤2)中获得已得到估计高程值的起始像素点为中心，将起始像素点4邻域中未解缠的像素点标记为待解缠像素点，按相位质量从大到小排列作为邻接列，如图2所示；

4)取出邻接列中质量最好的点，在以该点为中心3*3的窗口内，计算该窗口内已解缠像素点间的被估计点的高程值若最大高程差hdiff小于阈值hth，则认为该点保持连续，采用公式(4)计算；否则认为该点为突变点，采用公式(5)计算，并将该点4邻域中未解缠的点加入邻接列；4邻域中未解缠的点是指该点上下左右四个相邻点；被估计点的高程值计算公式为：

步骤1)中，利用生成相位质量图，将单基线相位解缠算法中的路径跟踪策略引入到多基线相位解缠算法中。

5)判断邻接列是否为空，若非空，则继续执行第四步；若为空，则算法结束，得到重建的高程图。

图3(f)为本实施例结合路径跟踪的多基线insar高程重建方法对图3(a)的高程重建结果，示意图的横纵坐标与图3(a)相同，图3(c)-图3(e)分别为现有的ml估计算法、map估计算法和基于tv模型的map估计算法对图3(a)的高程估计结果，图5(a)-图5(c)分别为map算法重建高程误差、tv+map算法重建高程误差以及本发明重建高程误差。

由图3(c)-图3(f)以及图5可以看出本发明实施例方法对多基线insar高程重建具有良好的效果，能够有效去除ml算法估计存在的毛刺现象，高程重建精度高；下述表1列出了在实施例1中本发明与ml算法、map算法以及tv+map算法在同一matlab计算环境下进行高程重建随消耗的时间和归一化均方误差。

表1不同算法的运算效率及精度对比表

从表1可以看出本发明在多基线insar高程重建中运行时间远小于map算法和tv+map算法，运算效率与ml算法相当，且高程重建精度远高于上述三种算法。

实施例2：

针对具有突变边缘的美国isola公园真实地形进行高程重建，由于存在地形突变，传统单基线算法无法对该地形进行有效的解缠。

图4(a)为isola公园真实地形高程图，图4(b)为含噪声的长基线干涉相位图；采用与实施例1相同的方法，对该地形进行高程重建，重建结果如图4(f)所示，另外图4(c)-图4(e)分别为ml算法估计高程、map算法估计高程、tv+map算法估计高程，图6(a)-图6(c)分别为map算法重建高程误差、tv+map算法重建高程误差以及本发明重建高程误差，下述表2列出了在实施例2中本发明与ml算法、map算法以及tv+map算法在同一matlab计算环境下进行高程重建随消耗的时间和归一化均方误差。

表2不同算法运算效率及精度对比表

由表2可知，本发明不仅可以有效消除ml算法估计中的毛刺现象，还能有效检测出突变边界，从而避免了由于引入邻域约束导致的误差传递。在提升高程重建精度的同时保证了算法的运算效率。

本发明提出的结合路径跟踪的多基线insar高程重建方法是在路径跟踪策略的指导下通过边缘检测判断不连续边界从而引入邻域约束的的高程重建方法，在模拟地形与真实地形的高程重建实验中都证明了其有效性。

上述实施例，仅为对本发明的目的、技术方案和有益效果进一步详细说明的具体个例，本发明并非限定于此。凡在本发明公开范围之内做任何修改、等同替换、改进等，均包含在本发明的保护范围之内。