基于D-InSAR与概率积分法的开采沉陷监测与预计

杨俊凯 陈炳乾 邓喀中 何 强 赵伟颖

(1.中国矿业大学环境与测绘学院，江苏 徐州 221116；2.矿山空间信息技术国家测绘地理信息局重点实验室,河南 焦作 454000)

基于D-InSAR与概率积分法的开采沉陷监测与预计

杨俊凯1陈炳乾2邓喀中1何 强1赵伟颖1

(1.中国矿业大学环境与测绘学院，江苏 徐州 221116；2.矿山空间信息技术国家测绘地理信息局重点实验室,河南 焦作 454000)

为掌握采动区附近高等级公路的变形趋势，针对传统开采沉陷监测方法存在诸多缺陷的问题，将合成孔径雷达技术(D-InSAR)与概率积分法相结合，提出了一种矿山开采沉陷监测与预计方法。该方法首先利用D-InSAR技术对兖州某矿9308工作面进行了监测和分析，并将监测结果与路面水准测量数据进行了对比，然后基于D-InSAR的监测值，求取了9308工作面非充分采动下的概率积分法预计参数，再根据三下采煤规程对参数进行修正得到工作面全采时的预计参数，最后利用修正后的参数预计出工作面全采对附近高等级公路的影响程度。结果表明D-InSAR技术高精度的优势可以为开采沉陷预计提供可靠的监测数据，该方法可以为实现矿山开采沉陷监测与预计一体化提供有力支持。

D-InSAR 概率积分法 开采沉陷 预计参数

煤炭是我国的主体能源，占据着我国能源消费结构中的70%左右，对于穿越高等级公路的采空区，煤层开采极易对公路造成严重的破坏，为了保护矿区附近重要建筑物，必须建立及时有效的开采沉陷监测与预计系统。D-InSAR技术作为近几年新兴的空间对地观测技术之一，已经在矿区开采沉陷监测中取得了广泛的应用，1996年Camec等人利用D-lnSAR技术监测了法国Gardanne附近煤矿的地表沉陷,证明D-lnSAR可以监测矿区地表沉陷[1]；其后国内外学者针对D-InSAR技术的不足做了许多改进,由此发展出了PS-InSAR、CPT、小基线集等技术[2-4],这些新技术往往在城市微小缓慢形变中能取得理想的效果，而在新兴的矿区开采沉陷监测手段中，D-InSAR技术仍然有着不可替代的作用。而在开采沉陷预计方面，目前最常用的预计方法有概率积分法，典型曲线法和负指数函数法[5-6]，但这些预计方法都是基于传统的测量手段获取监测资料，传统的测量手段虽然精度较高,却存在工作量大、成本高、测点稀疏、人身安全不易保证等缺点[7]。针对传统开采沉陷监测方法的诸多缺陷，为实现开采沉陷监测与预计一体化，本研究采用10景高分TerraSAR数据，利用D-InSAR技术获取了山东兖州某矿9308工作面上方时序沉降图，通过提取9308工作面走向和倾向2条观测线的下沉曲线，分析了采空区上方地表移动规律，然后对2条下沉曲线进行了概率积分法拟合求参，结合三下采煤规程中预计参数与回采尺寸的关系，对预计参数进行了修正，最后利用修正后的参数对工作面全采时进行了稳态预计，预计出了工作面全采对附近高等级公路的影响程度，从而实现了矿区开采沉陷监测与预计一体化。

1 矿区开采沉陷监测与预计模型的建立

1.1 D-InSAR技术原理及监测矿区开采沉陷可行性分析

合成孔径雷达差分干涉技术(D-InSAR)是以合成孔径雷达复数图像的相位信息获取地表变化信息的技术[8]。D-InSAR技术常用的处理方法主要有二轨法和三轨法，在偏远山区，外部DEM不易获取，常常使用三轨D-InSAR技术，但是三轨差分法容易引入解缠相位噪声[9]。二轨差分法是由Massonnet 等[10]提出的，它首先利用研究区域地表变化前后的2幅SAR影像生成干涉条纹图，再用外部DEM数据模拟地形相位的条纹图，从干涉条纹图中去除地形信息，从而获取地表形变信息。

D-InSAR采用重复轨道干涉测量模式，在其获取的2幅影像干涉纹图中，还包含了其他误差，干涉相位主要由以下几部分组成[11]：

(1)

式中，Φtop为地形相位，Φflt为平地相位，Φdef为视线方向(LOS)的形变相位，Φatm为大气延迟相位，Φniose为噪声相位，k为整周模糊度。其中Φtop和Φflt可以通过外部DEM来去除，Φatm可以用干涉图叠加的方法来提高形变信息和大气干扰信号之间的信噪比，降低其相位分量[12]。对于噪声相位Φniose，一般采用多视处理和滤波抑制干涉图的噪声。而整周模糊度可通过相位解缠来还原真实相位，常用的解缠方法有枝切法[13]和最小费用流法[14]。最后得出了沿视线方向的地表形变相位。

在矿区开采沉陷监测应用中，D-InSAR主要受制于矿区地表形变梯度大、速度快，地表相干性低造成的失相干问题。对于地表植被引起的失相干，选取冬季、初春植被相对稀少的影像，可提高干涉对相干性。而由于形变梯度造成的空间失相干，仍然是D-InSAR技术的难点。不同的SAR卫星的临界梯度主要受卫星的空间分辨率以及波长的影响。D-InSAR技术可监测到的临界形变量可表示为[15-16]

(2)

式中，Wmax为地表可监测的临界形变量；r为地表下沉主要影响半径;g为雷达影像的地面分辨率,λ为传感器波长。以TerraSAR卫星为例，其重访周期为11d，若主要影响半径r取100 m，根据式(2)可知，该卫星可监测最大沉降量为260 mm，可监测平均下沉速度为23.6 mm/d。本研究中煤层埋深620 m左右，属于较大深部开采[17]，加之采空区上覆岩层坚硬，地表下沉量相对较小，D-InSAR技术可以在本研究实验区内监测出完整的下沉盆地，可为概率积分法模型提供监测资料，从而实现矿区开采沉陷监测与预计集成化。

1.2 概率积分法原理

概率积分法是由我国学者刘宝琛、廖国华根据波兰学者李特威尼申(J.Litwiniszyn)的随机介质理论发展得来的[6]，对于半无限开采时，主断面内任意点地表移动的变形预计公式为

(3)

2 实 验

2.1 实验区域及数据

为了验证本研究所提出模型的可行性，以山东兖州某矿9308工作面为例，该工作面地表多被农田覆盖，并且工作面正南方有2010年开采完的9306工作面和某高等级公路，9308工作面的开采必将对附近公路产生一定的影响。9308工作面设计切眼面长140m，推进长度800m，设计开采时间为2011年11月到2012年7月。煤层平均采厚5m，平均采深620m，煤层倾角为6°，开采方向自东北向西南方向，如图1所示。

图1 工作面位置

为了有效地实现开采沉陷监测，选取实验区地表植被相对稀少时段的10景TerraSAR卫星影像，同时为了减小时间的失相干，将10景影像两两配对，共组成9组干涉对，具体影像参数见表1。

2.2 时序D-InSAR开采沉陷监测实验

表1 干涉对组合与对应参数

为了实现实验区域时序D-InSAR的开采沉陷监测，采用GAMMA软件分别对9组干涉对进行二轨差分处理，此外D-InSAR易受噪声的影响，采用自适应滤波方法多次滤波，来改善干涉条纹的质量。同时在去地形相位时选用SRTM90 m分辨率的高程数据，其垂直精度约为16 m，在二轨差分中，外部DEM10 m的垂直误差只会导致对形变量估计约为0.9 mm的偏差，所以外部DEM的误差可忽略[18]。然后将9组形变图进行叠加，得到时序D-InSAR形变图，最后经过地理编码转换成UTM坐标系的时序形变结果，如图2所示。

图2 D-InSAR时序形变图

图2(a)～图2(i)所示为实验区域2011年12月25日—2012年4月23日时序D-InSAR结果，在图2(a)中，工作面上方未有明显沉陷区域，表明虽然工作面起始开采时间为2011年11月，距起始监测时间工作面已经推进约100 m，但是由于工作面推进距离小于起动距(0.25～0.5)H0，所以图2(a)中未形成明显的下沉区域。随着工作面的推进，有明显的下沉盆地，并且下沉盆地影响范围慢慢变大，直至不再增大。在整个监测时段内下沉盆地中心最大下沉值达到了682 mm，如图2(i)，下沉盆地中心由于受工作面倾角和附近9306工作面重复采动的影响，明显的向下山方向和9306工作面偏移。

为了定量地分析矿山开采沉陷的地表形变规律，沿9308工作面走向、倾向各布置1条观测线，如图3所示，然后在时序D-InSAR形变图中分别提取出2条观测线上的下沉值，如图4。从图4(a)可以看出，走向始终处于非充分采动状态，从2011年12月25日至2012年4月2日走向累计最大下沉达到663 mm，在监测初期，下沉值普遍较小，说明未受提前开采的影响，可以不考虑其提前采动影响，下沉速度在2011-12-25—2012-01-16监测时段处于增加的状态，在2012-01-16—2012-01-27下沉速度明显减小，后期下沉速度出现了减小后增加，分析可能是由于受到工作面推进速度的影响。在图4(b)中看出，在监测时段内倾向线处于非充分采动状态，倾向整个监测时段累计最大的下沉达到682 mm，其下沉速度趋势与走向一致，随着工作面的开采，倾向下沉曲线的最大下沉值沿着工作面的推进方向向前移动。在走向方向下沉影响范围已经成形，不再有明显的变化，后期只是下沉量值的增加，如图4(a)所示。

图3 观测线A、B以及路面水准点位置图

图4 走向A和倾向B下沉曲线

2.3 D-InSAR监测精度评定

为了验证D-InSAR监测精度，取布设于高等级公路一侧的21个保存比较好的水准点进行验证，具体点位如图3所示中点L01-L21。这21个水准点于2012年1月5日至2012年4月3日进行了2次观测，其观测时间与D-InSAR监测时段2012年1月5日至2012年4月2日较为一致，以稳定点L01为参考点，求取了其余20个点相对于L01的沉降量，具体对比结果如图5所示。

图5 水准与D-InSAR监测值对比

由图5可看出，D-InSAR与水准测量值整体上保持一致，D-InSAR监测值相对于水准值的最大误差出现在点L21，为20.6 mm，最小误差在点L17，为0.1 mm，其他点误差分布较均匀，中误差为7.7 mm。经过以上验证，证明D-InSAR可以满足矿区开采沉陷监测的精度要求。

2.4 融合概率积分法的开采沉陷预计

为了有效地实现开采沉陷预计,选取我国应用最广泛的概率积分法进行开采沉陷预计，每隔20 m等距提取9308工作面整个监测时段内走向和倾向曲线上的下沉信息，共计69个监测点，如图6所示，利用开采沉陷软件(MSPS)进行预计参数，求得9308工作面开采到2012年4月2日时的下沉系数q=0.33，主要影响角正切tanβ=2.28，水平移动系数b=0.3，最大下沉角θ=78.25°，拐点偏移距S左=-37.18 m，S右=33.30 m，S上=-0.99 m，S下=-2.53 m。

图6 D-InSAR提取值与MSPS拟合值对比

由图4可知，9308工作面在整个监测时段走向和倾向均处于非充分采动，经过概率积分法求参，只得到了非充分采动情况下的预计参数。在《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规程》中指出,当回采尺寸较小时，要得到充分采动情况的预计参数，应对非充分采动预计参数进行修正，其中水平移动系数和开采影响传播角与回采尺寸关系不明显，可以不用修正。采用规程中下沉系数、主要影响角正切、拐点偏移距与回采尺寸的关系并结合该矿区移动覆岩岩性进行预计参数修正，修正后的参数q=0.47， tanβ=2.3，b=0.3，θ=78.25°，S左=-41.3 m，S右=37 m，S上=-1.1 m，S下=-2.5 m。最后采用修正后的参数进行9308工作面全采时概率积分法预计，如图7所示。

由图7可知，当工作面全采时，采用概率积分法预计出高等级公路受采动影响长度约为982 m，最大下沉出现在L8点，最大下沉值1 017 mm，最大倾斜在L3点，达到4.2 mm/m,最大曲率0.04 mm/m2,最大拉伸变形1.8 mm/m,最大压缩变形3.5 mm/m,倾斜和水平移动略微超过了参考文献[19]中规定的一级公路地基移动变形容许值，已经对公路造成了轻度影响，所以建议协调好路、矿有关部门及时对路面进行定期监测和修复，以保障行车安全。

图7 概率积分法预计的等值线图

3 结 论

本研究利用D-InSAR技术对9308工作面进行了开采沉陷监测，并且分析了采空区上方地表移动规律，同时将D-InSAR监测值与水准实测值，对比发现，D-InSAR监测值相对于水准值的最大误差为20.6 mm，最小误差仅为0.1 mm。在保证D-InSAR监测精度的前提下，获取了9308工作面非充分采动下的概率积分法预计参数，根据“三下”采煤规程对预计参数进行修正，最后预计出了工作面全采时对附近高等级公路的影响程度预计得出高等级公路最大下沉值达到1 017 mm，最大倾斜达到4.2 mm/m,最大曲率0.04 mm/m2,最大水平变形3.5 mm/m。结果证明融合D-InSAR技术和概率积分法可以有效地实现矿山开采沉陷监测与预计的集成化，对矿山安全生产具有重要意义。

[1] 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.

[2] Ferretti A,Claudio P,Rocca A.Nonliner subsidence rate estimation using permanent scatters in differential SAR interferometry[J].IEEE Transactions on Geoscience and Remote Sensing，2000,38(5):202-212.

[3] Mora O，Mallorqui J J，Broquetas A.Linear and nonlinear terrain deformation maps from a reduced set of interferometric SAR images[J].IEEE Transactions on Geoscience and Remote Sensing，2003，41(10)：2243-2253.

[4] Berardino P，Fornaro Q,Lanari R，et al.A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[J].IEEE Transactions on Geoscience and Remote Sensing,2002,40(11)：2375-2383.

[5] 邹友峰,邓喀中,马伟民.矿山开采沉陷工程[M].徐州：中国矿业大学出版社，2003:69-103. Zou Youfeng，Deng Kazhong，Ma Weimin.Mining Subsidence Engineering in Mine[M].Xuzhou：China University of Mining and Technology Press，2003：69-153.

[6] 何国清,杨 伦,凌赓娣,等.矿山开采沉陷学[M].徐州:中国矿业大学出版社,1995:72-90. He Guoqing，Yang Lun，Ling Gengdi，et al.Mining Subsidence Theory[M].Xuzhou：China University of Mining and Technology Press,1995:72-90.

[7] 范洪冬,邓喀中,祝传广,等.基于时序SAR技术的采空区上方高速公路变形监测及预测方法[J].煤炭学报，2012，37(11)：1841-1846. Fan Hongdong，Deng Kazhong，Zhu Chuanguang，et al.Deformation monitoring and prediction methods for expressway above goaf based on time series SAR technique[J].Journal of China Coal Society,2012，37(11)：1841-1846.

[8] 焦明连,蒋廷臣.合成孔径雷达干涉测量理论与应用[M].北京:测绘出版社,2009. Jiao Minglian，Jiang Tingchen.Synthetic Aperture Radar Interferometry Measurement Theory and Application[M].Beijing:Surveying and Mapping Press,2009.

[9] 周建民,李 震,李新武.基于ALOS /PALSAR雷达干涉数据的中国西部山谷冰川冰流运动规律研究[J].测绘学报,2009,38(4):341-347. Zhou Jianmin，Li Zhen，Li Xinwu.Research on rules of the valley glacier motion in western China based on ALOS/PALSAR interferometry[J].Acta Geodaetica et Cartographica Sinica.2009,38(4):341-347.

[10] Massonnet D，Rossi M，Carmona C，et al.The displacement field of the Landers earthquake mapped by radar interferometry[J].Nature，1993，364：138-142.

[11] 朱建军,邢学敏,胡 俊,等.利用InSAR技术监测矿区地表形变[J].中国有色金属学报，2011,21(10)：2564-2576. Zhu Jianjun，Xing Xuemin，Hu Jun，et al.Monitoring of ground surface deformation in mining area with InSAR technique[J].The Chinese Journal of Nonferrous Metals,2011,21(10):2564-2576.

[12] 闫建伟,汪云甲,陈国良,等.钱营孜煤矿地表沉陷的D-InSAR监测[J].金属矿山，2011(3):105-107. Yan Jianwei，Wang Yunjia，Chen Guoliang，et al.Ground subsidence monitoring with D-InSAR in Qianyingzi Coal Mine[J].Metal Mine，2011(3):105-107.

[13] Goldstein R M，Zebker H A，Werner C L.Satellite radar interferometry：two-dimensional phase unwrapping[J].Radio Science，1988，23(4)：713-720.

[14] Buckland J R，Huntley J M，Turner S R E.Unwrapping noisy phase maps by use of a minimum-cost-matching algorithm[J].Applied Optics，1995，34(23)：5100-5108.

[15] 武继峰,杨志强.InSAR技术应用于矿区开采沉陷监测的数据优选研究[J].测绘通报，2013(S):92-94. Wu Jifeng，Yang Zhiqiang.Data optimization of mining subsidence monitoring by InSAR technology[J].Bulletin of Surveying and Mapping，2013(S):92-94.

[16] 姚丹丹,吴 侃,何 强.基于D-InSAR技术的煤矿区沉陷监测[J].金属矿山，2014(11):151-155. Yao Dandan，Wu Kan，He Qiang.Coal mine subsidence monitoring based on D-InSAR technology[J].Metal Mine,2014(11):151-155.

[17] 李培现.深部开采地表沉陷规律及预测方法研究[D].徐州：中国矿业大学,2012. Li Peixian.Study on Regularity and Prediction Method of Surface Subsidence Due to Deep Mining[D].Xuzhou:China University of Mining and Technology,2012.

[18] 陈炳乾,邓喀中,范洪冬.基于D-InSAR技术和SVR算法的开采沉陷监测与预计[J].中国矿业大学学报,2014(5):880-886. Chen Bingqian，Deng Kazhong，Fan Hongdong.Combining D-InSAR and SVR for monitoring and prediction of mining subsidence[J].Journal of China University of Mining and Technology,2014(5):880-886.

[19] GB/T7714—2011JTG/TD31 采空区公路设计与施工技术明细[S].北京：人民交通出版社，2011. GB/T7714—2011JTG/TD31 Guidelines for Design and Construction of Highway Engineering in the Mined-out Area[S].Beijing:China Communications Press,2011.

(责任编辑 石海林)

Monitoring and Prediction of Mining Subsidence based on D-InSAR and Probability Integral Method

Yang Junkai1Chen Bingqian2Deng Kazhong1He Qiang1Zhao Weiying1

(1.SchoolofEnvironmentalScienceandSpatialInformatics，ChinaUniversityofMining&Technology，Xuzhou221116，China;2.KeyLaboratoryofMineSpatialInformationTechnologies，NationalAdministrationofSurveying，MappingandGeoinformation，Jiaozuo454000，China)

In order to obtain the deformation law of high grade highway near mining area and overcome the limitations of traditional mining subsidence monitoring methods，a mining subsidence monitoring and predicting method based on the combination of the Differential SAR Interferometry (D-InSAR) technique and probability integral method was presented.The D-InSAR technology was used to monitor and analyze the subsidence trend of the 9308 working face in Yanzhou，and the monitor results were compared with leveling data.Before modifying the prediction parameters into sufficiency mining condition，the probability integral method parameters of the 9308 working face under insufficiency mining were calculated using D-InSAR monitoring results.Then the influence of high grade highway due to coal mining was predicted using the modified prediction parameters.Experimental results verified that D-InSAR technology had high monitoring accuracy which can provide reliable data for the prediction of mining subsidence，and this method could provide strong support to realize the integration of monitoring and prediction of mining subsidence.

D-InSAR，Probability integral method，Mining subsidence，Prediction parameters

2015-02-12

测绘地理信息公益性行业科研专项(编号：201412016),江苏高校优势学科建设工程项目(编号：SZBF2011-6-B35),矿山空间信息技术国家测绘地理信息局重点实验室开放基金项目(编号：KLM201311)。

杨俊凯(1990—)，男，硕士研究生。

TD325，TD17

A

1001-1250(2015)-04-195-06