基于不同階次泰勒級數展開的含非本征衰減頻域振幅比平均的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）?？梢姡孩僭诘蛥⒖碱l段，只有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）?？梢姡弘S著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值估計結果?？梢姡簾o論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估計值?？梢钥闯觯孩佼敃r差較?。?00 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 數據，同時也適用于疊后地震記錄，而且無需事先進行非本征衰減補償，具有較強的適用性和實用性。