CN113779818A - 三维地质体其电磁场数值模拟方法、装置、设备及介质 - Google Patents

三维地质体其电磁场数值模拟方法、装置、设备及介质 Download PDF

Info

Publication number
CN113779818A
CN113779818A CN202111347154.4A CN202111347154A CN113779818A CN 113779818 A CN113779818 A CN 113779818A CN 202111347154 A CN202111347154 A CN 202111347154A CN 113779818 A CN113779818 A CN 113779818A
Authority
CN
China
Prior art keywords
electric field
dimensional
sub
matrix
interpolation
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
CN202111347154.4A
Other languages
English (en)
Other versions
CN113779818B (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 CN202111347154.4A priority Critical patent/CN113779818B/zh
Publication of CN113779818A publication Critical patent/CN113779818A/zh
Application granted granted Critical
Publication of CN113779818B publication Critical patent/CN113779818B/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
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

三维地质体其电磁场数值模拟方法、装置、设备及介质，包括：构建内部包含勘探目标的三维长方体模型；对三维长方体模型进行网格剖分，剖分成若干个小长方体单元，给每个小长方体单元的电阻率赋值，得到三维异常体模型；根据频率参数和三维异常体模型，构建对应的电场控制方程，并计算电场控制方程的系数矩阵和右端项；将三维异常体模型沿Z方向划分为多层子区域，在各层子区域分别构建对应的插值算子，并施加对应的插值算子给各层子区域；将各层子区域对应的插值算子合成后得到一个总的稀疏采样算子，并计算新的系数矩阵和右端项，由新的系数矩阵和右端项计算得到电磁场。该方法能够满足大规模电磁数据精细、快速正演成像的需求。

Description

三维地质体其电磁场数值模拟方法、装置、设备及介质
技术领域
本发明属于强磁体数值模拟技术领域，特别涉及一种三维地质体其电磁场数值模拟方法、装置、设备及介质。
背景技术
可控源电磁法是地球物理勘探中一种重要的勘探方法，具有勘探深度大、分辨率高、高阻屏蔽小、工作效率高等特点，目前已广泛应用于金属矿、油气、地热资源勘探等各个方面。其中可控源电磁正演技术对于可控源电磁勘探中的测线布置、可行性研究和数据解释等方面发挥着至关重要的作用，可有效地解释复杂的深部地质情况。目前，业界有大量的研究致力于发展电磁正演技术，其中有限差分模拟技术利用差分方程近似进行求解，理论简单且容易实现，在电磁正演中被广泛应用。
现阶段，由于模型精细离散化后形成的计算点数量的增加，正演所需的内存和时间急剧增长，在普通计算机上难以计算，已经成为三维电磁快速成像的难点问题。而数据解释具有实时性，因此正演模拟方法的计算效率直接影响反演解释的效果。
目前，八叉树网格算法由于其高精度，是国际上的主流算法。已有人基于八叉树思想提出一种多分辨率网格算法，仅对水平方向进行细化，通过混合粗-细网格代替原来的细网格，极大地提高了被动源电磁法正演的计算效率，但是不同粗糙度网格边界的处理是一项很困难的工作。为此，又有人提出了一种被动源电磁分级下采样算法，该算法直接减少离散电磁场采样点的数量，不改变网格单元的数量。但是分级下采样算法受频率影响很大，频率较低时计算效率没有提升。
加之，新世纪以来，电磁测量技术的发展大大促进了地球电磁场的研究，大尺度大规模的电磁测量技术已经蓬勃发展，传统的正演模拟算法难以满足大量高分辨率数据精细快速反演的要求，因此研究一种基于多分辨率采样的高效高精度数值模拟方法具有现实意义。
发明内容
针对目前基于可控源电磁的三维正演方法大多存在计算量大和计算效率低的问题，本发明旨在提出一种三维地质体其电磁场数值模拟方法、装置、设备及介质，以满足大规模电磁数据精细、快速正演成像的需求。
为实现上述技术目的，本发明提出的技术方案为：
一方面，本发明提供一种三维地质体其电磁场数值模拟方法，包括：
勘探目标为三维异常体，构建内部包含所述勘探目标的三维长方体模型；
对所述三维长方体模型沿x、y、z方向进行网格剖分，剖分成若干个小长方体单元，得到三维长方体模型的网格剖分参数，根据所述勘探目标的电阻率分布，给每个小长方体单元的电阻率赋值，每一个小长方体单元的电阻率为常值，不同小长方体单元的电阻率值不同，得到刻画任意电阻率分布的三维异常体模型；
根据频率参数和三维异常体模型，构建其对应的一次电场和二次电场间的电场控制方程AE a =b，其中E a 代表二次电场，A为电场控制方程的系数矩阵，
Figure 78186DEST_PATH_IMAGE001
，
Figure 427259DEST_PATH_IMAGE002
表示旋度算子，
Figure 396352DEST_PATH_IMAGE003
表示双旋度算子，
Figure 223493DEST_PATH_IMAGE004
，ω表示角频率，通过ω=2πf求取，f表示给定的频率，μ表示磁导率，其值为4π×10-7，ρ表示每个小长方体单元的电阻率，
Figure 227221DEST_PATH_IMAGE005
，
Figure 481616DEST_PATH_IMAGE006
表示异常体电阻率，E b 表示一次电场；
从三维异常体模型中每个小长方体单元的电阻率中分离出异常体电阻率，得到背景电阻率；根据三维异常体模型中每个小长方体单元的背景电阻率、所述三维长方体模型的网格剖分参数和频率参数计算三维异常体模型的一次电场，进而得到其对应的电场控制方程的右端项b；
将所述三维异常体模型沿Z方向划分为多层子区域，对各层子区域分别沿水平方向构建不同稀疏度的插值算子；
将各层子区域对应的插值算子合成后得到一个总的稀疏采样算子，电场控制方程的系数矩阵A和右端项b分别乘上总的稀疏采样算子得到新的稀疏矩阵A a,new 和修正后的右端项b a,new ，将新的稀疏矩阵A a,new 和修正后的右端项b a,new 带入电场控制方程，求解得到修正后的二次电场E a,new ；
对所述修正后的二次电场和一次电场求和得到总的电场E，由总的电场E即可计算得到磁场
Figure 938006DEST_PATH_IMAGE007
。
三维长方体模型的网格剖分参数包括三维长方体模型x、y、z方向网格剖分得到的小长方体的数量N x 、N y 、N z ，各个小长方体x、y、z方向的棱边长度，以及对各小长方体单元进行编号后得到的各小长方体单元的编号和坐标位置。
另一方面，本发明提供一种三维地质体其电磁场数值模拟装置，包括：
第一模块，用于构建内部包含所述勘探目标的三维长方体模型，所述勘探目标为三维异常体；
第二模块，用于对所述三维长方体模型沿x、y、z方向进行网格剖分，剖分成若干个小长方体单元，得到三维长方体模型的网格剖分参数，根据异常体的电阻率分布，给每个小长方体单元的电阻率赋值，每一个小长方体单元的电阻率为常值，不同小长方体单元的电阻率值不同，得到刻画任意电阻率分布的三维异常体模型；
第三模块，用于根据频率参数和三维异常体模型，构建其对应的一次电场和二次电场间的电场控制方程AE a =b，其中E a 代表二次电场，A为电场控制方程的系数矩阵，
Figure 568838DEST_PATH_IMAGE008
，
Figure 427073DEST_PATH_IMAGE002
表示旋度算子，
Figure 115018DEST_PATH_IMAGE009
表示双旋度算子，
Figure 58704DEST_PATH_IMAGE010
，ω表示角频率，通过ω=2πf求取，f表示预先给定的频率，μ表示磁导率，其值为4π×10-7，ρ表示每个小长方体单元的电阻率，
Figure 227648DEST_PATH_IMAGE011
，
Figure 674810DEST_PATH_IMAGE012
表示异常体电阻率，E b 表示一次电场；
第四模块，用于从三维异常体模型中每个小长方体单元的电阻率中分离出异常体电阻率，得到背景电阻率；根据三维异常体模型中每个小长方体单元的背景电阻率、所述三维长方体模型的网格剖分参数和频率参数计算三维异常体模型的一次电场，进而得到其对应的电场控制方程的右端项b；
第五模块，将所述三维异常体模型沿Z方向划分为多层子区域，对各层子区域分别沿水平方向构建不同稀疏度的插值算子；
第六模块，用于将各层子区域对应的插值算子合成后得到一个总的稀疏采样算子，电场控制方程的系数矩阵A和右端项b分别乘上总的稀疏采样算子得到新的稀疏矩阵A a,new 和修正后的右端项b a,new ，将新的稀疏矩阵A a,new 和修正后的右端项b a,new 带入电场控制方程，求解得到修正后的二次电场E a,new ；
第七模块，用于对所述修正后的二次电场和一次电场求和得到总的电场E，由总的电场E即可计算得到磁场
Figure 333324DEST_PATH_IMAGE013
。
另一方面，本发明提供一种计算机设备，包括存储器和处理器，存储器存储有计算机程序，处理器执行计算机程序时实现所述三维地质体其电磁场数值模拟方法中的步骤。
再一方面，本发明还提供一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现所述三维地质体其电磁场数值模拟方法中的步骤。
与现有技术相比，本发明的优点在于：
本发明所述主动源三维地质体其电磁场数值模拟方法，可以根据实际勘探任务，将整个研究区域灵活的分为多个子区域，在每个子区域水平方向施加不同“稀疏度”的插值算子，z方向计算点保持不变，通过分区域稀疏采样直接减少了离散电磁场的数量。因此，本发明大大提高了电磁正演的计算效率、减少了计算所需的内存；另外该方法不改变原有模型电阻率单元的数量，转变为降低离散计算点数量，与已有主流的方法相比，这更容易用于反演算法。
附图说明
图1是本发明一实施例中流程图；
图2是本发明一实施例中海洋油储体模型示意图；
图3为本发明一实施例得到的电场的参考解、数值解以及相对误差图，其中(a)为电场的实部的参考解，(b)为电场的实部的数值解，(c)为电场的实部的相对误差；（d)为电场的虚部的参考解，(e)为电场的虚部的数值解，(f)为电场的虚部的相对误差；
图4为一实施例中采用本发明方法相对传统方法降低的计算时间百分比的示意图；
图5为本发明一实施例的结构示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚明白，下面将以附图及详细叙述来清楚说明本发明所揭示内容的精神，任何所属技术领域技术人员在了解本发明内容的实施例后，当可由本发明内容所教示的技术，加以改变及修饰，其并不脱离本发明内容的精神与范围。本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。
参照图1，本发明一实施例中，提供一种三维地质体其电磁场数值模拟方法,包括：
(S1) 勘探目标为三维异常体，构建内部包含所述勘探目标的三维长方体模型；
(S2) 对所述三维长方体模型沿x、y、z方向进行网格剖分，剖分成若干个小长方体单元，得到三维长方体模型的网格剖分参数，根据所述勘探目标的电阻率分布，给每个小长方体单元的电阻率赋值，每一个小长方体单元的电阻率为常值，不同小长方体单元的电阻率值不同，得到刻画任意电阻率分布的三维异常体模型。需要说明的是，本发明中的网格剖分方式可以是均匀剖分，使得剖分成的各小长方体单元的尺寸相同，也可以是非均匀剖分，剖分成的各长方体单元可以具有不同的长宽高尺寸。
(S3) 根据频率参数和三维异常体模型，构建其对应的一次电场和二次电场间的电场控制方程AE a =b，其中E a 代表二次电场，A为电场控制方程的系数矩阵，
Figure 374093DEST_PATH_IMAGE014
，
Figure 471362DEST_PATH_IMAGE015
表示旋度算子，
Figure 913975DEST_PATH_IMAGE016
表示双旋度算子，
Figure 71287DEST_PATH_IMAGE017
，
Figure 661669DEST_PATH_IMAGE018
表示角频率，通过
Figure 906836DEST_PATH_IMAGE019
求取，f表示给定的频率，μ表示磁导率，其值为4π×10-7，ρ表示每个小长方体单元的电阻率，
Figure 594170DEST_PATH_IMAGE020
，
Figure 797749DEST_PATH_IMAGE021
表示异常体电阻率，E b 表示一次电场；
(S4) 从三维异常体模型中每个小长方体单元的电阻率ρ中分离出异常体电阻率
Figure 203322DEST_PATH_IMAGE021
，得到背景电阻率ρ b ；根据三维异常体模型中每个小长方体单元的背景电阻率ρ b 、所述三维长方体模型的网格剖分参数和频率参数f计算三维异常体模型的一次电场E b ，进而得到电场控制方程的右端项
Figure 520690DEST_PATH_IMAGE020
；
(S5) 将所述三维异常体模型沿Z方向划分为多层子区域，对各层子区域分别沿水平方向构建不同稀疏度的插值算子；
将各层子区域对应的插值算子合成后得到一个总的稀疏采样算子，电场控制方程的系数矩阵A和右端项b分别乘上总的稀疏采样算子得到新的稀疏矩阵A a,new 和修正后的右端项b a,new ，将新的稀疏矩阵A a,new 和修正后的右端项b a,new ，带入电场控制方程，求解得到修正后的二次电场E a,new ；
(S6) 对所述修正后的二次电场和一次电场求和得到总的电场E，由总的电场E即可计算得到磁场
Figure 62530DEST_PATH_IMAGE022
。
本发明的步骤(S1)中勘探目标为三维异常体，三维异常体的形状、大小、电阻率分布均不限，可以是任意形状、任意大小、任意电阻率分布的强磁性介质。
本发明的步骤(S2)中对所述三维长方体模型沿x、y、z方向进行网格剖分，其中沿x、y、z方向可以进行等间距的均匀剖分，也可以进行非等间距的剖分，具体的剖分方式不限，然后获取网格剖分参数，x、y、z方向剖分的小长方体单元个数分别用N x 、N y 、N z 表示，各个小长方体单元x、y、z方向的棱边长度（即每个小长方体单元的长宽高），如果采用均匀剖分，那么各个小长方体单元x、y、z方向的棱边长度分别为
Figure 233748DEST_PATH_IMAGE023
，最后将各小长方体单元进行编号，获取各小长方体单元的编号以及坐标位置等网格剖分参数。
本发明的步骤(S2)中，根据三维异常体的电阻率分布，给每个小长方体单元的电阻率赋值，每一个小长方体单元的电阻率均为常值，不同小长方体单元的电阻率值可以不同，以此刻画任意电阻率分布的三维异常体模型，其中所述空气部分中的各小长方体单元的电阻率为1010欧姆米，用来模拟电磁场响应。
本发明一实施例中的步骤（S3)中，电场控制方程的系数矩阵A通过以下方法获得：
三维异常体模型中各个小长方体单元x、y、z方向的棱边长度，构成长度元矩阵
Figure 470825DEST_PATH_IMAGE024
；
三维异常体模型中各个小长方体单元的表面积，构成面积元矩阵S；
三维异常体模型中各个小长方体单元的体积，构成体积元矩阵V；
双旋度算子
Figure 713588DEST_PATH_IMAGE025
，其中O为拓扑矩阵，由-1和1两种元素组成。拓扑矩阵O作用于由长度元与场值乘积组成的矩阵，实现积分的求和运算，几何映射关系为长度元L到面积元S，每行对应一个面积元。
最后得到电场控制方程的系数矩阵
Figure 250880DEST_PATH_IMAGE026
。
本发明一实施例中的步骤（S4)中，一次电场E b 的三个分量分别是E r 、E φ 、E z ，通过下式求解：
Figure 920895DEST_PATH_IMAGE027
其中：A 0 是矢量位，A r 、A φ 和A z 为A 0 的三个分量；ρ b表示背景电场率模型；r表示网格剖分参数产生的极矩。
Figure 645269DEST_PATH_IMAGE028
表示对矢量位的散度。对于各个小长方体单元，
Figure 426143DEST_PATH_IMAGE029
，
Figure 614679DEST_PATH_IMAGE030
,
Figure 330962DEST_PATH_IMAGE031
分别表示各个小长方体单元x、y方向的棱边长度。对于各个小长方体单元，φ表示接收点与小长方体单元中心点连线在xy平面投影与x轴的夹角，
Figure 932845DEST_PATH_IMAGE032
表示夹角φ的变化量，
Figure 392776DEST_PATH_IMAGE033
表示z方向棱边长度的变化量。
本发明一实施例中的步骤(S5)中，将所述三维异常体模型沿Z方向划分为五层子区域，由上至下分别为高空层子区域、空气-地表层子区域、浅地层子区域、中地层子区域和深地层子区域。可以理解，上述实施例中将三维异常体模型沿Z方向划分为五层子区域仅仅是本发明的一个优选实施例，在实际应用中本领域技术人员可以将三维异常体模型沿Z方向划分多层子区域，包括但不限于五层。
上述实施例中将三维异常体模型沿Z方向划分为五层子区域，接下来分别为5层子区域定义不同的“稀疏度”，并为每层子区域建立插值算子。根据多分辨率采样的思想，减少的场可以通过相邻的场进行插值，因此需要对电场控制方程的系数矩阵进行调整。由于多分辨率采样仅考虑在水平方向上进行，而垂直方向保持不变，因而构建的插值算子只是对水平方向做插值处理。
对于高空层子区域和浅地层子区域，这两个区域对应的插值算子为1阶稀疏采样算子，可表示为
Figure 763715DEST_PATH_IMAGE034
，s 1为一阶插值矩阵。
高空层子区域和浅地层子区域采用相同的方法分别构建一阶插值矩阵s 1。构建一阶插值矩阵s 1的方法是：先构建一个一阶初始矩阵，该矩阵大小为(n,n)，主对角线的元素为1；
对子区域（高空层子区域或浅地层子区域）中的各小长方体单元，构建x和y方向一阶插值算子系数为：
Figure 382390DEST_PATH_IMAGE035
Figure 205990DEST_PATH_IMAGE036
其中
Figure 469612DEST_PATH_IMAGE037
、
Figure 695057DEST_PATH_IMAGE038
分别表示子区域中x方向上第2j-1、2j+1个小长方体单元上计算点对应的插值算子系数，j用于计数，j为正整数，j=1,2,3...；
Figure 753143DEST_PATH_IMAGE039
、
Figure 595197DEST_PATH_IMAGE040
分别表示子区域中y方向第2j-1、2j+1个小长方体单元计算点对应的插值算子系数，而
Figure 131352DEST_PATH_IMAGE041
和
Figure 211303DEST_PATH_IMAGE042
分别代表第2j、2j-1、2j+1个小长方体单元沿x和y方向的棱边长度。
将所构建的一阶插值算子系数分别对应施加到子区域中的相应小长方体单元上，有
Figure 705870DEST_PATH_IMAGE043
Figure 769641DEST_PATH_IMAGE044
Figure 109486DEST_PATH_IMAGE045
和
Figure 309523DEST_PATH_IMAGE046
表示沿x方向上第2j-1、2j+1个小长方体单元上计算点对应的二次电场，该计算点仍然会保留，记录下其位置（2j-1，2j）和（2j+1，2j），然后在一阶初始矩阵里面找到对应的位置并分别赋值插值系数
Figure 974991DEST_PATH_IMAGE047
、
Figure 994900DEST_PATH_IMAGE048
；
Figure 404015DEST_PATH_IMAGE049
和
Figure 458559DEST_PATH_IMAGE050
表示沿y方向上第2j-1、2j+1个小长方体单元上计算点对应的二次电场，该计算点仍然保留，记录下其位置（2j，2j-1）和（2j，2j+1）然后在一阶初始矩阵里面找到对应的位置并分别赋值插值系数
Figure 20560DEST_PATH_IMAGE051
、
Figure 793344DEST_PATH_IMAGE052
；
Figure 740571DEST_PATH_IMAGE053
和
Figure 384042DEST_PATH_IMAGE054
表示x方向、y方向第2j个小长方体单元上计算点被插值的二次电场，该计算点将会被去除，记录下其位置（2j，2j），然后在一阶初始矩阵里面找到对应的位置并赋值0，以此规律，直至一阶初始矩阵中所有位置均赋值完毕，得到最终的一阶插值矩阵s 1。
对于空气-地表层子区域其插值算子为单位矩阵I，其矩阵大小为(n,n)，主对角线的元素为1。
对于中地层子区域，其插值算子为2阶稀疏采样算子，可表示为
Figure 656891DEST_PATH_IMAGE055
，s 2为二阶插值矩阵。
构建二阶插值矩阵s 2，方法是：先构建一个二阶初始矩阵，其矩阵大小为(n,n)，主对角线的元素为1；
对中地层子区域中的各小长方体单元，构建x和y方向二阶插值系数为：
Figure 651392DEST_PATH_IMAGE056
Figure 402311DEST_PATH_IMAGE057
其中
Figure 431446DEST_PATH_IMAGE058
、
Figure 344039DEST_PATH_IMAGE059
分别表示中地层子区域中x方向上第4j-3、4j+1个小长方体单元上计算点对应的插值系数；
Figure 91415DEST_PATH_IMAGE060
、
Figure 380445DEST_PATH_IMAGE061
分别表示中地层子区域中y方向第4j-3、4j+1个小长方体单元计算点对应的插值系数，
Figure 264087DEST_PATH_IMAGE062
和
Figure 347581DEST_PATH_IMAGE063
分别代表第4j-3、4j-2、4j-1、4j个小长方体单元沿x和y方向的棱边长度。
将所构建的二阶插值算子分别对应施加到中地层子区域中的相应小长方体单元上，有：
Figure 316674DEST_PATH_IMAGE064
Figure 409395DEST_PATH_IMAGE065
Figure 147544DEST_PATH_IMAGE066
和
Figure 664588DEST_PATH_IMAGE067
表示中地层子区域沿x方向上第4j-3、4j+1个小长方体单元上计算点对应的二次电场，该计算点仍然会保留，记录下其位置（4j-1，4j-3）和（4j-1，4j+1），然后在二阶初始矩阵里面找到对应的位置并分别赋值插值系数
Figure 120977DEST_PATH_IMAGE068
、
Figure 751810DEST_PATH_IMAGE069
；
Figure 610045DEST_PATH_IMAGE070
和
Figure 300920DEST_PATH_IMAGE071
表示沿y方向上第4j-3、4j+1个小长方体单元上计算点对应的电场，该计算点仍然保留记录下其位置（4j-3，4j-1）和（4j+1，4j-1），然后在二阶初始矩阵里面找到对应的位置并分别赋值插值系数
Figure 979026DEST_PATH_IMAGE072
、
Figure 413550DEST_PATH_IMAGE073
；
Figure 126291DEST_PATH_IMAGE074
和
Figure 988068DEST_PATH_IMAGE075
表示x方向、y方向第4j-1个小长方体单元上计算点被插值的电场，该计算点将会被去除，记录下其位置（4j-1，4j-1），然后在二阶初始矩阵里面找到对应的位置并赋值0，以此规律，直至二阶初始矩阵中所有位置均赋值完毕，得到最终的二阶插值矩阵s 2。
对于深地层子区域，其插值算子为3阶稀疏采样算子，可表示为
Figure 419049DEST_PATH_IMAGE076
，s 3为二阶插值矩阵。
构建三阶插值矩阵s 3，其方法是：先构建一个三阶初始矩阵，其矩阵大小为(n,n)，主对角线的元素为1；
对深地层子区域中的各小长方体单元，然后构建x和y方向三阶插值系数为：
Figure 391684DEST_PATH_IMAGE077
Figure 958932DEST_PATH_IMAGE078
其中
Figure 991610DEST_PATH_IMAGE079
、
Figure 644308DEST_PATH_IMAGE080
分别表示深地层子区域中x方向上第8j-7、8j+1个小长方体单元上计算点对应的插值系数；
Figure 420634DEST_PATH_IMAGE081
、
Figure 107967DEST_PATH_IMAGE082
分别表示深地层子区域中y方向第8j-7、8j+1个小长方体单元计算点对应的插值系数，
Figure 314476DEST_PATH_IMAGE083
和
Figure 454471DEST_PATH_IMAGE084
分别代表第8j-7、8j-6、8j-5、8j-4、8j-3、8j-2、8j-1、8j个小长方体单元沿x和y方向的棱边长度。
将所构建的三阶插值算子分别对应施加到深地层子区域中的相应小长方体单元上，有：
Figure 768909DEST_PATH_IMAGE085
Figure 310748DEST_PATH_IMAGE086
Figure 685229DEST_PATH_IMAGE087
和
Figure 312519DEST_PATH_IMAGE088
表示深地层子区域沿x方向上第8j-7、8j-1个小长方体单元上计算点对应的二次电场，该计算点仍然会保留，记录下其位置（8j-7，8j-3）和（8j-1，8j-3），然后在三阶初始矩阵里面找到对应的位置并分别赋值插值系数
Figure 430648DEST_PATH_IMAGE089
和
Figure 92574DEST_PATH_IMAGE090
；
Figure 637956DEST_PATH_IMAGE091
和
Figure 486963DEST_PATH_IMAGE092
表示沿y方向上第8j-7、8j+1个小长方体单元上计算点对应的电场，该计算点仍然保留，记录下其位置（8j-7，8j-3）和（8j-1，8j-3），然后在三阶初始矩阵里面找到对应的位置并分别赋值插值系数
Figure 408782DEST_PATH_IMAGE093
和
Figure 925214DEST_PATH_IMAGE094
；
Figure 641498DEST_PATH_IMAGE095
和
Figure 977801DEST_PATH_IMAGE096
表示x方向、y方向第8j-3个小长方体单元上计算点被插值的二次电场，该计算点将会被去除，记录下其位置（8j-3，8j-3），然后在三阶初始矩阵里面找到对应的位置并赋值0，以此规律，直至三阶初始矩阵中所有位置均赋值完毕，得到最终的三阶插值矩阵s 3
此时待求解的二次电场E a 被修正为E a,new ：
Figure 437732DEST_PATH_IMAGE097
右端项b被修正为b a,new ：
Figure 808671DEST_PATH_IMAGE098
其中
Figure 692926DEST_PATH_IMAGE099
表示第几阶稀疏采样。
将每个子区域的插值算子S k (k=0,1,2或3)总体合成后得到一个总的稀疏采样算子S，此时E a,new =SE a ，因而AE a,new =ASE a ，将矩阵AS作为新的系数矩阵A a,new 。
最终的电场控制方程表示为：
Figure 782105DEST_PATH_IMAGE100
MATLAB是一款求解大规模矩阵的数学软件，本发明通过调用MATLAB里面的稳定双共轭梯度法（bicgstab）求解器，对上述线性方程组进行求解，可以得到修正后的二次电场E a,new 。
下面对本发明提供的三维地质体其电磁场数值模拟方法的精度和效率进行检验。
对于如图2所示的复杂海洋油储体模型，模拟区域范围为：x和y方向均从-15 km到15 km，z方向均从-10 km到15 km；其中空气层的高度为-10 km到0 km，电阻率为1010欧姆米；海水层的高度为0 km到1 km，其电阻率为0.31欧姆米；油储体距离海水层2 km，其大小为 2 km × 2 km × 0.1 km，电阻率为20欧姆米；将该模型区域剖分为96 × 96 × 62个小长方体单元，其中划分整个模拟区域为高空层（包括最上方的6层小长方体单元）、空气-地表层（包括高空层下方的20层小长方体单元）、地层浅部层（包括空气-地表层下方的12层小长方体单元）、地层中部层（包括地层浅部层下方的10层小长方体单元）和地层底部层（包括地层中部层下方的14层小长方体单元）五个子区域，并为这五个子区域分别定义了1阶稀疏度、正常采样、1阶稀疏度、2阶稀疏度和3阶稀疏度。然后计算了该模型下的的电磁场。为了更好的体现本发明的创新性，将本发明和传统算法做了精度和计算时间的比较。传统算法为国际通用的正常采样算法，每条棱边上有一个采样点。对于本实施例中模型区域剖分为96 × 96 × 62个小长方体单元，x方向就会有96×97×63=586656条棱边，y方向就会有97×96×63=586656条棱边，z方向就会有97×97×62=583358条棱边，总共有1756670条棱边。即传统算法在x方向有586656个计算点，在y方向就会有586656个计算点，z方向就会有583358个计算点，总共约有175万个计算点；利用本发明的算法，只需要计算112万个计算点，将减小63万个计算点，这将有利于减小计算时间和优化内存。
本发明实施的三维地质体其电磁场数值模拟方法利用MATLAB语言编程实现，运行程序所用的个人电脑配置为：CPU-Inter Core i7-8700，主频为3.4GHz，运行内存为36GB。图3为本发明所述三维地质体其电磁场数值模拟方法计算的电场的参考解、数值解以及相对误差等值线图，其中(a)为电场的实部的参考解，(b)为电场的实部的数值解，(c)为电场的实部的相对误差；(d)为电场的虚部的参考解，(e)为电场的虚部的数值解，(f)为电场的虚部的相对误差，从图3中可以看出本发明的数值解和参考解高度一致，从相对误差图中可以看出实部最大相对误差为0.004%，虚部最大误差为0.04%，可以看出本发明方法精度很高。
图4为采用本发明所述方法和采用传统方法随着频率的改变所降低计算时间的百分比图，可以看出测试的7个频率，本发明计算速度更快，可以降低35%到65%的计算时间，特别是在0.25Hz和1Hz，本发明计算速度是传统方法的2倍。值得注意的是随着研究区域大小的增加，本发明方法的速度优势会更加明显。
本发明一实施例中提供一种三维地质体其电磁场数值模拟装置，包括：
第一模块，用于构建内部包含所述勘探目标的三维长方体模型，所述勘探目标为三维异常体；
第二模块，用于对所述三维长方体模型沿x、y、z方向进行网格剖分，剖分成若干个小长方体单元，得到三维长方体模型的网格剖分参数，根据所述勘探目标的电阻率分布，给每个小长方体单元的电阻率赋值，每一个小长方体单元的电阻率为常值，不同小长方体单元的电阻率值不同，得到刻画任意电阻率分布的三维异常体模型；
第三模块，用于根据频率参数和三维异常体模型，构建其对应的一次电场和二次电场间的电场控制方程AE a =b，其中E a 代表二次电场，A为电场控制方程的系数矩阵，
Figure 780148DEST_PATH_IMAGE101
，
Figure 271172DEST_PATH_IMAGE102
表示旋度算子，
Figure 657154DEST_PATH_IMAGE103
表示双旋度算子，
Figure 843416DEST_PATH_IMAGE104
，ω表示角频率，通过ω=2πf求取，f表示预先给定的频率，μ表示磁导率，其值为4π×10-7，ρ表示每个小长方体单元的电阻率，
Figure 769783DEST_PATH_IMAGE105
，
Figure 990680DEST_PATH_IMAGE106
表示异常体电阻率，E b 表示一次电场；
第四模块，用于从三维异常体模型中每个小长方体单元的电阻率中分离出异常体电阻率，得到背景电阻率；根据三维异常体模型中每个小长方体单元的背景电阻率、所述三维长方体模型的网格剖分参数、频率参数计算三维异常体模型的一次电场，进而得到其对应的电场控制方程的右端项b；
第五模块，将所述三维异常体模型沿Z方向划分为多层子区域，对各层子区域分别沿水平方向构建不同稀疏度的插值算子；
第六模块，用于将各层子区域对应的插值算子合成后得到一个总的稀疏采样算子，电场控制方程的系数矩阵A和右端项b分别乘上总的稀疏采样算子得到新的稀疏矩阵A a,new 和修正后的右端项b a,new ，将新的稀疏矩阵A a,new 和修正后的右端项b a,new 带入电场控制方程，求解得到修正后的二次电场E a,new ；
第七模块，用于对所述修正后的二次电场和一次电场求和得到总的电场E，由总的电场E即可计算得到电磁场
Figure 344301DEST_PATH_IMAGE107
。
上述各模块功能的实现方法，可以采用前述各实施例中相同的方法实现，在此不再赘述。
在本实施例中，提供了一种计算机设备，该计算机设备可以是服务器，其内部结构图可以如图5所示。该计算机设备包括通过系统总线连接的处理器、存储器、网络接口和数据库。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统、计算机程序和数据库。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的数据库用于存储样本数据。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现上述三维地质体其电磁场数值模拟方法。
本领域技术人员可以理解，图5中示出的结构，仅仅是与本申请方案相关的部分结构的框图，并不构成对本申请方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。
在一个实施例中，提供了一种计算机设备，包括存储器和处理器，该存储器存储有计算机程序，该处理器执行计算机程序时实现上述实施例中三维地质体其电磁场数值模拟方法的步骤。
在一个实施例中，提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述实施例中三维地质体其电磁场数值模拟方法的步骤。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器（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、y、z方向进行网格剖分，剖分成若干个小长方体单元，得到三维长方体模型的网格剖分参数，根据所述勘探目标的电阻率分布，给每个小长方体单元的电阻率赋值，每一个小长方体单元的电阻率为常值，不同小长方体单元的电阻率值不同，得到刻画任意电阻率分布的三维异常体模型；
根据频率参数和三维异常体模型，构建其对应的一次电场和二次电场间的电场控制方程AE a =b，其中E a 代表二次电场，A为电场控制方程的系数矩阵，
Figure 68478DEST_PATH_IMAGE001
，
Figure 251197DEST_PATH_IMAGE002
表示旋度算子，
Figure 703039DEST_PATH_IMAGE003
表示双旋度算子，
Figure 629406DEST_PATH_IMAGE004
，
Figure 850303DEST_PATH_IMAGE006
表示角频率，通过
Figure 203924DEST_PATH_IMAGE007
求取，f表示预先给定的频率，μ表示磁导率，其值为4π×10-7，ρ表示每个小长方体单元的电阻率，
Figure 877482DEST_PATH_IMAGE008
，
Figure 341961DEST_PATH_IMAGE009
表示异常体电阻率，E b 表示一次电场；
从三维异常体模型中每个小长方体单元的电阻率中分离出异常体电阻率，得到背景电阻率；根据三维异常体模型中每个小长方体单元的背景电阻率、所述三维长方体模型的网格剖分参数和频率参数计算三维异常体模型的一次电场，进而得到其对应的电场控制方程的右端项b；
将所述三维异常体模型沿Z方向划分为多层子区域，对各层子区域分别沿水平方向构建不同稀疏度的插值算子；
将各层子区域对应的插值算子合成后得到一个总的稀疏采样算子，电场控制方程的系数矩阵A和右端项b分别乘上总的稀疏采样算子得到新的稀疏矩阵A a,new 和修正后的右端项b a,new ，将新的稀疏矩阵A a,new 和修正后的右端项b a,new 带入电场控制方程，求解得到修正后的二次电场E a,new ；
对所述修正后的二次电场和一次电场求和得到总的电场E，由总的电场E即可计算得到磁场
Figure 417365DEST_PATH_IMAGE010
。
2.根据权利要求1所述的三维地质体其电磁场数值模拟方法，其特征在于，电场控制方程的系数矩阵A通过以下方法获得：
三维异常体模型中各个小长方体单元x、y、z方向的棱边长度，构成长度元矩阵
Figure 613991DEST_PATH_IMAGE011
；
三维异常体模型中各个小长方体单元的表面积，构成面积元矩阵S；
三维异常体模型中各个小长方体单元的体积，构成体积元矩阵V；
双旋度算子
Figure 899479DEST_PATH_IMAGE012
，其中O为拓扑矩阵，由-1和1两种元素组成；
最后得到电场控制方程的系数矩阵
Figure 777436DEST_PATH_IMAGE013
。
3.根据权利要求1所述的三维地质体其电磁场数值模拟方法，其特征在于，一次电场E b 的三个分量分别是E r 、E φ 、E z ，通过下式求解：
Figure 97559DEST_PATH_IMAGE014
其中：A 0 是矢量位，A r 、A φ 和A z 为A 0 的三个分量；ρ b表示背景电阻率模型；r表示网格剖分参数产生的极矩。
4.根据权利要求1、2或3所述的三维地质体其电磁场数值模拟方法，其特征在于，将所述三维异常体模型沿Z方向划分为五层子区域，由上至下分别为高空层子区域、空气-地表层子区域、浅地层子区域、中地层子区域和深地层子区域；
对于高空层子区域和浅地层子区域，其插值算子为1阶稀疏采样算子；
对于空气-地表层子区域其插值算子为单位矩阵I；
对于中地层子区域其插值算子为2阶稀疏采样算子；
对于深地层子区域其插值算子为3阶稀疏采样算子。
5.根据权利要求4所述的三维地质体其电磁场数值模拟方法，其特征在于，高空层子区域和浅地层子区域采用相同的方法分别构建1阶稀疏采样算子S1，S1=s1，s 1为一阶插值矩阵;
构建一阶插值矩阵s 1，方法是：先构建一个一阶初始矩阵，该矩阵大小为(n,n)，主对角线的元素为1；
对子区域中的各小长方体单元，构建x和y方向一阶插值算子系数为：
Figure 933928DEST_PATH_IMAGE015
Figure 441133DEST_PATH_IMAGE016
其中
Figure 119851DEST_PATH_IMAGE017
、
Figure 294481DEST_PATH_IMAGE018
分别表示子区域中x方向上第2j-1、2j+1个小长方体单元上计算点对应的插值算子系数；
Figure 301751DEST_PATH_IMAGE019
、
Figure 561831DEST_PATH_IMAGE020
分别表示子区域中y方向第2j-1、2j+1个小长方体单元计算点对应的插值算子系数，而
Figure 47170DEST_PATH_IMAGE021
和
Figure 810727DEST_PATH_IMAGE022
分别代表第2j、2j-1、2j+1个小长方体单元沿x和y方向的棱边长度；
将所构建的一阶插值算子系数分别对应施加到子区域中的相应小长方体单元上，有
Figure 988898DEST_PATH_IMAGE023
Figure 736274DEST_PATH_IMAGE024
Figure 25304DEST_PATH_IMAGE025
和
Figure 908947DEST_PATH_IMAGE026
表示沿x方向上第2j-1、2j+1个小长方体单元上计算点对应的二次电场，该计算点仍然会保留，记录下其位置（2j-1，2j）和（2j+1，2j），然后在一阶初始矩阵里面找到对应的位置并分别赋值插值系数
Figure 258020DEST_PATH_IMAGE027
、
Figure 227113DEST_PATH_IMAGE028
；
Figure 54254DEST_PATH_IMAGE029
和
Figure 57982DEST_PATH_IMAGE030
表示沿y方向上第2j-1、2j+1个小长方体单元上计算点对应的二次电场，该计算点仍然保留，记录下其位置（2j，2j-1）和（2j，2j+1）然后在一阶初始矩阵里面找到对应的位置并分别赋值插值系数
Figure 577957DEST_PATH_IMAGE031
、
Figure 34346DEST_PATH_IMAGE032
；
Figure 668108DEST_PATH_IMAGE033
和
Figure 260763DEST_PATH_IMAGE034
表示x方向、y方向第2j个小长方体单元上计算点被插值的二次电场，该计算点将会被去除，记录下其位置（2j，2j），然后在一阶初始矩阵里面找到对应的位置并赋值0，以此规律，直至一阶初始矩阵中所有位置均赋值完毕，得到最终的一阶插值矩阵s 1。
6.根据权利要求5所述的三维地质体其电磁场数值模拟方法，其特征在于，对于中地层子区域其2阶稀疏采样算子，表示为
Figure 951639DEST_PATH_IMAGE035
，s 2为二阶插值矩阵；
构建二阶插值矩阵s 2，方法是：先构建一个二阶初始矩阵，其矩阵大小为(n,n)，主对角线元素为1；
对中地层子区域中的各小长方体单元，构建x和y方向二阶插值系数为：
Figure 629745DEST_PATH_IMAGE036
Figure 329848DEST_PATH_IMAGE037
其中
Figure 42589DEST_PATH_IMAGE038
、
Figure 373207DEST_PATH_IMAGE039
分别表示中地层子区域中x方向上第4j-3、4j+1个小长方体单元上计算点对应的插值系数；
Figure 804188DEST_PATH_IMAGE040
、
Figure 776824DEST_PATH_IMAGE041
分别表示中地层子区域中y方向第4j-3、4j+1个小长方体单元计算点对应的插值系数，
Figure 609650DEST_PATH_IMAGE042
和
Figure 376749DEST_PATH_IMAGE043
分别代表第4j-3、4j-2、4j-1、4j个小长方体单元沿x和y方向的棱边长度；
将所构建的二阶插值算子分别对应施加到中地层子区域中的相应小长方体单元上，有：
Figure 295027DEST_PATH_IMAGE044
Figure 71353DEST_PATH_IMAGE045
Figure 493107DEST_PATH_IMAGE046
和
Figure 696686DEST_PATH_IMAGE047
表示中地层子区域沿x方向上第4j-3、4j+1个小长方体单元上计算点对应的二次电场，该计算点仍然会保留，记录下其位置（4j-1，4j-3）和（4j-1，4j+1），然后在二阶初始矩阵里面找到对应的位置并分别赋值插值系数
Figure 836680DEST_PATH_IMAGE048
、
Figure 413768DEST_PATH_IMAGE049
；
Figure 955608DEST_PATH_IMAGE050
和
Figure 330088DEST_PATH_IMAGE051
表示沿y方向上第4j-3、4j+1个小长方体单元上计算点对应的电场，该计算点仍然保留，记录下其位置（4j-3，4j-1）和（4j+1，4j-1），然后在二阶初始矩阵里面找到对应的位置并分别赋值插值系数
Figure 957379DEST_PATH_IMAGE052
、
Figure 75507DEST_PATH_IMAGE053
；
Figure 471854DEST_PATH_IMAGE054
和
Figure 17236DEST_PATH_IMAGE055
表示x方向、y方向第4j-1个小长方体单元上计算点被插值的电场，该计算点将会被去除，记录下其位置（4j-1，4j-1），然后在二阶初始矩阵里面找到对应的位置并赋值0，以此规律，直至二阶初始矩阵中所有位置均赋值完毕，得到最终的二阶插值矩阵s 2。
7.根据权利要求6所述的三维地质体其电磁场数值模拟方法，其特征在于，对于深地层子区域其3阶稀疏采样算子，表示为
Figure 131822DEST_PATH_IMAGE056
，
Figure 788063DEST_PATH_IMAGE058
为二阶插值矩阵；
构建三阶插值矩阵
Figure 304495DEST_PATH_IMAGE058
，其方法是：先构建一个三阶初始矩阵，其矩阵大小为(n,n)，主对角线的元素为1；
对深地层子区域中的小长方体单元，然后构建x和y方向三阶插值系数为：
Figure 20778DEST_PATH_IMAGE059
Figure 622660DEST_PATH_IMAGE060
其中
Figure 82592DEST_PATH_IMAGE061
、
Figure 453530DEST_PATH_IMAGE062
分别表示深地层子区域中x方向上第8j-7、8j+1个小长方体单元上计算点对应的插值系数；
Figure 340715DEST_PATH_IMAGE063
、
Figure 164314DEST_PATH_IMAGE064
分别表示深地层子区域中y方向第8j-7、8j+1个小长方体单元计算点对应的插值系数，
Figure 419147DEST_PATH_IMAGE065
和
Figure 644592DEST_PATH_IMAGE066
Figure 702678DEST_PATH_IMAGE067
分别代表第8j-7、8j-6、8j-5、8j-4、8j-3、8j-2、8j-1、8j个小长方体单元沿x和y方向的棱边长度；
将所构建的三阶插值算子分别对应施加到深地层子区域中的相应小长方体单元上，有：
Figure 279153DEST_PATH_IMAGE068
Figure 80887DEST_PATH_IMAGE069
Figure 160839DEST_PATH_IMAGE070
和
Figure 655405DEST_PATH_IMAGE071
表示深地层子区域沿x方向上第8j-7、8j-1个小长方体单元上计算点对应的二次电场，该计算点仍然会保留，记录下其位置（8j-7，8j-3）和（8j-1，8j-3），然后在三阶初始矩阵里面找到对应的位置并分别赋值插值系数
Figure 453597DEST_PATH_IMAGE072
和
Figure 59021DEST_PATH_IMAGE073
；
Figure 993479DEST_PATH_IMAGE074
和
Figure 658947DEST_PATH_IMAGE075
表示沿y方向上第8j-7、8j+1个小长方体单元上计算点对应的电场，该计算点仍然保留，记录下其位置（8j-7，8j-3）和（8j-1，8j-3），然后在三阶初始矩阵里面找到对应的位置并分别赋值插值系数
Figure 944435DEST_PATH_IMAGE076
和
Figure 87971DEST_PATH_IMAGE077
；
Figure 142515DEST_PATH_IMAGE078
和
Figure 244463DEST_PATH_IMAGE079
表示x方向、y方向第8j-3个小长方体单元上计算点被插值的二次电场，该计算点将会被去除，记录下其位置（8j-3，8j-3），然后在三阶初始矩阵里面找到对应的位置并赋值0，以此规律，直至三阶初始矩阵中所有位置均赋值完毕，得到最终的三阶插值矩阵
Figure 751668DEST_PATH_IMAGE080
。
8.三维地质体其电磁场数值模拟装置，其特征在于，包括：
第一模块，用于构建内部包含所述勘探目标的三维长方体模型，所述勘探目标为三维异常体；
第二模块，用于对所述三维长方体模型沿x、y、z方向进行网格剖分，剖分成若干个小长方体单元，得到三维长方体模型的网格剖分参数，根据异常体的电阻率分布，给每个小长方体单元的电阻率赋值，每一个小长方体单元的电阻率为常值，不同小长方体单元的电阻率值不同，得到刻画任意电阻率分布的三维异常体模型；
第三模块，用于根据频率参数和三维异常体模型，构建其对应的一次电场和二次电场间的电场控制方程AE a =b，其中E a 代表二次电场，A为电场控制方程的系数矩阵，
Figure 695966DEST_PATH_IMAGE081
，
Figure 870595DEST_PATH_IMAGE082
表示旋度算子，
Figure 612286DEST_PATH_IMAGE083
表示双旋度算子，
Figure 872366DEST_PATH_IMAGE084
，
Figure 357705DEST_PATH_IMAGE086
表示角频率，通过
Figure 386841DEST_PATH_IMAGE087
求取，f表示给定的频率，μ表示磁导率，其值为4π×10-7，ρ表示每个小长方体单元的电阻率，
Figure 299434DEST_PATH_IMAGE088
，
Figure 46810DEST_PATH_IMAGE089
表示异常体电阻率，E b 表示一次电场；
第四模块，用于从三维异常体模型中每个小长方体单元的电阻率中分离出异常体电阻率，得到背景电阻率；根据三维异常体模型中每个小长方体单元的背景电阻率、所述三维长方体模型的网格剖分参数和频率参数计算三维异常体模型的一次电场，进而得到其对应的电场控制方程的右端项b；
第五模块，用于将所述三维异常体模型沿Z方向划分为多层子区域，对各层子区域分别沿水平方向构建不同稀疏度的插值算子；
第六模块，用于将各层子区域对应的插值算子合成后得到一个总的稀疏采样算子，电场控制方程的系数矩阵A和右端项b分别乘上总的稀疏采样算子得到新的稀疏矩阵A a,new 和修正后的右端项b a,new ，将新的稀疏矩阵A a,new 和修正后的右端项b a,new 带入电场控制方程，求解得到修正后的二次电场E a,new ；
第七模块，用于对所述修正后的二次电场和一次电场求和得到总的电场E，由总的电场E即可计算得到磁场
Figure 335840DEST_PATH_IMAGE090
。
9.一种计算机设备，包括存储器和处理器，存储器存储有计算机程序，其特征在于：处理器执行计算机程序时实现权利要求1所述三维地质体其电磁场数值模拟方法中的步骤。
10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于：计算机程序被处理器执行时实现权利要求1所述三维地质体其电磁场数值模拟方法中的步骤。
CN202111347154.4A 2021-11-15 2021-11-15 三维地质体其电磁场数值模拟方法、装置、设备及介质 Active CN113779818B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111347154.4A CN113779818B (zh) 2021-11-15 2021-11-15 三维地质体其电磁场数值模拟方法、装置、设备及介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111347154.4A CN113779818B (zh) 2021-11-15 2021-11-15 三维地质体其电磁场数值模拟方法、装置、设备及介质

Publications (2)

Publication Number Publication Date
CN113779818A true CN113779818A (zh) 2021-12-10
CN113779818B CN113779818B (zh) 2022-02-08

Family

ID=78873940

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111347154.4A Active CN113779818B (zh) 2021-11-15 2021-11-15 三维地质体其电磁场数值模拟方法、装置、设备及介质

Country Status (1)

Country Link
CN (1) CN113779818B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117538945A (zh) * 2024-01-10 2024-02-09 中南大学 三维大地电磁多分辨率反演方法、装置、设备及介质

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090006053A1 (en) * 2006-03-08 2009-01-01 Carazzone James J Efficient Computation Method for Electromagnetic Modeling
CN103412989A (zh) * 2013-08-01 2013-11-27 电子科技大学 基于参数化降阶模型周期结构的三维电磁场仿真模拟方法
CN108873084A (zh) * 2018-05-10 2018-11-23 中南大学 一种基于单位分解积分的直流电阻率无单元正演方法
CN110346834A (zh) * 2019-07-22 2019-10-18 中国科学院地球化学研究所 三维频率域可控源电磁的正演方法、系统
US20210096203A1 (en) * 2019-09-27 2021-04-01 Q Bio, Inc. Maxwell parallel imaging
CN113221393A (zh) * 2021-01-29 2021-08-06 吉林大学 一种基于非结构有限元法的三维大地电磁各向异性反演方法
CN113553748A (zh) * 2021-09-22 2021-10-26 中南大学 一种三维大地电磁正演数值模拟方法
CN113569447A (zh) * 2021-07-06 2021-10-29 武汉市市政建设集团有限公司 一种基于舒尔补方法的瞬变电磁三维快速正演方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090006053A1 (en) * 2006-03-08 2009-01-01 Carazzone James J Efficient Computation Method for Electromagnetic Modeling
CN103412989A (zh) * 2013-08-01 2013-11-27 电子科技大学 基于参数化降阶模型周期结构的三维电磁场仿真模拟方法
CN108873084A (zh) * 2018-05-10 2018-11-23 中南大学 一种基于单位分解积分的直流电阻率无单元正演方法
CN110346834A (zh) * 2019-07-22 2019-10-18 中国科学院地球化学研究所 三维频率域可控源电磁的正演方法、系统
US20210096203A1 (en) * 2019-09-27 2021-04-01 Q Bio, Inc. Maxwell parallel imaging
CN113221393A (zh) * 2021-01-29 2021-08-06 吉林大学 一种基于非结构有限元法的三维大地电磁各向异性反演方法
CN113569447A (zh) * 2021-07-06 2021-10-29 武汉市市政建设集团有限公司 一种基于舒尔补方法的瞬变电磁三维快速正演方法
CN113553748A (zh) * 2021-09-22 2021-10-26 中南大学 一种三维大地电磁正演数值模拟方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
LI JIAN，ET AL.: "A Comparison of Different Divergence-free Solutions for 3D Anisotropic CSEM Modeling Using Staggered Finite Difference Method", 《THE 9TH INTERNATIONAL CONFERENCE ON ENVIRONMENTAL AND ENGINEERING GEOPHYSICS》 *
XULONG WANGA, ET AL.: "Fast numerical simulation of 2D gravity anomaly based on nonuniform fast Fourier transform in mixed space-wavenumber domain", 《JOURNAL OF APPLIED GEOPHYSICS》 *
徐凌华等: "基于有限单元法的二维／三维大地", 《物探化探计算技术》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117538945A (zh) * 2024-01-10 2024-02-09 中南大学 三维大地电磁多分辨率反演方法、装置、设备及介质
CN117538945B (zh) * 2024-01-10 2024-03-26 中南大学 三维大地电磁多分辨率反演方法、装置、设备及介质

Also Published As

Publication number Publication date
CN113779818B (zh) 2022-02-08

Similar Documents

Publication Publication Date Title
CN112800657B (zh) 基于复杂地形的重力场数值模拟方法、装置和计算机设备
CN112287534B (zh) 基于nufft的二维磁异常快速正演模拟方法和装置
CN113553748B (zh) 一种三维大地电磁正演数值模拟方法
CN114065511A (zh) 起伏地形下大地电磁二维正演数值模拟方法、装置、设备及介质
CN111103627B (zh) 大地电磁tm极化模式对电场数据二维反演方法和装置
CN110346835B (zh) 大地电磁的正演方法、正演系统、存储介质及电子设备
CN113255230B (zh) 基于mq径向基函数的重力模型正演方法及系统
CN113420487B (zh) 一种三维重力梯度张量数值模拟方法、装置、设备和介质
CN113779818B (zh) 三维地质体其电磁场数值模拟方法、装置、设备及介质
CN113962077B (zh) 三维各向异性强磁场数值模拟方法、装置、设备及介质
CN111967169B (zh) 二度体重力异常积分解数值模拟方法和装置
CN114021408A (zh) 二维强磁场数值模拟方法、装置、设备及介质
CN111103628B (zh) 大地电磁te极化模式对电场数据二维反演方法和装置
CN113051779B (zh) 一种三维直流电阻率法数值模拟方法
CN116842813B (zh) 一种三维大地电磁正演数值模拟方法
CN111339688B (zh) 基于大数据并行算法求解火箭仿真模型时域方程的方法
CN114970289B (zh) 三维大地电磁各向异性正演数值模拟方法、设备及介质
CN110968930B (zh) 一种地质体变属性插值方法及系统
CN113673163B (zh) 一种三维磁异常数快速正演方法、装置和计算机设备
CN114036806A (zh) 基于热导率各向异性介质的三维地温场数值模拟方法
CN113204905B (zh) 一种接触式激发极化法有限单元数值模拟方法
CN114114438A (zh) 一种回线源地空瞬变电磁数据的拟三维反演方法
CN113806686B (zh) 大规模复杂地质体重力梯度快速计算方法、装置和设备
CN117538945B (zh) 三维大地电磁多分辨率反演方法、装置、设备及介质
CN110543611B (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