應用近似L0范數(shù)的稀疏脈沖反演

劉百紅 李建華 鄭四連

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

1 問題的提出

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

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

(1)

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

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

(2)

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

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

(3)

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

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

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

(4)

式中：k表示稀疏度；p可等于1，也可等于2。當p=1時，式(4)轉(zhuǎn)化為線性規(guī)劃問題，即可用線性規(guī)劃算法求解； 當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)

該式可用許多通用的基于導數(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，重復第二步；

四、最終獲得最優(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ù)大小。因此，在實際應用時，將地震資料歸一化以后，在第一層循環(huán)中，從大到小、依次按照一定間隔取一系列σ。只要設定最大值、最小值及間隔步長即可。

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

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

3 模型試驗

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

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

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

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

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

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

4 實際地震資料處理

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

利用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ì)量，在做好標定與子波提取的基礎上，還需在反演過程中加入低頻模型或者橫向約束。

圖5 實際過井地震剖面L1

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

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

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

圖9 實際地震剖面L2

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

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

5 結(jié)論

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

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