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

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

0 引言

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

1 三維正演方法

1.1 控制方程

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

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

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

(5)

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

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

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

(6)

和

(7)

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

(8)

1.2 有限差分離散

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

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

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

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

(13)

1.3 激勵(lì)源的施加

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

(14)

式中：Js為源電流密度。

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

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

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

(15)

(16)

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

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

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

(17)

為保證電磁場(chǎng)傳播，取得滿足要求的時(shí)間步長。

對(duì)于頂邊界條件(地面)，采用Wang & Hohmann向上延拓邊界條件[14]，側(cè)邊界和底邊界使用Dirichlet邊界條件，電場(chǎng)切向分量和磁場(chǎng)法向分量為“0”。本文使用的邊界條件要求模型剖分?jǐn)U展到足夠大的區(qū)域。一般在背景電阻率為100 Ω·m時(shí)，對(duì)301×301×100個(gè)網(wǎng)格的均勻模型進(jìn)行FDTD計(jì)算時(shí)發(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)格尺寸。時(shí)間步的選取需滿足上述穩(wěn)定性條件。

2.1 水平回線源

研究回線源為水平狀態(tài)，起伏地形下層狀介質(zhì)模型的瞬變電磁響應(yīng)信號(hào)特征。溝谷或山峰均位于回線區(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。溝谷的剖面相對(duì)高差分別為Δh=0m、Δh=50 m、Δh=70 m、Δh=100 m，測(cè)點(diǎn)位于溝谷底部中心。

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

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

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

H型地層初段畸變程度較大(大于40%)，在0.000 04 s～0.001 s之間地形影響特征復(fù)雜，相對(duì)畸變量呈震蕩特征，0.001 s以后呈先增大后減小趨勢(shì)，但 0.000 04 s之后相對(duì)畸變量均小于20%。對(duì)比整段模擬結(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。山峰的相對(duì)高程分別為Δh=0 m、Δh=50 m、Δh=70 m、Δh=100 m，測(cè)點(diǎn)位于山峰頂部中心。

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

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

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

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

圖8 不同深度溝谷影響電阻率斷面圖對(duì)比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ā)射框隨地形起伏，中間地臺(tái)高差H(50 m、100 m)。數(shù)據(jù)采集點(diǎn)位于地臺(tái)左側(cè)、地臺(tái)上、地臺(tái)右側(cè)。

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

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

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

3 地形影響實(shí)例

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

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

4 結(jié)論

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

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

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

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

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