定源瞬变电磁数据的时频转换反演技术研究

2021-11-16张艳龙吴燕清邓刘洋崔少北

矿业安全与环保 2021年5期

张艳龙，吴燕清，邓刘洋，崔少北

(1.重庆大学煤矿灾害动力学与控制国家重点实验室，重庆 400044； 2.中煤科工集团重庆研究院有限公司，重庆 400037)

瞬变电磁法通过不接地回线以脉冲电流为场源，激励探测体感应二次电流，在脉冲间歇期间测量二次场随时间变化的响应[1-2]。采集并处理该响应信号，根据衰减异常信息，分析目标体电性特征[3]。

采用大定源回线装置测量响应信号是常用的工作方式，受限于数据源与处理方式，在实际应用中多以一维反演及相对视电阻率成像为主要成果体现。要实现三维正反演并成熟地在地面勘查中应用仍需进行大量的研究工作[4-6]。频域电磁法三维正反演技术日趋成熟[7-8]，如何通过数据转换进行三维正反演已成为热点问题。以傅里叶变换为基础进行数据转换[9]，形成了时频转换理论，为瞬变电磁数据的三维正反演提供了一种新方法[10-11]。

基于时频转换理论，笔者通过正则化技术，结合有限差分法与NLCG反演算法，在采用大定源回线装置进行瞬变电磁探测时，进行了理想模型及实测瞬变电磁数据的频域正反演研究，论证了三维反演的可行性，同时验证了其高分辨能力。

1 时频转换

从时域信号中获取频域响应时，时域信号需满足采样定律。只有具备时域和频域电磁响应等价条件，才能实现相互转化。瞬变电磁仪器是依据时间进行数据采集的，其采样范围有限，数据不满足傅里叶变换要求，“混频效应”会使频域响应的可信度降低。为解决以上问题，直接对傅里叶变换分段积分进行离散，得到矩阵形式。采用最小二乘法，并对频域响应加入光滑约束，完成时域到频域的转换。

1.1 关断时间校正

时频转换理论以阶跃发射波形为基础，而野外实测的瞬变电磁数据的发射波通常是具有一定关断时间的斜阶跃波，因此在对时域信号进行频域转换之前要对关断时间进行校正，得到理想阶跃关断条件下的瞬变电磁数据。关断时间通过最小二乘拟合后沿估计矫正。

1.2 时频转换矩阵计算方程

瞬变电磁的时、频域响应可以相互转换，其具体步骤为：傅里叶变换→欧拉变换→折线逼近法离散→矩阵变换→正则化求解。

傅里叶变换：

(1)

式中：h(t)为时域电磁响应；t为时间；I(ω)为电流波形的频谱；ω为角频率，ω=2πf；H(ω)为频域电磁响应。

在瞬变电磁法中，最常用的是下阶跃波形，故式(1)可变换为式(2)：

(2)

通过欧拉公式，推导频域—时域电磁场的正弦关系：

(3)

式中：Re为取实符号，表示频率域谐变电磁场的实分量；h(0)为零时刻磁场的初始值[12]。

对式(3)用折线逼近法进行离散[13]得到式(4)：

[sin(ωi+1t)-sin(ωit)]

(4)

式中：Fs(ω)=Re[Hz(ω)]/ω；Hz(ω)为频域电磁响应。

将式(4)用矩阵形式表达，令N离散的频率响应Fs构成矢量D，由M个测道组成的时域响应矢量为d，对于第ti道响应di有：

(5)

令系数：

进一步令：

(6)

则式(5)可写为：

用矩阵表示为：

d=LD

(7)

式中：Dj为第j个离散频率响应构成的矢量；L为系数矩阵。

从式(5)和式(6)可以推出矩阵系数L：

(8)

其中矩阵M可以通过式(6)求出，只与余弦和角频率有关。

1.3 正则化技术

为了提高分辨率，采用有利于解决不适定问题且收敛性更好的正则化反演方法[13]。由于矩阵L的奇异特性导致不能直接求解式(5)，解决该问题的经典方案是采用正则化技术，通过对解加以适当约束以得到稳定的解[14]，将矩阵方程计算转化为最优正则化问题，见式(9):

(9)

式中:WD为模型光滑矩阵;λ为正则化参数。

最优正则化参数的选择方案很多，这里采用的是L-曲线法求取最优正则化参数。L-曲线法是一种基于数据误差水平未知的启发式选取正则化因子λ的方法[15]。

1.4 数据计算验证

为验证时频转换的精度与可行性，进行三维理论模型正演数据及工程勘察数据的转换，对比响应特征。

1)三维模型试验

在一个电阻率为150 Ω·m均匀半空间，发射回线边长分别为400、400 m,在框内中心位置下方 100 m 安设一个三维高阻异常体，异常体尺寸为 80 m ×80 m×80 m，发射电流3 A,三维模型布设如图1 所示。

图1 三维模型布设平面图

图1中外部虚线框为发射回线，中部实线框为接收范围，内部有剖面线的线框为高阻异常体边界，中心圆点为坐标原点，黑色圆点为瞬变电磁测点，线距为40 m,点距为20 m。

选取框内A、B两个测点的数据进行时频转换对比，其中A点位于高阻异常体附近，B点位于高阻异常体较远处。把在A、B点采集的时域信号通过时频转换矩阵计算方程进行时频转换、对比。三维模型数据时频转换对比如图2所示。

(a)A点时域信号图

(b)A点垂直磁场响应和时频转换对比图

(c)B点时域信号图

(d)B点垂直磁场响应和时频转换对比图

图2中A点、B点转换计算的正则化因子λ分别为1.15×10-4、2.97×10-3。

数值试验讨论了时域到频率的计算效果。对时域垂直磁场响应进行了计算，利用时频转换方法，将A、B两个测点的瞬变电磁数据转化为频域响应，由图2可知，不同形态、反映数据的时频转换曲线与理论值基本一致，表明利用正则化技术的时频转换理论是可行的，变换结果的正确性证明了时频转换方法的可靠性。

2)实测数据试验

本次瞬变电磁法试验采用40 m×20 m的网格，线距为40 m，点距为20 m。使用V8多功能电法仪中心回线装置，发射线框尺寸为280 m×280 m，发射电流为10 A，采样频率为25 Hz，单点采样时间 2 min，增益×4。接收区域位于发射线框中心120 m×120 m范围。选取试验点为L4测线160 m里程点，标记为实测点C。

实测数据时频转换对比如图3所示。

(a)感应电动势与时域磁场图

(b)感应电动势及其拟合曲线图

(c)感应电动势OCCAM反演模型图

(d)频域磁场对比图

图3中转换计算的正则化因子λ为4.87×103。

图3给出了28个时间窗口的感应电动势值，从图中可以看出实测曲线比较光滑，无需进行滤波等预处理，首个和最后一个时间道分别为2.320×10-5s和1.002×10-2s，在转换之前，得到了理想阶跃关断条件下瞬变电磁响应。对关断时间校正数据进行积分，得到时间域垂直磁场，原始数据和积分后的垂直磁场如图3(a)所示，拟合差小于4%。反演后得到层状模型如图3(c)所示，应用时频转换方法得到的频率域垂直磁场响应与图3(c)模型正演得到的频率域垂直磁场数据对比如图3(d)所示。

通过对图3 分析可知，时频转换曲线与理论值基本一致。理论模型数据及实测数据试验证明了本文采用的时频转换方法的有效性、稳定性及精确度，为瞬变电磁法由时域向频域转换，以及在频域的三维正反演提供了依据。

2 三维模型正演

三维瞬变电磁正演研究及算法较多，如积分方程法[16-17]、时域有限差分等，频域响应计算后转换成时域[18-20]。为验证时频转换技术的可行性，建立正演模型，开展正演研究，以期得到理论模型下的正演成果，为下一步反演提供数据。

2.1 基本原理

基于麦克斯韦方程组，地下空间内的电场、磁场信号可被视为离散参数，将三维模型解构为大型几何空间集合体，每一组响应构成一个单元。求解几何单元外棱交点上的电场信号，通过非线性插值方法及电磁感应原理，可求得模型中任何单点的频域电场、磁场信号参数，最后求得三维模型的电磁场响应。

基于有限差分法实现三维瞬变电磁正演时还需注意几个问题：

1)磁场的散度：约束计算与求解，约束条件—磁场散度等于0；

2)网格剖分：模型剖分、单元划分时，网格设置应满足狄利克雷边界条件。

2.2 三维模型验证

建立电阻率为150 Ω·m的均匀半空间，在深度100 m安设1个5 Ω·m三维低阻异常体，200 m的正方形发射线框敷设于水平地面，发射电流为3 A，设三维低阻异常体的中心点在水平地表的投影为三维坐标原点，z轴方向为沿地表垂直向下，主剖面为沿x轴的测线剖面。三维模型布设与观测方式如图4 所示。

图4 三维模型布设与观测示意图

设置三维网格剖分参数如下：沿坐标轴方向 50 m×50 m×30 m，空间网格大小30 m×30 m×30 m。主剖面上0.01、20 ms的正演数据时域响应(方形点)如图5所示。使用积分方程方法求得的时域信号用曲线绘制[21]。

图5 三维模型正演数据时域响应曲线

由图5可知，该模型的正演试验表明，基于有限差分法的三维模型正演数据时域响应与Newman等[12]的结论基本一致，证明笔者采用的三维正演算法可靠并具有高精度，能为三维模型反演的研究提供依据。

3 三维模型反演

频域电磁法数据三维反演主要有拟牛顿法、高斯牛顿法、共轭梯度法和非线性共轭梯度法(NLCG)[22-24]。NLCG算法是非线性反演中确定性、统计性综合的方法，可将复杂计算线性化。通过精密的网格剖分与地质建模，基于函数的梯度，不断迭代计算，求得三维模型最优解。

在视电阻率150 Ω·m的均匀半空间中设计2个10 Ω·m的低阻异常体，埋深为100 m，其三维形态为100 m正方体，地面敷设2个400 m方形发射线框，发射电流为3 A，测区边长为280 m，测线共8条，每条测线布置10个测点，测框布置如图6所示。

图6 测框布置示意图

三维模型反演网格剖分密度为30 m×25 m×20 m，共计反演迭代计算120次，三维模型反演成果切片如图7所示。

图7 三维模型反演成果切片示意图

由图7可知，三维模型反演成果的低阻异常体范围与理论模型高度一致，表明NLCG反演算法可实现1条测线上不同发射源、不同测区的数据反演，同时证明NLCG反演算法在三维模型反演时的可行性并验证了其精度。

4 工程数据三维反演

4.1 地质概况

勘探井田位于太行山脉西部沁水盆地南部边缘，为低山丘陵地带。井田内地表大部分被第四系黄土覆盖，局部基岩出露。石炭系上统太原组与二叠系下统山西组为本勘探区主要含煤地层。

根据勘探区内钻孔地质信息，得到该勘探区煤、岩层电性特征，如表1所示。

表1 地层电性特征

上覆风化砂泥岩，二叠系分布有砂岩、泥岩、煤层；石炭系分布有石灰岩、煤层、砂岩、铁(铝)岩；奥陶系以泥灰岩、石灰岩为主。

由表1可以看出，该区地层电性特征有较大区别，上覆风化层及泥岩视电阻率值较低；二叠系底部、石炭系顶部煤层地层表现为较高视电阻率；奥陶系又呈高阻反映。表明勘探区是典型的三层H 模型。

4.2 数据处理

瞬变电磁法试验采用40 m×20 m的网格，线距为40 m，点距为20 m。使用V8多功能电法仪中心回线装置，发射线框尺寸为280 m×280 m，发射电流为10 A，采样频率25 Hz，单点采样时间为2 min，增益×4。接收测框尺寸为120 m×120 m，接收面积小于1/3发射面积，测区内共布设南北向12条测线，工程布置如图8所示。

图8 工程布置图

对在第5线300 m里程点采集的数据进行单点反演，该试验点衰减曲线如图9所示(B为磁感应强度)。结合附近钻孔资料，此处3#煤层埋深为98 m。

图9 试验点衰减曲线

试验点一维反演成果如图10所示。

图10 试验点一维反演成果

由图10可知，随着埋藏深度增加二叠系(P)地层视电阻率逐步升高，进入石炭系(C)地层后视电阻率进一步升高，到达奥陶系(O)后视电阻率依然表现为高电阻率值特征。与钻孔资料对比可知，反演结果基本反映出了钻孔附近地层的电性特征。

试验点数据表明,该次瞬变电磁探测参数合理、数据质量较高。将原始数据滤波去噪及关断时间校正后，利用正则化技术进行了本区数据的时频转换。

4.3 频域三维反演

经数据处理得到频域数据，在反演前，结合地质信息建立约束模型。本次反演依据钻孔地层信息进行物性层状建模，深度方向设置不同物性层状岩层共计15层，网格剖分密度为20 m×20 m×20 m，边界扩充系数为1.5，满足狄利克雷边界条件。全部433个测点数据均参与三维反演，迭代拟合差RMS=1.25，反演深度为400 m。三维反演垂直切片成果如图11 所示。