地球物理学报  2019, Vol. 62 Issue (11): 4437-4450   PDF    
引用本文
李昆, 戴世坤, 陈轻蕊, 等. 2019. 空间波数混合域磁异常场积分解三维数值模拟. 地球物理学报, 62(11): 4437-4450, doi: 10.6038/cjg2019M0362.  
Li K, Dai S K, Chen Q R, et al. 2019. Three-dimensional modeling of magnetic anomaly integral solution in a mixed space-wavenumber domain. Chinese J. Geophys. (in Chinese), 62(11): 4437-4450, doi: 10.6038/cjg2019M0362.  
空间波数混合域磁异常场积分解三维数值模拟
李昆1,2,3, 戴世坤1,2,3, 陈轻蕊1,2,3, 张钱江5, 赵东东1,2,3, 王顺国4, 凌嘉宣1,2,3     
1. 有色资源与地质灾害探查湖南省重点实验室, 长沙 410083;
2. 中南大学有色金属成矿预测与地质环境监测教育部重点实验室, 长沙 410083;
3. 中南大学地球科学与信息物理学院, 长沙 410083;
4. 加州大学圣地亚哥分校Scripps海洋研究所, 美国 92092;
5. 桂林理工大学地球科学学院, 广西桂林 541004
收稿日期 2018-06-13, 2018-12-02 收修定稿
基金项目: 国家重点研发计划（2017YFC0601305，2018YFC0603602），国家自然科学基金项目（41574127），中南大学研究生自主探索创新项目（2017zzts176），湖南省自然科学基金（2018JJ3642），广西自然科学基金（2018GXNSFAA050070），桂林理工大学科研启动基金资助项目（GUTQDJJ2018034）联合资助
第一作者简介: 李昆, 男, 1989年生, 博士, 研究方向为重、磁、电、震三维正反演理论方法研究.E-mail:155001015@csu.edu.cn
通讯作者: 戴世坤, 男, 1963年生, 教授, 博士生导师, 研究方向为重、磁、电、震三维正反演理论方法研究.E-mail:dskgmes@csu.edu.cn
摘要：本文提出一种空间波数混合域磁异常场三维数值模拟方法.该方法利用磁位三维空间域积分为卷积的特点，沿水平方向进行二维傅里叶变换，把空间域磁位满足的三维积分问题转化为不同波数之间相互独立的垂向一维积分问题.保留垂向为空间域，优势之一在于便于浅层单元剖分可适当加密，随着深度增加，单元剖分适当稀疏，可以准确模拟任意复杂地形和磁性体的磁异常，兼顾了计算精度与计算效率；优势之二在于一维积分垂向可离散为多个单元积分之和，每个单元采用二次形函数表征磁化强度，可得出单元积分的解析表达式，计算精度高、效率高.该方法充分利用一维形函数积分的高效和高精度、快速傅里叶变换的高效性及算法高度并行性，实现了磁异常场高效、高精度的数值模拟.设计棱柱体模型，将模型解析解与空间波数混合域法的数值解对比，结果表明该方法计算精度高、效率高.设计了组合棱柱体复杂模型，对比分析了标准FFT扩边法与Gauss-FFT法的计算精度与计算效率，总结了标准FFT的扩边系数选取策略.针对任意复杂地形条件下的磁异常模拟问题，本文提出一种适用于起伏地形条件下的磁异常场快速计算方法，并对其有效性进行了验证.
关键词: 磁异常场三维数值模拟      空间波数混合域      傅里叶变换      形函数法     
Three-dimensional modeling of magnetic anomaly integral solution in a mixed space-wavenumber domain
LI Kun1,2,3, DAI ShiKun1,2,3, CHEN QingRui1,2,3, ZHANG QianJiang5, ZHAO DongDong1,2,3, WANG ShunGuo4, LING JiaXuan1,2,3     
1. Hunan Key Laboratory of Nonferrous Resources and Geological Hazards Exploration, Changsha 410083, China;
2. Key Laboratory of Metallogenic Prediction of Nonferrous Metals and Geological Environment Monitoring of Ministry of Education, Central South University, Changsha 410083, China;
3. School of Geosciences and Info-Physics of Central South University, Changsha 410083, China;
4. Scripps Institution of Oceanography, UC San Diego, San Diego 92092, USA;
5. School of College of Earth Sciences of Guilin University of Technology, Guilin Guangxi 541004, China
Abstract: The paper presents a three-dimensional (3D) modeling method, engaged with a mixed space-wavenumber domain, for the magnetic anomaly field. Based on the fact that the 3D space domain integral for magnetic potential is convolution, the 3D integral problem of magnetic potential in space domain can be transformed into one dimensional (1D) integral problems in vertical direction independently to each other using the two-dimensional (2D) Fourier transform along the horizontal directions. By the proposed method, the depth direction is kept in space domain, which can provide two advantages in the modeling. Firstly, the shallow grid can be as fine as necessary for topography, however, with the increase of depth, the grid becomes coarse in order to reduce calculation cost. Thus, the calculation accuracy and efficiency in modeling can be guaranteed simultaneously. Secondly, the one-dimensional integral can be discretized into the sum of multiple element integrals along the depth direction. Each unit uses the second-order shape function to calculate the magnetization with the analytical expression of element integral. This method makes full use of 1D shape function for high accuracy, high parallelization among different wavenumbers, and fast Fourier transform to achieve high efficiency and high accuracy modeling of magnetic anomaly field. A prismatic model is designed to verify the correctness and high accuracy of the method by the comparison between the analytical solution and the numerical solution obtained by the proposed method. A complex model is designed to compare the accuracy and efficiency of the standard FFT expansion method and of the Gauss-FFT method. The strategy of edge selection for standard FFT method was also summarized based on the complex model. The paper presents a fast magnetic anomaly field calculation method and verify its validity in order to solve the problem of magnetic anomaly simulation under any complex terrain conditions.
Keywords: Magnetic anomaly 3D modeling    Mixed space-wavenumber domain    Fourier transform    Shape function method    
0 引言

任意复杂条件(规模面积大，地形复杂，异常体复杂)下的磁异常场高效、高精度数值模拟在磁异常反演成像及人机交互解释中有着重要作用，目前关于磁异常场正演的计算方法主要分为空间域法与频率域法.

空间域方法是通过规则形体的磁场解析式直接求得空间任意点的精确解.空间域方法的优点是计算精确，缺点是对于复杂形体，解析式复杂或者不存在.在空间域，磁场正演先期研究侧重于简单、规则形状的长方体、直立棱柱体、直立圆柱体等解析表达式的推导(Bhattacharyya，1964；Plouff，1976；Singh and Sabina, 1978), 之后学者们对均匀介质的多面体模型进行了详细研究(Furness，1994；Ivan，1996；Guptasarma and Singh, 1999；Singh and Guptasarma, 2001).同时，对于任意三度体简单物性变化模型的重磁场解析式推导(Pohánka，1998；Chakravarthi et al., 2002；García-Abdeslem，2005)、解析积分奇异点处理(Kwok，1991；Nagy et al., 2000；郭志宏等，2004；骆遥和姚长利，2007)、磁梯度张量目标定位(Hu et al., 2019)也取得了一定进展.近年来部分学者将压缩存储技术、GPU并行技术等应用于重磁位场正演计算中来提高计算效率(姚长利，2003；陈召曦，2012).

磁场正演在频率域中的计算称为谱正演法.谱正演法的优点是频谱表达式比较简洁，有些在空间域不能得到解析式的模型，在频率域容易得到(例如物性参数随深度呈指数衰减模型)，借助快速傅里叶变换，频率域的计算效率较高；缺点是快速傅里叶变换数值方法存在边界截断效应问题，限制了频率域正演方法的推广.在频率域，Pedersen(1978)首先给出了任意三度体(以三角形组合为其表面的多面体模型)的重磁异常频率域表达式，为三度体重磁异常的频率域正演提供了理论基础.之后人们陆续推导了物性连续变化及截面为简单几何形状的频率域表达式(吴宣志, 1981, 1983；Pedersen，1985).熊光楚(1984)推导了偶极子全空间磁位的三维傅里叶变换，并在此基础上推导了长方体模型的磁位、磁场.冯锐(1986)推导了Parker公式介质物性在水平方向任意变化与深度方向指数或线性变化的情况.柴玉璞(1988, 1996)采用偏移抽样方法，大大提高了离散傅里叶变换的精度，减小了基于傅里叶变换方法的位场数值模拟的误差.Tontini(2005)在前人基础上推导了物性为三维高斯分布源体的磁场频率域表达式，并将该方法应用于反演中.Tontini等(2009)对任意三度物性分布频率域正演方法进行了研究，并对减小频率域正演误差的扩边法与误差大小的关系作了研究.Wu和Tian(2014)和Wu(2016)将Gauss-FFT法应用于重磁位场数值模拟中，该方法对于压制边界效应问题提供了一种有效途径，数值模拟精度较高.Ren等(2017)提出一种基于非结构化网格的自适应快速多极子重磁场及其梯度张量正演方法，借助并行计算技术，实现了快速、高精度三维正演.

目前任意复杂条件下磁异常场数值模拟，为了得到较高的计算精度，需要对模型进行精细剖分，导致计算量和存储量剧增，使得大规模复杂条件下的磁异常正演需要高性能计算机，计算时间长.针对这一问题，本文提出一种空间波数混合域磁异常场三维数值模拟方法，该方法沿水平方向做二维傅里叶变换，把空间域磁位满足的三维积分问题转化为不同波数相互独立的垂向一维积分问题.保留垂向为空间域，便于浅层单元剖分适当加密，随着深度增加，单元剖分适当稀疏，可以准确模拟复杂地形和磁性体的磁异常，兼顾了计算精度与计算效率；垂直方向一维积分离散为多个单元积分之和，每个单元采用二次形函数描述磁化强度，可得出单元积分的解析表达式，计算精度高、效率高.该方法充分利用一维形函数积分的高效和高精度、快速傅里叶变换的高效性及算法高度并行性, 实现了磁异常场高效、高精度的数值模拟.模型算例中，设计了棱柱体模型，模型解析解与空间波数混合域法数值解对比表明，该方法正确、精度高、效率高.设计了组合棱柱体复杂模型，对比分析了标准FFT扩边法与Gauss-FFT法的计算精度与计算效率，总结了标准FFT的扩边系数选取策略.针对复杂地形条件磁异常模拟问题，本文提出一种适用于起伏地形条件下的磁异常场快速计算方法，并对其有效性进行了验证.

1 理论与方法

弱磁性体磁位空间域三维积分公式(Blakely，1996)：

(1)

其中，格林函数g(r, r′)=1/(4πR)，R= ，(x, y, z)为观测点坐标，(x′, y′, z′)为磁源分布坐标，梯度算子为真空磁导率，μ=4π×10－7H·m-1.磁化强度M(r′)=Ji+Jr, Ji为感应磁化强度，Jr为剩余磁化强度，单位均为A·m-1.感应磁化强度Ji=χ(H0+Ha)，H0为背景磁场，Ha为异常磁场，单位均为A·m-1，为磁化率，单位SI.对于弱磁性体Ha≪H0，可忽略异常磁场Ha，则感应磁化强度表达式简化为Ji=χH0.

对式(1)中x, y两个方向做傅里叶变换得空间波数混合域一维积分公式：

(2)

其中，为空间波数混合域磁化强度三个方向分量，kx, ky分别为x, y方向波数，为空间波数混合域格林函数，即g(r, r′)的二维傅里叶变换

(3)

其中，波数

将式(3)代入式(2)得空间波数混合域一维积分式：

(4)

其中，

空间域磁异常场、梯度张量与磁位的关系分别为(Blakely，1996)

(5)

(6)

其中，B为空间域磁异常场，单位nT，Bt为空间域磁异常梯度张量，单位nT·m－1.

根据傅里叶变换微分性质，由式(5)与式(6)得空间波数混合域磁异常场、磁异常梯度张量与磁位的关系分别为

(7)

(8)

其中，为空间波数混合域磁异常场，为空间波数混合域磁异常梯度张量，磁张量矩阵对称，独立分量为6个.在计算磁异常场和梯度时，若在异常内部，∂U/∂z、∂Bz/∂z利用形函数求偏导数法计算.

式(4)为空间波数混合域垂向(z方向)一维积分表达式，从表达式可以看出，不同波数一维空间域积分是相互独立的，可以并行计算.保留垂向为空间域，优势之一在于垂向单元可以进行灵活剖分，为了适应任意复杂地形和磁性体的磁异常场高精度计算，浅层单元剖分可适当加密，随着深度增加，单元剖分适当稀疏，同时兼顾了计算精度与计算效率；优势之二在于一维积分深度方向可离散为多个单元积分之和，每个积分单元采用有限单元法中的形函数表征磁化强度，可得出单元积分的解析表达式，计算精度高、效率高.

2 用形函数法计算一维积分

本文采用基于二次插值的形函数法计算空间波数混合域一维积分.将一维积分沿深度方向剖分为m个单元，每个单元内磁化强度要求满足二次变化，采用有限元单元法中的二次形函数描述磁化强度.

对空间波数混合域磁位积分公式(4)进行离散：

(9)

其中，Ix(kx, ky, z)、Iy(kx, ky, z)和Iz(kx, ky, z)分别为

(10)

(11)

(12)

其中，z′i和z′i+2分别为第i个单元积分的下限和上限.

将磁化强度表示为节点磁化强度的二次插值函数：

(13)

(14)

(15)

其中，N1, N2, N3为二次形函数.

将式(13)—(15)代入式(10)—(12)中，整理得

(16)

(17)

(18)

其中，Wi1、Wi2和Wi3为形函数积分系数，单元积分解析表达式详细推导过程见附录A.

由式(16)—(18)，离散后的一维单元积分可推导出解析表达式，计算精度近似于解析解，计算精度高.对于磁化强度满足二次变化的复杂形体，垂向单元采样间隔可以较大，可根据需要进行灵活选择.

3 起伏地形条件下磁异常场快速算法

水平地形条件下，测点在同一平面上，磁异常场计算只需进行一次积分，通过反傅里叶变换得水平面磁异常场.在起伏地形条件下，需要计算地形最高点与地形最低点之间n个深度节点所在平面磁异常场，然后通过插值得到起伏地形上磁异常场，计算量是水平地形条件下的n倍，计算效率低.针对这一问题，本文提出一种起伏地形磁异常场数值模拟快速算法.起伏地形模型示意图如下图 1所示.

图 1 起伏地形示意图 Fig. 1 A model with topography

以Ix一维积分为例，积分表达式为

(19)

其中，zp为起伏区域内不同深度节点，p=1, 2, …, n.

一维积分Ix(kx, ky, zp)积分区域由三部分组成(图 1所示中的(1)、(2)、(3))，式(19)可表示为

(20)

其中，

(21)

上式中，m1、m2和m3分别为三部分积分区域单元个数.

分析式(21)可知该积分具有以下特点：(1)计算不同深度节点zp所在面磁异常时，Ix3(kx, ky, zp)是相同的，即对于不同的深度节点zp，地形最低点以下的积分只需计算一次；(2)随着深度节点zp从地形最低点z1向上移动，深度节点zp的积分Ix1(kx, ky, zp)、Ix2(kx, ky, zp)与深度节点zp-1的积分Ix1(kx, ky, zp-1)、Ix2(kx, ky, zp-1)存在积分区域重叠，重叠区域可只计算一次；(3)当计算深度节点z1时，积分区域变为两部分Ix2(kx, ky, z1)、Ix3(kx, ky, z1)，当计算深度节点zn时，积分区域变为两部分Ix1(kx, ky, zn)、Ix3(kx, ky, zn).

根据式(21)积分特点，当计算深度节点zp所在面磁异常场时，可利用上一个深度节点zp-1的积分计算结果，与每个深度节点zp均对整个积分区域直接积分累加算法相比，避免了重复计算，提高了计算效率.快速算法具体计算过程为：

(1) 计算起伏面最低点z1节点所在平面磁异常场：由公式(21)，计算z1节点上下两部分积分Ix2(kx, ky, z1)与Ix3(kx, ky, z1)，并对Ix3(kx, ky, z1)与Ix2(kx, ky, z1)每个单元积分结果进行存储;

(2) 计算下一个深度节点z2：由积分公式特点可知，z2节点的Ix3(kx, ky, z2)与z1节点的Ix3(kx, ky, z1)相同，可直接调用已计算存储结果；z2节点的Ix2(kx, ky, z2)积分区域包含在z1节点Ix2(kx, ky, z1)积分区域以内，每个单元计算结果已经计算存储，可直接调用；计算Ix1(kx, ky, z2)并存储;

(3) 其他深度节点依次类推，直到计算出最高点zn节点所在平面磁异常场，然后通过插值得到起伏地形上磁异常场.

4 算例分析

设计棱柱体模型，采用基于四个高斯点的Gauss-FFT法(Wu and Tian, 2014)计算磁异常场及张量数值解，并与模型解析解(管志宁，2005)进行对比，验证该算法的正确性；研究了基于标准FFT扩边法和Gauss-FFT法的数值模拟精度和计算效率，分析总结了标准FFT的扩边策略；对本文提出的起伏地形磁异常场快速算法的高效性进行验证.

4.1 正确性验证

模型如图 2所示，其中，模拟区域范围为40 km×40 km×10 km，棱柱形异常体位于模拟区域中心，异常体大小为20 km×20 km×5 km，顶面埋深为1.5 km.模拟区域网格剖分为201×201×101，空间水平采样间隔均为200 m，垂向采样间隔为100 m.异常体磁化率为0.01 SI，背景磁场45000 nT，磁倾角45°，磁偏角5°.

图 2 棱柱体模型示意图 Fig. 2 Sketch of cuboid model

图 3为磁异常场在z=0 m平面的数值解、解析解与绝对误差.图 4、图 5为磁异常梯度张量在z=0 m平面的数值解、解析解与绝对误差.从图 3可以看出，磁异常场的绝对误差最大约为0.3 nT，比磁异常场解析解低三个数量级；从图 4、图 5可以看出磁异常梯度张量的绝对误差最大约为6×10-5nT·m－1，比梯度张量解析解低三个数量级.计算结果表明该算法正确、精度高.

图 3 磁异常场在z=0 m平面的数值解、解析解及绝对误差 Fig. 3 Numerical solution, analytical solution and absolute errors of magnetic anomaly field at z=0 m
图 4 磁梯度张量在z=0 m面的数值解、解析解及绝对误差 Fig. 4 Numerical solution, analytical solution and absolute errors of magnetic gradient tensor at z=0 m
图 5 磁梯度张量在z=0 m面的数值解、解析解及绝对误差 Fig. 5 Numerical solution, analytical solution and absolute error of magnetic gradient tensor at z=0 m
4.2 计算精度与效率

Wu和Tian(2014)首次将Gauss-FFT引入到重磁位场的数值模拟中，与标准的FFT相比，该方法削弱了截断边界的影响，提高了计算精度，但同时降低了计算效率.设计组合棱柱体模型，对基于Gauss-FFT与标准FFT扩边法的数值计算精度与计算效率进行对比，并分析总结了标准FFT扩边系数的选取策略.

引入相对均方根误差(Wu，2016)，该误差统计方式能突出大异常值所占的比重，避免小异常和过零异常造成的误差失真，计算公式为

(22)

其中，rrms为相对均方根误差，Nx, Ny分别为x, y方向的节点数，Bij为磁异常场(磁异常梯度张量)数值解，为磁异常场(磁异常梯度张量)解析解.

设计如图 6所示模型，模型模拟区域为66 km×66 km×10 km.设计了一个位于模拟区域中心的大棱柱形异常体和四个关于原点中心对称的小棱柱形异常体.其中，模型中心异常体大小为9 km×9 km×3 km，顶面埋深为1.5 km.距离模拟区域边界6 km均匀分布四个6 km×6 km×3 km的异常体，其顶面埋深均为1.5 km.棱柱形异常体的磁化率均为0.01 SI，磁倾角为45°，磁偏角为5°，背景磁场为45000 nT.研究区域网格剖分为221×221×101，空间水平采样间隔均为300 m，垂向采样间隔为100 m，观测面为z=0 m平面.

图 6 模型水平面示意图 Fig. 6 Sketch of the horizontal plane of the model

定义扩边系数:

(23)

其中，S为从模拟区域边界向外扩边距离，d为模拟区边缘异常体中心埋深，L为扩边系数，是本文定义的衡量特定规模异常体扩边程度的重要参数，反映了外扩距离与模拟区边缘异常体中心埋深的关系.

图 7为基于标准FFT扩边法的磁异常场及磁梯度张量的相对均方根误差.从图 7中可以看出当扩边系数选取为磁异常场(L=10)和梯度张量(L=4)时，各分量的相对均方根误差基本均在0.5%以下，基本满足实际应用的计算精度需求.由于梯度张量衰减比场更快，在边界截断效应影响相对较少，因此张量的扩边系数相对更小. 图 8为基于不同高斯点的Gauss-FFT法磁异常场及磁梯度张量的相对均方根误差.从图 8中可以看出基于两个高斯点的Gauss-FFT法的相对均方根误差均在2%以上，基于四个高斯点的Gauss-FFT法有效地压制了边界截断效应，相对均方根误差均在0.5%以下，计算精度高.图 7与图 8对比可以看出，标准FFT扩边法(场(L≥10)和梯度张量(L≥4))和基于四个高斯点的Gauss-FFT法的相对均方根误差均在0.5%左右，两种方法的计算精度相近.

图 7 基于标准FFT扩边法的磁异常场及梯度张量相对均方根误差 (a)磁异常场相对均方根误差；(b)磁梯度张量相对均方根误差. Fig. 7 The relative mean square error of the standard FFT method for magnetic anomaly field and gradient tensor (a) Field component error; (b) Gradient tensor error.
图 8 基于Gauss-FFT法的磁异常场及梯度张量相对均方根误差 (a)磁异常场相对均方根误差；(b)磁梯度张量相对均方根误差. Fig. 8 The relative mean square error of the Gauss-FFT method for magnetic anomaly field and gradient tensor (a) Field component error; (b) Gradient tensor error.

图 9为基于标准FFT法不同扩边系数与基于Gauss-FFT法不同高斯点下磁异常场或张量计算时间.表 1为不同扩边系数与不同高斯点下每个磁异常场或张量计算时间与对应的网格大小.测试计算机配置为CPU-Inter Core i5-4590, 主频为3.3 GHz, 内存为16 GB.从图 9中可以看出随着扩边系数的增大，基于标准FFT扩边法的计算时间近似线性增长，基于Gauss-FFT法的计算时间随着高斯点数的增大同样近似线性增长.从表 1可以看出当扩边系数(场的L=10和梯度张量的L=4)时，基于标准FFT扩边法的计算时间分别约为4.79 s、2.36 s，与基于四个高斯点的Gauss-FFT法计算精度相近，计算时间分别约为后者的1/5和1/11.从上述计算时间与精度对比可知，基于标准FFT扩边法(场的L=10和梯度张量的L=4)同时兼顾了计算精度与效率，更加适用于大规模、高效、高精度的磁异常数值模拟.

图 9 磁异常场或张量在不同扩边系数和高斯点下的计算时间 (a)不同扩边系数计算时间；(b)不同高斯点计算时间. Fig. 9 The calculation time of a certain component with different expansion coefficients and different gauss points (a) Calculation time of different boundary expansion factor; (b) The calculation times of different Gaussian point.
表 1 磁异常场或张量在不同扩边系数和高斯点下的计算时间与网格大小 Table 1 The calculation time and grid size of magnetic anomaly field or tensor with different expansion coefficients and different gauss points

图 10、图 11、图 12分别为基于标准FFT扩边法扩边系数为(场的L=10和梯度张量的L=4)时，磁异常场和梯度张量的数值解、解析解及绝对误差.从图 10可以看出，磁异常场的绝对误差最大约为0.5 nT，从图 11、图 12可以看出磁异常梯度张量的绝对误差最大约为10-4nT·m－1，磁异常场与张量的绝对误差均比解析解低两个数量级以上，满足实际计算精度需求.计算结果表明选取扩边系数(场的L=10和梯度张量的L=4)的策略是正确的.

图 10 扩边系数L=10时磁异常场数值解、解析解及绝对误差 Fig. 10 Numerical solution, analytical solution and absolute errors of magnetic anomaly field at L=10
图 11 扩边系数L=4时磁梯度张量数值解、解析解及绝对误差 Fig. 11 Numerical solution, analytical solution and absolute errors of magnetic gradient tensor at L=4
图 12 扩边系数L=4时磁梯度张量数值解、解析解及绝对误差 Fig. 12 Numerical solution, analytical solution and absolute errors of magnetic gradient tensor at L=4

综合考虑计算精度和计算时间两个因素，基于标准FFT扩边法比基于Gauss-FFT法更有利于实现大规模、高效、高精度的磁异常数值模拟，特别在面向实际问题时，标准FFT扩边法可根据需要选择不同的扩边系数，选择性更加灵活，更能兼顾计算精度与计算时间.

4.3 起伏地形条件下磁异常场快速算法验证

起伏地形如图 13所示，地形起伏范围为(-584 m，74 m)，异常体模型采用图 6所示模型，模拟区域网格剖分为221×221×101，观测点个数为221×221.采用每个深度节点zp均对整个区域直接积分累加的算法和快速算法分别计算了11个深度节点、51个深度节点、101个深度节点所在水平面磁异常场，验证起伏地形条件下磁异常场快速算法.

图 13 起伏地形示意图 Fig. 13 A model with topography

图 14为采用51个深度节点插值计算起伏面磁场的数值解、解析解与绝对误差.从图中可以看出磁异常场的绝对误差最大约为0.3 nT，比解析解低两个数量级以上，具有较高的计算精度.

图 14 插值个数n=51时磁异常场数值解、解析解及绝对误差 Fig. 14 Numerical solution, analytical solution and absolute errors of magnetic anomaly field at n=51

图 15为快速算法与直接积分累加算法的串行计算时间对比，其中n为不同深度节点所在水平面个数.从图中可以看出两种不同算法的计算时间都随着n的增大近似线性增长.快速算法计算11个不同深度节点所在平面磁异常的时间为34 s，约为直接积分累加算法计算时间的1/2；快速算法计算101个不同深度节点所在平面磁异常的时间为57 s，约为直接积分累加算法计算时间的1/7，并且随着n的增大，快速算法的优势越明显.计算结果表明在起伏地形条件下磁异常场数值模拟中，该快速算法具有更高的计算效率.

图 15 快速算法与直接积分累加算法计算时间 Fig. 15 The calculate time of fast algorithm and direct integral accumulation algorithm
5 结论

(1) 本文提出一种高效、高精度的空间波数混合域磁异常场积分解三维数值模拟方法.该方法通过沿水平方向进行二维傅里叶变换，把空间域磁位的三维卷积问题转化为不同波数之间相互独立的垂向一维积分问题.保留垂向为空间域，垂向单元可灵活剖分，兼顾了计算精度与计算效率；采用二次形函数描述垂向一维积分单元内磁化强度变化，可得出单元积分的解析表达式，计算精度高、效率高.该方法充分利用一维形函数积分的高效和高精度、快速傅里叶变换的高效性及算法高度并行性, 实现了磁异常场高效、高精度的数值模拟.设计棱柱体模型，模型的解析解与本文方法的数值解对比表明，该方法正确、精度高、效率高.

(2) 设计棱柱体组合模型，模拟区域网格剖分规模为221×221×101，观测点个数为221×221，研究了基于标准FFT扩边法和Gauss-FFT法数值模拟的计算精度与计算效率.模型算例数值试验结果表明：当基于标准FFT扩边系数为(场的L=10、梯度张量的L=4)时，与基于四个高斯点Gauss-FFT法计算精度相近，基于标准FFT扩边法的计算时间分别为4.79 s与2.36 s，计算时间约为Gauss-FFT法的1/5与1/11.综合考虑计算精度与计算效率，基于标准FFT扩边法更有利于实现大规模、高效、高精度的磁异常数值模拟.

(3) 设计起伏地形模型，模拟区域网格剖分为221×221×101，观测点个数为221×221.快速算法串行计算101个不同深度节点所在水平面的时间为57 s，约为直接积分累加算法计算时间的1/7，且随着深度节点所在水平面个数的增加，快速算法的优势越明显.

(4) 本文算法不同波数之间并行度高、不同深度节点正反傅里叶变换并行度高、起伏地形最低点与最高点之间不同深度节点所在平面上磁异常场计算高度并行，文中算例均采用串行算法，计算效率比目前主流算法效率高，如果将串行算法改为并行算法，计算效率将会进一步提高.

附录A

分量积分表达式为

(A1)

其中z′i, z′i+2为第i个单元的上限与下限, m为单元个数.

根据观测点与积分单元的相对位置，需要分三种情况进行讨论.对于第i个单元积分，

(1) 当观测点在该单元上方z＜z′i时，式(A1)第i个单元积分公式为

(A2)

(2) 当观测点在该单元下方z>z′i+2时，式(A1)第i个单元积分公式为

(A3)

(3) 当观测点在该单元内部z′i≤z≤z′i+2时，式(A1)第i个单元积分公式为

(A4)

对于第一种情况式(A2)，令γ=k，式(A2)可表示为

(A5)

令，则式(A5)可表示为

(A6)

对于第二种情况式(A3)，令γ=－k，式(A3)可表示为

(A7)

令，式(A7)可表示为

(A8)

对于第三种情况式(A4)，在(zi, z)积分区域令γ=－k，；在(z, z′i+2)积分区域令γ=k，，则式(A4)可表示为

(A9)

对式(A6)或式(A8)单元积分，采用二次形函数来表示插值函数，插值公式为

(A10)

式中，N1、N2、N3为形函数，表达式为

(A11)

(A12)

(A13)

其中，

(A14)

(A15)

单元积分可表示为

(A16)

式中，

(A17)

(A18)

(A19)

式(A17)—(A19)中单元积分部分推导过程如下：

(A20)

(A21)

(A22)

其中，

(A23)

(A24)

(A25)

式(A9)的单元形函数展开与式(A6)类似，不同在于积分的上下限不同.

利用每个单元的插值函数，可将所求积分表示为

(A26)

其中z表示计算深度节点水平面位置，m表示积分单元个数.

同理，分量与上述计算过程相同.

References
Bhattacharyya B K. 1964. Magnetic anomalies due to prism-shaped bodies with arbitrary polarization. Geophysics, 29(4): 517-531. DOI:10.1190/1.1439386
Blakely R J. 1996. Potential theory in gravity and magnetic applications: Cambridge University Press.
Chai Y P. 1988. Algorithm investigation of DFT of potential field. Acta Geophysics Sinica (in Chinese), 31(2): 211-224.
Chai Y P. 1997. Shift sampling theory of Fourier transform computation. Science in China (Series E), 40(1): 21-27. DOI:10.1007/BF02916587
Chakravarthi V, Raghuram H M, Singh S B. 2002. 3-D forward gravity modeling of basement interfaces above which the density contrast varies continuously with depth. Computers & Geosciences, 28(1): 53-57.
Chen Z X, Meng X H, Liu G F, et al. 2012. The GPU-based parallel calculation of gravity and magnetic anomalies for 3D arbitrary bodies. Geophysical and Geochemical Exploration (in Chinese), 36(1): 117-121.
Feng R. 1986. A computation method of potential field with three-dimensional density and magnetization distributions. Acta Geophysics Sinica (in Chinese), 29(4): 399-406.
Furness P. 1994. A physical approach to computing magnetic fields. Geophysical Prospecting, 42(5): 405-416. DOI:10.1111/j.1365-2478.1994.tb00218.x
García-Abdeslem J. 2005. The gravitational attraction of a right rectangular prism with density varying with depth following a cubic polynomial. Geophysics, 70(6): J39-J42. DOI:10.1190/1.2122413
Guan Z N. 2005. Geomagnetic Field and Magnetic Exploration (in Chinese). Beijing: Geological Publishing House.
Guo Z H, Guan Z N, Xiong S Q. 2004. Cuboid ΔT and its gradient forward theoretical expressions without analytic odd points. Chinese Journal of Geophysics (in Chinese), 47(6): 1131-1138.
Guptasarma D, Singh B. 1999. New scheme for computing the magnetic field resulting from a uniformly magnetized arbitrary polyhedron. Geophysics, 64(1): 70-74.
Hu S, Tang J, Ren Z, et al. 2019. Multiple underwater objects localization with magnetic gradiometry. IEEE Geoscience and Remote Sensing Letters, 16: 296-300. DOI:10.1109/LGRS.2018.2870839
Ivan M. 1996. Optimum expression for computation of the magnetic field of a homogeneous polyhedral body. Geophysical Prospecting, 44(2): 279-288. DOI:10.1111/j.1365-2478.1996.tb00150.x
Kwok Y K. 1991. Gravity gradient tensors due to a polyhedron with polygonal facets. Geophysical Prospecting, 39(3): 435-443. DOI:10.1111/j.1365-2478.1991.tb00320.x
Luo Y, Yao C L. 2007. Theoretical study on cuboid magnetic field and gradient expression without singular point. Oil Geophysical Prospecting (in Chinese), 42(6): 714-719.
Nagy D, Papp G, Benedek J. 2000. The gravitational potential and its derivatives for the prism. Journal of Geodesy, 74(7-8): 552-560. DOI:10.1007/s001900000116
Pedersen L B. 1978. Wavenumber domain expressions for potential fields from arbitrary 2-, 21/2-, and 3-dimensional bodies. Geophysics, 43(3): 626-630.
Pedersen L B. 1985. The gravity and magnetic fields from ellipsoidal bodies in the wavenumber domain. Geophysical Prospecting, 33(2): 263-281. DOI:10.1111/j.1365-2478.1985.tb00434.x
Plouff D. 1976. Gravity and magnetic fields of polygonal prisms and application to magnetic terrain corrections. Geophysics, 41(4): 727-741. DOI:10.1190/1.1440645
Pohánka V. 1998. Optimum expression for computation of the gravity field of a polyhedral body with linearly increasing density. Geophysical Prospecting, 46(4): 391-404. DOI:10.1046/j.1365-2478.1998.960335.x
Ren Z Y, Tang J T, Kalscheuer T, et al. 2017. Fast 3-D large-scale gravity and magnetic modeling using unstructured gridsand an adaptive multilevel fast multipole method. Journal of Geophysical Research:Solid Earth, 122(1): 79-109. DOI:10.1002/2016JB012987
Singh B, Guptasarma D. 2001. New method for fast computation of gravity and magnetic anomalies from arbitrary polyhedra. Geophysics, 66(2): 521-526. DOI:10.1190/1.1444942
Singh S K, Sabina F J. 1978. Magnetic anomaly due to a vertical right circular cylinder with arbitrary polarization. Geophysics, 43(1): 173-178. DOI:10.1190/1.1440818
Tontini F C. 2005. Magnetic-anomaly Fourier spectrum of a 3D Gaussian source. Geophysics, 70(1): L1-L5.
Tontini F C, Cocchi L, Carmisciano C. 2009. Rapid 3-D forward model of potential fields with application to the Palinuro Seamount magnetic anomaly (southern Tyrrhenian Sea, Italy). Journal of Geophysical Research:Solid Earth, 114(B2): B02103. DOI:10.1029/2008JB005907
Wu L Y, Tian G. 2014. High-precision Fourier forward modeling of potential fields. Geophysics, 79(5): G59-G68. DOI:10.1190/geo2014-0039.1
Wu L Y. 2016. Efficient modelling of gravity effects due to topographic masses using the Gauss-FFT method. Geophysical Journal International, 205(1): 160-178. DOI:10.1093/gji/ggw010
Wu X Z. 1981. Computation of spectrum of potential field due to 3-dimensional bodies (Homogeneous models). Acta Geophysics Sinica (in Chinese), 24(3): 336-348.
Wu X Z. 1983. The computation of spectrum of potential field due to 3-D arbitrary bodies with physical Parameters varying with depth. Acta Geophysics Sinica (in Chinese), 26(2): 177-187.
Xiong G C. 1984. Some problems about the 3-D fourier transforms of the gravity and magnetic fields. Acta Geophysics Sinica (in Chinese), 27(1): 103-109.
Yao C L, Hao T Y, Guan Z N, et al. 2003. High-speed Computation and efficient storage in 3-D gravity and magnetic inversion based on genetic algorithms. Chinese Journal of Geophysics (in Chinese), 46(2): 252-258.
柴玉璞. 1988. 位场DFT算法研究. 地球物理学报, 31(2): 211-224. DOI:10.3321/j.issn:0001-5733.1988.02.010
柴玉璞. 1996. Fourier变换数值计算的偏移抽样理论. 中国科学(E辑), 26(5): 450-456.
陈召曦, 孟小红, 刘国峰, 等. 2012. 基于GPU的任意三维复杂形体重磁异常快速计算. 物探与化探, 36(1): 117-121.
冯锐. 1986. 三维物性分布的位场计算. 地球物理学报, 29(4): 399-406. DOI:10.3321/j.issn:0001-5733.1986.04.010
管志宁. 2005. 地磁场与磁力勘探. 北京: 地质出版社.
郭志宏, 管志宁, 熊盛青. 2004. 长方体ΔT场及其梯度场无解析奇点理论表达式. 地球物理学报, 47(6): 1131-1138. DOI:10.3321/j.issn:0001-5733.2004.06.029
骆遥, 姚长利. 2007. 长方体磁场及其梯度无解析奇点表达式理论研究. 石油地球物理勘探, 42(6): 714-719. DOI:10.3321/j.issn:1000-7210.2007.06.019
吴宣志. 1981. 三度体(均质模型)位场波谱的正演计算. 地球物理学报, 24(3): 336-348. DOI:10.3321/j.issn:0001-5733.1981.03.012
吴宣志. 1983. 三度体(物性随深度变化模型)位场波谱的正演计算. 地球物理学报, 26(2): 177-187. DOI:10.3321/j.issn:0001-5733.1983.02.009
熊光楚. 1984. 重、磁场三维傅里叶变换的若干问题. 地球物理学报, 27(1): 103-109. DOI:10.3321/j.issn:0001-5733.1984.01.011
姚长利, 郝天珧, 管志宁, 等. 2003. 重磁遗传算法三维反演中高速计算及有效存储方法技术. 地球物理学报, 46(2): 252-258. DOI:10.3321/j.issn:0001-5733.2003.02.020
空间波数混合域磁异常场积分解三维数值模拟
李昆, 戴世坤, 陈轻蕊, 张钱江, 赵东东, 王顺国, 凌嘉宣