基于L1范數的波阻抗反演

劉百紅 ， 李建華 ， 鄭四連

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

0 引言

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

1 基本原理

在反射地震資勘探中，疊后地震道的形成可以用下面的褶積模型來描述：

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

(1)

其中:S(t)表示地震道；r(t)表示反射系數序列；W(t)表示子波；n(t)表示噪音；t為時間；*表示褶積。當已獲得疊后數據體時，希望由式(1)得到與數據相應的地下介質的參數，即反射系數序列，進而得到波阻抗，這就是疊后波阻抗反演。在這個過程中通常都是假設反射系數序列是稀疏的，并通過井震標定確定出子波，然后利用某種最優化算法求解設定好的目標函數，來獲得反射系數序列。最基本的目標函數是下面的誤差函數：

E(r)=‖S-d0‖2

(2)

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

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

(3)

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

實際上，當p=2、q=1時，式(3)僅僅是無約束優化的一般形式?？紤]到反射系數的定義：

(4)

應該有約束：|r|≤1。

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

(5)

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

(6)

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

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

(7)

(8)

其中，ε是一個很小的數，例如可以將其設置為：ε=10-10。則目標函數變為：

(9)

這樣，式(7)中的目標函數就變為連續可微，從而可以用基于導數的最優化方法求解式(9)。

2 模型試驗

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 實際地震資料處理

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

圖11 某商業軟件由圖9的地震數據反演得到的波阻抗剖面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)。用本文的反演方法進行了反演后加入了低頻趨勢，獲得了波阻抗剖面如圖13所示，反演使用的子波是提取的34 Hz非零相位子波。將井旁道的反演結果與井上計算的波阻抗進行了對比(圖14)。從圖14可見，加入低頻趨勢以后，反演結果在總體趨勢上與井上的波阻抗吻合得很好，既可以反映反射界面，也可以獲得絕對波阻抗，從而可以用于儲層描述。

4 結論

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

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