CN103760614A - 一种适用于不规则发射波形的瞬变电磁正演方法 - Google Patents

Abstract

本发明公开了一种适用于不规则发射波形的瞬变电磁正演方法，其特征在于，包括：建立圆柱坐标系，并建立接收点处的垂直磁场分量频率域表达式；将所述频率域表达式表示为Hankel积分形式，求出垂直磁场分量的频率域响应大小；根据所述垂直磁场分量的频率域响应大小采用GST算法求取垂直磁场分量的时间域响应大小；根据所述垂直磁场分量的时间域响应大小建立不规则发射波形信号在大地层状介质上的瞬变电磁响应表达式；利用所述瞬变电磁响应表达式求取每个时间段内的瞬变电磁响应，并将每个时间段内的瞬变电磁响应进行累加得到不规则发射波形的瞬变电磁响应。本发明将可用于地面、海洋和航空瞬变电磁法的系统预研、数据处理和解释。

Description

一种适用于不规则发射波形的瞬变电磁正演方法

技术领域

本发明涉及电磁学技术领域，尤其涉及一种适用于不规则发射波形的一维正演瞬变电磁(TEM)方法。

背景技术

时域电磁学是当今电磁学领域非常活跃并在蓬勃发展中的一个前沿分支，在通信、计算机网络、电子工程、雷达、地物探测、电磁兼容和生物电磁学等领域发挥着重要作用。其中，层状介质上展射电源的瞬变电磁响应是该课题的重点之一，是雷达与低频地球物理探测方法的常用物理模型。

在地物探测的瞬变电磁法实际作业中一方面由于仪器内部器件功能限制，地面瞬变电磁(GTEM)系统的发射波形也并非理想的斜变信号，而是类指数衰减形式。但在传统的正反演和数据处理中，波形则被看作理想的斜变信号，然后通过对阶跃响应进行关断时间上的平均求得。另一方面随着航空瞬变电磁(ATEM)系统的研制和应用，发射波形不再局限于单纯的斜变信号，半正弦、锯齿、梯形信号等往往成为首选，因此发射信号的时间展源效应必须考虑在内。

在不规则波形瞬变电磁响应的求取方面，国内外的相关文献往往是先在频域卷积然后变换到时域而得到。但对于某些复杂信号，数学上得到其频谱信息是非常困难的，因此传统方法往往只适用于简单信号响应的求取。

可见，如何对于目前瞬变电磁系统的复杂多变的波形，高效而准确的求取到其响应，成为瞬变电磁法系统预研、正反演解释的重要和关键技术之一。

发明内容

本发明的目的是提出一种适用于不规则发射波形的一维瞬变电磁正演方法，在满足一定的响应精度的同时，提高运算速度，便于该算法在实际工程中的应用。

根据本发明一方面，其提供了一种适用于不规则发射波形的瞬变电磁正演方法，其特征在于，包括：

步骤1、以圆形发射回线中心点对应的地面处为原点，建立z轴朝向地下的圆柱坐标系，并建立接收点处的垂直磁场分量的频率域表达式；

步骤2、将所述垂直磁场分量的频率域表达式表示为Hankel积分形式，并将其转换为滤波器积分形式，并将转换后的滤波器积分形式内的贝塞尔函数使用预先设定的滤波系数代替，并求出垂直磁场分量的频率域响应大小；

步骤3、根据所述垂直磁场分量的频率域响应大小采用GST算法求取垂直磁场分量的时间域响应大小；

步骤4、根据所述垂直磁场分量的时间域响应大小建立不规则发射波形信号在大地层状介质上的瞬变电磁响应表达式；

步骤5、利用所述瞬变电磁响应表达式求取每个时间段内的瞬变电磁响应，并将每个时间段内的瞬变电磁响应进行累加得到不规则发射波形的瞬变电磁响应

本发明直接从时间域出发，以阶跃响应为基响应，通过阶跃响应叠加的方法求取了任意波形的响应。本发明提出的上述方法是一种在时间域内将发射波形非均匀采样，进而进行阶跃响应叠加来实现任何发射波形的瞬变电磁响应；它将适用于各种瞬变电磁系统的任意发射波形，具有运算简单、高效率和适应性强等优点。

附图说明

图1是本发明中一种适用于不规则形状发射波形的一维瞬变电磁正演方法流程图；

图2A是使用本发明的高指数信号对实施例的处理结果与未用本发明的斜变处理结果对比图；

图2B是使用本发明的高指数信号对实施例的瞬变电磁响应与未用本发明的瞬变电磁响应对比图。

具体实施方式

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

TEM仪器在实际工作中，由于发射机电子器件特性以及回线感性负载等原因，发射波形并非理想的斜阶跃而为类指数等形状。传统方法在处理斜阶跃响应时，采用下阶跃响应在关断时间上的平均得到，因此传统方法并不能求取指数信号的响应。

因此，本发明提出了一种适用于不规则形状发射波形的一维瞬变电磁正演方法。图1示出了本发明提出的上述方法的流程图。如图1所示，该方法包括：

步骤1：首先以圆形发射回线中心点对应的地面处为原点，建立z轴朝向地下的圆柱坐标系

其中，ρ表示在上述圆柱坐标系中接收点所在的径向分量，

表示方位角，z表示高度。由于大部分地质为沉积型的分层地质，因此在TEM正演计算中将大地视为只在x-y平面内变化的层状介质。由于本发明中方法是基于阶跃响应叠加原理设计的，首先需要计算发射回线在中心点处的阶跃响应。对于层状介质以上z=-h处半径为a的圆形水平发射回线，其中，-h为发射回线的高度，由于z轴正向为指向地下，而接收线圈位于地上，故其高度为负；当接收点位于(0，0，z)处时所测得的垂直磁场分量H_z的频率域表达式如下所示：

$H_z=\frac{\mathrm{Ia}}{2}{\∫}_0^{\∞}[e^{{-u}_0(z+h)}+r_{\mathrm{TE}}{e}^{u_0(z-h)}]\frac{{\mathrm{\λ}}^2}{u_0}{J}_1\left(\mathrm{\λa}\right)\mathrm{d\λ}$

其中I为发射电流大小，J₁(λa)为第一类贝塞尔函数，下标表示贝塞尔函数的阶数，此处为1。a表示在上述圆柱坐标系中圆形发射回线的半径，λ为上式中的积分量。k_i为各层介质中电磁波的波数，下标i代表层状介质的层号。空气所在的半空间设定为第0层，其它层依次类推。k_i的表达式为

其中ω为傅里叶变换中的角频率，μ_i，ε_i和σ_i分别为各层中的相对磁导率、相对介电常数和电导率。r_TE为电场矢量从地表垂直入射到层状介质所对应的反射系数，可由

算出，其中Y₀=u₀/iωμ₀为空气的本征导纳，

为z=0处的地表波导纳，对于N层大地，地表波导纳可由如下递推关系给出：

$\widehat{Y_1}=Y_1\frac{\widehat{Y_2}+Y_1\tanh \left(u_1{h}_1\right)}{Y_1+\widehat{Y_2}\tanh \left(u_1{h}_1\right)}$

$\widehat{Y_1}=Y_i\frac{{\hat{Y}}_{i+1}+Y_i\tanh \left(u_i{h}_i\right)}{Y_i+{\hat{Y}}_{i+1}\tanh \left(u_i{h}_i\right)}$

$\widehat{Y_n}=Y_n$

其中，

为每层介质和下面一层介质界面上的波导纳。Y_i=u_i/iωμ_i为每层介质的本征导纳。由于最后一层的波导纳和本征导纳是一样的。而本征导纳是可以根据该层介质参数算出，因此从底层往上递推，可以得到地表波导纳

进而求得反射系数r_TE。

步骤2：将步骤1中的涡旋电场分量H_z表示为一Hankel积分。所述Hankel积分为：

$H_z={\∫}_0^{\∞}K\left(\mathrm{\λ}\right)J_1\left(\mathrm{\λa}\right)\mathrm{d\λ},K\left(\mathrm{\λ}\right)=-\frac{\mathrm{Ia}{\mathrm{\λ}}^2}{{2u}_0}[e^{-u_0(z+h)}+r_{\mathrm{TE}}{e}^{u_0(z-h)}]$

然后做如下代换：a=e^m和λ=e^-n，并在公式两端乘以e^m可得到如下用于计算线性移不变系统输出的褶积积分形式：

$e^m{H}_z\left(e^m\right)={\∫}_{-\∞}^{\∞}K\left(e^{-n}\right)\left[e^{m-n}{J}_1\left(e^{m-n}\right)\right]\mathrm{dn}$

其中，e^mH_z(e^m)是输出，K(e^-n)是输入，而e^m-nJ₁(e^m-n)是线性系统的移位脉冲响应，可视为一滤波器。经过代换转换为一滤波器形式。其中积分内的贝塞尔函数J₁(e^m-n)用滤波系数代替，该系数需提前根据有解析解的Hankel积分式运算出来以供使用。求取滤波系数所用的表达式与TEM响应的表达式越接近则精度越准确，滤波系数的求解在本领域中属于公知技术，在此不再赘述。最优的滤波器系数的自变量和因变量分别为如下两个矩阵。带入变换之后的公式即可求出垂直磁场频率域响应大小。

自变量为

[0.0129068125804799

0.0149207860690678

0.0172490191153463

0.0199405486456498

0.0230520632872256

0.0266490973363555

0.0308074110327511

0.0356145862113719

0.0411718709390677

0.0475963119875303

0.0550232200564072

0.0636090196687718

0.0735345437630571

0.0850088423716996

0.0982735856043615

0.113608153670764

0.131335521148493

0.151829059429431

0.175520400616997

0.202908528502544

0.234570288093798

0.271172535045600

0.313486180882605

0.362402429832490

0.418951549247639

0.484324568955362

0.559898366565402

0.647264667078035

0.748263567578565

0.865022293110741

1

1.15603957026802

1.33642748802547

1.54496305895134

1.78603843075007

2.06473109996649

2.38691085352428

2.75936339737628

3.18993327611618

3.68768909370502

4.26311451516882

4.92832907211913

5.69734342267199

6.58635444201507

7.61408635877997

8.80218512218761

10.1756743060733

11.7634821519804

13.5990508518309

15.7210409028036

18.1741453694431

21.0100312028795

24.2884274430946

28.0783832238011

32.4597220758638

37.5247231596010

43.3800648358516

50.1490715110368

57.9743110789593

67.0205976663047

77.4784629252608]

因变量为

[805577.048295096

-3807964.19063375

8764526.49993092

-13359277.2707393

15491675.0504720

-14875765.1589781

12540792.9801153

-9679146.37297799

7046616.73184324

-4940267.17979456

3382093.10700642

-2281474.97169558

1525215.98229803

-1014095.27761571

672053.159204177

-444506.544717767

293660.863327308

-193870.969371208

127938.337131712

-84407.7907618566

55680.3003614340

-36726.8107200606

24223.8451710691

-15976.7821280637

10537.2906951626

-6949.63076084068

4583.49186815395

-3022.87604203566

1993.70983203624

-1314.81992220079

867.248145638312

-571.845522126850

377.301951486474

-248.646503141396

164.223976387471

-108.041920993654

71.5537012062834

-46.9072518067537

31.1762761602306

-20.4555846129012

13.4207063382831

-9.15348771490349

5.62974288333525

-4.02602312401216

2.76905975626239

-1.39007401887418

1.16144402844965

-1.19350148451135

0.680913826618168

-0.118445989901162

-0.111956156520833

0.116112885945262

-0.0666582808304246

0.0303928418744887

-0.0123717975430253

0.00473897304666750

-0.00173902121992693

0.000606690894504463

-0.000194224656452397

5.22999198087292e-05

-9.20781851811029e-06]

步骤3：由于H_z是频率响应，故采用GST算法求取H_z的时间域响应，GST算法形式如下：

$h_z=[\ln \left(2\right)/t]\underset{j=1}{\overset{J}{\mathrm{\Σ}}}d(j,J)H_z[j\ln \left(2\right)/t],$

其中ln(2)代表2的自然对数值，t代表时间值。d(j，J)为GST的系数，属于已知量。J选取12时达到最佳运算效果。此时d(j，J)为如下矩阵，

[-0.0166666666666667

16.0166666666667

-1247

27554.3333333333

-263280.833333333

1324138.70000000

-3891705.53333333

7053286.33333333

-8005336.50000000

5552830.50000000

-2155507.20000000

359251.200000000]。

将此系数带入上式即可求取垂直磁场的时间域响应大小。

步骤4：若t表示时间，h_z(t)表示垂直磁场在时间域上的阶跃响应，n表示采样点，i₀(0⁺)表示电流初始值，Δi_n表示电流采样间隔，Δt表示时间采样间隔，N为总采样点数。那么在t-nΔt>0条件下，任意不规则发射波形信号在层状介质上的瞬变电磁响应可表示为：

$\tilde{h}\left(t\right)=i_0\left(0^{+}\right)h_z\left(t\right)+\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{\Δ}{i}_n{h}_z(t-\mathrm{n\Δt})$

由于越靠近仪器初始接收时间点，信号采样密度越对结果精度影响大，距离达到一定程度的信号采样值几乎对结果是毫无影响的，因此我们可以修正上面的采样方式从而得到一个更加高效的信号采样方式：即时间采样间隔中，Δt₁代表距离所求时间与最近的采样点之间的间隔，其他的靠前的采样点位置按照Δt_n+1=κΔt_n进行计算得到。由于信号是从响应时间点往前算起，因此递推系数κ必须大于1，其中采样初始间隔Δt₁与递推系数κ值的大小由计算时间和精度之间的平衡关系所决定。

步骤5：当时间采样间隔Δt_i大于10^-7时，步骤4公式中的i_nh_z(t-nΔt)采用自适应Lobbato积分

求取每个时间段内的响应，否则会因离散叠加而造成响应震荡。当间隔小于此值时，采用直接相乘i_nh_z(t-nΔt)计算此时间段的响应。求得每个时间段的响应之后，按照步骤4的公式累加求和即可得到任意不规则发射波形的瞬变电磁响应。

图2A示出了使用本发明的高指数信号对该实施例的处理结果与未用本发明的斜变处理结果对比图；图2B示出了使用本发明的高指数信号对该实施例的瞬变电磁响应与未用本发明的瞬变电磁响应对比图。如图2A和2B所示，使用本发明处理的结果非常精确地反应出了斜变信号与类指数信号的不同，而如果未使用该方法处理结果效果不理想，会造成很大的系统误差，这些误差势必会影响后期的数据处理与解释。

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

Claims (9)

1.一种适用于不规则发射波形的瞬变电磁正演方法，其特征在于，包括：

步骤1、以圆形发射回线中心点对应的地面处为原点，建立z轴朝向地下的圆柱坐标系，并建立接收点处的垂直磁场分量的频率域表达式；

步骤2、将所述垂直磁场分量的频率域表达式表示为Hankel积分形式，并将其转换为滤波器积分形式，并将转换后的滤波器积分形式内的贝塞尔函数使用预先设定的滤波系数代替，并求出垂直磁场分量的频率域响应大小；

步骤3、根据所述垂直磁场分量的频率域响应大小采用GST算法求取垂直磁场分量的时间域响应大小；

步骤4、根据所述垂直磁场分量的时间域响应大小建立不规则发射波形信号在大地层状介质上的瞬变电磁响应表达式；

步骤5、利用所述瞬变电磁响应表达式求取每个时间段内的瞬变电磁响应，并将每个时间段内的瞬变电磁响应进行累加得到不规则发射波形的瞬变电磁响应。

2.如权利要求1所述的方法，其中，垂直磁场分量的频率域表达式如下表示：

$H_z=\frac{\mathrm{Ia}}{2}{\∫}_0^{\∞}[e^{{-u}_0(z+h)}+r_{\mathrm{TE}}{e}^{u_0(z-h)}]\frac{{\mathrm{\λ}}^2}{u_0}{J}_1\left(\mathrm{\λa}\right)\mathrm{d\λ}$

其中，J₁(λa)为第一类贝塞尔函数，a表示在上述圆柱坐标系中圆形发射回线的半径，λ为上式中的积分量；I为发射电流大小，-h为圆形发射回线在所述圆柱坐标系中z轴的坐标值，

k_i为大地层状介质中各层的波数，下标i代表大地层状介质的层号，空气所在的半空间设定为第0层；k_i的表达式为

其中ω为角频率，μ_i，ε_i和σ_i分别为各层中的相对磁导率、相对介电常数和电导率，r_TE为电场矢量从地表垂直入射到层状介质所对应的反射系数。

3.如权利要求2所述的方法，其中，所述电场矢量从地表垂直入射到层状介质所对应的反射系数r_TE根据下式计算获得：

$r_{\mathrm{TE}}=(Y_0-\widehat{Y_1})/(Y_0+{\hat{Y}}_1)$

其中，Y₀=u₀/iωμ₀为空气的本征导纳，处的地表波导纳，对于N层大地，地表波导纳可由如下递推关系给出：

$\widehat{Y_1}=Y_1\frac{\widehat{Y_2}+Y_1\tanh \left(u_1{h}_1\right)}{Y_1+\widehat{Y_2}\tanh \left(u_1{h}_1\right)}$

$\widehat{Y_1}=Y_i\frac{{\hat{Y}}_{i+1}+Y_i\tanh \left(u_i{h}_i\right)}{Y_i+{\hat{Y}}_{i+1}\tanh \left(u_i{h}_i\right)}$

$\widehat{Y_n}=Y_n$

其中，

为每层介质和下面一层介质界面上的波导纳，Y_i=u_i/iωμ_i为每层介质的本征导纳，其中最后一层的波导纳和本征导纳相同，故从最后一层往上递推，可以得到地表波导纳

进而求得r_TE。

4.如权利要求2所述的方法，其中，所述Hankel积分如下所示：

$H_z={\∫}_0^{\∞}K\left(\mathrm{\λ}\right)J_1\left(\mathrm{\λa}\right)\mathrm{d\λ}$

$K\left(\mathrm{\λ}\right)=-\frac{\mathrm{Ia}{\mathrm{\λ}}^2}{{2u}_0}[e^{-u_0(z+h)}+r_{\mathrm{TE}}{e}^{u_0(z-h)}].$

5.如权利要求2所述的方法，其中，对所述Hankel积分做如下代换：a=e^m和λ=e^-n，并在其两端乘以e^m得到如下用于计算线性移不变系统输出的褶积积分形式：

$e^m{H}_z\left(e^m\right)={\∫}_{-\∞}^{\∞}K\left(e^{-n}\right)\left[e^{m-n}{J}_i\left(e^{m-n}\right)\right]\mathrm{dn}$

其中，e^mH_z(e^m)是输出，K(e^-n)是输入，而e^m-nJ₁(e^m-n)是线性移不变系统的移位脉冲响应，其中积分内的贝塞尔函数J₁(e^m-n)用滤波系数代替后计算得到垂直磁场分量的频率域响应H_z。

6.如权利要求5所述的方法，其中，垂直磁场分量的时间域响应大小如下得到：

$h_z=[\ln \left(2\right)/t]\underset{j=1}{\overset{J}{\mathrm{\Σ}}}d(j,J)H_z[j\ln \left(2\right)/t],$

其中，t代表时间值，d(j，J)为GST的系数，属于已知量。

7.如权利要求6所述的方法，其中，所述d(j，J)中J为12。

8.如权利要求1所述的方法，其中，所述瞬变电磁响应表达式如下所示：

$\tilde{h}\left(t\right)=i_0\left(0^{+}\right)h_z\left(t\right)+\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{\Δ}{i}_n{h}_z(t-\mathrm{n\Δt})$

其中，t表示时间，h_z(t)表示垂直磁场在时间域上的阶跃响应，n表示采样点，N为总采样点数，i₀(0⁺)表示电流初始值，Δi_n表示电流采样间隔，Δt表示时间采样间隔。

9.如权利要求8所述的方法，其中，步骤5中当时间采样间隔Δt_i大于10^-7时，所述瞬变电磁响应表达式中的i_nh_z(t-nΔt)采用自适应Lobbato积分

求取每个时间段内的响应，当当时间采样间隔Δt_i小于10^-7时，采用直接相乘i_nh_z(t-nΔt)计算此时间段的响应。

Images