基于球坐標(biāo)系下有限差分的地磁測深三維正演

李建平，翁愛華，李世文，李大俊， 李斯睿，楊 悅，唐 裕，張艷輝

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

0 引言

低頻電磁法對地幔中的揮發(fā)物、流體和部分熔融特別敏感[1]，而地磁測深(GDS)是低頻電磁法中一種重要的方法，成為研究地球地幔深部物理化學(xué)性質(zhì)以及地球動力學(xué)的重要技術(shù)手段。全球不同地磁臺站數(shù)據(jù)的一維反演結(jié)果表明了地幔電性的局部不均勻性[2]，此外，全球和半全球地磁數(shù)據(jù)的三維反演也揭示了地幔結(jié)構(gòu)的非一維性[3-4]。例如，Semenovo等[5]利用全球地磁臺站數(shù)據(jù)反演的結(jié)果揭示了地幔410～1 600 km，特別是670～900 km的深度范圍內(nèi)存在廣泛的橫向不均勻性。Kelbert等[6]通過全球中緯度地區(qū)的地磁數(shù)據(jù)反演得到全球三維電導(dǎo)率模型，并結(jié)合Huang等[7]和Yoshino等[8]的相關(guān)實驗室測量結(jié)果，認(rèn)為大洋板塊的俯沖作用將大量的水帶入轉(zhuǎn)換帶，使得轉(zhuǎn)換帶的電導(dǎo)率增大，說明地幔電性結(jié)構(gòu)的橫向變化很大。可見，地磁數(shù)據(jù)的一維反演無法獲得地幔精細(xì)的電性結(jié)構(gòu)，而通過地磁數(shù)據(jù)三維反演可以獲得比較精細(xì)的全球三維電導(dǎo)率模型。然而，有多種因素影響全球三維電導(dǎo)率模型的可靠性，例如，地磁臺站空間分布的不均勻性、正演求解精度不夠等問題[9]。

1 正演問題

1.1 數(shù)學(xué)物理模型

在長周期范圍，空氣層電導(dǎo)率小于10-10S/m[28]，介電常數(shù)為10-11F/m，則傳導(dǎo)電流在空氣和地層中占主導(dǎo)地位，忽略位移電流[16]。設(shè)時諧因子為eiωt，頻率域Maxwell方程的積分形式為

∮H·dl=?J·ndS，

(1)

∮E·dl=-?iωμH·dS，

(2)

J=σE。

(3)

在圖1球坐標(biāo)下，采用正交的曲面坐標(biāo)對地球進行離散。單元節(jié)點編號(i，j，k)對應(yīng)的球坐標(biāo)為(φ，θ，r)，磁場切向分量Hφ和Hθ分別定義為經(jīng)度和緯度方向，且它們的正方向分別指向東和南，垂直分量Hr定義在半徑方向。

邊界條件的設(shè)置對于計算精度非常重要?？諝鈱拥纳线吔鐟?yīng)該離地球足夠遠(yuǎn)，以保證地球內(nèi)部產(chǎn)生的二次磁場在源位置產(chǎn)生的干擾衰減掉。本文將計算域的外邊界設(shè)置在距地心10倍地球半徑R(R= 6 371 km)處，并且將地球以外的“空氣層”電導(dǎo)率設(shè)為10-10S/m；計算域的內(nèi)邊界為核幔邊界(CMB)，距離地表2 871 km，并且假設(shè)地核的電導(dǎo)率為無窮大。在此假設(shè)條件下，地核內(nèi)磁場切向分量為零，但是核幔邊界處磁場切向分量不為零。

圖1 球坐標(biāo)系下的磁場分量Fig.1 Magnetic field defined in a spherical coordinate system

1.2 球坐標(biāo)系下交錯網(wǎng)格有限差分

與傳統(tǒng)的直角坐標(biāo)系交錯網(wǎng)格有限差分[29]不同的是，本文采用Uyeshima等[16]的交錯網(wǎng)格有限差分策略，將整個球體剖分為表面彎曲的六面體棱柱，網(wǎng)格單元如圖2所示。假設(shè)單元網(wǎng)格內(nèi)的電導(dǎo)率值均勻。對第(i,j,k)棱柱，磁場三分量分別被定義在棱邊上，因此這些棱柱、棱邊和節(jié)點分別稱為H棱柱、H棱邊和H節(jié)點。φ(i)、θ(j)和r(k)分別是第(i,j,k)個H節(jié)點處的經(jīng)度、余緯度和到地心的距離。用L、M、N分別代表φ、θ、r方向的H棱柱數(shù)，則i∈[1,L]；j∈[1,M+1]；k∈[1,N+1]。用Hφ(i,j,k)、Hθ(i,j,k)和Hr(i,j,k)分別表示第(i,j,k)個H棱邊的磁場分量。用同樣方法定義電場。但電場法向分量在電導(dǎo)率不同的單元邊界是不連續(xù)的，因此在邊界的兩側(cè)需要分別定義，例如，與磁場Hθ(i,j,k)對應(yīng)的電場定義為Eθ(i,j-,k)和Eθ(i,j+,k)。

定義H棱邊φ、θ、r三個方向的長度分別為lφ(i,j,k)、lθ(i,j,k)和lr(i,j,k)，具體的計算公式為：

lφ(i,j,k)=r(k)sinθ(j)[φ(i+1)-φ(i)],

lθ(i,j,k)=r(k)[θ(j+1)-θ(j)],

lr(i,j,k)=r(k)-r(k+1)。

(4)

同樣，H棱柱三個方向的表面積表示為：

Sφ(i,j,k)=

sinθ(j)[φ(i+1)-φ(i)],

Sr(i,j,k)=r2(k)[cosθ(j)-cosθ(j+1)]

[φ(i+1)-φ(i)]。

(5)

利用式(4)和式(5)對Maxwell方程(1)、(2)、(3)分別進行離散(詳細(xì)過程參考文獻(xiàn)[16])，最后得到磁場解的三個分量滿足矩陣方程：

(6)

式中，矢量lH包含磁場的線積分。加入源后，代表邊界場值的常數(shù)移到方程的右端項，即包含在矢量lH中。Aθφ與磁場φ分量的線積分以及磁場θ分量的邊界值有關(guān)，其余類推。由Aθφ等子單元矩陣組成最終的系數(shù)矩陣A。系數(shù)矩陣A是稀疏對稱的、并且矩陣A除了其對角線元素外都是實數(shù)[16]。最后，對方程(6)式進行求解可以得到磁場分量。與Uyeshima等[16]采用的最小殘差加速(MRA)迭代法求解方程不同的是，本文利用求解速度更快的PARDISO求解器求解，避免了迭代求解過程中的散度校正[30]，同時也保證了求解精度和穩(wěn)定性。

左方H棱柱(黑線)的棱邊分別定義磁場的緯度、經(jīng)度和半徑方向分量，節(jié)點位于H棱柱中心的交錯棱柱(灰線)代表E棱柱。(-)代表H棱柱的外部，(+)代表電導(dǎo)率均勻棱柱的內(nèi)部?；谖墨I(xiàn)[16]。圖2 球坐標(biāo)下交錯網(wǎng)格有限差分網(wǎng)格剖分Fig.2 Mesh division geometry for the staggered-grid FD formulation in spherical coordinate system

1.3 全球地磁響應(yīng)函數(shù)

假設(shè)源為帶狀環(huán)形源[2]，根據(jù)Schultz等[31]地磁響應(yīng)函數(shù)的定義，由地表某點的磁場垂直和緯度方向分量，定義地磁測深中C響應(yīng)為

(7)

式中：tanθ為源空間結(jié)構(gòu)的補償項，因此一維情況下C響應(yīng)在地球表面的任意位置都是相同的；在空間某點上，Hr是指向地心的磁場垂直分量；Hθ是指向地磁南極的磁場余緯度分量。一般情況下C響應(yīng)的實部為正值，虛部為負(fù)值，并且可以在一定程度上反映深度與電阻率的關(guān)系[32]。

為了更好地描述地下電性結(jié)構(gòu)的橫向不均勻性，F(xiàn)ujii等[33]提出D響應(yīng)，其定義為

(8)

式中，Hφ是指向東的磁場沿經(jīng)度方向分量。D響應(yīng)在一維結(jié)構(gòu)下為零，且橫向不均勻性越強，D響應(yīng)的值越大。

2 精度驗證

三維有限差分模擬分別采用了36×18×54、72×36×78、180×90×98、360×180×118不同的粗細(xì)網(wǎng)格，發(fā)現(xiàn)不同網(wǎng)格的C響應(yīng)相對誤差在1%以內(nèi)。綜合計算精度和速度，下文數(shù)值模擬采用的網(wǎng)格為180×90×98，即地球在經(jīng)緯度方向被剖分為2°×2°、半徑方向被剖分為98層，其中空氣層為17層。計算機主頻為3.2 GHz，Intel處理器，內(nèi)存8 G。實際單個頻點計算占用內(nèi)存約3 G，單頻計算時間大約20 min。

2.1 一維模型

C響應(yīng)是地磁測深應(yīng)用最廣泛的響應(yīng)函數(shù)，本文驗證周期范圍從數(shù)小時到3 a的C響應(yīng)，選取35個對數(shù)均勻分布的頻率。一維模型電性結(jié)構(gòu)如Medin等[34]，其參數(shù)如表1所示。一維解析解采用李世文等[35]的層狀地球地表C響應(yīng)算法，驗證對比所用三維有限元正演結(jié)果來自Ribaudo等[15]。圖3a給出了本文有限差分算法數(shù)值精度驗證情況，可據(jù)文獻(xiàn)[34]。

a. 不同周期1D模型C響應(yīng)，正值是C的實部，負(fù)值是C的虛部；b. 本文有限差分解和一維解析解的相對誤差。圖3 球坐標(biāo)系中交錯網(wǎng)格有限差分算法的一維模型精度驗證結(jié)果Fig.3 Precision validation result of the 1D model for the staggered-grid FD formulation in spherical coordinate system

范圍深度/kmσ/(S/m)上地幔0～4000．01轉(zhuǎn)換帶400～8000．10下地幔800～28711．00地核>2871500000．00

以看出，有限差分得到的響應(yīng)與解析解以及有限元模擬響應(yīng)基本吻合。圖3b給出了有限差分?jǐn)?shù)值解和一維解析解的相對誤差，從圖可見：C響應(yīng)實部誤差都小于2%；虛部相對誤差大于實部誤差，最大也僅為4.8%，且比較大的誤差出現(xiàn)在1～10 d的小周期范圍以及大于100 d的超長周期。

2.2 雙半球模型

雙半球模型在全球電磁感應(yīng)正演模擬精度驗證中被廣泛采用[16, 22]，本文以南北(NS)雙半球模型(圖4)為例進行驗證。半球模型由一維電性結(jié)構(gòu)模型和地表無限薄的非均勻電導(dǎo)層組成。對于NS半球模型，北半球(余緯度0°～90°)薄層電導(dǎo)為20 000 S，南半球(余緯度90°～180°)薄層電導(dǎo)為20 S，薄層厚度設(shè)為12.6 km。

據(jù)文獻(xiàn)[22]。圖4 南北雙半球薄層模型Fig.4 North/south hemisphere model

圖5給出了本文有限差分算法在周期為4 d的磁場水平南向分量和垂直分量曲線。由圖5可知，本文計算結(jié)果與Kelbert等[21]的結(jié)果一致，說明了本文三維正演的精確性和可靠性。

3 三維數(shù)值模擬

與傳統(tǒng)數(shù)值模擬中異常體邊界的電性突變不同，“棋盤模型”(checkerboard model)不僅可以模擬地球橫向電導(dǎo)率的三維不均勻性，并且能很好地模擬局部異常體電阻率的連續(xù)性變化，對電磁法數(shù)值模擬更具有挑戰(zhàn)性[22]。棋盤模型通過球諧函數(shù)(SHFE)

據(jù)公司負(fù)責(zé)人介紹，廣東新華印刷是全省教材印制基地，承擔(dān)著由廣東省出版集團出版、發(fā)行的中小學(xué)教材70%的印刷任務(wù)。根據(jù)廣東中小學(xué)生人數(shù)情況，我們簡單做了推算，這意味著平均每年至少有上千萬名中小學(xué)生使用的是由廣東新華印刷印制的教材。如果從其開始承印教材的八九十年代算起，近三十年積累下來，無論從服務(wù)群體的數(shù)量還是服務(wù)內(nèi)容的質(zhì)量，廣東新華印刷在基礎(chǔ)教育的推廣方面均發(fā)揮了獨有的社會價值。

圖5 周期為4 d的南北雙半球模型響應(yīng)Fig.5 Response of North/South hemispheres at period 4 d

(9)

本文在表1所示一維模型的第二層，即在400～800 km深度的中低緯度進行球諧擾動。選擇球諧系數(shù)l=12、m=8來模擬橫向空間不均勻性，形成的“棋盤模型”如圖6所示。

圖7為周期為4 d時“棋盤模型”地球表面的磁場三分量的數(shù)值模擬結(jié)果。由圖7可見：Hφ對模型的橫向分辨率能力最好，尤其是Hφ的虛部對棋盤模型的邊界分辨能力更強，但是，Hφ反映不出來異常體的位置且在幅值上小于Hθ和Hr；Hθ的虛部響應(yīng)和棋盤模型(尤其是低阻異常)的位置對應(yīng)關(guān)系非常好，但實部分辨能力很差，且虛部數(shù)值比實部小很多；Hr的分辨能力較差，但實部和虛部對異常體的分辨能力相當(dāng)。

圖6 l=12，m=8時的棋盤模型Fig.6 Synthetic degree 12 order 8 checkerboard model

此外，由C響應(yīng)和D響應(yīng)的定義式(7)和式(8)可知，Hθ和Hr決定C響應(yīng)的大小，而D響應(yīng)中包含了對三維分辨能力更強的Hφ分量，這也印證了D響應(yīng)的橫向分辨率要優(yōu)于C響應(yīng)。但實測數(shù)據(jù)的D響應(yīng)容易受干擾影響[36]，所以現(xiàn)在對D響應(yīng)的應(yīng)用較少。隨著儀器技術(shù)和信號處理技術(shù)的發(fā)展，D響應(yīng)未來可以被廣泛應(yīng)用?？傊?，磁場分量能夠分辨出棋盤模型的異常體的大小和準(zhǔn)確的位置。

4 討論與結(jié)論

分別通過一維模型和雙半球模型的精度驗證，說明了本文數(shù)值模擬的正確性和精度；三維“棋盤模型”結(jié)果表明地表磁場分量對異常體具有很好的分辨能力，Hφ分量對異常體的邊界分辨能力更強，而Hθ分量和異常體的位置對應(yīng)關(guān)系非常好。同時，D響應(yīng)比C響應(yīng)對三維異常體分辨能力更強，隨著技術(shù)的發(fā)展，D響應(yīng)同樣具有重要的研究價值。

針對超長周期的地磁測深三維反演，直角坐標(biāo)系無法適用，本文的球坐標(biāo)系三維正演為地磁測深三維反演提供了技術(shù)基礎(chǔ)，對利用地磁數(shù)據(jù)獲得中國或全球地幔電性結(jié)構(gòu)，從而研究地球演化和地球動力學(xué)過程具有重要意義。

致謝：美國地質(zhì)調(diào)查局的Anna Kelbert研究員在本文數(shù)值模擬研究工作中提供了熱情幫助，Intel Fortran數(shù)學(xué)庫MKL提供了PARDISO求解器，在此一并致謝。

[1] Karato S I. Remote Sensing of Hydrogen in Earth’s Mantle[J]. Reviews in Mineralogy & Geochemistry, 2006, 62(1):343-375.

[2] Schultz A. On the Vertical Gradient and Associated Heterogeneity in Mantle Electrical Conductivity[J]. Phys Earth Planet Int，1990, 64(1): 68-86.

[3] Utada H, Koyama T, Obayashi M, et al. A Joint Interpretation of Electromagnetic and Seismic Tomography Models Suggests the Mantle Transition Zone Below Europe is Dry[J]. Earth Planet Sci Lett, 2009, 281(3/4): 249-257.

[4] Shimizu H, Utada H, Baba K, et al. Three-Dimen-sional Imaging of Electrical Conductivity in the Mantle Transition Zone Beneath the North Pacific Ocean by a Semi-Global Induction Study[J]. Phys Earth Planet Int， 2010, 183(1/2)：252-269.

[5]Semenov A，Kuvshinov A. Global 3-D Imaging of Mantle Conductivity Based on Inversion of Observatory C-Responses: II: Data Analysis and Results[J]. Geophysical Journal International，2012, 191(3): 965-992.

[6] Kelbert A，Schultz A，Egbert G. Global Electromag-netic Induction Constraints on Transition-Zone Water Content Variations[J]. Nature, 2009, 460：1003-1006.

[7] Huang X，Xu Y，Karato S I. Water Content in the Transition Zone from Electrical Conductivity of Wadsleyite and Ringwoodite[J]. Nature, 2005, 434: 746-749.

[8] Yoshino T，Manthilake G，Matsuzaki T，et al. Dry Mantle Transition Zone Inferred from the Conductivity of Wadsleyite and Ringwoodite[J]. Nature, 2008, 451: 326-329.

[9] Sun J，Egbert G. A Thin-Sheet Model for Global Electromagnetic Induction[J]. Geophysical Journal International，2012, 189(1): 343-356.

[10]Schmucker U. A Spherical Harmonic Analysis of Solar Daily Variations in the Years 1964-1965: Response Estimates and Source Fields for Global Induction: I: Methods[J]. Geophysical Journal International, 1999, 136(2): 439-454.

[11] Koyama T, Khan A, Kuvshinov A. Three-Dimen-sional Electrical Conductivity Structure Beneath Australia from Inversion of Geomagnetic Observatory Data: Evidence for Lateral Variations in Transition-Zone Temperature, Water Content and Melt[J]. Geophysical Journal International, 2014, 196(3): 1330.

[12] Kuvshinov A. 3-D Global Induction in the Oceans and Solid Earth: Recent Progress in Modeling Magnetic and Electric Fields from Sources of Magnetospheric, Ionospheric and Oceanic Origin[J]. Surveys in Geophysics, 2008, 29(2): 139-186.

[13] Everett M E，Schultz A. Geomagnetic Induction in a Heterogeneous Sphere: Azimuthally Symmetric Test Computations and the Response of an Undulating 600-km Discontinuity[J]. J Geophys of Research Atmospheres, 1996, 101(B2)：2765-2783.

[14] Yoshimura R，Oshiman N. Edge-Based Finite Element Approach to the Simulation of Geoelectromagnetic Induction on a 3-D Sphere[J]. Geophys Res Lett, 2002, 29(3): 1-4.

[15] Ribaudo J T，Constable C G，Parker R L. Scripted Finite Element Tools for Global Electromagnetic Induction Studies[J]. Geophysical Journal International, 2012, 188(2)：435-446.

[16] Uyeshima M，Schultz A. Geoelectromagnetic Induc-tion in a Heterogeneous Sphere: A New Three-Dimensional Forward Solver Using a Conservative Staggered-Grid Finite Difference Method[J]. Geophysical Journal International，2000, 140(3)：636-650.

[17]Kelbert A, Egbert G D, Schultz A. Non-Linear Conjugate Gradient Inversion for Global EM Induction: Resolution Studies[J]. Geophysical Journal International, 2008, 173(2): 365-381.

[18] Weiss C J. Triangulated Finite Difference Methods for Global-Scale Electromagnetic Induction Simulations of Whole Mantle Electrical Heterogeneity[J]. Geochemistry Geophysics Geosystems，2010，11(11)：129-143.

[19] Tarits P, Mande A M. The Heterogeneous Electrical Conductivity Structure of the Lower Mantle[J]. Physics of the Earth & Planetary Interiors, 2010, 183(1/2): 115-125.

[20]Martinec Z. Spectral-Finite Element Approach to Three-Dimensional Electromagnetic Induction in a Spherical Earth[J]. Geophysical Journal International, 1999, 136(1): 229-250.

[22] Kelbert A, Kuvshinov A, VelímskJ, et al. Global 3-D Electromagnetic Forward Modelling: A Benchmark Study[J]. Geophysical Journal Internationa，2014, 197(2)：785-814.

[23] 陳伯舫. 中國東南部地幔高導(dǎo)層的埋藏深度[J]. 地震學(xué)報, 1986,8(2): 172-178.

Chen Bofang. Buried Depth of the High Conductivity Layer Beneath the Region of South-East China[J]. Acta Seismologica Sinica,1986, 8(2): 172-178.

[24] 杜興信, 魯秀玲. 用單臺Z/H法研究陜西地區(qū)深部電導(dǎo)率[J]. 內(nèi)陸地震, 1994,8(4): 333-338.

Du Xingxin, Lu Xiuling. Investigation on Deep Electroconductivity of Shaanxi Region by Means of Z/H Method with Single Station[J]. Inland Earthquake,1994, 8(4): 333-338.

[25] 范國華, 姚同起, 顧左文，等. 利用磁梯度法研究我國地幔導(dǎo)電率[J]. 地震學(xué)報, 1997,19(2): 164-173.

Fan Guohua, Yao Tongqi, Gu Zuowen, et al. Research on Mantle Conductivity of China with Gradient Method[J]. Acta Seismologica Sinica, 1997, 19(2): 164-173.

[26] 趙國澤，湯吉，梁競閣，等. 用大地電磁網(wǎng)法在長春等地探測上地幔電導(dǎo)率結(jié)構(gòu)[J]. 地震地質(zhì)，2001, 23(2)：143-152.

Zhao Guoze，Tang Ji，Liang Jingge，et al. Measurement of Network-MT in Two Areas of NE China for Study of Upper Mantle Conductivity Structure of the Back-Arc Region[J]. Seismology and Geology，2001, 23(2)：143-152.

[27] 徐光晶，湯吉，黃清華，等. 華北地區(qū)上地幔及過渡帶電性結(jié)構(gòu)研究[J]. 地球物理學(xué)報，2015, 58(2)：566-575.

Xu Guangjing，Tang Ji，Huang Qinghua，et al. Study on the Conductivity Structure of the Upper Mantle and Transition Zone Beneath North China[J]. Chinese Journal of Geophysics，2015, 58(2)：566-575.

[28] Rokityansky I I. Geoelectromagnetic Investigation of the Earth’s Crust and Mantle[M]. Berlin: Springer-Verlag，1982.

[29] 李大俊，翁愛華，楊悅，等. 地-井瞬變電磁三維交錯網(wǎng)格有限差分正演及響應(yīng)特性[J]. 吉林大學(xué)學(xué)報(地球科學(xué)版), 2017, 47(5): 1552-1561.

Li Dajun, Weng Aihua, Yang Yue,et al. Three-Dimension Forward Modeling and Characteristics for Surface-Borehole Transient Electromagnetic by Using Staggered-Grid Finite Difference Method[J]. Journal of Jilin University (Earth Science Edition), 2017, 47(5): 1552-1561.

[30] Smith J T. Conservative Modeling of 3-D Electromag-netic Fields: Part 1: Properties and Error Analysis[J]. Geophysics, 1996, 61(5):1308-1318.

[31] Schultz A，Larsen J C. On the Electrical Conductivity of the Mid-Mantle: I: Calculation of Equivalent Scalar Magnetotelluric Response Functions[J]. Geophysical Journal International, 1987, 88(3):733-761.

[32] Schmucker U. Substitute Conductors for Electromag-netic Response Estimates[J]. Pure and Applied Geophysics, 1987, 125(2): 341-367.

[33] Fujii I, Schultz A. The 3D Electromagnetic Response of the Earth to Ring Current and Auroral Oval Excitation[J]. Geophysical Journal International, 2002, 151(3): 689-709.

[34] Medin A E，Parker R L，Constable S. Making Sound Inferences from Geomagnetic Sounding[J]. Physics of the Earth & Planetary Interiors, 2007, 160(1): 51-59.

[35] 李世文，翁愛華，李建平，等. 淺部約束的地磁測深C-響應(yīng)一維反演[J]. 地球物理學(xué)報, 2017，60(3):1201-1210.

Li Shiwen, Weng Aihua, Li Jianping, et al. 1-D Inversion of C-Response Data from Geomagnetic Depth Sounding with Shallow Resistivity Constraint [J]. Chinese Journal of Geophysics, 201760(3): 1201-1210.

[36] Shimizu H, Koyama T, Baba K. Three-Dimensional Geomagnetic Response Functions for Global and Semi-Global Scale Induction Problems[J]. Geophysical Journal International, 2009, 178(1): 123-144.