基于自適應(yīng)字典學(xué)習(xí)的可控震源數(shù)據(jù)諧波噪聲壓制方法

毛海波,馬俊彥,王曉濤,蔣 立,王麗麗,張有平,王惠迎,劉達(dá)偉

(1.中國石油天然氣股份有限公司新疆油田分公司勘探開發(fā)研究院地球物理研究所,新疆烏魯木齊830013;2.中國石油集團(tuán)東方地球物理勘探有限責(zé)任公司研究院,新疆烏魯木齊830002;3.中國石油集團(tuán)東方地球物理勘探有限責(zé)任公司新疆物探處,新疆烏魯木齊830002;4.西安交通大學(xué)電子與信息工程學(xué)院,陜西西安710049)

在地震勘探中,傳統(tǒng)的炸藥震源激發(fā)方式對環(huán)境具有不可逆轉(zhuǎn)的污染,且可控性差。與炸藥震源相比,可控震源不僅安全環(huán)保、施工效率高、成本低,而且組織靈活,可按照特定的深層地質(zhì)目標(biāo)和工區(qū)地表?xiàng)l件人為地調(diào)整所激發(fā)的信號,提升能量利用率?；诳煽卣鹪吹纳鲜鰞?yōu)點(diǎn),目前它已經(jīng)成為陸上地震勘探的主要激發(fā)方式之一,在國內(nèi)外得到了大規(guī)模的應(yīng)用[1]。

可控震源滑動掃描采集方式極大地提高了地震數(shù)據(jù)采集效率,覆蓋次數(shù)較炸藥震源提高幾倍乃至幾十倍,但是因滑動掃描方式采集的地震數(shù)據(jù)在不同炮記錄時間上部分重疊在一起,在相關(guān)處理后本炮數(shù)據(jù)會受到后一炮數(shù)據(jù)的諧波干擾。滑動掃描的滑動時間越短,采集到的地震記錄受諧波干擾就越嚴(yán)重。諧波噪聲降低了地震資料的信噪比和分辨率,對有效信號造成了干擾,進(jìn)而會影響震資料的成像與解釋。因此,高效的諧波噪聲壓制方法必然會受到關(guān)注。對于海量地震數(shù)據(jù),高效的諧波噪聲壓制技術(shù)也是其工業(yè)化應(yīng)用的前提。

按照地震數(shù)據(jù)的不同,可控震源地震記錄諧波噪聲壓制方法主要分為相關(guān)前和相關(guān)后兩大類。針對相關(guān)前的地震數(shù)據(jù)諧波噪聲壓制問題,LI等[2]對掃頻信號進(jìn)行分析得到純相移濾波算子,再對單道地震數(shù)據(jù)進(jìn)行純相移濾波處理壓制諧波噪聲;黃建平等[3]通過統(tǒng)計(jì)分析獲得統(tǒng)一的地面力信號,然后設(shè)計(jì)各諧波分量對應(yīng)的相移濾波器,該方法可以對多炮數(shù)據(jù)進(jìn)行處理,但是它只對相關(guān)前地震數(shù)據(jù)的高次諧波具有較好的壓制效果。由于相關(guān)前地震數(shù)據(jù)量較大,計(jì)算速度較慢,難以滿足實(shí)際生產(chǎn)需求,相關(guān)后的地震數(shù)據(jù)諧波噪聲壓制問題是目前主要的研究方向。YU等[4]提出了分頻異常振幅壓制方法,根據(jù)諧波噪聲與有效信號小波系數(shù)的差異,對地震數(shù)據(jù)有效信號的小波系數(shù)進(jìn)行閾值處理,進(jìn)而達(dá)到壓制諧波噪聲的目的,這種方法處理速度較快,但是當(dāng)有效信號與諧波噪聲的小波系數(shù)差異不明顯時,有效信號成分不能有效還原,從而造成較大損傷。SICKING等[5]提出地面力信號濾波法,該方法首先使用地面力信號分解得到基波和諧波,然后利用傅里葉變換將基波和諧波變換到頻率域,最后將諧波頻譜與基波頻譜的比值作為濾波器的系數(shù)來設(shè)計(jì)諧波噪聲濾波器。鐘飛等[6]將單道地震信號進(jìn)行Hilbert-Huang變換,在變換域設(shè)置閾值濾除諧波。WANG等[7]設(shè)計(jì)了一種諧波預(yù)測算子,在干擾區(qū)域成功將諧波噪聲剔除。然而,與常規(guī)方法相同,該方法必須從最后一炮地震記錄開始處理,無法進(jìn)行單炮處理。韓文功等[8]認(rèn)為諧波噪聲是諧波畸變信號與反射系數(shù)的褶積,諧波噪聲能夠由精確的反射系數(shù)和地面力信號估算出來。曲英銘等[9]首先采用諧波壓制濾波器基于最小二乘法對預(yù)測的諧波干擾進(jìn)行修正,以實(shí)現(xiàn)對力信號內(nèi)含諧波的準(zhǔn)確壓制,然后分頻濾波壓制諧波噪聲。羅勇等[10]將時頻域稀疏優(yōu)化諧波噪聲壓制方法應(yīng)用于準(zhǔn)噶爾盆地瑪湖地區(qū)三維典型單炮記錄諧波噪聲壓制,該方法具有良好的保真性。LI等[11]將時頻域稀疏優(yōu)化方法引入到諧波噪聲壓制問題中,該方法按照有效信號與諧波噪聲在時頻域的形態(tài)特征差異,通過分析選擇連續(xù)小波變換對有效信號進(jìn)行稀疏表示,選擇Chirplet變換對諧波噪聲進(jìn)行稀疏表示,組成可以用分塊坐標(biāo)松弛算法求解的超完備字典,進(jìn)而迭代求解得到有效信號和諧波噪聲。陳建友等[12]在LI等[11]提出的時頻域稀疏優(yōu)化方法的基礎(chǔ)上,應(yīng)用振幅譜比值衡量諧波噪聲強(qiáng)弱,可以自適應(yīng)地選擇參數(shù)壓制不同強(qiáng)度的諧波噪聲。時頻域稀疏優(yōu)化諧波噪聲壓制方法在實(shí)際應(yīng)用中取得了良好的效果,但考慮到地震資料的復(fù)雜性,采用固定字典的稀疏優(yōu)化方法,在有些區(qū)域仍然存在諧波噪聲壓制效果不理想及有效信號損傷的問題。同時,與傳統(tǒng)的固定字典稀疏表示方法相比,自適應(yīng)學(xué)習(xí)字典方法能更好地適應(yīng)復(fù)雜數(shù)據(jù)[13-14],進(jìn)而取得更好的諧波噪聲壓制效果。為此,本文提出了基于自適應(yīng)學(xué)習(xí)字典的諧波噪聲壓制方法,該方法直接從相關(guān)后地震炮集記錄中學(xué)習(xí)稀疏表示字典,根據(jù)諧波噪聲與有效信號的波形形態(tài)差異,將學(xué)習(xí)得到的字典原子分類為諧波噪聲子字典與有效信號子字典,分別稀疏表示諧波噪聲和有效信號。然后應(yīng)用分塊坐標(biāo)下降算法,利用分類得到的子字典分別重構(gòu)諧波噪聲與有效信號,實(shí)現(xiàn)自適應(yīng)地分離諧波噪聲的目的。

1 方法原理

1.1 建立模型

假設(shè)地震單道信號s為有效信號sr、諧波噪聲sh和隨機(jī)噪聲sn的線性疊加:

s=sr+sh+sn

(1)

其中,隨機(jī)噪聲sn可使得該模型具有更加普遍的適用性。有效信號sr與諧波噪聲sh的分離需要解決如下優(yōu)化問題:

(2)

式中:xr為有效信號字典Dr對應(yīng)的有效信號稀疏表示系數(shù)向量;xh為諧波噪聲字典Dh對應(yīng)的諧波噪聲稀疏表示系數(shù)向量;ε是誤差門限,目的是使該優(yōu)化問題在隨機(jī)噪聲或者其它干擾下保持穩(wěn)定。

公式(2)的求解關(guān)鍵是得到可以分別稀疏表示有效信號與諧波噪聲的字典,本文根據(jù)兩種信號的時間域波形特征差異,首先利用K-SVD得到超完備字典,然后對超完備字典進(jìn)行分類。由于諧波噪聲的能量呈現(xiàn)高頻端集中的特點(diǎn),為了最大程度地逼近諧波噪聲的波形特征,將諧波噪聲原子選擇為學(xué)習(xí)得到的超完備字典中能量集中在高頻端的原子。高、低頻分界點(diǎn)可以通過計(jì)算振幅譜比確定,設(shè)置振幅譜比的閾值將超完備字典分類為有效噪聲的子字典Dr與諧波信號的子字典Dh。最后,利用分類所得的兩個子字典分別稀疏重構(gòu)諧波噪聲與有效信號,實(shí)現(xiàn)分離諧波噪聲的目的,其技術(shù)流程如圖1所示。

圖1 本文算法流程

1.2 字典學(xué)習(xí)

傅里葉變換[15]、小波變換[16]、Ridgelet變換[17]、Curvelet變換[18]、Shearlet變換[19]等使用固定的波形字典,因此無法完全匹配地震信號的復(fù)雜特征。與固定字典相比,學(xué)習(xí)字典可以自適應(yīng)地學(xué)習(xí)信號特征以匹配復(fù)雜信號的特征。ENGAN等[20]提出的最優(yōu)方向字典學(xué)習(xí)算法(MOD)需要計(jì)算矩陣的偽逆,導(dǎo)致算法具有很高的復(fù)雜度,實(shí)用性不強(qiáng)。AHARON等[21]提出的K-SVD字典學(xué)習(xí)算法逐列更新字典原子,克服了MOD算法需要計(jì)算矩陣偽逆的缺陷,有效地提高了字典學(xué)習(xí)的效率。

K-SVD算法通過逐列更新字典的方法來更新字典A,其目標(biāo)方程可以表示為:

(3)

求解公式(3),首先利用初始固定字典A求得信號的稀疏表示系數(shù)矩陣x,然后再根據(jù)系數(shù)矩陣x通過逐列更新字典原子獲得最終的字典A。

首先,使用正交匹配追蹤算法(OMP)[22]求解稀疏系數(shù)。MALLAT等[23]提出匹配追蹤算法(MP)求解信號的稀疏表示,但是MP算法是純貪婪算法,計(jì)算復(fù)雜度大,無法保證每次迭代的投影是正交投影,最終求解結(jié)果可能并不稀疏。PATI等[22]提出的OMP算法,在MP算法基礎(chǔ)上每次迭代后進(jìn)行正交化處理,克服了MP算法的缺點(diǎn)。信號y經(jīng)過l次分解后如下式:

(4)

式中:ai為第i次分解的原子;xi為第i次分解時的系數(shù);Rly為l次分解后的殘差。

為了滿足正交條件,OMP算法需要由重建原子集合中利用最小二乘算法得到第l次分解的稀疏表示系數(shù):

x′l=argmin‖y-Alx‖

(5)

然后更新殘差:

Rly=y-Ax′ll=l+1

(6)

使得OMP算法殘差逐漸減小,直到算法收斂。

字典更新時首先固定系數(shù)矩陣x和字典A,更新字典的第k列ak時,令x中字典原子ak對應(yīng)第k行為xk,則字典學(xué)習(xí)的目標(biāo)函數(shù)可改寫為:

(7)

式中:Ek為去掉原子ak影響的誤差。將Ax分解成L個秩為1的矩陣的和。逐列更新字典的第k列原子時,固定其它的L-1列。

利用奇異值分解(SVD)方法處理字典矩陣A以更新ak和xk,但是經(jīng)過SVD之后,xk將會滿向量,即SVD得到的xk中非零值的數(shù)量和位置會和原xk不同,al將不能滿足稀疏表示的條件。

因此僅保留xk中非零值(去掉xk中所有的零值),再應(yīng)用SVD更新al和xk,(7)式轉(zhuǎn)化為下式:

(8)

1.3 字典分類

本文通過計(jì)算每個超完備字典原子振幅譜比,選定一個高、低頻分界點(diǎn)閾值,將字典原子分類為諧波噪聲子字典與有效信號子字典,用于分別稀疏表示諧波噪聲與有效信號。高、低頻分界點(diǎn)的閾值選擇是一個經(jīng)驗(yàn)值,主要評估信噪分離后有效信號的保真程度及壓制掉的噪聲強(qiáng)度。如果這個閾值選擇過小,壓制的噪聲能量大,有可能會損傷一定的有效信號;如果選擇過大,則有效信號保真度高,但是噪聲殘余多。

字典原子振幅譜比計(jì)算步驟如下。

1) 確定高、低頻的分界點(diǎn)p,計(jì)算分界點(diǎn)以上的高頻帶能量與總頻帶能量,地震信號的頻率一般不會超過100Hz,所以計(jì)算能量時將100Hz以上的能量舍棄以減少運(yùn)算量,算式如下。

(9)

(10)

式中:D′[k]為字典原子Fourier振幅譜;Efh為高頻能量;Ef為總頻帶能量。

2) 計(jì)算振幅譜比值α,即分界點(diǎn)p以上的高頻帶能量與總頻帶能量的比值:

α=Efh/Ef

(11)

2 算例分析

2.1 合成資料

圖2為合成地震數(shù)據(jù),由有效信號與諧波噪聲疊加得到。該數(shù)據(jù)包含3個反射層,利用卷積模型生成,所用掃描信號為線性升頻掃描信號:

圖2 含諧波噪聲合成地震數(shù)據(jù)

(12)

掃描信號的最低頻率為3Hz,最高頻率為90Hz,諧波噪聲包含二次諧波和三次諧波。該合成地震數(shù)據(jù)共有300道,采樣點(diǎn)為3000,采樣間隔為2ms。

首先應(yīng)用K-SVD字典學(xué)習(xí)方法在樣本集上進(jìn)行字典學(xué)習(xí),得到超完備字典如圖3a所示,然后計(jì)算字典原子的振幅譜比,選擇40Hz作為計(jì)算高低頻能量的分界頻率點(diǎn),即利用40～100Hz頻率范圍計(jì)算高頻帶能量,0～100Hz頻率范圍計(jì)算總頻帶能量。圖3b為超完備字典原子振幅譜比計(jì)算結(jié)果,其中紅線值為0.45,是設(shè)置的閾值,字典原子振幅譜比值高于閾值的為諧波噪聲字典原子,低于閾值的為有效信號字典原子,分類得到的有效信號字典和諧波噪聲字典如圖4所示。

圖3 合成數(shù)據(jù)字典訓(xùn)練(a)與K-SVD訓(xùn)練得到的超完備字典的字典原子振幅譜比(b)

為了驗(yàn)證本文方法分類得到的有效信號子字典與諧波噪聲子字典對有效信號與諧波噪聲表示的稀疏性,隨機(jī)提取一段用于構(gòu)建圖2合成數(shù)據(jù)的有效信號(如圖5),計(jì)算圖4中有效信號子字典與諧波噪聲子字典對圖5所示有效信號的表示系數(shù),結(jié)果見圖6。對比可得,本文分類得到的有效信號子字典相比于諧波噪聲子字典對原始有效信號的表示更稀疏,僅用5個有效信號的原子即可實(shí)現(xiàn)對有效信號的表示,而采用諧波噪聲子字典對有效信號進(jìn)行表示時需要大約30個原子。進(jìn)一步隨機(jī)提取一段用于構(gòu)建圖2 合成數(shù)據(jù)的諧波噪聲(圖7),計(jì)算圖4中有效信號子字典與諧波噪聲子字典對其的表示系數(shù),結(jié)果如圖8所示。對比可得,本文分類得到的諧波噪聲子字典相比于有效信號子字典對諧波噪聲的表示更稀疏。上述結(jié)果說明了本文提出的分類方法得到的有效信號子字典與諧波噪聲子字典都更能夠稀疏表示各自的信號,而對另一種信號的表示更不稀疏,滿足了信號分離的條件。

圖8 有效信號字典(a)和諧波噪聲字典(b)分別對諧波噪聲的稀疏表示系數(shù)

圖4 合成數(shù)據(jù)有效信號字典(a)與諧波噪聲字典(b)

圖5 選取的有效信號

圖6 有效信號字典(a)和諧波噪聲字典(b)分別對有效信號的稀疏表示系數(shù)

圖7 諧波噪聲

應(yīng)用分類得到的子字典分別重構(gòu)合成信號的有效信號與諧波噪聲如圖9所示。對比圖2a合成的有效信號與圖9a使用本文方法壓制諧波噪聲得到的有效信號可見,本文方法能夠高保真地恢復(fù)合成的有效信號,并且有效壓制諧波噪聲。圖9b為本文方法得到的諧波噪聲剖面,可以看到,噪聲剖面幾乎沒有反射波的損傷,只對直達(dá)波有少許損傷。因此,本文方法成功地壓制了諧波噪聲。

圖9 合成地震數(shù)據(jù)處理結(jié)果

2.2 實(shí)際數(shù)據(jù)

含有諧波噪聲的實(shí)際地震數(shù)據(jù)如圖10所示,每炮數(shù)據(jù)共400道,采樣間隔為2ms,采樣點(diǎn)為3501,記錄長度為7s。圖11是實(shí)際地震數(shù)據(jù)第210道信號的時頻圖,圖中諧波噪聲與有效信號相比,主要表現(xiàn)為高頻,有效信號相對于諧波噪聲主要表現(xiàn)為低頻。

圖10 原始地震數(shù)據(jù)

圖11 實(shí)際地震數(shù)據(jù)第210道數(shù)據(jù)時頻圖

首先進(jìn)行單道循環(huán)滑動截取組成樣本數(shù)據(jù)集,截取長度為300采樣點(diǎn),每隔10個采樣點(diǎn)截取一個樣本。應(yīng)用K-SVD算法在樣本集上訓(xùn)練,初始字典選擇離散余弦變換(DCT),迭代次數(shù)為50,字典原子為300采樣點(diǎn),冗余度為10,得到的超完備字典如圖12所示。然后計(jì)算字典原子的振幅譜比,選擇40Hz作為計(jì)算高頻能量的分界頻率點(diǎn),即使用40～100Hz頻率范圍計(jì)算高頻帶能量,0～100Hz頻率范圍計(jì)算總頻帶能量。圖13為振幅譜比計(jì)算結(jié)果,其中灰色直線值為0.40是設(shè)置的閾值,字典原子振幅譜比值高于閾值的為諧波噪聲字典原子,低于閾值的為有效信號字典原子,將超完備字典分類為有效信號字典(圖14a)與諧波噪聲字典(圖14b)。需要指出的是,這里的閾值0.40是一個經(jīng)驗(yàn)值,在其附近都可以實(shí)現(xiàn)諧波噪聲壓制。本文選擇閾值0.40主要用于評估信噪分離后有效信號的保真程度及壓制掉的噪聲強(qiáng)度。

圖12 實(shí)際地震數(shù)據(jù)訓(xùn)練得到的超完備字典

圖13 實(shí)際地震數(shù)據(jù)訓(xùn)練得到的字典原子振幅譜比

圖14 實(shí)際地震數(shù)據(jù)訓(xùn)練得到的有效信號字典(a)和諧波噪聲字典(b)

采用本文方法重構(gòu)得到的有效信號與諧波噪聲分別如圖15a、圖15b所示。對比圖10原始地震資料與采用本文方法得到的有效信號圖15a可見,原始地震資料中的諧波噪聲得到了有效壓制,且對有效信號的損傷較小,結(jié)果驗(yàn)證了本文方法壓制諧波噪聲的有效性。同樣,采用自適應(yīng)稀疏優(yōu)化方法[12]對圖10數(shù)據(jù)進(jìn)行處理,得到的結(jié)果如圖16所示,其中,圖16a 為有效信號,圖16b為諧波噪聲。對比圖15b與圖16b可知,本文方法對近炮點(diǎn)數(shù)據(jù)有些許損傷,對遠(yuǎn)炮點(diǎn)數(shù)據(jù)幾乎沒有損傷,而自適應(yīng)稀疏優(yōu)化方法幾乎對整炮數(shù)據(jù)均有損傷。

圖15 采用本文方法得到的有效信號(a)和諧波噪聲(b)

為了進(jìn)一步觀察,我們提取近炮點(diǎn)強(qiáng)噪聲區(qū)域第220道數(shù)據(jù)與遠(yuǎn)炮點(diǎn)幾乎無噪聲的第60道數(shù)據(jù)進(jìn)行分析。圖17和圖18分別為近炮點(diǎn)第220道地震數(shù)據(jù)本文方法和稀疏優(yōu)化方法處理對應(yīng)的原始數(shù)據(jù)、有效信號、諧波噪聲。圖19和圖20分別為遠(yuǎn)炮點(diǎn)第60道地震數(shù)據(jù)本文方法和稀疏優(yōu)化方法處理對應(yīng)的原始數(shù)據(jù)、有效信號、諧波噪聲。對比圖17與圖18諧波噪聲(紅色方框區(qū)域有效信號損傷)可見,在近炮點(diǎn)強(qiáng)噪聲區(qū)域,本文方法對有效信號損傷比自適應(yīng)稀疏優(yōu)化方法小。對比圖19與圖20諧波噪聲(紅色方框區(qū)域有效信號損傷)可見,在遠(yuǎn)炮點(diǎn)區(qū)域本文方法對有效信號幾乎沒有損傷,而自適應(yīng)稀疏優(yōu)化方法對有效信號的損傷很大。綜合分析認(rèn)為本文方法雖然在強(qiáng)噪聲區(qū)域淺層有少量噪聲殘留,但對整個剖面而言,對有效信號損傷較小。因此,本文方法有效地壓制了諧波噪聲。

圖17 近炮點(diǎn)地震道(第220道)本文方法結(jié)果

圖18 近炮點(diǎn)地震道(第220道)自適應(yīng)稀疏優(yōu)化方法結(jié)果

圖19 遠(yuǎn)炮點(diǎn)地震道(第60道)本文方法結(jié)果

圖20 遠(yuǎn)炮點(diǎn)地震道(第60道)自適應(yīng)稀疏優(yōu)化方法結(jié)果

3 結(jié)論與討論

本文提出了一種基于自適應(yīng)學(xué)習(xí)字典的諧波噪聲壓制方法?；谟行盘柲芰恐饕杏诘皖l端與諧波噪聲能量主要集中于高頻端的頻率特征,提出了應(yīng)用字典原子振幅譜比將K-SVD方法學(xué)習(xí)得到的超完備字典分類為有效信號稀疏表示子字典和諧波噪聲稀疏表示子字典,最后利用兩個子字典分別重構(gòu)有效信號與諧波噪聲以達(dá)到壓制諧波噪聲的目的。合成地震數(shù)據(jù)及實(shí)際數(shù)據(jù)的應(yīng)用結(jié)果驗(yàn)證了本文方法的可行性。對比本文方法與自適應(yīng)稀疏優(yōu)化方法的應(yīng)用結(jié)果,表明本文方法壓制了絕大部分的噪聲,同時無論在近炮點(diǎn)還是遠(yuǎn)炮點(diǎn)對有效信號的損傷都小于自適應(yīng)稀疏優(yōu)化方法,說明了本文方法的有效性。

同時,本文方法還有進(jìn)一步提升的空間。首先,字典分類方法是本文方法的核心,本文僅使用振幅譜比單一指標(biāo)作為字典分類的依據(jù),閾值的選擇會影響字典分類的準(zhǔn)確性,尤其是對選定閾值附近的字典很難準(zhǔn)確區(qū)分。而通過多個指標(biāo)綜合決定字典的分類,會提升字典分類的準(zhǔn)確性,進(jìn)而提升諧波噪聲壓制的效果,提出多個字典分類方法,并綜合應(yīng)用于字典分類是本文方法未來的一個研究方向;其次,K-SVD方法獲得的字典原子沒有明確的物理含義,若能加上一些人為約束,使得獲得的字典原子具有明確的物理含義,這會提升對字典的理解,字典分類的準(zhǔn)確性會有較大提升,這也是本文方法未來的研究方向。