一种电性源地空瞬变电磁三维聚焦反演方法

技术领域

本发明属于地球物理时间域电磁勘探领域，涉及地面、半航空瞬变电磁数据三维反演方法技术，特别涉及一种电性源地空瞬变电磁三维聚焦反演方法。

背景技术

地空瞬变电磁法属于电性源瞬变电磁勘探方法的重要分支。长久以来，受瞬变电磁正演理论、数据量和计算时间限制等因素的影响，地面、半航空瞬变电磁实用化解释技术仍然以一维为主。自然介质中三维特性较为普遍，因此采用三维反演解释技术更合理。然而，在对数据资料解释技术三维化时，计算效率和模型解释精度不足。其中，模型结构还原以及电性界面识别是研究的难点。针对以上问题，本发明以三维正演理论基础为支撑，在现有常规三维反演研究的基础上，提升了计算效率以及适应大数据集的实际需求，改进了反演模型期望以提升电性界面分辨能力，这两项是目前瞬变电磁三维反演技术发展的关键所在，对反演算法的广泛实用化具有极大的意义。

发明内容

本发明的目的在于克服现有技术的不足，提供一种能够较好的还原模型真实电性特征，具有边界锐化的效果，弥补了当今地空瞬变电磁探测数据反演解释系统的缺陷的电性源地空瞬变电磁三维聚焦反演方法。

本发明的目的是通过以下技术方案来实现的：一种电性源地空瞬变电磁三维聚焦反演方法，包含以下步骤：

(1)选择工区数据文件，包括参数文件、观测数据文件d；

(2)读取参数文件，包括：接受线圈面积、响应时间道、线源长度、偏移距、电流峰值、背景电阻率、响应时间；

(3)根据参数设置初始模型，包括：三维网格剖分、物性填充；

(4)设置反演终止条件，包括：最大迭代次数I

(5)计算频率域三维正演响应，然后采用伴随正演法求解频域率雅可比矩阵，更新数据项海塞矩阵；并通过频率域响应转时间域响应获得时间域磁场瞬态响应，并计算均方根误差Rms1；

(6)通过G-S变换将频域率雅可比矩阵转换为时间域雅可比矩阵，通过高斯牛顿反演方法得到L1正则反演目标函数；

(7)采用迭代重加权最小二乘法计算L1正则反演目标函数中的L1正则项，然后采用共轭梯度法求解线性方程组获得模型修正量Δm；

(8)更新模型向量：m＝m+Δm；再次计算频率域三维正演响应，通过频率域响应转时间域响应获得时间域磁场瞬态响应，并计算均方根误差Rms2；

(9)计算均方根误差△Rms：△Rms＝|Rms1-Rms2|/Rms1，判断均方根误差是否小于最小拟合误差Rms，或者判断迭代次数i是否小于最大迭代次数I

(10)终止迭代，输出当前模型m。

所述三维正演响应的具体实现方法为：将频率域电场双旋度方程通过伽辽金法改写为加权余量方程，方程通过不完全LU分解的双稳定共轭梯度法求解方程组得到电场；最后利用差分的方法实现辅磁场的计算；

单元格加权余量方程为：

其中，

将式(1)左端第一项表示为：

K

将左端第二项表示为：

为一阶插值函数，上标e表示三维网格中的某一个单元格，下标p和q表示棱边编号，K

右端项表示为：

对所有单元网格排序和组装则能构成一个大型稀疏线性方程组，通过解方程组能够得到全网格棱边电场分量，其表达式如下：

KE

其中K为总体刚度矩阵，由所有单元左端项组成，通过公式(2)和(3)得到，E

最后通过对E

所述L1正则反演目标函数表示为：

式(6)展开后表示为：

式(7)右端第一项为数据项目，第二项为L1正则项，其中N

所述迭代重加权最小二乘法的方法为：式(7)中L1正则项按矩阵展开近似表示为：

加权矩阵V

其中ε为一个小值；

式(7)对模型增量求导，得到：

公式(10)左端中括号内为数据项海塞矩阵，J

通过公式(10)得到模型修正量Δm。

本发明的有益效果是：本发明提供的电性源瞬变电磁三维聚焦反演成像方法，针对电性源发射，地面接受或者空中无人机接收二次场信号的地空时间域电磁探测方法。本发明首次在半航空瞬变电磁三维反演领域引入迭代重加权最小二乘法解决L1正则零点不可导问题，同时反演能够较好的还原模型真实电性特征，具有边界锐化的效果，弥补了当今地空瞬变电磁探测数据反演解释系统的缺陷。

附图说明

图1为电性源地空瞬变电磁三维聚焦反演的流程图；

图2为测点观测区域网格剖分示意图；

图3为半航空瞬变电磁L1、L2正则三维反演立体图；

图4为某测线不同工作结果对比；

图5为墨矿接地线源瞬变电磁数据L1、L2正则三维反演结果对比图。

具体实施方式

地空瞬变电磁三维聚焦反演方法，内容包括目标函数正则项的改进以及正则化因子和参考模型选取。本发明在三维反演技术选型中采用了具二阶收敛性的高斯牛顿法，其总目标函数海塞逼近阵可视为数据拟合项目的逼近海塞和精确正则项海塞之和，有利于实现自动变权重的鲁棒反演或聚焦反演策略。为增强反演模型对电性突变面的分辨能力，本发明将目标函数正则项采用L1范数，在高斯牛顿迭代框架下引入迭代重加权最小二乘策略近似求解，其迭代稳定性优于直接对L1范数逼近式求极小或次梯度搜索类的最优化策略。数据拟合的高斯牛顿型海塞逼近利用了雅可比矩阵内积形式，采用伴随正演法先并行求解频率域雅可比矩阵变换到时间域。为尽量降低反演算法的初始模型依赖性，以大正则因子反演得到的光滑模型为小正则因子反演的初始模型策略。基于上述研究，对数据进行反演能够较好的逼近模型的真实电性结构。下面结合说明书附图和具体实施方式对本发明作进一步的具体说明。

如图1所示，一种电性源地空瞬变电磁三维聚焦反演方法，包含以下步骤：

(1)选择工区数据文件，包括参数文件、观测数据文件d；

(2)读取参数文件，包括：接受线圈面积、响应时间道、线源长度、偏移距、电流峰值、背景电阻率、响应时间；

(3)根据参数设置初始模型，包括：三维网格剖分、物性填充；

(4)设置反演终止条件，包括：最大迭代次数I

(5)计算频率域三维正演响应，然后采用伴随正演法求解频域率雅可比矩阵，更新数据项海塞矩阵；并通过频率域响应转时间域响应(G-S变换)获得时间域磁场瞬态响应，并计算均方根误差Rms1；

三维正演响应的具体实现方法为：将频率域电场双旋度方程通过伽辽金法改写为加权余量方程，方程通过不完全LU分解的双稳定共轭梯度法求解方程组得到主场(频率域电场)；最后利用差分的方法实现辅助场(频率域磁场)的计算，将得到的总的频率域磁场通过G-S变换可以得到时间域磁场瞬态响应dBz/dt；

单元格加权余量方程为：

其中，

将式(1)左端第一项表示为：

K

将左端第二项表示为：

为一阶插值函数，上标e表示三维网格中的某一个单元格，下标p和q表示棱边编号，K

右端项表示为：

对所有单元网格排序和组装则能构成一个大型稀疏线性方程组，通过解方程组能够得到全网格棱边电场分量，其表达式如下：

KE

其中K为总体刚度矩阵，由所有单元左端项组成(公式2和3)，E

(6)通过G-S变换将频域率雅可比矩阵转换为时间域雅可比矩阵，通过高斯牛顿反演方法得到L1正则反演目标函数，L1正则反演目标函数表示为：

式(6)展开后表示为：

式(7)右端第一项为数据项目，第二项为L1正则项，其中N

(7)采用迭代重加权最小二乘法计算L1正则反演目标函数中的L1正则项，然后采用共轭梯度法求解线性方程组获得模型修正量Δm；

迭代重加权最小二乘法的方法为：式(7)中L1正则项按矩阵展开近似表示为：

加权矩阵V

其中ε为一个小值，这里取1e-3。

加权矩阵V

式(7)对模型增量求导，得到：

公式(10)左端中括号内为数据项海塞矩阵，J

通过公式(10)得到模型修正量Δm。反演过程即将求解出的Δm和m相加后得到一个修正后的新模型向量，再通过反演步骤(9)来判定结果是否满足条件，直到对m修正后的结果满足步骤(9)的终止条件后，结束反演过程，输出最终得到的模型向量。

(8)更新模型向量：m＝m+Δm；再次计算频率域三维正演响应，通过频率域响应转时间域响应(G-S变换)获得时间域磁场瞬态响应，并计算均方根误差Rms2；

(9)计算均方根误差△Rms：△Rms＝|Rms1-Rms2|/Rms1，判断均方根误差是否小于最小拟合误差Rms，或者判断迭代次数i是否小于最大迭代次数I

(10)终止迭代，输出当前模型m。

配套上述方法技术的建立对应电性源地空瞬变电磁三维聚焦反演成像软件模块包括：系统功能模块和底层支撑模块，其中，系统功能模块包括数据文件管理模块、预处理模块、正演成库模块、电磁响应值搜索和成图模块，底层支撑模块包括数据文件IO模块、嵌入式数据库模块、通用数学库模块；底层支撑模块，用于提供通用的功能函数给系统功能模块。

实施例1：规则模型数据三维聚焦反演

理论模型的背景电阻率设置为100Ω·m，在其中嵌入20Ω·m的低阻异常体，该异常体大小为300m×300m×100m，顶面埋深为200m，测点观测区域网格剖分示意图(不含扩展区域)如图2所示。图中，(a)为XY切面的测点，(b)为XZ(YZ)切面的测点；浅色为100Ω·m，深色为20Ω·m。总网格剖分为23×23×21，其中目标区域水平面x、y方向分别9网格共计81个测点，扩展区域x、y正负方向分别设置7个网格扩展至5km以外，空中设置6个网格扩展至3km以上，地下设置6个网格扩展至6km以下。设计线源长度1000m,电流大小为1A，以坐标原点作为线源中心点，线源到观测区的最小偏移距为300m，接收线圈高度分别为地面0m(地面瞬变电磁)和空中50m(半航空瞬变电磁)。

图3为半航空瞬变电磁L1、L2正则三维反演立体图，其中，(a)为L1正则，切除表层的俯视角；(b)为L2正则，切除表层的俯视角；(c)为L1正则，切除一角的侧视角；(d)为L2正则，切除一角的侧视角。从图中可以看出，三维聚焦反演能够还原模型真实电性特性。

实施例2：某地矿产实测数据三维聚焦反演

该地区地面瞬变电磁测线布置总体规则，线源长度2km，共5条测线，测线间距200m，测点间距80m。三维模型目标区域根据测线布置，按照规则六面体进行剖分，外围扩展10km。

图4为自然电场法电位曲线特征、已知矿体、三维聚焦反演结果对照，其中，上：自然电位曲线；中：地质勘查结果；下：L1正则三维反演结果。图5为墨矿接地线源瞬变电磁数据L1、L2正则三维反演结果对比图，左：L1正则，右L2正则。从图中可以看出本发明所提出的反演算法结果在水平位置上与已知矿体位置以及传统电法勘探一致。

综上所述，自然电场资料也与瞬变电磁数据的三维反演结果具有较好的对应关系。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。