基于基函數(shù)逼近和卡爾曼濾波的溫度場(chǎng)重建*

顏 華，顧夢(mèng)楠，王伊凡

(沈陽(yáng)工業(yè)大學(xué) 信息科學(xué)與工程學(xué)院，沈陽(yáng) 110870)

基于聲學(xué)CT的溫度場(chǎng)重建技術(shù)[1]是根據(jù)被測(cè)區(qū)域內(nèi)多條路徑上的聲波飛行時(shí)間數(shù)據(jù)推算被測(cè)區(qū)域的溫度分布，屬于“由效果反求原因”的逆問(wèn)題研究.該技術(shù)具有非接觸不干擾被測(cè)溫度場(chǎng)、測(cè)溫范圍廣、環(huán)境適應(yīng)能力強(qiáng)及可在線(xiàn)測(cè)量等優(yōu)點(diǎn).

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

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

1 聲學(xué)CT溫度場(chǎng)重建原理

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

(1)

采用聲學(xué)CT技術(shù)測(cè)量溫度場(chǎng)需要在被測(cè)區(qū)域的邊界處設(shè)置一定數(shù)量的聲波收發(fā)器.各聲波收發(fā)器依次發(fā)聲，所有的收發(fā)器都可以接收該信號(hào)，從而在被測(cè)區(qū)域內(nèi)部形成一定數(shù)量的聲波傳播路徑.測(cè)量出聲波在各聲波傳播路徑上的飛行時(shí)間，利用某種重建算法以及聲速與溫度的對(duì)應(yīng)關(guān)系，即可計(jì)算重建出被測(cè)區(qū)域內(nèi)的溫度分布.

2 重建BFAKF算法

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

2.1 基于基函數(shù)逼近的正問(wèn)題建模

聲學(xué)CT的正問(wèn)題是在已知聲波收發(fā)器的位置并確定有效聲線(xiàn)后，根據(jù)溫度(聲速)分布，求解各路徑上的聲波飛行時(shí)間.

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

(2)

式中，N為有效聲波路徑數(shù).將被測(cè)區(qū)域均勻地劃分為M個(gè)網(wǎng)格，且M>N，利用M個(gè)Markov徑向基函數(shù)的線(xiàn)性組合來(lái)表示f(x，y，z)，則有

(3)

式中：βi為待定系數(shù)；xi、yi、zi為第i個(gè)網(wǎng)格的中心點(diǎn)坐標(biāo)；a為徑向基函數(shù)的形狀參數(shù).將式(3)代入式(2)得

(4)

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

t=Aβ

(5)

2.2 基于卡爾曼濾波的逆問(wèn)題求解

卡爾曼濾波[7-8]是一種遞歸的估計(jì)算法，在已知系統(tǒng)模型、噪聲統(tǒng)計(jì)信息和初始狀態(tài)的前提下，利用上一時(shí)刻的狀態(tài)估計(jì)值和當(dāng)前狀態(tài)的觀(guān)測(cè)值實(shí)現(xiàn)對(duì)當(dāng)前時(shí)刻狀態(tài)值的估計(jì).假設(shè)考慮噪聲的系統(tǒng)離散模型為

(6)

式中：Xk為k時(shí)刻的被估計(jì)狀態(tài)；Φk/k-1為k-1時(shí)刻到k時(shí)刻狀態(tài)轉(zhuǎn)移矩陣；Ψk為系統(tǒng)噪聲驅(qū)動(dòng)矩陣；Wk為系統(tǒng)噪聲；Hk為狀態(tài)觀(guān)測(cè)值；Zk為觀(guān)測(cè)矩陣；Ok為觀(guān)測(cè)噪聲.

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

(7)

式中：Kk為濾波增益矩陣；Rk為觀(guān)測(cè)噪聲協(xié)方差矩陣.卡爾曼濾波正是通過(guò)不斷調(diào)整對(duì)估計(jì)和觀(guān)測(cè)的依賴(lài)程度，快速而有效地獲得被估計(jì)量的最優(yōu)估計(jì)值.

將卡爾曼濾波用于聲學(xué)CT溫度場(chǎng)重建逆問(wèn)題求解時(shí)先要建立系統(tǒng)模型.假設(shè)測(cè)量聲波飛行期間被測(cè)區(qū)域的溫度(聲速)不發(fā)生變化，即認(rèn)為狀態(tài)轉(zhuǎn)移矩陣Φk/k-1為單位陣，系統(tǒng)噪聲為0，則由k-1時(shí)刻到k時(shí)刻系統(tǒng)狀態(tài)預(yù)測(cè)方程為

βk=βk-1

(8)

由于聲波飛行時(shí)間測(cè)量過(guò)程中觀(guān)測(cè)噪聲的方差通常是比較穩(wěn)定的，故假設(shè)觀(guān)測(cè)噪聲協(xié)方差矩陣Rk=R不隨時(shí)間而變，則式(5)所對(duì)應(yīng)的系統(tǒng)狀態(tài)觀(guān)測(cè)方程為

tk=Aβk+Ok

(9)

相應(yīng)卡爾曼濾波方程可寫(xiě)為

(10)

3 重建誤差評(píng)價(jià)指標(biāo)

本文采用重建溫度場(chǎng)的平均相對(duì)誤差Rave、均方根誤差Rrms以及熱點(diǎn)溫度相對(duì)誤差Rh來(lái)評(píng)價(jià)重建質(zhì)量，其定義分別為

(11)

(12)

(13)

4 重建算法仿真驗(yàn)證

本文將32個(gè)聲波收發(fā)器均勻地布置在12 m×12 m×12 m的立方體區(qū)域周?chē)?，形?72條有效聲波路徑.采用最小二乘法(LSM)、Markov徑向基函數(shù)Tikhonov正則化法(MTR)以及本文提出的基函數(shù)逼近卡爾曼濾波法(BFAKF)對(duì)四種模型溫度場(chǎng)進(jìn)行了仿真重建，各溫度場(chǎng)表達(dá)式為

1) 熱點(diǎn)位于(0，0，0)的單峰對(duì)稱(chēng)場(chǎng)

(14)

2) 熱點(diǎn)位于(-2，-2，0)的單峰偏置場(chǎng)

(15)

3) 熱點(diǎn)位于(-2，0，-2)和(2，0，2)的雙峰對(duì)稱(chēng)場(chǎng)

(16)

4) 熱點(diǎn)位于(-2，0，2)、(2，0，2)、(-2，0，-2)和(2，0，-2)的四峰對(duì)稱(chēng)場(chǎng)

(17)

采用LSM法重建時(shí)，被測(cè)區(qū)域均勻地劃分為4×4×4=64個(gè)網(wǎng)格，即原始像素總數(shù)為64；采用MTR和BFAKF法重建時(shí)，被測(cè)區(qū)域均勻地劃分為10×10×10=1 000個(gè)網(wǎng)格，即原始像素總數(shù)為1 000.MTR法的正則化參數(shù)μ=1×10-6，形狀參數(shù)a=1×10-4；BFAKF法的形狀參數(shù)a=1×10-4，估計(jì)誤差協(xié)方差矩陣初值P0=5I，迭代次數(shù)為50，飛行時(shí)間數(shù)據(jù)無(wú)噪聲時(shí)取觀(guān)測(cè)噪聲協(xié)方差矩陣R0=1×10-5I，飛行時(shí)間數(shù)據(jù)有噪聲(噪聲標(biāo)準(zhǔn)差δ=2×10-5)時(shí)取R0=7×10-5I.

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

圖1 誤差隨迭代次數(shù)變化曲線(xiàn)Fig.1 Curves of errors varying with iteration number

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

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

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

表3 模型溫度場(chǎng)與重建溫度場(chǎng)對(duì)比Tab.3 Comparison between model and reconstruction temperature fields

5 結(jié) 論

本文提出一種利用基函數(shù)逼近和卡爾曼濾波的聲學(xué)CT溫度場(chǎng)重建算法.采用BFAKF算法以及有代表性的LSM、MTR算法分別對(duì)四種典型的溫度場(chǎng)模型進(jìn)行了仿真數(shù)據(jù)重建，重建結(jié)果表明：BFAKF法能使重建誤差以較快的速度下降并趨于穩(wěn)定；與LSM、MTR法相比，BFAKF法對(duì)應(yīng)的重建誤差更小，重建溫度場(chǎng)與模型溫度場(chǎng)更接近，BFAKF算法具有更好的復(fù)雜溫度場(chǎng)重建能力.但BFAKF算法的重建時(shí)間較長(zhǎng).一方面是因?yàn)樵撍惴ú捎昧溯^多的原始像素導(dǎo)致正問(wèn)題模型中的系數(shù)矩陣維數(shù)較大；另一方面由于該算法是迭代算法，系數(shù)矩陣維數(shù)大必然會(huì)導(dǎo)致逆問(wèn)題求解時(shí)間長(zhǎng).

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

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