基于射線路徑的疊前高精度Q值估計(jì)方法

王德利，戴建芳

（吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，吉林長(zhǎng)春130026）

地震波傳播過程中存在著能量衰減。引起地震波衰減的因素很多，總的來說可分為兩大類［1］：一類與地震波傳播特性有關(guān)，如球面擴(kuò)散、界面透射損失、層間多次波、振幅調(diào)諧引起的衰減、波的相位轉(zhuǎn)換導(dǎo)致的衰減等；另一類衰減與儲(chǔ)層巖石物理參數(shù)的變化相關(guān)，液體的性質(zhì)（如粘度、可壓縮性）及其飽和度，孔隙的性質(zhì)（包括孔隙度、滲透率）以及裂隙的幾何參數(shù)都會(huì)對(duì)這類衰減產(chǎn)生影響。對(duì)于后一類衰減我們通常用品質(zhì)因子Q來衡量。正確估計(jì)Q值對(duì)于提高地震記錄的分辨率及對(duì)儲(chǔ)層的預(yù)測(cè)與描述具有重要意義。

提取Q值的方法有很多種［2］，包括振幅衰減法、子波模擬法、譜模擬法、上升時(shí)間法、解析信號(hào)法和質(zhì)心頻移法等。Tonn［3］認(rèn)為這幾種方法都有一定的適用性，至今沒有任何一種方法是普遍適用的?，F(xiàn)在實(shí)際應(yīng)用中最常用的方法是譜比法，這是一種頻率域方法，其基本思路是分別截取目標(biāo)地層頂、底面一個(gè)時(shí)窗內(nèi)的地震記錄，通過線性擬合頻率和振幅之間的關(guān)系得出Q值［4］。

許多因素使得衰減測(cè)定的準(zhǔn)確度降低，這主要源于地震譜中不同波至走時(shí)不同，計(jì)算出的衰減受相鄰層位影響很大。Dasgupta等［5］提出了一種利用振幅值隨炮檢距的變化從CMP道集中提取Q值的方法；Hackert等［6］對(duì)該方法進(jìn)行了改進(jìn)，提出了對(duì)振幅譜進(jìn)行校正后提取Q值的方法（LSR）；王小杰等［7］使用兩步擬合算法（QVO 方法）通過偏移距的歸零處理進(jìn)一步提高Q值提取的精度。這些方法的有效性已經(jīng)被證明，包括巖性的識(shí)別［8］和裂縫的方位識(shí)別［9］。但是，在 QVO方法的實(shí)現(xiàn)過程中，時(shí)窗內(nèi)的地震記錄會(huì)受調(diào)諧效應(yīng)和多次波反射的影響，使得Q值提取精度不高；另外，計(jì)算過程中所用到的時(shí)差并非僅僅是目標(biāo)層的旅行時(shí)。由于介質(zhì)的各向異性以及震源與檢波器的方向性的影響［10－11］，使得目標(biāo)層的衰減受相鄰層位的影響，所以QVO方法第2步擬合中的歸零處理求取的零偏移距上的頻譜比斜率并不精確。

考慮到疊前地震資料比疊后地震資料的振幅和旅行時(shí)信息更加豐富，更能反映地層的一些細(xì)微特征，我們?cè)赒VO方法的基礎(chǔ)上提出一種基于射線路徑的疊前地震數(shù)據(jù)高精度Q值估計(jì)方法，其核心是通過τ－p變換對(duì)同一水平慢度的道集進(jìn)行處理，采用廣義S變換計(jì)算地層頂、底面反射波的瞬時(shí)頻譜，從而克服傳統(tǒng)Q值估計(jì)方法的兩步擬合所帶來的統(tǒng)計(jì)性誤差，以提高Q值計(jì)算的精度。

1 基本理論

1.1 常規(guī)Q值計(jì)算方法

振幅譜為S（f）的地震波在層狀介質(zhì)中傳播時(shí)，不僅有地層的粘彈性所引起的e指數(shù)衰減，還有透射和反射的能量損失P和球面擴(kuò)散能量損失G，而且這兩項(xiàng)與頻率無關(guān)。對(duì)于一個(gè)初始時(shí)刻振幅譜為S0（f）的地震子波，傳播t時(shí)間后，其振幅譜S（f）表示為

式中：Q是傳播過程中的平均衰減效應(yīng)。

同一道來自同一地層頂、底面反射波的瞬時(shí)振幅譜可分別表示為

式中：S1（f）為接收到的地層頂面反射的頻譜；S2（f）為接收到的地層底面反射的頻譜；t21為接收到的地層底面反射在地層1中的旅行時(shí)；t22為其在地層2中的旅行時(shí)（圖1）。

圖1 同一道接收到的地層2頂、底面的反射

（3）式除以（2）式且兩邊取對(duì)數(shù)得

通常的處理方法是假設(shè)接收到的同一地層頂、底面反射在地層1中的旅行時(shí)是一樣的，即t21＝t1，從而得到

1.2 高精度Q值計(jì)算方法

為了更正確地計(jì)算走時(shí)，解決地層各向異性和震源、檢波器的方向性問題［12］，疊前高精度Q值計(jì)算方法是把數(shù)據(jù)變換到τ－p域，對(duì)同一水平慢度的道集進(jìn)行處理，以避免傳統(tǒng)方法的兩步擬合所帶來的統(tǒng)計(jì)性不精確；并且引入走時(shí)各向異性參數(shù)，使計(jì)算出的衰減僅僅由目標(biāo)層位引起。

在以上常規(guī)Q值計(jì)算方法中做了一個(gè)假設(shè)，即t21＝t1，然而由圖1可知事實(shí)并不是這樣的，故計(jì)算出的Q值有一定的誤差。而在τ－p域中，同一水平慢度的道集接收到的地層2的反射波路徑在地層1中的走時(shí)是相同的，如圖2所示，由此可知t21＝t1。

圖2 同一水平慢度道集接收到的地層2的反射波路徑

為了把計(jì)算過程中相鄰層位的干涉效應(yīng)減到最小，疊前高精度Q值計(jì)算方法對(duì)振幅譜的計(jì)算是利用時(shí)頻分布圖［13］恰當(dāng)?shù)赜?jì)算地震波到達(dá)地層頂、底面的瞬時(shí)頻譜。S變換是一種介于短時(shí)傅里葉變換和小波變換之間的對(duì)非平穩(wěn)信號(hào)處理的時(shí)頻分析技術(shù)［14］，廣義S變換是在S變換的基礎(chǔ)上發(fā)展而來［15］，其窗函數(shù)加入了調(diào)節(jié)系數(shù)，是一種分辨率比較高的時(shí)頻分析方法。利用廣義S變換的變化時(shí)窗對(duì)地層上、下界面對(duì)應(yīng)的瞬時(shí)頻譜進(jìn)行分析，能夠克服傳統(tǒng)Q值計(jì)算方法中時(shí)窗長(zhǎng)度和寬度選擇對(duì)計(jì)算結(jié)果的影響。

1.2.1 τ－p 變換

對(duì)于同一水平慢度的道集數(shù)據(jù)，Snell定律保證了在水平各向同性速度和衰減的假設(shè)下其在地層1中的反射沿著同一路徑。Behura等［16］指出，具有同一水平慢度的反射在t－x域中表示為具有相同瞬時(shí)斜率的不同反射。同時(shí)，在把數(shù)據(jù)變換到τ－p域的過程由于綜合效應(yīng)能壓制部分隨機(jī)噪聲。

二維空間的τ－p變換的數(shù)學(xué)表達(dá)式為［17－18］

對(duì)于離散數(shù)據(jù)有

Van der Baan等［19］指出，對(duì)于τ－p 域 一個(gè)給定的層，其校正量τi為

式中：τ0是垂直入射波的旅行時(shí)；υi是層速度；ηi是Alkhalifah等［20］定義的各向異性參數(shù)。一個(gè)水平層的總校正量為

τ－p校正曲線表示的僅是平面波的垂直傳播時(shí)間部分。然而，為了提取Q值，需知總的傳播時(shí)間。由τ－p變換公式t＝τ＋px和（9）式可知

將（8）式代入（10）式有

由此，對(duì)于τ－p域中任意給定的一個(gè)點(diǎn)，給定一個(gè)p值，運(yùn)用（11）式和（12）式就可以找到t－x域中對(duì)應(yīng)的點(diǎn)［21］，這樣就可以在t－x域追蹤到圖2中相同水平斜率的反射波。

（11）式和（12）式是針對(duì)各向異性走時(shí)方程建立的，其中含有各向異性參數(shù)ηi，所以，我們提出的方法對(duì)各向異性介質(zhì)具有適用性。如果僅對(duì)各向同性介質(zhì)數(shù)據(jù)進(jìn)行處理，可默認(rèn)ηi為0。

1.2.2 廣義S變換

廣義S變換是由S變換發(fā)展而來。地震波在地層中傳播時(shí)，高頻成分要比低頻成分衰減嚴(yán)重，主頻變低；S變換的窗函數(shù)固定不變，不能隨著信號(hào)本身的特點(diǎn)進(jìn)行調(diào)整。我們采用了一種改進(jìn)的廣義S變換方法，其變化的時(shí)窗可以滿足吸收衰減分析的需求，即

當(dāng)λa和p變化時(shí)，廣義S變換的窗函數(shù)形狀隨頻率而變化。具體來說，窗函數(shù)的寬度隨頻率呈反比變化，幅值與頻率呈正比變化。低頻段時(shí)窗較寬，幅值較低，具有較高的頻率分辨率；高頻段時(shí)窗較窄，幅值較高，可獲得較高的時(shí)間分辨率［22］。

1.2.3 算法實(shí)現(xiàn)步驟

圖3給出了疊前高精度Q值估計(jì)方法的計(jì)算流程。如圖3所示，對(duì)一個(gè)疊前CMP道集數(shù)據(jù)，在給定的水平慢度范圍內(nèi)依照設(shè)置的間隔選取p1，p2，…，pn，依照（11）式和（12）式，把一個(gè)恒定水平慢度pi產(chǎn)生的反射（pi，τi）映射到t－x 域（ti，xi），公式中的τ0i，υi，ηi由高精度速度分析得到［23］，也就是選取目標(biāo)層位頂、底面在相同的表層路徑（即恒定pi）下所對(duì)應(yīng)的道集xi1，xi2。這里僅對(duì)各向同性介質(zhì)進(jìn)行處理，即

圖3 疊前高精度Q值估計(jì)方法計(jì)算流程

2 模型數(shù)據(jù)測(cè)試

對(duì)圖4a所示的4層地質(zhì)模型（各層的Q值依次為100，50，100，150）采用射線追蹤法進(jìn)行正演模擬。模擬采用的道間距為20m，炮檢距為0～1 180m；時(shí)間采樣間隔為1ms；所用子波為50Hz的雷克子波。模擬出的含衰減零偏移距CMP地震道集記錄見圖4b。

將圖4b的CMP道集記錄變換到τ－p域（圖5），地震記錄中接收到的第1層的反射波衰減程度僅由第1層的走時(shí)決定，故可用t－x域中不同地震道之間的頻譜比斜率值反演求取該層Q值。對(duì)第2層、第3層和第4層分別選取的水平慢度范圍為0～0.3，0～0.2，0～0.1，間隔都為0.001；采用本文提出的高精度方法計(jì)算這3層的Q值。

圖4 地質(zhì)模型（a）及其正演模擬地震記錄（b）

圖5 轉(zhuǎn)換到τ－p域中的地震記錄

表1給出了采用本文方法與用QVO（兩步擬合算法）方法計(jì)算的結(jié)果及其誤差比較，可以明顯看出本文方法的精確性。

為了驗(yàn)證本文方法的抗噪性，給上述模型的理論數(shù)據(jù)加入30dB的隨機(jī)噪聲，結(jié)果如圖6所示。圖7是水平慢度為0.09的頻譜。圖8給出了第2層、第3層的振幅譜比值與頻率的擬合關(guān)系。

表2給出了加入噪聲后用本文方法與用QVO方法計(jì)算的結(jié)果及其誤差，可見加入噪聲后本文方法仍保持較高的精度。比較加入噪聲前（表1）、后（表2）的誤差可以看出，本文方法抗噪性較強(qiáng)，加入噪聲后計(jì)算的Q值波動(dòng)不大。

表1 兩種方法計(jì)算出的Q值及其誤差

表2 加入30dB噪聲后兩種方法計(jì)算出的Q值及其誤差

圖6 加入30dB噪聲后的t－x 域（a）和τ－p 域（b）地震記錄

圖7 p＝0.09時(shí)模型4個(gè)界面的反射波頻譜

圖8 模型第2層（a）和第3層（b）的振幅譜比值與頻率的擬合關(guān)系

3 實(shí)際資料處理

為了驗(yàn)證本文方法在實(shí)際應(yīng)用中的可行性和有效性，選取了一個(gè)實(shí)際CMP道集數(shù)據(jù)進(jìn)行Q值試算。圖9為經(jīng)過去噪預(yù)處理之后的CMP道集，經(jīng)高精度速度分析后掃描出t0和vNMO，解釋出反射層位，圖9中A，B，C，D，E為各層反射界面。

為了驗(yàn)證本文方法的精確性，用計(jì)算出的Q值對(duì)CMP道集數(shù)據(jù)作反Q濾波處理。圖10顯示了截取的一段反Q濾波前、后的數(shù)據(jù)，通過對(duì)比可以看出，處理后的數(shù)據(jù)縱向分辨率有了較明顯的提高，很多層間小薄層經(jīng)過反Q濾波后能明顯分辨。可見本文方法在實(shí)際數(shù)據(jù)應(yīng)用中的效果較好，可用來做精確的反Q濾波處理，以提高地震記錄的分辨率。

圖9 經(jīng)去噪預(yù)處理后的實(shí)際CMP道集

圖10 反Q濾波前（a）、后（b）的一段道集記錄

4 結(jié)束語

我們提出的疊前高精度Q值估計(jì)方法在理論上具有一定的先進(jìn)性，相對(duì)于傳統(tǒng)譜比法的兩步擬合計(jì)算，疊前高精度Q值估計(jì)方法作出了兩個(gè)重要的改進(jìn)：①為使參考信號(hào)和測(cè)量信號(hào)具有相同的表層傳播路徑，對(duì)同一水平慢度的道集進(jìn)行處理，避免了傳統(tǒng)方法的兩步擬合所帶來的統(tǒng)計(jì)性不精確；②采用變時(shí)窗的時(shí)頻變換方法，最大限度地減少相鄰層位的干擾效應(yīng)。

通過理論模型數(shù)據(jù)測(cè)試證明了本文方法的精確性，實(shí)際數(shù)據(jù)的應(yīng)用結(jié)果也證明了該方法的有效性。

［1］李生杰，施行覺，王寶善，等.地層衰減在地震記錄上的特征分析［J］.石油地球物理勘探，2002，37（3）：248－253 Li S J，Shi X J，Wang B S，et al.Characteristics analysis of formation attenuation in seismic record［J］.Oil Geophysical Prospecting，2002，37（3）：248－253

［2］魏文，王小杰，李紅梅.基于疊前道集小波域Q值求取方法研究［J］.石油物探，2011，50（4）：355－360 Wei W，Wang X J，Li H M.Analysis of Qcalculating method based on the prestack gather in wavelet domain［J］.Geophysical Prospecting for Petroleum，2011，50（4）：355－360

［3］Tonn R.The determination of the seismic quality factor Qfrom VSP data：a comparison of different computational methods［J］.Geophysical Prospecting，1991，39（1）：1－27

［4］付勛勛，徐峰，秦啟榮，等.基于改進(jìn)的廣義S變換求取地層品質(zhì)因子Q值［J］.石油地球物理勘探，2012，47（3）：451－467 Fu X X，Xu F，Qin Q R，et al.Study on quality factors Qby modified generalized Stransform［J］.Oil Geophysical Prospecting，2012，47（3）：451－467

［5］Dasgupta R，Clark R A.Estimation of Qfrom surface seismic reflection data［J］.Geophysics，1998，63（6）：2120－2128

［6］Hackert C L，Parra J O.Improving Qestimation from seismic reflection data using well－log－based localized spectral correct ion［J］.Geophysics，2004，69（6）：1521－1529

［7］王小杰，印興耀，吳國(guó)忱.基于疊前地震數(shù)據(jù)的地層Q 值估計(jì)［J］.石油地球物理勘探，2011，46（3）：423－428 Wang X J，Yin X Y，Wu G C.Estimation of quality factors Qfrom pre－stack data［J］.Oil Geophysical Prospecting，2011，46（3）：423－428

［8］Clark R A，Carter A J，Nevill P C，et al.Attenuation measurements from surface seismic data—azimuthal variation and timelapse case studies［J］.Expanded Abstracts of 63rdEAGE Annual Conference，2001，L28

［9］Clark R，Benson P M，Carter A J，et al.Anisotropic P－wave attenuation measured from a multi－azimuth surface seismic reflection survey［J］.Geophysical Prospecting，2009，57（5）：835－845

［10］Reine C，Clark R，van der Baan M.Robust prestack Q－determination using surface seismic data：part 1—method and synthetic examples［J］.Geophysics，2012，77（1）：R45－R56

［11］Reine C，Clark R，van der Baan M.Robust prestack Q－determination using surface seismic data：part 2—3Dcase study［J］.Geophysics，2012，77（1）：B1－B10

［12］Reine C.A robust prestack Q－inversion in theτ－p domain using variable window spectral estimates［D］.England：University of Leeds，2009

［13］Reine C，van der Baan M，Clark R.The robustness of seismic attenuation measurements using fixed－and variable－window timefrequency transforms［J］.Geophysics，2009，74（2）：WA123－WA135

［14］劉喜武.基于廣義S變換研究地震地層特征［J］.地球物理學(xué)進(jìn)展，2006，21（2）：215－218 Liu X W.Study on characteristics of seismic stratigrephy by generalized S transform［J］.Progress in Geophysics，2006，21（2）：215－218

［15］Pinnegar C R，Mansinha L.The S transform with windows of arbitrary and varing shape［J］.Geophysics，2003，68（1）：381－385

［16］Behura J，Tsvankin I.Estimation of interval anisotropic attenuation from reflection data［J］.Geophysics，2009，74（6）：A69－A74

［17］Chapman C H.Generalized radon transforms and slant stacks［J］.Geophysical Journal of the Royal Astronomical Society，1981，66（2）：445－453

［18］Schultz P S，Claerbout J F.Velocity estimation and downward continuation by wavefront synthesis［J］.Geophysics，1978，43（4）：691－714

［19］Van der Baan M，Kendall J M.Estimating anisotropy parameters and traveltimes in theτ－p domain［J］.Geophysics，2002，67（4）：1076－1086

［20］Alkhalifah T，Tsvankin I.Velocity analysis for transversely isotropic media［J］.Geophysics，1995，60（5）：1550－1566

［21］Van der Baan M.Processing of anisotropic data in theτ－p domain：I—geometric spreading and moveout corrections［J］.Geophysics，2004，69（3）：719－730

［22］Stockwell R G，Mansinha L，Lowe R P.Localization of the complex spectrum：the S transform ［J］.IEEE Transactions on Signal Processing，1996，44（4）：998－1001

［23］Tsvankin H D，Thomsen I.Nonhyperbolic reflection moveout in anisotropic media［J］.Geophysics，1994，59（8）：1290－1304