劉金濤 王 孝 王小衛(wèi) 蘇 勤 劉夢(mèng)麗 袁 煥

(中國(guó)石油勘探開(kāi)發(fā)研究院西北分院，甘肅蘭州 730030)

0 引言

眾所周知，地震波在地下介質(zhì)中傳播時(shí)，由于大地濾波作用，會(huì)導(dǎo)致地震波能量衰減、速度頻散，降低分辨薄層的能力。目前已提出了許多地層吸收衰減模型(Q模型)描述這一現(xiàn)象，其中常Q模型是地震勘探領(lǐng)域使用最多的地層吸收衰減模型[1]。地震波的衰減效應(yīng)破壞了地震子波在時(shí)間方向的穩(wěn)態(tài)特征，導(dǎo)致傳統(tǒng)的反褶積方法很難從非穩(wěn)態(tài)地震記錄中準(zhǔn)確地恢復(fù)反射系數(shù)序列。因此，基于吸收補(bǔ)償?shù)奶岣叻直媛史椒ň哂忻鞔_的物理意義和較大應(yīng)用潛力[2-8]。

反Q濾波是吸收補(bǔ)償?shù)闹饕椒?，有多種實(shí)現(xiàn)方式，其中，波場(chǎng)延拓和稀疏反演是兩類應(yīng)用最為廣泛的方法。由于波場(chǎng)延拓方法在外推過(guò)程中存在放大函數(shù)，極易放大高頻噪聲，結(jié)果不穩(wěn)定。因此，稀疏反演類反Q濾波方法應(yīng)運(yùn)而生[9-10]。Zhang等[11]將衰減補(bǔ)償看作基于反演的最小二乘問(wèn)題，并用貝葉斯理論進(jìn)行正則化約束；Braga等[12]在小波域?qū)崿F(xiàn)了L2范數(shù)約束的反Q濾波；Oliveira等[13]在頻率域?qū)崿F(xiàn)了L1范數(shù)約束的衰減補(bǔ)償，Wang等[14]對(duì)其進(jìn)行改進(jìn)，在時(shí)間域?qū)崿F(xiàn)了L1范數(shù)約束的衰減補(bǔ)償。以上方法都只考慮了地震道縱向上的信息，而沒(méi)有考慮相鄰道之間的相關(guān)性，即都是單道補(bǔ)償方法，因而結(jié)果往往會(huì)出現(xiàn)橫向不連續(xù)問(wèn)題。為了解決這一問(wèn)題，將相鄰道之間的相關(guān)性作為先驗(yàn)信息加入反演系統(tǒng)進(jìn)行多道反演補(bǔ)償?shù)姆椒ㄖ饾u發(fā)展起來(lái)。

Zhang等[15]提出了一種多道反演方法，該算法將橫向約束融入反演算法中，反演的結(jié)果較常規(guī)單道算法具有更好的橫向連續(xù)性。Kazemi等[16]假設(shè)反射系數(shù)序列稀疏分布，提出了一種非線性多道盲源反褶積方法，主要優(yōu)點(diǎn)是抗噪性較強(qiáng)且不需要子波的先驗(yàn)信息。Yuan等[17]先將數(shù)據(jù)轉(zhuǎn)化到變換域，然后在變換域數(shù)據(jù)上加一個(gè)約束項(xiàng)，實(shí)現(xiàn)了多道波阻抗反演，在一定程度上緩解了單道反演算法的橫向不連續(xù)性問(wèn)題；在此基礎(chǔ)上，Yuan等[18]提出了一種在頻率—波數(shù)域的多道衰減補(bǔ)償方法，增強(qiáng)了地震同相軸橫向連續(xù)性，提高了吸收補(bǔ)償?shù)姆€(wěn)定性和抗噪性。

f-x域預(yù)測(cè)濾波是描述線性事件空間預(yù)測(cè)特性的有用工具，主要被用于數(shù)據(jù)插值和去噪領(lǐng)域[19-22]。為提高地層吸收補(bǔ)償?shù)姆€(wěn)定性和抗噪性，本文提出了一種基于反演的全局約束吸收補(bǔ)償方法。該方法首先在f-x域計(jì)算信號(hào)預(yù)測(cè)算子[23-25]，用于表征地震信號(hào)在空間上的結(jié)構(gòu)關(guān)系和能量變化；然后，將該算子作為橫向正則化約束引入到縱向吸收補(bǔ)償系統(tǒng)，建立具有雙向(縱向和橫向)即全局約束的反演目標(biāo)函數(shù)；最后應(yīng)用模型數(shù)據(jù)和實(shí)際數(shù)據(jù)驗(yàn)證了方法的有效性。

1 地層吸收模型

1.1 Q的定義

當(dāng)平面波通過(guò)均勻的黏彈性介質(zhì)時(shí)，其振幅衰減和速度頻散可以用一個(gè)無(wú)量綱的參數(shù)Q描述

(1)

式中：ω為圓頻率；α(ω)和v(ω)分別為衰減系數(shù)和相速度。由于品質(zhì)因子Q必須為正數(shù)，所以

(2)

進(jìn)一步假設(shè)

(3)

則式(1)可以近似為

(4)

1.2 修正的Kolsky-Futterman模型

Kolsky[26]假設(shè)在測(cè)量頻帶范圍內(nèi)衰減系數(shù)α(ω)與頻率呈嚴(yán)格的線性關(guān)系，即

(5)

式中ω0為參考圓頻率。同時(shí)，Kolsky定義其相速度為

(6)

將式(5)和式(6)代入式(4)，可得

(7)

對(duì)于較大的品質(zhì)因子，即Q(ω0)?1，有

(8)

式中γ=1/[πQ(ω0)]。則相速度和品質(zhì)因子可近似為

(9)

(10)

其中自由選擇的參考頻率ω0小于最低的測(cè)量頻率。

Wang等[27]指出，式(6)相速度僅僅是對(duì)于ω?ω0的一個(gè)漸近表達(dá)式。由于地震勘探資料的頻率范圍相對(duì)較低(一般小于500Hz)，其修正的相速度公式為

(11)

式中：h是一個(gè)頻率不相關(guān)的常數(shù)；ωh為重新定義的調(diào)諧參數(shù)。

對(duì)于修正的相速度模型，其相應(yīng)的修正Q模型為

(12)

本文所用的地層吸收模型即為上述模型。

2 基于反演的縱向約束吸收補(bǔ)償方法

在完全彈性介質(zhì)中，可以用褶積模型簡(jiǎn)單描述地震數(shù)據(jù)采集過(guò)程，即地震記錄可由地震子波和反射系數(shù)序列褶積外加噪聲得到

d(t)=r(t)*w(t)+n(t)

(13)

式中：d(t)為地震記錄；r(t)為反射系數(shù)；w(t)為震源子波；n(t)為噪聲。

在黏彈性介質(zhì)中，由于地層的吸收衰減作用，導(dǎo)致震源子波會(huì)隨著傳播時(shí)間而變化。因此，需要將穩(wěn)態(tài)的褶積模型推廣到非穩(wěn)態(tài)

(14)

式中w(τ，t)為時(shí)變子波，且有w(0，t)=w(t)。

(15)

將式(14)離散可改寫為矩陣向量形式，有

(16)

式中N為地震道的樣點(diǎn)數(shù)。

式(16)的加入縱向約束的L2范數(shù)最優(yōu)化問(wèn)題可用目標(biāo)函數(shù)表示為

(17)

式中：W為時(shí)變子波矩陣；r為反射系數(shù)向量；d為地震記錄向量；μt為縱向約束因子，主要作用為調(diào)節(jié)數(shù)據(jù)殘差項(xiàng)(第一項(xiàng))與正則化項(xiàng)(第二項(xiàng))的權(quán)重，使得反演結(jié)果更合理。通過(guò)極小化目標(biāo)函數(shù)，式(17)的解可表示為

r=(WTW+μtI)-1WTd

(18)

3 基于反演的全局約束吸收補(bǔ)償方法

假設(shè)地震剖面S(t，x)只有一個(gè)振幅為常數(shù)、斜率為p的線性同相軸，轉(zhuǎn)換到頻率域，可表示為

S(f，x)=A(f)exp(j2πfxp)

(19)

式中A(f)是同相軸子波的頻譜。假設(shè)道間隔Δx為常數(shù)，有x=(k-1)Δx，k=1，2，…，K，K為地震剖面的道數(shù)。對(duì)于任意頻率f，有

Sk(f)=Aexp(jθk)

(20)

式中θ=2πfpΔx。則第k道可用第k-1道表示為

Sk(f)=a1(f)Sk-1(f)

(21)

式中a1(f)=exp(jθ)。式(21)遞歸方程是一階差分方程，也稱為一階自回歸(AR)模型。通過(guò)該模型，可以沿著空間方向遞歸地預(yù)測(cè)有效信號(hào)。

類似地，假設(shè)在時(shí)空域有M條線性同相軸，則可以用一個(gè)M階的差分方程表示為

(22)

式中ai(f)表示信號(hào)預(yù)測(cè)算子，主要用于預(yù)測(cè)無(wú)噪聲的疊加復(fù)諧波。式(22)也可寫成預(yù)測(cè)誤差形式

(23)

式中：g0=1；gi=-ai。

在實(shí)際數(shù)據(jù)中，噪聲是不可避免的。假設(shè)地震數(shù)據(jù)由有效信號(hào)和噪聲疊加而成，即

Dk=Sk+Nk

(24)

將上式代入式(23)，可得

(25)

上式為自回歸滑動(dòng)平均(ARMA)模型，其自回歸分量和滑動(dòng)平均分量相等。

(26)

為了更加準(zhǔn)確地預(yù)測(cè)有效信號(hào)，AR模型中的參數(shù)M一般需大于實(shí)際同相軸的條數(shù)，而且M的確定需要應(yīng)用統(tǒng)計(jì)準(zhǔn)則。

假設(shè)M=2，式(26)可以寫成線性方程組

(27)

其矩陣表達(dá)式為

D-Ya=N

(28)

為了最小化噪聲能量，建立目標(biāo)函數(shù)

(29)

對(duì)目標(biāo)函數(shù)求偏導(dǎo)，并令導(dǎo)數(shù)為零，有

YHYa=YHD

(30)

式中上標(biāo)“H”表示共軛轉(zhuǎn)置。需要注意的是，矩陣YHY為Toeplitz矩陣，因此可用Levinson遞歸法快速求解。為了獲得穩(wěn)定的濾波器系數(shù)，通常需要在Toeplitz矩陣的對(duì)角線上增加一個(gè)小的擾動(dòng)μ1，則式(30)改寫為

(YHY+μ1Ι)a=YHD

(31)

那么待求的信號(hào)預(yù)測(cè)算子為

a=(YHY+μ1I)-1YHD

(32)

每一個(gè)頻率的預(yù)測(cè)算子是相互獨(dú)立的，將每個(gè)頻率的一維預(yù)測(cè)算子依次排列，可得一維預(yù)測(cè)算子

(33)

式中I為頻點(diǎn)數(shù)。將每個(gè)頻率的一維預(yù)測(cè)算子作用于含隨機(jī)噪聲的數(shù)據(jù)，可預(yù)測(cè)出有效信號(hào)

(34)

(35)

可見(jiàn)，該信號(hào)預(yù)測(cè)誤差算子除了在輸出位置的系數(shù)為-1以外，其他系數(shù)與信號(hào)預(yù)測(cè)算子的系數(shù)完全相同。地震信號(hào)的空間可預(yù)測(cè)性主要來(lái)源于地震同相軸的空間展布特征(即構(gòu)造特征)，信號(hào)預(yù)測(cè)誤差算子直接反映了地震反射結(jié)構(gòu)特征。需要注意的是，預(yù)測(cè)誤差算子作用于地震數(shù)據(jù)的結(jié)果不是有效信號(hào)，而是預(yù)測(cè)誤差。

現(xiàn)有的縱向約束吸收補(bǔ)償并未考慮到地震道橫向之間的關(guān)系。為了彌補(bǔ)這一缺陷，將地震信號(hào)的反射結(jié)構(gòu)特征(即信號(hào)預(yù)測(cè)誤差算子)也作為一種約束引入到吸收補(bǔ)償算法?；诜囱莸娜旨s束吸收補(bǔ)償方法具體計(jì)算過(guò)程如下。

(1)通過(guò)解最小二乘的問(wèn)題

(36)

(2)由預(yù)測(cè)算子可得到預(yù)測(cè)誤差算子

P=(a-M，…，a-2，a-1，-1，a1，a2，…，aM)

(37)

(3)通過(guò)反演

(38)

計(jì)算去噪地震數(shù)據(jù)。式中μ2為權(quán)衡參數(shù)。通過(guò)多次迭代求解式(38)，可得頻率域去噪數(shù)據(jù)。

(39)

式中μx為橫向約束因子。新的目標(biāo)函數(shù)的解可表示為

(40)

本文提出的基于反演的全局約束多道吸收補(bǔ)償方法的橫向約束算子由原始的地震剖面求得，包含了原始數(shù)據(jù)的反射結(jié)構(gòu)信息，能夠較明顯地提高吸收補(bǔ)償結(jié)果的穩(wěn)定性和抗噪性，同時(shí)還保持了同相軸的橫向連續(xù)性。

4 應(yīng)用

4.1 Marmousi模型試算

為驗(yàn)證本文方法相對(duì)于縱向約束吸收補(bǔ)償方法的優(yōu)越性，利用Marmousi模型進(jìn)行測(cè)試。

圖1a是原始反射系數(shù)與主頻為30Hz的Ri-cker子波褶積得到的Marmousi模型模擬數(shù)據(jù)，可以看出，Marmousi模型構(gòu)造復(fù)雜，既有平滑的同相軸，也有大型的背斜構(gòu)造和一些小型斷裂。圖1b是Marmousi模型的衰減數(shù)據(jù)，其中Q設(shè)為50。圖1c和圖1d是信噪比(SNR)分別為5和20的加噪衰減數(shù)據(jù)，可以看出，由于地層吸收的影響，地震資料的分辨率降低，原有的構(gòu)造反射特征變模糊。另外，隨機(jī)噪聲污染了地震反射信號(hào)。

分別采用傳統(tǒng)的縱向約束方法和本文的全局約束方法對(duì)圖1b所示的無(wú)噪衰減數(shù)據(jù)進(jìn)行補(bǔ)償，結(jié)果如圖2所示。其中縱向約束吸收補(bǔ)償過(guò)程中μt=0.1，本文方法中μt=0.1，μx=1.0。比較圖2a與圖2b可以看出，在無(wú)噪情況下，縱向約束吸收補(bǔ)償方法和本文方法均能很好地恢復(fù)地震波能量，提高地震資料的分辨率。對(duì)比圖中紅色箭頭處可知，本文方法的結(jié)果橫向更加連續(xù)、平滑，構(gòu)造結(jié)構(gòu)更清晰，弱反射得到了較好的恢復(fù)，證明了本文方法的效果更好。圖3為無(wú)噪衰減數(shù)據(jù)兩種方法補(bǔ)償前、后單道頻譜對(duì)比，可以看出，兩種方法都有效補(bǔ)償了中高頻的振幅衰減，并且拓寬了原始數(shù)據(jù)的地震頻帶，但本文方法補(bǔ)償結(jié)果的頻帶略寬，包含更豐富的高頻有效信息。

圖1 Marmousi模型數(shù)據(jù)(a)原始模擬數(shù)據(jù)； (b)衰減模擬數(shù)據(jù)； (c)SNR=5的加噪衰減數(shù)據(jù)； (d)SNR=20的加噪衰減數(shù)據(jù)

圖2 無(wú)噪衰減數(shù)據(jù)兩種方法補(bǔ)償結(jié)果(a)縱向約束補(bǔ)償； (b)本文方法

對(duì)于含噪衰減數(shù)據(jù)，縱向約束吸收補(bǔ)償結(jié)果雖然恢復(fù)了地震波能量，提高了地震資料的分辨率，但是吸收補(bǔ)償結(jié)果的信噪比較低，地震同相軸的橫向連續(xù)性較差(圖4a)。本文方法的補(bǔ)償結(jié)果不僅能夠提高地震資料的分辨率，使得構(gòu)造反射結(jié)構(gòu)更清晰，同時(shí)也能保證補(bǔ)償結(jié)果的信噪比，橫向上更加連續(xù)(圖4b)，如圖中紅色箭頭所示，本文方法恢復(fù)出了更多的細(xì)節(jié)，表明本文方法比縱向約束吸收補(bǔ)償方法具有更強(qiáng)的穩(wěn)定性和抗噪性。

圖5為SNR=5含噪衰減數(shù)據(jù)利用兩種方法補(bǔ)償結(jié)果的單道對(duì)比，可以看出，本文方法補(bǔ)償后與原始數(shù)據(jù)更相似，無(wú)論振幅還是相位，縱向補(bǔ)償效果均較好(紅色箭頭所示)。圖6為圖5結(jié)果的頻譜對(duì)比，可以看出，在含噪情況下，兩種方法均能補(bǔ)償中高頻的振幅衰減，并且拓寬了補(bǔ)償前數(shù)據(jù)的地震頻帶。此外，縱向約束補(bǔ)償結(jié)果的高頻部分略微強(qiáng)于本文方法補(bǔ)償結(jié)果，有可能是前者對(duì)高頻噪聲的過(guò)度放大所致。

圖3 無(wú)噪數(shù)據(jù)兩種方法補(bǔ)償前、后單道頻譜對(duì)比(a)第100道； (b)第103道； (c)第123道

圖4 SNR=5(左)和SNR=20(右)含噪數(shù)據(jù)兩種方法補(bǔ)償結(jié)果對(duì)比(a)縱向約束方法； (b)本文方法

以均方根振幅為衡量地震剖面橫向連續(xù)性，圖7為SNR=5含噪衰減數(shù)據(jù)兩種方法補(bǔ)償結(jié)果的橫向連續(xù)性對(duì)比，可見(jiàn)本文方法補(bǔ)償結(jié)果橫向上與補(bǔ)償前數(shù)據(jù)更加接近，橫向上抖動(dòng)較小(紅色箭頭所示)，體現(xiàn)了本文方法在保持剖面橫向連續(xù)性和穩(wěn)定性方面的優(yōu)勢(shì)。

圖5 SNR=5含噪數(shù)據(jù)兩種方法補(bǔ)償結(jié)果單道對(duì)比(a)第100道； (b)第103道； (c)第123道

圖6 SNR=5含噪數(shù)據(jù)兩種方法補(bǔ)償結(jié)果單道頻譜對(duì)比(a)第100道； (b)第103道； (c)第123道

圖7 SNR=5含噪數(shù)據(jù)兩種方法補(bǔ)償結(jié)果的橫向連續(xù)性對(duì)比

4.2 實(shí)際數(shù)據(jù)應(yīng)用

圖8a為二維實(shí)際地震疊加剖面，采用縱向約束吸收補(bǔ)償方法和本文方法的處理結(jié)果分別如圖8b、圖8c所示。可以看出，縱向約束補(bǔ)償方法雖然可以提高地震資料的分辨率，但補(bǔ)償后地震剖面中的噪聲被放大，信噪比較低；經(jīng)本文方法處理后，地震資料的能量衰減得到有效補(bǔ)償，地震分辨率得到一定提升，而且補(bǔ)償后地震記錄信噪比較高，噪聲得到有效壓制，地震剖面的同相軸更連續(xù)。由此可見(jiàn)，本文提出的基于反演的全局約束吸收補(bǔ)償方法比傳統(tǒng)的縱向約束吸收補(bǔ)償方法具有更強(qiáng)的穩(wěn)定性和抗噪性，其結(jié)果更利于后續(xù)的解釋和反演。

圖9為圖8數(shù)據(jù)的局部放大對(duì)比，可以看出，相較于原始剖面(圖9a)，本文方法補(bǔ)償剖面(圖9c)同相軸增多、更加清晰，提高了分辨率和信噪比，構(gòu)造更加清楚可辨；相較于縱向約束方法補(bǔ)償結(jié)果(圖9b)，本文方法補(bǔ)償結(jié)果信噪比較高，橫向更連續(xù)。

圖10為實(shí)際地震剖面兩種方法補(bǔ)償后的單道頻譜對(duì)比，可以看出，兩種方法均能有效補(bǔ)償?shù)卣饠?shù)據(jù)中高頻振幅衰減，一定程度上拓寬了數(shù)據(jù)的地震頻帶，其中縱向約束方法處理結(jié)果高頻部分相對(duì)較強(qiáng)，是高頻噪聲放大效應(yīng)的體現(xiàn)。

圖8 實(shí)際地震剖面兩種方法補(bǔ)償結(jié)果對(duì)比(a)原始； (b)縱向約束方法； (c)本文方法

圖9 實(shí)際地震剖面兩種方法補(bǔ)償結(jié)果放大對(duì)比(a)原始剖面； (b)縱向約束補(bǔ)償剖面； (c)全局約束補(bǔ)償剖面

圖10 實(shí)際地震剖面兩種方法補(bǔ)償結(jié)果的單道頻譜對(duì)比(a)第100道； (b)第200道； (c)第300道

5 結(jié)論

(1)本文將傳統(tǒng)的單道吸收補(bǔ)償方法拓展到了多道吸收補(bǔ)償，由地震信號(hào)的可預(yù)測(cè)表達(dá)計(jì)算信號(hào)預(yù)測(cè)誤差算子，并將其作為橫向約束算子引入多道反演框架中，提出了一種基于反演的全局約束多道吸收補(bǔ)償方法。

(2)本文方法的目標(biāo)函數(shù)表達(dá)式包含了縱向和橫向兩個(gè)方向上的約束信息，既考慮了縱向上的地震道稀疏分布信息，又考慮了橫向上的地震道間的連續(xù)性。

(3)模型數(shù)據(jù)和實(shí)際數(shù)據(jù)補(bǔ)償結(jié)果表明，本文提出的基于反演的全局約束多道吸收補(bǔ)償方法具有更強(qiáng)的抗噪性和更高的穩(wěn)定性，可有效減小補(bǔ)償對(duì)于高頻噪聲的放大效應(yīng)，提高了地震記錄分辨率和信噪比，同相軸橫向上更連續(xù)。

(4)由于基于反演的全局約束多道吸收補(bǔ)償方法屬于多道反演算法，運(yùn)算量大，計(jì)算效率較低，以后可從數(shù)值解法和矩陣運(yùn)算等方面進(jìn)行優(yōu)化以提高計(jì)算效率。