0 引 言

旋转机械的运行状态信息十分丰富，在故障发生或者发展过程中存在着大量的复杂、非平稳信号。其关键部件转轴和轴瓦发生摩擦时，金属部件间会产生瞬态冲击力，从而出现周期性冲击信号。但是，故障早期的冲击幅值很小，并且不可避免地受到各种噪声干扰。传统小波虽具有一定的时频分析能力，但其平移不变性差，提取碰摩类高频信号能力弱，无法准确兼顾信号局部特征^[1]，常常导致故障信息遗漏。

双树复小波变换(Dual-tree complex wavelet transform，DTCWT)最早由英国剑桥大学Kingsbury教授提出，它不仅具有优良的时频局部细化功能，同时降噪性能^[2]和抗频率混叠特性均很突出^[3]。Y.X.Wang^[4]针对旋转机械的混合故障诊断开展研究，证明双树复小波可以取得比提升小波和快速峭度图更好的分析效果，且DTCWT的高效率使得它能胜任机械设备的在线诊断。吴定海等^[5]提出了基于双树复小波包变换和自适应分块阈值降噪的能量提取方法，并应用于柴油机缸盖振动信号的特征提取，效果明显。胥永刚^[6]、艾树峰^[7]和苏文胜等^[8]通过对双树复小波分解后的信号进行包络谱分析，实现了对滚动轴承信号的故障特征提取。

笔者提出将双树复小波应用于转子等关键部件的故障诊断，对振动信号进行多层分解、细化，并利用迭代奇异值分解降噪(Iterative singular value decom-position，ISVD)^[9]代替传统的软、硬阈值处理方式，从而将碰摩产生的微弱周期性冲击振动信息提取分析，以更好地凸显机械故障特征。仿真结果和试验结果表明，该方法可以有效提取转子微弱故障的特征信息。

DTCWT是基于实数小波实现的复数小波变换，它通过2个并行的、高度对称互补的实数小波滤波器组来进行信号的分解与重构，确保了波形幅值和相位信息能够完整补偿，避免了传统的离散小波分解时，由于操作引起的频率混叠和奇异点敏感问题^[10]。

现对长度为N的振动信号X(n)进行双树复小波分解。以a树、b树分别代表DTCWT的实部和虚部。应用时需要合理设计a树和b树的低通滤波器H₀(p)和G₀(p)，使其满足半采样延迟条件^[11]，即b树的采样位置点始终位于a树中间，从而能够综合利用2树的分解系数，并且振动信号经变换后将获得长度为2N的分解系数，信息是2倍冗余，不会造成有用信息缺失。

令x(t)表示振动信号，ψ_h(t)和ψ_g(t)分别表示2个实小波，相对应的尺度函数分别为δ_h(t)和δ_g(t)，则双树复小波可以表示为：

实部树分解：根据小波理论，DTCWT的实部树小波系数和尺度系数可通过下式计算：

虚部树分解：将式(2)中的ψ_h(t)和δ_h(t)分别用ψ_g(t)及δ_g(t)代替，即可得到虚部树分解的小波系数

由式(2)和式(3)中的实、虚部树的输出系数，联合计算可得到双树复小波分解的小波系数和尺度系数：

DTCWT的重构过程为分解的逆过程，由末层依次进行插值、滤波和加法运算即可实现。对式(4)中得到的各层小波系数和末层尺度系数分别进行单支重构，则对应得到DTCWT分解的各层细节信号及末层近似信号。各层小波系数单支重构：

末层尺度系数单支重构：

重构信号x—(t)等于所有单支重构信号之和：

合适的滤波器可有效提高双树复小波的分解精度，并能较大程度地减小重构误差。相反，若构造的滤波器光滑性偏差，则分解信号将达不到预定的滤波效果。传统Odd-even滤波器组长度奇偶交替，这种设计是为了保证线性相位要求采用双正交小波，造成2树具有不同的频响特性，

影响分解系数精度。为克服此不足，N.G.Kingsbury^[12]提出Q-shift滤波器组的计算方法，不再强调线性相位，且在第2层及更高层上的设计长度均为偶数。

Q-shift滤波器组具有 1/4 采样周期的群时延，a树为 1/4 采样时延(q)时，b树为 3/4 时延(3q)。通过这种时延来保证在每一层滤波时，b树的采样点刚好位于a树采样点中间，保证采样过程的对称性。具有 1/4 采样延迟的Q-shift滤波器组应满足以下Z变换公式^[12]：

H₀=[0.003 3-0.003 90.034 7-0.038 9-0.117 20.275 30.756 10.568 80.011 9-0.106 70.023 80.017 0-0.005 4-0.004 6]；

G₀=[-0.004 6-0.005 40.017 00.023 8-0.106 70.011 90.568 80.756 10.275 3-0.117 2-0.038 90.034 7-0.003 90.003 3]；

表1形象地表示了输入点数为15的原始信号在第1层和第2层上的Q-shift滤波器组的采样结构。

假设转子振动信号x(t)只包含4种频率成分(30 Hz，180 Hz，380 Hz，600 Hz)，表达式为：

为证明DTCWT的优良抗频率混叠特性，对该信号进行3层DTCWT分解，单支重构后得到的各层时域信号的对应频谱如图1所示。对应特征频率从低到高依次为30、180、380和600 Hz，可见各频率信号被准确、无混叠地分布在理想频带中。

图1 DTCWT分解信号频谱图 Fig.1 Decomposed signal spectrogram of DTCWT

图选项

而使用传统Db4小波对该信号进行3层分解，得到的频谱图出现了明显的频率混叠，具有明显瑕疵，如图2所示。其中第2层细节信号d2理论上只有380 Hz的信号，但它实际上还包含原信号中不存在的频率成分320、400和620 Hz，即由小波分解时其他频率段与此频率段重合调制所导致。

图2 DWT分解信号频谱图 Fig.2 Decomposed signal spectrogram of DWT

图选项

振动信号经DTCWT多层分解后，蕴含各类频率信息的频带信号就得到了有效区分，可更准确地提取故障特征，但是波形与频谱上仍会显示由噪声引起的瑕疵和无用频率信息，影响分析效果。目前，Donoho等^[13]提出的软阈值降噪方法应用广泛,表达式为：

式(10)中阈值的选择对降噪结果影响很大：阈值太小，重构后的信号中仍然含有噪声；而阈值选择过大，信号中的有用信息将被滤除，导致信号失真。实际操作中该值需要反复尝试和验算。

鉴于此，笔者选用迭代奇异值分解合理选择降噪阶次，滤除噪声并保留各信号有用信息。给定一个秩为r的m×n维矩阵X，则X的奇异值分解形式为^[14]：

对A进行奇异值分解，得到(λ₁，λ₂，λ₃，……，λ_r)为A的奇异值。矩阵A的奇异值可以反映信号和噪声能量集中的情况，并且已经证明噪声将会使所有的奇异值都大于0。若矩阵秩为r，则前r个奇异值要远大于其他奇异值。将较大的前k个奇异值(k为降噪阶次)保留，其余置0，然后利用奇异值分解的逆过程得到矩阵A的最佳逼近矩阵A_r。将矩阵A_r的各列对应相加取平均值，即可得到降噪后的时间序列。由于A_r仅是对A的一次逼近，得到的时间序列可能没有完全剔除噪声，可以重复迭代以上各步，直到满足要求为止。信号降噪阶次k根据奇异熵增量法确定^[15]。

下面给出基于DTCWT-ISVD的振动信号降噪具体步骤：

(1)利用DTCWT对原始信号进行多层分解，并对分解后的各层小波系数和末层尺度系数分别进行单支节点重构，则得到DTCWT分解的各层细节信号及末层近似信号。

(2)按照SVD规则，对DTCWT后的各个信号进行分解，并根据奇异熵增量曲线确定降噪阶次，一般选择当曲线开始下降并趋向于一个稳定值时的阶次，防止因“一刀切”的阶次造成的高频信息中有用部分被错误滤除。

(3)重复上述步骤，对各信号进行迭代分解降噪。

(4)将降噪后的细节信号和近似信号全部叠加，即可获得降噪后真实的振动信号。

为了验证降噪方法的优越性，将该方法用于2种不同仿真信号的噪声处理，并与SVD、Db4等方法进行对比。

信号1为一叠加白噪声的heavy sine信号，信噪比5.3 dB。信号2为一个包含3种频率成分(30、150和360 Hz)的余弦信号和白噪声叠加的多谐波信号，信噪比5.3 dB，表达式为：

为了对降噪效果进行定量评价，引入2个降噪评价指标，即信噪比R_SNR和均方差E_MSE：

4种不同方法的降噪效果如表2所示。由表可见，基于DTCWT-ISVD的降噪方法针对以上2种不同仿真信号均获得了最高的信噪比和最小的均方差，获得了较好的降噪效果，有利于后续的分析。

Bently RK4转子故障模拟试验系统由调速电机、转轴、双圆盘、滑动轴承和支架等组成。通过在支架上安装碰摩螺栓可制造出转轴与螺栓的摩擦效果。数据采集仪采用中国石油大学(北京)故障诊断实验室自主开发的MDES-5型故障诊断系统。试验中，将加速度传感器吸附在轴承座表面，用数据线将传感器、采集仪和电脑连接起来。

转子振动信号采集示意图如图3a所示，碰摩螺栓安装图如图3b所示。试验参数：转速3 000 r/min；采样频率16 kHz，采样点数为4 096。

图3 转子故障模拟及分析系统 Fig.3 Rotor fault simulation and analysis system

图选项

采集的故障信号时域波形与频谱如图4所示。由图可以看出，该信号时域波形与转子正常振动波形(见图5)相差较大；而频谱图上出现了基频50 Hz、5倍频250 Hz以及一些高次谐波分量，与碰摩和松动等故障特征接近。但是谱线混叠严重、辨识度低，无法准确判断转子到底发生何种故障，下面基于DTCWT-ISVD法对振动信号进行分析。

图4 转子故障振动信号 Fig.4 Rotor fault vibration signal

图选项

图5 转子正常振动信号 Fig.5 Normal rotor vibration signal

图选项

用笔者的方法对转子故障振动信号进行4层分解，并对分解得到的4个细节信号进行奇异值分解降噪。其中降噪后第2层重构细节信号如图6所示。从图可见，每隔0.02 s就发生1次冲击，与转动周期(1/50 s)相同，即认为转子每转1周，转轴与螺栓发生1次冲击摩擦。

图6 故障信号经DTCWT-ISVD后的第2层细节信号 Fig.6 The second layer of detail signal from decomposition of fault signal by DTCWT-ISVD

图选项

试验结束后，对Bently转子试验台的转轴进行检查，发现碰摩螺栓与转轴间隙极小，运转时会产生周期性的接触，端面也被不均匀地磨损，因此产生了图6中的周期性冲击信号。

作为对比，使用传统小波分析对转子故障振动信号进行2 层分解，小波基函数为Db4。对第2 层的细节信号进行硬阈值处理，结果如图7所示。从图中看不出明显的周期性冲击现象。这是由于冲击幅值较小，约为0.05 m/s²，且冲击周期与转动周期相同，因而被淹没在正常周期信号之中，冲击特征不明显。

图7 故障信号经DWT分解后第2层细节信号 Fig.7 The second layer of detail signal from decomposition of fault signal by DWT

图选项

(1)通过构造的基于Q-shift滤波器组的双树复小波对振动信号进行多层分解，并综合利用复小波分解系数信息，可获得频率特征划分清晰的频带信号，从而凸显故障信息。与传统小波相比，双树复小波构造巧妙，并且具有平移不变性、抗频率混叠和可完全重构等优点。

(2)双树复小波变换和迭代奇异值分解相结合的分析方法比传统阈值降噪方法具有更好的降噪效果。信号经处理后，不仅能获得最大信噪比和最小均方值，还能够保证有用冲击成分不会被错误地滤除。

(3)利用双树复小波变换故障诊断方法成功提取了转子碰摩时产生的微弱冲击性信号，表明该方法故障特征提取作用明显。研究结果为现场设备转子故障诊断奠定了基础。