基于阻尼LSQR-LMBC 的火焰三維溫度場(chǎng)重建*

單良 趙騰飛 黃薈云 洪波 孔明

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

2) (中國(guó)計(jì)量大學(xué)計(jì)量測(cè)試工程學(xué)院，杭州 310018)

光場(chǎng)相機(jī)可以解決輻射測(cè)溫多相機(jī)系統(tǒng)光路復(fù)雜、同步觸發(fā)難等問(wèn)題，在輻射成像三維溫度重建時(shí)有其獨(dú)特優(yōu)勢(shì).LSQR 是求解基于大型稀疏矩陣最小二乘問(wèn)題的經(jīng)典算法，該算法用于重建三維溫度場(chǎng)時(shí)對(duì)溫度初值依賴(lài)較大，在信噪比較低的情況下重建精度不理想.本文提出阻尼LSQR-LMBC 重建算法，通過(guò)在LSQR方法中添加阻尼正則化項(xiàng)，提高火焰三維溫度場(chǎng)重建的抗噪性能，并結(jié)合LMBC 算法，實(shí)現(xiàn)吸收系數(shù)和三維溫度場(chǎng)同時(shí)求解.在數(shù)值模擬部分，隨著信噪比逐漸降低，阻尼LSQR 的重建效果比LSQR 更加穩(wěn)定，在信噪比達(dá)到13.86 dB 時(shí)，重建精度大約提高30%.阻尼LSQR-LMBC 的平均重建誤差為6.63%.用丁烷火焰進(jìn)行了實(shí)驗(yàn)，重建的丁烷火焰三維溫度場(chǎng)分布符合輻射火焰燃燒的特征，和熱電偶的測(cè)溫?cái)?shù)據(jù)結(jié)果進(jìn)行對(duì)比，相對(duì)誤差在6.8%左右.

1 引言

溫度是火焰燃燒評(píng)價(jià)的重要參數(shù)之一，溫度的研究有助于探究燃燒的本質(zhì)和規(guī)律，推進(jìn)燃燒理論和技術(shù)的發(fā)展[1].對(duì)于三維溫度場(chǎng)的測(cè)量，接觸式測(cè)溫技術(shù)會(huì)破壞火焰的原本形狀和溫度分布，要得到整個(gè)燃燒場(chǎng)的溫度信息，就必須合理布置測(cè)量點(diǎn)，通過(guò)估算得到溫度場(chǎng)的近似分布，這樣會(huì)大大增加測(cè)量系統(tǒng)的復(fù)雜度[2].非接觸式測(cè)溫最顯著的特點(diǎn)就是不對(duì)流場(chǎng)造成擾動(dòng)，具有測(cè)溫范圍寬，動(dòng)態(tài)響應(yīng)靈敏等優(yōu)點(diǎn).目前的光電探測(cè)元件技術(shù)發(fā)展相當(dāng)成熟，CCD 探測(cè)器可以非常好地采集火焰輻射強(qiáng)度信息，反演火焰的二維和三維溫度場(chǎng)信息，具有非接觸、準(zhǔn)確度高、實(shí)時(shí)性好等優(yōu)點(diǎn)，受到國(guó)內(nèi)外科研人員的青睞.

王式民等[3]結(jié)合單CCD 相機(jī)與光學(xué)分層成像法對(duì)火焰的三維溫度場(chǎng)進(jìn)行重建，其通過(guò)調(diào)整相機(jī)采集不同聚焦平面的二維圖像進(jìn)行疊加，最終得到燭光火焰的三維溫度場(chǎng)，但是在調(diào)焦過(guò)程中出現(xiàn)的延時(shí)誤差無(wú)法實(shí)現(xiàn)瞬態(tài)火焰的溫度重建.黃群星等[4]將截?cái)嗥娈愔捣纸?Truncate singular value decomposition，TSVD)的正則化方法應(yīng)用于三維火焰溫度場(chǎng)重建中，較好地解決了重建問(wèn)題的不適定性.Gilabert 等[5]提出一種采用多CCD 結(jié)合光學(xué)層析技術(shù)對(duì)火焰進(jìn)行三維溫度重建的方法.多相機(jī)系統(tǒng)雖然可以獲取火焰不同方向的火焰輻射信息，但是相機(jī)系統(tǒng)的復(fù)雜度會(huì)隨著CCD 數(shù)量的增加而增加，同時(shí)也會(huì)帶來(lái)相機(jī)的同步問(wèn)題.Hossain等[6]將傳統(tǒng)CCD 相機(jī)采用基于投影的CT 層析方法對(duì)火焰的三維溫度場(chǎng)進(jìn)行重建，但是層析法實(shí)現(xiàn)三維溫度場(chǎng)測(cè)量時(shí)存在需求設(shè)備多、光路復(fù)雜、耗時(shí)長(zhǎng)等常見(jiàn)弊端.陸永剛等[7]通過(guò)電動(dòng)機(jī)改變傳統(tǒng)相機(jī)的拍攝位置對(duì)火焰進(jìn)行移焦拍攝，再通過(guò)光學(xué)分層技術(shù)實(shí)現(xiàn)了火焰的三維重建，但是火焰的形狀是動(dòng)態(tài)變化的，移焦拍攝會(huì)出現(xiàn)較大的延時(shí)誤差，難以采集到同一時(shí)刻下的火焰輻射圖像序列.

光場(chǎng)相機(jī)作為一種新型的圖像采集設(shè)備，逐漸在三維重建上得到應(yīng)用和發(fā)展[8-14].光場(chǎng)相機(jī)是利用微透鏡陣列得到圖像，可以獲得比傳統(tǒng)相機(jī)更豐富的四維數(shù)據(jù)信息(光線(xiàn)的強(qiáng)度信息和方向信息)，再結(jié)合成像理論和重建算法實(shí)現(xiàn)物體的三維重建.Fahringer 等[15]提出了一種基于單個(gè)全光相機(jī)的三維三分量粒子圖像測(cè)速技術(shù).Horstmeyer 等[16]利用光場(chǎng)架構(gòu)，將濾光片與光場(chǎng)相機(jī)相結(jié)合構(gòu)建出一種多模式成像方案，可以以數(shù)字的方式對(duì)特定光譜和偏振光等信息進(jìn)行重建.Prevedel 等[17]將光場(chǎng)成像技術(shù)引入醫(yī)學(xué)神經(jīng)系統(tǒng)的研究中，并成功獲得了被研究者腦中單個(gè)神經(jīng)元活動(dòng)圖像.聶云峰等[18，19]對(duì)光場(chǎng)成像技術(shù)的進(jìn)展進(jìn)行了探索，根據(jù)光場(chǎng)成像原理研制了光場(chǎng)成像系統(tǒng)，開(kāi)展實(shí)驗(yàn)對(duì)光場(chǎng)相機(jī)的數(shù)字對(duì)焦能力進(jìn)行驗(yàn)證，并對(duì)離焦對(duì)象進(jìn)行了光譜復(fù)原.也有一些研究團(tuán)隊(duì)將光場(chǎng)成像技術(shù)應(yīng)用于火焰的三維溫度場(chǎng)重建.例如，Yuan 等[20，21]對(duì)火焰輻射光場(chǎng)的成像過(guò)程開(kāi)展了深入的研究，并分析和評(píng)價(jià)了火焰圖像的成像特征，奠定了火焰光場(chǎng)成像三維溫度場(chǎng)重建的理論基礎(chǔ).東南大學(xué)的許傳龍團(tuán)隊(duì)[22-26]將火焰?zhèn)鬏斈Ｐ团c光場(chǎng)成像技術(shù)相結(jié)合，對(duì)半透明介質(zhì)的火焰光場(chǎng)成像特征展開(kāi)了一系列的研究，并構(gòu)建了基于光場(chǎng)相機(jī)的火焰成像系統(tǒng)，對(duì)火焰進(jìn)行了三維溫度場(chǎng)的重建.

本文在光場(chǎng)成像的基礎(chǔ)上，將阻尼LSQR(Damped least squares QR decomposition，LSQR)和LMBC (Levenberg-Marquardt with boundary constraint，LMBC)算法相結(jié)合，在火焰溫度場(chǎng)的反演重建算法上進(jìn)行了新的探索與嘗試.利用雙平面參數(shù)法對(duì)光線(xiàn)進(jìn)行追跡，通過(guò)光場(chǎng)重建策略進(jìn)行輻射光場(chǎng)的正向模擬實(shí)驗(yàn)，使用阻尼LSQRLMBC 的混合算法對(duì)模擬光場(chǎng)圖像進(jìn)行反向溫度場(chǎng)重建.對(duì)丁烷火焰的三維溫度場(chǎng)進(jìn)行重建實(shí)驗(yàn)，并使用熱電偶對(duì)重建結(jié)果進(jìn)行對(duì)比分析.

2 原理與方法

2.1 光場(chǎng)成像模型

本文采用的傳統(tǒng)一代光場(chǎng)相機(jī)，成像示意如圖1 所示，其主要由主透鏡、微透鏡陣列以及CCD圖像傳感器組成.光源點(diǎn)發(fā)射出不同方向的光線(xiàn)，經(jīng)過(guò)主透鏡后再穿過(guò)微透鏡陣列投射到CCD 面成像，同一物點(diǎn)出發(fā)的不同方向的光線(xiàn)投射到同一微透鏡下的不同像素點(diǎn)，因此光場(chǎng)圖像傳感器的像素點(diǎn)獲取的信息包含方向信息[27].

圖1 光場(chǎng)相機(jī)結(jié)構(gòu)示意圖Fig.1.Structure diagram of light field camera.

本文將三維火焰中的某一個(gè)體積很小的區(qū)域定義為虛擬微元體.在對(duì)火焰進(jìn)行三維溫度場(chǎng)重建時(shí)，由于虛擬微元體的體積很小，因此可以將其視為點(diǎn)光源，再根據(jù)雙平面參數(shù)法對(duì)光線(xiàn)進(jìn)行追跡定位.根據(jù)光場(chǎng)相機(jī)的內(nèi)部光學(xué)器件的結(jié)構(gòu)順序可知，微元體發(fā)射出的光線(xiàn)到達(dá)光場(chǎng)相機(jī)圖像傳感器面的追跡過(guò)程及對(duì)應(yīng)的光線(xiàn)追跡模型如下.

1)從微元體出發(fā)的光線(xiàn)經(jīng)過(guò)自由空間到達(dá)主透鏡平面:

2)光線(xiàn)經(jīng)過(guò)主透鏡:

3)透過(guò)主透鏡的光線(xiàn)經(jīng)過(guò)自由空間到達(dá)微透鏡平面:

4)光線(xiàn)經(jīng)過(guò)微透鏡:

5)透過(guò)微透鏡的光線(xiàn)經(jīng)過(guò)自由空間，最終投射在圖像傳感器面:

式中，x和y分別表示光線(xiàn)穿過(guò)某個(gè)平面時(shí)對(duì)應(yīng)的光線(xiàn)交點(diǎn)位置坐標(biāo)，下標(biāo)對(duì)應(yīng)不同的平面;s0和si分別對(duì)應(yīng)物距和像距;fm和fi分別表示主透鏡和微透鏡的焦距;sx和sy表示微透鏡中心與Z軸的偏移距離.對(duì)θ和φ的數(shù)值可以利用正余弦定理，通過(guò)光場(chǎng)相機(jī)內(nèi)部元器件的參數(shù)和器件之間的距離計(jì)算得到.

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

火焰的輻射強(qiáng)度值通過(guò)圖像傳感器接收到的光線(xiàn)輻照強(qiáng)度進(jìn)行分析，采用穩(wěn)態(tài)傳輸方程進(jìn)行計(jì)算[28].在較短的曝光時(shí)間下可以忽略火焰動(dòng)態(tài)變化，直接考慮傳輸過(guò)程中火焰輻射能量對(duì)于介質(zhì)的散射、發(fā)射和吸收特性.根據(jù)Mie 理論，由于炭黑顆粒的吸收能力遠(yuǎn)大于散射能力，因此可以忽略火焰的散射特性，僅考慮火焰的吸收特性，火焰的穩(wěn)態(tài)傳輸方程可以表示為

式中，T為微元體的溫度，單位為K;c1和c2為輻射常數(shù)，分別為3.7418×10—16W·m2和1.4388×10—2m·K;λ為火焰輻射光線(xiàn)的波長(zhǎng).

在進(jìn)行火焰的三維重建時(shí)，需要將被測(cè)對(duì)象劃分為W個(gè)微元體，并依次編號(hào)為 (1，2，w，···，W)，光場(chǎng)相機(jī)的圖像探測(cè)器一共有M個(gè)子像素，依次編號(hào)為 (1，2，m，···，M).像素m的火焰輻射光線(xiàn)穿過(guò)n個(gè)微元體，沿著光線(xiàn)的傳輸方向依次編號(hào)為(w1，w2，wi，wj，···，wn)，將 (6) 式離散化得到 (8)式和(9) 式:

所有輻射光線(xiàn)的光譜輻射強(qiáng)度計(jì)算用 (10) 式的矩陣形式表示:

2.3 阻尼LSQR-LMBC 算法

在已知吸收系數(shù)的條件下，光場(chǎng)火焰圖像的三維溫度場(chǎng)重建問(wèn)題屬于線(xiàn)性反問(wèn)題.然而，在進(jìn)行實(shí)際火焰的測(cè)量時(shí)，火焰的吸收系數(shù)等物性參數(shù)是未知的，溫度重建為線(xiàn)性問(wèn)題和非線(xiàn)性問(wèn)題的混合問(wèn)題.為對(duì)實(shí)際火焰進(jìn)行溫度場(chǎng)的計(jì)算，很多研究人員將非線(xiàn)性算法與線(xiàn)性算法相結(jié)合，如孫俊[29]提出非負(fù)最小二乘(Non-negative least squares，NNLS-LMBC)的混合算法，孫俊陽(yáng)[30]提出Tikhonov-LMBC 的混合算法，這兩種算法都可以在火焰三維溫度場(chǎng)的重建過(guò)程中同時(shí)計(jì)算溫度和吸收系數(shù).本文使用阻尼LSQR-LMBC 混合算法對(duì)吸收系數(shù)未知的火焰進(jìn)行三維溫度場(chǎng)重建.

阻尼LSQR 是在LSQR 的基礎(chǔ)上添加n個(gè)阻尼方程[31]，目的是為了解決LSQR 算法對(duì)初始值的依賴(lài).(10) 式所示的溫度重建線(xiàn)性問(wèn)題的目標(biāo)函數(shù)如式(11)所示:

其中，m為像素點(diǎn)的個(gè)數(shù)，n為火焰劃分后的微元體的數(shù)目，I為單位矩陣，λ表示阻尼因子，是一個(gè)實(shí)常數(shù).再轉(zhuǎn)換為 (12) 式，結(jié)合Lanczos 迭代法利用最小二乘QR 分解算法進(jìn)行求解[32]，當(dāng)λ為0，阻尼LSQR 算法為常規(guī)標(biāo)準(zhǔn)的LSQR 算法，rIλ(m×n)-A(m×n)Ibλ為殘差向量.

阻尼LSQR-LMBC 混合算法的目標(biāo)函數(shù)為

式中，函數(shù)均在單色(G)輻射強(qiáng)度的前提下進(jìn)行計(jì)算，Iλ(G)為火焰輻射光線(xiàn)的單色輻射強(qiáng)度觀測(cè)值，為火焰輻射光線(xiàn)的單色輻射強(qiáng)度估計(jì)值.該算法的具體步驟如下:

1)根據(jù)目標(biāo)函數(shù)采用的單色輻射強(qiáng)度，選擇性地將像素點(diǎn)的光線(xiàn)輻射強(qiáng)度Iλ(R)或者Iλ(G)輸入.選取單位向量或者隨機(jī)向量作為吸收系數(shù)的初值，這是因?yàn)長(zhǎng)M 算法對(duì)初值沒(méi)有依賴(lài)性.對(duì)中間參數(shù)設(shè)置迭代次數(shù)k0，最大迭代次數(shù)為kmax300，v2，η2，吸收系數(shù)的向量為κκ0.然后對(duì)以下各參數(shù)進(jìn)行計(jì)算，參數(shù)ω，r，g，Stp，μ分別為對(duì)稱(chēng)矩陣、觀測(cè)值和估計(jì)值的差值向量、目標(biāo)函數(shù)的梯度向量，循環(huán)終止指數(shù)和衰減系數(shù).式中參與計(jì)算的參數(shù)J和ωii分別表示函數(shù)關(guān)于變量κ的一階偏微分雅各比矩陣和矩陣ω對(duì)角元素:

2)根據(jù)設(shè)置的單位向量或者隨機(jī)向量的吸收系數(shù)初值κ對(duì)系數(shù)矩陣A進(jìn)行計(jì)算.再通過(guò)阻尼LSQR 算法對(duì)火焰微元體的單色輻射強(qiáng)度Ibλ(G)進(jìn)行計(jì)算:

3)利用已經(jīng)計(jì)算出的火焰溫度T和系數(shù)矩陣A，再通過(guò) (15) 式對(duì)火焰光線(xiàn)的單色輻射強(qiáng)度值進(jìn)行計(jì)算.

4)計(jì)算一階偏微分雅各比矩陣J，根據(jù)步驟1)中(14)式對(duì)ω、r、g進(jìn)行計(jì)算，通過(guò)計(jì)算出的參數(shù)求解得到優(yōu)化步長(zhǎng) Δ (κ):

5)判斷計(jì)算的結(jié)果是否滿(mǎn)足 (17) 式所示的終止條件，如果滿(mǎn)足，則將循環(huán)終止指數(shù)Stp置為1，同時(shí)將運(yùn)算進(jìn)程轉(zhuǎn)到步驟8).反之，將 Δ (κ)疊加于κ得到新的吸收系數(shù)κnew，再將新的吸收系數(shù)投射到吸收系數(shù)的約束空間Q:

其中，l([l1，l2，···，lm]) 和u([u1，u2，···，um])分別為約束區(qū)間Q的下邊界和上邊界.

6)對(duì)迭代指數(shù)ρ進(jìn)行計(jì)算，并判斷ρ是否大于零，當(dāng)大于零時(shí)，轉(zhuǎn)入(20)式和步驟8)，反之，繼續(xù)往下執(zhí)行:

7)根據(jù) (20) 式對(duì)原吸收系數(shù)進(jìn)行更新，再重復(fù)對(duì)ω，r，g進(jìn)行計(jì)算，對(duì)ν 取值2，繼續(xù)對(duì)Stp和μ進(jìn)行計(jì)算:

8)判斷終止指數(shù)是否滿(mǎn)足Stp0，同時(shí)判斷迭代次數(shù)是否滿(mǎn)足k＜kmax，如果滿(mǎn)足這些判斷條件，運(yùn)算進(jìn)程轉(zhuǎn)回步驟2)，反之，結(jié)束循環(huán)，進(jìn)一步采用阻尼LSQR 算法對(duì)火焰的單色輻射強(qiáng)度Ibλ(R)進(jìn)行計(jì)算，最后計(jì)算出火焰的溫度.

3 火焰溫度場(chǎng)重建實(shí)驗(yàn)

3.1 數(shù)值仿真

根據(jù)火焰輻射光場(chǎng)成像模型，利用光線(xiàn)追跡對(duì)火焰輻射圖像進(jìn)行成像模擬.本文采用 (21)和(22) 式火焰分布模型，其中，T為火焰溫度大小，單位為K，z為火焰的軸向坐標(biāo)，r為火焰的徑向坐標(biāo).該模擬火焰符合火焰溫度分布特點(diǎn)，火焰底面半徑R和高度Z均為40 mm，火焰的溫度分布滿(mǎn)足旋轉(zhuǎn)對(duì)稱(chēng)分布.得到的火焰溫度分布如圖2 所示，溫度最高能達(dá)到2100 K.光場(chǎng)相機(jī)的參數(shù)為:主透鏡焦距為50 mm，微透鏡焦距為0.6 mm，微透鏡陣列大小為60*60，每個(gè)微透鏡直徑上覆蓋6 個(gè)像素，像素尺寸為8 μm*8 μm，采用光場(chǎng)相機(jī)1.0模型，微透鏡放置在主透鏡一倍焦距處，圖像探測(cè)器放置在微透鏡一倍焦距處.火焰中心與主透鏡中心的距離為505 mm.微透鏡的數(shù)量決定了記錄的空間物點(diǎn)數(shù)，即空間分辨率，每個(gè)微透鏡覆蓋的像素?cái)?shù)決定了記錄的光線(xiàn)方向數(shù)，即角度分辨率.圖像探測(cè)器與微透鏡的距離決定光場(chǎng)相機(jī)獲得的有效采樣輻射信息，一倍焦距處最低，這個(gè)參數(shù)主要和本文使用的實(shí)驗(yàn)設(shè)備的相機(jī)模型一致.

圖2 模擬火焰中心切片溫度分布Fig.2.Temperature distribution of simulated flame center slice.

以相機(jī)主透鏡中心O作為原點(diǎn)，Y軸垂直于圖像探測(cè)器平面，X和Z軸平行于圖像探測(cè)器平面.沿Y方向?qū)δM火焰進(jìn)行切層，每層間隔4 mm.根據(jù)光場(chǎng)成像模型，遍歷圖像探測(cè)器上像素的坐標(biāo)，逆向定位火焰切片上的微元體坐標(biāo)，將微元體坐標(biāo)代入 (21) 式得到火焰不同切片的溫度值.根據(jù)(15) 式計(jì)算得到火焰微元體的單色輻射Ibλ，由模擬火焰的光線(xiàn)穿過(guò)微元體的數(shù)量和設(shè)定的吸收系數(shù)計(jì)算得到 (8) 式中的矩陣A，進(jìn)而得到圖像傳感器上的火焰單色輻射強(qiáng)度圖，疊加RGB 三個(gè)通道的火焰圖像即可得到模擬光場(chǎng)火焰圖，如圖3所示.

圖3 模擬火焰光場(chǎng)輻射圖像Fig.3.Simulated flame light field radiation image.

分別利用LSQR 算法和阻尼LSQR 算法對(duì)火焰溫度場(chǎng)進(jìn)行重建.重建結(jié)果如圖4 所示，重建時(shí)使用原始溫度疊加20%的白噪聲作為溫度初值，一至三行分別為原始火焰分層、使用LSQR 和阻尼LSQR 對(duì)火焰光場(chǎng)輻射圖像進(jìn)行溫度重建得到的火焰溫度分布圖.從圖4 可以看出，這兩種算法均可以完成火焰的三維溫度重建.

圖4 QR 和阻尼LSQR 算法的火焰溫度場(chǎng)重建結(jié)果Fig.4.Flame temperature field reconstruction results of LSQR and damped LSQR algorithms.

LSQR 算法是一個(gè)對(duì)初值有一定依賴(lài)的算法[33]，初值的誤差對(duì)大型稀疏矩陣的求解的影響較為顯著，使計(jì)算效率變低，甚至無(wú)法收斂.為了定量分析LSQR 與阻尼LSQR 的重建精度，在原始溫度值上疊加了1%，5%，10%，15%和20%噪聲作為初值，從重建結(jié)果中提取了120 個(gè)采樣溫度點(diǎn)，將其與相對(duì)應(yīng)的原始溫度進(jìn)行對(duì)比分析，如圖5 所示.

從圖5 可以看出，隨著初值疊加噪聲的增大，求解值與原始數(shù)據(jù)之間的差值也逐漸變大.在1%的噪聲水平下，LSQR 與阻尼LSQR 的重建誤差差距較小，而隨著噪聲的增大，二者的重建誤差差距開(kāi)始變大，阻尼LSQR 的優(yōu)點(diǎn)逐漸凸顯.當(dāng)在20%的噪聲水平下，LSQR 算法的重建誤差在130 K以?xún)?nèi)，而阻尼LSQR 算法的重建誤差在90 K 以?xún)?nèi).

圖5 不同噪聲水平下，LSQR 與阻尼LSQR 算法重建誤差對(duì)比 (a)初值疊加1%噪聲;(b) 初值疊加5%噪聲;(c) 初值疊加10%噪聲;(d) 初值疊加15%噪聲;(e) 初值疊加20%噪聲Fig.5.Comparison of reconstruction errors between LSQR and Damped LSQR algorithm under different noise levels:(a) Initial value superimposed 1% noise;(b) initial value superimposed 5% noise;(c) initial value superimposed 10% noise;(d) initial value superimposed 15% noise;(e) initial value superimposed 20% noise.

為進(jìn)一步從整體上進(jìn)行算法優(yōu)劣的對(duì)比，將得到的所有重建結(jié)果通過(guò) (23) 式進(jìn)行了相對(duì)誤差計(jì)算:

式中，T為使用重建算法重建出的微元體的溫度，Tset為使用 (21) 式設(shè)置的火焰的微元體的溫度，Δ T為相對(duì)誤差.得到的誤差對(duì)比結(jié)果如表1 所示.

表1 LSQR 與阻尼LSQR 溫度重建結(jié)果相對(duì)誤差對(duì)比表Table 1.Comparison of relative errors of temperature reconstruction results of LSQR and Damped LSQR.

通過(guò)重建結(jié)果可以看出，當(dāng)初值疊加的噪聲較小時(shí)，阻尼LSQR 的略?xún)?yōu)于LSQR.隨著疊加的噪聲逐漸增大，阻尼LSQR 比LSQR 算法更加穩(wěn)定，而LSQR 受初值的影響更大，因此阻尼LSQR 更適合大型稀疏矩陣的求解.孫俊[29]采取的NNLS算法和孫俊陽(yáng)[30]采取的Tikhonov 算法均與LSQR算法進(jìn)行了對(duì)比，雖然與本文的光場(chǎng)相機(jī)模型不同，其相對(duì)的性能提升仍具有參考價(jià)值.對(duì)于NNLS算法，隨著噪聲水平的提高，LSQR 求解結(jié)果中的負(fù)值越來(lái)越多，但NNLS 未出現(xiàn)負(fù)值，兩種算法重建結(jié)果的相對(duì)誤差十分接近，信噪比為20 dB 時(shí)，相對(duì)誤差在5%左右;信噪比為10 dB 時(shí)，相對(duì)誤差在20%左右，信噪比為5 dB 時(shí)，相對(duì)誤差在40%左右.NNLS 相對(duì)LSQR 提升在5%以?xún)?nèi)，但NNLS的計(jì)算時(shí)長(zhǎng)比LSQR 高了兩個(gè)數(shù)量級(jí).對(duì)于Tikhonov算法，噪聲較小時(shí)，Tikhonov 正則化算法與LSQR算法重建結(jié)果的最大相對(duì)誤差基本一致.當(dāng)信噪比達(dá)到10 dB 時(shí)，LSQR 算法重建結(jié)果的相對(duì)誤差在1.8%以?xún)?nèi)，Tikhonov 正則化算法重建結(jié)果的相對(duì)誤差在1.2%以?xún)?nèi)，相對(duì)提升了約50%.對(duì)于本文的阻尼LSQR 算法，當(dāng)信噪比為13.86 dB 時(shí)，LSQR重建結(jié)果的相對(duì)誤差為14.58%，阻尼LSQR重建結(jié)果的相對(duì)誤差為10.99%，相對(duì)提升了32.7%.從計(jì)算時(shí)間上看，本文方法和LSQR 處于同一量級(jí)，計(jì)算的實(shí)時(shí)性對(duì)火焰來(lái)說(shuō)也是一個(gè)比較重要的因素.

對(duì)于實(shí)際的火焰溫度測(cè)量，火焰的吸收系數(shù)等物性參數(shù)是未知的，因此利用阻尼LSQR-LMBC混合算法對(duì)模擬火焰輻射光場(chǎng)圖像的吸收系數(shù)和溫度同時(shí)重建.火焰以4 mm 的間隔對(duì)火焰進(jìn)行分層，重建溫度分布如圖6 所示.圖中的火焰溫度分布圖可以正確的反映火焰的溫度趨勢(shì)特征，但是重建溫度的精度有待提高.

圖6 阻尼LSQR-LMBC 算法下重建的火焰溫度分布圖Fig.6.Flame temperature distribution reconstructed by Damped LSQR-LMBC algorithm.

為了進(jìn)一步對(duì)重建結(jié)果進(jìn)行分析，本文對(duì)重建后的前五層溫度進(jìn)行誤差分析，將計(jì)算出的每一個(gè)微元體的溫度與相對(duì)應(yīng)的原始溫度進(jìn)行相對(duì)誤差的計(jì)算，其結(jié)果如圖7 所示.

從圖7 可以看出，該算法的總體重建誤差在1.92%—17.2%之間.所有微元體的平均誤差為6.63%，因此該重建算法可以滿(mǎn)足高達(dá)2100 K 的火焰溫度重建.

圖7 阻尼LSQR-LMBC 算法重建的火焰溫度的相對(duì)誤差 (a) —16 mm 位置;(b) 火焰—12 mm 位置;(c) —8 mm 位置;(d) —4 mm位置;(e) 0 mm 位置Fig.7.Relative errors of flame temperature reconstructed by Damped LSQR－LMBC algorithm:(a) —16 mm position;(b) —12 mm position;(c) —8 mm position;(d) —4 mm position;(e) 0 mm position.

3.2 輻射強(qiáng)度標(biāo)定

光場(chǎng)相機(jī)圖像傳感器在對(duì)火焰進(jìn)行輻射信息采集時(shí)，輻射光線(xiàn)會(huì)先經(jīng)過(guò)光電轉(zhuǎn)換、模數(shù)轉(zhuǎn)換等過(guò)程再以灰度值的形式展示.因此，需要對(duì)光場(chǎng)相機(jī)的圖像傳感器進(jìn)行輻射強(qiáng)度標(biāo)定，建立灰度值與光線(xiàn)輻射強(qiáng)度之間的關(guān)系，進(jìn)而保證重建溫度的準(zhǔn)確度.實(shí)驗(yàn)裝置如圖8 所示.使用的黑體爐是由英國(guó)愛(ài)松特公司生產(chǎn)的，型號(hào)為ISOTECH Cyclops Model 878，溫度最高可以達(dá)到1573.15 K，精度為0.1 K，輻射率為0.999.相機(jī)為一代傳統(tǒng)光場(chǎng)相機(jī)，型號(hào)為L(zhǎng)ytro-Illum[34，35]，主透鏡焦距為40.11 mm，微透鏡焦距為48.38 μm，角度分辨率為15×15，空間分辨率為434×625.通過(guò)Lytro-Illum 光場(chǎng)相機(jī)拍攝得到的光場(chǎng)圖像為L(zhǎng)FR 格式，其中包含了原始光場(chǎng)圖像和拍攝時(shí)的參數(shù)信息.解碼時(shí)首先根據(jù)相機(jī)自帶的白圖像標(biāo)定微透鏡陣列的位置坐標(biāo)，微透鏡下覆蓋的像素記錄了來(lái)自同一物點(diǎn)不同方向的光線(xiàn)，這些像素集合又被稱(chēng)為宏像素.每個(gè)宏像素中相同坐標(biāo)的像素點(diǎn)記錄相同方向的光線(xiàn)，提取記錄相同方向光線(xiàn)的像素，按照宏像素位置排列即可得到2D 單視角圖像.依次提取不同方向的單視角圖像，按照記錄的方向排列即可得到多視角圖像.

圖8 光場(chǎng)相機(jī)輻射強(qiáng)度標(biāo)定實(shí)驗(yàn)裝置圖Fig.8.Experimental deviceofradiation intensity calibration of light field camera.

為避免其他雜光對(duì)標(biāo)定實(shí)驗(yàn)產(chǎn)生影響，在黑暗條件下，通過(guò)黑體爐溫度從1098.15 K 開(kāi)始，以25 K 的溫度間隔采集了8 組黑體輻射光場(chǎng)圖像，如圖9 所示.

普朗克定律指出，在黑體的條件下，在指定溫度下的輻射強(qiáng)度值可以通過(guò) (24) 式計(jì)算:

式中，E(λ，T) 為波長(zhǎng)λ處的光譜輻射強(qiáng)度值，ε(λ)為波長(zhǎng)λ處的黑體發(fā)射率，取值為1，T為黑體的溫度值，c1，c2分別為第一輻射常數(shù)和第二輻射常數(shù)，其值大小分別為 3.74177×10-16W·m2和1.4387752×10-2m·K.

對(duì)光場(chǎng)圖像標(biāo)定區(qū)域的灰度值進(jìn)行提取，再分別將R，G，B通道的灰度值與輻射強(qiáng)度值進(jìn)行曲線(xiàn)擬合，得到的擬合結(jié)果如 (25) 式—(27)式所示:

式中，IR，IG，IB分別為光場(chǎng)圖像的R，G，B通道下像素對(duì)應(yīng)的光譜輻射強(qiáng)度值，R，G，B分別為光場(chǎng)圖像的R，G，B通道下像素對(duì)應(yīng)的灰度值.光場(chǎng)輻射強(qiáng)度標(biāo)定結(jié)果對(duì)應(yīng)的擬合曲線(xiàn)圖如圖10 所示.

圖10 光場(chǎng)圖像灰度值與輻射強(qiáng)度的標(biāo)定擬合曲線(xiàn)圖Fig.10.Calibration fitting curve of gray value and radiation intensity of light field image.

3.3 實(shí)際火焰溫度重建

實(shí)際火焰的三維測(cè)溫裝置主要有光學(xué)平臺(tái)、光場(chǎng)相機(jī)以及火焰發(fā)生裝置等.將實(shí)驗(yàn)器材固定在光學(xué)平臺(tái)上，通過(guò)光場(chǎng)相機(jī)對(duì)火焰進(jìn)行圖像的采集.實(shí)驗(yàn)裝置圖如圖11(a)所示，圖中的丁烷氣罐的燃燒火焰，火焰發(fā)生器的噴嘴直徑為20 mm.

在對(duì)火焰進(jìn)行圖像采集時(shí)，要記錄火焰與相機(jī)之間的距離，并且保證火焰處于空氣流動(dòng)很微弱的暗室實(shí)驗(yàn)環(huán)境中.同時(shí)，需要合理控制相機(jī)的曝光時(shí)間，確保采集到的火焰圖像像素灰度值不會(huì)出現(xiàn)曝光過(guò)度，降低曝光不足帶來(lái)的影響.通過(guò)光場(chǎng)相機(jī)采集的丁烷火焰輻射光場(chǎng)圖像如圖11(b)所示，火焰與光場(chǎng)相機(jī)之間的距離為15 cm.

圖11 燃燒火焰實(shí)驗(yàn)圖 (a) 實(shí)驗(yàn)裝置;(b) 丁烷火焰輻射光場(chǎng)Fig.11.Experimental diagram of combustion flame:(a) Experimental device;(b) light field of butane flame radiation.

根據(jù)相機(jī)與火焰之間的距離以及相機(jī)內(nèi)部的參數(shù)等已知數(shù)據(jù)，利用 (1) 式— (5) 式所示的雙平面參數(shù)法將光場(chǎng)輻射圖像的方向信息轉(zhuǎn)換為三維空間的火焰位置信息.得到火焰的三維位置信息后使用微元體分割方法對(duì)火焰進(jìn)行網(wǎng)格劃分.通過(guò)輻射標(biāo)定結(jié)果將圖像灰度值幻化為輻射強(qiáng)度值，再利用阻尼LSQR-LMBC 混合算法對(duì)火焰的吸收系數(shù)和溫度同時(shí)重建，結(jié)果如圖12 所示.

圖12 丁烷火焰三維溫度重建后的不同層溫度分布圖Fig.12.Temperature distributions of different layers after three-dimensional temperature reconstruction of butane flame.

由丁烷發(fā)生器的噴嘴大小可知火焰的最大直徑為20 mm，圖中將火焰在2 mm，4 mm，6 mm，8 mm，10 mm，12 mm，14 mm，16 mm，18 mm處進(jìn)行豎直切片得到的各層溫度，其中10 mm 處的溫度分布圖為火焰的中心層.由于該火焰近乎旋轉(zhuǎn)對(duì)稱(chēng)，因此火焰在12 mm，14 mm，16 mm，18 mm處的溫度分布圖分別與8 mm，6 mm，4 mm，2 mm處的溫度分布圖近乎一致.圖中重建火焰的最高溫度在1400 K 左右，不超過(guò)1500 K.從10 mm 層的溫度分布圖中可以看出，火焰滿(mǎn)足焰心、內(nèi)焰和外焰的溫度分布趨勢(shì).外圍的溫度會(huì)出現(xiàn)低溫度區(qū)域的原因是由于火焰的溫度出現(xiàn)流失，擴(kuò)散到了低溫的空氣中，這也是2 mm 和18 mm 處的溫度較低的原因.而8 mm 與12 mm 處的火焰切片的較大區(qū)域?yàn)閮?nèi)焰和外焰，火焰較為明亮，溫度也較高.6 mm 與14 mm 處的火焰溫度切片的面積已經(jīng)明顯變小，外焰和內(nèi)焰的混合使得火焰的溫度較高.4 mm 和16 mm 處的火焰溫度已經(jīng)逐步擴(kuò)散到空氣中，溫度也有所降低.綜上，重建的丁烷火焰的三維溫度場(chǎng)分布符合基本輻射火焰燃燒的特征.

為進(jìn)一步對(duì)重建結(jié)果進(jìn)行準(zhǔn)確度驗(yàn)證，使用熱電偶對(duì)火焰進(jìn)行定點(diǎn)測(cè)溫，裝置如圖13(a)所示.選取不同位置的特征點(diǎn)進(jìn)行測(cè)溫，測(cè)溫點(diǎn)的選取如圖13(b)所示.熱電偶在測(cè)量時(shí)，其測(cè)量點(diǎn)和周?chē)h(huán)境之間存在輻射換熱和對(duì)流換熱，從而導(dǎo)致產(chǎn)生熱損失，為了減小引起的測(cè)量誤差，需要使用 (28)式對(duì)測(cè)量結(jié)果進(jìn)行修正[36]:

圖13 使用熱電偶測(cè)量丁烷火焰溫度 (a) 熱電偶測(cè)溫實(shí)驗(yàn)圖;(b) 熱電偶定點(diǎn)測(cè)溫位置圖Fig.13.Measurement of butane flame temperature using thermocouples:(a) Thermocouple temperature measurement experiment diagram;(b) location diagram of thermocouple fixed-point temperature measurement.

式中，T表示修正后的溫度值;Tc和Ts分別為測(cè)量的原始溫度和實(shí)驗(yàn)環(huán)境溫度，實(shí)驗(yàn)時(shí)室溫控制在293.15 K 左右;ε表示熱電偶探頭發(fā)射率，實(shí)驗(yàn)中使用的是K 型熱電偶，數(shù)值為0.89;σ表示斯蒂芬-玻爾茲曼常數(shù)，數(shù)值為5.67×10—8W/(m2·K4);hc表示對(duì)流換熱系數(shù)[29].

將重建溫度和熱電偶修正溫度進(jìn)行比對(duì)，如表2 所示.

從表2 可以看出，重建溫度與熱電偶的測(cè)溫?cái)?shù)據(jù)基本一致.其中測(cè)量偏差的平均值為76 K 左右，最大偏差約為89 K.測(cè)溫誤差控制在6.8%左右，因此本研究的重建方案的測(cè)溫精度滿(mǎn)足火焰的測(cè)溫標(biāo)準(zhǔn).

表2 阻尼LSQR-LMBC 重建溫度與熱電偶修正溫度對(duì)比表Table 2.Comparison of DampedLSQR-LMBC reconstruction temperature and thermocouple correction temperature.

4 結(jié)論

1) 本文針對(duì)經(jīng)典算法LSQR 在求解火焰三維溫度場(chǎng)時(shí)對(duì)初值依賴(lài)大的問(wèn)題，通過(guò)添加阻尼正則化項(xiàng)，提高了火焰三維溫度場(chǎng)的抗噪性能，仿真結(jié)果表明:在較低的噪聲水平下，這兩種算法的重建誤差基本一致，當(dāng)噪聲較大時(shí)，阻尼LSQR 算法的重建精度比LSQR 算法高30%左右.鑒于實(shí)際火焰的吸收系數(shù)未知，使用阻尼LSQR-LMBC 算法對(duì)吸收系數(shù)該火焰進(jìn)行數(shù)值計(jì)算，實(shí)驗(yàn)結(jié)果表明:火焰重建的相對(duì)誤差在1.92%—17.2%之間，平均誤差低于7%.證明了阻尼LSQR-LMBC 算法可以同時(shí)重建火焰的吸收系數(shù)和溫度.

2) 開(kāi)展了實(shí)際火焰溫度重建實(shí)驗(yàn)研究.搭建了光場(chǎng)成像的火焰三維溫度場(chǎng)重建的實(shí)驗(yàn)平臺(tái)，標(biāo)定了光場(chǎng)相機(jī)圖像傳感器的輻射強(qiáng)度，對(duì)丁烷火焰的三維溫度進(jìn)行重建，并將熱電偶的測(cè)量結(jié)果作為參考火焰溫度，將其與重建后的結(jié)果進(jìn)行對(duì)比分析和精度評(píng)價(jià).結(jié)果表明:火焰重建的溫度與熱電偶測(cè)量的溫度基本一致，平均偏差均低于100 K，重建誤差在6.8%左右，