CN113933905B - 一种圆锥型场源瞬变电磁反演方法 - Google Patents

Abstract

本发明公开了一种圆锥型场源瞬变电磁反演方法，属于瞬变电磁探测技术领域。本发明基于OCCAM反演算法，提出了一种圆锥型场源瞬变电磁反演方法。OCCAM法是一种带平滑约束的最小二乘反演，受初始模型影响小，运算稳定收敛，是一种有效的数据反演处理方法。本发明的反演方法不需要依赖于初始模型。对初始模型没有特别要求，且计算效率高。本发明采用OCCAM法得到的反演电阻率更接近模型电阻率。本发明反演采用简化后的拉格朗日乘子算法使得反演仍稳定、速度更快、结果更真实。本发明方法使得反演深度更接近层界面，尤其对目标层底界面的分辨能力获得加强，表现出较好的分层效果。

Description

一种圆锥型场源瞬变电磁反演方法

技术领域

本发明涉及一种圆锥型场源瞬变电磁反演方法，具体涉及一种基于OCCAM反演算法的圆锥型场源瞬变电磁数据的反演方法，属于瞬变电磁探测技术领域。

背景技术

瞬变电磁法(TEM)逐渐被应用于矿山巷道、交通隧道和地下空间探测，受到有限地下空间的限制，只能采用小尺寸发射-接收装置，如多匝小回线装置。杨海燕等发明了一种圆锥型场源装置，如图1所示，与多匝小回线装置相比，圆锥型装置具有互感低、关断时间短和“盲区”范围小的优点。但是，用于处理圆锥型场源探测资料的方法较少。

目前，用于圆锥型场源瞬变电磁的资料处理方法，主要有如下几种：

(1)杨海燕,李锋平,岳建华,刘旭华,赵海娇.基于“烟圈”理论的圆锥型场源瞬变电磁优化反演.中国矿业大学学报,2016,45(06):1230-1237.该论文提出了将“烟圈”反演结果转换到频率域作为初始模型，采用最小二乘正则化方法对反演结果进行改进的反演方案。

(2)杨海燕,李锋平,Chen Shen-En,岳建华,郭福生,陈晓,张华.An inversion oftransient electromagnetic data from a conical source.Applied Geophysics,2018,15(Z1):545-555+570.该论文提出了根据“烟圈”反演结果设定初始模型的阻尼最小二乘反演方法。对瞬变电磁正演所获得的视电阻率进行“烟圈”反演，以“烟圈”反演结果设定的初始模型在频率域内进行了阻尼最小二乘反演。同时，提出了以“烟圈”反演曲线的极值点和拐点为基点，建立阻尼最小二乘反演初始模型的方法。

(3)黄刚,杨海燕,李锋平,彭峰,焦俊俊,刘建鹏.圆锥型场源瞬变电磁数据AWPSO算法优化反演.地球物理学进展:1-13[2021-09-08].http://kns.cnki.net/kcms/detail/11.2982.p.20201109.1151.072.html.该论文提出了采用粒子群优化算法(PSO)和神经网络算法(BP)对圆锥型场源探测数据进行反演的方案，改进了一种基于神经网络算法Sigmoid函数的自适应加权粒子群优化(AWPSO)算法。

但上述反演方法的分层能力比较弱，为了使反演深度更接近层界面，加强对目标层底界面的分辨能力，获得较好的分层效果，需要提出一种新的圆锥型场源瞬变电磁反演方法。

发明内容

本发明所要解决的技术问题是：提供一种圆锥型场源瞬变电磁反演方法，该方法使得反演深度更接近层界面，尤其对目标层底界面的分辨能力获得加强，表现出较好的分层效果。

本发明为解决上述技术问题采用以下技术方案：

一种圆锥型场源瞬变电磁反演方法，该方法包括如下步骤：

步骤1，读入初始模型m₁即初始电阻率，根据初始模型m₁构造粗糙度矩阵，设置初始拉格朗日乘子

步骤2，对于第k次迭代，计算当前模型m_k对应的正演数据F[m_k]、雅克比矩阵J_k以及数据拟合差X²；

步骤3，搜索使数据拟合差X²最小时的模型m_k+1，同时计算模型m_k+1对应的粗糙度矩阵R_k+1以及步长a_k；

步骤4，判断步骤3得到的最小数据拟合差是否小于反演所要求达到的拟合差若是，则进入步骤5，否则，判断步骤3得到的最小数据拟合差是否小于第k-1次迭代得到的最小数据拟合差；若是，则进入步骤6，否则使步长a_k减半并返回步骤3；

步骤5，确定当数据拟合差X²等于反演所要求达到的拟合差时的拉格朗日乘子及对应的模型m′_k+1，根据模型m′_k+1计算粗糙度矩阵R′_k+1，判断R′_k+1是否大于R_k+1，若是，则使步长a_k减半并返回步骤3，否则进入步骤6；

步骤6，保存步长a_k和模型m_k+1；

步骤7，判断是否达到最大迭代次数或是否达到反演所要求达到的拟合差若达到最大迭代次数或达到反演所要求达到的拟合差/>则迭代结束，输出反演结果，否则，令k＝k+1并返回步骤2。

作为本发明的一种优选方案，步骤1所述根据初始模型m₁构造粗糙度矩阵为：

其中，R₁表示根据初始模型m₁构造的粗糙度矩阵，为N×N的矩阵，定义如下：

N为地层模型分段得到的总层数。

作为本发明的一种优选方案，所述步骤2的具体过程如下：

对于第k次迭代，向圆锥型场源通入电流，对于圆锥型场源的第i匝圆形发射线圈，其中心处垂直方向磁场的频率域响应为：

式中，H(ω)表示频率域响应，r_i为第i匝圆形发射线圈的半径，I(ω)为发射电流的傅式谱，λ为积分变量，h_i为第i匝圆形发射线圈离地高度，u₀为真空中的磁导率，Z₀为表层介质的输入波阻抗，Z⁽¹⁾为第一层的输入阻抗，J₁为第一类零阶贝塞尔函数；Z¹＝-iωμ₀/u₁，k₁为第一层介质的波数，σ₁为第一层介质的电导率，ω为角频率，μ₀为真空中磁导率；

对频率域响应公式进行140点汉克尔变换得到：

式中，H_z,i(ω)表示140点汉克尔变换得到的频率域，m为抽样点的位置，b为回线半径，a＝-7.91001919000e+00，s＝8.79671439570e-02，为汉克尔变换滤波系数；

设定激励电流I(t)为单位阶跃电流：

将140点汉克尔变换得到的频率域，通过傅里叶变换求得时间域中响应：

式中，t为时间，f(t)为时间域电磁场分布，F(ω)为频率域的电磁场分布，j为虚数单位；

将频率域响应公式代入上式即得时间域的磁场响应：

式中，Im表示虚部，Re表示实部，为角频率，H_z(t)为时间域电磁场，I₀为通入回线的电流，B_z(t)为时间域磁感应强度；

根据法拉第电磁感应定律，得到二次场的感应电位为v_i，将圆锥型场源产生的瞬变电磁感应电位视为n个线圈感应电位的叠加，则圆锥型场源下层状模型的瞬变电磁感应电位V(t)为：

式中，n为圆锥型场源缠绕的圆形发射线圈匝数，Δt_i是第i匝线圈与第1匝线圈距离产生的延时，对于第1匝线圈，Δt₁＝0；采用全区视电阻率算法将感应电位V(t)转换成反映电性变化的的视电阻率，即当前模型对应的正演数据；

雅克比矩阵J_k为：

其中，表示求F[m_k]的梯度，F[m_k]表示当前模型m_k对应的正演数据；

数据拟合差X²为：

X²＝‖Wd-WF[m_k]‖²

其中，W＝diag{1/σ₁,1/σ₂,…,1/σ_M}，σ_j为第j个数据的标准差，j＝1,…,M，M为场源装置测得的数据个数，d为场源装置测得的数据集合。

作为本发明的一种优选方案，步骤3所述步长a_k的计算公式为：

其中，表示拉格朗日乘子函数的模型集，m_k、m_k+1分别第k、k+1次迭代的模型，/>为第k次迭代的拉格朗日乘子。

作为本发明的一种优选方案，所述模型m_k+1定义为：

其中，为第k+1次迭代的拉格朗日乘子，/>为N×N的矩阵，N为地层模型分段得到的总层数，W＝diag{1/σ₁,1/σ₂,…,1/σ_M}，σ_j为第j个数据的标准差，j＝1,…,M，M为场源装置测得的数据个数，J_k为雅克比矩阵，d_k＝d-F[m_k]+J_km_k，d为场源装置测得的数据集合，F[m_k]表示当前模型m_k对应的正演数据。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

1、本发明的反演方法不需要依赖于初始模型，对初始模型没有特别要求，且计算效率高，较好地反映了模型的地质结构，对模型具有平滑效果。

2、本发明采用OCCAM法得到的反演电阻率更接近各层介质的真实电阻率，对电阻率的改进效果良好。

3、本发明反演采用简化后的拉格朗日乘子算法使得反演仍稳定、速度更快、结果更真实。

4、本发明在一定步长下逐次递减拉格朗日乘子，使反演仍稳定，速度更快，结果更真实。

附图说明

图1是一种圆锥型场源装置的结构图。

图2是本发明一种圆锥型场源瞬变电磁反演方法的流程图。

图3是D型模型反演曲线。

图4是G型模型反演曲线。

图5是H型模型反演曲线。

图6是K型模型反演曲线。

图7是KH型模型反演曲线。

具体实施方式

下面详细描述本发明的实施方式，所述实施方式的示例在附图中示出。下面通过参考附图描述的实施方式是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。

本发明基于OCCAM反演算法，提出了一种圆锥型场源瞬变电磁反演方法，具体流程如图2所示。OCCAM法是一种带平滑约束的最小二乘反演，受初始模型影响小，运算稳定收敛，是一种有效的数据反演处理方法。

首先，定义粗糙度，代表了模型的光滑程度。粗糙度可表示为模型相对某一坐标的一阶或二阶导数的平方的积分，如对z方向，则：

式中，z为深度，m为模型相应深度下所对应的电性参数，一般是地层的电阻率，或者电阻率的对数。模型的最底部为均匀半空间，将地层模型连续分段处理得到：

m_(z)＝m_i,z_i-1<z<z_i+1,i＝1,2,…,N (2)

式中，z₀＝0，z_i-1/z_i取为小于1的常数，m_i是电阻率或电导率值。

将粗糙度的微分形式离散化，得到：

粗糙度用矩阵表示为：

式中是N×N的矩阵，定义如下：

根据加权最小二乘法原理，定义模型数据与实测数据的拟合差X为：

式中，M为观测数据个数，d_j是第j个观测数据，σ_j为第j个数据的标准差，F_j[m]为模型m对应的正演数据。

正演部分如下：

对于半径r_i，通入电流I的第i匝圆形发射线圈，其中心处垂直方向磁场的频率域响应为：

式中，λ为积分变量，I(w)为发射电流的傅式谱，h_i为第i匝线圈离地高度，Z⁽¹⁾为第一层的输入阻抗，J₁为第一类零阶贝塞尔函数。对公式(7)进行汉克尔变换得到：

式中，I为回线源中通入的电流，m为抽样点的位置，在140点J₁汉克尔变换的滤波器系数中，a＝-7.91001919000e+00，s＝8.79671439570e-02。然后进行频率域向时间域的转换。假设激励电流为单位阶跃电流：

公式(8)通过傅里叶变换求得时域中响应：

将公式(7)代入上式即可得到时间域的磁场响应：

式中Im为虚部，Re为实部，为角频率，t为时间。式(11)、(12)含有特殊函数，内层贝塞尔函数积分采用140点快速汉克尔变换计算。根据法拉第电磁感应定律，可得到二次场的感应电位为v_i。圆锥形场源产生的瞬变电磁感应电位可视为多个线圈感应电位的叠加，因此，得到该场源下层状模型的瞬变电磁感应电位为：

式中,Δt_i是是各匝线圈与第1匝线圈距离不同而产生的延时，对于第1匝线圈i＝1，Δt₁＝0。进而，采用全区视电阻率算法将感应电位V(t)转换成反映电性变化的的视电阻率ρ_a(t)。

这样求解的数学问题就转化为：对于给定的数据集d及其误差，寻求当X²达到可接受的值时R₁或R₂尽可能小的模型m。

数据拟合差写为：

X²＝‖Wd-WF[m]‖² (14)

式中,W＝diag{1/σ₁,1/σ₂,…,1/σ_M}。

引入拉格朗日乘子μ^-1，构造无约束的目标函数：

式中，R₁为粗糙度矩阵，为反演所要求达到的拟合差。

在反演迭代过程中，为使总目标函数最小，求取目标函数的梯度，并令这时模型向量m满足：

式中，J是M×N阶的雅克比矩阵：

在迭代过程开始前，给定一个初始模型m₁，如果F在m₁可微，则对于足够小的向量Δ有：

F[m₁+Δ]＝F[m₁]+J₁Δ (18)式中，Δ＝m₂-m₁，将式(18)代入式(15)，则在m₂会出现如下的线性问题：

上式右端第二项圆括号中的表达式是一种数据向量，记作d₁。

进而，定义模型m₂，使U取得最小值，有

选取μ以达到期望的拟合差，在迭代中m₂由系列m₁,m₂,…中后面的模型依次取代，系列中每一个前面的向量都是下个向量的初始近似值。

当第k次迭代完成后，定义向量：

用一系列的μ值计算模型m_k+1(μ)的真正拟合差：

X_k+1(μ)＝‖Wd-WF[m_k+1(μ)]‖ (22)

为防止拟合问题发散，定义拉格朗日乘子函数的新模型集，即：

G(μ)＝(1-a)m_k+am_k+1(μ) (23)

其中，a为步长，其将被连续地平分直到获得一个合适的模型。

层状模型更能检验反演算法的有效性，尤其能反映出算法对层界面的分辨能力。分别建立典型的2层(D型和G型)、3层(H型、K型)以及4层(KH型)介质模型来验证。图3、图4、图5、图6和图7分别为层状模型的反演结果。图中左上部分是预设的初始模型及各项参数，左下为模型的一维正演结果，横坐标为时间，纵坐标为电阻率；右边为反演结果，横坐标为视电阻率，纵坐标为深度。图中反演曲线的形态反映出对应介质模型的特征，且对层界面的分辨较好。

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。

Claims (2)

1.一种圆锥型场源瞬变电磁反演方法，其特征在于，该方法包括如下步骤：

步骤1，读入初始模型m₁即初始电阻率，根据初始模型m₁构造粗糙度矩阵，设置初始拉格朗日乘子

步骤2，对于第k次迭代，计算当前模型m_k对应的正演数据F[m_k]、雅克比矩阵J_k以及数据拟合差X²；

具体过程如下：

对于第k次迭代，向圆锥型场源通入电流，对于圆锥型场源的第i匝圆形发射线圈，其中心处垂直方向磁场的频率域响应为：

式中，H(ω)表示频率域响应，r_i为第i匝圆形发射线圈的半径，I(ω)为发射电流的傅式谱，λ为积分变量，h_i为第i匝圆形发射线圈离地高度，u₀为真空中的磁导率，Z₀为表层介质的输入波阻抗，Z⁽¹⁾为第一层的输入阻抗，J₁为第一类零阶贝塞尔函数；k₁为第一层介质的波数，σ₁为第一层介质的电导率，ω为角频率，μ₀为真空中磁导率；

对频率域响应公式进行140点汉克尔变换得到：

式中，H_z,i(ω)表示140点汉克尔变换得到的频率域，m为抽样点的位置，b为回线半径，a＝-7.91001919000e+00，s＝8.79671439570e-02，为汉克尔变换滤波系数；

设定激励电流I(t)为单位阶跃电流：

将140点汉克尔变换得到的频率域，通过傅里叶变换求得时间域中响应：

式中，t为时间，f(t)为时间域电磁场分布，F(ω)为频率域的电磁场分布，j为虚数单位；

将频率域响应公式代入上式即得时间域的磁场响应：

式中，Im表示虚部，Re表示实部，为角频率，H_z(t)为时间域电磁场，I₀为通入回线的电流，B_z(t)为时间域磁感应强度；

根据法拉第电磁感应定律，得到二次场的感应电位为v_i，将圆锥型场源产生的瞬变电磁感应电位视为n个线圈感应电位的叠加，则圆锥型场源下层状模型的瞬变电磁感应电位V(t)为：

式中，n为圆锥型场源缠绕的圆形发射线圈匝数，Δt_i是第i匝线圈与第1匝线圈距离产生的延时，对于第1匝线圈，Δt₁＝0；采用全区视电阻率算法将感应电位V(t)转换成反映电性变化的视电阻率，即当前模型对应的正演数据；

雅克比矩阵J_k为：

其中，表示求F[m_k]的梯度，F[m_k]表示当前模型m_k对应的正演数据；

数据拟合差X²为：

X²＝‖Wd-WF[m_k]‖²

其中，W＝diag{1/δ₁,1/δ₂,…,1/δ_M}，δ_c为第c个数据的标准差，c＝1,…,M，M为场源装置测得的数据个数，d为场源装置测得的数据集合；

步骤3，搜索使数据拟合差X²最小时的模型m_k+1，同时计算模型m_k+1对应的粗糙度矩阵R_k+1以及步长a_k；

所述模型m_k+1定义为：

其中，为第k+1次迭代的拉格朗日乘子，θ为N×N的矩阵，N为地层模型分段得到的总层数，J_k为雅克比矩阵，d_k＝d-F[m_k]+J_km_k，F[m_k]表示当前模型m_k对应的正演数据；

所述步长a_k的计算公式为：

其中，表示拉格朗日乘子函数的模型集，m_k、m_k+1分别第k、k+1次迭代的模型，/>为第k次迭代的拉格朗日乘子；

步骤4，判断步骤3得到的最小数据拟合差是否小于反演所要求达到的拟合差若是，则进入步骤5，否则，判断步骤3得到的最小数据拟合差是否小于第k-1次迭代得到的最小数据拟合差；若是，则进入步骤6，否则使步长a_k减半并返回步骤3；

步骤5，确定当数据拟合差X²等于反演所要求达到的拟合差时的拉格朗日乘子/>及对应的模型m′_k+1，根据模型m′_k+1计算粗糙度矩阵R′_k+1，判断R′_k+1是否大于R_k+1，若是，则使步长a_k减半并返回步骤3，否则进入步骤6；

步骤6，保存步长a_k和模型m_k+1；

步骤7，判断是否达到最大迭代次数或是否达到反演所要求达到的拟合差若达到最大迭代次数或达到反演所要求达到的拟合差/>则迭代结束，输出反演结果，否则，令k＝k+1并返回步骤2。

2.根据权利要求1所述圆锥型场源瞬变电磁反演方法，其特征在于，步骤1所述根据初始模型m₁构造粗糙度矩阵为：

R₁＝‖θm₁‖²

其中，R₁表示根据初始模型m₁构造的粗糙度矩阵，θ为N×N的矩阵，定义如下：

N为地层模型分段得到的总层数。