CN106405665B - 基于dbim的瞬变电磁电导率反演方法 - Google Patents

Abstract

基于DBIM的瞬变电磁电导率反演方法，涉及地球物理勘探。包括以下步骤：1)读取观测数据；2)建立初始模型；3)更新模型参数；4)计算模型场值；5)计算误差；6)计算Frechet导数；7)计算更新量；8)判断收敛条件。从瞬变电磁的频谱信息出发，只要提取接收信号的准确频谱信息，就可以进行反演。基于DBIM方法建立迭代反演过程，最后反演的结果能很好与实际数据相吻合，可以大大提高瞬变电磁系统的计算速度和反演精度。不仅适用于半航空瞬变电磁系统，而且适用于全航空瞬变电磁系统。

Description

基于DBIM的瞬变电磁电导率反演方法

技术领域

本发明涉及地球物理勘探，尤其是涉及一种基于DBIM的瞬变电磁电导率反演方法。

背景技术

瞬变电磁法(Transient Electromagnetic Methods)又称时间域电磁方法(TimeDomain Electromagnetic Methods)，简称TDEM或TEM。航空电磁(AEM,AirborneElectromagnetics)法是在地面电磁法基础上发展起来的一种空中测量的电磁法，被广泛用来进行资源勘探，并且其非常适合探测一些复杂地形和人员难以达到的地区。其中半航空瞬变电磁(Semi-Transient Electromagnetic Methods)是为了克服全航空瞬变电磁的探测深度有限所发展出的新型航空电磁法，其是在地面铺设几公里长度的电性源，然后利用无人机或者直升机在空中接收磁场，通过分析接收到磁场的信息，来反演出地下介质的情况。

航空电磁法是基于岩石电性及磁性差异，利用电磁感应原理，以固定翼飞机或直升机等飞行器作为运载工具，实时地球物理探测的勘探方法，其方法具有高效、经济、适应性强等特点，能够广泛应用于地面大面积的矿产普查、水工环境普查、详查和精细测量等工作。

发明内容

本发明的目的在于提供一种在瞬变电磁系统中能准确反演地层电导率，在频率域中直接进行反演，不需要进行频时转换，在数据处理过程比较简单，只需要准确提取相应频率的频谱即可，不仅适用于半航空瞬变电磁系统，同样适用于全航空瞬变电磁系统的基于DBIM的瞬变电磁电导率反演方法。

本发明包括以下步骤：

1)读取观测数据，具体步骤如下：

将实际接收机接收数据作为输入数据。

在步骤1)中，所述输入数据包括电压数据或磁场数据。

2)建立初始模型，具体步骤如下：

根据已有的信息建立初始模型，初步定义下列初始参数：

ε＝(ε₁,ε₂...ε_n)；σ＝(σ₁,σ₂...σ_n)；h＝(h₁,h₂...h_n)

其中，ε：为地下分层模型的介电常数矩阵，ε₁,ε₂...ε_n：为地下第一层、第二层…第n层模型的介电常数；σ：为地下分层模型的电导率矩阵，σ₁,σ₂...σ_n：为地下第一层、第二层…第n层模型的电导率；h：为地下分层模型的电导率矩阵，h₁,h₂...h_n：为地下第一层、第二层…第n层模型的厚度；

在步骤2)中，所述初始模型包括分层信息、每层的厚度、每层介质电导率、介电常数等信息。

3)更新模型参数，具体步骤如下：

为迭代过程中更新模型的参数过程，每一次迭代的前一次都会得到一个模型的更新量，这里是δσ＝(δσ₁,δσ₂...δσ_n)，其中δσ为模型电导率更新矩阵，δσ₁,δσ₂...δσ_n为地下第一层、第二层…第n层的电导率更新量；得到电导率更新量之后，有：

σ^k+1＝σ^k+δσ^k+v

其中，σ^k+1为第k+1次迭代所需要的电导率参数矩阵；σ^k为第k次迭代过程中所需要的迭代矩阵；δσ^k+1为第k次迭代过程中所得到的模型电导率更新矩阵，其中第一次迭代，模型的更新量为0；

4)计算模型场值，具体步骤如下：

得到模型后，对模型的场值进行计算，采用频率域的计算方法：

E＝<G^EJ；J>

H＝<G^HJ；J>

其中，E和H表示电场和磁场；G^EJ和G^HJ为电流源电场并矢Green函数和磁场Green函数，<；>为点乘积分和；

5)计算误差，具体步骤如下：

得到模型的场值后，对实际测量的场值与模型计算的场值做差，得到两者之间的误差；

6)计算Frechet导数，具体步骤如下：

Frechet导数矩阵由以下方程求出：

其中为对应于地下介质所产生的等效源，为由源J'所产生的场值的变化量，得到去谱域的解后对其进行Fourier逆变换得到其频域的解。如下为在y方向源的作用下所产生的磁场的x分量：

δH_x＝δH_x1+δH_x2+δH_x3

δH_x为磁场的x方向分量，δH_x1、δH_x2、δH_x3分别为其分量的第一项、第二项与第三项，具体表达式如下：

上式中J_y表示y方向的源的强度；y_-l与y_l表示沿y方向源的起点与终点；y为接收点的y坐标；x'表示源的x坐标，x、y表示接收点的x和y坐标；表示为沿着第n层的上表面积分到n层的下表面，即第n层厚度的积分；i_rip(r)V_rip(r)表示位于z'的1V串联脉冲电流源(或1A并联脉冲电压源)在深度z处的电流和电压，下标p＝e,h分别表示横磁波与横电波；J₀、J₁表示第0类Bessel函数和第1类Bessel函数；其中，k_x和k_y分别为x方向和y方向的波数；δσ_n为地下第n层的电导率参数变化量。

对于上述的场值变化δH_x＝F·δσ，其中F即为地层参数的Frechet导数。

7)计算更新量，具体步骤如下：

对于上述地层模型，建立如下成本函数，

其中，C表示成本函数(Cost Function)；δf表示场值的误差；F为Frechet导数矩阵；δσ^k+1为第k次迭代所得到的k+1次的电导率更新量；f^obs为测量到的场值；σ^k为第k次迭代所得到的电导率矩阵；γ²表示正则化系数；||||²表示矩阵的二范数。

为了使成本函数取得最小值，等效于解如下方程：

其中，为F矩阵的转置，若其为复数矩阵，则表示为共轭转置。

由上式可以得到每一次迭代过程中电导率的更新量δσ^k+1。

8)判断收敛条件，具体步骤如下：

若未达到收敛条件，则返回步骤3)更新数据的参数；若达到收敛条件，则结束反演过程。

本发明应用于瞬变电磁快速反演，是利用瞬变电磁系统接收到的信号的频谱信息在频率域基于变形波恩迭代方法(DBIM Distorted Born iterative methods)的一种快速反演方法。

在瞬变电磁系统中，发射信号为双极性方波或者是半正弦信号，同时利用接收线圈接收一次场或者是二次场的时域信息，本发明可以使用任意发射信号，只要准确记录发射信号和接收信号，就可以准确的求得对应发射信号下的理论接收波形；对于时域的接收信号对其进行简单的去直流，去噪等操作后做FFT运算，得到其对应频率的频谱信息，就可以用来进行反演运算。

本发明不仅适用于半航空瞬变电磁系统，而且适用于全航空瞬变电磁系统。本发明从瞬变电磁的频谱信息出发，只要提取接收信号的准确频谱信息，就可以用本发明进行反演。本发明基于DBIM方法建立迭代反演过程，最后反演的结果能很好与实际数据相吻合，从而可以大大提高瞬变电磁系统的计算速度和反演精度。

附图说明

图1为本发明的系统反演流程图。

图2为瞬变电磁地下模型分层图。

图3为本发明理论模型5层的反演效果图。

图4为本发明理论模型10层的反演效果图。

图5为本发明实际数据效果图。

具体实施方式

以下将结合附图对本发明的技术方案及其突出效果作进一步说明。

参见图1，本发明的具体实施步骤如下：

(1)读取观测数据

将实际接收机接收数据作为输入数据，这些数据可以是电压数据，也可以是接收到的磁场数据。

(2)建立初始模型

可根据已有的信息建立一个简单的初始模型，初始模型需要包括分层信息，每层的厚度，每层介质电导率、介电常数等信息，参见图2，在图2中，Z₁、Z₂、…Z_n-1为层界面的高度，线源可采用电流源。

初步定义下列初始参数：

ε＝(ε₁,ε₂...ε_n)；σ＝(σ₁,σ₂...σ_n)；h＝(h₁,h₂...h_n)

ε：为地下分层模型的介电常数矩阵，ε₁,ε₂...ε_n：为地下第一层、第二层…第n层模型的介电常数；σ：为地下分层模型的电导率矩阵，σ₁,σ₂...σ_n：为地下第一层、第二层…第n层模型的电导率；h：为地下分层模型的电导率矩阵，h₁,h₂...h_n：为地下第一层、第二层…第n层模型的厚度。

(3)更新模型参数

为迭代过程中更新模型的参数过程，每一次迭代的前一次都会得到一个模型的更新量，这里是δσ＝(δσ₁,δσ₂...δσ_n)，其中δσ为模型电导率更新矩阵，δσ₁,δσ₂...δσ_n为地下第一层、第二层…第n层的电导率更新量。得到电导率更新量之后，有：

σ^k+1＝σ^k+δσ^k+1

其中σ^k+1为第k+1次迭代所需要的电导率参数矩阵；σ^k为第k次迭代过程中所需要的迭代矩阵；δσ^k+1为第k次迭代过程中所得到的模型电导率更新矩阵。其中第一次迭代，模型的更新量为0。

(4)模型场值计算

得到模型后，就要对模型的场值进行计算，这里采用频率域的计算方法：

E＝<G^EJ；J>

H＝<G^HJ；J>

这里E和H表示电场和磁场；G^EJ和G^HJ为电流源电场并矢Green函数和磁场Green函数，<；>为点乘积分和。

(5)误差计算

得到模型的场值后，需要对实际测量的场值与模型计算的场值做差，得到两者之间的误差。

(6)Frechet导数计算

其中Frechet导数矩阵可由以下方程求出：

其中为对应于地下介质所产生的等效源，为由源J'所产生的场值的变化量，得到去谱域的解后对其进行Fourier逆变换得到其频域的解。如下为在y方向源的作用下所产生的磁场的x分量：

δH_x＝δH_x1+δH_x2+δH_x3

δH_x为磁场的x方向分量，δH_x1、δH_x2、δH_x3分别为其分量的第一项、第二项与第三项，具体表达式如下：

上式中J_y表示y方向的源的强度；y_-l与y_l表示沿y方向源的起点与终点；y为接收点的y坐标；x'表示源的x坐标，x、y表示接收点的x和y坐标；表示为沿着第n层的上表面积分到n层的下表面，即第n层厚度的积分；i_rip(r)V_rip(r)表示位于z'的1V串联脉冲电流源(或1A并联脉冲电压源)在深度z处的电流和电压，下标p＝e,h分别表示横磁波与横电波；J₀、J₁表示第0类Bessel函数和第1类Bessel函数；其中，k_x和k_y分别为x方向和y方向的波数；δσ_n为地下第n层的电导率参数变化量。

对于上述的场值变化δH_x＝F·δσ，其中F即为地层参数的Frechet导数。

(7)更新量计算

对于上述地层模型，建立如下成本函数，

其中C表示成本函数(Cost Function)；δf表示场值的误差；F为Frechet导数矩阵；δσ^k+1为第k次迭代所得到的k+1次的电导率更新量；f^obs为测量到的场值；σ^k为第k次迭代所得到的电导率矩阵；γ²表示正则化系数；||||²表示矩阵的二范数。

为了使成本函数取得最小值，等效于解如下的方程

其中：为F矩阵的转置，若其为复数矩阵，则表示为共轭转置。

由上式可以得到每一次迭代过程中电导率的更新量δσ^k+1。

(8)收敛条件判断

判断收敛条件，若未达到收敛条件，则返回步骤(3)更新数据的参数；若达到收敛条件，则结束反演过程。

本发明在理论上具有较高的反演精度。图3为反演5层结果对比图，从图3中可以看出，反演结果可以很准确反演出地下电导率参数，并且反演的初值可以任意设定。这里仅使用1个观察点和20个频率点对分层情况进行反演。图4为反演10层的结果对比图，假设地下为5层的分层模型，而实际分10层来反演，结果对比图如图3，可以看出在10层结果也是可以很好与地下分层情况相符合的。图5为山东昌邑地区的实测数据的效果图，为山东昌邑某地区的铁矿采集区的一条半航空瞬变电磁数据反演结果图。

Claims (3)

1.基于DBIM的瞬变电磁电导率反演方法，其特征在于包括以下步骤：

1)读取观测数据，具体步骤如下：

将实际接收机接收场值作为输入数据；

2)建立初始模型，具体步骤如下：

根据已有的信息建立初始模型，初步定义下列初始参数

ε＝(ε₁,ε₂...ε_n)；σ＝(σ₁,σ₂...σ_n)；z＝(z₁,z₂...z_n)

其中，ε为地下分层模型的介电常数矩阵，ε₁,ε₂...ε_n为地下第一层、第二层…第n层模型的介电常数；σ为地下分层模型的电导率矩阵，σ₁,σ₂...σ_n为地下第一层、第二层…第n层模型的电导率；z为地下分层模型的层界面深度矩阵，z₁,z₂...z_n为地下第一层、第二层…第n层模型层界面深度；

3)更新模型参数，具体步骤如下：

在迭代过程中更新模型参数，并在每次迭代的前一次得到一个模型的更新量δσ＝(δσ₁,δσ₂...δσ_n)，其中δσ为模型电导率更新矩阵，δσ₁,δσ₂...δσ_n为地下第一层、第二层…第n层的电导率更新量；得到电导率更新量之后，有：

σ^k+1＝σ^k+δσ

其中，σ^k+1为第k+1次迭代所需要的电导率参数矩阵；σ^k为第k次迭代过程中所需要的迭代矩阵；δσ为第k次迭代过程中所得到的模型电导率更新矩阵，第一次迭代时，模型的更新量为0；

4)计算模型场值，具体步骤如下：

得到模型后，对模型的场值进行计算，采用频率域的计算方法：

E＝<G^EJ；J>

H＝<G^HJ；J>

其中，E和H表示电场和磁场；G^EJ和G^HJ是电流源电场的并矢Green函数和磁场Green函数，均为点乘积分和；

5)计算误差，具体步骤如下：

得到模型的场值后，与实际接收机接收场值做差，得到两者之间的差值；

6)计算Frechet导数，具体步骤如下：

Frechet导数矩阵由以下方程求出：

其中为对应于地下介质所产生的等效源，为由等效源J'所产生的场值的变化量；在y方向等效源的作用下所产生的磁场的x分量：

δH_x＝δH_x1+δH_x2+δH_x3

其中，δH_x为磁场的x方向分量，δH_x1、δH_x2、δH_x3分别为其分量的第一项、第二项与第三项，具体表达式如下：

上式中J_y表示y方向的等效源的强度；y_-l与y_l表示沿y方向等效源的起点与终点坐标；y为接收点的y坐标；x'表示等效源的x坐标，x、y表示接收点的x和y坐标；表示沿着第n层的上表面积分到n层的下表面，即第n层厚度的积分；i_rie(r)表示位于等效源的Z坐标z'处的1V串联脉冲电流源在深度z处的电流，V_sie(s)表示位于等效源的Z坐标z'处的1A并联脉冲电压源在深度z处的电压；J₀、J₁表示第0类Bessel函数和第1类Bessel函数；其中，k_x和k_y分别为x方向和y方向的波数；δσ_n为地下第n层的电导率参数变化量；

场值变化由公式δH_x＝F_x·δσ计算得出，其中F_x即为地层参数对于x方向的Frechet导数；

7)计算更新量，具体步骤如下：

对于步骤2)中建立的地层模型，建立如下成本函数：

其中，C表示成本函数；δf表示场值的误差；F为Frechet导数矩阵；δσ^k为第k次迭代所得到的k次的电导率更新量；f^obs为测量到的场值；σ^k为第k次迭代所得到的电导率矩阵；γ²表示正则化系数；

为了使成本函数取得最小值，等效于求解如下方程：

其中，为F矩阵的转置，若其为复数矩阵，则表示为共轭转置；由上式得到每一次迭代过程中电导率的更新量δσ^k+1；

8)判断收敛条件，具体步骤如下：

若未达到收敛条件，则返回步骤3)更新模型参数；若达到收敛条件，则结束反演过程。

2.如权利要求1所述基于DBIM的瞬变电磁电导率反演方法，其特征在于在步骤1)中，所述输入数据包括电压数据或磁场数据。

3.如权利要求1所述基于DBIM的瞬变电磁电导率反演方法，其特征在于在步骤2)中，所述初始模型包括分层信息、每层的厚度、每层介质电导率、介电常数。