热电偶动态响应的带外部输入自回归模型*

金敏俊,李文军,郑永军,曾九孙

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

工业领域中热电偶被广泛应用在各种温度测量系统中[1]。随着测量技术的发展,热电偶也越来越多地应用在各种介质的瞬态温度测量中[2-3],尤其是流体介质瞬态温度的测量[4]。当热电偶被用来测量瞬态温度时,需要对热电偶进行动态响应特性评价和校准。随着动态测量技术的进步,产生了多种评价和校准方法。

近年CANUTO E等分别从理论和实验两个角度出发[5],提出了动态相对校准的方法,讨论了多个温度传感器组合测量和校准,包括产生精确温度梯度问题、保持均匀稳定温度场问题。赵学敏等采用火焰为温度激励源[6],对气体燃烧环境下的热电偶动态响应特性进行了分析。DINIZ A.C.G.C.A.等提出一种转动式校准装置和校准方法[7],提供气流环境下的温度阶跃,对热电偶进行动态校准并做不确定度分析,该方法计算了温度对时间二阶导数的极值,以极值对应的时间点为阶跃发生的时间。

在采用上述方法考察热电偶动态响应时,都需要确定温度激励产生的时间起点。用辅助测量手段确定时间起点,会引入热电偶输入端的误差。此外,目前的激励方法中,施加的激励温度被视为一个稳定的准确值,但是如果采用流体比如气体作为激励介质,气流温度在中心值附近存在波动。因此,当以温度激励作为热电偶的输入,以热电偶响应为输出,对热电偶动态响应进行统计建模时[8],由于输入带有误差[9],所构成的系统是一种变量带误差EIV(Error in Variables)系统[10],这种系统的辨识问题属于变量带误差辨识问题[11]。本文从统计建模角度,建立了双气流环境下的热电偶动态响应实验系统,考察热电偶在高低温两种气流之间摆动时的动态响应。给出了描述热电偶动态响应的状态空间方程,采用一种带外部输入自回归模型ARX(Autoregressive with exogenous,)对热电偶动态响应过程进行辨识[12]。以一种露端式镍铬镍硅(Nickel-chromium/Nickel-silicon)热电偶为实验对象,给出了一个算例。

1 热电偶动态响应的实验方法及实验系统

1.1 热电偶动态响应的实验方法及影响因素

热电偶动态响应实验方法主要包括热风洞法、电加热法、激波管法、投入法和激光加热法等。这些方法都是对热电偶产生一个近似于理想阶跃的温度激励,通过测量热电偶温度响应随时间的变化,然后评价热电偶动态响应能力。由于热电偶测量接点与介质之间的换热受到多种换热因素的影响,热电偶动态响应也相应地受到这些换热因素的影响[13]。

以热风洞法为例分析,热风洞法产生温度激励的介质为气体,影响热电偶与气体介质之间换热的主要因素包括:对流以及流动状态、介质物理性质以及换热面形状、接触面几何尺寸和相对位置[14]。如果温度激励幅值较大,则热电偶热端指向冷端的热传导也不可忽略。同时,热电偶与气体之间辐射换热、与管道内壁面之间也不可忽略。

如上所述,在气流环境下存在较多因素影响热电偶的动态响应。如果试图从物理过程建模,对这些因素进行选取和舍弃都比较困难。于是从实验数据角度建立模型成为一种合适的选择,即利用统计建模和数据分析来估计影响因素、以及影响因素的显著性。考虑到气体是测量动态温度的常见被测介质,因此采用常见的空气为激励介质,建立了一种空气环境下的热电偶动态响应实验系统。

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

根据上述分析,建立了一种气流环境下的热电偶动态响应实验系统。为了避免管道壁面辐射引起的影响,以热风枪气流出口为温度激励点。热风枪出口段加装了整流段,以使气流流动均匀。热电偶定位点如图1所示。

图1 热电偶在管道出口的换热模型

实验系统主要包括计算机、数据采集卡、可编程控制器、机械摆动驱动机构和热风枪。其示意图如图2所示。

图2 热电偶动态响应实验系统示意图

实验原理是利用热风枪A、热风枪B分别产生高低温气流,用可编程控制器控制机械摆动驱动机构,控制并驱动热电偶测量接点在两种气流之间摆动,由数据采集卡采集热电偶输出,由可编程控制器返回热电偶进入气流的各时间点。

实验中,数据采集卡采用MCC USB1616HS高速数据采集卡,可编程控制器采用OMRON CPM1A型,热风枪采用QUICK2008型。

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

以一种K型露端式镍铬镍硅热电偶为实验对象,利用图2所示实验系统,按照表1给出的实验配置参数进行实验。

表1 动态响应实验配置参数

按照上述配置参数进行实验,获得了三支K型露端式镍铬镍硅热电偶的响应数据集。图3给出了低温气流为326 K且高温气流为515 K时,这三支热电偶的输入和输出。图4给出了低温气流温度为326 K且高温气流温度为752 K时,这三支热电偶所对应的输入和输出。

图3 气流温度326 K和515 K时热电偶的响应

图4 气流温度326 K和752 K时热电偶的响应

通过实验所获得的响应数据都带有测量误差,如果分别用来和ym(t)表示热电偶的输入和输出,则热电偶动态响应过程如图5所示。

图5 输入和输出带有误差的过程

图5中,um(t)为热电偶输入的测量值并含有干扰eu(t),ym(t)为热电偶输出的测量值并含有干扰ey(t),yz(t)为作用在热电偶输出上的加性干扰。干扰ey(t)和yz(t)可以合并处理,但eu(t)的存在,使得上述热电偶动态响应过程构成一类变量带误差回归问题[15]。因此,把热电偶视为一种带有随机误差的确定性系统,建立其状态空间模型,利用带外部输入自回归方法进行辨识。为了表达简便,本文以下的描述中热电偶的输入um(t)和输出ym(t)仍分别用u(t)和表示。

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

2.1 热电偶的连续时间状态空间模型

热电偶动态响应的实验,是由激励介质给予热电偶温度激励,这种激励包括正向温度阶跃和负向温度阶跃。热电偶受到温度激励后,响应随着时间连续变化。把热电偶视为一个单输入单输出连续时间系统,其状态空间方程可以描述为:

(1)

y(t)=cx(t)+du(t)

(2)

式(1)为热电偶状态方程,式(2)为热电偶输出方程。式中,x(t)为热电偶的状态变量,y(t)为热电偶输出,u(t)为热电偶输入,A∈Rn×n为系数矩阵,b∈R1×n、c∈R1×n、d∈R1×n都是系数矩阵。

把式(1)和式(2)写成解的形式,其描述为[25]:

(3)

(4)

在热电偶动态响应的实验中,需要通过采样方法获得热电偶的输入和输出。考虑到施加的温度激励通常为正向阶跃或者负向阶跃,即输入近似为方波信号,采用零阶保持器,利用阶跃响应不变变换把连续系统转化为离散系统。设采样周期为T,离散时间变量为k,采用零阶保持器,即:

u(t)=u(kT)

(5)

kT≤t≤(k+1)T

(6)

再定义:

x(k)=x(kT)

(7)

u(k)=u(kT)

(8)

y(k)=y(kT)

(9)

这时,式(3)和式(4)所表示的热电偶连续时间状态空间模型转化为离散时间状态空间模型:

x(k+1)=Gx(k)+fu(k)

(10)

y(k)=cx(k)+du(k)

(11)

式中:G=eAt,f=A-1[G-I]b,I为单位阵。

式(10)和式(11)给出了在确定的采样周期下,采用零阶保持器时,热电偶连续时间状态空间模型所对应的离散时间状态空间模型的形式。在热电偶动态响应实验中,需要通过采样获得离散的输入输出数据,因此使用离散时间系统建立模型。

2.2 热电偶的离散时间状态空间模型

把热电偶响应过程视为离散时间系统,其中离散时间变量用k表示,这时热电偶的状态空间模型可以表示为:

x(k+1)=Ax(k)+bu(k)

(12)

y(k)=cx(k)+du(k)

(13)

式(12)和式(13)中,x(k)为热电偶的状态变量,y(k)为热电偶的输出序列,u(k)为热电偶的输入序列,A∈Rn×n为常系数矩阵,b∈R1×n、c∈R1×n、d∈R1×n都是系数矩阵。

引入移位算子,即单位前移算子z:

zy(k)=y(k+1)

(14)

以及单位后移算子z-1:

z-1y(k)=y(k-1)

(15)

式(12)和式(13)可以表示为:

zx(k)=Ax(k)+bu(k)

(16)

y(k)=cx(k)+du(k)

(17)

把式(16)、式(17)表示为解的形式,有:

x(k)=(zI-A)-1bu(k)

(18)

y(k)=[c(zI-A)-1b+ud(k)

(19)

式(18)和式(19)是描述热电偶动态响应过程的离散时间系统。由于输入和输出都存在测量误差,热电偶动态响应过程可用随机扰动与确定性模型结合的模型来描述,即利用随机系统模型对离散时间随机信号的统计性质进行评估。

2.3 热电偶的变量带误差自回归模型

对于带随机信号的确定性模型,其典型结构有自回归模型AR(Autoregressive model)、带外部输入自回归模型ARX(Autoregressive model with exogenous input)、带外部输入自回归滑动平均ARMAX(Autoregressive Moving Average model with exogenous input)、输出误差模型OE(Out Error model)和博克斯-詹金斯模型BJ(Box-Jenkins model)等多种。对于离散时间系统,用多项式表示这几种模型的一般表达式为:

(20)

式(20)中,A(z),B(z),C(z),D(z)和F(z)是用移位算子z表示的多项式,ui表示第i个输入,nu表示输入数目,nki表示第i个输入的时延,n(k)为随机干扰。

在上述几种模型结构中,ARX模型在工业中得到广泛应用,可用于变量带误差问题的处理[16-17]。ARX模型的一般表达式为:

(21)

对于热电偶动态响应过程,如果忽略热电偶进出气流的移动速度等次要影响因素,式(21)可简化为单输入单输出ARX模型[18]:

A(z)y(k)=B(z)u(k-nk)+n(k)

(22)

即:

(23)

式(23)中,y(k)为热电偶输出序列,A(z)和B(z)是用移位算子z表示的多项式,u(k)为热电偶输入序列,nk表示输入时延,n(k)表示干扰序列。定义na为输出时延长度,nb为输入时延长度,则多项式A(z)可表示为常系数多项式:

A(z)=1+a1z-1+a2z-2+…+anaz-na

(24)

多项式B(z)可表示为:

B(z)=b1z-1+b2z-2+…+bnbz-nb

(25)

对应的模型示意图如图6。

图6 带外部输入自回归模型

图7 变量带误差系统的带外部输入自回归模型

考虑热电偶输入u(k)带有干扰,即变量带误差情形时,式(23)应表示为:

(26)

u(k)=u0(k)+eu(k)

(27)

对应的ARX模型如图7所示。

对于图7所示的模型,当用线性模型时,定义参数向量为:

θ=[a1,a2,…,ana,b1,b2,…,bnb]T

(28)

定义数据向量为:

φ=[-y(k-1),-y(k-2),…,-y(k-na),u(k-nk),u(k-nk-1),…,u(k-nk-nb+1)]T

(29)

并定义:

(30)

则由式(26)和式(27)可得:

(31)

(32)

式(31)包含热电偶当前输出温度y(k)值,以及输出时延值、输入当前值和输入时延值。把当前输出温度y(k)的预报值记为yp(k),则由式(31)可构造当前值的一种预报值,其展开式可以表示为:

yp(k)=[-a1,-a2,…,-ana,b1,b2,…,bnb]× [y(k-1),y(k-2),…,y(k-na),u(k-nk),u(k-nk-1),…,u(k-nk-nb+1)]T

(33)

式(33)中,时延值y(k-1),y(k-2),…,y(k-na)和u(k-nk),u(k-nk-1),…,u(k-nk-nb+1)为回归量,-a1,-a2,…,-ana,b1,b2,…,bnb为回归量的系数,表示回归量的权重。

当用非线性函数f预报y(k),记:

yp(k)=f[y(k-1),y(k-2),…,y(k-na),u(k-nk),u(k-nk-1),…,u(k-nk-nb+1)]

(34)

式(34)构成了预报yp(k)的非线性ARX模型。

式(33)和式(34)可以统一表示为图8所示的模型结构,其中左侧矩形框表示计算估计量所采用的回归量,右侧矩形框表示计算估计量所采用的非线性函数和线性函数。

图8 ARX模型的回归量和估计量

图8也描述了预报yp(k)的步骤。其中第一个步骤是计算回归量。对于不同类型的热电偶,可以采用不同类型的回归量。如果采用形如u(k)、u(k-1)、y(k-1)的形式,是线性形式。如果采用形如u(k-1)×y(k-1)的形式,是u(k-1)和y(k-1)所构成的非线性形式。

第二个步骤是用回归量预报估计量。对于非线性函数,仍用φ表示第一个步骤中的回归量,图8中计算估计量的非线性函数的一般表达式表示如下[19]:

(35)

θ=[αk,βk,γk]k=1,2,…,d

(36)

式(35)和式(36)中,θ为参数向量,φ为向量形式的回归量,κ是非线性函数,αk、βk和γk为参数。本文算例中非线性函数κ采用了工业上常用的 Sigmoid 函数,函数形式为[20]:

κ(x)=1/(1+e-x)

(37)

2.4 准则函数和评价指标

采用式(31)所表示的线性模型和式(35)所表示的非线性模型,能够对热电偶动态响应数据做有限数据长度的离线辨识。对于线性模型或非线性模型,或者两者的混合模型,y的单步预报值都可表示为[21]:

y(k)=f(θ,φ)+e(k)

(38)

式(38)中,θ为参数向量,φ为数据向量,且有:

φ=[y(k-1),y(k-2),…,y(k-na),u(k-nk),u(k-nk-1),…,u(k-nk-nb+1)]T

(39)

由式(39),单步预报误差为:

e(k)=y(k)-f(θ,φ)

(40)

定义准则函数为:

J(θ)=y(k)-f(θ,φ)2

(41)

再对准则函数J(θ)极小化,有:

(42)

根据上述准则函数,利用最小二乘法计算参数θ。在收敛情形下,可获得最优的θ和对应的模型。

在获得模型后,用拟合度FIT对模型进行评价,具体形式如下:

(43)

同时计算模型的均方根误差RMSE,具体形式如下:

(44)

3 仿真

3.1 仿真方案

为了验证上述模型和算法,先通过仿真的方法对一种已知热电偶模型进行辨识。仿真分为两种辨识方法,第一种方法是用传统的自回归模型AR进行辨识[22],第二种方法是用带外部输入自回归模型ARX进行辨识。

假设已知热电偶模型用以下传递函数描述:

(45)

根据式(45),用伪二值序列作为输入生成无噪声输出,并用数据集zr表示。数据集zr的时间步长为0.2,长度为501,如图9所示。

图9 已知热电偶模型的无噪声输出和输入

上述伪二值序列增加服从高斯分布N(0,0.012)的随机噪声,作为含噪声的输入,并由式(45)产生对应的输出,作为含噪声的输出,并用数据集zm表示,如图10所示。

图10 已知热电偶模型的含噪声输出和输入

对于图9中的无噪声输出数据集zr,划分为两个部分,前334个数据为zr1,后167个数据为zr2;对于图10中的含噪声输出数据集zm,同样划分为两个部分,前334个数据为zm1,后167个数据为zm2。

3.2 自回归模型

首先利用自回归模型AR对含噪声数据zm1做辨识。通过计算,确定AR模型的阶数为4,对应的回归量为:y(k-1)、y(k-2)、y(k-3)、y(k-4)。图11 和图12给出了对应的计算结果,其中图11为自回归模型AR输出与数据zr1的比较,图12为自回归模型AR输出与数据zr2的比较。

图11 自回归模型AR输出与数据zr1

图12 自回归模型AR输出与数据zr2

从图11和图12中可以看出,AR模型对zr1的拟合度FIT为53.22%,均方根误差RMSE为0.075 0;AR模型对zr2的拟合度FIT为64.81%,均方根误差RMSE为0.063 5。

3.3 带外部输入自回归模型

在利用自回归模型AR辨识的基础上,再利用带外部输入自回归模型ARX对含噪声数据zm1做辨识。先通过线性ARX模型计算,得到na为5,nb为4,nk为1,确定回归量为:y(k-1)、y(k-2)、y(k-3)、y(k-4)、y(k-5)以及u(k-1)、u(k-2)、u(k-3)、u(k-4)。再根据线性模型所确定的回归量,采用线性和非线性ARX混合模型做辨识,其中κ函数采用κ(x)=1/(1+e-x)函数。在选取估计量时,线性部分和非线性部分使用相同的回归量,混合后的形式为:

f(φ)=LT(φ-r)+d+κ(Q(φ-r))

(46)

式(46)中,φ为作为回归量的数据向量,L为向量形式的线性回归系数,r为φ的均值,d为标量形式的偏移量,Q为投影矩阵。上述ARX模型的计算,模型参数包含了L,r,d,Q以及表征κ函数的参数。图13和图14给出了对应的计算结果。其中图13 为ARX模型输出与数据zr1的比较,图14为ARX模型输出与数据zr2的比较。

图13 带外部输入自回归模型ARX 输出与数据zr1

图14 带外部输入自回归模型ARX 输出与数据zr2

计算得到ARX模型对zr1的拟合度FIT为94.69%,均方根误差RMSE为0.008 5,ARX模型对zr2的拟合度FIT为95.96%,均方根误差RMSE为0.007 3。

比较AR模型和ARX模型的辨识结果可以看出,ARX模型的拟合度比AR模型高,并且ARX模型的均方根误差RMSE比AR模型小。

4 算例

在仿真的基础上,对图3和图4所示热电偶动态响应实验数据做ARX辨识。

把图3中的实验数据做归一化处理并平均,再把数据集划分为两个部分,取前3 500个数据z1做辨识,取后1 750个数据z2做验证。对数据集z1,采用线性模型计算得到na为5,nb为5,nk为1,确定回归量为:y(k-1)、y(k-2)、y(k-3)、y(k-4)、y(k-5)以及u(k-1)、u(k-2)、u(k-3)、u(k-4)、u(k-5)。根据线性模型所确定的回归量,再采用线性和非线性ARX混合模型计算,其中κ函数采用κ(x)=1/(1+e-x)函数。图15和图16给出了对应的计算结果,其中图15为ARX模型输出与数据z1的比较,图16为ARX模型输出与数据z2的比较。

图15 带外部输入自回归模型ARX输出与数据z1

图16 带外部输入自回归模型ARX输出与数据z2

计算得到ARX模型对z1的拟合度FIT为 89.07%,均方根误差RMSE为0.032 2,ARX模型对z2的拟合度FIT为88.92%,均方根误差RMSE为0.031 9。

类似地,对于图4中的实验数据做归一化处理并平均,把数据集划分为两个部分,取前3 500个数据z3做辨识,取后1 750个数据z4做验证。对数据集z3,采用线性模型计算得到na为5,nb为4,nk为1,确定回归量为:y(k-1)、y(k-2)、y(k-3)、y(k-4)、y(k-5)以及u(k-1)、u(k-2)、u(k-3)、u(k-4)。根据线性模型所确定的回归量,再采用线性和非线性ARX混合模型计算,其中κ函数采用κ(x)=1/(1+e-x)函数。图17和图18给出了对应的计算结果,其中图17为ARX模型输出与数据z3的比较,图18为ARX模型输出与数据z4的比较。

图18 带外部输入自回归模型ARX输出与数据z4

图17 带外部输入自回归模型 ARX输出与数据z3

计算得到ARX模型对z3的拟合度FIT为91.67%,均方根误差RMSE为0.024 4,ARX模型对z4的拟合度FIT为91.64%,均方根误差RMSE为0.025 0。

综合图15～图18可以看出,对于图3所示双气流之间温度差较小的情形,以及图4所示双气流之间温度差较大的情形,ARX模型都得到了较好的拟合度FIT以及较小的均方根误差RMSE。

5 结论

建立了一种双气流环境下热电偶动态响应实验系统,对热电偶产生正负温度激励并采集了热电偶响应数据。建立了热电偶动态响应的状态空间方程,把热电偶动态响应过程化简为随机加确定性系统,采用带外部输入自回归模型,利用线性函数和非线性函数逼近热电偶动态响应函数,处理了热电偶动态响应回归分析中存在的变量带误差问题。实验和计算结果表明,带外部输入自回归模型适用于K型露端式热电偶动态响应过程的回归分析。上述工作为工业热电偶动态响应能力评价提供了一种测试方法。这种方法还可应用于热电偶测量动态温度的递推估计中。