劉百紅 李建華 鄭四連

(①中國石化石油物探技術(shù)研究院，江蘇南京 211103； ②中國石油東方地球物理公司研究院地質(zhì)研究中心，河北涿州 072751)

1 問題的提出

在地震勘探中，疊后地震道的形成可用褶積模型描述

S(t)=r(t)*W(t)+n(t)

(1)

式中：S(t)表示地震道；r(t)表示反射系數(shù)序列；W(t)表示子波；n(t)表示噪聲；t為時間； “*”表示褶積。當(dāng)已獲得疊后數(shù)據(jù)體時，往往希望由式(1)得到與數(shù)據(jù)相應(yīng)的地下介質(zhì)的參數(shù)，即反射系數(shù)序列，進而得到波阻抗，這就是所謂的疊后波阻抗反演[1，2]。這個過程中通常都假設(shè)反射系數(shù)序列是稀疏的，并通過井震標(biāo)定確定子波，然后利用某種最優(yōu)化算法求解設(shè)定好的目標(biāo)函數(shù)，獲得反射系數(shù)序列。

根據(jù)反射系數(shù)序列是稀疏的這一假設(shè)，可將反演問題表述為如下形式的基追蹤問題

(2)

式中：d0表示觀測地震道；‖·‖0表示L0范數(shù)。

由于求L0范數(shù)的極值問題很困難，且它對噪聲十分敏感，因此在實際計算中常將式(2)近似表示為

(3)

式中‖·‖1表示L1范數(shù)?？勺C明在一定程度上式(3)的解就是式(2)的解，且式(3)比式(2)容易求解。其中，最常見的式(3)解法就是匹配追蹤算法。

最初的匹配追蹤算法由Mallat等[3]提出，是用于信號分解，隨后出現(xiàn)了精度更高的正交匹配追蹤，并廣泛應(yīng)用于地球物理領(lǐng)域[4-10]。同時也出現(xiàn)了基于匹配追蹤算法的反演，如周東勇等[11]、Wen等[12]提出基于雙極子分解匹配追蹤算法的反演方法，劉曉晶等[13]提出基于正交匹配追蹤算法的反演方法，Zhang等[14-18]和Deborah等[19]提出基于雙極子分解基追蹤算法的反演方法，印興耀等[20]則在雙極子分解的基礎(chǔ)上用梯度投影算法求解反演問題，并開展實際應(yīng)用[21，22]。與此同時，曹靜杰[23]從貝葉斯稀疏反演理論出發(fā)，建立了非凸Lp范數(shù)正則化的盲反褶積模型，并利用最小二乘法進行反演； 梁東輝等[24]在對反射系數(shù)進行L0范數(shù)約束的基礎(chǔ)上，利用迭代硬閾值法進行反射系數(shù)反演。

實際上，在稀疏假設(shè)下，對于地球物理反演問題，其表現(xiàn)形式應(yīng)該是

(4)

式中：k表示稀疏度；p可等于1，也可等于2。當(dāng)p=1時，式(4)轉(zhuǎn)化為線性規(guī)劃問題，即可用線性規(guī)劃算法求解； 當(dāng)p=2時，式(4)轉(zhuǎn)化為二次規(guī)劃問題，即可用通用二次規(guī)劃算法求解。本文提出一種求解方法，依然基于式(2)的形式，即求L0范數(shù)的極值問題，但并不直接計算‖r‖0，而是將‖r‖0用一個平滑函數(shù)來近似。

2 基本原理

‖r‖0實際上定義為反射系數(shù)序列r中非零值的個數(shù)。也就是說，如果定義如下函數(shù)

(5)

那么‖r‖0就可表示為

(6)

其思路即是將式(6)中的階躍函數(shù)s(ri)用連續(xù)的可微函數(shù)來近似。這里選用零均值高斯函數(shù)

(7)

則有

(8)

或者

(9)

則有

于是

‖r‖0≈n-fσ(r)

其中參數(shù)σ為控制精度。σ越小，近似程度越高，σ越大，則函數(shù)越平滑。由于n是一個定值，因此求‖r‖0的極小就轉(zhuǎn)化為求fσ(r)的極大值。于是式(2)的基追蹤問題就可表示為

(10)

該式可用許多通用的基于導(dǎo)數(shù)的優(yōu)化算法[25](如最速下降法)求解。不同之處是在這個優(yōu)化問題中要選擇σ。由于|r|<1，根據(jù)式(9)，可選擇一系列σ，其范圍為0≤σ≤2，如[σ1,σ2,…,σK]=[2.0,1.9,…,0.1]。然后進行如下迭代計算：

一、固定σ=σ1，選擇最小二乘解r0=(WTW+λI)-1·WTd0作為初始解，用最速下降法進行迭代求解：r←r+μ·fσ(r)， 并將其投影到可行域Wr=d0中， 投影方式為r←r-WT·(WWT)-1·(Wr-d0)，從而獲得第一次迭代后的最優(yōu)解r1；

二、令σ=σi，用上一次迭代得到的最優(yōu)解ri作為初始解，用最速下降法進行下一次迭代，并將其投影到可行域，獲得新的最優(yōu)解ri；

三、令i=i+1，重復(fù)第二步；

四、最終獲得最優(yōu)解r=rK。

該過程中，無論是初始化還是在最速下降法中均需要選擇一個參數(shù)λ，但是其意義不同。在初始化時，最小二乘解中λ的意義是使矩陣求逆穩(wěn)定，因此這個時候的λ取值不宜太大。而在最速下降法中，λ可通過線性搜索來確定。為了簡化計算過程也可指定一個值，但是這個定值需試驗確定。本文通過模型試驗，取λ為0.01，并且在實際資料計算時，將實際地震資料進行了歸一化處理，沒有改變λ的值。

對于σ的取值范圍可以根據(jù)實際選擇。在實際地震勘探中由于|r|<1，不妨分別取r=0.5和r=0.05，根據(jù)式(7)可得到近似函數(shù)的變化趨勢，如圖1所示。由圖可見，即使是反射系數(shù)較大的情況，例如r=0.5，取σ最大為2即可使近似函數(shù)接近1；而對于反射系數(shù)較小的情況，例如r=0.05,σ最大為0.1即可使近似函數(shù)接近1。似乎只要取較大的σ(如2或4)，即可很好地近似所有情況，但是實際上其分辨率卻不高，即不能區(qū)分反射系數(shù)大小。因此，在實際應(yīng)用時，將地震資料歸一化以后，在第一層循環(huán)中，從大到小、依次按照一定間隔取一系列σ。只要設(shè)定最大值、最小值及間隔步長即可。

圖1 近似函數(shù)隨σ的變化趨勢

對每一個固定的σ，接下來就是用最優(yōu)化算法求解式(10)。由于此時的目標(biāo)函數(shù)是連續(xù)可微函數(shù)，因此，其優(yōu)化算法就可有很多選擇。本文用最速下降法，其終止條件可用最大迭代次數(shù)、目標(biāo)函數(shù)的迭代終止誤差或者目標(biāo)函數(shù)變量的迭代終止誤差進行控制。

3 模型試驗

眾所周知，如果式(1)中的噪聲項為0，子波準(zhǔn)確且為已知，則可容易且準(zhǔn)確地反演反射系數(shù)。但實際情況是噪聲項未知，子波雖然能提取但通常不十分準(zhǔn)確，此時難以準(zhǔn)確地反演反射系數(shù)，這對反演算法提出了嚴(yán)峻的考驗。為此，本文首先從測井資料獲得了準(zhǔn)確的反射系數(shù)系列； 再用30Hz零相位Ricker子波與之褶積得到模擬地震道； 然后在模擬地震道中加入隨機噪聲作為觀測數(shù)據(jù)(圖2b，信噪比為233)，子波分別用準(zhǔn)確的30Hz和28Hz零相位Ricker子波進行反演，反演結(jié)果如圖3。

從圖3a可見，本文提出的反演方法不僅可獲得準(zhǔn)確的絕對波阻抗值，還可獲得較可靠的波阻抗變化趨勢及界面位置。同時，弱隨機噪聲對本文反演方法的影響并不是很大，且實際資料中的隨機噪聲一般也較小，尤其是疊后或者局部疊加數(shù)據(jù)中的隨機噪聲較小。對比圖3a與圖3b可以看出，子波準(zhǔn)確與否對反演結(jié)果的影響很大。雖然反演結(jié)果也能刻畫波阻抗變化趨勢以及界面形態(tài)，但界面位置及波阻抗值大小都偏離了真實值。這與目前常用的M商用反演軟件的情況類似，即子波或者標(biāo)定的準(zhǔn)確與否決定了最終反演結(jié)果的可靠性。

用時間域波阻抗楔形模型(圖4a)進行試驗。時間采樣間隔是1ms，共1500道。圖4b是用35Hz零相位Ricker子波進行褶積得到的模擬數(shù)據(jù)。以此模擬數(shù)據(jù)作為輸入，用本文提出的反演方法進行反演，反演時的子波采用準(zhǔn)確的子波，即35Hz零相位Ricker子波。圖4c是反演得到的波阻抗。對比圖4a與圖4c可見，當(dāng)子波準(zhǔn)確時，本文的反演方法能十分準(zhǔn)確地獲得反射界面的真實位置和波阻抗真實值，即使兩個反射界面距離非常近，即所謂的薄層情況下，也能準(zhǔn)確刻畫反射界面的真實位置和波阻抗真實值。

圖2 無噪合成道(a)和含噪合成道(b，信噪比為233)

圖3 含噪、無噪合成道分別與準(zhǔn)確(a)、不準(zhǔn)確(b)子波的反演結(jié)果與真實值的對比

圖4 楔形模型真實波阻抗剖面(a)及其合成地震數(shù)據(jù)(b)和反演波阻抗剖面(c)

4 實際地震資料處理

對圖5實際剖面L1進行處理，圖6和圖7是反演后的反射系數(shù)剖面及相對波阻抗剖面。由于反演所用子波是30Hz零相位Ricker子波，而不是經(jīng)過標(biāo)定提取出來的，因此子波未必準(zhǔn)確； 另外因未經(jīng)標(biāo)定，故波阻抗值也只能是相對的。

利用M商用軟件也對圖5地震資料進行了相對波阻抗(圖8)計算。對比圖8和圖7可見，本反演方法獲得的反射系數(shù)剖面和波阻抗剖面能夠保持反射系數(shù)的稀疏性和橫向連續(xù)性，從縱向上看，波阻抗的變化趨勢以及反射界面還是能夠得到清晰的反映，尤其在縱向上，其分辨率并不比商業(yè)軟件反演的結(jié)果低。但是本文的反演方法本質(zhì)上是稀疏脈沖波阻抗反演的一種實現(xiàn)，且在實施過程中是按道進行的，沒有在橫向上進行約束。即本文的反演僅僅利用了地震數(shù)據(jù)，而沒有利用包含測井、地質(zhì)等信息的模型，即沒有利用地震信號帶限外的高低頻信息，尤其是低頻信息，因此本文的反演方法得到的結(jié)果本質(zhì)上是帶限的相對波阻抗。同時，由于在具體實現(xiàn)時是根據(jù)式(4)按道進行的反演，沒有規(guī)則化項，因此反演結(jié)果(波阻抗)仍存在“掛面條”現(xiàn)象(圖7)。

圖9是實際地震剖面L2，采樣間隔為1ms(實際顯示0.30～1.80s)。用本文方法反演后加入低頻趨勢，獲得波阻抗剖面(圖10)； 再將井旁道反演結(jié)果與井口計算的波阻抗(圖11)進行對比。從圖11可見，加入低頻趨勢以后，反演結(jié)果在總體趨勢上與井口的波阻抗吻合很好，既可反映反射界面，也可獲得絕對波阻抗，有利于儲層描述。另外，從圖10可見，反演結(jié)果從縱向上反映出波阻抗變化趨勢，但橫向上不平滑。為了進一步提高反演質(zhì)量，在做好標(biāo)定與子波提取的基礎(chǔ)上，還需在反演過程中加入低頻模型或者橫向約束。

圖5 實際過井地震剖面L1

圖6 本文方法反演的L1反射系數(shù)剖面

圖7 本文方法反演的L1相對波阻抗剖面

圖8 利用M商用軟件反演的L1相對波阻抗剖面

圖9 實際地震剖面L2

圖10 反演后加入低頻趨勢的L2波阻抗剖面

圖11 反演的井旁波阻抗(藍色)與井口波阻抗(綠色)的對比

5 結(jié)論

針對地下地層反射系數(shù)是稀疏的這一假設(shè)，本文將地震波阻抗反演視為L0范數(shù)約束下的最優(yōu)化問題，并將L0范數(shù)用一個平滑函數(shù)近似。然后將L0范數(shù)約束下的基追蹤問題轉(zhuǎn)化為無約束最優(yōu)化問題，并利用基于導(dǎo)數(shù)的最優(yōu)化方法求解無約束優(yōu)化問題，求得反射系數(shù)，進而計算得到波阻抗。由于目標(biāo)函數(shù)中包含了L0范數(shù)約束項，因此反演的反射系數(shù)是稀疏的。同時將L0范數(shù)用一個平滑函數(shù)近似，就可利用基于導(dǎo)數(shù)的優(yōu)化方法，從而使得反演快速收斂。

模型試驗表明，當(dāng)子波準(zhǔn)確時，即使有噪聲，該方法也可十分準(zhǔn)確地獲得波阻抗值和反射界面。由于實際地震資料是帶限的，因此實際資料計算結(jié)果是相對波阻抗，但波阻抗變化趨勢及反射界面可清晰刻畫，并且在加入低頻信息以后波阻抗真值可得到恢復(fù)。為了進一步提升反演效果，還需要在具體實現(xiàn)時加入橫向約束。