基于PCA 降維和迭代正則化的溫度場重建

白 云，顏 華，魏元焜

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

聲學(xué)CT 溫度場重建技術(shù)[1-6]具有非接觸、環(huán)境適應(yīng)力強(qiáng)、維護(hù)方便、實時監(jiān)測等優(yōu)點[7-8]，備受工業(yè)界青睞。重建系統(tǒng)的性能很大程度上依賴于重建算法的精度。依正問題建模方式的不同重建算法主要有兩大類。第一類基于網(wǎng)格內(nèi)聲線長度建模[9]，典型算法有最小二乘法（簡稱LSM 法）[10]，網(wǎng)格數(shù)N 通常小于聲線數(shù)M。針對邊緣信息缺失等問題，文獻(xiàn)[11]利用均衡優(yōu)化算法得到粗網(wǎng)格下的溫度分布，再利用高斯過程回歸算法預(yù)測細(xì)密網(wǎng)格下的溫度分布。第二類采用徑向基函數(shù)逼近聲慢度分布的方式建模，允許N>M，用正則化等手段求解嚴(yán)重病態(tài)的逆問題[12]。例如文獻(xiàn)[13]采用截斷的廣義奇異值分解法。

重建區(qū)域劃分的網(wǎng)格數(shù)越多，對復(fù)雜溫度場的描述能力越強(qiáng)，但重建矩陣的奇異性顯著增加。針對此問題本文提出一種基于主成分分析（principal component analysis，PCA）降維和迭代正則化的重建算法。該算法采用第二類建模方式，用PCA 降維改善逆問題的病態(tài)性，用迭代正則化法求解逆問題，進(jìn)而得到溫度分布。

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

聲波在氣體介質(zhì)中的傳播速度c 與氣體介質(zhì)的絕對溫度間T 間的關(guān)系為

式中：z 是聲音常數(shù)，由氣體組成成分決定，對空氣而言有z=20.045。

用聲學(xué)CT 法重建溫度場需事先將多組聲波收發(fā)器布置在待測區(qū)域周圍。系統(tǒng)工作時依據(jù)收發(fā)器采集到的聲波數(shù)據(jù)，計算聲波在收發(fā)器間的飛行時間，通過預(yù)設(shè)的重建算法，計算被測層面的聲慢度（聲速的倒數(shù)）分布，進(jìn)而利用式（1）得到溫度分布。

2 基于PCA 降維和迭代正則化的重建算法

2.1 正問題建模

利用布置在被測層面周圍的收發(fā)器在被測層面形成M 條聲線，并將被測區(qū)域離散成N 個網(wǎng)格。設(shè)被測區(qū)域內(nèi)的聲慢度分布為s（x，y），第i 個網(wǎng)格中心點的坐標(biāo)為（xi，yi），則可用式（2）計算聲波在第k 條聲線lk上的飛行時間：

本文利用N 個徑基函數(shù)的線性組合逼近聲慢度分布，如式（3）所示：

式中：βi為待定系數(shù)；α＞0 為徑向基函數(shù)的形狀參數(shù)，根據(jù)所需測量范圍以及收發(fā)器位置擺放，利用數(shù)值實驗確定此參數(shù)。令：

將式（4）代入式（2），有：

考慮M 條路徑并令A(yù)=（Aki）k=1，…，M，i=1，…，N，t=（t1，…，tM）T，β=（β1，…，βN）T，建立正問題模型如式（6）所示：

式中：β為N 維待定系數(shù)向量；A為M×N 維重建矩陣；t為M 維聲波飛行時間向量。本文M=54，N=81。

2.2 逆問題求解

2.2.1 基于主成分分析的降維處理

為獲得良好的重建質(zhì)量，需要獲得盡可能多的反映溫（聲慢）度分布信息的聲波飛行時間數(shù)據(jù)。相應(yīng)的，重建矩陣也成為典型的高維度數(shù)據(jù)，有很強(qiáng)的行相關(guān)性。為減輕冗余特征帶來的多重共線性問題，本文通過對重建矩陣A的主成分分析，對正問題模型進(jìn)行降維處理。

主成分分析[14-16]通過對原始的M 維變量的正交變換，得到被稱為主成分的線性無關(guān)的新變量。第一個主成分具有最大的信息量（方差值），后續(xù)主成分包含的信息量依次降低。前K（K

（1）將矩陣A的每一個行向量零均值化處理，得到矩陣Y，而后求矩陣Y的協(xié)方差矩陣C：

（2）先求矩陣C的特征值和特征向量，而后由上到下排列前K 大的特征值所對應(yīng)的特征向量，形成K×N 維轉(zhuǎn)換矩陣P；

（3）令A(yù)dr=PA得到降維后的K×N 維重建矩陣Adr。

本文K 取47，P為47×54 維轉(zhuǎn)換矩陣。去掉的7 個特征值累積貢獻(xiàn)率僅為0.0001%，降維所帶來的信息損失完全可以忽略。重建矩陣A條件數(shù)為1.68×108，降維重建矩陣Adr的條件數(shù)為8.77×107，重建矩陣的條件數(shù)降低47.8%。

用變換矩陣P左乘式（6）兩端，并令tdr=Pt，可得到PCA 降維后的正問題模型如式（8）所示：

2.2.2 逆問題求解

降維處理可有效地改善重建矩陣的病態(tài)性，但為了避免信息損失所帶來的重建失真，降維后Adr的仍具有嚴(yán)重的病態(tài)性。為此，本文提出的重建算法采用了式（9）所示的迭代正則化獲得式（8）的穩(wěn)定解：

式中：α 為松弛因子；μ 為正則化參數(shù)；I為單位矩陣；k為迭代次數(shù)。求出β后，用式（3）可以計算出N 個網(wǎng)格中心點處的聲慢度。由于聲慢度是聲速的倒數(shù)，故利用式（1）可得到N 個網(wǎng)格中心點處的溫度值。

3 重建誤差評價

本文采用均方根誤差Rrms對重建質(zhì)量進(jìn)行評估：

4 重建算法的仿真數(shù)據(jù)驗證

本文的被測層面為1.3 m×1.3 m 的正方形，每條邊長的四等分點上布置一個聲波收發(fā)器，形成54 條穿越被測區(qū)域的聲線，如圖1 所示。采用本文提出的基于PCA 降維和迭代正則化的重建算法（簡稱pcaIR 法）對單峰、雙峰和三峰模型溫度場分別進(jìn)行仿真重建。單峰模型溫度場的熱點位于（0，0），雙峰模型溫度場的熱點位于（-0.25，-0.25）和（0.25，0.25），三峰模型溫度場的熱點位于（-0.35，-0.35），（-0.35，0.35）和（0.2，-0.1）。作為對比算法，基于奇異值分解的直接正則化算法（簡稱svdDR法）和最小二乘算法（LSM）的重建結(jié)果也在本文中給出。

圖1 聲波收發(fā)器及所形成的聲波路徑Fig.1 Acoustic transceivers and acoustic paths formed

LSM 法[9]用式（11）重建出N 維聲慢度向量S：

式中：L=（Lki）k=1，…，M，i=1，…，N，Lki表示第k（k=1，2，…，M）條聲波路徑在第i（i=1，2，…，N）個網(wǎng)格內(nèi)的路徑長度，t=（t1，…，tM）T為M 維聲波飛行時間向量。

svdDR 法[1，4]基于重建矩陣A的奇異值分解和Tikhonov 正則化獲得式（6）的穩(wěn)定解，如式（12）所示：

式中：uj、vj分別是A的左、右奇異值向量；σ1≥σ2≥…≥σp＞0 是A 的奇異值；p 是非零奇異值的數(shù)目；μ（μ＞0）是正則化參數(shù)。

用LSM 法重建時，被測區(qū)域均勻地劃分為5×5=25 個網(wǎng)格。采用svdDR 法和本文提出的pcaIR 法時，被測區(qū)域劃分成9×9=81 個網(wǎng)格，徑向基函數(shù)形狀參數(shù)為0.0001。svdDR 法的正則化參數(shù)為1×1010。pcaIR 法的正則化參數(shù)為1×109，松弛因子α 為0.9，迭代次數(shù)k 為13 次。為更好地表述復(fù)雜溫度場，本文給出的重建圖像均采用三次樣條插值，細(xì)化后的像素為51×51=2601，即N′＝2601。

對聲波飛行時間數(shù)據(jù)添加標(biāo)準(zhǔn)差為1×10-5的高斯白噪聲。三種算法對應(yīng)的重建誤差如表1 所示。模型溫度場以及3 種重建算法給出的重建圖像如圖2 所示。由表1 和圖2 可以看出：①三幅重建圖像都正確地體現(xiàn)了模型溫度場的基本特征；②本文算法獲得的重建圖像與模型溫度場最接近；③本文算法對應(yīng)的重建誤差最小，與LSM 法、svdDR 法相比，本文算法的重建誤差最高可降低86.62%和29.1%；④LSM 法、svdDR 法和本文算法的重建時間分別為0.38 ms、0.82 ms、3.87 ms，雖本文算法的重建時間最長，但也只有3.87 ms，即3 種算法都是實時性很好的重建算法。

表1 重建誤差比較Tab.1 Comparison of reconstruction errors

圖2 模型溫度場及重建溫度場Fig.2 Model temperature fields and reconstructed temperature fields

5 重建算法的實測數(shù)據(jù)驗證

采用自開發(fā)的聲學(xué)CT 溫度場重建系統(tǒng)對本文算法進(jìn)行了實測數(shù)據(jù)驗證。收發(fā)器（由傳聲器和揚(yáng)聲器組成）布局及形成的54 條有效穿越被測層面的聲波路徑，如圖1 所示。在被測層面下方放置1～2個電加熱器，模擬單熱點、偏置的單熱點和雙熱點溫度場。在每個完整測量周期內(nèi)，揚(yáng)聲器輪流發(fā)聲；傳聲器將接收到的聲波轉(zhuǎn)換為電信號；計算機(jī)通過同步數(shù)據(jù)采集卡采集調(diào)理后的電信號，用基于互相關(guān)的時延估計算法，計算出聲波在54 條聲波路徑上的傳播時間，獲得實測飛行時間數(shù)據(jù)。

采用LSM 法、svdDR 法和本文算法對實測數(shù)據(jù)進(jìn)行了溫度場重建。各算法參數(shù)設(shè)置同仿真數(shù)據(jù)驗證。單熱點、偏置的單熱點和雙熱點溫度場的布置照片以及用3 種重建算法獲得的重建溫度場，如圖3 所示?？紤]到照片形變，熱點位置示意圖也在圖3中給出。用Fluke 54 II 和K 型熱電偶實測的熱點溫度與各算法重建的熱點溫度對比，如表2 所示。由圖3 和表2 可以看出，與LSM、svdDR 法相比，本文算法能更準(zhǔn)確地重建出溫度場特征，熱點溫度也更接近于熱電偶的測量值。

表2 熱點溫度Tab.2 Hot-spot temperature

圖3 實測溫度場Fig.3 Measured temperature fields

6 結(jié)語

本文提出一種基于PCA 降維和迭代正則化的聲學(xué)CT 溫度場重建算法，即pcaIR 法。數(shù)值仿真和實驗驗證表明，與經(jīng)典的最小二乘法、基于奇異值分解的直接正則化法相比，pcaIR 法的重建誤差更低，重建的溫度場更接近于真實分布。pcaIR 法雖是迭代算法，但由于計算簡單，所需要的迭代次數(shù)也很少，所以pcaIR 法的重建時間不足4 ms，遠(yuǎn)小于獲取一幀投影數(shù)據(jù)的數(shù)據(jù)采集時間。因此本文所提pcaIR 法有望應(yīng)用于實際的溫度場重建。