基于Julia語言的大地電磁三維正演模擬

邢雯淋, 曹 輝, 曹禮剛, 李小翠, 魏統(tǒng)彪, 張朝暉

(成都理工大學(xué) 地球勘探與信息技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，成都 610059)

0 引言

大地電磁測深(Magnetotelluric Sounding，MT)是由20世紀(jì)50年代蘇聯(lián)學(xué)者Tikhonov[1]和法國學(xué)者L. Cagnird[2]提出的以天然的平面電磁波為場源，根據(jù)趨膚效應(yīng)的原理研究地球電性結(jié)構(gòu)的一種地球物理勘探方法。目前已被應(yīng)用于各種領(lǐng)域，在油氣勘探、構(gòu)造研究、礦產(chǎn)探測、地震預(yù)報(bào)、地?zé)峥碧降确较蚓邢鄳?yīng)的研究[3-5]。

隨著目前計(jì)算設(shè)備和運(yùn)算軟件的高速發(fā)展以及數(shù)值計(jì)算方法的進(jìn)步，三維正反演發(fā)展取得了巨大的進(jìn)展[6]。但目前大部分關(guān)于大地電磁三維正反演的文獻(xiàn)主要集中在數(shù)值方法的闡述和討論上，而很少關(guān)注如何有效地對數(shù)值方法進(jìn)行編程，僅有幾種基于不同數(shù)值方法并充分采用了現(xiàn)代計(jì)算技術(shù)的三維編碼，被學(xué)術(shù)界廣泛應(yīng)用于現(xiàn)場實(shí)測數(shù)據(jù)，如Siripunvaraporn等[7]采用數(shù)據(jù)空間的方法使矩陣的維數(shù)取決于數(shù)據(jù)集的大小，而不是模型參數(shù)的數(shù)量；Kelbert等[8]設(shè)計(jì)了基于電磁反演計(jì)算機(jī)代碼的模塊化系統(tǒng)(稱為ModEM)，這些代碼在效率、精度和處理多尺度結(jié)構(gòu)和復(fù)雜地形的能力方面表現(xiàn)出了出色的性能;Geraskin[9]提出，通常三維電磁代碼中只有一小部分專門用于密集計(jì)算(一般采用Fortran和C，即靜態(tài)語言)，大部分的瑣碎工作使用Python和MATLAB(即動(dòng)態(tài)語言)更加便捷。這說明一般的三維電磁代碼的編寫主要還是使用 “雙重編程語言”來完成，而一個(gè)嶄新語言的Julia正好可以簡化“雙重編程語言”帶來的復(fù)雜性，它是在一定的性能需求下，犧牲一些不那么重要的動(dòng)態(tài)性尋找一個(gè)更優(yōu)的平衡點(diǎn)，可以僅使用一種語言來進(jìn)行編程開發(fā)，因此我們可以大大減少編碼期間所需的工作，并使代碼更容易理解和維護(hù)[10]。

因此，筆者嘗試使用一個(gè)全新的語言，Julia來實(shí)現(xiàn)由有限體積法(MFV)建立的基于各向同性介質(zhì)的大地電磁測深法三維正演模型，并對其進(jìn)行數(shù)值模擬。

1 Julia語言簡介

由于程序員在編寫高級而復(fù)雜的算法邏輯時(shí)習(xí)慣使用簡單的動(dòng)態(tài)語言(如Python、MATLAB)，但為了追求計(jì)算速度，編寫完后經(jīng)常會(huì)將需要密集計(jì)算的部分轉(zhuǎn)換為靜態(tài)語言(如C、Fortran)。而編寫兩種語言的代碼需要去考慮兩種語言之間的類型轉(zhuǎn)換和內(nèi)存管理，其復(fù)雜性可能更甚于僅使用其中一種語言，為此一種新的高級編程語言Julia被開發(fā)出來了[11]。

Julia是一種高級的、高性能的、動(dòng)態(tài)的語言，在科學(xué)計(jì)算領(lǐng)域出類拔萃，其語法像MATLAB或Python一樣簡單易懂，并且與傳統(tǒng)的動(dòng)態(tài)語言相比Julia的核心語言很小，有著豐富的基礎(chǔ)類型，具有可選類型標(biāo)注和多重派發(fā)這兩個(gè)特性，還有接近C語言的性能。Bezanson等[11]將Julia語言與其他六種語言(C++、Python、MATLAB、Octave、R 和 JavaScript)的速度進(jìn)行了定量分析比較(表1)，結(jié)果表明，Julia語言在各項(xiàng)測試中的速度均較快，因?yàn)镴ulia語言在設(shè)計(jì)之初就考慮如何利用現(xiàn)有的技術(shù)去高效加速動(dòng)態(tài)語言，它有一個(gè)與其系統(tǒng)結(jié)合在一起的即時(shí)編譯器(JIT)，提供了比Python，MATLAB等更高效地交互式編程和動(dòng)態(tài)性的同時(shí)，也可以達(dá)到靜態(tài)編譯語言(如C和Fortran)一般地性能，相當(dāng)于在一定程度上對“雙重編程語言”問題進(jìn)行了優(yōu)化[10]。

表1 微基準(zhǔn)測試結(jié)果(相對于C++的時(shí)間)[11]Tab.1 Microbenchmark results (times relative to C++)[11]

正是Julia語言的這些特性和它在數(shù)值計(jì)算上的專業(yè)性，才使其能有效降低編程工作的復(fù)雜性，并且擁有寬泛的應(yīng)用范圍，能被筆者運(yùn)用于實(shí)現(xiàn)大地電磁三維正演模擬。

2 大地電磁三維正演算法

2.1 邊值問題

大地電磁法遵循麥克斯韋方程，取時(shí)間因子為e-iωt，可滿足如下方程[3]：

▽×E=iωμH

(1)

▽×H=σ-iωεE

(2)

其中：E為電場；H為磁場；ω=2πf為角頻率；f為頻率；ε為介電常數(shù)；μ為真空中的磁導(dǎo)率；σ是介質(zhì)的電導(dǎo)率。

將式 (1) 代入式 (2) 可導(dǎo)出電場E所應(yīng)滿足的微分方程為式(3)。

▽×(▽×E)-k2E=0

(3)

式中，k2=iωμσ+ω2εμ≈iωμσ，由于在大地電磁勘探低頻時(shí)位移電流可忽略不計(jì)。

為了確保控制方程有確定解，必須設(shè)定所謂的“自然”邊界條件，由于大地電磁場的源都來自地球的外部，可以認(rèn)為初始大地電磁場是平面波場，初始電場E的偏振方向沿X軸(圖1)。由于空氣中的電磁場受到地形和不均勻體的影響而不均勻的分布，因此在數(shù)值模擬中必須包括空氣層和地下介質(zhì)[6]。如圖1所示，模擬區(qū)域設(shè)為Ω，由空氣層和地下介質(zhì)兩部分組成，Ω=Ω1∪Ω2。在模擬區(qū)域中，設(shè)空氣層的電導(dǎo)率為10-5S/m，這時(shí)空氣和地下介質(zhì)的接觸面(地表)變成內(nèi)部邊界。

圖1 正演三維模型示意圖Fig.1 Schematic diagram of forward 3D model

取狄利克雷邊界條件：

E(x,y,z)|?Ω=g(x,y,z)|?Ω

(4)

式中：g為邊界上的矢量電場;Ω為模擬區(qū)域，可以使用一維或二維的大地電磁計(jì)算值[12]。聯(lián)立式(3)和式(4)建立了大地電磁三維正演算法的邊值問題。

2.2 矢量有限體積分析

這里使用有限體積法求解式(1)和式(2)的解，為了處理任意地形和復(fù)雜電導(dǎo)率結(jié)構(gòu)模型，我們將電場和磁場在Yee[13-14]提出的交錯(cuò)網(wǎng)格上進(jìn)行離散化 (圖2)。

圖2 Yee交錯(cuò)網(wǎng)格用于確定電場和磁場位置Fig.2 Yee staggered grid is used to determine the position of electric and magnetic field

參照Haber[15]使用有限體積法對式(1)和式(2)離散化完整詳細(xì)的處理，可得到如下離散的麥克斯韋方程組的類比：

CURLe+iωb=0

(5)

CURLTMf(μ)b-Me(σ)e=s

(6)

A(σ)e=(CURLTMf(μ)CURL+

iωMe(σ)e=-iωs=q

(7)

式中：Me(σ)是唯一包含電導(dǎo)率的項(xiàng)。Haber[16]給出了各向同性電導(dǎo)率單元的Me(σ)的顯式表達(dá)式為式(8)。

(8)

式中：矩陣Aej為在j方向上單元中心到邊緣的平均變量；v是單元體積的矢量；⊙是點(diǎn)方向的哈達(dá)瑪積。

2.3 大型線性方程組的求解

在大地電磁的數(shù)值模擬中，最終都會(huì)得出一個(gè)大型稀疏線性方程組，式(7)最終可以表示為式(9)[16]。

A(σ)e=K+iωMe(σ)e=q

(9)

式中，K是由▽×μ-1▽×算子離散化得到的對稱半正定矩陣。筆者使用直接求解法通過直接分解系數(shù)矩陣來求解該稀疏線性方程組，其得出的結(jié)果更具高精度和高穩(wěn)定性[17]。

近年來，在計(jì)算成本方面已經(jīng)對大型稀疏線性方程組的直接求解技術(shù)進(jìn)行了高度優(yōu)化，已有開源的MPI分布式并行求解器MUMPS程序包和OpenMP共享內(nèi)存式并行求解器Pardiso程序包等[18-19]。隨著計(jì)算機(jī)性能的不斷升級，目前已經(jīng)可以使用直接法在個(gè)人計(jì)算機(jī)上建立適度大小的模型進(jìn)行大地電磁的三維正演模擬[20]。綜上所述，筆者選擇采用MUMPS求解器，對大地電磁三維正演中大型線性方程組進(jìn)行求解。

3 Julia語言與大地電磁的結(jié)合

近年來，在三維地球物理電磁模擬和實(shí)用軟件的開發(fā)中Julia得到了一定的運(yùn)用，如Peng等[21]完成了基于非精確高斯牛頓法的海洋可控源電磁三維反演其在Julia中的數(shù)值實(shí)現(xiàn)；B?rner等[22]使用Julia語言完成了瞬變電磁數(shù)據(jù)的三維反演；Ruthotto等[23]為了解決從噪聲和間接測量中估計(jì)偏微分方程(PDEs)的參數(shù)而求解不適定逆問題所需的大量計(jì)算，他們提出了一個(gè)用Julia語言編寫的開源軟件jInv，用于解決具有多種測量值的參數(shù)估計(jì)問題的并行算法；Han 等[20]用Julia語言實(shí)現(xiàn)采用插值方法將電場的不同分量平均到同一位置的一般各向異性三維大地電磁正演的有限體積算法，編寫過程非常高效，并產(chǎn)生了具有良好可讀性、可維護(hù)性和可擴(kuò)展性的代碼。由此可知,Julia語言可高效便捷地應(yīng)用于三維地球物理電磁模擬。

Julia提供了用于所有基本運(yùn)算的簡單函數(shù)，并對模塊和復(fù)合類型有著良好的支持，這里參照Han 等[20]所使用的代碼模塊化用于三維大地電磁各向同性建模，其結(jié)構(gòu)如圖3所示。其中，線性求解模塊使用了第三方的Julia包 MUMPS.jl[24]，提供了MUMPS求解器。

圖3 三維大地電磁各向同性正演建模代碼的主要結(jié)構(gòu)[20]Fig.3 The main structure of the 3D magnetotelluric isotropic forward modeling code [20]

4 Julia語言正演準(zhǔn)確性驗(yàn)證

為檢驗(yàn)本文采用Julia語言編寫的大地電磁三維正演算法程序模擬結(jié)果的準(zhǔn)確性，選取由MTNet國際論壇上發(fā)布的全球公共模型——Dublin測試模型1 (DTM1)進(jìn)行正確性驗(yàn)證[25]。

DTM1是一個(gè)三維正演模型，圖4是在均勻的半空間中由3個(gè)各向同性的異常體組合而成，均勻半空間的電阻率為100 Ω·m，各個(gè)異常體的電阻率如表2所示。

表2 DTM1模型異常體相關(guān)數(shù)據(jù)Tab.2 Data related to abnormal bodies in DTM1 model

圖4 DTM1模型三維示意圖Fig.4 3D schematic diagram of DTM1 model(a)xz軸垂直橫截面；(b)yz軸垂直橫截面；(c)xy軸水平橫截面

將模型域離散為Nx×Ny×Nz=50×55×82網(wǎng)格單元(包括7個(gè)空氣層)，最小單元尺寸為2×1×0.5 km，選用0. 1和10 000 s之間21個(gè)對數(shù)間隔周期進(jìn)行計(jì)算，并對上述頻率在(0，0，0)觀測點(diǎn)的正演結(jié)果與其他研究者模擬的結(jié)果進(jìn)行分析比對(圖5)。

圖5 本文Julia代碼正演結(jié)果(My)與其他研究者正演結(jié)果在(0,0,0) 觀測點(diǎn)視電阻率和相位的對比Fig.5 Comparison of apparent resistivity and phase at (0,0,0) observation point between forward modeling results of Julia code (My) in this paper and those of other researchers(a)視電阻率_XY模式；(b)相位_XY模式；(c)視電阻率_YX模式;(d)相位_YX模式(e)視電阻率_XX模式；(f)相位_XX模式；(g)視電阻率_YY模式;(h)相位_YY模式

由圖5可知，本文方法計(jì)算出視電阻率和相位的曲線與其他研究者得出的結(jié)果基本相吻合，只有在XX和YY模式時(shí)的相位高頻段上擬合較差(由于視電阻率數(shù)值太小，相位容易受到精度影響)，其他模式的視電阻率和相位無論在低頻還是在高頻上擬合的效果都很好。

對基于 Julia語言編寫的大地電磁正演程序的結(jié)果進(jìn)行精度分析，利用視電阻率和相位誤差計(jì)算公式：

(10)

(11)

式中：ερ為視電阻率相對誤差值；ρJulia為Julia程序模擬視電阻率值；ρOther為其他研究者的模擬視電阻率值；n為頻點(diǎn)的個(gè)數(shù)；εφ為相位相對誤差值；φJ(rèn)ulia為Julia程序模擬相位值；φOther為其他研究者的模擬視相位值；n為頻點(diǎn)的個(gè)數(shù)。表3為本文Julia程序與其他研究者XY和YX視電阻率和相位的模擬結(jié)果的相對誤差值。

表3 與其他研究者正演結(jié)果在(0,0,0) 觀測點(diǎn)視電阻率和相位的相對誤差Tab.3 Relative error of apparent resistivity and phase at (0,0,0) observation point between forward modeling results and those of other researchers

由表3可以看出，除了與wsinv3dmt的模擬結(jié)果相對誤差較大，與其他結(jié)果的視電阻率相對誤差值均小于7%，相位相對誤差值均小于0.02°，符合電法勘探的精度要求。

5 模型試算

為了驗(yàn)證基于 Julia語言編寫的大地電磁正演程序的可行性，圖6設(shè)計(jì)了一個(gè)高低阻組合模型，其中模型1為2 km×2 km×1.7 km，電阻率為10 Ω·m，埋深為0.3 km；模型2為2 km×1 km×1.7 km，電阻率為1 Ω·m，埋深為0.3 km；模型3為4 km×2 km×1 km，電阻率為1 000 Ω·m，埋深為2 km；背景場電阻率為100 Ω·m，空氣層電阻率設(shè)為10-8Ω·m。并根據(jù)圖6模型設(shè)計(jì)了圖7的兩種模型與之對比。

圖6 高低阻組合模型三維示意圖Fig.6 3D schematic diagram of combination model of high and low resistance(a)正視圖；(b)俯視圖

圖7 對比模型Fig.7 Comparison model(a)低阻模型;(b)低阻-高阻模型

先將模型離散為Nx×Ny×Nz= 64×40×47網(wǎng)格單元(包括7個(gè)空氣層)，然后在地表布置測線范圍為X=[-4 000 m，4 000 m] 和Y=[-4 000 m，4 000 m]，測點(diǎn)和測線之間的距離均為 200 m，最后對 0.1 Hz、1 Hz、10 Hz 和 100 Hz 的不同頻率的正演結(jié)果進(jìn)行計(jì)算，選取結(jié)果見圖8和圖9。對比模型0.1 Hz的結(jié)果如圖10所示。

圖8 XY和YX模式在Y=0測線上各頻率的對比圖Fig.8 Comparison of each frequencies of XY and YX modes on the survey line Y=0(a)視電阻率_XY模式；(b)相位_XY模式;(c)視電阻率_YX模式;(d)相位_YX模式

圖9 XY和YX模式的視電阻率和相位的等值面圖Fig.9 Isosurface plots of apparent resistivity and phase for XY and YX modes(a)0.1 Hz_XY_Resistivity/Ω·m;(b)0.1 Hz_XY_Phase/°;(c)0.1 Hz_XY_Resistivity/Ω·m;(d)0.1 Hz_XY_Phase/°;(e)1 Hz_XY_Resistivity/Ω·m;(f)1 Hz_XY_Phase/°;(g)1 Hz_XY_Resistivity/Ω·m;(h)1 Hz_XY_Phase/°;(i)10 Hz_XY_Resistivity/Ω·m;(j)10 Hz_XY_Phase/°;(k)10 Hz_XY_Resistivity/Ω·m;(l)10 Hz_XY_Phase/°(m)100 Hz_XY_Resistivity/Ω·m;(n)100 Hz_XY_Phase/°;(o)100 Hz_XY_Resistivity/Ω·m;(p)100 Hz_XY_Phase/°

圖10 對比模型0.1 Hz時(shí)XY 和YX模式的視電阻率和相位的等值面圖Fig.10 Isosurface of apparent resistivity and phase of XY and YX modes at 0.1 Hz for the comparison model(a)低阻模型;(b)低阻-高阻模型

從圖8與圖9可以看出，0.1 Hz時(shí)XY模式和YX模式均對模型中淺部低阻異常體在水平面上的數(shù)值大小和位置顯示較準(zhǔn)確(其中YX模式的位置顯示更為準(zhǔn)確)，并且兩個(gè)低阻體之間較小的電阻差異也能體現(xiàn)出來，且視電阻率值也與異常體實(shí)際電阻率較為接近；而對深部高阻體時(shí)，XY模式可以較好地顯示出水平面上淺部低阻體和深部高阻體重疊的綜合反映，而YX模式的顯示容易對高阻體的位置產(chǎn)生誤判。1 Hz時(shí)XY模式和YX模式的顯示與0.1 Hz較相近，均對模型淺部低阻體數(shù)值大小和位置顯示較好，XY模式視電阻率對深部高阻體的數(shù)值大小顯示較0.1 Hz略差，但也能明顯看出有高阻體的存在。

10 Hz時(shí)XY和YX模式對模型淺部低阻體數(shù)值大小和位置均有適度的顯示，但除了XY模式相位圖對深部高阻體略有顯現(xiàn)，其余視電阻率和相位圖均已完全無法顯示出深部高阻體的存在；100 Hz時(shí)XY和YX模式不僅對模型淺部低阻體數(shù)值大小顯示較差，并且完全也未能顯示出深部高阻體的存在，但對淺部低阻體的位置略有顯示。

對比模型在0.1 Hz時(shí)XY與YX結(jié)果可知，低阻模型與低阻-高阻模型，均對模型中淺部低阻異常體在水平面上的數(shù)值大小和位置顯示較準(zhǔn)確，其中低阻-高阻模型的結(jié)果，在深部高阻體處只有XY模式與低阻模型在數(shù)值上有所差異，與高低阻組合模型的結(jié)果情況基本一致。

綜上所述，對本文設(shè)計(jì)的高低阻組合模型進(jìn)行正演模擬，不同頻率下的正演結(jié)果對于淺部低阻異常體的位置反映均較好，但低頻時(shí)對淺部低阻異常體的模擬結(jié)果均較高頻時(shí)準(zhǔn)確，表明頻率越低異常體的分辨率越高，且由于低頻的勘探深度大，只有在低頻時(shí)的正演模擬結(jié)果才能顯示出對深部高阻體的響應(yīng)。

5 總結(jié)

這里通過Julia語言編寫的正演代碼實(shí)現(xiàn)由有限體積法的大地電磁各向同性的三維正演，并進(jìn)行了正確性驗(yàn)證，得到如下結(jié)論：

1)筆者參考全球公共模型DTM1進(jìn)行正確性驗(yàn)證，并對高低阻組合模型進(jìn)行正演試算，結(jié)果表明，在正確性驗(yàn)證中與多個(gè)研究者的正演結(jié)果進(jìn)行對比分析結(jié)果擬合較好；在模型試算中各頻率對淺層低阻體的正演效果較好，而對深部高阻體只有在低頻時(shí)才能得到較好的正演結(jié)果，這與大地電磁的相關(guān)理論相符，并且正演結(jié)果基本能將異常體的形態(tài)和位置特征體現(xiàn)出來。綜上所述兩個(gè)模型的正演結(jié)果,均證明了使用基于Julia語言的大地電磁正演程序是準(zhǔn)確、可靠、可行的。

2)對于通常需要進(jìn)行正反演模擬的研究者和學(xué)生來說，Julia作為一種高性能性和高動(dòng)態(tài)性的新型編程語言，提供了跟Python，MATLAB等相比更加高效地交互式編程和動(dòng)態(tài)性，并也達(dá)到了一般的靜態(tài)編譯語言(C、Fortran)的性能，它語言結(jié)構(gòu)簡單明了容易掌握，并在在數(shù)值計(jì)算上有著極強(qiáng)的專業(yè)性，可以更快更高效更簡易的進(jìn)行我們正反演代碼的編寫工作。