基于采空区深度变化磁梯度测量磁场信息的仿真

康 东，李庭筠，刘 伟，白成武，杨小峰，张丽娟

（1.成都同创众益科技有限公司 四川成都 610000；2.成都理工大学 四川成都 610000；3.四川煤矿安全监察局技术中心 四川成都610000）

0 引 言

随着经济的快速发展，我国对能源的需求量不断扩大。我国是煤炭生产大国，煤炭的大规模开发为国民经济的发展提供了保证，然而多年来对煤炭的频繁开采在矿区产生了许多人为挖掘的地下空洞（即采空区），采空区突水是导致矿井突水等灾害的主要因素之一[1]，采空区带来了诸多安全隐患，影响了矿业工程建设的进度和质量。

为了保障煤矿生产的安全与效率，解决采空区带来的安全隐患十分关键。目前国内外常用的探测采空区和地下水的方法有井下物探、化学探测、钻探方法[2]，其中物探方法主要包括高密度电法、瞬变电磁法、地质雷达法、微重力法、地震法、大地电磁法等[3]。但是传统的方法存在着探测距离小、反演精度不高的缺点[4]。超导量子干涉器（Superconducting Quantum Interference Device，简称SQUID）高精度磁测仪的发展使地下空洞的探测距离、探测精度和反演精度得到极大提高[5]，到目前为止，LTS-SQUID（Low Tc Superconductor-SQUID）的灵敏度可达 1fT，磁梯度测量方面，CSIRO和 DSTO研发的高温超导旋转梯度仪灵敏度可达 0.05nT/m，未来灵敏度还将达到1pT/m[6]。基于 SQUID 的磁梯度仪[7]和全张量磁梯度仪可以更加有效地反映出磁异常源特征[8]，并通过反演计算得到地下空洞和水源的位置、深度等特征参数，有效发现地下空洞存在导致的磁场异常，实现地下空洞、地下水的探测和定位。

磁法测量包括磁场总场强度BT、总场梯度（总场空间变化率）、总场三分量及磁梯度张量（总场分量空间变化率）。本文基于 COMSOL软件对采空区建立数学物理模型进行仿真计算，探讨地下空洞对该区域磁场的影响，对不同磁测量量进行对比分析。

1 计算模型

假设空间尺寸为 10km×10km×7km，建立几何模型。设定 2km厚度的地层和 5km厚度的大气层，假设地层和大气层为均匀介质。在空间的水平面中心处放置 1个三轴为 500m×50m×50m的椭球体，用来模拟人工挖掘产生的地下空洞，如图 1所示。空洞的延长方向选择东西向（x方向），可以更好地表现东西向磁强分量和梯度的变化。

图1 模型的几何构造（y方向为正北方）Fig.1 The geometric construction of the modely showing due north

假设空间处于均匀磁场，地表面的海拔为 0m。设地磁场磁倾角 Incl＝48.067°，磁偏角 Decl＝-2.017°，地磁场强度为 50763.7nT，地磁场随着纬度的变化存在细微的变化，约为 6.3×10-11T/m。地磁的强度和方向参考北纬 30°40'12"，东经 104°4'16"的地磁要素（模型坐标轴原点），数据来自美国国家地理数据中心（National Geophysical Data Center of the U.S.Government）[9]。仅考虑空洞中只存在空气，所以设空洞和大气层的相对磁导率μr1＝1。地层含有岩层或矿物，假设岩石和矿物均匀分布，含量占土层的1%，平均相对磁导率μr2＝3.5，地层不含磁性物质的相对磁导率μr0＝1，则土层的相对磁导率为μr＝μr0×99%+μr2×1%＝1.025。假设岩石和矿物存在剩余磁场，且剩余磁场方向与地磁场方向相同，平均剩余磁场强度 Mr＝10A/m，则岩石和矿物的平均剩余磁通密度 Br＝Mr×μ0，地层的剩余磁通密度为 Br0＝Br×1%＝125.7nT。

2 探测对象及原理

对搭建的模型进行多种磁强信息的测量，包含磁场强度总量 BT，总量的 3个梯度 BTx、BTy、BTz，磁场强度分量 Bx、By、Bz，分量和总场的关系为：

对于磁场强度三分量在x、y、z方向的梯度，即梯度张量，公式为[10]：

根据梯度张量的公式，可以推导出只有5个要素相互独立，主对称轴相关性为[11]：

将地下空洞假设为由若干磁偶极子构成的磁异常源，异常源的磁通密度B可由2个分量来表示[11]：

式中：?为磁偶极子指向观测点的位置矢量，θ为磁矩?与?的夹角，μ为磁导率。

3 模拟计算与分析

以地下空洞顶端与模型地表的距离d为变量，对选取的物理量进行测量。选取地层上表面与大气层下表面及地面作为测量平面（如图 2），测量该表面磁场强度及梯度，将测量信息以二维图像呈现，磁场强度大小用灰阶表示。选择 d1＝0.01m，d2＝100m，d3＝500m和d4＝1000m分析。

建立3条平行于y轴的地表测线，用于模拟地面磁力仪的测量，如图 3所示。从左到右依次为测线l1（x＝4750m）、测线 l2（x＝5000m）和测线 l3（x＝5250m）。

在 d1＝0.01m 条件下，地表测量的磁场总量、分量和梯度均能很好地呈现，此时空洞接近地表，图像可以看作磁异常源与背景场大小的对比。同时，图像很好地显示了不同磁场强度量在地表的分布情况。磁场△BT、△BTx、△BTy、△BTz的场强异常均在 10-7T量级，总场梯度和分量梯度的磁异常在 10-8～10-9T/m 量级。在 d2＝100m 下，随着磁异常源深度的增加，在测量平面的磁异常强度减弱，由式（1）和式（2）推断，随着空洞深度增加磁异常在地表的面积扩大，如图4所示。

图2 测量平面，箭头指向磁场强度方向Fig.2 Measurement plane arrows pointing to the direction of magnetic field intensity

图3 测线示意图Fig.3 From left to right in the schematic diagram of survey line showing the following—Survey Line l1（x＝4750m），Survey Line l2（x＝5000m）and Survey Line l3（x＝5250m）

图4 d2＝100m下磁场强度总量及分量（Y方向为正北）Fig.4 d2＝100m total and component of magnetic field intensity

磁场强度存在10-9T磁异常量，总场梯度和分量梯度均有明显的辨识度，信号幅值在 10-11～10-12T/m，可以有效地反映磁异常源信息。

d4＝1000m 时，BT的信号几乎与地磁场一致，总场梯度的 3个要素异常均在 10-12T/m量级，只有BTx在空间上的变化较为明显，磁强分量异常水平在10-10T，但只有 x分量在空间上具有较好的辨识度。通过图 5看出，对比梯度张量 9要素，以图像（见图 6）的方式表现了式（3），Bxy＝Byx，Bxz＝Bzx，Bzy＝Byz，梯度张量各要素除 Byy和 Bzz，在空间上均具有较好的辨识度，但磁异常幅值量级十分微小，介于0.5～10pT/m之间，接近磁梯度仪灵敏度的极限。

图5 d4＝1000m下磁场强度总量及分量Fig.5 d4＝1000m total and component of magnetic field intensity

图6 d4＝1000m下梯度张量Fig.6 d4＝1000m gradient tensor

4 结 论

准确地发现并且定位采空区是排除矿区突水安全隐患的一个重要因素，本文根据建模仿真计算，对比了磁法测量中多个测量对象在磁异常源不同深度下的强度，显示出总磁场强度及其分量探测采空区的局限性。预测总场梯度和梯度张量在采空区探测距离上有潜在的优势，梯度信号削弱了背景磁场的影响，增强了异常源变化的幅值，但在 1000m 深度下异常源只有 pT级的强度，需要使用极高精度的磁梯度仪（如 SQUID）配合测量。这对未来在矿区实地测量具有指导意义。

矿区的磁场环境具有多样性和不确定性，为进一步研究磁梯度在采空区勘探应用中的可行性和可靠性，还需改变模型中地下空洞和地层的性质参数，做更多的验证性仿真计算。