复杂形体重磁位场三维高效高精度数值模拟

2020-10-17周印明戴世坤凌嘉宣胡晓颖

石油地球物理勘探 2020年5期

周印明戴世坤李昆凌嘉宣胡晓颖熊彬

（①中南大学地球科学与信息物理学院，湖南长沙410083；②中南大学有色金属成矿预测与地质环境监测教育部重点实验室，湖南长沙410083；③东方地球物理公司综合物化探处，河北涿州072751；④桂林理工大学地球科学学院，广西桂林541006）

0 引言

重磁位场正演方法主要分为空间域方法和频率域方法。空间域方法又可以分为解析法和数值法。解析法是通过重磁异常场解析表达式直接求取空间任意点的精确异常场，其优点是原理简单、计算结果精确。解析法研究早期侧重于简单、规则模型的重磁场解析表达式的推导，如水平薄板、直立棱柱体、直立圆柱体等［1-4］。后来学者们逐步对均匀介质的多面体模型［5-9］和简单物性变化的任意三度体模型的重磁场解析式推导［10-11］、解析积分奇异点处理［9，12］、边界场连续性问题［12］等开展了大量的研究。一般来说，解析法表达式比较复杂、适应性较差、计算量较大。数值法主要通过有限差分法［13］、有限体积法［14］、有限单元法［15-16］求解重磁异常满足的泊松方程，但是当计算大量位置点的异常场时，耗时较长。

频率域方法的原理是采用数值方法对场源在某条测线（或平面／三维空间）产生的重磁异常频谱（傅里叶变换解析表达式）进行计算，并通过反傅里叶变换得到异常场［17］。其优点是表达式简单，计算效率相对较高，但是傅里叶变换的截断效应限制了该方法的普遍适用性。Bhattacharyya［3］推导了任意磁化方向的长方体磁场频谱表达式；Pedersen［18］给出了直立圆柱体和多面体模型的重磁异常频率域表达式；熊光楚［19］给出了长方体模型重力位的三维傅里叶变换式；后来一些学者对物性连续变化模型的频率域解析式［20-24］、Parker模型频率域表达式［25-26］、偏移抽样方法提高反傅里叶变换数值精度［27］开展了大量研究。Tontini等［28］实现了重磁异常三维全空间频率域正演，并研究了FFT 扩边法与误差的关系；Wu等［29-30］利用高斯FFT 提高了反傅里叶变换的数值精度，降低了由于FFT 引起的强制周期化、边界震荡效应等影响，并研究了地形对重磁异常的影响，分析了物性连续变化模型的梯度场；Ren等［31］提出一种基于非结构化网格的自适应快速多极子重磁场及其梯度张量正演方法，借助并行计算实现了快速、高精度的三维正演。

目前，随着计算机技术的不断发展和重磁异常正演算法的不断改进，频率域方法以其简捷、准确和高效等优势在数值模拟中的应用越来越广泛。但是，对于面向实际应用的超大规模复杂三维模型数值模拟问题，由于网格剖分单元数目剧增，计算量和存储量巨大，常规的空间域和频率域重磁异常正演方法已很难同时满足精度与效率上的要求。为此，本文提出一种重磁位场三维数值模拟方法，充分利用一维形函数积分的高精度性、不同波数之间高度并行性及傅里叶变换的高效性，实现高效、高精度、大规模的重磁位场数值模拟。

1 理论与方法

弱磁性体的磁位V 及重力位φ 的空间域积分表达分别为［32］

式中：M 为空间域磁化强度矢量；ρ为剩余密度；μ为真空磁导率，取值4π×10－7H／m；D 为万有引力常量，取值6.674×10－11m3kg－1s－2；（x，y，z）为观测点坐标；（x′，y′，z′）为异常点坐标；为哈密顿算符；G 为格林函数。

对式（1）沿x、y 方向做二维傅里叶变换，分别得到磁力位和重力位一维积分公式

式（2）的推导需利用波数域格林函数，其表达式为

磁位和重力位的一阶导数分别为

对应的二阶导数分别为

式中下标“B”和“g”分别代表磁力场和重力场。

利用傅里叶变换的微分性质，对式（4）和式（5）做二维傅里叶变换，得到波数域重磁异常场和张量的表达式分别为

式（2）为重磁位垂向一维积分表达式。从该表达式可以看出，不同波数一维空间域积分是相互独立的，可以并行计算；该一维积分深度方向可离散为多个单元积分之和，每个积分单元采用形函数法插值计算，得到单元积分的解析表达式，可提高计算精度和效率。

2 基于形函数的一维积分

本文采用基于二次插值的形函数法［31］计算傅里叶变换后的一维积分。将一维积分沿垂直方向剖分为m 个单元，采用二次形函数描述物性参数变化。从式（2）可以看出，重力位与磁位的积分表达式相似，只是系数不同；重力异常与磁异常积分中，不同方向的积分表达式也相似，只是系数不同。因此，下面以x 方向磁位积分为例，介绍形函数法计算一维积分的具体过程。

对积分区域进行离散

式中i表示积分单元。

对磁化强度在节点处采用二次形函数插值

式中：N1、N2、N3为二次形函数；i1、i2、i3分别为积分单元i的3个插值节点。

将式（9）代入式（8）中，整理得

由式（10）离散后的一维单元积分可推导出解析表达式，该表达式计算精度较高，接近于解析解。对于物性参数满足二次变化的复杂形体，垂向单元剖分间隔可根据需要灵活选择。

3 算法流程

本文按照如下流程进行重磁位场数值模拟计算。

（1）网格剖分。给定计算区域，确定沿x、y、z方向的剖分单元个数C1、C2、C3。水平方向均匀剖分，垂直方向可任意剖分，得到三维离散坐标点（xm，yn，zl），给定计算的观测平面z。

（2）模型参数设定。给所有的剖分单元赋予剩余密度或磁化率值。

（3）离散波数计算。由FFT 法的波数计算公式得出水平离散坐标（xm，yn）的离散波数（kx，ky）。

（4）重磁空间域二维傅里叶变换。对空间域重磁位场剩余密度或者磁化强度进行水平二维离散傅里叶变换。

（5）积分计算。采用形函数法计算垂向一维积分，得到波数谱。

（6）重磁波数谱二维傅里叶反变换。利用FFT算法，对波数谱场或张量进行二维离散傅里叶反变换，得到空间域重磁位场或张量。

4 算例分析

4.1 正确性验证

设计一个棱柱体模型，采用基于四个高斯点的高斯—FFT法［28］计算重磁异常场及张量的数值解，并与模型解析解［32］进行对比，验证该算法的正确性。

模型如图1所示。异常体大小为200m×200m×200m，顶面埋深h为200m。观测面为z＝0的水平面，模拟区域为800m×800m×400m，观测点个数为201×201。棱柱异常体与模拟区域水平中心点重合。模拟区域网格节点个数为201×201×201，空间采样间隔均为4m。设定异常体的磁化强度为0.01A／m，背景磁场为4500n T，磁倾角为45°，磁偏角为5°；设定其剩余密度为1000kg／m3。

图1 棱柱体模型示意图

图2为磁异常场在z＝0平面的数值解、解析解与绝对误差。可以看出，磁异常场的绝对误差最大约为5×10－3n T。图3为磁异常梯度张量在z＝0平面的数值解、解析解与绝对误差，可以看出，磁异常梯度张量的绝对误差最大约为5×10－5n T／m。磁异常场和磁异常梯度张量的绝对误差均比相应解析解小三个数量级。

图4为重力异常场在z＝0平面的数值解、解析解与绝对误差。可以看出，重力异常场的绝对误差最大约为2×10－4mGal。图5为重力异常梯度张量在z＝0平面的数值解、解析解与绝对误差，可以看出，重力异常梯度张量的绝对误差最大约为0.02E。同样地，重力异常场和重力异常梯度张量的绝对误差值也均比相应解析解低三个数量级。重磁场数值模拟结果表明，本文算法正确，且精度高。

图2 z＝0平面磁异常场的数值解（左）、解析解（中）及绝对误差（右）

4.2 复杂形体数值模拟适用性研究

对波数域的一维积分采用二次形函数法可得到每个单元的近似解析解，该方法计算精度较高，适用于复杂介质的数值模拟计算。以磁异常数值模拟为例设计一个复杂形体模型。异常体的水平方向磁化率不变，垂直方向磁化率呈二次变化。垂直方向分别采用均匀剖分与形函数插值两种方法获得数值解。通过数值解与解析解［25］的对比分析，可验证本文方法对于复杂模型的适用性。为了研究统计整个观测面的误差情况，引入相对均方根误差［30］

图3 z＝0平面磁梯度张量的数值解（左）、解析解（中）及绝对误差（右）

图4 z＝0平面重力异常场的数值解（左）、解析解（中）及绝对误差（右）

式中：P 和N 分别是x 和y 方向的节点数；Bij是磁场（张量）数值解；^Bij是磁场（张量）解析解。该误差统计方式能突出大异常值所占的比重。设计两个不同模型，水平方向磁化率均不变，垂向磁化率从0递增到0.02SI，其中模型1（图6）垂直方向磁化率变化为0.02／60002（z－2000）2（2000≤z≤8000），模型2（图7）磁化率变化为0.02／20002（z－4000）2（4000≤z≤6000），水平方向采样间距均为100m。两个模型的背景磁感应强度均为45000n T，磁倾角为45°，磁偏角为5°。

图8和图9分别为模型1磁场分量和磁张量随着采样间距变化的数值解与解析解的相对均方根误差变化图。可以看出，两种方法在垂直采样间距为25m 时的计算精度相近，在积分单元内磁化率均匀剖分的误差随着采样间距的增大而增大，而形函数二次插值的计算精度基本保持不变。

图10和图11分别为模型2磁场分量和磁张量随着采样间距变化的数值解与解析解的相对均方根误差变化图。

与模型算例1计算结果对比可以看出，两个模型在两种剖分方法下的误差变化趋势基本一致。但是与模型1相比，模型2中单元内磁化率均匀剖分的计算精度有所降低，这是由于模型2中模型介质磁化率的变化更剧烈，而单元内磁化率形函数二次插值的计算精度变化较小。这再次验证了该方法对于复杂介质模型具有较好的适用性。

对于模型区域剖分网格单元数为100×100×100、观测点数为101×101的三维正演计算，传统空间域累加算法耗时约800s，波数域Parker公式计算结果约耗时8.4s，本文算法约耗时5.1s，加速比分别为156.8与1.65，验证了本算法的高效性。