简介这是一份依据西安交通大学已故屈梁生院士全息谱理论编写的Matlab源程序包面向从事旋转机械故障诊断研究的工程师、研究生及相关领域学习者。全息诊断技术融合幅值、相位与频率信息可有效识别转子不平衡、碰摩等常见故障。压缩包共39个文件核心为M脚本及配套的DV、DH系列数据文件另有DH1/DV1等扩展数据与HLV/HLH辅助文件整体仅58KB便于快速部署与运行。程序内附典型故障仿真数据并提供二维/三维全息谱绘图入口方便对照学习已有868人学习使用。通过源码与示例数据读者可掌握全息谱图的构建思路、二维/三维全息谱绘制方法并可将程序迁移至自己的振动信号中节省算法复现时间。 做旋转机械振动诊断这些年我接触过不少同行大家普遍的感受是光看频谱图很多时候只能知道“哪个频率的振动大”却说不清楚“转子在这一频率下到底是怎么运动的”。频谱把幅值留下了却把两路信号之间的相位关系丢掉了。全息谱这个工具恰恰是把同一截面上两个互相垂直测点的幅值和相位融合起来在频域里画出一个椭圆轨迹让人一眼看出转子的进动方向和振型特征。这套原创的经典全息谱Matlab源程序就是围绕“全息诊断”这条主线写的从数据预处理、FFT幅值相位提取到椭圆参数合成、二维和三维全息谱绘制再到典型故障的特征判别全部打通。不管你是刚开始接触全息谱的研究生还是现场做故障诊断的工程师或者想把振动分析算法落地的开发人员这套程序和这篇文章都能提供一个可以“抄作业”的参考。1. 全息谱到底解决了什么问题——从“看频谱”到“看轨迹”1.1 单通道频谱分析的天然局限转子系统在运行中我们通常在轴承座的同一截面安装两个电涡流位移传感器方向互成90度习惯上叫X方向和Y方向。对这两个通道分别做FFT能得到各自的幅值谱和相位谱。问题在于单通道谱图只描述了“这个方向上有多少振动能量”它丢失了关键的空间信息。举个例子同一台机器1X分量幅值都是80微米X方向和Y方向的相位差是30度还是150度在两张独立的频谱图上是看不出本质区别的。但实际上前者可能对应正进动为主的椭圆轨迹后者可能已经接近反进动故障性质完全不同。这就是我常和同事说的频谱是“一维视角”而转子在空间里的实际运动是“二维轨迹”少看一个维度就少了一半信息。1.2 全息谱的“全息”到底在哪里全息谱的核心思路并不复杂把X、Y两个通道的时域信号分别做FFT在每一个关注的频率分量比如转频的1倍频、2倍频上取出该频率对应的幅值和相位然后把这两组幅值和相位合成为一个椭圆。为什么合成出来一定是椭圆因为同一频率下X方向的振动是A_x·cos(ωtφ_x)Y方向是A_y·cos(ωtφ_y)消去时间变量t之后这两个简谐振动的轨迹在数学上就是一个李萨如椭圆。这个椭圆的形状、大小、倾斜角度、偏心率以及进动方向包含了两个通道幅值和相位的全部信息。所以全息谱的“全息”本质上是把两路信号的幅值和相位信息压缩进一个几何图形里信息没有丢失反而更直观了。1.3 这套源程序能输出什么、适合谁来用这套Matlab源程序不是只画一个椭圆就完了它围绕全息诊断做了完整的链路一维全息谱在选定倍频处绘制该频率下的振动椭圆输出长半轴、短半轴、倾角、偏心率、进动方向等参数。这是最核心的诊断依据。二维全息谱把多个倍频0.5X、1X、2X、3X等的椭圆按频率轴排列在同一张图上横轴是频率纵轴上每个椭圆代表该频率下的振动形态。三维全息谱在二维全息谱的基础上把不同转速下的椭圆轨迹层叠起来形成“转速-频率-振动形态”的三维图谱启停机过程的分析特别有用。故障特征判别程序内置了几类典型故障的全息谱特征库比如不平衡、不对中、油膜涡动、动静碰磨可以辅助人工判断。适合的使用者很明确做转子动力学研究的在校学生省去从零啃理论到写代码的过程现场诊断工程师拿现场数据跑一遍几分钟就能看到全息谱图还有做状态监测系统开发的工程师这套程序的模块化写法可以直接嵌入自己的平台。2. 数据流的起点从信号读取到整周期截断的工程细节2.1 程序输入数据的组织方式源程序里约定了一套输入格式建议你按这个来整理数据调试起来会顺手很多。第一路是X方向位移信号第二路是Y方向位移信号第三路是键相脉冲或者转速信号。如果现场采集系统没有单独记录键相也可以用转速计通道程序里会自动换算成转频。采样率这个参数必须单独传入不能从数据里猜。很多现场数据标称采样率是2560Hz或者5120Hz但实际采集卡存在时钟偏差如果你在程序里写死采样率而实际值偏了0.1%整周期截断之后相位误差可能已经有几度了。全息谱对相位很敏感所以我的习惯是采集时同步记录采样率实测值或者用键相脉冲反推实际采样率。2.2 整周期截断为什么是“生死线”全息谱算法最基础也最要命的一步是截取整周期数据段。所谓整周期就是数据段长度包含整数个转频周期比如10转、20转。这样做的好处是FFT后转频及其倍频正好落在某根谱线上没有频谱泄漏幅值和相位都能准确还原。如果截取的不是整周期哪怕只差一个采样点转频能量就会泄漏到相邻谱线提取出来的幅值偏小相位还会出现偏移。更麻烦的是泄漏产生的旁瓣可能掩盖真实的小幅值故障特征比如油膜涡动的0.5X分量本来就不大被泄漏一污染诊断结论可能直接翻车。程序里我是这样处理的先根据转速信号计算转频f0然后根据采样率fs和目标周期数T计算截取点数N round(fs * T / f0)。比如转速1500rpm转频25Hz采样率2560Hz想截10个周期点数就是round(2560*10/25)1024点。然后从数据流里取连续N个点。转速波动小的机组按固定点数截取就够了。2.3 转速波动时怎么办现场工况不会那么理想特别是启停机过程中转频是变化的。这时候固定点数截取整周期就不再可靠。程序里预留了一个处理接口可以按相邻两个键相脉冲之间的采样点数作为“当前周期长度”然后拼出N个完整周期。严格一点的做法是阶比跟踪重采样把等时间间隔采样的信号转成等角度间隔的角域信号再做FFT。这个功能在程序里是作为一个可选模块实现的默认不开启因为阶比跟踪对键相信号质量要求高键相毛刺多的话反而引入新误差。我的建议是稳态工况用固定整周期截断程序默认路径简单稳定变转速工况再启用重采样模块同时把转速通道的滤波做好。3. 核心算法实现FFT幅值相位提取与椭圆参数合成3.1 从FFT谱线中准确提取幅值和相位整周期截断之后对X、Y通道分别做FFT然后找到关注的频率索引。假设采样率fsFFT点数N频率分辨率dffs/N那么第k倍频对应的谱线索引是round(k*f0/df)1。提取幅值和相位的Matlab代码很简洁但有几个细节必须注意。相位提取时用的是atan2(imag, real)这样能得到(-pi, pi]区间的相位值不会出现arctan象限判断错误。同时要把FFT结果的幅值换算成真实幅值单边谱要乘以2再除以N。function [A, phi] extract_component(data, fs, f_target) N length(data); df fs / N; k round(f_target / df) 1; spec fft(data) / N; % 取单边谱乘2恢复真实幅值 A 2 * abs(spec(k)); % 相位单位弧度相对数据段起点 phi atan2(imag(spec(k)), real(spec(k))); end相位基准是所有全息谱程序最容易出错的地方。FFT算出来的相位是相对于整段数据起点的如果X通道和Y通道的数据起点没有严格对齐或者数据采集时两通道之间存在时间延迟相位差就会失真。所以程序在读取数据时就强制校验两个通道的采样点数一致起止时刻一致。现场同步采集一般没问题但如果是把离线数据拼凑起来的这一步必须手动检查。3.2 椭圆方程推导与几何参数计算拿到X通道的幅值A_x、相位φ_x和Y通道的幅值A_y、相位φ_y之后把这两个方向的简谐振动合成轨迹满足如下椭圆方程[ \frac{x^2}{A_x^2} \frac{y^2}{A_y^2} - \frac{2xy\cos(\Delta\varphi)}{A_x A_y} \sin^2(\Delta\varphi) ]其中Δφφ_y-φ_x。只要Δφ不是0或±π轨迹就是一个有面积的真椭圆Δφ接近0或±π时椭圆退化为一根直线说明该频率下振动方向比较固定。Matlab里要算椭圆几何参数最稳妥的做法是用参数采样加主成分分析比纯公式推导更容易实现而且代码不容易出错。思路是在0到2π之间均匀取360个角度点用参数方程算出椭圆轨迹点云然后对点云做主成分分析第一主成分方向就是长轴方向第二主成分方向就是短轴方向两个方向上的投影极差就是长短半轴。function ell fit_ellipse(Ax, phix, Ay, phiy) theta linspace(0, 2*pi, 360); x Ax * cos(theta phix); y Ay * cos(theta phiy); pts [x, y]; % 椭圆中心在原点直接做主成分分析 coeff pca(pts); proj pts * coeff; a max(proj(:,1)) - min(proj(:,1)); b max(proj(:,2)) - min(proj(:,2)); % 长半轴和短半轴 if a b ell.a a / 2; ell.b b / 2; ell.angle atan2(coeff(2,1), coeff(1,1)); else ell.a b / 2; ell.b a / 2; ell.angle atan2(coeff(2,2), coeff(1,2)); end ell.ecc sqrt(1 - (ell.b / ell.a)^2); % 进动方向根据相位差判断 dphi phiy - phix; dphi mod(dphi pi, 2*pi) - pi; if dphi 0 ell.precession forward; else ell.precession reverse; end end这里进动方向的判断逻辑值得展开说。Δφ在(0, π)区间时Y方向振动在相位上超前X方向轨迹沿正方向旋转定义为正进动Δφ在(-π, 0)区间时为反进动。正进动常见于不平衡激励反进动往往和转子碰磨、某些流体激励有关所以这个参数对故障定性非常有价值。但要注意当Δφ接近0而幅值又差不多时椭圆趋近直线长轴方向不稳定进动方向意义也不大程序会输出退化提示。3.3 二维和三维全息谱的绘图逻辑二维全息谱画起来不复杂但要画得清楚需要处理好坐标缩放。不同倍频分量的幅值可能差很多比如1X是100微米2X只有10微米如果所有椭圆都用同一个坐标尺度2X椭圆会小到看不见。程序里默认对每个椭圆做独立归一化把长半轴归一化到同一显示尺寸同时标注真实长半轴数值。这样图上每个倍频的椭圆形态都能看明白数值也不丢失。三维全息谱本质上是把不同转速下的二维全息谱沿转速轴堆叠。视觉上用surf或者patch都可以画关键是控制透明度和视角。转频升速时正进动和反进动椭圆如果用同一种颜色三维图很容易糊成一片。我的程序里按进动方向着色正进动用暖色系反进动用冷色系这样图谱的辨识度明显提升。4. 调试这套程序时最容易翻车的四个环节4.1 频谱泄漏伪装成“故障特征”第一次用仿真数据验证时我最开始图省事没有加整周期截断直接拿一长段数据做FFT。结果1X旁边冒出一串边带2X幅值虚高我当时一度以为信号里存在调制故障排查了半天才发现是截断位置不对。这个问题最坑的地方在于它不会报错程序正常跑完输出看起来“很有内容”实际上全是泄漏造成的假象。解决方案就一句话整周期截断比任何窗函数都有效。理论上加汉宁窗也能抑制泄漏但窗函数会改变相位而全息谱恰恰对相位敏感。所以在这套程序里我坚持用矩形窗加整周期截断的组合而不是加窗。你如果非要用加窗的方式务必测试窗函数对相位差的影响不要默认“加窗不影响相位差”。4.2 相位基准不统一导致椭圆变形除了FFT本身的相位定义还要警惕一个坑很多采集系统对每个通道都有抗混叠滤波器不同通道的滤波器如果器件参数有偏差相位响应就不一致引入通道间附加相移。这个相移跟频率相关1X处差两三度2X处可能差十几度。现场验证方法是把同一个信号同时接入两个采集通道看看全息谱上是否出现了一个“本不应该存在”的椭圆。理想情况下同相信号合成的椭圆应该退化成一条直线如果画出来不是直线说明通道间存在相位失配。程序里提供了一个标定通道相位差的模块可以在分析前对相位差做补偿。4.3 转速波动导致倍频位置错位转速波动不光影响整周期截断还影响FFT谱线定位。比如目标转频25Hz频率分辨率2.5Hz转速波动0.5%那就是0.125Hz看起来不大但对应到谱线上已经偏移了半根谱线。这时候你用固定索引去提取幅值相位得到的幅值会偏小相位会漂移。程序里解决的办法是提取复数值时不是只看目标谱线而是在目标谱线附近做一个幅度加权重心计算相当于频域内插。这个办法在信噪比高的场合下效果很好计算量也不大。如果转速波动超过1%就别硬用定转速算法了老老实实开阶比跟踪模块。4.4 预处理滤波给相位“埋雷”很多人在数据分析前习惯先做带通滤波去掉高频和低频干扰。这个习惯本身没错但普通IIR滤波器是因果系统不同频率分量经过滤波器后相位延迟不一样这就改变了X、Y通道之间的相位关系。全息谱对相位差的敏感度是“度”级别的滤波器引入的相移可能达到几十度画出来的椭圆完全变形。如果必须在全息谱前做滤波我建议用零相位滤波。Matlab里就是filtfilt函数它把信号正着滤一遍再反着滤一遍相位偏移互相抵消。当然代价是计算量翻倍而且要求数据是离线处理。在线系统中没法用filtfilt那就在滤波之后做相位补偿但这需要对滤波器相位响应有精确的模型复杂度高不少能不做就不做。5. 源程序怎么验证、怎么改造成自己的诊断工具5.1 用一组仿真信号快速验证程序正确性拿到源程序后第一步不是拿现场数据跑而是用仿真信号做“盲测”验证算法本身的正确性。构造一组已知参数的信号转频25HzX方向幅值80微米、相位0.3弧度Y方向幅值60微米、相位1.2弧度采样率2560Hz截取10个周期。fs 2560; f0 25; N round(fs * 10 / f0); t (0:N-1) / fs; Ax 80; phix 0.3; Ay 60; phiy 1.2; x Ax * cos(2*pi*f0*t phix); y Ay * cos(2*pi*f0*t phiy) 5 * cos(2*pi*2*f0*t 0.8);程序跑完后1X椭圆提取结果应该与设定的参数高度一致。Y方向加了2X小分量全息谱的2X位置也应该出现一个小椭圆。如果这个基础验证都过不了先回头检查整周期截断和幅值相位提取别急着分析现场数据。5.2 现场数据的典型非理想情况与判断原则现场数据比仿真脏得多噪声大、存在不对中引起的多倍频耦合、传感器安装角度可能不是严格90度、两个通道灵敏度标定不一致。这些因素都会让全息谱椭圆出现变形甚至“畸形椭圆”。遇到畸形椭圆我的原则是先不急着下结论。首先检查输入数据是否有毛刺、缺包和磁干扰其次对比两个通道的时域峰峰值确认没有通道信号饱和最后看多个稳定工况下的椭圆是否重复出现。如果椭圆形态在不同转速下稳定复现才认为它反映了真实的转子运动特征。单凭一次测点数据画出的某个椭圆就下诊断结论很容易被偶然因素带偏。5.3 典型故障的全息谱特征对照程序里内置的故障判别逻辑主要依据下表这些经验特征。需要提醒的是全息谱给出的是“故障的可能性方向”不是绝对的判据现场诊断要结合工艺量和历史趋势一起看。故障类型主导倍频全息谱典型特征进动方向转子不平衡1X1X椭圆长轴稳定幅值随机组转速变化明显正进动轴系不对中2X2X椭圆明显常伴随1X椭圆增大椭圆偏心率大正进动为主油膜涡动0.4X~0.5X椭圆出现在半频附近形态不稳定常见反进动动静碰磨1X及高倍频1X椭圆相位不稳定多倍频椭圆出现可能出现段进动反向支承松动1X、2X、3X高倍频椭圆成分增多整体形态杂乱方向随机5.4 把程序改造成批量诊断工具的几点建议这套源程序本身是面向单个数据分析的实际工程使用中我建议做几个改造。第一把数据读取函数参数化改成输入是数据目录、测点编号、转速信息的批量接口。第二把椭圆参数、进动方向、主导倍频这些特征导出成结构化数据方便做趋势统计。第三增加报告生成模块自动输出全息谱图和特征表减少人工截图整理的工作量。值得多说一句的是不要把全息谱和传统频谱对立起来。实际诊断时我是先用频谱看全貌锁定异常频率再用全息谱对异常频率做进动分析和振型判断。两者配合比单独用任何一个都稳得多。在我实际使用中体会最深的一点是全息谱这个工具真正有价值的地方不在于把椭圆画得多好看而在于它逼着你把X、Y两个方向的振动信息放在同一个坐标系里去思考。幅值大不大是一回事振动方向怎么走、相位关系稳不稳定是另一回事。建议你拿到源程序后先花半天时间用仿真信号把每个环节跑熟体会到相位差对椭圆形态的影响之后再拿现场数据实战那时候你对全息诊断的理解就不只是会调用函数画图了。本文还有配套的精品资源点击获取