康 東，李庭筠，劉 偉，白成武，楊小峰，張麗娟

（1.成都同創(chuàng)眾益科技有限公司 四川成都 610000；2.成都理工大學(xué) 四川成都 610000；3.四川煤礦安全監(jiān)察局技術(shù)中心 四川成都610000）

0 引 言

隨著經(jīng)濟的快速發(fā)展，我國對能源的需求量不斷擴大。我國是煤炭生產(chǎn)大國，煤炭的大規(guī)模開發(fā)為國民經(jīng)濟的發(fā)展提供了保證，然而多年來對煤炭的頻繁開采在礦區(qū)產(chǎn)生了許多人為挖掘的地下空洞（即采空區(qū)），采空區(qū)突水是導(dǎo)致礦井突水等災(zāi)害的主要因素之一[1]，采空區(qū)帶來了諸多安全隱患，影響了礦業(yè)工程建設(shè)的進度和質(zhì)量。

為了保障煤礦生產(chǎn)的安全與效率，解決采空區(qū)帶來的安全隱患十分關(guān)鍵。目前國內(nèi)外常用的探測采空區(qū)和地下水的方法有井下物探、化學(xué)探測、鉆探方法[2]，其中物探方法主要包括高密度電法、瞬變電磁法、地質(zhì)雷達法、微重力法、地震法、大地電磁法等[3]。但是傳統(tǒng)的方法存在著探測距離小、反演精度不高的缺點[4]。超導(dǎo)量子干涉器（Superconducting Quantum Interference Device，簡稱SQUID）高精度磁測儀的發(fā)展使地下空洞的探測距離、探測精度和反演精度得到極大提高[5]，到目前為止，LTS-SQUID（Low Tc Superconductor-SQUID）的靈敏度可達 1fT，磁梯度測量方面，CSIRO和 DSTO研發(fā)的高溫超導(dǎo)旋轉(zhuǎn)梯度儀靈敏度可達 0.05nT/m，未來靈敏度還將達到1pT/m[6]?；?SQUID 的磁梯度儀[7]和全張量磁梯度儀可以更加有效地反映出磁異常源特征[8]，并通過反演計算得到地下空洞和水源的位置、深度等特征參數(shù)，有效發(fā)現(xiàn)地下空洞存在導(dǎo)致的磁場異常，實現(xiàn)地下空洞、地下水的探測和定位。

磁法測量包括磁場總場強度BT、總場梯度（總場空間變化率）、總場三分量及磁梯度張量（總場分量空間變化率）。本文基于 COMSOL軟件對采空區(qū)建立數(shù)學(xué)物理模型進行仿真計算，探討地下空洞對該區(qū)域磁場的影響，對不同磁測量量進行對比分析。

1 計算模型

假設(shè)空間尺寸為 10km×10km×7km，建立幾何模型。設(shè)定 2km厚度的地層和 5km厚度的大氣層，假設(shè)地層和大氣層為均勻介質(zhì)。在空間的水平面中心處放置 1個三軸為 500m×50m×50m的橢球體，用來模擬人工挖掘產(chǎn)生的地下空洞，如圖 1所示?？斩吹难娱L方向選擇東西向（x方向），可以更好地表現(xiàn)東西向磁強分量和梯度的變化。

圖1 模型的幾何構(gòu)造（y方向為正北方）Fig.1 The geometric construction of the modely showing due north

假設(shè)空間處于均勻磁場，地表面的海拔為 0m。設(shè)地磁場磁傾角 Incl＝48.067°，磁偏角 Decl＝-2.017°，地磁場強度為 50763.7nT，地磁場隨著緯度的變化存在細微的變化，約為 6.3×10-11T/m。地磁的強度和方向參考北緯 30°40'12"，東經(jīng) 104°4'16"的地磁要素（模型坐標軸原點），數(shù)據(jù)來自美國國家地理數(shù)據(jù)中心（National Geophysical Data Center of the U.S.Government）[9]。僅考慮空洞中只存在空氣，所以設(shè)空洞和大氣層的相對磁導(dǎo)率μr1＝1。地層含有巖層或礦物，假設(shè)巖石和礦物均勻分布，含量占土層的1%，平均相對磁導(dǎo)率μr2＝3.5，地層不含磁性物質(zhì)的相對磁導(dǎo)率μr0＝1，則土層的相對磁導(dǎo)率為μr＝μr0×99%+μr2×1%＝1.025。假設(shè)巖石和礦物存在剩余磁場，且剩余磁場方向與地磁場方向相同，平均剩余磁場強度 Mr＝10A/m，則巖石和礦物的平均剩余磁通密度 Br＝Mr×μ0，地層的剩余磁通密度為 Br0＝Br×1%＝125.7nT。

2 探測對象及原理

對搭建的模型進行多種磁強信息的測量，包含磁場強度總量 BT，總量的 3個梯度 BTx、BTy、BTz，磁場強度分量 Bx、By、Bz，分量和總場的關(guān)系為：

對于磁場強度三分量在x、y、z方向的梯度，即梯度張量，公式為[10]：

根據(jù)梯度張量的公式，可以推導(dǎo)出只有5個要素相互獨立，主對稱軸相關(guān)性為[11]：

將地下空洞假設(shè)為由若干磁偶極子構(gòu)成的磁異常源，異常源的磁通密度B可由2個分量來表示[11]：

式中：?為磁偶極子指向觀測點的位置矢量，θ為磁矩?與?的夾角，μ為磁導(dǎo)率。

3 模擬計算與分析

以地下空洞頂端與模型地表的距離d為變量，對選取的物理量進行測量。選取地層上表面與大氣層下表面及地面作為測量平面（如圖 2），測量該表面磁場強度及梯度，將測量信息以二維圖像呈現(xiàn)，磁場強度大小用灰階表示。選擇 d1＝0.01m，d2＝100m，d3＝500m和d4＝1000m分析。

建立3條平行于y軸的地表測線，用于模擬地面磁力儀的測量，如圖 3所示。從左到右依次為測線l1（x＝4750m）、測線 l2（x＝5000m）和測線 l3（x＝5250m）。

在 d1＝0.01m 條件下，地表測量的磁場總量、分量和梯度均能很好地呈現(xiàn)，此時空洞接近地表，圖像可以看作磁異常源與背景場大小的對比。同時，圖像很好地顯示了不同磁場強度量在地表的分布情況。磁場△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）的方式表現(xiàn)了式（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 結(jié) 論

準確地發(fā)現(xiàn)并且定位采空區(qū)是排除礦區(qū)突水安全隱患的一個重要因素，本文根據(jù)建模仿真計算，對比了磁法測量中多個測量對象在磁異常源不同深度下的強度，顯示出總磁場強度及其分量探測采空區(qū)的局限性。預(yù)測總場梯度和梯度張量在采空區(qū)探測距離上有潛在的優(yōu)勢，梯度信號削弱了背景磁場的影響，增強了異常源變化的幅值，但在 1000m 深度下異常源只有 pT級的強度，需要使用極高精度的磁梯度儀（如 SQUID）配合測量。這對未來在礦區(qū)實地測量具有指導(dǎo)意義。

礦區(qū)的磁場環(huán)境具有多樣性和不確定性，為進一步研究磁梯度在采空區(qū)勘探應(yīng)用中的可行性和可靠性，還需改變模型中地下空洞和地層的性質(zhì)參數(shù)，做更多的驗證性仿真計算。