本发明公开了一种考虑介质各向异性的三维天然源电磁场计算方法，采用基于二次场方法的偏微分方程的求解方法，引入洛伦兹规范，使最终得到的控制方程里仅有矢量位为未知量；在求解偏微分方程时，引入水平方向二维傅里叶变换，将三维偏微分方程转化为一维常微分方程，并利用伽辽金法最终合成三个五对角矩阵。本发明中最后合成的对角方程能用追赶法快速求解，计算效率高。使得一个复杂的各向异性介质三维天然源电磁场数值模拟转化为很多个常微分对角方程组的求解问题，简化了计算，又考虑了介质的各向异性。并提供了相应的计算设备，为大规模正演计算和反演成像提供了新的技术支持。
1.一种考虑介质各向异性的三维天然源电磁场计算方法，其特征在于，包括以下步骤：步骤一，建立直角坐标系中待测位置处的长方体状目标区域，再将目标区域基于立体坐标轴的三个方向进行网格剖分，得到三个方向上的网格节点数及网格坐标；并给定天然场源电磁场的计算频率和极化方向；
步骤二，给网格节点上的导电率赋值，其中导电率为各向异性导电率张量σ；
步骤四，根据天然场源的极化方向、计算频率和背景导电率这三个参数计算网格节点上的背景电场和背景磁场；
步骤五，将关于电磁场的Maxwell方程组基于洛伦兹规范，转换为关于矢量位的控制方程组，然后基于二次场方法得到二次场矢量位的控制方程，再进行二维水平方向傅里叶变换，得到二维波数一维空间域矢量位常微分方程组即控制方程，然后求解得到二维波数一维空间域二次场矢量位的解；
步骤六，利用二维波数一维空间域二次场矢量位与二次电磁场之间的关系求得二维波数一维空间域二次电磁场，再进行二维反傅里叶变换求得二次电磁场；
步骤七，以步骤六得到的二次电磁场，重新代入至步骤五中的二次场矢量位的控制方程组，然后继续迭代求解二次电磁场直到满足预设的误差要求；再将背景电磁场与二次电磁场相加得到总电磁场，从而得到天然场源在设定极化方向条件下的电磁场的解；
将天然源电磁场满足的Maxwell方程组从时间域转化到频率域的表达式为：其中 为拉普拉斯算子，E为电场强度，H 为磁场强度，阻抗率 ，i表示虚数单位， 为角频率，数值上 ，f为计算频率，μ为磁导率，各向异性导纳率，σ为各向异性电导率，ε为介电常数；
引入洛伦兹规范矢量位 ，在洛伦兹规范下， 与 之间的关系式为：式中背景各向异性导纳率 ， 表示背景导纳率， 为背景电导率，为背景介电常数；
基于二次场方法，将场分为背景场和二次场，其中二次场矢量位满足的偏微分方程即控制方程为：
其中波数 ， 为二次场的散射电流密度， 为异常导纳率；
将矢量位 展开为三个方向矢量位： ， 展开为三个方向散射电流：，并将展开后的偏微分方程组进行水平方向二维傅里叶变换，得到二维波数一维空间域的控制方程组为：
式中， 表示偏导， 分别为 对应的二维波数一维空间域二次矢量位， 分别为 对应的二维波数一维空间域散射电流密度，， 为波数域的电磁波传播系数；
采用伽辽金法将上式转化为有限元方程：其中Ne为垂直方向上剖分出的网格个数，Ni（Nj, Np, Nm）为第 个单元中第个节点的二次插值形函数，下标j, p, m分别表示二次插值单元的单元初始节点（j）、中间节点（p ）和末尾节点（m）的符号， 分别为第e个单元的z方向首尾坐标；上式中包含三种类型的积分有：第一类积分形式
将每一项单元积分进行总体合成，得到带宽为5的三个对角方程组：式中， 均为5对角矩阵， 分别为右端项，下标1，2，3分别代表待求的二维波数一维空间域矢量位 相关的系数；求解这三个对角方程即得到二维波数一维空间域矢量位。
2.根据权利要求1所述的方法，其特征在于，所述的步骤一中，天然场源的极化方向包括X极化方向和Y极化方向两种。
3.根据权利要求1所述的方法，其特征在于，所述的步骤二中，各向异性导电率张量σ为
3×3的矩阵，基于对角阵 =diag(σ1，σ2，σ3)和三个欧拉角度αL，αD，αS计算得到；即：其中
αL、αD或αS，Rz和Rx为角度变换矩阵，矩阵中的各个σxx、σxy、σxz、σyx、σyy、σyz、σzx、σzy、σzz为各向异性导电率张量的九个分量；
当介质导电率为各向同性时， 中的σ1=σ2=σ3，为相同的一个常数，αL=αD=αS=0，三个欧拉角均为0；当介质导电率为单轴各向异性时，σ中的σ1=σ2≠σ3，αL=αD=αS=0；当介质导电率为双轴各向异性时， 中的σ1≠σ2≠σ3，αL=αD=αS=0；当介质导电率为任意各向异性时， 中的σ1≠σ2≠σ3，αL、αD、αS不同时为0。
4.根据权利要求3所述的方法，其特征在于，所述的步骤三中，拆分出来的背景导电率记为σb，且同一z坐标上的所有网格节点的σb在数值上是相等的，σb为各向同性或各向异性，σb同样根据对角阵diag(σb1，σb2，σb3)和三个欧拉角度αbL，αbD，αbS计算得到；当背景导电率σb为各向同性时，σb中的σb1=σb2=σb3，αbL=αbD=αbS=0；当不满足各向同性条件时即为各向异性介质；拆分出来的异常导电率记为σa，且σa=σ‑σb。
5.根据权利要求1所述的方法，其特征在于，所述的步骤六中，二次场矢量位 与二次电磁场 之间的关系为：
6.根据权利要求1所述的方法，其特征在于，所述的步骤七包括：先将二次场与背景场相加得到总场：然后采用迭代求解二次电磁场：首先将二次电场设为0，然后利用上式计算总电场，将总电场带入二次场矢量位满足的偏微分方程中，然后再执行步骤五和六得到新的二次电场，将二次电场与背景电场相加得到总电场，并基于下式计算迭代误差 ：其中， 表示编号为 的剖分节点坐标，； 为x，y，z三个方向剖分
的节点数；式中 为本次计算得到的总电场， 为上一次的总电场，当迭代误差时，再将 带入入二次场矢量位满足的偏微分方程中，并重新执行步骤五和六再次得到新的二次电场，再计算总电场，从而循环直至当迭代误差 时，迭代结束，并同时计算二次电场和磁场，再加上背景电磁场，得到最终的总电磁场。
考虑介质各向异性的三维天然源电磁场计算方法
技术领域
[0001] 本发明涉及一种考虑介质各向异性的三维天然源电磁场计算方法。
背景技术
[0002] 自然界地层中，介质的各向异性普遍存在，特别是当岩层经过压实变质作用后会表现出很强的各向异性，因此采用各向同性介质模型进行野外资料的数据处理会导致错误的地质解释。此外，天然源电磁法由于其成本低、计算频率范围广，对地球深部结构了解、地质解释和区域电性层划分等方向上有广泛应用。但现有的天然场源电磁数据均是基于电导率各向同性进行处理，导致最终获得的天然场源电磁场并不准确。
发明内容
[0003] 为了解决目前天然场源电磁数据表现出各向异性，但在计算时却因基于各向同性进行计算而导致出现电磁场结果误差的技术问题，本发明提供一种考虑介质各向异性的三维天然源电磁场计算方法及设备。
[0004] 为了实现上述技术目的，本发明的技术方案是，
[0005] 一种考虑介质各向异性的三维天然源电磁场计算方法，包括以下步骤：
[0006] 步骤一，建立直角坐标系中待测位置处的长方体状目标区域，再将目标区域基于立体坐标轴的三个方向进行网格剖分，得到三个方向上的网格节点数及网格坐标；并给定天然场源电磁场的计算频率和极化方向；
[0007] 步骤二，给网格节点上的导电率赋值，其中导电率为各向异性导电率张量 ；
[0008] 步骤三，将导电率拆分为背景导电率和异常导电率；
[0009] 步骤四，根据天然场源的极化方向、计算频率和背景导电率这三个参数计算网格节点上的背景电场和背景磁场；
[0010] 步骤五，将关于电磁场的Maxwell方程组基于洛伦兹规范，转换为关于矢量位的控制方程组，然后基于二次场方法得到二次场矢量位的控制方程，再进行二维水平方向傅里叶变换，得到二维波数一维空间域矢量位常微分方程组即控制方程，然后求解得到二维波数一维空间域二次场矢量位的解；
[0011] 步骤六，利用二维波数一维空间域二次场矢量位与二次电磁场之间的关系求得二维波数一维空间域二次电磁场，再进行二维反傅里叶变换求得二次电磁场；
[0012] 步骤七，以步骤六得到的二次电磁场，重新代入至步骤五中的二次场矢量位的控制方程组，然后继续迭代求解二次电磁场直到满足预设的误差要求；再将背景电磁场与二次电磁场相加得到总电磁场，从而得到天然场源在设定极化方向条件下的电磁场的解。
[0013] 所述的方法，所述的步骤一中，天然场源的极化方向包括X极化方向和Y极化方向两种。
[0014] 所述的方法，所述的步骤二中，各向异性导电率张量 为3×3的矩阵，基于对角阵和三个欧拉角度 ， ， 计算得到；即：
[0015]
[0016] 其中
[0017] ， ，式中 为三个欧拉角度 ， 或 ， 和 为角度变换矩阵，矩阵中的各个 、 、 、 、
、 、 、 、 为各向异性导电率张量的九个分量；
[0018] 当介质导电率为各向同性时， 中的 ，为相同的一个常数， = ==0，三个欧拉角均为0；当介质导电率为单轴各向异性时， 中的 ， = =
=0；当介质导电率为双轴各向异性时， 中的 ， = = =0；当介质导电
率为任意各向异性时， 中的 ， ， ， 不同时为0。
[0019] 所述的方法，所述的步骤三中，拆分出来的背景导电率记为 ，且同一z坐标上的所有网格节点的 在数值上是相等的， 为各向同性或各向异性， 同样根据对角阵diag(σb1，σb2，σb3)和三个欧拉角度 ， ， 计算得到；当背景导电率 为各向同性时， 中的 ， ；当不满足各向同性条件时即为各向异性介质；拆分出来的异常导电率记为 ，且 。
[0020] 所述的方法，所述的步骤五包括：
[0021] 将天然源电磁场满足的Maxwell方程组从时间域转化到频率域的表达式为：
[0022]
[0023] 其中 为拉普拉斯算子， 为电场强度， 为磁场强度，阻抗率 ，i表示虚数单位， 为角频率，数值上 ， 为计算频率， 为磁导率，各向异性导纳率 ， 为各向异性电导率，为介电常数；
[0024] 引入洛伦兹规范矢量位 ，在洛伦兹规范下， 与 之间的关系式为：
[0025]
[0026] 式中 ， 表示背景导纳率， 为背景电导率， 为背景介电常数；
[0027] 基于二次场方法，将场分为背景场和二次场，其中二次场矢量位满足的偏微分方程即控制方程为：
[0028]
[0029] 其中波数 ， 为二次场的散射电流密度， 为异常导纳率 。
[0030] 将矢量位 展开为三个方向矢量位： ， 展开为三个方向散射电流： ，并将展开后的偏微分方程组进行水平方向二维傅里叶变换，得到二维波数一维空间域的控制方程组为：
[0031]
[0032] 式中， 表示偏导， 分别为 对应的二维波数一维空间域二次矢量位， 分别为 对应的二维波数一维空间域散射电流密度，
， 为波数域的电磁波传播系数；
[0033] 采用伽辽金法将上式转化为有限元方程：
[0034]
[0035] 其中Ne为垂直方向上剖分出的网格个数， 为第 个单元中第个节点的二次插值形函数，下标j, p, m分别表示二次插值单元的单元初
始节点（j ）、中间节点（p）和末尾节点（m）的符号， 分别为第e个单元的z方向首尾坐标；上式中包含三种类型的积分有：
[0036] ① 第一类积分形式
[0037]
[0038] ② 第二类积分形式
[0039]
[0040] ③ 第三类积分形式
[0041]
[0042] 其中，
[0043]
[0044] 将每一项单元积分进行总体合成，得到带宽为5的三个对角方程组：
[0045]
[0046] 式中， 均为5对角矩阵， 分别为右端项，下标1，2，3分别代表待求的二维波数一维空间域矢量位 相关的系数；求解这三个对角方程即
得到二维波数一维空间域矢量位。
[0047] 所述的方法，所述的步骤六中，二次场矢量位 与二次电磁场 之间的关系为：
[0048]
[0049] 从而根据上式求出一维空间域二次电磁场。
[0050] 所述的方法，所述的步骤七包括：
[0051] 先将二次场与背景场相加得到总场：
[0052]
[0053] 然后采用迭代求解二次电磁场：首先将二次电场设为0，然后利用上式计算总电场，将总电场带入二次场矢量位满足的偏微分方程中，然后再执行步骤五和六得到新的二次电场，将二次电场与背景电场相加得到总电场，并基于下式计算迭代误差ε：
[0054]
[0055] 其中， 表示编号为 的剖分节点坐标，i=1,2,⋯Nx，j=1,2,⋯,Ny，k=1,2,⋯,Nz； 为x，y，z三个方向剖分的节点数；式中 为本次计算得到的总电
场， 为上一次的总电场，当迭代误差 时，再将 带入入二次场矢量位满足的偏微分方程中，并重新执行步骤五和六再次得到新的二次电场，再计算总电场，从而循环直至当迭代误差 时，迭代结束，并同时计算二次电场和磁场，再加上背景电磁场，得到最终的总电磁场。
[0056] 本发明的技术效果在于，在对天然源电磁场进行计算时，考虑了介质导电率的各向异性，更符合实际情况；求解天然源电磁场时，采用了基于二次场方法的偏微分方程的求解方法，引入洛伦兹规范，使方程更加简化，最终得到的控制方程里仅有矢量位为未知量；在求解偏微分方程时，引入水平方向二维傅里叶变换，将三维偏微分方程转化为一维常微分方程，计算量和存储需求大大减少，利用伽辽金法最终合成三个五对角矩阵，传统有限元得到的是一个稀疏大型线性方程组，虽然也有并行求解器等快速求解方法，但本发明中最后合成的对角方程能用追赶法快速求解，计算效率高。使得一个复杂的各向异性介质三维天然源电磁场数值模拟转化为很多个常微分对角方程组的求解问题，简化了计算，又考虑了介质的各向异性，采用这种各向异性介质天然源电磁场计算系统能快速准确地计算任意复杂模型，分析不同模型的异常响应特征。为大规模天然源三维各向异性电磁场数值模拟及其反演成像提供了新的技术支持。
[0057] 下面结合附图对本发明作出进一步说明。
附图说明
[0058] 图1为本发明的目标区域示意图，其中(a)为XOY平面示意图，(b)为XOZ剖面示意图；
[0059] 图2为地面y=0m测线上天然源产生的电场本发明数值解和参考值的对比，其中(a)为电场幅值Ex参考解和本发明数值解的对比图，(b)为电场幅值Ey参考解和本发明数值解的对比图，(c)为电场幅值Ez参考解和本发明数值解的对比图；
[0060] 图3为地面y=0m测线上天然源产生的磁场本发明数值解和参考值的对比，其中(a)为磁场幅值Hx参考解和本发明数值解的对比图，(b)为磁场幅值Hy参考解和本发明数值解的对比图，(c)为磁场幅值Hz参考解和本发明数值解的对比图。
具体实施方式
[0061] 本实施例所提供的方法包括以下步骤：
[0062] 1）建立直角坐标系中的长方体目标区域，将目标区域x，y和z方向进行网格剖分，得到三个方向网格节点数为Nx，Ny，Nz。其中一个网格上有两个网格节点，而由于相邻的两个网格会共有一个网格节点，所以这里的网格节点总数量是比网格总数量多1个。然后给定天然场源电磁场的计算频率和极化方向。其中计算频率根据具体的计算需求自行设定即可，极化方向是根据天然源电磁场的情况设定，其中极化方向有X极化方向和Y极化方向两种。
[0063] 2）给网格节点上的导电率赋值，节点导电率为各向异性导电率张量σ，是一个3×3的矩阵，用对角阵 和三个欧拉角度 ， ， 计算得到，引用自文献“Pek Josef, and F.A.M. Santos, 2002, Magnetotelluric impedances and parametric sensitivities for 1‑d anisotropic layered media: Computers &Geoences, 28, 939‑950”，其具体的关系可以写为
[0064] （1）
[0065] 式中
[0066] ， ，式中 为某一角度，在本发明中， 为三个欧拉角度 ， 或 。 和 为角度变换矩阵，矩阵中的各个
、 、 、 、 、 、 、 、 为各向异性导电率张量的九个分量。
[0067] 当介质导电率为各向同性时， 中的 ，为相同的一个常数， = ==0，三个欧拉角均为0；当介质导电率为单轴各向异性时， 中的 ， = =
=0；当介质导电率为双轴各向异性时， 中的 ， = = =0；当介质导电
率为任意各向异性时， 中的 ， ， ， 不同时为0。
[0068] 3）从网格节点赋值的导电率中拆分出一个背景导电率和异常导电率。其中拆分出来的背景导电率记为 ，特征是同一z坐标上的所有网格节点的 在数值上是相等的，可以是各向同性，也可以是各向异性， 同样根据对角阵 和三个欧拉角度 ， ， 计算得到。当背景导电率 为各向同性时， 中的
， 。当不满足各向同性条件时即为各向异性介质。拆分
出来的异常导电率记为 ，数值上 。
[0069] 4）根据天然场源的极化方向、计算频率和背景导电率这三个参数，计算目标区域所有网格节点上的背景电场 和磁场 。当背景导电率为各向同性介质时，背景场计算方法详见文献“陈乐寿, 王光锷. 大地电磁测深法[M]. 地震出版社, 1987.”。当背景导电率为各向异性介质时，背景场计算方法详见文献“陈燊年, 洪清泉, 王建成. 介质为各向异性的电磁场[M]. 科学出版社, 2012.”。
[0070] 5）引入洛伦兹规范，将关于电磁场的Maxwell方程组转换为关于矢量位的控制方程组，基于二次场方法，可得二次场矢量位的控制方程组，进行二维水平方向傅里叶变换，可得二维波数一维空间域矢量位常微分方程组，求解方程可以求得二维波数一维空间域二次场矢量位的解。
[0071] 具体来说，设时谐因子为 ， 为角频率， 为频率，为虚数单位，t为时间。将天然源电磁场满足的Maxwell方程组从时间域转化到频率域的表达式为：
[0072] （2）
[0073] （3）
[0074] （4）
[0075] （5）
[0076] 式中为拉普拉斯算子， 为电场强度， 为磁场强度，阻抗率 ，i表示虚数单位， 为角频率，数值上 ， 为磁导率，各向异性导纳率 ，为各向异性电导率，为介电常数。
[0077] 引入洛伦兹规范矢量位 ，洛伦兹规范下， 与 之间的关系式为：
[0078] （6）
[0079] 式中背景各向异性导纳率 ， 为背景电导率， 为背景介电常数。
[0080] 基于二次场方法，将场分为背景场和二次场，背景场由步骤4）求得，二次场从Maxwell方程组出发，通过推导和简化，可得二次场矢量位满足的偏微分方程为：
[0081] （7）
[0082] 式中波数 ， 为二次场的散射电流密度，式中 为异常导纳率， 。
[0083] 将矢量位 展开为三个方向矢量位（ ）， 展开为三个方向散射电流（ ），并将展开后的偏微分方程组进行水平方向二维傅里叶变换，可得二维波数一维空间域的控制方程组为：
[0084] （8）
[0085] 式中， 分别为 对应的二维波数一维空间域二次矢量位，分别为 对应的二维波数一维空间域散射电流密度，
， 为波数域的电磁波传播系数。
[0086] 采用伽辽金法将式（8）转化为有限元方程：
[0087] （9）
[0088] 式（9）中Ne为垂直方向剖分单元个数即网格个数，N（i Nj, Np, Nm）为第 个单元中第 个节点的二次插值形函数，下标j, p, m分别表示二次插值单元的单元初始节点（j）、中间节点（p）和末尾节点（m）的符号，下面出现的符号说明相同，不再赘述。
分别为第e个单元的z方向首尾坐标。关于二次插值函数表达式等详细的说明见文献“徐世浙. 地球物理中的有限单元法[M]. 科学出版社, 1994.”。式（9）中包含三种类型的积分有：
[0089] 第一类积分形式：
[0090]
[0091] 第二类积分形式：
[0092]
[0093] 第三类积分形式：
[0094]
[0095] 其中，
[0096]
[0097] 将每一项单元积分进行总体合成，可以得到带宽为5的三个对角方程组，[0098] （13）
[0099] 式中， 均为5对角阵， 分别为右端项，下标1，2，3分别代表待求的二维波数一维空间域矢量位 相关的系数。求解这三个方程即可求
得二维波数一维空间域矢量位。
[0100] 6）利用二维波数一维空间域二次场矢量位与二次电场之间的关系求得二维波数一维空间域二次电磁场，进行二维反傅里叶变换求得二次电磁场；
[0101] 其中二次场矢量位与二次电磁场之间的关系为
[0102] （14）
[0103] 利用式（11）在二维波数域一维空间域中的关系式，可得二维波数域一维空间域电磁场，进行二维反傅里叶变换，可得空间域二次场电磁场。
[0104] 所述步骤7）中，二次场与背景场相加得到总场：
[0105] （15）
[0106] 二次电磁场的迭代求解的具体过程为：首先将二次电场设为0，然后利用式（12）计算总电场，将总电场带入式（7）中，然后利用步骤5）和6）得到新的二次电场，将二次电场与背景电场相加得到总电场，迭代误差 的计算公式为：
[0107] 其中， 表示编号为 的剖分节点坐标，； 为x，y，z三个方向剖分
的节点数。式中 为本次计算得到的总电场， 为上一次的总电场，当迭代误差时，将 带入式（7）中，利用步骤5）和6）再次得到新的二次电场，再计算总电场，循环往复，直至当迭代误差 时，迭代结束，并同时计算二次电场和磁场，再加上背景电磁场，得到最终的总电磁场。
[0108] 7）用迭代解法求解二次电磁场，然后将背景电磁场与二次电磁场相加得到总电磁场，得到天然场源在设定极化方向条件下的电磁场的解。
[0109] 下面通过模型实例，验证本发明提供的考虑介质各向异性的一种三维天然源电磁场计算方法的正确性和有效性。测试的计算机为Intel(R) Core(TM) i7‑6700HQ CPU 主频为2.60GHz ,内存为16GB、64位win10系统。
[0110] 模型XOY平面投影如图1所示，背景为均匀半空间介质，上半空间为空气，空气导电‑12率 0=10 S/m，下半空间的背景导电率为 =0.01 S/m，频率为10Hz，计算X方向极化的天然源在地面的电磁场响应，用文献“Liu Chang‑sheng, Ren Zheng‑yong, Tang Jing‑tian, et al.Three‑dimensional magnetotellurics modeling using edgebased finite‑element unstructured meshes[J]. Applied Geophysics, 2008(03):1855‑‑ikz
1859.”中的方法的计算结果为参照，验证方法的正确性。地下背景电场的计算公式为e ，地面电场幅度为1 V/m。计算范围x方向‑1000 1000m，y方向‑1000 1000m，z方向计算范围设~ ~
为0 1000m，剖分网格节点个数51×51×51，三个方向均匀剖分，△x、△y均为40m，△z为~
20m，异常体范围x方向‑100 100m，y方向‑200 200m，z方向200 400m，异常体导电率为任意~ ~ ~
各向异性，参数为： =diag(0.1，0.2，0.5)；三个欧拉角度 =5°、 =10°、 =15°。
[0111] 图2 和图3分别为X方向极化的天然场源在地面y=0m测线上电场和磁场本发明数值解和参考值的对比。由两图可以看出，Ez分量存在误差，这是由于Ez在地面为0，而图中两种方法解的Ez数值都很小，因此两种方法的误差可以忽略不计，认为都是正确的；除Ez外其他五个电场分量和磁场分量的吻合程度很好，计算误差小，验证了本实施例方法的正确性。本实施例的计算节点数为51×51×51=132651，本发明算法迭代终止时一共迭代了13次，每次迭代耗时约1.2 s, 计算总时间为16.2s，占用内存为120.5 MB，传统有限元算法占用内存约3.5 GB，迭代一次耗时约20.1s，正演计算总耗时约220.5 s。从迭代一次耗时和占用内存对比可以看出，本发明采用傅里叶变换，求解方程组的未知量由三维转换为一维，占用内存大幅度减小，采用追赶法求解方程组，迭代一次计算速度快，耗时短，总耗时速度提升了
10倍，表明本发明与传统有限元算法相比效率优势明显，证明本发明考虑介质各向异性的这种三维天然源电磁场计算方法准确且效率高。
[0112] 以上所述仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是在本发明的发明构思下，利用本发明说明书及附图内容所作的等效结构变换，或直接/间接运用在其他相关的技术领域均包括在本发明的专利保护范围内。
