2019, Vol. 62
2. 中国科学院地质与地球物理研究所, 北京 100029;
3. 中国科学院页岩气与地质工程重点实验室, 北京 100029;
4. 中国科学院地球科学研究院, 北京 100029
2. Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
3. Key Laboratory of Shale Gas and Geoengineering, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
4. Institute of Earth Science, Chinese Academy of Sciences, Beijing 100029, China
随着我国经济的不断发展，对矿产资源的需求越来越大，而我国目前浅层或者近地表的大型矿床越来越少，需要向深部“第二空间”进军，寻找深部隐伏矿已成为地质勘查工作的重点(李生虎，2017；杨政伟，2017；才立璇，2018；车春林，2018；陈光椿，2018).由于深部矿相对于浅部矿具有埋深大、信息弱和干扰大等特点，传统的矿产勘查技术在探测深度、探测分辨率以及抗干扰等方面需要有较大的改进和提高(曹新志等，2009).
可控源音频大地电磁法(简称CSAMT)作为一种人工源电磁法，具有采集信号强，分辨率高等特点，已逐步应用于金属矿勘查中，并在矿产勘查中取得了较好的效果(Zonge and Hughesh, 1991; Wannamaker, 1997; Unsworth et al., 2000; Takahashi T, 2004; Pellerin and Wannamaker, 2005).但深部矿相较于浅部矿有更厚的地表作为屏蔽层，在地表放置的接收仪器，接收到的有效信息弱，分辨率低，同时容易受到地表强电磁干扰，可能出现虚假的结果，造成勘探经费的浪费(曹新志等，2009；王艳等，2009；祁晓雨，2015；杨政伟，2017).本文尝试提出一种新的观测方式，将接收机放置在巷道中，接收机更加靠近深部目标体，对于深部矿体的分辨率会有所提高，同时在巷道中采集信号，地面产生的人为噪声对接收到的信号影响较小.为此，开展了地面激发巷道接收模式的CSAMT法对目标体分辨能力的研究.
目前，CSAMT三维正演模拟的方法主要包括有限差分，有限元和积分方程法等.有限差分和有限元法在进行正演模拟时需要对整个空间进行离散，而积分方程法只需要对于异常体进行剖分，也不需要处理边界问题，所以运算速度快，占用内存少(底青云等，2006；Tang et al., 2011；穆纳尔丁·托合提等，2016；关振玮等，2016；王若等, 2009, 2014；Key et al., 2011).因此，本文基于积分方程法，首先对简单异常体进行CSAMT正演模拟，分析其计算结果，之后，结合实际地电结构，设计了一个复杂模型，利用积分方程正演程序对其进行CSAMT正演计算，分析了在异常体的埋深不变时不同接收深度的CSAMT频率响应特征，并对比了不同接收深度的结果，同时为了验证二维模拟结果的正确性，利用MARE2DEM软件进行了结果对比.
1 方法原理积分方程法是以非均匀传导介质中的电磁场的不同类型的积分表达式为基础的数值方法.积分方程法和有限差分以及有限元法相比最大的优点是对层状介质背景条件下二维和三维地质体的电磁场响应能够精确而快速地模拟，其基本原理示于图 1.
|
图 1 积分方程法基本原理示意图 Fig. 1 Schematic diagram showing the basic principle of the integral equation method |
图 1为一个背景介质电导率为σb的三维地电模型，存在一个局部异常区域D，电导率分布为σ=σa+σb.在这个模型中，电磁场是由区域Q中的外部电流密度分布je产生.当电磁场为时谐变化时，时谐因子取e－iωt，采用积分方程法求解这个模型中的电磁场问题.
在模拟区域中，电磁场可以表示成正常场(背景场)和异常场(散射场)之和.
|
(1) |
|
(2) |
Eb和Hb分别表示由给定源在背景介质中产生的背景电场和背景磁场，Ea和Ha分别表示由于异常电导率分布σa而产生的异常电场和异常磁场，总场为二者之和.
根据麦克斯韦方程组，场强满足波动方程：
|
(3) |
|
(4) |
引入格林函数，(3)、(4)式等价于：
|
(5) |
根据格林函数，可以推导出基本的三维积分方程公式：
|
(6) |
|
(7) |
这里
|
(8) |
|
(9) |
最后得到了对于异常区域D上的散射电流进行积分的电磁场的表达式
|
(10) |
|
(11) |
其中rj是观测点的位置，如果在异常区域内的电场已知，则可求得在任意观测点rj的电磁场.式(10)和(11)式为积分方程法中三维电磁正演的基本方程.
在积分方程法中，给定源在异常区域产生的总电场强度很难精确求解，因此方程(10)和(11)需要利用近似求解方法对总电场强度进行近似.一般情况下异常场Ea比背景场Eb量级小很多，可以采用背景场代替总场的方法，即Born近似，但是这样求解精度较低.本文研究在巷道接收的问题，计算地下的场强，接收位置有可能穿过异常体，所以本文采用精度更高的连续迭代算法.
本文利用积分方程法可以求出空间中任意位置处的电场和磁场强度，为了更加清楚地分析在不同深度接收的电磁场对于深部异常体的分辨能力，计算了卡尼亚视电阻率：
|
(12) |
式中ω=2πf为角频率，f为频率，μ为磁导率，本文采用真空中的磁导率μ0.
2 数值模拟为了充分了解巷道接收时CSAMT对于异常体的分辨能力.首先设计了地下含有低阻球状异常体的模型，通过对比分析地面和巷道接收时的CSAMT响应特征，说明巷道接收时对异常体的分辨能力；其次将积分方程计算结果同开源软件MARE2DEM的模拟结果进行对比，验证利用积分方程法计算地下含有二维异常体时电磁响应的正确性；最后根据实际地质资料建立相应的地电模型，并进行CSAMT二维正演，分析其电磁响应特征，并将地面和巷道接收的结果进行了对比.
2.1 低阻球状异常体CSAMT正演模拟我国对于深部矿和浅部矿是以地表以下500 m为分界线，埋深小于500 m定义为浅部矿，埋深大于500 m定义为深部矿(杨政伟，2017).为了验证在巷道中接收可以提高CSAMT方法对深部异常体的分辨能力，设计了如图 2所示的低阻球体，异常体埋深为800 m.图 2a中第一层代表泥质层，电阻率ρ1=100 Ωm，层厚度H1=50 m，第二层为背景层，电阻率ρ2=5000 Ωm，在该层中存在一个低阻球体，电阻率ρ0=100 Ωm，球体顶部到地面的距离H=800 m，球体的直径为d=600 m.观测平面图如图 2b所示，发射偶极子长1.2 km，电流大小为20 A，源的中心距测线的最短距离为30 km，在目标体上方设置了11条测线，线距为100 m，每条测线上存在101个测点，点距为10 m.测线的埋深分别设计在地表 0 m、地下200 m、地下400 m进行接收，采用的观测频率为1、2、4、8、16、32、64、128、256 Hz、300、400、512、600、700、800、900、1024、2048、4096 Hz、8192 Hz.采用积分方程法进行正演计算，正演过程中仅对异常体进行网格剖分，大大节省计算量，异常体的网格剖分如图 3所示.
|
图 2 低阻球体模型示意图 (a)模型断面图; (b)观测装置平面图. Fig. 2 Schematic diagrams of a low-resistance sphere model (a) Cross section of the model; (b) Plane view of the observation device. |
|
图 3 异常体的网格剖分图 Fig. 3 Grid subdivision of the anomaly body |
首先，对20个频点进行筛选，选取在观测点处进入远区的频率所对应的场强和视电阻率进行对比.图 4所示为距发射源最近的测点(0，30000，0)处的电场总场、磁场总场及视电阻率的CSAMT频率响应曲线.从图中我们可以发现，在收发距为30 km的情况下，当频率低于32 Hz时，视电阻率曲线不满足远区场特征，因此后文中计算的都是大于32 Hz时的视电阻率.
|
图 4 低阻球体单点Ex、Hy、ρaxy频率响应曲线 (a)电场分量Ex的频率响应；(b)磁场分量Hy的频率响应；(c)视电阻率ρaxy频率响应 Fig. 4 Single-point Ex, Hy, ρaxy frequency response curves of low-resistance sphere (a) Frequency response of the electric field component Ex; (b) Frequency response of magnetic field component Hy; (c) Apparent resistivity ρaxy frequency response. |
图 5为过低阻球体正中间距地面不同深度接收的测线的视电阻率拟断面图，其中测线坐标为y=30500 m.从图中可以看出，在相同频带内，在地面接收的拟断面图上只能反映浅地表的低阻泥质层和深部高阻背景层，高频(2048~8192 Hz)对浅地表的低阻泥质层和低频(64~1024 Hz)对深部高阻背景层，而对于埋深较深的低阻球体则基本没有任何反映.在地下200 m、400 m接收的信号计算的视电阻率拟断面图则能较为明显地看出在高阻背景层中有一个埋深较深的低阻体，同时，我们还可以看出，地下400 m接收时得到的拟断面图中反映出的低阻异常更加明显，以上数值计算结果表明利用CSAMT在巷道中接收相比于地面接收对于深部异常体具有更好的探测效果.
|
图 5 中心测线不同接收深度的视电阻率拟断面图 (a)地面接收视电阻率拟断面图；(b)地下200 m接收视电阻率拟断面图; (c)地下400 m接收视电阻率拟断面图. Fig. 5 Psudo-section contours of apparent resistivity at different receiving depths of the central line (a) Ground surface; (b) 200 m below surface; (c) 400 m below surface. |
为了便于更加直观地分析，对比测线坐标y=30500 m，频率为64 Hz时不同深度接收时电场Ex、磁场Hy及视电阻率ρaxy曲线(图 6).从图 6(a, b)图中可以看出，在地下接收的场强的信号较地表弱，分析原因为地下介质对于电磁场能量的吸收衰减而导致，图 6c为视电阻率曲线，从图中可以看出，在地下200 m和400 m深度接收视电阻率曲线基本重合，同时，对异常体的反映较地表明显.在地表接收的信号受浅地表低阻覆盖层影响，对地下异常体感应的二次场吸收衰减较大，视电阻率表现为整体的体积效应，因此电阻率值在1800 Ωm左右，而地下200 m和400 m深度接收则对异常体感应出的二次场以及背景高阻层有较强的反应，同时由于高阻体对于电磁信号的衰减较弱，因此，地下200 m和400 m接收计算的视电阻率曲线基本重合.为了进一步对结果进行分析，定义了一个差视电阻率δρa，即总视电阻率减去背景场视电阻率.图 6f为差视电阻率曲线，从中可以明显的看出，在地下接收时对于深部异常体的异常响应明显大于地表接收时的异常响应.图 6(d, e)分别为电场Ex异常场、磁场Hy异常场响应曲线，由于在地下接收更加靠近异常体，所以纯异常响应也较地表明显.
|
图 6 低阻球体积分方程正演模拟CSAMT探测曲线 (a)电场总场；(b)磁场总场；(c)视电阻率；(d)电场纯异常；(e)磁场纯异常；(f)差视电阻率. Fig. 6 CSAMT sounding curves from forward modeling on low-resistance sphere by integral equation method (a) Total electric field; (b) Total magnetic field; (c) Apparent resistivity; (d) Pure anomalous electric field; (e) Pure anomalous magnetic field; (f) Differential apparent resistivity. |
为了验证利用积分方程法在计算复杂二维模型CSAMT响应的正确性，设计了如图 7所示的二维模型，发射偶极沿x轴在点(0, 0, 0)处，采用的是单位偶极子发射源，将积分方程法计算的结果同开源软件MARE2DEM的计算结果进行对比.测点设计为沿x轴5 km至7 km处，测点间隔为100 m.计算结果如图 8所示.可以看出两种方法的计算相对误差在3%以内，所以采用积分方程计算二维模型CSAMT响应的结果是可靠的.说明本文采用的积分方程法是可以应用到复杂地质模型的CSAMT数值计算中.
|
图 7 简单二维模型图 Fig. 7 Simple 2D model |
|
图 8 积分方程和MARE2DEM结果对比图 (a)电场对比；(b)电场相对误差；(c)磁场对比；(d)磁场相对误差. Fig. 8 Comparison of results by integral equation and MARE2DEM (a) Contrast of electric field; (b)Electric field relative error; (c) Contrast of magnetic field; (d) Magnetic field relative error. |
实际勘探区的地质条件比较复杂，为了验证巷道接收观测系统在实际野外地质条件下对目标体的探测精度，设计了如图 9所示的二维地质模型，9a为实际的模型图，9b为模型的切面图，电阻率为5000 Ωm的地质体为实际的矿藏，采用积分方程法进行CSAMT正演模拟.观测系统平面图如图 9c所示，发射偶极长为1.2 km，电流为20 A，源的中心距测线的距离为15 km，测线设计为从-5 km到5 km，测点间距为50 m，共201个测点，分别在地面、地下200 m、400 m进行观测接收，采用的观测频率为1、2、4、8、16、32、64、128、256、512、1024 Hz、2048、4096、8196 Hz.
|
图 9 复杂二维模型示意图 (a)模型图; (b)模型截面图; (c)观测装置平面图. Fig. 9 Schematic diagram of complex two-dimensional model (a) Model diagram; (b) Cross section of the model; (c) Plane view of the observation device. |
在收发距为15 km的情况下，当频率低于128 Hz时，此时场进入近场，所以后面视电阻率的计算都是在频率大于128 Hz进行的.图 10为该复杂二维模型正演的在频率为512 Hz时CSAMT响应曲线.从图 10(a, b)中看出，在地下接收的场强比地表弱，埋深每向下200 m，电场场强约降低2×10-7V·m-1, 磁场也是呈现这种规律，但是对于模型的电磁响应比地表更为清晰.图 10c为不同深度接收时频率为512 Hz的视电阻率曲线，可以看出地下400 m接收的曲线很明显地反应了在-1 km到2 km处存在一个高阻侵入体，这与真实模型的地质界限基本上一致.
|
图 10 二维模型积分方程正演模拟CSAMT探测曲线 (a)不同接收深度电场对比图；(b)不同接收深度磁场对比图；(c)不同接收深度视电阻率对比图. Fig. 10 CSAMT sounding curves of 2D forward modeling by integral equation method (a) Comparison of electric field at different receiving depths; (b) Comparison of magnetic field at different receiving depths; (c) Comparison of apparent resistivity at different receiving depths. |
地表接收的电磁信号，在浅部存在一个低阻的覆盖层，所以对于地下异常体感应的二次场的吸收衰减较多，因此视电阻率曲线对于浅地表的低阻体反映比较明显，整体表现为电阻率较低的波动不太大的水平线.而在地下进行接收，由于地下背景层是高阻体，接收点距异常体较近，异常体感应的二次场衰减缓慢，所以其视电阻率曲线能够更加清晰地反映出地下介质的情况.
图 11为二维模型不同接收深度的视电阻率拟断面图，从图中可以看出，在地面接收的信号，在观测的频段内，只能反应出浅地表的低阻层，对于高阻侵入体基本上没有反应，而在地下200 m、地下400 m接收的视电阻率拟断面图则能较为明显地看出有一个高阻侵入体，同时，其侵入位置以及侵入形态也得到了很好的反应.数值计算结果表明CSAMT在巷道中进行接收对于探测断裂以及侵入体相比于在地面接收有更好的效果.
|
图 11 二维模型不同接收深度的视电阻率拟断面图 (a)地面；(b)地下200 m；(c)地下400 m. Fig. 11 Psudo-section contours of forward modeling results of apparent resistivity at different receiving depths (a) Ground; (b) Underground 200 m; (c) Underground 400 m. |
本文利用积分方程法对设计的二维、三维模型进行了CSAMT正演模拟.数值结果表明在巷道中接收电磁信号能够提高对于地下异常体的分辨能力，但在计算时并没有考虑噪声的影响，而在实际工作中，往往会受到噪声的影响，所以在后续的工作中需要研究考虑地面噪声时巷道接收电磁信号对地下目标体的分辨能力.同时本文开展的CSAMT正演模拟可以为反演计算提供理论支持，但要对野外实际指导作用还要进行相关的反演计算，对比不同接收装置时的应用效果.
4 结论通过计算不同异常体模型开展基于积分方程法的CSAMT正演响应，分析在地下不同位置设置接收点的电磁响应特征，可以得到以下结论：
(1) 本文在对二维模型进行CSAMT响应正演模拟时，将积分方程法和MARE2DEM软件的计算结果进行了对比.对比结果表明使用积分方程法进行二维模型CSAMT正演模拟，采用精度较高的迭代算法时，其计算结果是可靠的.
(2) 对于埋深大于500m的简单矿藏，由于地下介质对电磁信号的吸收，在巷道中接收的CSAMT信号强度比地面弱.然而在相同的探测频段内，地面接收电磁场受浅地表覆盖层的影响，其对埋深较大的异常体的分辨能力是小于巷道接收的.
(3) 对于在低阻浅地表覆盖层下的高阻侵入体和地下断裂，从CSAMT正演的结果看出，在巷道中接收，电场和磁场曲线已经可以反映出地下不同地质体的大致形态，其视电阻率曲线和拟断面图则清楚的反映了异常体的位置，而在地表接收只能反映出背景高阻层，因此巷道接收对于探测地下断裂和侵入体的分辨能力更强.
(4) 本文利用积分方程法，采用巷道接收方式对埋深较深(大于500 m)的异常体，以及不同地层组合下的二维模型进行了正演模拟，并对这些模型的CSAMT响应特征进行分析，结果表明在巷道中接收可以提高对深部异常体的分辨力，这对于寻找深部矿藏，覆盖在地表下的侵入体以及断层具有一定的理论指导意义.
致谢 感谢犹他大学(University of Utah) Zhdanov教授CEMI研究组提供的三维积分方程正演程序.
Cai L X. 2018. Physical prospecting methods and their application in deep metal exploration. World Nonferrous Metal (in Chinese), (16): 137-139. |
Cao X Z, Zhang W S, Sun H S. 2009. Progress in the study of deep exploration in China. Geological Science and Technology Information (in Chinese), 28(2): 104-109. |
Che C L. 2018. Geological prospecting methods and problems in deep prospecting. World Nonferrous Metal (in Chinese), (12): 117-119. |
Chen G C. 2018. Effective application of comprehensive geophysical prospecting method in deep prospecting. World Nonferrous Metal (in Chinese), (11): 74-75. |
Di Q Y, Unsworth M, Wang M Y. 2006. 2006.2.5 D finite element CSAMT numerical inversion. . Oil Geophysical Prospecting (in Chinese), 41(1): 100-106. DOI:10.3321/j.issn:1000-7210.2006.01.020 |
Guan Z W, Li T L, Zhu C, et al. 2016. CSAMT simulation of three-dimensional geological model. Progress in Geophysics (in Chinese), 31(3): 1073-1079. DOI:10.6038/pg20160320 |
Key K, Ovall J. 2011. A parallel goal-oriented adaptive finite element method for 2.5-D electromagneticmodelling. Geophysical Journal International, 186(1): 137-154. DOI:10.1111/j.1365-246X.2011.05025.x |
Li S H. 2017. The geological exploration method and the problems of deep ore prospecting are discussed. World Nonferrous Metal (in Chinese), (14): 218-220. |
Munirdin T, Di Q Y, Zhang W W. 2016. CSAMT response of 3D geological structure model using ntegral equation method. Progress in Geophysics (in Chinese), 31(3): 1145-1151. DOI:10.6038/pg20160330 |
Pellerin L, Wannamaker P E. 2005. Multi-dimensional electromagnetic modeling and inversion with application to near-surface earth investigations. Computers and Electronics in Agriculture, 46(1-3): 71-102. DOI:10.1016/j.compag.2004.11.017 |
Qi X Y. 2015. The application and research of CSAMT field work methods. Journal of Railway Engineering Society (in Chinese), 32(6): 15-19. |
Takahashi T. 2004. ISRM Suggested methods for land geophysics in rock engineering. International Journal of Rock Mechanics and Mining Sciences, 41(6): 885-914. DOI:10.1016/j.ijrmms.2004.02.009 |
Tang J, Xiong B. 2018. The simulation of unstructured grid two-dimensional magnetotelluric forward modeling based on COMSOL software. Chinese Journal of Engineering Geophysics (in Chinese), 15(3): 347-356. |
Tang XG, Hu WB, Yan LJ. 2011. Topographic effects on long offset transient electromagnetic response. Applied Geophysics, 8(4): 277-284. DOI:10.1007/s11770-011-0297-x |
Unsworth M J, Lu X Y, Watts M D. 2000. CSAMT exploration at Sellafield:Characterization of a potential radioactive waste disposal site. Geophysics, 65(4): 1070-1079. DOI:10.1190/1.1444800 |
Wang R, Di Q Y, Wang M Y, et al. 2009. Research on the effect of 3D body between transmitter and receivers on CSAMT response using Integral Equation method. Chinese Journal of Geophysics (in Chinese), 52(6): 1573-1582. DOI:10.3969/j.issn.0001-5733.2009.06.019 |
Wang R, Wang M Y, Di Q Y, et al. 2014. 3D1C CSAMT modeling using finite element method. Progress in Geophysics (in Chinese), 29(2): 839-845. DOI:10.6408/pg20140249 |
Wang Y, Lin J, Zhou F D, et al. 2009. Study of resolution and method of CSAMT for the deep low-resistivity. Journal of China University of Mining & Technology (in Chinese), 38(1): 86-90. |
Wannamaker P E. 1997. Tensor CSAMT survey of the Sulphur Springs thermal area, Valles Caldera, New Mexico, PartsⅠ & Ⅱ:impilications for structure of the western caldera and for CSAMT methology. Geophysics, 62: 451-476. DOI:10.1190/1.1444156 |
Yang Z W. 2017. Discussion on deep prospecting method in geological and mineral exploration. World Nonferrous Metal (in Chinese), (10): 93-94. |
Zonge K L, Hughesh L J.1991. Controlled source audio frequency magnetotellurics. //Nabighian M N ed.Electromagnetic Methods in Applied Geophysics.Society for Exploratory Geophysics, 713-908.
|
Society of Exploration Geophysics P. O. Box 702740/Tulsa, Oklahoma.
|
才立璇. 2018. 深部金属矿勘查中常用的物探方法及运用. 世界有色金属, (16): 137-139. DOI:10.3969/j.issn.1002-5065.2018.16.082 |
曹新志, 张旺生, 孙华山. 2009. 我国深部找矿研究进展综述. 地质科技情报, 28(2): 104-109. DOI:10.3969/j.issn.1000-7849.2009.02.019 |
车春林. 2018. 地质勘查方法及深部找矿存在问题. 世界有色金属, (12): 117-119. DOI:10.3969/j.issn.1002-5065.2018.12.065 |
陈光椿. 2018. 综合物探方法在深部找矿中的有效应用. 世界有色金属, (11): 74-75. DOI:10.3969/j.issn.1002-5065.2018.11.047 |
底青云, Unsworth M, 王妙月. 2006. 2.5维有限元法CSAMT数值反演. 石油地球物理勘探, 41(1): 100-106. DOI:10.3321/j.issn:1000-7210.2006.01.020 |
关振玮, 李桐林, 朱成, 等. 2016. 三维地质模型的CSAMT模拟研究. 地球物理学进展, 31(3): 1073-1079. DOI:10.6038/pg20160320 |
李生虎. 2017. 论述地质勘查方法及深部找矿存在问题. 世界有色金属, (14): 218-220. |
穆纳尔丁·托合提, 底青云, 张文伟. 2016. 利用积分方程法研究复杂三维地质结构CSAMT响应特征. 地球物理学进展, 31(3): 1145-1151. DOI:10.6038/pg20160330 |
祁晓雨. 2015. CSAMT野外工作方法的研究与应用. 铁道工程学报, 32(6): 15-19. DOI:10.3969/j.issn.1006-2106.2015.06.004 |
唐杰, 熊彬. 2018. 基于COMSOL软件的非结构化网格下的二维大地电磁正演模拟. 工程地球物理学报, 15(3): 347-356. |
王若, 底青云, 王妙月, 等. 2009. 用积分方程法研究源与勘探区之间的三维体对CSAMT观测曲线的影响. 地球物理学报, 52(6): 1573-1582. DOI:10.3969/j.issn.0001-5733.2009.06.019 |
王若, 王妙月, 底青云, 等. 2014. CSAMT三维单分量有限元正演. 地球物理学进展, 29(2): 839-845. DOI:10.6408/pg20140249 |
王艳, 林君, 周逢道, 等. 2009. CSAMT法深部低阻分辨能力及方法研究. 中国矿业大学学报, 38(1): 86-90. |
杨政伟. 2017. 深部找矿方法在地质矿产勘查的探讨. 世界有色金属, (10): 93-94. |