CN110908000B - 基于变维贝叶斯的隧道瞬变电磁数据解释方法 - Google Patents

Abstract

本发明涉及一种基于变维贝叶斯的隧道瞬变电磁数据解释方法，该方法包括使用瞬变电磁仪器，放置于隧道掌子面处进行探测，获取隧道掌子面含水构造的观测数据；对获得的观测数据进行瞬变电磁数据预处理；使用变维贝叶斯算法对预处理后的数据反演，获得掌子面前方所有可能的地质电阻率分布及层界面位置信息，提取其中最大概率对应的电阻率及层界面位置分布信息；根据最大概率地质电阻率分布及层界面位置信息，给出隧道掌子面前方含水构造异常情况，并进行不确定度分析，得出准确的含水构造分布范围。在反演解释隧道掌子面前方低阻异常(含水构造)位置的同时，给出反演结果的不确定度分析，最后给出电阻率分布范围，为隧道工程安全开发提供了预警指导。

Description

基于变维贝叶斯的隧道瞬变电磁数据解释方法

技术领域

本发明涉及地球物理信号处理与分析技术研究领域，特别是一种基于变维贝叶斯的隧道瞬变电磁数据解释方法。

背景技术

隧道工程在施工过程中，常面临地下突涌水、塌方等地质灾害，这些灾害不仅会导致施工延误，更会造成隧道坍塌威胁施工人员的生命安全。瞬变电磁法由于具有探测便捷，对地下水体响应灵敏等特征，近年来越来越多被应用在隧道灾害超前地质预测中。

隧道瞬变电磁数据的正确解释十分重要，其决定着施工人员对前方地质情况的正确判断，瞬变电磁法的反演是一个非线性问题，其解通常也是不唯一的，即对于地下模型参数的选择，存在多个或者无穷个地电参数模型能够与实测数据所吻合，而常规瞬变电磁反演方法均只能够给出单个最优解，并不能获得当前反演模型参数的不确定度信息。通常而言，估计参数不确定性或模型非唯一性与估计参数值本身一样重要，尤其对于隧道瞬变电磁法而言，这种缺陷会造成较大的安全威胁。因此，对隧道瞬变电磁数据进行反演解释的同时评估反演结果的不确定性，掌握隧道掌子面前方地质可能性电阻率分布范围，从而准确预测突水层的位置，对隧道工程安全开发具有重要意义。

发明内容

本发明的目的就是针对上述现有技术的不足，提供基于变维贝叶斯的隧道瞬变电磁数据解释方法，在反演解释隧道掌子面前方低阻异常(含水构造)位置的同时，给出反演结果的不确定度分析，最后给出电阻率分布范围，为隧道工程安全开发提供了预警指导。

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

一种基于变维贝叶斯的隧道瞬变电磁数据解释方法，该方法包括以下步骤：

S1、使用瞬变电磁仪器，放置于隧道掌子面处进行探测，获取隧道掌子面含水构造的观测数据；

S2、对获得的观测数据进行瞬变电磁数据预处理；

S3、使用变维贝叶斯算法对预处理后的数据反演，获得掌子面前方所有可能的地质电阻率分布及层界面位置信息，提取其中最大概率对应的电阻率及层界面位置分布信息；

S4、根据最大概率地质电阻率分布及层界面位置信息，给出隧道掌子面前方含水构造异常情况，并进行不确定度分析，得出准确的含水构造分布范围。

进一步地，所述步骤S1中的瞬变电磁仪器是适用于隧道探测的小型瞬变电磁系统，采用重叠回线配置，发射系统贴靠掌子面发射双极性梯形波，接收系统同步采集观测信号数据。

进一步地，所述的步骤S1中的隧道掌子面含水构造的探测范围为掌子面前方0至80米以内。

进一步地，所述步骤S2中瞬变电磁数据预处理的具体步骤包括：

1)、将相同测点位置的瞬变电磁观测数据进行叠加，叠加次数与该测点位置上瞬变电磁仪的发射次数一致；

2)、瞬变电磁发射波形为双极性梯形波，对于叠加后的数据将正负极性半周期发射关断后对应的数据相减；

3)、对完成步骤2)操作的测点数据做对数等间隔数据抽道。

进一步地，所述步骤S3中的变维贝叶斯算法的包括如下具体步骤；

①设置反演初始模型及约束参数

在反演初始需要设置模型参数的初始值与范围作为先验信息，包括：初始模型电阻率ρ、初始地层数k、初始层界面位置z、最小地层层数k_min和最大地层层数k_max、最小层界面位置z_min和最大层界面位置z_max、最小和最大电阻率值(ρ_min,ρ_max)，定义最小层厚度h_min为：

h_min＝(z_max-z_min)/(2·k_max)；

②根据当前模型参数定义后验概率分布:

其中，p(m)为先验概率；p(d_obs|m)为似然函数；d_obs为观测数据；m为模型参数；p(k)、p(z|k)和p(ρ|k,z)分别为层数、界面位置、电阻率的先验概率分布，p(d_obs|k,z,ρ)为表征数据拟合程度的似然函数，计算公式如下：

其中，Δz代表新界面的可用空间；

其中，N_d为数据个数；F(m)为正演模型响应；C_d为数据方差矩阵；

其中，ρ₀为全空间最佳匹配电阻率，通过最小二乘法进行数据拟合得到；

其中ρ_r表示设置的电阻率变化范围；

③计算建议分布产生候选模型参数

在变维贝叶斯反演中，候选模型是由当前模型的建议分布函数生成，与之前的反演模型无关，让贝叶斯反演受初始模型的影响较小，候选模型的建议分布由下式定义：

q(m'|m)＝q(k'|k)q(z'|k',z)q(ρ'|k',z',ρ)

理想的建议分布与后验分布相同，以便加快模型的收敛速度，采用可跳跃马尔科夫链蒙特卡洛方法进行采样，该采样过程包括以下四种基本状态：

a、新层生成，在满足不超过最大设置层数且大于最小层厚度的前提下，在最大和最小层界面位置间随机生成一个新的界面，层数加1；

b、旧层灭亡，在现有的层界面中，随机删除一个界面，层数减1；

c、扰动更新，在层界面数保持不变的条件下，随机选择以层界面，将其所在的位置在(-h_min,h_min)范围内扰动；

d、保持不变，层数与界面位置都不变，只对电阻率进行反演。

上述四种基本状态的概率满足：

相应层界面的位置按照均匀分布给出满足：

电阻率的建议分布定义为以当前模型电阻率为均值的多维正态分布，满足：

其中，ρ_k'与

分别为候选模型的电阻率及对应候选模型界面位置的当前电阻率，C_f为后验电阻率方差，可由当前电阻率的正演模拟值与观测数据的二范数来定义：C_f＝λ·||F(m)-d_obs||₂，其中，λ为缩放系数；

④根据后验概率分布与候选模型建议分步计算接受概率

候选模型产生后，按照接受概率判断是否接受候选模型，定义式子为：

其中，|J|为候选模型对当前模型导数构成的雅克比矩阵，取1，接受概率化简为：

在0到1间随机产生一个数，比较该数与α的大小，如果α较大，则候选模型被接受，否则被拒绝接受，继续按建议分布生成新的候选模型，直至模型收敛。所有被接受的模型共同构成反演结果。

进一步地，所述步骤④中的模型收敛通过设置最大模型采样数来判断，当模型采样数达到该设定值时停止搜索。

进一步地，所述步骤S4中不确定度的分析根据变维贝叶斯反演结果，反演结果由所有被接受的模型构成，包含层界面位置信息和电阻率信息，根据层界面位置信息得到层厚度，将层界面厚度信息、对应的电阻率信息分别存入矩阵，建立模型层界面厚度与模型电阻率的关系图，根据关系图得出层界面厚度与对应电阻率的分布范围。

有益效果：在对隧道瞬变电磁数据进行预处理后，运用变维贝叶斯算法对隧道数据进行反演解释，该反演从全局对所有可能符合的模型进行搜索，从概率学上得到反演结果分布。该方法可以保证反演结果的全局最优性，并能够对反演结果进行不确定性分析，从而掌握隧道掌子面前方地质所有可能性电阻率分布范围，准确预测突水层的位置，为隧道安全施工提供保障，具有较大的实际应用价值。

附图说明

图1为隧道瞬变电磁探测方式示意图；

图2为变维贝叶斯算法流程图；

图3为隧道瞬变电磁含噪声接收信号仿真图；

图4为基于变维贝叶斯算法隧道瞬变电磁数据反演结果；

图5为变维贝叶斯反演地层维数概率分布图；

图6为变维贝叶斯反演结果层电阻率与层厚度对应关系图。

具体实施方式

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

如图1所示，将瞬变电磁仪器放置于隧道掌子面处进行探测，获取记录隧道掌子面处的瞬变电磁观测数据，其中瞬变电磁仪器是适用于隧道探测的小型瞬变电磁系统，采用重叠回线配置，发射系统贴靠掌子面发射双极性梯形波，接收系统同步采集观测信号数据；

⑵、对采集到的瞬变电磁数据进行数据预处理，包含下列步骤：

a、将相同测点位置的瞬变电磁观测数据进行叠加，叠加次数与该测点位置上瞬变电磁仪的发射次数一致；

b、瞬变电磁发射波形为双极性梯形波，对于叠加后的数据将其正负极性半周期发射关断后对应的数据相减；

c、对完成上述操作的测点数据做对数等间隔数据抽道；

⑶、本实施例将隧道模型设置为，隧道掌子面前方80米电阻率模型为掌子面前方20米到50米间存在低阻含水构造，对应电阻率为20Ω.m，其余地方为高阻，对应电阻率为300Ω.m，使用发射边长为6米，匝数为5匝，发射电流为10A的上述瞬变电磁系统对隧道进行探测，加入均值为10nV的随机噪声，模型正演观测信号d_obs如图3所示；

⑷、采用变维贝叶斯算法对预处理后的数据进行反演解释，其算法流程图如图2所示，包括以下具体步骤：

①设置反演初始模型及约束参数

在反演初始需要设置模型参数的初始值与范围作为先验信息，包括：初始模型电阻率ρ＝100Ω.m、初始地层数k＝1、最小最大地层层数(k_min＝1,k_max＝4)、最小最大层界面位置(z_min＝0m,z_max＝80m)、最小和最大电阻率变化范围(ρ_min＝10,ρ_max＝400)，定义最小层厚度h_min为：

h_min＝(z_max-z_min)/(2·k_max)＝10m；

②根据当前模型参数定义后验概率分布:

其中，p(m)为先验概率；p(d_obs|m)为似然函数；d_obs为观测数据；m为模型参数；k为模型层数；z为层界面深度；ρ为模型电阻率；p(k)、p(z|k)和p(ρ|k,z)分别为层数、界面位置、电阻率的先验概率分布，p(d_obs|k,z,ρ)为表征数据拟合程度的似然函数，它们计算公式如下：

其中，k_min与k_max分别代表设置的地层层数的最小值与最大值，这里设置为k_min为1，k_max为4；

其中，Δz代表新界面的可用空间；

其中，N_d为数据个数，本实施例为30个；F(m)为正演模型响应；C_d为数据方差矩阵；

其中，ρ₀为全空间最佳匹配电阻率，可以通过最小二乘法进行数据拟合得到；

其中ρ_r表示设置的电阻率变化范围，本实施例设置最小电阻率为10Ω.m，最大为400Ω.m。

③计算建议分布产生候选模型参数

在变维贝叶斯反演中，候选模型是由当前模型的建议分布函数生成的，与之前的反演模型无关，让贝叶斯反演受初始模型的影响较小，候选模型的建议分布由下式定义：

q(m'|m)＝q(k'|k)q(z'|k',z)q(ρ'|k',z',ρ)

理想的建议分布应与后验分布相同，以便加快模型的收敛速度，采用可跳跃马尔科夫链蒙特卡洛方法进行采样，该采样过程包括以下四种基本状态：

e、新层生成，在满足不超过最大设置层数且大于最小层厚度的前提下，在最大和最小层界面位置间随机生成一个新的界面，层数加1；

f、旧层灭亡，在现有的层界面中，随机删除一个界面，层数减1；

g、扰动更新，在层界面数保持不变的条件下，随机选择以层界面，将其所在的位置在(-h_min,h_min)范围内扰动；

h、保持不变，层数与界面位置都不变，只对电阻率进行反演。

上述四种基本状态的概率满足：

相应层界面的位置按照均匀分布给出满足：

电阻率的建议分布定义为以当前模型电阻率为均值的多维正态分布，满足：

其中，ρ_k'与

分别为候选模型的电阻率及对应候选模型界面位置的当前电阻率，C_f为后验电阻率方差，可由当前电阻率的正演模拟值与观测数据的二范数来定义：C_f＝λ·||F(m)-d_obs||₂，其中，λ为缩放系数；本实施例选取缩放系数λ为10000。

④根据后验概率分布与候选模型建议分步计算接受概率

候选模型产生后，按照接受概率判断是否接受候选模型，定义式子为：

其中，|J|为候选模型对当前模型导数构成的雅克比矩阵，取1，接受概率化简为：

在0到1间随机产生一个数，比较该数与α的大小，如果α较大，则候选模型被接受，否则被拒绝接受，继续按建议分布生成新的候选模型，直至模型收敛。所有被接受的模型共同构成反演结果。模型收敛通过设置最大模型采样数来判断，当模型采样数达到该设定值时停止搜索。本实施例设置最大模型采样数达到500万次时即为收敛，设置模型接收约束值为0.001；

⑸、隧道模型反演结果如图4所示，地层层数概率分布如图5所示，可以得出反演结果中最大概率对应的电阻率及层界面的位置信息与设置的模型基本一致，所有接收的模型都在设置模型的附近波动，进一步根据变维贝叶斯反演结果，将层界面厚度信息、对应的电阻率信息分别存入矩阵，建立模型层界面厚度与模型电阻率的对应关系，如图6所示；

⑹、根据变维贝叶斯反演结果，可以较为准确的判断出电阻率分层的层界面所在位置的信息，能够得到每一地层所在位置的全部电阻率分布范围，随着地层深度的增加，反演结果的不确定度也增大，体现在深层地层的电阻率分布范围与散度逐渐增大，因此采用变维贝叶斯反演算法，能够对隧道瞬变电磁数据做相对全面的反演解释，给出电阻率与层界面位置的分布情况。

Claims (4)

1.一种基于变维贝叶斯的隧道瞬变电磁数据解释方法，其特征在于，该方法包括以下步骤：

S1、使用瞬变电磁仪器，放置于隧道掌子面处进行探测，获取隧道掌子面含水构造的观测数据；

S2、对获得的观测数据进行瞬变电磁数据预处理；

S3、使用变维贝叶斯算法对预处理后的数据反演，获得掌子面前方所有可能的地质电阻率分布及层界面位置信息，提取其中最大概率对应的电阻率及层界面位置分布信息；

S4、根据最大概率地质电阻率分布及层界面位置信息，给出隧道掌子面前方含水构造异常情况，并进行不确定度分析，得出准确的含水构造分布范围；

所述步骤S1中的瞬变电磁仪器是适用于隧道探测的小型瞬变电磁系统，采用重叠回线配置，发射系统贴靠掌子面发射双极性梯形波，接收系统同步采集观测信号数据；

所述步骤S2中瞬变电磁数据预处理的具体步骤包括：

1)、将相同测点位置的瞬变电磁观测数据进行叠加，叠加次数与该测点位置上瞬变电磁仪的发射次数一致；

2)、瞬变电磁发射波形为双极性梯形波，对于叠加后的数据将正负极性半周期发射关断后对应的数据相减；

3)、对完成步骤2)操作的测点数据做对数等间隔数据抽道；

所述步骤S3中的变维贝叶斯算法的包括如下具体步骤；

①设置反演初始模型及约束参数

在反演初始需要设置模型参数的初始值与范围作为先验信息，包括：初始模型电阻率ρ、初始地层数k、初始层界面位置z、最小地层层数k_min和最大地层层数k_max、最小层界面位置z_min和最大层界面位置z_max、最小和最大电阻率值ρ_min,ρ_max，定义最小层厚度h_min为：

h_min＝(z_max-z_min)/(2·k_max)；

②根据当前模型参数定义后验概率分布:

其中，p(m)为先验概率；p(d_obs|m)为似然函数；d_obs为观测数据；m为模型参数；p(k)、p(z|k)和p(ρ|k,z)分别为层数、界面位置、电阻率的先验概率分布，p(d_obs|k,z,ρ)为表征数据拟合程度的似然函数，计算公式如下：

其中，Δz代表新界面的可用空间；

其中，N_d为数据个数；F(m)为正演模型响应；C_d为数据方差矩阵；

其中，ρ₀为全空间最佳匹配电阻率，通过最小二乘法进行数据拟合得到；

其中ρ_r表示设置的电阻率变化范围；

③计算建议分布产生候选模型参数

在变维贝叶斯反演中，候选模型是由当前模型的建议分布函数生成，与之前的反演模型无关，让贝叶斯反演受初始模型的影响较小，候选模型的建议分布由下式定义：

q(m'|m)＝q(k'|k)q(z'|k',z)q(ρ'|k',z',ρ)

理想的建议分布与后验分布相同，以便加快模型的收敛速度，采用可跳跃马尔科夫链蒙特卡洛方法进行采样，该采样过程包括以下四种基本状态：

a、新层生成，在满足不超过最大设置层数且大于最小层厚度的前提下，在最大和最小层界面位置间随机生成一个新的界面，层数加1；

b、旧层灭亡，在现有的层界面中，随机删除一个界面，层数减1；

c、扰动更新，在层界面数保持不变的条件下，随机选择以层界面，将其所在的位置在(-h_min,h_min)范围内扰动；

d、保持不变，层数与界面位置都不变，只对电阻率进行反演；

上述四种基本状态的概率满足：

相应层界面的位置按照均匀分布给出满足：

电阻率的建议分布定义为以当前模型电阻率为均值的多维正态分布，满足：

其中，ρ_k'与

分别为候选模型的电阻率及对应候选模型界面位置的当前电阻率，C_f为后验电阻率方差，可由当前电阻率的正演模拟值与观测数据的二范数来定义：C_f＝λ·||F(m)-d_obs||₂，其中，λ为缩放系数；

④根据后验概率分布与候选模型建议分步计算接受概率

候选模型产生后，按照接受概率判断是否接受候选模型，定义式子为：

其中，|J|为候选模型对当前模型导数构成的雅克比矩阵，取1，接受概率化简为：

在0到1间随机产生一个数，比较该数与α的大小，如果α较大，则候选模型被接受，否则被拒绝接受，继续按建议分布生成新的候选模型，直至模型收敛，所有被接受的模型共同构成反演结果。

2.按照权利要求1所述的方法，其特征在于，所述的步骤S1中的隧道掌子面含水构造的探测范围为掌子面前方0至80米以内。

3.按照权利要求1所述的方法，其特征在于，所述步骤④中的模型收敛通过设置最大模型采样数来判断，当模型采样数达到该设定值时停止搜索。

4.按照权利要求1所述的方法，其特征在于，所述步骤S4中不确定度的分析根据变维贝叶斯反演结果，反演结果由所有被接受的模型构成，包含层界面位置信息和电阻率信息，根据层界面位置信息得到层厚度，将层界面厚度信息、对应的电阻率信息建立模型层界面厚度与模型电阻率的关系图，根据关系图得出层界面厚度与对应电阻率的分布范围。

Images