宾川主动源地震波走时变化分析方法研究*

刘自凤，苏有锦，王宝善，王　彬，杨　军，李孝宾

(1.云南省地震局，云南 昆明 650224；2.中国地震局地球物理研究所，北京 100081)

宾川主动源地震波走时变化分析方法研究*

刘自凤1，苏有锦1，王宝善2，王彬1，杨军1，李孝宾1

(1.云南省地震局，云南 昆明 650224；2.中国地震局地球物理研究所，北京 100081)

摘要：挑选了2013年以来信噪比较高、数据较为完整的6个宾川主动源接收台站的气枪记录数据，采用互相关时延检测、近场校正等技术进行了高精度波速变化测量，得到以下初步结果：(1)接收台站下方介质存在δv/v≈10-5～10-2的相对波速变化；(2)2014年9月初至10月底，初至波走时变化较平稳，集中在-0.01～0.01 s之间变化，且6个台的变化趋势一致。

关键词：主动源；互相关时延检测；气枪信号；走时变化

0引言

地震波在地球介质中传播，带来了丰富的地下介质物性信息，为我们了解地球内部结构及运动变化提供了可能。地震波的很多性质(如衰减、各向异性等)也都可以用来研究地下介质应力状态的变化，但地震波速变化仍然是测量精度最高的研究方法。因此，地球内部地震波速的差异是确定地下地层结构和横向不均匀性的重要物理参数，地震波速的测量一直是地震学中一个令人注目的领域(徐荟，2013)。

波速测量是一个很经典的物理问题，如光速的测量、声速的测量等。在实验室内，测量岩石弹性波速度的方法基本上有3种：共振法、脉冲法(pulse transmission)和超声干涉法(ultrasonic interferometry)(王彬等，2012)。共振法特别适用于测量高温下介质的横波速度；脉冲法至今仍然在实验室内被广泛采用，野外小尺度在波速测量中也有所应用；干涉法是目前几种常用波速测量方法中精度最高的一种方法，可进一步细分为相位比较法和脉冲叠加法。由于观测技术的发展，人们越来越重视对地震波速度变化的精确测量与研究。由间接测量波速变化的原理可知，相对波速变化可以通过相对走时变化间接测量(徐荟，2013)。近年来，互相关检测技术在信号检测领域获得越来越多的应用，主要包括从噪声中提取信号时延估计、速度检测、距离检测和系统动态特性识别等。基于这一原理的地震波速度测量和尾波干涉(CWI)方法近年来已取得了一些突出的成果(Snieder，2006；Stehlyetal.，2007；罗桂纯等，2008)，在地震物理预报、监测地下应力变化、研究应力和波速变化的对应关系等方面都有重要的应用(罗桂纯等，2008)。

2011年中国地震局地球物理研究所和云南省地震局在多震的滇西北宾川大银甸水库顺利建成了全球首个常规观测的大容量气枪人工震源发射台——宾川主动震源发射台。发射台地处红河断裂带和程海断裂带交汇区中部。作为川滇菱形块体的西侧边界，红河断裂以西是著名的三江地槽褶皱系，以东则被金河—程海—宾川断裂一分为二。大致以弥渡一带为界北部称甘孜—丽江断块，南部称石棉—楚雄断块。红河断裂北部、南部在大地构造与深部特征、地热与地球化学场、断裂结构与力学性质、断层构造岩发育特征以及断裂活动等方面都存在着一系列显著异常。这些不同时间尺度的差异特征最终都集中体现在沿红河断裂带现代地震活动显著的南北差异上(张建国等，1993)。张建国等(1993)的统计结果表明，该断裂北部的地震活动明显强于南部。程海断裂是滇西北地区东缘一条十分醒目的多期活动性大断裂，主断裂及其分支断裂控制着永胜、程海、期纳、宾川、弥渡等盆地的发育。沿着断裂带除了发育诸多晚新生代构造盆地外，还出现了山脊、水系、冲沟的左行走滑位错，局部地区还发育了断层崖，表明程海断裂具有左行走滑兼正断层的性质。程海断裂早期阶段走向总体上近 NS向，倾向西，向南在弥渡附近与红河断裂相交，晚期阶段走向 NNE向(唐渊，刘俊来，2010)。在滇西地震预报实验场核心区域、活动强烈的红河断裂北段与程海断裂之间建立宾川主动源发射台，旨在利用人工震源主动向地下发射地震波，进行地下介质监测，变“被动监测”为“主动探察”，通过连续激发实现动态跟踪地震波走时变化，并进一步分析地震波速的时间演化特征，最终为该地区地震活动后续发展趋势判定提供一定判据。

宾川主动震源发射台于2012年11月投入正常运行，并在其周围布设了由40套短周期地震仪组成的密集观测台阵，形成了一套完整的人工主动震源发射—观测系统，并执行常规观测至今。该激发系统由4条容量为2 000 in3的单枪组合而成，在平均水深为20 m的大银甸水库内保持枪与水面的距离为12 m进行激发，用枪阵高压气体的容量和压力计算(Ronen，2002)得到激发一次能产生8.91×106J的能量，相当于一次ML0.7天然地震。为了跟踪分析宾川气枪源地震波走时差动态变化特征，本文对2013年以来宾川主动源实验数据进行处理，主要介绍互相关检测、近场校正等技术在测量宾川主动源地震波速度变化中的运用及初步处理结果。

1互相关时延检测技术原理

1.1相关函数

如果x(t)、y(t)是能量有限信号，它们的相关函数可定义为

(1)

(2)

式中，*表示复共轭。相关函数是两信号之间时延为τ的函数。若x(t) 与y(t) 不是同一信号， 则它们的相关函数Rxy(τ) 或Ryx(τ) 称为互相关函数。 如果x(t) 与y(t) 是同一信号， 即y(t)=x(t)， 此时， 它们的相关函数Rxx(τ) (或简写作R(τ) ) 称为自相关函数。即

(3)

在实际应用中，信号x(t)、 y(t)一般是实函数而不是虚函数，上述各式仍然适用。在[-∞，+ ∞]的区间里，功率函数是不可积的。通常把这类信号的相关函数定义为

(4)

(5)

(6)

此时有Rxy(τ)=Ryx(-τ)；Rxx(τ)=Rxx(-τ)，即实函数的自相关函数是时延τ的偶函数，τ是所研究两点间的时间间隔，即两信号的时延(Rxx(τ)是自相关函数， Rxy(τ)是互相关函数)。

1.2相似信号的时延检测

利用相关性计算两个记录信号的互相关函数，便可以较准确地求得两个记录信号的时间延迟，从而求得速度分布。其原理是拾震器的信号是由震源信号激振所引起的，源信号在传播过程中由于路径等因素的影响，会叠加许多干扰噪声，但由于这种噪声通常不与源信号同频，即源信号与这些干扰波的信号互相关函数很小，因而将源信号与拾震器信号进行互相关计算时，与源同频的激振信号便可以识别出来，从而得到其从源到拾震器之间的传播时间(Takeshietal.，2005；Clifford，1987；Azizul，1981)。即地震波形记录往往是功率有限信号，如果两个实数波形x(t)和y(t)波形高度相似， 只是时间上存在延迟， 即y(t)是x(t)时延τ得到的， 则我们可以通过计算两个波形的相关系数函数来得到时延τ。

x(t)和y(t)的时延相关系数定义为

(7)

当相关系数取得最大值时，对应的τ就是两个信号的时延。通过计算地震波记录中不同震相的走时变化，就可以确定不同类型波的波速变化。一般都是对于P波、S波以及面波震相来进行计算(徐荟，2013)。

2资料及数据处理过程介绍

2.1资料介绍

从目前收集到的激发记录可知(图1)，2013年以来分3个集中激发时段，分别为2013-01-07～2013-05-08、2013-09-09～2014-05-08、2014-09-08～2014-10-24。根据气枪源附近的参考台记录可知，2013-01-07～2014-10-24一共开展了100组激发实验，共记录到1 798条气枪信号。

气枪源附近布设有40个接收台站，受客观因素影响，并不是每个台站都能记录到每一次激发信号，因此本研究选用气枪源10 km范围内信噪比较高、数据记录缺失较少的6个台站开展地震波走时变化的跟踪分析，气枪源、所选台站的空间分布情况如图2a所示，图2b为2013年1月7日22时40分6个台站接收到的气枪信号，由图可见，6个台记录信号清晰、信噪比高。

图1　2013-01-07～2014-10-24宾川主动源

图2　台站分布示意图(a)及6个台的原始波形(b)

2.2数据处理过程

原始记录波形包含很多无效信息，不能直接使用。因此，在做相关干涉前，需要对原始波形数据进行去均值、去趋势、频谱分析及滤波、数据筛选、反卷积等预处理，处理流程如图3所示。

图3　数据处理流程图

2.2.1去均值和去趋势

为了避免一些非介质变化信号(如信号直流量、趋势项)的影响，在对信号进行处理之前，通常需要对信号进行去均值和去趋势。因为机械振动信号通过传感器拾取、信号调理和 A/D 采样转换成数字信号，最终都含有一定的直流量。如果积分，会在频谱处出现一个冲激，并影响其左右的频谱形状，产生较大误差。

另外由于传感器的零点漂移和现场环境的影响等，振动信号包含一定的线性项或非线性项的慢变成分，这些信号成分积分后被放大为缓慢变化的趋势项。且趋势项叠加在我们所要分析的信号上，对低频谱产生影响(徐荟，2013)。因此，我们分别使用 SAC 命令 Rmean 和 Rtrend 进行去均值和去趋势处理。图4a为2013年1月7日21时59分53264台站接收到的气枪原始记录，图4b为该条记录经过去均值和去趋势处理后的波形，经过处理有效消除了零漂和趋势项影响。

图4　波形预处理过程

2.2.2频谱分析及滤波

地震信号从震源激发并传播到最后被检波器接收，传播过程中会受到噪声的干扰，在低信噪比(SNR)情况下，可能会被噪声掩盖。为了尽可能消除噪声给观测带来的影响，需要对采集到的信号进行滤波处理以提高信噪比。滤波之前，首先要对信号进行频谱分析，判断波的优势能量集中在哪个频段，以此来选择滤波范围。林建民等(2008)的研究表明气枪震源激发的地震波具有丰富的低频能量，其优势频率范围为4～6 Hz。因此，本研究选4～6 Hz带宽进行带通滤波。图4c为2013年1月7日21时59分53264台站接收到的气枪记录经过滤波预处理后的波形。

2.2.3反卷积

为了对齐到时达到压缩子波、消除震源影响的目的，笔者对数据进行反卷积运算求其格林函数。作反卷积前，还需要对经过去均值、去趋势、滤波预处理后的数据进行筛选，一方面剔除不正常、有断记或记录到天然地震的波形，另一方面剔除信噪比低的记录。这里我们将水库边布设的参考台接收记录作为子波(图5a)，用远台接收到的相同时刻的气枪记录(图5b)与参考台对应记录作反卷积求其格林函数(图5c)，由图5可见，消除震源影响后震相更清晰。

图5　反卷积过程

2.2.4互相关时延检测

重复上述数据预处理过程，求得处理时段内各台站每次激发记录的格林函数，分别将该时段内各台站的所有格林函数叠加作为进行各台互相关计算的参考模板，例如处理53264台2013年以来的数据，得到1 604条记录的格林函数，将这1 604条记录叠加作为53264台站的相关检测参考模板(图6a)。为了进一步提高信噪比，可根据具体情况将每天、每周、或每月的格林函数叠加，本研究选择叠加每天的格林函数，图6b为53264台2013年1月7日14条记录的格林函数叠加结果；最后将每一天的叠加结果与参考模板作互相关计算。具体做法如下：从零时刻开始以一定窗长从两波形信号中选取出相似波形窗口，然后保持一窗口不动，以不同时延移动另一窗口，并计算不同时延情况下两相似波形窗口的互相关系数(图6c)，互相关系数Cm(t)最大时所对应的时间延迟就是这两段相似波形窗口的走时差t1(图6d)，再以一定步长向后整体移动两相似波形窗口，重复上述计算过程。为了提高精度，将t1作cos插值得到更高精度的走时差t2(图6e)，最后结合震相提取初至波的走时数据，要求选择相关系数最大、t1变化平稳的时间窗，如图6中矩形框所示。

3处理结果

用上述数据处理方法处理了2013年以来6个台站的记录，测得初至波10-5～10-2的相对走时变化，与已有的人工震源波速变化研究结果吻合。Furumoto和Tsuda(2001)用重复爆炸震源连续监测日本Kanto-Tokai地区P波速度随时间的变化，观测到的波速年相对变化量达到10-3；Silver等(2007)在California 的两口井中同样用压电超声波震源进行了波速测量实验，测得量级为10-6s 的波走时变化；Wang等(2008)在云南小江断裂带开展了主动震源监测介质波速变化实验，利用电动重锤并用尾波干涉时延检测方法测得10-3～10-2的相对波速变化，精度为10-4。

图7a为2013年年初至2014年10月底全时段的初至波走时变化曲线，但由于每年5～9月水位变化，不便于开展气枪激发实验造成缺数，图7b～d为2013年以来3个集中激发时段的地震波走时差曲线放大图。图7b显示，2013年上半年所选6个台走时差变化趋势基本一致，其中53278台变化幅度较其它5个台要稍大；图7c为2013-09～2014-05的走时变化曲线图，由图可见，在2013-09～2014-02期间53278台的走时变化趋势与其它5个台有所差异；2014年9月以来，6个台站的走时变化较平稳，集中在-0.01～0.01s之间变化，且变化趋势一致(图7d)。

4结论与讨论

对挑选的6个主动源接收台站的气枪信号，采用近场校正、互相关检测等技术，提取初至波走时差，初步结果与认识如下：

(1)本研究用互相关检测技术测得10-5～10-2的相对波速变化，对比可知，波速变化在已有相关研究得到的波速变化范围内。

图7　6个台站不同时段的初至波走时变化曲线

(2)2013-01～2014-05的初至波走时变化较2014年下半年的大，但变化趋势基本一致；2014年9月初至10月底，初至波走时变化较平稳，集中在-0.01～0.01 s之间变化，且6个台的变化趋势一致。

(3)用互相关检测技术能精确测定地震波走时变化，短期内能测量到10-5～10-4的相对波速变化，但是长期的相对波速变化范围较大(10-5～10-2)，这可能是由于在长时间监测过程中地质体或实验条件的突变引起地震波形畸变，用来互相关计算的随后波形与参考波形相似性降低所造成。因此，要深入分析走时变化与地下介质应力状态的关系还需要充分考虑地质体、实验条件等(包括水库水位、接收台站气温、气压等)诸多因素的影响。

本文在撰写过程中得到中国地震局地球物理研究所王宝善研究员；云南省地震局苏有锦研究员、王彬副局长；中国地震局滇西地震预报实验场金明培高级工程师、杨军工程师、李孝宾助理工程师、叶泵工程师的指导和帮助，本文使用了中国地震局地球物理研究所王宝善研究员提供的互相关时延计算程序，在此表示衷心感谢。

参考文献：

林建民，王宝善，葛洪魁，等.2008.大容量气枪震源特征及地震波传播的震相分析[J].地球物理学报，51(1)：206-211.

罗桂纯，葛洪魁，王宝善，等.2008.利用相关检测进行地震波速变化精确测量研究进展[J].地球物理学进展，23(1)：56-62.

唐渊，刘俊来.2010.川滇西部上新世以来构造地貌：断裂控制的盆地发育及对于远程陆内构造过程的约束[J].岩石学报，26(6)：1925-1934.

王彬，杨润海，王宝善，等.2012.地震波走时变化精确测量的实验研究[J].云南大学学报(自然科学版)，34(S2)：15-20.

徐荟.2013.利用主动震源初至波研究小江断裂带浅层波速变化[D].北京：中国地震局地球物理研究所，23-38.

张建国，汪良谋，徐煜坚，等.1993.红河断裂深部震源环境介质力学性质分析[J].地震地质，15(2)：131-137.

Azizul H.Q..1981.An Overview on the time delay estimate in active and passive syste ms for target localization[J].Ieee Transations on Acoustics，Speech，and Signal Processing，229(3)：527-533.

Clifford G..1987.Coherence and time delay estimation[J].Proceedings of the Ieee，75(2)：236-255.

Furumoto J.，Tsuda T..2001.Characteristics of Energy Dissipation Rate and Effect of Humidity on Turbulence Echo Power Revealed by MU radar-RASS Measurements[J].J.Atmos.Solar-Terr.Phys.，63(2-3)：285-294.

Ronen S..2002.Psi，pascal，bars，and decibels[J].Lead Edge，21(1)：60-62.

Silver P.G.，Daley T.M.，Niu F.，etal..2007.Active source monitoring of cross well seismic travel time for stress induced changes[J].Bull.Seism.Soc.Amen，97(1B)：281-293.

Snieder R..2006.The theory of coda wave interferometry[J].Pure Appl.Geophys，163：455-473.

Stehly L.，Campillon M.，Shapiro N.M..2007.Travel time measurements from noise correlation：stability and detection of instrumental time-shifts[J].Geophys.J.Int.，171(1)：223-230.

Takeshi N.，Satoru T.，Teruo Y..2005.Temporal changes in seismic velocity of t he crust around Iwate volcano，Japan，as inferred from analyses of repeated active seismic experiment data from 1998 to 2003[J].Earth Planets Space，57(6)：491-505.

Wang B.S.，Zhu P.，Chen Y.，etal..2008.Continuous subsurface velocity measurement with coda wave interferometry[J].J.Geophys. Res.， 113(B12313)：1-12.

Study on Analysis Method of Travel Time Variations of Seismic

Wave of Active Source in Binchuan

LIU Zi-fen1，SU You-jin1，WANG Bao-shan2，WANG Bin1，YANG Jun1，LI Xiao-bin1

(1. Earthquake Administration of Yunnan Province，Kunming 650224，Yunnan，China)

(2. Institute of Geophysics，CEA，Beijing 100081，China)

Abstract

Firstly，we chose the complete data with high signal-to-noise ratio recorded by six receiving stations of active source in Binchuan County since 2013. Secondly，we measured the velocity change with high precision by using cross-correlation delay detection，near-field correction technique etc.. Finally，the conclusions are obtained as follow：(1)The relative velocity variation wasδv/v≈10-5～10-2in the medium beneath the stations. (2)From early Sep.，2014 to the end of Oct. ，2014，the travel time variations of initial wave was relatively stable which was between -0.01 s and 0.01 s，and the change trend of waveform recorded by six selected stations was same.

Key words：active source；cross correlation delay technique；airgun signal；travel time variation

中图分类号：P313.22

文献标识码：A

文章编号：1000-0666(2015)04-0591-07

*收稿日期：2015-04-07.基金项目：云南省青年地震科学基金项目(201303)、2015年度震情跟踪定向任务(2015010107)和云南省陈颙院士工作站联合资助.