CN112949134B - 基于非结构有限元方法的地-井瞬变电磁反演方法 - Google Patents

Abstract

本发明公开了一种基于非结构有限元方法的地‑井瞬变电磁反演方法，包括以下步骤：采用非结构有限元方法为探测地质结构构建初始地‑井模型，同时根据观测电磁数据、正演得到的预测电磁数据以及对地‑井模型的反演参考模型构建正则化反演的目标函数；对初始地‑井模型进行正演计算获得预测电磁数据，依据该预测电磁数据和正则化反演的目标函数实现地‑井瞬变电磁的迭代反演，不断更新反演参考模型，直到满足反演终止条件，获取符合探测地质结构的最终反演参考模型。提升了地‑井瞬变电磁反演的速度和准确度。

Description

基于非结构有限元方法的地-井瞬变电磁反演方法

技术领域

本发明属于地质探测领域，具体涉及一种基于非结构有限元方法的地-井瞬变电磁反演方法。

背景技术

矿产资源为我国经济的快速发展提供了很大的帮助，但是随着大量的开采，地下浅部的矿产资源越来越少，已经基本开发殆尽。为了满足进一步发展的需要，必须将深部找矿作为当前能源开发的重点研究方向。

瞬变电磁(transient electromagnetic,TEM)作为一种常用的地球物理探测技术，目前已经广泛应用于矿产资源和能源勘探、地下水资源调查、地下构造特征研究等领域。TEM的基本原理是法拉第电磁感应定律。通过接地线源或不接地回线源向地下激发一次脉冲磁场，利用接地电极或不接地回线测量由导电大地产生的二次电磁场，并通过反演解释技术实现对地下目标体的刻画。

传统TEM的发射和接收均部署在地面上，由于离异常体较远导致分辨率受限，对深部异常体的探测能力较低。相比之下，地-井瞬变电磁(surface-to-borehole transientelectromagnetic,SBTEM)作为一项新技术，将接收机放入钻孔中，通过减少接收机与异常体的距离增强探测灵敏度，异常信号的强度和反演成像分辨率均可得到极大提高。

迄今为止，地-井瞬变电磁三维正演的研究已经相对成熟，理论体系也较为完善。目前在地-井瞬变电磁中应用的主要方法是有限差分法(杨海燕，徐正玉，2015；李大俊，翁爱华，2017)。也有少数学者将频域有限元法(毛立峰，王绪本,2006；李建慧,2011)应用到地-井瞬变电磁中。

有限差分法原理直观、容易在计算机上实现、实用性强，可以快速实现仿真模拟；但对物性参数复杂、几何特征不规则的模型适应性差表，严重影响仿真精度，且网格剖分的尺寸会受到很多条件的限制，很难精确模拟。相比之下，有限元法对复杂结构的建模更为灵活，空间离散也更加自由，对目标介质组成特点和复杂几何结构的模拟更为准确，能够更加有效地解决有限差分法对复杂研究对象建模的局限性问题，且有限元法比有限差分法稳定性更好。但不论是应用何种方法，都局限于正演层面，反演的相关研究十分匮乏。

发明内容

鉴于上述，本发明的目的是提供一种基于非结构有限元方法的地-井瞬变电磁反演方法，提升了地-井瞬变电磁反演的速度和准确度。

为实现上述发明目的，本发明提供以下技术方案：

一种于非结构有限元方法的地-井瞬变电磁反演方法，包括以下步骤：

1).构建初始地-井模型；

2).筛选观测地-井瞬变电磁数据中的可用数据用于反演；

3).将步骤2)的观测数据和对步骤1)中地-井模型(或迭代模型)进行非结构有限元正演得到的预测数据做数据拟合，并结合模型粗糙度构建地-井瞬变电磁三维正则化反演的目标函数；

4).对正演方程两侧求导，其中预测数据对模型参数的偏导数为灵敏度矩阵(雅可比矩阵)；对伴随正演求解方程得到灵敏度矩阵的转置与向量的乘积；

5).计算模型粗糙度，将灵敏度矩阵的转置与向量的乘积结合预测数据与观测数据的拟合差，可计算目标函数的梯度；

6).采用拟牛顿法或高斯-牛顿方法计算模型更新量，得到新的迭代模型；

7).重复执行步骤3)～步骤6)，直到达到最大迭代次数或者符合终止条件，得到最终的反演参考模型。

与现有技术相比，本发明具有的有益效果至少包括：

本发明实施例提供的基于非结构有限元方法的地-井瞬变电磁反演方法中，通过为探测地质结构构建初始地-井模型，并使用四面体网格进行空间离散，能够精确拟合复杂起伏地形和地下异常体；对初始地-井模型进行非结构矢量有限元正演计算获得预测电磁数据，可以有效避免散度条件引起的误差；依据预测电磁数据和观测数据构建正则化反演的目标函数；采用伴随正演技术对灵敏度信息进行计算和存储，采用加权平方和方法实现模型粗糙度的近似离散，具有节省计算空间提高计算效率的优点；同时基于拟牛顿方法或高斯-牛顿法实现了地-井瞬变电磁的迭代反演。获取符合探测地质结构的最终方案参考模型，依次大大提升了地-井瞬变电磁反演的效率和准确度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图做简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动前提下，还可以根据这些附图获得其他附图。

图1是本发明实施例提供的基于非结构有限元方法的地-井瞬变电磁反演方法的流程图；

图2是本发明实施例提供的均匀半空间模型示意图；

图3是本发明实施例提供的利用半空间半解析解验证本文算法的精度结果图，其中，(a)为Ex响应；(b)为Ex相对误差；(c)为Ez响应；(d)为Ez相对误差；(e)为dBy/dt响应；(f)为dBy/dt相对误差；(g)为dBz/dt响应；(h)为dBz/dt相对误差；(a)-(d)为测点1的电磁响应，而(e)-(h)为测点2的电磁响应，图中虚线代表负响应；

图4是本发明实施例提供的磁性源发射地-井瞬变电磁三维模型示意图；

图5是本发明实施例提供的磁性源地-井瞬变电磁三维反演真实模型与反演模型切片对比图；

图6是本发明实施例提供的磁性源地-井瞬变电磁三维反演参数随迭代变化曲线图。

具体实施方式

为使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例对本发明进行进一步的详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不限定本发明的保护范围。

图1是本发明实施例提供的基于非结构有限元方法的地-井瞬变电磁反演方法的流程图。如图1所示，实施例提供的基于非结构有限元方法的地-井瞬变电磁反演方法，包括以下步骤：

步骤1，构建初始地-井模型。

在构建初始地-井模型时，为探测地质结构构建几何模型，并将几何模型剖分成非结构的四面体网格单元，同时设置几何模型信息、传递数据信息以及控制参数，得到初始地-井模型。

实施例中，可以采用Delaunay方法将计算区域剖分为非结构的四面体单元。设置几何模型信息包括电导率信息、钻孔位置等，传递数据信息包括地面发射源以及接收井的个数与位置、井中接收点间隔和深度等。控制参数包括最大迭代次数和正则化因子λ等。

图2是实施例提供的均匀半空间模型示意图。构建的均匀半空间模型作为初始地-井模型，设置空气电导率为1×10⁸Ω·m，两个测试点分别在坐标(500m,0m,Z)和(-200m,500m,Z)处，Z表示深度，半空间电导率为1Ω·m/10Ω·m/100Ω·m。

实施例中，采用非结构四面体网格建模，在建模时具有较强的灵活性，能够解决无法精确拟合复杂起伏地形和地下异常体的问题。

步骤2，筛选观测地-井瞬变电磁数据中的可用数据用于反演。

实施例中，将灵敏度过大的观测数据的对应误差设置极大，可再反演中使该数据不发挥作用；

若进行反演的观测数据是实测数据，还需要叠加，去噪等预处理过程。

步骤3，构建正则化反演的目标函数。

实施例中，构建基于L2范数的目标函数，该目标函数用于正则化反演过程来优化反演参数模型，具体根据观测电磁数据、正演得到的预测电磁数据以及对地-井模型的反演参考模型构建正则化反演的目标函数为：

其中，m₀和m分别是M维的反演参考模型和解向量。d_pre和d_obs是N_d×n_channel维的预测电磁数据和观测电磁数据，N_d为接收点的个数，n_channel为时间道个数，λ为正则化因子，W_d为数据方差矩阵，W_m模型粗糙度算子；

其中在对初始地-井模型进行正演计算时，时间域电磁场满足如下Maxwell方程：

其中，E和H分别是电场和磁场，μ是磁导率，ε是介电常数，σ是电导率。j_s(r,t)是源电流密度。由(2)和(3)式，可以推导出电场满足的全波方程为：

由于初始地-井模型被剖分成非结构的四面体单元，则将模型的电场全波方程表示成有限元形式，可以有效避免散度条件引起的误差，具体为：

其中，t表示时间，矢量

表示各四面体网格单元棱边的待求解电场，介电常数矩阵

电导率矩阵

磁导率矩阵元素

以及源矢量分别为：

其中，上标e表示四面体单元，i、j分别表示四面体单元索引，ε^e表示介电常数，σ^e表示电导率，

分别表示四面体单元内电场矢量在棱边i和棱边j的插值基函数，V表示体积，

表示哈密顿算子，j_s(r,t)表示源电流密度；

对于瞬变电磁探测，当忽略位移电流时，有限元形式的电场全波方程简化并将整体单元进行组装得到：

为求解方程式(10)，需要对其中时间偏导数项进行离散。本实施例中采用阶跃波作为发射波形，由于空气和地下介质存在巨大的电性差异，如果选用显式时间离散，需要有足够短的时间步长Δt以满足稳定性条件，这对于晚期道响应计算非常不利。因此，本实施例中采用隐式后推欧拉方法进行时间离散，该方法无论时间步长Δt如何选择均无条件稳定。即将一阶后推欧拉方法应用于方程式(10)，可得：

Duⁿ⁺¹＝Buⁿ-Δtsⁿ⁺¹ (11)

其中，上标n和n+1表示递推过程，Δt表示时间步长，D＝(B+ΔtC)；

在正演过程中，还需要设置边界条件和初始条件。实施例中，采用非齐次狄利克雷边界条件(Dirichlet boundary condition)，将矩阵方程(11)在外边界处的电场设置为零，另外，对于相对比较复杂的下阶跃发射波形，将直流电场可利用节点有限元求解如下三维泊松问题得到，即

其中，j_s(r)表示源电流密度，

表示电势；

在求解得到方程式(12)中的电势后，利用导数关系求解出直流电场。

将边界条件和初始电场带入方程(11)并求解，即可获得预测电磁数据。

采用上述正演过程，通过合理加密测点、源和异常体附近的网格，稀疏扩边区域的网格，将网格数量控制在一定范围内，在保证计算精度的同时提高了计算效率。

为了检验正演过程的精度，设计如图2所示的半空间模型。利用提出的三维正演算法进行计算，并将计算结果与一维半解析解进行对比。模型参数如下：偶极子源沿x方向，长度为100m，中心坐标为(-200m,0m,0m)，电流强度为1A，空气电阻率取1×10⁸Ω·m，半空间电阻率为100Ω·m。

分别选取位于(500m,0m,500m)测点处的Ex、Ez，以及位于(-200m,500m,500m)处的dBy/dt和dBz/dt结果验证本文算法的准确性，验证结果如图2所示，其中，(a)为Ex响应；(b)为Ex相对误差；(c)为Ez响应；(d)为Ez相对误差；(e)为dBy/dt响应；(f)为dBy/dt相对误差；(g)为dBz/dt响应；(h)为dBz/dt相对误差；(a)-(d)为测点1的电磁响应，而(e)-(h)为测点2的电磁响应，图中虚线代表负响应；由图2可见，二者的相对误差基本小于5％，除|dBy/dt|零值带以外。

步骤4，伴随正演-求解灵敏度矩阵的转置与向量的乘积；

在对正则化反演的目标函数进行求解时，需要计算观测电磁数据的灵敏度和模型的粗糙度，

具体地，计算观测电磁数据的灵敏度时，将正演方程两端对第k^th个模型参数求导，可得：

其中，m_k表示模型参数，s_TEM表示源项，K表示系数矩阵，在地-井瞬变电磁正演中，预测电磁数据d从电场解e插值得到，即d＝Le,L＝L_tL_d，L_d将电场解插值到指定接收点位置的插值函数，L_t是将计算时间道的解插值到观测时间道上；

观测电磁数据的灵敏度写为：

其中，j^k是整体灵敏度矩阵中第k^th列元素，设定一个N×M维的矩阵G，N＝N_edge×n,n计算时间道的个数，即：

G＝{g₁,g₂,g₃,…,g_M} (15)

长度为N的向量g_k的表达式为：

由于正演过程中采用的是总场的方法，故有

但对于下阶跃波的电性源，初始场与模型介质的电导率相关，需要计算De_DC和Ce_DC关于模型参数的导数，从公式(14)和(15)可得：

J^T＝G^TK^-TL^T (17)

从目标函数梯度的表达式可知，需要计算灵敏度矩阵的转置J^T和一个向量的乘积，设相乘的向量为u，并定义：

w＝K^-TL^Tu (18)

则得到如下的伴随正演方程：

K^Tw＝L^Tu (19)

求解方程(19)后可得灵敏度矩阵的转置J^T和向量u的乘积：

J^Tu＝G^Tw. (20)

实施例中，采用伴随正演技术对灵敏度信息进行计算和存储，节省了计算空间提高了计算效率。

步骤5，计算目标函数的梯度。

将方程式(1)两侧对模型m^n-1求导，获得第n^th迭代中目标函数的梯度gⁿ：

其中，J是雅克比矩阵，而

求解(21)式需要计算观测电磁数据的灵敏度和模型的粗糙度，步骤4计算了灵敏度矩阵的转置与向量的乘积，而常用的模型粗糙度算子定义为模型梯度的L₂范数，对于非结构网格往往难以直接进行近似离散。本实施例中，在计算模型粗糙度时，采用加权平方和方法对非结构的四面体单元网格进行近似离散，表示式为：

式中，

其中，Vi是第i个四面体的体积，N(i)是与四面体i共顶点的所有四面体网格个数。

如此给出灵敏度矩阵的转置与向量的乘积，以及模型的粗糙度，可以计算目标函数的梯度。

步骤6，计算模型更新量，得到新的迭代模型。

依据预测电磁数据和正则化反演的目标函数实现地-井瞬变电磁的迭代反演，采用拟牛顿或高斯-牛顿方法，不断更新反演参考模型，重复步骤3～步骤6直到满足反演终止条件或达到最大迭代次数，获取符合探测地质结构的最终反演参考模型。

实施例中给定的终止条件是RMS<0.98，其中RMS的计算公式为

其中，Ndata是观测数据个数，d_pre和d_obs是预测电磁数据和观测电磁数据，W_d为数据方差矩阵。

实施例中，为了使反演结果更加接近真实情况，为反演过程中施加电导率上下限约束能够减少虚假异常。具体地，本发明基于对数参数的电导率转换函数将模型参数变化到转换域内(称为对数变换)，对反演结果进行电导率上下限约束，以减少了反演的多解性。

实施例中，依据该预测电磁数据和正则化反演的目标函数实现地-井瞬变电磁的迭代反演时，通过对正则化反演的目标函数进行求解以获得预测电磁数据和观测电磁数据的拟合差，将拟合差小于设置阈值作为反演迭代约束条件，在不满足约束条件时，依据求解正则化反演的目标函数获得的模型解向量更新方案参考模型。

实施例还提供了一种反演算例，以验证反演算法的正确性。图4是本发明实施例提供的磁性源发射地-井瞬变电磁三维模型示意图。如图4所示，起伏地形下埋藏一低阻月牙形异常体。发射线圈沿地形敷设，大小为500m×500m，发射电流1A，发射波形为阶跃波。9个接收井位于发射导线内，间距100m米等距分布。井深300m。低阻异常体位于-200m至-250m深度，长度约为400m，宽度约为200m，电阻率为5Ω·m。测点位于井中200m至300m深度，点距20m。测量时间道为1.2e-2ms至10.356ms对数等间隔分布，共29道数据。此算例使用dBz/dt作为反演数据，并加入5％高斯随机噪声模拟真实情况。空气电阻率取1×10⁸Ω·m，地下背景电阻率为100Ω·m。

图5是本发明实施例提供的磁性源地-井瞬变电磁三维反演真实模型与反演模型切片对比图，从图5可以看出三维反演结果与真实模型吻合较好且起伏地形并没有对反演结果有明显的影响，说明本发明的地-井瞬变电磁三维反演算法是可靠的，并且可以很好地解决起伏地形问题，具有一定的实用性。图6为磁性源地-井瞬变电磁三维反演参数随迭代变化曲线图。从图6可以看出，反演过程是稳定收敛的，进一步证明反演算法的正确性。

以上所述的具体实施方式对本发明的技术方案和有益效果进行了详细说明，应理解的是以上所述仅为本发明的最优选实施例，并不用于限制本发明，凡在本发明的原则范围内所做的任何修改、补充和等同替换等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于非结构有限元方法的地-井瞬变电磁反演方法，包括以下步骤：

1).构建初始地-井模型；

2).筛选观测地-井瞬变电磁数据中的可用数据用于反演；

3).将步骤2)的观测数据和对步骤1)中地-井模型进行非结构有限元正演得到的预测数据做数据拟合，根据观测电磁数据、正演得到的预测电磁数据以及对地-井模型的反演参考模型构建正则化反演的目标函数为：

其中，m₀和m分别是M维的反演参考模型和解向量，d_pre和d_obs是N_d×n_channel维的预测电磁数据和观测电磁数据，N_d为接收点的个数，n_channel为时间道个数，λ为正则化因子，W_d为数据方差矩阵，W_m模型粗糙度算子；

将方程式(9)两侧对模型m^n-1求导，获得第n^th迭代中目标函数的梯度gⁿ：

其中，J是雅克比矩阵，而

4).对正演方程两侧求导，其中预测数据对模型参数的偏导数为灵敏度矩阵；对伴随正演求解方程得到灵敏度矩阵的转置与向量的乘积；

5).将灵敏度矩阵的转置与向量的乘积结合预测数据与观测数据的拟合差，可计算目标函数的梯度；

6).采用拟牛顿法或高斯-牛顿方法计算模型更新量，得到新的迭代模型；

7).重复执行步骤3)～步骤6)，直到达到最大迭代次数或者符合终止条件，得到最终的反演参考模型；

依据该预测电磁数据和正则化反演的目标函数实现地-井瞬变电磁的迭代反演时，通过对正则化反演的目标函数进行求解以获得预测电磁数据和观测电磁数据的拟合差，将拟合差小于设置阈值作为反演迭代约束条件，在不满足约束条件时，依据求解正则化反演的目标函数获得的模型解向量更新方案参考模型。

2.如权利要求1所述的基于非结构有限元方法的地-井瞬变电磁反演方法，采用非结构有限元方法为探测地质结构构建初始地-井模型，包括：

为探测地质结构构建几何模型，并将几何模型剖分成非结构的四面体网格单元，同时设置几何模型信息、传递数据信息以及控制参数，得到初始地-井模型。

3.如权利要求1所述的基于非结构有限元方法的地-井瞬变电磁反演方法，在对初始地-井模型进行正演计算时，将模型的电场全波方程表示成有限元形式为：

其中，t表示时间，矢量

表示各四面体网格单元棱边的待求解电场，介电常数矩阵元素

电导率矩阵元素

磁导率矩阵元素

以及源矢量元素分别为：

其中，上标e表示四面体单元，i、j分别表示四面体单元索引，ε^e表示介电常数，σ^e表示电导率，

分别表示四面体单元内电场矢量在棱边i和棱边j的插值基函数，V表示体积，

表示哈密顿算子，j_s(r,t)表示源电流密度；

对于电场全波方程，忽略位移电流可得到电场扩散方程，有限元形式的电场扩散方程简化并将整体单元进行组装得到：

采用隐式后推欧拉方法进行时间离散，即将一阶后推欧拉方法应用于方程式(6)，可得：

Duⁿ⁺¹＝Buⁿ-Δtsⁿ⁺¹ (9)

其中，上标n和n+1表示递推过程，Δt表示时间步长，D＝(B+ΔtC)；

采用齐次狄利克雷边界条件，将矩阵方程(7)在外边界处的电场设置为零，同时直流电场可利用节点有限元求解如下三维泊松问题得到，即

其中，j_s(r)表示源电流密度，

表示节点电势；

在求解得到方程式(8)中的电势后，利用导数关系求解出直流电场；

将边界条件和初始电场带入方程并求解即可得到正演结果。

4.如权利要求1所述的基于非结构有限元方法的地-井瞬变电磁反演方法，其特征在于，在对正则化反演的目标函数进行求解时，需要计算观测电磁数据的灵敏度和模型的粗糙度，同时基于对数参数的电导率转换函数将模型参数变化到转换域内，对反演结果进行电导率上下限约束。

5.如权利要求4所述的基于非结构有限元方法的地-井瞬变电磁反演方法，其特征在于，计算观测电磁数据的灵敏度时，将正演方程两端对第k^th个模型参数求导，可得：

其中，m_k表示模型参数，s_TEM表示源项，K表示系数矩阵，在地-井瞬变电磁正演中，预测电磁数据d从电场解e插值得到，即d＝Le,L＝L_tL_d，L_d将电场解插值到指定接收点位置的插值函数，L_t是将计算时间道的解插值到观测时间道上；

预测电磁数据的灵敏度写为：

其中，j^k是整体灵敏度矩阵中第k^th列元素，设定一个N×M维的矩阵G，N＝N_edge×n,n计算时间道的个数，即：

G＝{g₁,g₂,g₃,…,g_M} (13)

长度为N的向量g_k的表达式为：

由于正演过程中采用的是总场的方法，故有

但对于下阶跃波的电性源，初始场与模型介质的电导率相关，需要计算De_DC和Ce_DC关于模型参数的导数，从公式(12)和(13)可得：

J^T＝G^TK^-TL^T (15)

从目标函数梯度的表达式可知，需要计算灵敏度矩阵的转置J^T和一个向量的乘积，设相乘的向量为u，并定义：

w＝K^-TL^Tu (16)

则得到如下的伴随正演方程：

K^Tw＝L^Tu (17)

求解方程(17)后可得灵敏度矩阵的转置J^T和向量u的乘积：

J^Tu＝G^Tw。 (18)

6.如权利要求4所述的基于非结构有限元方法的地-井瞬变电磁反演方法，其特征在于，在计算模型粗糙度时，采用加权平方和方法对非结构的四面体单元网格进行近似离散，表示式为：

式中，

其中，Vi是第i个四面体的体积，N(i)是与四面体i共顶点的所有四面体网格个数。

Images