一种基于广义地震子波的地层Q值计算方法

摘要

本发明涉及一种基于广义地震子波的地层Q值计算方法，包括以下步骤：1)对采集到的零偏VSP资料进行预处理，进行波场分离得到VSP下行波记录，然后利用窗函数提取每一道的直达波；2)以井中相邻或者预先设定的检波器的距离，将地下介质划分为多个地层；3)计算每一个地层的Q值。本发明基于广义地震子波假设，从粘弹介质中单程波传播理论出发，推导得到了一种利用子波包络峰值处瞬时频率的变化来计算地层Q值的新方法，由于广义地震子波具有参考频率和分数阶系数两个灵活的参数，可以选取最佳的参数来匹配实际地震子波的振幅谱，从而使得估计方法具有良好的子波适应性，并且提高了Q值估计的精度。

权利要求

1.一种基于广义地震子波的地层的品质因子Q值计算方法，包括以下步骤：

1)对采集到的零偏VSP资料进行预处理，利用中值滤波进行波场分离得到VSP下行波记录，然后利用窗函数提取每一道的直达波；

2)以井中相邻或者预先设定的检波器的距离，将地下介质划分为多个地层；

3)计算每一个地层的品质因子Q值，单个地层的品质因子Q值的计算方法如下：①对第n个地层的顶部和底部接收到的直达波，分别利用广义地震子波进行拟合，得到顶部直达波的分数阶系数 和参考频率 底部直达波的分数阶系数 和参考频率②对第n个地层顶部和底部接收到的信号，分别计算其包络，拾取包络峰值处的时间位置，分别为 和 从而地震波在此地层中的传播时间为：③对第n个地层顶部和底部接收到的信号，分别计算其瞬时频率和幅度谱加权平均频率，分别拾取包络峰值处 和 时刻对应的瞬时频率值，记为 和④计算第n个地层的品质因子Qn为：

其中，L(u)是关于广义地震子波分数阶系数的函数，其表达式为：式中的Γ()是Gamma函数；

衰减前后的包络峰值处瞬时频率为

均值参考频率 和均值分数阶系数 分别为顶部和底部相应参数的均值：

2.如权利要求1所述的一种基于广义地震子波的地层的品质因子Q值计算方法，其特征在于：所述步骤3)的④中，地层的品质因子Q值的计算公式的建立过程如下：根据粘弹介质中单程波的传播理论，对于在品质因子为Q的衰减介质中传播的震源子波，设其振幅谱为S(f)，经过时间τ，接收到的子波振幅谱A(f,τ)满足其中，f是频率；用广义地震子波来表示归一化的震源子波振幅谱：式中，u是广义地震子波的分数阶系数；f0是广义地震子波的参考频率；

广义地震子波包络峰值处的瞬时频率等于它的幅度谱加权平均频率：其中，Γ()是Gamma函数：

式中的x表示函数自变量，t表示积分变量；

衰减后的振幅谱A(f,τ)改写为

衰减 后 子波的 包 络峰 值 处瞬 时 频率 等 于其 幅度 谱 加权 平 均频 率在上式中，作变量替换，令则有

在上式中， 是一个小量，引入如下的近似

从而得到

从而，衰减后子波的包络峰值处瞬时频率为：对上式中的分子分母同时除以 并引入近似

衰减后子波的包络峰值处瞬时频率进一步化简为：注意到衰减前子波的包络峰值处瞬时频率为

那么，衰减后子波的包络峰值处瞬时频率改写为反过来，由经过传播时间τ前后，子波包络峰值处的瞬时频率来估计地层的品质因子Q值为：其中，系数函数L(u)是关于广义地震子波分数阶系数的函数，表示为：应用时，用广义地震子波分别来拟合衰减前子波与衰减后的子波，但一般地，衰减前后子波的参考频率和分数阶系数并不相同，因此，对衰减前后子波的参考频率和分数阶系数求均值：其中， 和 分别是衰减前后子波的参考频率，us和ur分别是衰减前后子波的分数阶系数；从而，地层的品质因子Q值计算公式修正为：

3.如权利要求2所述的一种基于广义地震子波的地层的品质因子Q值计算方法，其特征在于：所述步骤3)的①中，利用广义地震子波拟合得到地层顶部直达波和底部直达波的分数阶系数、参考频率的过程如下：对待拟合的实际地震子波，利用快速傅里叶变换计算其幅度谱为|H(f)|，则功率谱为|H(f)|2，利用功率谱计算得到子波的均值频率fm和偏差fσ：理论上，广义地震子波的均值频率为

那么

其中，与Gamma函数相关的项展开为

从而

这是一个关于 的多项式求根问题，用迭代法求得最优的u；进一步计算得到广义地震子波的参考频率：对于第n个地层而言，求得的衰减前的分数阶系数与参考频率即为地层顶部直达波的分数阶系数 和参考频率 求得的衰减后的分数阶系数与参考频率即为地层底部直达波的分数阶系数 和参考频率

说明书

一种基于广义地震子波的地层Q值计算方法

技术领域

[0001] 本发明涉及一种基于广义地震子波的地层Q值计算方法，属于地震勘探技术领域。

背景技术

[0002] 地震波在地下吸收介质中的传播会经历衰减和频散，这种粘弹性效应一般由品质因子Q来定量表示。在地震勘探中，准确估计Q值是非常重要的。Q可以作为烃类检测和储层描述的工具，也可以作为反Q滤波，非平稳反褶积或逆时偏移中实现吸收补偿的重要参数。

[0003] 在实际中，Q常常假设为在地震频段内是与频率无关的，也就是所谓的常Q模型。当地震波在吸收介质中传播，其高频分量较低频分量衰减得快。因此，可以通过地震波的频率成分的变化来计算地震衰减，如对数谱比法(Logarithm spectral ratio,LSR)、中心频率偏移法(Centroid frequency shift,CFS)和峰值频率偏移法(Peak frequency shift,PFS)等。

[0004] 对数谱比法对震源子波没有先验假设，但是其对噪音特别敏感，需要选择合适的带宽范围来拟合斜率。中心频率偏移法假设震源子波振幅谱是高斯形状的，由于其使用幅度谱的统计特征来估计衰减，因而估计结果较稳定，但如果震源子波不能满足高斯形状，那么此方法有较大的系统误差。峰值频率偏移法则假设震源子波是Ricker子波，由于振幅谱峰值位置容易受到噪音影响，因而峰值频率偏移法不够稳健。对CFS和PFS而言，其基于简单的震源子波假设可能导致不适应性问题。

[0005] Yang和Gao提出了一种利用包络峰值处瞬时频率(Envelope peak instantaneous frequency，EPIF)偏移计算Q值的方法，在计算时可以避免加时窗问题，提高了Q值估计的纵向分辨率。然而，此方法的前提是假设子波谱可以用一个常相位的高斯函数表示。实际上，子波谱常常并不是高斯的，而且随着波在吸收介质中的传播，子波谱的形状也会被显著的改变。

[0006] 综上，需要具有更好的子波适应性的Q值计算方法，从而提高Q值计算精度。

发明内容

[0007] 针对上述问题，本发明的目的是提供一种基于广义地震子波的地层Q值计算方法，该方法具有更好的子波适应性，从而使得Q值计算精度更高。

[0008] 为实现上述目的，本发明采用以下技术方案：一种基于广义地震子波的地层Q值计算方法，包括以下步骤：

[0009] 1)对采集到的零偏VSP资料进行预处理，利用中值滤波进行波场分离得到VSP下行波记录，然后利用窗函数提取每一道的直达波；

[0010] 2)以井中相邻或者预先设定的检波器的距离，将地下介质划分为多个地层；

[0011] 3)计算每一个地层的Q值，单个地层Q值的计算方法如下：

[0012] ①对第n个地层的顶部和底部接收到的直达波，分别利用广义地震子波进行拟合，得到顶部直达波的分数阶系数 和参考频率 底部直达波的分数阶系数 和参考频率[0013] ②对第n个地层顶部和底部接收到的信号，分别计算其包络，拾取包络峰值处的时间位置，分别为 和 从而地震波在此地层中的传播时间为：

[0014] ③对第n个地层顶部和底部接收到的信号，分别计算其瞬时频率或者加权瞬时频率，分别拾取包络峰值处 和 时刻对应的瞬时频率值，记为 和

[0015] ④计算第n个地层的品质因子Qn为：

[0016]

[0017] 其中，L(u)是关于广义地震子波分数阶系数的函数，其表达式为：

[0018]

[0019] 式中的Γ()是Gamma函数；

[0020] 衰减前后的包络峰值处瞬时频率为

[0021] 均值参考频率 和均值分数阶系数 分别为顶部和底部相应参数的均值：

[0022]

[0023] 进一步地，所述步骤3)的④中，Q值的计算公式的建立过程如下：

[0024] 根据粘弹介质中单程波的传播理论，对于在品质因子为Q的衰减介质中传播的震源子波，设其振幅谱为S(f)，经过时间τ，接收到的子波振幅谱A(f,τ)满足[0025]

[0026] 其中，f是频率；用广义地震子波来表示归一化的震源子波振幅谱：

[0027]

[0028] 式中，u是广义地震子波的分数阶系数；f0是广义地震子波的参考频率；

[0029] 广义地震子波包络峰值处的瞬时频率等于它的幅度谱加权平均频率：

[0030]

[0031] 其中，Γ()是Gamma函数：

[0032]

[0033] 式中的x表示函数自变量，t表示积分变量；

[0034] 衰减后的振幅谱A(f,τ)改写为

[0035]

[0036] 衰减后子波的包络峰值处瞬时频率等于其幅度谱加权的平均频率[0037]

[0038] 在上式中，作变量替换，令

[0039]

[0040] 则有

[0041]

[0042] 在上式中， 可以认为是一个小量，引入如下的近似

[0043]

[0044]

[0045] 从而得到

[0046]

[0047]

[0048] 从而，衰减后子波的包络峰值处瞬时频率为：

[0049]

[0050] 对上式中的分子分母同时除以 并引入近似

[0051]

[0052] 衰减后子波的包络峰值处瞬时频率可以进一步化简为：

[0053]

[0054] 注意到衰减前子波的包络峰值处瞬时频率为

[0055]

[0056] 那么，衰减后子波的包络峰值处瞬时频率可以改写为

[0057]

[0058] 反过来，可以由经过传播时间τ前后，子波包络峰值处的瞬时频率来估计地层的品质因子Q为：

[0059]

[0060] 其中，系数函数L(u)是关于广义地震子波分数阶系数的函数，表示为：

[0061]

[0062] 应用时，可用广义地震子波分别来拟合衰减前子波与衰减后的子波，但一般地，衰减前后子波的参考频率和分数阶系数并不相同，因此，对衰减前后子波的参考频率和分数阶系数求均值：

[0063]

[0064] 其中， 和 分别是衰减前后子波的参考频率，us和ur分别是衰减前后子波的分数阶系数；从而，地层品质因子Q计算公式修正为：

[0065]

[0066] 进一步地，所述步骤3)的①中，利用广义地震子波拟合得到地层顶部直达波和底部直达波的分数阶系数、参考频率的过程如下：

[0067] 对待拟合的实际地震子波，利用快速傅里叶变换计算其幅度谱为|H(f)|，则功率谱为|H(f)|2，利用功率谱计算得到子波的均值频率fm和偏差fσ：

[0068]

[0069]

[0070] 理论上，广义地震子波的均值频率为

[0071]

[0072]

[0073] 那么

[0074]

[0075] 其中，与Gamma函数相关的项可以展开为

[0076]

[0077] 从而

[0078]

[0079] 这是一个关于 的多项式求根问题，可以用迭代法求得最优的u；进一步，可以计算得到广义地震子波的参考频率：

[0080]

[0081] 对于第n个地层而言，求得的衰减前的分数阶系数与参考频率即为地层顶部直达波的分数阶系数 和参考频率 求得的衰减后的分数阶系数与参考频率即为地层底部直达波的分数阶系数 和参考频率

[0082] 本发明由于采取以上技术方案，其具有以下优点：本发明基于广义地震子波假设，从粘弹介质中单程波传播理论出发，推导得到了一种利用子波包络峰值处瞬时频率的变化来计算地层Q值的新方法，由于广义地震子波具有参考频率和分数阶系数两个灵活的参数，可以选取最佳的参数来匹配实际地震子波的振幅谱，从而使得估计方法具有良好的子波适应性，并且提高了Q值估计的精度。

[0083] 计算得到的地层Q值，将有利于含气性预测，或者用于衰减效应的补偿。

附图说明

[0084] 图1是本发明方法的流程示意图；

[0085] 图2是零偏VSP资料地层划分示意图；

[0086] 图3是震源子波分数阶系数变化时，分别采用本发明方法、包络峰值频率法、中心频率偏移法和峰值频率偏移法的Q值计算结果的相对误差的变化曲线比较图；

[0087] 图4是实际零偏VSP资料的下行波与直达波；其中，图(a)是下行波，图(b)是直达波；

[0088] 图5是本发明方法计算得到的地层Q值曲线。

具体实施方式

[0089] 下面结合附图和实施例对本发明进行详细的描述。

[0090] 如图1所示，本发明提出了一种基于广义地震子波的地层Q值计算方法，包括以下步骤：

[0091] 1)对采集到的零偏VSP(Vertical Seismic Profiling，垂直地震剖面)资料进行预处理，利用中值滤波进行波场分离得到VSP下行波记录，然后利用窗函数提取每一道的直达波；

[0092] 2)以井中相邻或者预先设定的检波器的距离，将地下介质划分为多个地层(如图2所示)；

[0093] 3)计算每一个地层的Q值，单个地层Q值的计算方法如下：

[0094] ①对第n个地层的顶部和底部接收到的直达波，分别利用广义地震子波进行拟合，得到顶部直达波的分数阶系数 和参考频率 底部直达波的分数阶系数 和参考频率[0095] ②对第n个地层顶部和底部接收到的信号，分别计算其包络，拾取包络峰值处的时间位置，分别为 和 从而地震波在此地层中的传播时间为：

[0096] ③对第n个地层顶部和底部接收到的信号，分别计算其瞬时频率或者加权瞬时频率，分别拾取包络峰值处 和 时刻对应的瞬时频率值，记为 和

[0097] ④计算第n个地层的品质因子Qn为：

[0098]

[0099] 其中，L(u)是关于广义地震子波分数阶系数的函数，其表达式为：

[0100]

[0101] 式中的Γ()是Gamma函数；

[0102] 衰减前后的包络峰值处瞬时频率为

[0103] 均值参考频率 和均值分数阶系数 分别为顶部和底部相应参数的均值：

[0104]

[0105] 上述实施例中，步骤3)的④中，Q值的计算公式的建立过程如下：

[0106] 根据粘弹介质中单程波的传播理论，对于在品质因子为Q的衰减介质中传播的震源子波，设其振幅谱为S(f)，经过时间τ，接收到的子波振幅谱A(f,τ)满足[0107]

[0108] 其中，f是频率。假设归一化的震源子波振幅谱可以用广义地震子波来表示：

[0109]

[0110] 式中，u是广义地震子波的分数阶系数；f0是广义地震子波的参考频率。

[0111] 根据Barnes，广义地震子波包络峰值处的瞬时频率等于它的幅度谱加权平均频率：

[0112]

[0113] 其中，Γ()是Gamma函数：

[0114]

[0115] 式中的x表示函数自变量，t表示积分变量；

[0116] 衰减后的振幅谱A(f,τ)改写为

[0117]

[0118] 衰减后子波的包络峰值处瞬时频率等于其幅度谱加权的平均频率[0119]

[0120] 在上式中，作变量替换，令

[0121]

[0122] 则有

[0123]

[0124] 在上式中， 可以认为是一个小量，引入如下的近似

[0125]

[0126]

[0127] 从而得到

[0128]

[0129]

[0130] 从而，衰减后子波的包络峰值处瞬时频率为：

[0131]

[0132] 对上式中的分子分母同时除以 并引入近似

[0133]

[0134] 衰减后子波的包络峰值处瞬时频率可以进一步化简为：

[0135]

[0136] 注意到衰减前子波的包络峰值处瞬时频率为

[0137]

[0138] 那么，衰减后子波的包络峰值处瞬时频率可以改写为

[0139]

[0140] 反过来，可以由经过传播时间τ前后，子波包络峰值处的瞬时频率来估计地层的品质因子Q为：

[0141]

[0142] 其中，系数函数L(u)是关于广义地震子波分数阶系数的函数，表示为：

[0143]

[0144] 实际中，可以用广义地震子波分别来拟合衰减前子波与衰减后的子波，但一般地，衰减前后子波的参考频率和分数阶系数并不相同。因此，可以对衰减前后子波的参考频率和分数阶系数求均值：

[0145]

[0146] 其中， 和 分别是衰减前后子波的参考频率，us和ur分别是衰减前后子波的分数阶系数。从而，地层品质因子Q计算公式修正为：

[0147]

[0148] 通过上述公式可知，为了估计地层Q值，需要已知衰减前后的子波包络峰值处的瞬时频率fep(0)和fep(τ)，旅行时间τ，以及用广义地震子波拟合衰减前后子波得到的平均参考频率和平均分数阶系数。

[0149] 上述实施例中，所述步骤3)的①中，利用广义地震子波拟合得到地层顶部直达波和底部直达波的分数阶系数、参考频率的过程如下：

[0150] 对待拟合的实际地震子波，利用快速傅里叶变换计算其幅度谱为|H(f)|，则功率谱为|H(f)|2。利用功率谱计算得到子波的均值频率fm和偏差fσ：

[0151]

[0152]

[0153] 理论上，广义地震子波的均值频率应为

[0154]

[0155]

[0156] 那么

[0157]

[0158] 其中，与Gamma函数相关的项可以展开为

[0159]

[0160] 从而

[0161]

[0162] 这是一个关于 的多项式求根问题，可以用迭代法求得最优的u。进一步，可以计算得到广义地震子波的参考频率：

[0163]

[0164] 利用上述方法，可以分别求得衰减前后的分数阶系数与参考频率。对于第n个地层而言，求得的衰减前的分数阶系数与参考频率即为地层顶部直达波的分数阶系数 和参考频率 求得的衰减后的分数阶系数与参考频率即为地层底部直达波的分数阶系数 和参考频率

[0165] 上述实施例中，步骤3)的②中，计算信号的包络的过程如下：

[0166] 对于待分析的实地震道x(t)，其对应的复地震道为：

[0167] z(t)＝x(t)+iy(t)＝a(t)exp[iθ(t)]

[0168] 其中， 是正交道， 表示Hilbert变换，θ(t)是瞬时相位，定义为[0169]

[0170] a(t)是信号的包络，也称为瞬时振幅，定义为

[0171]

[0172] 根据上述方法，可以分别计算衰减前后子波的包络。对于第n个地层而言，求得的衰减前子波的包络即为地层顶部接收到的信号包络，求得的衰减后子波的包络即为地层底部接收到的信号包络。

[0173] 上述实施例中，步骤3)的③中，对接收到的信号计算顺时频率、加权顺时频率的过程如下：

[0174] 对待分析的信号，瞬时频率定义为

[0175]

[0176] 其中，ε2是阻尼因子，选取为a2(t)最大值的0.00001～0.001倍。为进一步提高瞬时频率的抗噪性能，可以利用包络的平方作为权系数，计算加权瞬时频率：

[0177]

[0178] 其中W(t)＝a2(t)是加权系数，T是加权窗口的大小。

[0179] 下面利用合成数据测试本发明方法(GSW)的技术效果：

[0180] 设地层品质因子为Q＝50，地层旅行时为0.3秒。考虑当震源子波变化时，各种Q值计算方法估计结果的准确性，其中用相对误差来度量准确性。假设震源子波的参考频率为60Hz，分数阶系数从1到2.5变化，变化的步长为0.1，也即有16种情况。分别利用本发明的方法(GSW)，包络峰值频率法(EPIF)，中心频率偏移法(CFS)和峰值频率偏移法(PFS)计算Q值，相对误差(Relative error)随分数阶系数(Fractional value)的变化曲线如图3所示。从图中可以看到，对于EPIF和CFS，它们假设要求震源子波振幅谱是高斯的，因此估计结果存在明显的系统误差。对于PFS，它假设要求振幅谱是Ricker子波，注意到当分数阶系数为u＝
2时广义地震子波等价于Ricker子波，此时PFS方法的相对误差很小，但分数阶系数为其他值时，PFS计算结果的误差较大。本发明提出的方法具有最小的相对误差，表明方法的准确性和子波适应性。

[0181] 本发明的方法用于实际零偏VSP资料的地层Q值计算：

[0182] 此零偏VSP资料由炸药震源激发，在深度1500-3480米共记录了100道数据，道间距为20米，数据的采样时间步长是1毫秒。图4展示的是实际零偏VSP数据的下行波和拾取的直达波。利用本发明提出的方法计算地层Q值，得到的Q值曲线如图5所示。在3395-3460米范围(图中箭头标记)的低Q值区域，在实际钻井中显示为含气砂岩地层。因此，利用本发明提出的方法计算得到地层Q值，可以帮助解释人员进行含气性预测。

[0183] 上述各实施例仅用于说明本发明，其中方法的实施步骤等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。