孫雷鳴， 曾維輝， 方中于

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

0 引言

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

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

1 譜反演基本原理

1.1 譜反演目標函數的建立

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

圖1 反射系數奇偶分解示意圖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]對奇反射系數對的模擬研究表明，當厚度逐漸減小時，振幅先是增大(調諧)然后減?。欢鴮ε挤瓷湎禂祵?，則正好相反，隨著厚度減小振幅先減后增。因此在實際的大多數情況中反射系數都包含著這二個奇偶反射系數，而地震響應則與這二個反射系數對的大小有關。

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

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

(2)

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

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

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

(3)

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

(4)

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

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

(5)

同理，對于分析點位于時窗中點的多個地層模型，可推導頻率域反射系數表達式如下

jro(i)sin(πfTi)]

(6)

式中N表示奇偶反射系數對的個數；Ti表示為第i層反射系數位置到分析點的2倍時間距離。

圖2 兩個反射系數地層模型Fig．2 Two-layer reflectivity modle

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

(7)

(8)

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

1.2 譜反演目標函數線性化

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

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

(9)

式中f1=fL；fM=fL；M為頻率采樣個數，一般為等間隔采樣，采樣率與傅氏變換時窗長度有關。

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

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

2 線性譜反演存在的問題

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

(10)

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

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

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

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

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

(11)

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

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

3 線性譜反演約束方法

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

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

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

(12)

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

Ar=b

(13)

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

(14)

(ATA+Q)r=ATb

(15)

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

(16)

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

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

(17)

式中μ表示阻尼因子。

4 理論模型試驗分析

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

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

5 實際資料應用

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

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

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

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

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

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

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

6 結論

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

參考文獻：

[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].石油物探,2011,50(2):124-158.

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

[20] 王宇，韓立國，周家雄.L1-L2范數聯合約束稀釋脈沖反演的應用[J].地球科學——中國地質大學學報，2009,34(5):836-840.