基于D-InSAR技术和改进GM（1，1）模型的矿区沉降监测与预计

2020-10-26石晓宇魏祥平杨可明姚树一

金属矿山 2020年9期

石晓宇魏祥平杨可明王剑姚树一

（1.中国矿业大学（北京）地球科学与测绘工程学院，北京100083；2.淮北矿业股份有限公司，安徽淮北235000）

矿产资源大规模开采所引起的地面沉降会危及矿区及周边范围内建构（筑）物安全，甚至会影响矿区生态环境与居民生活安全，因此对矿区地表沉降进行监测预警很有必要。相较于传统监测手段，合成孔径雷达差分干涉测量（Differential Interferometry Synthetic Aperture Radar，D-InSAR）技术具有全天候全天时、高精度、大范围连续空间覆盖、成本低、零风险等优势。该技术通过雷达干涉图差分获得视线向微小形变［1-2］，被广泛应用于火山运动［3］、山体滑坡［4］、地震形变［5-6］、采矿区塌陷［7-8］以及城市地面沉降［9］等环境与灾害监测中。就开采沉陷预计而言，由开采引起的地表沉陷量会随着时空的变化而变化［10］，是一个关于多元变化的因变值。常规预计方法主要为概率积分法，但该方法无法较好地反映开采沉陷随时间变化的动态过程，在非充分采动或极不充分采动时，预计结果往往适用性不强［11-13］。灰色模型（Grey Model，GM（1，1））以部分信息已知、部分信息未知的小数据和贫信息不确定系统为研究对象［14］，通过对已有的时间序列数据进行累加、累减等变换处理，生成具有变化规律的新时间序列，可以弱化原始时间变量的随机性。因而可以通过建立GM（1，1）方程预测未来某一时刻的特征量，即可以有效演化地表沉降量随时间的变化规律。已有研究表明：矿区变形监测数据具有一定的灰色特性［15-17］，然而该类研究往往是基于观测站的水准或GPS测量结果，且与后续的演变预计环节彼此独立，不能形成一个完整且系统化的监测与预测模型。本研究联合D-In-SAR技术与GM（1，1）模型，建立描述沉降量与时间关系的模型方程，实现开采沉陷监测与预计的一体化，为实现矿区地面沉陷的动态监测与实时预测预警以及后续治理提供参考。

1 理论与方法

1.1 D-InSAR基本原理

D-InSAR是利用同一地区不同时相SAR影像间的差分干涉，获取该地区的地表形变信息的一种监测技术［18］。其主要思路是借助外部DEM数据以及利用形变前后两景SAR影像间雷达成像具有几何微小差异的干涉信号，获得干涉相位。通过去除其余相位的干扰影响，获取视线（LOS）方向的地表形变信息。两次雷达卫星成像记录的相位差值φ可近似表示为

式中，φflat为无高程变化时平坦地形引起的相位变化；φtopography为应用干涉或外部DEM数据可消除的地形相位；φdeformation为在LOS方向上地面形变引起的相位变化；φatmosphere为大气延迟引起的相位偏移；φnoise为噪声相位。

由式（1），经相位解缠提取地表LOS向形变相位φdeformation'，并将其转化为SAR卫星沿视线方向的形变值dLOS，即：

式中，k为相位整周模糊度；λ为雷达波波长。

由于垂直向形变是井下开采引起地表形变的主要贡献分量，故忽略水平向形变对雷达视线方向形变的影响，经过解三角形可得到形变前后SAR影像间的地表垂直形变值w为

1.2 经典GM（1，1）模型理论

式中，a,b为待估参数。

求解式（7）可得：

1.3 改进的GM（1，1）模型

本研究提出了2种改进的GM（1，1）模型，分别为基于补偿最小二乘法估计半参数GM（1，1）模型（BGM（1，1））和赋予相对权重的补偿最小二乘法估计半参数GM（1，1）模型（WGM（1，1））。

1.3.1 BGM（1，1）模型

BGM（1，1）相较于经典GM（1，1）模型，在方程中加入了非参数分量S，较好地顾及了模型中可能存在的系统误差，则模型方程为

式（9）中加入的S也是未知量，使得BGM（1，1）变为秩亏模型，若按最小二乘方法求解，无法得到唯一解，所以必须添加约束条件即补偿最小二乘估计准则，从而得到唯一解，即求解方法由VTPV=min转变为

根据螺旋CT成像分析和颈动脉内膜狭窄测量标准来确定患者的狭窄或闭塞程度:狭窄在30%之内为轻度狭窄,范围在31~69%之间为中度狭窄,70%~99%为重度狭窄,100%为闭塞。

式中，α为平滑参数［20］；R为正则矩阵。

1.3.2 WGM（1，1）模型

WGM（1，1）模型是在BGM（1，1）模型的基础上，基于平滑参数α求取方法对补偿最小二乘法进行相应改进，求解方法为

WGM（1，1）模型相较于BGM（1，1）模型，将平滑参数α的范围由[0 ,+∞]区间限定在[0 ,1]区间，在很大程度上，使平滑参数α的搜索范围减小，从而节省了搜索时间。

2 D-InSAR监测分析

2.1 研究区概况及数据源

淮北矿区袁二煤矿位于安徽省涡阳县境内，采区地面纵横分布有较多人工沟渠和村庄等建筑。研究区7221工作面为俯斜开采，于2017年12月6日开始回采，2018年6月3日停采。当工作面进行顶板完全垮落式回采时，工作面上方地表将形成塌陷区并对地面公共设施和村庄等建筑造成一定程度破坏。为了精确监测地表形变，分别在地表公路或建筑体上设置了一些水准测量与变形监测点，监测地表及其重要建筑的变形情况。由于该类监测点的数量与监测范围有限，且常常被人为破坏，无法对重要建筑物形变进行精确监测，进而无法对该类建筑物的损害情况进行有效评估。

研究区影像采用Sentinel-1A数据，外部数字高程模型选用分辨率为90 m的SRTM（Shuttle Radar Topography Mission，航天飞机雷达地形测绘使命）数据。从欧空局网站（https://scihub.copernicus.eu/）下载了2017年11月16日—2018年1月27日期间共6幅SLC数据，形成5个干涉对，干涉对影像的详细信息如表1所示，其中VV为同向垂直极化。

2.2 D-InSAR处理方法与监测结果

差分干涉常有“两过差分”、“三过差分”和“四过差分”［21］3种方法。本研究采用较为成熟的“两过差分”法，利用监测时间段内相邻时间获取的2幅SAR影像形成干涉对和对应的外部DEM数据，处理得到仅包含地表形变信息的干涉图，并对得到的5期相继时段内地表沉降量进行时序叠加处理，最终得到7221工作面2017年11月16日—2018年1月27日地表沉降量的时间序列形变图，如图1所示。图1（a）至图1（e）是以2017年11月16日为参考时间得出的地表形变，可以看出，在回采前期以及回采初期（图1（a）和图1（b）），地表下沉量及影响范围均较小，随着开采时间推移，开始形成下沉盆地（图1（c）至图1（e）），自西南向东北方向扩张（与7221工作面推进方向一致），平面范围逐渐增大。至2018年1月27日，最大下沉量为71 mm，位于7221工作面下方区域，7221工作面上方村庄范围内未发生明显沉降，说明所采用的注浆充填开采工艺（注浆孔位于工作面东北部）有效发挥了局域减沉减变作用。

2.3 D-InSAR监测成果验证与分析

为了验证D-InSAR技术监测成果的可靠性，从D-InSAR面状监测结果中提取了各水准点的下沉值，与已有水准监测数据进行对比。由于水准数据监测时段为2017年10月8日—2018年1月30日，而D-InSAR研究时段为2017年11月16日—2018年1月27日，两者时间基准不一致。为了尽可能减少因时间间隔不统一而带来的误差，以地表下沉速率为依据将水准数据起止日期内插至与研究时段相同的日期，即为2017年11月16日—2018年1月27日。采区7221工作面水准点布设位置如图2所示，水准监测成果与D-InSAR监测结果对比如图3所示。

结合图2和图3可知：D-InSAR监测结果与常规水准测量成果相比，下沉趋势基本一致，D-InSAR监测的工作面各观测线上分别在H3和I2点处（离工作面较近）下沉量最大；虽二者下沉趋势一致，但就下沉量而言，存在一定差异，在I观测线上表现得最为明显，H和C观测线上沉降趋势则大致相同；在I观测线上，I2、I3、I4点的水准测量下沉值与D-InSAR监测下沉值相差较大，下沉值之差最大约40 mm，推测可能是由于水准测量中人为误差传递或特别时间段的失相干导致；在F观测线上，F1点的水准数据与D-In-SAR监测数据相差最大，差值为25 mm。

3 改进的模型监测与验证分析

3.1 改进的GM（1，1）模型拟合试验

由D-InSAR技术获取的7221工作面开采阶段性地表沉降量，对该工作面上监测精度较高的H观测线上的4个观测点（H1、H2、H3、H4点）数据进行后续的建模与预测。从各形变图中提取2017年12月10日—2018年1月27日期间各观测点对应的像元值作为原始时间序列值，分别采用GM（1，1）、BGM（1，1）和WGM（1，1）模型进行拟合，拟合结果采用残差绝对值均值、相对精度、后验差等指标进行精度分析，3个模型与D-InSAR反演值的拟合精度分析见表2。

表2中，经典GM（1，1）模型加入系统误差之后，各点的拟合精度均有一定程度提高，就残差绝对值均值而言，减少了0.074～1.588 mm，平均相对误差以及相对精度最大分别减少和提高了2.51%和1.69%。对3个模型拟合结果进行综合分析，发现BGM（1，1）与WGM（1，1）模型均优于经典GM（1，1）模型；BGM（1，1）模型与WGM（1，1）模型虽然在拟合精度上相差很小，但在平滑因子α的选取方法上，WGM（1，1）模型确定了搜索上界k（k=0.5、0.8、1、5、10，步长均为0.01），因而节省了搜索时间。该模型在求解平滑参数α的过程中，最大搜索上界为1，因此可直接设置k等于1。两者模型平滑参数求取结果对比如表3所示。

注：残差绝对值=|实测监测值-模型拟合值|；相对误差=残差绝对值与实测值的百分比；关联度表示拟合曲线与原始曲线的相似程度；相对精度=1-平均相对误差；后验差为残差方差与原始数据方差之比。

根据上述分析，本研究采用WGM（1，1）模型进行沉降预测。

3.2 WGM（1，1）模型预计试验

采用WGM（1，1）模型预测4个像元第7期（2018年2月20日）沉降数据，并采用同期的水准监测成果进行对比验证，第7期预测值与水准数据对比见表4。由表4可知：H3和H4点第7期预测值与水准监测值基本吻合，相对误差在4%以下，精度较高，而H1和H2点预测相对误差最大达26.64%。原因为：利用D-In-SAR技术获取沉降预测样本的同时，也将沉降误差带入模型；时间序列周期较短，仅为12 d，监测点下沉规律复杂，6期数据中包含了沉降3个阶段（初始期、活跃期、衰退期），使得原序列随机性增强（如监测点H1），所以原序列经过模型处理后，仍然存在较大的随机性；H2点的原序列值虽然随机性较弱，但前4期沉降速度较快，形变至第5，6、7期阶段，沉降反而变为由缓慢至稳定阶段，原始序列值趋势呈“S”型曲线状态，而灰色模型适用于呈现指数增加的下凸曲线，致使模型方程不能准确地进行沉降预测。