黃興，齊宏，*，牛志田，任亞濤，阮立明

（1.哈爾濱工業(yè)大學(xué) 能源科學(xué)與工程學(xué)院，哈爾濱150001； 2.空天熱物理工業(yè)和信息化部重點(diǎn)實(shí)驗(yàn)室，哈爾濱150001）

高溫燃燒廣泛存在于航空航天、能源動(dòng)力、鋼鐵冶金、化工等領(lǐng)域，發(fā)動(dòng)機(jī)、燃?xì)廨啓C(jī)、內(nèi)燃機(jī)等高溫設(shè)備中均存在燃燒現(xiàn)象。火焰作為燃料燃燒時(shí)生成的劇烈的光熱現(xiàn)象，其溫度場直接表征了燃燒的狀態(tài)，如燃燒的穩(wěn)定性、燃燒效率、反應(yīng)速度、火焰三維結(jié)構(gòu)等。準(zhǔn)確測量燃燒過程中火焰的溫度分布，不僅可以更好地認(rèn)識燃燒產(chǎn)物的生成機(jī)理與火焰中的熱量傳遞和吸收過程，還可以為生產(chǎn)設(shè)備的安全運(yùn)行、燃料熱能的充分利用及燃燒污染物的有效控制提供數(shù)據(jù)支撐。如在航空發(fā)動(dòng)機(jī)的設(shè)計(jì)中，燃燒室出口溫度是一項(xiàng)很重要的設(shè)計(jì)參數(shù)，可作為衡量發(fā)動(dòng)機(jī)性能的一項(xiàng)重要指標(biāo)，對其溫度進(jìn)行測量可為發(fā)動(dòng)機(jī)設(shè)計(jì)提供參考。因此，對火焰空間測溫技術(shù)的研究有十分顯著的現(xiàn)實(shí)意義。

目前，常用的測溫手段主要分為接觸式測量和非接觸式測量2類。對于接觸式測量，測溫元件需要接觸火焰，這會(huì)導(dǎo)致火焰流場受到干擾，且接觸式測量為點(diǎn)測量方式，無法實(shí)現(xiàn)整個(gè)溫度場的測量。近年來，基于輻射成像的火焰測溫技術(shù)受到了國內(nèi)外學(xué)者的廣泛關(guān)注［1-3］。該方法具有非侵入式、無需外加激勵(lì)源等優(yōu)點(diǎn)，而且可以通過圖像記錄整個(gè)火焰輻射場的信息，進(jìn)而實(shí)現(xiàn)場參數(shù)的重建［4-7］。根據(jù)成像元件的數(shù)量，輻射成像測溫法可以分為單相機(jī)和多相機(jī)2類。單相機(jī)系統(tǒng)結(jié)構(gòu)簡單、易于實(shí)現(xiàn)，但由于成像角度受限，難以實(shí)現(xiàn)復(fù)雜形狀火焰溫度場的測量；多相機(jī)系統(tǒng)是在火焰周圍不同位置和角度布置多個(gè)相機(jī)進(jìn)行同時(shí)成像，以獲得多個(gè)角度的輻射圖像［8-10］，該系統(tǒng)空間分辨率高，但是系統(tǒng)復(fù)雜成本高，而且操作更為繁瑣。2005年，斯坦福大學(xué)的華裔學(xué)者Ng等［11］將光場成像的理論和技術(shù)應(yīng)用于商業(yè)化，也為火焰輻射成像測溫法提供了新思路。光場相機(jī)內(nèi)部特有的微透鏡陣列結(jié)構(gòu)，可以通過單相機(jī)一次曝光記錄全場的輻射光場信息，包括輻射強(qiáng)度的大小、位置及角度。該技術(shù)具有高的時(shí)間和空間分辨率，可以解決單相機(jī)采樣角度有限及多相機(jī)系統(tǒng)復(fù)雜操作困難等一系列問題，因此基于光場成像技術(shù)的高溫火焰溫度場測量技術(shù)是一種非常有應(yīng)用前景的測溫手段［12-14］。

基于輻射成像的高溫發(fā)光火焰溫度場重構(gòu)問題，從本質(zhì)上是根據(jù)邊界輻射強(qiáng)度分布求解輻射傳輸體系中溫度場的問題，這是一個(gè)典型的輻射反問題。許多學(xué)者開展了有關(guān)重構(gòu)模型與方法的研究。2002年，周懷春等［4］利用改進(jìn)的Tikhonov正則化方法在CCD相機(jī)獲取的火焰圖像中重建出三維溫度場，并證明了改進(jìn)的Tikhonov正則化方法具有強(qiáng)的抗測量誤差能力。2009年，黃群星等［5］采用立體適配器在一個(gè)高分辨率相機(jī)成像面上獲得不同探測角度的火焰圖像，充分利用了相機(jī)的分辨率，采用最小二乘QR（LSQR）分解算法和雙色法實(shí)現(xiàn)了非穩(wěn)態(tài)火焰斷面二維溫度和煙黑顆粒體積分?jǐn)?shù)同時(shí)重建計(jì)算，得出了火焰斷面的重建結(jié)果。2010年，劉冬和岑可法等［15］采用LSQR算法結(jié)合反向蒙特卡羅法實(shí)現(xiàn)二維、三維爐內(nèi)溫度場的反演，基于該方法，進(jìn)一步采用8臺(tái)CCD相機(jī)對300MW 煤粉爐爐膛的火焰溫度分布進(jìn)行了重建，獲得了不同高度截面上溫度分布規(guī)律［16］。2012年，劉冬等［17］在LSQR算法與截?cái)嗥娈愔捣纸馑惴ǎ═SVD）的基礎(chǔ)上，結(jié)合遺傳算法（Genetic Algorithm，GA）構(gòu)造了2種混合優(yōu)化算法，利用GA的隨機(jī)搜索能力減弱非最優(yōu)正則化參數(shù)對于重建精度的影響。2014年，婁春等［18］利用Tikhonov正則化和廣義截?cái)嗥娈愔捣纸獾幕旌纤惴ǎ═R-GSVD）重建了爐內(nèi)的三維溫度分布，并指出在重建過程中需要更新正則化參數(shù)。2016年，周懷春等［10］利用修正的Tikhonov正則化方法重建出了超臨界鍋爐爐膛內(nèi)部的溫度分布。目前，雖然大量學(xué)者針對輻射成像溫度場重構(gòu)問題開展了工作，但仍少有關(guān)于基于光場成像的火焰溫度場重建反問題的報(bào)道。本文主要在光場成像及火焰輻射傳輸模型的基礎(chǔ)上，開展吸收散射性火焰溫度場重建算法研究，將Landweber算法應(yīng)用到基于光場成像的火焰溫度場重建，以具有較高計(jì)算精度的廣義源項(xiàng)多流法作為正問題模型，采用Landweber算法作為反問題求解算法，根據(jù)光場成像模型模擬獲得吸收散射性介質(zhì)的輻射強(qiáng)度光場圖像，據(jù)此重構(gòu)出其三維溫度分布。同時(shí)將在輻射成像測溫中廣泛應(yīng)用的LSQR算法引入到本文研究中，作為對比以檢驗(yàn)Landweber算法的性能。

1 光場成像及火焰輻射傳輸模型

1.1 光場成像模型

對于基于光場成像的高溫火焰溫度場重建問題，首要問題是建立準(zhǔn)確可靠的成像模型以從拍攝得到的光場圖像中提取全場的輻射強(qiáng)度分布。光場相機(jī)的內(nèi)部結(jié)構(gòu)如圖1所示。與傳統(tǒng)相機(jī)不同，在光場相機(jī)中，主透鏡和傳感器之間加裝了微透鏡陣列?；鹧嫔弦稽c(diǎn)發(fā)出的若干采樣光線在經(jīng)過主透鏡之后匯聚于主透鏡像面上的一點(diǎn)，再被微透鏡陣列分散成若干束投射到不同的傳感器像素上。每個(gè)像素上的采樣光線方向信息可由像素和與之對應(yīng)的微透鏡確定，而其強(qiáng)度信息可通過圖像的灰度值確定。

圖1 光場成像的射線追蹤示意圖Fig.1 Schematic diagram of ray tracing of light-field imaging

根據(jù)式（1）和式（2）所示的透鏡成像公式，可從像素點(diǎn)開始對采樣光線進(jìn)行追蹤。

式中：sv為透鏡的物距，s為透鏡的像距，F(xiàn)為透鏡的焦距。對于主透鏡來說，sv為主透鏡與虛擬物面之間的距離，s為主透鏡與虛擬像面之間的距離，xv為光源點(diǎn)的坐標(biāo)，x為對應(yīng)的虛擬像點(diǎn)的坐標(biāo)，F(xiàn)為主透鏡的焦距；對于微透鏡來說，sv為微透鏡陣列與虛擬像面之間的距離，s為微透鏡陣列與探測面之間的距離，xv為虛擬像點(diǎn)的坐標(biāo)，x為與之對應(yīng)的像素點(diǎn)的坐標(biāo)，F(xiàn)為微透鏡的焦距。因此，當(dāng)光場相機(jī)的結(jié)構(gòu)參數(shù)已知時(shí)，可以根據(jù)式（1）和式（2）計(jì)算得到光線的位置信息，即光線經(jīng)過各個(gè)平面的坐標(biāo)。

在獲得各點(diǎn)的坐標(biāo)后，可以根據(jù)各點(diǎn)的連線得到光線的方向信息。對于圖1所示的采樣光線，其在火焰內(nèi)的傳輸方向可以根據(jù)點(diǎn)4和點(diǎn)5的連線得到。

式中：θ和φ分別為天頂角和圓周角。

基于上述光場成像模型，可以從光場相機(jī)拍攝到的光場圖像中獲得火焰邊界上的出射輻射強(qiáng)度的位置、方向和強(qiáng)度信息，為基于光場成像的火焰溫度場重建提供測量信號。

1.2 火焰輻射傳輸模型及算法

在通過光場圖像獲得火焰出射輻射強(qiáng)度信號之后，還需要建立出射輻射強(qiáng)度與火焰內(nèi)部輻射強(qiáng)度之間的聯(lián)系，才能實(shí)現(xiàn)火焰內(nèi)三維溫度場重建?；鹧媸且环N典型的參與性介質(zhì)，邊界上的出射輻射強(qiáng)度是其內(nèi)部各點(diǎn)輻射強(qiáng)度的累積?；鹧鎯?nèi)的輻射傳輸過程可以用如下輻射傳輸方程來描述：

式中：Iλ（ξ，Ω）為介質(zhì)在ξ位置處Ω 方向上的光譜輻射強(qiáng)度，W／（m2·μm·sr）；Ibλ（ξ）為介質(zhì)在ξ位置處的自身黑體光譜輻射強(qiáng)度，W／（m2·μm·sr）；κeλ（ξ）為衰減系數(shù)，m-1；κaλ（ξ）為吸收系數(shù)，m-1；κsλ（ξ）為散射系數(shù)，m-1；Φ（Ω′，Ω）為Ω′方向入射并向Ω 方向散射的散射相函數(shù)。等號右側(cè)第二項(xiàng)和第三項(xiàng)均可以使輻射強(qiáng)度增強(qiáng)，因此將二者之和定義為廣義源項(xiàng)Sλ（ξ，Ω），即

此時(shí)式（5）可以簡化為可以看出，若源項(xiàng)已知，則可將式（5）所示的積分微分方程簡化為式（7）所示的純微分方程，進(jìn)而可以通過沿程積分求解出火焰某一方向上的出射輻射強(qiáng)度。廣義源項(xiàng)多流法就是在這種求解思路的基礎(chǔ)上建立的。強(qiáng)度求解的關(guān)鍵在于對廣義源項(xiàng)的求解。本文利用有限體積法（Finite Volume Method，F(xiàn)VM）計(jì)算火焰的廣義源項(xiàng)。首先采用FVM求解較少離散方向上的輻射強(qiáng)度，從而根據(jù)式（6）獲得火焰內(nèi)的廣義源項(xiàng)分布，根據(jù)廣義源項(xiàng)的不變性［19］，增加了離散角度之后廣義源項(xiàng)的數(shù)值不發(fā)生改變，因此可以利用計(jì)算得到的廣義源項(xiàng)計(jì)算出更多離散方向上的輻射強(qiáng)度。如圖1所示，將Ω方向上探測線經(jīng)過的火焰區(qū)域離散為N個(gè)網(wǎng)格，對式（7）在微元控制體內(nèi)積分，可以得到

式中：κa為介質(zhì)吸收系數(shù)，m-1；κs為介質(zhì)散射系數(shù)，m-1；Ibi為控制體i自身黑體輻射強(qiáng)度，W／（m2·μm·sr）；τi為探測線穿過網(wǎng)格的光學(xué)厚度；為控制體i在離散方向sj上的輻射強(qiáng)度值；ΔΩj為控制立體角；Ii＋1（Ω）為穿過第i個(gè)網(wǎng)格后的輻射強(qiáng)度。因?yàn)榛鹧鏈囟冗h(yuǎn)高于周圍環(huán)境溫度，所以忽略環(huán)境輻射影響。假設(shè)每個(gè)網(wǎng)格內(nèi)溫度均勻分布，則Ω方向上的火焰邊界出射輻射強(qiáng)度可由積分得到，其離散形式可表示為

沿不同的探測方向?qū)鹧鎯?nèi)的輻射傳輸過程進(jìn)行積分，可以得到火焰邊界上的出射輻射強(qiáng)度分布，將其寫成矩陣形式可表示為

式中：Aλ為衰減系數(shù)矩陣；Sλ為廣義源項(xiàng)向量；Iλ為出射輻射強(qiáng)度向量，即光場相機(jī)最終接收到的輻射強(qiáng)度分布。

2 溫度場重建算法

基于光場成像的火焰溫度場重建，即根據(jù)由光場相機(jī)拍攝到的火焰光場圖像獲得火焰邊界上的出射輻射強(qiáng)度分布，進(jìn)而利用該出射輻射強(qiáng)度分布由火焰內(nèi)的輻射傳輸模型重建出火焰的三維溫度場。當(dāng)輻射特性參數(shù)已知時(shí)，式（11）中的未知數(shù)僅為廣義源項(xiàng)Sλ，可以將其視為線性方程組，通過求解式（11），可得到火焰內(nèi)的源項(xiàng)分布。由式（6），火焰的自身黑體輻射強(qiáng)度包含在廣義源項(xiàng)中，當(dāng)源項(xiàng)已知時(shí)，利用FVM 可以計(jì)算出火焰內(nèi)各點(diǎn)的輻射強(qiáng)度，進(jìn)而計(jì)算出火焰的自身黑體輻射強(qiáng)度項(xiàng)，根據(jù)普朗克定律，待重建的溫度場可以由火焰的自身黑體光譜輻射強(qiáng)度Ibλ計(jì)算得到，據(jù)此可實(shí)現(xiàn)火焰三維溫度場的重建。

對于式（11），考慮到其系數(shù)矩陣Aλ是一個(gè)大型的稀疏矩陣，且問題本身具有非適定性，本文利用Landweber算法來求解方程組。對于式（11）所示的方程組，Landweber算法的求解格式如下：式中：α為松弛因子；A*為A的轉(zhuǎn)置矩陣。Landweber算法適用于求解大規(guī)模不適定問題，方法較為穩(wěn)定，即使在測量值存在噪聲的情況下利用該算法也可以得到合理的解。但是傳統(tǒng)的Landweber算法迭代序列收斂速度很慢，為了加速收斂，本文采用了一種改進(jìn)的迭代格式［20］：

式中：a≥2為給定的正整數(shù)；E為單位矩陣；M 為過渡矩陣；n為迭代步數(shù)。雖然式（13）迭代一步的計(jì)算量要大于式（12），但是計(jì)算過程中總的迭代步數(shù)會(huì)大大減少。式（13）迭代k步，相當(dāng)于式（12）迭代ak步，因此該迭代格式會(huì)大大提高算法的效率。在本文中取a＝4。

基于光場成像與Landweber算法的火焰溫度場重建流程如圖2所示。

此外，為了衡量Landweber算法性能的優(yōu)劣，本文還利用LSQR算法重建火焰的三維溫度場以作為比較。該方法已被廣泛應(yīng)用于基于輻射成像的火焰溫度場重建問題中，并被證明具有很高的重建精度?；贚SQR算法的火焰溫度場重建流程與基于Landweber算法相似，在此不做贅述。

圖2 基于光場成像與Landweber算法的火焰溫度場重建流程Fig.2 Flowchart of flame temperature field reconstruction based on light-field imaging and Landweber algorithm

3 數(shù)值模擬結(jié)果與分析

在第1節(jié)成像理論模型的基礎(chǔ)上，開展數(shù)值模擬計(jì)算來檢驗(yàn)算法的性能?？紤]一圓柱型介質(zhì)，介質(zhì)半徑為0.05 m，高度為0.4 m，假設(shè)溫度分布為軸對稱分布。

式中：z和r分別為軸向和徑向坐標(biāo)。假設(shè)火焰的吸收系數(shù)為3m-1，散射系數(shù)為3m-1。在圓周角、半徑和軸向高度方向上將重建區(qū)域劃分為Nφ×NR×NL＝1×20×20個(gè)網(wǎng)格。

為了使模擬的輻射強(qiáng)度信號更貼合實(shí)際情況，在模擬的強(qiáng)度信號基礎(chǔ)上添加隨機(jī)測量噪聲，如下：

式中：為一個(gè)滿足標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)數(shù)；σ為添加噪聲的標(biāo)準(zhǔn)差；Imea和Iexa分別為強(qiáng)度信號的測量值和精確值；γ為測量誤差。

式中：εrel，i為單元網(wǎng)格溫度重建值與真值的相對誤差；為單元網(wǎng)格重建值；為單元網(wǎng)格的真實(shí)溫度值。溫度重建結(jié)果的平均相對誤差為

利用上述光場成像模型模擬構(gòu)建了介質(zhì)邊界輻射強(qiáng)度分布的光場圖像，如圖3所示。模擬研究中采用的光場相機(jī)結(jié)構(gòu)參數(shù)見文獻(xiàn)［21］。

圖3 模擬光場圖像Fig.3 Simulated light-field image

分別利用Landweber算法，LSQR算法對火焰溫度場進(jìn)行重建。為了研究測量誤差對重建精度的影響，分別添加了1%、3%、5%的測量誤差γ。不同測量誤差情況下2種算法的重建結(jié)果如圖4和圖5所示（圖中：R為半徑，Z為高度），縱截面上的重建相對誤差分布如圖6和圖7所示，平均重建相對誤差及計(jì)算時(shí)間對比如表1所示。

從圖4和圖5中可以看出，在添加測量誤差很?。ㄐ∮?%）的情況下，采用2種算法重建出的溫度分布均十分接近于假定的溫度場分布；當(dāng)測量誤差逐漸增大時(shí)，重建結(jié)果中出現(xiàn)了抖動(dòng)，溫度場輪廓逐漸變得不規(guī)則，尤其是在添加了5%測量誤差的情況下，在介質(zhì)頂部位置有了較明顯的波動(dòng)。但從整體上看，重建的溫度分布仍然清晰可見，且與真值相近，因而結(jié)果仍然是合理的。

圖4 LSQR算法三維溫度場重建結(jié)果Fig.4 Reconstructed 3D temperature field distribution using LSQR algorithm

圖5 Landweber算法三維溫度場重建結(jié)果Fig.5 Reconstructed 3D temperature field distribution using Landweber algorithm

圖6 LSQR算法重建相對誤差分布Fig.6 Reconstruction relative error distribution of LSQR algorithm

圖7 Landweber算法重建相對誤差分布Fig.7 Reconstruction relative error distribution of Landweber algorithm

表1 LSQR算法與Landweber算法計(jì)算結(jié)果對比Tab le 1 Com parison of calcu lation resu lts between LSQR algorithm and Landweber algorithm

進(jìn)一步從定量的角度分析重建結(jié)果的優(yōu)劣。從重建相對誤差中可以看出，無論有無測量誤差，利用2種算法都可以得到比較合理的溫度場重建結(jié)果。在沒有測量誤差時(shí)，2種算法溫度場的重建結(jié)果精度很高，平均重建相對誤差均在10-9量級；隨著測量誤差的增大，重建相對誤差也隨之增大；但即使在添加5%的測量誤差的條件下，其溫度場重建相對誤差平均值也分別只有0.91%和0.92%。因此可以看出，2種算法都具有很高的計(jì)算精度，且計(jì)算結(jié)果精度相當(dāng)。

從圖6和圖7中可以看出，火焰溫度場重建相對誤差較大的區(qū)域主要發(fā)生在火焰徑向中心處的軸向兩端位置，尤其是頂端區(qū)域誤差很明顯。對于溫度場軸對稱分布的火焰，由于內(nèi)部強(qiáng)度的沿程衰減作用，其中心位置的網(wǎng)格單元的自身輻射信號到達(dá)相機(jī)過程中衰減得最嚴(yán)重，從而導(dǎo)致信息的缺失；并且由于探測線的分布導(dǎo)致穿過軸向兩端位置的單元的探測線數(shù)量少，從而測量數(shù)據(jù)中攜帶軸向兩端位置單元的信息相對較少。此外，結(jié)合圖4和圖5所示的溫度場真值及圖3所示的光場圖像可以看出，頂端位置處的介質(zhì)單元溫度較低，其本身測量信號強(qiáng)度較弱，在添加了測量噪聲之后，測量信號更容易受到噪聲的干擾而失真，因此在這些位置的介質(zhì)單元的溫度場重建相對誤差較大。

計(jì)算效率是衡量算法優(yōu)劣的另一個(gè)重要指標(biāo)。本文比較了2種算法在重建溫度場時(shí)所消耗的時(shí)間，如表1所示?？梢钥闯?，基于Landweber算法重建溫度場所需時(shí)間均短于LSQR算法所需時(shí)間。當(dāng)2種算法達(dá)到相同的計(jì)算精度時(shí)，利用LSQR算法需要24 s左右的計(jì)算時(shí)間，而利用Landweber算法僅需要2.5 s，計(jì)算效率差距明顯。這主要是由于改進(jìn)的迭代格式使得Landweber算法的迭代步數(shù)大大減少，從而縮短了計(jì)算時(shí)間。事實(shí)上，對于本算例，在存在測量誤差的情況下，Landweber算法僅計(jì)算了10步，而LSQR算法計(jì)算了5 000步，要遠(yuǎn)多于Landweber算法。因此，在計(jì)算效率方面Landweber算法更勝一籌。綜合考慮計(jì)算效率與計(jì)算精度，Landweber算法更適用于基于光場成像的吸收散射性火焰三維溫度場重建問題。

4 結(jié) 論

1）本文提出了基于光場成像技術(shù)的吸收散射性火焰三維溫度場重建算法，構(gòu)建了光場成像模型，采用廣義源項(xiàng)多流法作為正問題計(jì)算方法，將Landweber算法引入到火焰三維溫度場重建問題中，利用Landweber算法計(jì)算廣義源項(xiàng)分布進(jìn)而重建出火焰的三維溫度分布。同時(shí)引入LSQR算法作為對比，以衡量Landweber算法的性能。

2）在理論模型的基礎(chǔ)上開展數(shù)值模擬研究，計(jì)算結(jié)果表明無論有無測量誤差，利用2種算法均可以很好地重建出介質(zhì)的三維溫度場。隨著測量誤差的增大，2種算法的重建相對誤差均有增大，但使在添加5%的測量誤差的條件下，溫度場平均重建相對誤差也分別只有0.91%和0.92%，仍在可接受范圍內(nèi)。

3）對比結(jié)果表明，2種算法具有相當(dāng)?shù)挠?jì)算精度，但Landweber算法的計(jì)算效率要明顯優(yōu)于LSQR算法，在相同重建精度條件下，Landweber算法的計(jì)算時(shí)間為2.5 s，而LSQR算法需要24 s，因而Landweber算法要更適用于火焰三維溫度場重建問題。結(jié)果表明，提出的基于光場成像與Landweber算法的吸收散射性火焰溫度場重建方法是可行的，在輻射成像測溫領(lǐng)域具有很廣闊的應(yīng)用前景。