基于EKF和UKF算法非均勻介質熱物性參數重建

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

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

引 言

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

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

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

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

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

1 正問題模型

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

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

直角坐標系中一維導熱方程可以表示為

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

相應的初始條件可以表示為

相應的邊界條件可以表示為

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

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

2 反問題算法

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

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

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

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

2.1 擴展卡爾曼濾波

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

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

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

①時間更新 根據tk-1時刻的最優估計對tk時刻的狀態量進行初步預測。

預測誤差的協方差矩陣可以表示為

②測量更新 根據tk時刻的測量信息修正初步預測結果。

估計誤差的協方差矩陣可以表示為

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

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

由于擴展卡爾曼濾波技術通過一階近似，而達到對非線性方程進行線性化處理的目的，重建結果的精度和穩定性與狀態向量和誤差協方差的初值選取有關。

2.2 無跡卡爾曼濾波

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

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

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

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

①利用狀態方程傳遞采樣點

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

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

④計算測量信息均值和協方差矩陣

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

⑤計算無跡卡爾曼濾波增益，并更新狀態向量和誤差協方差矩陣

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

3 結果與討論

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

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

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

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

圖3 熱導率重建結果Fig.3 Reconstructed results of thermal conductivity

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

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

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

圖4 不同測量誤差協方差下熱導率重建結果Fig.4 Retrieval results of thermal conductivity with different measurement error covariance

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

在現有工業問題中，除多層結構的材料外，材料的熱物性一般不隨材料的位置變化。但是，絕大多數介質的熱物性會隨著溫度的不同而發生改變，因此重構介質內部隨溫度變化的熱導率是一項非常有意義的工作。假設介質內部的熱導率為溫度的函數，且其可以被表示為式(31)

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

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

4 結 論

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

(1)UKF 算法可以較為準確重建介質中隨位置變化的熱導率。

(2) EKF 算法可以較為準確重建介質中隨時間變化的熱導率。

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

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

符 號 說 明

a——熱導率的系數，W·(m·K)-1

b——熱導率的系數，W·(m·K2)-1

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

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

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

K(k)——濾波增益

k——第k時刻

Lx——介質厚度，m

P,PX^y^,Py^——分別為狀態估計向量的誤差協方差，狀態-測量交叉協方差矩陣，預測測量協方差矩陣

Q——過程噪聲協方差矩陣

Q——過程噪聲協方差

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

R——測量噪聲協方差矩陣

R——測量噪聲協方差

T,T0,Ts,Tw——分別為溫度，初始溫度，環境溫度，壁面溫度，K

t——時間，s

v——測量噪聲，K

W——標量權重

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

X——狀態向量

x——x軸坐標，m

Z——測量信息，K

α——與權重相關參數

β——與權重相關參數

κ——比例參數

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

ρ——密度，kg·m-3

Φ,Φk-1——分別為狀態轉移矩陣和線性化后的狀態轉移矩陣

φ——預測測量

χ——Sigma點集