基于波恩迭代法的航空瞬变电磁电导率三维反演方法与流程

本发明涉及电磁法勘探，尤其是涉及到基于波恩迭代法(borniterativemethod)的航空瞬变电磁电导率三维反演方法。

背景技术：

电磁法勘探是利用天然场源或者人工发射电流源，在大地激发感应电磁场，通过分析和处理感应的电磁场，对地下良导电矿体的分布情况进行解释的一种地球物理方法。瞬变电磁法(transientelectromagneticmethods，简称tem)是一种在地球物理中常用到的电磁勘探方法。瞬变电磁法通过人工发射脉冲电流源或其他发射波形，在大地中产生交变电磁场，并接收地下异常体随时间变化的感应信号，从接收到的信号中推断和分析地下异常结构的分布情况。航空电磁法(airborneelectromagnetic，简称aem)是电磁法勘探中重要的一项分支技术，是一套基于航空系统的时间域勘探方法，主要应用于地面直接勘探难度大、开采困难的地区或者大量森林、湖泊、沼泽覆盖的区域。

航空电磁法可根据发射源的不同分为全航空电磁法(airbornetransientelectromagneticmethods，简称atem)和半航空电磁法(groundedairbornetransientelectromagneticmethods，简称greatem)。全航空电磁法将发射线圈和接受线圈放置于无人机或者固定翼飞机上，用于发射感应电流和接收大地的异常响应。全航空电磁法的优势在于覆盖面积大、快速高效且能勘探地面难以到达的区域，缺点在于勘探深度较浅。半航空电磁法的发射源通常选择2至3公里长的电性线源铺设于地表，将线圈悬挂于无人机下方接收感应信号。这样可以得到较大的收发距，从而提高探测深度。

技术实现要素：

本发明的目的在于提供对地下异常体的电导率的三维反演算法，并将算法用于瞬变电磁系统的数据的基于波恩迭代法的航空瞬变电磁电导率三维反演方法。

本发明包括以下步骤：

1)数据预处理；

在步骤1)中，所述数据预处理的具体方法可为：在航空电磁法中，接收机在空中接收的电压信号为总场响应信号，并且接收过程中会受到一定的噪声干扰；通过对接收机接收数据进行平均，滤波和去噪等处理，获取较为理想的总场响应信号；再将计算的理论均匀半空间响应从总响应信号中扣除，得到反演所需要的二次场响应信号。

2)三维初始模型设定；

在步骤2)中，所述三维初始模型设定的具体方法可为：根据具体实验飞行的信息建立三维分层初始模型，初始模型的参数包括大地的分层信息，每层的分界面位置，每层介质的介电常数、电导率、磁导率等信息以及异常体的位置和电导率、介电常数等。

3)模型参数的更新；

在步骤3)中，所述模型参数的更新的具体方法可为：在电导率的三维反演算法中，主要是对计算得到的二次场值进行拟合，其中二次场的数值是通过地下介质和异常体的对比度函数体现的，对比度函数χ(r)的定义如下：

其中表示复介电常数，定义为r表示空间位置，j表示虚数单位，ω表示角频率。在每一次迭代计算中，通过求解成本函数得到模型参数的更新量δχn+1(r)，χn+1(r)即为对比度函数在空间的分布；计算得到第n+1次的更新量之后，通过上一次的对比度函数的分布，计算得到下一次迭代时需要带入的对比度函数，即χn+1(r)＝χn(r)+δχn+1(r)；χn(r)为第n次迭代时对比度函数的分布，χn+1(r)为第n+1迭代时对比度函数的分布；在一般反演情况下，通常采用均匀半空间为初始模型，即第一次迭代时对比度函数为0。

4)二次场值的计算；

在步骤4)中，所述二次场值的计算的具体方法可为：在建立初始模型之后，在频率域计算模型在任意接收点的二次场的场值；采用的是体积分方程法对电磁场进行求解，在第m层中，计算区域内总的电场em(r)、入射电场和散射电场满足以下表达式：

计算区域内的散射场可以通过磁矢势的分层格林函数表示，即

为第q层的发射源r'在第m层产生的电场并矢格林函数，将散射电场的表达式带入公式(1)中，并用电场通量表示总场，可以得到电场积分方程的最终计算表达式如下：

求解上述方程得到计算区域内的电场分布，即可通过公式(2)求解任意位置的电场或者磁场的场值，其中磁场的二次场值可用磁场并矢格林函数表示：

5)误差计算；

在步骤5)中，所述误差计算的具体方法可为：通过理论计算得到二次场值f^cal之后，可以实验数据或仿真数据fo^bs进行做差，得到两者之间的误差为：

6)模型更新量的计算；

在步骤6)中，所述模型更新量的计算的具体方法可为：通过对计算量和误差的定义，根据具体计算的模型，得到离散之后理论计算的二次场的改变量，即：

上式可以简化的表示为δf＝m·δχ，定义计算更新量的成本函数如下：

其中，f表示成本函数，δf表示电场或者磁场的场值误差，mn表示公式(3)中简化的矩阵，δχn+1表示第n+1次迭代的对比度的更新量，f^obs表示实验观测值或者仿真数据的场值，γ表示求解成本函数时用来约束的正则化系数，χn表示第n次迭代时对比度函数的分布，||||²表示矩阵的二范数，令成本函数的值最小，得到等效的矩阵表达式：

等效表达式中表示mn矩阵的复共轭转置，i表示单位矩阵，在迭代过程中求解等效表达式的矩阵方程即可得到每一次迭代中对比度函数分布的更新量δχn+1。

7)收敛条件判断；

在步骤7)中，所述收敛条件判断的具体方法可为：迭代收敛的条件为理论计算的二次场场值与实测数据或仿真数据之间的差值满足收敛条件即终止迭代过程，即δf＜1％。当收敛条件未到达时，重复步骤4)～7)的过程，更新模型参量；当收敛条件满足时，结束迭代过程，输出模型参量。

本发明是基于波恩迭代方法(borniterativemethod)的一种频率域快速反演方法，能应用于瞬变电磁系统的勘探技术中，对地下良导体的电导率进行反演。在瞬变电磁勘探系统中，接收到的感应信号通常包括电磁噪声、地质噪声和接收与发射机产生的系统噪声等，在实际应用中这些噪声会对异常体的感应信号产生干扰。为了保证反演结果的可靠性，需要对数据进行去噪，滤波和校正等处理。另外，本发明的反演算法利用二次场感应信号进行反演，而实际接收的信号通常为包含了直达波、大地响应和异常体响应的总场响应，所以需要对数据进行二次场提取。通过快速hankel变换和分层格林函数可以计算得到时间域的均匀半空间响应。用接收到的总场响应扣除半空间响应，提取关于地下异常体的二次场响应。再将得到的二次场响应进行快速傅立叶变换，就能获得反演所需的二次场的频谱信息。利用二次场频谱信息进行三维电导率的反演，就能获得地下异常体的分布情况。

本发明的优点在于：

一、算法在频率域中进行反演，计算速度比时域中更快；

二、本发明利用的是二次场响应的数据，这样就能去除直达波和大地响应对数据的影响。

附图说明

图1为本发明实施例信息建立三维分层初始模型分界图。

图2为本发明实施例用于仿真计算的理论模型俯视图。

图3为本发明实施例用于仿真计算的理论模型剖面图。

图4为本发明实施例的瞬变电磁系统的总场响应和背景场响应计算对比图。

图5为本发明实施例的瞬变电磁系统的二次场提取结果图。

图6为本发明实施例的仿真数据三维反演效果图。

图7为本发明实施例的理论模型仿真数据反演效果x-y剖面图。

图8为本发明实施例的理论模型仿真数据反演效果x-z剖面图。

图9为本发明实施例的理论模型仿真数据反演效果y-z剖面图。

具体实施方式

以下结合附图对本发明的技术方案作进一步说明。

具体实施步骤如下：

(1)数据预处理

在航空电磁法中，接收机在空中接收的电压信号为总场响应信号，并且接收过程中受到一定的噪声干扰。通过对接收数据的平均，滤波和去噪处理，可以获得较为理想的总场响应信号。再将计算的理论均匀半空间响应从总响应中信号扣除，得到反演所需要的二次场响应信号。

(2)三维初始模型设定

根据具体实验飞行的信息建立三维分层初始模型，初始模型的参数包括大地的分层信息，如每层的分界面位置，每层介质的介电常数、电导率，磁导率等信息以及异常体的位置和电导率、介电常数等，如图1所示，其中具体的初始参数定义如下：ε1,ε2...εn：依次为地下第一层至第n层模型的介电常数；σ1,σ2...σn：依次为地下第一层至第n层模型的电导率；z1,z2…zn依次为地下分层介质的位置信息；其中地下异常体完全嵌入其中某一层，异常体的介电常数、电导率依次为εr,σr。

(3)模型参数的更新

在电导率的三维反演计算中，主要是利用二次场值进行拟合。其中二次场的数值是通过地下介质和异常体的对比度函数体现的。其中对比度函数χ(r)的定义如下：

表示复介电常数，定义为r表示空间位置，j表示虚数单位，ω表示角频率。每一次迭代计算中，通过求解成本函数得到模型参数的更新量δχn+1(r)，χn+1(r)即为对比度函数在空间的分布。计算得到第n+1次的更新量之后，通过上一次的对比度函数的分布，计算得到下一次迭代时需要带入的对比度函数，即χn+1(r)＝χn(r)+δχn+1(r)。χn(r)为第n次迭代时对比度函数的分布，χn+1(r)为第n+1迭代时对比度函数的分布。在一般反演情况下，通常采用均匀半空间为初始模型，即第一次迭代时对比度函数为0。

(4)二次场值的计算

在建立了初始模型之后，在频率域计算模型在任意接收位置的二次场的场值。这里我们采用的是体积分方程法对电磁场进行求解，在第m层中，计算区域内总的电场em(r)、入射电场和散射电场满足以下表达式：

计算区域内的散射场可以通过磁矢势的分层格林函数表示，即

为第q层的发射源r'在第m层产生的电场并矢格林函数，将散射电场的表达式带入公式(5)中，并用电场通量表示总场，可以得到电场积分方程的最终计算表达式如下：

求解上述方程得到计算区域内的电场分布，即可通过公式(6)求解任意位置的电场或者磁场的场值。

(5)误差计算

通过理论计算得到二次场值f^cal之后，可以实验数据或仿真数据f^obs进行做差，得到两者之间的误差为：

(6)模型更新量的计算

通过对计算量和误差的定义，根据具体计算的模型，可以得到离散之后理论计算的二次场的改变量，即：

上式可以简化的表示为δf＝m·δχ。定义计算更新量的成本函数如下：

f表示成本函数，δf表示电场或者磁场的场值误差，mn表示公式(7)中简化的矩阵，δχn+1表示第n+1次迭代的对比度的更新量，f^obs表示实验观测值或者仿真数据的场值，γ表示求解成本函数时用来约束的正则化系数，χn表示第n次迭代时对比度函数的分布，||||²表示矩阵的二范数。求解成本函数的方程，令成本函数取值最小，可以得到等效表达式：

等效表达式中表示mn矩阵的复共轭转置。在迭代过程中求解等效表达式的矩阵方程即可得到每一次迭代中对比度函数分布的更新量δχn+1。

(7)收敛条件判断

迭代收敛的条件为理计算的论二次场场值与实测数据或仿真数据之间的差值满足收敛条件即终止迭代过程，即error＜1％。当收敛条件未到达时，则重复步骤(4)到步骤(7)的过程，更新模型参量；当收敛条件满足时，结束迭代过程，输出模型参量。

本发明实施例用于仿真计算的理论模型俯视图参见图2，本发明实施例用于仿真计算的理论模型剖面图参见图3，本发明实施例的瞬变电磁系统的总场响应和背景场响应计算对比图参见图4，本发明实施例的瞬变电磁系统的二次场提取结果图参见图5，本发明实施例的仿真数据三维反演效果图参见图6，本发明实施例的理论模型仿真数据反演效果x-y剖面图参见图7，本发明实施例的理论模型仿真数据反演效果x-z剖面图参见图8，本发明实施例的理论模型仿真数据反演效果y-z剖面图参见图9。