尚曉榮 岳明鑫* 楊曉冬 吳小平 李 勇

(①中國(guó)科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院，安徽合肥 230026； ②中國(guó)科學(xué)技術(shù)大學(xué)安徽蒙城地球物理國(guó)家野外科學(xué)觀測(cè)研究站，安徽蒙城 233500； ③中國(guó)地質(zhì)科學(xué)院地球物理地球化學(xué)勘查研究所自然資源部地球物理電磁法探測(cè)技術(shù)重點(diǎn)實(shí)驗(yàn)室，河北廊坊 065000)

0 引言

可控源電磁法(Controlled-source electromagnetic method)采用人工場(chǎng)源激發(fā)交流電磁波信號(hào)(10-3～105Hz)，通過(guò)測(cè)量目標(biāo)區(qū)交變電磁場(chǎng)達(dá)到研究地下介質(zhì)電性結(jié)構(gòu)的目的。由于成本低、勘探深度大、抗干擾能力強(qiáng)等特點(diǎn)，可控源電磁法在油氣和礦產(chǎn)資源勘探[1-2]、地?zé)豳Y源開(kāi)發(fā)[3-4]及環(huán)境工程[5-6]等領(lǐng)域發(fā)揮了重要作用。可控源電磁法包括時(shí)間域電磁法和頻率域電磁法。時(shí)間域電磁法利用不接地回線源或接地線源向地下發(fā)射一次脈沖磁場(chǎng)，在一次場(chǎng)間歇期間，利用線圈或接地電極觀測(cè)二次渦流場(chǎng)[7-8]。頻率域電磁法常見(jiàn)的是可控源音頻大地電磁法，主要使用接地水平電偶源作為人工信號(hào)發(fā)射源[9]，本文研究的即是頻率域可控源電磁法。近年來(lái)，為適應(yīng)“深地”國(guó)家重大探測(cè)戰(zhàn)略中電磁精細(xì)結(jié)構(gòu)探測(cè)的需求，可控源電磁儀器的研發(fā)及數(shù)據(jù)采集技術(shù)獲得了迅猛發(fā)展，數(shù)據(jù)的反演解釋正逐步邁入三維實(shí)用化階段[10-13]。與實(shí)際需求相比，三維可控源電磁數(shù)據(jù)的反演技術(shù)尚未成熟，計(jì)算效率和準(zhǔn)確性均有待提高，特別是考慮起伏地形的三維反演尚不多見(jiàn)[14-15]。正演是反演的基礎(chǔ)，反演解釋的準(zhǔn)確性和計(jì)算效率很大程度上依賴(lài)于正演算法，因此研究可控源電磁數(shù)據(jù)的三維高精度正演方法十分必要。

三維可控源電磁法常用的數(shù)值模擬方法主要包括積分方程法[16-18]、有限差分法[19-20]、有限體積法[21-22]以及有限單元法[23-26]等。任政勇等[27]基于有限元方法提出一種新的非結(jié)構(gòu)化網(wǎng)格局部節(jié)點(diǎn)加密方法，實(shí)現(xiàn)了復(fù)雜模型完全非結(jié)構(gòu)化四面體全自動(dòng)剖分； Puzyrev等[28]基于節(jié)點(diǎn)有限元方法，采用復(fù)雜迭代求解器和基于稀疏近似逆的有效預(yù)調(diào)節(jié)器對(duì)大型稀疏線性方程組進(jìn)行求解，實(shí)現(xiàn)了各向異性介質(zhì)中三維可控源電磁場(chǎng)的正演模擬； 蔡紅柱等[29]提出一種基于非結(jié)構(gòu)化網(wǎng)格的海洋電磁有限單元正演算法，該算法適用于井中電磁、航空電磁、環(huán)境地球物理等非均勻且各向異性介質(zhì)的電磁感應(yīng)基礎(chǔ)研究； 湯文武等[30]比較了基于節(jié)點(diǎn)有限元與矢量有限元的三維可控源電磁正演，結(jié)果表明相同條件下矢量有限元法的計(jì)算精度更高，但速度較慢； He等[31]基于棱邊有限元方法實(shí)現(xiàn)了各向異性介質(zhì)中三維張量可控源電磁法的正演。

常規(guī)可控源電磁數(shù)據(jù)解釋多是基于平坦地形，這會(huì)導(dǎo)致反演異常體的位置與形態(tài)常發(fā)生畸變，甚至產(chǎn)生虛假異常。隨著近年來(lái)地質(zhì)勘查目標(biāo)和野外實(shí)地勘測(cè)場(chǎng)景越來(lái)越復(fù)雜，學(xué)者們開(kāi)始研究起伏地形下的三維電磁正、反演。Nam等[32]采用棱邊有限元方法研究了大地電磁法在三維地形下的響應(yīng)特征； Lin等[33]分析了可控源音頻大地電磁測(cè)深數(shù)據(jù)的三維地形效應(yīng)； 殷長(zhǎng)春等[34]采用自適應(yīng)非結(jié)構(gòu)有限元法進(jìn)行了帶地形的時(shí)間域海洋電磁法正演。積分方程法、有限差分法和有限體積法大多是基于結(jié)構(gòu)化網(wǎng)格實(shí)現(xiàn)的，其中有限單元法對(duì)網(wǎng)格剖分比較靈活，非結(jié)構(gòu)化網(wǎng)格更適于模擬真實(shí)地形起伏[35]?；诜墙Y(jié)構(gòu)化網(wǎng)格有限元法的二維和三維海洋電磁問(wèn)題已經(jīng)取得了一些研究成果[36-38]。相比之下，起伏地形下的陸地三維可控源電磁數(shù)值模擬研究相較滯后，目前關(guān)于起伏地形對(duì)各個(gè)場(chǎng)分量的影響特征的研究尚不多見(jiàn)。

因此，本文基于矢量有限元算法開(kāi)展了起伏地形下的三維可控源電磁法多分量數(shù)值模擬研究。采用非結(jié)構(gòu)化四面體網(wǎng)格對(duì)計(jì)算區(qū)域進(jìn)行剖分，以準(zhǔn)確模擬起伏地形，并基于棱邊基函數(shù)實(shí)現(xiàn)電磁場(chǎng)的插值和求解。首先，構(gòu)建層狀介質(zhì)模型，將其三維矢量有限元計(jì)算結(jié)果與一維解析解[39]進(jìn)行對(duì)比，驗(yàn)證算法的正確性；其次，構(gòu)建簡(jiǎn)單異常體模型，驗(yàn)證算法對(duì)探測(cè)低阻異常體的有效性； 然后，分別討論了測(cè)點(diǎn)處為山脊地形和山谷地形、發(fā)射源處于山脊地形、發(fā)射源與測(cè)點(diǎn)間為山脊地形等情況下頻率域可控源電磁法的響應(yīng)特征； 最后，研究了復(fù)雜三維地質(zhì)模型的電磁響應(yīng)特征，驗(yàn)證該算法的實(shí)用性。

1 數(shù)值模擬理論

從麥克斯韋方程組出發(fā)，假設(shè)時(shí)諧因子為e-iω t的平面電磁波入射到均勻介質(zhì)中，則計(jì)算區(qū)域電磁場(chǎng)滿(mǎn)足

?×E=iωμH

(1)

?×H=σE-iωεE+J

(2)

式中：E為電場(chǎng)強(qiáng)度；H為磁場(chǎng)強(qiáng)度；ω為角頻率；μ為介質(zhì)磁導(dǎo)率；σ為電導(dǎo)率；ε為相對(duì)介電常數(shù)；J為外加電流密度。

由式(1)和式(2)可得整個(gè)計(jì)算區(qū)域內(nèi)可控源電磁法滿(mǎn)足的電場(chǎng)波動(dòng)方程

(3)

利用Garlerkin方法，有

(4)

式中：k2=εμω2+iσμω；V表示積分空間。

對(duì)式(4)采用基于棱邊的非結(jié)構(gòu)四面體矢量有限元法進(jìn)行求解，電場(chǎng)矢量場(chǎng)可表示為六個(gè)矢量棱邊插值基函數(shù)的線性組合。非結(jié)構(gòu)四面體單元如圖1所示。

圖1 非結(jié)構(gòu)四面體單元e示意圖e1～e6為棱邊編號(hào)

圖1中每一個(gè)四面體單元的電場(chǎng)強(qiáng)度可表示為[39]

(5)

由式(4)和式(5)并應(yīng)用Garlerkin方法，可得電場(chǎng)強(qiáng)度關(guān)于單元e的有限元方程

AeEe-iωWeEe-ω2CeEe=iωSe

(6)

式中：Ae、We、Ce是單元e的剛度矩陣；Se為單元e的質(zhì)量矩陣。矩陣Ae、We、Ce的第(n，j)個(gè)元素和矩陣Se的第n個(gè)元素可分別表示為

(7)

(8)

(9)

(10)

式中Je為四面體單元e內(nèi)電偶極子的電流密度。為計(jì)算總場(chǎng)方法下的單元源項(xiàng)矩陣Se，將水平線源置于四面體網(wǎng)格的棱邊上，假設(shè)線源的中點(diǎn)坐標(biāo)為(x0,y0,z0)，則單元e內(nèi)的電流密度Je可表示為

Je=Iδ(x-x0)δ(y-y0)δ(z-z0)

(11)

式中：I表示穿過(guò)四面體單元的源電流強(qiáng)度； δ為Dirac-δ函數(shù)。

對(duì)起伏地形模型的正演方程組(式(3))求解時(shí)，右端項(xiàng)J不僅包含地下電阻率異常體引起的散射電流密度，還包括起伏地形產(chǎn)生的散射電流密度。使用非結(jié)構(gòu)有限元法剖分模型，同時(shí)對(duì)地形起伏區(qū)域進(jìn)行網(wǎng)格加密，以更好地對(duì)地形進(jìn)行擬合。

將單元矩陣進(jìn)行組裝，可得到大型稀疏復(fù)線性方程組

KE=S

(12)

式中：K=Ae-iωWe-ω2Ce，為正演系數(shù)矩陣；S=iωSe，反映源的作用。這里的E即為待求電場(chǎng)強(qiáng)度矢量。通常，K是一大型對(duì)稱(chēng)稀疏矩陣，采用CSR(Compressed Sparse Row，稀疏矩陣按行壓縮)格式進(jìn)行存儲(chǔ)以減小存儲(chǔ)空間。

最后，使用直接求解器Pardiso求解控制方程，可得到單元網(wǎng)格棱邊上的電場(chǎng)強(qiáng)度和磁感應(yīng)強(qiáng)度分量。

2 數(shù)值模擬算例

2.1 層狀模型

為驗(yàn)證算法的正確性，構(gòu)建一個(gè)三層地層的層狀模型。對(duì)模型的電、磁響應(yīng)特征進(jìn)行分析，并根據(jù)模擬結(jié)果與一維解析解的相對(duì)誤差(計(jì)算公式為(數(shù)值解-解析解)/解析解)分析本文算法的精度。

層狀模型的空間尺寸為100km×100km×100km，如圖2a所示。發(fā)射源位于y=-10km處，測(cè)點(diǎn)位于原點(diǎn)，發(fā)射頻率為0.1～5000Hz，共10個(gè)頻點(diǎn)：0.1、0.5、1、5、10、50、100、500、1000、5000Hz。用非結(jié)構(gòu)四面體網(wǎng)格剖分生成的最終網(wǎng)格分布見(jiàn)圖2b。采用Key[40]提出的一維可控源電磁正演算法與本文基于非結(jié)構(gòu)矢量有限元算法進(jìn)行模擬，并對(duì)比響應(yīng)曲線，結(jié)果見(jiàn)圖3。由圖可見(jiàn)，本文算法模擬結(jié)果與一維解析解吻合較好，電場(chǎng)的分量Ex實(shí)部和虛部的數(shù)值解與解析解的相對(duì)誤差的絕對(duì)值分別小于1%和2%，精度較高，驗(yàn)證了本文算法的正確性。

圖2 層狀模型幾何結(jié)構(gòu)(a)和網(wǎng)格剖分示意圖(b)

圖3 層狀模型Ex分量實(shí)部和虛部數(shù)值解與解析解對(duì)比(左)及其相對(duì)誤差(右)

在保證精度的前提下，分別測(cè)試了層狀地層模型不同剖分網(wǎng)格單元數(shù)情況下數(shù)值模擬所需要的時(shí)間，并統(tǒng)計(jì)了數(shù)值解與解析解的誤差，結(jié)果見(jiàn)表1。

從表1可以看出，網(wǎng)格數(shù)增加約1倍時(shí)，數(shù)值模擬耗時(shí)大約增加2倍； 網(wǎng)格數(shù)大于33萬(wàn)后，誤差下降較慢。因此，在綜合考慮耗時(shí)與精度的情況下，本文選擇33萬(wàn)網(wǎng)格數(shù)進(jìn)行后文模型剖分。

表1 層狀模型數(shù)值模擬時(shí)間和精度統(tǒng)計(jì)

2.2 簡(jiǎn)單組合異常體模型

為了分析本文算法對(duì)異常體電阻率的敏感性，構(gòu)建兩組簡(jiǎn)單組合異常體模型，背景為均勻半空間。模型包含4個(gè)相同的低阻(10Ω·m)/高阻(1000Ω·m)異常體，異常體尺寸為500m×500m×200m(圖4)。模型空間為100km×100km×100km； 發(fā)射源點(diǎn)位于y=-3km處； 測(cè)量點(diǎn)位于x軸的-1～1km、y軸-1～1km所圍區(qū)域的地面，x、y方向的點(diǎn)距均為200m； 發(fā)射頻率同前文層狀模型。對(duì)接收點(diǎn)附近進(jìn)行了網(wǎng)格加密，在遠(yuǎn)離觀測(cè)點(diǎn)處，網(wǎng)格剖分比較稀疏，模型最終網(wǎng)格剖分如圖5所示。

圖4 簡(jiǎn)單組合異常體模型

圖5 組合異常體模型網(wǎng)格剖分示意圖左：x=0剖面； 右：z=100m剖面。黃色方塊是異常體

分別對(duì)低阻異常體和高阻異常體兩個(gè)模型進(jìn)行模擬，結(jié)果見(jiàn)圖6。從圖中可以看出，低阻異常體產(chǎn)生的電場(chǎng)強(qiáng)度異常較高阻異常體明顯，且不同頻率時(shí)的相對(duì)異常(絕對(duì)值)也不同。

圖6 簡(jiǎn)單組合異常體模型的Ex振幅曲線(左)及其相對(duì)異常曲線(右)

為了進(jìn)一步分析背景模型、包含高阻異常體和低阻異常體模型的電場(chǎng)特征，繪制了0.5Hz的Ex振幅圖(圖7)?？梢钥闯觯惓ｓw的賦存區(qū)域在有、無(wú)異常體的情況下差異較大； 對(duì)比低阻和高阻異常體的情況可知，頻率域可控源電磁法對(duì)于低阻異常體更敏感。

圖7 0.5Hz時(shí)組合異常體模型地表Ex分量振幅圖(a)不存在異常體(背景模型)； (b)存在低阻異常體； (c)存在高阻異常體； (d)圖b與圖a的相對(duì)異常； (e)圖c與圖a的相對(duì)異常

2.3 三維起伏地形模型

2.3.1 測(cè)點(diǎn)處山脊地形模型

建立一個(gè)含有山脊地形的三維簡(jiǎn)單異常體模型(圖8)。山脊頂部長(zhǎng)、寬均為800m，底部長(zhǎng)、寬均為1600m，高度為400m，其中心點(diǎn)位于(0，0，200m)。假設(shè)大地背景電阻率為100Ω·m，異常體為不規(guī)則形狀(圖8中的藍(lán)色區(qū)域)，厚度為200m，長(zhǎng)、寬均為500m，電阻率為1Ω·m，中心點(diǎn)位于(0，0，750m)。采用非結(jié)構(gòu)化四面體單元對(duì)模型進(jìn)行網(wǎng)格剖分，對(duì)發(fā)射源和接收點(diǎn)處進(jìn)行網(wǎng)格加密，模型網(wǎng)格剖分結(jié)果見(jiàn)圖8。采用線源，源中點(diǎn)水平位置為(3km，-3km)，發(fā)射頻率同前文層狀模型采用的頻點(diǎn)； 接收點(diǎn)位于地面-1～1km(x)、-1～1km(y)區(qū)域內(nèi)，x、y方向測(cè)點(diǎn)間距均為200m。

圖9為圖8所示三個(gè)模型在測(cè)點(diǎn)(0，0，0)處的電場(chǎng)強(qiáng)度分量(Ex)和磁感應(yīng)強(qiáng)度分量(By)振幅計(jì)算結(jié)果。由圖可見(jiàn)，模型M02與模型M03的正演Ex和By振幅曲線較一致，而M01模型與二者相差較大，說(shuō)明山脊地形會(huì)削弱異常體的存在，同時(shí)也說(shuō)明正演時(shí)如果忽略山脊地形，會(huì)對(duì)正演結(jié)果產(chǎn)生較大影響。對(duì)比圖9b與圖9d可以發(fā)現(xiàn)，山脊地形對(duì)Ex振幅的影響大于對(duì)By振幅的影響。

圖8 山脊地形低阻異常體模型網(wǎng)格剖分示意圖(a)水平地形的低阻異常體模型(M01)； (b)含山脊地形的低阻異常體模型(M02)； (c)含山脊地形的無(wú)異常體模型(M03)

圖9 山脊地形模型Ex和By振幅正演曲線及相對(duì)異常(a)Ex振幅； (b)Ex振幅相對(duì)異常； (c)By振幅;(d)By振幅相對(duì)異常

圖10為圖8中三個(gè)模型在頻率為1Hz時(shí)沿x方向(y=0)的正演分量Ex和By的振幅及其相對(duì)異常。可以看出，山脊根部位置的正演Ex和By振幅均高于水平地形下的場(chǎng)分量幅值，山脊頂部位置的正演Ex和By振幅均低于水平地形下的場(chǎng)分量幅值，這是由于在地形高處(即山脊頂部)，供電電流是發(fā)散的，會(huì)呈現(xiàn)低電流密度，而在地形相對(duì)低處(即山脊根部)，供電電流是收斂的，呈現(xiàn)高電流密度。此現(xiàn)象說(shuō)明山脊地形會(huì)削弱異常體所產(chǎn)生的異常。

圖10 1Hz時(shí)山脊地形模型Ex和By振幅曲線(a)Ex振幅； (b)Ex振幅相對(duì)異常； (c)By振幅;(d)By振幅相對(duì)異常

2.3.2 測(cè)點(diǎn)處山谷地形模型

建立一個(gè)含有山谷地形的簡(jiǎn)單低阻異常體模型(圖11)。山谷頂部長(zhǎng)、寬均為1600m，底部長(zhǎng)、寬均為800m，高度為400m，其中心點(diǎn)位于(0，0，200m)。異常體參數(shù)與前文的山脊模型低阻體相同。采用與山脊地形模型相同的加密剖分網(wǎng)格，在相同位置布設(shè)發(fā)射源、接收點(diǎn)，采用同樣的發(fā)射頻率。模型網(wǎng)格剖分結(jié)果見(jiàn)圖11。

圖11 山谷地形模型網(wǎng)格剖分示意圖(a)低阻異常體模型(M04)； (b)無(wú)異常體模型(M05)

圖12為水平地形的低阻異常體模型(M01)和圖11中的山谷地形模型M04和M05的Ex和By正演振幅曲線及相對(duì)異常曲線。由圖可見(jiàn)，M01、M04和M05三個(gè)模型的電、磁場(chǎng)強(qiáng)度振幅曲線差異較大； 與山脊地形的By正演振幅曲線(圖9c和圖9d)相比，山谷地形低阻體模型的差異更大，可見(jiàn)頻率域可控源電磁法的電場(chǎng)強(qiáng)度、磁感應(yīng)強(qiáng)度分量受山谷地形的影響更大。

圖12 山谷地形模型Ex和By振幅正演曲線及相對(duì)異常(a)Ex振幅； (b)Ex分量振幅相對(duì)異常； (c)By振幅; (d)By振幅相對(duì)異常

圖13為山谷模型在頻率為1Hz時(shí)沿x方向(y=0)的正演Ex和By振幅及振幅相對(duì)異常?？梢钥闯觯谏焦葍蓚?cè)Ex和By振幅均低于水平地形下的場(chǎng)分量幅值，在山谷底部則均相反。這一特征與山脊地形的計(jì)算結(jié)果相似，這是由于在地形高處，地表下供電電流發(fā)散，呈現(xiàn)低電流密度； 在地形相對(duì)低處(即山谷底部)，供電電流收斂，呈現(xiàn)高電流密度。此現(xiàn)象說(shuō)明山谷地形會(huì)增強(qiáng)異常體產(chǎn)生的異常。

圖13 1Hz時(shí)山谷地形模型Ex和By振幅沿x方向變化曲線(a)Ex振幅； (b)Ex振幅相對(duì)異常； (c)By振幅; (d)By振幅相對(duì)異常

2.3.3 發(fā)射源處起伏地形模型

構(gòu)建一個(gè)發(fā)射源處為起伏地形的模型(圖14)。源的中點(diǎn)位于起伏地形區(qū)域的中心位置(0，-8000m，-395m)，山脊頂部長(zhǎng)、寬均為1000m，底部長(zhǎng)、寬均為1600m，高度為400m。其他參數(shù)設(shè)置均與測(cè)點(diǎn)處山脊地形模型參數(shù)相同。發(fā)射源處地形平坦模型(M06)和地形起伏模型(M07)的網(wǎng)格剖分結(jié)果見(jiàn)圖14。

圖14 發(fā)射源處地形平坦模型(M06，上)和地形起伏模型(M07，下)網(wǎng)格剖分示意圖

圖15為模型M06和M07在測(cè)點(diǎn)(0，0，0)處的電、磁場(chǎng)分量Ex和By振幅隨頻率變化曲線及相對(duì)異常。從圖15a可以看出，頻率高于100Hz時(shí)，模型M07的Ex振幅高于模型M06； 頻率低于100Hz時(shí)，情況相反。這是因?yàn)轭l率越低，探測(cè)深度越大，發(fā)射源處的地形為山脊地形，相對(duì)增加了介質(zhì)厚度； 而且，對(duì)于高阻異常體，在相同的頻率下，電場(chǎng)衰減較快，振幅衰減也更快。高頻信號(hào)的探測(cè)深度較小，發(fā)射源處向下傳導(dǎo)的電流在相同深度接收點(diǎn)處的電場(chǎng)信號(hào)強(qiáng)度會(huì)更大。另外，從圖15b可以看出，兩個(gè)模型的Ex相對(duì)振幅差較大，最大值達(dá)48%。由圖15b和圖15d可見(jiàn)兩個(gè)模型的By分量相對(duì)異常與Ex分量有相似特征。

圖15 發(fā)射源處地形起伏模型Ex和By振幅正演曲線及相對(duì)異常(a)Ex振幅； (b)Ex分量振幅相對(duì)異常； (c)By振幅; (d)By振幅相對(duì)異常

2.3.4 發(fā)射源與測(cè)點(diǎn)間地形起伏模型

建立一個(gè)源點(diǎn)與測(cè)點(diǎn)間地形起伏模型(圖16)，分析這類(lèi)地形對(duì)可控源電、磁場(chǎng)響應(yīng)的影響特征。起伏地形山脊頂部長(zhǎng)、寬均為800m，底部長(zhǎng)、寬均為1600m，高度為400m，中心點(diǎn)為(0，-2000m，0)。發(fā)射源中點(diǎn)位于地表(8km，-8km)處。其他參數(shù)設(shè)置同2.3.1節(jié)模型。圖16為發(fā)射源和測(cè)點(diǎn)間地形平坦模型(M06)和地形起伏模型(M08)網(wǎng)格剖分圖。

圖16 發(fā)射源與測(cè)點(diǎn)間地形平坦模型(M06，上)和地形起伏模(M08，下)網(wǎng)格剖分示意圖

圖17為模型M06和M08在測(cè)點(diǎn)(0，0，0)處的電、磁場(chǎng)分量Ex和By振幅曲線及振幅相對(duì)異常曲線。從圖17a和圖17b可見(jiàn)，兩個(gè)模型的Ex振幅較接近，差異較小。另外，發(fā)射源與起伏地形的距離、地形起伏程度、地形與異常體的距離等因素也會(huì)影響電、磁場(chǎng)。圖17c和圖17d是兩個(gè)模型By分量振幅曲線和相對(duì)異常曲線，可見(jiàn)By分量振幅相對(duì)異常曲線與Ex分量有相似規(guī)律。

圖17 發(fā)射源與測(cè)點(diǎn)間地形起伏模型Ex和By振幅正演曲線及相對(duì)異常(a)Ex振幅； (b)Ex分量振幅相對(duì)異常； (c)By振幅; (d)By振幅相對(duì)異常

從上述分析可知，由于地形的存在，尤其是發(fā)射源之下起伏地形和測(cè)點(diǎn)處的起伏地形，采用頻率域可控源電磁法對(duì)地下異常體進(jìn)行探測(cè)時(shí)，電、磁場(chǎng)分量均發(fā)生了畸變，出現(xiàn)了虛假異常。因此在實(shí)際數(shù)據(jù)的處理解釋中，需要考慮地表地形對(duì)可控源電磁法電、磁場(chǎng)的影響。

3 復(fù)雜金屬礦應(yīng)用實(shí)例

為了驗(yàn)證本文算法在實(shí)際應(yīng)用中的可靠性，以湘西某銅鉛鋅多金屬礦床的成礦模式為基礎(chǔ)建立圖18所示地質(zhì)模型，所包含的主要巖(礦)石的電阻率信息如表2所示。

圖18 湘西銅鉛鋅多金屬礦床成礦模式圖

表2 湘西銅鉛鋅多金屬礦床主要巖(礦)石電阻率統(tǒng)計(jì)

根據(jù)表中巖(礦)石電阻率范圍，取其平均值，給出圖19所示的金屬礦床模型二維電阻率填圖結(jié)果及其非結(jié)構(gòu)四面體模型網(wǎng)格剖分結(jié)果。對(duì)原模型四周區(qū)域進(jìn)行延拓，其中空—地邊界沿原模型邊界兩端分別向兩側(cè)延伸，而延拓區(qū)域的地層電阻率設(shè)為與原模型相鄰區(qū)域的電阻率一致。延拓模型空間范圍為40km×40km×40km，包含四個(gè)長(zhǎng)度為1km的線源，發(fā)射電流為1A，源的中點(diǎn)坐標(biāo)分別為(0，-3000m，105m)、(-3000m，0，115m)、(0，3000m，162m)和(3000m，0，115m)(圖20)。模型中異常體分布區(qū)域?yàn)?750m≤x≤750m，-750m≤y≤750m，0≤z≤400m，低阻目標(biāo)體為圖19中填圖顏色為藍(lán)紫色的小規(guī)模片狀塊體及青綠色的塊狀異常體。測(cè)點(diǎn)沿地形起伏均勻分布，地表范圍為-750m≤x≤750m，-750≤y≤750m，測(cè)點(diǎn)間距為100m。發(fā)射頻率為0.1～5000Hz，頻點(diǎn)分布同2.1節(jié)層狀模型。

圖19 湘西銅鉛鋅多金屬礦床模型二維電阻率填圖(上)及其非結(jié)構(gòu)四面體單元剖分局部示意圖(下)

圖20 湘西銅鉛鋅多金屬礦床延拓模型發(fā)射源布置示意圖

分別正演計(jì)算有、無(wú)異常體情況下，金屬礦床模型中心點(diǎn)在x=0、y=-200m處Ex振幅隨頻率的變化曲線，結(jié)果見(jiàn)圖21a，二者的相對(duì)異常見(jiàn)圖21b。由于模型深部電阻率均一，因此低頻時(shí)Ex振幅較穩(wěn)定，且相對(duì)異常基本保持為常數(shù)； 在高頻段，由于地質(zhì)異常體分布復(fù)雜，圖21b中的相對(duì)異常曲線出現(xiàn)了明顯的波動(dòng)。

圖21 金屬礦模型偏移距x=0、y=-200m處Ex振幅隨頻率變化曲線(a)及相對(duì)異常曲線(b)

選取正演計(jì)算的頻率1Hz時(shí)的Ex和By振幅數(shù)據(jù)，繪制圖22所示的等值線圖?？梢园l(fā)現(xiàn)，Ex和By振幅異常區(qū)域均與低阻地層的平面位置大致吻合(圖22中虛線橢圓區(qū)域)； 同時(shí)，Ex和By振幅相對(duì)異常圖中相對(duì)異常較大的區(qū)域(圖22中實(shí)線橢圓區(qū)域)也與實(shí)際低阻異常體的賦存區(qū)域存在很強(qiáng)的相關(guān)性。

圖22 金屬礦模型頻率1Hz時(shí)正演計(jì)算的地表位置Ex(上)和By(下)振幅等值線圖及相對(duì)異常圖(a)包含異常體； (b)不包含異常體； (c)圖b與圖a的相對(duì)異常

4 結(jié)論

本文基于非結(jié)構(gòu)四面體網(wǎng)格和矢量有限元法實(shí)現(xiàn)面向復(fù)雜地質(zhì)模型的可控源電磁法三維多分量數(shù)值模擬。通過(guò)對(duì)比層狀模型的數(shù)值解與一維解析解，以及組合異常體模型的電、磁場(chǎng)強(qiáng)度響應(yīng)特征分析，驗(yàn)證了本文算法的正確性和有效性，說(shuō)明了模型剖分方案的合理性。通過(guò)測(cè)點(diǎn)處山脊和山谷地形模型、發(fā)射源下山脊地形模型、發(fā)射源和測(cè)點(diǎn)間山脊地形模型的數(shù)值模擬，說(shuō)明了地形起伏對(duì)可控源電、磁場(chǎng)會(huì)產(chǎn)生顯著的影響。因此，在地形復(fù)雜的區(qū)域開(kāi)展可控源電磁工作必須考慮地形的影響。實(shí)際復(fù)雜帶地形金屬礦模型的計(jì)算結(jié)果也進(jìn)一步驗(yàn)證了該算法的實(shí)用性和可靠性。