金子奇,孫贊東
(中國石油大學(xué)(北京),北京102249)
根據(jù)應(yīng)用空間的不同,求取層間Q值的方法可以分成兩類,即時間域Q值估計方法和頻率域Q值估計方法。時間域Q值估計方法很難區(qū)分吸收衰減和其它類型的衰減,不建議使用。依據(jù)固有衰減作用與頻率相關(guān)的特點(diǎn),頻率域Q值估計方法能夠很好地解決這個問題。通過拾取不同時間位置的兩個子波并進(jìn)行傅里葉變換,頻率域Q值估計方法能夠根據(jù)子波振幅譜的變化來估計Q值。DASGUPTA等[1]首次利用譜比法估計Q值,在疊前反射波道集中通過建立子波對數(shù)譜和Q值間的關(guān)系估計Q值。隨后,許多學(xué)者對其做了改進(jìn),陳文爽等[2]、王小杰等[3]、王德利等[4]以S變換的頻譜作為譜比法的輸入,避免了通常計算地層吸收參數(shù)時的平均效應(yīng)。張繁昌等[5]將自適應(yīng)分解得到的子波應(yīng)用譜比法解決薄層和反射波干涉問題。李偉娜等[6]在微測井記錄中應(yīng)用譜比法并借鑒QVO(Q versus offset)的思想以適應(yīng)Q值橫向變化的情況。QUAN等[7]在VSP資料中,最先利用質(zhì)心頻率位移量估算Q值,建立了衰減后的子波峰值頻移量與Q值之間的關(guān)系。與質(zhì)心頻移法類似,ZHANG等[8]在疊前CMP道集中采用峰值頻率法估計Q值。隨后,許多學(xué)者對上述兩類頻移方法做了改進(jìn),王宗俊[9]利用譜模擬的子波對質(zhì)心頻移法進(jìn)行約束;張立彬等[10]引入積分中值參變量避免了質(zhì)心頻移法中對震源子波的先驗假設(shè);高靜懷等[11]利用特征結(jié)構(gòu)法提高地震子波峰值頻率的估計精度。
然而,許多因素影響頻域Q值估計方法的精度。震源強(qiáng)度和帶寬會影響Q值估計結(jié)果,在主頻范圍之外估計得到的Q值會變得不穩(wěn)定[12]。Q值估計結(jié)果隨著信噪比的降低而顯著畸變。對于一些頻率成分,遠(yuǎn)偏移距檢波器記錄的振幅有時會比近偏移距檢波器記錄的振幅強(qiáng),這種效應(yīng)會引起Q值估計結(jié)果的誤差[12],最嚴(yán)重的缺點(diǎn)是這些方法計算得到的Q值與路徑有關(guān)[13-14],Q值估計結(jié)果隨觀測系統(tǒng)的變化而變化。而真實地下介質(zhì)Q值分布應(yīng)僅與構(gòu)造、巖性和流體有關(guān),與觀測系統(tǒng)無關(guān)。另外,頻域Q值估計方法得到的是大套地層層間的等效Q值,無論是縱向上還是橫向上,分辨率都較低。BRZOSTOWSKI等[15]建立了不同傳播距離地震波振幅的變化與Q值之間的關(guān)系,采用聯(lián)合迭代重建技術(shù)SIRT求解方程實現(xiàn)振幅衰減Q值層析。振幅衰減層析方法在時間域?qū)崿F(xiàn),由于地震波振幅受多種因素(幾何擴(kuò)散、儀器響應(yīng)、震源/檢波器耦合特性、反射/透射等)的影響,因此其穩(wěn)定性較差。LIAO等[16]和嚴(yán)又生等[17]利用質(zhì)心頻移法結(jié)合層析反演方法計算Q值。該方法利用波傳播路徑下Q值與質(zhì)心位移量的對應(yīng)關(guān)系來迭代反演速度和Q值信息,對絕對振幅信息不敏感,只適應(yīng)較窄頻帶范圍內(nèi)的數(shù)據(jù)。
天然地震學(xué)中,NOWACK等[18]將衰減旅行時應(yīng)用于地震折射波數(shù)據(jù),反演淺層地殼內(nèi)部衰減及速度結(jié)構(gòu)。樊計昌等[19]利用三維Q值層析成像方法對岫巖坑的地震波衰減結(jié)構(gòu)進(jìn)行了研究。隨后,勘探地震學(xué)中也開始利用衰減旅行時層析方法反演地層Q值,CAVALCA等[20]、XIN等[21]和JIN[22]將衰減層析方法應(yīng)用于地震反射波數(shù)據(jù),利用對數(shù)譜比信息擬合衰減旅行時,并利用非線性反演方法計算地下真實Q值分布。衰減旅行時層析方法反演地下真實Q值分布,能適應(yīng)Q值剖面縱橫向變化,且對絕對振幅信息不敏感。衰減旅行時層析法,分別在實際資料和初始Q值模型中,求取實際衰減旅行時和模型衰減旅行時,建立二者之差與Q值模型更新量之間的關(guān)系,通過非線性反演的方法,計算Q模型更新量,更新初始模型直到接近真實地下Q值分布[22]。因此,衰減旅行時的計算是該方法的核心。傳統(tǒng)方法在計算實際衰減旅行時的時候,需要人工拾取最優(yōu)頻帶,人為誤差對Q值反演結(jié)果造成很大影響,并且,模型衰減旅行時的計算精度受射線追蹤精度影響,有偏差的旅行時和不合理的傳播路徑會導(dǎo)致反演結(jié)果產(chǎn)生誤差。魏文等[23]和王小杰等[24]采用時頻分析結(jié)合譜比法解決人為拾取最優(yōu)頻帶的問題,但該方法仍屬于頻率域Q值估計方法,得到的是大套地層的等效Q值。
因此,為了得到更為可靠和準(zhǔn)確的Q值估計結(jié)果,本文采用改進(jìn)后的衰減旅行時層析方法反演地下Q值分布。該方法在實際資料中拾取衰減旅行時時,采用對數(shù)譜比反演方法,利用多道反射波信息同時反演多道衰減旅行時,避免了傳統(tǒng)方法中人為選擇最優(yōu)頻帶可能引入的誤差,反演結(jié)果更穩(wěn)定。在計算模型衰減旅行時過程中,為提高衰減旅行時計算精度,采用方向梯度射線追蹤方法計算旅行時。改進(jìn)的算法用方向速度梯度代替?zhèn)鹘y(tǒng)的速度梯度,考慮了不同射線在不同傳播方向上的速度差異,得到更加準(zhǔn)確的旅行時和傳播路徑。而且,通過在反演過程中考慮反射波在不同界面內(nèi)的旅行時之差,可以有效解決上覆地層的影響。改進(jìn)后的衰減旅行時層析方法,經(jīng)模型數(shù)據(jù)和實際資料測試,Q值反演結(jié)果更為準(zhǔn)確,為后續(xù)吸收補(bǔ)償和成像提供可靠的衰減信息。
在地震勘探中,震源激發(fā)的地震波經(jīng)過地下介質(zhì)傳播,經(jīng)反射層反射回地面,再由地面檢波器接收。地震波本身攜帶了振幅和旅行時信息,即地震波的動力學(xué)信息和運(yùn)動學(xué)信息。我們利用旅行時信息進(jìn)行層析反演求取地下Q值分布。衰減旅行時層析反演求取Q值流程如圖1所示,具體算法可以分為以下幾步:①對數(shù)頻譜法反演方法計算炮集記錄上的反射波衰減旅行時;②建立初始Q值模型,利用方向梯度射線追蹤模擬地震波傳播,計算模型中的反射波衰減旅行時;③建立實際衰減旅行時與模型衰減旅行時之差與Q值模型更新量的關(guān)系,采用聯(lián)合迭代重建方法求解大型稀疏矩陣,更新Q值模型。
圖1 衰減旅行時層析反演求取Q值流程
圖2為單層二維介質(zhì)模型,由震源激發(fā)的射線經(jīng)反射后到達(dá)地表被接收。衰減旅行時層析方法,首先將地下介質(zhì)劃分成網(wǎng)格(圖2),并在每一網(wǎng)格內(nèi)定義獨(dú)立的速度和Q值,速度和Q值均為常量。一條由震源激發(fā)經(jīng)模型底界面反射并由地表檢波器接收的射線,在網(wǎng)格中被分為多段。
圖2 衰減旅行時層析反演中的網(wǎng)格劃分示意
第k條射線的衰減旅行時定義為[22]:
(1)
式中:lkl,tkl,vkl分別為第k條射線在第l個網(wǎng)格中的傳播距離、傳播時間和傳播速度。
將方程(1)改寫為矩陣形式:
(2)
式中:矩陣T包含由射線追蹤得到的每一網(wǎng)格內(nèi)的傳播時間tkl;Q為包含Q值的向量;t*向量為衰減旅行時。
(3)
(4)
公式(4)為m個方程n個未知數(shù)組成的大型稀疏矩陣,應(yīng)用正交分解最小二乘法(LSQR)求解該方程組。最終Q值結(jié)果Qf由初始模型和模型更新量共同表示:
(5)
傳統(tǒng)算法中,利用中心差分公式計算網(wǎng)格內(nèi)的速度梯度:
式中:下標(biāo)i,j分別表示在x方向和y方向的網(wǎng)格索引;v為定義的網(wǎng)格內(nèi)速度參數(shù)。
當(dāng)射線穿過分界面時,由于不準(zhǔn)確的速度梯度定義,傳統(tǒng)算法計算得到的射線方向會偏離真實射線路徑,尤其是當(dāng)網(wǎng)格間的速度變化超過10%時,偏離誤差更加明顯,使得傳播時間的計算也隨之偏離真實值。射線傳播距離s后的傳播時間定義為:
(7)
式中:向量n0為射線出射位置s=0處單位向量,指示射線傳播方向;向量λ=(λx,λy)為沿射線方向的速度梯度;v0為原點(diǎn)處的速度;o(λ3)為包含了λ的三階及更高階項。將旅行時結(jié)果代入公式(1)可得衰減旅行時,旅行時的誤差最終導(dǎo)致衰減旅行時的計算誤差。
改變傳統(tǒng)算法中對速度梯度的定義,考慮速度隨傳播方向的變化,用方向速度梯度代替公式(7)中的速度梯度。圖3展示了4種不同的射線傳播方向,重新定義速度梯度如下。
圖3 不同射線傳播方向示意a 射線沿西北方向傳播; b 射線沿東北方向傳播; c 射線沿東南方向傳播; d 射線沿西南方向傳播
1) 射線沿西北方向傳播時(圖3a):
2) 射線沿東北方向傳播時(圖3b):
3) 射線沿東南方向傳播時(圖3c):
4) 射線沿西南方向傳播時(圖3d):
根據(jù)對不同射線出射方向的方向梯度定義,利用公式(7)計算網(wǎng)格內(nèi)射線旅行時。
(12)
式中:f為地震信號頻率;A(t1,f),A(t2,f)分別為地震波在傳播時間t1和t2處的振幅譜;G為頻率獨(dú)立因子,代表與頻率無關(guān)的振幅衰減;B為獨(dú)立因子項;Δt*為兩個反射波衰減旅行時之差。傳統(tǒng)層析方法通過擬合對數(shù)譜比與頻率得到斜率,建立斜率與衰減旅行時的關(guān)系式,求得衰減旅行時。擬合斜率時,由于無法準(zhǔn)確選取最優(yōu)擬合頻帶,且無法消除噪聲等影響,通常會引入人為誤差。為解決上述問題,在共炮點(diǎn)道集中采用多道反射波同時計算多道衰減旅行時。同時反演頻率獨(dú)立因子G和Q值,可以壓制非吸收衰減作用(如球面擴(kuò)散等)的干擾。而且,傳統(tǒng)方法得到的是衰減旅行時之差。由于拾取不同傳播時間位置的反射波,若拾取的子波在上覆地層具有不同的傳播路徑,上覆地層內(nèi)的衰減影響則不能相互抵消,從而影響Q值估計結(jié)果。在反演過程中考慮反射波在不同界面內(nèi)的旅行時之差,可以有效解決上覆地層的影響。
以兩層介質(zhì)模型為例,沿偏移距方向,拾取t1時刻多個反射波A1,A3,…,A2n-1和t2時刻多個反射波A2,A4,…,A2n,如圖4所示。
圖4 兩層介質(zhì)模型內(nèi)不同界面反射波射線路徑示意(這些反射波均來自同一炮集記錄)
拾取的反射波經(jīng)傅里葉變換得到其振幅譜,整理并分別沿縱向和橫向取對數(shù)可得:
(13)
其中,b1,b2,…,b2n中包含對數(shù)譜比信息。將公式(13) 表示為d=Fm的矩陣形式:
(14)
(15)
得到的衰減旅行時。求解(15)式可得模型更新量,代入公式(5)求得最終Q值分布。
建立4層介質(zhì)模型,測試改進(jìn)后的層析反演方法估計Q值的準(zhǔn)確性,模型參數(shù)如表1所示。震源子波選擇主頻為60Hz的雷克子波,射線追蹤計算傳播路徑和旅行時,并計算振幅衰減項和相位畸變項,得到合成共炮點(diǎn)道集(為避免動校正拉伸作用對子波拾取的影響,正演過程中未考慮正常時差的變化),如圖5所示。
以第1層和第2層為例,說明改進(jìn)后的射線追蹤算法在射線路徑計算精度上的改善。將模型劃分網(wǎng)格,每一網(wǎng)格大小為10m×10m,模型尺度為80×40。對比改進(jìn)算法與傳統(tǒng)算法得到的折射角,結(jié)果如表2所示。圖6為兩層介質(zhì)模型下采用改進(jìn)前、后的射線追蹤算法計算得到的射線路徑對比結(jié)果。從圖6 中可以看出,兩種傳播路徑有明顯差異,射線傳播至200m深處界面時發(fā)生折射,且改進(jìn)后的算法由于考慮方向梯度的變化,得到的結(jié)果更接近真實射線路徑,提高了射線路徑的計算精度,尤其對于大角度入射的射線改進(jìn)更加明顯。與傳統(tǒng)算法誤差相比,在小入射角時(小于20°),改進(jìn)算法得到的折射角誤差小于傳統(tǒng)算法誤差的2%,在大入射角時(大于50°),誤差之比也不會超過10%。
表1 水平層狀介質(zhì)模型參數(shù)
圖5 水平層狀模型合成炮集記錄
表2 兩層介質(zhì)模型中改進(jìn)前后射線追蹤算法得到折射角對比
圖6 兩層介質(zhì)模型下改進(jìn)前、后射線追蹤算法計算得到的射線路徑對比
選取一炮地震記錄,拾取來自不同界面的反射波同相軸,分別采用多道同時反演法和傳統(tǒng)方法計算衰減旅行時,以第1層反射為例,計算結(jié)果如圖7所示,可以看出,與傳統(tǒng)方法相比,改進(jìn)后的方法更加穩(wěn)定且計算結(jié)果更加精確。
圖7 衰減旅行時計算結(jié)果
圖8 水平層狀模型中Q值反演結(jié)果
表3 Q值橫向變化介質(zhì)模型參數(shù)
圖9 Q值橫向變化模型中Q值反演結(jié)果
將衰減旅行時層析法應(yīng)用于實際資料求取Q值分布。我們將某一地區(qū)深度超過8000m的地層作為目標(biāo)進(jìn)行層析計算,地震剖面中對應(yīng)同向軸位置位于時間4.0s。抽取共炮點(diǎn)道集,選擇反射波同向軸明顯的位置3.2s和4.0s計算衰減旅行時(圖10),同時
在初始Q值剖面上射線追蹤得到相應(yīng)的正演衰減旅行時信息。初始Q值模型中Q設(shè)置為常數(shù)100,目標(biāo)區(qū)域以80m×80m進(jìn)行網(wǎng)格劃分,劃分網(wǎng)格后的模型橫向上包含80網(wǎng)格、縱向上包含80網(wǎng)格。最終反演Q值分布如圖11所示。
利用層析反演法求得的Q值分布,對地震數(shù)據(jù)的衰減做Q偏移補(bǔ)償。圖12a和圖12b分別以傳統(tǒng)的和改進(jìn)后的衰減旅行時層析求得的Q值作為輸入進(jìn)行Q補(bǔ)償偏移得到的CRP道集。對比可以看出,圖12b的道集質(zhì)量得到進(jìn)一步提高,振幅能量得到更好恢復(fù),相位畸變得到更合理校正,尤其3.5s位置改善十分明顯;同相軸更清晰,連續(xù)性增強(qiáng),原本無法識別的薄層得以識別。
圖10 共炮點(diǎn)道集記錄(拾取反射軸較明顯位置作為Q值層析反演輸入)
圖11 改進(jìn)前(a)、后(b)衰減旅行時Q值層析法計算得到的Q值剖面
圖12 以改進(jìn)前(a)、后(b)衰減旅行時層析法計算的Q值為輸入進(jìn)行Q補(bǔ)償偏移后的道集對比
對CRP道集做疊加得到偏移疊加剖面如圖13所示。從頻譜上看,改進(jìn)后的衰減旅行時層析方法的頻帶更寬,主頻由30Hz提高到40Hz,分辨率得到提高。從剖面上看,對于一些由于吸收衰減作用能量較弱、模糊不清的層位,經(jīng)吸收衰減補(bǔ)償后得以正確顯示,地下介質(zhì)結(jié)構(gòu)刻畫更加清晰。
圖13 改進(jìn)前(a)、后(b)衰減旅行時層析法計算的Q值進(jìn)行Q偏移補(bǔ)償?shù)寞B加剖面
本文采用改進(jìn)的衰減旅行時Q層析方法計算地下縱、橫向Q值分布。衰減旅行時層析法的關(guān)鍵在于準(zhǔn)確求取衰減旅行時,衰減旅行時上的較小誤差都會導(dǎo)致最終Q值反演的較大偏差。計算模型衰減旅行時,射線追蹤的精度是關(guān)鍵,在追蹤過程中考慮速度梯度隨傳播方向的變化,能夠更準(zhǔn)確地追蹤射線路徑,獲得更準(zhǔn)確的傳播時間。從實際資料中拾取衰減旅行時,最優(yōu)頻帶的選擇通常會引入人為誤差,對數(shù)譜比反演方法不僅無需人為選擇參與反演頻率成分,反演結(jié)果也更穩(wěn)定。而且,通過在反演過程中考慮反射波在不同界面內(nèi)的旅行時之差,可以有效解決上覆地層的影響。因此得到Q值反演結(jié)果更加符合地下介質(zhì)的吸收衰減分布。
采用改進(jìn)后的方法在數(shù)值模型和實際資料應(yīng)用中取得了較好的Q值估計結(jié)果。經(jīng)Q偏移補(bǔ)償后,提高了地震資料分辨率,恢復(fù)了衰減造成的能量衰減和相位畸變,為后續(xù)AVO分析和儲層預(yù)測提供更可靠資料。
[1]DASGUPTA R,CLARK R A.Estimation ofQfrom surface seismic reflection data[J].Geophysics,1998,63(6):2120-2128
[2]陳文爽,管路平,李振春,等.基于廣義S變換的疊前Q值反演方法研究[J].石油物探,2014,53(6):706-712
CHEN W S,GUAN L P,LI Z C,et al.PrestackQ-inversion based on generalized S transform[J].Geophysical Prospecting for Petroleum,2014,53(6):706-712
[3]王小杰,印興耀.基于零相位子波地層Q值估計[J].地球物理學(xué)進(jìn)展,2011,26(6):2090-2098
WANG X J,YIN X Y.Estimation of layer quality factors based on zero-phase wavelet[J].Progress in Geophysics,2011,26(6):2090-2098
[4]王德利,戴建芳.基于射線路徑的疊前高精度Q值估計方法[J].石油物探,2013,52(5):475-481
WANG D L,DAI J F.Prestack high accuracyQestimation based on ray path[J].Geophysical Prospecting for Petroleum,2013,52(5):475-481
[5]張繁昌,張汛汛,張立強(qiáng),等.基于自適應(yīng)子波分解的品質(zhì)因子Q提取方法[J].石油物探,2016,55(1):41-48
ZHANG F C,ZHANG X X,ZHANG L Q,et al.Extraction method ofQestimation based on adaptive wavelet decompesation[J].Geophysical Prospecting for Petroleum,2016,55(1):41-48
[6]李偉娜,云美厚,黨鵬飛,等.基于微測井資料的雙線性回歸穩(wěn)定Q估計[J].石油物探,2017,56(4):483-490
LI W N,YUN M H,DANG P F,et al.StabilityQestimation by dual linear regression based on uphole survey data[J].Geophysical Prospecting for Petroleum,2017,56(4):483-490
[7]QUAN Y,HARRIST J M.Seismic attenuation tomography using the frequency shift method[J].Geophysics,1997,62(3):895-905
[8]ZHANG C,ULRYCH T J.Estimation of quality factors from CMP records[J].Geophysics,2002,67(10):1542-1547
[9]王宗俊.基于譜模擬的質(zhì)心法品質(zhì)因子估算[J].石油物探,2015,54(3):267-273
WANG Z J.Quality factor estimation by centroid frequency shift of spectrum fitting[J].Geophysical Prospecting for Petroleum,2015,54(3):267-273
[10]張立彬,王華忠,馬在田.基于積分中值參變量法的質(zhì)心頻移Q值估算[J].石油物探,2010,49(3):214-223
ZHANG L B,WANG H Z,MA Z T.Quality factor estimation by centroid frequency shift of integration median[J].Geophysical Prospecting for Petroleum,2010,49(3):214-223
[11]高靜懷,楊森林.利用零偏移VSP 資料估計介質(zhì)品質(zhì)因子方法研究[J].地球物理學(xué)報,2007,50(4):1198-1209
GAO J H,YANG S L.On the method of guaiity factors estimation from zero-offset VSP data[J].Chinese Journal of Geophysics,2007,50(4):1198-1209
[12]PARRA J O,XU P C,HACKERT C L.A borehole-model-derived algorithm for estimatingQPlogs from full-waveform sonic logs[J].Geophysics,2007,72(4):E107-E117
[13]LIU Y,WEI X.IntervalQinversion from CMP gathers:Part I absorption equation[J].Expanded Abstracts of 75thAnnual Internat SEG Mtg,2005:1717-1720
[14]LIU Y,WEI X.IntervalQinversion from CMP gathers:Part II absorption equation[J].Expanded Abstracts of 75thAnnual Internat SEG Mtg,2005:1721-1724
[15]BRZOSTOWSKI M A,MCMECHAN G A.3-D tomographic imaging of near-surface seismic velocity and attenuation[J].Geophysics,1992,57(3):396-403
[16]LIAO Q,MCMECHAN G A.Tomographic imaging of velocity andQ,with application to crosswell seismic data from the Gypsy Pilot Site,Oklahoma[J].Geophysics,1997,62(6):1804-1811
[17]嚴(yán)又生,宜明理,魏新,等.井間地震速度和Q值聯(lián)合層析成像及應(yīng)用[J].石油地球物理勘探,2001,36(1):9-17
YAN Y S,YI M L,WEI X,et al.Joint tomographic imaging for cross-hole seismic velocity andQvalue[J].Oil Geophysical Prospecting,2001,36(1):9-17
[18]NOWACK R L,MICHAEL P M.Inversion of seismic attributes for velocity and attenuation structure[J].Geophysics of Journal International,1997,128(3):689-700
[19]樊計昌,劉明軍,趙成彬,等.岫巖隕石坑三維Q值層析成像[J].地球物理學(xué)報,2010,53(10):2367-2375
FAN J C,LIU M J,ZHAO C B,et al.Three-dimensionalQtomography for Xiuyuan meteorite impact center[J].Chinese Journal of Geophysics,2010,53(10):2367-2375
[20]CAVALCA M,MOORE I,ZHANG L.Ray-based tomography forQestimation andQcompensation in complex media[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:3989-3992
[21]XIN K F,HE Y,XIE Y.RobustQtomographic inversion through adaptive extraction of spectral features[J].Expanded Abstracts of 84thAnnual Internat SEG Mtg,2014:3726-3730
[22]JIN Z Q.An improvedQtomography inversion and its application[J].Expanded Abstracts of 85thAnnual Internat SEG Mtg,2015:1721-1724
[23]魏文,王小杰,李紅梅.基于疊前道集小波域Q值求取方法研究[J].石油物探,2011,50(4):356-360
WEI W,WANG X J,LI H M.Qestimation in wavelet domain using prestack data[J].Geophysical Prospecting for Petroleum,2011,50(4):356-360
[24]王小杰,印興耀,吳國忱.基于變換的吸收衰減技術(shù)在含氣儲層預(yù)測中的應(yīng)用研究[J].石油物探,2012,51(1):38-42
WANG X J,YIN X Y,WU G S.Application on reservoir prediction using attenuation technique based on transform[J].Geophysical Prospecting for Petroleum,2012,51(1):38-42
[25]TONN R.The determination of seismic quality factorQfrom VSP data:a comparison of different computational techniques[J].Geophysical Prospecting,1991,45(1):1-27