2. 中国地质调查局发展研究中心, 北京 100037;
3. 黑龙江省水文地质工程地质勘察院, 哈尔滨 150036
2. Development and Research Center, China Geological Survey, Beijing 100037, China;
3. Hydrogeology and Engineering Geology Prospecting Institute of Heilongjiang Province, Harbin 150036, China
大地电磁测深(Simpson and Bahr, 2005; Chave and Jones, 2012)是一种被动源电磁测深方法.它利用在地表测量随着时间变化的电场和磁场来推断地下电阻率分布的电磁方法,在油气勘探(Vozoff,1972; 王鑫等,2010)、矿产勘查(Livelybrooks et al., 1996; Saraev et al., 2010;汤井田等,2012),地热勘查(Stanley et al., 1977; Mogi and Nakama, 1993)和地壳深部探测中(Patro et al., 2005; Kelbert et al., 2008; 张乐天等,2012)得到广泛应用.在三维地质调查中,相对重力、磁力方法,大地电磁(MT)在深部构造探测上起到关键作用.MT已经从一维走向二维,三维正在发展之中.由于受计算机内存,运行时间,数据量庞大的制约,三维MT反演没有达到普遍应用(魏文博,2002;胡祖志和胡祥云,2005;陈向斌等,2011; Siripunvaraporn,2012; Miensopust et al., 2013).尽管二维MT达到普遍应用,但是在资料处理与反演中,还存在一些问题没有得到很好解决.例如二维MT反演的数据旋转方向选择(陈小斌,2008;赵维俊和孙中任,2013),数据极化模式选择(蔡军涛和陈小斌,2010;贺春艳和郭秋峰,2009)等.二维MT常见的三种反演方法(陈向斌等,2011;韩波等,2012)是非线性共轭梯度法(NLCG)(Rodi and Mackie, 2001)、奥卡姆法(Occam)法(Constable et al., 1987; deGroot-Hedin and Constable, 1990;吴小平和徐果明,1998)和数据子空间法(Siripunvaraporn and Egbert, 2000).这些方法都有各自优缺点,但是由于NLCG方法速度快,内存占用少而得到广泛应用.该方法已经应用在三维MT反演上(Newman and Alumbaugh, 2000; 胡祖志等,2006;张昆等,2013).该方法理论介绍比较详细,然而在实践中如何反演MT数据、如何进行参数设置的介绍却很少.所以利用实测MT数据,使用不同参数来进行NLCG反演试验对于大地电磁实践者具有一定指导意义.
本文首先介绍二维NLCG反演方法的原理.尽管Rodi与Mackie做了详细介绍,但是有必要在文章中阐述一下,以便读者更好理解后面的参数试验.在中国地质调查局项目“松辽盆地外围深部地质调查”支持下,2012年,我们在松辽盆地西部的扎鲁特盆地采集5条大地电磁剖面.选择其中D-D’剖面,进行初始模型电阻率,正则化参数,底板误差,反演模式选择试验.不同的参数选择,对反演有一定的影响.在特定参数下,个别地方反演结果会造成很大差别,容易形成错误解释.特别在反演模式选择上,要比较分析反演结果.否则容易造成错误的地质解释.
1 NLCG反演方法逆问题通常可以表示为



因为大地电磁反演具有非唯一性,所以很多模型都能够满足拟合条件.为了从很多模型中选取一个合适模型,所以需要惩罚函数-模型粗糙度(roughness)来约束反演.


在松辽盆地外围深部地质调查项目支持下,2012年,我们在内蒙古扎鲁特盆地采集5条NW向MT剖面.剖面测量方向是从西北到东南.MT设计点距为1 km,有效频率带宽为0.01 Hz到320 Hz.采用加拿大凤凰公司V8 MTU-A和RXU大地电磁记录仪采集五分量数据.电场X正方向与测线方向一致.选择质量好、具有代表性的106 km长D-D’剖面做二维MT反演试验.该剖面TM模式视电阻率和阻抗相位在图 1中显示,而TE模式视电阻率和阻抗相位在图 2中显示.注意测点160由于质量问题排除在外.从TM和TE电阻率拟断面图可以看出,西北向电阻率相当高,符合地质上火山岩与侵入岩电阻率特征.中部进入扎鲁特盆地西部边缘,显示大块低阻区.东南方向,浅部是第四系覆盖区,显示为低阻区,中深部显示电阻率为中高阻区.
2.1 初始模型电阻率试验在二维大地电磁反演之前,首先建立一个电阻率网格模型.初始模型一般是建立一个恒定电阻率的半空间网格.非线性共轭梯度法是局部最优化方法,受到初始电阻率模型的影响.纵观图 1和2中两种模式视电阻率从零点几欧姆到一万多欧姆,范围比较大.设置电阻率设为50 Ωm,网格数为61×452;设置电阻率设为100 Ωm,网格数为55×289. 设置电阻率设为200 Ωm,网格数为51×247.因为网格大小主要是由初始模型电阻率和参加反演的频率范围决定的.在初始电阻率模型网格设置后,其他参数保持默认,进行TM+TE模式的共轭梯度联合反演后,得到图 3的反演结果.从整体上看,三个不同电阻率初始模型的反演结果基本相同.然而,在X=22到25 km,褶皱在三个不同初始电阻率模型下形态不同的.初始电阻率50 Ωm模型的反演结果比较细致,而初始电阻率为100和200 Ωm的反演模型比较粗糙.此外,在x=111到120 km,初始模型为50 Ωm的结果显示高阻地质体顶界面更深一些,界线相对清晰一些.当然电阻率由于涉及内存存贮和计算时间限制不能设置非常小.尽管视电阻率范围跨度大,但是初始模型电阻率不能设置更大,因为过大会失去很多细节.
![]() | 图 1 大地电磁剖面D-D’TM模式视电阻率(上)和阻抗相位(下)拟断面图. 视电阻率上面4-430表示剖面D-D’上430号测点 Fig. 1 Apparent resistivity(upper) and impedance phase pseudo-sections of TM mode for MT profile D-D’.The label 4-430 denotes MT station 430 |
![]() | 图 2 大地电磁剖面D-D’TE模式视电阻率(上)和阻抗相位(下)拟断面图. 视电阻率上面4-430表示剖面D-D’上430号测点 Fig. 2 Apparent resistivity(upper) and impedance phase pseudo-sections of TE mode for MT profile D-D’.The label 4-430 denotes MT station 430 |
![]() | 图 3 大地电磁D-D’剖面 TM+TE模式二维联合反演电阻率网格模型.所有反演的初始模型为恒定电阻率 均匀半空间.图中(a),(b)与(c)反演结果使用的初始模型电阻率分别为50,100和200 Fig. 3 2-D joint TM+TE mode inversion resisitivity models for MT profile D-D’.All the inversion initial models are resistivity-constant half-spaces. The initial model resistivities in(a),(b) and (c)are 50,100 and 200,respectively |
在NLCG反演中,正则化参数λ(见方程5)作用是平衡模型拟合程度与光滑程度的参数.λ参数如果过大,模型拟合程度很低,也就是均方根误差(RMS)很大.反之,模型程度拟合很高,但是模型光滑程度很低.使用不同的λ值,其它参数保持默认值,TM+TE 模式联合反演的电阻率网格模型在图 4中显示.图 4中a,b,c,和d的值分别为3,6,12,24. 结果显示整体电阻率分布是非常相似的.西北部大块高阻岩石,间断地被低阻断层隔断,中部低阻区,进入扎鲁特盆地西部边缘,然后出现大块侵入高阻岩体.然而值得注意,在x=24~31 km区间,随着λ值增大,断层范围越来越大,边界越来越模糊.当λ等于3和6时,位于x=58 km处断层能够清晰辨别,而当λ等于12和24时,基本不能识别.此外,在x=112~120 km,随着λ值增大,高阻体电阻率变得越来越低,与上覆沉积地层电阻率越来越接近.λ的选取应尽量使反演的均方根(RMS)误差小,然而噪音和三维数据影响反演达到很小的RMS误差,例如1.5.可以适当降低λ值来减小RMS误差.但是RMS误差越小,并不等于反演效果越好.在MT数据反演中,λ值在3到300之间是很典型的(WinGlink User Guide,2011).所以对D-D’剖面,λ取值为6是一个合理的选择.
![]() | 图 4 应用不同正则化参数,大地电磁剖面D-D’TM+TE模式二维共轭梯度法反演电阻率网格图比较. 图中(a),(b),(c)与(d)中使用的正则化参数分别为3,6,12和24. Fig. 4 2-D joint TM+TE NLCG inversion resistivity model comparison with different regularization parameters for MT profile D-D’. The regularization parameters in(a),(b),(c) and (d)are 3,6,12 and 24,respectively. |
![]() | 图 5 大地电磁D-D’剖面在不同底板误差参数设置下的TM+TE模式联合二维反演电阻率网格模型.除底板误差 外,其它反演参数保持相同.视电阻率对数误差底板在(a),(b),(c)和(d)中分别为为0.5%,1%,2%和3%. Fig. 5 2-D joint TM+TE mode inversion resistivity models with distinct error floors for MT profile D-D’. Apart from the error floor,other parametes remain unchanged. Error floors in(a),(b),(c) and (d)are 0.5%,1%,2% and 3%,respectively. |
采用二维NLCG反演法,大地电磁剖面D-D’的TM模式,TE模式和TM+TE模式反演结果分别在图 6a,b和c中显示.这三种反演迭代80次,得到的RMS误差分别为2.98%,6.06%与5.16%.可见TE模式拟合误差最大.TE+TM联合模式处于中间.TM模式a反演结果显示剖面西北部有零散的深部高阻体出现,中间到东南部都是大片低阻体,除在x=88 km和x=96 km处存在两个高阻体之外.此反演结果与TM视电阻率拟断面图(图 1)非常吻合.x=60 m附近是扎鲁特盆地西部边缘.从此处向东南部进入扎鲁特盆地腹部.TE模式b反演结果与TM差别很大的.在X=14~56 km区间,出现大块连续高阻体.在中部及东南部,盆地内沉积岩电阻率比TM反演的要高一些.TM+TE模式联合反演结果c综合TM和TE单独反演的特征.明显特征是在x=30,48和58 km处的断层清晰可见.在x=92 km处,TM和TE单独反演中,似断层低阻体延伸超出海拔-4 km,而在TM+TE联合反演中,这个低阻体延伸到海拔-1 km多.在地质解释中,主要参考联合反演外,也应该参考TM和TE单独反演,因为有时候地下结构违反二维的假设.
![]() | 图 6 大地电磁D-D’剖面TM模式(a),TE模式(b)和TM+TE模式(c)二维NLCG法 反演后电阻率网格模型.反演参数保持相同. Fig. 6 TM mode(a),TE mode(b) and TM+TE mode(c)NLCG inversion resistivity models. Inversion parameters remain the same. |
为了增强TM,TE,TM+TE电阻率剖面地质解释的信心,我们从图 7到9列出三组视电阻率和阻抗相位曲线图.从测站340到350(x=49到50 km),红蓝曲线位置相对调换,明显具有断层的特征.TM,TM+TE模式反演结果中出现断层,而TE反演结果没有出现.同理,测站460到470(x=59到61 km),曲线位置颠倒,也是具有断层特征.TE,TM+TE模式反演结果中出现此断层,而TM反演结果没有体现.然而测站760到790,红蓝曲线没有整个对调,只有在大于100 Hz的范围内,部分调换,所以说明这个不太可能是深大断裂,只有浅部电阻率变化.所以这段既不符合TM模式反演结果,也不符合TE模式反演结果,符合TM+TE模式反演结果.
![]() | 图 7 连续4个测站330,340,350,360的视电阻率(上)与阻抗相位(下)曲线图.红色代表yx模式,蓝色代表xy模式 Fig. 7 Apparent resistivity(upper) and impedance phase(lower)plots of contiguous four MT stations 330,340, 350 and 360. The red denotes the YX mode,while the blue denotes the XY mode |
![]() | 图 8 连续4个测站440,450,460,470的视电阻率(上)与阻抗相位(下)曲线图.红色代表yx模式,蓝色代表xy模式 Fig. 8 Apparent resistivity(upper) and impedance phase(lower)plots of contiguous four MT stations 440,450, 460 and 470. The red denotes the YX mode,while the blue denotes the XY mode |
![]() | 图 9 连续4个测站760,770,780,790的视电阻率(上)与阻抗相位(下)曲线图.红色代表yx模式,蓝色代表xy模式 Fig. 9 Apparent resistivity(upper) and impedance phase(lower)plots of contiguous four measurement stations 760, 770,780 and 790. The red denotes the YX mode,while the blue denotes the XY mode |
3.1 本文采用内蒙古扎鲁特盆地采集的一条大地电磁剖面数据为例,对初始模型电阻率,正则化参数,误差底板做了一系列试验.在通常取值范围内,不同参数选择,对反演结果有一定影响,但是影响不大.在反演模式选择上,一定要谨慎.TM和TE模式反演结果差异很大,分辨的细节不同.TM+TE模式联合反演结果综合TM与TE单独反演后特征,相对更准确表达地下结构.在地质解释时,除了主要参考联合反演外,也应该参考独立反演.
3.2 在地质解释时,应该参考视电阻率和阻抗相位的曲线图、拟断面图.因为地下结构不能保证在整个剖面是二维的,有些地方是三维的,而且有的测站噪音大,所以二维反演的地下结构不一定准确.另一方面,二维MT反演的解是非唯一的,应该参考地面、井中地球物理资料、地质资料等信息进行地质解释.
致 谢 衷心感谢R and y L. Mackie先生和谭悍东教授的关于底板误差的讨论.同时,对在内蒙古扎鲁特盆地采集大地电磁数据的野外工作人员表示衷心感谢.作者感谢中国地质调查局松辽盆地外围深部地质调查项目的支持.第一作者感谢沈阳地质矿产研究所自设项目的支持.作者感谢地球物理进展杂志匿名审稿人提供的审稿意见和建议.| [1] | Cai J T, Chen X B. 2010. Refined techniques for data processing and two-dimensional inversion in magnetotelluricⅡ: Which data polarization mode should be used in 2D inversion[J]. Chinese Journal of Geophysics (in Chinese), 53(11): 2703-2714. |
| [2] | Chave A D, Jones A G. 2012. The Magnetotelluric Method: Theory and Practice[M].Cambridge: Cambridge University Press. |
| [3] | Chen X B, Zhao G Z, Ma X. 2008. Discussion of selection of data roration direction in 2-D magnetotelluric inversion[J]. Oil Geophysical Prospecting (in Chinese), 43(1): 113-118, 128. |
| [4] | Chen X B, Lü Q T, Zhang K. 2011. Review of magnetotelluric data inversion methods[J]. Progress in Geophyiscs (in Chinese), 26(5): 1607-1619. |
| [5] | Constable S C, Parker R L, Constable C G. 1987. occam's inversion: a practical algorithm for generating smooth models from electromagnetic sounding data[J].Geophysics, 52(3): 289-300. |
| [6] | deGroot-Hedin C, Constable S. 1990. Occam's inversion to generate smooth, two-dimensional models from magnetotelluric data[J]. Geophysics, 55(12): 1613-1624. |
| [7] | Egbert G D, Booker J R. 1986. Robust estimation of geomagnetic transfer function[J]. Geophysical Journal International, 87(1): 173-194. |
| [8] | Egbert G D, Livelybrooks D W. 1996. Single station magnetotelluric impedance estimation: Coherence weighting and the regression M-estimate[J]. Geophysics, 61(4): 964-970. |
| [9] | Han B, Hu X Y, He Z X, et al. 2012. Mathematical classification of magnetotelluric inversion methods[J]. Oil Geophysical Prospecting (in Chinese), 47(1): 177-187. |
| [10] | He C Y, Guo Q F. 2009. A study on magnetotelluric polarization mode identification by using WinGLink software[C]/Liu Y S, Wu S G, Wang H H, et al. Geophysics in Shandong province over Six Decades (in Chinese). Qingdao: China Ocean University Press, 557-564. |
| [11] | Hu Z Z, Hu X Y. 2005. Review of three dimensional magnetotelluric inversion methods[J]. Progress in Geophysics (in Chinese), 20(1): 214-220. |
| [12] | Hu Z Z, Hu X Y, He Z X. 2006. Pseudo-Three-Dimensional magnetotelluric inversion using nonlinear conjugate gradients[J]. Chinese Journal of Geophysics (in Chinese), 49(4): 1226-1234. |
| [13] | Kelbert A, Egbert G D, Schultz A. 2008. Non-linear conjugate gradient inversion for global EM induction: resolution studies[J]. Geophysical Journal International, 173(2): 365-381. |
| [14] | Larsen J C, Mackie R L, Manzella A, et al. 1996. Robust smooth magnetotelluric transfer functions[J]. Geophysical Journal International, 124(3): 801-819. |
| [15] | Livelybrooks D, Mareschal M, Blais E, et al. 1996. Magnetotelluric delineation of the Trillabelle massive sulfide body in Sudbury, Ontario[J]. Geophysics, 61(4): 971-986. |
| [16] | Miensopust M P, Queralt P, Jones A G. 2013. Magnetotelluric 3-D inversion—a review of two successful workshops on forward and inversion code testing and comparison[J]. Geophysical Journal International, 193(3): 1216-1238. |
| [17] | Mogi T, Nakama S. 1993. Magnetotelluric interpretation of the geothermal system of the Kuju volcano, southwest Japan[J]. Journal of Volcanology and Geothermal Research, 56(3): 297-308. |
| [18] | Newman G A, Alumbaugh D L. 2000. Three-dimensional magnetotelluric inversion using nonlinear conjugate gradients[J]. Geophysical Journal International, 140(2): 410-424. |
| [19] | Patro B P K, Sastry R S, Rao M, et al. 2005. Electrical imaging of Narmada-Son Lineament Zone, Centreal India from magnetotellurics[J]. Physics of the Earth and Planetary Interiors, 148(2-4): 215-232. |
| [20] | Rodi W, Mackie R L. 2001. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion[J]. Geophysics, 66(1): 174-187. |
| [21] | Saraev A K, Antaschuk K M, Nikiforov A B, et al. 2010. Audiomagnetotelluric soundings for the diamond exploration[J]. Chinese Journal of Geophysics, 53(3): 657-676. |
| [22] | Simpson F, Bahr K. 2005. Practical Magnetotellurics[M]. Cambridge: Cambridge University Press. |
| [23] | Siripunvaraporn S. 2012. Three-dimensional magnetotelluric inversion: An introductory guide for developers and users[J].Surveys in Geophysics, 33(1): 5-27. |
| [24] | Siripunvaraporn W, Egbert G. 2000. An efficient data-subspace inversion method for 2-D magnetotelluric data[J]. Geophysics, 65(3): 791-803. |
| [25] | Stanley W D, Boehl J E, Bostick F X, et al. 1977. Geothermal significance of magnetotelluric sounding in the Eastern Snake River Plain-Yellowstone Region[J]. Journal of Geophysical Research, 82(17): 2501-2514. |
| [26] | Tang J T, Xu Z M, Xiao X, et al. 2012. Effect rules of strong noise on magnetotelluric (MT) sounding in the Luzong ore cluster area[J]. Chinese Journal of Geophysics (in Chinese), 55(12): 4147-4159. |
| [27] | Vozoff K. 1972. The magnetotelluric method in the exploration of sedimentary basins[J]. Geophysics, 37(1): 98-141. |
| [28] | Wang X, Shan Y, Zhao G Z, et al. 2010. Deep electric structure beneath the northern section of the western margin of the Ordos basin[J]. Chinese Journal of Geophysics (in Chinese), 53(3): 595-604. |
| [29] | Wei W B. 2002. New advance and prospect of magnetotelluric sounding (MT) in China[J]. Progress in Geophysics (in Chinese), 17(2): 245-254. |
| [30] | WinGlink User Guide[M]. 2011.GEOSYSTEM SRL. |
| [31] | Wu X P, Xu G M. 1998. Improvement of OCCAM's inversion for MT data[J]. Chinese Journal of Geophysics (in Chinese), 41(4): 547-554. |
| [32] | Zhang K, Dong H, Yan J Y, et al. 2013. A NLCG inversion method of magnetotellurics with parallel structure[J]. Chinese Journal of Geophysics (in Chinese), 56(11): 3922-3931. |
| [33] | Zhang L T, Jin S, Wei W B, et al. 2012. Electrical structure of crust and upper mantle beneath the eastern margin of the Tibetan plateau and the Sichuan basin[J]. Chinese Journal of Geophysics (in Chinese), 55(12): 4126-4137. |
| [34] | Zhao W J, Sun Z R. 2013. A comparative study of magnetotelluric impedance tensor rotation and curve smoothing methods[J]. Geophysical and Geochemical Exploration (in Chinese), 37(6): 1125-1132. |
| [35] | 蔡军涛, 陈小斌. 2010. 大地电磁资料精细处理和二维反演解释技术研究(二)—反演数据极化模式选择[J]. 地球物理学报, 53(11): 2703-2714. |
| [36] | 陈小斌, 赵国泽, 马宵. 2008. 关于MT二维反演中数据旋转方向的选择问题初探[J]. 石油地球物理勘探, 43(1): 113-118, 128. |
| [37] | 陈向斌, 吕庆田, 张昆. 2011. 大地电磁测深反演方法现状与评述[J]. 地球物理学进展, 26(5): 1607-1619. |
| [38] | 韩波, 胡祥云, 何展翔,等. 2012. 大地电磁反演方法的数学分类[J]. 石油地球物理勘探, 47(1): 177-187. |
| [39] | 贺春艳, 郭秋峰. 2009. 应用Winglink进行大地电磁测深极化模式识别的研究[C]//刘元生, 吴时国, 王怀洪等. 山东地球物理六十年. 青岛: 中国海洋大学出版社, 557-564 |
| [40] | 胡祖志, 胡祥云. 2005. 大地电磁三维反演方法综述[J]. 地球物理学进展, 20(1): 214-220. |
| [41] | 胡祖志, 胡祥云, 何展翔. 2006. 大地电磁非线性共轭梯度拟三维反演[J]. 地球物理学报, 49(4): 1226-1234. |
| [42] | 汤井田, 徐志敏, 肖晓,等. 2012. 庐枞矿集区大地电磁测深强噪声的影响规律[J]. 地球物理学报, 55(12): 4147-4159. |
| [43] | 王鑫, 詹艳, 赵国泽,等. 2010. 鄂尔多斯盆地西缘构造带北段深部电性结构[J]. 地球物理学报, 53(3): 595-604. |
| [44] | 魏文博. 2002. 我国大地电磁测深新进展及瞻望[J]. 地球物理学进展, 17(2): 245-254. |
| [45] | 吴小平, 徐果明. 1998. 大地电磁数据的Occam反演改进[J]. 地球物理学报, 41(4): 547-554. |
| [46] | 张昆, 董浩, 严加永,等. 2013. 一种并行的大地电磁场非线性共轭梯度三维反演方法[J]. 地球物理学报, 56(11): 3922-3931. |
| [47] | 张乐天, 金胜, 魏文博,等. 2012. 青藏高原东缘及四川盆地的壳幔导电性结构研究[J]. 地球物理学报, 55(12): 4126-4137. |
| [48] | 赵维俊, 孙中任. 2013. 大地电磁阻抗张量旋转方法和曲线圆滑方法的比较[J]. 物探与化探, 37(6): 1125-1132. |
2014, Vol. 29










