基于三參數小波的頻譜分解方法

朱振宇 高佳倫 姜秀娣 孫文博 薛東川 王清振

(中海油研究總院有限責任公司,北京100028； ②中國石油大學(北京),北京 102200)

1 引言

近年來，隨著油氣勘探工作的不斷深入，頻譜分解方法已成為油氣儲層預測中的常用手段。Widess[1]利用楔形模型揭示了薄層反射與薄層厚度的關系。隨后，Neidell等[2]引入了調諧厚度的概念，即當層厚為地震子波波長的四分之一時，合成地震記錄的振幅最大，此時的地層厚度即為調諧厚度，這是利用地震反射振幅解釋薄層厚度的基礎。Partyka[3]發現，薄層反射在短時窗傅里葉變換頻譜上會出現陷頻現象，利用這種現象可以估算薄層厚度；同時，根據地震薄層反射的調諧理論，不同的調諧頻率對應不同的薄層厚度[4]，因此通過頻譜分解獲得的單頻體可以反映該頻率成分對應的地質現象[5]，挖掘更豐富的信息，還可以突出關鍵的地質目標，從而提高地震資料的解釋精度[6]。

利用頻譜分解技術進行地震成像時，需要選擇合適的頻譜分解方法[7]。Morlet提出的小波變換可以很好地刻畫非平穩信號的局部特征，長期以來一直是頻譜分解的有力工具[8]。小波變換的效果取決于基小波的選擇，由于Morlet小波具有較高的時間、頻率分辨率[9]，一直被廣泛使用。但是當Morlet小波的中心頻率較小時，其修正項就不可忽略。為此，Harrop等[10]提出了一種改進的Morlet小波，該小波在中心頻率較小時仍然可以滿足小波的允許條件。在Harrop等的研究基礎上，高靜懷等[11]提出一種靈活性高、具有三個可調參數的三參數小波，主要研究了σ(調制頻率)和τ(能量衰減因子)兩個參數。

本文研究了三參數小波的每個參數及其參數組合，重點分析了參數β(能量延遲因子)的影響。通過改變參數組合可以獲得多種形態的基小波，從而滿足不同的處理需求。首先闡述了頻譜分解方法的基本原理；然后介紹了三參數小波，重點研究了每個參數對小波的影響；隨后設計了正演模型，利用三參數小波進行時頻分析；最后利用基于三參數小波的頻譜分解方法預測儲層，取得了很好的效果。

圖1 楔形地質模型及其地震響應特征(a)楔形地質模型； (b)合成地震記錄；(c)薄層濾波器的振幅譜；(d)調諧曲線楔形地質模型的蓋層和底層具有相同的波阻抗

2 頻譜分解方法

頻譜分解方法的理論基礎來源于薄層反射的調諧現象。由楔形地質模型(圖1a)的正演地震記錄(圖1b)可見，薄層的存在對入射波而言相當于一個濾波器，由薄層濾波器響應的傅里葉振幅譜可以看到周期性頻陷現象(圖1c)。提取層厚度為25、50ms對應的兩條調諧曲線(圖1d)，曲線的第一個峰值頻率即為該厚度對應的調諧頻率，曲線的陷頻頻率Pf與薄層厚度T(雙程旅行時)存在如下對應關系[3]

(1)

利用式(1)就可以計算薄層厚度。

通過上述分析(圖1)可知，不同的薄層厚度對應不同的調諧頻率，換言之，不同的頻率成分包含的地質信息也會有所差異。所以在頻譜分解時可分頻解釋[12](圖2)，在不同的單頻體上突顯關鍵目標。利用頻譜分解可以開展針對性的解釋工作，包括刻畫砂體的展布和河道的邊界等[13-16]，其工作流程如下：

(1)選定目的層段；

(2)選擇恰當的時頻分析方法進行頻譜分解，借助分頻處理手段獲得離散的單頻體；

(3)在每個單頻體上制作沿層切片，利用頻率切片預測砂體展布、刻畫河道邊界；

(4)將不同的單頻切片進行RGB混色融合，提高儲層的解釋精度(此步驟可以靈活取舍)。

圖2 分頻解釋

3 三參數小波變換

在頻譜分解時需要選擇合適的時頻分析方法。目前已發展的時頻分析方法種類繁多[17]，小波變換是其中的一種方法。20世紀80年代初期，Grossmann 等[18]發現傳統的傅里葉變換在分析地震信號的局部特征時遇到了瓶頸，無法刻畫信號在某種頻率出現的時間位置。為此，引入小波變換的概念彌補這一缺陷。定義信號s(t)的小波變換(WT)為

(2)

式中：ψ(t)為基小波，上角“*”表示取共軛；a≠0為尺度因子；b為平移因子。小波變換的特點是既考慮頻率分辨率，也考慮時間分辨率，在分析信號時具有“變焦“功能，可很好地刻畫信號的局部特征[19]。小波變換結果與基小波的選擇關系緊密[20]。在實際的應用中，Morlet小波具有良好的聯合時頻分辨率，使用較為廣泛。Morlet小波的時間域表達式為

(3)

式中ω0≥5rad/s。由式(3)可以看出，Morlet小波只有一個可調參數，因此在使用時不夠靈活。并且當ω0取較小值時，其時域局部化效果不能令人滿意。三參數小波可以解決上述問題，其表達式為

ψ(t,Γ)= e-τ(t-β)2{p(Γ)[cos(σt)-k(Γ)]+

iq(Γ)sin(σt)}

(4)

其中

(5)

(6)

(7)

Γ=(σ,τ,β)

(8)

式中：σ為小波的調制頻率；τ為能量衰減因子；β為能量延遲因子[11]。

通過改變三參數小波的3個參數，可以得到不同形態的基小波。因此三參數小波的變化形式多樣，選擇靈活[21]，可以滿足不同的處理需求。下面著重討論每個參數對小波形態的影響。分析式(4)可以看出，三參數小波是由衰減函數和復變函數相乘得到的，這個復變函數的實部和虛部均為三角函數，而衰減函數相當于一個窗口。所以三參數小波使時間軸上無限分布的函數變為緊支撐集函數，從而符合小波“小”的概念。σ作為調制頻率，控制三角函數的頻率，影響小波的震蕩程度，即σ越大，小波的震蕩越劇烈。

圖3為不同σ的三參數小波。由圖可見，隨著σ的不斷增大，小波旁瓣增多，震蕩加快。τ作為能量衰減因子，控制衰減函數的衰減速度，即τ越大，衰減越快，窗口越窄，小波的波形就越窄。圖4為不同τ的三參數小波。由圖可見，隨著τ的不斷增大，小波在時間域的展布變窄，波形變窄。值得注意的是，σ對小波波形的影響程度與τ的取值有關，當τ取值較大(≥3)時，σ對小波的影響變小。這是因為τ較大時，窗口過窄，改變σ，就不易體現小波的周期變化，即震蕩程度變化不明顯。圖5為τ較大時不同σ的三參數小波，與τ較小時不同σ的三參數小波(圖3)相比，當τ較大時，改變σ，小波震蕩程度變化不明顯。

圖3 不同σ的三參數小波(a)Γ=(1,1,0)；(b)Γ=(3,1,0)； (c)Γ=(6,1,0)

圖4 不同τ的三參數小波(a)Γ=(2.5,1,0)；(b)Γ=(2.5,2,0)； (c)Γ=(2.5,3,0)

圖5 τ較大時不同σ的三參數小波(a)Γ=(1,3,0)；(b)Γ=(3,3,0)； (c)Γ=(6,3,0)

圖6 不同β的三參數小波(a)Γ=(π,3,0)； (b)Γ=(π,3,0.6)； (c)Γ=(π,3,1)； (d)Γ=(π,3,2)

最后討論β對小波的影響，β作為能量延遲因子對小波形態的影響較復雜。如果改變β，會使衰減函數產生時間延遲，即窗口在時間軸滑動。當β取三角函數周期的整數倍時，窗口就按周期的整數倍移動，則小波只會發生時移；當β為周期的非整數倍時，窗口內被改造的函數成分發生變化，因此小波不僅發生時移，還產生變形。由不同β的三參數小波(圖6)可見：①固定σ為π(周期為2)，當β=0.6時，小波不僅發生時移，還產生變形(圖6b)，這種變形即為相位延遲?？梢岳斫鉃橛捎讦赂淖?，導致窗口平移(也可以認為窗口不變，三角函數平移β個單位)，則相位延遲σβ個單位。②當β=1.0(為周期的一半)時，小波實部與虛部都發生反轉，且產生1s的時移(圖6c)。③當β=2.0(圖6d)時，小波只產生2s的時移，波形并沒有發生改變。由此可見：三參數小波在β為周期的整數倍時，可以匹配零相位子波；在β為周期的非整數倍時，可以匹配非零相位子波。

4 模型測試及實際地震資料應用

4.1 模型測試

通過模型測試三參數小波和Morlet小波時頻分析的差異，測試模型是反射系數變化的薄互層地質模型(圖7a)。圖7為測試模型及不同小波的時頻譜。由圖可見： ①對比ω0=5rad/s的Morlet小波(圖7c)和Γ=(5,2,0)的三參數小波時頻譜(圖7d)發現，前者只有一個強能量團，看不到任何薄層信息，后者顯示高頻能量呈鋸齒狀分布，即利用三參數小波做基小波可反映部分薄層信息，而Morlet小波不能精確地反映地下的地質細節。②對比ω0=1rad/s的Morlet小波(圖7e)和Γ=(1,2,0)的三參數小波時頻譜(圖7f)發現，當ω0較小時提高了Morlet小波的時間分辨率，但這是以犧牲頻率分辨率為代價的，其頻率分辨率幾乎為零(圖7e)，也就失去了時頻域聯合解釋的意義； 在三參數小波時頻譜(圖7f)中，主頻帶仍然在50Hz附近，且能量團發生分裂，出現多個錐形譜，可識別薄互層地質結構，即三參數小波的靈活性高，能更好地刻畫復雜地質現象。③對比Γ=(5,2,0)(圖7d)和Γ=(1,2,0)的三參數小波時頻譜(圖7f)發現，后者(即當σ取較小值時)的時間分辨率更高，刻畫薄層的能力更強。

圖7 測試模型及不同小波的時頻譜(a)反射系數序列； (b)圖a數據與50Hz雷克子波褶積得到的地震記錄； (c)ω0=5rad/s的Morlet小波時頻譜；(d)Γ=(5,2,0)的三參數小波時頻譜； (e)ω0=1rad/s的Morlet小波時頻譜； (f)Γ=(1,2,0)的三參數小波時頻譜圖a中的反射系數正、負相間，其絕對值從0.1增大到0.5，再從0.5減小至0.1，共計10個反射系數，相鄰兩個反射系數的時間間隔為10ms

4.2 實際地震資料應用

4.2.1 三參數小波頻譜分解

實際資料來源于渤海某靶區，目的層砂體薄厚不均，疊置關系復雜，砂體橫向分布不穩定，砂、泥巖交替頻繁。靶區內河道眾多，多期河道縱、橫分布，橫向變化快。對于河流相儲層而言，刻畫河道的展布及邊界，弄清河道之間的疊置關系是儲層預測中非常重要的工作。

三參數小波頻譜分解的關鍵在于參數的選取，結合前文的研究以及實際地震資料，參數選取時需要注意以下問題： 第一，σ取較小值更利于刻畫地質細節。第二，當τ值較大(≥3)時，由于σ對小波形態的影響變小，所以不同σ的頻譜分解結果差異不大(圖5)。如由三參數小波頻譜分解切片(圖8)可見， 即使σ存在很大差異，處理結果也沒有太大區別。第三，前文研究表明，令β=0可匹配零相位地震子波，β≠0雖然可以更好地匹配非零相位子波，但由于β非零造成基小波時移，則頻譜分解后的數據體也會產生時移，導致對地質體埋深的錯誤解釋。

圖9為不同β的三參數小波頻譜分解切片。由圖可見：當β=0時，在目的層下方38ms的切片(圖9a)上可見河道信息；當β=5時，在目的層下方38ms的切片(圖9b)上見不到任何河道信息，而在目的層下方92ms的切片(圖9c)上可見河道信息。可見，β非零會造成對河道埋深的錯誤解釋。

圖8 三參數小波頻譜分解切片(40Hz)(a)Γ=(1,3,0)； (b)Γ=(6,3,0)

圖9 不同β的三參數小波頻譜分解切片(40Hz)(a)Γ=(2.5,0.5,0)(目的層之下38ms)； (b)Γ=(2.5,0.5,5)(目的層之下38ms)； (c)Γ=(2.5,0.5,5)(目的層之下92ms)

針對該區目的層段的實際情況提取地震子波，并與不同的三參數小波做相關，優選相關性大的參數組合作為該目的層段的小波[7]。本文選擇Γ=(1,0.5,0)進行頻譜分解，并與Morlet小波的頻譜分解效果進行對比。對于河流相儲層而言，主要從河道展布及邊界刻畫是否清晰作為評價標準。圖10為三參數小波和Morlet小波的頻譜分解結果。由圖可見：在20Hz切片上，Morlet小波的信噪比低，對河道邊界(紅色圓圈內)的刻畫不如三參數小波清晰；在80Hz切片上，在Morlet小波展示的河道(藍色圓圈內)展布不連續，邊界不清晰，而三參數小波的頻譜分解效果較好。綜上所述，在該區的河流相儲層預測中，三參數小波頻譜分解的效果優于Morlet小波。

4.2.2 RGB混色融合

利用頻譜分解預測儲層時，可采用分頻解釋技術獲得離散的單頻體，但由于某個地質體在時頻域的響應不僅僅對應一種頻率，所以每個單頻體只反映地質體的一個時間厚度信息，為了提高解釋精度，需要綜合這些單頻信息，RGB混色融合技術的基本思想即源于此[22]。這里的R代表紅色，G代表綠色，B代表藍色，將三種顏色按照不同比例疊加，可以產生不同的顏色。本文用20Hz數據體作為紅色分量，40Hz數據體作為綠色分量，80Hz數據體作為藍色分量，即R、G、B分別代表低、中、高頻。圖11為20、40和80Hz切片的RGB融合顯示。由圖可見，通過RGB融合顯示，河道的能量更強，邊界更清晰，展布更連續。

圖10 三參數小波(PWT)和Morlet小波(WT)的頻譜分解結果

圖11 20Hz、40Hz和80Hz切片的RGB融合顯示

5 結論

(1)三參數小波的靈活性高，σ影響小波的震蕩程度。τ影響小波的寬窄，當τ較大時，σ對小波的影響減弱。β的影響頗為復雜，當其為周期的整數倍時，小波只產生時移，可以匹配零相位子波；當其為周期的非整數倍時，小波不僅產生時移，而且產生相位延遲，可以匹配非零相位子波。同時需要注意，β非零造成河道埋深解釋錯誤。

(2)三參數小波變換較Morlet小波變換具有較好的時頻分辨率，能更好地刻畫薄互層內細微的地質沉積結構。三參數小波變換處理結果的信噪比高，展示的河道邊界清晰、砂體展布連續。因此基于三參數小波的頻譜分解方法可以高精度地預測儲層。

(3)實際處理中可將從目的層段提取的地震子波與三參數小波做相關，優選相關性大的參數組合，從而獲得最佳的三參數小波。