融合D-InSAR与GIS技术的矿区开采沉陷形变监测及预测方法

陈银翠，徐良骥，余礼仁

(1. 安徽矿业职业技术学院，安徽 淮北 235000； 2. 安徽理工大学，安徽 淮南 232001)

1 数据处理方法和流程

1.1 D-InSAR二轨法数据处理流程

二轨法利用形变前后的两景SAR影像差分干涉处理再导入外部DEM可得到地表形变相位[7-9]；干涉数据处理是基于ENVI SARscape平台，其数据处理流程如图1所示。

1.2 GIS软件中数据处理方法及流程

D-InSAR获得的栅格结果和GIS中矢量数据进行叠加分析，在ArcMAP中加载出D-InSAR栅格格式形变图与矿区开采工作面、矿区界限矢量图并显示在一个界面范围。在完成空间坐标系的转换和D-InSAR栅格形变图与GIS矢量数据融合后，对融合数据可进行沉降区的提取、剖面分析和叠加分析[11]。

2 基于D-InSAR监测值的三次样条插值与GM(1,1)组合预测模型的建立

2.1 灰色理论预测模型

经典的GM(1,1)模型为等间距序列预测模型[12-13]。设X0为大于零的原始数据序列，将X0进行一次累加得到AGO序列X1，Z1是X1的邻均值等权序列。即

X(0)={X(0)(1),X(0)(2),X(0)(3),…,X(0)(k)}

(1)

X(1)={X(1)(1),X(1)(2),X(1)(3),…,X(1)(k)}

(2)

其中

X(1)(n)=X(0)(1)+X(0)(2)+X(0)(3)+…+X(0)(n)

(3)

则GM(1,1)模型相应的白化微分方程为

酒店认知实践是酒店管理专业大一学生的入学教育。课程通过学生参观酒店、访谈酒店管理层、酒店经理人讲座等具体形式来开展实践课程。学校安排学生参观优质校企合作单位，使学生设身处地地接触实践教学基地的真实对客服务场景，让学生主动调研高端酒店的基本情况，增加同学对酒店行业的认同感。

(4)

式中，a为发展系数；b为内生控制灰数。

若令c=[a,b]T,那么c可以表示为：c=(BTB)-1BTYn

其中

t=1,2,…,n

(5)

Yn=[x(0)(2),x(0)(3),…,x(0)(n)]T

(6)

白化微分方程的解为

(7)

式(7)求出的是一次累加后的生成值，对其进行累减还原，就能够得到预测值，即

(8)

2.2 结合三次样条插值的非等间距GM(1,1)模型

在开采沉陷监测中，相邻两次监测数据时间间隔一般不相等，因而不能使用经典GM(1,1)模型对矿区沉降结果进行有效的预测分析。为了克服单一模型算法所带来的局限性，提出将三次样条插值与经典GM(1,1)模型相结合的模型算法[14]。模型算法的具体步骤为：

(1) 将原始影像形变图中获取的非等时间间隔沉降数据进行三次样条插值。

(2) 由三次样条插值算法获得每个雷达重访周期内的下沉值，然后处理生成等时间间距序列，计算级比，确定所有样本值级比都落在区间[exp(-2/(n+1))，exp(2/(n+1))]内；如果不在，则需对数据做平滑处理。

(3) 使用经典GM(1,1)模型对处理后的数据进行计算得到白化微分方程。

(4) 求解白化微分方程得到预测序列并进行累减，对预测数据进行残差检验。

2.3 精度评定

表1 后验方差模型精度分级标准

3 D-InSAR融合GIS技术数据处理与结果分析

3.1 研究区概况及数据选择

袁店二矿位于安徽省亳州市内，该地区地势平坦，地面平均标高为30 m。7225工作面长630～637 m，宽144～145 m，标高为-411～-446 m。煤层平均厚度为3.81 m，煤层倾角为2～17.5°。由于工作面北部有村庄，采用注浆减沉方式进行回采。本文采用2景C波段Radarsat-2和6景Sentinel-1卫星影像数据对研究区进行分析，采用90 m分辨率的SRTM DEM数据。6个干涉对参数见表2。

表2 Radarsat-2和Sentinel-1数据干涉对参数

3.2 D-InSAR融合GIS技术干涉差分处理结果分析

采用二轨法对沉陷区范围影像数据进行差分干涉处理获取工作面形变图。将6组形变图和地表移动观测站图导入ArcGIS进行叠加分析，得到6个工作面的时序形变图如图2所示。

从图2(a)中可知，随着7225工作面开采工作的推进，在工作面的中间部分出现了下沉盆地。由图2(b)—(f)可知，从2015年7月至2016年1月，7225工作面上方每一个干涉对时间段内都形成了下沉区域，而且下沉较大的区域随着地下煤炭的开采逐渐上前推进。

4 基于三次样条插值与GM(1,1)组合模型的矿区沉降预测分析

选取2015年4月19日至12月21日6个监测点的沉降数据为研究对象，分别为倾向线和走向线工作面上方的A14、A15、A16和B10、B11、B12，见表3。将2015年4月1日间隔24 d的下沉值作为首次监测数据，通过三次样条插值获得2015年4月1日至12月21日每个运行周期内下沉值。将每个周期内下沉值逐次累加得到相对于首幅影像日期的累计下沉值见表4。

将此等间距数据作为GM(1,1)模型的原始数据前9期数据(2015-04-25—2015-11-03)建模预测后两期(2015-11-27和2015-12-21)累计下沉量，并与原始等间序列的后两期数据进行对比，验证预测模型的精度，结果见表5。由表5可知：到2015年11月27日的累计下沉值最大相对误差出现在B10点处，其值为11.5%，其余都在10%之内；到2015年12月21日的累计下沉值最大相对误差出现在B10点和B12点其值为11.2%。最大方差比为0.097，预测模型的精度能够满足D-InSAR定量分析的精度要求，因此可以进行沉陷预计分析。

m

表4 等间距下沉序列 m

表5 6个监测点预测值与实测值对比分析

由于7225工作面采用注浆减沉开采，因此预计地表的累计下沉量对安全开采有着重要的作用。由概率积分经验公式可知[15]，注浆开采的下沉系数一般在0.06～0.2之间。结合7225工作面地下采煤地质条件及工作面回采情况可以计算出地表充分采动时的最大下沉量。根据下沉系数公式得到地表最大下沉值在227～756 mm之间。因此根据位于中心盆地的B12点下沉值430 mm可知，7225工作面注浆减沉效果明显，在2015年5月到下沉稳定期间一般不会发生重大的塌陷事故。

5 结 论

(1) 本文采用融合D-InSAR和GIS技术对淮北矿业集团袁店二矿7225工作面进行分析，获得6个工作面的时序形变图。结果表明工作面下沉区域随着工作面开采工作的推进逐渐上前推进，秋冬季节与夏秋季节相比相干性较好。

(2) 利用三次样条插值与GM(1,1)组合模型，以7225工作面上6个监测点的沉降数据为研究对象进行预测。D-InSAR预测值与实测值最大相对误差和方差比分别为11.5%和0.097，该精度能够满足D-InSAR定量分析的精度要求；通过下沉系数计算最大下沉值430 mm处在227～756 mm之间，一般不会发生重大的塌陷事故。