一种直线导轨与滑块间刚度和阻尼的辨识方法

《Precision Engineering》：An identification method for stiffness and damping between linear guideway and carriage

　　

• 提出一种基于加速度信号的质心位移重构方法。
• 建立滑块-导轨接触的连续分布弹性阻尼模型。
• 使用拉格朗日方程推导五自由度动力学模型。
• 同步识别刚度与阻尼参数。
• 通过多次实验验证识别参数的精度。

1. 引言
线性导轨是数控机床进给系统中的关键部件，其力学性能直接影响机床的结构动态特性和最终加工精度 [1]。循环滚珠通常布置在导轨与滑块之间，通过滚动接触实现高精度、低摩擦的直线运动 [[2], [3], [4], [5]]。与传统滑动导轨相比，该结构不仅显著减少了磨损和能量损失，而且以更低的维护成本确保了稳定的运动质量 [6]。然而，当滑块内部的滚珠受到复杂的多方向外部载荷时，滚珠与滚道之间的局部接触区域会产生高表面压强，导致弹性压缩变形 [7]。这种变形激发了滑块在非进给方向上的耦合五自由度固有振动 [8,9]。这些振动导致定位偏差和姿态误差，从而影响高精度加工场景中的轨迹跟踪能力 [10]。因此，理解滚珠-滚道界面的动态行为、建立其动力学模型以及识别刚度和阻尼参数，对于预测运动响应和抑制振动至关重要。

关于线性导轨的振动特性，现有研究在多个方面取得了进展。在动态建模和仿真方法方面，研究逐渐从静态到动态、从刚体到柔性体、从单自由度到多自由度的方法演进。早期研究主要关注静态或准静态建模：Sun 等人 [11] 结合数值方法和有限元方法建立了线性导轨的单自由度静态模型；Liu 等人 [12,13] 引入几何误差开发了经实验验证的静态刚度模型；Kong 等人 [14] 进行了早期动态探索，建立了线性导轨的非线性动态模型，重点关注垂直振动和稳定性。随着研究的深入，多自由度和柔性体建模成为关键领域：在传统刚性模型的基础上，学者们相继开发了线性导轨的多自由度静态和动态模型 [15,16]。Yang 等人 [17] 分析了方向刚度及其在复杂载荷下的变化；Wang 等人 [18] 在 Kong 的工作基础上，提出了工作台系统的三自由度动态模型和振动抑制方法；Liu 等人 [19] 进一步引入装配误差建立了五自由度动态模型。同时，柔性体建模逐渐受到关注：Ohta 和 Tanaka [20] 将滑块视为柔性体，研究预加载导轨的垂直刚度并通过实验验证；Ma 等人 [21] 考虑滑块柔性的建立了三自由度静态刚度模型；Kwon 和 Tong 等人 [22,23] 提出了五自由度刚度模型，系统比较了柔性与刚性模型之间的差异；Cheng 等人 [24] 基于能量平衡原理，开发了柔性滚轮导轨的刚度矩阵模型，并讨论了外部载荷和预载荷对刚度的影响。在仿真方法与应用方面，数值分析和有限元方法得到了广泛应用：Zou 和 Wang [25] 结合这些方法研究了预载荷、热变形和磨损综合作用下的刚度演化；Pawe?ko 和 Dunaj 等人 [26,27] 使用有限元方法开发了简化刚度模型并通过实验验证了有效性；Majda [28] 通过有限元建模系统地研究了导轨几何误差对工作台运动精度的影响。此外，特定工作条件下的动态特性也引起了关注：Soleimanian 等人 [29] 建立了考虑润滑效应的摩擦动力学模型；Yu 等人 [30] 探索了磨损对导轨全生命周期动态可靠性的影响。这些研究共同推动了线性导轨动态建模和仿真方法的发展和完善。这些研究不断完善了线性导轨动态建模的理论框架，并增强了对系统动态特性的描述能力。然而，大多数现有多自由度动态模型都是根据牛顿运动定律逐自由度建立的，其中平动和转动自由度之间的几何关系通常通过额外的坐标变换来引入，导致动力学的独立建模和运动学描述。另一方面，虽然柔性模型可以更准确地描述滚珠-滚道界面的局部接触变形，但它们通常需要大量结构参数，对边界条件敏感，并且计算复杂度高 [31,32]。因此，对于面向工程应用的系统动态分析和参数识别，刚体多自由度模型在计算效率、模型简洁性和工程实用性方面仍具有显著优势 [[33], [34], [35]]。

随着正向动力学建模研究的推进，准确的参数识别在模型构建和验证中变得越来越重要。Zeng [36] 早期采用了数据处理和参数识别方法来提取与轴承刚度和阻尼密切相关的模态频率和阻尼比。Mao 等人 [37] 基于频率响应函数矩阵与动态刚度矩阵之间的逆关系，提出了主轴系统动态参数识别的解析方法。在磨削领域，Leonesio 等人 [38] 通过实验频率响应函数拟合识别了特定的磨削能和法向力比值，随后推断出磨削过程的动态刚度和阻尼特性。Qiao 等人 [39] 整合了最小二乘估计和蒙特卡洛模拟，使用运动误差数据实现了对主轴轴承动态刚度和阻尼的系统识别。然而，这些方法主要面向轴承和主轴等旋转机械系统，识别的参数大多是整体等效刚度和阻尼系数。关于线性导轨滚珠-滚道接触界面处连续分布刚度和阻尼参数识别的研究仍然相对有限。此外，现有研究主要使用频率响应函数或位移响应进行参数识别，而结合冲击试验加速度响应和信号处理的线性导轨系统动态参数系统化识别程序的开发仍需要进一步完善。

基于上述研究现状，本文提出了一种五自由度线性导轨系统的动态参数识别方法。首先，使用拉格朗日方程建立了滑块的五自由度动力学模型。在统一的广义坐标框架下，动力学方程与坐标变换有机融合，从而对平动和转动之间的耦合关系进行了一致描述。滑块与导轨之间的滚动接触被建模为等效零长度弹簧系统，并进一步近似为四条具有单位刚度和单位阻尼的连续分布支撑线，以表征滚珠-滚道接触界面的动态行为。随后，进行冲击激发实验以获得滑块的多方向加速度响应。通过结合频域积分和多节点三次样条插值重构位移信息。然后，使用 Newmark-β 数值积分方法和最小二乘估计同步识别四条支撑线的单位长度刚度和阻尼参数。最后，将识别的参数应用于正向动力学仿真，并通过实验验证模型的预测精度。该方法不仅获得了具有明确物理意义的接触参数，还揭示了线性导轨接触刚度和阻尼的分布特性以及滑块在激发下的姿态变化。它为线性导轨系统动态响应的高精度预测、振动抑制和结构优化提供了可靠的理论基础。

本文其余部分安排如下：第2节介绍理论模型的建立，包括线性导轨动力学方程的推导以及刚度和阻尼参数的识别方法。第3节介绍滑块质心位移重构方法。第4节介绍仿真验证。第5节报告实验研究，第6节进行讨论，第7节总结本文的主要结论。

2. 动态参数识别方法
2.1 线性导轨的动态建模
如图1所示，本研究关注一个线性导轨系统，在滑块顶部添加了一个金属块作为负载单元。这种布置便于施加多方向载荷，从而激发系统的五自由度振动并模拟实际负载条件。基于使用软件进行的三维建模和质心分析，确定了滑块本体和所加金属块的质心位置，具体分布如图2所示。

图1. 线性导轨的三维模型。
图2. 滑块质心位置示意图。

图3显示了线性导轨的典型横截面结构。滑块包含四排钢球，每排沿着周期性循环路径运行以确保在任何安装方向上都能稳定运行。图3进一步说明了具有升高结构的滑块中上下排滚珠循环路径的示意图 [40]，其中 L 代表直线沟槽工作区的长度。

图3. 线性导轨结构示意图。

为建立该系统的理论动力学模型，做出以下基本假设：
(1) 滑块和导轨均视为刚体，不发生宏观结构变形；
(2) 系统变形仅发生在滚动接触界面和滚动元件的局部接触区域；
(3) 滚动接触界面表面为理想几何面，初始接触状态完全由滚道和滚动元件半径决定；
(4) 导轨以理想直线状态安装，滑块在运行过程中任何姿态变化都将直接影响接触条件；
(5) 导轨与滑块之间的接触可简化为零长度弹簧模型；
(6) 四排滚珠与导轨之间的接触可等效表示为沿四条支撑线连续分布的弹性支撑单元，每条支撑线具有单位长度刚度和阻尼。这四条支撑线位于滚珠与滑块的接触点。每条支撑线的分布范围由第一个和最后一个受载滚珠的位置决定。对于第 j 条支撑线（j = 1, 2, 3, 4），起始和结束坐标分别记为 x_aj 和 x_bj，如图3所示。

单位长度的刚度和阻尼参数定义如下：
(1) { dk = ζ ds; dc = ψ ds }
其中 k 和 c 分别代表支撑线的总刚度和总阻尼；s 表示支撑线方向上的单位长度；ζ 和 ψ 分别对应单位长度的等效刚度系数和等效阻尼系数。

基于滑块严格结构对称和均一的假设，本研究采用分布参数模型对其进行描述。假设对于第 j 条支撑线，沿纵向（Y方向）单位长度横向刚度为 ζ_yj，垂直方向（Z方向）单位长度刚度为 ζ_zj。相应地，单位长度横向阻尼系数为 ψ_yj，垂直方向阻尼系数为 ψ_zj，条件为 ζ_y1 = ζ_y2 = ζ_y3 = ζ_y4，ζ_z1 = ζ_z2 = ζ_z3 = ζ_z4，ψ_y1 = ψ_y2 = ψ_y3 = ψ_y4，以及 ψ_z1 = ψ_z2 = ψ_z3 = ψ_z4。该模型具有连续分布的刚度和阻尼，更符合实际系统中滚珠的连续支撑特性。图4说明了线性导轨的动力学模型。

图4.

直线导轨动力学模型。图 5 说明了直线导轨的坐标系：惯性坐标系 OXYZ，原点 O 位于滑座系统质心，X 轴与水平横向方向一致，Y 轴与水平纵向方向一致，Z 轴与导轨垂直方向一致；平移坐标系 Gx'y'z'，G 为滑座系统质心，无振动时与点 O 重合，并随工作台移动，其轴与 OXYZ 平行；以及旋转坐标系 Gx"y"z"，它也随工作台移动。
下载：下载高分辨率图像 (333KB)
下载：下载全尺寸图像
图 5. 直线导轨的坐标系。
由于滑座不沿导轨随激振移动，激振下滑座质心的坐标为：
(2) ??=[0,??,??]T
在小角度 θx, θy, θz 的假设下，从 Gx'y'z' 到 Gx"y"z" 的旋转矩阵为
(3) ??=??????????????=?? ? ??1??????????????1?????????????1?? ? ??
(4) ??T=???1=?? ? ??1?????????????1??????????????1?? ? ??
其中：
(5) ????=?? ??10001????0?????1?? ???????=?? ? ??10?????010????01?? ? ???????=?? ??1????0?????10001?? ??
表 1 显示了支撑线两端弹簧对应的坐标值。这些坐标定义了滚珠支撑系统沿直线导轨长度的边界位置，为后续分布参数模型中刚度和阻尼的连续分布提供了几何参考。
表 1. 惯性坐标系中等效弹簧的坐标值。

| 弹簧排号 j | 第一列弹簧坐标 | 最后一列弹簧坐标 |
| :--- | :--- | :--- |
| | xaj, yaj, zaj | xbj, ybj, zbj |
| 1 | lax, ly, luzz | lbx, ly, luzz |
| 2 | lax, -ly, luzz | lbx, -ly, luzz |
| 3 | lax, -ly, ldz | lbx, -ly, ldz |
| 4 | lax, ly, ldz | lbx, ly, ldz |

当系统处于静平衡位置（即无振动）时，弹簧在坐标系 OXYZ 和 Gx"y"z" 中重合：
(6) ???????=??″?????=[???????,???????,???????]T, ???????=??″?????=[???????,???????,???????]T, ??=1,2,3,4
其中 xaj 是第 j 排第一列弹簧的坐标；xbj 是第 j 排最后一列弹簧的坐标。
系统运行期间，支撑点在参考系 Gx"y"z" 中的坐标保持不变，而它们在 Gx'y'z' 和 OXYZ 坐标系中的坐标值，分别记为 ??′?????, ??′????? 和 ????????, ????????，如下所示：
(7) ??′?????=??T???″?????, ??′?????=??T???″?????
(8) ????????=??+??T???″?????, ????????=??+??T???″?????
其中：
(9) ????????=?? ? ???????????????????????????? ? ??=?? ??0?????? ??+?? ? ??1?????????????1??????????????1?? ? ????? ? ????????????????????????? ? ??=?? ? ??????????????????????+??????????????+????????????+???????????????????????????????????+????????????+????????? ? ??
(10) ????????=?? ? ???????????????????????????? ? ??=?? ??0?????? ??+?? ? ??1?????????????1??????????????1?? ? ????? ? ????????????????????????? ? ??=?? ? ??????????????????????+??????????????+????????????+???????????????????????????????????+????????????+????????? ? ??
支撑线的方向由一个单位向量描述，可表示为：
(11) ????=(???????????????????,???????????????????,???????????????????)
以支撑线的起点 ???????? 为参考，引入参数 s 来表示从起点沿线段到任意位置的有向距离 (0 < s < L)。然后，振动后线段上任意点的位置 Pj(s) 可以参数化描述为：
(12) ?{ { {?{ { {??????(??)=????????+??·????????????????????????(??)=????????+??·????????????????????????(??)=????????+??·???????????????????
振动前两条接触线段上任意点的位置 Pj0(s) 可以参数化为：
(13) ?{ {?{ {??????0?(??)=???????+??·??????????????????????0?(??)=????????????0?(??)=???????
利用拉格朗日方程建立运动方程：
(14) dd?????(?????)?˙????(?????)???+????˙??=??
其中 q 是振动系统的广义坐标；Q 是系统广义力。
(15) ??=[??,??,????,????,????]T
振动系统的动能 T 可表示为：
(16) ??=12????˙??2+12????˙??2+12??????˙??2??+12??????˙??2??+12??????˙??2??
系统的势能 V 表示为：
(17) ??=12?4∑??=1∫??0(?????(??)??????0?(??))T??????(?????(??)??????0?(??))????
其中 Kj 是第 j 排支撑线的单位长度刚度矩阵，????=diag?(0,???????,???????)。基于均匀性假设，每排支撑线的刚度分布可得到：ζyj = ζy, ζzj = ζz。
因此，系统的最终势能表示为：
(18) ??=(??+?????????????????????????)2????????+(??+?????????????????????????)?????????????2+13?????2????????3?+(??+?????????????????????????)2????????+(??+?????????????????????????)?????????????2+13?????2????????3?+(???????????????+?????????)2?????????(???????????????+?????????)?????????????2+13?????2????????3?+(?????????????????????????)2?????????(?????????????????????????)?????????????2+13?????2????????3
系统的粘性耗散函数 D 可描述为：
(19) ??=12?4∑??=1∫??0˙?????(??)T??????˙?????(??)????
其中 Cj 是第 j 排支撑线的单位长度阻尼矩阵，其表达式为 ????=diag?(0,???????,???????)。基于均匀性假设，每排支撑线的阻尼分布可得到：ψyj = ψy, ψzj = ψz。
因此，系统的最终耗散能表示为：
(20) ??=(˙??+????????˙?????????????˙????)2????????+(˙??+????????˙?????????????˙????)?˙????????????2+13?˙????2????????3?+(˙??+????????˙?????????????˙????)2????????+(˙??+????????˙?????????????˙????)?˙????????????2+13?˙????2????????3?+(˙???????????˙????+?????˙????)2?????????(˙???????????˙????+?????˙????)?˙????????????2+13?˙????2????????3?+(˙???????????˙??????????˙????)2?????????(˙???????????˙??????????˙????)?˙????????????2+13?˙????2????????3
本研究中施加到系统上的外部激振力为 ??=[???????,???????,???????,???????,???????]T，如图 6 所示，其分量的分解结果如下：
(21) ?{ { { {?{ { { {????????=?????cos?45°???????=?????sin?45°???????=?????cos?45°·??1??????sin?45°·??2???????=?????sin?45°·??3???????=??????cos?45°·??3
下载：下载高分辨率图像 (434KB)
下载：下载全尺寸图像
图 6. 外部激振分解示意图。
当五自由度耦合激振施加到直线导轨时，其运动微分方程为：
(22) ???¨??+???˙??+?????=????
其中：
??=?? ? ? ? ???????????????????? ? ? ? ??
????=?? ? ? ? ? ???????????????????????? ? ? ? ? ??
??=?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??4????????0?2?????????(???????+???????)02?????????(2????????+??)04????????0?2?(2????????+??)????????0?2?(???????+???????)????????02?[(??2?????+??2?????)?????+2???2???????]???0?(???????+???????)?(2????????+??)????????0?2?(2????????+??)????????04?(??2?????+??????????+13???2)????????02?(2????????+??)????????0?(???????+???????)?(2????????+??)????????04?(??2?????+??????????+13???2)?????????? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??
??=?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??4????????0?2?(???????+???????)????????02?(2????????+??)????????04????????0?2?(2????????+??)????????0?2?(???????+???????)????????02?[(??2?????+??2?????)?????+2???2???????]???0?(???????+???????)?(2????????+??)????????0?2?(2????????+??)????????04?(??2?????+??????????+13???2)????????02?(2????????+??)????????0?(???????+???????)?(2????????+??)????????04?(??2?????+??????????+13???2)?????????? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??

2.2. 参数辨识方法
这里应用纽马克 -β 方法对公式 (22) 进行离散化。在时间区间 [t, t + Δt] 内，假设：
(23) ˙????+1=˙????+(1???)?Δ????¨????+???Δ????¨????+1
(24) ????+1=????+Δ????˙????+(12???)?Δ???2?¨????+???Δ???2?¨????+1
(25) ¨????+1=1???Δ???2?(????+1?????)?1???Δ????˙?????(12????1)?¨????
(26) ˙????+1=?????Δ????(????+1?????)+(1?????)?˙????+(1???2???)?Δ????¨????
其中 ??≥1/2, ??≥1/4?(1/2+??)2。这里取 ??=1/2 和 ??=1/4。
将公式 (25) 和 (26) 代入公式 (22) 并整理得到：
(27) ?{ { { { {?{ { { { {???(??)11?????+??(??)13?????=??(??)1??(??)22?????+??(??)24?????=??(??)2??(??)31?????+??(??)32?????+??(??)33?????+??(??)34?????=??(??)3??(??)42?????+??(??)44?????=??(??)4??(??)51?????+??(??)53?????=??(??)5
公式 (27) 中系数的表达式见附录 A。
公式 (27) 可以表示为以下矩阵形式：
(28) ??(??)???=??(??)
其中：
(29) ??(??)=?? ? ? ? ? ? ????(??)110??(??)1300??(??)220??(??)24??(??)31??(??)32??(??)33??(??)340??(??)420??(??)44??(??)510??(??)530?? ? ? ? ? ? ??
??=?? ? ? ???????????????????? ? ? ??
??(??)=?? ? ? ? ? ? ????(??)1??(??)2??(??)3??(??)4??(??)5?? ? ? ? ? ? ??
通过对所有时间点 i = 1, 2, …n 的方程进行积分，得到：
(30) ?????=??
其中：??=[??(1),??(2),…,??(??)]T, ??=[??(1),??(2),…,??(??)]T。
向量 η 利用最小二乘法估计，其解析解可得到为：
(31) ??=(??T???)?1???T???

3. 滑座姿态确定
3.1. 加速度信号积分
3.1.1. 频域积分
传统时域数值积分方法通过数值计算直接对加速度信号进行积分。虽然计算简单，但它们本质上均匀处理信号的整个频段，难以有效抑制低频漂移和高频噪声 [41,42]。因此，本研究采用频域积分方法处理数据。该方法的主要优势在于能够充分利用频域分析的灵活性。通过将加速度信号转换到频域，可以对不同频率分量进行差异化处理 [43]。
加速度振动信号的离散值为 {???(??)}?(??=0,1,2,?,??)。单次积分的数值公式为：
(32) ???(??)=???1∑??=01???2???????Δ???????(??)????(??)??????2?????????/??
二次积分的数值公式为：
(33) ???(??)=???1∑??=0?1(2???????Δ???)2????(??)????(??)??????2?????????/??
其中：
(34) ???(??)={1(????≤???Δ???≤????)0
其中 fd 和 fu 分别是截止频率的下限和上限；A(k) 是 a(r) 的傅里叶变换；Δ?? 是频率分辨率。

3.1.2. 初始位移处理
通过频域积分获得的位移信号已具有良好的精度和稳定性。然而，它仍然存在一个显著问题：初始位移不为零，这在物理上是不合理的。为了纠正这个问题，通过采用多节点三次样条插值方法进一步优化了处理过程。
多节点三次样条插值在连续性、光滑度和局部形状保持能力方面具有优势 [44,45]。与线性插值或多项式拟合方法相比，该方法可以确保相邻节点间一阶和二阶导数的连续性，从而能够更准确地描述振动信号中常出现的非线性和非平稳动态特征。同时，它可以避免高阶多项式可能出现的振荡现象 [46]。通过可调整的节点放置，该方法可以灵活适应局部信号波动，同时保持整体光滑度，从而提高从加速度重构位移积分路径的准确性。因此，多节点三次样条插值特别适合涉及分段平滑特征的工程信号处理任务，如结构振动响应。
具体实施如下：为了在时间 t0 强制执行零位移基准，同时确保向第一个显著峰值的平滑自然过渡，我们在区间 [t0, tp] 上构建多节点三次样条插值。插值节点 τj 选择为该区间内的 M 个等间距时间点，其对应值取自原始位移 ????=???(τ??)，其中第一个节点值明确设为零。插值函数 S(t) 在每个子区间 [????,????+1] 上由三次多项式表示，满足插值条件并保持一阶和二阶导数的连续性。
具体步骤如下：
(1). 区间和节点定义
过渡区间定义为：
(35) [??,??]=[??0,????]
在区间 [a, b] 内均匀选择 M 个节点，时间点由下式给出：
(36) ????=??0+???1???1?(???????0), ??=1,…,??
这些节点处的函数值从原始位移 ???(??) 中采样，其中第一个节点明确设为零：
(37) ??1=0, ????=???(????), ??=2,…,??

(2). 三次样条构造
在区间 [a, b] 上构建分段三次样条插值函数 S(t)。在每个子区间 [tj, tj+1] 上，样条表示为：
(38) ?????(??)=????+?????(???????)+?????(???????)2+?????(???????)3, ??=1,…,???1
每个多项式段的系数由以下条件确定：
1). 插值条件
(39) ?????(????)=????, ?????(????+1)=????+1, ??=1,…,???1
2). 一阶导数连续性
(40) ??′???(????+1)=??′??+1?(????+1), ??=1,…,???2
3). 二阶导数连续性
(41) ??″???(????+1)=??″??+1?(????+1), ??=1,…,???2
4). 边界条件
应用 MATLAB 中实现的非 knot 边界条件：
(42) ???1?(??2)=???2?(??2), ??????2?(?????1)=??????1?(?????1)
通过满足上述条件，得到 S(t)，它在区间 [a, b] 上二阶连续可微。

(3). 调整后位移
由此得到最终调整后的位移：
(43) ??adj?(??)=?{ {?{ {?0,??????

3.2. 滑座质心位移计算
为了确定滑座系统的质心位移，本研究采用三轴加速度传感器在滑座上的多个位置测量动态信号。然后通过最小二乘法进行反分析来计算质心的位移响应。具体测量装置如图 7 所示。在滑座表面共建立了 5 个测点。每个点都配备了一个三轴加速度传感器，以同时捕获三个正交方向上的加速度信号。这些测点的空间坐标详见表 2。基于测点 1 - 4 的测量数据和空间坐标信息，采用最小二乘拟合法重建滑座运动过程中质心位移的时间历程曲线。最终使用测点 5 的测量结果验证重建精度。
下载：下载高分辨率图像 (347KB)
下载：下载全尺寸图像
图 7. 测点位置示意图。
表 2. 惯性坐标系中测点的坐标。

| 测点 | 坐标 | 数值 (mm) |
| :--- | :--- | :--- |
| 1 | ltx1 × lty1 × ltz | (35, ?75, ?38.76) |
| 2 | ltx2 × lty2 × ltz | (35, 75, ?38.76) |
| 3 | ltx3 × lty3 × ltz | (-35, 75, ?38.76) |
| 4 | ltx4 × lty4 × ltz | (-35, ?75, ?38.76) |
| 5 | ltx5 × lty5 × ltz | (0, 0, ?38.76) |

当系统处于静平衡位置（即无振动）时，测点 1 到 5 在坐标系 OXYZ 和 Gx"y"z" 中重合：
(44) ????=??″??=[??????????,??????????,???????]T, ??=1,2,3,4,5
当系统运行时，测点 1 到 5 在参考系 Gx"y"z" 中的坐标 ??″?? (i = 1, 2, …,5) 保持不变。它们在坐标系 Gx'y'z' 和 OXYZ 中的坐标值 ??′?? 和 ????? 分别如下：
(45) ??′??=??T???″??
(46) ?????=??+??T???″??
其中：
(47) ?????=?? ??0?????? ??+?? ? ??1?????????????1??????????????1?? ? ????? ??????????????????????????????? ??=?? ? ????????????????????????????+??????????????+???????????????+?????????????????????????????????????????+???????????????+????????? ? ??
运行期间，每个测点的位移可表示为：
(48) ????=??????????=[???????,???????,???????]T
此外，每个测点的位移可以用矩阵形式表示为：
(49) ????=?? ??0?????? ??+?? ? ??0?????????????0??????????????0?? ? ????? ??????????????????????????????? ???=?? ? ??000??????????????????10????????0??????????01?????????????????????0?? ? ????? ? ? ? ???????????????????? ? ? ? ???=?????????
因此，所有测点的位移可以统一表示为：
(50) ??=???????
其中 ??=[??1,??2,??3,??4]T, ??=[??1,??2,??3,??4]T。
残差方程矩阵表达式为：
(51) ??=??????????
最小二乘准则是最小化：
(52) min?????‖??‖2=(??????????)T?(??????????)
通过对 XG 求导，得到以下正规方程矩阵形式：
(53) ??T??????=??T???
最终的质心位移为：
(54) ????=???1???T???
其中 ??=??T???。
用于验证质心位移精度的测点 5 的位移变化，在振动期间满足以下方程：
(55) ?{ {?{ {??????5=?????????5?????+?????????????????5=??+????????5???????????????????????5=???????????5?????+????????5?????

4. 仿真验证
为了验证所提出的辨识算法的数值实现和参数恢复能力，首先预设一组单位长度刚度和阻尼参数，并用于计算直线导轨系统的动态响应。随后将模拟响应视为观测数据，引入辨识过程以恢复预设参数。该数值研究在匹配模型条件下对参数辨识能力提供了初步验证，其中响应生成模型和辨识模型是一致的。响应生成模型和辨识模型之间由于理想化模型和实际系统之间的差异而产生的模型形式不确定性的影响，将在第 6.3 节中通过引入响应生成模型和辨识模型之间的非对称支撑条件进一步研究。
在仿真验证中，采用高斯脉冲作为激振力，其数学表达式如下：
(56) ???(??)=??·exp?(?(?????0)2/(2???2))
其中 F(t) 是脉冲幅度；t0 是峰值时间，即脉冲达到最大值 A 的时间；σ 是高斯函数的标准差，控制脉冲的时间宽度。

在本节所展示的仿真中，参数设置为A=500N，t?=0.005s，σ=0.0001s，激励力的施加位置如图6所示。本研究使用的线性导轨的基本参数详见表3。指定的单位长度刚度和单位长度阻尼值见表4。

表3. 线性导轨的基本参数。

参数 | 参数值
---|---
滑块和金属块的质量 m | 16.35 kg
ly | 17 mm
绕X轴的转动惯量 Ix | 0.08 kg·m2
lax | ?52.68 mm
绕Y轴的转动惯量 Iy | 0.07 kg·m2
lbx | 52.68 mm
绕Z轴的转动惯量 Iz | 0.03 kg·m2
luz×ldz | 37.99×52.48 mm
每列滚珠长度 L | 105.36 mm
l?×l?×l? | 33.76×110×50 mm

表4. 指定参数值。

参数 | 参数值
---|---
ψy (N·s/m2) | 10000
ψz (N·s/m2) | 10000
ζy (N/m2) | 500000000
ζz (N/m2) | 500000000

然而，在实际实验中，获取的信号不可避免地受到周围环境中各种内部和外部干扰的污染。由于所提出的识别框架使用位移响应作为观测数据，本节直接在位移信号中加入高斯白噪声来模拟测量噪声，并评估所提出框架对噪声污染的鲁棒性。含噪声的位移信号可表示为：

(57) $x_n(t) = x(t) + n(t)$

其中x(t)为无噪声位移，n(t)为均值为零的高斯白噪声。

噪声水平通过信噪比（SNR）来控制，其定义为：

(58) $SNR = 10\log_{10}\left(\frac{P_{signal}}{P_{noise}}\right) = 20\log_{10}\left(\frac{A_{signal}}{A_{noise}}\right)$

其中$P_{signal}$表示信号功率，$P_{noise}$为噪声功率，$A_{signal}$为信号幅度，$A_{noise}$为对应的噪声幅度。

本节选取了三种典型噪声水平进行分析，分别为SNR=20 dB、15 dB和10 dB，分别对应轻、中、重噪声条件，以全面评估识别方法在不同干扰水平下的性能。含噪声与无噪声位移信号的时间域对比结果如图8所示。

下载：下载高分辨率图片（2MB）| 下载：下载完整尺寸图片

图8. 不同信噪比水平下的位移响应对比：(a) Y方向；(b) Z方向；(c) θx方向；(d) θy方向；(e) θz方向。

将动态响应结果代入公式(31)，可以识别出线性导轨的单位长度刚度和阻尼参数，具体数值列于表5中。从表6中列出的识别结果的相对误差可以看出，所提出的方法在低至中等信噪比水平下能有效识别结构的刚度和阻尼参数，所有相对误差均控制在5%以内。然而，在高噪声条件下，ψy和ζz的识别结果不尽如人意，这表明尽管该方法具备一定的抗噪声能力和工程适用性，但其在强噪声环境下的识别性能仍有待改进。因此，如果要将采集到的位移信号直接用于识别，建议先进行降噪预处理。

表5. 参数识别结果。

参数 | 无噪声 | SNR=20 dB | SNR=15 dB | SNR=10 dB
---|---|---|---|---
ψy (N·s/m2) | 9898.77 | 10544.54 | 9815.68 | 11374.61
ψz (N·s/m2) | 10071.76 | 10546.07 | 9710.64 | 9348.19
ζy (N/m2) | 495229821.83 | 495957316.31 | 514244210.93 | 532383452.06
ζz (N/m2) | 500790825.96 | 508829991.85 | 521410386.79 | 595867394.81

表6. 参数识别结果的相对误差。

参数 | 无噪声 | SNR=20 dB | SNR=15 dB | SNR=10 dB
---|---|---|---|---
ψy (N·s/m2) | 1.012% | 5.45% | 1.84% | 13.75%
ψz (N·s/m2) | 0.718% | 5.46% | 2.89% | 6.52%
ζy (N/m2) | 0.954% | 0.81% | 2.85% | 6.48%
ζz (N/m2) | 0.158% | 1.77% | 4.28% | 19.17%

为进一步验证所识别参数的有效性，将理论位移曲线与无噪声条件下的仿真位移曲线进行了对比。五个自由度（即y、z、θx、θy、θz）的对比结果如图9所示。从曲线对比中可以看出，理论位移与仿真位移吻合良好。

下载：下载高分辨率图片（1MB）| 下载：下载完整尺寸图片

图9. 理论质心位移与仿真质心位移对比。

为了更全面、系统地评估所提出的识别方法在动态参数逆估计中的准确性和可靠性，本研究引入了四个互补的误差指标：RMSE和相对误差用于量化绝对误差和归一化幅值偏差，MAE用于表征计算结果与实验结果之间的平均绝对偏差，从而反映模型的整体误差水平，R2用于评估趋势一致性。

均方根误差（RMSE）反映了整个时间段内计算值与测量值之间误差的平均波动范围，其表达式如文献[47]所示：

(59) $RMSE = \sqrt{\frac{1}{n}\sum_{k=1}^{n}(x(t_k) - x_s(t_k))^2}$

其中$x(t_k)$表示测量值，$x_s(t_k)$表示计算值。RMSE的最小值为零，当所有计算值与对应的测量值完全匹配时出现。

为进一步评估误差相对于幅值的相对大小，可以引入以峰值$x_{peak}$为参考的相对误差，其表达式为：

(60) $Relative\ Error = \frac{RMSE}{x_{peak}} \times 100\%$，$x_{peak} = \max_k |x(t_k)|$

然而，RMSE对数据量级敏感，仅依赖此指标可能无法全面评估模型的预测精度。因此，有必要引入额外的误差指标以实现更全面的模型性能评估。平均绝对误差（MAE）评估预测值与真实值之间的平均绝对偏差，提供了一个客观的整体预测误差度量，不受正负误差相互抵消的影响。因此，本研究采用MAE作为补充评估指标，其表达式如下：

(61) $MAE = \frac{1}{n}\sum_{k=1}^{n}|x(t_k) - x_s(t_k)|$

决定系数R2用于表征计算值与测量值之间的一致性。其取值范围通常为[0, 1]，其表达式如下：

(62) $R^2 = 1 - \frac{\sum_{k=1}^{n}(x(t_k) - x_s(t_k))^2}{\sum_{k=1}^{n}(x(t_k) - \bar{x})^2}$

其中$\bar{x}$表示$x(t_k)$的均值。当R2趋近于1时，表明计算结果非常精确；反之，R2值越低则表明计算性能越差。

图10展示了五个方向的位移对比结果的误差评估指标。从图中可以观察到，在无噪声条件下，RMSE、相对误差和MAE的值均相对较小，R2值非常接近1，表明理论位移与仿真位移之间的偏差极小。随着SNR从20 dB降至10 dB，噪声水平逐渐增大，理论位移与仿真位移之间的误差呈现递增趋势。具体而言，当SNR=20 dB时，各方向的位移相对误差范围约为0.1%至1.2%，R2值非常接近1，表明在轻微噪声条件下理论位移与仿真位移吻合良好。

下载：下载高分辨率图片（1MB）| 下载：下载完整尺寸图片

图10. 解算精度和动态一致性评估指标示意图：(a) Y方向；(b) Z方向；(c) θx方向；(d) θy方向；(e) θz方向。

当SNR=15 dB时，θy方向的误差略有增加，但总体上仍保持在5%以内，R2值范围为0.87至0.98，表明该方法在中等噪声条件下仍保持较高的识别精度。然而，当SNR进一步降低至10 dB时，某些方向（如θy）的相对误差增加到17.545%，反映出强噪声干扰对特定参数识别精度有显著影响。在所有噪声水平下，RMSE和MAE值均保持较小，表明尽管噪声水平增加导致识别精度下降，但整体误差仍处于工程应用可接受的范围。

综上所述，所提出的方法在低至中等噪声条件（SNR≥15 dB）下表现出优异的鲁棒性，即使在强噪声环境下仍能提供合理的参数估计，显示出在工程应用中具有良好的前景。

## 5. 实验结果

### 5.1. 质心位移解算结果与检测结果

被测线性导轨为THK HSR35型号，预紧力为623 N。在本研究所采用的激励水平下，该预紧条件支持滚珠列与滚道之间连续接触的假设。因此，滚珠与滚道的接触行为在预紧平衡状态附近局部近似为等效线性和刚度与阻尼模型。

实验验证方案如图11所示。脉冲激励由PCB 086C01冲击锤施加，该冲击锤内置压电传感器，能够同步测量冲击力信号。试样的振动响应通过三轴加速度传感器采集。力信号和加速度信号由同一数据采集系统同步采集和记录，确保数据的时间域对齐。三轴加速度传感器的规格列于表7中。

下载：下载高分辨率图片（565KB）| 下载：下载完整尺寸图片

图11. 冲击锤测试示意图。

表7. 加速度传感器系统规格。

传感器 | 灵敏度(X/Y/Z) [mV/g] | 输入电压[V] | 带宽[kHz] | 采样频率[kHz]
---|---|---|---|---
1 | 19.08/19.47/20.17 | 10 | 51.2 | 8
2 | 20.08/20.09/19.96 | 10 | 51.2 | 8
3 | 20.28/19.52/20.24 | 10 | 51.2 | 8
4 | 19.08/19.47/20.17 | 10 | 51.2 | 8
5 | 20.19/20.37/20.38 | 10 | 51.2 | 8

三轴加速度传感器以12.8 kHz的采样频率采集测点1至5的时间域加速度信号，每次冲击的录制时长为0.5 s。对应的波形如图12所示。通过对所有测点采集的加速度信号进行快速傅里叶变换（FFT），获得了图13所示的频率谱分布。频谱分析清晰地揭示了各测点振动信号在频域中的能量集中特性和分布规律，为后续通过频域积分方法计算速度和位移响应提供了关键的频率结构信息。

下载：下载高分辨率图片（979KB）| 下载：下载完整尺寸图片

图12. 实验加速度曲线：(a)-(c) 测点1的X、Y、Z方向加速度；(d)-(f) 测点2的X、Y、Z方向加速度；(g)-(i) 测点3的X、Y、Z方向加速度；(j)-(l) 测点4的X、Y、Z方向加速度；(m)-(o) 测点5的X、Y、Z方向加速度。

下载：下载高分辨率图片（1MB）| 下载：下载完整尺寸图片

图13. 加速度频谱：(a) 测点1；(b) 测点2；(c) 测点3；(d) 测点4；(e) 测点5。

在频域中，使用第3节推导的频域积分公式进行二次积分。为抑制低频漂移和高频噪声干扰，施加了截止频率为80 Hz和4000 Hz的带通截断，去除低频成分以减轻趋势引起的漂移，并衰减高频成分以减少噪声影响。完成积分后，通过逆傅里叶变换将结果转换回时域。为处理积分初期可能出现的非物理初始偏移，采用多节点三次样条插值方法进行修正。选取振动起始时刻t?=0.0064 s作为基线重置点，将该时刻之前的位移统一设为零。在t?到第一个峰值的区间内选取四个等距节点作为约束，构建三次样条插值函数以重构位移序列。该方法有效消除了初始漂移，同时最大程度地保留了原始位移曲线的波形特征。通过上述处理，最终获得了各测点完整可靠的位移时程曲线，如图14所示。

下载：下载高分辨率图片（1MB）| 下载：下载完整尺寸图片

图14. 积分位移曲线：(a)-(c) 测点1的X、Y、Z方向位移曲线；(d)-(f) 测点2的X、Y、Z方向位移曲线；(g)-(i) 测点3的X、Y、Z方向位移曲线；(j)-(l) 测点4的X、Y、Z方向位移曲线；(m)-(o) 测点5的X、Y、Z方向位移曲线。

所获得的位移数据不仅为后续的质心运动分析提供了直接输入，其处理过程本身也验证了该信号处理方法在多测点振动测试中的适用性。图14中五个测点的位移曲线表现出良好的规律性和协调性，反映了激励下滑块系统的整体运动特征，从而为后续动态参数识别和模型验证奠定了数据基础。

通过将图14中测点1-4的位移数据代入第3节建立的质心位移解算公式(54)，可以反解出系统质心在五个自由度上的位移响应，结果如图15所示。该方法充分利用了多点位移信息，通过数学建模将局部观测数据与整体运动状态相关联，有效重构了质心的完整运动学行为，实现了从局部观测到系统整体运动特征的数学转换。该过程不仅提高了质心运动状态估计的精度和鲁棒性，也为后续动态分析和系统状态评估提供了关键的底层信息。

下载：下载高分辨率图片（808KB）| 下载：下载完整尺寸图片

图15.

质心五方向位移曲线

为验证所提出的质心位移求解方法的准确性，选取图7中的测点5作为独立验证基准。利用第3节中的公式(55)将质心位移响应反投影至该测点位置，获得理论位移曲线，随后将该理论曲线与测点5的实际测量位移数据进行对比分析。图16中的对比结果表明，理论曲线与实测数据在时程趋势、相位对应关系和幅值范围方面表现出高度一致性。这一结果充分验证了本文所提出的质心位移求解方法在数学上的合理性和工程上的可靠性，为后续基于质心运动的动态参数识别提供了有效的数据基础。

在参数识别之前，需验证线性导轨系统在冲击激励下的线性特性，以确保线性动态模型的适用性。为验证线性导轨系统在冲击激励下的线性特性，采用三种不同的冲击力水平（272 N、356 N和528 N）进行了对比实验。对每种工况下的响应信号进行频谱分析，结果如图17所示。如图17所示，三种工况下的主导频率相同，随着激励幅值的变化，没有出现频率偏移或新的频率成分出现。这表明系统的固有特性不随激励水平的变化而改变，符合线性时不变系统的行为特征。在此基础上，可以得出结论：所研究的线性导轨系统在施加的冲击范围内保持线性运行。

需要注意，所提出的线性化模型的有效性取决于线性导轨的预紧状态。当预紧力足够大，能够保持受载滚珠列与滚道在激励过程中持续接触时，所提出的方法适用。在预紧力不足、预紧损失、外载荷过大或大幅振动导致部分接触丧失的条件下，滚珠与滚道之间的相互作用可能表现出显著的非线性特征，此时应考虑采用非线性接触模型。

基于Newmark-β方法对动态方程进行离散化，根据公式(31)利用最小二乘法求解系统方程，识别线性导轨的单位长度刚度和阻尼参数。最终的识别结果汇总于表8中。

表8列出了系统参数的识别值及其相应的不确定性指标，包括标准差、95%置信区间和变异系数。识别得到的阻尼系数ψy和ψz分别为20030.3 N·s/m2和18947.1 N·s/m2，对应的变异系数分别为5.84%和5.80%。估算得到的刚度参数ζy和ζz分别为387248829.8 N/m2和990193737.3 N/m2，对应的变异系数分别为2.96%和1.43%，低于阻尼系数。

不确定性分析表明，刚度参数的变异系数均低于3%，显示出较强的统计稳定性和较高的工程应用可靠性。相比之下，阻尼参数的变异系数约为5.8%，略高于刚度参数，但仍处于可接受的工程范围内。这一差异主要源于阻尼项对系统动态响应的贡献相对较小，使得识别得到的阻尼参数对测量噪声和建模误差更为敏感。

总体而言，所有识别参数的标准差和95%置信区间均表明具有良好的统计一致性和较高的估计精度，证实了所提出的识别方法在工程应用中既有效又可靠。

5.2 重复冲击试验的可重复性分析

尽管单次瞬态响应中包含足够的系统动态信息用于参数识别，但在相同的实验条件下进行了三次独立的重复冲击试验，以进一步验证所提方法的鲁棒性和可重复性。试验过程中，冲击位置、激励方向、传感器布置和采样参数保持不变，以最大限度减少外部实验条件对识别结果的影响。每次冲击试验获得的响应均独立地利用所提出的频域双重积分和质心位移重建流程进行处理。三次重复冲击试验所获得的相应冲击力和重建的质心位移响应如图18所示。可以观察到，不同试验中的冲击载荷曲线和重建响应表现出相似的时间特性，尽管由于手动冲击激励不可避免的差异，幅值存在轻微变化。这证实了实验过程的一致性，并为后续的参数识别提供了可靠的响应数据。

在此基础上，独立地对三次重复冲击试验的响应数据进行动态参数识别。以第5.1节中原始冲击试验获得的识别结果作为参考基准，计算三次重复试验识别参数的相对偏差。详细的识别结果列于表9中，相应的与参考结果的相对偏差汇总于表10中。可以观察到，所有识别参数与参考值之间的偏差均较小。所有参数中的最大相对偏差约为3%，表明所提出的识别方法在重复冲击条件下具有良好的可重复性。

5.3 实验验证

将公式(25)和(26)代入公式(22)，可得qt+1的方程如下：

(63) K? qt+1 = R?t+1

其中：

(64) K? = K + (γ/(βΔt)) C + (1/(βΔt2)) M

(65) R?t+1 = Fp + M[(1/(βΔt2)) qt + (1/(βΔt)) q?t + (1/(2β) ? 1) q?t] + C[(γ/(βΔt)) qt + (γ/β ? 1) q?t + Δt(γ/(2β) ? 1) q?t]

根据公式(25)和(26)，可分别求解q?t+1和q?t+1。

将识别得到的刚度和阻尼参数代入动态微分方程中，并沿Y轴正方向和Z轴正方向对系统施加外部激励。

第一次冲击试验中，实验装置如图19(a)所示。Y方向各测点的加速度响应信号如图19(b)所示，施加的外部激励时程曲线如图19(c)所示。通过求解动态微分方程获得理论加速度响应。图20对比了测点5处的理论加速度与实测加速度信号，结果表明两者在时域波形、幅值和相位方面均表现出良好的一致性。此外，图21(a)和(b)分别对比了速度和位移响应，进一步验证了识别参数的准确性。

表11列出了加速度、速度和位移响应的误差评价指标。由表可知，加速度、速度和位移响应的相对误差分别为2.665%、4.450%和3.037%，均在4.5%以内。对应的R2值分别达到0.751、0.796和0.948，RMSE和MAE均保持在较低水平，表明理论计算结果与实验结果具有良好的相关性。

此外，对Y方向的加速度响应信号进行了频域分析，仿真结果与实验结果的频谱对比如图22所示。从图22可以看出，在Y方向，系统的主导频率成分主要集中在144 Hz附近。仿真结果与实验结果在主导频率位置方面吻合良好，整体频谱分布趋势基本一致，表明该模型能够准确捕捉该方向的主导模态特性。然而，在高频段，仿真与实验频谱之间仍存在一定差异，这表明所建立的5自由度模型在表征系统高阶模态和高频动态行为方面存在一定局限性。

第二次冲击试验中，实验装置如图23(a)所示。对应的Z方向加速度响应信号如图23(b)所示，外部激励时程曲线如图23(c)所示。图24对比了测点5处的理论数据与实测数据，加速度、速度和位移曲线在波形、幅值和相位方面均表现出高度一致性，充分验证了识别参数的准确性。

表12汇总了加速度、速度和位移响应的误差评价结果。数据显示，加速度、速度和位移的相对误差分别为2.365%、2.093%和2.090%，均低于2.5%。对应的决定系数R2分别为0.783、0.837和0.827，RMSE和MAE均保持在较低水平。这些结果表明，理论计算与实验测量结果吻合良好。

随后，对Z方向的加速度响应信号进行了频域分析，仿真结果与实验结果的频谱对比如图25所示。频谱分析结果表明，在Z方向，系统的主导振动频率主要集中在920 Hz附近。仿真结果准确预测了这一主导频率位置，并在主频范围内与实验结果保持良好的一致性，有效反映了系统的主要动态行为。然而，在高频成分处，两者之间仍存在一定偏差，这表明当前模型难以充分描述导轨系统的高频动态响应。

综上所述，所建立的5自由度动态模型不仅能够以合理的精度描述系统的时间域响应特性，还能有效预测Y和Z方向的主导共振频率，表明该模型能够较好地反映系统的低频主导动态特性。由于模型中进行的合理简化，在预测高阶模态和高频动态行为方面仍存在一些误差。因此，该模型更适合用于分析导轨系统的低频主导动态特性，而高频复杂动态行为的研究则需要开发更加精细的动态模型。

需要指出的是，本文所采用的等效线性刚度和阻尼假设与线性导轨的预紧状态密切相关。所提出的方法适用于能够保持受载滚珠列与滚道在激励过程中持续接触的预紧力水平下。在预紧力不足、激励幅值过大或出现部分接触丧失的情况下，滚珠与滚道之间的相互作用可能表现出显著的非线性特征，此时线性化模型可能无法准确表征系统的动态行为。

6. 讨论
6.1 识别参数的敏感性分析

为研究识别参数对冲击响应特性的影响，基于所建立的5自由度动态模型进行敏感性分析。采用高斯脉冲载荷来模拟施加在滑块上的瞬态冲击激励。

响应指标定义为各自由度瞬态位移响应的第一个主峰，以排除后续自由振动和振铃的影响，从而为冲击载荷下的直接动态响应提供更准确的表征。图26展示了系统参数在20%扰动下五个自由度（y、z、θx、θy、θz）峰值响应的相对变化。总体而言，结果表明响应对刚度参数（ζy、ζz）的变化比对阻尼参数（ψy、ψz）的变化更为敏感。具体而言，ζy的变化会导致y、θx和θz方向峰值响应出现约10%的变化。这表明Y方向相关的刚度及其耦合刚度项在控制系统的瞬态动态响应中起主导作用。类似地，ζz的变化对应z和θy方向峰值响应约8%–9%的变化，表明横向刚度主要控制相应耦合振动模态的动态特性。下载：下载高清图片（477KB）下载：下载全尺寸图片图26. 辨识参数的敏感性分析。相比之下，ψy和ψz的变化仅导致峰值响应出现相对较小的变化，量级在0.5%–1.5%之间。这表明在当前低阻尼条件下，阻尼变化主要影响能量耗散过程，而对瞬态峰值响应幅值的影响有限。这一观察与冲击载荷下系统的一般行为一致，即峰值响应主要由惯性与刚度之间的相互作用决定，而阻尼主要影响后续振动的衰减特性。总之，在冲击载荷条件下，所考虑的结构表现出以刚度为主导的动态行为。瞬态峰值响应主要随刚度相关参数的变化而变化，而阻尼参数的变化则导致相对较小的变化。这些发现为后续的结构优化和参数整定提供了定量依据。6.2 三次样条插值对辨识结果的影响为定量评估三次样条插值处理对参数辨识结果的影响，本节仅去除三次样条插值修正步骤，而保持所有其他处理程序不变。基于相同的原始加速度信号和频域积分框架，重新进行位移重建和参数辨识，以获取对照组的辨识结果。图27展示了未经三次样条插值处理时五个测点的重建位移和质心位移。下载：下载高清图片（2MB）下载：下载全尺寸图片图27. 测点1–5和重建质心的位移曲线。在未进行三次样条插值的情况下，最终辨识的参数为：ψy = 38101.5 N·s/m2，ψz = 22819.5 N·s/m2，ζy = 471978829.9 N/m2，ζz = 1018749829.9 N/m2。将上述辨识参数代入系统动力学方程求解系统响应。对两种条件（有插值处理和无插值处理）下辨识结果所得的加速度响应误差进行比较分析，如图28所示。下载：下载高清图片（742KB）下载：下载全尺寸图片图28. 解算精度和动态一致性评价指标示意图。比较结果表明，在对初始位移未应用三次样条插值修正的情况下，通过积分得到的位移结果违反了物理约束，从而导致参数辨识结果出现显著偏差。因此，在完成频域积分后，必须在参数辨识之前对重建位移进行初始位移修正，以确保辨识结果的准确性和可靠性。6.3 支撑不对称性和模型差异对辨识结果的影响为评估第2节中采用的全对称假设对参数辨识精度的影响，进一步考虑了线性导轨压缩侧和拉伸侧之间的刚度和阻尼不对称性。具体而言，在保持拉伸侧参数不变的情况下，将压缩侧的刚度和阻尼参数分别增加5%、10%和20%，以模拟实际运行中因载荷变化导致接触角变化而引起的参数偏差。对于每种不对称条件，首先基于非全对称的五自由度动力学方程计算系统动态响应，并将其作为参考响应。随后，将所得响应代入基于全对称五自由度模型建立的参数辨识流程中，以获得相应的等效对称刚度和阻尼参数。最后，通过比较重建响应与参考响应来评估全对称假设对模型预测能力的影响。本节仅考虑三种典型不对称条件，即5%、10%和20%。因此，与第4节的匹配模型数值验证不同，本分析中响应生成模型和辨识模型被有意设置为不同，从而可以直接评估模型差异对辨识结果的影响。不同不对称条件下压缩侧和拉伸侧的单位长度刚度和阻尼参数列于表13。拉伸侧参数与标称辨识结果保持一致，而压缩侧参数按照规定比例增加。表13. 不同不对称条件下线性导轨的参数设置。参数 | 5%不对称条件 | 10%不对称条件 | 20%不对称条件 | ψy1, ψy2, ψy3, ψy4 | 20030.3 | 20030.3 | 20030.3 | ψz1, ψz2 | 19894.5 | 20841.8 | 22736.5 | ψz3, ψz4 | 18947.1 | 18947.1 | 18947.1 | ζy1, ζy2, ζy3, ζy4 | 387248829.8 | 387248829.8 | 387248829.8 | ζz1, ζz2 | 1039703424.2 | 1089213111.0 | 1188232484.8 | ζz3, ζz4 | 990193737.3 | 990193737.3 | 990193737.3 | 如表14中的辨识结果所示，从全对称模型获得的等效参数随着压缩侧参数的增加呈现出明显的变化趋势。具体而言，Z方向参数ψz和ζz随着不对称度的增加逐渐增大，而Y方向参数ψy和ζy则几乎保持不变。表14. 不同不对称条件下辨识的等效参数。参数 | 5%不对称条件 | 10%不对称条件 | 20%不对称条件 | ψy (N·s/m2) | 20472.5 | 20472.7 | 20473.4 | ψz (N·s/m2) | 20241.5 | 20754.9 | 21784.6 | ζy (N/m2) | 378514504.0 | 378506083.6 | 378487875.3 | ζz (N/m2) | 976095909.2 | 999069652.0 | 1044917436.0 | 这种行为在物理上是合理的。由于不对称扰动仅引入在Z方向，Y方向的动力学方程保持原有的对称特性，导致相应辨识参数的变化可忽略不计。然而对于Z方向，实际系统包含压缩侧和拉伸侧之间的参数差异。因此，基于全对称模型的辨识过程确定了一组在对称约束下的等效参数来近似整体动态响应。因此，辨识得到的ψz和ζz值并不代表任何一侧的实际参数，而是表征两侧的综合动态行为，故随压缩侧参数的增加而增大。不同不对称条件下的响应预测误差如图29所示。总体而言，在5%、10%和20%的三种不对称条件下，所有自由度的RMSE和MAE值均保持在相对较低的水平，表明使用全对称模型辨识的等效参数能够有效重建不对称系统的动态响应。具体而言，平移响应的RMSE和MAE值大约保持在10??量级，而旋转响应的RMSE和MAE值大约保持在10??量级。此外，MAE值始终低于相应的RMSE值，表明预测误差主要分布在一个较小的范围内。下载：下载高清图片（732KB）下载：下载全尺寸图片图29. 解算精度和动态一致性评价指标示意图。相对误差和决定系数进一步证明了所提方法的有效性。Y和Z方向平移响应的相对误差分别约为1.34%–1.40%和4.44%–4.60%，而绕X、Y和Z轴旋转响应的误差分别为2.08%–2.18%、3.19%–3.54%和1.48%–1.49%。所有自由度的R2值均保持在0.85以上，最大值达到0.98，表明基于对称假设获得的等效参数仍能有效捕捉系统的主导动态特性。值得注意的是，响应误差并未随不对称度的增加而显著增大。这是因为所提辨识方法并不旨在独立恢复压缩侧和拉伸侧的参数，而是通过响应匹配辨识等效参数来表征系统的整体动态行为。因此，尽管随不对称度的增加，辨识得到的刚度和阻尼参数偏离实际的单侧参数，但重建响应仍与不对称参考模型保持高度一致。总之，不对称性分析表明所提辨识方法不依赖于响应生成模型与辨识模型之间的精确匹配。即使压缩侧刚度和阻尼参数与对称辨识模型中假设的参数相差高达20%，辨识得到的等效参数仍能提供满意的主导动态响应重建。这一结果表明，适度的模型差异主要将辨识参数的物理解释从局部接触参数转变为等效系统参数，而非使基于响应的辨识框架失效。除了支撑不对称性外，所提模型对滑块和直线导轨采用刚体假设，忽略了滑块、导轨和安装结构的弹性变形。虽然在实际工作条件下结构柔度可能对系统总变形有所贡献，但由于辨识过程基于实测系统响应，其影响部分地被纳入了辨识的等效参数中。因此，辨识参数代表直线导轨系统在特定安装、预载荷和载荷条件下的整体等效动态特性，而不仅仅是滚珠接触界面的局部参数。在本研究的实验范围内，在不同激励条件下系统主导频率几乎不变，表明系统可近似视为线性时不变系统。结果表明，刚体五自由度模型在小振幅振动、常规预载荷和有限载荷变化条件下，在计算效率、参数可辨识性和预测精度之间取得了合理平衡。对于大载荷条件、显著预载荷变化或结构柔度起主导作用的情况，未来研究应纳入结构柔度和载荷相关的不对称接触模型，以进一步提高预测能力。7. 结论本文提出了一种基于质心位移的直线导轨刚度和阻尼辨识方法。通过整合理论建模、参数辨识、仿真验证和实验验证，该系统地研究了五自由度激励下滚动直线导轨的动态特性和参数辨识问题。主要结论如下：（1）提出了一种基于频域积分和最小二乘算法的质心位移重建方法。该方法通过测量多点加速度响应，利用频域积分和多节点三次样条插值获得位移时间历程。此外，通过结合坐标变换和最小二乘拟合，实现了滑块质心位移的高精度重建。（2）建立了基于对称性和均匀性假设的线性五自由度动力学模型。该模型将滑块与导轨之间的接触简化为零长度弹簧，并进一步等效为两排连续分布单位刚度和阻尼的支撑线。基于此，利用拉格朗日方程推导了滚动直线导轨系统的运动方程。（3）提出了一种线性刚度和阻尼参数辨识方法。该方法采用Newmark-β方法进行时间域离散化，并结合最小二乘优化算法，同时辨识横向和竖直方向的单位长度等效刚度和阻尼参数。全面的实验验证、比较分析和参数敏感性研究表明，所提方法能够准确预测五自由度激励下直线导轨系统的动态响应，同时表现出良好的参数稳定性和模型适用性。此外，参数敏感性分析和不对称性分析验证了系统动态响应主要由刚度参数控制，并确认了基于全对称假设建立的五自由度模型在小振幅振动、常规预载荷和有限载荷变化条件下的合理性。总之，所提方法能够有效辨识直线导轨系统的等效刚度和阻尼参数，并实现高精度的动态响应预测，为直线导轨动态分析、振动抑制策略设计和机床进给系统动态性能优化提供了理论基础和技术支持。CRediT作者贡献声明李卓阳：撰写—初稿，验证，方法学，数据管理，概念化。曹书杰：可视化，验证，软件，形式分析。

Jinsong Zhao：可视化、软件、形式化分析、数据整理。Changyou Li：可视化、验证、资源、数据整理。Chunyu Zhao：写作 – 审阅与编辑、监督、资源、项目管理、资金获取。Guangyong Song：软件、调查、形式化分析。Hang Lu：软件、资源、数据整理。Zhuoyang Li|Shujie Cao|Jinsong Zhao|Changyou Li|Chunyu Zhao|Guangyong Song|Hang Lu