CN116755146A - 基于时频域振幅相位联合的弹性波全波形反演方法与煤田陷落柱建模 - Google Patents

Abstract

本发明公开了一种基于时频域振幅相位联合的弹性波全波形反演方法与煤田陷落柱建模，涉及一种煤田陷落柱弹性参数建模方法。首先，对地震信号进行时频变换，构建时频域振幅相位联合弹性波全波形反演目标函数；其次，引入振幅调节因子，增强弹性波全波形反演目标函数中相位信息，同时减弱地震数据振幅信息；然后，在目标函数中引入全变分正则化约束，提高弹性参数反演精度；最后，推导基于全变分约束的时频域振幅相位联合弹性波全波形反演的梯度算子，并利用局部优化算法对弹性参数进行更新迭代。提高了煤田陷落柱的成像精度和识别准确率。

Description

基于时频域振幅相位联合的弹性波全波形反演方法与煤田陷 落柱建模

技术领域

本发明涉及一种煤田陷落柱弹性参数建模方法，具体说是一种全变分约束的时频域振幅相位联合弹性波全波形反演方法，以及该方法在煤田陷落柱弹性参数建模中的应用。

背景技术

煤田地震勘探经历60多年的发展，从二维地震勘探发展到三维地震勘探，从原先折射波勘探到转变为多波多分量勘探，技术的进步使得煤田地震探测精度有了大幅度提高(彭苏萍等，2002；彭苏萍，2020；王超越和董守华，2020)。煤田隐伏陷落柱诱发矿井突水和瓦斯突出等灾害性事故，一直是煤矿安全开采的重大威胁，同时陷落柱探测也一直是煤田地震勘探的一项重要任务(徐佩芬等，2009；Wu等，2020)。由于陷落柱多呈陡倾角锥形，且柱面反射信号弱，导致常规反射波地震成像方法无法对陷落柱高陡侧翼进行有效地成像。为此，王贇等(2000)分析了煤层陷落柱的地震波场特征，指出陷落柱难以形成连续的反射波同相轴，并指出陷落柱柱璧绕射波是陷落柱地震响应的重要表示形式。曹志勇等(2012)利用共散射点道集，通过判断散射点的空间位置，进而间接地勾画陷落柱高陡侧翼形状，并在模拟数据和实际数据中得到了验证。朱广辉等(2020)利用纵波和转换波信息来对地下煤层陷落柱进行成像，数值测试结果表明利用PS波成像结果相对于PP波具有更高的分辨率，并能够更好地识别煤层陷落柱的真实位置。刘健等(2022)指出可以将地震记录中的绕射波提取出来，单纯利用绕射波来对陷落柱进行成像，证实了利用绕射波对陷落柱的成像明显优于反射波成像结果，可在一定程度上提高陷落柱的识别准确率。

此外，矿井地震也可以有效地对煤层陷落柱进行探测，且具有离目标体近、探测精度高等优势(滕吉文等，2013)。朱培民和蒋锦朋(2019)将弹性波全波形反演方法应用到煤层槽波信号反演中，并与常规的煤层槽波成像方法进行对比，发现槽波全波形反演能够更好地识别煤层中的小异常体，同时具有更高的成像分辨率。方金伟等(2022)利用弹性波全波形反演方法在回采工作面中对复杂陷落柱进行反演，指出横波相对于纵波信息能够更好地反演陷落柱的速度结构。但是，井下槽波只能在煤层内部传播，致使槽波地震勘探方法难以可靠识别煤层底板下部的隐伏陷落柱(程建远等，2018)。

时间域全波形反演方法首先由Lailly(1983)和Tarantola(1984)提出，其能够利用全波场的振幅和相位信息获取地下复杂构造的速度场信息。在此基础上，Pratt等(1998)将全波形反演方法推广到频率域，并指出只需几个离散频率的地震数据便可以达到时间域全波形反演精度。但是，要想获得煤田深部复杂构造的高精度速度反演结果，必须缓解全波形反演方法对初始模型的依赖，并提高煤田深部弱地震反射信号的可反演精度。地震数据的相位信息相比于振幅信息与地下介质有着更好的线性对应关系，利用相位信息构建速度反演目标函数，有助于缓解全波形反演方法对初始模型的依赖。Fichtner等(2008)在时频域构建相位目标函数，利用地震台站数据实现对全球尺度的全波形反演。等(2011)利用复地震道数据，提出瞬时相位目标函数，指出相位信息相比振幅信息可以获得更好的反演结果。Alkhalifah等(2012)针对地震数据的相位缠绕问题，提出利用解缠以后的相位信息能够有效缓解全波形反演对初始模型依赖。因此，可以利用振幅相位联合目标函数，来解决煤田陷落柱弹性参数的高精度速度建模问题。

发明内容

本发明针对以上问题，提出了一种基于时频域振幅相位联合的弹性波全波形反演方法与煤田陷落柱建模，利用弹性波数据的时频域局部化特征，采用了基于时频域目标函数的煤田陷落柱弹性参数建模方法，提高了煤田陷落柱的成像精度和识别准确率。

本发明的技术方案为：包括以下步骤：

步骤1、利用地震数据处理软件Geoeast对地震数据进行预处理，并将预处理后的地震数据作为时频域振幅相位联合的弹性全波形反演的数据输入。

步骤2、构建煤田陷落柱弹性参数初始模型，作为弹性波方程正演模拟的模型输入。

步骤3、定义观测系统，设定震源子波；利用初始弹性参数模型进行正演模拟，存储正传波场，其中正演模拟对应的弹性波方程为：

其中x,z为坐标轴，λ和μ为拉梅常数，ρ为密度，u＝[u_x,u_z]表示位移分量，σ＝[σ_xx,σ_zz,σ_xz]表示应力分量，t为时间；综合计算效率和模拟精度，优选用有限差分法进行正演模拟。

步骤4、利用Gabor变换获得观测数据和模拟数据的时频域振幅相位信息，其中观测数据和模拟数据的Gabor变换为：

其中h(τ-t)表示高斯窗函数；τ为中间变量；ω为角频率，u(t)和d(t)分别表示时间域模拟数据和观测数据；和/>表示时频域模拟数据和观测数据，F_h[·]表示对地震数据作用的Gabor变换算子。

步骤5、构建基于全变分约束的时频域振幅相位联合的弹性波全波形反演目标函数：

其中i表示虚数；ns为震源数目；nr为检波器数目；和/>表示模拟数据和观测数据在时频域的相位信息；/>和/>表示观测数据和模拟数据在时频域的振幅信息；/>表示时频域地震数据的位移分量，m＝[v_p,v_s]表示弹性波的P波速度和S波速度；||·||_TV表示全变分约束；λ_TV是用来调节数据拟合与模型约束的权重，其中当λ_TV＝0，上式不考虑模型约束；η＝max(|g_data|)/max(|g_TV|)是归一化权重因子，使数据拟合部分的更新量与模型约束保持同一个数量级；ε∈[0,1]表示振幅权重因子，用来控制目标函数中波形振幅和相位权重的比例。

步骤6、构建弹性波全波形反演梯度算子：

其中表示正向传播波场，/>表示逆时外推波场。其中拉梅常数[λ，μ]和地震波速度[v_p,v_s]有如下关系：

最终，目标函数对弹性波速度参数[v_p,v_s]的梯度可以表示为：

步骤7、计算目标函数关于模型参数的偏导数：

其中*为复共轭；Re[·]取数据的实部；/> 是数据的残差。其中时频域模拟数据的绝对值对模型参数的偏导数可以表示为：

其中全变分约束对模型参数的偏导数，可以表示为：

其中对应的差分格式可以表示为：

因此，目标函数关于模型参数的偏导数可以表示为：

步骤8、定义弹性波全波形反演伴随震源，并将伴随震源反传至模型空间，获得反传波场，用于计算梯度；对应的伴随震源为：

其中对应着水平分量和垂直分量的伴随震源。

步骤9、利用L-BFGS优化算法计算模型扰动的更新方向，迭代公式如下：

m_k+1＝m_k-α_kH_kg_k

其中m_k为第k步模型扰动的更新量，α_k为步长，H_k为近似海森矩阵的逆，为模型扰动的更新梯度。在L-BFGS优化更新中，只需要保存少数的向量对，用于更新Hessian矩阵，其更新公式如下：

H_k+1＝V_k ^TH_kV_k+ρ_ks_ks_k ^T

V_k＝I-ρ_ky_ks_k ^T

s_k＝m_k+1-m_k，y_k＝g_k+1-g_k

其中H_k+1是根据向量对{s_k,y_k}和H_k计算得到；H_kg_k的乘积可以通过梯度g_k与向量对{s_k,y_k}之间一系列向量的内积与向量的和来获得。其中的近似Hessian矩阵的逆矩阵H_k需满足以下更新公式：

步骤10、判断是否满足反演精度和迭代次数的终止条件，若满足则输出煤田陷落柱时频域振幅相位联合的弹性波全波形反演结果，获得煤田陷落柱的弹性参数模型；若不满足终止条件，将当前的反演结果继续作为下一次循环的初始弹性参数模型，直到满足终止条件。

本发明充分利用煤田地震信号的时频域振幅相位信息，提高煤田陷落柱的弹性参数建模精度。首先，对地震信号进行时频变换，构建时频域振幅相位联合的弹性波全波形反演目标函数；其次，引入振幅调节因子，增强时频域目标函数中相位信息的权重，减弱地震数据振幅信息对深部反演结果影响，提高煤田陷落柱深部弹性参数建模精度；最后，推导时频域振幅相位联合的弹性波全波形反演目标函数对应的梯度算子，并利用优化算法对弹性参数模型进行更新迭代。通过煤田模型数据测试，验证了本发明方法可以获得煤田陷落柱的高精度弹性参数反演结果。

附图说明

图1是本发明基于全变分约束的时频域振幅相位联合弹性波全波形反演方法，及其在煤田陷落柱速度建模上的应用测试对应的流程框图；

图2是煤田陷落柱模型，(a)真实Vp模型；(b)初始Vp模型；(b)真实Vs模型；(d)初始Vs模型；

图3是煤田陷落柱的常规多尺度弹性波全波形反演结果，(a)低频段Vp反演结果(6-8Hz)；(b)高频段Vp反演结果(6-12Hz)；(c)低频段Vs反演结果(6-8Hz)；(d)高频段Vs反演结果(6-12Hz)；

图4是煤田陷落柱的全变分约束多尺度弹性波全波形反演结果，(a)低频段Vp反演结果(6-8Hz)；(b)高频段Vp反演结果(6-12Hz)；(c)低频段Vs反演结果(6-8Hz)；(d)高频段Vs反演结果(6-12Hz)；

图5煤田陷落柱的全变分约束时频域振幅相位联合多尺度弹性波全波形反演结果，(a)低频段Vp反演结果(6-8Hz)；(a)低频段Vs反演结果(6-8Hz)；(c)高频段Vp反演结果(6-12Hz)；(d)高频段Vs反演结果(6-12Hz)。

具体实施方式

为能清楚说明本发明的技术特点，下面通过具体实施方式，并结合其附图，对本发明进行详细阐述。

由于煤层与围岩之间存在较强的波阻抗分界面，对地震波的传播具有很强的屏蔽作用，因此很难获得煤层下部构造的高精度成像结果。此外，由于煤田陷落柱陷是由灰岩经过溶蚀塌陷形成的，且多为高陡构造，内部充填介质差异性大，这给常规弹性波全波形反演带来了巨大挑战。为了解决这个问题，本发明使用基于全变分约束的时频域振幅相位联合弹性波全波形反演方法，充分利用弱地震信号的时频域相位信息，提高煤田复杂构造的成像精度，具体流程如图1所示。真实的煤田陷落柱Vp和Vs模型如图2a和2b所示，初始Vp和Vs模型如图2c和2d所示。在地表均匀分布40个震源，每个震源对应着300个检波器。地震数据记录时间为4.4s，时间间隔为2ms。震源采用主频为12Hz的雷克子波，同时利用L-BFGS优化算法来计算下降方向，并对模型扰动进行更新迭代。

本发明使用模型数据进行测试，首先利用真实的煤田陷落柱Vp和Vs模型(图2a、2b)，结合弹性波波动方程正演获得观测数据，并将其视为实际生产过程使用的地震数据。根据本案的流程：首先对地震数据进行预处理，构建初始Vp和Vs模型(2c、2d)，定义观测系统，以及提取地震数据对应的震源子波信息。在Gabor时频变换的框架下获得时频域模拟数据和观测数据，建立时频域振幅相位联合弹性波全波形反演目标函数，并求取目标函数相对于模型扰动的偏导数，进而获得模型扰动的更新量。此外，本发明还结合了全变分正则化方法来提高煤层和陷落柱的反演精度。最后利用L-BFGS优化算法在初始模型的基础上不断更新模型扰动。

本发明对比了常规多尺度弹性波全波形反演方法(图3)，基于全变分约束的多尺度弹性波全波形反演方法(图4)，以及基于全变分约束的时频域振幅相位联合多尺度弹性波全波形反演方法(图5)。

模型参数如下：

表1：基于全变分约束的时频域振幅相位联合多尺度弹性波全波形反演方法测试参数

将图3与图4对比可以看出，基于全变分约束的弹性波全波形反演相对于常规弹性波全波形反演的速度结构更均匀，尤其是在陷落柱和大尺度均匀地层内部的速度反演结果明显优于常规常规弹性波全波形反演方法。因此，基于全变分约束的弹性波全波形反演结果更接近真实模型扰动。将图3与图4和图5相对比可以看出，基于全变分约束的时频域振幅相位联合多尺度弹性波全波形反演结果在深部区域速度建模方面有了明显地优势。由于时频域振幅相位信息具有与地下速度模型更好的线性对应关系，极大程度上减弱了弹性波全波形反演方法对初始速度模型的依赖。此外，在时频域目标函数中引入权重因子，不仅增强了反演过程的稳定性，而且缓解扰动波场匹配过程中出现的周波跳跃现象。数值试验结果表明，基于全变分约束的时频域振幅相位联合多尺度弹性波全波形反演方法在煤田陷落柱复杂构造的高精度弹性参数建模方面具有一定的优势。

本发明利用地震信号的时频域局部化特征和弱地震信号的时频域振幅和相位信息，提高煤田陷落柱深部构造的反演精度。首先，对地震信号进行时频变换，构建时频域振幅相位联合弹性波全波形反演目标函数；其次，引入振幅调节因子，增强弹性波全波形反演目标函数中相位信息，同时减弱地震数据振幅信息；然后，在目标函数中引入全变分正则化约束，提高弹性参数反演精度；最后，推导基于全变分约束的时频域振幅相位联合弹性波全波形反演的梯度算子，并利用局部优化算法对弹性参数进行更新迭代。通过煤田陷落柱模型的数值测试，验证了本发明方法可以获得高精度的煤田陷落柱弹性参数建模结果。

本发明具体实施途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进，这些改进也应视为本发明的保护范围。

Claims (8)

1.一种基于时频域振幅相位联合的弹性波全波形反演的煤田陷落柱建模方法，其特征在于，包括以下步骤：

步骤1、利用地震数据处理软件Geoeast对地震数据进行预处理，并将预处理后的地震数据作为时频域振幅相位联合的弹性全波形反演方法的数据输入；

步骤2、构建弹性参数初始模型，作为弹性波方程正演模拟的模型输入；

步骤3、定义观测系统，设定震源子波；利用初始弹性参数模型进行正演模拟，存储正传波场；

步骤4、利用Gabor变换获得观测数据和模拟数据的时频域振幅相位信息；

步骤5、构建基于全变分约束的时频域振幅相位联合的弹性波全波形反演目标函数；

步骤6、构建弹性波全波形反演梯度算子；

步骤7、计算目标函数关于模型参数的偏导数；

步骤8、定义弹性波全波形反演伴随震源，并将伴随震源反传至模型空间，获得反传波场，用于计算梯度；

步骤9、利用L-BFGS优化算法计算模型扰动的更新方向；

步骤10、判断是否满足终止条件，若满足则输出煤田陷落柱时频域振幅相位联合的弹性波全波形反演结果，获得煤田陷落柱的弹性参数模型；若不满足终止条件，将当前的反演结果继续作为下一次循环的初始弹性参数模型，直到满足终止条件。

2.根据权利要求1所述的一种基于时频域振幅相位联合的弹性波全波形反演的煤田陷落柱建模方法，其特征在于，步骤3中正演模拟对应的弹性波方程为：

其中x,z为坐标轴，λ和μ为拉梅常数，ρ为密度，u＝[u_x,u_z]表示位移分量，σ＝[σ_xx,σ_zz,σ_xz]表示应力分量，t为时间；综合计算效率和模拟精度，优选用有限差分法进行正演模拟。

3.根据权利要求1所述的一种基于时频域振幅相位联合的弹性波全波形反演的煤田陷落柱建模方法，其特征在于，步骤4中观测数据和模拟数据的Gabor变换为：

其中h(τ-t)表示高斯窗函数；τ为中间变量；ω为角频率，u(t)和d(t)分别表示时间域模拟数据和观测数据；和/>表示时频域模拟数据和观测数据，F_h[·]表示对地震数据作用的Gabor变换算子。

4.根据权利要求1所述的一种基于时频域振幅相位联合的弹性波全波形反演的煤田陷落柱建模方法，其特征在于，步骤5中构建的函数为：

其中i表示虚数；ns为震源数目；nr为检波器数目；和表示模拟数据和观测数据在时频域的相位信息；/>和表示观测数据和模拟数据在时频域的振幅信息；/>表示时频域地震数据的位移分量，m＝[v_p,v_s]表示弹性波的P波速度和S波速度；||·||_TV表示全变分约束；λ_TV是用来调节数据拟合与模型约束的权重，其中当λ_TV＝0，上式不考虑模型约束；η＝max(|g_data|)/max(|g_TV|)是归一化权重因子，使数据拟合部分的更新量与模型约束保持同一个数量级；ε∈[0,1]表示振幅权重因子，用来控制目标函数中波形振幅和相位权重的比例。

5.根据权利要求1所述的一种基于时频域振幅相位联合的弹性波全波形反演的煤田陷落柱建模方法，其特征在于，步骤6中构建的梯度算子为：

其中表示正向传播波场，/>表示逆时外推波场；其中拉梅常数[λ，μ]和地震波速度[v_p,v_s]有如下关系：

最终，目标函数对弹性波速度参数[v_p,v_s]的梯度可以表示为：

6.根据权利要求1所述的一种基于时频域振幅相位联合的弹性波全波形反演的煤田陷落柱建模方法，其特征在于，步骤7计算目标函数关于模型参数的偏导数，具体为：

其中*为复共轭；Re[·]取数据的实部；/> 是数据的残差；其中时频域模拟数据的绝对值对模型参数的偏导数可以表示为：

其中全变分约束对模型参数的偏导数，可以表示为：

其中对应的差分格式可以表示为：

因此，目标函数关于模型参数的偏导数可以表示为：

7.根据权利要求1所述的一种基于时频域振幅相位联合的弹性波全波形反演的煤田陷落柱建模方法，其特征在于，步骤8对应的伴随震源为：

其中对应着水平分量和垂直分量的伴随震源。

8.根据权利要求1所述的一种基于时频域振幅相位联合的弹性波全波形反演的煤田陷落柱建模方法，其特征在于，步骤9所用的迭代公式如下：

m_k+1＝m_k-α_kH_kg_k

其中m_k为第k步模型扰动的更新量，α_k为步长，H_k为近似海森矩阵的逆，为模型扰动的更新梯度；在L-BFGS优化更新中，只需要保存少数的向量对，用于更新Hessian矩阵，其更新公式如下：

H_k+1＝V_k ^TH_kV_k+ρ_ks_ks_k ^T

V_k＝I-ρ_ky_ks_k ^T

s_k＝m_k+1-m_k，y_k＝g_k+1-g_k

其中H_k+1是根据向量对{s_k,y_k}和H_k计算得到；H_kg_k的乘积可以通过梯度g_k与向量对{s_k,y_k}之间一系列向量的内积与向量的和来获得；其中的近似Hessian矩阵的逆矩阵H_k需满足以下更新公式：