一种基于质量图导引法和枝切法的干涉相位图解缠方法

摘要

本发明公开了一种基于质量图导引法和枝切法的干涉相位图解缠方法，首先采用枝切法对干涉相位图进行初步解缠，得到解缠相位图和解缠相位的二值分布图，再基于枝切法中的枝切线分布图和解缠相位的二值分布图找到未解缠区域与解缠区域边界线上的未解缠像素点，根据质量图找出与枝切线上未解缠像素点相邻的已解缠像素点中质量最好的点，将该点作为起始参考点用质量图导引法解缠与其相邻的未解缠区域中的像素，即用枝切线上的像素相位使两侧不连通的区域的相位连续起来，然后利用质量图导引法解缠被枝切线隔开的小部分非连通区域，大大降低了要解缠图像的大小，即可降低质量图导引法的计算量，进而减少了运算时间。

说明书

技术领域

本发明涉及，具体涉及一种基于质量图导引法和枝切法的干涉相位图解缠 方法。

背景技术

利用干涉技术进行数字高程提取在科研及测量领域有着广泛的应用，如合 成孔径雷达成像、表面微观形貌测量等。相位解缠是采用干涉技术进行数字高 程模型提取的重要步骤。

在干涉法测量高程过程中，由干涉图像提取的相位值反映了物体表面形貌 的高低变化，但通过相位提取算法得到的相位值均是用反正切函数表示的，相 位值被压缩在[-π,π]范围内，而物体表面形貌的相位变化范围通常远远超出一 个周期的相位变化范围。因此为了得到真正的高程信息，需要对缠绕相位进行 解缠处理，将阶跃分布的不连续的相位分布图拼接成连续的相位分布图。目前 存在的相位解缠算法主要有三大类：路径跟踪法、最小范数法和网络流法。

细胞大小和形态的检测对于生物研究及临床医学都有重大意义，细胞的高 度信息更是直接关系到细胞迁移过程中的变化情况，是研究其内部机理非常重 要的指标。目前多个国家都在对不同的细胞高度测量方法进行研究，相移干涉 法因其精度高、无损测量的优点，受到广泛关注。在利用相移干涉法测量细胞 高度的过程中，相位解缠对细胞高度的重构起着非常重要的作用。

理想情况下，图像的采样率满足Nyquist定律，采样频率大于信号最高频率 的两倍，即相邻相位差不可能超过半个周期（π）。记为真实相位，为 缠绕相位。根据Nyquist定律，真实相位应满足：

对相位φ(m)的解缠过程为：

k(m+1)为使上式成立的整数，

令$\mathrm{\Delta \omega }\left(m\right)=\mathrm{\Delta \phi }(m+1)-2\mathrm{\pi k}(m+1),$得到：

即理想情况下，通过对相邻缠绕相位之差积分就可以实现相位解缠。但当 相位图受到噪声及欠采样的影响时，就会使一些点不满足相邻缠绕相位差的绝 对值小于π的条件，这些点的解缠相位会发生错误，并且错误会向后传播，导致 之后所有相位的解缠结果与真实相位相差甚远。因此需要采用其他一些解缠算 法来解缠相位，避免误差传递对结果的影响。

常用的解缠算法——枝切法和质量图导引算法属于路径跟踪法，它们都分 别有明显的优点和不足。

Goldstein枝切法从格林公式中曲线积分与路径无关的条件出发，首先找出 缠绕相位数据上的不可微的奇点（残差点），然后用枝切线将正负残差点相连， 并使整个枝切线电荷平衡。最后用洪水法绕过枝切线解缠相位。这种算法运算 速度快，在噪声或欠采样影响不是很大时能够得到很好地解缠结果。但当图像 质量较差时，枝切线分布过密会将相位图分割为许多非连通区域，由于解缠时 不能跨过枝切线，因此这些孤立的区域得不到解缠。

质量图导引法依据生成的相位质量图按照像元质量从高到低的顺序解缠， 首先选取质量最高的像元为种子像元，从种子像元开始相位解缠，按照质量高 低以种子像元为中心向周围生长。这种算法抗噪声性能很强，又由于没有枝切 线的限制，所以不会存在得不到解缠的孤立区域。但这种算法的运算速度比较 慢，运算时间随图像尺寸的增加成指数增长，不利于解缠面积较大的图像。

发明内容

有鉴于此，本发明提供了一种基于质量图导引法和枝切法的干涉相位图解 缠方法，能够先用枝切法解缠大部分包裹相位，然后利用质量图导引法解缠被 枝切线隔开的小部分非连通区域，在保证解缠效果的前提下，降低运算量，节 省运算时间。

本发明的一种基于质量图导引法和枝切法的干涉相位图解缠方法，包括如 下步骤：

步骤一，计算干涉相位图f_O(x，y)的残差点，根据残差点设置枝切线，得到枝 切线分布图f_B(x，y)；

步骤二，选择枝切法解缠的开始点，采用枝切法解缠干涉相位图，得到解 缠的相位图f_u(x，y)以及该图对应的二值分布图f_b(x，y)

步骤三，采用质量图导引法求出干涉相位图的质量图f_q(x，y)

步骤四，找出所有既位于所述步骤二的枝切线分布图f_B(x，y)中枝切线上又位 于所述步骤二的二值分布图f_b(x，y)中未解缠区域内的未解缠像素点M_i；

步骤五，根据所述步骤二的二值分布图f_b(x，y)和步骤四中得到的未解缠像素 点M_i，找到与每一个未解缠像素点M_i相邻的已解缠像素点，将所有未解缠像素 点的相邻的已解缠像素点组成相邻像素点集合；所述相邻的已解缠像素点为该 未解缠像素点所在未解缠区域外，且与该未解缠像素点在水平或在竖直方向上 相邻的像素点；

步骤六，利用步骤三得到的质量图f_q(x，y)找出相邻像素点集合中质量最小的 点M_q，将该点作为参考点，该点的相位值为参考相位值；

步骤七，以参考点M_q为中心，利用质量图导引法解缠与该点相邻的未解缠 的区域中的像素点的相位，在解缠过程中更新解缠的相位图f_u(x，y)，得到更新后 解缠的相位图f_u′(x,y)，然后求解该更新后解缠的相位图f_u′(x,y)对应的二值分布图 f_b′(x,y)，用该二值分布图f_b′(x,y)更新原来的二值分布图f_b(x，y)，即： f_b(x,y)＝f_b′(x,y)；

步骤八，对步骤七的二值分布图f_b(x，y)的像素值进行求和运算，并判断和值 是否为0：如果和值不为0，执行步骤四至步骤八；如果和值为0，结束干涉相 位解缠过程，执行步骤九；

步骤九、输出步骤七的更新后解缠的相位图f_u′(x,y)。

所述步骤一具体包括如下步骤：

S201、扫描干涉相位图中的像元直到找到一个没有被枝切线连接的残差点， 将此残差点作为新枝切线的起点；然后再以该残差点为中心，在大小为3×3的 窗口内搜索另一个残差点；

S202、如果在所述窗口内搜索到另一个残差点，转到S203；如果搜索到干 涉相位图的图像边界，用枝切线连接残差点和图像边界将使枝切线极性置0，转 到S201；如果既没搜索到残差点，又没有搜索到图像边界，则执行S204；

S203、如果被搜索到的残差点未被枝切线连接，则连接该残差点并求连接 该残差点的枝切线的极性和，执行S205；否则执行S204；

S204、如果此枝切线上所有残差点都被作为搜索中心搜索过，将搜索窗口 的大小增加2，转至S202；否则将此枝切线上的其中一个未被作为搜索中心的 残差点设为搜索中心，转到S202；

S205、判断该枝切线的极性和是否为0：

如果为0，则标识整个枝切线极性平衡，然后再判断是否将整个干涉相位图 扫描完毕，如果是，结束枝切线设置过程；如果否，转到S201；

如果不为0，转至S204。

所述步骤三具体包括如下方法：

首先采用传统解缠算法分别沿x方向和y方向解缠相位图，得到相位图 f_x(x，y)和f_y(x，y)；

分别计算相位图f_x(x，y)和f_y(x，y)中像素(i，j)相位与相邻像素相位的平均值 和

$\overline{f_x}(i,j)=\frac{f_x(i,j)+f_x(i+1,j)+f_x(i-1,j)+f_x(i,j+1)+f_x(i,j-1)}{5};$

$\overline{f_y}(i,j)=\frac{f_y(i,j)+f_y(i+1,j)+f_y(i-1,j)+f_y(i,j+1)+f_y(i,j-1)}{5};$

然后计算各点相位与平均值的差：

${\Delta }_{i,j}^x=f_x(i,j)-\overline{f_x}(i,j),$${\Delta }_{i+1,j}^x=f_x(i+1,j)-\overline{f_x}(i,j),$${\Delta }_{i-1,j}^x=f_x(i-1,j)-\overline{f_x}(i,j),$${\Delta }_{i,j+1}^x=f_x(i,j+1)-\overline{f_x}(i,j),$${\Delta }_{i,j-1}^x=f_x(i,j-1)-\overline{f_x}(i,j)$和${\Delta }_{i,j}^y=f_y(i,j)-\overline{f_y}(i,j),$${\Delta }_{i+1,j}^y=f_y(i+1,j)-\overline{f_y}(i,j),$${{\Delta }_{i-1,j}^y=f}_y(i-1,j)-\overline{f_y}(i,j),$${\Delta }_{i,j+1}^y=f_y(i,j+1)-\overline{f_y}(i,j),$${\Delta }_{i,j-1}^y=f_y(i,j-1)-\overline{f_y}(i,j);$

最后，像素点(i，j)处的质量Q(i，j)为：

$\left(\begin{array}{c}Q(i,j)=\\ \sqrt{{{\Delta }_{i+1,j}^x}^2+{{\Delta }_{i-1,j}^x}^2+{{\Delta }_{i,j+1}^x}^2+{{\Delta }_{i,j-1}^x}^2+{{\Delta }_{i,j}^x}^2}+\sqrt{{{\Delta }_{i+1,j}^y}^2+{{\Delta }_{i-1,j}^y}^2+{{\Delta }_{i,j+1}^y}^2+{{\Delta }_{i,j-1}^y}^2+{{\Delta }_{i,j}^y}^2}\end{array}\right)$

。

所述步骤七具体包括如下步骤：

S701、在参考点M_q相邻的未解缠区域内，以M_q作为相位展开的参考点，把 该点的所有相邻的未解缠像素点分别进行相位展开，将相位展开像素点标记为 已展开点；同时把已展开点各自相邻的未解缠点按质量值的大小排序存放到一 个毗邻队列Q中；

S702、取出毗邻队列Q中质量值最小的点，以该点的其中一个相邻的相位 展开像素点为参考点进行相位展开，并标记该质量值最小点为已展开点，把该 质量值最小点的相邻的未解缠点放入毗邻队列Q中，并按质量值大小对队列Q 重新排序；

S703、执行S701，直到毗邻队列Q为空。

本发明的一种基于质量图导引法和枝切法的干涉相位图解缠方法，具有如 下有益效果：

本发明首先采用枝切法对干涉相位图进行初步解缠，得到解缠相位图和解 缠相位的二值分布图，再基于枝切法中的枝切线分布图和解缠相位的二值分布 图找到未解缠区域与解缠区域边界线上的未解缠像素点，根据质量图找出与枝 切线上未解缠像素点相邻的已解缠像素点中质量最好的点，将该点作为起始参 考点用质量图导引法解缠与其相邻的未解缠区域中的像素，即用枝切线上的像 素相位使两侧不连通的区域的相位连续起来。本发明的方法继承了质量图导引 法在噪声比较严重的区域保持的良好鲁棒性的优点，同时由于先使用枝切法解 缠大部分像元，保证了相位解缠的效果，然后利用质量图导引法解缠被枝切线 隔开的小部分非连通区域，大大降低了要解缠图像的大小，即可降低质量图导 引法的计算量，进而减少了运算时间。

附图说明

图1为四相邻相位示意图。

图2为加入椒盐噪声的模拟细胞干涉相位图。

图3为加入椒盐噪声的模拟细胞干涉相位图（图2）的残差点分布图，图中 白点为正残差点，黑点为负残差点。

图4为枝切线分布图，根据残差点分布图（图3）得到。

图5为由枝切法得到的图2的解缠相位立体图。

图6为由枝切法得到的解缠相位灰度图。

图7为由枝切法得到的解缠相位的二值分布图，白色区域为已解缠区域， 黑色区域为未解缠区域。

图8为图2的质量图。

图9为由质量图导引法得到的解缠相位灰度图。

图10为由质量图导引法得到的解缠相位立体图。

图11为本发明的方法得到的解缠相位立体图。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

本发明提供了一种基于质量图导引法和枝切法的干涉相位图解缠方法，包 括如下步骤：

步骤一，如图3所示，计算干涉相位图f_O(x，y)的残差点，根据残差点设置枝 切线，得到如图4所示的枝切线分布图f_B(x，y)；

Goldstein枝切法中残差点的定义是：如图1所示，将四相邻像素的相位差 顺时针加和，加和的结果成为残差，如图3所示，残差为正的点为正残差点， 残差为负的点为负残差点。

相位差的定义为：将两相邻像素相位相减，若差小于-π正或大于则通过加减 n倍的2π使其在[-π，π]之间即得两像素相位差，n为整数。

Δp((i，j)，(i，j+1))=p(i，j)-p(i，j+1)+n*2π，n=…-3，-2，-1，0，1，2，3，...

Δp((i，j)，(i，j+1))为点(i，j)和(i，j+1)的相位差，p(i，j)和p(i，j+1)分别为点 (i，j)和(i，j+1)的相位，Δp((i，j)，(i，j+1)满足-π＜Δp((i，j)，(i，j+1))＜π。

其中，枝切线的设置步骤如下：

S201、扫描干涉相位图中的像元直到找到一个没有被枝切线连接的残差点， 将给残差点作为新枝切线的起点；然后再以该残差点为中心，在大小为3×3的 窗口内搜索另一个残差点；

S202、如果在所述窗口内搜索到另一个残差点，转到S203；如果搜索到干 涉相位图的图像边界，用枝切线连接残差点和图像边界将使枝切线极性置0，转 到S201；如果既没搜索到残差点，又没有搜索到图像边界，则执行S204；

S203、如果被搜索到的残差点未被枝切线连接，则连接该残差点并求连接 该残差点的枝切线的极性和，执行S205；否则执行S204；

S204、如果此枝切线上所有残差点都被作为搜索中心搜索过，将搜索窗口 的大小增加2，转至S202；否则将此枝切线上的其中一个未被作为搜索中心的 残差点设为搜索中心，转到S202；

S205、判断该枝切线的极性和是否为0：

如果为0，则标识整个枝切线极性平衡，然后再判断是否将整个干涉相位图 扫描完毕，如果是，结束枝切线设置过程；如果否，转到S201；

如果不为0，转至S204。

步骤二，选择枝切法解缠的开始点，采用枝切法解缠干涉相位图，得到解 缠的相位图f_u(x，y)以及该图对应的二值分布图f_b(x，y);其中的二值图分布图由以 下方法得到：在解缠的相位图f_u(x，y)中，将已解缠像素点的值置为1，将未解缠 像素点置为0，则形成二值分布图f_b(x，y)。

步骤三，采用质量图导引法求出干涉相位图的质量图f_q(x，y)，具体方法如下：

首先采用传统解缠算法分别沿x方向和y方向解缠相位图，得到相位图 f_x(x，y)和f_y(x，y)；

分别计算相位图f_x(x，y)和f_y(x，y)中像素(i，j)相位与相邻像素相位的平均值 和

$\overline{f_x}(i,j)=\frac{f_x(i,j)+f_x(i+1,j)+f_x(i-1,j)+f_x(i,j+1)+f_x(i,j-1)}{5};$

$\overline{f_y}(i,j)=\frac{f_y(i,j)+f_y(i+1,j)+f_y(i-1,j)+f_y(i,j+1)+f_y(i,j-1)}{5};$

然后计算各点相位与平均值的差：

${\Delta }_{i,j}^x=f_x(i,j)-\overline{f_x}(i,j),$${\Delta }_{i+1,j}^x=f_x(i+1,j)-\overline{f_x}(i,j),$${\Delta }_{i-1,j}^x=f_x(i-1,j)-\overline{f_x}(i,j),$${\Delta }_{i,j+1}^x=f_x(i,j+1)-\overline{f_x}(i,j),$${\Delta }_{i,j-1}^x=f_x(i,j-1)-\overline{f_x}(i,j)$和${\Delta }_{i,j}^y=f_x(i,j)-\overline{f_y}(i,j),$${\Delta }_{i+1,j}^y=f_y(i+1,j)-\overline{f_y}(i,j),$${\Delta }_{i-1,j}^y=f_y(i-1,j)-\overline{f_y}(i,j),$${\Delta }_{i,j+1}^y=f_y(i,j+1)-\overline{f_y}(i,j),$${\Delta }_{i,j-1}^y=f_y(i,j-1)-\overline{f_y}(i,j);$

最后，像素点(i,j)处的质量Q(i,j)为：

$\left(\begin{array}{c}Q(i,j)=\\ \sqrt{{{\Delta }_{i+1,j}^x}^2+{{\Delta }_{i-1,j}^x}^2+{{\Delta }_{i,j+1}^x}^2+{{\Delta }_{i,j-1}^x}^2+{{\Delta }_{i,j}^x}^2}+\sqrt{{{\Delta }_{i+1,j}^y}^2+{{\Delta }_{i-1,j}^y}^2+{{\Delta }_{i,j+1}^y}^2+{{\Delta }_{i,j-1}^y}^2+{{\Delta }_{i,j}^y}^2}\end{array}\right)$

，Q值越小对应相位图中的像素质量越好。

采用上述枝切法对干涉相位图进行解缠后，仍存在一些未解缠的区域，因 此本发明考虑采用抗噪声性能更好的质量图导引法进一步对未解缠区域进行处 理。为了使已解缠区域的相位能够与之后要采用的质量图导引法来解缠区域的 相位连接起来，所以首先要从未解缠区域边界线上的像素点选取参考点作为质 量图导引法的解缠起始点。

步骤四，找出所有既位于枝切线分布图f_B(x,y)中枝切线上又位于解缠相位的 二值分布图f_b(x,y)中未解缠区域内的未解缠像素点Mi_;

步骤五，根据解缠相位的二值分布图f_b(x,y)和步骤三中得到的未解缠像素点 M_i,找到与每一个未解缠像素点M_i相邻的已解缠像素点；所述相邻的已解缠像 素点为该未解缠像素点所在未解缠区域外，且与该未解缠像素点在水平或在竖 直方向上相邻的像素点；所有未解缠像素点的相邻的已解缠像素点组成相邻像 素点集合；

步骤六，利用质量图f_q(x,y)找出相邻像素点集合中质量最小的点M_q,将该点 作为参考点，该点的相位值为参考相位值；

步骤七，以参考点M_q为中心，利用质量图导引法解缠与该点相邻的未解缠 的区域中的像素点的相位，具体步骤为：

S701、在参考点M_q相邻的未解缠区域内，以M_q作为相位展开的参考点，把 该点的所有相邻的未解缠像素点分别进行相位展开，将相位展开像素点标记为 已展开点；同时把已展开点各自相邻的未解缠点按质量值的大小排序存放到一 个毗邻队列Q中；

S702、取出毗邻队列Q中质量值最小的点，以该点的其中一个相邻的相位 展开像素点为参考点进行相位展开，并标记该质量值最小点为已展开点，把该 质量值最小点的相邻的未解缠点放入毗邻队列Q中，并按质量值大小对队列Q 重新排序；

S703、执行S701，直到毗邻队列Q为空。

在解缠过程中同时更新解缠相位的二值分布图f_b(x,y)和解缠的相位图 f_u(x,y);

按照上述方法，完成参考点M_q相邻未解缠区域的相位解缠，为对其他区域 继续解缠，执行步骤八；

步骤八，对解缠相位的二值分布图f_b(x,y)的像素值进行求和运算： S＝Σ-f_b(x,y)，并判断和值S是否为0：如果和值S不为0，执行步骤四；如果 和值S为0，结束干涉相位解缠过程，执行步骤九；

步骤九、输出解缠的相位图f_u(x,y)。

从图5中可以看出，用枝切法得到解缠相位图中有若干空洞，这些空洞对 应于图6中的白色区域和图7中的黑色区域，是枝切法没有解缠的区域。从图9 和图10可以看出，由质量图导引法得到的解缠相位图中不存在大片的未解缠区 域，说明质量图导引法的鲁棒性更好。从图11中可以看出，通过本发明得到的 解缠相位图没有图5、6、7中的大片未解缠区域，本发明的解缠效果接近质量 图导引法解缠效果，并且用时更短。

表1对使用枝切法、质量图导引法和本算法解缠包裹相位（像素200×200） 的结果做了对比，t表示运算时间，η为像元解缠率（已解缠像素数与总像素数 的比）。本算法的解缠率比枝切法的解缠率提高了13%，而运算时间只有质量图 导引法的17%。从解缠效果（图5、10、11）上看，三种算法都能够较好地还原 原细胞的相位，虽然由于噪声的存在使一些点处没有得到平滑的解缠结果，但 并没有出现误差传递。

表1三种算法性能对比（图像尺寸：200×200）

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保 护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等， 均应包含在本发明的保护范围之内。