于和放, 胡小梅*

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

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

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

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

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

1.1 TDLAS測溫原理

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

(1)

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

(2)

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

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

(3)

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

(4)

1.2 非均勻溫度場重建原理

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

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

(5)

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

A=La

(6)

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

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

(7)

式中：k為迭代次數(shù);λ為迭代步長。

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

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

2.1 正則代數(shù)重建算法

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

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

(8)

式中：δ為正則化參數(shù)。

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

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

(9)

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

(10)

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

(11)

式中：μ為單位正則化權(quán)參數(shù);li為第i條投影所穿過的待測區(qū)域的路徑長度。

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

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

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

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

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

(12)

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

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

3 仿真驗證

基于上述溫度場重建原理,應用MATLAB軟件編寫了仿真重建程序,具體輸入?yún)?shù)和仿真重建結(jié)果如下。

3.1 仿真數(shù)據(jù)的設(shè)置

3.1.1 仿真溫度場的建立

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

圖1 仿真溫度場模型

3.1.2 投影分布設(shè)置

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

圖2 投影分布示意圖

3.1.3 待測區(qū)域設(shè)置和網(wǎng)格劃分

待測區(qū)域大小設(shè)置為100 mm×100 mm,網(wǎng)格劃分為19×19,如圖3所示。

圖3 待測區(qū)域設(shè)置和網(wǎng)格劃分示意圖

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

3.1.4 氣體吸收分子的選取

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

3.1.5 吸收譜線的選取

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

表1 譜線參數(shù)

3.2 仿真結(jié)果及誤差分析

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

(13)

(14)

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

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

表2 扇形投影分布

表3 平行投影分布

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

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

4 結(jié)束語

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