2011, Vol. 54
2. 长安大学地质工程与测绘学院, 西安 710054
2. College of Geology Engineering and Geomatics, Chang'an University,Xi'an 710054, China
可控源音频大地电磁法(CSAMT)于20 世纪70 年代提出并被应用到地球物理领域，在矿产普查，油气勘探，水文环境等各个方面发挥了巨大的作用.该方法最大的特点是采用人工场源，大大增加了电磁信号强度，弥补了天然场源信号微弱、不易观测等缺点.此外，它还具有工作效率高，利用一个偶极发射，可以在两侧很大的扇形区域内测量；勘探深度范围大；水平和垂直分辨率较高，高阻屏蔽作用小等优点[1].由于该方法的突出优点，近几年，可控源音频大地电磁测深法在南水北调西线工程深埋藏隧洞进行岩性分类、探测断层深部发育情况、指导深部钻孔布置等方面发挥了较大作用[2].
CSAMT 方法虽然有着不可比拟的优点，但是由于场源的存在，也有其固有的不足.可控源大地电磁法由于场源的影响，其理论要复杂很多，问题也最多，如场源附加效应，近区效应，场源阴影效应，过渡带效应等等.这些可以说是该方法本身固有的缺陷，目前的解释方法和技术水平很难对这些问题有一个很好的解释.尤其是对于南水北调西线工程这样沟谷深切、断崖纵横、地质构造复杂、地层陡倾角的高原地区，依据视电阻率拟断面图解决1000 m 以内的构造、岩溶、岩层划分等问题，必然存在一些异常的识别技术难题，比如异常的位置、埋深、边界、延伸、高低阻异常的差异、假异常的识别等，对于这些问题如果没有具体的认识，仅仅依据视电阻率拟断面图来进行地质解译，必然产生较大的误差甚至是严重的错误.因此，根据南水北调西线工程岩性的地球物理特征，设置了不同的三维模型，并对其进行数值计算分析，从三维空间中模拟场的规律，探索了不同地质异常体的特征，为提高可控源音频大地电磁法在南水北调西线工程深埋藏隧洞探测地质异常体的精度是十分必要的.
当前对于可控源电磁法的研究主要集中于一维层状地质体的模拟或沿走向无限延伸地质体的二维近似模拟，但严格来讲，地球物理电磁场问题都应该在三维空间中进行讨论[3].随着计算机硬件技术的不断发展，三维可控源电磁法正演逐渐变得可行.自20世纪70年代，就相继有许多学者对有源电磁法的数值模拟进行了研究.K.H.Lee 和H.F.Morrison[4]研究了二维地电结构、任意电导率分布的电磁散射的有限元数值解，并对结果进行了比较.Eugene A. Badea, Mark E. Everettz等[5]从基于库伦规范下的矢量势出发，采用二次场算法，针对井中回线源下的电磁响应进行了模拟.Yuji Mitsuhata[6~8]对频域CSEM 数据采用2.5 维方法进行了正演和反演的数值计算.国内底青云、王若等[9， 10]就二维线源和2.5维可控源电磁法进行了有限元数值模拟，闫述等[11]采用矢量有限元模拟了电偶源频率电磁测深，但其模型相对简单.Li等[12， 13]提出了基于后验误差估计的自适应有限元算法进行二维海洋可控源电磁法模拟，节省了计算成本，提高了数值精度.三维可控源因为源的加入使其理论和数值模拟非常复杂，它不像大地电磁研究的是平面波，只需在研究区域的上边界赋予电场的一个分量为常数，按照平面波的传播规律，忽略异常体的影响，在下边界施加一维吸收边界条件即可；也不像瞬变电磁方法研究的是纯二次场，根据场的因果关系，只需把初始时刻均匀大地表面的解析解赋予地表以上一个步长的网格.三维可控源数值模拟必须既要考虑空气区域，又要考虑大地区域，这将大大增加内存消耗和计算成本.此外，如何处理源的奇异性问题和外边界条件的施加都是三维可控源电磁法需要解决的难题.
2 CSAMT 法有限元三维数值模拟基本理论设大地为分区均匀、线性、各向同性、非色散介质的导电介质，假设角频率为ω，时间因子为e-iωt，那么在准静态近似下，忽略位移电流，在电偶极子源情况下频域麦克斯韦方程组如下:
|
(1) |
|
(2) |
其中，E为电场强度，H为磁场强度，ω 为角频率，μ为磁导率，σ 为电导率，Js 为电偶极子源电流密度.
将(1)式两边取散度后，再把(2)代入(1)式可导出电场E所满足的矢量波动方程:
|
(3) |
其中k为波数，k2 =iωμσ，方程(3)即为三维可控源电磁法中电场所满足的矢量波动方程，虽然和大地电磁的控制微分方程只差一个源项，但源的加载使得可控源电磁法的理论和数值模拟比大地电磁法复杂许多.
去掉源项，由广义变分原理可得到(3)式的泛函:
|
(4) |
其中，k是准静态近似条件下的波数，
南水北调西线工程引水线路地球物理勘查的主要任务是对引水隧洞附近地层结构和物性特征进行岩性分类和分组、基本查明引水线路区内主要构造线(尤其是断层)的宽度和走向，并对断层构造的赋水性进行初步评价.因此，在方法选择上，要结合西线工程的地理环境和地质条件，一方面探测深度足够大，探测效率足够高，另一方面探测效果足够好.基于以上原因，确定了以可控源音频大地电磁法为主，配合其他物探方法的技术思路.
探测区地层基本为2层结构:覆盖层、基岩.覆盖层主要为第四纪松散层，电阻率在几十欧姆米.基岩主要为砂岩、板岩或砂板岩互层，砂岩电阻率为800~5000Ωm, 板岩为20~500Ωm, 砂岩板岩互层为400~2000Ωm, 板岩夹砂岩为100~500Ωm.板岩表现为低阻，砂岩表现为相对高阻.断层或破碎带的电阻率则为300Ωm 以下，表现为相对低阻.
3.2 模型的设置及模拟结果根据南水北调西线工程岩性的特征及CSAMT法探测的目的，分别设置了以下7个模型.
3.2.1 单个低阻体异常在测区设计一个低阻体，该低阻体沿x方向长200m, 沿y方向长200m, 沿z方向长200m, 距地面分别为100m、300m、500m、800m.低阻异常体的电阻率为10Ωm、50Ωm、100Ωm, 背景电阻率为2000Ωm.模型结构如图 1.
|
图 1 模型结构图 (a)平面图；(b)剖面图. Fig. 1 Model structural drawing (a)Plain； (b)Profile. |
测线沿x方向，主测线距电偶源5100 m, 其接收点x坐标如下:
-750，-650，-550，-450，-350，-250，-150，-50，50，150，250，350，450，550，650，750.
所采用的频率为:256，512，1024，2048，4096，8192Hz.
图 2 是埋深d=100 m, 异常电阻率ρ异常=10Ωm, 背景电阻率ρ1=10Ωm 计算结果.
|
图 2 埋深d=100m, ρ异常=10Ωm视电阻率拟断面图 Fig. 2 d=100m, ρa=10Ωm resistivity drafting cross section |
图 3(a~c)是埋深d=300m, 异常电阻率分别为ρ异常=10 Ωm、50 Ωm、100 Ωm, 背景电阻率为ρ1=2000Ωm时的计算结果.从拟断面图可以看出，在埋深d=300 m 时，对不同低阻异常的响应反映都比较明显，但是低阻异常的值越高，视深度与异常体的埋深相比越偏深.异常体边界等值线分布密，说明异常体和围岩之间发生了电阻率突变，由此可以确定低阻异常体边界位置.
|
图 3 视电阻率拟断面图 (a)(b)(c)埋深d=300m；(d)(e)(f)d=500m ；(a)(d)ρ异常=10Ωm；(b)(e)ρ异常=50Ωm；(c)(f)ρ异常=100Ωm. Fig. 3 Resistivity drafting cross section |
图 3(d~f)是埋深d=500m, 异常电阻率分别ρ异常=10Ωm、50Ωm、100Ωm, 背景电阻率为ρ1=2000Ωm 时的计算结果.从拟断面图可以看出，在埋深d=500m 时，对不同低阻异常的响应反映较明显，仍可以确定低阻异常体边界位置，但是在低阻体两旁存在明显的由边界效应引起的对称假异常.
3.2.2 单个高阻体在测区设计一个高阻体，该高阻体沿x方向长200m, 沿y方向长200m, 沿z方向长200m, 距地面分别为100m、300m、500m、800m.高阻异常体的电阻率为3000Ωm, 背景电阻率为1000Ωm.进行数值模拟.
图 4是埋深d=100 m, 异常电阻率为ρ异常=3000Ωm, 背景电阻率为ρ1=1000Ωm时的计算结果.从图 4的视电阻率拟断面图中可以看到，高阻体的异常范围要小于低阻体异常范围，而且在异常体边界位置视电阻率等值线变化不像低阻体那样剧烈，视电阻率异常反映也稍偏上，两边出现对称的低阻体假异常.
|
图 4 埋深d=100m, ρ异常=3000Ωm视电阻率拟断面图 Fig. 4 d=100m, ρa=3000Ωm resistivity drafting cross section |
图 5(a, b, c)分别为埋深d=300m、500m、800m, ρ异常=3000Ωm 视电阻率拟断面图.从图可以看出，模型对不同埋深的高阻异常体均有反映，但是分辨率较低.
|
图 5 埋深d=300m(a)、d=500m(b)和d=800m(c)，ρa=3000Ωm 视电阻率拟断面图 Fig. 5 d=300m(a)、d=500m(b)和d=800m(c)，ρa=3000Ωm resistivity drafting cross section |
在测区设计两个相同大小的高阻体，相距500m, 高阻体的尺寸为沿x方向200m, 沿y方向200m, 沿z方向200m, 距地面为100m.高阻异常体的电阻率为3000Ωm, 背景场电阻率为500Ωm, 见图 6.
|
图 6 模型结构图 (a)平面图；(b)剖面图. Fig. 6 Model structural drawing (a)Plain；(b)Profile |
图 7a是埋深d=100 m, 异常电阻率为ρ异常=3000Ωm, 背景电阻率为ρ1=500 Ωm 时的计算结果.可以看出，两个高阻异常体的视电阻率响应都反映了出来，形态完全对称且一致，异常宽度和异常体基本相同，位置相符.
|
图 7 埋深d=100m 视电阻率拟断面图 (a)ρa=3000Ωm；(b)ρa1=30Ωm, ρa2=200Ωm；(c)ρa1=30Ωm, ρa2=2000Ωm Fig. 7 Resistivity drafting cross sectionwith d=100m |
在测区设计两个相同大小的低阻体，相距500m, 低阻体的尺寸为沿x方向200m, 沿y方向200m, 沿z方向200m, 距地面为100m.低阻异常体的电阻率为分别为30 Ωm, 200 Ωm, 背景场电阻率为2000Ωm.图 7b是埋深d=100m 时该模型的计算结果.从图 7b可以看出，两个低阻异常呈不对称分布，异常体的电阻率越低，异常宽度和异常深度影响越大，但异常的中心基本反映了异常体的深度，由此可见，可控源电磁法的横向分辨率要远大于纵向分辨率.
3.2.5 横向一个高阻异常体和一个低阻异常体模拟在测区设计两个相同大小的低阻体和高阻体，相距500m, 尺寸均为沿x方向200 m, 沿y方向200m, 沿z方向200m, 距地面为100m.低阻异常体的电阻率为分别为30Ωm, 高阻为2000Ωm, 背景场电阻率为1000 Ωm.图 7c是埋深d=100m时该模型的计算结果.从图 7c可以看出，低阻异常反映非常灵敏，低阻异常宽度和深度影响大，对高阻异常体基本没有反映.
3.2.6 单个倾斜高阻体在测区设计一个高阻体，该高阻体沿x方向长200m, 沿y方向长200m, 沿z方向长600m, 距地面为100m, 倾斜45°.高阻异常体的电阻率为2000Ωm, 背景电阻率为300Ωm, 见图 8，对该模型进行数值模拟.
|
图 8 模型结构图 (a)平面图；(b)剖面图. Fig. 8 Model structural drawing (a)Plain；(b)Profile. |
图 9a是埋深d=100m 时该模型的计算结果.从图 9a可以看出，高阻倾斜异常反映较明显，深部延伸情况不清晰.
|
图 9 埋深d=100m, 倾斜45°，ρ异常=2000Ωm (a)和ρ异常=50Ωm(b)视电阻率拟断面图 Fig. 9 d=100m, tiltangle45°，ρa=2000Ωm (a)andρa=50Ωm (b)resistivity drafting cross section |
在测区设计一个高阻体，该高阻体沿x方向长200m, 沿y方向长200m, 沿z方向长600m, 距地面分别为100 m, 倾斜45°.低阻异常体的电阻率为50Ωm, 背景电阻率为300Ωm.进行数值模拟.
图 9b是埋深d=100m 时该模型的计算结果.从图 9b可以看出，低阻倾斜异常反映非常明显，低阻异常宽度和深度影响大，低阻异常的两边呈现不对称的高阻假异常.
3.3 模型数值模拟结果分析根据3.2节的数值模拟结果，对7个模型模拟结果分析如表 1.
|
表 1 模拟结果分析 Table 1 The analysis of numerical simulation result |
本文根据南水北调西线工程岩性的地球物理特征，设置了7个不同的三维模型，并对其进行了有限元数值计算，探索地质体异常的规律.
对单个的低阻异常的响应反映明显，视电阻率响应宽度和异常体宽度基本相同，而视深度也基本符合异常体的埋深，但是在低阻体两旁会有对称的高阻假异常，在实际中要引起特别的注意；单个高阻体的异常范围要小于低阻体异常范围，而且在异常体边界位置视电阻率等值线变化不像低阻体那样剧烈，两边出现对称的低阻假异常；两个高阻异常体和两个低阻异常体的视电阻率响应都能反映出来，形态完全对称且一致，异常宽度和异常体基本相同，位置相符；对于高阻和低阻异常体组合，低阻异常反映非常灵敏，高阻异常体基本没有反映；对于倾斜异常体，低阻反映非常明显，高阻倾斜异常反映较明显，深部延伸情况不清晰；异常体的埋深和电阻率的大小对分辨率影响较大，埋深越大，分辨率越低，异常的电阻率越高，分辨率越低.
| [1] | 何继善编译. 可控源音频大地电磁法. 长沙: 中南大学出版社, 1990. He J S. Controlled Source Audio-Frequency Magnetotellurics Method . Changsha: Central South University Press, 1990 |
| [2] | 薛云峰, 何继善, 郭玉松. 南水北调西线工程深埋隧道地质超前预报系统研究的思考. 地球物理学进展 , 2006, 21(3): 993–997. Xue Y F, He J S, Guo Y S. Advance geological prediction system research used in deeply buried tunnel project of South-to-North Transfer Water. Progress in Geophysics (in Chinese) , 2006, 21(3): 993-997. |
| [3] | 闫 述. 基于三维有限元数值模拟的电和电磁探测研究. 西安: 西安交通大学, 2003. Yan S. Studies on the electrical and electromagnetic prospecting based on three-dimensional finite element numerical modeling. Xi'an: Xi'an Jiaotong University, 2003 |
| [4] | Lee K H, Morrison H F. A numerical solution for the electromagnetic scattering by a two-dimensional inhomogeneity. Geophysics , 1985, 50(3): 466-472. DOI:10.1190/1.1441924 |
| [5] | Badea E A, Everett M E, Newman G A, et al. Finite-element analysis of controlled-source electromagnetic induction using Coulomb-gauged potentials. Geophysics , 2001, 66(3): 786-799. DOI:10.1190/1.1444968 |
| [6] | Mitsuhata Y, Uchida T. 2.5-D modeling and inversion of CSEM data. Expanded Abstracts, 69th Ann. Int. Mtg., Soc. Expl. Geophys. 1999 |
| [7] | Mitsuhata Y. 2-D electromagnetic modeling by finite-element method with a dipole source and topography. Geophysics , 2000, 65(2): 465-475. DOI:10.1190/1.1444740 |
| [8] | Mitsuhata Y, Uchida T. 2.5-D inversion of frequency-domain CSEM data based on quasi-linearized ABIC. Calgary: 70th SEG Meeting, 2000 |
| [9] | 底青云, UnsworthM, 王妙月. 复杂介质有限元法2.5维可控源音频大地电磁法数值模拟. 地球物理学报 , 2004, 47(4): 723–730. Di Q Y, Unsworth M, Wang M Y. 2.5-D CSAMT modeling with the finite element method over 2-D complex earth media. Chinese J. Geophys. (in Chinese) , 2004, 47(4): 723-730. |
| [10] | 底青云, UnsworthM, 王妙月. 有限元法2.5维CSAMT数值模拟. 地球物理学进展 , 2004, 19(2): 317–324. Di Q Y, Unsworth M, Wang M Y. 2.5 D CSAMT modeling with the finite element method. Progress in Geophysics (in Chinese) , 2004, 19(2): 317-324. |
| [11] | 闫述, 陈明生. 电偶源频率电磁测深三维地电模型有限元正演. 煤田地质与勘探 , 2000, 28(3): 50–56. Yan S, Chen M S. Finite element solution of three-dimensional geoelectric models frequency electromagnetic sounding excited by a horizontal electric dipole. Coal Geology & Exploration (in Chinese) , 2000, 28(3): 50-56. |
| [12] | Li Y G, Key K. 2D marine controlled-source electromagnetic modeling: Part 1—An adaptive finite-element algorithm. Geophysics , 2007, 72(2): 51-62. |
| [13] | Li Y G, Constable S. 2D marine controlled-source electromagnetic modeling: Part 2—The effect of bathymetry. Geophysics , 2007, 72(2): 63-71. DOI:10.1190/1.2430647 |
| [14] | 张继锋, 汤井田, 喻言, 等. 基于电场矢量波动方程的三维可控源电磁法有限单元法数值模拟. 地球物理学报 , 2009, 52(12): 3132–3141. Zhang J F, Tang J T, Yu Y, et al. Three dimensional controlled source electromagnetic numerical simulation based on electric field vector wave equation using finite element method. Chinese J. Geophys. (in Chinese) , 2009, 52(12): 3132-3141. |
| [15] | 徐世浙. 地球物理中的有限单元法. 北京: 科学出版社, 1994 . Xu S Z. The Finite Element Method in Geophysics (in Chinese). Beijing: Science Press, 1994 . |
| [16] | 金建铭. 电磁场有限单元法. 西安: 西安电子科技大学出版社, 1998 . Jin J M. Finite Element Method of Electromagnetic Field (in Chinese). Xi'an: Xi'an Electronic University Press, 1998 . |