融合梯度向量流条纹探测与图切法相位解缠

刘成洲, 李 斌, 许艺腾

(中交天津港湾工程研究院有限公司岩土工程研究所， 天津 300202)

相位解缠是一个将缠绕相位恢复到绝对相位的过程，是合成孔径雷达干涉测量(interferometric synthetic aperture radar, InSAR)处理过程中至关重要的一环。解缠结果的正确与否直接影响InSAR在矿区变形监测的可靠性和精度[1]。大多数的相位解缠算法是建立在理想的缠绕条件下，即邻近像元的相位差分绝对值小于π，而实际的相位表面往往是不连续，或是缠绕相位含有噪声，这类算法很难正确解缠出真实相位。目前，相位解缠算法大致可以分为四类：路径跟踪[2]、最小Lp范数[3]、贝叶斯正则化[4]和参数模型[5]。路径跟踪法，是通过设置合理的积分路径，以积分相邻像元上的缠绕相位梯度的方式进行相位解缠。比较常见的路径跟踪法有Goldstein枝切法[6]和质量图引导的枝切法[7]。最小Lp范数法，是一种满足全局最优条件的相位解缠算法，该算法使解缠图像的误差达到相对最小，其解缠的结果具有连续性，不会出现孤岛现象。比较常见最小Lp范数法有最小二乘算法[8]、图切法[9]。贝叶斯正则化是依靠数据观察机制，以及相位先验知识建立的模型，常见的贝叶斯正则化方法有随机非线性相位重构法、基于局部平滑的自适应正则化相位估计。参数模型算法将解缠相位约束到一个参数表面。在文献[10]中，对图像进行分块，每一个块应用不同的参数模型进行解缠，最后将拟合后的图像融合成解缠相位。此外，还有遗传算法[11](利用遗传算法优化相位解缠参数，使得缠绕和解缠绕相位梯度差别位置总数最小)、蚁群算法[12](通过蚁群算法求取“枝切线”上残差点的最短路径，使得枝切线不易形成封闭路径，避免了部分区域无法解缠)、卡尔曼滤波[13](其将相位解缠问题转化为状态估计问题，建立相位的动态方程和观测方程来求解真实相位，在解缠的同时也消除了相位噪声)、区域生长法[14](利用区域增长法对缠绕相位进行编码标记，通过叠加的方式求取解缠相位)等应用到相位解缠。解缠方法很多，很难说哪种方法具有绝对的优势，因此，需要根据不同的数据特点选取不同的解缠方法。

针对矿区在沉降变形是出现的漏斗状变形，同时矿区地表由于覆盖农田、水域和树木等散射性强的地物，使其在干涉相位上表现为同心环状，不连续，高噪声的特点，利用本文提出的融合GVF-Snake条纹探测与马尔科夫随机场(Markov random field,MRF)图切法能够实现矿区变形精细化解缠。

1 算法的基本原理

1.1 基于GVF-Snake模型边界探测

Xu等[15]针对Snake模型无法收敛至凹型边界以及对初始化轮廓较为敏感等问题，提出了GVF-Snake模型。该模型重新定义了外部能量场，即外部能量场用全局梯度向量场来表示。具体实现流程：首先利用Sobel算子提取边界图f(x,y),然后通过极小化能量方程，得到整个图像域的梯度矢量流场V(x,y)=[u(x,y),v(x,y)]。则能量E方程可表示为

(1)

GVF场即式(1)中V(x,y)达到最小时E的值，并用它来代替Snake中Eext[Φ(s)]。为求能量函数的极小值，通常采用变分法，式(1)的解满足

(2)

(3)

引入时间变量t,把u、v看成时间t的函数，采用有限差分离散化迭代求解，解得u、v为

(4)

(5)

1.2 MRF图切法基本原理

Bioucas-Dias等[9]将缠绕相位当做目标函数，并用一阶马尔科夫随机场来表示。通过求取目标函数的最大后验估计来估计边界，并用图论中的最大流最小割对图像进行分割达到解缠的目的。

定义能量函数为

(6)

式(6)中：i、j表示像元的行列数，(i,j)∈G0≡{(k,l):k=1,2，…,M,l=1,2，…,N}(G0是图像二维格网行列索引号)；V(·)是由多个相互毗邻集合定义出的集合函数；hij表示沿着水平不连续参数一阶邻近，vij是沿着竖直不连续参数一阶邻近，hij、vij∈{0,1}，当hij、vij=0时，信号是不连续的，像元不连续。

式(6)各分量计算式为

(7)

(8)

(9)

(10)

式中：k≡{kij∈Z:(i,j)∈G0}是标记场；ψ≡{ψij∈[-π,π):(i,j)∈G0}是观测得到的缠绕相位； 上标h、v表示水平向和竖直向。

1.3 融合条纹探测与MRF图切法相位解缠

算法的实现包括以下步骤：

(1)利用基于GVF-Snake模型边界探测的相位解缠算法进行分步解缠，假设解缠到第n步。

(2)判断n步模型检测边界是否为真。

(3)对“不为真”的情况，边界探测的相位解缠已经实现了n-1步解缠，定义n-1步探测条纹外侧为解缠外部块，内侧为缠绕内部块。由于第n步条纹探测失败，需要对缠绕内部块进行MRF图切法相位解缠，将其解缠结果按照解缠准则增加或减少2π(n-1)后，再与缠绕外部块进行融合(具体的解缠准则为由外到内，相位从π向-π变化，为加，否则为减)，得到含有边界孤立点的粗解缠图像。

(4)高通滤波插值法消除边界上的孤立点，获得最终的解缠图像。

其流程如图1所示。

图1 融合条纹探测与MRF图切法相位解缠Fig.1 Fusion fringe detection and phase unwrapping by MRF graph cuts

流程中通过距离判别条纹是否为真，首先目视法手动选取初始边界点(与真实边界相差很小，可看作真实边界)，寻找每一个初始边界点到与之对应的GVF-Snake模型探测边界的最小距离的点，即

(11)

式(11)中：Xm、Ym为GVF-Snake模型探测边界的任意点；xm、ym为第m个初始化边界点坐标；m是小于等于初始化边界点个数的任意正整数，m=1,2,3,…。

通过多次试验获得的最小距离阈值，判定依据是通过真实条纹与探测边界的对比，探测条纹真假很容易分辨。详情如式(12)和图2所示。

(12)

式(12)中：R为判别率，规定R=0时，探测边界线“为真”，否则，“为假”；τ为最小距离阈值，如图2所示。

图2 判别探测边界是否为真Fig.2 Discriminate whether the detection boundary is true or not

图2是判断GVF-Snake探测边界是否为真的示意图。图2中红线表示其实际探测曲线，黑色线表示通过多次试验获得最小距离阈值曲线，蓝色点表示手动选取得初始化边界点，黑色箭头表示初始化边界点到探测边界得最小距离，当黑色箭头线在两条阈值边界之间，表示条纹探测线与真实边界相符，否则条纹探测失败。

1.4 高通插值滤波

图像粗解缠后会在其边界形成很多孤立点，这些孤立点相对高频，因此可以使用高频滤波将这些边界点标记出来，进行去除，然后再用双三次内插法将去除的点进行插值。实现步骤如图3所示。

图3 高通滤波插值法流程图Fig.3 Flow chart of interpolation method for high-pass filtering

从图3可知，该滤波算法的核心是将图像转换到频率域，通过设置高通滤波器来检测高频点，然后通过傅里叶逆变换转换到图像域。傅里叶变换与反变换，可以表示为

(13)

(14)

式中：M、N分别是图像的行数和列数；u=x=1,2,…,M,v=y=1,2,…,N，傅里叶变换后的相位谱为

(15)

将高频变换滤波器设置为

(16)

式(16)中：D0为滤波器的阻带半径；D(u,v)是点到滤波器中央的距离。

最后，将检测出的高频点进行剔除，并用双三次内插法对空值进行插值达到滤波目的。该算法的缺点是对不是边界的高频值也具有滤波作用，滤波后整幅图会更加平滑,滤波效果如图4所示。

图4 高通插值滤波效果图Fig.4 High-pass interpolation filtering effect map

图4(b)反映了高通滤波剔除高频噪声后的图像，白色表示空值，为高频噪声点，通过对图4(b)进行双三次内插法插值填补空区域，得到连续的解缠图像。

2 实验结果与分析

2.1 研究区域概况

研究区域为菏泽市巨野矿区郭屯煤矿某工作面，工作面走向全长(平距)1 795 m，面长(平距)125 m，开采标高-727～-815 m，表土层厚度约572.1 m，地面标高+42.9～+45.4 m。研究区域及布线情况如图5所示。

图5 研究区域与布线情况Fig.5 Research area and location

2.2 数据选取

选择sentinel-1A Interferometric Wide (IW) 模式，这种模式下SAR卫星的空间分辨率平均13 m(空间分辨率5～20 m),SAR卫星的波长为5.546 cm，选取的两组数据的基本信息如表1所示。

表1 选取的SAR数据的基本信息Table 1 Basic information of selected SAR data

2.3 数据处理预处理

实验预处理数据通过SNAP ESA软件进行干涉生成、去平地效应、去地形相位，然后经过地理编码得到，如图6所示。

图6 干涉数据获取流程图Fig.6 Flow chart of interference data acquisition

2.4 解缠算法分析

GVF-Snake模型解缠过程如图7(a)～图7(c)及其解缠结果如图7(d)与本文的解缠结果如图7(e)所示。可以看出，由于干涉条纹内部不连续性情况严重，致使探测边界“不为真”，使得解缠失败。通过本文方法的改进后，能够正确进行解缠，得到粗解缠图像。这些粗解缠图像条纹边界上分布少量孤立的散点,可以通过高通滤波插值法进行消除，得到最终的解缠图像。

图7 GVF-Snake模型解缠过程及融合前后解缠结果Fig.7 The unwrapping process of GVF-Snake model and the result of unwrapping before and after fusion

Bioucas-Dias等[16]将PEARLS算法用于相位重构，利用重构后绝对相位与真实地形进行比对，其误差众数集中在零上，且误差上限为+0.6 mm,误差下限-0.8 mm,在±0.3 mm之间的误差占比为0.936，这说明该方法对地形变形重构是很有效的，可以作为检核其他相位解缠算法的依据。在面分析上，本文列出了原始的相位图像[图8(a)]，选取最小费用流、高金斯枝切法等5种经典的相位解缠方法[图8(b)～图8(f)]与本文方法[图8(h)]进行对比分析，其中以图8(g)PEARLS相位重构为参照相位。

图8 多种相位解缠图像进行对比Fig.8 Comparisons of various phase unwrapping images

从图8中可看出各种相位解缠算法都能够解缠出相位的变化趋势，但部分算法缺点表现得很明显。其中最小费用流容易产生局部最优解，致使条纹不连续情况严重时，局部地区无法正确解缠；Lp最小范数法，对于整体解缠最优，忽略了细节的表达;高金斯枝切法由于受到间断相位的干扰，在图像下部出现解缠错误。解缠效果最好的三种算法，分别为本文方法、最小不连续法和Snaphu法相位解缠。

为定量评价本文相位解缠的精度，使用了5种评价指标来综合评定解缠相位的精度。

(1)平均绝对误差(average absolute error，MAE)。

(17)

(2)均方误差(mean square error,MSE)和均方根误差(root mean square error,RMSE)。

(18)

(19)

(3)绝对误差中值(median absolute error different,MED)。

(20)

(4)相位的最大绝对误差(max absolute error,MAX)。

(21)

式中：M是相位的行数；N是相位的列数；A为标准相位；B为待评价相位，median为取中值函数；max为取最大值函数。以PEARLS相位重构作为参照相位，将6种解缠方法的解缠图像作为待评价相位，计算出面评价精度指标，如表2所示。

由表2可知，最小不连续法是仅次于本文方法的解缠方法，从各项指标上看均优于其他四种；Lp最小范数法虽然直观上看上去解缠效果不佳，但其面分析指标上要比高金斯枝切法和最小费用流法，这也说明了其整体参数最优的解缠特性；其次比较好

表2 面分析评价指标Table 2 Evaluation index of surface analysis

的是Snaphu相位解缠方法，各项指标均优于最小费用流法和高金斯直切法；最好的解缠方法是本文方法，除平方根误差未达到最佳，其他指标均优于其他五种；最差劲的是高金斯枝切法。按照面性精度排名：本文方法>最小不连续法>Snaphu>Lp最小范数>最小费用流>高金斯枝切法。

也对解缠算法进行了线分析。主要步骤如下：

(1)将解缠后的绝对相位变形方向与视线向(LOS方向)转换到竖直变形，并去除大气相位，使绝对相位只含有地面变形信息。

(2)将相位图像上提取工作面以及布点信息。

(3)与实际水准测量数据进行对比。

绝对相位由视线向变形转换到竖直向变形，即

(22)

式(22)中：v为竖直形变；d为视线向形变；λ为哨兵1Å的波长；θ为入射角。

去除大气相位，选取Yu等[17-18]提出的迭代对流层分解模型(iterative tropospheric decomposition, ITD),此法是对海拔与涌动对流层延迟部分进行耦

合所得到的模型。将ITD修正相位转化到绝对相位下，然后通过解缠后的图像与之相减，获得大气改正图形。图9是大气相位改正以及改正后的相位评估。由图9可知，大气相位去除本质上是剔除解缠相位中的大气差分相位，从图9(d)中可以提取绝对相位的工作面信息和工作面布点的值。从而可以将实际的水准测量数据与相位解缠数据进行对比，如图10所示。

图9 大气相位改正及改正后的相位评估Fig.9 Atmospheric phase correction and phase assessment after correction

图10 水准数据与各解缠相位提取的变形对比Fig.10 Comparisons between leveling data and deformations extracted from unwrapping phases

由图10可知，6种解缠方式获取的变形，都能够定性地反映出工作面变形情况，但是在变形量级上存在差距。其中，Lp最小范数法大部分不符合实际变形，将A、B线上所有的点作为定量分析的目标，以实际水准测量值为标准，可评价出6种相位解缠的精度指标，如表3所示。

表3 线分析指标Table 3 Line analysis index

由表3可以看出，Lp最小范数法获得的绝对相位各向指标均为最大，因此可以确定此法不适合矿区变形监测。其他五种算法都能很好地反映线性变形情况，其中本文方法要明显优于另外四种算法，Snaphu解缠线性性能要高于最小不连续算法，总体的线性精度排名为：本文方法>Snaphu>最小不连续>高金斯枝切法>最小费用流>Lp最小范数法。

3 结论

通过上述阐述，可以得到以下结论。

(1)融合GVF-Snake条纹探测和MRF图切两种方法对相位进行解缠，提出了边界探测“为真”的判别方法，将边界探测错误的内部缠绕块用MRF图切法相位解缠，最后按照解缠规则将外部解缠块与内部解缠块相融合，得到粗解缠相位(边界处含有孤立点)。

(2)将高通滤波插值法应用于解缠边界孤立点去除，得到了可观的滤波结果。

(3)将ITD大气相位模型应用到InSAR的大气去噪，通过分析显示，大气去除效果比较显著。

(4)解缠精度分析上应用了面分析指标和线分析指标，充分说明了各个解缠方法在矿区变形相位解缠方面精度性能方面的优劣性，并就其精度性能进行了排名。在面分析上，其精度性能排名为：本文方法>最小不连续法>Snaphu>Lp最小范数>最小费用流>高金斯枝切法；在线分析上，其精度排名为：本文方法>Snaphu>最小不连续>高金斯枝切法>最小费用流>Lp最小范数法。