馬炳鎮(zhèn)， 郭建磊

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

0 引言

回線源瞬變電磁法在煤礦水文地質(zhì)勘查領(lǐng)域得到了廣泛應(yīng)用和發(fā)展[1-3]，但目前對于地形影響的響應(yīng)認(rèn)識仍需研究。隨著計算機(jī)和計算技術(shù)的發(fā)展，瞬變電磁數(shù)值模擬得到了快速的發(fā)展[4-6]，在三維時域有限差分正演方面， Wang T等[7]提出了一種三維電磁場的時間域有限差分方法，其采用均勻半空間瞬變電磁場解析解作為初始條件加入到迭代方程中；肖懷宇[8]采用Yee[9]的交錯網(wǎng)格方案和改進(jìn)的Du Fort-Frankel有限差分形式[10]來構(gòu)建三維正演，該算法適合于模擬任意電導(dǎo)率變化和地形起伏的地電模型;孫懷鳳等[11-12]將電流密度加入麥克斯韋方程組的安培環(huán)路定理方程，實現(xiàn)了源的加載，在計算區(qū)域內(nèi)通過電流源的激發(fā)和電磁場在有耗媒質(zhì)中的傳播特性形成激發(fā)一次場，這樣便為研究地表復(fù)雜，且回線源隨地形起伏的三維正演模擬提供了前提。這里是在以上三維時域有限差分正演研究基礎(chǔ)上，模擬并分析了回線源及隨地形起伏時的瞬變電磁響應(yīng)特征。

1 三維正演方法

1.1 控制方程

無源區(qū)域麥克斯韋方程組為

式中:E為電場強(qiáng)度;H為磁場強(qiáng)度;B為磁通量密度;σ為電導(dǎo)率;ε為介電常數(shù);t為時間。

瞬變電磁勘探中通常忽略位移電流，忽略位移電流后缺少電場對時間的導(dǎo)數(shù)，為了滿足三維運(yùn)算過程中對時間的迭代步要求，將式(2)介入虛擬介電常數(shù)得到式(5)。

(5)

式中:γ具有介電常數(shù)的量綱，稱其虛擬介電常數(shù)。

引入虛擬位移電流后的麥克斯韋方程組時域有限差分離散與方程(1)的離散方式類似。

在直角坐標(biāo)系中對Maxwell方程組中的式(1)和式(5)寫成分量的形式為

(6)

和

(7)

為了保證低頻條件下磁場計算結(jié)果的準(zhǔn)確性，必須將式(2)顯式的包含在迭代過程中。因此對于磁場的計算，可以通過電場先求解磁場的Hx和Hy兩個分量，然后通過這兩個分量以及Maxwell方程組中的式(4)來求解磁場的Hz分量。將方程(7)變形并由B=μH可得如下的磁場分量表達(dá)式方程(6)和方程(8)即無源區(qū)域電磁場計算的基本方程。

(8)

1.2 有限差分離散

直角坐標(biāo)系下的瞬變電磁場方程的有限差分離散使用Yee網(wǎng)格電磁場采樣格式和坐標(biāo)系統(tǒng)進(jìn)行網(wǎng)格離散(圖1)。

圖1 計算采用的Yee晶胞格式Fig.1 Calculation of the Yee unit cell format

以差商代替微分就可以對基本方程進(jìn)行求解，由于Euler前向差分對離散時間步的要求比較嚴(yán)格，故空間離散采用后向差分，時間離散采用中心差分，以f(x,y,z,t)表示電場或磁場在直角坐標(biāo)系中的某一分量，得到電磁場各分量一階偏導(dǎo)數(shù)的差分表達(dá)式為:

fn(i-1/2,j,k)]/Δx}+

O(Δx)

(9)

fn(i,j-1/2,k)]/Δy}+

O(Δy)

(10)

fn(i,j,k-1/2)]/Δz}+

O(Δz)

(11)

fn-1/2(i,j,k)/Δt]+O(Δt)}

(12)

磁場的z分量采用了式(13)的離散形式，n表示迭代步，其構(gòu)成了瞬變電磁場在有耗媒質(zhì)中傳播的電場和磁場的時域有限差分格式。

(13)

1.3 激勵源的施加

與大多學(xué)者采用初始條件來實現(xiàn)激勵源加載方法不同[7-8，13]，本文采用回線中加入階躍電流的方式進(jìn)行加載，源項為階躍源，且位于差分方程之中，計算過程中考慮激發(fā)源的上升沿、持續(xù)時間和下降沿，在有源區(qū)域，Maxwell方程中式(2)修改為式(14)。

(14)

式中：Js為源電流密度。

如圖2所示，將回線源施加在的網(wǎng)格棱邊，通過將矩形回線源電流密度加入麥克斯韋方程組的安培環(huán)路定理方程，實現(xiàn)回線源的加載。

圖2 回線源與網(wǎng)格位置示意圖Fig.2 Schematic diagram of loop source and grid location

對于地面瞬變電磁法，激勵源主要在水平方向，因此有源區(qū)域方程的直角坐標(biāo)分量形式為以中心差分格式離散方程，考慮源的加入方式將發(fā)射回線所在網(wǎng)格仍然按照差分格式進(jìn)行正常迭代，由于迭代格式中包含了源電流項，對源所在單元網(wǎng)格進(jìn)行特殊處理。根據(jù)法拉第電磁感應(yīng)定律和安培環(huán)路定理，根據(jù)法拉第電磁感應(yīng)定律和安培環(huán)路定理有：

(15)

(16)

源所在單元的電磁場值可以由上述積分求解得到，進(jìn)行源的處理時將其施加在Yee元胞的棱邊上，具體實現(xiàn)形式已有詳細(xì)的論文陳述[11]。

1.4 穩(wěn)定性條件及邊界條件

時域有限差分(FDTD)計算過程中需要滿足時間域和空間域的穩(wěn)定性條件：

(17)

為保證電磁場傳播，取得滿足要求的時間步長。

對于頂邊界條件(地面)，采用Wang & Hohmann向上延拓邊界條件[14]，側(cè)邊界和底邊界使用Dirichlet邊界條件，電場切向分量和磁場法向分量為“0”。本文使用的邊界條件要求模型剖分?jǐn)U展到足夠大的區(qū)域。一般在背景電阻率為100 Ω·m時，對301×301×100個網(wǎng)格的均勻模型進(jìn)行FDTD計算時發(fā)現(xiàn)可以迭代幾十萬步而不發(fā)散，能夠滿足精度需求。

2 地形影響特征正演模擬分析

模型采用非均勻剖分的形式，即在發(fā)射源和異常體處等關(guān)鍵區(qū)域采用均勻網(wǎng)格，按照10 m的立方體網(wǎng)格剖分，均勻網(wǎng)格往外以1.1倍網(wǎng)格擴(kuò)大，最大網(wǎng)格為30倍的均勻網(wǎng)格尺寸。時間步的選取需滿足上述穩(wěn)定性條件。

2.1 水平回線源

研究回線源為水平狀態(tài)，起伏地形下層狀介質(zhì)模型的瞬變電磁響應(yīng)信號特征。溝谷或山峰均位于回線區(qū)域內(nèi)。發(fā)射回線為600 m×600 m，發(fā)射電流為10 A。

2.1.1 溝谷地形

如圖3所示溝谷層狀模型(K型，H型)，K型模型參數(shù)為h1=100 m、ρ1=100 Ω·m，h2=200m、ρ2=1 Ω·m，ρ3=100 Ω·m；H型模型參數(shù)為h1=100 m、ρ1=10 Ω·m，h2=200 m、ρ2=1 000 Ω·m，ρ3=10 Ω·m。溝谷的剖面相對高差分別為Δh=0m、Δh=50 m、Δh=70 m、Δh=100 m，測點(diǎn)位于溝谷底部中心。

圖3 溝谷層狀介質(zhì)模型示意圖Fig.3 Schematic diagram of layered medium model in gully(a)溝谷模型;(b)溝谷與測點(diǎn)相對關(guān)系

圖4為不同高差下溝谷地形感應(yīng)電動勢二次場響應(yīng)曲線及對應(yīng)相對畸變量曲線(與水平地形相比)，從圖4中可看出,溝谷對早期信號影響嚴(yán)重，且溝谷高差越大，地形影響越嚴(yán)重，至晚期地形影響逐漸減弱。早期較水平地形表現(xiàn)為信號強(qiáng)度“減小”的特征，溝谷越深，信號強(qiáng)度越小；中、晚期信號較水平地形表現(xiàn)為信號強(qiáng)度“增強(qiáng)”的特征，溝谷越深信號強(qiáng)度越大，轉(zhuǎn)換為視電阻率后表現(xiàn)為“低阻異?！碧卣?。K型地層畸變曲線特征呈隨觀測時間增加，初段畸變量最大、后減小、再增大、最后減小的特征；在0.01 s之內(nèi)大部分區(qū)段相對畸變量均大于20%。

圖4 溝谷地形模擬結(jié)果Fig.4 Valley terrain simulation results(a)K型感應(yīng)電動勢曲線;(b)H型感應(yīng)電動勢曲線;(c)K型相對畸變量曲線;(d)H型相對畸變量曲線

H型地層初段畸變程度較大(大于40%)，在0.000 04 s～0.001 s之間地形影響特征復(fù)雜，相對畸變量呈震蕩特征，0.001 s以后呈先增大后減小趨勢，但 0.000 04 s之后相對畸變量均小于20%。對比整段模擬結(jié)果，H型地層溝谷地形影響畸變程度較小于K型地層。

2.1.2 山峰地形

如圖5所示山峰層狀模型，K型模型參數(shù)為h1=100 m、ρ1=100 Ω·m，h2=200 m、ρ2=1 Ω·m，ρ3=100 Ω·m；H型模型參數(shù)為h1=100 m、ρ1=10 Ω·m，h2=200 m、ρ2=1 000 Ω·m，ρ3=10 Ω·m。山峰的相對高程分別為Δh=0 m、Δh=50 m、Δh=70 m、Δh=100 m，測點(diǎn)位于山峰頂部中心。

圖5 山峰層狀介質(zhì)模型示意圖Fig.5 Schematic diagram of mountain layered medium model(a)山峰模型;(b)山峰與測點(diǎn)相對關(guān)系

圖6為不同高差下山峰地形感應(yīng)電動勢二次場響應(yīng)曲線及對應(yīng)相對畸變量曲線，從圖6中可以看出K及H型模型，早期較水平地形表現(xiàn)為信號“增強(qiáng)”的特征，山峰越高，信號強(qiáng)度越大；中、晚期信號較水平地形表現(xiàn)為“減弱”的特征，山峰越高，信號強(qiáng)度越弱，轉(zhuǎn)換為視電阻率后表現(xiàn)為“高阻異?！碧卣鳌Ｏ鄬兞壳€初段0.000 4 s之前呈先增大、后減小的特征(100 m山峰最大相對畸變量大于50%，50 m山峰最大相對畸變量大于20%)，0.000 4 s之后相對畸變量曲線呈先增大后減小的特征(100 m山峰H型最大相對畸變量約40%，K型最大相對畸變量約30%；50 m山峰H型最大相對畸變量約20%，K型最大相對畸變量約10%)，但H模型在0.000 1 s之后畸變特征更復(fù)雜，減小過程中有震蕩現(xiàn)象。從全時段模擬結(jié)果對比來看，山峰地形H型地層畸變程度亦較小于K型地層。

圖6 山峰地形模擬結(jié)果Fig.6 Mountain terrain simulation results(a)K型感應(yīng)電動勢曲線;(b)H型感應(yīng)電動勢曲線;(c)K型相對畸變量曲線;(d)H型相對畸變量曲線

2.2 起伏回線源

2.2.1 溝谷地形

如圖7所示模型介質(zhì)的電性ρ=100 Ω·m，發(fā)射源水平投影大小亦為600 m×600 m，圖8為溝谷示意圖，發(fā)射框隨地形起伏，中間溝谷凹槽高差H=50 m、100 m。數(shù)據(jù)采集點(diǎn)位于溝谷凹槽左側(cè)、凹槽內(nèi)、凹槽右側(cè)。

圖7 溝谷模型示意圖Fig.7 Schematic diagram of gully model

圖8為高差H=50 m、100 m時溝谷地形影響下深度視電阻率擬斷面圖(通過“煙圈效應(yīng)”等效電流環(huán)來計算深度[15])，從圖8中可以看出，發(fā)射線框隨地形起伏進(jìn)行布設(shè)時，地形影響使得整體信號“畸變”較為嚴(yán)重，100 m高差的溝谷最大相對畸變量達(dá)44.5%(不考慮較早期段)，表現(xiàn)為溝谷使得整體信號增強(qiáng)，呈“假低阻異?！保粶瞎仍缴?，地形影響約嚴(yán)重，晚期地形影響逐漸減弱，相對畸變量趨于10%以下(100 m高差)。

圖8 不同深度溝谷影響電阻率斷面圖對比Fig.8 Comparison of apparent resistivity sections affected by gullies of different depths(a)50 m;(b)100 m

2.2.2 山峰均勻介質(zhì)模型

如圖9所示模型介質(zhì)的電性ρ=100 Ω·m，水平發(fā)射框大小為600 m×600 m，發(fā)射框隨地形起伏，中間地臺高差H(50 m、100 m)。數(shù)據(jù)采集點(diǎn)位于地臺左側(cè)、地臺上、地臺右側(cè)。

圖9 山峰模型示意圖Fig.9 Schematic diagram of mountain model

圖10為高差H=50 m、100 m時山峰地形影響下深度視電阻率擬斷面圖(深度計算方式同圖8)，從圖10中亦可以看出，發(fā)射線框隨地形起伏布設(shè)時，地形影響使得整體信號發(fā)生“畸變”，100 m高差的山峰最大相對畸變量為25.6%(不考慮較早期段)，表現(xiàn)為山峰使得整體信號減弱，呈“假高阻異常”；山峰越高，地形影響約嚴(yán)重；晚期地形影響逐漸減弱，相對畸變量趨于10%(100 m正高差)。

圖10 不同高程山峰影響電阻率斷面圖Fig.10 Apparent resistivity section of different elevation peaks(a)50 m;(b)100 m

3 地形影響實例

勘探區(qū)位于山西河津某礦井，目標(biāo)層2號煤層最大埋深約600 m，采用邊長為600 m×600 m的瞬變電磁法大回線裝置，勘查主采煤層上覆含水層的富水異常范圍?？碧讲杉脑紨?shù)據(jù)生成的視電阻率擬斷面如圖11所示，紅色三角示意地表標(biāo)高，地形呈波狀起伏、溝谷連續(xù)分布，斷面兩側(cè)溝谷規(guī)模較大，深度約150 m，中間的溝谷較小，深度約100 m。溝谷位置的中淺部“低阻異?！迸c地形之間對應(yīng)關(guān)系密切，溝谷越深對應(yīng)的異常幅值越強(qiáng)，其異常形態(tài)與地表標(biāo)高起伏有一定的對應(yīng)性；隨著深度的增加，視電阻率等值線與地形的相關(guān)性明顯減弱，受地形影響程度減小。實測數(shù)據(jù)反映特征與正演結(jié)果基本一致，說明本次三維地形正演模擬結(jié)果正確。

圖11 實測數(shù)據(jù)電阻率擬斷面Fig.11 Apparent resistivity pseudo section

4 結(jié)論

采用時域有限差分模擬回線源瞬變電磁響應(yīng)特征，其結(jié)果與實測結(jié)果相一致?？偨Y(jié)如下：

1)通過三維時域有限差分正演模擬研究，在地形復(fù)雜勘探區(qū)，地形起伏超50 m時，大回線源瞬變電磁場應(yīng)考慮地形影響。

2)地形起伏是造成瞬變電磁響應(yīng)結(jié)果發(fā)生畸變的主導(dǎo)因素，不考慮較早期時間段，溝谷會導(dǎo)致“假低阻異常”，山峰會導(dǎo)致“假高阻異常”。

3)從模擬結(jié)果來看，K型模型地形影響畸變程度強(qiáng)于H型，但H型模型在0.000 1 s之后相對畸變量特征復(fù)雜，幅值有多次震蕩現(xiàn)象。

4)地形效應(yīng)早期較為嚴(yán)重，在本文的算例中，晚期響應(yīng)的影響減小。故在大回線TDEM的實測數(shù)據(jù)反演中，如采用一維反演技術(shù)，建議首先可對晚期數(shù)據(jù)做反演已獲取測區(qū)可靠的宏觀地電斷面，作為判斷全時間道反演結(jié)果是否合理的依據(jù)，并便于識別淺層可能由于地形畸變導(dǎo)致的假異常，避免錯誤的地質(zhì)解釋。