基于正則代數重建算法的二維溫度場重建

于和放, 胡小梅*

(1.上海大學 機電工程與自動化學院,上海 200444; 2.上海市智能制造及機器人重點實驗室,上海 200444)

可調諧半導體激光吸收光譜斷層診斷技術(Tunable Diode Laser Absorption Tomography,TDLAT)是可調諧二極管激光吸收光譜技術 (Tunable Diode Laser Absorption Spectroscopy,TDLAS)與計算機斷層成像(Computed Tomography,CT)技術相結合的一種新型的流場參數測量技術,可以實現對非均勻溫度場的測量。因其具有非侵入式、靈敏度高、抗噪聲能力強、實時性好、測量環境適應性強等優勢,目前,已廣泛應用于燃燒流場的診斷研究[1-3]中。

在基于TDLAT的溫度場重建中,重建算法對重建結果的精度起著關鍵作用。目前常用的重建算法主要包括解析類重建算法和迭代類重建算法。解析類重建算法在完備投影的條件下,可獲得較好的重建結果[4],但是由于硬件條件的限制,往往無法實現完備投影。迭代類重建算法可以在迭代過程中加入先驗信息,可以應用在不完備投影的溫度場重建中[5-6]。夏暉暉等[7]以水汽作為氣體吸收分子,使用代數重建算法(Algebraic Reconstruction Technique,ART)對燃燒溫度場和濃度場分布進行了仿真重建和實驗研究,但是網格劃分數量較少,重建結果的分辨率較低。Jeon等[8-11]提出了3種改進的ART,并進行了仿真驗證,給出了不同算法的適用范圍,并在后續幾年中將算法應用在實驗研究中。Liu等[12]在扇形投影分布下使用改進的Landweber算法對燃燒溫度場進行了仿真重建和實驗研究,實驗結果表明重建溫度場與預期分布吻合較好,表明了該算法檢測燃燒火焰的有效性。

上述算法本質上都是求最小二乘解或近似最小二乘解,在重建的過程中會引入噪聲和離散誤差。為了能夠有效抑制噪聲,平衡離散誤差,本文在ART的基礎上提出了基于留一交叉驗證的正則ART,并使用典型溫度場和典型投影分布對所提出的算法的重建精度和穩定性進行了仿真驗證。

1 TDLAS測溫和溫度場重建原理

1.1 TDLAS測溫原理

TDLAS的基本原理是激光穿過流場后的衰減程度不同,根據Beer-Lambert定律[13]可獲得激光衰減前后光強的關系[14],即

(1)

式中：v為激光頻率(cm-1);a(v)為激光在頻率v處的吸光度;l為激光穿過待測區域的路徑長度(cm);I0和I1分別為入射光強和出射光強;P為待測區域的壓強(atm,1 atm=101 325 Pa);X為氣體分子的濃度;T(l)為待測區域的溫度(K);S(T) 為氣體分子對激光的吸收強度(cm-2/atm);φ(v)為吸收譜線的線性函數。當溫度小于2 500 K、激光波長小于 2.5 μm時,線強S(T)可表示為

(2)

式中：T0為參考溫度296 K;Q(T)為配分函數;h為普朗克常數(J·S);c為光速(cm/s);E″為低狀態能級能量(cm-1);k為玻爾茲曼常數(J/K)。以上這些參數可以通過查詢HITRAN 光譜數據庫獲取。

TDLAS測溫通常采用線強比值法,首先根據硬件條件選取2條合適的吸收譜線,然后可采用掃描頻率波長直接吸收法獲得激光在整個選定吸收譜線的波長范圍的積分吸收面積：

(3)

線強比值法就是在同一條激光投影下的一個測量周期中,獲取該氣體的2個積分吸收面積。一次測量周期的時間很短,所以默認同一測量周期下,溫度、濃度和壓強不變,因此可以通過求解2個積分吸收面積的比值消除濃度和壓強的影響。經過變換后可以獲取待測區域的平均溫度為

(4)

1.2 非均勻溫度場重建原理

上述方法計算的溫度結果為投影上氣體溫度的均值。而實際待測區域多為非均勻溫度場,為了能獲取非均勻溫度場的溫度分布,需要將TDLAS和CT技術相結合。TDLAT的原理是獲取待測溫度場多個方向的投影數據,然后使用相應的重建算法重建溫度場?；贏RT的溫度場重建過程如下。

首先將待測區域離散成N個網格,并且認為每個網格內的溫度、壓力和濃度都是定值,離散化后積分吸收面積即投影值可以表示為

(5)

式中：i和j為投影和網格的編號;aj為第j個網格內的積分吸光度;Lij為第i條激光束在第j個網格內的路徑長度。表示成矩陣的形式為

A=La

(6)

式中：A為投影向量;L為系數矩陣;a為網格積分吸光度向量,是未知量。

然后采用ART對上述方程組進行求解。其基本思想是：先給定初始解,在沒有先驗信息的條件下一般初始解為0;然后計算當前計算投影值和測量投影值的殘差,沿著投影的傳播路徑進行反投影修正;之后不斷迭代直到殘差值滿足重建精度要求為止。ART的特點是每次迭代修正時只使用一個投影,因此只對該投影穿過的網格單元值進行修正。使用第i條投影對第j個網格進行修正的迭代公式為

(7)

式中：k為迭代次數;λ為迭代步長。

求解出每個網格的積分吸光度后,再使用線強比值法求解出每個網格的平均溫度,即可獲取待測場的溫度分布。

2 基于留一交叉驗證的正則代數重建算法

2.1 正則代數重建算法

在實際測量投影值時,由于設備和環境的影響,不可避免地會存在一定的噪聲。并且根據非均勻溫度場的重建原理,在進行重建時,需要將待測區域劃分成N個網格,因此基于網格計算出的投影值與測量投影值之間存在一定的離散誤差。所以在使用ART對方程組進行求解時,離散誤差和投影噪聲會隨著迭代修正而引入。為了減小重建算法對噪聲的敏感度,平衡離散誤差對結果帶來的影響,本文對ART進行了改進,提出正則代數重建算法。

ART實質上求解的是近似最小二乘解,所以對噪聲和離散誤差比較敏感,一般的解決辦法是在最小二乘準則的基礎上增添正則化項,將極小化目標函數改為

(8)

式中：δ為正則化參數。

正則化項的作用就是對解做一些限制,使解更加平滑,δ的作用是控制平滑的程度。采用梯度下降法對目標函數進行求解,將負梯度方向作為迭代方向,可以推導出正則Landweber算法的迭代格式為

ak+1=ak+λ[LT(A-Lak)+δak]

(9)

由式(9)可知,正則Landweber算法對投影的噪聲水平和離散誤差進行統一衡量,但是在實際情況下,每條投影穿過的重建區域都不相同,因此每條投影的噪聲水平和離散誤差的大小也不一致。所以有必要對每條投影的噪聲水平和離散誤差分別進行衡量。ART的特點就是每次迭代修正時只使用一個投影,可以實現對投影噪聲水平和離散誤差的分別衡量。由式(9)可以推導出正則代數重建算法的格式為

(10)

一條投影穿過待測區域的路徑越長,該投影路徑上的溫度梯度大概率越大,即離散誤差越大。所以每條投影的正則化權重由投影穿過的待測區域路徑的長度和單位正則化參數共同決定,因此正則代數重建算法的迭代公式可更改為

(11)

式中：μ為單位正則化權參數;li為第i條投影所穿過的待測區域的路徑長度。

正則代數重建算法的求解質量依賴于單位正則化參數的選擇,當正則化項的權重過小時不能很好地抑制噪聲和離散誤差,當正則化項權重較大時,會使得重建結果過于平滑。為了選取合適的正則化參數,本文采用留一交叉驗證法對正則化參數進行選取。

2.2 基于留一交叉驗證的正則化參數選取

交叉驗證也稱為循環估計,是一種可以評估參數的統計學方法,交叉驗證法評估參數的基本思想是將原始樣本數據集分割成2個子集,其中一個子集數據用于求解,另一個子集的數據用于驗證。對于本文來講,首先將投影數據分成2個數據集,設置不同的正則化參數,使用其中一個數據集重建溫度場,然后使用重建后的結果對另一個數據集的投影值進行計算,將計算投影值和測量投影值進行比較,獲取殘差值,選取使殘差值最小的正則化參數。

由于投影數量較少,在使用普通交叉驗證法選取正則化參數時,根據數據集拆分方式的不同,重建結果的差異性較大,使得求解的正則化參數不穩定。

為了減少普通交叉驗證方法的不穩定性,可以將數據集劃分為m個子集,每次取一個集合用于驗證,剩余m-1個集合用于求解,總共須進行m次交叉驗證,對這m次交叉驗證結果進行累加再取平均值作為最終的驗證結果。這種方法稱為m折交叉驗證法。為了保證在交叉驗證中求解的結果與使用整個數據集求解的結果最為接近,本文在劃分數據集時,使m等于數據集中投影數據的個數,即每次交叉驗證時,驗證集只有1個數據用于驗證,剩余的數據用于重建溫度場,這種方法被稱為留一交叉驗證法。該方法的極小化目標函數為

(12)

雖然留一交叉驗證法所需的計算成本較大,但可以最大限度地利用原始投影數據選出合適的正則化參數。

至此,對正則化參數的選取問題變成了目標函數為單變量函數且函數曲線為單峰的一維最優化問題,但是目標函數無法求導,所以本文采用黃金分割法求解近似最優值,其基本思想是首先確定最優點所在的區間,以固定區間縮短率0.618不斷縮小最優點所在的區間范圍,設置相應的閾值,區間長度小于閾值時,停止分割,取區間的中點作為近似最優點,將此時正則化權重值作為近似最優正則化參數。

3 仿真驗證

基于上述溫度場重建原理,應用MATLAB軟件編寫了仿真重建程序,具體輸入參數和仿真重建結果如下。

3.1 仿真數據的設置

3.1.1 仿真溫度場的建立

仿真模型采用高斯單峰對稱和高斯單峰偏置這2種溫度模型,如圖1所示,溫度范圍為500～1 000 K。

圖1 仿真溫度場模型

3.1.2 投影分布設置

為了驗證算法的普適性,分別在扇形投影和平行投影[16]的前提條件下對溫度場進行仿真重建。扇形投影分布和平行投影分布的示意圖如圖2所示,激光器從5個角度發射激光穿過待測區域,每個角度被6個間距相等的探測器探測,共可獲取30條有效投影數據。

圖2 投影分布示意圖

3.1.3 待測區域設置和網格劃分

待測區域大小設置為100 mm×100 mm,網格劃分為19×19,如圖3所示。

圖3 待測區域設置和網格劃分示意圖

從圖3可以看出,有部分網格沒有被投影穿過,而采用正則代數重建算法進行迭代修正的過程中,只會針對有投影穿過的網格進行修正,為求解沒有投影穿過的網格的溫度值,本文在迭代的過程中加入了高斯濾波。

3.1.4 氣體吸收分子的選取

H2O作為主要的燃燒產物,在1～3 μm波段有大量吸收,并且該波段的激光器較為成熟,所以H2O常作為測量流場溫度的媒介,因此在進行仿真溫度場重建時采用H2O作為測量媒介,并默認待測區域H2O的濃度為均勻分布。

3.1.5 吸收譜線的選取

在選取吸收譜線時首先要考慮硬件條件,選擇的譜線要有利于測量,盡量減少測量誤差,并且要保證線強比值R和溫度值T是一一對應的關系,即線強比是溫度的單調函數。其次要考慮所選擇的譜線組對待測區域的溫度區間有較高的測溫靈敏度[15]。根據仿真溫度場的設置,要保證譜線組在500～1 000 K的溫度范圍內有較高的測溫靈敏度。本文仿真測試所使用的一組吸收譜線的參數如表1所示。

表1 譜線參數

3.2 仿真結果及誤差分析

為表征溫度場的重建精度,對重建結果計算平均相對誤差Rave和最大相對誤差RE,max。

(13)

(14)

式中：Tm(j)和T(j)分別為模型溫度場和重建溫度第j個網格的的溫度值。

在常壓條件下,基于上述參數分別采用ART、正則Landweber算法和正則代數重建算法對溫度場進行重建,重建結果的誤差如表2和表3所示。

表2 扇形投影分布

表3 平行投影分布

從表2和表3中的數據可以看出,對上述2種溫度場而言,相比于ART,正則Landweber算法和正則代數重建算法的重建精度都有一定的提高。其中正則Landweber算法對投影的噪聲水平和離散誤差進行了統一衡量,優化效果并不明顯;正則代數重建算法根據投影穿過待測區域路徑的長度和單位正則化參數,動態地調整每條投影的正則化權重,實現了對每條投影的噪聲水平和離散誤差的分別衡量,使用該算法在扇形投影和平行投影下對上述2種溫度場的重建精度都有較明顯的提高。

在計算機處理器為AMD Ryzen 7 6800H的硬件條件下,基于正則代數重建算法的溫度場重建平均耗時為110.4 s,其中正則化參數的選取耗時占比約為99.68%,原因是采用了留一交叉驗證法,計算量較大。

4 結束語

本文依據最優化理論,對ART進行了改進,提出了正則代數重建算法,并使用該算法對單峰對稱和單峰偏置溫度場進行了仿真重建。在不同投影分布情況下的重建結果表明,相比ART和正則Landweber算法,該算法的重建精度有較為明顯的提高,并且穩定性較好。但是該算法的正則化參數選取采用了留一交叉驗證法,該方法計算成本較大,不能滿足實時重建的要求,因此須進一步研究快速選擇正則化參數的方法。