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

白 云，顏 華，魏元焜

（沈陽工業大學 信息科學與工程學院，沈陽 110870）

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

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

1 聲學CT 溫度場重建原理

聲波在氣體介質中的傳播速度c 與氣體介質的絕對溫度間T 間的關系為

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

用聲學CT 法重建溫度場需事先將多組聲波收發器布置在待測區域周圍。系統工作時依據收發器采集到的聲波數據，計算聲波在收發器間的飛行時間，通過預設的重建算法，計算被測層面的聲慢度（聲速的倒數）分布，進而利用式（1）得到溫度分布。

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

2.1 正問題建模

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

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

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

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

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

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

2.2 逆問題求解

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

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

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

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

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

（3）令Adr=PA得到降維后的K×N 維重建矩陣Adr。

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

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

2.2.2 逆問題求解

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

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

3 重建誤差評價

本文采用均方根誤差Rrms對重建質量進行評估：

4 重建算法的仿真數據驗證

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

圖1 聲波收發器及所形成的聲波路徑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）個網格內的路徑長度，t=（t1，…，tM）T為M 維聲波飛行時間向量。

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

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

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

對聲波飛行時間數據添加標準差為1×10-5的高斯白噪聲。三種算法對應的重建誤差如表1 所示。模型溫度場以及3 種重建算法給出的重建圖像如圖2 所示。由表1 和圖2 可以看出：①三幅重建圖像都正確地體現了模型溫度場的基本特征；②本文算法獲得的重建圖像與模型溫度場最接近；③本文算法對應的重建誤差最小，與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 重建算法的實測數據驗證

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

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

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

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

6 結語

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