0.   引言

轴承在航空发动机等机械系统中起着至关重要的作用，确保其长寿命的使用。为提升轴承热处理的品质，有必要对热处理过程中的炉温进行监测与解析。传统炉温测量一般采用热电偶布点的接触式监测技术，通过定期校准可获得炉内环境大致状况，拆装流程复杂，而且测量流程破坏被测环境，通常仅依赖于有限的热电偶装置无法达到对整个炉内状况的掌握和分析；传统的非接触测温技术通常需要直接对准目标。然而，大多数轴承热处理炉的炉壁上都没有视窗，无法对内部测温区域进行精准瞄准。因此，有必要研究相关传热反演分析算法以协助实现炉内温度环境的分析研究。

关于轴承热处理炉传热反演分析，属于一种典型的红外信号分析反问题。针对问题的研究始于数学与物理的反问题，这是一门新兴的研究领域，其中热传导反问题是这类问题中的关键部分^[1]。一个完整的热传导正问题应当涵盖以下几个方面的已知信息中的一种或多种^[2]：①描述热传导现象的偏微分控制方程，指描述热传导过程的数学模型，用来描述物质内部温度分布随时间变化的规律的方程，基于能量守恒和热传导定律，能够反映热量在物体内部从高温区域向低温区域传递的过程；②用于描述实验材料热物理性质的参数，包括热导率、热容、热扩散率、热膨胀系数、辐射系数、发射率、辐射角系数等，这些参数反映了材料在热传导、热存储和热响应等方面的特性；③热传递过程的起始状态以及边界约束条件，前者指热传递过程开始时物体或系统的温度分布状态，后者指热传递过程中物体边界上的温度或热流密度的分布情况，他们是热传导问题中的两个核心要素，共同决定了热量如何在物体内部传递以及温度如何随时间变化；④物体的几何形态及其内部热源的分布情况，这两个信息共同决定了热量在物体内部的传递方式和温度分布特性，在设计和优化热传导系统时，需要充分考虑这两个因素的影响。

在热传导反问题中，我们的目标是基于已获得的物体温度数据，去推算在上述已知条件中尚未明确的任何参数或状态。然而，热传导反问题领域面临的一项主要挑战，就在于如何有效地削弱噪声对反问题求解准确性的干扰。

目前，针对热传导反问题的解决方法主要可以分为确定性算法和随机性算法。在确定性算法领域，其中，共轭梯度法因其高精度和适用性而受到广泛应用。例如，清华大学的研究者们通过共轭梯度法来反演计算管道内壁温度^[3]，以实现对计算管道内壁温度的模拟。与此同时，北京空间飞行器总体设计部的研究人员则采用共轭梯度法反演航天器在轨运行时的瞬态外部热流^[4]。随机性算法中，常用的方法包括但不限于：QPSO优化算法和遗传算法。例如，江南大学的研究者采用QPSO优化算法来估算热传导系数^[5]；而大连理工大学的研究人员则运用遗传算法来反演计算热传导边界条件^[6]。

共轭梯度法在迭代计算过程中，对初始值的选取标准较高。若初始值的选取不当，解的稳定性将难以保证。随机性算法的初始值是任意生成的，这意味着算法在选取初始值时无严格规定。然而，由于随机性算法在计算过程中不采用正则化策略，因此其对噪声的抵抗力相对较弱。

在计算机视觉领域，一种广泛应用的基础原理是小波变换的滤噪技术——小波滤噪。小波变换将输入信号映射到小波域，根据噪声的频率系数在不同尺度上呈现出不同的属性和机理，对含有噪声的小波系数进行处理。小波去噪在一定程度上可视为低通滤波器，去噪后仍能保留信号特性在这一方面胜过传统的低通滤波器。

本研究采用量子粒子群随机优化算法（quantum particle swarm stochastic optimization algorithm, QPSO）以解决热传导反问题，同时为了提高算法的抗噪性，引入了小波滤噪机制。这一机制有助于提高算法对噪声的抵抗力。

1.   热传导反问题数学模型

传热过程示意图如图 1所示，物体为“无限大”物体，其厚度为$\delta $，物体的长、宽为无限大，在物体的一侧表面有温度为U_∞的气体进行对流加热，物体的另一侧表面满足理想绝热条件。

图  1  传热过程示意图

Figure  1.  Schematic diagram of heat transfer process

下载: 全尺寸图片 幻灯片

传热过程可采用下面的数学表达式来描述：

式中：K为热导系数，也就是表征材料导热能力大小的物理量；h(τ)为物体对流面的传送系数，通常指的是对流传导系数（也称为对流传热系数或表面传热系数），它反映了热传导过程中物体表面传热的效率和速度；a为被测物体的热交换系数，亦称为传热系数或热导率，是衡量单位时间内单位面积的物体表面通过单位温差传递热量能力的物理量；U_∞为物体初始温度，是物体在开始某个过程、实验或分析时所处的温度状态；U₀为进行对流的气体温度，其是一个动态变化的参数，它受到多种因素的影响，包括气体的性质、对流条件、外部环境等；U(x, τ)为物体内部温度场。

定义反演算法的控制函数如下：

之后我们运用QPSO来寻找一个恰当的h值（即表面传热系数）。这个h值的选取标准是使得通过计算得到的误差小于我们预先设定的反演误差。一旦找到了满足这一条件的h值，我们就可以认为已经成功地确定了材料表面的传热系数。

2.   QPSO随机优化算法

2.1   算法介绍

量子粒子群随机优化算法，缩写为QPSO，是一种在PSO优化算法的基础上进行提升的算法^[7-8]。在传统的PSO优化方法中，粒子的位置取决于其位置及移动速度共同影响，并且根据Marcelo J. Colaco和Helcio R. B.的研究结果，当粒子不断收敛，粒子的收敛区域逐渐缩小，收敛能力不足^[9]。

为了克服传统粒子群优化算法（PSO）存在的局限性，QPSO引入了一种新的概念，即将算法中的粒子视为具有量子特性的实体。这些量子粒子在量子空间中，根据一定的概率分布，在由某种势场决定的位置出现。这种量子行为使得粒子的搜索不再局限于局部区域，而是能够在全局范围内进行探索。因此，QPSO算法通过利用量子粒子的这一特性，实现了对解空间更广泛、更深入的搜索，从而提高了找到全局最优解的可能性，这种全局搜索能力有助于算法在复杂的优化问题中找到更好的解。

在QPSO算法中，粒子在量子空间中的位置由势场决定，这种势场引导机制使得粒子能够更高效地探索解空间，从而提高算法的收敛速度和精度。与PSO算法相比，QPSO算法需要设置的参数更少，在PSO算法中，通常需要设置惯性权重、加速常数等参数，而在QPSO算法中，主要需要设置的参数是创新参数a，这种参数设置的简化使得算法在实际应用中更加易于实现和调试，并且QPSO算法对参数的敏感性较低，这就代表即使参数设置不是最优的，算法仍然能够取得较好的性能。

图 2展示了QPSO优化算法的计算流程。

图  2  QPSO优化算法实现流程图

Figure  2.  QPSO optimization algorithm flowchart

下载: 全尺寸图片 幻灯片

在该计算流程中值得注意的是确定粒子群的规模至关重要。在优化问题中，每个粒子代表一个候选解，经过计算最终得到位于群体最优位置的粒子成为最优解。

开展QPSO实验前需要先定义几个量的含义：

① 粒子数N；

② 某个粒子的顺位序号排列为i；

③ 第i个粒子的坐标位置X_i；

④ 历史最佳坐标值P_i；

⑤ 全局最佳坐标值P_g。

粒子的势中心点位置θ_i为：

在式(3)中，任意随机数的分布应与均匀分布的随机数接近，取值范围通常设定在[20,30]区间内较为适宜。${r_1}$和${r_2}$在初始阶段，每个粒子的初始值被随机生成，接下来会对这些粒子进行精确的适应度计算。对于当前粒子位置的评估，其数值越小表示粒子与最优位置之间的距离越近，这种评估方法可以应用在式(3)中。接着根据各个粒子位置的适应程度来寻找粒子的历史最佳位置，其计算公式如下：

该表达式表示了粒子适应度函数F(x)，它反映了当前第i个粒子的历史最优位置P_i(t)，以及第i个粒子在更新后的历史最优位置P_i(t+1)。群体的全局最佳点是根据每一个案例的总体最优点确定的，其确定方法如下：

当已知历史最佳位置和总体最佳位置时，根据下式来调整粒子的分布：

式中：u为在(0, 1)上均匀分布的随机数。

L_i(t)的定义如下：

式(7)的系数$\alpha $，其数值一般根据下式由1降低到0.5：

式(7)中C(t)为平均最优位置，其计算公式如下：

将步骤重复进行，当达到迭代次数时，系统会终止迭代。在此时，群体中的全局最优位置即为优化问题的最优解。

2.2   算法的实现

采取随机的方式对粒子进行初始化，其中初始值的范围在[h_min, h_max]内：

接着，对各个粒子X_i(t)代入方程式(2)中的各个位置。因此，如果粒子适应度越小，则表明其与测量值之间的接近程度越高，反之亦然。通过计算每个粒子的历史最优值P_i和群体中粒子的全局最优值P_g，根据式(7)来更新粒子的最新位置。在继续进行迭代操作时，若达到预定的迭代次数或粒子已收敛，那么迭代过程将终止并退出。

QPSO反演算法的完整解方程表达式如下：

① 设置N个粒子的初始状态和粒子的历史最佳状态值X_i(0)＝P_i(0), (i＝1, …, N)；

② 初始化迭代次数t＝1和收缩系数α＝1；

③ 对各个例子进行适应度评估J(X_i(t))。若J(X_i(t))＜P_i(t)，更新粒子的历史最优值为P_i(t+1)＝X_i(t)；若J(P_i(t))＜P_g(t)，更新群体的全局最优值P_g(t+1)＝P_i(t)；

④ 根据式(9)计算群体平均最优值C(t)；

⑤ 根据式(4)和式(7)对每个粒子的数值进行更新；

⑥ 根据式(8)进行收缩系数α的优化迭代计算，为下一次预迭代做准备；

⑦ 若t＜t_max，则代表公式流程需要进行转换，需要把流程步骤跳转到③；否则，此刻应该退出迭代计算流程。

相较于其他PSO优化方法，QPSO在全球搜索能力方面表现出优于其他算法的优势，特别适用于解决空间广阔且复杂的问题^[10-11]。通常，QPSO只需进行较少的参数调整，这样其简洁性和易于使用的特点使得其在解决局部最优问题方面具有竞争优势。另一方面，由于量子行为特性，QPSO在避免局部最优方面具备明显优势。

3.   小波阈值滤噪算法

采用小波阈值滤波算法对测量值进行降噪处理，以降低噪声干扰，并提升反演算法的精确性。

检测到的电信号如下式所示：

式中，s(t)为方差为σ²的高斯白噪声n(t)，噪声则遵循正态分布N(0, σ²)。

为了从物理信号f(t)中提取原始信号s(t)，难度较大，需要通过其他变换方法来辅助。

在实际测定时，测得数据f(t)可变换成w_{j, k}：

通过以下论述将原始信号与小波系数之间的关系表述清晰。在对白噪声进行分析时，我们发现其对应的尺度下，小波系数呈现出均衡分布。在某些范围内小波系数的值会逐渐减小。依据信号特性，在这个公式中选取适当的阈值，将小于阈值的波形系数置零，或者保留或缩小大于阈值的波形系数。通过这种措施，可以实现原始信号与噪声的分离。

对于阈值的确定方法这一问题，Donoho提出了如下的公式来确定阈值^[12]：

式中，N为对应的小波系数长度；

$\sigma $为噪声的标准方差。

对于$\sigma $可通过经验公式(14)计算得出：

在这个过程中，我们需要计算小波系数的个数，相应尺度下小波系数的个数记为$N$。

对数据的滤波方法通常可以使用硬阈值法与软阈值法两种，而硬阈值法含义可如下：

软阈值法的公式组定义形式如下：

两个滤波方式的图如图 3所见^[13]，由图中的函数曲线可以知道：硬阈值法由于在阈值$ \pm \lambda $附近会产生不连续点，所以处理后根据得到的${\widehat w_{j, k}}$来重构信号时有可能产生振荡；在采用软阈值滤波的过程中，${\widehat w_{j, k}}$连续性得到了显著改善。虽然，运用小波变换来分离信号中的噪声分量是一种有效手段，不过同一时间应该注意的是在某些场合下这种方法也存在一定的局限性。

图  3  软、硬阈值法滤波函数

Figure  3.  Soft and hard thresholding filtering functions

下载: 全尺寸图片 幻灯片

利用小波数学代换进行噪声过滤^[13]，示意图如图 4所示。

图  4  小波滤噪计算过程示意图

Figure  4.  Schematic diagram of wavelet denoising calculation process

下载: 全尺寸图片 幻灯片

小波阈值滤噪基本计算过程如下：

① 小波变换：选择适合的小波基对含噪信号进行变换，得到小波系数；

② 阈值处理：选择合适的阈值方法确定阈值。应用硬阈值或软阈值函数对小波系数进行量化处理，去除噪声成分；

③ 小波重构：利用处理后的小波系数进行逆变换重构得到去噪后的信号，通过计算信噪比（SNR）或均方根误差（RMSE）等指标，评估去噪效果的好坏。

这个过程通过分解信号、处理噪声成分和重构信号，实现了对含噪信号的有效去噪。

4.   算法仿真

本章利用仿真实验评估QPSO在求解逆热传导问题中的有效性。在噪声条件下，我们进一步采用了两种不同的处理方法：一种是直接应用QPSO算法进行计算；第二种方法是首先使用小波技术去除噪声，然后使用QPSO算法进行计算。在实验中，我们选择40Cr13轴承钢作为测试样本，通过逆有限差分法获得热传导正问题的温度场数据，评估算法的性能。算法平台为MATLAB2023b，实验材料采用40Cr13轴承钢。

4.1   直接仿真

直接仿真实验中没有噪声的加入。其主要目的是评估QPSO算法在理想条件下的性能，即当数据纯净、不受外界干扰时，算法能否准确、高效地解决热传导反问题。这样的实验设计有助于了解算法的基本特性和能力，通过对比算法输出与预期结果（通常是基于已知条件通过正问题计算得出的结果），可以验证算法在无噪声环境下的计算准确性。

在无噪声条件下，算法的运行时间、迭代次数等性能指标能够更直接地反映其计算效率，这有助于确定算法在最优条件下的收敛速度和稳定性。

直接实验的结果可以作为后续有噪声实验的基准，用于比较和分析噪声对算法性能的具体影响。这有助于深入理解噪声对算法精度和稳定性的干扰机制，并为算法的改进和优化提供方向，这里分不同情况进行了3组实验。

在第一组实验过程里K固定为80 W/(m²⋅K)，Δw＝0.004 m，Δt＝0.2 s，结果如图 5所示。

图  5  K=80 W/(m²⋅K)反演结果

Figure  5.  Inversion result when K=80 W/(m²⋅K)

下载: 全尺寸图片 幻灯片

图 5中的两条线分别代表真实值和反演值，其中红色的线为真实值，因为热导真值一直为K＝80 W/(m²⋅K)没变所以其为一条直线，而蓝色的线就是由反演计算出来的反演值，可以看出在直接实验中其虽然由于代数方程组在求解过程中产生的迭代误差而有所波动但是紧紧依附于真值的红线之上，可以看出两者很接近。

在这个实验中的迭代误差是无法消除的。迭代法的基本思想是从一个初始估计出发，通过不断迭代更新解向量，直至达到预定的精度要求，在这个过程中，初始估计的不准确性、迭代公式的近似性、计算机浮点运算的误差等因素都会产生迭代误差，为了减小迭代误差，需要采取合适的策略来优化迭代过程和提高计算精度，相关方法的比较与分析在本文将不再赘述。

图 6进行了相应的误差统计，其中最大绝对误差为7.69×10^－12W/(m²⋅K)。

图  6  K=80 W/(m²⋅K)误差统计

Figure  6.  Error statistics of K=80 W/(m²⋅K)

下载: 全尺寸图片 幻灯片

在第二组实验过程中K有100 W/(m²⋅K)和150 W/(m²⋅K)两种值，并且以方波的形式出现，如图 7所示。

图  7  第二组反演结果

Figure  7.  Inversion results of group 2

下载: 全尺寸图片 幻灯片

结果显示与第一组情况显示的结论类似。图 8显示了反演值的绝对误差分布情况，其中最大绝对误差为1.39×10^－11W/(m²⋅K)。

图  8  第二组误差统计

Figure  8.  Error statistics of group 2

下载: 全尺寸图片 幻灯片

方波是一种非正弦曲线的波形，具有“高”和“低”两个明确的值，理想方波的占空比为50%，即高电平在一个波形周期内占有的时间比值等于低电平所占的时间比值。激励为方波的实验在通信系统、信号处理、控制系统、物理实验、生物医学以及电路分析与测试等多个领域中具有重要的意义。通过合理的实验设计和注意事项的遵循，可以充分发挥方波激励信号在各个领域中的优势和作用。

在第三组实验过程采用斜波，反演结果如图 9所示。

图  9  第三组反演结果

Figure  9.  Inversion results of group 3

下载: 全尺寸图片 幻灯片

观察图 9，同样可以发现结果显示与第一组情况显示的结论相似。图 10显示了反演值的绝对误差分布情况，其中最大绝对误差为5.26×10^－11 W/(m²⋅K)。

图  10  第三组误差统计

Figure  10.  Error statistics of group 3

下载: 全尺寸图片 幻灯片

根据上述实验，整体的反演效果较好，反演值与真值非常接近。

4.2   经随机噪声处理的仿真

加入随机噪声的仿真实验的主要目的是评估系统或算法在噪声环境下的性能，通过模拟真实世界中的随机干扰，可以测试系统或算法的鲁棒性、稳定性和准确性。本节将对4.1节的数据进行±0.4%、±0.8%、±1.2%、±1.6%和±2.0%的随机噪声处理并开展相应的仿真实验，从而得到该实验的大致抗噪性能。

在不采用滤噪算法的情况下，直接进行经±0.4%随机噪声数据处理的反演计算。反演结果如图 11所示。

图  11  无滤噪处理反演结果

Figure  11.  Non filtered noise processing inversion results

下载: 全尺寸图片 幻灯片

根据图 11所示的数据可知无滤噪处理时反演效果较差，红色基线数值为80，而最大相对误差已经超过了500%。由此可见，反演算法在实际应用中已经失去了原本应有的功能。误差统计分布情况如图 12所示。

图  12  无滤噪的反演结果误差统计

Figure  12.  Error statistics of unfiltered inversion results

下载: 全尺寸图片 幻灯片

从图 12可以看出，误差统计结果较差，完全无法作为有效处理的信号，由此可见随机优化算法的QPSO虽然在选取初始值方面表现出色，在全局寻优能力、维数限制、参数设置敏感性和算法稳定性等方面存在一定的局限性，这主要是因为QPSO算法在更新粒子位置时，依赖于当前种群的最佳位置和个体最佳位置，从而可能忽略了其他潜在的更优解。

为了提升QPSO算法对噪声的抵抗力，我们引入了基于小波分析的滤波技术，结果如图 13所示。

图  13  小波分析反演结果

Figure  13.  Inversion results after wavelet analysis

下载: 全尺寸图片 幻灯片

根据图 13中的数据，我们可以推断，当采用滤波器噪声修正后的温度数据进行反演时，反演结果的准确度得到了显著的提升。

图 14展示了经过小波去噪处理分析后的信号误差统计图。

图  14  小波分析后的误差统计

Figure  14.  Error statistics after wavelet analysis

下载: 全尺寸图片 幻灯片

从图像中可以看出，经过滤波后的信号与原始信号之间的差异十分小，相对误差小于3.2%，这进一步佐证了小波去噪算法的效果良好。

进行相同的实验，其他实验同样采用小波去噪去处理数据，用来处理的随机噪声数值分别取为±0.8%、±1.2%、±1.6%和±2.0%，他们的相对最大误差结果如图 15所示。

图  15  ±0.4%、±0.8%、±1.2%、±1.6%和±2.0%的结果

Figure  15.  Inversion results of ±0.4%, ±0.8%, ±1.2%, ±1.6%, and ±2.0%

下载: 全尺寸图片 幻灯片

根据图 15所示可以看出小波分析后反演结果相对真值来对比的偏差不大，在0.4%~±2.0%的区间内反演相对误差大小与加入的随机噪声大小呈正相关，随机噪声越大则反演相对误差越大，其中±0.4%、±0.8%、±1.2%的结果较好，均小于5%，完全可以作为有效实验。

5.   结语

研究采用QPSO随机优化算法来解决轴承热处理炉传热过程中传热系数反演估计问题。此外，通过对实际测量过程中的噪声干扰，采用基于小波分析的阈值去噪算法以减轻噪声对最终结果的影响。从算法的仿真数据可以看出：直接实验时反演值与真值非常接近，即使在目标值变化复杂情况下，仍然能取得较好的反演准确性；在进行叠加随机噪声模拟时，我们分别对每个组合±0.4%、±0.8%、±1.2%、±1.6%和±2.0%的随机噪声进行了测试，结果显示当随机噪声小于等于1.2%时反演结果比较准确。本文的研究有助于帮助我们了解轴承热处理炉表面热辐射信息的反演分析。