CN113051779B - 一种三维直流电阻率法数值模拟方法 - Google Patents

一种三维直流电阻率法数值模拟方法 Download PDF

Info

Publication number
CN113051779B
CN113051779B CN202110596835.8A CN202110596835A CN113051779B CN 113051779 B CN113051779 B CN 113051779B CN 202110596835 A CN202110596835 A CN 202110596835A CN 113051779 B CN113051779 B CN 113051779B
Authority
CN
China
Prior art keywords
electric field
resistivity
current
dimensional
space
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN202110596835.8A
Other languages
English (en)
Other versions
CN113051779A (zh
Inventor
凌嘉宣
戴世坤
陈轻蕊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN202110596835.8A priority Critical patent/CN113051779B/zh
Publication of CN113051779A publication Critical patent/CN113051779A/zh
Application granted granted Critical
Publication of CN113051779B publication Critical patent/CN113051779B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Measurement Of Resistance Or Impedance (AREA)

Abstract

本发明公开了一种三维直流电阻率法数值模拟方法，通过设置直流电法勘探区域和背景电阻率参数、三维异常体展布范围和电阻率参数、电流密度计算、二维离散傅里叶变换、波数域异常电位计算、二维反傅里叶变换、迭代收敛判断等步骤，实现了三维直流电法勘探高效、高精度的数值模拟。解决了目前直流电法勘探数值模拟中计算大规模模型时数据量大、存储要求高、计算时间慢，导致无法满足大规模直流电法数据精细反演的问题，有助于实现野外复杂地形或地下复杂结构下直流电法实测数据的精细反演和解释。

Description

一种三维直流电阻率法数值模拟方法
技术领域
本发明涉及数值模拟技术领域，具体是一种三维直流电阻率法数值模拟方法。
背景技术
直流电阻率法作为地球物理勘探基础方法之一，因其仪器操作简单、高效、高精度辨识、探测深度大、获取地下信息全面等优势而被广泛用于矿产资源勘探、工程环境勘探、考古、深海探测等方面。随着野外勘探深度和难度的增加，在复杂地形和地质结构条件下实现可靠、高效、高精度的直流电阻率法勘探成为研究的重点。正演是反演的基础，正演模拟的计算精度和效率决定着直流电阻率法勘探反演的效率与精度。
目前已有的直流电阻率法数值模拟方法中，为了在复杂条件下获得较高的计算精度，需要对模型进行精细剖分，导致计算量和存储量剧增，计算时间长，无法兼顾计算精度和计算效率，难以实现大规模复杂模型的数值模拟。因此，寻找一种能在复杂条件下兼顾计算精度和计算效率的数值模拟方法，对实现直流电阻率法高效反演成像具有重要的现实意义。
发明内容
为了解决常规方法对复杂模型数值模拟计算时无法兼顾计算精度和计算效率，无法满足直流电阻率法大规模实测数据精细反演成像的问题，本发明提供一种三维直流电阻率法数值模拟方法，有助于实现野外复杂地形或地下复杂结果下直流电法实测数据的精细反演和解释。
为实现上述目的，本发明提供一种三维直流电阻率法数值模拟方法，包括如下步骤：
步骤1，根据地下结构电阻率分布情况，构建目标区域，使三维电阻率异常体完全包含在该目标区域，并将目标区域沿x、y、z方向剖分为若干空间单元；
步骤2，对目标区域中剖分节点的电阻率进行赋值，并得到各个剖分节点的背景电场与背景电位，并将背景电场作为初始电场；
步骤3，基于各个剖分节点的电阻率与初始电场得到各个剖分节点的散射电流，进一步得到各个剖分节点的空间波数域散射电流；
步骤4，基于三维电阻率异常体内各个剖分节点的空间波数域散射电流，得到观测面的空间波数域异常电位，进一步得到观测面的空间波数域异常电场；
步骤5，基于观测面的空间波数域异常电场得到空间域异常电场，进一步得到观测面当前的总电场；
步骤6，判断观测面当前的总电场是否满足收敛条件：
若是则将观测面的波数域异常电位进行二维傅里叶变换后与背景电位累加，并将累加结果作为输出值；
否则将观测面当前的总电场作为初始电场，并重复步骤3-6。
在其中一个实施例中，步骤1中，所述将目标区域沿x、y、z方向剖分为若干空间单元，具体为：
建立三维直角坐标系，并将坐标原点置于目标区域中心位置；
对目标区域沿x、y、z方向剖分，将该区域分为若干空间单元，x、y、z方向的剖分节点数量分布为Nx、Ny、Nz；
其中，沿x、y方向空间单元间距相同，沿z方向空间单元间距的大小可根据目标区域不同深度电流密度的变化灵活剖分，对电流密度变化大的区域加密剖分，电流密度变化小的区域稀疏剖分。
在其中一个实施例中，步骤2具体包括：
对目标区域中剖分节点的电阻率进行赋值，其中，存在三维电阻率异常体区域的剖分节点赋予电阻率
Figure 477446DEST_PATH_IMAGE001
，其他区域的剖分节点赋予背景电阻率
Figure 290376DEST_PATH_IMAGE002
；
基于地下供电的电流大小I和背景电阻率
Figure 831079DEST_PATH_IMAGE002
，计算各剖分节点在x，y，z方向的背景电场
Figure 943391DEST_PATH_IMAGE003
三分量
Figure 606454DEST_PATH_IMAGE004
、
Figure 600954DEST_PATH_IMAGE005
、
Figure 883031DEST_PATH_IMAGE006
以及背景电位
Figure 646588DEST_PATH_IMAGE007
，其中，
Figure 496863DEST_PATH_IMAGE008
表示编号为
Figure 713081DEST_PATH_IMAGE009
的剖分节点坐标，
Figure 64428DEST_PATH_IMAGE010
，
Figure 744808DEST_PATH_IMAGE011
，
Figure 952935DEST_PATH_IMAGE012
；
将背景电场作为初始电场，即：
Figure 594132DEST_PATH_IMAGE013
、
Figure 280329DEST_PATH_IMAGE014
、
Figure 595641DEST_PATH_IMAGE015
，其中，
Figure 912353DEST_PATH_IMAGE016
、
Figure 103163DEST_PATH_IMAGE017
、
Figure 655367DEST_PATH_IMAGE018
为初始电场
Figure 982443DEST_PATH_IMAGE019
在x，y，z方向的三分量。
在其中一个实施例中，步骤3具体包括：
基于各个剖分节点的电阻率与初始电场得到各个剖分节点的散射电流，为：
Figure 470056DEST_PATH_IMAGE020
Figure 148162DEST_PATH_IMAGE021
Figure 254790DEST_PATH_IMAGE022
式中，
Figure 905214DEST_PATH_IMAGE023
、
Figure 360466DEST_PATH_IMAGE024
、
Figure 588185DEST_PATH_IMAGE025
表示剖分节点
Figure 419875DEST_PATH_IMAGE026
的散射电流在x，y，z方向的三分量；
采用二维傅里叶变换算法，计算各个剖分节点的空间波数域散射电流，为：
Figure 924805DEST_PATH_IMAGE027
Figure 550959DEST_PATH_IMAGE028
Figure 780821DEST_PATH_IMAGE029
式中，
Figure 353884DEST_PATH_IMAGE030
、
Figure 510059DEST_PATH_IMAGE031
、
Figure 635010DEST_PATH_IMAGE032
表示剖分节点
Figure 509425DEST_PATH_IMAGE033
的空间波数域散射电流在x，y，z方向的三分量，
Figure 886180DEST_PATH_IMAGE034
、
Figure 162440DEST_PATH_IMAGE035
分别表示x、y方向的波数，
Figure 209025DEST_PATH_IMAGE036
表示电流密度在空间波数域的坐标，
Figure 570736DEST_PATH_IMAGE037
表示虚数单位。
在其中一个实施例中，步骤4中，所述基于三维电阻率异常体内各个剖分节点的空间波数域散射电流，得到观测面的空间波数域异常电位，具体为：
将三维电阻率异常体的电流密度使用电流密度的二次插值函数表示，计算三维电阻率异常体内第m个空间单元在观测面
Figure 485603DEST_PATH_IMAGE038
下的空间波数域积分，其中，第m个空间单元即第i至第 i+2个剖分节点：
Figure 944266DEST_PATH_IMAGE039
Figure 348702DEST_PATH_IMAGE040
Figure 135393DEST_PATH_IMAGE041
式中，
Figure 650688DEST_PATH_IMAGE042
、
Figure 213125DEST_PATH_IMAGE043
、
Figure 788463DEST_PATH_IMAGE044
分别表示观测面为
Figure 62449DEST_PATH_IMAGE045
时，三维电阻率异常体内第m个空间单元积分在x，y，z方向上的空间波数域积分；
Figure 709331DEST_PATH_IMAGE046
表示空间波数域的坐标，k表示波数，
Figure 814691DEST_PATH_IMAGE047
、
Figure 498613DEST_PATH_IMAGE048
表示第i、i+1、i+2个剖分节点在x方向上的空间波数域电流密度，
Figure 56633DEST_PATH_IMAGE049
、
Figure 257938DEST_PATH_IMAGE050
、
Figure 155487DEST_PATH_IMAGE051
表示第i、i+1、i+2个剖分节点在y方向上的空间波数域电流密度，
Figure 72628DEST_PATH_IMAGE052
、
Figure 445840DEST_PATH_IMAGE053
、
Figure 106629DEST_PATH_IMAGE054
表示第i、i+1、i+2个剖分节点在z方向上的空间波数域电流密度；
Figure 124263DEST_PATH_IMAGE055
表示第m个空间单元中三个插值点，即第i、i+1、i+2个剖分节点的二次形函数，
Figure 795328DEST_PATH_IMAGE056
和
Figure 327941DEST_PATH_IMAGE057
分别
Figure 730103DEST_PATH_IMAGE058
第i和i+2个剖分节点z方向的长度范围，另外有：
Figure 398982DEST_PATH_IMAGE059
；
将
Figure 251400DEST_PATH_IMAGE060
、
Figure 208992DEST_PATH_IMAGE061
、
Figure 211583DEST_PATH_IMAGE062
累加，得到累加结果即为第m个空间单元在观测面
Figure 813597DEST_PATH_IMAGE063
下的异常电位；
将地下异常区域每个空间单元在观测面
Figure 509020DEST_PATH_IMAGE063
下的异常电位累加，即得到观测面的空间波数域异常电位：
Figure 688329DEST_PATH_IMAGE064
式中，
Figure 822507DEST_PATH_IMAGE065
为观测面
Figure 731557DEST_PATH_IMAGE066
的空间波数域异常电位。
在其中一个实施例中，步骤4中，所述得到观测面的空间波数域异常电场，具体为：
利用二维傅立叶变换的性质，基于观测面的空间波数域异常电位得到空间域异常电场，为：
Figure 269986DEST_PATH_IMAGE067
Figure 264487DEST_PATH_IMAGE068
Figure 920465DEST_PATH_IMAGE069
式中，
Figure 684021DEST_PATH_IMAGE070
、
Figure 393351DEST_PATH_IMAGE071
为观测面的空间波数域异常电场
Figure 875148DEST_PATH_IMAGE072
在x，y，z方向的三分量，
Figure 85550DEST_PATH_IMAGE073
表示虚数单位。
在其中一个实施例中，步骤5具体为：
对观测面的波数域异常电场
Figure 906875DEST_PATH_IMAGE074
进行二维反傅里叶变换，得到观测面的空间域异常电场
Figure 115003DEST_PATH_IMAGE075
，为
Figure 631566DEST_PATH_IMAGE076
Figure 317762DEST_PATH_IMAGE077
Figure 993594DEST_PATH_IMAGE078
基于观测面的空间域异常电场与背景电场得到观测面当前的总电场，为：
Figure 434940DEST_PATH_IMAGE079
式中，
Figure 625750DEST_PATH_IMAGE080
为观测面当前的总电场。
在其中一个实施例中，步骤6中，所述收敛条件为：
Figure 53320DEST_PATH_IMAGE081
式中，
Figure 645975DEST_PATH_IMAGE082
为期望收敛时的数值精度，
Figure 773069DEST_PATH_IMAGE083
为观测面的初始电场。
与现有的直流电阻率法数值技术相比，本发明提供的一种三维直流电阻率法数值模拟方法具有如下有益技术效果：
（1）将一个三维问题分解成多个一维问题进行求解，在计算大规模复杂形态时，可有效减少计算量和存储需要，计算速度快；
（2）垂直方向可根据地下电流密度变化灵活剖分，电流密度变化大的区域加密剖分，电流密度变化小的区域稀疏剖分，这样能很好地刻画地下复杂结构和地形，计算精度高；
（3）实现了直流电阻率法快速、高精度数值模拟，可满足大规模实测数据三维精细反演的需求。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图示出的结构获得其他的附图。
图1为本发明实施例中三维直流电阻率法数值模拟方法的流程示意图；
图2为本发明实施例中沿z方向空间单元灵活剖分的示意图；
图3为本发明实施例中示例的目标区域示意图；
图4为本发明实施例中示例的磁场计算值和理论值对比图；
图5为本发明实施例中示例的磁场计算值和理论值绝对误差图。
本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
需要说明，本发明实施例中所有方向性指示（诸如上、下、左、右、前、后……）仅用于解释在某一特定姿态（如附图所示）下各部件之间的相对位置关系、运动情况等，如果该特定姿态发生改变时，则该方向性指示也相应地随之改变。
另外，在本发明中如涉及“第一”、“第二”等的描述仅用于描述目的，而不能理解为指示或暗示其相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括至少一个该特征。在本发明的描述中，“多个”的含义是至少两个，例如两个，三个等，除非另有明确具体的限定。
在本发明中，除非另有明确的规定和限定，术语“连接”、“固定”等应做广义理解，例如，“固定”可以是固定连接，也可以是可拆卸连接，或成一体；可以是机械连接，也可以是电连接，还可以是物理连接或无线通信连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通或两个元件的相互作用关系，除非另有明确的限定。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。
另外，本发明各个实施例之间的技术方案可以相互结合，但是必须是以本领域普通技术人员能够实现为基础，当技术方案的结合出现相互矛盾或无法实现时应当认为这种技术方案的结合不存在，也不在本发明要求的保护范围之内。
如图1所示为本实施例公开的一种三维直流电阻率法数值模拟方法，包括如下步骤：
步骤1，根据地下结构电阻率分布情况，构建目标区域，使三维电阻率异常体完全包含在该目标区域，并将目标区域沿x、y、z方向剖分为若干空间单元。具体为：
建立三维直角坐标系，并将坐标原点置于目标区域中心位置；
对目标区域沿x、y、z方向剖分，将该区域分为若干空间单元，x、y、z方向的剖分节点数量分布为Nx、Ny、Nz；
其中，沿x、y方向空间单元间距相同，沿z方向空间单元间距的大小可根据目标区域不同深度电流密度的变化灵活剖分。
本实施例中，沿z方向空间单元间距的大小可根据目标区域不同深度电流密度的变化灵活剖分的具体实施过程为：电流密度变化大的区域加密剖分，电流密度变化小的区域稀疏剖分，如图2所示。
步骤2，对目标区域中剖分节点的电阻率进行赋值，并得到各个剖分节点的背景电场与背景电位，并将背景电场作为初始电场。其具体实施过程为：
对目标区域中剖分节点的电阻率进行赋值，其中，存在三维电阻率异常体区域的剖分节点赋予电阻率
Figure 388858DEST_PATH_IMAGE084
，其他区域的剖分节点赋予
Figure 416857DEST_PATH_IMAGE085
；
基于地下供电的电流大小I和背景电阻率
Figure 926336DEST_PATH_IMAGE085
，计算各剖分节点在x，y，z方向的背景电场
Figure 647167DEST_PATH_IMAGE086
三分量
Figure 750252DEST_PATH_IMAGE087
、
Figure 581942DEST_PATH_IMAGE088
、
Figure 962239DEST_PATH_IMAGE089
，为：
Figure 791655DEST_PATH_IMAGE090
Figure 444353DEST_PATH_IMAGE091
Figure 876471DEST_PATH_IMAGE092
以及背景电位
Figure 298225DEST_PATH_IMAGE093
，为：
Figure 298542DEST_PATH_IMAGE094
其中，
Figure 750121DEST_PATH_IMAGE095
表示编号为
Figure 189193DEST_PATH_IMAGE096
的剖分节点坐标，
Figure 403136DEST_PATH_IMAGE097
，
Figure 636672DEST_PATH_IMAGE098
，
Figure 60700DEST_PATH_IMAGE099
；
Figure 975566DEST_PATH_IMAGE100
表示点源供电的电流大小；
Figure 371912DEST_PATH_IMAGE101
表示圆周率；
Figure 589398DEST_PATH_IMAGE102
表示编号为
Figure 438406DEST_PATH_IMAGE103
剖分节点到点源
Figure 891384DEST_PATH_IMAGE104
的距离，为：
Figure 142236DEST_PATH_IMAGE105
式中，
Figure 779891DEST_PATH_IMAGE106
表示编号为
Figure 116195DEST_PATH_IMAGE107
的剖分节点坐标，
Figure 372864DEST_PATH_IMAGE108
表示点源在x，y，z方向上的坐标。将背景电场与背景电位作为初始电场与初始电位，即：
Figure 72965DEST_PATH_IMAGE109
、
Figure 819204DEST_PATH_IMAGE110
、
Figure 314907DEST_PATH_IMAGE111
，其中，
Figure 437584DEST_PATH_IMAGE112
、
Figure 459767DEST_PATH_IMAGE113
、
Figure 376907DEST_PATH_IMAGE114
为初始电场
Figure 891065DEST_PATH_IMAGE115
在x，y，z方向的三分量。
步骤3，基于各个剖分节点的电阻率与初始电场得到各个剖分节点的散射电流，进一步得到各个剖分节点的波数域散射电流。其具体实施过程为：
基于各个剖分节点的电阻率与初始电场得到各个剖分节点的散射电流，为：
Figure 364903DEST_PATH_IMAGE116
Figure 179275DEST_PATH_IMAGE117
Figure 470579DEST_PATH_IMAGE118
式中，
Figure 3192DEST_PATH_IMAGE119
、
Figure 264409DEST_PATH_IMAGE120
、
Figure 198867DEST_PATH_IMAGE121
表示剖分节点
Figure 926651DEST_PATH_IMAGE122
的散射电流在x，y，z方向的三分量；
采用二维傅里叶变换算法，计算各个剖分节点的波数域散射电流，为：
Figure 992565DEST_PATH_IMAGE123
Figure 260735DEST_PATH_IMAGE124
Figure 987383DEST_PATH_IMAGE125
式中，
Figure 417227DEST_PATH_IMAGE126
、
Figure 986749DEST_PATH_IMAGE127
、
Figure 793031DEST_PATH_IMAGE128
表示剖分节点
Figure 639764DEST_PATH_IMAGE129
的波数域散射电流在x，y，z方向的三分量，
Figure 319138DEST_PATH_IMAGE130
、
Figure 48060DEST_PATH_IMAGE131
分别表示x、y方向的波数，
Figure 595716DEST_PATH_IMAGE132
表示电流密度在波数域的坐标，
Figure 359272DEST_PATH_IMAGE133
表示虚数单位。
步骤4，基于各个剖分节点的波数域散射电流，得到观测面的波数域异常电位，进一步得到观测面的空间波数域异常电场。其具体实施过程为：
将三维电阻率异常体的电流密度使用电流密度的二次插值函数表示，计算异常体第m个空间单元（即第i至第 i+2个剖分节点）在观测面
Figure 193236DEST_PATH_IMAGE134
下的空间波数域积分：
Figure 675033DEST_PATH_IMAGE134
下的空间波数域积分：
Figure 760801DEST_PATH_IMAGE135
Figure 956028DEST_PATH_IMAGE136
Figure 898576DEST_PATH_IMAGE137
式中，
Figure 539773DEST_PATH_IMAGE138
、
Figure 491548DEST_PATH_IMAGE139
、
Figure 26435DEST_PATH_IMAGE140
分别表示观测面为
Figure 405464DEST_PATH_IMAGE141
时，三维电阻率异常体内第m个空间单元积分在x，y，z方向上的空间波数域积分；
Figure 533957DEST_PATH_IMAGE142
表示空间波数域的坐标，k表示波数，
Figure 23844DEST_PATH_IMAGE143
、
Figure 429548DEST_PATH_IMAGE144
表示第i、i+1、i+2个剖分节点在x方向上的空间波数域电流密度，
Figure 182741DEST_PATH_IMAGE145
、
Figure 595268DEST_PATH_IMAGE146
、
Figure 951163DEST_PATH_IMAGE147
表示第i、i+1、i+2个剖分节点在y方向上的空间波数域电流密度，
Figure 398324DEST_PATH_IMAGE148
、
Figure 56839DEST_PATH_IMAGE149
、
Figure 222241DEST_PATH_IMAGE150
表示第i、i+1、i+2个剖分节点在z方向上的空间波数域电流密度。
Figure 365515DEST_PATH_IMAGE151
表示第m个空间单元中三个插值点，即第i、i+1、i+2个剖分节点的二次形函数，
Figure 667184DEST_PATH_IMAGE152
和
Figure 496599DEST_PATH_IMAGE153
分别表示第i和i+2个剖分节点z方向的长度范围，另外有：
Figure 211614DEST_PATH_IMAGE154
；
将
Figure 846995DEST_PATH_IMAGE155
、
Figure 206432DEST_PATH_IMAGE156
、
Figure 269066DEST_PATH_IMAGE157
累加，得到累加结果即为第m个空间单元在观测面
Figure 956531DEST_PATH_IMAGE158
下的异常电位；
将地下异常区域每个空间单元在观测面
Figure 130023DEST_PATH_IMAGE158
下的异常电位累加，即得到观测面的空间波数域异常电场：
Figure 609546DEST_PATH_IMAGE159
式中，
Figure 905398DEST_PATH_IMAGE160
为观测面的空间波数域异常电位。
本实施例中，观测面的空间波数域异常电场的获取过程具体为：
利用二维傅立叶变换的性质，基于观测面
Figure 267109DEST_PATH_IMAGE161
的空间波数域异常电位得到空间域异常电场，为：
Figure 181976DEST_PATH_IMAGE162
Figure 578322DEST_PATH_IMAGE163
Figure 769044DEST_PATH_IMAGE164
式中，
Figure 680368DEST_PATH_IMAGE165
、
Figure 461242DEST_PATH_IMAGE166
为观测面的空间波数域异常电场
Figure 649778DEST_PATH_IMAGE167
在x，y，z方向的三分量，
Figure 225116DEST_PATH_IMAGE168
表示虚数单位。
步骤5，基于观测面的空间波数域异常电位得到空间域异常电场，进一步得到观测面当前的总电场，其具体实施过程为：
对观测面的波数域异常电场
Figure 108889DEST_PATH_IMAGE169
进行二维傅里叶反变换，得到观测面的空间域异常电位
Figure 427875DEST_PATH_IMAGE170
，为
Figure 736497DEST_PATH_IMAGE171
Figure 545053DEST_PATH_IMAGE172
Figure 368652DEST_PATH_IMAGE173
基于观测面的空间域异常电场与背景电场得到观测面当前的总电场，为：
Figure 163433DEST_PATH_IMAGE174
式中，
Figure 388878DEST_PATH_IMAGE175
为观测面当前的总电场。
步骤6，判断观测面当前的总电场是否满足以下收敛条件：
Figure 883182DEST_PATH_IMAGE176
式中，
Figure 866182DEST_PATH_IMAGE177
为期望收敛时的数值精度，
Figure 526970DEST_PATH_IMAGE178
为观测面的初始电场；
若满足，则将观测面的波数域异常电位
Figure 669239DEST_PATH_IMAGE179
进行二维反傅里叶变换得到
Figure 960543DEST_PATH_IMAGE180
，将
Figure 493155DEST_PATH_IMAGE181
与背景电位
Figure 770684DEST_PATH_IMAGE182
累加，并将累加结果
Figure 908404DEST_PATH_IMAGE183
作为输出值；
否不满足，将观测面当前的总电场
Figure 167347DEST_PATH_IMAGE184
作为初始电场
Figure 249573DEST_PATH_IMAGE185
，即令
Figure 517743DEST_PATH_IMAGE186
，重复步骤3-6，直至观测面当前的总电场满足收敛条件。
下面结合具体的示例，验证本实施例所提出的方法用于计算三维电阻体的电位计算方法的效果进行验证。
目标区域存有一低阻异常体，如图3所示。异常体沿x，y方向长度为1000m，沿z方向250m，顶面埋深50m，电阻率为
Figure 244390DEST_PATH_IMAGE187
，背景电阻率为
Figure 674235DEST_PATH_IMAGE188
。令异常体中心位置在地面的投影为坐标原点，供电电极A和B坐标分别为(-1000m,0, 0)和(1000m,0,0)，电流大小为1A。模型空间单元剖分节点数为101×101×51，x，y方向剖分间隔都为10m，z方向剖分间隔为5m。采用中梯装置测量，MN=20m。期望收敛的相对残差为
Figure 758603DEST_PATH_IMAGE189
。
本实施例方法利用Fortran语言编程实现，运行程序所用的个人台式机配置为：CPU-Inter Core i5-4590，主频为3.3GHz，内存为12GB。本实施例的三维电阻率低阻体的电位方法计算达到期望收敛相对残差为10-4时，需要迭代8次，计算总时间为6.62s，由此可见本实施例具有很高的计算效率。图4为沿着x方向从-500m至500m测量得到的视电阻率和相对误差图，具有本实施例计算得到的视电阻图与使用美国犹他大学IETEM3D软件计算得到的视电阻率图，图5为图4中两种算法之间的相对误差图。从图5中可看出两种算法的相对误差小于0.35%，验证了本实施例方法的正确性，说明本实施例方法具有很高的计算精度。
以上所述仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是在本发明的发明构思下，利用本发明说明书及附图内容所作的等效结构变换，或直接/间接运用在其他相关的技术领域均包括在本发明的专利保护范围内。

Claims (7)

1.一种三维直流电阻率法数值模拟方法，其特征在于，包括如下步骤：
步骤1，根据地下结构电阻率分布情况，构建目标区域，使三维电阻率异常体完全包含在该目标区域，并将目标区域沿x、y、z方向剖分为若干空间单元；
步骤2，对目标区域中剖分节点的电阻率进行赋值，并得到各个剖分节点的背景电场与背景电位，并将背景电场作为初始电场；
步骤3，基于各个剖分节点的电阻率与初始电场得到各个剖分节点的散射电流，进一步得到各个剖分节点的空间波数域散射电流；
步骤4，基于三维电阻率异常体内各个剖分节点的空间波数域散射电流，得到观测面的空间波数域异常电位，进一步得到观测面的空间波数域异常电场；
步骤5，基于观测面的空间波数域异常电场得到空间域异常电场，进一步得到观测面当前的总电场；
步骤6，判断观测面当前的总电场是否满足收敛条件：
若是则将观测面的空间波数域异常电位进行二维傅里叶变换后与背景电位累加，并将累加结果作为输出值；
否则将观测面当前的总电场作为初始电场，并重复步骤3-6；
步骤4中，所述基于三维电阻率异常体内各个剖分节点的空间波数域散射电流，得到观测面的空间波数域异常电位，具体为：
将三维电阻率异常体的电流密度使用电流密度的二次插值函数表示，计算三维电阻率异常体内第m个空间单元在观测面
Figure 13490DEST_PATH_IMAGE001
下的空间波数域积分，其中，第m个空间单元即第i至第 i+2个剖分节点：
Figure 717004DEST_PATH_IMAGE002
Figure 403200DEST_PATH_IMAGE003
Figure 875770DEST_PATH_IMAGE004
式中，
Figure 254798DEST_PATH_IMAGE005
、
Figure 180029DEST_PATH_IMAGE006
、
Figure 669916DEST_PATH_IMAGE007
分别表示观测面为
Figure 262572DEST_PATH_IMAGE008
时，三维电阻率异常体内第m个空间单元积分在x，y，z方向上的空间波数域积分；
Figure 546922DEST_PATH_IMAGE009
表示空间波数域的坐标，k表示波数，
Figure 959449DEST_PATH_IMAGE010
，
Figure 987448DEST_PATH_IMAGE011
、
Figure 434610DEST_PATH_IMAGE012
分别表示x、y方向的波数，
Figure 155441DEST_PATH_IMAGE013
、
Figure 55264DEST_PATH_IMAGE014
、
Figure 886954DEST_PATH_IMAGE015
表示第i、i+1、i+2个剖分节点在x方向上的空间波数域电流密度，
Figure 188622DEST_PATH_IMAGE016
、
Figure 814776DEST_PATH_IMAGE017
、
Figure 201895DEST_PATH_IMAGE018
表示第i、i+1、i+2个剖分节点在y方向上的空间波数域电流密度，
Figure 571696DEST_PATH_IMAGE019
、
Figure 993450DEST_PATH_IMAGE020
、
Figure 790505DEST_PATH_IMAGE021
表示第i、i+1、i+2个剖分节点在z方向上的空间波数域电流密度；
Figure 664920DEST_PATH_IMAGE022
表示第m个空间单元中三个插值点，即第i、i+1、i+2个剖分节点的二次形函数，
Figure 850131DEST_PATH_IMAGE023
为背景电阻率，
Figure 860812DEST_PATH_IMAGE024
表示虚数单位；
Figure 94348DEST_PATH_IMAGE025
和
Figure 190480DEST_PATH_IMAGE026
分别表示第i和i+2个剖分节点z方向的长度范围，另外有：
Figure 167663DEST_PATH_IMAGE027
；
将
Figure 564009DEST_PATH_IMAGE028
、
Figure 702866DEST_PATH_IMAGE029
、
Figure 286295DEST_PATH_IMAGE030
累加，得到累加结果即为第m个空间单元在观测面
Figure 801589DEST_PATH_IMAGE001
下的异常电位；
将地下异常区域每个空间单元在观测面
Figure 52442DEST_PATH_IMAGE001
下的异常电位累加，即得到观测面的空间波数域异常电位：
Figure 627780DEST_PATH_IMAGE031
式中，
Figure 432925DEST_PATH_IMAGE032
为观测面
Figure 751911DEST_PATH_IMAGE001
的空间波数域异常电位。
2.根据权利要求1所述三维直流电阻率法数值模拟方法，其特征在于，步骤1中，所述将目标区域沿x、y、z方向剖分为若干空间单元，具体为：
建立三维直角坐标系，并将坐标原点置于目标区域中心位置；
对目标区域沿x、y、z方向剖分，将该区域分为若干空间单元，对应的x、y、z方向的剖分节点数量分布为Nx、Ny、Nz；
其中，沿x、y方向空间单元间距相同，沿z方向空间单元间距的大小可根据目标区域不同深度电流密度的变化灵活剖分，对电流密度变化大的区域加密剖分，电流密度变化小的区域稀疏剖分。
3.根据权利要求2所述三维直流电阻率法数值模拟方法，其特征在于，步骤2具体包括：
对目标区域中剖分节点的电阻率进行赋值，其中，存在三维电阻率异常体区域的剖分节点赋予电阻率
Figure 591691DEST_PATH_IMAGE033
，其他区域的剖分节点赋予背景电阻率
Figure 337930DEST_PATH_IMAGE034
；
基于地下供电的电流大小I和背景电阻率
Figure 161530DEST_PATH_IMAGE034
，计算各剖分节点在x，y，z方向的背景电场
Figure 753048DEST_PATH_IMAGE035
三分量
Figure 978493DEST_PATH_IMAGE036
、
Figure 630054DEST_PATH_IMAGE037
、
Figure 675371DEST_PATH_IMAGE038
以及背景电位
Figure 336159DEST_PATH_IMAGE039
，其中，
Figure 150531DEST_PATH_IMAGE040
表示编号为
Figure 238573DEST_PATH_IMAGE041
的剖分节点坐标，
Figure 505606DEST_PATH_IMAGE042
，
Figure 970086DEST_PATH_IMAGE043
，
Figure 904544DEST_PATH_IMAGE044
；
将背景电场作为初始电场，即：
Figure 897907DEST_PATH_IMAGE045
、
Figure 917816DEST_PATH_IMAGE046
、
Figure 920407DEST_PATH_IMAGE047
，其中，
Figure 709372DEST_PATH_IMAGE048
、
Figure 139216DEST_PATH_IMAGE049
、
Figure 380841DEST_PATH_IMAGE050
为初始电场
Figure 187123DEST_PATH_IMAGE051
在x，y，z方向的三分量。
4.根据权利要求3所述三维直流电阻率法数值模拟方法，其特征在于，步骤3具体包括：
基于各个剖分节点的电阻率与初始电场得到各个剖分节点的散射电流，为：
Figure 565015DEST_PATH_IMAGE052
式中，
Figure 431340DEST_PATH_IMAGE053
、
Figure 160261DEST_PATH_IMAGE054
、
Figure 504655DEST_PATH_IMAGE055
表示剖分节点坐标
Figure 2633DEST_PATH_IMAGE056
的散射电流在x，y，z方向的三分量；
采用二维傅里叶变换算法，计算各个剖分节点的空间波数域散射电流，为：
Figure 774279DEST_PATH_IMAGE057
式中，
Figure 990497DEST_PATH_IMAGE058
、
Figure 138582DEST_PATH_IMAGE059
、
Figure 491066DEST_PATH_IMAGE060
表示剖分节点
Figure 699193DEST_PATH_IMAGE061
的空间波数域散射电流在x，y，z方向的三分量，
Figure 137128DEST_PATH_IMAGE062
、
Figure 88903DEST_PATH_IMAGE063
分别表示x、y方向的波数，
Figure 561473DEST_PATH_IMAGE064
表示电流密度在空间波数域的坐标，
Figure 674922DEST_PATH_IMAGE065
表示虚数单位。
5.根据权利要求4所述三维直流电阻率法数值模拟方法，其特征在于，步骤4中，所述得到观测面的空间波数域异常电场，具体为：
利用二维傅立叶变换的性质，基于观测面
Figure 865732DEST_PATH_IMAGE001
的空间波数域异常电位得到空间域异常电场，为：
Figure 90040DEST_PATH_IMAGE066
式中，
Figure 682696DEST_PATH_IMAGE067
、
Figure 967046DEST_PATH_IMAGE068
为观测面的空间波数域异常电场
Figure 379573DEST_PATH_IMAGE069
在x，y，z方向的三分量。
6.根据权利要求5所述三维直流电阻率法数值模拟方法，其特征在于，步骤5具体为：
对观测面的波数域异常电场
Figure 673151DEST_PATH_IMAGE070
进行二维反傅里叶变换，得到观测面的空间域异常电场
Figure 120313DEST_PATH_IMAGE071
，为
Figure 575565DEST_PATH_IMAGE072
基于观测面的空间域异常电场与背景电场得到观测面当前的总电场，为：
Figure 475388DEST_PATH_IMAGE073
式中，
Figure 307078DEST_PATH_IMAGE074
为观测面当前的总电场。
7.根据权利要求6所述三维直流电阻率法数值模拟方法，其特征在于，步骤6中，所述收敛条件为：
Figure 608746DEST_PATH_IMAGE075
式中，
Figure 500479DEST_PATH_IMAGE076
为期望收敛时的数值精度，
Figure 622019DEST_PATH_IMAGE077
为观测面的初始电场。
CN202110596835.8A 2021-05-31 2021-05-31 一种三维直流电阻率法数值模拟方法 Expired - Fee Related CN113051779B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110596835.8A CN113051779B (zh) 2021-05-31 2021-05-31 一种三维直流电阻率法数值模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110596835.8A CN113051779B (zh) 2021-05-31 2021-05-31 一种三维直流电阻率法数值模拟方法

Publications (2)

Publication Number Publication Date
CN113051779A CN113051779A (zh) 2021-06-29
CN113051779B true CN113051779B (zh) 2021-08-10

Family

ID=76518732

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110596835.8A Expired - Fee Related CN113051779B (zh) 2021-05-31 2021-05-31 一种三维直流电阻率法数值模拟方法

Country Status (1)

Country Link
CN (1) CN113051779B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113962077B (zh) * 2021-10-20 2022-08-09 中南大学 三维各向异性强磁场数值模拟方法、装置、设备及介质
CN113779816B (zh) * 2021-11-10 2022-02-08 中南大学 一种基于微分法的三维直流电阻率法数值模拟方法
CN114065585B (zh) * 2021-11-22 2024-05-10 中南大学 一种基于库伦规范的三维电性源数值模拟方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107742015A (zh) * 2017-09-30 2018-02-27 中南大学 基于任意偶极‑偶极装置的直流激电法三维数值模拟方法
CN108287371A (zh) * 2018-01-31 2018-07-17 中南大学 直流电阻率无单元法中的背景网格自适应剖分方法
CN108710156A (zh) * 2018-03-13 2018-10-26 中南大学 一种直流电阻率无单元法模拟的支持域快速构造方法
CN111291316A (zh) * 2020-01-21 2020-06-16 山东大学 一种基于小波变换的多尺度电阻率反演方法及系统

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101334484B (zh) * 2008-07-22 2011-01-26 江苏大学 三维高分辨电阻率勘探及直接成像方法
US8756017B2 (en) * 2011-02-17 2014-06-17 Yangtze University Method for detecting formation resistivity outside of metal casing using time-domain electromagnetic pulse in well
US8700372B2 (en) * 2011-03-10 2014-04-15 Schlumberger Technology Corporation Method for 3-D gravity forward modeling and inversion in the wavenumber domain
US8688423B2 (en) * 2012-01-31 2014-04-01 Willowstick Technologies, Llc Subsurface hydrogeologic system modeling
CN108108579B (zh) * 2018-01-31 2020-04-14 中南大学 直流电阻率无单元法中耦合有限单元法的边界处理方法
CN108388707B (zh) * 2018-02-05 2021-07-13 三峡大学 一种三维非对称结构土壤模型下基于场路耦合的直流偏磁计算方法
CN108873084B (zh) * 2018-05-10 2019-10-08 中南大学 一种基于单位分解积分的直流电阻率无单元正演方法
CN110457777B (zh) * 2019-07-23 2023-06-09 南方电网科学研究院有限责任公司 土壤电阻率反演方法、系统及可读存储介质
CN111323830B (zh) * 2020-01-14 2021-06-25 东华理工大学 一种基于大地电磁和直流电阻率数据的联合反演方法
CN112327376B (zh) * 2020-10-13 2022-06-24 长江大学 一种探测射孔金属套管外地层电阻率的井中时域电磁法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107742015A (zh) * 2017-09-30 2018-02-27 中南大学 基于任意偶极‑偶极装置的直流激电法三维数值模拟方法
CN108287371A (zh) * 2018-01-31 2018-07-17 中南大学 直流电阻率无单元法中的背景网格自适应剖分方法
CN108710156A (zh) * 2018-03-13 2018-10-26 中南大学 一种直流电阻率无单元法模拟的支持域快速构造方法
CN111291316A (zh) * 2020-01-21 2020-06-16 山东大学 一种基于小波变换的多尺度电阻率反演方法及系统

Also Published As

Publication number Publication date
CN113051779A (zh) 2021-06-29

Similar Documents

Publication Publication Date Title
CN113051779B (zh) 一种三维直流电阻率法数值模拟方法
CN107742015B (zh) 基于任意偶极-偶极装置的直流激电法三维数值模拟方法
CN112949134B (zh) 基于非结构有限元方法的地-井瞬变电磁反演方法
CN110058315B (zh) 一种三维各向异性射频大地电磁自适应有限元正演方法
CN108509693B (zh) 三维频率域可控源数值模拟方法
Politi Quantifying the dynamical complexity of chaotic time series
CN111103627B (zh) 大地电磁tm极化模式对电场数据二维反演方法和装置
CN105426339A (zh) 一种基于无网格法的线源时域电磁响应数值计算方法
CN117538945B (zh) 三维大地电磁多分辨率反演方法、装置、设备及介质
CN104360404A (zh) 基于不同约束条件的大地电磁正则化反演方法
CN113553748A (zh) 一种三维大地电磁正演数值模拟方法
CN114065586A (zh) 一种三维大地电磁空间-波数域有限元数值模拟方法
CN107810432B (zh) 模型压缩
CN114065585A (zh) 一种基于库伦规范的三维电性源数值模拟方法
CN113792445B (zh) 一种基于积分方程法的三维大地电磁数值模拟方法
CN114065511A (zh) 起伏地形下大地电磁二维正演数值模拟方法、装置、设备及介质
CN115238550A (zh) 自适应非结构网格的滑坡降雨的地电场数值模拟计算方法
CN109343134A (zh) 一种矿井瞬变电磁探测数据分析解释方法及系统
CN108873084B (zh) 一种基于单位分解积分的直流电阻率无单元正演方法
CN110119586B (zh) 轴向电导率各向异性瞬变电磁三分量三维fdtd正演方法
Yan et al. Adaptive finite element modeling of direct current resistivity in 2-D generally anisotropic structures
CN111611737A (zh) 一种三维任意各向异性介质的海洋可控源电磁正演方法
Persova et al. Geometric 3-D inversion of airborne time-domain electromagnetic data with applications to kimberlite pipes prospecting in a complex medium
CN114970289B (zh) 三维大地电磁各向异性正演数值模拟方法、设备及介质
CN114047554B (zh) 大地电阻率模型建模方法、装置、计算机设备和存储介质

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20210810