CN101710187A - 一种时间域航空电磁高度校正方法 - Google Patents

Abstract

本发明公开了一种时间域航空电磁高度校正方法，该方法利用x、z两个分量的感生电动势数据联合计算一种等效均匀半空间电阻率序列，再按等效均匀半空间电阻率为均匀半空间模型的电阻率参数值，在校正高度下正演计算均匀半空间模型的各瞬时感生电动势响应数据，即为最终的校正结果。其中，利用等效均匀半空间电阻率值实现观测感生电动势数据校正是本发明的创新部分。本发明计算速度快，校正结果的相对误差较小，符合航空电磁资料处理要求。

Description

一种时间域航空电磁高度校正方法

技术领域

本发明涉及一种时间域航空电磁高度校正方法，属于时间域航空电磁资料处理方法范畴。

背景技术

时间域航空瞬变电磁方法是一种快速、高效的地球物理勘探方法，在地质填图、矿产资源调查、地下水勘查等等方面发挥着重要作用。与地面方法相比，航空电磁法中，飞机飞行速度快，工作效率高，但另一方面，空中测量方法中飞机震动又易导致飞行高度的波动，而高度计数数据误差都会导致记录的高度数据与真实高度有较大的偏差。高度数据的偏差会导致剖面上出现假异常，而根据错误的高度数据进行反演解释，会得到错误的结果。所以需要对航空瞬变电磁资料进行高度校正处理，将不同飞行高度的资料校正为同一飞行高度的资料，或对存在的高度误差数据进行改正，以降低剖面解释的难度，方便反演成像的进行，减少由于飞行高度波动而造成的假异常。

与航磁、航放资料等其它航空地球物理资料的高度校正相比，航空电磁高度校正绝非一件容易的事，这是因为，校正量与飞机下方地下介质的电阻率构造之间有着复杂的函数关系，到目前为止，还没有一个较好的方法来完成航空电磁高度校正。人们只好寻找近似的校正方法，Green A.(1998)则提出基于偶极源的镜像深度的视深度校正方法，他使用了Maxwell的薄片镜像扩散近似的方法，将地下介质近似为一电导S的薄片，由各采样时刻t对应的B场数据计算镜像源的深度值，从而获得在一定的校正高度时对应的航空电磁响应，使得资料得以校正。Wolfgram(1993)也以摘要的形式提供这种校正思路，但没有公布详细的处理方法。Green认为均匀半空间的结果是误差最大的，所以仅提供了均匀半空间的理论校正结果，事实上，虽然用薄片模型来近似层状介质比均匀半空间更合适，但前提是需要得到较准确的电导值S和其深度值d，否则，得到的结果误差会更大，所以Green认为层状介质模型响应的高度校正误差是最大的这一说法是有条件的。如果校正模型用电阻率来描述的话，均匀半空间模型校正方法却是最合适的校正模型，而且仅有电阻率一个参数，而薄片电导模型需要获得薄片电导和薄片深度两个参数，准确地得到这两个值又困难得多。所以，本发明使用电阻率模型，根据响应数据确定地下介质等效均匀半空间模型参数，从而按等效均匀半空间模型的响应来做高度校正。

发明内容

本发明的目的在于提供一种校正速度快、校正结果较准确的时间域航空电磁高度校正方法，该方法用各校正时刻的等效均匀半空间模型来做高度校正，计算简单，不需要反演迭代过程，如果地下介质是均匀半空间的特殊情况时，得到的等效均匀半空间电阻率值为真值，校正的误差就是0。

为达到上述目的，本发明采用的方法包括以下步骤：

(1)对实测的两感生电动势分量数据进行转换

在实际飞行高度下，由两个感生电动势分量实测数据Vx、Vz计算对应各校正时刻T_k＝(t_k+t_k+nk)/2的幅值比值数据

$A\left(T_k\right)=\frac{V_x^2\left(t_k\right)+V_x^2\left(t_{k+\mathrm{nk}}\right)}{V_z^2\left(t_k\right)+V_z^2\left(t_{k+\mathrm{nk}}\right)},$ 所述k＝1，2，…，n-nk，n为采样时间个数；

(2)计算幅度比值数据查询表

给定m个离散化查询电阻率数据ρ_i ⁰(i＝1，2，...，m)，分别作为等效均匀半空间的电阻率值，先计算它的两分量航空瞬变电磁响应V_x ⁰、V_z ⁰，再计算理论幅度比值查询数据体

$A^0({\mathrm{\ρ}}_i^0,T_k)=\frac{{\left[V_x^0({\mathrm{\ρ}}_i^0,t_k)\right]}^2+{\left[V_x^0({\mathrm{\ρ}}_i^0,t_{k+\mathrm{nk}})\right]}^2}{{\left[V_z^0({\mathrm{\ρ}}_i^0,t_k)\right]}^2+{\left[V_z^0({\mathrm{\ρ}}_i^0,t_{k+\mathrm{nk}})\right]}^2},$ 并将其作为计算等效均匀半空间电阻率的查询数据表保存在一个二维数组中，所述i＝1，2，…，m，k＝1，2，…，n-nk，n为采样时间个数；

(3)计算等效均匀半空间电阻率值

以各校正时刻T_k对应的两分量幅度比值数据A⁰(ρ_i，T_k)作为电阻率的单值函数关系，用局部样条插值的方法计算实测两分量数据的幅度比值数据A(T_k)对应的电阻率值，即为所求的第k个校正时刻的等效均匀半空间电阻率值，逐一计算n-nk个校正时刻的等效均匀半空间电阻率值，即完成了等效均匀半空间电阻率的计算，所述i＝1，2，…，m；k＝1，2，…，n-nk，n为采样时间个数；

(4)计算各校正时刻的感生电动势数据

根据各校正时刻的等效均匀半空间电阻率值，在给定的校正飞行高h1下，逐一计算全部n-nk个校正时刻对应的各分量感生电动势数据，即为最终的校正结果。

本发明利用均匀半空间x、z分量的感生电动势响应数据定义一种幅度比值，并以该幅度比值作为均匀半空间电阻率的单值函数关系，将观测数据作为变量，通过插值方法计算对应的均匀半空间电阻率，即得到瞬时等效均匀半空间电阻率，再经正演计算新的校正高度下的感生电动势响应数据，即得到最终的校正结果。

本发明的具体原理如下：

如图1所示均匀大地上空的航空电磁系统，发射垂直磁偶极距为M，距地面高度为h0，接收线圈高度为hr，水平收发距为r，接收线圈有效面积为Sr，n个采样时间分别为t_k，k＝1，2，…，n。

假定离散化的均匀半空间电阻率为ρ_i ⁰，i＝1，2，...，m，则各电阻率值对应的均匀半空间模型的x、z分量航空瞬变感生电动势下阶跃脉冲响应为：

$V_x({\mathrm{\ρ}}_i^0,t_k)=-{\mathrm{\μ}}_0{S}_r\frac{M}{4\mathrm{\π}}{L}^{-1}\left\{{\∫}_0^{\∞}\right[e^{-\mathrm{\λ}(\mathrm{ht}-\mathrm{hr})}-\frac{\mathrm{\λ}-u}{\mathrm{\λ}+u}{e}^{-\mathrm{\λ}(h0+\mathrm{hr})}\left]{\mathrm{\λ}}^2{J}_1\left(\mathrm{\λr}\right)\mathrm{d\λ}\right\}---\left(1\right)$

$V_z({\mathrm{\ρ}}_i^0,t_k)=-{\mathrm{\μ}}_0{S}_r\frac{M}{4\mathrm{\π}}{L}^{-1}\left\{{\∫}_0^{\∞}\right[e^{-\mathrm{\λ}(\mathrm{ht}-\mathrm{hr})}+\frac{\mathrm{\λ}-u}{\mathrm{\λ}+u}{e}^{-\mathrm{\λ}(h0+\mathrm{hr})}\left]{\mathrm{\λ}}^2{J}_0\left(\mathrm{\λr}\right)\mathrm{d\λ}\right\}---\left(2\right)$

其中，μ₀为真空磁导率， $u=\sqrt{{\mathrm{\λ}}^2-{\mathrm{\κ}}^2},$ κ为波数，满足 ${\mathrm{\κ}}^2=-s{\mathrm{\μ}}_0/{\mathrm{\ρ}}_i^0,$ L^-1为逆拉普拉斯变换，J₀和J₁分别为第一类0阶和1阶贝塞尔函数。计算过程中，采用数值滤波算法(Guptasarm，1997)计算其中的Hankel变换；其中的逆拉普拉斯变换，借助Gaver-Stehfest概率变换算法(简称G-S变换法，Knight，1982)进行计算。

给定航空电磁的x、z两分量感生电动势响应时间序列ε_x、ε_z后，定义时刻T_k＝(t_k+t_k+nk)/2，(k＝1，2，...，n-nk)时的幅度比值函数，它是电阻率和T_k的函数：

$F_A(\mathrm{\ρ},T_k)=\frac{{\mathrm{\ϵ}}_x^2(\mathrm{\ρ},t_k)+{\mathrm{\ϵ}}_x^2(\mathrm{\ρ},t_{k+\mathrm{nk}})}{{\mathrm{\ϵ}}_z^2(\mathrm{\ρ},t_k)+{\mathrm{\ϵ}}_z^2(\mathrm{\ρ},t_{k+\mathrm{nk}})}---\left(3\right)$

按上式的定义，对两分量感生电动势观测数据Vx、Vz做变换，计算各校正时刻T_k时的幅度比值数据：

$A\left(T_k\right)=\frac{V_x^2\left(t_k\right)+V_x^2\left(t_{k+\mathrm{nk}}\right)}{V_z^2\left(t_k\right)+V_z^2\left(t_{k+\mathrm{nk}}\right)},(k=\mathrm{1,2},\·\·\·,n-\mathrm{nk})---\left(4\right)$

将A作为计算等效均匀半空间电阻率时的观测数据。由(1)、(2)和(4)知，比值计算导致本发明得到的等效均匀半空间电阻率值与发射磁矩、接收线圈面积无关，该方法的优势明显。

同样，计算幅度比理论查询数据，用于等效均匀半空间电阻率数据的计算。它由不同电阻率的均匀半空间的航空瞬变x、z分量响应数据按(3)式计算得到。其中的电阻率值选择在一般大地电阻率值范围内，并采用对数等间隔地离散为m个值，依次计算这m个电阻率对应的两分量感生电动势V_x ⁰、V_z ⁰，并计算n-nk个校正时刻的理论幅度比值查询数据：

$A^0({\mathrm{\ρ}}_i^0,T_k)=\frac{{\left[V_x^0({\mathrm{\ρ}}_i^0,t_k)\right]}^2+{\left[V_x^0({\mathrm{\ρ}}_i^0,t_{k+\mathrm{nk}})\right]}^2}{{\left[V_z^0({\mathrm{\ρ}}_i^0,t_k)\right]}^2+{\left[V_z^0({\mathrm{\ρ}}_i^0,t_{k+\mathrm{nk}})\right]}^2},(i=\mathrm{1,2},\·\·\·,m;k=\mathrm{1,2},\·\·\·,n-\mathrm{nk})---\left(5\right)$

将A⁰保存在二维数组中，可用它来求整个工区或区域的资料的等效均匀半空间电阻率值，从而节省大量计算时间。

由于单个感生电动势分量响应是均匀半空间电阻率的双值函数关系，所以由单个感生电动势数据得到的等效均匀半空间电阻率值有两个，无法得到准确结果，为克服这一难题，本发明提出的幅度比值定义的函数与电阻率值之间呈单值函数关系，这样可以得到唯一的等效均匀半空间电阻率值，从而可将其作为地下介质的等效模型进行飞行高度校正。如下分别以单个感生电动势数据V_x ⁰、V_z ⁰和用两分量联合数据A⁰，研究它们与电阻率值的函数关系。假定地下介质为三层介质，电阻率分别为100欧姆米、50欧姆米和100欧姆米，上两层厚度均分别为40米和20米，飞机飞行高度为120米，发射下阶跃脉冲信号，发射磁矩为30万Am²，接收线圈有效面积为100平方米，距地高度为70米，收发线圈水平距为100米，延时0.01ms开始采样，共28个采样时间，计算等效均匀半空间电阻率时采用的时间间隔数目nk取为3。图2为将各校正时刻的感生电动势数据与均匀半空间电阻率值之间的函数曲线，其中图(a)、图(b)和图(c)分别为各校正时刻的V_x ⁰、V_z ⁰和A⁰与电阻率的函数关系曲线集，可见，单个感生电动势数据V_x ⁰、V_z ⁰与电阻率之间是双值函数关系，而A⁰与电阻率之间是单值函数关系。所以，用两分量感生电动势幅度比值数据进行视电阻率查询得到的解是单值的，从而可根据观测数据单值地得到确定的一个等效均匀半空间电阻率值。

根据上述思路，由实测两分量感生电动势数据得到幅度比值观测数据A和理论查询表数据A⁰，对各校正时刻对应的数据，依次采用查询和局部样条插值方法得到对应时刻的等效均匀半空间电阻率值。先查询A⁰数据表对应的各电阻率值对应的两分量幅度比值理论值，与数据A中的对应校正时刻的幅度比值观测值比较，如果两者相等，就将对应的电阻率值作为该校正时刻的等效均匀半空间电阻率值，如果没有相等的，则搜索最接近观测数据的理论查询数据所对应的电阻率值ρ_p ⁰，用ρ_p ⁰及其附近若干个查询电阻率值序列和它们对应的理论幅度比值查询数据序列构成离散函数对应关系，进行局部样条搜索，得到幅度比值观测数值对应的等效均匀半空间电阻率值ρ_p ^s。用计算公式表示如下：

其中，ρ_p ^s是根据查询电阻率值与校正时刻T_i时的理论查询数据之间的局部离散函数关系，用样条插值得到观测数据对应的等效均匀半空间电阻率值，具体地说，就是用包含ρ_p0 ⁰的nn＝k2-k1+1个查询电阻率局部序列(ρ_k1 ⁰，ρ_k1+1 ⁰，…，ρ_p0 ⁰，…，ρ_k2 ⁰)，把它看作相对应的两分量比值理论查询值序列(A_k1 ⁰，A_k1+1 ⁰，…，A_k2 ⁰)的离散函数，再用样条插值得到的电阻率值ρ_p ^s，这里， ${\mathrm{\ρ}}_{p0}^0=\underset{l}{\min }\left(\right|A\left(T_i\right)-A^0({\mathrm{\ρ}}_l^0,T_i)\left|\right),$

Np为给定的参数，表示计算局部样条插值的电阻率值序列位于查询电阻率ρ_p0 ⁰值前后的个数，它不能选的过大，否则会导致计算量增加，对计算结果也无意义，在本发明使用的计算例子中均取为11。

最后，对各个校正时刻，用上述计算的等效均匀半空间电阻率值为等效地下模型电阻率值，在新的高度(即校正高度)下，按(1)和(2)计算均匀半空间的各分量感生电动势瞬时响应，即为对应校正时刻下的校正结果。全部校正时刻计算完毕后，即完成高度校正处理。

作为特例，如果被观测的实际地下介质为均匀半空间，则观测数据是由均匀半空间模型的响应数据转换得到，由上述的查询方法中幅度比值数据与电阻率的单值函数关系，反算得到的等效均匀半空间电阻率必然是真实的均匀半空间模型的电阻率值。显然，按得到的等效均匀半空间电阻率值正演计算校正高度下的响应数据必然是真实的，没有误差，所以对均匀半空间模型的响应数据而言，按本发明方法，无论校正高度值是多少，校正结果是精确的，校正误差均为0。

图3为按上述方法由上述的三层模型的响应数据计算得到的均匀等效半空间电阻率曲线。由图可见，早期反映第一层介质电阻率100欧姆米，随后，受第二层电阻率为50欧姆米的低阻层的影响，均匀等效半空间电阻率变小，达到最低值82.95欧姆米后，又受100欧姆米的基底层的影响，均匀等效半空间电阻率值又渐渐增高，在校正时刻7.3208ms时对应地下介质的等效均匀半空间电阻率值达到103.11欧姆米，最终维持在第3层介质电阻率值100欧姆米左右，符合真实情况。图4为校正高度分别为60m、90m、150m和180m时校正后的感生电动势分量的相对误差曲线对比图，其中相对误差的计算公式为(校正后的感生电动势-理论模型在校正高度时的感生电动势)/理论模型在校正高度时的感生电动势*100％。可见，飞行高度误差分别为-60m、-30m、30m和60m时，校正结果与真值的相对误差均在-10％～12％之内。实际航空瞬变电磁飞行高度波动一般不会超过30m，误差均在早期和中晚期部分，晚期后，误差逐渐变小。这是固定翼航空电磁常用的系统几何参数的校正模拟实验，结果表明校正误差均很小，可以满足一般情况下的处理要求。

现按Green A.设计的低拉与高拉系统两种极端情况下的几何参数做高度校正模拟实验，比较几种高度校正的校正效果。图5是低拉(-90，-70)系统(即接收线圈落后发射线圈90米，且在其下70米)分别从90米、180米和240米校正到中等高度120米的x、z分量感生电动势数据的校正量相对误差对比图；图6是高拉(-145，-30)系统分别从90米、180米和240米校正到中等高度120米的x、z分量感生电动势数据的校正量相对误差对比图。可以看出，无论是低拉系统还是高拉系统，两个感生电动势分量的相对误差曲线变化规律类似，飞行高度为90米向高处校正时，相对误差较大，变化范围在-20％～12％之间；反之，由高的飞行高度向低飞行高度校正时，相对误差要小一些，变化范围在-6％～12％之间。实际工作可利用这些特征选择工作参数，利用本发明的校正方法可以达到较好的处理效果。

到目前为止，仅有Green A.(1998)公开发表过时间域航空电磁高度校正的文章。Green A.使用的是B场数据，但是绝大多数瞬变电磁仪器是测量感生电动势数据的，测量B场对仪器的要求苛刻。与B场数据是电阻率的单值函数的优点不同，单个分量感生电动势是电阻率的双值函数，不能由传统的单个分量感生电动势数据反算得到唯一的等效均匀半空间电阻率值。所以本发明显然可应用B场数据进行高度校正，更重要的是提出应用两分量感生电动势数据幅度比值来进行电阻率值反算，可得到唯一的等效均匀半空间电阻率值，从而能够用两分量感生电动势资料完成高度校正。另外，Green A.的薄片近似方法高度校正仅给出简单的均匀半空间模型的校正结果，本发明则给出多层模型的校正结果，对均匀半空间模型，根据我们的算法原理，计算的等效均匀半空间模型的电阻率值就是真值，因而校正结果是完全正确的，相对误差为0。对一般常用的固定翼航空电磁系统几何参数条件下，本发明给出的校正误差均在10％～12％以内，结果稳定、可靠，符合航空电磁要求。再次，Green A.的薄板近似高度校正方法，需要计算薄板电导S和深度d，需要进行反演迭代(Tartaras，2000；Zhdanov，2002)或CDI(Guimin Liu，1993)的过程，计算时间长，而本发明采用数据查询的方法和对校正高度的一次均匀半空间响应的正演计算过程，速度很快，符合航空电磁的快速处理要求。

本发明的核心在于：(1)使用x、z两个分量的感生电动势计算等效均匀半空间电阻率；(2)用间隔为nk的时间道的观测数据来计算等效均匀半空间电阻率；(3)按两个分量的两个时间道响应数据的幅度之比定义幅度比值，依该定义分别由观测数据计算幅度比观测数据和由一系列电阻率查询数据对应的两分量幅度比值理论查询数据；(4)根据各校正时刻的理论查询数据与电阻率之间的单值函数关系，用局部样条插值法计算各校正时刻对应的等效均匀半空间电阻率值；(5)由等效均匀半空间电阻率瞬时值计算新的校正高度下均匀半空间模型的各分量感生电动势，即为最终的高度校正结果。基于上述特点，本发明具有如下优点：

(1)由于充分利用了两个分量的航空电磁测量信息，因此比任何一个单独感生电动势分量获得的地下介质信息更充分，从而使得最终计算得到的地下介质的等效电阻率值更接近最佳值；

(2)由于采用了以间隔为nk的两个时间道的感生电动势幅度比值数据查询等效均匀半空间电阻率的方法，因此相对于直接用一个时刻点的感生电动势数据计算视电阻率的方法，其抗噪性能高，计算过程也更稳定；

(3)可直接使用感生电动势数据进行校正，不必进行代价更高的测量磁场方法或将感生电动势转化为磁场的方法进行校正。本发明采用联合两分量感生电动势数据的方法定义的幅度比值函数是电阻率的单值函数，可以克服单分量感生电动势为电阻率的双值函数的缺点，使计算等效均匀半空间电阻率变得更为容易，确保了校正方法的可行性和校正结果的正确性；

(4)计算过程仅有简单的查询、插值和一次正演计算，计算速度快，符合航空瞬变电磁要求，使得航空电磁在飞行中的实时处理成为可能。

附图说明

图1为均匀大地上空的航空电磁系统示意图。

图2为各校正时刻下Vx、Vz和幅度比值A与查询电阻率的关系曲线。其中图(a)表示x方向感生电动势分量与查询电阻率的关系曲线；图(b)表示z方向感生电动势分量与查询电阻率的关系曲线；图(c)表示由x、z方向感生电动势分量联合计算的幅度比值与查询电阻率的关系曲线。

图3为由电阻率为100-50-100欧姆米、前两层厚度分别为40米和20米的三层模型的理论响应数据计算的等效均匀半空间电阻率曲线。

图4为将飞行高度为130m时，三层模型(100-50-100欧姆米、前两层厚度分别为40米和20米)的响应校正到60m、90m、150m和180m高度时的相对误差曲线对比图。其中图(a)为x分量感生电动势数据校正结果，图(b)为z分量感生电动势数据校正结果。

图5为低拉系统(-90，-70)系统在三层模型(100-50-100欧姆米、前两层厚度分别为40米和20米)上、飞行高度分别为90m、180m和240m的响应校正到120m时的相对误差曲线对比图。其中图(a)为x分量感生电动势数据校正结果，图(b)z分量感生电动势数据校正结果。

图6为高拉系统(-145，-30)系统在三层模型(100-50-100欧姆米、前两层厚度分别为40米和20米)上、飞行高度分别为90m、180m和240m的响应校正到120m时的相对误差曲线对比图。其中图(a)为x分量感生电动势数据校正结果，图(b)为z分量感生电动势数据校正结果。

具体实施方式

本发明是基于两分量感生电动势数据的一种幅度比值与电阻率的单值函数关系，由观测数据查询、局部样条插值的方法得到各校正时刻的航空瞬变电磁资料的一个等效均匀半空间电阻率序列，再按瞬时等效均匀半空间电阻率值计算校正高度下的均匀半空间响应，从而得到新的飞行高度的校正结果。具体包括如下步骤：

(1)输入计算参数：查询电阻率个数m、对数等间隔的查询电阻率序列ρ_i ⁰，i＝1，2，...，m，局部样条插值参数Np；航空瞬变电磁系统参数：飞行高度h0、校正高度h1、发射磁矩M、接收线圈有效面积Sr、接收线圈的高度hr、收发线圈的水平距r，接收信号的延时采样个数n及采样时间序列t_i(i＝1，2，...，n)；计算幅度比值用的两个采样数据的时间间隔个数nk；接收x、z分量感生电动势信号数据Vx、Vz；

(2)计算幅度比值观测数据：对n-nk个校正时刻T_k＝(t_k+t_k+nk)/2进行循环，计算各校正时刻时的幅度比值观测数据： $A\left(T_k\right)=\frac{V_x^2\left(t_k\right)+V_x^2\left(t_{k+\mathrm{nk}}\right)}{V_z^2\left(t_k\right)+V_z^2\left(t_{k+\mathrm{nk}}\right)},(k=\mathrm{1,2},\·\·\·,n-\mathrm{nk});$

(3)计算幅度比值理论查询数据表，计算步骤如下：

(a)对查询电阻率序列ρ_i ⁰，(i＝1，2，...，m)进行循环计算，令i＝1；

(b)计算电阻率为ρ_i ⁰的均匀半空间模型的理论航空瞬变电磁两分量感生电动势响应序列

${\hat{V}}_x^0({\mathrm{\ρ}}_i^0,t_k)=-{\mathrm{\μ}}_0{S}_r\frac{M}{4\mathrm{\π}}{L}^{-1}\left\{{\∫}_0^{\∞}\right[e^{-\mathrm{\λ}(\mathrm{ht}-\mathrm{hr})}-\frac{\mathrm{\λ}-u}{\mathrm{\λ}+u}{e}^{-\mathrm{\λ}(\mathrm{ht}+\mathrm{hr})}\left]{\mathrm{\λ}}^2{J}_1\left(\mathrm{\λr}\right)\mathrm{d\λ}\right\}$

${\hat{V}}_z^0({\mathrm{\ρ}}_i^0,t_k)=-{\mathrm{\μ}}_0{S}_r\frac{M}{4\mathrm{\π}}{L}^{-1}\left\{{\∫}_0^{\∞}\right[e^{-\mathrm{\λ}(\mathrm{ht}-\mathrm{hr})}+\frac{\mathrm{\λ}-u}{\mathrm{\λ}+u}{e}^{-\mathrm{\λ}(\mathrm{ht}+\mathrm{hr})}\left]{\mathrm{\λ}}^2{J}_0\left(\mathrm{\λr}\right)\mathrm{d\λ}\right\}$

其中，μ₀为真空磁导率， $u=\sqrt{{\mathrm{\λ}}^2+s{\mathrm{\μ}}_0/{\mathrm{\ρ}}_i^0},$ s为拉普拉斯因子。依次计算n个采样时刻下的

和

响应数据，其中时间循环指标k＝1，2，…，n；

(c)计算电阻率值为ρ_i ⁰对应各校正时刻T_k＝(t_k+t_k+nk)/2时的幅度比值理论查询值 $A^0({\mathrm{\ρ}}_i^0,T_k)=\frac{{\left[{\hat{V}}_x^0({\mathrm{\ρ}}_i^0,t_k)\right]}^2+{\left[{\hat{V}}_x^0({\mathrm{\ρ}}_i^0,t_{k+\mathrm{nk}})\right]}^2}{{\left[{\hat{V}}_z^0({\mathrm{\ρ}}_i^0,t_k)\right]}^2+{\left[{\hat{V}}_z^0({\mathrm{\ρ}}_i^0,t_{k+\mathrm{nk}})\right]}^2},$ 其中时间循环指标k＝1，2，…，n-nk；

(d)令i＝i+1，如果i＜m，则返回(a)，继续循环计算；如果i＝m时，全部m个电阻率值对应的幅度比值数据A⁰(ρ_i ⁰，T_k)，(i＝1，2，…m；k＝1，2，…，n-nk)计算完毕，保存在二维数组中作为数据查询表，进入步骤(4)；

(4)计算各校正时刻的等效均匀半空间电阻率值ρ_a(i)，i＝1，2，…，n-nk：

(a)令i＝1， $d=|A\left(T_1\right)-A^0({\mathrm{\ρ}}_1^0,T_1)|,$ j＝1，nj＝1；

(b)计算 $d_j=|A\left(T_i\right)-A^0({\mathrm{\ρ}}_j^0,T_i)|,$ 如果d_j＝0.0， ${\mathrm{\ρ}}_a\left(i\right)={\mathrm{\ρ}}_j^0,$ 转到(g)；否则，

如果d＜dj，d＝dj，nj＝j，转到(c)；

(c)j＝j+1；

(d)如果j≤m，返回(b)；如果j＞m，转到(e)；

(e)构造插值函数，如果nj＞Np，k1＝nj-Np，否则k1＝1；如果nj＜m-Np，k2＝nj＜m-Np，否则，k2＝m；则用于插值的nn＝k2-k1+1个局部电阻率序列为(ρ_k1 ⁰，ρ_k1+1 ⁰，…，ρ_p0 ⁰，…，ρ_k2 ⁰)，对应的理论查询数据为(A_k1 ⁰(ρ_k1 ⁰，T_i)，A_k1+1 ⁰(ρ_k1+1 ⁰，T_i)，…，A_k2 ⁰(ρ_k2 ⁰，T_i))；

(f)由电阻率序列(ρ_k1 ⁰，ρ_k1+1 ⁰，…，ρ_p0 ⁰，…，ρ_k2 ⁰)与查询数据(A_k1 ⁰(ρ_k1 ⁰，T_i)，A_k1+1 ⁰(ρ_k1+1 ⁰，T_i)，…，A_k2 ⁰(ρ_k2 ⁰，T_i))对应的函数关系，采用样条插值计算观测数据A(T_i)对应的电阻率值ρ_c，即得到第i个校正时刻的等效均匀半空间电阻率ρ_a(i)＝ρ_c；

(g)i＝i+1；

(h)如果i≤n-nk，返回(a)；否则，进入步骤(5)；

(5)对n-nk个校正时刻点循环计算校正高度h1下的各分量高度校正结果V′_x、V′_z：

${\hat{V}}_x^{\′}({\mathrm{\ρ}}_a\left(k\right),T_k)=-{\mathrm{\μ}}_0{S}_r\frac{M}{4\mathrm{\π}}{L}^{-1}\left\{{\∫}_0^{\∞}\right[e^{-\mathrm{\λ}(h1-\mathrm{hr})}-\frac{\mathrm{\λ}-u}{\mathrm{\λ}+u}{e}^{-\mathrm{\λ}(h1+\mathrm{hr})}\left]{\mathrm{\λ}}^2{J}_1\left(\mathrm{\λr}\right)\mathrm{d\λ}\right\}$

${\hat{V}}_z^{\′}({\mathrm{\ρ}}_a\left(k\right),T_k)=-{\mathrm{\μ}}_0{S}_r\frac{M}{4\mathrm{\π}}{L}^{-1}\left\{{\∫}_0^{\∞}\right[e^{-\mathrm{\λ}(h1-\mathrm{hr})}+\frac{\mathrm{\λ}-u}{\mathrm{\λ}+u}{e}^{-\mathrm{\λ}(h1+\mathrm{hr})}\left]{\mathrm{\λ}}^2{J}_0\left(\mathrm{\λr}\right)\mathrm{d\λ}\right\}$

其中， $u=\sqrt{{\mathrm{\λ}}^2+s{\mathrm{\μ}}_0/{\mathrm{\ρ}}_a\left(k\right)},$ 循环变量k＝1，2，…，n-nk。V′_x、V′_z即为高度校正结果，全部n-nk个校正时刻的高度校正结果完成后，计算完毕，程序结束。

Claims (1)

1.一种时间域航空电磁高度校正方法，其特征在于该方法包括以下步骤：

(1)对实测的两感生电动势分量数据进行转换

在实际飞行高度下，由两个感生电动势分量实测数据Vx、Vz计算对应各校正时刻T_k＝(t_k+t_k+nk)/2的幅值比值数据

所述k＝1，2，…，n-nk，n为采样时间个数；

(2)计算幅度比值数据查询表

给定m个离散化查询电阻率数据ρ_i ⁰(i＝1，2，...，m)，分别作为等效均匀半空间的电阻率值，先计算它的两分量航空瞬变电磁响应V_x ⁰、V_z ⁰，再计算理论幅度比值查询数据体

并将其作为计算等效均匀半空间电阻率的查询数据表保存在一个二维数组中，所述i＝1，2，…，m，k＝1，2，…，n-nk，n为采样时间个数；

(3)计算等效均匀半空间电阻率值

以各校正时刻T_k对应的两分量幅度比值数据A⁰(ρ_i，T_k)作为电阻率的单值函数关系，用局部样条插值的方法计算实测两分量数据的幅度比值数据A(T_k)对应的电阻率值，即为所求的第k个校正时刻的等效均匀半空间电阻率值，逐一计算n-nk个校正时刻的等效均匀半空间电阻率值，即完成了等效均匀半空间电阻率的计算，所述i＝1，2，…，m；k＝1，2，…，n-nk，n为采样时间个数；

(4)计算各校正时刻的感生电动势数据

根据各校正时刻的等效均匀半空间电阻率值，在给定的校正高度h1下，逐一计算全部n-nk个校正时刻对应的各分量感生电动势数据，即为最终的校正结果。

Images