基于EKF和UKF算法非均勻介質(zhì)熱物性參數(shù)重建

文爽，齊宏，劉少斌，任亞濤，阮立明

（1 哈爾濱工業(yè)大學(xué)能源科學(xué)與工程學(xué)院，黑龍江哈爾濱150001； 2 工業(yè)和信息化部航空航天熱物理重點實驗室，黑龍江哈爾濱150001）

引 言

熱傳導(dǎo)是自然界中存在的最基本的換熱方式之一，它廣泛存在于各種換熱器和熱管理系統(tǒng)中，且傳熱過程受到熱導(dǎo)率和比熱容等熱物性的控制。在工程應(yīng)用中為了保證機械的運轉(zhuǎn)效率、使用壽命和減少事故的發(fā)生，精確地掌握材料熱物性對于傳熱分析非常重要。在大多數(shù)實際問題中由于溫度的不均勻性會引起熱導(dǎo)率和比熱容的非均勻，當(dāng)材料內(nèi)部的溫度梯度增加時熱導(dǎo)率和比熱容的變化將加劇。除此之外，材料的熱導(dǎo)率和比熱容有時也會隨著位置的不同而不同。因此，確定與溫度或位置相關(guān)的熱導(dǎo)率和比熱容往往比較困難。

在過去的幾十年里，隨溫度或位置變化的熱物性參數(shù)重建引起了越來越多的關(guān)注，并且大量的數(shù)值算法被提出以解決熱傳導(dǎo)反問題[1-14]。Sawaf 等[1]使用Levenberg-Marquardt(L-M)方法重建了二維各向同性固體介質(zhì)中恒定的熱導(dǎo)率和比熱容。Ukrainczyk[2]利用L-M 方法反演了瓊脂水凝膠、甘油和渥太華石英砂等復(fù)雜介質(zhì)的熱擴散率，并獲得精度較高的重建結(jié)果。Huang 等[3-6]利用共軛梯度法重構(gòu)了一維和二維介質(zhì)中隨溫度變化的熱導(dǎo)率和比熱容，當(dāng)使用共軛梯度法求解熱傳導(dǎo)反問題時不需要給出重構(gòu)物理量的函數(shù)形式便可以獲得所有離散點的熱物性。上述算法都是基于梯度的算法，它們具有計算精度高、收斂速度快的優(yōu)點；但是它們的計算過程復(fù)雜且對初值依賴嚴(yán)重，若想獲得精度較高的重建結(jié)果必須給定較好的初始分布。

近幾年，學(xué)者們提出了大量的智能算法并且成功地將其引入到傳熱領(lǐng)域。Raudensky 等[15]將遺傳算法引入到傳熱領(lǐng)域，并利用它求解了一維鋼板中的熱導(dǎo)率和比熱容。Ardakani 等[16]利用微粒群算法(PSO)識別了二維固體介質(zhì)中內(nèi)含物的熱導(dǎo)率和幾何形狀，并且研究了測量數(shù)據(jù)和內(nèi)含物位置對重建結(jié)果的影響。Vakili 等[17]利用基本PSO、改進的排斥微粒群算法(RPSO)和完全改進排斥微粒群算法(CRPSO)對不同維度下的熱傳導(dǎo)反問題進行了研究，并指出當(dāng)利用CRPSO 求解導(dǎo)熱反問題時能獲得最好的重建結(jié)果。智能算法具有較好的全局收斂性，算法原理簡單易于編程實現(xiàn)；但是它的計算效率較低并且重建的參數(shù)個數(shù)較少。

Kalman[18]于1960 年提出了卡爾曼濾波算法，它最早被用于阿波羅登月計劃中的軌道預(yù)測問題，之后被廣泛地用于慣性導(dǎo)航、雷達追蹤和圖像識別等領(lǐng)域。近幾年，Scarpa 等[19]將卡爾曼濾波算法引入到傳熱領(lǐng)域。Ji 等[20-23]將卡爾曼濾波算法與遞歸最小二乘結(jié)合，以邊界溫度為測量信號，重建了一維線性導(dǎo)熱問題中的瞬態(tài)邊界熱流。LeBreux 等[24-26]雇傭基本卡爾曼濾波、擴展卡爾曼濾波和無跡卡爾曼濾波，對相變過程中的固液交界面進行了實時重建，他們指出無跡卡爾曼濾波的性能最好。由于采用非嵌入式的測量方式和實時在線重建的特點，卡爾曼濾波算法和其改進算法在熱傳導(dǎo)領(lǐng)域得到了廣泛的應(yīng)用。

當(dāng)前對卡爾曼濾波的研究主要集中在利用它實時估計介質(zhì)邊界上的瞬態(tài)熱流。本文主要采用擴展卡爾曼濾波算法和無跡卡爾曼濾波算法對介質(zhì)中隨時間和位置變化的熱導(dǎo)率進行重建。

1 正問題模型

本文重點研究一維平板介質(zhì)的熱傳導(dǎo)問題，并假設(shè)介質(zhì)各向同性且最初的溫度分布是均勻的。熱導(dǎo)率隨著介質(zhì)溫度或位置改變。介質(zhì)的左邊界受到高強度的熱流加熱，右邊界處于自然對流換熱的環(huán)境中，一維導(dǎo)熱模型如圖1所示。

圖1 一維介質(zhì)物理模型Fig.1 Schematic of physical model

直角坐標(biāo)系中一維導(dǎo)熱方程可以表示為

其中，ρ表示介質(zhì)密度；cp表示介質(zhì)比熱容；T(x,t)表示t時刻介質(zhì)x點的溫度；λ表示介質(zhì)熱導(dǎo)率。

相應(yīng)的初始條件可以表示為

相應(yīng)的邊界條件可以表示為

其中，qin(t)表示左邊界上的瞬態(tài)熱流；h 表示對流傳熱系數(shù)；Ts和Tw分別表示環(huán)境溫度和右壁面溫度。

如果所有的初始條件、邊界條件都已知，將能獲得任意時刻下介質(zhì)的溫度分布。該溫度場被視為測量信息以重建介質(zhì)的熱物性。

2 反問題算法

標(biāo)準(zhǔn)卡爾曼濾波算法是一種適用于線性時變系統(tǒng)的最優(yōu)遞推估計算法，主要用于狀態(tài)量的實時估計。隨后研究者們提出了改進的擴展卡爾曼濾波(EKF)和無跡卡爾曼濾波算法(UKF)，將卡爾曼濾波推廣到非線性時變系統(tǒng)的狀態(tài)估計。同時郝曉靜等[27-28]發(fā)現(xiàn)EKF 和UKF 可以用于參數(shù)辨識。當(dāng)EKF 和UKF 用于參數(shù)辨識問題的研究時，需引入狀態(tài)增廣矩陣，增廣項即為被識別的參數(shù)。利用卡爾曼濾波求解實際問題時，該問題的數(shù)學(xué)模型可以用隨機微分方程描述；除此之外，還要求系統(tǒng)的過程噪聲和測量噪聲必須是零均值的高斯白噪聲。

在一個非線性動態(tài)系統(tǒng)中，過程方程和測量方程可以被表示為式（5）～式（6）

其中，Φ 和H 表示非線性向量；w(k)表示過程噪聲，即模型誤差；v(k)表示測量噪聲，即熱電偶的測量誤差。協(xié)方差矩陣可以用如式（7）～式（8）表示

其中，Q 和Q 分別表示過程噪聲的協(xié)方差矩陣和協(xié)方差；同理，R 和R 分別表示測量噪聲的協(xié)方差矩陣和協(xié)方差。

2.1 擴展卡爾曼濾波

相對于僅適用于線性系統(tǒng)的標(biāo)準(zhǔn)卡爾曼濾波算法，EKF 被廣泛應(yīng)用于非線性系統(tǒng)的狀態(tài)估計和參數(shù)識別。EKF 通過在參考狀態(tài)附近對狀態(tài)方程和測量方程進行泰勒展開并忽略高階項，來達到線性化的目的，其線性化示意圖如圖2 所示。但是，EKF 在線性化的過程中忽略了高階非線性信息，因此它是一種次優(yōu)濾波。

擴展卡爾曼濾波算法的遞歸過程可以分為時間更新和測量更新兩部分[19]。

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

①時間更新 根據(jù)tk-1時刻的最優(yōu)估計對tk時刻的狀態(tài)量進行初步預(yù)測。

預(yù)測誤差的協(xié)方差矩陣可以表示為

②測量更新 根據(jù)tk時刻的測量信息修正初步預(yù)測結(jié)果。

估計誤差的協(xié)方差矩陣可以表示為

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

線性化狀態(tài)方程和測量方程的過程中會產(chǎn)生雅克比矩陣Φk-1和Hk，其表達式如式(14)～式(15)

由于擴展卡爾曼濾波技術(shù)通過一階近似，而達到對非線性方程進行線性化處理的目的，重建結(jié)果的精度和穩(wěn)定性與狀態(tài)向量和誤差協(xié)方差的初值選取有關(guān)。

2.2 無跡卡爾曼濾波

2000 年Julier 等[29]將卡爾曼濾波與無跡變換(UT)相結(jié)合，提出了無跡卡爾曼濾波技術(shù)(UKF)。與EKF 技術(shù)通過對非線性方程進行近似處理來達到線性化的目的不同，UKF 技術(shù)通過非線性函數(shù)的概率密度分布進行線性化處理而實現(xiàn)線性化。UKF方法基本思想是：在保證變量的統(tǒng)計特性，即誤差協(xié)方差和均值不變的情況下，確定一組點集(Sigma點集)，將非線性變換用于每個Sigma點集，得到變換后的Sigma 點集。計算變換后Sigma 點集的統(tǒng)計特性，即可近似看作原狀態(tài)變量經(jīng)過非線性變換后得到的統(tǒng)計特性，其線性化示意圖如圖2所示。

將比例修正算法和對稱采樣策略相結(jié)合，獲得對稱比例采樣策略。最終采樣公式和權(quán)重可以表示為如式(16)～式（17）

其中，λ=α2(n+κ)-n為調(diào)節(jié)參數(shù)，用于控制Sigma點與均值之間的距離。參數(shù)α決定了采樣點在均值附近的分布程度，其取值范圍一般為10-4≤α ≤1；參數(shù)β 描述了新信息的分布情況，對于正態(tài)高斯分布情況，β=2時能獲得最優(yōu)的結(jié)果；κ為比例參數(shù)，其取值通常為0或者3-n。

確定好采樣點和權(quán)值大小之后，無跡卡爾曼濾波的詳細(xì)計算步驟如下所示。

①利用狀態(tài)方程傳遞采樣點

②通過預(yù)測的采樣點信息χi(k+1),權(quán)值和計算協(xié)方差矩陣 P[ k + 1 k ] 和預(yù)測均值[ k + 1 k ]

③通過測量方程計算各預(yù)測采樣點的測量信息

④計算測量信息均值和協(xié)方差矩陣

其中，Py^[k + 1]表示測量向量協(xié)方差矩陣，PX^y^[k + 1]為狀態(tài)向量與測量向量的協(xié)方差矩陣。

⑤計算無跡卡爾曼濾波增益，并更新狀態(tài)向量和誤差協(xié)方差矩陣

利用無跡卡爾曼濾波算法求解非線性問題時，能獲得至少精確到二階矩的協(xié)方差和后驗均值。

3 結(jié)果與討論

本節(jié)設(shè)置了三個不同的算例，用于驗證擴展卡爾曼濾波和無跡卡爾曼濾波算法在熱物性參數(shù)重建中的可行性。算例1假設(shè)介質(zhì)的熱導(dǎo)率為一恒定值，既不隨時間變化也不隨空間變化；算例2中假設(shè)介質(zhì)的熱導(dǎo)率為空間位置的函數(shù)；算例3 中假設(shè)熱導(dǎo)率為溫度的函數(shù)。分別利用UKF 或EKF 算法對上述三種情況的熱導(dǎo)率進行重建。

采用正問題模型求得的介質(zhì)左右邊界或內(nèi)部的溫度信息作為測量信號，重建介質(zhì)的熱物性參數(shù)。作用于介質(zhì)左邊界時變的熱流可以通過式（28）獲得

以介質(zhì)右邊界的溫度信息作為測量信號，采用UKF 和EKF 算法分別對算例1 進行反演研究，重建結(jié)果如圖3所示。其中介質(zhì)的幾何尺寸、比熱容、測量噪聲、測量噪聲協(xié)方差、過程噪聲協(xié)方差、時間步長和空間節(jié)點數(shù)分別為：Lx= 0.02、cp= 7980、σR=10、R=100、Q=0.0001、Δt=0.1 和E=11。從圖3 中可以得到，UKF 和EKF 算法均獲得了較為準(zhǔn)確的熱導(dǎo)率，證明了當(dāng)前兩種算法求解導(dǎo)熱反問題的可行性。

深入分析算法中相關(guān)參數(shù)對重建結(jié)果準(zhǔn)確性和穩(wěn)定性的影響，有利于實驗的開展。采用不同量級的測量噪聲協(xié)方差，測試當(dāng)前算法的準(zhǔn)確性和穩(wěn)定性。不同測量噪聲協(xié)方差下的重構(gòu)結(jié)果如圖4所示，可以看出，隨著測量噪聲協(xié)方差的減小，熱導(dǎo)率能夠更快地逼近真值，即得到合理重建結(jié)果的時間減小。從式(25)和式(29)可以看出，隨著測量噪聲協(xié)方差的減小，卡爾曼濾波增益增加，這意味著更多的有效信息來自于當(dāng)前時刻的真實測量Z(k+1)而不是先驗信息Z^ [k + 1]。因此，當(dāng)使用UKF 求解導(dǎo)熱反問題時，為了獲得比較理想的重建結(jié)果，建議選取較小的測量誤差協(xié)方差。

圖3 熱導(dǎo)率重建結(jié)果Fig.3 Reconstructed results of thermal conductivity

在實際問題中，介質(zhì)的熱物性參數(shù)不是恒定不變的，它們往往隨著介質(zhì)的位置變化，或者是溫度的函數(shù)。因此，首先對隨位置變化的熱導(dǎo)率進行反演，其表達式如式（29）

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

介質(zhì)左邊界上的時變熱流可以通過式(30)獲得

圖4 不同測量誤差協(xié)方差下熱導(dǎo)率重建結(jié)果Fig.4 Retrieval results of thermal conductivity with different measurement error covariance

以介質(zhì)左右邊界和中間點的溫度信息作為測量信號，采用UKF 算法對算例2 進行反演研究。其中介質(zhì)的測量噪聲協(xié)方差、過程噪聲協(xié)方差和空間節(jié)點數(shù)分別為：R=100、Q=0.0001 和E=21，其他的參數(shù)設(shè)置與算例1相同。添加不同測量誤差的情況下重建介質(zhì)中隨位置變化的熱導(dǎo)率，其重建結(jié)果如圖5 所示，可以看出，當(dāng)添加不同的測量噪聲時，隨著測量誤差的增加，重建結(jié)果的波動性增強；但是，基于UKF 算法均能獲得較為合理的熱導(dǎo)率重建結(jié)果。

在現(xiàn)有工業(yè)問題中，除多層結(jié)構(gòu)的材料外，材料的熱物性一般不隨材料的位置變化。但是，絕大多數(shù)介質(zhì)的熱物性會隨著溫度的不同而發(fā)生改變，因此重構(gòu)介質(zhì)內(nèi)部隨溫度變化的熱導(dǎo)率是一項非常有意義的工作。假設(shè)介質(zhì)內(nèi)部的熱導(dǎo)率為溫度的函數(shù)，且其可以被表示為式(31)

圖5 不同測量誤差下隨位置變化的熱導(dǎo)率重建結(jié)果Fig.5 Reconstructed results of space-dependent thermal conductivity with different measurement errors

圖6 測量信息量對隨溫度變化的熱導(dǎo)率重建研究的影響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 進行重建研究，其中測量噪聲協(xié)方差R=100、狀態(tài)向量中增廣項所對應(yīng)的過程噪聲協(xié)方差分別為1 和0.1，其他的參數(shù)設(shè)置與算例1相同。考慮分別以介質(zhì)左右邊界以及中間點的溫度信息作為測量信號(即M=3)，和僅僅以介質(zhì)左右邊的溫度信息作為測量信號(M=2)，重建介質(zhì)中隨溫度變化的熱導(dǎo)率。從圖6 中可以看出，基于EKF 算法均可獲得較為合理的重建結(jié)果；但是，當(dāng)僅采用介質(zhì)左右邊界的溫度信息作為測量信號時，反演結(jié)果與真實熱物性間的誤差較大，且重建結(jié)果產(chǎn)生較為嚴(yán)重的時間滯后。這是由于EKF 算法是一種實時在線重建算法，每一時刻的重建結(jié)果主要與當(dāng)前時刻的測量信息有關(guān)，前面時刻的測量信息對其影響較小；因此僅以左右邊界的溫度作為測量信息時，可用的測量信息較少，對于多個參數(shù)同時重建問題其不能準(zhǔn)確地捕捉參數(shù)變化，最終導(dǎo)致熱物性參數(shù)不能精確重建。

4 結(jié) 論

當(dāng)前研究采用EKF、UKF 算法研究了一維介質(zhì)中熱物性參數(shù)重構(gòu)問題。采用有限差分法求解一維熱傳導(dǎo)問題，根據(jù)正問題獲得的介質(zhì)左右邊界和中心點溫度信息，分別采用EKF 和UKE 算法反演了介質(zhì)中均勻的熱導(dǎo)率。最后研究了介質(zhì)中隨位置和溫度變化的熱導(dǎo)率重構(gòu)問題。分析了測量誤差、測量噪聲協(xié)方差和測量信息多少對重建結(jié)果影響。得出以下主要結(jié)論。

(1)UKF 算法可以較為準(zhǔn)確重建介質(zhì)中隨位置變化的熱導(dǎo)率。

(2) EKF 算法可以較為準(zhǔn)確重建介質(zhì)中隨時間變化的熱導(dǎo)率。

(3)為獲得精度較高的反演結(jié)果，建議選取較小的測量噪聲協(xié)方差。

(4)UKF算法具有較強的魯棒性，及時添加較大的測量誤差時，也能獲得較為精確的重建結(jié)果。

符 號 說 明

a——熱導(dǎo)率的系數(shù)，W·(m·K)-1

b——熱導(dǎo)率的系數(shù)，W·(m·K2)-1

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

H,Hk——分別為測量矩陣和線性化后的測量矩陣

h——對流傳熱系數(shù)，W·(m2·K)-1

K(k)——濾波增益

k——第k時刻

Lx——介質(zhì)厚度，m

P,PX^y^,Py^——分別為狀態(tài)估計向量的誤差協(xié)方差，狀態(tài)-測量交叉協(xié)方差矩陣，預(yù)測測量協(xié)方差矩陣

Q——過程噪聲協(xié)方差矩陣

Q——過程噪聲協(xié)方差

qin——邊界時變熱流，W·m-2

R——測量噪聲協(xié)方差矩陣

R——測量噪聲協(xié)方差

T,T0,Ts,Tw——分別為溫度，初始溫度，環(huán)境溫度，壁面溫度，K

t——時間，s

v——測量噪聲，K

W——標(biāo)量權(quán)重

w——過程噪聲，W·m-2

X——狀態(tài)向量

x——x軸坐標(biāo)，m

Z——測量信息，K

α——與權(quán)重相關(guān)參數(shù)

β——與權(quán)重相關(guān)參數(shù)

κ——比例參數(shù)

λ——熱導(dǎo)率，W·(m·K)-1

ρ——密度，kg·m-3

Φ,Φk-1——分別為狀態(tài)轉(zhuǎn)移矩陣和線性化后的狀態(tài)轉(zhuǎn)移矩陣

φ——預(yù)測測量

χ——Sigma點集