AAR6轨道谱MATLAB实现：从PSD公式到时域不平顺样本生成与验证

简介：这是一份面向铁路工程、车辆动力学与 MATLAB 信号处理学习者的轨道谱资源，聚焦美国 AAR 六级谱不平顺标准与实现，以 MATLAB 脚本为核心模拟并生成符合六级谱要求的轨道激励数据，可用于课程设计、科研预研与轨道质量评估。压缩包共 3 个文件，包括两个 .m 源码文件和一个 .fig 图形界面文件，整体大小 2.58MB。其中 .m 文件分别承担 AAR 谱定义、参数设置与随机不平顺数据生成，可调整级别和频段以适配不同仿真需求；.fig 文件提供可视化交互界面，便于直接观察和检查轨道不平顺谱结果。已有 2141 人学习下载。借助这份资源，读者能掌握傅里叶变换、滤波器设计、随机信号生成等方法在轨道谱仿真中的实际应用，并可将生成的激励数据用于车辆动力学响应、乘客舒适度及运行安全性研究，是理解 AAR 六级谱不平顺、开展相关仿真的实用工具。

1. 从白噪声激励到AAR6轨道谱：为什么仿真先换激励源

做车辆动力学仿真时，最容易被低估的不是悬挂参数，而是轮轨激励源。很多人习惯用白噪声或几段正弦波代替轨道不平顺，结果垂向加速度响应在某个频率上被异常放大，排查半天发现是激励谱压根不对。AAR6_轨道谱_美国六级谱不平顺matlab，实际就是北美AAR标准里的六级轨道谱，用MATLAB把功率谱密度还原成时域不平顺样本的标准处理流程。它在轨道车辆动力学、构架疲劳、悬挂参数优化里是默认激励源之一。这篇文章写给需要自己从PSD公式生成样本的工程师，而不是只调工具箱接口的人。下面按公式、参数、MATLAB实现、FFT验证和收尾工程化的顺序讲清楚。

2. AAR6美国六级谱的PSD公式与参数选型

2.1 六级谱在AAR谱系中的位置：五级、六级、七级该用谁

AAR是美国铁路协会（Association of American Railroads）体系下定义的轨道几何不平顺谱，按照线路等级从低到高分为多个级别。六级谱是其中使用频率很高的一档，对应的是较高速度干线客货共线线路。常见对应关系如表所示。

选择六级谱做设计输入时，主要理由是它比五级谱幅值低，能反映较优轨道状态；又比七级谱更容易在工程中留出安全裕量。我一般做车辆悬挂前期选型时会直接用六级谱校核，等参数收敛后再用七级谱复核疲劳应力。这里要说明一点：轨道谱等级和车辆运行速度不是强绑定，六级谱在低于其设计速度下使用完全合理，因为激励幅值偏低，结果偏保守。

2.2 高低、轨向、水平、轨距四通道PSD公式

AAR谱采用空间角频率Ω作为自变量，Ω = 2π/λ，λ是轨道不平顺波长，单位是m。四个方向的功率谱密度常用解析式如下。

高低不平顺：

S_v(Ω) = A_v·Ω_c² / ((Ω² + Ω_r²)·(Ω² + Ω_c²))

轨向不平顺：

S_a(Ω) = A_a·Ω_c² / ((Ω² + Ω_r²)·(Ω² + Ω_c²))

水平不平顺：

S_c(Ω) = A_c·Ω_c² / ((Ω² + Ω_r²)·(Ω² + Ω_c²))

轨距不平顺：

S_g(Ω) = A_g·Ω_c² / ((Ω² + Ω_r²)·(Ω² + Ω_c²))

四个通道的差异集中在幅值系数A_v、A_a、A_c、A_g上。六级谱一组文献常用参数为：Ω_c = 0.8246，Ω_r = 0.0206，Ω_s = 0.4380，A_v = 0.0339，A_a = 0.0145，A_c = 0.0145，A_g = 0.00667。注意不同标准文本里A_v可能写为0.0339 cm²·rad/m，在MATLAB中参与计算时要换算成m²量级，也就是A_v取3.39e-6，A_a取1.45e-6，A_g取6.67e-7。如果直接抄0.0339进公式，生成序列幅值会大两个量级，这是最常见的数量级事故。

提示：具体项目里到底采用哪一组AAR参数，要以你手上购置或委托方的线路谱验收报告为准。不同文献的A值差异属于正常，不要拿本文数值对标准验收结果。

2.3 车速如何把空间谱变成时间谱

车辆以速度v运行时，轮对感受到的激励不再以空间波长为主，而是以时间频率为主。空间角频率Ω和时间圆频率ω的映射关系是ω = v·Ω。由于功率谱密度在频段内的积分代表不平顺方差，空间谱和时间谱之间存在换算系数：

S_t(ω) = S_Ω(Ω) / v

也就是说，车速越高，同样的空间不平顺被压缩到更高时间频率上，时间谱的幅值会降低。这个换算在做车辆频响分析时很关键，但生成位置域样本时并不需要；只在最后把位置序列变成时间序列、供Simulink等仿真器使用时才涉及到。

2.4 用结构体把谱参数和数据管线串起来

在MATLAB里，我习惯先把谱参数固化为一个结构体，而不是散落在脚本各处。这样后续更换五级谱、七级谱，或者调整截止波长，只需要改一处。

% AAR6谱参数与生成控制量 aar6.av = 3.39e-6; % 高低不平顺幅值系数, m^2*rad/m aar6.aa = 1.45e-6; % 轨向不平顺幅值系数, m^2*rad/m aar6.ac = 1.45e-6; % 水平不平顺幅值系数, m^2*rad/m aar6.ag = 6.67e-7; % 轨距不平顺幅值系数, m^2*rad/m aar6.omega_c = 0.8246; % 截断角频率, rad/m aar6.omega_r = 0.0206; % 低频滚降角频率, rad/m aar6.omega_s = 0.4380; % 轨距谱特征角频率, rad/m aar6.lambda_min = 2; % 最短波长, m aar6.lambda_max = 200; % 最长波长, m

参数说明：Omega_c和Omega_r决定了谱形状中高频衰减和低频抬升的转折点，通常不需要按线路等级修改；真正需要按等级换的是A_v、A_a、A_c、A_g四个幅值系数。lambda_min和lambda_max不来自标准，是合成时由用户定义的波长截断范围，建议和车辆模型关注的频率范围对应，不要盲目拉长到几百米。

3. 用MATLAB按AAR6谱合成时域不平顺的完整脚本

3.1 谐波叠加法：够用的精度和可解释性

把功率谱密度转成时域序列，工业界常用谐波叠加法，也叫三角级数法。原理是对PSD曲线按频率网格离散，在每一小段上把谱能量折算成正弦分量的幅值，再随机赋相位后叠加。公式是：

s(x) = Σ sqrt(2·S(Ω_k)·ΔΩ)·cos(Ω_k·x + φ_k)

其中φ_k是[0, 2π)均匀分布的随机相位。这个方法的优点是实现直观、代码量小、谱形误差容易控制，对一个固定PSD做几千个谐波叠加在MATLAB里几秒内完成。缺点是需要人为设置波长范围，超出截断范围的谱内容会丢失。另一种方案是频域白噪声滤波法，利用FFT做频域整形，速度更快，但边界效应和随机相位控制不如谐波叠加法直接。做动力学激励输入，谐波叠加法足够，还能逐个分量检查能量分布，因此后面代码采用这个方案。

3.2 aar6_vertical_psd_series：一个能直接跑的MATLAB函数

下面给出生成AAR6六级谱高低不平顺的位置域样本函数。代码采用空间域直接生成，以x为里程位置，输出单位为m。

function [s, x] = aar6_vertical_psd_series(L_total, dx, seed) % 生成AAR6美国六级谱高低不平顺位置域序列 % 输入: % L_total - 总里程, m % dx - 空间采样间隔, m % seed - 随机种子 % 输出: % s - 高低不平顺, m % x - 里程位置, m arguments L_total (1,1) double = 1000 dx (1,1) double = 0.1 seed (1,1) double = 20260601 end % 波长范围与空间角频率网格 lambda_min = 2; lambda_max = 200; Omega_min = 2*pi/lambda_max; Omega_max = 2*pi/lambda_min; N_psd = 2000; Omega = linspace(Omega_min, Omega_max, N_psd)'; dOmega = Omega(2) - Omega(1); % AAR6谱参数(已换算成m为单位) Av = 3.39e-6; Omega_c = 0.8246; Omega_r = 0.0206; % 高低不平顺功率谱密度 S_v = Av * Omega_c^2 ./ ((Omega.^2 + Omega_r^2) .* (Omega.^2 + Omega_c^2)); % 设置随机相位 rng(seed); phi = 2*pi*rand(N_psd, 1); % 位置向量 N = ceil(L_total / dx); x = (0:N-1)' * dx; % 谐波叠加 s = zeros(N, 1); for k = 1:N_psd s = s + sqrt(2 * S_v(k) * dOmega) * cos(Omega(k) * x + phi(k)); end end

逻辑说明：函数先按波长2m到200m生成空间角频率网格，再计算每个频点上的PSD值。sqrt(2*S_v(k)*dOmega)是该频点对应正弦分量的幅值，物理含义是该频段内不平顺方差贡献对应的等效幅值。随机相位phi保证各谐波相互独立，避免出现周期性叠加。N_psd取2000意味着在2m到200m波长范围内有2000个频率点，相邻谐波间隔足够密，生成序列不会出现明显的周期性重复感。dx取0.1m时，2m最短波长有20个采样点，满足空间采样定理，不会出现高频混叠。若把dx放宽到0.5m，最短波长必须同步提高到至少1m，否则生成序列高频段失真。

提示：seed=20260601是示例固定种子。固定种子是为了同一组参数能复现同一组不平顺序列，方便对比悬挂参数迭代前后的响应差异。

3.3 关键参数表：dx、波长范围、频域点数、随机相位

这组参数直接影响生成结果的谱形和计算开销，我把常用组合整理如下。

N_psd和dOmega的乘积决定了可模拟的最大波长范围，如果Omega_min取得过低而N_psd不够，最低频段能量会稀疏，低频PSD拟合会变差。我一般让N_psd保持在2000以上，这样pwelch验证时低频段不会出现明显锯齿。

3.4 从位置域到时间序列：插值与采样率选择

车辆仿真器里通常需要按固定时间步长读取输入。位置域生成后，用车速v把位置映射到时间，再插值即可。时间步长dt对应的等效空间步长是v·dt，因此dt必须满足v·dt < lambda_min/10，否则会丢失短波激励。比如车速176 km/h，v约48.9m/s，dt取0.001s时空间步长0.0489m，16.5倍的2m波长采样密度，结果可靠。

v = 176 / 3.6; % 车速, m/s t = (0:round(L_total/v/dt))' * dt; x_t = v * t; s_t = interp1(x, s, x_t, 'linear', 'extrap');

插值说明：interp1使用线性插值足够，因为位置域样本的采样密度远高于目标时间步长下的空间采样密度，线性插值不会引入额外高频误差。使用extrap是为了处理最后一点因取整造成的位置越界。如果后续要做FFT分析时间响应，建议用resample而不是插值，因为resample会做抗混叠滤波，但interp1此时更快且对激励影响极小。

4. 实战排错：用FFT验证AAR6生成结果

4.1 用pwelch把生成序列和理论谱对齐

合成完不能直接用，先验证生成序列的功率谱密度是否贴着理论PSD。验证脚本用pwelch做周期图估计，然后和理论曲线画在同一张图上。

Nfft = 4096; [pxx, f] = pwelch(detrend(s), hann(Nfft), Nfft/2, Nfft, 1/dx); Omega_est = 2*pi*f; % pwelch返回的是单边谱, 折到角频率域需要除以2*pi S_est = pxx / (2*pi); S_theory = Av * Omega_c^2 ./ ... ((Omega_est.^2 + Omega_r^2) .* (Omega_est.^2 + Omega_c^2)); loglog(Omega_est, S_est, '.', Omega_est, S_theory, '-'); legend('仿真序列PSD', 'AAR6理论谱', 'Location', 'best'); xlabel('空间角频率 \Omega (rad/m)'); ylabel('PSD (m^3/rad)');

逻辑说明：pwelch对位置序列做功率谱估计时，fs参数填1/dx，代表每米采样点数；得到的是按循环频率f（cycles/m）分布的单边谱。由于理论公式以Ω为自变量，而Ω=2πf，需要把f轴的PSD折算到角频率域，也就是除以2π。如果不做这步，低频段两条曲线会整体偏到理论值的2π倍附近。图里低频段出现明显离散点是正常的，因为波长200m对应的频段只有几个FFT bin，估计方差大；高频段应和理论线贴合。

4.2 轨距谱奇点和不平顺通道合成时的3个坑

轨距不平顺的公式形式比较特殊，当Ω趋近Omega_s时，分子分母会出现一个接近零的因子，导致PSD在特定波长上异常大。常见做法是把频率网格刻意避开Omega_s附近的点，或在合成前检查理论值不为Inf/NaN。我一般会在频点生成后加一句：

S_g(isinf(S_g)) = 0; S_g(isnan(S_g)) = 0;

要把Inf和NaN置0，而不是直接排除频点，否则该频段能量贡献完全丢失。除了这处奇点，实践里还有两个高频坑。

第一个坑是左右轨不平顺做完全独立随机序列。这样会丢失左右轨的相干性，水平不平顺的功率谱会异常放大。正确做法是用一组公共随机相位先生成中心线高低，再对右轨叠一个独立小扰动，左轨取公共部分减扰动，从而保留相干性。

第二个坑是没有做直流去除。谐波叠加法的每个分量cos(Ω_k·x+phi_k)在一个完整周期内均值为零，但如果x从0开始且周期没对齐，序列整体会有一个固定的直流偏置。动力学仿真中直流分量不会影响动态响应，却会在统计均值和轨距校核时造成干扰。合成后执行s = detrend(s)即可。

4.3 参数抄错引发的数量级错误排查表

把常见错误现象和对应原因列成表，适合快速定位问题。

排查时先看数量级，再看谱形。数量级错基本是单位问题，谱形错基本是波长范围和频点密度问题。这两类问题分别对应2分钟和10分钟内的修复成本。

5. 把AAR6不平顺变成可复现激励文件的3个收尾动作

5.1 用固定随机种子烘焙左轨、右轨样本

生成左轨和右轨的高低不平顺时，我推荐用公共随机数流构造，而不是分别调两次函数。做法是保存一组主相位，再为右轨叠加一个独立的次扰动，两者共用同一个seed。

rng(2026); phi_common = 2*pi*rand(N_psd, 1); phi_right_extra = 0.2 * 2*pi*rand(N_psd, 1); s_left = zeros(N, 1); s_right = zeros(N, 1); for k = 1:N_psd s_left = s_left + sqrt(2*S_v(k)*dOmega) * cos(Omega(k)*x + phi_common(k)); s_right = s_right + sqrt(2*S_v(k)*dOmega) * cos(Omega(k)*x + phi_common(k) + phi_right_extra(k)); end

参数说明：主相位phi_common决定左右轨的公共不平顺成分，phi_right_extra幅度取0.2rad，表示右轨在公共相位差上的小幅偏移。这个做法比生成两条完全独立序列更贴近真实轨道，因为左右轨共享轨枕、道床和路基的沉降趋势，只在钢轨自身偏差上有差异。固定seed后，任何同事重跑脚本都能得到同一对左右轨序列。

5.2 导出CSV/Mat进Simulink或外部求解器

生成完成后，把数据落成文件比直接留在工作区更可靠。导出CSV用writetable即可，导出给Simulink用timeseries。

T = table(x, s_left, s_right, 'VariableNames', {'x_m', 'left_mm', 'right_mm'}); writetable(T, 'AAR6_class6_irregularity.csv'); ts_left = timeseries(s_left, x / v, 'Name', 'vertical_left'); ts_right = timeseries(s_right, x / v, 'Name', 'vertical_right');

表格里把单位从m换算成mm更适合工程报表阅读，导出前乘以1000再写入。timeseries的时间轴用x/v把位置映射成时间序列，v是车辆速度，这样Simulink中直接从workspace引用变量名即可。若外部求解器只要位置序列，保留x_m列就行，注意时间序列和空间序列不要混淆。

5.3 留给下一版的自检函数

最后建议把验证写成可复用的独立函数，放进仿真工程里随激励生成一起跑。一个最小自检要覆盖三点：序列长度是否等于ceil(L_total/dx)+1，PSD估计在1m到100m波段的百分误差是否在20%以内，左右轨互相关在零延时的系数是否落在0.7到0.95之间。互相关系数过低说明左右轨独立随机成分太多，需要调小phi_right_extra；过高说明缺乏轨距随机分量，需要反向调整。这个自检函数跑一遍，整个AAR6六级谱激励文件的可靠性就有底了。

本文还有配套的精品资源，点击获取