2. 东华理工大学放射性地质与勘探技术国防重点学科实验室, 南昌 330013
2. Key Laboratory of Radioactive Geology and Exploration Technology Fundamental Science for National Defense, East China Institute of Technology, Nanchang 330013, China
可控源音频大地电磁测深法(Controlled-source audio-frequency magnetotellurics,缩写为CSAMT)是在大地电磁法和音频大地电磁法基础上发展起来的一种人工场源频率域测深方法.该方法具有抗干扰能力强、勘探深度范围大、观测效率高等特点(汤井田和何继善,2005),在石油、天然气、地热、固体矿产以及水文和环境工程勘查等领域得到了越来越广泛的应用.
由于CSAMT张量场源复杂的三维特性,为提高数据处理质量和解释水平,迫切需要开展CSAMT多维正反演研究.底青云等(2004)结合直立异常、倾斜异常和断陷模型研究了复杂介质2.5维CSAMT有限元法正演模拟.雷达(2010)使用三角网格剖分建立了起伏地形条件下基于Occam方法的CSAMT二维反演,有效地消除了地形对CSAMT的影响.邓居智等(2011)采用三维交错网格有限差分求解二次电场方程,实现CSAMT三维正演快速迭代求解.韩波等(2015)基于并行化直接求解器MUMPS和交错网格有限体积法实现频率域可控源三维正演,证明了直接求解对于模拟多场源可控源问题的有效性.Hu等(2015)基于六面体网格的矢量有限元法研究了带激电效应的三维张量CSAMT正演,并分析了典型模型的阻抗和视电阻率响应特征.王亚璐等(2017)采用非结构网格离散异常电场方程,引入散度条件约束实现CSAMT有限元数值模拟.
传统的CSAMT正演模拟主要基于各向同性模型,各向异性电阻率的CSAMT问题研究相对较少,大多局限于一维正反演.Li和Pedersen(1992)使用矢量势方法计算层状各向异性介质中的CSAMT响应,数值模拟结果表明即使是适中的各向异性也会产生和各向同性明显不同的结果.Yin和Maurer(2001)研究了CSAMT层状介质任意各向异性正演,分析了各向异性对CSAMT电场和视电阻率的影响.Li等(2000)使用非线性最小二乘法实现层状各向异性地层中CSAMT反演,实测数据反演表明各向异性模型比各向同性模型能更好地拟合数据.Zhou等(2014)研究了层状横向各向同性介质中的CSAMT数据正则化反演,发现CSAMT数据对于横向电阻率更敏感.曹强等(2016)运用Occam方法实现直立各向异性层状介质中CSAMT数据反演.三维正演方面,陈桂波(2009)利用积分方程算法研究了各向异性对CSAMT的场源效应与静态效应的影响.
考虑到在地层层理发育地区各向异性现象非常普遍,人们逐渐认识到电各向异性对大地电磁、海洋电磁和航空电磁等存在很大的影响并积极开展相关领域各向异性模拟研究(Yin and Fraser, 2004; Yin, 2006; Liu and Yin, 2014).相比之下,三维CSAMT各向异性正反演研究相对比较薄弱,各向异性对三维CSAMT的影响尚没有明确的认识.为此,本文研究各向异性介质中三维CSAMT正演模拟算法,通过引入3×3的对称正定矩阵表征各向异性电导率张量,使用基于非结构四面体网格的矢量有限元法离散电场Helmholtz方程,结合矩阵直接求解技术实现三维任意各向异性介质中CSAMT正演模拟.本文首先使用层状各向异性模型验证有限元算法的计算精度,然后分析各向异性异常体和各向异性围岩的CSAMT响应特征,最后提出使用视电阻率极性图识别地下介质各向异性电阻率的主轴方向.
1 控制方程和有限元离散从微分形式的麦克斯韦方程组出发,取时谐因子为eiωt,在有源区域电磁场满足的方程为(Jin, 2014)
(1) |
(2) |
(3) |
式(1)—(3)中,E表示电场强度,H表示磁场强度,Js表示外加电流源项,J表示总电流密度,i是虚数单位,ω=2πf是角频率,f是频率,σ是各向异性的电导率张量,μ是磁导率,ε是介电常数.为了将电导率张量从参考坐标系变换到计算坐标系,对主电导率张量σr进行三次欧拉旋转得到任意各向异性介质的电导率张量(Yin, 2000),即
(4) |
其中主电导率张量σr为
(5) |
式中σx、σy和σz称为主电导率,ρx、ρy和ρz称为主电阻率.(4)式中D=DxDyDz,Dx、Dy和Dz分别为绕不同坐标轴的旋转矩阵(Yin, 2000).
由(1)和(2)式,可得在整个计算区域Ω上电场强度E满足的方程为
(6) |
式中
为了保证电磁场存在唯一解,假设计算区域Ω的外边界Γ距离源足够远,电场强度满足Dirichlet边界条件:
(7) |
其中n是外边界面上单位法向量.
定义
(8) |
选取矢量测试函数V∈H(curl, Ω),对(6)式的第一项运用第一矢量格林定理并进行分部积分,其弱形式(Ren et al., 2013, 2014)表述为
(9) |
其中,
(10) |
(11) |
本文运用矢量有限元法求取电场强度E.为此,将整个计算区域Ω剖分成四面体单元组成的网格
(12) |
其中,
选取矢量测试函数V为矢量形函数N,将(12)式代入(10)和(11)式得到单元矩阵:
(13) |
其中,
(14) |
将(13)和(14)式中的单元矩阵组装可得到如下大型对称稀疏复线性代数方程组:
(15) |
由于张量CSAMT使用两个分离或重叠场源,需要求解两个极化方向的电磁场,直接求解仅需对系数矩阵K进行一次矩阵分解,然后回代两个极化方向的源项即可.因此,本文采用多波前直接求解器MUMPS(Amestoy et al., 2006)对方程组(15)进行求解,得到棱边上的电场Ẽ,进而利用(12)式进行插值求出观测点处的Ẽ,并利用(1)式计算观测点处的磁场Ẽ.
(14)式中的单元系数矩阵Ke的计算可参考(Jin, 2014).为计算使用总场方法时的单元源项矩阵Se,首先将水平导线发射源置于四面体网格中,计算每段导线与每个四面体单元的交点,并将各单元中所包含的线源部分近似为一个电偶极子.假设电偶极子的中点位置位于(x0, y0, z0),则通过单元的电流密度Jse可表示为(齐彦福等,2017)
(16) |
式中Ie表示穿过该单元的源电流,dl表示该段偶极子长度,δ为Dirac-δ函数.由此,(14)式中的源项矩阵可简化为
(17) |
对于任意形状、姿态和位置的电性源,均可利用电偶极子近似,进而由(17)式计算得到源项,(17)式中是矢量点乘,需要计算每个矢量形函数和外加电流源的点乘积.
为使单元电流源项满足电偶极子假设,剖分四面体网格时,在源所在的线段上嵌入一系列等距离点,通过施加网格棱边长度约束实现源附近网格的细化,如图 1a所示.为保证测点附近电磁场计算精度,同时尽可能减少网格单元个数,对测点附近的单元施加棱边长度约束实现网格的局部细化,如图 1b所示.
张量CSAMT通常在远区观测,测区场源可近似为垂直入射的平面波.可根据两次不同极化方向的电场和磁场进行转换得到阻抗张量,阻抗与电场和磁场分量关系式(Yin and Maurer, 2001)为
(18) |
(18)式中上标表示极化模式.由(18)式可得,阻抗张量为
(19) |
式中,ζ=Hy(2)Hx(1)-Hy(1)Hx(2).进一步根据阻抗张量可定义视电阻率和相位:
(20) |
式中i和j可分别表示x或y.
2 算例分析 2.1 层状各向异性模型精度验证首先,设计一个两层各向异性模型验证有限元算法的精度.如图 2所示,第一层为各向同性地层,电阻率为10 Ωm,层厚度为100 m;第二层为各向异性地层,x、y和z方向电阻率分别为200 Ωm、50 Ωm和50 Ωm.以坐标原点为中心沿x轴和y轴方向布设两个互相垂直的源,发射频率为0.1 Hz到10000 Hz,按对数等间隔取21个频点,测点位于(3000, 3000, 0)m处.
将计算结果与Yin和Maurer(2001)的一维层状各向异性半解析解比较,如图 3所示,本文有限元算法结果与Yin和Maurer(2001)一维结果吻合很好,对于等效视电阻率和相位(Yin和Maurer, 2001),两种算法的相对差均在4%以内,验证了本文提出的数值算法的可靠性.
下面研究各向同性半空间中存在各向异性异常体的张量CSAMT响应特征.均匀半空间电阻率为100 Ωm,其中含有一个大小为600 m×600 m×600 m的异常体,异常体中心位于(6500, 13500, 700)m处,其主轴电阻率ρx、ρy、ρz分别为1 Ωm、50 Ωm和1 Ωm.将异常体主电导率张量绕x和z轴分别旋转0°、45°和90°,计算视电阻率ρxya和ρyxa,并分析异常体各向异性对CSAMT视电阻率的影响.
作为对比,首先在图 4中给出电阻率分别为1 Ωm和50 Ωm的各向同性异常体对应的视电阻率响应.由于异常体电阻率较低,根据电流密度在介质分界面的连续性,异常体内电场幅值较小,导致图 4中视电阻率在低阻异常体上方出现极小值,清晰地指示了异常体的水平位置;且异常体电阻率值越小,视电阻率值越低.根据(19)式,ρxya主要由Ex决定,反映x方向电阻率变化.对于图 4a的情况,沿x方向的电流穿过电性分界面时,电流密度连续,而电导率存在很大差异,导致法向电场发生严重的不连续性,因此ρxya等值线沿x方向非常密集. ρyxa的分布特征可作出相似的解释,只是由于此时Ey起决定作用,ρyxa等值线分布特征与ρxya正交.
图 5给出各向异性异常体主电导率张量绕x轴旋转时的张量CSAMT视电阻率响应.从图 5(a-c)可以看出,随着异常体主电导率张量绕x轴旋转角α的增大,ρxya基本不变,且ρxya和图 4a中电阻率为1 Ωm的各向同性模型结果基本一致.这是由于ρxya主要由x方向电场决定,异常体电阻率变化对电场的影响较大而对磁场的影响较小,当异常体主电导率张量绕x轴旋转时,主轴电阻率ρx为1 Ωm保持不变,不影响x方向电场分布,因此ρxya对异常体主电导率张量绕x轴旋转不敏感.从图 5(d-f)可以看出,异常体各向异性对ρyxa有着很大影响,随着旋转角α增大,ρyxa逐渐减小.当α=0°时,ρyxa的幅值和分布规律和图 4d中电阻率为50 Ωm的各向同性模型结果相似;当α=90°时,ρyxa的幅值和分布规律和图 4c中电阻率为1 Ωm的各向同性模型结果相似.这是由于异常体主电导率绕x轴旋转时,主轴电阻率ρy由50 Ωm减小到1 Ωm,进而ρyxa分布特征由50 Ωm的各向同性异常体向1 Ωm的各向同性异常体的视电阻率分布特征转变.
从图 6可以看出,异常体主电导率张量绕z轴旋转对ρxya和ρyxa的幅值和分布规律都产生了很大影响.当旋转角χ=0°时,ρxya与图 4a中电阻率为1 Ωm的各向同性模型结果基本一致,ρyxa与图 4d中电阻率为50 Ωm的各向同性模型结果基本一致.当旋转角χ=90°时,ρxya和ρyxa分别和图 4b中电阻率为50 Ωm和图 4c中电阻率为1 Ωm的各向同性模型结果基本一致.当旋转角在0°到90°(比如图中χ=45°)变化时,ρxya和ρyxa分布特征均随旋转角度发生改变,然而对不同旋转角度的极性图研究发现,ρxya和ρyxa极轴的平均方向指示了异常体各向异性主轴方向.
在层理明显发育地区,围岩也可能呈现各向异性.本节研究各向异性围岩对CSAMT响应的影响特征.模型几何参数与2.2节相同.假设异常体为各向同性、电阻率为1 Ωm;围岩为各向异性,其主轴电阻率ρx、ρy、ρz分别为50 Ωm、100 Ωm和50 Ωm.将围岩主电导率张量分别绕x和z轴旋转0°、45°和90°,计算视电阻率ρxya和ρyxa,并分析围岩不同各向异性时的CSAMT视电阻率响应.
2.3.1 各向异性主轴绕x轴旋转由图 7可以看出,低阻异常体正上方,ρxya和ρyxa的形态和各向同性情况相似,分别在垂直极化电场方向(其中,ρxya在y方向,而ρyxa在x方向)被拉伸.从图 7(a-c)可进一步看出,随着围岩主电导率绕x轴旋转角α的增大,异常体上方视电阻率ρxya增大,围岩上方视电阻率ρxya为50 Ωm不变.这是因为当围岩主电导率张量绕x轴旋转时,围岩主轴电阻率ρx保持不变(50 Ωm),x方向电场分布受影响很小,因此ρxya对围岩主电导率张量绕x轴旋转不敏感;而视电阻率是真电阻率的综合反应,围岩主轴电阻率ρx由50 Ωm变为100 Ωm,导致视电阻率ρxya增大.相比之下,由图 7(d-f)可以看出,随着旋转角α的增大,围岩上方和异常体上方视电阻率ρyxa均减小.其中的物理原因与ρxya相似,这里不做赘述.
从图 8(a-c)可以看出,随着围岩主电导率绕z轴旋转角χ的增大,异常体和围岩上方视电阻率ρxya均增大.这是因为当围岩主电导率张量绕z轴旋转时,围岩主轴电阻率ρx由50 Ωm变为100 Ωm,围岩上方电场增大,视电阻率ρxya增大;视电阻率是异常体和围岩真电阻率的综合反映,围岩主轴电阻率ρx增大,导致异常体上方视电阻率ρxya也增大.同理,从图 8(d-f)可以看出,随着旋转角χ的增大,围岩上方和异常体上方视电阻率ρyxa均减小.当旋转角为γ=45°时,ρxya和ρyxa均沿异常中心发生了偏转,且偏转方向和绕z轴旋转方向相反.和前面讨论得出类似的结论,ρxya和ρyxa极性轴的平均方向指示围岩主电导率的旋转方向.
参考Yin(2006)关于海洋大地电磁的研究,我们将阻抗张量和视电阻率从直角坐标系转换到极坐标系,使用极坐标系下的视电阻率分布(极性图)来识别地下介质的各向异性电阻率主轴方向.设极轴
(21) |
式中,Zrφ和Zφr等价于使用极坐标下径向电场和切向磁场或者切向电场和径向磁场计算的阻抗分量.考虑到Zrφ和Zφr的正交性,这里仅给出ρrφa的计算结果.
图 9展示针对上节三维各向异性异常体模型计算得到的视电阻率ρrφa.从图 9a可以看出,主电导率张量绕x轴由0°旋转到90°时,主轴电阻率ρx保持不变(1 Ωm),因此ρrφa短轴位于x方向,且沿该方向ρrφa保持不变.相比之下,随着旋转角α增大,主轴电阻率ρy由50 Ωm变为1 Ωm,视电阻率长轴随之减小.从图 9b可以看出,当异常体各向异性主轴绕z轴旋转时,ρrφa发生旋转,其长轴方向指示各向异性异常体电阻率较大的主轴方向.
图 10展示针对上节各向异性围岩模型计算的视电阻率响应ρrφa.从图 10a可以看出,当各向异性主轴绕x轴旋转时,主轴电阻率ρx保持不变(50 Ωm),由此ρrφa短轴位于x方向,且该方向上ρrφa保持不变;相比之下,主轴电阻率ρy则由100 Ωm变到50 Ωm,导致视电阻率ρrφa长轴随着旋转角度增大明显变小.进一步从图 10b可以看出,当各向异性主轴绕z轴旋转时,ρrφa发生旋转,但与各向异性异常体的情况不同,此时ρrφa长轴方向垂直于电阻率较大的主轴方向.这可以给出如下解释:当围岩为各向异性时,沿层理方向的电流被低阻异常体吸引产生通道效应,导致垂直层理方向电流向异常体集中,因此围岩中电流密度减小,体现高阻特征.
本文基于非结构四面体网格和矢量有限元方法成功实现任意各向异性介质中张量CSAMT三维高效数值模拟.和层状各向异性模型的半解析解对比验证了本文提出的算法有效性和精确性.各向异性异常体和各向异性围岩模型的数值模拟结果表明,各向异性对视电阻率幅值及其地表分布特征均产生很大影响.极坐标下视电阻率ρrφa极性图能很好地识别各向异性主轴方向.因此,在各向异性效应较为显著的地区开展CSAMT工作必须考虑各向异性的影响.通过采用张量数据采集并结合相应的基于极坐标的视电阻率转换技术可有效识别地下目标体的各向异性特征.本文的研究结果可为认识CSAMT各向异性效应提供依据.如何利用反演技术从张量CSAMT数据中求解出各向异性参数将是我们未来的重要研究方向.
致谢作者向齐彦福博士和吉林大学电磁“千人计划”研究团队其他成员在文章准备过程中提供的帮助表示感谢
Amestoy P R, Guermouche A, L'Excellent J Y, et al. 2006. Hybrid scheduling for the parallel solution of linear systems. Parallel Computing, 32(2): 136-156. DOI:10.1016/j.parco.2005.07.004 |
Cao Q, Lin C H, Tan H D, et al. 2016. Occam inversion of CSAMT data in layered media with azimuthal anisotropy. Geophysical and Geochemical Exploration (in Chinese), 40(3): 594-602. DOI:10.11720/wtyht.2016.3.23 |
Chen G B. 2009. Integral equation method for 3-D electromagnetic modeling in layered anisotropic earth and its applications[Ph. D. thesis] (in Chinese). Changchun: Jilin University.
|
Deng J Z, Tan H D, Chen H, et al. 2011. CSAMT 3D modeling using staggered-grid finite difference method. Progress in Geophysics (in Chinese), 26(6): 2026-2032. DOI:10.3969/j.issn.1004-2903.2011.06.017 |
Di Q Y, Unsworth M, Wang M Y. 2004. 2.5-D CSAMT modeling with the finite element method over 2-D complex earth media. Chinese Journal of Geophysics (in Chinese), 47(4): 723-730. |
Han B, Hu X Y, Huang Y F, et al. 2015. 3-D frequency-domain CSEM modeling using a parallel direct solver. Chinese Journal of Geophysics (in Chinese), 58(8): 2812-2826. DOI:10.6038/cjg20150816 |
Hu Y C, Li T L, Fan C S, et al. 2015. Three-dimensional tensor controlled-source electromagnetic modeling based on the vector finite-element method. Applied Geophysics, 12(1): 35-46. DOI:10.1007/s11770-014-0469-1 |
Jin J M. 2014. The Finite Element Method in Electromagnetics. 3rd ed. Hoboken, NJ: Wiley-IEEE Pres.
|
Lei D. 2010. Studies and applications of 2-D CSAMT modeling and inversion with a dipole source and topography. Chinese Journal of Geophysics (in Chinese), 53(4): 982-993. DOI:10.3969/j.issn.0001-5733.2010.04.023 |
Li X B, Pedersen L B. 1992. Controlled-source tensor magnetotelluric responses of a layered earth with azimuthal anisotropy. Geophysical Journal International, 111(1): 91-103. DOI:10.1111/gji.1992.111.issue-1 |
Li X B, Oskooi B, Pedersen L B. 2000. Inversion of controlled-source tensor magnetotelluric data for a layered earth with azimuthal anisotropy. Geophysics, 65(2): 452-464. DOI:10.1190/1.1444739 |
Liu Y H, Yin C C. 2014. 3D anisotropic modeling for airborne EM systems using finite-difference method. Journal of Applied Geophysics, 109: 186-194. DOI:10.1016/j.jappgeo.2014.07.003 |
Nédélec J C. 1980. Mixed finite elements in R3. Numerische Mathematik, 35(3): 315-341. DOI:10.1007/BF01396415 |
Qi Y F, Yin C C, Liu Y H, et al. 2017. 3D time-domain airborne EM full-wave forward modeling based on instantaneous current pulse. Chinese Journal of Geophysics (in Chinese), 60(1): 369-382. DOI:10.6038/cjg20170131 |
Ren Z Y, Kalscheuer T, Greenhalgh S, et al. 2013. A goal-oriented adaptive finite-element approach for plane wave 3-D electromagnetic modelling. Geophysical Journal International, 194(2): 700-718. DOI:10.1093/gji/ggt154 |
Ren Z Y, Kalscheuer T, Greenhalgh S, et al. 2014. A finite-element-based domain-decomposition approach for plane wave 3D electromagnetic modeling. Geophysics, 79(6): E255-E268. DOI:10.1190/GEO2013-0376.1 |
Wang Y L, Di Q Y, Wang R. 2017. Three-dimensional modeling of controlled-source audio-frequency magnetotellurics using the finite element method on an unstructured grid. Chinese Journal of Geophysics (in Chinese), 60(3): 1158-1167. DOI:10.6038/cjg20170326 |
Yin C, Maurer H M. 2001. Electromagnetic induction in a layered earth with arbitrary anisotropy. Geophysics, 66(5): 1405-1416. DOI:10.1190/1.1487086 |
Yin C C. 2000. Geoelectrical inversion for a one-dimensional anisotropic model and inherent non-uniqueness. Geophysical Journal International, 140(1): 11-23. DOI:10.1046/j.1365-246x.2000.00974.x |
Yin C C, Fraser D C. 2004. The effect of the electrical anisotropy on the response of helicopter-borne frequency-domain electromagnetic systems. Geophysical Prospecting, 52(5): 399-416. DOI:10.1111/gpr.2004.52.issue-5 |
Yin C C. 2006. MMT forward modeling for a layered earth with arbitrary anisotropy. Geophysics, 71(3): G115-G128. DOI:10.1190/1.2197492 |
Zhou J M, Wang J X, Shang Q L, et al. 2014. Regularized inversion of controlled source audio-frequency magnetotelluric data in horizontally layered transversely isotropic media. Journal of Geophysics and Engineering, 11(2): 025003. DOI:10.1088/1742-2132/11/2/025003 |
曹强, 林昌洪, 谭捍东, 等. 2016. 可控源音频大地电磁法方位非各向同性层状介质Occam反演. 物探与化探, 40(3): 594-602. DOI:10.11720/wtyht.2016.3.23 |
陈桂波. 2009. 各向异性地层中电磁场三维数值模拟的积分方程算法及其应用[博士论文]. 长春: 吉林大学.
|
邓居智, 谭捍东, 陈辉, 等. 2011. CSAMT三维交错采样有限差分数值模拟. 地球物理学进展, 26(6): 2026-2032. DOI:10.3969/j.issn.1004-2903.2011.06.017 |
底青云, UnsworthM, 王妙月. 2004. 复杂介质有限元法2.5维可控源音频大地电磁法数值模拟. 地球物理学报, 47(4): 723-730. |
韩波, 胡祥云, 黄一凡, 等. 2015. 基于并行化直接解法的频率域可控源电磁三维正演. 地球物理学报, 58(8): 2812-2826. DOI:10.6038/cjg20150816 |
雷达. 2010. 起伏地形下CSAMT二维正反演研究与应用. 地球物理学报, 53(4): 982-993. DOI:10.3969/j.issn.0001-5733.2010.04.023 |
齐彦福, 殷长春, 刘云鹤, 等. 2017. 基于瞬时电流脉冲的三维时间域航空电磁全波形正演模拟. 地球物理学报, 60(1): 369-382. DOI:10.6038/cjg20170131 |
汤井田, 何继善. 2005. 可控源音频大地电磁法及其应用. 长沙: 中南大学出版社.
|
王亚璐, 底青云, 王若. 2017. 三维CSAMT法非结构化网格有限元数值模拟. 地球物理学报, 60(3): 1158-1167. DOI:10.6038/cjg20170326 |