多源頻率域地空系統(tǒng)三維電磁響應分析

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

1 長安大學地質工程與測繪學院， 西安 710054 2 中南大學地球科學與信息物理學院， 長沙 410083

0 引言

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

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

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

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

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

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

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

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

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

2.1 基于總場的公式

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

(1)

(2)

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

由此便可導出電場雙旋度方程：

(3)

(4)

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

(5)

2.2 有限元分析

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

圖2 非結構化網格剖分Fig.2 Unstructured grid

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

(6)

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

(7)

那么地空多源正演問題轉化為求解該變分方程的極值問題，該方程可看成多元函數，多元函數取極值應當滿足下列方程(徐世浙，1994):

(8)

其中i為節(jié)點編號，e為單元編號.

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

(9)

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

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

(10)

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

2.3 模型降階方法

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

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

(11)

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

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

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

(12)

2.4 垂直磁場的求取

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

(13)

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

(14)

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

3 多源三維電磁響應分析

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

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

為了進一步驗證程序的正確性，我們設計了一個兩層模型，第一層電阻率為200 Ωm，厚度為100 m，第二層電阻率為50 Ωm.第一個源沿x軸正方向放置，電偶源的中心位置為(0 m，-4000 m，0 m)，第二個源沿x軸負方向放置，電偶源的中心坐標為(0 m，2000 m，0 m)，設置測點位置為(1600 m，-2000 m，-100 m).如圖5所示，磁場相對誤差最大不超過3%，說明基于模型降階的有限元計算結果滿足精度要求.

3.2 三維電磁響應分析

3.2.1 低阻體模型

(1)單源和多源比較

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

圖4 均勻半空間3D-Forward、1D-Forward的磁場Hz及對應的誤差曲線(a) 磁場Hz； (b) 相對誤差曲線.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的磁場Hz及對應的誤差曲線(a) 磁場Hz； (b) 相對誤差.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為多個源的位置示意圖，我們分析幾種組合方式：①單源，①②組合，①③組合，①②③組合，①②③④組合，①⑤⑥組合.各個源的中心坐標已在圖7中給出.

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

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

圖6 低阻體模型(a) 剖面圖； (b) 平面圖； (c) 非結構化網格剖分圖.紅色代表低阻體，綠色代表背景.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.

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

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

圖9 中心測點全域視電阻率及相對誤差(a) 全域視電阻率； (b) 相對誤差.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單源和多源的總場幅值分布圖(a) ①； (b) ①②； (c) ①③； (d) ①②③； (e) ①②③④； (f) ①⑤⑥.Fig.10 Total magnetic field Hs at 256 Hz of single and multiple sources

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

圖12 256 Hz單源和多源的異常相對幅度分布圖(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單源和多源全域視電阻率、總場幅值和二次場響應(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

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

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

3.2.2 高阻體模型

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

3.2.3 兩個低阻體模型

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

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

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

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

4 結論

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

圖17 模型示意圖(a) 剖面圖； (b) 平面圖； (c) 源的加載方式； (d) 網格剖分. 紅色代表低阻體，綠色代表背景.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兩個低阻體的單源和多源全域視電阻率、總場幅值和二次場響應(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 不同頻率兩條中心測線的垂直磁場變化曲線圖(a)、(b) ①②③④； (c)、(d) ①⑤⑥； (e)、(f) ①.Fig.19 Vertical magnetic field curves of different frequencies in two central survey lines

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