热电偶动态响应模型的贝叶斯推断*

李文军,孙宏健,郑永军

(中国计量大学计量测试工程学院,杭州 310018)

热电偶被广泛应用于工业领域。当热电偶用于测量动态温度时,需要对热电偶动态响应特性进行评价[1-2]。热电偶动态响应特性评价方法通常遵循温度传感器动态特性评价方法,简化为热电偶时间常数的测试[3]。

热电偶时间常数的测试方法有多种。根据温度激励产生方式的不同,常用方法有热风洞法、电加热法、激波管法、投入法和激光加热法[4]。在这些测试方法中,都存在一个比较困难的环节,即如何确定温度激励起始时间点。对于投入法,通过激光对射开关等辅助测量手段,测量投入介质时的时间点并作为起始时间[5]。对于激光加热法,是用响应更快的红外探测设备辅助测量确定时间点[6]。但是这两种辅助测量手段本身又引入了误差,当评价小时间常数快响应热电偶时,这种误差影响变大[7-8]。

因此需要寻找其他的热电偶动态特性评价方法,比如采用神经网络技术对校准数据进行建模[9]。类似地,如果热电偶测温视为一个动态过程,则可以采用辨识和回归技术建模以考察热电偶的动态响应。建模可以从理论建模入手,也可以从实验建模入手,也可以从两者的结合入手。本文从实验建模的角度,建立了热电偶动态响应实验系统,采集了热电偶动态响应数据。把热电偶响应数据视为时间序列,采用贝叶斯回归方法估计参数分布。考虑到实验数据样本容量的不足,采用了马尔科夫链蒙特卡罗法MCMC(Markov Chain Monte Carlo)抽样。

1 热电偶动态响应及其实验评价

1.1 热电偶动态响应的影响因素

热电偶动态响应受到多种因素影响。从实验现象看,热电偶动态特性受热流密度影响较大[10-11]。另外,热电偶在使用中本身也存在温度漂移,以裸露式K型热电偶为例,在室温到600 ℃区间内,塞贝克系数会产生变化,变化引起的温度误差在0.5%[12]。以铠装式K型热电偶为例,在室温到600 ℃区间内,温度的漂移在2.5 ℃左右[13]。

从物理过程分析,当热电偶测量流体温度时,影响热电偶与流体之间非稳态换热的因素有:对流以及流动状态、流体的物理性质以及换热面形状、几何尺寸和相对位置[14]。如果流体为气体,则还存在热电偶与气体之间的辐射换热,以及热电偶的热端指向冷端的热传导[15],如图1所示。

图1 热电偶测气流温度时的换热模型

当热电偶测量固体表面温度时,影响热电偶与固体表面之间非稳态换热的因素有:接触压力、接触面间隙介质、接触面表面粗糙度和接触面形变[16]。

因此,无论是测量流体温度或者测量固体温度,都存在多种因素影响热电偶与被测物体之间的非稳态换热过程,进而影响到热电偶的动态响应。从物理过程建模时,选取和舍弃影响因素以及建立条件假设都比较困难。于是考虑从测量数据入手建模,为此建立了热电偶测量气体温度时的实验评价系统。

图2 实验评价系统组成示意图

1.2 热电偶动态响应的实验评价系统

热电偶测量气体温度的实验评价系统如图2所示,主要包括计算机、数据采集卡、可编程控制器、往复运动驱动机构、风道和工业调温风机。实验中,利用工业调温风机产生热气流,利用可编程控制器控制往复运动驱动机构,驱动热电偶进出风道,以获得温度激励。热电偶的输出信号由数据采集卡返回给计算机。

1.3 镍铬镍硅热电偶的响应数据集

以一种K型镍铬镍硅热电偶为对象,利用实验评价系统取得了几种温度激励下的响应数据,表1给出了一组实验评价时的设置参数。

表1 实验评价的参数设置

按照上述设置中激励温度进行实验,获得了热电偶响应数据集。图3给出了镍铬镍硅热电偶在3种温度激励下的响应。

图3 镍铬镍硅热电偶温度响应散点图

由于通过实验获得各种温度阶跃下的测试数据并不容易,因此评价热电偶响应能力要在测试数据不充分的条件下进行。当测试数据数量充分时,一般可采用矩法、极大似然法等经典推断方法,当测试数据不充分时,可采用小样本的推断方法[17]。把图3中热电偶温度响应视为时间序列,可采用贝叶斯推断方法,分析和评价热电偶响应能力[18]。

2 热电偶动态响应模型的贝叶斯推断

2.1 贝叶斯推断

热电偶的动态响应可以视为一个动态系统。对于动态系统,可以基于一组测量值或观察值,用推断方法来对参数做估计。推断的具体描述是,把实验中获得的热电偶动态响应数据T视为观测样本,对热电偶动态响应的模型参数θ进行估计。经典参数估计方法如矩法估计和最大似然估计,把参数θ视为参数空间Θ的一个未知但是存在的固定值。而从随机过程视角看,参数θ是参数空间Θ中取值的随机变量,可以采用贝叶斯推断方法来估计参数θ的分布。

贝叶斯推断的核心工具是贝叶斯准则,按照贝叶斯准则,参数和温度响应数据服从条件概率:

p(θ,T)=p(T|θ)p(θ)

(1)

p(θ,T)=p(θ|T)p(T)

(2)

按照贝叶斯准则有:

(3)

(4)

式(1)～(4)中:θ为模型参数;T为热电偶响应数据;p(θ,T)为模型参数和数据的联合概率密度函数;p(T|θ)为似然函数,即在给定一个模型时描述这个数据的似然;p(θ)为模型参数的先验分布,可以视为看到数据之前对模型参数分布的一种描述;p(θ|T)为后验分布;p(T)为边界似然,是对所有似然概率的积分。

根据贝叶斯准则,式(3)和式(4)给出了用热电偶动态响应模型的先验知识对模型进行推断的一种框架。按照这种框架,不直接估计参数值,而是考虑参数服从一定概率分布。贝叶斯推断热电偶响应模型参数的过程可以概括如下:

首先确定模型参数的先验分布p(θ);其次获取一系列观测T;再次求出参数θ的后验分布p(θ|T),求出θ的期望值作为其最终值。推断过程中,还需要定义参数的一个方差来评估推断的准确程度。

2.2 热电偶动态响应的回归模型

对于图3所示的热电偶动态响应数据集,为了计算方便采用过余温度进行转换,令:

(5)

图4 镍铬镍硅热电偶无量纲温度响应散点图

图4给出了通过转换后镍铬镍硅热电偶无量纲温度响应的散点图,以时间序列的形式列出了3种激励温度下的响应。

图4中的无量纲温度时间序列包含了温度激励未产生时即热电偶处于环境温度时的状态值。当评价的数据集包含这部分状态值时,数据处理中可以避免确定温度激励产生时间点。根据图4中散点的分布形状,采用如下形式的函数做回归:

(6)

式中:t为时间,单位为s;α,β为回归系数。

式(6)是一个非线性回归模型,可以先做线性化变化,令:

(7)

t*=t

(8)

则式(6)转化为:

T*=α+βt*

(9)

式(9)是一个线性回归函数,其中α为截距,β为回归系数。为简便起见,下文中省略*号,仍用T表示经过式(5)无量纲化和经过式(6)线性化后的温度,于是回归函数表示为:

T=α+βt

(10)

再设噪声方差为σ2,并假设温度为独立高斯分布[19],即:

Ti|θ～N(μi(θ),σ2)

(11)

式中:μi表示高斯分布均值。

把μi表示为时间预报值ti和参数α和β的函数,有:

(12)

假设α和β的先验分布都是高斯分布,即:

(13)

(14)

为方便,记噪声方差σ2的对数为待估计参数,并假设其为高斯分布,即:

logσ2～N(κ,ω2)

(15)

于是要推断的参数可以表示为:

θ=θ(α,β,logσ2)

(16)

式(10)确定了回归函数形式和待估计参数,式(13)～式(15)确定了待估计参数的先验。

由于后验分布通常不可解,所以一般采用各种近似贝叶斯推理方法,比如变分方法和MCMC方法。变分方法的思想是把一个待解问题,引入一个变量转变成一个优化问题。变分方法在优化问题中引入了约束,让问题简化以达到快速求解,但是如果引入的约束比较严格,近似的程度就会变差。而MCMC方法的思想是基于随机模拟,即构造一个随机过程来逐渐逼近所求的分布,通过随机过程不断采样达到刻画目标分布的目的。本文选择了MCMC方法来估计后验分布。

2.3 杂交蒙特卡罗抽样方法

在贝叶斯估计中,可以采用MCMC方法来估计后验分布。MCMC方法的思想是通过生成随机样本并用样本来估计后验分布。MCMC生成样本的有多种,其中杂交蒙特卡罗HMC(Hamiltonian Monte Carlo)方法是一种快速抽样方法。杂交蒙特卡罗方法的基本思想是在MCMC抽样中利用分子动力学模拟MD(Molecular Dynamics)方法来产生马尔科夫链移动[20]。用这个方法所产生的建议分布,往往与目标分布的动力学机理更接近[21]。

从方法上看,分子动力学模拟是对汉密尔顿方程做积分求解来估计复杂物理系统的典型特征,这种典型特征通常是组态函数(一般是时间的函数)关于时间的平均。蒙特卡罗模拟则是从由系统势能和温度所决定的玻尔兹曼分布中抽取典型组态的样本。分子动力学模拟和蒙特卡罗模拟之间的联系在于,物理系统中时间的平均会收敛到状态的平均。因此,分子动力学模拟产生的结果往往和蒙特卡洛模拟产生的估计通常会有较好的吻合。

把杂交蒙特卡罗抽样方法应用于热电偶响应分布的抽样,其具体过程描述如下。

假设从分布π(T)∝exp[-U(T)]中抽取样本,首先引入和T共轭的虚拟动量变量p,其次引入一个引导汉密尔顿函数H′(T,p):

H′(T,p)=U′(T)+k(p)

(17)

它用来产生一个建议状态。

再定义接收汉密尔顿函数H(T,p):

H(T,p)=U(T)+k(p)

(18)

它用来决定是接受还是拒绝。

在定义了引导汉密尔顿函数和接收汉密尔顿函数后,杂交蒙特卡罗抽样通过迭代过程实现,迭代过程用算法1表示。

算法1:

假设t时刻在组态空间处于T即T(t)=T,然后在t+1时刻执行以下步骤:

第1步 从高斯分布中产生一个新动量p,即p～∅(p)∝exp{-k(p)};

第2步 从状态(T,p)出发,运行Leap-frog算法L步,在相空间中得到一个新状态(T′,p′);

第3步 以概率min{1,exp{-H(T′,p′)+H(T,p)}}接收状态(T′,p′),以剩余概率令T(t+1)=T。用上述过程可获得样本。具体抽样过程中,还需要确定抽样的起始点和步长。

3 算例

3.1 杂交蒙特卡罗法抽样的起始点

对于图4中的无量纲温度时间序列,选择式(10)为回归公式,选择参数为式(13)、式(14)和式(15)形式的先验分布,并预置参数的均值和标准差,预置值如表2所示。

表2 参数先验分布均值和标准差的预置值

为计算方便,对式(3)两边取自然对数,有:

logp(θ|T)=const.+logp(T|θ)+logp(θ)

(19)

再引入辅助函数g(θ):

g(θ)=logp(T|θ)+logp(θ)

(20)

对于所建立的线性回归方程(10),利用上述辅助函数g(θ)计算似然函数与先验乘积,并求乘积的对数。

在进行杂交蒙特卡罗抽样之前,需要确定抽样起始点。考虑到抽样的效率,通常使用最大后验估计值MAP(Maximum a Posteriori Estimation)为起始点。为此,还需要先计算θmap,即:

(21)

计算θmap结合了式(20)和式(21),即通过计算g(θ)的梯度得到,同步得到其所对应的参数估计值。表3给出了算例中3个参数的最大后验估计值所对应的参数值。

表3 参数最大后验估计值对应的参数值

在确定了起始点后,抽样过程中还需要调整步长。调整步长通过一个自适应抽样过程完成。

3.2 杂交蒙特卡罗法抽样的步长调整

在热电偶响应的参数后验分布抽样过程中,为了能提高抽样效率,通过自适应方法调整杂交蒙特卡罗抽样的步长以及相应的虚拟动量向量p。通常把算法1中第3步的目标接受率初始值设为0.65,并把调整步长的起始点也定在最大后验估计值这点。当迭代过程收敛时可以得到目标步长。图5给出了算例中通过一个自适应过程确定步长的收敛过程。其中接受率终值为0.64。

图5 自适应抽样中步长的收敛过程

3.3 热电偶动态响应模型的回归结果

在确定了选择起始点和步长的方法后,使用独立的马尔科夫链从后验分布中抽取样本。抽取样本时选择互异的马尔科夫链初始点,且初始点随机分布在最大后验估计值θmap周围。利用老化(burn-in)过程,舍弃开始阶段的老化样本,获得所使用样本。在得到热电偶响应模型的参数后验分布后,再利用标准MCMC诊断给出评价[22]。表4给出了算例中热电偶响应参数回归结果和结果检验数据。

表4 回归结果和结果检验

在表4中,Mean为参数的后验均值;MCSE为Monte Carlo standard error,即蒙特卡罗标准误,是后验均值的标准差;SD为Standard error,即后验标准差;Q5为边际后验分布的5分位;Q95为边际后验分布的95分位;RHat为Gelman-Rubin 收敛诊断值,用来刻画平行的马尔科夫链的链内与链间方差的差异,当收敛诊断值小于等于1.1则说明抽样算法收敛[23]。在表5中看到,算例中的RHat值都小于1.1,可认为收敛到了期望的分布。

图6 参数的轨迹图

3.4 抽样样本的分布

在回归收敛的情形下,考察算例中的抽样样本是否构成一个目标分布上的合理集合。图6分别给出了算例中3个参数的第1个马尔科夫链轨迹,其中已经舍弃了老化样本。从图6可以看出,样本与样本之间没有长期相关性,样本混合比较合理。

图7给出了算例中3个参数的马尔科夫链混合后的散点图和直方图。其中直方图给出了参数的边际后验分布。

图7 参数的边际后验分布

如果采用参数的后验均值来表示,算例中的回归公式可以表示如下:

(22)

式(22)是在气流环境下,环境温度为29 ℃,激励温度区间为95 ℃～200 ℃时,一种K型镍铬镍硅热电偶无量纲温度动态响应的回归方程。

4 结论

针对热电偶动态响应能力的评价,建立了用贝叶斯推断方法计算参数分布的一种框架。对一种镍铬镍硅热电偶,用贝叶斯推断方法计算了其动态响应模型的参数分布。上述工作为热电偶动态响应能力的评价提供了一种测试方法。同时,利用上述框架,可以根据热电偶先验信息和所建回归模型进行预测,然后利用测量值和模型的后验信息更新预测值,再把后验信息置为新的先验信息再次进行预测。因此,这种推断过程还可应用于热电偶测量动态温度的递推估计中。

[1] 吴朋,林涛. 基于QGA_SVM的铠装热电偶传感器辨识建模研究[J]. 仪器仪表学报,2014,35(2):343-349.

[2] 吴德会,赵伟,黄松岭,等. 传感器动态建模FLANN方法改进研究[J]. 仪器仪表学报,2009,30(2):362-367.

[3] Doghmane M Y,Lanzetta F,Gavignet E. Dynamic Characterization of a Transient Surface Temperature Sensor[J]. Procedia Engineering,2015,120:1245-1248.

[4] 马勤弟,雷敏. 薄膜热电偶的动态校准及辨识建模[J]. 仪器仪表学报,1999,20(3):300-302.

[5] 赵学敏,王文廉,李岩峰,等. 火焰温度场测试中的传感器动态响应研究[J]. 传感技术学报,2016,29(3):368-372.

[6] 郝晓剑,李科杰,刘健,等. 基于CO2激光器的温度传感器可溯源动态校准[J]. 兵工学报,2009,30(2):156-159.

[7] Diniz A C G C A,Almeida M E K D,Vianna J N S,et al. Methodology for Estimating Measurement Uncertainty in the Dynamic Calibration of Industrial Temperature Sensors[J]. Journal of the Brazilian Society of Mechanical Sciences and Engineering,2017,39:1053-1060.

[8] 赵学敏,王文廉,李岩峰,等. 火焰温度场测试中的传感器动态特性补偿[J]. 传感技术学报,2017,30(5):735-741.

[9] Danisman K,Dalkiran I,Celebi F V. Design of A High Precision Temperature Measurement System Based on Artificial Neural Network for Different Thermocouple Types[J]. Measurement,2006,39:695-700.

[10] 成红刚,陈雄,鞠玉涛,等. 基于热电偶动态特性的温度预估方法实验研究[J]. 固体火箭技术,2012,35(6):842-846.

[11] 曾磊,桂业伟,贺立新,等. 镀层式同轴热电偶数据处理方法研究[J]. 工程热物理学报,2009,30(4):661-664.

[12] Webster E S. Drift in Type K Bare-Wire Thermocouples from Different Manufacturers[J]. International Journal of Thermophysics,2017,38(5):1-14.

[13] Webster E S. Thermal Preconditioning of MIMS Type K Thermocouples to Reduce Drift[J]. International Journal of Thermophysics,2017,38(1):1-14.

[14] Jaremkiewicz M. Accurate Measurement of Unsteady State Fluid Temperature[J]. Heat Mass Transfer,2017,53:887-897.

[15] 杨永军,蔡静. 特殊条件下的温度测量[M]. 北京:中国计量出版社,2009:130-132.

[16] Verma N N,Mazumder S. Extraction of Thermal Contact Conductance of Metal-Metal Contacts from Scale-Resolved Direct Numerical Simulation[J]. International Journal of Heat and Mass Transfer,2016,94:164-173.

[17] 姚继涛,王旭东. 小样本条件下可变作用代表值的贝叶斯推断方法[J]. 西南交通大学学报,2014,49(6):995-1001.

[18] Kantz H,Schreiber T. Nonlinear Time Series Analysis[M]. UK:Cambridge University Press,2004:22-145.

[19] Chakraborty S,Das P K. Application of Bayesian Inference Technique for the Reconstruction of An Isothermal Hot Spot inside A Circular Disc from Peripheral Temperature Measurement—A Critical Assess-ment[J]. International Journal of Heat and Mass Transfer,2015,88:456-469.

[20] 刘军. 科学计算中的蒙特卡罗策略[M]. 北京:高等教育出版社,2009:140-141.

[21] 康崇禄. 蒙特卡罗方法理论和应用[M]. 北京:科学出版社,2015:115-116.

[22] 刘金山. 基于MCMC算法的贝叶斯统计方法[M]. 北京:科学出版社,2016:49-53.

[23] Brooks S P,Gelman A. General Methods for Monitoring Convergence of Iterative Simulations[J]. Journal of Computational and Graphical Statistics,1998,7(4):434-455.