2011, Vol. 54
2. 中国地质大学地下信息探测技术与仪器教育部重点实验室,北京 100083;
3. 中国地质大学(北京)地球物理与信息技术学院,北京 100083
2. Key Laboratory of Geo-detection (China University of Geosciences), Ministry of Education, Beijing 100083, China;
3. School of Geophysics and Information Technology, China University of Geosciences, Beijing 100083, China
在利用天然电磁场进行勘探的现代地球物理方法中，采集垂直磁场分量的主要有两种:①同时采集电场和磁场信号的大地电磁测深法(MTS):观测两个水平方向的电场，两个水平方向和一个垂直方向的磁场；②只采集磁场信号的地磁测深法(GDS):观测两个水平方向和一个垂直方向的磁场[1， 2].在大地电磁测深理论和地磁测深理论中，倾子(也称垂直磁场转换函数)描述了磁场的垂直分量与水平分量之间的关系.在以z方向垂直向下的右手坐标系中，倾子矢量T的定义为
|
(1) |
式中，Hz，Hx和Hy分别为垂直磁场分量、x和y方向的磁场分量；Tzx和Tzy分别表示倾子T在x方向和y方向的分量.
大地电磁测深法的主要参数是张量阻抗参数(即水平电场分量和水平磁场分量之间的关系转换函数)和由阻抗整理出来的视电阻率和相位参数.当大地电磁测深工作同时采集三分量磁场信号时，可根据(1)式整理出倾子参数.相对视电阻率、相位和阻抗参数而言，倾子参数在大地电磁测深法中只是作为辅助参数使用.原因在于:倾子的异常值很小，且容易受到环境噪声的干扰，因而一般认为倾子资料信噪比低，利用倾子资料进行定量解释困难.因此，大地电磁测深和地磁测深所获得的倾子资料多被用于定性判别地下构造的维度和走向[3].
随着数字电子技术的发展和仪器灵敏度的提高，现代大地电磁野外工作，可以获得质量较高的垂直磁场信号并通过垂直磁场整理出倾子资料，特别是使用三分量磁探头进行长周期电磁测量的时候常采集垂直磁场分量.同时，现代地震地磁观测技术迅速发展，使得地磁观测资料对固体地球物理学、空间物理、地球内部电导率结构和地球动力学研究，以及空间天气学、地震学、地震监测预报研究等具有重要意义.多个国家在世界各地建立了很多采集三个正交磁场分量的地磁观测台阵.如何充分利用这些地磁观测资料获得地下电性结构并在地质勘探、环境监测、航空航海导航等生产实践和国防建设等领域予以应用是地球物理学家们需要考虑的问题[4-6].因此，研究者们研究倾子数据并发展倾子数据的定量解释方法.陈小斌等对倾子矢量的图示分析和应用进行了研究[7].陈清礼等对埋藏球体的倾子响应特征进行了分析[8].倾子或视倾子资料被应用于二维大地电磁资料的断裂解释[9， 10].Ledo等将倾子资料应用于大地电磁静态效应问题的研究[11].Berdichevsky等通过对二维地质模型倾子资料的分析研究，得出结论:与视电阻率、阻抗数据相比，倾子数据受大地电磁静态效应的影响小，使用倾子资料参与反演可以充分地提高地质模型的可靠性[1].许多二维反演代码[12， 13]实现了倾子资料的二维反演，并且将倾子数据用于大地电磁二维资料的反演解释[14-16].
对大地电磁视电阻率、相位或阻抗数据进行三维反演的研究有很多[17-28].在对大地电磁倾子数据和共轭梯度算法深入分析的基础上，我们实现了倾子资料三维共轭梯度反演算法.与大地电磁测深类似，地磁测深(地磁观测台阵观测的三分量磁场)也可以根据公式(1)整理出相同的倾子资料.因此，该反演算法也适用于地磁测深数据的解释.本文的目的是基于倾子资料的三维共轭梯度反演研究，探讨利用大地电磁测深或地磁测深的倾子资料进行三维反演获得地下三维电性结构的定量解释方法.文中，我们将阐述一种使用倾子资料进行反演的三维共轭梯度反演算法.首先讨论计算倾子响应的三维正演问题，其次介绍反演理论，包括目标函数和“拟正演"问题的计算，最后给出理论模型合成数据的三维反演结果.
2 三维正演由于三维反演中我们使用倾子参加反演计算，因此三维正演问题需要考虑如何从地下复杂的三维地质体中数值模拟得到精确的地表倾子响应，Tzx和Tzy.
参考谭捍东的思路，对于大地电磁场源S 的作用，总可以分解成两个正交的场源SX 和场源SY 等效作用的结果[29].在场源SX 的作用下，地表处产生的磁场分量分别记为Hx1、Hy1、Hz1，在场源SY 的作用下，地表处产生的磁场分量分别记为Hx2、Hy2、Hz2，则由它们计算倾子的表达式为
|
(2) |
从公式(2)可以看出，要想得到地表的倾子响应，首先需要获得地表的磁场值.对于给定的一种场源极化模式，地表观测点j处的磁场值Hj可以通过磁场值列向量H计算得到[19， 22]:
|
(3) |
它们之间通过内插向量hgjT转化.对于Hj的三个方向分量Hjx、Hjy和Hjz，方程(3)分别对应为
|
|
(4) |
从公式(4)可以看出，要想得到地表的磁场值，还需要先获得地下磁场值列向量H.而地下磁场值列向量H则可以通过采用交错采样有限差分法[19， 30]数值模拟得到.用交错采样网格(图 1)将麦克斯韦方程组的积分形式:
|
(5) |
|
图 1 编号为(i，j，k)的网格单元 Fig. 1 The (i，j，k) grid cell |
进行离散化，可得到关于地下各网格单元采样点处的磁场值H的正演方程:
|
(6) |
其中，K为对称的大型稀疏系数矩阵；H为待求解地下各网格单元采样点处的磁场三分量组成的列向量；s为与源及边界场值有关的列向量.如果，把场源SX 的作用下待求解的磁场值和方程的右端列向量分别记为H1 和s1 ，把场源SY 的作用下待求解的磁场值和方程的右端列向量分别记为H2 和s2 ，那么可以形成两个正演方程:KH1 =s1 ，KH2 =s2.解两个正演方程可分别得到地下各网格单元采样点处的两种极化模式磁场值H1 和H2.
将得到的两种极化模式磁场值H1 和H2 代入方程(3)，可求出地表观测点处两种极化模式的磁场值Hx1、Hy1、Hz1、Hx2、Hy2、Hz2.再将地表磁场值代入计算倾子的表达式(2)，可求得三维模型的倾子响应.
为了验证三维正演算法的正确性，我们设计了一个二维棱柱体地电模型.二维棱柱体的走向为y方向，几何尺寸为6km×3km, 电阻率为10Ωm, 顶面埋深为3km, 围岩电阻率为100Ωm.用三维正演算法计算出该棱柱体的倾子响应，再用二维有限差分法对该棱柱体进行数值模拟，得到倾子响应.图 2对比了频率分别为1Hz和0.1Hz时两种数值模拟计算结果.图 2a显示倾子Tzx的实部，图 2b显示倾子Tzx的虚部.“￮"表示频率为1Hz的三维数值解，“Δ "表示频率为0.1 Hz的三维数值解，“×"表示频率为1Hz的二维有限差分数值解，“+"表示频率为0.1Hz的二维有限差分数值解.由图可见，无论是倾子的实部还是倾子的虚部，两者的结果都拟合得非常好.该结果表明三维正演算法的计算结果是准确可靠的.
|
图 2 二维棱柱体模型倾子响应的二维有限差分数值解和三维数值解的对比 Fig. 2 The tipper response comparison between the 2D finite difference modeling results and the 3D modeling results generated from the 2D prism model |
倾子三维共轭梯度反演算法的目标函数定义为
|
(7) |
其中，Tobs表示观测倾子数据；F(m)为求取倾子响应的正演函数；V为数据方差；λ 为正则化参数；L为简单的二次差分拉普拉斯算子，m0 为先验模型.目标函数的梯度相应可表示为
|
(8) |
这里，A表示雅可比矩阵，数据误差向量:e=Tobs-F(m).
3.2 雅可比矩阵的计算---“拟正演"问题反演问题中，最关键的部分在于求取雅可比矩阵.在计算出雅可比矩阵后，便可获得每次模型更新所需的模型修正量，实现模型的更新迭代.倾子三维共轭梯度的反演过程避开了实际计算和存储雅可比矩阵的每个元素，而是通过计算雅可比矩阵的转置与一个向量的乘积以及雅可比矩阵与另一个向量的乘积直接获得模型修正量，其算法流程与阻抗三维共轭梯度反演的算法流程[19]类似，整个反演过程的主要计算量在于计算A(mref )TV-1e和A(mref )p的值.这里，mref 表示参考模型，p表示搜索方向.下面，我们介绍倾子三维共轭梯度反演如何求取A(mref )TV-1e和A(mref )p.
正演方程(6)两端同时对模型参数mk求偏导数，则有
|
(9) |
将公式(4)代入方程(9)并由∂hgjT/∂mk= 0可得:
|
(10) |
方程(2)两端同时对第k个模型参数mk求偏导数，可得:
|
(11) |
将方程(10)代入方程(11)，则第k个模型参数mk的扰动所引起的第j个观测点处倾子的改变量∂Tzxj/∂mk和∂Hzyj/∂mk可以写成如下形式:
|
其中，
|
(12) |
由于使用倾子资料参与反演，因此A(mref )TV-1e相应地可表示为
|
(13) |
其中，N为倾子数据的总个数:测点数× 频率数×4；k=1，2，…，M表示模型参数.因此，方程(12)中的∂Tzxj/∂mk、∂Hzyj/∂mk所求的就是雅可比矩阵A的元素.将式(12)代入式(13)并由K为对称阵可得:
|
(14) |
其中，1gn、2gn由式(12)确定.例如:如果Hzyj为第n个数据，则n1gT=jty1gT，2gnT=2gjtyT.
令
|
(15) |
|
(16) |
则
|
(17) |
|
(18) |
|
(19) |
把v1 视为解向量，$\sum\limits_{n = 1}^N {{g_n}\left( {{V^{ - 1}}e} \right)} n$ 视为方程(17)右端向量，式(17)是一个和式(6)类似的正演方程，我们称之为“拟正演"方程.通过解一次“拟正演"方程，可以得到v1.同样对于式(18)，通过解一次“拟正演"方程，可以得到v2.把v1 和v2 代入式(19)，可以计算出A(mref )TV-1e.
同理，对于A(mref )p，经过推导可得:
|
(20) |
令
|
(21) |
|
(22) |
则有
|
(23) |
|
(24) |
|
(25) |
把t1 视为解向量，$\sum\limits_{k = 1}^M {\left( {\partial K/\partial {m_k}{H_1}} \right)} $视为方程(23)右端向量，式(23)也是一个和式(6)类似的正演方程，我们称之为“拟正演"方程.通过解一次“拟正演"方程，可以得到t1.同样对于式(24)，通过解一次“拟正演"方程，可以得到t2.把t1 和t2 代入式(25)，可以计算出A(mref )p.这样通过解两次“拟正演"问题，也可以得到A(mref )p.
4 合成算例基于上面的理论和公式，我们编程开发了倾子资料三维共轭梯度反演程序.为了验证三维反演算法的正确性，我们设计了一个单棱柱体地电模型.
电阻率为10Ωm, 大小为6km×6km×3km, 顶面埋深为3km 的棱柱体埋藏于电阻率为100Ωm 的地下均匀半空间.三维网格剖分为38×38×35(含10个空气层).用倾子资料三维共轭梯度反演程序的正演代码部分计算出单棱柱体模型在地表所有剖分网格单元中心点处产生的频率为3.3 Hz、1 Hz、0.33Hz和0.1Hz的倾子数据.图 3给出了频率为0.1Hz时单棱柱体模型在地表处产生的倾子数据的平面等值线图，该图显示了Tzx和Tzy响应的平面分布特征.从图中可以看出倾子响应Tzx和Tzy的异常极值都位于棱柱体边界在地表的投影处附近:Tzx响应以x=0 处为分界线，在棱柱体边界x=-3km和x=3km 的附近分别出现极大值或极小值，并且异常值从边界处往两侧逐渐减弱，在x=0和x=∞处Tzx=0；而Tzy响应以y=0处为分界线，在棱柱体边界y=-3km 和y=3km 的附近分别出现极大值或极小值，并且异常值从边界处往两侧逐渐减弱，在y=0和x=∞处Tzy=0.
|
图 3 频率为0.1 Hz时单棱柱体模型产生的Tzx和Tzy响应的平面等值线图 黑色实线表示棱柱体的边界，从左到右分别对应Tzx实部、Tzx虚部、Tzy实部和Tzy虚部. Fig. 3 Tipper response Tzx and Tzy at 0.1 Hz generated from one prism model The black solid lines show the prism margins.Lett to right are the real part of ，the imaginary part of ，the real part of Tzy and the imaginary part of Tzy. |
在地表 784个测点处的倾子数据中加入3%高斯随机误差后用倾子资料三维共轭梯度反演程序在PC 机上进行反演.PC 机的配置为:Intel(R)core(TM)2Duo处理器，主频2.2GHz, 内存2.0G.经过116次反演迭代，耗时111h2 min, 数据的拟合方差从初始值16.69收敛到4.22迭代结束(见图 4).反演的结果见图 5的第二行.图 5 中同时给出了真实模型(图 5的第一行)和对相同测点处的张量阻抗数据进行三维反演的反演结果(图 5 的第三行)，其中张量阻抗数据也加入3% 的高斯随机误差.与阻抗数据的反演结果(电阻率最低值达到22 Ωm)相比，倾子数据的反演结果虽然电阻率值偏高(最低值为33.8Ωm)，棱柱体的中心位置也稍微往地表方向偏离，但基本上恢复出了棱柱体的形状、大小和位置.相对而言，在限制棱柱体横向边界方面，倾子数据的反演结果要优于阻抗数据的反演结果.
|
图 4 数据拟合方差随迭代次数收敛图 Fig. 4 Chi-square error versus iteration number |
|
图 5 对单棱柱体模型产生的阻抗和倾子数据进行三维反演的结果 图中第一行表示真实模型，第二行为倾子数据的反演结果，第三行为阻抗数据的反演结果.黑色虚线表示棱柱体的边界.第一列为深度4.25km 处的水平截面图，第二列为x=0km 处沿y方向的垂直断面图，第三列为y=0km 处沿x方向的垂直断面图. Fig. 5 Results from the 3D inversion of impedance and tipper data generated from one prism model The top row shows the test model, the second row shows the results from the inversion of tipper data, and the third row shows the results from the inversion of impedance data.The black dashed lines show the prism margins.The first column shows the horizontal slices at 4.25 km depth, the second column shows the vertical slices at x = 0 km along the y axis, and the third column shows the vertical slices aty=0 km along the x axis. |
频率为0.33Hz时输入的倾子数据和反演模型响应倾子数据的比较见图 6.从图中可以看出，输入的倾子数据和反演模型响应倾子数据吻合得比较好.
|
图 6 频率为0.33 Hz时输入的倾子数据和反演模型响应倾子数据的比较图 黑色实线表示棱柱体的边界.从左到右分别对应Tzx实部、Tzx虚部、Tzy实部和Tzy虚部.第一行(Calculated)表示反演模型响应，第二行(Observed)表示输入倾子数据. Fig. 6 The comparison between input tipper data at 0.33 Hz and the tipper response from the inversion result model The black solid lines show the prism margins.Lett to right are the real part of ’frzx, the imaginary part of ， the real part of and theimaginary part of Wzy.The first row (Calculated) shows the tipper response from the inversion result model and the second row(Observed) shows the input tipper data. |
我们实现了倾子资料三维共轭梯度反演算法.该反演算法可用于对大地电磁测深和地磁测深所整理出的倾子资料进行三维定量反演，获得地下三维模型的电阻率结构.倾子资料的三维反演过程避开了实际计算和存储雅可比矩阵的每个元素，而是通过计算雅可比矩阵的转置与一个向量的乘积以及雅可比矩阵与另一个向量的乘积直接获得每次模型更新所需的模型修正量，实现模型的更新迭代，有效地减少了反演的计算量，提高了反演效率.合成数据的反演算例验证了所实现的倾子资料三维共轭梯度反演算法的有效性和稳定性.
在对合成的倾子数据进行三维反演的过程中，我们发现初始模型的选择对反演结果具有影响:当初始模型的电阻率值接近真实模型的背景值时，能得到较好的三维反演结果；当初始模型的电阻率值偏离真实模型的背景值较多时，目标体和围岩的电阻率值的还原效果会差一些，但还是可以较好还原出目标体的大小、形状和位置等.因此，建议在对实测倾子数据进行三维反演时可先采取其他方法(例如，钻孔资料、地质资料、阻抗或视电阻率的反演结果等)获得工区的背景电阻率值，或结合其他资料对倾子数据的反演结果进行解释.
| [1] | Berdichevsky M N, Dmitriev V I, Golubtsova N S, et al. Magnetovariational sounding: new possibilities, Izvestiya. Physics of the Solid Earth , 2003, 39: 701-727. |
| [2] | Egbert G D, Booker J R. Robust estimation of geomagnetic transfer functions. J. R. astr. Soc. , 1986, 87: 173-194. DOI:10.1111/gji.1986.87.issue-1 |
| [3] | 胡文宝, 苏朱刘, 陈清礼, 等. 倾子资料的特征及应用. 石油地球物理勘探 , 1997, 32(2): 202–213. Hu W B, Su Z L, Chen Q L, et al. The characteristics and application of tipper data. Oil Geophysical Prospecting (in Chinese) , 1997, 32(2): 202-213. |
| [4] | 张贵宾, 申宁华. 地磁三分量测深方法与长春、通化壳幔电性结构. 地球物理学报 , 1994, 37(S2): 296–303. Zhang G B, Shen N H. Geomagnetic three-components depth sounding and the electrical conductivity of the crust and upper mantle in Changchun and Tonghua. Chinese J. Geophys. (in Chinese) , 1994, 37(S2): 296-303. |
| [5] | 腾云田.现代地震地磁观测技术研究.北京:中国地震局地球物理研究所,2001. Teng Y T. Study on modern technique of seismic and geomagnetic observation (in Chinese). Beijing: Institute of Geophysics, China Earthquake Administration, 2001 |
| [6] | 杜兴信, 鲁秀玲. 陕西地区单台Z/H 地磁测深研究. 地震地磁观测与研究 , 1995, 16(1): 27–34. Du X X, Lu X L. Research on the geomagnetic deep sounding in Shanxi region by single station Z/H method. Seismological and Geomagnetic Observation and Research (in Chinese) , 1995, 16(1): 27-34. |
| [7] | 陈小斌, 赵国泽, 詹艳, 等. 磁倾子矢量的图示分析及其应用研究. 地学前缘 , 2004, 11(4): 626–636. Chen X B, Zhao G Z, Zhan Y, et al. Analysis of tipper visual vectors and its application. Earth Science Frontiers (in Chinese) , 2004, 11(4): 626-636. |
| [8] | 陈清礼, 胡文宝, 李金铭, 等. 埋藏球体的倾子响应特征分析. 石油天然气学报 , 2007, 29(3): 75–78. Chen Q L, Hu W B, Li J M, et al. Characteristics of tipper response of buried sphere. Journal of Oil and Gas Technology (in Chinese) , 2007, 29(3): 75-78. |
| [9] | 余年, 王绪本, 阚瑷珂, 等. 倾子和视倾子的研究及在断裂解释中的应用. 工程地球物理学报 , 2007, 4(4): 275–281. Yu N, Wang X B, Han A K, et al. Study on tipper and apparent tipper & application in fault interpretation. Chinese Journal of Engineering Geophysics (in Chinese) , 2007, 4(4): 275-281. |
| [10] | 于鹏, 吴健生, 王家林, 等. 利用长周期MT数据研究沪浙地区深部断裂结构. 同济大学学报 , 2008, 36(4): 549–554. Yu P, Wu J S, Wang J L, et al. Long period magnetotelluric data-based study on deep fault structure in Shanghai and Zhejiang area. Journal of Tongji University (in Chinese) , 2008, 36(4): 549-554. |
| [11] | Ledo J, Gabas A, Marcuello A. Static shift leveling using geomagnetic transfer functions. Earth Planets Space , 2002, 54: 493-498. DOI:10.1186/BF03353040 |
| [12] | Siripunvaraporn W, Egbert G. An efficient data-subspace inversion method for 2-D magnetotelluric data. Geophysics , 2000, 65: 791-803. DOI:10.1190/1.1444778 |
| [13] | Rodi W, Mackie R L. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion. Geophysics , 2001, 66: 174-187. DOI:10.1190/1.1444893 |
| [14] | Gabas A, Marcuello A. The relative influence of different types of magnetotelluric data on joint inversions. Earth Planets Space , 2003, 55: 243-248. DOI:10.1186/BF03351755 |
| [15] | Tuncer V, Unsworth M J, Siripunvaraporn W, et al. Exploration for unconformity type uranium deposits with audio-magnetotelluric data: A case study from the McArthur River Mine, Saskatchewan (Canada). Geophysics , 2006, 71(6): 201-209. DOI:10.1190/1.2348780 |
| [16] | Becken M, Ritter O, Burkhardt H. Mode separation of magnetotelluric responses in three-dimensional environments. Geophys. J. Int , 2008, 172: 64-86. |
| [17] | Avdeev D B, Avdeeva A D. A rigorous three-dimensional magnetotelluric inversion. Progress in Electromagnetics Research , 2006, 62: 41-48. DOI:10.2528/PIER06041205 |
| [18] | 胡祖志, 胡祥云, 何展翔. 大地电磁非线性共轭梯度拟三维反演. 地球物理学报 , 2006, 49(4): 1226–1234. Hu Z Z, Hu X Y, He Z X. Pseudo-three-dimensional magnetotelluric inversion using nonlinear conjugate gradients. Chinese J. Geophys. (in Chinese) , 2006, 49(4): 1226-1234. |
| [19] | Lin C H, Tan H D, Tong T. Three-dimensional conjugate gradient inversion of magnetotelluric sounding data. Applied Geophysics , 2008, 5(4): 314-321. DOI:10.1007/s11770-008-0043-1 |
| [20] | Lin C H, Tan H D, Tong T. Parallel rapid relaxation inversion of 3D magnetotelluric data. Applied Geophysics , 2009, 6(1): 77-83. DOI:10.1007/s11770-009-0010-5 |
| [21] | Mackie R L, Madden T R. Three-dimensional magnetotelluric inversion using conjugate gradients. Geophys. J. Int. , 1993, 115: 215-229. DOI:10.1111/gji.1993.115.issue-1 |
| [22] | Newman G A, Alumbaugh D L. Three-dimensional magnetotelluric inversion using non-linear conjugate gradients. Geophys. J. Int , 2000, 140: 410-424. DOI:10.1046/j.1365-246x.2000.00007.x |
| [23] | Spichak V, Popova. Artificial neural network inversion of magnetotelluric data in terms of three-dimensional earth macroparameters. Geophys. J. Int , 2000, 142: 15-26. DOI:10.1046/j.1365-246x.2000.00065.x |
| [24] | Siripunvaraporn W, Egbert G, Lenbury Y, et al. Three-dimensional magnetotelluric inversion: data-space method. Physics of the Earth and Planetary Interiors , 2005, 150(1-3): 3-14. DOI:10.1016/j.pepi.2004.08.023 |
| [25] | Smith J T, Booker J R. Rapid inversion of two- and three-dimensional magnetotelluric data. J. Geophys. Res , 1991, 96: 3905-3922. DOI:10.1029/90JB02416 |
| [26] | 谭捍东, 余钦范, JohnBooker, 等. 大地电磁法三维快速松弛反演. 地球物理学报 , 2003, 46(6): 850–855. Tan H D, Yu Q F, Booker J, et al. Three-dimensional rapid relaxation inversion for the magnetotelluric method. Chinese J. Geophys. (in Chinese) , 2003, 46(6): 850-855. |
| [27] | Zhdanov M S, Fang S, Hursan G. Electromagnetic inversion using quasi-linear approximation. Geophysics , 2000, 65(5): 1501-1513. DOI:10.1190/1.1444839 |
| [28] | 林昌洪, 谭捍东, 佟拓. 利用大地电磁三维反演方法获得二维剖面附近三维电阻率结构的可行性. 地球物理学报 , 2011, 54(1): 245–256. Lin C H, Tan H D, Tong T. The possibility of obtaining nearby 3D resistivity structure from magnetotelluric 2D profile data using 3D inversion. Chinese J. Geophys. (in Chinese) , 2011, 54(1): 245-256. |
| [29] | 谭捍东, 魏文博, 邓明, 等. 大地电磁法张量阻抗通用计算公式. 石油地球物理勘探 , 2004, 39(1): 114–116. Tan H D, Wei W B, Deng M, et al. General use formula in MT tensor impedance. Oil Geophysical Prospecting (in Chinese) , 2004, 39(1): 114-116. |
| [30] | 谭捍东, 余钦范, JohnBooker, 等. 大地电磁法三维交错采样有限差分数值模拟. 地球物理学报 , 2003, 46(5): 705–711. Tan H D, Yu Q F, Booker J, et al. Magnetotelluric three-dimensional modeling using the staggered-grid finite difference method. Chinese J. Geophys. (in Chinese) , 2003, 46(5): 705-711. |