基于MT-InSAR和灰色预测模型的金属矿区地表变形研究*

王 华 ，李素敏,2,3，袁利伟，张 玮，马相松，管庆丹

(1.昆明理工大学 国土资源工程学院，云南 昆明650093；2.云南省高校高原山区空间信息测绘技术应用工程研究中心，云南 昆明650093；3.中国有色金属工业协会智慧矿山地理空间信息集成创新重点实验室，云南 昆明650093；4．昆明理工大学 公共安全与应急管理学院，云南 昆明650093；5.玉溪大红山矿业有限公司，云南 玉溪 653400；6.昆明有色冶金设计研究院股份公司，云南 昆明 650000)

0 引言

金属矿山由于高强度和大面积的地下开采，导致岩层结构遭到破坏,地表沉陷现象频发，易对矿区基础设施和周边生态环境造成不良影响[1]。掌握金属矿山重大风险灾害总体空间分布特征，明确其发展动态，并进行遥感监测和模型预测分析，对矿山安全监管与预测预警研究具有重要的现实意义。因此，近年来学者针对矿区地表变形监测和预测模型做了大量研究[2-3]。

矿区变形监测传统方法有水准测量、GPS测量等，虽然这些方法精度高，但是由于技术的局限性，仅能监测地面“离散点”，很难反映地表大面积的“面域”变形信息和时空演化特征，且需耗费大量的人力、物力和财力[4]。近些年，在微波遥感领域，合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar, InSAR)已成为对地观测的研究热点之一，InSAR技术具有全天时、全天候的工作能力和高精度、高分辨率、覆盖范围广等特点，并能对地连续观测获取地表“面域”的变形信息[5]。时间序列InSAR技术(Multiple Temporal InSAR, MT-InSAR)是在D-InSAR的基础上发展而来的。小基线集SABS技术最先由BERARDINO等[6]在2002年提出。李达等[7]利用SBAS和D-InSAR两种方法获取了某矿区的累计沉降值，并对其进行差值分析和量化分析，为矿区地表的变形监测与分析提供了新手段。目前对于矿区地表的监测研究成果已有很多，如：胡海斌等[8]运用灰色模型对沉陷区的沉降量进行了预测，经实例验证,灰色预测模型能满足工程预测精度要求；李路等[9]加入权值思想获得了改进后的灰色系统模型，并将其运用于实际沉降预测中；李金超等[10]提出了一种基于灰色支持向量机的组合预测模型并对矿区变形进行了预测，结果表明，该模型可实现对矿区变形的快速监测。

如何利用现有的贫信息、小样本监测数据[11]对矿区开展沉降预测是目前亟待解决的问题。本文以大红山铁矿为研究对象，利用MT-InSAR技术获取矿区时间序列沉降量及变形速率，根据矿区地质结构和降雨量分析矿区地表变形机理，结合地面水准监测数据验证MT-InSAR技术在矿区变形监测中的可行性和可靠性，并引入灰色沉降预测模型，以实现对小样本数据等劣势条件下的变形趋势短期预测，为矿区安全生产提供保障，预防重大灾害的发生。

1 研究区概况及数据来源

1.1 研究区概况

大红山铁矿是我国主要的铁矿石生产基地之一，位于云南省玉溪市新平县城西119 km的新化乡尺莫村，研究区地理位置与地形如图1所示。矿区位于滇中台坳南端以及红河断裂与绿汁江断裂夹持的三角地带，矿区西侧出露有变质较深、混合岩化较强的太古代哀牢山群，该区为北西、东西和南北向3组构造线交会地带。矿区属亚热带气候，干湿季分明，雨量较为充沛，年降雨量700～1 200 mm，多年平均降雨量930.8 mm，60%以上的雨量集中在6-9月。矿区目前有排土场、露天开采和地下开采区域，其中露天采场基本工作面高程在800～1 000 m，地下开采强度大且开采模式复杂，地表裂缝和滑坡不断出现，传统测量方法实施难度大、成本高。

图1 研究区地理位置与地形

1.2 数据来源

采用欧洲航天局拥有12 d重访周期的C波段传感器的Sentinel-1A卫星数据，选择干涉宽幅模式下的斜距单视复数影像，极化方式为垂直极化。采用30 m分辨率的SRTM DEM描述地形。影像数据参数见表1。

表1 Sentinel-1A影像数据参数

2 研究方法及原理

2.1 SBAS-InSAR数据处理

InSAR技术可以提供基于时间序列雷达数据的开采沉陷信息和地表变形的时空演化信息[12]。SBAS-InSAR技术是通过获取N景雷达数据设置短时间基线和空间基线阈值(本文为80%和10%)再将其组合，用Deluanary MCF最小费用流解缠方法得到662对短基线差分干涉图，根据相干系数选择稳定的高相干点目标，利用奇异值分解法(SVD)[13]解算估计线性变形速率和高程误差方程，再利用残余相位估计大气效应等相位，最终通过地理编码将相位转变为变形从而获得相干点的变形速率及累计时间序列沉降。本文通过如图2所示的SBAS-InSAR数据处理流程依次去除平地相位、地形相位、轨道误差、大气相位等，最终提取地面变形相位。通过获取累计InSAR时间序列变形并将其作为灰色模型的输入值，以实现对矿区变形趋势的准确预测。

图2 数据处理流程

2.2 预测模型原理

灰色预测模型是通过对原始数据进行累加处理生成一组灰色序列，减小其原始序列的随机性，以此凸显其指数增长规律，然后建立一阶灰色微分方程，利用最小二乘法原理求取灰参数，最后经累减还原得到预测值，其原理[14]叙述如下。

1)数据的检验和处理

原始数据序列：

x(0)=[x(0)(1),x(0)(2),…,x(0)(n)]。

(1)

2)建立灰色模型

生成x(0)的累计序列：

x(1)=[x(1)(1),x(1)(2),…,x(1)(n)]

=[x(0)(1),x(0)(1)+x(0)(2),…,

x(0)(1)+…+x(0)(n)] 。

(2)

生成x(1)的均值序列：

z(1)(k)=0.5x(1)(k)+0.5x(1)(k-1),k=2,3,…,n。

(3)

建立白化微分方程：

(4)

建立灰色微分方程：

x(0)(k)+az(1)(k)=b。

(5)

最小二乘法计算发展系数a和灰作用量b，得到预测方程：

(6)

3 结果分析

3.1 精度验证

为验证MT-InSAR技术的监测精度，将水准点数据与InSAR数据作对比分析。InSAR测量结果与真实三维空间位移的函数关系[15]为

DLOS=DV·cosθ+DN·sinφsinθ-DE·cosφsinθ，

(7)

式中，DLOS表示沿卫星雷达视线方向的位移；DV表示地面点沿垂直方向的位移；DN表示地面点沿南北方向的位移；DE表示地面点沿东西方向的位移；θ表示雷达卫星入射角度；φ表示雷达卫星飞行的方位角。

由于地表的水平位移不能忽略，且只有降轨数据，所以只能将获取的12个地表观测点实测水准数据三向位移(2020年4月8日至2020年10月8日)代入式(7)获得归一化的D′，并与由MT-InSAR观测得到的DLOS进行对比，以验证InSAR监测的可靠性，结果见表2、表3。

根据地面观测点的实测三向数据(北向、东向、垂向)经归一化后得到的D′与InSAR取得的DLOS比较，求得绝对误差和相对误差(见表3)。由表3可知，两者的绝对误差值均小于2 mm，相对误差在15%以内。这表明MT-InSAR技术不仅能获取地面变形情况，且其观测精度比传统测量技术的高。

表2 地面观测点实测三向位移 单位：mm

表3 归一化后的D′与DLOS的误差

3.2 InSAR结果分析

通过获取研究区138景(2014年10月至2020年10月)降轨数据，利用SBAS-InSAR技术获得研究区的变形场，从而掌握矿区的变形空间分布和变形特征。研究发现矿区出现了3个明显的沉降中心，即A、B、C区域，沉降面积分别为0.37、0.20、0.10 km2。其中，A沉降区域已经发育成较成熟的沉降集中带。另外，矿区还存在局部沉降区域D、E、F(见图3)。A、B、E沉降区域主要与堆排土体有关，C、D沉降区域主要是地下开采所致，F沉降区域主要受露天开采影响。为进一步探讨矿区的变形机理和变形特征，分别在6个沉降区域选取9个特征点，并截取A、B、C沉降区域中心的剖面图。

p-特征取样点；JC-水准点。

A、B区域是矿区的南部排土场，此处已经发育形成2个典型的沉降集中带，累计变形量在-555～-70 mm，平均变形速率在-90～-12 mm/a。对A区域经过两个沉降最大中心处沿a-a′的剖面(见图4b)分析可知，这两个沉降中心呈西北至东南方向的走势，且该排土场的海拔在800～1 200 m，落差较大，沉降边缘不断向两个方向扩散，且存在明显的不均匀沉降。

图4 A、B区域特征点时序曲线及a-a′剖面图

获取p1、p2、p3、p4四个点的时间序列变形，研究发现4个取样点呈近乎一致的沉降规律，沉降曲线逐渐减缓，表明变形速率不断减小；在每年的6-8月雨季期间，沉降曲线会发生小幅波动，这是因为该排土场虽然不存在地下开采，但堆排土质疏松(见图5)，且土壤中存在大量积水，弱化了排土场的稳定性，从而导致沉降的发生。

图5 排土场现状

C、D区域位于矿区地下开采的上覆地表，研究区变形场和变形速率图显示，该区域地表也形成了不同程度的沉陷，这两处区域累计变形量在-348～-70 mm，平均变形速率在-55～12 mm/a。该区域长期存在地下开采，开采深度达1 200 m，规模为400万t/a，地下高强度和大范围的采动引起地表大面积沉降。分析C、D区域的p5、p6点的时间序列沉降(见图6a)发现，雨季(6、7、8月)前后，两处的变形曲线有明显的波动，尤其是在2016年和2018年，且有加速下沉的趋势，这可能与该时期的降雨量较大有关。分析沿两个下沉中心的d-e剖面图(见图6b)发现，C、D区域依旧存在不均匀沉降，表明降雨也是造成地面塌陷的原因，且降雨量与地面下沉速率呈正相关，但此变形主要由地下开采所致。

图6 p5、p6点时间序列沉降曲线及d-e剖面图

由E区域的p7点的变形时间序列曲线(见图7a)可知，该区域的变形速率逐渐减小，与实际情况相符。F区域是露天采场，边坡上部坡度较大，由MT-InSAR获取的变形(见图7b)情况来看，下沉量较小，曲线呈周期性波动，这与露天边坡长期的矿体开采有关。

图7 p7、p和s点时间序列沉降曲线

3.3 预测分析

选取矿区水准点作为建模数据，以CD01、JC04两个实测点为研究对象，分别获取12期的水准监测数据，提取同时间段SBAS-InSAR值，将其作为输入值获取沉降预测方程，绘制预测曲线，通过已获取的小样本数据预测矿区变形趋势。将CD01、JC04点的第10、11、12期的预测值分别与InSAR值和水准数据对比，以验证模型的预测精度和可靠性，结果见图8、表4。

图8 CD01、JC04点的沉降曲线对比

表4 CD01、JC04点的沉降量相对误差

由图8、表4可知：监测点CD01、JC04的灰色模型预测曲线与水准及InSAR曲线较为吻合，模型拟合效果较好；CD01、JC04点在第10、第11期的监测相对误差在5%以内、预测相对误差在9%以内，而第12期的监测相对误差在3%以内、预测相对误差在12%以内。但随着监测周期的拉长，CD01、JC04点的监测相对误差和预测相对误差会增大，这与文献[16]的结论一致。该灰色模型计算简单快捷。

4 结论

a.利用时序SBAS-InSAR技术获取研究区的变形场，研究发现矿区内存在4个发育成熟的沉降集中带A、B、C、D及局部沉降区域E、F。区域A和B主要是由排土场局部的不稳定和雨季降水综合影响所致。区域C和D主要是由长期频繁的地下开采及雨季强降雨综合影响所致。区域E是废弃的排土场，其变形情况与实际相符。区域F是由南向北开采的露天采场，沉降量较小，曲线呈周期性波动，这与露天边坡长期的矿体开采有关。

b.比较地面观测点实测三向位移经归一化后得到的位移与InSAR取得的LOS向位移后发现，绝大部分点的相对误差在15%以内，表明MT-InSAR在观测精度上是可靠的。

c.灰色模型获得的预测曲线与水准及InSAR曲线较为吻合，预测模型的相对误差在12%以内。