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

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

1.氣象信息與信號處理四川省高校重點實驗室(成都信息工程大學)，四川 成都 610225

2.數學地質四川省重點實驗室(成都理工大學)，四川 成都 610059

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

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

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

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

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

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

(1)

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

lnX(t，f)=lnk-af

(2)

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

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

(3)

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

式(3)在頻域中變為：

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

(4)

對式(4)兩邊取對數得：

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

(5)

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

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

(6)

式中：u表示非線性和非平穩信號；ck表示第k個本征模態函數(IMF)；ωk是角頻率。

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

(7)

式中：n為IMF的個數。

對于對數域地震信號lnX(f)，經過VMD分解后，可以得到一系列從低頻到高頻的IMF。從而，可以找到一個分離點，在分離點之前的信號主要反映低頻的子波函數序列，分離點后的信號主要反映高頻的反射系數序列。為此，引入ALSD作為有用的相似性度量以有效地區分子波函數主導的IMF和反射系數序列主導的IMF。ALSD是兩個對數頻譜之間的距離度量，被廣泛應用于語音編碼領域中用于測量量化或內插性能[28-30]。通過最大化ALSD來確定子波函數主導的IMF和反射系數序列lnR(ω)主導的IMF之間的分離點：

(8)

其中，對于兩個IMF之間的ALSD定義為：

(9)

式中：m是ck的長度。然后，t1時刻地震子波函數W(ω)主導的子信號可以通過下式重構：

(10)

提取到地震子波函數主導的子信號后，結合式(2)進一步可以獲得衰減梯度估計式：

(11)

式中：f1為地震譜峰值頻率。考慮到處理后的窄頻譜特性，f2按下式進行動態選取：

(12)

式中：fm為頻譜主要范圍截止點。

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

(13)

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

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

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

2 實際地震數據應用

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

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

表1 IMF分量的ALSD

圖2 含氣水井A的測井解釋結果

圖3 過井道衰減分析

2.2 二維地震剖面處理

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

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

圖5 含氣井B的測井解釋結果

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

3 討論

3.1 VMD分解次數選取

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

此外，由于VMD的強局域分解能力和抗噪聲性能[27]，噪聲信息往往會被提取為單獨的IMF信號，因此本文的方法也可以適用于低信噪比地震數據的處理。

3.2 衰減梯度估計的高頻優勢頻段選取

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

3.3 含氣性檢測結果分析

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

4 結束語

本文提出了一種基于VMD和ALSD的抑制反射系數影響的頻率衰減分析算法，該算法避免了常規頻率衰減分析算法中反射系數的影響，能實現更準確地估計頻率衰減信息。為了抑制反射系數影響，該算法在時頻域對地震道的對數譜進行VMD分解，結合ALSD判斷子波函數主導的IMF和反射系數序列主導的IMF之間的分離點，通過重構獲得抑制了反射系數影響的頻譜。另一方面，為了保證提取的頻率衰減信息的準確性，對重構的抑制了反射系數影響的頻譜采用了動態選取頻段的方式計算頻率衰減信息，實現窄頻帶情況下高頻衰減信息計算頻段的優勢選取。四川盆地川西坳陷某氣田的疊后偏移地震數據處理結果表明，該算法可以有效識別出一些儲層特征不明顯的隱蔽儲層，檢測出弱烴類信號引起的能量異常，給出符合測井解釋結果的烴類解釋。