基于L1范數(shù)的全變分地震信號反褶積優(yōu)化算法*

盧明德

(遼河油田公司勘探開發(fā)研究院，遼寧 盤錦 124010)

0 引言

地震波在地層傳播過程中會受到大地低通濾波的影響，造成高頻信息嚴(yán)重衰減，使地震信號的分辨率大大降低。為了使地震數(shù)據(jù)能夠滿足構(gòu)造解釋和油藏勘探等要求，需要對地震資料進行高分辨率處理。反褶積算法是提高地震資料分辨率最常用的方法，在高分辨率地震資料處理中起著至關(guān)重要的作用。到目前為止，多種反褶積方法被提出，如傳統(tǒng)的維納(Wiener)濾波(Treitel，1974)、f-k濾波(Stewart，Schieck，1993)，還有以假設(shè)地層反射信號由稀疏脈沖信號組成的常規(guī)反褶積：L1范數(shù)稀疏脈沖反褶積(Tayloretal，1979)、參數(shù)化稀疏脈沖反褶積(Lei，Terence，2000；Danilo，2006)、L1-L2范數(shù)聯(lián)合反演法(王宇等，2009)、預(yù)條件共軛梯度法(朱振宇，劉洪，2005)等。L1范數(shù)稀疏脈沖反褶積對地震信號中缺失的高頻信息與低頻信息重建較好，但同相軸連續(xù)性較差；預(yù)條件共軛梯度法對單道處理效果不錯，但執(zhí)行多道反褶積時其剖面的連續(xù)性不好?？傊@些方法在提高地震信號分辨率和信噪比上具有相應(yīng)的效果，但均存在一定的局限性(Hennenfentetal，2005；馮志強等，2011)。

近些年，一些新的方法被提出。如Wang等(2016)考慮反射系數(shù)的稀疏性，建立了L1范數(shù)正則化反褶積模型；潘樹林等(2019)提出自適應(yīng)步長的快速迭代閾值收縮法(Beck，Teboulle，2009)，解決稀疏脈沖反褶積問題，具有超線性收斂速率；Pan等(2019)提出了自適應(yīng)的頻率域Bregman稀疏脈沖反褶積算法，在抗噪、自適應(yīng)與計算效率方面體現(xiàn)了較好的優(yōu)勢；馬濤等(2020)實現(xiàn)了利用可控震源力信號反褶積方法提取可控震源單炮記錄；杜鑫等(2021)通過優(yōu)選常規(guī)地震勘探反褶積處理方法，改進常規(guī)地震數(shù)據(jù)處理流程，使地震勘探的分辨率可以滿足城市地下空間調(diào)查的精度要求；張聯(lián)海等(2021)提出由數(shù)據(jù)驅(qū)動的深度卷積神經(jīng)網(wǎng)絡(luò)模型來解決地震反射信號的稀疏反褶積問題。然而，目前大多反褶積方法最主要的限制在于假設(shè)條件的需要，以及提高分辨率的同時也會提高噪聲的能量，降低了地震信號的信噪比。因此，研究出能同時提高分辨率、衰減噪聲，并保持剖面連續(xù)性的反褶積方法仍具挑戰(zhàn)性。

目前，基于全變分(Total Variation，TV)的理論(Rudinetal，1992)在信號去噪、醫(yī)學(xué)信號重建、運動成像等信號處理領(lǐng)域得到了廣泛應(yīng)用(牛和明等，2011；石明珠等，2011；Chenetal，2010)。對于全變分模型范數(shù)的選取，L2范數(shù)一般只針對高斯噪聲有效，然而實際地震信號同時會遭受到異常值等非高斯噪聲的攻擊。針對此種情況，L1范數(shù)能表現(xiàn)出較好的性能，并易于求解(Rudinetal，1992；Yangetal，2009b)。因此，基于數(shù)學(xué)范數(shù)最優(yōu)化理論，本文提出一種L1范數(shù)的全變分地震信號反褶積優(yōu)化算法(L1TV-SSD)，并將該算法對合成信號和野外采集地震數(shù)據(jù)進行實驗。

1 L1范數(shù)全變分地震信號反褶積(L1TV-SSD)

1.1 L1TV地震信號重建模型構(gòu)建

假設(shè)地震記錄的道數(shù)為n，各地震道包含的采樣點數(shù)為m，那么總采樣點數(shù)為n×m個元素。地震數(shù)據(jù)矩陣r可表示為：

(1)

令x*=vec(r)為矩陣r的拉伸向量，即x*∈Rnm；w∈Rnm是地震數(shù)據(jù)中的隨機噪聲，K∈Rnm×nm是褶積算子，f∈Rnm是實際地震信號，滿足如下關(guān)系：

f=Kx*+w

(2)

在K給定的前提下，從含噪聲的實際地震信號f中重建初始地震信號x*，即重建地震信號可視為反褶積和去噪同時處理。若K為單位算子，則此操作僅為去噪處理。依據(jù)Rudin等(1992)，Vogel和Oman(1998)研究，=Kx-f=1≈=w=1，即使用重建初始地震信號x*可對公式(3)進行最小化處理：

Φfid(x，f)==kx-f=1

(3)

式中：=·=1為L1范數(shù)；Φfid(x,f)是噪聲信號的保真項。

本文目標(biāo)是從f中重建x*，但此過程對噪聲w過于敏感，重建x*存在不穩(wěn)定性。為解決此問題，需融入TV正則項(Rudinetal，1992；Wangetal，2008)，即TV的離散形式：

(4)

式中：=·=2為L2范數(shù)；Di為局部有限差分算子；Dix∈R2為x在采樣點i處垂直方向與水平方向的一階有限差分。地震信號的重建過程為求解式(5)的最小化(Yangetal，2009a)：

(5)

即地震信號的重建模型為：

(6)

式中：μ>0為平衡二者的參數(shù)。注意：x*是理論上的真實地震信號數(shù)據(jù)；x為被恢復(fù)后的地震數(shù)據(jù)。本文通過交替方向乘子法(ADMM)求解式(6)。

1.2 交替方向乘子法(ADMM)

ADMM的基本思想需追溯到Gabay和Mercier(1976)的研究，即令θ1(·)和θ2(·)為凸函數(shù)，A為連續(xù)的線性算子，則最小化能量函數(shù)為：

(7)

通過引入輔助變量v，式(7)可等價轉(zhuǎn)換為：

(8)

式(8)的增廣Lagrangian函數(shù)為(Gabay，Mercier，1976)：

(9)

式中：λ和β為增廣Lagrangian求解最優(yōu)化的引入?yún)?shù)。因此，交替方向乘子法的迭代解如下：

(10)

式(10)的交替最小化迭代思想為：首先由初始化(vk，λk)獲取uk+1，再使用uk+1和λk獲取vk+1，從而可得到λk+1。最后，利用交替方向最小化的多次優(yōu)化來獲得恢復(fù)的地震信號。

1.3 L1TV地震信號重建及算法

首先，將式(6)轉(zhuǎn)化成式(8)形式，可表示為：

(11)

式中：y=(y1；y2)∈R2nm，y1和y2是長度為n×m維向量， 滿足((y1)i；(y2)i)=yi∈R2。利用二次罰函數(shù)法把式(11)轉(zhuǎn)成無約束優(yōu)化函數(shù)形式：

(12)

式中：β?0為罰參數(shù)。通過ADMM對式(12)進行求解(初始x=xk，λ=λk)，可得到：

(13)

式中：LA(x，y，λ)為式(12)的增廣Lagrangian函數(shù)：

(14)

式中：λi為R2空間中的向量。根據(jù)Wang等(2008)和Yang等(2009a)研究，式(13)中argminyLA(xk,y,λk)等價于如下形式：

(i=1,…,n·m)

(15)

式(15)通過二維收縮來求解(Wangetal，2008；Yangetal，2009a)：

(16)

在固定yk+1和λk的基礎(chǔ)上，通過最小二乘法求解式(13)中argminxLA(x,yk+1,λk)：

(17)

最后，由式(18)更新λ：

λk+1=λk-γβ(yk+1-Dxk+1)

(18)

2 實驗結(jié)果與分析

2.1 合成地震記錄實驗

為了對L1TV-SSD算法進行有效性評價，將其應(yīng)用于合成地震記錄上進行實驗，并與維納濾波和拉普拉斯濾波結(jié)果進行比較。實驗中選取參數(shù)為：γ=1.618，ε=10-3，罰參數(shù)β=10；為了使地震信號重建獲得最好的結(jié)果，不同的地震信號數(shù)據(jù)，μ值取值不同，范圍一般在(0，60]。

采用信噪比(Signal-Noise Ratio，簡稱SNR)來評價地震信號重建的質(zhì)量，SNR定義為：

(19)

式中：x*為原始地震記錄；E(x*)為x*的均值強度；x是重建后的地震信號數(shù)據(jù)；SNR單位為dB。

(a)原始合成地震記錄

原始合成地震記錄主頻為35 Hz、1 ms采樣、共150道，每道采樣點數(shù)為600，地震波為雷克子波(圖1a)。使用高斯核函數(shù)(標(biāo)準(zhǔn)差σ=2)與原始地震記錄褶積得到數(shù)據(jù)的SNR為7.69 dB。從圖1b中可以看出，箭頭指向處主頻降低。使用L1TV-SSD算法對圖1b數(shù)據(jù)進行處理(μ=5)，箭頭指向處的分辨率明顯提高，SNR值提高至13.88 dB(圖1c)。圖1d為SNR、優(yōu)化目標(biāo)函數(shù)值與L1TV-SSD算法迭代數(shù)的關(guān)系曲線，在迭代48次后SNR持續(xù)增加，優(yōu)化目標(biāo)函數(shù)值隨迭代次數(shù)的增加一直逐漸降低，在近140次時降到最低。

圖1e與f分別是使用維納濾波與拉普拉斯濾波的處理結(jié)果，SNR分別為10.42 dB與9.13 dB。維納濾波的處理效果提升不大，拉普拉斯濾波處理結(jié)果其主頻有所提高，但多了噪聲。綜上，L1TV-SSD算法獲得了較好的處理效果。

對原始地震記錄(圖2a)進行加噪處理后進行實驗。使用高斯核函數(shù)和原始地震記錄進行褶積，再通過Matlab中“randn”函數(shù)增加隨機噪聲的數(shù)據(jù)如圖2b所示。高斯核函數(shù)的標(biāo)準(zhǔn)差σ=1.5，隨機噪聲系數(shù)設(shè)置為0.02。使用L1TV-SSD算法對圖2b數(shù)據(jù)進行處理(μ=4)，從圖2c可明顯看出，L1TV-SSD算法有效地衰減了隨機噪聲，提高了分辨率，SNR由5.79 dB增至11.98 dB。從圖2d可以看出，隨著迭代次數(shù)的增加，SNR不斷提升，而目標(biāo)函數(shù)值不斷降低。對于維納濾波結(jié)果(圖2e)和拉普拉斯濾波結(jié)果(圖2f)，盡管二者的同相軸較圖2b清晰，但增加了較多的噪聲，使得信噪比下降，這表明這兩種結(jié)果極易受噪聲的影響。圖3為加噪聲數(shù)據(jù)、維納濾波方法和本文方法處理結(jié)果的第125道頻譜對比圖。從圖中可以看到，L1TV-SSD算法提高了子波的主頻，拓寬了有效頻帶，起到了反褶積的作用，高頻隨機噪聲同時也得到了較好的壓制；而維納濾波方法引入了較多的隨機噪聲，降低了子波的主頻。

圖3 加噪聲數(shù)據(jù)、維納濾波方法和本文方法處理結(jié)果的第125道頻譜對比Fig.3 Comparison between the spectrums of the 125thtrace respectively obtained by the noisy data method，the Wiener filtering method and the proposedmethod in this paper

對原始地震記錄(圖4a)增加較強的隨機噪聲，其數(shù)據(jù)如圖4b所示，隨機噪聲系數(shù)為0.15，用于褶積的高斯核函數(shù)的標(biāo)準(zhǔn)差σ=1.5。由于維納濾波和拉普拉斯濾波受噪聲影響太大，因此，本文不做討論。圖4c為L1TV-SSD算法的處理結(jié)果(μ=4)，從圖中明顯可看出，即使在信號受到較強噪聲污染時，也可以獲得不錯的處理結(jié)果，

SNR由-8.48 dB增至7.60 dB。同樣，從圖4d可以看到，隨著迭代次數(shù)的增加，SNR不斷提升，而優(yōu)化目標(biāo)函數(shù)值不斷降低。

2.2 實際資料處理

為了更好地說明L1TV-SSD算法的適用性，將其應(yīng)用于遼河某地區(qū)三維地震資料處理中。圖5a為一初始疊加剖面，信噪比較低；圖5b為L1TV-SSD算法處理結(jié)果；圖5c和d分別為應(yīng)用GeoEast軟件的徑向預(yù)測濾波(RadiPredicFilt)和隨機噪聲衰減(RNA)兩個模塊的處理結(jié)果。3種方法的對比剖面分析表明：L1TV-SSD算法對原始疊加剖面的背景噪聲具有更好的壓制效果，并且對模糊的同相軸有更好的處理效果，使其連續(xù)性增強，同相軸更加清晰，層間信息更加豐富，整體的分辨率明顯提高。從RadiPredicFilt和RNA兩個模塊來看，RNA在噪聲衰減以及細節(jié)的清晰度方面要優(yōu)于RadiPredicFilt。

為了進一步驗證L1TV-SSD算法的有效性，選取實際的CDP道集(圖6a)進行處理，從圖6a可以看出，實際資料存在很強的隨機噪聲，對有效信號造成了強烈的干擾。圖6b～d分別為使用L1TV-SSD算法、徑向預(yù)測濾波(RadiPredicFilt)和隨機噪聲衰減(RNA)模塊對實際的CDP道集處理的結(jié)果。從處理效果來看，3種方法均有效地壓制了隨機噪聲，L1TV-SSD算法有其更好的壓制效果，突顯出來的細節(jié)信息要優(yōu)于其它兩種方法。

3 結(jié)論

地震信號質(zhì)量的提高可為成像的精確性奠定基礎(chǔ)。本文提出一種提升地震信號分辨率和信噪比的新算法——L1TV-SSD算法，為了對該算法進行有效性評價，將其應(yīng)用于合成地震記錄和遼河凹陷野外采集數(shù)據(jù)進行實驗，主要結(jié)論如下：

(1)基于L1范數(shù)全變分理論構(gòu)建地震信號重建模型，通過交替方向乘子法解決轉(zhuǎn)化后最小化問題，設(shè)計地震信號優(yōu)化處理算法，同時實現(xiàn)反褶積與去噪處理。

(2)使用L1TV-SSD算法、維納濾波和拉普拉斯濾波對合成數(shù)據(jù)處理后發(fā)現(xiàn)：維納濾波的處理效果提升不大，拉普拉斯濾波處理結(jié)果其主頻有所提高，但多了噪聲，而L1TV-SSD算法處理效果分辨率明顯提高。

(3)從對含噪聲的合成地震處理結(jié)果發(fā)現(xiàn)：對于遼河凹陷野外采集的數(shù)據(jù)，L1TV-SSD算法提高了子波的主頻，有效提高了地震資料的分辨率和信噪比。

(4)本文提出的方法避免了常規(guī)地震信號反褶積和去噪相互影響的不足，提供了新的處理思路。