石戰(zhàn)戰(zhàn) ，夏艷晴 ，周懷來 ，王元君

（1.成都理工大學(xué)工程技術(shù)學(xué)院，四川樂山614000；2.成都理工大學(xué)地球物理學(xué)院，成都610059）

0 引言

地下空間的高分辨率成像是地震采集、處理的最終目標(biāo)［1-2］。由于地表?xiàng)l件和地下介質(zhì)的復(fù)雜性及速度分析、靜校正、動校正等處理方法的限制，處理后的疊前地震道集仍然存在剩余時(shí)差，表現(xiàn)為同相軸校正不平，這種資料難以滿足疊前偏移、疊前反演等處理、解釋的需求，會影響地震資料的成像精度。曲壽利［3］指出動校正方法本身的局限性、速度分析誤差和薄層調(diào)諧均是引起剩余時(shí)差的主要原因，并提出了一種基于統(tǒng)計(jì)方法的剩余時(shí)差校正方法。周鵬等［4］對國內(nèi)外相關(guān)方法技術(shù)進(jìn)行總結(jié)，將剩余時(shí)差校正方法歸納為精細(xì)速度分析、地震道集互相關(guān)技術(shù)、相位匹配技術(shù)和統(tǒng)計(jì)法等4類。石戰(zhàn)戰(zhàn)等［5］提出一種基于動態(tài)時(shí)間規(guī)整（Dynamic Time Warping，DTW）和形狀上下文（Shape Context，SC）的剩余時(shí)差校正方法，通過SC取代DTW算法中的距離算子，但該算法采用逐點(diǎn)搬家法進(jìn)行剩余時(shí)差校正，直接使用容易造成波形畸變，同時(shí)DTW算法對噪聲比較敏感，影響規(guī)整路徑計(jì)算的精確性。

由地震子波構(gòu)造過完備字典，對地震道集進(jìn)行稀疏表示，其結(jié)果（系數(shù)矩陣）相當(dāng)于地下介質(zhì)的單位沖激響應(yīng)，類似于地層反射系數(shù)，脈沖對應(yīng)于地層界面，消除了子波影響，而地層信息沒有損失，再對系數(shù)矩陣作時(shí)差校正，就能避免波形畸變。稀疏表示的核心思想是高維空間中的信號可以由低維信號模型來描述［6］。近年來受壓縮感知及相關(guān)數(shù)學(xué)理論發(fā)展的推動，產(chǎn)生了一系列稀疏表示方法，如匹配追蹤［7-8］、基追蹤（Basis Pursuit，BP）［9］和稀疏貝葉斯學(xué)習(xí)（Sparse Bayesian Learning，SBL）等，其中SBL是常用的稀疏表示方法，尤其適用于高相干字典條件下信號的稀疏表示［10］，該方法最早是由Tipping［6］提出并證明是最有效的回歸和分類方法。Wipf等［10-11］將SBL用于一維信號稀疏恢復(fù)，隨后將其拓展為基于多測量向量模型（Multiple Measurement Vectors，MMV）的多變量稀疏貝葉斯學(xué)習(xí)（Sparse Bayesian Learning for MMV，MSBL）算法。與BP算法相比，SBL和MSBL的優(yōu)勢在于其全局最優(yōu)解同時(shí)也是最稀疏解，并且具有較少的局部最小值，在實(shí)際信號處理中，稀疏解常具有一定的結(jié)構(gòu)，SBL框架則能夠充分利用這些結(jié)構(gòu)，提高稀疏表示的性能。Zhang等［12-13］針對塊稀疏信號的塊內(nèi)相關(guān)性，提出了塊稀疏貝葉斯學(xué)習(xí)（Block Sparse Bayesian Learning，BSBL），隨后通過將 MMV 模型轉(zhuǎn)化為單測量向量模型（Single Measurement Vectors，SMV）問題，提出了基于 TSBL（Time Varying Sparse Bayesian Learning）算法和 TMSBL（Sparse Bayesian Learning for Multiple Measurement Vector Model Exploiting Temporal Correlation）算法的聯(lián)合稀疏表示（Joint Sparse Representation，JSR）方法。TSBL算法和TMSBL算法能夠有效地處理時(shí)變稀疏（Time-Varying Sparsity）模型［14-15］。

DTW是一種基于動態(tài)規(guī)劃的模板匹配算法，具有簡單高效的特點(diǎn)，用于對系數(shù)矩陣進(jìn)行剩余時(shí)差校正。由于原始的DTW算法計(jì)算效率較低，出現(xiàn)了多種改進(jìn)算法，眾多學(xué)者［16-20］通過約束或抽象等方法降低計(jì)算量，但DTW算法存在固有缺陷，認(rèn)為模板匹配的2個(gè)信號之間的差異僅由時(shí)間畸變引起，估計(jì)出的規(guī)整路徑不夠精確，Compton等［21］提出平滑動態(tài)時(shí)間規(guī)整（Smooth Dynamic Time Warping，SDTW），能夠估計(jì)出更加平滑和精確的規(guī)整路徑，并將其應(yīng)用于多波匹配。Cui［22］將SDTW用于自動井震標(biāo)定。

單純使用DTW算法難以取得良好的校正效果，將SBL和DTW方法結(jié)合，提出一種聯(lián)合TMSBL和SDTW的疊前道集剩余時(shí)差校正方法，該方法從2個(gè)方面進(jìn)行改進(jìn)：①利用TMSBL對疊前道集進(jìn)行稀疏表示，對計(jì)算出的系數(shù)矩陣進(jìn)行SDTW時(shí)差校正，再利用校正后的系數(shù)矩陣重構(gòu)地震數(shù)據(jù)；②將原算法中的SC-DTW改為SDTW，以期避免波形畸變，同時(shí)實(shí)現(xiàn)高保真剩余時(shí)差校正和隨機(jī)噪聲壓制。

1 方法原理

1.1 稀疏貝葉斯學(xué)習(xí)

基于MMV模型的信號JSR的基本數(shù)學(xué)模型為

式中：S為測量矩陣；W為字典矩陣；X為源矩陣；N為噪聲矩陣。

設(shè)X具有稀疏性，JSR的目的是通過優(yōu)化算法恢復(fù)出系數(shù)矩陣X。受壓縮感知和相關(guān)數(shù)學(xué)理論研究推動，產(chǎn)生了多種JSR算法，其中TMSBL是常用的JSR方法之一，該方法適合于地震信號稀疏表示，原因是：①動、靜校正后，疊前地震道集中各道反映了同一地質(zhì)構(gòu)造的信息，如CMP道集反映了地下同一點(diǎn)的地質(zhì)信息，因此，系數(shù)矩陣具有聯(lián)合稀疏性，并且有一定的時(shí)間、空間結(jié)構(gòu)，TMSBL算法可以較好地利用這些信息，提高稀疏表示的精度；②處理后地震數(shù)據(jù)常存在剩余時(shí)差，系數(shù)矩陣的支撐是變

化，TMSBL可以較好對這種時(shí)變稀疏性進(jìn)行建模；③由地震子波構(gòu)造的字典矩陣，相干度較高，對于這類高相干性字典，與其他稀疏表示算法相比，SBL算法具有較好的適應(yīng)性。

SBL理論框架下，系數(shù)矩陣X的每一行都可以用參數(shù)化Gaussian分布建模

式中：p（·）概率密度函數(shù)；N（·）表示 Gaussian 分布；γi為非負(fù)超參數(shù)，控制系數(shù)矩陣X的行稀疏性。

當(dāng)γi=0時(shí)，對應(yīng)的Xi（X的第i行）為零，算法迭代過程中，大部分的γi都會變成0，實(shí)現(xiàn)了稀疏約束作用，Bi為協(xié)方差矩陣，描述Xi時(shí)序結(jié)構(gòu)。因此，參數(shù)γi和Bi的估計(jì)是SBL算法的核心，可以通過二型最大似然法（Type-ⅡMaximum Likelihood）進(jìn)行估計(jì)。

直接優(yōu)化求解式（2）較為困難，實(shí)際采用的優(yōu)化思路是，令 s=vec（ST）∈ RNL×1，D=W ? IL，x=vec（XT）∈ RML×1，n=vec（NT）∈ RNL×1，（? 為 Kronecker算子；vec（·）表示將矩陣按行展開為列向量；ST表示對矩陣S轉(zhuǎn)置），將MMV模型（1）轉(zhuǎn)化為塊SMV模型

可以看出，式中的X非零行正好對應(yīng)著中x的非零塊，可以采用BSBL算法進(jìn)行優(yōu)化。

TMSBL算法和詳細(xì)推導(dǎo)過程見文獻(xiàn)［12］，這里僅介紹基本思路。

考慮到式（2）中X每行對應(yīng)不同的Bi會引起算法過擬合，后續(xù)公式推導(dǎo)均采用統(tǒng)一的協(xié)方差矩陣B，將公式（1）轉(zhuǎn)化為SMV模型（3）后，向量x的先驗(yàn)概率可以寫為

假設(shè)地震噪聲為隨機(jī)噪聲，滿足Gaussian分布，條件概率密度為Gaussian分布

式中：參數(shù)λ表示噪聲方差；I為單位矩陣。

通過Bayesian方法求取后驗(yàn)概率密度也為Gaussian分布

當(dāng)參數(shù)γ和B確定后，系數(shù)矩陣的最優(yōu)解可以通過最大后驗(yàn)概率法（Maximum A Posterior，MAP）進(jìn)行估計(jì)。

地震道集中存在剩余時(shí)差，系數(shù)矩陣的支撐是變化的，須要采用時(shí)變稀疏模型進(jìn)行建模。時(shí)變稀疏模型可看作是多個(gè)MMV模型串聯(lián)結(jié)果［20-21］，因此，TMSBL可以便利地處理時(shí)變稀疏信號JSR［14］。

1.2 平滑動態(tài)時(shí)間規(guī)整算法原理

式中：ei，j為誤差矩陣E的第i行，第j列元素。

第2步，計(jì)算誤差累積矩陣C∈RM×N

第3步，通過C，計(jì)算出累積誤差最小的規(guī)整路徑u，并同時(shí)滿足3個(gè)約束條件：①邊界條件：u1=（1，1）和 uL=（M，N），要求規(guī)整路徑分別以 c1，1和cM，N為起點(diǎn)和終點(diǎn)；②單調(diào)條件：規(guī)整路徑的縱橫坐標(biāo)均是單調(diào)不減的，即當(dāng)相鄰2個(gè)點(diǎn)為ui=（a，b）和ui-1=（a'，b'），則滿足 a-a'≥0 和 b-b'≥0；③連續(xù)性條件：規(guī)整路徑中相鄰2點(diǎn)必須是C中的相鄰點(diǎn)格，當(dāng)相鄰 2 個(gè)點(diǎn)為 ui=（a，b）和 ui-1=（a'，b'），則滿足a-a'≤1和b-b'≤1，即保證f和g的每一個(gè)坐標(biāo)都會在u中出現(xiàn)。

第4步，計(jì)算流程第2步中所構(gòu)建累積誤差矩陣相當(dāng)于逐點(diǎn)搜索3種可能斜率（-1，0和1）的局部規(guī)整路徑，但估計(jì)出的規(guī)整路徑不夠精確和平滑，與實(shí)際地震資料不符。SDTW將該步驟改為間隔h點(diǎn)搜索2h+1種可能斜率（斜率范圍為-1～1）的局部路徑，計(jì)算出一個(gè)粗網(wǎng)格的累積誤差矩陣，這2種方法構(gòu)建累積誤差矩陣對比如圖1所示，其中圖1（a）為傳統(tǒng)方法搜索3種可能的局部路徑，圖1（b）以h=2為例，說明SDTW搜索2h+1種可能斜率的局部路徑，所構(gòu)造出的規(guī)整路徑精確度必然提高。

圖1 累積誤差矩陣構(gòu)建方法對比Fig.1 Comparison of construction methodsof cumulativeerror matrix

1.3 疊前道集剩余時(shí)差校正方法原理

疊前道集剩余時(shí)差校正方法假設(shè)：①地下地層為層狀介質(zhì)，各層介質(zhì)均為均勻、各向同性彈性介質(zhì)；②地震正演滿足褶積模型（Robison模型），表示為

式中：s為地震記錄；W為子波核矩陣；x為反射系數(shù)序列；n為噪聲。③疊前道集中各道信號反映相同的地下地層信息，因而具有相干性，從信號與系統(tǒng)的角度看，地震地質(zhì)模型等價(jià)于一個(gè)線性時(shí)不變系統(tǒng)，輸入為子波，輸出為地震記錄，反射系數(shù)表達(dá)了層狀地質(zhì)模型。鑒于褶積模型的精度和采集誤差，將x看作地層系統(tǒng)的單位沖擊響應(yīng)更為合適。

褶積模型［式（10）］和稀疏表示數(shù)學(xué)模型［式（1）］具有相同的數(shù)學(xué)形式，其計(jì)算流程如下：

（1）利用地震子波構(gòu)造過完備字典。隨著油氣勘探和開發(fā)工作的不斷深入發(fā)展，薄層儲層預(yù)測已經(jīng)成為一個(gè)研究熱點(diǎn)和難點(diǎn)問題，由于地層厚度薄，薄層的頂、底面?；ハ喔缮姣B加，從式（10）中難以精確反演出地下地層的單位沖擊響應(yīng)。針對這一問題，采用Zhang等［23］提出的過完備子波字典構(gòu)造方法，薄層頂、底面的單位沖激響應(yīng)構(gòu)成一個(gè)薄層反射，通過雙極子分解可以將其分解為一個(gè)偶序列和一個(gè)奇序列之和，其中奇序列具有更高的分辨率，通過將不同深度、不同厚度的薄層反射分解為一系列奇分量序列和偶分量序列，再與地震子波褶積就能構(gòu)造出過完備字典。

（2）采用SBL對地震道集進(jìn)行稀疏表示，構(gòu)造出系數(shù)矩陣，相當(dāng)于地下介質(zhì)的單位沖擊響應(yīng)，消除了子波影響。

（3）利用DTW對系數(shù)矩陣進(jìn)行剩余時(shí)差校正，就能避免直接使用DTW對地震道進(jìn)行處理造成的波形畸變。

（4）用校正后系數(shù)矩陣重構(gòu)地震記錄。

2 數(shù)值模擬

為了驗(yàn)證聯(lián)合SBL和DTW的疊前道集剩余時(shí)差校正技術(shù)的有效性，建立3組不同層厚的4層地質(zhì)模型，模型參數(shù)如表1所列。3組模型均采用30 Hz雷克子波進(jìn)行地震正演，無剩余時(shí)差地震記錄，由于剩余時(shí)差常由速度誤差引起，通過對模型地層速度加入隨機(jī)誤差來模擬剩余時(shí)差，剩余時(shí)差不僅引起同相軸時(shí)間深度差異，還會造成波形畸變（圖 2）。

表1模型參數(shù)Table 1 Model parameters

圖2 地質(zhì)模型正演結(jié)果Fig.2 Forward modeling of geological models

圖3 SBL與BP稀疏表示及重建精度對比Fig.3 Comparison of sparserepresentation and reconstruction accuracy of SBL and BP

對模型3第1道加入不同強(qiáng)度的隨機(jī)噪聲，對比SBL和BP這2種算法稀疏表示精度和重建精度（圖3），圖3中從上向下，噪聲信號的信噪比（Signal to Noise Ratio，SNR）分別為10 dB，5 dB，1 dB，結(jié)果顯示SBL計(jì)算出的系數(shù)向量具有較好的脈沖性，與地層單位沖擊響應(yīng)特征吻合［圖 3（a）—（c）］。SBL重建地震波形曲線光滑，消除了噪聲干擾［圖3（d）—（f）］，強(qiáng)噪聲條件下（SNR=1 dB），仍然能夠達(dá)到較高的重建精度［圖3（f）］。BP算法對噪聲敏感，計(jì)算出的系數(shù)向量中存在較多毛刺，與實(shí)際數(shù)據(jù)誤差較大［圖 3（g）—（i）］。BP 重建地震記錄與SBL算法相比誤差較大［圖 3（k）—（m）］。隨著噪聲變強(qiáng)，SNR逐漸降低，2種算法的稀疏表示精度和重建精度均有所降低，但SBL算法2種相對誤差均小于BP算法，在強(qiáng)噪聲條件下，BP算法2種相對誤差分別達(dá)到了19.057%和47.941%，而SBL算法僅為于5.982%和7.940%，遠(yuǎn)低于BP算法，可見SBL算法具有明顯的計(jì)算精度，用于地震數(shù)據(jù)稀疏表示是有效的。

聯(lián)合SBL和DTW算法（新方法）的優(yōu)勢是同時(shí)實(shí)現(xiàn)了高精度剩余時(shí)差校正和隨機(jī)噪聲壓制。在SNR=1 dB條件下，將新方法與DTW算法［3］的處理效果進(jìn)行對比（圖4），DTW算法對模型2應(yīng)用效果良好［圖4（e）］，而模型1由于地層較薄，速度誤差會使波形干涉疊加，DTW不能使疊加波形解調(diào)［圖4（d）中箭頭所示］，模型3中，上下2個(gè)同相軸存在互相干涉，DTW剩余時(shí)差校正算法使同相軸產(chǎn)生畸變［圖4（f）中箭頭所示］，而聯(lián)合SBL和DTW的新算法能夠使互相干涉的同相軸解調(diào)，剩余時(shí)差校正效果良好，波形畸變遠(yuǎn)小于DTW［圖4（g）—（i）］。由此可知，DTW適用于厚地層剩余時(shí)差校正，薄地層應(yīng)用效果較差，并且容易產(chǎn)生波形畸變，新方法實(shí)現(xiàn)剩余時(shí)差校正的同時(shí)壓制了地震噪聲。

圖4 低噪聲（SNR=10 d B）條件下新算法與DTW算法處理效果對比分析Fig.4 Comparative analysis of processing results of two residual moveout correction algorithms under high SNR（SNR=10 d B）

在SNR=5 dB和SNR=1 dB條件下能夠得到同樣的結(jié)論。同時(shí)，噪聲強(qiáng)度增加，信噪比降低，2種算法處理效果都會變差，DTW處理結(jié)果產(chǎn)生明顯畸變（圖5，圖6中箭頭所示）。強(qiáng)噪聲條件下，DTW處理結(jié)果噪聲干擾嚴(yán)重，同相軸被噪聲淹沒，識別困難，而新算法具有明顯的去噪能力，強(qiáng)噪聲條件下仍然能夠得到高信噪比處理結(jié)果，處理后剩余時(shí)差得到校正，且同相軸畸變不明顯。由此可知，結(jié)合SBL和DTW的剩余時(shí)差校正方法對不同品質(zhì)的地震資料都是有效的。

為了說明聯(lián)合SBL和DTW的剩余時(shí)差校正方法的原理，進(jìn)一步驗(yàn)證算法的有效性，以模型1第8道為例，對比分析強(qiáng)噪聲（SNR=1 dB）條件下2種算法處理結(jié)果（圖7）。取2個(gè)強(qiáng)噪聲污染地震信號，以第1道［圖7（a）］作為參考道，對第7道［圖7（b）］用新算法和DTW算法分別進(jìn)行剩余時(shí)差校正。單獨(dú)使用DTW算法對原始地震道進(jìn)行剩余時(shí)差校正處理后仍然存在明顯的波形畸變［圖7（d）中箭頭所示］。

圖5 中等噪聲（SNR=5 d B）條件下新算法與DTW算法處理效果對比分析Fig.5 Comparative analysis of processing results of two residual moveout correction algorithms under moderate noise contaminated（SNR=5 d B）

圖6 強(qiáng)噪聲（SNR=1 d B）條件下新算法與DTW算法處理效果對比分析Fig.6 Comparativeanalysisof processing resultsof two residual moveout correction algorithmsunder trong noisecontaminated（SNR=1 dB）

圖7 模型1第8道新算法與DTW算法處理結(jié)果對比分析Fig.7 Comparative analysis of processing results of two algorithms in the 8th trace of model 1

采用新方法，先對地震道進(jìn)行稀疏表示，由于存在剩余時(shí)差，同一地層的單位沖擊響應(yīng)存在明顯的時(shí)差［圖 7（e）—（f）中箭頭所示］，再對稀疏表示結(jié)果［圖 7（f）］作 DTW 剩余時(shí)差校正［圖 7（g）］，地層單位沖擊響應(yīng)校正到正確位置，由處理后的系數(shù)向量重構(gòu)地震道［圖7（c）］，剩余時(shí)差得到校正，地震波形沒有發(fā)生明顯的畸變。

3 應(yīng)用實(shí)例

為了驗(yàn)證聯(lián)合SBL和DTW的疊前道集剩余時(shí)差校正方法的有效性，以加拿大Nova Scotia省的Penobscot工區(qū)（數(shù)據(jù)由dGBEarth Sciences提供）為例進(jìn)行試算。研究區(qū)位于Scotian盆地Abenaki凹陷與Sable凹陷的過渡帶，主力儲層為侏羅系上統(tǒng)Abenaki組。工區(qū)面積88.6 km2，包括600條主測線和481條聯(lián)絡(luò)測線，地震面元為12.5 m×25.0 m，60次覆蓋，4 ms采樣，記錄長度為6 000 ms。

原始共中心點(diǎn)角道集，資料噪聲干擾嚴(yán)重，同相軸連續(xù)性差，受資料品質(zhì)影響，前期地震數(shù)據(jù)處理不能拉平同相軸，處理后數(shù)據(jù)存在明顯剩余時(shí)差，如圖 8（a），圖 9（a）中虛線標(biāo)注的同相軸所示，其中圖9為圖8在1.5～1.7 s部分的放大圖。經(jīng)DTW處理后，剩余時(shí)差得到部分校正，但同相軸仍未校平，受DTW算法制約，校正后角道集中存在明顯的波形畸變［圖 8（b），圖 9（b）中箭頭所示］。結(jié)合 SBL和DTW的新算法處理后剩余時(shí)差得到校正，同相軸拉平，同相軸連續(xù)性得到增強(qiáng)［圖 8（c），圖 9（c）中虛線所示］，同時(shí)噪聲得到壓制，提高了地震資料信噪比，利用這一方法處理后的地震數(shù)據(jù)進(jìn)行疊前反演、AVO分析等處理流程，其結(jié)果的精度和可信度也得以提高。

圖8 Penobscot工區(qū)新算法與DTW算法處理效果對比分析Fig.8 Comparative analysis of processing results of two residual moveout correction algorithms in Penobscot project

4 結(jié)論

（1）受地震資料品質(zhì)和處理方法的制約，處理后地震資料常存在剩余時(shí)差，提出一種聯(lián)合SBL和DTW的疊前道集剩余時(shí)差校正方法，能夠同時(shí)實(shí)現(xiàn)噪聲壓制和剩余時(shí)差校正，處理后疊前道集同相軸被拉平，連續(xù)性得到增強(qiáng)，滿足疊前反演、AVO分析等流程對資料品質(zhì)要求。

（2）DTW算法采用逐點(diǎn)搬家的方法進(jìn)行時(shí)差校正，容易引起波形畸變，在DTW時(shí)差校正前，利用SBL對地震數(shù)據(jù)進(jìn)行稀疏表示，其結(jié)果相當(dāng)于地下介質(zhì)的單位沖擊響應(yīng)，能夠消除地震子波影響，再對計(jì)算出的系數(shù)矩陣進(jìn)行剩余時(shí)差校正和重建計(jì)算，就能夠避免波形發(fā)生畸變。

（3）SBL算法能夠適用于高相干度字典，適用于地震數(shù)據(jù)稀疏表示，SBL算法的計(jì)算效率是限制其應(yīng)用于大規(guī)模地震數(shù)據(jù)剩余時(shí)差校正的瓶頸，可采用壓縮感知對地震數(shù)據(jù)進(jìn)行隨機(jī)采樣，降低處理方法的計(jì)算量。