基于L1范數(shù)的波阻抗反演

劉百紅 ， 李建華 ， 鄭四連

(1.中國石化 石油物探技術(shù)研究院，南京 211103；2.東方地球物理公司研究院 地質(zhì)研究中心，涿州 072751)

0 引言

目前的油氣勘探中最主要的方法是地震反射波法(反射地震勘探)。反射地震勘探的首要條件是地下介質(zhì)中存在反射界面(波阻抗界面)。而地震反演是一個相反的過程，即由地震資料來獲得地層的波阻抗。由于反演可以將地震數(shù)據(jù)轉(zhuǎn)化為地下地層的波阻抗、密度、速度、孔隙度、滲透率、砂泥巖百分比、壓力等信息，因此反演在目前的儲層預(yù)測中起到越來越重要的作用。反演技術(shù)也從疊后反演走向疊前反演，而實現(xiàn)方法則從直接反演走向基于模型反演。無論哪種反演都是基于褶積模型，而反演的基本思想也主要是數(shù)據(jù)擬合。同時為了減少反演的多解性，在反演的過程往往會加入測井、地震構(gòu)造解釋、地質(zhì)等信息來構(gòu)建初始模型對反演過程進(jìn)行約束。在這個過程除了要求較精確的子波外，還要求初始模型接近真實模型，才能達(dá)到可靠的結(jié)果。雖然基于全局優(yōu)化算法的反演方法對初始模型的依賴程度不高，但是其計算效率還是不如局部優(yōu)化算法。而局部優(yōu)化算法除了對初始模型依賴程度高以外，還要求目標(biāo)函數(shù)可微，因此，其目標(biāo)函數(shù)以及約束都是以L2范數(shù)為主。近年來，基于L1范數(shù)約束的反演開始被廣泛研究。相應(yīng)的算法主要是以匹配追蹤為主。最初的匹配追蹤算法由Mallat等[1]提出，是用于信號分解，隨后出現(xiàn)了精度更高的正交匹配追蹤，并在地球物理領(lǐng)域廣泛應(yīng)用[2-8]。同時也出現(xiàn)了基于匹配追蹤算法的反演，周東勇等[9]提出了基于雙極子分解匹配追蹤算法的反演方法;劉曉晶等[10]提出了基于正交匹配追蹤算法的反演方法；Zhang等[11]提出了基于雙極子分解基追蹤算法的反演方法。匹配追蹤算法需要首先建立完備的字典集，而且是一種貪婪迭代算法，因此其迭代次數(shù)多，并且有可能陷入局部最優(yōu)解。筆者將L1范數(shù)用一個平滑函數(shù)來近似，然后將L1范數(shù)約束下的基追蹤問題轉(zhuǎn)化為無約束最優(yōu)化問題，并利用基于導(dǎo)數(shù)的最優(yōu)化方法求解無約束優(yōu)化問題?；趯?dǎo)數(shù)的最優(yōu)化方法收斂速度快，并且當(dāng)初是解接近真解時可以得到真解的近似解。

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ù)序列，進(jìn)而得到波阻抗，這就是疊后波阻抗反演。在這個過程中通常都是假設(shè)反射系數(shù)序列是稀疏的，并通過井震標(biāo)定確定出子波，然后利用某種最優(yōu)化算法求解設(shè)定好的目標(biāo)函數(shù)，來獲得反射系數(shù)序列。最基本的目標(biāo)函數(shù)是下面的誤差函數(shù)：

E(r)=‖S-d0‖2

(2)

其中:d0表示觀測到的地震道。在實際反演中加入先驗信息時，應(yīng)在目標(biāo)函數(shù)中加入規(guī)則化項，從而使得目標(biāo)函數(shù)為式(3)。

F(r)=‖S-d0‖p+ε‖r‖q

(3)

其中：ε為加權(quán)系數(shù)；p、q為范數(shù)。一般取：p=q=2，則可以用基于導(dǎo)數(shù)的最優(yōu)化算法求解式(3)，獲得最小平方解。最小平方解是一個相對平滑的序列，為了提高反演的分辨力，近年來基于q=1的反演開始被廣泛研究，而q=1與反射系數(shù)序列是稀疏的這一假設(shè)更吻合。

實際上，當(dāng)p=2、q=1時，式(3)僅僅是無約束優(yōu)化的一般形式。考慮到反射系數(shù)的定義：

(4)

應(yīng)該有約束：|r|≤1。

另外，當(dāng)p=2、q=1時，式(3)的無約束優(yōu)化問題等價于如下約束優(yōu)化問題：

(5)

由于式(5)的優(yōu)化問題的解也是如下形式的優(yōu)化問題的解：

(6)

因此實際上常用匹配追蹤算法解式(6)，得到反射系數(shù)。匹配追蹤算法是一種貪婪算法，需要完備的原子空間，而該原子空間則是通過對反射系數(shù)進(jìn)行奇偶分解來建立的。Zhang[11]等在反射系數(shù)奇偶分解的基礎(chǔ)上，利用梯度投影算法解目標(biāo)函數(shù)。筆者將L1范數(shù)用一個平滑函數(shù)來近似；然后將L1范數(shù)約束下的基追蹤問題轉(zhuǎn)化為無約束最優(yōu)化問題，并利用基于導(dǎo)數(shù)的最優(yōu)化方法求解無約束優(yōu)化問題。其基本思想如下：

無論式(5)或者式(6)的最優(yōu)化問題都可以轉(zhuǎn)化為如下形式的優(yōu)化問題：

(7)

(8)

其中，ε是一個很小的數(shù)，例如可以將其設(shè)置為：ε=10-10。則目標(biāo)函數(shù)變?yōu)椋?/p>

(9)

這樣，式(7)中的目標(biāo)函數(shù)就變?yōu)檫B續(xù)可微，從而可以用基于導(dǎo)數(shù)的最優(yōu)化方法求解式(9)。

2 模型試驗

我們知道，如果式(1)中的噪聲項為“0”，子波準(zhǔn)確并且已知，那么反射系數(shù)是可以十分容易而且準(zhǔn)確地反演得到。但是實際情況是噪聲項未知，子波雖然提取出來了但也不是十分準(zhǔn)確，這種情況下，準(zhǔn)確地反演出反射系數(shù)就不那么容易了，這對反演算法提出了嚴(yán)峻的考驗。為此，首先從測井資料獲得了準(zhǔn)確的反射系數(shù)系列，用35 Hz零相位Ricker子波與之褶積得到模擬地震道，然后在模擬地震道中加入了隨機(jī)噪聲作為觀測數(shù)據(jù)(圖1)，子波分別用準(zhǔn)確的35 Hz零相位Ricker子波和37 Hz零相位Ricker子波進(jìn)行反演，反演結(jié)果見圖2～圖5。

圖1 合成道和含有隨機(jī)噪聲的合成道Fig.1 Synthetic record without noise and with random noise(a)合成道;(b)含有隨機(jī)噪聲的合成道

圖2 用不含噪聲的合成道和準(zhǔn)確的子波進(jìn)行反演得到的結(jié)果與真值對比Fig.2 Inversed impedance using synthetic record without noise and accurate wavelet and true impedance

圖3 用含噪的合成道和準(zhǔn)確的子波進(jìn)行反演得到的結(jié)果與真值對比Fig.3 Inversed impedance using synthetic record with random noise and accurate wavelet and true impedance

圖4 用不含噪的合成道和不準(zhǔn)確的子波進(jìn)行反演得到的結(jié)果與真值對比Fig.4 Inversed impedance using synthetic record without noise and inaccurate wavelet and true impedance

圖5 用含噪的合成道和不準(zhǔn)確的子波進(jìn)行反演得到的結(jié)果與真值對比Fig.5 Inversed impedance using synthetic record with noise and inaccurate wavelet and true impedance

圖6 楔形模型真實波阻抗剖面Fig.6 Impedance section of wedge model

圖7 楔形模型合成地震數(shù)據(jù)Fig.7 Synthetic seismic section of wedge model

從圖2可以看出，本文提出的反演方法在一些個別點上還是無法獲得準(zhǔn)確的絕對波阻抗值，但是基本上可以獲得比較可靠的波阻抗變化趨勢以及界面位置。

對比圖2和圖3可以看出，隨機(jī)噪聲對本文的反演方法的影響并不是很大，而且實際資料中的隨機(jī)噪聲一般也比較小，尤其是疊后或者局部疊加數(shù)據(jù)中的隨機(jī)噪聲是比較小的。

對比圖3、圖4、圖5可以看出，子波準(zhǔn)確與否對本文的反演方法的影響還是很大的。雖然反演結(jié)果也能刻畫波阻抗變化趨勢以及界面，但是界面位置以及波阻抗值大小都偏離了真值。

用楔形模型進(jìn)行了試驗。時間域波阻抗楔形模型如圖6所示，時間采樣間隔是1 ms，共1 500道。圖7是用35 Hz零相位Ricker子波與反射系數(shù)進(jìn)行褶積得到的模擬數(shù)據(jù)。我們用此模擬數(shù)據(jù)作為輸入，用本文提出的反演方法進(jìn)行了反演，反演時的子波采用準(zhǔn)確的子波，即35 Hz零相位Ricker子波。圖8是反演得到的波阻抗。對比圖6與圖8，反演能獲得反射界面的真實位置以及真實的波阻抗值。即使兩個反射界面距離非常近，即所謂的薄層情況下，也能準(zhǔn)確刻畫出反射界面的真實位置和波阻抗真實值。

圖8 反演得到的楔形模型波阻抗剖面Fig.8 Inversed impedance section of wedge model

圖9 實際地震剖面Fig.9 A real seismic section

圖10 由圖9的地震數(shù)據(jù)用本文方法反演得到的波阻抗剖面Fig.10 The inversed impedance section from the seismic section of Figure 9 by the proposed method

3 實際地震資料處理

我們對實際資料進(jìn)行了處理。圖9是地震剖面，圖10是相對波阻抗剖面。由于反演時的子波是30 Hz零相位Ricker子波，而不是經(jīng)過標(biāo)定提取出來的，另外由于沒有經(jīng)過標(biāo)定，因此波阻抗的值也只能是相對的。從圖10可見，本反演方法獲得的波阻抗剖面能夠保持反射系數(shù)的稀疏性和橫向連續(xù)性，從縱向上看，波阻抗的變化趨勢以及反射界面還是能夠得到清晰的反映。由于本文的反演僅僅利用了地震數(shù)據(jù)本身的貢獻(xiàn)，而沒有利用包含了測井、地質(zhì)等信息的模型，即沒有利用地震信號帶限外的高低頻信息，尤其是低頻信息，因此本文的反演方法得到的結(jié)果本質(zhì)上是帶限的相對波阻抗,即本文的反演方法本質(zhì)上是無約束稀疏脈沖波阻抗反演的一種實現(xiàn)。為了進(jìn)一步提高反演質(zhì)量，尤其是為了獲得絕對波阻抗值，除了與商業(yè)軟件(圖11)一樣，做好標(biāo)定與子波提取以外，還需要在反演過程中加入低頻模型或者橫向約束。

圖11 某商業(yè)軟件由圖9的地震數(shù)據(jù)反演得到的波阻抗剖面Fig.11 The inversed impedance section from the seismic section of Figure 9 by commercial software

圖12 實際地震剖面Fig.12 Another real seismic section

圖13 由圖12所示地震資料反演后加入低頻趨勢的波阻抗剖面Fig.13 The inversed impedance section from the seismic section of Figure 12

圖14 反演后的井旁道波阻抗與井波阻抗的對比Fig.14 Comparing the inversed impedance with the impedance at well

圖12是另一個實際資料的地震剖面，道長為3 s，采樣間隔1 ms(實際顯示300～1 800)。用本文的反演方法進(jìn)行了反演后加入了低頻趨勢，獲得了波阻抗剖面如圖13所示，反演使用的子波是提取的34 Hz非零相位子波。將井旁道的反演結(jié)果與井上計算的波阻抗進(jìn)行了對比(圖14)。從圖14可見，加入低頻趨勢以后，反演結(jié)果在總體趨勢上與井上的波阻抗吻合得很好，既可以反映反射界面，也可以獲得絕對波阻抗，從而可以用于儲層描述。

4 結(jié)論

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

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