地震薄層反射系數(shù)譜反演算法研究及應(yīng)用

孫雷鳴， 曾維輝， 方中于

(中海油能源發(fā)展股份有限公司鉆采工程研究院 地球物理研究所， 湛江 524057)

0 引言

目前，提高地震資料分辨率一直是人們討論的熱門話題，地球物理學(xué)家也一直在探索如何拓寬地震資料的頻帶寬度以獲得更高的分辨率。眾所周知，地震分辨率是從地震資料中提取信息的關(guān)鍵，根據(jù)Widess理論[1]，在地震上時(shí)間厚度小于八分之一波長的薄層是無法分辨的，如果存在噪音或者是隨著子波在地下傳播范圍增大，這個(gè)結(jié)論將變?yōu)樗姆种徊ㄩL[2]。如果地震薄層是潛在的儲層或者是重要的流動(dòng)單元，受限于現(xiàn)有的地震分辨能力，薄層的識別及其厚度的計(jì)算將變得非常困難。

譜反演是一種采用頻譜分解技術(shù),來提高小于調(diào)諧厚度的薄層成像的處理技術(shù)。最初是由John P. Castagna等人[2]提出的，在其發(fā)表的論文中提出了反射系數(shù)序列奇偶分解可以提高分辨率的理論，給出了譜反演的目標(biāo)函數(shù)，指出譜反演能夠分辨薄層反射系數(shù)的位置、極性和大小；John P. Castagna等[2-8,11-14]利用頻譜分解獲得的局部頻譜資料信息來進(jìn)行薄層反射系數(shù)譜反演，指出該方法具有不需要任何先驗(yàn)?zāi)Ｐ停恍枰魏畏瓷湎禂?shù)的數(shù)學(xué)假設(shè)，不需要任何層位約束，也不需要任何測井資料強(qiáng)制約束等優(yōu)點(diǎn)，可用來分辨低于調(diào)諧厚度的地震薄層；John P. Castagna等人提出的譜反演目標(biāo)函數(shù)為非線性方程，同時(shí)求解反射系數(shù)位置及大小，存在多解性和不確定性，直接求解難度較大。作者將譜反演目標(biāo)函數(shù)轉(zhuǎn)化線性方程組[17]并求解，發(fā)現(xiàn)要想求得精確解，需要0到Niquist頻率的所有頻率成分，而地震資料往往是帶限的，袁三一[10]在文獻(xiàn)中也提到了這一觀點(diǎn)。如何解決這個(gè)問題，作者對此進(jìn)行了探討。

1 譜反演基本原理

1.1 譜反演目標(biāo)函數(shù)的建立

譜反演根據(jù)地震記錄和地震子波的頻譜信息，從地震記錄中去除地震子波的影響，從而得到反射系數(shù)序列。根據(jù)函數(shù)奇偶分解理論[2]，對于任意反射系數(shù)r(t)，都可以將其分解成奇反射系數(shù)對ro(t)和偶反射系數(shù)對re(t)(圖1)。

圖1 反射系數(shù)奇偶分解示意圖Fig．1 Odd Even Decomposition of refection coefficients

r(t)=ro(t)+re(t)

(1)

其中ro(t)=[r(t)-r(t)]/2；re(t)=[r(t)+r(t)]/2。

John P. Castagna[2-8]對奇反射系數(shù)對的模擬研究表明，當(dāng)厚度逐漸減小時(shí)，振幅先是增大(調(diào)諧)然后減小；而對偶反射系數(shù)對，則正好相反，隨著厚度減小振幅先減后增。因此在實(shí)際的大多數(shù)情況中反射系數(shù)都包含著這二個(gè)奇偶反射系數(shù)，而地震響應(yīng)則與這二個(gè)反射系數(shù)對的大小有關(guān)。

假定水平疊加后的零炮檢距自激自收地震記錄可以用褶積模型來表示，即s(t)=w(t)*r(t)，在頻率域可表示為

S(f)=W(f)×R(f)

(2)

其中 “*”表示褶積運(yùn)算;“×”表示乘積運(yùn)算;S(f)為地震信號頻譜;W(f)為地震子波頻譜;R(f)為反射系數(shù)頻譜。

下面簡要給出譜反演的目標(biāo)函數(shù)的推導(dǎo)過程。圖2為一個(gè)兩個(gè)反射系數(shù)薄層模型，反射系數(shù)時(shí)間位置為t1和t2，大小分別為r1和r2，時(shí)間厚度為T，反射系數(shù)可表示為

r(t)=r1δ(t-t1)+r2δ(t-t1-T)

(3)

若將分析點(diǎn)設(shè)置于地層中點(diǎn)，則式(3)又可表示為

(4)

對式(4)進(jìn)行傅里葉變換，并根據(jù)反射系數(shù)奇偶分解理論、脈沖函數(shù)δ(t)性質(zhì)以及歐拉公式可得

R(re,ro,T,f)=2recos(πfT)+j2rosin(πfT)

(5)

同理，對于分析點(diǎn)位于時(shí)窗中點(diǎn)的多個(gè)地層模型，可推導(dǎo)頻率域反射系數(shù)表達(dá)式如下

jro(i)sin(πfTi)]

(6)

式中N表示奇偶反射系數(shù)對的個(gè)數(shù)；Ti表示為第i層反射系數(shù)位置到分析點(diǎn)的2倍時(shí)間距離。

圖2 兩個(gè)反射系數(shù)地層模型Fig．2 Two-layer reflectivity modle

由于分析點(diǎn)設(shè)置在時(shí)窗中點(diǎn)，因此需要對地震信號進(jìn)行時(shí)移，假如時(shí)窗長度為L，則時(shí)移量Δω=L/2，則式(6)等價(jià)于(7)

(7)

(8)

式中 Re[…]和lm[…]分別表示頻譜的實(shí)部和虛部；fi表示采樣頻率；α和β分別表示反射系數(shù)的偶反射系數(shù)對和奇反射系數(shù)對所占的比例；fL和fH分別表示為最低頻率和最高頻率。

1.2 譜反演目標(biāo)函數(shù)線性化

從式(8)可以看出，需求解的未知向量包括re、ro、T，即反射系數(shù)的大小及其位置，且這些未知量之間相互影響，存在多解性和不確定性，精確求解非常依賴于算法的全局尋優(yōu)性能[9,10]。如果將方程轉(zhuǎn)化一個(gè)線性方程組，問題可能會變得相對簡單。

假如在反射系數(shù)位置已知的情況下，即向量T已知，待求解未知量僅為re和ro，設(shè)α=β=1，根據(jù)式(6)和式(7)，可得到以下線性方程組，即

(9)

式中f1=fL；fM=fL；M為頻率采樣個(gè)數(shù)，一般為等間隔采樣，采樣率與傅氏變換時(shí)窗長度有關(guān)。

式(9)是在假設(shè)反射系數(shù)位置T=(T1,T2,…,TN)已知的情況下得到的結(jié)果，實(shí)際上，它是未知的。若假設(shè)地震道集每個(gè)采樣點(diǎn)都看做是一個(gè)反射系數(shù)，其所在的時(shí)間位置便是一個(gè)已知值，并且假設(shè)頻率分析時(shí)窗采樣點(diǎn)個(gè)數(shù)L為偶數(shù)，采樣率為Δt，則在式(9)中，N=L/2，Ti=(2i+1)Δt(i=1,2,…,N)，這樣就避免了反射系數(shù)位置的求解，只需求解re和ro即可。

通過求解方程組(9)便可得到反射系數(shù)的奇反射系數(shù)和偶反射系數(shù)，然后對其進(jìn)行重構(gòu)便可得到反射系數(shù)的解。求解方法可采用SVD、共軛梯度法、LSCG[2]、最小二乘QR分解[16-17]等方法，本文采用最小二乘QR分解法對其進(jìn)行求解。

2 線性譜反演存在的問題

由公式(2)可知，在不考慮噪音的情況下，反射系數(shù)頻譜可以表示為

(10)

式中f∈[0,fN]；fN表示Nyquist頻率。理論上基于(10)式的方程，旨在得到反射系數(shù)的反演，需要地震數(shù)據(jù)0到Nyquist頻率的所有頻率成分[10]。在假設(shè)每個(gè)時(shí)間采樣點(diǎn)都存在一個(gè)反射系數(shù)值的情況下，所求未知量的個(gè)數(shù)遠(yuǎn)遠(yuǎn)超過真實(shí)反射系數(shù)的個(gè)數(shù)，這種要求顯得更為苛刻。

圖3是用30 Hz雷克子波合成地震記錄在兩不同頻率范圍內(nèi)反演的結(jié)果，反射系數(shù)如圖3(b)所示，從圖3可以看出，該反射系數(shù)較為復(fù)雜，含有薄互層。合成地震記錄采樣間隔為 4 ms，圖3(d)是在0到Nyquist頻率(125 Hz)范圍內(nèi)反演的結(jié)果，從圖3(d)可以出，反演效果非常理想，薄層得到了很好地反演；圖3(f)是在 5 Hz 到 95 Hz 頻率反演內(nèi)反演的結(jié)果，從圖3(f)可以看出，反演結(jié)果能量分散，真實(shí)反射系數(shù)位置周圍存在 “抖動(dòng)”，并且隨著反演頻帶寬度的減小，反演效果會越來越差，能量越來越分散，分辨率也越來越低。

圖3 合成地震記錄不同頻率范圍反演的結(jié)果Fig．3 Inversion results of different spectrum of the same synthetic seismogram(a)Ricker子波；(b)反射系數(shù);(c)合成地震記錄；(d)[0,125Hz]反演結(jié)果；(e)[5,95 Hz]反演結(jié)果

在采樣間隔為 4 ms 的情況下，Nyquist頻率為 125 Hz，在f∈[0,125Hz]的情況下反演能取得好的效果，但在采樣間隔為 2 ms 的情況下，Nyquist為 250 Hz，在f∈[0,250Hz]的情況下是否也能取得好的效果呢？采用圖3的合成地震記錄，但采樣間隔為 2 ms。在頻譜范圍f∈[0,250Hz]反演得到的結(jié)果如圖4(b)所示，反演結(jié)果并不理想，存在較大誤差。將真實(shí)反射系數(shù)和反演得到的反射系數(shù)做傅里葉變換得到其頻譜，如圖4(d)和圖4(c)所示，從圖4可以看出，反射系數(shù)的高頻區(qū)域沒有得到正確反演。

事實(shí)上，地震子波是帶限的，其高頻成分的幅值接近于零值，在求解S(f)/W(f)時(shí)存在不確定性。基于式(10)的求解反射系數(shù)的方程組可以由式(11)來計(jì)算：[10]

(11)

其中相應(yīng)振幅值接近于零的采樣頻率f，f∈[fhigh,fN]，rand表示一個(gè)隨機(jī)復(fù)數(shù)。這便解釋了采樣間隔為 2 ms 時(shí)，基于所有頻率成分沒有取得好的反演結(jié)果的原因。

實(shí)際上所有的地震資料都是帶限的，不可能利用所有的頻率成分去反演反射系數(shù)，從而必須舍棄部分低頻成分和高頻成分。綜上所述，對于線性譜反演，僅僅利用帶限的地震資料去反演得不到精確解。 作者認(rèn)為，要解決這個(gè)問題，必須假定反射系數(shù)是稀疏的，即地層介質(zhì)波阻抗分界面處的反射系數(shù)為非零的，非地層介質(zhì)波阻抗分界面的反射系數(shù)為零，而線性譜反演方法認(rèn)為地震記錄每個(gè)采樣點(diǎn)都是一個(gè)反射界面，即反射系數(shù)是非稀疏的，因此需要對其進(jìn)行稀疏約束。

3 線性譜反演約束方法

由于地震資料是帶限的，不可能利用所有頻率去進(jìn)行反演，而線性譜反演假定反射系數(shù)在每個(gè)時(shí)間采樣點(diǎn)均存在值，要想取得精確解就需要所有頻率成分。在地震資料帶限的情況下，如何提高線性譜反演的精度，作者對線性反演進(jìn)行約束，成功地解決了這個(gè)問題。該方法分為兩步，先在方程(9)加入Cauchy約束并反演，在反演結(jié)果中篩選若干反射系數(shù)(個(gè)數(shù)應(yīng)不少真實(shí)反射系數(shù)的個(gè)數(shù))，這些反射系數(shù)是稀疏的，時(shí)間位置已知，然后再利用式(9)對這些反射系數(shù)進(jìn)行反演，反演時(shí)加入 L2 模約束。

圖4 合成地震記錄反演的結(jié)果及頻譜對比Fig．4 Inversion results of synthetic seismogram and spectrum comparison(a)反射系數(shù)；(b)反演結(jié)果；(c)反演誤差；(d)真實(shí)反射系數(shù)頻譜；(e)反演得到的頻譜

Cauchy 約束即為稀疏約束，研究表明，Cauchy 約束準(zhǔn)則在反射系數(shù)位置檢測和抗噪等方面較為適用，反射系數(shù)的 Cauchy 分布可表示為[18-19]。

(12)

線性方程式(9)可簡寫為

Ar=b

(13)

式中A為方程式(9)的系數(shù)矩陣；r為奇偶反射系數(shù)序列，即r=(re(1),…,re(N),ro(1),…,ro(N))T，b=(Re[Ψ(f1)],…,Re[Ψ(fM)],lm[Ψ(f1)],…,lm[Ψ(fM)])T。方程式(13)引入 Cauchy 約束準(zhǔn)則，則式(13)可表示為

(14)

(ATA+Q)r=ATb

(15)

式中Q表示約束對角矩陣。

(16)

式(15)為非線性隱式格式，采用最小二乘 QR 分解法對其進(jìn)行迭代求解。通過 Cauchy 約束求得的解并不是非常精確的，其值較真實(shí)值要偏小。為了提高反射系數(shù)反演的精度，從 Cauchy 約束反演得到的稀疏反射系數(shù)中篩選若干反射系數(shù)，這些反射系數(shù)的位置可以代表真實(shí)反射系數(shù)所在的位置，然后再對這些反射系數(shù)進(jìn)行反演，并且在反演過程中加入 L2 模約束，保證算法的穩(wěn)定性，收斂到真實(shí)解。此時(shí)反演的目標(biāo)函數(shù)可表示為

Obj2=‖Ar-b‖2+μ‖r‖

(17)

式中μ表示阻尼因子。

4 理論模型試驗(yàn)分析

圖5模型試驗(yàn)同樣采用圖3的反射系數(shù)及合成地震記錄，時(shí)間采樣間隔為 2 ms。僅對反演添加 Cauchy 約束，在頻率范圍[5,150 Hz]內(nèi)得到的反演結(jié)果如圖5(d)所示，與不加約束的反演結(jié)果(圖5(b))進(jìn)行對比可以看出，添加約束后反射系數(shù)使能量更加聚集到反射界面上，而不加約束的反演在真實(shí)反射系數(shù)周圍存在明顯的“抖動(dòng)”，能量分散。盡管如此，加約束后反演得到的結(jié)果與真實(shí)反射系數(shù)相比還是存在較大的誤差，如圖 5(f) 所示，這是由于 Cauchy 約束本身對反射系數(shù)也進(jìn)行壓制而不可避免的。

圖6是對 Cauchy 約束反演后篩選 80 個(gè)反射系數(shù)并再加 L2 模約束反演求得的結(jié)果，圖6(b)是在頻率范圍[5,130 Hz]內(nèi)反演得到的結(jié)果，從圖6可以看出，得到的反射系數(shù)非常接近于真實(shí)反射系數(shù)，并且誤差非常小；圖6(d)是在頻率范圍[10,100 Hz]內(nèi)反演求得的結(jié)果，從圖6(d)可以看出，減少部分高頻和低頻成分，依然能求得好的反演結(jié)果，誤差也非常小，但是相對于在頻率范圍[5,130 Hz]求得的反演結(jié)果誤差要偏大，這說明參與反演的頻率成分越多，反演效果就越精確。

5 實(shí)際資料應(yīng)用

將本文提出的譜反演方法應(yīng)用到實(shí)際疊后地震數(shù)據(jù)中，圖7(左)是來自某一個(gè)工區(qū)的實(shí)際地震剖面，采樣率是 1 ms，地震道長度為 8 000 ms，圖7(右)是利用譜反演方法反演的反射系數(shù)剖面。反演子波是采用統(tǒng)計(jì)法計(jì)算得到的雷克子波，沒有用到井?dāng)?shù)據(jù)約束，并且由于地震波在傳播過程存在高頻衰減現(xiàn)象，在不同時(shí)間段提取了不同的子波。從圖7可以看出，反演得到的反射系數(shù)剖面與真實(shí)地層相對應(yīng)，地震薄層也得到了很好地反演，如圖7(右)方框區(qū)域，這些薄層用肉眼在地震剖面是無法識別的，同時(shí)反射系數(shù)剖面給出了比原始剖面更多的地層細(xì)節(jié)和斷層細(xì)節(jié)，這對精細(xì)解釋來說是非常有益的。

事實(shí)上譜反演不僅能保留中低頻成分，還能較好地補(bǔ)償高頻成分。圖8 是反演得到反射系數(shù)剖面(方框區(qū))某一道與測井?dāng)?shù)據(jù)進(jìn)行對比，利用聲波測井曲線得到的真實(shí)反射系數(shù)與反演子波進(jìn)行褶積，得到合成地震記錄，地震道用藍(lán)色標(biāo)記，反演得到的反射系數(shù)與反演子波褶積得到的合成地震記錄用紅色標(biāo)記，實(shí)際地震道用黑色標(biāo)記。從圖8可以看出，兩者基本上是一致的，這說明反演結(jié)果真實(shí)反映了實(shí)際地層在地震已有的有限帶寬內(nèi)的信息。若將利用聲波測井曲線得到的真實(shí)反射系數(shù)與一帶寬較大的子波進(jìn)行褶積，如圖9所示，得到合成記錄用藍(lán)色標(biāo)記，反演得到的反射系數(shù)與這個(gè)子波褶積的得到合成地震記錄用紅色標(biāo)記。從圖9可以看出，地震比較微弱的高頻信息通過譜反演得到了很好補(bǔ)償，并且與測井?dāng)?shù)據(jù)對應(yīng)關(guān)系良好。

圖5 合成地震記錄[5,150 Hz]范圍內(nèi)加Cauchy約束前后的反演結(jié)果Fig．5 The comparison of inversion results before and after Cauchy restriction in synthetic seismogram's [5,150 Hz] frequency range (a)反射系數(shù)；(b)不加約束反演結(jié)果；(c)圖(b)反演誤差；(d)加約束反演結(jié)果；(e)圖(d)反演誤差

圖6 不同頻率范圍的反演結(jié)果Fig．6 Inversion results of different frequency rang(a)反射系數(shù)；(b)[5,150 Hz]反演結(jié)果；(c)圖(b)反演誤差；(d)[10,110 Hz]反演結(jié)果；(e)圖(d)反演誤差

圖7 原始地震剖面與和反演得到的反射系數(shù)對比Fig．7 The comparison of raw seismic section and inversed reflection coefficients

圖8 反演得到的反射系數(shù)和實(shí)際子波的褶積與測井合成記錄的對比Fig．8 The synthetic seismogram comparison of inversed reflectivity convolved with actual wavelet and well-log

圖9 反射系數(shù)和寬頻子波的褶積與測井合成記錄的對比Fig．9 The synthetic seismogram comparison of inversed reflectivity convolved with broad spectrum wavelet and well-log

6 結(jié)論

譜反演是一種高分辨率反射系數(shù)成像技術(shù)，能分辨低于調(diào)諧厚度的薄層，并且不需要任何初始模型和先驗(yàn)信息約束。作者對譜反演線性方法進(jìn)行了研究，揭示了假定所有時(shí)間采樣點(diǎn)均存在反射系數(shù)的需要所有頻率成分與地震資料帶限之間的矛盾，直接求解不能求得精確解。針對此問題，作者認(rèn)為只有在稀疏假設(shè)條件下的才能正確求解，通過在線性譜反演過程中添加Cauchy約束求得稀疏反射系數(shù)，然后挑選出這些稀疏反射系數(shù)，通過 L2 模約束再次求解，該方法在模型試算和實(shí)際資料應(yīng)用過程中取得了好的效果。結(jié)果表明，本文提出的方法能有效分辨薄層，并且拓寬地震頻帶提高分辨率。該方法在實(shí)際應(yīng)用僅考慮了地震子波的時(shí)變特征，若綜合考慮子波的空變特征，反演效果會更佳。

參考文獻(xiàn)：

[1] WIDESS M B. How thin is a thin bed? [J]. Geophysics, 1973, 38(6):1176-1180.

[2] PURYEAR C. I., CASTAGNA. J. P. Layer-thickness determination and stratigraphic interpretation using spectral inversion:Theory and application [J].Geophysics, 2008,73(2): 37- 48.

[3] CHOPRA S, CASTAGNA J. P., PORTNIAGUINE. O. Thin-bed reflectivity inversion [J].76th Annual International Meeting, SEG, Expanded Abstracts, 2006, 31(1): 2057-2061.

[4] CASTAGNA J. P. Spectral decomposition and high resolution reflectivity inversion [M]. Presented at the Oklahoma Section Meeting, SEG， 2004.

[5] CHOPRA S, CASTAGNA J P, XU Y, et al. Thin-bed reflectivity inversion and seismic interpretation [J].76th Annual International Meeting, SEG, Expanded Abstracts， 2007: 1923-1927.

[6] PORTNIAGUINE O, CASTAGNA J P. Spectral inversion: Lessons from modeling and Boonesville case study [J]. 75th Annual International Meeting, SEG, Expanded Abstracts, 2005:1638-1641.

[7] CHOPRA S, CASTAGNA J P, XU Y. Relative acoustic impedance application for thin-bed reflectivity inversion [J]. Expanded Abstracts of 79thAnnual Internat SEG Mtg, 2009:3554-3558.

[8] PORTNIAGUINE O, CASTAGNA J P. Inverse spectral decomposition [J], Expanded Abstracts of 74thAnnual Internat SEG Mtg, 2004,:1786-1789.

[9] YUAN S Y, WANG S X. A fast hybrid algorithm for spectral inversion [J]. Expanded Abstracts of 79thAnnual Internat SEG Mtg, 2009: 2447-2451.

[10] YUAN S Y, WANG SX, Xu Y C, et al. Ill-posed analysis for spectral inversion [J]. Expanded Abstracts of 79thAnnual Internat SEG Mtg, 2009:2442-2446.

[11] PURYEAR C I, CASTAGNA J P. An algorithm for calculation of bed thickness and reflection coefficients from amplitude spectrum [J]. Expanded Abstracts of 76th Annual International SEG Meeting, 2006: 1767-1770.

[12] CHOPRA S., CASTAGNA J. P., PORTNIAGUINE O. Seismic resolution and thin-bed reflectivity inversion[J]. CSEG Recorder, 2006:31(1):19-25.

[13] CASTANO K P, OJEDA G. Optimizing thin-layer mapping through Spectral Inversion: Performance of Genetic Algorithms and Simulated Annealing[J]. Expanded Abstracts of 80th Annual International SEG Meeting, 2010:1610-1614.

[14] CHOPRA S., CASTAGNA J. P., XU Y. Thin-bed reflectivity inversion and some applications [J]. First Break, 200905, 27(5): 17-24.

[15] KALLWEIT R S, WOOD L C. The limits of resolution of zero-phase wavelets [J]. Geophysics, 1982, 47(7):1035-1046.

[16] PAIGE C C, SAUNDERS M A .LSQR: an algorithm for sparse linear equations and sparse least squares [J]. ACM Transactions on Mathematical software, 1982,8(1):43-71.

[17] 柴新濤,李振春,韓文功. 基于LSQR算法的譜反演方法研究[J].石油物探, 2012, 51(1): 11-18.

[18] 邸海濱,郭玉倩,劉喜武. 基于稀疏約束貝葉斯估計(jì)的相對波阻抗反演[J].石油物探,2011,50(2):124-158.

[19] 曾銳，劉洪，陳世軍.柯西約束盲反褶積技術(shù)在井間地震的應(yīng)用[J].地球物理學(xué)進(jìn)展，2004,19(1):166-172.

[20] 王宇，韓立國，周家雄.L1-L2范數(shù)聯(lián)合約束稀釋脈沖反演的應(yīng)用[J].地球科學(xué)——中國地質(zhì)大學(xué)學(xué)報(bào)，2009,34(5):836-840.