融合多种差分干涉测量的矿区地表形变监测

乐颖， 夏元平*， 钱文龙

(1.东华理工大学测绘工程学院， 南昌 330013； 2.福州市勘测院， 福州 350108)

中国煤矿资源非常丰富，在开采过程中往往会出现高强度开采现象。而长期过度开采煤矿资源，将导致开采区域地表发生沉降，甚至可能对周边的居民区、道路和农田等造成严重影响[1-2]。目前，随着社会科技的进步，监测技术手段越来越多。合成孔径雷达干涉测量(interferometric synthetic aperture radar，InSAR)技术由于其全天时、全天候的特点，已经广泛应用于地表形变监测[3-4]、滑坡监测[5-7]、地形测绘、冰川测量[8]和地质灾害监测[9]等领域。PS-InSAR技术对数据数量要求较多，且监测结果与研究区域地貌密切相关，在植被较少的区域可取得较好的效果，而当地区覆盖较多植被时，取得的监测效果并不理想。煤矿开采区域往往位于山区，植被覆盖较多，且具有较大的开采量与较为明显的形变，易引起失相干问题。因此，PS-InSAR技术难以准确估算煤矿沉降中心形变值，但当PS-InSAR技术中的第一步反演所使用的控制点大多选择在形变小且相干性高的区域，符合SBAS-InSAR中参考点的选择要求，可将其直接应用于SABS-InSAR中作为GCP点从而提高轨道精炼的精度。而D-InSAR、SBAS-InSAR这两种技术方法在矿山形变监测领域应用较多。

为了避免单一方法的局限性，结合多种技术方法进行监测已经成为一种新的趋势。徐小波等[10]联合D-InSAR与偏移量追踪技术(offset-tracking)技术对高强度开采沉陷区进行了监测，并取得了良好效果。陈艳青等[11]利用D-InSAR和SBAS-InSAR技术监测大桥历史形变趋势，发掘雷达卫星的监测形变能力，为以后的监测预警提供了宝贵的经验。王勇等[12]融合D-InSAR与精密水准测量数据，利用线性回归构建校正模型，并应用于滑坡形变监测，从而提高了形变结果的精度。付波霖等[13]利用时序InSAR技术提取了西藏一滑坡点的形变信息，并对形变结果进行了交叉验证。柴华彬等[14]为解决因植被旺盛引起的InSAR失相干问题，提出了一种融合实测数据的SBAS-InSAR方法，通过对失相干区域采用反距离加权插值法，从而获取地表沉降信息。但是，大部分的学者更多的是根据不同技术得到的结果进行对比，或在缺少水准数据时通过交叉验证说明结果的准确性。并且，传统D-InSAR方法会受到时间失相关、空间失相关和大气效应等因素影响，从而限制了其在矿山形变中的应用；SBAS-InSAR可以弥补这些因素的制约，将两种技术融合或可取得更优的监测精度。鉴于此，现提出一种融合D-InSAR和SBAS-InSAR技术的矿区地表形变监测，将D-InSAR结果进行Kriging插值以补充SBAS-InSAR在矿区沉陷中心空缺部分，使得最终的矿区形变结果更加完整，并利用水准数据进行验证。

1 时序InSAR技术

1.1 基于相干性系数的 D-InSAR 地表形变信息提取方法

D-InSAR作为一门极具潜力的技术，基本思想是通过干涉相干性分析提取地表变化信息，并利用二次差分方法从干涉相位图中去除地形与其他因素的影响，从而达到提取形变信息的目的。“二轨法”差分干涉技术表达式[15-16]为

φ=φflat+φtopo+φdefo+φatm+φnoise

(1)

式(1)中：φflat为平地相位；φtopo为地形起伏造成的地形相位；φdefo为地表发生形变；φatm为由于大气发生延迟产生的大气相位；φnoise为观测时伴随的噪声相位。

总之，该技术最终目的就是将式(1)中φflat、φtopo、φatm与φnoise等相位信息剔除，最后只保留φdefo形变相位的技术。

相干系数于1993年最早由普拉蒂提出，该系数取值范围为[0,1]，取值越大说明相干性越高，理论模型表达式[17]为

(2)

式(2)中：E[·]为数学期望；u*为共轭复数；u1、u2为主、副图像的信号。然而，在实际的InSAR处理过程中，同一像素内很难获取计算所需数量的采样值。因此，根据式(2)的理论模型，基于SAR影像复数数据计算获得相干系数的标准表达式[18]为

(3)

式(3)中：u1(n,m)、u2(n,m)分别为主、副影像数据块内某个坐标(n,m)处的复数值；|·|2为对应的二阶范数；M、N为计算相干性的数据块尺寸大小；m、n为数据块内对应的行列号。

1.2 结合PS-InSAR的SBAS-InSAR技术

SBAS-InSAR干涉技术是由Ferretti等[19]提出的一种用于获取地表形变信息的时间序列分析方法，其基本原理是对相干目标进行相位分析，最终获取时序形变。通过选择合适的时空基线阈值构成差分干涉对，选定相干目标点，并利用线性相位变化构建模型与解算，通过时空滤波方法去除大气延迟，同时减弱了高程、大气误差以及D-InSAR处理过程中存在的失相关问题，从而获取地表的时间序列形变信息。

若SAR影像数据集中包含N+1景覆盖同一区域的影像，则影像获取时间序列[20]为

t=[t0,t1,…,tN]T

(4)

根据三基线和最小法选定其中一幅影像为主影像进行配准后，在全部自由组合差分干涉对中选取合适的时间基线和空间基线阈值干涉对，若获得M幅差分干涉图，则有

(5)

假设t0为参考时刻，则任意时刻ti(i=1,2,…,N)相对于t0时刻的差分相位φ(ti)为未知数，数据处理过程中所获取的差分干涉相位δφ(tk)(i=1,2,…,M)为观测量。若忽略失相关、高程误差以及大气延迟等因素影响，则第i幅(i=1,2,…,M)差分干涉图中像元(r,x)相位计算公式为

δφi(r,x)=φ(tA,r,x)-φ(tB,r,x)

(6)

式(6)中：λ为雷达波长；d(tA,r,x)和d(tB,r,x)分别为像元在主影像成像时刻tA和副影像成像时刻tB沿雷达视线方向的形变。若假设d(tA,r,x)=0，那么时间序列上的相位值为

(7)

式(7)中：d(ti,r,x)为任意时刻ti(i=1,2，…,N)相对于t0时刻沿雷达视线方向的形变。

图1 研究区域地理位置分布图Fig.1 Geographical distribution of the study area

2 融合D-InSAR和SBAS-InSAR技术的矿区地表形变监测

2.1 研究区域概况

峰峰矿区位于114°3′E～114°16′E、36°20′N～36°34′N[21]。峰峰矿区包括九龙矿、梧桐矿、新三矿、薛村矿、孙庄矿、牛儿庄矿、万年矿、通二矿、四矿等矿区，图1绘制了研究区域具体的地理位置分布情况。

现主要研究九龙矿区和万年矿区，从图1可以看出，这两大矿区的具体地理位置分别位于峰峰矿区的新坡镇和武安市的伯延镇。九龙矿隶属于冀中能源峰峰集团，是该集团有限公司名下的主力矿井之一。九龙矿位于太行山与华北平原的过渡地带，属于缓倾斜山前冲洪积地貌类型。矿区内地势西低东高，南低北高，总体形态为四周高中间低的盆地。万年矿区所在区域，四季分明，冬季寒冷干燥，夏季炎热多雨，地势低平，板块主体由南向北倾斜。值得注意的是，万年矿区和九龙矿区均位于农田下方，为了避免农田因测量而遭到破坏，同时考虑到观测站的稳定性，将水准测量观测站布设在附近的道路上[22]。

2.2 研究方法

融合D-InSAR和SBAS-InSAR技术的矿区地表形变监测，其目的主要是对D-InSAR和SBAS-InSAR这两种方法所得结果进行Kriging值，从而使得监测结果更加完整，以弥补SBAS-InSAR技术在监测矿区沉降中心存在空缺值的不足。在雷达影像数据处理过程中需要借助外部DEM，用来提供参考地形或者参考地理坐标系，常见的SRTM系统主要是30 m分辨率的SRTM1和90 m分辨率的SRTM3(https://srtm.csi.cgiar.org/)，本文选取的是SRTM 30 m数据。图2为了融合D-InSAR和SBAS-InSAR技术的矿区形变监测流程图。本文研究方法主要改进部分及具体操作步骤如下。

图2 融合D-InSAR和SBAS-InSAR技术的矿区形变监测Fig.2 Deformation monitoring in mining area by integrating DInSAR and SBAS-InSAR technology

(1)D-InSAR方法中的轨道精炼和重去平：由于轨道精炼时缺少形变先验认识及地面控制点信息，盲目选择控制点可能会引入误差。本文研究在选取D-InSAR方法中轨道精炼控制点时，引入相干系数法，该方法通过设置相干性系数阈值以剔除相干性低的像元点，从而保证选取的控制点相干性很高；再结合相位解缠结果，避免选取残余地形条纹上并远离形变区。

(2)SBAS-InSAR方法中主影像的选取：主影像选取对于差分干涉、相位解缠、估计平均速率等都有影响。因此，为了选取效果更好的主影像，本文研究采用三基线和最小法来选取主影像[23]。通过计算每景影像的时间基线、空间基线和多普勒基线，将影像集中这三种基线和最小的影像作为主影像，表达式为

(8)

(3)SBAS-InSAR方法中的轨道精炼：轨道精炼过程中需要选取若干参考点，该步骤主要由人为主观判断选择，并需要远离形变区和位于没有残余地形条纹上。若盲目选择控制点可能会引入误差，故此步骤对于操作人员的专业知识和经验能力要求极高。由于PS-InSAR技术中的第一步反演所使用的控制点大多选择在形变小且相干性高的区域，符合SBAS-InSAR中参考点的选择要求，因此可将其直接应用于SABS-InSAR中作为轨道精炼的控制点。利用地理编码将PS-InSAR技术中的参考点文件从地理坐标系转换为SAR坐标系，并将其作为SBAS-InSAR技术的轨道精炼控制点引入数据处理中，用于去除轨道数据所引起的误差。

表1 雷达影像数据的各个基线参数Table 1 Baseline parameters of radar image data

(4)D-InSAR结果和SBAS-InSAR结果的融合：本文采用的是Kriging插值法，将DInSAR方法得到的形变结果用于填补SBAS-InSAR结果在矿区中的空缺值，充分利用两种方法的互补性，以此获取最终的矿区地表形变结果。此外，Kriging法可对周围的测量值进行加权，从而得出未测量位置的预测，其表达式由数据的加权总和组成：

(9)

式(9)中：Z(si)为第i个位置的测量值；λi为第i个位置测量值的未知权重；S0为预测位置；N为测量值数。需要注意的是，在使用Kriging方法时，权重不仅取决于测量点之间的距离和预测位置，还取决于基于测量点的整体空间排列。若需要在权重中使用空间排列，必须量化空间自相关。因此，在普通Kriging法中，权重λi取决于测量点和预测位置的距离以及预测位置周围测量值之间空间关系的拟合模型。

2.3 监测结果评价指标

本文研究以标准差和均方根差作为精度指标，对研究区域形变结果精度进行定量评价[24]。标准差(standard deviation，SD)，是离均差平方的算术平均数的算术平方根，是测绘领域经常用来统计分布程度测量依据，用σ表示，σ越小表明数据离散程度越小，稳定性越高，其计算公式为

图3 2015-06-17—2015-10-03间的DInSAR形变结果Fig.3 DInSAR deformation results between 2015-06-17 and 2015-10-03

(10)

均方根差(root mean square error，RMSE)是一种常用的测量数值之间差异的量度，值越小说明观测值与真值之间误差越小，所得结果精度越高，其计算公式为

(11)

许诺是在电视上看到海岛发生地震的新闻后，匆忙放下手头的工作赶过来的，他提着一颗心飞了一路，看到丁小慧毫发无损地站在他面前，他才松了口气，紧紧地抱了抱她，“你是想吓死我啊？”

3 形变结果分析

3.1 融合D-InSAR和SBAS-InSAR技术的矿区地表形变监测结果

3.1.1 基于相干性系数的D-InSAR地表形变结果

图4 2015-06-29—2015-10-27间的DInSAR形变结果Fig.4 D-InSAR deformation results between 2015-06-29 and 2015-10-27

以2015年6月17日—10月3日和2015年6月29日—10月27日数据为例，将25景影像分别进行基于相干性系数的D-InSAR处理，图3和图4分别为两个时间段的数据所得形变结果图。为了更清晰地体现研究区域的形变结果，将万年矿区与九龙矿区形变结果进行局部放大，两者局部放大形变结果分别见图3和图4右侧的上方和下方。从图3和图4中可以明显看出，2015年6月17日—10月3日期间的形变量变化范围为-81.51～72.88 mm；2015年6月29日—10月27日期间的形变量变化范围为-130.46～77.06 mm。此外，通过对比分析图3和图4，能够发现两者地表形变量变化范围具有细微差异，主要是由于两者使用不同时间段的影像所致，但最终形变结果的沉降分布情况仍基本保持一致。从图3和图4的局部放大图中可以明显看出，万年矿区有两处沉降区域，九龙矿区有一处沉降区域。

图5 矿区的形变速率图Fig.5 Deformation rate diagram of the mining area

3.1.2 结合PS-InSAR的SBAS-InSAR技术的矿区形变结果

SBAS-InSAR技术相比于D-InSAR技术具有时序分析的优势，不仅可获得地表形变量结果，还可得到年均形变速率。将PS-InSAR技术与SBAS-InSAR技术相结合，并对25景SAR数据进行处理，可得到研究区域2015年6月—2017年3月间的年均形变速率，图5绘制了研究区域的形变速率结果。

从图5的左侧整体图可以看出，峰峰矿区较多煤矿开采区域存在地表沉降现象，年均形变速率最高达到了35 mm/a；而图5右侧局部放大区域分别绘制的是万年矿区和九龙矿区的形变速率变化情况，能够明显看出由煤矿开采引起的沉降特征呈漏斗形状，位于沉陷中心位置的沉降量最大，沉降速率也最快，且年均沉降速率最大可达-103.82 mm/a，符合矿区的沉降特征。

3.1.3 融合D-InSAR和SBAS-InSAR技术的矿区地表形变结果

图6 万年矿区形变时序结果Fig.6 Time series results of mine deformation in Wannian

将25景Sentinel-1影像根据3.1.1节与3.1.2节中的两种方法分别得到矿区形变结果，两者形变结果中研究区域的沉降趋势变化和沉降位置分布基本一致。利用Kriging插值法将这两种方法的形变结果进行融合得到矿区最终的地表形变结果，综合分析可反映两个矿区的沉降分布情况和沉降变化特征，图6和图7分别绘制了万年矿区和九龙矿区在2015年6月—2017年3月间的形变时序结果，其中正值代表地表向上发生抬升，负值代表地表向下发生沉降。从图6和图7中可以看出，地表形变在空间上的分布情况是随时间变化的，且累积沉降量随时间逐渐增加，并逐渐呈现漏斗形状，最大累积沉降量可达176.39 mm。总体而言，万年矿区和九龙矿区的地表沉降变化特征符合煤矿开采的变化特征，且融合后的累积沉降量和SBAS-InSAR得到的年均形变速率中的地表沉降分布情况吻合。

图7 九龙矿区的形变时序结果Fig.7 Time series results of mine deformation in Jiulon

3.2 融合形变结果分析

相干系数图反映了两幅影像的相关程度，可以用于评价SAR影像精确配准后所得到干涉条纹图的质量，主要是通过计算干涉条纹图的相干系数实现。并有研究表明，相干系数图质量的好坏还与季节的变化有关，主要是由于植被覆盖的影响，夏季的相干系数图质量普遍偏低，冬天相干系数图质量通常较高。根据1.1节中相干系数理论知识可知，相干系数的值越接近1，则两幅影像的相干性越高。矿区开采往往伴随着大范围的地表沉降变化，可以将其与光学影像相结合进行露天矿区的开采识别[25]。图8为研究区域的相干系数图，从图8中可以看出，除了部分水体和矿区大量级开采区域的相干性较低，其余整体的相干性都较高。从而说明，相干系数图的质量较优，所得地表形变结果准确性较高。

图8 研究区域的相干性系数图Fig.8 Coherence coefficient diagram of the study area

3.2.2 典型区域形变结果剖面分析

(1)万年矿区的剖面分析。为了研究矿区沉降漏斗的演化过程，将万年矿区沿着开采面的横纵剖面分别进行形变结果分析，图9绘制了地表形变横、纵剖面曲线结果。从图9中可以明显看出万年矿区横纵剖面的沉降情况，2016年10月9日、2016年12月20日和2017年3月2日这三景SAR数据的沉降趋势高度一致，在132610工作面开采区出现了明显的沉降漏斗。其中，图9(a)中，横向剖面结果显示有一个沉降中心，最大沉降量达到了88 mm左右；图9(b)中，纵向剖面结果显示有两个沉降中心，这与时序形变结果显示的沉降漏斗基本一致，最大沉降量达到了98 mm左右。

(2)九龙矿区的剖面分析。与万年矿区研究沉降漏斗的演化过程一样，将九龙矿区沿着开采面的横纵剖面分别进行形变结果分析，图10绘制了地表形变横、纵剖面曲线结果。从图10中可以观察到九龙矿区横纵剖面的沉降情况，2016年10月9日、2016年12月20日和2017年3月2日这三景SAR数据的沉降趋势高度一致，在15235工作面开采区出现了明显的沉降漏斗。其中，图10(a)中，横向剖面结果显示有一个沉降中心，最大沉降量达到了126 mm左右；图10(b)中，纵向剖面结果显示有一个沉降中心，最大沉降量达到了110 mm左右。

图9 万年矿区地表形变剖面曲线图Fig.9 Surface deformation profile curve of Wannian mining area

图10 九龙矿区地表形变剖面曲线图Fig.10 Surface deformation profile curve of Jiulong mining area

3.3 融合形变结果评估

表2 水准数据和形变结果误差对比Table 2 Comparison of errors between leveling data and deformation results

4 结论

为了确保结果的准确性，采用了多种InSAR技术对研究区域进行形变监测。首先根据PS-InSAR方法的GCP点改进了SBAS-InSAR方法中轨道精炼控制点选取的问题，并利用Kriging插值法将基于相干系数的D-InSAR方法处理结果填补SBAS-InSAR技术中的空缺值，使得最终得到的矿区形变结果更加全面、有效。

(1)通过基于相干性系数的D-InSAR地表形变结果对矿区进行形变监测，结果表明，2015年6月17日—10月3日期间的形变量变化范围为-81.51～72.88 mm；2015年6月29日—10月27日期间的形变量变化范围为-130.46～77.06 mm。并测得万年矿区有两处沉降区域，九龙矿区有一处沉降区域。

(2)利用Kriging插值法融合D-InSAR和SBAS-InSAR技术的矿区地表形变结果，综合分析可反映万年矿区和九龙矿区这两个典型矿区的沉降分布情况和沉降变化特征。结果表明，地表形变在空间上的分布情况是随时间变化的，且累积沉降量随时间逐渐增加，并逐渐呈现漏斗形状，最大累积沉降量可达176.39 mm。

(3)将本文结果与水准测量数据进行对比分析。结果表明，通过本文方法得到的形变结果和水准处理数据基本一致，且误差均小于13 cm，监测精度满足煤矿开采相关测量技术要求。

融合D-InSAR方法和SBAS-InSAR方法可以避免单一方法所产生的偏差，综合利用两种方法的优势可弥补另一种方法中的不足，从而提高形变监测结果的准确性。但是对于大形变量沉降的矿区监测，由于C波段可监测形变梯度的限制，将在今后的研究利用其他波段的卫星影像对矿区进行形变监测。