基于PCA與K?SVD的地震數據去噪方法

關鍵詞：稀疏表示，主成分分析，降維處理，K‐SVD，去噪

0 引言

地震數據對地質勘探和地震監測至關重要，提供了了解地球內部結構和地質活動的寶貴信息[1]。地震數據的質量和可靠性常受多種噪聲影響，因此，地震數據去噪技術在地球物理研究領域備受關注[2]。目前，地震數據去噪方法可以分為以下幾類：基于濾波的方法，通過頻域分析或時域分析[3]過濾特定頻率或特定時間段內的噪聲信號，如帶通濾波[4]、陷波濾波[5]等；基于時頻分析的方法，利用小波變換[6]、短時傅里葉變換[7]等時頻分析技術，捕捉地震信號的時頻特征，精確定位和去除噪聲；基于降秩的方法，利用奇異值分解（SVD）[8]、主成分分析（PCA）[9]等降秩技術，降低數據的維度和復雜度，從而去除數據中的噪聲；基于學習的方法，如深度學習和字典學習等。

基于深度學習的地震數據去噪方法主要是利用卷積神經網絡（CNN）[10]、循環神經網絡（RNN）[11]等深度學習模型，通過學習地震數據的時空特征和時序關系，捕獲地震波形中復雜的空間和時間特征，在復雜噪聲環境下實現精準去噪，提高地震數據的信噪比。與之不同的是，基于字典學習的地震數據去噪方法則是將地震數據表示為稀疏線性組合，利用學習到的字典和稀疏系數去除噪聲[12]。蘇遠超等[13]在遙感圖像的研究中對字典學習進行了詳細論述。

嚴春滿等[14]提出的K‐SVD算法是一種基于迭代優化的字典學習算法，用于學習數據的稀疏表示模型，其基本思想是通過迭代更新字典中的原子（基向量）和稀疏表示系數，使稀疏表示能夠最好地逼近原始信號。通過這種方式，K‐SVD能夠將信號表示為一組稀疏的基向量線性組合，進而實現高效的信號表示和去噪處理。在后續的研究中，封常青等[15]提出了并行K‐SVD技術，將字典學習過程中的字典更新操作并行化，加快了算法的運行速度，此技術尤其適用于多核處理器和分布式計算環境；Mairal等[16]提出了增量K‐SVD算法，允許在有新數據時對字典進行增量學習，而不需要重新訓練整個字典，從而減少了計算開銷，更適用于動態數據場景。

傳統的K‐SVD算法在地震數據去噪中表現出色，但其計算復雜度較高，在每次迭代過程中需要進行SVD。SVD是一種數值穩定且可靠的矩陣分解方法，計算成本較高，尤其是在處理大規模數據時，會導致算法的時間復雜度呈指數增加[17]。傳統的K‐SVD算法在每次迭代中，需要對每個訓練樣本進行稀疏編碼，并更新字典中的原子。這個過程需要對數據矩陣進行奇異值分解，得到稀疏編碼的系數矩陣和更新后的字典。隨著數據規模增大，奇異值分解的計算成本也隨之增長，導致算法的效率急劇下降。

為了克服以上問題，本文提出了一種基于PCA與K‐SVD結合的地震數據去噪方法。首先，利用PCA技術對地震數據進行降維處理，將數據從高維空間映射到低維空間，降低數據的維度和復雜度[18]。由于PCA的計算成本遠低于SVD，在降維后的數據上進行K‐SVD的稀疏表示和去噪處理極大提高了算法的計算效率。經PCA降維后，地震數據的維度顯著減小，后續的K‐SVD算法在低維空間中進行稀疏表示時的計算量大幅減少，使算法更加適用于處理大規模地震數據。本文主要研究字典學習方法在地震數據去噪領域的應用，創新性地將PCA與K‐SVD相結合，通過降低數據維度和復雜度，有效解決了傳統K‐SVD算法中SVD計算復雜度高的問題，并探討了其在不同噪聲環境下的性能表現以及與傳統去噪方法的對比分析。

1 K?SVD理論

1.1 字典學習

字典學習是機器學習[19]的一種，旨在從數據中學習一組“字典”，以便將數據表示為這組字典中的基向量線性組合。假設有一個由N個樣本組成的數據集Y=[y1，y2，…，yN]，其中每個樣本yi是一個d維向量，表示數據的一個特征。字典學習的目標是找到一個字典矩陣D，使得數據集Y能夠最優地用字典中的基向量表示，字典學習可以表示為以下優化問題

式中：X是稀疏編碼矩陣；λ是稀疏正則化參數；xi是X中的第i列，它是數據集Y中的第i個樣本在字典矩陣D上的稀疏表示。

優化問題的第一項是重構誤差，表示原始數據與稀疏編碼重構的近似差異；第二項是稀疏性約束，用于促使編碼系數盡可能稀疏，因為L0難以求解，所以使用L1正則項替代近似。通過求解上述優化問題，可以得到最優的D和X，實現數據的稀疏表示和重構。

1.1.1 SVD

奇異值分解建立在特征值分解的基礎上，而特征值分解能夠從一個矩陣中提取其特征信息[20]。定義一個矩陣A

式中：η是矩陣A的特征值；α是A在特征值為η時所對應的特征向量。當矩陣能夠相似對角化時，A可被分解為

式中：M表示由A的特征向量構成的矩陣；Λ是一個對角矩陣，其對角線上的元素是A的特征值。

一個矩陣實質上代表了一個線性變換，在將矩陣作用于向量時，實際上是對該向量進行線性變換。特征值分解將矩陣分解為一個對角矩陣，其對角線上的元素按照特征值大小排列，這些特征值所對應的特征向量描述了矩陣變換的主要和次要方向。通過特征值分解得到的前N個特征向量對應矩陣最主要的N個變換方向。利用這些主要變換方向可以提取矩陣的特征。特征向量代表特征，特征值代表特征的重要程度，可在應用中作為權重。雖然這種分解方法十分有效，但它僅適用于方陣。為了處理一般的矩陣，需要引入更適用的奇異值分解[21]。

式中：U是一個m×m階的正交矩陣，其列向量被稱為左奇異向量；V是一個n×n階的正交矩陣，其列向量被稱為右奇異向量；Σ是一個對角矩陣，對角線上的元素為奇異值，并且按照從大到小的順序排列。

1.1.2 稀疏表示

在字典學習中，稀疏表示[22]是指利用盡可能少的非零系數表示一個信號或數據，其基本思想是假設數據在某個特定的基下是稀疏的，只有少數系數是非零的。在給定字典矩陣D和待表示的數據樣本Y的條件下，稀疏表示問題則是通過求解式（1），得到最優的稀疏編碼系數，從而實現數據的稀疏表示。

1.1.3 字典更新

字典學習的核心在于將原始數據Y分解為字典矩陣D和稀疏編碼矩陣X[23]。D存儲了Y中的特征，而X則表示如何利用這些特征進行組合。為了得出D和X，提出一個優化問題

這個優化問題的目標是最小化原始矩陣與字典矩陣乘以稀疏編碼矩陣后差異的F范數，同時還要確保X中的非零元素盡可能少。

通過迭代更新D的每一列和X的每一行，最小化原始數據與重構數據之間的差異。將D和X表示為列向量乘以行向量，式（5）可重新表達為

式中：s是字典原子的數量，它決定了字典D的規模；j是一個索引變量，用于遍歷字典原子。

把D中第i列的元素di和X中第i行的元素xi抽取出來，優化目標公式被重新表達為

式中：Ei是一個誤差矩陣；di和xi具有相同的維度，在該式里以列向量形式m×1階的矩陣參與運算，共同影響著最終的重構誤差；dixi的結果也是一個m×1階的矩陣，表示字典矩陣乘以稀疏碼矩陣，代表所恢復的結果，Y-DX為原始矩陣與恢復矩陣之間的誤差。而減去di和xi，則是在不考慮字典矩陣的第i列和稀疏碼矩陣的第i行的情況下的誤差。

當更新完di和xi進入下一次迭代時，Ei也會隨之改變，且無法保證二次更新時X的第i行是稀疏的。因此，把之前的xi中提取出的非零元素形成新的行向量xi′，從Ei中提取與xi′對應的列，形成新的誤差矩陣Ei'，優化目標公式為

該步驟的意義在于更新X的第i行時，繼承之前該行稀疏碼的稀疏性。通過提取非零元素與對應的非零列，保證X向稀疏的方向發展。無論最終得到的di和xi′結果如何，與Ei的誤差仍然存在，只能通過增加更多的迭代次數來彌補這一損失。求解式（9）是一個最小二乘問題，通過找到dixi′來逼近Ei'。

2 PCA理論

PCA是一種常用的數據分析技術，用于降低數據維度并發現數據內在的模式和結構[24]。其核心思想是通過線性變換將高維數據投影到低維空間，同時盡可能保留數據信息。在降維過程中，PCA通過尋找數據中的主要方差方向（主成分）實現降維。

PCA對數據降維重構的步驟為：

（1）去平均值；

（2）計算協方差矩陣；

（3）計算協方差矩陣的特征值與特征向量；

（4）對特征值從大到小排序，選擇其中最大的K個特征值，然后將其對應的K個特征向量分別作為行向量組成新的特征向量矩陣；

（5）利用新的特征向量矩陣重構原始圖像。

3 基于PCA與K?SVD的地震數據去噪方法

在之前的研究中，傳統的K‐SVD算法在地震數據去噪領域中雖然取得了顯著成就，但隨著地震數據量的大規模增加，傳統的K‐SVD算法在矩陣分解時需要耗費大量時間和計算資源，限制了方法的實用性。于是，本文結合PCA的數據降維和KSVD的稀疏表示學習，在保證數據質量的前提下，顯著降低了計算復雜度，提高了計算效率。

圖1為本文方法的實現流程，涵蓋了從數據預處理到PCA、K‐SVD以及最終地震數據去噪的完整過程。首先，使用PCA提取地震數據特征，通過線性變換將原始數據投影到一個新的坐標系中，使數據在新坐標系中的方差最大化，從而實現數據的降維；其次，使用經過降維重構的地震數據作為輸入，初始化字典，其中的原子代表可能出現在地震信號中的各種模式或特征；再次，利用稀疏編碼算法（如正交匹配追蹤，OMP）計算每個地震信號的稀疏系數，這些系數表示每個原子在地震信號中的重要程度，通過SVD等方法，利用稀疏系數更新字典中的原子，以最小化地震信號與稀疏編碼的重建誤差，同時保持字典的稀疏性；然后，在迭代過程中檢查原子是否更新完畢或者是否達到最大迭代次數，以此為條件來結束算法；最后，得到更新后的字典和每個地震信號對應的稀疏編碼，并利用稀疏編碼去除地震信號中的噪聲。

4 數值試算

在地震數據處理領域，去噪方法的通用性與魯棒性至關重要。文本主要測試基于PCA與K‐SVD的地震數據去噪方法在多樣化數據類型及復雜地質結構環境下的表現，全面評估其性能與優勢。選取具有代表性的疊前與疊后模擬地震數據并加入傳統去噪技術進行比較。

為了量化評估方法的去噪效果，定義信噪比（SNR）為衡量信號中噪聲水平相對量的重要指標，數學表達式為

式中：Ps是不含噪聲的原始地震數據的功率；Pn是噪聲數據的功率。

圖2為一個模擬的疊后地震記錄，該記錄共有510個采樣點，20道地震數據，采樣頻率為1ms，在加入標準差為0.3的隨機高斯噪聲后使用本文方法去噪。圖3顯示的是本文方法的去噪效果，可以看到地震同相軸較噪聲數據更清晰。在采樣時間50～250ms間同相軸交錯處仍可看到比較明顯的噪聲殘留，對其他結構的去噪效果較好。

圖4為加噪地震數據去噪后的殘差，沒有有效信號殘留，這表明本文方法能夠有效去除噪聲，保留有效信號。為了體現方法的優越性，在接下來的試驗中，將本文方法與傳統f?x反褶積方法對比，通過數值模擬單炮地震記錄并添加噪聲，以展示不同方法的去噪性能。

圖5a為數值模擬生成的單炮地震記錄，模擬四層均勻水平層狀介質。該模型共有500個采樣點，60道地震數據，采樣間隔為1ms。在去噪前對地震數據歸一化處理，其目的是將數據映射到統一的范圍內，消除不同尺度數據對算法的影響。通過歸一化處理，可以更好地處理具有不同尺度特征的數據。

圖5b為在數據歸一化后，加入標準差為0.3高斯白噪聲后的效果。由圖可見，在加入噪聲后地震數據的同相軸變得模糊，直接影響了地震數據中地下結構的辨識，給數據的分析和解釋帶來了困難。圖5c和圖5d是f?x反褶積和本文方法去噪結果的對比，f?x反褶積能夠降低噪聲強度，一定程度上恢復同相軸的結構，但仍有較多噪聲殘留。而本文方法處理后，在地震同相軸處只有輕微的噪聲，能夠有效去除地震數據中的噪聲成分，使地震信號更清晰，減少了地震數據中的干擾，提高了地震數據的信噪比。

圖6是三種方法在加入不同強度噪聲下的信噪比曲線。由于本文方法與傳統K‐SVD方法去噪性能幾乎相同，區別在于計算效率，所以著重與f?x反褶積方法比較。整體可以看出在不同標準差的噪聲強度下，使用本文方法去噪后，地震數據的信噪比明顯高于使用f?x反褶積方法。隨著噪聲強度的增加，兩種方法的信噪比都有所下降，但本文方法的信噪比下降速度相對較慢，表現出更強的魯棒性。通過對比，在一定噪聲強度下，本文提出的去噪方法明顯優于f?x反褶積算法，這表明本文方法能夠更有效地去除不同強度的噪聲，從而獲得更清晰、更準確的地震數據。

圖7為f?x反褶積方法和本文方法去噪后的殘差對比。圖7a中可以觀察到明顯的有效信號殘留，表明f?x反褶積在去除噪聲的同時也造成了有效信號的缺失。為了更加全面地評估不同去噪方法的效果，引入局部相似性[25]概念。

圖8為f?x反褶積和本文方法處理后原始數據與去噪數據的局部相似性對比。為了比較結果的公平性，統一了相似性的色標范圍，能量越強，表明數據的相似性越好，去噪性能越強。在相同采樣時間和相同地震道上，使用f?x反褶積去噪后的能量強度明顯低于本文方法，這表明本文方法處理后的地震數據受到噪聲干擾小于f?x反褶積，反映出更優秀的去噪性能。

圖9a為一組隨機選擇的64個訓練塊，展示了字典原子在更新前、后的顯著變化。在未更新前，字典原子呈現出一種剛性特征，結構單一，缺乏對數據間復雜線性關系的準確反映。這意味著原字典結構的局限性，無法充分捕捉到數據中的細微差異和多樣性。更新后的字典原子之間展現出更豐富的變化，使字典更能適應數據的多樣化表現特征。這種靈活性的增加意味著更新的字典能夠更好地捕捉數據之間的線性關系，從而提高了字典在數據表示方面的性能，具有更強的適應性，能夠更好地反映數據中的細微變化和復雜關系。

圖10為兩種方法對單炮記錄去噪后的f?k譜[26]對比，圖10a、圖10b為原始數據與噪聲數據對應的f?k譜，圖10c、圖10d為f?x反褶積和本文方法去噪后對應的f?k譜。由圖可見f?x反褶積處理后能量強度較低，對有效信號的恢復較弱，本文方法的去噪效果更好。

5 實際數據試驗

圖11為二維地震數據在不同信息保留率下的可視化結果，清晰展示了PCA在數據降維中的效果。從表1可以看出，在保留不同百分比信息的情況下，降維后得到的主特征數量也存在較大的差異，特別是在保留99.99%的地質信息時，可以觀察到主特征數量相對減少，這意味著在保留極大比例數據重構相似性的同時，舍棄了大量的非主要特征。具體來說，在保留99.99%信息的情況下，成功地舍棄166個非主要特征，這些特征不再參與SVD的計算過程，從而顯著減少計算的復雜度，提高計算效率。這種降維后的特征選擇和信息保留策略在大規模數據處理中提供了更高效和可行的解決方案。

將使用PCA降維重構后的地震數據保存為numpy數組并對其去噪，從計算效率和去噪性能兩方面比較兩種方法的區別。圖12a為實際采集的疊前地震數據，該數據集共有100道，每道1500個采樣點，采樣間隔為1ms。受采集現場條件影響，噪聲淹沒了部分有效信號。對該數據使用f?x反褶積、傳統K‐SVD方法以及本文方法去噪。

由圖12可見，使用f?x反褶積方法處理噪聲后，地震同相軸的清晰度有所提高，有效信號也在一定程度上得到恢復。傳統K‐SVD方法和本文方法在去噪性能上比f?x反褶積方法表現更優秀，能夠在提高信噪比的同時保留地震數據中的有效信號。而傳統K‐SVD方法和本文方法之間的主要區別在于數據降維的實現方式。因此，需要對這兩種方法在計算性能方面進行深入比較，以便更全面地評估它們在地震數據處理中的表現。

從表2可知本文方法的去噪效果略好于傳統KSVD，所以并不足以體現出本文方法與傳統K‐SVD的區別，下面將從字典學習中的字典原子數以及最大迭代次數等方面評價兩種方法的計算效率。

表3是固定迭代次數的情況下比較不同字典原子數。字典原子數增加會增加算法的復雜度，算法的復雜度是影響計算效率的主要因素。數據規模的大小能直接影響程序的計算量和內存占用，處理大規模數據集可能會導致內存不足、計算速度變慢等問題。同樣的，增加字典原子個數可能會導致字典過度擬合訓練數據的風險增加。如果字典原子數量過多，導致字典過于靈活，從而捕捉到數據中的噪聲或不必要的細節，會影響字典的泛化能力，并在應用新數據時產生較差的效果。

表4是固定字典原子的情況下對不同迭代次數計算時間進行比較。字典更新是字典學習過程中最耗時的部分，每次對字典更新都涉及到大量的矩陣運算和優化操作，增加迭代次數會導致算法執行更多次更新操作，增加總體的計算時間。字典更新過程中還需要存儲和處理字典矩陣、稀疏表示等字典結構，這會占用一定的內存空間，增加迭代次數會導致需要存儲更多的中間結果和臨時變量，從而增加內存的占用。由此可知，在不斷增加字典原子數的情況下，本文方法的計算速度明顯快于傳統K‐SVD，這說明本文方法在使用PCA階段有效地實現了數據降維的工作且保留了有效信息，也切實地說明了在K‐SVD數據梳理階段添加PCA的可行性。

在地震數據去噪研究中，頻率—波數譜被視為評估不同去噪方法效果的重要指標。圖13為不同數據的f?k譜對比。在原始數據上幾乎看不到有效信號的能量，這表明有效信號受噪聲干擾嚴重，經過f?x反褶積去噪后，在中心軸附近的能量強度得到了恢復，但效果比傳統K‐SVD和本文方法要差。

為了更清晰地展示f?x反褶積與本文方法在復雜地質數據上去噪效果的差異，比較兩種方法的去噪結果與去除的噪聲之間的局部相似性（圖14），相似性的振幅值越低，表示對原始信號的損失越小。可以觀察到，使用f?x反褶積方法去噪后的振幅值高于本文方法，這表明去噪之后有效信號的損失較大；本文方法去噪后局部相似性能量較弱，表明使用本文方法去噪后有效信號的損失較小，更好地保留了完整的地震信號。

6 結論

地震數據去噪直接影響數據質量和后續分析結果的準確性。本文提出的基于PCA與K‐SVD的地震數據去噪方法在三種方法的數據對比中取得了顯著效果。K‐SVD作為一種基于字典學習的去噪方法，其計算復雜度較高，特別是在處理大規模數據時更明顯，通過引入PCA進行數據降維，能夠有效減少計算量，節約大量計算資源，提高算法的實用性和適用性。數據試驗結果驗證了本文方法的可行性，為地震數據處理提供了更高效的解決方案。盡管信噪比的提升有限，但與大幅降低的計算成本相比，這是一個值得考慮的權衡結果。

從試驗結果來看，該方法仍存在缺陷，去噪后提升信噪比不夠顯著，特別是在處理復雜的地震數據時，與傳統K‐SVD方法的去噪效果幾乎相同。在未來的研究中，可以在本文方法的基礎上從字典的初始化入手，不再采用隨機初始化方式，而是使用預訓練的深度學習模型，如自編碼器等來生成初始字典，以達到獲得更高信噪比的去噪結果。