引言

北京时间2015年8月12日23点34分左右，位于天津市滨海新区天津港的瑞海公司危险品仓库发生化学爆炸事故，导致165人死亡、8人失踪、798人受伤，并造成了爆炸附近大量建筑、车辆和集装箱受损，直接经济损失达66.86亿元(凤凰资讯，2016)。京津冀地区的数字测震台网和强震台网的部分台站记录到了该爆炸的清晰波形。波形记录显示，这起爆炸事故包含了前后2次较大的明显爆炸，本文记为爆炸-1和爆炸-2。Li等(2015)根据爆炸附近23个数字测震台的波形记录，采用波形互相关技术确定了各台站接收到的2次爆炸的时间差，根据到时差的方位差异，利用相对定位方法确定出爆炸-2位于爆炸-1西北方向55—70m。Zhao等(2016)基于天津市地震局某工程师利用爆炸附近GPS装置得到的爆炸-1的精确位置(117.7496°E，39.0439°N)，进一步分析了爆炸附近32个台站的测震资料，得出了第1次爆炸发生在北京时间23时34分4.68秒，第2次爆炸发生在23时34分36.94秒，间隔约为32.3s；并基于互相关得到2次事件到时差随方位的变化，推断出第2次爆炸发生在第1次爆炸西北方向65m附近；此外基于强震台网的资料，分析了2次爆炸的烈度分布。章鑫等(2016)分析了京津冀地区地电、地磁及地温等多种地球物理场观测对该爆炸事件的同震和震后响应，但是并没有涉及地震学的观测研究。李晔等(2017)通过对爆炸形成的深坑进行实地考察，认为实际爆炸地点为117.7510°E、39.0465°N，与Zhao等(2016)得到的位置相差约300m；并进一步利用118个测震台的资料，分析了爆炸周边Pg、Pn波的波速，其工作主要基于传统震相识别，未能给出2次爆炸发生的时间差及相对位置。

爆炸激发的声波在垂直和水平方向上分别以410m/s和560m/s的速度传播(Zhao等，2016)，且衰减较慢，导致爆炸-2的P波波头被淹没在爆炸-1的声波中，在震中距较大的记录中，无法清晰识别出爆炸-2的P波到时，影响了对爆炸-2的分析。2次爆炸的位置几乎重合，机理相同，因此波形具有一定的相似性，通常采用互相关检测方法(Schaff等，2005)来确定到时差。除了互相关检测，倒谱方法也是一种检测波形相似度和提取信号延时的方法，并且在地震学和声学等领域得到广泛的应用(Cohen，1970；Bakun等，1973；陈运泰等，2000；魏富胜等, 2003, 2010；郭祥云等，2010；高伟等，2013；解滔等，2016)。本文首先假定2次爆炸发生在同一位置，采用倒谱叠加法来确定2次爆炸的时间差；并根据到时差的方位差异，初步估算出爆炸-2相对爆炸-1的方位；同时基于倒谱的前3s特征，探讨爆炸-1前存在更小爆炸的可能性。

1.   数据和方法

1.1   数据

为方便相关人员进行研究，在天津港“8·12”爆炸发生后第2天，中国地震局国家测震台网备份中心(郑秀芬等，2009)立即向社会发布了记录到这2次爆炸100km范围内的32个台站的三分量波形数据。其中，29个台站位于天津市内，3个位于河北省境内(图 1)；除BAD、DOH、SJT、SJZ、ZTZ为宽频带记录台站外，其余27个台站均为短周期记录，数据采样率均为100sps。

图 1  爆炸(红色星号)附近地震台站(黑色三角形)分布

Figure 1.  The distribution of seismic stations (black triangles) near the explosion site (red star)

下载: 全尺寸图片 幻灯片

1.2   方法

在频率域中，角频率为ω的不同幅度的信号沿不同路径传播至接收点，若忽略衰减，则接收点接收到的信号可以表示为：

其中，$ {A_i}$与$ {\tau _i}$分别为第i条路径上信号的幅度和传播时间。功率谱密度为：

其中，$ {\tau _{ij}}$是第i条与第j条路径的传播延时差，“*”表示复共轭。$ \sum\limits_i {{A_i}A_i^*} $为非相干项；$\sum\limits_{i \ne j} {{A_i}A_j^*\cos (\omega {\tau _{ij}})} $为相干项，表征了多路径信号的相互干涉。对公式(2)取对数得到：

对式(3)中的$ {\rm{ln}}\left[ {1 + \frac{{\sum\limits_{i \ne j} {{A_i}A_j^*} \cos(\omega {\tau _{ij}})}}{{\sum\limits_i {{A_i}A_j^*} }}} \right]$做泰勒展开，得到近似表达式：

其中，$ O({\varepsilon ^2})$表示省略掉泰勒展开后的2阶以上小量。对公式(4)做傅里叶反变换，得到接收信号的倒谱：

其中，$ \delta $为狄拉克函数，当自变量为0时，函数值为1，否则为0。由公式(5)可以看出，倒谱$C(\tau) $自变量具有时间的量纲，在$\tau = 0$及$ \tau = {\tau _{ij}}$处存在峰值；在$\tau = 0$处的峰值对应于自相关函数，而在$ \tau = {\tau _{ij}}$处为延时峰值。因此，通过检测倒谱的延时峰值可以估算出路径延迟。作为以上推导的1种特殊情形，即路径相同而信号激励的时刻不同，则延时峰值对应于不同信号激励时刻的时间差。本文假定2次爆炸发生在同一地点，则台站接收到的关于这2次爆炸记录的倒谱包含了2次爆炸发生的时间差信息。以上分析还可以看出，与相关方法相比，倒谱对功率谱取对数是1种同态解卷积处理，有利于滤除激励源频谱起伏的影响(Cohen，1970；高伟等，2013)，且借助于快速傅里叶变换极容易实现。

2.   数据处理及结果

首先按照同一参考时间(北京时间8月12日23时34分0.68秒)截取32个台站的三分量记录，然后去掉线性趋势和均值。对于各台记录的每1个分量，采用2—10Hz的4阶巴特沃斯滤波器进行滤波，并基于快速傅里叶变换，采用周期图法计算功率谱密度函数，最后对功率谱密度函数的对数做傅里叶反变换得到倒谱。图 2(a)给出了DOH台垂直分量记录，由于震中距较大(约93km)，波形信噪比较低。图 2(b)展示了该记录的倒谱前100s的波形，图 2(c)、(d)分别为图 2(b)中0—3s和31—34s的放大波形。由图 2(b)、(c)可以看出，倒谱振幅在0s最大，在0.6s后衰减至0附近，反映了爆炸与天然地震的区别(魏富胜等, 2003, 2010；郭祥云等，2010)；由图 2(d)可以看出，在32—32.5s出现了较为明显的峰值，对应2次爆炸的时间差。考虑到2次爆炸位置并不完全重合，爆炸激发的波在各方向的传播速度有微小差异，同时各记录信噪比不同，峰值对应的时刻也有微小的差异，我们采用叠加来求2次爆炸发生的平均时差。计算中为了减小几何扩散和非弹性衰减的影响，对各记录进行了归一化处理；同时，为了压制0时刻的峰值，对倒谱两端采用了1s的汉宁窗。图 3(a)给出了32个台站垂向记录的倒谱及叠加的结果，图 3(c)为其31—34s的放大，可以看出在32—33s峰值明显增强，噪声得到明显压制，可以判定32—33s间的峰值对应于2次爆炸的时间差，准确值为32.32s。将同样的方法应用到东西分量和南北分量，得到的时间差均为32.31s，忽略掉百分位，可以判定2次爆炸发生的时间间隔为32.3s，这与Li等(2015)、Zhao等(2016)的结果完全一致。Li等(2015)与Zhao等(2016)采用的互相关检测主要利用的是2—10Hz的窗长为4s的回音壁波(whispering gallery wave)，而本文利用的是全波形，采用不同的滤波频段进行处理结果几乎相同，说明相比于Li等(2015)与Zhao等(2016)采用的互相关检测方法，本文采用的倒谱叠加法更为稳健可靠。

图 2  DOH台垂向分量波形记录及其倒谱

Figure 2.  Vertical component waveform and its cepstrum from station DOH

(a)DOH台垂向分量波形记录；(b)DOH台垂向分量波形记录的倒谱；(c)图(b)中0—3s的放大；(d)图(b)中31—34s的放大

下载: 全尺寸图片 幻灯片

图 3(b)给出了0—3s的叠加结果，可以看出在0.5s附近存在1个明显的峰值。李晔等(2017)根据爆炸-1前DAG台与BET台波形的微小起伏，认为在爆炸-1之前存在1个微小的爆炸。如果认为这个微小的爆炸与随后2次爆炸机制相同且发生的位置相近，则理论上也可以通过本文的倒谱叠加方法来获得它与其后2次爆炸的时间间隔。但是，0.5s附近的峰值是否对应爆炸-1与其前面微小爆炸的时间差值得商榷。理由之一是虽然由公式(5)可知倒谱在0s时刻出现最大的脉冲值，但因为序列的长度有限，“脉冲”有一定的宽度且会出现旁瓣，因此无法确定该峰值反映的是爆炸-1与其前面微小爆炸的时间延迟，还是因为信号有限长度造成的假象。加窗操作能在一定程度上抑制这种影响，但同时也会破坏原始振幅。由频率域车厢函数(box car function)的傅里叶反变换可知，旁瓣宽度与频带宽度成反比，衰减与时间成反比(Papoulis，1962)，当旁瓣峰值与延时峰值对应则出现极大值，当旁瓣极小值与延时峰值相对应则可能出现相互抵消的情形，因此很难进行定量分析。图 4分别给出了20—40Hz和10—20Hz滤波后32个台站三分量记录的倒谱单独叠加和总体叠加的结果，在0.5s附近并没有发现明显的峰值。理由之二是假设0.5s附近的峰值对应于爆炸-1前微小爆炸与爆炸-1的时间差，那么该微小爆炸与爆炸-2有约32.8s的时间差，在32.8s附近也应该出现1个峰值，但是图 3显示在32.8s附近并不存在峰值。需要指出的是，以上2个理由并不能否认在爆炸-1前存在微小爆炸的可能，由于微小爆炸能量太小，导致对它的记录信噪比极低，削弱了它与后续2次爆炸记录的相干性；同时，受到倒谱0值处峰值旁瓣的影响，故而难以出判定。

图 3  32条垂向记录的倒谱及其叠加

Figure 3.  The cepstra records of all the 32 vertical component waveforms and the stacking of them

(a)32条垂向波形记录的倒谱(黑色)及其叠加(红色)，其中Δ为震中距(km)；(b)图(a)中0—3s的放大；(c)图(a)中31—34s的放大

下载: 全尺寸图片 幻灯片

图 4  20—40Hz(a)与10—20Hz(b)后0—3s的倒谱叠加波形

Figure 4.  The cepstral-stacking waveforms of 0—3s filtered with 20—40Hz (a) and 10—20Hz (b)

下载: 全尺寸图片 幻灯片

倒谱叠加获得的32.3s处的峰值对应于爆炸-2与爆炸-1的延时，由于2次爆炸的位置可能略有差异以及地壳介质的各向异性，各台站峰值对应的延时也略有差异。对于每1条倒谱记录，在32—33s的窗内搜索倒谱幅值的最大值，其对应的时间为各台站的延迟时间。图 5分别给出了各台站延时与震中距和方位角(相对于爆炸-1)的关系，可见延时与震中距没有明显的关系。由于渤海海域内没有地震台站，致使在方位角90°—190°范围内没有延时观测，但依然可见延时随方位角变化明显。各台站延时与方位角的关系可以用如下函数拟合：

图 5  2次爆炸延迟时间随震中距(a)和方位角(b)的变化

Figure 5.  The observed time intervals between the two explosions varied with epicentral distances (a) and the station azimuths (b)

下载: 全尺寸图片 幻灯片

其中，t为延时(s)；t₀为平均延时(s)；$ \alpha $为方位角(rad)；$ \beta $为初始相角(rad)；$ \Delta t$为延时幅度(s)；$k = 2{\rm{ \mathsf{ π} }}/360 $，为已知常数，具有弧度倒数的量纲(1/rad)。当k已知时，式(6)的拟合可以利用三角函数的展开转化成二元函数线性回归问题，利用线性最小二乘方法进行求解，最终得到公式(6)的具体表达式为：

忽略各向异性的影响，认为这些延时完全是由爆炸-2与爆炸-1的位置不同引起的，则爆炸-2位于爆炸-1延时最小的方向上，对公式(7)求最小值，得到$\alpha $= -7°，即爆炸-2位于爆炸-1北西侧353°处，这与Li等(2015)、Zhao等(2016)估算的352°几乎完全一致。

3.   结论

本文通过对记录到的天津港“8·12”爆炸的32个台站波形数据进行倒谱叠加分析，得到以下结论：①2次明显爆炸发生的时间间隔为32.3s；②依靠倒谱叠加无法确定在爆炸-1前是否存在1个微小的爆炸事件；③爆炸-2发生在爆炸-1西北侧约353°处；④与相关检测相比，倒谱方法在检测识别发生在同一地点的多次爆炸及类似事件方面具有一定的优势。