功率谱密度PSD实战：随机振动试验谱与故障诊断校核-平芜编程栈

功率谱密度，也就是大家嘴里常说的 PSD，我第一次把它当成日常工具而不是课本上的一行公式，是从一次台架试验的扯皮开始的。当时供应商甩过来一条随机振动谱，要求按这条谱做几个小时的耐久考核，现场有人问这套试验的总量级是多少 g，另一位同事张口就来：峰值 0.04 g²/Hz，带宽 20 到 2000 Hz，乘一下大概 79，开根号差不多 8.9 g。听起来挺像那么回事，但这套算法的前提是谱线在整个带宽里是平的，而那条谱从 20 Hz 爬到 80 Hz 是上升段，从 500 Hz 往后又是一路往下掉的下降段，直接乘带宽得到的数字虚高了一大截。说白了，PSD 描述的是“能量在频率上怎么分布”，要拿到总量级，得沿着频率轴做积分，而不是拿一个点去乘一段宽度。

这件事之后我把 PSD 从头捋了一遍，发现它其实横跨了两个看起来完全不搭界的场景：一个是产品环境适应性这一侧的随机振动分析与试验谱制定，另一个是设备状态监测这一侧的故障诊断。前者关心的是“我这条谱能不能代表真实环境、试验台能不能复现出来”，后者关心的是“机器某个部位是不是已经开始不行了”。两边的数学内核是同一套，但工程上的坑完全不一样。下面这些内容适合刚接触振动与信号分析的朋友打底，也适合已经会用软件点两下出图、但对自己画出来的曲线心里没底的人对照着看，我会把参数怎么定、结果怎么校核、哪里最容易翻车都摊开讲。

1. 功率谱密度到底在描述什么，为什么随机振动非它不可

1.1 从一个“谱线不平”的现场问题说起

接着上面那个例子往下拆。随机振动试验谱通常用一张双对数坐标图表示，横轴是频率，纵轴是加速度谱密度，单位 g²/Hz，图上由若干段直线连成折线，每段直线的斜率用 dB/oct 标注。这张图表达的意思很朴素：在某个频率附近，每 1 Hz 带宽里分摊了多少加速度的均方值。注意是均方值，不是幅值，也不是峰值，这是理解 PSD 的第一道门槛。

为什么随机振动不能像正弦振动那样直接给一个“加速度幅值”？因为随机振动的瞬时值一直在变，你没法用“这一刻是多少 g”来描述它，只能用统计量。同一台设备在同一个台上跑两次，得到的时域波形完全不一样，但如果你把两次波形各自算一遍 PSD，两条曲线在统计意义上会高度重合。这个性质非常关键——随机信号的具体波形不可重复，但它的功率谱密度是可重复的，所以工程上用 PSD 来定义随机振动环境，而不是用时域波形。

注意：PSD 是统计量，它描述的是“大量样本平均之后的能量分布”。单次短时数据的 PSD 起伏会很大，看着像一堆毛刺，那不是噪声，那是估计方差，后面第 2 节会讲怎么压下去。

1.2 确定性信号看频谱，随机信号看密度

把一个确定性信号比如单频正弦做 FFT，得到的是两根尖峰，幅值有明确物理意义，可以读成“这里有 1 g 的加速度”。但对随机信号做 FFT，每一帧的结果都不一样，幅值本身没有稳定意义，只有取模平方再按时间或按样本做平均，才收敛出一个稳定的量。这就是 PSD 与普通幅值谱最本质的区别：幅值谱的纵轴是幅值，PSD 的纵轴是方差密度，前者可以逐点读，后者必须积分才有物理意义。

生活化的类比：一条高速公路上有多少车，这是“总量”；每公里路段上平均有多少车，这是“密度”。你只知道密度，得乘以路段长度才知道总量。PSD 就是“每赫兹有多少加速度方差”，乘上带宽才是方差，开根号才是我们平时说的 Grms。很多人算错就错在把密度当总量用了。

从数学上说，单边 PSD 满足一个非常实用的关系：信号均方值等于 PSD 在正频率轴上的积分。写出来就是 σ² = ∫S(f)df，两头开根号就是 Grms = sqrt(∫S(f)df)。这个关系在工程上至少有三种用法：一是从试验谱反算总量级，判断试验台推力够不够；二是用实测数据算出的 PSD 反推时域均方值，作为算法正确性的自检；三是在故障诊断里比较不同频带贡献了多少能量。

1.3 量纲这件事，比想象中更容易出错

g²/Hz 这个单位，拆开看是“加速度单位的平方，除以频率单位”。换成国际单位就是 (m/s²)²/Hz，两者相差 9.80665² ≈ 96.17 倍。这个系数在跨部门协作时经常出事：试验方给的谱用的是 g²/Hz，分析人员按 m/s² 直接读进去算，最后算出来的 Grms 差了将近 10 倍，然后互相甩锅说对方数据有问题。

还有一种情况是同一份数据被标注成 ASD（加速度谱密度）和 PSD 两个名字，其实是一回事，只是叫法习惯不同。速度谱密度、位移谱密度也类似，单位分别是 (m/s)²/Hz 和 m²/Hz，在低频段分析时用位移谱更直观，因为低频加速度幅值小，画在图上很容易被噪声淹没，转成位移谱能看得更清楚。转换关系靠的是对频率的幂次乘法：加速度谱乘 1/(2πf)² 得速度谱，再乘一次得位移谱。转换的时候别忘了，这个操作在 f 接近 0 的地方会爆炸，所以低频端一定要设一个截止频率，不然算出来的位移谱全是奇奇怪怪的尖峰。

2. 算法实现路线怎么选，参数怎么定才不吃亏

2.1 直接法（周期图）为什么不能拿来直接用

最朴素的做法是把整段数据做一次 FFT，取模平方，除以采样率和点数，得到一条 PSD 曲线。这个方法叫周期图法，理论上完全正确，实践上几乎不能用。原因在于它是一个“不一致估计”：随着数据长度增加，频率分辨率越来越好，但每个频点上的估计方差并不会变小，标准差始终维持在均值同一量级。换句话说，你采再多的数据，图上那堆毛刺还是那么高。

我最早踩这个坑是在分析一段几十秒的实测振动数据时，画出来的曲线上下跳动了好几倍，同一台设备换一段时间重算，峰值位置都能挪一点。当时以为是设备状态不稳定，折腾了一整天才反应过来是估计方法的问题。分辨率高不等于结果可信，这是两个维度的事。

解决办法是分段平均：把长数据切成若干段，每段算一个周期图，然后对所有段做平均。段数越多，方差越小，代价是每段变短，频率分辨率变差。这就是经典的 Welch 法，本质是用分辨率换稳定性。

2.2 三个关键参数：段长、重叠率、窗函数

Welch 法的参数就三个，但每一个都直接影响结果能不能用。

段长（nperseg）决定频率分辨率，Δf = fs/nperseg。要分辨间隔 5 Hz 的两根谱线，段长至少要 fs/5 个点。工程上的经验是：先确定你需要分辨的最小频率间隔，再倒推段长，而不是拿一个 1024 点当默认值到处套。低频分析尤其要注意，20 Hz 附近的问题用 1024 点、采样率 2000 Hz 去算，分辨率接近 2 Hz，勉强够用；如果要看 5 Hz 以下的成分，段长得拉到 8192 甚至更长。

重叠率（noverlap）是为了在段数有限的情况下多榨出一些平均次数。重叠 50% 是行业默认值，因为再往上增益递减，计算量却线性增长。实际效果是：不重叠时，1024 点数据能切出 20 段，重叠 50% 能切出将近 40 段，方差直接降下来一截。不过重叠段之间不再完全独立，所以有效自由度不会按段数线性增长，一般打个九折左右。

窗函数的选择要看你关心什么。汉宁窗（Hann）是最通用的，旁瓣衰减快，泄漏小，适合绝大多数随机信号分析；平顶窗（Flattop）幅值精度高，适合做标定和读绝对幅值，但主瓣宽，分辨率差；矩形窗分辨率最高，但泄漏严重，只在信号严格整周期截断时用。随机振动分析一律建议加汉宁窗，因为随机信号本身不存在整周期截断的可能，不加窗必然泄漏。

提示：加窗会降低信号能量，不同窗的能量修正系数不同。汉宁窗的幅值修正系数约 2.0，能量修正系数（ENBW）约 1.5。用现成的库函数时这些通常已经内部处理了，但如果你自己手写 FFT，一定要记得除，不然算出来的 PSD 会整体偏低。

2.3 手算一遍：采样率、段长、分辨率、自由度

举个例子把参数定死。某设备振动数据，采样率 fs = 2000 Hz，总时长 60 s，总点数 120000。关心的是 10 到 500 Hz 范围内的谱形，希望能分辨间隔约 5 Hz 的谱峰。

先定段长：Δf = fs/N ≤ 5，所以 N ≥ 400，取 2 的整数次幂，N = 512，对应 Δf = 3.9 Hz，够用。想更细一点取 N = 1024，Δf = 1.95 Hz，也完全可以，代价是段数变少。

再算段数：重叠 50%，每段前进 N/2 个点，段数 K ≈ (120000 - 1024)/512 + 1 ≈ 233 段。这个数量很充裕，方差会压得很低。

算自由度：不重叠时 dof ≈ 2K；50% 重叠打折，取 dof ≈ 1.9K ≈ 443。标准不确定度约为 1/sqrt(443) ≈ 4.8%。也就是说，你算出来的每个频点上的 PSD 值，误差大概在 5% 以内，这个精度做工程判断完全够了。

反过来看，如果总数据只有 5 s，即 10000 点，同样参数只能切出 18 段左右，dof ≈ 34，误差约 17%，图上就会明显毛糙。这时候要么接受毛糙、要么拉长采集时间、要么牺牲分辨率用更短的段长。这三者只能取其二，没有第四条路，这是做谱分析必须接受的基本约束。

2.4 单边谱、双边谱，以及那个容易漏掉的 2 倍

理论上功率谱密度定义在 −∞ 到 +∞ 上，是双边谱，正负频率各占一半能量。但实测信号是实数序列，正负频率的谱是共轭对称的，只有一半信息量，所以工程上习惯只画正频率部分，并把负频率那部分的能量折算过来，这就是单边谱。

折算的方法很简单：除直流分量和奈奎斯特频率点之外，其余频点的值乘以 2。直流和奈奎斯特不能乘 2，因为它们在双边谱里没有对称的“另一半”。这个细节在写代码时特别容易漏，漏了的结果是 PSD 整体偏低一半，Grms 偏低到原来的 0.707 倍，看起来好像只是“差一点”，实际上错了 30%，做量级判定时会得出完全相反的结论。

现成的库函数默认输出单边谱并且已经处理好了这个系数，但如果你拿到的数据是从别的系统导出的、或者是别人手写的脚本算的，第一件事就是拿一段已知均方值的白噪声去验一下，看积分结果对不对得上，这比看代码快得多。

3. 从一段加速度时程到一条能交付的 PSD 曲线

3.1 预处理：不做这几步，后面全是无用功

原始加速度数据拿到手，别急着算谱，先把三件事做掉。

第一是去趋势和去均值。传感器零点漂移、电荷放大器积分漂移，都会在数据里留下一个缓变的基线，这条基线在低频段会形成巨大的能量，把真正关心的低频成分完全压住。处理方法很简单，减去最小二乘拟合的一次或二次趋势项，再减去均值。第二是剔除异常段。设备启停瞬间的冲击、传感器被碰一下、数据采集卡丢包，都会在时程里留下远高于背景的尖峰，这些尖峰在做平方运算时会被放大得极其夸张，一个点就能污染整条谱。工程做法是设定阈值，超过中位数若干倍标准差的区间直接标记并剔除。第三是确认单位。前面提过 g 和 m/s² 差 96 倍，拿到数据先看清楚标签，别猜。

注意：剔除异常段之后不要简单地把前后拼接，那会在接缝处引入不连续，反而制造出新的高频成分。正确做法是把异常区间单独砍掉，剩余的有效段分别做谱，最后按段数加权平均。

3.2 完整实现：一段可以直接照抄的代码

下面这段用 Python 的 scipy 走一遍完整流程，参数就用上一节算出来的那组。

import numpy as np from scipy import signal fs = 2000.0 # 采样率 Hz nperseg = 1024 # 段长 noverlap = 512 # 50% 重叠 # x 是已经去过趋势、去过均值、剔除过异常段的加速度序列，单位 g x = x - x.mean() f, pxx = signal.welch( x, fs=fs, window='hann', nperseg=nperseg, noverlap=noverlap, detrend='constant', scaling='density', # 输出功率谱密度而不是功率谱 return_onesided=True # 单边谱 ) # 由 PSD 积分得到 Grms grms = np.sqrt(np.trapezoid(pxx, f)) print(f"分辨率 {f[1]-f[0]:.3f} Hz, Grms = {grms:.4f} g")

几个参数值得单独说明。scaling='density'表示输出的是密度（单位 g²/Hz），如果写成'spectrum'输出的是每个频点的功率（单位 g²），后者乘上分辨率才等于前者。这两个模式在画图时看起来形状一模一样，只是纵轴差一个常数，但拿去和别人的数据对比时就会对不上号。detrend='constant'表示每段内部只去均值，如果原始数据低频漂移严重，可以改成'linear'，但要注意这会连真实的低频趋势一起去掉，做低频分析时慎用。

3.3 结果校核：三种独立手段互相印证

算完不校核，等于没算。我一般用三招交叉验证。

第一招是 Parseval 校核。直接对原始时程求标准差，再对 PSD 积分开根号，两个数应该非常接近。因为 PSD 积分理论上就等于信号均方值。如果差得超过 10%，说明某个环节出问题了，优先怀疑单位、单双边系数、窗函数修正。第二招是白噪声标定。生成一段已知方差的白噪声，方差设成 1，走一遍流程，积分结果应该是 1.0 左右，偏差在统计误差范围内。这个方法能一次性验证整条链路，比逐行读代码高效得多。第三招是看曲线的物理合理性。如果谱图在 50 Hz 出现一根特别干净的尖峰，那大概率是工频干扰或者设备本身的转频，不是随机振动成分，这时候要考虑是否需要先做陷波处理，或者至少在报告里说明。

还有一个容易被忽略的校核点：看你算出来的 PSD 峰值数量级是否合理。一般工业现场的振动，宽带随机部分在 1e-4 到 1e-2 g²/Hz 这个量级，如果是 1 g²/Hz 这种数字，八成是单位错了或者没做归一化。

3.4 反向操作：由谱形生成时程，以及 Grms 的解析积分

工程上经常需要反过来做：手上有一条试验谱，要生成对应的时域驱动信号，或者单纯想算这条谱的总量级。这两件事都不难，但都有细节。

先说解析积分。试验谱在双对数坐标上是分段直线，每段的斜率用 dB/oct 表示，换算成幂指数是 m = 斜率 / (10·log10(2)) ≈ 斜率 / 3.0103。所以 −3 dB/oct 对应 m = −1，恰好是最常见的下降段。给定一段从 f1 到 f2、起始谱值 S1，末值 S2 = S1·(f2/f1)^m，这一段的面积有解析解：

当 m ≠ −1 时，A = S1·f1/(m+1) · [(f2/f1)^(m+1) − 1]
当 m = −1 时，A = S1·f1·ln(f2/f1)

把所有段加起来开根号就是 Grms。注意 −3 dB/oct 那一段必须用第二个公式，直接用第一个会除零，程序里如果不判这个分支，会输出 inf 或者 NaN，然后整条曲线报废。这个坑我见过不止一次，包括一些商用软件的老版本里也出现过。

代码上可以这样写，用对数比来推斜率，避免手工换算出错：

def grms_from_psd(f, s): f = np.asarray(f, dtype=float) s = np.asarray(s, dtype=float) m = np.log(s[1:] / s[:-1]) / np.log(f[1:] / f[:-1]) r = f[1:] / f[:-1] area = np.where( np.isclose(m, -1.0), s[:-1] * f[:-1] * np.log(r), s[:-1] * f[:-1] / (m + 1) * (r ** (m + 1) - 1) ) return np.sqrt(area.sum())

再说时程生成。思路是在频域按 PSD 造一个复频谱，幅值按 sqrt(S·fs·n/2) 给定，相位在 0 到 2π 均匀随机，然后做逆 FFT。这样生成的时程，其统计特性会收敛到目标谱。

def psd_to_time(f, s, fs, T, seed=0): rng = np.random.default_rng(seed) n = int(fs * T) fr = np.fft.rfftfreq(n, 1.0 / fs) S = np.interp(fr, f, s, left=0.0, right=0.0) amp = np.sqrt(S * fs * n / 2.0) ph = rng.uniform(0, 2 * np.pi, fr.size) spec = amp * np.exp(1j * ph) spec[0] = 0.0 spec[-1] = np.abs(spec[-1]) # 奈奎斯特点必须为实数 return np.fft.irfft(spec, n)

生成之后一定要回头算一遍 PSD 和 Grms，和目标谱对比。实测下来，单次生成的数据因为自由度有限，Grms 会有几个百分点的偏差，多生成几组取平均就稳了。另外，生成的信号峰均比天然偏高，如果要拿去驱动振动台，还得做削峰或者按峰均比截断，不然台子会在个别尖峰上触发保护。

4. 换一个战场：PSD 在故障诊断里的二次利用

4.1 先看谱，再看包络，别一上来就做花活

设备状态监测场景里，PSD 常常被当成“看一眼就过”的中间结果，大家更爱直接上包络谱、倒频谱、谱峭度这些听起来高级的方法。但我的习惯是无论如何先把 PSD 画出来，因为它是判断“有没有问题、问题大致在哪个频带”最快的工具，很多结论根本不需要更复杂的算法。

举个例子，滚动轴承早期损伤的冲击成分很弱，被齿轮啮合、电机电磁噪声盖住，在 PSD 上往往看不出来，只能看到一圈圈宽带底噪抬升。但如果 PSD 在某个特定频带明显鼓起来，并且随着运行时间缓慢变化，那说明那个频带对应的部件状态在变，这时候再针对那个频带做包络分析，目的性和成功率都会高很多。反过来，如果不看 PSD 直接对全频带做包络，解调出来的东西基本是一锅粥。

包络谱的做法其实不复杂：先用带通滤波把包含故障冲击的共振频带选出来，再对滤波后的信号做希尔伯特变换取包络，最后对包络做谱分析。关键在于那个带通频带怎么选，选窄了怕漏掉能量，选宽了怕混进其他干扰。谱峭度可以自动帮你找那个最“冲击性强”的频带，但算出来的结果需要人判断，不是无脑信。

4.2 特征频率手算一遍，别全靠软件输参数

轴承故障的特征频率有标准公式，但参数输错导致对不上号的情况非常多，还不如手算一遍心里踏实。以常见的滚动体轴承为例，设滚动体数量 n，滚动体直径 d，节圆直径 D，接触角 α，转频 fr：

外圈通过频率 BPFO = (n/2)·fr·(1 − (d/D)·cosα)
内圈通过频率 BPFI = (n/2)·fr·(1 + (d/D)·cosα)
滚动体自转频率 BSF ≈ (D/(2d))·fr·(1 − ((d/D)·cosα)²)
保持架频率 FTF = (fr/2)·(1 − (d/D)·cosα)

算完拿去做比对时有几个经验值可以参考：外圈和保持架的频率比较稳定，就摆在那儿；内圈频率因为随轴一起转，会受载荷区调制，包络谱上往往在 BPFI 两侧出现间隔等于转频的边带；滚动体故障的频率因为滚动体自转还公转，往往表现为 BSF 及其两倍频。这些细节不写在标准里，但实际看图时非常有用。

提示：包络谱上的特征频率对转速非常敏感。转速波动 1%，BPFO 就跟着漂 1%，如果分辨率不够，谱峰会糊成一片。所以做包络分析之前，一定要确认转速信息可靠，或者用键相/转速脉冲做角度域重采样。

4.3 三个反复踩到的坑

第一个坑是“传递路径污染”。传感器装在外壳上，轴承的冲击要经过壳体、法兰、螺栓层层传递，高频部分衰减很快，你看到的频谱形状里混着结构共振，不是纯粹的故障特征。解决办法一是尽量靠近测点，二是做传递函数归一化，三是别死磕某一个频带，多看几个方向。第二个坑是“只看谱不看趋势”。单次 PSD 里某根谱线高，可能就是这一批次的装配差异，不代表要坏。真正有意义的是同一测点、同一工况下随时间的演变，所以数据管理比算法重要。第三个坑是“把转速倍频当故障”。很多现场出现的新谱峰，追下去其实是转速变化导致某阶次扫过去了，或者负载变了导致啮合频率偏移，跟损伤毫无关系。遇到异常先排除工况变化，再怀疑设备，这条能省掉大量无效排查。

5. 常见问题与排查技巧实录

5.1 一张速查表，覆盖大部分翻车现场

现象	最可能的原因	排查动作
PSD 整体偏低约一半	单边谱没乘 2，或窗函数能量没补偿	用已知方差白噪声走一遍流程
低频段一片巨大的翘起	趋势项未去除、传感器漂移	加去趋势，检查线缆与供电
曲线毛刺极多，峰位每次算都不同	段数太少，自由度不足	加大重叠、延长采集时长
谱峰被“抹平”，看不出细节	段长太短，分辨率不够	增大 nperseg，接受方差变大
Grms 与目测量级明显不符	单位 g 与 m/s² 混用	统一换算系数 9.80665
−3 dB/oct 段积分报 inf 或 NaN	斜率 m = −1 未单独处理	加等值判断分支
转频及其倍频异常突出	不平衡、对中不良、工况变化	先查工况，再看动平衡与对中
包络谱上看不到特征频率	带通频带选错、转速不准	用谱峭度选带，核对转速