自适应分块的改进最小费用网络流解缠算法

张金翼，王 辉，吴思利，郑世超，顾约翰

（1.上海卫星工程研究所，上海 201109；2.上海市毫米波空天信息获取及应用技术重点实验室，上海 201109）

0 引言

近年来，地下资源的过度开采、水土流失等问题愈发严重，地质灾害（如地震、泥石流、矿区塌陷、城市地面沉降）频发，对我国人民群众的生命及财产安全造成威胁。在该类灾害中典型的表征为地表形变，因此监测地表形变是防灾减灾的重要手段。传统的测量地表形变技术手段为水准测量、GPS 测量等［1］，但上述方法存在费时、费力、数据点孤立的缺陷，无法及时、高效地预测和评估地质灾害。差分干涉合成孔径雷达（Differential Interferometric Synthetic Aperture Radar，DInSAR）是天基遥感中较为常用的测量地表形变的技术手段，因为全天时、全天候工作，速度快、范围广、测量密度高等特点逐渐受到人们的重视［2-3］，差分干涉合成孔径雷达形变测量主要利用相位的变化，直接测得的干涉相位被包裹在［-Π，Π）)[-π，π 内，要获得可用的相位信号须对其进行解缠绕，其效果对测量的精度有较大的影响［4-5］。

常用的解缠方法按照原理分为以下3 类：1）基于路径跟踪解缠，选择不同的路径对相邻像元进行差值积分恢复真实相位，较有代表性的是GOLSTEIN［6］于1988 年提出的枝切法，可在计算出残差点后，快速有效地寻找连接相邻残差点的最优路径；2）基于最小范数思想，通过拟合函数求得缠绕相位和解缠相位的离散偏微分差，进而求解相位解缠的整体最佳估计值，其中较为常用的是最小二乘法、Jacobi 迭代法、Gauss-Seidel 迭代法和逐次超松弛（Successive Over Relaxation，SOR）迭代法，这些算法的缺陷在于数据量较大时收敛速度慢，不适用于大规模数据［7-8］；3）基于网络规划的算法，由COSTANTINI［9］于1996 年提出，网络规划算法 在相位解缠中引入最小代价流的概念，通过搜索全区域最短路径生成最优枝切线进行解缠，可以将误差限制在相干质量差的区域内，保证其余区域结果的准确性，精度较高，因此常被用作其他算法的精度参照，但由于最小费用网络流（Minimum Cost Network Flow，MCF）的算法特性，当相位差图像相干性差、残差点数量多时，算法的复杂度较高，效率较低［10-12］。CHEN［13］提出改进的MCF 算法，分割整个图像，与传统的MCF 算法相比，提高了算法效率，但分块的规模为人工设定，没有具体的标准，不同的分块大小会对该算法的精度产生较大影响。

本文提出一种基于自适应分块的改进MCF 解缠算法，该方法根据质量图所表现的残差信息，将分块大小进行自适应寻优，更多地将高相干区域划分在一个子块。实验结果表明，该算法在提升算法效率的同时保证了准确度，有效地拓展了MCF 算法的应用场景。

1 MCF 算法解缠原理

相位解缠的基本思路是将已获取的干涉图中的缠绕相位，通过叠加缠绕数恢复相位的真实值，进而反映真实的形变情况［14-15］。

定义缠绕相位为ψi，j，真实相位为ϕi，j，对于一个M×N的矩阵，可得方程式如下：

式中：k为整数，0 ≤i≤M-1，0 ≤J≤N-1。

定义缠绕算子W，对相邻像素间的缠绕相位进行±2kπ 的操作，使其变为真实相位。定义x方向和y方向上的绝对相位差分别为，其方程式如下：

式中：ψi+1，j为M×N矩阵中(i+1，j)处点的缠绕相位，rad；ψi，j+1为M×N矩阵中(i，j+1)处点的缠绕相位，rad。

由式（2）变换可得方程式如下：

对于干涉图中大部分点，沿着点(i，j)、(i+1，j)、(i+1，j+1)、(i，j+1)的缠绕相位差环路积分为0，而实际应用中，由于噪声等因素的存在导致出现残差点，残差点处的缠绕相位差环路积分为±2π，即定义残差为

式 中：ei，j为残差；为(i+1，j)点在y方向上的绝对相位差，单位为rad；为(i，j+1)点 在x方向上的绝对相位差，单位为rad。

由于存在残差点，无法确定干涉图全图的积分，其变成与积分路径有关的变量，因此需考虑残差点的影响，使其降至最小。在最小费用流解缠算法中，将正/负残差点定义为供应/需求节点，残差为供应/需求量，将其他点定义为转运节点。相邻的节点之间通过带有流的弧连接，以调节供需平衡，最小费用流解缠的目的是在所有节点供需平衡的基础上，确保总费用取最小。达到供需平衡即［16-18］：

式 中：为 点(i，j)到点(i+1，j)的流量；为点(i，j)到点(i，j+1)的流量，1 ≤i≤M-1，1 ≤j≤N-1。

总费用流为

式中：为点(i，j)到点(i+1，j)对应的费用；为点(i，j)到点(i，j+1)的费用。

式中：ψ0，0为（0，0）位置处点的缠绕相位，rad。

综上所述，应用MCF 算法的解缠步骤如下。

1）对相干图做预处理，选定阈值，提取相干系数高于阈值的相位。

2）在相位的集合中建立Delaunay 三角网，进而构建对偶图，将残差映射到对偶网络中。

3）在对偶网络中，应用最小费用流法连接正负残差点对，计算最小费用流集合。

4）根据流的大小和方向对相位矩阵积分，得到解缠结果，再从高质量区域向低质量区域积分［12］。

2 自适应分块的改进MCF 解缠算法

MCF 算法一经推出，其准确度便得到肯定，但随着残差点的增多，该算法的效率有所下降，因此提高算法效率是优化该算法的一个方向［16］。CHEN［13］提出改进的MCF 算法，将整个相干图划分为不重叠的矩形块，将矩形块作为完全独立的图像进行相位解缠，降低解缠过程中所需的内存资源，提高效率。

原始相干图使用MCF 算法解缠的算法复杂度可表示为O(N2)，其中N为图中残点的数量，将相干图分为i个矩形块，则分割后算法复杂度可表示为O(N12+N22+N32+…+Ni2)，其 中N1+N2+N3+…+Ni=N，后者的算法复杂度远小于前者。

该种方法虽然提升了相位解缠的效率，但由于分割大小的不同会导致解缠精度的不稳定。在本文中引入一种基于质量图的自适应分块改进MCF算法，切割时，根据质量图自适应选取分割大小，保留相干性较强的区域，在提高效率的同时保证精度。

本文提出的算法流程如下。

1）对相干图做预处理，选定阈值，提取相干系数高于阈值的相位。

2）根据相干系数的计算方式计算质量图。

在干涉SAR 中残差点的存在是由于SAR 叠掩或去相关噪声引起的，因此相关系数较低的区域其残差点较多，相关系数较高的区域残差点较少。一般将表征干涉图像各点相位的特征图称为质量图，质量图可划分为相干系数图、伪相干系数图、相位导数变化图及最大相位梯度图4 种，本文采用使用范围较广的相干系数图［20］。

相干系数图中各个区域之间的相干性通过相干系数表示，相干系数值较高表示2 个区域之间的相关性较好，即相位信息准确；相干系数值较低说明2 个区域之间相关性较差，即相位信息有误。A、B2 张复图像的相干系数可表示为

式中：E为求期望；γ为相干系数；B*为B的共轭。

3）根据质量图选定应进行分割的块的尺寸。

质量图表示相干系数的大小，在最小费用流解缠方法中，影响解缠效率与准确度的是干涉图中集中存在的残差点。因此在分割时，将相干性较强的区域集中在一起，进而使得残差点被零散的分割到各个子块内，提高每个分割块的解缠效率。

在实际操作过程中，为滤去噪声求整张质量图的均值，将其作为门限滤波，表达式为

式中：T为门限值；M、N为质量图规模；f(i，j)为(i，j)位置点的相干系数。

在选择分块尺寸的过程中，选用密度峰值聚类算法。该算法是目前较为先进的聚类算法，需要较少的输入参数，能检测聚类中心，具有简单、聚类速度快、稳定性好、聚类高效的优点，且对噪声不敏感，是较为理想的聚类算法［21-22］。

对于质量图中的数据，采用欧氏距离计算两点之间的相似性：

式中：(xi，xj)为第i个点和第j个点；l为维度；m为维度的数量。

计算数据集中每个数据点的密度，表达式如下：

式中：S(m)为判断函数；dc为截断距离，设定为5。

当m≥0 时S(m)=1；当m≤0 时S(m)=0。

依照式（11），即可计算得到以每一点为中心的与其距离小于5 的点的个数，再求取每一数据点与其他密度比其大的数据点间的最小距离，如该点密度最大，则取与其他数据点距离的最大值为

式中：di，j为第i个点与第j个点的距离。

求得密度及最小距离后，将横轴表示数据点的密度，纵轴表示数据点的最小距离，选取密度高且距离大的点作为聚类中心点，再将其他点按密度值降序排列，归于聚类中心点［23-25］。

得到聚类结果后，将聚类点所在的矩形区域作为参考的干涉图分块区域，如聚类后产生多个子块，将其大小取均值作为最终的分块参考大小。

依照此方法得到的分块可将相干性较强的点集中在同一块区域，并依照大小划分其他区域，提高算法的效率和准确度。

4）在相位的集合中建立Delaunay 三角网，进而构建对偶图，将残差映射到对偶网络中。

5）在对偶网络中应用最小费用流法连接正负残差点对，计算最小费用流集合。

6）根据流的大小和方向对相位矩阵积分，得到解缠结果，再从高质量区域向低质量区域积分［12］。

算法流程如图1 所示。

图1 自适应分块的改进MCF 解缠算法流程Fig.1 Flow chart of the improved MCF unwrapping algorithm with adaptive chunking

3 仿真实验及结果分析计算机模拟结果

仿真实验在一台配备Intel Core i7-10750H 处理器（2.6 GHz，6 核12 线程）、NVIDIA GeForce GTX 1660 Ti GPU（6 GB 显存）和16 GB DDR4 内存的计算机上进行，运行于Windows 10 操作系统，使用Matlab R2019b 仿真，在该平台生成干涉条纹图，并按照本文方法对其解缠绕，解缠前干涉如图2所示。

图2 解缠前干涉Fig.2 Interference image before phase unwrapping

3.1 仿真实验

上文中所提到的方法对干涉图像做处理，经自适应算法选取可知高相干点位于（80～170，50～180），块大小为90×130，将上图保留高相干块并依照其大小划分其他区域，进而对其解缠绕，与从小到大随机选取其他分块大小做解缠对比，解缠后实验结果如图3 所示。图3（a）～图3（f）分别为寻优后分块的解缠结果、选取分割块大小为50×50 的解缠结果、选取分割块大小为100×100 的解缠结果、选取分割块大小为150×150 的解缠结果、选取分割块大小为200×200 的解缠结果和选取分割块大小为256×256 的解缠结果。图3（a）～图3（f）的解缠全流程处理时间分别为4.07、2.18、3.32、4.55、4.51 和8.92 s。

图3 不同分块方式解缠后图像Fig.3 Images after unwrapping with different chunking methods

不同分块解缠结果误差如图4 所示。其中图4（a）～图4（f）的解缠后均方根误差（Root Mean Square Error，RMSE）分别为4.002 8、5.564 1、5.574 6、5.232 5、5.232 5 和3.962 9 rad。

图4 不同分块方式解缠图像后均方根误差Fig.4 RMSE results after unwrapping the images with different chunking methods

仿真试验结果统计见表1。

表1 仿真试验结果Tab.1 Simulation results

3.2 结果分析

由上述仿真结果可知，采用文中算法选取到的分割大小可在4.07 s 内完成解缠，误差为4.002 8 rad，对比其他随机选取的分块大小，本算法在保证解缠效率的同时，其准确度没有大幅下滑。对比误差最小的分块方式（256×256），解缠误差恶化了1.01%，但效率提升了54.37%。证明本算法在平衡效率和准确度方面具有一定的优势，对提升改进MCF 算法的应用场景与计算效率有一定的贡献。

4 结束语

针对改进的基于分割的最小费用网络流算法准确率及效率受分割块大小影响的问题，提出利用质量图表征相干性优劣的特性自适应确定分割块大小的方法，进而将相干性较好的区域聚集在一起，将残差点散落地分布在不同分割块中。仿真结果表明，本方法通过自适应寻优方式选取合适的分割块，使得准确度在不受过多影响的情况下，有效地提高算法效率。不足之处是在分割干涉图进行的过程中，存在将其分块过小的可能，进而导致解缠后图像出现较为严重的马赛克现象。在后续研究中，将探究如何对自适应分块后的图像做针对性的解缠、融合处理。