基于虚拟波域法的瞬变电磁建模和反演方法与流程

本发明涉及瞬变电磁法的建模和反演技术，具体为基于虚拟波域法的瞬变电磁建模和反演方法。

背景技术：

目前，瞬变电磁法在固体矿物勘探、油气勘探、地下水探测、地质测绘和精确超前检测含水层结构等方面被广泛应用。特别是，随着高精度、智能化、灵活可靠的野外勘探仪器的出现，大大提高了其工作效率和观测精度。最近较为流行的瞬变电磁建模和反演方法是，将满足扩散方程的瞬变电磁场转化为满足波动方程的虚拟波场，然后借助地震学中已经发展起来的和正在发展的一些较为成熟的方法技术，来求解物性和几何参数。

地球物理学家在对层状介质中的大地电磁问题建模的研究中，发现大地电磁场及其满足的扩散方程，与地震波场及其满足的波动方程之间存在对应关系；随后，前苏联学者对这一对应原理进行了数学关系推导。陈本池等人在题目为“瞬变电磁波场的波场变换研究”和“瞬变电磁场的波场变换算法”的论文中采用了一种奇异值分解算法对方程系数矩阵进行分解。该方法的难点在于采样时间窗口宽度和阻尼系数的选择，这一过程需要通过大量的试算来完成，有较大的不确定性。最后，通过有限差分逆时深度偏移实现反演计算。

由于波场变换式是典型的欠定方程组问题，具有严重的不适定性。李貅等人在题目“从瞬变电磁场到波场的优化算法”的论文中，在正演计算中，采用两步最优算法来减少积分系数和离散采样点个数，应用特殊函数简化积分方程。在反演计算中，采用对瞬变电磁场时间信号分段的处理技术，通过偏差原理和牛顿迭代选择正则化参数，并运用正则化方法实现虚拟波场值得求解。

但是，奇异值分解算法的计算过程是一个冗长繁复的过程，特别是在系数矩阵比较大的时候，计算成本非常高，并且超大型矩阵的奇异值求解并不是那么容易。在有限差分法中对于边界问题的处理是一个较为棘手的问题。虽然对瞬变电磁时间信号分段的处理方式可以有效降低系数矩阵的病态程度，但是同样也引入了段与段衔接问题。

技术实现要素：

针对上述存在的问题，本发明的目的是提供一种高效的基于虚拟波域法的瞬变电磁建模和反演方法。

实现本发明目的的技术方案是：首先基于扩散场的准静态麦克斯韦方程，采用拉普拉斯变化和矢量恒等式，将瞬变电磁场从时域转换到虚拟波域，为改善线性方程组系数矩阵的谱性质，降低问题的病态程度，采用对称超松弛预处理技术和正则化共轭梯度算法求解线性方程组，得到虚拟波场值。

本发明基于虚拟波域法的瞬变电磁建模和反演方法，包括如下步骤：

(1)对于典型的地电导率，传导电流比位移电流大几个数量级，忽略位移电流，扩散场的准静态麦克斯韦方程如下：

其中，e为电场强度，h为磁场强度，μ为磁导率，σ为电导率，为在t方向的导数算子，外部源为电流密度j和磁流密度k；

(2)对以上准静态麦克斯韦方程两边求旋度有：

其中，

(3)引入虚拟波域的电场强度e^u和磁场强度h^u：

其中，和为虚拟波场的源项，为在虚拟时间t'方向的二阶导数算子，这里e^u和h^u表现为一个速度为(μσ)^-1/2，单位为的传播波；

(4)替代本构关系b＝μ0h，使用拉普拉斯变换法，实现瞬变电磁场从真实扩散域到虚拟波域的映射，如下：

其中，和分别为虚拟波域和实际扩散域中的感应电动势；

(5)将模型写成数值积分形式：

其中，f(t)为瞬变电磁场的电场分量、磁场分量或感应电动势，hj为积分系数，a(t,t'j)为核函数；

(6)方程的矩阵形式为au＝f；要得到虚拟波场值，即求解线性方程组：

ax＝b

其中，a＝[a·hj]，矩阵a包含核函数和积分系数hj，待求量x为虚拟波场u的值，f为接收到的瞬变电磁场时间信号，将f写成右端项常数b＝(b1,b2,lbn)；

(7)步骤(6)所列公式为典型的不适定问题，离散化后所得到的线性方程组是严重病态的；为共轭梯度迭代求解，首先在线性方程组两边同时左乘a^t：

a^tax＝a^tb

只要满足a是列满秩矩阵，a^ta就是对称正定矩阵；

(8)为改善系数矩阵的谱性质，降低方程系数矩阵的条件数，采用对称超松弛预处理技术：

构造预处理矩阵

m＝[(k+ωcl)^-1k^-1(k+ωcu)]/[ω(2-ω)]

其中，k、cl和cu分别为矩阵a^ta的对角阵、下三角阵和上三角阵，超松弛因子ω为(0，2)内的参数；

(9)将线性方程组转化为

m^-1a^tax＝m^-1a^tb

其中，m^-1a^ta接近单位阵，可显著改善求解的收敛速度；

(10)对于严重病态的线性方程组，还需进一步做正则化处理；上述问题的正则化方程为:

m^-1(a^ta+vi)x＝m^-1(a^tb+vx)

其中，v≥0，为正则化参数，i∈r^n×n，为单位阵。

(11)通过共轭梯度迭代求解虚拟波场的值x，共轭梯度法的基本步骤为：

给定一组初值x⁽⁰⁾∈rⁿ，计算残差r⁽⁰⁾＝a^tb-a^tax，取搜索方向d⁽⁰⁾＝r⁽⁰⁾；对k＝0,1,l，计算

x^(k+1)＝x^(k)+αkd^(k)

r^(k+1)＝r^(k)+αka^tad^(k)

d^(k+1)＝r^(k+1)+βkd^(k)

当r^(k)＝0或(d^(k),a^tad^(k))＝0，则停止计算，x^(k)＝x^*；其迭代通式如下：

m^-1(a^ta+vi)x^(k+1)＝m^-1(a^tb+vx^(k))。

上述瞬变电磁建模和反演方法的优选步骤如下：

步骤ⅰ：构造方程：m^-1(a^ta+vi)x^(k+1)＝m^-1(a^tb+vx^(k))；

步骤ⅱ：输入瞬变电磁场时间信号f；给定初值x⁽⁰⁾，超松弛因子ω和正则化参数v；

步骤ⅲ：对系数矩阵做超松弛预条件处理；

步骤ⅳ：正则化共轭梯度迭代求解线性方程组，得到全时域的虚拟波场值u。

所述步骤ⅱ中初值x⁽⁰⁾为单位向量，超松弛因子ω＝0.3，正则化参数v＝3.0e-08。

所述正则化是用一族相邻近的适定问题的解去逼近原问题的解。

本发明提高了基于虚拟波域法的瞬变电磁建模和反演方法的效率，从而改进了瞬变电磁法的性能；同时，该方法不再需要对瞬变电磁时间信号进行分段处理，避免了段与段的衔接问题，使其简单、易行，便于推广。

附图说明

图1为本基于虚拟波域法的瞬变电磁建模和反演方法实施例中h-型二维三层地电模型；

图2为本基于虚拟波域法的瞬变电磁建模和反演方法实施例中磁场强度的波场正变换结果；

图3为本基于虚拟波域法的瞬变电磁建模和反演方法实施例中磁场强度的波场反变换结果；

图4为本基于虚拟波域法的瞬变电磁建模和反演方法实施例中感应电动势的波场正变换结果；

图5为本基于虚拟波域法的瞬变电磁建模和反演方法实施例中感应电动势的波场反变换结果。

具体实施方式

实施例为二维三层地电模型的正反演模型，如图1所示，其中三层地的电导率关系为σ1>σ2<σ3，称为h-型地电模型。本实施例将h-型地电模型和均匀地电模型正演值，作为反演的输入信号：瞬变电磁场时间信号f，即方程的右端项常数b。让t在[10^-5,10^-2]范围内等间隔采样101个点，将它们代入正演模型

可直接求得瞬变电磁场时间信号f；反演的具体实施步骤如下：

步骤ⅰ：构造方程：m^-1(a^ta+vi)x^(k+1)＝m^-1(a^tb+vx^(k))；

步骤ⅱ：输入瞬变电磁场时间信号f；给定初值x⁽⁰⁾为单位向量，超松弛因子ω＝0.3和正则化参数v＝3.0e-08；

步骤ⅲ：对系数矩阵做超松弛预条件处理；

步骤ⅳ：正则化共轭梯度迭代求解线性方程组，得到虚拟波场值u。

图2为磁场在h型地电截面的衰减曲线与磁场在均匀半空间的衰减曲线对比。我们可以看到，由于二次磁场通过两个不同导电层之间的界面，h型地电截面的衰减曲线分别有一个向上和一个向下的微小波动。图3为将图2的瞬变电磁信号作为输入得到的对应虚拟波场反演值与理论值对比。在图3中，这种能够反映两个不同导电层之间的界面的波动现象更为明显。另外，虚拟波场反演值与理论值很好的吻合，表明我们的方法是有效的。在该计算中，总共需要迭代35567个时间步长，计算时间为5.235s，平均绝对误差为0.0568。

图4为感应电动势在h型地电截面的衰减曲线与感应电动势在均匀半空间的衰减曲线对比。负号表示感应电动势的方向与磁场的方向相反。图5为将图4的瞬变电磁信号作为输入得到的对应虚拟波场反演值与理论值对比。在该计算中，总共需要迭代24136个时间步长，计算时间为3.0623s，平均绝对误差为0.0741。

上述实施例，仅为对本发明的目的、技术方案和有益效果进一步详细说明的具体个例，本发明并非限于此。凡在本发明公开的范围之内所做的任何修改、同等替换、改进等，均包含在本发明的保护范围之内。