基于電磁場合矢量的陷落柱三維模擬及其探測

易洪春， 袁永榜， 牟義， 常江浩

(1.瓦斯災害監控與應急技術國家重點實驗室， 重慶 400039； 2.中煤科工集團重慶研究院有限公司， 重慶 400037； 3.煤炭科學技術研究院有限公司安全分院， 北京 100013； 4.河北地質大學地球科學學院， 石家莊 050031)

陷落柱是一種埋藏型巖溶在動力作用下發生坍塌和剝落后形成的錐狀塌陷體。陷落柱除自身常儲聚大量地下水外，還常成為溝通其他水源的導水通道，故陷落柱的存在容易造成礦井突水事故。突水性強、規模大、破壞嚴重一直是陷落柱突水的突出特點。

瞬變電磁法(transient electromagnetic method，TEM)是利用不接地回線或接地線源向地下發射一次脈沖磁場，在一次脈沖磁場間歇期間利用線圈或接地電極觀測地下介質中引起的二次感應渦流場，從而探測介質電阻率的一種方法，近年來被廣泛應用到煤礦井下水害探測中，取得了大量精準有效的探測成果。而瞬變電磁場在含水陷落柱中的傳播機理及響應特征一直是利用瞬變電磁法探測含水陷落柱需要弄清的問題。時域有限差分(finite-difference time-domain，FDTD)算法已被許多學者借鑒到電磁場的數值計算問題上來[1-3]，它不僅可以研究瞬變電磁場的產生機理，擴散過程及響應特征，同時也為瞬變電磁法探測含水低阻體提供了大量理論支撐。岳建華等[4]研究了層狀介質中地下瞬變電磁場的全空間效應；楊海燕等[5]對礦用多匝小回線源瞬變電磁場的分布規律進行了詳細分析；李建慧等[6]對基于電場方程的回線源產生的瞬變電磁場進行了三維正演研究；孟慶鑫等[7]對井中磁源瞬變電磁場進行了三維時域有限差分數值模擬；席振銖等[8]首次對瞬變電磁法三分量測量方法進行了研究；常江浩等[9]對礦井全空間磁源瞬變電磁場的擴散規律進行了詳細研究；程久龍等[10]對礦井全空間三維主軸各向異性介質瞬變電磁場響應特征進行了深入研究，解釋了探測結果存在偏差的原因。

目前，關于瞬變電磁場的研究主要集中于三分量中的某一分量，而對某兩個或三個分量的合矢量的研究相對較少。鑒于此，在前人研究的基礎之上，現采用基于非均勻網格的時域有限差分算法，研究全空間瞬變電磁場的擴散，并首次基于電場合矢量，詳細分析電場的傳播規律及響應特征，為瞬變電磁法超前探測掘進頭前方的含水區提供理論支撐，并為用合適量來進行瞬變電磁數據處理解釋奠定基礎。最后通過數值模型算例以及實際應用驗證瞬變電磁法探測含水陷落柱的有效性和適用性。

1 瞬變電磁場的FDTD算法

1.1 FDTD方程

時域有限差分法(FDTD)是求解三維時域電磁場問題的主要方法之一，其不需要求解大型方程，因此求解速度快且能夠較好地實行并行運算。基于此，現采用時域有限差分法來進行瞬變電磁場三維模擬計算。FDTD算法是一種直接基于時域電磁場方程的數值算法，算法的核心思想是將電磁場的時域微分方程的求解轉化為差分方程的迭代求解，近年來被廣泛用于瞬變電磁場的數值計算。經常用到的Yee網格其磁場分量位于面的中心，而電場分量位于棱的中間[11-12]。 在計算瞬變電磁場時，使用準靜態條件下的麥克斯韋方程[13]：

(1)

(2)

(3)

(4)

式中：H(r，t)為磁場強度，A/m；E(r，t)為電場強度，V/m；σ(r)為介質電導率，S/m；μ為介質磁導率，H/m；γ為位移介電常數；r為電磁波源點到場點的位置矢量；t為電磁波傳播時間。電磁場的差分方程可以通過對微分方程[式(1)～式(4)]的離散來獲得[14-15]。

1.2 源的加載

全空間條件下磁偶源中電流被突然切斷所產生的瞬變磁場[16-17]為

(5)

1.3 穩定條件

電磁場時域有限差分法采用顯式差分格式，其空間網格大小和時間步長需要滿足一定的條件才能保證迭代的穩定。電磁場時域有限差分迭代的數值穩定性由Courant條件確定[18-19]。

(6)

式(6)中：ν為電磁波傳播速度；Δx、Δy、Δz分別為x、y、z方向的空間最小網格間距。

1.4 邊界條件

由于計算機模擬有限空間，截斷邊界會對計算結果產生影響，需采用合適的邊界條件進行消除，本文研究中采用適合瞬變電磁場計算的卷積完全匹配層(CPML)邊界進行求解。

2 低阻陷落柱影響下礦井全空間瞬變電磁場的擴散及響應特征

2.1 低阻陷落柱影響下的全空間瞬變電磁場的擴散

陷落柱是煤田廣泛存在的一種地質構造，其導水性強，易于充水，已成為礦井水害的元兇之一。建立煤系地層中含水陷落柱的井下全空間地球物理模型，模擬陷落柱對瞬變電磁的響應特征及規律，把控瞬變電磁探測陷落柱的機理，可有效預防陷落柱突水事故的發生[20-21]。

圖1 陷落柱模型Fig.1 Collapse column model

圖1為煤系地層含水陷落柱的井下全空間地球物理模型。含水陷落柱發育于煤層底板以下，其頂部距煤層20 m，底部距離煤層70 m，長20 m、寬20 m，頂部中心坐標為(60，60，-20)，電阻率設為10 Ω·m，巷道位于煤層中，其截面大小為6 m×6 m。煤層厚10 m，其電阻率及圍巖的電阻率為100 Ω·m。發射回線大小為2 m×2 m，位于煤層中，回線的中心為坐標系原點，其法向沿x軸方向，滿足右手坐標系，發射電流為1 A。

圖2為柱狀低阻陷落柱影響下的不同時刻x方向感應磁場分量Hx在xy平面(z=-50 m)的等值線平面圖，單位A/m。如圖2所示，磁場的最大值始終保持在xy平面的中心位置(z軸上)。隨著感應磁場從發射回線附近不斷向周圍空間擴散，其擴散范圍到達陷落柱位置時，感應磁場的分布受到陷落柱的影響。

圖3為柱狀低阻陷落柱影響下的不同時刻y方向感應電場分量Ey在yz平面(x=60 m)的等值線剖面圖，單位V/m。如圖3所示，感應電場的正負最大值分別位于z=0平面的上下兩側，隨著感應電場向上下空間不斷地擴散，感應電場的正負最大值也向上下空間移動，當電場擴散到陷落柱體內，電場分布受到陷落柱的顯著影響，而電場等值線畸變的范圍恰好與陷落柱的位置和深度吻合(圖3中虛線框所示)。當t=500 μs時，感應電場正負最大值已經擴散到周圍較遠空間，此時電場等值線近似水平分布。

2.2 基于合矢量的瞬變電磁場的響應特征

某一方向的場強(電場或磁場)分量的擴散過程，一定程度上反映了介質及異常體的特性，但很難體現出異常體的整體形態，而實際上，場強是x、y、z三個方向的合矢量，用合矢量來描述電磁場的擴散過程，更能體現低阻異常體的電性、形態、位置等響應特征。下面以電場為例，進行電磁場擴散特性的描述。

圖4 感應電場分量Ey、Ez及其合矢量Ey+z的等值線分布圖Fig.4 Contour distribution of induced electric field components Ey，Ez and their resultant vector Ey+z

圖5 感應電場分量Ex、Ey及其合矢量Ex+y的等值線分布圖Fig.5 Contour distribution of induced electric field components Ex，Ey and their resultant vector Ex+y

圖6為不同深度(從上至下分別為z=0 m、z=-50 m和z=-100 m)的x和y方向的電場合矢量Ex+y的在xy平面的等值線平面圖，圖7為不同橫向位置(從左到右分別為x=40 m、x=60 m和x=80 m)的y和z方向電場合矢量Ey+z在yz剖面的等值線剖面圖。如圖6所示，在初始階段，感應電場主要分布在z=0平面附近，其分布只受發射回線周圍介質電阻率的影響。隨著時間的推移，感應電場的擴散范圍慢慢擴大，其分布受到低阻陷落柱的影響。在t=500 μs時，感應電場已經進入陷落柱體內，它的等值線分布受到陷落柱的顯著影響，等值線在陷落柱內發生畸變。在不同深度的畸變的范圍與陷落柱模型的位置十分吻合，如圖6(b)所示，特別是在z=-50 m的平面上，能清晰地看到等值線畸變的范圍恰好是陷落柱在z=-50 m平面所在的位置。在不同橫向位置的yz剖面的電場合矢量Ey+z的等值線分布的畸變區域的位置和高度也恰好與陷落柱模型的位置和深度十分吻合，如圖7所示。當陷落柱完全包含于磁場畸變的中心時，此時的陷落柱已經成為一個“二次源”，向空間輻射電磁場。此時，觀測到的信號主要反映了瞬變電磁場在低阻體影響下的衰減情況。

通過對電磁場的擴散過程進行數值模擬，表明瞬變電磁場的分布會受到低阻體的顯著影響，并且可以通過電磁場的變化來確定低阻體的形態和位置。因此，可以通過觀測瞬變電磁場的變化來探測低阻異常體。

圖6 不同深度的x和y方向的電場合矢量Ex+y的等值線平面圖Fig.6 Plan views of contours of the vector Ex+y in the electrical field in the x and y directions at different depths

3 應用實例

基于前文研究的理論基礎，將瞬變電磁法應用到煤礦含水陷落柱的實際探測中，通過應用實例，進一步驗證了瞬變電磁法探測低阻陷落柱的有效性。

3.1 實例1

山西官地礦23511工作面掘進巷道沿3#煤層布置，距2#煤底板平均5 m，平均煤厚3.0 m，煤層傾角平均4°。頂板為泥巖、細砂巖及砂質泥巖互層，底板為泥巖。工作面整體為褶皺構造，陷落柱較發育。據2#煤采掘資料，該巷道預計揭露正斷層4條，其中斷距大于1 m的2條，均不存在導水現象；工作面內存在1個陷落柱，積水情況不明。為保障巷道的安全掘進，需提前查明前方100 m范圍內的陷落柱富含水情況。

圖7 不同橫向位置的y和z方向的電場合矢量Ey+z的等值線剖面圖Fig.7 Section views of contours of the vector Ey+z in the electric field in the y and z directions at different positions

為了能夠準較直觀地反映巷道掘進頭前方低阻異常體的方位和距離，常使用扇形探測技術[18]來描述前方介質的地電特性，即在視電阻率等值線圖上將每個探測方向上的視電阻率值按照時深轉換結果繪制到該方向對應的位置上，形成二維的扇形圖。其中由感應電動勢求取視電阻率參考文獻[19]，并通過時深轉換來獲得對應深度[20]。

如圖8所示，根據探測任務，每隔15°進行一個方向探測，共13個方向，形成從左幫到右幫180°范圍的扇形掃描。掃描結果經過視電阻率計算，時深轉換和坐標換算，形成的二維扇形視電阻率等值線擬斷面圖，如圖9所示。據圖9可知，在探測位置前方30～90 m、左幫10～80 m范圍視電阻率值相對正前方及右幫區域明顯偏低，且為明顯的視電阻率等值線低阻“圈閉”異常，結合水文地質資料解釋此區域為陷落柱含水造成的低阻異常。經礦方鉆探驗證為陷落柱含水水區，出水點位于掘進前方33 m處頂板上方4.5 m，積水范圍自巷道至左側幫75 m，出水量約20 m3/h，采取防治水措施后安全掘過了此區域。

圖8 瞬變電磁超前探測掃描方向示意圖Fig.8 Schematic diagram of the scanning direction of transient electromagnetic advance detection

圖9 瞬變電磁超前探測視電阻率等值線擬斷面圖Fig.9 Pseudo-section diagram of the apparent resistivity contour of the transient electromagnetic advanced detection

3.2 實例2

山西晉城礦區王臺鋪礦地質條件較復雜，為保障巷道的安全掘進，需要準確查明前方采空區及陷落柱含水情況。探測方式除同實例1中水平方向掃描巷道迎頭左幫到右幫180°范圍外(每15°一個方向，共13個方向)，如圖10(a)所示，還在剖面方向探測頂板向上+90°、+60°、+30°、水平0°、底板向下-30°、-60°、-90°共7個方向，如圖10(b)所示。

圖10 探測方式布置示意圖Fig.10 Schematic diagram of detection mode layout

探測結果經過視電阻率計算，時深轉換和坐標轉換，并利用Voxler成像軟件形成三維超前探測成果圖。圖中坐標原點為探測位置，+x方向沿巷道右幫方向，+y方向沿巷道掘進方向，+z方向為巷頂底板方向，形成的三維數據體立體圖和三維等值面圖分別如圖11和圖12所示，圖11和圖12中從藍色到紅色(冷色調到暖色調)表示視電阻率不斷升高。可以看出圖中存在3處較大的異常范圍，其中在迎頭正前方存在1處高阻范圍，推測為不含水區；在迎頭兩幫各存在1處低阻范圍；推斷為含水陷落柱，而在巷道兩幫都存在低阻響應推測為瞬變電磁的全空間效應引起的。

后經過礦方打鉆驗證，正前方50 m之內的高阻異常區，為無水區；結合其他地質資料，向右側幫的低阻異常區打鉆，于+x方向35 m處出水，水量較大，揭露為含水陷落柱。探測結果說明，瞬變電磁法對探測含水陷落柱十分有效。

圖11 三維數據體立體圖圖Fig.11 Three-dimensional view of a 3D data body

圖12 三維數據體等值面圖Fig.12 Isosurface map of 3D data body

4 結論

通過全空間三維時域有限差分數值算法，模擬了在低阻陷落柱影響下的全空間瞬變電磁場的擴散過程，并對電磁場某一方向上的分量及其合矢量進行了詳細分析，最后將結果應用于實際探測。得出如下結論。

(1)FDTD法是模擬瞬變電磁場特征規律的有效方法，該方法在求解瞬變電磁場微分方程時將時空都離散化，把空間劃分為網格，時間上按步進的方式迭代，在計算過程中的每一時間步都對全場所有的網格進行一次計算，因而能夠展示出場隨時間的波動在空間的傳播過程，并給出清晰的物理圖像。這種模擬過程可以為瞬變電磁法探測低阻陷落柱提供理論支撐。

(2)低阻體會嚴重影響瞬變電磁場的擴散過程。某一方向的場強分量的擴散過程，一定程度上反映了介質及異常體的特性，但很難體現出異常體的整體形態，而用合矢量來描述電磁場的擴散過程，更能體現異常體的電性、形態、位置等特征。

(3)多方向探測方法為掘進頭前方的目標定位提供了一種解決方案，該方法不僅可以檢測到低阻體，還可以確定其方位和距離。數值模擬和實際應用都證明了瞬變電磁法探測含水陷落柱的有效性。

(4)對于瞬變電磁的分支——鉆孔瞬變電磁或者地井瞬變電磁，目前常規的分析方法是針對不同單一分量進行曲線分析，確定異常方位和形狀，該方法分析起來十分繁瑣，而用合矢量來分析異常響應特征可能會找到一種簡便的數據處理解釋方法，這需要進一步研究。