CN115292973A - 一种任意采样的空间波数域三维磁场数值模拟方法及系统 - Google Patents

一种任意采样的空间波数域三维磁场数值模拟方法及系统 Download PDF

Info

Publication number
CN115292973A
CN115292973A CN202211225001.7A CN202211225001A CN115292973A CN 115292973 A CN115292973 A CN 115292973A CN 202211225001 A CN202211225001 A CN 202211225001A CN 115292973 A CN115292973 A CN 115292973A
Authority
CN
China
Prior art keywords
domain
dimensional
magnetic
magnetic field
wave number
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202211225001.7A
Other languages
English (en)
Other versions
CN115292973B (zh
Inventor
戴世坤
张莹
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN202211225001.7A priority Critical patent/CN115292973B/zh
Publication of CN115292973A publication Critical patent/CN115292973A/zh
Application granted granted Critical
Publication of CN115292973B publication Critical patent/CN115292973B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/141—Discrete Fourier transforms
    • G06F17/142—Fast Fourier transforms, e.g. using a Cooley-Tukey type algorithm
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Operations Research (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Measuring Magnetic Variables (AREA)

Abstract

本发明涉及磁法勘探技术领域，公开了一种任意采样的空间波数域三维磁场数值模拟方法及系统，包括：构建目标区域的三维目标模型；基于磁化强度构建磁化强度与空间域磁场异常场磁位的三维泊松方程，并通过任意采样的二维傅里叶正变换将三维泊松方程转换为空间波数混合域的一维常微分方程，求解得到波数域异常场磁位；根据波数域异常场磁位与波数域异常场磁场强度的关系求解得到异常场磁场强度，对波数域异常场磁场强度做任意采样的二维傅里叶反变换，得到空间域异常场磁场强度；通过空间域异常场磁场强度与空间域磁感应强度的关系，得到空间域磁感应强度；本发明解决了现有的傅里叶变换算法对于振荡谱变换的效果不理想的问题，提高了计算精度。

Description

一种任意采样的空间波数域三维磁场数值模拟方法及系统
技术领域
本发明涉及磁法勘探技术领域，尤其涉及一种任意采样的空间波数域三维磁场数值模拟方法及系统。
背景技术
磁法勘探是地球物理勘探的重要手段之一。高效率、高精度的磁异常数值模拟对于磁测资料处理、精细的反演成像具有重要作用。
频率域数值模拟方法在磁测数值模拟中应用非常广泛，而基于傅里叶变换的频率域方法的精度和效率受到所采用的傅里叶变换方法的制约。标准快速傅里叶变换存在截断效应，影响数据处理效果。文献（Wu L Y，Tian G. 2014. High-precision Fourierforward modeling of potential fields. Geophysics， 79(5) G59–G68.）为解决快速傅里叶变换的截断效应，在移位采样方案基础上提出Gauss-FFT方法，对于计算精度大幅提高，并削弱了截断效应影响，但是其仍需要等间隔采样，且内存需求高于快速傅里叶变换。为实现非均匀采样，相关学者在兼顾快速傅里叶变换的快速性和非均匀采样的灵活性的基础上，提出一系列非均匀采样傅里叶变换方法。其中基于最小二乘误差插值的NUFFT算法不考虑指数项的快速振荡，因此对于振荡谱变换的效果不理想。
发明内容
本发明提供了一种任意采样的空间波数域三维磁场数值模拟方法及系统，以解决现有的数值模拟算法对于振荡谱变换的效果不理想的问题。
为了实现上述目的，本发明通过如下的技术方案来实现：
第一方面，本发明提供一种任意采样的空间波数域三维磁场数值模拟方法，包括：
构建含有异常体的目标区域三维目标模型，对所述三维目标模型进行非均匀剖分得到一系列节点，并根据磁化率分布数据对每个节点进行磁化率赋值，得到每个节点的磁化率，其中，包含异常体的节点根据异常体的磁化率分布数据进行赋值，不包含异常体的节点赋值为0，其中，非均匀剖分包括以下任一方式：
方式一：对预设的第一区域采取非均匀剖分，进行加密；其中，所述第一区域满足以下公式：
Figure DEST_PATH_IMAGE001
；
Figure DEST_PATH_IMAGE002
为所述第一区域对应节点的剩余密度，
Figure DEST_PATH_IMAGE003
为第一区域对应节点周围的第j个节点的剩余密度，
Figure DEST_PATH_IMAGE004
为第一区域周围节点数量；
Figure DEST_PATH_IMAGE005
为权重，取值范围为（0,1）；
方式二：对预设的第二区域进行稀疏采样，其中，所述第二区域满足以下公式：
Figure DEST_PATH_IMAGE006
；
Figure DEST_PATH_IMAGE007
为所述第二区域对应节点的剩余密度，
Figure DEST_PATH_IMAGE008
为第二区域对应节点周围的第j个节点的剩余密度，
Figure 544125DEST_PATH_IMAGE004
为第一区域周围节点数量；
Figure 523582DEST_PATH_IMAGE005
为权重，取值范围为（0,1）；
根据地球主磁场模型计算每个节点处的地球主磁场强度，并将所述地球主磁场强度作为空间域背景场磁场强度计算空间域总场磁场强度，根据空间域总场磁场强度和磁化率之间的关系得到磁化强度的计算模型，通过磁化强度的计算模型得到磁化强度；
基于所述磁化强度构建磁化强度与空间域磁场异常场磁位的三维泊松方程，并通过任意采样的二维傅里叶正变换将所述三维泊松方程转换为空间波数混合域的一维常微分方程，求解所述一维常微分方程得到波数域异常场磁位；
根据所述波数域异常场磁位与波数域异常场磁场强度的关系构建计算方程，并对方程求解得到波数域异常场磁场强度，对所述波数域异常场磁场强度做任意采样的二维傅里叶反变换，得到空间域异常场磁场强度；
通过所述空间域异常场磁场强度与空间域磁感应强度的关系，得到空间域磁感应强度。
可选的，所述磁化强度的计算模型为：
Figure DEST_PATH_IMAGE009
；
其中，M表示磁化强度，
Figure DEST_PATH_IMAGE010
表示磁化率，H表示空间域总场磁场强度，空间域总场磁场强度为空间域背景场磁场强度和空间域异常场磁场强度之和；
可选的，所述三维泊松方程为：
Figure DEST_PATH_IMAGE011
；
其中，
Figure DEST_PATH_IMAGE012
表示空间域磁场异常场磁位，M表示磁化强度，式中，
Figure DEST_PATH_IMAGE013
，
Figure DEST_PATH_IMAGE014
，i，j，k分别为x，y，z方向的单位向量。
上式展开为
Figure DEST_PATH_IMAGE015
；
其中，
Figure DEST_PATH_IMAGE016
、
Figure DEST_PATH_IMAGE017
、
Figure DEST_PATH_IMAGE018
分别为磁化强度M在x，y，z方向的分量，
Figure DEST_PATH_IMAGE019
为偏导数符号。
可选的，所述任意采样的二维傅里叶正变换公式为：
Figure DEST_PATH_IMAGE020
；
其中，
Figure DEST_PATH_IMAGE021
表示x方向的波数，
Figure DEST_PATH_IMAGE022
表示y方向上的波数，
Figure DEST_PATH_IMAGE023
表示空间域函数，
Figure DEST_PATH_IMAGE024
表示波数谱；
所述任意采样的二维傅里叶正变换公式通过两次一维傅里叶正变换得到，两次一维傅里叶正变换分别为；
对
Figure 479818DEST_PATH_IMAGE023
进行x方向一维傅里叶正变换，变换式为：
Figure DEST_PATH_IMAGE025
；
其中，x，y表示两个相互垂直的方向；
Figure 609448DEST_PATH_IMAGE021
表示x方向的波数，
Figure DEST_PATH_IMAGE026
表示空间域函数，
Figure DEST_PATH_IMAGE027
为对
Figure 151435DEST_PATH_IMAGE026
做x方向做一维傅里叶变换后的波数谱；
对
Figure 646002DEST_PATH_IMAGE027
进行y方向一维傅里叶正变换，变换式为：
Figure DEST_PATH_IMAGE028
；
其中，
Figure 991664DEST_PATH_IMAGE022
表示y方向上的波数，
Figure DEST_PATH_IMAGE029
为对
Figure 659405DEST_PATH_IMAGE026
做二维傅里叶变换后的波数谱。
可选的，所述一维傅里叶正变换具体变换方法为：
设连续一维傅里叶正变换分别表示为：
Figure DEST_PATH_IMAGE030
；
其中，
Figure 282279DEST_PATH_IMAGE021
表示波数，
Figure DEST_PATH_IMAGE031
表示空间域函数，
Figure DEST_PATH_IMAGE032
表示波数谱；
对上述连续一维傅里叶正变换进行离散得到：
Figure DEST_PATH_IMAGE033
；
其中，
Figure DEST_PATH_IMAGE034
表示单元个数，
Figure DEST_PATH_IMAGE035
表示第j个单元，i为虚数，
Figure 839424DEST_PATH_IMAGE021
表示x方向上的波数；
用二次形函数对
Figure DEST_PATH_IMAGE036
进行插值：
单元内采用二次插值形函数拟合时，设任一单元内部三个节点坐标分别为
Figure DEST_PATH_IMAGE037
，
Figure DEST_PATH_IMAGE038
为中点，满足
Figure DEST_PATH_IMAGE039
，每个节点上值分别为
Figure DEST_PATH_IMAGE040
，
Figure DEST_PATH_IMAGE041
采用二次形函数进行表示可得：
Figure DEST_PATH_IMAGE042
；
其中，N1、N2、N3表示二次插值函数，分别为：
Figure DEST_PATH_IMAGE043
；
上式写为：
Figure DEST_PATH_IMAGE044
Figure DEST_PATH_IMAGE045
；
令
Figure DEST_PATH_IMAGE046
，
Figure DEST_PATH_IMAGE047
，
Figure DEST_PATH_IMAGE048
为单元内傅里叶变换节点系数，即W1、W2、W3分别表示每个节点对应的傅里叶变换系数，则上式写为：
Figure DEST_PATH_IMAGE049
；
当波数
Figure DEST_PATH_IMAGE050
不为0时，将N 1、N 2、N 3代入
Figure DEST_PATH_IMAGE051
中，可得单元内傅里叶变换节点系数W1、W2、W3分别为：
Figure DEST_PATH_IMAGE052
Figure DEST_PATH_IMAGE053
Figure DEST_PATH_IMAGE054
Figure DEST_PATH_IMAGE055
；
Figure DEST_PATH_IMAGE056
Figure DEST_PATH_IMAGE057
Figure DEST_PATH_IMAGE058
积分核函数都包括
Figure DEST_PATH_IMAGE059
，其在
Figure DEST_PATH_IMAGE060
上单元积分解析解为：
Figure DEST_PATH_IMAGE061
Figure DEST_PATH_IMAGE062
Figure DEST_PATH_IMAGE063
；
Figure DEST_PATH_IMAGE064
因此，可得
Figure 382969DEST_PATH_IMAGE050
不为0时的
Figure DEST_PATH_IMAGE065
半解析解为：
Figure DEST_PATH_IMAGE066
Figure DEST_PATH_IMAGE067
Figure DEST_PATH_IMAGE068
；
Figure DEST_PATH_IMAGE069
Figure DEST_PATH_IMAGE070
当波数
Figure 394746DEST_PATH_IMAGE050
为0时，
Figure DEST_PATH_IMAGE071
，
Figure DEST_PATH_IMAGE072
，
Figure DEST_PATH_IMAGE073
，
Figure DEST_PATH_IMAGE074
分别表示波数为0时的傅里叶变换系数，进行简单积分即可求得零波数下傅里叶变换节点系数为：
Figure DEST_PATH_IMAGE075
；
将不同单元解析表达式进行累加，可得最终一维傅里叶正变换结果。
可选的，所述一维常微分方程通过将所述三维泊松方程进行水平方向二维傅里叶正变换得到，所述一维常微分方程为：
Figure DEST_PATH_IMAGE076
；
其中，
Figure DEST_PATH_IMAGE077
表示波数域异常场磁位，
Figure DEST_PATH_IMAGE078
、
Figure DEST_PATH_IMAGE079
、
Figure DEST_PATH_IMAGE080
表示波数域磁化强度
Figure DEST_PATH_IMAGE081
的x分量、y分量和z分量，
Figure 452219DEST_PATH_IMAGE021
、
Figure 163954DEST_PATH_IMAGE022
分别表示x、y方向的波数，
Figure DEST_PATH_IMAGE082
为偏导数符号；
所述求解所述一维常微分方程得到波数域异常场磁位，包括：
在笛卡尔坐标系下，取Z轴垂直向下为正向，计算区域取水平地面为上边界Z min ，取地下离异常体足够远处为下边界Z max ，其上下边界条件满足：
上边界：
Figure DEST_PATH_IMAGE083
；
下边界：
Figure DEST_PATH_IMAGE084
；
其中，
Figure DEST_PATH_IMAGE085
；
联立一维常微分方程和上下边界，得到：
Figure DEST_PATH_IMAGE086
；
运用变分法，得到与边值问题等价的变分问题：
Figure DEST_PATH_IMAGE087
沿z方向进行单元剖分，在每个单元内采用二次插值函数，得到每个节点处的波数域异常场磁位
Figure 500520DEST_PATH_IMAGE077
。
可选的，所述根据所述波数域异常场磁位与波数域异常场磁场强度的关系构建计算方程，包括：
波数域异常场磁位
Figure 57534DEST_PATH_IMAGE077
与波数域异常场磁场强度
Figure DEST_PATH_IMAGE088
满足如下关系式：
Figure DEST_PATH_IMAGE089
；
其中，i为虚数。
可选的，所述二维傅里叶反变换公式为：
Figure DEST_PATH_IMAGE090
；
其中，
Figure 779634DEST_PATH_IMAGE021
、
Figure 662270DEST_PATH_IMAGE022
表示波数，
Figure 63296DEST_PATH_IMAGE023
为空间域函数，
Figure 345372DEST_PATH_IMAGE024
表示波数谱；
所述空间域异常场磁场强度与空间域磁感应强度的关系为：
Figure DEST_PATH_IMAGE091
；
其中，
Figure DEST_PATH_IMAGE092
表示介质的绝对磁导率，单位为H/m，
Figure 390820DEST_PATH_IMAGE092
与
Figure DEST_PATH_IMAGE093
的关系满足下式：
Figure DEST_PATH_IMAGE094
；
其中，
Figure DEST_PATH_IMAGE095
表示真空中磁导率，
Figure DEST_PATH_IMAGE096
H/m。
第二方面，本申请实施例提供一种任意采样的空间波数域三维磁场数值模拟系统，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述第一方面中任一所述方法的步骤。
有益效果：
本发明提供的任意采样的空间波数域三维磁场数值模拟方法，通过任意采样的二维傅里叶变换，将三维问题降至一维，只保留z方向在空间域，运用一维有限单元法，单元内部采用形函数二次插值，从而对微分方程进行求解，之后再通过任意采样的傅里叶反变换回空间域，极大的提高了计算精度和计算效率，算法并行度好，占用内存少。
此外，本发明中应用的基于形函数二次插值的任意采样傅里叶变换方法（AS-FT），能够提前计算好傅里叶变换系数，根据场和谱的分布灵活设定采样间隔，根据需求适量稀疏和加密采样点，同时在积分区间内能够得到傅里叶振荡算子
Figure DEST_PATH_IMAGE097
的解析解，兼顾精度和效率，并且不存在截断效应，将该傅里叶变换方法应用于偏微分方程解法中，能够完美解决边界问题，且计算效率较高。
附图说明
图1为本发明优选实施例的任意采样的空间波数域三维磁场数值模拟方法的流程图；
图2为本发明优选实施例的任意采样的空间波数域三维磁场数值模拟方法的三维目标模型结构示意图；
图3为本发明优选实施例的任意采样的空间波数域三维磁场数值模拟方法的剖分示意图之一；
图4为本发明优选实施例的任意采样的空间波数域三维磁场数值模拟方法的剖分示意图之二；
图5为本发明优选实施例的任意采样的空间波数域三维磁场数值模拟方法的单元节点结构示意图；
图6为本发明优选实施例的任意采样的空间波数域三维磁场数值模拟方法的边界条件示意图。
图7为本发明优选实施例的任意采样的空间波数域三维磁场数值模拟方法的结果与解析解对比图。
具体实施方式
下面对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
除非另作定义，本发明中使用的技术术语或者科学术语应当为本发明所属领域内具有一般技能的人士所理解的通常意义。本发明中使用的“第一”、“第二”以及类似的词语并不表示任何顺序、数量或者重要性，而只是用来区分不同的组成部分。同样，“一个”或者“一”等类似词语也不表示数量限制，而是表示存在至少一个。“连接”或者“相连”等类似的词语并非限定于物理的或者机械的连接，而是可以包括电性的连接，不管是直接的还是间接的。“上”、“下”、“左”、“右”等仅用于表示相对位置关系，当被描述对象的绝对位置改变后，则该相对位置关系也相应地改变。
应理解，本申请的一种任意采样的空间波数域三维磁场数值模拟方法可以应用于磁法勘探中，例如地下矿脉勘探，石油、天然气勘查，地质构造推断等等，此处只做示例，不做限定。
实施例1，请参见图1：
本申请实施例提供一种任意采样的空间波数域三维磁场数值模拟方法，包括：
构建含有异常体的目标区域三维目标模型，对所述三维目标模型进行剖分得到一系列节点，并根据磁化率分布数据对每个节点进行磁化率赋值，得到每个节点的磁化率，其中，包含异常体的节点根据异常体的磁化率分布数据进行赋值，不包含异常体的节点赋值为0；
根据地球主磁场模型计算每个节点处的地球主磁场强度，并将所述地球主磁场强度作为空间域背景场磁场强度计算空间域总场磁场强度，根据空间域总场磁场强度和磁化率之间的关系得到磁化强度的计算模型，通过磁化强度的计算模型得到磁化强度；
基于所述磁化强度构建磁化强度与空间域磁场异常场磁位的三维泊松方程，并通过任意采样的二维傅里叶正变换将所述三维泊松方程转换为空间波数混合域的一维常微分方程，求解所述一维常微分方程得到波数域异常场磁位；
根据所述波数域异常场磁位与波数域异常场磁场强度的关系构建计算方程，并对方程求解得到波数域异常场磁场强度，对所述波数域异常场磁场强度做任意采样的二维傅里叶反变换，得到空间域异常场磁场强度；
通过所述空间域异常场磁场强度与空间域磁感应强度的关系，得到空间域磁感应强度。
在上述实施例中，通过对目标区域进行模型建立，得到目标区域的三维目标模型，并在三维目标模型中确定异常体，异常体可以是任意复杂形状，接着对包含有异常体的三维目标模型进行剖分，得到一系列节点，并根据磁化率分布数据对每个节点进行赋值，然后根据赋值结果计算节点的磁化强度，并通过任意采样的二维傅里叶正变换得到波数域异常场磁位，再通过二维傅里叶反变换得到空间域异常场磁场强度，最后通过空间域异常场磁场强度与空间域磁感应强度的关系得到空间域磁感应强度结果，并完成计算。
其中，空间域总场磁场强度为空间域背景场磁场强度与空间域异常场磁场强度之和，由于空间域异常场磁场强度数值非常小可以忽略不计，因此，此处空间域总场磁场强度等于空间域背景场磁场强度。
上述实施例中对三维目标模型的剖分方法除了非均匀剖分外，还可以对三维目标模型进行均匀剖分，均匀剖分方法包括：对三维目标模型的空间域x、y、z三个方向都均匀剖分，其中，x、y、z为三个分别垂直的方向；同时，上述实施例中在进行完任意采样的二维傅里叶反变换后，求解得到的空间域异常场磁场强度，不需要再进行紧算子迭代以及其他判断步骤，来对空间域异常场磁场强度进行判断，可以直接将求解得到的空间域异常场磁场强度应用在空间域异常场磁场强度与空间域磁感应强度的关系式上求解空间域磁感应强度。
实施例2，请参见图2-6：
基于目前三维频率域磁异常数值模拟方法基本还未实现任意采样，且存在边界效应和效率不高等问题，本发明提出一种任意采样的空间波数域三维磁场数值模拟方法。
本发明的研究方案为:
本发明提出的任意采样的空间波数域三维磁场数值模拟方法包括以下步骤：
步骤一：模型建立
对数值模拟的计算区域完成地质建模工作。首先确定整个计算区域的大小，之后再确定异常体的分布，异常体可以是任意复杂条件和任意复杂形状，异常体要在计算区域内部。一个简单模型的示意图如图2所示，图中异常体为球体。
步骤二：模型剖分
空间域建模：
建立模型后，对模型进行剖分，其x，y，z三个方向的采样点数分别为Nx,Ny,Nz。本发明的优势之一便是模型剖分在x,y,z三个方向都是任意的，可以在模型异常体变化较快的地方采取非均匀剖分，进行加密，在异常体变化较慢的地方或者无变化的地方进行稀疏采样。也可以三个方向都是均匀采样，如图4所示。对图2所示模型，为了对球体更好的拟合，在水平方向的剖分和采样可以如图3所示，同样的，在z方向也可以进行非均匀剖分。
依据空间域剖分，确定波数
Figure DEST_PATH_IMAGE098
的截止频率（
Figure 850882DEST_PATH_IMAGE098
的最大正值和最小负值）和
Figure 614570DEST_PATH_IMAGE098
的采样方式。
截止频率与空间域对应方向的最小剖分间隔有关，设x方向最小剖分间隔为
Figure DEST_PATH_IMAGE099
，y方向最小剖分间隔为
Figure DEST_PATH_IMAGE100
,则对应截止频率为：
Figure DEST_PATH_IMAGE101
；
Figure DEST_PATH_IMAGE102
；
在截止频率内采样能够保证所有的频谱信息都被采样。确定截止频率后，再确定采样个数，假设
Figure 716650DEST_PATH_IMAGE098
的采样个数分别为Nkx，Nky。
可选择均匀采样，即
Figure 131450DEST_PATH_IMAGE098
排列的间隔相同；也可选择在对数域均匀采样，对数域采样对于磁法数值模拟的波数选取较为适合。
对数区间采样时，设波数选取范围为
Figure DEST_PATH_IMAGE103
，波数域采样点数为
Figure DEST_PATH_IMAGE104
，在对数域采样进行等间隔采样，采样间隔为：
Figure DEST_PATH_IMAGE105
；
其中，
Figure DEST_PATH_IMAGE106
为一个小数，一般为
Figure DEST_PATH_IMAGE107
。
波数排列在
Figure DEST_PATH_IMAGE108
上为：
Figure DEST_PATH_IMAGE109
；
波数排列在
Figure DEST_PATH_IMAGE110
上为：
Figure DEST_PATH_IMAGE111
；
Figure 146371DEST_PATH_IMAGE098
的对数域采样都可用公式（1）和公式（2）的采样方式，由此给定了空间域x, y, z和波数域k x , k y 的排列。
步骤三：磁化率赋值
给图3或者图4中的节点进行磁化率赋值。异常体部分根据异常体的磁化率值赋给对应的每个节点，无异常部分的节点上磁化率为0，磁化率用
Figure DEST_PATH_IMAGE112
表示，为标量，单位为SI。
步骤四：计算节点对应的磁化强度M
根据地球主磁场模型IGRF，计算每个节点处的地球主磁场强度H0，为数值模拟中的背景场，即无异常时的磁场，单位为A/m，其三个方向的分量分别表示为H 0x 、H 0y 、H 0z 。节点处由异常体产生的磁场强度为H a ，为数值模拟中的异常场，即由磁化率异常所产生的磁场，单位为A/m，其三个分量分别为H ax 、H ay 、H az 。总场H为背景场与异常场之和。本发明只考虑弱磁条件下，即
Figure 725251DEST_PATH_IMAGE112
＜0.01SI的情况，H a 一般远小于H0，因此忽略不计。
背景场三个分量由下式计算，其中
Figure DEST_PATH_IMAGE113
表示背景场H0的L2范数，α为研究区域磁倾角，β为研究区域磁偏角。
Figure DEST_PATH_IMAGE114
；
Figure DEST_PATH_IMAGE115
；
Figure DEST_PATH_IMAGE116
；
因此在得到了每个节点的H0后，磁化强度M由下式计算：
Figure DEST_PATH_IMAGE117
；
步骤五：通过任意采样的二维傅里叶变换得到波数域异常场磁位
Figure DEST_PATH_IMAGE118
满足的一维常微分方程。
空间域磁场异常场磁位U a 和磁化强度M满足的方程如下所示：
Figure DEST_PATH_IMAGE119
；
对上式进行二维傅里叶变换。
该处任意采样的二维傅里叶正变换原理如下：
二维傅里叶正变换公式为：
Figure DEST_PATH_IMAGE120
；
式中
Figure DEST_PATH_IMAGE121
表示波数，
Figure DEST_PATH_IMAGE122
为空间域函数，
Figure DEST_PATH_IMAGE123
表示波数谱。
该二维变换是通过两次一维
Figure DEST_PATH_IMAGE124
傅里叶正变换完成的，首先对一维傅里叶正变换原理进行介绍。
一维傅里叶正变换可分别表示为：
Figure DEST_PATH_IMAGE125
；
其中表示波数，
Figure DEST_PATH_IMAGE126
为空间域函数，
Figure DEST_PATH_IMAGE127
为波数谱。
对上式中的正变换积分进行离散可得：
Figure DEST_PATH_IMAGE128
；
其中，N表示单元个数，
Figure DEST_PATH_IMAGE129
表示第j个单元，其中，i为虚数。
用二次形函数对
Figure DEST_PATH_IMAGE130
进行插值。单元内采用二次插值形函数拟合时，设任一单元内部三个节点坐标分别为
Figure 352060DEST_PATH_IMAGE037
，
Figure 903258DEST_PATH_IMAGE038
为中点，满足
Figure 79025DEST_PATH_IMAGE039
，单元内节点如图5所示。
每个节点上值分别为
Figure DEST_PATH_IMAGE131
，
Figure 20567DEST_PATH_IMAGE041
采用二次形函数进行表示可得：
Figure DEST_PATH_IMAGE132
；
其中，
Figure DEST_PATH_IMAGE133
；
上式可写为：
Figure DEST_PATH_IMAGE134
Figure DEST_PATH_IMAGE135
；
令
Figure DEST_PATH_IMAGE136
，
Figure DEST_PATH_IMAGE137
，
Figure DEST_PATH_IMAGE138
为单元内傅里叶变换节点系数，则上式简写为：
Figure DEST_PATH_IMAGE139
；
当波数
Figure 339815DEST_PATH_IMAGE050
不为0时，将N1、N2、N3代入
Figure DEST_PATH_IMAGE140
中，可得单元内傅里叶变换节点系数
Figure DEST_PATH_IMAGE141
Figure DEST_PATH_IMAGE142
Figure DEST_PATH_IMAGE143
Figure DEST_PATH_IMAGE144
；
Figure DEST_PATH_IMAGE145
Figure DEST_PATH_IMAGE146
Figure DEST_PATH_IMAGE147
积分核函数都包括
Figure DEST_PATH_IMAGE148
，其在
Figure DEST_PATH_IMAGE149
上单元积分解析解为：
Figure DEST_PATH_IMAGE150
Figure DEST_PATH_IMAGE151
Figure DEST_PATH_IMAGE152
；
Figure DEST_PATH_IMAGE153
因此，可得
Figure 810766DEST_PATH_IMAGE050
不为0时的
Figure 157434DEST_PATH_IMAGE147
半解析解为：
Figure DEST_PATH_IMAGE154
Figure DEST_PATH_IMAGE155
Figure DEST_PATH_IMAGE156
；
Figure DEST_PATH_IMAGE157
Figure DEST_PATH_IMAGE158
当波数
Figure 196060DEST_PATH_IMAGE050
为0时，
Figure DEST_PATH_IMAGE159
，
Figure DEST_PATH_IMAGE160
，
Figure DEST_PATH_IMAGE161
，进行简单积分即可求得零波数下傅里叶变换节点系数为：
Figure DEST_PATH_IMAGE162
；
将不同单元解析表达式进行累加，可得最终一维傅里叶正变换结果。易知，当空间域和频率域剖分不变时，傅里叶变换节点系数
Figure DEST_PATH_IMAGE163
和
Figure DEST_PATH_IMAGE164
均不变，提前将傅里叶变换系数算出并存储，可减少重复计算，提高算法效率，这也是本算法的优势之一。
二维傅里叶正变换是在对x做完一维傅里叶变换：
Figure DEST_PATH_IMAGE165
；
之后对
Figure DEST_PATH_IMAGE166
进行
Figure DEST_PATH_IMAGE167
方向一维傅里叶变换：
Figure DEST_PATH_IMAGE168
；
两次一维傅里叶变换的原理与过程完全相同，因此不再赘述。
通过任意采样的傅里叶变换，得到空间波数混合域一维常微分方程，保留z方向为空间域：
Figure 390501DEST_PATH_IMAGE076
（3）；
上式为波数域异常场磁位
Figure 634401DEST_PATH_IMAGE077
在所满足的一维常微分方程，其中
Figure 637123DEST_PATH_IMAGE077
表示波数域异常场磁位，
Figure 599263DEST_PATH_IMAGE078
、
Figure 244002DEST_PATH_IMAGE079
、
Figure 607987DEST_PATH_IMAGE080
表示波数域磁化强度
Figure DEST_PATH_IMAGE169
的x分量、y分量和z分量，
Figure 578348DEST_PATH_IMAGE021
、
Figure 778517DEST_PATH_IMAGE022
分别表示x、y方向的波数。保留了垂直方向在空间域，因此垂向也可以任意剖分。
与现有傅里叶变换所不同的是，现有的快速傅里叶变换在计算的过程中，边界效应影响较大必须扩边进行计算，扩边计算的加入增加了计算步骤，导致计算效率降低，而高斯函数的傅里叶变换通过高斯积分计算傅里叶变换积分公式，这种方式虽然削弱了截断效应的影响，但同时也牺牲了计算效率，同时标准的快速傅里叶变换和高斯函数的傅里叶变换只能均匀采样，不适用于需要非均匀采样的情况；非均匀快速傅里叶变换算法虽然实现了非均匀采样，但其内核仍然是传统傅里叶变换，仍存在边界效应。
因此提出上述基于形函数二次插值的任意采样的傅里叶变换方法，该方法基本不存在边界效应，且能够非均匀采样，并且具有较高效率，在偏微分方程求解中具有很强的适应性；此处任意采样的傅里叶变换是指，该方法的傅里叶变换采样点能够进行均匀采样，也能够进行非均匀采样。
采用上述基于形函数二次插值的任意采样的傅里叶变换方法与现有傅里叶变换相比，能够对傅里叶变换积分进行完整积分计算，正变换空间的采样相比于现有傅里叶变换更加准确，反变换空间的波数选取相比于现有傅里叶变换更加准确，进而在计算过程中不会存在频谱泄露，因此，与现有傅里叶变换相比不会存在边界效应，同时，上述基于形函数二次插值的任意采样的傅里叶变换方法在进行积分区间计算时，只需要准确计算傅里叶变换积分即可，与传统的傅里叶变换采样点需要均匀采样不同，上述方法采样点的选择对计算的结果不会产生影响，因此，上述方法在采用过程中能够实现非均匀采样。
步骤六：运用一维形函数法对波数域异常场磁位
Figure 210635DEST_PATH_IMAGE077
进行求解。
为了得到控制方程（3）的定解，需给出合适的边界条件，边界条件示意图如图6所示，在笛卡尔坐标系下，取Z轴垂直向下为正向，计算区域取水平地面为上边界Z min ，取地下离异常体足够远处为下边界Z max 。其上下边界条件满足：
上边界：
Figure DEST_PATH_IMAGE170
下边界：
Figure DEST_PATH_IMAGE172
；
其中，
Figure DEST_PATH_IMAGE173
。
联立得空间波数混合域异常场磁位满足的边值问题：
Figure 789963DEST_PATH_IMAGE086
；
运用变分法，得出与边值问题等价的变分问题：
Figure DEST_PATH_IMAGE174
在图6所示的笛卡尔坐标系中，沿z方向进行单元剖分，在每个单元内采用二次插值函数。对方程每一次求解的问题都是变分问题，对变分问题项进行单元分析、总体合成，得到一个由全体节点组成的五对角方程，采用追赶法可实现快速高效的求解，得到每个节点处的波数域异常场磁位
Figure 8586DEST_PATH_IMAGE077
。
步骤七：基于波数域异常场磁位
Figure 260707DEST_PATH_IMAGE077
求得波数域异常场磁场强度。
波数域异常场磁场强度
Figure DEST_PATH_IMAGE175
与波数域异常场磁位
Figure 146754DEST_PATH_IMAGE077
满足以下关系：
Figure DEST_PATH_IMAGE176
；
其中，i为虚数。
步骤八：运用任意采样的二维傅里叶反变换方法求得空间域异常场强度
Figure DEST_PATH_IMAGE177
。
此处运用任意采样二维傅里叶反变换也是本发明专利的一大创新，保证了本发明专利磁场数值模拟时反变换回空间域时可以任意采样，从而提高精度和效率。
二维任意采样傅里叶反变换公式为：
Figure DEST_PATH_IMAGE178
；
式中
Figure 599864DEST_PATH_IMAGE121
表示波数，
Figure 23892DEST_PATH_IMAGE122
为空间域函数，
Figure 548545DEST_PATH_IMAGE123
表示波数谱。该反变换公式同正变换公式形式完全相同，原理也相同，不再赘述。
步骤九：求解空间域磁感应强度B a ，结束数值模拟。
由异常场磁感应强度B a （单位为T）与异常场磁场强度H a 的关系可求得磁感应强度B a ，进而得到B a 的三个分量B ax ，B ay ，B az 。
Figure DEST_PATH_IMAGE179
；
其中，μ为介质的绝对磁导率，单位为H/m。μ与
Figure DEST_PATH_IMAGE180
的关系满足下式：
Figure DEST_PATH_IMAGE181
；
其中
Figure DEST_PATH_IMAGE182
为真空中磁导率，
Figure DEST_PATH_IMAGE183
。
实施例3，请参见图7：
对任意采样的空间波数域三维磁场数值模拟方法的精度和效率进行检验；
设计球体模型，背景磁场强度50000nT，磁倾角45º，磁偏角5.9º。模型计算大小为500m×500m×500m，范围：x方向-250 m~250 m，y方向-250 m~250 m，z方向0~500 m。异常体球体模型中心位于（0m，0m，250m），球半径为100m，球体磁化率为0.01SI，模型示意图如图1所示。将水平方向进行非均匀剖分，水平方向剖分方式如图2所示，最小间隔为1m，最大间隔为32m，在异常体以外采样间隔逐渐从1m增大至32m。z方向采用等间隔剖分方式。且三个方向节点个数均为101。波数域采样范围为-0.1~0.1，且采取对数采样方式，最小数为10-4，在对数域进行等间隔采样。波数域采样k x 和k y 都为101个。地面场值
Figure DEST_PATH_IMAGE184
的数值解与弱磁球体解析解对比，相对均方根误差分别为0.05%、0.05%、0.06%，如图7所示，占用内存0.8GB，耗时0.62s。
其中图7中（a）为B ax 数值解，（b）为B ax 解析解，(c)为B ax 数值解与解析解的绝对误差；（d）为B ay 数值解，（e）为B ay 解析解，(f)为B ay 数值解与解析解的绝对误差；（g）为B az 数值解，（h）为B az 解析解，(i)为B az 数值解与解析解的绝对误差。
本申请实施例提供一种任意采样的空间波数域三维磁场数值模拟系统，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述任意采样的空间波数域三维磁场数值模拟方法中任一所述方法的步骤。
上述的任意采样的空间波数域三维磁场数值模拟系统，可以实现上述的任意采样的空间波数域三维磁场数值模拟方法的各个实施例，且能达到相同的有益效果，此处，不做赘述。
以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。

Claims (9)

1.一种任意采样的空间波数域三维磁场数值模拟方法，其特征在于，包括：
构建含有异常体的目标区域三维目标模型，对所述三维目标模型进行非均匀剖分得到一系列节点，并根据磁化率分布数据对每个节点进行磁化率赋值，得到每个节点的磁化率，其中，包含异常体的节点根据异常体的磁化率分布数据进行赋值，不包含异常体的节点赋值为0，其中，非均匀剖分包括以下任意方式：
方式一：对预设的第一区域采取非均匀剖分，进行加密；其中，所述第一区域满足以下公式：
Figure 626682DEST_PATH_IMAGE001
；
Figure 463182DEST_PATH_IMAGE002
为所述第一区域对应节点的剩余密度，
Figure 284507DEST_PATH_IMAGE003
为第一区域对应节点周围的第j个节点的剩余密度，
Figure 695897DEST_PATH_IMAGE004
为第一区域周围节点数量；
Figure 851941DEST_PATH_IMAGE005
为权重，取值范围为（0,1）；
方式二：对预设的第二区域进行稀疏采样，其中，所述第二区域满足以下公式：
Figure 6979DEST_PATH_IMAGE006
；
Figure 682811DEST_PATH_IMAGE007
为所述第二区域对应节点的剩余密度，
Figure 999523DEST_PATH_IMAGE008
为第二区域对应节点周围的第j个节点的剩余密度，
Figure 144327DEST_PATH_IMAGE004
为第一区域周围节点数量；
Figure 837477DEST_PATH_IMAGE005
为权重，取值范围为（0,1）；
根据地球主磁场模型计算每个节点处的地球主磁场强度，并将所述地球主磁场强度作为空间域背景场磁场强度计算空间域总场磁场强度，根据空间域总场磁场强度和磁化率之间的关系得到磁化强度的计算模型，通过磁化强度的计算模型得到磁化强度；
基于所述磁化强度构建磁化强度与空间域磁场异常场磁位的三维泊松方程，并通过任意采样的二维傅里叶正变换将所述三维泊松方程转换为空间波数混合域的一维常微分方程，求解所述一维常微分方程得到波数域异常场磁位；
根据所述波数域异常场磁位与波数域异常场磁场强度的关系构建计算方程，并对方程求解得到波数域异常场磁场强度，对所述波数域异常场磁场强度做任意采样的二维傅里叶反变换，得到空间域异常场磁场强度；
通过所述空间域异常场磁场强度与空间域磁感应强度的关系，得到空间域磁感应强度。
2.根据权利要求1所述的任意采样的空间波数域三维磁场数值模拟方法，其特征在于，所述磁化强度的计算模型为：
Figure 898974DEST_PATH_IMAGE009
；
其中，M表示磁化强度，
Figure 386587DEST_PATH_IMAGE010
表示磁化率，H表示空间域总场磁场强度，空间域总场磁场强度为空间域背景场磁场强度和空间域异常场磁场强度之和。
3.根据权利要求1所述的任意采样的空间波数域三维磁场数值模拟方法，其特征在于，
所述三维泊松方程为：
Figure 251643DEST_PATH_IMAGE011
；
其中，
Figure 482905DEST_PATH_IMAGE012
表示空间域磁场异常场磁位，M表示磁化强度，式中，
Figure 664487DEST_PATH_IMAGE013
，
Figure 588581DEST_PATH_IMAGE014
，i，j，k分别为x，y，z方向的单位向量；
上式展开为
Figure 176819DEST_PATH_IMAGE015
；
其中，
Figure 211771DEST_PATH_IMAGE016
、
Figure 716702DEST_PATH_IMAGE017
、
Figure 342856DEST_PATH_IMAGE018
分别为磁化强度M在x，y，z方向的分量，
Figure 182504DEST_PATH_IMAGE019
为偏导数符号。
4.根据权利要求1所述的任意采样的空间波数域三维磁场数值模拟方法，其特征在于，所述任意采样的二维傅里叶正变换公式为：
Figure 755568DEST_PATH_IMAGE020
；
其中，
Figure 380585DEST_PATH_IMAGE021
表示x方向的波数，
Figure 177639DEST_PATH_IMAGE022
表示y方向上的波数，
Figure 740470DEST_PATH_IMAGE023
表示空间域函数，
Figure 851645DEST_PATH_IMAGE024
表示波数谱；
所述任意采样的二维傅里叶正变换公式通过两次一维傅里叶正变换得到，两次一维傅里叶正变换分别为；
对
Figure 65589DEST_PATH_IMAGE023
进行x方向一维傅里叶正变换，变换式为：
Figure 751654DEST_PATH_IMAGE025
；
其中，x，y表示两个相互垂直的方向；
Figure 316628DEST_PATH_IMAGE021
表示x方向的波数，
Figure 497073DEST_PATH_IMAGE026
表示空间域函数，
Figure 96682DEST_PATH_IMAGE027
为对
Figure 704381DEST_PATH_IMAGE026
做x方向做一维傅里叶变换后的波数谱；
对
Figure 241804DEST_PATH_IMAGE027
进行y方向一维傅里叶正变换，变换式为：
Figure 225940DEST_PATH_IMAGE028
；
其中，
Figure 945635DEST_PATH_IMAGE022
表示y方向上的波数，
Figure 458655DEST_PATH_IMAGE029
为对
Figure 247489DEST_PATH_IMAGE026
做二维傅里叶变换后的波数谱。
5.根据权利要求4所述的任意采样的空间波数域三维磁场数值模拟方法，其特征在于，所述一维傅里叶正变换具体变换方法为：
设连续一维傅里叶正变换分别表示为：
Figure 504158DEST_PATH_IMAGE030
；
其中，
Figure 343938DEST_PATH_IMAGE021
表示波数，
Figure 293439DEST_PATH_IMAGE031
表示空间域函数，
Figure 805454DEST_PATH_IMAGE032
表示波数谱；
对上述连续一维傅里叶正变换进行离散得到：
Figure 600235DEST_PATH_IMAGE033
；
其中，
Figure 763363DEST_PATH_IMAGE034
表示单元个数，
Figure 680503DEST_PATH_IMAGE035
表示第j个单元，i为虚数，
Figure 912771DEST_PATH_IMAGE021
表示x方向上的波数；
用二次形函数对
Figure 776821DEST_PATH_IMAGE036
进行插值：
单元内采用二次插值形函数拟合时，设任一单元内部三个节点坐标分别为
Figure 794456DEST_PATH_IMAGE037
，
Figure 836492DEST_PATH_IMAGE038
为中点，满足
Figure 306788DEST_PATH_IMAGE039
，每个节点上值分别为
Figure 505688DEST_PATH_IMAGE040
，
Figure 377829DEST_PATH_IMAGE041
采用二次形函数进行表示得：
Figure 89302DEST_PATH_IMAGE042
；
其中，N1、N2、N3表示二次插值函数，分别为：
Figure 46894DEST_PATH_IMAGE043
；
上式写为：
Figure 252747DEST_PATH_IMAGE044
Figure 776133DEST_PATH_IMAGE045
；
令
Figure 148253DEST_PATH_IMAGE046
，
Figure 858720DEST_PATH_IMAGE047
，
Figure 337106DEST_PATH_IMAGE048
为单元内傅里叶变换节点系数，即W1、W2、W3分别表示每个节点对应的傅里叶变换系数，则上式写为：
Figure 449418DEST_PATH_IMAGE049
；
当波数
Figure 784585DEST_PATH_IMAGE050
不为0时，将N1、N2、N3代入
Figure 966036DEST_PATH_IMAGE051
中，得单元内傅里叶变换节点系数W1、W2、W3分别为：
Figure 248113DEST_PATH_IMAGE052
Figure 214932DEST_PATH_IMAGE053
Figure 940574DEST_PATH_IMAGE054
Figure 891212DEST_PATH_IMAGE055
；
Figure 242559DEST_PATH_IMAGE056
Figure 798305DEST_PATH_IMAGE057
Figure 458963DEST_PATH_IMAGE058
积分核函数都包括
Figure 100159DEST_PATH_IMAGE059
，其在
Figure 520777DEST_PATH_IMAGE060
上单元积分解析解为：
Figure 462188DEST_PATH_IMAGE061
Figure 529632DEST_PATH_IMAGE062
Figure 658125DEST_PATH_IMAGE063
；
Figure 351274DEST_PATH_IMAGE064
因此，得
Figure 412771DEST_PATH_IMAGE065
不为0时的
Figure 149652DEST_PATH_IMAGE066
半解析解为：
Figure 31020DEST_PATH_IMAGE067
Figure 262282DEST_PATH_IMAGE068
Figure 443864DEST_PATH_IMAGE069
；
Figure 853111DEST_PATH_IMAGE070
Figure 956196DEST_PATH_IMAGE071
当波数
Figure 725569DEST_PATH_IMAGE072
为0时，
Figure 745347DEST_PATH_IMAGE073
，
Figure 105921DEST_PATH_IMAGE074
，
Figure 696302DEST_PATH_IMAGE075
，
Figure 534945DEST_PATH_IMAGE076
分别表示波数为0时的傅里叶变换系数，进行简单积分即求得零波数下傅里叶变换节点系数为：
Figure 379536DEST_PATH_IMAGE077
；
将不同单元解析表达式进行累加，得最终一维傅里叶正变换结果。
6.根据权利要求1所述的任意采样的空间波数域三维磁场数值模拟方法，其特征在于，所述一维常微分方程通过将所述三维泊松方程进行水平方向二维傅里叶正变换得到，所述一维常微分方程为：
Figure 645432DEST_PATH_IMAGE078
；
其中，
Figure 254268DEST_PATH_IMAGE079
表示波数域异常场磁位，
Figure 631022DEST_PATH_IMAGE080
、
Figure 110545DEST_PATH_IMAGE081
、
Figure 531031DEST_PATH_IMAGE082
表示波数域磁化强度
Figure 830426DEST_PATH_IMAGE083
的x分量、y分量和z分量，
Figure 10871DEST_PATH_IMAGE021
、
Figure 95633DEST_PATH_IMAGE022
分别表示x、y方向的波数，
Figure 234490DEST_PATH_IMAGE084
为偏导数符号；
所述求解所述一维常微分方程得到波数域异常场磁位，包括：
在笛卡尔坐标系下，取Z轴垂直向下为正向，计算区域取水平地面为上边界Z min ，取地下离异常体足够远处为下边界Z max ，其上下边界条件满足：
上边界：
Figure 286760DEST_PATH_IMAGE085
；
下边界：
Figure 739738DEST_PATH_IMAGE086
；
其中，
Figure 177541DEST_PATH_IMAGE087
；
联立一维常微分方程和上下边界，得到：
Figure 487300DEST_PATH_IMAGE088
；
运用变分法，得到与边值问题等价的变分问题：
Figure 26866DEST_PATH_IMAGE089
沿z方向进行单元剖分，在每个单元内采用二次插值函数，得到每个节点处的波数域异常场磁位
Figure 283535DEST_PATH_IMAGE079
。
7.根据权利要求1所述的任意采样的空间波数域三维磁场数值模拟方法，其特征在于，所述根据所述波数域异常场磁位与波数域异常场磁场强度的关系构建计算方程，包括：
波数域异常场磁位
Figure 342889DEST_PATH_IMAGE079
与波数域异常场磁场强度
Figure 26811DEST_PATH_IMAGE090
满足如下关系式：
Figure 319252DEST_PATH_IMAGE091
；
其中，i为虚数。
8.根据权利要求1所述的任意采样的空间波数域三维磁场数值模拟方法，其特征在于，所述二维傅里叶反变换公式为：
Figure 645191DEST_PATH_IMAGE092
；
其中，
Figure 792008DEST_PATH_IMAGE021
、
Figure 912410DEST_PATH_IMAGE022
表示波数，
Figure 957727DEST_PATH_IMAGE023
为空间域函数，
Figure 556198DEST_PATH_IMAGE024
表示波数谱；
所述空间域异常场磁场强度与空间域磁感应强度的关系为：
Figure 58986DEST_PATH_IMAGE093
；
其中，
Figure 615869DEST_PATH_IMAGE094
表示介质的绝对磁导率，单位为H/m，
Figure 86165DEST_PATH_IMAGE094
与
Figure 285065DEST_PATH_IMAGE095
的关系满足下式：
Figure 157206DEST_PATH_IMAGE096
；
其中，
Figure 868679DEST_PATH_IMAGE097
表示真空中磁导率，
Figure 91850DEST_PATH_IMAGE098
H/m。
9.一种任意采样的空间波数域三维磁场数值模拟系统，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现上述权利要求1至8任一所述方法的步骤。
CN202211225001.7A 2022-10-09 2022-10-09 一种任意采样的空间波数域三维磁场数值模拟方法及系统 Active CN115292973B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211225001.7A CN115292973B (zh) 2022-10-09 2022-10-09 一种任意采样的空间波数域三维磁场数值模拟方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211225001.7A CN115292973B (zh) 2022-10-09 2022-10-09 一种任意采样的空间波数域三维磁场数值模拟方法及系统

Publications (2)

Publication Number Publication Date
CN115292973A true CN115292973A (zh) 2022-11-04
CN115292973B CN115292973B (zh) 2023-01-24

Family

ID=83834769

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211225001.7A Active CN115292973B (zh) 2022-10-09 2022-10-09 一种任意采样的空间波数域三维磁场数值模拟方法及系统

Country Status (1)

Country Link
CN (1) CN115292973B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116029110A (zh) * 2022-12-26 2023-04-28 长江岩土工程有限公司 一种考虑退磁和剩磁的磁场数值模拟方法和系统
CN116774301A (zh) * 2023-06-26 2023-09-19 中国人民解放军火箭军工程大学 重磁位场正则化向下延拓方法、系统、电子设备及介质

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102440778A (zh) * 2010-09-30 2012-05-09 株式会社东芝 磁共振成像装置
US20130275102A1 (en) * 2010-07-19 2013-10-17 Terje Graham Vold Computer simulation of electromagnetic fields
CN112287534A (zh) * 2020-10-21 2021-01-29 中南大学 基于nufft的二维磁异常快速正演模拟方法和装置
CN113656750A (zh) * 2021-10-20 2021-11-16 中南大学 基于空间波数域的强磁介质的磁感应强度计算方法
CN113962077A (zh) * 2021-10-20 2022-01-21 中南大学 三维各向异性强磁场数值模拟方法、装置、设备及介质
CN114004127A (zh) * 2021-11-05 2022-02-01 中南大学 二维主轴各向异性强磁场数值模拟方法、装置、设备及介质
CN114065586A (zh) * 2021-11-22 2022-02-18 中南大学 一种三维大地电磁空间-波数域有限元数值模拟方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130275102A1 (en) * 2010-07-19 2013-10-17 Terje Graham Vold Computer simulation of electromagnetic fields
CN102440778A (zh) * 2010-09-30 2012-05-09 株式会社东芝 磁共振成像装置
CN112287534A (zh) * 2020-10-21 2021-01-29 中南大学 基于nufft的二维磁异常快速正演模拟方法和装置
CN113656750A (zh) * 2021-10-20 2021-11-16 中南大学 基于空间波数域的强磁介质的磁感应强度计算方法
CN113962077A (zh) * 2021-10-20 2022-01-21 中南大学 三维各向异性强磁场数值模拟方法、装置、设备及介质
CN114004127A (zh) * 2021-11-05 2022-02-01 中南大学 二维主轴各向异性强磁场数值模拟方法、装置、设备及介质
CN114065586A (zh) * 2021-11-22 2022-02-18 中南大学 一种三维大地电磁空间-波数域有限元数值模拟方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
SHIKUN DAI ET AL.: "Three-dimensional magnetotelluric modeling in a mixed space-wavenumber domain", 《GEOPHYSICS》 *
戴世坤等: "空间-波数域三维大地电磁场积分方程法数值模拟", 《地球物理学报》 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116029110A (zh) * 2022-12-26 2023-04-28 长江岩土工程有限公司 一种考虑退磁和剩磁的磁场数值模拟方法和系统
CN116774301A (zh) * 2023-06-26 2023-09-19 中国人民解放军火箭军工程大学 重磁位场正则化向下延拓方法、系统、电子设备及介质
CN116774301B (zh) * 2023-06-26 2024-05-14 中国人民解放军火箭军工程大学 重磁位场正则化向下延拓方法、系统、电子设备及介质

Also Published As

Publication number Publication date
CN115292973B (zh) 2023-01-24

Similar Documents

Publication Publication Date Title
CN115292973B (zh) 一种任意采样的空间波数域三维磁场数值模拟方法及系统
CN112287534B (zh) 基于nufft的二维磁异常快速正演模拟方法和装置
CN106646645B (zh) 一种重力正演加速方法
CN108197389A (zh) 二维强磁性体磁场的快速、高精度数值模拟方法
CN113962077B (zh) 三维各向异性强磁场数值模拟方法、装置、设备及介质
CN109254327B (zh) 三维强磁性体的勘探方法及勘探系统
Abbas et al. An exact and fast computation of discrete Fourier transform for polar and spherical grid
Wu Comparison of 3-D Fourier forward algorithms for gravity modelling of prismatic bodies with polynomial density distribution
CN114021408A (zh) 二维强磁场数值模拟方法、装置、设备及介质
Piip 2D inversion of refraction traveltime curves using homogeneous functions
Ying et al. The phase flow method
CN114004127A (zh) 二维主轴各向异性强磁场数值模拟方法、装置、设备及介质
Brandt et al. Renormalization multigrid (RMG): Statistically optimal renormalization group flow and coarse-to-fine Monte Carlo acceleration
Stepanova et al. Analytical models of the physical fields of the Earth in regional version with ellipticity
CN115017782B (zh) 考虑介质各向异性的三维天然源电磁场计算方法
CN115659579A (zh) 基于3d as-ft的三维重力场数值模拟方法及系统
CN115113286B (zh) 基于多分量频率域航空电磁数据融合三维反演方法
CN115795231A (zh) 一种空间波数混合域三维强磁场迭代方法及系统
Bello et al. Forward plane‐wave electromagnetic model in three dimensions using hybrid finite volume–integral equation scheme
Lan et al. High-precision 3D reconstruction of multiple magnetic targets based on center weighting method
CN116244877B (zh) 基于3d傅里叶变换的三维磁场数值模拟方法及系统
Yin et al. A fast 3D gravity forward algorithm based on circular convolution
Schnizer et al. Magnetic field analysis for superferric accelerator magnets using elliptic multipoles and its advantages
CN116911146B (zh) 三维重力场全息数值模拟及cpu-gpu加速方法
Nengem et al. A Study of the Effects of the Shape Parameter and Type of Data Points Locations on the Accuracy of the Hermite-Based Symmetric Approach Using Positive Definite Radial Kernels

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant