矿区沉降监测的InSAR质量图相位解缠方法

刘万利 张秋昭 胡 江

（1.徐州工程学院数学与统计学院，江苏 徐州 221018；2.中国矿业大学环境与测绘学院，江苏 徐州 221116）

二维相位解缠是合成孔径雷达干涉测量技术（Interferometric Synthetic Aperture Radar，InSAR）中的一个关键环节，解缠的成功与否直接影响到目标高程信息的准确性，因此一直是InSAR技术研究的热点和难点［1-2］。近年来，许多学者对相位解缠进行了深入研究，主要可分为两大类：一类是基于积分的路径跟踪算法，另一类是最小范数算法［3］。质量图在这两类算法中起着非常重要的作用，直接影响其解缠精度［4］。尤其对于路径跟踪算法，一个好的质量指标通常可以避免或减少积分路径穿过残差点，可以使得高质量区域获得比较高的精度［5］。

基于此，许多学者就如何选择合适的质量指标展开了相关研究［6-14］。GHIGLIA和PRITT在文献［12］中比较了相干系数（Correlation Coefficient，CC）、伪相干系数（Pseudo Correlation Coefficient，PCC）、相位导数变化（Phase Derivative Variation，PDV）、最大相位梯度（Maximum Phase Gradient，MPG）4种质量指标，并指出相干系数是最有效的一个。OSMANOGLU等［13］描述了6种质量指标，试验结果显示Fisher距离的性能总体上优于相干系数、相位导数变化、结合枝切法的相位导数变化、二阶导数可靠度及线扫描指标等5种质量指标。LIU等［15］提出一种基于灰度共生矩阵（Gray Level Co-occurrence Matrix，GLCM）的熵差质量指标，该指标对分级级数要求比较高，需要事先确定级数值。刘万利等［5］对几种常见的质量指标做了比较深入的分析研究，并提出了一种改进的Fisher距离指标，试验结果表明，Fisher和改进的Fisher信息质量指标总体上比较稳定。

综上分析，基于相干系数和相位导数变化的质量指标简单易行，相干系数强调时间相干性，相位导数变化强调空间相似性［5］；并且上述质量图指标大多单方面偏重时间相关、空间相关、噪声等某一个指标。InSAR技术在监测山区矿区地表形变时，获得的干涉图通常具有噪声大、条纹密且复杂等特点，单一质量图指标用于山区矿区等复杂区域干涉图相位解缠时，适应性较差。为此，本研究提出一种融合相位导数变化和相关系数两种质量指标的新的质量图指标。该指标包含了空间相关性、时间相关性等更多的质量信息，因此更适用于地物覆盖复杂地区、条纹比较密集干涉图的相位解缠（尤其是山区煤矿区沉降监测干涉图）。采用模拟数据和实测数据试验验证了所提质量指标的可靠性和适用性。

1 相位解缠理论

基于滤波的相位解缠方法在展开密条纹高噪声干涉图相位时具有独特的优势［5，15］，本研究采用容积卡尔曼滤波相位解缠方法［16］进行分析。

1.1 相位解缠模型

容积卡尔曼滤波相位解缠模型在继承卡尔曼滤波相位解缠模型的基础上，采用容积卡尔曼滤波器处理非线性方程，其状态方程和观测方程为

式中，y(k)和x(k)分别为像元k的观测值和解缠相位；φ(k)为相位梯度估计值，通过最大似然估计器获得；h[⋅]表示观测函数；v(k)为k像元的观测误差向量；w(k)为梯度估计误差，被认为高斯白噪声且满足如下条件：

式中，E[⋅]表示期望；D[⋅]表示方差；为梯度估计误差方差，其计算参考文献［5］。v1(k)和v2(k)分别为复观测值实部和虚部的观测误差，被认为是高斯白噪声且满足如下条件：

式中，R(k)为观测误差协方差阵；，SkSNR为像元k的信噪比。

容积卡尔曼滤波相位解缠是沿着某一特定路径进行的，而路径的选择是基于某一质量指标从高质量像元到低质量像元。质量指标选择得当与否直接关系到解缠效果。

1.2 相位解缠步骤

本研究相位解缠实施流程如图1所示。具体步骤为：①根据某一质量指标生成质量图；②选取质量图中的一个高质量点A1作为解缠始点；③检测像素点A1的4个邻接像素点；④利用容积卡尔曼滤波解缠模型估计4个邻接像素点中最高质量的像素点A2；⑤选取像素点A1、A2邻接的质量最高的未解缠像素点A3并解缠。重复以上步骤，直到相位图中的全部像素点都被解缠完毕。

2 基于相干系数与相位导数变化的新质量指标

相干系数和相位导数变化是两种较常用的衡量干涉图像质量的质量指标，分别介绍如下。

2.1 相干系数

相干系数是最直接、最常用的衡量［12］指标，计算公式为

式中，s1=c1+n1，s2=c2+n2，c1、c2为图像的信号部分，n1、n2为图像的噪声部分；[⋅]*为共轭运算。

对于实测干涉图，CC是综合多种失相干因素得出的结果，表征的是两幅参与干涉的影像的相干性。通常情况下，反射能力强的地物对应像元的CC值较高。本质上是用相位的一阶导数来描述相位质量，对相位梯度比较敏感。因此，对于地形陡峭的区域，通过计算相位梯度，会将其标识为低质量数据。

2.2 相位导数变化

对于二维干涉图数据，某像元的PDV值计算公式为［12］

由式（5）可以看出，PDV是以当前像元为中心的一定窗口内的像元在水平方向和竖直方向上梯度方差大小之和，反映的是目标窗口内梯度数值的离散程度，是梯度空间相似性好坏的评价指标。本质上是用相位的二阶导数来描述相位质量，对梯度变化比较敏感。因此在一些相干性较低的同质地物覆盖区，该指标会表现出高质量结果。

2.3 PDV和CC融合质量指标

由2.1及2.2节可以看出CC可以很好地反映像元的个体特征，但不能很好地反映空间相似性，对于一些存在孤立高相干点的干涉图，容易造成解缠路径穿过低质量区，去解缠那些高质量点然后再穿过低质量区去寻找次高质量点，从而造成误差过早累积。而PDV会将一些空间上匀质但相干性不高的区域误判为高质量区。基于此，本研究提出了一种融合PDV和CC的质量指标。运算步骤如下：

（1）数据标准化。由式（4）和式（5）获得的CC和PDV值量纲和数量级不同，首先对其进行标准化处理可得：

由式（7）可知，PDV_CC值越小，质量越高，该指标融合了像元相位的空间相似性和时间相干性两个特征，是一种信息比较全面的质量指标。

3 试验与分析

我国多数煤矿区地表下沉速度快、植被覆盖多，尤其是山区矿区地形起伏大，因此InSAR干涉图通常呈现出条纹密集复杂、噪声大、残差点多等特点。如果质量指标选用不当，会造成解缠路径穿过大量残差点，进而导致沉降监测精度不理想。本研究采用一组多山场景模拟数据（干涉图具有条纹密集、噪声大的特点）和矿区实测数据对所提出的质量指标进行试验分析。为详细比较试验结果，整个相位解缠过程只有质量指标不同，解缠算法采用容积卡尔曼滤波相位解缠方法。采用解缠误差d和偏差χ2两个指标进行定量评价，两者计算公式为

式中，φu为解缠相位值；φr为相位先验值；σ为标准差；N为样本数。

3.1 模拟数据试验

多山场景地形模拟数据如图2（a）所示。为了详细分析地形起伏对相位解缠的影响，模拟数据的垂直基线设为150 m，像元个数为250×250。数据具有地势起伏大、陡峭度较高等特点。此场景的模拟成像加噪干涉图如图2（c）所示。本试验的含噪干涉图体现的是几何失相干程度，几何失相干严重的区域噪声强度大，由干涉图可以看出多个区域条纹比较密集。

图3给出了模拟数据的3种质量图。从大致的纹理来看，CC质量图与其余两种质量图差别较大，高低质量区区分比较明显。这是因为模拟数据CC值的计算仅体现了几何失相干，条纹密集的区域质量低，稀疏的区域质量高。PDV_CC与PDV的大致纹理比较接近，但颜色有所区别，前者颜色介于CC和PDV之间。这是由于PDV_CC是融合了CC和PDV两种质量指标所致。另外，由图3还可以看出在条纹稀疏区，PDV_CC和PDV也有少量低质量散点，是因为含噪干涉图的噪声添加具有随机特性，而这两者的计算都是基于含噪干涉相位的。

图4给出了3种质量指标引导解缠路径的相位解缠结果图。相对于真实相位图（图2（b）），可以看出CC引导的解缠失败（图4（a）），PDV（图4（b））和PDV_CC（图4（c））解缠结果总体保持了原样，但在条纹非常密集的区域解缠结果表现出了不连续现象，即解缠相位误差较大。这是因为这些区域噪声强度大、图像质量低、条纹密集且变化较快，进而导致条纹梯度的估计误差较大。进一步观察可以发现，在条纹密集的区域，PDV_CC较PDV解缠结果更连续光滑（主要体现在椭圆和矩形区域）。说明兼顾空间相似性和时间相干性的PDV_CC质量指标优于仅考虑空间相似性的PDV质量指标。

为了进一步讨论所提质量指标的可靠性，图5给出了3种指标的解缠误差图及其对应的误差直方图。从误差图可以明显看出CC解缠误差图右侧部分出现了大面积的大误差区域，PDV和PDV_CC解缠误差图只是在条纹密集噪声较大的区域存在较大误差，这说明CC指标最不可靠。对比图5（b）和图5（c），可以看出PDV_CC方法解缠图中椭圆和矩形区域的不连续区域明显小于PDV方法，这说明PDV_CC质量指标优于PDV。从误差直方图也可以看出PDV_CC的误差分布图形中心峰（落在0附近部分）高度最高且宽度最窄，说明该指标引导的解缠结果误差方差最小。因此PDV_CC的解缠结果优于其他两种指标。

表1给出了模拟数据3种质量指标性能的定量评价指标取值，分别为误差均值、误差方差以及偏差均值。从解缠误差和偏差表达式可知，以上3个定量评价指标的值越接近0说明解缠结果越好。因此从表1可以看出CC指标最不可靠，PDV_CC优于CC，是一种比较鲁棒可靠的质量指标。

分析偏差值随解缠路径的变化情况可以很好地解释误差从哪里产生、误差随路径如何传播等问题。一个好的质量指标应该是较大的偏差值出现得越晚越好，这样可以避免误差较早累积。模拟数据解缠结果的偏差值随解缠路径的走向曲线图如图6所示。可以看出，CC的偏差值最大，PDV次之，PDV_CC最小；CC方法的大误差出现最早（在3 700个像元左右），PDV次之，PDV_CC方法的大误差出现最晚（在56 000个像元左右）。这说明PDV_CC方法出现误差累积现象较晚，有更多的像元解缠结果较可靠。

3.2 矿区实测数据试验

研究区为神东矿区大柳塔煤矿52304工作面，像元个数为300×250，区域内地形略有缓坡，有荒草沙石覆盖。植被类型有低矮灌木、沙生植被、荒草等。干涉图由2013年4月2日及2013年4月24日两景TerraSAR-X 影像生成（图7（a）），外部DEM采用课题组通过收集的航拍地形图资料并结合GPS控制点生成的Relief-DEM。干涉图的生成采用Gamma软件。由于地表覆盖导致的失相干较严重，故干涉图采用的是2视处理结果。3种质量指标对应的质量图见图7（b）至图7（d）。由于试验采用的是沉降监测干涉图，整个区域的真实沉降值难以获得，因此解缠结果主要从连续性角度进行分析。

由图7（b）至图7（d）可以看出，CC质量图的高质量区比较小，而PDV和PDV_CC显现出大面积高质量区，这是因为后两个指标考虑了空间相似性，而研究区地面覆盖空间相似性较高。PDV_CC与PDV质量图相似度较高，但比后者颜色稍淡，与公式所体现的特征相符。由图7（e）至图7（g）解缠相位图可以看出，CC质量指标解缠结果在形变区右下角出现了较大的相位突变，相位连续性最差。在此区域，PDV和PDV_CC两种方法解缠效果相当，均优于CC解缠方法。结合干涉图和质量图可以看出，此区域噪声大并且失相干比较严重，说明在低质量区PDV和PDV_CC指标鲁棒性较好。在形变区右上方，信噪比较高，CC和PDV_CC取得了比较连续的解缠结果，而PDV方法出现了相位跳变。这说明PDV指标在引导解缠路径时出现了从低质量到高质量的情况，导致误差提前大量累积从而造成高质量区解缠结果不理想。从3种质量指标引导的解缠顺序图（图7（h）至图7（j），尤其是图7（j）椭圆区域）可以看出，解缠结果的不同是由于质量指标引导的路径不同导致的误差累积不同所致。由此试验可以看出本研究提出的PDV_CC质量指标鲁棒性较好。

综合以上试验结果可以看出，无论是模拟数据还是实测数据，PDV_CC都取得了优于PDV和CC的解缠结果。这说明融合指标可靠性好，是一种比较适合引导复杂地物覆盖、沉降较快工作面干涉图相位解缠的指标。此外，对于模拟数据，CC指标解缠结果极其不可靠，是因为模拟的地形在小范围内起伏变化较快，存在高质量孤岛；CC指标引导的路径出现了从低质量到高质量情形，从而造成了误差累积传播。对于实测数据，在沉陷区PDV_CC和PDV取得了相似的解缠结果，较符合实际沉降情况，但在非沉陷区PDV出现了解缠相位不连续现象，而CC和PDV_CC结果相似且较符合实际情况。两种试验结果有所差别，主要是因为模拟数据的噪声产生原因比较单一，主要由几何失相干引起，而实测数据噪声原因比较复杂。

4 结论

（1）提出了一种融合PDV和CC的新质量指标，该指标不仅考虑了时间相干性还兼顾了空间相似性，更适用于噪声大、条纹密集的复杂干涉图相位解缠。

（2）高山地形模拟数据和矿区实测数据试验均表明，新的质量指标与单一质量指标相比，具有较好的鲁棒性和实用性。

（3）相位梯度估计噪声与质量指标的相互影响机制是下一步研究的方向。