一种锐边界模型的瞬变电磁反演方法与流程

本发明涉及地球物理探测技术领域，尤其地涉及一种锐边界模型的瞬变电磁反演方法。

背景技术：

瞬变电磁法作为一种地球物理勘探的有力手段，能够解决矿产勘查、能源、工程、水文、环境地质调查、考古探测等多种地球物理问题。它利用阶跃波或其他脉冲电流地下发射一次脉冲磁场，地下导电地质体在一次场激发下产生感应涡流，进而激发出二次磁场。当发射电流关断时，一次场消失，地下涡流随时间衰减，其衰减时间与导电地质体的电性参数(体积、结构、电阻率、埋深)有关。通过观测二次场值，可以重建地下不均匀体的电性参数，这个过程称为反演。反之，利用已知电性参数的地下结构计算二次场变化的过程，称为正演。

常用的瞬变电磁一维反演方法有Marquardt方法和Occam方法。Marquardt方法通过假设模型的正演数据与实测数据直接进行拟合，迭代求出大地模型的每层厚度及其电阻率。该方法计算速度快，原理简单，结果能反映大地电阻率的突变边界，但反演过程不稳定，易陷入局部最优解。Occam方法在Marquardt方法约束函数的基础上施加了纵向平滑约束，使纵向电阻率梯度的二范数最小，采用固定的大地模型厚度，利用高斯牛顿法迭代求解大地每层电阻率。该方法反演过程稳定，能够找到符合纵向平滑约束的全局最优解，但对某些电阻率突变的地电界面分辨率不佳，无法准确刻画大地的分层结构。

技术实现要素：

(一)要解决的技术问题

本发明提供了一种锐边界模型的瞬变电磁反演方法，用以解决Marquardt方法反演不稳定、Occam方法尽可能对大地电阻率变化进行平滑，对电阻率突变的分层边界分辨率不佳，无法准确刻画大地的分层结构以及上述两种方法混合使用在实际工程应用中不够简便的问题。

(二)技术方案

本发明提供了一种锐边界模型的瞬变电磁反演方法，所述反演为通过实测二次场数据，重建不均匀体的电性参数，该瞬变电磁反演方法包括构建反演初始模型，其特征在于还包括如下步骤：

根据所述反演初始模型的待反演电性参数构建最小梯度支撑泛函；以及

根据最小梯度支撑泛函构建反演目标函数。

上述方案中，还包括：

正演所述反演初始模型，得到理论响应值；

根据所述实测二次场数据与所述理论响应值构建不匹配泛函。

上述方案中，所述构建反演目标函数还包括：

根据所述不匹配泛函构建反演目标函数。

上述方案中，使用如下公式构建所述最小梯度支撑泛函为：

式中，P_s为最小梯度支撑泛函，β称为聚焦因子，是一个远小于1的正数，m为待反演的电性参数，R为描述梯度变化的矩阵：N为反演层数。

上述方案中，使用如下公式构建所述反演目标函数为：

P＝P_t+P_s

式中P为反演目标函数，P_t为不匹配泛函，α为正则化因子，用于调整反演过程中不匹配泛函和最小梯度支撑泛函之间的比重，P_s为最小梯度支撑泛函。

上述方案中，还包括：

以初始模型为初始条件，求解反演目标函数取最小值时对应的电性参数，得到反演结果。

(三)有益效果

本发明提供的方法，可以得到以下有益效果：

1、本发明方法通过使用正则化目标函数，继承了Occam反演稳定的优点，反演结果不依赖初始模型，能够找到全局最优解。

2、本发明方法选取最小梯度支撑泛函作为电阻率约束条件，产生类似Marquardt反演的结果，能突出电阻率突变的地层边界，弥补了Occam反演对电性界面反映不清楚的缺陷。

3、本发明方法给出了使用过程中每一步骤的具体公式，操作简便，能一步实现Occam反演与Marquardt反演混合使用的反演结果。

附图说明

图1示意性示出了根据本发明实施例的锐边界模型的瞬变电磁反演方法的流程图。

图2示意性示出了根据本发明实施例的瞬变电磁二次场相应数据。

图3示意性示出了根据本发明实施例的锐边界模型反演结果。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明作进一步的详细说明。

图1示意性示出了根据本发明实施例的锐边界模型的瞬变电磁反演方法的流程图。如图所示，

在步骤S1，输入反演参数，反演参数包括发射电流I，发射线圈半径R，发射线圈匝数N_t，观测偏移距r，接收线圈有效面积S。根据本发明实施例，发射线圈匝数为N_t＝1，发射电流为I＝5A，发射半径R＝50m，观测点偏移距r＝0m，接收线圈有效面积为2000m²。

在步骤S2，构建反演初始模型，反演初始模型包括反演层数N，各层厚度h＝[h₁，h₂，...，h_N-1]，反演初始电阻率m＝[m₁，m₂，...，m_N]。由于随深度变化反演分辨率降低，通常令h_i+1/h_i＞1。根据本发明实施例，反演模型为30层，厚度增长比例h_i+1/h_i＝1.12，第一层厚度h₁＝5m，每层的初始电阻率为100Ωm。

在步骤S3，正演初始模型得到理论数据。根据本发明实施例，正演所述反演初始模型，得到理论响应F(m)。

在步骤S4，输入实测的二次场数据，输入由瞬变电磁法观测的二次场数据，包括采集时间T及其对应的场值d_obs。其中，场值d_obs要满足衰减特性，需剔除早期包含一次场的部分和晚期由于噪声产生震荡的部分。根据本发明实施例，如图2所示，该数据由三层大地模型正演求得，其厚度分别为100m，200m，均匀半空间。相应的电阻率分别为300Ωm，90Ωm，300Ωm。

在步骤S5，构建实测的二次场数据与理论响应值的不匹配泛函。根据本发明实施例，不匹配泛函可写为：

P_t＝||W_d(F(m)-d_obs)||²

式中，W_d为描述数据比重的权值矩阵，F(m)是理论响应值，W_d可选择为：

式中，σ_i为各时间测道的噪声，M为时间测道总数。

根据本发明实施例，在此假设各时间测道数据所占比重一致，即权值矩阵为单位对角阵。

在步骤S6，构建最小梯度支撑泛函。根据本发明实施例，构建约束电阻率变化的最小梯度支撑泛函，可写为：

式中，β称为聚焦因子，是一个远小于1的正数。R为描述梯度变化的矩阵，可写为：

根据本发明实施例，为了突出大地电阻率的突变边界，β应尽量小，取β＝5×10^-11。

在步骤S7，由不匹配泛函和最小梯度支撑泛函构建反演目标函数，根据本发明实施例，反演目标函数可写为：

P＝P_t+αP_s

式中，α称为正则化因子，用于调整反演过程中不匹配泛函和最小梯度支撑泛函之间的比重，根据本发明的实施例，设置α的搜索范围为10^-4-10⁰。

在步骤S8，以初始模型为初始条件，采用高斯牛顿法求解目标函数取最小时对应的电阻率，得到反演结果。将理论响应值F(m)近似写为：

F(m)＝F(m₀)+J₀(m₁-m₀)

式中，m₀为某步迭代的初始条件，m₁为本次迭代的待求电阻率，J₀为某步迭代时正演函数对初始条件的导数矩阵，其每个元素为：

将F(m)＝F(m₀)+J₀(m₁-m₀)代入目标函数P＝P_t+αP_s中，令得迭代公式：

每步迭代后用新求得的结果更新初始条件，当拟合差小于期望拟合差时，迭代终止。每步迭代过程需在一定区间内动态选择α值，使目标函数达到最小。本发明实施例选择但不仅限于高斯牛顿法。

在步骤S9，得到最终的反演结果，即不均匀体的电性参数，结果如图3所示，可以看到，锐边界反演结果与理论模型有较高程度的拟合，对电阻率突变边界有较高的分辨率。

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