基于多重水平集的三维磁性体边界反演

肖 晓 段雅婷 胡双贵 汤井田 谢 勇 刘长生

(①中南大学有色金属成矿预测与地质环境监测教育部重点实验室，湖南长沙 410083； ②有色资源与地质灾害探查湖南省重点实验室，湖南长沙 410083； ③中南大学地球科学与信息物理学院，湖南长沙 410083； ④自然资源部覆盖区深部资源勘查工程技术创新中心，湖南长沙 410083； ⑤长沙航空职业技术学院，湖南长沙 410124)

0 引言

卫星、航空、地面、井中等多维空间的磁场数据测量技术迅速发展并逐渐成熟，配套的数据处理和解释技术也在不断更新和完善，这势必拓宽磁法勘探的应用领域，提高磁法的探测能力[1-4]。磁场数据三维反演不仅能够提供磁异常体的空间位置，还能获得磁性异常体的形状、物性分布等定量信息，在磁场数据处理和解释中具有重要意义，一直是该领域的研究热点问题。

目前的磁异常反演包括磁化率和磁性界面两方面的内容。磁化率反演首先固定反演网格，假设反演单元内磁化率为常数或者线性分布，建立反演目标函数，然后最小化反演目标函数，从而获得目标的磁化率分布[5-6]。Li等[6]利用一个或多个适当的加权函数将先验信息引入目标函数；Pilkington[7]对目标函数进行平滑及深度加权;Miguel[8]通过地质统计学模型实现了岩性约束下的重磁数据联合反演；Fregoso等[9]实现了重磁数据的交叉梯度联合反演。上述方法可定量描述地质构造的磁化率分布。磁性界面反演方法一般假定已知地下异常体的磁化率，通过反演得到异常体的边界和位置。管志宁等[10]提出了频率域磁性单界面反演方法，推导了常磁性与变磁性单界面反问题解的近似表达式；Pilkington等[11]提出了利用重磁资料确定地壳界面起伏形态；Wang等[12]对多面体的顶点位置进行反演；刘沈衡等[13]利用欧拉算法反演磁性界面；Fullagar等[14]通过混合参数化反演得到目标体磁化率的分布及垂直界面的位置。赵文举等[15]通过BP神经网络预测磁性体顶面埋深；郑强等[16]基于磁梯度张量特征值成功地识别了磁性体边界。然而，由于缺乏地下异常体拓扑结构和数量信息，现有的磁边界反演算法无法灵活描述复杂地下磁性目标体的几何形状。

因此，需要一种无需过多关于异常体拓扑结构和数量的先验信息就可以直接反演地下磁异常体空间位置和边界的算法。基于水平集方法[17-21]的磁边界反演算法正好符合这一要求，该方法仅需已知目标体的平均磁化率就可以反演得到地下异常体的边界及几何形状。Isakov等[22]将水平集方法用于重力数据反演；Zheglova等[23]将水平集用于地震旅行时层析成像，实现了边界的二维重建；Li等[24]在磁性数据的反演中引入水平集，用两个水平集对地下磁性体的边界进行了三维反演；在矿产勘查中，Li等[25]在地表数据与钻孔磁资料的联合反演中引入水平集，实现了目标圈定。基于水平集的反演可以得到异常体明确的边界，将水平集引入磁边界反演，是改善磁边界反演效果的重要途径。然而，目前已发表的基于水平集的磁性目标体边界反演算法仅使用了两个水平集，而实际地球物理勘探的目标往往是多个具有不同磁化率的磁性体。

针对上述问题，本文在Li等[24]的基础上提出了一种基于多重水平集的磁边界反演算法，该算法具有抗噪性强、精度高、多目标识别的优点。首先基于多个水平集函数建立磁界面反演目标函数，并详细推导了在多重水平集反演构架下磁边界反演的实现过程；然后采用任意四面体单元磁法解析解算法[24]进行高精度正演计算，构建正演网格与反演网格的映射关系，并基于加权基本无振荡(WENO)格式对水平集函数进行更新及重新初始化；最后，通过多个理论模型算例验证本文反演算法的有效性，并分析初始模型对反演结果的影响。

1 方法与原理

1.1 基于磁化率的多重水平集表示

假设非磁性背景下，区域Ω中包含N个磁性体，且每个磁性体具有恒定的磁化率，则磁化率分布可表示为

(1)

上述模型的磁化率可用水平集公式表示为

(2)

式中：φ(·)表示水平集函数；H(·)表示Heaviside函数[25]。其表达式分别为

(3)

(4)

(5)

(6)

式中δ(·)是Dirac delta函数。

1.2 基于多重水平集的磁边界反演算法

反演的总目标函数Ei可表示为数据拟合差泛函Ed和正则化项Er的线性组合

Et=Ed+αEr

(7)

1.2.1 数据失配项及其导数

(8)

式中：M是观测点数；σk是与第k个数据相关的误差标准偏差，数据拟合差的最优值一般为0.5。

(9)

(10)

(11)

(12)

(13)

令

(14)

则

(15)

1.2.2 正则化项及其导数

为了进一步降低反演的多解性，本文在水平集反演目标函数中引入Tikhonov正则化项

(16)

式中P代表水平集个数。Er的Fréchet导数为

(17)

(18)

1.3 水平集函数的更新及重新初始化

对均匀网格，利用WENO格式对水平集方程和重新初始化方程进行空间离散[29-34]，并利用三阶精确总变差递减龙格—库塔格式(TVD-RK3)[32]进行时间离散化。

1.3.1 水平集方程的近似

对于式(18)，使用TVD-RK3格式及Godunov法对H(φi)=Fi|φi|进行空间离散化，得到

(19)

式中：

使用TVD-RK3格式进行时间离散化。例如,对水平集方程φt+Fi|φ|=0，首先对φ(l)执行一次欧拉计算，得到t(l+1)时刻的解

(20)

(21)

(22)

通过线性平均更新φ(l+1)

(23)

并建立一个Courant-Friedrichs-Lewy(CFL)条件[34]

(24)

保持更新的稳定性。在本文反演中，正则化参数α远小于max|Fi|，所以设

(25)

1.3.2 重新初始化方程的近似

一般情况下，水平集函数在通过式(18)进行迭代时，并不保留其符号距离属性，所以Sussman等[19]引入了重新初始化方程

(26)

(27)

2 反演算例

2.1 反演流程

本文采用任意四面体单元磁法解析解算法[24]进行高精度正演计算，利用正、反演网格之间的物性映射实现正演网格与反演网格相互独立运行，并基于WENO格式对水平集函数进行更新及重新初始化，具体反演流程如下。

(1)初始化水平集函数φi；

(4)基于WENO格式计算总目标函数Et的负梯度方向-∂Et/∂φi，并更新水平集函数φi；

(5)重新初始化水平集函数φi，使其保持带符号距离属性；

(6)返回步骤(2)，直至迭代收敛或达到最大迭代次数。

2.2 两个水平集模型

设计一个由具有不同磁化率的两个倾斜棱镜体组成的模型，分析基于两个水平集公式的磁边界反演效果。首先根据先验信息确定两个不同的磁化率，并使用基于两个水平集公式的磁边界反演算法确定异常体数目并恢复磁性体边界。

设计模型如图1a所示，由κ1=0.04SI和κ2=0.08SI的两个倾斜棱镜体组成。图1b为加噪总磁异常图。初始猜测如图1c所示。反演中引入两个水平集函数φ1和φ2，初始猜测分别为

(28)

φ2=1-

(29)

图1d是使用多水平集反演得到的模型，其数据拟合差Ed=0.5022，可见反演结果与原模型的形状和位置吻合较好。图1e是磁场强度B0=50000nT、磁倾角为75°、磁偏角为25°的背景感应场下，反演模型(图1d)在观测面正演的磁异常图。图1f为反演模型的正演磁异常(图1e)与原模型加噪磁异常(图1b)之差。反演结果的更多细节见图2。显然，使用两个水平集的反演结果很好地恢复了每个棱镜体的倾角、空间位置及形态，且反演误差与预期的一样，呈随机分布特征。

图1 两个不同磁化率的倾斜棱镜模型及两个水平集计算结果(a)模型示意图； (b)加噪总磁异常图； (c)初始猜测模型； (d)两个水平集的反演模型； (e)图d正演磁异常图； (f)图e与图b的差值

图2 两个不同磁化率的倾斜棱镜体模型使用两个水平集的反演结果细节展示(a)y=300m截面； (b)y=700m截面； (c)z=150m平面图； (d)z=250m平面图虚线表示模型的真实边界，实线表示模型的恢复边界

将本文反演结果与Li等[24]的反演结果做对比，见图3。可以看出，本文采用基于任意四面体单元的磁法解析解算法[26]进行高精度正演，并引入HJ-WENO格式对水平集方程进行更新及重新初始化，反演得到的磁异常体边界与模型吻合更好。

图3 本文反演结果与Li等[24]的反演结果对比(a)y=300m截面(过磁性体①)； (b)y=700m截面(过磁性体②)虚线为模型的轮廓线，实线为本文反演磁性体边界，*线表示Li等反演的磁性体边界

2.3 三个水平集模型

实际地下地质情况是非常复杂的，需分析多个磁性体存在的情况下，本文方法是否有效。对于具有不同磁化率值的磁性体，需要引入更多的水平集函数，对其进行磁边界反演。

假设模型由三个具有不同磁化率的直立长方体组成(图4a)，且磁化率值已知。图4b为加噪总磁异常图。根据图4b，在反演算法中引入三个水平集函数φ1、φ2和φ3，初始猜测(图4c)分别为

φ1=1-

(30)

φ2=1-

(31)

φ3=1-

(32)

图4 三个不同磁化率模型及三个水平集计算结果(a)模型示意图； (b)加噪总磁异常图； (c)初始猜测模型； (d)三个水平集的反演模型；(e)图d正演磁异常图； (f)图e与图b的差值图a中磁性体①、②、③的长方体中心坐标分别为(500m,150m,175m)、(500m,475m,350m)和(500m,825m,300m)，大小分别为400m×200m×200m、200m×150m×200m和400m×150m×200m，磁化率分别为0.04、0.08、0.12SI

图5 图4d细节展示(a)y=200m截面； (b)y=500m截面； (c)y=800m截面； (d)z=100m平面图虚线为真实模型的轮廓线，实线为本文反演结果

2.4 深部异常体模型

前面两个模型中的磁异常体深度较小，说明了本反演算法对浅层磁性矿藏具有很好的圈定作用。为了进一步验证该方法对深部目标体边界的圈定效果，设计了图6a所示的模型。

模型包括两个磁性长方体，其中心坐标分别为 (2000m,1500m,900m)、(2950m,2550m,825m)，大小分别为1000m×500m×400m、700m×300m×350m,极化率已知，分别为0.16、0.30SI。基于加噪总磁异常数据(图6b)，反演中引入两个水平集函数φ1和φ2，其初始猜测(图6c)分别为

φ1=1-

(33)

φ2=1-

(34)

图6d为反演结果，其数据拟合差Ed=0.5916。图6e和图6f分别为反演模型的正演磁异常图及误差。图7是图6d的细节展示。由图可见，反演结果在各个方向上对目标体的边界都吻合得很好，展现了真实磁性体的形状和空间位置。

图6 深部磁异常体模型及反演结果(a)模型示意图； (b)加噪总磁异常图； (c)初始猜测模型； (d)基于两个水平集的反演结果； (e) 图d正演磁异常图； (f)图e与图b的差值

图7 图6d细节展示(a)z=800m平面图； (b)x=2000m截面；(c)x=3000m截面图。虚线表示原始模型轮廓线，实线表示反演模型边界

本算例验证了本文反演算法对深部矿产的圈定效果很好。

2.5 水平集数目对反演结果的影响

基于多重水平集方法的磁界面反演算法的关键假设是地下异常体的磁化率已知。实际反演中，一般可以基于已有的地质知识或其他先验信息得到研究区磁性体确切的磁化率值，但无法确切知道某一区域具体对应哪个磁化率值。复杂的地质环境中，这个不可预知的因素可能会极大地限制多重水平集方法的应用。本节针对图4a所示模型，假设仅知道一个磁化率值，分析水平集数目的选择对反演结果的影响。

假设图4a中仅知道②号异常体的磁化率(κ=0.08SI)，并据此形成初始猜测(图8a)，并在反演中引入一个水平集函数

φ=1-

(35)

图8b为引入单水平集的反演结果。可见即使初始模型只假设了一个磁源体，反演结果也能准确地反映有三个独立的磁源体，说明水平集的个数对反演结果影响不大。

图9是图8b的细节展示。对比图9a与图6a可以清晰地看出，当κ=0.04SI的①号长方体被赋予大于其实际磁化率值(κ=0.08SI)时，反演得到的磁性边界范围更小、深度更大；对比图9b与图6a可以看出，当磁化率为κ=0.12SI的②号长方体被赋予低于其实际磁化率值(κ=0.08SI)时，反演得到的边界范围更大，且深度小于实际深度；图9c中，反演的③号磁性体边界与模型(图6a)几乎完全一致，这是因为反演时设定的磁化率值与实际情况完全一致。

图8 图6a模型单水平集反演结果(a)初始猜测； (b)反演模型

图9 图8b细节展示(a)y=200m截面(过①号磁性体)； (b)y=500m截面(过②号磁性体)； (c)y=800m截面(过③号磁性体)虚线表示原模型磁性体轮廓线，实线表示反演模型边界

该模型反演结果说明，基于水平集方法反演时模型的磁化率与真实值差别太大时，恢复的磁性体边界会出现较大偏差：如果反演模型的磁化率赋值低于实际值，反演结果中磁性源的体积会大于实际体积，且位置偏深；反之则会导致反演的磁源体积小于实际体积，且位置偏浅。

由本节反演结果可知，反演初始模型的各个子域磁化率值是否准确，对反演结果影响较大，而反演初始模型水平集的数目对反演结果的影响较小。因此，在实际应用中，需要利用先验信息尽量准确地估计各个区域的磁化率值，以各区域的平均磁化率作为各子域的磁化率进行多水平集反演，才可得到比较可靠的磁边界位置。

3 结论

本文基于多重水平集的磁界面反演算法，实现了对磁化率已知的磁性目标体空间位置和几何形状的反演。

(1)本文算法中，磁化率值由先验信息确定，采用任意四面体单元磁法解析解算法[24]进行高精度正演计算，并基于WENO格式对水平集函数进行更新及重新初始化。

(2)具有两个水平集和三个水平集的反演算例验证了本算法的有效性与可靠性；深部异常体模型反演算例验证了本算法对深部矿产资源的圈定具有良好效果。

(3)针对实际生产中，难以准确获取地下目标磁化率的情况，讨论了水平集数目的确定对反演结果的影响。即便水平集数目与实际不符，反演结果仍能准确划分磁性体区域，但各反演磁性体的深度和边界存在一定误差。若子域被赋予偏大的磁化率值时，反演得到的磁性体范围小于实际情况，且深度偏小；反之，若子域被赋予较小磁化率值时，则反演磁性体的范围大于实际情况，且深度偏大。