混合驱动框架下主轴-轴承系统热特性分析：热网络、粒子滤波与SINDy协同

简介：面向机械工程、热力学与数据科学交叉领域的技术资源包，聚焦主轴-轴承系统在多变工况下的热特性建模。基于论文复现思路，系统介绍数据驱动与模型驱动相结合的混合框架，解决物理机制复杂、温度传感器部署受限时的热源计算与参数估计难题；内容涵盖热网络法温度预测模型、SIAN与Sobol全局敏感性分析、状态空间方程构建、结构可识别性与可估计性评估、多层粒子滤波（MLPF）及稀疏识别，并附详细可运行Python代码与逐段解释。包体为单个PDF文档，共929KB，既有论文解读，也有完整可运行代码与注释。适合有一定编程基础、从事主轴-轴承系统设计优化的研发人员拓展实践，目前已有131人浏览学习。读者可据此掌握少传感器条件下温度场的准确预测与热参数辨识思路，并结合具体应用场景扩展算法，提高预测鲁棒性并减少对大量温度传感器的依赖。

1. 主轴-轴承系统热特性分析：混合驱动框架凭什么用少量传感器稳住温度场

机床主轴高速运转时，轴承温度每升高10°C，润滑脂寿命就缩短一半，热伸长直接改变加工精度——做过主轴热特性分析的人都清楚这个痛点。传统热网络模型依赖专家经验定热阻，纯数据驱动模型又缺乏物理可解释性，两边都有硬伤。这份资源复现的是一篇基于混合驱动框架的主轴-轴承系统热特性分析论文，把热网络法、Sobol全局敏感性分析、多层粒子滤波(MLPF)和SINDy稀疏识别串成完整链路，用少量温度传感器就能对不同工况下的温度场做有效预测。适合有机床设计、热误差补偿、轴承故障诊断需求的工程师，也适合正在复现论文的研究生。代码用Python写成，热源计算、热阻网络、粒子滤波、稀疏识别都有可运行实现和逐段解释。

2. 混合驱动框架的底层逻辑：热网络模型与参数可识别性怎么衔接

2.1 热网络法为什么够用又不够用

热网络法的思路很直白：把主轴-轴承系统离散成若干节点，轴承、轴、箱体、环境各占一个，节点之间用热阻连接，每个节点挂一个热容，然后按热平衡方程迭代求解。好处是物理意义明确、计算量小，四节点模型跑几百步也就是毫秒级，可以直接嵌进数控系统的热误差补偿逻辑里实时运行，这是它到今天仍然是工程主流的原因。

问题出在参数上。热阻不是查手册能查准的，接触热阻跟表面粗糙度、装配预紧力、润滑油膜状态都有关，而这些量在机床运行中是动态变化的。论文里专门点名了几个常被忽略的因素——热变形耦合效应、制造和装配误差、接触压力变化，这些都会让固定参数的热网络模型在不同工况之间迁移时误差放大。说穿了就是一句话：模型结构是对的，但参数不可信。

混合驱动框架的思路就是承认这一点：物理模型提供结构骨架，数据驱动负责修正不确定性。先通过可识别性分析筛出哪些参数值得估计，再用粒子滤波把这些参数和节点温度一起在线估计，最后用稀疏识别从数据里校验模型结构是不是对。三步走下来，既保住了物理可解释性，又获得了纯数据驱动方法的适应能力，这也是标题里"混合驱动"四个字的核心含义。

2.2 SIAN与Sobol全局敏感性分析：从一堆热参数里筛出关键变量

框架的第一步不是上来调参，而是先做参数可识别性分析，这里用了两个工具配合。SIAN解决的是"这个参数在数学上能不能被观测数据唯一确定"的问题——比如两个接触热阻串联在同一个测点路径上，它们对温度的影响方式几乎一样，单独估计任何一个都是病态的，强行辨识只会得到一堆不稳定的伪解。Sobol全局敏感性分析解决的是"哪个参数对温度输出影响最大"的问题，它不假设参数和输出之间是线性关系，能在整个参数空间上评估每个参数的独立贡献和交互贡献。

代码里用SALib实现Sobol分析的套路是这样：

from SALib.sample import saltelli from SALib.analyze import sobol # 定义参数空间：两个参数，采样区间都在[0,1] problem = { 'num_vars': 2, 'names': ['param1', 'param2'], 'bounds': [[0, 1], [0, 1]] } # 生成Saltelli采样序列 param_values = saltelli.sample(problem, 1000) outputs = np.array([example_model(*params) for params in param_values]) # 计算Sobol一阶指数S1和总效应指数ST Si = sobol.analyze(problem, outputs) print(f"一阶指数: {Si['S1']}") print(f"总效应指数: {Si['ST']}")

注意saltelli.sample(problem, 1000)里的1000不是总样本数，而是每个参数的采样基数，实际生成的样本数量是N * (2D + 2)，D是参数个数。参数一多，总样本数会暴涨，所以我一般先设1024跑一次，看ST指数是否稳定再决定要不要降基数。筛选参数时看ST比看S1更稳妥：ST包含参数自身贡献和它与其它参数的交互作用，ST明显偏小的参数固定成经验值就行，不必放进粒子滤波的状态向量；ST大的参数才是需要在线估计的对象。

代码里的SIAN部分做了简化处理，直接打印"假设所有参数在数学上都是可识别的"。实际用的时候要小心这个假设，SIAN需要基于雅可比矩阵的符号计算，工程上我常用数值扰动法近似：对每个参数加微小扰动，看输出温度轨迹是否产生可区分的变化，变化不可区分的参数直接合并或固定，省得给粒子滤波留一堆不可观测的维度。

2.3 状态空间方程：把热网络改造成适合粒子滤波的形式

筛选完参数，下一步是把热网络模型改写成状态空间方程，这是连接物理模型和粒子滤波的接口。热网络的离散化迭代本身就是天然的递推形式：x(k+1) = f(x(k), u(k), θ) + w(k)，其中x是各节点温度向量，u是转速、载荷这些工况输入，θ是筛选出来的不确定参数，w是过程噪声。观测方程y(k) = h(x(k)) + v(k)把节点温度映射到传感器实际安装位置的读数。

代码里ThermalNetworkModel的update_temperatures就是f的具体实现，每调用一次完成一步状态递推，输出下一时刻的节点温度向量。粒子滤波做的事情，就是每来一个温度观测值，用这个递推模型做一步预测，再按预测和实测的偏差调整粒子权重，最终得到节点温度和不确定参数的联合后验估计。这也是为什么写状态空间方程时，θ必须显式地放进状态向量里，而不是藏在模型的类属性中——滤波算法只认状态向量，不认外部变量，这个顺序搞反了，后面所有代码都会绕远路。

3. 热源计算与热阻网络建模：轴承发热量算不准后面全白做

3.1 Palmgren摩擦热公式与热变形耦合修正

热网络模型的输入是热源，主轴-轴承系统里最主要的热源就是轴承摩擦热。代码用的是Palmgren公式的工程简化版，把摩擦扭矩分成黏滞摩擦和载荷摩擦两项：

class BearingHeatSource: def __init__(self, bearing_params): """ bearing_params: 必须包含 viscosity: 润滑油运动黏度(mm²/s) dm: 轴承节圆直径(mm) f0: 载荷相关摩擦系数 """ self.params = bearing_params def calculate_heat(self, speed, load, temp): """speed: 转速(rpm), load: 载荷(N), temp: 轴承温度(°C)""" friction_torque = self._calculate_friction_torque(speed, load) # 摩擦扭矩与转速折算成发热功率(W) heat_basic = 1.047e-4 * friction_torque * speed # 热变形修正：温度升高会改变游隙和接触状态 thermal_deformation_factor = self._thermal_deformation_effect(temp) return heat_basic * thermal_deformation_factor def _calculate_friction_torque(self, speed, load): # 黏滞摩擦项：与转速的1.5次方成正比，高速时占主导 viscous = (10e-7) * (self.params['viscosity']**0.5) * (speed**1.5) * (self.params['dm']**3) # 载荷摩擦项：与载荷成正比，低速重载时占主导 load_component = self.params['f0'] * load * self.params['dm'] return viscous + load_component def _thermal_deformation_effect(self, temp): # 相对室温(25°C)的温差，每升高1°C热源增加0.5% delta_temp = temp - 25 return 1 + 0.005 * delta_temp

1.047e-4 * friction_torque * speed是把摩擦扭矩和转速折算成发热功率的常用系数，这个常数值得背下来，做轴承热分析到处用。黏滞项里speed**1.5说明高速工况下发热是超线性增长的，电主轴高速空转时的温升就是这么来的；标准Palmgren公式里黏滞项其实是(ν·n)^(2/3)的形式，代码做了简化，趋势一致但标定时要注意。_thermal_deformation_effect对应论文强调的热变形耦合：温度升高→轴承游隙变化→摩擦力矩变化→发热量变化，代码用每升高1°C热源增加0.5%的线性修正近似，工程上够用，但这个系数在不同预紧力下要重新标定，不能一个值走天下。

3.2 热阻网络模型的Python实现与节点设计

有了热源，核心就是热网络主体，代码实现是这样的：

class ThermalNetworkModel: def __init__(self, nodes, connections, initial_temp=25.0): """ nodes: 节点名列表，如['bearing', 'shaft', 'housing', 'ambient'] connections: 连接关系列表，元素为(node1, node2, thermal_resistance) """ self.nodes = nodes self.connections = connections self.temperatures = {node: initial_temp for node in nodes} self.heat_sources = {node: 0.0 for node in nodes} def update_temperatures(self, dt=1.0): """按热平衡方程递推一步，dt为时间步长(秒)""" new_temps = {} for node in self.nodes: total_heat = self.heat_sources[node] # 累加所有相邻节点通过热阻传来的热流 for conn in self.connections: if node in conn[:2]: other_node = conn[0] if conn[1] == node else conn[1] resistance = conn[2] temp_diff = self.temperatures[other_node] - self.temperatures[node] total_heat += temp_diff / resistance heat_capacity = 1000.0 # 温升 = 总热流 * 时间步长 / 热容 temp_change = total_heat * dt / heat_capacity new_temps[node] = self.temperatures[node] + temp_change self.temperatures.update(new_temps) return self.temperatures

update_temperatures的核心逻辑就是一句话：每个节点的温升等于"自身热源加上所有相邻节点通过热阻传过来的热流"乘以时间步长再除以热容。temp_diff / resistance是傅里叶导热定律的离散形式，热流方向由温差符号决定，温度高的节点向温度低的节点传热。连接关系里('bearing', 'shaft', 0.5)表示轴承和轴之间热阻0.5 K/W，热阻定义是温差除以热流，数值越小导热越强。

拿到代码后我会立刻做两处改造。第一，热容1000.0写死会导致所有节点温升速率一样，实际轴承座热容和轴段热容差一个量级，不改的话温升时序全错。第二，环境节点必须固定温度不参与更新，否则整个系统温度会整体漂移，这一步我在第5章会展开讲。

3.3 热容、热阻初值怎么定：量级比精度重要

节点和热阻的初值设计，我的经验是"量级对，精度可以不那么讲究"——因为后面粒子滤波会在线修正这些参数，初值只要保证滤波能从正确的搜索空间起步就行。四个节点是起步配置，要做前后轴承温差分析就加到六个节点，但节点越多参数辨识的欠定问题越严重，别贪多。

热容按节点代表区域的质量乘比热容估算，钢的比热容约460 J/(kg·°C)，轴承节点质量小所以热容小、温升快，箱体反过来。热阻量级的区分很明显：固定接触面的接触热阻通常在0.001~0.01 K/W，对流换热热阻在0.1~1 K/W，代码里0.5、0.8、1.2这几个值基本落在这个分布里。初值不用精确，但要确认方向——热源只加在轴承节点，轴段和箱体之间不能有直连热阻，否则热量路径就错了。

4. 多层粒子滤波与稀疏识别：数据驱动部分怎么落地

4.1 粒子滤波为什么适合热参数估计

热参数估计的难点在于系统强非线性，热阻到温度之间的映射不是线性的。标准卡尔曼滤波要求线性高斯系统，直接用不了；EKF靠一阶泰勒展开近似，在热网络这种强非线性下误差累积快，而且热参数的后验分布经常是非高斯的，EKF强行用高斯分布拟合会丢信息。粒子滤波用一组带权重的粒子近似后验分布，不假设分布形态，天然适配这个场景。

论文的MLPF在标准粒子滤波上加了分层结构：一层做状态估计（节点温度），一层做参数估计（热阻、热容），中间用SINDy稀疏识别的结果做层间信息交换。分层的理由是状态和参数的收敛速率不一样，状态几秒内就能稳定，参数需要长期观测才能收敛，混在一个滤波器里容易出现状态收敛了参数还在飘的"跷跷板"现象。分层之后各层可以用不同的粒子数和噪声设置，整体鲁棒性明显更好。

4.2 MLPF代码实现：预测、更新、重采样与关键参数

class MultiLayerParticleFilter: def __init__(self, n_particles, state_dim, obs_dim, process_noise, obs_noise): self.n_particles = n_particles self.state_dim = state_dim self.obs_dim = obs_dim self.process_noise = process_noise self.obs_noise = obs_noise # 初始化：从标准正态分布采样n_particles个粒子，权重均匀 self.particles = np.random.randn(n_particles, state_dim) self.weights = np.ones(n_particles) / n_particles def predict(self, transition_func, control_input=None): """预测步骤：用状态转移函数推进每个粒子，添加过程噪声""" for i in range(self.n_particles): if control_input is not None: self.particles[i] = transition_func(self.particles[i], control_input) else: self.particles[i] = transition_func(self.particles[i]) self.particles[i] += np.random.multivariate_normal( np.zeros(self.state_dim), self.process_noise) def update(self, observation, observation_func): """更新步骤：高斯似然更新权重，归一化后重采样""" for i in range(self.n_particles): pred_obs = observation_func(self.particles[i]) error = observation - pred_obs # 马氏距离形式的高斯似然 self.weights[i] *= np.exp(-0.5 * error.T @ np.linalg.inv(self.obs_noise) @ error) self.weights += 1e-300 # 防下溢保底 self.weights /= np.sum(self.weights) self.resample() def resample(self): """多项式重采样：按权重抽取粒子，抽中的粒子被复制""" indices = np.random.choice( range(self.n_particles), size=self.n_particles, p=self.weights) self.particles = self.particles[indices] self.weights = np.ones(self.n_particles) / self.n_particles def estimate(self): """加权平均得到状态估计值""" return np.average(self.particles, weights=self.weights, axis=0)

predict里的transition_func在热网络场景下就是ThermalNetworkModel.update_temperatures的封装，过程噪声协方差process_noise代表模型误差——热网络本身是简化模型，这个噪声不能设成零，否则滤波会过度信任模型，观测数据修正不了偏差。update里error.T @ np.linalg.inv(self.obs_noise) @ error是马氏距离，相当于把观测误差按噪声协方差归一化。观测噪声协方差obs_noise按传感器真实精度设，我用±0.1°C的PT100时一般取方差0.02左右。

重采样是粒子滤波的命门。1e-300这行不能删，连续几轮更新后权重乘积会下溢成零，没这个保底直接出NaN。粒子数1000对四节点热网络够用，状态维度超过6就往3000加，但计算量跟着翻倍，实时性要求高的场合要权衡。资源里还带了一个EnhancedMLPF扩展类，在标准滤波基础上加了sparse_identification_layer，把稀疏识别结果接入滤波循环做层间交互，这正是论文里"多层"的落点。

4.3 SINDy稀疏识别与MLPF的分层协作

SINDy（稀疏识别非线性动力学）的作用是从温度时间序列里自动发现状态变量之间的动力学关系，不靠人工推导，是数据驱动部分最核心的一环。

import pysindy as ps def sparse_identification(X, t, derivatives=None): """X: 状态时间序列，t: 对应时间点""" if derivatives is None: derivatives = np.gradient(X, t, axis=0) # 数值求导 model = ps.SINDy( optimizer=ps.STLSQ(threshold=0.1), # 稀疏化阈值 feature_library=ps.PolynomialLibrary(degree=2) # 候选函数库 ) model.fit(X, t=t, x_dot=derivatives) model.print() return model

STLSQ(threshold=0.1)会把系数绝对值低于0.1的项剔除，得到最简方程。阈值太小方程项太多，失去"稀疏识别"的意义；阈值太大把真实动力学项也删了。我通常从0.1起步跑一遍看项数，再在0.01到0.5之间做网格搜索，用"拟合残差 + 非零项数"做交叉验证。PolynomialLibrary(degree=2)限定候选函数库为二次多项式，热网络里温度项之间基本是线性加弱非线性关系，二次足够，加高次项只会助长过拟合。

在MLPF框架里，稀疏识别不是独立跑的，而是和粒子滤波分层协作：粒子滤波收敛后，把估计出的温度轨迹喂给SINDy，如果识别出的方程结构和当前热网络不一致——比如多出一项非线性耦合——说明模型结构错了而不是参数错了。这时候该回头改热网络结构，而不是继续调滤波参数。这个判断帮我省过好几次白调参的冤枉时间，遇到温度预测误差压不下去时，先跑一次稀疏识别看结构，比闷头改协方差矩阵高效得多。

5. 避坑指南：热参数估计过程中的五个高频翻车现场

5.1 热网络阶段的坑：环境节点、热容初值与热阻连接

翻车现场一：环境节点温度跟着系统一起漂移。

现象：仿真跑几百步，环境节点从25°C飘到40°C，整个温度场整体抬升，怎么调热阻都压不回来。

原因：update_temperatures对所有节点一视同仁地做热平衡更新，没有把环境节点设置为固定温度边界，散热路径变成"假散热"——热量传出去又通过环境节点的温升传回来一部分。

解决：更新前把环境节点从迭代列表里剔除，或者每步更新后强制temperatures['ambient'] = initial_temp。我习惯在类里加一个fixed_nodes参数把这层逻辑固化进去，一劳永逸，避免每次建模都重复踩。

翻车现场二：所有节点共用同一个热容，温升曲线全长得一样。

现象：轴承节点和箱体节点温升速率几乎相同，跟实测曲线对不上，"轴承先热、箱体后热"的时序出不来。

原因：代码里heat_capacity = 1000.0是写死的，没有区分节点代表区域的质量差异，热惯性全被抹平了。

解决：改成节点字典，按每个节点代表区域的质量乘比热容算。轴承节点热容小、温升快，箱体热容大、温升慢，这个先后顺序是热网络的灵魂，省掉这一步等于白建模。

翻车现场三：热阻连接写错节点对，热量串流。

现象：热量从低温节点流向高温节点，明明只有轴承加了热源，轴段温度反而比轴承还高。

原因：连接列表写错，比如把('shaft', 'bearing', 0.8)这样的对写进去，让轴段和轴承直接互通，绕过了本该经过的热阻路径，热量按错误路径传导。

解决：建模后做一次邻接矩阵可视化，确认每个节点只连线物理上相邻的部件。热源只挂在轴承节点上，轴段和箱体之间不要有直连热阻，检查完再跑仿真。

5.2 粒子滤波与敏感性分析阶段的坑

翻车现场四：粒子权重退化，有效粒子数骤减。

现象：滤波跑到几十步，estimate()结果开始剧烈跳动，重采样后粒子几乎全是同一个值的复制品，状态估计完全失去意义。

原因：观测噪声设得太小，似然函数过于尖锐，少数粒子独占了几乎所有权重；或者过程噪声太大，粒子早早散掉失去了多样性。这个权衡本身就是粒子滤波里最玄学的地方。

解决：实时算有效粒子数1 / np.sum(self.weights**2)，降到粒子总数的三分之一就提前触发重采样。观测噪声方差不要小于传感器的物理精度，过程噪声从0.01 * np.eye(state_dim)起步，再根据残差逐步调整，一次调到位很难。

翻车现场五：Sobol分析两次结果对不上，参数筛选结论变来变去。

现象：同一组参数范围，两次运行ST指数差异巨大，这次筛掉的热阻下次又变成敏感参数，参数筛选结论完全不可信。

原因：saltelli.sample(problem, 1000)的样本数跟参数维度直接挂钩，参数一多基数1000覆盖不了高维参数空间，而且没有固定随机种子，两次采样完全不同。

解决：固定随机种子np.random.seed(42)，把基数提到2048以上，先跑一次确认ST稳定再取结果。参数维度过高时，先做一轮粗筛把明显不敏感的参数固定掉，再对剩余参数做精细分析，别指望一次到位。

6. 把预测误差压下来的验证技巧：从仿真数据到实测数据的过渡

前面所有代码跑通，只代表框架在仿真数据上闭环了，真正落地到主轴-轴承系统还得按顺序验证，这个顺序错了很容易被误差来源误导。第一步是残差分析：拿仿真输出和实测温度做差，画温度残差曲线。残差有趋势性，比如前100秒持续为正，说明模型结构有问题，不是噪声问题；残差是零均值白噪声，才继续往下走。第二步是工况外推测试：用低速工况数据训练，用高速工况数据验证。这是混合驱动框架最有价值的地方，纯数据驱动模型在这一步通常直接崩掉，而物理模型加数据修正的结构能扛住一部分外推。

第三步是关键参数的物理合理性检验。粒子滤波估计出的热阻值，要跟材料手册和经验公式对比数量级。我遇到过一次滤波收敛后轴承热阻估计出0.0001 K/W，比理论值低两个量级，但温度预测误差仍然很小——这说明参数过拟合了，模型在用错误的参数拟合正确的输出。从那以后我每次跑完滤波都强制走一遍参数合理性检查：把每个估计值和经验范围逐一对比，超出一个数量级就直接判模型失败，重新审查热网络结构，绝不看温度误差小了就放行。

还有一个容易漏的细节：传感器安装位置要跟热网络节点严格对应。热网络里的bearing节点代表轴承内圈发热点，但传感器只能装在轴承座上，中间隔着接触热阻和热容延迟。我一般会在观测函数observation_func里加一个固定延迟项，而不是让传感器读数直接等于节点温度，这个修正能让滤波收敛速度提升一个量级，值得优先试。数据驱动和模型驱动不是二选一，而是像论文里那样分层协作——模型提供结构骨架，数据负责修正不确定性，希望这份拆解能帮到你，少走几个月的弯路。

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