基于不同阶次泰勒级数展开的含非本征衰减频域振幅比平均的Q值估计方法

张瑾，王彦国，兰慧田，张国书，潘耶莉

（1.江西省防震减灾与工程地质灾害探测工程研究中心，江西南昌 330013；2.东华理工大学地球物理与测控技术学院，江西南昌 330013；3.东华理工大学核资源与环境国家重点实验室，江西南昌 330013；4.大庆油田勘探开发研究院，黑龙江大庆 163712；5.东华理工大学核科学与工程学院，江西南昌 330013）

0 引言

现今油气勘探目标日益复杂，对地震勘探提出了更高要求。实际地下介质是黏弹性介质，地震波在其中传播将发生两类衰减：一类是与地层吸收性质相关的衰减，称为本征衰减；另一类是与地震波传播性质有关的衰减，称为非本征衰减，如球面扩散、透射损失等[1]。

品质因子Q值作为描述本征衰减的关键参数，在提高岩性预测精度及地震资料纵向分辨率方面具有重要作用。一方面，Q值与地层岩性、饱和度及渗透率等参数密切相关，可用于储层预测和油气识别[2-5]；另一方面，Q是能量补偿反Q滤波方法的必要参数，可间接提高地震资料成像精度[6-8]。

自20世纪70年代起，Q值估计方法备受人们关注。其中，对数谱比（LSR）法是常用的传统算法，具有广泛的适用性，可估计透射波数据（如VSP 资料[9-11]）及地面反射波数据（如叠后地震资料[12-13]和叠前CMP数据[14]）的Q值。但LSR 法受频段、子波干涉和噪声等因素的影响较大[15-17]，限制了其在实际地震资料中的应用。

为了克服这些因素的影响，人们一直在寻求优化算法。Wang 等[17]基于积分运算弱化噪声的思想，利用对数谱面积差（LSAD）线性拟合对数谱斜率求取Q值，提高了算法的抗噪性，但Q值估计结果受非本征衰减的影响，仅适用于VSP数据。为消除非本征衰减的影响，An 等[18]、Cheng 等[19]利用高、低参考频段对数谱面积差（LSAD_LH）法及分段频率对数谱面积双差值（LSADD）法求取Q值，进一步改善了算法的实用性。许李囡等[20]在S 变换域利用变分法估计Q值，受震源子波类型及时窗宽度的影响较小。杨登峰等[21]基于LSR 法加入高斯函数权重算子求取斜率，降低了Q值估计算法对频段选择的依赖性。Zhang 等[22-23]基于振幅谱积分差和振幅谱积分比值稳定地估计Q值。

目前，Q值估计方法发展迅速，但易受频段选择及非本征衰减等因素影响，Q值估计精度及稳定性不高。为此，本文提出基于不同阶次泰勒级数展开的含非本征衰减频域振幅比平均的Q值估计方法。首先，利用高、低双参考频段消除非本征衰减因素对Q值估计的影响，提出考虑非本征衰减的1～4阶频域振幅比单频点（FARS）Q值估计算法。在此基础上，采用频点平均模式构建多频点平均（FARA）Q值估计方法。最后通过模型试验和实例应用测试方法的效果和精度。

1 方法原理

1.1 FARS 法的基本原理

基于常Q值的Futterman 体波衰减模型[24]，在考虑非本征衰减影响时，地震波在黏弹性介质中的传播规律可用不同时刻的振幅变化表述[25]

式中：P（t）为t时刻的非本征衰减项[17-18]；Δt为时间延迟；A（0,f）为震源子波的振幅，f为频率。将式（1）与式（2）相比，得

为了消除非本征衰减对Q值估计的影响，定义一个相对参考频段内振幅比的连乘，即

求解以上四式，分别得到x在1～4 阶泰勒级数展开的正解

其中

而

则不同阶次i（i=1,2,3,4）泰勒级数展开的Q值为

式（14）为不同阶次FARS法估计Q值基本公式。为削弱参考频段选取的影响，可选取低参考频段（如10～20 Hz）和高频段（如（2fm－10）～2fm，fm为地震子波的主频）的组合模式求取Q值。

1.2 FARA 法的基本原理

由于本文研究不随频率变化的Q估计值，因此可在振幅谱主值频段[fmin，fmax]采用FARA 法提高Q值估计的稳定性

2 模型试验

2.1 参考频段对Q 估计值的影响

图1为两层介质模型在200、300 ms时刻的含非本征衰减地震记录。鉴于本文算法由子波振幅谱获取Q估计值，故展示图1中4个衰减地震子波的振幅谱（图2）。可以看出，随着旅行时间增加及Q值减小，主频向低频方向移动，不同时刻衰减子波的振幅关于主频呈近似对称分布。基于这一思想，本文将在主值频段[22]估计Q值。

图1 两层介质模型在200、300 ms时刻的含非本征衰减地震记录

图2 图1 中4 个衰减地震子波的振幅谱

选取不同参考频段，利用FARS 法估计Q值（图3）。由图可见：①在低参考频段，随着频率增加，Q估计值先迅速增大，后缓慢下降，在20 Hz之后Q估计值逐渐趋于真值；在参考频段出现突变（图3a 的14 Hz 处）。②在高参考频段，在超低频率处Q估计值随着频率增加迅速增大，随后趋于平缓；在参考频段也存在突变（图3b 的88 Hz 及70 Hz 处）。因此，图3a、图3b的Q估计值在大部分频段分布在真实Q值两侧，其中Q1估计值最明显。③将低、高参考频段组合，提高了Q值估计精度（图3c），Q估计值在大部分频段更接近真值。

图3 真实Q 值为160（左）、40（右）的Q 值估计结果

2.2 Q 估计值影响因素分析

为了探究参考频段、时差、时窗、噪声等因素对FARA 法的影响，进行了相应的测试分析。此外，LSADD 法[19]可直接由含非本征衰减的地震记录估计Q值，与FARA 方法的理论基础较相似。因此，对比、分析了FARA 法与LSADD 法的效果。

2.2.1 参考频段的影响

FARA 法在估计Q值时需给定参考频段，基于图1，分别对不同衰减子波（Q=160、120、80、40）由FARA 法估计Q值。鉴于FARA 法在高、低参考频段附近的Q估计值不稳定，故在利用FARA 法计算Q值时，选择主值频段为计算频段。如：在低参考频段10～(10+K)（模型中频率采样率为1 Hz，为简便表达，下文直接将与参考频点K有关的频率值表示为K），计算频段为(10+K)～2fm[22]；在高参考频段(2fm-K)～2fm，计算频段为10～（2fm-K）。这里选取低参考频段10～20 Hz 及其计算频段20～2fm、高参考频段（2fm-10）～2fm及其计算频段10～(2fm-10)估计高、低参考频段Q值，并求高、低频段Q估计值的平均值（表1）。可见：①在低参考频段，只有1 阶FARA 法的Q估计值明显低于理论值，2～4 阶FARA 法的Q估计值与理论值相差不大。②在高参考频段，只有1 阶FARA 法的Q估计值明显偏大，2～4 阶FARA 法的Q估计值基本上都略大于理论值，相对误差（估计值和真值之差与真值的比值）略高于低参考频段。③高、低频段Q估计值的平均值表明，1 阶FARA 法有效提高了Q值估计精度；无论是低参考频段、还是高参考频段，由于2～4阶FARA 法的Q估计值均略大于理论值，因此高、低频段Q估计值的平均值不能明显提高Q值估计精度。不同阶次FARA 法获得的高、低频段Q估计值的平均值最大相对误差约为3%（Q= 40，1阶FARA 法），说明由FARA 法求取Q值的精度较高。需要说明的是，下文均采用高、低频段Q估计值的平均值作为FARA 法的Q估计值。

表1 不同参考频段的FARA 法Q 估计值

为了探究参考频带宽度（即参考频点个数K）对FARA 法的影响，基于图1 采用FARA 法获取Q估计值，并将高、低频的K均从1 增至40，由此得到Q估计值随K的变化曲线（图4）。可见：随着K增加，所有阶次的Q估计值基本呈现先逐渐下降再缓慢上升的趋势；虽然不同K的Q估计值都大于理论值，但最大偏差均未超过3%，可见选取不同的K基本不影响Q值估计精度，但适当增加K，可提高Q值估计的稳定性。

图4 Q 估计值随K 的变化曲线

试验结果表明：对于FARA 法来说，选择低参考频段计算Q值的精度高于高参考频段，采用高、低频段Q估计值的平均值可以提高Q值的估计精度，对于低阶次Q值更是如此；参考频点个数K对Q值估计结果影响不大，但适当增加K可以在一定程度上提高算法的稳定性。

2.2.2 旅行时差的影响

由于FARA 法是一种基于泰勒级数展开的近似算法，其误差将随着旅行时与Q的比值增大而增大。为了测试FARA 法对旅行时差的敏感程度，基于图1，固定参考子波的旅行时为200 ms，选择衰减子波的旅行时为250、300、400 ms，子波截取的时窗宽度为50 ms，高、低频的K均为40。图5 为不同旅行时差的LSADD 法和FARA 法的Q值估计结果。可见：无论Q值大小，LSADD 法的Q估计值均大于理论值，且随着旅行时差的增加而增大；不同阶次FARA 法的Q估计值基本也随旅行时差的增加而增大，但增大幅度整体小于LSADD 法；无论旅行时差大小，各阶次FARA法的Q值估计精度均高于LSADD 法。

图5 不同旅行时差的LSADD 法和FARA 法的Q 值估计结果

试验结果表明：FARA 法由于采用高、低频段Q估计值的平均值，即使1阶Q估计值的精度也较高，但误差随着旅行时差的增大而略微增大，且任意阶次的Q值估计精度均高于LSADD 法。

2.2.3 时窗宽度的影响

为测试截取子波时窗宽度（子波不完整性）对FARA 法的影响，基于图1，得到Q估计值随时窗宽度的变化曲线（图6）。可见，LSADD 法、各阶次FARA法的Q估计值随时窗宽度的变化曲线基本一致，即随着时窗宽度增加先迅速降低、再逐渐增大、再缓慢下降并逐渐趋于真值。显然，当旅行时差较小时，Q估计值波动较大，明显偏离真值。当截取子波宽度大于25 ms（完整子波延续宽度约为50 ms）时，FARA 法的Q估计值的波动小于LSADD 法，估计精度更高。

图6 Q 估计值随时窗宽度的变化曲线

试验结果表明，FARA 法和LSADD 法在一定程度上均受子波不完整性的影响，但当时窗宽度大于子波宽度一半时，两种方法均能获得稳定、准确的Q估计值，但FARA 法的Q估计值受时窗宽度影响更小。

2.2.4 噪声的影响

为测试FARA 法的抗噪性，将图1 数据添加随机噪声（图7），图8 为加噪前、后衰减地震子波的振幅谱。可见，加噪前、后振幅谱的主要能量均集中于2倍主频之内；且在主值频段（10～2fm），加噪前、后子波的整体一致性较好，说明选取主值频段估计Q值更合理。

图7 图1 数据添加噪声

图8 加噪前、后衰减地震子波的振幅谱

图9～图12 为不同Q值、噪声的LSADD 法和FARA 法的Q估计值概率分布图。由图可见：①当理论Q值较大且噪声较弱时（图9），不同阶次的FARA法的Q估计值几乎完全吻合，较LSADD 法的Q估计值更接近真值（LSADD 法、FARA 法的Q估计值相对误差分别为17.7%、小于5.5%）。②当理论Q值较小且噪声较弱时（图10），LSADD 法与FARA 法均能获得准确的Q估计值，但后者的精度更高。③当噪声较强时（图11、图12），无论理论Q值大小，两种算法的Q估计值波动均增大（LSADD 法、FARA 法的Q估计值相对误差分别为32%、21%），且在理论Q值较大时（图11）得到的Q估计值明显偏小，在理论Q值较小时（图12）得到的Q估计值精度较高，但FARA 法的Q值估计精度更高。④不同阶次FARA 法的Q估计值离散性更小，稳定性更强。

图9 LSADD 法和FARA 法的Q 估计值概率分布图（理论Q 值为160，信噪比为10 dB）

图10 LSADD 和FARA 法的Q 估计值概率分布图（理论Q 值为 40，信噪比为10dB）

图11 LSADD 和FARA 法的Q 估计值概率分布图（理论Q 值为160，信噪比为5 dB）

图12 LSADD 和FARA 法的Q 估计值概率分布图（理论Q 值为40，信噪比为5 dB）

试验结果表明，FARA 法的效果在一定程度上受噪声影响，Q值估计误差随着理论Q值的增大、噪声强度的增加而增大。由于FARA 法采用高、低频段Q估计值的平均值作为FARA 法的Q估计值，因此FARA 法受噪声的影响程度低于LSADD 法，且稳定性更高。

2.2.5 理论Q值与时差的综合影响

由于FARA 法是一种估计Q值的近似算法，其精度受Δt/Q影响，Δt/Q越大，泰勒级数展开的误差越大，因此有必要研究理论Q值、不同时差对FARA 法的综合影响。表2为不同理论Q值、不同时差的FARA法Q估计值。可以看出：①当时差较小（100 ms）、理论Q值较大时，FARA 法的Q估计值相对误差较小，不同阶次FARA 法的Q估计值无明显差别；2 阶FARA 法Q估计值的精度略高。②当时差较大、理论Q值较小时，1 阶、3 阶、4阶FARA 法的Q估计值精度仍较高，Q估计值受Δt/Q影响有限。③整体而言，3 阶、4 阶FARA 法的Q估计值精度更高，主要是由于泰勒级数展开的余项误差更小所致。④对比表1 和表2 可见，无论理论Q值大小，低参考频段的Q估计值均小于高参考频段，高、低频段Q估计值的平均值可弱化泰勒级数展开误差的影响，1 阶FARA 的Q估计值更是如此。试验结果表明，即使时差较大、Q值较小，导致泰勒级数截断误差较大，但也并不会引起低阶次FARA 法的Q值估计失效，这主要是由于低、高参考频段的泰勒级数展开误差恰好抵消所致，进一步证明本文方法能较好地估计Q值。

表2 不同理论Q 值、不同时差的FARA 法Q 估计值

3 实例应用

为了测试新算法的实际应用效果，选取中国南海A 区海洋叠前地震数据进行测试。图13 为中国南海A 区CMP 道集、第1 道振幅谱及第1 段截取的不同道地震子波的振幅谱。由图可见：由于不同主频子波的振幅相互叠加，CMP 道集（图13a）的第1道振幅谱（图13b）曲线波动较大，但振幅谱能量主要集中分布在10～47 Hz 范围。根据子波分布特征，在旅行时分别为532、744、936、1188、1356、1800 及1976 ms处，将整个地震记录道分割为7 段，然后利用LSADD 法和本文提出的FARA 法求取Q值，再求取Q估计值的平均值作为该段的平均等效Q值[17,19]。在计算过程中选用的时窗宽度为48 ms，主值频段仍为10～2fm，FARA 法参考频段的频点数K为40。为了了解处理流程，以第1 段（起点为图13a 红色箭头处）为例详述实际资料处理过程。由第1 段截取的不同道地震子波的振幅谱（图13c）可见：①第1 道子波主频约为29 Hz，其余各道子波主频均降低，由于旅行时差较小，因此主频降低不明显。②随着旅行时间增加，整体上振幅最大值不断减小，所有子波的振幅谱分布于10 Hz 至2 倍主频之间，与图2相似。③第1道～第5道子波主频分别为28.9、26.8、26.4、26.4 及28.6 Hz，说明由于激发、接收条件等因素的影响，导致部分接收子波的主频并没有随着传播时间增加而降低，若采用的参考子波振幅（或主频）小于衰减子波，则Q估计值会出现负值[26-27]。因此在计算平均等效Q值时，需要剔除负值。

图13 中国南海A 区CMP 道集（a）、第1 道振幅谱（b）及第1 段截取的不同道地震子波的振幅谱（c）

图14 为由不同方法估计的平均等效Q值。由图可见：①整体上，不同阶次的FARA 法Q估计值的一致性较好，可相互印证，从而提高Q值估计精度；②FARA 法得到的等效Q值整体上大于LSADD 法，这与由图9～图12 所得认识一致，表明实际地震记录中的噪声导致LSADD 法Q估计值整体偏小。

图14 由不同方法估计的平均等效Q值

4 结论

准确估计品质因子Q值对提高岩性预测及地震成像精度均具有重要作用。本文从含非本征衰减的地震波振幅表达式出发，利用参考频段内振幅比的连乘消除非本征衰减对Q估计值的影响，然后在参考频点处再对振幅衰减因子进行1～4阶泰勒级数展开，由此推导不受非本征衰减影响的振幅比单频点Q值计算公式。此外，利用高、低频段Q估计值的平均值削弱参考频段的影响，并提高算法的稳定性。由模型试验及实例应用得到以下认识：

（1）本文提出的FARA 法可在主值频段内稳定估计Q值，与LSADD 法相比，具有受时差及时窗等因素影响小和抗干扰能力强的优点；

（2）选取主值频段估计Q值，可有效降低计算频段对Q估计值的影响，从而提高小Q值估计的稳定性；

（3）FARA 法包含1～4 阶的Q值估计结果，这些结果互为印证，进一步提高了方法的实用性；

（4）实例应用中，不同阶次FARA 法获得的等效Q值一致性较好，整体上大于LSADD 法的Q估计值，因此由FARA 法得到的Q值更可靠；

（5）FARA 法不仅适用于叠前CMP 道集、VSP 数据，同时也适用于叠后地震记录，而且无需事先进行非本征衰减补偿，具有较强的适用性和实用性。