可控源音频大地电磁法方位非各向同性层状介质Occam反演
曹强1,2,3, 林昌洪1,3, 谭捍东1,3, 龚应丽1,3
1. 中国地质大学(北京) 地下信息探测技术与仪器教育部重点实验室,北京 100083
2. 中国中铁二院工程集团有限责任公司,四川 成都 610031
3. 中国地质大学(北京) 地球物理与信息技术学院,北京 100083

作者简介: 曹强(1989-),男,硕士,中国中铁二院工程集团有限责任公司,主要从事地电学和铁路隧道超前地质预报研究工作。 E-mail: didacaoqiang@163.com

摘要

目前,可控源音频大地电磁法测深资料反演主要基于各向同性介质,在复杂的地电体中,由于岩石裂隙或矿物形状等原因使得在各个方向的电阻率不同。本文主要探讨更为符合实际地质体情况的方位非各向同性层状介质在正交水平电偶源激发条件下的可控源音频大地电磁数据的反演。运用Occam反演方法,实现了可控源音频大地电磁法方位非各向同性层状介质反演。构建理论地电模型,通过理论模型合成数据的反演算例检验了所实现的可控源音频大地电磁法方位非各向同性层状介质反演算法的有效性和稳定性。同时,不同数据的反演结果表明:方位非各向同性介质的反演应同时使用视电阻率 ρs xy ρs yx和相位 ϕxy ϕyx数据。

关键词: CSAMT; 正交水平电偶源; 方位非各向同性; 层状介质; Occam反演
中图分类号:P631 文献标志码:A 文章编号:1000-8918(2016)03-0594-09 doi: 10.11720/wtyht.2016.3.23
Occam inversion of CSAMT data in layered media with azimuthal anisotropy
CAO Qiang1,2,3, LIN Chang-Hong1,3, TAN Han-Dong1,3, GONG Ying-Li1,3
1. Key Laboratory of Geo-detection (China University of Geosciences), Ministry of Education, Beijing 100083, China
2. China Railway ErYuan Engineering Group CO.LTD, Chengdu 610031,China
3. School of Geophysics and Information Technology, China University of Geosciences, Beijing 100083, China
Abstract

At present, the inversion of controlled source audio-frequency magnetotelluric data is mainly based on the isotropic media. Because of the rock fracture or the shape of mineral, the resistivity is different in each direction in the complicated geoelectric body. To fit actual geological conditions better, this paper mainly probes into the inversion of controlled source audio-frequency magnetotelluric data in layered media with azimuthal anisotropy induced by ground orthogonal horizontal electric dipoles. We develop an Occam inversion algorithm to invert controlled source audio-frequency magnetotelluric data for azimuthal anisotropy in layered media. Designing the theoretical geoeletric models, the synthetic examples demonstrate the validity and stability of the Occam inversion algorithm. At the same time, the inversion results of different data suggest that in order to get better inversion results of azimuthally anisotropic layered media the apparent resistivities ρs xy, ρs yx, and phases ϕxy, ϕyx, should be used.

Keyword: CSAMT; orthogonal horizontal electric dipole source; azimuthal anisotropy; layered media; Occam inversion

可控源音频大地电磁法(CSAMT)是在音频大地电磁法(AMT)基础上发展起来的一种人工源频率域测深方法。近年来, CSAMT在寻找金属矿、石油、煤、地热和解决工程地质等问题上有广泛的应用[1, 2]。目前, 实测CSAMT数据主要采用各向同性的反演技术, 然而, 实际地下电性结构是非常复杂的, 许多情况下, 地下介质是方位非各向同性的。例如, 含片状或者针状良导矿物的网脉状或细脉状金属矿石, 沿网脉状或细脉状方向的电阻率值明显低于同等金属矿物含量的浸染状矿石的电阻率值, 而含片状、树枝状高阻矿物(如石英脉)的岩石, 垂直于岩脉方向的电阻率值往往很高[3]。因此, 为了使CSAMT获得更为符合实际地质体情况的结果, 需要展开对地下方位非各向同性介质的反演研究。

对于非各向同性介质, 一维非各向同性介质的研究相对较多, 二维和三维的非各向同性介质的研究 48相对较少。20世纪四五十年代, 前苏联学者就对非各向同性大地中的稳定电流场作了研究[9]; 国内学者付良魁在1964年对层状非各向同性大地的稳定电流场作了理论探讨[10]; 阮爱国、毛桐恩等在2002年推导了直流电法层状方位各向异性介质中的电位分布、边界条件及视电阻率计算公式[11]; 沈金松、郭乃川等对各向异性层状介质中直流电视电阻率响应进行过探讨[12, 13]。在大地电磁法(MT)中, 早期F. Abramovici和Y. Shoham将广义反演应用到一维各向异性模型的MT资料[14]; 林长佑, 武玉霞等在1996年对水平层状对称各向异性大地电磁资料进行了反演研究[15]; 林长佑, 杨长福等对大地电磁法各向异性层状介质做了进一步研究[16, 17], 并运用于实测数据; 2003年, Changchun Yin研究了一维层状各向异性模型MT反演的内在非唯一性[18] ; 2006年, Josef Pek 和Fernando A. M. Santos利用一维Occam反演方法, 实现MT各向异性层状介质反演[19]; 秦林江、杨长福采用广义逆法对一维MT水平层状各向异性介质模型反演进行了探索研究, 实现了MT全张量响应的一维各向异性反演[20] ; 霍光谱、胡祥云等利用带加权因子的马奎特反演对层状各向异性介质大地电磁视电阻率和相位进行反演研究[21]。在可控源音频大地电磁法(CSAMT)中, Li等在1992年研究分析了可控源张量大地电磁法方位非各向同性半空间和层状介质的阻抗张量和倾斜因子[22, 23]; 同年, Yu和R. N. Edwards开发了计算海底任意有限源多层横向各向异性介质电磁场的算法[24]; Yin和H.-M. Maurer在2001年对CSAMT各向异性层状介质的电磁感应和视电阻率响应进行了研究[25], 说明地表视电阻率和相位对各向异性具有很好的灵敏性; 2007年, L. O. Loseth和B. Ursin研究并推导了偶极子源下电磁波在各向异性层状介质下的传播[26]; Li等在2000年利用阻尼最小二乘法, 实现了对一维方位非各向同性层状介质可控源张量大地电磁阻抗张量数据和倾子矢量数据的反演[27]。笔者在正交水平电偶源激发条件下, 将讨论一种使用Occam方法[28]反演卡尼亚视电阻率和相位资料, 得到方位非各向同性层状地电模型的反演方法。

1 正演问题
1.1 理论

考虑有M层的方位非各向同性的地层(如图1), 水平地层两个正交方向的电阻率不同, 根据此正交方向建立笛卡尔坐标系, 定义准静态条件(即忽略位移电流的影响)下的张量电阻率在第i层地层中为

ρi=ρtiρniρti,

式中, ρ tix方向的电阻率, ρ niy方向的电阻率, 并定义非各向同性系数λ i= ρni/ρti

i层层状非各向同性和非磁性介质的频率域麦克斯韦方程组为

式中, Ei是电场总矢量, Hi是磁场总矢量, ρ i是张量电阻率, ω 是角频率, μ 0为真空磁导率。引入矢量势A和标量势Ψ , 借助二维傅里叶变换对将矢量势A的各个分量从空间域转换到波数域, 利用快速汉克尔变换数值算法, 可求解得到x方向电偶源产生的地表总电场Ex和地表总磁场Hx, 同时可求解得到y方向电偶源产生的地表总电场Ey和地表总磁场 Hy23

图1 M层方位非各向同性介质模型

利用张量阻抗通用计算公式[29]

Zxx=ExxHyy-ExyHyxHxxHyy-HxyHyx, Zxy=ExyHxx-ExxHxyHxxHyy-HxyHyx, Zyx=EyxHyy-EyyHyxHxxHyy-HxyHyx, Zyy=EyyHxx-EyxHxyHxxHyy-HxyHyx,

式中, ExxEyxExyEyy分别表示xy方向电偶源产生的地表x方向和y方向总电场, HxxHyxHxyHyy分别表示xy方向电偶源产生的地表x方向和y方向总磁场。参考卡尼亚视电阻率的定义, 可得到计算地表视电阻率和相位的表达式

ρsxx=1μ0ω|Zxx|2,  ρsxy=1μ0ω|Zxy|2, ρsyx=1μ0ω|Zyx|2,  ρsyy=1μ0ω|Zyy|2, ϕxx=Arg(Zxx),  ϕxy=Arg(Zxy), ϕyx=Arg(Zyx),  ϕyy=Arg(Zyy)

1.2 方位非各向同性层状介质的阻抗正演响应

三层方位非各向同性层状介质地层参数见表1, 正交水平电偶源位于坐标原点, 接收点坐标为(7 070 m, 7 070 m), 接收频率13个(21~213 Hz)。三层方位非各向同性层状介质的阻抗正演响应如图2所示, 由图可知, ZxxZyy的实部值、虚部值都接近于零, 故笔者在反演时没有利用ZxxZyy对应的视电阻率和相位数据。

图2 三层方位非各向同性层状介质的阻抗正演响应

表1 三层方位非各向同性层状介质地层参数
1.3 方位非各向同性层状介质与各向同性层状介质的正演响应对比

三层方位非各向同性层状介质地层参数见表1, 三层各向同性层状介质地层参数见表2。准静态条件下, 在远区的各向同性层状介质的CSAMT视电阻率值的关系是ρ sxysyx, 且ρ sxxsyy=0。三层方位非各向同性层状介质与三层各向同性层状介质的视电阻率和相位正演响应对比结果, 如图3所示。由图3可知, 三层非各向同性层状介质的视电阻率 ρsxyanisotropic和相位 ϕxyanisotropic响应曲线与三层各向同性层状介质的视电阻率 ρsxyisotropic和相位 ϕxyisotropic响应曲线基本重合, 没有明显差异, 但是, ρsyxanisotropicϕyxanisotropic的响应曲线与 ρsyxisotropicϕyxisotropic的响应曲线相比, 存在明显的幅值差。从以上正演响应对比结果可以看出, 对于地下层状介质, 如果实测的视电阻率ρ sxyρ syx, 以及二者的相位存在明显差异, 应考虑地下介质可能存在方位非各向同性。

表2 三层各向同性层状介质地层参数

图3 三层方位非各向同性层状介质与三层各向同性层状介质的正演响应对比结果

2 反演问题
2.1 目标函数

反演采用Occam方法, Occam反演中除了要把拟合数据差控制到期望值以内, 反演的模型应尽可能简单、光滑, 使模型粗糙度最小。考虑模型粗糙度, 同时考虑数据的拟合差, 引入拉格朗日乘子μ -1, 得到目标函数:

ψ(m)=m2+μ-1Wd-WF(m)2-X* 2}

式中, ‖ ∂ m2= (dm/dz)2dz是粗糙度, 是模型m对深度z的一阶导数平方的积分的范数形式; d是数据向量; F(m)是正演响应; m是模型参数向量; WM× M的对角矩阵, 它的作用是利用数据标准差对数据进行归一化; X* 为数据拟合差的期望值。

文中反演输入数据使用了3种数据向量, 具体形式为

d1=(lgρsxyi, φxyi)T, d2=(lgρsyxi, φyxi)T, d3=(lgρsxyi, φxyi, lgρsyxi, φyxi)T,

式中, i是观测频点数。三种反演数据得到不同反演结果, 并对其结果进行对比。

2.2 灵敏度矩阵计算

在反演问题中, 最为关键的部分在于灵敏度矩阵Jij的计算。考虑模型正演响应Fi(m), 对其在mj处进行泰勒级数一阶展开

Fi(mj+1)=Fi(mj+δm)=Fi(mj)+Jij(mj+1-mj),

式中, Jij是灵敏度矩阵Jij=∂ Fi(mj)/∂ mj , 本文采用中心差分法计算灵敏度矩阵, 即

Jij=Fi[m]mj=Fi(mj+δm)-Fi(mj-δm)2δm

2.3 反演流程

可控源音频大地电磁法方位非各向同性层状介质Occam反演算法大致流程(见图4)为:

1) 构建初始模型, 输入反演数据。

2) 计算模型的正演响应F(m)灵敏度矩阵J和模型粗糙度‖ ∂ m2

3) 计算目标函数ψ (m)和其关于m的微分, 得到

mk+1=[μT+(WJk)T(WJk)]-1(WJk)TWdk,

其中,

dk=d-F(mk)+Jkmk

4) 采用一维线性搜索方法搜索拉格朗日乘子μ , 搜索范围1~105, 最小化拟合差X2

5) 拟合差是否满足期望值 X* 2, 如果满足, 选择最大μ 的模型保存, 退出程序; 否则, 选择X2最小的模型保存并继续执行程序。

6) 是否达到最大迭代次数, 如果是, 退出程序; 否则, 继续执行程序。

7) 将步骤5保存得到的模型作为新的初始模型, 返回到步骤2继续执行程序。

图4 可控源音频大地电磁法方位非各向同性层状 介质Occam反演算法流程

3 理论模型合成算例
3.1 三层模型

3.1.1 各向同性正演模型合成数据的反演结果

为了说明算法可行性, 建立一个各向同性正演模型, 模型地层参数见表2, 并对其合成数据进行非各向同性反演。反演输入数据使用d3=(lg ρsxyi, φxyi, lg ρsyxi, φyxi)T进行反演。反演初始模型共20层, 每层x方向电阻率ρ t=100 Ω · m, 方位非各向同性系数λ =1, 即各向同性地层, 地层厚度25 m, 正交水平电偶源位于坐标原点, 接收点坐标为(7 070 m, 7 070 m), 接收频率13个(21~213 Hz)。在理论模型的视电阻率和相位中加入1%的高斯随机误差, 用可控源音频大地电磁法方位非各向同性层状介质Occam反演程序在PC 机上进行反演。

各向同性正演模型合成数据的非各向同性反演结果如图5所示, 每层x方向电阻率ρ t的反演结果与各向同性理论模型吻合, 每层非各向同性系数λ 也接近于1, 是各向同性地层。反演结果模型的响应与反演输入数据的拟合结果, 如图6所示, 视电阻率ρ sxysyx, 相位ϕ xyyx, 二者拟合良好, 说明笔者提出的算法是正确、可行的。

图5 三层各向同性正演模型合成数据的反演结果

图6 三层各向同性反演结果模型的数据拟合结果

3.1.2 方位非各向同性正演模型合成数据的反演结果

理论方位非各向同性层状介质模型地层参数见表1。初始模型20层, 每层x方向电阻率ρ t=100 Ω · m, 方位非各向同性系数λ =1, 即各向同性地层, 地层厚度25 m, 正交水平电偶源位于坐标原点, 接收点坐标(7 070 m, 7 070 m), 接收频率13个(21~213 Hz)。在理论模型的视电阻率和相位中加入1%的高斯随机误差后, 用可控源音频大地电磁法方位非各向同性层状介质Occam反演程序在PC 机上进行反演。

对合成数据进行方位非各向同性反演的反演结果如图7所示。由图7a可知, 当反演输入数据是d1=(lg ρsxyi, φxyi)T时, x方向的电阻率ρ t反演结果和理论模型吻合, 非各向同性系数λ 反演结果第一层偏大, 第二层偏小, 第三层偏大; 由图7b可知, 当反演输入数据是d2=(lg ρsyxi, φyxi)T时, x方向的电阻率ρ t反演结果第一层偏小, 第二层和第三层偏大, 非各向同性系数λ 反演结果第一层偏大, 第二层偏小, 第三层先大后小; 由图7c可知, 当反演输入数据是d3=(lg ρsxyi, φxyi, lg ρsyxiφyxi)T时, x方向的电阻率ρ t及非各向同性系数λ 反演结果均与理论模型吻合, 且两者的反演结果明显好于利用单个卡尼亚视电阻率和相位数据的反演结果。

反演结果模型响应与反演输入数据拟合结果, 如图8所示。二者拟合良好, 说明反演结果可靠。

图7 三层方位非各向同性不同正演模型合成数据的反演结果

图8 三层方位非各向同性反演结果模型的数据拟合结果

3.2 五层模型

理论模型地层参数见表3, 初始模型40层, 每层x方向电阻率ρ t=100 Ω · m, 非各向同性系数λ =1, 即为各向同性地层, 地层厚度25 m, 正交水平电偶源位于坐标原点, 接收点坐标(7 070 m, 7 070 m), 接收频率13个(21~213 Hz)。在理论模型的视电阻率和相位中加入1%的高斯随机误差后, 用可控源音频大地电磁法方位非各向同性层状介质Occam反演程序在PC 机上进行反演。

表3 五层理论模型地层参数

五层方位非各向同性层状介质反演结果如图9所示。由图9a可知, 当反演输入数据是d1=(lg ρsxyi, φxyi)T时, x方向的电阻率ρ t反演结果与理论模型吻合, 非各向同性系数λ 反演结果却很差, 无法分辨层位, 且幅值整体偏大; 由图9b可知, 当反演输入数据是d2=(lg ρsyxi, φyxi)T时, x方向的电阻率ρ t反演结果第一、二层偏小, 第四层偏大, 第三层和第五层与理论模型吻合; 非各向同性系数λ 反演结果整体趋势和理论模型一样, 但是反演结果的深度位置比理论模型的深度位置整体偏低。由图9c可知, 当反演输入数据是d3=(lg ρsxyi, φxyi, lg ρsyxi, φyxi)T时, x方向的电阻率ρ t反演结果与理论模型吻合, 非各向同性系数λ 反演结果也与理论模型吻合。

反演结果模型响应与反演输入数据拟合结果, 如图10所示。二者拟合良好, 说明反演结果可靠。

图9 五层方位非各向同性不同正演模型合成数据的反演结果

图10 五层反演结果模型的数据拟合结果

4 结论

利用Occam反演方法, 实现了可控源音频大地电磁法方位非各向同性层状介质视电阻率和相位资料的反演, 通过理论模型合成数据的反演算例检验了可控源音频大地电磁法方位非各向同性层状介质反演算法的正确性和有效性。方位非各向同性层状介质的正演响应与各向同性层状介质的正演响应存在明显差异。对比正交水平电偶源激发条件下所产生的不同数据的反演结果发现, 对于方位非各向同性层状介质的反演, 同时利用视电阻率ρ sxyρ syx和相位ϕ xyϕ yx的反演结果最可靠, 最接近理论模型。这说明, 在考虑地下介质方位非各向同性情况下, 为了得到更为可靠的地下电性结构, 应在可控源音频大地电磁法实际工作中采用正交水平电偶源, 同时对地表电场和磁场进行多分量的观测, 并在反演中尽可能地使用视电阻率ρ sxyρ syx和相位ϕ xyϕ yx数据。

The authors have declared that no competing interests exist.

参考文献
[1] Routh P S, Oldenburg D W. Inversion of controlled source audio-frequency magnetotellurics data for a horizontally layered earth[J]. Geophysics, 1999, 64(6): 1689-1697. [本文引用:1]
[2] Lu X Y, Unsworth J M, Booker J. Rapid relaxation inversion of CSAMT data[J]. Geophysical journal International, 1999, 138(2): 381-392. [本文引用:1]
[3] 傅良魁. 应用地球物理教程—电法放射性地热 [M]. 北京: 地质出版社, 1991. [本文引用:1]
[4] 杨长福. 大地电磁二维对称各向异性介质有限元数值模拟[J]. 西北地震学报, 1997, 19(2): 27-33. [本文引用:1]
[5] Li Y G, Dai S K. Finite element modelling of marine controlled-sourceelectromagnetic responses in two-dimensional dipping anisotropicconductivity structures[J]. Geophys. J. Int. , 2011, 185(1): 622-636. [本文引用:1]
[6] 王威, 吴小平. 电阻率任意各向异性三维有限元快速正演[J]. 地球物理学进展, 2010, 25(4): 1365-1371. [本文引用:1]
[7] Plotkin V V. Inversion of heterogeneous anisotropic magnetotelluric responses[J]. Russian Geology and Geophysics, 2012, 53: 829-836. [本文引用:1]
[8] 陈桂波, 汪宏年, 姚敬金, . 各向异性海底地层海洋可控源电磁响应三维积分方程法数值模拟[J]. 物理学报, 2009, 58(6): 3848-3857. [本文引用:1]
[9] 克维亚特可夫斯基. 电法勘探[M] . 北京: 地质出版社, 1958. [本文引用:1]
[10] 付良魁. 论在岩石非各向同性条件下电测资料解释的理论基础[G]//北京地质学院科学研究论文集(二), 北京地质学院, 1964. [本文引用:1]
[11] 阮爱国, 毛桐恩, 李清河, . 层状方位各向异介质的视电阻率计算[J]. 地震学报, 2002, 24(5): 502-509. [本文引用:1]
[12] 沈金松, 郭乃川. 各向异性层状介质中视电阻率与磁场响应研究[J]. 地球物理学报, 2008, 51(5): 1608-1619. [本文引用:1]
[13] 沈金松, 郭乃川, 苏本玉. 直流电测深中各向异性层状地层的视电阻率响应特征[J]. 中国石油大学学报, 2009, 33(3): 59-65. [本文引用:1]
[14] Abramovici F, Shoham Y. Inversion of anisotropic magnetotelluric data[J]. Geophys. J. R. Astr. Soc. , 1977, 50(1): 55-74. [本文引用:1]
[15] 林长佑, 武玉霞, 杨长福, . 水平层状对称各向异性大地电磁资料反演[J]. 地球物理学报, 1996, 39(S): 342-348. [本文引用:1]
[16] 林长佑, 王书明, 孙崇赤, . 天祝—永登地区大地电磁资料的一维各向异性反演[J]. 西北地震学报, 2004, 26(1): 72-77. [本文引用:1]
[17] 杨长福, 林长佑, 孙崇赤, . 二维对称各向异性介质大地电磁反演[J]. 地震学报, 2005, 27(3): 339-345. [本文引用:1]
[18] Yin C C. Inherent nonuniqueness in magnetotelluric inversion for 1D anisotropic models[J]. Geophysics, 2003, 68(1): 138-146. [本文引用:1]
[19] Pek J, Santos F A M. Magnetotelluric inversion for anisotropic conductivities in layered media[J]. Physics of the Earth and Planetary Interiors, 2006, 158(2-4): 139-158. [本文引用:1]
[20] 秦林江, 杨长福, 孙崇赤, . 大地电磁全张量响应的一维各向异性反演[J]. 地球物理学报, 2012, 55(2): 693-701. [本文引用:1]
[21] 霍光谱, 胡祥云, 方慧, . 层状各向异性介质大地电磁联合反演[J]. 物理学报, 2012, 61(12): 129101. [本文引用:1]
[22] Li X B, Pedersen L B. The electromagnetic response of an azimuthally anisotropic half-space[J]. Geophysics, 1991, 56(9): 1462-1473. [本文引用:1]
[23] Li X B, Pedersen L B. Controlled-source tensor magnetotelluric responses of layered earth with azimuthal anisotropic[J]. Geophys. J. Int. , 1992, 111(1): 91-103. [本文引用:1]
[24] Yu L M, Edwards R N. Algorithms for the computation of the electromagnetic response of amultilayered, laterally anisotropic sea floor to arbitrary finite sources[J]. Geophys. J. Int. , 1992, 111(1): 185-189. [本文引用:1]
[25] Yin, Maurer H-M. Electromagnetic induction in a layered earth with arbitrary anisotropy[J]. Geophysics, 2001, 66(5): 1405-1416. [本文引用:1]
[26] Loseth L O, Ursin B. Electromagnetic fields in planarly layered anisotropic media[J]. Geophys. J. Int. , 2007, 170(1): 44-80. [本文引用:1]
[27] Li X B, Oskooi B, Peserent L B. Inversion of controlled-source tensor magnetotelluric data for a layered earth with azimuthal anisotropy[J]. Geophysics, 2000, 65(2): 452-464. [本文引用:1]
[28] Constable S C, Parker R L, Constable C G. Occam's inversion : A practical algorithm for generating smooth models form electromagnetic sounding data[J]. Geophysics, 1987, 52(3): 289-300. [本文引用:1]
[29] 谭捍东, 魏文博, 邓明, . 大地电磁法张量阻抗通用计算公式[J]. 石油地球物理勘探, 2004, 39(1): 113-116. [本文引用:1]
[30] 刘鸿洲, 安亚婷. 倾斜电各向异性介质大地电磁正演研究[J]. 物探与化探, 2014, 38(6): 1241-1245. [本文引用:1]