基于非结构网格的带地形MT二维Occam反演及应用

doi:10.11720/wtyht.2022.1397

基于非结构网格的带地形MT二维Occam反演及应用

姜奋勇^,, 叶益信^,, 陈海文, 杨烁健

东华理工大学地球物理与测控技术学院,江西南昌 330013

Application of 2D inversion of magnetotelluric data bearing terrain information based on an unstructured mesh

JIANG Fen-Yong^,, YE Yi-Xin^,, CHEN Hai-Wen, YANG Shuo-Jian

School of Geophysics and Measurement-Control Technology, East China University of Technology,Nanchang 330013,China

通讯作者: 叶益信(1983-),男,博士,副教授,硕士生导师,主要从事电法勘探正反演研究及其应用工作。Email:yixinye321@126.com

责任编辑: 王萌

收稿日期: 2021-07-16 修回日期: 2021-09-25

基金资助:

国家自然科学基金项目“带复杂海底地形的海洋可控源电磁法三维正反演研究”(41774078)
江西省自然科学基金青年重点项目“基于非结构双网格的频率域可控源电磁法三维快速反演研究”(20202ACBL211006)

Received: 2021-07-16 Revised: 2021-09-25

作者简介 About authors

姜奋勇(1996-),男,研究方向为电法勘探正反演。Email: stujiangfy@163.com

摘要

本文对基于自适应非结构三角形网格的带地形MT数据二维Occam反演进行应用研究。自适应非结构三角形网格能够准确地模拟起伏地形和复杂地质构造,正演网格由有限元解的后验误差估计指导自动细化,保证了模型响应的准确性,反演网格在反演目标区域采用精细网格剖分,在模型边界区域采用粗网格剖分,在满足反演精度的前提下减少了不必要的反演参数。基于快速Occam正则化反演方法,采用伴随方程法推导灵敏度矩阵并结合加权平方法计算模型粗糙度。通过对陆地起伏地形模型和起伏海底模型反演试算,验证了算法的精确性和适用性,能够对陆地及海底起伏地形下的多尺度结构进行成像。然后,将该方法应用于克拉玛依后山区域的实测数据反演处理中,反演得到的电阻率结构与地质资料以及采用非线性共轭梯度反演所得结果相吻合。

关键词： 非结构网格; MT; 二维Occam反演; 地形; 克拉玛依

Abstract

This paper focuses on the application of the 2D inversion of magnetotelluric data that bear terrain information based on an unstructured adaptive triangular mesh. An adaptive unstructured triangular grid can be used to accurately simulate undulating terrain and complex geological structures. The adaptive unstructured triangular grids for magnetotelluric forward modeling are automatically refined using the a posteriori error estimation with the finite element solution, which ensures the accuracy of the model response. As for adaptive unstructured triangular grids for magnetotelluric inversion, fine mesh generation is adopted for the inversion target areas, while coarse grid generation is utilized for the boundary areas of the model, thus reducing unnecessary inversion parameters on the premise of satisfying the inversion accuracy. According to the inversion of an undulating-terrain model of land and an undulating-seabed model, the accuracy and applicability of the algorithm are verified and the algorithm can be used to image the multi-scale structures under the undulating terrain of land and seabed. Then, the method was applied to the inversion of the measured data of the Houshan area in Karamay. As a result, the resistivity structure obtained through the inversion was consistent with the geological data and the results obtained through the nonlinear conjugate gradient inversion.

Keywords： unstructured mesh; MT; 2D inversion; terrain; Karamay

PDF (2133KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

姜奋勇, 叶益信, 陈海文, 杨烁健. 基于非结构网格的带地形MT二维Occam反演及应用[J]. 物探与化探, 2022, 46(2): 482-489 doi:10.11720/wtyht.2022.1397

JIANG Fen-Yong, YE Yi-Xin, CHEN Hai-Wen, YANG Shuo-Jian. Application of 2D inversion of magnetotelluric data bearing terrain information based on an unstructured mesh[J]. Geophysical and Geochemical Exploration, 2022, 46(2): 482-489 doi:10.11720/wtyht.2022.1397

0 引言

在过去的几十年中,大地电磁法被广泛应用于深部构造研究以及油气和矿产勘查中^[1⇓-3]。常规的MT数据处理和解释往往忽略地形因素的影响,而在实际MT数据采集中,难以避免的会遇到地形起伏的问题,而起伏地表会引起MT响应的畸变^[4-5]。因此,将地形因素引入到大地电磁数值模拟中是非常必要的,许多学者利用有限元法对起伏地形下二维大地电磁响应进行了数值模拟研究,并对畸变后的大地电磁响应做了进一步的分析^[3,6]。常规基于规则网格的算法对复杂起伏地形的模拟精度不高,而且使用起来不方便,因此越来越多的学者采用非结构三角网格有限元方法计算电磁场^[7⇓-9]。非结构三角网格能够准确地模拟起伏地形和复杂地质构造,并且可以通过自适应有限元法自动调节网格,以保证正演响应的精度^{[10⇓⇓⇓⇓-15]}。已有研究将非结构三角网格自适应有限元法用于求解大地电磁的响应,取得了很好的效果,并结合Occam反演算法成功用于带地形大地电磁数据二维反演中,证明了其算法具有较强的实用性和较高的精度^[16-17]。在一些研究区内地形比较复杂,水平地形的反演往往会产生一些不准确的结果,要实现高质量的探测,就要对研究区域采集的MT数据进行带地形的反演,基于非结构网格自适应有限元的带地形反演结果将会提供更加准确的电阻率信息,这对以后的工程勘查和矿产勘探是很有帮助的。

本文对基于自适应非结构网格有限元的MT二维带地形反演做了应用研究,对起伏地形模型进行了MT二维反演试算,得到了比较好的结果。最后将本算法应用于克拉玛依后山区域实测MT数据的反演,同样取得了较好的结果。

1 反演理论

1.1 Occam反演

MT数据反演方法采用的是Occam算法^[18-19],该算法是基于正则化约束的最小二乘法,其泛函U表达式为:

(1)

U = μ ‖ R m ‖^{2} + ‖ W (d - F (m)) ‖^{2}

式中:m为N维的模型参数矢量,一般为电阻率值;R为粗糙度算子矩阵;μ为拉格朗日乘子,用以平衡模型粗糙度和数据拟合误差,当 μ 取较大值时,反演以搜索光滑模型为主,反之则以搜索最小拟合误差为主;W为关联的对角加权矩阵;d为观测数据矢量,F(m)为模型m对应的正演响应。

给定初始模型m_k的函数,通过以下迭代方法实现目标函数的最小化:

(2)

m_{k + 1} = [μ R^{T} R + (W J_{k})^{T} W J_{k}]^{- 1} \times (W J_{k})^{T} W d

其中,修改数据向量

(3)

\dot{d} = d - f (m_{k}) + J_{k} m_{k}

雅可比矩阵J为M×N大小的矩阵,其中每个分量的表达式为:

(4)

J_{i j} = \frac{\partial F_{i} (m)}{\partial m_{j}} 。

1.2 灵敏度计算

灵敏度矩阵的计算对于电磁反演是非常重要的步骤,能够精确求解灵敏度矩阵,对于最终的电磁数据反演分辨率至关重要。求解灵敏度矩阵的方法有很多,最经典的方法是将问题线性化,然后迭代求解^[20]。而在我们的反演问题中,往往模型参数的数量要远大于数据的数量,所以这种方法比较低效,本文采用的是伴随方程法推导灵敏度矩阵^[21-22]。

频率域麦克斯韦方程组表示为:

(5)$\left\{\begin{array}{l}\nabla \times E=-\mathrm{i} \omega \mu H+J_{m}, \\ \nabla \times H=\sigma E+J_{e},\end{array}\right.$

式中:ε为介电常数;μ为磁化率;σ为可变电导率;E为电场强度;H为磁场强度;J_e和J_m分别为电源和磁源。

设置合适的边界条件,将ϕ_j(r)作为一个基函数,则电导率被表示为一组基函数的线性组合:

(6)

σ (r) = \overset{N}{\sum_{j = 1}} σ_{j} ϕ_{j} (r)

式中:σ_j为系数。将式(6)带入方程组(5),然后对σ_j求偏导,得到灵敏度方程:

(7)$\left\{\begin{array}{l}\nabla \times \frac{\partial E}{\partial \sigma_{j}}=-\mathrm{i} \omega \mu \frac{\partial H}{\partial \sigma_{j}} \\ \nabla \times \frac{\partial H}{\partial \sigma_{j}}=\sigma \frac{\partial E}{\partial \sigma_{j}}+\phi_{j} E_{\circ}\end{array}\right.$

式中: $\frac{\partial E}{\partial σ_{j}}$ 和 $\frac{\partial H}{\partial σ_{j}}$ 是电场和磁场对第j个单元电导率的灵敏度。

引入伴随电源 $J_{e}^{*}$ 和伴随磁源 $J_{m}^{*}$ ,得到如下麦克斯韦方程组:

(8)$\left\{\begin{array}{l}\nabla \times E^{*}=-\mathrm{i} \omega \mu H^{*}+J_{m}^{*}, \\ \nabla \times H^{*}=\sigma E^{*}+J_{e}^{*} \circ\end{array}\right.$

考虑边界问题并带入其中,得到下式:

(9)

\int_{D} (J_{m}^{*} \cdot \frac{\partial H}{\partial σ_{j}} + J_{e}^{*} \cdot \frac{\partial E}{\partial σ_{j}}) d v = \int_{D} E^{*} \cdot E ϕ_{j} d v

方程(9)便是求解电磁场灵敏度的基本公式,通过选择合适的伴随场源计算灵敏度矩阵。如果需要获得观测点r₀处E_x的灵敏度矩阵,令 $J_{e}^{*}$ =σ(r-r₀) $\dot{e}$ 和 $J_{m}^{*}$ =0,则得到:

(10)

{\frac{\partial E_{x}}{\partial σ_{j}}|}_{r_{0}} = \int_{D} E^{*} \cdot E ϕ_{j} d v

式(10)即为E_x分量的灵敏度计算公式,为了得到 $\frac{\partial E_{x}}{\partial σ_{j}}$ 可以通过求解原始电磁场得到E,然后计算得到伴随场源情况下的电场值E^*,最后点积E^*·E即可。

1.3 粗糙度计算

模型粗糙度算子R主要用于反演结果的稳定,避免产生虚假的电性结构,目前比较常用的是将模型梯度的L₂范数作为模型的粗糙度算子^[23-24],其表达式为:

(11)$\|R(m)\| \int_{\Omega} \nabla m^{2} \cdot{ }_{\mathrm{d}} \Omega_{\circ}$

基于结构网格的对积分的离散近似较为容易,而在非结构网格中情况却较为复杂,因此本文采用加权平方和的方法^[17]处理梯度点积的积分:

(12)

‖ R (m) ‖^{2} = \overset{m}{\sum_{i = 1}} A_{i} [\overset{N (i)}{\sum_{j = 1}} ω_{j} {(\frac{Δ m_{i j}}{Δ r_{i j}})}^{2}]

其中:

(13)

\begin{array}{l} Δ m_{i j} = m_{i} - m_{j}, \end{array}

(14)

\begin{array}{l} Δ r_{i j} = \sqrt{(y_{i} - y_{j})^{2} + (z_{i} - z_{j})^{2}}, \end{array}

(15)

\begin{array}{l} ω_{j} = \frac{V_{j}}{\overset{N (i)}{\sum_{k = 1}} V_{k}}, \end{array}

式中:V_i表示第i个三角单元的面积;N_i表示与第i个非结构三角单元共享一个顶点所有三角单元的集合;Δr_ij表示三角形之间质心的距离。

2 模型实例

2.1 陆地起伏地形模型

为了验证该方法反演的效果,本文首先设计了一个起伏地形的复杂模型。模型结构如图1所示,在不同的地层中分别设置了两个低阻异常体,1号低阻异常体电阻率为3 Ω·m,2号低阻异常体电阻率为0.1 Ω·m。MT接收点从左往右依次设置了21个测点,如图1中白色小三角形所示,每个测点间距800 m。MT的频率范围为0.001~1 000 Hz,共选取了48个频点,取对数等间隔分布。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 陆地起伏地形模型

Fig.1 Land complex topography model

对上述起伏地形模型进行自适应有限元正演计算,得到模型正演响应与自适应细化网格,如图2所示,图2a为第一次自适应细化网格;图2b为第五次自适应细化网格,表1为细化次数及每次细化生成的节点和单元数量。可以看出,由于测点附近的网格剖分对测点接收到的电磁场分量影响较大,单元的误差估计值较大,所以在自适应网格优化过程中,对浅地表和测点附近的网格进行了加密,而边界和深部的网格单元误差很小,所以并没有进行加密。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 陆地起伏地形模型不通细化次数网格剖面

Fig.2 Sections of different mesh refinement times of land complex topography model

表1 陆地起伏地形模型网格自适应细化次数

Table 1 Mesh adaptive refinement of land complex topography model

细化次数	节点数	单元数
1	2137	4246
2	2977	5910
3	4493	8933
4	6644	13228
5	9526	18986
6	13448	26823
7	18815	37548
8	26028	51965
9	35471	70849

新窗口打开| 下载CSV

模型进行正演计算后,得到了TE、TM模式的视电阻率和相位值,然后添加4%的随机噪声,则生成了反演计算所需的观测数据。反演的初始模型包含空气层和均匀地层,其中空气不参与反演,给定均匀地层的初始电阻率为1.0 Ω·m。计算区域所划分的单元越多,所需要的计算量越大,因此只对目标反演区域采用较精细的网格剖分,其余部分则采用粗网格以减少不必要的计算量。如图3所示,在y范围为-10~10 km及地表到10 km深的目标区域,调用Triangle程序对该区域进行精细网格剖分,共生成19 514个三角单元,而在其余区域共剖分产生了1 130个三角形单元。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 反演初始模型及其网格剖分

Fig.3 The initial model and its grid subdivision used in the inversion

利用并行机进行反演试算,如表2所示,反演经过69次迭代,耗时7 h,占用内存约为12 GB,最终模型均方根拟合差(RMS)值为4.91,反演结果如图4所示。从图中可以看出,该算法能够准确地反演出两个低阻异常体在复杂模型中的埋深、大小和形态,同时带地形的MT反演能够较好地区分不同电性的地层,这也说明MT无论是对局部低阻异常体还是对较大尺度的电性结构都能有效地探测。

表2 反演参数

Table 2 Inversion parameters

种类	数值
迭代次数/次	69
消耗时间/h	7
RMS	4.91
占用内存/GB	12

新窗口打开| 下载CSV

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 陆地起伏地形模型反演结果

Fig.4 Inversion result of land complex topography model

图5为真实模型响应数据与反演最终模型响应数据的对比,由图5可以看出,真实模型的响应数据与反演模型得到的视电阻率和相位信息高度吻合,验证了反演结果的可靠性和准确度。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 真实模型响应(a)与反演模型响应(b)拟断面

Fig.5 Quasi-sectional view of MT response with real model (a) and inverse model (b)

2.2 海底起伏地形模型

设计了一个起伏海底模型,模型结构如图6所示,海水深度在2 242~3 012 m之间变化,设海水的电阻率为0.3 Ω·m,在不同的地层中分别设置了两个低阻异常体,1号异常体电阻率为0.6 Ω·m,2号异常体电阻率为3 Ω·m。MT接收点从左往右依次设置了21个测点每个测点间距800 m。MT的频率范围为0.001~1 Hz,共选取了21个频点,取对数等间隔分布。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 起伏海底模型

Fig.6 Undulating seabed model

对上述起伏海底模型进行自适应有限元正演计算,得到模型正演响应与自适应细化网格,表3为起伏海底模型细化次数,从表3中可以看出,相比于陆地反演模型,海底模型的剖分次数更多,产生的节点数以及单元数也比陆地模型的数量多。

表3 起伏海底模型网格自适应细化次数

Table 3 Mesh adaptive refinement times of undulating seabed model

细化次数	节点数	单元数
1	2299	4549
2	3098	6142
3	4304	8548
4	5965	11869
5	8556	17048
6	11651	23230
7	15885	31693
8	21834	43579
9	30020	59933
10	41090	82060
11	56278	112421
12	76582	153009
13	79516	158768

新窗口打开| 下载CSV

在并行机上进行反演试算,如表4所示,反演经过21次迭代,耗时91.5 min,最终模型的均方根拟合差(RMS)值为1.151 8,占用内存约为3.35 GB。反演结果如图7所示,可以看出反演结果能够准确呈现出两个异常体所处位置、体积形态等特征,并且具有较高的分辨率,验证了算法对于起伏海底模型的适用性。图8为真实模型响应数据与反演最终模型响应数据的对比,由图8可以看出,真实模型的响应数据与反演模型得到的视电阻率和相位信息高度吻合,体现了反演结果的真实性和可靠性。

表4 反演参数

Table 4 Inversion parameters

种类	数值
迭代次数/次	21
消耗时间/min	91.5
RMS	1.1518
占用内存/GB	3.35

新窗口打开| 下载CSV

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 起伏海底模型反演结果

Fig.7 Inversion result of undulating seabed model

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 真实模型响应(a)与反演模型响应(b)拟断面

Fig.8 Quasi-sectional view of MT response with real model (a) and inverse model (b)

3 实测数据反演

3.1 研究区概况

研究区域地处西伯利亚板块和哈萨克斯坦板块交界地带,其具体位置处于克拉玛依后山一带,邻近克拉玛依市区。综合该区域地质构造及前人研究资料可知,研究区在地层划分上属于北疆—兴安岭地层大区北疆地层区,位于哈图断裂以南的克拉玛依地层小区,以石炭系为主,未见泥盆系。区内的地层分区在二叠纪之前表现明显,早二叠世晚期的卡拉岗组陆相火山岩及火山碎屑岩的出现代表了地层分区性消失,进入统一的陆内演化阶段。测区内中新生界地层除了古近系外,其余均有出露。

图9为研究区地质及测线布置图,MT测点采用约2 km的间距布置12个测点,测线方向为135°,测线长度约为22 km。从图9中可以看出达尔布特断裂横穿测线而过,在测线上方则有哈图断裂横穿而过,而在测线下方含有大量二叠纪侵入体。

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 研究区地质与测线布置

Fig.9 Research area geological and line arrangement map

3.2 数据反演与解释

对实测数据进行处理,选取其中42个频点的观测数据进行反演,频点选取范围是1×10^-2~1×10² Hz。为了减少数据的噪声对反演结果带来的偏差,在反演进行之前,对实测数据进行了圆滑处理。数据在经过36次迭代,耗时118.2 min,最终得到反演结果。

图10为测线已知地质剖面情况和MT数据非线性共轭梯度反演结果的电阻率剖面,图11为测线浅部剖面Occam二维反演结果断面。从图11中可以看出,浅部剖面的反演结果具有清晰的电性特征,在剖面中段(10~14 km)存在明显的近似垂直低阻特征,为达尔布特断裂,断裂向西北方向倾,最大深度可达2 km左右,并且有向深部延伸的趋势,这与地质资料相匹配,验证了将本文算法程序应用于MT数据得到的反演结果具有较高的准确性。

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 测线地质剖面(a)和MT数据非线性共轭梯度反演结果(b)

Fig.10 Geological section (a)and nonlinear conjugate gradient inversion results (b) of MT data

图11

新窗口打开| 下载原图ZIP| 生成PPT

图11 MT数据Occam反演电阻率浅部断面

Fig.11 Shallow section of Occam inversion resistivity of MT data

图12为测线反演结果的深部断面,其中黑色倒三角为测点位,测线方向为135°,结合浅部数据,可以观察到深度2 km以下的部分在断裂的北侧有向左下方延伸的低阻带,并且在0 km附近的深部区域出现了明显的低阻异常,与资料中哈图断层相吻合。而作为一个大断裂,达尔布特断裂两侧均分布有太勒古拉组中段,在断裂带北侧电阻率偏低,南侧偏高;在达尔布特断裂南侧剖面图上(14~20 km)存在很高的电阻率特征,为克拉玛依花岗岩岩体的侵入;在达尔布特断裂北侧高阻特征不明显,这和地质图对应得很好,进一步验证了算法的准确性。

图12

新窗口打开| 下载原图ZIP| 生成PPT

图12 MT数据Occam反演电阻率深部断面

Fig.12 Deep section of Occam inversion resistivity of MT data

4 结论

对基于非结构自适应有限元的带地形二维MT反演算法进行了应用研究,通过二维带地形陆地及海底模型反演计算,表明该算法可以有效压制静态效应,反演对TE数据的拟合要比对TM数据的拟合更好,同时带地形模型反演对异常体反演效果较好,对异常体的大小和形态还原度很高,验证了该算法对于陆地及海底MT数据的共同实用性,证明带地形MT反演能够准确地揭示地下电性结构和探明低阻异常体。

同时,应用该方法对克拉玛依后山区域实测数据进行反演,最终得到电阻率结构清晰的反演结果,并将反演结果与已知地质资料以及非线性共轭梯度法的反演结果进行了对比验证,能够清楚显示出地质资料中表明的达尔布特断裂等地质构造。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Huang

X Y

, Deng

J Z

, Chen

, et al.

Magnetotelluric extremum boundary inversion based on different stabilizers and its application in a high radioactive waste repository site selection

[J]. Applied Geophysics, 2019, 16(3) : 367-377.

DOI:10.1007/s11770-019-0755-4 URL [本文引用: 1]

[2]

李磊

湘南骑田岭锡铅锌多金属矿区岩矿石电性研究

[J]. 物探与化探, 2007(S1) : 77-80,93.