盧明德

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

0 引言

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

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

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

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

1.1 L1TV地震信號重建模型構建

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

(1)

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

f=Kx*+w

(2)

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

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

(3)

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

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

(4)

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

(5)

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

(6)

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

1.2 交替方向乘子法(ADMM)

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

(7)

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

(8)

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

(9)

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

(10)

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

1.3 L1TV地震信號重建及算法

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

(11)

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

(12)

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

(13)

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

(14)

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

(i=1,…,n·m)

(15)

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

(16)

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

(17)

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

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

(18)

2 實驗結果與分析

2.1 合成地震記錄實驗

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

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

(19)

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

(a)原始合成地震記錄

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

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

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

圖3 加噪聲數據、維納濾波方法和本文方法處理結果的第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)增加較強的隨機噪聲，其數據如圖4b所示，隨機噪聲系數為0.15，用于褶積的高斯核函數的標準差σ=1.5。由于維納濾波和拉普拉斯濾波受噪聲影響太大，因此，本文不做討論。圖4c為L1TV-SSD算法的處理結果(μ=4)，從圖中明顯可看出，即使在信號受到較強噪聲污染時，也可以獲得不錯的處理結果，

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

2.2 實際資料處理

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

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

3 結論

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

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

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

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

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