基于VMD的高分辨率頻率衰減分析算法

薛雅娟，陳輝，劉哲哿，楊佳，鐘劍丹，周娟，李東芳

1.氣象信息與信號(hào)處理四川省高校重點(diǎn)實(shí)驗(yàn)室(成都信息工程大學(xué))，四川 成都 610225

2.數(shù)學(xué)地質(zhì)四川省重點(diǎn)實(shí)驗(yàn)室(成都理工大學(xué))，四川 成都 610059

地震波的能量衰減是地震波在地下介質(zhì)中傳播期間發(fā)生的能量損耗，是介質(zhì)的一種內(nèi)在屬性。除了巖層的不均勻性和幾何擴(kuò)散等造成的彈性損耗會(huì)導(dǎo)致地震波傳播過程中的能量衰減以外，巖石的黏彈性性質(zhì)也會(huì)導(dǎo)致地震波能量衰減，使地震波的彈性能轉(zhuǎn)化為熱能然后耗散掉，造成地震波的總能量衰減[1-6]。由地震波彈性勢(shì)能轉(zhuǎn)化為熱能形成的地震波本征衰減與各種介質(zhì)和流體特性(如流體含量，滲透率，黏度)有關(guān)[7]，使其成為潛在的儲(chǔ)層監(jiān)測(cè)工具[8-9]。在相對(duì)穩(wěn)定的垂直和水平巖性幾乎沒有變化的地層結(jié)構(gòu)中，地震波的衰減主要是由流體特性引起的[10-13]。實(shí)驗(yàn)室和實(shí)際數(shù)據(jù)測(cè)量均表明，在大多數(shù)頻率帶寬上黏性流體飽和巖石的地震波衰減比干燥巖石的地震波衰減更為明顯，且地震波的高頻能量衰減大于低頻能量衰減[10-14]。

能量吸收分析方法[15]是一種通過計(jì)算吸收系數(shù)來估計(jì)介質(zhì)高頻能量吸收的有效的衰減估計(jì)方法，已被廣泛用于油氣探測(cè)中[13-21]。能量吸收分析方法的關(guān)鍵是準(zhǔn)確計(jì)算吸收系數(shù)或者是衰減梯度。常規(guī)能量吸收分析方法使用兩點(diǎn)斜率法計(jì)算吸收系數(shù)，僅適用于旁瓣很少或幾乎沒有旁瓣的平滑頻譜或具有高信噪比(SNR)的地震數(shù)據(jù)[17]。近年來為了提高衰減估計(jì)的準(zhǔn)確性，具有更高時(shí)頻分辨率的時(shí)頻分析方法，如小波變換等被引入能量吸收分析方法中，同時(shí)，為了提高吸收系數(shù)計(jì)算的準(zhǔn)確性，最小二乘法也被引入替代傳統(tǒng)的二點(diǎn)斜率法[15-16]，為避免不同頻率地震波衰減不同造成的相互影響而采用自適應(yīng)經(jīng)驗(yàn)?zāi)B(tài)分解及其衍生方法[13,17-19]、量子信息處理算法[20]的高精度衰減估計(jì)方法也引入到能量吸收分析中。

對(duì)于地震波反射數(shù)據(jù)，能量吸收分析方法通常是通過地震波的譜信息來計(jì)算。然而，目前依賴于地震波譜信息分析的能量吸收分析方法中很少考慮反射系數(shù)的影響。實(shí)際上，在許多沉積地層中，地震波反射譜的主要成分來自相互干擾的近距離反射，即眾所周知的調(diào)諧效果[22]，地震數(shù)據(jù)譜的形狀(包括幅度、頻率和相位)不可避免地受到調(diào)諧效應(yīng)的影響[23-25]。為了抑制反射系數(shù)的影響，校正地震波頻譜形狀，給出更為準(zhǔn)確可靠的衰減估計(jì)結(jié)果，本文引入變分模態(tài)分解(VMD)方法，結(jié)合平均對(duì)數(shù)譜距離(ALSD)發(fā)展一種自適應(yīng)魯棒的能量吸收分析方法。

1 基于VMD的抑制反射系數(shù)影響的頻率衰減分析算法

能量吸收分析方法主要利用指數(shù)衰減函數(shù)exp(-af)在時(shí)頻域中估計(jì)高頻能量吸收的數(shù)量[14,17]：

X(t，f)=ke(-af)

(1)

式中：X(t，f)為地震波x(t)的時(shí)頻分布；k為常數(shù)；f為頻率；a為吸收系數(shù)或稱作衰減梯度，被廣泛應(yīng)用于烴類指示。對(duì)式(1)兩邊取以e為底的對(duì)數(shù)，得：

lnX(t，f)=lnk-af

(2)

目前衰減梯度估計(jì)中，通常對(duì)式(2)直接利用兩點(diǎn)斜率法或者最小二乘法估計(jì)衰減梯度。這里，為了抑制反射系數(shù)影響，校正地震波頻譜，首先考慮地震波的平穩(wěn)卷積模型：

x(t)=w(t)*r(t)

(3)

式中：w(t)為地震子波函數(shù)；r(t)為反射系數(shù)序列。

式(3)在頻域中變?yōu)椋?/p>

X(f)=W(f)R(f)

(4)

對(duì)式(4)兩邊取對(duì)數(shù)得：

lnX(f)=lnW(f)+lnR(f)

(5)

通常情況下，反射系數(shù)序列r(t)具有類似隨機(jī)噪聲特性[26]，相對(duì)于子波函數(shù)，它的頻譜變化非常快，而子波函數(shù)相對(duì)于反射系數(shù)序列通常變化很慢，通常是一個(gè)低頻信號(hào)。地震信號(hào)在對(duì)數(shù)域可以看作是低頻子波函數(shù)序列和相對(duì)高頻反射系數(shù)序列的和。

接下來，引入VMD方法進(jìn)行分析。由于VMD分解是在時(shí)頻域中進(jìn)行處理的，因此首先利用具有較高分辨率的廣義S變換將地震信號(hào)轉(zhuǎn)換到時(shí)頻域。此處VMD方法可以看作是一個(gè)自適應(yīng)、多頻帶濾波器，VMD是一個(gè)約束變分問題，表示為[27]：

(6)

式中：u表示非線性和非平穩(wěn)信號(hào)；ck表示第k個(gè)本征模態(tài)函數(shù)(IMF)；ωk是角頻率。

ck在ωk周圍幾乎是緊湊的。通常，使用二次懲罰和拉格朗日乘數(shù)來解決問題式(6)。 通過采用乘數(shù)的交替方向方法(ADMM)，可以獲得每個(gè)篩選過程中的每個(gè)模式和中心頻率。VMD方法的更多細(xì)節(jié)及主要步驟詳見文獻(xiàn)[27]。經(jīng)過VMD分解后，原始信號(hào)u被分解為一系列窄帶本征模態(tài)函數(shù)ck的和：

(7)

式中：n為IMF的個(gè)數(shù)。

對(duì)于對(duì)數(shù)域地震信號(hào)lnX(f)，經(jīng)過VMD分解后，可以得到一系列從低頻到高頻的IMF。從而，可以找到一個(gè)分離點(diǎn)，在分離點(diǎn)之前的信號(hào)主要反映低頻的子波函數(shù)序列，分離點(diǎn)后的信號(hào)主要反映高頻的反射系數(shù)序列。為此，引入ALSD作為有用的相似性度量以有效地區(qū)分子波函數(shù)主導(dǎo)的IMF和反射系數(shù)序列主導(dǎo)的IMF。ALSD是兩個(gè)對(duì)數(shù)頻譜之間的距離度量，被廣泛應(yīng)用于語音編碼領(lǐng)域中用于測(cè)量量化或內(nèi)插性能[28-30]。通過最大化ALSD來確定子波函數(shù)主導(dǎo)的IMF和反射系數(shù)序列l(wèi)nR(ω)主導(dǎo)的IMF之間的分離點(diǎn)：

(8)

其中，對(duì)于兩個(gè)IMF之間的ALSD定義為：

(9)

式中：m是ck的長度。然后，t1時(shí)刻地震子波函數(shù)W(ω)主導(dǎo)的子信號(hào)可以通過下式重構(gòu)：

(10)

提取到地震子波函數(shù)主導(dǎo)的子信號(hào)后，結(jié)合式(2)進(jìn)一步可以獲得衰減梯度估計(jì)式：

(11)

式中：f1為地震譜峰值頻率?？紤]到處理后的窄頻譜特性，f2按下式進(jìn)行動(dòng)態(tài)選?。?/p>

(12)

式中：fm為頻譜主要范圍截止點(diǎn)。

為了給出更為準(zhǔn)確的衰減梯度估計(jì)值，此處利用最小二乘法擬合獲得衰減梯度a。設(shè)擬合多項(xiàng)式φ(f)=a0+a1f，則最小二乘法是對(duì)于給定的數(shù)據(jù)y=s(f)的數(shù)據(jù)點(diǎn)(fi，lnW(t，fi)；i=1，2，…，m)求近似曲線φ(f)，且使φ(f)與原曲線y=s(f)偏差最小，可以表示為：

(13)

式中：δi為近似曲線在點(diǎn)(fi，ln(W(t，fi))處的偏差。

對(duì)式(13)右邊求對(duì)于ai的偏導(dǎo)數(shù)，列出系數(shù)矩陣求解，確定擬合系數(shù)a1。則該時(shí)刻的衰減梯度a=|a1|?；赩MD的抑制反射系數(shù)影響的能量吸收分析算法的處理流程如圖1所示。

圖1 基于VMD的抑制反射系數(shù)影響的頻率衰減分析算法的流程圖

2 實(shí)際地震數(shù)據(jù)應(yīng)用

利用四川盆地川西坳陷某氣田的疊后偏移地震數(shù)據(jù)進(jìn)行分析。該數(shù)據(jù)前期處理中經(jīng)過了地表一致性振幅補(bǔ)償提高振幅保真性的處理，但不包括自適應(yīng)反褶積或時(shí)變頻譜白化、幾何擴(kuò)散補(bǔ)償、自動(dòng)增益AGC等幅度補(bǔ)償處理。該氣田研究區(qū)河道砂發(fā)育，儲(chǔ)層厚度較薄，物性較差且致密，平均孔隙度8.66%，平均滲透率0.21 mD。受非均質(zhì)性的影響，不同砂組及砂體含氣性差異較大，氣水分布十分復(fù)雜，在構(gòu)造高、低部位均有天然氣及地層水分布。該研究區(qū)存在的一類儲(chǔ)層從地球物理特征角度來看具有中等波阻抗、弱地震響應(yīng)的特征，筆者認(rèn)為其屬于一類隱蔽儲(chǔ)層，在疊后數(shù)據(jù)處理中傳統(tǒng)方法較難識(shí)別。

2.1 過井道基于VMD的抑制反射系數(shù)影響的衰減特征及算法分析

表1 IMF分量的ALSD

圖2 含氣水井A的測(cè)井解釋結(jié)果

圖3 過井道衰減分析

2.2 二維地震剖面處理

圖4為研究區(qū)一條過含氣井B的疊后偏移地震剖面。含氣井B為斜井，儲(chǔ)層位置如圖4中橢圓區(qū)域所示，目標(biāo)層地震響應(yīng)較弱，橢圓區(qū)域地震信號(hào)與周圍地震信號(hào)差異不明顯，含氣性檢測(cè)困難，屬于隱蔽儲(chǔ)層。含氣井B的測(cè)井解釋結(jié)果如圖5所示。該剖面基于VMD的抑制反射系數(shù)影響的衰減梯度剖面如圖6(a)所示，可以看出，橢圓區(qū)域有強(qiáng)振幅異常響應(yīng)，由于該段內(nèi)巖性穩(wěn)定，排除巖性、地層等因素影響，認(rèn)為該強(qiáng)振幅異常特征表明該橢圓區(qū)域存在烴類，符合測(cè)井解釋結(jié)果的含氣性檢測(cè)結(jié)果。而在基于小波變換的衰減梯度剖面(見圖6(b))中，在相同橢圓區(qū)域位置并沒有強(qiáng)振幅異常響應(yīng)，說明該算法不能很好地給出符合測(cè)井解釋結(jié)果的含氣性檢測(cè)結(jié)果，表明針對(duì)隱蔽儲(chǔ)層的含氣性檢測(cè)問題，傳統(tǒng)方法無法有效檢測(cè)到微弱地震響應(yīng)信號(hào)的強(qiáng)振幅異常特征。

圖4 過含氣井B的疊后偏移地震剖面

圖5 含氣井B的測(cè)井解釋結(jié)果

圖6 基于VMD抑制反射系數(shù)影響的衰減梯度剖面與基于小波變換的衰減梯度剖面對(duì)比

3 討論

3.1 VMD分解次數(shù)選取

基于VMD的抑制反射系數(shù)影響的頻率衰減分析方法中，重構(gòu)抑制反射系數(shù)影響的頻譜是該算法的一個(gè)關(guān)鍵步驟。合適的VMD分解次數(shù)影響著重構(gòu)抑制反射系數(shù)影響的頻譜的準(zhǔn)確性。盡管通常情況下，選擇IMF的數(shù)量(即分解級(jí)數(shù))是VMD的關(guān)鍵問題，VMD分解次數(shù)設(shè)置過大，雖然可以盡量地避免子波信息和反射系數(shù)信息的重疊，但是會(huì)導(dǎo)致分解產(chǎn)生的IMF分量過分解，捕獲額外的噪聲或者導(dǎo)致模態(tài)重復(fù)；而如果VMD分解次數(shù)設(shè)置過小，會(huì)導(dǎo)致數(shù)據(jù)分割不足，導(dǎo)致子波信息和反射系數(shù)信息重疊較多，產(chǎn)生模態(tài)混疊。目前雖然關(guān)于VMD分解次數(shù)選取有很多研究并發(fā)展了很多方法[31-34]，但在本次研究的算法中選擇并不困難。在本次研究的算法中，主要通過計(jì)算各個(gè)IMF分量的ALSD來確定子波函數(shù)主導(dǎo)的IMF和反射系數(shù)序列主導(dǎo)的IMF之間的分離點(diǎn)，而不是獲得每個(gè)IMF的精確表達(dá)，因此，由過多模式引起的模式重復(fù)的傳統(tǒng)難題在本次研究中不需要考慮，只需要避免生成太少的模式導(dǎo)致的模態(tài)混疊問題，即子波信息和反射系數(shù)信息重疊太多的情況。這種情況很容易解決，對(duì)于實(shí)際地震數(shù)據(jù)，通常使用過井道結(jié)合測(cè)井解釋結(jié)果來獲得適當(dāng)?shù)腣MD分解級(jí)數(shù)。此外，需要說明的是，盡管本次研究使用VMD算法來分離子波信息和反射系數(shù)信息，但是實(shí)際上子波信息和反射系數(shù)信息在對(duì)數(shù)譜中并不是完全分開的，由于篩選出的信息主要是子波信息，且反射系數(shù)信息相對(duì)子波信息頻率更高，體現(xiàn)在頻譜上主要是高頻細(xì)節(jié)波動(dòng)信息，而子波信息決定了總體頻譜曲線的趨勢(shì)，個(gè)別細(xì)節(jié)微小波動(dòng)信息不會(huì)影響到總體頻譜曲線趨勢(shì)，從而對(duì)衰減梯度計(jì)算結(jié)果的影響很小。圖6(a)的結(jié)果也表明了該方法可以給出準(zhǔn)確度很高的衰減梯度估計(jì)值。

此外，由于VMD的強(qiáng)局域分解能力和抗噪聲性能[27]，噪聲信息往往會(huì)被提取為單獨(dú)的IMF信號(hào)，因此本文的方法也可以適用于低信噪比地震數(shù)據(jù)的處理。

3.2 衰減梯度估計(jì)的高頻優(yōu)勢(shì)頻段選取

基于VMD的抑制反射系數(shù)影響的頻率衰減分析方法中，影響結(jié)果精度的另一個(gè)因素是高頻優(yōu)勢(shì)頻段的選取。傳統(tǒng)衰減梯度估計(jì)算法中，通常選取總能量65%到總能量85%的頻段進(jìn)行衰減梯度計(jì)算。本次研究中，針對(duì)采用了廣義S變換后時(shí)頻譜分辨率較高、頻帶較窄，且經(jīng)過IMF選擇重構(gòu)后的抑制反射系數(shù)頻譜頻帶又變窄的特性，采用了動(dòng)態(tài)頻段選取方式，頻段選取方式如式(12)所示。該式中采取了高頻優(yōu)勢(shì)頻段最大頻段計(jì)算范圍為40 Hz的處理方式。在實(shí)際地震數(shù)據(jù)處理中，也可以根據(jù)實(shí)際地震數(shù)據(jù)特征，適當(dāng)調(diào)整該參數(shù)。

3.3 含氣性檢測(cè)結(jié)果分析

含烴類區(qū)域通常在衰減梯度剖面中表現(xiàn)為強(qiáng)振幅異常特征，排除巖性、地層等其他因素影響，通常認(rèn)為衰減梯度剖面中的強(qiáng)振幅異常特征可以給出含烴類信息的解釋結(jié)果。本次研究中針對(duì)四川盆地川西坳陷某氣田的疊后偏移地震數(shù)據(jù)進(jìn)行分析，選取了隱蔽儲(chǔ)層段進(jìn)行研究。從過含氣水井A的地震道分析及過含氣井B的地震剖面分析中可以看到，常規(guī)基于小波變換的衰減梯度計(jì)算方法檢測(cè)不到隱蔽儲(chǔ)層中的烴類信息，不能給出很好的烴類解釋結(jié)果，而基于VMD的抑制反射系數(shù)影響的衰減分析方法可以檢測(cè)到這類弱含氣響應(yīng)信號(hào)的強(qiáng)振幅異常信息，給出符合測(cè)井解釋結(jié)果的含烴類檢測(cè)結(jié)果。

4 結(jié)束語

本文提出了一種基于VMD和ALSD的抑制反射系數(shù)影響的頻率衰減分析算法，該算法避免了常規(guī)頻率衰減分析算法中反射系數(shù)的影響，能實(shí)現(xiàn)更準(zhǔn)確地估計(jì)頻率衰減信息。為了抑制反射系數(shù)影響，該算法在時(shí)頻域?qū)Φ卣鸬赖膶?duì)數(shù)譜進(jìn)行VMD分解，結(jié)合ALSD判斷子波函數(shù)主導(dǎo)的IMF和反射系數(shù)序列主導(dǎo)的IMF之間的分離點(diǎn)，通過重構(gòu)獲得抑制了反射系數(shù)影響的頻譜。另一方面，為了保證提取的頻率衰減信息的準(zhǔn)確性，對(duì)重構(gòu)的抑制了反射系數(shù)影響的頻譜采用了動(dòng)態(tài)選取頻段的方式計(jì)算頻率衰減信息，實(shí)現(xiàn)窄頻帶情況下高頻衰減信息計(jì)算頻段的優(yōu)勢(shì)選取。四川盆地川西坳陷某氣田的疊后偏移地震數(shù)據(jù)處理結(jié)果表明，該算法可以有效識(shí)別出一些儲(chǔ)層特征不明顯的隱蔽儲(chǔ)層，檢測(cè)出弱烴類信號(hào)引起的能量異常，給出符合測(cè)井解釋結(jié)果的烴類解釋。