简介这份MATLAB程序包面向铁路工程、车辆动力学及轨道不平顺仿真研究人员用于生成符合美国AAR六级谱标准的轨道激励数据帮助评估轨道质量对列车运行性能、乘客舒适度与安全性的影响。压缩包仅3个文件共2.58MB包含2个M脚本AAR6.m为核心实现m_G2Xm.m负责谱计算与随机信号生成和1个FIG图形界面文件用于展示和交互查看轨道不平顺谱结果结构精简适合直接运行和二次开发。目前已有2141人学习使用。通过阅读和运行源码可以掌握傅里叶变换、滤波器设计与随机过程建模在轨道谱生成中的应用思路理解AAR六级谱的数学模型与参数设置并在此基础上拓展不同等级不平顺仿真为车辆动力学仿真提供输入激励。1. 先从美国六级谱说起搞车辆动力学仿真的人第一道绕不开的坎就是轨道不平顺激励。现场测试数据拿不到时最常用的替代方案是标准轨道谱而 AAR6 就是其中代表美国六级谱对应高等级、较平顺的干线线路条件。这篇文章把 AAR6 轨道谱从数学表达式到 MATLAB 时域生成串成一条完整链路包含单位换算、三角级数/IFFT 两种生成思路、谱密度还原验证以及一个能直接复用的小函数。适合正在做车辆-轨道耦合仿真、或第一次把标准谱转成时域不平顺样本的 MATLAB 工程师。六级谱参数干净、高频衰减平缓比一级谱更容易验证算法对错是练手和工程入门的理想起点。2. AAR6 轨道谱的数学模型与 MATLAB 里的解析式定义2.1 轨道不平顺分类美国谱的对应关系轨道不平顺按几何形态分四类高低垂向、方向横向、水平、轨距。前两类对车辆垂向和横向动力学影响最大美国谱也主要对垂向和横向给出解析公式。美国 FRA/AAR 谱的表达式统一为S(Ω) A · Ωc² / ( (Ω² Ωr²) · (Ω² Ωc²) )其中Ω是空间角频率单位rad/mΩc和Ωr是截断频率A是粗糙度常数。六级谱只是把参数取成了高等级线路对应的值。这个表达式是典型的“低通高通”组合Ωr控制低频拐点Ωc控制高频衰减起始位置两个截断频率把能量限制在中间频带物理意义清晰。2.2 AAR6 的参数表与量纲陷阱六级谱的经典参数如下参数垂向 (A_v)横向 (A_l)说明A0.03390.0339单位 cm²·rad/mΩc0.82450.8245单位 rad/mΩr0.02060.0206单位 rad/m参数本身不难记容易出错的是单位。文献里A常写成cm²·rad/m也就是 PSD 单位为cm²/(rad/m)这是厘米克秒制。MATLAB 里不平顺数据我们一般用米所以必须先统一单位1 cm² 1e-4 m²代入公式前把A乘1e-4。2.3 MATLAB 函数句柄直接定义 AAR6 谱在 MATLAB 里最直接的做法是定义一个函数句柄后续所有生成和验证都复用这个句柄% aar6_psd.m function S aar6_psd(Omega, type) % Omega: 空间角频率, rad/m % type: vert 垂向, lateral 横向 A 0.0339 * 1e-4; % cm^2*rad/m - m^3 OmegaC 0.8245; OmegaR 0.0206; if strcmpi(type, lateral) A 0.0339 * 1e-4; end S A * OmegaC^2 ./ ((Omega.^2 OmegaR^2) .* (Omega.^2 OmegaC^2)); end注意S(Ω)的单位是m²/(rad/m)约分后就是m³。很多新手在这里只换了厘米到米忘了rad/m对分母的影响最后生成的时域样本方差差好几个数量级。写完函数后先算一下Ω 1时的谱密度量级六级谱大约在1e-6 ~ 1e-5数量级才正常。这个函数句柄在后面生成和验证环节会被反复调用。3. 用 MATLAB 生成 AAR6 轨道不平顺时域样本3.1 从空间谱到时间谱的换算轨道谱描述的是空间域的统计特性但动力学仿真需要的是时域输入序列。假定列车以速度v匀速通过空间角频率Ω与时间角频率ω的关系为ω v · Ω按能量守恒作变量替换S_ω(ω) S_Ω(Ω) / v。这里S_ω是时间角频率谱单位变成m²·s/rad。速度越大同一空间波长在时域里出现得越快谱密度幅值越低这是合理的物理趋势。这个换算关系是后续一切代码的地基。3.2 三角级数叠加法生成不平顺三角级数法是工程中最常用的做法把连续谱按频率离散每根谱线赋予随机相位然后叠加余弦波。代码如下function x gen_aar6_tri(v, fs, trackLen, type) % 输入 % v: 车速, m/s % fs: 采样率, Hz % trackLen: 需要的轨道长度, m % type: vert 或 lateral % 输出 % x: 轨道不平顺时域序列, m lambda_min 0.5; % 最短波长 m lambda_max 200; % 最长波长 m F_min 1 / lambda_max; F_max 1 / lambda_min; nF 8192; % 频率离散点数 dF (F_max - F_min) / nF; F F_min (0:nF-1) * dF; % cycles/m Omega 2 * pi * F; % rad/m S_Omega aar6_psd(Omega, type); % m^3 S_w S_Omega / v; % m^2*s/rad, 时间角频率谱 dw 2 * pi * v * dF; % rad/s T_total trackLen / v; N floor(T_total * fs); t (0:N-1) / fs; x zeros(N, 1); for k 1:nF phase 2 * pi * rand; x x sqrt(2 * S_w(k) * dw) * cos(omega_w(k) * t phase); end x x - mean(x); % 去直流 end代码逻辑分四步先把连续谱按空间频率离散成nF根谱线再由ω v·Ω转到时间频率然后用sqrt(2·S_w·dw)计算每根谱线的幅值最后叠加随机相位余弦。幅值公式里的系数2来自单边谱处理S_w是双边定义只取正频率叠加余弦时要乘 2。nF取 8192 是为了让频率分辨率足够细避免生成序列出现周期感如果序列太长导致循环慢可以改为向量化累加。3.3 IFFT 法生成不平顺IFFT 法适合批量生成多条不平顺原理是在频域构造共轭对称的复谱再逆变换回时域function x gen_aar6_ifft(v, fs, T_total, type) N 2^nextpow2(T_total * fs); df fs / N; f (0:N-1) * df; omega 2 * pi * f; % 只取正频率部分构造单边随机谱 nPos floor(N / 2) - 1; S_Omega aar6_psd(2 * pi * f(2:nPos1) * v, type); S_w S_Omega / v; dw 2 * pi * df; X zeros(N, 1); X(2:nPos1) sqrt(2 * S_w * dw) .* exp(1i * 2 * pi * rand(nPos, 1)); X(N:-1:N-nPos1) conj(X(2:nPos1)); % 负频率共轭对称 x real(ifft(X)) * N; % ifft 自带 1/N这里乘 N 恢复幅值 x x - mean(x); end这段代码的关键在共轭对称部分X(k)和X(N-k2)必须共轭逆变换结果才是实数。x real(ifft(X)) * N里的N是很多新手容易漏的归一化系数去掉它幅值会小N倍频谱验证时对不上。提示三角级数法物理意义清楚、循环可读性好适合理解原理IFFT 法一次生成整条序列计算效率高适合 Monte Carlo 批量仿真。两种方法在N足够大时统计结果一致但 IFFT 的归一化更容易出错。3.4 两种方法的适用边界方法代码复杂度计算效率适用场景三角级数低中等单条不平顺生成、教学验证IFFT中高批量生成、随机样本集构造实际工程中如果只生成一条用于仿真三角级数就够了如果要做随机振动统计IFFT 更划算。两者生成结果的方差理论上都应等于∫ S_Ω dΩ可以互相校验。4. AAR6 生成结果怎么验证参数怎么调4.1 用 pwelch 还原 PSD 并与目标谱对比生成完不平顺序列第一件事不是直接用而是验证频谱是否和 AAR6 目标谱吻合。最常用的做法是 Welch 法谱估计v 160 / 3.6; % 160 km/h - m/s fs 500; x gen_aar6_tri(v, fs, 3000, vert); nfft 4096; [pxx, f] pwelch(x, hann(nfft), nfft / 2, nfft, fs, psd); % 目标谱: 注意 pwelch 输出是 per Hz, 所以要乘 2*pi Omega_target 2 * pi * f / v; S_target 2 * pi * aar6_psd(Omega_target, vert) / v; figure; loglog(f, pxx, b); hold on; loglog(f, S_target, r--, LineWidth, 1.5); legend(Welch估计, AAR6目标谱); xlabel(频率 Hz); ylabel(PSD m^2/Hz);pwelch默认返回单边功率谱密度单位是m²/Hz而我们构造的S_w是角频率谱单位是m²·s/rad两者差一个2π因子所以目标谱要乘2π。这一步不做低频段看起来会差一个常数倍。验证时重点看中频带0.1 ~ 20 Hz是否贴合高频端因为谱估计窗函数泄漏会有轻微下偏属正常现象。4.2 车速、采样率、轨道长度的匹配车速决定时间频率范围最短波长λ_min对应最高时间频率f_max v / λ_min。采样率fs至少要取2·f_max实际工程建议 4 倍以上车速 km/hλ_min 0.5 m 时 f_max Hz建议 fs Hz7240200160895003501941000轨道长度决定最低可分辨频率T trackLen / v最低时间频率约为1 / T。要生成 0.01 Hz 以下的低频成分需要几百米以上的轨道长度。仿真时如果只关心 0.2 Hz 以上轨道长度可以适当缩短。4.3 常见错误与排查手段生成结果 RMS 不对优先查A是否完成cm² - m²换算频谱形状在高频段翘起多半是随机相位构造时共轭对称做错低频段能量异常高检查是否忘记去直流分量。还有一个隐蔽问题nF取得太少时三角级数法的序列会出现明显的周期性重复把nF提到 4096 以上基本能消除。排查时先用rms(x)和方差积分sqrt(sum(S_w * dw))对比差 5% 以内说明能量对得上再去看频谱形状。5. 落一个能用的轨道谱生成函数5.1 函数封装与参数校验把前面代码封装成一个接口清晰、带参数校验的函数后续项目直接调用function x genAAR6Track(vKmh, fs, trackLen, type, method) % genAAR6Track 生成 AAR6 六级谱轨道不平顺时域序列 % x genAAR6Track(vKmh, fs, trackLen, type, method) % vKmh: 车速 km/h % fs: 采样率 Hz % trackLen:轨道长度 m % type: vert 垂向 / lateral 横向默认 vert % method: tri 三角级数 / ifft 逆傅里叶默认 tri validateattributes(vKmh, {numeric}, {scalar, positive}, 1); validateattributes(fs, {numeric}, {scalar, positive}, 2); v vKmh / 3.6; switch nargin case 4 method tri; case 5 % user specified otherwise error(需要 4 或 5 个输入参数); end if strcmpi(method, tri) x gen_aar6_tri(v, fs, trackLen, type); else x gen_aar6_ifft(v, fs, trackLen / v, type); end end这个函数同时保留了三角级数和 IFFT 两种内核method参数切换生成策略。参数校验放在入口处避免后期仿真脚本传错单位造成静默错误。5.2 自动校验脚本每次生成后自动跑一遍频谱对比把目标频带内的平均偏差打到命令行x genAAR6Track(160, 500, 3000, vert, tri); nfft 4096; [pxx, f] pwelch(x, hann(nfft), nfft/2, nfft, fs, psd); valid f 0.1 f 20; S_target 2 * pi * aar6_psd(2 * pi * f / v, vert) / v; err mean(abs(log10(pxx(valid)) - log10(S_target(valid)))); fprintf(RMS %.4f mm, 平均对数谱偏差 %.2f dB\n, ... rms(x) * 1000, 10 * err); x genAAR6Track(160, 500, 3000, vert, ifft);这个校验脚本的价值在于改速度、采样率、波长范围后只用一条命令就能确认生成质量。平均对数谱偏差小于 1 dB 说明生成正确超过 2 dB 就回去检查频率范围或归一化系数。把这个函数和校验脚本放到同一个目录之后做车辆动力学仿真时输入不平顺就成了两条命令的事。本文还有配套的精品资源点击获取
AAR6轨道谱的MATLAB时域生成：从单位换算到PSD验证
简介这份MATLAB程序包面向铁路工程、车辆动力学及轨道不平顺仿真研究人员用于生成符合美国AAR六级谱标准的轨道激励数据帮助评估轨道质量对列车运行性能、乘客舒适度与安全性的影响。压缩包仅3个文件共2.58MB包含2个M脚本AAR6.m为核心实现m_G2Xm.m负责谱计算与随机信号生成和1个FIG图形界面文件用于展示和交互查看轨道不平顺谱结果结构精简适合直接运行和二次开发。目前已有2141人学习使用。通过阅读和运行源码可以掌握傅里叶变换、滤波器设计与随机过程建模在轨道谱生成中的应用思路理解AAR六级谱的数学模型与参数设置并在此基础上拓展不同等级不平顺仿真为车辆动力学仿真提供输入激励。1. 先从美国六级谱说起搞车辆动力学仿真的人第一道绕不开的坎就是轨道不平顺激励。现场测试数据拿不到时最常用的替代方案是标准轨道谱而 AAR6 就是其中代表美国六级谱对应高等级、较平顺的干线线路条件。这篇文章把 AAR6 轨道谱从数学表达式到 MATLAB 时域生成串成一条完整链路包含单位换算、三角级数/IFFT 两种生成思路、谱密度还原验证以及一个能直接复用的小函数。适合正在做车辆-轨道耦合仿真、或第一次把标准谱转成时域不平顺样本的 MATLAB 工程师。六级谱参数干净、高频衰减平缓比一级谱更容易验证算法对错是练手和工程入门的理想起点。2. AAR6 轨道谱的数学模型与 MATLAB 里的解析式定义2.1 轨道不平顺分类美国谱的对应关系轨道不平顺按几何形态分四类高低垂向、方向横向、水平、轨距。前两类对车辆垂向和横向动力学影响最大美国谱也主要对垂向和横向给出解析公式。美国 FRA/AAR 谱的表达式统一为S(Ω) A · Ωc² / ( (Ω² Ωr²) · (Ω² Ωc²) )其中Ω是空间角频率单位rad/mΩc和Ωr是截断频率A是粗糙度常数。六级谱只是把参数取成了高等级线路对应的值。这个表达式是典型的“低通高通”组合Ωr控制低频拐点Ωc控制高频衰减起始位置两个截断频率把能量限制在中间频带物理意义清晰。2.2 AAR6 的参数表与量纲陷阱六级谱的经典参数如下参数垂向 (A_v)横向 (A_l)说明A0.03390.0339单位 cm²·rad/mΩc0.82450.8245单位 rad/mΩr0.02060.0206单位 rad/m参数本身不难记容易出错的是单位。文献里A常写成cm²·rad/m也就是 PSD 单位为cm²/(rad/m)这是厘米克秒制。MATLAB 里不平顺数据我们一般用米所以必须先统一单位1 cm² 1e-4 m²代入公式前把A乘1e-4。2.3 MATLAB 函数句柄直接定义 AAR6 谱在 MATLAB 里最直接的做法是定义一个函数句柄后续所有生成和验证都复用这个句柄% aar6_psd.m function S aar6_psd(Omega, type) % Omega: 空间角频率, rad/m % type: vert 垂向, lateral 横向 A 0.0339 * 1e-4; % cm^2*rad/m - m^3 OmegaC 0.8245; OmegaR 0.0206; if strcmpi(type, lateral) A 0.0339 * 1e-4; end S A * OmegaC^2 ./ ((Omega.^2 OmegaR^2) .* (Omega.^2 OmegaC^2)); end注意S(Ω)的单位是m²/(rad/m)约分后就是m³。很多新手在这里只换了厘米到米忘了rad/m对分母的影响最后生成的时域样本方差差好几个数量级。写完函数后先算一下Ω 1时的谱密度量级六级谱大约在1e-6 ~ 1e-5数量级才正常。这个函数句柄在后面生成和验证环节会被反复调用。3. 用 MATLAB 生成 AAR6 轨道不平顺时域样本3.1 从空间谱到时间谱的换算轨道谱描述的是空间域的统计特性但动力学仿真需要的是时域输入序列。假定列车以速度v匀速通过空间角频率Ω与时间角频率ω的关系为ω v · Ω按能量守恒作变量替换S_ω(ω) S_Ω(Ω) / v。这里S_ω是时间角频率谱单位变成m²·s/rad。速度越大同一空间波长在时域里出现得越快谱密度幅值越低这是合理的物理趋势。这个换算关系是后续一切代码的地基。3.2 三角级数叠加法生成不平顺三角级数法是工程中最常用的做法把连续谱按频率离散每根谱线赋予随机相位然后叠加余弦波。代码如下function x gen_aar6_tri(v, fs, trackLen, type) % 输入 % v: 车速, m/s % fs: 采样率, Hz % trackLen: 需要的轨道长度, m % type: vert 或 lateral % 输出 % x: 轨道不平顺时域序列, m lambda_min 0.5; % 最短波长 m lambda_max 200; % 最长波长 m F_min 1 / lambda_max; F_max 1 / lambda_min; nF 8192; % 频率离散点数 dF (F_max - F_min) / nF; F F_min (0:nF-1) * dF; % cycles/m Omega 2 * pi * F; % rad/m S_Omega aar6_psd(Omega, type); % m^3 S_w S_Omega / v; % m^2*s/rad, 时间角频率谱 dw 2 * pi * v * dF; % rad/s T_total trackLen / v; N floor(T_total * fs); t (0:N-1) / fs; x zeros(N, 1); for k 1:nF phase 2 * pi * rand; x x sqrt(2 * S_w(k) * dw) * cos(omega_w(k) * t phase); end x x - mean(x); % 去直流 end代码逻辑分四步先把连续谱按空间频率离散成nF根谱线再由ω v·Ω转到时间频率然后用sqrt(2·S_w·dw)计算每根谱线的幅值最后叠加随机相位余弦。幅值公式里的系数2来自单边谱处理S_w是双边定义只取正频率叠加余弦时要乘 2。nF取 8192 是为了让频率分辨率足够细避免生成序列出现周期感如果序列太长导致循环慢可以改为向量化累加。3.3 IFFT 法生成不平顺IFFT 法适合批量生成多条不平顺原理是在频域构造共轭对称的复谱再逆变换回时域function x gen_aar6_ifft(v, fs, T_total, type) N 2^nextpow2(T_total * fs); df fs / N; f (0:N-1) * df; omega 2 * pi * f; % 只取正频率部分构造单边随机谱 nPos floor(N / 2) - 1; S_Omega aar6_psd(2 * pi * f(2:nPos1) * v, type); S_w S_Omega / v; dw 2 * pi * df; X zeros(N, 1); X(2:nPos1) sqrt(2 * S_w * dw) .* exp(1i * 2 * pi * rand(nPos, 1)); X(N:-1:N-nPos1) conj(X(2:nPos1)); % 负频率共轭对称 x real(ifft(X)) * N; % ifft 自带 1/N这里乘 N 恢复幅值 x x - mean(x); end这段代码的关键在共轭对称部分X(k)和X(N-k2)必须共轭逆变换结果才是实数。x real(ifft(X)) * N里的N是很多新手容易漏的归一化系数去掉它幅值会小N倍频谱验证时对不上。提示三角级数法物理意义清楚、循环可读性好适合理解原理IFFT 法一次生成整条序列计算效率高适合 Monte Carlo 批量仿真。两种方法在N足够大时统计结果一致但 IFFT 的归一化更容易出错。3.4 两种方法的适用边界方法代码复杂度计算效率适用场景三角级数低中等单条不平顺生成、教学验证IFFT中高批量生成、随机样本集构造实际工程中如果只生成一条用于仿真三角级数就够了如果要做随机振动统计IFFT 更划算。两者生成结果的方差理论上都应等于∫ S_Ω dΩ可以互相校验。4. AAR6 生成结果怎么验证参数怎么调4.1 用 pwelch 还原 PSD 并与目标谱对比生成完不平顺序列第一件事不是直接用而是验证频谱是否和 AAR6 目标谱吻合。最常用的做法是 Welch 法谱估计v 160 / 3.6; % 160 km/h - m/s fs 500; x gen_aar6_tri(v, fs, 3000, vert); nfft 4096; [pxx, f] pwelch(x, hann(nfft), nfft / 2, nfft, fs, psd); % 目标谱: 注意 pwelch 输出是 per Hz, 所以要乘 2*pi Omega_target 2 * pi * f / v; S_target 2 * pi * aar6_psd(Omega_target, vert) / v; figure; loglog(f, pxx, b); hold on; loglog(f, S_target, r--, LineWidth, 1.5); legend(Welch估计, AAR6目标谱); xlabel(频率 Hz); ylabel(PSD m^2/Hz);pwelch默认返回单边功率谱密度单位是m²/Hz而我们构造的S_w是角频率谱单位是m²·s/rad两者差一个2π因子所以目标谱要乘2π。这一步不做低频段看起来会差一个常数倍。验证时重点看中频带0.1 ~ 20 Hz是否贴合高频端因为谱估计窗函数泄漏会有轻微下偏属正常现象。4.2 车速、采样率、轨道长度的匹配车速决定时间频率范围最短波长λ_min对应最高时间频率f_max v / λ_min。采样率fs至少要取2·f_max实际工程建议 4 倍以上车速 km/hλ_min 0.5 m 时 f_max Hz建议 fs Hz7240200160895003501941000轨道长度决定最低可分辨频率T trackLen / v最低时间频率约为1 / T。要生成 0.01 Hz 以下的低频成分需要几百米以上的轨道长度。仿真时如果只关心 0.2 Hz 以上轨道长度可以适当缩短。4.3 常见错误与排查手段生成结果 RMS 不对优先查A是否完成cm² - m²换算频谱形状在高频段翘起多半是随机相位构造时共轭对称做错低频段能量异常高检查是否忘记去直流分量。还有一个隐蔽问题nF取得太少时三角级数法的序列会出现明显的周期性重复把nF提到 4096 以上基本能消除。排查时先用rms(x)和方差积分sqrt(sum(S_w * dw))对比差 5% 以内说明能量对得上再去看频谱形状。5. 落一个能用的轨道谱生成函数5.1 函数封装与参数校验把前面代码封装成一个接口清晰、带参数校验的函数后续项目直接调用function x genAAR6Track(vKmh, fs, trackLen, type, method) % genAAR6Track 生成 AAR6 六级谱轨道不平顺时域序列 % x genAAR6Track(vKmh, fs, trackLen, type, method) % vKmh: 车速 km/h % fs: 采样率 Hz % trackLen:轨道长度 m % type: vert 垂向 / lateral 横向默认 vert % method: tri 三角级数 / ifft 逆傅里叶默认 tri validateattributes(vKmh, {numeric}, {scalar, positive}, 1); validateattributes(fs, {numeric}, {scalar, positive}, 2); v vKmh / 3.6; switch nargin case 4 method tri; case 5 % user specified otherwise error(需要 4 或 5 个输入参数); end if strcmpi(method, tri) x gen_aar6_tri(v, fs, trackLen, type); else x gen_aar6_ifft(v, fs, trackLen / v, type); end end这个函数同时保留了三角级数和 IFFT 两种内核method参数切换生成策略。参数校验放在入口处避免后期仿真脚本传错单位造成静默错误。5.2 自动校验脚本每次生成后自动跑一遍频谱对比把目标频带内的平均偏差打到命令行x genAAR6Track(160, 500, 3000, vert, tri); nfft 4096; [pxx, f] pwelch(x, hann(nfft), nfft/2, nfft, fs, psd); valid f 0.1 f 20; S_target 2 * pi * aar6_psd(2 * pi * f / v, vert) / v; err mean(abs(log10(pxx(valid)) - log10(S_target(valid)))); fprintf(RMS %.4f mm, 平均对数谱偏差 %.2f dB\n, ... rms(x) * 1000, 10 * err); x genAAR6Track(160, 500, 3000, vert, ifft);这个校验脚本的价值在于改速度、采样率、波长范围后只用一条命令就能确认生成质量。平均对数谱偏差小于 1 dB 说明生成正确超过 2 dB 就回去检查频率范围或归一化系数。把这个函数和校验脚本放到同一个目录之后做车辆动力学仿真时输入不平顺就成了两条命令的事。本文还有配套的精品资源点击获取
