基于層次聚類多道奇異譜分析的地震數(shù)據(jù)同時(shí)重建與去噪方法

曹靜杰，許昌昊，朱躍飛

（1.自然資源部京津冀城市群地下空間智能探測(cè)與裝備重點(diǎn)實(shí)驗(yàn)室，河北石家莊 050031； 2.河北省戰(zhàn)略性關(guān)鍵礦產(chǎn)資源重點(diǎn)實(shí)驗(yàn)室，河北石家莊 050031； 3.河北地質(zhì)大學(xué)數(shù)理學(xué)院，河北石家莊 050031； 4.中國(guó)礦業(yè)大學(xué)（北京）地球科學(xué)與測(cè)繪工程學(xué)院，北京 100083）

0 引言

地震勘探數(shù)據(jù)處理的最終目標(biāo)是獲得高信噪比、高保真度、高分辨率的反映地下結(jié)構(gòu)的數(shù)據(jù)，而提高地震數(shù)據(jù)信噪比是獲得高保真和高分辨率地震圖像的基礎(chǔ)。在地震數(shù)據(jù)采集過(guò)程中，由于環(huán)境及壞道等因素，經(jīng)常會(huì)出現(xiàn)地震道缺失及隨機(jī)噪聲干擾等問(wèn)題，降低了地震數(shù)據(jù)的信噪比。地震數(shù)據(jù)重建與去噪是獲得高信噪比完備地震數(shù)據(jù)的主要手段，對(duì)于地震數(shù)據(jù)處理的后續(xù)環(huán)節(jié)，如AVO 分析［1］、地震屬性分析、波動(dòng)方程偏移成像［2］等，都有非常重要的影響。

地震數(shù)據(jù)重建與去噪方法可分為四類。最早應(yīng)用的方法是基于波動(dòng)方程，通過(guò)速度模型正演地震道記錄，代替缺失的地震數(shù)據(jù)［3-4］。這類方法的重建精度依賴于速度的準(zhǔn)確性，且計(jì)算量大，難以在實(shí)際應(yīng)用中普及。

第二類是基于稀疏變換的方法。該類方法不依賴地下介質(zhì)速度等信息，以地震數(shù)據(jù)在某些變換域的稀疏性為先驗(yàn)約束，構(gòu)建反演模型并求解得到重建的數(shù)據(jù)。壓縮感知理論出現(xiàn)后，該類方法迅速應(yīng)用于地震勘探領(lǐng)域。常用的稀疏變換包括Fourier 變換［5］、Radon 變換［6-7］、Curvelet 變換［8-10］及Seislet 變換［11-12］等。Hennenfent等［8］提出一種基于非規(guī)則采樣Curvelet 變換和閾值法對(duì)地震數(shù)據(jù)進(jìn)行重建與去噪； 仝中飛［13］將Curvelet 變換與閾值迭代法相結(jié)合，提出Curvelet 變換閾值迭代法，把地震數(shù)據(jù)隨機(jī)噪聲衰減和地震缺失道重建問(wèn)題描述為L(zhǎng)1模最優(yōu)化問(wèn)題，利用閾值迭代法求解以達(dá)到去噪與重建的目的； Cao等［14］將低冗余曲波變換引入地震數(shù)據(jù)重建，大幅度提高了曲波變換在數(shù)據(jù)重建中的計(jì)算速度； 葛子建［15］將Fourier變換法與凸集投影法相結(jié)合，建立了一套有效的針對(duì)不規(guī)則缺失、含假頻地震數(shù)據(jù)和含噪數(shù)據(jù)的重建方法； 曹靜杰等［16］由迭代閾值法推導(dǎo)出加權(quán)凸集投影法，且證明它是解無(wú)約束優(yōu)化問(wèn)題的一種有效方法，加權(quán)因子可看作擬合誤差項(xiàng)系數(shù)，進(jìn)而提出一種改進(jìn)的凸集投影法，與原始凸集投影法相比不需增加任何計(jì)算量，只需通過(guò)選取閾值進(jìn)行重建與去噪。由于采用的變換不同，因此各種方法的特點(diǎn)也不同?；贔ourier變換的方法計(jì)算效率高，但基于線性同相軸假設(shè)，需將數(shù)據(jù)分塊； 基于Radon變換的方法需采用反演方法求解Radon變換系數(shù)，計(jì)算效率低； 基于Curvelet變換的方法，冗余度高，計(jì)算效率低，但重建質(zhì)量高。

第三類是基于矩陣降秩的方法［17］。矩陣降秩方法是根據(jù)相鄰地震道之間的相干性壓制隨機(jī)噪聲，減少去噪時(shí)有效信號(hào)的損失，且通過(guò)迭代重建方法補(bǔ)全由于空間采樣不均勻造成的地震道缺失［18-19］。在缺失的地震數(shù)據(jù)矩陣中，缺失的列會(huì)導(dǎo)致解不穩(wěn)定，一般需將缺失地震數(shù)據(jù)排成Hankel 矩陣，再通過(guò)對(duì)Hankel矩陣降秩實(shí)現(xiàn)數(shù)據(jù)重建與噪聲消除。該類方法基于奇異譜分解（SVD），最早應(yīng)用于核磁共振圖像處理，稱為Cadzow 濾波［20］，后被Trickett［21］引入地震數(shù)據(jù)處理領(lǐng)域，并發(fā)展成多道奇異譜分析（MSSA）方法［18，22］。黃建平等［23］提出一種基于SVD 的聯(lián)合去噪、規(guī)則化方法，迭代時(shí)在地震道缺失位置進(jìn)行重建，在不缺失地震道位置壓制噪聲。Huang 等［24］將阻尼因子引入MSSA 去噪方法，提出阻尼多道奇異譜分析（DMSSA）方法，以阻尼因子減弱隨機(jī)噪聲帶來(lái)的奇異值的增量，將信號(hào)子空間與噪聲子空間進(jìn)行更徹底的分離，但并未討論地震數(shù)據(jù)重建問(wèn)題。馬繼濤等［25］基于頻率域SVD 矩陣降秩運(yùn)算，利用凸集投影迭代方法，實(shí)現(xiàn)了地震數(shù)據(jù)去噪和重建的同步處理。Chen 等［26］將阻尼算子引入5D 地震數(shù)據(jù)去噪和重建，使矩陣降秩方法在地震數(shù)據(jù)重建中更好地壓制了隨機(jī)噪聲，但并未對(duì)奇異值數(shù)量的選取進(jìn)行相關(guān)研究。

矩陣降秩方法同樣基于線性同相軸假設(shè)，一般需將數(shù)據(jù)分塊，當(dāng)重排的Hankel矩陣維數(shù)較大時(shí)，SVD需耗費(fèi)較長(zhǎng)時(shí)間，且須首先選取每個(gè)數(shù)據(jù)塊的奇異值個(gè)數(shù)。朱躍飛等［27］探討了自動(dòng)確定奇異值個(gè)數(shù)的方法，利用一種Akaike信息準(zhǔn)則自動(dòng)地確定地震信號(hào)的奇異值個(gè)數(shù)，然后基于DMSSA 方法的框架去噪，避免了通過(guò)人工選取每個(gè)數(shù)據(jù)塊的奇異值個(gè)數(shù)。

第四類是基于人工智能的方法。由于近年來(lái)計(jì)算機(jī)硬件的不斷更新，人工智能得到快速發(fā)展，特別是深度學(xué)習(xí)理論在各領(lǐng)域都產(chǎn)生了極大影響。字典學(xué)習(xí)方法［28-29］是其中的一個(gè)分支，該方法與基于稀疏變換法的區(qū)別在于，基函數(shù)不固定和對(duì)復(fù)雜構(gòu)造有更好的表達(dá)能力，但其計(jì)算量大、泛化能力較差，難以實(shí)施大規(guī)模數(shù)據(jù)計(jì)算。深度學(xué)習(xí)方法，如卷積神經(jīng)網(wǎng)絡(luò)［30］、卷積自編碼器［31］、殘差學(xué)習(xí)［32］等，在地震數(shù)據(jù)去噪、重建、屬性分析、斷層解釋、初值拾取等領(lǐng)域，已得到一些應(yīng)用。針對(duì)地震數(shù)據(jù)重建，江金生［33］提出一種多層、多模塊的卷積自編碼器，改進(jìn)后的卷積自編碼器可更有效地提取地震數(shù)據(jù)特征，訓(xùn)練后的卷積自編碼器深度神經(jīng)網(wǎng)絡(luò)可有效地完成地震數(shù)據(jù)中隨機(jī)噪聲、線性噪聲、面波等的壓制及地震數(shù)據(jù)重建。Liu等［34］提出基于部分卷積方法提高重建結(jié)果中的模糊問(wèn)題，使同相軸更清晰。Siahkoohi等［35］通過(guò)炮—檢互易原理，利用共檢波點(diǎn)頻率—波數(shù)域隨機(jī)缺失數(shù)據(jù)訓(xùn)練卷積網(wǎng)絡(luò)，可實(shí)現(xiàn)很大比例缺失數(shù)據(jù)的重建。王峰［36］在地震數(shù)據(jù)去噪中采用DnCNN 模型，讓卷積神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)數(shù)據(jù)中噪聲的特征，建立統(tǒng)一用于二維和三維地震數(shù)據(jù)、同時(shí)去噪和重建的模型。除了卷積網(wǎng)絡(luò)，Wang 等［37］將ResNet 用于地震數(shù)據(jù)重建，生成對(duì)抗網(wǎng)絡(luò)GAN［38］和U-net 網(wǎng)絡(luò)［39-40］也被用于數(shù)據(jù)重建。

MSSA 法可通過(guò)矩陣降秩法有效壓制噪聲，并重建缺失地震道，其依據(jù)為隨機(jī)噪聲的存在會(huì)引起地震數(shù)據(jù)在頻率域矩陣秩的增加，通過(guò)降秩法壓制隨機(jī)噪聲。即首先將地震數(shù)據(jù)進(jìn)行Fourier 變換，再將每一頻率分量變?yōu)閴KHankel 矩陣，然后對(duì)該Hankel 矩陣進(jìn)行奇異值截?cái)?，通過(guò)多次迭代達(dá)到去噪及重建的效果。如何讓計(jì)算機(jī)自動(dòng)確定奇異值的個(gè)數(shù)，目前尚無(wú)很好的解決方案，但這對(duì)該類方法能否具有工業(yè)應(yīng)用價(jià)值至關(guān)重要。

針對(duì)上述問(wèn)題，本文基于DMSSA 的地震數(shù)據(jù)重建與去噪框架，提出一種層次聚類法自動(dòng)確定奇異值的個(gè)數(shù)，即通過(guò)聚類法將奇異值分為兩類，設(shè)定其中一類對(duì)應(yīng)有效信號(hào)，另一類對(duì)應(yīng)噪聲，且在有效信號(hào)頻率范圍內(nèi)取最大值以確定奇異值個(gè)數(shù)，從而實(shí)現(xiàn)地震數(shù)據(jù)同時(shí)重建和去噪。數(shù)值實(shí)驗(yàn)表明，該方法對(duì)模擬數(shù)據(jù)和實(shí)際數(shù)據(jù)均能取得較好效果。

1 方法原理

1.1 MSSA 去噪原理

設(shè)大小為M×N×O的三維含噪地震數(shù)據(jù)為D(x，y，t)，且x=x1，x2，…，xM；y=y1，y2，…，yN；t=t1，t2，…，tO。通過(guò)以下五步實(shí)現(xiàn)MSSA 法去噪。

第一步，對(duì)D(x，y，t)做離散Fourier 變換，使地震數(shù)據(jù)從時(shí)間域變換到頻率域，即F(x，y，ω)，其中ω=ω1，ω2，…，ωJ，J為信號(hào)的頻率個(gè)數(shù)。

第二步，對(duì)F(x，y，ω)每一頻率分量做Hankel變換。顯然

是F(x，y，ω)中頻率為ωi的分量。將F(x，y，ωi)的每行構(gòu)成Hankel矩陣

式中：Aq為F(x，y，ωi)第q行構(gòu)成的Hankel 矩陣；；Ky=N-Kx+1； 符號(hào)表示下取整運(yùn)算。將所有Aq排列成一個(gè)塊Hankel矩陣

第三步，對(duì)Hωi進(jìn)行SVD，得到

式中：Uωi為L(zhǎng)x×Kx階酉矩陣；為L(zhǎng)y×Ky階酉矩陣；Σ=diag(σ1，σ2，…，σs)為對(duì)角矩陣（s為非零奇異值個(gè)數(shù)），且σ1＞σ2＞…＞σs。

第四步，對(duì)Hωi做降秩處理，確定Σ中合理的奇異值個(gè)數(shù)。假設(shè)地震信號(hào)中有P個(gè)不同斜率的線性同相軸，對(duì)Σ中的奇異值做截?cái)嗵幚?，保留前P個(gè)奇異值，其余奇異值設(shè)為0，可達(dá)到降秩的效果。其過(guò)程可表示為

式中：Γ-1表示Hankel 矩陣變換的逆；Ψ-1為Fourier反變換算子；FMSSA(x，y，ω)為(x，y，ωi)的集合；(x，y，t)為輸出的去噪結(jié)果。

1.2 DMSSA 去噪原理

MSSA 只能將數(shù)據(jù)分解為噪聲子空間和信號(hào)加噪聲子空間，并不能完全去除隨機(jī)噪聲。當(dāng)?shù)卣饠?shù)據(jù)中只有有效信號(hào)時(shí)，Σ'中包含與地震圖像中線性同相軸數(shù)量相同個(gè)數(shù)的非零奇異值。當(dāng)?shù)卣饠?shù)據(jù)混入隨機(jī)噪聲，所有的奇異值大小都會(huì)發(fā)生改變，非零奇異值個(gè)數(shù)將會(huì)增加。傳統(tǒng)的MSSA 方法去噪保留前P個(gè)奇異值，其余奇異值置0，但保留的前P個(gè)奇異值中依然殘留隨機(jī)噪聲引起的變化量，其去噪結(jié)果中依然殘存隨機(jī)噪聲。

Huang 等［24］提出的DMSSA 方法，將阻尼因子引入MSSA 去噪方法，通過(guò)阻尼因子減弱隨機(jī)噪聲造成的奇異值增量，有效壓制殘存的隨機(jī)噪聲。

式（5）和式（6）被替換為下列兩式

式中：I是單位矩陣；D表示阻尼因子；T表示阻尼算子；P為預(yù)估的奇異值個(gè)數(shù)。通過(guò)式（9）和式（10）兩式可利用阻尼因子D對(duì)第P+1 個(gè)奇異值進(jìn)行放大或縮小，用前P個(gè)奇異值與其相減，達(dá)到壓制前P個(gè)奇異值中與隨機(jī)噪聲相關(guān)部分的目的。

值得一提的是，D值越小，阻尼效果越強(qiáng)，其去噪效果越好，但會(huì)對(duì)有效信號(hào)造成輕微損害。當(dāng)D→∞時(shí)，式（10）退化為式（6）。根據(jù)地震數(shù)據(jù)的復(fù)雜度及噪聲的強(qiáng)弱，D的選取范圍通常是1～5。本文研究重點(diǎn)在估計(jì)有效信號(hào)對(duì)應(yīng)的奇異值，未討論阻尼因子的選擇，但經(jīng)過(guò)實(shí)驗(yàn)選擇阻尼因子為2，該阻尼因子能得到可靠的結(jié)果。將上述過(guò)程用DMSSA 表示。

1.3 DMSSA 同時(shí)重建和去噪原理

缺失的地震道通常在地震數(shù)據(jù)中以0 的形式存在，在對(duì)地震數(shù)據(jù)進(jìn)行重建時(shí)，地震道缺失部分可看作有效信號(hào)與隨機(jī)噪聲相加為0 的情況，因此三維地震數(shù)據(jù)重建法的步驟與去噪原理的步驟相似，只是在重建空白地震道過(guò)程中采用了類似加權(quán)凸集投影的方式。具體重建流程如表1。

表1 DMSSA 同時(shí)重建和去噪流程

表1 中：F0(ω)=F(ω)，an是一個(gè)隨循環(huán)線性減少的標(biāo)量，a1=1，anmax=0，為采樣算子。

上述方法的含義是，對(duì)于頻率ω1～ωmax范圍內(nèi)的每個(gè)頻率ω，設(shè)定最大迭代次數(shù)為nmax，每次迭代都執(zhí)行表1 中Fn(ω)的計(jì)算過(guò)程，如果迭代次數(shù)未達(dá)到nmax，但是滿足，則返回Fn(ω)為數(shù)據(jù)結(jié)果，將這些頻率范圍都執(zhí)行上述運(yùn)算，就得到頻率域的同時(shí)重建和去噪結(jié)果，再通過(guò)Fourier反變換就能得到時(shí)間域的重建與去噪結(jié)果。

該方法通過(guò)循環(huán)迭代的方式完成地震數(shù)據(jù)的重建與去噪，每次迭代中，在地震道缺失位置處由ADMSSA 方法進(jìn)行重建，在地震道沒(méi)有缺失的位置由anF(ω)+(1-an)·DMSSA[Fn-1(ω) ]對(duì)數(shù)據(jù)進(jìn)行去噪，F(xiàn)(ω)為原始觀測(cè)數(shù)據(jù)，F(xiàn)n-1(ω)為前一輪奇異譜分析后得到的數(shù)據(jù)，an控制兩者所占的比率，由1 線性遞減至0。該方法可有效降低直接重建時(shí)由于原始觀測(cè)數(shù)據(jù)的噪聲對(duì)數(shù)據(jù)相干性造成的影響，實(shí)現(xiàn)地震數(shù)據(jù)的同時(shí)重建和去噪。

1.4 層次聚類法自動(dòng)確定奇異值個(gè)數(shù)

DMSSA 類方法實(shí)現(xiàn)同時(shí)重建和去噪時(shí)，奇異值個(gè)數(shù)P的選取至關(guān)重要。若奇異值個(gè)數(shù)選取過(guò)大，則重建后殘留的噪聲過(guò)大。若奇異值個(gè)數(shù)選取過(guò)少，則會(huì)損害有效信號(hào)。在計(jì)算機(jī)運(yùn)算過(guò)程中，由于實(shí)際地震資料的數(shù)據(jù)量過(guò)大，常需對(duì)地震數(shù)據(jù)進(jìn)行分塊處理； 而每塊的同相軸數(shù)量并不相同，因此需要人工確定P，無(wú)法實(shí)現(xiàn)自動(dòng)化重建與去噪。本節(jié)通過(guò)對(duì)奇異值的分布特點(diǎn)進(jìn)行分析，采用層次聚類方法使計(jì)算機(jī)自動(dòng)確定合適的奇異值個(gè)數(shù)P，進(jìn)而更好地完成地震數(shù)據(jù)的重建。

圖1a 是采用主頻為30 Hz 的雷克子波正演得到一塊大小為40×40×300的模擬數(shù)據(jù)，時(shí)間采樣率為2 ms，該數(shù)據(jù)有3 個(gè)線性同相軸。圖1b 為隨機(jī)抽取50%地震道且增加10%隨機(jī)噪聲后的數(shù)據(jù)，信噪比為-3.9 dB。圖2為頻率為40 Hz時(shí)塊Hankel矩陣的奇異值分布，其中綠線為圖1a 數(shù)據(jù)對(duì)應(yīng)的奇異值，紅線為待重建地震數(shù)據(jù)（圖1b）對(duì)應(yīng)的奇異值。由圖2中看出，圖1a地震數(shù)據(jù)第4 項(xiàng)及之后的奇異值遠(yuǎn)小于前3 項(xiàng)奇異值，且絕大多數(shù)為0。隨機(jī)噪聲的增加與地震道的缺失使奇異值序列中前3 項(xiàng)奇異值變小，非零項(xiàng)增多，但前3 項(xiàng)奇異值依然與3 項(xiàng)之后的奇異值有較大區(qū)別。上述含噪聲地震數(shù)據(jù)的奇異值可分為兩類：一類是與有效信號(hào)有關(guān)的奇異值，數(shù)值較大； 另一類是與噪聲和道缺失有關(guān)的奇異值，數(shù)值較小。本文引入一種層次聚類的方法，將多個(gè)元素按照相似程度自下而上的進(jìn)行分類，從而確定奇異值個(gè)數(shù)。

圖1 模擬數(shù)據(jù)一

圖2 圖1a 和圖1b 頻率為40 Hz 時(shí)的奇異值分布

定義類與類之間的距離

式中：Gp、Gq為兩個(gè)不同的類；zi和zj均為二維向量，分別屬于Gp和Gq；dij為zi與zj之間的歐式距離。

將式（4）中Σ的奇異值排列成為二維數(shù)組

使用(z1，z2，…，zs)進(jìn)行表示。

層次聚類的具體過(guò)程如下：

（1）將(z1，z2，…，zs)中每個(gè)數(shù)組各自分為一類，計(jì)算各類之間的距離；

（2）合并類間距最小的兩個(gè)類，構(gòu)成一個(gè)新類；

（3）計(jì)算新類與當(dāng)前各類的距離；

（4）若類的個(gè)數(shù)為2，即終止計(jì)算； 否則，返回步驟（2）。

上述過(guò)程將奇異值分為兩類，由于奇異值個(gè)數(shù)P遠(yuǎn)遠(yuǎn)小于奇異值的總數(shù)量s，故選取數(shù)量較少的一類，計(jì)算其包含奇異值的數(shù)量即可得到P。

研究發(fā)現(xiàn)，使用DMSSA 方法時(shí)只需對(duì)有效信號(hào)頻率范圍進(jìn)行聚類分析，對(duì)有效頻率范圍以外可取值為零。在有效頻率范圍內(nèi)，每個(gè)頻率分量的聚類結(jié)果也有可能發(fā)生變化，為了使該方法結(jié)果穩(wěn)定，選擇有效信號(hào)頻率范圍內(nèi)聚類，以結(jié)果的最大值作為最終的奇異值個(gè)數(shù)。

為驗(yàn)證上述觀點(diǎn)，圖3給出了圖1數(shù)據(jù)1～150 Hz層次聚類法測(cè)得的奇異值數(shù)量。對(duì)圖1 數(shù)據(jù)而言，10～90 Hz 在信號(hào)的有效頻率范圍內(nèi)，層次聚類法確定的奇異值個(gè)數(shù)穩(wěn)定在3； 超出該范圍時(shí)測(cè)得的奇異值數(shù)量不斷上下波動(dòng)，其原因是超出該頻率范圍的為噪聲。根據(jù)MSSA 理論，可驗(yàn)證在10～90 Hz 內(nèi)測(cè)得的奇異值數(shù)量是準(zhǔn)確的。

圖3 模擬數(shù)據(jù)一層次聚類的結(jié)果

2 模擬數(shù)據(jù)實(shí)驗(yàn)

首先用模擬數(shù)據(jù)驗(yàn)證方法的有效性。圖4a 為采用層次聚類法確定的奇異值個(gè)數(shù)P=3時(shí)圖1b的重建結(jié)果，阻尼因子D取2，重建后信噪比為19.9 dB；圖4b為重建誤差。從圖4可看出，重建后線性同相軸清晰可見(jiàn)，噪聲殘留較少，說(shuō)明該方法預(yù)測(cè)的P值準(zhǔn)確，能重建出較高信噪比的地震圖像。

圖4 模擬數(shù)據(jù)一重建效果

為了驗(yàn)證方法的穩(wěn)定性，下面采用另一組模擬數(shù)據(jù)進(jìn)行測(cè)試。圖5a 是對(duì)主頻為30 Hz 的雷克子波合成得到的一塊維數(shù)為40×40 ×300 的模擬地震數(shù)據(jù)，時(shí)間采樣率為2 ms，該數(shù)據(jù)有5 個(gè)線性同相軸。圖5b 為隨機(jī)抽取50%地震道且增加10%隨機(jī)噪聲后的數(shù)據(jù)，其信噪比為-0.7 dB。

圖5 模擬數(shù)據(jù)二

圖6 為通過(guò)層次聚類測(cè)得的奇異值個(gè)數(shù)，可見(jiàn)在10～90 Hz頻率范圍內(nèi)，測(cè)得的P值為5。

圖6 模擬數(shù)據(jù)二層次聚類的結(jié)果

圖7a為模擬數(shù)據(jù)二基于層次聚類法的重建結(jié)果，阻尼因子取2，重建后信噪比為18.6 dB，圖7b 為重建誤差。圖4與圖7兩組模擬地震數(shù)據(jù)的重建效果都相對(duì)較好，線性同相軸恢復(fù)清晰，噪聲殘留較少，證明了本文方法預(yù)測(cè)奇異值數(shù)量的可行性與準(zhǔn)確性。

圖7 模擬數(shù)據(jù)二重建效果

3 實(shí)際地震數(shù)據(jù)實(shí)驗(yàn)

為了驗(yàn)證本文方法對(duì)實(shí)際數(shù)據(jù)的處理效果，選用了兩塊三維疊后地震數(shù)據(jù)。

圖8a是一塊大小為40×40×276的三維實(shí)際地震數(shù)據(jù)，時(shí)間上有276 個(gè)采樣點(diǎn)，時(shí)間采樣率為2 ms。主測(cè)線方向與聯(lián)絡(luò)測(cè)線方向上各有40 個(gè)采樣點(diǎn)。對(duì)其加入20%的隨機(jī)噪聲（圖8b），然后進(jìn)行隨機(jī)30%地震道抽?。▓D8c）?？梢?jiàn)該數(shù)據(jù)噪聲能量較強(qiáng)，信噪比低，地震道缺失嚴(yán)重。

圖8 三維實(shí)際地震數(shù)據(jù)一及其加噪和缺失數(shù)據(jù)

在計(jì)算大規(guī)模實(shí)際地震數(shù)據(jù)去噪和重建時(shí)需做分塊處理。通過(guò)滑動(dòng)分塊將大塊地震數(shù)據(jù)分成各個(gè)小塊，對(duì)每個(gè)小塊通過(guò)DMSSA 方法進(jìn)行重建。通過(guò)多次比較，發(fā)現(xiàn)空間分塊采用10×10 左右時(shí)能取得穩(wěn)定的效果。

圖9 為對(duì)其中一個(gè)小塊通過(guò)層次聚類法確定奇異值的結(jié)果?？梢?jiàn)實(shí)際數(shù)據(jù)相比于模擬數(shù)據(jù)，有較大的差異性，奇異值波動(dòng)更大；超出有效頻率（10～90 Hz）范圍后，測(cè)得的奇異值個(gè)數(shù)出現(xiàn)嚴(yán)重波動(dòng)；在有效頻率范圍內(nèi)，奇異值測(cè)得的數(shù)量相對(duì)較為穩(wěn)定，但與模擬地震數(shù)據(jù)不同的是，在有效頻率范圍內(nèi)測(cè)得的奇異值數(shù)量并不完全一樣。如何選擇恰當(dāng)?shù)钠娈愔敌柽M(jìn)行細(xì)致探討。若選取的奇異值過(guò)大，則會(huì)影響去噪效果； 若選取的奇異值過(guò)小，則會(huì)損害有效信號(hào)。因此，本文對(duì)各頻率測(cè)得的奇異值個(gè)數(shù)進(jìn)行統(tǒng)計(jì)（表2），并按測(cè)得數(shù)量從大到小排列，取其前面90%，選取其中最大值作為該小塊最終確認(rèn)的奇異值，這樣選取奇異值可避免損傷大部分有效信號(hào)，且不會(huì)對(duì)去噪效果產(chǎn)生較大影響。

圖9 實(shí)際數(shù)據(jù)一某一分塊基于層次聚類法確定的奇異值個(gè)數(shù)

表2 圖9 中測(cè)得奇異值個(gè)數(shù)統(tǒng)計(jì)

對(duì)圖9 在有效頻率范圍內(nèi)（10～90 Hz）的各頻率所測(cè)奇異值的個(gè)數(shù)進(jìn)行統(tǒng)計(jì)（表2）。可見(jiàn)測(cè)得的奇異值個(gè)數(shù)為1、2、3 的頻率的累計(jì)總數(shù)為76，占所有頻率個(gè)數(shù)的93.76%，因此取奇異值個(gè)數(shù)1、2、3 中的最大值3 可避免對(duì)有效信號(hào)造成損傷，且會(huì)得到較好的去噪效果，顯然此小塊在有效頻率范圍內(nèi)測(cè)得的奇異值的個(gè)數(shù)為3。

對(duì)圖8 實(shí)際地震數(shù)據(jù)采用滑動(dòng)分塊重建方法，每小塊大小為10×10，在主測(cè)線方向與聯(lián)絡(luò)測(cè)線方向均取步長(zhǎng)為2，共分得256塊。圖10a為基于層次聚類法確定奇異值的重建結(jié)果，阻尼因子取2，其每塊測(cè)得的奇異值數(shù)量如圖11 所示。圖10b 為去噪后殘留的噪聲。可見(jiàn)該實(shí)際數(shù)據(jù)的重建與去噪效果（圖10a）良好，有效信號(hào)損失較少，剖面的同相軸輪廓恢復(fù)清晰，噪聲清除干凈，證明基于本文方法在實(shí)際地震數(shù)據(jù)中的有效性與可行性。

圖10 實(shí)際數(shù)據(jù)一重建效果

圖11 對(duì)所有分塊測(cè)得的奇異值

對(duì)圖8實(shí)際地震重建中，由圖11可見(jiàn)大部分?jǐn)?shù)據(jù)塊中測(cè)得的奇異值在1～3。為了驗(yàn)證層次聚類法所測(cè)奇異值的準(zhǔn)確性，指定所有數(shù)據(jù)塊的奇異值全部為1、2、3，所得到的重建結(jié)果與本文每個(gè)數(shù)據(jù)塊都自動(dòng)確定奇異值個(gè)數(shù)的方法所得重建結(jié)果進(jìn)行對(duì)比。

從圖12 可見(jiàn)，當(dāng)對(duì)每個(gè)數(shù)據(jù)塊的奇異值直接取1時(shí)，其圖像左下方的斷層變得不再明顯。對(duì)每個(gè)數(shù)據(jù)塊的奇異值取2 和3 時(shí)有較好重建效果，但仍然殘留少量噪聲，且存在線性同相軸恢復(fù)不完整的情形。采用層次聚類確定奇異值的方法每個(gè)數(shù)據(jù)塊自動(dòng)確定適合的奇異值個(gè)數(shù)，重建圖像的噪聲去除干凈，線性同相軸清晰可見(jiàn)，圖像重建效果良好。

圖12 指定奇異值與通過(guò)層次聚類法取確定奇異值重建結(jié)果對(duì)比

圖13a 為三維實(shí)際地震數(shù)據(jù)實(shí)例二，大小為150×150× 276，時(shí)間上有276 個(gè)采樣點(diǎn)，時(shí)間采樣率為2 ms，主測(cè)線方向與聯(lián)絡(luò)測(cè)線方向各有150 個(gè)采樣點(diǎn)。對(duì)其加入20%隨機(jī)噪聲（圖13b），再對(duì)其隨機(jī)抽稀30%地震道（圖13c）。該數(shù)據(jù)同樣存在噪聲能量強(qiáng)、信噪比低、地震道缺失嚴(yán)重等特點(diǎn)。

圖13 三維實(shí)際地震數(shù)據(jù)二及其缺失的數(shù)據(jù)

圖14為針對(duì)圖13基于層次聚類法確定奇異值的重建結(jié)果，此時(shí)阻尼因子取值為2，得到去噪后的殘留噪聲（圖14b）。可見(jiàn)該重建與去噪結(jié)果中，同相軸輪廓恢復(fù)清晰，噪聲清除干凈，有效信號(hào)損失較少。再次證明了本文方法對(duì)實(shí)際地震數(shù)據(jù)的有效性與可行性。

圖14 三維實(shí)際地震數(shù)據(jù)二重建效果

4 討論

本文方法需要對(duì)地震信號(hào)有效頻率進(jìn)行分析，選定有效信號(hào)的范圍。對(duì)地震數(shù)據(jù)進(jìn)行分塊處理和拼接處理也是必須的，分塊的大小對(duì)本文方法有一定影響。經(jīng)驗(yàn)表明：當(dāng)空間分塊大約為10×10時(shí)，能取得穩(wěn)定的結(jié)果； 拼接方法的選擇對(duì)于重建與去噪結(jié)果同樣會(huì)產(chǎn)生影響； 另外，阻尼因子也應(yīng)根據(jù)噪聲的強(qiáng)弱程度而選取。

5 結(jié)論

矩陣降秩方法是地震數(shù)據(jù)去噪與重建的一種重要方法，秩的選擇對(duì)于矩陣降秩效果的影響巨大。當(dāng)秩選取過(guò)小時(shí)，將損害有效信息；但當(dāng)秩選取過(guò)大時(shí)，會(huì)在去噪及重建后留下明顯的殘留噪聲。面對(duì)海量地震數(shù)據(jù)，人工選取奇異值的方法具有一定的主觀性，且不利于去噪和重建的工業(yè)化實(shí)現(xiàn)。本文通過(guò)對(duì)缺道且含有噪聲的地震數(shù)據(jù)的奇異值進(jìn)行分析，引入了層次聚類方法，將奇異值自動(dòng)分為有效信號(hào)和噪聲兩類，選擇在有效頻帶范圍內(nèi)實(shí)現(xiàn)聚類，使得多道奇異譜分析方法能自適應(yīng)且穩(wěn)定地確定奇異值的個(gè)數(shù)，并通過(guò)模擬數(shù)據(jù)與實(shí)際數(shù)據(jù)驗(yàn)證了方法的優(yōu)越性。實(shí)驗(yàn)表明，利用層次聚類法確定的奇異值數(shù)量可靠穩(wěn)定，重建后線性同相軸清晰可見(jiàn)，無(wú)太多殘留噪聲，地震數(shù)據(jù)重建效果良好。