基于L1/2正則化理論的地震稀疏反褶積

康治梁,張雪冰

(1.成都理工大學(xué)地球物理學(xué)院,四川成都610059;2.加拿大多倫多大學(xué)物理學(xué)院,安大略多倫多,ON M5S 1A7;3.吉林建筑大學(xué)測繪與勘查工程學(xué)院,吉林長春130118)

隨著油氣勘探開發(fā)的不斷深入,小尺度(薄互層)巖性油氣藏的識別已成為愈發(fā)重要的問題[1],在褶積模型的假設(shè)下,要求能夠從地震信號中反褶積出準(zhǔn)確的反射系數(shù)信息,包括它的位置和振幅[2]。有限頻地震信號的反褶積常常是一個多解問題,而在實(shí)際應(yīng)用中,地質(zhì)模型的成層性決定了我們需要的是具有稀疏性質(zhì)的解[2]。壓縮感知理論的信號恢復(fù)問題就是求解信號在某個基底下的稀疏表示系數(shù),這與地震反褶積有類似之處。VELIS[2]將模擬退火算法和線性最小二乘法相結(jié)合,求取反射系數(shù)序列中脈沖的時間位置和振幅,然而模擬退火算法非常耗時,應(yīng)用于大規(guī)模數(shù)據(jù)的計(jì)算成本較高。匹配追蹤(MP)算法可用于稀疏問題的求解,YANG等[3]利用MP算法實(shí)現(xiàn)了壓縮感知框架下的重力信號恢復(fù),但MP算法通常在信號超稀疏且噪聲不顯著時才有優(yōu)勢[4]。目前廣泛使用的稀疏問題求解辦法是L1正則化約束。ZHANG等[5]認(rèn)為地震反射系數(shù)可以用一系列代表楔狀模型的奇偶字典線性表示,求取表示系數(shù)時使用L1范數(shù)進(jìn)行稀疏性約束,實(shí)現(xiàn)了模型約束意義下的反演。OLIVEIRA等[6]使用L1正則化約束實(shí)現(xiàn)了衰減介質(zhì)中的稀疏反褶積。余慧敏等[7]利用L1正則化約束實(shí)現(xiàn)了壓縮感知框架下的探地雷達(dá)成像。蔣川東等[8]實(shí)現(xiàn)了基于L1范數(shù)的低場核磁共振T2譜稀疏反演。然而,壓縮感知理論的一些實(shí)踐表明,在稀疏恢復(fù)能力和適應(yīng)噪聲方面,L1正則化不是最佳的稀疏正則化[9],這也啟示我們必須在地震稀疏反褶積中尋找更好的稀疏約束方法。

XU等[10-12]對非凸的Lq(0

1 方法原理

1.1 稀疏反褶積

在自激自收條件下,ROBINSON[17]提出利用褶積模型描述地震響應(yīng),地震記錄可近似看作是地震子波與地層反射系數(shù)的褶積結(jié)果,即:

s(t)=w(t)*r(t)+n(t)

(1)

式中:s(t)表示合成地震記錄;w(t)表示地震子波;r(t)表示反射系數(shù)序列;n(t)表示隨機(jī)噪聲。在矩陣格式下(1)式可以表達(dá)為:

s=Wr+n

(2)

式中:W為卷積核矩陣。求取反射系數(shù)問題通常存在多解性,然而層狀地層的假設(shè)決定了反射系數(shù)r是一個稀疏脈沖序列,于是反褶積問題可以描述為在約束條件下的稀疏求解問題。

min‖r‖0s.t.s=Wr+n

(3)

式中:‖r‖0是L0范數(shù),代表向量的非零項(xiàng)系數(shù)的個數(shù)。類似(3)式的稀疏問題也可以從壓縮感知理論導(dǎo)出。由壓縮感知理論可知,當(dāng)信號在某個基底下能稀疏表示時(如余弦變換基、小波變換基),可以用低于奈奎斯特采樣率的方式感知信號,且仍能高概率地恢復(fù)出信號[18]。壓縮感知的信號恢復(fù)問題可以表示為:

min‖x‖0s.t.y=φψx+ε

(4)

式中:ψ表示基底;φ表示測量矩陣;y表示測量值(維度小于原始信號維度);x表示原始信號的稀疏表示系數(shù);ε表示噪聲。可令A(yù)=φψ。(3)式和(4)式都可以直接建模為L0正則化問題[12],以反褶積為例:

(5)

式中:λ是正則化參數(shù)。我們面臨的問題是,L0正則化問題是一個NP-hard問題,無法直接求解。已有的許多關(guān)于壓縮感知和稀疏反演的文獻(xiàn)中,問題((5)式)通常凸松弛為L1范數(shù)問題[12]:

(6)

L1范數(shù)問題是凸優(yōu)化問題,常見算法包括梯度投影法、同倫法、迭代收縮閾值法、近似梯度法和增廣拉格朗日乘子法[19]。本文使用迭代收縮閾值算法(iterative shrinkage-thresholding algorithm,ISTA),BECK等[20]給出了該算法的詳細(xì)推證。ISTA算法源自一般化的最小化問題:

(7)

(8)

式中:tk為步長。對于L1范數(shù)問題((6)式),f(x)=‖s-Wr‖2,g(x)=‖r‖1,(8)式進(jìn)一步轉(zhuǎn)化為:

(9)

問題(9)由于變量可分,rk的每一個分量可以通過求解一個簡單的軟閾值函數(shù)Tsoft得到:

rk=Tsoft[(rk-1-2tkWT(Wrk-1-s)]

(10)

同時算法收斂要求步長tk滿足tk∈(0,1/‖WT·W‖)。

盡管L1范數(shù)約束已被廣泛使用,但壓縮感知領(lǐng)域的研究結(jié)果表明,使用L1范數(shù)約束不能保證有最稀疏的結(jié)果,也不能確保在較小采樣率的情況下精確恢復(fù)信號[12],此外,從地震稀疏反褶積得到的結(jié)果也能發(fā)現(xiàn),L1范數(shù)約束得到的結(jié)果在幅值上與真實(shí)反射系數(shù)的擬合度也有不足[21]。

1.2 L1/2正則化理論

針對L1范數(shù)存在的問題,XU等[10-12]系統(tǒng)地研究了形如(11)式的Lq正則化問題,提出了“L1/2正則化理論”。

(11)

XU等[10-12]通過壓縮感知研究方面的大量計(jì)算證實(shí),q∈(0,1/2]時,Lq正則化性能沒有顯著差異,q∈(1/2,1]時,Lq正則化性能遞減,他們的工作揭示了L1/2正則化在各種稀疏約束中的突出的代表性現(xiàn)象。L1/2正則化是非凸的,常規(guī)的凸優(yōu)化工具難以求解,XU等[12]基于嚴(yán)格的數(shù)學(xué)分析給出了求解L1/2正則化問題的特定閾值迭代算法,使得該問題能快速高效地求解。該算法的迭代格式為:

xn+1=Hλ(B(xn))

(12)

Hλ(x)=(hλ(x1),hλ(x2),…,hλ(xn))T

(13)

其中,hλ(·)定義為:

(14)

XU等[12]有針對性地對比了L1/2正則化和L1正則化在壓縮感知信號恢復(fù)問題研究方面的性能差異,在有噪聲和無噪聲條件下,L1/2正則化的恢復(fù)精度都更高,且對欠采樣的穩(wěn)健性更強(qiáng)。這些工作為我們提供了另一種地震稀疏反褶積的思路,即將L1/2范數(shù)作為反射系數(shù)的稀疏約束項(xiàng),并用該特定閾值迭代算法進(jìn)行求解。

2 數(shù)值實(shí)驗(yàn)

本文假設(shè)子波已知,實(shí)驗(yàn)中選取主頻為30Hz的雷克子波(圖1)。在實(shí)際處理中,子波可以從地震記錄中通過盲提取獲得,如利用高階統(tǒng)計(jì)量的方法[22]。

2.1 簡單模型

在各類正則化問題中,正則化參數(shù)λ對反演效果有顯著的影響。通常依賴經(jīng)驗(yàn)來設(shè)定λ,為了對比在常規(guī)取值范圍內(nèi)、不同λ條件下L1/2范數(shù)約束和L1范數(shù)約束的反演效果,構(gòu)建了一個簡單的反射系數(shù)序列(圖2a),時間采樣間隔為2ms,合成對應(yīng)的地震數(shù)據(jù),并加入15%的噪聲,結(jié)果見圖2b所示。本文使用SRE(signal-to-reconstruction error)來評價反演效果,SRE定義為:

(15)

圖1 實(shí)驗(yàn)中使用的子波

圖5a是構(gòu)建的另一個反射系數(shù)模型,圖5b是對應(yīng)的加入了15%噪聲的合成地震記錄。分別用L1和L1/2范數(shù)約束進(jìn)行反褶積,結(jié)果如圖6a和圖6b所示,圖7a和圖7b是對應(yīng)的反演結(jié)果與真實(shí)反射系數(shù)的殘差。綜合圖6和圖7可以看出,這兩種稀疏約束均能得到較準(zhǔn)確的反射系數(shù)模型形態(tài)特征,但L1/2范數(shù)約束的殘差更小,能更好地保護(hù)反射系數(shù)幅值。

圖2 簡單的稀疏反射系數(shù)模型(a)及其含噪的合成地震數(shù)據(jù)(b)

圖3 SRE值與λ的關(guān)系(橫坐標(biāo)為對數(shù)坐標(biāo))

圖4 SRE值與噪聲水平的關(guān)系

圖5 簡單反射系數(shù)模型(a)及其含噪的合成地震記錄(b)

圖6 采用L1(a)和L1/2(b)正則化約束得到的反褶積結(jié)果

圖7 采用L1(a)和L1/2(b)正則化約束的殘差

2.2 復(fù)雜模型

利用Marmousi2 P波速度模型的局部(圖8紅色方框中)構(gòu)建了一個反射系數(shù)模型,如圖9a所示。圖9b和圖9c分別為對應(yīng)的不含噪聲的和含20%噪聲的地震記錄,共150道,時間采樣間隔為2ms。

圖10是采用不含噪聲的數(shù)據(jù)反演得到的結(jié)果,其中圖10a是采用L1范數(shù)約束得到的結(jié)果,圖10b是采用L1/2范數(shù)約束得到的結(jié)果。大體來看,采用這兩種方法都能得到反射系數(shù)模型主要格架,但從標(biāo)注紅色方框的部分可以看出,采用L1范數(shù)約束得到的結(jié)果中存在反射系數(shù)缺失的情況,不利于橫向連續(xù)性的保持和弱反射系數(shù)的保護(hù),而采用L1/2范數(shù)約束能較好地保存反射系數(shù)信息,橫向連續(xù)性好,刻畫地層尖滅的能力更強(qiáng)。

圖11是采用含噪聲的數(shù)據(jù)反演得到的結(jié)果,其中圖11a是采用L1范數(shù)約束得到的結(jié)果,圖11b是采用L1/2范數(shù)約束得到的結(jié)果。由于噪聲的存在,此時采用這兩種方法反演的質(zhì)量相比無噪聲情形下的結(jié)果都有所下降。從紅色方框指示的位置可以看到,采用L1范數(shù)約束得到的反射系數(shù)模型存在明顯的層位缺失的現(xiàn)象,有的弱反射層被整體破壞;采用L1/2范數(shù)約束得到的反射系數(shù)模型雖然也存在缺失或振幅削弱的地方,但尚能保存大致輪廓,對地層格架特征的描述依然優(yōu)于L1范數(shù)約束。

圖8 Marmousi2 P波速度模型(空間采樣間隔為1.25m)

圖9 反射系數(shù)模型(a)和不含噪聲的合成記錄(b)及含噪聲的合成記錄(c)

圖10 采用不含噪聲數(shù)據(jù)的反演結(jié)果a L1正則化約束的結(jié)果; b L1/2正則化約束的結(jié)果

圖11 采用含噪聲數(shù)據(jù)的反演結(jié)果a L1正則化約束的結(jié)果; b L1/2正則化約束的結(jié)果

3 實(shí)際數(shù)據(jù)測試

為了驗(yàn)證上述數(shù)值實(shí)驗(yàn)結(jié)果,對某實(shí)際二維地震數(shù)據(jù)進(jìn)行了應(yīng)用測試,其疊加剖面如圖12所示。圖13 是用統(tǒng)計(jì)性方法從實(shí)際地震記錄中提取的子波。該區(qū)域薄層發(fā)育,縱向分辨率的提升十分必要。

圖14是L1/2正則化約束反演結(jié)果,從標(biāo)注紅色方框的部分和箭頭指示位置來看,L1/2正則化約束能夠較好地揭示某些薄層結(jié)構(gòu)的稀疏反射系數(shù),有較強(qiáng)的刻畫薄層和微小透鏡體(箭頭指示的部分)的能力,也能較好地保持原始地層模型的縱向展布特征。圖15 展示了反褶積前、后地震數(shù)據(jù)頻譜的變化,可以看出,采用本文方法反褶積之后,子波影響得到了有效消除,地震資料有效頻帶寬度顯著提高,剖面的地層結(jié)構(gòu)分辨率得到了明顯提升。

圖12 實(shí)際地震數(shù)據(jù)疊加剖面

圖13 從實(shí)際地震記錄中提取的子波

圖14 采用L1/2正則化約束反演結(jié)果

圖15 反褶積前、后的地震數(shù)據(jù)頻譜

4 結(jié)束語

本文將用于壓縮感知領(lǐng)域的L1/2正則化理論引入到地震稀疏反褶積中,采用了它的特定閾值迭代算法進(jìn)行求解,形成了一種新的非凸正則化約束的稀疏反褶積方法。理論合成數(shù)據(jù)測試結(jié)果表明,L1/2正則化的表現(xiàn)優(yōu)于L1正則化,適應(yīng)噪聲的能力更強(qiáng),能更好地保護(hù)弱反射系數(shù)振幅。實(shí)際資料應(yīng)用結(jié)果表明,L1/2正則化突出的稀疏表達(dá)能力,使它能夠很好地刻畫薄層結(jié)構(gòu)和透鏡體,為地震資料高分辨率處理提供了有力的工具。當(dāng)然,作為一種數(shù)據(jù)處理手段,無論采用何種算法,反褶積的結(jié)果都需要其它手段驗(yàn)證,如測井信息的標(biāo)定,這是實(shí)際工作中不能忽視的。

稀疏問題在地震信號處理中普遍存在,如基于原子庫的信號稀疏分解,將信號用一系列具有特定時頻特征的原子稀疏表示,可以達(dá)到去噪的目的,用L1/2正則化理論求取稀疏表示系數(shù),有助于分離有效信息和噪聲,這也是下一步研究方向之一。