黃亦磊+運乃丹+周仕勇

摘要：為研究噪音對主動源信號走時變化提取精度的影響，采用對比數(shù)值實驗的方法。按照常規(guī)流程，首先用反褶積去除震源激發(fā)時間的不確定性干擾，壓縮子波以得到更加尖銳的震相，由于發(fā)現(xiàn)反褶積會降低信號的信噪比，于是提出了臺差法來代替反褶積矯正氣槍激發(fā)的不確定性，再結合干涉法提取走時變化，并定量對比了兩者的走時變化提取精度，最后嘗試用臺差法提取主動源走時的日變化。結果表明，臺差法可以有效矯正震源激發(fā)時刻的影響，且不降低信噪比，從而也能獲得較高走時變化提取的精度。

關鍵詞：主動源；走時變化；走時臺差法；信噪比

中圖分類號：P315.31文獻標識碼：A文章編號：1000-0666（2017）04-0595-10

0引言

地震波信號是照亮地球內部結構的一盞明燈。復雜的地震波信號是由斷層的破裂過程和介質格林函數(shù)耦合而成。傳統(tǒng)的地震學分為研究震源運動學和地球結構研究（Lay，Wallace，1995）。從地震信號研究震源S（t）或者介質結構G（t）就必須使兩者解耦。因此，大多數(shù)震源或者介質的研究中，都必須假設另一者有先驗的模型。

由于主動源具有精確的激發(fā)時刻、激發(fā)地點，為地下介質的研究提供了良好的可重復震源，從而使精細獲取介質結構變化成為可能。Niu 等（2008）利用壓電陶瓷作為主動源，發(fā)現(xiàn)在實驗場附近2個小地震發(fā)生前，壓電陶瓷產(chǎn)生地震波信號走時有顯著的延遲。在物理上，由于地震發(fā)生之前都會有裂隙和孔隙壓力的變化，這個變化在介質波速上有所體現(xiàn)（Simmons，1964；Yukutake，1988），當這種精細變化被主動源信號所捕獲時，將有助于我們認識地震孕育的物理過程。壓電陶瓷具有頻率高、衰減快、傳播距離十分有限的特點，不能在大區(qū)域的研究中作為廣泛應用的主動源。因此，陳颙院士在其所主導建設的賓川地震信號發(fā)射臺中，引入了4槍組合的氣槍激發(fā)源作為主動源，從而克服了傳播距離有限的問題（羅桂純等，2006；王寶善等，2011；王偉濤等，2009；楊微等，2013）。

目前，賓川發(fā)射臺氣槍實驗已順利進行多年，積累了大量的實驗數(shù)據(jù)。如何提取高精度的主動源信號走時變化是眾多學者關心和研究的問題。筆者首先介紹目前提取走時變化常用的干涉法和處理流程，之后對提取走時變化流程中反褶積這一步驟提出可能存在放大噪音而影響走時變化提取精度的問題，并進行相應的數(shù)值實驗來驗證這一想法，進而提出一種可能替代反褶積的臺差法來矯正震源的影響，并進行了相應的實證。

1干涉法提取走時變化簡介

波速測量關鍵的問題在于走時的獲取，對于相對走時變化在0.01%～1%量級的精細變化（Wang et al，2008；劉自鳳等，2015）探測中，傳統(tǒng)手動拾取震相到時的方法已經(jīng)不再適用。由于主動源氣槍激發(fā)的重復性，波形上也表現(xiàn)出一定的相似性，基于互相關提取相對走時變化的方法比較適用（劉自鳳等，2015；王彬等，2012）。本文以Wang 等（2008）所提出的干涉法為基礎，開展信號走時提取。

1.1互相關函數(shù)

對于兩個能量有限、定義在[t1，t2]上的有限長地震實信號x（t）和y（t），相關函數(shù)定義為：Rxy（τ）=∫t2t1x（t）y（t-τ）dt（1）其中：τ為時間延遲。物理意義上，相關函數(shù)表征了將一個信號相對另一個信號移動之后的相似性，為了方便比較，引入無量綱歸一化的相關系數(shù)ρxy：ρxy（τ）=Rxy（τ）∫t2t1x2（t）∫t2t1y2（t）（2）互相關系數(shù)在[-1，1]間變化，互相關系數(shù)值越大說明2個信號越相像。2個相同的信號在τ=0時，互相關系數(shù)取1。

1.2干涉法

這里的干涉法特指為氣槍設計的一種基于互相關提取走時變化的方法。在氣槍主動源的走時變化探測過程中，存在以下2個問題，必須設計與之相應、適合氣槍主動源提取走時變化的方法：

（1）氣槍信號由很多震相組成，不同震相代表了地震波在地下介質中傳播的不同路徑，因此形成了對地下不同的采樣空間。不同路徑上的走時變化可能不同，所以有必要截取出一定長度的時窗，使之只包含一個特定的震相，以保證窗內所提取的走時變化是穩(wěn)定而一致的。

地震研究40卷第4期黃亦磊等：臺差法提高低信噪比下走時變化探測精度的研究（2）地震儀的采樣率為100 Hz，對應的時間采樣間隙0.01 s，然而介質速度變化引起的走時變化可能小于時間采樣間隙，因此我們有必要對記錄的地震信號進行有效插值，將走時測量精度大幅提高至0.001 s以上。

干涉法的原理如圖1所示。首先我們需要得到一個參考波形，將參考波形和待測信號分成若干個小時窗，每個小時窗可以左右移動一定的步長，將移動過程中最大的互相關系數(shù)（圖1c）對應的時間移動定為該時窗的相對時移（圖1d）。每一個時窗都對應一個時間延遲值，我們需要選定一個波形顯著、互相關系數(shù)較大且平穩(wěn)的時間段，取這段時間中波形的互相關系數(shù)（圖1c）作為走時變化的權重，加權求平均值，作為該道信號的走時變化。

從時移量得到速度變化有很簡單的數(shù)學關系（Schaff，Beroza，2004）。由于信號的傳播距離s（臺源距）可以確定，假設介質速度原來是v，對應的走時為t，速度變化為dv，走時變化為dt，則有：

s=v·t=（v+dv）（t+dt）（3）忽略二階小量dv·dt，則有：dvv=-dtt（4）從干涉法的實現(xiàn)過程可以看出，信號走時變化提取或多或少會受到信號處理流程的影響?，F(xiàn)列舉如下，以提示讀者之后研究的注意：

（1）地震信號的前期處理。包括去均值、去除線性趨勢、濾波和反褶積。其中去均值、去除線性趨勢、濾波都是提高信號信噪比的有效手段。而反褶積則能有效矯正震源發(fā)震時刻不確定性的影響，并壓縮子波得到更加尖銳的震相序列，但同時會降低信噪比。

（2）干涉法中時窗大小選取以及滑動步長選取。時窗大小以及滑動步長的選取帶有經(jīng)驗性和探索性，也依據(jù)不同研究問題的目標而定。endprint

（3）選取哪一段波形對應的走時變化做平均得到最后的走時變化。

如何選取時窗滑動步長、滑動步長選多大以及最后算平均走時變化的時窗，具有經(jīng)驗性和探索性，這是比較主觀的影響因素，并非本文所要研究的問題。在之后的研究中，我們盡可能控制使窗口長度以及滑動步長等主觀因素一致。主要定量研究噪聲及濾波和反褶積等步驟對走時變化提取精度的影響，并嘗試提出一種保持原信號信噪比，但又能一定程度上矯正震源時間不確定性的臺差法，并對這種方法的有效性開展了實證研究。

2干涉法提取走時變化可能影響因素探究

影響一個物理量準確測量的原因有兩種：一是客觀存在的誤差和噪音；二是主觀人為因素。

儀器的誤差無法消除，自然的噪音也客觀存在，因此要在提取走時的過程中，盡量降低噪音，并使用合適的信號處理流程。筆者將重點研究在數(shù)據(jù)處理流程中可能會影響走時變化提取的因素。

劉自鳳等（2015）提出了賓川主動源氣槍定時變化提取的一個流程，如圖2所示。這個流程是目前大多數(shù)走時變化提取所遵循的流程，這個流程看起來是比較合理的，但是劉自鳳等（2015）并沒有做數(shù)值方面的測試，也沒有做主要步驟對提升結果的可靠性的分析。接下來，筆者將從理論測試出發(fā)，研究在不同的噪音水平和處理流程中（主要是濾波和反褶積）對走時變化提取結果的影響。

2.1噪音背景下反褶積對走時變化提取的影響

在與地震信號同頻率的噪音存在的背景下（濾波無法消除噪音），實際地震信號對震源時間函數(shù)反褶積，噪音的影響可能會被放大。為了研究這一問題，需如下數(shù)值實驗。

選擇2013年前80天53281臺（位置如圖3所示）數(shù)據(jù)作為研究對象，將80天數(shù)據(jù)疊加，如圖4a 所示。并把參考臺CKT0（位置如圖3所示）相同時間段的數(shù)據(jù)疊加作為震源時間函數(shù)（圖4b）。對53281臺與CKT0臺的疊加信號加入信噪比為20、15、8、5、-1、-5的噪音高斯白噪音，如圖5所示。

進行以下4組實驗，每組實驗我們都先進行去均值、去除線性趨勢：

（1）對加入噪音之后的信號不做任何處理。

（2）對加入噪音之后的信號做2～8 Hz的濾波，如圖5b所示。

（3）對加入噪音之后的信號直接反褶積，如圖5c所示。

（4）對加入噪音之后的信號線做2～8 Hz的濾波，之后反褶積，再做2～8 Hz濾波，如圖5d所示。

在4組實驗中，把圖5中4幅子圖的最底下一道信號作為組內信號的參考模板，以保證模板和待檢測信號的一致性。將模板和每條信號做干涉（干涉允許最大的左右移動步長為0.2 s），提取每一個時點的走時變化，并把5～15 s的時窗內每一個時點對應的走時變化展示在圖6中。由于這些信號都是原始信號加噪音而來，理論上走時變化應該為0，走時變化偏移0的程度代表了噪音對走時變化提取的影響。對比4組實驗數(shù)據(jù)，總結如下：

（1）組內實驗結果對比，可以發(fā)現(xiàn)噪音越高，對走時變化提取影響越大。符合常規(guī)預期。

（2）濾波可以有效提高信噪比，降低走時變化提取的不確定性。符合常規(guī)預期。

（3）反褶積降低信噪比，反褶積之后提取走時變化的結果都差于不做反褶積，反褶積實質上起到了一個放大噪音的效果。

在此基礎上，我們以53281臺80天疊加信號為基礎，添加噪音，使其信噪比從-10到15連續(xù)變化，對每一個信噪比下的信號濾波之后，做反褶積（嘗試0.000 1，0.001和0.01的反褶積waterlevel水平）和不反褶積（只濾波）處理，之后提取走時變化，仍然以原始數(shù)據(jù)疊加而未添加噪音的信號作為參考模板。取每條信號的信號段（8～12s）提取的走時變化的平均值，作為該信號的走時變化。結果如圖7所示，可以看出，所有反褶積處理之后的信號，走時變化提取結果都比只濾波處理的結果差，而且，隨著反褶積waterlevel水平的增加，反褶積之后提取的走時變化結果也會更加好且穩(wěn)定。

綜上所述，我們在圖7中定量展示了在存在不同水平噪音時，反褶積這一處理流程可能帶來的走時提取誤差。反褶積是一把雙刃劍，它可以消除震源的影響，壓縮子波，但是也會降低信號的信噪比，給走時變化提取帶來不確定性。且高噪音水平要求反褶積時用更高的反褶積水準值。

3替代方案——臺差法及其有效性實證

3.1反褶積的可能替代方法——臺差法

提出用反褶積處理信號的初衷之一，是矯正氣槍源的不穩(wěn)定影響，比如激發(fā)時刻不同步以及水位變化等。反褶積可以去除震源的影響，留下介質的格林函數(shù)。但是，上述的數(shù)值實驗說明，在存在噪音的情況下，反褶積會降低信噪比，從而導致變化提取不準確。另外，由于參考臺得到的信號并非真實的震源時間函數(shù)，可能還需考慮方位帶來的輻射花樣問題，有一定的局限性。

因此，我們提出用臺差法來代替反褶積以消除震源的影響，從而避免由于信噪比低而引起提取走時變化的精度問題。臺差法的想法比較簡單：震源如果存在激發(fā)時刻不準確或水位變動的影響，這種走時變化是一種整體時移，在每一個臺站的記錄都是近乎相同的，因此利用非常靠近震源的參考臺進行自身干涉而提取到的走時變化，可以認為就是源的變化所導致的，然后再用其他臺站提取的走時變化結果減去參考臺提取的走時變化結果，這樣便可矯正源的影響?？偨Y臺差法的步驟如下：（1）對信號去均值，去除線性趨勢，濾波；（2）參考臺用干涉法提取走時變化；（3）待研究臺站用干涉法提取走時變化；（4）待研究臺站走時變化減去參考臺走時變化得到最終的走時變化。

為了驗證臺差法的正確性，我們做了如下實驗：選取2013年前80天的53281臺和CKT0臺數(shù)據(jù)，分2組做對比實驗：第一組用53281臺單槍數(shù)據(jù)對CKT0同一槍數(shù)據(jù)做反褶積（反褶積水準值0.001），然后單天數(shù)據(jù)疊加，提取走時變化（走時變化用8.5～10.5 s時段走時變化平均）；第二組，對CKT0臺和53281臺數(shù)據(jù)單天疊加，然后分別用干涉法提取CKT0臺和53281臺走時變化，然后用臺差法（53281臺走時變化用9～11 s時段走時變化平均，CKT0臺走時變化用0.5～1.5 s時段的走時變化平均）處理數(shù)據(jù)。上述2組實驗在提取走時變化時的模板均為2臺在這80天內的所有數(shù)據(jù)的疊加，結果如圖8所示。可以發(fā)現(xiàn)，CKT0臺和53281臺提取的走時變化有很好的相關性，說明源的誤差確實在這兩個臺都有所體現(xiàn)，相減之后，得到的走時變化非常平緩，在這個階段內，并無大的環(huán)境因素改變，我們有理由認為這段時間的介質速度基本不變。對比反褶積提取的走時變化結果，可以發(fā)現(xiàn)臺差法提取的走時變化結果比反褶積提取的走時變化結果更加穩(wěn)定。endprint

3.2實際資料處理

為了進一步驗證本文提出的臺差法的有效性，本文研究了氣壓以及地震可能對介質速度的影響（Spane，2002）。首先，選取了2015年11月20—28日的集中實驗階段作為研究對象，在此階段，氣槍激發(fā)比較密集，且有氣壓記錄，適合研究氣壓對走時變化的影響。首先選擇53264臺、53274臺、53281臺和53266臺作為研究對象（圖3），

首先用臺差法對單槍數(shù)據(jù)提取走時變化，進行加權平均的走時變化的時間范圍分別為3～5 s、3～5 s、8～10 s和4～6 s，模板為9日內所有數(shù)據(jù)的疊加，4個臺的走時變化和氣壓變化如圖9所示?？梢园l(fā)現(xiàn)，在誤差范圍內，這4個臺的走時變化都與氣壓有非常一致的變化，但是走時變化似乎領先氣壓變化一些相位。為了進一步觀察他們的相關關系，本文細化了實驗，選取53274和53281臺作為研究對象，將這兩個臺在集中實驗期間，每天同一小時、同一分鐘的信號疊加，然后提取走時變化，進行加權平均的走時變化的時間范圍仍為3～5 s和8～10 s，模板同樣為9日內所有數(shù)據(jù)的疊加，得到的將是一天之內的走時變化。從圖10可以看出，這兩個臺的走時變化隨著氣壓都有正向相關的趨勢，但是走時變化都領氣壓變化一些相位。

除此之外還計算了一天之內氣壓對于走時變化的具體數(shù)值影響，從走時變化結果中取出最大值和最小值做差，得到dt；再取平均一天氣壓中最高氣壓和最低氣壓，做差得到dp。則氣壓對走時變化的最大影響約為：Re=dtt×dP（5）其中：t為干涉過程中對走時變化進行加權平均時所選時間段內的地震波平均到時。53274臺和53281臺的結果分別為Re74=5.3×10-7/ Pa和Re81=1.2×10-6/ Pa，這與前人的研究基本一致（Silver et al，2007；Yamamura et al，2003）。

實驗室進行的巖石實驗表明，氣壓的變化會改變巖石內部裂隙中物質的密度，從而改變地震波在巖石中的傳播速度，使波的走時發(fā)生改變（Nur，1971；Walsh，1965），但實際地球內部的巖體結構更加的復雜，氣壓與走時變化之間的具體關系仍待進一步探究。另外，除了氣壓對走時變化有影響之外，固體潮和溫度（Leary et al，1979）等因素也可能通過影響介質速度而使走時變化產(chǎn)生波動，目前我們還沒有進行這方面的探究，這也是接下來的研究內容之一。

4討論

研究地震孕育過程中介質可能發(fā)生的變化，對于我們認識、理解地震發(fā)生過程具有重要意義。如果能有效地提取出介質的變化，并研究清楚這些變化與地震之間的關系，則我們在地震前兆的認識上將更進一步。

反褶積理論上能去除源的干擾。但實際上，由于氣槍水下激發(fā)產(chǎn)生的地震波的機理尚未研究透徹，只能用距離氣槍西面50 m處的一個參考臺來近似源，在現(xiàn)有的條件下是一個不錯的選擇。另在一個問題，即參考臺雖然距離氣槍很近，但是受到震源輻射花樣方位角的影響，因此方位角偏離參考臺比較大的臺站，用參考臺信號反褶積的影響程度仍有待考究。

至于臺差法能否準確地提取出介質的走時變化，有待進一步的驗證。臺差法目前來看有一個缺點，就是不能很好地將參考臺信號對應到待研究的臺站。由于參考臺P波、S波混疊在一起，而在待研究臺站可以明確選擇某一震相，這樣做差可能難以解釋他們的對應關系，但是震源的不確定性在各個震相應該都有體現(xiàn)，所以我們相信可以部分消除源的影響。臺差法的第二個問題在于射線路徑的問題，臺差法可以明確的矯正氣槍的觸發(fā)誤差、水位變化帶來的誤差，但是對于不同的臺站，射線路徑不同，在從震源到參考臺部分的走時，我們沒有考慮不同臺站是不一樣的，這一部分誤差不能保證完全消除，有待進一步研究。但是鑒于CKT0臺離源非常近，我們認為影響不是很大。

進一步，我們將研究走時變化提前于氣壓變化的原因，并將用理論計算得出固體潮的理論值進行對比。另外，下一步將回溯地震，觀察在其發(fā)生前是否有走時變化的突變。

5結論

在目前無法理論推導震源函數(shù)的情況下，前人想到用離震源非常近的參考臺的垂向分量作為震源時間函數(shù)，來進行反褶積，在高信噪比的情況下，是一種可行之舉。但是我們通過實踐發(fā)現(xiàn)在存在較高噪音的情況下反褶積可能導致測量誤差加大，使得走時變化的提取不準確，進而設計了理論實驗和實際數(shù)據(jù)的實驗，發(fā)現(xiàn)反褶積這個過程在走時提取中確實會降低信號的信噪比，從而影響走時變化的精確提取。然后，我們定量研究了噪音水平在走時變化提取中的影響。結果表明，反褶積只有在高信噪比的情況下（8以上）才不會影響走時變化提取的精度。而且原始信號的信噪比越低時，反褶積可能將信號的噪音放大越嚴重。所以當原始信號的信噪比相對較低時，應該提高反褶積的水準值。實驗發(fā)現(xiàn)，不進行反褶積的信號與反褶積之后的信號提取走時變化相比，在信號信噪比低于8時，未反褶積的信號提取的走時變化更加準確和平穩(wěn)，且對噪音的抵抗能力更強。

為了解決這個問題，筆者提出了一種可能的解決方法，也就是用參考臺自身信號的走時變化來矯正遠臺的走時變化，這種方法，從目前的結果來看，也不失為一種矯正源效應的方法。另外，進一步進行了該方法的實證，結果發(fā)現(xiàn)，氣壓變化與氣槍走時變化有很強的相關性，但是氣槍走時變化總是領先氣壓變化，猜想可能和氣溫、固體潮等因素相關，這個負載的物理問題有待進一步研究。

參考文獻：

劉自鳳，蘇有錦，王寶善，等.2015.賓川主動源地震波走時變化分析方法研究[J].地震研究，38（4）：591-597.

羅桂純，王寶善，葛洪魁，等.2006.氣槍震源在地球深部結構探測中的應用研究進展[J].地球物理學進展，21（2）：400-407.

王寶善，楊微，王偉濤，等.2011.賓川氣槍信號發(fā)射臺性能分析[C]//中國地球物理學會年會第二十七屆年會論文集.endprint

王彬，楊潤海，王寶善，等.2012.地震波走時變化精確測量的實驗研究[J].云南大學學報：自然科學版（增刊2）：15-20.

王偉濤，王寶善，葛洪魁，等.2009.利用主動震源檢測汶川地震余震引起的淺層波速變化[J].中國地震，25（3）：223-233.

楊微，王寶善，葛洪魁，等.2013.大容量氣槍震源主動探測技術系統(tǒng)及試驗研究[J].中國地震，29（4）：399-410.

LAY T，WALLACE T C.1995.Modern global seismology[M].Manhattan Academic press.

LEARY P C，MALIN P E，PHINNEY R A，et al.1979.Systematic monitoring of millisecond travel time variations near Palmdale，California[J].Journal of Geophysical Research Atmospheres，84（B2）：659-666.

NIU F，SILVER P G，DALEY T M，et al.2008.Preseismic velocity changes observed from active source monitoring at the Parkfield SAFOD drill site[J].Nature，454（7201）：204-208.

NUR A.1971.Effects of stress on velocity anisotropy in rocks with cracks[J].Journal of Geophysical Research，76：2022-2034.

SCHAFF D P，BEROZA G C.2004.Coseismic and postseismic velocity changes measured by repeating earthquakes[J].Journal of Geophysical Research：Solid Earth，109（B10）：67-85.

SILVER P G，DALEY T M，NIU F，et al.2007.Active Source Monitoring of Cross-Well Seismic Travel Time for Stress-Induced Changes[J].Bulletin of the Seismological Society of America，97（1b）：281-293.

SIMMONS G.1964.Velocity of shear waves in rocks to 10 kilobars，1[J].Journal of Geophysical Research Atmospheres，69（6）：1123-1130.

SPANE F A.2002.Considering barometric pressure in groundwater flow investigations[J].Water Resources Research，38（6）：14-1—14-18.

WALSH J B.1965.The Effect of Cracks on the Compressibility of Rock[J].Journal of Geophysical Research Atmospheres，70（2）：381-389.

WANG B，ZHU P，CHEN Y，et al.2008.Continuous subsurface velocity measurement with coda wave interferometry[J].Journal of Geophysical Research：Solid Earth，113（B12）：36-44.

YAMAMURA K，SANO O，UTADA H，et al.2003.Long-term observation of in situ seismic velocity and attenuation[J].Journal of Geophysical Research Solid Earth，108（B6）：1-5.

YUKUTAKE H.1988.In situ measurements of elastic wave velocity in a mine，and the effects of water and stress on their variation[J].Tectonophysics，149（1-2）：165-175.endprint