CN102478670B - 一种利用地震衰减属性进行储层流体性质预测的方法 - Google Patents

Abstract

本发明实施例提供一种利用地震衰减属性进行储层流体性质预测的方法，所述方法包括：确定储层饱和度与地层品质因子关系，所述储层饱和度与地层品质因子关系包括储层饱和度与地层品质因子变化曲线、衰减峰值出现的饱和度位置；对地震成果数据进行储层反射界面构造解释，确定储层在地震剖面中的反射时间位置；确定储层品质因子的大小；沿反射界面对地层品质因子进行平滑与处理，对地震品质因子数据采用微分计算方法确定微分衰减属性大小；利用储层饱和度与地层品质因子变化曲线分析，确定储层流体性质分布特征。本发明实施例可以提高地震储层性质预测的可靠性和精度。

Description

一种利用地震衰减属性进行储层流体性质预测的方法

技术领域

本发明涉及石油地球物理勘探地震储层预测方法，尤其涉及一种利用地震衰减属性进行储层流体性质预测的方法。

背景技术

利用地震波衰减属性进行储层流体性质检测已成为近年来研究的热点。当地震波在含流体地层中传播时，波动能量会发生衰减。引起地震波能量衰减的主要因素包括地层界面散射引起的能量衰减和因地层本身(地层岩性、孔隙及其充填流体性质等)性质差异引起的衰减。通常将后者称之为地层本征衰减。因地层界面散射引起的衰减，例如：波在地层中的微曲多次、几何扩散等，在一定平面范围内其作用横向变化不大。地层本征衰减是我们利用地震反射记录的衰减属性进行储层流体预测的主要依据。

一般而言，地层中流体性质不同，其对地震波传播能量衰减程度不一样，含气地层吸收地震波能量较其它地层强，含油水地层对地震波能量的衰减作用相对较弱。通常采用地层品质因子来度量地层本征衰减大小。含流体地层品质因子的大小主要与地层岩石性质、孔渗条件以及孔隙流体性质密切相关。

利用地震反射记录的衰减属性进行地层流体性质检测主要涉及由地震记录提取地层品质因子(或吸收系数)、提取参数的处理和解释、地层流体性质检测检测分析。目前常用的主要方法包括：(1)地层品质因子或吸收系数提取与分析：由地震记录提取品质因子的方法主要包括谱比法、振幅衰减法、解析信号法、信号(子波、相位、频率等)模拟法、脉冲振幅法、上升时间法等，通过计算地层品质因子，采用定性的方法比较品质因子相对大小的分布特征，进行地层流体性质推测；(2)地震记录频谱分析：即采用傅里叶变换或小波变换等，计算得到(储层时间段)地震记录频谱，利用谱中能量分布特征、谱的几何形态，定性确定地震记录的衰减差异变化(高频成分衰减多，则地层品质因子相对较小)，进而利用衰减差异特征推测储层流体性质；(3)时频分析技术，即采用数学变换或匹配追踪等方法将地震记录分解为不同谐波分量记录，比较不同分量记录中能量变化关系，推测储层流体性质。

现行利用地震衰减属性进行储层流体预测方法存在重要不足，即利用地震波衰减属性进行储层流体性质预测的结果中存在许多不确定性因素。其主要原因有：(1)地层本征品质因子具有明确的物理意义，即品质因子与地层性质、孔渗参数以及孔隙流体性质有关。但由地震记录计算品质因子其精度受到地震记录质量的约束。(2)使用地震记录谱信息或时频分析，因其缺乏明确的物理关系，这些属性参数仅能起到粗略推测地层性质的作用。(3)实际储层中完全饱和一种性质流体的情况非常少见，实际储层中多以两相或多相态流体共存方式赋存在岩石孔隙中。(4)由地震记录计算得到品质因子与储层性质之间缺乏定量标定，是导致地震储层流体检测结果多解性的重要因素。

发明内容

本发明实施例提供一种利用地震衰减属性进行储层流体性质预测的方法，以提高地震储层性质预测的可靠性和精度。

一方面，本发明实施例提供了一种利用地震衰减属性进行储层流体性质预测的方法，所述方法包括：确定储层饱和度与地层品质因子关系，所述储层饱和度与地层品质因子关系包括储层饱和度与地层品质因子变化曲线、衰减峰值出现的饱和度位置；对地震成果数据进行储层反射界面构造解释，确定储层在地震剖面中的反射时间位置；确定储层品质因子的大小；沿反射界面对地层品质因子进行平滑与处理，对地震品质因子数据采用微分计算方法确定微分衰减属性大小；利用储层饱和度与地层品质因子变化曲线分析，确定储层流体性质分布特征。

上述技术方案具有如下有益效果：因为采用利用地层品质因子与衰减峰值饱和度特征，通过提取地震记录衰减参数，采用求极值方法寻找地震数据体中衰减峰值位置，并对衰减峰值进行标定，进而实现利用地震衰减属性进行储层流体性质预测的方法，相对于现有技术，可以提高地震储层性质预测的可靠性和精度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例一种利用地震衰减属性进行储层流体性质预测的方法流程图；

图2为本发明实施例参考样品与岩石样品声波测试信号图；

图3为本发明实施例参考样品与岩石样品振幅谱示意图；

图4为本发明实施例振幅谱比法计算地层品质因子示意图；

图5为本发明实施例地层品质因子与饱和度关系曲线图；

图6为本发明实施例含油储层地质剖面图；

图7为本发明实施例含油储层地震记录剖面图；

图8为本发明实施例地震品质因子剖面图；

图9为本发明实施例含油储层地震微分衰减属性参数剖面；

图10为本发明实施例各地层单元地层衰减模型示意图；

图11为本发明实施例地层衰减正演模型记录示意图；

图12为本发明实施例含油砂层在地震品质因子剖面中特征示意图；

图13为本发明实施例微分衰减属性参数剖面图；

图14为本发明实施例衰减因子预测有利含油砂层分布区域示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

大量的实验测试结果证明：地层品质因子不仅与流体性质有关，而且与流体饱和度有关。一般，地层品质因子随着含油饱和度增加逐渐减小(衰减程度增强)，并在即将完全饱和时出现衰减最小值-我们称之为衰减峰值饱和特征参数；在此饱和度之后，流体饱和度继续增加，反而使得地震波能量衰减急剧减小。为了克服现行地震衰减属性进行储层流体性质预测的不足，提高地震储层性质预测的可靠性和精度，利用实际储层衰减峰值饱和特征参数，本发明实施例提出了基于储层岩石物理的地震衰减峰值进行储层流体性质确定性分析新方法。

为实现上述目的，如图1所示，为本发明实施例一种利用地震衰减属性进行储层流体性质预测的方法流程图，所述方法包括：

101、确定储层饱和度与地层品质因子关系，所述储层饱和度与地层品质因子关系包括储层饱和度与地层品质因子变化曲线、衰减峰值出现的饱和度位置；

102、对地震成果数据进行储层反射界面构造解释，确定储层在地震剖面中的反射时间位置；

103、确定储层品质因子的大小；

104、沿反射界面对地层品质因子进行平滑与处理，对地震品质因子数据采用微分计算方法确定微分衰减属性大小；

105、利用储层饱和度与地层品质因子变化曲线分析，确定储层流体性质分布特征。

可选的，通过对储层进行试验分析，确定储层饱和度与地层品质因子关系。

可选的，对储层反射界面的地震记录进行傅里叶变换，确定储层品质因子的大小。

可选的，对储层反射界面的地震记录进行傅里叶变换，采用振幅谱比法确定储层品质因子的大小。

可选的，所述微分衰减属性为对地震品质因子进行微分运算所产生的微分品质因子参数。

可选的，所述利用储层饱和度与地层品质因子变化曲线分析，确定储层流体性质分布特征，包括：利用已知钻测井确定的地层信息，对微分衰减参数进行标定，结合储层饱和度与地层品质因子变化曲线分析，确定储层流体性质分布特征。

可选的，所述利用储层饱和度与地层品质因子变化曲线分析，确定储层流体性质分布特征，包括：确定地震品质因子的极值点位置，或微分衰减属性的最小绝对值位置；根据所述极值点位置或所述最小绝对值位置附近微分衰减参数的变化特征，微分衰减参数缓慢变化一侧为低含流体饱和度方向，微分衰减参数变化急剧一侧为接近完全饱和流体的一侧。

本发明实施例中采取的技术路线是通过对储层进行试验分析，确定储层饱和度与地层品质因子关系，特别是衰减随饱和度变化关系曲线、衰减峰值出现的饱和度位置；对地震成果数据进行必要储层反射界面构造解释，确定储层在地震剖面中的反射时间位置；对储层反射界面的较小时窗内地震记录进行傅里叶变换，采用振幅谱比法确定储层品质因子的大小；沿反射界面对地层品质因子进行适当平滑与处理，对品质因子数据采用微分计算方法确定微分衰减属性(即对地震品质因子进行微分运算所产生的微分品质因子参数)大小，最后利用已知井信息标定，结合储层饱和度与品质因子变化曲线分析，确定储层流体性质分布特征。本发明所采用的主要技术手段如下：

(1)含流体地层衰减特征分析

利用地震波衰减属性进行储层流体性质检测需要了解含流体储层品质因子变化规律，确定所研究储层地层吸收与流体性质及其饱和度之间关系，为地震储层流体检测技术提供必要的基础依据。地震波的衰减与地层中孔隙流体性质、粘滞性以及内摩擦性质有关。含流体地层具有较大的粘滞性和内摩擦性。地震波在穿过含流体地层时，地震波高频能量衰减较快，低频能量衰减相对较慢。因此，在含流体地层位置地震波出现能量显著减少。

引起地震波衰减的主要因素包括：地层性质、孔隙及其流通性质。当地层性质变化不大时，一般地层孔隙度大小也比较稳定。但孔隙流体性质改变或饱和度发生变化，对地震波能量衰减有重要影响，其作用主要表现为流体粘滞性和内摩擦引起的高频能量损失。地层品质因子是度量含流体地层中地震波能量损耗大小的物理量，该参数的大小依赖于波动频率、应变幅度、地层性质(特别是孔隙流体状态)等物理参数。

通常，确定储层品质因子与孔隙流体性质关系的方法主要以实验测试为主。具体方法包括实验样品制备、声波特性测试以及地层品质计算等。

(a)实验样品制备：选择储层段具有代表性岩石，将其加工成直径为25.4mm，长度为50.0mm的圆柱形标准样，样品两端需用研磨仪打磨，使之光滑、两端保持高度平整，端面与圆柱轴线垂直。采用高精度计量器确定样品几何尺寸，精度保持在0.001mm以内。根据样品尺寸，采用均质的金属铝材料，制作一个与岩石样品尺寸完全相同的铝样品，将其作为参考样品。

进行测试前，先将岩石样品放入烘干箱处理48h，并在干燥箱保持24h，然后取出样品，使用高精度电子天平秤(精度为0.001g)测定样品在干燥条件下的重量，并记录样品重量。

将样品放置在盛有油(或水)的真空容器中进行抽真空饱和，真空度大约为0.001τ，真空保持时间在24h以上，完全饱和流体样品在容器中看不到气泡从样品中逸出。取出完全饱和样品后，即刻称重，确定完全饱和样品的重量；利用干燥样品重量和完全饱和样品重量之差dW，并按下面方程确定样品总体积。

V_sat＝dW·ρ_f

其中ρ_f是孔隙流体的密度。

(b)声波特性测试

将样品放置在岩芯夹持器中，使岩芯受到岩芯夹持器中围压和孔隙压力的作用，以模拟实际地层条件。同时，采用不同流体作为孔隙压力介质(如气体等)，改变样品流体饱和度变化，具体做法如下：通过精确控制孔隙压端注入流体量和出口端驱替出来的原孔隙流体重量，确定饱和度状态下岩石样品的重量，并利用干燥样品重量和部分饱和样品重量之差dW，计算岩石样品在不同时刻的饱和度，计算方法如下。

记部分饱和状态下，岩石体积为：V_part＝dW′·ρ_f

则岩石样品饱和度为：

S_fluid＝V_part/V_sat＝dW′/dW

同时，在每一个饱和度点状态下，采用脉冲透射法测试样品声波特性，测试方法简述如下：在某一个稳定的饱和度状态下，利用岩芯夹持器中的声波探头，通过发射电脉冲激励加载在样品一端的声波探头压电陶瓷震源产生振动信号，该振动信号由样品的一端传播至另一段段，并被加载在样品另一端的检测探头记录，最终所记录的信号转换为离散数字信号，并存储在磁盘介质中。

通过对每一个饱和度状态下样品进行声学参数测试，获得岩石样品在不同饱和度条件下声波传播曲线。

采用脉冲透射法对参考金属铝样品进行声学参数测试，并记录铝样品的声波传播记录。

(c)地层品质因子计算

地层品质因子的定义为：

$Q=\frac{\mathrm{\ωE}}{-\mathrm{dE}/\mathrm{dt}}=\frac{2\mathrm{\πW}}{\mathrm{\ΔW}}$

其中，E为所研究系统中的瞬时能量，-dE/dt是能量损失的变化率，W为介质在一个完整的应力-应变周期内所存储的最大能量，ΔW为每个振动周期内波动能量损失大小。

对于声波传播，其能量损耗常用衰减系数表示：

A＝A₀exp(-αx)

其中，A为声波传播了距离x后的能量，A₀为原始(未衰减)声波能量；x为声波传播距离；α为衰减系数，其大小不仅与波动频率有关，而且与介质品质因子和速度相关，衰减系数与品质因子关系可表示为：

$\mathrm{\α}=\frac{\mathrm{\πf}}{\mathrm{QV}}$

其中，f为声波振动频率；Q为地层品质因子；V为地层速度。

利用实验测试方法确定地层品质因子的方法有很多，本方法采用谱比法由声波实验记录确定地层品质因子。与其它方法相比，该方法具有实验计算结果稳定、计算精度高等特点。具体方法如下：

假定由激发探头产生的震源为一平面波S(ω)，该平面波由样品一端传播至另一端。接收探头记录到的振幅谱为Ar(ω)。平面波衰减方程可表示为：

A_r(ω)＝S(ω)exp(-α(ω)L)

其中，L为样品的长度；α平面波传播过程中能量衰减系数。

通常，采用脉冲透射法进行声波测量时，所接受的信号与样品几何尺寸(样品长度与半径)、波的传播频率或波长有关以及传播介质与声波探头直接的透射系数有关。假定，样品与探头之间透射系数大小为T，样品几何尺寸与性质对波形的影响为G(λ，r，L)，则声波测量中发射探头能量S(ω)与接收探头能量Ar(ω)之间存在如下关系：

A_r(ω)＝S(ω)·T₁·T₂·G(λ，r，L)·exp(-α(ω)L)

其中，T1表示探头与样品接触面的透射系数；T2表示样品与探头接触面的透射系数；

实际上存在T₁＝1-R₁，T₂＝1-R₂，且有R₁＝-R₂，故有：

T₁·T₂＝1-R²

对于不同饱和度点情况下岩石样品声波曲线的振幅与频率关系可表示为：

$A_r^s\left(\mathrm{\ω}\right)=S\left(\mathrm{\ω}\right)\·(1-R_s^2){\·G}_s(\mathrm{\λ},r,L)\·\exp (-{\mathrm{\α}}_s\left(\mathrm{\ω}\right)L)$

其中，

为被测样品所接收的声波振幅谱；

为样品与探头之间反射系数，可根据声波速度计算得到；G_s(λ，r，L)为样品几何衍射因子；α_s(ω)所求衰减系数。

由该方程可知，欲求衰减系数需要确定因样品G_s(λ，r，L)参数，该参数大小仅是样品长度、半径等几何尺寸的函数，直接定量计算G_s(λ，r，L)非常困难。

铝样品的品质因子非常大(衰减系数非常小，可以近似为零)，通过制作一个几何尺寸与实际样品完全一样的参考样品，并采集其传播振幅谱。

$A_r^a\left(\mathrm{\ω}\right)=S\left(\mathrm{\ω}\right){\·G}_a(\mathrm{\λ},r,L)\·\exp (-{\mathrm{\α}}_a\left(\mathrm{\ω}\right)L)$

其中，

为金属铝样品所接收的声波振幅谱；因声波探头由铝制作，故铝样品与探头之间反射系数为0；G_a(λ，r，L)为样品几何衍射因子；α_a(ω)所铝样品衰减系数(近似为零)。

通过使用岩石样品振幅谱与铝样品振幅相比方法，可以消除样品几何尺寸对接收声波振幅谱的影响，进而利用振幅信息计算得到声学衰减系数，即将铝样品振幅谱与岩石样品振幅谱相除，两边分别取对数，整理得：

${\mathrm{\α}}_s\left(\mathrm{\ω}\right)=\frac{20}{L}{\log }_{10}[\frac{A_r^a\left(\mathrm{\ω}\right)}{A_r^s\left(\mathrm{\ω}\right)}\·\left({1-R}^2\right)]-\frac{20}{L}{\log }_{10}\left[\frac{G_a(\mathrm{\λ},r,L)}{G_s(\mathrm{\λ},r,L)}\right]$

注意到，声波方程中几何衍射因子G(λ，r，L)只是样品尺寸的函数，故上述方程中右边第二项中G_a(λ，r，L)与G_s(λ，r，L)比值为1，即第二项为零，计算方程可写为：

${\mathrm{\α}}_s\left(\mathrm{\ω}\right)=\frac{20}{L}{\log }_{10}[\frac{A_r^a\left(\mathrm{\ω}\right)}{A_r^s\left(\mathrm{\ω}\right)}\·(1-R^2)]$

利用上式中频率与振幅谱比关系，拟合频率-振幅谱比曲线，其斜率为衰减系数。利用衰减系数与品质因子关系：

$\frac{1}{Q}=\frac{{\mathrm{\α}}_s(\mathrm{dB}/\mathrm{\λ})}{8.68\mathrm{\π}}$

我们可以得到测试岩石样品的品质因子。

具体计算如下：

(1)确定参考样品与岩石样品所采集的记录，如图2所示，为本发明实施例参考样品与岩石样品声波测试信号图。

(2)采用傅里叶变换，分别计算参考样品和岩石样品声波信号的频谱。图3为本发明实施例参考样品与岩石样品振幅谱示意图。

(3)分别对参考样品与岩石样品振幅谱取对数，并进行线性拟合，计算振幅谱比值的斜率，确定地层品质因子，参见图4所示，为本发明实施例振幅谱比法计算地层品质因子示意图。

(4)根据不同饱和度点岩石样品测试结果，计算得到品质因子与流体饱和度变化关系曲线，如图5所示，为本发明实施例地层品质因子与饱和度关系曲线图，使用该曲线可对地震衰减属性进行标定和分析。

影响地层品质因子大小的因素有很多，但孔隙流体的存在是最重要且作用最大的因素。在地层沉积和埋深条件变化较小时，地层孔隙流体性质及其赋存状态决定了地层品质因子的大小。一般认为，地震波通过干燥岩石时，地震波能量衰减非常小；当地层含流体后，地震波能量衰减显著。含流体地层引起地震波能量衰减的原因主要是孔隙流体与岩石骨架的相对运动差异造成，当部分饱和流体岩石受到应力波的作用，岩石骨架与孔隙流体都会产生位移变化，岩石骨架因其刚性较强，其应变大小直接与所受应力成正比，但孔隙流体位移则需要一定的能量维持孔隙流体压力场平衡。因此，含流体地层的品质因子通常随着流体饱和度的增加而减小(地震波衰减更加大)，可能是随着饱和度增加，岩石中孔隙流体体积增加，需要驱动流体运动的能量更大，因此品质因子会减小。当岩石完全饱和流体后，地层品质因子会显著增大(地震波衰减变小)，产生这种现象的原因主要是地层完全饱和后，岩石孔隙中充满了流体，岩石骨架与孔隙流体一起运动，孔隙流体的存在增强了岩石刚性，故地层品质会显著增加。

(2)地震衰减品质因子提取方法

地震衰减属性参数的提取方法较多，M.Bath(1974)提出了振幅频谱比率方法，该方法取两个深度或两个时间上子波，进行频谱分析后得到两个频谱比的斜率，这个斜率就是品质因子的函数。Gladwin和Stacy等根据衰减过程地震波脉冲增宽现象，提出了上升时间原理计算地震衰减参数的方法。Biai进一步研究上升时间的特征，并利用上升时间(或脉冲宽度)进行Q值估测的方法，称之为上升时间法，该方法比频率比率方法的条件容易满足，但困难在于寻找最大斜率点的位置和确定斜率都有误差。Kan(1981)等对频谱比率方法进行改进，试图解决散射衰减分离的问题。Quan(1997)提出了一种基于质心频移的方法从VSP资料中估算地震的吸收衰减，并且导出了Q与质心频率的关系。Shengjie(2000年)提出了一种利用VSP资料通过散射迭代逼近方法计算Q值，试图解决散射衰减与本征衰减分离的困难。

这些方法在无噪声情形下，使用解析信号法计算结果最好，频谱模拟法和频谱比法也能获得很好的结果；在较小噪声的情形下，除相位模拟法和振幅衰减法外，其它方法生成的结果一样好；在强噪声情形下，没有一种方法可靠；一般而言，波形建立时间法和脉冲振幅法强烈地依赖于采样率和资料的质量，计算的品质因子精度不高；振幅衰减法计算精度很差。

为了克服地震记录计算Q值精度差的问题，一般采用振幅谱比方法。该方法可适用于含噪声地震记录的能量衰减分析。采用振幅谱比法计算品质因子的具体实现方法如下：

对于叠后地震记录可将其视为自激自收的一维(双程时间)地震波传播结果，其波动方程可表示为：

$\frac{{\∂}^{2}U(x,\mathrm{\ω})}{{\∂}^2{x}^2}=k^{2}U(x,\mathrm{\ω})$

其中U(x，ω)为平面波位移，x为平面波传播距离，ω为角频率，k为波数。

该方程的解析解可表示为：

U(x+Δx，ω)＝U(x，ω)exp(-ikΔx)

其中i为单位虚数。

对于双程反射波的传播距离增量Δx，可用地震波的传播速度v(f)和传播时间增量Δt表示为：

Δx＝v(f)Δt

对于粘滞性介质，地震波的传播将受到地层吸收作用。大量实验和实际地震数据表面，地层吸收作用呈指数衰减特征，由物理因果原理可知，采用波动方程表征这种地层衰减特征需将波数或角频率中引入虚数分量。例如，将波数用虚数表示，可写为：

$k=\frac{\mathrm{\ω}}{v\left(\mathrm{\ω}\right)}(1-\frac{j}{2Q})$

其中，v(ω)为在频率为ω条件下的地震波相速度。Q为介质的品质因子。

为了将距离转换为时间变量t，可定义一个相对参考频率ω₀，在该频率对应的速度为v(ω₀)。将上述复波数代入波动方程中，可得：

$U(t+\mathrm{\Δt},\mathrm{\ω})=U(t,\mathrm{\ω})\exp (-\frac{\mathrm{\ωv}\left({\mathrm{\ω}}_0\right)}{2\mathrm{Qv}\left(\mathrm{\ω}\right)}\mathrm{\Δt})\exp \left(\frac{\mathrm{i\ωv}\left({\mathrm{\ω}}_0\right)}{v\left(\mathrm{\ω}\right)}\mathrm{\Δt}\right)$

其中，v(ω₀)和v(ω)分别为地震波在参考频率ω₀和某一传播频率ω时的传播速度。

上述方程中有两个指数项，第一个指数项其作用主要是衰减地震波的传播能量，第二个指数项是由于地震波传播速度与频率的关系而引起的相位改变项。

根据Kjartansson频散模型，地震波速度频散关系可表示为：

$v\left(\mathrm{\ω}\right)=v\left({\mathrm{\ω}}_0\right){\left|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_{\text{0}}}\right|}^{\mathrm{\γ}}$

其中：

当地震波速度频散非常小的情况下，有v(ω)＝v(ω₀)

此时，考虑了地层吸收作用的波动方程的解可表示为：

$U(t+\mathrm{\Δt},\mathrm{\ω})=U(t,\mathrm{\ω})\exp (-\frac{\mathrm{\πf}}{Q}\mathrm{\Δt})\exp \left(\mathrm{i\ω\Δt}\right)$

假设，地震波在时间增量Δt₁的频率响应为：

$U(t+\mathrm{\Δ}{t}_1,\mathrm{\ω})=U(t,\mathrm{\ω})\exp (-\frac{\mathrm{\πf}}{Q}\mathrm{\Δ}{t}_1)\exp \left(\mathrm{i\ω\Δt}\right)$

地震波在时间增量Δt₂的频率响应为：

$U(t+\mathrm{\Δ}{t}_2,\mathrm{\ω})=U(t,\mathrm{\ω})\exp (-\frac{\mathrm{\πf}}{Q}\mathrm{\Δ}{t}_2)\exp \left(\mathrm{i\ω\Δt}\right)$

对不同时刻频率响应相除，得：

$\frac{U(t+{\mathrm{\Δt}}_1,\mathrm{\ω})}{U(t+{\mathrm{\Δt}}_2,\mathrm{\ω})}=\exp (-\frac{\mathrm{\πf}}{Q}{\mathrm{\Δt}}_1)/\exp (-\frac{\mathrm{\πf}}{Q}{\mathrm{\Δt}}_2)$

两边取对数，可得：

$\ln \left(\frac{U(t+{\mathrm{\Δt}}_1,\mathrm{\ω})}{U(t+{\mathrm{\Δt}}_2,\mathrm{\ω})}\right)=-\frac{\mathrm{\πf}}{Q}\mathrm{\Δt}$

通过计算传播时间Δt前后振幅谱比与频率域的斜率k

$k=\frac{1}{f}\ln \frac{U(t+{\mathrm{\Δt}}_1,\mathrm{\ω})}{U(t+{\mathrm{\Δt}}_2,\mathrm{\ω})}$

以及传播时间Δt

Δt＝t₂-t₁

其中，t1和t2分别是计算窗口(分析目的层顶-底界面)传播时间。

目的层段地层品质因子则可提取得到：

$Q_i=-\frac{\mathrm{\π\Δt}}{k}$

滑动计算窗口可以得到不同地层的地震品质因子参数。对于三维地震数据，逐条侧线进行计算可得地震品质因子数据体。

(3)地震衰减峰值进行储层流体性质预测方法

随着地震勘探对象的日益复杂，对地震储层预测方法的精度要求显得尤为重要。油气储层的非均质粘弹性质及其对地震波吸收特征成为近年来研究的热点。一般而言，地震波在非均质粘弹介质中传播，其能量会发生衰减。地震波在地层中的衰减主要由两方面因素促成，其一是地震波在地层边界引起散射效应而产生的能量耗散；其二是地层本身孔渗及其流体性质差异引起的本征能量衰减。因此，油气储层所固有的本征衰减是人们利用地震衰减属性进行储层流体检测的基础参数。大量国内外实验测试和实践结果表明，地层本征衰减或品质因子与储层性质、孔渗分布、孔隙流体性质及饱和度等参数密切相关。与地震波其它参数(走时，速度等)相比，地震波能量衰减对孔隙流体更加敏感。因此，研究如何利用地震衰减属性参数进行储层流体检测、合理解释地震波振幅与储层关系是一项重要的基础研究工作。

地震波衰减属性近年来逐渐受到人们的重视，该参数常用于地震储层预测和流体检测。因此，寻求有效的地震衰减属性检测储层流体方法显得尤为重要，具有重要实际意义。地震波的衰减与储层骨架性质、孔渗分布、流体性质及其饱和度等相关，通过对衰减属性研究，不仅有助于地震资料的解释，而且还可利用地层品质因子进行储层性质预测、流体检测以及油气藏动态描述。实际应用表明，储层中若孔隙发育，孔隙中饱和流体(油、气、水)后，地震波能量衰减加剧，特别是高频成分衰减严重，含油气地层衰减参数与非油气层差异非常大。目前，利用地震波在油气储层中能量衰减特征进行地震储层流体检测的方法比较多，主要分为地震记录频谱分析技术、时频分析技术以及地震品质因子反演分析技术等。频谱分析技术根据油气储层对地震信号的衰减作用而形成的频谱特征，提取各种反映频谱变化特征的衰减属性参数，如频率衰减梯度、峰值频率、总能量、频谱能量百分比等参数，进行定性油气储层地震检测研究；时频分析技术则通过对地震记录进行数学变换，将地震记录分解为各种频率分量构成的单频或部分频率分量地震记录，利用油气储层吸收特性，对比不同频率分量地震数据，进行油气储层定性预测的方法；地震品质因子反演分析技术采用不同理论模型(振幅谱比值法、上升时间法、各类模拟方法等)由不同时刻地震记录直接反演求取品质因子参数(称之为地震品质因子)，形成地震品质因子剖面或数据体，根据已知井油气储层位置的地震品质因子变化特点，应用地震品质因子剖面(或数据体)进行地震储层预测方法。

上述利用地震衰减属性参数进行储层预测方法存在不足。其一，由于地震记录质量等原因，上述各类方法很难从地震记录中提取到准确的品质因子参数，如频谱分析与时频分解等方法仅依据地震记录频率成分变化(严格意义，频率变化与品质因子是不同的概念，品质因子的变化可引起地震波的高频衰减程度不同，但地震波频率变化并非一定是因品质因子改变而引起)；其二，这些预测方法是一类利用地震衰减属性进行油气储层定性的预测方法，地震检测结果中存在许多不确定因素，地震储层预测结果精度不高，多解性强。

为克服上述地震衰减属性进行储层预测的不足，本发明提出利用储层衰减峰值饱和度这一标志，根据地震品质因子参数的相对值变化特征，寻求地震衰减剖面中储层位置的衰减峰值区域，结合储层品质因子测量结果，合理检测储层孔隙流体的分布范围。本发明利用地震品质因子参数的相对变化值，通过分析衰减曲线变化特征，进行储层流体检测。因此，地震记录质量对本方法的影响较小，且直接使用品质因子相对变化关系检测储层流体可提高储层预测精度。

根据本方法，使用地震品质因子进行储层流体性质检测时，需要结合储层品质因子和流体饱和度关系，寻找地震品质因子剖面中品质因子的极值位置和变化率特征位置，以便确定含流体位置。寻求地震品质因子剖面中品质因子的极值位置的关键方法是对地震品质因子进行微分运算，形成新的微分衰减属性参数，具体实现方法如下：

设地震品质因子剖面中某地层品质因子为Q(j)，则根据定义，地震品质因子求导可写为：

$Q_j^{\′}\≈\frac{Q(j+\mathrm{\Δx})-Q\left(\text{j}\right)}{\mathrm{\Δx}};$ j＝1，2，...N；

其中，j表示第j道位置；Δx为道间距；N为分析时窗中地震记录个数。

由上述方程直接计算品质因子导数会产生较大误差。因为，Q(j+Δx)的泰勒级数展开在应用时必然存在截尾误差与舍入误差的影响，即有：

$Q(j+\mathrm{\Δx})=Q\left(j\right)+Q^{\′}\left(j\right)\mathrm{\Δx}+\frac{1}{2}{\mathrm{\Δx}}^2{Q}^{\′\′}\left(j\right)+\frac{1}{6}{\mathrm{\Δx}}^3{Q}^{\′\′\′}\left(j\right)+\·\·\·$

通常，将一阶倒数项以上省去，产生截尾误差。

克服上述误差影响的方法之一是采用中心均差方程进行计算：

$Q_j^{\′}=\frac{Q_{j+1}-Q_{j-1}}{2\mathrm{\Δx}}-\frac{1}{6}{\mathrm{\Δx}}^3{Q}^{\′\′\′}\left(\mathrm{\ζ}\right),x_{j-1}\≤\mathrm{\ζ}{\≤x}_{j=1,}$ j＝1，2，…N-1

上式中

为高阶截尾项，即误差项。

为了提高计算结果的精度，我们采用三点近似计算方法，计算方程如下：

$Q_j^{\′}=\frac{-3Q_j+4Q_{j+1}-Q_{j+2}}{2\mathrm{\Δx}}+\frac{1}{3}{\mathrm{\Δx}}^2{Q}^{\′\′\′}\left(\mathrm{\ζ}\right)\text{,}$ j＝1，2…，N-2

对于分析时窗中最后一个点的导数计算可采用下面公式：

$Q_N^{\′}=\frac{Q_{j-2}-4Q_{j-1}+{3Q}_j}{2\mathrm{\Δx}}+\frac{1}{3}{\mathrm{\Δx}}^2{Q}^{\′\′\′}\left(\mathrm{\ζ}\right)\text{,}$ j＝3，4…，N

上述方程中Δx²相关项为误差项。

根据该方程，可实现地震品质因子微分运算的程序化计算。

由地层品质因子随饱和度变化关系曲线可知，利用地震品质因子参数检测储层流体性质时，首先需要确定地震品质因子的极值点位置，或微分衰减属性的最小绝对值位置；然后，根据该位置附近微分衰减参数的变化特征，微分衰减参数缓慢变化一侧为低含流体饱和度方向，微分衰减参数变化急剧一侧为接近完全饱和流体的一侧。实际应用时，还可利用已知钻测井确定的地层信息，对微分衰减参数进行标定，进一步提高微分衰减属性参数检测储层流体性质的可靠性。

本发明的效果结合如下实例进行说明：

有如下五个地层构成的地质剖面，见图6所示，为本发明实施例含油储层地质剖面图。第一层为Q＝200(吸收能量较小)的致密岩石构成，平均厚度约为200米；第二层与第四层为Q＝164的泥岩层，第二层厚度约为250米，第四层厚度约为100米；其间夹有一储层，厚度为80米。泥岩层下覆为Q＝176的非储层。储层孔隙中充填了油气，含油饱和度由上而下逐渐增大，并在储层底部完全饱和油，储层段地层品质因子的变化趋势类似于图5所示，衰减峰值饱和度层段(图中红色区域)Q＝80；完全饱油层段品质因子为100，储层段最小含油饱和度为50％。

根据含油储层构造特征，采用二维粘弹波动方程进行了正演计算，并对数值模拟的道集进行叠前处理，处理后叠加地震记录见图7所示，为本发明实施例含油储层地震记录剖面图。该剖面中椭圆标注区域为储层位置，因该储层较厚，储层与围岩存在较大的弹性差异，故其在地震剖面中反射能量与特征明显。采用上述地震品质因子计算方法，求取了该地震剖面的品质因子参数，计算结果见图8所示，为本发明实施例地震品质因子剖面图。

由地震记录求取的地震品质因子是一个相对品质因子变化属性参数，它反映了地层品质因子的相对变化关系。由计算结果可知，储层位置地震品质因子存在一定的异常，但该异常范围比较大，与其上参数分布特征相近，难以确定储层边界的准确位置。常规的地震衰减属性分析方法需要结合类似图6的已知储层位置及地层品质因子等信息，对衰减属性剖面标定，检测衰减属性参数在剖面中异常区域和边界。但因地震品质因子的分辨率很低，使用该参数进行油气检测只能是一种定性分析方法。

为了克服上述地震品质因子储层检测的分辨率低、可靠性差的不足，采用本发明所述微分衰减参数计算方法，求取了微分衰减属性参数，并对计算结果进行归一化处理，计算结果见图9所示，为本发明实施例含油储层地震微分衰减属性参数剖面。

对比地震品质因子剖面与微分衰减参数剖面可知，在微分衰减参数剖面中含油储层衰减峰值位置出现微分衰减参数异常，其特征比较明显，衰减峰值饱和度位置易于分辨(即蓝色与红色相交位置)。与此相比，地震品质因子剖面中，虽然油气储层位置地震品质因子明显降低，但地震品质因子的变化较小，与模型所示储层品质因子特征有一定差别。

上述模型计算结果表明，采用本发明所述微分衰减属性预测方法与传统衰减属性检测方法相比具有显著的效果，本方法可以提高地震储层流体检测的可靠性和精度。

实例为我国西部某油田三维地震衰减属性流体性质检测分析，实现方法如下：

首先，根据研究区储层地层品质因子测试结果，确定各个地层地层品质因子的变化范围和均值，建立研究区各地层单元衰减模型，如图10所示，为本发明实施例各地层单元地层衰减模型示意图。并针对主要目的层，开展储层品质因子与含流体饱和度关系曲线研究。

该地层衰减模型中各层系衰减量反映了地层平均品质因子的大小。总体而言，研究区地层品质因子随着地层深度增加而增大，较浅地层的衰减量变化较大，且地层吸收作用强；而较深地层的衰减量变化相对缓慢，地层吸收作用较之浅层比较弱。地层衰减量的变化特征一方面与岩石致密程度有关，另一方面也与各层系地层所受压实作用有关。

地震波的衰减信息中包括了地层本征衰减信息和其它引起地震波衰减的因素。为了分析地层本征衰减引起地震波衰减的特征，依据已知钻井W1井的地层模型及品质因子测试结果，分别计算了地层无衰减和有衰减情况下地震记录，如图11所示，为本发明实施例地层衰减正演模型记录示意图。在计算衰减地层的地震响应时，分别依据图10中侏罗系地层衰减量及含油砂层测试的地层品质因子，以侏罗系地层为目标进行正演计算。其中含油砂层品质因子是根据相应层段试验测量结果确定的，图中各模型记录之间所插入曲线为W1井声波速度曲线。所使用的地层品质因子分别为208和105。

正演模型记录中，为便于对比，将反射振幅强度进行了归一化处理。由该图可知，随着地层品质因子减小，模型记录中的反射振幅减弱。与无衰减记录相比，当地层品质因子Q＝208时，储层位置反射振幅减小了75％左右。而当Q＝105时，反射振幅约减小95％左右。

利用地层品质因子进行储层横向预测时，图11中Q＝208的正演记录可视为衰减因子的背景变化趋势，而将Q＝105的正演记录作为含油层衰减异常区域的一种标志。两者之间反射能量相差约60％。

衰减数据中的品质因子的变化范围为50-550左右。将地震品质因子参数与测井曲线换算的品质因子曲线相比，地震品质因子剖面中的品质因子变化范围大，且普遍高于测井计算的品质因子曲线(这可能说明地震品质因子参数中包含了非本征衰减信息)。由于地震波的频带及波长等因素，使得地震品质因子参数剖面无法反映不同地层的界面特征。使用Q＝208作为衰减背景趋势分析，可突出含油砂层的位置。如图12所示，为本发明实施例含油砂层在地震品质因子剖面中特征示意图。

地震品质因子剖面中含油砂层的基本轮廓比较明显。含油砂层的底界位置如图中箭头所指位置。显然衰减量最小的位置并不在含油砂层所处位置处(略高于含油砂层反射位置)。因此，只能定性地根据正演模型中衰减量的背景值与含油砂层关系对地震品质因子剖面进行储层流体检测分析。

为了提高地震储层检测的可靠性和精度，结合该储层峰值衰减饱和度曲线(图5所示)，对地震品质因子数据体进行微分计算，求取微分衰减属性参数记录。如图13所示，为本发明实施例微分衰减属性参数剖面图。

根据衰减峰值饱和度特征，地震流体检测可利用峰值饱和度衰减前、后(特别是峰值饱和度以后衰减因子急剧减小)品质因子的快速变化特征，分析含流体地层的空间位置。

由地震品质因子数据计算的微分衰减属性参数剖面中(图13所示)，已知含油砂层位置正是处于品质因子变化最大的位置(如图13中箭头所示)。这与上述利用峰值饱和度概念预测储层流体性质的思路完全一致，微分衰减参数在横向上分布与含油砂层的构造形态比较接近。通过对微分衰减属性与测井计算品质因子标定，微分衰减属性异常区域与已知井油层位置一致。

根据上述分析，将微分衰减参数最低点和最大变化率区域作为预测含油砂层的有利位置，对三维数据体进行了储层检测分析工作。为便于对比，分别绘制了地震品质因子预测储层平面图和采用微分衰减属性参数检测的含油储层平面图。如图14所示，为本发明实施例衰减因子预测有利含油砂层分布区域示意图，其中，图14中(a)为微分衰减属性参数的；图14中(b)为地震品质因子参数的。

研究区内侏罗系三工河组底部砂层为主要目的层，研究区内共有钻井三口，其中，CHE71和CHE26在三工河组见工业油流，CHE29井为干井。对比微分衰减属性参数与地震品质因子的预测平面图可知，地震品质因子较小的区域(图14中(b)中椭圆所示区域)与微分衰减属性参数异常区域密切相关(如图14中(a)箭头所指位置)，而在图14中(a)中其它区域，微分衰减属性参数相对比较高，但微分衰减属性参数分布区域比较集中，其边界易识别。对比研究区内已知钻井油气显示结果，说明利用微分衰减属性参数预测的含油砂层平面分布是可靠的。

尽管利用地震品质因子预测结果和微分衰减属性预测的有利区域在平面上比较相似，但微分衰减属性参数预测油气有利区域的边界清楚、有利区域分布合理。而地震品质因子预测结果较为散乱。根据已知钻井的油气显示结果可知，微分衰减属性参数进行流体检测具有较高的可靠性和有效性。

本发明首次提出根据利用储层流体饱和度与品质因子关系进行地震储层预测方法，详细阐述了含流体储层品质因子测试方法、地震数据提取地震品质因子技术以及利用地震衰减峰值饱和度特征进行储层预测的方法。其中含流体岩石衰减峰值饱和度及微分衰减属性概念是本发明首创。模型数据与实际地震数据应用表明，本发明较之传统地震衰减预测方法具有显著的有效性和可靠性。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分步骤是可以通过程序来指令相关硬件来完成，所述的程序可以存储于一计算机可读取存储介质中，该程序在执行时，包括上述全部或部分步骤，所述的存储介质，如：ROM/RAM、磁盘、光盘等。

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种利用地震衰减属性进行储层流体性质预测的方法，其特征在于，所述方法包括：

确定储层饱和度与地层品质因子关系，所述储层饱和度与地层品质因子关系包括储层饱和度与地层品质因子变化曲线、衰减峰值出现的饱和度位置；

对地震数据进行储层反射界面构造解释，确定储层在地震剖面中的反射时间位置；

确定储层品质因子的大小；

沿反射界面对地层品质因子进行平滑与处理，对地震品质因子数据采用微分计算方法确定微分衰减属性大小；

利用储层饱和度与地层品质因子变化曲线分析，确定储层流体性质分布特征；

其中，所述确定储层饱和度与地层品质因子关系，具体包括：

（1）确定参考样品与岩石样品所采集的记录；

（2）采用傅里叶变换，分别计算参考样品和岩石样品声波信号的频谱；

（3）分别对参考样品与岩石样品振幅谱取对数，并进行线性拟合，计算振幅谱比值的斜率，确定地层品质因子；

（4）根据不同饱和度点岩石样品测试结果，计算得到品质因子与流体饱和度变化关系曲线；

其中，所述确定储层品质因子的大小，具体包括：

对于叠后地震记录能将其视为自激自收的一维地震波传播结果，其波动方程表示为：

$\frac{{\∂}^{2}U(x,\mathrm{\ω})}{{\∂}^2{x}^2}=k^{2}U(x,\mathrm{\ω})$

其中U(x,ω)为平面波位移，x为平面波传播距离，ω为角频率，k为波数；

该方程的解析解表示为：

U(x+Δx,ω)=U(x,ω)exp(-ikΔx)

其中i为单位虚数；

对于双程反射波的传播距离增量⊿x，用地震波的传播速度v(f)和传播时间增量⊿t表示为：

Δx=v(f)Δt

对于粘滞性介质，地震波的传播将受到地层吸收作用，大量实验和实际地震数据表明，地层吸收作用呈指数衰减特征，由物理因果原理可知，采用波动方程表征这种地层衰减特征需将波数或角频率中引入虚数分量，将波数用虚数表示，写为：

$k=\frac{\mathrm{\ω}}{v\left(\mathrm{\ω}\right)}(1-\frac{j}{2Q})$

其中，v(ω)为在频率为ω条件下的地震波相速度，Q为介质的品质因子；

为了将距离转换为时间变量t，定义一个相对参考频率ω₀，在该频率对应的速度为v(ω₀)，将复波数代入波动方程中，可得：

$U(t+\mathrm{\Δt},\mathrm{\ω})=U(t,\mathrm{\ω})\exp (-\frac{\mathrm{\ωv}\left({\mathrm{\ω}}_0\right)}{2\mathrm{Qv}\left(\mathrm{\ω}\right)}\mathrm{\Δt})\exp \left(\frac{\mathrm{i\ωv}\left({\mathrm{\ω}}_0\right)}{v\left(\mathrm{\ω}\right)}\mathrm{\Δt}\right)$

其中，v(ω₀)和v(ω)分别为地震波在参考频率ω₀和某一传播频率ω时的传播速度；

上述方程中有两个指数项，第一个指数项其作用主要是衰减地震波的传播能量，第二个指数项是由于地震波传播速度与频率的关系而引起的相位改变项；

根据Kjartansson频散模型，地震波速度频散关系表示为：

$v\left(\mathrm{\ω}\right)=v\left({\mathrm{\ω}}_0\right){\left|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_0}\right|}^{\mathrm{\γ}}$

其中： $\mathrm{\γ}\≈\frac{1}{\mathrm{\πQ}}$

当地震波速度频散非常小的情况下，有v(ω)=v(ω₀)

此时，考虑了地层吸收作用的波动方程的解表示为：

$U(t+\mathrm{\Δt},\mathrm{\ω})=U(t,\mathrm{\ω})\exp (-\frac{\mathrm{\πf}}{Q}\mathrm{\Δt})\exp \left(\mathrm{i\ω\Δt}\right)$

假设，地震波在时间增量Δt₁的频率响应为：

$U(t+\mathrm{\Δ}{t}_1,\mathrm{\ω})=U(t,\mathrm{\ω})\exp (-\frac{\mathrm{\πf}}{Q}\mathrm{\Δ}{t}_1)\exp \left(\mathrm{i\ω\Δt}\right)$

地震波在时间增量Δt₂的频率响应为：

$U(t+\mathrm{\Δ}{t}_2,\mathrm{\ω})=U(t,\mathrm{\ω})\exp (-\frac{\mathrm{\πf}}{Q}\mathrm{\Δ}{t}_2)\exp \left(\mathrm{i\ω\Δt}\right)$

对不同时刻频率响应相除，得：

$\frac{U(t+\mathrm{\Δ}{t}_1,\mathrm{\ω})}{U(t+\mathrm{\Δ}{t}_2,\mathrm{\ω})}=\exp (-\frac{\mathrm{\πf}}{Q}\mathrm{\Δ}{t}_1)/\exp (-\frac{\mathrm{\πf}}{Q}\mathrm{\Δ}{t}_2)$

两边取对数，可得：

$\ln \left(\frac{U(t+\mathrm{\Δ}{t}_1,\mathrm{\ω})}{U(t+\mathrm{\Δ}{t}_2,\mathrm{\ω})}\right)=-\frac{\mathrm{\πf}}{Q}\mathrm{\Δt}$

通过计算传播时间⊿t前后振幅谱比与频率域的斜率k；

$k=\frac{1}{f}\ln \frac{U(t+\mathrm{\Δ}{t}_1,\mathrm{\ω})}{U(t+\mathrm{\Δ}{t}_2,\mathrm{\ω})}$

以及传播时间⊿t

Δt=t₂-t₁

其中，t₁和t₂分别是计算窗口分析目的层顶界面及层底界面的传播时间；

目的层段地层品质因子则能提取得到：

$Q_i=-\frac{\mathrm{\π\Δt}}{k}$

滑动计算窗口能得到不同地层的地震品质因子参数；

其中，所述利用储层饱和度与地层品质因子变化曲线分析，确定储层流体性质分布特征，具体包括：

确定地震品质因子的极值点位置，或微分衰减属性的最小绝对值位置；

根据所述极值点位置或所述最小绝对值位置附近微分衰减参数的变化特征，微分衰减参数缓慢变化一侧为低含流体饱和度方向，微分衰减参数变化急剧一侧为接近完全饱和流体的一侧；

其中，所述确定地震品质因子的极值点位置，或微分衰减属性的最小绝对值位置，具体为：

设地震品质因子剖面中某地层品质因子为Q(j),则根据定义，地震品质因子求导写为：

$Q_j^{\′}\≈\frac{Q(j+\mathrm{\Δx})-Q\left(j\right)}{\mathrm{\Δx}};j=\mathrm{1,2},...N;$

其中，j表示第j道位置；Δx为道间距;N为分析时窗中地震记录个数；

由上述方程直接计算品质因子导数会产生较大误差；因为，Q(j+Δx)的泰勒级数展开在应用时必然存在截尾误差与舍入误差的影响，即有：

$Q(j+\mathrm{\Δx})=Q\left(j\right)+Q^{\′}\left(j\right)\mathrm{\Δx}+\frac{1}{2}{\mathrm{\Δx}}^2{Q}^{\′\′}\left(j\right)+\frac{1}{6}{\mathrm{\Δ}{x}^{3}Q}^{\′\′\′}\left(j\right)+...$

通常，将一阶倒数项以上省去，产生截尾误差；

克服上述误差影响的方法之一是采用中心均差方程进行计算：

$Q_j^{\′}=\frac{Q_{j+1}-Q_{j-1}}{2\mathrm{\Δx}}-\frac{1}{6}{\mathrm{\Δx}}^3{Q}^{\′\′\′}\left(\mathrm{\ζ}\right),x_{j-1}\≤\mathrm{\ζ}\≤x_{j=1},j=\mathrm{1,2},...N-1$

上式中

为高阶截尾项，即误差项；

为了提高计算结果的精度，我们采用三点近似计算方法，计算方程如下：

$Q_j^{\′}=\frac{-3Q_j+{4Q}_{j+1}-Q_{j+2}}{2\mathrm{\Δx}}+\frac{1}{3}{\mathrm{\Δx}}^2{Q}^{\′\′\′}\left(\mathrm{\ζ}\right),j=\mathrm{1,2}...,N-2$

对于分析时窗中最后一个点的导数计算采用下面公式：

$Q_N^{\′}=\frac{Q_{j-2}-{4Q}_{j-1}+{3Q}_j}{2\mathrm{\Δx}}+\frac{1}{3}{\mathrm{\Δx}}^2{Q}^{\′\′\′}\left(\mathrm{\ζ}\right),j=\mathrm{3,4}...,N$

上述方程中Δx²相关项为误差项；

根据该方程，可实现地震品质因子微分运算的程序化计算。

Images