利用主動震源直達波互相關時延檢測技術監測小江斷裂帶淺層地震波波速變化*

徐 薈，劉學軍，王 彬，王寶善

0 引言

地震是由在地應力作用下活動斷裂失穩錯動引起的，監測斷裂帶介質應力狀態的變化對研究斷裂帶的演化規律、理解地震的孕育發生過程具有重要意義，但地下介質應力狀態變化的測量是非常困難的。斷裂帶應力狀態的改變會導致地下介質屬性的變化，故可以通過監測介質屬性從而間接監測地下介質應力狀態。地震學家們很早就致力于地殼介質屬性時變性質的研究，提出了包括測量波速比的變化 (Semenov，1969)、剪切波分裂的變化 (Crampin et al．，1990)、尾波Q值的變化 (Jin，Aki，1986)、衰減系數的變化及利用重復人工 (Reasenberg，Aki，1974)或天然地震(Poupinet et al．，1984)測量波速的變化等地震學方法。波速變化研究是其中一個很重要的方法。大量巖石實驗結果顯示波速與作用于巖石的正壓力有關，壓力通過改變巖石中微裂隙的開合來影響波速 (Walsh，1965;Nur，1971)，因此，理論上如果我們能夠連續精確監測波速變化，波速變化可以作為輔助判斷斷裂帶地下應力狀態的指示計。

隨著高重復性人工主動震源和高精度走時測量技術的發展，波速變化的精確測量取得了一些進展 (王偉濤等，2009)。Yamamura等 (2003)用壓電超聲波震源連續監測了日本Miura Bay海岸線附近波速變化，發現波速隨著海潮的應力加載變化而規律地變化。Silver等 (2007)在California的兩口井中同樣用壓電超聲波震源進行了波速測量實驗，測得量級為10-6s的波走時變化，其中一口井測得的走時變化與大氣壓呈正相關，而另一口井的走時變化與大氣壓呈負相關，其原因為多孔彈性介質的近場和遠場效應，兩口井處相對波速變化對大氣壓的敏感系數分別約為10-6Pa-1和10-7Pa-1。Niu等 (2008)通過對SAFOD鉆井內波速測量，發現波速變化與大氣壓呈負相關，并且發現實驗期間的兩次地震前波走時都有異常變化，可能是由于震前地下介質應力發生了變化。Wang等 (2008)在云南小江斷裂帶開展了主動震源監測介質波速變化實驗，并用尾波干涉時延檢測方法測得10-3～10-2的相對波速變化，精度為10-4，相對波速變化與大氣壓有很好的負相關性，且對大氣壓的敏感系數為10-6Pa-1。尾波由介質內部非均勻性顆粒對彈性波的多次反射、散射等產生，傳播路徑復雜，所測波速變化反映的是存在強散射區域介質的平均波速變化，而很難確定變化發生的相對位置。與尾波相比，直達波傳播規律很清楚，其波速變化可以反映其傳播路徑上的波速變化，因此筆者使用直達波互相關時延檢測方法監測云南小江斷裂帶附近地震波波速的連續變化。

1 實驗簡介

實驗地點位于云南省昆明市小哨臺站附近，實驗時間從2006年4月7日06時30分至5月8日23時30分，近一個月。實驗選址有如下3方面考慮:(1)觀測點所處的大構造背景是小江斷裂帶，小江斷裂帶形成與發育歷史久遠，無論是斷塊垂直差異運動還是走向滑移運動都非常突出，地震活動十分頻繁;(2)測線布設在小哨臺站東北方一塊較平緩的山坡上，其覆蓋層厚度不超過2 m，地表基巖出露較多;(3)實驗地點附近有多種觀測數據可以用來和波速變化進行比對標定。小哨臺 (XS)有一口地下水位觀測井，井深2 156 m，在小哨臺以西33 km處的昆明臺 (KM)有連續重力觀測系統，以東32 km處的嵩明臺(SM)有氣象三要素 (氣壓、溫度、降水量)觀測系統，以南36 km處的宜良臺 (YL)有水管儀形變觀測系統。這些觀測臺站都在實驗地點40 km范圍內。綜合上述3點，此次實驗地點是一個非常理想的野外實驗場。Wang等 (2008)給出了本次實驗更詳細的描述，實驗位置及場地觀測系統布設如圖1所示。

實現地下介質波速變化的精確測量有兩個要素:有性質單一、重復性好、穩定性好的震源;有高精度的數據記錄系統。

震源采用了GISCO公司生產的ESS200電動落錘，重錘質量為99．8 kg，完全由電力控制自動激發，落錘提升高度由鉸鏈決定，每次提升高度一致，約為1 m，激發能量相同，約為30 kJ，以保證有很好的震源重復性。由于實驗期間正處小哨地區的風季，為避開風擾，同時考慮到電落錘電源的充電時間 (充滿電的電落錘可以連續擊震120次左右，充滿電一次需4～6 h)，選擇在每天風動較小的6個固定時段 (00:30，01:30，06:30，07:30，22:30，23:30)啟動電落錘。每次啟動，電落錘在12 min內完成30次垂直擊震。

測線由48道地震檢波器組成，電壓靈敏度為22 V/m·s，主頻率為60 Hz，頻帶為10～200 Hz，道間距5 m，測線總長度235 m，第1道距震源15 m，共48道，道間距5 m，如圖1中實直線所示。為保證高數據記錄精度，數據采集器方面使用了Geometrix NZ StraII高采樣淺勘儀，采樣率為32 000 Hz。數據采集器會將震源每次啟動產生的30條擊震信號自動線性疊加。

由于第1道距震源很近，可近似認為第1道的記錄波形就是震源信號波形。圖2a為兩次激發記錄的震源信號波形，圖2b為所有震源信號波形的互相關系數，從圖中可看出互相關系數大部分在0．99以上，說明震源具有高度重復性，從震源方面保證了此次實驗的可靠性。

2 方法原理

2． 1 互相關時延檢測方法

假設地震波沿某一固定路徑傳播的基準走時為t0，走時變化為τ，則地震波速度的相對變化與走時的相對變化關系為由式 (1)可以看出，相對波速變化可以通過相對走時變化間接測量。按照一般方法，精確測量走時變化的關鍵在于準確判斷波初動到時。由于各種干擾噪聲及能量衰減等因素，實際記錄波形初至的精確識別難度大且精度不高。互相關時延檢測技術通過計算兩個記錄波形的相關函數，利用有效信號確定性和噪聲信號隨機性的特點，可以精確測量兩波形信號走時差，從而解決了波速測量中精確拾取初至這一難題。

用互相關時延檢測方法計算兩個波形信號之間走時差，首先選用適當窗長從兩波形信號中選擇相似波形窗口，其中一窗口保持不動，以不同時延移動另一窗口，并計算其互相關系數，互相關系數最大時對應的時間延遲即是此兩波形信號走時差。

計算波的連續走時變化有兩種方法。一種是普通互相關時延檢測方法，另一種是前后互相關時延檢測方法。普通互相關時延檢測方法是將第一個波形指定為參考波形，然后用互相關時延檢測方法分別計算其他波形相對于第一個波形的走時差。前后互相關時延檢測方法不設定唯一參考波形，而是隨時間變化將參考波形逐一向后輪換，這樣計算得到的是前后兩地震波形的走時差，最后將這些走時差進行累加，便可得到某一波形相對于第一個波形的走時差。理論上，這兩種方法都可以計算得到波形連續走時變化，但是前后互相關時延檢測方法的結果更可靠一些，其原因可能是普通互相關時延檢測方法用來互相關計算的隨后的波形和第一條波形相似性越來越差 (圖3a)，不利于提取最大相關系數點而得到精確到時，而前后互相關時延檢測方法用來互相關計算的相鄰波形相似性較高 (圖3b)。

2． 2 近場校正消除震源同步觸發誤差

直達波互相關時延檢測技術的誤差來源于兩方面:互相關時延檢測方法自身誤差和定時中的時間同步誤差。時間同步誤差σ又包括震源同步觸發誤差σsource和道間同步采集誤差σsample(李宜晉，2011)，則

假設遠近場波測量走時為T1，T2，則分別表示為

其中，t1，t2分別為遠近場波實際走時，σsample1為遠場道間同步采集誤差，量級為10-12s，可近似忽略。

如果進行近場校正

可得到完全去除震源同步觸發誤差的遠近場之間走時差。

3 波速變化精確測量

實驗共獲得48道每道221條自動疊加后的垂直分量波形，每條波形的記錄長度為0．512 s。筆者挑選2006年4月7日22時30分第8、18、28道的波形，并對波形進行頻譜分析 (圖4a)，從圖中可看出信號頻率主要集中在10～80 Hz。再將波形信號經過10～80 Hz的4階Butterworth帶通濾波以降低噪聲，然后手動剔除信噪比差和記錄明顯錯誤的個別波形，用于下一步時延計算。圖4b為第8、18、28道波形濾波后波形圖，從圖中可以看出，濾波后波形信噪比明顯增加，可以清晰看到波形中的直達波部分 (圖4b中陰影部分)。

筆者選用前后互相關時延檢測方法計算波的連續走時變化。在計算前后兩波形的走時差時，首先從零時刻開始以0．05 s窗長，從兩波形信號中選取出相似波形窗口，然后保持其中一窗口不動，以不同時延移動另一窗口，并計算不同時延情況下兩相似波形窗口的互相關系數，互相關系數Cm(t)最大時所對應的時間延遲就是這兩段相似波形窗口的走時差。再以步長0．001 s向后整體移動兩相似波形窗口，重復上述計算過程。由于此次實驗是用直達波的走時變化來表征淺層介質中地震波的走時變化，直達波形寬度為0．035 s。將中心點位于直達波形范圍內的所有相似波形窗口的走時差求平均，便得到前后兩直達波波形的平均走時差，即前后兩波形走時差。最后將這些走時差進行累加，得到后面波形相對第一個波形的走時差，即得到了波的連續走時變化。

將除第1道的其它道與第1道進行近場校正，完全去除震源同步觸發誤差，得到第2～48道近場校正后連續走時變化。選取的直達波波形的中心點時刻即為直達波平均走時，走時變化除以平均走時可得到相對走時變化，再利用相對波速變化與相對走時變化的關系得到相對波速變化。

第2～48道連續相對波速變化趨勢基本一致，但幅值不同，故將每10道進行平均，圖5a即為此次實驗的760 h期間第2～9、10～19、20～29、30～39道近場校正后平均連續相對波速變化。相對波速變化隨時間的長周期趨勢在下雨期后明顯上升 (實際為一小段時間，這里近似為時間點500 h。根據長周期趨勢的斜率，將一個月的實驗期間大致劃分為兩個階段:階段Ⅰ下雨前 (0～500 h)和階段Ⅱ下雨后 (500～760 h)。從圖中看到，每10道平均連續相對波速變化的走勢基本一致，但第2～9道的相對波速變化比其他道更劇烈。各道相對波速變化的標準差如圖5b所示，前10道標準差明顯高于后30道，后30道標準差大小基本相近，其原因是離震源較近的道所測得的波速變化反映的是淺層介質情況，離震源較遠的道測得的波速變化反映的是深層介質情況，相比于深層介質，淺層介質強度較小，大氣壓變化相同時，淺層介質中裂紋和孔隙開合程度的變化更大，引起的波速變化更大，所以近道的波速變化更為劇烈。

互相關時延檢測方法產生的誤差在理論上存在一個下限值，被稱作Cramer-Rao下限值 (Silver et al．，2007;Wang et al．，2008;Wang et al．，2009)，表示的是波形時延估算的最小標準差

其中，f0為信號中心頻率，T為所選相關窗口長度，B為頻寬比，ρ為波形的相關系數，SNR為信噪比。由式 (5)可以看出，信噪比是決定測量誤差的重要因素，提高信噪比可以大大減小走時測量的誤差。

此實驗中 f0≈50 Hz，T≈0．05 s，B≈0．4。筆者計算了各道的信噪比SNR，最低信噪比約為100，SNR?1。ρ用所有連續相似波形窗口分別對應的最大互相關系數Cm(t)的平均值Cm(圖3b)來表示，大部分的相關系數平均值Cm都大于 0．99，近似認為ρ=1。故式 (5)簡化為

通過計算，各道相對波速變化的平均標準差為4．6×10-3，而我們的測量精度達到10-4，所以可以由我們的測量系統進行精確測量。

4 波速變化原因探討

Silver等 (2007)，Sens-Sch?nfelder和 Wegler(2006)研究表明波速變化與地下水位聯系緊密。因此，筆者對波速變化與地下水位的相關關系進行了研究，發現波速變化與地下水位有一定的相關性，但相關性不是很好。而大氣壓和固體潮是導致水位變化的主要原因 (Spane，2002)，所以筆者進一步探究大氣壓和固體潮與波速變化的相關性。研究后發現波速變化和大氣壓相關性很好，而和固體潮的相關性不明顯，這與 Wang等(2008)對該地區的研究結果一致。圖6中顯示的是實驗期間第30～39道平均連續相對波速變化與對應時段內小哨臺 (XS)記錄的地下水位、嵩明臺 (SM)記錄的大氣壓以及用軟件MAPSIS計算得到的理論的固體潮引起的面應變 (eNS+eEW)。從圖中可以看出，地下水位主要受大氣壓和固體潮影響，大氣壓和波速變化相關性很好，但固體潮和波速變化的相關性不明顯，所以地下水位與波速變化呈一定相關性但相關性不會特別好。

大氣壓和波速變化的相關關系很明顯，第10～19(圖7c)、20～29(圖7e)、30～39(圖7g)道在階段Ⅰ下雨前 (0～500 h)相對波速變化與氣壓呈正線性相關，第2～9(圖7b)、10～19(圖7d)、20～29(圖7f)、30～39(圖7h)道階段Ⅱ下雨后 (500～760 h)相對波速變化與氣壓呈負線性相關，只有第2～9道在階段Ⅰ下雨前相對波速變化與氣壓相關關系離散 (圖7a)，其原因可能為實驗點處覆蓋有薄風化層，第2～9道接收到的波是經風化層傳播的，而遠道接收到的波則穿過了風化層到達下方巖石層中傳播 (圖8)。風化層性質結構比較復雜，故下雨前第2～9道波速變化與大氣壓相關關系離散，下雨后風化薄層連通性變得非常好，所以波速變化與大氣壓關系也變得很明顯，且對大氣壓敏感度很高。

應用一元最小二乘線性擬合大氣壓變化對波速變化的影響關系，線性擬合表達式為

其中，α1為相對波速變化對大氣壓的敏感系數，α2為常數。根據波速變化隨時間的長周期趨勢的斜率，筆者將實驗期間大致劃分為兩個階段，所以也分別對這兩個階段進行線性回歸擬合段 (圖7中直線所示)，擬合所得的α1值如表1所示。實驗點處相對波速變化對大氣壓的敏感系數估計值約為10-6Pa-1，這與 Yamamura等 (2003)、Silver等 (2007)、Wang等 (2008)實驗結果一致。

表1 相對波速變化對大氣壓敏感系數Tab．1 Sensitive coefficient of the relative velocity variation on barometric pressure

地下介質由骨架和填充了氣液體的微裂隙組成，微裂隙間通過更細微的毛細管道連通。介質中波速大小主要和有效應力Pe有關 (Pe=Pc-Pp，其中Pc為圍壓，Pp為孔隙壓)。大氣壓對有效壓力的影響比較復雜，會同時影響圍壓和孔隙壓:大氣壓增大導致圍壓增大，同時也通過擠壓微裂隙中氣液體增大了孔隙壓。下雨前，大氣壓增大，圍壓增大占主導，有效應力增大，孔隙率減小，波速增大，于是表現出波速和大氣壓呈正相關;下雨后，大氣壓增大，由于下雨導致孔隙連通性大大增強導致孔隙壓增大占主導，有效應力減小，孔隙率增大，波速減小，于是表現出波速和大氣壓呈負相關。Wang等 (2008)則發現該實驗點處相對波速變化與大氣壓在下雨前后一直呈很好的負相關。因為尾波所測波速變化反映的是存在強散射區域介質的平均波速變化，波經過強散射之后還能從地面被接收到說明強散射區域深度比較淺，介質孔隙連通性很好，下雨后大氣壓增大時孔隙壓增大占主導，有效應力減小，波速減小，于是表現出波速與大氣壓呈負相關。

5 結論

筆者在云南小江斷裂帶附近進行了一個月的淺層地震波波速變化監測實驗，并利用互相關時延檢測技術計算得到了精度為10-4的波速相對變化10-3～10-2。實驗點處地下介質中不同深度的波速變化趨勢相同，但變化劇烈程度不同。離震源較近道相對波速變化與較遠道相比更劇烈，原因是近道反映的是淺層介質中波速變化情況，淺層介質強度較小，大氣壓變化相同時，淺層介質中裂紋和孔隙開合程度的變化更大，引起的波速變化更大，所以近道的波速變化更為劇烈。

波速變化和大氣壓有很好的相關性，下雨前(0～500 h)相對波速變化和氣壓的變化呈正線性相關關系，而下雨后 (500～760 h)相對波速變化和氣壓的變化呈負線性相關關系。下雨前，大氣壓增大時圍壓增大占主導，有效應力增大，波速增大，故波速變化與大氣壓表現為正相關關系。下雨后，介質孔隙連通性增強，大氣壓增大時孔隙壓增大占主導，有效應力減小，波速減小，故波速變化與大氣壓表現為負相關關系。該地區相對波速變化對大氣壓的敏感系數估計值約為10-6Pa-1。

李宜晉，辛維，王寶善，等．2011．利用高采樣率數采實現地震波速變化高精度測量［J］．地震地磁觀測與研究，32(2):34-40．

李宜晉．2011．小尺度人工震源地震波速變化觀測系統的技術研究［D］．北京:中國地震局地球物理研究所．

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

王彬．2009．利用多種震源測量介質波速變化的實驗研究［D］．合肥:中國科學技術大學．

王偉濤，王寶善，葛洪魁，等．2009．利用主動震源檢測汶川地震余震引起的淺層波速變化［J］．中國地震，25(3):223-233．

楊微，葛洪魁，王寶善，等．2010．由精密控制人工震源觀測到的綿竹5．6級地震前后波速變化［J］．地球物理學報，53(5):1149-1157．

Crampin S．，Booth D．C．，Evans R．，et al．．1990．Changes in Shear Wave Splitting at Anza near the Time of the North Palm Springs Earthquake［J］．J．Geophys．Res．，95(B7):11197-11212．

Jin A．，Aki K．．1986．Temporal Change in Coda Q before the Tangshan Earthquake of 1976 and the Haicheng Earthquake of 1975 ［J］．J．Geophys．Res．，91(B1):665-673．

Niu F．L．，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］．J．Geophys．Res．，76(8):2022-2034．

Poupinet G．，Ellsworth W．L．，Frechet J．．1984．Monitoring Velocity Variations in the Crust using Earthquake Doublets:An Application to the Calaveras Fault，California［J］．J．Geophys．Res．，89(B7):5719-5731．

Reasenberg P．，Aki K．．1974．A Precise，Continuous Measurement of Seismic Velocity for Monitoring in Situ Stress［J］．J．Geophys．Res．，79(2):399-406．

Semenov A．N．．1969．Variation in the Travel Time of Transverse and Longitudinal Waves before Violent Earthquakes［J］．Bull．Acad．Sci．USSR，Phys．Solid Earth，3:245-248．

Sens-Sch? nfelder C．，Wegler U．．2006．Passive Image Interferometry and Seasonal Variations of Seismic Velocities at Merapi Volcano，Indonesia［J］．Geophys．Res．Lett．，33(21)，L21302，doi:10．1029/2006 GL027797．

Silver P．G．，Daley T．M．，Niu F．L．，et al．．2007．Active Source Monitoring of Cross-well Seismic Travel Time for Stress-induced Changes［J］．BSSA，97(1B):281-293．

Silver P．G．，Daley T．M．，Niu F．L．，et al．．2007．Active source monitoring of cross-well seismic travel time for stress-induced changes［J］．BSSA，97(1B):281-293．

Spane F．A．．2002．Considering Barometric Pressure in Groundwater Flow Investigations［J］．Water Resour．Res．，38(6):14-1-14-18．

Walsh J．B．．1965．The effect of Cracks on the Compressibility of Rocks［J］．J．Geophys．Res．，70(2):381-389．

Wang B．S．，Zhu P．，Chen Y．，et al．．2008．Continuous Subsurface Velocity Measurement with Coda Wave Interferometry［J］．J．Geophys．Res．，113(B12)，B12313，doi:10．1029/2007 JB005023．

Yamamura K．，Sano O．，Utada H．，et al．．2003．Long-term Observation of in Situ Seismic Velocity and Attenuation［J］．J．Geophys．Res．，108(B6)，B62317，doi:10．1029/2002JB002005．