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

孫雷鳴， 曾維輝， 方中于

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

0 引言

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

譜反演是一種采用頻譜分解技術(shù),來(lái)提高小于調(diào)諧厚度的薄層成像的處理技術(shù)。最初是由John P. Castagna等人[2]提出的，在其發(fā)表的論文中提出了反射系數(shù)序列奇偶分解可以提高分辨率的理論，給出了譜反演的目標(biāo)函數(shù)，指出譜反演能夠分辨薄層反射系數(shù)的位置、極性和大?。籎ohn P. Castagna等[2-8,11-14]利用頻譜分解獲得的局部頻譜資料信息來(lái)進(jìn)行薄層反射系數(shù)譜反演，指出該方法具有不需要任何先驗(yàn)?zāi)Ｐ停恍枰魏畏瓷湎禂?shù)的數(shù)學(xué)假設(shè)，不需要任何層位約束，也不需要任何測(cè)井資料強(qiáng)制約束等優(yōu)點(diǎn)，可用來(lái)分辨低于調(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è)問(wèn)題，作者對(duì)此進(jìn)行了探討。

1 譜反演基本原理

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

譜反演根據(jù)地震記錄和地震子波的頻譜信息，從地震記錄中去除地震子波的影響，從而得到反射系數(shù)序列。根據(jù)函數(shù)奇偶分解理論[2]，對(duì)于任意反射系數(shù)r(t)，都可以將其分解成奇反射系數(shù)對(duì)ro(t)和偶反射系數(shù)對(duì)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]對(duì)奇反射系數(shù)對(duì)的模擬研究表明，當(dāng)厚度逐漸減小時(shí)，振幅先是增大(調(diào)諧)然后減??；而對(duì)偶反射系數(shù)對(duì)，則正好相反，隨著厚度減小振幅先減后增。因此在實(shí)際的大多數(shù)情況中反射系數(shù)都包含著這二個(gè)奇偶反射系數(shù)，而地震響應(yīng)則與這二個(gè)反射系數(shù)對(duì)的大小有關(guān)。

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

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

(2)

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

下面簡(jiǎn)要給出譜反演的目標(biāo)函數(shù)的推導(dǎo)過(guò)程。圖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)

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

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

(5)

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

jro(i)sin(πfTi)]

(6)

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

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

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

(7)

(8)

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

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

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

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

(9)

式中f1=fL；fM=fL；M為頻率采樣個(gè)數(shù)，一般為等間隔采樣，采樣率與傅氏變換時(shí)窗長(zhǎng)度有關(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即可。

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

2 線性譜反演存在的問(wèn)題

由公式(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)超過(guò)真實(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)”，并且隨著反演頻帶寬度的減小，反演效果會(huì)越來(lái)越差，能量越來(lái)越分散，分辨率也越來(lái)越低。

圖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ū)域沒(méi)有得到正確反演。

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

(11)

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

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

3 線性譜反演約束方法

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

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

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

(12)

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

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表示約束對(duì)角矩陣。

(16)

式(15)為非線性隱式格式，采用最小二乘 QR 分解法對(duì)其進(jìn)行迭代求解。通過(guò) Cauchy 約束求得的解并不是非常精確的，其值較真實(shí)值要偏小。為了提高反射系數(shù)反演的精度，從 Cauchy 約束反演得到的稀疏反射系數(shù)中篩選若干反射系數(shù)，這些反射系數(shù)的位置可以代表真實(shí)反射系數(shù)所在的位置，然后再對(duì)這些反射系數(shù)進(jìn)行反演，并且在反演過(guò)程中加入 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。僅對(duì)反演添加 Cauchy 約束，在頻率范圍[5,150 Hz]內(nèi)得到的反演結(jié)果如圖5(d)所示，與不加約束的反演結(jié)果(圖5(b))進(jìn)行對(duì)比可以看出，添加約束后反射系數(shù)使能量更加聚集到反射界面上，而不加約束的反演在真實(shí)反射系數(shù)周圍存在明顯的“抖動(dòng)”，能量分散。盡管如此，加約束后反演得到的結(jié)果與真實(shí)反射系數(shù)相比還是存在較大的誤差，如圖 5(f) 所示，這是由于 Cauchy 約束本身對(duì)反射系數(shù)也進(jìn)行壓制而不可避免的。

圖6是對(duì) 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é)果，誤差也非常小，但是相對(duì)于在頻率范圍[5,130 Hz]求得的反演結(jié)果誤差要偏大，這說(shuō)明參與反演的頻率成分越多，反演效果就越精確。

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

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

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

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

圖9 反射系數(shù)和寬頻子波的褶積與測(cè)井合成記錄的對(duì)比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)信息約束。作者對(duì)譜反演線性方法進(jìn)行了研究，揭示了假定所有時(shí)間采樣點(diǎn)均存在反射系數(shù)的需要所有頻率成分與地震資料帶限之間的矛盾，直接求解不能求得精確解。針對(duì)此問(wèn)題，作者認(rèn)為只有在稀疏假設(shè)條件下的才能正確求解，通過(guò)在線性譜反演過(guò)程中添加Cauchy約束求得稀疏反射系數(shù)，然后挑選出這些稀疏反射系數(shù)，通過(guò) L2 模約束再次求解，該方法在模型試算和實(shí)際資料應(yīng)用過(guò)程中取得了好的效果。結(jié)果表明，本文提出的方法能有效分辨薄層，并且拓寬地震頻帶提高分辨率。該方法在實(shí)際應(yīng)用僅考慮了地震子波的時(shí)變特征，若綜合考慮子波的空變特征，反演效果會(huì)更佳。

參考文獻(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ì)的相對(duì)波阻抗反演[J].石油物探,2011,50(2):124-158.

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

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