一种基于堆排序的质量图相位解缠方法

技术领域

本发明属于InSAR图像处理技术领域，更为具体地讲，涉及一种基于堆排 序的质量图相位解缠方法。

背景技术

InSAR(interferometric Synthetic Aperture Radar，干涉合成孔径雷达)技术 最初的发展目的是为了进行地形测绘，最为成功的应用就是获取高精度DEM (Digital Elevation Model，数字高程模型)数据。而随着技术的发展，InSAR 应用领域在不断拓展，主要包括：地表形变探测、动目标检测、海洋测绘、森 林制图、洪涝监测、交通监测以及冰川研究等。星载InSAR系统是一个涉及领 域宽、层次多、高度复杂的大系统，其对应的处理环节也较为复杂。星载InSAR 系统数据处理主要步骤包括：InSAR复图像配准、干涉相位滤波、干涉相位解 缠以及DEM等。由于处理过程环环相扣，任意一个环节处理的误差都将直接被 带入下一步处理中去。

InSAR的基本原理就是通过两部天线同时观测或者经过两次的平行的观测， 获取地面呈现在SAR(Synthetic Aperture Radar，合成孔径雷达)图像上的相位 差，从而反演得到每点的高程值。而在实际的物理过程中，呈现在图像上的相 位差是真实相位差的卷叠。也就是说，实际上干涉图中和地面相关联的像点是 以为模2π的。为了能够反演得到高程，就必须对SAR图像的每一点加上2π的 整数倍，从而求得其真实的相位差值。因而相位解缠是InSAR成像处理的关键 步骤之一，其准确程度直接影响到DEM或者地表形变的精度。

当前，已有的干涉相位解缠方法主要分为两大类：路径跟踪法、最小范数 法。路径跟踪法的基本策略是选择合适的积分解缠路径，阻止相位误差传递。 这种算法的代表是Goldstein于1988年提出的枝切法。枝切法通过设置枝切线对 隔离的残差点进行解缠。此后，有许多学者改进了这种算法。Huntley利用最近 邻算法连接不同极性的残差点，从而得到最短长度的枝切线。Buckland利用整 数规划中的匈牙利算法连接正负残差点，获取配对枝切线的全局最小解。Flynn 提出了Mask-Cut算法，通过质量图引导设置枝切线，提高了枝切线设置的准确 度。Xu提出区域生长算法用于解缠，将相位区域按照质量高高低分开，并从高 质量的区域向低质量区域解缠。最小范数法实质上是将相位解缠的问题转化为 数学上的最小范数问题，通过寻找使相位梯度和缠绕相位的误差的L^p达到最小 来解缠相位。Takajo最先提出采用最小二乘法进行相位解缠，并利用FFT求解， 该方法没有引入质量图。随后，Ghigilia提出了利用质量图的最小二乘法，提高 了解缠精度。

值得注意的是，目前关于InSAR和D-InSAR(Differential InSAR差分干涉 合成孔径雷达)数据处理的相位解缠算法种类虽然较多，但是普适性较差。每 一种方法都有自身的局限性。枝切法解缠速度快，但是“枝切线”设置不当会导致 误差传播。最小二乘法容易产生全局性误差，且计算速度较慢。而由于噪声和 欠采样导致相位出现不一致的问题一直是相位解缠的难点。因而，必须对现有 的算法进行改进，找到既能满足精度要求、又能快速解缠的算法。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于堆排序的质量图相 位解缠方法，以实现解缠过程中兼顾速度与精度的目的。

为实现上述发明目的，本发明一种基于堆排序的质量图相位解缠方法，其 特征在于，包括以下步骤：

(1)、SAR图像预处理

将两幅SAR图像共轭相乘，得到干涉图像，再对干涉图像依次进行去平地 效应和中值滤波，得到待解缠相位图像A；

(2)、从待解缠相位图像A中提取相位质量图P

(2.1)、根据待解缠的相位图像A，分别获取到行缠绕相位梯度图像矩阵 Ax和列缠绕相位梯度图像矩阵Ay；

(2.2)、遍历行缠绕相位梯度图像矩阵Ax中所有点，计算出待解缠相位图 像A距离向的每点的相位质量m_i,j；

其中，γ表示窗口的大小；符号表示取整数；Ax_i,j表示在行缠绕相位梯 度图像矩阵Ax中，第i行、第j列个点的像素值，i＝1,2,…,m，m为行缠绕相位 梯度图像矩阵Ax的行数，j＝1,2,…,n，n为行缠绕相位梯度图像矩阵Ax的列数；

再将所有点的相位质量m_i,j构成行缠绕相位质量图Ax^*；

(2.3)、遍历列缠绕相位梯度图像矩阵Ay中所有点，计算出待解缠相位图 像A方位向的每点的相位质量n_i,j；

其中，Ay_i,j表示在列缠绕相位梯度图像矩阵Ay中，第i行、第j列个点的 像素值；

再将所有点的相位质量n_i,j构成列缠绕相位质量图Ay^*；

(2.4)、计算相位质量图P

P＝α₁Ax^*+α₂Ay^*；

其中，α₁和α₂为影响因子，且满足α₁+α₂＝1；

(3)、建立最大堆

(3.1)、在待解缠图像A中任选一点作为解缠的起始点，并将该起始点标 记为已解缠点Q；

(3.2)、以已解缠点Q为中心，寻找距离已解缠点Q最近的点，即Q点的 邻域，再判断该邻域中某一点是否解缠，如果该点解缠，则舍弃该点；如果该 点没有解缠，则放入堆栈中，直到判断完邻域中的所有点；

(3.3)、利用计算机科学堆排序法对堆栈中所有点按照质量图P中质量大小 生成完全二叉树，建立最大堆；

(4)、对最大堆的根节点进行解缠

(4.1)、将最大堆的根节点与最后一个叶节点交换，同时将该根节点移出最 大堆，并对该根节点在待解缠图像A中对应的点进行相位解缠；

(4.2)、将最大堆中剩余的节点保留在堆栈中，再以步骤(4.1)中最大堆 的根节点作为步骤(3.1)中选取的解缠起始点，并返回到步骤(3.2)，重复上 述步骤，直到待解缠图像A中所有的点解缠完成。

本发明的发明目的是这样实现的：

本发明基于堆排序的质量图相位解缠方法，先通过去平地效应以及滤波得 到待解缠相位图像，然后从待解缠相位图像中提取相位质量图，并使用堆排序 法对相位质量图进行排序，即确定出待解缠相位图像中待解缠点的解缠顺利， 最后，根据解缠顺序对待解缠点解缠，直至整幅待解缠相位图像解缠完毕。本 发明利用待解缠相位图像的质量信息引导像待解缠点进行积分解缠，同时又添 加堆排序的法对相位质量图进行排序，从而缩短相位解缠所需时间，最终快速、 高精度的完成相位解缠，生成相位解缠图像。

同时，本发明基于堆排序的质量图相位解缠方法还具有以下有益效果：

(1)、通过方位向和距离向影响因子的引入，能够在确定相位质量图时更 好的与实际成像过程中距离向和方位向的噪声影响相匹配，从而能够更好的指 导相位解缠，提升相位解缠的准确度；

(2)、堆排序算法是一种稳定的计算机排序算法；堆排序的方法与具体的 数据无关，因而在获取SAR图像的相位质量图后，不必担心相位质量图数据的 不规则情况；虽然堆排序在数据小的情况下优势不明显，但是却非常适用于较 大数据的排序，而SAR图像基本都是数据量极大的图像，因而使用堆排序的方 法能够更好的提升SAR图像的处理速度；

(3)、由于最大堆满足根节点大于左右叶子节点的性质，因而其根节点始 终是最大值，每次选取相位质量最高点时只需要取出根节点即可，并且，由于 堆排序自身的特点，在每次加入新的邻域点进入堆栈中时，都能够很快的生成 最大堆。

附图说明

图1是本发明基于堆排序的质量图相位解缠方法的流程图；

图2是对伊朗巴姆地区进行预处理后的待解缠相位图像；

图3是求取行缠绕相位梯度图像矩阵的说明图；

图4是待解缠相位图像中距离向的相位质量计算示意图；

图5是从伊朗巴姆地区的待解缠相位图像中提取的质量图；

图6是选取解缠的起始点的示意图；

图7是以已知点邻域四个点为例的说明图；

图8是对图7中四个邻域点的最大堆排序过程示意图；

图9是加入新邻域点并且建立最大堆的示意图；

图10是图2所示基于质量图堆排序方法的相位解缠结果图；

图11是菲尼克斯亚利桑那州凤凰城的两幅SAR图像截图；

图12是图11所示的SAR图像经过共轭相乘后得到的的干涉图；

图13是图12所示的干涉图经过处理得到的质量图；

图14是图13所示的质量图经过相位解缠后的结图；

图15是菲尼克斯亚利桑那州凤凰城地区的高程图。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更 好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设 计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

图1是本发明基于堆排序的质量图相位解缠方法的流程图。

在本实施例中，如图1所示，本发明基于堆排序的质量图相位解缠方法， 包括以下步骤：

S1、获取干涉图像

将两幅SAR图像共轭相乘，即分别将两幅SAR图像对应的像素共轭相乘， 得到干涉图像；

在本实施例中，通过星载SAR对同一目标区域成像两次，得到两幅SAR 图像。

S2、对干涉图像进行去平地效应，即得到相干图像；

S3、对相干图像进行中值滤波，得到待解缠相位图像A；

本实施例中，在伊朗巴姆地区(BAM)，对尺寸为2823×1973的ASAR数 据图经过预处理后，即依次经过步骤S1～步骤S3的处理，仍有干涉条纹存在， 如图2所示，这种干涉条纹产生的原因即为相位缠绕的结果。为了能够进行高 程差的重建提取，就必须对此图像进行相位解缠，恢复真实的相位差。

S4、从待解缠相位图像A中提取相位质量图P

S4.1)、根据待解缠的相位图像A，分别获取到行缠绕相位梯度图像矩阵Ax 和列缠绕相位梯度图像矩阵Ay；

具体方法为：

将待解缠相位图像A中所有像素点的像素值从[0,2π]映射到[-0.5，+0.5]， 的到映射图像A₁；

遍历映射图像A₁，求取映射图像A₁中每行内相邻两个像素点的像素差值 Δ(i,j)，即：

Δ(i,j)＝A(i,j+1)-A(i,j)

判断Δ(i,j)是否属于[-0.5,+0.5]；

${\Delta }^{\prime }(i,j)=\left(\begin{array}{c}\Delta (i,j)+1,\Delta (i,j)<-0.5\\ \Delta (i,j),-0.5\le \Delta (i,j)\le 0.5\\ \Delta (i,j)-1,\Delta (i,j)>0.5\end{array}\right)$

其中，i为行坐标，i＝1,2,…,m，m为待解缠图像的行数；j为列坐标， j＝1,2,…,n，n为待解缠图像的列数；

在本实施例中，如图3所示，遍历此图像，计算Δ'(i,j)。

a)、计算Δ'(1,1)：先计算出Δ(1,1)＝0.2-0.1＝0.1，因为-0.5≤0.1≤0.5，因 此Δ'(1,1)＝Δ(1,1)＝0.1。

b)、计算Δ'(1,2)：先计算出Δ(1,2)＝-0.4-0.2＝-0.6，因为-0.6＜-0.5，因 此Δ'(1,2)＝Δ(1,2)+1＝0.4。

c)、计算Δ'(1,3)：先计算出Δ(1,3)＝0.5-(-0.4)＝0.9，因为0.9≥0.5，因此 Δ'(1,3)＝Δ(1,3)-1＝-0.1。

将所有的像素差值Δ'(i,j)构成m行n-1列的矩阵，再将矩阵中第n-1列的值 赋给第n列，构成行缠绕相位梯度图像矩阵Ax；

获取列缠绕相位梯度图像矩阵Ay的方法为：

遍历映射图像A₁，求取映射图像A₁中每列内相邻两个像素点的像素差值 即：

$\tilde{\Delta }(i,j)=A(i+1,j)-A(i,j)$

判断是否属于[-0.5,+0.5]；

${\tilde{\Delta }}^{\prime }(i,j)=\left(\begin{array}{c}\tilde{\Delta }(i,j)+1,\Delta (i,j)<-0.5\\ \tilde{\Delta }(i,j),-0.5\le \Delta (i,j)\le 0.5\\ \tilde{\Delta }(i,j)-1,\Delta (i,j)>0.5\end{array}\right)$

此处的像素点的像素差值的计算方法与行缠绕相位梯度图像矩阵相同，仅 仅是求每列中像素点的差值，因此不再赘述。

将所有的像素差值构成m-1行n列的矩阵，再将矩阵中第m-1行的值 赋给第m行，构成列缠绕相位梯度图像矩阵Ay。

S4.2)、遍历行缠绕相位梯度图像矩阵Ax中所有点，计算出待解缠相位图 像A距离向的每点的相位质量m_i,j；

其中，γ表示窗口的大小；符号表示取整数；Ax_i,j表示在行缠绕相位梯 度图像矩阵Ax中，第i行、第j列个点的像素值，i＝1,2，…,m，m为行缠绕相位 梯度图像矩阵Ax的行数，j＝1,2,…,n，n为行缠绕相位梯度图像矩阵Ax的列数；

再将所有点的相位质量m_i,j构成行缠绕相位质量图Ax^*；

在本实施例中，如图4所示，每个方格代表矩阵Ax中的一个点，实线方格 构成矩阵Ax，虚线方格是矩阵Ax边界点外补充的点，且虚线方格对应的数值 为0。

下面以点1、点2和点8三个点分为三种情况为例加以说明。

a)、求取边界点：方格1(即第1行第1列)的距离向相位质量，则以方格 1为中心，γ＝3*3时，周边包含五个虚线方格和2、7、8三个实线方格，利用上 述的计算公式可以计算方格1在3*3窗口中的相位质量m_1,1。

b)、求取边界点：方格2(即第1行第2列)的距离向相位质量，则以方格 2为中心，γ＝3*3时，周边包含三个虚线方格和1、3、7、8、9五个实线方格， 同样可以计算出方格2在3*3窗口中的相位质量的m_1,2。

c)、求取非边界点：方格8(即第2行第2列)的距离向相位质量，则以方 格8为中心，γ＝3*3时，周边只包含1、2、3、7、9、13、14、15八个实线方 格，同样也可以计算出方格8在3*3窗口中的相位质量的m_2,2。

S4.3)、遍历列缠绕相位梯度图像矩阵Ay中所有点，计算出待解缠相位图 像A方位向的每点的相位质量n_i,j；

其中，Ay_i,j表示在列缠绕相位梯度图像矩阵Ay中，第i行、第j列个点的 像素值；在此，方位向每个点的相位质量n_i,j的计算方法与距离向计算相同，在 此不再赘述；

再将所有点的相位质量n_i,j构成列缠绕相位质量图Ay^*；

S4.4)、计算出相位质量图P

P＝α₁Ax^*+α₂Ay^*；

其中，α₁和α₂为影响因子，且满足α₁+α₂＝1；在本实施例中，α₁＝α₂＝0.5， 如图5所示，对BAM地区的图像经过上述步骤处理后，即可得到相位质量图P；

S5、确定待解缠图像A中的解缠的起始点

在待解缠图像A中任选一点作为解缠的起始点，该起始点尽量位于相位质 量高的区域，并将该起始点标记为已解缠点Q；

S6、将Q点邻域中的未解缠点放入堆栈

以已解缠点Q为中心，寻找距离已解缠点Q最近的点，即Q点的邻域，再 判断该邻域中某一点是否解缠，如果该点解缠，则舍弃该点；如果该点没有解 缠，则放入堆栈中，直到判断完邻域中的所有点；

在本实施例中，如图6所示，每个正方形的小方格代表图像中的一个点； 图像中已经被解缠的点标记浅灰色，未被解缠的点标记深灰色，选择其中质量 最大的点作为解缠的起始点，并将此点标记为浅灰色，即图5中箭头所指的点；

S7、建立最大堆

利用计算机科学堆排序法对堆栈中所有点按照质量图P中质量大小生成完 全二叉树，建立最大堆；

下面对其具体步骤描述如下：

S7.1)建立最大堆

设堆栈中邻接列的数组为{r₁,r₂,...,r_n}，n为数组中元素的个数，最大堆需满 足如下性质：

$\left(\begin{array}{c}{r}_i\le r_{2i}\\ r_i\le r_{2i+1}\end{array}\right)$或$\left(\begin{array}{c}{r}_i\ge r_{2i}\\ r_i\ge r_{2i+1}\end{array}\right);$其中$i=\mathrm{1,2},...,\frac{n}{2};$

在本实施例中，邻接列的数组为：a[4]＝{0.89,0.93,0.81,0.85}，其建立最大 堆步骤为：

S7.1.2)、根据数组a[4]＝{0.89,0.93,0.81,0.85}，随机生成一个二叉树；

在本实施例中，如图7(c)所示，以起始点邻域中四个点为例，随机生成 一个二叉树，二叉树结构如图7(a)所示，那么这四个点在计算中的存储结构， 即邻接列，如图7(b)所示；

S7.2.1)、从二叉树的最后一个非叶子节点开始调整；其调整规则为：取此 非叶子节点、左子节点和右子节点中最大值点与此非叶子节点相交换，直到所 有的非叶子节点都满足最大堆的定义为止，即生成完全二叉树；

在本实施例中，如图8所示，选取8(a)中的最后一个非叶子节点为0.93， 其左子节点为0.85，无右子节点。取其中最大值点与该非叶子节点交换后，仍 为8(a)所示的二叉树。继续向前寻找非叶子节点为0.89，其左子节点为0.93， 右子节点为0.81，取其中最大值点0.93与该非叶子节点交换后，生成完全二叉 树，如图8(b)所示；

S7.2)、如图8(c)所示，将最大堆的根节点0.93与最后一个叶节点0.85 相交换，并将根节点0.93移出最大堆；

S7.3)、如图8(d)所示，将剩余的节点保留在堆栈中；

S7.4)、如图9(a)所示，当新的节点加入堆栈时，即节点0.9、0.7、0.6加 入堆栈时，则并重复步骤S7.1)～步骤S7.3)，直到步骤S7.3)中没有剩余节点 为止，如图9(b)所示；

S8、对最大堆的根节点进行解缠

S8.1)、将最大堆的根节点与最后一个叶节点交换，同时将该根节点移出最 大堆，并对该根节点在待解缠图像A中对应的点进行相位解缠；

S8.2)、将最大堆中剩余的节点保留在堆栈中，再以步骤S8.1)中最大堆的 根节点作为步骤S5中选取的解缠起始点，并返回到步骤S5，重复上述步骤，直 到待解缠图像A中所有的点解缠完成。

在本实施例中，以图2所示的待解缠相位图经过本发明的相位解缠处理后， 其结果如图10所示，图像中已经无明显的干涉条纹存在。

实例

一般而言，本发明需要依据地形没有发生变化的两幅SAR数据进行一系列 的处理得到的。而由于伊朗巴姆地区获取的两幅SAR图像分别是在地震前后， 因而其解缠结果不能够与实际的地形相比较。为此，需要找寻更为合适的星载 数据。

图11是菲尼克斯亚利桑那州凤凰城的两幅SAR图像截图；

图12是图11所示的SAR图像经过共轭相乘后得到的的干涉图；

图13是图12所示的干涉图经过处理得到的质量图；

图14是图13所示的质量图经过相位解缠后的结果图；

图15是菲尼克斯亚利桑那州凤凰城地区的高程图。

在本实施例中，如图11所示，该地区的两幅RADARSAT-2图像分别获取于 2008年5月4号和2008年5月28号，图像尺寸为11471×14206。并且该地区 有三座山，能够更直观的验证，如图15所示。

好的解缠结果应该显示出三座山的基本轮廓。由于计算机内存的限制，实验截 取了图像中坐标(6000,5000)到坐标(9000,11000)的3000×6000图像，该图 像刚好覆盖整幅图像中的三座山。图12为这两幅图像经过共轭相乘后的干涉图； 图13为该图像经过处理后得到的质量图，图14为最终的解缠相位图；通过对 比图14和图15，可以清楚看到图15中三座山的轮廓能够在图14中清晰可见。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的 技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本 技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的 本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明 创造均在保护之列。