復(fù)數(shù)域約束最小二乘拓頻

楊 培 杰

(中國石化勝利油田分公司勘探開發(fā)研究院，山東東營 257015)

0 引言

地震拓頻技術(shù)大致可分為兩類，一類是時間域拓頻，主要包括反褶積[1]、盲反褶積[2]、反Q濾波[3]、盲源分離[4]、壓縮感知[5]等。反褶積通過壓縮地震子波提高時間分辨率，計算過程中往往需要假設(shè)地震子波為最小相位，以及反射系數(shù)為高斯白噪，或是需要子波的參與。為了克服這些問題，盲反褶積應(yīng)運而生，盲反褶積不需要或者是弱化了對子波和反射系數(shù)的先驗假設(shè)，因而可以應(yīng)用于時變系統(tǒng)和非最小相位系統(tǒng)，但是由于其對噪聲較為敏感，因此并沒有得到較好的應(yīng)用。反Q濾波通過補償?shù)貙拥酿椥运p，從而提高地震資料的分辨率，但易受噪聲、Q值誤差等多種因素影響。盲源分離利用圍巖反射在各地震道的相似性和目標(biāo)儲層弱信號的差異性分離地震弱信號，以達到提高地震資料分辨率的目的，往往需要對相鄰的地震道做出一定的假設(shè)。

另一類是頻率域拓頻，主要包括譜白化[6]、譜藍(lán)化[7]、有色反演[8]、譜反演[9-10]等。譜白化通過對地震數(shù)據(jù)頻譜進行補償，以達到拓寬地震頻帶的目的。但現(xiàn)實中很難得到完全白化的頻譜，因此譜藍(lán)化技術(shù)得到了更為廣泛的應(yīng)用。有色反演實際上是一種測井約束頻率域拓頻，其核心是用地震的頻譜和井的波阻抗頻譜相匹配，從而拓寬地震數(shù)據(jù)的頻譜。譜反演是在譜分解的基礎(chǔ)上，通過反演使頻率域目標(biāo)函數(shù)達到極小而得到反射系數(shù)的方法，由于譜反演法能分辨調(diào)諧厚度以下的薄層，因而得到了較為廣泛的應(yīng)用。

最小二乘法在地震勘探中的應(yīng)用十分廣泛。約束最小二乘地震反演[11]使用最小二乘方法構(gòu)建目標(biāo)函數(shù)，反演過程穩(wěn)定、計算效率高；最小二乘偏移[12-13]把成像問題當(dāng)作一個反問題處理，以尋求更接近于真實的地下反射；最小二乘橫波估計[14]通過最小二乘將橫波估計目標(biāo)函數(shù)的最優(yōu)化問題轉(zhuǎn)化為求解線性矩陣方程組問題，從而大大提高了橫波估計的效率和穩(wěn)定性；最小二乘時頻分析方法[15-17]將時頻分析轉(zhuǎn)換成一個求解反問題的過程，通過求解矩陣方程組進而得到信號的高分辨率時頻分析結(jié)果。

本文開發(fā)了復(fù)數(shù)域約束最小二乘譜藍(lán)化拓頻(Constrained Complex-domain Least-squares Spectrum Blueing，簡記為Y-SpecB)技術(shù)，其核心思想是設(shè)計一個寬頻的約束目標(biāo)譜[18]，在復(fù)數(shù)域最小二乘[19-21]的思路下，將原始地震頻譜向約束目標(biāo)譜靠近，進而達到提高地震資料主頻、拓寬頻帶的目的。

1 數(shù)學(xué)模型及其解

考慮復(fù)數(shù)域矩陣方程

E(ω)D(ω)+N(ω)=F(ω)

(1)

式中：E(ω)為頻率域拓頻算子；D(ω)為地震數(shù)據(jù)的譜；F(w)為寬頻約束目標(biāo)譜；N(ω)為白噪聲的譜；ω為頻率。其中D(ω)和F(ω)為已知項，E(ω)為未知待求項。

式(1)即為Y-SpecB的正演數(shù)學(xué)模型，其意義為，用拓頻算子去拓寬原始地震譜，加上白噪聲譜，從而逼近寬頻約束目標(biāo)譜。

復(fù)數(shù)矩陣E(ω)、D(ω)、F(ω)的具體形式為

(2)

(3)

(4)

式中：下標(biāo)“R”表示復(fù)數(shù)的實部；下標(biāo)“I”表示復(fù)數(shù)的虛部；n為待拓頻地震道的頻率域樣點數(shù)。

將E(ω)定義為下式的最優(yōu)解

(5)

顯然，該式為復(fù)數(shù)域的最優(yōu)化問題。令

(6)

則式(1)可以寫為

(7)

用矩陣方程組式(7)中的系數(shù)矩陣的轉(zhuǎn)置左乘該矩陣方程，有

(8)

進一步可表示為

(9)

式中

(10)

則式(9)的解可以寫為

(11)

實際計算中，為了提高解的穩(wěn)定性，需要對目標(biāo)函數(shù)的最優(yōu)化過程進行約束[22-24]，則目標(biāo)函數(shù)式(5)變?yōu)?/p>

(12)

式中Regu[E(ω)]為正則化項。由于L2范數(shù)正則化可以有解析形式的解，易于求解，因此在求解式(12)時選用L2范數(shù)正則化約束，則式(11)變?yōu)?/p>

(13)

式中：μ≥0，為控頻因子；I為單位對角矩陣。

將式(13)的解表示為

EO(ω)=ER(ω)+iEI(ω)

(14)

將求解的拓頻算子EO(ω)與原始地震頻譜D(ω)相乘，可得拓頻后的地震數(shù)據(jù)頻譜

DY(ω)=EO(ω)D(ω)

(15)

對其進行傅里葉反變換，可得拓頻后的時間域地震數(shù)據(jù)

dY=IFFT[DY(ω)]

(16)

式中IFFT表示傅里葉反變換。

2 寬頻約束目標(biāo)譜的設(shè)計

寬頻約束目標(biāo)譜的設(shè)計是本文方法的關(guān)鍵點之一，常用兩種目標(biāo)譜的設(shè)計方法，一是理論目標(biāo)譜，二是實際目標(biāo)譜。

理論目標(biāo)譜基于各種窗函數(shù)進行設(shè)計，常用的理論目標(biāo)譜如圖1所示。一般來說，海寧窗和高斯窗具有窄頻的特征，梯形窗和廣義高斯窗具有寬頻的特征，而梯形窗函數(shù)具有頻率突變點，不利于最小二乘的擬合效果。理論目標(biāo)譜的優(yōu)點是設(shè)計靈活，可控性高。

為了更好地拓展地震資料的高低頻，一般會采用廣義高斯窗形狀的理論目標(biāo)譜(圖1d)，廣義高斯目標(biāo)譜Gg定義為

圖1 常用理論目標(biāo)譜(a)海寧窗；(b)梯形窗；(c)高斯窗；(d)廣義高斯窗

(17)

式中：ωL表示低截頻率；ωH表示高截頻率；σL表示低截標(biāo)準(zhǔn)差；σH表示高截標(biāo)準(zhǔn)差。

實際目標(biāo)譜的設(shè)計來自具體研究區(qū)實際測井資料，多采用波阻抗的頻譜作為實際目標(biāo)譜，這和有色反演思路較為類似。圖2為M地區(qū)一口井波阻抗的頻譜(去掉部分低頻信息后)，可以看出，井資料頻譜很寬，即使到了500Hz仍然有能量，并且振幅譜的變化較大，因此，實際目標(biāo)譜往往不太穩(wěn)定和統(tǒng)一。

圖2 實際井資料的波阻抗頻譜

上述約束目標(biāo)譜的設(shè)計主要是針對目標(biāo)譜的振幅譜。Y-SpecB方法不改變拓頻后地震數(shù)據(jù)的相位譜，具體做法是：首先對地震道進行傅里葉變換，然后保留相位譜，將相位譜所對應(yīng)的振幅譜修正為所設(shè)計的目標(biāo)譜的振幅譜，即可得到包含了振幅譜和相位譜的約束目標(biāo)譜。此時，通過反演求得的拓頻算子為零相位，在時間域?qū)ΨQ。

需要說明的是，不論采用理論目標(biāo)譜還是實際目標(biāo)譜，對拓頻結(jié)果的影響并不大，這主要是因為本文方法并不是讓拓頻后的頻譜去擬合目標(biāo)譜，而是在控頻因子μ的控制下向目標(biāo)譜靠近，拓展原始地震數(shù)據(jù)微弱的低頻和高頻信息，進而拓寬頻帶，提高地震資料的分辨率。

實際應(yīng)用中，可以嘗試不同的理論目標(biāo)譜和實際目標(biāo)譜分別進行拓頻，從中選擇效果更好的目標(biāo)譜，進而開展全區(qū)的拓頻處理。

3 模型驗證

通過模型數(shù)據(jù)說明本文方法的準(zhǔn)確性和有效性。設(shè)計如圖3所示的抽象前積體地質(zhì)模型[25]，紅色和黃色表示前積砂體，綠色表示泥巖，巖性地層單元的時間厚度約20ms，前積體時間厚度約5ms，共11期。

圖3 前積體地質(zhì)模型

分別采用不同主頻的Ricker子波進行地震正演，結(jié)果如圖4所示?？梢钥闯?，30Hz合成記錄主頻較低，頻帶較窄，無法分辨巖性界面內(nèi)部的11期前積砂體，而只能反映大套的巖性地層界面，50Hz合成記錄主頻較高，頻帶較寬，可以較好地分辨這11期前積砂體。

圖4 前積體模型不同主頻子波地震正演結(jié)果(左)及其頻譜(右)(a)30Hz；(b)50Hz

選擇廣義高斯目標(biāo)譜，用本文方法對30Hz合成記錄進行拓頻處理，并分析不同的控頻因子μ對拓頻結(jié)果的影響。經(jīng)過反復(fù)測試，發(fā)現(xiàn)當(dāng)ωL=18Hz、ωH=110Hz、σL=10Hz、σH=30Hz時，拓頻效果最佳。圖5為使用不同控頻因子μ拓頻結(jié)果，可以看出，本文方法不僅能增加高頻端的能量，同時也能拓展低頻段的能量，隨著μ變小，拓頻后地震記錄的主頻逐漸變高，頻帶逐漸變寬，前積砂體也逐漸變得清晰可識別；當(dāng)μ=0.0001時，拓頻結(jié)果能夠較為真實地恢復(fù)50Hz的合成記錄，巖性界面中的10期砂體得以重建，如圖6所示。

圖5 前積體模型不同控頻因子μ的拓頻后剖面(左)及CDP=180單道頻譜(右)的對比(a)μ=0.01；(b)μ=0.001；(b)μ=0.0001

圖6 50Hz記錄(上)與30Hz記錄拓頻后(下)的對比

頻率域拓頻算子類似于一個帶通濾波器，對其進行傅里葉反變換，可以得到時間域拓頻算子，如圖7所示。該算子關(guān)于零時刻對稱，相位為零，保證了拓頻不會改變原始地震數(shù)據(jù)的相位信息，進而保證拓頻后的地震數(shù)據(jù)的準(zhǔn)確性。

圖7 時間域拓頻算子

進一步分析拓頻前、后數(shù)據(jù)與高頻合成記錄的相關(guān)性。30Hz合成地震道在經(jīng)過拓頻后同相軸增多，與50Hz合成地震道的同相軸具有很好的一致性(圖8)；30Hz合成地震道與50Hz合成地震道互相關(guān)最大值為0.45，而30Hz拓頻后地震道與50Hz合成地震道互相關(guān)最大值為0.87，相關(guān)性提高了0.93倍(圖9)。由圖10可以看出，30Hz合成地震道相位譜和拓頻后地震道相位譜在有效頻帶范圍內(nèi)并沒有改變，進而保證了拓頻結(jié)果的可靠性。

圖8 30Hz合成記錄拓頻前、后與50Hz合成記錄的波形對比

圖9 30Hz合成記錄拓頻前、后與50Hz合成記錄相關(guān)性分析(CDP165)

圖10 30Hz合成記錄拓頻前(黑色)、后(紅色)相位譜對比

關(guān)于約束目標(biāo)譜參數(shù)ωL、ωH、σL、σH及控頻因子μ如何設(shè)置，并沒有一個具體的要求，一般的原則是，約束目標(biāo)譜要完全包含原始地震道的頻譜，為保證拓頻效果，ωH可適當(dāng)高一些；在此基礎(chǔ)上，μ越小，拓頻后的頻譜越接近目標(biāo)譜。缺省值為ωL=18Hz、ωH=100Hz、σL=10Hz、σH=30Hz、μ=0.0001，在此基礎(chǔ)上，可以不斷調(diào)試這些參數(shù)，直到達到滿意的拓頻效果。

模型試算結(jié)果表明，本文方法可以在基本不改變相位譜的基礎(chǔ)上，通過調(diào)節(jié)控頻因子μ，有效地補償?shù)卣饠?shù)據(jù)的低頻信息，增強高頻信息，保證了拓頻結(jié)果的客觀、準(zhǔn)確性。

4 實際應(yīng)用

4.1 薄互層識別

濟陽坳陷CD油田古近系東營組4砂組(Ed4)以濁積沉積為主，儲量豐富，是尋找優(yōu)質(zhì)儲量的重點攻關(guān)層系。該層較深，和泥巖呈互層分布，總體較薄，由部分井測井曲線可知，平均速度約為3600m/s。在地震剖面上目的層主頻約為25Hz，地震數(shù)據(jù)的分辨率較低，大約為18m，同相軸出現(xiàn)相互疊置現(xiàn)象，導(dǎo)致大多砂體無法有效分辨。

在開展Y-SpecB拓頻前，首先要設(shè)計目標(biāo)譜。對地震數(shù)據(jù)的主頻、帶寬和有效信號占比(ESR)進行分析，最終選擇廣義高斯窗的理論目標(biāo)譜，如圖11所示，參數(shù)ωL=15Hz、ωH=80Hz、σL=10Hz、σH=30Hz。對μ的取值進行試驗，反復(fù)設(shè)定不同值并觀察拓頻后地震數(shù)據(jù)的保真度和ESR，最終確定μ=0.0001。

圖11 CD油田約束目標(biāo)譜、拓頻算子、原始單道及拓頻后的頻譜對比

圖12為Line1783拓頻前、后的剖面及頻譜對比，可以看出，拓頻后地震數(shù)據(jù)主頻由28.5Hz提高到51Hz，帶寬由53Hz增加到102Hz，由該區(qū)地層速度計算，拓頻后的地震資料對于地層的絕對分辨率由18m左右提高到8.7m左右，大大提高了地震資料的分辨能力，拓頻后的地震剖面能量均衡，中、高頻信息豐富，分辨率明顯提高(箭頭所示)，總體效果較好。

圖12 拓頻前(a)、后(b)Line1783地震剖面(左)及其頻譜(右)的對比

圖13為過cb81井地震剖面拓頻前、后對比，其中黑色井曲線為波阻抗，目的層段速度具有砂高泥低的特點，黑色箭頭處為幾套薄層砂體或砂泥巖薄互層的組合。從上往下的第一、第三箭頭處砂體組合，由于儲層之間的泥巖隔層較薄，在原始剖面上兩套砂巖與泥巖隔層的同相軸出現(xiàn)了相互干涉的現(xiàn)象，而在拓頻后的地震剖面上，這兩套砂體被清楚地分辨了出來。中間箭頭處為一套砂泥巖薄互層，在原始剖面上表現(xiàn)為一個地震同相軸，而在拓頻后的地震剖面上，隱約可以看到多個同相軸，分辨率得到了明顯提高。

圖13 過cb81井地震剖面拓頻前(a)、后(b)的對比

進一步通過井震標(biāo)定分析本文方法拓頻結(jié)果的準(zhǔn)確性。根據(jù)拓頻后數(shù)據(jù)的主頻范圍，選擇主頻45Hz的Ricker子波制作合成記錄，進行井震標(biāo)定。從圖14中的紅色箭頭處可以看出，在拓頻后的剖面上可以較為清晰的看到原始地震剖面上沒有的地震同相軸。將合成記錄和拓頻前、后的地震剖面進行對比，合成記錄和拓頻后的地震剖面的相關(guān)性更好，相關(guān)系數(shù)由0.64提高到了0.86，即用拓頻后的地震數(shù)據(jù)進行井震標(biāo)定，可提高標(biāo)定結(jié)果的準(zhǔn)確性，進一步說明本文方法拓頻結(jié)果的客觀性、準(zhǔn)確性。

圖14 cb81井拓頻前、后井震標(biāo)定對比

4.2 超覆尖滅點識別

濟陽坳陷GB洼陷南斜坡古近系沙三段(Es3)以發(fā)育地層超覆油氣藏為主，已發(fā)現(xiàn)多套含油層，是該區(qū)最主要的儲量層系。在構(gòu)造研究基礎(chǔ)上，開展地層展布規(guī)律分析，通過北東向的井間地層對比，發(fā)現(xiàn)在古地貌坡折處，井間存在地層超覆尖滅現(xiàn)象。

通過分析和試算，選擇Zh243井的實際目標(biāo)譜作為約束目標(biāo)普。由于實際目標(biāo)譜的變化較大且不穩(wěn)定(圖2)，因此對原始的實際目標(biāo)譜進行了去直流和平滑處理，結(jié)果如圖15所示。經(jīng)過反復(fù)實驗，取μ=0.001。

圖15 GB洼陷原始地震單道頻譜、約束目標(biāo)譜、拓頻算子及拓頻后頻譜對比

為驗證本文方法的客觀性、準(zhǔn)確性，對拓頻后地震剖面進行梯形帶通濾波，濾波器參數(shù)設(shè)置為0～20～30～50Hz，帶通濾波后的地震數(shù)據(jù)(圖17)和原始數(shù)據(jù)(圖16a)非常接近。

圖16 拓頻前(a)、后(b)連井剖面(左)及其頻譜(右)的對比

圖17 拓頻結(jié)果的帶通濾波后地震剖面

5 結(jié)束語

本文建立了復(fù)數(shù)域約束最小二乘拓頻方法的正演數(shù)學(xué)模型，在寬頻目標(biāo)譜的約束下，通過復(fù)數(shù)域最小二乘方法，實現(xiàn)了地震數(shù)據(jù)頻率域的拓頻處理。模型驗證表明，本文方法可以在不改變地震信號相位譜的情況下，有效增強地震數(shù)據(jù)的低頻和高頻信息，并可以通過設(shè)定不同的控頻因子來調(diào)節(jié)拓頻后地震數(shù)據(jù)的頻譜和目標(biāo)頻譜之間的距離，以得到不同分辨率的拓頻結(jié)果。實際應(yīng)用結(jié)果表明，本文方法可以提高地震主頻、拓寬頻帶，有效地提高了地震資料對薄互層以及尖滅點的識別能力，可以為地震地質(zhì)綜合研究提供高分辨率的地震資料。