熱電偶動態響應的帶外部輸入自回歸模型*

金敏俊,李文軍,鄭永軍,曾九孫

(中國計量大學計量測試工程學院,杭州 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型露端式熱電偶動態響應過程的回歸分析。上述工作為工業熱電偶動態響應能力評價提供了一種測試方法。這種方法還可應用于熱電偶測量動態溫度的遞推估計中。