0 引言
煤炭是目前我国主要的供给能源[1]。地下煤炭资源的大面积开采使得采区周围岩层原始应力平衡被打破[2]，引起岩层移动和变形，不仅直接造成矿区地表沉陷，还会带来一系列潜在的灾害问题。对矿区地表沉陷规律的研究是采矿工程中一个重大课题，其中如何更高效、更系统、更安全地对矿区地表进行沉陷监测与预计备受关注。
地下采煤引起的岩层移动和变形是一个复杂的时空过程，传统的地表沉陷预计方法(如概率积分法、剖面函数法等)仅能获得地表沉陷最终稳定后的空间变化，不能较好地反映开采沉陷随时间变化的动态过程[3]。众多学者综合考虑因岩体本身复杂性而难以全面研究采煤工作面上覆岩层移动及地表沉陷动态过程的情况，以局部观测点为切入点，建立沉陷区单点沉陷量时间序列模型，从而研究矿区地表动态沉陷规律[4-6]。GM(1,1)与灰色Verhulst模型为灰色系统理论中2个基础的时间序列模型，可在小样本、贫信息、不确定性等劣势情况下，鉴别系统因素间发展趋势的相异程度，寻找系统变化规律，并对事物发展趋势作出模糊性的短期、中长期预测[7-8]。GM(1,1)是用一阶微分方程对1个变量建立的模型[9]，主要解决具有指数变化规律的生成序列，只能描述单调变化过程；灰色Verhulst模型是由德国生物学家Verhulst在研究物种繁殖趋势时提出的一阶非线性动态模型[10]，其基本思想是生物繁殖数量受到周围环境的影响，繁殖曲线会以“S”形呈现，最终达到一个平衡稳定的状态。将这2种灰色模型(Grey Model，GM)应用于矿区地表沉陷预计，建立沉陷量时间序列模型，可较好地揭示地表单点连续动态变化的过程，实现地表沉陷动态预计。然而，现有文献只针对2种GM在矿区沉陷预计方面的应用进行研究，并没有结合模型自身特点，分析和比较二者在矿区地表沉陷预计中的适用性。此外，目前矿区地表单点动态沉陷预计方法主要基于传统的水准测量数据[11]，监测方法单一，成本高，且观测点易破坏，不能保证地表形变信息的实时性。合成孔径雷达差分干涉测量(Differential Interferometry Synthetic Aperture Radar，D-InSAR)技术通过雷达干涉图的差分获得地表形变信息，以高精度、全天候、连续空间覆盖、实时性等优势成为新兴的形变监测手段，被广泛应用于矿区地表沉陷监测中[12-15]。将D-InSAR技术与GM相结合，避免监测与预计独立进行，可较好地形成系统化及模块化流程，实现矿区地表沉陷监测与动态预计一体化。
本文以淮北矿业(集团)股份有限公司袁店二矿7221工作面为试验区域，融合D-InSAR技术和GM(1,1)与灰色Verhulst模型，通过建立沉陷量与时间的GM进行地表沉陷动态预计，实现监测与预计一体化，并针对2种GM的特点，结合预计结果分析和比较其在矿区地表沉陷预计中的适用性。
1 矿区地表沉陷D-InSAR监测
1.1 试验区域与相关数据
袁店二矿7221工作面上方地表主要覆盖物为村庄和农田，四邻无其他采掘活动。该工作面为俯斜开采，工作面宽约110 m，走向长约544 m，煤层倾角为10～20°，基岩厚度为132～162 m，采高为4.5 m，于2017-12-06开始回采，2018-06-03停采。为了控制地表采动损害在建筑物等安全保护等级范围之内，分别在地表公路或建筑体上设置水准测量与变形监测点，以监测地表及重要建筑物变形情况。但该区域内人工沟渠纵横分布，监测点稀少、空间配置不合理且常被人为破坏，无法准确对重要建筑物进行损害评估。因此，有必要在工作面开采过程中对该区域进行地表沉陷D-InSAR监测。
为了尽可能减少D-InSAR技术相位失相干对地表沉陷监测结果精度的影响，选取冬季和初春时段(该时段植被萧疏且地物波谱特征无明显变化)的9景哨兵卫星(Sentinel-1A)雷达影像作为D-InSAR监测的数据源。影像时间为2017-11-04—2018-02-08，相邻时间间隔(12 d)的影像可以形成8个干涉对，见表1，其中VV为同向垂直极化方式，IW为干涉测量宽幅模式。外部数字高程模型选用分辨率为90 m的SRTM(Shuttle Radar Topography Mission，航天飞机雷达地形测绘使命)数据。
表1 Sentinel-1A雷达影像干涉对
Table 1 Interference pairs of Sentinel-1A radar
序号主影像辅影像极化方式工作模式垂直基线/m时间基线/d12017-11-042017-11-16VVIW-21.9281222017-11-162017-11-28VVIW-69.6481232017-11-282017-12-10VVIW42.0991242017-12-102017-12-22VVIW92.7101252017-12-222018-01-03VVIW17.2421262018-01-032018-01-15VVIW-60.4071272018-01-152018-01-27VVIW-38.8091282018-01-272018-02-08VVIW-62.15612
1.2 D-InSAR监测结果及分析
根据表1数据，利用D-InSAR技术二轨法获取8组7221工作面时间序列沉陷分布。为方便分析，以2017-11-04雷达影像为参考对象，将各组沉陷量依次累加，即可得到各组时间序列沉陷分布。由于2017-11-04—28地形无明显变化，所以仅给出后6组时间序列沉陷分布，如图1所示。可看出在回采初期(2017-12-06—10)，沉陷量为-10～0 mm。随着工作面推进，进入沉陷发展阶段(图1(b)—图1(f))，在2017-12-22开始形成下沉盆地，与7221工作面推进方向一致。下沉盆地自西南向东北方向扩张，沉陷量和平面范围随之增大。至2018-02-08，最大沉陷量为-85 mm，位于7221工作面切眼附近，离注浆孔(注3、注4、注5，位于工作面东北部)较近的村庄未发生明显沉陷，说明注浆充填开采工艺有很好的局域减沉效果。
(a) 2017-11-04—2017-12-10 (b) 2017-11-04—2017-12-22 (c) 2017-11-04—2018-01-03
(d) 2017-11-04—2018-01-15 (e) 2017-11-04—2018-01-27 (f) 2017-11-04—2018-02-08
图1 7221工作面时间序列沉陷分布
Fig.1 Subsidence distribution of 7221 working face in time series
为更好地对D-InSAR监测结果进行定量研究，提取下沉盆地中沿7221工作面倾向l1和走向l2(l1，l2位置如图2所示)的沉陷剖面线，如图3所示。l1长807 m，方向为自西向东，设其最西端为零点；l2长1 089 m，方向为自南向北，设其最南端为零点。从图3(a)可看出，2017-11-04—2018-02-08，l1最大沉陷量为-81 mm，2017-11-28—2018-01-03沉陷剖面线变陡，沉陷量显著增大，后期进入沉陷衰减阶段，开始微小沉陷。从图3(b)可看出，2017-11-04—2018-02-08，l2最大沉陷量为-84 mm，2017-11-28—2018-01-03沉陷速度较快，动态最大沉陷量不断增加，2018-01-03—02-08沉陷速度开始衰减，但平面影响范围不断增大。对比图3(a)和图3(b)可看出，开采对倾向沉陷影响范围小于走向，且倾向沉陷阶段发育早于走向，下沉盆地近似出现平底状态。
2 GM(1,1)与灰色Verhulst模型动态预计
2.1 模型建立与预计
基于D-InSAR监测获取的2017-11-16—2018-02-08 7221工作面地表沉陷量，从各期沉陷分布(以表1中干涉对组合为参考， 2017-11-16—28为第1期观测数据，2017-11-16—12-10为第2期观测数据，以此类推至2018-02-08，共7期数据)中提取43个水准点的像元信息，并绘制沉陷量曲线，如图4所示。选取H2，H3，E2，G5，F4为观测点，对7期数据建立GM(1,1)与灰色Verhulst模型方程，获得前7期沉陷量拟合值,并预计第8，9期(2018-02-20，2018-03-04)沉陷量。1—7期各观测点沉陷实际值与2种GM拟合值对比如图5所示。
图2 l1，l2位置
Fig.2 l1 and l2 positions
(a) 倾向l1
(b) 走向l2
图3 倾向l1和走向l2沉陷剖面线
Fig.3 Subsidence profile curves of tendency l1 and trend l2
2.2 预计结果精度检验
GM精度检验一般分为事中检验和事后检验[16]。事中检验包含残差检验、关联度检验和后验差检验。残差检验属于算术检验，定义为模型值与实际值之间的差，是一种较为客观的逐点检验法；关联度检验属于几何检验，表征2个事物之间的关联程度，关联值越大(越接近1)则相关性越好，模型等级越高；后验差检验属于统计检验，是按照均方差比值C和小误差概率P来评估灰色模型优劣，C越小越好，相反P越大越好。事后检验即预测检验，通常采用残差检验和相对误差检验来了解预测误差。本文对于GM预计结果的精度检验分为2步：采用残差检验、关联度检验和后验差检验对GM(1,1)和灰色Verhulst模型进行拟合精度检验，结果见表2；对预计的第8,9期沉陷值进行残差和相对误差检验，结果见表3和表4。
(a) 水准点位置
(b) 各水准点沉陷量曲线
图4 水准点沉陷量曲线
Fig.4 Subsidence value curves of leveling points
(a) H2
(b) H3
(c) E2
(d) G5
(e) F4
图5 观测点沉陷量曲线
Fig.5 Subsidence value curves of observation points
表2 GM(1,1)和灰色Verhulst模型拟合精度检验
Table 2 Fitting precision checking of GM(1,1) and
grey Verhulst model
观测点GM(1,1)灰色Verhulst模型平均残差/mm关联度后验差(C/P)平均残差/mm关联度后验差(C/P)H25.1250.9910.189/0.8333.7660.9930.148/1H38.1890.9910.162/15.1830.9870.157/1E21.3420.9930.139/11.0220.9950.187/1G51.8130.9950.058/13.0620.9630.111/1F41.4760.9980.081/11.0220.9930.069/1
表3 第8期沉陷量预计精度检验
Table 3 Prediction precision checking of subsidence value
in phase 8
观测点实际值/mmGM(1,1)灰色Verhulst模型预计值/mm残差/mm相对误差/%预计值/mm残差/mm相对误差/%H2-57-671018-54-35H3-73-1103751-881521E2-22-21-15-19-314G5-113-95-1816-89-2421F4-52-51-12-60815
2.3 预计结果分析
(1) 从图5(a)、图5(b)可看出，实际沉陷量曲线后期呈平底饱和状态，此时GM(1,1)和灰色Verhulst模型前期(1,2,3期)拟合结果差别不大，但后期(4—7期)灰色Verhulst模型较GM(1,1)收敛快，拟合结果与实际沉陷量吻合较好。图5(c)中实际沉陷量曲线后期包含非单调摆动过程，使用GM(1,1)拟合误差较大。图5(d)、图5(e)中实际沉陷量曲线呈近似单峰型，GM(1,1)拟合曲线更接近实际值，且后期GM(1,1)和灰色Verhulst模型拟合结果相差不大。
表4 第9期沉陷量预计精度检验
Table 4 Prediction precision checking of subsidence value
in phase 9
观测点实际值/mmGM(1,1)灰色Verhulst模型预计值/mm残差/mm相对误差/%预计值/mm残差/mm相对误差/%H2-57-792238-5812H3-75-1315675-901520E2-26-24-28-20-623G5-113-11744-104-98F4-63-62-12-862337
(2) 从表2可看出，观测点H2,H3使用灰色Verhulst模型拟合结果明显优于GM(1,1)，平均残差分别减小了约26.5%,36.7%。根据后验差精度等级(表5)，观测点H2的GM(1,1)拟合结果只能达到合格，而其灰色Verhulst模型拟合结果精度等级达到良好；观测点G5采用GM(1,1)拟合时平均残差较灰色Verhulst模型减小了1.249 mm，关联度提高了约3.2%；观测点E2，F4采用2种GM拟合的效果相近，精度等级均达到良好，且平均残差之差均小于1 mm。
表5 后验差精度等级
Table 5 Precision level of posterior error
等级均方差比值C小误差概率P良好C<0.35P≥0.95合格0.35≤C<0.50.8≤P<0.95勉强0.5≤C<0.650.7≤P<0.8不合格C≥0.65P<0.7
(3) 从表3、表4可看出，采用GM(1,1)和灰色Verhulst模型预计第8，9期5个观测点的沉陷量时精度与表2中拟合精度较一致，证明可根据各点拟合精度(残差、关联度和后验差)选择沉陷预计适用模型。GM(1,1)对观测点H3的预计误差最大，最大残差绝对值为56 mm，显然该点沉陷量曲线更适宜选用灰色Verhulst模型来描述。灰色Verhulst模型对第8期观测点G5的预计精度最低，最大相对误差为21%，而使用GM(1,1)预计时残差减少了6 mm，且相对误差减小了5%，预计效果更佳。
可见对于单点沉陷曲线后期呈平底饱和状态的观测点(如H2,H3)，地表进入沉陷衰退阶段时选用灰色Verhulst模型进行沉陷动态预计，精度明显优于GM(1,1)，最大相对误差减少了55%；地表进入沉陷发展阶段后，对于单点沉陷量曲线呈近似单峰型的观测点(如G5，F4)，选用GM(1,1)预计沉陷量更合适，最大相对误差可减少35%。
3 结论
(1) 井下采煤地表移动是一个先加速后减速的过程，即随着工作面推进和开采时间推移，地表点会经历开始移动、剧烈移动、缓慢移动及停止移动4个阶段，其移动轨迹是一个随机变化的过程。
(2) 矿山开采沉陷开始至活跃前期，地表单点沉陷量曲线呈近似单峰型，采用GM(1,1)进行短期预计效果更佳；当矿区地表单点沉陷量曲线呈平底饱和状态，沉陷进入衰退阶段，此时采用灰色Verhulst模型预计的结果与实际值误差较小，适用于中长期预计。
(3) GM(1,1)和灰色Verhulst模型较传统地表沉陷预计方法而言，是假设地表沉陷历程符合已知时间序列所描述的动态过程，通过先验信息来研究或预测未知领域的时间函数模型，而开采引起的地表沉陷是一个复杂的时空变化过程，仅依靠时间序列模型并不能实现准确预计，后续会将时间序列模型与某种常规静态方法(如概率积分法)相结合，建立更精确的地表任意点动态沉陷预计模型。
[1] 王诺,张进,吴迪,等.世界煤炭资源流动的时空格局及成因分析[J].自然资源学报,2019,34(3):487-500.
WANG Nuo,ZHANG Jin,WU Di,et al.The temporal and spatial patterns and causes of coal resource flow in the world[J].Journal of Natural Resources,2019,34(3):487-500.
[2] 高永格,牛矗,张强,等.厚松散层下采煤地表沉陷特征研究[J].煤炭科学技术,2019,47(6):192-198.
GAO Yongge,NIU Chu,ZHANG Qiang,et al.Study on surface subsidence characteristics of coal mining under thick loose layer[J].Coal Science and Technology,2019,47(6):192-198.
[3] 陈博.开采沉陷动态预计方法研究[D].淄博:山东理工大学,2013.
CHEN Bo.Research on dynamic prediction method of mining subsidence[D].Zibo:Shandong University of Technology,2013.
[4] 刘玉成,戴华阳.近水平煤层开采沉陷预计的双曲线剖面函数法[J].中国矿业大学学报,2019,48(3):676-681.
LIU Yucheng,DAI Huayang.Hyperbolic function model for predicting the main section surface deformation curve due to approximate horizontal coal seam underground longwall mining[J].Journal of China University of Mining and Technology,2019,48(3):676-681.
[5] 闫伟涛,陈俊杰,阎跃观.倾斜煤层开采沉陷预计模型的构建[J].河南理工大学学报(自然科学版),2018,37(3):12-16.
YAN Weitao,CHEN Junjie,YAN Yueguan.Subsidence prediction model construction in inclined coal seam mining[J].Journal of Henan Polytechnic University(Natural Science),2018,37(3):12-16.
[6] 陈海燕,戎晓力,林阳.矿区开采沉陷预计的改进BP神经网络模型[J].金属矿山,2017(4):119-122.
CHEN Haiyan,RONG Xiaoli,LIN Yang.Mining subsidence prediction based on improved BP neural network model[J].Metal Mine,2017(4):119-122.
[7] 邓聚龙.灰色理论基础[M].武汉:华中科技大学出版社,2002.
DENG Julong.Grey theoretical basis[M].Wuhan: Huazhong University of Science and Technology Press,2002.
[8] 潘申运.补偿最小二乘法改进灰色预测模型的应用分析[J].测绘科学,2014,39(11):104-107.
PAN Shenyun.Application of penalized least squares method in gray forecasting model[J].Science of Surveying and Mapping,2014,39(11):104-107.
[9] 姜刚,杨志强,张贵钢.卡尔曼滤波算法的灰色理论模型在变形监测中的应用[J].测绘科学,2011,36(4):19-21.
JIANG Gang,YANG Zhiqiang,ZHANG Guigang.Application of gray model in deformation monitoring based on Kalman filter method[J].Science of Surveying and Mapping,2011,36(4):19-21.
[10] 肖海平,陈兰兰.灰色理论模型在矿山变形监测中的应用[J].金属矿山,2009(1):154-155.
XIAO Haiping,CHEN Lanlan.Application of grey theory model in monitoring mine deformation[J].Metal Mine,2009(1):154-155.
[11] 王磊,张鲜妮,陈元非,等.基于D-InSAR LOS向变形的开采沉陷预计参数反演方法研究[J].中国矿业大学学报,2017,46(5):1159-1165.
WANG Lei,ZHANG Xianni,CHEN Yuanfei,et al.Method of mining subsidence prediction parameters inversion based on D-InSAR LOS deformation[J].Journal of China University of Mining & Technology,2017,46(5):1159-1165.
[12] HE Liming,WU Lixin,LIU Shanjun,et al.Mapping two-dimensional deformation field time-series of large slope by coupling D-InSAR-SBAS with MAI-SBAS[J].RemoteSensing,2015,7(9):12440-12458.
[13] CARNEC C,MASSONNET D,KING C.Two examples of the use of SAR interferometry on displacement fields of small spatial extent[J].Geophysical Research Letters,1996,23(24):3579-3582.
[14] 朱建军,杨泽发,李志伟.InSAR矿区地表三维形变监测与预计研究进展[J].测绘学报,2019,48(2):135-144.
ZHU Jianjun,YANG Zefa,LI Zhiwei.Recent progress in retrieving and predicting mining-induced 3D displacements using InSAR[J].Acta Geodaetica et Cartographica Sinica,2019,48(2):135-144.
[15] GABRIEL A K,GOLDSTEIN R M,ZEBKER H A.Mapping small elevation changes over large areas:differential radar interferometry[J].Journal of Geophysical Research,1989,94(B7):9183-9191.
[16] 杨俊凯,范洪冬,赵伟颖,等.基于D-InSAR技术和灰色Verhulst模型的矿区沉降监测与预计[J].金属矿山,2015(3):143-147.
YANG Junkai,FAN Hongdong,ZHAO Weiying,et al.Monitoring and prediction of mining subsidence based on D-InSAR and gray Verhulst model[J].Metal Mine,2015(3):143-147.