随着小波变换研究的日渐深入，小波分析作为一种有效的时频分析方法在地震学研究和应用中得到了越来越广泛的应用(杜兴信，1997；于锦海等，1999；刘希强，2000，2002；徐义贤，王家映，2000；严尊国等，2000；杨文采等，2001；李杰等，2005)，已经成为一种不可或缺的数学处理工具. 实际工作中，利用Matlab时频分析工具箱可以很好地实现信号的小波时频分析，而对于通过二进小波变换和小波包变换实现信号的带通滤波，则需要明确小波包分解树结点(node)与信号子空间频带的对应关系. 但通过小波包变换得到的信号子空间频带内的频率并非按照分解树结点编号的大小顺序排列，各个结点重构信号的频率范围不易判定，不便于实现信号的带通滤波.

对于小波包分解树结点与信号子空间频带对应关系的研究，目前尚未见到相关文献. 本文通过分析小波包变换的Mallat分解算法与分解滤波器的关系，揭示分解树结点与信号子空间频带的对应关系，从而确定分解树每个结点所对应的重构信号的频带范围. 随着测震台网和地震前兆观测网络“十五”建设项目的完成和“十一五”建设项目的开展，将产生海量的测震和前兆数字化观测数据，本文可为通过小波包变换实现数字信号指定带通滤波提供依据，对小波包理论实际应用于数字化观测具有一定意义.

小波包{W_n(x)}_n∈Z+定义如下(李弼程, 罗建书，2003):

信号f(x)在子空间Ω_n^j中的小波包变换系数(李弼程，罗建书, 2003)为

小波包变换的Mallat分解算法(李弼程, 罗建书，2003)为式中，C_k^{2n, j-1}和C_k^{2n+1, j-1}分别表示信号在子空间Ω_2n^j-1和Ω_2n+1^j-1中的小波包变换系数.

由式(3)可知，小波包变换的Mallat分解算法可通过滤波器组电路来实现. 首先通过低通分解滤波器(LP)和高通分解滤波器(HP)将信号f(x)的频带0—ω平均分成低频0—ω/2和高频ω/2—ω两个频带，然后分别对低频变换系数和高频变换系数1/2抽样后，再分别通过滤波器组LP和HP，这样频带0—ω/2即可平均划分成频带0—ω/4和ω/4—ω/2，频带ω/2—ω可平均划分为ω/2—3ω/4和3ω/4—ω两个频带. 如此重复下去，最后得到相等带宽且频带互不重叠的分解信号. Mallat算法的塔式结构如图 1所示.

图 1(Fig. 1)

图 1 小波包变换的Mallat分解算法LP和HP分别表示低通分解滤波器和高通分解滤波器 Fig. 1 Mallat decomposition algorithm of wavelet packet transform LP and HP respectively represent low-pass and high-pass decomposition filter

对被分析信号进行某尺度的小波包分解，如果某信号子空间频带是经过偶数次HP而得到的，则对该频带进行再分解时，得到的两个均匀带宽的频带会按照低频带到高频带的顺序排列；反之如果子空间频带是经过奇数次HP而得到的，则对该频带进行再分解后所得到的两个均匀带宽的频带会按照高频带到低频带的顺序排列. 信号子空间频带排列顺序与LP无关.

图 2中，假定某信号的采样率为64 Hz，根据Nyquist采样定理，可检测信号频带为0—32 Hz. 对该信号进行尺度j=4的小波包分解，信号频带则被均匀分成2⁴个频带，每个频带的宽度为2 Hz. 按照上述频带排列规则，结点[1, 1]对应的频带16—32 Hz是经过1次HP得到的，对该频带再分解得到的两个频带按照24—32 Hz和16—24 Hz排列，它们分别对应结点[2, 2]和[2, 3]；结点[3, 5]对应的频带24—28 Hz是依次经过HP，LP和HP得到的，共2次HP，对该频带再分解得到的两个频带按照24—26 Hz和26—28 Hz排列，它们分别对应结点[4, 10]和[4, 11]. 依据以上频带排列规则，其余结点对应的频带见图 2.

图 2(Fig. 2)

图 2 小波包分解树结点与信号子空间频带的对应关系 Fig. 2 Corresponding relationships between nodes and frequency-bands

对某信号进行尺度为j的小波包分解，则该信号频带0—ω将被均匀划分成2^j个互不重叠的频带，带宽为ω/2^j，各个频带的范围为nω/2^j—(n+1)ω/2^j(j=1, 2, …； n=0, 1, 2, …, 2^j). 由于n与频带下限(或上限)一一对应，因此用n表示频带编号. 图 2中方括号内的数字分别表示小波包分解树的分解尺度j和该分解尺度下的结点编号m.

经分析发现，当m值和n值均用二进制数表示，通过设定二者之间的运算规则，可以实现n到m的转化. 如果n值的二进制表示有k个连续的1，对k个连续的1中最高位1后面的k个数均取非，得到的数则为m值的二进制表示. 例如，n=11，换算为二进制数为1011，按照以上的运算规则转化为二进制数1110，换算为十进制数14，即m=14. 如果ω=32 Hz，j=4，则频带22—24 Hz对应结点[4, 14]，与图 2中所示的对应关系一致. 假定该信号的采样率为64 Hz，可检测信号频带为0—32 Hz. 通过该运算规则，不同分解尺度下信号子空间频带与结点的对应关系见表 1. 与图 2相比，二者是一致的.

表1(Table 1)

表1 结点与信号子空间频带的对应关系Table 1 Corresponding relationships between nodes and bands

表1 结点与信号子空间频带的对应关系 Table 1 Corresponding relationships between nodes and bands

设模拟信号由频率分别为1, 5, 13 Hz和27 Hz的正弦信号叠加而成，解析表达式为S=sin(2πt)+sin(10πt)+sin(26πt)+sin(54πt)， t∈[0, 10]，取信号采样率为64 Hz. 采用DMeyer小波基对该模拟信号进行4层小波包分解，其分解树见图 2. 分解树最底层结点变换系数灰度图见图 3. 图中横轴代表时间，纵轴代表分解树最底层互不重叠的 分解频带，自下到上按照分解频带由低到高排列，纵轴数字代表与分解频带对应的结点编号，灰 度深浅代表系数绝对值的大小.由图 3可以看出，模拟 信号S主要包含4种频率成分，且各成分的信号均为显著 的周期信号. 同时，图 3中纵轴的结点编号和分解频带的对应关系与图 2中小波包分解树最底层结点和信号子空间频带的对应关系是一致的. 这说明以上总结的小波包分解信号子空间频带的排列规则是正确的. 通过对分解树最底层结点变换系数进行重构，得到了各个结点对应的不同分解频带的重构信号(图 4). 由图 4可以看出，模拟信号S由4种频率成分的信号叠加而成，分别为1, 5, 13 Hz和27 Hz的正弦信号，它们对应的结点编号分别为[4, 0]，[4, 3]，[4, 5]和[4, 11]. 这种对应关系与图 2中小波包分解树最底层结点和信号子空间频带的对应关系是一致的，再次说明了以上总结的小波包分解信号子空间频带的排列规则是正确的.

图 3(Fig. 3)

图 3 分解树最底层结点变换系数灰度图 Fig. 3 Gray image of transform coefficients of nodes at the most bottom of decomposition tree

图 4(Fig. 4)

图 4 分解树最底层结点变换系数重构信号 Fig. 4 Reconstruction signals from transform coefficients of nodes at the most bottom of decomposition tree

运用频带能量比方法(杨选辉等，2005)，对宁夏及邻区震中距在55 km左右的地震与爆破进行识别，这样的震中距足以满足所取数据均在P波段内. 本文事件参数取自宁夏地震观测报告，事件记录选取信噪比较高的台站记录. 记录采样率为50 Hz，选取P波段垂直向记录的256个数据点. 如表 2所示，序号1-5的事件为地震事件，序号6-10的事件为经过落实后确认的爆破事件.

表2(Table 2)

表2 宁夏及邻区的地震事件和爆破事件目录Table 2 Catalogue of earthquakes and explosions locating in Ningxia and its adjacent area 序号 日期 发震时刻 φN/° λE/° 地点 深度/km ML 台站 震中距/km 年-月-日 时:分 1 2003-10-26 21:20 38.15 106.07 青铜峡 26 1.7 YCH 53 2 2004-01-24 09:02 39.10 106.00 阿拉善左旗 15 2.0 YCH 54 3 2004-07-23 06:29 38.15 106.20 永宁 28 1.8 YCH 56 4 2004-07-24 14:52 38.15 106.23 灵武 16 1.7 YCH 57 5 2005-01-13 03:45 38.15 106.25 灵武 14 2.2 YCH 59 6 2005-11-16 11:31 39.13 106.08 阿拉善左旗 0 2.7 YCH 60 7 2005-11-16 14:22 39.10 106.17 阿拉善左旗 0 2.5 YCH 58 8 2005-11-24 12:50 39.12 106.08 阿拉善左旗 0 1.8 YCH 54 9 2005-11-24 16:25 39.10 106.08 阿拉善左旗 0 2.3 YCH 56 10 2005-11-25 13:20 37.13 105.37 同心 0 2.1 ZHW 58 注：YCH代表银川台，ZHW代表中卫台.

表2 宁夏及邻区的地震事件和爆破事件目录 Table 2 Catalogue of earthquakes and explosions locating in Ningxia and its adjacent area

在此采用11阶Daubechies小波基对数据段进行尺度j=2的小波包分解. E₀，E₁，E₂和E₃分别表示频带0—6.25 Hz, 6.25—12.5 Hz, 12.5—18.75 Hz和18.75—25 Hz内的能量. 根据某频带内的小波包变换系数的平方和等于该频带内的信号能量(曾宪伟等，2008)

图 5给出了表 2中所列事件不同频带能量比的自然对数值. 爆破信号的ln(E₀/E_i)均大于地震信号的ln(E₀/E_i)(i=1, 2, 3)，且爆破信号频率越高能量越低，该能量随频率衰减的速率远大于地震信号衰减的速率. 由此说明，地震信号和爆破信号在不同频带上的能量比是能够分开的，且爆破记录频率越高信号越弱. 这是因为爆破事件发生于地表，信号在传播过程中有较长的路程在浅层，其频率越高衰减越快.

图 5(Fig. 5)

图 5 不同频带能量比值图 Fig. 5 Energy ratio of different frequency bands

通过分析小波包变换的Mallat分解算法与分解滤波器的关系，并由滤波器实现电路联想到设定频带编号与结点编号间进行二进制转化的运算规则，得到了小波包分解树结点与信号子空间频带的对应关系，从而揭示出小波包信号子空间频带的排列规则. 通过模拟信号给出检验，结果表明本文总结的小波包信号子空间频带的排列规则是正确的. 利用这一规则，可以确定小波包分解树每个结点对应的重构信号的频带范围. 通过将频带间的能量比转化为相应结点变换系数的平方和之比，对宁夏及邻区的地震与爆破进行了有效识别. 由此，本研究可为通过小波包变换实现数字信号指定带通滤波提供依据，也可为小波包理论应用于测震和前兆数字化观测提供支持.