利用主动震源直达波互相关时延检测技术监测小江断裂带浅层地震波波速变化*

徐 荟，刘学军，王 彬，王宝善

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．