張博, 殷長春, 劉云鶴, 蔡晶

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院， 長春　130026

起伏地表頻域/時(shí)域航空電磁系統(tǒng)三維正演模擬研究

張博, 殷長春*, 劉云鶴, 蔡晶

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院， 長春130026

摘要由于航空電磁系統(tǒng)具有工作頻率低、時(shí)間延遲短等特點(diǎn)，地形對航空電磁響應(yīng)有很大影響，忽略地形影響會(huì)給航空電磁數(shù)據(jù)解釋造成很大誤差.本文將基于非結(jié)構(gòu)化網(wǎng)格的矢量有限元法應(yīng)用于模擬起伏地表?xiàng)l件下頻域/時(shí)域(FD/TD)三維航空電磁系統(tǒng)響應(yīng).該方法由于采用非結(jié)構(gòu)網(wǎng)格，與傳統(tǒng)的結(jié)構(gòu)化網(wǎng)格電磁正演算法相比，能更好地?cái)M合地形和地下不規(guī)則異常體，提高對不規(guī)則地形和地下介質(zhì)航空電磁響應(yīng)的計(jì)算精度.通過將計(jì)算結(jié)果與半空間模型的半解析解及已發(fā)表的結(jié)果進(jìn)行對比，檢驗(yàn)了本文算法的精度.通過對典型山峰和山谷地形航空電磁響應(yīng)分析對比，總結(jié)了地形對航空電磁響應(yīng)的影響特征.研究結(jié)果對航空電磁地形效應(yīng)的識別和校正具有指導(dǎo)意義.關(guān)鍵詞航空電磁； 三維正演； 地形效應(yīng)； 矢量非結(jié)構(gòu)有限元法

Starting from the frequency-domain Maxwell′s equations, we obtain a vector Helmholtz equation for a stable solution when air is present. We use the Galerkin method to discretize the Helmholtz equation and obtain the final finite-element equations. To accelerate the calculation speed, the source response in a free-space is used as the primary field, while the direct solver is used to deal with the multiple-source problem. After calculating the frequency-domain AEM responses, a Hankel′s transform is applied to obtain the time-domain EM responses.

To check the accuracy of the presented algorithm for AEM topographic modeling,we compare our results with both the semi-analytical results and those from published literatures.After that, we calculate the EM responses for 3D models with only topography and with both topography and anomalous targets.

From profiled AEM responses, we draw similar conclusions to those of Yin (2015), while from area AEM responses, we find that: 1) for frequency-domain AEM systems, the real part of AEM response contains more information to deep earth than the imaginary part; 2) for time-domain AEM system, the magnetic induction dB/dtreveals the underground conductivity distribution better than theBfield. These features provide the theoretical basis for identification and correction of topographic effect from the AEM measurements.

1引言

航空電磁法由于其高效、經(jīng)濟(jì)、且無需地面人員接近勘查區(qū)等特點(diǎn)，在地質(zhì)填圖、油氣勘探、礦產(chǎn)普查及工程和環(huán)境調(diào)查等方面得到了廣泛的應(yīng)用(Wolfgram and Goiden, 2001; Paine and Collins, 2003;Smith et al., 2006; Pfaffhuber et al., 2009;Tan et al., 2009; Yang and Oldenburg, 2013).國內(nèi)外眾多學(xué)者研究表明，地形效應(yīng)對航空電磁響應(yīng)的影響巨大(Liu and Becker, 1992; Annetts et al.，1998；Mitsuhata, 2000; Baba and Seama, 2002; Sasaki and Nakazato, 2003; Nam et al.，2007；劉云鶴和殷長春，2013；蔡晶等，2014).忽略這種影響會(huì)給航空電磁數(shù)據(jù)解釋帶來很大誤差(殷長春等,2015).然而人們對帶地形條件下頻/時(shí)域三維航空電磁響應(yīng)的研究還十分有限.傳統(tǒng)的有限差分、積分方程和基于規(guī)則網(wǎng)格的有限元法已經(jīng)被廣泛應(yīng)用于解決三維航空電磁正演模擬問題(Newman and Alumbaugh，1995；Avdeev et al.，1998)，然而這些方法往往采用矩形或長方體對模型進(jìn)行剖分，在處理不規(guī)則地形或地下復(fù)雜異常體時(shí)往往需要將網(wǎng)格剖分得很細(xì)，浪費(fèi)計(jì)算資源，且擬合地形效果不佳.近年，一些學(xué)者使用非結(jié)構(gòu)網(wǎng)格對電磁正演問題進(jìn)行的研究主要集中在海洋電磁領(lǐng)域(Um et al., 2010;Ren et al.，2014; Fu et al., 2015)，采用非結(jié)構(gòu)化網(wǎng)格對起伏地表?xiàng)l件下三維頻/時(shí)域航空電磁響應(yīng)進(jìn)行模擬的學(xué)者較少.由于航空電磁已經(jīng)越來越廣泛地應(yīng)用于地球物理勘探，地形對航空電磁響應(yīng)的影響十分巨大，本文將基于非結(jié)構(gòu)網(wǎng)格的矢量有限元法應(yīng)用于帶地形條件下三維頻/時(shí)域航空電磁系統(tǒng)響應(yīng)正演模擬.

由于非結(jié)構(gòu)網(wǎng)格能夠有效剖分任意復(fù)雜地質(zhì)模型，其在電磁正演模擬中的應(yīng)用比傳統(tǒng)的結(jié)構(gòu)化網(wǎng)格和二叉樹網(wǎng)格更加廣泛.本文采用體積和質(zhì)量約束的Delauany四面體剖分技術(shù)對計(jì)算區(qū)域進(jìn)行網(wǎng)格剖分.由于對Maxwell方程進(jìn)行直接求解會(huì)出現(xiàn)解不穩(wěn)定的問題，我們使用電場的Helmholtz方程進(jìn)行問題求解.一次/二次場分離算法與直接計(jì)算總場算法相比具有更高的計(jì)算精度和更好的穩(wěn)定性，因此本文采用分離算法對航空電磁響應(yīng)進(jìn)行模擬.我們使用Galerkin方法對Helmholtz方程進(jìn)行離散，得到三維航空電磁正演問題對應(yīng)的數(shù)學(xué)方程.由于航空電磁具有多源性，我們采用直接分解技術(shù)對矩陣方程進(jìn)行求解.在得到不同頻率的航空電磁響應(yīng)后，我們對頻率域響應(yīng)進(jìn)行漢克爾變換，從而得到時(shí)間域航空電磁響應(yīng).

我們首先通過對比本文計(jì)算結(jié)果與半空間模型的解析解以及其他已發(fā)表的結(jié)果檢驗(yàn)算法的精度.然后，我們針對典型山峰和山谷地形模型計(jì)算航空電磁響應(yīng)，分析純地形模型航空電磁響應(yīng)特征.為了進(jìn)一步分析地形對航空電磁系統(tǒng)響應(yīng)的影響，我們還計(jì)算和分析了地形模型中埋藏異常體模型的電磁響應(yīng).

2正演算法

電磁場在發(fā)射源附近變化非常劇烈，使用直接計(jì)算總場方法進(jìn)行正演模擬時(shí)，為保證正演精度，必須將源附近的網(wǎng)格劃分得很細(xì)，浪費(fèi)大量的計(jì)算資源.一次場二次場分離算法將異常體看成虛擬源以代替真實(shí)源，從而避免了真實(shí)發(fā)射源附近場值變化劇烈的問題.然而，理論研究表明，一次場二次場分離方法不適合處理背景場難以獲得的復(fù)雜地電模型(Mitsuhata, 2000; Jahandari and Farquharson, 2014)，這是由于獲得背景場耗費(fèi)的計(jì)算資源十分巨大.此外，該方法必須將起伏地表看成異常體，當(dāng)發(fā)射源或接收點(diǎn)位于起伏地表附近時(shí)，發(fā)射源和接收點(diǎn)附近的起伏地表必須被劃分成很密的網(wǎng)格以保證求解精度.由于航空電磁收發(fā)裝置距離地面有一定的飛行高度，因此無需將起伏地表網(wǎng)格劃分得很密.另外，由于該方法對于任意復(fù)雜模型均可使用源在全空間的響應(yīng)作為背景場，能大幅度減少背景場計(jì)算工作量.從以上分析可以看出，一次場二次場分離算法在處理航空電磁正演問題時(shí)保留了計(jì)算精度高、網(wǎng)格劃分少的優(yōu)點(diǎn)；同時(shí)避免了復(fù)雜模型背景場計(jì)算量大，對發(fā)射源和接收點(diǎn)附近的起伏地表網(wǎng)格要求剖分密等缺點(diǎn).本文采用一次場二次場分離算法對三維航空電磁響應(yīng)進(jìn)行模擬.二次場滿足如下形式的Maxwell方程：

(1)

(2)

式中σ是電導(dǎo)率，σs=σ-σp(σp是背景電導(dǎo)率)，Ep是背景電場，Es是待求解的二次場，μ是磁導(dǎo)率.直接求解(1)式和(2)式會(huì)出現(xiàn)解不穩(wěn)定問題.由此，三維航空電磁通常對Helmholtz方程進(jìn)行求解.從(1)和(2)式我們可以得到電、磁場滿足的Helmholtz方程.我們利用電場的Helmholtz方程進(jìn)行求解，即

(3)

為確保電場具有唯一解，本文對計(jì)算區(qū)域的外邊界添加了Dirichlet邊界條件：

(4)

式中Ω表示計(jì)算區(qū)域外邊界，g表示二次場在邊界處的場值.本文令計(jì)算區(qū)域的邊界遠(yuǎn)離異常體，因此可以假設(shè)g=0.(3)和(4)式構(gòu)成了二次電場對應(yīng)的邊值問題.

本文使用有限元法求解邊值問題包括以下幾個(gè)步驟：

1) 將計(jì)算區(qū)域離散成小單元；

2) 計(jì)算(3)式與形函數(shù)N的L2內(nèi)積以得到邊值問題對應(yīng)的泛函形式(L2內(nèi)積的定義見附錄B)；

3) 計(jì)算上述邊值問題的泛函，得到整體剛度矩陣；

4) 求解有限元方程，得到問題的解.

圖1　矢量有限元四面體單元內(nèi)電場分布Fig.1　Electrical fields in a tetrahedron element for vector FEM

我們首先使用Delaunay方法對計(jì)算區(qū)域進(jìn)行四面體網(wǎng)格剖分.Delaunay網(wǎng)格的具體性質(zhì)及剖分方法見附錄A.為了保證網(wǎng)格剖分質(zhì)量，對網(wǎng)格剖分添加了體積和質(zhì)量約束.使用節(jié)點(diǎn)有限元法進(jìn)行電磁場數(shù)值模擬存在嚴(yán)重的缺陷.首先,由于未強(qiáng)加散度條件，節(jié)點(diǎn)有限元法求解出來的場可能出現(xiàn)非物理解；其次，由于節(jié)點(diǎn)有限元法在介質(zhì)分界面處求解出的場不具有連續(xù)性，在介質(zhì)分界面處強(qiáng)加邊界條件十分不便.矢量有限單元法由于將電場賦于棱邊上, 并且其基函數(shù)在每個(gè)單元內(nèi)自然滿足散度條件，確保了電場在介質(zhì)分界面處的連續(xù)性，從而克服了節(jié)點(diǎn)有限元的缺陷.本文選擇的矢量有限單元每個(gè)單元中電場在棱邊上的分布情況如圖1所示.

對(3)式與基函數(shù)在每個(gè)小單元進(jìn)行L2內(nèi)積，并使用格林公式和分部積分，可以得到邊值問題對應(yīng)的變分形式.對變分形式的邊值問題使用有限元方法進(jìn)行離散，并將所有小單元對應(yīng)的系數(shù)矩陣進(jìn)行組合，我們可以得到形式如下的整體矩陣方程：

(5)

式中Es是棱邊上的矢量電場，系數(shù)矩陣K是由每個(gè)單元對應(yīng)的小矩陣K1和K2組成的對稱稀疏非Hermitian矩陣，右端項(xiàng)向量b由每個(gè)單元對應(yīng)的源項(xiàng)向量be=K3Ep組成.公式(5)的推導(dǎo)過程見附錄B.

由于航空電磁的移動(dòng)平臺特征具有多源性，傳統(tǒng)的迭代求解算法處理該問題時(shí)需要對每個(gè)源進(jìn)行一次單獨(dú)正演，計(jì)算速度慢.本文采用直接分解算法對(5)式進(jìn)行求解.該方法針對模型不變的情況只需一次分解就能得到所有源對應(yīng)的響應(yīng)(Oldenburgetal., 2013)，求解效率較之于傳統(tǒng)的迭代算法有很大提高.由式(5)求得棱邊上的電場Es后，空間任意位置的磁場可由Faraday′s定律得到：

(6)

由(5)和(6)式計(jì)算出多個(gè)頻率對應(yīng)的航空電磁系統(tǒng)響應(yīng)，針對本文發(fā)射的階躍波信號，利用正弦變換(等價(jià)于0.5階漢克爾變換)對頻率域航空電磁響應(yīng)進(jìn)行反傅里葉變換：

(7)

(8)

式中t表示時(shí)間，J表示貝塞爾函數(shù)，ω為角頻率.經(jīng)過以上變換，我們可以由頻率域航空電磁響應(yīng)得到時(shí)間域航空電磁響應(yīng).

3精度驗(yàn)證

為了檢驗(yàn)本文算法對三維航空電磁響應(yīng)的模擬精度，我們將本文計(jì)算結(jié)果與Sasaki和Nakazato(2003)使用有限差分計(jì)算的起伏地表頻率域航空電磁響應(yīng)進(jìn)行對比；對于時(shí)間域航空電磁響應(yīng)，我們將本文的結(jié)果與半空間模型的半解析解進(jìn)行對比.Sasaki和Nakazato(2003)計(jì)算的梯形山峰模型(簡稱模型一)及其網(wǎng)格劃分情況如圖2a所示.山峰高50m，頂部和底部分別寬20m和220m；模型的電阻率為100Ωm；收發(fā)裝置采用水平共面裝置，收發(fā)距10m，發(fā)射頻率為16kHz，飛行高度固定為距離起伏地表30m.圖2b給出了本文模擬結(jié)果與Sasaki和Nakazato(2003)模擬結(jié)果之間的對比，圖2c給出兩者之間的相對誤差.從圖中可以看出，兩者吻合較好，最大相對誤差不超過4%.

圖2　本文模擬結(jié)果與Sasaki和Nakazato(2003)模擬結(jié)果對比(a) 山峰模型； (b) 模擬結(jié)果對比； (c) 相對誤差.Fig.2　Comparison of results from this paper with those from Sasaki and Nakazato(2003) (a) A trapezoid hill model; (b) Comparison of the FE solution of this paper with those from Sasaki and Nakazato(2003); (c) Relative errors.

圖3　本文模擬結(jié)果與半解析解對比(a) Bz模擬結(jié)果對比； (b) Bz相對誤差曲線； (c)dBz/dt模擬結(jié)果對比； (d) dBz/dt相對誤差曲線.Fig.3　Comparison of FE results from this paper with the semi-analytical results from Yin et al. (2013)(a) Comparison of time-domain Bz; (b) Relative errors for Bz; (c) Comparison of dBz/dt; (d) Relative errors for dBz/dt.

圖5　純地形模型頻率域航空電磁測線響應(yīng)(a) 模型二響應(yīng)； (b) 模型三響應(yīng).Fig.5　Frequency-domain AEM responses in profile(a) EM responses for model 2； and (b) for model 3.

由于帶地形模型時(shí)間域航空電磁響應(yīng)發(fā)表的正演模擬結(jié)果較少，本文將計(jì)算結(jié)果與半空間模型的半解析解(殷長春等, 2013)進(jìn)行對比.收發(fā)裝置參考FUGRO系統(tǒng)，參數(shù)設(shè)置如下：發(fā)射高度30m，接收高度50m，發(fā)射、接收線圈的水平距離為10m，發(fā)射偶極矩為615000Am2，發(fā)射波形為階躍波(下文所有時(shí)間域系統(tǒng)參數(shù)均相同).正演模擬結(jié)果對比如圖3a和3c所示，相對誤差如圖3b和3d所示.從圖中可以看出，B場最大誤差不超過2%，dB/dt僅個(gè)別點(diǎn)相對誤差在3%附近，整體相對誤差不超過3%.由此我們得出結(jié)論：本文的正演模擬方法無論對頻率域還是時(shí)間域均具有較高的計(jì)算精度.

4地形影響特征分析

下面我們通過模擬山峰和山谷地形頻率/時(shí)間域三維航空電磁響應(yīng)，分析地形對航空電磁響應(yīng)的影響特征.山峰和山谷地形模型如圖4a和4b所示(簡稱模型二和模型三).山峰高度和山谷深度均為20m，在頂部和底部的寬度分別為60m和220m，模型的電阻率為100Ωm.對于頻率域航空系統(tǒng)，本文模擬了水平共面(HCP)裝置380Hz, 1600Hz, 6300Hz, 25kHz, 120kHz五個(gè)頻率的電磁響應(yīng)，裝置距離地表30m，收發(fā)距10m(下文所有頻率域系統(tǒng)均與此相同).對于時(shí)間域航空電磁響應(yīng)，本文模擬了FUGRO系統(tǒng)在10-5～10-2s之間12個(gè)時(shí)間道的航空電磁響應(yīng).由于在低頻段和晚期，航空電磁正演響應(yīng)受地形影響很小，本文的正演響應(yīng)結(jié)果平面分布僅給出120kHz(高頻)和1600Hz(中低頻)以及10-5s(早期時(shí)間道)和3.2×10-4s(中晚期道)的結(jié)果，而對測線正演響應(yīng)我們給出了全部12個(gè)時(shí)間道以及5個(gè)頻率的計(jì)算結(jié)果.圖5a和5b分別給出了模型二和模型三對應(yīng)5個(gè)頻率的電磁響應(yīng)信號，測線位于模型中心剖面；圖6和圖7分別給出了模型二和模型三頻率域航空電磁響應(yīng)的平面分布；圖8給出了模型二和模型三12個(gè)時(shí)間道電磁信號的響應(yīng)，測線位于模型中心剖面；圖9和10分別給出了模型二和模型三時(shí)間域航空電磁B和dB/dt的平面響應(yīng).

從圖5、8的模擬結(jié)果和響應(yīng)特征分析我們可以得出如下結(jié)論： (1) 地形對頻率域和時(shí)間域航空電磁響應(yīng)均有較大影響，航空電磁數(shù)據(jù)處理解釋不能忽略地形效應(yīng)； (2) 頻率越高或時(shí)間道越早，電磁響應(yīng)受地形的影響越大.這是因?yàn)楦哳l段電磁信號穿透能力較弱，反映的信息主要來自距離航空電磁收發(fā)裝置較近的地表介質(zhì)；而低頻段電磁信號穿透能力較強(qiáng)，反映的信息主要來自距離航空電磁收發(fā)裝置較遠(yuǎn)的地下深部介質(zhì).由此也揭示了地形影響主要出現(xiàn)在高頻段電磁信號中.與之相對應(yīng)，早期道電磁信號反映的信息來自距離航空電磁收發(fā)裝置較近的地表介質(zhì)，而晚期道信號反映的信息主要來自離航空電磁收發(fā)裝置較遠(yuǎn)的地下深部介質(zhì).因此，地形對早期時(shí)間道電磁信號影響較大； (3) 帶地形航空電磁響應(yīng)的曲線形態(tài)與地形之間存在鏡像關(guān)系.在地形的拐點(diǎn)位置，電磁響應(yīng)顯示出劇烈的變化.從圖6、7、9、10的模擬結(jié)果和響應(yīng)特征分析我們可以進(jìn)一步看出： (1) 對于同一頻率的電磁信號，航空電磁響應(yīng)的實(shí)部受地形影響小于虛部，而虛部則能夠更好地刻畫地形的形態(tài)，對精細(xì)結(jié)構(gòu)的反應(yīng)能力更強(qiáng)； (2) 對于同一時(shí)間道的電磁信號，B場受到地形的影響小于dB/dt場，而dB/dt場能夠更好地刻畫地形的形態(tài)，對精細(xì)結(jié)構(gòu)的反映能力更強(qiáng).這是由于磁場B是磁感應(yīng)dB/dt的積分，在一定程度上平滑了異常電磁響應(yīng)，使得地下(包括地形)精細(xì)結(jié)構(gòu)的刻畫能力受到減弱.

5地形對異常體響應(yīng)的影響特征

為分析地形對異常體響應(yīng)的影響，我們設(shè)計(jì)了山峰和山谷地形下埋藏低阻異常體的模型(簡稱模型四和模型五).模型四和模型五如圖11a和11b所示.山峰和山谷模型參數(shù)與模型二和模型三參數(shù)相同，異常體埋在距離地表30m的位置，大小40m×

圖6　模型二頻率域航空電磁響應(yīng)平面分布

圖7　模型三頻率域航空電磁響應(yīng)平面分布

圖8　時(shí)間域航空電磁測線響應(yīng)(a) 模型二Bz響應(yīng)； (b) 模型三Bz響應(yīng)； (c) 模型二dBz/dt響應(yīng)； (d) 模型三dBz/dt響應(yīng).Fig.8　Time-domain AEM responses in profile(a) Bz for model 2; (b) Bz for model 3; (c) dBz/dt for model 2; (d)dBz/dt for model 3.

40 m×30 m，電阻率為1 Ωm，圍巖電阻率為100 Ωm.圖12a和12b分別給出模型四和五5個(gè)頻率的電磁響應(yīng)信號.測線位于模型中心剖面，而圖13和14分別給出了模型四和模型五頻率域航空電磁的平面響應(yīng).圖15給出了模型四和模型五12個(gè)時(shí)間道的電磁響應(yīng)，測線位于模型中心剖面，而圖16和17分別給出了模型四和模型五時(shí)間域航空電磁B和dB/dt的平面響應(yīng).

從圖12、15的模擬結(jié)果和響應(yīng)特征分析我們可以得到如下結(jié)論：純地形與帶地形異常體的響應(yīng)曲線在高頻段和早期時(shí)間道差異較?。蝗欢?，隨著頻率的降低或時(shí)間的延長，兩者差異逐漸變大，地下電性不均勻體的影響突出表現(xiàn)出來.這從物理上很容易解釋：在高頻段和早期時(shí)間道，電磁信號穿透淺，主要反映地表信息，此時(shí)地形響應(yīng)占主導(dǎo)地位，而異常體的響應(yīng)微弱、難以識別；隨著頻率降低，或時(shí)間向晚期道推移，電磁信號逐步穿透地下介質(zhì)，地表和地形響應(yīng)減弱，而異常體的響應(yīng)增強(qiáng).由此，我們可以得出結(jié)論，地形效應(yīng)主要體現(xiàn)在頻域電磁信號的高頻段或時(shí)域電磁信號的早期道，而有一定埋深目標(biāo)體的電磁響應(yīng)主要體現(xiàn)在中低頻和中晚期時(shí)間道(取決于良導(dǎo)體導(dǎo)電性和埋深).

從圖13、14、16、17中的模擬結(jié)果和響應(yīng)特征分析我們可以進(jìn)一步看出： (1) 對于頻率域結(jié)果，高頻信號的虛部刻畫地形細(xì)節(jié)的能力強(qiáng)于實(shí)部，而低頻信號的實(shí)部對具有一定埋深的異常體反映能力更強(qiáng).以上現(xiàn)象表明，對于本文設(shè)計(jì)的模型，無論是高頻信號還是中低頻信號，頻率域響應(yīng)的實(shí)部包含的深部信息都較虛部豐富，而虛部則包含更多的淺層信息； (2) 對于時(shí)間域響應(yīng)，無論是早期時(shí)間道信號刻畫地形精細(xì)結(jié)構(gòu)的能力，還是晚期時(shí)間道信號反映異常體的能力，dB/dt都強(qiáng)于B.這一現(xiàn)象再次說明dB/dt能夠更好地刻畫介質(zhì)的電性變化特征.

圖9　模型二時(shí)間域航空電磁響應(yīng)平面分布(a) 10-5s的Bz響應(yīng)； (b) 3.2×10-4s的Bz響應(yīng)； (c) 10-5s的dBz/dt響應(yīng)； (d) 3.2×10-4s的dBz/dt響應(yīng).Fig.9　Plane view of time-domain AEM response for model 2(a) Bz field at 10-5s; (b) Bz field at 3.2×10-4s; (c) dBz/dt at 10-5s; (d) dBz/dt at 3.2×10-4s.

圖10　模型三時(shí)間域航空電磁響應(yīng)平面分布(a) 10-5s Bz響應(yīng)； (b) 3.2×10-4s Bz響應(yīng)； (c) 10-5s dBz/dt響應(yīng)； (d) 3.2×10-4s dBz/dt響應(yīng).Fig.10　Plane view of time-domain AEM response for model 3(a) Bz field at 10-5s; (b) Bz field at 3.2×10-4s; (c) dBz/dt at 10-5s; (d) dBz/dt at 3.2×10-4s.

圖11　帶地形和異常體模型(a) 山峰模型(模型四)； (b) 山谷模型(模型五).Fig.11　An abnormal body under topography(a) An hill model (model 4); (b) A valley model (model 5).

圖12　帶地形和異常體模型頻率域航空電磁響應(yīng)(a) 模型四響應(yīng)； (b) 模型五響應(yīng).Fig.12　Frequency-domain AEM responses in profile for (a) model 4 and (b) model 5

6結(jié)論

本文基于有限元方法能夠模擬復(fù)雜模型，而非結(jié)構(gòu)網(wǎng)格可以很好地模擬起伏地表的特點(diǎn)，成功將非結(jié)構(gòu)矢量有限元法應(yīng)用于模擬帶地形條件下三維航空電磁響應(yīng)正演模擬.通過與均勻半空間的半解析解及已發(fā)表模型結(jié)果進(jìn)行對比，驗(yàn)證了本文算法具有較高的精度.通過選擇使用源在全空間的解析解作為背景場，并采用直接分解求解器對矩陣方程進(jìn)行求解，大大提高了正演計(jì)算速度.本文研究成果一定程度上為有效識別地形效應(yīng)和異常體響應(yīng)提供可能，對航空電磁數(shù)據(jù)解釋在一定程度上提供了理論支撐.

通過大量理論模型的計(jì)算結(jié)果，我們得出如下結(jié)論：

(1) 對于頻率域航空電磁系統(tǒng)，高頻電磁信號虛部刻畫地形細(xì)節(jié)的能力強(qiáng)于實(shí)部，而低頻電磁信號實(shí)部反映具有一定埋深異常體的能力更強(qiáng).由此我們可以合理推測，無論高頻信號還是中低頻信號，頻率域航空電磁響應(yīng)的實(shí)部包含的深部信息都較虛部豐富；

(2) 對于時(shí)間域航空電磁系統(tǒng)，早期時(shí)間道磁感應(yīng)dB/dt刻畫起伏地表細(xì)節(jié)的能力強(qiáng)于磁場B，而在晚期時(shí)間道dB/dt對異常體的反映能力較之于磁場B更強(qiáng).因此，dB/dt對地質(zhì)體電性分布變化更為敏感，在刻畫異常體精細(xì)結(jié)構(gòu)方面更具優(yōu)勢.

(3) 然而，由于地形和異常體響應(yīng)處于相同數(shù)量級，為有效實(shí)現(xiàn)航空電磁數(shù)據(jù)反演解釋，對地形影響進(jìn)行模擬和校正十分必要.這將是我們未來的研究內(nèi)容.

圖13　模型四頻率域航空電磁響應(yīng)平面分布

圖14　模型五頻率域航空電磁響應(yīng)平面分布

圖15　時(shí)間域航空電磁測線響應(yīng)(a) 模型四Bz響應(yīng)； (b) 模型五Bz響應(yīng)； (c) 模型四dBz/dt響應(yīng)； (d) 模型五dBz/dt響應(yīng).Fig.15　Time-domain AEM responses in profile(a) Bz for model 4; (b) Bz for model 5; (c) dBz/dt for model 4; (b) dBz/dt for model 5.

附錄ADelaunay非結(jié)構(gòu)網(wǎng)格性質(zhì)及剖分方法

Delaunay四面體剖分方法能夠得到高質(zhì)量的四面體網(wǎng)格，目前已經(jīng)成為電磁數(shù)值模擬中常用的一種網(wǎng)格剖分方法.Delaunay四面體剖分方法剖分網(wǎng)格具有如下性質(zhì)：

(1) 空圓特性，即確保任意Delaunay四面體網(wǎng)格中任意四面體外接球范圍內(nèi)不會(huì)有其他點(diǎn)存在；

(2) 最大化最小角特性，即在散點(diǎn)集可能形成的四面體網(wǎng)格中，Delaunay四面體剖分所形成的四面體的最小角最大.性質(zhì)1確保了劃分網(wǎng)格的唯一性，并為網(wǎng)格剖分提供了一種依據(jù)，性質(zhì)2確保了劃分網(wǎng)格具有較高的質(zhì)量.圖A1比較了四個(gè)節(jié)點(diǎn)情況下Delaunay和非Delaunay網(wǎng)格剖分結(jié)果.圖A1(a)中三角形BCD的外接圓將點(diǎn)A包含在內(nèi)，因此剖分結(jié)果為非Delaunay網(wǎng)格；圖A1(b)中三角形ABC和ACD的外接圓未將任何不屬于三角形的節(jié)點(diǎn)包含在內(nèi)，剖分結(jié)果為Delaunay網(wǎng)格.從圖中我們可以看出Delaunay網(wǎng)格剖分比非Delaunay網(wǎng)格剖分具有更好的質(zhì)量.

圖A1　二維網(wǎng)格剖分示例(a) 非Delaunay網(wǎng)格； (b) Delaunay網(wǎng)格.Fig.A1　An example for 2D grids(a) General grids； (b) Delaunay grids.

圖16　模型四時(shí)間域航空電磁響應(yīng)平面分布(a) 10-5s的Bz響應(yīng)； (b) 3.2×10-4s的Bz響應(yīng)； (c) 10-5s的dBz/dt響應(yīng)； (d) 3.2×10-4s的dBz/dt響應(yīng).Fig.16　Plane view of time-domain AEM responses for model 4(a) Bz at 10-5s; (b) Bz at 3.2×10-4s; (c) dBz/dt at 10-5s; (d) dBz/dt at 3.2×10-4s.

圖17　模型五時(shí)間域航空電磁響應(yīng)平面分布(a) 10-5s的Bz響應(yīng)； (b) 3.2×10-4s的Bz響應(yīng)； (c) 10-5s的dBz/dt響應(yīng)； (d) 3.2×10-4s的dBz/dt響應(yīng).Fig.17　Plane view of time-domain AEM responses for model 5(a) Bz at 10-5s; (b) Bz at 3.2×10-4s; (c) dBz/dt at 10-5s; (d) dBz/dt at 3.2×10-4s.

基于性質(zhì)1，我們給出一種Delaunay網(wǎng)格剖分方法如下(以二維情況為例)：

1) 構(gòu)造一個(gè)大型輔助四邊形，將空間所有節(jié)點(diǎn)包含在內(nèi)；

2) 連接輔助四邊形的一條對角線，形成兩個(gè)大型三角形；

3) 將節(jié)點(diǎn)逐一插入輔助四邊形內(nèi)，若新節(jié)點(diǎn)的插入使某個(gè)三角形不再遵守性質(zhì)1，則記錄這個(gè)三角形，所有上述不遵守性質(zhì)1的三角形構(gòu)成了Delaunay腔；

4) 刪除Delaunay腔內(nèi)的已有三角形，將新插入的節(jié)點(diǎn)與Delaunay腔的頂點(diǎn)連接，得到腔內(nèi)重新劃分的網(wǎng)格；

5) 重復(fù)3—4步驟，直到所有節(jié)點(diǎn)被插入四邊形內(nèi)；

6) 刪除與大輔助四邊形有關(guān)的三角形單元，完成網(wǎng)格剖分.

附錄B單元?jiǎng)偠染仃囃茖?dǎo)

為方便下文推導(dǎo)，重寫公式(3)：

(B1)

對(B1)式與基函數(shù)在每個(gè)小單元進(jìn)行L2內(nèi)積：

(B2)

并使用格林公式和分部積分，可以得到邊值問題對應(yīng)的變分形式：

?ΩiωμσsNe·EpdΩ，

(B3)

(B4)

(B5)

式中L為標(biāo)量有限元基函數(shù)，l表示棱邊長度.從式(B5)中我們可以看出，矢量基函數(shù)自然滿足散度為0條件.將(B4)式中的電場插值函數(shù)代入(B3)式，對于每個(gè)四面體單元我們可以得到

K1Es+K2Es=K3Ep，

(B6)

式中

(B7)

(B8)

(B9)

(B10)

將所有小單元對應(yīng)的系數(shù)矩陣進(jìn)行組合，我們可以得到形式如下的整體矩陣方程：

K Es=b.

(B11)

References

Annetts D, Sugeng F, Raiche A. 1998. The effect of topography on airborne electromagnetic response.∥ SEG Expanded Abstracts. New Orleans, Louisiana.Avdeev D B, Kuvshinov A V, Pankratov O V, et al. 1998. Three-dimensional frequency-domain modeling of airborne electromagnetic responses.ExplorationGeophysics, 29(2): 111-119. Baba K, Seama N. 2002. A new technique for the incorporation of seafloor topography in electromagnetic modelling.GeophysicalJournalInternational, 150(2): 392-402,doi:10.1046/j.1365-246X.2002.01673.x. Cai J, Qi Y F, Yin C C. 2014. Weighted Laterally-constrained inversion of frequency-domain airborne EM data.ChineseJournalofGeophysics(in Chinese), 57(3): 953-960, doi: 10.6038/cjg2014032.

Fu H H, Wang Y Q, Um E S, et al. 2015. A parallel finite-element time-domain method for transient electromagnetic simulation.Geophysics, 80(4): E231-E224, doi: 10.1190/geo2014-0067.1. Jahandari H, Farquharson C G. 2014. A finite-volume solution to the geophysical electromagnetic forward problem using unstructured grids.Geophysics, 79(6): E287-E302, doi: 10.1190/geo2013-0312.1. Jin J M. 1998. Electromagnetic Finite-Element Method (in Chinese). Xi′an: Xidian University Press. Liu G M, Becker A. 1992. Evaluation of terrain effects in AEM surveys using the boundary element method.Geophysics, 57(2): 272-278, doi: 10.1190/1.1443240.

Liu Y H, Yin C C. 2013. 3D inversion for frequency-domain HEM data.ChineseJ.Geophys. (in Chinese), 56(12): 4278-4287, doi: 10.6038/cjg20131230.

Mitsuhata Y. 2000. 2-D electromagnetic modeling by finite-element method with a dipole source and topography.Geophysics, 65(2): 465-475, doi: 10.1190/1.1444740.

Nam M J, Kim H J, Song Y, et al. 2007. 3D magnetotelluric modelling including surface topography.GeophysicalProspecting, 55(2): 277-287, doi: 10.1111/j.1365-2478.2007.00614.x. Newman G A, Alumbaugh D L. 1995. Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences.GeophysicalProspecting, 43(8): 1021-1042.

Oldenburg D W, Haber E, Shekhtman R. 2013. Three dimensional inversion of multisource time domain electromagnetic data.Geophysics, 78(1): E47-E57, doi: 10.1190/geo2012-0131.1.

Paine J G, Collins E W. 2003. Applying airborne electromagnetic induction in groundwater salinization and resource studies, West Texas.∥ Symposium on the Application of Geophysics to Engineering and Environmental Problems. San Antonio, Texas, U.S.A., 722-738. Pfaffhuber A A, Monstad S, Rudd J. 2009. Airborne electromagnetic hydrocarbon mapping in Mozambique.ExplorationGeophysics, 40(3): 237-245, doi: 10.1071/ASEG2009ab040. Ren Z Y, Kalscheuer T, Greenhalgh S, et al. 2014. A finite-element-based domain-decomposition approach for plane wave 3D electromagnetic modeling.Geophysics, 79(6): E255-E268, doi: 10.1190/geo2013-0376.1.

Sasaki Y, Nakazato H. 2003. Topographic effects in frequency-domain helicopter-borne electromagnetics.ExplorationGeophysics, 34(2): 24-28, doi: 10.1071/EG03024. Smith B D, Thamke J N, Cain M J, et al. 2006. Helicopter Electromagnetic and Magnetic Survey Maps and Data, East Poplar oil field area, August 2004, Fort Peck Indian Reservation, northeastern Montana. Open-File Report.Tan K P, Munday T, Halas L, et al. 2009. Utilizing airborne electromagnetic data to map groundwater salinity and salt store

at Chowilla, SA.∥ASEG Extended Abstracts 2009: 20th Geophysical Conference. ASEG, 1-6.Um E S, Harris J M, Alumbaugh D L. 2010. 3D time-domain simulation of electromagnetic diffusion phenomena: A finite-element electric-field approach.Geophysics, 75(4): F115-F126, doi: 10.1190/1.3473694. Wolfgram P, Goiden H. 2001. Airborne EM applied to sulphide nickel-examples and analysis.ExplorationGeophysics, 32(4): 136-140.

Yang D K, Oldenburg D W. 2013. 3D conductivity models of Lalor Lake VMS deposit from ground loop and airborne EM data sets.∥23rd International Geophysical Conference and Exhibition. Melbourne, Australia: ASEG, 1-4. Yin C C, Huang W, Ben F. 2013. The full-time electromagnetic modeling for time-domain airborne electromagnetic systems.ChineseJ.Geophys. (in Chinese),56(9):3153-3162,doi:10.6038/cjg20130928. Yin C C, Zhang B, Liu Y H, et al. 2015. 2.5-D forward modeling of the time-domain airborne EM system in areas with topographic relief.ChineseJ.Geophys. (in Chinese), 58(4): 1411-1424, doi: 10.6038/cjg20150427.

附中文參考文獻(xiàn)蔡晶, 齊彥福, 殷長春. 2014. 頻率域航空電磁數(shù)據(jù)的加權(quán)橫向約束反演. 地球物理學(xué)報(bào), 57(3): 953-960, doi: 10.6038/cjg2014032. 金建銘. 1998. 電磁場有限元方法. 西安: 西安電子科技大學(xué)出版社. 劉云鶴, 殷長春. 2013. 三維頻率域航空電磁反演研究. 地球物理學(xué)報(bào), 56(12): 4278-4287, doi: 10.6038/cjg20131230. 殷長春, 黃威, 賁放. 2013. 時(shí)間域航空電磁系統(tǒng)瞬變?nèi)珪r(shí)響應(yīng)正演模擬. 地球物理學(xué)報(bào), 56(9): 3153-3162, doi: 10.6038/cjg20130928. 殷長春, 張博, 劉云鶴等. 2015. 2.5維起伏地表?xiàng)l件下時(shí)間域航空電磁正演模擬. 地球物理學(xué)報(bào), 58(4): 1411-1424, doi: 10.6038/cjg20150427.

(本文編輯胡素芳)

3D modeling on topographic effect for frequency-/time-domain airborne EM systems

ZHANG Bo, YIN Chang-Chun*, LIU Yun-He, CAI Jing

CollegeofGeo-explorationSciencesandTechnology,JilinUniversity,Changchun130026,China

AbstractAirborne electromagnetic (AEM) exploration is an effective geophysical tool, especially suitable for survey in ridged mountain areas. Theoretical research shows that topography has serious effect on airborne EM responses. Negligence of this effect may create large errors in AEM interpretations. However, until now little study has been done on 3D topographic effect on frequency-/time-domain(FD/TD) airborne EM systems. Most scholars dealt with this problem using a structured grid that can′t simulate topography accurately. We present analgorithm using edge-based unstructured finite-element method(FEM).

KeywordsAirborne EM; 3D forward modeling; Topographic effect; Unstructured vector finite-element method

基金項(xiàng)目國家自然科學(xué)基金重點(diǎn)和面上項(xiàng)目(41530320，41274121)，國家青年基金項(xiàng)目(41404093)，國家重大科研裝備研究項(xiàng)目(ZDYZ2012-1-03，20130523MTEM05)聯(lián)合資助.

作者簡介張博，1989年生，男，博士生，主要從事航空電磁正反演方面的研究. E-mail:em_zhangbo@163.com *通訊作者殷長春， 1965年生，男，教授，主要從事電磁勘探理論，特別是航空和海洋電磁方面的研究. E-mail:yinchangchun@jlu.edu.cn

doi:10.6038/cjg20160431 中圖分類號P631

收稿日期2015-09-11，2015-11-06收修定稿

張博, 殷長春, 劉云鶴等. 2016. 起伏地表頻域/時(shí)域航空電磁系統(tǒng)三維正演模擬研究.地球物理學(xué)報(bào)，59(4):1506-1520,doi:10.6038/cjg20160431.

Zhang B, Yin C C, Liu Y H, et al. 2016. 3D modeling on topographic effect for frequency-/time-domain airborne EM systems.ChineseJ.Geophys. (in Chinese),59(4):1506-1520,doi:10.6038/cjg20160431.