一种基于积分方程法的三维大地电磁数值模拟方法与流程

1.本发明涉及数值模拟技术领域，具体是一种基于积分方程法的三维大地电磁数值模拟方法。


背景技术：

2.由于大地电磁测深法（mt）野外工作成本低、操作简单，具有勘探深度大、不受高阻屏蔽和对高导层敏感等优点，已成为矿产资源探查、工程环境探查和深部地质构造研究的一种有效的勘探方法。
3.大地电磁三维数值模拟方法主要有有限差分法、有限单元法、有限体积法和积分方程法。其中，积分方程法由于仅需对异常体进行剖分，占用内存小，计算量小等优点对于开发高效高精度的正演算法具有一定的研究价值。而现有的数值模拟方法均基于空间域的微分或积分方程进行计算，最终合成空间域线性方程组。大规模复杂条件下最终的方程组维数巨大，导致现有方法计算量大、所需存储量大，影响正演效率。而三维大地电磁正演计算速度是影响反演成像效率与人机交互定量解释的关键因素，因此，高效、高精度的大规模复杂地质条件下三维数值模拟方法仍然是当前正演研究的主要焦点之一。


技术实现要素：

4.针对上述现有技术中的不足，本发明提供一种基于积分方程法的三维大地电磁数值模拟方法，能适用于各频段大地电磁数值模拟，且精确度高。
5.为实现上述目的，本发明提供一种基于积分方程法的三维大地电磁数值模拟方法，包括如下步骤：步骤1，建立大地电磁积分方程，并在直角坐标系中建立目标区域，使三维导电率异常体包含在目标区域内，将目标区域沿x、y、z方向进行剖分，得到空间域离散的采样节点，并根据地下地质条件导电率的分布，给定目标区域采样节点的导电率，设置计算频率；步骤2，将目标区域的导电率分布拆分为背景导电率和异常导电率，设置大地电磁场源的极化方向，然后根据背景导电率和计算频率计算每个采样节点的背景电场和背景磁场，将背景电场设置为初始电场；步骤3，将目标区域采样节点的异常导电率与初始电场相乘得到采样节点的散射电流，再将散射电流进行水平方向二维傅里叶变换，得到空间波数域散射电流；步骤4，将大地电磁满足的二次电场的积分方程进行傅里叶变换，得到空间波数域二次电场积分方程，将空间波数域散射电流代入空间波数域二次电场积分方程中，计算空间波数域二次电场；步骤5，将空间波数域二次电场进行二维傅里叶反变换，得到空间域二次电场，将二次电场与背景电场相加得到空间域总电场；步骤6，判断得到的空间域总电场是否满足给定的迭代收敛条件，若满足迭代收敛条件，则输出总电场，并利用总电场求导计算总磁场，否则将计算得到的总电场替代初始电
场后返回步骤3；步骤7，分别计算x和y方向极化的大地电磁场源正演得到的电磁场，利用大地电磁视电阻率计算公式求取测点处的视电阻率和相位。
6.作为上述技术方案的进一步改进，步骤1中，所述将目标区域沿x、y、z方向进行剖分，具体为：将目标区域沿x、y、z方向进行剖分，x、y、z方向的采样节点数量分别设为n
x
、n
y
、n
z
。
7.作为上述技术方案的进一步改进，步骤2中，所述目标区域的导电率；其中，为总导电率；为背景导电率，在水平方向不变化，仅在垂直z方向上变化；为异常导电率，当异常体不存在时，，当异常体存在时，等于给定节点的导电率减去背景导电率。
8.作为上述技术方案的进一步改进，步骤4中，所述空间波数域二次电场积分方程为：式中，为空间波数域二次电场，k
x
、k
y
分别表示x、y方向傅里叶变换的波数，(k
x
、k
y
、z
k
)表示空间波数域二次电场密度在空间波数域的坐标，为空间波数域格林函数，为空间波数域散射电流，表示垂向坐标，z1，z2表示垂向积分的长度范围。
9.作为上述技术方案的进一步改进，步骤6中，给定的空间域总电场的迭代收敛条件为：式中，为期望收敛时的数值精度，为坐标为的节点正演得到的新的总电场，为坐标为的节点的初始电场。
10.本发明提供的一种基于积分方程法的三维大地电磁数值模拟方法，利用水平方向二维傅里叶变换，将三维空间域积分方程拆分成很多波数下的一维积分方程，计算量和存储需求大大减小，能提高数值模拟的计算效率；且该数值模拟方法能适用于大地电磁研究的频段104~10
‑
4 hz内的任意频率，采用迭代法求解电磁场，计算速度快，占用内存小，使得该发明方法计算效率高，计算精度高，适用于大规模复杂地质条件的大地电磁数值模拟，能有效地提高三维大地电磁数据反演的精度和效率。
附图说明
11.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图示出的结构获得其他的附图。
12.图1为本发明实施例中基于积分方程法的三维大地电磁数值模拟方法的流程示意图；图2为本发明实施例中示例的目标区域示意图，其中，（a）为目标区域的侧视图，（b）为目标区域的俯视图；图3为本发明实施例中示例的大地电磁视电阻率和相位的计算值对比示意图，其中，(a)为大地电磁视电阻率ρ
xy
的示意图，(b)为大地电磁视电阻率ρ
yx
的示意图，(c)为大地电磁相位φ
xy
的示意图，(d)为大地电磁相位φ
yx
的示意图；图4为本发明实施例中示例的大地电磁视电阻率和相位的计算值和参考值的相对误差图，其中，(a)为两种方法大地电磁视电阻率ρ
xy
的相对误差示意图，(b)为两种方法大地电磁视电阻率ρ
yx
的相对误差示意图，(c)为两种方法大地电磁相位φ
xy
的相对误差示意图，(d)为两种方法大地电磁相位φ
yx
的相对误差示意图。
13.本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。
具体实施方式
14.下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
15.需要说明，本发明实施例中所有方向性指示（诸如上、下、左、右、前、后
……
）仅用于解释在某一特定姿态（如附图所示）下各部件之间的相对位置关系、运动情况等，如果该特定姿态发生改变时，则该方向性指示也相应地随之改变。
16.另外，在本发明中如涉及“第一”、“第二”等的描述仅用于描述目的，而不能理解为指示或暗示其相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括至少一个该特征。在本发明的描述中，“多个”的含义是至少两个，例如两个，三个等，除非另有明确具体的限定。
17.在本发明中，除非另有明确的规定和限定，术语“连接”、“固定”等应做广义理解，例如，“固定”可以是固定连接，也可以是可拆卸连接，或成一体；可以是机械连接，也可以是电连接，还可以是物理连接或无线通信连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通或两个元件的相互作用关系，除非另有明确的限定。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。
18.另外，本发明各个实施例之间的技术方案可以相互结合，但是必须是以本领域普通技术人员能够实现为基础，当技术方案的结合出现相互矛盾或无法实现时应当认为这种技术方案的结合不存在，也不在本发明要求的保护范围之内。
19.如图1所示为本实施例公开的一种基于积分方程法的三维大地电磁数值模拟方
法，包括如下步骤：步骤1，建立大地电磁积分方程，并在直角坐标系中建立目标区域，使三维导电率异常体包含在目标区域内，将目标区域沿x、y、z方向进行剖分，得到空间域离散的采样节点，并根据地下地质条件导电率的分布，给定目标区域采样节点的导电率，设置计算频率。
20.将目标区域沿x、y、z方向进行剖分，具体为：将目标区域沿x、y、z方向进行剖分，将目标区域分为若干空间单元，x、y、z方向的采样节点数量分别设为n
x
、n
y
、n
z
。
21.步骤2，将目标区域的导电率分布拆分为背景导电率和异常导电率，设置大地电磁场源的极化方向，然后根据背景导电率和计算频率计算每个采样节点的背景电场和背景磁场，将背景电场设置为初始电场。其中，目标区域的导电率；其中，为总导电率；为背景导电率，在水平方向不变化，仅在垂直z方向上变化；为异常导电率，当异常体不存在时，，当异常体存在时，等于给定节点的导电率减去背景导电率。计算得到的各采样节点在x，y，z方向的背景电场的三分量、、，其中，表示编号为的采样节点坐标，。至于根据背景导电率和计算频率计算每个采样节点的背景电场和背景磁场的具体过程则为所属领域的常规手段，本实施例中不再赘述。
22.步骤3，将目标区域采样节点的异常导电率与初始电场相乘得到采样节点的散射电流，再将散射电流进行水平方向二维傅里叶变换，得到空间波数域散射电流。其中，将散射电流进行水平方向二维傅里叶变换的具体过程则为所属领域的常规手段，本实施例中不再赘述。
23.步骤4，将大地电磁满足的二次电场的积分方程进行傅里叶变换，得到空间波数域二次电场积分方程，将空间波数域散射电流代入空间波数域二次电场积分方程中，计算空间波数域二次电场。其中，空间波数域二次电场积分方程为：式中，为空间波数域二次电场，k
x
、k
y
分别表示x、y方向傅里叶变换的波数，(k
x
、k
y
、z
k
)表示空间波数域二次电场密度在空间波数域的坐标，为空间波数域格林函数，为空间波数域散射电流，z
s
表示垂向坐标，z1，z2表示垂向积分的长度范围。
24.步骤5，将空间波数域二次电场进行二维傅里叶反变换，得到空间域二次电场，将二次电场与背景电场相加得到空间域总电场。其中，将空间波数域二次电场进行二维傅里叶反变换的具体过程则为所属领域的常规手段，本实施例中不再赘述。
25.步骤6，判断得到的空间域总电场是否满足给定的迭代收敛条件，若满足迭代收敛条件，则输出总电场，并利用总电场与背景磁场求导计算总磁场，否则将计算得到的总电场替代初始电场后返回步骤3。给定的空间域总电场的迭代收敛条件为：
式中，为期望收敛时的数值精度，为坐标为的节点正演得到的新的总电场，为坐标为的节点的初始电场。
26.步骤7，分别计算x和y方向极化的大地电磁场源正演得到的电磁场，利用大地电磁视电阻率计算公式求取测点处的视电阻率和相位。
27.下面结合设计模型实例，对本实施例提供的基于积分方程法的大地电磁三维数值模拟方法进行精度检验，测试的计算机为intel(r) core(tm) i7
‑
6700hq cpu 主频为2.60ghz ,内存为16gb、64位win10系统，算法在microsoft visual studio 2015开发平台上运行。
28.模型xoy平面投影如图2所示，背景为均匀半空间介质，上半空间为空气，空气导电率σ0=10
‑
12 s/m，下半空间的背景导电率为=0.01 s/m，频率为1hz，进行大地电磁场三维数值模拟。计算范围x方向
‑
1000~1000m，y方向
‑
1000~1000m，z方向计算范围设为0~1000m，剖分节点个数101
×
101
×
101，三个方向均匀剖分，x方向上的节点间隔δx、y方向上的节点间隔δy和z方向上的节点间隔δz均为10m，异常体范围x方向
‑
100~100m，y方向
‑
200~200m，z方向200~400m，异常体导电率σ=0.1 s/m。用美国犹他大学开发的基于积分方程法的三维正演软件intem3d的计算结果为参照，验证方法的正确性。
29.本实施例的大地电磁三维积分方程数值模拟计算达到期望收敛相对残差为10
‑4时，需要迭代17次，计算总时间为20.35s，占用内存为832.5 mb，由此可见本实施例方法计算速度快，占用内存小，具有很高的计算效率。图3为z=0平面本实施例方法和intem3d软件计算得到的大地电磁视电阻率和相位图，图4为二者之间的相对误差。从图4中可看出，两种算法视电阻率和相位的相对误差均小于0.4%，验证了本实施例方法的正确性，说明本实施例方法具有很高的计算精度。
30.以上所述仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是在本发明的发明构思下，利用本发明说明书及附图内容所作的等效结构变换，或直接/间接运用在其他相关的技术领域均包括在本发明的专利保护范围内。