零相位化和非线性叠加对提取气枪震源地震信号波形的影响研究
姚佳琪+曹文忠+叶泵+张伟
摘要：使用零相位化处理和非线性叠加方法提取低信噪比记录中的有效信号，通过定量评估结果波形的走时和振幅等信息的精准程度，探究此处理流程是否可提高提取气枪震源有效信号波形的准确性。对合成波形和云南宾川主动源实验中实际资料数据进行零相位化处理，选用相似性加权、时间域相位加权、时频域相位加权、改进的时频域相位加权4种非线性叠加方法提取有效信号。结果表明，基于上述处理方法可有效提高提取波形的准确性：（1）零相位处理可以提高波形数据的分辨率，使有效信号到时点变为波形峰值点，有利于在非线性叠加后拾取出准确的到时位置；（2）零相位化数据的非线性叠加结果卷积震源子波可基本恢复地震信号波形记录，可降低非线性叠加对有效信息的损失；（3）使用时间域相位加权叠加和相似性加权叠加方法可获得较好结果，但可能会压制低信噪比小振幅信号；基于时频域的相位加权类叠加方法对有效波形成分影响较大，但对小振幅信号保幅较好。
关键词：主动源；气枪震源；零相位化处理；非线性叠加；信号提取
中图分类号：P315.3文献标识码：A文章编号：1000-0666（2017）04-0581-14
0引言
人工震源技术可以克服传统天然震源时空分布不均的制约，为获得高时空分辨率的区域尺度地下结构和介质变化信息提供可能（陈颙，朱日祥，2005；王宝善等，2016）。与天然震源相比，人工震源单次激发能量较低、远距离传播易受噪声干扰，需要大量重复激发再从中提取有效信号（张军华，2011；杨微等，2013；倪宇东，2014；顾庙元等，2016；武安绪等，2016）。提取结果的可靠性将直接影响后续地震资料处理结果的准确性。因此，如何从低信噪比的重复激发地震记录中获取高质量有效信号，是人工震源实验中的重点研究问题之一。
非线性叠加类方法是从低信噪重复激发记录中提取有效信号最为常用的方法之一。此类方法利用了有效信号的一致性，从而实现保留有效随机信号、压制噪声的目的（Kennett，1987；Schimmel，Paulssen，1997；Schimmel，Gallart，2007；Li，Gao，2014）。相比线性叠加，非线性叠加可大幅提高信噪比，但会造成信噪比过低的信号发生波形畸变和相位偏移等（武安绪等，2016）。在主动源实验中，非线性叠加方法也被应用于大容量气枪震源、低频可控震源等地震资料有效信号提取中。武安绪等（2016）使用线性叠加、时间域相位加权叠加、时频域相位加权叠加方法对接收到的波形記录直接处理并定量评估提取结果的可靠性。经合成测试和实际数据处理后，认为加权叠加会影响有效信号提取的波形质量和时间精度。与大容量气枪震源不同，低频可控震源选用震动时间较长的扫频信号激发且激发能量更低（陶知非等，2011）。为了降低不同震相的地震信号重叠干扰并提高数据分辨率，需要对接收到的原始波形做零相位化处理，通常是将地震记录与震源子波做互相关，将地震记录转化为地下介质脉冲响应与震源信号自相关子波的卷积记录（熊翥，2008；倪宇东，2014）。顾庙元等（2016）使用线性叠加、非线性叠加、直接分离类多种方法对主动源实验中零相位化的低频可控震源数据进行处理，获得了较好的有效信号信息。本文将零相位化处理流程引入气枪震源资料有效信号提取中，通过定量评估提取结果波形的走时和振幅等信息的准确性，探究使用零相位化数据是否可降低非线性叠加对低信噪比记录中有效信号波形的影响。
本文首先使用相似性加权叠加、时间域相位加权叠加、时频域相位加权叠加、改进的时频域相位加权叠加4种非线性叠加方法提取未经零相位化处理的加噪合成低信噪比信号，提出了非线性叠加直接提取气枪震源信号所存在的问题。然后，引入零相位化处理流程，分别使用互相关和反褶积两种常用的零相位化处理方法对气枪震源子波零相位化，对比两种处理手段得到的零相位化子波形态特征来确定适合气枪震源资料的零相位化处理方法。最后，通过对合成地震剖面数据和云南宾川主动源实验实际数据的处理，探究并验证先对数据做零相位化处理再叠加的可行性及优势，为主动源实验的信号提取提供参考。
1非线性叠加
1.1方法原理
非线性叠加是在线性叠加基础上发展起来的，是基于记录中有效信息波形或相位等的一致性来计算各采样点的叠加权重以提高信噪比，又称加权叠加。此类叠加方法既可在时间域中实现，又可在时频域中实现。时间域非线性叠加可表示为：yws（t）=w（t）γ yls（t）（1）式中：yws（t）为时间域加权叠加结果；w（t）表示时间域的加权系数；γ为加权指数因子，用于控制加权强度，一般在1.5～2.5之间；yls（t）为线性叠加结果，即：yls（t）=1N∑Nj=1xj（t）（2）式中：N为叠加次数；xj（t）为第j次激发所接收的地震记录；yls（t）为线性叠加结果。而时频域非线性叠加利用了有效信号时频域分布的一致性，可表示为：yws（t）=T-1wγ（t，f）1N∑Nj=1T{xj }（t，f）（3）式中：T和T-1分别表示某种时频变换及其逆变换算子，通常使用Stockwell变换及其逆变换；w（t，f）为时频域下的加权系数。下面，介绍本文所用非线性叠加方法的具体原理。
1.1.1相似性加权叠加
相似性加权叠加是以波形相似性系数（Semblance）作为一致性测量准则而实现非线性叠加的（Kennett，1987）。所谓波形相似性系数是多道记录叠加后瞬时总能量与各单道记录瞬时能量总和之比，用于测量多道记录一致性程度的常用参数（Taner，Koehler，1969；Neidell，Taner，1971），可直接用于加权系数计算：w（t）=∑n/2i=-n/2wG（i）∑Nj=1sj（t+i）2N∑n/2i=-n/2wG（i）∑Nj=1s2j（t+i）（4）式中：wG（i）为高斯窗函数exp[-（τ-t）22k2]，这里k为高斯窗宽度因子，一般应大于震源子波宽度，合理选取可使加权系数光滑稳定；n为窗函数的个数。
1.1.2时间域相位加权叠加
时间域相位加权叠加是以波形记录的瞬时相位作为衡量有效信号一致性依据的叠加方法（Schimmel，Paulssen，1997），其加权系数计算公式如下：w（t）=1N∑Nj=1exp[iΦj（t）]（5）式中：Φj（t）第j次激发记录随时间变化的瞬时相位，可通过Hilbert变换得到。由接收到的地震记录x（t）和其Hilbert变换形式H[x（t）]可构建出解析信号：X（t）=x（t）+iH[x（t）]（6）式中：X（t）表示解析信号。再将解析信号用指数形式表达：X（t）=A（t）exp[iΦ（t）]（7）式中：A（t）为局部振幅；Φ（t）为瞬时相位，是用于判断有效信号一致性、计算加权系数的重要参数。
1.1.3时频域相位加权叠加
时频域相位加权叠加是时间域相位加权在时频域中的拓展，这种方法是以地震记录的时频域瞬时相位的一致性作为评估依据（Schimmel，Gallart，2007）。该方法是通过Stockwell变换将信号变换到时频域中（Stockwell et al，1996，1999），其加权系数的表达式为：w（τ，f）=1N∑Nj=1Sj（τ，f）exp（i2πfτ）Sj（τ，f）（8）式中：Sj（τ，f）表示信号经Stockwell变换的时频谱。Stockwell变换的表达式可写为：Sj（τ，f）=∫∞∫-∞xj（t） wG（t-τ，f）exp（-i2πfτ）dt（9）这里，高斯窗函数wG（t-τ，f）随频率变化，其表达式为：wG（t，f）=f·（2π）-1/2·exp（-t2f2/2）（10）在時频域中实现相位叠加后，再通过Stcokwell逆变换后可得到叠加后的波形记录，如（3）式所示。
1.1.4改进的时频域相位加权叠加
改进的时频域相位加权叠加起初是为解决CMP道集叠加中的时移问题（Li，Gao，2014）而提出的，这里将其用于提取低信噪比信号。此方法是在基于Stockwell变换的自适应滤波叠加基础上发展而来的（Pinnegar，Eaton，2003），它不仅以时频域的瞬时相位为衡量一致性依据，同时也考虑了时频域中信号局部振幅的一致性程度。与以上加权叠加方法不同，改进的时频域相位加权叠加需要计算每一道记录的加权系数，以第k次重复激发记录的权重系数为例：wk（τ，f）=Uk（τ，f）-Vk（τ，f）maxτ，f{Uk（τ，f）-Vk（τ，f）}（11）其中，Uk（τ，f）=1N-1∑Nj=1，j≠kSj（τ，f）maxτ，f{Sj（τ，f） }v（12）Vk（τ，f）=minτ，f{Uk（τ，f）}（13）第k次记录的加权系数是通过其他各记录叠加、归一化等计算得到的。然后对每条记录乘以相应的权重系数后叠加，最终叠加结果为：y（t）=S-11N∑Nk=1Sk（τ，f）wk（τ，f）（14）该方法在实现过程中进行了两次叠加，因此很大程度上压制了随机噪声，保留并突出有效信号，但计算成本较高。
1.2非线性叠加提取有效信号存在的问题
这里使用简单的合成波形来对比评估上述4种非线性叠加方法的信号提取效果。首先使用云南宾川实验CKT台记录到的气枪震源激发信号作为震源子波（图1）构建时长10 s、采样间隔0.01 s、有效信号到时3 s、最大绝对振幅为1的合成波形（图2a）。再向合成波形中加入标准差为1、均值为0的高斯随机噪声（图2b）得到了加噪波形记录。重复上述操作，构建出300道加噪记录。
对上述加噪波形带通滤波2～8 Hz后（图2c）做叠加，结果如图3a所示。由图可知，非线性叠加相比于线性叠加可有效压制噪声、提高信噪比，但会对波形等有效信息造成影响：（1）相似性加权叠加和时间域相位加权叠加方法在压制噪声的同时，也会压制有效信号中小振幅波形使得小振幅波形信息难以完好保留；（2）基于Stockwell变换的时频域相位加权叠加和改进的时频域相位加权叠加方法虽然可以提取出小振幅波形信息，但会造成相位偏移产生虚假波形；（3）由于震相到时为波形起跳点，基于此从结果中拾取到时，从拾取结果（图3中深灰色竖线所标）可看出震相的到时难以准确拾取；（4）线性叠加虽不会对信号波形造成影响，但随机噪声难以压制完全，且当叠加次数有限时仍难以从结果中拾取出准确的到时位置。
由以上结果可看出，直接对气枪震源记录非线性叠加会对波形造成影响，甚至可能导致初至起跳时间拾取错误。之所以初至时间会被非线性叠加破坏，是因为对于因果信号（最小相位信号），初至时间位于能量较弱的起跳点。如果地震震源采用零相位子波激发，则初至时间位于零相位子波的峰值，有可能在非线性叠加过程中保留下来。本文尝试先对气枪信号进行零相位化处理，再进行非线性叠加，考察非线性叠加对零相位化后的信号走时和振幅等波形信息的定量影响，进而对能否应用非线性叠加提取气枪震源激发的有效信号提供理论依据。
2零相位化处理
地震波的传播可以看作是一个线性系统，符合褶积模型：x（t）=r（t）*s（t）（15）式中：x（t）为台站接受到的地震记录；r（t）为地下介质的脉冲响应（天然地震学中称为格林函数）；s（t）为震源子波信号。由此可见，地震记录可近似认为是地层脉冲响应与震源子波卷积的结果，各震相实际是混合相位的（熊翥，2008）。对于接收到的记录，震相的到时位置为波形起跳点，但往往因噪声的影响而难以准确拾取走时（倪宇东等，2010）。零相位化处理可将混合相位波形转化为零相位子波，从而使震相到时点变为零相位子波的峰值点，易拾取有效信号信息。此外，零相位化处理可将震源子波压缩为近似δ函数的零相位子波，提高数据分辨率，避免了多震相相互叠加、干扰等。下面介绍互相关、反褶积2种常用的零相位化处理方法原理。
2.1方法原理
2.1.1互相关
互相关方法是把地震记录与震源子波作互相关，将记录中的震源子波信号转化为零相位的互相关子波。互相关方法常用于可控震源资料处理中（图4a-1、a-2）。地震记录与震源子波做互相关后可表示为：x（t）s（t）=r（t）*s（t）s（t）（16）式中：表示互相关算子。对上式简化可表示为：x（t）s（t）=r（t）*k（t）（17）式中：k（t）表示震源信号互相关子波，这是零相位化的子波，此时震相的起跳点变为峰值点。
2.1.2反褶积
反褶积方法是将地震记录中的震源子波去除掉，从而得到具有较高分辨率的地下介质脉冲响应。反褶积处理被广泛应用于地震学偏移处理等中（熊翥，2008；Shearer，2009），气枪震源资料处理中也常用此种手段提取传播路径的格林函数。反褶积的方法有很多，最直接的处理方法是频率域反褶积，即将式（15）作Fourier变换将其转换至频率域中：X（ω）=R（ω）S（ω）（18）式中：X（ω）为地震记录x（t）的频谱；R（ω）为地下介质脉冲响应r（t）的频谱；S（ω）为震源子波的频谱。将X（ω）与S（ω）相除可得到较高分辨率的地下介质信息：R（ω）=X（ω）S（ω）（19）因实际数据中混合有随机噪声，频率域直接相除常由于分母频谱出现极小值或零点，造成结果不稳定。为了避免此种现象，这里选用计算稳定且效率高的频率域水准反褶积方法（Helmberger et al，1971；翟秋实等，2016）进行后续基于反褶积方法的零相位化处理。
2.2气枪震源的零相位化处理方法
根据以上常用的零相位化处理方法，本文将探究适合气枪震源的零相位化方式。首先分别使用互相关和反褶积方法获得气枪震源的零相位子波波形，对比2个子波形态以选择适合气枪震源资料的零相位化处理方法。为了更好地对比评估，选择适合气枪震源的零相位化方法，这里同时对云南宾川实验气枪震源ckt台波形（图4b-1）和长江计划安徽实验中低频可控震源理论子波（图4a-1）做零相位化处理，然后對比相应零相位化子波波形特征。在使用主动源进行探测试验中，可控震源资料必须经零相位化处理后获得较高分辨率数据再进行后续处理（Brittle et al，2000；倪宇东，2014；顾庙元等，2016），而气枪震源较少使用此流程来提高数据质量。此外，可控震源由于扫频信号自身的特性，其互相关子波主瓣明显、旁瓣噪声低可近似为δ函数，且互相关方法实现简单且稳定（特别是对于低信噪比记录），所以对于此类数据最常用的零相位化方法是互相关方法。
经互相关后得到的震源自相关子波波形分别如图4a-2、b-2所示。使用反褶积方法进行零相位化处理后的子波波形如图4a-3、b-3所示。
从图4结果中可看出，气枪震源的自相关子波（图4b-2）不具备可控震源自相关子波（图4a-2）高清晰度的特征，且时长过长而难以达到压缩震源子波提高分辨率的目的，而利用反褶积得到的子波波形（图4b-3）具有很好的分辨率。基于此结果，本文认为对于气枪震源数据的零相位处理使用反褶积方法可获得更好的零相位化处理效果，后续使用频率域水准反褶积来实现气枪震源数据的零相位化。
为了探究零相位化方法是否可以降低非线性叠加对波形信息的影响，这里对1.2中的300道加噪合成数据进行基于反褶积方法的零相位化处理，然后再进行非线性叠加。单道加噪合成记录（图2c）的零相位化处理结果如图2d所示（经2～8 Hz滤波），零相位化后300次叠加结果如图3b所示。对图3b各结果设置2～4 s波形对比时窗定量评估提取波形的准确度，如图中文字所标。由结果可看出：（1）非线性叠加可以很好地压制随机噪声，并且较直接叠加结果（图3a）相比可更为准确地保留波形信息；（2）经零相位化处理后震相拾取位置由拾取起跳点变为拾取波形峰值点，经非线性叠加后峰值点位置不受影响，因此震相到时仍可准确拾取；（3）对比4种非线性叠加结果，可以看出相似性加权叠加和时间域相位加权叠加可较好保留零相位子波信息并压制子波旁瓣，而基于Stockwell变换的时频域相位加权和改进的时频域相位加权方法会导致子波波形相位变化而产生虚假波形，但波形峰值位置即有效信号的到时不受影响。由此可见，数据经零相位化后可降低非线性叠加对结果的影响，有利于保证后续处理的准确性。
3合成地震剖面数据处理测试
为了更准确地评估零相位化处理是否可降低非线性叠加对有效信号的影响，本文使用一维速度模型合成地震剖面数据来模拟大容量气枪震源信号处理过程。
3.1数据合成
首先在PREM全球平均速度模型（Dziewonski，Anderson，1981）的基础上（图5a），仿照真实地下介质情况构建一维速度结构用于合成仿真剖面数据，即在0～0.5 km近地表深度范围内加入了强衰减的低速层，衰减系数设置为QP=80、QS=40，并对加入了0.5～24.4 km地壳中设置速度逐渐变化的薄层，构建出的一维速度模型分布，如图5b所示。然后使用FK程序（Zhu，Rivera，2002）计算上述一维速度结构（图5b）的理论波形。气枪震源的震源机制复杂，这里将大容量气枪简化为作用于地表的垂直单力源，震源子波采用云南宾川主动源实验场ckt台接收单次激发记录（图1）。为了便于研究体波信息，这里从正演所得波形中去除面波成分，最终得到了该速度结构的理论地震剖面记录（图6a）。并使用Taup震相走时计算程序（Crotwell et al，1990）计算对应的Pg、Sg、Pn、Sn震相走时，如图6中所示。
然后，模仿真实噪声情况向理论地震剖面中加入随机噪声，构建合成地震剖面数据。本文通过测量震中距为5 km左右的台站接收的实际记录中信噪比来确定所加随机噪声强度。这里所定义的信噪比为有效信号时段内有效信号的最大振幅与背景噪声均方根振幅之比（林建民等，2008）。图7为云南宾川主动源实验中震中距为5.009 6 km的53260台大容量气枪震源单次激发记录，从记录中设置Pg波到时开始至10 s区间范围内的信号窗和时窗长度为5 s的噪声窗，分别如图中深灰色阴影和浅灰色阴影区域所示，信噪比为244.437 4。根据测量出的信噪比，向震中距为5 km处的理论波形中加入相应强度的高斯随机噪声。由于接收台站的噪声水平基本一致，因此对其他震中距的理论合成记录中加入与5 km台站相同标准差的高斯随机噪声，合成出的单次激发仿真地震剖面如图6c所示。重复上述加噪操作，构造出300次重复激发仿真数据。
3.2数据处理与分析
首先，对加噪合成数据进行去均、去势、去尖灭、2～8 Hz带通滤波等预处理。然后，对预处理后的数据直接叠加得到线性及非线性叠加结果，分别如图8a、c所示；而对于零相位化数据的叠加，则先使用频率域水准反褶积对预处理后数据作零相位化处理，然后进行2～8 Hz带通滤波处理再叠加，最后得到线性及非线性叠加结果分别如图8b、d所示。
这里抽出震中距分别为20、50、100 km合成数据的叠加结果（图9）以更好地观察提取有效信号效果。根据Taup震相走时计算程序计算出的震相到时，可以从结果中看出直接叠加和零相位化后叠加的波形特征。对于直接叠加走时点位置为波形起跳点，但非线性叠加后初至时间破坏；而对于零相位化后叠加的结果，走时点位置为该震相时段内波形出现的第一个振幅明显的峰值点处。因此本文认为：（1）低振幅的震相到时起跳点在叠加过程中振幅发生变化，导致很难从结果中拾取出准确到时；（2）而零相位化后数据的叠加结果可以避免上述问题，较好地保留有效信号信息；（3）将未零相位化数据的叠加结果和零相位化数据的叠加结果相结合可用于震相走时拾取。
4实际数据处理分析
将前文所用的先零相位化后非线性叠加处理流程应用于2013年云南宾川主动源实验大容量气枪震源地震资料实际数据处理中，验证零相位化处理是否可以提高提取信号的可靠程度。本文选取了该实验中的53280、53259、53270台（震中距分别为18.108 3 km、46.545 2 km、91.609 4 km）记录。对此实际数据进行了预处理、零相位化、叠加等流程处理，其中预处理流程中除去均、去势、去尖灭及带通滤波处理外，还需进行去仪器响应和分量旋转。根据合成测试中的零相位化子波波形和震相到时点位置波形等特征，本文将气枪震源波形记录直接线性叠加结果和零相位化后的波形记录非线性叠加结果相结合拾取震相到时，本文认为的到P、S波时位置分别如图10中灰色实、虚线所示。
5讨论
为定量评估各方法提取有效信号的效果，本文选用最大互相关系数R和时间延迟Td来评估提取波形的准确程度。最大互相关系数是所选取的波形对比时窗内的叠加结果与对应的理论波形的互相关系数最大值，时间延迟是互相关系数最大值所对应的时间偏移大小（王彬等，2012；肖卓，高原，2015；武安绪等，2016）。然后分别从图3a、b两组处理结果选取2～6 s和2～4 s的波形对比时窗，计算出的R、Td大小如图中文字所标。
5.1非线性叠加对各频段成分的影响
选取1.2和1.3中的非线性叠加结果，以时间域相位加权叠加结果为例分析非线性叠加结果的各频率成分的影响。从频率范围为2～8 Hz的直接叠加结果（图3a）和零相位化后叠加结果波形（图3b）中分离出2～4、4～6、6～8 Hz频段波形，如图11红线所示。将分离出的各频段成分与对应频段的理论信号做比较（图11文字标出最大互相关系数R的大小）结果显示，非线性叠加对4～6 Hz频率成分影响最小，零相位化数据叠加结果各频率成分保持良好。而大容量气枪震源主频在5 Hz左右，可以看出主频范围的信号成分保留较好。
5.2非线性叠加对振幅的影响
为了探究各叠加方法是否可以保留有效信号振幅，这里使用1.2中的合成波形（图2a）构建一组振幅变化的合成数据，最大振幅从1至0.1，间隔0.1变化，共10道数据，如图12a所示。再向数据中加入标准差为1的高斯随机噪声，产生300组数据，然后进行叠加。
从叠加结果的最大振幅变化情况（图12b、c）中可看出：（1）非线性叠加可以基本保持未零相位化和零相位化数据的振幅信息，相似性加权和时间域相位加权叠加振幅结果与线性叠加基本一致，基本吻合真实振幅；（2）而时频域相位加权叠加类方法虽基本拟合真实振幅，但是所得振幅偏小；（3）此外对于最大振幅为0.1的记录，零相位化处理后数据的相似性加权叠加和时间域相位加权叠加结果的振幅保持效果较基于时频域的相位加权叠加结果差，此测试结果与图9b中远震中距叠加结果相似。本文认为當信噪比过低时，相似性加权叠加和时间域相位加权叠加仅利用了时间域上的波形和相位信息，在压制噪声提高信噪比的同时，也将小振幅的有效信号信息压制而无法获得准确的波形信息；而时频域相位加权叠加充分利用了有效信号的时频分布特征，可从低信噪比记录中提取出有效信息。
5.3零相位化非线性叠加结果恢复波形能力的讨论由图3所标文字可知，对于气枪震源合成波形直接叠加，当叠加次数足够时，非线性叠加的提取效果较线性叠加差；而对于零相位化的气枪震源合成波形叠加，相似性加权叠加和时间域相位加权叠加的提取效果优于线性叠加方法；时频域相位加权和改进的时频域相位加权叠加会导致有效信号相位偏移而产生虚假波形，所提取的结果与真实波形差异较大。这里将零相位化数据叠加结果（图3b）卷积气枪震源子波波形（图2a）测试零相位化的叠加结果是否可用于恢复原始信号并提高提取结果的可靠性，结果如图13所示。设置2～6 s波形对比时窗与气枪震源子波波形（图3a）进行比较发现（图13文字所标），零相位化数据非线性叠加结果的最大互相关系数高于数据直接进行非线性叠加结果的，其中零相位化波形叠加结果中相似性加权和时间域相位加权叠加所恢复的气枪震源波形结果略优于波形直接线性叠加结果。这是由于波形经零相位化处理后原记录中波形的起跳点转化为了零相位子波的峰值点，峰值点的振幅较大，在叠加过程中不易受噪声干扰，因此可以保留较好的走时信息和地下介质的脉冲响应信息；此外，相似性加权叠加和时间域相位加权叠加可以一定程度上压制零相位子波的旁瓣噪声，从而进一步提高数据的分辨率，而基于时频域的相位加权类叠加方法由于Stockwell变换分辨率等问题造成相位偏移等现象，使得最终经逆变换至时间域后产生虚假波形而影响数据质量。
5.4以近台站记录近似震源子波对零相位化结果的影响由于气枪震源的震源机制复杂，尚未有很好的理论震源子波，目前是在距震源0.05 km处放置参考台站（CKT台），使用参考台接收记录作为震源子波做后续处理。下面探究采用近台站记录作为理论震源子波是否会影响零相位化结果的准确性。使用FK程序计算3.1合成测试中所用一维速度模型下震中距为0.05 km处的理论波形，如图14所示。以此近台站理论波形作为零相位化处理中所用的震源子波与合成数据（图6a）的波形做反褶积。图15为震中距分别为20、50、100 km零相位化波形结果对比，图中红线波形为以近台记录近似震源子波得到的零相位化结果，从图中可看出与以真实子波的零相位化结果（图中黑线所示）存在的差异较小，但仍有必要对主动源实验中气枪震源机制进行深入研究。
6結论
从合成数据和实际数据的处理结果和分析中可得到以下结论：
（1）零相位处理可以提高波形数据的分辨率，零相位化后有效信号到时点变为波形峰值点，从而减小了叠加过程中噪声的干扰，经非线性叠加后可提高震相到时的准确性。在比较气枪震源的零相位子波后，本文认为以反褶积方法作为零相位化处理手段更适用于气枪震源有效信号的提取。且在实际数据处理中，使用近震中距的ckt台站记录近似气枪震源的真实子波对零相位化处理结果影响很小。
（2）非线性叠加结果卷积震源子波可获得与原始波形记录较为一致的波形记录，较大程度降低了非线性叠加对有效信息的损失。
（3）使用时间域相位加权和相似性加权叠加方法可获得较好结果，而基于时频域的相位加权类叠加方法对有效波形成分影响较大。基于时间域的相位加权和相似性加权叠加方法对于低信噪比记录中小振幅波形的保幅能力较差，可将小能力信号波形压制而损伤部分有效信息；基于时频域的相位加权叠加类方法可较好保留叠加波形中的小振幅信息（保幅能力与线性叠加相比较差），但会造成信号相位偏移而产生虚假波形，进而影响提取信号的准确性。
参考文献：
陈颙，朱日祥.2005.设立“地下明灯研究计划”的建议[J].地球科学进展，20（5）：485-489.
顾庙元，姚佳琪，张伟，等.2016.地学长江计划安徽实验中低频可控震源地震波信号提取方法评估[J].中国地震，32（2）：356-378.
林建民，王宝善，葛洪魁，等.2008.大容量气枪震源特征及地震波传播的震相分析[J].地球物理学报，51（1）：206-212.
倪宇东，祖云飞，李海翔，等.2010.可控震源地震数据初至时间拾取方法[J].石油地球物理勘探，45（6）：793-796.
倪宇东.2014.可控震源地震勘探采集技术[M].北京：石油工业出版社.
陶知非，赵永林，马磊.2011.低频地震勘探与低频可控震源[J].物探装备，21（2）：71-76.
王宝善，葛洪魁，王彬，等.2016.利用人工重复震源进行地下介质结构及其变化研究的探索和进展[J].中国地震，32（2）：168-179.
王彬，杨润海，王宝善，等.2012.地震波走时变化精确测量的实验研究[J].云南大学学报（自然科学版），34（S2）：15-20.
武安绪，叶泵，李红，等.2016.气枪震源弱信号提取的可靠性初步分析[J].中国地震，32（2）：319-330.
肖卓，高原.2015.尾波干涉原理及其应用研究进展综述[J].地震学报，（3）：516-526.
熊翥.2008.地震数据处理应用技术[M].北京：石油工业出版社.
杨微，王宝善，葛洪魁，等.2013.大容量气枪震源主动探测技术系统及试验研究[J].中国地震，29（4）：399-410.
翟秋实，姚华建，王宝善.2016.气枪震源资料反褶积方法及处理流程研究[J].中国地震，32（2）：295-304.
张军华.2011.地震资料去噪方法[M].青岛：石油大学出版社.
BRITTLE K F，LINES L R，DEY A K.2001.Vibroseis deconvolution：a comparison of cross correlation and frequency domain sweep decnvolution[J].Geophysical Prospecting，49（6）：675-686.
CROTWELL H P，OWENS T J，RITSEMA J.1990.The TauP Toolkit：Flexible Seismic Travel-time and Ray-path Utilities[J].Seismological Research Letters，70（2）：154-160.
DZIEWONSKI A M，ANDERSON D L.1981.Preliminary reference Earth model[J].Physics of the Earth & Planetary Interiors，25（4）：297-356.
HELMBER D，WIGGINS R A.1971.Upper mantle structure of Mid western united states[J].Journal of Geophysical Resaearch，76（14）：3229-3245.
KENNETT B L N.1987.Observational and theoretical constraints on crustal and upper mantle heterogeneity[J].Physics of the Earth & Planetary Interiors，47（1）：319-332.
LI Q，GAO J.2014.Application of Seismic Data Stacking in Time–Frequency Domain[J].Geoscience & Remote Sensing Letters IEEE，11（9）：1484-1488.
NEIDELL N S，TANER M T.1971.Semblance and other coherency measures for multichannel data[J].Geophysics，36（3）：482-497.
PINNEGAR C R，EATON D W.2003.Application of the S transform to prestack noise attenuation filtering[J].Journal of Geophysical Research：Solid Earth，108（B9），doi：10.1029/2002JB002258.
SCHIMMEL M，GALLART J.2007.Frequency-dependent phase coherence for noise suppression in seismic array data[J].Journal of Geophysical Research：Solid Earth，112（B4）：1-10.
SCHIMMEL M，PAULSSEN H.1997.Noise reduction and detection of weak，coherent signals through phase-weighted stacks[J].Geophysical Journal of the Royal Astronomical Society，130（2）：497-505.
SHEARER P M.2009.Introduction to seismology[M].Oxford city：Cambridge University Press.
STOCKWELL R G，MANSINHA L，LOWE R P.1996.Localization of the complex spectrum：the S transform[J].IEEE Transactions on Signal Processing，44（4）：998-1001.
STOCKWELL R G.1999.S-transform analysis of gravity wave activity from a small scale network of airglow imagers[J].Thesis，4662.
TANER M T，KOEHLER F.1969.Velocity spectra-digital computer derivation and applications of velocity functions[J].Geophysics，34（6）：859-881.
ZHU L，RIVERA L A.2002.A note on the dynamic and static displacements from a point source in multilayered media[J].Geophysical Journal International，148（3）：619-627.
摘要：使用零相位化处理和非线性叠加方法提取低信噪比记录中的有效信号，通过定量评估结果波形的走时和振幅等信息的精准程度，探究此处理流程是否可提高提取气枪震源有效信号波形的准确性。对合成波形和云南宾川主动源实验中实际资料数据进行零相位化处理，选用相似性加权、时间域相位加权、时频域相位加权、改进的时频域相位加权4种非线性叠加方法提取有效信号。结果表明，基于上述处理方法可有效提高提取波形的准确性：（1）零相位处理可以提高波形数据的分辨率，使有效信号到时点变为波形峰值点，有利于在非线性叠加后拾取出准确的到时位置；（2）零相位化数据的非线性叠加结果卷积震源子波可基本恢复地震信号波形记录，可降低非线性叠加对有效信息的损失；（3）使用时间域相位加权叠加和相似性加权叠加方法可获得较好结果，但可能会压制低信噪比小振幅信号；基于时频域的相位加权类叠加方法对有效波形成分影响较大，但对小振幅信号保幅较好。
关键词：主动源；气枪震源；零相位化处理；非线性叠加；信号提取
中图分类号：P315.3文献标识码：A文章编号：1000-0666（2017）04-0581-14
0引言
人工震源技术可以克服传统天然震源时空分布不均的制约，为获得高时空分辨率的区域尺度地下结构和介质变化信息提供可能（陈颙，朱日祥，2005；王宝善等，2016）。与天然震源相比，人工震源单次激发能量较低、远距离传播易受噪声干扰，需要大量重复激发再从中提取有效信号（张军华，2011；杨微等，2013；倪宇东，2014；顾庙元等，2016；武安绪等，2016）。提取结果的可靠性将直接影响后续地震资料处理结果的准确性。因此，如何从低信噪比的重复激发地震记录中获取高质量有效信号，是人工震源实验中的重点研究问题之一。
非线性叠加类方法是从低信噪重复激发记录中提取有效信号最为常用的方法之一。此类方法利用了有效信号的一致性，从而实现保留有效随机信号、压制噪声的目的（Kennett，1987；Schimmel，Paulssen，1997；Schimmel，Gallart，2007；Li，Gao，2014）。相比线性叠加，非线性叠加可大幅提高信噪比，但会造成信噪比过低的信号发生波形畸变和相位偏移等（武安绪等，2016）。在主动源实验中，非线性叠加方法也被应用于大容量气枪震源、低频可控震源等地震资料有效信号提取中。武安绪等（2016）使用线性叠加、时间域相位加权叠加、时频域相位加权叠加方法对接收到的波形記录直接处理并定量评估提取结果的可靠性。经合成测试和实际数据处理后，认为加权叠加会影响有效信号提取的波形质量和时间精度。与大容量气枪震源不同，低频可控震源选用震动时间较长的扫频信号激发且激发能量更低（陶知非等，2011）。为了降低不同震相的地震信号重叠干扰并提高数据分辨率，需要对接收到的原始波形做零相位化处理，通常是将地震记录与震源子波做互相关，将地震记录转化为地下介质脉冲响应与震源信号自相关子波的卷积记录（熊翥，2008；倪宇东，2014）。顾庙元等（2016）使用线性叠加、非线性叠加、直接分离类多种方法对主动源实验中零相位化的低频可控震源数据进行处理，获得了较好的有效信号信息。本文将零相位化处理流程引入气枪震源资料有效信号提取中，通过定量评估提取结果波形的走时和振幅等信息的准确性，探究使用零相位化数据是否可降低非线性叠加对低信噪比记录中有效信号波形的影响。
本文首先使用相似性加权叠加、时间域相位加权叠加、时频域相位加权叠加、改进的时频域相位加权叠加4种非线性叠加方法提取未经零相位化处理的加噪合成低信噪比信号，提出了非线性叠加直接提取气枪震源信号所存在的问题。然后，引入零相位化处理流程，分别使用互相关和反褶积两种常用的零相位化处理方法对气枪震源子波零相位化，对比两种处理手段得到的零相位化子波形态特征来确定适合气枪震源资料的零相位化处理方法。最后，通过对合成地震剖面数据和云南宾川主动源实验实际数据的处理，探究并验证先对数据做零相位化处理再叠加的可行性及优势，为主动源实验的信号提取提供参考。
1非线性叠加
1.1方法原理
非线性叠加是在线性叠加基础上发展起来的，是基于记录中有效信息波形或相位等的一致性来计算各采样点的叠加权重以提高信噪比，又称加权叠加。此类叠加方法既可在时间域中实现，又可在时频域中实现。时间域非线性叠加可表示为：yws（t）=w（t）γ yls（t）（1）式中：yws（t）为时间域加权叠加结果；w（t）表示时间域的加权系数；γ为加权指数因子，用于控制加权强度，一般在1.5～2.5之间；yls（t）为线性叠加结果，即：yls（t）=1N∑Nj=1xj（t）（2）式中：N为叠加次数；xj（t）为第j次激发所接收的地震记录；yls（t）为线性叠加结果。而时频域非线性叠加利用了有效信号时频域分布的一致性，可表示为：yws（t）=T-1wγ（t，f）1N∑Nj=1T{xj }（t，f）（3）式中：T和T-1分别表示某种时频变换及其逆变换算子，通常使用Stockwell变换及其逆变换；w（t，f）为时频域下的加权系数。下面，介绍本文所用非线性叠加方法的具体原理。
1.1.1相似性加权叠加
相似性加权叠加是以波形相似性系数（Semblance）作为一致性测量准则而实现非线性叠加的（Kennett，1987）。所谓波形相似性系数是多道记录叠加后瞬时总能量与各单道记录瞬时能量总和之比，用于测量多道记录一致性程度的常用参数（Taner，Koehler，1969；Neidell，Taner，1971），可直接用于加权系数计算：w（t）=∑n/2i=-n/2wG（i）∑Nj=1sj（t+i）2N∑n/2i=-n/2wG（i）∑Nj=1s2j（t+i）（4）式中：wG（i）为高斯窗函数exp[-（τ-t）22k2]，这里k为高斯窗宽度因子，一般应大于震源子波宽度，合理选取可使加权系数光滑稳定；n为窗函数的个数。
1.1.2时间域相位加权叠加
时间域相位加权叠加是以波形记录的瞬时相位作为衡量有效信号一致性依据的叠加方法（Schimmel，Paulssen，1997），其加权系数计算公式如下：w（t）=1N∑Nj=1exp[iΦj（t）]（5）式中：Φj（t）第j次激发记录随时间变化的瞬时相位，可通过Hilbert变换得到。由接收到的地震记录x（t）和其Hilbert变换形式H[x（t）]可构建出解析信号：X（t）=x（t）+iH[x（t）]（6）式中：X（t）表示解析信号。再将解析信号用指数形式表达：X（t）=A（t）exp[iΦ（t）]（7）式中：A（t）为局部振幅；Φ（t）为瞬时相位，是用于判断有效信号一致性、计算加权系数的重要参数。
1.1.3时频域相位加权叠加
时频域相位加权叠加是时间域相位加权在时频域中的拓展，这种方法是以地震记录的时频域瞬时相位的一致性作为评估依据（Schimmel，Gallart，2007）。该方法是通过Stockwell变换将信号变换到时频域中（Stockwell et al，1996，1999），其加权系数的表达式为：w（τ，f）=1N∑Nj=1Sj（τ，f）exp（i2πfτ）Sj（τ，f）（8）式中：Sj（τ，f）表示信号经Stockwell变换的时频谱。Stockwell变换的表达式可写为：Sj（τ，f）=∫∞∫-∞xj（t） wG（t-τ，f）exp（-i2πfτ）dt（9）这里，高斯窗函数wG（t-τ，f）随频率变化，其表达式为：wG（t，f）=f·（2π）-1/2·exp（-t2f2/2）（10）在時频域中实现相位叠加后，再通过Stcokwell逆变换后可得到叠加后的波形记录，如（3）式所示。
1.1.4改进的时频域相位加权叠加
改进的时频域相位加权叠加起初是为解决CMP道集叠加中的时移问题（Li，Gao，2014）而提出的，这里将其用于提取低信噪比信号。此方法是在基于Stockwell变换的自适应滤波叠加基础上发展而来的（Pinnegar，Eaton，2003），它不仅以时频域的瞬时相位为衡量一致性依据，同时也考虑了时频域中信号局部振幅的一致性程度。与以上加权叠加方法不同，改进的时频域相位加权叠加需要计算每一道记录的加权系数，以第k次重复激发记录的权重系数为例：wk（τ，f）=Uk（τ，f）-Vk（τ，f）maxτ，f{Uk（τ，f）-Vk（τ，f）}（11）其中，Uk（τ，f）=1N-1∑Nj=1，j≠kSj（τ，f）maxτ，f{Sj（τ，f） }v（12）Vk（τ，f）=minτ，f{Uk（τ，f）}（13）第k次记录的加权系数是通过其他各记录叠加、归一化等计算得到的。然后对每条记录乘以相应的权重系数后叠加，最终叠加结果为：y（t）=S-11N∑Nk=1Sk（τ，f）wk（τ，f）（14）该方法在实现过程中进行了两次叠加，因此很大程度上压制了随机噪声，保留并突出有效信号，但计算成本较高。
1.2非线性叠加提取有效信号存在的问题
这里使用简单的合成波形来对比评估上述4种非线性叠加方法的信号提取效果。首先使用云南宾川实验CKT台记录到的气枪震源激发信号作为震源子波（图1）构建时长10 s、采样间隔0.01 s、有效信号到时3 s、最大绝对振幅为1的合成波形（图2a）。再向合成波形中加入标准差为1、均值为0的高斯随机噪声（图2b）得到了加噪波形记录。重复上述操作，构建出300道加噪记录。
对上述加噪波形带通滤波2～8 Hz后（图2c）做叠加，结果如图3a所示。由图可知，非线性叠加相比于线性叠加可有效压制噪声、提高信噪比，但会对波形等有效信息造成影响：（1）相似性加权叠加和时间域相位加权叠加方法在压制噪声的同时，也会压制有效信号中小振幅波形使得小振幅波形信息难以完好保留；（2）基于Stockwell变换的时频域相位加权叠加和改进的时频域相位加权叠加方法虽然可以提取出小振幅波形信息，但会造成相位偏移产生虚假波形；（3）由于震相到时为波形起跳点，基于此从结果中拾取到时，从拾取结果（图3中深灰色竖线所标）可看出震相的到时难以准确拾取；（4）线性叠加虽不会对信号波形造成影响，但随机噪声难以压制完全，且当叠加次数有限时仍难以从结果中拾取出准确的到时位置。
由以上结果可看出，直接对气枪震源记录非线性叠加会对波形造成影响，甚至可能导致初至起跳时间拾取错误。之所以初至时间会被非线性叠加破坏，是因为对于因果信号（最小相位信号），初至时间位于能量较弱的起跳点。如果地震震源采用零相位子波激发，则初至时间位于零相位子波的峰值，有可能在非线性叠加过程中保留下来。本文尝试先对气枪信号进行零相位化处理，再进行非线性叠加，考察非线性叠加对零相位化后的信号走时和振幅等波形信息的定量影响，进而对能否应用非线性叠加提取气枪震源激发的有效信号提供理论依据。
2零相位化处理
地震波的传播可以看作是一个线性系统，符合褶积模型：x（t）=r（t）*s（t）（15）式中：x（t）为台站接受到的地震记录；r（t）为地下介质的脉冲响应（天然地震学中称为格林函数）；s（t）为震源子波信号。由此可见，地震记录可近似认为是地层脉冲响应与震源子波卷积的结果，各震相实际是混合相位的（熊翥，2008）。对于接收到的记录，震相的到时位置为波形起跳点，但往往因噪声的影响而难以准确拾取走时（倪宇东等，2010）。零相位化处理可将混合相位波形转化为零相位子波，从而使震相到时点变为零相位子波的峰值点，易拾取有效信号信息。此外，零相位化处理可将震源子波压缩为近似δ函数的零相位子波，提高数据分辨率，避免了多震相相互叠加、干扰等。下面介绍互相关、反褶积2种常用的零相位化处理方法原理。
2.1方法原理
2.1.1互相关
互相关方法是把地震记录与震源子波作互相关，将记录中的震源子波信号转化为零相位的互相关子波。互相关方法常用于可控震源资料处理中（图4a-1、a-2）。地震记录与震源子波做互相关后可表示为：x（t）s（t）=r（t）*s（t）s（t）（16）式中：表示互相关算子。对上式简化可表示为：x（t）s（t）=r（t）*k（t）（17）式中：k（t）表示震源信号互相关子波，这是零相位化的子波，此时震相的起跳点变为峰值点。
2.1.2反褶积
反褶积方法是将地震记录中的震源子波去除掉，从而得到具有较高分辨率的地下介质脉冲响应。反褶积处理被广泛应用于地震学偏移处理等中（熊翥，2008；Shearer，2009），气枪震源资料处理中也常用此种手段提取传播路径的格林函数。反褶积的方法有很多，最直接的处理方法是频率域反褶积，即将式（15）作Fourier变换将其转换至频率域中：X（ω）=R（ω）S（ω）（18）式中：X（ω）为地震记录x（t）的频谱；R（ω）为地下介质脉冲响应r（t）的频谱；S（ω）为震源子波的频谱。将X（ω）与S（ω）相除可得到较高分辨率的地下介质信息：R（ω）=X（ω）S（ω）（19）因实际数据中混合有随机噪声，频率域直接相除常由于分母频谱出现极小值或零点，造成结果不稳定。为了避免此种现象，这里选用计算稳定且效率高的频率域水准反褶积方法（Helmberger et al，1971；翟秋实等，2016）进行后续基于反褶积方法的零相位化处理。
2.2气枪震源的零相位化处理方法
根据以上常用的零相位化处理方法，本文将探究适合气枪震源的零相位化方式。首先分别使用互相关和反褶积方法获得气枪震源的零相位子波波形，对比2个子波形态以选择适合气枪震源资料的零相位化处理方法。为了更好地对比评估，选择适合气枪震源的零相位化方法，这里同时对云南宾川实验气枪震源ckt台波形（图4b-1）和长江计划安徽实验中低频可控震源理论子波（图4a-1）做零相位化处理，然后對比相应零相位化子波波形特征。在使用主动源进行探测试验中，可控震源资料必须经零相位化处理后获得较高分辨率数据再进行后续处理（Brittle et al，2000；倪宇东，2014；顾庙元等，2016），而气枪震源较少使用此流程来提高数据质量。此外，可控震源由于扫频信号自身的特性，其互相关子波主瓣明显、旁瓣噪声低可近似为δ函数，且互相关方法实现简单且稳定（特别是对于低信噪比记录），所以对于此类数据最常用的零相位化方法是互相关方法。
经互相关后得到的震源自相关子波波形分别如图4a-2、b-2所示。使用反褶积方法进行零相位化处理后的子波波形如图4a-3、b-3所示。
从图4结果中可看出，气枪震源的自相关子波（图4b-2）不具备可控震源自相关子波（图4a-2）高清晰度的特征，且时长过长而难以达到压缩震源子波提高分辨率的目的，而利用反褶积得到的子波波形（图4b-3）具有很好的分辨率。基于此结果，本文认为对于气枪震源数据的零相位处理使用反褶积方法可获得更好的零相位化处理效果，后续使用频率域水准反褶积来实现气枪震源数据的零相位化。
为了探究零相位化方法是否可以降低非线性叠加对波形信息的影响，这里对1.2中的300道加噪合成数据进行基于反褶积方法的零相位化处理，然后再进行非线性叠加。单道加噪合成记录（图2c）的零相位化处理结果如图2d所示（经2～8 Hz滤波），零相位化后300次叠加结果如图3b所示。对图3b各结果设置2～4 s波形对比时窗定量评估提取波形的准确度，如图中文字所标。由结果可看出：（1）非线性叠加可以很好地压制随机噪声，并且较直接叠加结果（图3a）相比可更为准确地保留波形信息；（2）经零相位化处理后震相拾取位置由拾取起跳点变为拾取波形峰值点，经非线性叠加后峰值点位置不受影响，因此震相到时仍可准确拾取；（3）对比4种非线性叠加结果，可以看出相似性加权叠加和时间域相位加权叠加可较好保留零相位子波信息并压制子波旁瓣，而基于Stockwell变换的时频域相位加权和改进的时频域相位加权方法会导致子波波形相位变化而产生虚假波形，但波形峰值位置即有效信号的到时不受影响。由此可见，数据经零相位化后可降低非线性叠加对结果的影响，有利于保证后续处理的准确性。
3合成地震剖面数据处理测试
为了更准确地评估零相位化处理是否可降低非线性叠加对有效信号的影响，本文使用一维速度模型合成地震剖面数据来模拟大容量气枪震源信号处理过程。
3.1数据合成
首先在PREM全球平均速度模型（Dziewonski，Anderson，1981）的基础上（图5a），仿照真实地下介质情况构建一维速度结构用于合成仿真剖面数据，即在0～0.5 km近地表深度范围内加入了强衰减的低速层，衰减系数设置为QP=80、QS=40，并对加入了0.5～24.4 km地壳中设置速度逐渐变化的薄层，构建出的一维速度模型分布，如图5b所示。然后使用FK程序（Zhu，Rivera，2002）计算上述一维速度结构（图5b）的理论波形。气枪震源的震源机制复杂，这里将大容量气枪简化为作用于地表的垂直单力源，震源子波采用云南宾川主动源实验场ckt台接收单次激发记录（图1）。为了便于研究体波信息，这里从正演所得波形中去除面波成分，最终得到了该速度结构的理论地震剖面记录（图6a）。并使用Taup震相走时计算程序（Crotwell et al，1990）计算对应的Pg、Sg、Pn、Sn震相走时，如图6中所示。
然后，模仿真实噪声情况向理论地震剖面中加入随机噪声，构建合成地震剖面数据。本文通过测量震中距为5 km左右的台站接收的实际记录中信噪比来确定所加随机噪声强度。这里所定义的信噪比为有效信号时段内有效信号的最大振幅与背景噪声均方根振幅之比（林建民等，2008）。图7为云南宾川主动源实验中震中距为5.009 6 km的53260台大容量气枪震源单次激发记录，从记录中设置Pg波到时开始至10 s区间范围内的信号窗和时窗长度为5 s的噪声窗，分别如图中深灰色阴影和浅灰色阴影区域所示，信噪比为244.437 4。根据测量出的信噪比，向震中距为5 km处的理论波形中加入相应强度的高斯随机噪声。由于接收台站的噪声水平基本一致，因此对其他震中距的理论合成记录中加入与5 km台站相同标准差的高斯随机噪声，合成出的单次激发仿真地震剖面如图6c所示。重复上述加噪操作，构造出300次重复激发仿真数据。
3.2数据处理与分析
首先，对加噪合成数据进行去均、去势、去尖灭、2～8 Hz带通滤波等预处理。然后，对预处理后的数据直接叠加得到线性及非线性叠加结果，分别如图8a、c所示；而对于零相位化数据的叠加，则先使用频率域水准反褶积对预处理后数据作零相位化处理，然后进行2～8 Hz带通滤波处理再叠加，最后得到线性及非线性叠加结果分别如图8b、d所示。
这里抽出震中距分别为20、50、100 km合成数据的叠加结果（图9）以更好地观察提取有效信号效果。根据Taup震相走时计算程序计算出的震相到时，可以从结果中看出直接叠加和零相位化后叠加的波形特征。对于直接叠加走时点位置为波形起跳点，但非线性叠加后初至时间破坏；而对于零相位化后叠加的结果，走时点位置为该震相时段内波形出现的第一个振幅明显的峰值点处。因此本文认为：（1）低振幅的震相到时起跳点在叠加过程中振幅发生变化，导致很难从结果中拾取出准确到时；（2）而零相位化后数据的叠加结果可以避免上述问题，较好地保留有效信号信息；（3）将未零相位化数据的叠加结果和零相位化数据的叠加结果相结合可用于震相走时拾取。
4实际数据处理分析
将前文所用的先零相位化后非线性叠加处理流程应用于2013年云南宾川主动源实验大容量气枪震源地震资料实际数据处理中，验证零相位化处理是否可以提高提取信号的可靠程度。本文选取了该实验中的53280、53259、53270台（震中距分别为18.108 3 km、46.545 2 km、91.609 4 km）记录。对此实际数据进行了预处理、零相位化、叠加等流程处理，其中预处理流程中除去均、去势、去尖灭及带通滤波处理外，还需进行去仪器响应和分量旋转。根据合成测试中的零相位化子波波形和震相到时点位置波形等特征，本文将气枪震源波形记录直接线性叠加结果和零相位化后的波形记录非线性叠加结果相结合拾取震相到时，本文认为的到P、S波时位置分别如图10中灰色实、虚线所示。
5讨论
为定量评估各方法提取有效信号的效果，本文选用最大互相关系数R和时间延迟Td来评估提取波形的准确程度。最大互相关系数是所选取的波形对比时窗内的叠加结果与对应的理论波形的互相关系数最大值，时间延迟是互相关系数最大值所对应的时间偏移大小（王彬等，2012；肖卓，高原，2015；武安绪等，2016）。然后分别从图3a、b两组处理结果选取2～6 s和2～4 s的波形对比时窗，计算出的R、Td大小如图中文字所标。
5.1非线性叠加对各频段成分的影响
选取1.2和1.3中的非线性叠加结果，以时间域相位加权叠加结果为例分析非线性叠加结果的各频率成分的影响。从频率范围为2～8 Hz的直接叠加结果（图3a）和零相位化后叠加结果波形（图3b）中分离出2～4、4～6、6～8 Hz频段波形，如图11红线所示。将分离出的各频段成分与对应频段的理论信号做比较（图11文字标出最大互相关系数R的大小）结果显示，非线性叠加对4～6 Hz频率成分影响最小，零相位化数据叠加结果各频率成分保持良好。而大容量气枪震源主频在5 Hz左右，可以看出主频范围的信号成分保留较好。
5.2非线性叠加对振幅的影响
为了探究各叠加方法是否可以保留有效信号振幅，这里使用1.2中的合成波形（图2a）构建一组振幅变化的合成数据，最大振幅从1至0.1，间隔0.1变化，共10道数据，如图12a所示。再向数据中加入标准差为1的高斯随机噪声，产生300组数据，然后进行叠加。
从叠加结果的最大振幅变化情况（图12b、c）中可看出：（1）非线性叠加可以基本保持未零相位化和零相位化数据的振幅信息，相似性加权和时间域相位加权叠加振幅结果与线性叠加基本一致，基本吻合真实振幅；（2）而时频域相位加权叠加类方法虽基本拟合真实振幅，但是所得振幅偏小；（3）此外对于最大振幅为0.1的记录，零相位化处理后数据的相似性加权叠加和时间域相位加权叠加结果的振幅保持效果较基于时频域的相位加权叠加结果差，此测试结果与图9b中远震中距叠加结果相似。本文认为當信噪比过低时，相似性加权叠加和时间域相位加权叠加仅利用了时间域上的波形和相位信息，在压制噪声提高信噪比的同时，也将小振幅的有效信号信息压制而无法获得准确的波形信息；而时频域相位加权叠加充分利用了有效信号的时频分布特征，可从低信噪比记录中提取出有效信息。
5.3零相位化非线性叠加结果恢复波形能力的讨论由图3所标文字可知，对于气枪震源合成波形直接叠加，当叠加次数足够时，非线性叠加的提取效果较线性叠加差；而对于零相位化的气枪震源合成波形叠加，相似性加权叠加和时间域相位加权叠加的提取效果优于线性叠加方法；时频域相位加权和改进的时频域相位加权叠加会导致有效信号相位偏移而产生虚假波形，所提取的结果与真实波形差异较大。这里将零相位化数据叠加结果（图3b）卷积气枪震源子波波形（图2a）测试零相位化的叠加结果是否可用于恢复原始信号并提高提取结果的可靠性，结果如图13所示。设置2～6 s波形对比时窗与气枪震源子波波形（图3a）进行比较发现（图13文字所标），零相位化数据非线性叠加结果的最大互相关系数高于数据直接进行非线性叠加结果的，其中零相位化波形叠加结果中相似性加权和时间域相位加权叠加所恢复的气枪震源波形结果略优于波形直接线性叠加结果。这是由于波形经零相位化处理后原记录中波形的起跳点转化为了零相位子波的峰值点，峰值点的振幅较大，在叠加过程中不易受噪声干扰，因此可以保留较好的走时信息和地下介质的脉冲响应信息；此外，相似性加权叠加和时间域相位加权叠加可以一定程度上压制零相位子波的旁瓣噪声，从而进一步提高数据的分辨率，而基于时频域的相位加权类叠加方法由于Stockwell变换分辨率等问题造成相位偏移等现象，使得最终经逆变换至时间域后产生虚假波形而影响数据质量。
5.4以近台站记录近似震源子波对零相位化结果的影响由于气枪震源的震源机制复杂，尚未有很好的理论震源子波，目前是在距震源0.05 km处放置参考台站（CKT台），使用参考台接收记录作为震源子波做后续处理。下面探究采用近台站记录作为理论震源子波是否会影响零相位化结果的准确性。使用FK程序计算3.1合成测试中所用一维速度模型下震中距为0.05 km处的理论波形，如图14所示。以此近台站理论波形作为零相位化处理中所用的震源子波与合成数据（图6a）的波形做反褶积。图15为震中距分别为20、50、100 km零相位化波形结果对比，图中红线波形为以近台记录近似震源子波得到的零相位化结果，从图中可看出与以真实子波的零相位化结果（图中黑线所示）存在的差异较小，但仍有必要对主动源实验中气枪震源机制进行深入研究。
6結论
从合成数据和实际数据的处理结果和分析中可得到以下结论：
（1）零相位处理可以提高波形数据的分辨率，零相位化后有效信号到时点变为波形峰值点，从而减小了叠加过程中噪声的干扰，经非线性叠加后可提高震相到时的准确性。在比较气枪震源的零相位子波后，本文认为以反褶积方法作为零相位化处理手段更适用于气枪震源有效信号的提取。且在实际数据处理中，使用近震中距的ckt台站记录近似气枪震源的真实子波对零相位化处理结果影响很小。
（2）非线性叠加结果卷积震源子波可获得与原始波形记录较为一致的波形记录，较大程度降低了非线性叠加对有效信息的损失。
（3）使用时间域相位加权和相似性加权叠加方法可获得较好结果，而基于时频域的相位加权类叠加方法对有效波形成分影响较大。基于时间域的相位加权和相似性加权叠加方法对于低信噪比记录中小振幅波形的保幅能力较差，可将小能力信号波形压制而损伤部分有效信息；基于时频域的相位加权叠加类方法可较好保留叠加波形中的小振幅信息（保幅能力与线性叠加相比较差），但会造成信号相位偏移而产生虚假波形，进而影响提取信号的准确性。
参考文献：
陈颙，朱日祥.2005.设立“地下明灯研究计划”的建议[J].地球科学进展，20（5）：485-489.
顾庙元，姚佳琪，张伟，等.2016.地学长江计划安徽实验中低频可控震源地震波信号提取方法评估[J].中国地震，32（2）：356-378.
林建民，王宝善，葛洪魁，等.2008.大容量气枪震源特征及地震波传播的震相分析[J].地球物理学报，51（1）：206-212.
倪宇东，祖云飞，李海翔，等.2010.可控震源地震数据初至时间拾取方法[J].石油地球物理勘探，45（6）：793-796.
倪宇东.2014.可控震源地震勘探采集技术[M].北京：石油工业出版社.
陶知非，赵永林，马磊.2011.低频地震勘探与低频可控震源[J].物探装备，21（2）：71-76.
王宝善，葛洪魁，王彬，等.2016.利用人工重复震源进行地下介质结构及其变化研究的探索和进展[J].中国地震，32（2）：168-179.
王彬，杨润海，王宝善，等.2012.地震波走时变化精确测量的实验研究[J].云南大学学报（自然科学版），34（S2）：15-20.
武安绪，叶泵，李红，等.2016.气枪震源弱信号提取的可靠性初步分析[J].中国地震，32（2）：319-330.
肖卓，高原.2015.尾波干涉原理及其应用研究进展综述[J].地震学报，（3）：516-526.
熊翥.2008.地震数据处理应用技术[M].北京：石油工业出版社.
杨微，王宝善，葛洪魁，等.2013.大容量气枪震源主动探测技术系统及试验研究[J].中国地震，29（4）：399-410.
翟秋实，姚华建，王宝善.2016.气枪震源资料反褶积方法及处理流程研究[J].中国地震，32（2）：295-304.
张军华.2011.地震资料去噪方法[M].青岛：石油大学出版社.
BRITTLE K F，LINES L R，DEY A K.2001.Vibroseis deconvolution：a comparison of cross correlation and frequency domain sweep decnvolution[J].Geophysical Prospecting，49（6）：675-686.
CROTWELL H P，OWENS T J，RITSEMA J.1990.The TauP Toolkit：Flexible Seismic Travel-time and Ray-path Utilities[J].Seismological Research Letters，70（2）：154-160.
DZIEWONSKI A M，ANDERSON D L.1981.Preliminary reference Earth model[J].Physics of the Earth & Planetary Interiors，25（4）：297-356.
HELMBER D，WIGGINS R A.1971.Upper mantle structure of Mid western united states[J].Journal of Geophysical Resaearch，76（14）：3229-3245.
KENNETT B L N.1987.Observational and theoretical constraints on crustal and upper mantle heterogeneity[J].Physics of the Earth & Planetary Interiors，47（1）：319-332.
LI Q，GAO J.2014.Application of Seismic Data Stacking in Time–Frequency Domain[J].Geoscience & Remote Sensing Letters IEEE，11（9）：1484-1488.
NEIDELL N S，TANER M T.1971.Semblance and other coherency measures for multichannel data[J].Geophysics，36（3）：482-497.
PINNEGAR C R，EATON D W.2003.Application of the S transform to prestack noise attenuation filtering[J].Journal of Geophysical Research：Solid Earth，108（B9），doi：10.1029/2002JB002258.
SCHIMMEL M，GALLART J.2007.Frequency-dependent phase coherence for noise suppression in seismic array data[J].Journal of Geophysical Research：Solid Earth，112（B4）：1-10.
SCHIMMEL M，PAULSSEN H.1997.Noise reduction and detection of weak，coherent signals through phase-weighted stacks[J].Geophysical Journal of the Royal Astronomical Society，130（2）：497-505.
SHEARER P M.2009.Introduction to seismology[M].Oxford city：Cambridge University Press.
STOCKWELL R G，MANSINHA L，LOWE R P.1996.Localization of the complex spectrum：the S transform[J].IEEE Transactions on Signal Processing，44（4）：998-1001.
STOCKWELL R G.1999.S-transform analysis of gravity wave activity from a small scale network of airglow imagers[J].Thesis，4662.
TANER M T，KOEHLER F.1969.Velocity spectra-digital computer derivation and applications of velocity functions[J].Geophysics，34（6）：859-881.
ZHU L，RIVERA L A.2002.A note on the dynamic and static displacements from a point source in multilayered media[J].Geophysical Journal International，148（3）：619-627.
