基于响应面法及改进粒子群算法的重载铁路大机捣固作业参数优化

Optimizing Tamping Parameters for Heavy-Haul Railway Based on Response Surface Method and Improved Particle Swarm Optimization

摘要

大型机械作业是改善有砟轨道线路几何形位的有效手段，科学合理的捣固参数是确保铁路安全运营的关键因素。当前捣固参数的选择主要基于道床恢复质量作为评价指标，而忽视了捣固对道砟受力特性的影响。为此，采用离散元与多体动力学耦合方法，构建有砟道床与捣固装置的耦合模型，并通过现场试验验证模型的可靠性。在此基础上，结合响应面法和改进粒子群算法，以道床支承刚度和道砟间接触力超过临界值的比率作为响应值，对捣固频率、捣固振幅及捣镐移速等关键参数进行优化。结果表明：最优捣固参数组合为捣固频率39 Hz、捣固振幅8 mm和捣镐移速1 m · s^－1，采用该优化参数后，在道床支承刚度保持基本不变的前提下，道砟间作用力超过临界值的比例明显降低，降幅达33.33%；捣固参数优化后道砟受力减小，不仅提升了作业质量，还可延长道砟的使用寿命。

Abstract

Large-scale mechanical operation plays a pivotal role in improving the geometric alignment of ballasted railway tracks, and precise tuning of tamping parameters are key factors ensuring operational safety. Currently, tamping parameters are selected mainly based on ballast bed restoration quality as the evaluation metric, often neglecting the influence of ballast force characteristics during tamping. This study introduces a coupled model that integrates ballasted track dynamics with tamping machinery by combining the discrete element method and multibody dynamics, and model reliability was validated through field tests. Building on this validated framework, tamping frequency, amplitude, and pick speed were optimized by combining response surface method and improved particle swarm algorithm, with ballastbed support stiffness and the ratio of the inter-ballast contact forces exceeding the critical threshold as response values. The results reveal that optimal performance is achieved at a tamping frequency of 39 Hz, an amplitude of 8 mm, and a tamping pick speed of 1 m · s^-1. Using these parameters, while ballast bed support stiffness remains basically unchanged, the proportion of inter-ballast contact forces exceeding the critical value decreases by 33.33%. The optimized parameters reduce detrimental forces on the ballast, thereby enhancing work quality and extending its service life.

Graphical abstract

关键词

重载铁路 / 大机捣固 / 响应面法 / 粒子群算法 / 捣固参数优化

Key words

Heavy-haul railway / Machine tamping / Response surface method / Particle swarm optimization / Tamping parameter optimization

列车荷载长期作用下，有砟道床结构易受损，导致轨道服役性能下降，并引发影响行车安全的轨道不平顺等问题^［1］。为确保有砟轨道的安全运营，大型机械（以下简称大机）捣固作业已成为改善道床不平顺的重要手段^［2］。因此，深入研究捣固参数对捣固效果的影响，特别是从宏观与细观多角度进行综合分析，具有重要的理论和实践意义。

在大机捣固作业中，捣固参数的变化对捣固效率和质量具有显著影响。国内外学者通过试验和仿真手段，深入研究了大机捣固对有砟道床宏观与细观特征的影响。Zhou等^［3］采用离散元法（Discrete Element Method，DEM）构建了道砟箱捣固模型，研究表明夹持力和捣镐振动频率对道床密实度具有显著影响。Aingaran等^［4］通过室内三轴试验模拟列车荷载及捣固作业过程，揭示了捣固过程中道砟颗粒间力的传递机制，并指出应力反转会削弱原有道床经碾压逐渐形成的竖向承载结构，进而引起道床垂向承载能力的降低。Liu等^［5］基于现场试验，探讨了捣固次数对道床纵横向阻力及支承刚度的影响规律。Zhang等^［6-7］建立了离散元与多体动力学（Multi Body Dynamics，MBD）耦合的捣固模型，分析了捣固作业对不同道床区域密实度的扰动效应，并基于该模型提出了捣镐夹持力的优化方法，研究表明当夹持力为7.8 kN时道床刚度可达到最大值41.43 kN · mm^－1。

另一方面，尽管捣固作业能够有效恢复轨道几何状态，但不合理的捣固参数可能导致道砟破碎率升高，进而威胁列车运行安全。Saussine等^［8］通过数值仿真分析了捣固频率对道床密实度的影响，提出了合理的捣固频率范围。Descantes等^［9］通过室内单枕捣固试验，研究了不同粒径道砟的破碎行为，分析了捣固过程中道砟颗粒的典型破碎模式。Wang等^［10］基于数值方法探讨了捣固次数对道床扰动及破坏的影响，提出最佳捣固次数应控制在4次以内。高亮等^［11］研究了捣固频率对捣镐冲击力和道床力学质量的影响，发现当频率处于30~40 Hz时，道砟受力最小且破碎率最低。Shi等^［12-13］综合考虑捣固作业对有砟道床的扰动破坏及养护效果，采用离散元-多体动力学耦合方法，分析了捣固参数对道床宏观与细观特征的影响，并提出当捣固深度设置为20~30 mm时养护效果最佳。

从上述研究可以看出，当前大型机械捣固参数的研究主要集中于单一因素的影响，而对多因素多目标共同作用的探讨相对较少。然而，在实际捣固作业中，各参数间的相互影响不容忽视。近年来，响应面法（Response Surface Method，RSM）因其能够综合考虑多因素相互作用，已广泛应用于铁路运输领域。占玉林等^［14］基于数值仿真结果，利用响应面法和粒子群算法（Particle Swarm Optimization，PSO）优化了大跨钢箱拱桥整体临时提升支架的结构参数。杨宏印等^［15］提出了一种结合贝叶斯法和响应面法的桥梁有限元模型修正方法，显著优化了模型计算结果的准确性。响应面法与数值仿真的结合不仅能够有效减小试验规模，还能充分考虑不同因素间的相互作用。在大型机械捣固参数优化方面，Chi等^［16］利用响应面法研究了道岔区域道床的最优捣固参数，提出最佳振动频率为42 Hz，最佳振动幅度为0.8°；Zhou等^［17］通过响应面法确定了捣固后最小孔隙率的大型机械捣固参数；Shi等^［18］结合响应面法与遗传算法（Genetic Algorithm，GA）提出了一种捣固参数优化方法，确定了最优捣固参数为振动频率42 Hz、振动幅度8.4 mm、挤压力10.47 kN。与仅考虑单一影响因素的方法相比，响应面法能够综合考虑更多因素的影响，从而获取更为全面的优化捣固参数。

由上述研究可知，现有优化指标多集中于道床阻力和支承刚度等因素，忽略了捣固参数对道砟受力的影响，从而制约了捣固质量的提升。为此，结合捣固作业的实际动态过程，基于离散元与多体动力学方法，构建有砟道床与捣固装置的耦合模型。采用Box-Behnken设计法开展数值仿真试验，通过分析捣固参数与道床状态间的响应关系，以道床支承刚度及道砟间作用力超过临界值的比率作为响应值构建响应面函数，并进一步利用改进粒子群算法（Improved Particle Swarm Optimization，IPSO）进行寻优计算。最终获得使道床支承刚度和道砟破碎概率达到最优解的捣固作业参数。该方法为提升捣固作业效率和质量提供理论支持，可为现场维修作业提供科学指导。

1 捣固装置-有砟道床耦合模型

1.1　模型建立

道砟复杂的颗粒形状对维持轨道稳定性具有重要作用，同时能够有效吸收轮轨相互作用产生的冲击与振动。科学合理地精确模拟道砟外形尺寸，对于复现道砟间复杂的接触关系至关重要。目前，道砟颗粒模型主要有3种：引入滚动阻力矩的球颗粒模型、镶嵌颗粒模型和黏结颗粒模型。相较于使用球形颗粒简化道砟外形，精细化的道砟颗粒模型能够更准确地复现道砟间的力学特性^［19］。基于此，本文采用镶嵌颗粒模型，道砟真实颗粒及颗粒模型如图1所示。

依据规范要求，本文选取图1所示6种不同尺寸的道砟，以重载铁路I级道砟级配为标准^［20］，通过落雨法构建符合图2所示级配曲线的有砟道床离散元模型。

捣固装置-有砟道床耦合模型由有砟道床模型及捣镐模型组成。在离散元软件中，采用墙体模拟Ⅲ型混凝土枕。其中，道砟-道砟、道砟-捣镐及道砟-轨枕间的接触均选用赫兹-明德林无滑移接触模型^［21-22］，该模型不考虑相互间的黏结效应。

捣固装置及轨枕属于典型的连续介质，且在捣固作业过程中涉及平动、转动及振动等一系列复杂机械运动。因此，本文采用MBD方法模拟轨枕及捣固装置。DEM-MBD耦合方法构建的大机捣固装置-有砟道床耦合模型如图3所示。

捣固装置-有砟道床耦合模型参数取值见表1。其中，道砟参数参考文献［23］，轨枕及捣固装置的参数依据文献［11，24］，并结合道床横向阻力测试结果确定。

数值仿真时参照现场实际捣固过程，将捣固作业分为起道、插镐、夹持和撤镐提升4个阶段。起道阶段，轨枕被大机抬升至指定高度；插镐阶段，捣镐竖直下插至指定深度，此阶段捣固装置与道砟发生强烈相互作用；夹持阶段，捣镐进行夹持运动，将轨枕侧部较深位置的道砟挤压至轨枕底部，以提高轨枕底部道床的密实度；撤镐提升阶段，捣镐逐渐恢复至初始位置。本文以单次捣固作业为分析对象，完成1个捣固循环内捣固参数的优化分析。

1.2　模型可靠性验证

现有研究^［25］表明，道床横向阻力能够有效反映道床状态。本文在构建捣固装置-有砟道床耦合模型后，选取捣固前后的道床横向阻力作为验证模型可靠性的指标。采用文献［26］所述的加载及数据采集方法开展试验，获取道床横向阻力曲线。

在现场试验中，利用某重载铁路捣固现场，测试捣固前后的道床横向阻力。数值仿真时，为减少边界效应的影响，采用推动中间轨枕的方式获取道床横向阻力。现场试验与数值仿真的道床横向阻力变化曲线如图4所示。

由图4可知：现场测试与数值仿真获得的道床横向阻力曲线具有相似的变化规律；道床横向阻力随轨枕位移的增加先增大，在位移达到2 mm后趋于平缓，道床整体力学响应呈现明显的弹塑性特征；无论是数值仿真还是现场试验，道床横向阻力在道床受力进入弹塑性阶段时都会出现波动现象，这是道砟间互锁状态在受到扰动时发生接触状态的改变导致的；从数值对比结果来看，捣固前数值仿真与试验结果在屈服阶段的横向阻力差值为0.47 kN，捣固后差值为0.34 kN，表明数值仿真与实测结果在趋势和数值上均具有较高的一致性；此外，重载铁路捣固后的道床横向阻力较捣固前显著降低，这是因为捣固后道床密实度下降，颗粒间接触变得松散，数值仿真与现场试验结果规律相一致。

综上所述，捣固作业前后的道床横向阻力试验结果与仿真结果规律一致，表明本文构建的捣固装置-有砟道床耦合模型具有良好的可靠性与准确性。相较于现场试验，数值仿真在捣固参数优化设计中更易于控制，并能获取更为全面的道床宏观与细观状态数据。

2 基于响应面法与改进粒子群算法优化方法

针对捣固参数的优化，在建立数值模型的基础上，提出了1种基于响应面法和改进粒子群算法的优化方法。该方法首先建立基于DEM-MBD的捣固数值模型，并确定捣固作业参数作为研究变量。在此基础上，采用Box-Behnken设计试验获取不同参数组合下的试验数据，并构建目标响应值的响应面函数，以描述参数与目标函数之间的关系。最后，引入改进粒子群算法对响应面函数进行优化求解，获得捣固作业的最优解，从而为捣固参数的优化配置提供理论依据与方法支持。计算流程如图5所示。

2.1　响应面法基本原理

响应面法是一种渐进近似的优化方法，常用于分析多个试验变量与指标间的关系。该方法通过一系列合理有效的试验，计算或采集所需试验数据，并基于这些数据构建响应面函数，以拟合设计变量与响应值之间的数学关系，最终通过求解该函数获得响应值的最优解。其过程包括筛选试验、试验设计、方程回归、响应面函数检验和优化求解。

响应面函数如式（1）所示。

y={\beta }_{\mathrm{0}}+{\displaystyle \sum _{i=\mathrm{1}}^k}{\beta }_i{x}_i+{\displaystyle \sum _{i=\mathrm{1}}^k}{\beta }_{ii}{x}_i^{\mathrm{2}}+{\displaystyle \sum _{i=\mathrm{1}}^{k-\mathrm{1}}}{\displaystyle \sum _{j=i+\mathrm{1}}^k}{\beta }_{ij}{x}_i{x}_j+\epsilon（1）

式中：y为响应面函数的值；{\beta }_{\mathrm{0}}为常数项；x_i和x_j为设计因子；{\beta }_i为因子x_i的1次项系数；{\beta }_{ii}为因子x_i的2次项系数；{\beta }_{ij}为因子x_i和因子x_j之间的线性交互效应系数；\epsilon为随机误差项；i和j为下标索引，范围为1—k，k为设计因子数目。

响应面函数拟合模型能否反映目标和变量的关系，需要通过决定系数R^{\mathrm{2}}和调整决定系数R_{\mathrm{a}}^{\mathrm{2}}检验。R^{\mathrm{2}}和R_{\mathrm{a}}^{\mathrm{2}}计算式如下。

R^{\mathrm{2}}=\mathrm{1}-\frac{S_{\mathrm{S}\mathrm{E}}}{S_{\mathrm{S}\mathrm{T}}}（2）

R_{\mathrm{a}}^{\mathrm{2}}=\mathrm{1}-\frac{n-\mathrm{1}}{n-k-\mathrm{1}}(\mathrm{1}-R^{\mathrm{2}})（3）

其中，

S_{\mathrm{S}\mathrm{T}}={\displaystyle \sum _{l=\mathrm{1}}^n}(y_l{-\overline{y})}^{\mathrm{2}}

S_{\mathrm{S}\mathrm{E}}={\displaystyle \sum _{l=\mathrm{1}}^n}(y_l{-y)}^{\mathrm{2}}

式中：S_{\mathrm{S}\mathrm{T}}为总平方和；S_{\mathrm{S}\mathrm{E}}为残差平方和；y_l，\overline{y}和y分别为第l点响应的试验值、平均值及响应面函数值；n为试验数目。

数据结果拟合后若决定系数R^{\mathrm{2}}和调整决定系数R_{\mathrm{a}}^{\mathrm{2}}的值大于0.9，则考虑该响应面函数足够逼近实际模型，越趋近于1则响应面函数就越能反映真实问题。经过检验后的响应面函数可用于预测捣固参数。考虑到Box-Behnken设计试验次数少、效率较高的优点，本研究中的响应面试验的组别设计采用Box-Behnken方法。

2.2　Box-Behnken捣固参数响应面试验设计

响应面函数的建立需遵循正交设计原则，以选择合理且有效的试验方案。在众多试验设计中，Box-Behnken试验设计是一种常见的基于3水平的2阶试验方案，能够揭示各因素与响应值之间的非线性关系。该方法相较于其他响应面设计方法的优势在于，能够在预测精度与运行经济性之间实现平衡。

根据大型机械捣固作业的实际需求，本研究选取捣固频率、捣固振幅和捣镐移速作为影响因子，以道床支承刚度和道砟间作用力超过临界值的比例作为响应变量，设计了3因素3水平的响应面试验。在Box-Behnken设计中，3水平的因素取值分别位于2个端点和中点位置，各因素水平的取值见表2。表2中，x_{\mathrm{1}}表示捣固频率，x_{\mathrm{2}}表示捣固振幅，x_{\mathrm{3}}表示捣镐移速。

在3因素3水平的Box-Behnken设计中，试验数量n按下式计算。

n=\mathrm{2}k(k-\mathrm{1})+C_{\mathrm{0}}（4）

式中：C_{\mathrm{0}}为中心点的数量。

本文采用3因素试验，选择3个中心点作为中心试验，包含12个角点及3个中心点，共计15组仿真试验，试验设计如图6所示。其中，1，0和－1为3水平取值，分别代表所选捣固参数的最大值、中间值和最小值。

2.3　改进粒子群算法

粒子群算法通过个体粒子在群体中的协同进化过程来生成与优化解。其算法假设范围内的每个粒子都具有独立的速度与位置，计算开始后会被随机赋予初速度与初始位置，并在限定的搜索空间内自由飞行，搜索当前最优粒子位置。在迭代计算过程中，由每个粒子个体极值与全局极值对比后，完成粒子的位置和速度的迭代更新。PSO算法的进化方程如下所示。

\begin{array}{l}{v}_{m,t}+\mathrm{1}=w{v}_{m,t}+c_{\mathrm{1}}{r}_{\mathrm{1}}(p_{m,t}-x_{m,t})+\\ c_{\mathrm{2}}{r}_{\mathrm{2}}(g_t-x_{m,t})\end{array}（5）

x_{m,t}+\mathrm{1}=x_{m,t}+v_{m,t}+\mathrm{1}（6）

式中：m为粒子编号；t为迭代次数，取值范围为\mathrm{0}\le t\le t_{\mathrm{m}\mathrm{a}\mathrm{x}}-\mathrm{1}，其中t_{\mathrm{m}\mathrm{a}\mathrm{x}}为最大迭代次数；v_{m,t}为第t次迭代时粒子m的速度；w为计算时的惯性权重；c_{\mathrm{1}}和c_{\mathrm{2}}分别为个体认知因子和社会认知因子；r_{\mathrm{1}}\mathrm{和}{r}_{\mathrm{2}}均为0~1之间的随机数；p_{m,t}为第t次迭代时粒子m的最佳位置；g_t为种群内截至第t次迭代时的全局最优分量；x_{m,t}为第t次迭代时粒子m的位置。

PSO算法在处理高维复杂优化问题时易陷入局部最优解。尽管粒子通过全局信息和局部信息调整自身位置，但在某些情况下，粒子群可能过早收敛至局部最优解而无法跳出。在PSO算法的寻优过程中，算法性能高度依赖于权重w、学习因子c_{\mathrm{1}}\mathrm{和}{c}_{\mathrm{2}}的取值，其大小直接影响全局搜索能力与局部搜索能力的平衡。不合理的权重值和学习因子都可能导致粒子群结果集中于边界区域，收敛速度过快，从而陷入局部最优解，难以找到全局最优解。

为解决这一问题，本文采用动态调整惯性权重和学习因子的方法，确保在迭代过程中粒子能够覆盖更广泛的搜索空间，从而有效平衡全局搜索与局部搜索能力。动态权重系数和学习因子调整后的表达式分别为

w=w_{\mathrm{m}\mathrm{a}\mathrm{x}}-\frac{t}{t_{\mathrm{m}\mathrm{a}\mathrm{x}}}(w_{\mathrm{m}\mathrm{a}\mathrm{x}}-w_{\mathrm{m}\mathrm{i}\mathrm{n}})（7）

c_{\mathrm{1}}=(c_{\mathrm{1}\mathrm{m}\mathrm{a}\mathrm{x}}-c_{\mathrm{1}\mathrm{m}\mathrm{i}\mathrm{n}})\frac{t_{\mathrm{m}\mathrm{a}\mathrm{x}}-t}{t_{\mathrm{m}\mathrm{a}\mathrm{x}}}+c_{\mathrm{1}\mathrm{m}\mathrm{i}\mathrm{n}}（8）

c_{\mathrm{2}}=(c_{\mathrm{2}\mathrm{m}\mathrm{a}\mathrm{x}}-c_{\mathrm{2}\mathrm{m}\mathrm{i}\mathrm{n}})\frac{t}{t_{\mathrm{m}\mathrm{a}\mathrm{x}}}+c_{\mathrm{2}\mathrm{m}\mathrm{i}\mathrm{n}}（9）

式中：w_{\mathrm{m}\mathrm{a}\mathrm{x}}和w_{\mathrm{m}\mathrm{i}\mathrm{n}}分别为最大及最小的权重因子；c_{\mathrm{1}\mathrm{m}\mathrm{a}\mathrm{x}}，c_{\mathrm{1}\mathrm{m}\mathrm{i}\mathrm{n}}，c_{\mathrm{2}\mathrm{m}\mathrm{a}\mathrm{x}}\mathrm{和}{c}_{\mathrm{2}\mathrm{m}\mathrm{i}\mathrm{n}}分别为算法运算时个体认知因子、社会认知因子的最大值及最小值。

为验证PSO算法的正确性与准确性，采用经典的多峰测试函数Rastrigin对IPSO算法、PSO算法、GA算法进行验证。算法种群数量为30，最大迭代次数设为500。其中，IPSO算法的w_{\mathrm{m}\mathrm{a}\mathrm{x}}= 0.9，w_{\mathrm{m}\mathrm{i}\mathrm{n}}=\mathrm{0.4}，c_{\mathrm{1}\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{2}，c_{\mathrm{2}\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{2}，c_{\mathrm{1}\mathrm{m}\mathrm{i}\mathrm{n}}=\mathrm{0.5}，c_{\mathrm{2}\mathrm{m}\mathrm{i}\mathrm{n}}=\mathrm{0.5}。3种优化算法测试结果见表3。

由表3可知，IPSO算法的平均迭代次数较少，其最优值标准差也相对较小。图7为3种优化算法运行时的收敛曲线。综合表3和图7中的相关数据也可看出，IPSO算法具有更好的全局搜索性能和稳定性。相较于标准PSO算法，改进后的权重及学习因子动态变化的方式提高了PSO算法的收敛速度和搜索可靠性，改善了PSO算法易陷入局部最优解的问题。

3 捣固质量影响因素

3.1　道床支承刚度

道床支承刚度直接影响道床对轨道结构的支承能力。捣固后的道床较为松散，密实度下降，因此研究不同捣固参数组合下的道床支承刚度变化，有助于从宏观角度分析捣固后道床安全性的差异。采用Box-Behnken设计方法设计的试验组别进行数值仿真，并提取了15组试验后的道床支承刚度，结果如图8所示。

由图8可知，不同条件下支承能力存在较大差异，角点组的支承刚度普遍低于中心组，且差异较为明显。中心组支承刚度平均为64.71 kN · mm^－1，高于其余角点组。将图8所示的道床支承刚度各组试验结果采用响应面函数模型拟合数据，并利用最小二乘法进行方差和统计分析，得到道床支承刚度响应面函数Q如下。

Q=-\mathrm{222.581}\mathrm{8}+\mathrm{1.993}\mathrm{2}{x}_{\mathrm{1}}+\mathrm{3.539}\mathrm{7}{x}_{\mathrm{2}}+\mathrm{408.656}\mathrm{3}{x}_{\mathrm{3}}-\mathrm{0.015}\mathrm{1}{x}_{\mathrm{1}}^{\mathrm{2}}-\mathrm{0.191}\mathrm{3}{x}_{\mathrm{2}}^{\mathrm{2}}-\mathrm{167.942}\mathrm{7}{x}_{\mathrm{3}}^{\mathrm{2}}-\mathrm{0.005}\mathrm{3}{x}_{\mathrm{1}}{x}_{\mathrm{2}}-\mathrm{0.637}\mathrm{5}{x}_{\mathrm{1}}{x}_{\mathrm{3}}-\mathrm{0.038}\mathrm{5}{x}_{\mathrm{2}}{x}_{\mathrm{3}}（10）

获取道床支承刚度响应面函数后，需要对该函数进行方差分析，判断数据的拟合精度效果。方差分析结果见表4。其中，平方和能够衡量数据变异程度；均方是平方和除以相应的自由度，代表平均变异；F为检验响应面函数显著性的参数，其值越大，则对该目标的影响越大；P为响应面函数拒绝原假设的值，若小于0.05则表明具有统计学意义，函数拟合结果较好^［27］。

由表4可知，Q的显著性检验参数P=0.000 4，小于0.05，模型拟合精度较好。另一方面，由式（2）和式（3）计算可知，道床支承刚度的响应面函数的决定系数R^{\mathrm{2}}为0.99，调整决定系数R_{\mathrm{a}}^{\mathrm{2}}为0.96，响应面函数具有较好的拟合精度。

从显著性参数F来看，在单因素影响下，捣固频率、捣固振幅及捣镐移速的F分别为5.20，0.58和40.51。因此，影响道床支承刚度的主次顺序为：捣镐移速>捣固频率>捣固振幅。此外，捣固频率与捣镐移速的交互影响即x_{\mathrm{1}}{x}_{\mathrm{3}}的F为11.31，大于其他交互影响的F取值，在道床支承刚度的影响中最为显著。

基于道床支承刚度响应面函数，得到道床支承刚度与捣固频率、捣固振幅及捣镐移速的等值线图，如图9所示。

由图9可知：道床支承刚度等值线图中捣镐移速方向的等值线分布更为密集，相比之下，捣固频率与振幅方向的变化较为平缓，表明捣镐移速对道床支承刚度的影响最为显著，这与方差分析结果一致；无论是捣固频率、捣固振幅还是捣镐移速，对道床支承刚度的影响均呈现先增大后减小的趋势。因此，各影响因素的取值应选择适中值，更有利于捣固后道床的稳定性。

3.2　道砟间作用力临界值超限分数

根据文献［12］可知，当道砟间作用力超过10 kN时，道砟破碎的概率显著增加。基于此，以10 kN作为临界值，分析不同捣固参数条件下道砟间接触力超过10 kN的比例。文献［11］指出，下插阶段的道砟间作用力通常大于夹持阶段，且出现较大作用力的概率较高。因此，本文以下插阶段的道砟间作用力作为分析对象。

由于单个道砟间的作用力无法全面反映道床的力学特性，本文选取道砟间作用力超过临界值f_{\mathrm{c}}的数目与道砟间作用力总数的比值作为分析指标，定义该比值为道砟间作用力临界值超限分数\eta（后续简称为临界值超限分数），其计算式如下。

\eta =\frac{D_{f>f_{\mathrm{c}}}}{D_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}\times \mathrm{100}\mathrm{\%}（11）

式中：D_{f>f_{\mathrm{c}}}为道砟间作用力f超过临界值f_{\mathrm{c}}的数目；D_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}为道砟间作用力的总数目。

根据设计试验得到的临界值超限分数如图10所示。

由图10可知，临界值超限分数在不同设计组别中的变化规律与道床支承刚度并不一致。中心组试验结果并未表现为最优，其平均值仅为0.76，高于部分试验组。

将图10所示的临界值超限分数各组试验结果采用响应面函数模型拟合数据，并采用最小二乘法进行方差和统计分析，得到临界值超限分数响应面函数如下。

\begin{array}{l}\eta =\mathrm{3.076}\mathrm{22}-\mathrm{0.031}\mathrm{45}{x}_{\mathrm{1}}-\mathrm{0.132}\mathrm{71}{x}_{\mathrm{2}}-\\ \mathrm{4.292}\mathrm{19}{x}_{\mathrm{3}}+\mathrm{0.000}\mathrm{58}{x}_{\mathrm{1}}^{\mathrm{2}}+\mathrm{0.005}\mathrm{85}{x}_{\mathrm{2}}^{\mathrm{2}}+\\ \mathrm{2.518}\mathrm{23}{x}_{\mathrm{3}}^{\mathrm{2}}-\mathrm{0.000}\mathrm{42}{x}_{\mathrm{1}}{x}_{\mathrm{2}}-\\ \mathrm{0.006}\mathrm{56}{x}_{\mathrm{1}}{x}_{\mathrm{3}}+\mathrm{0.066}\mathrm{67}{x}_{\mathrm{2}}{x}_{\mathrm{3}}\end{array}（12）

得到上述响应面函数后，需对上述的临界值超限分数响应面函数进行方差分析，判断该函数的拟合精度效果，结果见表5。

由表5可知，\eta的显著性检验参数P<0.000 1，说明该模型具有较好的精度。另一方面，由式（2）和式（3）计算可知，\eta的决定系数R^{\mathrm{2}}为0.99，调整决定系数R_{\mathrm{a}}^{\mathrm{2}}为0.99，表明该模型具有较高的拟合度，可用于后续分析。

从显著性参数F来看，在单因素的影响下，捣固频率、捣固振幅以及捣镐移速的F分别为3.43，83.00和1 155.30，表明对临界值超限分数的影响程度依次为：捣镐移速>捣固振幅>捣固频率。此外，从交互性影响的显著性参数可以看出，捣固振幅与捣镐移速的交互影响即x_{\mathrm{2}}{x}_{\mathrm{3}}的F为20.53，大于其他交互影响的F取值，临界值超限分数的影响最为显著。

基于临界值超限分数响应面函数，得到临界值超限分数与捣固频率、捣固振幅与捣镐移速的等值线图，如图11所示。

由图11可知：捣镐移速方向的数值变化梯度更大，变化区域更为集中，相较之下，捣固频率与振幅方向的变化幅度较小，说明捣镐移速对临界值超限分数的影响作用更为突出，这与方差分析结果一致，且其影响程度较道床支承刚度作为响应时更为明显。

与道床支承刚度响应面函数不同的是，捣镐移速对临界值超限分数的影响几乎呈持续增加趋势，即捣镐移速越大，道砟受破坏的概率越高；相比之下，捣固频率对临界值超限分数的影响呈现先增大后减小的趋势，选择适中的捣固频率更有利于降低临界值超限分数。

4 基于响应面法与改进粒子群算法的捣固参数优化设计

本文基于响应面法与改进粒子群算法的优化分析，旨在确保道床支承刚度尽可能大的同时，降低道砟间作用力达到破坏临界值的概率，从而满足有砟轨道系统维护的可持续性要求，并延长有砟道床的使用寿命。基于数值仿真结果及响应面函数分析，应用前述优化策略对捣固作业参数进行设计，以获得最优捣固方案。

由于道床支承刚度与临界值超限分数的数值差异较大，为便于求解，对道床支承刚度和临界值超限分数的响应面函数取值进行归一化处理，计算式如下。

Q_{\mathrm{1}}=\frac{Q-Q_{\mathrm{m}\mathrm{i}\mathrm{n}}}{Q_{\mathrm{m}\mathrm{a}\mathrm{x}}-Q_{\mathrm{m}\mathrm{i}\mathrm{n}}}（13）

{\eta }_{\mathrm{1}}=\frac{\eta -{\eta }_{\mathrm{m}\mathrm{i}\mathrm{n}}}{{\eta }_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\eta }_{\mathrm{m}\mathrm{i}\mathrm{n}}}（14）

式中：Q_{\mathrm{1}}为归一化后道床支承刚度响应面函数取值；{\eta }_{\mathrm{1}}为归一化后临界值超限分数响应面函数取值；Q_{\mathrm{m}\mathrm{a}\mathrm{x}}和Q_{\mathrm{m}\mathrm{i}\mathrm{n}}分别为取值范围内道床支承刚度的最大值和最小值；{\eta }_{\mathrm{m}\mathrm{a}\mathrm{x}}和{\eta }_{\mathrm{m}\mathrm{i}\mathrm{n}}分别为取值范围内临界超限分数的最大值和最小值。

利用响应面函数Q及\eta偏导数等于0的方法求解Q_{\mathrm{m}\mathrm{a}\mathrm{x}}，Q_{\mathrm{m}\mathrm{i}\mathrm{n}}，{\eta }_{\mathrm{m}\mathrm{a}\mathrm{x}}和{\eta }_{\mathrm{m}\mathrm{i}\mathrm{n}}，可得Q_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{65.62}，Q_{\mathrm{m}\mathrm{i}\mathrm{n}}=\mathrm{13.94}，{\eta }_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{2.99}，{\eta }_{\mathrm{m}\mathrm{i}\mathrm{n}}=\mathrm{0.30}。为了能够更好地平衡道床支承刚度Q最大化以及临界值超限分数\eta的最小化，增加权重系数\lambda，利用归一化函数后的Q_{\mathrm{1}}值和{\eta }_{\mathrm{1}}值，构造加权函数\mu，即

\mu =-\lambda Q_{\mathrm{1}}+(\mathrm{1}-\lambda ){\eta }_{\mathrm{1}}（15）

利用IPSO算法对构造后的加权函数寻优分析，粒子数目设定为30，迭代次数设为500，自变量取值范围与响应面函数时一致。由于单次寻优结果具有一定的随机性，本节考虑权重系数\lambda从0到1.0共11种组合，分别计算不同系数下的道床支承刚度和临界值超限分数，并求解不同权重系数下的最优解，结果如图12所示。其中，当\lambda =\mathrm{1}时表示仅考虑道床支承刚度最大化，\lambda =\mathrm{0}时则表示仅考虑临界值超限分数最小化。

由图12可知：随着权重系数\lambda的增加，道床支承刚度逐渐增大，并在权重系数达到0.4后趋于平稳，临界值超限分数则随权重系数的增加而持续上升；仅考虑道床支承刚度最大（\lambda =\mathrm{1}）与\lambda =0.4时，Q与\eta分别由65.62 kN · mm^－1及0.63降低至63.69 kN · mm^－1及0.46；在道床支承刚度基本保持不变的前提下，临界值超限分数明显降低。

进一步考虑现场可操作性，将权重系数为0.4时的最优捣固参数（捣固频率39.09 Hz、捣固振幅7.77 mm、捣镐移速1.04 m · s^－1）按整数取为捣固频率39 Hz、捣固振幅8 mm、捣镐移速1 m · s^－1，此时的Q和\eta分别为62.15 kN · mm^－1和0.42。相较于仅考虑道床支承刚度最优值的解，Q仅降低5.29%，而\eta可降低33.33%。因此，基于目标函数的优化分析可实现道床支承刚度变化较小，而临界值超限分数明显下降的实际需求，从而降低道砟破坏概率。

5 结论

（1）通过数值仿真试验，利用响应面法确定了道床支承刚度Q和临界值超限分数\eta的响应函数，并建立了与捣固频率、捣固振幅及捣镐移速的回归关系，响应面函数的R_a^{\mathrm{2}}分别为0.96和0.99，模型预测精度较高。

（2）影响道床支承刚度的捣固参数按显著性从大到小排序依次为捣镐移速、捣固频率和捣固振幅；影响临界值超限分数的捣固参数按显著性从大到小排序依次为捣镐移速、捣固振幅及捣固频率。

（3）基于改进粒子群算法与响应面函数结合的方法，实现了捣固参数的优化求解。综合分析优化结果与大机作业可行性，最优捣固工艺方案为：捣固频率39 Hz，捣固振幅8 mm，捣镐移速1 m · s^－1。与仅考虑道床支承刚度最大化的方法相比，采用优化方法选择的捣固参数可以在道床支承刚度数值上基本不变的情况下，将临界值超限分数占比降低33.33%。该优化方案提升了有砟道床在捣固后的稳定性与耐久性。

参考文献

[1] 高亮，徐旸，杨国涛，等.铁路有砟道床劣化研究进展综述［J］.铁道学报，2022，44（8）：78-92．

[2] GAO Liang， XU Yang， YANG Guotao， et al. Review of Research of Railway Ballast Bed Deterioration ［J］. Journal of the China Railway Society， 2022， 44 （8）： 78-92. in Chinese

[3] GUO Y L， MARKINE V， JING G Q. Review of Ballast Track Tamping： Mechanism， Challenges and Solutions ［J］. Construction and Building Materials， 2021， 300： 123940.

[4] ZHOU T Y， HU B， SUN J F， et al. Discrete Element Method Simulation of Railway Ballast Compactness during Tamping Process ［J］. The Open Electrical & Electronic Engineering Journal， 2013， 7 （1）： 103-109.

[5] AINGARAN S， LE PEN L， ZERVOS A， et al. Modelling the Effects of Trafficking and Tamping on Scaled Railway Ballast in Triaxial Tests ［J］. Transportation Geotechnics， 2018， 15： 84-90.

[6] LIU J X， WANG P， LIU G Z， et al. Influence of a Tamping Operation on the Vibrational Characteristics and Resistance-Evolution Law of a Ballast Bed ［J］. Construction and Building Materials， 2020， 239： 117879.

[7] ZHANG Z H， XIAO H， ZHU Y J， et al. Macro-Meso Mechanical Properties of Ballast Bed during Three-Sleeper Tamping Operation ［J］. International Journal of Rail Transportation， 2023， 11 （6）： 886-911.

[8] 张智海，肖宏，刘光鹏，等.三枕捣固作业中捣镐夹持力对有砟道床作用效应影响研究［J］.铁道学报，2024，46（10）：132-143.

[9] ZHANG Zhihai， XIAO Hong， LIU Guangpeng， et al. Effect of Squeezing Force of Tamping Picks on Ballast Bed during Three-Sleeper Tamping Operation ［J］. Journal of the China Railway Society， 2024， 46 （10）： 132-143. in Chinese

[10] SAUSSINE G， AZÉMA E， PERALES R， et al. Compaction of Railway Ballast during Tamping Process： a Parametric Study ［C］∥ Powders and Grains 2009： Proceedings of the 6th International Conference on Micromechanics of Granular Media. Golden （Colorado）， USA. Colorado： American Institute of Physics， 2009， 1145 （1）： 469-472.

[11] DESCANTES Y， PERALES R， SAUSSINE G， et al. On the Damaging Effects of the Ballast Tamping Operation ［C］∥ 9th World Congress on Railway Research. Lille， France. UK： Rail Safety and Standards Board Ltd.， 2011： 1-10.

[12] WANG C Y， SHI S W， GAO L， et al. Influence of Consecutive Tamping Operations on Failures of Ballast Bed and Tamping Machine ［J］. Engineering Failure Analysis， 2024， 163： 108583.

[13] 高亮，石顺伟，蔡小培，等.捣固装置与有砟道床相互作用及捣固频率的影响研究［J］.工程力学，2024，41（12）：30-41.

[14] GAO Liang， SHI Shunwei， CAI Xiaopei， et al. Research on Interaction of Tamping Machine-Ballast Bed and Influence of Tamping Frequency ［J］. Engineering mechanics， 2024， 41 （12）： 30-41. in Chinese

[15] SHI S W， GAO L， CAI X P， et al. Mechanical Characteristics of Ballasted Track under Different Tamping Depths in Railway Maintenance ［J］. Transportation Geotechnics， 2022， 35： 100799.

[16] SHI S W， GAO L， XIAO H， et al. Research on Ballast Breakage under Tamping Operation Based on DEM-MBD Coupling Approach ［J］. Construction and Building Materials， 2021， 272： 121810.

[17] 占玉林，许江辉，许俊，等.基于响应面法和粒子群算法的桥梁高耸临时提升支架优化［J］.中国铁道科学，2022，43 （6）：39-46.

[18] ZHAN Yulin， XU Jianghui， XU Jun， et al. Optimization of Bridge High Temporary Lifting Support Based on Response Surface Method and Particle Swarm Algorithm ［J］. China Railway Science， 2022， 43 （6）： 39-46. in Chinese

[19] 杨宏印，吴楠昊，王波，等.基于贝叶斯法和响应面法的重载铁路桥梁模型修正［J］.中国铁道科学，2022，43（6）：47-54.

[20] YANG Hongyin， WU Nanhao， WANG Bo， et al. Model Updating of Heavy-Haul Railway Bridge Based on Bayesian Method and Response Surface Method ［J］. China Railway Science， 2022， 43 （6）： 47-54. in Chinese

[21] CHI Y H， XIAO H， ZHANG Z H， et al. Analysis of the Influence of Vibration Frequency and Amplitude on Ballast Bed Tamping Operation in Railway Turnout Areas ［J］. Computational Particle Mechanics， 2024， 11 （2）： 771-788.

[22] ZHOU T Y， HU B， PENG Y X， et al. Research on Vibration Parameters of Tamping Operation ［J］. Applied Mechanics and Materials， 2015， 724： 283-286.

[23] SHI S W， WANG J， GAO L， et al. Intelligent Optimization of Tamping Parameters Using Discrete Element and Radial Basis Function-Multi-Objective Genetic Algorithm （RBF-MOGA） ［J］. Automation in Construction， 2024， 168： 105752.

[24] 边学成，李伟，李公羽，等.基于颗粒真实几何形状的铁路道砟剪切过程三维离散元分析［J］.工程力学，2015，32（5）：64-75.

[25] BIAN Xuecheng， LI Wei， LI Gongyu， et al. Three-Dimensional Discrete Element Analysis of Railway Ballast’s Shear Process Based on Particles’ Real Geometry ［J］. Engineering Mechanics， 2015， 32 （5）： 64-75. in Chinese

[26] 国家铁路局. TB 10082—2017 铁路轨道设计规范 ［S］.北京：中国铁道出版社，2017.

[27] National Railway Administration of the People’s Republic of China. TB 10082—2017 Code for Design of Railway Track ［S］. Beijing： China Railway Publishing House， 2017. in Chinese ）

[28] MINDLIN R D. Compliance of Elastic Bodies in Contact ［J］. Journal of Applied Mechanics， 1949， 16 （3）： 259-268．

[29] MINDLIN R D， DERESIEWICZ H. Elastic Spheres in Contact under Varying Oblique Forces ［J］. Journal of Applied Mechanics， 1953， 20 （3）： 327-344．

[30] XIAO J L， LIU G Z， LIU J X， et al. Parameters of a Discrete Element Ballasted Bed Model Based on a Response Surface Method ［J］. Journal of Zhejiang University： Science A （Applied Physics ＆ Engineering）， 2019， 20 （9）： 685-700．

[31] SHI S W， GAO L， CAI X P， et al. Effect of Tamping Operation on Mechanical Qualities of Ballast Bed Based on DEM-MBD Coupling Method ［J］. Computers and Geotechnics， 2020， 124： 103574.

[32] 刘浩，杨国涛，江成，等.高速铁路有砟道床状态评定参数关联关系分析［J］.铁道学报，2021，43（6）：128-134.

[33] LIU Hao， YANG Guotao， JIANG Cheng， et al. Analysis on Correlation between State Evaluation Parameters of Ballasted Bed in High-Speed Railway ［J］. Journal of the China Railway Society， 2021， 43 （6）： 128-134. in Chinese

[34] 井国庆，付豪，贾文利，等.高速铁路有砟轨道横向阻力宏细观分析［J］，铁道工程学报，2018，35（9）：21-25

[35] JING Guoqing， FU Hao， JIA Wenli， et al. Macro-Micro Analysis of Lateral Resistance for High Speed Railway Ballasted Track ［J］. Journal of Railway Engineering Society， 2018， 35 （9）： 21-25. in Chinese

[36] 赵新元，杨科，何祥，等.基于RSM-BBD的源煤基固废胶结体配比及性能研究［J］.材料导报，2024，38（9）：153-159.

[37] ZHAO Xinyuan， YANG Ke， HE Xiang， et al. Study on Proportioning and Performance of Cemented Body from Multi-Source Coal-Based Solid Waste Based on RSM-BBD Experiment ［J］. Materials Reports， 2024， 38 （9）： 153-159. in Chinese

基金资助

国家重点研发计划项目(2022YFB2602900)

国能铁路装备有限责任公司技术开发项目(TZKY22-12)