FUTTERMAN^[1]在介质对地震波吸收衰减效应的研究中指出, 岩石的吸收衰减特征是地层的基本属性之一, 因此地层的吸收衰减参数的精确求取问题一直是研究的热点。而与深部地层相比, 近地表地层对地震波的吸收衰减作用更为强烈, 但在处理过程中, 很难准确地求取地层吸收衰减参数, 尤其是近地表Q值参数。目前的研究已经给出了多种Q值估计方法, 且各种方法都有优缺点和相应的应用范围。基于数据来源的不同, 可将已有的方法分为VSP数据和地面反射地震数据两大类, 而微测井可以归为VSP数据类。基于VSP数据类的算法主要包括谱比法、质心频移法及峰值频移法等。BATH^[2]提出了谱比法, 即提取两个时间或深度处的子波频谱, 使用两振幅谱谱比的斜率来计算Q值。GLADWIN等^[3]基于衰减时脉冲宽度增加的现象提出了上升时间原理, KJARTANSSON^[4]依据上升时间原理提出了采用上升时间法求Q值的方法。JANNSEN等^[5]提出子波模拟法, 通过改变Q值进行正演, 选择正演结果与实际信号最接近的模拟Q值作为近似, 之后依照近似原理又出现了频率域近似的频谱模拟法。TONN^[6]使用合成VSP资料对不同Q值估计方法(包括谱比法、解析信号法、振幅衰减法、子波模拟法、拟合技术等)进行了对比分析, 认为大多数方法求取的Q值精度在很大程度上取决于地震资料的质量, 且方法大都建立在特定假设的条件下, 实际应用时缺乏普适性。

近年来, 诸多学者在前人研究的基础上提出了一些改进算法, 其中, 赵秋芳等^[7]梳理了现有近地表品质因子估算方法, 将其初步划分为岩石样本测试Q值估算和地层原位测量Q值估算两大类; 金子奇等^[8]采用对数谱比反演方法, 利用多道反射波信息同时反演多道衰减旅行时, 从而避免了传统方法中人为选择最优频带带来的误差; 陈文爽等^[9]、王小杰等^[10]和王德利等^[11]基于S变换的时频谱作为谱比法的输入消除传统谱比法在计算地层吸收参数时的平均效应; 张繁昌等^[12]采用自适应分解的子波通过谱比法来分析薄层干涉问题; 李伟娜等^[13]借鉴微测井分层速度回归分析思想, 提出了一种双线性回归稳定Q值估计算法提高了精度; 王宗俊^[14]利用谱模拟获得的子波谱约束质心频移法进行Q值估算; 张立彬等^[15]引入积分中值参变量避免了对质心频移法中震源子波的先验假设; 高静怀等^[16]利用特征结构法改善地震子波峰值频率提取精度提高Q值估算精度; 冯玮等^[17]采用时间域子波匹配法通过子波褶积矩阵引入衰减响应, 从根本上避免谱估计的误差; 郭锐等^[18]引入Capon2D Q值估计新方法, 对其求解过程进行了加权改进; 宋吉杰等^[19]提出一种基于信息融合的近地表介质Q值估计方法以及稳定的反Q滤波处理技术, 提高了地震资料的分辨率; 王静等^[20]利用零井源距VSP资料, 采用改进的频谱拟合法进行稳定的深层Q值估算; 苏勤等^[21]基于地表一致性原理, 引入表层相对衰减系数的概念, 在共炮检域迭代计算相对衰减系数并实现稳定的表层Q值空变补偿。上述方法重点是在Q值求取算法中进行了各种改进, 本文重点则是要在微测井地震数据初至波分离、走时自动提取和频谱求取方面进行改进。MALLAT^[22]首次提出了一种使用Gabor小波作为时频原子库的匹配追踪方法。LIU等^[23]进一步提出了使用Ricker子波库和Morlet子波库的快速匹配追踪分解方法, 并在Morlet原子库中将子波尺度参数作为常数。WANG^[24]将尺度参数改为可变参数, 并对Morlet原子库进行了改进, 使残差下降速度更快, 完备库中的原子数更多。SONG^[25]提出了一种基于快速匹配追踪的近地表Q值估计方法, 提高了Q值估计的精度和抗噪能力, 但使用零相位Ricker小波原子在匹配追踪提取初至时会产生一定误差。在无噪情况下, 常规的谱比法是可靠的方法, 当信噪比较高时, 谱比法依赖于诸多因素, 如时窗长度、时窗形状、求取斜率的起止频率等, 当信噪比较低时, 谱比法可靠性差, 因此, 需要深入研究谱比法解的稳定性。本文基于非零相位雷克子波库的复数域快速匹配追踪算法, 实现微测井资料初至波分离、走时自动拾取和频谱求取, 并通过引入整形正则化算子和最小平方反演方法求取谱比值, 提出了一种新的近地表Q值估计方法。

微测井观测点在整个地震工区比较稀疏, 使得利用稀疏采样点信息构建三维近地表模型较为困难, 传统的双线性、双三次样条及克里金插值建模的方法很难获得较好的建模效果。近年来, 机器学习和深度学习被广泛应用于多个领域, 在石油勘探、开发领域也取得较好的应用效果。RÖTH等^[27]通过神经网络反演了地下速度, 证明了机器学习方法能够解决反演问题。MOHAMED等^[28]使用神经网络完成了地震数据的叠后反演。KORJANI等^[29]通过深度学习算法预测了地下岩石的物理参数。DAHLKE等^[30]提出使用深度学习方法直接根据原始地震记录进行断层预测, 在模型数据中实验准确性较高。基于深度学习的多元非线性回归方法, 利用工区范围内和周边工区大量的微测井数据作为样本进行深度学习, 在一个近地表条件相对稳定的地区, 可以利用大量的已知样本数据建立统一的多元非线性回归算子作为初始回归模型, 再通过迁移学习方法, 将目标工区内新获得的微测井数据作为新的学习标签对多元非线性回归进行调优, 不断优化多元非线性回归算子, 以保证近地表模型的建模精度。本文将工区内不同位置的微测井数据估计得到的Q值函数作为已知样本标签, 并通过训练深度神经网络来构建多元非线性回归算子, 实现了三维近地表品质因子建模, 并将建立的三维近地表模型用于地震资料处理, 以验证方法的有效性。

1 方法原理 1.1 复数域快速匹配追踪分解

常规的匹配追踪算法要求在每次匹配中扫描全部的原子库, 但原子库是过完备的, 扫描过程计算量较大, 计算速度较慢。为了减少匹配中的扫描工作量, 可以使用快速复数域匹配追踪(Complex domain fast matching pursuit decomposition, CFMP)方法来提高计算效率。CFMP方法需要根据实地震记录获取复地震记录, 从复信号中获取先验信息(如振幅包络最大值处的时间、瞬时频率以及瞬时相位), 然后在先验信息的约束下扫描一小部分的原子库得到最佳匹配, 完成复信号的重构后再将结果返回实数域, 从而在一定程度上减小匹配中的计算量, 提高算法效率。匹配追踪算法的基本思想是将信号投影到一系列时频原子上, 选取的时频原子能很好地匹配信号的局部特征, 利用这些时频原子精确地表示原始信号。

设H是Hilbert空间, 满足$\left\| f \right\| = \int_{ - \infty }^{ + \infty } {{{\left| {f(t)} \right|}^2}} {\rm{d}}t < + \infty {\rm{ 。}}G = \left\{ {{g_r}\mid \;\left\| {{g_r}} \right\| = 1,r \in \mathit{\Gamma }} \right\}$是H中的基元函数, 每一个g_r都经过标准化, Γ是索引。基元函数G是一个时频原子库, 其中每一个g_r都是一个时频原子。通过对复数域信号进行匹配迭代计算, 每一次匹配都得到复时频原子。假设一个复信号F∈H, CFMP算法的目的是将F表示成从G选择的复时频原子的线性组合。信号的复信号可以通过希尔伯特变换获得:

$ F(t) = f(t) + {\rm{i}} \cdot H(f(t)) $

(1)

其中, f(t)为实信号, 虚部H(f(t))为实信号的Hilbert变换, 其实质是实部的90°相移。

同理, 假设g(t)为时频原子, H(g(t))表示时频原子的Hilbert变换, 则复时频原子的表达式为:

$ G(t)=g(t)+\mathrm{i} \cdot H(g(t)) $

(2)

匹配追踪从初始残量R⁰F=F开始, 首先求取残量的瞬时振幅包络、瞬时频率和瞬时相位, 将振幅包络最大值处的时间、瞬时频率以及瞬时相位分别表示为t_j, f_j和φ_j。通过n次匹配计算, 得到:

$ {R^n}F = \left\langle {{R^n}F,{G_{{\gamma _n}}}} \right\rangle {G_{{\gamma _n}}} + {R^{n + 1}}F $

(3)

RⁿF为复残量。其中, G_γn为复时频原子, 且‖G_γn‖=1, 〈RⁿF, G_γn〉表示信号与残差的内积。很容易得出G_γn正交于Rⁿ⁺¹F, 因此有:

$ {\left\| {{R^n}F} \right\|^2} = {\left| {\left\langle {{R^n}F,{G_{{\gamma _n}}}} \right\rangle } \right|^2} + {\left\| {{R^{n + 1}}F} \right\|^2} $

(4)

为了使‖Rⁿ⁺¹F‖最小, 则所选择的复子波G_γn∈G使得|〈RⁿF, G_γn〉|最大, 即最佳原子的选择条件为:

$ {G_{{\gamma _n}}} = \arg \mathop {\sup }\limits_{\gamma \in \varGamma } \left| {\left\langle {{R^n}F,{G_{{\gamma _n}}}} \right\rangle } \right| $

(5)

通过(5)式选择的时频原子能使每次迭代计算所得的残量最小。当信号的局部与匹配原子性质完全相同或者完全匹配时, 复信号在原子上的投影值〈RⁿF, G_γn〉为一实数, 即有real(〈RⁿF, G_γn〉)=〈RⁿF, G_γn〉, 故匹配条件改进为:

$ {G_{{\gamma _n}}} = \arg \mathop {\sup }\limits_{\gamma \in \varGamma } \left| {{\mathop{\rm real}\nolimits} \left( {\left\langle {{R^n}F,{G_{{\gamma _n}}}} \right\rangle } \right)} \right| $

(6)

通过(6)式可计算出最佳匹配子波原子, 即计算出最佳的时移、主频、相位等参数。为了使匹配效果最佳, 本文在构建时频原子库时, 对前面所求参数(时移、主频、相位)的一定邻域内构建时频原子字典库W, W={W_γ|‖W_γ‖=1}, 其中:指标集γ是集合Γ中的元素, 复时频原子都经过单位化。经过匹配条件(6)式, 得到首次匹配的结果为:

$ {R^0}F = {\mathop{\rm real}\nolimits} \left( {\left\langle {{R^0}F,{W_{{\gamma _1}}}} \right\rangle } \right){W_{{\gamma _1}}} + {R^1}F $

(7)

其中, R¹F是首次迭代的复残量, 可直接用于下一次迭代计算, 省去了Hilbert变换所带来的计算量, 计算过程与前一次相同。通过设置迭代次数或误差阈值, 复匹配追踪最后的结果表示为:

$ F(t) = \sum\limits_{n = 0}^{m - 1} {{\mathop{\rm real}\nolimits} } \left( {\left\langle {{R^n}F,{W_{{\gamma _n}}}} \right\rangle } \right) + {R^m}f $

(8)

由于信号中存在各种噪声以及其它干扰因素, 使得使用瞬时信息求得的匹配参数存在误差, 本文使用复数域的最小二乘法对求得的匹配参数进行校正, 使匹配参数所确定的复雷克子波能更接近地震信号的局部信息。设RⁿF=A_jW_γ+Rⁿ⁺¹F, 其中W_γ为利用上面瞬时参数确定所得的复时频原子, A_j=|A_j|·exp(iθ_j)表示对子波原子的校正项, 其中包括振幅校正|A_j|与相位校正θ_j。为使校正后的原子子波能最佳匹配地震信号, 故‖Rⁿ⁺¹F‖取最小值, 使用复数域最小二乘法, 求得:

$ {A_j} = {\left( {\overline {W_\gamma ^{\rm{H}}} \cdot {W_\gamma } + \varepsilon \mathit{\boldsymbol{I}}} \right)^{ - 1}} \cdot \overline {W_\gamma ^{\rm{H}}} \cdot {R^n}F $

(9)

式中: $\overline {W_\gamma ^{\rm{H}}} $表示W_γ的共轭转置; ε为阻尼因子; I为单位矩阵。

通过(9)式求得的A_j得到振幅和相位校正项后的匹配子波原子变为G_γ=|A_j|W_γ(t), 其中匹配参数为γ={u_j, f_j, φ_j}, φ_j=α_j+θ_j为校正之后的相位。由(9)式可计算出最佳匹配子波原子, 即计算出最佳的时移、主频、相位等参数。对复信号求取其先验信息, 包括:振幅包络、瞬时频率和瞬时相位。复信号的振幅包络、瞬时频率和瞬时相位的计算公式分别为:

$ {S(t) = \sqrt {{f^2}(t) + {f^{*2}}(t)} } $

(10)

$ {{f_j} = {{\left. {\frac{1}{{2{\rm{ \mathsf{ π} }}}}\frac{{\rm{d}}}{{{\rm{d}}t}}\arg (F(t))} \right|}_{t = {u_j}}}} $

(11)

$ {{\theta _j} = {{\left. {\arg \frac{{{f^*}(t)}}{{f(t)}}} \right|}_{t = {u_j}}}} $

(12)

其中, f^*(t)为信号的虚部。n次迭代后的振幅系数为:

$ \left| {{A_j}} \right| = \frac{{\mid \langle {R^n}F,W({t_j},{f_j},{\theta _j}\rangle |}}{{{{\left| {W\left( {{t_j},{f_j},{\theta _j}} \right)} \right|}^2}}} $

(13)

另外, 匹配追踪采用非零相位雷克子波构成时频原子库, 具体构建非零相位雷克子波库分两步。首先给出零相位雷克子波W的解析表达式如下:

$ W = \left[ {1 - 2{{({\rm{ \mathsf{ π} }}ft)}^2}} \right]{{\rm{e}}^{ - ({\rm{ \mathsf{ π} }}ft)2}} $

(14)

然后, 对(14)式零相位雷克子波W进行相位旋转, 得到相位旋转φ的雷克子波W_φ表达式为:

$ W_{\varphi}=W \cos \varphi+\operatorname{Imag}(H(W)) \sin \varphi $

(15)

其中, 符号Imag表示取虚部。进一步地, 以振幅包络最大值处的瞬时频率和瞬时相位为基准, 按一定的邻域范围将瞬时频率和瞬时相位左右扩展, 再基于(14)式和(15)式, 形成局部非零相位雷克子波时频原子库。

1.2 整形正则化谱比反演

常规的对数谱比法就是两式作比值并求对数, 即

$ \ln \left| {\frac{{{A_1}(f)}}{{{A_2}(f)}}} \right| = {\rm{ \mathsf{ π} }}\frac{{\Delta t}}{Q}f + c $

(16)

式中:Δt为地震波在某地层中传播的双程旅行时; A₁(f)为衰减后的子波频谱值; A₂(f)为衰减前的子波频谱值; Q为地层的品质因子; c为常量。显然, 谱比对数值为频率f的线性函数, 其斜率为k=π(Δt/Q), 通过线性回归求出斜率k值即可确定Q值。其中关键是求取稳定的对数谱比值。

为了提高谱比计算的稳定性和抗噪能力, 可采用优化稳定反演算法来求取谱比参数。匹配追踪方法可用于信号的分解, 能够提取不同时刻的子波, 并获得相应的幅度谱。子波的谱比表达式可写为:

$ O(f) = \frac{{{A_{i + 1}}W\left( {{t_{i + 1}},f} \right)}}{{{A_i}W\left( {{t_i},f} \right)}} $

(17)

式中:A_iW(t_i, f)是在t_i处匹配得到的最佳原子振幅谱。当分母数值很小时, 除法运算的结果会出现不稳定, 简单的方法是在分母上加一个白噪系数(ξ)来防止不稳定, 即:

$ O(f)=\frac{A_{i+1} W\left(t_{i+1}, f\right)}{A_{i} W\left(t_{i}, f\right)+\xi} $

(18)

其中, ξ的选取和信号的信噪比有关。FOMEL^[26]利用整形正则化来稳定反演问题, 通过引入整形正则化的方法来获得较为稳定的谱比结果。通过引入正则化算子, 将(18)式转化成最小平方问题:

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{O(f)} \sum\limits_f {{{\left\| {{A_i}\left( {{t_i},f} \right) - O(f){A_{i + 1}}\left( {{t_{i + 1}},f} \right)} \right\|}^2}} + }\\ {\xi D(O(f))} \end{array} $

(19)

式中:D是Tikhonov正则化算子; t_i, t_i+1分别为微测井不同深度位置的初至波旅行时。采用共轭梯度迭代来完成正则化反演, 谱比的线性模型可以写成:

$ \boldsymbol{m}=\boldsymbol{L}^{-1} \boldsymbol{d} $

(20)

式中:m是O(f)的矩阵形式; d是A_i+1W(t_i+1, f)的矩阵形式; L是A_iW(t_i, f)的矩阵形式。模型的理论解写为:

$ \mathit{\boldsymbol{\hat m}} = {\left( {{\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{L}} + {\mathit{\boldsymbol{\xi }}^2}{\mathit{\boldsymbol{D}}^{\rm{T}}}\mathit{\boldsymbol{D}}} \right)^{ - 1}}{\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{d}} $

(21)

引入整形正则化算子S=(I+ξ²D^TD)^－1, 则正则化的理论解可以写为:

$ \mathit{\boldsymbol{\hat m}} = {\left[ {{\lambda ^2}\mathit{\boldsymbol{I}} + S\left( {{\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{L}} + {\lambda ^2}\mathit{\boldsymbol{I}}} \right)} \right]^{ - 1}}\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{d}} $

(22)

其中, λ为阻尼调整因子。

1.3 深度学习Q值的建模方法

神经网络是学习算法的一种, 由多层互连节点组成“网络”。受动物神经系统的启发, 神经网络的网络节点类似于神经元, 边缘类似于突触, 如图 1所示。在每个边缘都有一个相关权重, 网络定义了输入数据从网络的输入层传递到输出层的计算规则。用神经网络的网络函数表征输入和输出层之间的关系, 由权重进行参数化。适当定义网络函数后, 可以通过最小化网络权重的代价函数来执行学习任务。多层的神经网络可以进行特征学习, 网络在隐藏层学习输入, 随后在输出层进行分类或回归。

微测井资料是精确求取近地表模型的关键, 但是通常情况下由于成本的原因, 微测井的空间采样密度一般是在1~4 km²内布设一口井, 在一块100 km²的三维工区内, 通常微测井点有50~100个, 样本点过少很难得到好的深度学习模型, 因此在开展多元非线性三维回归模型训练过程中, 尽可能多地收集工区附近的微测井数据, 作为初始模型训练的样本, 在此基础上, 基于迁移学习原理, 利用本工区的样本点深化调优, 提高模型预测的精度。深度学习训练分为4个步骤:数据集整理、定义模型、训练/学习和预测/评估。数据集是本工区用前述谱比法估算的每个观测点的Q值函数, 由测点坐标(X, Y), 测点深度(Depth)和Q值构成, 其中X、Y和深度构成训练集的输入X_train, 对应的Q值构成训练集的标签label_train。我们需要对输入数据进行标准化, 使其满足正态分布, 开源的机器学习算法包Scikit-learn中包含的StandardScalar函数可以对任何输入数据进行标准化^[31-32]。这种数据标准化方法经过处理后数据符合标准正态分布, 即均值为0, 标准差为1, 转化函数为:

$ x=\frac{x-\mu}{\sigma} $

(23)

式中:x为输入数据; μ为均值; σ为方差或标准差。深度学习模型定义采用Keras框架的序贯模型(Sequential), 首先在模型中添加一个全连接层(Dense), 然后添加5个隐藏层, 最后一层是输出层, 神经网络结构为:3-200-200-100-100-100-1, 如图 2所示, 即输入层为3个神经元, 对应变量X、Y和深度(Depth), 5个隐藏层分别有200, 200, 100, 100, 100个神经元, 加入隐藏层可以拟合更加复杂模型, 输出层为1个神经元, 对应预测的Q值, 这是典型的多元回归网络模型设置, 用来预测单个连续值。激活函数则采用线性整流函数指代数学中的斜坡函数, 即Relu函数:f(x)=max(0, x), 如图 3所示, 加入激活函数来拟合非线性模型。模型定义后, 则需要对模型进行训练和学习。模型学习过程中, 最优化算子选用自适应矩估计(adaptive moment estimation, ADAM)算子, 损失函数选用最小平方误差函数(mean squared error, MSE)。模型训练达到精度要求后, 就形成了一个三维多元非线性回归模型, 只要输入工区内任意点的X、Y和深度, 就可以预测该点的Q值, 从而实现三维Q值建模。算法流程如图 4所示。

2 模型分析

为了便于开展微测井正演分析, 结合生产实际, 建立了如图 5所示的4层水平层状介质模型。正演模拟计算采用Futteman方程频率域Q衰减算法实现。初至波采用非零相位雷克子波, 相位为30°, 主频为40Hz。依次在图 5中的1, 2, 3三个位置激发, 地面接收。为使模拟更加符合实际情况, 在初至后面添加一个主频为60Hz, 相位为60°的子波, 和初至子波呈半叠置关系, 同时在初至地震信号添加了随机噪声。

图 6a中蓝色波形是图 5中炮点3激发地面接收得到的地震记录, 红色波形是利用CFMP方法匹配得到的初至波形。图 6b是采用CFMP方法提取的初至信号(红色)和在炮点3直接激发不考虑后续波形和噪声获得的初至波(蓝色)叠合显示。图 6c中蓝色波形是图 5中在炮点2激发地面接收得到的地震记录, 红色波形是利用CFMP方法匹配得到的初至波形。图 6d是采用CFMP方法提取的初至信号(红色)和在炮点2直接激发不考虑后续波形和噪声获得的初至波(蓝色)叠合显示。可以看到, 尽管有噪声和续至波形的影响, 采用CFMP方法提取的初至波和实际初至波基本一致, 表明该算法在初至波提取方面具有较好的抗干扰能力。

图 7a中的曲线是在炮点2(蓝色)和炮点3(红色)激发地震波地面接收的初至波通过加窗截取一段信号后初至信号求取的频谱, 受噪声和续至波的影响, 频谱呈毛刺状变化。图 7b是用这两炮的频谱计算的谱比值曲线, 可以看到, 该曲线抖动严重。图 7c中蓝虚线是谱比值曲线取对数得到的对数谱比曲线, 同样曲线抖动严重; 红实线是基于最小二乘线性拟合的结果, 进一步利用拟合斜率, 由(16)式可以求得第2层的Q值。

图 8a中的曲线是在炮点2(蓝色)和炮点3(红色)激发地面接收的地震波通过CFMP提取初至信号求取的频谱, 由于抗噪性好, 计算的频谱较为光滑。图 8b是用这两炮的频谱计算的谱比值曲线, 可以看到, 该曲线基本呈线性。图 8c中蓝虚线是谱比值曲线取对数得到的对数谱比曲线, 红实线是基于最小二乘线性拟合的结果, 基于(16)式可以求得第2层的Q值。对比分析可以看出, 与传统的加窗谱分析方法相比, 匹配追踪得到的频谱能够分离干扰信号, 得到的频谱更加平滑, 频谱分析精度更高, 同时加入正则化的谱比相较于传统的谱比法运算更加稳定。

表 1给出了模型Q值、估计Q值及其误差。从表 1可以看出, CFMP谱比法和常规谱比法相比, 抗噪能力强, Q值估算累计误差小, 精度高。

3 近地表Q值估计与三维建模 3.1 近地表Q值估计

实测微测井数据来自中国石油华北油田。探区地表平缓, 但表层结构松散且复杂多变、速度低, 对地震波的吸收、衰减严重, 地震资料分辨率低。岩石物理分析发现, 介质的吸收作用与介质孔隙度、固结程度以及孔隙充填物关系密切, 所以表层未固结岩层的吸收规律与深层固结岩层的规律不同。图 9为由近地表取样实验室测定的泥、黄土和沙3种岩样Q值与速度交会分析结果。由于表层通常是砂、泥混合的, 可以看出, 近地表的Q值为2~15, 速度为300~800m/s。

图 10给出了探区内微测井测点位置, 探区面积约60km², 区内共有42个微测井测量点。图 11是采用CFMP对数谱比法估算的图 10中不同点所在位置的Q值函数。

3.2 深度学习三维非线性Q值建模

以探区内42口微测井资料估计的近地表Q值函数中的36口井的数据构成训练集作为深度学习模型的输入样本, 即训练标签, 其余6口井作为测试集。进一步按照前述深度学习三维Q值建模原理, 通过深度神经网络的多元非线性回归, 可得到一个三维近地表Q值回归数学模型。图 12给出了图 10中6号井点微测井资料估算的Q值函数(测试集中的一个样本)与深度学习Q值回归数学模型所得预测值的对比结果, 可以看出, 预测值与微测井估算Q值之间的绝对误差值均介于±0.1内, 基本可以忽略, 说明二者具有良好的一致性。

图 13为基于深度学习多元非线性回归所得近地表Q值数学模型建立的研究区三维近地表Q值模型。图 14a和图 14b分别为10m和20m的近地表Q值等深度切片。由图 14可见, Q值空间变化规律性较强, 尽管控制样本点很少, 但是训练获得的多元非线性回归模型预测的Q值变化较为平滑, 证明了该算法在数据量较为稀疏的条件下可以很好地构建三维模型。图 14中红线是三维工区的边界, 边界线外预测的结果可靠性要差一些。

4 地震资料近地表衰减补偿效果分析

地层非弹性吸收引起的地震波能量衰减的量可通过反Q滤波方法对衰减量进行补偿校正, 以提高地震剖面的分辨率。在反Q滤波的过程中, 外推步长是地震道的时间采样率, 并要求在每一个采样点上输出结果, 这样计算量巨大, 限制了其在叠前处理过程中的应用。根据地震波在地下介质中传播的衰减规律, 可将近地表层Q值转化为等效Q值:

$ {Q_{0N}}(\tau ) = \frac{{{\tau _{N + 1}} - {\tau _1}}}{{\sum\limits_{j = 1}^N {\frac{{{\tau _{j + 1}} - {\tau _j}}}{{{Q_j}}}} }} $

(24)

式中:Q_j为近地表第j层的Q值; Q_0N(τ)表示从地表到第N层的等效Q值, 描述了地面激发点与地下相应点间的Q值关系。最后将单程旅行时概念下的等效Q值转化为双程旅行时概念下的等效Q(τ)值。基于近地表等效Q(τ)值的衰减公式表示如下:

$ \begin{array}{l} A(\tau ,\omega ) = A(0,\omega )\exp \left[ { - \int_0^\tau {\frac{\omega }{{2Q(\tau )}}} {{\left( {\frac{\omega }{{{\omega _r}}}} \right)}^{\mu (\tau )}}{\rm{d}}\tau } \right] \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \exp \left[ {{\rm{i}}\int_0^\tau {{{\left( {\frac{\omega }{{{\omega _r}}}} \right)}^{\mu (\tau )}}} \omega {\rm{d}}\tau } \right] \end{array} $

(25)

式中:μ(τ)=1/πQ(τ); A(0, ω)是地面接收的地震信号振幅; A(τ, ω)是传过等效Q层后的地震信号振幅; ω是圆频率; ω_r是参考圆频率; i是虚数单位。该方法可以对振幅和相位同时进行校正, 且等效Q值的引入更好地适应了近地表的衰减补偿处理, 更有利于提高地震资料分辨率。按上述方法可以将随深度变化的近地表层Q模型分为多个等效Q模型, 在每一等效Q层内实施反Q滤波。将波场从地面延拓到第N层的顶部, 可以递归实现, 即:

$ A\left(\tau_{0}, \omega\right) \rightarrow A\left(\tau_{1}, \omega\right) \rightarrow \cdots \rightarrow A\left(\tau_{n-1}, \omega\right) $

(26)

在近地表反Q滤波中, 则是将地面地震记录按反Q滤波的形式延拓到近地表Q值模型的底部, 在近地表模型底部以下, 给定一个较大的Q值(例如200), 以保证反Q滤波只受近地表参数的影响。

图 15为近地表Q值补偿滤波前、后的单炮记录, 显然, 近地表Q值补偿后地震记录的分辨率有明显的提高。图 16是图 15所示的单炮记录近地表Q值补偿前、后的频谱, 可以看出, 无论是1500ms以上还是以下, 高频成分都得到了有效补偿, 频带展宽明显。图 17为近地表Q值补偿前、后的叠前时间偏移剖面, 其中, 红色波形是插入的合成地震记录, 可以看出, 近地表Q值补偿后, 目标层分辨率明显提高, 与合成地震记录的一致性更好。

5 结论

本文提出了一种基于非零相位雷克子波的复数域快速匹配追踪分解并结合对数谱比法估算微测井数据近地表Q值的方法, 并基于深度学习构建了多元非线性回归算子, 实现了近地表Q值三维建模, 进而将Q值模型用于实际资料处理, 得出以下结论。

1) 采用基于非零相位雷克子波原子库改进的快速复数域匹配追踪方法, 可自动提取微测井初至波的频谱信息和初至时间, 提高了初至波频谱求取的精度和抗噪能力, 并提高了微测井资料处理的工作效率。

2) 通过引入整形正则化算子和最小平方反演方法求取谱比值, 提高了谱比算法的准确性和稳定性。

3) 将工区内不同位置的微测井数据估计得到的Q值函数作为已知样本标签, 基于迁移学习原理, 通过深度学习训练形成近地表Q值的多元非线性回归算子, 可建立较传统方法更精确的三维近地表Q值模型。

4) 华北油田实际资料应用结果表明, 本文近地表Q值估算方法具有合理性和准确性, 且高效、抗噪能力强。利用建立的三维近地表Q值模型进行近地表Q值补偿, 可展宽目的层反射波频带, 提高地震资料分辨率。

研究发现尽管可以通过匹配追踪获取精确的初至波频谱和旅行时, 但是在谱比计算中, 需要选取计算谱比值的频带宽度参数, 而估算的Q值对谱比值比较敏感, 因此往往会造成估算结果的不稳定。可以采用理论精度相对差一点的质心频移法代替谱比法进行Q值估算, 以提高算法的稳定性。