InSAR相位解缠算法在矿区沉降监测中的对比与分析

王志勇，孟庆颖

(山东科技大学 测绘科学与工程学院，山东 青岛 266590)

1 引 言

随着地下矿产资源的开采，矿区地表沉陷问题日益严重，在矿区形成大量的塌陷坑、塌陷槽，引起农田破坏、房屋损毁等，严重影响人民的生命财产安全和经济的可持续发展。雷达干涉测量(Interferometric Sythetic Aperture Radar，InSAR)技术提供了一种全新的地表形变监测方法，具有全天候、全天时、低成本等优势[1-6]。相位解缠是InSAR数据处理及应用的一个关键技术和处理步骤，特别是在InSAR矿区沉降监测数据处理中，它直接影响InSAR地表形变监测的可靠性及精确性。

InSAR技术已经初步被应用于矿区地面沉降监测中[7]，但标志性的成果还很少，尚未在矿区沉降监测中得到大规模应用，这主要是由于矿区环境的复杂性以及InSAR技术的局限性决定的。矿区形变比较特殊[7]，主要表现为沉降速率不均匀、矿区形变量非常大(如在活跃期内，一个月的沉降量可达半米，甚至更大)、沉降范围分散；其次，沉降范围很小(多为半径几百米的沉降漏斗)；矿区农田多、地表植被覆盖率高，雷达数据相干性很差，致使在干涉条纹图中出现很多不连续的区域。基于这些特点，部分InSAR相位解缠在矿区沉降监测中不能正确解缠，致使InSAR监测的可靠性存在很大问题。

由于不同的相位解缠算法都具有一定的适应性，并不是所有的相位解缠算法都能解决矿区大形变沉降监测时的相位解缠问题。因此，针对矿区形变的特殊性，对现有的相位解缠算法进行对比分析，探讨适合矿区沉降特点的InSAR相位解缠算法，对于矿区沉降监测的可靠性及精度方面都具有十分重要的意义。

本文选取沉降特点非常明显的济宁市梁宝寺矿区为试验区域，采用真实的ALOS PALSAR(日本雷达卫星，2006年发射)雷达数据，生成差分干涉图。分别用经典的相位解缠方法进行矿区大形变区域的相位解缠实验，并对解缠结果进行分析评价，找出适于矿区环境的、精度较高、可靠性较好的解缠方法。

2 相位解缠算法分析

所谓的相位解缠就是将相位由主值恢复为真实值的过程[2]。在InSAR数据处理中，针对相位解缠已经进行了大量的研究，提出了许多算法，主要可以概括为以下三大类[3，8-9，14]：路径跟踪类相位解缠算法、最小范数类相位解缠算法、基于最优估计的相位解缠算法。

路径跟踪类算法主要采用不同的策略寻找最优化的积分路径，避开残差点造成的误差传递，以满足相位梯度闭合路径积分为零的条件，它一般是通过识别残差点，设置正确的枝切线阻止积分路径穿过；或者是在相位质量图的帮助下，从高质量数据开始积分。主要包括：枝切法[10]、掩膜切线法[9]、Flynn的最小不连续法等。

最小范数类相位解缠算法是一种全局求解算法，以缠绕相位的离散偏微分与解缠相位的离散偏微分的差最小为准则建立全局的数学模型来求解相位解缠的估计值。主要包括：最小LP范数法[2-3]、高斯-赛德尔迭代法、PCG预解共轭梯度法、FFT/DCT最小二乘法、多级格网法[2-3]。

基于最优估计的相位解缠算法将最优估计的算法应用到相位解缠中，主要包括：网络规划法[11]、Kalman滤波法[12]、遗传算法[13]等。

表1列出了几种经典的相位解缠算法的核心思想以及其优缺点。

表1 相位解缠方法及特点分析

3 研究区域及InSAR相位解缠实验

3.1 研究区域及实验数据

选取了我国矿区地面沉降比较严重的济宁矿区的梁宝寺煤矿作为实验区，选用2009年01月10日和2009年02月25日获取的两景ALOS PALSAR数据组成干涉对，其path号为449，Frame号为700，其工作模式为FBS(Fine Beam Single Polarization)模式，L波段(波长23.6cm)，产品级别为Level1.1(即单视复影像数据)，HH极化，视角为34.3°，方位向像元大小为3.148m，距离向像元大小为4.684m。干涉对的垂直基线为268.542m，时间间隔46d。

除此之外，为了去除地形相位的影响，本文采用了SRTM DEM数据。

3.2 InSAR技术生成差分干涉图

采用双轨法差分干涉测量[7]的方法生成差分干涉图，并采用改进的Goldstein滤波算法滤除了部分相位噪声。为使实验更具针对性，对原始单视复影像(SLC)进行区域裁剪，获得单个矿区的差分干涉图，采用1∶2(距离向：方位向)多视处理，生成的差分干涉图的大小为：470行×481列。图1(a)为梁宝寺矿区的滤波后的差分干涉图，图1(b)为相干图。

图1 差分干涉图及相干图

3.3 相位解缠对比实验

对滤波后的差分干涉相位图分别采用枝切法、区域生长法(Region Grow)、质量引导掩膜切线法、最小不连续法、LP最小范数法、无权多级网络法、加权多级网络法、PCG预解共轭梯度法、最小费用流法(minimum cost flow algorithm，MCF)进行解缠实验，解缠结果如图2所示。

图2 不同相位解缠方法解缠结果

4 相位解缠结果分析及评价

为了对比各种相位解缠算法在矿区沉降监测中的解缠结果，从剖面图、统计分析、运行时间、形变量等几个方面进行对比分析。

4.1 直接分析

根据图2，分析不同相位解缠方法得到的解缠结果发现，基于路径跟踪的相位解缠算法整体效果较差，只有Flynn最小不连续法解缠效果相对较好，这是因为枝切法在相位不连续点较多且比较密集的时候，无法进行解缠，区域生长法需要从质量较好的区域开始，质量引导掩膜切线需要提供质量较好的品质图像，矿区环境复杂，相干性较差，相位噪声严重，很难满足上述要求，只有Flynn最小不连续法可以在没有相位质量图指导的情况下较好地完成相位解缠；基于最小范数的相位解缠算法总体解缠效果不错，但是加权多级网络由于权重选择导致其在噪声较大的区域无法实现正确解缠；基于网络规划的相位解缠算法，最小费用流法(MCF)出现了较多不可解缠的区域，解缠效果较差。

4.2 剖面图

实验选取了其中的一行数据(位于差分干涉图的第210行)，对比分析每种相位解缠算法的结果图的剖面图，所生成的剖面图如图3所示。

根据图3，通过分析不同解缠方法下的剖面图，发现枝切法、最小不连续法 、LP最小范数、加权多级网络剖面线走势大致相同；区域生长法、最小费用流(MCF)剖面线出现大量不连续点；质量引导掩膜切线出现突变值，形成两个波峰，不符合整体走势。无权多网络、PCG预解共轭梯度虽然整体光滑，但是与其他剖面线有一定差别，可能在解缠过程中引入了误差，导致解缠结果出现问题。

图3 不同解缠算法的剖面图比较

4.3 定量统计比较及运行效率对比

对各相位解缠算法得到的解缠结果进行了统计分析，分别统计了其最大值、最小值、均值及标准差等，表2列出了各种相位解缠方法的统计结果。

结果显示加权多级网络和LP最小范数法的平均值和标准差较大，其他方法大致接近。在运行速度上，枝切法运行时间最短，LP最小范数运行时间最长。枝切法运行速度最快，但是稳定性较差；区域生长法稳定性较好，但是运行时间过长；无权多级网络运行较快且稳定性较好；最小不连续法在解缠效果、标准差、运行时间、直方图形状和剖面线走势上比较均衡，没有突出优势，但是整体较好；LP最小范数稳定性一般、运行最慢，但是解缠效果较好。

本文采用的9种相位解缠算法均运行在Windows XP环境下，电脑配置为Intel(R) Core(TM)2 Duo CPU2.80GHz、4G内存。各种解缠算法运行时间如表2所示。从表2中可以看出，不同的相位解缠算法的运行时间有很大差异，其中枝切法最快，而LP最小范数法由于是一种全局最优解，其运行时间相对较长，是所有算法中最慢的一种方法。

4.4 相位解缠对形变结果的比较

以上的统计及分析还不能确定相位解缠算法在矿区大形变量解缠时的优劣，还需要考虑解缠算法可靠性的问题。

在差分干涉图中，一个条纹代表了半个波长的形变量，因此，根据条纹的数量可以得到最大形变量的大小，通过不同相位解缠算法得到的矿区最大形变量与条纹数目直接反应的最大形变量的差值应该在较小的范围内，表3列出了各相位解缠方法得到的最大形变量与理论值之间的差值。

表2 相位解缠算法统计比较

表3 各种解缠方法最大形变量与理论值的差异

从表3中可以看出，尽管所有的相位解缠算法都能恢复相位的主值并且得到形变量，但在得到的最大形变量方面却存在较大的差异。并且有的算法与理论值相差较大，在这些算法中，区域生长法、无权多网络算法、PCG预解共轭梯度方法得到的形变量与理论值相差较小。

通过分析各种解缠方法解缠后的形变区域的直方图，发现区域生长法、最小费用流法(MCF)的直方图较平滑，分布较好；枝切法、质量引导掩膜切线、最小不连续法分布虽然较好，但是零值较高；LP最小范数、加权多级网络、PCG预解共轭梯度直方图分布不好且有毛刺；无权多网络出现两个峰值。

5 结束语

对梁宝寺矿区的差分干涉图，采用多种方法进行了解缠实验，并对解缠结果进行了定性和定量分析，通过分析发现：在矿区环境中，最小费用流法没有显著优势，反而出现了较多不可解缠的区域，而最小不连续法、LP最小范数、无权多级网络和PCG预解共轭梯度在解缠效果相对不错。其中PCG在解缠过程中引入误差导致直方图和剖面线与整体不一致，LP最小范数运行时间过长。综合分析得出最适于用于矿区的相位解缠方法是最小不连续法和无权多级网络法，它们解缠效率较高、精度较高，稳定性和可靠性较好，在矿区解缠实验中，受噪声影响较小，能较好地完成相位解缠。本文较全面地分析了在矿区环境中各种相位解缠方法的优缺点，为矿区相位解缠方法的选择提供了依据，对矿区沉降监测精度的提高具有重要意义。

致谢：本研究还得到了山东省泰山学者建设工程专项以及山东科技大学科研创新团队支持计划项目(2011KYTD103)的资助，在此表示感谢。

参考文献：

[1] HANSSEN R.Radar interferometry-data interpretation and analysis[M].Kluwer Academic Publisher，2001.

[2] 李平湘，杨杰.雷达干涉测量原理与应用[M].北京：测绘出版社，2006.

[3] 王超，张红.星载合成孔径雷达干涉测量[M].北京：科学出版社，2002.

[4] 张永红，张继贤，龚文瑜.基于SAR干涉点目标分析技术的城市地表形变监测[J].测绘学报，2009，38(6)：482-487.

[5] 程璞，许才军，王华.InSAR相位解缠算法研究[J].大地测量与地球动力学，2007，27(3)：50-55.

[6] 黄国满，张继贤，赵争，等.机载干涉SAR测绘制图应用系统研究[J].测绘学报，2008，37(3)：277-279.

[7] 王志勇，张继贤，黄国满.基于InSAR技术的济宁矿区沉降精细化监测与分析[J].中国矿业大学学报，2014，43(1)：169-174.

[8] 靳国旺，徐国华，余懋勋，等.基于瞬时频率估计的InSAR相位解缠[J].测绘科学技术学报，2009，26(1)：33-35.

[9] GHIGLIA D C.Two-dimensional phase unwrapping：Theory，algorithms and software[M].New York：John Wiley & Sons，Inc，1998.

[10] GOLDSTEIN R M，ZEBKER H A，WERNER C L.Satellite radar interferometry：Two-dimensional phase unwrapping[J].Radio Science，1988，23(4)：713-720.

[11] CHEN C W.Statistical-cost network-flow approaches to two-dimensional phase unwrapping for radar interferometry[D].California：Standford University，2001.

[12] 刘国林，郝华东，陶秋香.卡尔曼滤波相位解缠及其与其他方法的对比分析[J].武汉大学学报(信息科学版)，2010，35(10)：1174-1178.

[13] 赵争，张继贤，张过.遗传算法在InSAR相位解缠中的应用[J].测绘科学，2002，27(3)：37-39.

[14] 刘志敏，张景发，罗毅，等.InSAR相位解缠算法的实验对比研究[J].遥感信息，2012，27(2)：71-76.