一種聯合SBL和DTW的疊前道集剩余時差校正方法

石戰戰 ，夏艷晴 ，周懷來 ，王元君

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

0 引言

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

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

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

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

1 方法原理

1.1 稀疏貝葉斯學習

基于MMV模型的信號JSR的基本數學模型為

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

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

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

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

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

當γi=0時，對應的Xi（X的第i行）為零，算法迭代過程中，大部分的γi都會變成0，實現了稀疏約束作用，Bi為協方差矩陣，描述Xi時序結構。因此，參數γi和Bi的估計是SBL算法的核心，可以通過二型最大似然法（Type-ⅡMaximum Likelihood）進行估計。

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

可以看出，式中的X非零行正好對應著中x的非零塊，可以采用BSBL算法進行優化。

TMSBL算法和詳細推導過程見文獻［12］，這里僅介紹基本思路。

考慮到式（2）中X每行對應不同的Bi會引起算法過擬合，后續公式推導均采用統一的協方差矩陣B，將公式（1）轉化為SMV模型（3）后，向量x的先驗概率可以寫為

假設地震噪聲為隨機噪聲，滿足Gaussian分布，條件概率密度為Gaussian分布

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

通過Bayesian方法求取后驗概率密度也為Gaussian分布

當參數γ和B確定后，系數矩陣的最優解可以通過最大后驗概率法（Maximum A Posterior，MAP）進行估計。

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

1.2 平滑動態時間規整算法原理

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

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

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

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

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

1.3 疊前道集剩余時差校正方法原理

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

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

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

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

（2）采用SBL對地震道集進行稀疏表示，構造出系數矩陣，相當于地下介質的單位沖擊響應，消除了子波影響。

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

（4）用校正后系數矩陣重構地震記錄。

2 數值模擬

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

表1模型參數Table 1 Model parameters

圖2 地質模型正演結果Fig.2 Forward modeling of geological models

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

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

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

圖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條件下能夠得到同樣的結論。同時，噪聲強度增加，信噪比降低，2種算法處理效果都會變差，DTW處理結果產生明顯畸變（圖5，圖6中箭頭所示）。強噪聲條件下，DTW處理結果噪聲干擾嚴重，同相軸被噪聲淹沒，識別困難，而新算法具有明顯的去噪能力，強噪聲條件下仍然能夠得到高信噪比處理結果，處理后剩余時差得到校正，且同相軸畸變不明顯。由此可知，結合SBL和DTW的剩余時差校正方法對不同品質的地震資料都是有效的。

為了說明聯合SBL和DTW的剩余時差校正方法的原理，進一步驗證算法的有效性，以模型1第8道為例，對比分析強噪聲（SNR=1 dB）條件下2種算法處理結果（圖7）。取2個強噪聲污染地震信號，以第1道［圖7（a）］作為參考道，對第7道［圖7（b）］用新算法和DTW算法分別進行剩余時差校正。單獨使用DTW算法對原始地震道進行剩余時差校正處理后仍然存在明顯的波形畸變［圖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 強噪聲（SNR=1 d B）條件下新算法與DTW算法處理效果對比分析Fig.6 Comparativeanalysisof processing resultsof two residual moveout correction algorithmsunder trong noisecontaminated（SNR=1 dB）

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

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

3 應用實例

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

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

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

4 結論

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

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

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