CN110346834B - 三维频率域可控源电磁的正演方法、系统 - Google Patents

三维频率域可控源电磁的正演方法、系统 Download PDF

Info

Publication number
CN110346834B
CN110346834B CN201910663566.5A CN201910663566A CN110346834B CN 110346834 B CN110346834 B CN 110346834B CN 201910663566 A CN201910663566 A CN 201910663566A CN 110346834 B CN110346834 B CN 110346834B
Authority
CN
China
Prior art keywords
electric field
field
determining
simulation
simulation area
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.)
Expired - Fee Related
Application number
CN201910663566.5A
Other languages
English (en)
Other versions
CN110346834A (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.)
Institute of Geochemistry of CAS
Original Assignee
Institute of Geochemistry of CAS
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 Institute of Geochemistry of CAS filed Critical Institute of Geochemistry of CAS
Priority to CN201910663566.5A priority Critical patent/CN110346834B/zh
Publication of CN110346834A publication Critical patent/CN110346834A/zh
Application granted granted Critical
Publication of CN110346834B publication Critical patent/CN110346834B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/08—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation operating with magnetic or electric fields produced or modified by objects or geological structures or by detecting devices
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/38—Processing data, e.g. for analysis, for interpretation, for correction

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Electromagnetism (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本申请提供一种三维频率域可控源电磁的正演方法、系统、存储介质及电子设备，方法包括：确定模拟区域并将模拟区域剖分为多个单元；获取每个单元的电导率张量和磁导率张量；确定出模拟区域内的一次场电场；根据模拟区域内的一次场电场及电导率张量和磁导率张量，确定出模拟区域内的二次场电场；进一步确定出模拟区域的电场和磁场；并由此确定出模拟区域内地表处的视电阻率及相位。通过在三维频率域可控源电磁的正演方法中考虑电导率和磁导率的各向异性，可以定量地分析电导率与磁导率各向异性对三维频率域可控源电磁响应的影响，这有利于更加准确地了解和分析实际中复杂的地质结构的电磁响应，因此适用于更多类型的地质结构，应用领域也更加广阔。

Description

三维频率域可控源电磁的正演方法、系统
技术领域
本申请涉及勘探地球物理领域，具体而言，涉及一种三维频率域可控源电磁的正演方法、系统、存储介质及电子设备。
背景技术
频率域可控源电磁法是一类非常重要的地球物理方法，具有横向分辨能力强、工作效率高及信噪比高等优点，目前被广泛应用于矿产资源勘查、油气勘探及工程勘察等领域。在实际的勘查中，通常对观测数据进行数据模拟，以确定地下介质的电性结构。但目前的正演方法及正演系统的应用范围有限，而基于现有的正演系统对具有复杂电磁响应的地下介质进行的反演，难以有效地应用在更复杂的环境中。
发明内容
本申请实施例的目的在于提供一种三维频率域可控源电磁的正演方法、系统、存储介质及电子设备，以使基于频率域可控源电磁法的正演方法和正演系统应用于更复杂的环境中。
为了实现上述目的，本申请的实施例通过如下方式实现：
第一方面，本申请的实施例提供一种三维频率域可控源电磁的正演方法，包括：确定模拟区域，所述模拟区域包括表征异常体的异常体区域及与所述异常体区域互补的周边区域，所述异常体区域内的电导率和磁导率不同于所述周边区域的电导率和磁导率；将所述模拟区域剖分为多个单元；获取每个单元的电导率张量
Figure BDA0002139012770000011
和磁导率张量
Figure BDA0002139012770000012
确定出模拟区域内的一次场电场，其中，所述一次场电场表示排除所述异常体时所述模拟区域的电场；根据每个单元的一次场电场、
Figure BDA0002139012770000013
和
Figure BDA0002139012770000014
确定出模拟区域内的二次场电场，其中，所述二次场电场为所述异常体产生的电场；根据模拟区域内的一次场电场和二次场电场，确定出所述模拟区域的电场和磁场；根据所述模拟区域的电场和磁场，确定出所述模拟区域内地表处的视电阻率及相位。
通过在三维频率域可控源电磁的正演方法中考虑电导率和磁导率的各向异性，可以定量地分析电导率与磁导率各向异性对三维频率域可控源电磁响应的影响，这有利于更加准确地了解和分析实际中复杂的地质结构的电磁响应。并且，本方法也可以为三维频率域可控源电磁电导率与磁导率各向异性反演提供基础。以及，由于本方法同时考虑了电导率和磁导率的各向异性，适用于更多类型的地质结构，因此，可以在更复杂的环境中应用。
结合第一方面，在第一方面的第一种可能的实现方式中，所述将所述模拟区域剖分为多个单元，包括：采用稀疏网格将所述异常体区域均匀剖分，以及，沿远离所述异常体的方向采用逐步增大的网格对所述周边区域进行剖分，以将所述模拟区域剖分为多个单元。
通过这样的方式对模拟区域进行剖分，能尽可能减小异常体区域在模拟区域的外边界处的影响，甚至可以忽略不计。另外，采用稀疏网格配置方法，可以降低矩阵方程维度，节约求解过程花费的时间，从而提高效率。
结合第一方面，在第一方面的第二种可能的实现方式中，所述获取每个单元的电导率张量
Figure BDA0002139012770000021
和磁导率张量
Figure BDA0002139012770000022
包括：获取所述模拟区域内每个单元的主轴各向异性坐标系中x”主轴方向的电导率σ1、y”’主轴方向的电导率σ2、z’主轴方向的电导率σ3，并获取三次欧拉旋转过程中每次旋转的角度αS、αD和αL，以及，获取所述模拟区域内每个单元的主轴各向异性坐标系中x”主轴方向的磁化率χ1、y”’主轴方向的磁化率χ2、z’主轴方向的磁化率χ3，并获取三次欧拉旋转过程中每次旋转的角度βS、βD和βL；根据每个单元的σ1、σ2、σ3、αS、αD和αL，确定出每个单元的
Figure BDA0002139012770000023
以及，根据每个单元的χ1、χ2、χ3、βS、βD和βL，确定出每个单元的
Figure BDA0002139012770000024
主轴各向异性坐标系可以体现电导率的各向异性和磁导率的各向异性，是实现将电导率和磁导率的各向异性考虑进正演方法中的基础。
结合第一方面，在第一方面的第三种可能的实现方式中，所述根据模拟区域内的一次场电场、
Figure BDA0002139012770000025
和
Figure BDA0002139012770000026
确定出模拟区域内的二次场电场，包括：根据模拟区域内的一次场电场、
Figure BDA0002139012770000027
和
Figure BDA0002139012770000028
结合麦克斯韦方程组，确定所述模拟区域的电场积分方程；根据所述电场积分方程，确定出所述模拟区域的总体系数矩阵；根据预设的边界条件和所述总体系数矩阵，确定出所述模拟区域内的二次场电场。
采用这种方式可以实现基于模拟区域内的一次场电场、
Figure BDA0002139012770000031
和
Figure BDA0002139012770000032
确定出模拟区域内的二次场电场，使得二次场电场中能够体现异常区域电导率和磁导率的各向异性。
结合第一方面的第三种实现方式，在第一方面的第四种可能的实现方式中，所述根据模拟区域内的一次场电场、
Figure BDA0002139012770000033
和
Figure BDA0002139012770000034
结合麦克斯韦方程组，确定所述模拟区域的电场积分方程，包括：根据模拟区域内的一次场电场，确定出每个单元的电场双旋度方程；由根据每个单元的电场双旋度方程及
Figure BDA0002139012770000035
和
Figure BDA0002139012770000036
确定出的方程点乘电场强度的变分后对整个模拟区域积分，得到中间方程；利用矢量恒等式和散度定理，基于所述中间方程确定出所述模拟区域的电场积分方程。
采用这种方式可以准确地得到对模拟区域积分的电场积分方程。
结合第一方面的第三种实现方式，在第一方面的第五种可能的实现方式中，在每个单元为六面体时，电场赋在棱边上，所述根据所述电场积分方程，确定出所述模拟区域的总体系数矩阵，包括：根据所述电场积分方程，确定出每个单元的单元方程；根据每个单元的单元方程中包含的系数矩阵，确定出所述模拟区域的总体系数矩阵。
采用这种方式，可以允许电场法向分量在电性不连续界面处突变，且考虑了电导率和磁导率的各向异性，更符合实际情况，从而使应用此正演方法具有更加广泛的应用范围和更高的准确性。
结合第一方面，在第一方面的第六种可能的实现方式中，所述根据模拟区域内的一次场电场和二次场电场，确定出所述模拟区域的电场和磁场，包括：基于模拟区域内的一次场电场和二次场电场，确定出模拟区域内的电场；每条棱边上的磁场可根据该棱边及其邻近棱边上的电场确定出，其中，所述模拟区域内所有单元的电场和磁场的总体表示所述模拟区域的电场和磁场。
第二方面，本申请的实施例提供一种三维频率域可控源电磁的正演系统，包括：区域确定单元，用于确定模拟区域，所述模拟区域包括表征异常体的异常体区域及与所述异常体区域互补的周边区域，所述异常体区域内的电导率和磁导率不同于所述周边区域的电导率和磁导率；区域剖分单元，用于将所述模拟区域剖分为多个单元；张量获取单元，用于获取每个单元的电导率张量
Figure BDA0002139012770000037
和磁导率张量
Figure BDA0002139012770000038
一次场确定单元，用于确定出模拟区域内的一次场电场，其中，所述一次场电场表示排除所述异常体时所述模拟区域的电场；二次场确定单元，用于根据模拟区域内的一次场电场、
Figure BDA0002139012770000041
和
Figure BDA0002139012770000042
确定出模拟区域内的二次场电场，其中，所述二次场电场为所述异常体产生的电场；正演处理单元，用于根据模拟取悦的一次场电场和二次场电场，确定出所述模拟区域的电场和磁场；以及，根据所述模拟区域的电场和磁场，确定出所述模拟区域内地表处的视电阻率及相位。
第三方面，本申请的实施例提供一种存储介质，所述存储介质包括存储的程序，其中，在所述程序运行时控制所述存储介质所在设备执行第一方面或第一方面的第一种至第六种可能的实现方式中任一所述的三维频率域可控源电磁的正演方法。
第四方面，本申请的实施例提供一种电子设备，包括存储器和处理器，所述存储器用于存储包括程序指令的信息，所述处理器用于控制程序指令的执行，所述程序指令被处理器加载并执行时实现第一方面或第一方面的第一种至第六种可能的实现方式中任一所述的三维频率域可控源电磁的正演方法的步骤。
为使本申请的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。
附图说明
为了更清楚地说明本申请实施例的技术方案，下面将对本申请实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本申请的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。
图1为本申请实施例提供的一种三维频率域可控源电磁的正演方法的流程图。
图2示出了本申请实施例提供的一种模拟区域的示意图。
图3示出了本申请实施例提供的一种六面体单元的结构框图。
图4示出了本申请实施例提供的一种欧拉旋转的方法。
图5示出了本申请实施例提供的一种矢量有限元法应用于六面体单元的示意图。
图6示出了本申请实施例提供的一种电阻率任意各向异性模型的示意图。
图7示出了本申请实施例提供的电阻率任意各向异性模型的测试结果收敛状况示意图。
图8示出了本申请实施例提供的有限元的结果和有限差分的结果的对比图。
图9示出了本申请实施例提供的一种磁导率层状介质模型的示意图。
图10示出了本申请实施例提供的磁导率层状介质模型的解析结果对比的示意图。
图11示出了本申请实施例提供的一种三维频率域可控源电磁的正演系统的结构框图。
图12示出了本申请实施例提供的一种电子设备的结构框图。
图示：10-三维频率域可控源电磁的正演系统；11-区域确定单元；12-区域剖分单元；13-张量获取单元；14-一次场确定单元；15-二次场确定单元；16-正演处理单元；20-电子设备；21-存储器；22-通信模块；23-总线；24-处理器。
具体实施方式
下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行描述。
频率域可控源电磁法在实际勘察的应用中，由于发射源的强度相对于天然场源信号强度明显增强，对针对目标体的勘探分辨能力提高，观测电磁场值所需的时间大大减少，所以该方法具有较高的生产效率。但目前的频率域可控源电磁的正演系统多数是基于地球介质是各向同性的这一假设建立起来的，例如基于可控源音频大地电磁法的正演系统、基于海洋可控源电磁法的正演系统等。
由于各种宏观微观因素的影响，如地球构造应力场、地球介质形变带及晶体定向排列等，地球介质的各向异性是客观存在的，因此，地球介质是各向同性的这一假设的解释能力有限。而固体地球介质物理性质的各向异性现象已被国内外地球物理学界证实和接受，以这一物理现象为基础发展的地球物理理论与方法技术也已在一些领域获得了初步应用，但距离系统性的大范围应用还存在诸多问题。
本申请的发明人发现，基于现有的频率域可控源电磁正演系统，在其基础上对复杂地质的实际观测数据进行数值模拟时，如果忽略各向异性，可能给反演结果带来偏差甚至错误。基于此，本申请的发明人经过创造性的劳动提出本申请。
其中，可控源音频大地电磁法：一种在大地电磁法和音频大地电磁法的基础上发展起来的人工源电磁勘探方法。海洋可控源电磁法：海洋电磁探测法中的一种，用于对海洋进行探测，利用不同介质的电导率差异来探测海底内部的地质结构。正演：是指对于一个给定的模型，在一定激发源的作用下，根据一定的物理原理求其响应的过程。反演：是指根据实测的数据来反推产生这些数据的系统内在信息的一种数学物理过程。
如图1所示，本申请提供的频率域可控源电磁的正演方法可以包括：步骤S10、步骤S20、步骤S30、步骤S40、步骤S50、步骤S60和步骤S70。
下面对本申请实施例提供的频率域可控源电磁的正演方法进行详细完整的介绍。
步骤S10：确定模拟区域，所述模拟区域包括表征异常体的异常体区域及与所述异常体区域互补的周边区域，所述异常体区域内的电导率和磁导率不同于所述周边区域的电导率和磁导率。
在本实施例中，可以在应用本申请实施例提供的三维频率域可控源电磁正演方法的正演系统中，以获取的模拟数据为基础建立模拟区域；当然，模拟区域也可以是正演系统中预设的模拟区域，此处不作为限定。模拟区域中包括表征异常体的异常体区域，以及与异常体区域互补的周边区域。其中，异常体：是指引起物探异常的地质体，如引起重力异常的不同密度体、引起磁异常的磁性体、引起电法异常的具有电性差异的地质体等。
如图2所示出的模拟区域在正演系统中的示意图，在模拟区域中，包括空气、地面和地下介质，而地下介质中存在异常体。异常体所在的区域表示异常体区域，而模拟区域中与异常体区域互补的区域为周边区域，周边区域可以包括空气、地面和部分地下介质。需要说明的是，在本实施例中，异常体区域的电导率和磁导率可以与周边区域的电导率和磁导率不同，而异常体区域内各部分的电导率和磁导率可以相同，也可以由于各部分的不同而可能不同，此处不作为限定。
步骤S20：将所述模拟区域剖分为多个单元。
在本实施例中，为了尽可能减小甚至消除异常体区域在模拟区域中与周边区域的交界的外边界处的影响，可以采用稀疏网格将异常体区域进行均匀剖分；以及，可以沿远离异常体的方向采用逐步增大的网格对周边区域进行剖分，从而将整个模拟区域剖分为多个单元。
示例性的，可以采用规则六面体剖分网格对异常体区域进行均匀剖分，并从异常体区域开始向周边区域沿x、y及z三个方向不断增大剖分网格的大小，直至异常体区域在模拟区域的外边界处的影响可以忽略不计。通过规则六面体剖分网格对模拟区域的剖分而得到的矩形块如图3所示，矩形块状的单元，也便于分析。当然，同样可以采用其他方式例如四面体剖分切割方法，对模拟区域进行剖分，得到多个四面体单元。因此，此处不应视为对本申请的限定。
通过对模拟区域进行剖分，能尽可能减小异常体区域在模拟区域的外边界处的影响，甚至可以忽略不计。另外，采用稀疏网格配置方法，可以降低矩阵方程维度，节约求解过程花费的时间，从而提高效率。
步骤S30：获取每个单元的电导率张量
Figure BDA0002139012770000071
和磁导率张量
Figure BDA0002139012770000072
在本实施例中，对于电导率和磁导率各向异性的异常体，其电导率张量
Figure BDA0002139012770000073
磁导率张量
Figure BDA0002139012770000074
都可以由主轴各向异性通过欧拉旋转得到。
示例性的，以图4为例，可以为每个单元确定一个测量坐标系(x，y，z)。之后将测量坐标系沿z轴顺时针方向旋转，旋转角度为αS，得到坐标系(x’，y’，z)；然后将坐标系(x’，y’，z)绕x’轴顺时针旋转，旋转角度为αD，得到坐标系(x’，y”，z’)；再将坐标系(x’，y”，z’)绕z’轴旋转，旋转角度为αL，得到主轴各向异性坐标系(x”，y”’，z’)。
确定出主轴各向异性坐标系(x”，y”’，z’)后，可以获取主轴各向异性坐标系中x”主轴方向的电导率σ1、y”’主轴方向的电导率σ2、z’主轴方向的电导率σ3，以及，获取三次欧拉旋转过程中每次旋转的角度αS、αD和αL；其中，αS：各向异性走向角，αD：各向异性倾角，αL：各向异性倾向角。
类似的，对测量坐标系进行三次欧拉旋转，得到主轴各向异性坐标系(x”，y”’，z’)。可以获取主轴各向异性坐标系中x”主轴方向的磁化率χ1、y”’主轴方向的磁化率χ2、z’主轴方向的磁化率χ3，以及，获取三次欧拉旋转过程中每次旋转的角度βS、βD和βL，其中，βS：各向异性走向角，βD：各向异性倾角，βL：各向异性倾向角。需要说明的是，同一单元的βS、βD、βL与αS、αD、αL之间相互独立，旋转角度可以不同。
参数(σ1、σ2、σ3、μ1、μ2、μ3、αS、αD、αL、βS、βD和βL)可以通过预设的方式进行设定，也可以依靠用户输入进行设定。确定出上述参数后，可以根据每个单元的σ1、σ2、σ3、αS、αD和αL，确定出每个单元的电导率张量
Figure BDA0002139012770000081
以及，根据每个单元的χ1、χ2、χ3、βS、βD和βL，确定出每个单元的磁导率张量
Figure BDA0002139012770000082
具体地，在本实施例中，根据每个单元的σ1、σ2、σ3、αS、αD和αL，确定出每个单元的
Figure BDA0002139012770000083
的方式，可以是通过每次旋转的旋转变换矩阵。其中，将测量坐标系沿z轴顺时针方向旋转，旋转角度为αS，得到坐标系(x’，y’，z)，对应的旋转变换矩阵记为R1(αS)；将坐标系(x’，y’，z)绕x’轴顺时针旋转，旋转角度为αD，得到坐标系(x’，y”，z’)，对应的旋转变换矩阵记为R2(αD)；将坐标系(x’，y”，z’)绕z’轴旋转，旋转角度为αL，得到主轴各向异性坐标系(x”，y”’，z’)，对应的旋转变换矩阵记为R3(αL)。其中，R1(αS)、R2(αD)和R3(αL)分别为：
Figure BDA0002139012770000084
Figure BDA0002139012770000085
Figure BDA0002139012770000086
基于此，可以确定出总的旋转变换矩阵：
R＝R1(αS)R2(αD)R3(αL)····(4)，
而总的旋转变换矩阵R满足：
RRT＝I························(5)，
其中，RT为R的转置矩阵，I为单位矩阵。
由此可知，R是正交矩阵，故RT＝R-1，R-1为R的逆矩阵。
在各向异性介质中，电导率是个大小为3×3阶的张量。在得到了总的旋转变换矩阵之后，便可进一步得到测量坐标系中各物理量与主轴各向异性坐标系中各物理量之间的关系：
J′＝RJ························(6)，
Figure BDA0002139012770000091
Figure BDA0002139012770000092
E′＝RE······················(9)，
其中，J、E和
Figure BDA0002139012770000093
分别为测量坐标系下的电流密度、电场和电导率张量；J′、E′和
Figure BDA0002139012770000094
分别为主轴各向异性坐标系下的电流密度、电场和电导率张量。
进一步，由式(6)-(9)可确定出：
Figure BDA0002139012770000095
其中，
Figure BDA0002139012770000096
从而，确定出
Figure BDA0002139012770000097
为：
Figure BDA0002139012770000098
而
Figure BDA0002139012770000099
由三个电导率σ1、σ2、σ3及三个欧拉旋转角αS、αD、αL共同决定，可以表示任意各向异性。其中每个元素分别为：
Figure BDA00021390127700000910
Figure BDA00021390127700000911
σzz＝sin2αD(sin2αLσ1+cos2αLσ2)+cos2αDσ3········(15)，
Figure BDA0002139012770000101
Figure BDA0002139012770000102
Figure BDA0002139012770000103
而
Figure BDA0002139012770000104
是对称正定的，此外，电阻率张量
Figure BDA0002139012770000105
为
Figure BDA0002139012770000106
的逆：
Figure BDA0002139012770000107
根据麦克斯韦方程组，并忽略位移电流，可以确定：
Figure BDA0002139012770000108
与之相似地，磁导率张量
Figure BDA0002139012770000109
也可由类似方式确定出，此处不再赘述。确定出的
Figure BDA00021390127700001010
和
Figure BDA00021390127700001011
可分别表示如下：
Figure BDA00021390127700001012
Figure BDA00021390127700001013
其中，
Figure BDA00021390127700001014
为电导率张量，σ为电导率；
Figure BDA00021390127700001015
为磁导率张量；μ0为空气中的磁导率，χ为磁化率，
Figure BDA00021390127700001016
为磁阻率张量。
需要说明的是，本申请中采用欧拉旋转，以获取各向异性的
Figure BDA00021390127700001017
和
Figure BDA00021390127700001018
并且以z轴、x轴、z轴的方式旋转，在一些可实现的方式中，也可以通过其他的旋转轴和旋转顺序进行旋转，以获取各向异性的
Figure BDA00021390127700001019
和
Figure BDA00021390127700001020
此处不作为限定。
步骤S40：确定出每个单元的一次场，其中，所述一次场表示排除所述异常体时所述模拟区域的电场。
在本实施例中，关于一次场的值，可以通过计算一次场模型(通常可以为均匀半空间或水平层状介质模型)时确定出来，由于一次场对应的是简单的均匀版空间或者层状介质，可通过模型获取，因此，其计算过程此处不作详细描述。示例性的，一次场电场的计算过程可参考：一维电阻率各向异性对海洋可控源电磁响应的影响研究(罗鸣，李予国，地球物理学报，2015年8月，第58卷第8期)。
需要说明的是，所述一次场电场表示排除所述异常体时所述模拟区域的电场是指：在模拟区域不存在异常体时模拟区域的电场，即在存在异常体区域时，将模拟区域中的异常体区域替换为与周边区域电导率和磁导率一致的区域而形成的模拟区域，此时整个模拟区域为均匀半空间或层状介质模型，而计算此时的模拟区域的电场，即为一次场电场。
步骤S50：根据模拟区域内的一次场电场、
Figure BDA0002139012770000111
和
Figure BDA0002139012770000112
确定出模拟区域内的二次场电场，其中，所述二次场电场为所述异常体产生的电场。
在本实施例中，将对存在电导率和磁导率各向异性时的模拟区域中每个单元的二次场电场的计算过程进行详细的描述。
由于电场双旋度方程中只包含电磁场中的电场，而磁场双旋度方程只包含电磁场中的磁场，而求得电场后就可求得磁场，或求得磁场后也可求得电场。因此，本实施例中，对采用求解电场后确定磁场的方式进行介绍。可以理解，在其他可实现的方式中，也可以类似地采用求解磁场后确定电场的方式，此处不作限定。
在本实施例中，可以根据模拟区域内的一次场电场、
Figure BDA0002139012770000113
和
Figure BDA0002139012770000114
结合麦克斯韦方程组，确定模拟区域的电场积分方程，具体的过程如下所述。
取时谐因子e-iωt，当只存在电性源且外电流密度矢量为JS，且背景层磁导率为μ0，忽略位移电流，根据麦克斯韦方程组，可以得到关于电场的双旋度方程：
Figure BDA0002139012770000115
其中，
Figure BDA0002139012770000116
和
Figure BDA0002139012770000117
分别为电导率张量和磁导率张量，
Figure BDA0002139012770000118
为哈密顿算子，E是电场，i是
Figure BDA0002139012770000119
ω是角频率。
在本实施例中，也可以将总电场E分解为一次场EP和二次场ES：
E＝EP+ES·················(24)，
从而可以得到关于一次场的电场双旋度方程：
Figure BDA00021390127700001110
其中，
Figure BDA0002139012770000121
为计算一次场模型(通常可以为均匀半空间或水平层状介质模型)时的电导率，背景场相对磁导率为1(即磁化率为0)。
而式(23)减去(25)可以得到二次场的电场双旋度方程：
Figure BDA0002139012770000122
其中，I是单元矩阵，
Figure BDA0002139012770000123
为异常体的电阻率。
而关于
Figure BDA0002139012770000124
有：
Figure BDA0002139012770000125
由此，可以确定出二次场的电场双旋度方程。并且，由于采用的是二次场计算，因此，当边界条件距离异常体足够远时，异常体产生的二次场在边界上可以认为等于零。
在本实施例中，根据每个单元的电场双旋度方程及
Figure BDA0002139012770000126
和
Figure BDA0002139012770000127
确定出的方程，即电场双旋度方程，可以点乘电场强度的变分后对整个模拟区域积分，得到中间方程，具体过程如下所述。
可以用二次场Es的变分δES点乘式(26)，并对整个区域积分得到：
Figure BDA0002139012770000128
式中，
δEs＝δEs，x+δEs，y+δEs，z·············(29)，
其中，δEs，x、δEs，y和δEs，z分别表示δEs在x轴、y轴和z轴方向的分量。
确定出中间方程后，可以利用矢量恒等式和散度定理，确定出模拟区域的电场积分方程，具体过程如下所述。
在本实施例中，可以根据矢量恒等式：
Figure BDA0002139012770000129
可将式(28)中的第一项变为：
Figure BDA00021390127700001210
而根据散度定理(或高斯定理)：
Figure BDA0002139012770000131
则可将式(31)的右边第一项变为：
Figure BDA0002139012770000132
式中，U代表外边界面。
当采用第一类边界条件时，边界上的场值为零，同样根据矢量恒等式和散度定理，式(28)最后一项也可简化，可以得到：
Figure BDA0002139012770000133
由此，就可以确定出模拟区域的电场积分方程，即式(34)。
通过这样的方式，即使用伽里金加权余量法，能够快速准确地确定出模拟区域的电场积分方程，提高效率。需要说明的是，在其他一些实现方式中，也可以采用变分方程推导出模拟区域的电场积分方程，因此，此处不应视为对本申请的限定。确定出模拟区域的电场积分方程后，可以对模拟区域中的一个单元进行分析，可以据此确定出模拟区域的总体系数矩阵，具体过程如下所述。
在本实施例中，为了方便描述，将以规则六面体的单元为例，对本实施例的正演方法进行介绍。
如图5所示出的一个六面体单元，x、y和z三个方向的边长分别记为a、b和c，其原点坐标为(x0，y0，z0)。可以将切向量场赋在棱边上，则三个方向上的场分向量可表示为：
Figure BDA0002139012770000134
Figure BDA0002139012770000135
其中，N为Whitney插值基函数；i，j，k分别表示x分量，y分量，z分量的序号，e表示第e个单元。
具体地：
Figure BDA0002139012770000141
Figure BDA0002139012770000142
Figure BDA0002139012770000143
Figure BDA0002139012770000144
Figure BDA0002139012770000145
Figure BDA0002139012770000146
Figure BDA0002139012770000147
Figure BDA0002139012770000148
Figure BDA0002139012770000149
Figure BDA00021390127700001410
Figure BDA00021390127700001411
Figure BDA00021390127700001412
而单元上的E可表示为whitney插值基函数与棱边上场值乘积的和：
Figure BDA00021390127700001413
Figure BDA00021390127700001414
其中，i＝1，2，3，4；带x、y、z脚标的分别表示在x轴、y轴、z轴上的分量。
将整个模拟区域积分分解为各单元积分：
Figure BDA0002139012770000151
以下，将对式(51)中的各项分别处理，此处有：
Figure BDA0002139012770000152
对式(51)的第一项进行处理，在单元内有：
Figure BDA0002139012770000153
对式(51)的第二项进行处理，在单元内有：
Figure BDA0002139012770000154
对式(51)的第三项进行处理，在单元内有：
Figure BDA0002139012770000155
对式(51)的第四项进行处理，在单元内有：
Figure BDA0002139012770000161
其中，
Figure BDA0002139012770000162
式中，I表示3×3的单元矩阵。
在一个单元内，将式(53)减去式(54)，再减去式(55)，最后再减去式(57)，便可确定出一个单元的系数矩阵。
在本实施例中，确定出一个单元内的系数矩阵后，可以将一个单元内的系数矩阵扩展到所有单元，对整个模拟区域中所有单元的系数矩阵进行整合，便得到模拟区域的总体系数矩阵：
Figure BDA0002139012770000163
而由于
Figure BDA0002139012770000164
的任意性，可以确定出：
(K1-K2)ES＝(K3+K4)EP············(59)。
将所有棱边上的电场一次场值EP及边界上的二次场值ES代入式(54)，解之便可以得到每条棱边上的二次场电场值。
具体地，确定出模拟区域的总体系数矩阵后，可以采用无损解精度的直接法加入预设的第一类边界条件，可得到方程组：
Figure BDA0002139012770000171
如果已知φ3＝p3，φ5＝p5，φ6＝p6，那么求解式(60)等效于求解下列方程组：
Figure BDA0002139012770000172
通过直接法添加边界条件，消除了方程中第一类边界条件处的量，降低了方程的维数，尤其是在已知量比较多的情况下，能有效地降低矩阵方程的维度，提升处理速度，从而提高三维频率域可控源电磁的正演方法的运行效率。
加入边界条件后，总体系数矩阵则变为：
KE＝P·············(62)，
其中，K是最后形成的总体系数矩阵，即式(59)中的K1-K2，P则为式(59)中的(K3+K4)Ep。解方程组，便可得到所有单元的每条棱边上的二次场值。此处，可以使用对称超松弛迭代法预处理，并采用稳定双共轭梯度法求解器进行解方程，求得模拟区域内所有单元上的二次场值。
通过采用从一个单元进行分析后扩散到整个模拟区域的方式确定出模拟区域的总体系数矩阵，而矢量有限元法的应用，可以允许二次场法向分量在电性不连续界面处突变，且考虑了电导率和磁导率的各向异性，更符合实际情况，从而使应用此正演方法的正演系统具有更加广泛的应用范围和更高的准确性。
需要说明的是，在本实施例中，根据一次场电场、
Figure BDA0002139012770000173
和
Figure BDA0002139012770000174
确定二次场电场的过程，可以为本实施例中介绍的方式。在实际的计算过程中，由一次场电场、
Figure BDA0002139012770000175
和
Figure BDA0002139012770000176
确定二次场电场、再确定出总体系数矩阵的过程，不需要在中间过程中计算出二次场的值、中间方程的值等，只需在最后过程中将各个量的实际值代入最后的总体系数矩阵计算出所有单元的二次场值即可。这样在最后计算，从理论上说，可以避免由中间过程的计算带来的误差，从而提高结果的准确性。
步骤S60：根据模拟区域内的一次场电场和二次场电场，确定出所述模拟区域的电场和磁场。
在本实施例中，可以结合模拟区域内的一次场电场值和二次场电场值，求得每个单元的总电场值。而确定出每个单元的总电场值后，可以确定出模拟区域内每个单元的磁场H。具体地：
Figure BDA0002139012770000181
Figure BDA0002139012770000182
Figure BDA0002139012770000183
确定出模拟区域内每个单元的电场和磁场后，就可得到模拟区域的电场和磁场了。
步骤S70：根据所述模拟区域的电场和磁场，确定出所述模拟区域内地表出的使电阻率及相位。
在本实施例中，可以基于模拟区域内的电场和磁场，确定出地表处的张量阻抗Z：
Figure BDA0002139012770000184
Figure BDA0002139012770000185
Figure BDA0002139012770000186
Figure BDA0002139012770000187
其中，Zxx、Zxy、Zyx、Zyy分别为：xx、xy、yx、yy模式下的张量阻抗。
根据每个单元的张量阻抗，可以确定出模拟区域内地表处的视电阻率ρij及相位
Figure BDA0002139012770000188
Figure BDA0002139012770000191
Figure BDA0002139012770000192
(i＝x，y；j＝x，y)
在本实施例中，模拟区域内地表处的视电阻率及相位，就可以表征异常体的电磁响应结果，即三维频率域可控源电磁对异常体区域的正演结果。其中，地表指代模拟区域中空气与地下介质交界处的面。
应用本申请提供的三维频率域可控源电磁的正演方法，通过在三维频率域可控源电磁的正演方法中考虑电导率和磁导率的各向异性，可以定量地分析电导率与磁导率各向异性对三维频率域可控源电磁响应的影响，这有利于更加准确地了解和分析实际中复杂的地质结构的电磁响应，因此适用于更多类型的地质结构，应用领域也更加广阔。并且，本方法也可以为三维频率域可控源电磁电导率与磁导率各向异性反演提供基础，能够进一步扩展本方法的应用领域。
为了尽可能全面合理地对本申请实施例提供的三维频率域可控源电磁的正演方法进行验证，本实施例中，对一个电阻率任意各向异性模型进行了计算，并将其与有限差分的结果进行对比；另外，还计算了一个磁导率层状介质模型，并与解析解进行了对比。具体如下：
电阻率任意各向异性模型：如图6所示，在一个电阻率为100Ω·m磁化率为0的均匀半空间中有一个三维异常体，其尺寸为350m×350m×350m，埋深为150m，三个主轴电阻率分别为50/25/20Ω·m，三个欧拉角分别为20°/45°/30°。计算频点为100Hz。张量源中心坐标为(0m，-8000m，1m)，单个源长度为300m，电流为1A。其中，图6(a)为电阻率任意各向异性模型的侧视图，图6(b)为电阻率任意各向异性模型的平面图。
需要说明的是，在本实施例中，除了特殊说明外，采用的都是此种张量源，但并不限定于此。本实施例中所用的张量源限定为同一种，只是为了方便说明，在其他一些实施例中，也可以采用其他的张量源。
经计算，比较可控源电磁的视电阻率和大地电磁的视电阻率，最大相对误差小于0.5％，二者结果吻合很好。
请参阅图7，本实例采用对称超松弛迭代法预处理，并采用准最小余量迭代求解器进行解方程，以求得方程的解。收敛情况如图7所示，可以看到二者收敛较好。
图8所示为本实例中有限元的结果和有限差分的结果的对比，第一行对应xy模式的视电阻率，第二行对应yx模式的视电阻率，第一列对应有限元的两种模式的视电阻率，第二列对应有限差分的两种模式的视电阻率，第三列对应相对误差，可以看到有限元和有限差分的结果较一致，最大相对误差小于0.5％。
磁导率层状介质模型：磁导率模型如图9所示，在100Ω·m的均匀半空间中有一个磁化率为2的地层。本实施例中计算三个频点(100Hz、40Hz和20Hz)，数值解和解析解对比如图10所示，其中第二列表示解析解，第三列表示数值解，可以看出数值解和解析解吻合很好，最大相对误差小于0.5％。
通过以上两个模型的计算对比和验证，证明了该正演方法的准确性。而将本申请实施例提供的三维频率域可控源电磁的正演方法的计算结果分别与电导率模型和磁导率模型进行对比，可以全面地验证本方法的准确性。
而本申请实施例提供的三维频率域可控源电磁的正演方法，可以用于验证多种电导率和磁导率各向异性的异常体介质的电磁响应，以及定性研究电导率和磁导率对介质的电磁响应的影响，以便更好地了解各种介质的电磁响应。而基于本正演方法建立的反演模型，也能够适用于更多类型的地质，以及，提供更准确的反演结果，对地质勘测起到积极的促进作用。
请参阅图11，基于同一发明构思，本申请实施例中还提供一种三维频率域可控源电磁的正演系统11。包括：
区域确定单元11，用于确定模拟区域，所述模拟区域包括表征异常体的异常体区域及与所述异常体区域互补的周边区域，所述异常体区域内的电导率和磁导率不同于所述周边区域的电导率和磁导率；区域剖分单元12，用于将所述模拟区域剖分为多个单元；张量获取单元13，用于获取每个单元的电导率张量
Figure BDA0002139012770000201
和磁导率张量
Figure BDA0002139012770000202
一次场确定单元14，用于确定出区域内的一次场电场，其中，所述一次场电场表示排除所述异常体时所述模拟区域的电场；二次场确定单元15，用于根据模拟区域内的一次场电场、
Figure BDA0002139012770000203
知
Figure BDA0002139012770000204
确定出模拟区域内的二次场电场，其中，所述二次场电场为所述异常体产生的电场；正演处理单元16，用于根据模拟区域内的一次场电场和二次场电场，确定出所述模拟区域的电场和磁场；以及，根据所述模拟区域的电场和磁场，确定出所述模拟区域内地表处的视电阻率及相位。
在本实施例中，所述区域剖分单元12，还用于采用稀疏网格将所述异常体区域均匀剖分，以及，沿远离所述异常体的方向采用逐步增大的网格对所述周边区域进行剖分，以将所述模拟区域剖分为多个单元。
在本实施例中，所述张量获取单元13，还用于获取所述模拟区域内每个单元的主轴各向异性坐标系中x”主轴方向的电导率σ1、y”’主轴方向的电导率σ2、z’主轴方向的电导率σ3，并获取三次欧拉旋转过程中每次旋转的角度αS、αD和αL，以及，获取所述模拟区域内每个单元的主轴各向异性坐标系中x”主轴方向的磁化率χ1、y”’主轴方向的磁化率χ2、z’主轴方向的磁化率χ3，并获取三次欧拉旋转过程中每次旋转的角度βS、βD和βL；根据每个单元的σ1、σ2、σ3、αS、αD和αL，确定出每个单元的
Figure BDA0002139012770000211
以及，根据每个单元的χ1、χ2、χ3、βS、βD和βL，确定出每个单元的
Figure BDA0002139012770000212
在本实施例中，所述二次场确定单元15，还用于根据模拟区域内的一次场电场、
Figure BDA0002139012770000213
和
Figure BDA0002139012770000214
结合麦克斯韦方程组，确定所述模拟区域的电场积分方程；根据所述电场积分方程，确定出所述模拟区域的总体系数矩阵；根据预设的边界条件和所述总体系数矩阵，确定出所述模拟区域内的二次场电场。
在本实施例中，所述二次场确定单元15，还用于根据模拟区域内的一次场电场，确定出每个单元的电场双旋度方程；由根据每个单元的电场双旋度方程及
Figure BDA0002139012770000215
和
Figure BDA0002139012770000216
确定出的方程点乘电场强度的变分后对整个模拟区域积分，得到中间方程；利用矢量恒等式和散度定理，基于所述中间方程确定出所述模拟区域的电场积分方程。
在本实施例中，在每个单元为六面体时，电场赋在棱边上，所述二次场确定单元15，还用于根据所述积分方程，确定出每个单元的单元方程；根据每个单元的单元方程中包含的系数矩阵，确定出所述模拟区域的总体系数矩阵。
在本实施例中，所述正演处理单元16，还用于基于模拟区域内的一次场电场和二次场电场，确定出模拟区域内的电场；每条棱边上的磁场可根据该棱边及其邻近棱边上的电场确定出，其中，所述模拟区域内所有单元的电场和磁场的总体表示所述模拟区域的电场和磁场。
如图12所示，本申请的实施例还提供一种电子设备20，本申请实施例提供的正演方法，可以由电子设备20执行；而应用本申请提供的正演方法的正演系统，可以基于电子设备20实现其功能。电子设备20可以为终端，例如智能手机、平板电脑、个人电脑、个人数字助理等；电子设备20也可以为服务器，例如网络服务器、云服务器、服务器集群、数据服务器等，此处不作限定。
请参阅图12，电子设备20可以包括：通过网络与外界连接的通信模块22、用于执行程序指令的一个或多个处理器24、总线23和不同形式的存储器21，例如，磁盘、ROM、或RAM，或其任意组合。
示例性的，存储器21中存储有程序。处理器24可以从存储器21调用并运行这些程序，从而便可以通过运行程序而执行大地电磁的正演方法。处理器24通过对三维频率域可控源电磁的正演方法的执行，处理器24可以对多种不同的地质类型进行正演，以完成对三维频率域可控源电磁的正演方法的执行。
本申请的实施例还提供一种存储介质，存储介质中包含存储的程序，其中，在该程序运行时控制存储介质所在设备执行本实施例的三维频率域可控源电磁的正演方法。
综上所述，本申请的实施例提供一种三维频率域可控源电磁的正演方法、系统、存储介质及电子设备，通过在三维频率域可控源电磁的正演方法中考虑电导率和磁导率的各向异性，可以定性地分析电导率与磁导率各向异性对三维频率域可控源电磁响应的影响，这有利于更加准确地了解和分析实际中复杂的地质结构的电磁响应。并且，本方法也可以为三维频率域可控源电磁电导率与磁导率各向异性反演提供基础。以及，由于本方法同时考虑了电导率和磁导率的各向异性，适用于更多类型的地质结构，因此，可以在更复杂的环境中应用。
在本申请所提供的实施例中，应该理解到，所揭露装置和方法，可以通过其它的方式实现。以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式。
另外，作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。
以上所述仅为本申请的实施例而已，并不用于限制本申请的保护范围，对于本领域的技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。

Claims (9)

1.一种三维频率域可控源电磁的正演方法，其特征在于，包括：
确定模拟区域，所述模拟区域包括表征异常体的异常体区域及与所述异常体区域互补的周边区域，所述异常体区域内的电导率和磁导率不同于所述周边区域的电导率和磁导率；
将所述模拟区域剖分为多个单元；
获取每个单元的电导率张量
Figure FDA0002606560170000011
和磁导率张量
Figure FDA0002606560170000012
确定出模拟区域内的一次场电场，其中，所述一次场表示排除所述异常体时所述模拟区域的电场；
根据模拟区域内的一次场电场、
Figure FDA0002606560170000013
和
Figure FDA0002606560170000014
确定出模拟区域内的二次场电场，其中，所述二次场电场为所述异常体产生的电场；
根据模拟区域内的一次场电场和二次场电场，确定出所述模拟区域的电场和磁场；
根据所述模拟区域的电场和磁场，确定出所述模拟区域内地表处的视电阻率及相位；
其中，所述获取每个单元的电导率张量
Figure FDA0002606560170000015
和磁导率张量
Figure FDA0002606560170000016
包括：
获取所述模拟区域内每个单元的主轴各向异性坐标系中x”主轴方向的电导率σ1、y”’主轴方向的电导率σ2、z’主轴方向的电导率σ3，并获取三次欧拉旋转过程中每次旋转的角度αS、αD和αL，以及，获取所述模拟区域内每个单元的主轴各向异性坐标系中x”主轴方向的磁化率χ1、y”’主轴方向的磁化率χ2、z’主轴方向的磁化率χ3，并获取三次欧拉旋转过程中每次旋转的角度βS、βD和βL；根据每个单元的σ1、σ2、σ3、αS、αD和αL，确定出每个单元的
Figure FDA0002606560170000017
以及，根据每个单元的χ1、χ2、χ3、βS、βD和βL，确定出每个单元的
Figure FDA0002606560170000018
确定出的
Figure FDA0002606560170000019
及
Figure FDA00026065601700000110
分别为：
Figure FDA0002606560170000021
Figure FDA0002606560170000022
其中，
Figure FDA0002606560170000023
为电导率张量，σ为电导率；
Figure FDA0002606560170000024
为磁导率张量；μ0为空气中的磁导率，χ为磁化率，
Figure FDA0002606560170000025
为磁阻率张量。
2.根据权利要求1所述的正演方法，其特征在于，所述将所述模拟区域剖分为多个单元，包括：
采用稀疏网格将所述异常体区域均匀剖分，以及，沿远离所述异常体的方向采用逐步增大的网格对所述周边区域进行剖分，以将所述模拟区域剖分为多个单元。
3.根据权利要求1所述的正演方法，其特征在于，所述根据模拟区域内的一次场电场、
Figure FDA0002606560170000026
和
Figure FDA0002606560170000027
确定出模拟区域内的二次场电场，包括：
根据模拟区域内的一次场电场、
Figure FDA0002606560170000028
和
Figure FDA0002606560170000029
结合麦克斯韦方程组，确定所述模拟区域的电场积分方程；
根据所述电场积分方程，确定出所述模拟区域的总体系数矩阵；
根据预设的边界条件和所述总体系数矩阵，确定出所述模拟区域内的二次场电场。
4.根据权利要求3所述的正演方法，其特征在于，所述根据模拟区域内的一次场电场、
Figure FDA00026065601700000210
和
Figure FDA00026065601700000211
结合麦克斯韦方程组，确定所述模拟区域的电场积分方程，包括：
根据模拟区域内的一次场，确定出模拟区域内的电场双旋度方程；
由根据每个单元的电场双旋度方程及
Figure FDA00026065601700000212
和
Figure FDA00026065601700000213
确定出的方程点乘电场强度的变分后对整个模拟区域积分，得到中间方程；
利用矢量恒等式和散度定理，基于所述中间方程确定出所述模拟区域的电场积分方程。
5.根据权利要求3所述的正演方法，其特征在于，在每个单元为六面体时，电场赋在棱边上，所述根据所述电场积分方程，通过有限元分析，确定出所述模拟区域的总体系数矩阵，包括：
根据所述电场积分方程，确定出每个单元的单元方程；
根据每个单元的单元方程中包含的系数矩阵，确定出所述模拟区域的总体系数矩阵。
6.根据权利要求1所述的正演方法，其特征在于，所述根据模拟区域内的一次场电场和二次场电场，确定出所述模拟区域的电场和磁场，包括：
基于模拟区域内的一次场电场和二次场电场，确定出模拟区域内的电场；
每条棱边上的磁场可根据该棱边及其邻近棱边上的电场确定出，其中，所述模拟区域内所有单元的电场和磁场的总体表示所述模拟区域的电场和磁场。
7.一种三维频率域可控源电磁的正演系统，其特征在于，包括：
区域确定单元，用于确定模拟区域，所述模拟区域包括表征异常体的异常体区域及与所述异常体区域互补的周边区域，所述异常体区域内的电导率和磁导率不同于所述周边区域的电导率和磁导率；
区域剖分单元，用于将所述模拟区域剖分为多个单元；
张量获取单元，用于获取每个单元的电导率张量
Figure FDA0002606560170000031
和磁导率张量
Figure FDA0002606560170000032
一次场确定单元，用于确定出模拟区域内的一次场电场，其中，所述一次场表示排除所述异常体时所述模拟区域的电场；
二次场确定单元，用于根据模拟区域内的一次场电场、
Figure FDA0002606560170000033
和
Figure FDA0002606560170000034
确定出模拟区域内的二次场电场，其中，所述二次场电场为所述异常体产生的电场；
正演处理单元，用于根据每个单元的一次场电场和二次场电场，确定出所述模拟区域的电场和磁场；以及，根据所述模拟区域的电场和磁场，确定出所述模拟区域内地表处的视电阻率及相位；
其中，所述张量获取单元，还用于获取所述模拟区域内每个单元的主轴各向异性坐标系中x”主轴方向的电导率σ1、y”’主轴方向的电导率σ2、z’主轴方向的电导率σ3，并获取三次欧拉旋转过程中每次旋转的角度αS、αD和αL，以及，获取所述模拟区域内每个单元的主轴各向异性坐标系中x”主轴方向的磁化率χ1、y”’主轴方向的磁化率χ2、z’主轴方向的磁化率χ3，并获取三次欧拉旋转过程中每次旋转的角度βS、βD和βL；根据每个单元的σ1、σ2、σ3、αS、αD和αL，确定出每个单元的
Figure FDA0002606560170000041
以及，根据每个单元的χ1、χ2、χ3、βS、βD和βL，确定出每个单元的
Figure FDA0002606560170000042
确定出的
Figure FDA0002606560170000043
及
Figure FDA0002606560170000044
分别为：
Figure FDA0002606560170000045
Figure FDA0002606560170000046
其中，
Figure FDA0002606560170000047
为电导率张量，σ为电导率；
Figure FDA0002606560170000048
为磁导率张量；μ0为空气中的磁导率，χ为磁化率，
Figure FDA0002606560170000049
为磁阻率张量。
8.一种存储介质，其特征在于，所述存储介质包括存储的程序，其中，在所述程序运行时控制所述存储介质所在设备执行权利要求1至6中任意一项所述的三维频率域可控源电磁的正演方法。
9.一种电子设备，包括存储器和处理器，所述存储器用于存储包括程序指令的信息，所述处理器用于控制程序指令的执行，其特征在于：所述程序指令被处理器加载并执行时实现权利要求1至6任意一项所述的三维频率域可控源电磁的正演方法的步骤。
CN201910663566.5A 2019-07-22 2019-07-22 三维频率域可控源电磁的正演方法、系统 Expired - Fee Related CN110346834B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910663566.5A CN110346834B (zh) 2019-07-22 2019-07-22 三维频率域可控源电磁的正演方法、系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910663566.5A CN110346834B (zh) 2019-07-22 2019-07-22 三维频率域可控源电磁的正演方法、系统

Publications (2)

Publication Number Publication Date
CN110346834A CN110346834A (zh) 2019-10-18
CN110346834B true CN110346834B (zh) 2020-11-17

Family

ID=68179684

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910663566.5A Expired - Fee Related CN110346834B (zh) 2019-07-22 2019-07-22 三维频率域可控源电磁的正演方法、系统

Country Status (1)

Country Link
CN (1) CN110346834B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111812732B (zh) * 2020-06-29 2024-03-15 中铁二院工程集团有限责任公司 一种基于卷积神经网络的大地电磁非线性反演方法
CN113109876B (zh) * 2021-02-22 2022-06-10 长江勘测规划设计研究有限责任公司 基于三维磁电阻率法张量观测模式的解译方法及装置
CN113779818B (zh) * 2021-11-15 2022-02-08 中南大学 三维地质体其电磁场数值模拟方法、装置、设备及介质
CN116341332A (zh) * 2023-03-30 2023-06-27 重庆大学 基于电导率分块连续变化的大地电磁三维有限元正演方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7894989B2 (en) * 2005-06-09 2011-02-22 Exxonmobil Upstream Research Co. Method for determining earth vertical electrical anisotropy in marine electromagnetic surveys
CN106547031A (zh) * 2016-10-09 2017-03-29 煤炭科学技术研究院有限公司 一种基于场变换确定地质电性界面的瞬变电磁探测法
CN106842343A (zh) * 2017-02-14 2017-06-13 中国科学院地质与地球物理研究所 一种电性源瞬变电磁电场响应成像方法
CN106980736A (zh) * 2017-04-11 2017-07-25 吉林大学 一种各向异性介质的海洋可控源电磁法有限元正演方法
CN109407160A (zh) * 2018-12-19 2019-03-01 中国科学院地质与地球物理研究所 一种基于wem矢量有限元法进行三维正演的方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7894989B2 (en) * 2005-06-09 2011-02-22 Exxonmobil Upstream Research Co. Method for determining earth vertical electrical anisotropy in marine electromagnetic surveys
CN106547031A (zh) * 2016-10-09 2017-03-29 煤炭科学技术研究院有限公司 一种基于场变换确定地质电性界面的瞬变电磁探测法
CN106842343A (zh) * 2017-02-14 2017-06-13 中国科学院地质与地球物理研究所 一种电性源瞬变电磁电场响应成像方法
CN106980736A (zh) * 2017-04-11 2017-07-25 吉林大学 一种各向异性介质的海洋可控源电磁法有限元正演方法
CN109407160A (zh) * 2018-12-19 2019-03-01 中国科学院地质与地球物理研究所 一种基于wem矢量有限元法进行三维正演的方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
二维大地电磁各向异性参数对视电阻率的影响研究;蔡义宇 等;《地球物理学进展》;20190308;1-11 *

Also Published As

Publication number Publication date
CN110346834A (zh) 2019-10-18

Similar Documents

Publication Publication Date Title
CN110346834B (zh) 三维频率域可控源电磁的正演方法、系统
CN110346835B (zh) 大地电磁的正演方法、正演系统、存储介质及电子设备
Li et al. A finite-element time-domain forward solver for electromagnetic methods with complex-shaped loop sources
CN113221393B (zh) 一种基于非结构有限元法的三维大地电磁各向异性反演方法
Ren et al. 3D direct current resistivity modeling with unstructured mesh by adaptive finite-element method
Jahandari et al. A finite-volume solution to the geophysical electromagnetic forward problem using unstructured grids
CN110058315B (zh) 一种三维各向异性射频大地电磁自适应有限元正演方法
Čuma et al. Large‐scale 3D inversion of potential field data
Wu et al. High-precision Fourier forward modeling of potential fields
Schaa et al. PDE-based geophysical modelling using finite elements: examples from 3D resistivity and 2D magnetotellurics
CN113962077B (zh) 三维各向异性强磁场数值模拟方法、装置、设备及介质
CN116842813B (zh) 一种三维大地电磁正演数值模拟方法
CN106019394A (zh) 海洋大地电磁场非线性共轭梯度三维并行反演方法
Ren et al. Recursive analytical formulae of gravitational fields and gradient tensors for polyhedral bodies with polynomial density contrasts of arbitrary non-negative integer orders
Ouyang et al. Iterative magnetic forward modeling for high susceptibility based on integral equation and Gauss-fast Fourier transform
Yu et al. Two-dimensional anisotropic magnetotelluric inversion using a limited-memory quasi-Newton method
CN114970289B (zh) 三维大地电磁各向异性正演数值模拟方法、设备及介质
Yin et al. A finite volume method for global electromagnetic induction forward modeling on collocated unstructured grids
Xia et al. Numerical investigations of blockiness of fractured rocks based on fracture spacing and disc diameter
Muravyev et al. Approximation of anomalous magnetic field using array of singular sources (magnetized rods)
Jia et al. Modeling of complex geological body and computation of geomagnetic anomaly
Ivanov et al. Three-dimensional inversion of the single-profile magnetotelluric data
Liu et al. Three-Dimensional Inversion of Time-Domain Electromagnetic Data Using Various Loop Source Configurations
Prihantoro et al. Numerical solution of 3-D magnetotelluric using vector finite element method
Bomba et al. Identifying the structure of soil massifs by numerical quasiconformal mapping methods

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
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20201117

Termination date: 20210722