基于磁异常ΔT相关成像约束的非结构网格磁化强度反演

杨景旭; 安国强; 李柏森; 鲁宝亮

doi:10.12404/j.issn.1671-1815.2502820

科学技术与工程 ›› 2026, Vol. 26 ›› Issue (9) :3686-3695 DOI: 10.12404/j.issn.1671-1815.2502820

天文学、地球科学

基于磁异常ΔT相关成像约束的非结构网格磁化强度反演

杨景旭 ¹ ,
安国强 ² ,
李柏森 ¹ ,
鲁宝亮 ¹^,³^,⁴^,^*

作者信息 +

Unstructured Grid Magnetization Inversion Based on Magnetic Anomaly ΔT Correlation Imaging Constraints

YANG Jing-xu ¹ ,
AN Guo-qiang ² ,
LI Bo-sen ¹ ,
LU Bao-liang ¹^,³^,⁴^,^*

Author information +

文章历史 +

摘要

针对复杂地形与地质体的情况下地下三维磁性结构分布高精度反演难题,提出基于磁异常ΔT相关成像约束的非结构网格磁化强度反演方法。采用四面体非结构网格以提升复杂地形地质拟合精度;借助相关成像识别重点区域进而实现模型降维,大幅压缩解空间规模,提升反演效率;对降维后的重点区域,引入深度加权约束与物性范围约束进行正则化反演,能有效聚焦目标区域,提高反演精度;运用预条件共轭梯度法求解目标函数,并引入OpenMP并行算法,进一步提升反演效率。通过模型试验验证了该方法在提高磁化强度反演精度和效率方面的优势,将其应用于加拿大安大略省McFaulds Lake地区实测航空磁异常数据,也进一步证实了方法的正确性与有效性。该方法为复杂地质体的磁化强度反演提供了全新技术手段,具有显著的创新价值与应用潜力。

Abstract

To address the challenge of high-precision inversion for 3D subsurface magnetic structure distribution under complex terrain and geological body conditions, an unstructured grid magnetization inversion method constrained by magnetic anomaly ΔT correlation imaging was proposed. Tetrahedral unstructured grids were adopted to enhance the fitting accuracy of complex topography and geology. Correlation imaging was utilized to identify key regions, thereby enabling model dimensionality reduction, significantly compressing the solution space scale, and improving inversion efficiency. For the key regions after dimensionality reduction, depth-weighted constraints and physical property range constraints are introduced for regularized inversion, which can effectively focus on the target regions and improve the inversion accuracy. The preconditioned conjugate gradient method was adopted to solve the objective function, and the OpenMP parallel algorithm was introduced to further boost the inversion efficiency. The advantages of this method in improving the accuracy and efficiency of magnetization inversion were verified by model tests, and when it is applied to the measured aeromagnetic anomaly data in the McFaulds Lake area, Ontario, Canada, the correctness and effectiveness of the method are further confirmed. This method provides a new technical means for the magnetization inversion of complex geological bodies and possesses significant innovative value and application potential.

Graphical abstract

关键词

磁异常 / 非结构网格 / 相关成像 / 模型降维 / 三维磁化强度反演

Key words

magnetic anomalies / unstructured grid / correlation imaging / model dimension reduction / 3D magnetization inversion

引用本文

引用格式 ▾

杨景旭, 安国强, 李柏森, 鲁宝亮. 基于磁异常ΔT相关成像约束的非结构网格磁化强度反演[J]. 科学技术与工程, 2026, 26(9): 3686-3695 DOI:10.12404/j.issn.1671-1815.2502820

YANG Jing-xu, AN Guo-qiang, LI Bo-sen, LU Bao-liang. Unstructured Grid Magnetization Inversion Based on Magnetic Anomaly ΔT Correlation Imaging Constraints[J]. Science Technology and Engineering, 2026, 26(9): 3686-3695 DOI:10.12404/j.issn.1671-1815.2502820

登录浏览全文

4963

注册一个新账户忘记密码

磁法勘探是地球物理勘探中最重要的勘探手段之一,其核心在于通过观测地表的磁异常来推断地下地质结构和矿产资源的分布,具有勘探深度大、成本低、效率高的优点,被广泛应用于构造单元划分、矿产资源勘探等领域。三维磁化强度反演是揭示地球内部物质结构最常用的方法之一,该方法通过将观测区域地下空间离散化为大量的网格单元,经过反演得到每个网格单元的物性值,进而确定不均匀磁性体的空间位置与物性分布,准确获取复杂地形和地质体条件下的地下三维磁性结构分布。

在三维磁化强度反演中,结构化直立六面体网格剖分是使用最广泛的一种方式^[1-3],这种网格剖分方法简单,但是对于起伏地形和复杂的地质体边界的拟合效果较差,直接影响反演的精度。非结构化四面体网格的剖分方法近年来受到了广泛的关注^[4-8]。Okabe^[9]通过应用散度定理和格林公式,推导了在任意方向上由多边形面组成的均匀多面体重磁异常的完整解析表达式,并分析了可能存在的奇点问题。非结构化四面体网格可以更精细地拟合起伏地形和复杂的地质体边界,从而可以获得更加准确的反演结果,然而,相比于直立六面体网格,四面体网格的正演计算公式更加复杂,计算量更大。

由于地球物理场的等效性以及观测数据的有限性,地球物理反演面临着多解性和不稳定性的问题,提升反演算法的精度与效率已成为近年来的研究热点。目前,Tikhonov正则化作为一种有效减少多解性的方法,被广泛应用于地球物理反演中^[1,7,10-13]。为了更准确地反映地下真实异常体的位置与物性分布,前人提出了多种约束条件,主要包括深度加权约束^[1,11-14]、光滑约束^[1]、聚焦约束^[15-16]、物性范围约束^[17-18]以及参考模型约束^[19]等。然而在地球物理反演中,地下剖分单元过多会导致矩阵过大,进而引发内存不足、运算时间过长等问题,而通过减少模型数量以及采用并行计算的方式,可较为有效地解决这些问题。例如,借助边缘识别方法选取重点反演区域,进而实现模型降维,降低反演的多解性^[8,20-22]。

在前人的研究基础上,针对复杂地形和地质体条件下地下三维磁性结构分布问题,现提出一种基于磁异常ΔT相关成像约束的非结构网格磁化强度反演方法,通过使用四面体非结构网格来拟合复杂地质体,利用相关系数进行模型降维,并使用正则化方法对降维后的重点区域进行反演,引入深度加权约束和物性范围约束提高反演的精度,最后利用预条件共轭梯度法对目标函数进行求解,同时为了提高反演效率引入OpenMP并行算法,得到高精度的磁化强度反演结果。最后,通过模型试验和加拿大安大略省McFaulds Lake地区实测航空磁异常数据反演验证所提方法的正确性和有效性。

1 方法原理

1.1 非结构网格磁异常正演

采用基于Delaunay三角剖分的非结构网格剖分^[23]。Delaunay三角剖分依据的准则是空外接圆(球)准则,即在单个四面体的外接球内不包含除本身的第5个点。为了能简单、快速地将复杂的模型剖分为规则的四面体,采用开源软件TetGen,它是一款优良的四面体网格生成软件,可以对复杂的形体或地形进行剖分^[24]。TetGen剖分的输出文件中每个四面体的三角面的顶点坐标在数组中的顺序不满足正演理论的排序要求,会影响到线积分的结果。因此,需要对顶点进行排序,使得从外表面方向观察,顶点坐标在数组中必须按照逆时针排序,且各个面的排序方向相同。

对一个任意的四面体ABCD(图1),Okabe^[9]给出了以三角面逼近任意多面体的重力梯度计算方法,其基本思想是根据散度定理将多面体的体积分转化为面积分,再根据格林公式转化为线积分,从而完成重力梯度的计算。即分别计算4个三角面的面积分,将面积分又转化为线积分,则完成了一个四面体的体积分计算。

在直角坐标系中,地下地质体$Q(\xi,\eta,\zeta )$在计算点P(x, y, z)产生的重力位V的计算公式为

(1)$V=G\underset{V}{\int \int \int }\frac{\rho }{[{(\xi -x)}^{2}+{(\eta -y)}^{2}+{(\zeta -z)}^{2}{]}^{1/2}}\mathrm{d}V$

式(1)中:G为万有引力常数;ρ为场源密度;(x, y, z)为观测点坐标;$(\xi,\eta,\zeta )$为场源坐标。

对于将任意形体用非结构网格四面体进行剖分,计算所有四面体在场点产生的重力梯度值然后求和,即得到任意形体的重力梯度。一个四面体在计算点(x, y, z)产生的沿kl方向的重力梯度V_kl为

(2)${V}_{kl}=G\rho \stackrel{4}{\sum _{i}}\stackrel{3}{\sum _{j}}{l}^{\mathrm{T}}{n}_{i}{L}_{j}\left(i\right)$

式(2)中:l^T=(l_x,l_y,l_z)为l方向的方向向量;n_i为第i个三角形的外法向量;L_j(i)为四面体第i个三角形第j边的积分值。

(3)$\begin{array}{l}{L}_{j}\left(i\right)=\left\{\right({M}_{y}\mathrm{c}\mathrm{o}\mathrm{s}\psi -{M}_{x}\mathrm{s}\mathrm{i}\mathrm{n}\psi \left)\mathrm{l}\mathrm{n}\right[\xi +\\ ({\xi }^{2}+{\eta }^{2}+{Z}^{2}{)}^{1/2}]+\\ {M}_{z}\mathrm{t}\mathrm{a}{\mathrm{n}}^{-1}\frac{-\xi \eta +({\eta }^{2}+{Z}^{2})\mathrm{t}\mathrm{a}\mathrm{n}\psi }{Z({\xi }^{2}+{\eta }^{2}+{Z}^{2}{)}^{1/2}}\}{|}_{{\xi }_{i}}^{{\xi }_{i+1}}\end{array}$

式(3)中:ψ为旋转角度;Z为旋转后的z坐标。M_x、M_y、M_z的表达式如下。

$\begin{array}{l}{M}_{x}={k}_{x}\mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{c}\mathrm{o}\mathrm{s}\varphi +{k}_{y}\mathrm{s}\mathrm{i}\mathrm{n}\theta \mathrm{c}\mathrm{o}\mathrm{s}\varphi -{k}_{z}\mathrm{s}\mathrm{i}\mathrm{n}\varphi ;\\ {M}_{y}=-{k}_{x}\mathrm{s}\mathrm{i}\mathrm{n}\theta +{k}_{y}\mathrm{c}\mathrm{o}\mathrm{s}\theta ;\\ {M}_{z}={k}_{x}\mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{s}\mathrm{i}\mathrm{n}\varphi +{k}_{y}\mathrm{s}\mathrm{i}\mathrm{n}\theta \mathrm{s}\mathrm{i}\mathrm{n}\varphi +{k}_{z}\mathrm{c}\mathrm{o}\mathrm{s}\varphi 。\end{array}$

根据重磁位泊松公式即可正演出磁异常,表达式为

(4)$\begin{array}{l}\mathrm{\Delta }T={H}_{\mathrm{a}x}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{t}}+{H}_{\mathrm{a}y}\mathrm{c}\mathrm{o}\mathrm{s}{\beta }_{\mathrm{t}}+{Z}_{\mathrm{a}}cos{\gamma }_{\mathrm{t}}\\ =\frac{{\mu }_{0}M}{4\mathrm{\pi }G\rho }\left[\right({V}_{xx}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{s}}+{V}_{xy}\mathrm{c}\mathrm{o}\mathrm{s}{\beta }_{\mathrm{s}}+{V}_{xz}\mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_{\mathrm{s}})\times \\ \mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{t}}+({V}_{xy}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{s}}+{V}_{yy}\mathrm{c}\mathrm{o}\mathrm{s}{\beta }_{\mathrm{s}}+{V}_{yz}\mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_{\mathrm{s}})\times \\ \mathrm{c}\mathrm{o}\mathrm{s}{\beta }_{\mathrm{t}}+({V}_{xz}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{s}}+{V}_{yz}\mathrm{c}\mathrm{o}\mathrm{s}{\beta }_{\mathrm{s}}+{V}_{zz}\mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_{\mathrm{s}})\times \\ \mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_{\mathrm{t}}]\end{array}$

式(4)中:cosα_s、cosβ_s、cosγ_s为磁化强度的方向余弦;cosα_t、cosβ_t、cosγ_t为正常地磁场的方向余弦;H_ax、H_ay、Z_a分别为水平磁异常x分量、y分量和垂直磁异常;μ₀为真空中的磁导率;M为磁化强度;V_kl为kl方向的重力梯度,$k,l\in (x,y,z)$。

1.2 磁异常ΔT相关成像

磁异常相关成像作为一种半定量的方法, 通过计算地下每个磁偶极子在地面产生的磁异常与实测磁异常之间的归一化互相关系数, 用来表示各个位置的概率分布。

定义测网内测点的实测磁异常ΔT与在地下三维空间内第q个磁偶极子的磁异常ΔT_q的互相关系数^[25]为

(5)${C}_{q}=\frac{\stackrel{{N}_{\mathrm{s}}}{\sum _{i=1}}\mathrm{\Delta }T({x}_{i},{y}_{i},{z}_{i}){B}_{q}({x}_{i},{y}_{i},{z}_{i})}{\sqrt{\stackrel{{N}_{\mathrm{s}}}{\sum _{i=1}}\mathrm{\Delta }{T}^{2}({x}_{i},{y}_{i},{z}_{i})\stackrel{{N}_{\mathrm{s}}}{\sum _{i=1}}{B}_{q}^{2}({x}_{i},{y}_{i},{z}_{i})}}$

式(5)中:N_s为测区观测点总数;B_q(x_i,y_i,z_i)为地下三维测区内第q个磁偶极子的基函数,表达式为

(6)$\begin{array}{l}{B}_{q}({x}_{i},{y}_{i},{z}_{i})=\frac{1}{{r}^{5}}\left[A\right({x}_{q}-{x}_{i}{)}^{2}+B({y}_{q}-{y}_{i}{)}^{2}+\\ C({x}_{q}-{x}_{i})({z}_{q}-{z}_{i})+\\ D({y}_{q}-{y}_{i})({z}_{q}-{z}_{i})+\\ E({x}_{q}-{x}_{i})({y}_{q}-{y}_{i})+\\ F({z}_{q}-{z}_{i}{)}^{2}]\end{array}$

式(6)中:各物理量的表达式如下。

$\begin{array}{l}r=\sqrt{({x}_{q}{-x)}^{2}+({y}_{q}{-y)}^{2}+({z}_{q}{-z)}^{2}};\\ l=\mathrm{c}\mathrm{o}\mathrm{s}{I}_{0}\mathrm{c}\mathrm{o}\mathrm{s}{A\text{'}}_{0};m=\mathrm{c}\mathrm{o}\mathrm{s}{I}_{0}\mathrm{s}\mathrm{i}\mathrm{n}{A\text{'}}_{0};n=\mathrm{s}\mathrm{i}\mathrm{n}{I}_{0};\\ L=\mathrm{c}\mathrm{o}\mathrm{s}I\mathrm{c}\mathrm{o}\mathrm{s}A\text{'};M=\mathrm{c}\mathrm{o}\mathrm{s}I\mathrm{s}\mathrm{i}\mathrm{n}A\text{'};N=\mathrm{s}\mathrm{i}\mathrm{n}I\text{'};\\ A=2lL-nN-mM;\\ B=2mM-lL-nN;\\ C=3(lN+nL);\\ D=3(mN+nM);\\ E=3(lM+mL);\\ F=2nN-mM-lL。\end{array}$

相关系数C_q表明了地下某个磁偶极子引起异常的概率大小,取值范围为[0,1],其值的意义是实测磁异常ΔT是由第q个磁偶极子产生的可能性。C_q越趋向于1,表明该磁偶极子磁性高的可能性越大;反之,其值越接近于0,表明该磁偶极子磁性低的可能性越大。

在计算所有网格单元的相关系数后,选择相关系数较大的网格单元进行反演,即可达到模型降维的目的。经过模型降维处理,将圈定的异常源区域作为靶区进行反演,认为其他区域不存在异常源从而不参与反演,可以有效地解决运算时间长、内存不足的问题。

1.3 正则化反演

基于Tikhonov正则化反演方法,构建目标函数为

(7)$\varphi ={\varphi }_{d}+\lambda {\varphi }_{\mathrm{m}}$

式(7)中:${\varphi }_{\mathrm{d}}=\Vert {W}_{\mathrm{d}}(Gm-d){\Vert }^{2}$为数据约束项,W_d为数据加权对角矩阵,G为核函数矩阵,m为模型参数向量,d为观测数据向量; ${\varphi }_{\mathrm{m}}=\Vert {W}_{\mathrm{m}}m{\Vert }^{2}$为模型约束项,W_m为模型加权对角矩阵; λ为正则化参数,控制着数据拟合以及模型约束项在反演计算过程中的相对重要程度。一般使用多次试验、广义交叉验证法^[11,26]或L曲线法^[27]对正则化参数进行取值。本文研究中,正则化参数根据多次试验的反演效果确定。

在正则化反演过程中,由于数据约束项与模型约束项两者的物理量纲不一致,正则化参数λ很难确定,故采用最大幅值差进行归一化处理。

(8)${\gamma }_{d}={d}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{d}_{\mathrm{m}\mathrm{i}\mathrm{n}}$

(9)${\gamma }_{m}={m}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{m}_{\mathrm{m}\mathrm{i}\mathrm{n}}$

式中:γ_d为数据约束项最大幅值差归一化因子;d_max、d_min分别为观测数据的最大、最小值;γ_m为模型约束项最大幅值差归一化因子;m_max、m_min分别为磁化强度的最大、最小值。此时目标函数为

(10)$\varphi ={‖{W}_{\mathrm{d}}(Gm-d)/{\mathrm{\gamma }}_{\mathrm{d}}‖}^{2}+\lambda {‖{W}_{\mathrm{m}}m/{\mathrm{\gamma }}_{\mathrm{m}}‖}^{2}$

磁数据反演中,核函数G随深度的增大迅速衰减,使得反演的结果集中于地表附近,不能够准确地反映地下空间真实的物性分布状态。为了克服上述问题,利用Li等 ^[1]提出的深度加权函数,其表达式为

(11)$w\left(z\right)=(z+{z}_{0}{)}^{-\beta /2}$

式(11)中:z为剖分单元的中心点深度,取四面体各顶点深度的平均值; z₀和β为常数,z₀取决于观测面的高度,β为深度加权因子,根据磁场随深度的衰减规律,一般在磁反演中β取3.0,将深度加权约束矩阵W_z加入目标函数后形式为

(12)$\varphi ={‖{W}_{\mathrm{d}}(Gm-d)/{\gamma }_{d}‖}^{2}+\lambda {‖{W}_{\mathrm{m}}{W}_{\mathrm{z}}m/{\gamma }_{m}‖}^{2}$

对应的求解方程为

(13)$\begin{array}{l}\left\{\right({W}_{\mathrm{d}}{G)}^{\mathrm{T}}\left({W}_{\mathrm{d}}G\right)/{\gamma }_{d}^{2}+\lambda ({W}_{\mathrm{m}}{W}_{\mathrm{z}}{)}^{\mathrm{T}}({W}_{\mathrm{m}}{W}_{\mathrm{z}})/{\gamma }_{m}^{2}\}m=\\ \left({W}_{\mathrm{d}}G\right)T{W}_{\mathrm{d}}d/{\gamma }_{d}^{2}\end{array}$

在求解方程的过程中,反演得到的物性结果很可能不具有实际物理意义,给后续解释带来困难,可以在反演过程中引入物性范围对反演解进行物理意义上的约束。主要的物性范围约束方式包括转换函数法^[28]、对数转换法^[1]、罚函数法^[11]。本文研究使用强制约束的方法,将反演磁化强度值约束在给定的范围。

共轭梯度法具有计算速度快、内存需求少、能够避免大型矩阵的乘积计算的特点,是求解线性方程组的有效方法之一^[29-31]。为了改善迭代法的收敛性,增强解的稳定性,选用预条件共轭梯度法(PCG)进行求解,将原方程Ax=b替换为SAx=Sb,用深度加权函数即式(11)作为预条件因子S。

另外,为了进一步提高反演的效率,采用并行方法进行计算加速。OpenMP是共享存储的并行编程模型,具有可移植和易编程性等特点,适用于循环级并行等任务,对程序代码改动小,在科学计算上占据优势。因此采用OpenMP并行计算方法来提升三维磁化强度的反演效率。

模型试验和实际资料处理均在工作站上完成,工作站的配置为:2 颗INTEL(R) XEON(R) PLATINUM 8575C 2.80 GHz CPU,每个CPU 48个内核,576 GB内存。

2 理论模型试验

2.1 模型试验1(单个模型)

构建单个直立长方体模型,模型参数如表1所示。x方向的范围为0~1 000 m,y方向为0~800 m,测点间距为25 m,平面测网高度-10 m,反演深度设定为0~800 m。该模型的位置示意图如图2(a)所示。在不考虑剩磁的情况下,对模型进行正演计算,获得其磁异常数据如图2(b)所示。

首先对整个研究区域进行非结构网格剖分,共剖分为30 495个单元。将正演得到的磁异常数据作为观测数据,开展反演模型试验。其中选择正则化参数λ为50,深度加权参数z₀为10 m,不进行模型降维的三维反演结果如图3(a)所示,反演结果无法准确反映模型的真实位置和形状,磁性体的分布较为分散;重建的磁异常、残差如图3(b)、图3(c)所示,均方根误差为4.49 nT,表明反演结果与观测数据的拟合效果较差。计算效率方面,反演的迭代次数为10 000次,所用时间为10 026.72 s,表明不降维反演结果的精度和效率均较低。

为了提高反演的精度和效率,提出了基于磁异常ΔT相关成像约束的非结构网格磁化强度反演方法。即通过分析相关成像的结果进行模型降维,选择异常体所在的重点区域参与反演,降维反演的结果如图4(a)所示,准确地反映了模型的真实位置和形状。数据重建的结果、重建残差结果如图4(b)、图4(c)所示,均方根误差为0.62 nT,表明反演结果与观测数据的拟合效果较好,反演精度较高。同时,降维后参与反演的网格单元为2 211个,迭代次数为1 284次,反演时间为6.69 s,表明降维反演提高了反演的效率,同时获得了更优的反演结果。

2.2 模型试验2(组合模型)

为了进一步验证所提方法的有效性,构建了一个包含两个直立长方体的组合模型,其具体参数如表2所示。观测网和地形剖分与模型试验1相同。该模型的位置示意图如图5(a)所示。在不考虑剩磁的情况下,对模型进行正演计算,获得其磁异常数据如图5(b)所示。

将正演得到的磁异常数据作为观测数据,开展反演模型试验。其中选择正则化参数λ为50,深度加权参数z₀为10 m,不进行模型降维的三维反演结果如图6(a)所示,反演结果无法准确反映模型的真实位置和形状,磁性体的分布较为分散且浮于地表,与实际模型的两个直立长方体的分布存在较大差异;重建的磁异常、残差如图6(b)、图6(c)所示,均方根误差为1.19 nT,表明反演结果与观测数据的拟合效果较差。计算效率方面,反演的迭代次数为10 000次,所用时间为9 478.41 s,表明不降维反演结果的精度和效率均较低。

通过分析相关成像的结果进行模型降维,选择异常体所在的重点区域参与反演,降维反演的结果如图7(a)所示,准确地反映了模型的真实位置和形状。数据重建的结果、重建残差结果如图7(b)、图7(c)所示,均方根误差为0.67 nT,表明反演结果与观测数据的拟合效果较好,反演精度较高。同时,降维后参与反演的网格单元为4 409个,迭代次数为1 466次,反演时间为55.67 s,表明降维反演提高了反演的效率,同时获得了更优的反演结果。

3 实际资料处理

McFaulds Lake位于加拿大安大略省西北部,在此地发现了鹰巢镍、铜和铂族元素(Ni-Cu-PGE)矿床。该区域属于源自地幔的、具有高磁性的超镁铁质侵入体的一部分。由于该侵入体在航磁数据上形成了一个直径约60 km的半圆环形强磁异常,在勘探早期被称为“火环”。图8展示了“火环地区”已知矿化情况的地质简图^[32]。该地区已勘探出多个经济矿床,包括:岩浆型Ni-Cu-PGE矿床、岩浆型铬铁矿矿床、火山块状硫化物矿床以及由金伯利岩承载的钻石矿床。“火环”底部的主要层状侵入岩体中,蕴藏着品位极高的Ni-Cu-PGE矿床,而在其层状侵入岩体地层结构更靠东、更高处则有层状铬铁矿矿床^[33-34]。铬铁矿床通常出现在层状超镁铁质侵入岩中,常与磁铁矿和蛇纹石共生。因此,人们可以通过磁性追踪铬铁矿的可能分布。2010年在McFaulds Lake地区开展了航空地球物理调查,采集了磁力数据^[33]。

如图8所示,选取了存在3个铬铁矿床Big Daddy、Black Thor、Black Label的矩形范围为反演目标区域,目标区域的长×宽×深为8 km×8 km×2 km,选取大地水准面为垂直方向零点,将起伏地形高程数据[图9(a)]组合为反演区域,并对该区域进行四面体剖分,剖分数量为99 886个,观测面x方向范围为548~556 km,y方向范围为5 843~5 851 km,点距为50 m,总观测点数25 921个。

目标区域总磁异常ΔT数据如图9(b)所示,通过相关成像选取重点区域进行降维,降维后的模型数量为82 341个,对矿区进行磁化强度反演,深度加权的参数z₀为100 m,参数β为3.0,物性约束范围为0~20 A/m。

相关成像降维的三维磁化强度反演结果(图10)显示了McFaulds Lake地区存在的高磁性岩体。在x=552 km[图10(a)]、y=5 846 km[图10(b)]处的垂直切片也显示出同样的高磁性岩体分布特征和z=0.35 km[图10(c)]处的水平切片,在垂直切片中,高磁性区域主要集中在深度为300~700 m,呈现出条带状分布,水平切片图进一步表明高磁性岩体的深度范围和分布形态。图10(d) 、图10(e)分别为重建的磁异常图和重建残差图,重建磁异常的均方根误差仅为1.94 nT,表明反演模型能够有效拟合实际观测数据,验证了反演结果的可靠性。与前人的反演结果^[34-35]比较,形态基本吻合,进一步验证了本文方法的正确性和有效性。

4 结论

针对复杂地形和地质体条件下地下三维磁性结构分布问题,提出了一种基于磁异常相关成像约束的非结构网格磁化强度反演方法,得到以下结论。

(1)模型降维显著提高了反演精度和效率。通过相关成像进行模型降维,有效减少了反演计算中的网格单元数量,提高了反演结果的精度和效率。模型试验和实际资料处理结果进一步验证了本文方法的正确性和有效性。

(2)非结构网格剖分能够更好地拟合复杂地质体边界。相比于传统的结构化网格,非结构四面体网格能够更精细地拟合起伏地形和复杂的地质体边界,从而获得更加精细和准确的反演结果。这为复杂地质条件下的磁化强度反演提供了一种有效的技术手段。

综上所述,本文方法在提高反演精度、拟合复杂地质体以及提升计算效率方面均表现出显著优势,为复杂地质体的三维磁化强度反演提供了一种可靠的技术手段,具有广阔的应用前景。

参考文献

原文顺序 | 出版日期 | 本文引用

[1]	Li Y G, Oldenburg D W. 3-D inversion of magnetic data[J]. Geophysics, 1996, 61(2): 394-408.

[2]	孟小红, 刘国峰, 陈召曦, 等. 基于剩余异常相关成像的重磁物性反演方法[J]. 地球物理学报, 2012, 55(1): 304-309.

[3]	Meng Xiaohong, Liu Guofeng, Chen Zhaoxi, et al. 3-D gravity and magnetic inversion for physical properties based on residual anomaly correlation[J]. Chinese Journal of Geophysics, 2012, 55(1): 304-309.

[4]	高秀鹤, 黄大年. 基于共轭梯度算法的重力梯度数据三维聚焦反演研究[J]. 地球物理学报, 2017, 60(4): 1571-1583.

[5]	Gao Xiuhe, Huang Danian. Research on 3D focusing inversion of gravity gradient tensor data based on a conjugate gradient algorithm[J]. Chinese Journal of Geophysics, 2017, 60(4): 1571-1583.

[6]	曹书锦, 朱自强, 鲁光银, 等. 基于单元细分H-自适应有限元全张量重力梯度正演[J]. 地球物理学进展, 2010, 25(3): 1015-1023.

[7]	Cao Shujin, Zhu Ziqiang, Lu Guangyin, et al. Forward modelling of full gravity gradient tensors based H-adaptive mesh refinement[J]. Progress in Geophysics, 2010, 25(3): 1015-1023.

[8]	Jahandari H, Farquharson C G. Forward modeling of gravity data using finite-volume and finite-element methods on unstructured grids[J]. Geophysics, 2013, 78(3): 69-80.

[9]	Ren Z Y, Tang J T, Kalscheuer T, et al. Fast 3-D large-scale gravity and magnetic modeling using unstructured grids and an adaptive multilevel fast multipole method[J]. Journal of Geophysical Research: Solid Earth, 2017, 122(1): 79-109.

[10]	Ma G Q, Niu R X, Gao T, et al. High-efficiency gravity data inversion method based on locally adaptive unstructured meshing[J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60: 5910507.

[11]	高新宇, 鲁宝亮, 安国强, 等. 起伏地形下基于非结构网格模型降维的重力三维密度反演[J]. 地球物理学报, 2024, 67(1): 379-393.

[12]	Gao Xinyu, Lu Baoliang, An Guoqiang, et al. 3D density inversion of gravity based on unstructured grid and model dimension reduction under undulating terrain[J]. Chinese Journal of Geophysics, 2024, 67(1): 379-393.

[13]	Okabe M. Analytical expressions for gravity anomalies due to homogeneous polyhedral bodies and translations into magnetic anomalies[J]. Geophysics, 1979, 44(4): 730-741.

[14]	Li Y G, Oldenburg D W. 3-D inversion of gravity data[J]. Geophysics, 1998, 63(1): 109-119.

[15]	Li Y G, Oldenburg D W. Fast inversion of large-scale magnetic data using wavelet transforms and a logarithmic barrier method[J]. Geophysical Journal International, 2003, 152(2): 251-265.

[16]	秦朋波, 黄大年. 重力和重力梯度数据联合聚焦反演方法[J]. 地球物理学报, 2016, 59(6): 2203-2224.

[17]	Qin Pengbo, Huang Danian. Integrated gravity and gravity gradient data focusing inversion[J]. Chinese Journal of Geophysics, 2016, 59(6): 2203-2224.

[18]	李泽林, 姚长利, 郑元满. 基于L_p范数稀疏优化算法的重力三维反演[J]. 地球物理学报, 2019, 62(10): 3699-3709.

[19]	Li Zelin, Yao Changli, Zheng Yuanman. 3D inversion of gravity data using L_p-norm sparse optimization[J]. Chinese Journal of Geophysics, 2019, 62(10): 3699-3709.

[20]	朱杰君, 姜登峰, 周文纳. 基于深度加权的重力异常广义线性反演成像方法研究及应用[J]. 科学技术与工程, 2018, 18(24): 22-27.

[21]	Zhu Jiejun, Jiang Dengfeng, Zhou Wenna. The imaging method of generalized linear inversion for gravity data and its application based on depth weighting function[J]. Science Technology and Engineering, 2018, 18(24): 22-27.

[22]	Zhdanov M S, Ellis R, Mukherjee S. Three-dimensional regularized focusing inversion of gravity gradient tensor component data[J]. Geophysics, 2004, 69(4): 925-937.

[23]	鄢文强, 程天君, 陈晓, 等. 基于指数型聚焦稳定器和惩罚函数约束的重力三维反演[J]. 科学技术与工程, 2025, 25(3): 893-899.

[24]	Yan Wenqiang, Cheng Tianjun, Chen Xiao, et al. Three-dimensional gravity inversion based on exponential focus stabilizer and penalty function constraints[J]. Science Technology and Engineering, 2025, 25(3): 893-899.

[25]	Portniaguine O, Zhdanov M S. Focusing geophysical inversion images[J]. Geophysics, 1999, 64(3): 874-887.

[26]	Portniaguine O, Zhdanov M S. 3-D magnetic inversion with data compression and image focusing[J]. Geophysics, 2002, 67(5): 1532-1541.

[27]	陈石, 王青华, 王谦身, 等. 云南鲁甸M_S6.5地震震源区和周边三维密度结构及重力场变化[J]. 地球物理学报, 2014, 57(9): 3080-3090.

[28]	Chen Shi, Wang Qinghua, Wang Qiansheng, et al. The 3D density structure and gravity change of Ludian M_S6.5 Yunnan epicenter and surrounding regions[J]. Chinese Journal of Geophysics, 2014 57(9): 3080-3090.

[29]	Yang M, Wang W Y, Welford J K, et al. 3D gravity inversion with optimized mesh based on edge and center anomaly detection[J]. Geophysics, 2019, 84(3): G13-G23.

[30]	安国强, 鲁宝亮, 高新宇, 等. 边缘特征和深度加权约束的重力三维相关成像反演[J]. 物探与化探, 2024, 48(1): 113-124.

[31]	An Guoqiang, Lu Baoliang, Gao Xinyu, et al. Correlation tomography inversion of gravity anomalies constrained by edge features and depth weighting[J]. Geophysical and Geochemical Exploration, 2024, 48(1): 113-124.

[32]	安国强, 鲁宝亮, 李柏森, 等. 边缘特征和相关成像联合约束自适应网格剖分的重力三维反演方法[J]. 地球物理学报, 2025, 68(3): 1102-1115.

[33]	An Guoqiang, Lu Baoliang, Li Bosen, et al. 3D gravity inversion based on adaptive meshing constrained by edge features and correlation imaging[J]. Chinese Journal of Geophysics, 2025, 68(3): 1102-1115.

[34]	Lawson C L. Properties of n-dimensional triangulations[J]. Computer Aided Geometric Design, 1986, 3(4): 231-246.

[35]	Hang S. TetGen, a Delaunay based quality tetrahedral mesh generator[J]. ACM Transactions on Mathematical Software, 2015, 41(2): 1-36.

[36]	郭良辉, 孟小红, 石磊. 磁异常ΔT三维相关成像[J]. 地球物理学报, 2010, 53(2): 435-441.

[37]	Guo Lianghui, Meng Xiaohong, Shi Lei. 3D correlation imaging for magnetic anomaly ΔT data[J]. Chinese Journal of Geophysics, 2010, 53(2): 435-441.

[38]	Golub G H, Heath M, Wahba G. Generalized cross-validation as a method for choosing a good ridge parameter[J]. Technometrics, 1979, 21(2): 215-223.

[39]	Hansen P C. Analysis of discrete ill-posed problems by means of the L-curve[J]. SIAM Review, 1992, 34(4): 561-580.

[40]	Commer M. Three-dimensional gravity modelling and focusing inversion using rectangular meshes[J]. Geophysical Prospecting, 2011, 59(5): 966-979.

[41]	Wilson G, Čuma M, Zhdanov M S. Massively parallel 3D inversion of gravity and gravity gradiometry data[J]. Preview, 2011, 2011(152): 29-34.

[42]	陈少华, 朱自强, 鲁光银, 等. 重力梯度张量的预条件共轭梯度法反演[J]. 中南大学学报(自然科学版), 2013, 44(2): 619-625.

[43]	Chen Shaohua, Zhu Ziqiang, Lu Guangyin, et al. Inversion of gravity gradient tensor based on preconditioned conjugate gradient[J]. Journal of Central South University (Science and Technology), 2013, 44(2): 619-625.

[44]	Hou Z L, Wei X H, Huang D N, et al. Full tensor gravity gradiometry data inversion: performance analysis of parallel computing algorithms[J]. Applied Geophysics, 2015, 12(3): 292-302.

[45]	Mungall J E, Harvey J D, Balch S J, et al. Eagle’s nest: a magmatic Ni-sulfide deposit in the James Bay Lowlands, Ontario, Canada[J]. Society of Economic Geologists Special Publication, 2010, 15: 539-557.

[46]

Ontario Geological Survey and Geological Survey of Canada. Ontario airborne geophysical surveys, gravity gradiometer and magnetic data, grid and profile data (ASCII and Geosoft R􀳳 formats) and vector data, McFaulds Lake area: Geophysical Data Set 1068[R]. Ontario: Ontario Geological Survey and Geological Survey of Canada, 2011.

[47]	Zhu Y, Zhdanov M S, Cuma M. Gramian constraints in the joint inversion of airborne gravity gradiometry and magnetic data[C]// SEG Annual Meeting. Houston, Texas: SEG, 2013: 1166-1170.

[48]	Wei L, Zhdanov M S. Joint multinary inversion of gravity and magnetic data using Gramian constraints[J]. Geophysical Journal International, 2018, 215(3): 1540-1557.