張繼鋒， 劉寄仁， 馮兵*， 董巖， 鄭一安

1 長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院， 西安 710054 2 中南大學(xué)地球科學(xué)與信息物理學(xué)院， 長(zhǎng)沙 410083

0 引言

傳統(tǒng)的地面電磁法(包括可控源音頻大地電磁法、瞬變電磁法等)是資源探查的重要手段，特別是在礦產(chǎn)資源及環(huán)境工程領(lǐng)域發(fā)揮了重要的作用(底青云等, 2019; Xue et al.，2020).隨著“十三五”國(guó)家“四深”發(fā)展戰(zhàn)略的提出，地球深部礦產(chǎn)資源的精細(xì)化探測(cè)提上日程，迫切需要發(fā)展新的技術(shù)方法.多道瞬變電磁法、廣域電磁法、短偏移距瞬變電磁以及地空電磁法等一批新的探測(cè)技術(shù)應(yīng)運(yùn)而生(陳衛(wèi)營(yíng)等，2016；薛國(guó)強(qiáng)等，2020；Ji et al., 2016).地空電磁法融合了地面電磁法大功率發(fā)射和無(wú)人機(jī)空中快速非接觸式采集的雙重優(yōu)點(diǎn)，能夠進(jìn)入地形復(fù)雜區(qū)域開(kāi)展大深度資源勘探，近年來(lái)成為地球物理電磁法研究的熱點(diǎn)(Allah et al., 2013, 2014; 李貅等, 2015；張瑩瑩和李貅，2017；Xue et al., 2018；Wu et al., 2019).該方法與無(wú)人機(jī)飛行平臺(tái)相結(jié)合，有利于在中小區(qū)域開(kāi)展大深度地下結(jié)構(gòu)精細(xì)探測(cè).相對(duì)于地面和航空電磁法，更加經(jīng)濟(jì)、安全和便捷，因此具有廣闊的市場(chǎng)前景和應(yīng)用價(jià)值.

地空電磁法的發(fā)展最早可追溯到20世紀(jì)70年代，加拿大TURAIR系統(tǒng)研制成功并用于安大略湖北部金屬硫化礦的探測(cè)(Pemberton et al., 1970; Bosschart and Seigel, 1972).該系統(tǒng)采用地面大回線作為發(fā)射源，在近空區(qū)域測(cè)量特定頻率下兩個(gè)水平線圈或正交線圈的磁場(chǎng)比值以及相位差，發(fā)射功率可達(dá)15 kW，探測(cè)深度可達(dá)到200 m.90年代后期，一系列新的地空儀器設(shè)備相繼出現(xiàn)，包括澳大利亞的FLAIRTEM系統(tǒng)，加拿大的TerraAir系統(tǒng)和日本GREATEM系統(tǒng)f621(Elliott,1998; Mogi et al., 1998 ).這幾種設(shè)備采集方式都是觀測(cè)一次場(chǎng)關(guān)斷后瞬變電磁響應(yīng)，解釋仍然采用地面回線源的解釋方法，只是前兩種系統(tǒng)仍然采用大回線發(fā)射源激勵(lì)，GREATEM系統(tǒng)采用地面電性源激勵(lì)方式.進(jìn)入21世紀(jì)，Smith等(2001)對(duì)地面、地空和航空瞬變電磁系統(tǒng)進(jìn)行了對(duì)比研究，系統(tǒng)分析了三種方法的噪聲水平和探測(cè)效果.結(jié)果表明：地空電磁系統(tǒng)的探測(cè)深度和數(shù)據(jù)質(zhì)量均優(yōu)于航空電磁系統(tǒng)，僅次于地面電磁系統(tǒng).此后，電性源地空系統(tǒng)GREATEM先后用于火山區(qū)和海岸帶地電結(jié)構(gòu)探測(cè)，深度可達(dá)800 m，成為一種重要的深部電磁勘探方法(Mogi et al., 2009; Ito et al., 2011, 2014).德國(guó)BGR研究所(Nittinger et al.，2017)也推出了一套頻率域地空電磁探測(cè)系統(tǒng)，并成功應(yīng)用于金屬礦探測(cè)中.目前地空電磁系統(tǒng)的研究主要以單輻射源為主，大量的研究集中在儀器研制和信號(hào)的去噪方面(嵇艷鞠等，2013；Wang et al., 2013；Qu et al., 2017; Li et al., 2017; 劉富波等，2017；Yuan et al., 2019).由于地空電磁法接收信號(hào)干擾較大，單源發(fā)射功率有限，探測(cè)深度和測(cè)量范圍都會(huì)受到限制，因此，迫切需要發(fā)展地空多源電磁勘探系統(tǒng)和方法.

多源勘探概念最早應(yīng)用在直流電阻率法中，Bibby(1986)提出多源測(cè)量模式并引入電阻率張量，該方法克服了視電阻率對(duì)目標(biāo)體位置和方向的依賴，減少了單源測(cè)量結(jié)果的“假異?！爆F(xiàn)象，證明多源測(cè)量能夠獲得更可靠的信息.Caldwell和Bibby(1998)在長(zhǎng)偏移距瞬變電磁(LOTEM)勘探中提出了多源瞬時(shí)電阻率張量，保留了直流電阻率張量的許多有用性質(zhì)，特別是在地下三維電阻率分布情況下，許多解釋技術(shù)可以應(yīng)用于時(shí)間域電磁數(shù)據(jù).90年代初，阻抗張量和傾子矢量被引入到地面可控源電磁測(cè)量中(Li and Pedersen,1991)，該方法一般采用兩個(gè)相互垂直的激勵(lì)源，分別從不同的方向激勵(lì)，兩個(gè)源單獨(dú)工作，能夠獲取更加豐富的地電信息，解釋結(jié)果更加可信，后來(lái)進(jìn)一步得到發(fā)展(Caldwell et al., 2002).席振銖等(2011)研究了正交磁偶源電磁場(chǎng)分布規(guī)律，實(shí)現(xiàn)了可控源高頻大地電磁張量測(cè)量；王顯祥等(2014)等研究了“L”型源的各分量表達(dá)式，設(shè)計(jì)了一種新的信號(hào)發(fā)射模式，實(shí)現(xiàn)了全張角范圍內(nèi)測(cè)量，顯示了張量可控源方法的優(yōu)勢(shì).

近年來(lái)，多源地空電磁系統(tǒng)的研究也逐漸受到了重視，李貅等(2015)提出了多輻射場(chǎng)源的全域視電阻率定義方式，不僅提供了一種定性直觀的解釋方法，也為三維多源地空瞬變電磁的發(fā)展提供了理論依據(jù).Zhou等(2018a，b)年提出頻率域多源電磁探測(cè)方法，主要從儀器的研制、激勵(lì)信號(hào)的同步以及多源之間的相互影響等方面進(jìn)行了研究，表明多源地空電磁勘探具有廣闊的發(fā)展前景.因此，需要進(jìn)一步從三維角度對(duì)多源電磁響應(yīng)進(jìn)行模擬和分析.目前地球物理電磁法三維數(shù)值模擬主要集中在單激勵(lì)源或被動(dòng)源三維模擬方面(李建慧等, 2016; Ren et al., 2013, 2014; 殷長(zhǎng)春等, 2017)，多源頻率域地空三維電磁響應(yīng)特征研究比較少，特別是全域視電阻率分布特征和規(guī)律.

本文采用多源三維有限元數(shù)值模擬方法，從多源輻射角度出發(fā)，設(shè)置了雙源、三源以及四源輻射方式，建立了低阻體、高阻體模型以及兩個(gè)低阻體模型，分別從總場(chǎng)、二次場(chǎng)以及視電阻率參數(shù)等方面對(duì)不同發(fā)射源模式進(jìn)行分析，總結(jié)每種源的電磁響應(yīng)和全域視電阻率特征，為實(shí)際選擇多源輻射方式進(jìn)行地空電磁勘探提供理論基礎(chǔ).

1 多源頻率域地空系統(tǒng)

地空頻率域電磁探測(cè)方法是采用地面布設(shè)發(fā)射源，空中測(cè)量垂直磁場(chǎng)響應(yīng)的工作模式.該方法由吉林大學(xué)教授林君提出，其儀器(GAFEM)已經(jīng)研制成功，采用偽隨機(jī)發(fā)射波形，在一些老礦區(qū)進(jìn)行了野外試驗(yàn)，取得了良好效果.單輻射源能量有限，信號(hào)強(qiáng)度弱，加之受運(yùn)動(dòng)噪聲的影響，信噪比和探測(cè)深度都會(huì)受到一定的影響.多場(chǎng)源激勵(lì)可改變空間磁場(chǎng)分布，提高信號(hào)強(qiáng)度，擴(kuò)大觀測(cè)區(qū)域，提高信噪比.通過(guò)多場(chǎng)源不同角度和不同方向激勵(lì)，可加強(qiáng)發(fā)射源與目標(biāo)體的耦合程度，改善對(duì)目標(biāo)體的分辨率(圖1).頻率域地空電磁法測(cè)量空中磁場(chǎng)垂直分量，不再測(cè)量電場(chǎng)分量，場(chǎng)的特征和解釋方法完全不同，地面的處理解釋方法不再適用.因此，必須基于空中磁場(chǎng)分量，研究三維地電模型的電磁場(chǎng)空間分布特征和新的解釋方法.

圖1 頻率域多源地空電磁采集方式示意圖Fig.1 Schematic diagram showing data acquisition for multiple-source frequency domain grounded-airborne electromagnetics

2 多源頻率域地空有限元分析

2.1 基于總場(chǎng)的公式

在準(zhǔn)靜態(tài)極限條件下的電場(chǎng)E和磁場(chǎng)H所滿足的頻率域方程為(湯井田和何繼善，2005):

(1)

(2)

其中假設(shè)電磁波的諧變因子為eiω t，式中μ為磁導(dǎo)率，本文μ選取真空磁導(dǎo)率μ0，σ為電導(dǎo)率，ω為圓頻率，Js為傳導(dǎo)電流密度.

由此便可導(dǎo)出電場(chǎng)雙旋度方程：

(3)

(4)

其中，V為研究區(qū)域的體積，在進(jìn)行節(jié)點(diǎn)有限元求解時(shí)，該方程不滿足散度條件，產(chǎn)生的偽解直接影響計(jì)算精度.為了提高計(jì)算精度，可在方程中增加一個(gè)罰項(xiàng)來(lái)削弱偽解產(chǎn)生的影響(金建銘，1998)：

(5)

2.2 有限元分析

本文采用非結(jié)構(gòu)化網(wǎng)格剖分方法，利用TetGen生成高質(zhì)量的四面體網(wǎng)格，將源和異常體附近場(chǎng)值變化劇烈的區(qū)域進(jìn)行網(wǎng)格加密，網(wǎng)格邊界大于10倍的趨膚深度，以保證場(chǎng)值在邊界上已衰減為零，如圖2所示.

圖2 非結(jié)構(gòu)化網(wǎng)格剖分Fig.2 Unstructured grid

圖3 四面體單元Fig.3 Tetrahedron element

(6)

若將整個(gè)研究區(qū)域剖分為n個(gè)單元，離散化后的變分方程為：

(7)

那么地空多源正演問(wèn)題轉(zhuǎn)化為求解該變分方程的極值問(wèn)題，該方程可看成多元函數(shù)，多元函數(shù)取極值應(yīng)當(dāng)滿足下列方程(徐世浙，1994):

(8)

其中i為節(jié)點(diǎn)編號(hào)，e為單元編號(hào).

由此可得到所有單元所滿足的矩陣方程：

(9)

按照單元節(jié)點(diǎn)的編號(hào)，可將每個(gè)單元生成的子矩陣組合成總體剛度矩陣，得到最終的離散矩陣方程：

(C+iωM)·E=-iωb.

(10)

由于產(chǎn)生的剛度矩陣是大型的對(duì)稱稀疏矩陣，直接存儲(chǔ)將會(huì)浪費(fèi)大量存儲(chǔ)空間，因此本文采用雙向鏈表組裝矩陣，按CSR方式存儲(chǔ).由于多頻可控源正演模擬的頻點(diǎn)之間是相互獨(dú)立的，解方程將會(huì)浪費(fèi)大量的時(shí)間，嚴(yán)重影響運(yùn)算效率，而經(jīng)前人研究發(fā)現(xiàn)，可以通過(guò)模型降階方法，將一個(gè)復(fù)雜的系統(tǒng)簡(jiǎn)化，只要經(jīng)過(guò)一次模型降階，便可實(shí)現(xiàn)多頻點(diǎn)的計(jì)算，提高運(yùn)算速率(蔣耀林，2010；B?rner et al., 2008；周建美等，2018).

2.3 模型降階方法

對(duì)于矩陣方程(10)，我們可以做如下變形：

(M-1C+iωI)·E=-iωM-1b，

(11)

令A(yù)=M-1C，X=M-1b，gω(A)=-iω(A+iωI)-1，解得E(ω)=-iω(A+iωI)-1X=gω(A)X.那么求解E(ω)的關(guān)鍵便是求解gω(A)，因此我們采用Krylov子空間(張繼鋒等，2009, 2020)投影算法，可降低矩陣運(yùn)算的維數(shù)，大大提高運(yùn)算效率.

模型降階算法實(shí)際上就是尋找一個(gè)有效逼近rm(A)，使得E≈rm(A)成立，求解rm(A)等價(jià)于求解矩陣A在m維Krylov子空間Qm+1=qm(A)-1κm+1(A,X)上的投影(周建美等，2018).由于有限元離散矩陣的特點(diǎn)，矩陣M并非單位矩陣，為了避免因傳統(tǒng)正交算法造成運(yùn)算的復(fù)雜性，我們引入一種M正交運(yùn)算(蔣耀林，2010; B?rner, et al., 2015)，進(jìn)而得到列正交向量Vm+1，下面給出模型降階后的最終表達(dá)式：

E=‖X‖MVm+1gω(Am+1)e1，

(12)

2.4 垂直磁場(chǎng)的求取

求解方程(10)可得到整個(gè)區(qū)域所有節(jié)點(diǎn)的電場(chǎng)三分量值，而在任意單元e中，任意點(diǎn)(x0,y0,z0)處的電場(chǎng)大小可通過(guò)插值基函數(shù)求得：

(13)

在地空系統(tǒng)中，電偶源在地面發(fā)射電流，線圈在空中接收磁場(chǎng)，因此本文通過(guò)研究z方向的磁場(chǎng)值Hz來(lái)研究大地的電性特征.由式(1)可得：

(14)

其中相應(yīng)場(chǎng)在某一點(diǎn)的導(dǎo)數(shù)值可通過(guò)式(13)代入式(14)求得.在實(shí)際數(shù)據(jù)處理和資料解釋中，為了能夠直觀的反映地下介質(zhì)的電性特征，我們需要將磁場(chǎng)轉(zhuǎn)化為視電阻率.張瑩瑩等(2015)定義了時(shí)間域瞬變電磁場(chǎng)的全域視電阻率，我們將其引入到頻率域，利用泰勒級(jí)數(shù)展開(kāi)，求得視電阻率的近似表達(dá)式，然后通過(guò)迭代算法得到全域視電阻率.

3 多源三維電磁響應(yīng)分析

3.1 有限元解和解析解的比較

圖4a、b分別給出了有限元解和解析解的垂直磁場(chǎng)分量以及相對(duì)誤差.可以看出，三維有限元計(jì)算結(jié)果同一維正演結(jié)果吻合的很好，磁場(chǎng)最大誤差為1.8%左右.

為了進(jìn)一步驗(yàn)證程序的正確性，我們?cè)O(shè)計(jì)了一個(gè)兩層模型，第一層電阻率為200 Ωm，厚度為100 m，第二層電阻率為50 Ωm.第一個(gè)源沿x軸正方向放置，電偶源的中心位置為(0 m，-4000 m，0 m)，第二個(gè)源沿x軸負(fù)方向放置，電偶源的中心坐標(biāo)為(0 m，2000 m，0 m)，設(shè)置測(cè)點(diǎn)位置為(1600 m，-2000 m，-100 m).如圖5所示，磁場(chǎng)相對(duì)誤差最大不超過(guò)3%，說(shuō)明基于模型降階的有限元計(jì)算結(jié)果滿足精度要求.

3.2 三維電磁響應(yīng)分析

3.2.1 低阻體模型

(1)單源和多源比較

為了說(shuō)明頻率域多源地空電磁方法的優(yōu)勢(shì)，設(shè)計(jì)一個(gè)低阻體模型，針對(duì)垂直磁場(chǎng)分量,比較幾種多源組合模式的電磁場(chǎng)分布.低阻體模型如圖6所示，異常體電阻率為10 Ωm，大小為500 m×500 m×300 m，異常體中心坐標(biāo)為(0 m，0 m，250 m), 圍巖電阻率為100 Ωm.圖6c為非結(jié)構(gòu)化網(wǎng)格剖分圖.為便于對(duì)比多源組合方式，假設(shè)所有的電性源偶極矩為1000 Am，飛行高度30 m.

圖4 均勻半空間3D-Forward、1D-Forward的磁場(chǎng)Hz及對(duì)應(yīng)的誤差曲線(a) 磁場(chǎng)Hz； (b) 相對(duì)誤差曲線.Fig.4 Magnetic field Hz and relative error for homogenous half space from 3D- and 1D-forward modeling(a) Magnetic field Hz; (b) Relative error.

圖5 兩層模型3D-Forward與1D-Forward的磁場(chǎng)Hz及對(duì)應(yīng)的誤差曲線(a) 磁場(chǎng)Hz； (b) 相對(duì)誤差.Fig.5 Magnetic field Hz and relative error of two-layer model from 3D- and 1D-forward modeling(a) Magnetic field Hz； (b) Relative error.

圖7為多個(gè)源的位置示意圖，我們分析幾種組合方式：①單源，①②組合，①③組合，①②③組合，①②③④組合，①⑤⑥組合.各個(gè)源的中心坐標(biāo)已在圖7中給出.

我們首先給出256 Hz各種組合源作用下的均勻半空間電磁響應(yīng).圖8為各種組合源作用下的垂直磁場(chǎng)切片圖，可以看出，靠近源的一側(cè)磁場(chǎng)強(qiáng)，而遠(yuǎn)離源的一側(cè)磁場(chǎng)逐漸衰減，同時(shí)我們也注意到，這幾種組合源都將使得背景場(chǎng)值增大，信號(hào)強(qiáng)度增強(qiáng)，這是多源場(chǎng)之間的相互疊加作用導(dǎo)致的.圖9給出了中心測(cè)線全域視電阻率測(cè)深曲線及相對(duì)誤差，可以看出,全域視電阻率基本接近于背景電阻率100 Ωm，誤差不超過(guò)2.5%，完全滿足計(jì)算精度要求，也說(shuō)明本文全域視電阻率算法的正確性.

圖10給出了單源和多源作用下，含異常體的總場(chǎng)幅值分布切片圖，圖11為相應(yīng)的二次場(chǎng)分布圖，圖12是單源和多源情況下異常幅度的分布圖，圖13給出了單源及多源作用下的全域視電阻率分布圖.從中可以看出，若在單源作用下測(cè)量垂直磁場(chǎng)時(shí)，在靠近源的一側(cè)和遠(yuǎn)離源的一側(cè)由于在界面上電阻率突變，導(dǎo)致界面產(chǎn)生積累電荷，磁場(chǎng)響應(yīng)增強(qiáng)，異常體的中心與異常響應(yīng)的極值并不對(duì)應(yīng)，在靠近源的一側(cè)邊界將會(huì)使得場(chǎng)值增大，而在遠(yuǎn)離源的一側(cè)邊界將會(huì)使場(chǎng)值減小，如圖10a和圖11a所示，而磁場(chǎng)增大將會(huì)使得全域視電阻率增大，因而在靠近源的一側(cè)邊界處出現(xiàn)出高阻極大值，而在遠(yuǎn)離源的一側(cè)邊界處出現(xiàn)了低阻極小值，如圖13a所示.可見(jiàn)，對(duì)于三維地電模型的電磁響應(yīng)分布比較復(fù)雜，異常體的位置和磁場(chǎng)及全域視電阻率形態(tài)并不一定吻合，需要進(jìn)行仔細(xì)分辨.磁場(chǎng)分布不僅與異常體的形態(tài)有關(guān)，也與異常體與源的相對(duì)位置有很大關(guān)系.

圖6 低阻體模型(a) 剖面圖； (b) 平面圖； (c) 非結(jié)構(gòu)化網(wǎng)格剖分圖.紅色代表低阻體，綠色代表背景.Fig.6 Model of a conductive body(a) Cross-section； (b) Plane view； (c) Unstructured grid.Red shows the conductive body and green shows background.

分析①②組合源和①②③組合源的電磁響應(yīng)，發(fā)現(xiàn)總場(chǎng)和二次場(chǎng)的幅度都在一定程度上有所增強(qiáng)，如圖10b、圖10d、圖11b、圖11d所示.這也是多源的優(yōu)勢(shì)之一，可增強(qiáng)觀測(cè)區(qū)域的信號(hào)強(qiáng)度.為了更好的反映異常響應(yīng)分布，需要分析其異常響應(yīng)的相對(duì)幅值，如圖12b、d所示，我們發(fā)現(xiàn)在邊緣處的極大極小值的相對(duì)幅度在減小，這是不同方向上的源產(chǎn)生的場(chǎng)在邊界區(qū)域進(jìn)行疊加的結(jié)果，導(dǎo)致異常場(chǎng)的分布形態(tài)發(fā)生了變化，而圖13b、d視電阻率分布圖也說(shuō)明了這一點(diǎn)，在邊緣處的電阻率極大值在減小，而極小值在增大.圖10c、圖11c是①③組合的總場(chǎng)幅值和二次場(chǎng)大小，這與單源情況下相比，也提高了信號(hào)強(qiáng)度，而從相對(duì)幅度分布圖12c及全域視電阻率圖13c來(lái)看，在異常體的中心出現(xiàn)了低阻圈閉，在邊緣部位出現(xiàn)了高阻響應(yīng).這種源的組合方式在異常形態(tài)分布上優(yōu)于前面分析的三種方式，能更好的反映異常的中心位置，削弱了異常體邊界處的假極值，但是在源的方向上，出現(xiàn)了拉長(zhǎng)的低阻，而在垂直源的方向上出現(xiàn)了高阻.

圖8 256 Hz單源和多源的垂直磁場(chǎng)分布圖(a) ①； (b) ①②； (c) ①③； (d) ①②③； (e) ①②③④； (f) ①⑤⑥.Fig.8 Vertical magnetic field of single and multiple sources at 256 Hz

圖9 中心測(cè)點(diǎn)全域視電阻率及相對(duì)誤差(a) 全域視電阻率； (b) 相對(duì)誤差.Fig.9 Full-domain apparent resistivity and relative error curves of central measuring point(a) Full-domain apparent resistivity; (b) Relative error.

圖10 256 Hz單源和多源的總場(chǎng)幅值分布圖(a) ①； (b) ①②； (c) ①③； (d) ①②③； (e) ①②③④； (f) ①⑤⑥.Fig.10 Total magnetic field Hs at 256 Hz of single and multiple sources

圖11 256 Hz單源和多源的二次場(chǎng)響應(yīng)分布圖(a) ①； (b) ①②； (c) ①③； (d) ①②③； (e) ①②③④； (f) ①⑤⑥.Fig.11 Secondary magnetic field at 256 Hz of single and multiple sources

圖12 256 Hz單源和多源的異常相對(duì)幅度分布圖(a) ①； (b) ①②； (c) ①③； (d) ①②③； (e) ①②③④； (f) ①⑤⑥.Fig.12 Relative percentage of anomaly at 256 Hz of single and multiple sources

圖13 256 Hz單源和多源的全域視電阻率分布圖(a) ①； (b) ①②； (c) ①③； (d) ①②③； (e) ①②③④； (f) ①⑤⑥.Fig.13 Full-domain apparent resistivity at 256 Hz of single and multiple sources

圖14 單源和多源的主軸剖面的廣域視電阻率圖(a)、(g) ①； (b)、(h) ①②； (c)、(i) ①③； (d)、(j) ①②③； (e)、(k) ①②③④； (f)、(l) ①⑤⑥.Fig.14 Full-domain apparent resistivity of single-source and multiple sources spindle sections

圖15 256 Hz單源和多源全域視電阻率、總場(chǎng)幅值和二次場(chǎng)響應(yīng)(a)—(c) ①； (d)—(f) ①②③④； (g)—(i) ①⑤⑥.Fig.15 Full-domain apparent resistivity, total field and secondary filed responses at 256 Hz of single and multiple sources

圖16 單源和多源的主軸剖面的廣域視電阻率圖(a)、(b) ①； (c)、(d) ①②③④； (e)、(f) ①⑤⑥.Fig.16 Full-domain apparent resistivity of spindle sections of single source and multiple sources

從①③組合方式來(lái)看，當(dāng)源關(guān)于異常體的中心出現(xiàn)對(duì)稱分布時(shí)，對(duì)異常體的分辨較好，因而本文還設(shè)計(jì)了①②③④組合和①⑤⑥組合，而①⑤⑥組合三個(gè)源相互呈120°對(duì)稱，從總場(chǎng)幅值和二次場(chǎng)分布可看出，在測(cè)網(wǎng)中心出現(xiàn)了低阻圈閉，且①②③④組合的總場(chǎng)極小值比其余組合方式的極小值更大，說(shuō)明①②③④組合在增強(qiáng)信號(hào)強(qiáng)度方面是最強(qiáng)的，而在二次場(chǎng)分布方面，由于場(chǎng)源和異常體的對(duì)稱性，這兩種方式基本消除了異常體邊界的極值，突出了異常中心的電性特征，而異常相對(duì)幅度和全域視電阻率分布更好的說(shuō)明了這一點(diǎn)，在沒(méi)有異常體的區(qū)域，受到異常體的影響很小，幾乎接近于背景值，而在異常體中心形成了很明顯的低阻異常圈閉，因而這兩種方式都能夠提高磁場(chǎng)的信號(hào)強(qiáng)度，同時(shí)使全域視電阻率分布形態(tài)更接近于真實(shí)的地電分布.若不考慮其他因素，①②③④組合是相對(duì)較優(yōu)的多源組合方式，若考慮野外施工效率，也可選擇①⑤⑥多源組合方式.

接下來(lái)我們分析每個(gè)組合方式的主測(cè)線剖面圖，包括x=0測(cè)線和y=0測(cè)線，如圖14所示，圖14a—f表示單源和多源組合的x=0測(cè)線剖面，圖14g—l表示單源和多源組合的y=0測(cè)線剖面.從中可看出，單源激勵(lì)下，兩條測(cè)線僅平行于源的測(cè)線出現(xiàn)了低阻圈閉，另一條垂直測(cè)線靠近源的一側(cè)呈現(xiàn)高阻，遠(yuǎn)離源的一側(cè)呈現(xiàn)低阻，這和上述切片圖是一致的(圖14a、 g).組合源①②的兩條測(cè)線均沒(méi)有在異常體位置出現(xiàn)低阻圈閉.對(duì)于①③組合，在切片圖中已展現(xiàn)出較寬的低阻圈閉，而兩條中心測(cè)線的剖面圖也能較好的體現(xiàn)異常體的位置，但平行源的測(cè)線低阻圈閉橫向范圍寬(圖14i)，而其x=0測(cè)線在邊界處呈現(xiàn)了高阻響應(yīng)，中心低阻圈閉較小(圖14c).①②③組合的切片圖形與單源形式是類(lèi)似的，因而其剖面圖的形態(tài)也類(lèi)似，只是出現(xiàn)低阻圈閉的測(cè)線相反.①②③④組合和①⑤⑥組合的兩條中心測(cè)線都出現(xiàn)了明顯的低阻圈閉，而且與異常體的位置對(duì)應(yīng)的較好，周?chē)鷩鷰r呈現(xiàn)出背景電阻率，這兩種組合更容易識(shí)別異常體的分布.

3.2.2 高阻體模型

同樣，我們也分析了高阻體的地空三維電磁響應(yīng)分布，圖15a—c為單源作用下的高阻異常體響應(yīng)，可以看出在靠近源的一側(cè)表現(xiàn)出了低阻性質(zhì)，而在遠(yuǎn)離源的一側(cè)表現(xiàn)出了高阻性質(zhì)，這與低阻體的響應(yīng)是相反的.從①⑤⑥組合和①②③④組合的響應(yīng)來(lái)看，異常響應(yīng)呈現(xiàn)高阻圈閉，消除了邊緣的假極值.而圖16給出了相應(yīng)的主軸剖面，從中看出，單源情況下的x=0測(cè)線剖面和y=0測(cè)線剖面的形態(tài)和低阻體單源情況下的剖面形態(tài)是一致的，只是高阻和低阻發(fā)生了倒轉(zhuǎn)，如圖16a、b所示，也可以看出，y=0測(cè)線的高阻圈閉形態(tài)出現(xiàn)誤差，這是因?yàn)榈乜枕憫?yīng)對(duì)高阻體不敏感造成的，產(chǎn)生的異常響應(yīng)和圍巖的電阻率較為接近，因而對(duì)高阻體的靈敏度不如低阻體.從①②③④組合和①⑤⑥組合的中心測(cè)線剖面可以看出，每種組合的兩條中心測(cè)線都出現(xiàn)了高阻圈閉，雖然異常幅度較低，但是能夠很好的反映異常體的位置.

3.2.3 兩個(gè)低阻體模型

建立兩個(gè)大小相同的低阻體，大小均為400 m×500 m×300 m，電阻率都為10 Ωm，關(guān)于y軸對(duì)稱，左邊異常體中心為(-400 m，0 m，250 m),右側(cè)異常體中心為(400 m，0 m，250 m)，背景電阻率為100 Ωm，圖17a為模型剖面圖，圖17b為模型平面圖，圖17c為源的位置示意圖，圖17d為非結(jié)構(gòu)化網(wǎng)格剖分圖.

共設(shè)計(jì)了三種場(chǎng)源的組合方式，在單源①的作用下，靠近源的一側(cè)呈現(xiàn)兩個(gè)了高阻，而遠(yuǎn)離源的一側(cè)呈現(xiàn)了兩個(gè)低阻圈閉，且對(duì)稱分布，如圖18a—c所示，這與前文的結(jié)果是吻合的.而①⑤⑥組合和①②③④組合的響應(yīng)也呈對(duì)稱分布，但是較單源而言，低阻異常的位置同實(shí)際位置更加接近，兩側(cè)呈對(duì)稱分布的高阻異常也已大致趨于背景值，而異常響應(yīng)在兩個(gè)低阻體相互靠近的一側(cè)邊緣出現(xiàn)了低阻極小值，雖然和異常體并沒(méi)有完全吻合，但比單源作用下的響應(yīng)形態(tài)更簡(jiǎn)單，更能反映異常體的電性分布特征.

為了更好的比較，我們給出了不同頻率下兩條中心測(cè)線(x=0和y=0)的垂直磁場(chǎng)變化曲線,如在圖19, 可以發(fā)現(xiàn)，在y=0測(cè)線上，多源激發(fā)下垂直磁場(chǎng)都出現(xiàn)了兩個(gè)低值凹陷，但明顯四源和三源情況下的凹陷位置更接近于異常體的中心，而單源激發(fā)下，在中間位置磁場(chǎng)響應(yīng)卻增大了，與實(shí)際模型電性不一致(圖19e).從剖面圖整體上可以看出，多源的信號(hào)強(qiáng)度明顯強(qiáng)于單源的信號(hào)，這對(duì)于實(shí)際測(cè)量非常有利，在相同噪聲背景下能夠提高信噪比，如圖19a、c、e所示.而對(duì)于x=0測(cè)線，四源和三源在中心位置出現(xiàn)了低阻凹陷，單源在中心位置兩側(cè)則出現(xiàn)了一個(gè)突起和一個(gè)凹陷，這個(gè)與平面圖的異常特征一致，多源激勵(lì)不僅可以增強(qiáng)信號(hào)的強(qiáng)度，而且可以改變目標(biāo)體的異常響應(yīng)特征，如圖19b、d、f所示.

為了方便說(shuō)明，在圖20中給出了512 Hz的中心測(cè)線(x=0和y=0)全域視電阻率變化曲線，對(duì)于y=0測(cè)線，四源和三源均出現(xiàn)了兩個(gè)低阻凹陷，且剛好對(duì)應(yīng)異常體的中心位置，而單源在異常體邊緣處出現(xiàn)了極小值，中間位置則出現(xiàn)了極大值(圖20a).而對(duì)于x=0測(cè)線，四源和三源情況下均出現(xiàn)了低阻凹陷，這是因?yàn)殡m然測(cè)線沒(méi)有經(jīng)過(guò)異常體，但旁側(cè)兩個(gè)低阻異常體在多源激勵(lì)下出現(xiàn)視電阻率低值，而單源激勵(lì)下兩側(cè)分別出現(xiàn)了較強(qiáng)的極大極小值分布，不利于異常體電性特征及平面位置的定性分析(圖20b).

4 結(jié)論

本文采用基于電場(chǎng)雙旋度方程的非結(jié)構(gòu)化節(jié)點(diǎn)有限元方法，成功實(shí)現(xiàn)了多源地空電磁響應(yīng)特征的模擬和分析.通過(guò)Krylov子空間模型降階技術(shù)大大降低了大型稀疏矩陣的維度，然后直接求逆矩陣便可得到傳遞函數(shù)，進(jìn)而求得多頻電磁方程的解.采用均勻半空間和層狀模型，通過(guò)和一維解析解進(jìn)行對(duì)比，驗(yàn)證了本文算法的正確性和有效性.通過(guò)泰勒級(jí)數(shù)展開(kāi)取線性主部的迭代算法計(jì)算多源激勵(lì)的全域視電阻率，并將其應(yīng)用于三維地下介質(zhì)的定性分析解釋.三維地電模型的多源電磁響應(yīng)分析表明：對(duì)于單個(gè)常體，一般單源激勵(lì)會(huì)在靠近源的一側(cè)出現(xiàn)高視電阻率，遠(yuǎn)離源的一側(cè)出現(xiàn)低視電阻率；而對(duì)稱的三源或四源激勵(lì)模式可削弱場(chǎng)在邊界上的響應(yīng)，其全域視電阻率特征相對(duì)簡(jiǎn)單，和地下介質(zhì)電性分布較吻合，這對(duì)于尋找淺地表獨(dú)立的目標(biāo)體，如未爆炸的炮彈、地雷等是有利的.對(duì)于鄰近兩個(gè)低阻異常體，多源激勵(lì)模式不僅能夠增加信號(hào)的強(qiáng)度，而且能夠提高橫向分辨率.實(shí)際地下介質(zhì)的電性變化非常復(fù)雜，多源激勵(lì)的電磁響應(yīng)及全域視電阻率可為定性解釋提供參考，精細(xì)化的解釋將進(jìn)一步依賴于地空三維電磁反演方法的發(fā)展.

圖17 模型示意圖(a) 剖面圖； (b) 平面圖； (c) 源的加載方式； (d) 網(wǎng)格剖分. 紅色代表低阻體，綠色代表背景.Fig.17 Model of two conductive bodies(a) Cross-section； (b) Plane view； (c) Loading manner of source； (d) Unstructured grid. Red shows the conductive body and green shows background.

圖18 256 Hz兩個(gè)低阻體的單源和多源全域視電阻率、總場(chǎng)幅值和二次場(chǎng)響應(yīng)(a)—(c)①; (d)—(f)①⑤⑥; (g)—(i)①②③④.Fig.18 Full-domain apparent resistivity, total and secondary vertical magnetic field at 256 Hz of single and multiple sources

附錄A

(A1)

(A2)

(A3)

(A4)

(A5)

圖19 不同頻率兩條中心測(cè)線的垂直磁場(chǎng)變化曲線圖(a)、(b) ①②③④； (c)、(d) ①⑤⑥； (e)、(f) ①.Fig.19 Vertical magnetic field curves of different frequencies in two central survey lines

圖20 512 Hz下中心測(cè)線的全域視電阻率曲線(a) y=0； (b) x=0.Fig.20 Full-domain apparent resistivity curves at 512 Hz