CN111967169A - 二度体重力异常积分解数值模拟方法和装置 - Google Patents

二度体重力异常积分解数值模拟方法和装置 Download PDF

Info

Publication number
CN111967169A
CN111967169A CN202011127944.7A CN202011127944A CN111967169A CN 111967169 A CN111967169 A CN 111967169A CN 202011127944 A CN202011127944 A CN 202011127944A CN 111967169 A CN111967169 A CN 111967169A
Authority
CN
China
Prior art keywords
wave number
domain
sampling
value
gravity anomaly
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
CN202011127944.7A
Other languages
English (en)
Other versions
CN111967169B (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 CN202011127944.7A priority Critical patent/CN111967169B/zh
Publication of CN111967169A publication Critical patent/CN111967169A/zh
Application granted granted Critical
Publication of CN111967169B publication Critical patent/CN111967169B/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/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

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

Abstract

本发明提出了一种二度体重力异常积分解数值模拟方法、装置、计算机设备和存储介质。本发明通过重力模型表示、波数域重力异常公式推导、波数域非均匀采样、空间域重力异常场计算等步骤，利用非均匀采样快速傅里叶变换实现了重力异常场高效、高精度数值模拟。在本发明中，将非均匀采样傅里叶变换应用于重力异常模拟，充分融合了快速傅里叶变换的效率优势和高斯傅里叶变换的精度优势，有效兼顾了重力异常数值模拟的计算精度与计算效率，解决了波数域重力异常数值模拟只适合规则测线、数值模拟方法不能同时兼顾计算效率和计算精度的问题。

Description

二度体重力异常积分解数值模拟方法和装置
技术领域
本申请涉及二度体重力异常计算技术领域，特别是涉及一种二度体重力异常积分解数值模拟方法、装置、计算机设备和存储介质。
背景技术
重力勘探是一种以地下介质密度差异作为物性基础的一种地球物理勘探方法。该方法施工方便、成本低、效率高、勘探深度大，已广泛应用于地球深部构造研究、大地及区域地质构造研究、构造单元划分、隐伏岩体或岩层及断裂的探测、成矿远景区的划分、石油、天然气或煤等远景盆地圈定等方面。在实际生产中，常把走向方向尺度远比垂直走向方向尺度大的地质体用走向方向无限延伸的二度体代替。针对走向方向无限延伸的二度体模型，高效、高精度数值模拟作为反演成像与定量解释的基础，一直是研究的焦点。
目前，有关二度体的数值模拟方法主要分为空间域方法和波数域方法。空间域方法是通过异常场解析式直接求取空间任意点的精确异常场，其优点是原理简单、结果精确，缺点是解析式比较复杂，推导过程繁琐，且当计算大量位置点异常场数据时，速度较慢。波数域方法是根据场源产生异常场的傅里叶变换解析表达式，通过数值方法计算该异常频谱的反傅里叶变换得到空间域异常场。相对于空间域方法而言，其优点是表达式简洁，计算效率较高。随着快速傅里叶变换扩边法和高斯快速傅里叶变换法Gauss-FFT在重力异常正演数值模拟中的广泛应用，频率域方法以其简洁性、准确性和高效性成为处理大规模复杂模型正演的首选方法。文献(Tontini, F.C., L. Cocchi, C. Carmisciano. Rapid 3-Dforward model of potential fields with application to the Palinuro Seamountmagnetic anomaly (southern Tyrrhenian Sea, Italy). Journal of GeophysicalResearch, 2009. 114.) 利用傅里叶变换，详细推导了全空间的波数域表达式，借助快速傅里叶变换算法，实现了重力异常快速数值模拟，但在空间域和波数域只能均匀采样，并且由于截断边界效应的影响，采用扩边的方式使得数值模拟精度相对较低)；文献(Wu, L. ,Tian, G. High-precision Fourier forward modeling of potential fields.Geophysics, 2014, 79(5): G59-G68.) 提出了一种重力异常正演模拟的 Gauss-FFT 方法，该方法有效克服了快速傅里叶变换方法的边界截断效应问题，提高了数值模拟精度，但在空间域和频率域只能均匀采样，并且增加了成倍的计算量，效率相对较低。文献(Lee J Y, Greengard L . The type 3 nonuniform FFT and its applications. Journal ofComputational Physics, 2005, 206(1):1-5.) 将非均匀快速傅里叶变换(NUFFT)应用于核磁共振模拟成像中，取得了很好的效果。
目前，不管是传统快速傅里叶变换扩边法还是Gauss-FFT法，在空间域和波数域只能均匀采样，只适用于规则测线，且难以兼顾计算精度与计算速度。
发明内容
基于此，有必要针对上述技术问题，提供一种能够兼顾计算精度与计算速度的二度体重力异常积分解数值模拟方法、装置、计算机设备和存储介质。
一种二度体重力异常积分解数值模拟方法，所述方法包括：
根据观测点在笛卡尔坐标系下的位置信息，构建包含目标区域的二度体模型；所述观测点坐标包括x轴方向和z轴方向；所述二度体模型是将所述目标区域分成若干个矩形，所述二度体的截面形状用于确定每个矩形的空间域异常体剩余密度值；
根据所述二度体模型，确定对应的空间域重力异常表达式；所述空间域重力异常表达式包含空间域异常体剩余密度值；
对所述空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式；其中，所述波数域重力异常表达式包含波数域异常体剩余密度值；所述波数域异常体剩余密度值是对所述空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的；
根据所述矩形沿x轴方向的尺寸信息，确定截止频率，根据所述截止频率确定波数域波数采样值范围，根据所述波数采样值范围得到波数采样值；
将所述波数采样值代入所述波数域重力异常表达式，得到波数域重力异常场值；
通过对所述波数域重力异常场值进行一维非均匀快速傅里叶反变换，得到目标区域内任一点的空间域重力异常场值。
在其中一个实施例中，还包括：空间域重力异常表达式为：
Figure 259620DEST_PATH_IMAGE001
式中：
Figure 966545DEST_PATH_IMAGE002
表示空间域重力异常场；
Figure 157354DEST_PATH_IMAGE003
表示万有引力常数；
Figure 850504DEST_PATH_IMAGE004
表示x方向剖分的矩形个数；
Figure 505476DEST_PATH_IMAGE005
表示z方向剖分的矩形个数；
Figure 258669DEST_PATH_IMAGE006
表示所述观测点坐标；
Figure 999091DEST_PATH_IMAGE007
表示编号为
Figure 495932DEST_PATH_IMAGE008
的矩形中心坐标；
Figure 943094DEST_PATH_IMAGE009
表示编号为
Figure 991821DEST_PATH_IMAGE010
的矩形的空间域异常体剩余密度值，
Figure 94906DEST_PATH_IMAGE011
表示矩形的x轴方向尺寸，
Figure 254492DEST_PATH_IMAGE012
矩形的z轴方向尺寸。
在其中一个实施例中，还包括：波数域重力异常表达式为：
Figure 556161DEST_PATH_IMAGE013
式中：
Figure 385576DEST_PATH_IMAGE014
表示波数域重力异常场；k表示波数；
Figure 100591DEST_PATH_IMAGE015
表示波数域异常体剩余密度值；
Figure 939234DEST_PATH_IMAGE016
表示符号函数：
Figure 360989DEST_PATH_IMAGE017
。
在其中一个实施例中，还包括：非均匀采样快速傅里叶变换为：
Figure 751519DEST_PATH_IMAGE018
式中：i为虚数单位，
Figure 563617DEST_PATH_IMAGE019
为给定离散采样点
Figure 737109DEST_PATH_IMAGE020
对应采样点的值，
Figure 341266DEST_PATH_IMAGE021
为计算的离散采样点
Figure 512484DEST_PATH_IMAGE022
傅里叶变换频谱，N表示采样点总数；
其中，所述非均匀快速傅里叶变换的实现步骤为：
根据
Figure 202092DEST_PATH_IMAGE023
临近q个均匀采样点的傅里叶变换基，得到近似非均匀傅里叶变换基
Figure 382537DEST_PATH_IMAGE024
为：
Figure 778883DEST_PATH_IMAGE025
式中：m表示过采样因子，
Figure 511216DEST_PATH_IMAGE026
表示权重因子，
Figure 32327DEST_PATH_IMAGE027
为精度因子，
Figure 813201DEST_PATH_IMAGE028
表示对
Figure 391950DEST_PATH_IMAGE029
取整；
根据采样点的值和权重因子，计算新傅里叶变换基下对应的傅里叶变换系数
Figure 170550DEST_PATH_IMAGE030
：
Figure 303592DEST_PATH_IMAGE031
采用均匀的快速傅里叶变换，得到：
Figure 825840DEST_PATH_IMAGE032
式中，
Figure 931199DEST_PATH_IMAGE033
表示傅立叶变换之后的频谱。
在其中一个实施例中，还包括：将所述目标区域分成若干个矩形，每个矩形的x轴方向矩形尺寸
Figure 5334DEST_PATH_IMAGE034
、z轴方向矩形尺寸
Figure 32196DEST_PATH_IMAGE035
可以相同，也可以不同；在重力异常场变化快的地方加密网格，在重力异常场变化慢的地方稀疏网格。
在其中一个实施例中，还包括：根据所述矩形沿x轴方向的尺寸信息，确定截止频率为：
Figure 217190DEST_PATH_IMAGE036
式中，
Figure 380318DEST_PATH_IMAGE037
表示截止频率；
Figure 297458DEST_PATH_IMAGE038
表示x轴方向矩形尺寸的最小值。
在其中一个实施例中，还包括：
根据所述截止频率确定波数域波数采样值范围为
Figure 670671DEST_PATH_IMAGE039
；
根据所述波数采样值范围和预先设置的采样点总数，在波数域进行非均匀采样，得到波数采样值，包括：
在
Figure 269142DEST_PATH_IMAGE040
上波数依次取值为：
Figure 676990DEST_PATH_IMAGE041
在
Figure 765032DEST_PATH_IMAGE042
上波数依次取值为：
Figure 500907DEST_PATH_IMAGE043
式中，
Figure 27703DEST_PATH_IMAGE044
，N表示预先设置的采样点总数。
一种二度体重力异常积分解数值模拟装置，所述装置包括：
二度体模型构建模块，用于根据观测点在观测点坐标中的位置信息，构建包含目标区域的二度体模型；所述观测点坐标包括x轴方向和z轴方向；所述二度体模型是将所述目标区域分成若干个矩形，所述二度体的截面形状用于确定每个矩形的空间域异常体剩余密度值；
空间域重力异常表达式确定模块，用于根据所述二度体模型，确定对应的空间域重力异常表达式；所述空间域重力异常表达式包含空间域异常体剩余密度值；
波数域重力异常表达式确定模块，用于对所述空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式；其中，所述波数域重力异常表达式包含波数域异常体剩余密度值；所述波数域异常体剩余密度值是对所述空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的；
波数采样值确定模块，用于根据所述矩形沿x轴方向的尺寸信息，确定截止频率，根据所述截止频率确定波数域波数采样值范围，根据所述波数采样值范围得到波数采样值；
波数域重力异常场值确定模块，用于将所述波数采样值代入所述波数域重力异常表达式，得到波数域重力异常场值；
空间域重力异常场值确定模块，用于通过对所述波数域重力异常场值进行一维非均匀快速傅里叶反变换，得到目标区域内任一点的空间域重力异常场值。
一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现以下步骤：
根据观测点在笛卡尔坐标系下的位置信息，构建包含目标区域的二度体模型；所述观测点坐标包括x轴方向和z轴方向；所述二度体模型是将所述目标区域分成若干个矩形，所述二度体的截面形状用于确定每个矩形的空间域异常体剩余密度值；
根据所述二度体模型，确定对应的空间域重力异常表达式；所述空间域重力异常表达式包含空间域异常体剩余密度值；
对所述空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式；其中，所述波数域重力异常表达式包含波数域异常体剩余密度值；所述波数域异常体剩余密度值是对所述空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的；
根据所述矩形沿x轴方向的尺寸信息，确定截止频率，根据所述截止频率确定波数域波数采样值范围，根据所述波数采样值范围和预先设置的采样点总数，在波数域进行非均匀采样，得到波数采样值；
将所述波数采样值代入所述波数域重力异常表达式，得到波数域重力异常场值；
通过对所述波数域重力异常场值进行一维非均匀快速傅里叶反变换，得到目标区域内任一点的空间域重力异常场值。
一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现以下步骤：
根据观测点在笛卡尔坐标系下的位置信息，构建包含目标区域的二度体模型；所述观测点坐标包括x轴方向和z轴方向；所述二度体模型是将所述目标区域分成若干个矩形，所述二度体的截面形状用于确定每个矩形的空间域异常体剩余密度值；
根据所述二度体模型，确定对应的空间域重力异常表达式；所述空间域重力异常表达式包含空间域异常体剩余密度值；
对所述空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式；其中，所述波数域重力异常表达式包含波数域异常体剩余密度值；所述波数域异常体剩余密度值是对所述空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的；
根据所述矩形沿x轴方向的尺寸信息，确定截止频率，根据所述截止频率确定波数域波数采样值范围，根据所述波数采样值范围和预先设置的采样点总数，在波数域进行非均匀采样，得到波数采样值；
将所述波数采样值代入所述波数域重力异常表达式，得到波数域重力异常场值；
通过对所述波数域重力异常场值进行一维非均匀快速傅里叶反变换，得到目标区域内任一点的空间域重力异常场值。
上述二度体重力异常积分解数值模拟方法、装置、计算机设备和存储介质，通过构建包含目标区域的二度体模型，确定对应的空间域重力异常表达式，其中空间域重力异常表达式包含空间域异常体剩余密度值，通过对空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式，其中波数域异常体剩余密度值是对空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的，根据矩形沿x轴方向的尺寸信息，确定截止频率，根据截止频率确定波数域波数采样值范围，根据波数采样值范围和预先设置的采样点总数进行波数域非均匀采样，得到波数采样值，根据波数采样值，代入波数域重力异常表达式中，得到波数域重力异常场值，最后通过对波数域重力异常场值进行一维非均匀快速傅里叶反变换，得到目标区域内任一点的空间域重力异常场值。采用本方法进行二度体重力异常积分解数值模拟，通过借助非均匀采样快速傅里叶变换，能够实现空间域和波数域非均匀采样，可以有效克服快速傅里叶变换边界效应导致计算精度低，又能保证精度的前提下解决Gauss-FFT方法随着高斯点数增加计算效率低的问题。有效兼顾了重力异常数值模拟的计算精度与计算效率，解决了波数域重力异常数值模拟只适合规则测线、数值模拟方法不能同时兼顾计算效率和计算精度的问题。
附图说明
图1为一个实施例中二度体重力异常积分解数值模拟方法的流程示意图；
图2为一个实施例中二度体重力异常积分解数值模拟方法的算法流程图；
图3为一个实施例中实现非均匀快速傅里叶变换的流程示意图；
图4为一个实施例中截面为圆形二度体模型示意图；
图5为一个实施例中波数域非均匀采样模拟结果与理论值对比图；
图6为一个实施例中波数域非均匀采样模拟结果与理论值的相对误差；
图7为另一个实施例中空间域和波数域都非均匀采样模拟结果与理论值对比图；
图8为另一个实施例中空间域和波数域都非均匀采样模拟结果与理论值的相对误差；
图9为另一个实施例中二度体重力异常积分解数值模拟装置的结构框图；
图10为一个实施例中计算机设备的内部结构图。
具体实施方式
为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。
本申请提供的二度体重力异常积分解数值模拟方法，可以应用于如下应用环境中。根据观测点构建包含目标区域的二度体模型，确定对应的空间域重力异常表达式，其中空间域重力异常表达式包含空间域异常体剩余密度值，通过对空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式，其中波数域异常体剩余密度值是对空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的，根据矩形沿x轴方向的尺寸信息，确定截止频率，根据截止频率确定波数域波数采样值范围，根据波数采样值范围和预先设置的采样点总数进行波数域非均匀采样，得到波数采样值，根据波数采样值，代入波数域重力异常表达式中，得到波数域重力异常场值，最后通过对波数域重力异常场值进行一维非均匀快速傅里叶反变换，得到目标区域内任一点的空间域重力异常场值。
在一个实施例中，如图1所示，提供了一种二度体重力异常积分解数值模拟方法，包括以下步骤：
步骤102，根据观测点在观测点坐标中的位置信息，构建包含目标区域的二度体模型。
观测点坐标包括x轴方向和z轴方向；二度体模型是将目标区域分成若干个矩形，矩形的截面形状用于确定每个矩形的空间域异常体剩余密度值。
具有一定走向，且沿走向方向异常值变化不明显的异常，称为二度异常，它是二维空间坐标
Figure 962161DEST_PATH_IMAGE045
的函数，引起二度异常的地质体称为二度体，如层状矿体、脉状矿体、向斜、背斜等。
构建包含目标区域的二度体模型，将目标区域分成若干个矩形，每个小矩形的几何尺寸
Figure 424366DEST_PATH_IMAGE046
,
Figure 506592DEST_PATH_IMAGE047
可以相同，也可以不同，剖分矩形的网格剖分灵活，可以在重力异常场变化快的地方加密网格，在重力异常场变化慢的地方适当稀疏网格。根据二度体的截面形状，对矩形密度进行赋值，每个矩形的密度是一个常数，不同矩形密度值不同，以此表示任意截面形状、任意密度分布的二度体模型。
步骤104，根据二度体模型，确定对应的空间域重力异常表达式；空间域重力异常表达式包含空间域异常体剩余密度值。
由于实际地球内部的物质密度分布非常不均匀，因而实际观测重力值与理论上的正常重力值总是存在着偏差，这种在排除各种干扰因素影响之后，仅仅是由于物质密度分布不匀而引起的重力的变化，就称为重力异常。重力场正演是指根据密度分布计算重力场，反演是指根据观测重力值计算密度分布。正演是反演的基础，正演计算的效率直接影响反演计算的效率，而正演计算精度直接影响反演结果的质量。在二度体重力异常积分解数值模拟方法中，可以根据构建的二度体模型得到对应的空间域重力异常表达式。
步骤106，对空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式；其中，波数域重力异常表达式包含波数域异常体剩余密度值；波数域异常体剩余密度值是对空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的。
由于构建的二度体模型中剖分的矩形尺寸不一定一样，每个矩形对应的密度值不一定一样，因此，空间域的采样可以是非均匀的。通过对空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式，当剖分的矩形尺寸不一的时候，波数域异常体剩余密度值是对空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的。
步骤108，根据矩形沿x轴方向的尺寸信息，确定截止频率，根据截止频率确定波数域波数采样值范围，根据波数采样值范围和预先设置的采样点总数，在波数域进行非均匀采样，得到波数采样值。
理论上，采样点选取越多计算精度越高，但同时计算量也相应增加，效率下降。因此，需要在满足计算精度的前提下尽可能减少采样点总数。通过对不同模型重力异常频谱分析发现，波数采样值绝对值比较小的时候重力异常频谱比较大，随着波数采样值绝对值的增大，重力异常频谱逐渐减小，直到波数采样值绝对值足够大的时候重力异常频谱几乎衰减为零，重力异常频谱几乎衰减为零说明对应频率的重力异常场能量是极低的，因此不需要选取绝对值过大的波数采样值。根据矩形沿x轴方向的尺寸信息，确定截止频率，根据截止频率确定波数域波数采样值范围，可以在满足精度的前提下尽可能减少采样点总数，提高方法的计算效率。
步骤110，将波数采样值代入所述波数域重力异常表达式，得到波数域重力异常场值。
相对于空间域重力异常表达式，波数域重力异常表达式的优点是表达式简洁，计算效率较高。
步骤112，通过对波数域重力异常场值进行一维非均匀快速傅里叶反变换，得到目标区域内任一点的空间域重力异常场值。
通过一维非均匀快速傅里叶反变换，将波数域重力异常场值转换到空间域重力异常场值，可以得到目标区域内任一点的空间域重力异常场值。
上述二度体重力异常积分解数值模拟方法中，通过构建包含目标区域的二度体模型，确定对应的空间域重力异常表达式，其中空间域重力异常表达式包含空间域异常体剩余密度值，通过对空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式，其中波数域异常体剩余密度值是对空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的，根据矩形沿x轴方向的尺寸信息，确定截止频率，根据截止频率确定波数域波数采样值范围，根据波数采样值范围和预先设置的采样点总数进行波数域非均匀采样，得到波数采样值，根据波数采样值，代入波数域重力异常表达式中，得到波数域重力异常场值，最后通过对波数域重力异常场值进行一维非均匀快速傅里叶反变换，得到目标区域内任一点的空间域重力异常场值。采用本方法进行二度体重力异常积分解数值模拟，波数域非均匀采样可以根据重力异常场的频谱，对应波数域采样的间隔，实现根据波数域异常场变化规律进行采样，在波数域异常场变化较快的地方多采样，在场变化缓慢的地方少采样，以最少的采样个数达到数值模拟精度要求，进一步提高计算效率，有效兼顾了重力异常数值模拟的计算精度与计算效率，解决了波数域重力异常数值模拟只适合规则测线、数值模拟方法不能同时兼顾计算效率和计算精度的问题。
在其中一个实施例中，空间域重力异常表达式为：
Figure 978024DEST_PATH_IMAGE048
式中：
Figure 766989DEST_PATH_IMAGE049
表示空间域重力异常场；
Figure 259150DEST_PATH_IMAGE050
表示万有引力常数；
Figure 969617DEST_PATH_IMAGE051
表示x方向剖分的矩形个数；
Figure 775899DEST_PATH_IMAGE052
表示z方向剖分的矩形个数；
Figure 747266DEST_PATH_IMAGE053
表示所述观测点坐标；
Figure 551274DEST_PATH_IMAGE054
表示编号为
Figure 608092DEST_PATH_IMAGE055
的矩形中心坐标；
Figure 218065DEST_PATH_IMAGE056
表示编号为
Figure 919304DEST_PATH_IMAGE057
的矩形的空间域异常体剩余密度值。
对本实施例中空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式为：
Figure 18847DEST_PATH_IMAGE058
式中：
Figure 438327DEST_PATH_IMAGE059
表示波数域重力异常场；k表示波数；
Figure 586412DEST_PATH_IMAGE060
表示波数域异常体剩余密度值；
Figure 266792DEST_PATH_IMAGE061
表示符号函数：
Figure 678182DEST_PATH_IMAGE062
其中，波数域异常体剩余密度值
Figure 381696DEST_PATH_IMAGE063
是对空间域异常体剩余密度值
Figure 661367DEST_PATH_IMAGE064
进行非均匀采样快速傅里叶变换得到的。
由于构建的二度体模型中，矩形尺寸是根据重力场的分布剖分的，不一定均匀剖分，本具体实施例中的波数域重力异常表达式可以实现在波数域进行非均匀采样，提高计算效率。
在其中一个实施例中，如图2所示，一种二度体重力异常积分解数值模拟方法主要包括四个部分：二度体模型表示、波数域重力异常公式推导、波数域非均匀采样和空间域重力异常计算。二度体模型构建包括设定目标区域及矩形几何尺寸，设定矩形个数和根据二度体形状对每个矩形进行赋值，以此刻画不同截面二度体；波数域重力异常公式推导是对空间域重力异常表达式进行一维傅里叶变换；波数域非均匀采样是根据重力异常频谱变换规律进行非均匀采样；空间域重力异常计算是对计算的非均匀采样频谱进行一维NUFFT反变换，得到空间域重力异常。
在其中一个实施例中，非均匀采样快速傅里叶变换为：
Figure 71620DEST_PATH_IMAGE065
式中：i为虚数单位，
Figure 778545DEST_PATH_IMAGE066
为给定离散采样点
Figure 969355DEST_PATH_IMAGE067
对应采样点的值，
Figure 396925DEST_PATH_IMAGE068
为计算的离散采样点
Figure 317476DEST_PATH_IMAGE069
傅里叶变换频谱，N表示采样点总数；
如图3所示，其中，实现非均匀快速傅里叶变换的步骤为：
步骤302：根据
Figure 601827DEST_PATH_IMAGE070
临近q个均匀采样点的傅里叶变换基，得到近似非均匀傅里叶变换基
Figure 483196DEST_PATH_IMAGE071
为：
Figure 839091DEST_PATH_IMAGE072
式中：m表示过采样因子，
Figure 489515DEST_PATH_IMAGE073
表示权重因子，
Figure 210346DEST_PATH_IMAGE074
为精度因子，
Figure 438065DEST_PATH_IMAGE075
表示对
Figure 207438DEST_PATH_IMAGE076
取整；
步骤304：根据采样点的值和权重因子，计算新傅里叶变换基下对应的傅里叶变换系数
Figure 102582DEST_PATH_IMAGE077
：
Figure 994314DEST_PATH_IMAGE078
步骤306：采用均匀的快速傅里叶变换，得到：
Figure 584696DEST_PATH_IMAGE079
式中，
Figure 305831DEST_PATH_IMAGE080
表示傅立叶变换之后的频谱。
快速傅里叶变换扩边法和Gauss-FFT在重力异常正演数值模拟中有着广泛应用，但这两种方法在重力异常正演数值模拟的应用中只能均匀采样，难以实现计算精度和计算效率的统一。而本实施例中的非均匀快速傅里叶变换是基于将局部插值与FFT相结合来实现非均匀采样的快速傅里叶变换，将本实施例中的非均匀快速傅里叶变换应用于重力异常正演数值模拟中，可以实现波数域的非均匀采样，进而在满足精度的前提下提高计算效率。
在其中一个实施例中，将目标区域分成若干个矩形，每个矩形的x轴方向尺寸
Figure 930847DEST_PATH_IMAGE081
、z轴方向尺寸
Figure 993481DEST_PATH_IMAGE082
可以相同，也可以不同；在重力异常场变化快的地方加密网格，在重力异常场变化慢的地方稀疏网格。采用这种方式可以在空间域中根据重力异常场的分布情况决定采样间隔，在重力异常场变化快的地方多采样，在重力异常场变化慢的地方少采样，从而提高计算效率。
在其中一个实施例中，根据矩形沿x轴方向的尺寸信息，确定截止频率为：
Figure 195792DEST_PATH_IMAGE083
式中，
Figure 306968DEST_PATH_IMAGE084
表示截止频率；
Figure 911124DEST_PATH_IMAGE085
表示x轴方向矩形尺寸的最小值。
根据截止频率确定波数域波数采样值范围为
Figure 144660DEST_PATH_IMAGE086
，再根据波数采样值范围得到波数采样值，包括：
在
Figure 444054DEST_PATH_IMAGE087
上波数依次取值为：
Figure 749133DEST_PATH_IMAGE088
在
Figure 879900DEST_PATH_IMAGE089
上波数依次取值为：
Figure 487599DEST_PATH_IMAGE090
式中，
Figure 398923DEST_PATH_IMAGE091
，N表示采样点总数。
本实施例中的波数采样方式是将波数范围
Figure 383060DEST_PATH_IMAGE092
在对数轴上划分为N个区间，每个区间中波数范围的对数距离为1，采用该方式可以实现在波数绝对值较小的区域密集采样，随着波数绝对值增大，采样逐渐稀疏。通过这种方式可以在满足精度的前提下减少采样点的数量，进而提高计算效率。
在一个具体实施例中，目标区域有一个截面为圆形的二度体如图4所示，目标区域范围为：x轴方向从-128km到128km，z轴方向从0km到100km，
Figure 368334DEST_PATH_IMAGE093
和
Figure 271567DEST_PATH_IMAGE094
均为1km，可知x轴方向剖分矩形个数为256个，z轴方向剖分矩形个数为100个，圆心坐标为(0km,80km)，圆的半径为5km，剩余密度为10kg/m3，计算测线为水平地面上的重力异常。
利用Fortran语言编程实现，运行程序所用的个人电脑配置为：CPU-Inter Corei7-8700，主频为3.2GHz，运行内存为8.00GB。
图5为空间域均匀剖分无限长圆柱体产生的重力异常波数域非均匀采样模拟结果与理论值对比图，其中，NUFFT表示非均匀快速傅里叶变换，N表示采样点总数。从图5中可以看出采用非均匀快速傅里叶变换选取150和256个波数域采样点总数，模拟结果与理论值吻合的很好，精度较高。图6表示本实施例中波数域非均匀采样模拟结果与理论值的相对误差，可以看出，波数域采样点总数为256时比波数域采样点总数为100时，模拟结果与理论值的相对误差更小，计算精度更高。
对于同样的二度体模型，统计了快速傅里叶变换扩边法（FFT扩边法）、高斯快速傅里叶变换法（Gauss-FFT方法）与非均匀快速傅里叶变换法（NUFFT方法）计算的相对均方根误差及计算时间，对比如表1，可以看出在达到相同精度的情况下，非均匀快速傅里叶变换法计算所用时间更短，计算效率更高。
Figure 811133DEST_PATH_IMAGE095
其中，4L、8L为快速傅里叶变换扩边法的扩边比，2、3分别为Gauss-FFT中采用高斯点的个数，rrms(%)表示相对均匀方根误差。
Figure 192436DEST_PATH_IMAGE096
Figure 501058DEST_PATH_IMAGE097
表示解析解，
Figure 247297DEST_PATH_IMAGE098
表示数值解。
在另一个具体实施例中，目标区域范围为：x轴方向从-128km到128km，z轴方向从0km到100km，在空间域进行非均匀采样，在(0km,15km)范围内，剖分的矩形
Figure 133213DEST_PATH_IMAGE099
为0.5km，在(15km,98km) 范围内，矩形
Figure 193573DEST_PATH_IMAGE099
为1km，在(98km,128km)范围内，矩形
Figure 153439DEST_PATH_IMAGE099
为2km，-128km到0km关于（0,0）对称。采用非均匀快速傅里叶变换选取150和256个波数域采样点总数，计算结果如图7所示，可以看出，模拟结果与理论值吻合的很好，精度较高。图8表示本实施例中波数域非均匀采样模拟结果与理论值的相对误差，可以看出，波数域采样点总数为256时比波数域采样点总数为150时，模拟结果与理论值的相对误差更小，计算精度更高。
应该理解的是，虽然图1的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，图1中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。
在一个实施例中，如图9所示，提供了一种二度体重力异常积分解数值模拟装置，包括：二度体模型构建模块902、空间域重力异常表达式确定模块904、波数域重力异常表达式确定模块906、波数采样值确定模块908、波数域重力异常场值确定模块910和空间域重力异常场值确定模块912，其中：
二度体模型构建模块902，用于根据观测点在观测点坐标中的位置信息，构建包含目标区域的二度体模型；观测点坐标包括x轴方向和z轴方向；二度体模型是将目标区域分成若干个矩形，矩形的截面形状用于确定每个矩形的空间域异常体剩余密度值；
空间域重力异常表达式确定模块904，用于根据二度体模型，确定对应的空间域重力异常表达式；空间域重力异常表达式包含空间域异常体剩余密度值；
波数域重力异常表达式确定模块906，用于对空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式；其中，波数域重力异常表达式包含波数域异常体剩余密度值；波数域异常体剩余密度值是对空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的；
波数采样值确定模块908，用于根据矩形沿x轴方向的尺寸信息，确定截止频率，根据截止频率确定波数域波数采样值范围，根据波数采样值范围和预先设置的采样点总数，在波数域进行非均匀采样，得到波数采样值；
波数域重力异常场值确定模块910，用于将波数采样值代入波数域重力异常表达式，得到波数域重力异常场值；
空间域重力异常场值确定模块912，用于通过对波数域重力异常场值进行一维非均匀快速傅里叶反变换，得到目标区域内任一点的空间域重力异常场值。
二度体模型构建模块902还包括将目标区域分成若干个矩形，每个矩形的x轴方向尺寸
Figure 664055DEST_PATH_IMAGE100
、z轴方向尺寸
Figure 647054DEST_PATH_IMAGE101
可以相同，也可以不同；在重力异常场变化快的地方加密网格，在重力异常场变化慢的地方稀疏网格。
波数采样值确定模块908还包括根据所述矩形沿x轴方向的尺寸信息，确定截止频率为：
Figure 635739DEST_PATH_IMAGE102
式中，
Figure 715690DEST_PATH_IMAGE103
表示截止频率；
Figure 6994DEST_PATH_IMAGE104
表示x轴方向矩形尺寸的最小值。
波数采样值确定模块908还包括根据截止频率确定波数域波数采样值范围为
Figure 601924DEST_PATH_IMAGE105
；
根据波数采样值范围得到波数采样值，包括：
在
Figure 66403DEST_PATH_IMAGE106
上波数依次取值为：
Figure 204123DEST_PATH_IMAGE107
在
Figure 525383DEST_PATH_IMAGE108
上波数依次取值为：
Figure 748554DEST_PATH_IMAGE109
式中，
Figure 751145DEST_PATH_IMAGE110
，N表示采样点总数。
关于二度体重力异常积分解数值模拟装置的具体限定可以参见上文中对于二度体重力异常积分解数值模拟方法的限定，在此不再赘述。上述二度体重力异常积分解数值模拟装置中的各个模块可全部或部分通过软件、硬件及其组合来实现。上述各模块可以硬件形式内嵌于或独立于计算机设备中的处理器中，也可以以软件形式存储于计算机设备中的存储器中，以便于处理器调用执行以上各个模块对应的操作。
在一个实施例中，提供了一种计算机设备，该计算机设备可以是终端，其内部结构图可以如图10所示。该计算机设备包括通过系统总线连接的处理器、存储器、网络接口、显示屏和输入装置。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现一种二度体重力异常积分解数值模拟方法。该计算机设备的显示屏可以是液晶显示屏或者电子墨水显示屏，该计算机设备的输入装置可以是显示屏上覆盖的触摸层，也可以是计算机设备外壳上设置的按键、轨迹球或触控板，还可以是外接的键盘、触控板或鼠标等。
本领域技术人员可以理解，图10中示出的结构，仅仅是与本申请方案相关的部分结构的框图，并不构成对本申请方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。
在一个实施例中，提供了一种计算机设备，包括存储器和处理器，该存储器存储有计算机程序，该处理器执行计算机程序时实现上述方法实施例中的步骤。
在一个实施例中，提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述方法实施例中的步骤。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器（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 (10)

1.一种二度体重力异常积分解数值模拟方法，其特征在于，所述方法包括：
根据观测点在笛卡尔坐标系下的位置信息，构建包含目标区域的二度体模型；所述观测点坐标包括x轴方向和z轴方向；所述二度体模型是将所述目标区域分成若干个矩形，所述二度体的截面形状用于确定每个矩形的空间域异常体剩余密度值；
根据所述二度体模型，确定对应的空间域重力异常表达式；所述空间域重力异常表达式包含空间域异常体剩余密度值；
对所述空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式；其中，所述波数域重力异常表达式包含波数域异常体剩余密度值；所述波数域异常体剩余密度值是对所述空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的；
根据所述矩形x轴方向的尺寸信息，确定截止频率，根据所述截止频率确定波数域波数采样值范围，根据所述波数采样值范围和预先设置的采样点总数，在波数域进行非均匀采样，得到波数采样值；
将所述波数采样值代入所述波数域重力异常表达式，得到波数域重力异常场值；
通过对所述波数域重力异常场值进行一维非均匀快速傅里叶反变换，得到目标区域内任一点的空间域重力异常场值。
2.根据权利要求1所述的方法，其特征在于，所述空间域重力异常表达式为：
Figure 380630DEST_PATH_IMAGE001
式中：
Figure 942061DEST_PATH_IMAGE002
表示空间域重力异常场；
Figure 241455DEST_PATH_IMAGE003
表示万有引力常数；
Figure 953059DEST_PATH_IMAGE004
表示x方向剖分的矩形个数；
Figure 700739DEST_PATH_IMAGE005
表示z方向剖分的矩形个数；
Figure 42859DEST_PATH_IMAGE006
表示所述观测点坐标；
Figure 891866DEST_PATH_IMAGE007
表示编号为
Figure 735057DEST_PATH_IMAGE008
的矩形中心坐标；
Figure 923593DEST_PATH_IMAGE009
表示编号为
Figure 498931DEST_PATH_IMAGE010
的矩形的空间域异常体剩余密度值，
Figure 897551DEST_PATH_IMAGE011
表示矩形的x轴方向尺寸，
Figure 419800DEST_PATH_IMAGE012
矩形的z轴方向尺寸。
3.根据权利要求2所述的方法，其特征在于，所述波数域重力异常表达式为：
Figure 525159DEST_PATH_IMAGE013
式中：
Figure 599294DEST_PATH_IMAGE014
表示波数域重力异常场；k表示波数；
Figure 94997DEST_PATH_IMAGE015
表示波数域异常体剩余密度值；
Figure 217674DEST_PATH_IMAGE016
表示符号函数：
Figure 505436DEST_PATH_IMAGE017
。
4.根据权利要求1所述的方法，其特征在于，所述非均匀采样快速傅里叶变换为：
Figure 360260DEST_PATH_IMAGE018
式中：i为虚数单位，
Figure 671155DEST_PATH_IMAGE019
为给定离散采样点
Figure 659840DEST_PATH_IMAGE020
对应采样点的值，
Figure 411895DEST_PATH_IMAGE021
为计算的离散采样点
Figure 765516DEST_PATH_IMAGE022
傅里叶变换频谱，N表示采样点总数；
其中，所述非均匀快速傅里叶变换的实现步骤为：
根据
Figure 626025DEST_PATH_IMAGE023
临近q个均匀采样点的傅里叶变换基，得到近似非均匀傅里叶变换基
Figure 762608DEST_PATH_IMAGE024
为：
Figure 697066DEST_PATH_IMAGE025
式中：m表示过采样因子，
Figure 283905DEST_PATH_IMAGE026
表示权重因子，
Figure 241497DEST_PATH_IMAGE027
为精度因子，
Figure 509667DEST_PATH_IMAGE028
表示对
Figure 626528DEST_PATH_IMAGE029
取整；
根据采样点的值和权重因子，计算新傅里叶变换基对应的傅里叶变换系数
Figure 994055DEST_PATH_IMAGE030
：
Figure 501260DEST_PATH_IMAGE031
采用均匀的快速傅里叶变换，得到：
Figure 369859DEST_PATH_IMAGE032
式中，
Figure 216592DEST_PATH_IMAGE033
表示傅立叶变换之后的频谱。
5.根据权利要求1所述的方法，其特征在于，将所述目标区域分成若干个矩形包括：
将所述目标区域分成若干个矩形，每个矩形的x轴方向尺寸
Figure 82917DEST_PATH_IMAGE034
、z轴方向尺寸
Figure 139735DEST_PATH_IMAGE035
可以相同，也可以不同；在重力异常场变化快的地方加密网格，在重力异常场变化慢的地方稀疏网格。
6.根据权利要求1所述的方法，其特征在于，根据所述矩形沿x轴方向的尺寸信息，确定截止频率，包括：
根据所述矩形沿x轴方向的尺寸信息，确定截止频率为：
Figure 687391DEST_PATH_IMAGE036
式中，
Figure 778843DEST_PATH_IMAGE037
表示截止频率；
Figure 550490DEST_PATH_IMAGE038
表示x轴方向矩形尺寸的最小值。
7.根据权利要求6所述的方法，其特征在于，根据所述截止频率确定波数域波数采样值范围，根据所述波数采样值范围和预先设置的采样点总数，在波数域进行非均匀采样，得到波数采样值，包括：
根据所述截止频率确定波数域波数采样值范围为
Figure 235550DEST_PATH_IMAGE039
；
根据所述波数采样值范围得到波数采样值，包括：
在
Figure 445951DEST_PATH_IMAGE040
上波数依次取值为：
Figure 64014DEST_PATH_IMAGE041
在
Figure 209825DEST_PATH_IMAGE042
上波数依次取值为：
Figure 975655DEST_PATH_IMAGE043
式中，
Figure 130693DEST_PATH_IMAGE044
，N表示采样点总数。
8.一种二度体重力异常积分解数值模拟装置，其特征在于，所述装置包括：
二度体模型构建模块，根据观测点在笛卡尔坐标系下的位置信息，构建包含目标区域的二度体模型；所述观测点坐标包括x轴方向和z轴方向；所述二度体模型是将所述目标区域分成若干个矩形，所述二度体的截面形状用于确定每个矩形的空间域异常体剩余密度值；
空间域重力异常表达式确定模块，用于根据所述二度体模型，确定对应的空间域重力异常表达式；所述空间域重力异常表达式包含空间域异常体剩余密度值；
波数域重力异常表达式确定模块，用于对所述空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式；其中，所述波数域重力异常表达式包含波数域异常体剩余密度值；所述波数域异常体剩余密度值是对所述空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的；
波数采样值确定模块，用于根据所述矩形沿x轴方向的尺寸信息，确定截止频率，根据所述截止频率确定波数域波数采样值范围，根据所述波数采样值范围和预先设置的采样点总数，在波数域进行非均匀采样，得到波数采样值；
波数域重力异常场值确定模块，用于将所述波数采样值代入所述波数域重力异常表达式，得到波数域重力异常场值；
空间域重力异常场值确定模块，用于通过对所述波数域重力异常场值进行一维非均匀快速傅里叶反变换，得到目标区域内任一点的空间域重力异常场值。
9.一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至7中任一项所述方法的步骤。
10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1至7中任一项所述的方法的步骤。
CN202011127944.7A 2020-10-21 2020-10-21 二度体重力异常积分解数值模拟方法和装置 Active CN111967169B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011127944.7A CN111967169B (zh) 2020-10-21 2020-10-21 二度体重力异常积分解数值模拟方法和装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011127944.7A CN111967169B (zh) 2020-10-21 2020-10-21 二度体重力异常积分解数值模拟方法和装置

Publications (2)

Publication Number Publication Date
CN111967169A true CN111967169A (zh) 2020-11-20
CN111967169B CN111967169B (zh) 2021-01-01

Family

ID=73387104

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011127944.7A Active CN111967169B (zh) 2020-10-21 2020-10-21 二度体重力异常积分解数值模拟方法和装置

Country Status (1)

Country Link
CN (1) CN111967169B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113076678A (zh) * 2021-04-15 2021-07-06 中南大学 一种频率域二度体重力异常快速数值模拟方法和装置
CN115659579A (zh) * 2022-09-05 2023-01-31 中南大学 基于3d as-ft的三维重力场数值模拟方法及系统
CN116244877A (zh) * 2022-09-05 2023-06-09 中南大学 基于3d as-ft的三维磁场数值模拟方法及系统

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2004090575A1 (en) * 2003-04-09 2004-10-21 Norsar Method for simulating local prestack depth migrated seismic images
US20110082666A1 (en) * 2009-10-02 2011-04-07 Wen-Chih Chen Numerical Engineering Model System And Accessing Method Thereof
CN106855904A (zh) * 2017-01-10 2017-06-16 桂林理工大学 一种二度体重力异常计算方法
CN107402409A (zh) * 2017-09-26 2017-11-28 西南石油大学 一种三维不规则地层起伏界面重力正演方法
CN108197389A (zh) * 2018-01-04 2018-06-22 中南大学 二维强磁性体磁场的快速、高精度数值模拟方法
CN110428498A (zh) * 2019-07-12 2019-11-08 南京理工大学 一种将三维地质模型转化为数值计算模型的方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2004090575A1 (en) * 2003-04-09 2004-10-21 Norsar Method for simulating local prestack depth migrated seismic images
US20110082666A1 (en) * 2009-10-02 2011-04-07 Wen-Chih Chen Numerical Engineering Model System And Accessing Method Thereof
CN106855904A (zh) * 2017-01-10 2017-06-16 桂林理工大学 一种二度体重力异常计算方法
CN107402409A (zh) * 2017-09-26 2017-11-28 西南石油大学 一种三维不规则地层起伏界面重力正演方法
CN108197389A (zh) * 2018-01-04 2018-06-22 中南大学 二维强磁性体磁场的快速、高精度数值模拟方法
CN110428498A (zh) * 2019-07-12 2019-11-08 南京理工大学 一种将三维地质模型转化为数值计算模型的方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
DAI SHI-KUN: "three-dimensional numerical modeling of gravity anomalies based on poisson equation in space-wavenumber mixed domain", 《APPLIED GEOPHYSICS》 *
李颖梅: "复杂条件二维重力场及重力张量场空间波数域正演方法", 《物探化探计算技术》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113076678A (zh) * 2021-04-15 2021-07-06 中南大学 一种频率域二度体重力异常快速数值模拟方法和装置
CN113076678B (zh) * 2021-04-15 2022-04-19 中南大学 一种频率域二度体重力异常快速数值模拟方法和装置
CN115659579A (zh) * 2022-09-05 2023-01-31 中南大学 基于3d as-ft的三维重力场数值模拟方法及系统
CN116244877A (zh) * 2022-09-05 2023-06-09 中南大学 基于3d as-ft的三维磁场数值模拟方法及系统
CN116244877B (zh) * 2022-09-05 2023-11-14 中南大学 基于3d傅里叶变换的三维磁场数值模拟方法及系统

Also Published As

Publication number Publication date
CN111967169B (zh) 2021-01-01

Similar Documents

Publication Publication Date Title
CN111967169B (zh) 二度体重力异常积分解数值模拟方法和装置
CN112287534B (zh) 基于nufft的二维磁异常快速正演模拟方法和装置
CN112800657B (zh) 基于复杂地形的重力场数值模拟方法、装置和计算机设备
Manolis et al. Elastic waves in continuous and discontinuous geological media by boundary integral equation methods: A review
Minniakhmetov et al. High-order spatial simulation using Legendre-like orthogonal splines
Chen et al. A robust method of thin plate spline and its application to DEM construction
CN113420487B (zh) 一种三维重力梯度张量数值模拟方法、装置、设备和介质
CN109459787B (zh) 基于地震槽波全波形反演的煤矿井下构造成像方法及系统
CN114021408A (zh) 二维强磁场数值模拟方法、装置、设备及介质
CN114065511A (zh) 起伏地形下大地电磁二维正演数值模拟方法、装置、设备及介质
Godin Wentzel–Kramers–Brillouin approximation for atmospheric waves
CN109490978B (zh) 一种起伏地层的频率域快速高精度正演方法
US20150134308A1 (en) Method and device for acquiring optimization coefficient, and related method and device for simulating wave field
CN112666612A (zh) 基于禁忌搜索的大地电磁二维反演方法
CN113221409B (zh) 一种有限元和边界元耦合的声波二维数值模拟方法和装置
CN113076678B (zh) 一种频率域二度体重力异常快速数值模拟方法和装置
CN113268702B (zh) 一种频率域磁梯度张量变换方法、装置和计算机设备
CN113673163B (zh) 一种三维磁异常数快速正演方法、装置和计算机设备
CN111158059A (zh) 基于三次b样条函数的重力反演方法
Sun et al. Probabilistic analysis of width‐limited 3D slope in spatially variable soils: UBLA enhanced with efficiency‐improved discretization of horn‐like failure mechanism
Yuen et al. Geophysical applications of multidimensional filtering with wavelets
Bucha et al. Cap integration in spectral gravity forward modelling up to the full gravity tensor
Desmars et al. Reconstruction of ocean surfaces from randomly distributed measurements using a grid-based method
Shargatov Dynamics and stability of air bubbles in a porous medium
CN114970289A (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