脈沖響應字典架構(gòu)的薄層譜分解方法

李雪英，王福霖，萬喬升

(1.東北石油大學 地球科學學院， 黑龍江 大慶 163318； 2.黑龍江省油氣藏形成機理與資源評價重點實驗室， 黑龍江 大慶 163318)

0 引 言

譜分解，也稱時頻分析或譜估計，是一種基本的信號分析工具，其目的在于揭示地震記錄的時頻特征[1]。譜分解通過將一維時域信號變換到二維的時頻譜中，描述頻率隨時間的變化規(guī)律。短時傅里葉變換、連續(xù)小波變換、S變換和魏格納威利分布等傳統(tǒng)的譜分解方法均受到海森堡測不準原理的限制[2]。

Chen等[3]提出基追蹤算法，通過稀疏約束反演，在譜分解中表現(xiàn)出優(yōu)異的精度和穩(wěn)定性。Portniaguine等[4]提出將譜分解作為一種反演問題去解決，比較了不同反演約束條件下的時頻譜，得出稀疏約束具有最高時頻分辨率的結(jié)論。目前，基于基追蹤的譜分解已廣泛應用于微地震數(shù)據(jù)處理[5]、衰減估計[6]和薄層分析等領(lǐng)域[7]，在地震勘探中，薄層及薄互層的研究一直是一個難點問題。薄層是指厚度小于入射子波λ/8的層狀介質(zhì)，考慮噪音、地震波衰減等因素，薄層被重新定義為小于λ/4的層狀介質(zhì)[8]。薄層是地層中最簡單的一種地質(zhì)結(jié)構(gòu)，但卻是研究薄互層地震特征響應的重要基礎(chǔ)，深入了解薄層和薄互層的地震響應特征，尤其是時頻域響應特征，對提高薄層分辨率、開展薄層預測具有重要意義。

利用基追蹤進行譜分解的過程中，尤其研究薄層方面，其字典原子多采用Ricker或Mallet子波，通過褶積模型進行實驗[9-13]，但褶積模型所能表征的最小薄層厚度為一個采樣間隔，尚不能與地層厚度精確匹配，同時該模型忽略地震波在傳播過程中的頻散與衰減效應，不能完全模擬地震波在地下介質(zhì)中真實的傳播情況[14-18]。因此，以波動方程為基礎(chǔ)，基于波場延拓開展薄層、薄互層正演模擬，研究薄層反射復合波可能更具實際意義[8]。

筆者基于波場延拓理論，利用深度域相移法對不同厚度單砂層和等厚薄互層進行正演模擬，研究各個主頻Ricker子波激發(fā)的層界面脈沖響應，截取其有效部分，歸一化處理后作為原子，構(gòu)建脈沖響應字典，利用基追蹤算法對理論合成的薄層反射復合波進行譜分解，驗證該方法對薄層時頻分辨率的改善情況。

1 基追蹤原理

2004年，Portniaguine等[4]提出地震記錄可以表示成一系列小波和相應的系數(shù)褶積之和：

(1)

式中：Nb——小波的數(shù)量；

t——時間；

n——小波的尺度，該參數(shù)控制著小波頻率。

式(1)改寫為矩陣：

(2)

式中：ψN——主頻的卷積矩陣；

Φ——字典，其每一列被稱為原子；

η——隨機噪聲；

α——字典對應的系數(shù)序列向量，其中包含了地震記錄s(t)的時頻分布信息。

將上述α系數(shù)序列向量按頻率重列成t×n的矩陣，即為時頻譜。求解最稀疏系數(shù)序列α本質(zhì)上是尋找地震信號最稀疏表達，該問題可以表示為

(3)

式中，‖α‖0——系數(shù)序列向量中非零元素的數(shù)目。

由于字典Φ是過完備的，因此使用式(3)求稀疏分解的結(jié)果是NP-Hard問題。Chen等提出基追蹤算法[3]來解決信號稀疏表示的問題，基追蹤基于全局優(yōu)化算法來實現(xiàn)稀疏表示的問題，其求解模型為

(4)

用l1范數(shù)代替l0范數(shù)的微妙差異改變了問題的性質(zhì)，將NP難問題轉(zhuǎn)化成基追蹤問題，而基追蹤問題可以變形為線性規(guī)劃問題從而求解?；粉櫵惴ǖ闹饕獎?chuàng)新于基于線性規(guī)劃理論，通過最小化l1范數(shù)求解系數(shù)向量，將稀疏分解轉(zhuǎn)化為約束極值問題。 Donoho進行的實驗證明，基追蹤算法不僅可以使用盡可能少的原子實現(xiàn)原始信號的表示，而且可以獲得比匹配追蹤更穩(wěn)定的結(jié)果。

將基追蹤算法轉(zhuǎn)化為線性規(guī)劃問題即可得到稀疏約束下的精確解，進而求得系數(shù)序列α。此外，投影梯度法、對偶對數(shù)障礙規(guī)劃算法等都是有效的求解方法。除算法外，預定義的字典也很大程度上影響算法求解精度[2]。

采用波動方程正演模擬的單界面脈沖響應信號，制備了全新的脈沖響應字典，以式(4)為稀疏反演框架，采用線性規(guī)劃算法獲取薄層、薄互層的系數(shù)序列矩陣。

2 脈沖響應字典的制備與地質(zhì)模型

2.1 脈沖響應字典制備

基于脈沖響應字典架構(gòu)的譜分解方法是用人為構(gòu)建的脈沖響應字典實現(xiàn)地震信號的譜分解，其核心在于脈沖響應字典的構(gòu)建。三種不同子波的時域波形如圖1所示。

圖1 三種不同子波的時域波形 Fig. 1 Time-domain waveforms of three different wavelets

建立單反射界面模型，界面上下的厚度均為3 000 m，其中砂巖的速度為2 918 m/s，密度為2.14 g/cm3；泥巖的速度為3 180 m/s，密度為2.32 g/cm3(全文模型參數(shù)與此一致)。

基于波場延拓理論，應用深度相移法正演模擬獲取單界面地震響應。采用零相位雷克子波作為激發(fā)震源，激發(fā)主頻以1 Hz為間隔，從1 Hz變化到100 Hz，得到100個不同激發(fā)主頻的地震反射記錄。

此時激發(fā)得到的地震記錄不是復合波，而是單一界面的脈沖響應，提取中間道的脈沖信號，截取信號的有效部分(圖1a)進行歸一化處理，以消除信號振幅大小對譜分解結(jié)果的影響。將該信號作為脈沖響應字典中的一系列原子，即式(1)中的ψ(t,n)。將其沿著方陣主對角線方向平移構(gòu)建卷積矩陣ψN。將不同頻率的卷積矩陣ψN按頻率變化依次排列，組成的矩陣即構(gòu)成了脈沖響應字典。

ψ1={(ψ(t,1))(ψ(t-1,1))…(ψ(1,1))},ψ2={(ψ(t,2))(ψ(t-1,2))…(ψ(1,2))},…,ψN={(ψ(t,N))(ψ(t-1,N))…(ψ(1,N))},

(5)

式中：()——由脈沖信號構(gòu)成的列向量；

{}——不同頻率的卷積矩陣；

ψ(t-τ,N)——時移量為τ、主頻為N的脈沖信號。

需要說明的是，因為脈沖信號都經(jīng)過歸一化處理，所以界面反射系數(shù)的大小不影響譜分解結(jié)果。文中根據(jù)實驗需要，構(gòu)建了頻率范圍足夠?qū)挼拿}沖響應字典，在應用中，應根據(jù)實際需要設(shè)計合適的頻率范圍。

2.2 地質(zhì)模型的建立

建立單砂層(雙界面)模型，該模型由均勻泥巖背景中包含單砂層組成，單砂層厚度以1 m為間隔，從1 m遞增到30 m，以考察不同厚度下，脈沖響應字典的時頻譜反演精度。薄層上下泥巖厚度各為3 000 m，使單砂層位于地質(zhì)模型的中心位置。使用了正演模擬方法和觀測系統(tǒng)設(shè)置方式[19]，共模擬出30套數(shù)據(jù)體。

建立等厚薄互層模型，該模型由均勻泥巖背景中包含互層數(shù)分別為2、3的砂泥巖等厚地質(zhì)體組成，用以考察不同互層數(shù)對脈沖響應字典下時頻譜反演效果的影響。模型參數(shù)與觀測系統(tǒng)設(shè)置方式與上文保持一致。單層厚度以1 m為間隔，從1 m變化至30 m，共模擬出60套數(shù)據(jù)體。

3 仿真實驗

單砂層結(jié)構(gòu)簡單，是構(gòu)成薄互層的基本單元之一，其時頻譜簡單明了，有助于脈沖響應字典與復Ricker子波字典、實Ricker子波字典時頻分析效果的對比。同時，分析單砂層的時頻特征，總結(jié)規(guī)律，對薄互層時頻特征的分析有指導作用。

等厚薄互層較單砂層結(jié)構(gòu)復雜，時頻譜也相對復雜，有助于判斷上述字典在互層構(gòu)造下分解的精度與穩(wěn)定性，并進行對比。

文中分別用上述字典對單砂層和等厚薄互層正演模擬信號進行譜分解，對比時頻譜，分析脈沖響應字典的優(yōu)勢、劃定時頻譜頻率和位置精度并分析各自的時頻特征。

3.1 單砂層模型的時頻特征

基于波動方程正演模擬的單砂層合成地震記錄時域波形如圖2所示。當厚度在1到9 m(λ/8)之間時，波形振幅隨厚度的增加線性增大；當厚度在9到19 m(λ/4，地層厚度整數(shù)值與對波長的等分不完全對應)之間時，來自頂、底的兩個反射波發(fā)生相消干涉，復合反射波的振幅隨地層厚度減薄而減小，成非線性變化；當厚度在19到37 m(λ/2，圖中未畫出)時，自頂、底的兩個反射波發(fā)生相長干涉，隨厚度增大，振幅略有減小[25]。

圖2 不同厚度的單砂層時域波形 Fig. 2 Time-domain waveforms of single sand layer with different thicknesses

圖3為使用復Ricker子波字典對圖2時域波形進行譜分解的時頻譜，小圖內(nèi)左上方的標識代表層厚度。無論厚度如何變化，頻譜均為雙峰。受薄層升頻特性影響，隨著厚度的不斷增大，雙峰的峰值頻率從36 Hz附近下降到30 Hz附近，同時隨著相互干涉的兩個反射波逐漸分離，基追蹤算法求得的解越來越稀疏，使得其頻率范圍從32～42 Hz縮短到26～32 Hz，頻帶寬度從10 Hz減小到6 Hz。

當薄層厚度小于λ/4，即厚度小于19 m時，反射波相互干涉程度加劇，波峰與波谷間的頻率大于波形主頻，能量團雙峰向中間高頻處聚攏，呈“八”字型分布；當厚度等于λ/4時，能量團近于平行，當厚度大于λ/4，即厚度大于19 m時，干涉程度減弱，波峰與波谷間的頻率小于波形主頻，能量團雙峰向中間低頻處聚攏，呈倒“八”字型分布。

圖3 復Ricker子波字典分解的單砂層時頻譜Fig. 3 Time-frequency spectrum of single sand layer decomposed by complex Ricker wavelet dictionary

圖4為使用實Ricker子波字典對圖2時域波形進行譜分解的時頻譜。與圖3相比，時頻譜時間分辨率更高；受薄層升頻特性影響，隨著厚度的不斷增大，雙峰的峰值頻率從36 Hz附近下降到30 Hz附近；大部分能量團的頻帶寬度穩(wěn)定在6 Hz左右。

當薄層厚度小于λ/4，即厚度小于19 m時，能量團呈“八”字型分布；當厚度等于λ/4時，能量團近于平行，當厚度大于λ/4，即厚度大于19 m時，能量團呈倒“八”字型分布。其形成原因與圖3相同，但能量的分布形態(tài)相較于圖3不明顯。

圖5為使用脈沖響應字典對圖2時域波形進行譜分解的時頻譜。與圖4相比，該圖的能量團更為聚焦，出現(xiàn)的位置更為精確，分解出的頻率范圍更小。無論厚度如何變化，頻譜均為雙峰(干擾能量團除外)，當厚度在1～19 m之間時，雙峰的間距無明顯變化，當厚度大于19 m時，雙峰的間距隨時域復合波的分離而逐漸增加。受薄層升頻特性影響，隨厚度的不斷增大，雙峰的峰值頻率從34 Hz附近下降到30 Hz附近；單個能量團的頻帶寬度在2 Hz左右，頻率聚焦性比圖3和4時頻譜提升3倍。由此說明，相比于使用復Ricker子波和實Ricker字典，使用脈沖響應字典對薄層進行譜分解能更精確地識別復合波頻率。

當厚度在1～14 m(3/16λ)之間時，時頻譜中出現(xiàn)一些干擾能量團，這是薄層厚度過小，復合波時域波形和原子的時域波形不匹配，基追蹤算法難以求解出稀疏解所致；當厚度在14～19 m之間時，干擾能量團逐漸消失，時域波形得到了較為精確地時頻分解；當厚度大于19 m時，干擾能量團消失，時域波形得以精確分解。能量團發(fā)育位置反映時域波形波峰、波谷位置，對應的頻率反映時域波形的峰值頻率。

當厚度在1～18 m之間時，能量團成正“八”字型分布，當厚度在19 m時，能量團近于平行，當厚度在20～24 m之間時能量團逐漸成倒“八”字型分布。該現(xiàn)象的出現(xiàn)與圖3原因相同。此外，時頻譜中發(fā)育等頻的雙能量團可局部指示出單砂層結(jié)構(gòu)，這與單砂層是薄互層單元之一的概念是統(tǒng)一的，且當單層厚度大于19m 時，首尾能量團可精確指示單砂層頂?shù)捉缑嫖恢谩?/p>

圖4 實際Ricker子波字典分解的單砂層時頻譜Fig. 4 Time-frequency spectrum of single sand layer decomposed by actual Ricker wavelet dictionary

圖5 脈沖響應字典分解的單砂層時頻譜Fig. 5 Time-frequency spectrum of single sand layer decomposed by impulse response dictionary

3.2 等厚薄互層模型的時頻特征

基于波動方程正演模擬的薄互層合成地震記錄時域波形如圖6所示。

圖6 不同互層數(shù)、不同厚度的等厚薄互層時域波形Fig. 6 Time-domain waveforms of equal thickness and thin interlayers with different numbers of interlayers and different thicknesses

由圖6可見，等厚薄互層時域波形振幅的變化規(guī)律與圖3相同。

圖7和圖8為使用復Ricker子波字典對圖6時域波形進行譜分解的時頻譜。圖9和10為使用實Ricker子波字典對圖6時域波形進行譜分解的時頻譜。

圖7 復Ricker子波字典分解的互層數(shù)為2的等厚薄互層時頻譜Fig. 7 Time-frequency spectrum of equal thickness and thin interlayer with 2 interlayers decomposed by complex Ricker wavelet dictionary

由圖7、8可見，當互層數(shù)為2，厚度在1～9 m之間時，時頻譜表現(xiàn)為雙峰；厚度在9～14 m之間時，時頻譜峰數(shù)由二向四過渡；厚度大于或等于14 m時，時頻譜表現(xiàn)為四峰。當互層數(shù)為3，厚度在1～5 m之間時，時頻譜表現(xiàn)為雙峰；厚度在5～14 m之間時，時頻譜峰數(shù)由二向六過渡；厚度大于或等于14 m時，時頻譜表現(xiàn)為六峰。

圖8 復Ricker子波字典分解的互層數(shù)為3的等厚薄互層時頻譜Fig. 8 Time-frequency spectrum of equal thickness and thin interlayer with 3 interlayers decomposed by complex Ricker wavelet dictionary

總體上看，受薄層升頻特性影響，隨著厚度的不斷增大，時頻譜峰值頻率由37 Hz下降到23 Hz，單個能量團的頻帶寬度大部分在6 Hz左右。

當薄層厚度小于λ/4，厚度小于19 m時，反射波相互干涉程度加劇，波形中間位置的波峰與波谷頻率較兩側(cè)更高，能量團向高頻處移動，呈“n”型分布(5 m例外)；當厚度大于λ/4，即厚度大于19 m時，干擾能量團出現(xiàn)，反射波相互干涉程度減弱，波形中間位置的波峰與波谷頻率較兩側(cè)更低，能量團向低頻處移動，呈“u”型分布。

圖9 實Ricker子波字典分解的互層數(shù)為2的等厚薄互層時頻譜Fig. 9 Time-frequency spectrum of equal thickness and thin interlayer with 2 interlayers decomposed by real Ricker wavelet dictionary

薄層厚度小于λ/4，厚度小于19 m，能量團呈“n”型分布(5 m例外)；即厚度大于19 m時，能量團呈“u”型分布。其形成原因與圖7和8相同。

圖10 實Ricker子波字典分解的互層數(shù)為3的等厚薄互層時頻譜Fig. 10 Time-frequency spectrum of equal thickness and thin interlayer with 3 interlayers decomposed by real Ricker wavelet dictionary

圖9、10與圖7、8相比，時頻譜時間分辨率更高。時頻譜峰數(shù)變化規(guī)律與圖7和8相同?？傮w上看，受薄層升頻特性影響，隨著厚度的不斷增大，時頻譜峰值頻率由35 Hz下降到23 Hz，單個能量團的頻帶寬度大部分在6 Hz左右。

圖11和12為使用脈沖響應字典對圖6時域波形進行譜分解的時頻譜。

圖11 脈沖響應字典分解的互層數(shù)為2的等厚薄互層時頻譜Fig. 11 Time-frequency spectrum of equal thickness and thin interlayer with 2 interlayers decomposed by impulse response dictionary

圖11、12與圖9、10相比，能量團更為聚焦，觀測到的波形頻率更為精確。通過觀察可見，當互層數(shù)為2，厚度在1～9 m之間時，時頻譜表現(xiàn)為雙峰；厚度在9～14 m之間時，時頻譜峰數(shù)由二向四過渡；厚度大于或等于14 m時，時頻譜表現(xiàn)為四峰。當互層數(shù)為3，厚度在1～14 m之間時，時頻譜表現(xiàn)為四峰，并逐漸向六峰過渡；厚度大于或等于14 m時，時頻譜表現(xiàn)為六峰。受薄層升頻特性影響，隨厚度的不斷增大，峰值頻率從33 Hz下降到23 Hz，單個能量團的頻帶寬度在2 Hz左右，頻率聚焦性比圖7和8時頻譜提高3倍。由此說明，相比于使用復Ricker子波和實Ricker子波字典，使用脈沖響應字典對薄互層進行譜分解能更精確地識別復合波頻率。

圖12 脈沖響應字典分解的互層數(shù)為3的等厚薄互層時頻譜Fig. 12 Time-frequency spectrum of equal thickness and thin interlayer with 3 interlayers decomposed by impulse response dictionary

當厚度在1～14 m(3/8λ)之間時，時頻譜中出現(xiàn)一些干擾能量團，原因與圖5相同；當厚度在14～19 m之間時，干擾能量團逐漸消失，時域波形得到了較為精確地時頻分解；當厚度大于19 m時，干擾能量團消失，時域波形得以精確分解。能量團發(fā)育位置反映時域波形波峰、波谷位置，對應的頻率反映時域波形峰頻。

當厚度小于19 m時，中間兩個能量團向高頻處移動，時頻譜呈“n”型分布；當厚度大于19 m時，中間兩個能量團向高頻處移動，時頻譜成“u”型分布。其形成原因與圖7和8相同。

4 結(jié) 論

(1)提出的基于脈沖響應字典架構(gòu)的基追蹤譜分解方法相比于傳統(tǒng)的Ricker子波字典譜分解方法具有相同的時域分辨率，且頻率聚焦性提升3倍，有利于薄砂層及等厚薄互層時頻特性分析。

(2)隨著單砂層厚度減薄，時頻譜峰值頻率逐漸升高，高頻能量對薄層分辨能力增強，說明薄的單砂層具有升頻特性。

(3)等厚薄互層時頻譜能量團由“n”型分布向“u”型過渡，隨著單層厚度的減薄，峰值頻率逐漸升高，峰值頻率數(shù)量及間距可清楚展現(xiàn)等厚薄互層內(nèi)部結(jié)構(gòu)特征。