煤矿采空区地表重力异常效应模拟研究_汤伏全.pdf

第 43 卷第 4 期煤炭学报Vol．43No．4 2018 年4 月 JOURNAL OF CHINA COAL SOCIETYApr．2018 汤伏全， 李庚新， 原一哲．煤矿采空区地表重力异常效应模拟研究［J］ ．煤炭学报， 2018， 43 4 945－950．doi 10．13225/j．cnki．jccs． 2017．1578 TANG Fuquan， LI Gengxin， YUAN Yizhe．Simulation and study on ground gravity anomaly effect in goaf of coal mine［ J］ ．Journal of China Coal Society， 2018， 43 4 945－950．doi 10．13225/j．cnki．jccs．2017．1578 煤矿采空区地表重力异常效应模拟研究 汤伏全， 李庚新， 原一哲 西安科技大学 测绘科学与技术学院， 陕西 西安710054 摘要 地下大规模采煤形成的空洞效应导致重力场发生变化， 通过地面重力变化来反演采空区上 方空洞特征及其演变对于煤矿安全生产和地面保护具有实际意义。长期以来， 煤矿地下采空导致 的地表重力异常效应因量级较小而被忽视。基于布格重力理论采用数值模拟方法得到采空区地表 的重力异常分布特征， 建立采空区深度、 宽度、 厚度及开采边界与地表重力异常分布的关系， 并使用 三阶紧致差分逼近算法得到采空区地表重力异常的一次差分曲线。通过加入随机误差模拟得到采 空区地表重力异常值， 采用不同拟合方法进行曲线平滑处理， 依据其一次差分曲线的峰值来推断采 空区的边界位置， 为现有技术条件下利用重力变化反演煤矿采空区演变提供了可行性。 关键词 煤矿; 采空区; 重力异常; 重力测量; 反演; 紧致差分格式 中图分类号 TD82文献标志码 A 文章编号 0253－9993 2018 04－0945－06 收稿日期 2017－11－12修回日期 2018－01－25 责任编辑 常明然 基金项目 国家自然科学基金资助项目 51674195 ; 陕西省自然科学基金资助项目 2016JM5048 作者简介 汤伏全 1966 ， 男， 湖南湘潭人， 教授， 博士。E－mail 2504557922 qq. com Simulation and study on ground gravity anomaly effect in goaf of coal mine TANG Fuquan， LI Gengxin， YUAN Yizhe College of Geomatics， Xi’ an University of Science and Technology， Xi’ an710054， China Abstract The goaf effect ed by the large- scale underground coal mining will change ground gravity field．There- fore， it is of practical significance to inverse the evolution and characteristics of variation of the hollow above the coal mine goaf by ground gravity survey for the safety of coal mine and ground protection．For a long time， due to its small magnitude， this kind of gravity anomaly effect caused by coal mine goaf is neglected．Based on the Bouguer gravity the- ory， the author uses numerical simulation to obtain ground gravity anomaly distribution characteristics of the goaf in longwall face， and establish the corresponding relationship between depth， width， thickness， mining boundary of working face and the distribution of gravity anomalies， use third- order compact difference scheme to acquire its first de- rivative curve．Different kinds of fitting s are adopted to smooth the distribution curve of the ground gravity anomaly in the goaf， which is simulated by adding random errors of different sizes．Also， according to the peak point of the first order difference curve， the boundary position of the goaf is deduced， which provides the feasibility of using the gravity to inverse the change of coal mine goaf under the existing technical conditions． Key words coal mine; coal mine goaf; gravity anomaly; gravity survey; inversion; compact difference scheme 我国煤炭资源丰富， 其开采方式主要为地下开 采。多年来的大规模高强度开采形成了大量的地下 采空区 ［1 ］， 在上覆岩层的自重作用下顶板将逐层垮 落导致空洞效应向上传递， 直至地表形成塌陷盆地。 ChaoXing 煤炭学报 2018 年第 43 卷 地表塌陷的程度取决于煤层开采深度、 工作面大小、 覆岩特性等多种因素。由于受岩层垮落碎胀影响和 基岩关键层的控制作用， 多数采空区上覆岩层并未充 分垮落和压实， 不同层位上 包括采空区 存在大量 的空洞， 尤其是一些建筑物下的条带开采和地方小煤 窑， 以及在巨厚坚硬岩层下进行采煤时， 采空区的垮 落基本上被上覆岩层所阻隔， 并未传递到地表。一些 小煤窑的盗采和滥采使得老采空区的空间位置亦无 法确定。这种长期存在的地下空洞效应导致老采空 区地面发生持续的残余变形或突然塌陷等安全隐患， 严重影响了矿区的安全生产和各种工程的选址与建 设 ［2 ］。目前， 对于老采空区的探测一般采用钻探方 法或物探方法。常用物探技术有电阻率法、 探地雷达 法、 电磁法、 三维地震法、 放射性法、 重力法等 ［3－8 ］。 其中， 重力法是通过重力场变化判定和解读地球内部 结构构造、 矿产资源分布规律等的一种重要手段。老 采空区的空洞效应会引起地面重力场发生变化， 但由 于这种重力异常较小， 一般只有 50 ～ 500 μGal 的量 级。因此多数学者认为， 通过地表获得的重力异常数 据来反演地下采空区的空间分布和变化特征在技术 上难以实现。而近年来， 随着重力测量仪的精度不断 提高， 如 CG － 5 型相对重力仪的测量随机误差在 5 μGal 左右， 展开上述问题研究已具有现实可能性。 为此， 本文以现有的布格重力理论模型为基础， 将地 下采空区简化为立方体空洞， 分析采空区地表的重力 异常与采空区特征之间的关系， 通过数值模拟方法探 讨基于重力异常效应进行采空区反演的可行性， 为煤 矿采空区探测和地面长期稳定性监测提供一种新的 技术途径。 1煤矿采空区地表重力异常计算模型 地下岩层、 矿产密度分布不均匀或因矿产开发等 质量亏损会导致地表重力发生变化， 称为重力异 常 ［9 ］。本文采用布格重力模型来模拟计算规则采空 区地表的重力变化 ［10 ］。如图 1 所示， 以采空区左边 界中央在地表的垂直投影作为地表坐标系原点， 地水 平面内沿采空区延伸方向为 X 轴， 沿采空区宽度方 向为 Y 轴， 垂直向下为 Z 轴正向。采空区长、 宽、 高 分别 A， B， C， 埋深为 H。 根据重力场的计算公式， 该采空区的重力场可表 示为 Δg x， y， z Gρ∫ A 0∫ B 2 － B 2 ∫ －H －H－C ζ － z R3 dξdηdζ 1 式中， G 为万有引力常数; ρ 为场源密度; x， y， z 为 观测点的坐标; ξ， η， ζ 为采空区场源点的坐标。 图 1计算模型坐标系 Fig. 1Coordinate system of calculation R { ξ－ x 2 η－ y 2 ζ－ z 2} 1 2 对式 1 的积分求解， 许多学者都进行了研 究 ［11 ］， 本文使用式 2 来计算该积分。 Δg － Gρ‖ ξ－ x ln［ η－ y R］ η－ y ln［ ξ－ x R］－ ζ－ z arctan ξ－ x η － y ζ－ z R A 0 B 2 － B 2 －H －H－C 2 2采空区地表的重力异常分布特征 为分析采空区特征对于地表重力异常的影响， 以 矩形工作面开采条件建立模拟计算的基准模型， 设采 空区长度 A、 宽度 B、 高度 C 及埋深 H 分别为 2 000， 200， 4， 300 m。剩余密度为 ρ 2 680 g/cm 3， 在其它 尺度参数不变的情况下， 地表重力异常随采空区长度 A 变化的分布曲线如图 2 图中地表点位置坐标零点 在采空区左边界正上方， 向右为正 所示。当采空区 长度达到 1 000 m 后， 重力异常峰值大小趋于不变， 在采空区中央地表达到最大重力异常值， 在该基准模 型条件下约为90 μGal。重力异常曲线呈现以采空区 中心为对称的分布形态， 这对于通过重力异常分布来 反演采空区边界具有参考意义。 图 2不同采空区长度 A 的重力异常分布曲线 Fig. 2Distribution curves of surface gravity anomalies at different length 当采空区深度 H 从 300 m 变化到 100 m 时， 地 表重力异常分布曲线如图 3 所示。各曲线变化形 态基本相似， 在采空区边界附近的地表有所不同。 随着空洞深度变小， 重力异常峰值由 59 μGal 增大 到 156 μGal， 其发生位置由采空区中央向边界方向 649 ChaoXing 第 4 期汤伏全等 煤矿采空区地表重力异常效应模拟研究 图 3不同采深 H 的地表重力异常分布曲线 Fig. 3Distribution curves of surface gravity anomalies at different mining depth 偏移。 当采空区宽度 B 从100 m 变化到300 m 时， 地表 重力异常分布曲线如图 4 所示， 随着开采宽度的增 大， 重力异常峰值由 49 μGal 增大到 122 μGal， 各曲 线分布形态基本呈一致。 图 4不同采宽 B 的地表重力异常分布曲线 Fig. 4Distribution curves of surface gravity anomalies at different mining width 当采空区高度 C 从3 m 变化到7 m 时， 地表重力 异常分布曲线如图 5 所示。重力异常峰值由 68 μGal 增大到 158 μGal， 其发生位置位于采空区内侧地表， 各曲线分布形态基本呈一致。 图 5不同采高 C 的地表重力异常分布曲线 Fig. 5Distribution curves of surface gravity anomalies at different mining depths 分析图 4 和图 5 可知， 在基准模型条件下地表重 力异常峰值与采宽、 采高均呈线性正相关。分析图 3 可知， 地表重力异常峰值与采空区埋深呈非线性的负 相关， 如图 6 所示。 图 6采空区埋深与重力异常峰值的关系 Fig. 6Relationship between different mining depths and maximum gravity anomalies 在基准模型条件下， 采空区重力异常分布曲 线 图 7 a 的形态较为简单， 在采空区边界位置的 地表重力异常未见明显特征。所以， 仅依据重力异常 曲线难以反演确定重力变化与采空区边位置的关系。 为此， 对重力异常曲线 图 7 a 进行一次差分处理 来分析其变化特征 ［12－13 ］。由于传统的有限差分格式 在步长较大的情况下精度偏低， 若要提高精度则需要 增加网格点。因此， 这里采用计算物理学中广泛应用 的紧致差分格式 ［14－15 ］逼近算法来进行一次差分处 理。该方法所需的网格点较少， 且截断误差系数较 小， 其基本算法如下 对于区间［ a， b］上的离散函数 F x ， 节点为 x1， x2， ， xn， 步长 h xi －x i－1， i 1， 2， ， n。节点处的函数值为 f i ， 一阶导数值为 df i ， 对于节点 xi， 设一阶导数与函数值之间满足式 3 ∑ l≥0 al df i l df i － l 2 ∑ l ＞ 0 al f i l f i － l 2l ∑ l≥0 al∑ l ＞ 0 al 1 3 式中， al为差分格式的待定系数， l1， 2， 3， ， n。 图 7基准模型下地表重力异常及其一阶导数曲线 Fig. 7Ground gravity anomaly and its first derivative curve under the normal model 749 ChaoXing 煤炭学报 2018 年第 43 卷 则函数 F x 具有三阶逼近精度的一次差分。即 其一阶导数 df i 可通过对式 3 采用 Taylor 级数展 开， 根据精度需要得到待定系数 al的代数方程， 用追 赶法求解该代数方程得到相应的待定系数。本文采 用计算步长为 h10 m， 并对边界点使用 5 点偏差逼 近 i1 和 i n 时 ， 求得三阶紧致差分逼近算法满 足 df i － 1 12［ 25f i － 48f i 1 36f i 2 － 16f i 3 3f i 4 ］ ， i 1 3 2 t i － 1 3 df i － 1 [] ， i ∈［ 2 ∶ n － 1］ 1 12［ 25f i － 48f i － 1 36f i － 2 － 16f i － 3 3f i － 4 ］ ， i n              4 式中， t i 5 6 ［ f i － f i － 1 ］ 1 6 ［ f i 1 － f i ］ ， i ∈［ 2 ∶ n － 1］ ， n b － a h ， h 为步长。 对图 7 a 的曲线函数值利用上述算法计算其一 次差分， 得到基准模型下地表重力异常分布的一次差 分曲线， 如图 7 b 所示。 由图 7 b 可以看出， 重力异常一次差分曲线的 两个峰值正好位于采空区边界正上方附近。两峰值 之间的距离与采空区长度相等。因此在理论模型下， 可根据采空区地表重力异常的一次差分曲线的峰值 位置来判断采空区边界位置并确定采空区范围。 3基于重力异常模拟数据反演采空区边界位 置 由于实地观测的重力异常数据中不可避免会存 在一定的随机误差， 使得重力异常曲线偏离图 7 所示 的标准形态。在大地测量及岩土工程领域， 常通过数 值模拟和加入随机误差的方法来研究实际工程问 题 ［16 ］。因此， 本文在重力异常曲线中加入一定的随 机误差后， 来分析基于重力异常反演采空区边界位置 的可行性。 首先， 根据目前高精度相对重力仪 如 CG－5 型 仪器 的测量精度情况， 在理论重力异常曲线中加入 5 μGal 的随机测量误差， 得到地表重力 g 异常模拟 测量值的散点分布， 分别采用有理数拟合法、 正弦曲 线逼近法、 平滑样条拟合法和傅立叶变换拟合法对模 拟测量值进行拟合和平滑处理， 得到拟合曲线及其残 差分布， 如图 8 所示。将图 8 的拟合平滑曲线利用三 阶紧致差分逼近算法， 计算其一阶导数得到地表重力 异常的一阶差分曲线， 如图 9 所示。 为了衡量上述拟合方法的优劣程度， 采用确定系 数 R－square 和回归拟合标准差 RMSE 来评价拟 合优度。R－square 为预测数据与原始数据均值之差 的平方和除以原始数据与均值之差的平方和。R－ square 的变化范围为［ 0， 1］ ， 其值越大表明拟合程度 越好。RMSE 为预测数据和原始数据对应点误差的 平方和之均值的平方根。上述 4 种拟合方法的拟合 优度系数见表 1。 图 8不同方法拟合的重力异常模拟曲线 Fig. 8Fitting curves and residuals of gravity anomaly values by different s 849 ChaoXing 第 4 期汤伏全等 煤矿采空区地表重力异常效应模拟研究 图 9加5 μGal 随机误差的重力异常一阶差分曲线 Fig. 9First order difference curve of gravity anomaly with 5 μGal random error 表 1 4 种拟合方法的拟合优度参数 Table 1Goodness of fit parameters for four fitting s 拟合参数 拟合方法 有理数 拟合法 正弦曲线 逼近法 平滑样条 拟合法 傅里叶 变换法 R－square0. 982 20. 992 70. 995 10. 992 7 RMSE5. 463 33. 526 12. 937 03. 517 7 对图 8 和表 1 的拟合效果参数分析可得 在加入 5 μGal 随机误差情况下， 采用平滑样条拟合方法得 到最佳拟合效果， 傅里叶变换法次之， 正弦曲线逼近 法和有理数拟合法法拟合效果最差。对比图 9 和图 7 可知， 平滑样条拟合法， 有理数拟合法和正弦曲线 逼近法的一次差分曲线与没有随机误差的标准差分 曲线相比， 其峰值位置相差 10 m。即在5 μGal 的随 机误差影响下， 若采用这 3 种方法获取一次差分曲 线， 其反演的采空区位置偏离正确位置 10 m。而傅 里叶变换法的一次差分曲线， 其分布特征和峰值位置 与标准曲线 图 7 一致。 将随机误差增大到10 μGal， 利用傅里叶变换法 对上述重力异常模拟数值进行拟合后并得到其一阶 差分曲线， 如图 10 所示。由图 10 可知， 其一阶差分 曲线的峰值位置与采空区边界位置相一致。 为探讨随机误差对于采空区边界位置反演效果 的影响， 将随机误差不断增大到50 μGal， 通过上述 方法得到空洞边界反演的最大位置偏差与随机测量 误差的关系曲线， 如图 11 所示。图 11 中当测量随机 误差为15 uGal 时， 所反演的空洞边界位置偏差达 10 m， 反演效果变差。 图 10加入10 μGal 随机误差后重力异常一阶差分曲线 Fig. 10First order difference curve of gravity anomaly with 10 μGal random error 图 11最大反演偏差与测量随机误差的关系曲线 Fig. 11Relationship between the maximum inversion and measurement random error 4结论 1 采空区地表的重力异常与采宽、 采高呈线性 正相关， 与开采深度呈非线性的负相关。基于理论模 949 ChaoXing 煤炭学报 2018 年第 43 卷 型的重力异常一阶差分曲线的峰值位置与采空区边 界位置相一致。 2 对于随机误差不超过10 μGal 的模拟重力 异常数据， 使用傅里叶拟合方法进行平滑并利用三阶 紧致差分逼近得到其一阶差分曲线， 可通过该曲线峰 值位置来推断地下空洞的边界位置。因此， 在一定的 测量误差情况下， 利用地表重力异常数据反演地下空 洞仍然具有可行性。 应该指出， 在实际情形中当地下空洞呈非规则形 状和离散分布时， 也可利用数值叠加方法得到地表重 力异常值。因此， 在一定的边界条件和简化假设下， 通过重力异常数据可反演得到诸如条带开采或小煤 窑老采空区等空洞在埋深、 边界、 大小等方面特征的 演变过程。 参考文献 References ［ 1］杜坤， 李夕兵， 刘伟科， 等．采空区危险性评价的综合方法及工程 应用［J］．中南大学学报 自然科学版 ， 2011， 42 9 2803－ 2804． DU Kun， LI Xibing， LIU Kewei， et al．Comprehensive uation of underground goaf risk and engineering application［J］ ． Journal of Central South University Science and Technoogy ， 2011， 42 9 2803－2804． ［ 2］魏晓刚， 麻凤海， 刘书贤， 等．煤矿采空区地震安全防护的若干问 题［J］．地震研究， 2016， 39 1 151－152． WEI Xiaogang， MA Fenghai， LIU Shuxian， et al．Problems of safe- ty control of seismic dynamic disasters in doal mine goaf［ J］．Journal of Seismological Research， 2016， 39 1 151－152． ［ 3］刘菁华， 王祝文， 朱士， 等．煤矿采空区及塌陷区的地球物理探查 ［J］．煤炭学报， 2005， 30 6 716－717． LIU Jinghua， WANG Zhuwen， ZHU Shi， et al．The geophysical explo- ration about exhausted area and singking area in coal mine［J］．Jour- nal of China Coal Society， 2005， 30 6 716－717． ［ 4］程久龙， 潘冬明， 李伟， 等．强电磁干扰区灾害性采空区探地雷达 精细探测研究［J］．煤炭学报， 2010， 35 2 227－228． CHENG Jiulong， PAN Dongming， LI Wei， et al．Study on the detec- ting of hazard abandoned workings by ground penetrating radar on strong electromagnetic interference area ［J］ ． Journal of China Coal Society， 2010， 35 2 227－228． ［ 5］覃思， 程建远， 胡继武， 等．煤矿采空区及巷道的井地联合地震超 前勘探［J］．煤炭学报， 2015， 40 3 636－637． QIN Si， CHENG Jianyuan， HU Jiwu， et al．Coal- seam- ground- seismic for advance detection of goaf and roadway［J］ ． Journal of China Coal Society， 2015， 40 3 636－637． ［ 6］姜志海， 杨光．浅埋特厚煤层小窑采空区瞬变电磁探测技术研究 及应用［J］．采矿与安全工程学报， 2014， 31 5 770－771． JIANG Zhihai， YANG Guang．Research and application of TEM de- tection technology for small mine gob area in shallowly- buried and extremely thick seam［J］．Journal of Mining ＆ Safety Engineering， 2014， 31 5 770－771． ［ 7］李文．煤矿采空区地面综合物探方法优化研究［J］．煤炭科学技 术 2017， 45 1 194－195． LI Wen．Optimization study of surface comprehensive geophysical de- tection s of coal mine goafs［J］．Coal Science and Technology 2017， 45 1 194－195． ［ 8］付天光．综合物探方法探测煤矿采空区及积水区技术研究［J］． 煤炭科学技术， 2014， 42 8 90－91． FU Tianguang． Study on technology of comprehensive geophys- ic exploration of mine goaf and water accumulated area［J］． Coal Science and Technology， 2014， 42 8 90－91． ［ 9］吴亮， 刘长弘， 王庆宾， 等．重力异常对地下异常体边界的识别算 法［J］．测绘科学， 2015， 40 12 34－35． WU Liang， LIU Changhong， WANG Qingbin， et al． An identifica- tion of underground anomaly boundary using gravity anomaly ［J］．Science of Surveying and Mapping， 2015， 40 12 34－35． ［ 10］SU Yongjun， CHENG Lizhen， MICHEL Chouteau， et al．New im- proved ulas for calculating gravity anomalies based on a cylin- der model［J］．Journal of Applied Geophysics， 2012， 86 2 36－ 43． ［ 11］骆遥．两种新的长方体重力场正演表达式及其理论推导［J］．工 程地球物理学报， 2008， 5 2 211－212． LUO Yao．New Expressions for gravitational attraction of a homoge- neous rectangular prism and its deduction［J］．Chinese Journal of Engineering Geophysics， 2008， 5 2 211－212． ［ 12］ZUO Boxin， HU Xiangyun．Edege detection of gravity field using ei- genvalue analysis of gravity gradient tensor［J］． Journal of Ap- plied Geophysics， 2015， 114 1 263－270． ［ 13］Wu Heyu， Li Lu， Xing Congcong， et al．A new of edge de- tection based on the total horizontal derivative and the modulus of full tensor gravity gradient［J］． Journal of Applied Geophysics， 2017， 139 1 239－245． ［ 14］孙建安， 贾伟， 吴广智．一种非均匀网格上的高精度紧致差分格 式［J］．西北师范大学学报 自然科学版 ， 2014， 50 4 31－33． SUN Jian’an， JIA Wei， WU Guangzhi． A higher accurate com- pact difference scheme on non- uni grid［J］．Jounral of North- west Normal University Natural Science ， 2014， 50 4 31－33． ［ 15］李佳， 罗纪生．一阶、 二阶导数耦合的紧致差分格式及其应用 ［J］．计算机工程与应用， 2012， 48 1 13－15． LI Jia， LUO Jisheng． Derivation and application compact scheme of coupling of first and second derivative［J］．Computer Engineering and Applications， 2012， 48 1 13－15． ［ 16］刘惠敏， 刘繁明， 荆心， 等．浅地表地下质量异常重力与磁探测 数值模拟方法［J］．中国惯性技术学报， 2016， 24 4 491－492． LIU Huimin， LIU Fanming， JING Xin， et al． Numerical simula- tion for gravity and magnetic observations of mass anomaly in shallow subsurface［J］．Journal of Chinese Inertial Technology， 2016， 24 4 491－492． 059 ChaoXing