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

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

0 引言

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

1 正演問(wèn)題

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

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

∮H·dl=?J·ndS，

(1)

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

(2)

J=σE。

(3)

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

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

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

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

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

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

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棱柱三個(gè)方向的表面積表示為：

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)對(duì)Maxwell方程(1)、(2)、(3)分別進(jìn)行離散(詳細(xì)過(guò)程參考文獻(xiàn)[16])，最后得到磁場(chǎng)解的三個(gè)分量滿(mǎn)足矩陣方程：

(6)

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

左方H棱柱(黑線)的棱邊分別定義磁場(chǎng)的緯度、經(jīng)度和半徑方向分量，節(jié)點(diǎn)位于H棱柱中心的交錯(cuò)棱柱(灰線)代表E棱柱。(-)代表H棱柱的外部，(+)代表電導(dǎo)率均勻棱柱的內(nèi)部?；谖墨I(xiàn)[16]。圖2 球坐標(biāo)下交錯(cuò)網(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ù)的定義，由地表某點(diǎn)的磁場(chǎng)垂直和緯度方向分量，定義地磁測(cè)深中C響應(yīng)為

(7)

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

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

(8)

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

2 精度驗(yàn)證

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

2.1 一維模型

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

a. 不同周期1D模型C響應(yīng)，正值是C的實(shí)部，負(fù)值是C的虛部；b. 本文有限差分解和一維解析解的相對(duì)誤差。圖3 球坐標(biāo)系中交錯(cuò)網(wǎng)格有限差分算法的一維模型精度驗(yàn)證結(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ù)值解和一維解析解的相對(duì)誤差，從圖可見(jiàn)：C響應(yīng)實(shí)部誤差都小于2%；虛部相對(duì)誤差大于實(shí)部誤差，最大也僅為4.8%，且比較大的誤差出現(xiàn)在1～10 d的小周期范圍以及大于100 d的超長(zhǎng)周期。

2.2 雙半球模型

雙半球模型在全球電磁感應(yīng)正演模擬精度驗(yàn)證中被廣泛采用[16, 22]，本文以南北(NS)雙半球模型(圖4)為例進(jìn)行驗(yàn)證。半球模型由一維電性結(jié)構(gòu)模型和地表無(wú)限薄的非均勻電導(dǎo)層組成。對(duì)于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的磁場(chǎng)水平南向分量和垂直分量曲線。由圖5可知，本文計(jì)算結(jié)果與Kelbert等[21]的結(jié)果一致，說(shuō)明了本文三維正演的精確性和可靠性。

3 三維數(shù)值模擬

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

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

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

(9)

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

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

圖6 l=12，m=8時(shí)的棋盤(pán)模型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)中包含了對(duì)三維分辨能力更強(qiáng)的Hφ分量，這也印證了D響應(yīng)的橫向分辨率要優(yōu)于C響應(yīng)。但實(shí)測(cè)數(shù)據(jù)的D響應(yīng)容易受干擾影響[36]，所以現(xiàn)在對(duì)D響應(yīng)的應(yīng)用較少。隨著儀器技術(shù)和信號(hào)處理技術(shù)的發(fā)展，D響應(yīng)未來(lái)可以被廣泛應(yīng)用?？傊?，磁場(chǎng)分量能夠分辨出棋盤(pán)模型的異常體的大小和準(zhǔn)確的位置。

4 討論與結(jié)論

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

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

致謝：美國(guó)地質(zhì)調(diào)查局的Anna Kelbert研究員在本文數(shù)值模擬研究工作中提供了熱情幫助，Intel Fortran數(shù)學(xué)庫(kù)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] 陳伯舫. 中國(guó)東南部地幔高導(dǎo)層的埋藏深度[J]. 地震學(xué)報(bào), 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] 杜興信, 魯秀玲. 用單臺(tái)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] 范國(guó)華, 姚同起, 顧左文，等. 利用磁梯度法研究我國(guó)地幔導(dǎo)電率[J]. 地震學(xué)報(bào), 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] 趙國(guó)澤，湯吉，梁競(jìng)閣，等. 用大地電磁網(wǎng)法在長(zhǎng)春等地探測(cè)上地幔電導(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ū)上地幔及過(guò)渡帶電性結(jié)構(gòu)研究[J]. 地球物理學(xué)報(bào)，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] 李大俊，翁愛(ài)華，楊?lèi)偅? 地-井瞬變電磁三維交錯(cuò)網(wǎng)格有限差分正演及響應(yīng)特性[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(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] 李世文，翁愛(ài)華，李建平，等. 淺部約束的地磁測(cè)深C-響應(yīng)一維反演[J]. 地球物理學(xué)報(bào), 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.