台差法提高低信噪比下走时变化探测精度的研究
黄亦磊+运乃丹+周仕勇
摘要：为研究噪音对主动源信号走时变化提取精度的影响，采用对比数值实验的方法。按照常规流程，首先用反褶积去除震源激发时间的不确定性干扰，压缩子波以得到更加尖锐的震相，由于发现反褶积会降低信号的信噪比，于是提出了台差法来代替反褶积矫正气枪激发的不确定性，再结合干涉法提取走时变化，并定量对比了两者的走时变化提取精度，最后尝试用台差法提取主动源走时的日变化。结果表明，台差法可以有效矫正震源激发时刻的影响，且不降低信噪比，从而也能获得较高走时变化提取的精度。
关键词：主动源；走时变化；走时台差法；信噪比
中图分类号：P315.31文献标识码：A文章编号：1000-0666（2017）04-0595-10
0引言
地震波信号是照亮地球内部结构的一盏明灯。复杂的地震波信号是由断层的破裂过程和介质格林函数耦合而成。传统的地震学分为研究震源运动学和地球结构研究（Lay，Wallace，1995）。从地震信号研究震源S（t）或者介质结构G（t）就必须使两者解耦。因此，大多数震源或者介质的研究中，都必须假设另一者有先验的模型。
由于主动源具有精确的激发时刻、激发地点，为地下介质的研究提供了良好的可重复震源，从而使精细获取介质结构变化成为可能。Niu 等（2008）利用压电陶瓷作为主动源，发现在实验场附近2个小地震发生前，压电陶瓷产生地震波信号走时有显著的延迟。在物理上，由于地震发生之前都会有裂隙和孔隙压力的变化，这个变化在介质波速上有所体现（Simmons，1964；Yukutake，1988），当这种精细变化被主动源信号所捕获时，将有助于我们认识地震孕育的物理过程。压电陶瓷具有频率高、衰减快、传播距离十分有限的特点，不能在大区域的研究中作为广泛应用的主动源。因此，陈颙院士在其所主导建设的宾川地震信号发射台中，引入了4枪组合的气枪激发源作为主动源，从而克服了传播距离有限的问题（罗桂纯等，2006；王宝善等，2011；王伟涛等，2009；杨微等，2013）。
目前，宾川发射台气枪实验已顺利进行多年，积累了大量的实验数据。如何提取高精度的主动源信号走时变化是众多学者关心和研究的问题。笔者首先介绍目前提取走时变化常用的干涉法和处理流程，之后对提取走时变化流程中反褶积这一步骤提出可能存在放大噪音而影响走时变化提取精度的问题，并进行相应的数值实验来验证这一想法，进而提出一种可能替代反褶积的台差法来矫正震源的影响，并进行了相应的实证。
1干涉法提取走时变化简介
波速测量关键的问题在于走时的获取，对于相对走时变化在0.01%～1%量级的精细变化（Wang et al，2008；刘自凤等，2015）探测中，传统手动拾取震相到时的方法已经不再适用。由于主动源气枪激发的重复性，波形上也表现出一定的相似性，基于互相关提取相对走时变化的方法比较适用（刘自凤等，2015；王彬等，2012）。本文以Wang 等（2008）所提出的干涉法为基础，开展信号走时提取。
1.1互相关函数
对于两个能量有限、定义在[t1，t2]上的有限长地震实信号x（t）和y（t），相关函数定义为：Rxy（τ）=∫t2t1x（t）y（t-τ）dt（1）其中：τ为时间延迟。物理意义上，相关函数表征了将一个信号相对另一个信号移动之后的相似性，为了方便比较，引入无量纲归一化的相关系数ρxy：ρxy（τ）=Rxy（τ）∫t2t1x2（t）∫t2t1y2（t）（2）互相关系数在[-1，1]间变化，互相关系数值越大说明2个信号越相像。2个相同的信号在τ=0时，互相关系数取1。
1.2干涉法
这里的干涉法特指为气枪设计的一种基于互相关提取走时变化的方法。在气枪主动源的走时变化探测过程中，存在以下2个问题，必须设计与之相应、适合气枪主动源提取走时变化的方法：
（1）气枪信号由很多震相组成，不同震相代表了地震波在地下介质中传播的不同路径，因此形成了对地下不同的采样空间。不同路径上的走时变化可能不同，所以有必要截取出一定长度的时窗，使之只包含一个特定的震相，以保证窗内所提取的走时变化是稳定而一致的。
地震研究40卷第4期黄亦磊等：台差法提高低信噪比下走时变化探测精度的研究（2）地震仪的采样率为100 Hz，对应的时间采样间隙0.01 s，然而介质速度变化引起的走时变化可能小于时间采样间隙，因此我们有必要对记录的地震信号进行有效插值，将走时测量精度大幅提高至0.001 s以上。
干涉法的原理如图1所示。首先我们需要得到一个参考波形，将参考波形和待测信号分成若干个小时窗，每个小时窗可以左右移动一定的步长，将移动过程中最大的互相关系数（图1c）对应的时间移动定为该时窗的相对时移（图1d）。每一个时窗都对应一个时间延迟值，我们需要選定一个波形显著、互相关系数较大且平稳的时间段，取这段时间中波形的互相关系数（图1c）作为走时变化的权重，加权求平均值，作为该道信号的走时变化。
从时移量得到速度变化有很简单的数学关系（Schaff，Beroza，2004）。由于信号的传播距离s（台源距）可以确定，假设介质速度原来是v，对应的走时为t，速度变化为dv，走时变化为dt，则有：
s=v·t=（v+dv）（t+dt）（3）忽略二阶小量dv·dt，则有：dvv=-dtt（4）从干涉法的实现过程可以看出，信号走时变化提取或多或少会受到信号处理流程的影响。现列举如下，以提示读者之后研究的注意：
（1）地震信号的前期处理。包括去均值、去除线性趋势、滤波和反褶积。其中去均值、去除线性趋势、滤波都是提高信号信噪比的有效手段。而反褶积则能有效矫正震源发震时刻不确定性的影响，并压缩子波得到更加尖锐的震相序列，但同时会降低信噪比。
（2）干涉法中时窗大小选取以及滑动步长选取。时窗大小以及滑动步长的选取带有经验性和探索性，也依据不同研究问题的目标而定。
（3）选取哪一段波形对应的走时变化做平均得到最后的走时变化。
如何选取时窗滑动步长、滑动步长选多大以及最后算平均走时变化的时窗，具有经验性和探索性，这是比较主观的影响因素，并非本文所要研究的问题。在之后的研究中，我们尽可能控制使窗口长度以及滑动步长等主观因素一致。主要定量研究噪声及滤波和反褶积等步骤对走时变化提取精度的影响，并尝试提出一种保持原信号信噪比，但又能一定程度上矫正震源时间不确定性的台差法，并对这种方法的有效性开展了实证研究。
2干涉法提取走时变化可能影响因素探究
影响一个物理量准确测量的原因有两种：一是客观存在的误差和噪音；二是主观人为因素。
仪器的误差无法消除，自然的噪音也客观存在，因此要在提取走时的过程中，尽量降低噪音，并使用合适的信号处理流程。笔者将重点研究在数据处理流程中可能会影响走时变化提取的因素。
刘自凤等（2015）提出了宾川主动源气枪定时变化提取的一个流程，如图2所示。这个流程是目前大多数走时变化提取所遵循的流程，这个流程看起来是比较合理的，但是刘自凤等（2015）并没有做数值方面的测试，也没有做主要步骤对提升结果的可靠性的分析。接下来，笔者将从理论测试出发，研究在不同的噪音水平和处理流程中（主要是滤波和反褶积）对走时变化提取结果的影响。
2.1噪音背景下反褶积对走时变化提取的影响
在与地震信号同频率的噪音存在的背景下（滤波无法消除噪音），实际地震信号对震源时间函数反褶积，噪音的影响可能会被放大。为了研究这一问题，需如下数值实验。
选择2013年前80天53281台（位置如图3所示）数据作为研究对象，将80天数据叠加，如图4a 所示。并把参考台CKT0（位置如图3所示）相同时间段的数据叠加作为震源时间函数（图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组实验数据，总结如下：
（1）组内实验结果对比，可以发现噪音越高，对走时变化提取影响越大。符合常规预期。
（2）滤波可以有效提高信噪比，降低走时变化提取的不确定性。符合常规预期。
（3）反褶积降低信噪比，反褶积之后提取走时变化的结果都差于不做反褶积，反褶积实质上起到了一个放大噪音的效果。
在此基础上，我们以53281台80天叠加信号为基础，添加噪音，使其信噪比从-10到15连续变化，对每一个信噪比下的信号滤波之后，做反褶积（尝试0.000 1，0.001和0.01的反褶积waterlevel水平）和不反褶积（只滤波）处理，之后提取走时变化，仍然以原始数据叠加而未添加噪音的信号作为参考模板。取每条信号的信号段（8～12s）提取的走时变化的平均值，作为该信号的走时变化。结果如图7所示，可以看出，所有反褶积处理之后的信号，走时变化提取结果都比只滤波处理的结果差，而且，随着反褶积waterlevel水平的增加，反褶积之后提取的走时变化结果也会更加好且稳定。
综上所述，我们在图7中定量展示了在存在不同水平噪音时，反褶积这一处理流程可能带来的走时提取误差。反褶积是一把双刃剑，它可以消除震源的影响，压缩子波，但是也会降低信号的信噪比，给走时变化提取带来不确定性。且高噪音水平要求反褶积时用更高的反褶积水准值。
3替代方案——台差法及其有效性实证
3.1反褶积的可能替代方法——台差法
提出用反褶积处理信号的初衷之一，是矫正气枪源的不稳定影响，比如激发时刻不同步以及水位变化等。反褶积可以去除震源的影响，留下介质的格林函数。但是，上述的数值实验说明，在存在噪音的情况下，反褶积会降低信噪比，从而导致变化提取不准确。另外，由于参考台得到的信号并非真实的震源时间函数，可能还需考虑方位带来的辐射花样问题，有一定的局限性。
因此，我们提出用台差法来代替反褶积以消除震源的影响，从而避免由于信噪比低而引起提取走时变化的精度问题。台差法的想法比较简单：震源如果存在激发时刻不准确或水位变动的影响，这种走时变化是一种整体时移，在每一个台站的记录都是近乎相同的，因此利用非常靠近震源的参考台进行自身干涉而提取到的走时变化，可以认为就是源的变化所导致的，然后再用其他台站提取的走时变化结果减去参考台提取的走时变化结果，这样便可矫正源的影响。总结台差法的步骤如下：（1）对信号去均值，去除线性趨势，滤波；（2）参考台用干涉法提取走时变化；（3）待研究台站用干涉法提取走时变化；（4）待研究台站走时变化减去参考台走时变化得到最终的走时变化。
为了验证台差法的正确性，我们做了如下实验：选取2013年前80天的53281台和CKT0台数据，分2组做对比实验：第一组用53281台单枪数据对CKT0同一枪数据做反褶积（反褶积水准值0.001），然后单天数据叠加，提取走时变化（走时变化用8.5～10.5 s时段走时变化平均）；第二组，对CKT0台和53281台数据单天叠加，然后分别用干涉法提取CKT0台和53281台走时变化，然后用台差法（53281台走时变化用9～11 s时段走时变化平均，CKT0台走时变化用0.5～1.5 s时段的走时变化平均）处理数据。上述2组实验在提取走时变化时的模板均为2台在这80天内的所有数据的叠加，结果如图8所示。可以发现，CKT0台和53281台提取的走时变化有很好的相关性，说明源的误差确实在这两个台都有所体现，相减之后，得到的走时变化非常平缓，在这个阶段内，并无大的环境因素改变，我们有理由认为这段时间的介质速度基本不变。对比反褶积提取的走时变化结果，可以发现台差法提取的走时变化结果比反褶积提取的走时变化结果更加稳定。
3.2实际资料处理
为了进一步验证本文提出的台差法的有效性，本文研究了气压以及地震可能对介质速度的影响（Spane，2002）。首先，选取了2015年11月20—28日的集中实验阶段作为研究对象，在此阶段，气枪激发比较密集，且有气压记录，适合研究气压对走时变化的影响。首先选择53264台、53274台、53281台和53266台作为研究对象（图3），
首先用台差法对单枪数据提取走时变化，进行加权平均的走时变化的时间范围分别为3～5 s、3～5 s、8～10 s和4～6 s，模板为9日内所有数据的叠加，4个台的走时变化和气压变化如图9所示。可以发现，在误差范围内，这4个台的走时变化都与气压有非常一致的变化，但是走时变化似乎领先气压变化一些相位。为了进一步观察他们的相关关系，本文细化了实验，选取53274和53281台作为研究对象，将这两个台在集中实验期间，每天同一小时、同一分钟的信号叠加，然后提取走时变化，进行加权平均的走时变化的时间范围仍为3～5 s和8～10 s，模板同样为9日内所有数据的叠加，得到的将是一天之内的走时变化。从图10可以看出，这两个台的走时变化随着气压都有正向相关的趋势，但是走时变化都领气压变化一些相位。
除此之外还计算了一天之内气压对于走时变化的具体数值影响，从走时变化结果中取出最大值和最小值做差，得到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）。
实验室进行的岩石实验表明，气压的变化会改变岩石内部裂隙中物质的密度，从而改变地震波在岩石中的传播速度，使波的走时发生改变（Nur，1971；Walsh，1965），但实际地球内部的岩体结构更加的复杂，气压与走时变化之间的具体关系仍待进一步探究。另外，除了气压对走时变化有影响之外，固体潮和温度（Leary et al，1979）等因素也可能通过影响介质速度而使走时变化产生波动，目前我们还没有进行这方面的探究，这也是接下来的研究内容之一。
4讨论
研究地震孕育过程中介质可能发生的变化，对于我们认识、理解地震发生过程具有重要意义。如果能有效地提取出介质的变化，并研究清楚这些变化与地震之间的关系，则我们在地震前兆的认识上将更进一步。
反褶积理论上能去除源的干扰。但实际上，由于气枪水下激发产生的地震波的机理尚未研究透彻，只能用距离气枪西面50 m处的一个参考台来近似源，在现有的条件下是一个不错的选择。另在一个问题，即参考台虽然距离气枪很近，但是受到震源辐射花样方位角的影响，因此方位角偏离参考台比较大的台站，用参考台信号反褶积的影响程度仍有待考究。
至于台差法能否准确地提取出介质的走时变化，有待进一步的验证。台差法目前来看有一个缺点，就是不能很好地将参考台信号对应到待研究的台站。由于参考台P波、S波混叠在一起，而在待研究台站可以明确选择某一震相，这样做差可能难以解释他们的对应关系，但是震源的不确定性在各个震相应该都有体现，所以我们相信可以部分消除源的影响。台差法的第二个问题在于射线路径的问题，台差法可以明确的矫正气枪的触发误差、水位变化带来的误差，但是对于不同的台站，射线路径不同，在从震源到参考台部分的走时，我们没有考虑不同台站是不一样的，这一部分误差不能保证完全消除，有待进一步研究。但是鉴于CKT0台离源非常近，我们认为影响不是很大。
进一步，我们将研究走时变化提前于气压变化的原因，并将用理论计算得出固体潮的理论值进行对比。另外，下一步将回溯地震，观察在其发生前是否有走时变化的突变。
5结论
在目前无法理论推导震源函数的情况下，前人想到用离震源非常近的参考台的垂向分量作为震源时间函数，来进行反褶积，在高信噪比的情况下，是一种可行之举。但是我们通过实践发现在存在较高噪音的情况下反褶积可能导致测量误差加大，使得走时变化的提取不准确，进而设计了理论实验和实际数据的实验，发现反褶积这个过程在走时提取中确实会降低信号的信噪比，从而影响走时变化的精确提取。然后，我们定量研究了噪音水平在走时变化提取中的影响。结果表明，反褶积只有在高信噪比的情况下（8以上）才不会影响走时变化提取的精度。而且原始信号的信噪比越低时，反褶积可能将信号的噪音放大越严重。所以当原始信号的信噪比相对较低时，应该提高反褶积的水准值。实验发现，不進行反褶积的信号与反褶积之后的信号提取走时变化相比，在信号信噪比低于8时，未反褶积的信号提取的走时变化更加准确和平稳，且对噪音的抵抗能力更强。
为了解决这个问题，笔者提出了一种可能的解决方法，也就是用参考台自身信号的走时变化来矫正远台的走时变化，这种方法，从目前的结果来看，也不失为一种矫正源效应的方法。另外，进一步进行了该方法的实证，结果发现，气压变化与气枪走时变化有很强的相关性，但是气枪走时变化总是领先气压变化，猜想可能和气温、固体潮等因素相关，这个负载的物理问题有待进一步研究。
参考文献：
刘自凤，苏有锦，王宝善，等.2015.宾川主动源地震波走时变化分析方法研究[J].地震研究，38（4）：591-597.
罗桂纯，王宝善，葛洪魁，等.2006.气枪震源在地球深部结构探测中的应用研究进展[J].地球物理学进展，21（2）：400-407.
王宝善，杨微，王伟涛，等.2011.宾川气枪信号发射台性能分析[C]//中国地球物理学会年会第二十七届年会论文集.
王彬，楊润海，王宝善，等.2012.地震波走时变化精确测量的实验研究[J].云南大学学报：自然科学版（增刊2）：15-20.
王伟涛，王宝善，葛洪魁，等.2009.利用主动震源检测汶川地震余震引起的浅层波速变化[J].中国地震，25（3）：223-233.
杨微，王宝善，葛洪魁，等.2013.大容量气枪震源主动探测技术系统及试验研究[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.
摘要：为研究噪音对主动源信号走时变化提取精度的影响，采用对比数值实验的方法。按照常规流程，首先用反褶积去除震源激发时间的不确定性干扰，压缩子波以得到更加尖锐的震相，由于发现反褶积会降低信号的信噪比，于是提出了台差法来代替反褶积矫正气枪激发的不确定性，再结合干涉法提取走时变化，并定量对比了两者的走时变化提取精度，最后尝试用台差法提取主动源走时的日变化。结果表明，台差法可以有效矫正震源激发时刻的影响，且不降低信噪比，从而也能获得较高走时变化提取的精度。
关键词：主动源；走时变化；走时台差法；信噪比
中图分类号：P315.31文献标识码：A文章编号：1000-0666（2017）04-0595-10
0引言
地震波信号是照亮地球内部结构的一盏明灯。复杂的地震波信号是由断层的破裂过程和介质格林函数耦合而成。传统的地震学分为研究震源运动学和地球结构研究（Lay，Wallace，1995）。从地震信号研究震源S（t）或者介质结构G（t）就必须使两者解耦。因此，大多数震源或者介质的研究中，都必须假设另一者有先验的模型。
由于主动源具有精确的激发时刻、激发地点，为地下介质的研究提供了良好的可重复震源，从而使精细获取介质结构变化成为可能。Niu 等（2008）利用压电陶瓷作为主动源，发现在实验场附近2个小地震发生前，压电陶瓷产生地震波信号走时有显著的延迟。在物理上，由于地震发生之前都会有裂隙和孔隙压力的变化，这个变化在介质波速上有所体现（Simmons，1964；Yukutake，1988），当这种精细变化被主动源信号所捕获时，将有助于我们认识地震孕育的物理过程。压电陶瓷具有频率高、衰减快、传播距离十分有限的特点，不能在大区域的研究中作为广泛应用的主动源。因此，陈颙院士在其所主导建设的宾川地震信号发射台中，引入了4枪组合的气枪激发源作为主动源，从而克服了传播距离有限的问题（罗桂纯等，2006；王宝善等，2011；王伟涛等，2009；杨微等，2013）。
目前，宾川发射台气枪实验已顺利进行多年，积累了大量的实验数据。如何提取高精度的主动源信号走时变化是众多学者关心和研究的问题。笔者首先介绍目前提取走时变化常用的干涉法和处理流程，之后对提取走时变化流程中反褶积这一步骤提出可能存在放大噪音而影响走时变化提取精度的问题，并进行相应的数值实验来验证这一想法，进而提出一种可能替代反褶积的台差法来矫正震源的影响，并进行了相应的实证。
1干涉法提取走时变化简介
波速测量关键的问题在于走时的获取，对于相对走时变化在0.01%～1%量级的精细变化（Wang et al，2008；刘自凤等，2015）探测中，传统手动拾取震相到时的方法已经不再适用。由于主动源气枪激发的重复性，波形上也表现出一定的相似性，基于互相关提取相对走时变化的方法比较适用（刘自凤等，2015；王彬等，2012）。本文以Wang 等（2008）所提出的干涉法为基础，开展信号走时提取。
1.1互相关函数
对于两个能量有限、定义在[t1，t2]上的有限长地震实信号x（t）和y（t），相关函数定义为：Rxy（τ）=∫t2t1x（t）y（t-τ）dt（1）其中：τ为时间延迟。物理意义上，相关函数表征了将一个信号相对另一个信号移动之后的相似性，为了方便比较，引入无量纲归一化的相关系数ρxy：ρxy（τ）=Rxy（τ）∫t2t1x2（t）∫t2t1y2（t）（2）互相关系数在[-1，1]间变化，互相关系数值越大说明2个信号越相像。2个相同的信号在τ=0时，互相关系数取1。
1.2干涉法
这里的干涉法特指为气枪设计的一种基于互相关提取走时变化的方法。在气枪主动源的走时变化探测过程中，存在以下2个问题，必须设计与之相应、适合气枪主动源提取走时变化的方法：
（1）气枪信号由很多震相组成，不同震相代表了地震波在地下介质中传播的不同路径，因此形成了对地下不同的采样空间。不同路径上的走时变化可能不同，所以有必要截取出一定长度的时窗，使之只包含一个特定的震相，以保证窗内所提取的走时变化是稳定而一致的。
地震研究40卷第4期黄亦磊等：台差法提高低信噪比下走时变化探测精度的研究（2）地震仪的采样率为100 Hz，对应的时间采样间隙0.01 s，然而介质速度变化引起的走时变化可能小于时间采样间隙，因此我们有必要对记录的地震信号进行有效插值，将走时测量精度大幅提高至0.001 s以上。
干涉法的原理如图1所示。首先我们需要得到一个参考波形，将参考波形和待测信号分成若干个小时窗，每个小时窗可以左右移动一定的步长，将移动过程中最大的互相关系数（图1c）对应的时间移动定为该时窗的相对时移（图1d）。每一个时窗都对应一个时间延迟值，我们需要選定一个波形显著、互相关系数较大且平稳的时间段，取这段时间中波形的互相关系数（图1c）作为走时变化的权重，加权求平均值，作为该道信号的走时变化。
从时移量得到速度变化有很简单的数学关系（Schaff，Beroza，2004）。由于信号的传播距离s（台源距）可以确定，假设介质速度原来是v，对应的走时为t，速度变化为dv，走时变化为dt，则有：
s=v·t=（v+dv）（t+dt）（3）忽略二阶小量dv·dt，则有：dvv=-dtt（4）从干涉法的实现过程可以看出，信号走时变化提取或多或少会受到信号处理流程的影响。现列举如下，以提示读者之后研究的注意：
（1）地震信号的前期处理。包括去均值、去除线性趋势、滤波和反褶积。其中去均值、去除线性趋势、滤波都是提高信号信噪比的有效手段。而反褶积则能有效矫正震源发震时刻不确定性的影响，并压缩子波得到更加尖锐的震相序列，但同时会降低信噪比。
（2）干涉法中时窗大小选取以及滑动步长选取。时窗大小以及滑动步长的选取带有经验性和探索性，也依据不同研究问题的目标而定。
（3）选取哪一段波形对应的走时变化做平均得到最后的走时变化。
如何选取时窗滑动步长、滑动步长选多大以及最后算平均走时变化的时窗，具有经验性和探索性，这是比较主观的影响因素，并非本文所要研究的问题。在之后的研究中，我们尽可能控制使窗口长度以及滑动步长等主观因素一致。主要定量研究噪声及滤波和反褶积等步骤对走时变化提取精度的影响，并尝试提出一种保持原信号信噪比，但又能一定程度上矫正震源时间不确定性的台差法，并对这种方法的有效性开展了实证研究。
2干涉法提取走时变化可能影响因素探究
影响一个物理量准确测量的原因有两种：一是客观存在的误差和噪音；二是主观人为因素。
仪器的误差无法消除，自然的噪音也客观存在，因此要在提取走时的过程中，尽量降低噪音，并使用合适的信号处理流程。笔者将重点研究在数据处理流程中可能会影响走时变化提取的因素。
刘自凤等（2015）提出了宾川主动源气枪定时变化提取的一个流程，如图2所示。这个流程是目前大多数走时变化提取所遵循的流程，这个流程看起来是比较合理的，但是刘自凤等（2015）并没有做数值方面的测试，也没有做主要步骤对提升结果的可靠性的分析。接下来，笔者将从理论测试出发，研究在不同的噪音水平和处理流程中（主要是滤波和反褶积）对走时变化提取结果的影响。
2.1噪音背景下反褶积对走时变化提取的影响
在与地震信号同频率的噪音存在的背景下（滤波无法消除噪音），实际地震信号对震源时间函数反褶积，噪音的影响可能会被放大。为了研究这一问题，需如下数值实验。
选择2013年前80天53281台（位置如图3所示）数据作为研究对象，将80天数据叠加，如图4a 所示。并把参考台CKT0（位置如图3所示）相同时间段的数据叠加作为震源时间函数（图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组实验数据，总结如下：
（1）组内实验结果对比，可以发现噪音越高，对走时变化提取影响越大。符合常规预期。
（2）滤波可以有效提高信噪比，降低走时变化提取的不确定性。符合常规预期。
（3）反褶积降低信噪比，反褶积之后提取走时变化的结果都差于不做反褶积，反褶积实质上起到了一个放大噪音的效果。
在此基础上，我们以53281台80天叠加信号为基础，添加噪音，使其信噪比从-10到15连续变化，对每一个信噪比下的信号滤波之后，做反褶积（尝试0.000 1，0.001和0.01的反褶积waterlevel水平）和不反褶积（只滤波）处理，之后提取走时变化，仍然以原始数据叠加而未添加噪音的信号作为参考模板。取每条信号的信号段（8～12s）提取的走时变化的平均值，作为该信号的走时变化。结果如图7所示，可以看出，所有反褶积处理之后的信号，走时变化提取结果都比只滤波处理的结果差，而且，随着反褶积waterlevel水平的增加，反褶积之后提取的走时变化结果也会更加好且稳定。
综上所述，我们在图7中定量展示了在存在不同水平噪音时，反褶积这一处理流程可能带来的走时提取误差。反褶积是一把双刃剑，它可以消除震源的影响，压缩子波，但是也会降低信号的信噪比，给走时变化提取带来不确定性。且高噪音水平要求反褶积时用更高的反褶积水准值。
3替代方案——台差法及其有效性实证
3.1反褶积的可能替代方法——台差法
提出用反褶积处理信号的初衷之一，是矫正气枪源的不稳定影响，比如激发时刻不同步以及水位变化等。反褶积可以去除震源的影响，留下介质的格林函数。但是，上述的数值实验说明，在存在噪音的情况下，反褶积会降低信噪比，从而导致变化提取不准确。另外，由于参考台得到的信号并非真实的震源时间函数，可能还需考虑方位带来的辐射花样问题，有一定的局限性。
因此，我们提出用台差法来代替反褶积以消除震源的影响，从而避免由于信噪比低而引起提取走时变化的精度问题。台差法的想法比较简单：震源如果存在激发时刻不准确或水位变动的影响，这种走时变化是一种整体时移，在每一个台站的记录都是近乎相同的，因此利用非常靠近震源的参考台进行自身干涉而提取到的走时变化，可以认为就是源的变化所导致的，然后再用其他台站提取的走时变化结果减去参考台提取的走时变化结果，这样便可矫正源的影响。总结台差法的步骤如下：（1）对信号去均值，去除线性趨势，滤波；（2）参考台用干涉法提取走时变化；（3）待研究台站用干涉法提取走时变化；（4）待研究台站走时变化减去参考台走时变化得到最终的走时变化。
为了验证台差法的正确性，我们做了如下实验：选取2013年前80天的53281台和CKT0台数据，分2组做对比实验：第一组用53281台单枪数据对CKT0同一枪数据做反褶积（反褶积水准值0.001），然后单天数据叠加，提取走时变化（走时变化用8.5～10.5 s时段走时变化平均）；第二组，对CKT0台和53281台数据单天叠加，然后分别用干涉法提取CKT0台和53281台走时变化，然后用台差法（53281台走时变化用9～11 s时段走时变化平均，CKT0台走时变化用0.5～1.5 s时段的走时变化平均）处理数据。上述2组实验在提取走时变化时的模板均为2台在这80天内的所有数据的叠加，结果如图8所示。可以发现，CKT0台和53281台提取的走时变化有很好的相关性，说明源的误差确实在这两个台都有所体现，相减之后，得到的走时变化非常平缓，在这个阶段内，并无大的环境因素改变，我们有理由认为这段时间的介质速度基本不变。对比反褶积提取的走时变化结果，可以发现台差法提取的走时变化结果比反褶积提取的走时变化结果更加稳定。
3.2实际资料处理
为了进一步验证本文提出的台差法的有效性，本文研究了气压以及地震可能对介质速度的影响（Spane，2002）。首先，选取了2015年11月20—28日的集中实验阶段作为研究对象，在此阶段，气枪激发比较密集，且有气压记录，适合研究气压对走时变化的影响。首先选择53264台、53274台、53281台和53266台作为研究对象（图3），
首先用台差法对单枪数据提取走时变化，进行加权平均的走时变化的时间范围分别为3～5 s、3～5 s、8～10 s和4～6 s，模板为9日内所有数据的叠加，4个台的走时变化和气压变化如图9所示。可以发现，在误差范围内，这4个台的走时变化都与气压有非常一致的变化，但是走时变化似乎领先气压变化一些相位。为了进一步观察他们的相关关系，本文细化了实验，选取53274和53281台作为研究对象，将这两个台在集中实验期间，每天同一小时、同一分钟的信号叠加，然后提取走时变化，进行加权平均的走时变化的时间范围仍为3～5 s和8～10 s，模板同样为9日内所有数据的叠加，得到的将是一天之内的走时变化。从图10可以看出，这两个台的走时变化随着气压都有正向相关的趋势，但是走时变化都领气压变化一些相位。
除此之外还计算了一天之内气压对于走时变化的具体数值影响，从走时变化结果中取出最大值和最小值做差，得到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）。
实验室进行的岩石实验表明，气压的变化会改变岩石内部裂隙中物质的密度，从而改变地震波在岩石中的传播速度，使波的走时发生改变（Nur，1971；Walsh，1965），但实际地球内部的岩体结构更加的复杂，气压与走时变化之间的具体关系仍待进一步探究。另外，除了气压对走时变化有影响之外，固体潮和温度（Leary et al，1979）等因素也可能通过影响介质速度而使走时变化产生波动，目前我们还没有进行这方面的探究，这也是接下来的研究内容之一。
4讨论
研究地震孕育过程中介质可能发生的变化，对于我们认识、理解地震发生过程具有重要意义。如果能有效地提取出介质的变化，并研究清楚这些变化与地震之间的关系，则我们在地震前兆的认识上将更进一步。
反褶积理论上能去除源的干扰。但实际上，由于气枪水下激发产生的地震波的机理尚未研究透彻，只能用距离气枪西面50 m处的一个参考台来近似源，在现有的条件下是一个不错的选择。另在一个问题，即参考台虽然距离气枪很近，但是受到震源辐射花样方位角的影响，因此方位角偏离参考台比较大的台站，用参考台信号反褶积的影响程度仍有待考究。
至于台差法能否准确地提取出介质的走时变化，有待进一步的验证。台差法目前来看有一个缺点，就是不能很好地将参考台信号对应到待研究的台站。由于参考台P波、S波混叠在一起，而在待研究台站可以明确选择某一震相，这样做差可能难以解释他们的对应关系，但是震源的不确定性在各个震相应该都有体现，所以我们相信可以部分消除源的影响。台差法的第二个问题在于射线路径的问题，台差法可以明确的矫正气枪的触发误差、水位变化带来的误差，但是对于不同的台站，射线路径不同，在从震源到参考台部分的走时，我们没有考虑不同台站是不一样的，这一部分误差不能保证完全消除，有待进一步研究。但是鉴于CKT0台离源非常近，我们认为影响不是很大。
进一步，我们将研究走时变化提前于气压变化的原因，并将用理论计算得出固体潮的理论值进行对比。另外，下一步将回溯地震，观察在其发生前是否有走时变化的突变。
5结论
在目前无法理论推导震源函数的情况下，前人想到用离震源非常近的参考台的垂向分量作为震源时间函数，来进行反褶积，在高信噪比的情况下，是一种可行之举。但是我们通过实践发现在存在较高噪音的情况下反褶积可能导致测量误差加大，使得走时变化的提取不准确，进而设计了理论实验和实际数据的实验，发现反褶积这个过程在走时提取中确实会降低信号的信噪比，从而影响走时变化的精确提取。然后，我们定量研究了噪音水平在走时变化提取中的影响。结果表明，反褶积只有在高信噪比的情况下（8以上）才不会影响走时变化提取的精度。而且原始信号的信噪比越低时，反褶积可能将信号的噪音放大越严重。所以当原始信号的信噪比相对较低时，应该提高反褶积的水准值。实验发现，不進行反褶积的信号与反褶积之后的信号提取走时变化相比，在信号信噪比低于8时，未反褶积的信号提取的走时变化更加准确和平稳，且对噪音的抵抗能力更强。
为了解决这个问题，笔者提出了一种可能的解决方法，也就是用参考台自身信号的走时变化来矫正远台的走时变化，这种方法，从目前的结果来看，也不失为一种矫正源效应的方法。另外，进一步进行了该方法的实证，结果发现，气压变化与气枪走时变化有很强的相关性，但是气枪走时变化总是领先气压变化，猜想可能和气温、固体潮等因素相关，这个负载的物理问题有待进一步研究。
参考文献：
刘自凤，苏有锦，王宝善，等.2015.宾川主动源地震波走时变化分析方法研究[J].地震研究，38（4）：591-597.
罗桂纯，王宝善，葛洪魁，等.2006.气枪震源在地球深部结构探测中的应用研究进展[J].地球物理学进展，21（2）：400-407.
王宝善，杨微，王伟涛，等.2011.宾川气枪信号发射台性能分析[C]//中国地球物理学会年会第二十七届年会论文集.
王彬，楊润海，王宝善，等.2012.地震波走时变化精确测量的实验研究[J].云南大学学报：自然科学版（增刊2）：15-20.
王伟涛，王宝善，葛洪魁，等.2009.利用主动震源检测汶川地震余震引起的浅层波速变化[J].中国地震，25（3）：223-233.
杨微，王宝善，葛洪魁，等.2013.大容量气枪震源主动探测技术系统及试验研究[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.
