CN113806686B - 大规模复杂地质体重力梯度快速计算方法、装置和设备 - Google Patents

大规模复杂地质体重力梯度快速计算方法、装置和设备 Download PDF

Info

Publication number
CN113806686B
CN113806686B CN202111372974.9A CN202111372974A CN113806686B CN 113806686 B CN113806686 B CN 113806686B CN 202111372974 A CN202111372974 A CN 202111372974A CN 113806686 B CN113806686 B CN 113806686B
Authority
CN
China
Prior art keywords
gravity gradient
calculation
abnormal body
abnormal
kernel function
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.)
Active
Application number
CN202111372974.9A
Other languages
English (en)
Other versions
CN113806686A (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 CN202111372974.9A priority Critical patent/CN113806686B/zh
Publication of CN113806686A publication Critical patent/CN113806686A/zh
Application granted granted Critical
Publication of CN113806686B publication Critical patent/CN113806686B/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
    • 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
    • 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
    • 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

Landscapes

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

Abstract

本申请涉及一种大规模复杂地质体重力梯度快速计算方法、装置和设备。所述方法包括：对计算区域和异常体分别进行网格剖分，异常体网格间隔与计算区域一致，计算任一观测高度的重力梯度核函数积分系数，并进一步得到任一观测高度上任一层异常体对应的核函数系数矩阵，与异常体密度值扩展矩阵进行二维快速离散卷积，得到任一观测高度上任一层异常体对应的重力梯度空间域卷积结果；将各层卷积结果累加，得到计算区域中任一观测高度上重力梯度各个分量的值。本发明解决了重力梯度张量计算过程中二维离散卷积计算占用内存大、计算速度低的问题，为大规模复杂地质体的模拟提供了技术支撑。

Description

大规模复杂地质体重力梯度快速计算方法、装置和设备
技术领域
本申请涉及计算机技术领域，特别是涉及一种大规模复杂地质体重力梯度快速计算方法、装置和计算机设备。
背景技术
随着航空重力仪器的发展，单一的重力垂直分量的测量已经不能满足精细化勘探的需求，因此航空重力梯度仪应运而生，航空重力梯度测量技术由于其分辨率高且不受地形条件的影响等优势，在矿产资源勘查、地质科学研究及国防建设方面逐渐受到广泛应用。
重力梯度张量正演计算是重力梯度张量反演和重力梯度张量分量转换的基础，根据模型研究可以了解不同物性及几何参数的异常体产生的重力梯度张量异常的空间特征，更能在缺乏实测数据的情况为其他研究提供对比数据及模型约束。
目前，针对重力梯度张量正演计算的研究相对较少，并且现有的计算方法存在占用内存大计算效率低等问题。
对于以上问题，文献（陈涛, 张贵宾, 索奎,等. 不同小波基函数在重力梯度异常正演计算中的应用研究. 地球物理学报物探与化探, 2015, 39: 91-97.）采用基于小波变换理论，给出了小波域中的重力梯度张量正演方程，并针对不同小波基函数对计算效率进行了对比分析，实现了重力梯度张量的计算，但随着网格剖分数量的增加其矩阵相乘时间及占用内存均会快速增加，计算效率低。文献(Zhang S , Meng X C , Zhang M , et al.The improved residual node density based gravity forward method and itsapplication[J]. Journal of Applied Geophysics, 2018, 159:765-772.)提出了一种新的空间域重力梯度张量正演方法，但是对于任意密度分布的复杂地质体模型，该方法计算效率低。
文献（Wu, L., Chen, ChenL. Fourier forward modeling of vector andtensor gravity fields due to prismatic bodies with variable density contrast.Geophysics: Journal of the Society of Exploration Geophysicists, 2016, 81(1):G13-G26.）提出了一种Gauss-FFT方法，有效压制常规快速傅里叶变换的边界效应，提高了重力场及梯度张量的精度，但是随着高斯点的增加，计算量和计算时间成倍增加，不能实现计算精度和效率的统一，无法实现大规模快速反演的要求。
现有的重力梯度张量计算方法存在计算效率低、占用内存大的问题。
发明内容
基于此，有必要针对上述技术问题，提供一种能够减少计算内存占用，提高计算效率的大规模复杂地质体重力梯度快速计算方法、装置和计算机设备。
一种大规模复杂地质体重力梯度快速计算方法，所述方法包括：
根据待探测的异常体的分布构建三维坐标系，并确定计算区域，将所述计算区域沿x，y方向进行等间隔剖分，沿z方向进行灵活剖分，得到计算区域模型；
对所述异常体进行网格剖分，得到异常体网格模型，并在所述异常体网格模型中设置异常体密度值，得到异常体密度分布源函数；其中，异常体网格模型与所述计算区域模型对应部分的网格间隔一致；
根据所述计算区域模型、所述异常体密度分布源函数，得到任一观测高度的重力梯度核函数积分系数；
根据所述重力梯度核函数积分系数，得到任一观测高度上任一层异常体对应的核函数系数矩阵；
根据所述核函数系数矩阵和预设的异常体密度值扩展矩阵，通过二维傅里叶变换实现快速离散卷积计算，并进一步进行傅里叶反变换，得到任一观测高度上任一层异常体对应的重力梯度空间域卷积结果；
将各层异常体对应的重力梯度空间域卷积结果进行累加，得到所述计算区域中任一观测高度上重力梯度各个分量的值。
在其中一个实施例中，还包括：根据待探测的异常体的分布构建三维坐标系，x轴方向指向正东，y轴方向指向正北，z轴垂直向下；
确定计算区域范围；
将所述计算区域沿x，y方向进行等间隔剖分，沿z方向根据场变化的快慢进行灵活剖分，在场变化快的区域减小间隔，在场变化慢的区域增大间隔，得到计算区域模型。
在其中一个实施例中，还包括：根据所述计算区域模型、所述异常体密度分布源函数，得到观测高度
Figure 613677DEST_PATH_IMAGE001
上的重力梯度核函数积分系数为：
Figure 301010DEST_PATH_IMAGE002
Figure 379956DEST_PATH_IMAGE003
Figure 519950DEST_PATH_IMAGE004
其中，
Figure 959022DEST_PATH_IMAGE005
分别表示
Figure 704124DEST_PATH_IMAGE006
的核函数积分系数；
Figure 203238DEST_PATH_IMAGE007
表示万有引力常数；
Figure 82726DEST_PATH_IMAGE008
表示编号为
Figure 591068DEST_PATH_IMAGE009
的计算区域模型单元的中心坐标，
Figure 987414DEST_PATH_IMAGE010
表示编号为
Figure 657430DEST_PATH_IMAGE011
的异常体网格模型单元的中心坐标，
Figure 975278DEST_PATH_IMAGE012
，
Figure 772464DEST_PATH_IMAGE013
,
Figure 23317DEST_PATH_IMAGE014
，
Figure 864234DEST_PATH_IMAGE015
,
Figure 731696DEST_PATH_IMAGE016
，
Figure 565529DEST_PATH_IMAGE017
和
Figure 936467DEST_PATH_IMAGE018
分别为所述计算区域模型在x，y方向上的网格数，
Figure 151548DEST_PATH_IMAGE019
、
Figure 240727DEST_PATH_IMAGE020
和
Figure 363403DEST_PATH_IMAGE021
分别为所述异常体网格模型在x，y，z方向上的网格数，
Figure 605160DEST_PATH_IMAGE022
，
Figure 53459DEST_PATH_IMAGE023
，
Figure 364354DEST_PATH_IMAGE024
，
Figure 290722DEST_PATH_IMAGE025
、
Figure 839515DEST_PATH_IMAGE026
、
Figure 439474DEST_PATH_IMAGE027
、
Figure 237666DEST_PATH_IMAGE028
、
Figure 967724DEST_PATH_IMAGE029
，
Figure 433341DEST_PATH_IMAGE030
、
Figure 708595DEST_PATH_IMAGE031
和
Figure 994083DEST_PATH_IMAGE032
分别为x，y，z方向上的网格间隔，
Figure 262254DEST_PATH_IMAGE033
。
在其中一个实施例中，还包括：根据所述重力梯度核函数积分系数，得到观测高度
Figure 520060DEST_PATH_IMAGE034
上第n层异常体对应的核函数系数矩阵为：
Figure 481062DEST_PATH_IMAGE035
Figure 768693DEST_PATH_IMAGE036
Figure 574975DEST_PATH_IMAGE037
。
在其中一个实施例中，还包括：获取第n层异常体的密度值
Figure 749604DEST_PATH_IMAGE038
，表示为矩阵
Figure 615929DEST_PATH_IMAGE039
；
将矩阵
Figure 79272DEST_PATH_IMAGE039
填零扩展，得到异常体密度值扩展矩阵：
Figure 439977DEST_PATH_IMAGE040
通过二维傅里叶变换实现快速离散卷积计算，并进一步进行傅里叶反变换，得到任一观测高度上任一层异常体对应的重力梯度空间域卷积结果。
在其中一个实施例中，还包括：通过二维傅里叶变换实现快速离散卷积计算，即所述核函数系数矩阵与所述异常体密度值扩展矩阵的快速相乘：
Figure 469113DEST_PATH_IMAGE041
Figure 37497DEST_PATH_IMAGE042
Figure 519294DEST_PATH_IMAGE043
其中，
Figure 136220DEST_PATH_IMAGE044
表示二维离散傅里叶正变换算子；
通过二维离散傅里叶反变换，得到观测高度
Figure 272060DEST_PATH_IMAGE045
上第n层异常体对应的重力梯度空间域卷积结果：
Figure 745767DEST_PATH_IMAGE046
Figure 980439DEST_PATH_IMAGE047
Figure 932214DEST_PATH_IMAGE048
其中，
Figure 686675DEST_PATH_IMAGE049
表示二维离散傅里叶反变换算子。
在其中一个实施例中，还包括：将各层异常体对应的重力梯度空间域卷积结果进行累加，得到所述计算区域中观测高度
Figure 534545DEST_PATH_IMAGE050
上重力梯度各个分量的值为：
Figure 990934DEST_PATH_IMAGE051
Figure 746401DEST_PATH_IMAGE052
Figure 339056DEST_PATH_IMAGE053
其中，
Figure 669412DEST_PATH_IMAGE054
为所述计算区域在z方向上剖分的网格数。
一种大规模复杂地质体重力梯度快速计算装置，所述装置包括：
应用场景设置模块，用于根据待探测的异常体的分布构建三维坐标系，并确定计算区域，将所述计算区域沿x，y方向进行等间隔剖分，沿z方向进行灵活剖分，得到计算区域模型；对所述异常体进行网格剖分，得到异常体网格模型，并在所述异常体网格模型中设置异常体密度值，得到异常体密度分布源函数；其中，异常体网格模型与所述计算区域模型对应部分的网格间隔一致；
重力梯度核函数积分系数确定模块，用于根据所述计算区域模型、所述异常体密度分布源函数，得到任一观测高度的重力梯度核函数积分系数；
核函数系数矩阵确定模块，用于根据所述重力梯度核函数积分系数，得到任一观测高度上任一层异常体对应的核函数系数矩阵；
二维快速离散卷积模块，用于根据所述核函数系数矩阵和预设的异常体密度值扩展矩阵，通过二维傅里叶变换实现快速离散卷积计算，并进一步进行傅里叶反变换，得到任一观测高度上任一层异常体对应的重力梯度空间域卷积结果；
重力梯度张量各分量输出模块，用于将各层异常体对应的重力梯度空间域卷积结果进行累加，得到所述计算区域中任一观测高度上重力梯度各个分量的值。
一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现以下步骤：
根据待探测的异常体的分布构建三维坐标系，并确定计算区域，将所述计算区域沿x，y方向进行等间隔剖分，沿z方向进行灵活剖分，得到计算区域模型；
对所述异常体进行网格剖分，得到异常体网格模型，并在所述异常体网格模型中设置异常体密度值，得到异常体密度分布源函数；其中，异常体网格模型与所述计算区域模型对应部分的网格间隔一致；
根据所述计算区域模型、所述异常体密度分布源函数，得到任一观测高度的重力梯度核函数积分系数；
根据所述重力梯度核函数积分系数，得到任一观测高度上任一层异常体对应的核函数系数矩阵；
根据所述核函数系数矩阵和预设的异常体密度值扩展矩阵，通过二维傅里叶变换实现快速离散卷积计算，并进一步进行傅里叶反变换，得到任一观测高度上任一层异常体对应的重力梯度空间域卷积结果；
将各层异常体对应的重力梯度空间域卷积结果进行累加，得到所述计算区域中任一观测高度上重力梯度各个分量的值。
一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现以下步骤：
根据待探测的异常体的分布构建三维坐标系，并确定计算区域，将所述计算区域沿x，y方向进行等间隔剖分，沿z方向进行灵活剖分，得到计算区域模型；
对所述异常体进行网格剖分，得到异常体网格模型，并在所述异常体网格模型中设置异常体密度值，得到异常体密度分布源函数；其中，异常体网格模型与所述计算区域模型对应部分的网格间隔一致；
根据所述计算区域模型、所述异常体密度分布源函数，得到任一观测高度的重力梯度核函数积分系数；
根据所述重力梯度核函数积分系数，得到任一观测高度上任一层异常体对应的核函数系数矩阵；
根据所述核函数系数矩阵和预设的异常体密度值扩展矩阵，通过二维傅里叶变换实现快速离散卷积计算，并进一步进行傅里叶反变换，得到任一观测高度上任一层异常体对应的重力梯度空间域卷积结果；
将各层异常体对应的重力梯度空间域卷积结果进行累加，得到所述计算区域中任一观测高度上重力梯度各个分量的值。
上述大规模复杂地质体重力梯度快速计算方法、装置、计算机设备和存储介质，构建三维坐标系确定计算区域后，将计算区域沿x，y方向进行等间隔剖分，沿z方向进行灵活剖分，得到计算区域模型；对异常体进行网格剖分，其与计算区域模型的网格间隔一致，并得到异常体密度分布源函数；计算任一观测高度的重力梯度核函数积分系数，并进一步处理得到任一观测高度上任一层异常体对应的核函数系数矩阵；根据核函数系数矩阵和预设的异常体密度值扩展矩阵，通过二维傅里叶变换实现快速离散卷积计算，并进一步进行傅里叶反变换，得到任一观测高度上任一层异常体对应的重力梯度空间域卷积结果；将各层异常体对应的重力梯度空间域卷积结果进行累加，得到计算区域中任一观测高度上重力梯度各个分量的值。本发明通过对计算区域和异常体分别进行网格剖分，避免了传统二维离散卷积算法核函数系数矩阵中计算区域和异常体剖分网格数必须一致，并且需要用零元素填充以形成核函数系数矩阵的问题，不仅可以处理更一般的卷积情况，并且可以减少核函数系数矩阵元素数存储与计算。
附图说明
图1为一个实施例中大规模复杂地质体重力梯度快速计算方法的流程示意图；
图2为另一个实施例中大规模复杂地质体重力梯度快速计算方法的流程示意图；
图3为一个具体实施例中模型示意图；
图4为一个具体实施例中本发明方法计算的地表重力梯度分量gxx平面等值线图；
图5为一个具体实施例中解析解计算的地表重力梯度分量gxx平面等值线图；
图6为一个具体实施例中解析解与本发明方法计算的相对误差结果图；
图7为一个具体实施例中本发明方法计算的地表重力梯度分量gyy平面等值线图；
图8为一个具体实施例中解析解计算的地表重力梯度分量gyy平面等值线图；
图9为一个具体实施例中解析解与本发明方法计算的相对误差结果图；
图10为一个具体实施例中本发明方法计算的地表重力梯度分量gzz平面等值线图；
图11为一个具体实施例中解析解计算的地表重力梯度分量gzz平面等值线图；
图12为一个具体实施例中解析解与本发明方法计算的相对误差结果图；
图13为一个具体实施例中本发明方法和传统二维离散卷积算法的计算时间随网格剖分数目变化的曲线图；
图14为一个实施例中大规模复杂地质体重力梯度快速计算装置的结构框图；
图15为一个实施例中计算机设备的内部结构图。
具体实施方式
为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。
在一个实施例中，如图1所示，提供了一种大规模复杂地质体重力梯度快速计算方法，包括以下步骤：
步骤102，根据待探测的异常体的分布构建三维坐标系，并确定计算区域，将计算区域沿x，y方向进行等间隔剖分，沿z方向进行灵活剖分，得到计算区域模型。
z方向可以根据场变化的快慢灵活剖分，在场变化快的地方网格间隔Δz适当变小，远离场的地方垂向网格间隔适当稀疏，在保证计算精度的前提下提高了计算效率。
步骤104，对异常体进行网格剖分，得到异常体网格模型，并在异常体网格模型中设置异常体密度值，得到异常体密度分布源函数。
其中，异常体网格模型与计算区域模型对应部分的网格间隔一致。
步骤106，根据计算区域模型、异常体密度分布源函数，得到任一观测高度的重力梯度核函数积分系数。
步骤108，根据重力梯度核函数积分系数，得到任一观测高度上任一层异常体对应的核函数系数矩阵。
通过对计算区域和异常体分别进行网格剖分，并且计算区域网格和异常体网格间隔一致，本发明网格的构建方式使得在得到核函数系数矩阵时，避免了传统二维离散卷积算法核函数系数矩阵中计算区域和异常体剖分网格数必须一致，并且需要用零元素填充以形成核函数系数矩阵的问题。由于异常体网格数通常比计算区域网格数少很多，本发明的方法大大减少了核函数系数矩阵元素数的存储与计算。
步骤110，根据核函数系数矩阵和预设的异常体密度值扩展矩阵，通过二维傅里叶变换实现快速离散卷积计算，并进一步进行傅里叶反变换，得到任一观测高度上任一层异常体对应的重力梯度空间域卷积结果。
步骤112，将各层异常体对应的重力梯度空间域卷积结果进行累加，得到计算区域中任一观测高度上重力梯度各个分量的值。
上述大规模复杂地质体重力梯度快速计算方法中，构建三维坐标系确定计算区域后，将计算区域沿x，y方向进行等间隔剖分，沿z方向进行灵活剖分，得到计算区域模型；对异常体进行网格剖分，其与计算区域模型的网格间隔一致，并得到异常体密度分布源函数；计算任一观测高度的重力梯度核函数积分系数，并进一步处理得到任一观测高度上任一层异常体对应的核函数系数矩阵；根据核函数系数矩阵和预设的异常体密度值扩展矩阵，通过二维傅里叶变换实现快速离散卷积计算，并进一步进行傅里叶反变换，得到任一观测高度上任一层异常体对应的重力梯度空间域卷积结果；将各层异常体对应的重力梯度空间域卷积结果进行累加，得到计算区域中任一观测高度上重力梯度各个分量的值。本发明通过对计算区域和异常体分别进行网格剖分，避免了传统二维离散卷积算法核函数系数矩阵中计算区域和异常体剖分网格数必须一致，并且需要用零元素填充以形成核函数系数矩阵的问题，不仅可以处理更一般的卷积情况，并且可以减少核函数系数矩阵元素数存储与计算。
在其中一个实施例中，还包括：根据待探测的异常体的分布构建三维坐标系，x轴方向指向正东，y轴方向指向正北，z轴垂直向下；确定计算区域范围；将计算区域沿x，y方向进行等间隔剖分，沿z方向根据场变化的快慢进行灵活剖分，在场变化快的区域减小间隔，在场变化慢的区域增大间隔，得到计算区域模型。
在其中一个实施例中，还包括：根据计算区域模型、异常体密度分布源函数，得到观测高度
Figure DEST_PATH_725185DEST_PATH_IMAGE001
上的重力梯度核函数积分系数为：
Figure 172255DEST_PATH_IMAGE056
Figure 884996DEST_PATH_IMAGE057
Figure 356560DEST_PATH_IMAGE058
其中，
Figure 990803DEST_PATH_IMAGE059
分别表示
Figure 88072DEST_PATH_IMAGE060
的核函数积分系数；
Figure 920899DEST_PATH_IMAGE061
表示万有引力常数；
Figure 78211DEST_PATH_IMAGE062
表示编号为
Figure 266264DEST_PATH_IMAGE063
的计算区域模型单元的中心坐标，
Figure 167224DEST_PATH_IMAGE064
表示编号为
Figure 792240DEST_PATH_IMAGE065
的异常体网格模型单元的中心坐标，
Figure 120453DEST_PATH_IMAGE066
，
Figure 260448DEST_PATH_IMAGE067
,
Figure 981410DEST_PATH_IMAGE068
，
Figure 523250DEST_PATH_IMAGE069
,
Figure 22364DEST_PATH_IMAGE070
，
Figure 649655DEST_PATH_IMAGE071
和
Figure 95680DEST_PATH_IMAGE072
分别为计算区域模型在x，y方向上的网格数，
Figure 741293DEST_PATH_IMAGE073
和
Figure 411309DEST_PATH_IMAGE074
分别为异常体网格模型在x，y，z方向上的网格数，
Figure 791475DEST_PATH_IMAGE075
，
Figure 572349DEST_PATH_IMAGE076
，
Figure 839513DEST_PATH_IMAGE077
，
Figure 883693DEST_PATH_IMAGE078
、
Figure 485575DEST_PATH_IMAGE079
、
Figure 70141DEST_PATH_IMAGE080
、
Figure 706658DEST_PATH_IMAGE081
、
Figure 970674DEST_PATH_IMAGE082
，
Figure 794273DEST_PATH_IMAGE083
和
Figure 182529DEST_PATH_IMAGE084
分别为x，y，z方向上的网格间隔，
Figure 611237DEST_PATH_IMAGE085
。
在其中一个实施例中，还包括：根据重力梯度核函数积分系数，得到观测高度
Figure 793956DEST_PATH_IMAGE086
上第n层异常体对应的核函数系数矩阵为：
Figure 386743DEST_PATH_IMAGE087
Figure 313111DEST_PATH_IMAGE088
Figure 393062DEST_PATH_IMAGE089
。
给出的核函数系数矩阵的大小为
Figure 12262DEST_PATH_IMAGE090
，而并非传统二维离散卷积算法核函数系数矩阵的
Figure 13716DEST_PATH_IMAGE091
和
Figure 993042DEST_PATH_IMAGE092
，
Figure 927500DEST_PATH_IMAGE093
和
Figure 983181DEST_PATH_IMAGE094
必须大小相同，并且在形成核函数系数矩阵必须用零元素填充为大小为
Figure 268669DEST_PATH_IMAGE095
的矩阵，该方法不仅可以处理更一般的卷积情况，并且可以减少核函数系数矩阵元素数存储与计算。
在其中一个实施例中，还包括：获取第n层异常体的密度值
Figure 287572DEST_PATH_IMAGE096
，表示为矩阵
Figure 545378DEST_PATH_IMAGE097
；
将矩阵
Figure 506380DEST_PATH_IMAGE097
填零扩展，得到异常体密度值扩展矩阵：
Figure 13585DEST_PATH_IMAGE098
通过二维傅里叶变换实现快速离散卷积计算，并进一步进行傅里叶反变换，得到任一观测高度上任一层异常体对应的重力梯度空间域卷积结果。
在其中一个实施例中，还包括：通过二维傅里叶变换实现快速离散卷积计算，即核函数系数矩阵与异常体密度值扩展矩阵的快速相乘：
Figure 351026DEST_PATH_IMAGE099
Figure 771993DEST_PATH_IMAGE100
Figure 638318DEST_PATH_IMAGE101
其中，
Figure 898398DEST_PATH_IMAGE102
表示二维离散傅里叶正变换算子；
通过二维离散傅里叶反变换，得到观测高度
Figure 711633DEST_PATH_IMAGE103
上第n层异常体对应的重力梯度空间域卷积结果：
Figure 740769DEST_PATH_IMAGE104
Figure 528727DEST_PATH_IMAGE105
Figure 541683DEST_PATH_IMAGE106
其中，
Figure 955347DEST_PATH_IMAGE107
表示二维离散傅里叶反变换算子。
在其中一个实施例中，还包括：将各层异常体对应的重力梯度空间域卷积结果进行累加，得到计算区域中观测高度
Figure 838989DEST_PATH_IMAGE108
上重力梯度各个分量的值为：
Figure 515958DEST_PATH_IMAGE109
Figure 734319DEST_PATH_IMAGE110
Figure 686094DEST_PATH_IMAGE111
其中，
Figure 955401DEST_PATH_IMAGE112
为计算区域在z方向上剖分的网格数。
第n层二维离散卷积的前
Figure 334430DEST_PATH_IMAGE113
行和
Figure 541552DEST_PATH_IMAGE114
列，即为第n层长方体组合模型重力梯度
Figure 297018DEST_PATH_IMAGE115
，将各层组合模型的重力梯度进行累加可以得到整个模型的重力梯度张量。
在其中一个实施例中，如图2所示，提供了一种大规模复杂地质体重力梯度快速计算方法，包括：
应用场景设置；
核函数积分系数计算：给出重力梯度张量核函数计算公式计算积分系数；
核函数系数矩阵预处理：对某一观测平面核函数系数矩阵进行处理，使其适用于一般情况；
二维快速离散卷积：采用二维快速离散卷积实现核函数与密度的快速相乘；
输出重力梯度各分量：对各层模型的重力梯度累加得到整个模型的重力梯度张量。
应该理解的是，虽然图1、2图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，图1、2少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。
在一个具体实施例中，如图3所示，计算区域内有一个棱柱型异常体，计算区域范围为：x和y方向均从-1000m到1000m，z方向从0m到1000m(z轴垂直向下为正)，网格剖分间隔均为10m，整个计算区域剖分为200×200×200个单元，异常体分布范围为：x和y方向均从-100m到100m，z方向从150m到350m，剩余密度为2000kg/m3，计算整个三维网格上的重力异常。
利用Fortran语言编程实现，运行程序所用的个人电脑配置为：CPU-Inter Corei7-8700，主频为3.2GHz。计算了地面重力梯度分量gxx、gyy和gzz异常值。图4、图5分别为本发明方法和解析解计算的重力梯度分量gxx结果图，从两幅图可以看出解析解与本发明方法计算结果形态吻合很好；图6为解析式法与本发明方法计算的相对误差，可以看出整个平面相对误差都小于4.0×10-9；图7、图8分别为本发明方法和解析解计算的重力梯度分量gyy的平面等值线图；图9为相对误差图，计算结果与gxx分量类似。图10、图11分别为本发明方法和解析解计算的重力梯度分量gzz结果图，从两幅图可以看出解析解与本发明方法计算结果形态吻合很好；图12为解析式法与本发明方法计算的相对误差，可以看出整个平面相对误差都小于2.0×10-10，可以看出本发明方法计算精度较高。
图13给出了本发明方法和传统二维离散卷积算法的计算时间随网格剖分数目变化的曲线。从图中可以看出，随着网格剖分数目的增加，传统二维离散卷积算法计算时间急剧增加，而本发明方法的计算时间增长缓慢，总体可以看出随着网格剖分数目的增加，本发明方法的优势越明显，进一步表明本方法更适用于大规模地质体的快速计算。
在一个实施例中，如图14，提供了一种大规模复杂地质体重力梯度快速计算装置，包括：应用场景设置模块1402、重力梯度核函数积分系数确定模块1404、核函数系数矩阵确定模块1406、二维快速离散卷积模块1408和重力梯度张量各分量输出模块1410，其中：
应用场景设置模块1402，用于根据待探测的异常体的分布构建三维坐标系，并确定计算区域，将计算区域沿x，y方向进行等间隔剖分，沿z方向进行灵活剖分，得到计算区域模型；对异常体进行网格剖分，得到异常体网格模型，并在异常体网格模型中设置异常体密度值，得到异常体密度分布源函数；其中，异常体网格模型与计算区域模型对应部分的网格间隔一致；
重力梯度核函数积分系数确定模块1404，用于根据计算区域模型、异常体密度分布源函数，得到任一观测高度的重力梯度核函数积分系数；
核函数系数矩阵确定模块1406，用于根据重力梯度核函数积分系数，得到任一观测高度上任一层异常体对应的核函数系数矩阵；
二维快速离散卷积模块1408，用于根据核函数系数矩阵和预设的异常体密度值扩展矩阵，通过二维傅里叶变换实现快速离散卷积计算，并进一步进行傅里叶反变换，得到任一观测高度上任一层异常体对应的重力梯度空间域卷积结果；
重力梯度张量各分量输出模块1410，用于将各层异常体对应的重力梯度空间域卷积结果进行累加，得到计算区域中任一观测高度上重力梯度各个分量的值。
应用场景设置模块1402还用于根据待探测的异常体的分布构建三维坐标系，x轴方向指向正东，y轴方向指向正北，z轴垂直向下；确定计算区域范围；将计算区域沿x，y方向进行等间隔剖分，沿z方向根据场变化的快慢进行灵活剖分，在场变化快的区域减小间隔，在场变化慢的区域增大间隔，得到计算区域模型。
重力梯度核函数积分系数确定模块1404还用于根据计算区域模型、异常体密度分布源函数，得到观测高度
Figure 879238DEST_PATH_725185DEST_PATH_IMAGE001
上的重力梯度核函数积分系数为：
Figure 174024DEST_PATH_IMAGE117
Figure 117709DEST_PATH_IMAGE118
Figure 663485DEST_PATH_IMAGE119
其中，
Figure 641805DEST_PATH_IMAGE120
分别表示
Figure 628216DEST_PATH_IMAGE121
的核函数积分系数；
Figure 262459DEST_PATH_IMAGE122
表示万有引力常数；
Figure 359728DEST_PATH_IMAGE123
表示编号为
Figure 677708DEST_PATH_IMAGE124
的计算区域模型单元的中心坐标，
Figure 100599DEST_PATH_IMAGE125
表示编号为
Figure 18877DEST_PATH_IMAGE126
的异常体网格模型单元的中心坐标，
Figure 654258DEST_PATH_IMAGE127
，
Figure 544853DEST_PATH_IMAGE128
,
Figure 122334DEST_PATH_IMAGE129
，
Figure 262328DEST_PATH_IMAGE130
,
Figure 232558DEST_PATH_IMAGE131
，
Figure 774398DEST_PATH_IMAGE132
和
Figure 24245DEST_PATH_IMAGE133
分别为计算区域模型在x，y方向上的网格数，
Figure 854798DEST_PATH_IMAGE134
和
Figure 831981DEST_PATH_IMAGE135
分别为异常体网格模型在x，y，z方向上的网格数，
Figure 493906DEST_PATH_IMAGE136
，
Figure 429501DEST_PATH_IMAGE137
，
Figure 790426DEST_PATH_IMAGE138
，
Figure 571300DEST_PATH_IMAGE139
、
Figure 87732DEST_PATH_IMAGE140
、
Figure 131911DEST_PATH_IMAGE141
、
Figure 733794DEST_PATH_IMAGE142
、
Figure 334671DEST_PATH_IMAGE143
，
Figure 705609DEST_PATH_IMAGE144
和
Figure 717428DEST_PATH_IMAGE145
分别为x，y，z方向上的网格间隔，
Figure 541027DEST_PATH_IMAGE146
。
核函数系数矩阵确定模块1406还用于根据重力梯度核函数积分系数，得到观测高度
Figure 132545DEST_PATH_IMAGE147
上第n层异常体对应的核函数系数矩阵为：
Figure 607258DEST_PATH_IMAGE148
Figure 789978DEST_PATH_IMAGE149
Figure 632032DEST_PATH_IMAGE150
。
二维快速离散卷积模块1408还用于获取第n层异常体的密度值
Figure 558399DEST_PATH_IMAGE151
，表示为矩阵
Figure 389083DEST_PATH_IMAGE152
；
将矩阵
Figure 211546DEST_PATH_IMAGE152
填零扩展，得到异常体密度值扩展矩阵：
Figure 9738DEST_PATH_IMAGE153
通过二维傅里叶变换实现快速离散卷积计算，并进一步进行傅里叶反变换，得到任一观测高度上任一层异常体对应的重力梯度空间域卷积结果。
二维快速离散卷积模块1408还用于通过二维傅里叶变换实现快速离散卷积计算，即核函数系数矩阵与异常体密度值扩展矩阵的快速相乘：
Figure 474217DEST_PATH_IMAGE154
其中，
Figure 939833DEST_PATH_IMAGE155
表示二维离散傅里叶正变换算子；
通过二维离散傅里叶反变换，得到观测高度
Figure 982132DEST_PATH_IMAGE156
上第n层异常体对应的重力梯度空间域卷积结果：
Figure 2041DEST_PATH_IMAGE157
Figure 739053DEST_PATH_IMAGE158
Figure 793596DEST_PATH_IMAGE159
其中，
Figure 489020DEST_PATH_IMAGE160
表示二维离散傅里叶反变换算子。
重力梯度张量各分量输出模块1410还用于将各层异常体对应的重力梯度空间域卷积结果进行累加，得到计算区域中观测高度
Figure 278115DEST_PATH_IMAGE161
上重力梯度各个分量的值为：
Figure 349977DEST_PATH_IMAGE162
Figure 259027DEST_PATH_IMAGE163
Figure 390931DEST_PATH_IMAGE164
其中，
Figure 588694DEST_PATH_IMAGE165
为计算区域在z方向上剖分的网格数。
关于大规模复杂地质体重力梯度快速计算装置的具体限定可以参见上文中对于大规模复杂地质体重力梯度快速计算方法的限定，在此不再赘述。上述大规模复杂地质体重力梯度快速计算装置中的各个模块可全部或部分通过软件、硬件及其组合来实现。上述各模块可以硬件形式内嵌于或独立于计算机设备中的处理器中，也可以以软件形式存储于计算机设备中的存储器中，以便于处理器调用执行以上各个模块对应的操作。
在一个实施例中，提供了一种计算机设备，该计算机设备可以是终端，其内部结构图可以如图15所示。该计算机设备包括通过系统总线连接的处理器、存储器、网络接口、显示屏和输入装置。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现一种大规模复杂地质体重力梯度快速计算方法。该计算机设备的显示屏可以是液晶显示屏或者电子墨水显示屏，该计算机设备的输入装置可以是显示屏上覆盖的触摸层，也可以是计算机设备外壳上设置的按键、轨迹球或触控板，还可以是外接的键盘、触控板或鼠标等。
本领域技术人员可以理解，图15中示出的结构，仅仅是与本申请方案相关的部分结构的框图，并不构成对本申请方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。
在一个实施例中，提供了一种计算机设备，包括存储器和处理器，该存储器存储有计算机程序，该处理器执行计算机程序时实现上述方法实施例中的步骤。
在一个实施例中，提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述方法实施例中的步骤。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器（ROM）、可编程ROM（PROM）、电可编程ROM（EPROM）、电可擦除可编程ROM（EEPROM）或闪存。易失性存储器可包括随机存取存储器（RAM）或者外部高速缓冲存储器。作为说明而非局限，RAM以多种形式可得，诸如静态RAM（SRAM）、动态RAM（DRAM）、同步DRAM（SDRAM）、双数据率SDRAM（DDRSDRAM）、增强型SDRAM（ESDRAM）、同步链路（Synchlink） DRAM（SLDRAM）、存储器总线（Rambus）直接RAM（RDRAM）、直接存储器总线动态RAM（DRDRAM）、以及存储器总线动态RAM（RDRAM）等。
以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求为准。

Claims (8)

1.一种大规模复杂地质体重力梯度快速计算方法，其特征在于，所述方法包括：
根据待探测的异常体的分布构建三维坐标系，并确定计算区域，将所述计算区域沿x，y方向进行等间隔剖分，沿z方向进行灵活剖分，得到计算区域模型；
对所述异常体进行网格剖分，得到异常体网格模型，并在所述异常体网格模型中设置异常体密度值，得到异常体密度分布源函数；其中，异常体网格模型与所述计算区域模型对应部分的网格间隔一致；
根据所述计算区域模型、所述异常体密度分布源函数，得到观测高度
Figure 725185DEST_PATH_IMAGE001
上的重力梯度核函数积分系数为：
Figure 237069DEST_PATH_IMAGE002
Figure 983308DEST_PATH_IMAGE003
Figure 806908DEST_PATH_IMAGE004
其中，
Figure 70530DEST_PATH_IMAGE005
分别表示
Figure 295975DEST_PATH_IMAGE006
的核函数积分系数；
Figure 354061DEST_PATH_IMAGE007
表示万有引力常数；
Figure 664956DEST_PATH_IMAGE008
表示编号为
Figure 466690DEST_PATH_IMAGE009
的计算区域模型单元的中心坐标，
Figure 546642DEST_PATH_IMAGE010
表示编号为
Figure 41208DEST_PATH_IMAGE011
的异常体网格模型单元的中心坐标，
Figure 839400DEST_PATH_IMAGE012
，
Figure 303879DEST_PATH_IMAGE013
，
Figure 113703DEST_PATH_IMAGE015
Figure 903805DEST_PATH_IMAGE016
和
Figure 85167DEST_PATH_IMAGE017
分别为所述计算区域模型在x，y方向上的网格数，
Figure 87758DEST_PATH_IMAGE018
、
Figure 17668DEST_PATH_IMAGE019
和
Figure 978670DEST_PATH_IMAGE020
分别为所述异常体网格模型在x，y，z方向上的网格数，
Figure 220296DEST_PATH_IMAGE021
，
Figure 167523DEST_PATH_IMAGE022
，
Figure 342153DEST_PATH_IMAGE023
，
Figure 83844DEST_PATH_IMAGE024
、
Figure 78345DEST_PATH_IMAGE025
、
Figure 563684DEST_PATH_IMAGE026
、
Figure 592820DEST_PATH_IMAGE027
，
Figure 364467DEST_PATH_IMAGE028
Figure 987209DEST_PATH_IMAGE029
和
Figure 400873DEST_PATH_IMAGE030
分别为x，y，z方向上的网格间隔，
Figure 159881DEST_PATH_IMAGE031
；
根据所述重力梯度核函数积分系数，得到任一观测高度上任一层异常体对应的核函数系数矩阵；
根据所述核函数系数矩阵和预设的异常体密度值扩展矩阵，通过二维傅里叶变换实现快速离散卷积计算，并进一步进行傅里叶反变换，得到任一观测高度上任一层异常体对应的重力梯度空间域卷积结果；
将各层异常体对应的重力梯度空间域卷积结果进行累加，得到所述计算区域中任一观测高度上重力梯度各个分量的值。
2.根据权利要求1所述的方法，其特征在于，根据待探测的异常体的分布构建三维坐标系，并确定计算区域，将所述计算区域沿x，y方向进行等间隔剖分，沿z方向进行灵活剖分，得到计算区域模型，包括：
根据待探测的异常体的分布构建三维坐标系，x轴方向指向正东，y轴方向指向正北，z轴垂直向下；
确定计算区域范围；
将所述计算区域沿x，y方向进行等间隔剖分，沿z方向根据场变化的快慢进行灵活剖分，在场变化快的区域减小间隔，在场变化慢的区域增大间隔，得到计算区域模型。
3.根据权利要求2所述的方法，其特征在于，根据所述重力梯度核函数积分系数，得到任一观测高度上任一层异常体对应的核函数系数矩阵，包括：
根据所述重力梯度核函数积分系数，得到观测高度
Figure 725185DEST_PATH_IMAGE001
上第n层异常体对应的核函数系数矩阵为：
Figure 212468DEST_PATH_IMAGE034
Figure 164243DEST_PATH_IMAGE036
Figure 902392DEST_PATH_IMAGE038
。
4.根据权利要求3所述的方法，其特征在于，根据所述核函数系数矩阵和预设的异常体密度值扩展矩阵，通过二维傅里叶变换实现快速离散卷积计算，并进一步进行傅里叶反变换，得到任一观测高度上任一层异常体对应的重力梯度空间域卷积结果，包括：
获取第n层异常体的密度值
Figure 153858DEST_PATH_IMAGE039
，表示为矩阵
Figure 610247DEST_PATH_IMAGE040
；
将矩阵
Figure 975500DEST_PATH_IMAGE041
填零扩展，得到异常体密度值扩展矩阵：
Figure 833735DEST_PATH_IMAGE043
通过二维傅里叶变换实现快速离散卷积计算，并进一步进行傅里叶反变换，得到任一观测高度上任一层异常体对应的重力梯度空间域卷积结果。
5.根据权利要求4所述的方法，其特征在于，所述通过二维傅里叶变换实现快速离散卷积计算，并进一步进行傅里叶反变换，得到任一观测高度上任一层异常体对应的重力梯度空间域卷积结果，包括：
通过二维傅里叶变换实现快速离散卷积计算，即所述核函数系数矩阵与所述异常体密度值扩展矩阵的快速相乘：
Figure 524610DEST_PATH_IMAGE045
Figure 202716DEST_PATH_IMAGE047
Figure 496294DEST_PATH_IMAGE049
其中，
Figure 84402DEST_PATH_IMAGE050
表示二维离散傅里叶正变换算子；
通过二维离散傅里叶反变换，得到观测高度
Figure 725185DEST_PATH_IMAGE001
上第n层异常体对应的重力梯度空间域卷积结果：
Figure 111580DEST_PATH_IMAGE052
Figure 208849DEST_PATH_IMAGE053
Figure 917042DEST_PATH_IMAGE054
其中，
Figure 74354DEST_PATH_IMAGE055
表示二维离散傅里叶反变换算子。
6.根据权利要求5所述的方法，其特征在于，所述将各层异常体对应的重力梯度空间域卷积结果进行累加，得到所述计算区域中任一观测高度上重力梯度各个分量的值，包括：
将各层异常体对应的重力梯度空间域卷积结果进行累加，得到所述计算区域中观测高度
Figure 725185DEST_PATH_IMAGE001
上重力梯度各个分量的值为：
Figure 237799DEST_PATH_IMAGE057
Figure 800499DEST_PATH_IMAGE058
Figure 863133DEST_PATH_IMAGE059
其中，
Figure 3127DEST_PATH_IMAGE060
为所述计算区域在z方向上剖分的网格数。
7.一种大规模复杂地质体重力梯度快速计算装置，其特征在于，所述装置包括：
应用场景设置模块，用于根据待探测的异常体的分布构建三维坐标系，并确定计算区域，将所述计算区域沿x，y方向进行等间隔剖分，沿z方向进行灵活剖分，得到计算区域模型；对所述异常体进行网格剖分，得到异常体网格模型，并在所述异常体网格模型中设置异常体密度值，得到异常体密度分布源函数；其中，异常体网格模型与所述计算区域模型对应部分的网格间隔一致；
重力梯度核函数积分系数确定模块，用于根据所述计算区域模型、所述异常体密度分布源函数，得到任一观测高度的重力梯度核函数积分系数；
根据所述计算区域模型、所述异常体密度分布源函数，得到观测高度
Figure 580215DEST_PATH_IMAGE062
上的重力梯度核函数积分系数为：
Figure 856475DEST_PATH_IMAGE063
其中，
Figure 230956DEST_PATH_IMAGE064
分别表示
Figure 592667DEST_PATH_IMAGE065
的核函数积分系数；
Figure 569850DEST_PATH_IMAGE066
表示万有引力常数；
Figure 107142DEST_PATH_IMAGE067
表示编号为
Figure 511578DEST_PATH_IMAGE068
的计算区域模型单元的中心坐标，
Figure 501531DEST_PATH_IMAGE070
表示编号为
Figure 282405DEST_PATH_IMAGE072
的异常体网格模型单元的中心坐标，
Figure 533258DEST_PATH_IMAGE073
，
Figure 249541DEST_PATH_IMAGE074
,
Figure 851424DEST_PATH_IMAGE075
，
Figure 311355DEST_PATH_IMAGE076
,
Figure 416715DEST_PATH_IMAGE077
，
Figure 303899DEST_PATH_IMAGE078
和
Figure 127499DEST_PATH_IMAGE079
分别为所述计算区域模型在x，y方向上的网格数，
Figure 391121DEST_PATH_IMAGE080
和
Figure 350987DEST_PATH_IMAGE081
分别为所述异常体网格模型在x，y，z方向上的网格数，
Figure 533706DEST_PATH_IMAGE082
，
Figure 982618DEST_PATH_IMAGE083
，
Figure 643406DEST_PATH_IMAGE084
，
Figure 598724DEST_PATH_IMAGE085
、
Figure 217924DEST_PATH_IMAGE086
、
Figure 750537DEST_PATH_IMAGE087
、
Figure 90382DEST_PATH_IMAGE027
、
Figure 24840DEST_PATH_IMAGE088
，
Figure 690308DEST_PATH_IMAGE089
和
Figure 710216DEST_PATH_IMAGE090
分别为x，y，z方向上的网格间隔，
Figure 119332DEST_PATH_IMAGE092
；
二维快速离散卷积模块，用于根据所述核函数系数矩阵和预设的异常体密度值扩展矩阵，通过二维傅里叶变换实现快速离散卷积计算，并进一步进行傅里叶反变换，得到任一观测高度上任一层异常体对应的重力梯度空间域卷积结果；
重力梯度张量各分量输出模块，用于将各层异常体对应的重力梯度空间域卷积结果进行累加，得到所述计算区域中任一观测高度上重力梯度各个分量的值。
8.一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至6中任一项所述方法的步骤。
CN202111372974.9A 2021-11-19 2021-11-19 大规模复杂地质体重力梯度快速计算方法、装置和设备 Active CN113806686B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111372974.9A CN113806686B (zh) 2021-11-19 2021-11-19 大规模复杂地质体重力梯度快速计算方法、装置和设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111372974.9A CN113806686B (zh) 2021-11-19 2021-11-19 大规模复杂地质体重力梯度快速计算方法、装置和设备

Publications (2)

Publication Number Publication Date
CN113806686A CN113806686A (zh) 2021-12-17
CN113806686B true CN113806686B (zh) 2022-02-08

Family

ID=78938412

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111372974.9A Active CN113806686B (zh) 2021-11-19 2021-11-19 大规模复杂地质体重力梯度快速计算方法、装置和设备

Country Status (1)

Country Link
CN (1) CN113806686B (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107577641A (zh) * 2017-08-21 2018-01-12 吉林大学 一种基于gpu并行的重力梯度张量数据快速密度反演方法
CN112800657A (zh) * 2021-04-15 2021-05-14 中南大学 基于复杂地形的重力场数值模拟方法、装置和计算机设备
CN113138426A (zh) * 2021-04-12 2021-07-20 中南大学 一种基于多变量核密度欧拉解概率密度成像方法
CN113420487A (zh) * 2021-08-25 2021-09-21 中南大学 一种三维重力梯度张量数值模拟方法、装置、设备和介质
CN113514899A (zh) * 2021-07-12 2021-10-19 吉林大学 一种重力梯度的自适应剖分反演方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10242126B2 (en) * 2012-01-06 2019-03-26 Technoimaging, Llc Method of simultaneous imaging of different physical properties using joint inversion of multiple datasets

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107577641A (zh) * 2017-08-21 2018-01-12 吉林大学 一种基于gpu并行的重力梯度张量数据快速密度反演方法
CN113138426A (zh) * 2021-04-12 2021-07-20 中南大学 一种基于多变量核密度欧拉解概率密度成像方法
CN112800657A (zh) * 2021-04-15 2021-05-14 中南大学 基于复杂地形的重力场数值模拟方法、装置和计算机设备
CN113514899A (zh) * 2021-07-12 2021-10-19 吉林大学 一种重力梯度的自适应剖分反演方法
CN113420487A (zh) * 2021-08-25 2021-09-21 中南大学 一种三维重力梯度张量数值模拟方法、装置、设备和介质

Also Published As

Publication number Publication date
CN113806686A (zh) 2021-12-17

Similar Documents

Publication Publication Date Title
CN112287534B (zh) 基于nufft的二维磁异常快速正演模拟方法和装置
CN112800657B (zh) 基于复杂地形的重力场数值模拟方法、装置和计算机设备
Jahandari et al. Forward modeling of gravity data using finite-volume and finite-element methods on unstructured grids
CN113420487B (zh) 一种三维重力梯度张量数值模拟方法、装置、设备和介质
CN113962077B (zh) 三维各向异性强磁场数值模拟方法、装置、设备及介质
CN111967169B (zh) 二度体重力异常积分解数值模拟方法和装置
Larios et al. Continuous data assimilation with a moving cluster of data points for a reaction diffusion equation: A computational study
CN114021408A (zh) 二维强磁场数值模拟方法、装置、设备及介质
D'Urso et al. On a generalized Love's problem
CN116842813B (zh) 一种三维大地电磁正演数值模拟方法
CN113806686B (zh) 大规模复杂地质体重力梯度快速计算方法、装置和设备
CN113673163B (zh) 一种三维磁异常数快速正演方法、装置和计算机设备
CN113076678B (zh) 一种频率域二度体重力异常快速数值模拟方法和装置
CN113779818B (zh) 三维地质体其电磁场数值模拟方法、装置、设备及介质
CN113268702B (zh) 一种频率域磁梯度张量变换方法、装置和计算机设备
CN113221409B (zh) 一种有限元和边界元耦合的声波二维数值模拟方法和装置
CN113642189B (zh) 一种基于积分解的重力梯度张量快速数值模拟方法和装置
CN113656976B (zh) 一种二维磁梯度张量快速数值模拟方法、装置和设备
Yin et al. A fast 3D gravity forward algorithm based on circular convolution
CN114036806A (zh) 基于热导率各向异性介质的三维地温场数值模拟方法
CN111324972A (zh) 一种基于gpu并行的探地雷达电磁波数值模拟计算方法
CN114065528A (zh) 一种二度体重力异常快速数值模拟方法、装置和设备
Browne et al. Nonlinear solution techniques for solving a Monge-Amp\ere equation for redistribution of a mesh
CN114002749A (zh) 二度体井中磁场梯度张量计算方法、装置、设备及介质
CN114004092A (zh) 二度体井中重力梯度张量计算方法、装置、设备及介质

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