一种任意位置高精度接收器响应的计算方法、装置和设备

1.本技术涉及数值模拟领域，特别是涉及一种任意位置高精度接收器响应的计算方法、装置和计算机设备。


背景技术：

2.大地电磁测深法是以天然交变电磁场为源，假设平面波垂直入射大地，通过测量地表电磁场来研究地下电性结构的一种地球物理方法。利用宽频带的天然电磁场作为场源来探测地下介质的电性变化，具有探测深度范围大和深部高导结构非常敏感的优势。
3.大地电磁测深正演模拟就是通过解析或者数值模拟的方法，在给定地下介质分布及激发源的情况下，计算相应的地球物理响应的过程。我们可以基于地球物理正演来研究不同物理模型下的响应的分布规律，从而可以指导实际勘探工作。
4.通过正演求解有限差分形成的线性方程组后，可以得到剖分网格各棱边上的电场分量值。但为了获得任意接收器位置的电场，通常还需要在选定区域内进行插值。由于实际研究的地质环境比较复杂，如：复杂的地形、电阻率的剧烈变化和各向异性介质等。目前，比较常用的是移动最小二乘法进行插值，这种方法精度比较高，使用比较广泛，但是存在一些问题：如目标函数的构建复杂和多个未知参数求解繁琐。因此，现有技术存在效率较低、适应性不佳的问题。


技术实现要素：

5.基于此，有必要针对上述技术问题，提供一种能够提高接收器响应正演效率和适应性的任意位置高精度接收器响应的计算方法、装置、计算机设备和存储介质。
6.一种任意位置高精度接收器响应的计算方法，所述方法包括：
7.构建内部包含勘探目标的三维长方体模型，对所述三维长方体模型进行网格剖分，形成若干个长方体单元，给每个长方体单元的电阻率赋值，得到三维异常体模型；
8.根据所述三维长方体模型网格剖分参数和预设的频率参数得到电场控制方程的系数矩阵；
9.根据所述三维长方体模型网格剖分参数和所述频率参数得到所述三维异常体模型边界面上的电场，根据三维大地电磁场源和边界条件得到电场控制方程的右端向量；
10.根据所述系数矩阵和所述右端向量确定所述三维异常体模型对应的电场控制方程；
11.以所述三维长方体模型为坐标原点进行坐标转换，确定每个长方体单元的坐标，通过人机交互界面获取每个接收器的坐标位置，根据每个长方体单元的坐标和接收器的坐标位置确定接收器所处的长方体单元，并进一步根据所述电场控制方程确定所述长方体单元各个顶点的电场分量，得到电场对角矩阵；
12.根据接收器的坐标位置和接收器所在长方体单元的坐标，确定接收器测量点x、y和z三个方向棱边长度的权重值，并进一步确定接收器所在长方体单元各个顶点的总权重
值；
13.根据所述电场对角矩阵和所述总权重值确定接收器上的电场计算量，若所述电场计算量在周围点电场最小值和最大值之间，则判定接收器的响应是精确的，输出接收器响应，否则，细化所述三维长方体模型的网格剖分并重新计算，直到能够判定接收器的响应是精确的，输出接收器响应。
14.在其中一个实施例中，还包括：根据所述三维长方体模型网格剖分参数确定双旋度算子；
15.根据所述双旋度算子和预设的频率参数得到电场控制方程的系数矩阵为：
[0016][0017]
其中，a表示系数矩阵，表示旋度算子，表示双旋度算子，ω表示角频率，通过ω＝2πf求取，f表示给定的频率，μ表示磁导率，其值为4π
×
10-7
，ρ表示每个长方体单元的电阻率。
[0018]
在其中一个实施例中，还包括：当三维大地电磁场源x极化时，空气层顶部x方向的电场e
x
为场源，取其电场值e
x
＝1，通过二维正演控制方程求取三维异常体模型所在区域的左侧面、右侧面的边界条件，此时前侧面和后侧面边界上的电场取0，通过插值左右两测边界最底下的电场值得到底面边界条件；其中，表示对应长方体单元y方向长度值δy，表示表示对应长方体单元z方向长度值δz；
[0019]
当三维大地电磁场源y极化时，空气层顶部y方向的电场ey为场源，取其电场值ey＝1，通过二维正演控制方程求取三维异常体模型所在区域的前侧面、后侧面的边界条件，此时左侧面和右侧面边界上的电场取0，通过插值前后两测边界最底下的电场值得到底面边界条件；其中，表示对应长方体单元x方向长度值δx；
[0020]
由所述三维大地电磁场源和边界条件构成电场控制方程的右端向量。
[0021]
在其中一个实施例中，还包括：根据所述系数矩阵和所述右端向量确定所述三维异常体模型对应的电场控制方程为：
[0022]
ae＝b
[0023]
其中，b表示右端向量，e表示电场。
[0024]
在其中一个实施例中，还包括：确定所述长方体单元各个顶点的电场分量，其中，沿x方向的电场为沿x方向的电场为和沿y方向的电场沿y方向的电场和沿z方向的电场沿z方向的电场和
[0025]
得到电场对角矩阵f
x
,fy,fz。
[0026]
在其中一个实施例中，还包括：根据接收器的坐标位置和接收器所在长方体单元的坐标，确定接收器测量点x、y和z三个方向棱边长度的权重值为：
[0027][0028]
其中，α
x1
，α
x2
，α
y1
，α
y2
，α
z1
和α
z2
分别为接收器测量点x、y和z三个方向棱边长度的权重值，(l
x
,ly,lz)表示接收器的坐标位置，(xi,yi,zi)表示第i个长方体单元的坐标；
[0029]
进一步确定接收器所在长方体单元各个顶点的总权重值s为：
[0030][0031]
其中，λ为顶点的权重值。
[0032]
在其中一个实施例中，还包括：根据所述电场对角矩阵和所述总权重值确定接收器上的电场计算量为：
[0033][0034]
其中，e
rx
、e
ry
、e
rz
分别为接收器在x、y和z方向上的电场计算量。
[0035]
在其中一个实施例中，还包括：若e
rx
、e
ry
、e
rz
满足下式：
[0036]erxmin
≤e
rx
≤e
rxmax
且e
rymin
≤e
ry
≤e
rymax
且e
rzmin
≤e
rz
≤e
rzmax
[0037]
则判定接收器的响应是精确的，输出接收器响应；其中，e
rxmin
,e
rymin
和e
rzmin
和e
rxmax
,e
rymax
和e
rzmax
分别为接收器所在长方体单元周围各个顶点电场各个分量的最小值和最大值。
[0038]
一种任意位置高精度接收器响应的计算装置，所述装置包括：
[0039]
模型构建模块，用于构建内部包含勘探目标的三维长方体模型，对所述三维长方体模型进行网格剖分，形成若干个长方体单元，给每个长方体单元的电阻率赋值，得到三维异常体模型；
[0040]
电场控制方程确定模块，用于根据所述三维长方体模型网格剖分参数和预设的频率参数得到电场控制方程的系数矩阵；根据所述三维长方体模型网格剖分参数和所述频率
参数得到所述三维异常体模型边界面上的电场，根据三维大地电磁场源和边界条件得到电场控制方程的右端向量；根据所述系数矩阵和所述右端向量确定所述三维异常体模型对应的电场控制方程；
[0041]
电场对角矩阵确定模块，用于以所述三维长方体模型为坐标原点进行坐标转换，确定每个长方体单元的坐标，通过人机交互界面获取每个接收器的坐标位置，根据每个长方体单元的坐标和接收器的坐标位置确定接收器所处的长方体单元，并进一步根据所述电场控制方程确定所述长方体单元各个顶点的电场分量，得到电场对角矩阵；
[0042]
总权重值确定模块，用于根据接收器的坐标位置和接收器所在长方体单元的坐标，确定接收器测量点x、y和z三个方向棱边长度的权重值，并进一步确定接收器所在长方体单元各个顶点的总权重值；
[0043]
响应输出模块，用于根据所述电场对角矩阵和所述总权重值确定接收器上的电场计算量，若所述电场计算量在周围点电场最小值和最大值之间，则判定接收器的响应是精确的，输出接收器响应，否则，细化所述三维长方体模型的网格剖分并重新计算，直到能够判定接收器的响应是精确的，输出接收器响应。
[0044]
一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现以下步骤：
[0045]
构建内部包含勘探目标的三维长方体模型，对所述三维长方体模型进行网格剖分，形成若干个长方体单元，给每个长方体单元的电阻率赋值，得到三维异常体模型；
[0046]
根据所述三维长方体模型网格剖分参数和预设的频率参数得到电场控制方程的系数矩阵；
[0047]
根据所述三维长方体模型网格剖分参数和所述频率参数得到所述三维异常体模型边界面上的电场，根据三维大地电磁场源和边界条件得到电场控制方程的右端向量；
[0048]
根据所述系数矩阵和所述右端向量确定所述三维异常体模型对应的电场控制方程；
[0049]
以所述三维长方体模型为坐标原点进行坐标转换，确定每个长方体单元的坐标，通过人机交互界面获取每个接收器的坐标位置，根据每个长方体单元的坐标和接收器的坐标位置确定接收器所处的长方体单元，并进一步根据所述电场控制方程确定所述长方体单元各个顶点的电场分量，得到电场对角矩阵；
[0050]
根据接收器的坐标位置和接收器所在长方体单元的坐标，确定接收器测量点x、y和z三个方向棱边长度的权重值，并进一步确定接收器所在长方体单元各个顶点的总权重值；
[0051]
根据所述电场对角矩阵和所述总权重值确定接收器上的电场计算量，若所述电场计算量在周围点电场最小值和最大值之间，则判定接收器的响应是精确的，输出接收器响应，否则，细化所述三维长方体模型的网格剖分并重新计算，直到能够判定接收器的响应是精确的，输出接收器响应。
[0052]
一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现以下步骤：
[0053]
构建内部包含勘探目标的三维长方体模型，对所述三维长方体模型进行网格剖分，形成若干个长方体单元，给每个长方体单元的电阻率赋值，得到三维异常体模型；
[0054]
根据所述三维长方体模型网格剖分参数和预设的频率参数得到电场控制方程的系数矩阵；
[0055]
根据所述三维长方体模型网格剖分参数和所述频率参数得到所述三维异常体模型边界面上的电场，根据三维大地电磁场源和边界条件得到电场控制方程的右端向量；
[0056]
根据所述系数矩阵和所述右端向量确定所述三维异常体模型对应的电场控制方程；
[0057]
以所述三维长方体模型为坐标原点进行坐标转换，确定每个长方体单元的坐标，通过人机交互界面获取每个接收器的坐标位置，根据每个长方体单元的坐标和接收器的坐标位置确定接收器所处的长方体单元，并进一步根据所述电场控制方程确定所述长方体单元各个顶点的电场分量，得到电场对角矩阵；
[0058]
根据接收器的坐标位置和接收器所在长方体单元的坐标，确定接收器测量点x、y和z三个方向棱边长度的权重值，并进一步确定接收器所在长方体单元各个顶点的总权重值；
[0059]
根据所述电场对角矩阵和所述总权重值确定接收器上的电场计算量，若所述电场计算量在周围点电场最小值和最大值之间，则判定接收器的响应是精确的，输出接收器响应，否则，细化所述三维长方体模型的网格剖分并重新计算，直到能够判定接收器的响应是精确的，输出接收器响应。
[0060]
上述任意位置高精度接收器响应的计算方法、装置、计算机设备和存储介质，通过构建内部包含勘探目标的三维长方体模型，对三维长方体模型进行网格剖分，给每个长方体单元的电阻率赋值，得到三维异常体模型；分别计算电场控制方程的系数矩阵和右端向量，确定三维异常体模型对应的电场控制方程；通过人机交互界面获取每个接收器的坐标位置，确定接收器所处的长方体单元，并进一步根据电场控制方程确定长方体单元各个顶点的电场分量，得到电场对角矩阵；根据接收器的坐标位置和接收器所在长方体单元的坐标，确定接收器测量点x、y和z三个方向棱边长度的权重值，并进一步确定接收器所在长方体单元各个顶点的总权重值；根据电场对角矩阵和总权重值确定接收器上的电场计算量，若电场计算量在周围点电场最小值和最大值之间，则判定接收器的响应是精确的，输出接收器响应，否则，细化三维长方体模型的网格剖分并重新计算，直到能够判定接收器的响应是精确的，输出接收器响应。本发明可以避免主流方法存在的目标函数复杂和参数求解繁琐等问题，同时又解决了线性插值在电阻率剧烈变化时候精度不高的问题，能够满足任意位置接收器响应的高精度计算需求。
附图说明
[0061]
图1为一个实施例中任意位置高精度接收器响应的计算方法的流程示意图；
[0062]
图2为另一个实施例中任意位置高精度接收器响应的计算方法的流程示意图；
[0063]
图3为一个实施例中高-低阻混合矿体模型示意图；
[0064]
图4为一个实施例中得到的电场的参考解、数值解以及相对误差图，其中(a)为电场的实部的参考解，(b)为电场的实部的数值解，(c)为电场的实部的相对误差；(d)为电场的虚部的参考解，(e)为电场的虚部的数值解，(f)为电场的虚部的相对误差；
[0065]
图5为一个实施例中任意位置高精度接收器响应的计算装置的结构框图；
[0066]
图6为一个实施例中计算机设备的内部结构图。
具体实施方式
[0067]
为了使本技术的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本技术进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本技术，并不用于限定本技术。
[0068]
在一个实施例中，如图1所示，提供了一种任意位置高精度接收器响应的计算方法，包括以下步骤：
[0069]
步骤102，构建内部包含勘探目标的三维长方体模型，对三维长方体模型进行网格剖分，形成若干个长方体单元，给每个长方体单元的电阻率赋值，得到三维异常体模型。
[0070]
本发明中的网格剖分方式可以是均匀剖分，使得剖分成的各长方体单元的尺寸相同，也可以是非均匀剖分，剖分成的各长方体单元可以具有不同的长宽高尺寸。
[0071]
根据勘探目标的电阻率分布，给每个小长方体单元的电阻率赋值，每一个长方体单元的电阻率为常值，不同长方体单元的电阻率值不同，得到刻画任意电阻率分布的三维异常体模型。
[0072]
步骤104，根据三维长方体模型网格剖分参数和预设的频率参数得到电场控制方程的系数矩阵。
[0073]
具体地，根据三维长方体模型网格剖分参数确定双旋度算子，双旋度算子的计算方式如下：
[0074]
三维异常体模型中各个小长方体单元x、y、z方向的棱边长度，构成长度元矩阵；
[0075]
三维异常体模型中各个小长方体单元的表面积，构成面积元矩阵s；
[0076]
三维异常体模型中各个小长方体单元的体积，构成体积元矩阵v；
[0077]
双旋度算子，其中o为拓扑矩阵，由-1和1两种元素组成。拓扑矩阵o作用于由长度元与场值乘积组成的矩阵，实现积分的求和运算，几何映射关系为长度元l到面积元s，每行对应一个面积元。
[0078]
根据双旋度算子和预设的频率参数得到电场控制方程的系数矩阵为：
[0079][0080]
其中，a表示系数矩阵，表示旋度算子，表示双旋度算子，ω表示角频率，通过ω＝2πf求取，f表示给定的频率，μ表示磁导率，其值为4π
×
10-7
，ρ表示每个长方体单元的电阻率。
[0081]
步骤106，根据三维长方体模型网格剖分参数和频率参数得到三维异常体模型边界面上的电场，根据三维大地电磁场源和边界条件得到电场控制方程的右端向量。
[0082]
求得的三维异常体模型边界面上的电场即为边界条件。
[0083]
具体地，右端向量b由三维大地电磁场源和边界条件构成的向量，其中场源指在空气层垂直向下入射的平面波。三维大地电磁场源x极化时，三维大地电磁场源位于空气层顶端，其电场值e
x
为1，三维异常体模型所在区域的左侧面、右侧面和底端为边界条件，其中左
侧面和右侧面边界条件的电场值e
x
通过调用二维正演控制方程求取，表示对应长方体单元y方向长度值δy，表示表示对应长方体单元z方向长度值δz，此时前侧面和后侧面边界上的电场取0，底部边界电场值e
x
通过插值左右两测边界最底下的电场值得到。y极化时，三维大地电磁场源位于空气层顶端，其电场值ey为1，三维异常体模型所在区域的前侧面、后侧面和底端为边界条件，其中前侧面和后侧面边界条件的电场值ey通过调用二维正演控制方程求取，表示对应长方体单元x方向长度值δx，此时左侧面和右侧面边界上的电场取0，底部边界电场值ey通过插值前后两测边界最底下的电场值得到。
[0084]
由三维大地电磁场源和边界条件构成电场控制方程的右端向量。
[0085]
步骤108，根据系数矩阵和右端向量确定三维异常体模型对应的电场控制方程。
[0086]
具体地，根据系数矩阵和右端向量确定三维异常体模型对应的电场控制方程为：
[0087]
ae＝b
[0088]
其中，b表示右端向量，e表示电场。
[0089]
步骤110，以三维长方体模型为坐标原点进行坐标转换，确定每个长方体单元的坐标，通过人机交互界面获取每个接收器的坐标位置，根据每个长方体单元的坐标和接收器的坐标位置确定接收器所处的长方体单元，并进一步根据电场控制方程确定长方体单元各个顶点的电场分量，得到电场对角矩阵。
[0090]
坐标转换是将坐标原点由研究区域西北方向角点转换到研究区域地表的中心处，转换方式如下：
[0091]
结合长方体单元x、y和z方向的棱边长度δx,δy和δz，先通过求和公式得到各个测量点的坐标位置(xi,yj,zk)：
[0092]
(xi,yj,zk)＝(sum(δxi),sum(yj),sum(zk))
[0093]
其中(xi,yj,zk)表示以研究区域西北方向角点为原点的坐标位置，xi表示x方向第i个网格单元距离原点的距离，yj表示y方向第j个网格单元距离原点的距离和zk表示z方向第k个网格单元的测量点坐标位置。通过坐标转换将坐标原点转换到研究区域地表的中心处：
[0094]
(xi,yj,zk)＝(x
i-sum(δx)/2,y
j-sum(y)/2,z
k-sum(z)/2)
[0095]
此时(xi,yj,zk)表示以研究区域西北方向角点为原点的坐标位置。
[0096]
由电场控制方程可以确定长方体单元棱边上的电场分量，进而可以确定长方体单元顶点的电场分量。
[0097]
根据实际勘探需要，在控制面板上输入接收器的个数n及每个接收器对应的坐标位置(l
x
,ly,lz)。假设接收器的坐标位置范围为{xi≤x≤x
i 1
；yj≤y≤y
j 1
；zk≤z≤z
k 1
}，搜索并记录距离接收器位置最近的八个测量点各方向的电场，八个测量点可认为是接收器所处长方体单元的各个顶点，其中沿x方向的电场为处长方体单元的各个顶点，其中沿x方向的电场为和
沿y方向的电场沿y方向的电场和和沿z方向的电场和并记录下对应的坐标位置(i,j,k)，(i 1,j,k)，(i,j 1,k)，(i,j,k 1)，(i,j 1,k 1)，(i 1,j,k 1)，(i 1,j 1,k)和(i 1,j 1,k 1)，然后分别存储上述测量点各个方向的电场于对角矩阵f
x
,fy和fz。
[0098]
步骤112，根据接收器的坐标位置和接收器所在长方体单元的坐标，确定接收器测量点x、y和z三个方向棱边长度的权重值，并进一步确定接收器所在长方体单元各个顶点的总权重值。
[0099]
具体地，根据接收器的坐标位置和接收器所在长方体单元的坐标，确定接收器测量点x、y和z三个方向棱边长度的权重值为：
[0100][0101]
其中，α
x1
，α
x2
，α
y1
，α
y2
，α
z1
和α
z2
分别为接收器测量点x、y和z三个方向棱边长度的权重值，(l
x
,ly,lz)表示接收器的坐标位置，(xi,yi,zi)表示第i个长方体单元的坐标；进一步确定接收器所在长方体单元各个顶点的总权重值s为：
[0102][0103]
其中，λ为顶点的权重值。
[0104]
步骤114，根据电场对角矩阵和总权重值确定接收器上的电场计算量，若电场计算量在周围点电场最小值和最大值之间，则判定接收器的响应是精确的，输出接收器响应，否则，细化三维长方体模型的网格剖分并重新计算，直到能够判定接收器的响应是精确的，输出接收器响应。
[0105]
具体地，根据电场对角矩阵和总权重值确定接收器上的电场计算量为：
[0106][0107]
其中，e
rx
、e
ry
、e
rz
分别为接收器在x、y和z方向上的电场计算量。
[0108]
若e
rx
、e
ry
、e
rz
满足下式：
[0109]erxmin
≤e
rx
≤e
rxmax
且e
rymin
≤e
ry
≤e
rymax
且e
rzmin
≤e
rz
≤e
rzmax
[0110]
则判定接收器的响应是精确的，输出接收器响应；其中，e
rxmin
,e
rymin
和e
rzmin
和e
rxmax
,e
rymax
和e
rzmax
分别为接收器所在长方体单元周围各个顶点电场各个分量的最小值和最大值。
[0111]
否则，说明该模型下，无法精确得到该接收器的响应，需要细化模型剖分。
[0112]
求解周围八个测量点的平均电阻率、细化后的单元长度和该局部区域的坐标位置：
[0113][0114]
和通过坐标位置确定该局部区域{xi≤x≤x
i 1
；yi≤y≤y
i 1
；zi≤z≤z
i 1
}。
[0115]
返回步骤102,并输入上述平均电阻率、细化后的单元长度和该局部区域范围重新细化构建内部包含勘探目标的三维长方体模型，再次进行计算，直至满足条件为止。
[0116]
上述任意位置高精度接收器响应的计算方法中，通过构建内部包含勘探目标的三维长方体模型，对三维长方体模型进行网格剖分，给每个长方体单元的电阻率赋值，得到三维异常体模型；分别计算电场控制方程的系数矩阵和右端向量，确定三维异常体模型对应的电场控制方程；通过人机交互界面获取每个接收器的坐标位置，确定接收器所处的长方体单元，并进一步根据电场控制方程确定长方体单元各个顶点的电场分量，得到电场对角矩阵；根据接收器的坐标位置和接收器所在长方体单元的坐标，确定接收器测量点x、y和z三个方向棱边长度的权重值，并进一步确定接收器所在长方体单元各个顶点的总权重值；根据电场对角矩阵和总权重值确定接收器上的电场计算量，若电场计算量在周围点电场最小值和最大值之间，则判定接收器的响应是精确的，输出接收器响应，否则，细化三维长方体模型的网格剖分并重新计算，直到能够判定接收器的响应是精确的，输出接收器响应。本发明可以避免主流方法存在的目标函数复杂和参数求解繁琐等问题，同时又解决了线性插值在电阻率剧烈变化时候精度不高的问题，能够满足任意位置接收器响应的高精度计算需求。
[0117]
应该理解的是，虽然图1的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，图1中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，
而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。
[0118]
在其中一个实施例中，还包括：由接收器上的电场er即可计算得到磁场
[0119]
在另一个实施例中，如图2所示，提供了一种任意位置高精度接收器响应的计算方法，包括：构建三维地质体模型、构建及离散计算电场的控制方程，计算任意位置接收器上的响应。构建三维地质体模型包括根据勘探目标体的形状、大小及电阻率分布构建包含勘探目标的三维长方体模型，将上述模型分解成若干个小的长方体模型，获取网格剖分参数，并对每个小的单元赋值电阻率。构建及离散计算电场的控制方程包括根据背景模型，网格剖分参数，频率参数计算边界面上的电场，从而得到右端向量b，根据三维模型、网格剖分参数、频率参数得到控制方程的系数矩阵a。计算任意位置接收器上的响应包括输入接收器的个数及每个接收器对应的坐标位置，搜索距离接收器最近的八个测量点上的电场、电阻率及坐标位置，构建局部单元测量点三个方向棱边长度的权重值，进而计算接收器总的权重值，最后求解任意接收器上的响应。若模型不能满足计算需求，需要细化模型，包括输入平均电阻率、细化后的单元长度和该局部区域的坐标位置。
[0120]
在一个具体实施例中，对于如图3所示的高-低阻混合矿体模型，模拟区域范围为：x和y方向均从-120km到120km，z方向均从-100km到160km；其中空气层的高度为-100km到0km，电阻率为10
10
欧姆米；高阻体电阻率为2000欧姆米，低阻体电阻率为20欧姆米，其大小为3km
×
3km
×
3km，距离地表1km；将该模型区域剖分为80
×
80
×
70个小长方体单元，先计算了该模型下的的电场，然后输入接收器的个数81，以及每个接收器的坐标位置(在本次实例中坐标位置地表间隔200米)。为了更好的体现本发明的创新性，将本发明和传统算法做了精度和计算时间的比较。传统算法为国际通用的插值算法；国家通用插值算法需要时间8秒，本发明算法需要5秒；利用本发明的算法，在保证计算精度的同时，也在一定程度提高了效率。
[0121]
本发明实施的任意位置高精度的接收器其响应计算方法利用matlab语言编程实现，运行程序所用的个人电脑配置为：cpu-inter core i7-8700，主频为3.4ghz，运行内存为36gb。图4为本发明一个实施例中任意位置高精度的接收器其响应计算方法计算的电场的参考解、数值解以及相对误差等值线图，其中(a)为电场的实部的参考解，(b)为电场的实部的数值解，(c)为电场的实部的相对误差；(d)为电场的虚部的参考解，(e)为电场的虚部的数值解，(f)为电场的虚部的相对误差，从图3中可以看出本发明的数值解和参考解高度一致，从相对误差图中可以看出实部最大相对误差为0.4％，虚部最大误差为0.5％，可以看出本发明方法精度很高。
[0122]
在一个实施例中，如图5所示，提供了一种任意位置高精度接收器响应的计算装置，包括：模型构建模块502、电场控制方程确定模块504、电场对角矩阵确定模块506、总权重值确定模块508和响应输出模块510，其中：
[0123]
模型构建模块502，用于构建内部包含勘探目标的三维长方体模型，对三维长方体模型进行网格剖分，形成若干个长方体单元，给每个长方体单元的电阻率赋值，得到三维异常体模型；
[0124]
电场控制方程确定模块504，用于根据三维长方体模型网格剖分参数和预设的频
率参数得到电场控制方程的系数矩阵；根据三维长方体模型网格剖分参数和频率参数得到三维异常体模型边界面上的电场，根据三维大地电磁场源和边界条件得到电场控制方程的右端向量；根据系数矩阵和右端向量确定三维异常体模型对应的电场控制方程；
[0125]
电场对角矩阵确定模块506，用于以三维长方体模型为坐标原点进行坐标转换，确定每个长方体单元的坐标，通过人机交互界面获取每个接收器的坐标位置，根据每个长方体单元的坐标和接收器的坐标位置确定接收器所处的长方体单元，并进一步根据电场控制方程确定长方体单元各个顶点的电场分量，得到电场对角矩阵；
[0126]
总权重值确定模块508，用于根据接收器的坐标位置和接收器所在长方体单元的坐标，确定接收器测量点x、y和z三个方向棱边长度的权重值，并进一步确定接收器所在长方体单元各个顶点的总权重值；
[0127]
响应输出模块510，用于根据电场对角矩阵和总权重值确定接收器上的电场计算量，若电场计算量在周围点电场最小值和最大值之间，则判定接收器的响应是精确的，输出接收器响应，否则，细化三维长方体模型的网格剖分并重新计算，直到能够判定接收器的响应是精确的，输出接收器响应。
[0128]
电场控制方程确定模块504还用于根据双旋度算子和预设的频率参数得到电场控制方程的系数矩阵为：
[0129][0130]
其中，a表示系数矩阵，表示旋度算子，表示双旋度算子，ω表示角频率，通过ω＝2πf求取，f表示给定的频率，μ表示磁导率，其值为4π
×
10-7
，ρ表示每个长方体单元的电阻率。
[0131]
电场控制方程确定模块504还用于当三维大地电磁场源x极化时，此时前侧面和后侧面边界上的电场取0，通过二维正演控制方程求取三维异常体模型所在区域的左侧面、右侧面的边界条件，通过插值左右两测边界最底下的电场值得到底面边界条件；其中，表示对应长方体单元y方向长度值δy，表示对应长方体单元z方向长度值δz，e
x
表示三维大地电磁场源在空气层顶端的电场值，e
x
＝1；当三维大地电磁场源y极化时，此时左侧面和右侧面边界上的电场取0，通过二维正演控制方程求取三维异常体模型所在区域的前侧面、后侧面的边界条件，通过插值左右两测边界最底下的电场值得到底面边界条件；其中，表示对应长方体单元x方向长度值δx；由三维大地电磁场源和边界条件构成电场控制方程的右端向量。
[0132]
电场控制方程确定模块504还用于根据系数矩阵和右端向量确定三维异常体模型对应的电场控制方程为：
[0133]
ae＝b
[0134]
其中，b表示右端向量，e表示电场。
[0135]
电场对角矩阵确定模块506还用于确定长方体单元各个顶点的电场分量，其中，沿
x方向的电场为x方向的电场为和沿y方向的电场沿y方向的电场和沿z方向的电场和得到电场对角矩阵f
x
,fy,fz。
[0136]
总权重值确定模块508还用于根据接收器的坐标位置和接收器所在长方体单元的坐标，确定接收器测量点x、y和z三个方向棱边长度的权重值为：
[0137][0138]
其中，α
x1
，α
x2
，α
y1
，α
y2
，α
z1
和α
z2
分别为接收器测量点x、y和z三个方向棱边长度的权重值，(l
x
,ly,lz)表示接收器的坐标位置，(xi,yi,zi)表示第i个长方体单元的坐标；
[0139]
进一步确定接收器所在长方体单元各个顶点的总权重值s为：
[0140][0141]
其中，λ为顶点的权重值。
[0142]
响应输出模块510还用于根据电场对角矩阵和总权重值确定接收器上的电场计算量为：
[0143][0144]
其中，e
rx
、e
ry
、e
rz
分别为接收器在x、y和z方向上的电场计算量。
[0145]
响应输出模块510还用于若e
rx
、e
ry
、e
rz
满足下式：
[0146]erxmin
≤e
rx
≤e
rxmax
且e
rymin
≤e
ry
≤e
rymax
且e
rzmin
≤e
rz
≤e
rzmax
[0147]
则判定接收器的响应是精确的，输出接收器响应；其中，e
rxmin
,e
rymin
和e
rzmin
和e
rxmax
,e
rymax
和e
rzmax
分别为接收器所在长方体单元周围各个顶点电场各个分量的最小值和
最大值。
[0148]
关于任意位置高精度接收器响应的计算装置的具体限定可以参见上文中对于任意位置高精度接收器响应的计算方法的限定，在此不再赘述。上述任意位置高精度接收器响应的计算装置中的各个模块可全部或部分通过软件、硬件及其组合来实现。上述各模块可以硬件形式内嵌于或独立于计算机设备中的处理器中，也可以以软件形式存储于计算机设备中的存储器中，以便于处理器调用执行以上各个模块对应的操作。
[0149]
在一个实施例中，提供了一种计算机设备，该计算机设备可以是终端，其内部结构图可以如图6所示。该计算机设备包括通过系统总线连接的处理器、存储器、网络接口、显示屏和输入装置。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现一种任意位置高精度接收器响应的计算方法。该计算机设备的显示屏可以是液晶显示屏或者电子墨水显示屏，该计算机设备的输入装置可以是显示屏上覆盖的触摸层，也可以是计算机设备外壳上设置的按键、轨迹球或触控板，还可以是外接的键盘、触控板或鼠标等。
[0150]
本领域技术人员可以理解，图6中示出的结构，仅仅是与本技术方案相关的部分结构的框图，并不构成对本技术方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。
[0151]
在一个实施例中，提供了一种计算机设备，包括存储器和处理器，该存储器存储有计算机程序，该处理器执行计算机程序时实现上述方法实施例中的步骤。
[0152]
在一个实施例中，提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述方法实施例中的步骤。
[0153]
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本技术所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器(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)等。
[0154]
以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
[0155]
以上所述实施例仅表达了本技术的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本技术构思的前提下，还可以做出若干变形和改进，这些都属于本技术的保护
范围。因此，本技术专利的保护范围应以所附权利要求为准。