时间序列InSAR技术在矿区地表形变监测中的应用——以内蒙古霍林河露天矿区为例

王洪明，李如仁，覃怡婷，刘竹青，顾 骏

（1.沈阳建筑大学交通工程学院,辽宁 沈阳 110168；2.扎鲁特旗扎哈淖尔煤业有限公司,内蒙古通辽 029100；3.沈阳市勘察测绘研究院,辽宁 沈阳 110004）

0 引言

露天矿的开采导致矿区生态遭到不同程度的破坏，矿区及其周边的地质灾害主要表现形式是地表形变，同时矿区的形变对城区地表也有不同程度的影响[1−2]。传统的监测手段（如水准测量、GPS 等）对于矿区形变提取有一定的局限性，水准测量方法过于依赖人工，在连续性和实时性方面略有不足；GPS 方法技术成本过高，且对于面域的沉降范围有较大的限制。合成孔径雷达可对地表进行全天时、全天候对地观测[3]。

雷达干涉测量技术可以对露天矿区及其周边地表形变进行监测，具有大范围和高精度等优点[3]。InSAR的形变监测发展较迅速，最初的D-InSAR技术只能获取地表短时间的形变，后来采用时序InSAR技术提取长时间跨度的形变[4−6]。传统的时序InSAR技术对于矿区及其周边地表的形变提取存在明显的不足，由于露天矿地表的相干性差，故只能提取分布式点目标，而城区地表具有足够相干性可以提取永久散射体目标。刘一霖等[7]利用时序研究开采进程中地表大量级沉陷的完整形变时间序列，得到开采工作面地表形变的时空演变规律；仝云霄等[8]利用经典的D-InSAR技术对矿区沉降进行了精细化分析，与最邻近时间的精密水准监测结果总体趋势吻合，为矿区整体规划、灾害预警和开采沉陷治理等提供参考；张学东等[9]以唐山市为例，使用相干目标短基线InSAR技术监测矿业城市地面沉降，查明了其地面沉降量及其空间分布特征，最大年沉降速率达到−46.8 mm/a，主城区沉降速率普遍低于−11 mm/a；任文静等[10]、蒋金雄等[11]利用小基线集（SBAS）方法和相干矩阵特征值分解（T-EVD）的DSInSAR技术分别对矿区进行地表沉降提取，为矿区地表形变监测提供了良好的应用前景；谢文斌等[12]利用时序InSAR技术对抚顺市地表形变信息进行提取，结合信息熵方法分析地表不均匀形变现状。总体而言，目前在矿区及其周边城区地表形变提取及分析研究较少。随着露天矿的逐年开采和城市的迅速发展，矿区沉陷及矿业城市地表沉降是亟待探究的复杂问题。

传统的时序分析方法误差较大，也无法整体评估出矿区形变对其周边城区地表形变产生的影响。鉴于此，基于露天矿区及其周边城区地表的相干性差异巨大的特点，文章提出改进的时序InSAR技术，此方法在点目标提取的基础上，使用狄龙尼三角网将分布式目标与永久散射体目标进行统一构网，进行城市与矿区的形变一体化反演。以霍林河矿区及其周边城区为试验案例，获得了露天矿区及其周边城区地表的形变速率。

1 研究区与数据

研究区霍林河煤田地处大兴安岭南段脊部，是一个山间盆地，四周为中低山峦所环抱（图1），海拔标高1 100 ～1 350 m。盆地内部地势较平坦，东北和西南两端为低丘陵，中间区段是开扩的平原，海拔标高930 ～980 m。研究区内部包括中部与北部两个露天矿和霍林郭勒市市区。为了较好覆盖研究区，采用29 景Sentinel 1 降轨数据及90 m 分辨率Version4 SRTM 数据，Sentinel 1 数据具体参数如表1。

表1 研究区SAR 数据参数表Table 1 SAR data parameter table in the research area

图1 研究区Google earth 影像图Fig.1 Google Earth image of the research area

2 时序InSAR 处理

由于露天矿区及其城区周边地表相干性差异巨大，为保持矿区周边具备一定的相干性，采用改进的时序InSAR技术，该方法通过严格筛选差分干涉数据，对差分干涉数据提取分布式散射体与永久散射体，作为形变分析的目标点，通过进行回归分析与模型精化，得到研究区地表形变速率。改进的时序InSAR技术既继承了传统时序InSAR 的技术优势，又很好的解决了目标研究区域相干性差异性较大的问题，具体技术流程如图2所示。

图2 改进时序InSAR技术流程图Fig.2 Flow chart of improved sequence InSAR technical

2.1 预处理与差分干涉图生成

通过GAMMA 预先设置时间基线为30 d，空间基线为临界基线的30% 进行干涉生成，为了更好的筛选干涉数据[11−13]，选用SRTM 数据除去差分干涉对原始的地形相位，通过对差分干涉数据的易发生形变区域进行目视分析比较剔除干涉效果差的像对，最后选择34 个差分干涉结果进行后续试验处理，且为了保证低相干区域结果的准确性对9 个相对组成的子集进行3D 解缠，低相干地区参照3D 解缠区域进行时序处理。筛选出部分的差分干涉结果如图3所示，6 张差分干涉图在矿区位置都有明显的形变漏斗产生。

图3 差分干涉结果图Fig.3 Differential interference result diagram

2.2 干涉点目标选取

由于矿区及其周边相干性较低，导致干涉结果较差，需要对矿区及其周边分布式目标点进行提取，文章提出经过参数优化的分布式散射点选取方法对研究区进行选点，主要包括以下三步：

（1）将同一区域分干涉结果像素数据进行时间维度平均，采用24×24 窗口进行计算，将其中心区视为均值即为参考值，邻域像素值作为待估计值，为了减小误差，考虑空间临近关系，只把临近点像素作为待估计值[14]。

（2）运用瑞丽分布函数获得置信域区间，函数如式（1）：

式中：γref——参考像素均值；

γarb——任意临近像素均值；

z1−α/2——标准正态分位点，服从标准正态分布，为了增加同质点样本数量，α预先设置值为50%；

P——概率函数模型；

L——多视视数，得到置信域区间并获得同质点[15−16]。

（3）为进一步提高样本选取精度，在小样本数据情况下，传统的高斯假设不能够逼近真实的分布，可以根据SAR 强度服从指数分布的性质，其累加服从伽马分布，可以得到更加精确的置信域区间。

式中，σ代表SAR 强度指数分布的期望，根据伽马分布性质可知其分位点为 2/α，在有N 个样本的情况下，则可得置信域区间为：

将所有样本均值带入计算，提高同质样本的选取精度，进而得到分布式目标点。分布式目标（DS 点）在干涉体系中占据大比例，包括裸漏的岩石、低矮植被等，对矿区选取DS 点可以弥补永久散射体数量不足的缺点，避免失相干，进而提高干涉处理精度。

通过时序技术分别提取永久散射体点目标，将分布式点目标与永久散射体目标进行叠加（图4），采用Delaunay 三角网法对两种相干目标点统一构网。差分干涉图中,相邻目标点相位差表示为：

图4 干涉点目标分布图Fig.4 Image of interference point target distribution

式中：Δυ——相邻点之间的形变速率；

Δε——高程误差；

Δω——残余相位（包括大气延迟相位，非线性形变相位和噪声相位）；

B⊥——干涉对的空间基线；

T——干涉对的时间基线；

λ——波长；

R——传感器到目标的距离；

θ——雷达入射角。

2.3 干涉点目标形变处理

对目标点进行第一次回归分析时可以得到其高程误差、线性形变速率以及残余相位，由于第一次回归分析主要是为了估算初始的形变速率范围与数据处理效率，并没有顾及大气的与非线性形变相位的影响，因此在回归分析时，需要通过设置相位标准差小于1.5 rad剔除质量较低的点，通过两次迭代可以使高程误差的趋近于0，标准差显著降低[17]。在虑性形变的条件下，形变速率增量以及高程误差增量通过整体相位相干系数最大化模型来进行模型优化，公式为：

其中，γ越大，表示模型估计与差分相位越相似误差越小，以 γ作为权重，使用带权的最小二乘方法，从已知参考点解算出稀疏格网上目标点的时序线性沉降速率 υ和高程误差ε。将解算的每个目标点的线性形变和高程误差从初始差分干涉图中减去,即得到残余相位 Δω，其包含了大气扰动相位、非线性形变相位以及噪声相位。

再对点目标进行第二次回归分析，目的是使得非线性形变相位于大气延迟相位分离，因为大气延迟相位在1 km2范围内具有较强的空间相关性，在空间域上表现为平滑的低频信号，在时间域上呈现高频信号可被当作随机噪声处理[16−18]。因此根据其各自的特征进行合适的空间滤波，可将非线性形变相位和大气延迟相位分离，得到精确的非线性形变量。最后，根据计算的线性形变量和非线性形变量，得到形变相位进行相位并将其转到为形变量，视线向形变分解为垂直向形变。

3 结果与分析

对时序InSAR 结果进行区域裁剪，去除失相干区域，得到矿区及其周边城区的形变速率图与部分形变累积量如图5、图6，可得到在2019年全年中部露天矿区靠近城市侧的最大沉降速率达到569 mm/a，北部露天矿区远离城市侧最大沉降速率达到630 mm/a，中部露天矿附近城区地表最大沉降速率达到23 mm/a，北部露天矿附近城区地表最大沉降速率达到26 mm/a。由于露天矿区周边地质条件不稳定，露天矿旁道路年形变速率最大达到71 mm/a，城区的西侧沉降速率明显高于东侧。

图5 地表年形变速率图Fig.5 Annual surface deformation rate map

图6 形变累积量图Fig.6 Deformation cumulant diagram

3.1 地表形变及影响机理分析

矿区与城区之间地质条件不够稳定，最大沉降速率达到71 mm/a，平均沉降速率达到28 mm/a。试验区在1—3月平均沉降量5.3 mm、4—6月平均沉降量14.7 mm、7—9月平均沉降量16.2 mm、10—12月平均沉降量12.5 mm，其中4—6月和7—9月两个时间阶段沉降累积量较大。三个月为一组时间段，其平均沉降量、最大沉降量见表2。

表2 研究区累积沉降量表Table 2 Cumulative settlement scale

4—6月形变量分析，中部露天矿最大沉降量达到57 mm，北部露天矿最大形变量达到54 mm，城区地表最大沉降量为4 mm。通过对降雨量的分析，霍林郭勒市4—6月降雨量与往年同期相比减少约30%，土地解冻是矿区沉降量剧增的主要因素。中部矿区西侧的地表形变量为正值，是由于矿区西侧地质条件不稳定和排土场堆弃排土所致。7—9月形变量分析，中部露天矿最大沉降量达到61 mm，北部露天矿最大沉降量达到59 mm，城市地表沉降量最大值达到12 mm。由于霍林郭勒市7—9月的降雨量与往年同期相比增加约50%，雨水可能导致地表沉降量增大，北部露天矿南侧地表形变量为正值，由于北部矿区开采排土场积土导致。对比上述两时间段形变量，可以看出中部与北部露天矿的南侧形变量均从由负值变为正值，均为排土堆砌所致，可以从矿区开采得到验证。采用Kenda 系数法相关系数验证雨季时两者之间关系，经计算得系数为0.600，在夏季与秋季侵蚀性降雨与地表形变速率有正向的影响关系，可知在雨季降雨侵蚀力对露天矿坑周边地表环境变化产生较大影响，图7 为侵蚀性降雨量与监测点形变量关系图。

由于矿区形变复杂，将像元值数目低于总数5%的剔除，并将形变速率图与强度影像进行叠加得到图8，可以看出，城区地表形变速率与其距离露天矿的距离成反比，即城区离露天矿越近地表的形变速率值越大，从客观上也符合矿业城市地表沉降的规律，城市并未出现大规模的不规则沉降，矿区与城市之间道路沉降值较大均在32 mm/a 以上，靠近露天矿一侧的道路沉降值远大于另一侧，道路倾斜严重。

图8 基于强度影像的地表形变速率图Fig.8 Surface deformation rate map based on intensity image

3.2 GPS 数据验证与分析

为了验证改进得时序InSAR技术监测结果的精度，采用同期GPS 结果进行比对。在图7 上提取6 个点并与同样的GPS 点监测结果进行比较（表3），最小互差0.76，说明该InSAR 监测方法是可行与可靠的。

图7 侵蚀性降雨量与监测点形变速率图Fig.7 Chart of erosive rainfall and deformation rate at monitoring points

表3 GPS 监测结果与同期InSAR 监测结果比较Table 3 Comparison of GPS monitoring results with InSAR monitoring results in the same period

4 结论

（1）考虑到露天矿区地表的低相干性，提出了一种改进的时序InSAR技术，采用24×24 窗口同质样本选取方法提取分布式目标，最大化的识别矿区的相干点目标，对于霍林郭勒市城区采用经典的永久散射体干涉测量方法提取永久散射体，对两种干涉点目标进行统一的构网分析和形变反演，与GPS 数据进行对比验证，证明改进的时序InSAR技术方法具有较好的监测效果，对矿业城市地表形变监测提出了一种新的InSAR 监测方法。

（2）将霍林郭勒市的月降雨量数据和干涉点月沉降量关联分析，Kenda 系数达到0.736 证明二者具有较高相关性，可知侵蚀性降雨对于露天矿区及其周边地表沉降具有正向的影响。

（3）通过本方法得到露天矿区地表沉降速率最大可达346 mm/a，城区地表形变速率最大可达31 mm/a，露天矿和城区之间道路沉降量巨大，且公路两侧形变速率值差异巨大，道路出现明显的倾斜。剔除形变速率极值，可知城区地表形变速率从客观上符合矿业城市地表沉降的规律，即城区形变速率值与其到露天矿的距离成反比，可得城区并未出现大规模的不规则沉降，矿区与城市之间道路沉降值较大均在32 mm/a 以上，靠近露天矿一侧的道路沉降值远大于另一侧，道路倾斜严重。

感谢欧洲空间局( ESA)为本文提供的SAR 数据和国家科学数据服务平台提供SRTM-Version4 数据。