張繁昌,張汛汛,張立強(qiáng),宗兆云

(中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院,山東青島266580)

?

基于自適應(yīng)子波分解的品質(zhì)因子Q提取方法

張繁昌,張汛汛,張立強(qiáng),宗兆云

(中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院,山東青島266580)

摘要:品質(zhì)因子Q和速度一樣是地層的固有參數(shù),其大小決定地震波吸收衰減的強(qiáng)弱。由于Q值的大小與地震子波的傳播時間和頻率的乘積有關(guān),因此,提出了一種基于自適應(yīng)子波分解的品質(zhì)因子提取方法。首先,利用自適應(yīng)子波分解算法對地震信號進(jìn)行子波分解,得到的匹配子波可認(rèn)為是不同時刻的時變子波;然后,將不同匹配子波的中心時間和中心頻率的乘積與相應(yīng)匹配子波的對數(shù)振幅投影到平面坐標(biāo)中,得到隨該乘積值變化的對數(shù)振幅散點分布;最后通過擬合這些離散點斜率提取Q值。為壓制離群值對擬合精度的影響,采用L1模線性回歸方法對離散點進(jìn)行線性擬合。理論模型測試和實際資料應(yīng)用的結(jié)果表明,該方法能比較準(zhǔn)確地提取品質(zhì)因子,有效補償?shù)卣鸩芰?提高了地震資料的分辨率。

關(guān)鍵詞：品質(zhì)因子;自適應(yīng)子波分解;匹配子波;L1模線性回歸;分辨率

地震波在地下介質(zhì)中傳播時會發(fā)生能量衰減,其中部分能量不可逆地轉(zhuǎn)化為其它形式能量,且高頻能量的耗散快于低頻,導(dǎo)致地震波的振幅衰減和波形畸變,降低地震資料的分辨率[1-2]。地層吸收衰減特性通常用品質(zhì)因子Q來度量[3]。FUTTERMAN[4]在1962年第一次提出品質(zhì)因子方面的理論及相關(guān)估計方法;隨后,KJARTANSSON[5]提出了通過研究脈沖寬度在時間域提取Q值的上升時間法;WARD等[6]提出了通過計算空間上兩個不同點的振幅衰減提取品質(zhì)因子的振幅衰減法;HAUGE等[7]提出了利用不同時刻的頻譜比值反演Q值的頻率域Q值計算方法;QUAN等[8]提出了質(zhì)心頻移法,該方法利用頻譜質(zhì)心頻率的偏移研究地層吸收衰減特性;GAO等[9]推導(dǎo)了頻譜峰值頻率與品質(zhì)因子的關(guān)系,并提出了峰值頻率法。宮同舉等[10]對品質(zhì)因子提取方法進(jìn)行了比較,指出沒有哪種方法適用于任何情況,且各種方法對地震資料品質(zhì)的要求都很高;針對譜比法和質(zhì)心頻移法的優(yōu)缺點,趙寧和曹思遠(yuǎn)等先后提出頻域統(tǒng)計屬性組合方法[11]和對數(shù)譜根式方法[12]提取Q值,提高了Q值的準(zhǔn)確性;王小杰等[13]將小波分頻技術(shù)運用到品質(zhì)因子提取中,通過分頻處理避免不同頻率成分地震信號的干擾,以提高求取Q的穩(wěn)定性。

自適應(yīng)子波分解與重構(gòu)算法是MALLAT等[14]提出的一種較好的信號表示方法,但龐大的計算量一直制約該方法的應(yīng)用;LIU等[15]提出了動態(tài)匹配追蹤算法,該方法以瞬時特征為先驗信息,提高了算法計算效率;為減小原子庫規(guī)模,張繁昌等[16-17]在LIU的算法基礎(chǔ)上提出了雙參數(shù)快速匹配追蹤算法,在復(fù)數(shù)域中將時頻原子控制參數(shù)減少為2個;張汛汛等[18]用相位展開法得到穩(wěn)定的瞬時頻率,通過掃描得到最優(yōu)匹配頻率,降低了計算量。分析品質(zhì)因子計算方法以及衰減機(jī)理后我們發(fā)現(xiàn),Q值的計算與地震子波的傳播時間和頻率的乘積有關(guān)。鑒于此,提出了一種基于自適應(yīng)子波分解的品質(zhì)因子Q提取方法。首先,基于自適應(yīng)子波分解算法將地震信號分解為具有不同中心時間t0，中心頻率f0以及振幅A的一系列匹配子波,將不同匹配子波t0和f0的乘積(t0·f0)作為橫坐標(biāo),相應(yīng)的振幅取對數(shù)作為縱坐標(biāo)投影到坐標(biāo)平面內(nèi),得到隨t0·f0變化的對數(shù)振幅離散點;然后,鑒于L1模線性回歸方法能較好地壓制離群值的影響[19],采用該方法對離散點進(jìn)行線性擬合,得到最優(yōu)斜率;最后利用斜率與品質(zhì)因子的關(guān)系提取Q值。

1方法原理

1.1自適應(yīng)子波分解

自適應(yīng)子波分解方法的核心思想是通過創(chuàng)建超完備匹配子波庫對信號進(jìn)行自適應(yīng)分解,最后信號可表示為一系列匹配子波的線性疊加。一個超完備的匹配子波庫是一系列匹配子波的集合,是信號自適應(yīng)子波分解的基礎(chǔ)。常用的匹配子波通過對基函數(shù)g(t)進(jìn)行伸縮、平移以及調(diào)制得到,并且滿足條件：g(t)是實的連續(xù)可微函數(shù),g(t)∈O[1/(t2+1)],‖g(t)‖=1,∫g(t)dt≠0且g(0)≠0。

令γ=(α,u,ζ),則匹配子波定義為：

(1)

設(shè)D={gγ(t)}γ∈Γ是由一系列匹配子波組成的超完備子波庫,從子波庫中搜索最優(yōu)匹配子波gγ0,使其滿足：

(2)

(2)式中,h是待分解信號,可以表示為：

(3)

式中：r1h是信號一次分解后的殘差項;“<·>”表示內(nèi)積。整個自適應(yīng)子波分解過程需要迭代進(jìn)行,直到殘差項滿足所設(shè)定的閾值。最后,經(jīng)過自適應(yīng)子波分解后,信號h(t)可表示為一系列最優(yōu)匹配子波和殘差項的線性疊加：

(4)

式中：Ai,ti,fi,φi分別表示各個最優(yōu)匹配子波的振幅、中心時間、中心頻率以及相位。經(jīng)過多次迭代后,殘差項的值已達(dá)到設(shè)定的閾值,則信號h(t)為：

(5)

圖1是用衰減合成地震記錄進(jìn)行自適應(yīng)子波分解的分解過程與重構(gòu)效果對比圖。從圖1a可以看出,經(jīng)過3次迭代分解后,衰減記錄的各子波被完全分解出來,分解得到的各個匹配子波可以認(rèn)為是不同時刻的衰減記錄時變子波;圖1b是重構(gòu)結(jié)果對比圖,可見自適應(yīng)子波分解算法能實現(xiàn)信號的完全重構(gòu)。

1.2衰減模型

當(dāng)一個平面波在粘彈性介質(zhì)中傳播時,在時間域可表示為：

(6)

其中,f是頻率,k是復(fù)波數(shù)[20]：

(7)

式中：α是衰減因子;v表示相速度。在KOLSKY衰減模型[21]中,相速度定義為：

(8)

衰減因子表示為：

(9)

式中：Q是品質(zhì)因子;v0和f0分別表示參考相速度和參考頻率。將(8)式和(9)式代入(7)式可得：

(10)

在頻率域求解波動方程平面波,可得頻域內(nèi)衰減信號H(f)與原始信號H0(f)的關(guān)系為：

(11)

式中：δ表示傳播距離。將(10)式代入(11)式得到：

(12)

從式(12)看出,當(dāng)僅考慮振幅衰減時,可通過振幅信息和傳播時間與頻率的乘積之間的關(guān)系提取Q值。

1.3L1模線性回歸提取Q值

對信號進(jìn)行自適應(yīng)子波分解后,得到一系列具有不同中心時間t0，中心頻率f0和振幅A的匹配子波。將中心時間和中心頻率的乘積t0·f0,以及相應(yīng)匹配子波的對數(shù)振幅lnA投影到坐標(biāo)平面內(nèi),得到隨t0·f0變化的對數(shù)振幅散點分布。研究發(fā)現(xiàn),對數(shù)振幅和t0·f0有以下關(guān)系：

(13)

式中：A0是未衰減振幅;χ=t0·f0。這樣Q值就可以通過公式(13)的斜率得到。為減小離群值對Q值提取的影響,采用基于L1模的線性回歸方法求取斜率。

將各離散點用(13)式表示,并變換為矩陣形式：

(14)

式中：a和b分別表示擬合曲線的斜率和截距。為求取斜率,在L1模下需滿足：

(15)

其中,

r是殘值振幅向量。對(15)式求導(dǎo)：

(16)

即E(m)的梯度為：

(17)

其中,R是對角矩陣,其對角元素為：

(18)

由(17)式可得L1模的解[22],即：

m=(GTRG)-1GTRd

(19)

這樣Q值就可以通過最優(yōu)斜率a提取：

(20)

2方法測試與應(yīng)用效果分析

2.1理論模型測試

圖2說明了Q值的求取過程,圖2a表示常Q衰減信號。首先利用自適應(yīng)子波分解算法對該信號進(jìn)行子波分解,將不同匹配子波中心時間和中心頻率的乘積和相應(yīng)對數(shù)振幅投影到坐標(biāo)平面,如圖2b 中離散點所示,然后通過L1模回歸方法擬合最優(yōu)直線(圖2b中的紅線),最后采用(20)式計算Q值。

圖3a和圖3b分別表示理論Q值為50和80的合成衰減地震記錄,可以看出,隨著旅行時的增加地震記錄衰減逐漸加大,并且Q值越小,衰減越明顯。為補償?shù)卣鹩涗浀乃p,利用本文方法進(jìn)行Q值提取,最后計算得到圖3a和圖3b的Q值分別為49.1和79.3,與模型值的誤差僅為0.9%和1.7%,證明本文方法能夠比較精確地提取Q值。分別利用所求Q值對圖3a和圖3b進(jìn)行補償,結(jié)果見圖3c和圖3d,與無衰減記錄(圖3e)對比發(fā)現(xiàn),補償后衰減振幅得到有效恢復(fù)。圖4是Q值為80衰減地震記錄補償前、后時頻譜對比圖,可以看出,經(jīng)過補償后中深部衰減能量得到有效恢復(fù),特別是高頻成分的能量。

圖5是對Q值為80的合成衰減地震記錄(圖3b)進(jìn)行處理后的分時間段0～0.7s,0～1.4s,0～2.0s離散振幅點圖以及擬合得到平均品質(zhì)因子Qa;然后,利用Qa計算得到不同時間段地層品質(zhì)因子Qn。表1為不同時間段的Qa和Qn,從中可以看出,計算得到的不同時間段的地層品質(zhì)因子具有較好的準(zhǔn)確度。

表1　不同時間段平均品質(zhì)因子Qa及地層品質(zhì)因子Qn比較

圖2　Q值的提取過程a 常Q衰減信號; b 利用對數(shù)線性擬合曲線求Q

圖3　合成地震記錄衰減補償a Q=50衰減記錄; b Q=80衰減記錄; c Q=50補償記錄; d Q=80補償記錄; e 無衰減記錄

圖4　Q值為80衰減模型補償前(a)和后(b)時頻譜對比

圖5　不同時間段對數(shù)線性擬合求平均品質(zhì)因子Qa

2.2實際資料應(yīng)用

我們從一個160道、采樣間隔為2ms的實際地震剖面中取一道地震數(shù)據(jù),然后分別采用譜比法和自適應(yīng)子波分解方法求取Q值，其結(jié)果如圖6所示。從圖6可以看出,兩條曲線趨勢一致,但基于自適應(yīng)子波分解方法提取的Q值比譜比法得到的Q值更穩(wěn)定。圖7為截取一段地震剖面(0.8～1.4s)進(jìn)行地震衰減補償處理前后的剖面,圖7a是補償前地震剖面,圖中深層能量被吸收衰減,分辨率受到極大影響。圖7b是利用本文方法對每道提取Q值后的補償后剖面,與圖7a 對比發(fā)現(xiàn),深層地震波能量得到了補償,波形具有更好的一致性,同相軸連續(xù)性增強(qiáng),弱反射得以恢復(fù),復(fù)合波得到分離,波場信息更加豐富。

圖8為圖7數(shù)據(jù)對應(yīng)的頻譜,可見原始地震剖面的地震波能量主要集中在中低頻,高頻能量較弱,因此時間域剖面的分辨率較低;經(jīng)過補償后的地震記錄頻譜,高頻能量得到增強(qiáng),頻帶拓寬,因而時間域地震記錄的分辨率得以提高。

圖6　提取的Q值曲線對比

圖7　實際地震資料吸收衰減補償前后結(jié)果a 補償前; b 補償后

圖8　補償前后振幅譜對比

3結(jié)束語

地震波吸收衰減補償是增強(qiáng)深層能量、提高分辨率的重要途徑之一。本文基于自適應(yīng)子波分解算法,將地震信號分解為不同的匹配子波,通過不同匹配子波中心時間和中心頻率的乘積與相應(yīng)對數(shù)振幅的線性關(guān)系提取Q值。本方法通過降維處理,將二維信息時間和頻率降低為一維信息時間與頻率的乘積,增加了算法的穩(wěn)健性;同時,在線性擬合過程中運用基于L1模的回歸方法,減小離群值對擬合精度的影響,確保Q值的求取精度。理論模型測試發(fā)現(xiàn),本方法能夠精確求取Q值,并且具有良好的穩(wěn)健性;實際地震資料應(yīng)用表明,經(jīng)過提取Q值的吸收衰減補償后,地震資料衰減能量得到有效恢復(fù),分辨率有了明顯提高。

參考文獻(xiàn)

[1]魏文,王小杰,李紅梅.基于疊前道集小波域Q值求取方法研究[J].石油物探,2011,50(4),355-360

WEI W,WANG X J,LI H M.Study on extraction method for Q based on pre-stack gather in wavelet domain[J].Geophysical Prospecting for Petroleum,2011,50(4),355-360

[2]王小杰,印興耀,吳國忱.粘彈性介質(zhì)地震波傳播特征及反射特征研究[J].物探化探計算技術(shù),2012,34(3):258-266

WANG X J,YIN X Y,WU G C.Characteristics of seismic wave propagation and reflection in visoelastic medium[J].Physical and Geochemical Exploration Technology,2012,34(3),258-266

[3]AKI K,RICHARDS P G.Quantitative seismology (second edition)[M].Sausalito:University Science Books,2002:1-700

[4]FUTTERMAN W I.Dispersive body waves[J].Journal of Geophysical Research,1962,67(18):5279-5291

[5]KJARTANSSON E.Constant Q wave propagation and attenuation[J].Journal of Geophysical Research,1979,84(B9):4737-4748

[6]WARD R E.Mapping seismic attenuation within geothermal systems using teleseisms with application to the Geysers2 Clear Lake region[J].Journal of Geophysical Research,1980,85(B8):5227-5236

[7]HAUGE P S.Measurements of attenuation from vertical seismic profiles[J].Geophysics,1981,46(10):1548-1558

[8]QUAN Y,HARRIS J M.Seismic attenuation tomography using the frequency shift method[J].Geophysics,1997,62(3):895-905

[9]高靜懷,楊森林.利用零偏移VSP資料估計介質(zhì)品質(zhì)因子方法研究[J].地球物理學(xué)報,2007,50(4):1198-1209

GAO J H,YANG S L.On the method of quality factors estimation from zero-offset VSP data[J].Chinese Journal of Geophysics,2007,50(4):1198-1209

[10]宮同舉,孫成禹,彭洪超,等.幾種提取品質(zhì)因子方法的對比分析[J].勘探地球物理進(jìn)展,2009,32(4):252-256

GONG T J,SUN C Y,PENG H C,et al.Comparison of several computational methods of quality factor[J].Progress in Exploration Geophysics,2009,32(4):252-256

[11]趙寧,曹思遠(yuǎn),王宗俊,等.頻域統(tǒng)計性屬性組合提取品質(zhì)因子Q[J].石油地球物理勘探,2013,48(4):545-552

ZHAO N,CAO S Y,WANG Z J,et al.Seismic Q estimation by combination of frequency statistics attributes[J].Oil Geophysical Petroleum,2013,48(4):545-552

[12]曹思遠(yuǎn),譚佳,高明,等.對數(shù)譜根式法Q值反演[J].石油地球物理勘探,2014,49(1):161-166

CAO S Y,TAN J,GAO M,et al.Seismic Q estimation with logarithmic spectrum equation root[J].Oil Geophysical Petroleum,2014,49(1):161-166

[13]王小杰,欒錫武.基于小波分頻技術(shù)的地層Q值提取方法研究[J].石油物探,2015,54(3):260-266

WANG X J,LUAN X W.Q value extraction method based on wavelet frequency division technology[J].Geophysical Prospecting for Petroleum,2015,54(3):260-266

[14]MALLAT S,ZHANG Z.Matching pursuit with time-frequency dictionaries[J].IEEE Transactions on Signal Processing,1993,41(12):3397-3415

[15]LIU J,WU Y,HAN D,et al.Time-frequency decomposition based on Ricker wavelet [J].Expanded Abstracts of 74thAnnual Internat SEG Mtg,2004:1937-1940

[16]張繁昌,李傳輝.基于正交時頻原子的地震信號快速匹配追蹤[J].地球物理學(xué)報,2012,55(1):277-283

ZHANG F C,LI C H.Seismic signal fast matching pursuit based on the orthogonal time-frequency atom[J].Chinese Journal of Geophysics,2012,55(1):277-283

[17]張繁昌,李傳輝.匹配追蹤保幅地震AVF剖面及同頻率剖面的構(gòu)建[J].石油物探,2013,52(3):234-239

ZHANG F C,LI C H.The construction of the preserved amplitude seismic AVF section and the same frequency section based on the matching pursuit[J].Geophysical Prospecting for Petroleum,2013,52(3):234-239

[18]張汛汛,張繁昌,劉漢卿.基于快速匹配追蹤算法的地震道集剩余時差校正[J].石油物探,2015,54(4):420-426

ZHANG X X,ZHANG F C,LIU H Q.Seismic gathers residual moveout correction based on fast matching pursuit algorithm[J].Geophysical Prospecting for Petroleum,2015,54(4):420-426

[19]ASTER R C,BORCHERS B,THURBER C H.Parameter estimation and inverse problems[M].New York:Elsevier Academic Press,2005,1-480

[20]WANG Y H.Seismic inverse Q filtering [M].New York:Wiley-Black well,2008:1-248

[21]KOLSKY H.The propagation of stress pulses in viscoelastic solids[J].Philosophical Magazine,1956,1(8):693-710

[22]ZHANG F C,DAI R H,LIU H Q.Seismic inversion based on L1-norm misfit function and total variation regularization[J].Journal of Applied Geophysics,2014,109(1):111-118

(編輯：朱文杰)

Extraction method for quality factorQbased on adaptive wavelet decomposition

ZHANG Fanchang,ZHANG Xunxun,ZHANG Liqiang,ZONG Zhaoyun

(SchoolofGeosciences,ChinaUniversityofPetroleum,Qingdao266580,China)

Abstract:Quality factor Q is the inherent parameter of strata,which determines the strength of seismic attenuation.Since the Q value is related with the product of travel time and frequency of seismic wavelet,we proposed a method to extract the value of Q based on the adaptive wavelet decomposition.Firstly,we decompose seismic signal using the adaptive wavelet decomposition algorithm.The decomposed pursuit wavelets can be considered as time-varying seismic wavelets at different time.Then,the products of central time and central frequency versus the corresponding logarithmic amplitudes of the pursuit wavelets are projected into a coordinate plane,to obtain the scatter distribution of logarithmic amplitudes with these products.Finally,the slope of these scatters is fitted to extract Q value.In order to decrease the influence of outliers,L1 norm linear regression method is used in the linear fitting.Theoretical and real data tests show this method can accurately extract quality factor,compensate seismic wave energy and improve the resolution of seismic data.

Keywords:quality factor,adaptive wavelet decomposition,pursuit wavelet,L1 norm linear regression,resolution

文章編號：1000-1441(2016)01-0041-08

DOI:10.3969/j.issn.1000-1441.2016.01.006

中圖分類號：P631

文獻(xiàn)標(biāo)識碼：A

基金項目：國家重點基礎(chǔ)研究發(fā)展計劃(973計劃)項目(2013CB228604)、國家科技重大專項(2016ZX0500804)、山東省中青年科學(xué)家獎勵基金(2014BSE28009)和中國石油大學(xué)(華東)研究生創(chuàng)新工程(YCX2015006)聯(lián)合資助。

作者簡介：張繁昌(1972—),男,教授,博士,主要從事地震儲層預(yù)測方法研究。通訊作者：張汛汛(1991—),男,碩士,主要從事地球物理探測方法與技術(shù)研究。

收稿日期：2015-06-19;改回日期：2015-09-22。

張繁昌,張汛汛,張立強(qiáng),等.基于自適應(yīng)子波分解的品質(zhì)因子Q提取方法[J].石油物探,2016,55(1):-48

ZHANG Fanchang,ZHANG Xunxun,ZHANG Liqiang,et al.Extraction method for quality factorQbased on adaptive wavelet decomposition[J].Geophysical Prospecting for Petroleum,2016,55(1):-48

This research is financially supported by the National Key Basic Research and Development Program of China (973 Program) (Grant No.2013CB228604),the National Science and Technology Major Project of China (Grant No.2016ZX0500804),Young and Middle-Aged Scientists Research Awards Foundation of Shangdong Province (Grant No.2014BSE28009) and Postgraduate Innovation Project of China University of Petroleum (Grant No.YCX2015006).