于和放, 胡小梅*

(1.上海大學(xué) 機(jī)電工程與自動(dòng)化學(xué)院,上海 200444; 2.上海市智能制造及機(jī)器人重點(diǎn)實(shí)驗(yàn)室,上海 200444)

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

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

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

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

1.1 TDLAS測溫原理

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

(1)

式中：v為激光頻率(cm-1);a(v)為激光在頻率v處的吸光度;l為激光穿過待測區(qū)域的路徑長度(cm);I0和I1分別為入射光強(qiáng)和出射光強(qiáng);P為待測區(qū)域的壓強(qiáng)(atm,1 atm=101 325 Pa);X為氣體分子的濃度;T(l)為待測區(qū)域的溫度(K);S(T) 為氣體分子對激光的吸收強(qiáng)度(cm-2/atm);φ(v)為吸收譜線的線性函數(shù)。當(dāng)溫度小于2 500 K、激光波長小于 2.5 μm時(shí),線強(qiáng)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測溫通常采用線強(qiáng)比值法,首先根據(jù)硬件條件選取2條合適的吸收譜線,然后可采用掃描頻率波長直接吸收法獲得激光在整個(gè)選定吸收譜線的波長范圍的積分吸收面積：

(3)

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

(4)

1.2 非均勻溫度場重建原理

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

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

(5)

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

A=La

(6)

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

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

(7)

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

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

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

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

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

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

(8)

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

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

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

(9)

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

(10)

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

(11)

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

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

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

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

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

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

(12)

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

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

3 仿真驗(yàn)證

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

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

3.1.1 仿真溫度場的建立

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

圖1 仿真溫度場模型

3.1.2 投影分布設(shè)置

為了驗(yàn)證算法的普適性,分別在扇形投影和平行投影[16]的前提條件下對溫度場進(jìn)行仿真重建。扇形投影分布和平行投影分布的示意圖如圖2所示,激光器從5個(gè)角度發(fā)射激光穿過待測區(qū)域,每個(gè)角度被6個(gè)間距相等的探測器探測,共可獲取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ù)重建算法進(jìn)行迭代修正的過程中,只會(huì)針對有投影穿過的網(wǎng)格進(jìn)行修正,為求解沒有投影穿過的網(wǎng)格的溫度值,本文在迭代的過程中加入了高斯濾波。

3.1.4 氣體吸收分子的選取

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

3.1.5 吸收譜線的選取

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

表1 譜線參數(shù)

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

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

(13)

(14)

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

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

表2 扇形投影分布

表3 平行投影分布

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

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

4 結(jié)束語

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