李金鍵，王丙寅，齊 琪，張 彪，許傳龍

(東南大學(xué) 大型發(fā)電裝備安全運(yùn)行與智能測(cè)控國(guó)家工程研究中心&能源與環(huán)境學(xué)院，江蘇 南京 210096)

引 言

固體火箭發(fā)動(dòng)機(jī)的燃燒產(chǎn)物經(jīng)噴管加速后噴出的火焰流簡(jiǎn)稱(chēng)為固發(fā)羽流。固發(fā)羽流具有溫度高、流速快、強(qiáng)輻射等特性[1-2]。羽流溫度的準(zhǔn)確測(cè)量能夠?yàn)榘l(fā)動(dòng)機(jī)流場(chǎng)仿真和燃料特性研究等提供重要依據(jù)，并指導(dǎo)工程熱防護(hù)實(shí)踐[3-4]。目前，接觸式測(cè)溫方法受限于探針材料的性能，難以滿(mǎn)足羽流測(cè)溫范圍需求?；诠鈭?chǎng)相機(jī)的非接觸測(cè)溫方法，憑借測(cè)溫范圍寬、響應(yīng)快、可靠性高等優(yōu)點(diǎn)被廣泛應(yīng)用于火焰溫度測(cè)量[6-8]。該方法基本思路是利用火焰自發(fā)輻射的光譜信息構(gòu)建關(guān)于溫度與物性參數(shù)的方程組，通過(guò)求解該方程組，進(jìn)而獲得火焰溫度分布。這類(lèi)問(wèn)題屬于輻射傳輸反問(wèn)題范疇，其主要難點(diǎn)是極小的測(cè)量噪聲可能會(huì)導(dǎo)致極大的求解誤差[9]。因此，選擇有效的重建算法至關(guān)重要。

目前，對(duì)于輻射傳輸反問(wèn)題的求解算法主要分為以下3類(lèi)：第一類(lèi)是正則化算法，例如吉洪諾夫正則化算法、截?cái)嗥娈愔捣纸馑惴?TSVD)[10]、最小二乘QR分解算法(LSQR)[11]、非負(fù)最小二乘算法(NNLS)等。正則化算法由于需要選取合適的正則化策略和預(yù)估輻射特性參數(shù)，難以適用于羽流三維溫度分布重建；第二類(lèi)算法為智能優(yōu)化算法，例如模擬退火算法(SA)[12]和代數(shù)迭代算法(SART)、微粒群算法[13]等。智能優(yōu)化算法隨著重建網(wǎng)格劃分?jǐn)?shù)量的增加，重建時(shí)間急劇增加，難以滿(mǎn)足快速、高分辨的三維溫度分布重建的需求[14]；第三類(lèi)算法為混合算法，該類(lèi)算法是將不同算法各自的優(yōu)勢(shì)有機(jī)結(jié)合，實(shí)現(xiàn)算法間的優(yōu)劣互補(bǔ)，從而提高算法整體性能。例如NNLS-LMBC混合算法[15]、LSQR-CG算法[16]等。當(dāng)前混合算法的研究中普遍存在重建精度較低、魯棒性差等問(wèn)題。因此，有必要構(gòu)建一種魯棒性強(qiáng)、求解精度高的重建算法。

Levenberg-Marquaredt(LM)算法是非線(xiàn)性?xún)?yōu)化算法的一種，其具備優(yōu)秀的全局搜尋能力，能夠快速地接近最優(yōu)解，因此被廣泛應(yīng)用于非線(xiàn)性問(wèn)題的求解[17]。Gauss-Newton(GAUSS)算法具備優(yōu)秀的局部搜尋能力和較快的收斂速度[18]。兩種算法的結(jié)合能夠較好地解決現(xiàn)有三維溫度求解算法重建精度低、魯棒性差等問(wèn)題。

因此，本研究建立了基于單光場(chǎng)相機(jī)的羽流輻射傳輸模型，并提出了一種基于LM-GAUSS混合算法的羽流三維溫度重建方法。通過(guò)在循環(huán)迭代過(guò)程中，對(duì)輻射特性參數(shù)進(jìn)行調(diào)整，進(jìn)而減弱求解問(wèn)題病態(tài)性對(duì)求解過(guò)程的影響。通過(guò)數(shù)值仿真和實(shí)驗(yàn)研究，在一定程度噪聲的條件下，對(duì)比分析了LSQR、NNLS算法與所提混合算法的重建性能，為重建固發(fā)羽流溫度分布提供高可靠性、高精度的重建算法。

1 基于光場(chǎng)成像的羽流測(cè)溫模型

該模型主要由光場(chǎng)相機(jī)、待測(cè)物體兩部組成，如圖1所示。光場(chǎng)相機(jī)由主透鏡、微透鏡陣列與傳感器3部分組成，羽流自發(fā)輻射光線(xiàn)通過(guò)光場(chǎng)成像系統(tǒng)后被傳感器接收，由于背景輻射強(qiáng)度相較于羽流輻射強(qiáng)度可以忽略，通過(guò)單次曝光即可獲得羽流自發(fā)輻射的光線(xiàn)位置、方向與強(qiáng)度等信息。將采集的光譜輻射信息與重建算法相結(jié)合，即可實(shí)現(xiàn)羽流三維溫度分布的重建[7, 8]。

羽流內(nèi)部的輻射傳輸過(guò)程由輻射傳輸方程進(jìn)行描述，方程如下所示：

(1)

式中：Iλ(r,Ω)為r位置處波長(zhǎng)為λ并沿Ω方向傳輸?shù)墓庾V輻射強(qiáng)度，W/(m3·sr)；k(r)為吸收系數(shù)，單位為m-1；δ(r)為散射系數(shù)，m-1；Ф(Ω′,Ω)為羽流在r位置處的散射相函數(shù)；Ibλ(r)為黑體輻射強(qiáng)度，W/(m3·sr)。

為簡(jiǎn)化研究過(guò)程，僅考慮其吸收特性[19-20]。結(jié)合光學(xué)厚度定義，如式(2)所示，將由上述方程推導(dǎo)得到相鄰網(wǎng)格間的輻射傳輸方程，如式(3)所示。

Δτ=kΔL

(2)

Iλ(r+Δr,Ω)=Iλ(r,Ω)·exp[-k(r)·ΔL]+

Ibλ(r)·(1-exp[-k(r)·ΔL])

(3)

式中：Δτ為光線(xiàn)沿Ω方向穿過(guò)網(wǎng)格的光學(xué)厚度，無(wú)量綱；ΔL為光線(xiàn)沿Ω方向穿過(guò)網(wǎng)格的距離，m；Δr為沿Ω入射方向上相鄰網(wǎng)格的距離，m。

若將羽流劃分為V個(gè)網(wǎng)格，并按一定排列方式將各網(wǎng)格編號(hào)(1，2，v，…，V)，光場(chǎng)相機(jī)圖像探測(cè)器共M個(gè)像素接收到輻射光線(xiàn)，并按一定排列方式將各像素編號(hào)(1，2，m，…，M)。對(duì)于任意一個(gè)像素m的輻射光線(xiàn)，其穿過(guò)的羽流網(wǎng)格總數(shù)記為n，則輻射光線(xiàn)方向?qū)⒏骶W(wǎng)格編號(hào)依次記為v1，v2，vi，(vj)，…，vn，結(jié)合式(3)可得到各像素輻射光線(xiàn)的光譜輻射強(qiáng)度，如式(4)與式(5)所示，將N條光線(xiàn)計(jì)算式聯(lián)立可得離散的輻射傳輸方程組，如式(6)所示。溫度可通過(guò)求解方程組后，由黑體輻射定律獲得，如式(7)所示。

(4)

(5)

(6)

A·Ibλ=Iλ,ccd

(7)

2 LM-GAUSS算法

LM-GAUSS混合算法的基本原理是在循環(huán)迭代過(guò)程中，通過(guò)比較當(dāng)前目標(biāo)函數(shù)與梯度向量的無(wú)窮范數(shù)來(lái)選擇該次迭代需使用的算法，并使用所選擇的算法計(jì)算求解向量該次迭代需前進(jìn)的步長(zhǎng)，直到目標(biāo)函數(shù)滿(mǎn)足終止條件或滿(mǎn)足最大迭代次數(shù)為止[21]。兩種算法的結(jié)合能夠保證求解精度的同時(shí)，提高了算法魯棒性。

結(jié)合羽流溫度重建問(wèn)題，目標(biāo)函數(shù)向量f和目標(biāo)函數(shù)F的定義如式(8)、式(9)所示，求解過(guò)程中算法的選擇策略如式(10)所示：

f(x)=A·Ibλ-Iλccd

(8)

(9)

(10)

式中：J為雅可比矩陣；hGAUSS為GAUSS算法計(jì)算的前進(jìn)步長(zhǎng)；hLM為L(zhǎng)M算法計(jì)算的前進(jìn)步長(zhǎng)；F′(x)為目標(biāo)函數(shù)的梯度向量；B為GAUSS算法的系數(shù)矩陣；μ為L(zhǎng)M算法的阻尼系數(shù)；I為單位矩陣。LM-GAUSS混合算法步驟與流程圖如下：

圖2 LM-GAUSS算法流程圖Fig.2 Flow chart of LM-Gauss algorithm

步驟一：對(duì)溫度T隨機(jī)賦值，將吸收系數(shù)k與溫度T組成初始迭代向量x，設(shè)定最大迭代次數(shù)itermax與算法截止值ε1，ε2。

步驟二：設(shè)定初始迭代的算法類(lèi)型Method=LM；迭代狀態(tài)數(shù)Better=false；滿(mǎn)足選擇GAUSS算法條件的次數(shù)計(jì)數(shù)器值count=0；對(duì)LM算法中阻尼系數(shù)的調(diào)整因子q=2；GAUSS算法初始系數(shù)矩陣B=I；迭代次數(shù)計(jì)數(shù)器值iter=0；GAUSS算法中解向量x前進(jìn)步長(zhǎng)調(diào)整因子Δ=0。

步驟三：計(jì)算當(dāng)前目標(biāo)函數(shù)F、雅可比矩陣J、函數(shù)向量f、混合算法終止條件數(shù)Found與LM算法的阻尼因子μ。其中，阻尼因子取轉(zhuǎn)置雅克矩陣與其自身相乘后矩陣對(duì)角線(xiàn)元素的最大值，如式(11)與式(12)所示。

Found=‖F(xiàn)′(x)‖∞≤ε1

(11)

u=max{diag(JTJ)}

(12)

步驟四：判斷Method是否為L(zhǎng)M，若是則繼續(xù)下一步驟，否則轉(zhuǎn)到步驟六。

步驟五：將相關(guān)參數(shù)Found、x、count、μ、q、J、f輸入LM算法模塊進(jìn)行計(jì)算，返回參數(shù)xnew、Found、Better、Method、count、μ、q、Δ后轉(zhuǎn)到步驟七。

步驟六：將相關(guān)參數(shù)Found、x、B、Δ輸入GAUSS算法模塊進(jìn)行計(jì)算并返回獲得的參數(shù)(xnew、Found、Better、Method、count、Δ)。

步驟七：根據(jù)BFGS策略對(duì)系數(shù)矩陣B進(jìn)行更新。

步驟八：判斷狀態(tài)數(shù)Better是否為true。若是，則更新當(dāng)前解向量x=xnew；

步驟九：判斷混合算法條件數(shù)是否為真或iter>itermax。若混合算法條件數(shù)為1，則將當(dāng)前解向量并轉(zhuǎn)化為溫度重建值T，繼續(xù)下一步驟，否則iter=iter+1，轉(zhuǎn)到步驟三。

步驟十：輸出重建溫度值T，結(jié)束計(jì)算。

3 仿真計(jì)算結(jié)果與討論

通過(guò)仿真計(jì)算對(duì)混合算法性能進(jìn)行研究，將半徑Rmax為8mm，高度Hmax為30mm的圓柱羽流沿軸向、徑向和周向進(jìn)行劃分，如圖3所示。同時(shí)，選取610nm和530nm為中心波段的兩組光譜輻射信息構(gòu)建離散化的輻射傳輸方程組，設(shè)定羽流溫度T(K)與吸收系數(shù)k(m-1)的分布函數(shù)，如式(13)、式(14)所示，圖3中AD為軸向方向，RD為半徑方向，CD為周向方向。

(13)

(14)

圖3 網(wǎng)格劃分與設(shè)定溫度分布圖Fig.3 Grid partition and setting temperature distribution

為研究噪聲對(duì)算法性能的影響，在添加2%至10%不同程度的高斯白噪聲后，使用混合算法對(duì)三維溫度進(jìn)行重建，并通過(guò)比較平均重建誤差對(duì)算法的準(zhǔn)確性與魯棒性進(jìn)行評(píng)估，噪聲添加公式與平均重建誤差的定義如式(15)～式(17)所示：

(15)

(16)

(17)

式中：In,ccd為第n根光線(xiàn)的輻射強(qiáng)度，W/(m3·sr)；In,ccd,exact為第n根光線(xiàn)準(zhǔn)確的輻射強(qiáng)度值；noiseI為對(duì)輻射強(qiáng)度的高斯白噪聲強(qiáng)度，無(wú)量綱；km為第m個(gè)網(wǎng)格的吸收系數(shù)；km,exact為網(wǎng)格準(zhǔn)確的吸收系數(shù)值；noisek為對(duì)吸收系數(shù)的高斯白噪聲強(qiáng)度；ζ為均值為0，標(biāo)準(zhǔn)差為1的正態(tài)分布隨機(jī)數(shù)；δT為溫度平均重建誤差。

3.1 輻射強(qiáng)度噪聲下算法性能對(duì)比

對(duì)輻射強(qiáng)度添加不同程度的噪聲后，使用LSQR、NNLS與LM-GAUSS這3種算法對(duì)溫度進(jìn)行重建，結(jié)果如圖4示。

圖4 輻射強(qiáng)度噪聲下算法重建誤差圖Fig.4 Reconstruction error under radiation intensity noise

由圖4可知，3種算法均能在無(wú)噪聲條件下準(zhǔn)確地重建羽流三維溫度分布。隨著輻射強(qiáng)度噪聲增加，由于LSQR算法在求解過(guò)程中無(wú)法保證溫度非負(fù)性的特征，因此，算法的平均重建誤差高于NNLS算法。當(dāng)輻射強(qiáng)度噪聲大于2%時(shí)，LSQR算法的平均重建誤差高于10%，NNLS算法的平均重建誤差高于7%。其主要原因是待求解問(wèn)題屬于病態(tài)問(wèn)題，系數(shù)矩陣過(guò)大的條件數(shù)導(dǎo)致在求解過(guò)程中極小的強(qiáng)度噪聲會(huì)造成極大的重建誤差。

隨著強(qiáng)度噪聲的增加，LM-GAUSS算法重建誤差略微增加，最大平均重建誤差為0.6%。其證明了在含有輻射強(qiáng)度噪聲重建時(shí)，LM-GAUSS算法具有比LSQR、NNLS算法更高的重建精度與更強(qiáng)的魯棒性。

3.2 吸收系數(shù)噪聲下算法性能對(duì)比

在求解離散化的輻射傳輸方程時(shí)，吸收系數(shù)會(huì)含有一定程度的噪聲。為討論該部分噪聲對(duì)算法性能的影響，使用LSQR、NNLS以及LM-GAUSS算法在添加吸收系數(shù)噪聲的條件下對(duì)羽流三維溫度進(jìn)行重建，結(jié)果如圖5所示。

圖5 吸收系數(shù)噪聲下算法重建誤差圖Fig.5 Reconstruction error under absorption coefficient noise

隨著吸收系數(shù)噪聲增加，LSQR算法的重建誤差快速增加至40%以上，NNLS算法的重建誤差快速增加到20%以上。主要原因是對(duì)吸收系數(shù)添加噪聲后，系數(shù)矩陣產(chǎn)生了較大偏差，加劇了問(wèn)題的病態(tài)程度，使得重建誤差遠(yuǎn)大添加輻射強(qiáng)度噪聲時(shí)的重建誤差，此時(shí)上述兩種算法已經(jīng)無(wú)法在含有吸收系數(shù)噪聲的條件下對(duì)溫度進(jìn)行求解。

LM-GAUSS算法在計(jì)算過(guò)程中，對(duì)含有噪聲的吸收系數(shù)進(jìn)行調(diào)整，能夠減弱吸收系數(shù)噪聲對(duì)問(wèn)題求解的影響。隨著吸收系數(shù)噪聲的增加，算法重建誤差略微增加，最大平均重建誤差為0.7%。上述結(jié)果表明，相比于LSQR、NNLS算法對(duì)吸收系數(shù)噪聲過(guò)于敏感的問(wèn)題，LM-GAUSS算法具有更好的重建性能，證明了LM-GAUSS算法是優(yōu)于NNLS、LSQR算法的選擇。

3.3 不同噪聲水平對(duì)LM-GAUSS算法的影響

在實(shí)際過(guò)程中，吸收系數(shù)噪聲與輻射強(qiáng)度噪聲同時(shí)存在，因此，設(shè)立以下4種工況，如表1所示，使用LM-GAUSS算法對(duì)溫度進(jìn)行重建，重建后的溫度分布如圖6所示。

表1 數(shù)值計(jì)算的工況設(shè)計(jì)表 Table 1 Working condition design example

圖6 不同噪聲條件下溫度重建分布圖Fig.6 Temperature distribution under different noise conditions

在不同的噪聲條件下，LM-GAUSS算法均能夠較準(zhǔn)確地重建羽流內(nèi)部的溫度分布，4種工況對(duì)應(yīng)的平均重建誤差分別為1.67%、2.1%、2.58%和4.66%。隨著噪聲的增加，靠近中軸線(xiàn)位置的重建誤差增加，但溫度分布與設(shè)定趨勢(shì)保持一致。羽流頭部的溫度與設(shè)定溫度相比偏高，可能與穿過(guò)頭部網(wǎng)格的光線(xiàn)數(shù)量較少有關(guān)。在吸收系數(shù)噪聲與輻射強(qiáng)度噪聲條件同時(shí)存在時(shí)，混合算法的重建誤差優(yōu)于5%，保持在可接受范圍內(nèi)。

3.4 網(wǎng)格數(shù)量對(duì)LM-GAUSS算法性能影響

在三維溫度測(cè)量過(guò)程中，網(wǎng)格劃分?jǐn)?shù)量越多，分辨率越高。因此，需要討論網(wǎng)格數(shù)量對(duì)LM-GAUSS算法重建性能的影響。本研究設(shè)計(jì)了3種不同的網(wǎng)格劃分形式，如表2所示，在分別含有吸收系數(shù)噪聲與輻射強(qiáng)度噪聲的條件下對(duì)溫度進(jìn)行重建，重建結(jié)果如圖7所示。

表2 網(wǎng)格劃分工況設(shè)計(jì)表 Table 2 Design of grid number example

圖7 輻射強(qiáng)度噪聲下網(wǎng)格數(shù)量對(duì)重建誤差的影響Fig.7 Influence of mesh number on the reconstruction error under radiation noise

圖8 吸收系數(shù)噪聲下網(wǎng)格數(shù)量對(duì)重建誤差的影響Fig.8 Influence of mesh number on the reconstruction error under absorption coefficient noise

由圖7與圖8可知，在噪聲不變的條件下，隨著網(wǎng)格數(shù)量增加，平均重建誤差略微增加。其主要原因是網(wǎng)格數(shù)量越多，對(duì)應(yīng)解向量的維度越大，問(wèn)題的病態(tài)性越強(qiáng)。對(duì)比兩種情況可發(fā)現(xiàn)LM-GAUSS算法在網(wǎng)格數(shù)量增加過(guò)程中，對(duì)輻射強(qiáng)度噪聲更為敏感，其重建誤差增加程度更大。當(dāng)網(wǎng)格劃分較多時(shí)，輻射強(qiáng)度噪聲下混合算法的重建誤差低于3%，吸收系數(shù)噪聲下的重建誤差低于1.5%，證明了LM-GAUSS算法能夠在高分辨率條件下重建三維溫度分布。

綜上所述，通過(guò)比較不同噪聲條件下傳統(tǒng)LSQR、NNLS算法與LM-GAUSS算法的性能，驗(yàn)證了LM-GAUSS算法具備更高的重建精度以及更強(qiáng)的魯棒性，并討論了網(wǎng)格數(shù)量對(duì)LM-GAUSS算法性能的影響，結(jié)果表明，隨著網(wǎng)格數(shù)量的增加，重建誤差略有增加，但仍保持較高的重建精度，可為后續(xù)固體火箭發(fā)動(dòng)機(jī)熱試車(chē)試驗(yàn)提供算法依據(jù)。

4 試驗(yàn)結(jié)果與分析

為了進(jìn)一步研究算法在重建實(shí)際火焰的性能，基于課題組的燃燒平臺(tái)，對(duì)乙烯層流擴(kuò)散火焰進(jìn)行溫度重建，實(shí)驗(yàn)裝置如圖9所示。

圖9 乙烯擴(kuò)散火焰拍攝裝置Fig.9 Ethylene diffusion flame temperature measurement device

燃燒工況參數(shù)如下，乙烯流率VC2H4為0.08L/min，空氣流率VO2為3L/min。光場(chǎng)相機(jī)由主鏡頭(Nikon NIKKOR 50mm f/1.8D)、微透鏡陣列(F數(shù)為4.2，透鏡單元尺寸為100μm×100μm)、相機(jī)(分辨率為3312×2488，像素大小為5μm)組成。使用光場(chǎng)相機(jī)對(duì)火焰進(jìn)行拍攝，通過(guò)光場(chǎng)相機(jī)拍攝到的光場(chǎng)圖像，如圖10(a)所示。再利用黑體爐進(jìn)行輻射強(qiáng)度標(biāo)定，獲得光線(xiàn)的輻射強(qiáng)度分布，將火焰區(qū)域劃分8×8×8的網(wǎng)格，重建得到溫度分布，如圖10(b)所示。將熱電偶(型號(hào)為Omega P13R-020)測(cè)量火焰中心位置的溫度作為參考點(diǎn)，驗(yàn)證算法重建精度。

圖10 光場(chǎng)相機(jī)的火焰圖像和重建溫度分布圖Fig.10 Flame image of light field camera and reconstructed temperature distribution

由圖10(b)可知，重建火焰溫度范圍在700～2200K范圍內(nèi)。火焰中心區(qū)域算法重建溫度為1431K，熱電偶測(cè)溫為1350K，相對(duì)誤差為6%?；鹧娴闹亟ㄗ罡邷囟葹?214K，最低溫度為700K。隨著高度的增加，火焰高溫區(qū)域的半徑逐漸減小，符合實(shí)際火焰溫度分布的特點(diǎn)，證明LM-GAUSS算法能夠較好地反應(yīng)實(shí)際火焰的三維溫度分布特點(diǎn)，驗(yàn)證了其重建實(shí)際火焰的可行性。

5 結(jié) 論

(1)對(duì)比LSQR、NNLS與LM-GAUSS算法的重建結(jié)果，驗(yàn)證了LM-GAUSS混合算法重建精度高、魯棒性強(qiáng)的特點(diǎn)。分析了網(wǎng)格數(shù)量對(duì)算法重建精度的影響，結(jié)果顯示隨著網(wǎng)格數(shù)量的增加，重建誤差略有上升，但仍?xún)?yōu)于5%。

(2)使用LM-GAUSS算法對(duì)乙烯擴(kuò)散火焰的溫度場(chǎng)進(jìn)行重建試驗(yàn)，證明了算法在重建實(shí)際火焰的可行性，為后續(xù)開(kāi)展固體火箭發(fā)動(dòng)機(jī)溫度場(chǎng)現(xiàn)場(chǎng)測(cè)量提供了算法支撐。