基于EKF和UKF算法非均匀介质热物性参数重建

文爽，齐宏，刘少斌，任亚涛，阮立明

（1 哈尔滨工业大学能源科学与工程学院，黑龙江哈尔滨150001； 2 工业和信息化部航空航天热物理重点实验室，黑龙江哈尔滨150001）

引 言

热传导是自然界中存在的最基本的换热方式之一，它广泛存在于各种换热器和热管理系统中，且传热过程受到热导率和比热容等热物性的控制。在工程应用中为了保证机械的运转效率、使用寿命和减少事故的发生，精确地掌握材料热物性对于传热分析非常重要。在大多数实际问题中由于温度的不均匀性会引起热导率和比热容的非均匀，当材料内部的温度梯度增加时热导率和比热容的变化将加剧。除此之外，材料的热导率和比热容有时也会随着位置的不同而不同。因此，确定与温度或位置相关的热导率和比热容往往比较困难。

在过去的几十年里，随温度或位置变化的热物性参数重建引起了越来越多的关注，并且大量的数值算法被提出以解决热传导反问题[1-14]。Sawaf 等[1]使用Levenberg-Marquardt(L-M)方法重建了二维各向同性固体介质中恒定的热导率和比热容。Ukrainczyk[2]利用L-M 方法反演了琼脂水凝胶、甘油和渥太华石英砂等复杂介质的热扩散率，并获得精度较高的重建结果。Huang 等[3-6]利用共轭梯度法重构了一维和二维介质中随温度变化的热导率和比热容，当使用共轭梯度法求解热传导反问题时不需要给出重构物理量的函数形式便可以获得所有离散点的热物性。上述算法都是基于梯度的算法，它们具有计算精度高、收敛速度快的优点；但是它们的计算过程复杂且对初值依赖严重，若想获得精度较高的重建结果必须给定较好的初始分布。

近几年，学者们提出了大量的智能算法并且成功地将其引入到传热领域。Raudenský 等[15]将遗传算法引入到传热领域，并利用它求解了一维钢板中的热导率和比热容。Ardakani 等[16]利用微粒群算法(PSO)识别了二维固体介质中内含物的热导率和几何形状，并且研究了测量数据和内含物位置对重建结果的影响。Vakili 等[17]利用基本PSO、改进的排斥微粒群算法(RPSO)和完全改进排斥微粒群算法(CRPSO)对不同维度下的热传导反问题进行了研究，并指出当利用CRPSO 求解导热反问题时能获得最好的重建结果。智能算法具有较好的全局收敛性，算法原理简单易于编程实现；但是它的计算效率较低并且重建的参数个数较少。

Kalman[18]于1960 年提出了卡尔曼滤波算法，它最早被用于阿波罗登月计划中的轨道预测问题，之后被广泛地用于惯性导航、雷达追踪和图像识别等领域。近几年，Scarpa 等[19]将卡尔曼滤波算法引入到传热领域。Ji 等[20-23]将卡尔曼滤波算法与递归最小二乘结合，以边界温度为测量信号，重建了一维线性导热问题中的瞬态边界热流。LeBreux 等[24-26]雇佣基本卡尔曼滤波、扩展卡尔曼滤波和无迹卡尔曼滤波，对相变过程中的固液交界面进行了实时重建，他们指出无迹卡尔曼滤波的性能最好。由于采用非嵌入式的测量方式和实时在线重建的特点，卡尔曼滤波算法和其改进算法在热传导领域得到了广泛的应用。

当前对卡尔曼滤波的研究主要集中在利用它实时估计介质边界上的瞬态热流。本文主要采用扩展卡尔曼滤波算法和无迹卡尔曼滤波算法对介质中随时间和位置变化的热导率进行重建。

1 正问题模型

本文重点研究一维平板介质的热传导问题，并假设介质各向同性且最初的温度分布是均匀的。热导率随着介质温度或位置改变。介质的左边界受到高强度的热流加热，右边界处于自然对流换热的环境中，一维导热模型如图1所示。

图1 一维介质物理模型Fig.1 Schematic of physical model

直角坐标系中一维导热方程可以表示为

其中，ρ表示介质密度；cp表示介质比热容；T(x,t)表示t时刻介质x点的温度；λ表示介质热导率。

相应的初始条件可以表示为

相应的边界条件可以表示为

其中，qin(t)表示左边界上的瞬态热流；h 表示对流传热系数；Ts和Tw分别表示环境温度和右壁面温度。

如果所有的初始条件、边界条件都已知，将能获得任意时刻下介质的温度分布。该温度场被视为测量信息以重建介质的热物性。

2 反问题算法

标准卡尔曼滤波算法是一种适用于线性时变系统的最优递推估计算法，主要用于状态量的实时估计。随后研究者们提出了改进的扩展卡尔曼滤波(EKF)和无迹卡尔曼滤波算法(UKF)，将卡尔曼滤波推广到非线性时变系统的状态估计。同时郝晓静等[27-28]发现EKF 和UKF 可以用于参数辨识。当EKF 和UKF 用于参数辨识问题的研究时，需引入状态增广矩阵，增广项即为被识别的参数。利用卡尔曼滤波求解实际问题时，该问题的数学模型可以用随机微分方程描述；除此之外，还要求系统的过程噪声和测量噪声必须是零均值的高斯白噪声。

在一个非线性动态系统中，过程方程和测量方程可以被表示为式（5）～式（6）

其中，Φ 和H 表示非线性向量；w(k)表示过程噪声，即模型误差；v(k)表示测量噪声，即热电偶的测量误差。协方差矩阵可以用如式（7）～式（8）表示

其中，Q 和Q 分别表示过程噪声的协方差矩阵和协方差；同理，R 和R 分别表示测量噪声的协方差矩阵和协方差。

2.1 扩展卡尔曼滤波

相对于仅适用于线性系统的标准卡尔曼滤波算法，EKF 被广泛应用于非线性系统的状态估计和参数识别。EKF 通过在参考状态附近对状态方程和测量方程进行泰勒展开并忽略高阶项，来达到线性化的目的，其线性化示意图如图2 所示。但是，EKF 在线性化的过程中忽略了高阶非线性信息，因此它是一种次优滤波。

扩展卡尔曼滤波算法的递归过程可以分为时间更新和测量更新两部分[19]。

图2 EKF和UKF的线性化原理Fig.2 Linearization principle of EKF and UKF technique

①时间更新 根据tk-1时刻的最优估计对tk时刻的状态量进行初步预测。

预测误差的协方差矩阵可以表示为

②测量更新 根据tk时刻的测量信息修正初步预测结果。

估计误差的协方差矩阵可以表示为

滤波增益K(k)可以用如下的方程表示

线性化状态方程和测量方程的过程中会产生雅克比矩阵Φk-1和Hk，其表达式如式(14)～式(15)

由于扩展卡尔曼滤波技术通过一阶近似，而达到对非线性方程进行线性化处理的目的，重建结果的精度和稳定性与状态向量和误差协方差的初值选取有关。

2.2 无迹卡尔曼滤波

2000 年Julier 等[29]将卡尔曼滤波与无迹变换(UT)相结合，提出了无迹卡尔曼滤波技术(UKF)。与EKF 技术通过对非线性方程进行近似处理来达到线性化的目的不同，UKF 技术通过非线性函数的概率密度分布进行线性化处理而实现线性化。UKF方法基本思想是：在保证变量的统计特性，即误差协方差和均值不变的情况下，确定一组点集(Sigma点集)，将非线性变换用于每个Sigma点集，得到变换后的Sigma 点集。计算变换后Sigma 点集的统计特性，即可近似看作原状态变量经过非线性变换后得到的统计特性，其线性化示意图如图2所示。

将比例修正算法和对称采样策略相结合，获得对称比例采样策略。最终采样公式和权重可以表示为如式(16)～式（17）

其中，λ=α2(n+κ)-n为调节参数，用于控制Sigma点与均值之间的距离。参数α决定了采样点在均值附近的分布程度，其取值范围一般为10-4≤α ≤1；参数β 描述了新信息的分布情况，对于正态高斯分布情况，β=2时能获得最优的结果；κ为比例参数，其取值通常为0或者3-n。

确定好采样点和权值大小之后，无迹卡尔曼滤波的详细计算步骤如下所示。

①利用状态方程传递采样点

②通过预测的采样点信息χi(k+1),权值和计算协方差矩阵 P[ k + 1 k ] 和预测均值[ k + 1 k ]

③通过测量方程计算各预测采样点的测量信息

④计算测量信息均值和协方差矩阵

其中，Py^[k + 1]表示测量向量协方差矩阵，PX^y^[k + 1]为状态向量与测量向量的协方差矩阵。

⑤计算无迹卡尔曼滤波增益，并更新状态向量和误差协方差矩阵

利用无迹卡尔曼滤波算法求解非线性问题时，能获得至少精确到二阶矩的协方差和后验均值。

3 结果与讨论

本节设置了三个不同的算例，用于验证扩展卡尔曼滤波和无迹卡尔曼滤波算法在热物性参数重建中的可行性。算例1假设介质的热导率为一恒定值，既不随时间变化也不随空间变化；算例2中假设介质的热导率为空间位置的函数；算例3 中假设热导率为温度的函数。分别利用UKF 或EKF 算法对上述三种情况的热导率进行重建。

采用正问题模型求得的介质左右边界或内部的温度信息作为测量信号，重建介质的热物性参数。作用于介质左边界时变的热流可以通过式（28）获得

以介质右边界的温度信息作为测量信号，采用UKF 和EKF 算法分别对算例1 进行反演研究，重建结果如图3所示。其中介质的几何尺寸、比热容、测量噪声、测量噪声协方差、过程噪声协方差、时间步长和空间节点数分别为：Lx= 0.02、cp= 7980、σR=10、R=100、Q=0.0001、Δt=0.1 和E=11。从图3 中可以得到，UKF 和EKF 算法均获得了较为准确的热导率，证明了当前两种算法求解导热反问题的可行性。

深入分析算法中相关参数对重建结果准确性和稳定性的影响，有利于实验的开展。采用不同量级的测量噪声协方差，测试当前算法的准确性和稳定性。不同测量噪声协方差下的重构结果如图4所示，可以看出，随着测量噪声协方差的减小，热导率能够更快地逼近真值，即得到合理重建结果的时间减小。从式(25)和式(29)可以看出，随着测量噪声协方差的减小，卡尔曼滤波增益增加，这意味着更多的有效信息来自于当前时刻的真实测量Z(k+1)而不是先验信息Z^ [k + 1]。因此，当使用UKF 求解导热反问题时，为了获得比较理想的重建结果，建议选取较小的测量误差协方差。

图3 热导率重建结果Fig.3 Reconstructed results of thermal conductivity

在实际问题中，介质的热物性参数不是恒定不变的，它们往往随着介质的位置变化，或者是温度的函数。因此，首先对随位置变化的热导率进行反演，其表达式如式（29）

其中，a1=20 W/(m·K)-1，b2=300 W/(m·K2)-1。

介质左边界上的时变热流可以通过式(30)获得

图4 不同测量误差协方差下热导率重建结果Fig.4 Retrieval results of thermal conductivity with different measurement error covariance

以介质左右边界和中间点的温度信息作为测量信号，采用UKF 算法对算例2 进行反演研究。其中介质的测量噪声协方差、过程噪声协方差和空间节点数分别为：R=100、Q=0.0001 和E=21，其他的参数设置与算例1相同。添加不同测量误差的情况下重建介质中随位置变化的热导率，其重建结果如图5 所示，可以看出，当添加不同的测量噪声时，随着测量误差的增加，重建结果的波动性增强；但是，基于UKF 算法均能获得较为合理的热导率重建结果。

在现有工业问题中，除多层结构的材料外，材料的热物性一般不随材料的位置变化。但是，绝大多数介质的热物性会随着温度的不同而发生改变，因此重构介质内部随温度变化的热导率是一项非常有意义的工作。假设介质内部的热导率为温度的函数，且其可以被表示为式(31)

图5 不同测量误差下随位置变化的热导率重建结果Fig.5 Reconstructed results of space-dependent thermal conductivity with different measurement errors

图6 测量信息量对随温度变化的热导率重建研究的影响Fig.6 Retrieval results of thermal conductivity with different quantities of measurement signals

其中,a2=20 W·(m·K)-1，b2=0.001 W·(m·K2)-1。

采用EKF 算法对算例3 进行重建研究，其中测量噪声协方差R=100、状态向量中增广项所对应的过程噪声协方差分别为1 和0.1，其他的参数设置与算例1相同。考虑分别以介质左右边界以及中间点的温度信息作为测量信号(即M=3)，和仅仅以介质左右边的温度信息作为测量信号(M=2)，重建介质中随温度变化的热导率。从图6 中可以看出，基于EKF 算法均可获得较为合理的重建结果；但是，当仅采用介质左右边界的温度信息作为测量信号时，反演结果与真实热物性间的误差较大，且重建结果产生较为严重的时间滞后。这是由于EKF 算法是一种实时在线重建算法，每一时刻的重建结果主要与当前时刻的测量信息有关，前面时刻的测量信息对其影响较小；因此仅以左右边界的温度作为测量信息时，可用的测量信息较少，对于多个参数同时重建问题其不能准确地捕捉参数变化，最终导致热物性参数不能精确重建。

4 结 论

当前研究采用EKF、UKF 算法研究了一维介质中热物性参数重构问题。采用有限差分法求解一维热传导问题，根据正问题获得的介质左右边界和中心点温度信息，分别采用EKF 和UKE 算法反演了介质中均匀的热导率。最后研究了介质中随位置和温度变化的热导率重构问题。分析了测量误差、测量噪声协方差和测量信息多少对重建结果影响。得出以下主要结论。

(1)UKF 算法可以较为准确重建介质中随位置变化的热导率。

(2) EKF 算法可以较为准确重建介质中随时间变化的热导率。

(3)为获得精度较高的反演结果，建议选取较小的测量噪声协方差。

(4)UKF算法具有较强的鲁棒性，及时添加较大的测量误差时，也能获得较为精确的重建结果。

符 号 说 明

a——热导率的系数，W·(m·K)-1

b——热导率的系数，W·(m·K2)-1

cp——比热容，J·(kg·K)-1

H,Hk——分别为测量矩阵和线性化后的测量矩阵

h——对流传热系数，W·(m2·K)-1

K(k)——滤波增益

k——第k时刻

Lx——介质厚度，m

P,PX^y^,Py^——分别为状态估计向量的误差协方差，状态-测量交叉协方差矩阵，预测测量协方差矩阵

Q——过程噪声协方差矩阵

Q——过程噪声协方差

qin——边界时变热流，W·m-2

R——测量噪声协方差矩阵

R——测量噪声协方差

T,T0,Ts,Tw——分别为温度，初始温度，环境温度，壁面温度，K

t——时间，s

v——测量噪声，K

W——标量权重

w——过程噪声，W·m-2

X——状态向量

x——x轴坐标，m

Z——测量信息，K

α——与权重相关参数

β——与权重相关参数

κ——比例参数

λ——热导率，W·(m·K)-1

ρ——密度，kg·m-3

Φ,Φk-1——分别为状态转移矩阵和线性化后的状态转移矩阵

φ——预测测量

χ——Sigma点集