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

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

0 引言

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

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

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

1 匹配追蹤強(qiáng)反射識別與分離方法

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

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

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

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

(1)

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

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

(2)

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

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

(3)

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

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

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

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

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

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

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

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

(4)

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

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

(5)

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

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

(6)

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

(7)

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

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

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

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

(8)

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

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

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

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

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

紅色虛線框?yàn)楹瑥?qiáng)反射的時(shí)窗，藍(lán)色虛線框?yàn)檫x取的背景能量時(shí)窗，黑線為原始地震記錄，綠線為根據(jù)強(qiáng)反射分析后拾取的最佳匹配子波，紅線為去強(qiáng)反射后的地震記錄

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

2 模型試算

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

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

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

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

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

圖11 模型原始正演道集

圖12 去除強(qiáng)反射后的地震道集

3 實(shí)際資料處理

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

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

圖13 原始地震記錄

圖14 去除強(qiáng)反射的地震記錄

圖15 原始地震道集

圖16 去除強(qiáng)反射的地震道集

4 結(jié)束語

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