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

侯思安，張峰，李向陽

中國石油大學(xué)(北京)油氣資源與工程國家重點(diǎn)實(shí)驗(yàn)室，北京 102249

0 引言

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

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

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

1 方法原理

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

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

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

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

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

因?yàn)殡S機(jī)噪音E也是滿足Gaussian分布的，所以：

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

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

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

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

因?yàn)閄是已知的，因此概率密度函數(shù)P(X)是一個(gè)常數(shù)，所以有：

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

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

根據(jù)公式(8)，在采集數(shù)據(jù)X已知時(shí)，求取基函數(shù)M和系數(shù)A的最大后驗(yàn)概率等價(jià)于求解如下的最優(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和隨機(jī)噪音E的方差；表示l2模規(guī)則化；?表示采集數(shù)據(jù)的觀測(cè)系統(tǒng)。

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

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

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

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

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

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

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

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

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

其中，

同理：

其中，

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

因?yàn)榇藭r(shí)采集數(shù)據(jù)X，系數(shù)矩陣A和超參數(shù)都是已知的，所以條件概率密度函數(shù)P(X| A,ξM,ΛM)是一個(gè)常數(shù)，那么：

又因?yàn)榛瘮?shù)矩陣M和系數(shù)矩陣A是無關(guān)的：

所以：

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

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

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

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

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

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

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

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

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

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

2 方法測(cè)試

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

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

2.2 合成地震記錄測(cè)試

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

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

圖5 合成地震記錄測(cè)試：(a)原始數(shù)據(jù)；(b)含隨機(jī)缺失數(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

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

由于地震數(shù)據(jù)的復(fù)雜性以及信號(hào)噪音和缺失的干擾，這些重建參數(shù)是難以直接獲取的，也沒有一個(gè)明確的選取準(zhǔn)則。但是通過這一組數(shù)據(jù)對(duì)比，當(dāng)k從15逐漸增加到25并且規(guī)則化參數(shù)λM和λA從0.03逐漸減少到0.01時(shí)，PMF算法數(shù)據(jù)重建的精度在不斷提高，但是也出現(xiàn)了較大幅度的波動(dòng)，這說明算法的穩(wěn)定性也有所降低，但是BPMF算法始終保持了一個(gè)較高的計(jì)算精度和穩(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對(duì) 比 測(cè) 試：(a)k =15測(cè) 試 結(jié) 果；(b)k =20測(cè)試結(jié)果；(c)k=25測(cè)試結(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 含噪音合成地震記錄測(cè)試：(a)原始數(shù)據(jù)；(b)含隨機(jī)缺失數(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 含噪音合成地震記錄全部測(cè)試Fig. 9 Overview of synthetic noisy data reconstruction test

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

2.3 實(shí)際地震記錄測(cè)試

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

圖10 實(shí)際地震數(shù)據(jù)測(cè)試：(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í)際地震數(shù)據(jù)測(cè)試單道記錄：(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í)際地震數(shù)據(jù)基函數(shù)Fig. 12 Basis matrix of real seismic data

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

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

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