賓川主動源地震波走時變化分析方法研究*

劉自鳳，蘇有錦，王寶善，王　彬，楊　軍，李孝賓

(1.云南省地震局，云南 昆明 650224；2.中國地震局地球物理研究所，北京 100081)

賓川主動源地震波走時變化分析方法研究*

劉自鳳1，蘇有錦1，王寶善2，王彬1，楊軍1，李孝賓1

(1.云南省地震局，云南 昆明 650224；2.中國地震局地球物理研究所，北京 100081)

摘要：挑選了2013年以來信噪比較高、數(shù)據(jù)較為完整的6個賓川主動源接收臺站的氣槍記錄數(shù)據(jù)，采用互相關時延檢測、近場校正等技術進行了高精度波速變化測量，得到以下初步結果：(1)接收臺站下方介質存在δv/v≈10-5～10-2的相對波速變化；(2)2014年9月初至10月底，初至波走時變化較平穩(wěn)，集中在-0.01～0.01 s之間變化，且6個臺的變化趨勢一致。

關鍵詞：主動源；互相關時延檢測；氣槍信號；走時變化

0引言

地震波在地球介質中傳播，帶來了豐富的地下介質物性信息，為我們了解地球內(nèi)部結構及運動變化提供了可能。地震波的很多性質(如衰減、各向異性等)也都可以用來研究地下介質應力狀態(tài)的變化，但地震波速變化仍然是測量精度最高的研究方法。因此，地球內(nèi)部地震波速的差異是確定地下地層結構和橫向不均勻性的重要物理參數(shù)，地震波速的測量一直是地震學中一個令人注目的領域(徐薈，2013)。

波速測量是一個很經(jīng)典的物理問題，如光速的測量、聲速的測量等。在實驗室內(nèi)，測量巖石彈性波速度的方法基本上有3種：共振法、脈沖法(pulse transmission)和超聲干涉法(ultrasonic interferometry)(王彬等，2012)。共振法特別適用于測量高溫下介質的橫波速度；脈沖法至今仍然在實驗室內(nèi)被廣泛采用，野外小尺度在波速測量中也有所應用；干涉法是目前幾種常用波速測量方法中精度最高的一種方法，可進一步細分為相位比較法和脈沖疊加法。由于觀測技術的發(fā)展，人們越來越重視對地震波速度變化的精確測量與研究。由間接測量波速變化的原理可知，相對波速變化可以通過相對走時變化間接測量(徐薈，2013)。近年來，互相關檢測技術在信號檢測領域獲得越來越多的應用，主要包括從噪聲中提取信號時延估計、速度檢測、距離檢測和系統(tǒng)動態(tài)特性識別等。基于這一原理的地震波速度測量和尾波干涉(CWI)方法近年來已取得了一些突出的成果(Snieder，2006；Stehlyetal.，2007；羅桂純等，2008)，在地震物理預報、監(jiān)測地下應力變化、研究應力和波速變化的對應關系等方面都有重要的應用(羅桂純等，2008)。

2011年中國地震局地球物理研究所和云南省地震局在多震的滇西北賓川大銀甸水庫順利建成了全球首個常規(guī)觀測的大容量氣槍人工震源發(fā)射臺——賓川主動震源發(fā)射臺。發(fā)射臺地處紅河斷裂帶和程海斷裂帶交匯區(qū)中部。作為川滇菱形塊體的西側邊界，紅河斷裂以西是著名的三江地槽褶皺系，以東則被金河—程海—賓川斷裂一分為二。大致以彌渡一帶為界北部稱甘孜—麗江斷塊，南部稱石棉—楚雄斷塊。紅河斷裂北部、南部在大地構造與深部特征、地熱與地球化學場、斷裂結構與力學性質、斷層構造巖發(fā)育特征以及斷裂活動等方面都存在著一系列顯著異常。這些不同時間尺度的差異特征最終都集中體現(xiàn)在沿紅河斷裂帶現(xiàn)代地震活動顯著的南北差異上(張建國等，1993)。張建國等(1993)的統(tǒng)計結果表明，該斷裂北部的地震活動明顯強于南部。程海斷裂是滇西北地區(qū)東緣一條十分醒目的多期活動性大斷裂，主斷裂及其分支斷裂控制著永勝、程海、期納、賓川、彌渡等盆地的發(fā)育。沿著斷裂帶除了發(fā)育諸多晚新生代構造盆地外，還出現(xiàn)了山脊、水系、沖溝的左行走滑位錯，局部地區(qū)還發(fā)育了斷層崖，表明程海斷裂具有左行走滑兼正斷層的性質。程海斷裂早期階段走向總體上近 NS向，傾向西，向南在彌渡附近與紅河斷裂相交，晚期階段走向 NNE向(唐淵，劉俊來，2010)。在滇西地震預報實驗場核心區(qū)域、活動強烈的紅河斷裂北段與程海斷裂之間建立賓川主動源發(fā)射臺，旨在利用人工震源主動向地下發(fā)射地震波，進行地下介質監(jiān)測，變“被動監(jiān)測”為“主動探察”，通過連續(xù)激發(fā)實現(xiàn)動態(tài)跟蹤地震波走時變化，并進一步分析地震波速的時間演化特征，最終為該地區(qū)地震活動后續(xù)發(fā)展趨勢判定提供一定判據(jù)。

賓川主動震源發(fā)射臺于2012年11月投入正常運行，并在其周圍布設了由40套短周期地震儀組成的密集觀測臺陣，形成了一套完整的人工主動震源發(fā)射—觀測系統(tǒng)，并執(zhí)行常規(guī)觀測至今。該激發(fā)系統(tǒng)由4條容量為2 000 in3的單槍組合而成，在平均水深為20 m的大銀甸水庫內(nèi)保持槍與水面的距離為12 m進行激發(fā)，用槍陣高壓氣體的容量和壓力計算(Ronen，2002)得到激發(fā)一次能產(chǎn)生8.91×106J的能量，相當于一次ML0.7天然地震。為了跟蹤分析賓川氣槍源地震波走時差動態(tài)變化特征，本文對2013年以來賓川主動源實驗數(shù)據(jù)進行處理，主要介紹互相關檢測、近場校正等技術在測量賓川主動源地震波速度變化中的運用及初步處理結果。

1互相關時延檢測技術原理

1.1相關函數(shù)

如果x(t)、y(t)是能量有限信號，它們的相關函數(shù)可定義為

(1)

(2)

式中，*表示復共軛。相關函數(shù)是兩信號之間時延為τ的函數(shù)。若x(t) 與y(t) 不是同一信號， 則它們的相關函數(shù)Rxy(τ) 或Ryx(τ) 稱為互相關函數(shù)。 如果x(t) 與y(t) 是同一信號， 即y(t)=x(t)， 此時， 它們的相關函數(shù)Rxx(τ) (或簡寫作R(τ) ) 稱為自相關函數(shù)。即

(3)

在實際應用中，信號x(t)、 y(t)一般是實函數(shù)而不是虛函數(shù)，上述各式仍然適用。在[-∞，+ ∞]的區(qū)間里，功率函數(shù)是不可積的。通常把這類信號的相關函數(shù)定義為

(4)

(5)

(6)

此時有Rxy(τ)=Ryx(-τ)；Rxx(τ)=Rxx(-τ)，即實函數(shù)的自相關函數(shù)是時延τ的偶函數(shù)，τ是所研究兩點間的時間間隔，即兩信號的時延(Rxx(τ)是自相關函數(shù)， Rxy(τ)是互相關函數(shù))。

1.2相似信號的時延檢測

利用相關性計算兩個記錄信號的互相關函數(shù)，便可以較準確地求得兩個記錄信號的時間延遲，從而求得速度分布。其原理是拾震器的信號是由震源信號激振所引起的，源信號在傳播過程中由于路徑等因素的影響，會疊加許多干擾噪聲，但由于這種噪聲通常不與源信號同頻，即源信號與這些干擾波的信號互相關函數(shù)很小，因而將源信號與拾震器信號進行互相關計算時，與源同頻的激振信號便可以識別出來，從而得到其從源到拾震器之間的傳播時間(Takeshietal.，2005；Clifford，1987；Azizul，1981)。即地震波形記錄往往是功率有限信號，如果兩個實數(shù)波形x(t)和y(t)波形高度相似， 只是時間上存在延遲， 即y(t)是x(t)時延τ得到的， 則我們可以通過計算兩個波形的相關系數(shù)函數(shù)來得到時延τ。

x(t)和y(t)的時延相關系數(shù)定義為

(7)

當相關系數(shù)取得最大值時，對應的τ就是兩個信號的時延。通過計算地震波記錄中不同震相的走時變化，就可以確定不同類型波的波速變化。一般都是對于P波、S波以及面波震相來進行計算(徐薈，2013)。

2資料及數(shù)據(jù)處理過程介紹

2.1資料介紹

從目前收集到的激發(fā)記錄可知(圖1)，2013年以來分3個集中激發(fā)時段，分別為2013-01-07～2013-05-08、2013-09-09～2014-05-08、2014-09-08～2014-10-24。根據(jù)氣槍源附近的參考臺記錄可知，2013-01-07～2014-10-24一共開展了100組激發(fā)實驗，共記錄到1 798條氣槍信號。

氣槍源附近布設有40個接收臺站，受客觀因素影響，并不是每個臺站都能記錄到每一次激發(fā)信號，因此本研究選用氣槍源10 km范圍內(nèi)信噪比較高、數(shù)據(jù)記錄缺失較少的6個臺站開展地震波走時變化的跟蹤分析，氣槍源、所選臺站的空間分布情況如圖2a所示，圖2b為2013年1月7日22時40分6個臺站接收到的氣槍信號，由圖可見，6個臺記錄信號清晰、信噪比高。

圖1　2013-01-07～2014-10-24賓川主動源

圖2　臺站分布示意圖(a)及6個臺的原始波形(b)

2.2數(shù)據(jù)處理過程

原始記錄波形包含很多無效信息，不能直接使用。因此，在做相關干涉前，需要對原始波形數(shù)據(jù)進行去均值、去趨勢、頻譜分析及濾波、數(shù)據(jù)篩選、反卷積等預處理，處理流程如圖3所示。

圖3　數(shù)據(jù)處理流程圖

2.2.1去均值和去趨勢

為了避免一些非介質變化信號(如信號直流量、趨勢項)的影響，在對信號進行處理之前，通常需要對信號進行去均值和去趨勢。因為機械振動信號通過傳感器拾取、信號調理和 A/D 采樣轉換成數(shù)字信號，最終都含有一定的直流量。如果積分，會在頻譜處出現(xiàn)一個沖激，并影響其左右的頻譜形狀，產(chǎn)生較大誤差。

另外由于傳感器的零點漂移和現(xiàn)場環(huán)境的影響等，振動信號包含一定的線性項或非線性項的慢變成分，這些信號成分積分后被放大為緩慢變化的趨勢項。且趨勢項疊加在我們所要分析的信號上，對低頻譜產(chǎn)生影響(徐薈，2013)。因此，我們分別使用 SAC 命令 Rmean 和 Rtrend 進行去均值和去趨勢處理。圖4a為2013年1月7日21時59分53264臺站接收到的氣槍原始記錄，圖4b為該條記錄經(jīng)過去均值和去趨勢處理后的波形，經(jīng)過處理有效消除了零漂和趨勢項影響。

圖4　波形預處理過程

2.2.2頻譜分析及濾波

地震信號從震源激發(fā)并傳播到最后被檢波器接收，傳播過程中會受到噪聲的干擾，在低信噪比(SNR)情況下，可能會被噪聲掩蓋。為了盡可能消除噪聲給觀測帶來的影響，需要對采集到的信號進行濾波處理以提高信噪比。濾波之前，首先要對信號進行頻譜分析，判斷波的優(yōu)勢能量集中在哪個頻段，以此來選擇濾波范圍。林建民等(2008)的研究表明氣槍震源激發(fā)的地震波具有豐富的低頻能量，其優(yōu)勢頻率范圍為4～6 Hz。因此，本研究選4～6 Hz帶寬進行帶通濾波。圖4c為2013年1月7日21時59分53264臺站接收到的氣槍記錄經(jīng)過濾波預處理后的波形。

2.2.3反卷積

為了對齊到時達到壓縮子波、消除震源影響的目的，筆者對數(shù)據(jù)進行反卷積運算求其格林函數(shù)。作反卷積前，還需要對經(jīng)過去均值、去趨勢、濾波預處理后的數(shù)據(jù)進行篩選，一方面剔除不正常、有斷記或記錄到天然地震的波形，另一方面剔除信噪比低的記錄。這里我們將水庫邊布設的參考臺接收記錄作為子波(圖5a)，用遠臺接收到的相同時刻的氣槍記錄(圖5b)與參考臺對應記錄作反卷積求其格林函數(shù)(圖5c)，由圖5可見，消除震源影響后震相更清晰。

圖5　反卷積過程

2.2.4互相關時延檢測

重復上述數(shù)據(jù)預處理過程，求得處理時段內(nèi)各臺站每次激發(fā)記錄的格林函數(shù)，分別將該時段內(nèi)各臺站的所有格林函數(shù)疊加作為進行各臺互相關計算的參考模板，例如處理53264臺2013年以來的數(shù)據(jù)，得到1 604條記錄的格林函數(shù)，將這1 604條記錄疊加作為53264臺站的相關檢測參考模板(圖6a)。為了進一步提高信噪比，可根據(jù)具體情況將每天、每周、或每月的格林函數(shù)疊加，本研究選擇疊加每天的格林函數(shù)，圖6b為53264臺2013年1月7日14條記錄的格林函數(shù)疊加結果；最后將每一天的疊加結果與參考模板作互相關計算。具體做法如下：從零時刻開始以一定窗長從兩波形信號中選取出相似波形窗口，然后保持一窗口不動，以不同時延移動另一窗口，并計算不同時延情況下兩相似波形窗口的互相關系數(shù)(圖6c)，互相關系數(shù)Cm(t)最大時所對應的時間延遲就是這兩段相似波形窗口的走時差t1(圖6d)，再以一定步長向后整體移動兩相似波形窗口，重復上述計算過程。為了提高精度，將t1作cos插值得到更高精度的走時差t2(圖6e)，最后結合震相提取初至波的走時數(shù)據(jù)，要求選擇相關系數(shù)最大、t1變化平穩(wěn)的時間窗，如圖6中矩形框所示。

3處理結果

用上述數(shù)據(jù)處理方法處理了2013年以來6個臺站的記錄，測得初至波10-5～10-2的相對走時變化，與已有的人工震源波速變化研究結果吻合。Furumoto和Tsuda(2001)用重復爆炸震源連續(xù)監(jiān)測日本Kanto-Tokai地區(qū)P波速度隨時間的變化，觀測到的波速年相對變化量達到10-3；Silver等(2007)在California 的兩口井中同樣用壓電超聲波震源進行了波速測量實驗，測得量級為10-6s 的波走時變化；Wang等(2008)在云南小江斷裂帶開展了主動震源監(jiān)測介質波速變化實驗，利用電動重錘并用尾波干涉時延檢測方法測得10-3～10-2的相對波速變化，精度為10-4。

圖7a為2013年年初至2014年10月底全時段的初至波走時變化曲線，但由于每年5～9月水位變化，不便于開展氣槍激發(fā)實驗造成缺數(shù)，圖7b～d為2013年以來3個集中激發(fā)時段的地震波走時差曲線放大圖。圖7b顯示，2013年上半年所選6個臺走時差變化趨勢基本一致，其中53278臺變化幅度較其它5個臺要稍大；圖7c為2013-09～2014-05的走時變化曲線圖，由圖可見，在2013-09～2014-02期間53278臺的走時變化趨勢與其它5個臺有所差異；2014年9月以來，6個臺站的走時變化較平穩(wěn)，集中在-0.01～0.01s之間變化，且變化趨勢一致(圖7d)。

4結論與討論

對挑選的6個主動源接收臺站的氣槍信號，采用近場校正、互相關檢測等技術，提取初至波走時差，初步結果與認識如下：

(1)本研究用互相關檢測技術測得10-5～10-2的相對波速變化，對比可知，波速變化在已有相關研究得到的波速變化范圍內(nèi)。

圖7　6個臺站不同時段的初至波走時變化曲線

(2)2013-01～2014-05的初至波走時變化較2014年下半年的大，但變化趨勢基本一致；2014年9月初至10月底，初至波走時變化較平穩(wěn)，集中在-0.01～0.01 s之間變化，且6個臺的變化趨勢一致。

(3)用互相關檢測技術能精確測定地震波走時變化，短期內(nèi)能測量到10-5～10-4的相對波速變化，但是長期的相對波速變化范圍較大(10-5～10-2)，這可能是由于在長時間監(jiān)測過程中地質體或實驗條件的突變引起地震波形畸變，用來互相關計算的隨后波形與參考波形相似性降低所造成。因此，要深入分析走時變化與地下介質應力狀態(tài)的關系還需要充分考慮地質體、實驗條件等(包括水庫水位、接收臺站氣溫、氣壓等)諸多因素的影響。

本文在撰寫過程中得到中國地震局地球物理研究所王寶善研究員；云南省地震局蘇有錦研究員、王彬副局長；中國地震局滇西地震預報實驗場金明培高級工程師、楊軍工程師、李孝賓助理工程師、葉泵工程師的指導和幫助，本文使用了中國地震局地球物理研究所王寶善研究員提供的互相關時延計算程序，在此表示衷心感謝。

參考文獻：

林建民，王寶善，葛洪魁，等.2008.大容量氣槍震源特征及地震波傳播的震相分析[J].地球物理學報，51(1)：206-211.

羅桂純，葛洪魁，王寶善，等.2008.利用相關檢測進行地震波速變化精確測量研究進展[J].地球物理學進展，23(1)：56-62.

唐淵，劉俊來.2010.川滇西部上新世以來構造地貌：斷裂控制的盆地發(fā)育及對于遠程陸內(nèi)構造過程的約束[J].巖石學報，26(6)：1925-1934.

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

徐薈.2013.利用主動震源初至波研究小江斷裂帶淺層波速變化[D].北京：中國地震局地球物理研究所，23-38.

張建國，汪良謀，徐煜堅，等.1993.紅河斷裂深部震源環(huán)境介質力學性質分析[J].地震地質，15(2)：131-137.

Azizul H.Q..1981.An Overview on the time delay estimate in active and passive syste ms for target localization[J].Ieee Transations on Acoustics，Speech，and Signal Processing，229(3)：527-533.

Clifford G..1987.Coherence and time delay estimation[J].Proceedings of the Ieee，75(2)：236-255.

Furumoto J.，Tsuda T..2001.Characteristics of Energy Dissipation Rate and Effect of Humidity on Turbulence Echo Power Revealed by MU radar-RASS Measurements[J].J.Atmos.Solar-Terr.Phys.，63(2-3)：285-294.

Ronen S..2002.Psi，pascal，bars，and decibels[J].Lead Edge，21(1)：60-62.

Silver P.G.，Daley T.M.，Niu F.，etal..2007.Active source monitoring of cross well seismic travel time for stress induced changes[J].Bull.Seism.Soc.Amen，97(1B)：281-293.

Snieder R..2006.The theory of coda wave interferometry[J].Pure Appl.Geophys，163：455-473.

Stehly L.，Campillon M.，Shapiro N.M..2007.Travel time measurements from noise correlation：stability and detection of instrumental time-shifts[J].Geophys.J.Int.，171(1)：223-230.

Takeshi N.，Satoru T.，Teruo Y..2005.Temporal changes in seismic velocity of t he crust around Iwate volcano，Japan，as inferred from analyses of repeated active seismic experiment data from 1998 to 2003[J].Earth Planets Space，57(6)：491-505.

Wang B.S.，Zhu P.，Chen Y.，etal..2008.Continuous subsurface velocity measurement with coda wave interferometry[J].J.Geophys. Res.， 113(B12313)：1-12.

Study on Analysis Method of Travel Time Variations of Seismic

Wave of Active Source in Binchuan

LIU Zi-fen1，SU You-jin1，WANG Bao-shan2，WANG Bin1，YANG Jun1，LI Xiao-bin1

(1. Earthquake Administration of Yunnan Province，Kunming 650224，Yunnan，China)

(2. Institute of Geophysics，CEA，Beijing 100081，China)

Abstract

Firstly，we chose the complete data with high signal-to-noise ratio recorded by six receiving stations of active source in Binchuan County since 2013. Secondly，we measured the velocity change with high precision by using cross-correlation delay detection，near-field correction technique etc.. Finally，the conclusions are obtained as follow：(1)The relative velocity variation wasδv/v≈10-5～10-2in the medium beneath the stations. (2)From early Sep.，2014 to the end of Oct. ，2014，the travel time variations of initial wave was relatively stable which was between -0.01 s and 0.01 s，and the change trend of waveform recorded by six selected stations was same.

Key words：active source；cross correlation delay technique；airgun signal；travel time variation

中圖分類號：P313.22

文獻標識碼：A

文章編號：1000-0666(2015)04-0591-07

*收稿日期：2015-04-07.基金項目：云南省青年地震科學基金項目(201303)、2015年度震情跟蹤定向任務(2015010107)和云南省陳颙院士工作站聯(lián)合資助.