熱電偶動態(tài)響應模型的貝葉斯推斷*

李文軍,孫宏健,鄭永軍

(中國計量大學計量測試工程學院,杭州 310018)

熱電偶被廣泛應用于工業(yè)領域。當熱電偶用于測量動態(tài)溫度時,需要對熱電偶動態(tài)響應特性進行評價[1-2]。熱電偶動態(tài)響應特性評價方法通常遵循溫度傳感器動態(tài)特性評價方法,簡化為熱電偶時間常數(shù)的測試[3]。

熱電偶時間常數(shù)的測試方法有多種。根據(jù)溫度激勵產(chǎn)生方式的不同,常用方法有熱風洞法、電加熱法、激波管法、投入法和激光加熱法[4]。在這些測試方法中,都存在一個比較困難的環(huán)節(jié),即如何確定溫度激勵起始時間點。對于投入法,通過激光對射開關等輔助測量手段,測量投入介質(zhì)時的時間點并作為起始時間[5]。對于激光加熱法,是用響應更快的紅外探測設備輔助測量確定時間點[6]。但是這兩種輔助測量手段本身又引入了誤差,當評價小時間常數(shù)快響應熱電偶時,這種誤差影響變大[7-8]。

因此需要尋找其他的熱電偶動態(tài)特性評價方法,比如采用神經(jīng)網(wǎng)絡技術(shù)對校準數(shù)據(jù)進行建模[9]。類似地,如果熱電偶測溫視為一個動態(tài)過程,則可以采用辨識和回歸技術(shù)建模以考察熱電偶的動態(tài)響應。建模可以從理論建模入手,也可以從實驗建模入手,也可以從兩者的結(jié)合入手。本文從實驗建模的角度,建立了熱電偶動態(tài)響應實驗系統(tǒng),采集了熱電偶動態(tài)響應數(shù)據(jù)。把熱電偶響應數(shù)據(jù)視為時間序列,采用貝葉斯回歸方法估計參數(shù)分布?？紤]到實驗數(shù)據(jù)樣本容量的不足,采用了馬爾科夫鏈蒙特卡羅法MCMC(Markov Chain Monte Carlo)抽樣。

1 熱電偶動態(tài)響應及其實驗評價

1.1 熱電偶動態(tài)響應的影響因素

熱電偶動態(tài)響應受到多種因素影響。從實驗現(xiàn)象看,熱電偶動態(tài)特性受熱流密度影響較大[10-11]。另外,熱電偶在使用中本身也存在溫度漂移,以裸露式K型熱電偶為例,在室溫到600 ℃區(qū)間內(nèi),塞貝克系數(shù)會產(chǎn)生變化,變化引起的溫度誤差在0.5%[12]。以鎧裝式K型熱電偶為例,在室溫到600 ℃區(qū)間內(nèi),溫度的漂移在2.5 ℃左右[13]。

從物理過程分析,當熱電偶測量流體溫度時,影響熱電偶與流體之間非穩(wěn)態(tài)換熱的因素有:對流以及流動狀態(tài)、流體的物理性質(zhì)以及換熱面形狀、幾何尺寸和相對位置[14]。如果流體為氣體,則還存在熱電偶與氣體之間的輻射換熱,以及熱電偶的熱端指向冷端的熱傳導[15],如圖1所示。

圖1 熱電偶測氣流溫度時的換熱模型

當熱電偶測量固體表面溫度時,影響熱電偶與固體表面之間非穩(wěn)態(tài)換熱的因素有:接觸壓力、接觸面間隙介質(zhì)、接觸面表面粗糙度和接觸面形變[16]。

因此,無論是測量流體溫度或者測量固體溫度,都存在多種因素影響熱電偶與被測物體之間的非穩(wěn)態(tài)換熱過程,進而影響到熱電偶的動態(tài)響應。從物理過程建模時,選取和舍棄影響因素以及建立條件假設都比較困難。于是考慮從測量數(shù)據(jù)入手建模,為此建立了熱電偶測量氣體溫度時的實驗評價系統(tǒng)。

圖2 實驗評價系統(tǒng)組成示意圖

1.2 熱電偶動態(tài)響應的實驗評價系統(tǒng)

熱電偶測量氣體溫度的實驗評價系統(tǒng)如圖2所示,主要包括計算機、數(shù)據(jù)采集卡、可編程控制器、往復運動驅(qū)動機構(gòu)、風道和工業(yè)調(diào)溫風機。實驗中,利用工業(yè)調(diào)溫風機產(chǎn)生熱氣流,利用可編程控制器控制往復運動驅(qū)動機構(gòu),驅(qū)動熱電偶進出風道,以獲得溫度激勵。熱電偶的輸出信號由數(shù)據(jù)采集卡返回給計算機。

1.3 鎳鉻鎳硅熱電偶的響應數(shù)據(jù)集

以一種K型鎳鉻鎳硅熱電偶為對象,利用實驗評價系統(tǒng)取得了幾種溫度激勵下的響應數(shù)據(jù),表1給出了一組實驗評價時的設置參數(shù)。

表1 實驗評價的參數(shù)設置

按照上述設置中激勵溫度進行實驗,獲得了熱電偶響應數(shù)據(jù)集。圖3給出了鎳鉻鎳硅熱電偶在3種溫度激勵下的響應。

圖3 鎳鉻鎳硅熱電偶溫度響應散點圖

由于通過實驗獲得各種溫度階躍下的測試數(shù)據(jù)并不容易,因此評價熱電偶響應能力要在測試數(shù)據(jù)不充分的條件下進行。當測試數(shù)據(jù)數(shù)量充分時,一般可采用矩法、極大似然法等經(jīng)典推斷方法,當測試數(shù)據(jù)不充分時,可采用小樣本的推斷方法[17]。把圖3中熱電偶溫度響應視為時間序列,可采用貝葉斯推斷方法,分析和評價熱電偶響應能力[18]。

2 熱電偶動態(tài)響應模型的貝葉斯推斷

2.1 貝葉斯推斷

熱電偶的動態(tài)響應可以視為一個動態(tài)系統(tǒng)。對于動態(tài)系統(tǒng),可以基于一組測量值或觀察值,用推斷方法來對參數(shù)做估計。推斷的具體描述是,把實驗中獲得的熱電偶動態(tài)響應數(shù)據(jù)T視為觀測樣本,對熱電偶動態(tài)響應的模型參數(shù)θ進行估計。經(jīng)典參數(shù)估計方法如矩法估計和最大似然估計,把參數(shù)θ視為參數(shù)空間Θ的一個未知但是存在的固定值。而從隨機過程視角看,參數(shù)θ是參數(shù)空間Θ中取值的隨機變量,可以采用貝葉斯推斷方法來估計參數(shù)θ的分布。

貝葉斯推斷的核心工具是貝葉斯準則,按照貝葉斯準則,參數(shù)和溫度響應數(shù)據(jù)服從條件概率:

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

(1)

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

(2)

按照貝葉斯準則有:

(3)

(4)

式(1)～(4)中:θ為模型參數(shù);T為熱電偶響應數(shù)據(jù);p(θ,T)為模型參數(shù)和數(shù)據(jù)的聯(lián)合概率密度函數(shù);p(T|θ)為似然函數(shù),即在給定一個模型時描述這個數(shù)據(jù)的似然;p(θ)為模型參數(shù)的先驗分布,可以視為看到數(shù)據(jù)之前對模型參數(shù)分布的一種描述;p(θ|T)為后驗分布;p(T)為邊界似然,是對所有似然概率的積分。

根據(jù)貝葉斯準則,式(3)和式(4)給出了用熱電偶動態(tài)響應模型的先驗知識對模型進行推斷的一種框架。按照這種框架,不直接估計參數(shù)值,而是考慮參數(shù)服從一定概率分布。貝葉斯推斷熱電偶響應模型參數(shù)的過程可以概括如下:

首先確定模型參數(shù)的先驗分布p(θ);其次獲取一系列觀測T;再次求出參數(shù)θ的后驗分布p(θ|T),求出θ的期望值作為其最終值。推斷過程中,還需要定義參數(shù)的一個方差來評估推斷的準確程度。

2.2 熱電偶動態(tài)響應的回歸模型

對于圖3所示的熱電偶動態(tài)響應數(shù)據(jù)集,為了計算方便采用過余溫度進行轉(zhuǎn)換,令:

(5)

圖4 鎳鉻鎳硅熱電偶無量綱溫度響應散點圖

圖4給出了通過轉(zhuǎn)換后鎳鉻鎳硅熱電偶無量綱溫度響應的散點圖,以時間序列的形式列出了3種激勵溫度下的響應。

圖4中的無量綱溫度時間序列包含了溫度激勵未產(chǎn)生時即熱電偶處于環(huán)境溫度時的狀態(tài)值。當評價的數(shù)據(jù)集包含這部分狀態(tài)值時,數(shù)據(jù)處理中可以避免確定溫度激勵產(chǎn)生時間點。根據(jù)圖4中散點的分布形狀,采用如下形式的函數(shù)做回歸:

(6)

式中:t為時間,單位為s;α,β為回歸系數(shù)。

式(6)是一個非線性回歸模型,可以先做線性化變化,令:

(7)

t*=t

(8)

則式(6)轉(zhuǎn)化為:

T*=α+βt*

(9)

式(9)是一個線性回歸函數(shù),其中α為截距,β為回歸系數(shù)。為簡便起見,下文中省略*號,仍用T表示經(jīng)過式(5)無量綱化和經(jīng)過式(6)線性化后的溫度,于是回歸函數(shù)表示為:

T=α+βt

(10)

再設噪聲方差為σ2,并假設溫度為獨立高斯分布[19],即:

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

(11)

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

把μi表示為時間預報值ti和參數(shù)α和β的函數(shù),有:

(12)

假設α和β的先驗分布都是高斯分布,即:

(13)

(14)

為方便,記噪聲方差σ2的對數(shù)為待估計參數(shù),并假設其為高斯分布,即:

logσ2～N(κ,ω2)

(15)

于是要推斷的參數(shù)可以表示為:

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

(16)

式(10)確定了回歸函數(shù)形式和待估計參數(shù),式(13)～式(15)確定了待估計參數(shù)的先驗。

由于后驗分布通常不可解,所以一般采用各種近似貝葉斯推理方法,比如變分方法和MCMC方法。變分方法的思想是把一個待解問題,引入一個變量轉(zhuǎn)變成一個優(yōu)化問題。變分方法在優(yōu)化問題中引入了約束,讓問題簡化以達到快速求解,但是如果引入的約束比較嚴格,近似的程度就會變差。而MCMC方法的思想是基于隨機模擬,即構(gòu)造一個隨機過程來逐漸逼近所求的分布,通過隨機過程不斷采樣達到刻畫目標分布的目的。本文選擇了MCMC方法來估計后驗分布。

2.3 雜交蒙特卡羅抽樣方法

在貝葉斯估計中,可以采用MCMC方法來估計后驗分布。MCMC方法的思想是通過生成隨機樣本并用樣本來估計后驗分布。MCMC生成樣本的有多種,其中雜交蒙特卡羅HMC(Hamiltonian Monte Carlo)方法是一種快速抽樣方法。雜交蒙特卡羅方法的基本思想是在MCMC抽樣中利用分子動力學模擬MD(Molecular Dynamics)方法來產(chǎn)生馬爾科夫鏈移動[20]。用這個方法所產(chǎn)生的建議分布,往往與目標分布的動力學機理更接近[21]。

從方法上看,分子動力學模擬是對漢密爾頓方程做積分求解來估計復雜物理系統(tǒng)的典型特征,這種典型特征通常是組態(tài)函數(shù)(一般是時間的函數(shù))關于時間的平均。蒙特卡羅模擬則是從由系統(tǒng)勢能和溫度所決定的玻爾茲曼分布中抽取典型組態(tài)的樣本。分子動力學模擬和蒙特卡羅模擬之間的聯(lián)系在于,物理系統(tǒng)中時間的平均會收斂到狀態(tài)的平均。因此,分子動力學模擬產(chǎn)生的結(jié)果往往和蒙特卡洛模擬產(chǎn)生的估計通常會有較好的吻合。

把雜交蒙特卡羅抽樣方法應用于熱電偶響應分布的抽樣,其具體過程描述如下。

假設從分布π(T)∝exp[-U(T)]中抽取樣本,首先引入和T共軛的虛擬動量變量p,其次引入一個引導漢密爾頓函數(shù)H′(T,p):

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

(17)

它用來產(chǎn)生一個建議狀態(tài)。

再定義接收漢密爾頓函數(shù)H(T,p):

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

(18)

它用來決定是接受還是拒絕。

在定義了引導漢密爾頓函數(shù)和接收漢密爾頓函數(shù)后,雜交蒙特卡羅抽樣通過迭代過程實現(xiàn),迭代過程用算法1表示。

算法1:

假設t時刻在組態(tài)空間處于T即T(t)=T,然后在t+1時刻執(zhí)行以下步驟:

第1步 從高斯分布中產(chǎn)生一個新動量p,即p～?(p)∝exp{-k(p)};

第2步 從狀態(tài)(T,p)出發(fā),運行Leap-frog算法L步,在相空間中得到一個新狀態(tài)(T′,p′);

第3步 以概率min{1,exp{-H(T′,p′)+H(T,p)}}接收狀態(tài)(T′,p′),以剩余概率令T(t+1)=T。用上述過程可獲得樣本。具體抽樣過程中,還需要確定抽樣的起始點和步長。

3 算例

3.1 雜交蒙特卡羅法抽樣的起始點

對于圖4中的無量綱溫度時間序列,選擇式(10)為回歸公式,選擇參數(shù)為式(13)、式(14)和式(15)形式的先驗分布,并預置參數(shù)的均值和標準差,預置值如表2所示。

表2 參數(shù)先驗分布均值和標準差的預置值

為計算方便,對式(3)兩邊取自然對數(shù),有:

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

(19)

再引入輔助函數(shù)g(θ):

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

(20)

對于所建立的線性回歸方程(10),利用上述輔助函數(shù)g(θ)計算似然函數(shù)與先驗乘積,并求乘積的對數(shù)。

在進行雜交蒙特卡羅抽樣之前,需要確定抽樣起始點?？紤]到抽樣的效率,通常使用最大后驗估計值MAP(Maximum a Posteriori Estimation)為起始點。為此,還需要先計算θmap,即:

(21)

計算θmap結(jié)合了式(20)和式(21),即通過計算g(θ)的梯度得到,同步得到其所對應的參數(shù)估計值。表3給出了算例中3個參數(shù)的最大后驗估計值所對應的參數(shù)值。

表3 參數(shù)最大后驗估計值對應的參數(shù)值

在確定了起始點后,抽樣過程中還需要調(diào)整步長。調(diào)整步長通過一個自適應抽樣過程完成。

3.2 雜交蒙特卡羅法抽樣的步長調(diào)整

在熱電偶響應的參數(shù)后驗分布抽樣過程中,為了能提高抽樣效率,通過自適應方法調(diào)整雜交蒙特卡羅抽樣的步長以及相應的虛擬動量向量p。通常把算法1中第3步的目標接受率初始值設為0.65,并把調(diào)整步長的起始點也定在最大后驗估計值這點。當?shù)^程收斂時可以得到目標步長。圖5給出了算例中通過一個自適應過程確定步長的收斂過程。其中接受率終值為0.64。

圖5 自適應抽樣中步長的收斂過程

3.3 熱電偶動態(tài)響應模型的回歸結(jié)果

在確定了選擇起始點和步長的方法后,使用獨立的馬爾科夫鏈從后驗分布中抽取樣本。抽取樣本時選擇互異的馬爾科夫鏈初始點,且初始點隨機分布在最大后驗估計值θmap周圍。利用老化(burn-in)過程,舍棄開始階段的老化樣本,獲得所使用樣本。在得到熱電偶響應模型的參數(shù)后驗分布后,再利用標準MCMC診斷給出評價[22]。表4給出了算例中熱電偶響應參數(shù)回歸結(jié)果和結(jié)果檢驗數(shù)據(jù)。

表4 回歸結(jié)果和結(jié)果檢驗

在表4中,Mean為參數(shù)的后驗均值;MCSE為Monte Carlo standard error,即蒙特卡羅標準誤,是后驗均值的標準差;SD為Standard error,即后驗標準差;Q5為邊際后驗分布的5分位;Q95為邊際后驗分布的95分位;RHat為Gelman-Rubin 收斂診斷值,用來刻畫平行的馬爾科夫鏈的鏈內(nèi)與鏈間方差的差異,當收斂診斷值小于等于1.1則說明抽樣算法收斂[23]。在表5中看到,算例中的RHat值都小于1.1,可認為收斂到了期望的分布。

圖6 參數(shù)的軌跡圖

3.4 抽樣樣本的分布

在回歸收斂的情形下,考察算例中的抽樣樣本是否構(gòu)成一個目標分布上的合理集合。圖6分別給出了算例中3個參數(shù)的第1個馬爾科夫鏈軌跡,其中已經(jīng)舍棄了老化樣本。從圖6可以看出,樣本與樣本之間沒有長期相關性,樣本混合比較合理。

圖7給出了算例中3個參數(shù)的馬爾科夫鏈混合后的散點圖和直方圖。其中直方圖給出了參數(shù)的邊際后驗分布。

圖7 參數(shù)的邊際后驗分布

如果采用參數(shù)的后驗均值來表示,算例中的回歸公式可以表示如下:

(22)

式(22)是在氣流環(huán)境下,環(huán)境溫度為29 ℃,激勵溫度區(qū)間為95 ℃～200 ℃時,一種K型鎳鉻鎳硅熱電偶無量綱溫度動態(tài)響應的回歸方程。

4 結(jié)論

針對熱電偶動態(tài)響應能力的評價,建立了用貝葉斯推斷方法計算參數(shù)分布的一種框架。對一種鎳鉻鎳硅熱電偶,用貝葉斯推斷方法計算了其動態(tài)響應模型的參數(shù)分布。上述工作為熱電偶動態(tài)響應能力的評價提供了一種測試方法。同時,利用上述框架,可以根據(jù)熱電偶先驗信息和所建回歸模型進行預測,然后利用測量值和模型的后驗信息更新預測值,再把后驗信息置為新的先驗信息再次進行預測。因此,這種推斷過程還可應用于熱電偶測量動態(tài)溫度的遞推估計中。

[1] 吳朋,林濤. 基于QGA_SVM的鎧裝熱電偶傳感器辨識建模研究[J]. 儀器儀表學報,2014,35(2):343-349.

[2] 吳德會,趙偉,黃松嶺,等. 傳感器動態(tài)建模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] 馬勤弟,雷敏. 薄膜熱電偶的動態(tài)校準及辨識建模[J]. 儀器儀表學報,1999,20(3):300-302.

[5] 趙學敏,王文廉,李巖峰,等. 火焰溫度場測試中的傳感器動態(tài)響應研究[J]. 傳感技術(shù)學報,2016,29(3):368-372.

[6] 郝曉劍,李科杰,劉健,等. 基于CO2激光器的溫度傳感器可溯源動態(tài)校準[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] 趙學敏,王文廉,李巖峰,等. 火焰溫度場測試中的傳感器動態(tài)特性補償[J]. 傳感技術(shù)學報,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] 成紅剛,陳雄,鞠玉濤,等. 基于熱電偶動態(tài)特性的溫度預估方法實驗研究[J]. 固體火箭技術(shù),2012,35(6):842-846.

[11] 曾磊,桂業(yè)偉,賀立新,等. 鍍層式同軸熱電偶數(shù)據(jù)處理方法研究[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算法的貝葉斯統(tǒng)計方法[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.