空间—波数域二维强磁场及其梯度张量数值模拟研究

戴世坤，冉应强*，张莹，陈轻蕊，凌嘉宣，贾金荣

(1.中南大学有色金属成矿预测与地质环境监测教育部重点实验室，湖南长沙 410083； 2.中南大学地球科学与信息物理学院，湖南长沙 410083)

0 引言

近年来，强磁性体的退磁效应逐渐成为磁场正、反演领域的研究热点。当磁性体磁化率低于0.1 SI时，可忽略退磁效应[1]；当磁性体磁化率高于0.1 SI时，退磁效应会显著影响磁异常的幅值、形态等，进而影响磁测数据的处理和解释结果，所以在强磁性体磁场正、反演过程中必须考虑退磁效应的影响。考虑到三维模型的计算量太大，且实际勘探过程中，当地质体沿走向的尺度远大于垂直走向的尺度时(如断层、接触带等)，这类地质体可用走向方向无限延伸的二度体近似，既能满足勘探要求，也能大大提高计算效率[2]。

目前，对于二度体磁异常的研究多是基于弱磁情况，可忽略退磁效应。Bhattacharyya等[3-4]、Peder-sen[5]、吴宣志[6]给出了频率域不同二度体模型的弱磁异常场表达式；Talwani等[7]、Shuey等[8]、Plouff[9]、Singh等[10]、Won等[11]、贾真[12]提出了空间域二度体数值模拟算法；对于频率域弱磁正演计算方法，柴玉璞[13]基于偏移抽样提高了反傅里叶变换数值精度；Tontinif等[14]研究了快速傅里叶变换扩边法与误差的关系；Wu等[15]在偏移采样的基础上提出了高斯傅里叶变换方法，有效提高了弱磁场数值模拟的精度和效率。

考虑退磁效应的情况下，强磁场数值模拟方法主要有以下三类：积分方程法、有限体积法和有限单元法。Eskola等[16]基于积分方程提出了考虑退磁效应的静磁方程，其解的准确性与代数方程组的规模有关；Ouyang等[17]基于强磁场积分方程，采用棱柱体剖分方法实现了强磁性体的正演计算，采用Gauss-FFT方法减小了截断效应；Kostrov[18]提出了一种采用三角单元的体积分方法求解考虑退磁效应的二维磁场正演问题；付文祥[19]基于体积分方法提出一种采用三角棱柱单元拟合二度体的方法，并给出了考虑退磁影响的无限长水平圆柱体的总磁场强度异常的理论解析表达式；欧洋等[20]利用有限体积法实现了同时考虑退磁和剩磁的正演模拟；王书惠[21]在不引入均匀磁化设定、考虑退磁效应影响的情况下，利用有限元计算了强磁性复杂形体的磁异常；刘双等[22]利用FlexPDE软件实现了二度圆柱体、板状体的有限元强磁场正演，总结了退磁作用对磁异常幅值影响的规律；另外，刘鑫磊[23]计算了具有各向异性、非均匀磁化率分布的任意形态地质体的磁场，计算过程考虑了退磁效应，但计算结果可能会不连续；Purssm等[24]采用迭代方法推导了复杂地形下高磁化率地质体的磁场表达式，但计算结果误差较大。

现有的强磁场数值模拟方法大多基于有限单元法、有限体积法等，这些方法最终都归结为大型线性方程组的求解，对于复杂条件下的大规模二维强磁场数值的模拟存在计算效率低的问题。此外，由重磁位场的研究结论可知，与总磁场相比，梯度张量对异常体边界的识别能力更强、分辨率更高，而现有关于强磁场梯度张量的数值模拟研究较少。因此，本文提出一种适用于大规模复杂条件(复杂磁化率分布和起伏地形)的二度体强磁场及其张量梯度的正演方法——空间—波数混合域二维强磁场及梯度张量数值模拟方法。该方法利用一维傅里叶变换将空间域二维强磁性体磁位满足的二维偏微分方程转换为空间—波数混合域磁位满足的一维常微分方程，利用有限单元法对每个波数下的空间—波数混合域磁位进行计算，并利用追赶法求解定带宽线性方程组，得到空间—波数混合域磁位，利用位场与张量梯度之间的关系得到相应的空间—波数域磁场分量，再经过反傅里叶变换求得空间域场。该算法通过迭代对强磁场进行数值逼近，确保了算法稳定和收敛。

1 方法原理

1.1 控制方程

强磁性体二维磁位方程[25]为

2U(x,z)=·M(x,z)

(1)

关于磁介质受到外部磁场的作用，可用磁化强度M衡量磁化程度，它与磁场强度H之间的关系为

(2)

对式(1)沿x方向进行一维傅里叶变换，得到空间—波数混合域磁位方程

(3)

1.2 边界条件

式(3)的通解为

(4)

式中a和b为常数。在笛卡尔坐标系中，令z轴向下为正，模型的上边界z=zmin，下边界z=zmax，根据上行波和下行波相关理论，上、下边界条件分别为

(5)

(6)

1.3 磁场和磁场梯度张量

对磁位求导可得到相应的磁异常场及其梯度张量，空间域磁异常场Ba、磁异常场梯度张量T和总强度磁异常ΔT[26]的计算公式分别为

(7)

(8)

ΔT=‖B0+Ba‖-‖B0‖

(9)

(10)

(11)

1.4 空间—波数域常微分方程的求解

联立式(3)、式(5)、式(6)，可求解空间—波数混合域磁位的边值问题。

对于二维强磁性体，磁位在空间—波数混合域满足边值问题，可采用基于二次插值的一维有限单元法对其进行求解。该方法保留垂向方向为空间域，可根据实际需求灵活改变垂向网格密度，兼顾了计算精度和计算效率；利用有限元法得到的五对角线性方程组，可利用追赶法[27]实现方程组的高效、高精度求解。

基于变分原理[25]，得到与边值问题等价的变分问题

(12)

(13)

(14)

由于δu≠0，故

Ku=P

(15)

1.5 紧算子

求解空间—波数混合域磁位的边值问题时，磁化强度M的表达式(式2)中存在未知项Ha，若直接求解，不能得到五对角线性方程组。本文采用迭代法对真解进行逐次逼近。在电磁法数值模拟中，Torres-Verdín等[28]提出了一种电磁场扩展的Bron近似法，这是一种基于内部场的新型近似法，通过将背景电场投影到散射张量上对电磁场进行近似表示，该方法适用于异常体体积较大、电阻率差异明显及宽频的情形。Zhdanov等[29]构建了一个收敛的Bron级数，Gao[30]进一步修改得到新的收敛Bron近似。Ouyang等[17]基于前人研究，根据电磁场紧算子推导过程构造了适用于强磁场稳定收敛的迭代格式，其迭代公式为

(16)

式中：i表示迭代次数;H(i)和H(i+1)分别表示第i次和第i+1次迭代得到的总场。

1.6 算法流程

本文提出的适用于大规模复杂条件(复杂磁化率分布和起伏地形)的二度体强磁场及其张量梯度的正演方法，即空间—波数混合域二维强磁场及梯度张量数值模拟方法流程见图1，具体步骤如下：

图1 空间—波数混合域二维强磁场及梯度张量数值模拟迭代算法流程

(1)输入背景场磁场强度H0，即第一次迭代时总磁场强度的初始值H(0)；

2 算例分析

首先设计一个二维圆柱体模型和一个棱柱体模型分别验证本文算法的正确性和高效性，然后设计一个带起伏地形的、包含一个顺磁性异常体和一个强磁性异常体的模型，验证正演算法对复杂条件的适应性。算例采用串行计算，算法中的傅里叶变换通过Gauss-FFT实现。算法基于Fortran95语言编程，电脑配置为：Intel Core i3-4150 CPU，主频为3.50 GHz，内存为12 GB。

2.1 正确性验证

设计图2所示二维圆柱体模型，圆柱体沿y方向无限延伸，顶部埋深为300 m。模型的研究区域为：x=[-2000 m，2000 m]，z=[0，1000 m]；网格数为800(x)×400(y)，水平采样间隔为5 m，垂向采样间隔为2.5 m；观测面为z=0。圆柱体磁异常场的x分量Bax、z分量Baz、梯度张量水平分量∂Bax/∂x和垂向分量∂Bax/∂z的解析表达式详见附录B。该模型的总强度磁异常ΔT、磁异常场Ba及其梯度张量T的解析解和数值解及观测面上各点的相对误差分别见图3、图4和图5。

图2 圆柱体模型示意图

图3 圆柱体模型总强度磁异常ΔT解析解和数值解(a)及其相对误差(b)

图4 圆柱体模型磁异常场Bax(上)和Baz(下)的解析解和数值解(a)及其相对误差曲线(b)

图5 圆柱体模型磁异常场梯度张量∂Bax/∂x(上)和∂Bax/∂z(下)的解析解和数值解(a)及其相对误差(b)

由图3～图5可以看出，ΔT的相对误差小于0.5%，Bax和Baz的相对误差小于0.5%，∂Bax/∂x和∂Bax/∂z的相对误差小于0.1%，即ΔT、Bax、Baz、∂Bax/∂x和∂Bax/∂z的解析解和数值解基本吻合，仅在极少数零值点附近相对误差较大。该圆柱体模型的正演结果验证了本算法的正确性。

图6为该模型的磁异常水平分量Bax和垂直分量Baz的相对误差随迭代次数的变化曲线。可以看出，经9次迭代后，计算结果即可满足迭代终止条件。

图6和图7表明本文迭代算法适应于不同磁化率的模型，收敛稳定。

图6 磁异常分量相对误差随迭代次数变化曲线

图7 迭代次数随圆柱体磁化率的变化曲线

2.2 算法效率分析

图8 棱柱体模型示意图

参照文献[31]中利用有限体积法计算得到的数值解，衡量COMSOL软件和本文算法的精度。引入相对均方根误差rrms[33]统计整条观测线上总强度磁异常ΔT的相对误差

(17)

式中：N表示观测总点数；Xi表示COMSOL或者本文算法计算得到的第i个观测点的总强度磁异常；Yi表示有限体积法计算得到的第i个观测点的总强度磁异常。

当网格数为200×100时，采用本文算法的rrms为2.84%，计算时间0.75 s，占用内存12.0 MB。表1为不同网格剖分方案下基于COMSOL软件的计算精度、耗时和内存占用。基于COMSOL计算时需对研究范围进行扩边处理，当网格扩边4倍时，rrms为2.92%，计算时间618 s，占用内存2387.1 MB。因此，与COMSOL软件相比，在获得相似精度的情况下，本文算法耗时更短，占用内存更少，计算效率更高。

表1 COMSOL计算时间与占用内存统计

网格数rrms%计算时间s占用内存MB200×100(不扩边)110.901480.5600×300(扩边2倍)9.57130887.41000×500(扩边4倍)2.926182387.1

有限体积法、COMSOL软件(扩边4倍)和本文算法的总强度磁异常见图9a，以有限体积法的数值解为参照，采用COMSOL软件和本文方法计算结果的相对误差见图9b。可以看出，采用COMSOL软件和空间—波数域方法计算得到的总强度磁异常与有限体积法的计算结果图形均大致吻合。

图9 棱柱体模型不同方法得到的总强度磁异常ΔT(a)及相对误差(b)

对于图8模型，改变模型网格剖分个数，分析本文算法计算时间和占用内存随网格节点数的变化，结果见图10和表2。对比表1与表2，可见本文算法计算时间更短、占用内存更少，占用内存和计算时间与网格节点数均呈现近似线性正相关变化趋势。

图10 棱柱体模型本文算法计算时间和占用内存随网格剖分节点数的变化

表2 棱柱体模型本文算法计算时间与占用内存统计

网格数计算时间/s占用内存/MB200×1000.7512.0400×2002.9443.1600×3005.96108.8800×4009.50183.81000×50014.39280.2

2.3 适应性验证

为了验证本文正演算法对复杂模型的适应性，设计一个起伏地形模型(图11)，模型同时包含顺磁性异常体和强磁性异常体。模型的背景磁场T0=50000 nT，磁倾角α=45°，磁偏角β=5°。研究区域：x=[-2000 m，2000 m]，网格数为400；z=[-320 m，820 m]，其中z方向区间[320 m，820 m]为起伏地形最低点以下的计算区域，网格数为100。起伏地形高程满足函数

图11 起伏地形模型示意图

x∈[-2000,2000]

(18)

基于该模型研究起伏地形条件下插值平面个数对数值解计算精度的影响。起伏地形在z方向的区域为[-320 m，320 m]，该方向采用均匀网格剖分，将插值平面个数分别设为5、9、13、17、21、25、29和33。采用式(17)计算整条起伏观测线上的总强度磁异常ΔT、磁异常场Ba及梯度张量T的解析解。

图12为该起伏地形模型的总强度磁异常、磁异常场及梯度张量的解析解和数值解的相对均方根误差随插值平面个数的变化曲线。可以看出，总强度磁异常、磁异常场与其梯度张量的rrms均随插值平面个数的增加而降低。采用5个平面进行插值时，起伏地形上总强度磁异常、磁异常场与其梯度张量的rrms均大于1%；当采用25个及以上不同高度平面进行插值时，起伏地形上总强度磁异常、磁异常场与其梯度张量的rrms均低于1%，约为0.9%，说明采用25个插值平面即可达到满意的计算精度。

图12 rrms随插值平面数的变化曲线

另外，采用25个插值平面时，耗时约2.330 s，占用内存52.8 MB，表明本文算法计算速度快、占用内存小，适用于复杂地形下二度体磁异常场及其梯度张量的计算，计算效率较高。

采用25个高度平面进行插值，基于本文方法得到的总强度磁异常、磁异常场及其梯度张量分别见图13、图14和图15。可以看出，解析解和数值解的吻合度很高，表明了本文算法对复杂地形模型的计算精度较高。

图13 起伏地形模型总强度磁异常解析解与数值解对比

图14 起伏地形模型磁异常场解析解和数值解对比

(a)Bax； (b)Baz

图15 起伏地形模型磁异常场梯度张量解析解和数值解对比

(a)∂Bax/∂x； (b)∂Bax/∂z

3 结论

针对强磁性体的退磁效应不能忽略的问题，本文提出了一种空间—波数混合域二度体磁异常场及其梯度张量正演方法。该算法充分利用了傅里叶变换的快速性、追赶法求解的高效性和迭代算法的稳定性的优势，实现了二维强磁性体磁异常场及其梯度张量的高效、高精度数值模拟。设计了水平地形圆柱体模型进行数值试验，通过对比解析解与数值解验证了本算法的正确性和较高的精度；以有限体积法数值解为参照，对比了本文算法和COMSOL软件的计算效率，证明了本算法的高效性；利用起伏地形下同时包含顺磁性异常体和强磁性异常体的复杂模型验证了本算法对复杂地质条件的适应性。另外，在本文算法的计算环节中，正、反傅里叶变换和求解常微分方程这两个环节均有很好的并行性，采用并行方式可进一步提高计算效率。本文为地下二度体磁异常场及其张量梯度的正演提供了一种高效、高精度算法，对于二度体磁异常的反演和人机交互正、反演解释具有重要意义。

附录A 三维强磁场变分问题求解

求解变分问题(式(12))时，可将其写为

(A-1)

式中

(A-2)

令

(A-3)

式中：j、m分别为单元两端节点坐标；p为单元中点坐标；积分区域q为单元范围。可得

(A-4)

式(A-1)中第一项单元积分为

(A-5)

其中

(A-6)

式中l为单元长度。

式(A-1)中第二项单元积分为

(A-7)

其中

(A-8)

式(A-1)中第三项单元积分为

(A-9)

利用形函数可将上式中的Jaz表示为

Jaz=JazjNj+JazpNp+JazmNm

(A-10)

其中

szq=(Jazj,Jazp,Jazm)T

(A-11)

则式(A-9)中右端项可表示为

(A-12)

其中

(A-13)

式(A-1)中，z=zmin、z=zmax分别为第一个和最后一个节点，可将其分别扩展为

(A-14)

其中

(A-15)

(A-16)

综上所述，式(A-1)可表示为

F(u)=uTKu-2uTP

(A-17)

附录B二维圆柱体强磁场及其梯度张量解析解

Ha=

(B-1)

(B-2)

沿走向无限长圆柱体的磁场梯度张量

(B-3)

其中

H0z(z-z0)[(z-z0)2-3(x-x0)2]]

(B-4)

H0z(x-x0)[(x-x0)2-3(z-z0)2]]

(B-5)

式中H0x、H0z分别为背景磁场H0的x分量和z分量。