地球物理学进展  2016, Vol. 31 Issue (3): 943-948   PDF    
引用本文
柴伦炜, 李桐林, 苏晓波, 朱成, 刘永亮, 关振玮. 2016. 尖锐边界反演在临江市六道沟可控源音频大地电磁测深数据反演中的应用[J]. 地球物理学进展, 31(3): 943-948, doi: 10.6038/pg20160302  
CHAI Lun-wei, LI Tong-lin, SU Xiao-bo, ZHU Cheng, LIU Yong-liang, GUAN Zhen-wei. 2016. Sharp boundary inversion used in Linjiang Liudaogou controlled source audio magnetotelluric sounding data processing. Progress in Geophysics, 31(3): 943-948, doi: 10.6038/pg20160302   
尖锐边界反演在临江市六道沟可控源音频大地电磁测深数据反演中的应用
柴伦炜, 李桐林, 苏晓波, 朱成, 刘永亮, 关振玮     
吉林大学 地球探测科学与技术学院, 长春 130026
收稿日期: 2015-12-02 修回日期: 2016-04-14.
基金项目: 中国地质调查局项目吉黑东部综合找矿方法研究(201301-060-004)资助.
作者简介: 柴伦炜,男,1990年生,硕士,主要从事电磁法勘探理论模拟与应用研究.(E-mail:18644988852@163.com)
通讯作者: 李桐林,男,1962年生,教授,博导,从事电磁法勘探理论问题及其应用和固体矿产地球物理综合地质解释等.(E-mail:litl@jlu.edu.cn)
摘要: 在实际大地电磁测深数据反演过程中,为了得到更加清晰的电性单元分界面,传统的Occam平滑反演方法遇到了挑战,本文讨论了sharp boundary inversion(SBI)反演的理论及其应用.通过模型试算,对比Occam和SBI的反演效果,证明了SBI反演划分层状介质的可靠性.本文对临江市六道沟地区实测可控源音频大地电磁测深(CSAMT)数据进行了全区视电阻率转换,然后对全区视电阻率进行Occam反演和SBI反演.最后对比已有地质资料,表明SBI反演能够用于CSAMT数据的处理中.
关键词: 可控源大地电磁     Occam反演     SBI反演    
Sharp boundary inversion used in Linjiang Liudaogou controlled source audio magnetotelluric sounding data processing
CHAI Lun-wei, LI Tong-lin, SU Xiao-bo, ZHU Cheng, LIU Yong-liang, GUAN Zhen-wei     
College for Geoexploration of Science and Technology, Jilin University, Changchun 130026, China
Abstract: In processing of magnetotelluric sounding data inversion,in order to obtain a clearer electrical interface unit, the traditional method of Occam inversion face challenges, this article discusses the sharp boundary inversion (SBI) theory and its applications. By comparing Occam inversion and the SBI inversion, SBI inversion proves to be more reliable in layered media. In this paper, we convert the measured controlled source audio magnetotelluric (CSAMT) data in Linjiang Liudaogou areas to the region apparent resistivity, and then we use region apparent resistivity to do OCCAM inversion and SBI inversion. Finally, comparing the existing geological data, it can show that the SBI inversion can be used to processing CSAMT data.
Key words: controlled source audio magnetotelluric     occam inversion     SBI inversion    
0 引 言

OCCAM反演是具有代表性的光滑模型反演方法，OCCAM反演方法通过使模型粗糙度最小来达到压制虚假的干扰的目的.对于二维构造，Constable等人在反演目标函数中添加了横向和纵向的粗糙度的约束，要使模型粗糙度达到最小，就会导致反演对地下地质体边界位置刻画模糊，从而使得反演界面位置不清晰.为了克服上述缺陷，de Groot-Hedlin和Constable(2004)提出了二维MT的尖锐边界反演(sharp boundary inversion)简称SBI反演方法，参与该反演的参数是各层的电阻率和所有层的层深度，这种方法把地质构造表示为一定数量的层.这些层电阻率是均匀的，层深度在横向是变化的，即不同的测点同一层的层深度不全是相同的.通过引入新的粗糙度限制条件，不仅控制地层电阻率之间的差异，同时也对层深度进行约束，避免了边界位置过多的冗余构造.

Smith等(1999)最先提出了反演二维MT数据的SBI方法，de Groot-Hedlin和Constable(2004)研究了用于反演二维数据的SBI算法.欧东新和王家林(2005)将SBI方法的构建与RRI(快速松弛反演)方法结合，提出了一种新的二维快速反演方法，对SBI反演进行的一定的改善.张罗磊等(2010)对该尖锐边界反演方法做了改进研究，通过引入对角梯度支撑改善了倾斜电性分界面的反演效果.叶益信等(2009)提出了将SBI与RRI结合，提出了先通过RRI反演得到一个光滑模型，再利用模糊聚类方法得到一个适合反演的初始模型，在对初始模型进行SBI反演.

吉林省矿产资源勘探深度平均300 m，大约有三分之一面积的陆地玄武岩覆盖区没有进行勘察，深部资源有巨大的潜力.因此了解长白山玄武岩的厚度以及中古生代沉积盆地的厚度，中古生代沉积盆地基底的起伏形态，对于探测玄武岩覆盖区深部资源十分必要.

本文利用在临江市六道沟多金属矿区内采集的可控源音频大地电磁测深(CSAMT)数据，使用电场正演迭代拟合卡尼亚视电阻率方法对CSAMT数据进行了全区视电阻率转换，获得了全区视电阻率，最后对全区视电阻率进行了Occam反演和SBI反演，比较了两种反演的结果，最后对反演的结果进行分析和解释.

1 理论与方法 1.1 反演理论 1.1.1 SBI反演理论

对于一个nr层的一个均匀半空间，含有nc个小单元块，那么模型的参数是一个长度为(nr+1)+nr×nc的变量，公式为

式中ρi代表第i个电性单元层的电阻率，dij代表第i层的第j个小块的深度.因为上面模型参数变化范围较大，因而对参数取对数值进行反演，以提高反演稳定性.

反演的目标函数为

其中μ-1是拉格朗日因子，d代表大地电磁测深数据，F(m)是通过模型m产生的非线性的正演响应函数，χ2*是期望拟合差，W是一个对角矩阵，其对角值是数据误差的倒数.

SBI反演目标函数与光滑模型反演的目标函数相似，不同之处在于它们对粗糙度的定义，SBI反演的粗糙度R=Sm2，S是一个奖惩矩阵，SBI法不仅对电阻率进行奖惩，而且对相邻地层进行奖惩.下面是一般情况下SBI反演中粗糙度的定义，公式为

其中Srho是一个作用于电阻率参数上的nr×(nr+1)阶矩阵， 式中，
其中c是用来控制是否对地层电阻率进行粗糙度的约束以及在粗糙度矩阵中地层电阻率所占有的权重，如果取常数c=0，表明不对层电阻率进行粗糙度的约束.
其中Si是(nc-1)×nc阶矩阵，表示层深度的粗糙度矩阵，

由于模型正演响应的函数F[m]是非线性的，利用泰勒公式将其线性化，公式为

Ji是灵敏度矩阵.

这样目标函数的迭代形式为

式中，$\hat{d}$=d-F(mi)+Jimi通过迭代形式就可以得到下一个模型修正值.

1.1.2 CSAMT全区视电阻率转换方法

电偶极子在均匀水平大地上的电场和磁场为

式中I0，K0是零阶的第一类和第二类修正贝塞尔函数，I1，K1为一阶的第一类和第二类修正贝塞尔函数，其宗量为$\frac{1}{2}$k1r；r为收发距；波数为电流偶极矩PE=I×AB，变量AB表示供电极距；φ为r方向与x轴正方向的夹角.

视电阻率的运算定义式为

式中μ代表磁导率，单位是H/m；ω为协变电流的圆频率，单位rad/s.本文转换全区视电阻率的程序过程如下：

① 给定一个很小的电阻率的初值ρ1；

② 将ρ1代入上面(1)、(2)、(3)中，求出ρs；

③ 判断ρs-ρ实测/ρs≤ε(ε是给定很小的正数)是否成立，如果成立，程序停止计算，ρs就是转换后的全区视电阻率；否则，ρ1=ρ1+Δρ(Δρ为搜索步长，为给定的值)，返回上一步重新计算，直到ρs-ρ实测/ρs≤ε式成立；

④ 输出转换的全区视电阻率ρs.

1.2 模型反演的效果 1.2.1 模型Ⅰ

如图 1所示为孤立体模型，一个均匀半空间存在一个异常体，异常体电阻率为10 Ω·m，均匀半空间电阻率是100 Ω·m，顶面埋深为3.2 km，大小为6.7 km×10 km，采样点数为9个，采样间距2 km，采集频点有10个，分别为(32，16，8，4，2，1，0.5，0.25，0.125，0.0625)Hz，利用有限元进行正演，正演网格是90×31.首先利用OCCAM进行反演，再进行尖锐边界反演，反演网格是138×40，网格横向间距500 m，

图 1 真实模型 Fig. 1 Real model

在边界位置网格间距扩大，纵向初始层间距50 m，随后以1.1倍递增.OCCAM反演(图 2)和尖锐边界反演(图 3)分别经过6次迭代和25次迭代，拟合差都达到了0.193，由图 2，图 3可以看出，OCCAM反演的刻画边界位置比较模糊，也没有准确反演出异常体确切电阻率，而尖锐边界反演结果能够较好的刻画出上边界的位置，反演电阻率与真实的电阻率也比较接近.由于SBI反演对模型粗糙度约束的是横向层深度变化，对纵向并没有进行约束，所以反演结果左右两边有部分凸凹.

图 2 OCCAM反演结果 Fig. 2 OCCAM inversion results

图 3 SBI反演结果 Fig. 3 SBI inversion results
1.2.2 模型Ⅱ

如图 4是一个阶梯状的模型，阶梯状上方电阻率是1000 Ω·m，阶梯电阻率10 Ω·m，下方电阻率为100 Ω·m，采样点数为9个，采样间距2 km，采集频点有10个，分别为(32，16，8，4，2，1，0.5，0.25，0.125，0.0625)HZ，利用有限元进行模型正演，正演网格大小是90×31.进行SBI反演，反演网格是138×40，网格横向间距500 m，在边界位置网格间距扩大，纵向初始层深度150 m，随后层间距以1.1倍递增.图 5是经过六次Occam反演的结果，图 6是经过33次SBI反演的结果，两种反演拟合差值都达到1.49.

图 4 真实模型 Fig. 4 Real model

图 5 Occam反演结果 Fig. 5 Occam inversion results

图 6 SBI反演结果 Fig. 6 SBI inversion results

由于Occam反演参数个数多于SBI反演参数的个数，Occam反演拟合差值可以继续减小，为了对比两者在相同拟合差情况下的拟合情况，两者期望拟合差均设置为1.49.SBI反演结果能够较好的将地层电性单元分界面区分开.对比Occam反演结果，SBI反演对于第一层位置以及左右边界位置的刻画更好，并且反演的电阻率也能与真实模型电阻率相一致.

1.3 实测CSAMT数据反演

图 7是吉林省临江市六道沟处于吉林省长白山玄武岩覆盖区，位于华北古陆北缘东段.研究区岩石地层单元发育齐全，发育的有太古宙、远古界、古生界、中新生界地层.研究区构造岩活动频繁强烈，以太古宙和中生代岩浆活动最为活跃，主要岩石类型是太古宙奥长花岗岩、花闪长岩、英云闪长岩.工作区内矿产资源丰富，区域内矿产以金、铁、铜铅锌为主要矿种.进行CSAMT测量时候用标量装置，发送机通过两个接地电极A，B向地下供交变电流，实际探测中供电极距AB=1200 m，供电电流是10 A.L2线是该地区所做可控源的测线位置，图 7中白色粗线就是本文所用实际数据的位置.

图 7 吉林省临江市六道沟地区地质图 Fig. 7 Jilin Linjiang Liudaogou region geological map

首先将实测的CSAMT数据进行了全区视电阻率转换，然后对全区视电阻率进行Occam和SBI二维反演，反演结果如下图所示：

分析Occam反演(图 8)和SBI反演(图 9)结果，频率范围是8192 Hz~0.125 Hz，17个频率，18个测点；前17个测点之间间距是200 m，最后两个测点间隔是71 m，最终拟合差RMS都达到1.0.图 8中能看出电阻率的变化范围，但边界处模糊，不容易分辨出地层界面.图 9是SBI反演结果，能够很明显看出地层是一个三层层状构造.根据尖锐边界反演结果分析，第一层的电阻率较低，层厚度约为0.36 km，研究区的踏勘资料显示主要是灰黑色玄武岩，判断其为橄榄玄武岩、安山玄武岩等粗面玄武岩.第二层，界面厚度大约0.3 km，电阻率稍高，再根据地质图，推断该层位为果松组的安山岩，局部地区有流纹岩.第三层电阻率较高，结合本地区地质演化规律，推测该地层是晚期侏罗世的闪长岩.

图 8 Occam反演结果 Fig. 8 Occam inversion results

图 9 SBI反演结果 Fig. 9 SBI inversion results

SBI反演结果显示出了该地区层状构造，地层分界面位置清晰，能够很好地与该地区地质构造吻合.说明了SBI反演能够应用于CSAMT数据，提高了CSAMT数据的勘探深度，对于实际数据处理有很好的指导意义.

2 结 论

2.1  尖锐边界大地电磁测深反演可以避免光滑模型反演中界面不清晰的问题，能够清晰显示出界面的位置.

2.2  尖锐边界大地电磁测深反演效果依赖于给定的初始模型，当给定的初始模型比较差的时候，反演会出现不收敛情况.

2.3  本文所用的全区视电阻率转换方法能够用于CSAMT数据的处理中，并且增加了CSAMT的反演深度，对于实际问题具有应用价值.

参考文献
[1] Constable S C, Parker R L, Constable C G. 2012. Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data[J]. Geophysics, 52(3): 289-300.
[2] de Groot-Hedlin C, Constable S C. 2004. Inversion of magnetotelluric data for 2D structure with sharp resistivity contrasts[J]. Geophysics, 69(1): 78-86.
[3] de Groot-Hedlin C, Constable S C. 2012. Occam's inversion to generate smooth, two-dimensional models from magnetotelluric data[J]. Geophysics, 55(12): 1613-1624.
[4] He M X, Hu X Y, Ye Y X, et al. 2011. 2.5 D controlled source audio-frequency magnetotellurics OCCAM inversion[J]. Progress in Geophysics (in Chinese), 26(6): 2163-2170, doi: 10.3969/j.issn.1004-2903.2011.06.033.
[5] Hu Z Z, Hu X Y. 2005. Review of three dimensional magnetotelluric inversion methods[J]. Progress in Geophysics (in Chinese), 20(1): 214-220, doi: 10.3969/j.issn.1004-2903.2005.01.037.
[6] Lei D. 2010. Studies and applications of 2-D CSAMT modeling and inversion with a dipole source and topography[J]. Chinese Journal Geophysics (in Chinese), 53(4): 982-993, doi: 10.3969/j.issn.0001-5733.2010.04.023.
[7] Li H, Li T L, Wu L, et al. 2015. Transformation of all-time apparent resistivity of CSAMT and analysis of its effect[J]. Progress in Geophysics (in Chinese), 30(2): 889-893, doi: 10.6038/pg20150255.
[8] Li J P, Li T L, Zhao X F, et al. 2007. Study on the TEM all-time apparent resistivity of arbitrary shape loop source over the layered medium[J]. Progress in Geophysics (in Chinese), 22(6): 1777-1780, doi: 10.3969/j.issn.1004-2903.2007.06.015.
[9] Liu G X, Zhang H Z, Han J T, et al. 2006. Features of the electric structure of the lithosphere beneath the Hinggan-Inner Mongolia and Jilin-Heilongjiang regions[J]. Geology in China (in Chinese), 33(4): 824-831.
[10] Ou D X, Wang J L. 2005. The fast magnetotelluric inversion of 2-D block structure[J]. Geophysical Prospecting for Petroleum (in Chinese), 44(5): 525-528.
[11] Shi M J, Xu S Z, Liu B. 1997. Finite element method using quadratic element in Mt forward modeling[J]. Acta Geophysica Sinica (in Chinese), 40(3): 421-430.
[12] Siripunvaraporn W, Egbert G. 2000. An efficient data-subspace inversion method for 2-D magnetotelluric data[J]. Geophysics, 65(3): 791-803.
[13] Smith T, Hoversten M, Gasperikova E, et al. 1999. Sharp boundary inversion of 2D magnetotelluric data[J]. Geophysical Prospecting, 47(4): 469-486.
[14] Tang J T, He J S. 1994. A new method to define the full-zone resistivity in horizontal electric dipole frequency soundings on a layered earth[J]. Acta Geophysica Sinica (in Chinese), 37(4): 543-552.
[15] Tong T G, Liu C M, He J S. 2009. Numerical simulation and application discussion of the CSAMT Full-zone resistivity method[J]. Progress in Geophysics (in Chinese), 24(5): 1855-1860, doi: 10.3969/j.issn.1004-2903.2009.05.041.
[16] Wang J Y. 1997. New development of magnetotelluric sounding in China[J]. Acta Geophysica Sinica (in Chinese), 40(S1): 206-216.
[17] Wannamaker P E, Stodt J A, Rijo L. 1987. A stable finite element solution for two-dimensional magnetotelluric modelling[J]. Geophysical Journal of the Royal Astronomical Society, 88(1): 277-296.
[18] Weng A H. 2007. Occam's inversion and its application to transient electromagnetic method[J]. Geology and Prospecting (in Chinese), 43(5): 74-76.
[19] Wu L, Li T L, Zhu C, et al. 2015. Research and inversion static effect in magnetotelluric[J]. Progress in Geophysics (in Chinese), 30(2): 840-846, doi: 10.6038/pg20150248.
[20] Wu X P, Xu G M. 1998. Improvement of Occam's inversion for Mt data[J]. Acta Geophysica Sinica (in Chinese), 41(4): 547-554.
[21] Ye T, Chen X B, Yan L J. 2013. Refined techniques for data processing and two-dimensional inversion in magnetotelluric (Ⅲ): Using the impressing method to construct starting model of 2D magnetotelluric inversion[J]. Chinese Journal of Geophysics (in Chinese), 56(10): 3596-3606, doi: 10.6038/cjg20131034.
[22] Ye Y X, Hu X Y, KimKangsop, et al. 2009. Application analysis of sharp boundary inversion of magnetotelluric data for 2D structure[J]. Progress in Geophysics (in Chinese), 24(2): 668-674, doi: 10.3969/j.issn.1004-2903.2009.02.041.
[23] Zhang L L, Yu P, Wang J L, et al. 2009. Smoothest model and sharp boundary based two-dimensional magnetotelluric inversion[J]. Chinese Journal of Geophysics (in Chinese), 52(6): 1625-1632, doi: 10.3969/j.issn.0001-5733.2009.06.025.
[24] Zhang L L, Yu P, Wang J L, et al. 2010. A study on 2D magnetotelluric sharp boundary inversion[J]. Chinese Journal of Geophysics (in Chinese), 53(3): 631-637, doi: 10.3969/j.issn.0001-5733.2010.03.017.
[25] Zhao X Y, Weng A H, Li Y B, et al. 2015. The vertical profiles analysis of magnetotelluric sounding response[J]. Progress in Geophysics (in Chinese), 30(3): 1185-1189, doi: 10.6038/pg20150324.
[26] Zhao Y, Li X, Wang Y P. 2015. Full-domain apparent resistivity definition for large-loop TEM[J]. Progress in Geophysics (in Chinese), 30(4): 1856-1863, doi: 10.6038/pg20150446.
[27] 何梅兴, 胡祥云, 叶益信,等. 2011. 2.5维可控源音频大地电磁法Occam反演理论及应用[J]. 地球物理学进展, 26(6): 2163-2170, doi: 10.3969/j.issn.1004-2903.2011.06.033.
[28] 胡祖志, 胡祥云. 2005. 大地电磁三维反演方法综述[J]. 地球物理学进展, 20(1): 214-220, doi: 10.3969/j.issn.1004-2903.2005.01.037.
[29] 雷达. 2010. 起伏地形下CSAMT二维正反演研究与应用[J]. 地球物理学报, 53(4): 982-993, doi: 10.3969/j.issn.0001-5733.2010.04.023.
[30] 李鹤, 李桐林, 伍亮,等. 2015. CSAMT全区视电阻率转换及其效果分析[J]. 地球物理学进展, 30(2): 889-893, doi: 10.6038/pg20150255.
[31] 李建平, 李桐林, 赵雪峰,等. 2007. 层状介质任意形状回线源瞬变电磁全区视电阻率的研究[J]. 地球物理学进展, 22(6): 1777-1780, doi: 10.3969/j.issn.1004-2903.2007.06.015.
[32] 刘国兴, 张志厚, 韩江涛,等. 2006. 兴蒙、吉黑地区岩石圈电性结构特征[J]. 中国地质, 33(4): 824-831.
[33] 欧东新, 王家林. 2005. 二维块状结构大地电磁快速反演[J]. 石油物探, 44(5): 525-528.
[34] 史明娟, 徐世浙, 刘斌. 1997. 大地电磁二次函数插值的有限元法正演模拟[J]. 地球物理学报, 40(3): 421-430.
[35] 汤井田, 何继善. 1994. 水平电偶源频率测深中全区视电阻率定义的新方法[J]. 地球物理学报, 37(4): 543-552.
[36] 佟铁钢, 刘春明, 何继善. 2009. CSAMT全区电阻率法数值模拟及应用探讨[J]. 地球物理学进展, 24(5): 1855-1860, doi: 10.3969/j.issn.1004-2903.2009.05.041.
[37] 王家映. 1997. 我国大地电磁测深研究新进展[J]. 地球物理学报, 40(S1): 206-216.
[38] 翁爱华. 2007. Occam反演及其在瞬变电磁测深中的应用[J]. 地质与勘探, 43(5): 74-76.
[39] 伍亮, 李桐林, 朱成,等. 2015. 大地电磁测深法中静态效应及其反演[J]. 地球物理学进展, 30(2): 840-846, doi: 10.6038/pg20150248.
[40] 吴小平, 徐果明. 1998. 大地电磁数据的Occam反演改进[J]. 地球物理学报, 41(4): 547-554.
[41] 叶涛, 陈小斌, 严良俊. 2013. 大地电磁资料精细处理和二维反演解释技术研究(三)——构建二维反演初始模型的印模法[J]. 地球物理学报, 56(10): 3596-3606, doi: 10.6038/cjg20131034.
[42] 叶益信, 胡祥云, 金钢燮,等. 2009. 大地电磁二维陡边界反演应用效果分析[J]. 地球物理学进展, 24(2): 668-674, doi: 10.3969/j.issn.1004-2903.2009.02.041.
[43] 张罗磊, 于鹏, 王家林,等. 2009. 光滑模型与尖锐边界结合的MT二维反演方法[J]. 地球物理学报, 52(6): 1625-1632, doi: 10.3969/j.issn.0001-5733.2009.06.025.
[44] 张罗磊, 于鹏, 王家林,等. 2010. 二维大地电磁尖锐边界反演研究[J]. 地球物理学报, 53(3): 631-637, doi: 10.3969/j.issn.0001-5733.2010.03.017.
[45] 赵祥阳, 翁爱华, 李亚彬,等. 2015. 大地电磁测深垂直剖面响应分析[J]. 地球物理学进展, 30(3): 1185-1189, doi: 10.6038/pg20150324.
[46] 赵越, 李貅, 王祎鹏. 2015. 大回线源瞬变电磁全域视电阻率定义[J]. 地球物理学进展, 30(4): 1856-1863, doi: 10.6038/pg20150446.