瞬變電磁法數(shù)值模擬實(shí)驗(yàn)教學(xué)軟件研究

昌彥君, 張 瑩, 肖明順

(1. 中國地質(zhì)大學(xué)(武漢) 地球物理與空間信息學(xué)院, 湖北 武漢 430074；2. 中國地質(zhì)大學(xué)(武漢) 地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室, 湖北 武漢 430074； 3. 中國冶金地質(zhì)調(diào)查局 中南地勘院, 湖北 武漢 430074)

瞬變電磁法數(shù)值模擬實(shí)驗(yàn)教學(xué)軟件研究

昌彥君1,2, 張 瑩1, 肖明順3

(1. 中國地質(zhì)大學(xué)(武漢) 地球物理與空間信息學(xué)院, 湖北 武漢 430074；2. 中國地質(zhì)大學(xué)(武漢) 地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室, 湖北 武漢 430074； 3. 中國冶金地質(zhì)調(diào)查局 中南地勘院, 湖北 武漢 430074)

瞬變電磁法是一種廣泛應(yīng)用于礦產(chǎn)資源、水資源、工程與環(huán)境地質(zhì)調(diào)查的地球物理方法。由于電磁場偏微分方程的求解十分復(fù)雜,使得其理論教學(xué)過程非常費(fèi)時。該文以提高教學(xué)效果為目的開發(fā)出瞬變電磁法數(shù)值模擬實(shí)驗(yàn)教學(xué)軟件。該教學(xué)軟件包括網(wǎng)格剖分模塊、有限元計(jì)算模塊、GS系數(shù)循環(huán)模塊,傅氏系數(shù)循環(huán)模塊等,通過按數(shù)值計(jì)算流程圖的演示,使得學(xué)生更容易理解。數(shù)值模擬實(shí)驗(yàn)教學(xué)軟件簡化了教學(xué)過程,節(jié)約了教學(xué)時間,具有生動形象和易于理解的教學(xué)效果。

瞬變電磁法； 數(shù)值模擬實(shí)驗(yàn)； 2.5維正演； 教學(xué)軟件

瞬變電磁法(transient electromagnetic method,TEM)是一種廣泛應(yīng)用于礦產(chǎn)資源、水資源、工程與環(huán)境地質(zhì)調(diào)查的地球物理方法[1-2]。由于TEM的應(yīng)用廣泛,目前許多大學(xué)都開設(shè)了關(guān)于TEM的基礎(chǔ)理論和應(yīng)用技術(shù)課程。在瞬變電磁法的理論教學(xué)中,電磁場偏微分方程的求解十分復(fù)雜[3],尤其是復(fù)雜模型條件下的數(shù)值模擬,教學(xué)過程非常費(fèi)時[4]。實(shí)驗(yàn)教學(xué)軟件不僅可以簡化教學(xué)過程,而且具有生動形象和易于理解的教學(xué)效果[5],因此,開發(fā)一套瞬變電磁法數(shù)值模擬實(shí)驗(yàn)教學(xué)軟件,對于提高瞬變電磁法課程的教學(xué)效果具有重要的意義。筆者基于多年科研成果[6-8],以實(shí)驗(yàn)教學(xué)為目的,開發(fā)出瞬變電磁法數(shù)值模擬實(shí)驗(yàn)教學(xué)軟件,通過數(shù)值計(jì)算流程圖的演示,使學(xué)生更容易理解復(fù)雜模型情況下的數(shù)值求解過程。

1 瞬變電磁法數(shù)理方程及其求解

瞬變電磁法的數(shù)值模擬內(nèi)容包括一維、二維、2.5維(二維地電模型和三維場源)和三維問題。由于一維問題相對簡單,二維問題不符合實(shí)際的三維源性質(zhì),三維問題涉及的計(jì)算量過大,因而在科研中常選擇介于二維和三維之間的所謂2.5維問題進(jìn)行正演理論的研究[9-11]。相對三維,2.5維的計(jì)算量顯著減少；相對二維,2.5維能更好地接近野外實(shí)際地質(zhì)情況。因此,本文探討的瞬變電磁法數(shù)值模擬實(shí)驗(yàn)教學(xué)軟件也是基于2.5維模型開發(fā)的。

在三維各向同性均勻?qū)щ娊橘|(zhì)中，電(磁)偶極子形成的電磁場滿足Maxwell方程。在時間域中，有如下方程組：

在介質(zhì)分界面上，滿足銜接條件：

其中,ε為介電常數(shù)，σ為電導(dǎo)率，μ為磁導(dǎo)率, Γ為介質(zhì)分界面，n為分界面的法向方向矢量。E為電場矢量，H為磁場矢量；Je為場源電流密度矢量，Jm為場源磁流密度矢量。以下標(biāo)1和2表示不同的介質(zhì)。

電(磁)偶極子源位于笛卡爾坐標(biāo)原點(diǎn)時，有

其中，δ(x)、δ(y)和δ(z)為狄拉克(Dirac)源函數(shù)；Pe、Pm分別為源的電偶極矩和磁偶極矩，且有Pe=IΔs，Δs為電流元長度，Pm=IS，S為小電流環(huán)的面積；μ0是真空中的磁導(dǎo)率；I為按階躍規(guī)律變化的供電電流強(qiáng)度，它滿足下式:

瞬變電磁法三維正演問題為時間(t)和空間(x,y,z)域的四維復(fù)雜問題。所謂的2.5維問題，是三維問題的一個特例，即二維地電模型(模型的電、磁性參數(shù)在地質(zhì)體走向方向無變換)和三維場源。求解時，利用拉普拉斯變換對時間(t)消維，將時空(t,x,y,z)域四維問題轉(zhuǎn)換為離散的拉普拉斯域的(x,y,z)三維問題；對走向(y)方向電磁場分量作傅里葉變換，進(jìn)一步降維得到離散的拉普拉斯和傅里葉域下的(x,z)二維問題；運(yùn)用有限元單元法求解二維問題，通過反傅里葉變換和逆拉普拉斯變換得到時空(t,x,y,z)域電磁場響應(yīng)，這便是時間域有源電磁法2.5維有限元數(shù)值模擬的計(jì)算思路[11-12]。

2 算法流程及其功能模塊

瞬變電磁法數(shù)值模擬算法流程圖對學(xué)生理解瞬變電磁法的數(shù)值計(jì)算過程有很大幫助，該流程圖反映了瞬變電磁法的核心計(jì)算原理和計(jì)算步驟，是整個數(shù)值模擬的精華和學(xué)習(xí)要點(diǎn)。圖1為軟件的計(jì)算流程圖，各主要功能模塊介紹如下。

圖1 瞬變電磁法數(shù)值模擬算法流程圖

(1) 網(wǎng)格剖分預(yù)處理模塊。網(wǎng)格剖分預(yù)處理模塊負(fù)責(zé)軟件的主要GUI交互。運(yùn)用有限單元算法求解離散的拉普拉斯和傅里葉域下的(x,z)二維問題，需要剖分地下計(jì)算區(qū)域網(wǎng)格，建立網(wǎng)格單元電性參數(shù)和幾何參數(shù)信息，提供給用于數(shù)值計(jì)算的FORTRAN動態(tài)計(jì)算庫。軟件使用者僅僅需要提供測點(diǎn)信息和采樣時間參數(shù)信息等，程序?qū)⒆詣悠史钟邢拊W(wǎng)格，并提供GUI交互供軟件使用者建立地電網(wǎng)格模型，劃分異常地質(zhì)體。軟件并未解決2.5維瞬變電磁的全自動反演，僅僅提供了一種粗糙的人機(jī)聯(lián)作試錯反演。因此，正演和試錯反演的網(wǎng)格剖分模塊是一致的。

(2) 數(shù)據(jù)參數(shù)傳遞模塊。數(shù)據(jù)參數(shù)傳遞模塊負(fù)責(zé)軟件向FORTRAN動態(tài)計(jì)算庫的數(shù)據(jù)傳輸，即傳輸?shù)叵掠?jì)算區(qū)域網(wǎng)格單元的電性和幾何參數(shù)以及其他信息。由于瞬變電磁法數(shù)值模擬需要考慮空氣因素的影響，實(shí)際計(jì)算的有限元網(wǎng)格需要包括向上拓展的空氣網(wǎng)格，因此數(shù)據(jù)參數(shù)傳遞模塊也負(fù)責(zé)擴(kuò)充地下網(wǎng)格，形成適用于計(jì)算的計(jì)算網(wǎng)格。

(3) 有限元計(jì)算模塊。有限元計(jì)算模塊是FORTRAN動態(tài)計(jì)算庫的核心模塊，即計(jì)算拉普拉斯和傅里葉域下不同測點(diǎn)位置的(x,z)二維電磁場響應(yīng)。FORTRAN動態(tài)計(jì)算庫采用了Gaver-Stehfest逆拉普拉斯算法進(jìn)行拉普拉斯逆變換[13]，一個時間點(diǎn)僅僅需要十幾個拉普拉斯域的計(jì)算值。因此，2.5維瞬變電磁場響應(yīng)的計(jì)算為時間、傅里葉系數(shù)、Gaver-Stehfest系數(shù)、場源測點(diǎn)位置的四重循環(huán)。

(4) 源點(diǎn)循環(huán)模塊。源點(diǎn)循環(huán)是時間、傅里葉系數(shù)、Gaver-Stehfest系數(shù)、場源測點(diǎn)位置四重循環(huán)的最里層循環(huán)。上述有限元計(jì)算到最后歸結(jié)于一個線性方程組的求解。在相同的地電模型參數(shù)、相同的采樣時間、發(fā)射線圈參數(shù)和一定的近似條件下,不同測點(diǎn)的有限元線性方程組有相同的計(jì)算部分。因此,將源點(diǎn)循環(huán)放在四重循環(huán)最里層,可減少重復(fù)計(jì)算量。

(5) GS系數(shù)循環(huán)模塊。GS系數(shù)循環(huán)即Gaver-Stehfest系數(shù)循環(huán)。由于Gaver-Stehfes逆拉普拉斯變換算法對精度要求很高[14],因此,將其放在第三層循環(huán)。

(6) 傅氏系數(shù)循環(huán)模塊。傅氏系數(shù)循環(huán)即傅里葉系數(shù)循環(huán)。通過Gaver-Stehfes逆變換得到不同的采樣時間、不同的測點(diǎn)位置、不同的傅里葉系數(shù)下的二維電磁場響應(yīng)之后,運(yùn)用數(shù)值計(jì)算實(shí)現(xiàn)傅里葉反變換,得到不同的采樣時間、不同的測點(diǎn)位置下的空間域電磁場響應(yīng)值。

(7) 數(shù)據(jù)傳遞輸出模塊。數(shù)據(jù)傳遞輸出模塊負(fù)責(zé)FORTRAN動態(tài)計(jì)算庫的數(shù)據(jù)傳輸,即傳輸瞬變電磁場響應(yīng)值,并調(diào)用相關(guān)計(jì)算模塊得到視電阻率、視深度等輔助信息。該模塊亦負(fù)責(zé)數(shù)據(jù)的圖形顯示,包括視電阻率擬斷面圖和磁感電動勢多測道剖面圖等。

3 實(shí)驗(yàn)教學(xué)軟件的部分演示

軟件的開發(fā)是基于VC6.0平臺進(jìn)行的,采取界面GUI開發(fā)和數(shù)值計(jì)算核心剝離開發(fā)的思路,以便于增加新的計(jì)算模型,以及根據(jù)教學(xué)需要不斷地完善與擴(kuò)充內(nèi)容。瞬變電磁法數(shù)值模擬實(shí)驗(yàn)教學(xué)軟件的部分功能演示如下:

(1) 導(dǎo)入測點(diǎn)數(shù)據(jù)；(2)選擇測點(diǎn)數(shù)據(jù)(見圖2)；(3)形成網(wǎng)格(見圖3)；(4)選擇矩形工具并用鼠標(biāo)左鍵畫矩形；(5)右鍵點(diǎn)擊賦予電性參數(shù)(見圖4)；(6)正演計(jì)算(見圖5)；(7)查看視電阻率斷面圖(見圖6)。

圖3 形成網(wǎng)格界面

圖4 賦予電性參數(shù)界面圖

圖5 正演計(jì)算頁面

圖6 視電阻率斷面圖

4 結(jié)束語

通過數(shù)值模型和計(jì)算流程圖,將復(fù)雜的電磁場偏微分方程的求解過程形象化,使得瞬變電磁法課程的教學(xué)過程得到了簡化。軟件中儲存的大量常用模型的數(shù)值計(jì)算結(jié)果,極大地豐富了教學(xué)內(nèi)容,縮短了授課時間。對于復(fù)雜且不容易理解的教學(xué)內(nèi)容,通過使用數(shù)值模擬實(shí)驗(yàn)教學(xué)軟件,不但提高了教學(xué)效果,也激發(fā)了學(xué)生的學(xué)習(xí)興趣,促使學(xué)生從被動學(xué)習(xí)轉(zhuǎn)變?yōu)橹鲃訉W(xué)習(xí),這種教學(xué)方法值得進(jìn)一步完善和推廣。

References)

[1] 陳載林,黃臨平,陳玉梁.我國瞬變電磁法應(yīng)用綜述[J].鈾礦地質(zhì),2010,26(1):51-54.

[2] 郭寧寧,昌彥君.瞬變電磁法在城區(qū)建筑基坑中的防空洞探測[J].工程地球物理學(xué)報(bào),2011,8(6):705-708.

[3] 蔣邦遠(yuǎn).實(shí)用近區(qū)磁源瞬變電磁法勘探[M].北京:地質(zhì)出版社,1998.

[4] Wannamaker P E,Hohmann G W,Sanfilipo W A.Electromagnetic modeling of three-dimensional bodies in layered earths using integral equations[J].Geophysics,1984,49(1):60-74.

[5] 昌彥君,張瑩,曹中林.探地雷達(dá)數(shù)值模擬實(shí)驗(yàn)研究[J].實(shí)驗(yàn)技術(shù)與管理,2009,26(4):69-72.

[6] 昌彥君,張桂青.電磁場從頻率域轉(zhuǎn)換到時間域的幾種算法比較[J].物探化探計(jì)算技術(shù),1995(3):25-29.

[7] 昌彥君,肖明順,武毅.瞬變電磁數(shù)據(jù)一維反演時的初值選取研究[J].石油地球物理勘探,2010,45(2),295-298.

[8] 昌彥君,肖明順,孟永良,等.帶地形的瞬變電磁2.5維有限元數(shù)值模擬研究[G]//地球物理科學(xué)研究十年進(jìn)展論文集.武漢:中國地質(zhì)大學(xué)出版社,2012:267-278.

[9] Mitsuhata Y.2-D electromagnetic modeling by finite-element method with a dipole source and topography[J].Geophysics,2000,65(2):465-475.

[10] 孟永良,羅延鐘,昌彥君.時間譜電阻率法的二維正演算法[J].地球科學(xué),2000,25(6):656-662.

[11] 王華軍,羅延鐘.中心回線瞬變電磁2.5維維有限元算法[J].地球物理學(xué)報(bào),2003,46(6):855-862.

[12] 徐世浙.地球物理中的有限單元法[M].北京:科學(xué)出版社,1994.

[13] Knight J H,Raiche A P.Transient electromagnetic calculation using the Gaver-Stehfest inverse Laplace transform method[J].Geophysics,1982,47(1):47-50.

[14] 羅延鐘,昌彥君.G-S變換快速算法[J].地球物理學(xué)報(bào),2000,43(5):684-690.

Research on teaching software of numerical simulation of transient electromagnetic method

Chang Yanjun1,2, Zhang Ying1, Xiao Mingshun3

(1. Institute of Geophysics and Geomatics,China University of Geosciences,Wuhan 430074,China; 2. Subsurface Multi-scale Imaging Laboratory,China University of Geosciences,Wuhan 430074,China; 3. Zhongnan Institute of Geological Exploration,China Metalurgical Geology Bureau,Wuhan 430074,China)

The transient electromagnetic method (TEM) is a kind of tool used widely in mineral and oil exploration as well as water exploration,engineering and environment geophysics exploration.Its theory teaching process is very time-consuming due to the complexity of the solution of electromagnetic field partial differential equations.This paper develops the numerical simulation experiment teaching software of the transient electromagnetic method for the purpose of experiment teaching based on scientific research.The teaching software modules,including grid subdivision finite element calculation module,GS coefficient circulation module,Fourier coefficient circulation module and so on,make students easy to understand by demonstrating the flow chart of numerical calculation. It may not only simplify the teaching process and save the teaching time, but also have good teaching effect of vivid image and easy to understand.

transient electromagnetic method; numerical simulation experiment; 2.5-dimensional forward; teaching software

2014- 11- 25 修改日期：2015- 01- 07

中國地質(zhì)大學(xué)創(chuàng)新實(shí)驗(yàn)和地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室基金項(xiàng)目資助

昌彥君(1964—),男,湖北仙桃,博士,教授,主要從事電磁法勘探的數(shù)值模擬研究.

P631

A

1002-4956(2015)7- 0131- 04