2012, Vol. 55
很多自然事件在孕育、发生和发展的过程中均有可能向大气中辐射低频的次声波，地震也不例外.这些次声波由于频率很低，可以在大气中远距离传播[1-2].因此，建立一个广域的次声监测网具有极其重要的实际意义.利用这个网络，我们可以对地震等自然事件产生的次声波进行监测，并通过对次声波信号的分析，揭示这些自然事件产生次声波的机理，研究传播规律，从而为有效预测这些自然事件提供有价值的信息.
鉴于此目的，国内外很多学者长期致力于自然事件次声信号的分析和研究.Milton A G 等从20世纪90年代开始对火山爆发产生的次声信号进行研究[3]，并在1999 年通过对日本Sakurajima火山爆发之前的次声信号的记录和分析，研究了火山爆发前兆次声波的相关特性[4].AssinkJD 和EversLG 等[5]在2008年研究了闪电产生次声波的相关特征.他们通过分析距离闪电区域50km 之外次声阵所接收到的信号，结果表明:闪电产生的次声波具有爆炸波的特征，其主频在1~5Hz之间.与地震相关联的次声信号的研究也是由来已久，早在1964年，BruceAB[6]就分析了当时Alaska大地震所产生的一种大气波，其研究表明:这种波是由地震所引起的震中附近地表突发性地垂直起伏和倾斜而产生的，可以称作地震空气波.同样，这种机制也可以向大气中辐射次声波[7].20世纪90年代中期，基于全面禁止核试验条约组织(CTBTO)在全球逐步建立起的60多个次声监测站，越来越多的与地震相关的次声波被观测和研究.MikumoT 等[8]观测到了产生于2004年Sumatra-Andaman9.2 级大地震的大气低频声-重力波，其周期范围是350~715s之间，波的群速度达到300~314 m/s.同时，他们假设了几种地震源参数，利用真实大气模型得到了与观测信号一致的合成波形，从而研究了信号源的相关特征.LePichon等[9]观测到了相同地震发生时来自于震中以及大陆块附近的次声波.随着地震次声信号观测研究的发展，一系列次声信号处理方法也得以建立和发展，例如PMCC[10-12]，F-K 分析法[13]等.但是，CTBTO次声站点主要用于监测的周期在0.05~100s范围内的由核爆炸产生的次声波，而在地震发生时，特别是地震发生前辐射的一些次声波的周期要大于100s[14-15].因此，为了更好地探测频率较低的次声波，中国科学院声学研究所设计并制造了InSYS2008型低频宽带次声测量传感器.该传感器可测量的次声波周期范围为10-4~102s, 适合于探测地震发生前或发生时所产生的次声信号.同时，在华北地区建立了一个广域次声监测网，阵元信号的数据都是在线实时传输，可用于实时探测低频的次声波信号.
本文介绍了发生于2011年10 月12 日北京海淀区的一次震级为1.8 级的小型地震前4 天(2011年10月8日)由北京地区的五个次声监测站所得到的异常次声信号，并对这五路信号的时域波形、时频特性以及信号间的相关性进行了分析，最后对声源位置进行了成像定位研究.
2 广域次声监测网为了监测和研究与自然事件相关联的低频次声波，中国科学院声学研究所设计并制造了InSYS2008型低频宽带次声测量传感器.该传感器测量的频带范围为0.0001~100 Hz, 动态范围为128 dB(0.00002~500Pa)，一致性误差小于1dB，温度漂移为0.01mV/℃，抗电磁干扰能力强，可实现对次声波长期不间断的测量.图 1 是该传感器和用于传递数据的网络传输仪的实物图.图 2 为传感器的灵敏度曲线.从图 2 中可以看到，该传感器的灵敏度高，频带宽，在1000s 的周期处，依然可以达到30mV/Pa以上的幅值.
|
图 1 低频宽带次声传感器与网络传输仪 Fig. 1 Broad-band and low-frequency sensor with network transmission device |
|
图 2 低频宽带次声传感器灵敏度曲线 Fig. 2 Infrasound sensors' sensitivity curve |
同时，声学所计划在全国范围内布设次声传感器，建立一个广域次声监测网.自2011年7月以来，先后在华北地区建立次声监测站点17个，其中北京地区设有5个监测站.图 3 是目前为止广域次声监测网的站点分布图，其中黑色实心圈代表着布设站点.
|
图 3 广域次声监测网站点分布图(北京设有5个站点) Fig. 3 Distribution of stations of the monitoring infrasound network, including 5 stations in Beijing |
2011年10月12日16时40分，在北京市海淀区(40°0'0″N，116°12'0″E)发生了里氏1.8级的小型地震，震源深度为5.1km.图 4给出了中国地震台网所给出的地震发生地点与北京地区的五个监测站点的相对位置.其中A-E 五个黑点所标识的是监测站点的位置，其经纬度依次为:A:39°59'38″N，116°25'23″E；B:39°59'30″N，116°12'39″E；C:39°59'19″N，116°19'44″E；D:40°3'44″N，116°17'25″E；E:39°49'52″N，116°21'07″E.红点标识的是地震发生位置，其靠近B点，位于B点的西北方向.由图 4可以看到，监测站点距离震源并不远.因此，虽然该地震的震级很小，但是传感器依然清晰地记录下相关信号.
|
图 4 海淀地震发生位置与北京地区五个次声监测站点的分布图 Fig. 4 Epicenter of the Haidian earthquake and locations of five stations for monitoring infrasound in Beijing |
图 5给出的是上述五个站点于2011年10月8日16:03:00到18∶49∶40共10000s时间内所记录到的特征次声信号.从时域波形图中可以看到，这五路信号的时域波形非常相似，都为“N"形的脉冲波，脉冲波峰值最高能达到10Pa.这五路信号波形从波峰到波谷的“下降"部分都不同程度地受到了高频噪声的影响，E路信号受到的高频噪声影响最为明显.
|
图 5 2011年10月8日16时03分A-E五个站点接收的低频次声信号 Fig. 5 Waveforms of five infrasound signals detected by A-E stations in Beijing, at 16 :03 :00，8 0ct, 2011 |
这种次声信号属于非平稳信号，其频谱是随着时间变化的.因此，在分析和处理这类信号时，单纯的Fourier分析方法已经不再适用.在此，采用了Wigner-Ville分布时频分析方法对这类信号频谱特性进行了研究.信号S(t)的Wigner-Ville 分布定义为[16]:
|
(1) |
这里，z(t)为S(t)的解析信号，可由S(t)的Hilbert变换而得到，即:
|
(2) |
对于离散时间信号x(n)，其Wigner-Ville分布可以表示为:
|
(3) |
图 6是图 5所示五路信号的时频分布图.从图 6中可以发现，这五路信号的能量都主要集中在0.025Hz以下，波形的波峰和波谷处，而高频噪声的能量很弱.并且，波峰区域的高能量持续时间要长于波谷区域高能量的持续时间.因此，根据时频分析的结果将它们进行0.025Hz以下的低通滤波处理.
|
图 6 A-E 五路信号的时频分布 Fig. 6 Spectrograms of the five infrasound signals detected by A一E stations |
图 7中的A-E 五幅图是低通滤波后的时域波形，其中标注的四个时间依次是信号的起始、波峰、波谷以及结束的时间.可以看到，经过滤波后五路信号信噪比高，五个波形相似性非常高.次声波最先到达的站点为D，到达时间为16∶39∶46，然后依次到达B、C、A、E 四个站点.A-E 中五个“N"型波的持续时间依次为55分41秒，57分30秒，54分56秒，57分14秒以及1小时0分59秒，持续时间基本是一致的，均在1 个小时左右.只有E 路信号的持续时间略长，这是因为该信号从波峰到波谷所经历的时间要更长，与上述E 路信号受到的高频噪声影响较强有关.
|
图 7 低通滤波后五路信号的时域波形 Fig. 7 Waveforms of five signals after low-pass filtering |
对于两路信号x(n)和y(n)，其信号间的相关系数定义为:
|
(4) |
图 8是根据(4)式而得到的B 路信号与其它四路信号的相关系数分布图.从图 8 中可以清晰地看到，B路信号与其它四路信号的相关系数分布的主瓣值均非常得高，且副瓣值相对于主瓣值要小很多.其中，A 和E 两路信号与B 路信号的相关系数最大值达到了0.8，而C 和D 两路信号与B 路信号的相关系数最大值均在0.8以上，达到了0.87.由于A、E 两个站点距离B 站点更远，并且这两个站点均位于市区内，因此次声波在传播过程中受到的背景噪声干扰要更大，从而造成了这两路信号在传播过程中更大的波形畸变.这也是A 和E 两路信号与B路信号的相关系数要略低于C 和D 两路信号与B 路信号相关系数的原因.五个信号之间的高相关性以及较高的主-旁瓣比值表明信号间的相似度很高，说明这一次声信号并非偶然被五个站点所接收到的，五路信号应该来自于同一波源.
|
图 8 A、C、D、E 四站点的次声信号与B站点次声信号的相关系数分布图 Fig. 8 Correlation coefficients between signals detected by B station and other four signals |
根据上述五路站点所接收到的高信噪比和高相似度的次声信号，利用波束形成方法对次声波源进行了成像定位研究.
首先，将阵列接收到的五路原始时域信号Sm(t)(m=1，2，…，5)进行低通滤波处理.然后，将滤波后的信号进行傅里叶变换，得到频域信号Sm(ω).接着，通过延迟求和波束形成方法得到所求区域范围(xi，yj)内的信号幅值:
|
(5) |
这里，i，j分别表示x，y方向所划分网格数.Δτm表示次声波到原点与其到五个站点的时延，这里选取B站点为原点.时延量由相关系数最大值所对应的时间值所决定，计算可得A，C，D，E 四路信号与B路信号间的时延量分别为:τBA =tA-tB =689.1s, τBC =tB -tC =314.2s, τBD =tD -tB =-284.2s, τBE =tE -tB =1410.9s.同时，根据时延量和站点间的距离求得波的传播速度约为25m/s.
图 9是实测次声波源的成像定位图.由图 9可以清晰地看到这五路信号来自阵列的西北方向，与震中方向基本相同.但是，波源位置处于西北偏北方向，更加靠近D站点，与震中位置有所偏移，偏移距离约为5km.
|
图 9 北京五点阵列实测次声波源的成像定位 Fig. 9 Imaging location of the infrasound source |
从图 7中可以看出，次声信号是最先到达D 站点的，这说明该次声波的声源位置更加靠近D 站点，与波源定位的结果是相符的.也就是说，该地震前异常次声信号源位置与地震发生时震中位置的方位基本相同，均位于阵列西北方向，但两点并不完全重合.文献[7]中指出，地震产生次声波有两种机制:(1)震中附近地面活动可能产生次声波，其声源位于震中附近；(2)地震表面波与特征地形(例如:山脉)的相互作用，通过衍射过程产生次声波，其声源位置距离震中要远.值得注意的是，此次地震发生地点靠近香山山脉，很有可能是第二种机制所产生的，所以，其声源与地震发生位置有些偏差也是合理的.
为了检验北京地区五点阵列对来自其西北方向的声波波源的定位能力，进行如下仿真研究.图 10是对该五点阵列定位特性进行的仿真计算结果.假设次声波源位于阵列西北方向，(-20km, 20km)处.根据监测信号的频率范围和传播速度，假设声源为频率为0.02Hz的单频正弦波信号，波的传播速度为25m/s.从图 10中可以看到，该阵列能够比较准确地对声源进行定位.由于波的频率很低以及阵元个数较少的原因，被定位区域的范围就显得比较大，但这并不影响阵列对次声波源方位的估计.
|
图 10 北京五点阵列对单点次声波源成像定位的仿真结果 Fig. 10 Imaging location simulation of a singleinfrasound source by the Beijing infrasound array composed of five stations |
本文通过对2011年10月12日北京的一次小震级地震发生之前所观测到的异常次声信号进行分析和研究后发现，五路信号的时域波形均为“N"形脉冲波，波形非常相似；时频分析的结果表明信号能量主要集中在0.025Hz以下，“N"波的波峰和波谷处.相比而言，高频噪声的能量很弱.低通滤波之后的五路信号的相关性很高---A 和E 两路信号与B路信号的相关系数最大值达到了0.8，而C 和D 两路信号与B 路信号的相关系数最大值均在0.8 以上，且副瓣值与主瓣值相比很低.由于此次接收的五个站点位于北京市区内，因此波传播过程中必然会受到环境噪声及地形等因素的影响，引起波形的部分畸变，从而造成了A、E 两路信号与B路信号的相关性要略差于C、D 两路信号与B路信号的相关性.五路信号之间的高相关性以及较高的主旁瓣比表明了信号间的高相似度，说明该次声信号并非偶然被接收到的大气波，这五路次声信号应该来自于同一波源.对次声波源进行了成像定位研究，其结果表明:震前异常次声信号源的方向与地震发生时的震中方向基本一致，附近山脉可能造成了声源位置与震中位置有所偏差.
该信号是在地震发生前4天接收到的，由于震级较小，仅为里氏1.8级，因此位于北京地区以外的站点并没有接收到明显的相关信号.尽管如此，依然可以表明地震之前可能会有异常次声信号的出现.未来将加大广域次声监测网的布设范围，继续加强对地震前后出现的次声波的观测研究，这将使我们得到更有规律性的认知，可望为地震预测技术提供更有价值的信息.
| [1] | Le Pichon A, Blanc E, Hauchecorne A. Infrasound Monitoring for Atmospheric Studies. New York: Springer. 2010 . |
| [2] | Evers L G. The inaudible symphony: on the detection and source identification of atmospheric infrasound. Netherlands: Delft University of Technology, 2008. |
| [3] | Garcés M A, Hansen R A. Waveform analysis of seismoacoustic signals radiated during the fall 1996 eruption of Pavlof volcano, Alaska. Geophys. Res. Lett. , 1998, 25(7): 1051-1054. DOI:10.1029/98GL00543 |
| [4] | Garcés M A, Iguchi M, Ishihara K, et al. Infrasonic precursors to a vulcanian eruption at Sakurajima volcano, Japan. Geophys. Res. Lett. , 1999, 26(16): 2537-2540. DOI:10.1029/1998GL005327 |
| [5] | Assink J D, Evers L G, Holleman I, et al. Characterization of infrasound from lightning. Geophys. Res. Lett. , 2008, 35(15): L15802. DOI:10.1029/2008GL034193 |
| [6] | Bolt B A. Seismic air waves from the great 1964 Alaskan earthquake. Nature , 1964, 202(4937): 1095-1096. DOI:10.1038/2021095a0 |
| [7] | Mutschlecner J P, Whitaker R W. Infrasound from earthquakes. J. Geophys. Res. , 2005, 110: D0108. DOI:10.1029/2004JD005067 |
| [8] | Mikumo T, Shibutani T, Le Pichon A, et al. Low-frequency acoustic-gravity waves from coseismic vertical deformation associated with the 2004 Sumatra-Andaman earthquake (Mw=9.2). J. Geophys. Res. , 2008, 113: B12402. DOI:10.1029/2008JB005710 |
| [9] | Le Pichon A, Herry P, Mialle P, et al. Infrasound associated with 2004-2005 large Sumatra earthquakes and tsunami. Geophys. Res. Lett. , 2005, 32(19): L19802. DOI:10.1029/2005GL023893 |
| [10] | Cansi Y. An automatic seismic event processing for detection and location: The P. M. C. C. method. Geophys. Res. Lett. , 1995, 22(9): 1021-1024. DOI:10.1029/95GL00468 |
| [11] | Havelock D, Kuwano S, Vorlander M. Handbook of Signal Processing in Acoustics (Vol.1). New York: Springer, 2008 . |
| [12] | Cansi Y, Klinger Y. An automated data processing method for mini-arrays. Newsletter of the European-Mediterranean Seismological Centre , 1997, 11: 2-4. |
| [13] | Smart E, Flinn E A. Fast frequency-wavenumber analysis and Fisher signal detection in real-time infrasonic array data processing. Geophys. J. Roy. Astr. Soc. , 1971, 26(1-4): 279-284. |
| [14] | 林琳, 杨亦春. 大气中一种低频次声波观测研究. 声学学报 , 2010, 35(2): 200–207. Lin L, Yang Y C. Observation & study of a kind of low-frequency atmospheric infrasonic waves. Acta Acustica (in Chinese) , 2010, 35(2): 200-207. |
| [15] | 胡心康, 石金瑞, 任建国. 大气中一种压力波动与地震关系的初步探索. 地球物理学报 , 1980, 23(4): 450–458. Hu X K, Shi J R, Ren J G. Preliminary study of the relation between a kind of pressure waves in the atmosphere and earthquakes. Chinese J. Geophys. (in Chinese) , 1980, 23(4): 450-458. |
| [16] | 张贤达, 保铮. 非平稳信号分析与处理. 北京: 国防工业出版社, 1998 . Zhang X D, Bao Z. Non-stationary Signal Analysis and Processing (in Chinese). Beijing: National Defence Industry Press, 1998 . |