盧明德

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

0 引言

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

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

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

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

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

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

(1)

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

f=Kx*+w

(2)

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

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

(3)

式中：=·=1為L(zhǎng)1范數(shù)；Φfid(x,f)是噪聲信號(hào)的保真項(xiàng)。

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

(4)

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

(5)

即地震信號(hào)的重建模型為：

(6)

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

1.2 交替方向乘子法(ADMM)

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

(7)

通過(guò)引入輔助變量v，式(7)可等價(jià)轉(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)化來(lái)獲得恢復(fù)的地震信號(hào)。

1.3 L1TV地震信號(hào)重建及算法

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

(11)

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

(12)

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

(13)

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

(14)

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

(i=1,…,n·m)

(15)

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

(16)

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

(17)

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

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

(18)

2 實(shí)驗(yàn)結(jié)果與分析

2.1 合成地震記錄實(shí)驗(yàn)

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

采用信噪比(Signal-Noise Ratio，簡(jiǎn)稱SNR)來(lái)評(píng)價(jià)地震信號(hào)重建的質(zhì)量，SNR定義為：

(19)

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

(a)原始合成地震記錄

原始合成地震記錄主頻為35 Hz、1 ms采樣、共150道，每道采樣點(diǎn)數(shù)為600，地震波為雷克子波(圖1a)。使用高斯核函數(shù)(標(biāo)準(zhǔn)差σ=2)與原始地震記錄褶積得到數(shù)據(jù)的SNR為7.69 dB。從圖1b中可以看出，箭頭指向處主頻降低。使用L1TV-SSD算法對(duì)圖1b數(shù)據(jù)進(jìn)行處理(μ=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次時(shí)降到最低。

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

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

圖3 加噪聲數(shù)據(jù)、維納濾波方法和本文方法處理結(jié)果的第125道頻譜對(duì)比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

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

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

2.2 實(shí)際資料處理

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

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

3 結(jié)論

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

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

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

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

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