基于LSMR算法的火焰三維溫度場重建

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

（1 中國計量大學信息工程學院，浙江省電磁波信息技術與計量檢測重點實驗室，浙江 杭州 310018；2 中國計量大學計量測試工程學院，浙江 杭州 310018）

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

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

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

1 測量原理與方法

1.1 光場成像

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

圖1 光場相機成像

1.2 火焰輻射光場成像模型

將通過光場相機采集到的火焰輻射信息和重建算法相結合，進行火焰三維溫度場的重建?；鹧孑椛涔鈭鰣D像的輻射強度計算采用的是穩(wěn)態(tài)傳輸方程，為簡化過程，僅考慮火焰的吸收特性[17]，火焰的穩(wěn)態(tài)傳輸方程可以表示為式(1)。

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

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

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

1.3 LSMR算法

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

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

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

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

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

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

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

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

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

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

2 數(shù)值模擬實驗結果與分析

2.1 實驗條件

為了對LSMR算法的重建性能進行研究，基于火焰輻射光場成像模型，進行了相關的數(shù)值模擬實驗?；鹧嫘螤钤O定為底面半徑R=20mm 和高度Z=40mm 的圓柱體，吸收系數(shù)設置為10m-1。火焰分布模型如式(12)和式(13)所示[17]。式中T為火焰的溫度大小，K；z為火焰的軸向坐標；r為火焰的徑向坐標。其溫度分布滿足旋轉對稱分布，該模擬火焰符合火焰的溫度分布特點?；鹧鏈囟确植既鐖D2所示，溫度范圍為900～2100K。

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

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

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

2.2 火焰輻射強度的非負性分析

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

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

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

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

圖3 中心層內焰區(qū)域輻射強度重建對比

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

2.3 火焰輻射強度的誤差分析

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

圖4 LSQR、NNLS 和LSMR算法計算火焰輻射強度的平均相對誤差

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

2.4 火焰輻射強度的計算效率分析

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

表2 不同噪聲下計算時間的對比

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

2.5 溫度場重建結果分析

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

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

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

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

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

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

3 結論

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

（1）使用LSMR 方法對火焰輻射強度求解時，不會出現(xiàn)負值。因此，LSMR 求解的精度和NNLS相當，比LSQR有所提高。

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

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