郭建磊，高小偉，侯彥威

(中煤科工集團西安研究院有限公司，陜西 西安 710077)

煤礦開采條件復雜，安全生產(chǎn)問題突出，其中，在掘進工作面發(fā)生重特大水災事故占比達51.16%[1]，因此，提高超前預報的準確性對保障煤礦安全生產(chǎn)至關重要。煤礦超前探測中一般采用鉆探和物探兩種方法，其中，鉆探方法結(jié)果可靠，但施工周期長、費用高，對巷道掘進生產(chǎn)影響較大。物探方法中，礦井瞬變電磁法[2]因無損、快速、信息豐富且對含水體反應靈敏等獨特優(yōu)點而被廣泛應用，但受巷道復雜環(huán)境影響較大。巷-孔瞬變電磁超前探測方法[3]有效結(jié)合了礦井瞬變電磁法和超前鉆孔的優(yōu)勢，該方法將發(fā)射線框布置在掘進工作面處，接收探頭放入鉆孔中，按間隔逐點采集三分量響應。孫懷鳳等[4]通過物理模擬試驗證明可以通過巷-孔瞬變電磁垂直分量信號在幅值和方向反號時間上存在的差異判斷掘進工作面前方是否存在異常構(gòu)造。陳丁等[5]采用積分方程法以煤層底板下方存在低阻體為模型研究了全空間條件下巷-孔瞬變電磁場的響應特征，結(jié)果表明，水平分量響應的延續(xù)時間優(yōu)于垂向分量，垂向分量響應的幅值強于水平分量，在實際觀測中不僅要觀測垂向分量，也要觀測水平分量。范濤[6-7]利用鉆孔瞬變電磁垂直分量實現(xiàn)二維擬地震反演，并進行了三維數(shù)值模擬、水槽物理模擬和井下現(xiàn)場模擬試驗；后續(xù)利用趨勢面提取技術(shù)提取了水平分量純異常場，通過異常場形態(tài)組合和幅值關系實現(xiàn)異常體中心方位角定位，并在此基礎上通過聚類算法對水平分量異常形態(tài)進行自動分類，完成異常體空間角度定位[8]。該方法達到“一孔多用”的目的，有效克服了電磁和人文干擾，具有觀測點距離異常體近、采集數(shù)據(jù)信噪比高等優(yōu)點，在全國多個煤礦獲得了廣泛應用[9-10]。

巷-孔瞬變電磁法在資料處理解釋過程中將地層假設為各向同性介質(zhì)，但隨著電磁勘探方法的發(fā)展，越來越多的學者認識到電導率各向異性在地下構(gòu)造中普遍分布[11-13]，并對介質(zhì)電導率各向異性特性與矢量電磁場關系做了大量研究。當前，關于頻率域電磁法，如大地電磁法[14-18]、海洋可控源電磁法[19-21]、直流電法[22-24]、井間電磁法[25-26]和航空電磁法[27-28]等方法的各向異性正反演研究均取得了較為豐富的成果。關于時間域電磁法的各向異性研究主要針對半空間介質(zhì)和電導率各向異性對垂直分量的影響。周建美等[29]基于有限體積法實現(xiàn)雙軸各向異性瞬變電磁三維正演，劉亞軍等[30]基于有限體積法實現(xiàn)任意各向異性瞬變電磁三維正演。在全空間介質(zhì)方面，程久龍等[31]基于時域有限差分法研究了不同主軸各向異性介質(zhì)對礦井瞬變電磁場垂直分量的影響特征，解釋了探測結(jié)果存在偏差的原因，并結(jié)合實例進行說明。上述研究并未涉及電導率各向異性對三分量響應的影響。

目前，巷-孔瞬變電磁法從硬件到軟件均實現(xiàn)了三分量信號的采集與解釋，通過將地層假設為各向同性介質(zhì)，依據(jù)水平分量形態(tài)組合及幅值關系判定孔旁異常體的分布象限，依據(jù)垂直分量反演計算異常體相關參數(shù)，最終實現(xiàn)巷-孔瞬變電磁立體成像解釋。由于資料解釋受制于對介質(zhì)電導率各向同性的假設，導致在存在明顯電導率各向異性的勘探區(qū)會產(chǎn)生較大的解釋誤差，因此，非常有必要研究全空間條件下電導率各向異性對巷-孔瞬變電磁響應的影響方式及程度，為進一步提高超前精度提供指導依據(jù)。本文研究采用時域有限差分算法(Finite Difference Time Domain，F(xiàn)DTD)[32]實現(xiàn)各向異性瞬變電磁三維正演，通過引入電導率各向異性張量構(gòu)建控制方程，以差分代替微分對控制方程進行離散，空間離散采用后向差分，時間離散采用中心差分，以電流密度的形式加入各向異性Maxwell 方程組的安培環(huán)路定理實現(xiàn)任意電流源的加載，進而推導出有源和無源區(qū)域電磁場迭代表達式。最后，構(gòu)建全空間模型、層狀模型和三維塊狀模型進行正演并分析電導率軸向各向異性對巷-孔瞬變電磁三分量響應的影響程度與方式。

1 正演理論基礎

1.1 控制方程

電導率各向異性、無源媒質(zhì)中的Maxwell 方程組為：

式中：E為電場強度，V/m；H為磁場強度，A/m；為張量電導率，S/m；μ為磁導率，Wb/(A·m)；ε為介電常數(shù)，F(xiàn)/m；t為時間，s。

地球物理低頻電磁勘探中一般忽略位移電流，為了滿足三維時域有限差分計算要求，因此，加入虛擬介電常數(shù)構(gòu)成顯式的時間迭代格式，式(1b)變?yōu)椋?/p>

式中：γ為虛擬介電常數(shù)。

采用經(jīng)典Yee 氏網(wǎng)格對計算區(qū)域進行剖分，電場分量位于棱邊中心點，磁場分量位于面中心點，每一個電場(磁場)分量均被相應的4 個磁場(電場)分量包圍，電場分量和磁場分量在時間軸上交替采樣。

采用非均勻網(wǎng)格形式[32]進行模型剖分，將式(2)在t(n+1/2)時刻進行空間和時間離散，以差分代替微分，其中，空間上采用后向差分、時間上采用中心差分，無源區(qū)域得到如下軸向各向異性的電場強度迭代公式：

在有源媒質(zhì)區(qū)域，式(2)必須包含源電流密度項，將式(2)修改為：

式中：Js為源電流密度。

巷-孔瞬變電磁法采用回線源進行發(fā)射，坐標系滿足“右手定則”，假定，巷道方向為z方向，上下為x方向，左右為y方向，發(fā)射源將被放置在xoy平面，因此，式(6)中僅存在x、y方向上的源電流密度項(Jsx、Jsy)?？紤]低頻近似后，在直角坐標系中將式(6)展開得到下式，按照差分格式正常離散，可以得到有源區(qū)域網(wǎng)格電場迭代公式。磁場迭代公式不包含電導率項，因此，各向異性條件下磁場迭代公式與各向同性條件下磁場迭代公式相同[32]。

1.2 算法驗證

構(gòu)建均勻全空間各向同性模型進行正演模擬，通過將數(shù)值解與解析解進行對比驗證算法精度。全空間介質(zhì)電導率為0.01 S/m，采用笛卡爾坐標系，回線源中心位于坐標系原點，邊長為3 m×3 m，電流為1 A，發(fā)射源上升沿和下降沿均為1 μs，持續(xù)時間為10 ms，二次場采樣時間為10 ms。最小網(wǎng)格尺寸為1 m，相鄰網(wǎng)格放大系數(shù)為1.05，剖分網(wǎng)格數(shù)為221×221×200，接收點位于(0.5 m，0.5 m，0 m)，在接收點處接收?Bx/?t、?By/?t、?Bz/?t響應。

數(shù)值解與解析解對比結(jié)果如圖1 所示。由圖1 可知，數(shù)值解和解析解的三分量響應曲線基本重合，當時間達到0.003 ms 時誤差控制在5%以內(nèi)，其中，?Bx/?t和?By/?t的相對誤差稍大于?Bz/?t的相對誤差，由于接收點位于坐標系對角線上，故?Bx/?t與?By/?t響應曲線重合且相對誤差一致。故本文算法滿足計算精度要求。

圖1 數(shù)值解與解析解對比結(jié)果Fig.1 Comparison of numerical solutions and analytical solutions

2 各向異性巷-孔瞬變電磁響應特征

2.1 各向異性全空間模型

巷-孔瞬變電磁法勘探環(huán)境為全空間介質(zhì)，因此，建立如圖2 所示模型研究全空間介質(zhì)電導率各向異性對巷-孔瞬變電磁三分量響應的影響方式與程度。如圖2 所示，采用上節(jié)的坐標系系統(tǒng)，巷道內(nèi)空氣電導率為0.000 1 S/m，空腔范圍為(-3～3 m,-3～3 m,0～-∞ m)。發(fā)射源平貼放置于巷道工作面且中心位于坐標原點，測點位于超前鉆孔中，鉆孔孔口位于(0.5 m,0.5 m,0 m)，鉆孔平行于z軸，長度范圍為0～80 m，按照1 m 的點間距從孔底至孔口逐點采集三分量響應數(shù)據(jù)。模型a為各向同性，電導率為σx=σy=σz=0.01 S/m；模型b 為x軸各向異性，電導率為σx=0.1 S/m、σy=σz=0.01 S/m；模型c 為y軸各向異性，電導率為σy=0.1 S/m、σx=σz=0.01 S/m；模型d 為z軸各向異性，電導率為σz=0.1 S/m、σx=σy=0.01 S/m。模型剖分與發(fā)射源參數(shù)與1.2 節(jié)一致。

圖2 全空間各向異性巷-孔瞬變電磁超前探測模型Fig.2 Full-space anisotropic tunnel-hole transient electromagnetic advance detection model

圖3 為4 個模型的三分量響應多測道圖，時間道的時間為0.06～0.21 ms。由圖3 可知，4 個模型的?Bx/?t、?By/?t多測道圖形態(tài)呈反“C”型、?Bz/?t多測道圖形態(tài)呈半“C”型；模型a 的?Bx/?t、?By/?t多測道圖形態(tài)和幅值一致，最大幅值位于25 m；?Bz/?t多測道圖最大幅值位于0 m，隨深度增加逐漸衰減；與模型a 相比，模型b 的?Bx/?t多測道圖形態(tài)和幅值在發(fā)射源處發(fā)生微弱變化，?By/?t多測道圖幅值整體增大且最大幅值位于10 m、至40 m 趨于0，?Bz/?t多測道圖幅值整體增大，至40 m 趨于0；與模型a 相比，模型c 的?Bx/?t多測道圖幅值整體增大且最大幅值位于10 m、至40 m 趨于0，?By/?t多測道圖形態(tài)和幅值在發(fā)射源處發(fā)生微弱變化；?Bz/?t多測道圖幅值整體增大、至40 m處趨于0。模型d 的三分量響應多測道圖特征與模型a 的基本一致。

圖3 三分量響應多測道圖Fig.3 Three-component responses multi-channel diagram

將4 個模型的0.06 ms 時刻的三分量響應曲線進行對比(圖4)發(fā)現(xiàn)，x軸各向異性主要影響?By/?t和?Bz/?t響應，y軸各向異性主要影響?Bx/?t和?Bz/?t響應，z軸各向異性對三分量響應基本不產(chǎn)生影響。根據(jù)“煙圈效應”，回線源產(chǎn)生的感應電流主要在水平方向，三分量響應受到水平方向電導率的影響較大，而受到垂直方向的電導率影響很??；x方向感應電流垂直穿過y軸電導率，y方向感應電流垂直穿過x軸電導率，因此，x軸電導率主要影響?Bx/?t響應，y軸電導率主要影響?By/?t響應。

圖4 三分量響應曲線對比(t=0.06 ms)Fig.4 Comparison of three-component responses (t=0.06 ms)

2.2 各向異性層狀模型

對于水平鉆孔而言，結(jié)合煤炭固有的地層特征，巷-孔瞬變電磁法勘探環(huán)境基本處于層狀介質(zhì)。當回線源位置固定后，地下不同深度地層受到一次場激發(fā)的程度不同，地層電導率各向異性對三分量響應的影響程度也不一樣。

以各向同性全空間模型為基礎，分別將距回線源不同距離的地層設置為各向異性(圖5)。如圖5 所示，模型a 為上部地層(40～∞ m)電導率分別為各向同性和各向異性、模型b 為中間地層(-40～40 m)電導率分別為各向同性和各向異性、模型c 為下部地層(-40～-∞ m)電導率分別為各向同性和各向異性。當?shù)貙訛楦飨蛲詴r，電導率σx=σy=σz=0.01 S/m；當x軸為各向異性時，電導率σx=0.1 S/m、σy=σz=0.01 S/m；當y軸為各向異性時，電導率σy=0.1 S/m、σx=σz=0.01 S/m；當z軸為各向異性時，電導率σz=0.1 S/m、σx=σy=0.01 S/m。模型剖分與發(fā)射源參數(shù)與1.2 節(jié)一致。

圖5 層狀模型Fig.5 Schematic diagram of the layered models

對比分析模型a-模型c 的不同電導率參數(shù)下0.06 ms 時刻的三分量響應曲線(圖6)。由圖6a 可知，對于回線源上部的地層，y軸各向異性對?Bx/?t響應影響明顯，且使?Bx/?t幅值由負變正，對?By/?t、?Bz/?t響應幾乎沒有影響；同時，x軸、z軸各向異性與各向同性的三分量響應曲線幾乎重合。由圖6c 可知，y軸各向異性對?Bx/?t響應影響明顯，與模型a 的?Bx/?t響應反映不同的是，模型c 的?Bx/?t響應幅值正負不變且幅值增大。因此，通過水平分量形態(tài)和幅值的相互關系可以判斷各向異性介質(zhì)所在方位及電導率主軸方向，這主要是因為模型a 的地層位于回線源上方，模型c的地層位于回線源下方，巷-孔瞬變電磁三分量響應具有方向性。由圖6b 可知，x軸、z軸各向異性和各向同性的?Bx/?t響應曲線重合且偏離y軸各向異性的?Bx/?t響應曲線；y軸、z軸各向異性和各向同性的?By/?t響應曲線重合且偏離x軸各向異性的?By/?t響應曲線；x軸和y軸各向異性的?Bz/?t響應曲線重合且與z軸各向異性和各向同性的?Bz/?t響應曲線偏離較大，說明距離發(fā)射源較近地層的電導率各向異性對三分量響應的影響較大，也再次說明巷-孔瞬變電磁三分量響應主要受水平方向電導率的影響。

圖6 三分量響應曲線對比(t=0.06 ms)Fig.6 Comparison of three-component responses (t=0.06 ms)

2.3 各向異性塊狀模型

煤炭資源廣泛分布于沉積巖中，沉積巖地區(qū)由于層理發(fā)育導致地下介質(zhì)電阻率隨電流方向發(fā)生變化[33]；煤炭開采產(chǎn)生不同形態(tài)、高度和密度的采空裂隙和裂縫，裂隙和裂縫充水后形成裂隙含水體[34]，上述情況均會導致地層介質(zhì)具有電導率各向異性特性。因此，本節(jié)主要研究各向同性地層中含有各向異性塊狀模型和各向異性地層中含有各向同性塊體模型。

1)各向同性地層中含有各向異性塊狀模型

建立如圖7 所示的模型，異常體1 的規(guī)模為20 m×20 m×20 m，異常體2 的規(guī)模為30 m×20 m×20 m，異常體3 的規(guī)模為20 m×30 m×20 m，3 個異常體的中心均位于(20 m，20 m，20 m)。將3 個異常體分別設置為各向同性和各向異性，為各向同性時電導率為σx=σy=σz=0.1 S/m；為x軸各向 異性時電導率 為σx=1 S/m、σy=σz=0.1 S/m；為y軸各向異性時電導率為σy=1 S/m、σx=σz=0.1 S/m；為z軸各向異性時電導率為σz=1 S/m、σx=σy=0.1 S/m。模型剖分與發(fā)射源參數(shù)與1.2 節(jié)一致。

圖7 塊狀異常體模型俯視圖Fig.7 Top view of block abnormal model

將3 個異常體的0.06 ms 時刻的三分量響應曲線進行對比(圖8)。由圖8 可知，當3 個異常體均為各向同性時，?Bx/?t響應大小為異常體3＞異常體2＞異常體1，說明y方向邊長改變影響?Bx/?t響應；?By/?t響應的大小為異常體2＞異常體3＞異常體1，說明改變x方向邊長改變影響?By/?t響應；?Bz/?t響應大小為異常體2=異常體3 ＞異常體1。當3 個異常體均為z軸各向異性時，其三分量響應形態(tài)、規(guī)律與各向同性時相同，但其水平分量響應幅度大于各向同性；當異常體均為x軸和y軸各向異性時，其?Bx/?t和?By/?t響應形態(tài)、規(guī)律與各向同性的相同，不同之處在于，當3 個異常體均為x軸各向異性時，異常體2 的?Bz/?t響應＞異常體3，當3 個異常體均為y軸各向異性時，異常體3 的?Bz/?t響應＞異常體2。

圖8 三分量響應曲線對比(t=0.06 ms)Fig.8 Comparison of three-component responses (t=0.06 ms)

將水平分量響應曲線進行進一步對比(圖9)。當3 個異常體均為x軸各向異性時，異常體1 的?By/?t響應＞?Bx/?t響應，異常體3 的?Bx/?t(?By/?t)響應＞異常體2 的?By/?t(?Bx/?t)響應。當3 個異常體均為y軸各向異性時，異常體1 的?Bx/?t響應＞?By/?t響應，異常體2 的?By/?t(?Bx/?t)響應＞異常體3 的?Bx/?t(?By/?t)響應。當3 個異常體均為z軸各向異性時，異常體1 的?Bx/?t響應和?By/?t響應曲線重合，異常體2 的?By/?t(?Bx/?t)響應與異常體3 的?Bx/?t(?By/?t)響應曲線重合，且異常體2 的?By/?t(異常體3 的?Bx/?t)響應＞異常體2 的?Bx/?t(異常體3 的?By/?t)響應。上述情況說明異常體邊長不同對水平分量響應的影響強弱不同，與模型長軸平行方向的電導率改變對三分量響應的影響較大。

圖9 水平分量響應對比(t=0.06 ms)Fig.9 Comparison diagram of horizontal component responses (t=0.06 ms)

2)各向異性地層中含有各向同性塊體模型

基于圖7 中的異常體1 模型，將異常體1 設置為各向同性，電導率為σx=σy=σz=0.1 S/m，將全空間介質(zhì)分別設置為各向同性和各向異性，為各向同性時電導率為σx=σy=σz=0.01 S/m；為x軸各向異性時電導率為σx=0.1 S/m、σy=σz=0.01 S/m；為y軸各向異性時電導率為σy=0.1 S/m、σx=σz=0.01 S/m；為z軸各向異性時電導率為σz=0.1 S/m、σx=σy=0.01 S/m。

將全空間介質(zhì)為各向同性和各向異性時0.06 ms時刻的三分量響應曲線進行對比(圖10)。由圖10 可知，全空間介質(zhì)為各向同性與z軸各向異性的三分量響應幾乎重合，?Bx/?t和?By/?t響應曲線在距離發(fā)射源40 m 后出現(xiàn)偏離；改變?nèi)臻g介質(zhì)x軸和y軸電導率對異常體的三分量響應均影響較大，使?Bx/?t和?By/?t異常幅值由正變負，且改變y軸電導率對?Bx/?t響應的影響＞改變x軸電導率、改變x軸電導率對?By/?t響應的影響＞改變y軸電導率，改變x軸和y軸電導率對?Bz/?t響應的影響一樣，這與上述小節(jié)的規(guī)律一致。全空間介質(zhì)為各向異性的三分量異常響應大于各向同性的三分量異常響應，說明異常體產(chǎn)生的異常響應被淹沒在全空間介質(zhì)電導率各向異性產(chǎn)生的異常響應中，如果此時不考慮介質(zhì)的各向異性特性會對巷-孔瞬變電磁解釋造成較大誤差。

圖10 三分量響應對比(t=0.06 ms)Fig.10 Comparison of three-component responses (t=0.06 ms)

3 結(jié)論

a.介質(zhì)距回線源的距離不同，其各向異性特性對巷-孔瞬變電磁三分量響應的影響強弱和方式不同，可以通過其三分量響應的形態(tài)和幅值的相互關系辨別各向異性介質(zhì)所在方位及主軸電導率方向，同時由于三維異常體的邊長不同對其三分量響應的影響強弱不同，進一步增加了各向異性條件下三分量解釋的復雜性。

b.回線源產(chǎn)生的感應電流主要是水平方向，并且x方向感應電流垂直穿過y軸電導率，y方向感應電流垂直穿過x軸電導率，因此，z軸各向異性對三分量響應基本沒有影響，x軸和y軸各向異性對巷-孔瞬變電磁三分量響應有較大影響，同時?Bx/?t響應主要受y軸電導率影響，?By/?t響應主要受x軸電導率影響。

c.當異常體所處介質(zhì)呈電導率各向異性時，異常體產(chǎn)生的異常響應被淹沒在全空間介質(zhì)電導率各向異性產(chǎn)生的異常響應中，介質(zhì)的電導率各向異性特征在巷-孔瞬變電磁法解釋過程中不可忽略，因此，本文研究內(nèi)容為進一步研究各向異性提供一定的基礎，也為后續(xù)進行各向異性反演提供參考。

d.未來的工作包括引入更復雜的電導率張量、進行亞網(wǎng)格剖分、添加吸收邊界條件等，實現(xiàn)任意各向異性的不規(guī)則復雜模型三維快速正演，為實際生產(chǎn)等提供重要的參考意義。