魏元焜,顏 華,周英鋼

(沈陽工業(yè)大學(xué)信息科學(xué)與工程學(xué)院,遼寧沈陽 110870)

0 引言

聲學(xué)CT溫度場重建技術(shù)[1]具有非接觸、不干擾被測溫度場、測量范圍廣、環(huán)境適應(yīng)能力強(qiáng)、維護(hù)方便等優(yōu)點(diǎn),作為一種新型工業(yè)測溫儀表的關(guān)鍵技術(shù),具有廣泛的應(yīng)用前景。該技術(shù)可用于聲學(xué)高溫計(jì),用來監(jiān)測鍋爐爐膛、煙道燃燒室的溫度分布[1-2]。

以聲學(xué)CT重建溫度場時(shí),需要在被測區(qū)域周圍設(shè)置多個(gè)聲波收發(fā)器,形成M條從不同方向穿越被測區(qū)域的聲波路徑(簡稱聲線),并測得聲波在各聲線上的飛行時(shí)間(time of flight,TOF)。根據(jù)溫度分布求TOF是聲學(xué)CT溫度場重建的正問題,由測得的TOF反推溫度分布是其逆問題。為求解逆問題,通常將測溫區(qū)域劃分成N個(gè)網(wǎng)格,利用正問題建立數(shù)學(xué)模型,然后利用重建算法根據(jù)TOF反推N個(gè)網(wǎng)格內(nèi)的聲慢(聲速倒數(shù))值,最后以插值算法獲得細(xì)密的溫度分布。經(jīng)典的重建算法有最小二乘法(least square method,LSA)[3-4]、聯(lián)合迭代重建技術(shù)(simultaneous iterative reconstruction technique,SIRT)[5-6]等。實(shí)際應(yīng)用中可獲得的TOF數(shù)據(jù)非常有限,因此N大于M的重建算法更適合于復(fù)雜聲慢度(溫度)分布的重建。但N大于M時(shí)逆問題是欠定的(即線性方程組的方程數(shù)少于未知數(shù)個(gè)數(shù)),且通常N值越大,聲學(xué)CT投影矩陣條件數(shù)越大,即重建結(jié)果對TOF的測量誤差越敏感。為解決這一問題,取得穩(wěn)定的逆問題解,經(jīng)典算法往往要求N

為提升聲學(xué)CT溫度場重建精度,應(yīng)對N>M導(dǎo)致的逆問題求解的不穩(wěn)定,本文提出了不同的思路,提出了基于稀疏重構(gòu)(sparse reconstruction)的聲學(xué)CT溫度場重建算法,簡稱SR算法。

1 SR算法

1.1 數(shù)學(xué)模型的建立

以二維聲學(xué)CT溫度場重建為例考慮聲學(xué)CT正問題,設(shè)測溫平面內(nèi)聲慢度分布函數(shù)為f(x,y),則被測區(qū)域內(nèi)第j(1

(1)

式中:fi為網(wǎng)格i內(nèi)的聲慢度值;aji為聲線pj在網(wǎng)格i內(nèi)的長度。

同時(shí)考慮M條聲線,定義t=(t1,…,tM)T,A=(aji),f=(f1,…,fN)T,j=1,…,M,i=1,…,N,則聲學(xué)CT系統(tǒng)獲得TOF的過程可表示為

t=Af

(2)

式中:f為N×1維向量,由N個(gè)網(wǎng)格內(nèi)的聲慢值構(gòu)成,可看作系統(tǒng)的觀測對象;A為M×N維矩陣,為描述系統(tǒng)的觀測矩陣,表征了聲學(xué)CT系統(tǒng)對聲慢度的觀測過程。

分析式(1)、式(2)可知,矩陣A的元素aji可視為聲線pj對網(wǎng)格i內(nèi)聲慢度信號(hào)的采樣;M×1維向量t為系統(tǒng)的觀測結(jié)果,由M條聲線上的聲波TOF構(gòu)成,可由聲學(xué)CT系統(tǒng)測得。收發(fā)器位置固定后,矩陣A為已知量。因此在式(2)所示模型下,聲學(xué)CT逆問題的實(shí)質(zhì)是根據(jù)已知的t和A求f。

由壓縮感知理論[15-17]可知,只要信號(hào)在某個(gè)變換域是稀疏的,那么就可以由一個(gè)與變換基不相干的觀測矩陣所獲得的少量觀測數(shù)據(jù),然后利用稀疏重構(gòu)算法極高概率地重構(gòu)出原信號(hào)。這一過程與聲學(xué)CT溫度場重建相似,可作為N>M的逆問題求解思路。但聲速不可能為無窮大,所以聲慢度不會(huì)是0,即聲慢度向量f不是稀疏的,因此需要按式(3)將f變換到稀疏域:

f=Ψθ

(3)

式中:Ψ為N×N維矩陣,通常Ψ是一個(gè)標(biāo)準(zhǔn)正交基;θ為N×1維向量,且θ中大部分元素為0或非常接近0。

稱Ψ為f的稀疏基,θ為f在Ψ域的稀疏表示,θ中非零元素的個(gè)數(shù)為其稀疏度。常用的稀疏基有離散傅里葉變換(DFT)基和離散余弦變換(DCT)基,本文研究發(fā)現(xiàn)聲慢信號(hào)在DCT域中具有更好的稀疏性,本文以DFT作為SR算法的稀疏基Ψ。

基于f在DFT域的稀疏性,將式(3)帶入式(2),即可得到基于稀疏重構(gòu)的聲學(xué)CT數(shù)學(xué)模型:

t=AΨθ=Uθ

(4)

式中:θ為邏輯上的被測稀疏信號(hào);U為等效觀測矩陣,U=AΨ。

(5)

式中z為由氣體組分決定的常數(shù),空氣中z=20.05。

最后利用三次樣條插值算法,即可獲得更加細(xì)致像素下的溫度場描述。

在稀疏重構(gòu)中有2個(gè)重要環(huán)節(jié)——觀測矩陣的優(yōu)化和稀疏重構(gòu),將在下文展開論述。

1.2 觀測矩陣的優(yōu)化

由壓縮感知理論可知,為實(shí)現(xiàn)稀疏信號(hào)θ的高精度重構(gòu),觀測矩陣A和稀疏基Ψ的相干性要盡可能低。高斯隨機(jī)矩陣因其隨機(jī)性,在很大概率上與任意稀疏基不相干,因此壓縮感知中常采用獨(dú)立同分布的高斯隨機(jī)矩陣作為觀測矩陣[18-19]。但聲學(xué)CT系統(tǒng)的結(jié)構(gòu)化特質(zhì)使其無法實(shí)現(xiàn)高斯隨機(jī)矩陣,因此需要對矩陣A進(jìn)行優(yōu)化,使其盡量接近高斯隨機(jī)矩陣的隨機(jī)特性。本文通過提高聲線穿過各個(gè)網(wǎng)格的“無差別化”和“無偏頗性”,使聲線在各網(wǎng)格內(nèi)的分布更充分、更均勻,從而等效提高觀測矩陣A對聲慢度信號(hào)采樣的隨機(jī)性。

為定量評價(jià)聲線在各網(wǎng)格內(nèi)布局的充分、均勻程度,定義N×1維向量Γ如式(6)所示:

(6)

式中:Γi為網(wǎng)格i內(nèi)所有聲線長度之和,i=1,2,…,N;當(dāng)測溫區(qū)域和網(wǎng)格劃分確定后,aji為各收發(fā)器位置坐標(biāo)構(gòu)成的向量X的函數(shù)。

定義聲線布局充分均勻度F(X)如式(7)所示,

F(X)=μF(Γ)-σF(Γ)

(7)

式中:X為各收發(fā)器位置坐標(biāo)構(gòu)成的向量;μF(Γ)為向量Γ各元素的均值,反映聲線布局的充分性,其數(shù)值越大充分性越好;σF(Γ)為向量Γ各元素的標(biāo)準(zhǔn)差,反映聲線布局的均勻性,其數(shù)值越小均勻性越好。

F(X)數(shù)值越大,聲線布局充分均勻性越好,F(X)數(shù)值越小,聲線充分均勻性越差。在給定被測區(qū)域及N值和M值的情況下,F(X)僅由收發(fā)器位置坐標(biāo)向量X決定。

本文以求F(X)最大值為優(yōu)化目標(biāo),用遺傳算法求解X。種群規(guī)模設(shè)置為200,最大代數(shù)為300,遷移率為0.2,雜交率為0.8。設(shè)置被測區(qū)域?yàn)檫呴L1.3 m的正方形,劃分N=15×15=225個(gè)網(wǎng)格。4個(gè)邊框上各布置3個(gè)聲波收發(fā)器,形成M=54條有效穿越被測區(qū)域的聲線。由于收發(fā)器固定在正方形的邊框上,因此每個(gè)收發(fā)器必有一個(gè)坐標(biāo)是常數(shù),要優(yōu)化的收發(fā)器位置坐標(biāo)X可縮減為12×1維向量。

圖1(a)給出了常用的收發(fā)器標(biāo)準(zhǔn)布局[3,4,10,13]以及對應(yīng)的聲線。該布局也被用作遺傳算法的初值。圖1(b)給出了用遺傳算法獲得的優(yōu)化布局。優(yōu)化前、后的聲線布局充分均勻度F(X)分別為0.122 7和0.138 7。相對于優(yōu)化前,優(yōu)化后的充分均勻度提高11.5%。

圖1 優(yōu)化前后的收發(fā)器布局

1.3 稀疏信號(hào)的重構(gòu)

根據(jù)式(4)由已知的t和U求θ的過程可歸結(jié)為式(8)所示的l0范數(shù)最小化問題[19]。

(8)

式(8)可由OMP算法[15-16]求解,具體流程如下:

(1)初始化,令殘差r0=t,索引向量Λ0為空,支撐矩陣V0為空,迭代計(jì)數(shù)m=1。

(2)列選擇,計(jì)算矩陣U各列與殘差rm-1的內(nèi)積。將內(nèi)積最大的列向量加入支撐矩陣Vm-1,構(gòu)成新的支撐集矩陣Vm;將選中的列向量在矩陣U中的列索引號(hào)加入索引向量Λm-1,構(gòu)成新的索引向量Λm;最后,以零向量替換U中被選中列。

(3)更新殘差,以最小二乘法求解式(9),計(jì)算θm并以式(10)更新殘差。

(9)

rm=t-Vmθm

(10)

(4)確定是否滿足迭代終止條件,若滿足,則停止迭代并跳轉(zhuǎn)到步驟(5);否則令m=m+1,返回步驟(2)。

2 溫度場重建實(shí)驗(yàn)

在自行設(shè)計(jì)的聲學(xué)CT溫度場重建裝置上實(shí)施了重建實(shí)驗(yàn),被測區(qū)域?yàn)檫呴L1.3 m的正方形。實(shí)驗(yàn)裝置的系統(tǒng)框圖如圖2所示,圖2中關(guān)鍵模塊實(shí)物如圖3所示。計(jì)算機(jī)控制12個(gè)揚(yáng)聲器(TC9FD-19-08)輪流發(fā)聲,12個(gè)傳聲器(CHZ 221)將接收到的聲波信號(hào)轉(zhuǎn)換成電信號(hào),經(jīng)信號(hào)調(diào)理器后送入數(shù)據(jù)采集器(PCI 6123),利用互相關(guān)時(shí)延估計(jì)算法可獲取各聲線上的TOF,重建程序調(diào)用重建算法重建測溫層面上的溫度分布。

圖2 硬件系統(tǒng)框圖

圖3 硬件實(shí)物圖

實(shí)驗(yàn)選取了4種聲學(xué)CT重建算法:本文提出的SR算法,其觀測矩陣對應(yīng)圖1(b)所示的優(yōu)化布局,記為SR-b;采用了標(biāo)準(zhǔn)收發(fā)器布局的SR算法,其觀測矩陣對應(yīng)圖1(a)所示的標(biāo)準(zhǔn)布局,記為SR-a算法;采用了標(biāo)準(zhǔn)收發(fā)器布局的LSA算法以及SIRT算法,其觀測矩陣對應(yīng)圖1(a)所示的標(biāo)準(zhǔn)布局。SR算法采用N=15×15的網(wǎng)格劃分,LSA和SIRT為避免逆問題求解的不穩(wěn)定,采用了N=5×5的網(wǎng)格劃分。SIRT和LSA算法按文獻(xiàn)[5,7]實(shí)現(xiàn),其中SIRT算法的松弛因子為0.3,當(dāng)相鄰2次迭代結(jié)果之差相對上次迭代結(jié)果之比小于1×10-4時(shí),停止SIRT迭代。SR算法估計(jì)σ參數(shù)的方法如下:數(shù)據(jù)采集卡PCI 6123的最大采樣率為500 kHz,由此可先確定采樣間隔所對應(yīng)的最大量化誤差σs=2 μs;然后在均溫場TF1下連續(xù)進(jìn)行了40次TOF測量,計(jì)算出所有聲線TOF標(biāo)準(zhǔn)差的平均值σm=1.1 μs;最后取σs和σm中的數(shù)值較大者作為

TOF所含噪聲標(biāo)準(zhǔn)差σ的估計(jì)值,即σ=2×10-6。

以4種算法分別重建301.2 K的均溫場(TF1)和單峰對稱溫度場(TF2)。為保證標(biāo)準(zhǔn)布局和優(yōu)化布局下TF2的近似性,電熱爐的位置保持不變。利用空調(diào)調(diào)節(jié)室內(nèi)溫度,在保證加熱前測溫區(qū)溫度均為301.2 K的情況下,開啟電熱爐。圖4給出了兩種收發(fā)器布局下的TF2溫度場。

(a)標(biāo)準(zhǔn)布局

(b)優(yōu)化布局圖4 2種收發(fā)器布局下的單峰溫度場TF2

圖5給出了4種算法對均溫場TF1、單峰場TF2的重建圖像。由圖5可以看出:4種算法的重建結(jié)果均能反映實(shí)際溫度場的主要特征;對于TF1而言,SR-b和SR-a算法重建的圖像呈現(xiàn)出非常好的均勻性,而SIRT和LSA算法重建的圖像有明顯的非均勻性;對于TF2,SR-b算法重建圖像中偽影的顏色與背景的顏色最接近,而其他3種算法在被測區(qū)域的四角或邊界附近偽影的顏色與背景顏色差異明顯。

(a)SR-b重建TF1

(b)SR-a重建TF1

(c)LSA重建TF1

(d)SIRT重建TF1

(e)SR-b重建TF2

(f)SR-a重建TF2

(g)LSA重建TF2

為對重建溫度場進(jìn)行誤差評價(jià),通常以熱電偶進(jìn)行接觸式測量。由于在測試平面布置大量的熱電偶會(huì)影響聲學(xué)CT的精確重建,因此重建誤差評價(jià)通常只針對關(guān)鍵點(diǎn)進(jìn)行[10-13]。本文將80PK-1型熱電偶安裝在高精度手持式數(shù)字溫度計(jì)(Fluke 54 II)上實(shí)現(xiàn)溫度的實(shí)時(shí)顯示。對于TF1溫度場,用熱電偶測量了被測層面的最高溫度Tmax和最低溫度Tmin。對于TF2溫度場,用熱電偶測量了被測層面上電熱爐上方的最高溫度。表1給出了熱電偶測量和算法重建溫度的對比。重建溫度相對于熱電偶測量值的誤差e也在表1中給出。為便于整體評價(jià),表1中的最后一行還給出了同一算法所對應(yīng)的|e|的平均值。由表1可以看出:SR-b算法的重建溫度與熱電偶的測量溫度最為接近,SR-a算法與LSA算法重建溫度與熱電偶測溫偏差逐漸增大;從|e|的平均值來看,SR-b的重建誤差平均值最小,SR-b、SR-a的重建誤差明顯低于SIRT、LSA算法,SR-b和SR-a的重建誤差低于LSA誤差最為明顯。

表1 熱電偶測量與算法重建溫度對比

溫度場TF2是在301.2 K的均溫場基礎(chǔ)上通過電熱爐局部加熱形成的,因此TF2溫度場中所有像素(三次樣條插值的細(xì)密像素37 pixel×37 pixel)的真實(shí)溫度都不應(yīng)出現(xiàn)低于301.2 K的異常結(jié)果。但受重建誤差的影響,4種算法的TF2重建圖像中都存在溫度值低于301.2K的“異常”像素。顯然,異常像素?cái)?shù)目越少,異常像素的溫度越接近301.2 K,重建溫度場越接近實(shí)際溫度場,則算法性能越好。表2給出了4種算法所對應(yīng)的異常像素的數(shù)目Nabn、異常像素溫度的平均值Tavgabn和最小值Tminabn。由表2可以看出,4種算法中SR-b算法的Nabn最小,Tavgabn和Tminabn與301.2 K也最接近。

表2 異常像素統(tǒng)計(jì)

綜合圖5、表1和表2可知,4種算法中SR-b算法的重建質(zhì)量最優(yōu),LSA算法的重建質(zhì)量最差,SR-a算法的重建質(zhì)量優(yōu)于SIRT、LSA算法。因此從重建誤差來看,本文提出的稀疏重構(gòu)算法優(yōu)于經(jīng)典的LSA和SIRT算法,且本文的觀測矩陣優(yōu)化方法可以有效降低重建誤差。

3 結(jié)束語

本文提出的SR算法利用聲慢度信號(hào)在DFT域的稀疏性,在稀疏重構(gòu)的框架中求解聲學(xué)CT逆問題,允許N>M的聲學(xué)CT溫度場重建。通過溫度場重建實(shí)驗(yàn)可以得出以下結(jié)論:

(1)SR算法重建溫度場的誤差低于LSA算法和SIRT算法,證實(shí)了稀疏重構(gòu)算法在溫度場重建應(yīng)用中的可行性與優(yōu)越性。

(2)SR-b算法重建誤差低于SR-a算法,說明通過優(yōu)化收發(fā)器布局,可以有效提高聲學(xué)CT采樣隨機(jī)性,通過滿足壓縮感知理論中對采樣隨機(jī)性的要求,可以有效提高基于稀疏重構(gòu)的聲學(xué)CT溫度場重建質(zhì)量。

(3)SR算法中采用的噪聲水平相關(guān)的OMP迭代終止條件在溫度場重建實(shí)驗(yàn)中,表現(xiàn)出了良好的穩(wěn)定性,表明SR算法在室內(nèi)實(shí)際應(yīng)用中可以穩(wěn)定地重建溫度場。