基于任意偶极‑偶极装置的直流激电法三维数值模拟方法与流程

本发明具体涉及一种基于任意偶极-偶极装置的直流激电法三维数值模拟方法。

背景技术：

地球物理学是通过定量的物理方法(如：地震弹性波、重力、地磁、地电、地热和放射能等方法)研究地球以及寻找地球内部矿藏资源的一门综合性学科，研究范围包括地球的地壳、地幔、地核和大气层等。近年来，地球物理学发展越来越快，也越来越受到人们的重视。

直流激电法是地球物理学中一个非常重要的探测方法。直流激电法根据激发和和观测条件的不同可以分为地面直流激电法、井中直流激电法和隧道直流激电法，该方法作为一种重要的电法勘探方法，在矿产资源、水资源探测、考古、环境和灾害监测、隧道超前预报、路基岩溶及煤田采空区探测等领域得到了大量的应用，并取得了良好的实际效果。

直流激电法数值模拟方法与直流电阻率法完全相同，极化率数据的反演也可以在电阻率反演的基础上得到。直流激电法与直流电阻率法一样，其观测方式与数值模拟具有直接的联系，直流电法供电装置主要包括单极供电装置和偶极供电装置，偶极供电装置测量数据由于极化信息丰富和分辨率高，获得了更广泛的应用。但是在三维勘探中，由于受电场零线的影响，偶极-偶极装置很难实现全区观测。考虑到野外施工环境和三维正反演方法的限制，基于任意布极的三维测量面临极大的困难，目前主要包括如下两种观测方式：①采用以e-scan为代表的规则三维测网进行电极的布设；②将多条二维测线组成三维测网。

e-scan方法采用的是单极-单极(偶极)装置，将大量的电极连接到一根多芯电缆上，由计算机对多路传输系统进行控制，采用电极转换的方式确定供电电极和测量电极。其基本思想是将地面上的电极按照二维网格布置，每一个网格节点上有一根电极，假设共有n根电极，则完成一次全测量所需的次数为n＝n(n-1)。具体操作过程为：从m＝1(m为电极编号)开始为供电电极，其余电极转换为测量电极，接着m＝2转换为供电电极，其余电极转换为测量电极，直到m＝n时为止。这种测量方法当电极个数n过大时，对于单通道的二维测量仪进行三维测量来说工作效率较低。为了提高野外施工工作效率，jonathan(2005)针对已有的e-scan三维测量方式做了改进，采用“十字-对角方向测量”，即以源点为中心以45°为间隔进行旋转，测量旋转方向上的测点，并通过模型和实测资料进行了验证，这种测量方式与以往测测量方式几乎不降低地下目标体的信息量。

但是，e-scan测网方式存在的缺点包括：

1)由于采用单极-单极或三极装置，很难向地下供大电流，只能够进行浅层勘探，此外由于“无穷远”电极的存在，野外施工不方便；

2)所有电极都连接到同一根电缆上，造成野外施工不够灵活，布极要求非常严格，在地形较为复杂的地区很难施工；

3)测量数据可能过于密集，造成施工资源的浪费。

而二维测线组成三维测网的施工方式虽然简单可行，但由于二维测线获取地下目标体的三维信息有限，难以获得高精度成像结果。

技术实现要素：

本发明的目的在于提供一种能够实现全区观测、适应性好、数据处理简单可靠且准确的基于任意偶极-偶极装置的直流激电法三维数值模拟方法。

本发明提供的这种基于任意偶极-偶极装置的直流激电法三维数值模拟方法，包括如下步骤：

s1.采用三维直流激电任意偶极-偶极装置对待测区域进行全区观测，从而获得待测区域测点的水平电场和

s2.根据步骤s1得到的待测区域测点的水平电场和计算得到实测数据电阻率合成电场和激发极化合成电场

s3.计算步骤s2得到的实测数据的归一化背景电场并计算实测数据归一化传输函数和极化率归一化传输函数

s4.根据步骤s1中观测系统参数对理论模型进行有限元数值模拟，获得理论模型测点电阻率合成电场和激发极化合成电场

s5.根据步骤s2中选取的背景电阻率模型和步骤s4中模型网格剖分参数，计算理论模型归一化背景电场

s6.计算理论模型电阻率归一化传输函数和极化率归一化传输函数至此，分别获得了实测数据和理论模型归一化传输函数，从而为三维直流激电归一化传输函数数据的反演提供条件。

步骤s1所述的采用三维直流激电任意偶极-偶极装置对待测区域进行全区观测，具体为采用如下规则进行全区观测：

r1.发电电极和测量电极分开；

r2.发射电极为偶极装置，在待测区域内可以沿任意方向布设，且可以布设于地表或井中；

r3.测量电极为偶极装置，测量方式为矢量观测方式，测点在待测区域内可以任意布设，且在任意小区域内可以对测点进行加密或稀疏；

r4.矢量测量时，先后沿着x方向和y方向观测测点的水平电场和

r5.对供电电极进行断电操作，观测激发极化电场和

步骤s2所述的实测数据电阻率合成电场和激发极化合成电场具体为采用如下算式进行计算：

步骤s3所述的计算实测数据的归一化背景电场具体为采用如下步骤计算：

a.选取计算归一化背景场的电阻率模型；

b.采用网络加密-收缩算法对步骤a选定的电阻率模型进行网格离散，从而提升归一化背景场的计算精度和效率；

c.根据步骤b对背景场电阻率模型进行网格离散后，按照野外观测系统参数进行有限单元法数值模拟，计算得到实测数据归一化背景合成电场并计算实测数据电阻率归一化传输函数和极化率归一化传输函数具体采用如下算式进行计算：

步骤s4所述的理论模型有限单元法数值模拟，具体为采用如下步骤进行计算：

a.设计电阻率模型和极化率模型，并采用结构化网格对模型进行离散；

b.按照步骤a对模型离散后，根据野外观测系统参数采用有限单元法求解三维直流激电微分方程组，计算出模型中每个节点的电位值；

c.根据步骤b得到的每个测点电位u后，采用下式计算得到测点水平电场分量和

e＝-▽u

式中▽为梯度运算，即电场e与电位u为负梯度关系。

d.根据步骤c计算得到的水平电场分量计算理论模型电阻率合成电场和激发极化合成电场具体为采用如下算式进行计算：

步骤s5所述的计算理论模型归一化背景电场，具体为采用如下步骤进行计算：

y1.选取计算理论模型归一化背景场的电阻率模型，该电阻率模型与步骤s3中实测数据归一化背景场电阻率模型相同；

y2.采用结构化网格对步骤y1中选取的电阻率模型进行网格离散，其网格离散方式与步骤s4中模型离散方法完全相同；

y3.按照步骤y2对背景电阻率模型离散后，根据野外观测系统参数采用有限单元法进行数值模拟，计算得到理论模型归一化背景电场合成背景电场的计算方式与步骤s4中相同。

步骤s6所述的计算理论模型电阻率归一化传输函数和极化率归一化传输函数，具体为采用如下算式进行计算：

本发明提供的这种基于直流激电法的三维数值模拟方法，能实现复杂条件下三维直流激电任意偶极-偶极装置全区观测三维数值模拟，并为该条件下的三维反演成像提供准备条件；矢量测量方式能够有效地探测出地下目标体的三维信息，并为以后各向异性介质的研究提供条件；偶极-偶极装置测量数据获取地下目标体信息更为丰富，分辨率更高；不需要计算视电阻率，不受电场零线的限制，其布极更为灵活，可以针对重点勘探区域任意加密测点密度，能够适应复杂环境的勘探要求；在数值模拟中，由于采用了归一化处理，场源附近节点误差抵消，大大提升了数值模拟的计算精度。

附图说明

图1为本发明方法的方法流程图。

图2为本发明中的任意偶极-偶极观测装置示意图。

图3为本发明实施算例中场源位置示意图。

图4为本发明实施例一中的双层介质隆起构造模型示意图。

图5为本发明实施例一中的不同背景模型正演数值模拟结果示意图。

图6为本发明实施例一中的不同激发方位的rρ和rη计算结果示意图。

图7为本发明实施例一中的归一化传输函数分辨率的示意图。

图8为本发明实施例二中的起伏地形模型的yz剖面图。

图9为本发明实施例二中的起伏地形条件下归一化传输函数平面图。

图10为本发明实施例二中的物性结构不一致模型的yz剖面图。

图11为本发明实施例二中的物性结构不一致时归一化传输函数的平面图。

具体实施方式

如图1所示为本发明的方法流程图：本发明提供的这种基于任意偶极-偶极装置的直流激电法三维数值模拟方法，包括如下步骤：

s1.采用三维直流激电任意偶极-偶极装置对待测区域进行全区观测，从而获得待测区域测点的水平电场和具体为采用如下规则进行全区观测：

r1.发电电极和测量电极分开；

r2.发射电极为偶极装置，在待测区域内可以沿任意方向布设，且可以布设于地表或井中；

r3.测量电极为偶极装置，测量方式为矢量观测方式，测点在待测区域内可以任意布设，且在任意小区域内可以对测点进行加密或稀疏；

r4.矢量测量时，先后沿着x方向和y方向观测测点的水平电场和

r5.对供电电极进行断电操作，观测激发极化电场和

s2.根据步骤s1得到的待测区域测点的水平电场和计算得到实测数据电阻率合成电场和激发极化合成电场合成电场的计算公式如下：

s3.计算步骤s2得到的实测数据的归一化背景电场并计算实测数据归一化传输函数和极化率归一化传输函数具体计算步骤如下：

a.选取计算归一化背景场的电阻率模型；当背景模型为简单模型(如均匀水平半空间模型，水平层状介质模型)时，归一化背景电场可以根据解析表达式或者数值滤波方法计算得到；当背景模型包含起伏地形或已知模型约束条件存在时，采用数值计算获得合成电场数据。由于数值模拟受场源奇异性影响，但实测数据没有这个问题，因此，实测数据归一化背景场值的计算至少需要满足如下三个方面的条件：测点不受源点奇异性的影响、数值计算精度高(地表测点)以及计算效率快；

b.采用网络加密-收缩算法对步骤a选定的电阻率模型进行网格离散，从而提升归一化背景场的计算精度和效率；具体为：采用网格加密-收缩算法对背景模型进行网格离散，该算法主要分为两个步骤：

(1)在常规结构网格基础上均匀加密计算区域中水平方向网格节点密度，从而降低场源影响范围，提升背景场计算精度；

(2)在垂直方向采用网格收缩方法降低网格节点的规模数，从而提升背景场计算效率；

上述的网格加密-收缩算法可以参见论文：finiteelementnumericalsimulationof2.5ddirectcurrentmethodbasedonmeshrefinementandrecoarsement.appliedgeophysics,2016,13(2):257-266。

c.根据步骤b对背景场电阻率模型进行网格离散后，按照野外观测系统参数进行有限单元法数值模拟，计算得到实测数据归一化背景合成电场并分别计算实测数据电阻率归一化传输函数和极化率归一化传输函数具体采用如下算式进行计算：

s4.根据步骤s1中观测系统参数对理论模型进行有限元数值模拟，获得理论模型测点电阻率合成电场和激发极化合成电场具体为采用如下步骤进行计算：

a.设计电阻率模型和极化率模型，并采用结构化网格对模型进行离散；

b.按照步骤a对模型离散后，根据野外观测系统参数采用有限单元法求解三维直流激电微分方程组，计算出模型中每个节点的电位值；

c.根据步骤b得到的每个测点电位u后，采用下式计算得到测点水平电场分量和

e＝-▽u

式中▽为梯度运算，即电场e与电位u为负梯度关系。

d.根据步骤c计算得到的水平电场分量计算理论模型电阻率合成电场和激发极化合成电场具体为采用如下算式进行计算：

s5.根据步骤s2中选取的背景电阻率模型和步骤s4中模型网格剖分参数，计算理论模型归一化背景电场具体计算步骤如下：

y1.选取计算理论模型归一化背景场的电阻率模型，该电阻率模型与步骤s3中实测数据归一化背景场电阻率模型相同；

y2.采用结构化网格对步骤y1中选取的电阻率模型进行网格离散，其网格离散方式与步骤s4中模型离散方法完全相同；

y3.按照步骤y2对背景电阻率模型离散后，根据野外观测系统参数采用有限单元法进行数值模拟，计算得到理论模型归一化背景电场合成背景电场的计算方式与步骤s4中相同。

s6.计算理论模型电阻率归一化传输函数和极化率归一化传输函数至此，分别获得了实测数据和理论模型归一化传输函数，从而为三维直流激电归一化传输函数数据的反演提供条件。步骤s6中理论模型归一化传输函数计算方式如下：

从上式可以看出：电阻率归一化传输函数需要选取背景模型，而单独计算极化率归一化传输函数时对背景模型和网格剖分不做特殊要求。对比步骤s3中实测数据归一化传输函数可以看出：电阻率归一化传输函数计算中最大的区别在于背景电场的不同，实测数据背景电场是精确的，而理论模型背景电场在场源附近存在数值误差，但因为理论模型在场源附近同样存在数值误差，并且这种误差具有相关性，因此归一化处理后能消除场源奇异性的影响。

以下结合两个实施例对本发明方法进行进一步说明：

如图3所示为模型算例中标量场源布设位置图，s1-s5为5对偶极子源，其中s1﹑s2﹑s3和s5布设于地表，s1﹑s2和s3沿y方向布设，s5沿x方向布设，偶极距为80m；s4布设于井中，位于隆起构造下部，埋深400m，偶极距为60m。在2km×2km的测区范围内，矢量测量测点可以任意布设，以计算区域中101×101节点作为测点，测点间距为20m。

实施例一：水平地形下隆起构造模型的三维ip正演

构造三维双层介质隆起构造电阻率模型，第一层介质电阻率为40ω·m，极化率为0，厚度为100m；第二层介质电阻率为100ω·m，极化率为0；隆起构造距离顶界面为150m，厚度为100m，顶界面沿x方向走向200m，沿y方向走向360m,构造为两层介质：上层介质电阻率为10ω·m，极化率为20％，厚度为40m，下层介质电阻率为100ω·m，极化率为0。

如图4所示为双层介质隆起构造模型剖面图，其中，图4中(a)为yz剖面图，图中隆起构造呈梯形，顶宽360m，底宽760m，厚度为100m；图4中(b)为xz剖面图，图中隆起构造呈矩形，宽为200m，厚度为100m。

背景模型对电阻率归一化传输函数的影响：

为了研究背景模型的选取对电阻率归一化传输函数的影响，选取了如下5个模型计算归一化背景电场

背景模型一，电阻率为50ω·m的均匀半空间模型；

背景模型二，电阻率为40ω·m的均匀半空间模型；

背景模型三，算例中的双层介质电阻率模型(不含隆起构造)；

背景模型四，算例中含隆起构造的双层介质电阻率模型，隆起构造中电阻率均为100ω·m。

分别计算不同背景模型时电阻率归一化传输函数：

其中为隆起构造中含低阻异常体的电阻率模型计算结果；

如图5所示为偶极源位于地表(s1)时正演数值模拟结果，从图中可以看出：

1)本文定义的电阻率归一化传输函数在整个测区形态光滑,不受传统视电阻率计算中电场零线的影响；

2)即使背景模型与实际模型存在差异时(如中背景模型给的50ω·m的均匀半空间模型，而实际模型第一层介质为40ω·m)，源点附近的计算误差均能在归一化中被抵消；

3)在中，由于背景模型电阻率与真实模型第一层介质电阻率相同，且偶极子源远离异常体，因此源点附近区域测点不受地下高阻层的影响，归一化传输函数数值为1，距离源点较远的测点，受地下高阻层介质和隆起构造的影响，归一化传输函数数值逐渐增大。对比可以发现，两个模型计算结果整体形态相似，但由于背景模型较真实模型第一层介质电阻率更高，在源点附近归一化传输函数数值低于1；

背景模型越接近真实模型，电阻率归一化传输函数中异常体特征越明显，且异常差异更大(如和)，从而为具有一定地质信息的测区采用约束反演提供条件。

场源位置对归一化传输函数的影响：

分别对观测系统(图3)中4对偶极子源(s1-s4)进行正演数值模拟，并计算电阻率归一化传输函数rρ和极化率归一化传输函数rη，其中计算rρ时背景模型为双层介质模型(上述算例中模型三)。令为隆起构造中含低阻异常体和极化率的等效电阻率模型计算结果，可得计算表达式:

其中k＝1,2,3,4分别对应偶极源s1-s4。

计算结果如图6所示：

1)场源位置发生改变时，由于激发方位和照射角度的不同，归一化传输函数具有形态差异很大﹑极化信息丰富和分辨率高的特点；

2)电阻率异常和极化率异常同源时，rρ和rη具有相似的形态分布图；

极化率归一化传输函数直接反映了地下极化异常体的形态和在地表的投影位置，因此，通过对大量理论模型的数值结果进行研究,总结出规律性认识后，可以根据归一化传输函数形态分布图大致推断出地下目标体和极化体的构造特征。

归一化传输函数分辨率：

如图3所示，隆起构造在x和y方向长度不等，沿x走向为200m，沿y走向为360m。因此，本文分别对观测系统(图3)中偶极源s1和s5进行正演数值模拟，通过对比两次正演计算结果分析归一化传输函数的水平分辨率。算例中，s1沿y方向布设，s5沿x方向布设,两对偶极子源在空间上的关系为s5为s1沿逆时针旋转90°，偶极距均为80m，供电电流强度均为20a。

计算结果如图7所示，图中白色方框对应隆起构造中低阻异常体和高极化率体在地表的投影位置。电阻率归一化传输函数计算中背景模型为双层隆起构造模型(背景模型四)，从而使得在计算rρ和rη时电阻率异常体和极化率异常体同源.从图中可以看出：当偶极子源沿y方向布设时(s1)，归一化传输函数在y方向分辨率较高，y方向主异常形态与异常体在空间位置上一一对应；当偶极子源沿x方向布设时(s2)，归一化传输函数在x方向分辨率较高，x方向异常形态与异常体在空间位置上一一对应。综上可以推断出归一化传输函数沿偶极子源走向方向具有很高的分辨率。在实际勘探过程中,综合地下目标体的构造走向，采用偶极子源进行多方位和多角度照射，可以大大提升勘探效果。

实施例二：复杂起伏地形条件下三维ip正演

设计如图8所示(模型yz剖面图)的三维凹嵌地形构造电阻率模型，模型电阻率为100ω·m，极化率为1％；凹嵌地形距离水平地表30m，地形底端为正方向，边长为800m；在地形下方置入一电阻率为20ω·m﹑极化率为20％的低阻异常体，异常体距离顶界面70m，厚度为40m，沿x和y方向长度为200m。

地形对归一化传输函数的影响：

以背景电阻率为100ω·m﹑极化率为0％的凹嵌地形模型为背景场计算模型，计算归一化电阻率传输函数rρ和极化率归一化传输函数rη，其中偶极源s1位于地表(如图3所示)。计算结果如图9所示，对比水平地形算例计算结果可以发现：起伏地形对rρ和rη没有任何影响，从而可以得出结论：归一化传输函数既可以解决源点奇异性问题又能抵消地形对数值计算结果的影响，以此为基础的三维正反演方法可以适用于复杂环境勘探中。

物性结构不一致：

对图8所示的极化率模型做如下改造：取模型低阻异常体中极化率等于背景模型极化率(1％)，沿y轴对称的方向置入一个与低阻异常体大小﹑埋深相同的高极化体，极化率为20％，电阻率为100ω·m。模型yz剖面图如图10所示。

采用与上一个算例中相同的计算方法：以背景电阻率为100ω·m﹑极化率为0％的凹嵌地形模型为背景场计算模型，计算归一化电阻率传输函数rρ和极化率归一化传输函数rη，其中偶极子源s5位于地表(如图3所示)。计算结果如图11所示：当物性结构不一致时，rρ和rη平面图均能反映出电阻率异常体和极化率异常体在地表的投影位置；对比图9可以发现：当异常体为低阻高极化时，极化率归一化传输函数rη具有更大的异常数值，这与直流电阻率法对低阻体更敏感的特性相符。