CN112882114B - 一种面向复杂地质的瞬变电磁自适应横向约束反演方法 - Google Patents

Abstract

本发明公开一种面向复杂地质的瞬变电磁自适应横向约束反演方法，该方法包括：对瞬变电磁测线包含的N_s个测点数据进行虚拟波场变换，计算并记录瞬变电磁相邻测点波场变换后数据的二范数；对得到的N_s‑1个二范数进行求倒数，进行归一化处理，并乘以缩放系数后，进行横向约束自适应加权矩阵的计算；将横向约束自适应加权矩阵应用于电阻率和厚度的横向约束、深度约束中，进行自适应横向约束反演，得到瞬变电磁自适应横向约束反演结果。能够在数据反演时自动识别数据连续性较好的测点及电性差异大的断面测点，对连续性好的测点施加正常横向约束，对连续性差的测点不施加或施加弱横向约束，在突出断面异常的同时，有效保留瞬变电磁数据本身良好的横向分辨率。

Description

一种面向复杂地质的瞬变电磁自适应横向约束反演方法

技术领域

本发明涉及地球物理信号处理与分析技术研究领域，特别是一种面向复杂地质的瞬变电磁自适应横向约束反演方法。

背景技术

瞬变电磁单点一维反演得到的电阻率和厚度经常会出现横向不连续性，给后续的地球物理解释工作带来困难，为此针对连续性较好的地质情况，如沉积岩地区，学者们提出了瞬变电磁横向约束反演方法，该方法通过在相邻测点间施加横向约束，保证了电阻率及厚度的横向连续性，是一种基于一维模型的拟二维反演技术。然而，对于一些复杂地质情况，如城市地质、岩溶地质，这种地质环境的特点是既有连续性较好的区域，又存在连续性较差区域，直接用横向约束反演方法会弱化这种差异，造成断面地带数据的错误解释，单纯采用单点反演一维反演会使得层界面不光滑。此外，在缺乏地质先验信息的探测区域，地下结构的整体连续性及断层分布情况未知，便难以选择最佳反演解释方案。因此，需要研究适用于各种地质情况下的自适应横向约束反演，能够在数据反演时自动识别数据连续性较好的测点及电性差异大的断面测点，对连续性好的测点施加正常横向约束，对连续性差的测点不施加或施加弱横向约束，在突出断面异常的同时，有效保留瞬变电磁数据本身良好的横向分辨率，对实现复杂地质情况下瞬变电磁数据的高精度解释具有重要意义。

发明内容

本发明的目的就是针对上述现有技术的不足，提供一种面向复杂地质的瞬变电磁自适应横向约束反演方法。

本发明的目的是通过以下技术方案实现的：

一种面向复杂地质的瞬变电磁自适应横向约束反演方法，该方法包括以下步骤：

S1、对瞬变电磁测线包含的N_s个测点数据进行虚拟波场变换，计算并记录瞬变电磁相邻测点波场变换后数据的二范数；

S2、对得到的N_s-1个二范数进行求倒数，进行归一化处理，并乘以缩放系数后，进行横向约束自适应加权矩阵的计算；

S3、将横向约束自适应加权矩阵应用于电阻率和厚度的横向约束、深度约束中，进行自适应横向约束反演，得到的反演结果即为瞬变电磁自适应横向约束反演结果。

所述的步骤S1中的虚拟波场变换，是指将采集到的瞬变电磁数据利用虚拟波场变换公式变换至波场域，虚拟波场变换公式的计算包括以下步骤：

T11、对经典瞬变电磁波场变换公式求导，反映波动场与扩散场关系的经典波场变换公式为：

其中B(t)为扩散场量，指的是电场或磁场分量，t为扩散场时间，q为虚拟域时间，U表示虚拟波动场量；对经典波场变换公式(1)进行求导变换可得到感应电压与波动场之间的转换关系式为：

T12、对转换方程(2)进行数值化离散将其转化为线性方程形式，即：

其中A为积分系数矩阵，积分系数矩阵中的元素计算满足：

其中，q域的选取采用对数进行离散，范围大于瞬变电磁探测时间范围。

T13、对公式(3)进行反变换求解虚拟波场，采用正则化算法寻找虚拟波动场量U的最优近似解：

其中，D是n阶差商算子；α是正则化因子，控制约束的权重；φ是正则解,令公式(5)取值最小，完成虚拟波动场量U的求解。

所述的步骤S2中的横向约束自适应加权矩阵包括电阻率和厚度约束及深度约束，具体包括：

T21、对步骤S1计算得到的相邻测点波场变换后的数据二范数进行求倒数：

ψ_i＝1/||U_obsi-U_obs(i+1)||₂，i＝1,2,...,N_s-1 (6)

其中U_obsi为第i个测点数据的虚拟波场变换数据，ψ_i为第i测点与第i+1测点的二范数的倒数，它表示了对相邻测点的差异性取倒数；

T22、求解相邻测点间的相关性系数rel_i，满足计算公式(7)：

其中ψ_max为ψ_i中的最大值；λ为缩放系数，为(0,1)之间的一个固定值；rel_i表示了第i测点与第i+1测点的相关度，其取值越接近1，相邻两测点之间的相关性越高；

T23、进行电阻率和厚度横约束自适应加权矩阵计算：

R_p'Δm＝Δr_p+e_rp (8)

R'_p＝W_p.*R_p，Δr_p＝-R'_pm₀ (9)

其中R_p'为电阻率和厚度横向约束自适应加权矩阵；m为全部反演模型参数；m₀为初始反演模型参数，Δm＝m-m₀，为模型修正量，e_rp表示同一测线上相邻测点间同一层模型参数的差异，R_p为横向约束矩阵，W_p为电阻率与厚度约束自适应相关加权因子，分别满足：

其中，

式中，S＝N_m*(N_s-1)，T＝N_m*N_s，N_m为每个测点的模型参数个数，N_s为测点个数,W_p为与相邻测点数据二范数差异性相关的加权矩阵，加权因子越大，横向约束越强，反演得到的参数横向连续性越好，反演结果越光滑。

T24、进行深度约束自适应加权矩阵计算：

R_d'Δm＝Δr_d+e_rd (13)

R'_d＝W_d·*R_d，Δr_d＝-R'_dm₀ (14)

式中，e_rd为测线上相邻测点之间层界面深度的差异；R_d为由厚度和深度表示得稀疏矩阵，W_d为电阻率与厚度约束自适应相关加权因子，分别满足：

其中：

式中，

M＝(N_s-1)×(N-1)；h_i,j为第i个测点第j层的厚度；d_i,j为第i个测点第j层下界面的深度。至此，完成横向约束自适应加权矩阵的计算。

所述的步骤S3中的自适应横向约束反演，是将数据拟合方程式、电阻率和厚度约束方程式及深度约束方程式进行总体集成，具体包括如下步骤：

T31、可得在横向约束条件下的反演方程：

其中，J为总体雅可比矩阵：

J_i为第i个测点的雅可比矩阵，将式(18)简化为式(20)：

A·Δm＝b+e (20)

T32、利用奇异值分解法求解式(20)，将式(20)两边乘以矩阵A的转置，得到方程(21)：

A^TA·Δm＝A^Tb (21)

求解方程(21)为

Δm＝(A^TA)^-1A^Tb (22)

为避免A^TA过小导致反演解出现不稳定性，在系数矩阵中引入阻尼因子λ²，将式(22)变为

Δm＝(A^TA+λ²I)^-1A^T·b (23)

式中，I为单位矩阵，式(23)可使用通用奇异值分解方法求解，在获得模型改正量后，可以进行模型更新。如此循环直到数据拟合差小于事先设定的阈值或者当迭代次数大于设定的迭代次数为止，即可输出反演结果。

有益效果：本发明提出了一种面向复杂地质的瞬变电磁自适应横向约束反演方法。首先将采集到瞬变电磁数据进行虚拟波场变换，提高对电性界面差异的识别能力，利用这种差异计算相邻测点的相关加权因子，根据相关加权因子生成横向约束自适应加权矩阵，进行自适应横向约束反演，能够在数据反演时自动识别数据连续性较好的测点及电性差异大的断面测点，对连续性好的测点施加正常横向约束，对连续性差的测点不施加或施加弱横向约束，在突出断面异常的同时，有效保留瞬变电磁数据本身良好的横向分辨率，具有较好的实际应用价值。

附图说明

图1为两组不同半空间电阻率下瞬变电磁正演结果；

图2为对瞬变电磁正演数据进行虚拟波场变换后结果；

图3为含有不连续性地质下瞬变电磁理论模型；

图4为对理论模型进行单点一维反演结果；

图5为对理论模型进行常规横向约束反演结果；

图6为瞬变电磁自适应横向约束反演结果。

具体实施方式

下面结合附图和实施例对本发明做进一步的详细说明：

一种面向复杂地质的瞬变电磁自适应横向约束反演方法的详细说明：

图1展示了在层状大地结构下，在背景电阻率为20欧姆每米的介质中，在地下10米处分别含有一个层厚度为10米，电阻率为50欧姆每米、200欧姆每米的高阻介质时瞬变电磁感应电压数据正演计算结果，可以看出两者曲线基本一致，仅有细微差别，因此难以直接看出电性界面存在明显差异；将采集到的瞬变电磁数据利用虚拟波场变换公式变换至波场域，能够提高数据对电性界面差异的识别能力，虚拟波场变换公式的计算包括以下步骤：

T11、对经典瞬变电磁波场变换公式求导，反映波动场与扩散场关系的经典波场变换公式为：

其中E为扩散场量(电场或磁场分量)，t，q分别为扩散场时间和虚拟域时间，U表示虚拟波动场量；由于常规瞬变电磁仪器采集的是衰减的感应电压数据，通过对经典波场变换公式进行求导变换可得到感应电压与波动场之间的转换关系式为：

T12、对转换方程进行数值化离散将其转化为线性方程形式，即：

其中A为积分系数矩阵，其计算满足：

其中，q域的选取采用对数进行离散，其范围需要大于瞬变电磁探测时间范围，本实施例中瞬变电磁探测时间范围为[10^-6～10^-3s],q域的选取范围为[10^-6～10^-2s]，采用对数等间隔进行离散，离散点数为50个。

T13、求解虚拟波场需要对式(3)进行反变换，考虑到方程的病态性，采用正则化算法寻找虚拟波动场量U的最优近似解：

其中，D是n阶差商算子；α是正则化因子，控制约束的权重，本实施例利用L曲线法进行求解；φ是正则解。令上式取值最小，即可完成虚拟波动场量U的求解。图2展示了对图1瞬变电磁感应电压数据进行虚拟波场转换后的图像，可以看出此时两者间电阻率的差异得到了明显放大；

图3展示了采用瞬变电磁正演得到的地下含断面带的理论模型，在地下5米与20处存在低阻异常带，但在横向上呈分段连续分布。在后续反演中，对理论模型数据加入了原始信号幅值5％的随机噪声，图4展示了采用传统SVD反演成像结果，可以看出反演结果横向连续性较差，层界面不光滑；图5展示了采用传统横向约束反演的成像结果，可以看出传统横向约束反演虽然使得层界面光滑，但整体反演结果会弱化断面带的存在，使得在断面处成像结果相对失真；下面进行自适应横向约束反演。

对图3中的理论测线模型数据进行虚拟波场变换后，求解其横向约束自适应加权矩阵，包括电阻率和厚度约束及深度约束矩阵，具体包括以下计算步骤：

T21、对步骤S1计算得到的相邻测点波场变换后的数据二范数进行求倒数：

ψ_i＝1/||U_obsi-U_obs(i+1)||₂，i＝1,2,...,N_s-1 (6)

其中U_obsi为第i个测点数据的虚拟波场变换数据，ψ_i为第i测点与第i+1测点的二范数的倒数，它表示了对相邻测点的差异性取倒数；本实施例中，理论模型包含了11个测点数据，即N_s＝11。

T22、求解相邻测点间的相关性系数rel_i，满足计算公式为：

其中ψ_max为ψ_i中的最大值；λ为缩放系数，为(0,1)之间的一个固定值；rel_i表示了第i测点与第i+1测点的相关度，其取值越接近1，相邻两测点之间的相关性越高，本实施例中λ＝0.3；

T23、进行电阻率和厚度横约束自适应加权矩阵计算：

R_p'Δm＝Δr_p+e_rp (8)

R'_p＝W_p.*R_p，Δr_p＝-R'_pm₀ (9)

其中R_p'为电阻率和厚度横向约束自适应加权矩阵；m为全部反演模型参数；m₀为初始反演模型参数，Δm＝m-m₀，为模型修正量，e_rp表示同一测线上相邻测点间同一层模型参数的差异，R_p为横向约束矩阵，W_p为电阻率与厚度约束自适应相关加权因子，分别满足：

其中，

式中，S＝N_m*(N_s-1)，T＝N_m*N_s，N_m为每个测点的模型参数个数，N_s为测点个数,W_p为与相邻测点数据二范数差异性相关的加权矩阵，加权因子越大，横向约束越强，反演得到的参数横向连续性越好，反演结果越光滑，本实施中采用三层模型进行反演，即N_m＝5；

T24、进行深度约束自适应加权矩阵计算：

R_d'Δm＝Δr_d+e_rd (13)

R'_d＝W_d·*R_d，Δr_d＝-R'_dm₀ (14)

式中，e_rd为测线上相邻测点之间层界面深度的差异；R_d为由厚度和深度表示得稀疏矩阵，W_d为电阻率与厚度约束自适应相关加权因子，分别满足：

式中，

M＝(N_s-1)×(N-1)；h_i,j为第i个测点第j层的厚度；d_i,j为第i个测点第j层下界面的深度，本实施例中采用三层模型进行反演，即N＝3。至此，完成横向约束自适应加权矩阵的计算。

进行自适应横向约束反演，是将数据拟合方程式、电阻率和厚度约束方程式及深度约束方程式进行总体集成，具体包括如下步骤：

T31、可得在横向约束条件下的反演方程：

其中，J为总体雅可比矩阵：

J＝diag(J₁,J₂,...,J_Ns) (19)

J_i为第i个测点的雅可比矩阵。式(18)可进一步简化为

A·Δm＝b+e (20)

T2、利用奇异值分解法求解式(20)。将式(20)两边乘以矩阵A的转置，可得到如下方程：

A^TA·Δm＝A^Tb (21)

其解为

Δm＝(A^TA)^-1A^Tb (22)

为避免A^TA过小导致反演解出现不稳定性，在系数矩阵中引入阻尼因子λ²，将式(22)变为

Δm＝(A^TA+λ²I)^-1A^T·b (23)

式中，I为单位矩阵，式(23)可使用通用奇异值分解方法求解，在获得模型改正量后，可以进行模型更新。如此循环直到数据拟合差小于事先设定的阈值或者当迭代次数大于设定的迭代次数为止，即可输出反演结果。本实施例中，最大迭代次数设置为50次，阻尼因子λ²采用将[10⁰～10^-5]区间对数离散为50份的值，在逐次迭代中不断缩小。

图6展示了对理论模型进行自适应横向约束反演的成像结果，相比于传统SVD反演及传统横向约束反演，自适应横向约束反演成像结果更接近真实理论模型，层界面光滑，且能准备呈现断面带信息，具有较好的应用价值。

Claims (3)

1.一种面向复杂地质的瞬变电磁自适应横向约束反演方法，其特征在于，该方法包括以下步骤：

S1、对瞬变电磁测线包含的N_s个测点数据进行虚拟波场变换，计算并记录瞬变电磁相邻测点波场变换后数据的二范数；

S2、对得到的N_s-1个二范数进行求倒数，进行归一化处理，并乘以缩放系数后，进行横向约束自适应加权矩阵的计算；

S3、将横向约束自适应加权矩阵应用于电阻率和厚度的横向约束、深度约束中，进行自适应横向约束反演，得到的反演结果即为瞬变电磁自适应横向约束反演结果；

其中，所述步骤S2中的横向约束自适应加权矩阵包括电阻率和厚度约束及深度约束，具体包括：

T21、对步骤S1计算得到的相邻测点波场变换后的数据二范数进行求倒数：

ψ_i＝1/||U_obsi-U_obs(i+1)||₂，i＝1,2,...,N_s-1

其中U_obsi为第i个测点数据的虚拟波场变换数据，ψ_i为第i测点与第i+1测点的二范数的倒数，它表示了对相邻测点的差异性取倒数；

T22、求解相邻测点间的相关性系数rel_i，满足：

其中ψ_max为ψ_i中的最大值；λ为缩放系数，为(0,1)之间的一个固定值；rel_i表示了第i测点与第i+1测点的相关度，其取值越接近1，相邻两测点之间的相关性越高；

T23、进行电阻率和厚度横约束自适应加权矩阵计算：

R_p'Δm＝Δr_p+e_rp

R'_p＝W_p.*R_p，Δr_p＝-R'_pm₀

其中R_p'为电阻率和厚度横向约束自适应加权矩阵；m为全部反演模型参数；m₀为初始反演模型参数，Δm＝m-m₀，为模型修正量，e_rp表示同一测线上相邻测点间同一层模型参数的差异，R_p为横向约束矩阵，W_p为电阻率与厚度约束自适应相关加权因子，分别满足：

其中，

式中，S＝N_m*(N_s-1)，T＝N_m*N_s，N_m为每个测点的模型参数个数，N_s为测点个数,W_p为与相邻测点数据二范数差异性相关的加权矩阵，加权因子越大，横向约束越强，反演得到的参数横向连续性越好，反演结果越光滑；

T24、进行深度约束自适应加权矩阵计算：

R_d'Δm＝Δr_d+e_rd

R'_d＝W_d·*R_d，Δr_d＝-R'_dm₀

式中，e_rd为测线上相邻测点之间层界面深度的差异；R_d为由厚度和深度表示得稀疏矩阵，W_d为电阻率与厚度约束自适应相关加权因子，分别满足：

其中：

式中，

M＝(N_s-1)×(N-1)；h_i,j为第i个测点第j层的厚度；d_i,j为第i个测点第j层下界面的深度；至此，完成横向约束自适应加权矩阵的计算。

2.按照权利要求1所述的方法，其特征在于，

所述步骤S1中的虚拟波场变换，是指将采集到的瞬变电磁数据利用虚拟波场变换公式变换至波场域，虚拟波场变换公式的计算包括以下步骤：

T11、对经典瞬变电磁波场变换公式求导，反映波动场与扩散场关系的经典波场变换公式为：

其中B(t)为扩散场量，指的是电场或磁场分量，t为扩散场时间，q为虚拟域时间，U表示虚拟波动场量；对经典波场变换公式(1)进行求导变换可得到感应电压与波动场之间的转换关系式为：

T12、对转换方程(2)进行数值化离散将其转化为线性方程形式，即：

其中A为积分系数矩阵，积分系数矩阵中的元素计算满足：

其中，q域的选取采用对数进行离散，范围大于瞬变电磁探测时间范围；

T13、对公式(3)进行反变换求解虚拟波场，采用正则化算法寻找虚拟波动场量U的最优近似解：

其中，D是n阶差商算子；α是正则化因子，控制约束的权重；φ是正则解,令公式(5)取值最小，完成虚拟波动场量U的求解。

3.按照权利要求1所述的方法，其特征在于，所述的步骤S3中的自适应横向约束反演，是将数据拟合方程式、电阻率和厚度约束方程式及深度约束方程式进行总体集成，具体包括如下步骤：

T31、可得在横向约束条件下的反演方程：

其中，J为总体雅可比矩阵：

J_i为第i个测点的雅可比矩阵，将反演方程简化为：

A·Δm＝b+e

T32、利用奇异值分解法求解简化后的反演方程，将简化后的反演方程两边乘以矩阵A的转置，得到方程：

A^TA·Δm＝A^Tb

求解方程为

Δm＝(A^TA)^-1A^Tb

为避免A^TA过小导致反演解出现不稳定性，在系数矩阵中引入阻尼因子λ²，变为

Δm＝(A^TA+λ²I)^-1A^T·b

式中，I为单位矩阵，使用通用奇异值分解方法求解，在获得模型改正量后，进行模型更新，如此循环直到数据拟合差小于事先设定的阈值或者当迭代次数大于设定的迭代次数为止，即可输出反演结果。

Images