焦敘明, 杜啟振, 趙 強(qiáng)

(1.中海油服物探事業(yè)部物探研究院,天津 300000;2.中國(guó)石油大學(xué)(華東)深層油氣重點(diǎn)實(shí)驗(yàn)室,山東青島 266580;3.青島海洋科學(xué)與技術(shù)國(guó)家實(shí)驗(yàn)室海洋礦產(chǎn)資源評(píng)價(jià)與探測(cè)技術(shù)功能實(shí)驗(yàn)室,山東青島 266237;4.中國(guó)石油大學(xué)(華東)CNPC物探重點(diǎn)實(shí)驗(yàn)室,山東青島 266580)

在實(shí)際地震資料采集過(guò)程中,受野外采集環(huán)境以及檢波器本身等因素的影響,接收到的地震數(shù)據(jù)往往存在缺道、壞道等問(wèn)題。地震道的缺失導(dǎo)致地震信號(hào)空間采樣的不完備,但常規(guī)地震資料處理往往基于信號(hào)完備的基本假設(shè)。缺失的信息作為干擾源,影響到多次波壓制、偏移成像以及道集分析等后續(xù)地震資料處理[1-7]。重建缺失地震數(shù)據(jù)的方法可大致劃分為3類(lèi):波動(dòng)方程法、預(yù)測(cè)濾波法以及信號(hào)變換法。波動(dòng)方程法通過(guò)先驗(yàn)的速度場(chǎng)等信息,對(duì)地震道缺失位置進(jìn)行補(bǔ)全。但波動(dòng)方程法往往會(huì)引入巨大的計(jì)算量,同時(shí)先驗(yàn)速度場(chǎng)信息的準(zhǔn)確性也會(huì)極大地影響重建地震道的精度。預(yù)測(cè)濾波法利用相鄰地震道的相似性以及可預(yù)測(cè)性,通過(guò)構(gòu)建合適的預(yù)測(cè)因子,以現(xiàn)有的地震道為輸入,對(duì)相鄰的缺失地震道進(jìn)行預(yù)測(cè)重建,但其預(yù)測(cè)精度往往取決于預(yù)測(cè)因子的準(zhǔn)確性以及相鄰地震道的相似性。信號(hào)變換法以數(shù)學(xué)變換作為基本手段,以變換域中信號(hào)的稀疏性以及缺道干擾的非稀疏性作為前提條件,對(duì)缺失信息進(jìn)行重建。針對(duì)如何利用信號(hào)稀疏性進(jìn)行地震缺失信息重建的問(wèn)題,國(guó)內(nèi)外學(xué)者進(jìn)行了大量的研究[8-12]。常規(guī)的信號(hào)變換法通常利用L2范數(shù)作為殘差度量準(zhǔn)則,可較穩(wěn)定地處理有隨機(jī)噪聲干擾的重建問(wèn)題[13]。實(shí)際資料中廣泛分布的異常振幅噪聲往往不滿足L2范數(shù)所需要的正態(tài)分布,這會(huì)嚴(yán)重影響該類(lèi)算法的性能,進(jìn)而影響地震道重建效果。一般來(lái)說(shuō),提高傳統(tǒng)信號(hào)變換法的手段主要有兩類(lèi):合適的異常振幅噪聲檢測(cè)方法[14-16]以及魯棒的度量標(biāo)準(zhǔn)[17-18]?，F(xiàn)階段基于魯棒度量標(biāo)準(zhǔn)的算法普遍需要求取噪聲的方差作為輸入,利用中位絕對(duì)離差通?？稍诋惓Ｕ穹肼暣嬖诘那闆r下對(duì)噪聲方差進(jìn)行大致的估計(jì)。然而在實(shí)際資料處理過(guò)程中,噪聲情況往往是未知的,這對(duì)中位絕對(duì)離差的估計(jì)造成了一定的不確定性。Liu等[19]通過(guò)引入高斯核函數(shù)對(duì)殘差進(jìn)行度量,提出了最大相關(guān)熵準(zhǔn)則(maximum correntropy criterion, MCC)。最大相關(guān)熵準(zhǔn)則以誤差的局部相關(guān)熵的最大化作為目標(biāo)函數(shù),利用高斯核函數(shù)對(duì)預(yù)測(cè)值與觀測(cè)值的誤差進(jìn)行正態(tài)化約束,使殘差結(jié)果趨向于零均值,有效地限制了異常振幅噪聲對(duì)于有效信號(hào)估計(jì)的影響?；谶@種相關(guān)熵準(zhǔn)則的相關(guān)算法,在異常振幅噪聲影響下仍可表現(xiàn)出較好的性能[18-20]。筆者從含異常振幅噪聲的缺失地震道重建問(wèn)題出發(fā),引入相關(guān)熵,構(gòu)建以最大相關(guān)熵為目標(biāo)函數(shù)的魯棒地震道重建模型。

1 基于壓縮感知理論的含噪缺失地震道重建

不同于香農(nóng)采樣定理,壓縮感知重建所需的樣本數(shù)不由信號(hào)本身的帶寬決定,而是由信號(hào)的稀疏結(jié)構(gòu)及其采樣方式的有限等距性質(zhì)決定[22]。在地震勘探過(guò)程中,受地理環(huán)境等因素的影響,容易出現(xiàn)不規(guī)則的地震道缺失、壞道的問(wèn)題。不規(guī)則的地震道缺失保證了壓縮感知理論應(yīng)用于地震道恢復(fù)的可能性。除此之外,野外采集信號(hào)往往受各種噪聲的干擾。干擾噪聲一般可分為兩種:相干噪聲以及非相干噪聲。相干噪聲指代面波、多次波等波形噪聲,非相干噪聲指代隨機(jī)噪聲以及奇異振幅噪聲等非波形噪聲。考慮到地震道缺失以及非相干噪聲的影響,地震資料采集可表示為

y=RΩ(x+n+e).

(1)

基于壓縮感知理論的信號(hào)x的恢復(fù)過(guò)程可表述為求解其最稀疏解α,即

(2)

(3)

式(3)所表示的最優(yōu)化問(wèn)題可利用迭代閾值算法[22]進(jìn)行求解

ak+1=Tλ(ak+ΦT(y-RΩΦak)),

(4)

其中

Tλ(a)=sign(a)max(0,|a|-λ).

式中,ΦT為逆稀疏變換基;ak和ak+1分別為第k、k+1次迭代的稀疏表示結(jié)果;λ為迭代過(guò)程的閾值;sign(a)表示取α的正負(fù)號(hào)。

理想情況下,

φ(ri)=ri.

(5)

式中,ri為恢復(fù)結(jié)果相對(duì)于觀測(cè)結(jié)果的殘差。式(5)中的敏感度指示了對(duì)應(yīng)的殘差在L2范數(shù)度量中的影響程度。圖1(a)為不同殘差ri所對(duì)應(yīng)的L2范數(shù)度量值曲線,圖1(b)為不同ri所對(duì)應(yīng)的導(dǎo)數(shù)值。可以看出,L2范數(shù)的導(dǎo)數(shù)隨差值的增加而遞增,說(shuō)明L2范數(shù)的敏感度也隨著差值的增大而更加明顯。利用隨機(jī)生成1 000個(gè)分布于[-10,10]之間的零均值正態(tài)分布數(shù)據(jù),并額外添加10個(gè)相同的異常值,得到的不同異常值對(duì)原始零均值正態(tài)分布數(shù)據(jù)的影響如圖1(c)所示。由于L2范數(shù)不是魯棒型的估計(jì),它的代價(jià)函數(shù)是基于最小均方誤差準(zhǔn)則的,其統(tǒng)計(jì)量可無(wú)限增加,觀測(cè)結(jié)果中較大的異常值的影響在L2范數(shù)估計(jì)中是不受限的,所造成的影響將會(huì)被放大。另一方面,如圖1(c)所示,均值的估計(jì)誤差受異常值的影響明顯,呈線性增加的趨勢(shì),即異常值數(shù)值越大,均值的估計(jì)偏差越大。對(duì)于式(3)而言,其求解過(guò)程是一個(gè)高度非線性的過(guò)程,異常值對(duì)最終求解結(jié)果的影響程度將會(huì)更大。當(dāng)異常值較大或者較多時(shí),會(huì)使算法性能退化,甚至出現(xiàn)不收斂的情況。

圖1 L2范數(shù)分布示意圖Fig.1 Distribution of L2 norm

2 基于最大相關(guān)熵準(zhǔn)則的壓縮感知地震道重建

2.1 非線性相關(guān)熵度量準(zhǔn)則

(6)

式中,κσ為帶寬為σ的核函數(shù);N為樣本總數(shù)。高斯核函數(shù)可以有效地將低維數(shù)據(jù)映射到高維空間[18],適合描述從低維采樣到高維數(shù)據(jù)恢復(fù)的過(guò)程。考慮高斯核函數(shù),式(6)變?yōu)?/p>

(7)

其中ri=yi-xi為恢復(fù)殘差。

令N=1,可以通過(guò)以下導(dǎo)數(shù)來(lái)評(píng)估式(6)對(duì)于不同殘差值ri的敏感度:

(8)

圖2為不同帶寬(w)所對(duì)應(yīng)的相關(guān)熵以及相關(guān)熵敏感度隨殘差值的變化趨勢(shì),可以看出,相關(guān)熵具有遠(yuǎn)端衰減的作用。對(duì)比圖1與圖2可以看出,異常值在L2范數(shù)計(jì)算過(guò)程中起重要作用,敏感度更高,而在相關(guān)熵計(jì)算中由于遠(yuǎn)端衰減的作用,其影響可得到一定程度的衰減。不同的帶寬對(duì)應(yīng)的相關(guān)熵及其敏感度曲線的作用范圍不同:帶寬越小,相關(guān)熵的有效計(jì)算范圍越小,其導(dǎo)數(shù)收斂趨向于零的速度越快,遠(yuǎn)端異常值的衰減越明顯。對(duì)于高斯核函數(shù),Liu等[19]指出合適的帶寬取值范圍為[0.8,3],并且在取值為1時(shí),固有誤差達(dá)到最小。本文中將帶寬固定為1。

圖2 相關(guān)熵分布示意圖Fig.2 Distribution of correntropy

2.2 基于最大相關(guān)熵準(zhǔn)則的壓縮感知地震道重建模型

在信息論中, 相關(guān)熵準(zhǔn)則可用來(lái)處理受到噪聲影響的信號(hào)分析[19-21]。在最大相關(guān)熵準(zhǔn)則下,可以通過(guò)最大化觀測(cè)值和恢復(fù)信號(hào)之間差值的相關(guān)熵來(lái)估量恢復(fù)結(jié)果。從圖2(a)可以看出,相關(guān)熵函數(shù)隨著帶寬σ的變化存在尺度問(wèn)題,為此對(duì)其進(jìn)行尺度歸一化處理。當(dāng)帶寬σ=1時(shí),歸一化后的最大相關(guān)熵JM(r)可以表示為

(9)

其中

r=y-RΩ(Φα).

常規(guī)的L1范數(shù)是求解問(wèn)題的凸最優(yōu)化極小值問(wèn)題。對(duì)于嚴(yán)格凸優(yōu)化問(wèn)題來(lái)說(shuō),局部最優(yōu)解就是全局最優(yōu)解。從相關(guān)熵的表達(dá)式及其示意圖(圖2(a))可以看出,相關(guān)熵是明顯的凹函數(shù)優(yōu)化問(wèn)題,難以直接利用凸函數(shù)優(yōu)化算法進(jìn)行求解。另一方面,同L2范數(shù)一樣,相關(guān)熵函數(shù)同樣滿足極值的唯一性,區(qū)別在于L2范數(shù)求取極小值,而相關(guān)熵函數(shù)求解極大值。一個(gè)簡(jiǎn)單直接的解決方案是對(duì)常規(guī)的相關(guān)熵函數(shù)乘以-1,便得到以相關(guān)熵為度量函數(shù)的凸函數(shù)極小值優(yōu)化目標(biāo)函數(shù)JCM(r)為

(10)

式(10)本質(zhì)上等價(jià)于求解最大化相關(guān)熵,因此仍將式(10)稱(chēng)為最大化相關(guān)熵。結(jié)合式(10)與式(3),可以得到基于最大化相關(guān)熵準(zhǔn)則的地震道重建模型

(11)

2.3 基于偽觀測(cè)數(shù)據(jù)的迭代閾值求解算法

雖然式(11)是一個(gè)凸優(yōu)化問(wèn)題,然而相關(guān)熵度量標(biāo)準(zhǔn)中指數(shù)函數(shù)的存在導(dǎo)致其求解過(guò)程更加復(fù)雜。為此,借鑒Zhao等[18]在解決非線性Huber函數(shù)計(jì)算時(shí)所引入的偽觀測(cè)數(shù)據(jù)的構(gòu)建方式,構(gòu)造如下第k次迭代時(shí)的偽觀測(cè)數(shù)據(jù)yP(αk):

yP(αk)=RΩΦαk-φ(y-RΩΦαk).

(12)

將相關(guān)熵的梯度表達(dá)式代入式(12)中,整理得到

yP(αk)=RΩΦαk+μ(y-RΩΦαk).

(13)

對(duì)于式(3)所表示的常規(guī)壓縮感知恢復(fù)模型,第k次迭代時(shí)關(guān)于殘差r的梯度為

(14)

將偽觀測(cè)數(shù)據(jù)yP替換式(13)中觀測(cè)數(shù)據(jù)y,并將式(8)表達(dá)式代入,得到

(15)

另外,基于最大化相關(guān)熵準(zhǔn)則的地震道重建模型在第k次迭代時(shí)關(guān)于殘差r的梯度為

(16)

比較式(15)及式(16)可以看出,當(dāng)基于L2范數(shù)的常規(guī)壓縮感知模型的輸入數(shù)據(jù)為偽觀測(cè)數(shù)據(jù)時(shí),其關(guān)于殘差的更新梯度等價(jià)于基于最大化相關(guān)熵準(zhǔn)則的地震道重建模型的梯度。利用偽觀測(cè)數(shù)據(jù),基于最大化相關(guān)熵準(zhǔn)則的地震道重建模型便可以借鑒式(4)的迭代閾值算法進(jìn)行求解,具體形式為

(17)

3 實(shí)驗(yàn)與方法對(duì)比

利用如圖3(a)所示的模擬地震數(shù)據(jù)進(jìn)行測(cè)試分析。該模型數(shù)據(jù)由有限差分方法對(duì)SEG/EAGE鹽丘模型進(jìn)行正演模擬得到。正演模擬子波為主頻30 Hz的雷克子波,空間及時(shí)間采樣步長(zhǎng)分別為10 m和1 ms。為驗(yàn)證所研究方法的重建效果,首先隨機(jī)地對(duì)地震記錄中40%地震道進(jìn)行充零,模擬地震道缺失的問(wèn)題;其次,對(duì)剩余的地震道添加隨機(jī)噪聲,并對(duì)20%的剩余地震道添加異常振幅噪聲,結(jié)果如圖3(b)所示?；贚2范數(shù)以及最大相關(guān)熵準(zhǔn)則的重建結(jié)果如圖3(c)、(d)所示,去除的噪聲剖面如圖3(e)、(f)所示。圖4為圖3(a)～(d)所對(duì)應(yīng)的f-k譜圖。對(duì)比兩種殘差度量的重建結(jié)果可以看出,基于L2范數(shù)的重建結(jié)果中缺失地震道基本得到恢復(fù),但重建結(jié)果受所添加異常振幅噪聲的影響明顯,存在局部振幅異常的現(xiàn)象,且去除的噪聲剖面中存在明顯的有效信號(hào)泄漏,頻譜圖中因噪聲及缺失道所引起的干擾去除不完全;基于最大相關(guān)熵準(zhǔn)則的重建模型對(duì)所添加的異常振幅噪聲具有更好的魯棒性,重建結(jié)果的信噪比更高,去除的噪聲剖面中有效信號(hào)相對(duì)更少。

圖3 模擬地震數(shù)據(jù)測(cè)試分析Fig.3 Test and analysis of simulated seismic data

圖4 模擬地震數(shù)據(jù)f-k譜圖Fig.4 The f-k spectrum of simulated seismic data

為了進(jìn)一步驗(yàn)證在不同噪聲含量及不同地震道缺失情況下的重建效果,分別固定噪聲含量或地震道缺失量,測(cè)試兩種方法的重構(gòu)質(zhì)量及收斂速度。測(cè)試結(jié)果如圖5所示,其中圖5(a)為固定地震道缺失比、不同異常噪聲含量的重構(gòu)質(zhì)量及收斂速度對(duì)比結(jié)果,圖5(b)為固定異常噪聲含量、不同地震道缺失的重構(gòu)質(zhì)量及收斂速度對(duì)比結(jié)果。可以看出,隨著噪聲含量與缺失道的增加,兩種方法均能逐漸收斂,驗(yàn)證了壓縮感知模型在隨機(jī)缺失地震道重建中的有效性。在初始迭代環(huán)節(jié),L2范數(shù)的收斂速度要快于最大相關(guān)熵準(zhǔn)則,這是因?yàn)樽铋_(kāi)始?xì)埐钪写嬖诖罅坑行盘?hào),高斯核函數(shù)除了對(duì)異常振幅噪聲具有壓制作用外,還間接地導(dǎo)致了有效信號(hào)能量的損失;隨著迭代次數(shù)的增加,殘差中的有效信號(hào)逐漸減少,高斯核對(duì)于有效信號(hào)的影響逐漸降低,而異常振幅噪聲仍可得到有效壓制。雖然最大相關(guān)熵準(zhǔn)則的初始收斂速度要慢于L2范數(shù),但在一定的迭代次數(shù)后,本文方法的重構(gòu)質(zhì)量要明顯高于L2范數(shù)。

圖5 高密集的噪聲以及欠采樣恢復(fù)收斂曲線Fig.5 Recovery convergence curve of high-density noise and under-sampling input

圖6為利用本文方法對(duì)本測(cè)試環(huán)節(jié)的最復(fù)雜輸入地震數(shù)據(jù)的測(cè)試結(jié)果。可以看出,對(duì)于高密集的異常振幅噪聲以及高度缺失的地震道數(shù)據(jù),本文方法仍然能實(shí)現(xiàn)魯棒重建。圖4(a)及圖6(d)所示的f-k譜圖對(duì)比表明,本文方法重建結(jié)果的高低頻信息均能得到較高質(zhì)量的恢復(fù)。由此可見(jiàn),本文方法對(duì)于嚴(yán)重異常振幅噪聲下的高缺失數(shù)據(jù)重構(gòu)具有較強(qiáng)的魯棒性。

為進(jìn)一步驗(yàn)證本文方法的有效性和實(shí)用性,利用如圖7(a)所示的某海上拖纜數(shù)據(jù)進(jìn)行實(shí)際資料測(cè)試驗(yàn)證。在該實(shí)際資料中受海浪波動(dòng)影響,約30%的檢波器脫離預(yù)定位置。對(duì)于脫離預(yù)定位置的檢波器,假定其采集信息缺失。除地震道缺失影響外,該拖纜數(shù)據(jù)存在明顯的單道異常振幅噪聲。分別對(duì)該地震數(shù)據(jù)應(yīng)用基于L2范數(shù)及最大相關(guān)熵準(zhǔn)則的重建方法進(jìn)行信號(hào)恢復(fù),恢復(fù)結(jié)果如圖7(b)、(c)所示?？梢钥闯?兩種方法均可以較好地重建缺失道位置的數(shù)據(jù),但L2范數(shù)的重建結(jié)果仍受異常振幅噪聲影響明顯;基于最大相關(guān)熵準(zhǔn)則的重建方法可較好地進(jìn)行異常振幅噪聲壓制以及缺失地震道重建。對(duì)比不同數(shù)據(jù)的f-k譜圖特征可以明顯看出,異常振幅噪聲的存在導(dǎo)致L2范數(shù)重建結(jié)果的f-k譜圖改善不明顯,而本文方法對(duì)于異常振幅噪聲的去除較徹底,可更好地恢復(fù)有效信號(hào)的頻譜信息。

圖7 實(shí)際地震數(shù)據(jù)測(cè)試分析Fig.7 Test and analysis of field seismic data

另外,該海上拖纜數(shù)據(jù)的深層信號(hào)能量較弱,為便于對(duì)比,圖8給出了圖7(a)～(c)中方框位置的信息。可以看出,涌浪的存在導(dǎo)致L2范數(shù)對(duì)信號(hào)的估計(jì)出現(xiàn)偏差,其求解結(jié)果不再是最優(yōu)解,重建結(jié)果上部分弱信號(hào)缺失;本文方法重建結(jié)果中弱信號(hào)能量清晰、連續(xù)性較好,在一定程度上體現(xiàn)了本文方法在異常振幅噪聲影響下的弱信號(hào)保護(hù)能力。

圖8 實(shí)際地震數(shù)據(jù)測(cè)試局部放大Fig.8 Magnified areas of field seismic data

4 結(jié)束語(yǔ)

在傳統(tǒng)壓縮感知地震信號(hào)重建方法的基礎(chǔ)上,通過(guò)引入信息論中的相關(guān)熵,構(gòu)建了基于最大相關(guān)熵準(zhǔn)則的壓縮感知地震道重建方法,實(shí)現(xiàn)對(duì)欠采樣地震記錄中異常振幅噪聲的有效衰減,從而保證地震道重建過(guò)程的穩(wěn)定性。實(shí)驗(yàn)結(jié)果表明,欠采樣數(shù)據(jù)中的異常振幅噪聲會(huì)嚴(yán)重影響傳統(tǒng)壓縮感知地震信號(hào)重建方法的重建精度,而最大相關(guān)熵準(zhǔn)則可有效抑制噪聲對(duì)重建過(guò)程的影響,可有效提高重建方法的精度及穩(wěn)定性,表明了本文方法在含異常振幅噪聲的非規(guī)則地震道缺失重建中的準(zhǔn)確性及有效性。