本发明公开了一种三维直流电阻率法数值模拟方法，通过设置直流电法勘探区域和背景电阻率参数、三维异常体展布范围和电阻率参数、电流密度计算、二维离散傅里叶变换、波数域异常电位计算、二维反傅里叶变换、迭代收敛判断等步骤，实现了三维直流电法勘探高效、高精度的数值模拟。解决了目前直流电法勘探数值模拟中计算大规模模型时数据量大、存储要求高、计算时间慢，导致无法满足大规模直流电法数据精细反演的问题，有助于实现野外复杂地形或地下复杂结构下直流电法实测数据的精细反演和解释。
1.一种三维直流电阻率法数值模拟方法，其特征在于，包括如下步骤：步骤1，根据地下结构电阻率分布情况，构建目标区域，使三维电阻率异常体完全包含在该目标区域，并将目标区域沿x、y、z方向剖分为若干空间单元；
步骤2，对目标区域中剖分节点的电阻率进行赋值，并得到各个剖分节点的背景电场与背景电位，并将背景电场作为初始电场；
步骤3，基于各个剖分节点的电阻率与初始电场得到各个剖分节点的散射电流，进一步得到各个剖分节点的空间波数域散射电流；
步骤4，基于三维电阻率异常体内各个剖分节点的空间波数域散射电流，得到观测面的空间波数域异常电位，进一步得到观测面的空间波数域异常电场；
步骤5，基于观测面的空间波数域异常电场得到空间域异常电场，进一步得到观测面当前的总电场；
步骤6，判断观测面当前的总电场是否满足收敛条件：若是则将观测面的空间波数域异常电位进行二维傅里叶变换后与背景电位累加，并将累加结果作为输出值；
否则将观测面当前的总电场作为初始电场，并重复步骤3‑6；
步骤4中，所述基于三维电阻率异常体内各个剖分节点的空间波数域散射电流，得到观测面的空间波数域异常电位，具体为：将三维电阻率异常体的电流密度使用电流密度的二次插值函数表示，计算三维电阻率异常体内第m个空间单元在观测面 下的空间波数域积分，其中，第m个空间单元即第i至第 i+2个剖分节点：
式中， 、 、 分别表示观测面为 时，三维电阻率异常体内第m个空间单元积分在x，y，z方向上的空间波数域积分； 表示空间波数域的坐标，k表示波数， ， 、 分别表示x、y方向的波数， 、 、 表示第i、i+1、i+2个剖分节点在x方向上的空间波数域电流密度， 、 、 表示第i、i+1、i+2个剖分节点在y方向上的空间波数域电流密度， 、 、 表示第i、i+1、i+2个剖分节点在z方向上的空间波数域电流密度； 表示第m个空间单元中三个插值点，即第i、i+1、i+2个剖分节点的二次形函数， 为背景电阻率，表示虚数单位； 和 分别表示第i和i+2个剖分节点z方向的长度范围，另外有：；
将 、 、 累加，得到累加结果即为第m个空间单元在观测面 下的异常电位；
将地下异常区域每个空间单元在观测面 下的异常电位累加，即得到观测面的空间波数域异常电位：
2.根据权利要求1所述三维直流电阻率法数值模拟方法，其特征在于，步骤1中，所述将目标区域沿x、y、z方向剖分为若干空间单元，具体为：建立三维直角坐标系，并将坐标原点置于目标区域中心位置；
对目标区域沿x、y、z方向剖分，将该区域分为若干空间单元，对应的x、y、z方向的剖分节点数量分布为Nx、Ny、Nz；
其中，沿x、y方向空间单元间距相同，沿z方向空间单元间距的大小可根据目标区域不同深度电流密度的变化灵活剖分，对电流密度变化大的区域加密剖分，电流密度变化小的区域稀疏剖分。
3.根据权利要求2所述三维直流电阻率法数值模拟方法，其特征在于，步骤2具体包括：对目标区域中剖分节点的电阻率进行赋值，其中，存在三维电阻率异常体区域的剖分节点赋予电阻率 ，其他区域的剖分节点赋予背景电阻率 ；
基于地下供电的电流大小I和背景电阻率 ，计算各剖分节点在x，y，z方向的背景电场三分量 、 、 以及背景电位 ，其中， 表示编号为 的剖分节点坐标， ， ， ；
将背景电场作为初始电场，即： 、 、，其中， 、 、 为初始电场在x，y，z方向的三分量。
4.根据权利要求3所述三维直流电阻率法数值模拟方法，其特征在于，步骤3具体包括：基于各个剖分节点的电阻率与初始电场得到各个剖分节点的散射电流，为：式中， 、 、 表示剖分节点坐标 的散射电流在x，y，z方向的三分量；
采用二维傅里叶变换算法，计算各个剖分节点的空间波数域散射电流，为：式中， 、 、 表示剖分节点 的空间波数域散射电流在x，y，z方向的三分量， 、 分别表示x、y方向的波数， 表示电流密度在空间波数域的坐标，表示虚数单位。
5.根据权利要求4所述三维直流电阻率法数值模拟方法，其特征在于，步骤4中，所述得到观测面的空间波数域异常电场，具体为：利用二维傅立叶变换的性质，基于观测面 的空间波数域异常电位得到空间域异常电场，为：
式中， 、 为观测面的空间波数域异常电场在x，y，z方向的三分量。
6.根据权利要求5所述三维直流电阻率法数值模拟方法，其特征在于，步骤5具体为：对观测面的波数域异常电场 进行二维反傅里叶变换，得到观测面的空间域异常电场 ，为
基于观测面的空间域异常电场与背景电场得到观测面当前的总电场，为：式中， 为观测面当前的总电场。
7.根据权利要求6所述三维直流电阻率法数值模拟方法，其特征在于，步骤6中，所述收敛条件为：
式中， 为期望收敛时的数值精度， 为观测面的初始电场。
一种三维直流电阻率法数值模拟方法
技术领域
[0001] 本发明涉及数值模拟技术领域，具体是一种三维直流电阻率法数值模拟方法。
背景技术
[0002] 直流电阻率法作为地球物理勘探基础方法之一，因其仪器操作简单、高效、高精度辨识、探测深度大、获取地下信息全面等优势而被广泛用于矿产资源勘探、工程环境勘探、
考古、深海探测等方面。随着野外勘探深度和难度的增加，在复杂地形和地质结构条件下实
现可靠、高效、高精度的直流电阻率法勘探成为研究的重点。正演是反演的基础，正演模拟
的计算精度和效率决定着直流电阻率法勘探反演的效率与精度。
[0003] 目前已有的直流电阻率法数值模拟方法中，为了在复杂条件下获得较高的计算精度，需要对模型进行精细剖分，导致计算量和存储量剧增，计算时间长，无法兼顾计算精度
和计算效率，难以实现大规模复杂模型的数值模拟。因此，寻找一种能在复杂条件下兼顾计
算精度和计算效率的数值模拟方法，对实现直流电阻率法高效反演成像具有重要的现实意
义。
发明内容
[0004] 为了解决常规方法对复杂模型数值模拟计算时无法兼顾计算精度和计算效率，无法满足直流电阻率法大规模实测数据精细反演成像的问题，本发明提供一种三维直流电阻
率法数值模拟方法，有助于实现野外复杂地形或地下复杂结果下直流电法实测数据的精细
反演和解释。
[0005] 为实现上述目的，本发明提供一种三维直流电阻率法数值模拟方法，包括如下步骤：
[0006] 步骤1，根据地下结构电阻率分布情况，构建目标区域，使三维电阻率异常体完全包含在该目标区域，并将目标区域沿x、y、z方向剖分为若干空间单元；
[0007] 步骤2，对目标区域中剖分节点的电阻率进行赋值，并得到各个剖分节点的背景电场与背景电位，并将背景电场作为初始电场；
[0008] 步骤3，基于各个剖分节点的电阻率与初始电场得到各个剖分节点的散射电流，进一步得到各个剖分节点的空间波数域散射电流；
[0009] 步骤4，基于三维电阻率异常体内各个剖分节点的空间波数域散射电流，得到观测面的空间波数域异常电位，进一步得到观测面的空间波数域异常电场；
[0010] 步骤5，基于观测面的空间波数域异常电场得到空间域异常电场，进一步得到观测面当前的总电场；
[0011] 步骤6，判断观测面当前的总电场是否满足收敛条件：
[0012] 若是则将观测面的波数域异常电位进行二维傅里叶变换后与背景电位累加，并将累加结果作为输出值；
[0013] 否则将观测面当前的总电场作为初始电场，并重复步骤3‑6。
[0014] 在其中一个实施例中，步骤1中，所述将目标区域沿x、y、z方向剖分为若干空间单元，具体为：
[0015] 建立三维直角坐标系，并将坐标原点置于目标区域中心位置；
[0016] 对目标区域沿x、y、z方向剖分，将该区域分为若干空间单元，x、y、z方向的剖分节点数量分布为Nx、Ny、Nz；
[0017] 其中，沿x、y方向空间单元间距相同，沿z方向空间单元间距的大小可根据目标区域不同深度电流密度的变化灵活剖分，对电流密度变化大的区域加密剖分，电流密度变化
小的区域稀疏剖分。
[0018] 在其中一个实施例中，步骤2具体包括：
[0019] 对目标区域中剖分节点的电阻率进行赋值，其中，存在三维电阻率异常体区域的剖分节点赋予电阻率 ，其他区域的剖分节点赋予背景电阻率 ；
[0020] 基于地下供电的电流大小I和背景电阻率 ，计算各剖分节点在x，y，z方向的背景电场 三分量 、 、 以及背景电位
，其中， 表示编号为 的剖分节点坐标， ，
， ；
[0021] 将 背 景 电 场 作 为 初 始 电 场 ，即 ： 、、 ，其中， 、
、 为初始电场 在x，y，z方向的三分量。
[0022] 在其中一个实施例中，步骤3具体包括：
[0023] 基于各个剖分节点的电阻率与初始电场得到各个剖分节点的散射电流，为：
[0024]
[0025]
[0026]
[0027] 式中， 、 、 表示剖分节点 的散射电流在x，y，z方向的三分量；
[0028] 采用二维傅里叶变换算法，计算各个剖分节点的空间波数域散射电流，为：
[0029]
[0030]
[0031]
[0032] 式中， 、 、 表示剖分节点 的空间波数域散射电流在x，y，z方向的三分量， 、 分别表示x、y方向的波数， 表示电
流密度在空间波数域的坐标，表示虚数单位。
[0033] 在其中一个实施例中，步骤4中，所述基于三维电阻率异常体内各个剖分节点的空间波数域散射电流，得到观测面的空间波数域异常电位，具体为：
[0034] 将三维电阻率异常体的电流密度使用电流密度的二次插值函数表示，计算三维电阻率异常体内第m个空间单元在观测面 下的空间波数域积分，其中，第m个空间单元即第i
至第 i+2个剖分节点：
[0035]
[0036]
[0037]
[0038] 式中， 、 、 分别表示观测面为 时，三维电阻率异常体内第m个空间单元积分在x，y，z方向上的空间波数域积分； 表示空间
波数域的坐标，k表示波数， 、 表示第i、i+1、i+2个剖分节点在x
方向上的空间波数域电流密度， 、 、 表示第i、i+1、i+2个剖分节点在y方向上的空
间波数域电流密度， 、 、 表示第i、i+1、i+2个剖分节点在z方向上的空间波数域电
流密度； 表示第m个空间单元中三个插值点，即第i、i+1、i+2个剖分节点的二次
形函数， 和 分别 第i和i+2个剖分节点z方向的长度范围，另外有：
[0039] ；
[0040] 将 、 、 累加，得到累加结果即为第m个空间单元在观测面 下的异常电位；
[0041] 将地下异常区域每个空间单元在观测面 下的异常电位累加，即得到观测面的空间波数域异常电位：
[0042]
[0043] 式中， 为观测面 的空间波数域异常电位。
[0044] 在其中一个实施例中，步骤4中，所述得到观测面的空间波数域异常电场，具体为：
[0045] 利用二维傅立叶变换的性质，基于观测面的空间波数域异常电位得到空间域异常电场，为：
[0046]
[0047]
[0048]
[0049] 式中， 、 为观测面的空间波数域异常电场 在x，y，z方向的三分量，表示虚数单位。
[0050] 在其中一个实施例中，步骤5具体为：
[0051] 对观测面的波数域异常电场 进行二维反傅里叶变换，得到观测面的空间域异常电场 ，为
[0052]
[0053]
[0054]
[0055] 基于观测面的空间域异常电场与背景电场得到观测面当前的总电场，为：
[0056]
[0057] 式中， 为观测面当前的总电场。
[0058] 在其中一个实施例中，步骤6中，所述收敛条件为：
[0059]
[0060] 式中，为期望收敛时的数值精度， 为观测面的初始电场。
[0061] 与现有的直流电阻率法数值技术相比，本发明提供的一种三维直流电阻率法数值模拟方法具有如下有益技术效果：
[0062] （1）将一个三维问题分解成多个一维问题进行求解，在计算大规模复杂形态时，可有效减少计算量和存储需要，计算速度快；
[0063] （2）垂直方向可根据地下电流密度变化灵活剖分，电流密度变化大的区域加密剖分，电流密度变化小的区域稀疏剖分，这样能很好地刻画地下复杂结构和地形，计算精度
高；
[0064] （3）实现了直流电阻率法快速、高精度数值模拟，可满足大规模实测数据三维精细反演的需求。
附图说明
[0065] 为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本
发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以
根据这些附图示出的结构获得其他的附图。
[0066] 图1为本发明实施例中三维直流电阻率法数值模拟方法的流程示意图；
[0067] 图2为本发明实施例中沿z方向空间单元灵活剖分的示意图；
[0068] 图3为本发明实施例中示例的目标区域示意图；
[0069] 图4为本发明实施例中示例的磁场计算值和理论值对比图；
[0070] 图5为本发明实施例中示例的磁场计算值和理论值绝对误差图。
[0071] 本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。
具体实施方式
[0072] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基
于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其
他实施例，都属于本发明保护的范围。
[0073] 需要说明，本发明实施例中所有方向性指示（诸如上、下、左、右、前、后……）仅用于解释在某一特定姿态（如附图所示）下各部件之间的相对位置关系、运动情况等，如果该
特定姿态发生改变时，则该方向性指示也相应地随之改变。
[0074] 另外，在本发明中如涉及“第一”、“第二”等的描述仅用于描述目的，而不能理解为指示或暗示其相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、
“第二”的特征可以明示或者隐含地包括至少一个该特征。在本发明的描述中，“多个”的含
义是至少两个，例如两个，三个等，除非另有明确具体的限定。
[0075] 在本发明中，除非另有明确的规定和限定，术语“连接”、“固定”等应做广义理解，例如，“固定”可以是固定连接，也可以是可拆卸连接，或成一体；可以是机械连接，也可以是
电连接，还可以是物理连接或无线通信连接；可以是直接相连，也可以通过中间媒介间接相
连，可以是两个元件内部的连通或两个元件的相互作用关系，除非另有明确的限定。对于本
领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。
[0076] 另外，本发明各个实施例之间的技术方案可以相互结合，但是必须是以本领域普通技术人员能够实现为基础，当技术方案的结合出现相互矛盾或无法实现时应当认为这种
技术方案的结合不存在，也不在本发明要求的保护范围之内。
[0077] 如图1所示为本实施例公开的一种三维直流电阻率法数值模拟方法，包括如下步骤：
[0078] 步骤1，根据地下结构电阻率分布情况，构建目标区域，使三维电阻率异常体完全包含在该目标区域，并将目标区域沿x、y、z方向剖分为若干空间单元。具体为：
[0079] 建立三维直角坐标系，并将坐标原点置于目标区域中心位置；
[0080] 对目标区域沿x、y、z方向剖分，将该区域分为若干空间单元，x、y、z方向的剖分节点数量分布为Nx、Ny、Nz；
[0081] 其中，沿x、y方向空间单元间距相同，沿z方向空间单元间距的大小可根据目标区域不同深度电流密度的变化灵活剖分。
[0082] 本实施例中，沿z方向空间单元间距的大小可根据目标区域不同深度电流密度的变化灵活剖分的具体实施过程为：电流密度变化大的区域加密剖分，电流密度变化小的区
域稀疏剖分，如图2所示。
[0083] 步骤2，对目标区域中剖分节点的电阻率进行赋值，并得到各个剖分节点的背景电场与背景电位，并将背景电场作为初始电场。其具体实施过程为：
[0084] 对目标区域中剖分节点的电阻率进行赋值，其中，存在三维电阻率异常体区域的剖分节点赋予电阻率 ，其他区域的剖分节点赋予 ；
[0085] 基于地下供电的电流大小I和背景电阻率 ，计算各剖分节点在x，y，z方向的背景电场 三分量 、 、 ，为：
[0086]
[0087]
[0088]
[0089] 以及背景电位 ，为：
[0090]
[0091] 其中， 表示编号为 的剖分节点坐标， ，， ；表示点源供电的电流大小；表示圆周率； 表示编号为 剖
分节点到点源 的距离，为：
[0092]
[0093] 式中， 表示编号为 的剖分节点坐标， 表示点源在x，y，z方向上的坐标。将背景电 场与背景电位作为初始电场与初始电位 ，即：
、 、 ，
其中， 、 、 为初始电场 在x，y，z方向的三分量。
[0094] 步骤3，基于各个剖分节点的电阻率与初始电场得到各个剖分节点的散射电流，进一步得到各个剖分节点的波数域散射电流。其具体实施过程为：
[0095] 基于各个剖分节点的电阻率与初始电场得到各个剖分节点的散射电流，为：
[0096]
[0097]
[0098]
[0099] 式中， 、 、 表示剖分节点 的散射电流在x，y，z方向的三分量；
[0100] 采用二维傅里叶变换算法，计算各个剖分节点的波数域散射电流，为：
[0101]
[0102]
[0103]
[0104] 式中， 、 、 表示剖分节点 的波数域散射电流在x，y，z方向的三分量， 、 分别表示x、y方向的波数， 表示电流密度
在波数域的坐标，表示虚数单位。
[0105] 步骤4，基于各个剖分节点的波数域散射电流，得到观测面的波数域异常电位，进一步得到观测面的空间波数域异常电场。其具体实施过程为：
[0106] 将三维电阻率异常体的电流密度使用电流密度的二次插值函数表示，计算异常体第m个空间单元（即第i至第 i+2个剖分节点）在观测面 下的空间波数域积分：
[0107] 下的空间波数域积分：
[0108]
[0109]
[0110]
[0111] 式中， 、 、 分别表示观测面为 时，三维电阻率异常体内第m个空间单元积分在x，y，z方向上的空间波数域积分； 表示空间
波数域的坐标，k表示波数， 、 表示第i、i+1、i+2个剖分节点在x
方向上的空间波数域电流密度， 、 、 表示第i、i+1、i+2个剖分节点在y方向上的空间
波数域电流密度， 、 、 表示第i、i+1、i+2个剖分节点在z方向上的空间波数域电流
密度。 表示第m个空间单元中三个插值点，即第i、i+1、i+2个剖分节点的二次
形函数， 和 分别表示第i和i+2个剖分节点z方向的长度范围，另外有：
[0112] ；
[0113] 将 、 、 累加，得到累加结果即为第m个空间单元在观测面 下的异常电位；
[0114] 将地下异常区域每个空间单元在观测面 下的异常电位累加，即得到观测面的空间波数域异常电场：
[0115]
[0116] 式中， 为观测面的空间波数域异常电位。
[0117] 本实施例中，观测面的空间波数域异常电场的获取过程具体为：
[0118] 利用二维傅立叶变换的性质，基于观测面 的空间波数域异常电位得到空间域异常电场，为：
[0119]
[0120]
[0121]
[0122] 式中， 、 为观测面的空间波数域异常电场 在x，y，z方向的三分量，表示虚数单位。
[0123] 步骤5，基于观测面的空间波数域异常电位得到空间域异常电场，进一步得到观测面当前的总电场，其具体实施过程为：
[0124] 对观测面的波数域异常电场 进行二维傅里叶反变换，得到观测面的空间域异常电位 ，为
[0125]
[0126]
[0127]
[0128] 基于观测面的空间域异常电场与背景电场得到观测面当前的总电场，为：
[0129]
[0130] 式中， 为观测面当前的总电场。
[0131] 步骤6，判断观测面当前的总电场是否满足以下收敛条件：
[0132]
[0133] 式中，为期望收敛时的数值精度， 为观测面的初始电场；
[0134] 若满足，则将观测面的波数域异常电位 进行二维反傅里叶变换得到，将 与背景电位 累加，并将累加结果
作为输出值；
[0135] 否不满足，将观测面当前的总电场 作为初始电场 ，即令，重复步骤3‑6，直至观测面当前的总电场满足收敛条件。
[0136] 下面结合具体的示例，验证本实施例所提出的方法用于计算三维电阻体的电位计算方法的效果进行验证。
[0137] 目标区域存有一低阻异常体，如图3所示。异常体沿x，y方向长度为1000m，沿z方向250m，顶面埋深50m，电阻率为 ，背景电阻率为 。令异常体中心位置在地面
的投影为坐标原点，供电电极A和B坐标分别为(‑1000m,0, 0)和(1000m,0,0)，电流大小为
1A。模型空间单元剖分节点数为101×101×51，x，y方向剖分间隔都为10m，z方向剖分间隔
为5m。采用中梯装置测量，MN=20m。期望收敛的相对残差为 。
[0138] 本实施例方法利用Fortran语言编程实现，运行程序所用的个人台式机配置为：CPU‑Inter Core i5‑4590，主频为3.3GHz，内存为12GB。本实施例的三维电阻率低阻体的电
位方法计算达到期望收敛相对残差为10‑4时，需要迭代8次，计算总时间为6.62s，由此可见
本实施例具有很高的计算效率。图4为沿着x方向从‑500m至500m测量得到的视电阻率和相
对误差图，具有本实施例计算得到的视电阻图与使用美国犹他大学IETEM3D软件计算得到
的视电阻率图，图5为图4中两种算法之间的相对误差图。从图5中可看出两种算法的相对误
差小于0.35%，验证了本实施例方法的正确性，说明本实施例方法具有很高的计算精度。
[0139] 以上所述仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是在本发明的发明构思下，利用本发明说明书及附图内容所作的等效结构变换，或直接/间接运用
在其他相关的技术领域均包括在本发明的专利保护范围内。
