0 引言

地震的发生是长期积累的应变在短期快速释放的过程，其孕育及破裂过程与断层区介质属性密切相关(Yang et al，2014；Yang，2015；Weng et al，2016)。研究地下介质属性及应力状态变化对深入了解地震的孕育、发生及震后恢复过程有重要作用(Li et al，2003)。然而，在孕震深度上测量应力的变化却十分困难。虽然应力状态可通过研究地震波的衰减、各向异性等获得，但对地震波速度变化的精确测量则是更为有效的研究方法(罗桂纯等，2008)。在地震发生前后，由于同震破坏和震后愈合效应，断层区的介质性质会随时间发生变化(Yang et al，2010、2011；Yang，2015)。而在孕震期间的应力积累及相应的裂隙变化，也会使断层区介质性质随时间产生变化。大量实验室研究表明，地壳岩石的地震波速度明显地随应力而变化，敏感度为10^-9~10^-8 Pa^-1(Birch，1960、1961)。这种依赖关系是由于作用在裂隙表面的应力变化导致微裂隙的张开/闭合所致(Nur，1971)。因此，应力的变化可以通过测量地震波速度的变化获得，而精确测量地下介质的波速变化则需要高质量的、可控的重复震源。由于天然地震的发生受地震重复性以及时空分布的制约(Wang et al，2012)，且不可控制，因此通常使用主动震源来监测地震波速度的变化，以期为孕震深度的岩石应力变化及岩石破裂愈合状态的研究提供信息(Li et al，2003)。目前，利用重复地震及主动源(压电陶瓷、气枪阵列、可控震源、炸药等)已经进行了许多探测地下速度结构变化的有意义的尝试(Schaffet al，2004；Chen et al，2007；Silver et al，2007；Niu et al，2008；Wang et al，2008)。

在所有人工震源中，气枪阵列以其可重复性高、探测距离远、能量高效、环保等优点被广泛应用于海洋及陆地的地下结构探测中(Chen et al，2007；Wang et al，2010；Xu et al，2014)。为了对地震频发的滇西地区，特别是对位于其中的红河、程海2个大断裂带地下的波速变化进行监测，2011年4月在云南省大理州宾川县大银甸水库建立了世界上第1个气枪地震信号发射台站(图 1)，该台可以连续不断地向地下发射传播距离远、重复性高的地震信号(Wang et al，2012；Chen et al，2014；杨微等，2013)。

相比于海洋中的气枪信号，陆地水库水深较浅，水域面积有限，且水位随季节变化明显，故气枪震源的信号并非完全重复(图 2)。因此，直接对各台站记录的气枪信号进行叠加得到的走时变化并不能反映真实的地下介质变化信息(图 2)。林建民等(2010)根据陆地水库中水听器记录到的气枪信号，对水库气枪的近场子波特性与沉放深度、工作压力等激发参数之间的关系进行了研究，但是没有给出更远处岸上地震台站记录到的气枪信号与这些激发参数之间的关系(陈蒙，2015)。远处岸上地震台站记录到的气枪信号受震源的影响程度与近场气枪信号受震源的影响间存在一定差别，因而有必要对震源变化的影响范围以及远场气枪信号的影响因素等进行分析，以了解该震源系统的特性。同时，对近场气枪信号与远场气枪信号分别选用合适的处理方法，以使气枪信号能够更准确地提供地下介质变化等信息。

1 主动源探测宾川实验站

宾川实验台站位于滇西南，该区域构造结构以印度板块与欧亚板块的碰撞为背景，地质构造十分复杂，是我国大陆内部地震活动最为强烈的地区之一。宾川实验站位于西南红河断裂、北部剑川断裂和东部程海断裂等围成的三角地块内(24.5°~28.0°N，98.5°~101.5°E)。地震信号发射台位于宾川县以西、距其7km的大银甸水库，水库汇水面积407.8km²。水库周围架设了40个流动地震台站(图 1)，使用了Guralp CMG -40T短周期流动地震计和Reftek 130B数据采集器，采样率为每秒100个数据点。

2013年1~12月间，在宾川气枪地震信号发射台站进行了1077次实验。每次由1个4条气枪组成的气枪阵列进行激发，激发压力为15MPa，枪阵的沉放深度为8m。大部分实验在14~22时进行，平均每天激发15次。由于降雨以及周围居民生活、灌溉用水，大银甸水库的水位会随季节发生明显的变化(>20m)。其中，水位较低的5~8月无法进行激发实验。

我们选用34个流动地震台站的数据进行分析，其中，最近的台站ckt距震源50m，最远的台站距震源151km。单次气枪激发信号记录可在20km的范围内清晰可见，在更远处的台站由于信噪比的降低则信号难以分辨(图 3(a))。但由于气枪信号重复性好，因此可通过对激发信号进行叠加来提高信噪比。本研究中我们采用了相位加权叠加方法(Schimmelet al，1997)。相较于线性叠加，相位加权叠加后的波形具有更高的信噪比，1000次叠加的波形在150km范围内非常清晰(图 3(b))。

2 气枪震源 2.1 气枪信号特征

首先对气枪地震信号进行去均值、去趋势处理。为了确定气枪信号的有效频率范围，对其信号利用S变换进行时频分析(Stockwell et al，1996)。信号h(t)的S变换为

$ {S_x}\left({\tau, f} \right) = \int_{_{ - \infty }}^{^\infty } {h\left(t \right)} \frac{{|f|}}{{\sqrt {2{\rm{\pi }}} }}{{\rm{e}}^{\frac{{ - {{\left({\tau - t} \right)}^2}{f^2}}}{2}}}{{\rm{e}}^{ - i2{\rm{\pi }}ft}}{\rm{d}}t $

(1)

其中，τ为时移参数；f为频率；窗函数w(t，f)为高斯窗函数，即

$ w\left({t, f} \right) = \frac{{|f|}}{{\sqrt {2{\rm{\pi }}} }}{{\rm{e}}^{\frac{{ - {t^2}{f^2}}}{2}}} $

(2)

气枪阵列单次激发(水库水位12.5m)产生的气枪原始地震信号及其时频图如图 4所示。首先以参考台(ckt)为例，距大容量非调制气枪阵列50m的ckt的原始记录可显示2~50Hz的能量，时间集中在0.5~1.0s，最大振幅约为12Hz。对于更远处的台站如53278台站(距震源2.3km)的时频分析结果显示，主要能量集中在2~5Hz，10~20Hz的能量幅度较小。随着震中距的增加，53281台站记录到的气枪地震信号主频为3~5Hz，5~10Hz有较小的P波能量。53257台站由于受不同频率噪声的影响，难以识别清晰震相(图 4(d))。以上分析显示，气枪信号主要能量集中在3~5Hz，因此，可在此频带采取带通滤波。

气枪激发的地震信号主要由2部分组成，即高压气体瞬间释放产生的压力脉冲和气泡振荡产生的气泡脉冲(林建民等，2010)。此外，在压力脉冲与气泡脉冲形成的同时，也会产生这些脉冲的虚反射，表现为负脉冲(陈浩林等，2008)。压力脉冲的频率较高，气泡脉冲的频率较低。如在参考台ckt，其压力脉冲主频为12Hz，持续约0.3s；气泡脉冲主频为4Hz，持续时间约1s(图 4(a))。在海洋石油勘探中为了提高探测的分辨率，主要使用压力脉冲，气泡脉冲是需要压制的噪声。但对于深部探测，为了提高探测的距离和深度，主要使用低频气泡脉冲。为此，下面我们主要对气枪阵列产生的气泡脉冲进行分析。

2.2 震源聚类与频率域反褶积

大银甸水库的水位由于季节性灌溉而变化，冬季水位可达20m，而夏季不足8m(图 5(a))。水位变化会导致气枪震源函数的变化，致使气枪震源波形并非完美重复(图 2、图 5(c))。我们对所有参考台的波形利用S变换进行时频分析后发现，水位变化会导致信号主频发生变化，随着水库水位的增加，气枪信号的主频也增大(图 5(b))。水位的变化也导致震源的波形不重复，其波形变化与水位变化相一致，水位较高时，气泡振荡周期变短；水位较低时，气泡振荡周期变长(图 5(c))。

为了系统地研究水位变化引起的震源函数的变化，我们对参考台ckt的波形利用互相关进行聚类分析(Yang et al，2009；王伟涛等，2012)。在此，未对远处台站进行聚类分析，这是因为由于受信噪比的影响，远处台站单天叠加波形之间的相关系数可能变动很大，而且不同台站之间的相关系数阈值不同，聚类结果也未必一致。而参考台ckt距震源很近，可以作为近似源函数，故对波形的聚类分析非常合适。首先，拾取参考台站的初至P波到时t₀，互相关窗口长度为5s，即P波之前0.5s、之后4.5s。将互相关方法应用于参考台ckt在2013年记录到的数据后发现，单天波形的互相关系数大于0.995，因此，仅对每天的单次激发进行相位加权叠加。由于单天实验激发次数最少为14次，故对单天数据进行14次叠加，再对单天叠加数据进行波形互相关，最后根据相关系数将参考台的波形分为7组(图 5(a))。同组波形的互相关系数大于0.99，如第7~38、333~360天。第47、252、261、267、272天由于与其他组的相关系数较低，因此不归为任何一组(图 5(a))。不同天之间的相关系数明显受水位变化的影响，但水位变化并非唯一的影响因素。如第1组与第272天，平均水位差为0.09m，相关系数为0.89，小于第272天与第2组的相关系数0.928(平均水位差为2.58m)，也小于第272天与第283天的相关系数0.973(平均水位差为1.89m)。这可能是由于受震源本身的影响，或季节变化导致水库水位及地下介质变化的结果。

为了研究水位变化引起的波形变化是否仅集中在参考台或某一距离范围内，我们对所有台站的波形根据参考台聚类的结果进了相位加权叠加，并对叠加后的波形进行比较。由于第4组数据包含第112、122、128天共42个事件的数据，因此，为了保证不同组的叠加次数相同，对不同台站记录气枪信号的同组波形均进行42次叠加。第6组数据只被个别台站记录到，因此仅对比其他6组波形。对比结果显示，水位变化不仅影响近场气枪信号，而且明显影响远场气枪信号(图 6(a))。不同水位下激发的气枪信号初至P波到时相互吻合，但后续震相到时差较大，因此，若要利用不同水位下激发信号来探测地下介质变化，则有必要去除震源变化的影响。

为了去除震源的影响，我们利用参考台记录的信号来近似震源函数s(t)，将远台记录u(t)与震源函数之间进行频率域水准反褶积(Helmberger et al，1971)

$ u\left(t \right) = s\left(t \right) * g\left(t \right) $

(3)

$ G\left(\omega \right) = \frac{{U\left(\omega \right)S * \left(\omega \right)}}{{\{ S\left(\omega \right)S * \left(\omega \right), {\rm{max}}\{ S\left(\omega \right)S * \left(\omega \right)\} c\} }} $

(4)

其中，*为褶积算子；g(t)为震源到台站间的格林函数；G(ω)为g(t)的频谱；U(ω)为u(t)的频谱；S(ω)为s(t)的频谱；S*(ω)为S(ω)的复共轭；c为水准因子，用于控制振幅最小值，以保证反褶积的稳定性。

首先根据聚类结果对各组波形进行42次叠加，然后分别对不同的震源函数(参考台记录)进行反褶积(图 6(b))。选取水准因子c为0.001，进行反褶积的数据窗口为从P波到时至P波到时之后16s。结果发现，当震中距小于10km时，各组反褶积后的波形互相吻合。当震中距大于10km时，反褶积的结果无法很好匹配。

3 讨论与结论

如上所述，云南宾川气枪震源信号主频为3~5Hz，建议在数据处理中先采用此频段的带通滤波。由于目前利用参考台的记录对远处台站叠加信号进行反褶积无法完全消除震源函数变化的影响，因此建议，首先对震源聚类分析，再基于聚类分析的结果对同一组震源函数的信号在所有台站进行叠加以提高信噪比。

由于每一天的震源重复性很好，相关系数大于0.995，故可对每一天的气枪信号作叠加。以2013年第1次激发实验(2013年1月7日)的叠加信号作为参考波形，将所有叠加信号对第7天进行互相关。若远处台站的后续信号与参考信号的相关系数变化小于参考台ckt的相关系数变化，则由目前的数据和处理手段暂时无法推断远处台站下或射线路径上是否有介质变化；反之，则可能存在一定的地下介质变化或者信号接收到其他噪声的干扰。我们首先拾取各台站记录的初至P波到时t₀，互相关窗口时长为5s，即P波到时之前0.5s和之后4.5s，各台站选取的互相关窗口长度一致，主要对比P波部分。由于第128~252天水库水位过低，无法进行激发实验，且第252~360天的波形只被部分台站记录到，因此选取第7~128天的记录进行分析。

图 7(a)为不同震中距的台站在不同天记录的波形与第7天记录的波形互相关的结果。由图 7(a)可见，参考台ckt不同天的波形与第7天波形的互相关系数变化与水位变化相一致。第47天，水位下降1.77m，互相关系数从0.999降至0.928。第56~80天，每天都进行激发实验，数据密集，互相关系数出现锯齿状起伏，但由于水位平稳，互相关系数维持在0.92左右。至第112天，水位下降了7.32m，互相关系数降至0.539。第128~252天，水位过低，无法进行激发实验。第272天，水位升至与第7天水位差仅为0.09m，互相关系数升为0.900。随着震中距的增加，波形受水位变化的影响逐渐减小，互相关系数的变化小于参考台ckt。在这些台站下方或者射线路径上即使有地下介质的变化，目前的数据及叠加方法也无法区分该变化是水位变化所致，还是真正的地下介质变化。

但是我们发现有4个台站(53265、53274、53267、53280) 的互相关系数变化大于参考台站(图 7(b))。53265台站距气枪震源4km，第7~59天的波形与第7天波形的互相关系数均大于参考台，若参考台站的相关系数变化可作为水位影响的基准，则53265台站受水位的影响应小于参考台，因此，互相关系数变化量小于参考台。第59~128天，53265台站的互相关系数变化趋势与参考台的吻合，但变化量更大，表明在射线路径上或该台站附近可能存在介质变化。另外3个台站，即53274、53267、53280，第7~112天的互相关系数变化量均大于参考台，表明该时段内这些台站之下可能存在介质变化。

波形互相关系数的变化除了由震源变化、地下介质变化等引起之外，也可能由信噪比引起。为了排除信噪比对波形互相关的影响，我们对上述4个台站(53265、53274、53267、53280) 计算了信噪比。首先对信号进行Hilbert变换，基于该变换计算出包络信号，计算信噪比时，取信号的最大值与噪声水平的平均值。由于对数据进行互相关的窗口为P波到时之前0.5s和之后4.5s，因此在计算信噪比时，选取同样窗口的数据，噪声部分选取0s到P波到时之前0.3s。由于已对数据进行了3~5Hz滤波处理及14次相位加权叠加，因此信号的信噪比非常高，以参考台ckt为例，不同天数据的平均信噪比达65dB。这4个台站中距震源最远的53280台站，不同天数据的平均信噪比也达到37dB。为了分析信噪比对数据互相关的影响，我们对每天波形的信噪比进行归一化，并对这4个台站的互相关系数以及在不同天记录的信噪比(经归一化处理)进行了对比(图 8)。第7~128天，归一化后的信噪比变化总体较互相关系数的变化小。第7~65天，53265台站的归一化信噪比小于互相关系数，表明信噪比对互相关系数起主要的影响作用。第65天后，信噪比远小于互相关系数，说明影响互相关结果的因素除了信噪比以外还有地下介质的变化。对于53274、53267、53280台站，信噪比小于互相关系数，说明这3个台站也记录到了地下介质变化的信息。

值得注意的是，53267、53280台站与53265台站方位角近似，从震源到台站的路径一致(图 1)，但53267、53280台站从第112天起的互相关系数即大于参考台站，并且大于距离气枪震源更近的53265台站。这表明53265台站下方在第112天起可能存在浅层介质变化，但这种变化可能未被53267、53280台站记录到。第80~112天，这3个台站记录信号的互相关系数相近，但在震源正西方的53274台站信号的互相关系数小于另外3个台站，这可能是由于波形受水位变化的影响小，互相关系数大，无法识别出该时段的地下介质变化。

为了消除气枪震源的影响，通常采用叠加与反褶积的方法。此时叠加与反褶积的顺序就显得尤为重要。先叠加再反褶积基于震源在叠加时段内没有变化的假设，因此，需要保证对互相关系数足够大的波形进行叠加，以减少水位变化对震源的影响。通过对震源聚类分析我们认为，需要对同组波形进行叠加，然后进行反褶积，这样可以在保证信噪比的前提下尽可能少地引入水位的影响。本研究的实验中，当震中距小于10km时，各组反褶积后的波形互相吻合，说明反褶积方法能有效消除水库水位变化导致的近场气枪子波的变化。但是在震中距大于10km时，由于远场气枪信号受震源的影响较小，且气枪信号除了受震源(主要是水位)的影响，还受地下介质变化及信噪比的影响，先叠加再反褶积的效果不如对相近水位气枪信号直接进行叠加的效果。

通过对云南宾川气枪震源数据处理的研究，我们建议采用如下气枪震源数据的处理流程：① 预处理，主要包括去均值，去趋势，进行带通滤波(范围应包含气泡脉冲主频3~5Hz)。对于有时钟误差的台站，需要进行钟差校正。② 对震源波形利用互相关进行聚类分析，对同组的地震波形进行相位加权叠加以提高信噪比。

我们对宾川主动源实验场周围不同台站2013年记录的波形进行如上处理，对同组波形进行了叠加，与第1次激发实验(2013年1月7日)记录的叠加波形进行互相关后发现，4个台站(53265、53274、53267、53280) 的波形互相关系数变化量大于参考台站。通过计算不同台站的信噪比发现，这4个台站信噪比的变化小于其互相关系数的变化，表明台站下方或由震源至台站的路径上可能存在介质的变化。

补充资料：

本文提供了对参考台ckt进行聚类分析时对不同天叠加数据的互相关结果，见附图。

本文进行相位加权叠加、互相关、反褶积处理前的数据挑选：

② 53277台站有很严重的时间零飘，因此未使用该台站数据。

② 2013034163959噪声很严重，故删除所有台站的此事件。

③ 2013086073808、2013086073459、2013086073237走时存在问题，可选择删除或校正走时。

④ 由于从第325天开始，在21：00左右将4个气枪变为1个气枪进行激发，波形与4个气枪有很大差别，因此将这些使用1个气枪激发的事件删除，见下表：

2013325210959 2013325211959 2013325212459 2013325212959 2013325213459 2013325213959 2013325214459 2013325214959 2013325215459 2013325215959 2013325220459 2013325220959 2013325221459

2013338211459 2013338211959 2013338212459 2013338212959 2013338213459 2013338213804 2013338214959 2013338215459 2013338215959 2013338220459 2013338220959 2013338221459 2013338221959 2013338222459

2013345201959 2013345202459 2013345202959 2013345203459 2013345203959 2013345204459 2013345204959 2013345205459 2013345205959 2013345210459 2013345210959 2013345211459 2013345211959 2013345212459 2013345213459 2013345213959 2013345214459 2013345214959 2013345215459 2013345215959 2013345220459 2013345220959

2013353201959 2013353202459 2013353202959 2013353203459 2013353203959 2013353204459 2013353204959 2013353205459 2013353205959 2013353210459 2013353210959 2013353211459 2013353211959 2013353212459 2013353212959 2013353213459 2013353213959 2013353214459 2013353214959 2013353215459 2013353215959 2013353220459

2013360210959 2013360211459 2013360211959 2013360212459 2013360212959 2013360213459 2013360213959 2013360214459 2013360214959 2013360215459 2013360215959 2013360220459

以参考台ckt为例，第325、338、345、353、360天不同时刻事件的波形对比图：

致谢: 感谢两位审稿人提出的中肯修改建议。