基于改进迭代法的地空时间域电磁系统高精度下延拓方法与流程

本发明涉及一种地球物理勘探领域的基于改进迭代法的地空时间域电磁系统高精度下延拓方法，尤其适用于地空时间域电磁系统探测方法。

背景技术

地空时间域电磁探测系统的发射装置位于地表，通常以电性源为主，向地下发射大功率电流。接收装置固定在无人机、飞艇等无人操纵的飞行器中，通过地面和空中协同测量方式，获取含有地下信息的z分量感应电动势，相比于航空电磁探测系统，它具有探测深度更深的优势，而相比于地面电磁探测系统，更适合在高山、海陆交互带、滩涂地区等人不易进入的区域开展大规模的电磁勘探工作，具有广泛的应用前景。

空间域向下延拓迭代法作为一种广泛使用的向下延拓方法，由著名地球物理学家徐世浙院士提出，他将重、磁资料解释中的向下延拓问题巧妙转化为向上延拓问题，通过采用差值对延拓面上的位场进行校正和反复迭代，有效解决了位场数据的向下延拓问题。

中国专利201010573168.3公开了一种时间域地空电磁探测系统及标定方法。测量时，通过去除大地背景场数据，提取含几何误差的纯闭合异常环电磁信号，输入数据质量评价系统后，与计算的闭合异常环理论值进行比较、拟合，确定系统误差、几何参数误差以及时间域地空电磁探测系统的探测分辨率，实现电磁探测系统性能的测试和标定。

中国专利201310322111.x公开了一种基于积分迭代算法的航空重力向下延拓方法，在对测区航空重力测量数据进行向下延拓时，将航空重力测量数据作为积分迭代的初始值进行向上延拓；将向上延拓得到的结果与航空重力测量数据求差后作为向下延拓结果的反馈量；将反馈量与初始值叠加得到新的迭代初始值进行下一次迭代，该方法可提高数据的精度与置信度。

中国专利201610005275.3公开了一种基于逆延拓理论的地空时域电磁数据高度校正方法，采用无人机或飞艇在空中进行飞行测量，将测量的数据进行电磁噪声去除、双极性叠加取样、插值后，应用逆延拓理论，将空中z＝-h高度处测量的磁场变化率进行延拓计算，获得地面z＝0处的磁场变化率数据，进而实现飞行高度的校正。

上述三个专利，其中专利201010573168.3公布了时间域地空电磁探测系统及标定方法，包括系统的构成、几何参数误差、探测分辨率的研究等，但是并未涉及空中感应电动势的向下延拓以及视电阻率-视深度计算及成像问题。专利201310322111.x公布了采用积分迭代算法解决航空重力测量数据的下延拓问题，航空重力测量数据为磁场数据，而地空系统测量的时间域电磁数据为感应电动势，数据本身存在差异，因此为了得到高精度的下延拓数据必须对迭代方法进行修正，使其适合地空系统测量数据的下延拓。专利201610005275.3发布的逆延拓方法将fft向上延拓中的延拓算子中的高度参数变号并加入了i(kxx+kyy)项，但是仍然无法改变e指数项在高频时的不稳定性，与未延拓时相比精度有所提高但是变化不大。

技术实现要素：

本发明所要解决的技术问题在于提供一种基于改进迭代法的地空时间域电磁系统高精度下延拓方法，通过理论计算数据，利用鱼群方法优化迭代公式中的三个参数，利用修正后的迭代式对地表场进行逐步校正，验证方法的正确性，最终将其应用在实测数据中，通过对数据处理后的空中感应电动势进行向下延拓、以及视电阻率-视深度成像，提高地空电磁数据的解释精度。

本发明是这样实现的，一种基于改进迭代法的地空时间域电磁系统高精度下延拓方法包括：

1)计算长接地导线源在地表z＝0平面的感应电动势，接收点之间需等间距分布，利用fft向上延拓理论计算z＝h(h>0)平面的空中感应电动势；

2)采用鱼群智能方法确定改进迭代法中的三个参数：α、β、γ；

3)通过改进迭代法，下延拓出z＝0平面感应电动势，与步骤1)计算出的地表理论感应电动势进行比较，验证方法的正确性；

4)采用接收线圈传感器测量z＝h平面感应电动势，为保证下延拓精度，测线条数最少为五条，且保持平行；

5)将测得数据进行叠加、采样处理，判断处理后的数据是否为等间距，否则进行插值；

6)采用改进迭代法，下延拓出z＝0平面感应电动势，计算视电阻率，视深度，并完成视电阻率-视深度成像。

进一步地，步骤2)中，通过鱼群方法选择最优的α、β、γ三个参数值，首先完成相对误差目标函数的构造：

η＝max{[uⁿ(x,y,z＝h)-f(x,y,z＝h)]/uⁿ(x,y,z＝h)×100}(1)

其中f(x,y,z＝h)为步骤1)通过fft向上延拓出的z＝h平面感应电动势，uⁿ(x,y,z＝h)为向下延拓的地表感应电动势再进行fft向上延拓得到的空中感应电动势，进行鱼群算法的初始化，其中包括α、β、γ三个参数值的初始化，通过预设最大计算次数m、精度值q、逐步迭代确定最优参数值，首先判断η＜q，若满足条件迭代结束，否则按照最大计算次数m完成最优参数确定。

进一步地，步骤3)中，采用改进迭代法实现空中感应电动势的下延拓计算，进一步提高地-空电磁系统数据解释的精度，将步骤1)中延拓出的z＝h平面感应电动势f(x,y,z＝h)作为地表z＝0平面的初始感应电动势，为了保证下延拓在高频时解的稳定性，不能采用fft向下延拓方法，因此采用改进迭代法实现空中场的向下延拓，并将此次延拓出的z＝h平面感应电动势uⁿ(x,y,z＝h)与f(x,y,z＝h)之间的差值作为误差项：

errorⁿ＝uⁿ(x,y,z＝h)-f(x,y,z＝h)(2)

对向下延拓出的地表感应电动势uⁿ(x,y,z＝0)进行修正：

其中n为迭代次数，增加积分项和微分项可以在相同迭代次数时增加计算精度；

对uⁿ⁺¹(x,y,z＝0)进行向上延拓、计算误差项、修正，直到errorⁿ小到可以忽略为止，将此时得到的uⁿ⁺¹(x,y,z＝0)作为最终的z＝0平面感应电动势，完成空中感应电动势的向下延拓。

本发明与现有技术相比，有益效果在于：本发明针对逆延拓理论的地空时域电磁数据高度校正方法延拓精度不高，并且对于电磁信号中存在的高频分量延拓后数据不具有稳定性，以及迭代积分法的参数手动选取、迭代次数高等不足，本专利提出采用修正后的迭代法实现空中感应电动势的下延拓，修正项由只有误差项一项增加为误差积分项、误差微分项共三项误差项，其中对应的三个参数由鱼群方法智能寻找，是一种迭代次数少精度高的下延拓方法。

附图说明

图1是地空电磁系统示意图；

图2基于改进迭代法的地空时间域电磁系统高精度下延拓方法流程图；

图3是基于鱼群算法的参数确定方法流程图；

图4是改进迭代法的正确性验证流程图；

图5是向下延拓精度验证图；

图6是视电阻率-视深度成像图；

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图和实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

实施例

参见图2结合图1所示，一种基于改进迭代法的地空时间域电磁系统高精度下延拓方法，包括：

1)计算长接地导线源在地表z＝0平面的感应电动势，接收点之间需等间距分布，利用fft向上延拓理论计算z＝h(h>0)平面的空中感应电动势；

如图1所示，步骤1)中，根据长接地导线源层状大地z方向磁场响应计算式：

其中l为长接地导线半长度，r＝[(x-x′)²+y²]^1/2，为反射系数，u0≈λ。再利用guptasarma频时变换方法将其从频率域变换到时间域，从而得到地表理论感应电动势f(x,y,z＝0)。在计算过程中为了实现向上延拓，x与y方向上的节点间距需为等间距。计算各节点响应后采用基于fft的向上延拓方法实现感应电动势的延拓：

f(kx，ky，z＝h)＝g(kx，ky，z＝h)·f(kx，ky，z＝0)(5)

其中f(kx,ky,z＝0)、f(kx,ky,z＝h)分别为f(x,y,z＝0)、f(x,y,z＝h)进行二维傅里叶变换后结果，其中延拓算子h为接收线圈距离地面高度，kx、ky为波数，最后再将f(kx,ky,z＝h)其变换到空间域。

2)采用鱼群智能方法确定改进迭代法中的三个参数：α、β、γ；

如图3所示，步骤2)中，通过鱼群算法选择最优的α、β、γ三个参数值，采用步骤1计算的长接地导线源在地表z＝0平面延拓出的感应电动势f(x,y,z＝h)，与u(x,y,z＝0)延拓出的感应电动势之间的相对误差趋于一个极小值，鱼群算法的具体步骤如下：

①完成相对误差目标函数的构造；

η＝max{[u(x,y,z＝h)-f(x,y,z＝h)]/uⁿ(x,y,z＝h)×100}(1)

②对鱼群个数n、感知距离visual、觅食行为最大尝试次数try_number、预设最大相对误差精度值q，最大计算次数m等参数进行赋值，实现鱼群的初始化；

③进行聚群行为和追尾行为，若未满足拥挤度和最优食物浓度的条件，则进行觅食行为，若未在视觉范围内找到更优食物浓度则进行反复尝试，若达到最大尝试次数仍未找到更优位置则进行随机行为，确定α、β、γ参数值；

④一次计算结束后，更新鱼群数据，判断是否达到预设精度或最大计算次数，若达到输出参数，否则转入③继续执行。

3)通过改进迭代法，下延拓出z＝0平面感应电动势，与步骤1)计算出的地表理论感应电动势进行比较，验证方法的正确性；

如图4所示，步骤3)中采用改进迭代法实现空中感应电动势的下延拓计算，进一步提高地-空电磁系统数据解释的精度，将步骤一中延拓出的z＝h平面感应电动势f(x,y,z＝h)作为地表z＝0平面的感应电动势uⁿ(x,y,z＝0)初值，为了保证下延拓在高频时解的稳定性，仍然采用fft向上延拓方法，将此次延拓出的z＝h平面感应电动势uⁿ(x,y,z＝h)与f(x,y,z＝h)之间的差值作为误差项：

errorⁿ＝uⁿ(x,y,z＝h)-f(x,y,z＝h)(2)

对地表感应电动势uⁿ(x,y,z＝0)进行修正：

其中n为迭代次数，增加积分项和微分项可以在相同迭代次数时增加计算精度。

对uⁿ⁺¹(x,y,z＝0)进行向上延拓、计算误差项、修正，直到errorⁿ小到可以忽略为止，将此时得到的uⁿ⁺¹(x,y,z＝0)作为最终的z＝0平面感应电动势，完成空中感应电动势的向下延拓，延拓结果如图5所示。

4)采用接收线圈传感器测量z＝h平面感应电动势，为保证下延拓精度，测线条数最少为五条，且保持平行；

步骤4)中，接地导线长度为2000m，电流为20a-40a，接收线圈距离地面高度为30m，无人机与接收线圈之间的连接线长度为10m-15m。不同测线之间平行，为保证精度测线条数最少为五条，采用接收线圈对测点进行感应电动势测量。

5)将测得数据进行叠加、采样处理，判断处理后的数据是否为等间距，否则进行插值；

步骤5)中，地空电磁信号中存在随机噪声，由于信号具有周期性而噪声具有随机性，因此可采用周期性叠加方法消除电磁数据中的随机噪声，叠加方法如下：

n为叠加次数，f(t)为地空电磁信号，本实施例中的叠加次数n＝5。

对于数据中的采样部分，将电磁数据按对数等间隔截取，在切开的每一个时间段内做算术平均，判断叠加与采样后的各测点数据是否为等间距，若非等间距采用双三次样条插值方法对其进行插值运算。

6)采用改进迭代法，下延拓出z＝0平面感应电动势，计算视电阻率，视深度，并完成视电阻率-视深度成像。

步骤6)中，z方向电偶极源的电磁响应公式为：

其中i为发射电流，ds为电偶极子长度，σ为电导率，μ为磁导率，r为收发距，t为采样时间。

长导线源的电磁响应可认为是电偶极子的叠加，考虑线圈的面积s，并将bz＝μhz代入公式(6)，可得长导线源的响应：

公式(6)为θ的函数，通过计算我们能够得到θ的值，视电阻率ρ和θ的关系可以通过下式表示出来：

本实施例中采用扩散深度计算电磁波的传播深度，考虑关断时间的影响，扩散深度可以表示为：

其中t是采样时间，toff是关断时间。图6为地空时间域电磁数据下延拓后的视电阻率-视深度成像结果。

图5为采用图2所示的本发明一个实施例的基于改进迭代法延拓出的地表感应电动势与理论计算出的电动势两者之间的相对误差，为0.01％，延拓后的结果与地表解析式计算结果可以忽略不计，充分验证了基于地空时间域电磁系统的高精度下延拓方法的正确性，为地空电磁数据的高精度解释提供了新方法。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。