單良，仰文淇，洪波，周榮幸，孔明

（1 中國(guó)計(jì)量大學(xué)信息工程學(xué)院，浙江省電磁波信息技術(shù)與計(jì)量檢測(cè)重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 310018；2 中國(guó)計(jì)量大學(xué)計(jì)量測(cè)試工程學(xué)院，浙江 杭州 310018）

燃燒和火焰作為自然界和人類生產(chǎn)、生活中普遍存在的現(xiàn)象，在內(nèi)燃機(jī)、電廠鍋爐、火箭推進(jìn)器、燃?xì)廨啓C(jī)和化學(xué)工程等工業(yè)生產(chǎn)領(lǐng)域扮演著重要角色[1-2]。及時(shí)并且準(zhǔn)確地獲得火焰溫度分布，可以更好地探究燃燒的本質(zhì)和規(guī)律，并在此基礎(chǔ)上提高燃燒效率，優(yōu)化燃燒過程[3]。目前圖像探測(cè)器的分辨率可以達(dá)到很高的像素級(jí)別，高集成度、高分辨率、高性能、低成本的圖像傳感器的迅速發(fā)展，推動(dòng)了火焰輻射成像測(cè)溫技術(shù)的發(fā)展?；鹧孑椛涑上駵y(cè)溫技術(shù)具有測(cè)量精度高、測(cè)量實(shí)時(shí)性好、測(cè)量范圍廣等優(yōu)點(diǎn)，逐漸成為目前火焰溫度場(chǎng)測(cè)量技術(shù)領(lǐng)域的研究熱點(diǎn)[4-5]。

閆勇等[6]使用單個(gè)CCD相機(jī)，結(jié)合圖像處理和雙色輻射測(cè)溫技術(shù)，獲得了燃燒系統(tǒng)中火焰的三維溫度信息。這里使用單個(gè)CCD相機(jī)拍攝火焰圖像，獲取的火焰信息有限，導(dǎo)致溫度場(chǎng)重建的效果會(huì)受到一定影響。如果要提高火焰溫度場(chǎng)重建的準(zhǔn)確性和效果，則需要使用多臺(tái)相機(jī)從不同的角度拍攝火焰圖像。周懷春等[7]使用多個(gè)CCD相機(jī)獲取火焰圖像，引入改進(jìn)的Tikhonov正則化方法重建出三維溫度場(chǎng)。岑可法等[8]基于多CCD相機(jī)系統(tǒng)，利用代數(shù)重建方法（ART）測(cè)量火焰斷面溫度場(chǎng)。這些團(tuán)隊(duì)系統(tǒng)研究了求解火焰輻射反問題的算法，并分析了這些方法的重建特性。但是多相機(jī)系統(tǒng)的復(fù)雜性較大，對(duì)不同空間位置的相機(jī)進(jìn)行耦合同步比較困難。光場(chǎng)相機(jī)可以通過單臺(tái)相機(jī)在一次曝光中采集四維光場(chǎng)信息[9]。此特點(diǎn)可以解決輻射測(cè)溫多相機(jī)系統(tǒng)光路復(fù)雜、同步觸發(fā)難等問題，在輻射成像三維溫度重建時(shí)有其獨(dú)特優(yōu)勢(shì)[10-15]。許傳龍等[16]將光場(chǎng)成像技術(shù)與火焰溫度重建結(jié)合，然后用LSQR對(duì)該數(shù)學(xué)模型進(jìn)行求解，重構(gòu)出火焰的溫度分布。LSQR 方法是一種迭代正則化方法，適用于求解大規(guī)模不適定問題。但是在求解時(shí)無(wú)法保證結(jié)果的非負(fù)性[17]，而且在重建的準(zhǔn)確性和效率方面仍有提升空間。孫俊等[17]提出了NNLS-LMBC 混合算法，可同時(shí)重建出火焰的溫度場(chǎng)和吸收系數(shù)，NNLS 算法能夠保證求解結(jié)果的非負(fù)性，且具有較好的計(jì)算穩(wěn)定性。但是計(jì)算效率低，計(jì)算時(shí)間長(zhǎng)。

LSMR 算法是一種用于求解稀疏最小二乘問題的迭代算法，其原理是將極小殘差法用來(lái)求解線性方程組和最小平方問題，已經(jīng)被應(yīng)用于圖像處理等方面的研究，證明了LSMR 算法在解決不適定問題上具有良好的性能[18-20]。針對(duì)LSQR 算法在火焰微元體輻射強(qiáng)度求解中出現(xiàn)負(fù)值和NNLS算法的計(jì)算效率低等問題，本文提出將LSMR 算法用于火焰光場(chǎng)成像三維溫度場(chǎng)重建。構(gòu)建單光場(chǎng)相機(jī)的輻射成像模型，從輻射強(qiáng)度重建結(jié)果的非負(fù)性、求解精度和求解效率三個(gè)方面分析LSMR 算法的求解性能。然后使用LSMR 算法對(duì)模擬光場(chǎng)火焰進(jìn)行溫度場(chǎng)重建，驗(yàn)證其進(jìn)行火焰三維溫度場(chǎng)重建的準(zhǔn)確性。

1 測(cè)量原理與方法

1.1 光場(chǎng)成像

光場(chǎng)相機(jī)由主透鏡、微透鏡陣列和CCD 陣列組成。傳統(tǒng)一代光場(chǎng)相機(jī)的成像如圖1所示。光源點(diǎn)發(fā)射出不同方向的輻射光線，經(jīng)過主透鏡后到達(dá)微透鏡，微透鏡根據(jù)方向分離這些光線，并在微透鏡后面的CCD 陣列上記錄下輻射強(qiáng)度信息，同時(shí)光線的位置和方向信息可以通過光場(chǎng)相機(jī)中CCD陣列和微透鏡陣列來(lái)確定[21]。

圖1 光場(chǎng)相機(jī)成像

1.2 火焰輻射光場(chǎng)成像模型

將通過光場(chǎng)相機(jī)采集到的火焰輻射信息和重建算法相結(jié)合，進(jìn)行火焰三維溫度場(chǎng)的重建。火焰輻射光場(chǎng)圖像的輻射強(qiáng)度計(jì)算采用的是穩(wěn)態(tài)傳輸方程，為簡(jiǎn)化過程，僅考慮火焰的吸收特性[17]，火焰的穩(wěn)態(tài)傳輸方程可以表示為式(1)。

式中，ILλ為火焰輻射光線L 的光譜輻射強(qiáng)度；κL為光線L 穿過火焰微元體的吸收系數(shù)；lL為該微元體的幾何長(zhǎng)度；Ibλ為火焰微元體的黑體輻射強(qiáng)度，W/(m3·sr)，根據(jù)普朗克定律，Ibλ表示為式(2)。

式中，T為微元體的溫度，K；c1為第一輻射常數(shù)，c1=3.7418×10-16W·m2；c2為第二輻射常數(shù)，c2=1.4388×10-2m·K；λ為火焰輻射光線的波長(zhǎng)。

將火焰劃分為V個(gè)微元體，依次編號(hào)為(1,2, …,v, …,V)，光場(chǎng)相機(jī)的圖像探測(cè)器共有M個(gè)像素接收到火焰輻射光線，依次編號(hào)為(1, 2, …,m, …,M)。像素m的火焰輻射光線穿過n個(gè)微元體，沿著光線的傳輸方向依次編號(hào)為(v1, v2, …,vi, …, vj, …, vn)，將式(1)離散化得到式(3)、式(4)。

1.3 LSMR算法

LSMR迭代算法可以用于求解線性系統(tǒng)Ax-b和最小平方問題min‖Ax-b‖2。LSMR 算法基于Golub-Kahan 雙對(duì)角化過程，等同于將最小殘差法應(yīng)用于正規(guī)方程ATAx=ATb。

對(duì)于光場(chǎng)相機(jī)圖像傳感器各像素對(duì)應(yīng)的火焰輻射光線的輻射強(qiáng)度向量Iλ，通過構(gòu)造系數(shù)矩陣A，可建立線性方程Iλ=AIbλ。利用LSMR 算法求解線性方程Iλ=AIbλ的過程等價(jià)于求解線性系統(tǒng)Ax-b的過程。

對(duì)于線性系統(tǒng)Ax-b，LSMR算法的具體實(shí)現(xiàn)是通過Golub-Kahan過程實(shí)現(xiàn)雙對(duì)角化，即式(6)。

選擇合適的標(biāo)量αk≥0和βk≥0，使得uk和vk的范數(shù)等于1。式(6)經(jīng)過k次迭代可得式(7)。

e表示單位矩陣；ek+1表示單位矩陣中第k+1 列向量。結(jié)合式(7)可以得到式(8)。

對(duì)于正規(guī)方程ATAx=ATb，rk=b-Axk為其每一步的迭代殘差，每一步迭代過程中通過選擇yk尋求方程的近似解，如式(9)所示。

LSMR算法在每一步迭代中選擇yk，使‖ATrk‖達(dá)到極小值，如式(10)所示。

而LSQR 算法的每一步迭代中選擇yk，只使得‖rk‖最小化，如式(11)所示。

在線性系統(tǒng)的求解過程中，LSQR 算法[22]和LSMR 算法都可以得到最小范數(shù)解。在LSQR 算法的迭代過程中，只有殘差范數(shù)‖rk‖單調(diào)遞減；但是在LSMR算法中，殘差范數(shù)‖rk‖和‖ATrk‖都是單調(diào)遞減的。LSMR 算法求解得到的結(jié)果比LSQR算法更加精確和穩(wěn)定。

相比于LSQR算法，LSMR算法求解速度更快，得到的結(jié)果更加安全。NNLS 算法[23]可以保證求解的非負(fù)性，但是其算法中存在矩陣求逆運(yùn)算，而LSQR 和LSMR 算法沒有這一過程，計(jì)算效率會(huì)顯著提高。

2 數(shù)值模擬實(shí)驗(yàn)結(jié)果與分析

2.1 實(shí)驗(yàn)條件

為了對(duì)LSMR算法的重建性能進(jìn)行研究，基于火焰輻射光場(chǎng)成像模型，進(jìn)行了相關(guān)的數(shù)值模擬實(shí)驗(yàn)?；鹧嫘螤钤O(shè)定為底面半徑R=20mm 和高度Z=40mm 的圓柱體，吸收系數(shù)設(shè)置為10m-1?；鹧娣植寄Ｐ腿缡?12)和式(13)所示[17]。式中T為火焰的溫度大小，K；z為火焰的軸向坐標(biāo)；r為火焰的徑向坐標(biāo)。其溫度分布滿足旋轉(zhuǎn)對(duì)稱分布，該模擬火焰符合火焰的溫度分布特點(diǎn)。火焰溫度分布如圖2所示，溫度范圍為900～2100K。

圖2 火焰數(shù)值模擬模型

以0mm 的位置為火焰的中心層，每隔5mm 取一層火焰溫度數(shù)據(jù)，將被測(cè)對(duì)象火焰分成7層，每層的大小為72×72。分別用LSQR、NNLS 和LSMR三種算法對(duì)火焰的輻射強(qiáng)度進(jìn)行重建，為了定量分析不同算法的重建效果，在原始火焰輻射強(qiáng)度分布數(shù)據(jù)上加入1%、5%、10%、15%和20%高斯白噪聲，再將重建結(jié)果與相對(duì)應(yīng)的原始數(shù)據(jù)進(jìn)行對(duì)比和分析。其中絕對(duì)誤差和相對(duì)誤差的定義如式(14)和式(15)所示。

式中，ΔIabs為重建的絕對(duì)誤差；ΔIrel為重建結(jié)果的相對(duì)誤差；I為使用算法重建出來(lái)的值；Iset為根據(jù)式(12)設(shè)定的火焰數(shù)據(jù)。

2.2 火焰輻射強(qiáng)度的非負(fù)性分析

在設(shè)定的火焰輻射強(qiáng)度分布中疊加1%、5%、10%、15%、20%的噪聲，再用LSMR、LSQR 和NNLS 三種算法分別求解火焰的輻射強(qiáng)度。統(tǒng)計(jì)7層中所有重建結(jié)果中出現(xiàn)負(fù)值的數(shù)量并計(jì)算所占的百分比，結(jié)果見表1。

表1 輻射強(qiáng)度重建時(shí)出現(xiàn)負(fù)值點(diǎn)的數(shù)量及百分比

如表1所示，在噪聲水平較低的情況下，例如1%噪聲時(shí)，LSQR、NNLS 和LSMR 三種算法對(duì)火焰微元體輻射強(qiáng)度的重建都未出現(xiàn)負(fù)值。但是隨著噪聲的增大，LSQR 算法無(wú)法保證求解結(jié)果的非負(fù)性，對(duì)火焰微元體輻射強(qiáng)度求解的結(jié)果中會(huì)出現(xiàn)越來(lái)越多的負(fù)值。相比之下，LSMR 和NNLS 求解的結(jié)果中未出現(xiàn)負(fù)值，保證了求解的非負(fù)性。

LSQR 算法對(duì)輻射強(qiáng)度的求解會(huì)收斂至負(fù)值，主要集中在火焰的內(nèi)焰部分。用LSQR、NNLS 和LSMR三種算法在不同噪聲水平下對(duì)火焰的輻射強(qiáng)度進(jìn)行重建后，選取中心切片層大小為21×21的內(nèi)焰區(qū)域，對(duì)比和分析此區(qū)域的重建效果，如圖3所示。

圖3 中心層內(nèi)焰區(qū)域輻射強(qiáng)度重建對(duì)比

圖3(a)為該區(qū)域的原始輻射強(qiáng)度分布，圖3(b)、(c)、(d)分別為L(zhǎng)SQR、NNLS 和LSMR 算法的重建結(jié)果。可以看出，LSQR 算法在1%噪聲時(shí)還呈現(xiàn)出較好的重建效果，當(dāng)噪聲增加到5%時(shí)，出現(xiàn)了少量的負(fù)值點(diǎn)。隨著噪聲的不斷增大，其求解的結(jié)果中出現(xiàn)越來(lái)越多的負(fù)值解。在5%、10%、15%、20%情況下此區(qū)域出現(xiàn)負(fù)值解的百分比分別為4.3%、14.7%、28.6%、34.7%，中心層內(nèi)焰區(qū)域負(fù)值解的百分比遠(yuǎn)高于所有層全部解的情況。LSMR 和NNLS 算法始終保持著求解結(jié)果的非負(fù)性。正因?yàn)榇?，在噪聲逐漸增大的情況下，LSQR 算法重建的偏差比LSMR 和NNLS 要顯著增大。LSMR 和NNLS 算法比LSQR 算法具有更好的求解穩(wěn)定性。

2.3 火焰輻射強(qiáng)度的誤差分析

計(jì)算LSMR、LSQR和NNLS三種算法在不同噪聲水平下對(duì)火焰輻射強(qiáng)度重建的平均相對(duì)誤差，如圖4所示。

圖4 LSQR、NNLS 和LSMR算法計(jì)算火焰輻射強(qiáng)度的平均相對(duì)誤差

由圖4 可以看出，在噪聲很低的情況下，LSQR、NNLS 和LSMR 三種算法的平均相對(duì)誤差十分接近，這是因?yàn)長(zhǎng)SQR算法求解的結(jié)果中還未出現(xiàn)負(fù)值解。隨著噪聲水平的增加，三種算法的平均相對(duì)誤差都在增大，但是由于LSQR算法求解結(jié)果中出現(xiàn)的負(fù)值解變多，平均相對(duì)誤差也一直高于LSMR 和NNLS 算法。在噪聲水平為5%、10%、15%、20%的情況下，LSMR算法重建的平均相對(duì)誤差比LSQR 算法分別要低11.4%、14.5%、13.3%。LSMR和NNLS算法的平均相對(duì)誤差都很接近，展現(xiàn)出比LSQR算法更好的抗噪性能。

2.4 火焰輻射強(qiáng)度的計(jì)算效率分析

計(jì)算效率是評(píng)估算法性能的另一個(gè)重要指標(biāo)。在噪聲水平為1%、5%、10%、15%、20%的情況下，分別統(tǒng)計(jì)LSQR、NNLS和LSMR三種算法多次計(jì)算的時(shí)間，并求得平均計(jì)算時(shí)間進(jìn)行對(duì)比，如表2所示。

表2 不同噪聲下計(jì)算時(shí)間的對(duì)比

在計(jì)算中，系數(shù)矩陣會(huì)因?yàn)樗砑拥脑肼暥a(chǎn)生偏差，使得方程組的病態(tài)程度加劇。所以在不同的噪聲水平下，方程組的計(jì)算復(fù)雜度也不一樣，計(jì)算的時(shí)間也會(huì)受到影響。隨著噪聲的增加，三種算法的求解時(shí)間都隨著計(jì)算復(fù)雜度的變化而發(fā)生波動(dòng)。如表2 所示，得到LSQR 的平均用時(shí)為2.90～3.05s，NNLS的平均用時(shí)為2201～2278s，LSMR的計(jì)算時(shí)間僅為0.15s 左右。可以看出，LSMR 算法的計(jì)算時(shí)間比NNLS 算法降低了4 個(gè)數(shù)量級(jí)。同時(shí)LSMR 算法的計(jì)算時(shí)間只有LSQR 算法用時(shí)的二十分之一，LSMR算法展現(xiàn)出在不同噪聲水平下進(jìn)行計(jì)算時(shí)高效的求解效率。

2.5 溫度場(chǎng)重建結(jié)果分析

為了驗(yàn)證LSMR 算法對(duì)火焰溫度場(chǎng)的重建效果，分別疊加1%、5%、10%、15%、20%的噪聲，重建后的溫度分布和絕對(duì)誤差分布如圖5 和圖6所示。

圖5 LSMR算法在不同噪聲下溫度的重建效果

圖6 LSMR算法在不同噪聲下溫度重建的絕對(duì)誤差分布

如圖5 所示，在不同的噪聲水平下，LSMR 算法重建出來(lái)的火焰溫度分布均可以正確反映火焰的溫度趨勢(shì)特征。如圖6所示，重建的誤差主要分布在火焰的高溫區(qū)，在1%、5%、10%、15%、20%噪聲水平下，LSMR 重建的最大誤差分別為4K、29K、68K、141K、215K，對(duì)于最高溫度在2100K的高溫區(qū)域，在合理的誤差范圍之內(nèi)。

LSMR算法對(duì)火焰溫度場(chǎng)進(jìn)行重建時(shí)的平均相對(duì)誤差和最大相對(duì)誤差見表3，兩者都隨著疊加噪聲的增加而增大，但是在不同的噪聲水平下，平均相對(duì)誤差和最大相對(duì)誤差都保持在一個(gè)較低的水平。LSMR 算法能夠準(zhǔn)確地重建出火焰的三維溫度場(chǎng)，表明了提出的LSMR算法是精確、合理、可靠的。

表3 LSMR算法溫度重建的誤差

3 結(jié)論

基于光場(chǎng)成像技術(shù)，采用LSMR方法對(duì)火焰的三維溫度場(chǎng)重建進(jìn)行了數(shù)值模擬研究，并與LSQR和NNLS方法對(duì)火焰輻射強(qiáng)度重建的效果進(jìn)行了比較，主要結(jié)論如下。

（1）使用LSMR 方法對(duì)火焰輻射強(qiáng)度求解時(shí)，不會(huì)出現(xiàn)負(fù)值。因此，LSMR 求解的精度和NNLS相當(dāng)，比LSQR有所提高。

（2）在計(jì)算時(shí)間方面，LSMR 的平均計(jì)算時(shí)間僅為0.15s 左右，低于LSQR 方法的2.90～3.05s 和NNLS方法的2201～2278s。

（3）在不同的噪聲水平下，LSMR 都能夠完成對(duì)火焰三維溫度場(chǎng)的重建。并且在噪聲水平較高的情況下，溫度重建的平均相對(duì)誤差能保持在1.2%以內(nèi)。