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

邢雯淋, 曹 輝, 曹禮剛, 李小翠, 魏統彪, 張朝暉

(成都理工大學 地球勘探與信息技術教育部重點實驗室，成都 610059)

0 引言

大地電磁測深(Magnetotelluric Sounding，MT)是由20世紀50年代蘇聯學者Tikhonov[1]和法國學者L. Cagnird[2]提出的以天然的平面電磁波為場源，根據趨膚效應的原理研究地球電性結構的一種地球物理勘探方法。目前已被應用于各種領域，在油氣勘探、構造研究、礦產探測、地震預報、地熱勘探等方向均有相應的研究[3-5]。

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

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

1 Julia語言簡介

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

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

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

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

2 大地電磁三維正演算法

2.1 邊值問題

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

▽×E=iωμH

(1)

▽×H=σ-iωεE

(2)

其中：E為電場；H為磁場；ω=2πf為角頻率；f為頻率；ε為介電常數；μ為真空中的磁導率；σ是介質的電導率。

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

▽×(▽×E)-k2E=0

(3)

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

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

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

取狄利克雷邊界條件：

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

(4)

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

2.2 矢量有限體積分析

這里使用有限體積法求解式(1)和式(2)的解，為了處理任意地形和復雜電導率結構模型，我們將電場和磁場在Yee[13-14]提出的交錯網格上進行離散化 (圖2)。

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

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

CURLe+iωb=0

(5)

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

(6)

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

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

(7)

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

(8)

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

2.3 大型線性方程組的求解

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

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

(9)

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

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

3 Julia語言與大地電磁的結合

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

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

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

4 Julia語言正演準確性驗證

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

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

表2 DTM1模型異常體相關數據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網格單元(包括7個空氣層)，最小單元尺寸為2×1×0.5 km，選用0. 1和10 000 s之間21個對數間隔周期進行計算，并對上述頻率在(0，0，0)觀測點的正演結果與其他研究者模擬的結果進行分析比對(圖5)。

圖5 本文Julia代碼正演結果(My)與其他研究者正演結果在(0,0,0) 觀測點視電阻率和相位的對比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可知，本文方法計算出視電阻率和相位的曲線與其他研究者得出的結果基本相吻合，只有在XX和YY模式時的相位高頻段上擬合較差(由于視電阻率數值太小，相位容易受到精度影響)，其他模式的視電阻率和相位無論在低頻還是在高頻上擬合的效果都很好。

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

(10)

(11)

式中：ερ為視電阻率相對誤差值；ρJulia為Julia程序模擬視電阻率值；ρOther為其他研究者的模擬視電阻率值；n為頻點的個數；εφ為相位相對誤差值；φJulia為Julia程序模擬相位值；φOther為其他研究者的模擬視相位值；n為頻點的個數。表3為本文Julia程序與其他研究者XY和YX視電阻率和相位的模擬結果的相對誤差值。

表3 與其他研究者正演結果在(0,0,0) 觀測點視電阻率和相位的相對誤差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的模擬結果相對誤差較大，與其他結果的視電阻率相對誤差值均小于7%，相位相對誤差值均小于0.02°，符合電法勘探的精度要求。

5 模型試算

為了驗證基于 Julia語言編寫的大地電磁正演程序的可行性，圖6設計了一個高低阻組合模型，其中模型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，空氣層電阻率設為10-8Ω·m。并根據圖6模型設計了圖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網格單元(包括7個空氣層)，然后在地表布置測線范圍為X=[-4 000 m，4 000 m] 和Y=[-4 000 m，4 000 m]，測點和測線之間的距離均為 200 m，最后對 0.1 Hz、1 Hz、10 Hz 和 100 Hz 的不同頻率的正演結果進行計算，選取結果見圖8和圖9。對比模型0.1 Hz的結果如圖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時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時XY模式和YX模式均對模型中淺部低阻異常體在水平面上的數值大小和位置顯示較準確(其中YX模式的位置顯示更為準確)，并且兩個低阻體之間較小的電阻差異也能體現出來，且視電阻率值也與異常體實際電阻率較為接近；而對深部高阻體時，XY模式可以較好地顯示出水平面上淺部低阻體和深部高阻體重疊的綜合反映，而YX模式的顯示容易對高阻體的位置產生誤判。1 Hz時XY模式和YX模式的顯示與0.1 Hz較相近，均對模型淺部低阻體數值大小和位置顯示較好，XY模式視電阻率對深部高阻體的數值大小顯示較0.1 Hz略差，但也能明顯看出有高阻體的存在。

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

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

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

5 總結

這里通過Julia語言編寫的正演代碼實現由有限體積法的大地電磁各向同性的三維正演，并進行了正確性驗證，得到如下結論：

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

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