一种无边缘效应的中心回线TEM全期真电阻率计算方法

摘要

本发明公开了一种无边缘效应的中心回线TEM全期真电阻率计算方法，从载流点微元的圆形回线电场公式出发，根据电场、磁场、感生电动势之间的关系、解决了中心点外场点解析求解困难的问题；根据电阻率勘探中的相对概念，应用大宗量Bessel函数的渐进式，解决了含Bessel函数积分的计算问题。本发明包括以下步骤：获得任意场点感生电动势解析表达式的步骤；获得任意场点单Bessel函数的感生电动势解析表达式

说明书

技术领域

本发明属于地球物理勘探领域，具体涉及一种电和电磁法勘探资料解释方法。。

背景技术

为了提高中心回线TEM（Transient Electro-Magnetic，TEM）的施工效率，在野外勘探中， 将观测点从中心点扩大到了中心1/3的区域如图1所示。但是，中心区域场的不均匀性，形 成如图2所示的边缘效应。由于采用了中心点公式定义的视电阻率，边缘效应在视电阻率- 深度剖面上，形成了与发射回线有关、与地质结构无关的韵律变化，导致对地下地质结构的 误判。为了解决这个问题，现有技术将大定源回线和中心回线理论公式进行了统一^[1-2]，获得 的仍然是视电阻率，且大定源回线的理论公式是从偶极子微元出发导出的，中心回线的场点 到源的距离与偶极子假设相差较远，误差仍然较大；对于真电阻率的计算，现有技术采用测 点归位的校正方法消除边缘效应^[3]，但是边缘效应是回线源场的本质表现，该方法不仅不能 消除边缘效应，还会引入新的误差。对地下地质结构的误判依然存在。

对比文件与参考文献

[1]李建平,李桐林,赵雪峰,梁太木.层状介质任意形状回线源瞬变电磁全区视电阻率的研究.地球物理 学进展,2007,22(6):1777-1780

[2]石显新,闫述,傅君眉,陈明生.瞬变电磁法中心回线装置资料解释方法的改进.地球物理学报,2009, 52(7):1931-1936

[3]http://www.phoenix-geophysics.com/

[4]Knight J H,Raiche A P.transient electromagnetic calculations using the Gaver-Stehfest inverse Laplace  transform method.Geophysics,1982,47(1):47-50

[5]Anderson W L.Numerical integration of related Hankel transforms of order 0 and 1 by adaptive digital  filtering.Geophysics,1979,44(7):1287-1305.

[6]Koefoed O,Ghoch D P，Polmen G J.Computation of type curves for electromagnetic depth sounding with a  horizontal transient coil by means of a digital linear filter.Geophysical Prospecting,1972,20:406-420.

[7]Verma R K,Koefoed O.A note on the linear filter method of computing electromagnetic sounding curves. Geophysical Prospecting,1973,21:70-76.

[8]陈明生,陈乐寿,王天生,白改先.用改进的广义逆矩阵方法解释大地电磁测深及电测深资料.地球物 理学报,1983,26(4):390-400.

发明内容

为了克服现有技术中消除边缘效应方法的缺陷，本发明提供一种无边缘效应的中心回线 TEM全期真电阻率计算方法，消除由边缘效应引起的对地下地质结构的误判。

要获得无边缘效应的中心回线TEM真电阻率的关键：第一，获得回线内任一点感生电动 势V(t)的解析表达式。但是，由于场分布的不均匀性，除中心点以外的感生电动势，并不能 通过绕接收线圈对电场E_θ的积分获得；第二，中心回线TEM理论公式中以不同形式出现的 双Bessel函数的数值积分计算问题。

为了解决以上技术问题，本发明所采用的技术方案如下。

一种无边缘效应的中心回线TEM全期真电阻率计算方法包括以下步骤：

步骤一，获得中心回线TEM任意场点垂直磁场解析表达式

在圆柱坐标系中，当回线中点与坐标原点重合时，大地表面上中心回线TEM电场E_θ的 频率域表达式为

$E_{\theta }(r,\omega )=-\mathrm{j\omega }{\mu }_{0}I\left(\omega \right)a\underset{0}{\overset{\infty }{\int }}{R}_n(\lambda ,\omega )J_1\left(\mathrm{\lambda a}\right)J_1\left(\mathrm{\lambda r}\right)\mathrm{d\lambda }---\left(1\right)$

式中r为地面上一点到坐标原点距离；ω＝2πf为圆频率，其中f为频率；μ₀＝4π×10^-7H/m 为非磁性大地磁导率；I为发射电流，a为发射回线半径；J₁为1阶Bessel函数，R_n层状大 地表面上的总反射系数；λ为Hankel变换的积分变量；

利用Maxwell旋度方程

$▿\times E=-\mathrm{j\omega }{\mu }_{0}H---\left(2\right)$

其中E为电场强度，H为磁场强度。由于电场仅有θ分量、且仅是r的函数，故有如下垂直 分量的磁场H_z

$H_z=\frac{1}{r}\frac{\partial }{\partial r}\left(r{E}_{\theta }\right)=\frac{I\left(\omega \right)a}{r}\underset{0}{\overset{\infty }{\int }}{R}_n(\lambda ,\omega )J_1\left(\mathrm{\lambda a}\right)[\lambda J_0\left(\mathrm{\lambda r}\right)-\frac{1}{r}{J}_1\left(\mathrm{\lambda r}\right)]\mathrm{d\lambda }---\left(3\right)$

式（3）中J₀是0阶Bessel函数。解析公式（3）中含有双Bessel函数，还需要进一步化成单 Bessel函数，才能应用现有滤波系数等算法求得积分的值。

步骤二，获得任意场点单Bessel函数的感生电动势解析表达式

对于普遍应用的a＝600m~800m的大发射回线，利用Bessel函数的渐进式

$J_1\left(x\right)\approx \sqrt{\frac{2}{\mathrm{\pi x}}}\cos (x-\frac{3\pi }{4})(x\to \infty )---\left(4\right)$

将公式（4）代入公式（3）

$h_z\left(t\right)\approx \frac{I\left(\omega \right)a}{r}\sqrt{\frac{2}{\mathrm{\pi a}}}\underset{0}{\overset{\infty }{\int }}{R}_n(\lambda ,\omega )\sqrt{\frac{1}{\lambda }}\cos (\mathrm{\lambda a}-\frac{3\pi }{4})[\lambda J_0\left(\mathrm{\lambda r}\right)-\frac{1}{r}{J}_1\left(\mathrm{\lambda r}\right)]\mathrm{d\lambda }---\left(5\right)$

对公式（5）做逆Laplace变换，得到时间域形式

实测感生电动势V(t)和h_z(t)的关系为

$V\left(t\right)=\frac{\partial }{\partial t}{h}_z\left(t\right)---\left(7\right)$

将公式（6）代入公式（7）后得任意场点单Bessel函数的感生电动势解析表达式为

步骤三，将公式（8）代入反演程序中，即获得无边缘效应的中心回线TEM的全期真电 阻率。

上述公式（8）的求值利用Laplace的微分性质、G-S算法、滤波系数算法。

所述的滤波系数算法为Koefoed等提供的47个滤波系数计算1阶Bessel函数，用Verma 提供的51个滤波系数计算0阶Bessel函数；或者利用Anderson提供的441个滤波系数计算 1阶和0阶Bessel函数。

本发明具有有益效果。本发明公开的无边缘效应中心回线TEM全期真电阻率解决方法， 从载流点微元的圆形回线电场公式出发，克服了偶极子微元的固有误差。根据电场、磁场、 感生电动势之间的关系、解决了中心点外场点解析求解困难的问题；根据电阻率勘探中的相 对概念，应用大宗量Bessel函数的渐进式，解决了含Bessel函数积分的计算问题，从根本上 消除了边缘效应；通过将本发明的任意场点单Bessel函数的感生电动势公式（8）代入反演程 序，即获得分层大地的全期真电阻率和厚度。提高了中心回线TEM探测地下地质结构的准确 度。全期真电阻率可以更好地反映地下地质结构；为回线中任意点的TEM响应研究提供了解 析形式的理论公式。

附图说明

图1是中心回线TEM实际勘探中测点的分布示意图，图中的矩形方框为发射回线，十 字标记为测点。

图2是具有边缘效应的实测数据测道图，观测时长30ms、时间道20道，后4道干扰较 大，仅取了前16道，发射回线600m×600m，中心测区200m×200m。

图3是无边缘效应和校正边缘效应的中心回线TEM电阻率-深度剖面对比图，其中(a)是加 拿大Phoenix公司V8仪器标配软件计算的、校正后仍有边缘效应影响的电阻率-深度剖面图， (b)是本发明的无边缘效应的全期中心回线TEM真电阻率-深度剖面图。

具体实施方式

下面结合附图，对本发明的具体实施方案作进一步详细说明。

以煤田水文地质勘探中获得的实测数据为例，发射回线600m×600m，观测在其中200m ×200m的中心区域内进行。图3中(a)是加拿大phoenix公司V8仪器标配校正边缘效应软件 计算的电阻率-深度剖面，可以看出采用测点归位的校正并没有消除边缘效应的影响。同样的 实测数据，用无边缘效应的全期真电阻率方法计算如下：

勘探施工时需要记录每个测点与发射回线的归属关系，及到发射回线中点的位置。考虑 到实际施工中发射回线一般为正方形，公式（8）中的发射半径由下式换算

$a=\frac{L}{\sqrt{\pi }}---\left(9\right)$

式中L为正方形发射回线的边长。公式（8）中场点r的取值范围在发射回线中部三分之一的 范围内，如图1所示。

应用Laplace的微分性质、G-S算法^[4]；滤波系数算法采用Anderson提供的441个滤波系数^[5]对上述公式（8）进行编程，然后与如图2所示的实测数据一起代入改进的广义逆矩阵^[8]反演 程序进行反演计算，得到全期真电阻率和地层厚度，用Surfer等绘图软件生成如图3中(b)所 示的电阻率-深度剖面。反演初始参数设置有多种方法，本实施例采用的是均匀半空间地电模 型，大地电阻率ρ₁用早期或晚期视电阻率公式估算，地层层数与时间道数相同为16与现有 的如图3中(a)所示的校正后仍有边缘效应影响的真电阻率-深度剖面^[3]相比，本发明消除了边 缘效应的影响。