Observation of Local Infrasound Coupled by Seismic Wave on Wide Spread Infrasound Network
-
摘要: 通过提出应用于广域次声传感器阵列的最小方差法信号源定位模型, 分析了阵列信号相关系数特征和本地次声波实时大气传播特性, 对阵列阵元数量、阵元组成结构引起的定位误差以及本地次声波的真实大气传播射线进行仿真, 并利用中国境内布置的广域次声传感器网络监测到了2013年4月20日四川芦山(雅安)地震的瑞利波激发的本地次声波, 验证了上述模型和仿真, 结合中国地震台网的国内的地震监测台站数据, 从信号走时、信号互相关系数、小波时频图、质点运动轨迹等方面进行了分析与对比, 并使用广域最小方差法搜索的算法对次声波和地震波进行定位, 结果显示: 各次声站点接收到由地震瑞利波引起次声站附近地表震动产生并垂直地表向上传播的次声波, 在地震瑞利波之后到达, 而且相关系数都达到0.6~0.9, 计算得到次声波源方位角为230°(以北京为原点), 距离震中小于150 km, 而且本地次声波受大气传播影响较小, 能够较容易的被广域次声阵列探测到, 因此地震本地次声波监测能够作为地震监测、研究地面起伏运动与大气波动关系的有效手段.Abstract: A kind of least-square-error localization algorithm applied on wide spread infrasound network is proposed in this article. Models of cross correlation between distant sensors and atmosphere infrasound propagation are analyzed. The localization error caused by quantity and distribution structure of network and ray tracing of local infrasound in real atmosphere are also calculated. Infrasound coupled by local seismic Rayleigh wave of Lushan (Ya'an) earthquake on April 20th, 2013 is detected by infrasound network and could prove the algorithm and analysis above. Comparing infrasound signals with seismic recording of IRIS global network, we find that they ware well correlated for the corresponding time period in signal travel time, signal correlation (0.6-0.9), particle motion trajectory analysis, etc.. The zone of infrasound source calculated by the least-square-error localizing algorithm is not compact but its center (minimum value determined by least-square-error method) is less than 150 km distant from the epicenter. Due to the less absorption and refraction in atmosphere propagation, local infrasound is easily detected and recognized and could be a possible and feasible way to monitor earthquake and relationship between ground motion and pressure perturbation in atmosphere.
-
火山、地震、泥石流等自然事件在发生过程中往往会释放出一些特殊的次声波.次声波具有频率低、传播距离远和易于检测的特点.通过观测地震活动伴生的次声波的发生过程，研究次声波的信号特征，有利于科研人员掌握地震活动的状态和研究其发生发展规律及产生次声波的机理.
国内外次声波观测研究中多次监测到地震伴生的次声波.近年来国际禁止核试验条约组织(CTBTO)在全球建立了数十个次声监测站(IMS)，在监测核试验信号的同时，也用于监测地震等自然活动的次声波，中科院声学所安装在香山的次声阵列也监测到2008年5月12日的汶川地震次声波信号(林琳和杨亦春，2010)和2011年3月11日的日本大地震.地震波检测通常使用地震传感器直接检测周边的振动，其中的瑞利波成分引起地表振动向空气中耦合出的次声波，被称作本地次声波(Mikumo and Watada, 2009).本地次声波是由本地的地震波直接耦合产生的，是地震发生后次声台站最早监测到的次声波.由于地震波传播具有速度快、衰减小的特点，因此地震动传感器能够在较远距离上探测到来自震中的快速、高信噪比的地震动信号.
1964年美国阿拉斯加地震产生的本地次声波被美国西海岸的次声站接收到，并与地震瑞利波具有较高的相关性(Donn and Posmentier, 2009; Mikumo and Watada, 2009).2001年中国昆仑山地震(Le Pichon et al., 2003)、秘鲁Arequipa地震(Le Pichon et al., 2002)和2002年阿拉斯加的Denali地震(Olson et al., 2003)都被IMS网络监测到地震引发的本地次声波.2003年日本北海道tokachi-oki地震，日本科学家Watada使用次声传感器监测到了主周期为15~20 s，持续时间20 min的微气压计波动(Watada et al., 2006)，同时，附近的地震传感器监测到的瑞利波与该次声信号具有较高的相关性，认为是瑞利波激发的次声波引起的微气压波动；中国邵长金等(2005)利用CC-1T型传感器接收到了一系列地震次声波，并通过特定时间识别波形幅值等方法分离出本地次声波.但是以上工作只对本地次声波进行了识别和机理研究，并没有联合地震与次声站网络对声源进行定位研究.
本地次声波从次声站附近地表直接耦合产生，受到大气影响较小，较容易观测和提取出来，并进行分析和定位.因此本地次声波能够作为一种地震波通过岩石层向大气层耦合，并影响上层大气和电离层波动的研究手段.中国科学院声学研究所建立了广域次声传感器网络(吕君等，2012)，以监测自然事件产生的大气次声波.在2013年4月20日的四川芦山地震发生后的10 min内，多个次声站都接收到较大幅度且信噪比较高的次声信号，经过与地震动数据的相关性分析，分离提取出有效的本地次声波信号，并得到与震中接近的定位结果.
1. 广域次声波阵列定位理论
广域次声传感器阵列是相对于传统的小尺度次声阵列，其阵列的孔径、阵元数量和覆盖面积有较大的不同：孔径达到1 000~2 000 km，阵元数量达到上百个，在重点监测区域布设为小孔径(1 km)三点阵，使用小尺度次声阵进行时延估计定向只能反映局部传播方向，且地震本地次声波在大气中近似垂直传播，在测量区域内无法确定实际的速度和方位角.因此，在使用广域次声传感器阵列进行地震等自然事件次声波测量时，需要提出适应于该阵列使用的传播模型和定位算法.
1.1 阵列定位模型
由于在地震波传播过程中会出现频散现象，从而影响到各个地点的地震波信号以及本地次声波的波形特征.如果采用传统小尺度次声阵列基于互相关算法估计得到的传感器之间时延量的方法，将会由于信号相关度的显著下降导致时延估计误差的明显增大，从而影响到地震波信号以及地震波耦合产生的本地次声波处理与走时分析精度(崔东源和和跃时，2003；孔选林和李录明，2008).
大气中传播的次声波广域传感器之间同一次声源信号的互相关系数为(Mack and Flinn, 1971)：
$$ \begin{array}{l} \;\;\;\;\;C\left( {\overrightarrow r ,\mathit{T}} \right) = \\ \sqrt {{{\left| {\frac{{\sin (2\pi x\sin \Delta \theta /cT)}}{{2\pi x\sin \Delta \theta /cT}}} \right|}^2} \cdot {{\left| {\frac{{\sin \{ 2\pi y\Delta c/[cT(c + \Delta c)]\} }}{{2\pi y\Delta c/cT(c + \Delta c)}}} \right|}^2}} , \end{array} $$ (1) 式中：T为信号的周期，C为传感器之间信号平均速度，(x, y)是两传感器之间的向量的坐标.Δc为传播过程中沿传播方向波速产生的误差，Δθ为沿波阵面方向出现的来波方位角误差.
对于远距离大气传播的同源次声波，随着传感器之间的信号相关系数随距离增大而减小，据公式(1)，当传感器之间的距离超过300 km时，任何方位角的次声源被传感器接收到信号的相关性都降低0.2以下，这种情况下通过信号互相关计算时延将不够准确，因此需要使用信号峰值检测方法来确定传感器之间的时延.另外由于地震本地次声波频率在0.03~0.20 Hz之间，相应波长λ为2~10 km，由于声源S与阵列的距离满足$d \le \frac{{2{L^2}}}{\lambda }$，则声源S处于近场范围，需要按照球面波来计算，广域传感器阵列中的任意两个传感器i、j之间满足公式(2).
$$ {D_{si}} - {D_{sj}} = \Delta {t_{ij}} \cdot v \gg, \left({i \ne j} \right), $$ (2) 其中:
$$ \begin{array}{l} {D_{si}} = R \cdot \arccos [\sin \left({90 - Lats} \right) \cdot \sin \left({90 - lati} \right) \cdot \\ \cos \left({Lons - Loni} \right) + \cos \left({90 - lats} \right) \cdot \cos \left({90 - Lati} \right)] \cdot \\ {\rm{pi/180}}{\rm{.}} \end{array} $$ 传感器地理坐标为ri=(Lati, Loni)，R为地球平均半径，由于芦山地震位置及其扫描点位于北半球，且由于距离较近(地球大圆距离小于15°)，能够忽略地球表面曲率带来的影响，因此只考虑水平位置之间的地球大圆距离.另一次声传感器与扫描点的距离为Dsj，并设信号传播速度估计值为v≫，其中Δt≫ij为根据信号峰值点估计得到的时延.
1.2 基于最小时延方差的广域次声源搜索定位法
$$ \varepsilon \left( {{r_s}} \right) = \sqrt {{{\sum\limits_{i = 1,j = 1}^{C_n^2} {\left( {\frac{{\left( {{D_{si}} - {D_{sj}}} \right)}}{{}} - \Delta t{ \gg _{ij}}} \right)} }^2}} , $$ (3) 式中：n为广域次声传感器阵列中传感器数量，rs =(Lats, Lons)为区域S上某个地理位置.
当扫描点rs在空间S中扫描，ε(rs) 取得极小值时所对应的rs则是声源最可能的位置，从而实现次声源的定位.对定位精度影响最大的因素是时延估计Δt≫ij的精度.当传播过程中各向异性介质对信号特征带来随机误差，而且随着广域次声传感器网络中阵元距离的增大，通过互相关计算获得时延的精度会降低(Christie and Campus, 2009)，只能通过信号峰值检测的方法计算时延，而信号传播中存在频散效应又使得提取信号峰值时刻的误差增大.这些因素使得信号之间的Δt≫ij也引入了随机误差，并影响到整体的定位精度.
1.3 有效目标信号识别算法
由于次声阵列周围环境背景噪声的存在，当目标平面波到达阵列时，通过互相关计算得到的信号不同频段的各个时延是一系列随机的数值，其中包含有平面波成分和背景噪声所造成的成分，因此需要建立判别方法，以判断接收到的本地次声波的有效性.以三点阵为例，传感器为S1、S2、S3，三角形各边长分别为la、lb、lc，S1、S2位置的夹角分别为α、β.当次声平面波信号波速为c，其与水平面入射角为θ，$\theta \in [ - \frac{\pi }{2}, \frac{\pi }{2}]$，方位角为φ，φ∈[0, 2π]，入射到阵列被接收到，此时各传感器在入射平面波径向的投影点分别为a、b、c.设次声波波速为c，则有公式(4)：
$$ \left\{ \begin{array}{l} la\;\cos \varphi /\cos \theta = c\Delta {t_{12}} = {d_{ab}}\\ lb\;\cos \left( {\varphi - \alpha } \right)/\cos \theta = - c\Delta {t_{31}} = - {d_{ac}}\\ lc\;\cos \left( {180° - \beta - \varphi } \right)/\cos \theta = c\Delta {t_{23}} = {d_{bc}} \end{array} \right.. $$ (4) 声程差之间存在关系：
$$ {d_{ab}} + {d_{bc}} - {d_{ac}} = 0. $$ (5) 理想情况下，三点阵信号之间的时延估计不存在误差的情况，当来自任意方位角和入射角的平面波到达时，传感器之间满足公式(5)的关系，传感器之间时延为Δt12、Δt23、Δt31的闭环和等于零，但实际情况下该闭环时延和只能趋向于零，即满足如下关系：
$$ {r_{123}} = \Delta {t_{12}} + \Delta {t_{23}} + \Delta {t_{31}} < \varepsilon \left({\varepsilon 可以根据阵列孔径得到} \right), $$ (6) 公式(6)可以作为判断信号是否到达的依据.这种方法能够避免次声三点阵附近接收到的随机信号产生的干扰.当次声阵的传感器数量时N≥3，有：
$$ \begin{array}{l} {c_n} = \sqrt {\frac{6}{{n\left({n - 1} \right)\left({n - 2} \right)}}\sum\limits_{_{i > j > k}}^M {r_{ijk}^2}, } i, j, k \in {R_n}, \\ n = 3, 4 \cdots, N, M = C_n^3 + C_n^4 + \cdots + C_n^{n - 1} + 1, \end{array} $$ (7) 式中：Rn为具有传感器数量为n的子阵，当cn或r123小于ε(ε可以根据阵列孔径得到)时，则认为该段次声波信号是有效的平面波信号，由于当地面振动引起的本地次声波以近似垂直的入射角被阵列中各传感器同时接收到，公式(6)和公式(7)能够有效判断本地次声波是否到达，而且有Δt12≈Δt23≈Δt31 < ε.
1.4 广域次声传播模型
广域次声传感器网络的阵元之间的次声波的传播时间和传播距离都大于一个周期和波长，在传播过程中，大气温度梯度和风剖面对次声波声线波阵面矢量的影响都将十分显著.在分析次声波大气传播过程时，使用大气射线追踪理论和大气分层理论模型，由于次声波能量主要在大气的平流层和对流层内传播，因此使用欧洲中期天气预报中心(ECMWF)的ERA-interim的大气分层模型，能够提供高度0~50 km、地理精度0.75°、每6 h更新一次的大气剖面数据.次声声线追踪计算(Drob and Meier, 2009)的公式(8)和公式(9)：
$$ \frac{{\rm{d}}}{{{\rm{d}}t}}\left[ \begin{array}{l} {r_x}\\ {r_y}\\ {r_z} \end{array} \right] = c\left[ \begin{array}{l} {n_x}\\ {n_y}\\ {n_z} \end{array} \right] + \left[ \begin{array}{l} {u_x}\\ {u_y}\\ {u_z} \end{array} \right], $$ (8) $$ \begin{array}{l} \frac{{\rm{d}}}{{{\rm{d}}t}}\left[ \begin{array}{l} {n_x}\\ {n_y}\\ {n_z} \end{array} \right] = n \cdot \nabla (c + \left[ \begin{array}{l} {n_x}\\ {n_y}\\ {n_z} \end{array} \right] \cdot \left[ \begin{array}{l} {u_x}\\ {u_y}\\ {u_z} \end{array} \right]) \cdot \left[ \begin{array}{l} {n_x}\\ {n_y}\\ {n_z} \end{array} \right] - \\ \nabla (\mathit{c} + \left[ \begin{array}{l} {n_x}\\ {n_y}\\ {n_z} \end{array} \right] \cdot \left[ \begin{array}{l} {u_x}\\ {u_y}\\ {u_z} \end{array} \right]). \end{array} $$ (9) 波阵面矢量能够通过公式(9)计算得到.即当地面某一点以一定仰角向大气辐射次声波，次声波声线波会由于大气分层造成的有效声速变化而在分层界面产生折射，并进一步产生弯折而从大气平流层顶向高空逸出或返回地面，在地面产生一系列次声波接收的影区和靶区，并最终影响到次声站的信号接收情况.
1.5 阵元位置误差计算
阵列的定位误差与声阵列的几何分布有关.阵列的几何分布和空间尺寸的差异，将会在某个方向上产生较大的方位角计算误差.在广域次声传感器阵列几何分布确定的情况下，通过估计时延误差，能够得到由于阵列几何分布所带来的方位角计算误差.
对于次声三点阵，其中τ12、τ23、τ31为3组时延估计量，σ02是方差，可以从该阵列接收到的信号大量观测值中统计得到.并且有公式(10)：
$$ \begin{array}{l} \;\;\;\;\;\;\;{\tau _{12}} \sim N\left( {{\tau _{12m}},\sigma _1^2} \right),{\tau _{23}} \sim N\left( {{\tau _{23m}},\sigma _2^2} \right),{\tau _{31}} \sim \\ N\left( {{\tau _{31m}},\sigma _3^2} \right). \end{array} $$ (10) 在广域阵列中，各通道信号时延量的方差较小，因此有σ12≈σ22≈σ32，根据大量阵列信号观察值可以得到σ < 0.05，std(σ12, σ22, σ32)≤0.05
因此将每组时延的方差(这里设置σi=2·σi，i=1, 2, 3)代入公式(3)计算得到确定由于阵列本身分布所造成的方位角计算误差.
2. 仿真分析
2.1 阵元数量对定位性能的影响
广域次声监测网络在使用最小方差法对次声源进行定位时，监测网络所包含的阵元数量会影响次声源定位的精度.因此这里依据公式(5)改变阵元数量的大小，来对阵元数量变化所导致的精度变化进行仿真.假设关于次声监测网络在全国每一个省会城市安装一台传感器，即共有32台传感器，然后将传感器数量减半，即为16台，分别计算两种情况下阵列使用最小方差法搜索次声源位置的精度.根据图 1所示，次声站数量越多，曲线在靠近次声源方位角时下降斜率越大，搜索速度更快且精度更高.
图 1 阵元数量对最小方差法方位角估计精度影响Fig. 1. Accuracy of azimuth estimation affected by number of array spot2.2 阵列排布对定位误差影响
使用基于广域传感器阵列的最小方差法搜索次声源时，除了阵列数量以外，阵列中阵元的相对位置也对声源方位角的计算误差存在较大影响(这里使用四边形阵列进行计算).
以北京、山东、山西、湖北的次声传感器组成的阵列，并且设波速为3 km/s，4 km/s和5 km/s进行仿真，并将获得不同角度上的标准差显示在如图 2所示的雷达图上.结果显示：当来波速度一定的情况下，在对于确定阵元排布的广域次声阵列的声源方位角计算误差随着方位角变化而不同，当声源位于200°~210°或30°时，引入的计算误差最大，当声源位于280°或145°时，引入的计算误差最小.
图 2 次声阵列分布对方位角估计产生的归一化误差Fig. 2. Normalized error of azimuth error caused by distribution of infrasound array2.3 传播特性影响
地震本地次声波由地面垂直振动产生，由于地表起伏，次声波可以认为是向整个空间上半球方向传播.根据公式(8)、公式(9)以及地震发震时刻2013年4月20日8时(UTC+8)的ERA-interim的大气剖面情况，来计算次声站附近地表震动产生的次声波沿地震波传播径向(平行于传播方向)和切向方向(垂直于传播方向)的射线轨迹.如图 3中北京和枝江两地的大气传播射线剖面，以蓝线表示从平流层顶向高空逸出的声线，红线表示能够反射回地面的声线.
图 3 枝江(a, b)与北京(c, d)本地次声波大气传播剖面声线仰角为-90°~+90°，仰角变化间隔为2°，azi为方位角，蓝色为从平流层顶逸出的声线，红色代表能够反射回地面的声线Fig. 3. Atmosphere propagation profile of local infrasound in Zhijiang (a, b) and Beijing (c, d)由于大气有效声速梯度所造成的声线折射作用，枝江次声站附近的本地次声波从-90°~+90°仰角的次声波都向上传播并从平流层顶逸出，北京次声站的径向仰角为0°~10°和切向仰角为0°~3°的本地次声波经过几十千米水平距离的传播返回地面，其他角度的声线都向上传播从平流层顶穿出.由于次声波的能量主要在平流层和对流层传播，在中间层和热层传播时能量衰减十分显著，因此次声站接收到的本地次声波都是由于次声站的地表震动直接产生的，而不是周边地区产生并传播过来的.因此，由次声传感器接收到的本地次声波的速度应与地震瑞利波速度是相同的，且传播方向为垂直地面向上的，所以某地区小尺寸次声阵列的各阵元信号之间是没有时延的.
3. 芦山地震波及本地次声波监测与分析
2013年4月20日8时2分46秒(UTC+8)，四川雅安芦山县境内发生Ms7.0地震，为了研究地震波与本地次声波的关系，本文使用了中国境内的地震和次声台站数据震后相同长度的数据进行波形、时频、相关性等分析.
3.1 地震和次声波监测台站及站点分布
本文中使用IRIS数据库中提供的中国境内的10个台站垂直分量数据，分别为北京(BJT)、西安(XAN)、上海(SSE)、泰安(TIA)、恩施(ENH)、昆明(KMI)、拉萨(LSA)、乌鲁木齐(WMQ)、海拉尔(HIA)、牡丹江(MDJ)和琼中(QIZ)，数据来自于Streckeisen STS-1VBB w/E300，JCZ-1V/VBB地震动传感器.
次声波监测数据来自于中国科学院声学研究所，在北京、山东、湖北等地建设的广域次声传感器网络，北京地区共有6个台站(BJ)，而湖北有襄阳(XY)、枝江(ZJ)两个台站，山东有济南(JN)一个台站，每个台站包含一台INSAS2008型次声传感器和配套的数字采集仪，传感器的3dB频带测量范围为0.008~40.000 Hz，满足对地震本地次声波的测量需求.配套的数字采集仪能够使用GPS信号完成对采集时间基准的校准，同时提供次声监测台站本地的经纬度和海拔高度数据，对次声传感器采集的信号的采样频率为100 Hz，并将AD转换后通过广域网传输到中心服务器上，从而实现对次声波的长期不间断的测量.由于INSAS2008型次声传感器属于电容式次声传感器，其声学系统工作在弹性控制区，即只对周边大气压力波动敏感，但是对传感器本身的振动不敏感(谢金来和陶中达，2003).因此，当地震波传播到次声台站附近地表时，该传感器只能够测量到由于地震动耦合产生的大气次声波，而测量不到由于地震动引起传感器本身的振动.
3.2 芦山地震全国地震台站和次声站接收信号
根据台站与震中距离和地震瑞利波的传播速度，选择地震和次声台站震后1 000 s内的记录数据.如图 4所示的地震波数据，在芦山地震发生(08∶02∶46 UTC+8)后150~700 s内，地震波被国内9个地震台站接收到.通过波形峰值对应的时刻点计算地震波走时和地震波的传播距离，得到结果如表 1所示.
图 4 芦山地震国内地震台站垂直分量监测信息(起始时间2013-04-20T08∶02∶00 UTC+8，发震时间为信号开始后的46 s)Fig. 4. Vertical-component of Lushan earthquake monitored by seismic station表 1 各地震台站地震走时Table Supplementary Table Seismic travel time of each station站点 走时(s) 震中距离(km) 速度(km/s) 地理坐标 KMI(昆明) 184 539 2.93 25.12°N, 102.74°E ENH(恩施) 220 625 2.84 30.28°N, 109.49°E XAN(西安) 234 717 3.06 34.03°N, 108.92°E LSA(拉萨) 384 1 142 2.97 29.70°N, 91.13°E QIZ(琼中) 484 1 398 2.89 19.03°N, 109.84°E TIA(泰安) 505 1 485 2.94 36.21°N, 117.12°E BJT(北京) 543 1 635 3.01 40.02°N, 116.17°E SSE(上海) 614 1 744 2.84 31.09°N, 121.19°E WMQ(乌鲁木齐) 664 2 044 3.07 43.81°N, 87.70°E HIA(海拉尔) 844 2 566 3.04 49.27°N, 119.74°E MDJ(牡丹江) 1 084 2 838 2.62 44.62°N, 129.59°E 如图 5，位于北京、济南、襄阳、枝江的次声传感器分别在震后550 s、505 s、330 s、280 s的时刻接收到了信噪比较大的次声信号，但是相比地震动信号，次声信号中的噪声干扰十分明显.北京是与震中距离最远的次声台站，而且在信号到达时刻附近有较大的噪声干扰，北京次声监测网络(三点阵，孔径600 m)的本地次声波信号互相关分析如图 6a所示，在震后500~750 s的时间段内，次声阵列中的各传感器的信号之间的互相关系数最大达到0.9，属于强相关信号，根据公式(6)的条件，在同一时间段使得三通道信号的各组时延量与闭合时延和同时约等于零，满足本地次声波有效信号到达条件，即本地次声波的传播方向为垂直于地表，与3.3节的仿真结论相同，如图 6b.通过观察信号到达时间，北京的次声信号与北京地震台站信号到达时刻十分接近，根据芦山地震到北京的距离和地震波走时，可以计算得到地震波与本地次声波传播速度都为2.89 km/s，其他次声站山东济南、湖北襄阳、枝江3个次声站的信号，也计算得到相似的传播速度(表 2)，该传播速度远远超过通过大气传播的震中次声波的传播速度，这也符合地震瑞利波传播速度特征.
图 5 芦山地震各次声站本地次声波波形(起始时间2013-04-20T08∶02∶00 UTC+8，发震时间为信号开始后的46 s)Fig. 5. Local infrasound waveform monitored by infrasound stations
图 6 北京次声台站三点阵次声相关系数(a)和闭合时延和与各阵元时延量(b)Fig. 6. Correlation coefficient of Beijing tripartite infrasound array (a) and consistence of time-delay between each element of array (b)表 2 本地次声波信号走时Table Supplementary Table Local infrasound travel time站点 信号传播时间(s) 与震中距离(km) 估计速度(km/s) 北京 565 1 635 2.89 济南 470 1 400 2.97 襄阳 310 910 2.92 枝江 280 821 2.93 3.3 地震波与本地次声波相关性分析
为了确定芦山地震后地震波与本地次声波的相关性，不仅需要验证地震质点位移轨迹是否具备瑞利波的特征(Poggi and Fäh, 2010)，还需要对信号分别进行地震波的小波频谱分析和互相关系数分析.
质点位移轨迹分析能够显示出地震波激发质点振动时所引起的沿着垂直方向和地震波传播方向的运动轨迹.如图 7所示，以恩施地震台接收到的地震波的峰值点的数据为例，质点运动轨迹以一个外径不断减小的近似圆形或椭圆形轨迹进行瑞利波的质点运动，其他各地震台站接收到的地震波的峰值点附近信号具备地震瑞利波的特征.因此可以证明在地震波垂直分量的峰值附近的波主要为瑞利波(Hobiger and Cornou, 2013)，并且瑞利波分布的时间范围与次声波相同.
图 7 芦山地震(恩施ENH)垂直分量信号峰值时刻起始50 s内的质点运动轨迹Fig. 7. Particle motion trajectory analysis (ENH) of 100 s from the time of peak amplitude in Lushan earthquake由于北京同时具有地震台站(BJT)和多个次声台站，能够进行次声传感器之间和次声与地震波之间的相关分析，由于恩施的地震台站(ENH)与位于湖北枝江的次声站相对于震中的方位角都在267°~269°范围内，而且两站的相对距离为200 km，因此恩施与枝江两地的地震动波形是相似的，另外由于济南次声站和泰安的地震台站距离大约60 km，与震中的方位角相差1°以内，所以能够对以上各地的地震与次声信号之间进行比较与分析.
图 10为北京地区和恩施地区的地震波信号和次声波信号的时频小波分析，对比分析了地震信号及其激发出的次声波信号的时频分析图中的信号频率、持续时间和主能量区域变化趋势等.北京地区的地震台站信号(图 8a)与次声台站信号(图 8b)，二者信号开始时间都在600 s前后，即震后大约550 s左右，信号主要能量频率大约为0.05~0.15 Hz之间，在700 s(震后650 s)左右二者信号主要能量的频率慢慢增大到0.10 Hz左右，在800 s左右信号停止.
图 10 最小方差法计算地震源与次声波源Fig. 10. Location of seismic source and infrasound source calculated by least-squared-error method
图 8 次声与地震台站信号小波分析结果a.北京地震台站；b.北京次声台站.起始时间2013-04-20T08∶02∶00 UTC+8，发震时间为信号开始后的46 sFig. 8. Wavelet analysis of infrasound and seismic signal of Beijing seismic infrasound恩施地区的地震台站信号与枝江地区次声台站信号，由于两地都位于相同方位角的地震波传播路径，二者信号开始时间分别为震后大约210 s和290 s左右，走时差80 s，频率大约为0.05~0.15 Hz之间，而且也是随着时间从0.05 Hz增大到0.15 Hz.由于恩施与枝江距离震中距离较近，其信号分别在350 s和420 s左右停止.因此两地的地震波和次声波信号在所研究的时间段内的时频特征是相似的.
图 9是对北京地震波和次声波信号的峰值区域进行的互相关分析，北京地区的地震与次声波信号的相关系数为0.6~0.7，由于湖北的台站距离震中较近，信号信噪比高，相关系数达到了0.7，山东的台站(济南—泰安)信号的互相关系数为0.5~0.6，3组信号都属于显著相关，信号相关系数分布的主瓣值远高于旁瓣值，因此表明地震波信号与次声信号之间的具有十分高的相似性.
图 9 芦山地震后北京地震台站(BJT)与北京次声台站信号互相关分析Fig. 9. Cross-correlation of signals between infrasound and seismic station (BJT) of Beijing地面垂直起伏产生本地次声波的声压可由公式(11)得到：
$$ {p_0} = {\rho _0}{c_0}{w_0}, $$ (11) 式中：p0是次声波压力幅值，ρ0是空气密度，c0是当前空气声速，w0是地面起伏速度.
根据恩施的同震地表位移速度峰值(图 4)约为0.001 m/s，空气密度为1.29 kg/m3，声速为340 m/s，计算可知地表起伏产生的次声波压力幅值0.43 Pa(图 5)，与附近的次声监测站监测到的0.40 Pa相符.
综合以上质点运动的运动轨迹分析、小波分析、信号相关性分析、同震地表位移速度与次声波声压等方面的分析，能够说明北京、枝江、济南地区探测到的芦山地震后最早到达的次声波不是偶然或其他次声源干扰造成，是由于本地地震波中的瑞利波成分引起地面震动耦合产生的.
3.4 地震波源与次声波源定位分析
利用广域阵列定位算法公式(5)分别对地震波源和次声波源进行方位角和地理位置估计.首先根据地震波主能量频率范围(0.03~0.20 Hz)，以北京为坐标原点，以到震中的距离为半径，进行方位角估计.
根据公式(5)，可以得到地震波的方位角估计值，以北京为中心，方位角0°~360°扫描获得在不同波速情况下各地震站点信号之间时延量的方差值，其中最小值指向的即为地震源的方位角.当波速大约在2.5~3.0 km/s之间时，在230°的位置(地震实际方位角为231.5°)得到方差的最小值，因此可以认为该方向为地震源方向、且速度大约为2.5~3.0 km/s，与地震走时估计到的速度相似.将该速度代入公式(5)，设置以北京为原点经纬度正负15°的地理区域，以0.5°的间隔进行搜索，估计得到地震波源的位置如图 10a所示.由于地震波信号接收台站数量多，信号信噪比高，因此得到的定位经过与实际的震中十分接近.
以北京为中心，方位角0°~360°扫描，获得在不同方位角和波速下各次声站信号之间时延量的方差值，其中最小值指向的方位角即为次声源概率最大的方位角.当波速大约在2.5~3.0 km/s之间时，238°的位置得到时延量方差的最小值，因此认为该方向是次声源的方向，且速度也大约为2.5~3.0 km/s.同样，使用该速度进行空间搜索计算，得到次声源估计位置如图 10b所示，计算得到声源区域呈现较大的分布，其中心在32°N、106°E附近，距离实际震中约150 km，且考虑到阵列阵型的空间分布和阵元数量引入次声源位置估计的误差，认为该次声波的声源来自于地震震中附近.
3.5 广域次声阵列定位误差分析
广域次声阵列的几何分布使得某个方向上出现较大的方位角估计误差.在广域次声传感器阵列几何分布确定的情况下，通过估计时延误差，能够得到由于阵列几何分布所带来的方位角计算误差, 方位角误差曲线在200°附近出现峰值，185°左右出现波谷，同样的情况存在于图 10b中.但该最小方差曲线是二次型的，不存在两个极小值，因此可以使用上述方法对该误差进行分析.由于襄阳枝江之间距离远小于阵列孔径，因此将襄阳枝江作为一个阵元，与北京和济南组成大尺度三点阵进行估计.
根据图 2得到归一化阵列方位角误差，图中的3组蓝色误差曲线分别代表不同速度下(2 km/s、3 km/s、4 km/s)方位角0°~360°的方位角计算误差.这种情况是由于4个传感器组成的阵列阵型的空间分布导致在方位角200°附近出现了方位角估计误差的极值，而且随着速度减小，方位角计算误差也会增大.这说明阵列的几何分布可以产生上述定位误差.由于面波传播速度是随其频率和传播途径的地质构造改变的，而大尺度阵列计算中使用整体速度代替局部速度，这也会导致定位出现误差，综上所述，能够说明该地震瑞利波可能来自于芦山地震震中位置.
4. 结论
本文对阐述了广域次声传感器阵列对地震本地次声波信号定位、传播、误差分析的模型，并通过2013年4月20日四川芦山地震后的地震波和本地次声波观测数据进行了分析和信号定位.本地次声波在次声记录数据上是到达最快的次声波，且与地震波垂直分量成分的峰值部分速度相同，而且从大气传播射线仿真和实际信号分析也证明了该次声波是垂直地表向上传播的.
在地震垂直分量最大振幅部分进行了地震动数据的质点运动轨迹分析，显示该部分特征与瑞利波相同.通过对该部分地震波与相应次声波的小波时频分析和相关性分析，二者表现出了较高的相似性，但由于北京市区建筑以及西山地形对次声波的干扰，造成了两地次声信号互相关系数的差异.可以说明，在北京与枝江接收到的次声波是四川芦山地震波中的瑞利波成分传播经过两地时引起地面振动耦合产生的，并非大气中偶然产生的.
地震波源与芦山地震震中位置基本相同，次声波源的方位角与地震波源的方位角指向芦山地震震中方向，且次声波源位置与震中位置相差150 km，该误差能够归因于次声站数量较少以及时延估计误差对计算精度的影响.
相对于地震后另外两种次声波(震中次声波、衍射次声波)，本地次声波更容易被监测到.本地次声波监测作为一种研究地表起伏引起的大气层和电离层波动的手段，也是一种通过监测大气次声波动来研究地震活动的方法.目前，地震前次声波(提前2~15 d)作为在大地震预测的一种可能性手段被研究；同时，由于地震后次声波(本地次声波)与地震具有明确的关联性，其作为验证次声监测网络的可靠性、次声波识别方法以及远距离(几千km)次声源定位算法的有效性是十分必要的；另外，研究地表起伏速度、位移等参数与激发的本地次声波的关系，也有助于认识研究震前次声波的产生机理以及与地面起伏的关联性，进而推动震前次声波的研究.根据芦山地震本地次声波研究中所表现出设备和算法方面的问题，下一步工作要增加广域次声监测网的次声站数量、优化阵型结构、改进定位算法、提高次声波接收及识别能力，为进一步研究震前次声波机理以及震前次声波与地震的关联性提供可靠的依据.
-
图 1 阵元数量对最小方差法方位角估计精度影响
Fig. 1. Accuracy of azimuth estimation affected by number of array spot
图 2 次声阵列分布对方位角估计产生的归一化误差
Fig. 2. Normalized error of azimuth error caused by distribution of infrasound array
图 3 枝江(a, b)与北京(c, d)本地次声波大气传播剖面
声线仰角为-90°~+90°，仰角变化间隔为2°，azi为方位角，蓝色为从平流层顶逸出的声线，红色代表能够反射回地面的声线
Fig. 3. Atmosphere propagation profile of local infrasound in Zhijiang (a, b) and Beijing (c, d)
图 4 芦山地震国内地震台站垂直分量监测信息(起始时间2013-04-20T08∶02∶00 UTC+8，发震时间为信号开始后的46 s)
Fig. 4. Vertical-component of Lushan earthquake monitored by seismic station
图 5 芦山地震各次声站本地次声波波形(起始时间2013-04-20T08∶02∶00 UTC+8，发震时间为信号开始后的46 s)
Fig. 5. Local infrasound waveform monitored by infrasound stations
图 6 北京次声台站三点阵次声相关系数(a)和闭合时延和与各阵元时延量(b)
Fig. 6. Correlation coefficient of Beijing tripartite infrasound array (a) and consistence of time-delay between each element of array (b)
图 7 芦山地震(恩施ENH)垂直分量信号峰值时刻起始50 s内的质点运动轨迹
Fig. 7. Particle motion trajectory analysis (ENH) of 100 s from the time of peak amplitude in Lushan earthquake
图 10 最小方差法计算地震源与次声波源
Fig. 10. Location of seismic source and infrasound source calculated by least-squared-error method
图 8 次声与地震台站信号小波分析结果
a.北京地震台站；b.北京次声台站.起始时间2013-04-20T08∶02∶00 UTC+8，发震时间为信号开始后的46 s
Fig. 8. Wavelet analysis of infrasound and seismic signal of Beijing seismic infrasound
图 9 芦山地震后北京地震台站(BJT)与北京次声台站信号互相关分析
Fig. 9. Cross-correlation of signals between infrasound and seismic station (BJT) of Beijing
表 1 各地震台站地震走时
Table 1. Seismic travel time of each station
站点 走时(s) 震中距离(km) 速度(km/s) 地理坐标 KMI(昆明) 184 539 2.93 25.12°N, 102.74°E ENH(恩施) 220 625 2.84 30.28°N, 109.49°E XAN(西安) 234 717 3.06 34.03°N, 108.92°E LSA(拉萨) 384 1 142 2.97 29.70°N, 91.13°E QIZ(琼中) 484 1 398 2.89 19.03°N, 109.84°E TIA(泰安) 505 1 485 2.94 36.21°N, 117.12°E BJT(北京) 543 1 635 3.01 40.02°N, 116.17°E SSE(上海) 614 1 744 2.84 31.09°N, 121.19°E WMQ(乌鲁木齐) 664 2 044 3.07 43.81°N, 87.70°E HIA(海拉尔) 844 2 566 3.04 49.27°N, 119.74°E MDJ(牡丹江) 1 084 2 838 2.62 44.62°N, 129.59°E 
下载: 导出CSV
表 2 本地次声波信号走时
Table 2. Local infrasound travel time
站点 信号传播时间(s) 与震中距离(km) 估计速度(km/s) 北京 565 1 635 2.89 济南 470 1 400 2.97 襄阳 310 910 2.92 枝江 280 821 2.93
下载: 导出CSV
-
[1] Christie, D.R., Campus, P., 2009. The IMS Infrasound Network: Design and Establishment of Infrasound Stations. In: Le Pichon, A., Blanc, E., Hauchecorne, A., eds., Infrasound Monitoring for Atmospheric Studies. Springer, Berlin, 29-75. [2] Cui, D.Y., He, Y.S., 2003. Attempt on Applying Azimuth of Multi-Station to Calculate the Earthquake Center. Seismological Research of Northeast China, 19(4): 30-34 (in Chinese with English abstract). http://en.cnki.com.cn/Article_en/CJFDTOTAL-DDYJ200304005.htm [3] Donn, W.L., Posmentier, E.S., 1964. Ground-Coupled Air Waves from the Great Alaskan Earthquake. Journal of Geophysical Research, 69(24): 5357-5361. doi: 10.1029/JZ069i024p05357 [4] Drob, D.P., Meier, R.R., et al., 2009. Inversion of Infrasound Signals for Passive Atmospheric Remote Sensing. In: Le Pichon, A., Blanc, E., Hauchecorne, A., eds., Infrasound Monitoring for Atmospheric Studies. Springer, Berlin, 701-731. [5] Hobiger, M., Cornou, C., et al., 2013. Ground Structure Imaging by Inversions of Rayleigh Wave Ellipticity: Sensitivity Analysis and Application to European Strong-Motion Sites. Geophysical Journal International, 192(1): 207-229. doi: 10.1093/gji/ggs005 [6] Kong, X.L., Li, L.M., Luo, S.X., et al., 2008. Seismic Wave Ray Forward in Anisotropy Medium. Computing Techniques for Geophysical and Geochemical Exploration, 30(3): 178-184 (in Chinese with English abstract). http://www.researchgate.net/publication/292920022_Seismic_wave_ray_forward_in_anisotropy_medium [7] Le Pichon, A., Guilbert, J., Vallée, M., et al., 2003. Infrasonic Imaging of the Kunlun Mountains for the Great 2001 China Earthquake. Geophysical Research Letters, 30(15). doi: 10.1029/2003GL017581 [8] Le Pichon, A., Guilbert, J., Vega, A., et al., 2002. Ground-Coupled Air Waves and Diffracted Infrasound from the Arequipa Earthquake of June 23, 2001. Geophysical Research Letters, 29(18): 33-1-33-4. doi: 10.1029/2002GL015052 [9] Lin, L., Yang, Y.C., 2010. Observation & Study of a Kind of Low-Frequency Atmospheric Infrasonic Waves. Acta Acustica. , 35(2): 200-207 (in Chinese with English abstract). [10] Lü, J., Guo, Q., Feng, H.N., et al., 2012. Anomalous Infrasonic Waves before a Small Earthquake in Beijing. Chinese J. Geophys. , 55(10) : 3379-3385 (in Chinese with English abstract). doi: 10.1002/cjg2.1751/full [11] Mack, H., Flinn, E.A., 1971. Analysis of the Spatial Coherence of Short-Period Acoustic-Gravity Waves in the Atmosphere. Geophysical Journal of the Royal Astronomical Society, 26(1-4): 255-269. doi: 10.1111/j.1365-246X.1971.tb03399.x [12] Mikumo, T., Watada, S., 2009. Acoustic-Gravity Waves from Earthquake Sources. In: Le Pichon, A., Blanc, E., Hauchecorne, A., eds., Infrasound Monitoring for Atmospheric Studies. Springer, Berlin, 263-279. [13] Olson, J.V., Wilson, C.R., Hansen, R.A., 2003. Infrasound Associated with the 2002 Denali Fault Earthquake, Alaska. Geophysical Research Letters, 30(23). doi: 10.1029/2003GL018568 [14] Poggi, V., Fäh, D., 2010. Estimating Rayleigh Wave Particle Motion from Three-Component Array Analysis of Ambient Vibrations. Geophysical Journal International, 180(1): 251-267. doi: 10.1111/j.1365-246X.2009.04402.x [15] Shao, C.J., Tang, L., Li, X.F., 2005. Characteristics of Infrasonic Waves Caused by the Ms8.0 Earthquake in 2003 in Hokkaido, Japan. Earthquake, 25(1): 74-80 (in Chinese with English abstract). http://www.researchgate.net/publication/289896019_Characteristics_of_infrasonic_waves_caused_by_the_Ms80_earthquake_in_2003_in_Hokkaido_Japan [16] Watada, S., Kunugi, T., Hirata, K., et al., 2006. Atmospheric Pressure Change Associated with the 2003 Tokachi-Oki Earthquake. Geophysical Research Letters, 33(24). doi: 10.1029/2006GL027967 [17] Xie, J.L., Tao, Z.D., Xie, Z.H., 2003. A High Sensitivity Wide-Band Infrasound Condenser Microphone. Nuclear Electronics & Detection Technology, 23(5): 428-432 (in Chinese with English abstract). http://www.researchgate.net/publication/292668839_High_sensitivity_wide-band_infrasound_condenser_microphone [18] 崔东源, 和跃时, 2003. 利用多台方位角计算震中的尝试. 东北地震研究, 19(4): 30-34. https://www.cnki.com.cn/Article/CJFDTOTAL-DDYJ200304005.htm [19] 孔选林, 李录明, 罗省贤, 2008. 各向异性介质中的地震波射线正演. 物探化探计算技术, 30(3): 178-184. doi: 10.3969/j.issn.1001-1749.2008.03.003 [20] 林琳, 杨亦春, 2010. 大气中一种低频次声波观测研究. 声学学报, 35(2): 200-207. https://www.cnki.com.cn/Article/CJFDTOTAL-XIBA201002017.htm [21] 吕君, 郭泉, 冯浩楠, 等, 2012. 北京地震前的异常次声波. 地球物理学报, 55(10): 3379-3385. doi: 10.6038/j.issn.0001-5733.2012.10.020 [22] 邵长金, 唐炼, 李相方, 2005.2003年日本北海道8.0级地震次声波特征研究. 地震, 25(1): 74-80. https://www.cnki.com.cn/Article/CJFDTOTAL-DIZN200501010.htm [23] 谢金来, 陶中达, 谢照华, 2003. 高灵敏度宽频带电容次声传感器. 核电子学与探测技术, 23(5): 428-432. doi: 10.3969/j.issn.0258-0934.2003.05.011 -
下载:
点击查看大图
计量
- 文章访问数: 3487
- HTML全文浏览量: 143
- PDF下载量: 510
- 被引次数: 0
