0 引言

利用人工震源产生的地震波对区域尺度的地下结构进行主动探测，是一种研究地球内部结构的重要手段(陈颙等，2007a；王宝善等，2011；杨微等，2013)。自20世纪80年代以来，随着气枪技术的迅猛发展，气枪震源已经在海洋地球物理勘探领域取得了巨大的成功(Lutter et al，1999；McIntosh et al，2005；赵明辉等，2004；丘学林等，2007)。近10年来，陆上水库中的气枪激发实验也广泛开展，已有的研究发现，大容量气枪具有激发能量强、可重复性好、安全性高和绿色环保等特点，应用前景广阔(陈颙等，2007b；林建民等，2008、2010；唐杰等，2009；杨微等，2013)。

气枪震源到接收台站间的格林函数反映了研究区域地下结构的信息，在从地震信号中提取格林函数的过程中，需要去除气枪震源子波信息，这是一个反褶积的过程(许力生等，1996)。前人已经提出了很多种反褶积方法(Bell et al，1995；Clayton et al，1976；Gurrola et al，1995；Helmberger et al，1971；Kauppinen et al，1981；Kikuchi，1982；Sheehan et al，1995；黄绪德，1992；王西文等，2006；吴庆举等，2007)，这些方法在用于处理高信噪比的数据时效果都很好，但用于信噪比一般的数据时则有所不同(Ligorría et al，1999)，且这些方法的计算效率也不相同，上述问题都影响着气枪震源资料处理中反褶积方法的选择。

我们以云南宾川主动源资料为例，对频率域水准反褶积和时间域迭代反褶积2种常用的反褶积方法处理气枪资料的效果进行了对比，并对反褶积操作在气枪资料处理流程中的顺序问题进行了研究。

1 气枪实验简介

云南省宾川县大银甸水库(25°48′N，100°30′E)实验选用的是Bolt公司生产的1500LL型气枪，单枪容量为2000 in³。大银甸水库的总库容为4085万m³，激发点的水深约为15m，气枪的沉放深度约为10m。实验时采用4条气枪同时组合激发，1次激发释放的总能量相当于1次M_L0.7的天然地震，气枪源的优势频带为3~7Hz(杨微等，2013)。

在本研究中，我们使用了30个短周期地震仪记录的共493次激发(2013年2月25日~3月27日)的地震资料(图 1)。图 1中台站的偏移距范围为60m(离气枪震源最近的参考台)至90.93km。

2 方法与流程

地震仪器记录到的地震波数据可以表示为

$ d\left( t \right) = s\left( t \right) * g\left( t \right) * i\left( t \right) + n\left( t \right) $

(1)

式中，d(t)为观测数据；s(t)为震源时间函数；g(t)为气枪源到接收点之间的格林函数；i(t)为仪器响应；n(t)为噪声；“*”表示时间域的褶积。由式(1) 可以获取只与地球内部结构有关的格林函数g(t)。在实际数据处理过程中，忽略噪声n(t)的影响，先将仪器响应i(t)从观测数据d(t)中去除，得到

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

(2)

由于参考台距气枪源非常近(60m)，可以将参考台的记录近似看作气枪子波信号即震源时间函数s(t)，这样就可以从u(t)中反褶积掉s(t)，最终得到气枪源到接收点之间的格林函数g(t)。

2.1 频率域水准反褶积

频率域水准反褶积是一种在频率域内直接相除的反褶积方法。将式(2) 作傅里叶变换，可得其在频率域的表达式

$ U\left( \omega \right) = S\left( \omega \right)G\left( \omega \right) $

(3)

式中，U(ω)为u(t)的频谱；S(ω)为s(t)的频谱；G(ω)为g(t)的频谱。将式(3) 变换后可得

$ G\left( \omega \right) = \frac{{U\left( \omega \right)}}{{S\left( \omega \right)}} $

(4)

在实际应用中，这种频率域内直接相除的做法常常会由于分母的频谱出现零点或者极小的值使得商的频谱发生畸变。为此，Helmberger等(1971)提出了频率域水准反褶积方法，将式(4) 改进为

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

(5)

式中，S*(ω)为S(ω)的复共轭；max为求最大值运算符；c为水准比例因子。该方法的核心在于，先根据S(ω)的频谱特征选取一个合理的水准阈值，然后将分母频谱中振幅低于该阈值的提升至阈值，振幅高于阈值的则保持不变。这样既保证了反褶积的稳定性，又不会影响原始数据的频带峰值。由于反褶积经常会引入高频噪声(Ammon et al，1993；Velasco et al，1995)，所以得到G(ω)之后，还需要对其进行低通滤波，然后再将其反傅里叶变换回时间域，最终得到气枪震源到台站之间的格林函数g(t)。

2.2 时间域迭代反褶积

时间域迭代反褶积是一种通过在时间域内迭代构建预测信号并与观测信号作最小二乘拟合的途径来获取反褶积结果的方法。该方法最初由Kikuchi等(1982)提出，并被用来估算大地震的震源时间函数，后来也被Ligorría等(1999)应用到接收函数的求解中，其数学推导过程已由Kikuchi等(1982)详尽阐述，在此仅简要介绍操作流程：

(1) 将s(t)和u(t)作互相关扫描，得到互相关函数绝对值的最大值对应的时间延迟t₁，并通过求解Kikuchi等(1982)列出的一个关于最小二乘的方程可得幅值m₁，这样对u(t)中振幅最大的信号的预测就可表示为

$ u_{_{{\rm{pre}}}}^{^1}\left( t \right) = s\left( t \right) * {m_1}\delta (t - {t_1}) $

(6)

(2) 用原始信号u(t)减去预测信号u_pre¹(t)得到新的原始信号u¹(t)；

(3) 迭代执行步骤(1)、(2) 直至达到预设的最大迭代次数或者新的原始信号的总能量不再大于原始信号总能量的某个预设的比例。假设一共执行了N次迭代，则原始信号u(t)可以近似表示为

$u\left( t \right) \approx {u_{{\rm{pre}}}}\left( t \right) = \sum\limits_{i = 1}^N {} u_{_{{\rm{pre}}}}^{^i}\left( t \right) = s\left( t \right) * \sum\limits_{i = 1}^N {{m_1}} \delta (t - {t_1}) $

(7)

结合式(2) 可得

$ g\left( t \right) \approx \sum\limits_{i = 1}^N {{m_i}} \delta (t - {t_i}) $

(8)

这样就可通过时间域迭代反褶积方法近似获得气枪震源到接收点之间的格林函数g(t)。另外，与频率域水准反褶积类似，在实际操作中需要对信号进行低通滤波以消除高频噪声的影响。

2.3 数据处理流程

除了反褶积方法的选择，各种数据处理方法在整体流程中的排序也是气枪震源资料处理中会影响最终格林函数质量的重要问题。为此，首先需要对气枪震源的原始地震波形资料进行截取、去均、去势和去仪器响应等预处理；然后为了便于SH波、P-SV波的解耦，将经过预处理的波形资料从U(垂直)、N(北)、E(东)3个分量旋转到Z(垂直)、R(径向)、T(切向)3个分量；除此之外，还需要将不同炮的波形资料进行线性叠加以压制随机噪声。

基于不同的假设，可以选择先反褶积再叠加，也可以选择先叠加再反褶积。先反褶积再叠加基于的假设是研究区域地下结构在叠加时间范围内没有发生变化；先叠加再反褶积基于的假设是研究区域地下结构在叠加时间范围内没有发生变化且气枪震源在叠加时间范围内也没有发生变化。从静态地震成像的角度出发，可以忽略研究区域地下结构随时间的变化；在水库水位基本不变的前提下，大容量气枪震源具有极高的可重复性(唐杰等，2009；杨微等，2013)，实验期间水库水位基本保持稳定，因此也可以忽略气枪震源随时间的变化。在先反褶积的过程中，我们采用同一次激发观测到的台站信号和参考台信号做反褶积；在先叠加的过程中，选取数目相同且具有一一对应关系的台站信号和参考台信号分别做线性叠加，即每个叠加后的台站信号都有一个与之对应的叠加后的参考台信号。

反褶积是一种稳定性较差的数据处理方式，所以在反褶积之前要对数据进行带通滤波，以减少无关频带的干扰，并且在反褶积之后也需要对结果进行带通滤波。还应注意的是，如果滤波前数据的两端不为0，则需对其两端做尖灭处理。综上所述，我们设计了4种气枪震源资料处理流程(图 2)，并对比了它们的优劣。

3 结果 3.1 数据频谱特征

基于前述提及的2点假设，可将多炮叠加后的数据看作信噪比较高的单炮数据，以这样的单炮数据为例，我们绘制了参考台和53281台(图 1)去除仪器响应后数据的频谱特征图(图 3)，2个台的波形数据如图 4(a)、4(b)所示。根据气枪源的优势频带为3~7Hz(杨微等，2013)和数据的频谱特征，将资料处理流程中第1次带通滤波的频带设置为2~8Hz(图 2)，第2次的设置为2.5~5.0Hz(图 2)，最终得到气枪信号优势频率的波形资料。

3.2 反褶积结果

对于经过叠加和带通滤波之后的数据，我们分别用时间域迭代反褶积方法和频率域水准反褶积方法处理，最后再将结果滤波到信号的核心频带2.5~5.0Hz(图 4、5)。用时间域迭代反褶积时，设定当迭代次数≥100或数据残差小于原始信号的千分之一时终止迭代，最终反褶积结果如图 4(c)所示；用频率域水准反褶积时，设定的水准比例因子c为0.0001。将2种方法得到的反褶积结果再与参考台信号做褶积得到预测信号(图 4(d))，将预测信号和原始信号(图 4(a))做互相关，得到的互相关系数均大于0.99。采用同样的参数设置，也获得了先时间域反褶积后叠加和先频率域反褶积后叠加的结果。

3.3 4种资料处理流程的结果

以53281台站(偏移距约21km)Z分量数据为例，我们得到了用4种流程处理气枪震源资料的结果(图 5)及其信噪比。本文采用的是振幅信噪比，即有效信号的最大振幅与背景噪声振幅的均方根之比(林建民等，2008)。以53281台站为例，设定的背景噪声时窗为0~2.5s，P波时窗为2.5~5.0s，S波时窗为5~10s。最终得到的先叠加后时间域反褶积结果的P波信噪比为515.03，S波信噪比为1206.64；先叠加后频率域反褶积结果的P波信噪比为33.70，S波信噪比为66.63；先时间域反褶积后叠加结果的P波信噪比为48.99，S波信噪比为116.79；先频率域反褶积后叠加结果的P波信噪比为36.26，S波信噪比为71.68。由此可见，先叠加后时间域反褶积的方法可以得到信噪比最高的结果，且用其它分量数据进行对比的结果与用Z分量的基本一致。

对比4种处理流程的相对计算效率发现，先叠加后反褶积的计算效率要高于先反褶积后叠加，频率域水准反褶积的计算效率要高于时间域迭代反褶积。

3.4 走时剖面

对所有台站的数据都采用相同的处理流程(先叠加后时间域反褶积)，将最终得到的气枪震源到台站之间的格林函数按偏移距排列，得到这些台站的走时剖面(图 6)。从先叠加后时间域迭代反褶积的Z分量走时剖面图(图 6(a))中可以看到清晰的P波初至信息，S波等后续震相在不同台站的表现有所不同；在T分量走时剖面图(图 6(c))中P波能量较弱，S波等后续震相的能量较强。

4 讨论 4.1 反褶积方法

在选取合适的信号频带的条件下，时间域迭代反褶积结果的信噪比远高于频率域水准反褶积。由于时间域反褶积每次迭代时都拾取振幅最大的信号，所以对于信噪比大于1的情况，时间域反褶积具有自动压制噪声的功效，可使得P波初至信息更加清晰，频率域水准反褶积则没有这个优势。并且在频率域反褶积时，水准比例因子的选取也是一个难点，该因子越大则反褶积结果越稳定，但失真度也越大。因此，水准比例因子需要根据数据质量进行多次尝试来进行选择。

但是由于时间域反褶积需要进行多次迭代，这直接导致该方法的耗时显著高于频率域水准反褶积，由时间域迭代反褶积获得的格林函数是由一些(本文中N=100) 经过平移缩放的δ(t)组合而成的(式(8))，因而采用该方法得到的频谱与原始信号的频谱可能差别较大，引入了很多其它频带的噪声。另外，再考虑到远处台站的高频衰减效应，因而，最后我们在对结果做滤波时需要根据原始信号的频谱特征选择1个相对较窄的核心频带范围，若频带过宽，则可能会出现反褶积后波形畸变的情况。

4.2 资料处理流程

从提高信噪比的角度出发，特别是对时间域反褶积方法而言，先叠加后反褶积的处理流程更具优势。这是因为反褶积是一种稳定性较差的操作，叠加之前的信号噪声较强，对反褶积的结果造成了很大的干扰，并且由于2种反褶积方法中一些重要的参数(如水准比例因子c和迭代次数N)是相对值，所以并不能通过后期叠加而将这种干扰完全消除，这使得叠加后的信号噪声较弱，以至于在时间域反褶积的前N次迭代过程中，有可能并没有对这些噪声信号作预测(式(8))，因而导致这些噪声被自动忽略，进一步提高了格林函数的信噪比。

但是先做叠加需要同时满足研究区域地下结构在叠加时间范围内没有发生变化及气枪震源在叠加时间范围内也没有发生变化等2点假设，这导致这种资料处理流程的应用范围具有一定的局限性。另外，在原始资料质量非常好或非常差的情况下，这种处理流程的优势均不明显，这是因为在这2种情况下，资料处理的流程问题已经降为次要问题。通常情况下，为了通过气枪信号得到地下介质速度随时间的变化信息，可以采用一定时间段内(例如几天)的信号先叠加再反褶积的处理方法，从而保证相对较高的信噪比。

4.3 走时剖面

Z分量走时剖面中的P波初至信息较为清晰，S波等后续震相的能量较高，但无法清晰分辨。由于S波在水中无法传播，所以剖面中的S波信号不是气枪直接激发产生的，可能是P波在水库边界上转换而来的(杨微等，2013)，也可能是由于水库地形不对称引起的。在一维均一各向同性介质中，T分量不应该含有P波成分，但在T分量走时剖面图中，还是可以看到个别台站的P波震相，这可能是由于三维结构的复杂性所致，并且这种三维结构的复杂性也可能是造成剖面中偏移距相近的台站之间波形存在较大差异的原因。所以，如何正确识别气枪震源震相仍然是关键的科学问题之一。

5 结论

通过对云南宾川大容量气枪震源资料处理方法的研究，得出以下结论：

(1) 在一般的气枪资料处理过程中，时间域迭代反褶积方法在结果的信噪比方面优于频率域水准反褶积，但是计算效率相对较低。

(2) 在满足一定假设的情况下，先叠加后反褶积的处理流程在结果的信噪比、计算效率等方面都优于先反褶积后叠加的处理流程。

综上所述，我们提出了从气枪震源资料中提取气枪源到台站之间的格林函数的一般流程，即：① 截取、去均、去势、去仪器响应、旋转；② 叠加、带通滤波；③ 时间域或频率域反褶积；④ 带通滤波。

致谢: 感谢中国地震局地球物理研究所和云南省地震局等多家单位相关研究人员的大力支持，同时感谢审稿人对文章内容提出的宝贵意见。