基于基函數逼近和卡爾曼濾波的溫度場重建*

顏 華，顧夢楠，王伊凡

(沈陽工業大學 信息科學與工程學院，沈陽 110870)

基于聲學CT的溫度場重建技術[1]是根據被測區域內多條路徑上的聲波飛行時間數據推算被測區域的溫度分布，屬于“由效果反求原因”的逆問題研究.該技術具有非接觸不干擾被測溫度場、測溫范圍廣、環境適應能力強及可在線測量等優點.

重建算法在聲學CT溫度場重建中起到至關重要的作用.為實現溫度場重建，通常要將被測區域劃分成若干個網格，每個網格重建一個溫度值，相當于一個原始溫度采樣點[2].典型的聲學CT溫度場重建算法按照正問題建模方式的不同主要有兩大類：一類利用網格內聲線路徑長度建模，以最小二乘法(LSM)[3-5]為代表；另一類利用基函數對溫度(聲速)分布逼近建模，以Markov徑向基函數Tikhonov正則化法(MTR)[6]為代表.LSM的優點是運算簡單，重建速度快，但該算法要求被測區域網格剖分總數小于可獲得的聲波數據數；MTR算法的優點是被測區域網格剖分數不受聲線數的限制，能夠更有效地重建復雜溫度場，但其重建質量仍然有較大的提升空間.卡爾曼濾波是一種實時遞推算法，其實質是以最小均方誤差為估計準則的最優估計，與最小二乘估計相比，卡爾曼濾波能夠有效提升估計精度，抗噪能力也更為突出[7-8].

本文提出一種聲學CT溫度場重建新算法，利用Markov徑向基函數逼近復雜的溫度(聲速)分布，建立正問題模型，利用卡爾曼濾波進行逆問題求解獲得聲速分布，進而利用溫度與聲速的關系得到溫度場分布.

1 聲學CT溫度場重建原理

聲波在氣體介質中的傳播速度c、氣體介質的絕對溫度T以及由氣體介質組成成分所決定的聲音常數B之間的關系[2]可以表示為

(1)

采用聲學CT技術測量溫度場需要在被測區域的邊界處設置一定數量的聲波收發器.各聲波收發器依次發聲，所有的收發器都可以接收該信號，從而在被測區域內部形成一定數量的聲波傳播路徑.測量出聲波在各聲波傳播路徑上的飛行時間，利用某種重建算法以及聲速與溫度的對應關系，即可計算重建出被測區域內的溫度分布.

2 重建BFAKF算法

聲學CT重建屬逆問題研究，通常要借助正問題模型來求解.本文提出的基于基函數逼近和卡爾曼濾波的重建算法(BFAKF)可分為正問題建模和逆問題求解兩個環節.

2.1 基于基函數逼近的正問題建模

聲學CT的正問題是在已知聲波收發器的位置并確定有效聲線后，根據溫度(聲速)分布，求解各路徑上的聲波飛行時間.

假設被測區域內聲速倒數的分布可以表示為f(x，y，z)，則聲波在第k條傳播路徑lk上的傳播時間tk可以表示為

(2)

式中，N為有效聲波路徑數.將被測區域均勻地劃分為M個網格，且M>N，利用M個Markov徑向基函數的線性組合來表示f(x，y，z)，則有

(3)

式中：βi為待定系數；xi、yi、zi為第i個網格的中心點坐標；a為徑向基函數的形狀參數.將式(3)代入式(2)得

(4)

令A=(aki)k=1，2，…，N；i=1，2，…，M，t=[t1，t2，…，tN]T，β=[β1，β2，…，βM]T，可以得到用基函數逼近方法建立的聲學CT溫度場重建正問題模型，即

t=Aβ

(5)

2.2 基于卡爾曼濾波的逆問題求解

卡爾曼濾波[7-8]是一種遞歸的估計算法，在已知系統模型、噪聲統計信息和初始狀態的前提下，利用上一時刻的狀態估計值和當前狀態的觀測值實現對當前時刻狀態值的估計.假設考慮噪聲的系統離散模型為

(6)

式中：Xk為k時刻的被估計狀態；Φk/k-1為k-1時刻到k時刻狀態轉移矩陣；Ψk為系統噪聲驅動矩陣；Wk為系統噪聲；Hk為狀態觀測值；Zk為觀測矩陣；Ok為觀測噪聲.

針對式(6)的卡爾曼濾波方程可表示為

(7)

式中：Kk為濾波增益矩陣；Rk為觀測噪聲協方差矩陣.卡爾曼濾波正是通過不斷調整對估計和觀測的依賴程度，快速而有效地獲得被估計量的最優估計值.

將卡爾曼濾波用于聲學CT溫度場重建逆問題求解時先要建立系統模型.假設測量聲波飛行期間被測區域的溫度(聲速)不發生變化，即認為狀態轉移矩陣Φk/k-1為單位陣，系統噪聲為0，則由k-1時刻到k時刻系統狀態預測方程為

βk=βk-1

(8)

由于聲波飛行時間測量過程中觀測噪聲的方差通常是比較穩定的，故假設觀測噪聲協方差矩陣Rk=R不隨時間而變，則式(5)所對應的系統狀態觀測方程為

tk=Aβk+Ok

(9)

相應卡爾曼濾波方程可寫為

(10)

3 重建誤差評價指標

本文采用重建溫度場的平均相對誤差Rave、均方根誤差Rrms以及熱點溫度相對誤差Rh來評價重建質量，其定義分別為

(11)

(12)

(13)

4 重建算法仿真驗證

本文將32個聲波收發器均勻地布置在12 m×12 m×12 m的立方體區域周圍，形成172條有效聲波路徑.采用最小二乘法(LSM)、Markov徑向基函數Tikhonov正則化法(MTR)以及本文提出的基函數逼近卡爾曼濾波法(BFAKF)對四種模型溫度場進行了仿真重建，各溫度場表達式為

1) 熱點位于(0，0，0)的單峰對稱場

(14)

2) 熱點位于(-2，-2，0)的單峰偏置場

(15)

3) 熱點位于(-2，0，-2)和(2，0，2)的雙峰對稱場

(16)

4) 熱點位于(-2，0，2)、(2，0，2)、(-2，0，-2)和(2，0，-2)的四峰對稱場

(17)

采用LSM法重建時，被測區域均勻地劃分為4×4×4=64個網格，即原始像素總數為64；采用MTR和BFAKF法重建時，被測區域均勻地劃分為10×10×10=1 000個網格，即原始像素總數為1 000.MTR法的正則化參數μ=1×10-6，形狀參數a=1×10-4；BFAKF法的形狀參數a=1×10-4，估計誤差協方差矩陣初值P0=5I，迭代次數為50，飛行時間數據無噪聲時取觀測噪聲協方差矩陣R0=1×10-5I，飛行時間數據有噪聲(噪聲標準差δ=2×10-5)時取R0=7×10-5I.

圖1以單峰對稱和四峰對稱模型溫度場為例，給出了有噪聲情況下BFAKF法的均方根誤差、平均相對誤差隨迭代次數的變化曲線.表1、2分別給出了采用無噪聲飛行時間和有噪聲飛行時間時LSM法、MTR法以及BFAKF法對應的重建誤差.由于LSM法無法正確地重建出四峰對稱溫度場特征，故本文不給出四峰對稱溫度場的LSM法熱點溫度相對誤差.算法的重建時間分別是：LSM算法0.023 61 s、MTR算法0.116 46 s、BFAKF算法2.157 2 s.LSM法只有64個原始像素，正問題模型中的系數矩陣A只有172×64個元素，MTR法有1 000個原始像素，矩陣A有172×1 000個元素，所以逆問題求解需要更長的時間.

圖1 誤差隨迭代次數變化曲線Fig.1 Curves of errors varying with iteration number

表1 采用無噪聲飛行時間和原始像素時的重建誤差Tab.1 Reconstruction errors with noise-free time-of-flight and original pixels %

表2 采用有噪聲飛行時間和原始像素時的重建誤差Tab.2 Reconstruction errors with noisy time-of-flight and original pixels %

表3給出了飛行時間數據無噪聲時LSM、MTR以及BFAKF法對應的重建溫度場與模型溫度場.為用盡可能少的篇幅顯示溫度場特征，單峰溫度場只給出了x=0 m、y=0 m、z=0 m對應的切片圖；雙峰和四峰溫度場只給出了y=0 m處切片圖.由于僅用64或者1 000個像素難以細致地描述一個復雜的溫度場，故而表4、5給出了細化像素對應的重建誤差.具體做法是先將重建出的溫度值賦予原始像素的中心點，然后將12 m×12 m×12 m的區域細化成31×31×31=29 791個像素，再用三次樣條插值的方法求出細化像素中心點的溫度值.LSM法的原始像素數M=4×4×4=64，細化像素數M′=23×23×23=12 167，細化像素的描述區域為9 m×9 m×9 m；MTR和BFAKF法的原始像素數M=10×10×10=1 000，細化像素數M′=27×27×27=19 683，細化像素的描述區域為10.8 m×10.8 m×10.8 m.

表3 模型溫度場與重建溫度場對比Tab.3 Comparison between model and reconstruction temperature fields

5 結 論

本文提出一種利用基函數逼近和卡爾曼濾波的聲學CT溫度場重建算法.采用BFAKF算法以及有代表性的LSM、MTR算法分別對四種典型的溫度場模型進行了仿真數據重建，重建結果表明：BFAKF法能使重建誤差以較快的速度下降并趨于穩定；與LSM、MTR法相比，BFAKF法對應的重建誤差更小，重建溫度場與模型溫度場更接近，BFAKF算法具有更好的復雜溫度場重建能力.但BFAKF算法的重建時間較長.一方面是因為該算法采用了較多的原始像素導致正問題模型中的系數矩陣維數較大；另一方面由于該算法是迭代算法，系數矩陣維數大必然會導致逆問題求解時間長.

表4 采用無噪聲飛行時間和細化像素時的重建誤差Tab.4 Reconstruction errors with noise-free time-of-flight and refined pixels %

表5 采用有噪聲飛行時間和細化像素時的重建誤差Tab.5 Reconstruction errors with noisy time-of-flight and refined pixels %