2. 中国地震局兰州地震研究所, 兰州市东岗西路450号 730000
2. Lanzhou Institute of Seismology, China Earthquake Administration, Lanzhou 730000, China
地震的孕育和发生往往伴随着地下介质应力状态的变化,而介质应力状态的变化会引起地震波波速的变化(Birch,1960、1961;Scholz,1968;Nur et al,1969),通过分析跟踪特定区域的波速变化,识别大震发生前具有明确物理意义的前兆异常信息,对于地震预测研究等具有重要意义。地下介质应力变化引起的波速变化很小,其动态监测对观测系统精度的要求很高,不仅需要有重复性和稳定性较好的震源,还需要高精度的数据记录系统(王伟涛等,2009)。天然地震定位精度有限,且分布不均匀,因此利用天然震源测量地下介质波速变化的时空精度不太高。人工震源的激发位置已知,激发可控,可重复性好,探测精度高,但与天然地震相比,人工震源激发的能量低,探测距离有限,无论是炸药爆破、锤击震源还是列车震源,其能量都比天然地震小(罗桂纯等,2006)。而高性能人工震源和高精度测量系统的发展,则为地震波波速变化的高精度检测提供了条件(杨微等,2010)。
近年来,利用人工震源向地下发射地震波,进行地下介质性质的四维地震观测已成为各国地球物理学家关注的研究方向。早在20世纪70年代,Reasenberg等(1974)就尝试利用主动震源监测地下介质的变化。目前,小尺度(数米至数十米)的高精度主动源监测已取得了重要进展。日本学者(Yamamura et al,2003)利用压电陶瓷震源,基于超声波监测的方法发现了由固体潮等引起的地下介质的波速变化。Silver等(2007)发展了一套井间高精度监测系统,Niu等(2008)将其应用于SAFORD钻井中,发现数千米外2个小地震前出现了明显的波速变化。Wegler等(2006)利用气枪震源进行了基线长度为3km的观测,研究了火山爆发前介质S波波速的变化。由此可见,利用高性能人工震源主动向地下发射地震波,通过高度同步的地震台站接收地震信号是实现高精度监测地下介质变化的可行手段(陈顒等,2007;王宝善等,2011)。初步研究结果显示,大容量气枪震源在区域尺度地震观测中具有独特的优势:①能够产生足够能量的低频(数Hz)地震信号;②震源可高度重复;③震源对近场破坏小,是进行介质性质监测的理想震源(陈颙等,2007)。
地震波走时是与波速密切相关的参数。本文基于新疆呼图壁气枪震源发射台激发的气枪信号,以距气枪发射台108km的2015年2月22日沙湾MS5.0地震、153km的2015年12月6日沙湾MS4.8地震为震例,利用震中附近石场、石河子和乌苏等台记录的气枪信号,选取反射P波震相为时间窗(即震相窗),初步研究了2次5级左右地震前的走时变化。
1 新疆呼图壁气枪震源发射台简介2013年6月新疆主动震源野外科学观测研究站投入运行,这标志着中国第1座基于人工水体的大容量气枪发射台站顺利建成。气枪发射台位于距乌鲁木齐约80km的北天山北麓呼图壁县境内的古河滩荒地,激发池近圆形,上直径100m,下直径20m,深度18m。发射台配备6条2000 in3的Bolt 1500LL型气枪,2台大型空气压缩机,气枪沉放深度为10m,工作压力为15MPa,可连续每15min激发1次,6条枪同时1次激发释放的能量与1次ML0.9天然地震释放的能量相当(苏金波等,2015;杨微等,2013)。2013年8月开始持续激发实验,每周4晚至周5早上连续激发,每20min激发1次,截至2016年2月29日积累了4000余次的激发数据。从前期实验来看,在距激发点380km的基岩台可识别到数次连续激发实验叠加后产生的气枪震源的地震波信号(苏金波等,2015)。
2 研究思路地震波在传播过程中,高频信号衰减快,传播距离近,若想传播距离远,则需要水平传播距离足够远且能量足够强的低频信号。气枪震源作为表面源,是通过枪体内高压气体在水中瞬间释放来激发地震波的,震源信号的低频成分源于气泡脉冲,这种低频地震波可以在地下介质中传播较远,非常适合区域性地下结构的探测和地下介质动态变化的监测。新疆呼图壁气枪信号的主要频率集中在2~6Hz(苏金波等,2015),气枪信号经过叠加等处理后,在最远380km处的台站可以识别到气枪信号,但并不是在380km范围内的所有台站都可以清晰识别出信号,这主要是受台站的台基、信噪比和传播路径等的影响。
新疆主动源科研团队的初步研究表明,气枪震源与天然地震二者激发出的体波的性质存在差异。新疆呼图壁气枪震源激发出的S波信号不强,但P波信号较强,与天然地震的震源特性截然不同。天然地震震源具有一定深度,Pg震相是来自结晶基底以下弱梯度层的折射波,Pn震相是上地幔顶部弱速度梯度层的折射波(林建民等,2008),其大部分在密度较大的岩层中传播,波速较快;而气枪震源是表面源,Pg大部分在松散的土层中传播,速度较慢。但是,传播到一定距离后,特别是首波出现后,天然地震和气枪震源的信号走时越来越接近,说明气枪震源激发的地震波信号与天然地震在传播路径上存在着一定的相似性。据此,本研究选取P波反射震相作为走时获取对象。由于反射波在特定界面以上的介质中传播,携带丰富的介质信息,若反射波在震中台站附近出射,则很大程度上是穿越震源区的,这可能对研究地震前孕震区介质特性的异常变化等有益。因此,我们以反射波震相走时为数据来源分析2次震例,以期发现震前走时变化的异常特征。
3 气枪震源反射波震相识别 3.1 地壳分层及速度模型选取要确定气枪震源信号的不同震相信息,需要找到合适的地壳速度模型。速度模型一直是地震学研究中的经典问题,全球的分层速度模型包括IASPEI91、AK135等模型,地壳速度模型与地震定位的精度密切相关。
前苏联И.Л.涅尔赛索夫等1964年发表了“3400km走时表”(以下简称“3400走时表”),其测线从帕米尔至贝加尔,全长3500km,穿过中亚、哈萨克斯坦、阿勒泰萨彦和贝加尔沿岸等几个地区(新疆地震局分析预报室,1982;唐明帅等,2010;孙安辉等,2011)。“3400走时表”的西部走时表也部分适用于新疆地区,“3400走时表”是一个大范围内的平均走时,震中距愈大,符合愈好,而在小区域范围内,结果偏差稍大。陈向军等(2014)、王俊等(2014)对“3400走时表”拟合得到了符合该走时表的速度结构模型,解析后发现,其速度结构为1个固定震源深度10km、双层结构的走时表(表 1),将其进一步转化为对应的单层速度模型(表 2)。由表 1、2可见,在不同地壳模型结构下地震波走时存在差异。
双层模型和单层模型下地震波的传播路径如图 1所示。地震反射波是指由震源出发,经过界面反射后的波,在近震范围内,记录较好的反射波是莫霍面反射波PmP和SmS,该波通常在震中距60~120km的范围内较为明显(陈向军等,2014)。双层地壳模型时,P波和S波在康纳德面外侧的反射波为Pc和Sc(刘瑞丰等,2014)。
根据单、双层地壳模型,设定震源深度为0km,计算得到新疆地区地震反射波震相的理论走时曲线(图 2)。
依据新疆地区地震波P波反射震相理论走时的计算结果,本文中的3个台站与气枪发射台的距离分别为:89(石河子台)、112(石场台)、186km(乌苏台)。双层模型计算得到的PmP理论走时分别为:23.45、25.91、35.35s;单层模型计算得到的PmP理论走时分别为:23.44、25.9、35.35s,单、双层地壳模型得到的理论走时差异很小。图 3给出了3个台经数据处理后记录到的气枪信号,由图 3可见,选取PmP震相的理论走时作为参考走时,在实际记录波形上可以识别出PmP震相,且实际走时与理论走时基本一致,这样就为选取震相窗口提供了方便。因此,参考理论走时选取P波反射震相是可靠的。
自2013年6月以来,气枪震源每周进行1组连续激发,选择环境噪声相对较小的夜间为激发时间,每周可以获得超过30次激发的数据。由于2014年1月中旬至4月下旬,空压机出现故障,因而仅有1天的激发记录;2014年7月,气管出现大面积漏气,直至8月下旬修理完成,此间未进行激发试验;2015年8月初至9月中旬,参考台GPS授时系统出现故障,不能提取准确的零时,时间数据无法使用。除这3段时间以外,其他时间数据基本连续。此外,2014年6~7月,还进行了为期1个月的、以30min为时间间隔的连续激发。目前,已累计激发4000余次,积累了一定量的数据。
观测实验期间,距气枪发射台108km处发生2015年2月22日沙湾MS5.0地震和153km处发生2015年12月6日沙湾MS4.8地震,这2次5级左右地震震中相距约50km(图 4),本文以这2次地震为震例样本,研究地震前走时变化的时空演化特征。
选取距离2次5级左右地震约80km范围内的3个台站记录的2013年8月~2016年2月气枪震源数据,3个台站分别为距气枪震源89km的石河子台、112km的石场台和186km的乌苏台(图 4),其台基噪声分别为Ⅲ、Ⅰ、Ⅰ级,考虑到台站台基及接收信号的能力,选择石河子台、石场台2周接收到的气枪震源数据进行叠加等处理来选取震相窗,而对于乌苏台则将1周接收的数据进行叠加。
4.2 原理尾波干涉法基于路径叠加理论,其认为波场是沿所有可能散射路径的波的叠加(Snieder,2006),通过对比扰动前、后时间窗口内的2个波形的走时偏移分布规律,进而求得波速变化。在测量介质波速变化方面,尾波干涉法精度很高,应用广泛(张金川等,2014)。利用互相关技术提取气枪激发时刻时,相关检测原理的一项重要应用就是测量地震波速的变化,相关检测技术则是利用信号周期性和噪声随机性的特点,通过自相关或者互相关运算,达到去除噪声的目的。计算不同时刻2个波形信号的互相关函数,可以获取2个信号的延时,从而获得介质速度的变化。如果介质中的波速变化是均匀的,则介质中的相对波速变化可以表示为
$ \frac{{{\rm{d}}v}}{v} = - \frac{{{\rm{d}}t}}{t} = {\rm{const}} $ | (1) |
其中,v为介质中的波速;dv为介质中的波速变化;t为尾波某部分的到时;dt为尾波该部分的到时变化。由式(1) 可见,尾波某时刻的到时变化dt随到时t线性变化,通过测量不同时刻到时变化dt随到时t变化的斜率,就可得到介质中的波速变化(Niu et al,2003;Wang et al,2008)。具体做法是,在2个相似波形中选取适当窗长,选取1个窗口为参考窗口,移动另1个波形窗口,计算互相关系数,互相关系数最大时对应的时间延迟即为2个波形信号的走时差。
4.3 数据处理 4.3.1 利用参考台数据提取激发时刻将激发池边的参考台数据作为提取激发时刻的源数据时,波形数据会存在非0的均值或存在长周期的线性趋势,为了不影响数据分析,必须在数据分析前进行去均值、去趋势等处理。因新疆呼图壁主动源气枪信号的主频为2~6Hz,故首先对波形信号进行2~6Hz的4阶Butterworth带通滤波,以降低噪声干扰。然后,选取参考波形窗口,利用互相关检测方法提取激发时刻及截取实验波形信号,同时需要手动剔除记录信号错误和信噪比较低的个别波形。
4.3.2 固定台数据处理依据参考台提取到的激发时刻,对选取的3个固定台记录波形截取记录信号时段波形。通常在信号处理之前,首先需要对信号进行去均值、去趋势、滤波及剔除质量差的数据等操作;之后,将激发池边的参考台接收的气枪信号作为参考波形,3个固定台站接收的相同时刻的气枪信号与参考台相应信号作反卷积得到格林函数,以消除震源的影响并获取清晰的震相。最后,将3个台得到的格林函数分别进行全部叠加作为参考模板。因3个台站距离气枪震源较远,1次激发信号的接收范围有限,并且受接收台站的台基、背景噪声水平等的影响,若要清晰地识别气枪信号,需要多条记录叠加。因此,乌苏台按1天数据记录进行叠加,得到相应的格林函数,石河子台和石场台按14天窗长、7天步长滑动叠加得到相应的格林函数,将叠加得到的每个格林函数与参考模板作互相关计算,获取相应震相窗的走时变化。
5 地震前后地震波走时变化分析利用前述原理处理了2013年8月~2016年2月底乌苏台、石河子台和石场台记录的数据,依据“3400走时表”得到震源深度为0km的理论PmP震相走时,选取震相窗长为0.2s的完整周期,每个台选取3组震相窗获取走时变化,重点分析2次5级左右地震前走时的变化。
图 5(a)为石河子台的3组PmP震相窗走时变化曲线。依据“3400走时表”,石河子台距主动源89km,距2015年2月22日沙湾5.0级地震震中约19km,距2015年12月6日沙湾4.8级地震震中约68km。在双层模型下,PmP理论走时为23.45s;而在单层模型下,PmP理论走时为23.44s,2种模型差异不大,故结合记录波形,选取了21~24s之间的3组震相窗(图 5(b))。沙湾5.0级地震前,2014年10月24日走时变化开始下降,低值异常持续近4个月,走时变化幅度最大达0.06s,相对变化幅度达0.26%,其后发生了沙湾5.0级地震。而沙湾4.8级地震前,石河子台记录的走时变化自2015年9月25日起出现低值起伏波动,持续近2个月之后发生了沙湾4.8级地震,走时变化幅度达0.03s,相对变化幅度达0.12%,但震后走时变化依然低值,说明震前波速增加,区域应力积累,震后应力释放并不完全,导致波速下降缓慢。
图 6(a)为石场台的3组PmP震相窗走时变化曲线。依据“3400走时表”,石场台距主动源112km,距2015年2月22日沙湾5.0级地震震中约22km,距2015年12月6日沙湾4.8级地震震中约42km。在双层模型下,PmP理论走时为25.91s,比单层模型下理论走时增加0.01s,参考理论走时,选取了23.7~26.0s之间的3组震相窗(图 6(b)),实际记录波形震相也可以清晰地看到PmP震相到时。沙湾5.0地震前,2014年10月17日走时变化开始下降,低值异常持续近4个月,走时变化幅度最大达0.057s,相对变化幅度达0.21%,随后发生了沙湾5.0级地震。而沙湾4.8级地震前石场台记录走时变化低值异常约2个月,走时变化幅度达0.053s,走时相对变化幅度最大达0.2%左右,走时变幅较大可能是由于石场台距本次地震震中较近。
图 7(a)为乌苏台的3组PmP震相窗走时变化曲线。依据“3400走时表”,乌苏台距主动源186km,距2015年2月22日沙湾5.0地震震中约80km,距2015年12月6日沙湾4.8级地震震中约54km。在双层模型和单层模型下,PmP理论走时均为35.35s,选取了34.6~36.0s之间的3组震相窗(图 7(b))。沙湾5.0地震前,2014年10月24日走时变化开始下降,低值起伏波动长达4个月,走时变化最大幅度约0.051s,相对变化幅度达0.15%。而沙湾4.8级地震前走时变化低值异常持续近2个月,走时相对变化幅度最大达0.1%。
2次5级左右地震前,3个台记录的PmP走时变化曲线均出现了低值异常。沙湾5.0级地震前,石河子台和石场台相对走时变化幅度达0.2%左右,而乌苏台相对走时变化幅度较小,这可能与台站所处位置有关。而沙湾4.8级地震前,3个台记录走时变化趋势整体较小,这可能与地震震级大小有关,孕震区域应力变化存在的差异,导致相对走时变化幅度较小。
6 结论与讨论本文基于尾波干涉法并利用互相关检测等技术,对3个台站接收到的主动源信号选取PmP震相窗以获取走时变化,通过分析沙湾2次5级左右地震前走时变化特征,得到如下初步认识:
(1) 利用P波反射震相走时作为震相窗选取的依据,因反射波携带了沉积层、上中地壳和壳幔边界等丰富的介质变化信息,可能有利于研究震源区介质的异常特性。
(2) “3400走时表”的西部走时表也部分适用于新疆地区,根据单、双层速度模型,参考理论PmP震相走时,从石河子台、石场台和乌苏台记录的气枪震源信号中识别到PmP震相,表明参考理论PmP走时选取震相窗口是可靠的。
(3) 取2013年8月~2016年2月实验数据,利用互相关检测等技术选取3组震相窗,获得每个台的走时变化曲线,分析2015年2月22日沙湾5.0级地震、2015年12月6日沙湾4.8级地震前地震波走时变化特征。3个台站获取的PmP走时变化曲线显示,2次5级左右地震前均出现了不同程度的低值异常,沙湾5.0级地震前均出现了近4个月的低值异常,且走时变化幅度较大;而沙湾4.8级地震前均出现近2个月的低值异常,走时变化幅度较小。2次5级左右地震前石河子台和石场台获取的走时异常变化特征显著于乌苏台。
(4) 2次5级左右地震前,3个台站获取的PmP走时变化存在不同程度的低值异常特征,这可能与地震震级等因素有关。一般来说,震级小,释放能量少,波速变化就小(杨微等,2010)。本研究中沙湾5.0级地震前走时变化幅度比沙湾4.8级地震前相对较大,可能与震级大小有关。徐荟等的(2015)研究也发现,距震源较近的道相波速变化比较远的道相更剧烈,石场台、石河子台距2次5级左右地震震中较近,PmP反射震相的传播路径可能经过2次地震的震源区后被石河子、石场台接收,而乌苏台距2次地震震中较远,PmP反射震相可能穿越震源区较少,所以走时异常变化特征没有其他2个台显著。3个台记录走时相对变化幅度为0.10%~0.26%,与王伟涛等(2009)观测到的2%的相对变化幅度相比偏小,但与以往研究(杨微等,2010;Ikuta et al,2004)得出的平均波速变化幅度0.1%~0.4%比较则相近;Wang等(2008)在云南开展了主动源监测介质波速变化实验,并用尾波干涉时延检测方法得到10-3~10-2的相对变化,精度为10-4;徐荟等(2015)利用主动源也得到了10-3~10-2的相对波速变化;刘自凤等(2015)利用云南气枪震源得到了初至波10-5~10-2相对走时变化;本文利用新疆气枪震源可以提取到10-3~10-2相对变化,与前人的研究结果相近。地震发生前,断裂带附近应力逐渐积累,波速增加;相反,地震后应力释放,波速相对下降,地震后一段时间应力调整,波速又相对增加(杨微等,2010)。因此,若要准确分析波速变化与地下介质应力变化间的关系还需充分考虑地质体、实验环境等因素的影响。综上所述,获取高精度的波速变化,不仅需要高精度的观测系统,也需要考量其它因素的影响,这样才有望提取到与地震相关的波速变化,才更加有利于我们认识地震的孕育过程。
(5) 利用气枪震源研究介质波速变化工作仍处于探索阶段,研究比较初步,震例选取较少,今后还需要积累更多的震例以待进一步深入研究。
致谢: 本文在撰写过程中得到了中国地震局地球物理研究所冀战波博士、新疆地震局向元助理工程师的大力支持和帮助,本文使用了中国地震局地球物理研究所王宝善研究员提供的互相关时延计算程序,在此表示衷心的感谢。陈向军, 上官文明, 宋秀青, 等. 2014, 新疆全区和分区地壳速度模型的分析. 中国地震, 30(2): 178–187. |
陈顒, 王宝善, 葛洪魁, 等. 2007, 建立地震发射台的建议. 地球科学发展, 22(5): 441–446. |
罗桂纯, 王宝善, 葛洪魁, 等. 2006, 气枪震源在地球深部结构探测中的应用研究进展. 地球物理学进展, 21(2): 400–407. |
林建民, 王宝善, 葛洪魁, 等. 2008, 大容量气枪震源特征及地震波传播的震相分析. 地球物理学报, 51(1): 206–212. |
刘瑞丰, 陈翔, 沈道康, 等. 2014, 宽频带数字地震记录震相分析: 26–27. |
刘自凤, 苏有锦, 王宝善, 等. 2015, 宾川主动源地震波走时变化分析方法研究. 地震研究, 38(4): 591–597. |
苏金波, 王宝善, 王海涛, 等. 2015, 利用大容量气枪震源资料研究北天山地区介质衰减特征. 地震研究, 38(4): 598–605. |
孙安辉, 陈棋福, 陈颐. 2011, 天山东北部地震的重新定位和一维地壳速度模型的改善. 中国地震, 27(3): 235–246. |
唐明帅, 王海涛, 段天山. 2010, 利用和田地震台阵数据对2008年于田7.3级地震序列重新定位. 内陆地震, 24(3): 227–235. |
王宝善, 王伟涛, 葛洪魁, 等. 2011, 人工震源地下介质变化动态监测. 地球科学进展, 26(3): 249–256. |
王俊, 宋秀青, 陈向军, 等. 2014, 新疆于田M7.3级主震精定位研究. 中国地震, 30(2): 188–197. |
王伟涛, 王宝善, 葛洪魁, 等. 2009, 利用主动源检测汶川地震余震引起的浅层波速变化. 中国地震, 25(3): 223–233. |
徐荟, 刘学军, 王彬, 等. 2015, 利用主动震源直达波互相关时延检测技术监测小江断裂带浅层地震波波速变化. 地震研究, 38(1): 7–15. |
新疆地震局分析预报室. 1982, "3400km走时表"及其应用、 新疆走时表工作的设想. 地震地磁观测与研究, 3(4): 25–30. |
张金川, 王勤彩, 薛兵, 等. 2014, 用尾波干涉法监测介质波速变化研究进展. 地震, 34(3): 62–73. |
杨微, 葛洪魁, 王宝善, 等. 2010, 由精密控制人工震源观测到的绵竹5.6级地震前后波速变化. 地球物理学报, 53(5): 1149–1157. |
杨微, 王宝善, 葛洪魁, 等. 2013, 大容量气枪震源主动探测技术系统及实验研究. 中国地震, 29(4): 399–410. |
Birth F. 1960, The velocity of compressional waves in rocks to 10 kilobars.part 1. J Geophys Res, 65: 1083–1102. DOI:10.1029/JZ065i004p01083. |
Birth F. 1961, The velocity of compressional waves in rocks to 10 kilobars.part 2. J Geophys Res, 66: 2199–2224. DOI:10.1029/JZ066i007p02199. |
Ikuta R, Yamaoka K. 2004, Temporal variation in the shear wave anisotropy detected using the Accurately Controlled Routinely Operated Signal System. J Geophys Res, 109: B09305. DOI:10.1029/2003JB002901. |
Silver P G, Daley T M, Niu F, et al. 2007, Active source monitoring of cross-well seismic travel time for stress-induced changes. Bull Seism Soc Am, 97(1B): 281–293. DOI:10.1785/0120060120. |
Niu F L, Silver P G, Daley T M, et al. 2008, Preseismic velocity changes observed from active source monitoring at the Parkfield SAFOD drill site. Nature, 454: 204–208. DOI:10.1038/nature07111. |
Nur A, Simmons G. 1969, The effect of saturation on velocity in low porosity rocks. Earth and Planetary Science Letters, 7(2): 183–193. DOI:10.1016/0012-821X(69)90035-1. |
Reasenberg P, Aki K. 1974, A precise continuous measurement of seismic velocity for monitoring in situ stress. J Geophys Res, 79(2): 399–406. DOI:10.1029/JB079i002p00399. |
Scholz C H. 1968, Microfracturing and the inelastic defomation of rock in compression. J Geophys Res, 73: 1417–1432. DOI:10.1029/JB073i004p01417. |
Snieder R. 2006, The Theory of Coda Wave Interferom etry. Pure Appl Geophys, 163(2-3): 455–473. DOI:10.1007/s00024-005-0026-6. |
Wegler U, Lühr B G, Snieder R, et al. 2006, Increase of shear wave velocity before the 1998 eruption of Merapi volcano(Indonesia). Geophys Res Lett, 33(9): L09303. |
Wang B S, Zhu P, Chen Y, et al. 2008, Continuous subsurface velocity measurement with coda wave interferometry. J Geophys Res, 113(B12): B12312. DOI:10.1029/2007JB005154. |
Yamamura K, Sano O, Utada H, et al. 2003, Long-term observation of in situ seismic velocity and attenuation. J Geophys Res, 108(B6): 2317. |