一种时间域电磁数据反演成像方法

技术领域

本发明涉及地球探测与信息技术领域，特别是一种时间域电磁数据反演成像方法。

背景技术

航空电磁数据成像经历了几十年的发展历程，成像方法多、成像理论较为成熟，在航空地球物理界获得广泛应用。航空电磁数据成像是将观测数据(电磁响应)转换为表征地下介质电性分布特征的中间参数，如视电导率、视深度等。成像算法速度快，能从海量航空电磁数据中快速提取地下电性主要信息，适用于现场快速数据处理，同时它也可为复杂的航空电磁反演提供初始模型。常用的航空电磁成像方法有Sengpiel差分视电阻率法、电导率深度成像、基于浮动薄板理论的成像。

这些算法在内存需求和计算速度上可以满足三维电磁反演的需要。但对于大规模的三维电磁数据反演，为节省内存和计算时间，通常灵敏度矩阵非显式地计算，即通过伴随方法计算灵敏度矩阵和向量的乘积。与高斯牛顿方法相比，拟牛顿法和非线性共轭梯度法每次迭代正演计算量少，但反演收敛速度较慢。

发明内容

本发明的目的在于克服现有技术的不足，提供一种时间域电磁数据反演成像方法，对航空电磁法瞬变电磁快熟反演成像效果十分明显，通过多次迭代，大大降低拟合误差，反演收敛速度较快。

本发明的目的是通过以下技术方案来实现的：

一种时间域电磁数据反演成像方法，包括以下步骤：

S1：将电磁测量数据进行CDI成像，得到地层初始模型；

S2：根据Tikhonov正则化方法建立非线性目标函数：

式中，J(m)为目标函数，m表示模型向量，m

S3：利用全波形灵敏度矩阵

m

式中，m

进一步，所述步骤S2包括以下子步骤：

S201：建立电场耦合势摄动方程：

式中，

S202：求解频率域电磁场Frechet导数：

(4)式中，第一个式中，δE(r

S203：进行频率-时间转换：

式中，b(t)是时间域函数，t为时间参数，公式等号右边为基于傅里叶变换的时频转换计算；

S204：结合电磁测量数据，建立非线性目标函数：

式中，J(m)为目标函数，m表示模型向量，m

进一步，所述步骤S3包括以下子步骤：

S301：建立线性化法方程：

式中，

S302：完成Cholesky分解的共轭梯度迭代计算：

式中，m

本发明的有益效果是：

本方法对航空电磁法瞬变电磁快熟反演成像效果十分明显，通过多次迭代，大大降低拟合误差，反演收敛速度较快。

附图说明

图1为本发明的流程图；

图2为梯形波dBz(t)/dt响应及其相对误差；

图3为不同测点上梯形波dBz/dt响应；

图4为地层原始模型；

图5为第1次迭代后的反演模型；

图6为第7次迭代后的反演模型；

图7为第14次迭代后的反演模型；

图8为迭代反演的拟合误差。

具体实施方式

以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需说明的是，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。

需要说明的是，以下实施例中所提供的图示仅以示意方式说明本发明的基本构想，遂图式中仅显示与本发明中有关的组件而非按照实际实施时的组件数目、形状及尺寸绘制，其实际实施时各组件的型态、数量及比例可为一种随意的改变，且其组件布局型态也可能更为复杂。

实施例一：

如图1至图8所示，一种时间域电磁数据反演成像方法，包括以下步骤：

S1：将三维航空电磁测量数据进行CDI成像，得到地层初始模型；

对三维航空电磁测量数据需进行预处理；

S2：根据Tikhonov正则化方法建立非线性目标函数：

式中，J(m)为目标函数，m表示模型向量，m

S3：利用全波形灵敏度矩阵

m

式中，m

这样将反演成像转化为含有稳定泛函的非线性泛函的极小化问题。

所述步骤S2包括以下子步骤：

S201：建立电场耦合势摄动方程：

式中，

S202：求解频率域电磁场Frechet导数：

(4)式中，第一个式中，δE(r

S203：进行频率-时间转换：

式中，b(t)是时间域函数，t为时间参数，公式等号右边为基于傅里叶变换的时频转换计算；

S204：结合电磁测量数据，建立非线性目标函数：

式中，J(m)为目标函数，m表示模型向量，m

所述步骤S3包括以下子步骤：

S301：建立线性化法方程：

式中，

S302：完成Cholesky分解的共轭梯度迭代计算：

式中，m

在迭代反演过程中，利用Morozov偏差原理、以及不完全Cholesky分解与共轭梯度(ICCG)法等优化方法，实现分区数据反演，并行全波高阶逼近技术，大大降低反演迭代次数，自适应正则化保证反演质量。

通过并行运算与信号开窗法对整个时间域电磁资料进行分区域，同时迭代反演，并对反演结果进行组合，利用OpenMP多核并行，实现了多频率宽带电磁信号散射的快速求解，得到三维反演成像数据。

采用“moving footprint”技术进行正演加速，即依次对每个测点提取其下方footprint大小的模型重建子模型，并在子模型内完成相应正演计算，形成局部反演模式，提高机载电磁复合效应三维快速反演成像计算的效率。

通过上述全波形自适应正则化迭代计算方法，实现三维电阻率反演成像。

本实施例基于三维时间域电磁高效仿真算法研究，建立一套二维、三维显式时间域航空电磁灵敏度矩阵与反演成像方法，实现机载时间域电磁数据反演成像解释。利用差异场理论与并行计算实现二维、三维全波形电磁响应与灵敏度矩阵快速计算，并结合全波形自适应正则化反演技术，用于开发相应的反演成像算法。

采用正演模拟结果进行验证，地层模型参数选择背景电阻率，异常体电阻率，异常体大小100m*200m*200m，异常体中心离地面深度200米。飞行高度30米。发射线圈为边长为6.66m的正八边形，接收点位于正八边形中心。基于面积相等原理，用半径为8.26m的圆形发射线圈近似边长6.66m的正八边形。梯形波上升下降沿时间均为0.2ms，稳定时间为3.6ms，最大电流强度。频率采样：[1.0e-3，1.0e8]，采样率：对数等间距67个采样点。计算时间：计算时间与网格数量有关。当网格数量为50*50*70时，计算441个采样点(21条测线、每条测线采样间距为10米)所需时间约为12小时；当网格数量为30*30*60时，大约时间约为1.5小时。结果网格数量为50*50*70，中心区域x,y,z方向步长为10米。图2所示，数值结果与文献结果(地球物理学报，2017年，Vol.60，No.1：369-382)对比，最大相对误差低于5％。

本方法对航空电磁法瞬变电磁快熟反演成像效果十分明显：图2是埋有不同深度的两个异常体的地层模型，电阻率均假定。其中，一个长方形异常体，埋深60米，长宽高分别为80米、80米、60米；另外一个圆柱体异常体，埋深100米，半径为40米，高度100米。图3是二维平面不同测点的磁场分布，从图像中可以看出，正演模拟结果能够明显体现出异常体的水平方向的大致位置，异常体埋深越浅，磁感应强度越明显。选定一个理论模型(图4)，利用正演模拟软件对其进行模拟，得到不同测点不同时刻的正演模拟曲线。设定半空间模型为初始模型，利用编制的三维模拟软件，地层模型进行反演，得到初步的反演结果(图5-图7)，经过14次的迭代反演，最终拟合误差下降至20％左右(图8)。

以上所述实施例仅表达了本发明的具体实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。