基于白噪声法和直接法（周期图法）的车辆振动功率谱密度（PSD）计算的实现方法

一、理论基础

1. 白噪声法建模原理

• 路面激励模型：根据ISO标准，路面不平度速度功率谱密度为常数（白噪声特性）

  \(G_v(ω)=\frac{G_0⋅u}{0.1}\)

  其中：

  • \(G_0\)：路面不平度系数（B级路面\(G_0=64×10^{−6}m^3\)）
  • \(u\)：车速（m/s）

• 车辆动力学模型：二自由度四分之一车体模型

  

  其中zg为路面激励。

2. 直接法（周期图法）原理

通过FFT直接计算信号功率谱密度：

• \(N\)：信号长度
• \(Δf\)：频率分辨率
• \(X(f)\)：FFT变换结果

---

二、MATLAB实现代码

1. 白噪声生成与车辆振动仿真

%% 参数设置
L = 2.5;        % 轴距(m)
u = 20;         % 车速(m/s)
G0 = 64e-6;     % B级路面
m = 1500;       % 整车质量(kg)
mr = 300;       % 悬架质量(kg)
k = 200000;     % 悬架刚度(N/m)
c = 2000;       % 阻尼系数(Ns/m)
kt = 150000;    % 轮胎刚度(N/m)

%% 生成白噪声路面激励
t = 0:0.01:10;  % 时间向量(10秒)
N = length(t);
zg = sqrt(G0*u/0.1) * randn(1,N);  % 白噪声生成

%% 建立振动微分方程
M = [m, mr; mr, mr];
C = [c+c, -c; -c, c];
K = [k+kt, -kt; -kt, kt];

%% 数值求解振动响应
dt = 0.01;
[y, yr] = ode45(@(t,z) [0,1; -K/M, -C/M]*z, [0,10], [0;0;0;0], [], zg);

%% 提取垂直加速度
z_accel = y(:,1);
yr_accel = yr(:,1);

2. 直接法功率谱密度计算

%% 直接法PSD计算
Fs = 1/dt;      % 采样频率
Nfft = 2^nextpow2(length(z_accel));  % 补零到2的幂次
window = hanning(Nfft);  % 汉宁窗numoverlap = Nfft/2;

[Pxx, f] = pwelch(z_accel, window, numoverlap, Nfft, Fs, 'onesided');

%% 绘制结果
figure;
subplot(2,1,1);
plot(t, zg*1e3);
xlabel('时间(s)'); ylabel('路面激励(mm)');
title('白噪声路面激励时域信号');

subplot(2,1,2);
plot(f, 10*log10(Pxx));
xlabel('频率(Hz)'); ylabel('PSD (dB/Hz)');
title('垂直加速度功率谱密度');
grid on;

---

三、关键参数说明

参数	物理意义	典型值	影响因素
G0	路面不平度基准幅值	64×10⁻⁶ m³	路面等级（A-F级）
u	车速	20-120 km/h	车辆行驶状态
c	悬架阻尼系数	1500-3000 Ns/m	减振器性能
k	悬架刚度	150000-250000 N/m	弹簧刚度

---

四、结果分析

1. 时域特征 路面激励呈现随机波动特性，符合高斯白噪声分布 垂直加速度幅值与车速呈正相关
2. 频域特征 主能量集中在0-30Hz范围内（人体敏感频段） 峰值频率对应车辆固有频率（约1-5Hz）
3. 参数敏感性 阻尼比增加 → 加速度PSD幅值降低（但悬架动行程增大） 车速提高 → 高频成分能量占比增加

---

五、验证方法

1. 理论对比

   计算系统固有频率验证PSD主峰位置：

   wn = sqrt(k/m);  % 固有频率(rad/s)
   fn = wn/(2*pi);  % 固有频率(Hz)
   

2. 实验验证 使用加速度传感器实测数据与仿真结果对比，要求： 均方根误差 < 15% 主频偏差 < 2Hz

参考代码 白噪声法求车辆振动，直接法的功率谱密度 www.youwenfan.com/contentcnn/84838.html

六、扩展应用

1. 多通道分析

   同时计算俯仰角加速度PSD：

   pitch_accel = gradient(z_accel, dt);
   [P_pitch, f_pitch] = pwelch(pitch_accel, window, numoverlap, Nfft, Fs);
   

2. 时频分析

   使用小波变换分析瞬态振动事件：

   [cfs, f] = cwt(z_accel, 1:128, 'morl', 'SamplingPeriod', dt);
   contourf(t, f, abs(cfs).^2); shading interp;
   

---

七、注意事项

1. 数据预处理 去除直流分量：z_accel = z_accel - mean(z_accel); 阻塞处理：对超过奈奎斯特频率的信号进行镜像折叠
2. 计算效率优化 使用FFT重采样减少计算量 采用多核并行计算（parfor）
3. 工程应用建议 采样频率 ≥ 200Hz（满足奈奎斯特定理） 信号长度 ≥ 10倍系统响应时间