郭 建 磊

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

0 引言

在沉積巖地區(qū)由于層理發(fā)育導致地下介質(zhì)電阻率隨電流方向發(fā)生變化，地層表現(xiàn)出很強的各向異性導電特性。地層電導率各向異性的研究最早可以追溯到1920年，Schlumberger在進行電法勘探時就發(fā)現(xiàn)了沉積地層的電性各向異性現(xiàn)象[1]。在近幾十年的歷程中，較多的地球物理學者對地層各向異性進行了研究[2-3]，目前學者們對電導率各向異性與矢量電磁場的研究取得了初步進展[4]。

對于地層電導率而言，各向同性與各向異性的區(qū)別在于地層電導率與方向的關系。地層各向同性指地層電導率不隨方向的改變而變化，地層各向異性指地層電導率隨方向的改變而變化。地層多發(fā)育為連續(xù)沉積地層[5]，當沉積層具有明顯層理時，導電大地不同方向的導電能力發(fā)生變化，地下巖層表現(xiàn)出宏觀電導率各向異性；另外地層在形成過程中有大量的孔隙和裂隙，孔隙和裂隙會被其他固體物質(zhì)或者流體物質(zhì)充填，最終導致地層具有強烈的各向異性特性[6-8]；地層開采會產(chǎn)生一定的擾動，改變原巖的構(gòu)造應力，引起上覆巖層的破壞和采空區(qū)裂隙分布的改變，產(chǎn)生不同形態(tài)、不同高度和密度的采空裂隙和裂縫[9]，經(jīng)過一定時間后裂隙和裂縫會充水形成裂隙含水體，由于構(gòu)造和含水充填的原因，裂隙含水體會產(chǎn)生電導率各向異性特性[10-11]。

瞬變電磁法目前是探測地層電性分布、地層富水性以及采空區(qū)積水的主要地球物理探測方法，在煤田、金屬礦等領域應用廣泛。隨著計算機的發(fā)展以及瞬變電磁理論的進步，眾多學者在瞬變電磁法各向異性三維正演、瞬變電磁三分量解釋及應用等方面做了大量工作，比如：周建美等基于有限體積法實現(xiàn)雙軸各向異性瞬變電磁三維正演，通過正演模擬發(fā)現(xiàn)水平方向電導率變化嚴重影響瞬變電磁磁場垂直分量響應，垂直方向電導率變化對磁場垂直分量響應幾乎沒有影響[12]；劉亞軍等基于有限體積法實現(xiàn)任意各向異性瞬變電磁三維正演，發(fā)現(xiàn)主軸水平方向電導率各向異性嚴重影響瞬變電磁垂直分量的磁場響應，且傾斜各向異性比水平各向異性的影響更為嚴重，通過水平各向異性信號能較清晰地判斷出各向異性主軸方向[13]；程久龍等基于時域有限差分法研究地層各向異性對礦井瞬變電磁的影響特征，發(fā)現(xiàn)相比HTI-X介質(zhì)，HTI-Y介質(zhì)和VTI介質(zhì)各向異性系數(shù)的變化對瞬變電磁場的影響顯著，并結(jié)合實例進行了說明[14]。

目前，瞬變電磁法無論從硬件到軟件均已經(jīng)實現(xiàn)了三分量數(shù)據(jù)采集與解釋[15-16]，并且模型試驗和實際資料處理結(jié)果表明三分量解釋方法能夠利用更多的瞬變電磁數(shù)據(jù)，獲得更高分辨率的處理結(jié)果[17]。瞬變電磁法的解釋方法主要有視電阻率解釋法[18]、浮動薄板解釋法[19]、煙圈理論解釋法[20]、一維反演解釋法[17]和擬地震成像法[21]，但上述解釋方法均是基于將地層假設為各向同性介質(zhì)為前提進行。前人研究結(jié)果表明介質(zhì)各向異性特性會嚴重影響瞬變電磁響應，尤其是三維異常體的各向異性影響會更為突出，因此，基于各向同性模型的傳統(tǒng)成像及反演成像技術(shù)不再適用于各向異性數(shù)據(jù)的處理解釋，這導致傳統(tǒng)的基于各向同性模型的成像方法和一維反演技術(shù)不再適用于各向異性數(shù)據(jù)的處理與反演解釋，并且隨著瞬變電磁技術(shù)的發(fā)展以及計算機計算性能的快速提升，三分量采集及三分量聯(lián)合反演解釋等成為目前實現(xiàn)瞬變電磁精細探測的前沿技術(shù)[15-17]。前人的研究工作更多的是探討電導率各向異性對瞬變電磁磁場垂直分量響應的影響，并未涉及對磁場水平分量響應的影響。研究地層各向異性和含水體各向異性對回線源瞬變電磁三分量磁場響應的影響特征，有利于地質(zhì)精細化勘探。該項研究最便捷的途徑就是三維正演計算。

時域有限差分算法經(jīng)過近些年的快速發(fā)展，其理論不斷完善并趨于成熟。孫懷鳳等[22]通過安培環(huán)路定理方程將矩形回線源電流密度加入麥克斯韋方程組，實現(xiàn)了回線源的加入，實現(xiàn)TEM全波場正演模擬；盧緒山[23]在孫懷鳳的基礎上利用Open-ACC Fortran語言進行了以異構(gòu)并行算法開發(fā)，加快了計算速度；李展輝等[24]將復頻率完全匹配層吸收邊界應用到三維瞬變電磁時域有限差分正演中，給出了CFS-PML在TEM正演中參數(shù)設置準則；孫懷鳳等[25]后續(xù)又在常規(guī)FDTD的迭代基礎上開發(fā)了Crank-Nicolson FDTD正演算法，在保證求解精度的同時計算效率大幅提高。

本文基于三維時域有限差分算法，引入軸向電導率張量構(gòu)建各向異性電磁場控制方程，進而推導無源和有源條件下軸向各向異性電場三分量迭代公式，其中磁場分量的迭代計算公式與各向同性磁場迭代計算公式一致，以電流密度的形式加入各向異性麥克斯韋方程組安培環(huán)路定理實現(xiàn)電流源的加載，實現(xiàn)各向異性三維瞬變電磁三分量正演。通過與各向同性及各向異性半空間模型解析結(jié)果對比驗證算法的準確性，分別建立各向異性半空間模型、層狀模型、含水體模型并進行正演計算，研究x、y、z軸各向異性分別對瞬變電磁?Bx/?t、?By/?t和?Bz/?t響應的影響。

1 軸向各向異性有限差分瞬變電磁正演算法

1.1 正演算法原理

各向異性、有耗、無源媒質(zhì)中的Maxwell方程組表示為[14]：

式中：E為電場強度，B為磁感應強度，H為磁場強度，ε為虛擬介電常數(shù)，t為時間，σ為張量電導率。

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

(2)

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

圖1 Yee晶胞網(wǎng)格剖分示意Fig.1 Schematic diagram of Yee unit cell meshing

圖2 電磁場時間采樣分布示意Fig.2 Schematic diagram of time sampling distribution of electromagnetic field

空間上采用后向差分，時間上采用中心差分，將式(2)在tn+1/2時刻進行空間和時間離散，整理后得到電場強度隨時間的推進公式：

(3)

為分析軸向各向異性瞬變電磁響應特征，將任意各向異性簡化為軸向各向異性，σ為各個節(jié)點上等效軸向各向異性電導率張量，通常表示為：

(4)

將σ分別代入矩陣P、Q，求得:

(5)

以差分代替微分，采用非均勻網(wǎng)格形式[22]進行模型剖分，將式(3)在各分量所在空間節(jié)點處進行空間和時間離散，得到電導率各向異性下電場強度迭代表達式：

式中：i、j、k分別表示x、y、z方向網(wǎng)格剖分的節(jié)點號，Δx、Δy、Δz分別表示剖分網(wǎng)格的長度，Δt表示時間間隔。

采用回線源中加入階躍電流的方式進行電流源加載，采用梯形波發(fā)射波形，考慮發(fā)射電流的上升沿、持續(xù)時間和下降沿，有源區(qū)域式(2)可修改為：

(9)

式中Js為源電流密度；因電流在水平方向，因此方程中僅存在Jsx和Jsy。在直角坐標系中加入電流密度展開后，將發(fā)射回線所在網(wǎng)格按照差分格式正常離散,得到有源區(qū)域網(wǎng)格電場迭代形式[22]：

磁場迭代公式不包含電導率項，各向異性磁場迭代公式與各向同性迭代公式相同，具體推導過程可以參看參考文獻[22,23]。

瞬變電磁法是寬頻帶勘探方法，加載梯形波作為激發(fā)源。選用第一類邊界條件，即在6個邊界面上將電場和磁場分量全部強制賦零，計算過程中就需要滿足時間域和空間域的穩(wěn)定性條件，即Courant穩(wěn)定性條件:

(12)

式中v表示電磁波的相速度，保證電磁場傳播的因果關系。

1.2 正演算法驗證

采用各向同性和各向異性半空間模型，通過與一維數(shù)值解的比較驗證本文算法的精度。采用笛卡爾正交坐標系，模型中發(fā)射回線采用190 m×190 m的方形回線，回線源放置于xoy平面且回線源中心位于坐標系中心，發(fā)射電流為1 A，電流方向呈順時針，上升沿和下降沿均為1 μs，脈沖持續(xù)時間10 ms，二次場采樣時間為10 ms，模型剖分均勻網(wǎng)格尺寸為10 m，模型剖分網(wǎng)格數(shù)為221×221×200。

首先，建立各向同性半空間模型。模型電導率為0.01 S/m，考慮到水平分量磁場響應在回線源中心點位置結(jié)果為0，故將接收點選擇在(-60 m,-60 m,0 m) 處，分別接收?Bx/?t、?By/?t、?Bz/?t響應，與一維解析結(jié)果[16]進行對比。對比結(jié)果(圖3)顯示：?Bx/?t、?By/?t、?Bz/?t響應的解析解與數(shù)值解的誤差均在10%以下，在0.02～10 ms時間范圍內(nèi)誤差均在5%以下，?Bz/?t的響應誤差稍大于?Bx/?t和?By/?t的響應誤差，接收點位于發(fā)射回線的對角線上，?Bx/?t與?By/?t響應曲線重合且響應誤差一致。

圖3 各向同性半空間模型垂直分量磁場響應(a)與相對誤差(b)Fig.3 Isotropic half-space model vertical component magnetic field response (a) and relative error (b)

然后，建立如圖4所示的各向異性半空間模型，模型軸向電導率參數(shù)如表1所示。在回線源中心點處接收?Bz/?t響應，圖5為各向異性半空間模型有限差分數(shù)值解與一維解析解[30]的對比結(jié)果。通過圖5可以發(fā)現(xiàn)：3個模型的?Bz/?t響應一維解析結(jié)果與數(shù)值解的誤差均在5%以下，滿足計算精度要求。

圖4 各向異性半空間模型Fig.4 Anisotropic half-space models

圖5 各向異性半空間模型垂直分量磁場響應(a)與相對誤差(b)Fig.5 Anisotropic half-space model vertical component magnetic field response (a) and relative error (b)

表1 圖4模型的軸向電導率Table 1 Axial conductivity of the Fig.4 model

通過各向同性和各向異性半空間模型對比，驗證了本文算法的準確性與精度。本文主要分析地層的軸向各向異性對磁場三分量的影響特征，因此為考慮成圖方便，下文中將接收點坐標點選擇在(-60 m,-60 m,0 m) 處。

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

2.1 各向異性半空間模型

基于上述各向異性半空間模型進行正演，研究回線源瞬變電磁三分量響應(圖6)。圖6中，由于在雙對數(shù)坐標系中無法表示負值，因此將感應電動勢為負值的部分取絕對值后繪制于雙對數(shù)坐標系中，后續(xù)繪圖采用如上約定。通過?Bx/?t響應曲線

圖6 各向異性半空間模型瞬變電磁磁場三分量響應Fig.6 Three-component TEM response of anisotropic half-space models

可以發(fā)現(xiàn)：增大x軸電導率后響應幅值整體增大；增大y軸電導率后曲線出現(xiàn)“變號”現(xiàn)象，“變號”時間約為0.1 ms，響應幅值增大且增大幅度大于增大x軸電導率；增大z軸電導率對響應幾乎沒有影響，與背景響應重合。通過?By/?t響應曲線可以發(fā)現(xiàn)：增大x軸電導率后曲線出現(xiàn)“變號”現(xiàn)象，“變號”時間為0.1 ms左右，響應幅值增大且增大幅度大于增大y軸電導率；增大y軸電導率后響應幅值整體增大；增大z軸電導率對響應幾乎沒有影響，與背景響應重合。通過?Bz/?t響應曲線可以發(fā)現(xiàn)：增大x、y軸電導率后響應幅值整體增大且響應曲線重合，此時接收點位于對角線上，x、y軸電導率各向異性對?Bz/?t響應的影響程度一致；增大z軸電導率對?Bz/?t響應幾乎沒有影響，與背景響應重合。

為研究軸向電導率各向異性對三分量磁場等值線整體分布的影響，基于上述模型，選取關斷后0.1 ms時刻的三分量磁場值繪制地面三分量磁場強度等值線(圖7)，并以各向同性半空間模型三分量磁場響應分布圖為參考背景。通過?Bx/?t等值線分布圖可知，增大x軸電導率后等值線圖左、右兩側(cè)幅值增強，等值線漸變?yōu)榻茩E圓形；增大y軸電導率后等值線圖左、右兩側(cè)幅值增強且增大幅度大于增大x軸電導率，等值線圖左、右兩側(cè)出現(xiàn)“變號”現(xiàn)象；增大z軸電導率后等值線圖幾乎沒有變化。通過?By/?t等值線圖可知，增大x軸電導率后等值線圖左、右兩側(cè)幅值增強且增大幅度大于增大y軸電導率，等值線圖左、右兩側(cè)出現(xiàn)“變號”現(xiàn)象；增大y軸電導率后等值線圖左、右兩側(cè)幅值增強，等值線漸變?yōu)闄E圓形；z軸電導率增大后等值線圖幾乎沒有變化。通過?Bz/?t等值線圖可知，增大x、y軸電導率后等值線圖幅值增強，等值線由同心圓向外擴散模式變?yōu)?個極值橢圓中心向外擴散模式，增大z軸電導率后等值線圖幾乎沒有變化。

圖7 地面磁場三分量0.1 ms時刻分布情況Fig.7 Ground three-component magnetic field distribution at 0.1 ms

圖6和圖7均證明水平方向電導率對磁場三分量的影響較大，垂直方向電導率對磁場三分量基本沒有影響，這主要是因為回線源產(chǎn)生的電流主要是水平方向的，因此三分量響應主要受水平方向電導率的影響，垂直方向電導率的影響很?。粂方向的電導率對?Bx/?t響應的影響大于x方向的，x方向的電導率對?By/?t響應的影響大于y方向的。

2.2 各向異性層狀地層

沉積巖具有層理結(jié)構(gòu)。當沉積層具有明顯層理時，導電大地不同方向的導電能力發(fā)生變化，地下巖層表現(xiàn)出宏觀電導率各向異性；另外地層在形成的過程中有大量的孔隙和裂隙，孔隙和裂隙會被其他固體物質(zhì)或者流體物質(zhì)充填，最終導致地層具有強烈的各向異性特性。因此，研究層狀地層各向異性對回線源瞬變電磁三分量磁場響應的影響特征具有現(xiàn)實意義。

以各向同性半空間模型為背景模型，設地層電導率為0.01 S/m，基于背景模型設置不同深度的地層為各向異性，分別改變上部地層(深度在0～100 m)、中間地層(深度在100～200 m)、下部地層(深度在100 m以下的地層)以及全部地層(地表以下地層)為各向異性地層，改變各向異性地層的3個方向電導率，電導率張量參數(shù)如表2所示。

表2 層狀模型的軸向電導率Table 2 Anisotropy parameters of axial conductivity

采用與前文1.2章節(jié)中正演算法相同的發(fā)射源和發(fā)射電流參數(shù)，經(jīng)過三維計算得到層狀各向異性地層瞬變電磁三分量響應結(jié)果(圖8)。圖8表明：增大x軸電導率，改變上部地層與改變?nèi)康貙拥?By/?t響應出現(xiàn)“變號”現(xiàn)象；增大y軸電導率，改變上部地層與改變?nèi)康貙拥?Bx/?t響應出現(xiàn)“變號”現(xiàn)象；增大水平方向電導率，改變上部地層與改變?nèi)康貙拥娜至宽憫€重合，改變中間地層和改變下部地層的三分量響應曲線重合，改變上部地層與改變?nèi)康貙拥娜至宽憫荡笥诟淖冎虚g地層和改變下部地層的三分量響應幅值；增大z軸電導率，所有模型的瞬變電磁三分量響應曲線全部重合。

通過圖8再次驗證了:水平方向電導率對磁場三分量的影響較大，垂直方向電導率對磁場三分量基本沒有影響;y方向電導率對?Bx/?t響應的影響大于x方向電導率，x方向電導率對?By/?t響應的影響大于y方向電導率;淺部地層的電導率各向異性對三分量磁場響應的影響占主導成分。

圖8 各向異性層狀模型瞬變電磁三分量響應Fig.8 Three-component TEM response of anisotropic layered models

2.3 各向異性含水體模型

設置各向異性含水體模型，研究含水體各向異性對瞬變電磁磁場三分量響應的影響。各向異性含水體模型如下：在電導率為0.01 S/m的各向同性半空間中放置一個頂部埋深50 m、規(guī)模為110 m×110 m×50 m的含水體。當含水體為各向同性時，軸向電導率均為1 S/m；分別增大含水體的x、y、z軸電導率,設置模型為各向異性，電導率張量如表3所示。

表3 含水體模型的軸向電導率Table 3 Anisotropy parameters of axial conductivity

采用與前文相同的發(fā)射源和發(fā)射電流參數(shù)，將含水體中心分別放于發(fā)射源正下方、x軸正向55 m、y軸正向55 m處(圖9)。圖10為含水體模型瞬變電磁磁場三分量響應結(jié)果，圖11為其x、y、z軸各向異性瞬變電磁水平分量響應對比。

圖9 含水體相對于回線源位置俯視圖Fig.9 Top view of the position of the hydrate relative to the loop source

由圖10a和圖11可知，當含水體位于源中心正下方時，增大x、y軸電導率后，?Bx/?t、?By/?t、?Bz/?t響應幅值大于增大z軸電導率和各向同性模型，且x軸各向異性模型與y軸各向異性模型的?Bx/?t、?By/?t、?Bz/?t響應曲線基本重合，z軸各向異性模型與各向同性模型的?Bx/?t、?By/?t、?Bz/?t響應曲線基本重合。此時，接收點x方向距含水體中心點距離與y方向相等，x軸各向異性模型的?By/?t響應幅值大于?Bx/?t響應幅值，y軸各向異性模型的?Bx/?t響應幅值大于?By/?t響應幅值，z軸各向異性模型的?Bx/?t響應曲線與?By/?t響應曲線重合;這再次說明了x軸各向異性對?By/?t分量響應的影響大于?Bx/?t，y軸各向異性對?Bx/?t分量響應的影響大于?By/?t。

將含水體沿x軸正向移動55 m并保持深度不變，此時接收點x方向距含水體中心點距離大于y方向。通過圖10b可知，增大y軸電導率后?Bx/?t、?By/?t、?Bz/?t響應幅值增大且幅度大于增大x軸電導率。這主要是因為異常體位于x軸，根據(jù)“煙圈理論”，通過異常體位置的磁力線主要與y軸電導率方向平行，受其影響較大。將含水體沿y軸正向移動55 m并保持深度不變，此時接收點x方向距含水體中心點距離小于y方向。通過圖10c可知，增大x軸電導率后?Bx/?t、?By/?t、?Bz/?t響應幅值增大且幅度大于增大y軸電導率，這主要是因為異常體位于y軸，根據(jù)“煙圈理論”，通過異常體位置的磁力線主要與x軸電導率方向平行，受其影響較大。通過圖11可知，模型位于回線源正下方時響應反映均大于模型中心位于偏x軸正向和偏y軸正向，根據(jù)“煙圈理論”，位于回線源正下方的異常體受回線源所激發(fā)的二次場能量最強。

圖10 各向異性含水體瞬變電磁三分量響應Fig.10 Transient electromagnetic three-component response of anisotropic aquifer

圖11 x、y、z軸各向異性瞬變電磁水平分量響應對比Fig.11 x,y,z axis anisotropic transient electromagnetic horizontal component response comparison diagram

3 結(jié)論

本文基于三維時域有限差分算法，引入軸向電導率張量構(gòu)建各向異性的電磁場控制方程，實現(xiàn)電導率軸向各向異性瞬變電磁三維正演。通過與各向同性及各向異性半空間模型一維解析結(jié)果對比,驗證算法的準確性。

建立各向異性半空間模型、層狀模型、含水體模型并進行正演計算，研究x、y、z軸各向異性分別對瞬變電磁?Bx/?t、?By/?t和?Bz/?t響應的影響,得到如下結(jié)論：

1)通過各向異性半空間模型發(fā)現(xiàn)：水平方向電導率對磁場三分量的影響較大，垂直方向電導率對磁場三分量基本沒有影響，三分量響應主要受水平方向電導率的影響，垂直方向電導率的影響很??；y方向電導率的?Bx/?t響應的影響大于x方向電導率，同理,x方向電導率的?By/?t響應的影響大于y方向電導率；

2)通過層狀各向異性模型發(fā)現(xiàn)：各向異性地層所處深度不同對瞬變電磁三分量磁場響應影響特征不同，認為淺部地層的電導率各向異性對三分量磁場響應的影響占主導成分；

3)通過各向異性含水體模型發(fā)現(xiàn)：接收點與各向異性含水體的相對位置不同會影響瞬變電磁三分量磁場響應，當接收點x方向距離含水體中心小于y方向時，x軸各向異性對三分量響應影響程度大于y軸各向異性，反過來，當接收點x方向距離三維地質(zhì)體中心距離大于y方向時，x軸各向異性對三分量響應影響程度小于y軸各向異性，并且模型位于回線源正下方時響應反映均大于模型中心位于偏x軸正向和偏y軸正向。

研究發(fā)現(xiàn)地層電導率軸向各向異性對瞬變電磁磁場三分量響應具有較大影響，為提高瞬變電磁精細勘探，需要進一步研究不同地電類型、不同探測方式的軸向各向異性瞬變電磁三分量磁場響應特征，并進一步研究軸向電導率對三分量磁場的物理影響機理，為三分量聯(lián)合探測及精細反演提供指導依據(jù)。

致謝：本研究在孫懷鳳等開發(fā)的TEM3DFDTD程序基礎上完成，在此表示感謝！