三维频率域可控源电磁的正演方法、系统与流程

本申请涉及勘探地球物理领域，具体而言，涉及一种三维频率域可控源电磁的正演方法、系统、存储介质及电子设备。

背景技术：

频率域可控源电磁法是一类非常重要的地球物理方法，具有横向分辨能力强、工作效率高及信噪比高等优点，目前被广泛应用于矿产资源勘查、油气勘探及工程勘察等领域。在实际的勘查中，通常对观测数据进行数据模拟，以确定地下介质的电性结构。但目前的正演方法及正演系统的应用范围有限，而基于现有的正演系统对具有复杂电磁响应的地下介质进行的反演，难以有效地应用在更复杂的环境中。

技术实现要素：

本申请实施例的目的在于提供一种三维频率域可控源电磁的正演方法、系统、存储介质及电子设备，以使基于频率域可控源电磁法的正演方法和正演系统应用于更复杂的环境中。

为了实现上述目的，本申请的实施例通过如下方式实现：

第一方面，本申请的实施例提供一种三维频率域可控源电磁的正演方法，包括：确定模拟区域，所述模拟区域包括表征异常体的异常体区域及与所述异常体区域互补的周边区域，所述异常体区域内的电导率和磁导率不同于所述周边区域的电导率和磁导率；将所述模拟区域剖分为多个单元；获取每个单元的电导率张量和磁导率张量确定出模拟区域内的一次场电场，其中，所述一次场电场表示排除所述异常体时所述模拟区域的电场；根据每个单元的一次场电场、和确定出模拟区域内的二次场电场，其中，所述二次场电场为所述异常体产生的电场；根据模拟区域内的一次场电场和二次场电场，确定出所述模拟区域的电场和磁场；根据所述模拟区域的电场和磁场，确定出所述模拟区域内地表处的视电阻率及相位。

通过在三维频率域可控源电磁的正演方法中考虑电导率和磁导率的各向异性，可以定量地分析电导率与磁导率各向异性对三维频率域可控源电磁响应的影响，这有利于更加准确地了解和分析实际中复杂的地质结构的电磁响应。并且，本方法也可以为三维频率域可控源电磁电导率与磁导率各向异性反演提供基础。以及，由于本方法同时考虑了电导率和磁导率的各向异性，适用于更多类型的地质结构，因此，可以在更复杂的环境中应用。

结合第一方面，在第一方面的第一种可能的实现方式中，所述将所述模拟区域剖分为多个单元，包括：采用稀疏网格将所述异常体区域均匀剖分，以及，沿远离所述异常体的方向采用逐步增大的网格对所述周边区域进行剖分，以将所述模拟区域剖分为多个单元。

通过这样的方式对模拟区域进行剖分，能尽可能减小异常体区域在模拟区域的外边界处的影响，甚至可以忽略不计。另外，采用稀疏网格配置方法，可以降低矩阵方程维度，节约求解过程花费的时间，从而提高效率。

结合第一方面，在第一方面的第二种可能的实现方式中，所述获取每个单元的电导率张量和磁导率张量包括：获取所述模拟区域内每个单元的主轴各向异性坐标系中x”主轴方向的电导率σ1、y”’主轴方向的电导率σ2、z’主轴方向的电导率σ3，并获取三次欧拉旋转过程中每次旋转的角度αs、αd和αl，以及，获取所述模拟区域内每个单元的主轴各向异性坐标系中x”主轴方向的磁化率χ1、y”’主轴方向的磁化率χ2、z’主轴方向的磁化率χ3，并获取三次欧拉旋转过程中每次旋转的角度βs、βd和βl；根据每个单元的σ1、σ2、σ3、αs、αd和αl，确定出每个单元的以及，根据每个单元的χ1、χ2、χ3、βs、βd和βl，确定出每个单元的

主轴各向异性坐标系可以体现电导率的各向异性和磁导率的各向异性，是实现将电导率和磁导率的各向异性考虑进正演方法中的基础。

结合第一方面，在第一方面的第三种可能的实现方式中，所述根据模拟区域内的一次场电场、和确定出模拟区域内的二次场电场，包括：根据模拟区域内的一次场电场、和结合麦克斯韦方程组，确定所述模拟区域的电场积分方程；根据所述电场积分方程，确定出所述模拟区域的总体系数矩阵；根据预设的边界条件和所述总体系数矩阵，确定出所述模拟区域内的二次场电场。

采用这种方式可以实现基于模拟区域内的一次场电场、和确定出模拟区域内的二次场电场，使得二次场电场中能够体现异常区域电导率和磁导率的各向异性。

结合第一方面的第三种实现方式，在第一方面的第四种可能的实现方式中，所述根据模拟区域内的一次场电场、和结合麦克斯韦方程组，确定所述模拟区域的电场积分方程，包括：根据模拟区域内的一次场电场，确定出每个单元的电场双旋度方程；由根据每个单元的电场双旋度方程及和确定出的方程点乘电场强度的变分后对整个模拟区域积分，得到中间方程；利用矢量恒等式和散度定理，基于所述中间方程确定出所述模拟区域的电场积分方程。

采用这种方式可以准确地得到对模拟区域积分的电场积分方程。

结合第一方面的第三种实现方式，在第一方面的第五种可能的实现方式中，在每个单元为六面体时，电场赋在棱边上，所述根据所述电场积分方程，确定出所述模拟区域的总体系数矩阵，包括：根据所述电场积分方程，确定出每个单元的单元方程；根据每个单元的单元方程中包含的系数矩阵，确定出所述模拟区域的总体系数矩阵。

采用这种方式，可以允许电场法向分量在电性不连续界面处突变，且考虑了电导率和磁导率的各向异性，更符合实际情况，从而使应用此正演方法具有更加广泛的应用范围和更高的准确性。

结合第一方面，在第一方面的第六种可能的实现方式中，所述根据模拟区域内的一次场电场和二次场电场，确定出所述模拟区域的电场和磁场，包括：基于模拟区域内的一次场电场和二次场电场，确定出模拟区域内的电场；每条棱边上的磁场可根据该棱边及其邻近棱边上的电场确定出，其中，所述模拟区域内所有单元的电场和磁场的总体表示所述模拟区域的电场和磁场。

第二方面，本申请的实施例提供一种三维频率域可控源电磁的正演系统，包括：区域确定单元，用于确定模拟区域，所述模拟区域包括表征异常体的异常体区域及与所述异常体区域互补的周边区域，所述异常体区域内的电导率和磁导率不同于所述周边区域的电导率和磁导率；区域剖分单元，用于将所述模拟区域剖分为多个单元；张量获取单元，用于获取每个单元的电导率张量和磁导率张量一次场确定单元，用于确定出模拟区域内的一次场电场，其中，所述一次场电场表示排除所述异常体时所述模拟区域的电场；二次场确定单元，用于根据模拟区域内的一次场电场、和确定出模拟区域内的二次场电场，其中，所述二次场电场为所述异常体产生的电场；正演处理单元，用于根据模拟取悦的一次场电场和二次场电场，确定出所述模拟区域的电场和磁场；以及，根据所述模拟区域的电场和磁场，确定出所述模拟区域内地表处的视电阻率及相位。

第三方面，本申请的实施例提供一种存储介质，所述存储介质包括存储的程序，其中，在所述程序运行时控制所述存储介质所在设备执行第一方面或第一方面的第一种至第六种可能的实现方式中任一所述的三维频率域可控源电磁的正演方法。

第四方面，本申请的实施例提供一种电子设备，包括存储器和处理器，所述存储器用于存储包括程序指令的信息，所述处理器用于控制程序指令的执行，所述程序指令被处理器加载并执行时实现第一方面或第一方面的第一种至第六种可能的实现方式中任一所述的三维频率域可控源电磁的正演方法的步骤。

为使本申请的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。

附图说明

为了更清楚地说明本申请实施例的技术方案，下面将对本申请实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本申请的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图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：获取每个单元的电导率张量和磁导率张量

在本实施例中，对于电导率和磁导率各向异性的异常体，其电导率张量磁导率张量都可以由主轴各向异性通过欧拉旋转得到。

示例性的，以图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，确定出每个单元的电导率张量以及，根据每个单元的χ1、χ2、χ3、βs、βd和βl，确定出每个单元的磁导率张量

具体地，在本实施例中，根据每个单元的σ1、σ2、σ3、αs、αd和αl，确定出每个单元的的方式，可以是通过每次旋转的旋转变换矩阵。其中，将测量坐标系沿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)分别为：

基于此，可以确定出总的旋转变换矩阵：

r＝r1(αs)r2(αd)r3(αl)····(4)，

而总的旋转变换矩阵r满足：

rr^t＝i························(5)，

其中，r^t为r的转置矩阵，i为单位矩阵。

由此可知，r是正交矩阵，故r^t＝r^-1，r^-1为r的逆矩阵。

在各向异性介质中，电导率是个大小为3×3阶的张量。在得到了总的旋转变换矩阵之后，便可进一步得到测量坐标系中各物理量与主轴各向异性坐标系中各物理量之间的关系：

j′＝rj························(6)，

e′＝re······················(9)，

其中，j、e和分别为测量坐标系下的电流密度、电场和电导率张量；j′、e′和分别为主轴各向异性坐标系下的电流密度、电场和电导率张量。

进一步，由式(6)-(9)可确定出：

其中，

从而，确定出为：

而由三个电导率σ1、σ2、σ3及三个欧拉旋转角αs、αd、αl共同决定，可以表示任意各向异性。其中每个元素分别为：

σzz＝sin²αd(sin²αlσ1+cos²αlσ2)+cos²αdσ3···········(15)，

而是对称正定的，此外，电阻率张量为的逆：

根据麦克斯韦方程组，并忽略位移电流，可以确定：

与之相似地，磁导率张量也可由类似方式确定出，此处不再赘述。确定出的和可分别表示如下：

其中，为电导率张量，σ为电导率；为磁导率张量；μ0为空气中的磁导率，χ为磁化率，为磁阻率张量。

需要说明的是，本申请中采用欧拉旋转，以获取各向异性的和并且以z轴、x轴、z轴的方式旋转，在一些可实现的方式中，也可以通过其他的旋转轴和旋转顺序进行旋转，以获取各向异性的和此处不作为限定。

步骤s40：确定出每个单元的一次场，其中，所述一次场表示排除所述异常体时所述模拟区域的电场。

在本实施例中，关于一次场的值，可以通过计算一次场模型(通常可以为均匀半空间或水平层状介质模型)时确定出来，由于一次场对应的是简单的均匀版空间或者层状介质，可通过模型获取，因此，其计算过程此处不作详细描述。示例性的，一次场电场的计算过程可参考：一维电阻率各向异性对海洋可控源电磁响应的影响研究(罗鸣，李予国，地球物理学报，2015年8月，第58卷第8期)。

需要说明的是，所述一次场电场表示排除所述异常体时所述模拟区域的电场是指：在模拟区域不存在异常体时模拟区域的电场，即在存在异常体区域时，将模拟区域中的异常体区域替换为与周边区域电导率和磁导率一致的区域而形成的模拟区域，此时整个模拟区域为均匀半空间或层状介质模型，而计算此时的模拟区域的电场，即为一次场电场。

步骤s50：根据模拟区域内的一次场电场、和确定出模拟区域内的二次场电场，其中，所述二次场电场为所述异常体产生的电场。

在本实施例中，将对存在电导率和磁导率各向异性时的模拟区域中每个单元的二次场电场的计算过程进行详细的描述。

由于电场双旋度方程中只包含电磁场中的电场，而磁场双旋度方程只包含电磁场中的磁场，而求得电场后就可求得磁场，或求得磁场后也可求得电场。因此，本实施例中，对采用求解电场后确定磁场的方式进行介绍。可以理解，在其他可实现的方式中，也可以类似地采用求解磁场后确定电场的方式，此处不作限定。

在本实施例中，可以根据模拟区域内的一次场电场、和结合麦克斯韦方程组，确定模拟区域的电场积分方程，具体的过程如下所述。

取时谐因子e^-iωt，当只存在电性源且外电流密度矢量为js，且背景层磁导率为μ0，忽略位移电流，根据麦克斯韦方程组，可以得到关于电场的双旋度方程：

其中，和分别为电导率张量和磁导率张量，为哈密顿算子，e是电场，i是ω是角频率。

在本实施例中，也可以将总电场e分解为一次场ep和二次场es：

e＝ep+es·················(24)，

从而可以得到关于一次场的电场双旋度方程：

其中，为计算一次场模型(通常可以为均匀半空间或水平层状介质模型)时的电导率，背景场相对磁导率为1(即磁化率为0)。

而式(23)减去(25)可以得到二次场的电场双旋度方程：

其中，i是单元矩阵，为异常体的电阻率。

而关于有：

由此，可以确定出二次场的电场双旋度方程。并且，由于采用的是二次场计算，因此，当边界条件距离异常体足够远时，异常体产生的二次场在边界上可以认为等于零。

在本实施例中，根据每个单元的电场双旋度方程及和确定出的方程，即电场双旋度方程，可以点乘电场强度的变分后对整个模拟区域积分，得到中间方程，具体过程如下所述。

可以用二次场es的变分δes点乘式(26)，并对整个区域积分得到：

式中，

δes＝δes，x+δes，y+δes，z·············(29)，

其中，δes，x、δes，y和δes，z分别表示δes在x轴、y轴和z轴方向的分量。

确定出中间方程后，可以利用矢量恒等式和散度定理，确定出模拟区域的电场积分方程，具体过程如下所述。

在本实施例中，可以根据矢量恒等式：

可将式(28)中的第一项变为：

而根据散度定理(或高斯定理)：

则可将式(31)的右边第一项变为：

式中，u代表外边界面。

当采用第一类边界条件时，边界上的场值为零，同样根据矢量恒等式和散度定理，式(28)最后一项也可简化，可以得到：

由此，就可以确定出模拟区域的电场积分方程，即式(34)。

通过这样的方式，即使用伽里金加权余量法，能够快速准确地确定出模拟区域的电场积分方程，提高效率。需要说明的是，在其他一些实现方式中，也可以采用变分方程推导出模拟区域的电场积分方程，因此，此处不应视为对本申请的限定。确定出模拟区域的电场积分方程后，可以对模拟区域中的一个单元进行分析，可以据此确定出模拟区域的总体系数矩阵，具体过程如下所述。

在本实施例中，为了方便描述，将以规则六面体的单元为例，对本实施例的正演方法进行介绍。

如图5所示出的一个六面体单元，x、y和z三个方向的边长分别记为a、b和c，其原点坐标为(x0，y0，z0)。可以将切向量场赋在棱边上，则三个方向上的场分向量可表示为：

其中，n为whitney插值基函数；i，j，k分别表示x分量，y分量，z分量的序号，e表示第e个单元。

具体地：

而单元上的e可表示为whitney插值基函数与棱边上场值乘积的和：

其中，i＝1，2，3，4；带x、y、z脚标的分别表示在x轴、y轴、z轴上的分量。

将整个模拟区域积分分解为各单元积分：

以下，将对式(51)中的各项分别处理，此处有：

对式(51)的第一项进行处理，在单元内有：

对式(51)的第二项进行处理，在单元内有：

对式(51)的第三项进行处理，在单元内有：

对式(51)的第四项进行处理，在单元内有：

其中，

式中，i表示3×3的单元矩阵。

在一个单元内，将式(53)减去式(54)，再减去式(55)，最后再减去式(57)，便可确定出一个单元的系数矩阵。

在本实施例中，确定出一个单元内的系数矩阵后，可以将一个单元内的系数矩阵扩展到所有单元，对整个模拟区域中所有单元的系数矩阵进行整合，便得到模拟区域的总体系数矩阵：

而由于的任意性，可以确定出：

(k1-k2)es＝(k3+k4)ep············(59)。

将所有棱边上的电场一次场值ep及边界上的二次场值es代入式(54)，解之便可以得到每条棱边上的二次场电场值。

具体地，确定出模拟区域的总体系数矩阵后，可以采用无损解精度的直接法加入预设的第一类边界条件，可得到方程组：

如果已知φ3＝p3，φ5＝p5，φ6＝p6，那么求解式(60)等效于求解下列方程组：

通过直接法添加边界条件，消除了方程中第一类边界条件处的量，降低了方程的维数，尤其是在已知量比较多的情况下，能有效地降低矩阵方程的维度，提升处理速度，从而提高三维频率域可控源电磁的正演方法的运行效率。

加入边界条件后，总体系数矩阵则变为：

ke＝p·············(62)，

其中，k是最后形成的总体系数矩阵，即式(59)中的k1-k2，p则为式(59)中的(k3+k4)ep。解方程组，便可得到所有单元的每条棱边上的二次场值。此处，可以使用对称超松弛迭代法预处理，并采用稳定双共轭梯度法求解器进行解方程，求得模拟区域内所有单元上的二次场值。

通过采用从一个单元进行分析后扩散到整个模拟区域的方式确定出模拟区域的总体系数矩阵，而矢量有限元法的应用，可以允许二次场法向分量在电性不连续界面处突变，且考虑了电导率和磁导率的各向异性，更符合实际情况，从而使应用此正演方法的正演系统具有更加广泛的应用范围和更高的准确性。

需要说明的是，在本实施例中，根据一次场电场、和确定二次场电场的过程，可以为本实施例中介绍的方式。在实际的计算过程中，由一次场电场、和确定二次场电场、再确定出总体系数矩阵的过程，不需要在中间过程中计算出二次场的值、中间方程的值等，只需在最后过程中将各个量的实际值代入最后的总体系数矩阵计算出所有单元的二次场值即可。这样在最后计算，从理论上说，可以避免由中间过程的计算带来的误差，从而提高结果的准确性。

步骤s60：根据模拟区域内的一次场电场和二次场电场，确定出所述模拟区域的电场和磁场。

在本实施例中，可以结合模拟区域内的一次场电场值和二次场电场值，求得每个单元的总电场值。而确定出每个单元的总电场值后，可以确定出模拟区域内每个单元的磁场h。具体地：

确定出模拟区域内每个单元的电场和磁场后，就可得到模拟区域的电场和磁场了。

步骤s70：根据所述模拟区域的电场和磁场，确定出所述模拟区域内地表出的使电阻率及相位。

在本实施例中，可以基于模拟区域内的电场和磁场，确定出地表处的张量阻抗z：

其中，zxx、zxy、zyx、zyy分别为：xx、xy、yx、yy模式下的张量阻抗。

根据每个单元的张量阻抗，可以确定出模拟区域内地表处的视电阻率ρij及相位

(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，用于获取每个单元的电导率张量和磁导率张量一次场确定单元14，用于确定出区域内的一次场电场，其中，所述一次场电场表示排除所述异常体时所述模拟区域的电场；二次场确定单元15，用于根据模拟区域内的一次场电场、知确定出模拟区域内的二次场电场，其中，所述二次场电场为所述异常体产生的电场；正演处理单元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，确定出每个单元的以及，根据每个单元的χ1、χ2、χ3、βs、βd和βl，确定出每个单元的

在本实施例中，所述二次场确定单元15，还用于根据模拟区域内的一次场电场、和结合麦克斯韦方程组，确定所述模拟区域的电场积分方程；根据所述电场积分方程，确定出所述模拟区域的总体系数矩阵；根据预设的边界条件和所述总体系数矩阵，确定出所述模拟区域内的二次场电场。

在本实施例中，所述二次场确定单元15，还用于根据模拟区域内的一次场电场，确定出每个单元的电场双旋度方程；由根据每个单元的电场双旋度方程及和确定出的方程点乘电场强度的变分后对整个模拟区域积分，得到中间方程；利用矢量恒等式和散度定理，基于所述中间方程确定出所述模拟区域的电场积分方程。

在本实施例中，在每个单元为六面体时，电场赋在棱边上，所述二次场确定单元15，还用于根据所述积分方程，确定出每个单元的单元方程；根据每个单元的单元方程中包含的系数矩阵，确定出所述模拟区域的总体系数矩阵。

在本实施例中，所述正演处理单元16，还用于基于模拟区域内的一次场电场和二次场电场，确定出模拟区域内的电场；每条棱边上的磁场可根据该棱边及其邻近棱边上的电场确定出，其中，所述模拟区域内所有单元的电场和磁场的总体表示所述模拟区域的电场和磁场。

如图12所示，本申请的实施例还提供一种电子设备20，本申请实施例提供的正演方法，可以由电子设备20执行；而应用本申请提供的正演方法的正演系统，可以基于电子设备20实现其功能。电子设备20可以为终端，例如智能手机、平板电脑、个人电脑、个人数字助理等；电子设备20也可以为服务器，例如网络服务器、云服务器、服务器集群、数据服务器等，此处不作限定。

请参阅图12，电子设备20可以包括：通过网络与外界连接的通信模块22、用于执行程序指令的一个或多个处理器24、总线23和不同形式的存储器21，例如，磁盘、rom、或ram，或其任意组合。

示例性的，存储器21中存储有程序。处理器24可以从存储器21调用并运行这些程序，从而便可以通过运行程序而执行大地电磁的正演方法。处理器24通过对三维频率域可控源电磁的正演方法的执行，处理器24可以对多种不同的地质类型进行正演，以完成对三维频率域可控源电磁的正演方法的执行。

本申请的实施例还提供一种存储介质，存储介质中包含存储的程序，其中，在该程序运行时控制存储介质所在设备执行本实施例的三维频率域可控源电磁的正演方法。

综上所述，本申请的实施例提供一种三维频率域可控源电磁的正演方法、系统、存储介质及电子设备，通过在三维频率域可控源电磁的正演方法中考虑电导率和磁导率的各向异性，可以定性地分析电导率与磁导率各向异性对三维频率域可控源电磁响应的影响，这有利于更加准确地了解和分析实际中复杂的地质结构的电磁响应。并且，本方法也可以为三维频率域可控源电磁电导率与磁导率各向异性反演提供基础。以及，由于本方法同时考虑了电导率和磁导率的各向异性，适用于更多类型的地质结构，因此，可以在更复杂的环境中应用。

在本申请所提供的实施例中，应该理解到，所揭露装置和方法，可以通过其它的方式实现。以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式。

另外，作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。

以上所述仅为本申请的实施例而已，并不用于限制本申请的保护范围，对于本领域的技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。