CN103777248A - 一种适用于不规则发射回线的tem一维正演方法 - Google Patents

Abstract

本发明公开了一种适用于不规则发射回线的瞬变电磁正演方法，包括：以接收线圈中心对应的地面点处为原点，建立圆柱坐标系，并建立接收点处的涡旋电场分量的频率域表达式；将所述涡旋电场分量的频率域表达式表示为Hankel积分形式，并求出涡流电场分量的频率域响应大小；求取涡流电场分量的时间域响应大小；记录不规则发射回线的拐点坐标，求出两个拐点之间分段回线的函数表达式；将所述涡流电场分量的时间域响应在所述分段回线上进行积分得到每个分段回线上的电压值，并将各个分段回线上的电压值求和得到所述不规则发射回线的瞬态响应。本发明提出的方法适用于不规则发射回线的地面、海洋和航空瞬变电磁法一维正演算法。

Description

一种适用于不规则发射回线的TEM一维正演方法

技术领域

本发明涉及地质探测领域，尤其涉及一种适用于不规则发射回线的瞬变电磁法(TEM)一维正演方法。

背景技术

瞬变电磁法是利用向地面发射电磁波而探测地下地质构造的地电特征的一种方法。瞬变电磁法一般采用大回线向地面发射电磁波，然后利用近似磁偶极子的接收机进行接收地下介质反射回波。通过后期的处理可以对地下地电特征甚至地质结构进行解释、瞬变电磁方法可用于贱金属矿、地下水和石油等资源的勘探，是国家能源勘探中非常重要的一种非地震勘探方法。

在瞬变电磁法的实际作业中，由于勘探处地形、人为因素等原因，发射回线多呈现不规则形；这些总体特征使得发射源成为任意形状的空间展源。在传统的正反演和数据处理中：发射回线常用等面积的圆形回线近似代替。但近年来，一方面随着计算机处理能力和精密采集技术的发展，使得发射回线空间展源效应的考虑成为可能。目前国内外瞬变电磁法的相关文献在空间展源方面主要涉及了圆形、矩形和一些简单多边形响应的求取，因其计算方法都是针对特定的形状，使得其在实际作业中往往面临很大的不足。

发明内容

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

为此，本发明提出了一种适用于不规则发射回线的瞬变电磁正演方法，其包括：

步骤1：以接收线圈中心对应的地面点处为原点，建立z轴朝向地下的圆柱坐标系

并建立接收点处的涡旋电场分量的频率域表达式；

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

步骤3：采用GST算法求取涡流电场分量的时间域响应大小；

步骤4：记录不规则发射回线的拐点坐标，并对整个发射回线进行分段三次样条插值，求出两个拐点之间分段回线的函数表达式；

步骤5：将所述涡流电场分量的时间域响应在所述分段回线上进行积分得到每个分段回线上的电压值，并将各个分段回线上的电压值求和得到所述不规则发射回线的瞬态响应。

本发明从发射回线的空间坐标出发。首先借助有源空间的互易定理求得层状大地上磁偶极子响应，进而将可视作磁偶极子的接收线圈与任意发射回线源进行数学运算上的收发互换，然后利用基于矢量运算的第二类曲线积分求取了任意形状回线的响应。因为该方法是针对发射源坐标构建的，使用者无需关心发射源的实际形状，只需用使用者确定回线坐标即可。这使得该算法在实际作业中具有很大的优越性。仿真对比也显示了所提出方法的正确性。该方法涉及地球物理电磁勘探领域中的各种常见平台瞬变电磁法方法与技术，包括，一维层状模型磁偶极子响应算法；收发互换算法和矢量积分算法。本发明方法可以有效解决一维层状模型上任意几何形状发射回线的瞬变电磁响应问题。本发明将可用于瞬变电磁法中的高精度系统预研、正反演计算和数据解译中。

附图说明

图1是本发明中适用于不规则发射回线TEM一维正演方法流程图；

图2A是一具有不规则形状的发射回线的示意图；

图2B是使用本发明的方案和传统等效圆形回线方案对图2A示出的不规则发射回线中的A点进行了处理后得到的瞬态响应对比图；

图2C是使用本发明的方案和传统的方案对图2A示出的不规则发射回线中的B点进行了处理后的结果对比图。

具体实施方式

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

在瞬变电磁的实际作业中，由于勘探处地形、人为因素等原因，发射回线多呈现不规则形，假设一种比较常见的情形：实际作业在铺设发射回线的过程中由于河流、树林或山丘等原因将方形回线铺设成了如图2所示的不规则形状。由于传统的处理方法一般用等面积的圆形回线代替该不规则回线计算其响应。

本发明直接从时间域出发，以磁偶极子的响应为基响应，通过矢量积分算法求取了任意形状回线的响应。为了提高运算效率，本文利用高等电磁场理论中的Carson互易定理，提出了一种基于坐标的矢量积分运算方法。仿真分析显示本文方法很好地解决了运算时间和结果精度之间的矛盾。

图1示出了本发明提出的一种适用于不规则发射回线TEM一维正演方法流程图。如图1所示，其包括以下步骤：

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

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

表示方位角，z表示高度。因为实际应用中接收线圈足够小且距离发射线圈足够远，一般将水平放置的接收线圈作为垂直磁偶极子进行分析计算。此外由于大部分地质为沉积型的分层地质。因此在TEM正演计算中将大地视为只在x-y平面内变化的层状介质。由于本发明中方法是基于电磁场Carson互易定理设计的，首先将接收线圈作为发射源，不规则回线作为接收回线处理。对于层状介质以上z=-h处磁偶极矩为m的垂直磁偶极子，其中，-h为垂直磁偶极子即接收线圈的高度，由于z轴正向为指向地下，而接收线圈位于地上，故其高度为负；当接收点位于

处时所测得的涡旋电场分量

的频率域表达式如下所示：

其中J₁(λρ)为第一类贝塞尔函数。下标表示贝塞尔函数的阶数，此处为1。ρ表示在上述圆柱坐标系中接收点所在的径向分量，λ为上式中的积分量。

k_i为各层介质中的波数，下标i代表层状介质的层号。空气所在的半空间设定为第0层，其它层依次类推。k_i的表达式为

其中ω为角频率，μ_i，ε_i和σ_i分别为各层中的相对磁导率、相对介电常数和电导率。

为各层介质的阻抗率。i为复数中的虚数单位。r_TE为电场矢量从地表垂直入射到层状介质所对应的反射系数，可由

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

${\hat{Y}}_1=Y_1\frac{{\hat{Y}}_2+Y_1\tanh \left(u_1{h}_1\right)}{Y_1+{\hat{Y}}_2\tanh \left(u_1{h}_1\right)}$

${\hat{Y}}_i=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)}$

${\hat{Y}}_n=Y_n$

其中，

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

进而求得反射系数r_TE。

步骤2：步骤1中的涡旋电场分量

表示为一Hankel积分。将其简写为

$K\left(\mathrm{\λ}\right)=-\frac{{\hat{z}}_{0}m{\mathrm{\λ}}^2}{4\mathrm{\π}{u}_0}[e^{-u_0(z+h)}+r_{\mathrm{TE}}{e}^{u_0(z-h)}].$ 做代换ρ=e^m和λ=e^-n，并在公式两端乘以e^m可得到如下用于计算线性移不变系统输出的褶积积分形式 $e^{m}E\left(e^m\right)={\∫}_{-\∞}^{\∞}K\left(e^{-n}\right)\left[e^{m-n}{J}_1\left(e^{m-n}\right)\right]\mathrm{dn}.$ 其中e^mE(e^m)是输出，K(e^-n)是输入，而e^m-nJ_n(e^m-n)是线性系统的移位脉冲响应，可视为一滤波器。经过代换转换为一滤波器形式。其中积分内的贝塞尔函数J₁(e^m-n)用如下滤波系数代替。该滤波器系数的自变量和因变量分别为如下两个矩阵。带入变换之后的公式即可求出涡流电场的频率域响应大小。

自变量为

[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：由于

是频率响应，故采用GST算法求取的时间域响应，GST算法形式如下：

其中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]。

将此系数带入上式即可求取涡流电场的时间域响应大小。

求得垂直磁偶极子的时间域涡流电场值之后，根据高等电磁场理论中的Carson互易定理可知，各向同性的线性介质中将发射回线和接收线圈功能互换，只要通入接收线圈中的电流和之前发射线圈中的电流值是相等的，那么发射回线中所测得的电压值与互换前所测电压值相等。因此将此涡流电场值沿发射回线进行积分所得到的电压值，就是不规则回线作为发射时在接收线圈中接收到的电压值。

步骤4：记录不规则发射回线的拐点坐标(x_i，y_i)，并对整个回线进行分段三次样条插值，用于求出两个拐点之间直线的表达式y(x)，为后面步骤在每段上的积分做准备。

步骤5：不规则回线每段上的电压可以由矢量积分

求取得出，其中v_i为每段回线的电压值，C_i代表每段回线的积分路径。然后利用公式

将此矢量积分转换为普通的一次积分形式进行积分运算。其中y′(x)为y(x)的导数。cosα和cosβ分别为

与所述圆柱坐标系对应的直角坐标系中x和y坐标轴的夹角的余弦值，表达式分别为 $\cos \mathrm{\α}=\frac{y\left(x\right)}{\sqrt{{\left(x\right)}^2+{\left[y\left(x\right)\right]}^2}}$ 和 $\cos \mathrm{\β}=\frac{x}{\sqrt{{\left(x\right)}^2+{[y\left(x\right)}^2}}.$ 求得每段回线的电压后，求代数和

即可得到任意不规则发射回线的瞬态响应。

图2A示出了一具有不规则形状的发射回线示意图；图2B示出了使用本发明的方案和传统等效圆形回线方案对图2A示出的不规则发射回线中的A点进行了处理后得到的瞬态响应对比图；图2C示出了使用本发明的方案和传统的方案对图2A示出的不规则发射回线中的B点进行了处理后的结果对比图。如图2A所示，其中发射回线由于需要绕开树林和房屋等因素造成发射回线呈现一种不规则五边形的形状。设定发射回线下方的地层参数为σ₁=0.1S／m，σ₂=0.001S／m和h₁=100m的两层大地上。接收点分别设定在中心A点和大收发距B点；如2B所示，其中使用本发明计算得到的瞬变电磁响应在前期低于等效圆形回线的响应，但在后期高于等效圆形回线的响应；如图2C所示，其中使用本发明就算得到的瞬变电磁响应全部高于等效圆形回线的响应，并且该误差不随时间变小。对比图2B和2C可见使用本发明处理的结果非常理想，这种精确的结果在大收发距的情况下将会更加明显。而未使用该方法处理结果效果并不理想，会造成很大的系统误差。这种误差将会影响后期的反演解释，甚至造成假异常。

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

Claims (8)

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

步骤1：以接收线圈中心对应的地面点处为原点，建立z轴朝向地下的圆柱坐标系

并建立接收点处的涡旋电场分量的频率域表达式；

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

步骤3：采用GST算法求取涡流电场分量的时间域响应大小；

步骤4：记录不规则发射回线的拐点坐标，并对整个发射回线进行分段三次样条插值，求出两个拐点之间分段回线的函数表达式；

步骤5：将所述涡流电场分量的时间域响应在所述分段回线上进行积分得到每个分段回线上的电压值，并将各个分段回线上的电压值求和得到所述不规则发射回线的瞬态响应。

2.如权利要求1所述的方法，其特征在于，所述涡旋电场分量的频率域表达式如下所示：

其中，J₁(λρ)为第一类贝塞尔函数，ρ表示在上述圆柱坐标系中接收点所在的径向分量，λ为上式中的积分量；m为垂直磁偶极子的磁偶极矩，其z轴坐标为垂直磁偶极子的高度-h，k_i为大地层状介质中各层的波数，下标i代表大地层状介质的层号，空气所在的半空间设定为第0层；k_i的表达式为其中ω为角频率，μ_i，ε_i和σ_i分别为各层中的相对磁导率、相对介电常数和电导率，

为各层介质的阻抗率；r_TE为电场矢量从地表垂直入射到层状介质所对应的反射系数。

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

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

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

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

${\hat{Y}}_1=Y_1\frac{{\hat{Y}}_2+Y_1\tanh \left(u_1{h}_1\right)}{Y_1+{\hat{Y}}_2\tanh \left(u_1{h}_1\right)}$

${\hat{Y}}_i=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)}$

${\hat{Y}}_n=Y_n$

其中，

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

进而求得r_TE。

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

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

5.如权利要求4所述的方法，其中，对所述Hankel积分中进行如下代换：ρ=e^m和λ=e^-n，并在其两端乘以e^m可得到如下滤波器形式：

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

6.如权利要求5所述的方法，其中，涡流电场分量的时间域响应大小如下获得：

其中，ln(2)代表2的自然对数值，t代表时间值，d(j，J)为GST的系数。

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

8.如权利要求1所述的方法，其中，所述涡流电场分量的时间域响应在各个分段回线上的电压值如下获得：

其中，v_i为每段回线的电压值，y′(x)为每个分段回线函数表达式y(x)的导数；cosα和eosβ分别为涡流电场分量的时间域响应与x和y坐标轴的夹角的余弦值；x_i和x_i+1为构成第i个分段回线的发射回线拐点的横坐标。

Images