基于貝葉斯概率矩陣分解的地震數據重建算法

侯思安，張峰，李向陽

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

0 引言

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

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

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

1 方法原理

1.1 基于概率矩陣分解的地震數據插值算法

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

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

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

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

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

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

根據概率密度函數的平移關系有：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

在進行缺失地震數據重建時，可以利用貝葉斯邊緣化處理(Marginalizing)來提高插值算法的精度。基本原理就是對任何可能出現的超參數{ΘM,ΘA}分別計算數據重建的結果，并對所有的結果進行加權求和，求和的權系數就是每個超參數對應的概率密度函數，因此基于BPMF的地震數據重建可以表示為如下的一個積分：

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

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

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

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

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

其中，

同理：

其中，

第三步，更新基函數M，根據貝葉斯公式：

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

又因為基函數矩陣M和系數矩陣A是無關的：

所以：

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

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

數據重建的流程如圖2所示。

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

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

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

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

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

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

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

2 方法測試

2.1 單道地震數據測試

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

2.2 合成地震記錄測試

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

圖4 單道地震記錄測試：(a)含有隨機缺失的單道數據，從左到右的子波依次為：主頻15 Hz振幅1.0，主頻30 Hz振幅1.0、主頻15 Hz振幅3.0和主頻30 Hz振幅3.0；(b)重建的單道數據；(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算法進行地震數據矩陣分解得到的基函數M，其中圖6a表示BPMF的基函數，圖6b表示離散余弦變換(Discrete Cosine Transform，DCT)的基函數，每一個單獨的數據塊代表基函數的一列。可以看出BPMF算法的基函數的每一個數據塊都非常類似于地震數據的同相軸，區(qū)別在于振幅和相位的不同。而DCT基函數是具有解析數學表達式的傳統(tǒng)數學變換，每個數據塊都是固定的，與輸入數據無關。通過對比可以看出矩陣分解算法可以根據輸入數據自適應的構建基函數M，具有特征提取的功能。

圖5 合成地震記錄測試：(a)原始數據；(b)含隨機缺失數據；(c)Curvelet重建數據；(d)PMF重建數據；(e)BPMF重建數據；(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算法，測試了隨機缺失道數從10道遞增到100道，基函數M的列數k分別為15、20和25，以及規(guī)則化參數λM和λA分別為0.01、0.02和0.03的情況，測試結果如圖7所示，峰值信噪比的計算公式為：

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

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

圖7 PMF和BPMF對 比 測 試：(a)k =15測 試 結 果；(b)k =20測試結果；(c)k=25測試結果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)原始數據；(b)含隨機缺失數據；(c)BPMF重建數據；(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展示了隨機缺失道數從10道遞增到100道時的全部重建結果。通過這一組測試可以看出當數據中存在隨機噪音時，BPMF算法依然可以穩(wěn)定高效的對缺失地震數據進行重建。

2.3 實際地震記錄測試

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

圖10 實際地震數據測試：(a)原始地震數據；(b)含缺失道地震數據；(c)BPMF重建地震數據；(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 實際地震數據測試單道記錄：(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 結論與展望

圖12 實際地震數據基函數Fig. 12 Basis matrix of real seismic data

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

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

[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.