融合分布式散射体时序InSAR技术在矿区形变调查中的应用

贾会会, 张海清, 李克达，张小朋

1.华北地质勘查局五一四地质大队，河北 承德 067000 2.河北省地质灾害监测预警技术创新中心，河北 承德 067000 3.河北大学建筑工程学院，河北 保定 071002

0 引言

河北省承德市滦平县张百湾镇周台子村铁矿区由于多年矿山开采导致历史遗留大量的浅埋采空区，其中部分采空区已发生明显的地面塌陷。其中：20170814，滦平县启星矿业集团启泰采区华兵选矿厂发生采空塌陷事件，形成深约40 m的塌陷坑，致使华兵选矿厂全部陷入地下，造成巨大的经济损失；201904—201905，距离周台子村约90 m的三采区连续发生两次采空区地面塌陷，塌陷坑深度约40 m，直径约50 m，造成矿路损毁断交。潜在的塌陷安全隐患会对周边建筑和居民区造成巨大威胁，因此对矿区进行形变调查和监测，对采空区塌陷进行灾害预警十分必要。

传统矿区形变监测的方法有水准测量、GPS测量、三维激光扫描[1-3]等。这些以“点”为基础的形变监测技术耗费大量人力、物力、财力，且监测结果空间和时间分辨率较低，因此无法进行大范围、高密度的矿区地表形变监测[4]。近几十年发展起来的InSAR(合成孔径雷达干涉测量技术)[5]具有监测精度高、全天时全天候、监测结果时间分辨率和空间分辨率较高等优势，其中监测精度可达到厘米甚至毫米级，为矿区形变监测和范围调查提供了新的有效监测手段[6-12]。传统时序InSAR技术采用永久散射体来获取地表形变信息，但这种方法的缺陷是自然地表测量点较少、空间采样率低，且很难确定出铁矿采空区的沉降范围和形变趋势。针对矿区内监测点数量不足的问题，本文采用DS(distributed scatterer，分布式散射体)和高相干点时间序列InSAR技术来对矿区进行形变解算。DS点是指分辨单元内具有相同散射特性的散射体，例如沙漠、草甸等。DS点处理流程包括同质像素点的提取和最优相位估计。提取同质像素点的现有方法是利用配准后的时间序列SAR(合成孔径雷达)振幅影像进行t检验、KS(Kolmogorov-Smirnov)检验、AD(Anderson-Darling)检验等[13]。这些方法虽然在效率和鲁棒性方面存在不足，但大多数都能有效地解决同质像素点的识别问题。现有的最优相位估计方法主要有MLE(极大似然估计)、PL(phase linking, 相位链接)算法和PT(phase triangulation,相位三角)算法[14-16]等。然而，这些方法不能在一个分辨单元中分解多个散射机制，否则可能会影响同质像素点相位的估计精度。由于中低分辨率SAR数据的分辨单元包括很多不同类型的地面目标；因此，可采用样本协方差矩阵或者相干矩阵EVD(特征值分解)方法进行单个分辨率单元中多个散射机制的区分，并进一步提高相位优化精度。此外，可基于分布式散射体的时间序列InSAR技术获得高密度监测点，提取铁矿开采区引起的地表形变不同时空特征，进而对矿区形变进行精确识别和监测。

为查明滦平县周台子村铁矿采空区空间分布，本文采取融合分布式散射体和相干散射体(distribute scatterer-coherent scatterer InSAR,DS-CSInSAR)的方法对研究区地表形变进行监测。即先采用AD检验的方法提取同质像素点，再采用基于协方差矩阵特征值分解算法来对DS点的最优相位进行估计，然后通过融合相干点来对矿区地表形变进行反演，并确定矿区形变的空间分布；以期为采空区塌陷灾害预警及治理提供技术支撑。

1 研究区和数据集

1.1 研究区介绍

研究区位于河北省承德市滦平县周台子村，铁矿采区位于滦河北岸。研究区的地理位置见图1a，范围为117.447°E—117.678°E，40.867°N—41.082°N，红色边框为研究区覆盖范围。该区地貌类型为低山侵蚀构造地貌(图1b)，高程范围为269～1 650 m；气候属大陆季风气候，冬长寒冷，夏短而炎热，多年平均气温9.1 ℃，最热月(7月)平均气温24.4 ℃，最冷月(1月)平均气温-9.4 ℃，极端最高气温41.5 ℃，极端最低气温-24.2 ℃，最大日温差23.8 ℃；历年最大降水量835.9 mm，最小降水量326.7 mm，平均降水量557.9 mm。研究区GPS监测点分别位于华兵矿业、JC33广场、JC34二采、JC35路边(图1c，d)，A—J铁矿区的位置见图1d。

a.研究区地理位置；b.研究区地质地貌图；c.研究区数字高程模型(DEM)图；d.铁矿和GPS监测点位置图。a图中，蓝色边框为Sentinel-1数据的覆盖范围，红色边框为研究区；b图中，红色边框为研究区；c图中，红色圆圈为该地区GPS监测点；d图中，A—J为铁矿区位置。

1.2 数据集

本文采用C波段48景250 km幅宽的Sentinel-1数据来对调查区进行形变反演。Sentinel-1数据集时间跨度为20170316—20201002，IW(interfe-rometric wide swath)TOPS(terrain observation with progressive scans，渐进扫描地形观测)模式，升轨轨道号为69，距离向和方位向分辨率分别为2.3 m和13.9 m。48景SAR图像的Sentinel-1数据时空基线图见图2，设置时间基线和空间基线的阈值分别为120 m和80 d，同时剔除一些受大气影响严重的干涉对，最终得到88个干涉对，其中最大的垂直基线和最大的时间基线分别为115.435 m和60 d。

图2 Sentinel-1数据时空基线图

2 技术方法

为了更精确地识别和监测调查区，首先提取出DS点并与相干点进行融合，得到研究区更精确的地表形变信息；然后使用沉降信息进行矿区的时空分析及沉降影响因素分析；最后确定出矿区的沉降区域。具体技术流程见图3，主要包括Sentinel-1 TOPS模式预处理、相干点和DS点处理、时序形变解算等。

γ.平均相干系数；γt.时间相干性；γth.平均相干系数的固定阈值。

2.1 Sentinel-1 TOPS模式预处理

进行预处理时，首先采用30 m分辨率的SRTM(shuttle radar topography mission，航天飞机雷达地形测绘使命)数据进行DEM(digital elevation model，数字高程模型)配准；其次，由于TOPS成像模式中多普勒中心频率的变化，因此需要采用ESD(enhanced spectral diversity，增强谱分集)配准方法来减少Burst间的相位跳变，并进一步优化方位向的偏移；然后，在将辅图像重新采样到主图像的坐标系中之前，需要对辅图像进行降噪去斜，同时对每个Burst都进行去斜、插值和反去斜的流程[17]，从而使SAR图像方位向上的配准精度优于1/1000像素[18]；之后，将每个条带的Burst进行拼接，因为研究区位于第2个条带，因此不需要进行条带拼接；最后按照研究区的范围进行ROI(感兴趣区)的提取。

2.2 相干点和DS点处理

本文利用研究区时间序列的相干系数图得到平均相干系数图，利用平均相干系数阈值法，即γ≥γth，提取平均相干系数γ在固定阈值(γth=0.7)以上的相干点作为高相干点。为了提取研究区的DS点，利用AD检验提取同质像素点，并利用特征值分解算法对同质像素点进行最优相位估计。

2.2.1 AD检验

相对于KS检验对尾部差异不敏感的缺点，AD检验在尾部增加了权因子，以有效减少第二类错误[19]。对配准后时间序列的SAR振幅图像，其中任意两个像素a和b的AD检验统计参数A2可表示为

(1)

式中：N为SAR图像个数；Fa(x)和Fb(x)分别为两个像素a和b的振幅的CDF(经验累积分布函数)；Fa,b(x)为a和b两像素点的CDF。当A2满足特定的阈值，则a和b为SHP(statistically homogeneous pixels，同质像素点)。

为了提取同质像素点，取13×13窗口作为本文的搜索窗口，使用上述的AD检验去比较该窗口内的中心像素点和其他像素点的检验统计量。本文中显著性水平取0.95，同质像素点阈值取20。当初始的同质像素点提取后，先判断其与中心像素点是否连接，舍弃与中心像素点没有连接的像素点，则剩下的像素点即为初始的DS点集。

2.2.2 同质像素点的最优相位估计

最优相位估计是SHP识别的关键步骤，本文基于协方差矩阵特征值分解算法对最优相位进行估计，该方法优势是不用进行相干矩阵求逆操作，避免了迭代耗时且计算效率高，且该方法适用于中低分辨单元内的多种散射机制[20]。

上述提取DS点的协方差矩阵C见下式：

(2)

式中：E为数学期望；p为时间序列像素x的复数向量，即p=[p1,p2,p3,…,pN]†；†为厄米特共轭转置；k为同质像素点集合的第k个点；Ω为Nx个同质像素点的集合。

因此基于协方差矩阵的特征分解可见下式：

(3)

(4)

由于分辨单元内主导散射机制数量n没有建立先验知识且不同特征值对应的特征向量相互正交；因此，求解最大特征值λi对应的特征向量vi的相位分量作为优化相位的估计值。

2.2.3 DS点与相干点融合

根据时间相干性γt对DS点进一步筛选，先将平均相干系数γ大于γt的DS点提取出来，即

(5)

式中：M为最优干涉相位值数目；φjk为原始干涉相位；θ为时间域滤波后求得的最优相位，j和k分别表示第j个和第k个干涉对;i为虚数单位。

然后，将上述提取的高相干点与DS点进行融合，同时将DS点最优相位估计与相干点的原始相位进行融合，进行时序形变解算。

2.3 时序形变解算

将DS点和相干点融合后，对这些点集构建Delanuary三角网。步骤如下：首先，构建线性形变模型进行邻近测量点的干涉相位差分，对三角网中每条边进行相对形变速率和相对高程误差解算，采用空间搜索方法对三角网中每条边进行求解，边的阈值取0.6，小于该阈值的边则被剔除，不参与形变解算，得到三角网每条边的相对形变速率和高程误差。然后，选择一个参考点，建立绝对形变速率与相对形变速率的系数矩阵，采用最小二乘算法对绝对形变速率和绝对高程误差进行求解。再将求出的形变速率和高程误差从干涉相位中减去，即可得到每个点的残余相位。最后，对残余相位进行MCF(minimum cost flow，最小费用流)相位解缠，得到其真实相位。残余相位包括大气相位、非线性形变相位和噪声相位，利用大气相位在时域上表现为高频信号、在空域上表现为低频信号的特性，采用时间域的高频和空间域的低频滤波对大气相位进行去除并估计出大气相位，将上述求得线性形变与非线性形变相加即可得到研究区的时序形变量结果。

3 结果与分析

3.1 相干点和DS点预处理结果

同质像素点的提取结果能够反映出同质区域同质点的分布情况，采用AD检验识别的同质像素点集见图4a。该区域大部分同质像素点的数目在100～255之间。图4b为利用相干性系数大于0.7筛选出来的高相干点，可见研究区由于位于自然地表区域，相干点非常少，仅为3 514个。图4c为利用时间相干性阈值为0.5筛选的最终的DS点，数量为57 766，图4d为DS点和相干点融合的结果。可见经过DS点加密，测量点的数量得到极大的提升，提高了时序形变解算的精度，且减小了相位解缠的误差。

a.同质像素点集分布图；b.高相干点目标分布图；c.DS点分布图；d.DS点与相干点融合图。b、c、d图中，红色点代表选出的点。

利用DS点所有干涉对的相位信息，采用基于协方差矩阵特征值分解算法对单一主图像干涉对的最优相位进行估计，以20181013—20181130和20181130—20181224两个干涉对为例，结果见图5。从图5a、5c可看出原始干涉相位去相干噪声较严重；对分布式散射体进行相位优化后，相位质量和整体的相干性得到明显的提高，在低相干地区，干涉相位的相干性也得到提高(图5b、5d)。

a.20181013—20181130原始相位；b.20181013—20181130最优相位；c.20181130—20181224原始相位；d.20181130—20181224最优相位。

3.2 铁矿区形变结果

通过本文的融合分布式和相干散射体时序InSAR处理后，得到铁矿区年平均形变速率(图6)，可知研究区的形变速率区间为-34.50～24.50 mm/a，形变趋势较大的地区主要集中在图6红色矩形框内。提取图6中红色边框的形变速率数据，结果见图7a。由该地区A—J铁矿区域红色曲线范围(图7b、c)可发现，这些地区的形变速率较大，且均出现不同程度的沉降。提取图7b和图7c中A—J铁矿区域中点目标(白色圆圈)进一步分析时序形变位移。图8为A—J形变区域提取的点目标时序形变位移曲线图，可发现：A—D区域在每年6—9月出现不同程度的沉降，沉降幅度在-33.40～32.70 mm，其中负值表示下沉，正值表示抬升；E区域在20200827沉降量最大，为50.76 mm；F—I区域在每年6—9月形变量变化区间为-35.02～26.03 mm，其中周台子村附近矿区F、G区域的最大沉降量分别为32.00 mm和34.00 mm，离路边最近的矿区H区域最大沉降量为33.00 mm；在J区域窑岭沟矿区出现明显的抬升现象，最大抬升量为24.03 mm。

红色矩形框内代表形变趋势较大的地区。

a.年平均形变速率图；b.铁矿区A—G形变区域；c.铁矿区H—J形变区域。

3.3 精度验证

为验证本文时序InSAR结果的准确性和精度，本研究获取滦平县矿区华兵矿业、JC33广场、JC34二采、JC35路边4个GPS监测点20200715—20201002期间24 h间隔的GPS实测数据和图8中GPS点对应的InSAR时序测量值进行对比验证。二者InSAR时序形变量与GPS的结果对比见图9，进而可计算出4个监测点与InSAR测量值之间的绝对误差分别为-1.70～3.70，-0.80～4.20，-3.42～3.77，-2.86～4.79 mm。InSAR与GPS监测相关性对比结果(图10)表明：JC33广场和JC35路边的相关系数分别为0.93和0.81，这是由于这2个GPS监测点位于建筑物附近，形变趋势稳定，且InSAR干涉图的相干性较好，在进行形变反演时受到的大气、地形误差影响较小，InSAR形变精度较高，因此与GPS测量值有着较高的相关性；华兵矿业和JC34二采相关系数分别为0.51和0.65，这是由于这2个GPS监测点位于矿区附近，属于自然地表，在InSAR处理过程中，易受到去相干的影响，形变解算精度没有建筑区的高。此外InSAR测量值与入射角和LOS(line-of-sight，雷达视线方向)、残余大气相位也有关系，因此与GPS测量值有着较低的相关性，但是它们的绝对误差范围较小，从而证明了InSAR监测滦平县矿区的地表形变结果与GPS野外观测结果一致性较好，结果可信度较高。

a.A—E区域；b.F—J区域。

a.华兵矿业；b.JC33广场；c.JC34二采；d.JC35路边。

a.华兵矿业；b.JC33广场；c.JC34二采；d.JC35路边。

为了对比分析DS-CSInSAR 技术结果的准确性，我们采用CSInSAR技术对该铁矿区进行了PSI(永久散射体雷达干涉)分析、Delaunay三角网构建、相对形变和绝对形变解算等步骤得到研究区的年平均形变速率，结果见图11。图11表明：两种方法得到的形变趋势基本一致，但与 CSInSAR 技术对比，本文提出的DS-CSInSAR方法极大地提高了点的密度；通过 DS 的空域滤波和最优相位估计，DS-CSInSAR方法提高了低相干铁矿区干涉图相位的质量和相干性，从而间接证明了该算法的可靠性。

图11 CSInSAR技术(a)与DS-CSInSAR技术(b)结果对比图

4 结论与建议

1)为弥补矿区监测点目标不足的情况，本文提出了DS-CSInSAR技术的铁矿区识别与监测的方法，利用双样本AD检验和协方差矩阵特征值分解方法进行相位优化，大大提高了矿区的测量点数量和形变解算的精度。与CSInSAR技术方法相比，DS-CSInSAR技术方法在监测点密度和干涉图的质量上均有明显改善。

2)通过对研究区地表形变的时空分析发现，研究区年平均形变速率为-34.50～24.50 mm/a。不同矿区在每年的6—9月出现不同程度的沉降，其中周台子村和路边附近F、G、H区域的最大沉降量分别为32.00、34.00、33.00 mm；窑岭沟矿区出现明显抬升现象，最大抬升量为24.03 mm。InSAR 监测结果与4个GPS站点结果一致，相关性较好，绝对误差为-3.42～4.79 mm。

3)本文利用高密度的DS点确定了铁矿区形变趋势，相对于其他技术具有空间分辨率高、监测范围大的优势。该方法较好地解决了在地形复杂铁矿区进行形变监测的问题，进一步提高了利用InSAR形变结果识别和监测铁矿区的能力；同时，根据沉降结果进行的分析也为预测铁矿区形变强度和位移趋势提供了重要信息。