馬繼濤,王建花,劉國昌

(1.中國石油大學(xué)(北京)地球物理與信息工程學(xué)院物探系,北京102249;2.中海油研究總院,海洋石油勘探國家工程實(shí)驗室,北京100028)

?

基于頻率域奇異值分解的地震數(shù)據(jù)插值去噪方法研究

馬繼濤1,王建花2,劉國昌1

(1.中國石油大學(xué)(北京)地球物理與信息工程學(xué)院物探系,北京102249;2.中海油研究總院,海洋石油勘探國家工程實(shí)驗室,北京100028)

摘要:奇異值分解(Singular Value Decomposition,SVD)是近幾年發(fā)展起來的一種地震數(shù)據(jù)隨機(jī)噪聲壓制方法?；陬l率域奇異值分解矩陣降秩運(yùn)算,利用凸集投影(Projection Onto Convex Sets,POCS)迭代算法,實(shí)現(xiàn)了地震數(shù)據(jù)去噪和插值的同步處理,給出了方法的實(shí)現(xiàn)步驟。實(shí)際資料處理時，采用分窗處理方式減少了算法對內(nèi)存的需求,降低了插值和去噪處理的運(yùn)算量,同時使有效信號的同相軸線性化,滿足方法的假定條件。模擬數(shù)據(jù)和實(shí)際資料測試均表明,頻率域奇異值分解方法可以在壓制地震數(shù)據(jù)噪聲的同時進(jìn)行插值處理,具有廣闊的應(yīng)用前景。

關(guān)鍵詞：奇異值分解(SVD);頻率域;插值;去噪

奇異值分解(SVD)是提高地震資料信噪比的有效手段之一,被廣泛應(yīng)用于VSP波場分離、地震信號增強(qiáng)等。SERGIO等[1]首先提出基于K-L變換的奇異值分解方法,通過特征圖像濾波分離VSP波場。高靜懷等[2]、徐伯勛等[3]指出,從圖像處理的角度看,奇異值分解相當(dāng)于K-L變換。針對奇異值分解僅能提高水平同相軸信噪比、不能有效確定奇異值分解濾波參數(shù)等問題,陳遵德等[4]、陸文凱等[5]對SVD濾波方法做了一些改進(jìn)。吳亞東等[6]改進(jìn)了奇異值分解濾波方法,使其對脈沖干擾、側(cè)反射以及其它反向干擾的壓制效果更好。李文杰等[7]提出利用奇異值分解濾波法衰減地震記錄中的直達(dá)波和折射波。2009年,SACCHI[8]研究了基于奇異譜分析(頻率域奇異值分解)的地震數(shù)據(jù)隨機(jī)噪聲壓制技術(shù)。魏小強(qiáng)等[9]基于多道奇異譜分析理論,對三維地震數(shù)據(jù)進(jìn)行了插值重建,但是沒有涉及地震數(shù)據(jù)去噪問題。黃建平等[10]基于奇異譜分析理論對三維地震數(shù)據(jù)進(jìn)行去噪和插值,但未對反對角線平均的問題加以闡述。

頻率域奇異值分解壓制地震數(shù)據(jù)隨機(jī)噪聲的原理相對簡單:如果不含噪聲的地震數(shù)據(jù)由k個線性同相軸組成,則其對應(yīng)頻率域的Hankel矩陣的秩應(yīng)為k;若地震數(shù)據(jù)中存在隨機(jī)噪聲,則該矩陣的秩會增加。因此,通過降秩處理可以壓制地震數(shù)據(jù)中的隨機(jī)噪聲。同樣,地震數(shù)據(jù)有缺失道也會使Hankel矩陣的秩增加,因此,通過降秩處理還可以對地震數(shù)據(jù)進(jìn)行插值。本文基于頻率域奇異值分解理論,研究了地震數(shù)據(jù)插值和去噪同步處理方法,并對奇異值分解降秩處理后的反對角線平均問題做了詳細(xì)闡述。理論數(shù)據(jù)和實(shí)際資料的測試結(jié)果表明方法具有可行性與有效性。

1方法原理

1.1基于頻率域奇異值分解的插值去噪原理

(1)一般取M=E(N/2)+1,該矩陣反對角線上的元素相同。若某一線性同相軸在時間域為:

(2)式中:t為時間,x為偏移距,p為線性同相軸的斜率,w為地震子波。經(jīng)過傅里葉變換,(2)式在頻率域可表示為:

(3)對于規(guī)則排列的地震數(shù)據(jù),偏移距x=nΔx,其中n為地震道序號。令α=ωpΔx,則公式(3)又可表示為:

(4)簡寫為:

(5)由Fn和Fn-1可以得到信號的一階差分方程預(yù)測形式:

(6)將F1～Fn代入矩陣H中(公式(1)),此時矩陣H的秩為1。同理,當(dāng)?shù)卣鹌拭嬷泻琸個線性同相軸時,存在方程:

(7)此時,矩陣H的秩為k。因此,可以用前k個奇異值對地震數(shù)據(jù)中的有效信號進(jìn)行重構(gòu),從而壓制不相關(guān)的隨機(jī)噪聲。

由上述推導(dǎo)可以看出,若地震記錄中只含有一個線性同相軸且不含噪聲,則此地震記錄對應(yīng)Hankel矩陣的秩為1;同樣,若地震記錄含有k個線性同相軸且不含噪聲,對應(yīng)Hankel矩陣的秩為k。地震數(shù)據(jù)中存在的隨機(jī)噪聲會使Hankel矩陣的秩增加,利用奇異值分解對Hankel矩陣做降秩處理可以濾除隨機(jī)噪聲。類似地,地震數(shù)據(jù)若有地震道缺失,則缺失的地震道在數(shù)據(jù)矩陣中對應(yīng)列的數(shù)值為0。將缺失的地震道變換到f-x域時,其對應(yīng)頻率域的樣點(diǎn)數(shù)值仍然為0。這些0值樣點(diǎn)變換到Hankel矩陣中后,會使矩陣的秩增加。因此,地震數(shù)據(jù)缺失時,降秩處理的目的是得到一個與原Hankel矩陣類似、沒有缺失道的低秩矩陣。

利用頻率域奇異值分解對地震記錄中的缺失道進(jìn)行插值與隨機(jī)噪聲壓制的過程是類似的,因此兩種處理可以同時進(jìn)行。本文隨機(jī)噪聲壓制和缺失地震道插值的步驟如下:

1) 通過傅里葉變換將地震數(shù)據(jù)轉(zhuǎn)換至f-x域;

2) 將每個頻率分量轉(zhuǎn)變?yōu)镠ankel矩陣;

3) 通過奇異值分解,得到奇異值譜;

4) 對Hankel矩陣進(jìn)行降秩處理;

5) 沿Hankel矩陣反對角線方向求取平均值并轉(zhuǎn)換為頻率域地震道數(shù)據(jù);

6) 通過傅里葉反變換將頻率域地震道數(shù)據(jù)轉(zhuǎn)換至?xí)r間域,完成去噪、插值處理。

1.2反對角線平均處理

在利用頻率域奇異值分解進(jìn)行去噪、插值的過程中,若降秩處理后的矩陣仍然是Hankel矩陣,則將Hankel矩陣轉(zhuǎn)變?yōu)轭l率域地震道數(shù)據(jù)之后,通過傅里葉反變換就可以得到去噪和插值后的時間域數(shù)據(jù)。假定s為降秩后的頻率域數(shù)據(jù),那么s(n)可以從降秩后Hankel矩陣H′的反對角線元素得到,反對角線上矩陣元素的序號i,j滿足i+j-1=n。如s(1)可以由降秩后Hankel矩陣的元素H′(1,1)得到,s(2)可以由降秩后Hankel矩陣的元素H′(1,2),H′(2,1)得到。但在實(shí)際數(shù)據(jù)處理中,降秩后得到的矩陣一般不再是Hankel矩陣,需要沿反對角線方向求取矩陣元素平均值,以得到規(guī)則的Hankel矩陣。GOLYANDINA等[11]引入了一個能在降秩后求取矩陣反對角線元素平均值的算子。為簡化對算子的描述,假設(shè)Hankel矩陣H(L,K)中,L≤K,令i+j-1=n,J=L+K-1,則降秩后的頻率域數(shù)據(jù)s(n)可以由以下公式計算得到:

(8)

利用公式(8),可由降秩處理后的矩陣恢復(fù)去噪和插值后的頻率域數(shù)據(jù)s。

1.3POCS迭代算法

若只對地震數(shù)據(jù)進(jìn)行去噪處理,則利用公式(8) 恢復(fù)頻率域數(shù)據(jù),并通過傅里葉反變換至?xí)r間域后,可得到去噪后的地震數(shù)據(jù);若對地震數(shù)據(jù)去噪的同時進(jìn)行插值處理,則利用公式(8)恢復(fù)頻率域數(shù)據(jù),并通過傅里葉反變換至?xí)r間域后,會在空地震道處出現(xiàn)一些降低矩陣秩的數(shù)值,這些數(shù)值比其它地震道振幅要小。為得到空地震道處的正確振幅,本文利用POCS算法進(jìn)行迭代處理[12]。在f-x域進(jìn)行POCS迭代處理,可以在恢復(fù)空地震道正確振幅的同時,保持原有數(shù)據(jù)的振幅特性。

假定輸入數(shù)據(jù)S是二維數(shù)據(jù),建立一算子T,其作用是在插值后識別空地震道的位置,并提取插值前正常的地震道,將其放在原處。正常地震道處T(i,j)=1,空地震道處T(i,j)=0。將算子T作用于數(shù)據(jù)S,得到數(shù)據(jù)Sobs,即Sobs=T·S?？梢钥闯?Sobs=T·Sobs。

2.3.5 非心臟外科手術(shù)：為減少外科手術(shù)圍術(shù)期心臟并發(fā)癥風(fēng)險，在術(shù)前應(yīng)首先評估外科手術(shù)的緊迫性、出血風(fēng)險和心血管事件的發(fā)生風(fēng)險。在充分權(quán)衡出血和血栓風(fēng)險的基礎(chǔ)上，圍手術(shù)期抗血小板治療應(yīng)由多學(xué)科(外科醫(yī)師、麻醉師、心內(nèi)科醫(yī)生)和患者共同決定：出血危險較低的患者，可繼續(xù)服用阿司匹林。如患者進(jìn)行小型牙科手術(shù)、皮膚科操作、白內(nèi)障手術(shù)等出血風(fēng)險低的手術(shù)；手術(shù)相關(guān)出血風(fēng)險高，應(yīng)術(shù)前停用抗血小板藥物，通常術(shù)前停用P2Y12受體拮抗劑至少5d，術(shù)前需停用所有抗血小板治療的患者，如遇到血栓風(fēng)險高患者，可給予靜脈抗血小板藥物GPⅡb/Ⅲa 受體拮抗劑或低分子肝素“橋接”。

對重建后的數(shù)據(jù)再次進(jìn)行頻率域奇異值分解降秩處理,并將原正常的地震道再次放回原處,如此迭代直至缺失地震道振幅滿足設(shè)定的條件為止。假定I是與T大小相同且元素全為1的矩陣,則利用POCS迭代和奇異譜分析進(jìn)行去噪、插值的算法如下:

forp=1∶N

(9)

end

end

因?qū)嶋H地震資料數(shù)據(jù)量較大,利用本算法進(jìn)行插值和去噪時,奇異值分解對內(nèi)存的需求及計算量均非常大;同時實(shí)際資料并不完全滿足本方法假定有效信號是線性且平穩(wěn)信號的條件。為減少算法對內(nèi)存的需求,降低插值和去噪的運(yùn)算量,同時使有效信號的同相軸線性化,本文對數(shù)據(jù)進(jìn)行分窗處理。一般而言,分窗可以加快運(yùn)算速度,使地震數(shù)據(jù)的同相軸局部線性化和平穩(wěn)化,但窗口選取過大,會影響信號的局部線性化和平穩(wěn)化,使窗口內(nèi)大部分同相軸不滿足方法的假定條件,從而導(dǎo)致處理效果變差。窗口選取過小,很明顯會增加運(yùn)算量。因此,窗口的空間寬度和時間長度選取要適中,需根據(jù)實(shí)際數(shù)據(jù)的復(fù)雜程度進(jìn)行最優(yōu)化選擇。

2模擬數(shù)據(jù)測試

2.1有效性驗證

圖1　模擬數(shù)據(jù)迭代6次插值去噪處理前后對比(信噪比為5)a 處理前的數(shù)據(jù); b 插值去噪后的結(jié)果; c 壓制掉的噪聲

圖2　模擬數(shù)據(jù)迭代處理過程中的奇異值分布a 第1次迭代; b 第2次迭代; c 第3次迭代; d 第4次迭代; e 第5次迭代; f 第6次迭代

2.2抗噪性分析

為測試本文插值去噪方法的適用性,對圖1所示模擬數(shù)據(jù)分別加入一定程度的隨機(jī)噪聲,加噪后數(shù)據(jù)的信噪比分別為2,1(有效信號最大振幅與隨機(jī)噪聲最大振幅比值分別為2,1)。利用本文方法對其進(jìn)行插值去噪處理,同樣迭代6次,信噪比為2的數(shù)據(jù)處理結(jié)果如圖3所示,信噪比為1的數(shù)據(jù)處理結(jié)果如圖4所示。從圖3和圖4可以看出,即使在數(shù)據(jù)信噪比為1的情況下,本文方法仍然能夠在重建出空缺地震道的同時,有效壓制隨機(jī)噪聲。但同時可以看出,數(shù)據(jù)中噪聲能量較強(qiáng)時(圖4),去噪并重建后的結(jié)果中仍殘存部分隨機(jī)噪聲。進(jìn)一步測試表明,增加迭代次數(shù)可改善噪聲壓制效果。因為增加迭代次數(shù)可進(jìn)一步降低矩陣的秩,進(jìn)一步壓制隨機(jī)噪聲。

2.3重建效果分析

圖3　模擬數(shù)據(jù)迭代6次插值去噪處理前后對比(信噪比為2)a 處理前的數(shù)據(jù); b 插值去噪后的結(jié)果; c 壓制掉的噪聲

圖4　模擬數(shù)據(jù)迭代6次插值去噪處理前后對比(信噪比為1)a 處理前的數(shù)據(jù); b 插值去噪后的結(jié)果; c 壓制掉的噪聲

圖5　隨機(jī)空缺20道的模擬數(shù)據(jù)迭代10次插值去噪處理前后對比(信噪比為5)a 處理前的數(shù)據(jù); b 插值去噪后的結(jié)果; c 壓制掉的噪聲

圖6　隨機(jī)空缺40道的模擬數(shù)據(jù)迭代10次插值去噪處理前后對比(信噪比為5) a 處理前的數(shù)據(jù); b 插值去噪后的結(jié)果; c 壓制掉的噪聲

圖7　隨機(jī)空缺60道的模擬數(shù)據(jù)迭代10次插值去噪處理前后對比(信噪比為5)a 處理前的數(shù)據(jù); b 插值去噪后的結(jié)果; c 壓制掉的噪聲

對含噪聲且地震道缺失程度不同的數(shù)據(jù)進(jìn)行重建和去噪效果對比測試。模擬數(shù)據(jù)總道數(shù)為80,對其加入一定量的隨機(jī)噪聲,并隨機(jī)刪除其中的20,40和60道,最大連續(xù)空道個數(shù)為7道,最小連續(xù)空道個數(shù)為1道。利用前述頻率域奇異值分解插值去噪方法對其進(jìn)行插值去噪處理,迭代次數(shù)均為10次,處理結(jié)果分別如圖5,圖6和圖7所示?？梢钥闯?即使在地震道缺失程度高達(dá)75%的情況下,本文方法仍然可以較好地完成對地震數(shù)據(jù)的重建和去噪處理。然而從圖7c可以明顯看出,空缺60道的模擬數(shù)據(jù)插值去噪結(jié)果與原始數(shù)據(jù)之差存在一些有效波的信息,說明重建后數(shù)據(jù)與原始數(shù)據(jù)之間存在一定程度的差異。將圖7中的模擬數(shù)據(jù)迭代次數(shù)增加到30次后,差異剖面中基本看不到有效波信息(圖8),說明增加迭代次數(shù)可以改善最終的插值去噪處理效果。

圖8　隨機(jī)空缺60道的模擬數(shù)據(jù)迭代30次插值去噪處理前后對比(信噪比為5)a 處理前的數(shù)據(jù); b 插值去噪后的結(jié)果; c 壓制掉的噪聲

3實(shí)際資料測試處理

用南海某區(qū)實(shí)際資料對本文方法進(jìn)行了測試。考慮實(shí)際資料數(shù)據(jù)量大,采用滑動分窗處理方式,每個窗口的大小相同,為橫向100個樣點(diǎn),縱向100個樣點(diǎn)。在每個窗口內(nèi)同相軸可近似視為線性的。數(shù)據(jù)除自身存在著隨機(jī)噪聲外,還被隨機(jī)刪除了50道,刪除的地震道中,最大連續(xù)空道個數(shù)為8道。圖9a為缺失50道的地震數(shù)據(jù),圖9b為應(yīng)用本文方法迭代8次插值去噪后的結(jié)果,圖9c為去除的隨機(jī)噪聲,運(yùn)算中選取奇異值的個數(shù)為12個。從圖9b和圖9c可以看出,本文方法較好地重建了缺失的地震道,壓制的噪聲除了隨機(jī)噪聲之外,還有部分繞射噪聲。為了更好地展示本文方法的插值去噪效果,對圖9中的局部區(qū)域進(jìn)行放大顯

示。圖10展示了圖9中虛線框內(nèi)數(shù)據(jù)的重建及去噪結(jié)果,圖11展示了圖9中實(shí)線框內(nèi)數(shù)據(jù)的重建及去噪結(jié)果。由圖10和圖11可以看出,重建后的信號同相軸連續(xù),并且具有更高的信噪比。

圖9　實(shí)際地震資料插值去噪處理前后對比a 處理前的剖面; b 插值去噪后的結(jié)果; c 壓制掉的噪聲

圖10　實(shí)際地震剖面局部放大(針對虛線框)a 處理前的剖面; b 插值去噪后的結(jié)果

圖11　實(shí)際地震剖面局部放大(針對實(shí)線框)a 處理前的剖面; b 插值去噪后的結(jié)果

4結(jié)論與建議

本文基于頻率域奇異值分解實(shí)現(xiàn)了地震數(shù)據(jù)插值和去噪的同步處理,并用模擬數(shù)據(jù)和實(shí)際資料對方法進(jìn)行了驗證。通過分析和研究,得到以下結(jié)論與建議:

1) 噪聲的存在和地震道的缺失會影響地震數(shù)據(jù)中奇異值的分布,在頻率域?qū)?shù)據(jù)矩陣進(jìn)行降秩可實(shí)現(xiàn)地震數(shù)據(jù)去噪和插值的同步處理;

2) 為重建空缺地震道數(shù)據(jù),同時保持原始數(shù)據(jù)的振幅特征,在降秩過程中需利用POCS算法進(jìn)行迭代處理;

3) 對于數(shù)據(jù)量較大的實(shí)際地震資料,采用分窗方式進(jìn)行處理可以減少算法對內(nèi)存的需求和運(yùn)算量,同時使有效信號的同相軸線性化。

參考文獻(xiàn)

[1]SERGIO L M F,TAD J U.Application of singular value decomposition to vertical seismic profiling[J].Geophysics,1988,53(6):778-785

[2]高靜懷,朱光明,王玉貴.奇異值分解在VSP中的應(yīng)用[J].西安地質(zhì)學(xué)院學(xué)報,1992,14(3):71-78

GAO J H,ZHU G M,WANG Y G.Application of singular value decomposition to VSP[J].Journal Xi’an College of Geology,1992,14(3):71-78

[3]徐伯勛,白應(yīng)甫,白旭濱,等.奇異值分解在垂直地震剖面中的應(yīng)用[J].石油物探,1993,32(3):60-72

XU B X,BAI Y F,BAI X B,et al.Application of singular value decomposition (SVD) to VSP[J].Geophysical Prospecting for Petroleum,1993,32(3):60-72

[4]陳遵德,段天友,朱廣生.SVD濾波方法的改進(jìn)與應(yīng)用[J].石油地球物理勘探,1994,29(6):783-792

CHEN Z D,DUAN T Y,ZHU G S.Improvement of singular-value-decomposition filtering and the application[J].Oil Geophysical Prospecting,1994,29(6):783-792

[5]陸文凱,牟永光.一種改進(jìn)的SVD濾波器[J].石油地球物理勘探,1996,31(5):736-741

LU W K,MOU Y G.An improved SVD filter[J].Oil Geophysical Prospecting,1996,31(5):736-741

[6]吳亞東,符溪,文鵬飛,等.奇異值分解壓制隨機(jī)噪聲的方法及應(yīng)用[J].新疆石油地質(zhì),2003,24(2):144-145

WU Y D,FU X,WEN P F,et al.SVD method and its application in attenuating Random noise[J].Xinjiang Petroleum Geology,2003,24(2):144-145

[7]李文杰,魏修成,劉洋,等.SVD濾波法在直達(dá)波和折射波衰減處理中的應(yīng)用[J].石油勘探與開發(fā),2004,31(5):71-73

LI W J,WEI X C,LIU Y,et al.Application of SVD filtering in the processing of decaying direct wave and refracted wave[J].Petroleum Exploration and Development,2004,31(5):71-73

[8]SACCHI M.FX singular spectrum analysis[R].Calgary,Alberta,Canada:CSPG CSEG CWLS Convention,2009:392-395

[9]魏小強(qiáng),雷秀麗,馬慶珍.基于多道奇異譜分析的三維地震數(shù)據(jù)規(guī)則化方法[J].石油地球物理勘探,2014,49(5):846-851

WEI X Q,LEI X L,MA Q Z.3D seismic data regularization based on multi-channel singular spectrum analysis[J].Oil Geophysical Prospecting,2014,49(5):846-851

[10]黃建平,李闖,李國磊,等.基于奇異譜分析的聯(lián)合去噪及規(guī)則化方法[J].地球物理學(xué)進(jìn)展,2014,29(4):1666-1671

HUANG J P,LI C,LI G L,et al.Simultaneous seismic data de-noising and regularization method based on singular spectrum analysis[J].Progress in Geophysics (in Chinese),2014,29(4):1666-1671

[11]GOLYANDINA N,NEKRUTKIN V,ZHIGLJAVSKY A.Analysis of time series structure:SSA and related techniques[M].Boca Raton,Fla:Chapman & Hall/CRC,2001:28-29

[12]ABMA R,KABIR N.3D interpolation of irregular data with a POCS algorithm[J].Geophysics,2006,71(6):E91-E97

(編輯：戴春秋)

Seismic data noise attenuation and interpolation using singular value decomposition in frequency domain

MA Jitao1,WANG Jianhua2,LIU Guochang1

(1.GeophysicsDepartment,CollegeofGeophysicsandInformationEngineering,ChinaUniversityofPetroleum,Beijing102249,China; 2.CNOOCResearchInstitute,OffshoreOilExplorationNationalEngineeringLaboratory,Beijing100028,China)

Abstract:Singular value decomposition (SVD) is a seismic random noise attenuation method developed in recent years.Using the rank reduction character of frequency singular value decomposition,and by a projection onto convex sets (POCS) iterative algorithm,a method that can simultaneously process seismic data interpolation and random noise reduction is realized,and the anti-diagonal averaging process is characterized in detail.Sliding window in real data processing can linearize the seismic events in the localized window,reduce the memory storage requirement,and decrease the computation cost.The method is tested with synthetic data and real datasets,and the results illustrate that the rank reduction method using singular value decomposition in the frequency domain can suppress seismic noise while interpolate seismic data well,which shows a broad application prospect.

Keywords:singular value decomposition (SVD),frequency domain,interpolation,noise attenuation

文章編號：1000-1441(2016)02-0205-09

DOI:10.3969/j.issn.1000-1441.2016.02.006

中圖分類號：P631

文獻(xiàn)標(biāo)識碼：A

基金項目：國家自然科學(xué)基金(41404099)、中國石油大學(xué)(北京)科研基金項目(2462015YQ0512)、海洋石油勘探國家工程實(shí)驗室“斜纜采集地震數(shù)據(jù)分析與處理技術(shù)研究”項目聯(lián)合資助。

作者簡介：馬繼濤(1983—),男,博士,講師,主要從事地震數(shù)據(jù)處理方法研究工作。

收稿日期：2015-05-07;改回日期：2015-10-25。

This research is financially supported by the National Natural Science Foundation of China (Grant No.41404099),the Science Foundation of China University of Petroleum,Beijing (Grant No.2462015YQ0512) and the Project from National Engineering Laboratory for Offshore Oil Exploration.