轨道不平顺分析的MATLAB程序毕业论文及实习报告1 精品

目 录

第1章 绪论„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„1 1．1 引言„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„1 1．2 研究背景„„„„„„„„„„„„„„„„„„„„„„„„„„„„2 1．3 本文主要研究工作„„„„„„„„„„„„„„„„„„„„„„„„4 第2章 轨道随机不平顺„„„„„„„„„„„„„„„„„„„„„„„„„5 2．1 轨道不平顺及其形成原因„„„„„„„„„„„„„„„„„„„„„5 2．2 轨道不平顺的分类„„„„„„„„„„„„„„„„„„„„„„„„6 第3章 轨道不平顺功率谱分析„„„„„„„„„„„„„„„„„„„„„„9 3．1 轨道不平顺状态的评估方法„„„„„„„„„„„„„„„„„„„„9 3．1．1 局部不平顺幅值超限评分法„„„„„„„„„„„„„„„„„9 3．1．2 轨道质量指数（简称TQI）评价方法„„„„„„„„„„„„„9 3．1．3 局部不平顺幅值超限评分法与轨道质量指数评价法的比较 „„„10 3．2 轨道谱研究概述„„„„„„„„„„„„„„„„„„„„„„„„ 11 3．2．1 国外铁路轨道谱的研究情况 „„„„„„„„„„„„„„„„11 3．2．2 国内铁路轨道谱的研究情况 „„„„„„„„„„„„„„„„16 3．3 国内外轨道谱比较分析 „„„„„„„„„„„„„„„„„„„„„21 3．3．1 普通线路轨道谱的比较 „„„„„„„„„„„„„„„„„„22 3．3．2 高速线路轨道谱的比较 „„„„„„„„„„„„„„„„„„24 3．3．3 结论 „„„„„„„„„„„„„„„„„„„„„„„„„„27 第4章 轨道谱估计 „„„„„„„„„„„„„„„„„„„„„„„„„„28 4．1 随机过程及其特征描述 „„„„„„„„„„„„„„„„„„„„„28 4．1．1 随机过程 „„„„„„„„„„„„„„„„„„„„„„„„28 4．1．2 平稳随机过程 „„„„„„„„„„„„„„„„„„„„„„28 4．1．3 随机信号的相关函数 „„„„„„„„„„„„„„„„„„„29 4．1．4 随机信号的功率谱 „„„„„„„„„„„„„„„„„„„„30

4．2 功率谱估计的各种方法 „„„„„„„„„„„„„„„„„„„„„31 4．2．1 古典谱估计 „„„„„„„„„„„„„„„„„„„„„„„31 4．2．2 最大熵谱估计„„„„„„„„„„„„„„„„„„„„„„36 4．2．3 谱估计方法的比较分析„„„„„„„„„„„„„„„„„„37 4．3 实测轨道谱与现有轨道谱的比较„„„„„„„„„„„„„„„„„40 4．4 线路平顺性趋势分析„„„„„„„„„„„„„„„„„„„„„„41 第5章 轨道不平顺数值模拟„„„„„„„„„„„„„„„„„„„„„„43 5．1 国内外常用的数值模拟方法„„„„„„„„„„„„„„„„„„„43 5．1．1 白噪声滤波法及二次滤波法„„„„„„„„„„„„„„„„43 5．1．2 三角级数法„„„„„„„„„„„„„„„„„„„„„„„44 5．2 逆傅氏变换法„„„„„„„„„„„„„„„„„„„„„„„„„45 5．2．1 估计功率谱的Blackman-Turkey（BT）法 „„„„„„„„„„46 5．2．2 逆傅氏变换法的计算步骤„„„„„„„„„„„„„„„„„46 5．2．3 算例„„„„„„„„„„„„„„„„„„„„„„„„„„48 结论与展望 „„„„„„„„„„„„„„„„„„„„„„„„„„„„„50 致谢 „„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„52 参考文献 „„„„„„„„„„„„„„„„„„„„„„„„„„„„„„53 附录 „„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„55

第1章 绪论

1.1 引言

铁路轨道是一种特殊的结构物，它大多支承在密实度和弹性都很不均匀的道床和路基上，其工作条件十分复杂。在线路建设和机车运行过程中，钢轨将不可避免的产生左右、高低、方向、轨距等不平顺，不仅对列车运行的稳定性和舒适度产生不良影响，同时作为机车车辆/轨道系统的激扰源，将引起轮轨接触的动作用力，对设备造成破坏。

国内外铁路运输系统的实践证明，即使轨道结构在强度方面完全满足要求，而当轨道的平顺性不良时，由轨道不平顺所引起的车辆振动和轮轨动作用力将随着行车速度的提高而成倍的急剧增大，增加了列车脱轨的危险。反之，若轨道的平顺性满足要求，列车的振动和轮轨动作用力不大，行车的安全和平稳舒适度就能够得到保证，轨道和车辆部件的寿命、维修周期也较长。轨道的平顺性是线路方面制约行车速度的主要因素，是铁路管理的核心问题。

我国目前铁路技术发展的目标是逐步实现客运高速、货运重载和行车高密度。铁路线路设备作为重要的基础设施，将面临快速和重载的双重压力。研究轨道不平顺，对于机车车辆、线路的设计，车/轨系统动力学研究以及轨道状态的科学评定都有重要意义。

轨道不平顺的产生和发展是很多因素共同作用的结果。受载荷的随机性、路基的不均匀沉降、养护水平不同等因素影响。轨检车是测量轨道几何形位的检测设备，用于检测轨道高低、轨向、水平、轨距等不平顺以及车体、轴箱的垂向及横向振动加速度、里程、速度、钢轨断面等参数，采用模拟数字混合处理方法，对检测信号进行预处理、解偏、修正补偿、超限值摘取、统计、评价、显示及存储。

经轨检车检测并进行处理后所得到的轨道不平顺数据是里程的随机函数，应当作为随机过程来研究。随机过程不能用确定的数学函数关系来描述，只能考察其统计分布，用集平均的方式加以描述。在工程实际中，通常使用统计均值、自相关函数、功率谱密度函数等来描述轨道不平顺的统计特征和频域结构。本文中着重采用轨道功率谱密度来描述轨道不平顺的幅值相对于空间频率的分布特征，这样得到函数图形就称为轨道谱图。

如上所述，轨道不平顺是一个随机过程，它可从实测数据选取某段采用，也可

以根据已知功率谱密度函数的拟合公式进行数值模拟，这种从频域向时域转换的方法在求解车辆-轨道系统随机振动响应时十分有用。目前国内外常用的轨道不平顺数值模拟方法主要有二次滤波法、三角级数法、白噪声法以及逆傅氏变换法。本文在数值模拟方面也作了不少工作，利用文献[1]中提出的逆傅氏变换法对美国六级线路谱做了数值模拟，再通过模拟谱线与解析谱线的比较，再次验证了该方法的正确性和可靠性。

1.2 研究背景

本文对轨道不平顺采样信号进行的谱分析和数值分析，均将轨道不平顺信号视作均匀平稳的随机信号，而且可以证明这一假设是合理的，所得到的分析结果精度是足够的。

轨道不平顺的统计特征只能依靠线路实地测量获得。英国铁路于1964年就开始了这项测试工作，是世界上开展这一研究最早的国家之一。目前，英国、日本、德国、美国、俄罗斯、印度、捷克等国家都测定了各自的轨道不平顺的谱密度和相关函数。我国也在这方面做了不少研究工作。1982年铁道科学研究院罗林等讨论了各种轨道不平顺的测量方法，用“惯性基准法”测量了轨道随机不平顺，并对大量轨检车实测的不平顺数据进行了分析处理，列举了平稳轨道不平顺的样本记录功率谱密度。1985年原长沙铁道学院随机振动研究室将轨道不平顺分为弹性和几何不平顺，对先后三次用地面测量方法在京广线测定的轨道不平顺进行分析处理得到了各种不平顺谱，并且统计了我国Ⅰ级干线轨道不平顺功率谱密度的解析表达式。但是应该认识到，两单位早期研究中所获得的轨道谱分辨率精度都不高，尤其是样本数据太少（当时长沙铁道学院测取的数据仅数百米，铁道科学研究院测取的数据也只有数十公里），所以都不足以代表我国铁路轨道不平顺的统计特征。有鉴于此，20世纪90年代末，铁道科学研究院对我国轨道不平顺进行了深入细致的研究，在我国东南西北各主要干线约四万公里轨检车检测数据和部分地面测量数据的基础上，经筛选、分类处理、统计分析，提出了我国主要干线高低、水平、轨向三种轨道不平顺和部分轨道长波长不平顺的功率谱密度，其中包括重载线、提速线、准高速线、高速试验线、不同轨道结构以及特大桥梁等各种情况下的轨道不平顺功率谱密度。

[2]

国际上对采样信号进行谱分析的研究始于18世纪，英国科学家牛顿最早给出了“谱”的概念。后来，1822年，法国工程师傅立叶提出了著名的傅立叶谐波分析

理论。该理论至今依然是进行信号分析和信号处理的理论基础。傅立叶级数提出后，首先在人们观测自然界中的周期现象时得到应用。19世纪末，Schuster提出用傅立叶系数的幅度平方作为函数中功率的度量，并将其命名为“周期图”（periodogram）。这是经典谱估计的最早提法，这种提法至今仍被沿用，只不过现在是用快速傅立叶变换（FFT）来计算离散傅立叶变换（DFT），用DFT的幅度平方作为信号中功率的度量。周期图法较差的方差性能促使人们研究另外的分析方法。1927年，Yule提出用线性回归方程来模拟一个时间序列，Yule的工作实际上成了现代谱估计中最重要的方法——参数模型法谱估计的基础。Walker利用Yule的分析方法研究了衰减正弦时间序列，得出了Yule-Walker方程，可以说Yule和Walker都是开拓自回归模型的先锋。1930年，著名控制理论专家Wiener在他的著作中首次精确定义了一个随机过程的自相关函数及功率谱密度，并把谱分析建立在随机过程统计特征的基础上，即“功率谱密度是随机过程二阶统计量自相关函数的傅立叶变换”，这就是Wiener-Khintchine定理，该定理把功率谱密度定义为频率的连续函数，而不再像以前定义为离散的谐波频率的函数。1949年，Turkey 根据Wiener-Khintchine定理提出了对有限长数据进行谱估计的自相关法。即利用有限长数据估计自相关函数，再对该自相关函数求傅立叶变换，从而得到谱的估计。1958年，Blackman（二阶升余弦窗的提出者）和Turkey在出版的有关经典谱估计的专著中讨论了自相关谱估计方法，所以经典谱估计的自相关法又叫BT（Blackman-Turkey）法。周期图法和自相关法是经典谱估计的两个基本方法。1948年，Bartlett首次提出了用自回归模型系数计算功率谱，自回归模型和线性预测都用到了1911年提出的Toeplitz矩阵结构，Levinson曾根据该矩阵的特点于1947年提出了解Yule-Walker方程的快速解法，这些都为现代谱估计的发展打下了良好的理论基础。1965年，Cooley和Tuekey提出的FFT算法，也促进了谱估计的迅速发展。现代谱估计的提出主要是针对经典谱估计（周期图法和自相关法）分辨率低和方差性能不好的问题。1967年，Burg提出的最大熵谱估计，就是朝着高分辨率谱估计所作的最有意义的努力。

[3]

对于轨道随机不平顺的数值模拟，国内外学者也作了大量的研究，现在国内外常用的数值模拟方法主要有二次滤波法、三角级数法、白噪声滤波法和逆傅氏变换法

[1，4~6]

，其中，逆傅氏变换法通用性强，数据处理速度快，精度较高，是一种简单

实用的方法。其它几种方法在不同程度上都存在问题，要么通用性差，要么处理速度慢，要么精度差，因此本文中运用逆傅氏变换法进行数值模拟，以便快速准确的得到结果。

1.3 本文主要研究工作

1、简要分析了轨道不平顺的形成原因和影响因素，介绍了轨道不平顺的分类。在不失一般性的前提下，将轨道不平顺视为平稳的各态历经随机过程，明确了研究对象的数学模型是随机过程均匀平稳采样信号。

2、分析比较了几种评价轨道状况的方法，提出建立轨道的功率密度谱才能最有效的反映轨道状况。针对我国三大干线和郑武线高速试验段轨道谱、美国五、六级线路谱和德国高、低干扰轨道谱，利用MATLAB绘制各种轨道谱的谱线并分析比较它们之间的差异，得到了我国轨道谱的优劣与适用条件。

3、根据已有的轨道不平顺实际检测数据（京广线K80~K90区间下行线左钢轨高低不平顺检测数据（已作预处理）），运用多种功率谱估计方法，利用MATLAB编程计算，得到相应的轨道不平顺功率谱密度曲线，即轨道谱图,并对这几种谱估计方法进行分析比较，得到各种谱估计方法的优劣。

4、针对上述轨道谱图，综合线路实际条件，分析了该路段轨道不平顺的发展变化趋势，提出了一些具有针对性的养护维修意见。

5、基于美国六级线路谱，利用MATLAB编程实现轨道随机不平顺时域样本的数值模拟（逆傅氏变换法），得到轨道不平顺时域的随机样本，再将时域样本反衍，与功率谱的解析解比较，验证该方法的正确性和可靠性。

第2章 轨道随机不平顺

2.1 轨道不平顺及其形成原因

铁路轨道由于列车车轮的作用，轨面会产生不均匀磨耗，轨头会被磨伤，在无缝线路时存在焊接接头，在列车和环境温度载荷作用下，轨道的几何形状会发生恶化，使轨道不再处于平顺状态。在线路的平直区段，钢轨并不是呈理想的平直状态，两根钢轨在高低和左右方向上相对于理想的平直轨道呈某种波状变化而产生偏差，这种几何参数的偏差就称为轨道的不平顺。轨道在没有车轮载荷作用时所呈现的不平顺称为静态不平顺；车辆沿轨道运行时，轨道在车轮载荷作用下沿长度方向每点呈现不均匀的弹性下沉，由此形成的不平顺称为动态不平顺。

图2-1 轨道不平顺的状态

轨道不平顺是轮轨系统的激扰源，是引起机车车辆产生振动和轮轨动作用力的主要原因，对行车安全、平稳、舒适性、车辆和轨道部件的寿命以及环境噪声等都有重要影响。在快速高速、重载行车条件下，轨道不平顺对行车的影响更大，是轨道方面直接限制提高行车速度的主要因素。轨道平顺性也是轨道结构综合性能和承载能力的重要体现。

轨道不平顺还是机车车辆动力学设计，确定悬挂减震参数不可缺少的主要输入函数，为了研究轮轨相互作用，研究车辆、轨道动力学性能，进行动力学试验、计算机仿真，都需要深入研究轨道不平顺。

轨道不平顺的形成和发展是诸多具有随机性的因素共同作用的结果，这些因素包括：钢轨的初始平直性，钢轨磨耗、损伤，钢轨间距不均、质量不一，线路施工高程偏差，道床的级配和强度不均、松动、脏污、板结、路基下沉不均匀、刚度变

化，道床、路基的不均匀残余变形积累，机车车辆时刻变化的动力作用，以及雨雪、气温、地震等自然环境因素，它们综合作用，造成了轨道不平顺的随机特性。实际运营中的轨道不平顺都是经常变化的，显得很不规则。通常不同位置的轨道不平顺幅值和波长都各不相同。所以轨道不平顺波形不能用单一的简谐、三角、指数或抛物线等规则的波形来描述，可以看作是由许多无法预知的不同频率、不同幅值、不同相位的简谐波叠加而成的复杂的随机波。从本质上讲，轨道不平顺是一个随机过程，是里程位置的随机函数，任一特定区段的轨道不平顺可看成随机过程的一个样本。轨道不平顺的随机性特征决定了轨道不平顺的描述不能用一个明确的数学表达式来表示，而只能用随机振动理论中描述随机数据的“均方差”、“方差”和“功率谱密度函数”等统计函数来表达轨道不平顺的特征，从时域、频域对轨道不平顺的幅值特征、波长结构以及是否包含周期性波形等作全面的描述。此外，对于轨道不平顺的局部特征，可以用幅值、半波长、四分之一波长变化率以及平均变化率等参 量来表述。

[7]

实测轨道不平顺105高低不平顺幅值/mm0-5-10-15-20010002000300040005000600070008000900010000里程/m 图2-2 轨道不平顺实测波形（高低不平顺）

2.2 轨道不平顺分类

轨道几何不平顺主要可以分为以下五种：

[8]

1、轨道高低不平顺

轨道高低变化是垂直方向的不平顺，是指钢轨表面在同一轮载作用下形成的沿长度方向的高低不平，是由于轨面不均匀的磨耗、低接头、弹性垫层和轨枕、道床、

路基的弹性不均、各扣件和部件间的扣紧程度和间隙不等、轨枕底部有暗坑、道床和路基的永久变形等原因所造成的。轨道高低不平顺如图2-3中zr所示。

图2-3 轨道高低不平顺

轨道高低不平顺是引起机车车辆竖向振动，使轮轨间产生巨大惯性力的主要因素。研究表明：波长大的不平顺对车体振动的影响较大，幅值及幅值相应的平均变化率较大也是大的不平顺，会引起车体强烈振动。波长较短或变化率较大的高低不平顺使轮轨间产生激烈的冲击，引起极大的相互作用和很大的簧下加速度。

高低不平顺的数值以左右轨高低不平顺平均值来表示：zr(zLzR)/2。 2、轨道水平不平顺

轨道水平不平顺是指左右轨对应点的高差所形成的沿长度方向的不平顺，是由轨道高低不平顺所派生的，是使机车车辆产生侧滚的主要原因。当水平不平顺幅值较大，并接连不断，频率与车辆侧滚的自振频率相近时，车辆将产生较大的左右倾斜振动，使车轮一侧轮载增大，另一侧减小，形成脱轨条件，影响行车安全。轨道水平不平顺如图2-4中zc所示。

水平不平顺的数值以两股钢轨顶面水平偏差沿轨道方向的变化率来表示：

zczLzR。

3、轨道轨距不平顺

轨距不平顺是指左右两根钢轨的轨距沿长度方向上的偏差，轨距不平顺对轮轨磨耗和车辆运行稳定性和安全性影响较大，轨距过大还会引起掉道。轨距不平顺如图2-4中yg所示。

轨距不平顺的数值用实际轨距与名义轨距之差来表示：ygy2y1g。

图2-4 轨道水平不平顺和轨距不平顺

4、轨道方向不平顺

轨道方向不平顺是指左右两根钢轨沿长度方向在横向水平面内呈现的弯曲不直，是由轨道铺设时的原始弯曲、养护和运用中积累的轨道横向弯曲变形等原因造成的。方向不平顺会引起机车车辆的横移-侧滚振动，使得轮对产生很大的横向水平力和侧滚力矩（尤其是强制机车车辆轮对的蛇形运动），对于车辆，很大的横向水平力容易造成脱轨系数超过允许值而导致脱轨。轨道方向不平顺如图2-5中ya所示。

图2-5 轨道方向不平顺

轨道方向不平顺的数值以实际轨道中心线相对于理论中心线的偏差来表示：

ya(y1y2)/2。

5、三角坑不平顺

三角坑不平顺是指轨道一定间隔的水平不平顺的变化量，它是轨道在平面上的扭曲状态。存在三角坑的地方，不仅车辆摇晃严重，而且会出现所谓三轮支撑、一轮悬空（或大大减载）容易出轨的危险状态。

第3章 轨道不平顺功率谱分析

3.1 轨道不平顺状态的评估方法

测得各种轨道不平顺的数据后，需对轨道平顺状态进行科学评价，诊断不平顺的严重程度，以确定列车运行的速度限值和指导维修、养护等工务作业。评定诊断轨道平顺状态好坏和恶化程度的依据，是轨道不平顺对机车车辆响应的影响和经验。我国对轨道不平顺状态的评价方法主要采用局部不平顺幅值超限评分法(即峰

[7]

值扣分法)和轨道质量指数法两种。

[7]

3.1.1 局部不平顺幅值超限评分法

局部不平顺幅值超限评分法从轨道几何尺寸指标、动力学指标的角度出发，根据轨道局部不平顺超限等级，以一公里为单位计算总扣分的方式来评价轨道的质量。检查评定项目包括轨距、水平、高低、轨向、三角坑、车体垂向振动加速度和横向振动加速度共七项。局部不平顺幅值超限评分法把轨道动态几何尺寸允许偏差管理值按线路允许速度分为四级：I级为保养标准，每处扣1分；II级为舒适度标准，每处扣5分；III级为临时补修标准，每处扣100分；Ⅳ级为限速标准，每处扣301分。每公里扣分总数S按下式计算：

SKiTjCij (3-1)

式中，S为每公里扣分数；Ki为各级扣分数，K11，K25，K3100，K4301，即一级、二级、三级、四级扣分分别为1分、5分、100分和301分；Tj为各项的权系数，目前T1~T7均取值为1；Cij为各项不平顺各级偏差的个数，i1~4，

j1~7。

按上式计算结果，按每公里总的扣分数的多少，把轨道质量状态分为如下三种：

A． 优良：S50分； B． 合格：51S300分； C． 失格：S300分。

3.1.2 轨道质量指数（简称TQI）评价方法

《铁道建筑》在1994年以讲座的形式连续刊出了四篇文章，分别介绍了TQI

在国外铁路中的应用概况、TQI的基本概念及计算机处理技术、TQI数据库管理软件的开发、TQI在济南铁路局线路维修工作中的实践和效果。这些成果为TQI技术在我国铁路中的应用打下了坚实的基础。此后，铁路相关科技人员在实践中又不断探索，总结TQI在我国干线中的应用情况，从检测数据的可靠性、TQI在工程中的应用等方面都提出了有益的研究成果，完善了我国铁路轨道状态的管理方法。

TQI从统计学（离散性），物理学（轨道质量均衡性）的角度（相对峰值扣分法）反映轨道状态的恶化程度，是衡量轨道区段整体质量状态的综合指标，可以作为工务部门编制轨道维修、养护计划，指导作业的依据，是对轨道质量状态进行宏观管理和质量控制的重要手段。以200m轨道区段作为单元区段，分别计算单元区段内左、右高低、左、右轨向、轨距、水平、三角坑七项几何参数的标准差。各单项几何不平顺幅值的标准差称为单项指数，七个单项指数之和作为评价该单元区段轨道平顺性综合质量状态的轨道质量指数。其计算公式为：

TQIi1(xijxi)2 (3-2) n式中，i为单项几何参数标准差，i1,2,7；n为单元区段中采样点的个数，200m单元区段中n800点；xij为各项参数在单元区段中采样点的幅值，j1,2,800点；

xi为各项参数在单元区段中连续点幅值xij的平均值。

3.1.3 局部不平顺幅值超限评分法与轨道质量指数评价法的比较

局部不平顺幅值超限评分法能够找出轨道的局部病害及病害的类型、程度和所在位置，作为指导现场紧急补修非常实用，但仅用超限点峰值的大小、超限的数量及扣分多少，不能全面地评价轨道区段的质量状态，比如不能反映周期性不平顺所产生的谐波的影响。轨道质量指数评价法能够判别轨道质量的均衡性，能做出更为符合实际情况的评价。但是这两种方法都是从轨道不平顺幅值的角度出发来评价轨道平顺状态的，因此都具有一些局限性。而轨道不平顺的功率谱密度能清楚地表明某一段轨道不平顺所包含的波长成份及各波长成分的均方值密度，能够提供轨道不平顺幅值和波长两方面的信息，可以对利用局部不平顺幅值超限评分法和轨道质量指数评价法对轨道平顺性进行评定时做出有益的补充。在我国，许多科技人员已经做了大量工作，但是还没有形成较为通用的轨道谱，铁路平顺状态的评定和管理应

用中，轨道谱的应用也十分有限。鉴于此，本文针对轨道谱做了进一步深入的研究，下面简单介绍轨道谱的研究概况。

3.2 轨道谱研究概述

3.2.1 国外铁路轨道谱的研究情况

轨道随机不平顺的统计特征只能依靠线路实地测量获得。自二十世纪60年代中期以来，国外很多国家已开始了实地测试工作，对轨道不平顺的功率谱进行研究。迄今为止，国外轨道谱的研究已经相当完善，极具参考价值。下面简单介绍国外比较典型的轨道不平顺功率谱密度。

[2]

3.2.1.1 美国轨道不平顺轨道谱

美国运输部联邦铁路总署（FRA）制定的铁路安全法规（安全标准）将美国铁路按平顺状态的安全限度和相应的允许速度分为六个等级（98年又增补三个高速等级变为9个等级），并公布了六个等级线路轨道不平顺的功率谱密度，其拟合曲线函数表达式如下（单位：cm2/(rad/m)）：

高低不平顺：

2kAvc (3-3) Sv()222(c)轨向不平顺：

2kAac (3-4) Sa()222(c)水平不平顺和轨距不平顺具有相同的谱密度表达式：

24kAvc (3-5) Sc()Sg()2222(s)(c)式中，S()为轨道不平顺功率谱密度[cm2/(rad/m)]；为轨道不平顺的空间频率

(rad/m)；Av、Aa是粗糙度常数(cm2rad/m)；c、s是截断频率(rad/m)；k是安

全系数，可根据要求在0.25~1.0之间选取，一般取为0.25。

上述各式中的诸参数值见表3-1，可见，由于参数值的不同，不同等级线路的同一名称的谱密度曲线也是不同的。表中还同时列出了根据安全标准制定的不同等级线路所允许的车辆最高运行速度。

表3-1 美国轨道谱的参数值 参数 线路等级 一级 二级 三级 四级 五级 六级 Av/(cm2rad/m) 1.2107 1.0181 0.6816 0.5376 0.2095 0.0339 Aa/(cm2rad/m) 3.3634 1.2107 0.4128 0.3027 0.0762 0.0339 s/(rad/m) 0.6046 0.9308 0.8520 1.1312 0.8209 0.4380 c/(rad/m) 允许最高速度 /（km/h） 货车 客车 0.8245 0.8245 0.8245 0.8245 0.8245 0.8245 16 24 40 48 64 96 96 128 128 144 176 176 美国六个级别线路的功率谱密度（PSD）是在约7万英里的各级线路状态数据库中，每级选取5-10个区段的轨道不平顺检测数据，经计算统计分析得出的。这些区段每个长8-16公里，广泛分布在整个美国，并反映了各铁路公司轨道运营情况和养修状态。图3-1、3-2、3-3为根据美国轨道谱公式(3-3)、(3-4)和(3-5)，经过单位变换（横坐标均变成波长/m，纵坐标为功率谱密度/cm2/(rad/m)），利用MATLAB编程绘出美国六级线路谱，分别为高低、方向和水平和轨距不平顺。（具体程序详见附录1）

106美国六级高低不平顺104功率谱密度/[mm2/(rad/m)]10210010-210-410310210110010-1波长/m 图3-1 美国六级高低不平顺轨道谱

106美国六级方向不平顺104功率谱密度/[mm2/(rad/m)]10210010-210-410310210110010-1波长/m 图3-2 美国六级方向不平顺轨道谱

103美国六级水平及轨距不平顺102功率谱密度/[mm2/(rad/m)]10110010-110-210-310-410310210110010-1波长/m 图3-3 美国六级水平及轨距不平顺轨道谱

3.2.1.2 德国轨道不平顺轨道谱

20世纪80年代初，德国在进行高速列车理论研究时采用了式(3-6)、(3-7)和(3-8)的轨道谱分析式：

高低不平顺：

2Avc (3-6) Sv()222(2)()rc轨向不平顺：

2Aac (3-7) Sa()2222(r)(c)水平不平顺：

2Avb2c2 (3-8) Sc()222222(r)(c)(s)这里，高低、方向不平顺功率谱密度Sv()、Sa()的单位为m2/(rad/m)；水平不平顺由于采用倾角度量，因而其功率谱密度Sc()的单位为1/(rad/m)；为轨道不平顺的空间频率(rad/m)；c、r、s是截断频率(rad/m)；Av、Aa是粗糙度常数

(m2rad/m)；b为左右滚动圆距离之半(m),一般可取0.75m。

德国不平顺轨道谱没有给出轨距不平顺的功率谱密度表达式，但规定轨距不平

顺在-3~3mm范围内变化。一般轨距不平顺与水平不平顺具有相同的功率谱密度表达式，据此，文献[2]提议了一种轨距不平顺功率谱密度的表达式：

Sg()22Agc2222(22)()(rcs) (3-9)

低干扰水平系数、高干扰水平系数和截断频率见表3-2，其中Ag是基于轨距不平顺在-3-3mm范围内变化时经试算得出的参考值。

表3-2 德国轨道谱粗糙度系数及截断频率 c 轨道级别 /(rad/m) r /(rad/m) s /(rad/m) Aa Av Ag /(m2rad/m)/(m2rad/m)/(m2rad/m) 4.03210-7 5.3210-8 低干扰 0.8246 0.0206 0.4380 2.11910-7 高干扰 0.8246 0.0206 0.4380 6.12510-7 1.0810-6 1.03210-7 图3-4、3-5、3-6为根据德国轨道谱公式(3-6)、(3-7)和(3-8)，经过单位变换（横坐标均变成波长/m，纵坐标为功率谱密度/cm2/(rad/m)），利用MATLAB编程绘出的德国低干扰线路谱，分别为高低、方向和水平不平顺谱。（具体程序见附录1）

104德国低干扰高低不平顺102功率谱密度/[mm2/(rad/m)]10010-210-410-610310210110010-1波长/m 图3-4 德国低干扰高低不平顺轨道谱

104德国低干扰方向不平顺102功率谱密度/[mm2/(rad/m)]10010-210-410-610310210110010-1波长/m 图3-5 德国低干扰方向不平顺轨道谱

108德国低干扰水平不平顺106功率谱密度/[mm2/(rad/m)]10410210010-210-410-610310210110010-1波长/m 图3-6 德国低干扰水平不平顺轨道谱

3.2.2 国内铁路轨道谱的研究情况

我国对轨道不平顺功率谱密度的应用较晚，但研究工作起步较早，60~70年代长沙铁道学院、铁道科学院就已用地面检测和轨检车检测等不同方法获得了数百米和数十公里高低、水平、轨向不平顺的轨道谱。此后，很多铁路科技人员对轨道谱展开了深入的研究，取得了一些重要成果，下面简单加以概述。

[2]

3.2.2.1 长沙铁道学院随机振动研究室建议的轨道谱

长沙铁道学院等单位分别于1965年9月、79年10月、82年11月，用地面测量法先后三次在京广线上对轨道不平顺进行了实测，并用傅氏变换法和最大熵谱估计法，求得了京广线郑州－五里堡区间（50kg/m普通线路），京广线猴子石－长沙南站（50kg/m普通线路），京广线桃林一范家园（50kg/m无缝线路）三段约数百米的轨道无轮载作用时，和轮载作用下高低、轨向、水平、轨距四种不平顺的功率谱密度，建议使用的一级铁路干线的轨道不平顺功率谱密度的表达式如下：

轨道高低不平顺：

f28.879101 (3-10) Sv(f)2.755104227f2.52410f9.61103轨道方向不平顺：

f29.701102 (3-11) Sa(f)9.404104f3.768102f22.6661053轨道水平不平顺：

f26.346103 (3-12) Sc(f)5.100104226f3.15710f7.791108轨道轨距不平顺：

f23.863102 (3-13) Sg(f)7.001104225f3.35510f1.464103以上各式中，轨道不平顺功率谱密度S(f)的单位为mm2/(1/m)；f是空间频率

(1/m)。

3.2.2.2 铁道科学研究院建议的轨道谱

1999年7月铁道科学院完成的铁道部重点课题“我国干线轨道不平顺功率谱”，对轨道谱的数据采集、处理、计算、分析进行了较全面深入地研究，提出了我国主要干线和不同轨道结构、质量状态以及曲线、桥梁、焊缝等特殊地段的轨道谱，轨道高低、水平、轨向不平顺功率谱密度采用系数不同的同一解析式表达：

A(f2BfC) (3-14) S(f)432fDfEfFfG式中，S(f)为轨道不平顺功率谱，f为空间频率，AG为轨道不平顺功率谱密度的特征系数，对不同线路和不同类型的轨道不平顺有不同的数值。表3-3和表3-4分别给出了我国京哈、京广、京沪三大提速干线和60kg/m钢轨超长无缝线路轨道的特征参数，后者实质上是基于1998年郑武（郑州—武昌）线200km/h以上高速试验段轨道状态的测试结果而拟合得出的，因而也可称之为郑武线高速试验段轨道谱。

表3-3 我国京沪、京广、京哈三大干线轨道谱的特征参数 参数 A B C D E F 左轨高低 1.1029 -1.4709 0.5941 0.8480 3.8016 -0.2500 右轨高低 0.8581 -1.4607 0.5848 0.0407 2.8428 -0.1989 左轨轨向 0.2244 -1.5746 0.6683 -2.1466 1.7665 -0.1506 右轨轨向 0.3743 -1.5894 0.7265 0.4353 0.9101 -0.0270 水平 0.1214 -2.1603 2.0214 4.5089 2.2227 -0.0396 注：测量时三大干线提速目标为160km/h。

图3-7、3-8、3-9为根据铁科院轨道谱公式(3-14)，经过单位变换（横坐标均变成波长/m，纵坐标为功率谱密度/cm2/(rad/m)），利用MATLAB编程绘出中国三大干线谱，分别为高低、方向和水平不平顺，谱线均用左轨参数作出。（具体程序详见附录1）

102G 0.0112 0.0094 0.0052 0.0031 0.0073 三大干线左轨高低功率谱密度/[mm2/(rad/m)]10110010-110-2101100波长/m 图3-7 三大干线左轨高低不平顺轨道谱

102三大干线左轨轨向功率谱密度/[mm2/(rad/m)]10110010-110-2101100波长/m 图3-8 三大干线左轨轨向不平顺轨道谱

102三大干线水平功率谱密度/[mm2/(rad/m)]10110010-110-2101100波长/m 参数 左轨高低 右轨高低 左轨轨向 右轨轨向 水平 图3-9 三大干线水平不平顺轨道谱

表3-4 郑武线高速试验段轨道谱的特征参数 A B C D E F 0.1270 -2.1531 1.5503 4.9835 1.3891 -0.0327 0.3326 -1.3757 0.5497 2.4907 0.4057 0.0858 0.0627 -1.1840 0.6773 2.1237 -0.0847 0.0340 0.1595 -1.3853 0.6671 2.3331 0.2561 0.0928 0.3328 -1.3511 0.5415 1.8437 0.3813 0.2068 G 0.0018 -0.0014 -0.0005 -0.0016 -0.0003 注：试验段最高试验速度达到240km/h。

图3-10、3-11、3-12为根据铁科院轨道谱公式(3-14)，经过单位变换（横坐标均变成波长/m，纵坐标为功率谱密度/cm2/(rad/m)），利用MATLAB编程绘出的郑武线高速试验段轨道谱，分别为高低、方向和水平不平顺。（具体程序详见附录1）

102郑武线高速试验段左轨高低不平顺101功率谱密度/[mm2/(rad/m)]10010-110-210-3101100波长/m 图3-10 郑武线高速试验段左轨高低不平顺轨道谱

102郑武线高速试验段左轨轨向不平顺功率谱密度/[mm2/(rad/m)]10110010-110-2101100波长/m 图3-11 郑武线高速试验段左轨轨向不平顺轨道谱

102郑武线高速试验段水平不平顺功率谱密度/[mm2/(rad/m)]10110010-110-2101100波长/m 图3-12 郑武线高速试验段水平不平顺轨道谱

3.2.2.3 轨道短波不平顺功率谱

以上的轨道不平顺功率谱的波长范围都在几米到几十米范围内，所以它们一般只能满足机车车辆和桥梁结构的低频随机振动分析，无法满足轨道结构随机振动研究的需要，因为簧下质量和轨下结构的振动主频可达数百到数千赫兹。为此，1988年铁科院王澜在博士论文中对我国石太（石家庄—太原）线的轨道垂向短波不平顺进行了实测，测量方法是采用地面测量方式，用Colmar钢轨磨耗测量仪进行测量，经回归分析，提出了我国50kg/m钢轨线路垂向短波不平顺的功率谱密度函数的表达式，其波长范围为0.01~1m，即

S(f)0.036f3.15 (3-15)

式中，S(f)的单位为mm2/(1/m)；f是空间频率(1/m)。

3.3 国内外轨道谱比较分析

由于我国至今没有轨道谱标准，3.2.2节介绍了根据实测轨道几何参数得出的我国几种轨道谱的研究结果，其中最具代表性的当属京沪、京广、京哈三大干线轨道谱和郑武线高速试验段轨道谱。然而，尚不清楚它们的基本特征是什么？与国外典型轨道谱有何差异？为了便于我们在动力学分析中正确选取轨道随机不平顺激扰，有必要探明这些问题。

3.3.1 普通线路轨道谱的比较

首先比较常规速度运行条件下的轨道谱。我国普通线路以京沪、京广、京哈三大提速干线为例，其轨道谱适用于160km/h以下速度；国外普通线路轨道谱选取美国AAR标准五级线路谱（适用于144km/h以下速度）及六级线路谱（适用于176km/h以下速度）。为了方便比较，必须先将它们换算成同一单位（此处均换算成波长(m)－功率谱密度[mm2/(rad/m)]），对于公式(3-3)、(3-4)的美国谱和公式(3-14)的中国干线谱换算过程如下：

设空间波长为L，则可得： 对于美国谱：

高低不平顺：

222()2kAvcL2LckAvSv(L)22 (3-16) 2222L(c)()2[()2()2]2[(c)21]LLLcLkAv21，f，根据能量守恒有S(L)dLS()d，则LL方向不平顺：

222()2kAacL2LckAaSa(L)22 (3-17) 2222L(c)()2[()2()2]2[(c)21]LLLcLkAacm2rad/m=cm2/m。 由此可得Sv(L)、Sa(L)的单位为

rad对于三大干线谱：

11A[()2BC]A(fBfC)1LLS(L)4fDf3Ef2FfG(1)4D(1)3E(1)2F1GL2 (3-18)

LLLLA(L2BL3CL4)2L(1DLEL2FL3GL4)2由于根据上式无法推倒出S(L)的单位，故利用S(L)S(f)dfS(f)2，由S(f)的dLL单位mm2/(1/m)可得到S(L)的单位为mm2m/m2=mm2/m。

经过单位换算后，在MATLAB中编写程序可以得到以下普通线路轨道谱比较图：（具体程序详见附录1）

103普通线路高低不平顺 102美国五级线路谱美国六级线路谱中国三大干线谱功率谱密度/[mm2/(rad/m)]10110010-110-210-3 101100波长/m 图3-13 普通线路轨道高低不平顺功率谱密度

103普通线路方向不平顺 102美国五级线路谱美国六级线路谱中国三大干线谱功率谱密度/[mm2/(rad/m)]10110010-110-210-3 101100波长/m 图3-14 普通线路轨道方向不平顺功率谱密度

图3-13、3-14分别给出了1~30m波长范围内（当时我国轨检车仅能测量1~30m波长范围），美国五、六级线路谱和中国三大干线谱的高低不平顺功率谱密度的对比结果和方向不平顺功率谱密度的对比结果。

由图3-13可见，我国三大干线轨道谱的高低不平顺总体上介于美国五级谱和六级谱之间。具体比较结果发现：在1~25m波长范围内，中国三大干线轨道谱高低不平顺数值小于美国五级谱而大于美国六级谱；在25~30m波长范围内，我国三大干线谱还略优于美国六级谱。由于长波不平顺主要影响车辆运行平稳性，短波不平顺对轮轨动力作用有重要影响，由此可推知，在我国三大干线谱激扰下，中低速列车垂向平稳性要优于美国五级谱下的舒适度，但较美国六级谱的舒适度要差。

由图3-14可以看出，我国三大干线轨道谱的方向不平顺总体大于美国五级谱和六级谱数值，这表明线路方向几何状态较差。具体比较结果如下：与美国六级谱相比，在1~25m波长范围内，我国三大干线谱密度值明显较大，而在25m以上的少许波长下，三大干线谱密度值要小些；与美国五级谱相比，在7~20m波长范围内，二者的功率谱密度非常接近，波长大于20m时，三大干线谱优于美国五级谱，波长小于7m时，前者差于后者。因此可以预计，车辆在我国三大干线上的横向平稳性与美国五级谱条件下的基本相当，但差于美国六级谱下的舒适性，而轮轨横向动力作用性能均不及美国五级、六级线路情形。

3.3.2 高速线路轨道谱的比较

关于高速运行条件下的轨道谱，国外最典型的标准要数德国轨道谱，包括低干扰（适合250km/h及其以上速度）和高干扰谱（适合250km/h以下速度）两种，我国尚无高速轨道谱标准，1998年曾在郑武线上按高速运行要求改造了一段所谓高速试验段（试验段从许昌至小商桥之间，曾创下当时中国铁路最高试验速度

240km/h），通过轨检车测量并拟合了轨道谱，这里就以此作为比较对象。同样对

德国轨道谱进行单位变换，得

高低不平顺：

222()22AvcLLcSv(L)2222222(2r)(c)(()2()2)(()2()2)LLrLLc (3-19)

AvAv11112L2L2(2)(22)c2LLrLLc方向不平顺：

222()2AacL2LcSa(L)2222222(2)()rc(()2()2)(()2()2) (3-20) LLrLLcAaAa11112L2L2c(22)(22)LLrLLcm2rad/m=m2/m=m。 由此可得Sv(L)、Sa(L)的单位为

11radm422mm经过单位换算后，在MATLAB中编写程序可以得到以下高速线路轨道谱比较图：（具体程序详见附录1）

104高速线路高低不平顺 10功率谱密度/[mm2/(rad/m)]3 德国低干扰轨道谱德国高干扰轨道谱中国高速试验段谱10210110010-110-210-3 102101100波长/m 图3-15 高速线路轨道高低不平顺功率谱密度

103高速线路方向不平顺 10功率谱密度/[mm2/(rad/m)]2 德国低干扰轨道谱德国高干扰轨道谱中国高速试验段谱10110010-110-210-310-4 102101100波长/m 图3-16 高速线路轨道方向不平顺功率谱密度

图3-15、3-16分别给出了德国低干扰谱、高干扰谱和郑武段高速轨道谱的高低不平顺及方向不平顺功率谱密度在1~30m波长范围内的对比结果。文献［2］中通过地面测量得到了30m以上的长波高低不平顺，这是因为长波不平顺对高速列车运行的舒适度影响很大。

根据图3-15和图3-16的结果可知：

对于高低不平顺，中国高速试验段轨道谱与德国高干扰谱互有优劣，在7~10m波长范围内，二者的功率谱密度值相差不多；在10~30m波长范围内，中国谱密度值要比德国高干扰谱的功率谱密度值小；而在1~7m波长范围内，中国谱密度值要比德国高干扰谱的功率谱密度值大一些。中国高速试验谱与德国低干扰谱相比，从波长1~30m范围内，中国高速试验谱明显要差。由于长波不平顺对高速列车运行舒适性有重要影响，短波不平顺对高速列车的动力特性影响较大，故可以推断，对于垂向舒适度，中国高速试验谱要优于德国高干扰谱，而比德国低干扰谱差；而对于轮轨垂向动力作用性能，中国高速试验谱不仅差于德国高干扰谱，更差于德国低干扰谱。

对于方向不平顺，在1~15m波长范围内，中国高速试验段轨道谱要差于德国高

干扰谱和低干扰谱，尤其对数米以下波长非常明显；在15~30m波长范围内，中国高速试验段轨道谱略优于德国高干扰轨道谱，但明显劣于德国低干扰谱。同样可以推断，在德国低干扰谱激扰下列车横向运行舒适性优于其它两种谱下的舒适性，而中国高速试验谱激扰下的车辆横向舒适性接近于德国高干扰谱激扰结果。

3.3.3 结论

① 我国三大干线谱高低不平顺在25m以下波长范围内，比美国五级线路谱要好，比美国六级线路谱要差。

② 我国三大干线轨道谱的方向不平顺在7~20m的波长范围内基本与美国五级线路谱相当，在1~7m的波长范围内尚不如美国五级线路谱好，在25m以下波长范围均比美国六级线路谱明显要差。

③ 中国高速试验段轨道谱的高低不平顺在7~30m波长范围内，其高低不平顺和方向不平顺基本与德国高干扰谱相当，而在1~7m波长范围内的高低和方向不平顺，尚不如德国高干扰谱，离德国低干扰谱相差很远。

第4章 轨道谱估计

在第1章已经提到，轨道不平顺影响因素的多样性和随机性决定了轨道不平顺的随机性，沿线路长度方向的轨道不平顺可以看作一个随机过程，轨道不平顺是里程的随机函数，任一特定区段的轨道不平顺可以看作随机过程的一个样本。所以轨道不平顺不能用一个确定的数学表达式来描述，只能应用随机信号理论加以描述。

[3]

4.1 随机过程及其特征描述

4.1.1 随机过程

随机过程通常用X(,t)表示，定义为两个自变量和tT的一个集合函数。其中，是随机变量族的样本空间，T是参数t的集合。随机过程X(,t)能以两种方法描述：

1、表示成随机变量族{X(,ti),tiT}，对任何固定的tiT，X(,ti)是随机变量。

2、表示成T上的一个函数集{X(i,t),i}，对任何固定的i，X(i,t)是一个样本函数，或称为过程的一个实现。

4.1.2 平稳随机过程

1．严平稳过程——从n维概率分布函数出发

设{X(t),tT}是一个随机过程，如果对于任意的T，过程X(t)与

X(t)(tT)有同样的分布，即

F(x1,x2,,xn;t1,t2,,tn)F(x1,x2,,xn;t1,t2,,tn) (4-1)

对一切有限集{ti}T和任意T都成立，则称为X(t)严平稳过程（也称狭义平稳过程或强平稳过程）。

严平稳过程描述的物理系统，其任意的有限维分布不随时间改变。注意，只是分布不随时间变化，并不涉及服从什么分布。

2．宽平稳过程——只考虑一阶矩和二阶矩 设{X(t),tT}是一个随机过程，如果有

① X(t),(tT)为二阶矩过程，即

E[X(t)]， tT (4-2)

2② E[X(t)]为常数，且相关函数只与时间间隔有关，而与起点无关，即

RX(t,t)RX() (4-3)

其中，RX()为随机过程的自相关函数，则称X(t)为一个宽平稳过程（也称广义平稳过程或弱平稳过程），简称平稳过程。

宽平稳过程不一定是严平稳过程；反过来，严平稳过程也不一定是宽平稳过程，因为宽平稳过程必须是二阶矩过程。

3．各态历经性

如果随机过程任一样本函数的时间平均等于过程的集合平均，则称随机过程是各态历经的或各态遍历的，即，如果有

At[x(t)]E[X(t)] (4-4)

则称X(t)为均值各态历经过程；如果有

At[x(t)x(t)]E[X(t)X(t)]RX() (4-5)

则称X(t)为自相关各态历经随机过程。

可见，对各态历经随机过程，可以用一个样本函数的时间平均计算随机过程的集合平均。这一点为实际操作带来很大的方便，因为对一个样本过程进行长时间统计比同时对许多样本进行统计要容易实现。本文计算所涉及的随机数据均看作广义平稳和各态历经的随机过程。

4.1.3 随机信号的相关函数

平稳随机信号的自相关函数Rx(m)定义为

Rx(m)E[x*(n)x(nm)] (4-6)

式中x(n)是平稳随机信号，\"*\"代表取共轭，如果是x(n)实随机信号，则上式成为

Rx(m)E[x(n)x(nm)] (4-7)

平稳随机信号的相关函数具有如下性质：

① 若x(n)为实信号，则Rx(m)为实偶函数，即

*Rx(m)Rx(m) (4-8)

Rx(m)Rx(m) (4-9)

若x(n)为复信号，则Rx(m)共轭偶对称，即

*Rx(m)Rx(m) (4-10)

本文所用到的随机信号均为实随机信号。

② 在m0时，Rx(m)取得最大值，即

Rx(0)Rx(m)

且Rx(0)就是序列的平均功率，即

Rx(0)E[x2(n)] (4-11)

③ 维纳-辛钦定理

随机信号自相关函数的傅立叶变换是信号的功率谱密度。如果用Sx(ejw)表示随机信号序列x(n)的功率谱密度，则有

Sx(e)F[Rx(m)]Rx(m)F1[Sx(ejw)]12jwmR(m)exjwm (4-12)

jSx(ejw)ejwmd (4-13)

在式中令m0，可得

1Rx(0)2S(ex)d (4-14)

因为Rx(0)E[x2(n)]为信号的平均功率，所以将Sx(ej)称为随机信号的功率谱密度(PSD)或功率谱密度谱，简称功率谱。

4.1.4 随机信号的功率谱

对确定性能量信号，可以用FFT（快速傅立叶变换）作频谱分析，得到频域特性。对于平稳随机信号，因为是无限能量的信号，故其傅立叶变换不存在（在z平面的单位圆上不满足绝对可和的条件）。由于自相关函数序列Rx(m)是一个能量有限的确定性序列，故能满足序列傅立叶变换绝对可和的条件。由维纳-辛钦定理可知，

对序列Rx(m)求傅立叶变换得到的就是序列的功率谱Sx(ej)。

平稳随机信号的功率谱具有如下性质：

① 不论x(n)是实序列还是复序列，Sx(ej)都是的实函数。 ② 如果x(n)是实序列，Sx(ej)具有偶对称性。

③ Sx(ej)对所有的都是非负的，且是的周期函数，周期为2。

4.2 功率谱估计的各种方法

根据维纳-辛钦定理，广义平稳随机过程的功率谱是自相关函数的傅立叶变换，它取决于无限多个自相关函数值。但对于许多实际应用问题，可资利用的观测数据往往是有限的，所以要准确计算功率谱通常是不可能的。比较合理的目标是设法得到功率谱的一个好的估计值，这就是功率谱估计。也就是说，功率谱估计是根据平稳随机过程一个实现的有限个观测值，来估计随机过程的功率谱密度（PSD）。这里涉及到两个问题，怎样评价一个估计是好的估计？怎样得到好的估计？

功率谱估计评价指标包括客观度量和统计度量。在客观度量中，谱分辨率特性是一个主要指标。谱分辨率是指估计谱对真实谱中两个靠的很近的谱峰的分辨能力。统计度量是指估计的偏差、方差、均方误差、一致性等评价指标。但需要注意的是，对统计特性的分析方法只适用于长数据记录。所以，利用统计度量对不同的谱估计方法进行比较是不妥当的，只能用来对某种谱估计方法进行描述，并且一般只用来描述古典谱估计方法，因为现代谱估计方法往往用于短数据记录情况。

至于怎样得到好的估计，这就是后面将要介绍的各种谱估计方法。这些方法主要分为两大类。通常，将以傅立叶分析为理论基础的谱估计方法叫做古典谱估计方法或经典谱估计方法；而把不同于傅立叶分析的新的谱估计方法叫做现代谱估计或近代谱估计。本文主要介绍了古典谱估计和最大熵谱估计，分述如下：

4.2.1 古典谱估计

古典谱估计主要有相关法（间接法）和周期图法（直接法）两种，以及由此派生出来的各种改进算法。 4.2.1.1 相关法谱估计

相关法谱估计是以相关函数为媒介来计算功率谱，所以又叫间接法。它的理论

基础是维纳-辛钦定理，因为是由Blackman和Turkey提出的，所以又叫BT法。其具体步骤如下：

第一步：由获得的N点数据构成的有限长序列xN(n)来估计自相关函数Rx(m)序列，即

1N1Rx(m)xN(n)xN(nm) (4-15)

Nn0这一步要注意三个问题：①Rx(m)是双边序列，自变量的取值范围应为

m(N1),1,0,1,N,(。对实数序列xN(n)，由自相关函数的偶对称性，只

需求出m0,1,,(N1)的Rx(m)，另一半也就知道了。②因为得到的只是x(n)的N个观测值xN(0),xN(1),xN(N1)，所以对nN时的x(n)的值只能假设为零。③为了较少地使用已知数据段之外的数据，通常取mN。所以，计算自相关函数时只需求出m0,1,,(m1)的Rx(m)（这里MN），然后利用Rx(m)Rx(m)求出另一半。

第二步：由自相关函数的傅立叶变换求功率谱，即

Sx(e)M1jm(M1)Rx(m)ejm (4-16)

以上两步中经历了两次截断，一次是估计Rx(m)时仅利用了x(n)的N个观测值

xN(n)，这相当于对x(n)加矩形窗截断。该窗是加在数据上的，一般称为加数据窗。另一次是估计Sx(e)时仅利用了从(M1)到(M1)的Rx(m)，这相当于对Rx(m)加矩形窗截断，将Rx(m)截成(2M1)长，这称为加延迟窗。式中的Rx(m)和Sx(ej)分别表示对Rx(m)和Sx(ej)的估值。由于MN，使得式(4-16)的运算量不是很大，因此，在FFT问世之前，相关法是最常用的谱估计方法。

当FFT问世后情况有所变化。因为截断后的xN(n)可视作能量信号，由相关卷积定理，可将式(4-15)写为

Rx(m)j1[xN(m)xN(m)] (4-17) N对上式两边取(2N1)点的DFT，并将时域卷积变为频域乘积，有

Sx(ej)112*[X2(k)X(k)]X(k) (4-18) N12N12N1NN于是可得用FFT求Rx(m)的完整方案如下：

①对N长的xN(n)充(N1)个零，成为(2N1)长的x2n1(n)。 ②求(2N1)点的FFT，得

2N2X2N1(k)xn02n1(n)W2nkN1 (4-19)

③求

12X2N1(k)。 N④求(2N1)点的IFFT，得

N1112Rx(m)X2N1(k)W2Nnk1 (4-20)

2N1k(N1)N⑤由式(4-16)便可求出Sx(ej)。

上面的相关运算中，充零的目的是为了能用圆周卷积代替线性卷积，以便采用快速卷积算法。

4.2.1.2 周期图法谱估计

周期图法又称直接法，其具体步骤如下：

第一步：由获得的N点数据构成的有限长序列xN(n)直接求傅立叶变换，得频谱XN(ej)，即

XN(e)xN(n)ejm (4-21)

jn0N1第二步：取频谱幅度的平方，并除以N，以此作为对x(n)真实功率谱Sx(ej)的估计，即

Sx(ej)21XN(ej) (4-22) N4.2.1.3 古典谱估计的改进

相关法和周期图法都是用获得的N个数据对随机过程进行谱估计，它隐含着对

无限长数据加了一个矩形窗截断。时域中与矩形窗函数相乘对应于频域中与矩形窗频谱相卷积，就这一点来说，估计谱相当于真实谱与矩形窗频谱相卷积的结果。矩形窗频谱等于矩形序列RN(n)的傅立叶变换，即

)jN11e2e2 (4-23) WR(e)RN(n)eej1enn0sin()2N)/sin()，对功率谱有影响的是矩形窗频谱的幅度谱sin(该函数与功率谱相卷积22jjnN1jnjNsin(N必然使所得的估计谱不同于真实的功率谱，因为函数不同于函数（任何函数与函数相卷积形状不变）。

为了使估计谱逼近真实谱，应设法使窗谱幅度函数逼近函数。矩形窗频谱的幅度函数为sin(N)/sin()，它与函数相比存在两方面的差异，一是主瓣不是无

22限窄，二是有旁瓣。由于主瓣不是无限窄，真实的功率谱与主瓣卷积后将使功率向附近频域扩散，使信号模糊，降低了谱分辨率，主瓣越宽分辨率越低。由于存在旁瓣，又会产生两个后果，一是功率谱主瓣内的能量“泄漏”到旁瓣使谱估计的方差增大，二是旁瓣卷积后得到的功率谱完全属于干扰。严重情况下，强信号与旁瓣的卷积可能大于弱信号与主瓣的卷积，使弱信号淹没在强信号的干扰中，无法检测出来，这是古典谱估计的主要缺点。

为了减少泄漏和提高谱估计的分辨率，改善窗的形状是必要的。对窗函数幅度谱的总要求是希望它的主瓣尽可能窄，同时旁瓣尽可能小，以减少能量外泄，改善方差，提高谱估计的分辨率。但实现时两者往往是矛盾的，任何减小旁瓣的努力都要以牺牲主瓣宽度为代价，反之亦然。两者只能互换，不可兼得。因此古典谱估计无法克服谱分辨率低的缺点。如果用统计量度量进行评价，可以得出相关法谱估计和周期图法谱估计都不是对Sx(ej)的一致估计，主要问题是方差大。于是人们研究各种改进方法以减少方差，尤其是周期图法。改进周期图法的主要途径是平滑和平均。

平滑是用一个适当的窗函数W(e)与算出的功率谱Sx(ej)进行卷积，使谱线

j平滑。这种方法得出的谱估计是无偏的，方差也小，但分辨率下降。平均就是将截取的数据段xN(n)再分成L个小段，分别计算功率谱后取功率谱的平均，这种方法

叫Bartlett方法。因为L个平均的方差比随机变量的单独方差小L倍，所以当

L时，L个平均的方差趋于零，可以达到一致谱估计的目的。

现在比较常用的方法是Welch法，又叫加权交叠平均法，简记为WOSA法。这种方法以加窗（加权）求取平均，以分段重叠求得平均，因此集平均与平滑的优点于一体，但同时也不可避免地带有两者的缺点。其主要步骤如下：

① 将N长的数据段分成L个小段，每小段M点，相邻小段间交叠M/2点，于是段数为

LNM/2 (4-24)

M/2② 对各小段加同样的平滑窗w(n)后求傅立叶变换，得

Xi(e)xi(n)w(n)ejn， i1,,L (4-25)

jn0M1③ 求各小段功率谱的平均，得

21L1Sx(e)Xi(ej) (4-26)

Li1MUj1这里，UMM1n0w(n)代表窗函数平均功率，所以MUw(n)是M长窗函数

22n0M1w(n)的能量。

Welch法对数据分段加非矩形窗，因为窗函数在其边沿为零，从而减小了分段时的截断效应。分段平均减小了由数据样本本身的随机性带来的方差，段数越多方差越小，但分辨率下降。分段多了每段的点数必然减少，分段时允许数据有部分重叠，可以在段数较多的情况下拉长每小段长度，有利于平滑过渡。Welch法得出的

Sx(e)是无偏谱估计，当段数L增大时方差减小，Sx(ej)趋于一致估计。

j古典谱估计出现较早，运算量较小，物理概念明确，便于工程实现，对长数据记录来说，还是比较实用的。但古典谱估计有一些难以克服的缺点，可以总结如下：

① 谱的分辨率较低。

② 加窗的坏影响不可避免。较宽的主瓣降低分辨率，较大的旁瓣有可能掩盖真实谱中较弱的成分，或是产生虚假的谱峰。没有一个窗函数能使谱估计在方差、偏差和分辨率各方面同时得到改善，使用窗函数只是一种折中的技巧，不是解决问

题的根本办法。

③ 古典谱不是真实谱的一致估计。

上述缺点促使人们寻求更好的谱估计方法，这就导致了现代谱估计的产生。现代谱估计重点研究包含短序列在内的各种高分辨率的、有效的谱估计方法。

4.2.2 最大熵谱估计

古典谱估计隐含着数据窗以外的序列值为0的假设，显然这是不合理的。如何利用有限的数据记录，尽可能得到PSD的良好估计？Burg提出外推给定的有限长自相关序列，使之变成无限长序列，再由维纳-辛钦定理计算功率谱。

假定已知自相关序列{Rx(0),Rx(1),,Rx(p)}，现通过合理的外推求得

Rx(p1),Rx(p2),，并要保证外推后的自相关矩阵是非负定的。一般有无限多种可能的外推方法，都能得到比较合适的自相关序列。Burg证明了如果外推后的自相关序列所对应的时间序列具有最大熵，那么这种外推方法才是最合理的。为了说明这种方法的合理性，先简单介绍一下熵的概念。

在香农（Shannon）信息论中，离散性随机变量X的熵定义为

H(X)P(xi)lnP(xi)E[lnP(xi)] (4-27)

i1N式中P(xi)表示Xxi这一事件发生的概率，lnP(xi)是信息量的定义，信息量是解除事件不确定性所需要信息的度量。E[]表示求数学期望，上式表明，熵是平均不确定性的度量，它在数值上等于平均信息量。对于必然事件，由于其发生的概率等于1，故其信息量等于零，对应的熵等于零。越是小概率事件，其信息量越大，对应的熵值越大。

根据上面的概念不难想象，如果外推后的自相关序列所对应的时间序列具有最大熵，则意味着在具有已知的p1个自相关值的所有时间序列中，该时间序列的不确定性最大，将是最随机或最不可预测的。统计学认为，最大熵是最合理、最自然、最无主观性的假定。这正是选择最大熵准则外推自相关序列的合理性所在。

在此思路下，应用解有约束优化问题的拉格朗日乘子法，可以推出高斯随机过程的最大熵功率谱为

S(e)j21akejkk1p2 (4-28)

式中的ak可以根据Yule-Walker方程由已知的p1个自相关函数值求得，也就是AR模型的参数。

于是可以得出结论，对高斯随机过程，在已知{Rx(0),Rx(1),,Rx(p)}的情况下，其最大熵谱与其AR模型谱是一致的。因为最大熵谱是建立在自相关函数外推基础上的，所以AR模型谱也等效于一个外推后的自相关函数的谱。需要注意的是，如果不是高斯型随机过程，其最大熵谱与AR模型谱是不一样的。

4.2.3 谱估计方法的比较分析

图4-1、4-2、4-3、4-4、4-5所示是用周期图法、自相关法、Bartlett法、Welch法和最大熵法对轨道实测信号数据（京广线K80~K90区间下行线左钢轨高低不平顺检测数据）进行谱估计的结果。（具体程序详见附录1）

由于我国轨检车只能检测1~50m波长范围的轨道不平顺，其对应的频率范围为0.83~41.7Hz，故在图中我们只比较该频率范围内的功率谱密度值。又因为频率20Hz之后的功率谱密度值基本不变，故可以认为频率从20~41.7Hz的数据段不可用。故只需比较0.83~20Hz频率段内的功率谱，由各图中曲线可以看出，随着频率的增大，功率谱密度值都在减小，频率为0.83Hz处的功率谱密度约为100mm2，频率为20Hz处的功率谱密度值约在10-1mm2处徘徊。因此可以判断各种方法功率谱值的范围基本一致。

另一方面，随着使用方法的优化，曲线变得越来越平滑，分辨率越来越高。周期图法和自相关法作为古典方法，其精度不高，曲线毛刺很多，难以分辨各频率对应幅值，但相比之下，自相关法的精度要稍高于周期图法。Bartlette法和Welch法所作曲线明显比前面两者平顺的多，二者曲线差异不是很明显，Welch法精度稍高，这是因为Bartlett法只是将数据分段，各数据段加矩形窗，而Welch法是将数据有重叠分段，各数据段加了汉明窗，数据被更进一步的平均和平滑。图4-5是用最大熵法作出的曲线，它是一种现代谱估计方法，毫无疑问，不论从分辨率和平滑度上都是很出色的，它是对数据进行功率谱估计最合理、最自然、最无主观性的 方法。

104周期图法谱估计102100功率谱密度(mm2)10-210-410-610-810-1010-1100101102频率(Hz) 图4-1 周期图法谱估计图

103自相关法102功率谱密度(mm2)10110010-110-110-2100101102频率(Hz) 图4-2 自相关法谱估计图

103bartlett法谱估计102功率谱密度(mm2)10110010-110-210-110-3100101102频率(Hz) 图4-3 Bartlett法谱估计图

1010103welch法谱估计21功率谱密度(mm2)10101010100-1-2-3-410-110-5100101102频率(Hz) 图4-4 Welch法谱估计图

1010103最大熵法谱估计21功率谱密度(mm2)10101010100-1-2-3-410-110-5100101102频率(Hz) 图4-5 最大熵法谱估计图

4.3 实测轨道谱与现有轨道谱的比较

为了分析该路段不平顺状态，可以用该路段的轨道谱与现有轨道谱进行比较，本文将实测高低不平顺轨道谱与美国六级高低不平顺线路谱进行了比较，分析了该路段的不平顺状况，提出了一些养护维修意见。

由于实测轨道谱的横、纵坐标分别为频率/Hz和功率谱密度/mm2，故必须将现有的美国高低不平顺谱（如式(3-3)所示）进行单位转化，转化过程如下：

因为V，2f，故得：

2f，根据能量守恒有S(f)dfS()d，可V22fc2()kAVfc2VVV (4-29) SV(f)2222f22f22fc2()[()()]2f(ffc)VVVkAV式中：为空间频率(rad/m)，为角频率(rad)，f时间频率(Hz)，V为车速(m/s)，

cm2rad1mms2s=cm2s=cm2/Hz，只需这里为取为100km/h。由此得到的SV(f)单位为

11rad22ss再给SV(f)乘以100f，则其单位就变为mm2，横坐标为频率/Hz。查表3-1可知Av、

c的值，经过上面的转换便可求得fc，利用MATLAB编程实现这一过程，就可以得到下面的比较图：

104实测高低不平顺轨道谱与美国六级高低不平顺谱的比较 实测高低不平顺轨道谱线美国六级高低不平顺谱线102功率谱密度(mm2)10010-210-410 -110-6100101频率(Hz) 图4-6 实测高低不平顺轨道谱与美国六级

高低不平顺谱的比较图

由4.2.3节的分析比较结果可知，只需考虑0.83~20Hz频率段内的功率谱，故在比较实测高低不平顺轨道谱与美国六级高低不平顺谱时，只需分析该频率段内的功率谱密度，由上图明显可以看出，在0.83~20Hz频率段内，实测的功率谱密度值总体上大于美国六级功率谱值，只有少数值在落在了美国六级高低不平顺谱线以下，因此可以得出结论，京广线K80~K90区间下行线左钢轨高低不平顺总体状况不及美国六级线路来的好，功率谱数值上相差一个数量级，轨道高低质量较差，为了使列车运行更快速、更平稳，还需进一步改进该路段的线路状况。

4.4 线路平顺性趋势分析

如上章所述，评价线路平顺性的方法很多，但是利用功率谱密度来分析是比较宏观、比较合理的方法。本节将运用实测数据得到的功率谱曲线来分析该线路在一定时间内平顺性的变化趋势。

针对京广线K80-K90区间下行线左钢轨高低不平顺04年到07年的部分检测数据，调入已有的功率谱估计函数（运用Welch法），将所有谱线画到一张图纸上，得到下图：

1010103趋势分析 2104050607功率谱密度(mm2)10101010100-1-2-3-410 -110-5100101102频率(Hz) 图4-7 04~07年京广线K80-K90区间下行线左钢轨

高低不平顺趋势分析

图中将各个年份轨道谱线画成不同颜色以便区分，由于轨检车所能测试的波长范围为1~50m，取列车行车速度为150km/h，则对应的频率范围为41.7~0.83Hz，故只需分析该频率段内的谱线。由图知，在频率0.83~4Hz范围内，即波长50~10.4m，功率谱密度值变小，轨道高低平顺状况逐年变好，也就是说，在较长波范围，高低平顺性趋于良好，车辆运行的平稳性逐年提高；在频率4~10.5Hz范围内，即波长10.4~3.97m，06年功率谱密度值最小，07年次之，05年最差，因此可得出在该波长范围内，07年较06年恶化，但比04、05年良好些；在频率10.5~41.7Hz范围内，即波长3.97~1m，07年功率谱密度值最小，较前三年都平顺些，但同时可以看出，06年部分月份轨道平顺性较差，在得到维修后，平顺性得到改善。因此，总体来讲，04~07年，该路段轨道高低不平顺整体状态是向好的方向发展的，但在波长10.4~3.97m范围里，不平顺幅度有所增大，这也需要铁路工务部门及时维护以改善这一波长范围的高低平顺度，使轨道状况得到改善，行车更加平稳、舒适。

第5章 轨道不平顺数值模拟

对铁路车辆—轨道系统进行深入的动力学分析以及对车辆进行振动台试验，都需要可靠的线路不平顺样本。考虑到轮轨接触等非线性因素，轨道不平顺随机过程的数值模拟也是车辆—轨道系统动力研究的重要工具，其模拟质量直接影响研究的有效性。国内外常用的轨道不平顺数值模拟方法有白噪声滤波法、二次滤波法、三角级数法以及逆傅氏变换法。本章针对各种模拟方法在理论上作了比较，最后从时间域谱密度函数出发，对逆傅氏变换模拟方法进行了探讨，给出了相应的时域模拟方法的原理及步骤。

5.1 国内外常用的数值模拟方法

国内外常用的轨道不平顺数值模拟方法有白噪声滤波法、二次滤波法、三角级

[6]

数法以及逆傅氏逆变换法，下面将介绍前三种方法，逆傅氏逆变换法将单独列为一节详细介绍。

5.1.1 白噪声滤波法及二次滤波法

白噪声滤波法将轨道不平顺随机过程抽象为满足一定条件的白噪声，并将其视为一个线性系统的输入，而输出就是满足特定功率谱要求的随机过程样本，这也是地震波人工模拟的一种常用方法；二次滤波法针对不同的轨道谱设计不同的滤波器，把白噪声经过两次滤波拟合为满足特定谱的随机样本。这两种方法的样本拟合精度都和滤波器的设计有关。为了提高样本的模拟质量，文献[9]在滤波器设计中采用了遗传算法。时域样本的二次滤波模拟方法，具体步骤如下：

（1）模拟产生一个均匀分布于（0，1）的随机数列ui（可以采用累乘法程序生成）；

（2）将ui转换成均值为0、均方差为的正态分布随机数列xi,即

xilnuicos(2ui1)xi1lnuisin(2ui1) (5-1)

式中: i1,3,5,n1，22kAvvt，t为时间采样间隔。

（3）通过下面的滤波过程,可以通过随机变量xi产生随机过程样本函数i(i的谱密度即美国六级垂向时域不平顺谱)

i(a1)i1ai2(1a)xi (5-2)

式中，aexp(ct)。

上述模拟方法产生的轨道不平顺随机样本的频率很宽。为了满足轨道频谱的实际测量范围,需要对样本进行带宽滤波,这样一个满足要求的样本就形成了。

5.1.2 三角级数法

设轨道不平顺为各态历经随机过程x(t)，则其功率谱密度函数为

TFx(,T)x(t)eitdt (5-3)

012Fx(,T) (5-4) TT轨道不平顺随机过程可以采用复数形式的傅立叶级数表达

Sx()limx(t)Ckeik0t (5-5)

k0N11Ckx(t)eik0tdt (5-6)

T0其离散表达式为

2nkNTx(nt)Ckek0N1i (5-7)

2nki1N1Ckx(nt)eN (5-8)

Nn0由式(5-3)、(5-4)与式(5-5)、(5-6)得

11Ckx(t)eik0tdtFx(k0,T) (5-9)

T0T将式(5-7)、(5-8)带入式(5-3)、(5-4)得 122Sx(k0)limFx(k0,T)limTCk (5-10)

TTTTCklimTSx(k0) (5-11) T三角级数法是工程上一种常用方法，它假设轨道不平顺是具有零均值的各态历经平稳高斯过程，并将其用三角级数表示。下面对三角级数法从原理及步骤上给予阐述。

将轨道不平顺展开成复数傅立叶级数

x(t)1CjTjkT2Ceijjk0t (5-12)

x(t)eij0tdtT2由于CjCj且C00，故式(5-12)可改写为

x(t)2Ckcos(j0tj)j1kCj1TT2 (5-13)

x(t)eij0tdtT2为描述系数Cj的随机性质，引入随机变量j并令其在0~2区间上服从均匀分布，于是可构造出如下轨道不平顺随机过程的三角级数表达式

kx(nt)j12Sx(j)ktcos(jn) (5-14)

jktk式中，0/kt。

由于三角级数法同功率谱的快速傅式变换数值计算方法间存在区别，从而三角级数模拟样本的功率谱与解析功率谱间存在误差。

5.2 逆傅氏变换法

逆傅氏变换法基于功率谱快速数值算法原理，但在计算方法上采用了傅氏逆变

[1]

换，使随机样本的模拟速度较三角级数法有大幅度的提高。二次滤波法需要进行滤波器的设计，对不同功率谱密度函数的轨道不平顺，均需设计出合理的滤波器，因此缺乏通用性。三角级数法和白噪声滤波法是将轨道不平顺看作平稳高斯随机过程，这显然不完全与实际情况相符，文献[10]证明轨道不平顺并不都是平稳随机过程，事实上，轨道不平顺功率谱是对时域采样信号通过周期图法估计而获得的，其计算核心是快速傅立叶变换，鉴于以上方法的各种缺点，而逆傅氏变换法误差较小，运算速度较快，通用性强，故本文采用逆傅氏变换法对功率谱函数进行数值模拟。

5.2.1 估计功率谱的Blackman-Turkey（BT）法

在介绍逆傅氏变换法之前，有必要先介绍一下估计功率谱的Blackman-Turkey法。设时间序列{xs}，s0,1,,(N1)，记录长度TN，为时间间隔，对于相关函数的时迟r，也为离散值，则

1Rxx()x(t)x(t)dt

T01N1RrRxx(r)xsxsr (5-15)

Ns0T式中：r0,1,,(N1)；xsx(s) ；Rxx、Rr为时间序列{xs}的自相关函数。则

r111N1ik2N (5-16) Sxx(f)Rxx()Sxx(k)Sxx(fk)RreTTNr0式中：Sxx(k)为时间序列{xs}的功率谱密度函数。所以

222N1N1ikriksik(sr)1N11N111Sxx(k){xsxsr}eN{xseN}{xsreN} (5-17)

Nr0Ns0Ns0Nr0令jrs，则

xr0N1sreik2(sr)NN1sjsxjeik2jN (5-18)

又因为对离散傅氏变换，时间序列{xs}已离散周期化，周期为N。所以

N1sjsxjeik2jNxjej0N1ik2jN (5-19)

将式(5-18)和式(5-19)带入式(5-17)得

1Sxx(k){Nxess0N1ik2sN1}{Nxejj0N1ik2jN} (5-20)

112*DFT[x][X(k)X(k)]sN2N2式中：X(k)为时间序列{xs}，s0,1,,(N1)的频谱；k0,1,,(N1)。

5.2.2 逆傅立叶变换的计算步骤

已有的轨道不平顺功率谱均是空间域谱，这给轨道不平顺样本的数值模拟及时域

评价带来不便，所以应先将空间频谱转化成时域频谱，再由上述估计功率谱的周期图法可知，功率谱密度值在离散的采样点上与信号的频谱有着一个确定的关系。如果能够在功率谱密度函数上离散采样，构造出频谱X(k)，然后再对其进行傅立叶逆变换，即可得到时域的模拟轨道不平顺激扰函数x(t)。实施的步骤如下：

（1）已有的轨道不平顺功率谱均是空间域谱，先将空间频谱化成时域谱，以美国六级线路轨道高低不平顺功率谱为例，利用4.3节提到的单位换算方法可以将其转化为时域空间的功率谱。

（2）得到时域功率谱后，由于其为单边谱，所以要将其转化为双边谱Sx(f)，如图5-1所示。设需模拟的时间序列总时间为Ts，时间间隔，则时域和频域采样点数为NrTs/，一般Nr不是正好为2的整数次幂，这时必须将其转化为2的整数次幂以保证快速傅立叶的计算速度足够快。频域采样间隔为f1/(Nr)。众所周知，由周期图法估计出的功率谱具有周期性，且为偶对称序列。如图5-1所示，设功率谱有效频率段上截止频率为fu和下截止频率fl，则有效频率段内的采样点数为

(01的)采样点值记为0。若Nf(fufl)/f。又设N0fl/f，则0~NN0NfNr/2，则Nf~Nr/2的采样点值记为0；若N0NfNr/2，则可增大Ts以满足N0NfNr/2。于是得到功率谱Sx(f)的Nr/2个离散采样点值Sx(fkf)，

k0,1,,Nr/2。最后再以此形成以Nr/2为对称中心的偶对称序列Sx(fkf)，k0,1,,Nr/2。

图5-1 周期功率谱密度采样图

（3）由式(5-16)和式(5-20)可知，时域序列的频谱模值为

X(k)DFT[x(n)]Nr2Sk(k)NrSk(k)NrSk(fkf)f (k0,1,,Nr1) (5-21)

（4）式(5-21)给定了序列x(n)的频谱X(k)的模值。由于时间序列x(n)为一随机过程，其频谱相位一定具有随机性。

设n为独立相位序列，它的各分量均值为零。又由于实序列的FFT为复序列（实部偶对称，虚部奇对称）。所以n应为复数，且有n1，故设

ncosnisinnexp(in) (5-22)

式中n服从0~2均匀分布。

又因为X(k)的实部关于Nr/2偶对称，虚部关于Nr/2奇对称，所以只需求出频谱X(k)，(k0,1,,Nr/2)，即由式(5-21)可得

X(k)(k)X(k)Nr(k)Sx(fkf)f (k0,1,,Nr/2) (5-23)

显然由对称条件容易得到X(k)，(k0,1,,Nr1)。

（5）将得到的复序列X(k)进行IFFT（快速傅立叶逆变换）可得

1x(n)NNr1k0X(k)ei2knNr (n0,1,,Nr1) (5-24)

5.2.3 算例

以美国六级线路轨道高低不平顺功率谱为例，其数学表达式如式(3-3)示，其中：

k0.25；V100 km/h27.78 m/s；

fccV0.824527.783.647 Hz；22AV0.0339 cm2rad/m，取空间波长为0.5~50 m。则

fuV/0.555.56 Hz flV/500.5556 Hz Ts10 s 0.0001 s 根据以上的计算步骤，利用MATLAB编程工具编写程序，可以得出数值模拟图线如下图示：（具体程序见附录 ）

模拟的时间序列0.80.6高低不平顺幅值/cm0.40.20-0.2-0.4-0.6012345678910时间/s 图5-1 模拟轨道不平顺的时间序列

100不平顺解析值与模拟值的比较 10-1模拟值解析值10-2功率谱密度/(cm2/Hz)10-310-410-510-610-710 -110-8100101102频率/Hz 图5-2 轨道不平顺功率谱解析值与模拟值的比较

由图5-1可以看出，在1~10s之间，时间序列呈十分随机的分布；模拟的时间序列幅值在-0.6~0.6cm之间，与文献[1]的结果相差不大，从而验证了模拟时间序列的正确性。

图5-2中解析值与模拟值的比较可以看出，二者基本重合，模拟精度较高，结果令人满意，这也从另一方面再次证明了模拟时间序列的正确性，这样的到的时间序列就可以用于动力学计算了。

结论与展望

轨道是铁路系统的主要技术部件，用于引导列车的运行，并将所承受的荷载传递给路基。轨道结构由钢轨、轨枕、联结零件、道床、防爬设备和道岔等主要部件构成，是一种非常复杂的组合体，具有组合性和散体性，并且所承受的荷载具有随机性和重复性，外加环境、养护维修等诸多因素的影响，致使轨道系统在其运营寿命中自始至终都存在不平顺。轨道的平顺状态，对行车的安全性、平稳性、舒适性，对车辆和轨道部件的寿命以及环境噪声等都有重要影响，是轨道方面直接限制提高行车速度的主要因素。因此，深入研究轨道平顺状态对于保证铁路安全运营具有十分重要的指导意义。本文从轨道谱和轨道不平顺的时频分析等方面对轨道不平顺进行了研究，主要研究成果和结论如下：

1、通过本次论文的撰写，进一步理解了轨道不平顺的概念及其形成原因，为以后深入研究轨道不平顺打下了坚实的基础。

2、在理解轨道不平顺概念的基础上，分析对比了国内外几种典型的轨道谱，通过分析比较，认识到我国在研究轨道谱方面的欠缺及我国轨道质量状况较差这一现实，从而明确了以后研究工作的方向和目标。

3、实测的轨道不平顺数据都是时域数据，通过学习研究，学会利用各种谱估计方法（本文中主要用到了周期图法、自相关法、改进的周期图法（Bartlett法和Welch法）和最大熵法）将其转化为轨道功率密度谱来分析评价线路的平顺性，从而在实现了由时域数据向频域数据的转化。

4、为满足动力学分析的需求，还要将频域数据转化为时域数据，这一过程就叫做数值模拟。在总结国内外常用方法的基础上，学会了利用逆傅氏变换法完成对轨道谱函数的数值模拟。

5、以上工作都是在MATLAB编程条件下完成的，因此通过此次对轨道不平顺的研究，熟练掌握了用MATLAB编程来完成数值计算的方法。

由于轨道系统自身和所承受荷载的特点，决定了轨道不平顺的研究是一个庞大的系统工程，需要大量数据资料的积累，涉及众多学科的综合应用，需要大量人力物力的投入，兼之作者水平、实践经验和研究时间的限制，因此，在已做工作的基础上，仍应在以下几个方面进行更深入细致的研究：

1、由于我国在轨道谱研究方面十分欠缺，只做出了少数线路的轨道谱，所以

在以后的研究中，应以建立我国通用轨道谱为目标，进一步完善轨道不平顺的研究工作。

2、在此次的研究中，主要应用经典谱估计方法进行了轨道功率谱的估计，而经典方法存在分辨率低、精度差等诸多弊端，因此有必要学习更多的现代谱估计方法来估计轨道的功率谱，以使估计的分辨率和精度提高，最终达到一致估计的目标。

3、对轨道不平顺的研究要为机车车辆的运行服务，因此在对轨道不平顺研究的基础上，还要深入到车辆动力学领域，研究在轨道不平顺信号激扰下车辆的动力状况，即车辆-轨道耦合动力学。

致 谢

在本文的研究和撰写过程中，得到了导师翟婉明教授的悉心指导、教诲、鼓励和关怀，翟老师在本文的选题、研究内容和研究方法上给予了具体的指导和严格的要求，翟老师在繁忙的科研工作中抽出时间，为学生指明方向，指点迷津，提点不足，一直给与了学生极大的鞭策和帮助，在此谨致以深深的谢意和真挚的敬意。

在论文的撰写过程中，得到了赵春发老师的无私帮助，在一些研究难题上提出了许多宝贵意见，学生在此谨致以诚挚的感谢。

在论文撰写过程中，得到了高建敏博士、王少林博士、徐鹏博士的无私帮助，在研究内容和研究方法上提出了许多宝贵的指导意见，同时还在百忙中抽出时间仔细审阅了论文，学生在此谨致以衷心的感谢。

对以上诸位和其它对本文工作给予了关心帮助的所有人员一并致以诚挚的谢意！

特别感谢含辛茹苦的父母，你们用浸满艰辛与汗水的默默支持，铺就了儿子的求学之路，无以为报，只期盼在我取得一点点进步之时，我深爱的家人脸上能够展露出一丝的欣慰笑容！

感谢其他亲朋好友给予我的帮助和支持！

最后，真诚的感谢参加本文评阅和答辩的专家和学者！

参考文献

[1] 陈果，翟婉明.铁路轨道不平顺随机过程的数值模拟[J].西南交通大学学报.1999，34(2)：

138~142

[2] 翟婉明.车辆-轨道耦合动力学[M].第三版.北京：科学出版社，2007，106~114 [3] 张玲华，郑宝玉.随机信号处理[M].北京：清华大学出版社.2003

[4] 王元丰，王颖，王东军.铁路轨道不平顺模拟的一种新方法[J].铁道学报.1997，19(6)：

110~115

[5] 柯在田.轨道不平顺的数值模拟[J].铁道部科学研究院.1990

[6] 张立伟，冯军和.轨道不平顺随机过程的数值模拟[J].世界地震工程.2008，24(1)：103~109 [7] 张宪麦.轨道不平顺时频域分析及预测方法的研究[D].北京：铁道科学研究院博士学位论

文.2006，1~19

[8] 练松良.轨道工程[M].上海：同济大学出版社，2006

[9] 钱雪军.轨道不平顺的时域模拟法.铁道学报[J].2000，22(4)：94~98 [10] 罗林.轨道随机干扰函数[J].中国铁道科学.1982，3(1)：74~111

[11] 姚武川，姚天任.经典谱估计方法的MATLAB分析[J].华中理工大学学报.2000，28(4)：

45~47

[12] 常巍，谢光军，黄朝峰.MATLABR2007基础与提高[M].北京：电子工业出版社.2007 [13] 肖先赐.现代谱估计-原理与应用[M].哈尔滨：哈尔滨工业大学出版社.1991 [14] Michael J. Roberts.信号与系统[M].胡剑凌译.北京：机械工业出版社.2006 [15] 傅广操，樊明捷.MATLAB在现代功率谱估计中的应用[J].电脑学习.2003

[16] 张洁.轨道不平顺非均匀采样信号谱分析[D].成都：西南交通大学博士学位论文.2004 [17] 陈桂林，张照明，戚红雨.应用MATLAB语言处理数字信号与数字图像[M].北京：科学

出版社.2000

[18] 郭仕剑，王宝顺，贺志国，杨可心.MATLAB7.X数字信号处理[M].北京：人民邮电出版

社，2006.

[19] 万永革.数字信号处理的MATLAB实现[M].北京：科学出版社.2007 [20] 童大陨.铁路轨道及路基[M].北京：中国铁道出版社.1982

[21] 魏鑫，张平.周期图法功率谱估计中的窗函数分析[J].测控技术.2004，25(3)：14~15

附录1

轨道不平顺分析的MATLAB程序

1.1 轨道不平顺功率谱分析的MATLAB程序

1.1.1 美国六级谱

function [q1,q2,q3]=usanew() %函数名 fs=1000; %采样频率 l=0.5:1/fs:300; %波长范围 Av=0.2095; Aa=0.0762; Wc=0.8245; Ws=0.8209; k=0.25; Lc=2*pi/Wc;

Ls=2*pi/Ws; %参数设置 %美国五级线路

S1=100*k*Av./(2*pi*(2*pi./(Wc*l)).^2+2*pi); %高低不平顺 S2=100*k*Aa./(2*pi*(2*pi./(Wc*l)).^2+2*pi); %方向不平顺

S3=100*2*k*Av./(pi*l.^2*Lc.^2.*(1./l.^2+1/Lc.^2).*(1./l.^2+1/Ls.^2)); %水平及轨距不平顺 %美国六级线路 Av2=0.0339; Ws2=0.4380;

Ls2=2*pi/Ws2; %参数设置

S4=100*k*Av2./(2*pi*(2*pi./(Wc*l)).^2+2*pi); %高低不平顺 S5=100*k*Av2./(2*pi*(2*pi./(Wc*l)).^2+2*pi); %方向不平顺

S6=100*2*k*Av2./(pi*l.^2*Lc.^2.*(1./l.^2+1/Lc.^2).*(1./l.^2+1/Ls2.^2)); %水平及轨距不平顺 figure;

loglog(l,l.^2.*S4,'r'); %绘制波长与功率谱密度的双对数坐标图

title('美国六级高低不平顺','fontsize',25); %题名，并设置字号为25 ylabel('功率谱密度/[mm2/(rad/m)]','fontsize',20); %y轴名称 grid on %加网格

xlabel('波长/m','fontsize',20); %x轴名称 figure

loglog(l,l.^2.*S5,'r');

title('美国六级方向不平顺','fontsize',25); ylabel('功率谱密度/[mm2/(rad/m)]','fontsize',20); grid on

xlabel('波长/m','fontsize',20); figure

loglog(l,l.^2.*S6,'r');

title('美国六级水平及轨距不平顺','fontsize',25); ylabel('功率谱密度/[mm2/(rad/m)]','fontsize',20); grid on

xlabel('波长/m','fontsize',20);

1.1.2 德国低干扰轨道谱

function [q1,q2,q3]=germannew() %函数名 fs=1000; %采样频率 l=0.5:1/fs:300; %波长范围 %参数设置

Av=4.032*10^(-7); Wc=0.8246; Wr=0.0206; Ws=0.4380; Aa=2.119*10^(-7); Lc=2*pi/Wc; Lr=2*pi/Wr; Ls=2*pi/Ws;

b=0.75;

S1=Av./(2*pi*l.^2*Lc^2.*(1./l.^2+1/Lr^2).*(1./l.^2+1/Lc^2)); %高低不平顺

S2=Aa./(2*pi*l.^2*Lc^2.*(1./l.^2+1/Lr^2).*(1./l.^2+1/Lc^2)); %水平不平顺

S3=Av./(2*pi*b^2*l.^4*Lc^2.*(1./l.^2+1/Lr^2).*(1./l.^2+1/Lc^2).*(1./l.^2+1/Ls^2)); %轨距不平顺 figure; %开新窗口

q1=loglog(l,S1.*l.^2*10^6); %绘制波长与功率谱密度的双对数坐标图 ylabel('功率谱密度/[mm2/(rad/m)]','fontsize',20); %y坐标名称 grid on %加网格

xlabel('波长/m','fontsize',20); %x坐标名称

title('德国低干扰高低不平顺','fontsize',25); %题名，并设置字号为25 figure;

q2=loglog(l,S2.*l.^2*10^6);

ylabel('功率谱密度/[mm2/(rad/m)]','fontsize',20); grid on

xlabel('波长/m','fontsize',20);

title('德国低干扰方向不平顺','fontsize',25); figure;

q3=loglog(l,S3.*l.^4*10^6);

ylabel('功率谱密度/[mm2/(rad/m)]','fontsize',20); grid on

xlabel('波长/m','fontsize',20);

title('德国低干扰水平不平顺','fontsize',25);

1.1.3 中国三大干线轨道谱

function [q]=china() %函数名称

t=fopen('BJ.txt','r'); %打开中国三大干线参数文件 p=fscanf(t,'%e'); %读取数据

fs=1000; %设置采样频率 l=1:1/fs:30; %波长范围

A=p(1,1);B=p(2,1);C=p(3,1);D=p(4,1);E=p(5,1);F=p(6,1);G=p(7,1); %给各参数赋值

S=A.*(l.^2+B*l.^3+C*l.^4)./(l.^2.*(G*l.^4+F*l.^3+E*l.^2+D*l+1)); %经单位转化后的轨道谱公式

q=loglog(l,l.^2.*S); %画波长与功率谱密度的双对数坐标图 title('三大干线左轨高低','fontsize',25) %题名，并设置字号为25 xlabel('波长/m','fontsize',20) %x轴名称 grid on %加网格

ylabel('功率谱密度/[mm2/(rad/m)]','fontsize',20) %y坐标名称 figure

A=p(8,1);B=p(9,1);C=p(10,1);D=p(11,1);E=p(12,1);F=p(13,1);G=p(14,1); %给各参数赋值

S=A.*(l.^2+B*l.^3+C*l.^4)./(l.^2.*(G*l.^4+F*l.^3+E*l.^2+D*l+1)); loglog(l,l.^2.*S);

title('三大干线左轨轨向','fontsize',25) xlabel('波长/m','fontsize',20) grid on

ylabel('功率谱密度/[mm2/(rad/m)]','fontsize',20) figure

A=p(15,1);B=p(16,1);C=p(17,1);D=p(18,1);E=p(19,1);F=p(20,1);G=p(21,1); %给各参数赋值

S=A.*(l.^2+B*l.^3+C*l.^4)./(l.^2.*(G*l.^4+F*l.^3+E*l.^2+D*l+1)); loglog(l,l.^2.*S);

title('三大干线水平','fontsize',25) xlabel('波长/m','fontsize',20) grid on

ylabel('功率谱密度/[mm2/(rad/m)]','fontsize',20)

1.1.4 中国郑武高速试验段轨道谱

function [q]=china2()

t=fopen('ZW.txt','r'); %打开参数文件 p=fscanf(t,'%e') fs=1000; %采样频率 l=1:1/fs:30; %波长范围

A=p(1,1);B=p(2,1);C=p(3,1);D=p(4,1);E=p(5,1);F=p(6,1);G=p(7,1); %参数赋值

S=A.*(l.^2+B*l.^3+C*l.^4)./(l.^2.*(G*l.^4+F*l.^3+E*l.^2+D*l+1)); loglog(l,l.^2.*S); %画波长与功率谱密度的双对数坐标图 title('郑武线高速试验段左轨高低不平顺','fontsize',25) xlabel('波长/m','fontsize',20) grid on

ylabel('功率谱密度/[mm2/(rad/m)]','fontsize',20) figure

A=p(8,1);B=p(9,1);C=p(10,1);D=p(11,1);E=p(12,1);F=p(13,1);G=p(14,1); S=A.*(l.^2+B*l.^3+C*l.^4)./(l.^2.*(G*l.^4+F*l.^3+E*l.^2+D*l+1)); q=loglog(l,l.^2.*S);

title('郑武线高速试验段左轨轨向不平顺','fontsize',25) xlabel('波长/m','fontsize',20) grid on

ylabel('功率谱密度/[mm2/(rad/m)]','fontsize',20) figure

A=p(15,1);B=p(16,1);C=p(17,1);D=p(18,1);E=p(19,1);F=p(20,1);G=p(21,1); S=A.*(l.^2+B*l.^3+C*l.^4)./(l.^2.*(G*l.^4+F*l.^3+E*l.^2+D*l+1)); loglog(l,l.^2.*S);

title('郑武线高速试验段水平不平顺','fontsize',25) xlabel('波长/m','fontsize',20) grid on

ylabel('功率谱密度/[mm2/(rad/m)]','fontsize',20)

1.1.5 中美轨道谱比较

function [q1,q2,q3]=usachina() %函数名 fs=1000; %采样频率 l=1:1/fs:30; %波长范围 %参数设置 Av=0.2095; Aa=0.0762; Wc=0.8245; k=0.25; %美国五级线路

S1=100*k*Av./(2*pi*(2*pi./(Wc*l)).^2+2*pi); %高低不平顺 S2=100*k*Aa./(2*pi*(2*pi./(Wc*l)).^2+2*pi); %方向不平顺 %美国六级线路 Av2=0.0339;

S4=100*k*Av2./(2*pi*(2*pi./(Wc*l)).^2+2*pi); %高低不平顺 S5=100*k*Av2./(2*pi*(2*pi./(Wc*l)).^2+2*pi); %方向不平顺 %高低不平顺对比 figure;

loglog(l,l.^2.*S1,'r:'); %绘制波长与功率谱密度的双对数坐标图（五级） title('普通线路高低不平顺','fontsize',25); %题名，并设置字号为25 ylabel('功率谱密度/[mm2/(rad/m)]','fontsize',20); %y轴名称 grid on

xlabel('波长/m','fontsize',20); %x轴名称 hold on %图形保持

loglog(l,l.^2.*S4,'r-.'); %六级 hold on

%中国三大干线谱 t=fopen('BJ.txt','r');

p=fscanf(t,'%e');

A=p(1,1);B=p(2,1);C=p(3,1);D=p(4,1);E=p(5,1);F=p(6,1);G=p(7,1); S=A.*(l.^2+B*l.^3+C*l.^4)./(l.^2.*(G*l.^4+F*l.^3+E*l.^2+D*l+1)); loglog(l,l.^2.*S);

legend('美国五级线路谱','美国六级线路谱','中国三大干线谱') %加图示 %方向不平顺对比 figure

q2=loglog(l,l.^2.*S2,'r:');

title('普通线路方向不平顺','fontsize',25); ylabel('功率谱密度/[mm2/(rad/m)]','fontsize',20); grid on

xlabel('波长/m','fontsize',20); grid on hold on

loglog(l,l.^2.*S5,'r-.'); hold on

A=p(8,1);B=p(9,1);C=p(10,1);D=p(11,1);E=p(12,1);F=p(13,1);G=p(14,1); S=A.*(l.^2+B*l.^3+C*l.^4)./(l.^2.*(G*l.^4+F*l.^3+E*l.^2+D*l+1)); loglog(l,l.^2.*S);

legend('美国五级线路谱','美国六级线路谱','中国三大干线谱')

1.1.6 中德轨道谱比较

function [q1,q2]=germanchinanew() fs=1000; %采样频率 l=1:1/fs:100; %波长范围 %参数设置

Av=4.032*10^(-7); Wc=0.8246; Wr=0.0206; Aa=2.119*10^(-7);

Lc=2*pi/Wc; Lr=2*pi/Wr; %德国低干扰轨道谱

S1=Av./(2*pi*l.^2*Lc^2.*(1./l.^2+1/Lr^2).*(1./l.^2+1/Lc^2)); %高低不平顺

S2=Aa./(2*pi*l.^2*Lc^2.*(1./l.^2+1/Lr^2).*(1./l.^2+1/Lc^2)); %方向不平顺

%德国高干扰轨道谱 Av2=10.8*10^(-7); Aa2=6.125*10^(-7);

S4=Av2./(2*pi*l.^2*Lc^2.*(1./l.^2+1/Lr^2).*(1./l.^2+1/Lc^2)); %高低不平顺

S5=Aa2./(2*pi*l.^2*Lc^2.*(1./l.^2+1/Lr^2).*(1./l.^2+1/Lc^2)); %方向不平顺 figure;

q1=loglog(l,S1.*l.^2*10^6,'r:'); %绘制波长与功率谱密度的双对数坐标图（低干扰）

ylabel('功率谱密度/[mm2/(rad/m)]','fontsize',20); grid on

xlabel('波长/m','fontsize',20);

title('高速线路高低不平顺','fontsize',25); hold on

loglog(l,S4.*l.^2*10^6,'r-.'); %绘制波长与功率谱密度的双对数坐标图（高干扰） hold on

%中国郑武高速试验段轨道谱 t=fopen('ZW.txt','r'); p=fscanf(t,'%e'); l=1:1/fs:30;

A=p(1,1);B=p(2,1);C=p(3,1);D=p(4,1);E=p(5,1);F=p(6,1);G=p(7,1); S=A.*(l.^2+B*l.^3+C*l.^4)./(l.^2.*(G*l.^4+F*l.^3+E*l.^2+D*l+1)); loglog(l,l.^2.*S);

legend(' 德国低干扰轨道谱','德国高干扰轨道谱','中国高速试验段谱') %加图示 figure; l=1:1/fs:100;

q2=loglog(l,S2.*l.^2*10^6,'r:');

ylabel('功率谱密度/[mm2/(rad/m)]','fontsize',20); grid on

xlabel('波长/m','fontsize',20);

title('高速线路方向不平顺','fontsize',25); hold on

loglog(l,S5.*l.^2*10^6,'r-.'); hold on l=1:1/fs:30;

A=p(8,1);B=p(9,1);C=p(10,1);D=p(11,1);E=p(12,1);F=p(13,1);G=p(14,1); S=A.*(l.^2+B*l.^3+C*l.^4)./(l.^2.*(G*l.^4+F*l.^3+E*l.^2+D*l+1)); loglog(l,l.^2.*S);

legend(' 德国低干扰轨道谱','德国高干扰轨道谱','中国高速试验段谱')

1.2 实测数据轨道谱估计MATLAB程序

1.2.1 周期图法的功率谱图(直接法)

function [w]=tt() %函数名称 l=1; %波长下限 v=150/3.6; %车速 fs=2*v/l; %采样频率

f=fopen('0412.txt','r'); %打开04年12月轨道高低不平顺实测数据 y=fscanf(f,'%e'); %读入原始数据

nfft=2^15; %FFT的点数

Xk=fft(y,nfft); %快速傅立叶变换

%save fft.txt Xk -ascii; %存储傅立叶变换值 p=abs(Xk).^2/length(Xk); %计算功率谱 %save psd.txt p -ascii; %存储功率谱值 index=0:round(nfft/2-1); %时间坐标取一半 k=index.*fs/nfft; %对应的频率 figure

w=loglog(k,p(index+1)); %画频率和功率谱密度的双对数坐标图 title('周期图法谱估计','fontsize',25) %标题名,字号为25 xlabel('频率(Hz)','fontsize',20); %x轴名称 grid on

ylabel('功率谱密度(mm2)','fontsize',20); %y轴名称

1.2.2 自相关法(间接法)

function [w]=indirect() %函数名 l=1; %波长下限 v=150/3.6; %车速 fs=2*v/l; %采样频率

f=fopen('0412.txt','r'); %打开04年12月轨道高低不平顺实测数据 y=fscanf(f,'%e'); %读取原始数据

nfft=2^15; %FFT的点数

R=xcorr(y,'unbiased'); %自相关函数 F=fft(R,nfft); %快速傅立叶变换 p=abs(F); %求解功率谱密度

index=0:round(nfft/2-1); %时间坐标取一半 k=index*fs/nfft; %对应的频率 figure %新窗口

w=loglog(k,p(index+1)); %画频率和功率谱密度的双对数坐标图

title('自相关法','fontsize',25) %标题名,字号为25 xlabel('频率(Hz)','fontsize',20); %x轴名称 grid on %加网格

ylabel('功率谱密度(mm2)','fontsize',20); %y轴名称

1.2.3 改进的直接法估计(bartlett法)

function [t]=Bartlett() l=1; %波长下限 v=150/3.6; %车速 fs=2*v/l; %采样频率

f=fopen('0412.txt','r'); %打开04年12月轨道高低不平顺实测数据 y=fscanf(f,'%e'); %读入原始数据

nfft=2^15; %FFT的点数 window=boxcar(2^8); %加矩形窗 noverlap=0; %不重叠 w=0.9; %置信度90%

p=psd(y,nfft,fs,window,noverlap,w); %求功率谱密度 index=0:round(nfft/2-1); k=index*fs/nfft; %确定频率 figure

t=loglog(k,p(index+1)); %画频率和功率谱密度的双对数坐标图 title('bartlett法谱估计','fontsize',25) %标题名,字号为25 xlabel('频率(Hz)','fontsize',20); %x轴名称 grid on

ylabel('功率谱密度(mm2)','fontsize',20); %y轴名称

1.2.4 改进的直接法估计(welch法)

function [t]=welch() %函数名 l=1; %波长下限 v=150/3.6; %车速

fs=2*v/l; %采样频率

f=fopen('0412.txt','r'); %打开04年12月轨道高低不平顺实测数据 y=fscanf(f,'%e'); %读入原始数据

nfft=2^15; %FFT的点数 window=hamming(2^8); %加汉明窗 noverlap=128; %重叠128点 w=0.9;%置信度90%

p=psd(y,nfft,fs,window,noverlap,w); %计算功率谱密度 index=0:round(nfft/2-1);

k=index*fs/nfft; %确定对应频率 figure

t=loglog(k,p(index+1)); %画频率和功率谱密度的双对数坐标图 title('welch法谱估计','fontsize',25) %标题名,字号为25 xlabel('频率(Hz)','fontsize',20); %x轴名称 grid on

ylabel('功率谱密度(mm2)','fontsize',20); %y轴名称

1.2.5 最大熵法的功率谱图

function [w]=MEM() %函数名称 l=1; %波长下限 v=150/3.6; %车速 fs=2*v/l; %采样频率

f=fopen('0412.txt','r'); %打开04年12月轨道高低不平顺实测数据 y=fscanf(f,'%e'); %读入原始数据

nfft=2^15; %FFT的点数

[p,f]=pmem(y,64,nfft,fs); %最大熵法计算功率谱密度 figure

w=loglog(f,p); %画频率和功率谱密度的双对数坐标图

title('最大熵法谱估计','fontsize',25) xlabel('频率(Hz)','fontsize',20); grid on

ylabel('功率谱密度(mm2)','fontsize',20);

1.2.6 实测轨道谱与美国六级谱高低不平顺对比

function [w]=ttyushice() %函数名称 l=0.5; %波长下限 v=100/3.6; %车速 fs=2*v/l; %采样频率 f=fopen('0412.txt','r'); y=fscanf(f,'%e'); %读入原始数据

nfft=2^15; %FFT的点数

Xk=fft(y,nfft); %快速傅立叶变换

%save fft.txt Xk -ascii; %存储傅立叶变换值 p=abs(Xk).^2/length(Xk); %计算功率谱 %save psd.txt p -ascii; %存储功率谱值 index=0:round(nfft/2-1); %时间坐标取一半 k=index.*fs/nfft; %对应的频率 figure

w=loglog(k,p(index+1)); %画频率和功率谱密度的双对数坐标图

title('实测高低不平顺轨道谱与美国六级高低不平顺谱的比较','fontsize',25) xlabel('频率(Hz)','fontsize',20); grid on

ylabel('功率谱密度(mm2)','fontsize',20); grid on hold on

v=100/3.6; %列车运行速度 delta=0.0001; %时间间隔

Nr=2^17; %时域和频域采样点数 df=1/(Nr*delta); %频域采样间隔 fu=v/0.5; %上截止频率 fl=v/50; %下截止频率

Nf=fix((fu-fl)/df); %有效频率段内的采样点数 f=linspace(fl,fu,Nf);

fc=82.45/(7.2*pi); %截断频率对应的时间频率

S=0.25*0.0339*(fc)^2*100/3.6./(2*pi*f.^2.*(f.^2+(fc)^2)); loglog(f,100*S.*f,'r--.'); %绘制解析谱线

legend('实测高低不平顺轨道谱线','美国六级高低不平顺谱线')

1.2.7 趋势比较分析(welch法)

function [t]=qs() %函数名 l=1; %波长下限 v=150/3.6; %车速 fs=2*v/l;%采样频率

f=fopen('0412.txt','r'); %打开04年12月轨道高低不平顺实测数据 y=fscanf(f,'%e'); %读入原始数据

nfft=2^15; %FFT的点数 window=hamming(2^8); %加汉明窗 noverlap=128; %重叠128点 w=0.9; %置信度90%

p=psd(y,nfft,fs,window,noverlap,w); %计算功率谱密度 index=0:round(nfft/2-1);

k=index*fs/nfft; %确定对应频率 figure

t=loglog(k,p(index+1)); %画频率和功率谱密度的双对数坐标图 title('趋势分析','fontsize',25) xlabel('频率(Hz)','fontsize',20);

ylabel('功率谱密度(mm2)','fontsize',20); %05年3月谱图 hold on %保持图形 f=fopen('0503.txt','r'); y=fscanf(f,'%e');

p=psd(y,nfft,fs,window,noverlap,w); loglog(k,p(index+1),'r-'); %06年1月谱图 hold on

f=fopen('0601.txt','r'); y=fscanf(f,'%e');

p=psd(y,nfft,fs,window,noverlap,w); loglog(k,p(index+1),'g-'); %07年1月谱图 hold on

f=fopen('0701.txt','r'); y=fscanf(f,'%e');

p=psd(y,nfft,fs,window,noverlap,w); loglog(k,p(index+1),'y-'); %04年9月谱图 hold on

f=fopen('0409.txt','r'); y=fscanf(f,'%e');

p=psd(y,nfft,fs,window,noverlap,w); loglog(k,p(index+1)); %05年12月谱图 hold on

f=fopen('0512.txt','r'); y=fscanf(f,'%e');

p=psd(y,nfft,fs,window,noverlap,w); loglog(k,p(index+1),'r-'); %05年3月谱图 hold on

f=fopen('0503.txt','r'); y=fscanf(f,'%e');

p=psd(y,nfft,fs,window,noverlap,w); loglog(k,p(index+1),'r-'); %06年1月谱图 hold on

f=fopen('0601.txt','r'); y=fscanf(f,'%e');

p=psd(y,nfft,fs,window,noverlap,w); loglog(k,p(index+1),'g-'); %06年2月谱图 hold on

f=fopen('0602.txt','r'); y=fscanf(f,'%e');

p=psd(y,nfft,fs,window,noverlap,w); loglog(k,p(index+1),'g-'); %06年3月谱图 hold on

f=fopen('0603.txt','r'); y=fscanf(f,'%e');

p=psd(y,nfft,fs,window,noverlap,w); loglog(k,p(index+1),'g-'); %06年4月谱图 hold on

f=fopen('0604.txt','r');

y=fscanf(f,'%e');

p=psd(y,nfft,fs,window,noverlap,w); loglog(k,p(index+1),'g-'); %06年6月谱图 hold on

f=fopen('0606.txt','r'); y=fscanf(f,'%e');

p=psd(y,nfft,fs,window,noverlap,w); loglog(k,p(index+1),'g-'); %06年9月谱图 hold on

f=fopen('0609.txt','r'); y=fscanf(f,'%e');

p=psd(y,nfft,fs,window,noverlap,w); loglog(k,p(index+1),'g-'); %06年12月谱图 hold on

f=fopen('0612.txt','r'); y=fscanf(f,'%e');

p=psd(y,nfft,fs,window,noverlap,w); loglog(k,p(index+1),'g-'); %07年1月谱图 hold on

f=fopen('0701.txt','r'); y=fscanf(f,'%e');

p=psd(y,nfft,fs,window,noverlap,w); loglog(k,p(index+1),'y-'); %07年3月谱图 hold on

f=fopen('0703.txt','r'); y=fscanf(f,'%e');

p=psd(y,nfft,fs,window,noverlap,w); loglog(k,p(index+1),'y-'); %07年5月谱图 hold on

f=fopen('0705.txt','r'); y=fscanf(f,'%e');

p=psd(y,nfft,fs,window,noverlap,w); loglog(k,p(index+1),'y-'); %07年6月谱图 hold on

f=fopen('0706.txt','r'); y=fscanf(f,'%e');

p=psd(y,nfft,fs,window,noverlap,w); loglog(k,p(index+1),'y-');

legend('04','05','06','07');%添加图例

1.3 数值模拟的MATLAB程序

function [Nf]=sm() v=100/3.6;%列车运行速度 delta=0.0001;%时间间隔 Nr=2^17;%时域和频域采样点数 df=1/(Nr*delta);%频域采样间隔 fu=v/0.5;%上截止频率 fl=v/50;%下截止频率

Nf=fix((fu-fl)/df);%有效频率段内的采样点数 No=round(fl/df);%下截止频率以下点数 t=zeros(1,No);%下截止频率以下点数置为零 m=zeros(1,Nr/2-Nf-No);%上截止频率到Nf置为零

f=linspace(fl,fu,Nf);%将上下截止频率之间等分Nf份 fc=82.45/(7.2*pi);%截断频率对应的时间频率

S=0.25*0.0339*(fc)^2*100/3.6./(2*pi*f.^2.*(f.^2+(fc)^2));%美国六级谱对应的公式

Sx=[t,S,m,0];%补全采样点 k=Nr/2+1;

dfai=rand(1,k)*2*pi;%角度在0~2pi间均匀分布 Y=exp(dfai*i);%独立相位序列 Xk1=Nr*Y.*sqrt(Sx*df);%频谱 Xk2=Xk1(2:k-1); Xk3=rot90(real(Xk2),2); Xk4=rot90(imag(Xk2),2); Xk5=Xk3-i*Xk4;

Xk=[Xk1,Xk5];%补全频谱

x=ifft(Xk,Nr);%傅立叶逆变换求时域样本 figure

plot(0.0001*(1:64:Nr),x(1:64:Nr));%画时域样本 title('模拟的时间序列','fontsize',25) xlabel('时间/s','fontsize',20)

ylabel('高低不平顺幅值/cm','fontsize',20) Sxn=(abs(Xk)/Nr).^2/df; Kn=No+1:No+Nf; f=Kn*df; figure;

loglog(f,Sxn(Kn),'b-.');%绘制由时域样本得到的功率谱曲线 hold on

f=linspace(fl,fu,Nf);

fc=82.45/(7.2*pi);%截断频率对应的时间频率

S=0.25*0.0339*(fc)^2*100/3.6./(2*pi*f.^2.*(f.^2+(fc)^2));

loglog(f,S);%绘制解析谱线 grid on

title('不平顺解析值与模拟值的比较','fontsize',25) xlabel('频率/Hz','fontsize',20)

ylabel('功率谱密度/(cm2/Hz)','fontsize',20) legend('模拟值','解析值')

附录2

毕业实习报告

实习时间：2008年5月7日至2008年5月9日

实习单位：东方电机厂、东方汽轮机厂

春暖花开，鸟语花香，在这样美丽的春季我们进行了本科阶段最后一次实习，实习的目的地是坐落在德阳的东方电机厂和东方汽轮机厂，这两大工厂实力十分雄厚，技术国际领先，在全国排名也是数一数二的。东方电机厂曾为三峡大坝制作过发电机组，东方汽轮机厂的设备也广泛分布在全国各个发电厂中。

德阳是一座十分漂亮的城市，驱车来到这里，一下子就感觉到小城的恬静与舒适。我们首先参观的是东方电机厂，在入厂之前，厂方派了一名我校的老校友给我们简介东方电机厂的辉煌历史，这位老校友讲的十分生动，谈笑间底气十足，他无时无刻不以作为一名东方电机厂职工而骄傲，听得我们对该厂也是心驰神往。听完这段精彩的介绍后，老校友带着我们前往厂区参观。一个企业的形象靠什么来反映，靠细节，进入厂房之前我就注意到了厂区内卫生十分整洁，上班时间基本没有几个员工在厂区游荡，当时我就感慨东方电机不愧是大厂，果然有大厂的风范。集体领过安全帽之后，我们分为4个组进入厂房，厂房内也十分干净，零件整齐堆放，工业垃圾也是集中处理的。我们十几个人被一个技术员带着在厂房内参观，每到一处，她都细致的讲解这个零件是如何制造的，用到了那些技术，还有这个零件的作用，虽然我们听得不太懂，但是我们也十分佩服她对专业知识掌握的如此熟练。介绍中也不乏关于企业文化的内容，她还无意间给我们讲了个小故事，讲的是在三峡大坝一台发电机组的竞标中，只有东方电机和通用公司入围，而经过多次分析比较，东方电机被选中，通用公司也声明以后再不涉足该发电机组的研究，这可是为中国人争了光啊！由此也可一窥东方电机技术力量的雄厚。技术员领着我们参观了主要的几个生产厂房，我们不仅了解了生产方面的一些技术问题，更重要的是了解了一个企业的文化，一个企业的精神，这才是企业长远发展的支撑力量。很快我们就结束了对东方电机厂的参观，我也期待着一睹东方汽轮机厂的风采。

东方汽轮机厂的总厂坐落在德阳绵竹县的汉旺镇上，也正因此汉旺镇被评为全国百强镇之一，可见一个大的企业对地方经济发展的巨大推动作用。汉旺镇不愧是全国百强镇之一，这里经济十分发达，环境十分优美，依山傍水，我们还打趣地说

老了以后在这里安享晚年呢！东方汽轮机厂坐落在紫岩山脚下，同属东方电气集团旗下，主要制造大型汽轮机，产品远销海内外。经过技术人员的简单介绍之后，我们赶往厂区进行参观。厂区内被分为好几个小的生产厂区，我们参观了几个主要的厂区。每个厂区都有技术人员给我们详细介绍，我们也不时提出问题，技术人员总是很耐心的回答，直到我们明白为止。厂区内卫生干净整洁，员工衣着统一，各个精神饱满，一派干劲十足的样子。各个厂区分工明确，生产效率很高，产品也是供不应求。花了整整一个上午我们结束了东汽之旅，这里先进的生产技术，优秀的企业文化，深深的打动着我，也算我见了一回大世面。

花了两天时间，我们结束了本次实习，带着不一样的心情回到了学校。这次实习不同于以往的实习，它让我们近距离的接触大型企业，接触它先进的生产技术，良好的管理机制和优秀的企业文化，使我们认识到了企业不仅靠技术，更靠管理和文化，我也在潜意识里形成了对管理和文化的重视。我们现在学的是技术，但我们还要学管理、学文化，这对以后我们的就业有长远的影响，我也深深感谢这次实习，给了我新颖的感受，给了我不一样的理念。