一种瞬变电磁小波伽辽金2.5维正演算法

1.本发明涉及地理勘探技术领域，具体涉及一种瞬变电磁小波伽辽金2.5维正演算法。


背景技术：

2.随着我国经济建设的快速持续发展，矿产资源的短缺形势日益凸显。众所周知，瞬变电磁法作为地球物理勘探中比较有效的勘探方法之一，早已广泛应用于矿产资源勘探，水文地质调查和工程勘察等领域。随着勘探复杂度的加大，如何提高勘探的精度已然成为国内外研究的热点问题。正演是反演的基础，因此，研究高精度的瞬变电磁2.5正演对实现高效的2.5维反演是具有重要的理论意义和实际价值。
3.目前，常用的直流电阻率数值模拟方法主要有积分方程法(iem)、有限差分法(fdm)和有限单元法(fem)等。然而，现有的模拟方法容易受到引入场源所带来的奇异性的影响，从而使得勘探的精度降低的情况。


技术实现要素：

4.针对现有技术存在的不足，本发明提出一种瞬变电磁小波伽辽金2.5维正演算法，以解决现有的模拟方法容易受到引入场源所带来的奇异性的影响，从而使得勘探的精度降低的问题。
5.为达到上述目的，本发明采用如下技术方案：一种瞬变电磁小波伽辽金2.5维正演算法，包括以下步骤：
6.读取模型区域的模拟文件和与模型区域中的地下介质为各向同性介质的背景层介质文件；
7.对模拟文件中的模型区域采用非均匀几何网格进行离散剖分，得到模型区域参数，所述模型区域参数包括网格单元、网格点和模型介质电导率；
8.读取背景层介质文件中的背景层介质和背景层介质电导率；
9.对所述背景层介质进行计算获得一次场电磁场值；
10.读取根据小波伽辽金法对偏微分方程进行推导得出的线性方程组；
11.通过模型介质电导率根据刚度矩阵公式得出刚度矩阵；
12.通过一次场电磁场值、模型介质电导率和背景层介质电导率根据右端项公式计算出线性方程组的右端项；
13.加载本质边界条件，根据所述刚度矩阵和线性方程组的右端项对线性方程组进行求解，得到模型区域中每个网格点的二次场电磁场值；
14.将所有网格点的二次场电磁场值与一次场电磁场值相加，得到所有网格点的总场电磁场值。
15.优选的，所述线性方程组为：
16.kfe＝beꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)；
17.kfh＝bhꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)；
18.其中，刚度矩阵公式为：
[0019][0020]
右端项公式为：
[0021][0022][0023]
其中：
[0024][0025][0026][0027][0028][0029]
其中：和分别为y方向的一次电磁场分量；ky是走向方向的波数；
[0030]
σ为模型介质电导率，σ
p
为背景层介质电导率，δσ为异常电导率；
[0031]
即δσ＝σ-σ
p
；s是拉普拉斯变量；
[0032]
为小波尺度函数与其导数的乘积积分，为2-term连接系数；k是刚度矩阵，区域网格单元e，fe表示为y方向的二次电磁场分量fe，b为右端项。
[0033]
优选的，所述偏微分方程满足拉式域和走向方向的瞬变电磁二次场，偏微分方程具体为：
[0034][0035][0036]
优选的，所述小波伽辽金法所采用的2.5维小波尺度函数可通过x，z方向的一维紧凑支撑小波基函数的乘积得到。
[0037]
优选的，所述2.5维小波尺度函数的具体推导步骤为：
[0038]
将daubechies小波尺度函数定义为:
[0039][0040]
得出相应的小波函数，定义为：
[0041][0042]
其中，l为小波的支集长度，pj为滤波器系数；
[0043]
由小波的多分辨率性质，得出：
[0044][0045]
即，通过x，z方向的一维紧凑支撑小波基函数的乘积得到的2.5维小波基函数为：
[0046][0047]
优选的，所述线性方程组的具体推导步骤为：
[0048]
二次场电电磁场(fe,fh)和一次电磁场分别采用小波的尺度函数进行展开为：
[0049][0050][0051]
其中，和分别是y方向的二次电磁场分量；
[0052]
将方程(17)和(18)代入方程(11)-(1)中可得：
[0053][0054][0055]
令
[0056]
方程(19)-(20)可变为：
[0057][0058][0059]
其中：
[0060][0061]
将方程(21)-(22)左右两边乘上和并对所有的网格点进行积分，得到线性方程组(1)-(10)。
[0062]
优选的，所述线性方程组采用并行直接求解器mumps进行求解。
[0063]
优选的，所述本质边界条件采用狄利克雷边界条件为：
[0064]
x
|γ
＝0 (24)。
[0065]
与现有技术相比，本方案产生的有益效果是：
[0066]
1、本发明采用不均匀网格对模型区域的研究区域进行离散剖分，能够构建复杂的地质结构及物性分布，避免了对全区域进行同等程度的加密，更好地模拟真实的地质情况，同时能够避免了全区域进行同样程度的加密，有效地减少了未知数个数，从而减少了计算量。
[0067]
2、与现有的迭代算法求解线性方程组相比，提出了采用并行直接求解器mumps对线性方程组直接求解，有效地提高了计算效率，同时进一步提高了数值解的精度
[0068]
3、本发明首次推导了基于走向方向的二次场电磁场的电导率各向同性条件下的瞬变电磁2.5维小波伽辽金的线性方程组，该线性方程组能够有效避免源点奇异性的影响，提高数值解的精度。
附图说明
[0069]
为了更清楚地说明本发明具体实施方式，下面将对具体实施方式中所需要使用的附图作简单地介绍。在所有附图中，各元件或部分并不一定按照实际的比例绘制。
[0070]
图1为本发明一种瞬变电磁小波伽辽金2.5维正演算法的流程框图；；
[0071]
图2为本发明一种瞬变电磁小波伽辽金2.5维正演算法中小波伽辽金法的尺度函数的波形图；
[0072]
图3为本发明一种瞬变电磁小波伽辽金2.5维正演算法中小波伽辽金法的小波函数的波形图；
[0073]
图4为试验例1中建立的水平层状电导率各向同性介质模型的示意图；
[0074]
图5为试验例1中各向同性层状模型拟解析解与小波伽辽金数值解、有限元数值解的对比图；
[0075]
图6为试验例1中各向同性层状模型小波伽辽金数值解和有限元数值解与拟解析解的相对误差的对比图；
[0076]
图7为试验例2中建立的各向同性介质中含有异常体模型；
[0077]
图8为建立的各向同性介质中含有异常体模型多时刻瞬变电磁x方向的水平电动势的剖面曲线；
[0078]
图9为建立的各向同性介质中含有异常体模型多时刻瞬变电磁y方向的水平电动势的剖面曲线；
[0079]
图10为建立的各向同性介质中含有异常体模型多时刻瞬变电磁垂直电动势的剖面曲线。
具体实施方式
[0080]
下面将结合附图对本发明技术方案的实施例进行详细的描述。以下实施例仅用于更加清楚地说明本发明的技术方案，因此只作为示例，而不能以此来限制本发明的保护范围。
[0081]
请参阅图1-图10，一种瞬变电磁小波伽辽金2.5维正演算法，包括以下步骤：
[0082]
步骤一、读取地电模型的模拟文件和与地电模型中的地下介质为各向同性介质的背景层介质文件。
[0083]
步骤二、对模拟文件中的模型区域采用非均匀几何网格进行离散剖分，得到模型区域参数，模型区域参数包括网格单元、网格点和模型介质电导率。
[0084]
步骤三、读取背景层介质文件中的背景层介质和背景层介质电导率。
[0085]
步骤四、对背景层介质进行计算获得一次场电磁场值。具体的，可采用均匀半空间或者水平层状介质作为背景介质通过现有的算法进行计算。
[0086]
读取根据小波伽辽金法对偏微分方程进行推导得出的线性方程组。
[0087]
通过模型介质电导率根据刚度矩阵公式得出刚度矩阵。
[0088]
通过一次场电磁场值、模型介质电导率和背景层介质电导率根据右端项公式计算出线性方程组的右端项。
[0089]
加载本质边界条件，根据刚度矩阵和线性方程组的右端项对线性方程组进行求解，得到模型区域中每个网格点的二次场电磁场值。
[0090]
具体的，本步骤中，具体的线性方程组为：
[0091]
kfe＝beꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)。
[0092]
kfh＝bhꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)。
[0093]
其中，刚度矩阵公式为：
[0094][0095]
右端项公式为：
[0096][0097][0098]
其中：
[0099][0100][0101][0102][0103][0104]
其中：和分别为y方向的一次电磁场分量；ky是走向方向的波数。
[0105]
σ为模型介质电导率，σ
p
为背景层介质电导率，δσ为异常电导率。
[0106]
即δσ＝σ-σ
p
；s是拉普拉斯变量。
[0107]
为小波尺度函数与其导数的乘积积分，为2-term连接系数；k是刚度矩阵，区域网格单元e，fe表示为y方向的二次电磁场分量fe，b为右端项。
[0108]
线性方程组的具体推导步骤为：
[0109]
偏微分方程满足拉式域和走向方向的瞬变电磁二次场，偏微分方程具体为：
[0110][0111][0112]
同时，小波伽辽金法所采用的2.5维小波尺度函数可通过x，z方向的一维紧凑支撑小波基函数的乘积得到。
[0113]
其中，2.5维小波尺度函数的具体推导步骤为：
[0114]
将daubechies小波尺度函数定义为:
[0115][0116]
得出相应的小波函数，定义为：
[0117][0118]
其中，l为小波的支集长度，pj为滤波器系数。
[0119]
由小波的多分辨率性质，得出：
[0120]
[0121]
即，通过x，z方向的一维紧凑支撑小波基函数的乘积得到的2.5维小波基函数为：
[0122][0123]
其中，线性方程组的具体推导步骤为：
[0124]
二次场电电磁场(fe,fh)和一次电磁场分别采用小波的尺度函数进行展开为：
[0125][0126][0127]
其中，和分别是y方向的二次电磁场分量；
[0128]
将方程(17)和(18)代入方程(11)-(1)中可得：
[0129][0130][0131]
令
[0132]
方程(19)-(20)可变为：
[0133][0134][0135]
其中：
[0136][0137]
将方程(21)-(22)左右两边乘上和并对所有的网格点进行积分，得到线性方程组(1)-(10)。
[0138]
其中，本质边界条件采用狄利克雷边界条件为：
[0139]
x
|γ
＝0 (24)。
[0140]
本实施例中，线性方程组采用并行直接求解器mumps进行求解。
[0141]
步骤五、将所有网格点的二次场电磁场值与一次场电磁场值相加，得到所有网格点的总场电磁场值。
[0142]
本发明提出一种新的技术来求解瞬变电磁2.5维正演问题，首先假设地下介质为各向同性介质。并且首次推导了基于走向方向的二次场电磁场的电导率各向同性条件下的瞬变电磁2.5维小波伽辽金方程，该算法可以有效避免源点奇异性的影响,提高数值解的精度。为了更好模拟复杂地质结构，本发明采用不均匀网格对研究区域进行离散剖分，可构建复杂的地质结构及物性分布，同时避免了对全区域进行同等程度的加密，减少了计算量.针对常规迭代算法求解线性方程组收敛慢的问题，提出了采用并行直接求解器mumps对线性方程组直接求解，有效地提高了计算效率。常规的迭代求解方法，进一步提高了数值解的精度。因此，本发明可以实现高精度、快速的电导率呈各向同性条件下的瞬变电磁2.5维数值模拟。
[0143]
验证试验：
[0144]
试验例1
[0145]
为了验证本发明提出的基于小波伽辽金法(wgm)瞬变电磁2.5维正演算法的正确性，如图3所示，设计一个地电模型为三层水平层状介质模型。第一层厚度为200m，电阻率100欧姆米，第二层的电阻率为10欧姆米，厚度为100m。第三层介质的电阻率为1000欧姆米。坐标原点位于地标中心，采用中心回线装置进行测量，电流大小为10a，线圈大小为500m*500m。沿x轴和z轴，每个网格点间隔20m，网格总数覆盖4000m的距离。
[0146]
验证方法：采用现有的数字滤波算法计算拟解析解，并将本发明中的线性方程组计算的数值解与解析解和分别采用有限单元法所计算的数值解进行对比。
[0147]
验证结论：从图4中的对比结果可知，本发明中基于小波伽辽金法的瞬变电磁2.5维正演的线性方程组计算出的数值解与解析解、有限元解基本吻合，验证了本发明中求解线性方程组的正确性。
[0148]
请一并参阅图5，为基于小波伽辽金法的瞬变电磁2.5维正演线性方程组和fem两种数值模拟算法的精度对比，可以明显看到，基于小波伽辽金法的瞬变电磁2.5维正演线性方程组解的精度稍微比有限元高。表明本发明提出的基于小波伽辽金法的瞬变电磁2.5维正演线性方程组是正确的，可行的。
[0149]
表1是有限单元法和基于小波伽辽金法的瞬变电磁2.5维正演线性方程组计算时消耗的时间对比，从表1可以看到，本发明提出的基于小波伽辽金法的瞬变电磁2.5维正演线性方程组计算速度比有限单元快。进一步表明基于小波伽辽金法的瞬变电磁2.5维正演线性方程组可作为瞬变电磁2.5维数值模拟新途径。
[0150]
methodcpu time(s)meshfem83.5200
×
200wgm64.9200
×
200
[0151]
表1不同数值模拟方法的计算时间对比
[0152]
试验例2
[0153]
为了进一步验证本文算法的有效性，如图5所示，设计一个均匀半空间中有一个低阻异常体模型。低阻异常体模型的长为120m，宽为50m。埋深为50m，中心位置为(0，0，75)。低阻异常体模型的电阻率为10欧姆米，围岩电阻率为1000欧姆米。坐标原点设置在模型的中
心位置。供电电流为1a，地表上的方环源尺寸为50m x 50m。接收回路的等效面积为43m^2；网格的参数与前一个模型中的参数相同。
[0154]
模型2的tem响应如图7所示。从图7中可以看出，在x方向的水平电动势出现了最大值和最小值，如图8所示。最小值和最大值的位置与异常体左右边缘在地表的投影一致。y方向的水平电动势，如图9所示，反映了异常体沿走向长方向延伸的信息。图10显示了垂直电动势的剖面，该剖面被认为是解释瞬变电磁数据最可靠的信息。在取样时间为1.8e-2s之前，确定垂直电动势曲线位于水平线上。这是因为当时的信号主要反映在覆盖层中。随着采样时间的延迟，垂直电动势图逐渐反映了异常体。然后将曲线从水平直线更改为具有两个最大值的曲线。此外，两个最大值的位置与水库左右边缘在地表的投影一致。随着采样时间的进一步延迟，垂直电动势曲线变为单峰曲线，然后变为水平线，这与物理模拟中获得的结果一致。表2为有限单元法和小波伽辽金法cpu时间对比.
[0155][0156]
表2不同数值模拟方法在不同的模型下的计算时间对比
[0157]
从表2可以发现，在各向同性条件下，基于小波伽辽金法的瞬变电磁2.5维正演线性方程组消耗的时间均比fem少。
[0158]
以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围，其均应涵盖在本发明的权利要求和说明书的范围当中。