一种直流电阻率小波伽辽金三维正演方法

1.本发明属于地球物理勘探领域，具体地而言为一种直流电阻率小波伽辽金三维正演方法。


背景技术：

2.随着经济建设的快速持续发展，矿产资源的短缺形势日益凸显，当前传统矿产勘探开发技术的理论研究不足和探测深度以及精度上的局限性以及目标体越隐蔽、复杂，导致油气勘探开发难度相应越来越大。因此，开展矿产资源勘探技术研究，提高勘探的精度，对于提高矿产资源勘探开发效益，具有重要的意义。
3.众所周知，直流电阻率法作为地球物理勘探中最古老的勘探方法之一，早已广泛应用于矿产资源勘探，水文地质调查，工程勘察等领域。随着勘探复杂度的加大，如何提高勘探的精度已然成为国内外研究的热点问题。正演是反演的基础，因此，研究高精度的电阻率三维正演对实现高效的三维反演是具有重要的理论意义和实际价值。目前为止，常用的直流电阻率数值模拟方法主要有积分方程法(iem),有限差分法(fdm)，有限单元法(fem)等。
4.直流电阻率法三维正演中存在计算速度慢，计算精度低等问题。


技术实现要素：

5.本发明所要解决的技术问题在于提供一种直流电阻率小波伽辽金三维正演方法，解决现有的方法存在的计算速度慢，计算精度低的问题。
6.本发明是这样实现的，
7.一种直流电阻率小波伽辽金三维正演方法，该方法包括：
8.对地电模型计算区域进行不均匀网格剖分为多个网格单元，得到网格单元编号，节点编号和坐标；
9.读取设计的地电模型参数，包括背景层参数、网格单元编号、节点编号和坐标；
10.计算设计的背景层的地下背景介质的一次场电位值；
11.对所有网格单元进行循环计算单元系数矩阵，接着总体合成所有网格单元的系数矩阵；
12.根据一次场电位值计算直流电阻率法小波伽辽金线性方程组的右端项；
13.加载本质边界条件，求解直流电阻率法小波伽辽金线性方程组，得到各个节点的二次场电位，进而得到总场电位值。
14.进一步地：
15.所述直流电阻率法小波伽辽金线性方程组为：
16.k
sus
＝b
17.其中：
[0018][0019][0020]ks
为系数矩阵，us为求解域中各个节点待求的场值，b为右端项，为求解域中各个节点待求的场值，b为右端项，均为2-term连接系数。
[0021]
进一步地：
[0022]
所述直流电阻率法小波伽辽金线性方程组通过下列方法得到：
[0023]
采用狄利克雷边界条件的二次场电位所满足的边值问题表达式：
[0024][0025]
其中us为待求解的二次场电位，u
p
为一次场电位，σ为介质的电导率，σ
p
为背景介质的电导率，δσ为异常电导率，δσ＝σ-σ
p
；
[0026]
二次场电位us和一次场电位u
p
分别采用小波尺度函数进行展开为：
[0027][0028][0029]
代入二次场电位所满足的边值问题表达式中得到：
[0030][0031]
令变为：
[0032][0033]
其中
[0034][0035]
左右两边乘上∫φ
m,p
‘
(x)dx∫φ
m,q
‘
(y)dy∫φ
m,r
‘
(z)dz得到：
[0036][0037]
其中
[0038][0039][0040][0041]
其中为小波尺度函数与其导数的乘积积分，称为2-term连接系数，将其表示如下所示：
[0042][0043][0044][0045]
根据小波尺度函数正交的性质得到：
[0046]
∫φ
m,p
‘
(x)φ
m,p
(x)dx＝δ
p,p
′
[0047]
∫φ
m,q
‘
(y)φ
m,q
(y)dy＝δ
q,q
′
[0048]
∫φ
m,r
‘
(z)φ
m,r
(z)dz＝δ
r,r
′
[0049]
代入方程得
[0050][0051]
其中
[0052][0053]
进一步化简得到：
[0054]ksus
＝b。
[0055]
进一步地：小波尺度函数为daubechies小波尺度函数。
[0056]
进一步地：采用非均匀几何网格对研究区域进行离散剖分。
[0057]
进一步地：采用均匀半空间或者水平层状介质作为背景介质来进行计算一次场电位值。
[0058]
进一步地：求解直流电阻率法小波伽辽金线性方程组采用并行直接求解器 mumps。
[0059]
进一步地：本质边界条件为狄利克雷边界条件，在边界处的电位值为：us|
γ0
＝0。
[0060]
本发明与现有技术相比，有益效果在于：
[0061]
本发明首次提出了一种基于daubechies小波函数的小波迦辽金法的直流电阻率三维数值模拟算法，该算法既可以有效回避引入场源所带来的奇异性的影响,同时采用不均匀网格对研究区域进行离散，在场源和异常体附近进行加密，远离场源和异常体区域采用稀疏网格，这样一方面有利于构建几何复杂地电模型，更好地模拟真实的地质情况,另一方面能够避免了全区域进行同样程度的加密，有效地减少了未知数个数，从而减少了计算量。因此，综合以上所述技术的应用，以实现高效、快速的直流电阻率法各向同性介质的数值模拟。
[0062]
本发明适用于几何复杂、物性分布的井中电法、海洋直流电法、激发极化等地球物理勘探方法数值模拟研究。
[0063]
本发明提出一种高效的数值模拟技术来求解直流电阻率三维正演问题，可以提高计算精度和减少时间。首先假设地下介质为各向同性介质。并且首次推导了基于二次场电位的电导率各向同性条件下的直流电阻率小波伽辽金方程，该算法可以有效避免源点奇异性的影响,提高数值解的精度。为了更好模拟复杂地质结构，本发明采用不均匀网格对研究区域进行离散剖分，可构建复杂的地质结构及物性分布，同时避免了对全区域进行同等程度的加密，减少了计算量. 针对常规迭代算法求解线性方程组收敛慢的问题，提出了采用并行直接求解器 mumps对线性方程组直接求解，有效地提高了计算效率。常规的迭代求解方法，进一步提高了数值解的精度。因此，本发明可以实现高精度、快速的电导率呈各向同性条件下的直流电阻率法数值模拟。
附图说明
[0064]
图1是本发明实施例提供的小波伽辽金法的尺度函数(a)和小波函数图(b)；
[0065]
图2是本发明实施例提供的方法流程图；
[0066]
图3是本发明实施例建立的水平层状电导率各向同性介质模型；
[0067]
图4是本发明实施例各向同性层状模型拟解析解与小波伽辽金数值解、有限元数值解的对比；
[0068]
图5是本发明实施例各向同性介质层状模型小波伽辽金数值解和有限元数值解与拟解析解的相对误差的对比；
[0069]
图6是本发明实施例建立的各向同性介质中含有异常体模型(a)各向同性介质模型的y＝0的垂直断面示意图；(b)各向同性介质模型的z＝5m的平面示意图；
[0070]
图7是本发明实施例建立的各向同性介质中含有异常体模型采用e-scan 法测量时的视电阻率平面图(a)供电点坐标为(0,0,0)m的z＝0处的视电阻率平面图；(b)供电点坐标为(0,-12,0)m的z＝0处的视电阻率平面图；(c)供电点坐标为(-12,-12,0)m的z＝0处的视电阻率平面图。
具体实施方式
[0071]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0072]
参见图2，直流电阻率小波伽辽金三维正演方法，包括如下的步骤：
[0073]
1)从拉普拉斯方程出发，采用小波伽辽金法得到电导率呈各向同性条件下的直流电阻率小波伽辽金线性方程组；
[0074]
2)对地电模型研究区域进行网格剖分，得到网格单元编号，节点编号和坐标等参数；
[0075]
3)读取设计的地电模型参数，包括背景层参数、电导率模型、网格单元编号、节点坐标等；
[0076]
4)计算背景层相应的一次场电位；
[0077]
5)对所有网格单元进行循环计算系数矩阵，接着总体合成所有单元的单元系数矩阵，同时，根据背景层计算出来的一次电位计算直流电阻率小波伽辽金线性方程组的右端项。
[0078]
6)加载本质边界条件，求解线性方程组，得到各个节点的二次场电位；
[0079]
7)将二次场电位与一次场电位相加，得到所有节点的总场电位值。
[0080]
其中，从电导率各向同性介质的直流电阻率法边值问题出发，
[0081]
为了避免源点奇异性的影响，引入二次场算法开展直流电阻率三维正演，在此直接给出考虑采用狄利克雷边界条件的二次场电位所满足的边值问题：
[0082][0083]
其中us为待求解的二次场电位，u
p
为一次场电位，可通过以均匀半空间或者全空间或者层状介质为背景介质计算得来(步骤4)(计算一次场电位)。σ为介质的电导率，σ
p
为背景介质的电导率，δσ为异常电导率，即δσ＝σ-σ
p
。
[0084]
二次场电位us和一次场电位u
p
分别采用小波尺度函数进行展开。小波尺度函数采用daubechies小波尺度函数。
[0085]
daubechies小波尺度函数的定义为:
[0086][0087]
相应的小波函数定义为
[0088][0089]
其中l为小波的支集长度，pj为滤波器系数。有关daubechies小波的尺度函数见图1，并且由小波的多分辨率性质，有
[0090][0091]
三维问题的小波基函数可通过x，y，z方向的一维紧凑支撑小波基函数的乘积而获得，如下所示：
[0092]
φ
r,p,q
(x,y,z)＝φr(x)φ
p
(y)φq(z)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0093]
其中φr(x),φ
p
(y),φq(z)为x，y，z三方向的小波尺度函数。
[0094]
本发明中，将研究区域剖分成有限多个六面体单元e，在各个单元中，电导率是固定不变的，单元中的二次场电位us和一次场电位u
p
分别采用小波的尺度函数进行展开为(对
应步骤2网格剖分)
[0095][0096][0097]
将方程(6)、(7)代入方程(1)中可得：
[0098][0099]
令
[0100]
方程(8)可变为
[0101][0102]
其中
[0103][0104]
为了进一步简化方程(9)，方程(9)左右两边乘上∫φ
m,p
‘
(x)dx∫φ
m,q
‘
(y)dy∫φ
m,r
‘
(z)dz,
[0105]
方程(10)进一步变化成
[0106][0107]
其中：
[0108][0109][0110][0111]
其中为小波尺度函数与其导数的乘积积分，一般称为 2-term连接系数。可将其表示如下所示：
[0112][0113][0114][0115]
根据小波尺度函数正交的性质
[0116]
∫φ
m,p
‘
(x)φ
m,p
(x)dx＝δ
p,p
′
ꢀꢀꢀꢀ
(18)
[0117]
∫φ
m,q
‘
(y)φ
m,q
(y)dy＝δ
q,q
′
ꢀꢀꢀ
(19)
[0118]
∫φ
m,r
‘
(z)φ
m,r
(z)dz＝δ
r,r
′
ꢀꢀꢀ
(20)
[0119]
将(18)-(20)代入方程(9)，可得
[0120][0121]
其中
[0122][0123]
则方程(21)可进一步化简为
[0124]ksus
＝b
ꢀꢀꢀꢀ
(23)
[0125]
其中(步骤5)(计算系数矩阵以及右端项)
[0126][0127][0128]
通过求解方程(23)，便可获得二次场电位，进而求出总场电位。
[0129]
为了求解方程(23)，还需要加入相应的边界条件。本发明采用狄利克雷边界条件，即在边界处的电位值为0：(步骤6)(加入本质边界条件)
[0130]us
|
γ0
＝0
ꢀꢀꢀ
(24)
[0131]
步骤6中线性方程组的求解：
[0132]
直接求解法相比于krlov子空间迭代算法，因求解精度高,而且稳定性好等优点而被广泛用于求解大型稀疏线性方程组，特别是能够一次性求解多右端项问题，能有效加快求解正演问题的速度。因此，本发明采用并行直接求解器mumps对线性方程组进行求解。
[0133]
为了验证本发明提出的基于小波伽辽金法(wgm)直流电阻率三维正演的正确性，设计一个两层水平层状介质模型。如图3所示，第一层厚度为100m，电阻率100欧姆米，第二层的电阻率为1000欧姆米。坐标原点位于地标中心，采用对称四级装置进行测量，电流大小为1a。采用数字滤波算法计算拟解析解，并将小波伽辽金法计算的数值解与解析解和分别采用有限单元法所计算的数值解进行对比。对比结果如图4所示。从图4可知，小波伽辽金法数值解与解析解、有限元解基本吻合，验证了本发明的正确性。图5为wgm和 fem两种数值模
拟算法的精度对比，可以明显看到，wgm解的精度比有限元高。表明本发明提出的直流电阻率wgm正演算法是正确的，可行的。表 1是有限单元法和wgm计算时消耗的时间对比，从表1可以看到，本发明提出的小波伽辽金算法计算速度比有限单元快。进一步表明小波伽辽金法可作为直流电阻率法数值模拟新途径。
[0134]
表1不同数值模拟方法的计算时间对比
[0135][0136]
为了进一步验证本发明方法的有效性，设计一个均匀半空间中有一个低阻异常体模型，如图6(a)和图6(b)所示。低阻异常体的长宽高均为6m，埋深为3m，中心位置为(0，0，6)。电阻率为50欧姆米,围岩电阻率为1000 欧姆米。坐标原点设置在模型的中心位置。供电电流为1a，设计三个供电点分别进行供电，具体坐标为(0,0,0)、(0,-12,0)、(-12,-12,0)。采用e-scan法进行测量，计算结果如图5所示。
[0137]
图7为采用e-scan法测量时的视电阻率平面图。图7(a)-(c)分别供电点为(0,0,0)、(0,-12,0)、(-12,-12,0)的视电阻率异常平面图。图7是本发明实施例建立的各向同性介质中含有异常体模型采用e-scan法测量时的视电阻率平面图(a)供电点坐标为(0,0,0)m的z＝0处的视电阻率平面图；(b)供电点坐标为(0,-12,0)m的z＝0处的视电阻率平面图；(c)供电点坐标为(-12,-12,0) m的z＝0处的视电阻率平面图。从图7(a)可以看到，供电点位于低阻体正上方时，视电阻率平面图呈现出一个规则对称的圆形异常,并且异常规模大小与异常体在地表上的投影位置相对应。当供电点远离异常体时，从图7(b)-(c)可以看到出现了高、低阻异常成对的现象，并且高、低阻异常中心与供电点处于同一条直线上。表2为有限单元法和小波伽辽金法cpu时间对比，从表2可以发现，在各向同性条件下，wgm消耗的时间均比fem少。
[0138]
表2不同数值模拟方法在不同的模型下的计算时间对比
[0139][0140]
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。