利用泰勒級數展開的振幅譜積分差值的Q值估計方法

張 瑾 張國書 王彥國*③ 李紅星

(①東華理工大學地球物理與測控技術學院，江西南昌 330013； ②東華理工大學核科學與技術學院，江西南昌 330013； ③東華理工大學核資源與環境國家重點實驗室，江西南昌 330013)

0 引言

地震波在黏彈性介質中傳播時會產生速度頻散和能量衰減[1-2]，而品質因子Q作為測量速度頻散和能量衰減的一個重要參數[3-4]，在地震數據處理和解釋中起重要作用。一方面，品質因子Q是反Q濾波必不可少的參數，可提高地震成像的縱向分辨率[5-6]。另一方面，品質因子Q與巖石類型、流體飽和度和滲透率等有關，常被應用于儲層預測和油氣識別[7-8]。

近年來，由于品質因子Q值的作用日趨顯著，Q值估計方法也得到了快速發展。最早的估計方法是時域估計方法，如上升時間[9]、脈沖振幅法[10]。由于該類算法僅適用于無噪聲的地震記錄，未能得到廣泛的應用。隨著快速傅里葉變換技術的發展，人們相繼提出了一系列的頻域方法，如對數譜比法(Logical Spectrum Ratio，LSR)[11]、中心頻移法(Centroid Frequency Shift，CFS)[12]、小波域包絡峰值瞬時頻率法(Wavelet Envelope Peak Instanta-neous Frequency，WEPIF)[13]。但這些方法的Q值估計準確程度通常受噪聲、時窗、帶寬、子波干涉等因素影響[14-16]。因此，針對Q估計的魯棒性和準確性，人們提出了一系列穩鍵的頻域方法。

趙寧等[17]采用一階、二階泰勒級數展開理論計算質心頻率，進一步提高了Q值估計的精度和抗噪性。Wang等[18]提出對數譜面積差法(Logarithmic Spectral Area Difference，LSAD)，提高了Q估計的穩定性。王靜等[19]基于頻譜擬合法實現VSP地震數據的準確Q值估計。劉國昌等[20]在譜比法的理論基礎上，利用局部斜率相同的反射波形進行Q值估計，得到了較好的Q值估計結果。蘇勤等[21]基于地表一致性原理，通過引入表層相對衰減系數，推導出表層相對品質因子Q值的計算公式，并將相應方法應用于實際數據Q值估計，取得了較好的應用效果。

本文對描述大地濾波作用的地震波振幅衰減因子進行二階泰勒級數展開，從不同時刻地震子波振幅譜的積分差值與Q值關系出發，推導出基于泰勒展開振幅譜積分差值(Amplitude Spectral Integral Difference，ASID)的Q值估計方法； 并通過頻帶寬度、時窗寬度及噪聲干擾試驗，測試了本文方法的有效性，同時與LSR、CFS和LSAD三種Q值估計方法進行對比分析以說明ASID法的優越性。

1 方法原理

忽略地震波在地下傳播過程中的球面擴散和透射損失作用，地層的“大地濾波作用”可用頻域一維波動方程的解析解描述[3]

U(r+Δr,f)=U(r,f)×

(1)

式中： i為虛數單位；f為頻率；r為傳播距離； Δr為傳播距離增量；U(r，f)和U(r+Δr，f)對應傳播距離分別為r和r+Δr的波場；Q是傳播距離r與r+Δr之間的地層品質因子；v(r，f)是相速度，可由Kjartansson相速度模型[4]表示為

(2)

(3)

式中：f0為地震波主頻；v(r,f0)是頻率f0的地震波在r處的相速度；γ是與Q有關的頻率指數項。另外，若主頻為f0的地震波在Δr之內的相速度v(r,f0)保持不變，則Δr可由v(r,f0)表示為

Δr=v(r,f0)Δt

(4)

式中Δt是地震波通過Δr的旅行時間。

將式(2)～式(4)代入式(1)，可得[22-23]

(5)

當頻率在f0附近時，有|f/f0|-(1/πQ)?1[18,24]，因此式(5)可簡化為

U(t+Δt,f)=U(t,f)exp(-2πifΔt)×

(6)

從該式發現Q只影響實數部分，因此不同時刻地震子波振幅譜的關系[16,25]可表示為

(7)

式中A(t+Δt,f)和A(t,f)分別是t+Δt和t時刻的振幅譜。定義旅行時t處振幅譜積分為

(8)

式中F為頻率積分區間，可預先給定，且該區間通常選取包含主頻的頻段。同理，t+Δt時刻的振幅譜積分為

(9)

基于二階泰勒級數展開，式(9)指數項可表示為

(10)

式中Rn(πΔtf/Q)是泰勒公式的余項，與旅行時t、頻率f和品質因子Q有關。當πΔtf/Q<1時，可認為該余項接近于0。因此，估算Q時選取的旅行時差Δt不宜過大。當忽略余項時，式(9)可改寫為

(11)

兩個接收點之間的振幅譜積分差值為

ΔG=G(t)-G(t+Δt)

(12)

上式可簡化為

(13)

其中

雖然式(13)存在兩個解，但通過實驗發現，其中一個解是一個非常小的數，即非Q值的真解，而Q值的有效解為

(14)

該式為本文ASID算法Q值估計的基本公式，即是利用同層位不同旅行時間的地震子波振幅譜之差實現Q值估計，一般適用于疊前多道地震數據Q值估計。

圖1為多層介質疊前地震數據ASID算法的流程圖。對于同層地層來說，可通過改變參考道和對比道的位置，利用其子波振幅譜信息獲取一系列Q值。但在實際資料Q值估計時，同層位上獲取的一系列Q值中，可能會存在一些Q值小于零或數值過大，這些不合理結果須剔除。然后再將有效Q值估計結果進行平均，獲得平均等效Q值，提高Q值估計的準確程度，弱化噪聲等環境因素對Q值估計的影響。

圖1 多層介質疊前地震數據ASID法Q值估計流程

2 模型試驗

2.1 影響因素分析

估算Q值時需選擇頻帶寬度和時窗寬度。另外，隨機干擾也會影響Q值估計。為此，對不同帶寬、不同時窗和不同噪聲環境進行了試驗，并將本文方法與LSR、CFS和LSAD三種常用方法做對比分析。

2.1.1 頻帶寬度試驗

為獲取準確的Q值，許多頻域Q值估計方法都需選取適當的頻帶寬度，過寬或過窄的頻帶都可能導致Q值估計效果不理想[18,25]。為了測試帶寬對ASID法的影響，選取主頻為50Hz的Ricker子波，分別模擬了Q=200和50時不同旅行時的一維衰減地震記錄(圖2)。相同Q值的地震子波旅行時分別為200ms和300ms。原始Ricker子波的寬度約為40ms，由于ASID算法是基于不同時刻地震子波振幅譜積分差值估算Q值的，這里給出了不同衰減程度的子波振幅譜(圖3)。可見隨著旅行時間的增大及Q值的減小，振幅譜的主頻向低頻方向移動，且振幅譜有效頻帶寬度減小。

由于頻率過低時振幅譜較小，因此選取起始頻率為10Hz，截止頻率從20Hz開始遞增，獲得了ASID法和其他三種常規方法的Q值估計結果隨截止頻率的變化(圖4)。Q值估計過程中，均使用60ms時窗拾取子波。從圖4可看出：截止頻率較小時，LSR法和CFS法的Q值估計誤差較大，隨著截止頻率增加，四種方法獲得Q值均逐漸接近理論值；在一定的中頻段內，Q值估計結果較穩定，即效果較好； 但截止頻率較大時，隨著截止頻率的增加，LSR法和LSAD法的估計效果明顯變差，而ASID法基本不受截止頻率的影響。由此可知： LSR法不適于截止頻帶較小和較大時的Q值估計，需選擇合適的頻段； CFS法不適于截止頻率較小時的Q值估計，不過該方法在高頻段時，其受頻率影響程度弱于LSR法和LSAD法； LSAD法則不適于高頻段的Q值估計； 而ASID法受頻率影響最小，尤其高頻處頻率變化基本不影響Q值估計。

為了探究ASID算法不易受頻率影響的原因，并了解泰勒級數展開對Q值估計的影響，圖5給出不同截止頻率下振幅譜積分差的理論值和泰勒級數展開近似、這兩者的差值，以及Q估算值。從圖5a和圖5b可見，隨著截止頻率的增加(即帶寬增大)，振幅譜積分差先迅速增加，后趨于平緩，這主要緣于高頻處不同時刻的地震子波振幅譜均趨于零(圖3)； 泰勒級數展開形式下的振幅譜積分差近似值與理論值基本一致，不存在明顯差值。從圖5c和圖5d可見，理論值與近似值之差也隨截止頻率的增加先迅速遞增，后平緩變化； 兩者雖存在一定差值，但幅度較小，且該差值也是Q值估算誤差的主要來源。從圖5e、圖5f可見，Q估計值與頻率關系不大，在截止頻率為30Hz時誤差最大，但相對誤差仍小于5%； 無論頻率多高，Q估算值基本無變化，相對誤差都在3%以內。上述分析表明，本文方法所用振幅譜積分差對高頻不敏感，高頻干擾對方法的影響不大； 雖使用二階泰勒級數展開，但它帶來的誤差較小；Q值估計與頻率關系不大，無論截止頻率設定得過高或過低，Q值估計誤差均較小。

圖2 不同Q值的合成衰減地震記錄

圖3 Q=200(a)和Q=50(b)的不同時刻地震子波的振幅譜

圖4 Q=200(a)和Q=50(b)時不同方法估計的Q值隨截止頻率的變化(起始頻率為10Hz)波浪線表示該段的值已超出取值范圍

圖5 不同時刻子波振幅譜積分差理論值與泰勒級數展開下的近似值(a)和(b)分別是Q=200和Q=50的振幅譜理論值與近似值； (c)和(d)分別是Q=200和Q=50的振幅譜理論值與近似值差值； (e)和(f)分別是Q=200和Q=50時ASID法的Q估算結果

2.1.2 時窗寬度試驗

為了厘清ASID算法受時窗寬度的影響程度，仍基于圖2的一維合成衰減地震記錄，選擇不同時窗寬度進行試驗。從不同時窗子波振幅譜(圖6)的整體上看，60ms與40ms的子波頻譜幾乎完全重合，即當選取完整子波時，其寬度對振幅譜影響較小。隨著時窗寬度的減小，振幅譜之間差異變大，低頻段能量不斷增加，且頻帶有所拓寬； 當時窗為10ms時，振幅譜峰值處于零頻點附近。在不同時窗寬度下，分別利用本文方法及LSR、CFS、LSAD方法估算Q值，為了使其他三種方法的Q估算值最佳，選定頻帶范圍均為10～100Hz，不同方法的Q估算值見圖7?？梢姰敃r窗寬度較小(10ms或20ms)時，LSR和CFS法估計的Q值遠大于真實值，LSAD和ASID法則能估算出更接近真實的Q值，且ASID法更精確(25ms、30ms)； 隨著時窗寬度的增加，四種方法均能較好地獲取真實Q值，尤其時窗寬度大于子波寬度時。該試驗表明，LSR法和CFS法并不適合小時窗(即非完整子波)下的Q值估計，而ASID法受時窗寬度影響最小。

2.1.3 噪聲干擾試驗

為了分析噪聲對Q值估計的影響，對圖2一維合成衰減地震道分別加入5%、10%和15%的高斯隨機噪聲(圖8)，對所得記錄均隨機運行1000次，得到三種噪聲環境下Q值估計概率分布圖(圖9和圖10)。在估算Q值過程中，選取的時窗寬度均為40ms，頻帶范圍仍為10～100Hz。需指出的是，在圖9和圖10中，Q估算值超出橫坐標取值范圍的，均統計在坐標軸的兩側，即小于0的統計在-40～0之內，超出橫坐標最大值的統計在最大值處，不過在計算平均值和均方差時使用的是實際Q估算值。

圖6 不同時窗(10～60ms)的子波振幅譜(Q=50)

圖7 不同時窗寬度下四種方法的Q值估計柱狀圖(Q=50)

從該兩圖可看出： LSR法受噪聲干擾影響最嚴重，Q估算值幾乎完全失真，即使噪聲含量較低(圖9a、圖10a，5%噪聲)時，其Q值估計誤差也很大； CFS法在噪聲含量較低、Q值較大(圖9b)時，Q估算值較可靠，但隨著噪聲增加、Q值減小，其計算誤差迅速增加，也難以有效估計Q值； LSAD法在Q值較大(圖9c、圖9g、圖9k)時，具有較強抗噪能力，能估算出較真實的Q值，在Q值較小(圖10c、圖10g、圖10k)時，其計算誤差隨噪聲環境的惡化而增加，已不能有效進行強噪聲、小Q值的估算(圖10k)； 相對于其他三種方法，ASID法抗干擾能力最強，能有效準確地估算Q值，且其概率分布更接近于正態分布，標準差也是最小，不過隨著噪聲含量的增加，誤差會有所增大。

圖8 分別添加5%(a)、10%(b)、15%(c)噪聲的一維地震記錄

該試驗表明，LSR法受噪聲影響最嚴重，即使少量噪聲也會導致較大Q值估算誤差； CFS法僅適用于Q值較大、噪聲含量較低時的Q值估計； LASD法不適用于強噪聲、小Q值的估計； ASID法受噪聲強度和Q值影響較小，能在不同噪聲環境中準確估算Q值。

圖9 Q=200時三種噪聲環境下四種方法Q值估計概率分布圖(a)～(d)5%噪聲下LSR法、CFS法、LSAD法、ASID法； (e)～(h)10%噪聲下LSR法、CFS法、LSAD法、ASID法；(i)～(l)15%噪聲下LSR法、CFS法、LSAD法、ASID法

2.2 含薄層的共炮點道集測試

為了驗證ASID法在薄層Q值估計中的效果，設計一個水平三層介質模型，模型參數如表1所示。

在合成地震記錄時，為了與實際相符，兼顧了反射系數和大地濾波作用的影響(圖11a)，震源采用主頻為50Hz的Ricker子波，道間距為20m，第一道檢波器位置與震源重合，共80道地震數據。需說明的是，在利用ASID算法和常規Q值估算法預測Q值之前，首先要獲得消除反射系數影響的地震記錄(圖11b)。設薄層(第二層)厚度為10m，確保小于薄層定義的最大厚度(子波寬約40ms，則薄層最大厚度為λa/2=VTa/2=1400m/s×40ms×0.5=28m)。

表1 三層水平介質模型參數

圖12是消除反射系數影響后的首道地震記錄，可見接收到的第二層頂、底界面的反射子波發生了明顯干涉，地震子波負邊瓣幅值相互疊加，幅值明顯增加，已與第二層地震子波的主峰幅值相當。實際地震記錄中常遇此情形，尤其對于油氣層。

由于此模型為多層模型，可采用圖1的流程完成ASID算法的Q值估計，其他算法則只需將圖1虛線部分處理步驟換成對應算法的Q值估計流程，即可實現不同算法的Q值估計。

圖10 Q=50時三種噪聲環境下四種方法Q值估計概率分布圖(a)～(d)5%噪聲下LSR法、CFS法、LSAD法、ASID法； (e)～(h)10%噪聲下LSR法、CFS法、LSAD法、ASID法； (i)～(l)15%噪聲下LSR法、CFS法、LSAD法、ASID法

圖11 合成衰減CSP道集(a)兼顧反射系數、大地濾波作用； (b)消除反射系數影響

為進一步提高方法的合理性，首先刪除LSR、CFS和LSAD法Q估算值中的負值和大于1000的數值，再將一系列Q值做平均化處理，得到等效Q值，最后利用常規層間Q值求取方法[11，18]獲取層間Q值(圖13)。與此同時，為避免子波干涉的影響，第一、第二層選取子波寬度為12ms，第三層選取整個子波寬度(40ms)，頻帶范圍仍選用10～100Hz。

從圖13可見，由于第一層選取非完整子波，導致LSR和CFS法在第一層上估計的Q值(分別為413.4和414.0)遠高于真實值，這與前述時窗寬度試驗結果相一致； 第二層同樣選取非完整子波，且第二層層間Q值計算需使用第一層Q值，由于第一層Q估算值已嚴重失真，導致第二層Q值為負值(分別為-91.4和-74.9)； 第三層Q值估計受到上兩層Q值估計不準確的影響，誤差同樣較大。LSAD法對含薄層的地震記錄Q值估計效果較好(圖13)，在三層介質上的Q估算值分別為216.1、34.3和145.9，相對誤差分別為8.06%、71.50%和27.08%，遠小于LSR和CFS法的估算結果。而ASID是受時窗寬度影響最小的方法(圖7)，故對薄層Q值的估計則更準確，它在三層介質上的Q估算值分別為200.6、20.4和206.0，相對誤差分別僅為0.28%、1.94%和2.98%，整體精度比LSAD法提高了近20倍。

圖12 圖11b的第一道地震數據

圖13 不同Q值估計方法計算結果

3 實例應用

選用M海域共中心點(CMP)地震道集(圖14)，分別應用ASID和LSAD法進行Q值估計試驗。時間采樣間隔為4ms。需補充說明的是，該組數據是已消除球面擴散及透射損失作用影響的疊前地震數據。通過對比同相軸一致的不同地震道，發現基本所有地層對應的地震子波橫向衰減均較快，表明地層Q值可能相對較小。

為了合理確定Q值估計時的頻段，從第一道地震記錄的振幅譜(圖15)可知，主頻約為21Hz，帶寬約為10～40Hz。在估算Q值時，選取時窗為28ms(包含完整子波)、頻帶處于10～40Hz?；趫D1的Q值估計流程，分別利用ASID和LSAD方法獲取等效Q值(圖16)，可見兩種方法估計的等效Q值相似度較高，說明所得等效Q值具有一定的可信度。該等效Q值在30～60區間，整體上偏小，反映了沉積層可能具有一定的黏度。

圖14 M海域CMP地震道集

圖15 第一道振幅譜

圖16 ASID(紅線)與LSAD(藍線)估算的等效Q值

4 結論

品質因子Q是用于提高地震記錄縱向分辨率和分析儲層特征的關鍵參數。本文在地震波振幅衰減項的二階泰勒級數展開基礎上，利用不同時刻地震子波振幅譜積分差值與Q值的關系，構建了一種Q值估計(ASID)方法。針對頻帶寬度、時窗寬度及噪聲干擾等的試驗結果表明，相對于LSR、CFS和LSAD法，ASID法受頻帶寬度和時窗寬度影響更小、抗干擾能力更強。含薄層的三層介質模型試驗表明ASID法具有更高的Q值估計精度，能更準確地估算薄層Q值。

還應說明，本文算法是一種基于疊前(CMP道集)數據Q值估計方法，在前期地震數據處理中，應預先做好振幅增益恢復、抽道集、透射損失和球面擴散補償等預處理，且這些預處理的質量直接影響Q值估計的準確程度。如增益恢復處理得不好，會導致無效的Q值估計結果；未做透射損失和球面擴散補償會導致Q估算值偏低。因此，在Q值估計之前，應充分做好數據處理的前期準備以獲取更準確的Q值。