侯思安，張峰，李向陽

中國石油大學(xué)(北京)油氣資源與工程國家重點實驗室，北京 102249

0 引言

地震數(shù)據(jù)重建是提高資料品質(zhì)和降低勘探綜合成本的關(guān)鍵技術(shù)。經(jīng)典的數(shù)據(jù)重建算法主要有預(yù)測濾波算法(Prediction Filter)和促稀疏反演算法(Sparsity Promoting Inversion)：預(yù)測濾波算法是將含有缺失道的地震數(shù)據(jù)從時間—空間域變換到頻率—空間域，然后根據(jù)數(shù)據(jù)在頻率—空間域中表現(xiàn)出來的線性特征進行數(shù)據(jù)重建[1];而促稀疏反演算法是假設(shè)原始的地震信號在變換域中是稀疏的，當(dāng)數(shù)據(jù)存在缺失時稀疏性會變?nèi)酰虼丝梢酝ㄟ^求解一個l1模規(guī)則化的最優(yōu)化問題來恢復(fù)信號的稀疏特征實現(xiàn)數(shù)據(jù)重建[2-3]。這些算法的一個共同點是都依賴于解析數(shù)學(xué)變換如快速傅里葉變換(Fast Fourier Transform)[4]，拉東變換(Radon Transform)[5]，曲波變換(Curvelet Transform)[6]和地震小波變換(Seislet Transform)[7]等。但是受限于計算機的性能和實際地震資料的復(fù)雜性，解析數(shù)學(xué)變換已經(jīng)很難滿足石油工業(yè)對數(shù)據(jù)處理精度的需求。因此學(xué)術(shù)界將研究的重心轉(zhuǎn)向了機器學(xué)習(xí)(Machine Learning)算法，希望能通過機器學(xué)習(xí)達到更好的數(shù)據(jù)重建效果。

就地震數(shù)據(jù)重建這個問題而言，主要的機器學(xué)習(xí)算法有稀疏字典學(xué)習(xí)算法(Sparsity Dictionary Learning)和低秩矩陣分解算法(Low-rank Matrix Factorization)。與數(shù)學(xué)變換算法相似，稀疏字典學(xué)習(xí)算法是通過輸入數(shù)據(jù)自適應(yīng)的生成一組變換基函數(shù)，在該基函數(shù)下原始的輸入數(shù)據(jù)具有稀疏特征。目前計算精度最高的稀疏字典學(xué)習(xí)算法是k-SVD算法[8-9]，在使用該算法進行地震數(shù)據(jù)插值時[10-11]：第一步，對原始數(shù)據(jù)進行Patch處理[12]并初始化一組基函數(shù)(在k-SVD算法中基函數(shù)被稱為字典)；第二步，求解l1模規(guī)則化問題計算稀疏系數(shù)；第三步，分別對基函數(shù)的每一列進行更新，更新每一列時要對殘差矩陣進行一次SVD分解，并用分解的第一列左特征向量矩陣更新基函數(shù)，用右特征向量矩陣的第一行和第一個特征值的乘積更新系數(shù)矩陣的對應(yīng)行；第四步，循環(huán)步驟二至步驟四若干次直至收斂為止?？梢钥闯霎?dāng)基函數(shù)具有k列時，該算法每一次迭代都要進行k次SVD分解，巨大的計算量限制了該算法在實際資料處理中的應(yīng)用。數(shù)據(jù)驅(qū)動的緊致框架(Data-Driven Tight Frame，DDTF)[13]則是一種改進的稀疏字典學(xué)習(xí)算法，該算法對基函數(shù)施加了緊致性約束條件，在每次迭代時僅需要一次SVD分解，極大地提高了數(shù)據(jù)處理的效率[14-15]。低秩矩陣分解算法則假設(shè)實際的地震信號具有低秩特征，但是當(dāng)數(shù)據(jù)中存在缺失或噪音時，信號的秩會提高，因此可以設(shè)計一個降秩的優(yōu)化算法實現(xiàn)地震數(shù)據(jù)重建[16]，常用于地震數(shù)據(jù)處理的減秩算法有阻尼減秩算法(Damped Rank-reduction Method)[12]，特征值收縮算法(Singular Value Shrinkage)[17]和多通道奇異譜分析算法(Multichannel Singular Spectrum Analysis, MSSA)[18]等。此外實際用于地震數(shù)據(jù)插值的都是4D或5D數(shù)據(jù)體，針對這種高維信號可以采用構(gòu)建Hankel矩陣對數(shù)據(jù)降維[17]或直接應(yīng)用張量分解算法進行求解[19-20]。

相比較而言低秩矩陣分解算法在求解最優(yōu)化問題時不用處理l1模規(guī)則化，也不是必須要用SVD分解，因此在計算效率上具有一定的優(yōu)勢。但是在求解該算法時需要考慮多個最優(yōu)化參數(shù)，多數(shù)情況下是通過設(shè)置不同的參數(shù)進行數(shù)據(jù)處理然后人工選取較優(yōu)的結(jié)果，這樣做會增大數(shù)據(jù)處理的計算量，且處理精度也難以保障。針對這一問題，Salakhutdinov和Mnih在求解推薦系統(tǒng)(Recommendation System)的矩陣分解問題時提出了貝葉斯概率矩陣分解(Bayesian Probabilistic Matrix Factorization，BPMF)算法[21]，該算法可以對基函數(shù)和系數(shù)的均值和方差等統(tǒng)計學(xué)參數(shù)進行隨機模擬，通過計算不同參數(shù)的概率密度函數(shù)自適應(yīng)的選取最優(yōu)結(jié)果，Net fl ix問題測試表明該方法具有良好的應(yīng)用效果。本文將該算法引入到地震數(shù)據(jù)重建問題中，大量的數(shù)值測試表明該算法可以提高數(shù)據(jù)處理的精度和穩(wěn)定性。本文首先介紹了地震數(shù)據(jù)矩陣分解的概率解釋和BPMF算法的原理；然后，通過單道子波、合成地震記錄和實際資料對該方法進行測試；最后，對算法進行總結(jié)和展望。

1 方法原理

1.1 基于概率矩陣分解的地震數(shù)據(jù)插值算法

基于概率矩陣分解(Probabilistic Matrix Factorization，PMF)[22]的地震數(shù)據(jù)插值算法通常假設(shè)實際地震信號矩陣是低秩的，并且可以表示為兩個矩陣的乘積形式：

其中，X ∈ Rn×m表示含有隨機缺失的地震數(shù)據(jù)，n和m分別表示矩陣的行數(shù)和列數(shù)；M ∈ Rn×k和 A ∈ Rk×m分別表示基函數(shù)和對應(yīng)的系數(shù)，理論上k等于數(shù)據(jù)X的秩；E表示隨機噪音矩陣。

由于基函數(shù)M和系數(shù)A都是未知的，很自然的假設(shè)的這兩個矩陣的元素都是符合Gaussian分布的：

其中，N(y|?, σ2)表示y符合均值為?，方差為σ2的Gaussian分布；μi,j和ai,j分別表示M和A的元素；和分別表示M和A的方差；∏表示乘積運算符。

因為隨機噪音E也是滿足Gaussian分布的，所以：

其中，σE2表示隨機噪音E的方差。

根據(jù)概率密度函數(shù)的平移關(guān)系有：

其中，表示重建地震數(shù)據(jù)X*= MA的元素；上標(biāo)Ii,j用于標(biāo)記數(shù)據(jù)的缺失信息，當(dāng)該元素缺失時Ii,j=0，否則Ii,j=1。

因此基于PMF的地震數(shù)據(jù)重建等價于在采集數(shù)據(jù)X已知時，求取M和A的最大后驗概率：

因為X是已知的，因此概率密度函數(shù)P(X)是一個常數(shù)，所以有：

對公式(7)求對數(shù)并帶入概率密度函數(shù)(2)、(3)、(5)和高斯分布函數(shù)，可以得到：

其中，c是用于平衡方程(8)左右兩端的一個常數(shù)。

根據(jù)公式(8)，在采集數(shù)據(jù)X已知時，求取基函數(shù)M和系數(shù)A的最大后驗概率等價于求解如下的最優(yōu)化問題：

圖1 貝葉斯概率矩陣分解示意圖Fig. 1 The graphical model for Bayesian Probabilistic Matrix Factorization

其中，xi,j和分別表示采集得到的地震數(shù)據(jù)和重建的地震數(shù)據(jù)；表示最優(yōu)化權(quán)系數(shù)，和表示基函數(shù)M、系數(shù)A和隨機噪音E的方差；表示l2模規(guī)則化；?表示采集數(shù)據(jù)的觀測系統(tǒng)。

公式(9)描述了基于PMF的地震數(shù)據(jù)插值算法，在應(yīng)用該算法進行數(shù)據(jù)重建時需要提前知道規(guī)則化參數(shù)λM和λA，而這兩個參數(shù)又和基函數(shù)M、系數(shù)A和隨機噪音E的方差和相關(guān)。其中隨機噪音方差可以通過分析地震數(shù)據(jù)進行估算，但是和通常是無法獲取的。一種普遍的做法是嘗試用不同的參數(shù)進行數(shù)據(jù)處理，然后選取其中效果最好的作為最后的處理結(jié)果，但是這樣做的計算量是比較大的。因此，本文將貝葉斯概率矩陣分解算法引入到地震數(shù)據(jù)插值問題中，基于該算法可以有效地解決規(guī)則化參數(shù)不確定的問題。

1.2 基于貝葉斯概率矩陣分解的地震數(shù)據(jù)插值算法

貝葉斯概率矩陣分解(Bayesian Probabilistic Matrix Factorization，BPMF)的核心原理是假設(shè)決定基函數(shù)M的系數(shù)A分布特征的超參數(shù)和是符合Gaussian-Wishart分布的，通過馬爾科夫蒙托卡羅(Markov Chain Monte Carlo，MCMC)方法對不同參數(shù)的結(jié)果進行隨機模擬，并計算不同參數(shù)對應(yīng)的重建結(jié)果的概率密度函數(shù)自適應(yīng)的選取最優(yōu)結(jié)果。BPMF的示意圖如圖1所示。

根據(jù)前文的闡述，在BPMF算法中基函數(shù)M和系數(shù)A滿足如下的高斯分布：

其中，Mi,:表示基函數(shù)矩陣M的第i行，A:,j表示系數(shù)矩陣A的第j列，和是決定基函數(shù)M和系數(shù)A分布特征的超參數(shù)，這兩個超參數(shù)都符合Gaussian-Wishart分布：

其中，ξ0等于0，v0等于基函數(shù)的列數(shù)k，W0為大小是k×k的單位矩陣；W(Λ|W0,v0)表示W(wǎng)ishart分布：

其中，c表示歸一化參數(shù)，Tr表示矩陣的跡。

在進行缺失地震數(shù)據(jù)重建時，可以利用貝葉斯邊緣化處理(Marginalizing)來提高插值算法的精度?；驹砭褪菍θ魏慰赡艹霈F(xiàn)的超參數(shù){ΘM,ΘA}分別計算數(shù)據(jù)重建的結(jié)果，并對所有的結(jié)果進行加權(quán)求和，求和的權(quán)系數(shù)就是每個超參數(shù)對應(yīng)的概率密度函數(shù)，因此基于BPMF的地震數(shù)據(jù)重建可以表示為如下的一個積分：

直接求解積分函數(shù)(15)幾乎是不可能的，因此Salakhutdinov和Mnih提出了基于馬爾科夫蒙托卡羅(Markov Chain Monte Carlo，MCMC)的求解方法。在MCMC的框架下，積分函數(shù)(15)可以近似如下的求和函數(shù)：

其中，表示馬爾科夫隨機采樣序列，序列的長度是K。

本文使用Gibbs方法對公式(16)進行隨機采樣。由于地震數(shù)據(jù)插值僅需要概率最大的基函數(shù)M和系數(shù)A，所以實際計算時不需要完整的求解概率密度函數(shù)(16)，而是通過多次迭代使得其達到穩(wěn)定即可，那么基于BPMF的地震數(shù)據(jù)重建流程如下：

第一步，初始化基函數(shù)M(0)和系數(shù)A(0)，初始化時可以采用其他算法求得的結(jié)果或直接采用隨機函數(shù)生成得到，對輸入地震數(shù)據(jù)進行Patch處理[12];

第 二 步， 更 新 超 參 數(shù)和由于基函數(shù)M和系數(shù)A在這一步中是已知的，超參數(shù)和的概率密度函數(shù)表示為：

其中，

同理：

其中，

第三步，更新基函數(shù)M，根據(jù)貝葉斯公式：

因為此時采集數(shù)據(jù)X，系數(shù)矩陣A和超參數(shù)都是已知的，所以條件概率密度函數(shù)P(X| A,ξM,ΛM)是一個常數(shù)，那么：

又因為基函數(shù)矩陣M和系數(shù)矩陣A是無關(guān)的：

所以：

第四步，更新系數(shù)A，與第三步同理可以得：

最后，重復(fù)步驟二至步驟四直至基函數(shù)M和系數(shù)A達到穩(wěn)定為止。

數(shù)據(jù)重建的流程如圖2所示。

1.3 貝葉斯概率矩陣分解算法的穩(wěn)定性分析

BPMF算法是一種隨機缺失數(shù)據(jù)的恢復(fù)算法，要對缺失數(shù)據(jù)進行精確的恢復(fù)除了要求原始數(shù)據(jù)是低秩的，還要滿足一定的采樣條件。Candes和Recht對這一問題進行了深入研究，當(dāng)隨機采樣滿足如下的條件時即可用低秩算法進行數(shù)據(jù)重建[16]，具體條件如下：

其中，e表示隨機采樣的數(shù)據(jù)點數(shù)；f =max(m, n)表示原始數(shù)據(jù)矩陣行數(shù)或列數(shù)的大值；r表示原始數(shù)據(jù)矩陣的秩；C表示一個數(shù)值常數(shù)。根據(jù)采樣條件，當(dāng)原始數(shù)據(jù)的秩r比較小時(Cf65r log f＜ m × n )，可以用隨機采樣的數(shù)據(jù)精確地恢復(fù)出原始數(shù)據(jù)；但是當(dāng)r比較大時(Cf65r log f＞ m × n )，即便知道了所有有效信息，低秩算法也會造成原始數(shù)據(jù)的缺失。

圖3表示用隨機生成的低秩矩陣進行的BPMF方法穩(wěn)定性測試。在測試中，原始數(shù)據(jù)矩陣是一個100× 100的方陣，矩陣的秩從1逐漸增加到41，采樣率從0%逐漸增加到100%，數(shù)據(jù)重建的效果使用峰值信噪比(公式25)進行評估，并假設(shè)信噪比大于15時為重建成功，否則為失敗。結(jié)果表明BPMF算法是一種精確穩(wěn)定的低秩恢復(fù)算法，僅需要不到20%的隨機采樣數(shù)據(jù)就可以對低秩數(shù)據(jù)進行恢復(fù)。

基于BPMF的地震數(shù)據(jù)插值算法，在理論上可以有效的減少數(shù)據(jù)處理過程中對優(yōu)化參數(shù)的依賴，提高數(shù)值計算的穩(wěn)定性和精度，并用一個隨機數(shù)矩陣驗證了方法的穩(wěn)定性。但是真正的地震數(shù)據(jù)并不是完全隨機缺失的(表現(xiàn)為在時間方向連續(xù)缺失，在空間方向隨機缺失)，并且地震數(shù)據(jù)也不是完全低秩數(shù)據(jù)，這些問題可能會影響基于BPMF算法的地震數(shù)據(jù)重建效果。本文將通過合成地震數(shù)據(jù)和實際資料對本方法進行測試。

圖2 貝葉斯概率矩陣分解算法計算流程Fig. 2 Workflow of Bayesian probabilistic matrix factorization algorithm

圖3 隨機數(shù)矩陣重建測試：白色表示重建成功；黑色表示重建失敗Fig. 3 Illustration of random matrix recovery：white and black represents the successful recovery and failed recovery,respectively

2 方法測試

2.1 單道地震數(shù)據(jù)測試

首先通過隨機缺失的單道地震子波數(shù)據(jù)對方法進行測試。單道記錄如圖4a所示，道集總共有401個采樣點，采樣間隔為1 ms，4個獨立的子波主頻分別為 15 Hz、30 Hz、15 Hz和 30 Hz，振幅分別為 1.0、1.0、3.0和3.0，數(shù)據(jù)總計有201道隨機缺失。數(shù)據(jù)重構(gòu)結(jié)果如圖4b所示，4條不同的線段分別表示Patch大小為4、8、12和16的重構(gòu)子波(Patch處理類似于地震數(shù)據(jù)的時窗概念，是一種常用的圖像處理手段，具體定義可參考文獻12)，圖4c和圖4d分別表示了0～200 ms和 200～400 ms的局部放大圖。根據(jù)Patch=4和Patch=8的重建結(jié)果，地震數(shù)據(jù)重建的Patch尺寸不能太小，否則會出現(xiàn)局部數(shù)據(jù)的采樣不足，無法實現(xiàn)數(shù)據(jù)重建；同時對比Patch=12和Patch=16的重建結(jié)果，當(dāng)Patch過大時，重建結(jié)果容易過度平滑；并且對比不同的子波結(jié)果發(fā)現(xiàn)，本文方法更容易使主頻較高、振幅較強的數(shù)據(jù)產(chǎn)生平滑，因此在實際數(shù)據(jù)處理時需要嘗試不同的Patch參數(shù)然后選取最優(yōu)的計算結(jié)果。

2.2 合成地震記錄測試

本文通過合成地震記錄對提出方法進行測試。測試數(shù)據(jù)如圖5a所示，原始數(shù)據(jù)具有201個地震道，每個地震道具有201個采樣點，空間和時間采樣步長分別為10 ms和7 ms，隨機選取其中的40道為缺失數(shù)據(jù)，如圖5b所示。在數(shù)據(jù)重建時選取基函數(shù)M的列數(shù)k等于20，規(guī)則化參數(shù)λM=λA=0.01。圖5c、圖5d和圖5e分別表示基于Curvelet方法[23]，PMF方法和BPMF方法的重建結(jié)果；圖5f、圖5g和圖5h分別表示對應(yīng)的重建數(shù)據(jù)殘差。通過這一組數(shù)據(jù)對比表明，基于矩陣分解原理的PMF算法和BPMF算法比經(jīng)典的Curvelet算法能更加有效地恢復(fù)地震數(shù)據(jù)，重建結(jié)果的殘差明顯減少。但是當(dāng)數(shù)據(jù)在一定的區(qū)域內(nèi)存在大量的缺失時，PMF算法的重建效果出現(xiàn)了較為嚴(yán)重的下降，如圖5d的右側(cè)所示。相比較而言，在相同的區(qū)域BPMF算法能較好的提高數(shù)據(jù)重建的效果，差剖面中沒有明顯的數(shù)據(jù)畸變，如圖5h所示。

圖4 單道地震記錄測試：(a)含有隨機缺失的單道數(shù)據(jù)，從左到右的子波依次為：主頻15 Hz振幅1.0，主頻30 Hz振幅1.0、主頻15 Hz振幅3.0和主頻30 Hz振幅3.0；(b)重建的單道數(shù)據(jù)；(c)0～200 ms局部放大圖；(d)200～400 ms局部放大圖Fig. 4 Single trace test：(a)Random missing single trace data, wavelets parameters(from left to right): 15 Hz amp=1.0, 30 Hz amp=1.0, 15 Hz amp=3.0 and 30 Hz amp=3.0; (b) Reconstructed single trace data; (c) Zoom of reconstructed data(0～200ms); (d)Zoom of reconstructed data(200～400 ms)

圖6展示了基于BPMF算法進行地震數(shù)據(jù)矩陣分解得到的基函數(shù)M，其中圖6a表示BPMF的基函數(shù)，圖6b表示離散余弦變換(Discrete Cosine Transform，DCT)的基函數(shù)，每一個單獨的數(shù)據(jù)塊代表基函數(shù)的一列?？梢钥闯鯞PMF算法的基函數(shù)的每一個數(shù)據(jù)塊都非常類似于地震數(shù)據(jù)的同相軸，區(qū)別在于振幅和相位的不同。而DCT基函數(shù)是具有解析數(shù)學(xué)表達式的傳統(tǒng)數(shù)學(xué)變換，每個數(shù)據(jù)塊都是固定的，與輸入數(shù)據(jù)無關(guān)。通過對比可以看出矩陣分解算法可以根據(jù)輸入數(shù)據(jù)自適應(yīng)的構(gòu)建基函數(shù)M，具有特征提取的功能。

圖5 合成地震記錄測試：(a)原始數(shù)據(jù)；(b)含隨機缺失數(shù)據(jù)；(c)Curvelet重建數(shù)據(jù)；(d)PMF重建數(shù)據(jù)；(e)BPMF重建數(shù)據(jù)；(f)Curvelet重建差剖面；(g)PMF重建差剖面；(h)BPMF重建差剖面Fig. 5 Synthetic data reconstruction test: (a) Original data; (d) Random missing data; (c) Reconstruction data via Curvelet; (d)Reconstruction data via PMF; (e) Reconstruction data via BPMF; (f) Residual of Curvelet reconstruction data; (g) Residual of PMF reconstruction data; (h) Residual of BPMF reconstruction data

為了更進一步的對比PMF算法和BPMF算法，測試了隨機缺失道數(shù)從10道遞增到100道，基函數(shù)M的列數(shù)k分別為15、20和25，以及規(guī)則化參數(shù)λM和λA分別為0.01、0.02和0.03的情況，測試結(jié)果如圖7所示，峰值信噪比的計算公式為：

由于地震數(shù)據(jù)的復(fù)雜性以及信號噪音和缺失的干擾，這些重建參數(shù)是難以直接獲取的，也沒有一個明確的選取準(zhǔn)則。但是通過這一組數(shù)據(jù)對比，當(dāng)k從15逐漸增加到25并且規(guī)則化參數(shù)λM和λA從0.03逐漸減少到0.01時，PMF算法數(shù)據(jù)重建的精度在不斷提高，但是也出現(xiàn)了較大幅度的波動，這說明算法的穩(wěn)定性也有所降低，但是BPMF算法始終保持了一個較高的計算精度和穩(wěn)定性。

圖6 基函數(shù)示意圖：(a)BPMF基函數(shù)示意圖；(b)DCT基函數(shù)示意圖Fig. 6 Illustration of basis matrix: (a) Basis matrix of BPMF; (b) Basis matrix of DCT

圖7 PMF和BPMF對 比 測 試：(a)k =15測 試 結(jié) 果；(b)k =20測試結(jié)果；(c)k=25測試結(jié)果Fig. 7 Comparison PMF with BPMF: (a) Test result of k=15; (a) Test result of k =20; (a) Test result of k=25

圖8 含噪音合成地震記錄測試：(a)原始數(shù)據(jù)；(b)含隨機缺失數(shù)據(jù)；(c)BPMF重建數(shù)據(jù)；(d)BPMF重建差剖面Fig. 8 Synthetic noisy data reconstruction test: (a) Original data; (b) Random missing data; (c) Reconstruction data via BPMF;(d) Residual of BPMF reconstruction data

圖9 含噪音合成地震記錄全部測試Fig. 9 Overview of synthetic noisy data reconstruction test

圖8和圖9則測試了在含有隨機噪音的情況下，BPMF算法的重建效果，其中圖8展示了隨機25道缺失的重建效果，圖9展示了隨機缺失道數(shù)從10道遞增到100道時的全部重建結(jié)果。通過這一組測試可以看出當(dāng)數(shù)據(jù)中存在隨機噪音時，BPMF算法依然可以穩(wěn)定高效的對缺失地震數(shù)據(jù)進行重建。

2.3 實際地震記錄測試

最后通過實際資料測試BPMF地震數(shù)據(jù)重建算法。測試使用的是海上OBC數(shù)據(jù)的一個剖面，該剖面有101個地震道，每個地震道具有2501個采樣點，時間采樣間隔為2 ms，截取其中1000 ms至3000 ms的一段數(shù)據(jù)用于測試，如圖10a所示；測試時隨機選取了其中的38道數(shù)據(jù)為缺失道，如圖10b所示；基于BPMF算法的重建結(jié)果(k=25)和差剖面分別如圖10c和圖10d所示；圖11則分別展示了第22道、第74道和第97道重建數(shù)據(jù)和原始數(shù)據(jù)的單道記錄。通過對比重建的剖面和單道記錄，大部分缺失的地震數(shù)據(jù)被精確的恢復(fù)出來了，特別是對于主頻較低的深層信號，差剖面中幾乎不含殘留的有效信號，這說明本文提出的BPMF算法可以有效的對缺失地震數(shù)據(jù)進行重建。但是對能量較強的信號，數(shù)據(jù)恢復(fù)的結(jié)果中出現(xiàn)了比較明顯的平滑，恢復(fù)信號的振幅能量降低了，這與單道記錄測試結(jié)果是一致的，也是今后研究需要解決的問題。圖12則展示了基于BPMF算法計算的基函數(shù)，該基函數(shù)依然體現(xiàn)出了非常明顯的地震數(shù)據(jù)同相軸特征。

圖10 實際地震數(shù)據(jù)測試：(a)原始地震數(shù)據(jù)；(b)含缺失道地震數(shù)據(jù)；(c)BPMF重建地震數(shù)據(jù)；(d)BPMF 重建差剖面Fig. 10 Real seismic data test: (a) Original seismic data; (b) Missing seismic data; (c) Reconstructed seismic data via BPMF;(d) Residual of BPMF reconstruction data

圖11 實際地震數(shù)據(jù)測試單道記錄：(a)第22道記錄；(b)第74道記錄；(c)第97道記錄Fig. 11 Trace data of real seismic data test: (a) No. 22 trace data; (b) No. 74 trace data; (c) No. 97 trace data

3 結(jié)論與展望

圖12 實際地震數(shù)據(jù)基函數(shù)Fig. 12 Basis matrix of real seismic data

矩陣分解算法是在地震數(shù)據(jù)信號處理中非常具有研究價值的工業(yè)應(yīng)用前景的機器學(xué)習(xí)算法，在應(yīng)用該算法時需要求解一個較為復(fù)雜的最優(yōu)化問題，影響最終處理結(jié)果的有基函數(shù)M、系數(shù)A和隨機噪音E的方差和以及基函數(shù)的列數(shù)k(理論上k應(yīng)該等于實際地震數(shù)的秩)等參數(shù)。本文針對和的選取問題開展了研究，通過引入馬爾科夫蒙托卡羅方法對不同的參數(shù)進行隨機模擬并計算相應(yīng)的概率密度函數(shù)自適應(yīng)的選取最優(yōu)結(jié)果，合成地震數(shù)據(jù)和實際資料測試表明該算法在計算精度和穩(wěn)定性方面均有所提高。

但是該算法也存在一些需要研究和解決的問題。首先，當(dāng)數(shù)據(jù)中能量較強的同相軸，本文算法會對結(jié)果產(chǎn)生一定的平滑作用，導(dǎo)致重建數(shù)據(jù)失真，初步分析認為這與Patch處理有直接關(guān)系，但是如何解決這個問題還需要深入研究；然后貝葉斯算法可以根據(jù)輸入數(shù)據(jù)自適應(yīng)的生成一個基函數(shù)，該基函數(shù)表現(xiàn)出了非常明顯的地震數(shù)據(jù)特征，是否可以利用該基函數(shù)以及如何利用該基函數(shù)還是一個需要深入研究的問題。

[1] SPITZ S. Seismic trace interpolation in the F-X domain[J]. Geophysics, 1991, 56(6): 785-794.

[2] BAI L S, LIU Y K, LU H Y, et al. Curvelet-domain joint iterative seismic data reconstruction based on compressed sensing[J]. Chinese Journal of Geophysics, 2014, 57(9): 2937-2945.

[3] HERRMANN F J, HENNENFENT G. Non-parametric seismic data recovery with curvelet frames[J]. Geophysical Journal of the Royal Astronomical Society, 2010, 173(1): 233-248.

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

[5] WANG J, NG M, PERZ M. Seismic data interpolation by greedy local Radon transform[J]. Geophysics, 2010, 75(6): WB225-WB234.

[6] CANDèS E, DEMANET L, DONOHO D, et al. Fast discrete curvelet transforms[J]. multiscale modeling & simulation, 2006, 5(3):861-899.

[7] FOMEL S, LIU Y. Seislet transform and Seislet frame[J]. Geophysics, 2010, 75(3): V25-V38.

[8] AHARON M, ELAD M, BRUCKSTEIN A. K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[J].IEEE Transactions on Signal Processing, 2006, 54(11): 4311-4322.

[9] ELAD M, AHARON M. Image denoising via sparse and redundant representations over learned dictionaries[J]. IEEE Transactions on Image Processing, 2006, 15(12): 3736-3745.

[10] BECKOUCHE S, MA J. Simultaneous dictionary learning and denoising for seismic data[J]. Geophysics, 2014, 79(3): A27-A31.

[11] CHEN Y. Fast dictionary learning for noise attenuation of multidimensional seismic data[J]. Geophysical Journal International, 2017,209(1): 21-31.

[12] MA J. Three-dimensional irregular seismic data reconstruction via low-rank matrix completion[J]. Geophysics, 2013, 78(5):V181-V192.

[13] CAI J F, JI H, SHEN Z, et al. Data-driven tight frame construction and image denoising[J]. Applied & Computational Harmonic Analysis, 2014, 37(1): 89-105.

[14] LIANG J, MA J, ZHANG X. Seismic data restoration via data-driven tight frame[J]. Geophysics, 2014, 79(3): V65-V74.

[15] YU S, MA J, ZHANG X, et al. Interpolation and denoising of high-dimensional seismic data by learning a tight frame[J]. Geophysics,2015, 80(5): V119-V132.

[16] CANDèS E J, RECHT B. Exact matrix completion via convex optimization[J]. Foundations of Computational Mathematics, 2008, 9(6):717.

[17] CHEN K, SACCHI M D. Robust reduced-rank fi ltering for erratic seismic noise attenuation[J]. Geophysics, 2015, 80(1): V1-V11.

[18] CHEN Y, ZHANG D, JIN Z, et al. Simultaneous denoising and reconstruction of 5-D seismic data via damped rank-reduction method[J].Geophysical Journal International, 2016, 206(3): 1695–1717.

[19] KREIMER N, SACCHI M D. A tensor higher-order singular value decomposition for prestack seismic data noise reduction and interpolation[J]. Geophysics, 2012, 77(3): V113-V122.

[20] GAO J, STANTON A, SACCHI M D. Parallel matrix factorization algorithm and its application to 5D seismic reconstruction and denoising[J]. Geophysics, 2015, 80(6): V173-V187.

[21] SALAKHUTDINOV R, MNIH A. Bayesian probabilistic matrix factorization using Markov chain Monte Carlo[C]. International Conference on Machine Learning, Helsinki, Fabianinkatu, ACM, 2008: 880-887.

[22] SALAKHUTDINOV R, MNIH A. Probabilistic matrix factorization[C]. International Conference on Neural Information Processing Systems. Vancouver, Canada, 2007: 1257-1264.

[23] DAUBECHIES I, DEFRISE M, DE MOL C. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[J]. Communications on Pure & Applied Mathematics, 2004, 57(11): 1413-1457.