基于約束Markov決策過程的初至自動識別技術

羅飛， 王華忠

同濟大學海洋與地球科學學院， 波現象與智能反演成像研究組， 上海 200092

0 引言

在地震勘探層析成像、靜校正、速度分析、AVO分析和地質解釋中地震數據走時信息起著重要的作用，眾多地球物理學家提出了許多算法，用以進行走時的拾取，這些方法都有各自的優缺點以及適用范圍.所以，走時自動拾取算法的研究有其現實意義，算法的穩定性和準確性在經濟上具有重要價值.傳統走時拾取算法大體可分為兩類：滑動時窗法和相干法(Molyneux and Schmitt，1999).在滑動時窗方法中，地震信號序列其屬性在連續或重疊的移動窗口中計算(Coppens，1985；Spagnolini，1991；Sabbione and Velis，2010).相干類方法依賴于使用一些相似度測量技術比較單個或多個波形(López and Aldana，2009).同時，近幾年來，人工智能技術迅猛發展，機器學習算法在地球物理勘探領域的應用越來越廣泛(Jia and Ma，2017；Shi et al.，2020).

走時拾取工作是地震數據處理的重要一步，也是人工智能算法應用較為廣泛的一個領域.比如，傳統的神經網絡(Artifical Neural Network，ANN)算法就已廣泛應用于地震(微震)事件的自動分類識別中(Turhan et al.，1988；Veezhinathan and Wagner，1990；Murat and Rudman，1992；McCormack et al.，1993；Scarpetta et al.，2005；Esposito et al.，2006，2013；Langer et al.，2006；AitLaasri et al.，2013；Vallejos and McKinnon，2013；Maity et al.，2014；Riggelsen and Ohrnberger，2014；Mousavi and Langston，2016).但早期用于地震(微震)信號識別的網絡結構受限于計算機能力一般設計比較淺，其泛化能力比較欠缺(Murat and Rudman，1992).此外，Mousavi和Langston(2016)指出傳統ANN算法常以工程特征(Engineered features)作為網絡的輸入，通過分析不同特征對自動化地震事件識別精度的影響，認為從原始數據中提取工程特征包含過多不確定性.卷積神經網絡(CNN)一般包括卷積、池化和全連通層，是一種能力強大的深度學習算法(LeCun et al.，1995；于子葉等，2018；李薇薇等，2021).CNN使用其多個卷積層直接從圖像或信號中提取不同的特征或屬性，然后通過完全連接的層對其進行分類，將特征提取與分類放在同一網絡結構中，這樣很大程度上降低了提取工程特征對走時拾取精度的影響.Yuan等(2018)將CNN直接應用于地震初至走時中，將二維原始地震數據作為輸入，相比單道輸入，考慮了波形的空間橫向連續性特征.為了進一步證明CNN算法在海量數據支持下，具有很強的分類能力，Loginov等(2019)以5000個訓練樣本訓練包含4個隱藏層的CNN網絡并用其完成了某3D地震數據(450萬道)的初至走時拾取工作，正確率達到了95%.神經網絡類方法，屬于有監督學習，大量的標簽樣本產生不僅耗時，同時也會引入人的先驗認識.無監督學習算法(比如模糊聚類分析、支持向量機等)直接根據特征屬性將地(微)震信號自動分為幾類，不但能夠完成走時拾取工作，還能夠為有監督學習提供標簽樣本(Chen，2000；蔣一然和寧杰遠，2019；許鑫等，2020).Ma等(2019)基于強化學習理論，在能量比譜上自動化全局尋優實現初至走時拾取，但該方法缺乏對獎勵函數和初始狀態選擇的詳細描述，難以適應復雜波形.還有部分研究人員在傳統走時拾取算法(Sabbione and Velis，2010)的基礎上，使用人工智能算法實現窗函數/拾取策略的自適應選取、質量監控等(Duan et al.，2018；Hollander et al.，2018；Mezyk and Malinowsk，2018).總之，人工智能算法的引入，提高了拾取地震走時信息的自動化程度和精度，為地震數據處理提供了很大幫助.

在本文中，我們將初至拾取看作高維特征屬性空間內帶約束的Markov決策過程(Constrained Markov Decision Processes，CMDPs)，在一定全局尋優準則的約束下獲得積累獎勵值最大的路徑，從而自動化的拾取地震數據的初至信息.文中首先介紹了約束Markov決策理論，通過引入折扣因子γ、受空間幾何信息約束的動作和轉移概率，CMDPs算法能夠自動獲取地震數據的初至信息.接著，基于合成地震數據，分析了文中算法的抗噪能力以及常規MDPs和約束MDPs算法對參數的敏感程度，說明CMDPs算法不僅降低了對起始狀態和折扣因子選擇的難度，同時考慮地震數據的空間橫向連續性，自動的回避地震數據中的壞道信息，具有一定QC功能，保證初至拾取更加準確和自動化.最后，實際地震數據的測試結果，證明CMDPs方法能準確地實現地震初至事件的自動拾取，特別是對弱初至信號或淺層相鄰復雜波形的拾取.

1 約束Markov決策理論

1.1 Markov決策過程

馬爾科夫鏈(Markov chains)，指在隨機過程中，下一個狀態僅取決于當前狀態，和之前(歷史)其他狀態無關，即：

P(st+1|st,…,s1)=P(st+1|st)，

(1)

其中，P(st+1|st)為狀態轉移概率，表示從狀態st到st+1的概率.

馬爾科夫決策過程(Markov Decision Processes，MDPs)即為在滿足馬爾科夫鏈前提下，一個智能體(Agent)采取行動(Action)從而改變自己的狀態(State)獲得獎勵(Reward)與環境(Environment)發生交互的循環學習過程.其一般由五元組〈S,A,T,P,r〉描述(Sigaud and Buffet，2013)：

·S：有限狀態集；

·A：控制狀態發生變化的所有可能動作的集合；

·T：隨機演變的時間(空間)序列；

·P：狀態轉移；

·r：狀態轉移時的獎勵函數(瞬時).

圖1 Markov決策過程示意圖(Sigaud and Buffet， 2013)Fig.1 Sketch showing Markov decision processes (Sigaud and Buffet, 2013)

圖1為Markov決策過程的示意圖.在每個時間步t和當前狀態s下，執行動作at后，狀態s有p(s′|s,a)的概率轉移到下一個狀態s′.r為實現該過程應獲得的瞬時獎勵值.

(2)

(3)

表1 值函數迭代法—折扣策略Table 1 Value iteration algorithm-discounted criterion

1.2 初至拾取中的MDPs

如何將初至識別問題，提煉成一個包含〈S,A,T,P,r〉五個元素的Markov決策過程？以圖2中合成2D地震數據為例，狀態空間S和時間序列T定義為

(4)

其中，i=1,…,nt;j=1,…,nx；sij表示時空域地震數據每一點的位置，ti和xj分別為時間和偏移距/道索引，i和j分別為時間和空間采樣點.

基于Markov決策理論，在獎勵值空間(高維特征屬性集合)中根據一定的準則全局尋優獲得一條累積獎勵值最大的路徑，便對應著初至走時的子集.為了降低Markov決策過程的非線性程度，一個合適的獎勵函數的選取是至關重要的.對于獎勵函數的設計，一般以“生成與初至走時相關性高的屬性”為目標，需對相應的地震數據進行分析，構造一個包含多個屬性的特征空間.因此，我們能夠定義瞬時獎勵函數為

(5)

由式(5)可知，獎勵值函數可以由多個屬性組成.fm(sij)表示第m種特征屬性，αm為加權因子，用于調節多個特征屬性之間的權重.對于初至拾取，其能量變化屬性一般可以認為是相關度較高的敏感屬性.如圖2中地震數據，采用STA/LTA算法(Trnkoczy，1999)得到對應的歸一化獎勵函數(能量比)譜(圖2b).接著，需選取合適起始狀態，設計合適的全局尋優準則.Ma等(2019)基于強化學習理論，將遠偏移距的初至到達時位置作為初始狀態，由下至上進行全局尋優，但當遇到地震數據信噪比較低的情況時，自動化的選取初始狀態位置則會出現困難.考慮到地震數據近偏移距的淺層信噪比往往高于遠偏移距，因此，本文采用由上至下的全局尋優策略并結合后續空間結構信息，以此弱化種子點選擇的困難，相應的動作集合函數A和執行某種動作的概率π(a|s)可以寫為

(6)

=(10%;10%;10%;70%).

(7)

式(7)說明，在狀態s執行不同動作a后能夠轉移到狀態s′的概率可能不同，且選取由上往下全局尋優策略，故向下移動的概率最大，為70%.根據上述公式，將初至拾取定義為一個特征屬性空間內的Markov決策尋優過程后，求解相應(3)式，便能夠獲得初至走時的子集(圖2中綠點所示).

圖2 常規MDPs方法拾取初至走時示意圖(a) 合成2D地震數據； (b) 自動拾取結果覆蓋于獎勵值空間上； (c) 自動拾取結果覆蓋于地震數據上.綠點，自動拾取結果.Fig.2 Schematic diagram of picking up first breaks by traditional MDPs method(a) Synthetic seismic data (2D); (b) Auto-picking result overlaid on the reward value function; (c) Auto-picking result overlaid on seismic data. Green dots: auto-picking result.

1.3 約束Markov決策過程

為了進一步的提高Markov理論自動拾取初至信息的精度，處理更加復雜的情況，本文在常規Markov決策理論的基礎上，結合圖像分割技術(Luo et al.，2018)，引入地震數據的空間結構信息dip(s)，式(6)重新寫為

(8)

(9)

其中，當a∈(0,-π)時，狀態發生轉移的概率較大，從而能夠保證從上而下的尋優準則.改進后的式(8)、(9)，使全局尋優過程不僅保留常規MDPs的特點，還能一定程度上受空間結構先驗信息的約束，后續數值試驗的結果也表明了約束Markov理論的優勢.值得注意的是，空間幾何信息的獲取，除了文中提到的圖像分割技術(Luo et al.，2018)，也能采用其他方法技術.

2 數值試驗

2.1 抗噪能力測試

如圖3所示，對比了不同初至拾取算法的抗噪能力.本試驗所采用合成地震數據共36道，縱向701個采樣點，采樣間隔為1 ms.對其添加不同級別的高斯白噪聲(White Gaussian Noise，WGN)，生成不同信噪比的模擬地震數據，所使用信噪比公式可表示為

圖3 不同信噪比下四種方法初至拾取結果示意圖(a) SNR=5 dB； (b) SNR=-2 dB； (c) SNR=-6 dB； (d) SNR=-10 dB.青色： STA/LTA；藍點：MER；紅點：常規MDPs；綠圈：約束MDPs.Fig.3 First-breaks picked by four methods with different SNR(a) SNR=5 dB； (b) SNR=-2 dB； (c) SNR=-6 dB； (d) SNR=-10 dB. Cyan dots: STA/LTA; Blue dots: MER; Red dots: traditional MDPs; Green circles: constrained MDPs.

(10)

其中，s(t)為原始有效信號，n(t)為高斯白噪聲.按公式(10)計算圖3a—d的信噪比依次為5 dB、-2 dB、-6 dB和-10 dB.圖3中，將四種不同的初至拾取算法的結果覆蓋在地震數據上，其中青色圓點為經典的STA/LTA方法(Trnkoczy，1999)；藍色圓點為MER(Modified Energy Ratio)算法(Wong et al.，2009)；紅色圓點為常規MDPs算法；綠色圓圈為本文改進的CMDPs方法.對比不同信噪比下，四種算法的拾取結果可知，傳統基于單道處理的算法，忽視了地震數據所包含的空間特征(橫向連續性)，從而導致難以識別信噪比低的弱初至信號，而常規MDPs和CMDPs算法，均基于Markov理論，在高維空間進行初至拾取，并考慮與距離相關的折扣因子γ，很自然的回避掉異常拾取結果，能夠抓住弱初至信號，使拾取結果更加符合物理，故基于Markov決策理論的兩種算法抗噪能力強于傳統單道拾取算法.同時，對比常規MDPs方法和本文改進的CMDPs算法，由于CMDPs算法引入了先驗結構約束，其結果比常規MDPs算法更加穩健.圖4定量顯示了不同信噪比下四種算法的性能，根據定義的拾取誤差函數：

圖4 不同信噪比下四種方法初至拾取結果誤差曲線圖Fig.4 Error curves of first break picking using four methods with different SNRs

(11)

式中，Error代表算法的拾取誤差，即所有N道的拾取誤差(測量值與真實值之差的絕對值)之和，單位為ms，第i道的真實值true(i)通過CMDPs算法在原始有效信號s(t)上拾取所得.如圖4所示，MER和STA/LTA兩種基于單道拾取的算法其拾取積累誤差明顯高于考慮地震數據橫向連續性的算法(MDPs和CMDPs).CMDPs(綠線)誤差增加最為緩慢，進一步說明在考慮空間結構先驗信息后，CMDPs算法的性能要好于常規MDPs算法.

2.2 折扣因子γ敏感性測試

當獎勵函數和起始狀態確定后，基于Markov決策理論拾取地震數據初至信息，需要考慮折扣因子γ，用以回避由壞道引起的錯誤初至拾取結果(異常值或奇異點).但是，常規MDPs和本文改進的約束MDPs算法對其敏感程度存在差異.如圖5a所示，為理論合成地震數據，共26道，其中第23～26道信噪比較低，初至能量弱于相鄰直達波能量.圖5b—d為考慮不同折扣因子γ情況下常規MDPs和CMDPs的初至拾取結果.當γ設置較大為0.6時，MDPs算法對淺層相鄰復雜波形識別能力較弱，出現拾取錯誤(圖5b中紅點)，這是因為γ越大，遠處的狀態獎勵值影響越大，全局尋優時自然難以忽視強直達波能量.隨著γ值逐漸減小到0.2，MDPs算法能獲得較準確的初至信息(圖5c中紅點，遠偏移距仍存在輕微抖動)；但是當γ值進一步減小為0.1時，根據公式(2)可知，當前狀態只受相鄰很近的狀態獎勵值影響，從而再次出現拾取偏差(圖5d中紅點).而圖5b—d中綠色圓圈為本文改進的約束Markov決策理論初至拾取結果，可見在不同折扣因子下，CMDPs由于空間趨勢的約束，拾取結果基本一致，說明其對γ的取值并不敏感.總體來說，折扣因子γ越小，常規MDPs拾取結果應該更加平緩，即相鄰初至走時變化小.當地震數據的淺層信噪比較高時，常規MDPs和CMDPs的折扣因子γ的選取還是比較容易的，兩者的拾取效果也相當.但是，當出現數值試驗中所展示的復雜波現象后，約束MDPs的優勢就得到了體現.

2.3 OBC數據測試

前面理論數據試驗，分析了約束Markov理論在初至拾取工作中的特點.在本節，選取南海某OBC實際地震數據(圖6)說明CMDPs算法在實際資料初至拾取中的效果.如圖6a所示，該地震數據為氣槍震源激發，道間距25 m，共計185道，最大偏移距3.5 km，時間采樣間隔2 ms，采樣時間為3 s，圖中黃框處顯示鄰近弱初至波的下方存在強能量的折射波.針對該數據，采用STA/LTA算法生成狀態獎勵值空間(圖6b)，選取折扣因子γ為0.3，分別采用常規MDPs算法和約束MDPs算法自上而下進行初至走時的拾取.圖6b為拾取結果覆蓋在獎勵值空間上，圖6c則是拾取結果覆蓋在地震數據上.對比兩種方法的拾取效果，可以發現如果采用常規MDPs方法進行全局尋優自動拾取初至走時，在經過黃框處時容易忽略較弱的初至波，得到錯誤的拾取結果(圖6b、c中紅色圓點).因此，在常規MDPs的基礎上，本文改進的CMDPs算法引入結構信息dip(s)，相當于在狀態s執行動作a后能夠轉移到狀態s′的概率同時受到dip(s)的約束.圖6b、c中綠色圓圈即為約束MDPs算法得到的拾取結果，其在黃框處仍然能拾取準確的弱初至信號，說明CMDPs方法包含了常規MDPs和空間幾何信息的優點，能夠更加精確的自動化拾取初至信息，特別是在處理相鄰較近復雜波現象時，存在明顯優勢.

圖5 常規MDPs與約束MDPs在不同折扣因子γ情況下初至拾取結果對比圖(a) 理論合成地震數據； (b) γ=0.6時初至拾取結果； (c) γ=0.2時初至拾取結果； (d) γ=0.1時初至拾取結果.綠圈，約束MDPs；紅點，常規MDPs.Fig.5 Comparison of first-breaks picking results between traditional MDPs and constrained MDPs under different discount factors(a) Noisy synthetic seismic data; (b) γ=0.6; (c) γ=0.2; (d) γ=0.1. Green circles: constrained MDPs; Red dots: traditional MDPs.

圖6 (a) OBC數據； (b) 自動拾取結果覆蓋于獎勵值空間上； (c) 自動拾取結果覆蓋于地震數據上綠圈，約束MDPs；紅點，常規MDPs.Fig.6 (a) Real data (from the OBC); (b) Auto-picking result overlaid on the reward value function; (c) Auto-picking result overlaid on seismic dataGreen circles: constrained MDPs; Red dots: traditional MDPs.

圖7 (a) 山前帶數據； (b) 自動拾取結果覆蓋于獎勵值空間上(圖7a黃框區域放大顯示)； (c) 自動拾取結果覆蓋于地震數據上(圖7a黃框區域放大顯示)綠圈，約束MDPs； 紅點，常規MDPs.Fig.7 (a) Data from complex foothills; (b) Auto-picking result overlaid on the reward value function enlarged view with the yellow box in Fig.7a; (c) The auto-picking result overlaid on the shot gather enlarged view with the yellow box in Fig.7aGreen circles: the constrained MDPs; Red dots: the traditional MDPs.

2.4 陸上山前帶數據測試

海上地震數據信噪比一般會比陸上復雜地區數據信噪比高一些，同時由于地表高程影響，海上數據的道間時差變化也相對緩慢一點.為了進一步說明本文算法在實際應用中的普適性，如圖7a所示，為西部某山前帶地區實際地震資料.該數據主要用于近地表建模，因此使用炸藥震源，小道間距(2 m)采集，共351道，最大偏移距為830 m，時間采樣點數為501，采樣間隔4 ms.由于山前帶地區地表存在起伏，地震數據道間時差相對平原地區變化更為劇烈，為了更好地適應陸上山前帶復雜地震數據，對于該實際資料，基于2.2節的分析，假如選擇較小的折扣因子，會導致常規MDPs算法難以適應道間時差變化大的情況，因此，這里常規MDPs和CMDPs算法相較于OBC數據(圖6)均選取更大的折扣因子γ(取值為0.5)使其拾取結果允許一定程度的抖動，狀態獎勵值空間仍然采用STA/LTA構建.通過放大圖7a中黃框區域，圖7b、c更加清晰的展示了兩種算法自上而下拾取的效果，其中紅色圓點為常規MDPs的拾取結果，綠色圓圈為CMDPs的拾取結果.相比于圖6b的狀態獎勵值空間，圖7b中狀態獎勵值空間能量變化更加劇烈，遠偏移距上信噪比更低，在采用較大折扣因子后，常規MDPs算法在信噪比低或壞道處(圖7b、c中紅色圓點)出現了少量錯誤拾取，而約束MDPs算法由于引入了空間構造信息的約束，受折扣因子的影響更小，拾取結果更加合理穩健，說明本文改進的CMDPs算法，考慮空間橫向連續性并受先驗空間信息約束，能夠適應道間時差變化大的情形，相比于常規MDPs，CMDPs對參數依賴程度更低，適用范圍更廣，拾取精度也更高.

3 討論

通過第2節的數值試驗分析，考慮空間橫向連續性的Markov決策過程比傳統基于單道初至拾取的算法抗噪能力更好，同時引入先驗空間幾何信息約束后，相較于常規MDPs算法，約束MDPs方法對折扣因子的依賴程度更低，其拾取結果更穩健.海上和陸上實際地震資料初至拾取結果也表明，CMDPs算法比常規MDPs算法適用性更廣，初至拾取精度更高.

在約束Markov決策理論進行初至拾取過程中，折扣因子的選取、狀態獎勵值空間構建以及空間幾何信息的獲取這幾個方面需要注意.折扣因子選取方面，通過第2.2節的數值試驗，說明CMDPs算法比常規MDPs算法對折扣因子的依賴程度低，一般根據地震數據的信噪比以及道間時差變化情形選擇一個適中的值即可.關于狀態獎勵值空間的構建，文中只提及了STA/LTA算法，這是因為能量的變化對初至比較敏感，文中所涉及的地震數據使用STA/LTA方法生成的獎勵值空間已經能夠滿足CMDPs需求.但是，理論上，構建獎勵值函數的方式有很多，并且選取的獎勵函數越合適，CMDPs算法的非線性程度越低，初至拾取的精度也就越高.對于空間幾何信息的提取，文中使用的是圖像分割技術，也能夠依靠其他更先進的構造提取算法.總之，CMDPs算法可以看作一套地震數據初至識別技術的流程框架，與其相關的一些重要組成部分，可以根據處理人員對數據的認識程度進行調整，以期獲得更加精確合理的初至信息.

任何初至拾取算法均有其適用性，CMDPs算法在遇到比文中信噪比更低，道間時差變化更劇烈的實際地震資料時，也需要進行必要的地震數據預處理，比如去噪、地表一致性校正以及靜校正等.同時，本文數值試驗均使用二維數據，從第1節的理論分析可知，CMDPs理論也能直接構建3D的狀態獎勵值空間、狀態轉移動作、轉移概率以及空間結構信息，從而很自然的將約束Markov理論拓展到3D情形，但是，3D情況下的全局尋優難度更高效率更低，因此在實際地震資料處理中，選擇在2D空間還是在3D空間進行地震數據的初至拾取工作，需要結合實際生產需求設計合理方案.

4 結論

致謝感謝中國石油勘探開發研究院及西北分院、中海油研究院和湛江分公司、中國石化物探技術研究院和勝利油田分公司對波現象與智能反演成像研究組(WPI)研究工作的資助與支持.感謝審稿專家提出的修改意見.