第27卷 第4期2006年 8月地震地磁观测与研究SEISMOLO GICAL AND GEOMA GN ETIC OBSERVA TION AND RESEARCH Vol 27 No 4Aug 2006
时频分析方法及在地震波谱研究中的应用张 帆1),2) 钟羽云2) 朱新运2) 熊 丹2) 周 昕2)
1)中国合肥230026中国科学技术大学地球与空间科学学院
2)中国杭州310013浙江省地震局
摘要 地震波是一种非平稳信号，傅立叶分析已不适合。时频分析能对非平稳信号进行分析，并给出地震波信号能量的时间2频率分布。但是在时频分析时，存在交叉干扰项，会带来一些错误的信息。本文讨论了波恩2约旦(Born 2Jordan )时频分布、乔伊2威廉斯(Choi 2Williams )时频分布、赵2阿特拉斯2马克斯(Zhao 2Atlas 2Marks )时频分布3者之间抑制和消除交叉干扰项的特点，认为赵2阿特拉斯2马克斯(Zhao 2Atlas 2Marks )时频分布抑制交叉干扰项的效果和时频聚集性更好，更适合于分析地震波。
关键词 非平稳信号；时频分布；数字化地震波；交叉干扰项；时变频谱
中图分类号:P315.3+1 文献标识码:A 文章编号:100323246(2006)0420017206
作者简介：张帆(1978～)，男，工程师，主要从事地震预报工作
本文收到日期:2006202227
引言
信号的各种统计量(如频率、相关函数、功率谱等)随时间变化的信号称为非平稳信号。经典的傅立叶分析已不适合于分析非平稳信号，它只能把信号分解成单个的频率分量，并建立起每个分量的相对强度，而不能得到什么时间存在哪些频率。因此人们建立一种分布，在时间和频率上同时表示信号的能量或者强度，即时频分布，其目的就是用这种分布来研究信号的频率或者频谱含量随时间的变化以及信号的时变频谱。自从Ville 推导出Wigner 在1932年研究量子统计力学时所得到的一种分布———Wigner 时频分布后，这种分布在时频分析的研究过程中起到了非常重要的作用，并由此演化出多种时频分布。
时频分布可分为线性和非线性时频表示。常用的线性时频表示有短时傅立叶变换(STF T )、小波变换(W T )和戈勃(Gabor )展开等。而非线性时频表示则是二次型的能量分布表示(胡昌华等,2002)，常用的非线性时频表示有WV 时频分布、马根诺2希尔时频分布、乔伊2威廉斯时频分布等。本文所分析的3种时频分布属于后者。
时频分布已经在很多领域(包括工程、物理、天文学等)得到了广泛的应用，取得了很有价值的成果。随着对地震波认识的深入，地震波被看作是一种非平稳信号，而现代化的数字化地震仪提供了采样率
和分辨率较高、记录频带宽、动态范围大的数字化地震波形资料(何少林等,1997)，利用计算机可以方便地对地震波进行时频分析。刘希强等(1998,2000,2003)已经用线性的小波分析方法来识别地震弱震相，区别爆破、塌方和地震。利用非线性时频分布来研究地
71
地 震 地 磁 观 测 与 研 究27卷震波还比较少见。本文讨论了波恩2约旦时频分布、乔伊2威廉斯时频分布、赵2阿特拉斯2马克斯时频分布3者之间抑制和消除交叉干扰项的特点，认为赵2阿特拉斯2马克斯时频分布抑制交叉干扰项的效果和时频聚集性更好，更适合于分析地震波，提高地震波时频分布的可读性。1 理论和方法
王者荣耀ar相机解析信号s (t )的Wigner 分布为
起诉欠钱不还的程序和费用W (t ,ω)=12π∫s 3(t -12τ)s (t +12τ)e -j τωd τ(1)
式中3表示解析信号复共轭。对于两个信号的和s (t )=s 1(t )+s 2(t )带入(1)式得到
W (t ,ω
)=W 11(t ,ω)+W 22(t ,ω)+W 12(t ,ω)+W 21(t ,ω)(2)式中W 11(t ,ω),W 22(t ,ω),W 12(t ,ω),W 21(t ,ω)分别是s 1(t ),s 2(t )的自项和他们之间的交叉干扰项。
所有的时频表示都可以由下面的方程得到
C (t ,ω)=1
4π2∫∫∫s 3(u -12τ)s (u +12τ)<(θ,τ)e -j θτ-j τωτ+j θu d τd u d θ(3)
<(θ,τ)是二维函数，叫做核(Claasen et al ,1998)。取不同的核函数可得到不同的时频分布。对于一个简单的多分量信号s (t )=s 1(t )+s 2(t )=A 1e j ω1t +A 2e j ω2t 分布为
C (t ,ω
)=C 11(t ,ω)+C 22(t ,ω)+C 12(t ,ω)+C 21(t ,ω)(4)
其中
C 11=|A 1|22π∫<(0,τ)e -j τ(ω-ω1) C 12=A 31A 2e -j (ω1-ω2)t K (ω
)而
K (ω)=∫<(ω1-ω2,τ)e -j τ(ω-ω12)d τ(5)其中ω12=(ω1+ω2)/2
假定ω2>ω1，那么K (ω
)的性质决定了交叉干扰项的大小和分布(Cohen ,1995)。选择不同的核得到不同的K (ω
破釜沉舟百二秦关终属楚)，不同大小和分布的交叉干扰项。对于波恩2约旦时频分布(Cohen ,1966)，取核
<BJ (θ,τ)=sin (θτ/2)θτ/2
(6)K (ω)=1/(ω2-ω1) (ω1≤ω≤ω2)(7)
此时，交叉干扰项完全地控制在两个频率之间。对于乔伊2威廉斯时频分布，取核<CW (θ,τ)=e -θ2τ2/σ(8)
C 12=2A 1A 2co s [(ω2-ω1)t ]η(ω,ω1,ω2,
σ)(9)其中η(ω,ω1,ω2,σ)=1
4π(ω2-ω1)2/σexp [-(ω-ω)2
4(ω2-ω1)2/σ](9)式说明，只要σ有限，交叉干扰项将有限，并且不集中在一点，而随着较低的最大强度延伸(白居宪译,1997)。对于赵2阿特拉斯2马克斯时频分布，取核
<ZAM (θ,τ)=g (τ)|τ|sin (αθτ)αθτ(10)
81
第4期张帆等：时频分析方法及在地震波谱研究中的应用在赵2阿特拉斯2马克斯的原著中,g (τ
)取1,α取1/2。此时K (ω)=14πj ∫2|τ|(ω2-ω1)[e -j τ(ω-ω1)-e -j τ(ω-ω2)]d τ(11)(11)式说明K (ω)是(ω-ω1)和(ω-ω2)的函数，交叉干扰项位于自项中间。2 理论信号的3种时频分布及Wigner 时频分布的比较
设计3个高斯元组成的理论信号，采样率1Hz ，折叠频率为0.5Hz 。取第1个高斯核的图1 3个高斯元信号及其时频图解
频率中心为0.3Hz ，时间中心为32s ，持续
时间为32s ，振幅为1；第2个高斯核的频
率中心为0.15Hz ，时间中心为56s ，持续
时间为48s ，振幅为1.22；第3个高斯核的
频率中心为0.41Hz ，时间中心为102s ，持
续时间为20s ，振幅为0.7。此信号及其时
陈美诗频图解见图1。图1中从左到右分别是高
斯元1、元2、元3。元1、元2和元3两两之小海绵的照片>新浪 评论
间在频率域无重叠；元2和元3、元1和元3
在时间域无重叠；元1和元2时间域有
重叠。
分别把(6)式、(8)式、(10)式代入(3)式，求出该信号的BJ 、CW
、ZAM 分布，并将它们同Wigner 时频分布对比(图2)。
图2 图1中信号的WV ,CW ,BJ ,ZAM 时频分布的比较
9
1
地 震 地 磁 观 测 与 研 究27卷
Wigner分布中3个高斯元中两两之间存在交叉干扰项，并且元之间的距离越近，幅度越大，交叉干扰项的幅度也就越大。乔伊2威廉斯时频分布中，虽然元1和元3及元2和元3之间的交叉干扰项已经被消除，但时间域有重叠的元1和元2之间，在整个重叠区域存在与频率轴平行的小幅交叉干扰项；波恩2约旦时频分布中，元1和元3及元2和元3之间的交叉干扰项已被消除，但元1和元2之间的交叉干扰项相比,CW分布已经有效地控制在元1和元2的时间重叠区域之间；赵2阿特拉斯2马克斯时频分布中，高斯元两两之间的交叉干扰项已被消除。相比之下,4种分布中时频聚集性较好的是Wigner分布和赵2阿特拉斯2马克斯分布。
根据理论信号的分析，说明3种时频分布中,ZAM时频分布消除交叉干扰项和时频聚集的特性更显著。因此，用ZAM时频分布表示地震波能更好地揭示地震波的频率和时间的变化关系。
3 资料处理和结果
选取2002年7月28日19时39分文成2泰顺交界M L3.5地震以及2002年9月25日
06
图3 永康台M L3.5地震记录垂直分量时43分台湾省台东县M S4.9地震永康地震台垂直分量进行处理并比较。永康台所用的数字地震仪是FSS23型短周期地震仪。该台离震中的距离分别约为140km 和694km。计算之前扣除了仪器响应(胡昌华等,2002)，根据累计梯形积分算法，利用Matlab中提供的cumt rapz函数(李海涛等,2002)，把速度记录积分转化为位移记录，并进行滤波等处理。M L3.5地震记录仪器校正后的结果见图3，从上到下分别是原始记录，速度记录，位移记录。
实信号的能量密度谱关于原点对称，一半对应于负频率，一半对应于正频率，此时，平均频率为零，频率扩展是正负频率谱尖峰之间的距离。物理上并不存在负频率，只需取正频率信息，应把实信号转化成解析信号，即
^s(t)=s(t)+j・H T[s(t)]
式中H T[s(t)]表示信号的Hilbert变换。
使用解析信号可以消除信号的傅立叶变换的负频值，使其负频率的那一半频谱值为零，从而得到正频率频谱的平均频率和频率扩展，同时使用解析信号还可以消除时频分布中的负频率区的成分和正负频率区之间的交叉干扰项。
将赵2阿特拉斯2马克斯时频分布核代入(4)式，并将(4)式中的信号s(t)换成^s便可计算出地震波的C ZAM(t,ω)。两次地震P波前12s垂直分量记录及ZAM分布见图4，给出了该记录95%的能量分布。图4中可以清楚的分辨出什么时刻存在着什么频率。文成地震的能量主要集中在1～2Hz之间，存在7Hz左右的频率成分，并且在2s时P波的振幅最大，此时的P波的主频率为1.5Hz左右。随后，地震波主要是1～2Hz的频率成分，随着时间的延续，地震的能量有所衰减，到10s左右能量基本衰减到最大能量的5%以下。平均频率约为2.9Hz，频率扩展约为7.9Hz。台东地震的能量集中在0.5～2Hz之间,2Hz以上频率成分能量很弱, 02
第4期张帆等：时频分析方法及在地震波谱研究中的应用振幅在约9s 处达到最大，平均频率约为0.98Hz ，频率扩展约为1.5Hz
。
图4 永康台记录P 波垂直分量及ZAM 时频分布
(a )文成地震;(b )台湾省台东地震
4 结论
比较3个高斯核构成信号的3种时频分布特点，得出以下结论。
(1)对于CW 时频分布，当信号在时域和频域有重叠时交叉干扰项比较突出，范围比较大，尤其是在时域重叠时。
(2)对于BJ 时频分布，信号的交叉干扰项都被限制在两个信号之间。
(3)对于ZAM 时频分布，信号的交叉干扰项被很好的抑制。
(4)ZAM 的时频聚集性优于其他两种分布，更接近于WV 时频分布。
从文成地震和台东地震的永康台P 波记录的时频分布看出了两者的特征区别，较近的地震其平均频率高于较远的地震，频率扩展也高于较远的地震。较近的文成地震的记录高频成分较多，主要能量的集
中频率比台东地震的高，这是由于台东地震的震中距大于文成地震，传播路径较长，介质过滤了台东地震高频成分的缘故。在记录的前5s ，文成地震的成分主要为1～3Hz ，台东地震主要为0.5Hz 。这是由于永康台记录到了台东地震的Pn 震相，文成地震的Pn 震相不发育，而理论上Pn 震相的频率要比Pg 小。
地震波是一种非平稳信号，它的频率等统计量随着时间而变化，经典的傅立叶分析不能反映出地震波的这种非平稳性。时频分析可以在时间和频率平面表示地震波的能量分布，解释地震波的非平稳性，并可研究地震波的时频谱。因此，时频分析在地震的模式识别等方面有很大的应用潜力。由于时频分布是信号的双线性函数，不可避免存在着交叉干扰项。在时频分析时，应尽可能抑制和消除这种干扰项。综上所述，认为ZAM 时频分布消除交叉干扰项和时频聚集能力较好，更能精确地表示地震波的时频特性。
12
发布评论