基于采空區深度變化磁梯度測量磁場信息的仿真

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

（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）配合測量。這對未來在礦區實地測量具有指導意義。

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