传统的矿区开采沉陷变形监测常采用精密水准测量、全站仪三角测量和GPS测量等方法[1-4]。但常规的监测方法存在费用高、需要大量人力物力、监测周期长等缺点，而D-InSAR技术作为在SAR基础上发展而来的新型空间对地观测技术，具有连续空间覆盖能力、高度自动化和高精度的优势，近年来在矿区开采沉陷中得到广泛的应用[5-6]。D-InSAR技术主要是获取监测区域内地表形变信息，地表形变是一个长期持续的过程，大多数研究方法是采用多幅影像进行差分干涉处理，得到地表时序形变信息。在此基础上提出了干涉叠加技术，主要包括永久散射体干涉测量(PS-InSAR)和小基线集技术(SBAS-InSAR)[7-10]。两种监测方法都要求获得连续时间序列的SAR影像数据，因此虽然D-InSAR技术在矿区开采沉陷变形监测方面具有很多成功且典型的工程实例，但在具体的实践中和后续数据处理中受到很多限制性约束，对于矿山开采沉陷监测而言，SAR影像数据不足的情况下主要影响因素是时间失相干，从而不能获得连续时间段的地表下沉值，因此SAR影像源数据的获取和影像费用高的问题成为影响该技术的发展前景之一[11]。为了克服这一缺陷，本文提出将合成孔径雷达差分干涉测量技术(D-InSAR)与三次样条插值法相融合的开采沉陷监测方法。该方法利用D-InSAR影像上对应监测点在已知雷达卫星重访周期内的监测值作为三次样条插值的节点，通过三次样条插值法来内插反演其他雷达重访周期内监测点的下沉值。
1 D-InSAR沉陷监测原理与内插反演模型的建立 1.1 D-InSAR测量技术原理合成孔径雷达差分干涉测量(D-InSAR)技术以合成孔径雷达(SAR)复数据所提供的相位及强度信息为信息源，利用同一目标区域的两幅或多幅干涉纹图获取SAR的相位信息，从而提取地表三维信息和高程变化信息，然后通过差分处理(除去地球表面、地形起伏等因素)来获取地表微小形变[12-13]。
目前，D-InSAR获取地表形变信息的实现方法主要有3种：二轨法、三轨法和四轨法[14-17]。二轨法是获取同一地区形变发生前后的两幅SAR影像，将其进行干涉处理，干涉结果与已有DEM数据进行差分，消除地形因素而获取地表形变信息的方法；三轨法需要3幅SAR影像，其中两幅进行差分干涉处理生成干涉图，作为DEM信息，再将另一幅影像与生成的干涉图进行干涉处理，最终得出地表形变信息；四轨法基本思想是获取形变前后的两对SAR影像，分别进行干涉处理，形成形变前后的两幅干涉相位图，形变前的主要用来生成DEM，形变后的主要用来反映形变信息，将二者进行差分，得到形变相位，再将相位转换为斜距从而计算出雷达视线向的形变量。
1.2 融合D-InSAR监测值的三次样条插值沉陷内插反演函数的建立三次样条插值函数是最常用的插值曲线拟合函数，其良好的收敛性、可靠的稳定性，以及具有二阶光滑度等优点使其在函数逼近、微积分和微分方程等科学计算中应用广泛。在矿区开采沉陷中，由于沉降是一个随时间平滑变化的过程，因此，三次样条插值在开采沉陷数据处理中起到非常重要的作用。
本文利用D-InSAR技术获得(矿区沉降面上)一系列与时间相关的离散数据{xj}＝{x0, x1, …, xn}，xj表示雷达卫星重访周期。使用三次样条插值进行沉陷监测，建立已有观测数据和未知预测数据之间的函数关系。以区间[xj, xj+1]为例，S (x)的二阶导数为S″(x)=Mj(j=0, 1, …, n)，则该区间上S (x)的三次样条插值函数为
(1) 若令
(2) 线性方程组的解为M，代入式(2) 中可以得到各个小区间上的三次样条插值函数式[18-19]。采用三次样条插值建立开采沉陷监测函数步骤如下：
(1) 为构建三次样条插值所需要的已知观测数据，首先对已有的雷达影像进行差分干涉处理，得到一系列影像周期内监测点的下沉值。以与首幅影像日期(2012-11-01) 间隔24 d为第一个周期开始起算，将得到的下沉值作为三次样条插值函数的已知数据点。
(2) 运用MATLAB软件进行编程，使用三次样条插值函数spline，建立内插模型
其中，x、y为已知雷达卫星重访周期及该时间段内下沉值(D-InSAR监测下沉值)；xx是插值节点(雷达卫星重访周期数)；Ppval (cs, xx)是插值节点函数(估测下沉值)。
2 D-InSAR数据处理和结果分析 2.1 研究区概况及相关数据本次试验采用淮北矿业集团袁店二矿作为研究区域。在2012年11月-2013年12月对7226工作面进行了开采工作，7226工作面长约768 m，工作面平均宽为158.9 m，工作面标高为-453~-416 m，煤层平均厚度3.89 m，煤层倾角为3°~14°。由于该工作面上方有村庄，进一步的地下采煤会引起地表建筑物的损坏，因此在该工作面采用注浆充填开采方式，从2013年7月开始进行注浆开采。2013年10月-2015年5月对7225工作面进行了开采工作，7225工作面长约960 m，工作面平均宽为176 m，工作面平均标高为-450 m，煤层平均厚度为4.07 m，煤层倾角为6°~19°。
本文采用9景C波段RADARSAT-2影像数据对试验矿区进行研究分析，影像分辨率为3 m，雷达卫星重访周期为24 d。本文选用的外部DEM是SRTM3 DEM数据，是美国奋进号航天飞机在2000年2月11-22日进行的为期11 d的航天飞机雷达地形测绘任务SRTM所获得的。为减少时间失相关，尽可能选取时间间隔最小的两幅影像进行干涉处理，9景影像共组成5组干涉对，其相关参数见表 1。
| 序号 | 主影像日期 | 辅影像日期 | 极化方式 | 时间基线/d | 空间基线/m |
| 1 | 2012-11-01 | 2012-12-19 | VV | 48 | 151 |
| 2 | 2013-11-13 | 2013-12-07 | HH | 24 | 148 |
| 3 | 2014-02-24 | 2014-03-20 | VV | 24 | 187 |
| 4 | 2015-04-01 | 2015-05-19 | HH | 48 | 317 |
| 5 | 2015-05-19 | 2015-12-21 | HH | 216 | 245 |
为了得到矿区地表下沉形变场，采用SARscape软件中D-InSAR Displacement Workflow工作流分别对5组干涉对进行二轨差分处理，经过多视和滤波抑制斑点噪声，相位解缠后进行地理编码和辐射定标得到形变图，将得到的5组形变图和在7225、7226工作面设立的地表移动观测站Cass图导入ArcMap中进行叠加分析，得到5个时间段内工作面地表形变图，如图 1所示。
|
| 图 1 7225、7226工作面不同时间段内沉降分布 |
(1) 由图 1(a)可知，形成的最大下沉区域位于距离7226工作面340 m的7221工作面，该工作面开采起止时间是2011年1月-12月，最大下沉区域是由于该工作面开采残余变形的影响，并非7226工作面开采引起，而7226工作面在开采初期沉降值和影响范围很小。
(2) 由图 1(b)可知，注浆减沉开采的影响导致注浆周边区域下沉分布不均匀，在2013年11月-12月期间，地表最大下沉区域和影响范围随着7226工作面开采逐渐向前推进，形成了明显的形变区沉降漏斗。
(3) 由图 1(c)可知，由于7225工作面开采与7226工作面残余变形叠加影响，造成下沉影响范围较大的区域分别出现在7226工作面的尾端和7225工作面的开切眼处，而在7226工作面注浆区域地表下沉值依然分布不均匀。
(4) 由图 1(d)可知，随着7225工作面采煤的推进，下沉影响范围增加，而下沉值最大区域不在7225工作面正上方，而是在两个工作面之间，这是由于受到7226工作面开采的残余变形的影响。
(5) 由图 1(e)可知，随着两个工作面采煤工作的结束，地表受残余变形的影响，在216 d内下沉范围增广明显，但下沉速率开始减缓，下沉值大于10 mm等值线的区域平均沉降速率为0.23 mm/d，属于下沉衰退期[20]。
3 结合三次样条插值函数的开采沉陷监测为了克服SAR影像数据不足而存在的时间失相干现象，只能得到7225、7226工作面5个时间段内的下沉值和影响范围，不能获取2012年11月-2015年12月时间内的连续下沉影响范围。本文采用三次样条插值方法，将D-InSAR监测值作为三次样条插值函数的初始值。根据开采沉陷相关知识可知，开采工作面的走向线和倾向线方向的下沉值最能代表地表的下沉规律，因此选取走向线上L系列监测点中的ML02、ML06、ML09、ML11和倾向线S系列监测点中的MS01、MS04、MS08、MS12为研究对象，如图 1所示。以雷达卫星一个重访周期的时间间隔为三次样条插值函数的自变量步长，将1号和4号干涉对中两个周期内的下沉值作平均值处理，所选取的8个监测点下沉信息见表 2。
| mm | ||||||||
| 雷达重访周期数(x步长) | 下沉值(y) | |||||||
| ML02 | ML06 | ML09 | ML11 | MS01 | MS04 | MS08 | MS12 | |
| 1 | 1.9 | 4.1 | 5.5 | 4.3 | 2.6 | 4.9 | 5.5 | 5.3 |
| 2 | 1.9 | 4.1 | 5.5 | 4.3 | 2.6 | 4.9 | 5.5 | 5.3 |
| 18 | 1.1 | 10.2 | 16.4 | 18.1 | 7.1 | 5.4 | 2.7 | 3.4 |
| 21 | 3.1 | 20.9 | 30.8 | 35.5 | 13.2 | 10.7 | 6.1 | 5.2 |
| 37 | 2.8 | 5.8 | 8.2 | 9 | 0.4 | 1 | 0.7 | 0.6 |
| 38 | 2.8 | 5.8 | 8.2 | 9 | 0.4 | 1 | 0.7 | 0.6 |
利用三次样条插值法求出2012年11月-2015年5月上述8个水准点每个雷达卫星重访周期内的下沉值，运用MATLAB软件绘制地表下沉曲线图，如图 2所示。由图 2可知，雷达卫星重访周期内该区域地表下沉趋势略有不同，符合实际地表下沉规律。最大下沉值出现在靠近中心盆地监测点ML11处，单个周期内最大下沉量为42.4 mm。
|
| 图 2 第11~26个雷达卫星重访周期内所选监测点下沉值 |
为了进一步验证精度，实地观测了2013年7月23日(对应SAR重访周期第11周)-2014年7月18日(对应SAR重访周期第26周)7225、7226工作面16个周期(共360 d)内的水准数据，将内插反演值与实测值对比，结果见表 3。由表 3可知，插值解算累计下沉值与实测结果基本一致，其中最大误差出现在监测点ML11处，为31.5 mm，其余监测点误差均在20 mm以内。最大相对误差出现在监测点MS01处，为17%，其余监测点处相对误差都不大于8.2%。
| 点号 | 沉降值 | 误差/mm | 相对误差/(%) | |
| 预计值/mm | 实测值/mm | |||
| ML02 | 32.7 | 34 | 1.3 | 3.8 |
| ML06 | 214.1 | 217 | 2.9 | 1.3 |
| ML09 | 324.7 | 312 | 12.7 | 4.1 |
| ML11 | 369.5 | 401 | 31.5 | 7.9 |
| MS01 | 135.7 | 116 | 19.7 | 17 |
| MS04 | 112.5 | 105 | 7.5 | 7.1 |
| MS08 | 61.7 | 57 | 4.7 | 8.2 |
| MS12 | 64.8 | 61 | 3.8 | 6.2 |
(1) 本文采用D-InSAR技术对淮北矿业集团袁店二矿7225、7226工作面开采沉陷进行了分析，获取了5个时期内地表沉陷的下沉影响范围，结果与实际矿区开采导致地表下沉影响范围符合得很好。
(2) 融合D-InSAR技术与三次样条插值的方法有效地解决了在SAR影像数据不足的情况下对袁店二矿沉陷区进行监测分析的问题。应用结果表明：该方法内插解算反演结果较实测水准监测结果最大误差和最大相对误差分别为31.5 mm和17%，内插反演结果与实测值比较吻合，从而有效地说明了该方法的可行性。
| [1] | 岳建平, 田林亚. 变形监测技术与应用[M]. 北京: 国防工业出版社, 2007. |
| [2] | WEBLEY P W, BINGLEY R M, DODSON A H, et al. Atmospheric Water Vapour Correction to InSAR Surface Motion Measurements on Mountains:Results from a Dense GPS Network on Mount Etna[J]. Physics and Chemistry of the Earth, 2002, 27(4): 363–370. |
| [3] | 陈炳乾, 邓喀中, 范洪冬. 基于D-InSAR技术和SVR算法的开采沉陷监测与预计[J]. 中国矿业大学学报, 2014, 43(5): 880–886. |
| [4] | CHEN Bingqian, DENG Kazhong, FAN Hongdong, et al. Large-scale Deformation Monitoring in Mining Area by D-InSAR and 3D Laser Scanning Technology Integration[J]. International Journal of Mining Science and Technology, 2013, 23(4): 555–561. DOI:10.1016/j.ijmst.2013.07.014 |
| [5] | BAYUAJI L, SUMANTYO J T S, KUZE H. ALOS PALSAR D-InSAR for Land Subsidence Mapping in Jakarta, Indonesia[J]. Canadian Journal of Remote Sensing, 2010, 36(1): 1–8. DOI:10.5589/m10-023 |
| [6] | PEI Liang, LI Wenjie, TAN Yang. Study of Monitoring Mining Subsidence in Coal Mining Area by D-InSAR Technology[J]. Journal of Coal Science and Engineering, 2008, 14(4): 591–593. DOI:10.1007/s12404-008-0417-2 |
| [7] | 杨成生, 刘媛媛, 敖萌. 基于SBAS时序分析的大同地面沉降与地下水活动研究[J]. 国土资源遥感, 2015, 27(1): 127–132. DOI:10.6046/gtzyyg.2015.01.20 |
| [8] | 侯安业, 张景发, 刘斌, 等. PS-InSAR与SBAS-InSAR监测地表沉降的比较研究[J]. 大地测量与地球动力学, 2012, 34(4): 125–128. |
| [9] | 李永生, 张景发, 罗毅, 等. 利用高分辨率聚束模式TerraSAR-X影像的PSInSAR监测地表变形[J]. 武汉大学学报(信息科学版), 2012, 27(12): 1452–1455, 1514. |
| [10] | 许才军, 何平, 温扬茂. 利用PSInSAR研究意大利Etna火山的地表形变[J]. 武汉大学学报(信息科学版), 2011, 26(9): 1012–1016. |
| [11] | 吴立新, 高均海, 葛大庆, 等. 基于D-InSAR的煤矿区开采沉陷遥感监测技术分析[J]. 地理与地理信息科学, 2004, 20(2): 22–25. |
| [12] | 范洪冬, 邓喀中, 祝传广, 等. 基于时序SAR技术的采空区上方高速公路变形监测及预测方法[J]. 煤炭学报, 2012, 37(11): 1841–1846. |
| [13] | ALLAN M L, MUTHAMA M N, KINYANJUI S Z. Time Series Analysis Model for the Rate of Influx of Refugees in Kenya[J]. The International Journal of Engineering and Science, 2013, 2(9): 7–18. |
| [14] | 魏长婧, 汪云甲, 闫建伟. D-InSAR技术二轨法监测矿区地表沉陷的方法研究[J]. 煤炭技术, 2012, 31(7): 129–130. |
| [15] | 王小兵. 基于DInSAR技术的矿山开采沉陷监测研究现状[J]. 金属矿山, 2015(S1): 65–71. |
| [16] | SHARIFIKIA M. A Comprehensive Interferometric Process for Monitoring Land Deformation Using ASAR and PALSAR Satellite Interferometric Data[J]. GIScience & Remote Sensing, 2015, 52(1): 58–77. |
| [17] | 刘国林, 郝华东, 陶秋香. 卡尔曼滤波相位解缠及其与其他方法的对比分析[J]. 武汉大学学报(信息科学版), 2010, 25(10): 1174–1178. |
| [18] | 陈文略, 王子羊. 三次样条插值在工程拟合中的应用[J]. 华中师范大学学报(自然版), 2004, 38(4): 418–422. |
| [19] | 张威, 杨月婷. 数值分析[M]. 5版.北京: 清华大学出版社, 2010: 41-44. |
| [20] | 高延法, 贾君莹, 李冰, 等. 地表下沉衰减函数与塌陷区稳定性分析[J]. 煤炭学报, 2009, 34(7): 892–896. |