基于局部頻率約束的動態(tài)匹配追蹤強反射識別與分離方法

許 璐 吳笑荷 張明振 印興耀* 宗兆云

(①中國石油大學(華東)地球科學與技術(shù)學院,山東青島 266580； ②海洋國家實驗室海洋礦產(chǎn)資源評價與探測技術(shù)功能實驗室,山東青島 266580； ③中國石油化工股份有限公司勝利油田分公司物探研究院,山東東營 257022)

0 引言

石油和天然氣為重要的戰(zhàn)略資源，隨著國民經(jīng)濟的發(fā)展，油氣需求量逐漸增大，油氣勘探目標逐漸從構(gòu)造油氣藏轉(zhuǎn)向巖性油氣藏[1]，薄層及薄互層逐漸成為勘探工作的重點。砂體厚度小、橫向連通性差、難以確定沉積尖滅位置等因素增大了薄層和薄互層的識別難度[2]，特別是當目標儲層上、下存在高阻層時，較強的波阻抗差在地震剖面上呈強反射[3]，屏蔽了目標砂層的弱反射信息，難以確定薄儲層的尖滅線位置。地層尖滅線是控制油氣成藏的重要因素，對于油氣勘探、開發(fā)評價具有重要意義[4]。在一定的地域范圍內(nèi)，通常地層的結(jié)構(gòu)特征變化不大，在地震剖面上同相軸的橫向連續(xù)性較好[5]，其反射子波的振幅、頻率、相位等屬性保持相對一致。因此根據(jù)這種特征利用匹配追蹤稀疏理論拾取和剝離強反射干擾，進而削弱強反射對薄層砂體有效反射信息的屏蔽作用。

Mallat等[6]提出的匹配追蹤算法(Matching Pursuit)根據(jù)信號局部結(jié)構(gòu)特征構(gòu)建超完備的字典原子庫，在一組能夠適應各種時頻特性的原子上將信號展開，最終利用有限個原子稀疏表示原始信號，實現(xiàn)信號的自適應分解。Liu等[7-8]提出了動態(tài)匹配追蹤算法的概念，在每次迭代過程中參與內(nèi)積運算的時頻原子隨信號的局部瞬時屬性發(fā)生變化，提高了信號的分解效率。匹配追蹤算法在各領(lǐng)域得到廣泛應用。宋維琪等[9]將復數(shù)子波匹配追蹤算法應用于儲層預測領(lǐng)域，可較好地識別薄層砂體信息；楊昊等[10]利用匹配追蹤算法自動提取薄層頂、底反射系數(shù)，定量預測薄層厚度；張繁昌等[11]研究了基于匹配追蹤瞬時譜計算的砂體尖滅線識別方法，并用于檢測三角洲沉積砂體的尖滅線；李海山等[12]采用貪婪匹配追蹤算法分離煤層強反射；李坤等[13]利用匹配追蹤算法識別流體、刻畫儲層邊界，減少了反演流程，且反演結(jié)果與測井資料保持一致。

基于以上研究背景，本文采用動態(tài)匹配追蹤算法識別強反射，通過得到的最佳匹配強反射的子波剝離強反射?？紤]到傳統(tǒng)的瞬時頻率運算結(jié)果存在負頻率現(xiàn)象且易受噪聲影響，即當計算結(jié)果出現(xiàn)負頻率時，只能進行時頻原子頻率屬性的全局搜索，并不是真正意義上的動態(tài)搜索[14]。為此，本文引進局部頻率計算模式用于約束動態(tài)匹配追蹤算法最優(yōu)原子的搜索范圍，其頻率計算結(jié)果合理，可以確定最佳匹配強反射的地震子波。通過構(gòu)建含高阻層的地質(zhì)模型和實際資料驗證了方法的可行性。

1 匹配追蹤強反射識別與分離方法

1.1 基于局部頻率約束的動態(tài)匹配追蹤算法

傳統(tǒng)的匹配追蹤算法是一種貪婪算法[15]，通過創(chuàng)建超完備的時頻原子庫，對時移、頻率和相位參數(shù)進行全局掃描，以尋找最佳相關(guān)的匹配原子，但該算法每一次迭代的內(nèi)積運算量大，對信號的分解效率極低。動態(tài)匹配追蹤算法將信號的瞬時特征作為先驗信息約束匹配原子的搜索范圍[16]，以先驗信息創(chuàng)建動態(tài)子波庫，使每一次迭代的內(nèi)積運算量大為減小，極大地提高了分解效率。傳統(tǒng)方法以瞬時頻率確定最優(yōu)原子搜索范圍，而瞬時頻率是瞬時相位對時間的導數(shù)，計算結(jié)果存在“負頻率”現(xiàn)象，且受噪聲干擾大，計算結(jié)果極不穩(wěn)定。當頻率計算結(jié)果不合理時，只能對時頻原子的頻率屬性進行全局搜索，并不是真正意義上的動態(tài)搜索。為此，本文將局部頻率算法用于動態(tài)匹配追蹤算法，根據(jù)整形光滑理論，利用差分算子綜合某數(shù)據(jù)點鄰域內(nèi)的頻率值，計算得到該數(shù)據(jù)點的頻率屬性信息，局部頻率計算結(jié)果合理、穩(wěn)定，抗噪性較強,因此可以直接用于確定最優(yōu)時頻原子的搜索鄰域。

連續(xù)信號x(t)的復地震道為[17]

c(t)=x(t)+ih(t)=A(t)eiφ(t)

(1)

式中：h(t)為x(t)的Hilbert變換；A(t)為地震道包絡；φ(t)為信號的瞬時相位。

傳統(tǒng)的瞬時角頻率ω(t)是瞬時相位φ(t)的時間變化率[18]，即

(2)

引入Fomel[19]的局部角頻率ωloc的計算策略

ωloc=[λ2I+S(V-λ2I)]-1SU

(3)

式中：λ為平滑尺度因子；I為單位向量；S為整形正則化算子[20]；U和V分別為u(t)和v(t)組成的向量矩陣。

為了驗證局部頻率計算結(jié)果的合理性，設計兩個用于對比的瞬時頻率與局部頻率信號(圖1)。圖2為線性調(diào)諧信號與實際地震信號的瞬時頻率。由圖可見，由于瞬時頻率是相位對時間的導數(shù)，計算結(jié)果出現(xiàn)很多“負頻率”現(xiàn)象，利用負頻率奇異點不能直接確定動態(tài)匹配追蹤頻率參數(shù)搜索范圍，需要進行全局掃描，因此不是真正意義上的動態(tài)搜索方式。圖3為線性調(diào)諧信號與實際地震信號的局部頻率。由圖可見，線性調(diào)諧信號的局部頻率計算結(jié)果與理論值一致(圖3a)，實際地震信號的局部頻率計算結(jié)果較為合理(圖3b)，在理論范圍內(nèi)不存在“負頻率”現(xiàn)象，可以用于確定最佳匹配原子的頻率參數(shù)。

圖1 線性調(diào)諧信號(a)與實際地震信號(b)圖a的頻率由70Hz向50Hz減小

圖2 線性調(diào)諧信號(a)與實際地震信號(b)的瞬時頻率

圖3 線性調(diào)諧信號(a)與實際地震信號(b)的局部頻率

首先確定最大包絡振幅的時間位置u0，即為此次迭代搜索的中心，進而得到該時間位置u0處的局部角頻率ω0和相位φ0，以此創(chuàng)建匹配搜索的先驗信息點α0(u0,ω0,φ0)，并設定相鄰控制參數(shù)的間隔Δα，得到最佳匹配原子的搜索范圍γ′=[α0-nΔα,α0+nΔα](n為正整數(shù),α為時頻原子的三個控制參數(shù)u、ω和φ,Δα為相鄰控制參數(shù)α的間隔)，即得到迭代搜索的動態(tài)字典原子庫D′={gγ′(t)}(D′∈D,gγ′(t)為時頻原子),D為過完備字典原子庫。

假定第1次迭代需要掃描動態(tài)字典原子庫中的m個時頻原子

φi∈U[φ(u0),δφ]}

(4)

式中：u0為每個最佳匹配子波的搜索中心位置處的時間；ω(u0)和φ(u0)分別為搜索中心位置處信號的瞬時角頻率和相位信息；U[ω(u0),δω]、U[φ(u0),δφ]分別為角頻率和相位的搜索鄰域,δω、δφ為相應參數(shù)搜索半徑。

第1次迭代與信號最佳相關(guān)的時頻原子gγ1(t)滿足

(5)

則信號f被表示為沿gγ1(t)方向的分量和殘差分量，即

f=〈f,gγ1〉gγ1+R1f

(6)

式中R1f為第一次迭代后的殘差信號。按上述方式依次迭代下去，直至達到預設的迭代次數(shù)和殘差閾值，則迭代停止。若經(jīng)過n次迭代后，完成對信號的分解過程，最終將信號表示為n個最佳匹配原子的線性組合

(7)

式中： Ri-1f為i-1次迭代后的殘差信號； Rnf為n次迭代后的殘差信號，當i=1時,R0f為原始信號。

利用9個Morlet小波合成理論信號，基于局部頻率約束的匹配追蹤算法迭代分解理論信號，圖4為無噪聲理論信號的迭代分解過程。由圖可見，每一次迭代優(yōu)選出一個最佳匹配的原子，經(jīng)過9次迭代后，重構(gòu)信號與理論信號高度吻合，殘差曲線為零值。為驗證本文提出的匹配追蹤算法的抗噪性，對理論合成信號進行加噪，在信噪比為5:1的情況下，對加噪信號迭代分解，圖5為加噪理論信號的迭代分解過程。由圖可見，基于局部頻率約束的匹配追蹤算法對信號的迭代分解過程穩(wěn)定，經(jīng)過9次迭代后，得到9個最佳匹配原子重構(gòu)理論信號，殘差曲線基本為噪聲信號，證明該算法具有一定抗噪性。

1.2 基于局部頻率約束的動態(tài)匹配追蹤算法

在對原始數(shù)據(jù)進行強反射識別與分離時，可以利用層位約束或一定的時窗范圍內(nèi)信號包絡振幅最大位置拾取強反射具體位置。通過求取強反射時間位置的復地震道的局部頻率、瞬時相位信息，構(gòu)建局部反射特征波形庫，利用式(5)搜索最佳相關(guān)的子波作為強反射子波。在對實際數(shù)據(jù)進行強反射剝離前，需要先確定強反射剝離的強度[21]，即

(8)

式中：dN為強反射時窗內(nèi)地震數(shù)據(jù)能量,dM為背景時窗內(nèi)地震數(shù)據(jù)能量,N和M分別代表其時窗內(nèi)地震數(shù)據(jù)的采樣點數(shù)；ξ為強反射壓制參數(shù)；dstrong為拾取的強反射地震數(shù)據(jù)。

圖4 無噪聲理論信號的迭代分解過程

圖5 加噪理論信號的迭代分解過程

圖6為去除強反射前、后記錄。根據(jù)壓制強反射的能量均衡原則(式8)，認為使去除強反射后的平均能量(紅色虛線框)與背景能量(藍色虛線框)基本保持一致的ξ為最佳強反射壓制系數(shù)，能夠達到較為理想的強反射壓制效果(圖6b)。圖7為M區(qū)不同ξ的單道處理結(jié)果。為了滿足能量均衡原則(式8)，即使去除強反射后的平均能量(紅色虛線框內(nèi))與地震數(shù)據(jù)背景能量(藍色虛線框內(nèi))保持一致，通過對比選取ξ=0.6為最佳強反射壓制系數(shù)，以此系數(shù)對實際地震數(shù)據(jù)進行強反射剝離，能夠較為理想地壓制強反射并突顯目標儲層弱反射信息(圖7d)。

圖6 去除強反射前(a)、后(b)記錄

紅色虛線框為含強反射的時窗，藍色虛線框為選取的背景能量時窗，黑線為原始地震記錄，綠線為根據(jù)強反射分析后拾取的最佳匹配子波，紅線為去強反射后的地震記錄

圖7 M區(qū)不同ξ的單道處理結(jié)果

2 模型試算

將本文提出的基于優(yōu)化局部頻率約束的動態(tài)匹配追蹤強反射識別與分離方法用于模型數(shù)據(jù)的處理，以驗證方法可靠性。構(gòu)建含高阻層的理論地質(zhì)模型(圖8)，利用28Hz的Morlet小波合成理論模型的正演地震記錄(圖9)，可見由于受高速層的影響，頂部薄砂體的有效反射信息被強反射淹沒，難以識別上覆砂體真正的形態(tài)和位置。圖10為由本文方法識別和剝離強反射后的地震記錄。由圖可見，明顯地削弱了強干擾的影響，能夠恢復上覆薄砂巖的地震響應。為了更清晰地看出強反射對砂體特征識別的影響，圖11和圖12分別為去除強反射前、后的地震道集。由圖可見，剝離強反射后突顯了砂體形態(tài)特征，有助于對儲層進行精細描述，為后續(xù)地震數(shù)據(jù)處理提供了有效的數(shù)據(jù)基礎(圖12)。

圖8 理論地質(zhì)模型

模型頂部(灰色)和底部(綠色)為大套泥巖(平均速度為4300m/s，密度為2.7g/cm3)，中部(黑色)為因泥巖脫水形成的高阻層(平均速度為4700m/s，密度為2.8g/cm3)，高阻層上覆薄層砂巖(紅色，平均速度為4350m/s，密度為2.5g/cm3)

圖9 理論模型的正演地震記錄

圖10 由本文方法識別和剝離強反射后的地震記錄

圖11 模型原始正演道集

圖12 去除強反射后的地震道集

3 實際資料處理

A區(qū)位于濟陽坳陷，該區(qū)存在多級不整合，利于發(fā)育與不整合有關(guān)的油藏。實際地震資料存在典型的不整合反射特征，為準確描述巖性圈閉提供了基礎數(shù)據(jù)。由于該區(qū)儲層傾角小、儲層厚度小且目的層上、下存在高阻抗體，地震反射特征不能真實反映不整合圈閉的地質(zhì)特征。為了準確地拾取砂體反射特征信息，有必要識別與分離強反射。

圖13為原始地震記錄。由圖可見，砂巖尖滅反射信息淹沒于上覆高速層強反射中，難以識別真正的砂體尖滅點。考慮到強反射記錄具有很好的橫向連續(xù)性，利用本文提出方法將T2和Tr強反射從工區(qū)實際數(shù)據(jù)中剝離，從剝離后的地震記錄(圖14)上可見，砂體的反射信息得到突顯。從去除強反射前、后的地震道集上也能明顯地看到：由于Tr強反射與砂體反射干涉、疊加，無法準確識別砂體尖滅位置(圖15)； 剝離Tr和T2強反射后砂體尖滅反射信息更清楚(圖16)。因此，高速層強反射識別與剝離技術(shù)為砂體尖滅識別提供了可靠的數(shù)據(jù)。

圖13 原始地震記錄

圖14 去除強反射的地震記錄

圖15 原始地震道集

圖16 去除強反射的地震道集

4 結(jié)束語

為削弱目標儲層上、下高阻層形成的強反射干擾影響，本文提出了基于局部頻率約束的動態(tài)匹配追蹤強反射識別與分離方法?？紤]到瞬時頻率計算結(jié)果存在“負頻率”現(xiàn)象，認為其不適用于作為時頻原子的頻率搜索參數(shù)，本文提出局部頻率計算模式，計算結(jié)果較為合理，將其作為構(gòu)建動態(tài)子波庫的控制參數(shù)。利用強反射位置處的局部頻率和瞬時相位信息作為先驗信息，構(gòu)建動態(tài)反射特征波形庫，以此搜索最佳匹配強反射子波，通過單道測試確定強反射壓制系數(shù)，以保證剝離強反射后的均衡能量與背景能量一致，有助于拾取目標儲層的反射信息，最終實現(xiàn)儲層的精細描述。通過構(gòu)建含高阻層的模型數(shù)據(jù)驗證了基于局部頻率約束的動態(tài)匹配追蹤強反射識別和分離方法的可行性，并將該方法應用于實際地震數(shù)據(jù)，削弱了強反射干擾影響，突顯了弱反射信息。