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

(①東華理工大學(xué)地球物理與測控技術(shù)學(xué)院，江西南昌 330013； ②東華理工大學(xué)核科學(xué)與技術(shù)學(xué)院，江西南昌 330013； ③東華理工大學(xué)核資源與環(huán)境國家重點(diǎn)實(shí)驗(yàn)室，江西南昌 330013)

0 引言

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

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

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

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

1 方法原理

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

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

(1)

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

(2)

(3)

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

Δr=v(r,f0)Δt

(4)

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

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

(5)

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

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

(6)

從該式發(fā)現(xiàn)Q只影響實(shí)數(shù)部分，因此不同時(shí)刻地震子波振幅譜的關(guān)系[16,25]可表示為

(7)

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

(8)

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

(9)

基于二階泰勒級(jí)數(shù)展開，式(9)指數(shù)項(xiàng)可表示為

(10)

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

(11)

兩個(gè)接收點(diǎn)之間的振幅譜積分差值為

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

(12)

上式可簡化為

(13)

其中

雖然式(13)存在兩個(gè)解，但通過實(shí)驗(yàn)發(fā)現(xiàn)，其中一個(gè)解是一個(gè)非常小的數(shù)，即非Q值的真解，而Q值的有效解為

(14)

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

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

圖1 多層介質(zhì)疊前地震數(shù)據(jù)ASID法Q值估計(jì)流程

2 模型試驗(yàn)

2.1 影響因素分析

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

2.1.1 頻帶寬度試驗(yàn)

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

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

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

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

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

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

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

2.1.2 時(shí)窗寬度試驗(yàn)

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

2.1.3 噪聲干擾試驗(yàn)

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

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

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

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

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

該試驗(yàn)表明，LSR法受噪聲影響最嚴(yán)重，即使少量噪聲也會(huì)導(dǎo)致較大Q值估算誤差； CFS法僅適用于Q值較大、噪聲含量較低時(shí)的Q值估計(jì)； LASD法不適用于強(qiáng)噪聲、小Q值的估計(jì)； ASID法受噪聲強(qiáng)度和Q值影響較小，能在不同噪聲環(huán)境中準(zhǔn)確估算Q值。

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

2.2 含薄層的共炮點(diǎn)道集測試

為了驗(yàn)證ASID法在薄層Q值估計(jì)中的效果，設(shè)計(jì)一個(gè)水平三層介質(zhì)模型，模型參數(shù)如表1所示。

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

表1 三層水平介質(zhì)模型參數(shù)

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

由于此模型為多層模型，可采用圖1的流程完成ASID算法的Q值估計(jì)，其他算法則只需將圖1虛線部分處理步驟換成對(duì)應(yīng)算法的Q值估計(jì)流程，即可實(shí)現(xiàn)不同算法的Q值估計(jì)。

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

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

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

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

圖12 圖11b的第一道地震數(shù)據(jù)

圖13 不同Q值估計(jì)方法計(jì)算結(jié)果

3 實(shí)例應(yīng)用

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

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

圖14 M海域CMP地震道集

圖15 第一道振幅譜

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

4 結(jié)論

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

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