地震动分析是研究地震危害控制的基础，也是控制地震危害的前提。小波包变换可将地震信号均匀分解到不同的频带上，从而准确分析地震信号各频带的特征。但在MATLAB平台中，由于算法的限制，小波包变换后得到的树结点顺序与频带顺序并不一致，即与树结点对应的频带并不是依次由低到高，而会出现“跳频”现象^[1-4]。若忽视此问题，则不能准确定位频带，进而影响后续时频分析的正确性。

现有方法^[3]虽然已经明确了上述的“跳频”现象并得到了树结点与频带结点的排列规律，但其只能按照所得规律手动进行序号对照，效率太低。本文借助小波包算法和异或运算定义，以唐山南北向地震波为例，基于MATLAB平台，编制相应程序，得到的地震频带序号与其小波包树结点序号的排列规律与现有方法^[3]得到的完全一致。对地震信号进行小波包高层次分解时，利用该转换算法能够快速高效地得到各个频带序号对应的结点序号，对快速高效地提取地震信号局部特征、有层次地展现地震信号等具有重要意义。

1 小波包理论

设小波包变换中的正交尺度函数为ϕ(t)，小波函数为ψ(t)，由它们产生的滤波器组中低通滤波器为l(n)，高通滤波器为h(n)，根据多分辨分析理论，ϕ(t)、ψ(t)、l(n)和h(n)之间满足双尺度方程^[5]：

(1)

小波包分解过程可以用分解树来描述，如图 1所示。图中U_{j, n}表示小波包分解的子空间编号，也可以称为结点编号，j表示分解层数，其对应的分解尺度为2^j，第j层一共有2^j个正交基，n表示第j分解层的子空间编号，第j层的第n个子空间由集合{U_{j, n}(t-k)}构成。

由小波包变换理论知识可知，同一层的各个子空间之间相互正交，且每一层上的所有子空间的直和构成这个信号空间。

2 地震信号频带划分参数

对地震动非平稳信号进行小波包变换时，主要参数有分解层数和分解路径。

2.1 分解层数

进行小波包分解时，除了确定小波基函数，还需要确定分解层数。选用字母j表示小波包分解层数，通常分解层数按式(2)选取^{[2, 6]}：

(2)

式中，L_s为输入信号的长度。对于地震非平稳信号，采样时间间隔一般为0.01 s或0.02 s，即采样频率为100 Hz或50 Hz。地震动信号一般持续时间为10~20 s，若采样频率为50 Hz，则信号长度为2⁹~2¹⁰。综合考虑频带宽度及分辨率，本文认为分解层数取7能满足地震信号的分析要求。

2.2 分解路径

如上所述的小波包分解树中，树结点从左至右依次编号为0, 1, …2^j-1，称作结点序号。地震信号的小波包分解实际上是把信号通过一组滤波器组合(低通或高通)后，到达底端的树结点。设定0表示低通滤波，1表示高通滤波，则各滤波器组合可由一个二进制数表示^[7]，称该二进制数为实际路径。例如，实际路径101表示信号经过一次高通滤波器，一次低通滤波器，最后又经过一次高通滤波器，如图 2中虚线所示。该二进制数对应的十进制数就是分解终端树结点的频带序号，因此实际路径又称频带路径。与频带路径相对应，把结点序号的二进制数称作结点路径或理想路径。

除了频带路径和结点路径，每个结点还对应一个通频带。通频带是指通过该结点的频带范围，与频带序号有关。若信号的有效频带为f，分解层数为j，则频带序号m对应的通频带为m·(f/2^j)-(m+1)·(f/2^j)，其中，m=0, 1, …, 2^j-1。

需要说明的是，分解路径的每个二进制数的位数为获得该频带所需要的分解层数，这样可以准确获取各个频带分解路径，尤其是以0开头的分解路径，同时为后续准确得到频带序号及对应的通频带作铺垫。

3 频带序号到结点序号的快速转换算法

当分解层数较少时，可以根据现有规律^[3]手动依次得到每个结点的频带路径和通频带，进而得到完整的小波包树。但当分解层数较多时，手动获取比较困难，需要找到频带序号和结点序号之间的快速换算关系。通过对多个结点的分析，得到频带序号到结点序号的快速换算过程如下：对频带序号进行二进制转换得到频带路径，即一个二进制数；从左至右依次对该二进制数的每位数与其左侧位数进行异或运算，得到结点路径，即一个新的二进制数；对该新的二进制数进行十进制转换，即可得到结点序号。比如，经过3层小波包分解得到的小波包树，频带序号为a，a的频带路径(二进制数)为ABC，对其进行异或运算后可得到结点路径为XYZ(二进制数)。其中，X=A、Y=A⊕B、Z=B⊕C，将二进制数XYZ转化为十进制数b，则b为频带a对应的结点序号。

上述“异或”是一种数学逻辑运算符，数学符号为“⊕”，计算机符号为“xor”，运算法则为：

(3)

根据上述换算过程和运算法则，在MATLAB平台上编制程序，可以实现小波包分解到任意分解层数时，地震信号子频带序号与结点序号的转换，表 1为7层小波包变换下的频带序号与结点序号的对应关系。

4 实例验证

选取唐山地震波南北向(北京观测)的强震记录进行实例分析。地震信号采样点数为2 000，采样的时间间隔为0.01 s，采用sym5小波函数，在MATLAB平台上编制程序，对唐山波进行7层小波包变换，得到各树结点的频谱图。本例中，采样频率为100 Hz，根据采样定理，有效频率为0~50 Hz，进行7层分解后每个频带宽度为f/2⁷，即50/2⁷ Hz。限于篇幅，这里给出频带序号分别为8、13、37、38、56、60的频谱图，如图 3所示。按照本文给出的变换算法，参照表 1的频带序号与结点序号的对应关系，取图 3中的各点如表 2所示。

5 结语

小波包变换可将地震信号均匀分解到不同的频带上，但由于算法限制，小波包变换后得到的树结点顺序与频带顺序并不一致，且现有方法不能快速高效地得到频带与小波包树结点的对应关系。本文针对这一问题，基于MATLAB平台，从小波包算法和异或运算的角度，剖析了滤波器组实现小波包变换算法，分析了频带均匀划分理论。借助分解路径，总结了小波包变换的频带顺序排列规律，并由此得到了小波包分解树的频带序号和树结点序号的转换算法以及各频带对应的通频带，为确定频带序号到小波包分解树结点序号的转换关系提供了一种有效的方法。

根据总结的小波包变换的频带顺序排列规律和频带序号到结点序号的转换算法，对实测地震记录进行分析。可以看出，小波包空间分解的子频带的结点序号并不是完全按照频率递增顺序排列的，直观上证实了小波包处理地震非平稳信号过程中存在的“跳频”现象。频带序号对应结点序号的结果与现有方法的排列规律完全一致，验证了本文所述频带序号到结点序号转换算法的正确性。若不考虑能量泄漏，各结点通频带范围与本文所提供的转换算法结果基本一致，再次验证了本文所述转换算法的正确性。该转换算法对快速高效地提取地震信号局部特征、有层次地展现地震信号等具有重要意义。