基于稀疏優(yōu)化的煙羽斷層重建方法*

鐘鳴宇 奚亮 司福祺? 周海金 王煜

1)(中國(guó)科學(xué)院安徽光學(xué)精密機(jī)械研究所,環(huán)境光學(xué)與技術(shù)重點(diǎn)實(shí)驗(yàn)室,合肥 230031)

2)(中國(guó)科學(xué)技術(shù)大學(xué)研究生院科學(xué)島分院,合肥 230031)

3)(安徽理工大學(xué)電氣與信息工程學(xué)院,淮南 232001)

1 引 言

燃煤發(fā)電在我國(guó)能源結(jié)構(gòu)中占比高達(dá)65%,為經(jīng)濟(jì)社會(huì)發(fā)展和人民生活提供堅(jiān)實(shí)的能源保障的同時(shí),排放SO2等也給環(huán)境帶來(lái)了巨大的壓力.最新的研究表明,大氣中的有機(jī)物與硫酸在大氣新粒子的形成和增長(zhǎng)中具有重要作用[1].而燃煤電廠排放的硫超過(guò)總量的40%,因此研究電廠SO2的排放情況及空間分布具有重要意義.

差分吸收光譜儀(differential optical absorption spectroscopy,DOAS)僅能測(cè)得氣體濃度沿光傳輸路徑的積分,而斷層圖像重建使這種大氣參數(shù)的間接測(cè)量成為實(shí)用技術(shù)[2],在大氣遙感領(lǐng)域得以廣泛應(yīng)用.用于斷層重建的光學(xué)遙感設(shè)備通常僅有2臺(tái)且固定在地面,氣體斷層重建是典型的不完全角度重建.早期的研究通過(guò)直接解方程來(lái)重建氣體分布,重建圖像質(zhì)量不高[3-5].一種改進(jìn)方法是對(duì)氣體建模作為先驗(yàn)信息,為方程組增加約束.常用的模型有Drescher 等[6]提出的高斯模型、Price等[7]提出的低三階導(dǎo)數(shù)(low third derivative method,LTD)模型和Olaguer等[8-10]提出的基于流體力學(xué)的歐拉方程模型.LTD模型假設(shè)氣體濃度對(duì)位置的三階導(dǎo)數(shù)值全部為零,相當(dāng)于隱含地假設(shè)了氣體濃度嚴(yán)格滿足空間位置的二階多項(xiàng)式.Johansson等[11]分別將該方法用于重建火山和煙囪煙羽分布.Kazahaya等[12]將LTD項(xiàng)乘以權(quán)重系數(shù)后和投影方程相加,并使用最小二乘法求解,但重建的結(jié)果仍然不夠理想.Casaballe等[13]利用LTD模型,對(duì)投影方程組進(jìn)行Tikhonov正則化,并使用cvx優(yōu)化工具箱進(jìn)行求解,取得了很好的數(shù)值模擬效果,但該方法實(shí)際抗誤差能力弱,無(wú)法用于外場(chǎng)實(shí)驗(yàn)數(shù)據(jù).總的來(lái)說(shuō),重建算法還很不成熟,重建圖像存在大量偽影.壓縮感知(compressed sensing,CS)理論為氣體分布重建提供了新的思路: 如果能找到氣體分布在某種變換下具有稀疏性,就可以在采樣數(shù)很少的情況下精確重建氣體分布[14-17].

本文使用成像差分吸收光譜儀(imaging differential optical absorption spectroscopy,IDOAS)采集數(shù)據(jù),提出了一種基于低三階模型的全變分(low third derivative total variation,LTD-TV)法重建煙羽分布.該方法基于LTD模型和壓縮感知理論對(duì)煙羽分布進(jìn)行重建,首先使用代數(shù)重建算法(algebraic reconstruction technique,ART)對(duì)重建數(shù)據(jù)進(jìn)行初始化,用對(duì)數(shù)障礙函數(shù)法[17]確定目標(biāo)函數(shù),全變分法確定下降梯度,最后使用Barzilai-Borwein(BB)算法[18]確定步長(zhǎng)后對(duì)重建結(jié)果進(jìn)行優(yōu)化.本文對(duì)該方法進(jìn)行了數(shù)值模擬,并進(jìn)行了外場(chǎng)實(shí)驗(yàn),重建了煙囪煙羽的斷層2維分布,并對(duì)重建結(jié)果進(jìn)行了分析.

2 測(cè)量原理與IDOAS斷層掃描系統(tǒng)

2.1 測(cè)量原理

DOAS技術(shù)基于Lambert-Beer定律

其中I0(λ)為大氣層外的太陽(yáng)光強(qiáng),I(λ)為經(jīng)過(guò)大氣層后,DOAS接收到的太陽(yáng)光強(qiáng),σm(λ)為第m種氣體的吸收截面,

其中Sm表示該種氣體的斜柱濃度,等于氣體的濃度cm對(duì)光程r積分,積分距離為L(zhǎng).通過(guò)最小二乘法求解(1)式,可以得到污染氣體的斜柱濃度[19,20].

2.2 斷層掃描系統(tǒng)

如圖1所示,掃描煙羽的過(guò)程中,假設(shè)風(fēng)沿Z軸方向,虛擬的掃描平面垂直于大地,煙羽被該虛擬平面截取了一個(gè)平面,兩臺(tái)IDOAS放置在掃描平面與大地的交線上.將IDOAS#1指向IDOAS#2的方向設(shè)為X軸的正方向,豎直向上設(shè)為Y軸正方向.掃描區(qū)域離散化如圖1中的虛線所示.兩臺(tái)IDOAS的視場(chǎng)角為30°,相鄰掃描線間隔0.625°,在2 s內(nèi)采樣并存儲(chǔ)48個(gè)點(diǎn)的柱濃度數(shù)據(jù).而采用以往常用的多軸差分吸收光譜儀(multi-axis differential optical absorption spectroscopy,MAX-DOAS),采集相同數(shù)量的數(shù)據(jù)需要5 min以上.IDOAS采集數(shù)據(jù)的時(shí)間分辨率比MAX-DOAS提高了160多倍.在正式掃描之前先進(jìn)行預(yù)掃描,調(diào)整IDOAS的仰角,使煙羽位于IDOAS的視場(chǎng)角內(nèi).將掃描區(qū)域劃分為多個(gè)網(wǎng)格,在這些網(wǎng)格上使用重建算法重建煙羽的濃度分布.

圖1 掃描區(qū)域離散化Fig.1.Scanning region discretization.

3 重建算法

在離散情況下,(2)式可表示為

S(i)為 第i條射線的路徑積分濃度,C(j)表示掃描區(qū)域第j個(gè)像素中的氣體平均濃度,H(i,j)表示第i條射線穿越第j個(gè)像素的長(zhǎng)度.把整個(gè)系統(tǒng)的投影系數(shù)和濃度寫(xiě)成矩陣與向量相乘的形式

其中H表示投影矩陣,C表示圖1中重建區(qū)域中像素中排列成的列向量,通常C中像素的數(shù)目遠(yuǎn)大于S中射線的數(shù)目,所以(4)式是一個(gè)欠定方程組.

3.1 傳統(tǒng)LTD算法

為了解方程組(4),LTD模型假設(shè)氣體濃度相對(duì)位置的三階導(dǎo)數(shù)等于0[7],于是有

式中c(k,l)是圖1所示重建區(qū)域中像素矩陣c的第k行l(wèi)列像素.把(5)式和(6)式改寫(xiě)為矩陣的形式

將(4)式和(7)式聯(lián)立得到一個(gè)過(guò)定方程組,求解就得到了濃度C的近似值.

3.2 LTD-TV算法

在傳統(tǒng)的LTD方法中,隱含地假設(shè)氣體的濃度是位置的二階多項(xiàng)式.然而,氣體的濃度不可能?chē)?yán)格地按二階多項(xiàng)式分布,(7)式左端直接設(shè)為零顯然會(huì)導(dǎo)致圖像像素間約束過(guò)強(qiáng),圖像邊緣出現(xiàn)大量偽峰,嚴(yán)重影響重建氣體的實(shí)際分布.而本文提出的LTD-TV算法,僅要求氣體的濃度大體按照二階多項(xiàng)式分布,也就是假設(shè)氣體濃度值在三階導(dǎo)數(shù)下是稀疏的,該假設(shè)顯然比(7)式合理得多.

LTD-TV法首先使用ART算法對(duì)重建圖像進(jìn)行初始化,然后使用基于全變分的優(yōu)化算法優(yōu)化目標(biāo)函數(shù).目標(biāo)函數(shù)為

式中C?即為所求的氣體濃度分布,表示氣體濃度分布的全變分的模,ε和鏡頭接收的光子數(shù)目的泊松分布有關(guān),這里簡(jiǎn)單設(shè)置為10—12,不作進(jìn)一步的深入討論.使用對(duì)數(shù)障礙函數(shù)法[17],方程(8)中的約束條件可以寫(xiě)入優(yōu)化問(wèn)題中

如果使用低三階導(dǎo)數(shù)模型,根據(jù)(5)式和(6)式,(9)式中氣體分布的全變分表達(dá)式為

對(duì)(9)式右邊括號(hào)中的部分求導(dǎo),得到(9)式的梯度為

(11)式中,全變分的梯度的計(jì)算公式為

使用優(yōu)化算法求解(13)式,αn是優(yōu)化算法的步長(zhǎng),pn是優(yōu)化算法的下降方向.

式中的gn(j)由(11)式確定,j和(3)式中j的含義相同,表示第j個(gè)像素,(14)式隱含地給氣體濃度添加了非負(fù)約束.(13)式中αn用BB算法得到

式中sn-1=Cn-Cn-1,yn-1=pn-pn-1,分別表示前后兩次迭代的氣體濃度分布差,和前后兩次迭代的梯度差.當(dāng)(13)式中兩次迭代間的圖像差別很小時(shí),說(shuō)明算法收斂.

式中N表示圖像中像素的總數(shù).當(dāng)σresidual小于某個(gè)閾值,或者迭代次數(shù)大于某個(gè)數(shù)時(shí)迭代停止.

4 數(shù)值模擬

假設(shè)真實(shí)的氣體濃度分布服從高斯模型,掃描系統(tǒng)如圖1所示,兩臺(tái)IDOAS分別位于X軸的左右兩端,Y軸代表氣體所在的高度,掃描線間隔為0.625°,重建圖像解析度為 20×20.為了比較方便,將高斯函數(shù)的濃度歸一化為1.在知道氣體真實(shí)分布的情況下,使用接近度σnearness作為重建效果的指標(biāo)

其中C?(j)是測(cè)試圖像第j個(gè)像素中的氣體濃度,C(j)是 重建圖像第j個(gè)像素的氣體濃度,是所有像素氣體濃度的平均值.接近度越小,重建效果越好.

4.1 LTD-TV算法與傳統(tǒng)LTD算法重建結(jié)果對(duì)比

圖2給出幾幅用傳統(tǒng)的LTD算法和LTD-TV算法重建氣體分布的數(shù)值模擬等高線圖.圖2(a)—圖2(c)分別是單高斯、雙高斯、三高斯氣體擴(kuò)散模型.

圖2 傳統(tǒng)LTD法與LTD-TV法比較(a),(b),(c)測(cè)試圖形;(d),(e),(f)傳統(tǒng)LTD法重建圖形;(g),(h),(i)LTD-TV法20000次迭代重建圖形Fig.2.Comparison between traditional LTD algorithm and LTD-TV algorithm:(a),(b),(c)Test distribution;(d),(e),(f)reconstruction of distribution using traditional LTD algorithm;(g),(h),(i)reconstruction of distribution using LTD-TV algorithm with 20000 iterations.

圖2(d)—圖2(f)使用傳統(tǒng)LTD法重建,接近度分別為0.6363,0.5930,0.5778; 圖2(g)—圖2(i)是在循環(huán)20000次的情況下用LTD-TV法重建,接近度分別為0.1127,0.1052,0.1995,分別比傳統(tǒng)的LTD算法接近度減小82.28%,82.25%,65.47%.從圖2中可以看出,LTD-TV法重建圖形相比傳統(tǒng)的LTD法有很大的改善,大大減小了重建的偽影.該算法不但能重建出煙羽擴(kuò)散的大致情況,還能重建出豐富的細(xì)節(jié).然而,LTD法假設(shè)氣體濃度是位置的二階多項(xiàng)式,重建過(guò)程其實(shí)是用二階多項(xiàng)式去擬合高斯函數(shù),導(dǎo)致其峰值濃度比真實(shí)濃度低10%—15%,擴(kuò)散范圍比測(cè)試圖形稍大,總體上有變圓的趨勢(shì).

4.2 測(cè)量誤差對(duì)重建結(jié)果的影響

外場(chǎng)實(shí)驗(yàn)中,溫度、濕度、氣溶膠等都會(huì)影響DOAS的精確性.SO2的測(cè)量值往往具有2%—20%的測(cè)量誤差,不好的重建算法可能會(huì)使誤差在整個(gè)圖形擴(kuò)散,造成重建結(jié)果偏離真實(shí)值,甚至完全淹沒(méi)在重建噪聲中.因此,重建算法的抗誤差能力是考察重建算法的一個(gè)重要方面.由于誤差來(lái)源的復(fù)雜性和不可重復(fù)性[22,23],一般采用疊加隨機(jī)數(shù)的方法模擬誤差對(duì)測(cè)量的影響[13].給圖2中測(cè)試圖形的路徑積分濃度Sk加入加性噪聲ΔSk

式中f表示誤差系數(shù),Rrand是方差為1的隨機(jī)數(shù).圖3給出了f從1%—20%變化的情況下,用LTD-TV重建圖2中的測(cè)試圖像得到的σnearness曲線.

從圖3可以看出,隨著誤差系數(shù)f變大,LTDTV算法的重建接近度也逐漸變大.當(dāng)誤差系數(shù)達(dá)到某個(gè)值時(shí),接近度的上升變慢,說(shuō)明誤差越大,LTD-TV算法的抗誤差效果越明顯.從圖3還可看出高斯數(shù)為2時(shí)重建誤差最小,對(duì)比圖2(a)—圖2(c)可以看出,高斯數(shù)為2的測(cè)試圖形的非零元素最少,圖像本身的稀疏性最強(qiáng),因此LTDTV算法的性能也發(fā)揮得更加充分.而在誤差系數(shù)f變化的情況下,傳統(tǒng)LTD算法重建接近度仍然為0.6363,0.5930,0.5778,比使用LTD-TV法重建的接近度更大.總的來(lái)說(shuō),LTD-TV算法比傳統(tǒng)的LTD算法的重建質(zhì)量好得多.

圖3 接近度隨誤差系數(shù)變化曲線Fig.3.Nearness as the functions of error factors.

4.3 大氣流動(dòng)對(duì)重建結(jié)果的影響

已有的研究表明,當(dāng)煙羽位于重建圖像的中心,且兩臺(tái)DOAS與煙羽中心的連線之間成90°夾角時(shí),重建效果最好[11,13].由于風(fēng)向的不確定性,該條件往往得不到滿足.圖4中給出了圖2(a)—圖2(c)中高斯煙羽模型偏離初始位置時(shí),接近度隨偏離距離變化的曲線,偏移距離小于0時(shí),表示煙羽位置左移,偏移距離大于0時(shí),表示煙羽位置右移.

圖4 接近度偏離距離變化曲線Fig.4.Nearness as the functions of offset distance.

從圖4中可以看出,當(dāng)煙羽偏離圖像中心時(shí),重建接近度和氣體分布的稀疏性仍然相關(guān).當(dāng)單高斯煙羽和雙高斯煙羽偏離初始位置時(shí),重建圖像的接近度發(fā)生振蕩,這種振蕩是重建圖像離散化造成的,與以往的算法相類似[11].除此之外,接近度曲線未表現(xiàn)出明顯的變化規(guī)律,說(shuō)明煙羽稀疏度較高的情況下,LTD-TV算法不受煙羽與IDOAS相對(duì)位置的影響.而當(dāng)三高斯煙羽模型偏離圖像中心時(shí),重建接近度發(fā)生變化.這是因?yàn)闊熡鹬亟ㄊ且环N極端的不完全角度重建,當(dāng)煙羽中氣體濃度分布較為復(fù)雜時(shí),不同的煙羽位置會(huì)對(duì)IDOAS采集數(shù)據(jù)造成較大的影響,從而造成重建效果發(fā)生起伏.此外,外場(chǎng)實(shí)驗(yàn)時(shí),特別是在湍流的影響下,要表示煙羽中氣體濃度的復(fù)雜分布,需要更多的像素?cái)?shù)目.從信息論的角度來(lái)看,在兩臺(tái)IDOAS信息采集能力有限的情況下,低三階導(dǎo)數(shù)模型增加的先驗(yàn)知識(shí)會(huì)將IDOAS采集的數(shù)據(jù)淹沒(méi),導(dǎo)致重建結(jié)果偏離真實(shí)分布.因此,如果需要進(jìn)一步提高LTDTV算法的重建圖像的分辨率,只能進(jìn)一步增加儀器的數(shù)量[5,24].

基于煙羽位置與形狀的復(fù)雜性,使用LTDTV法時(shí)仍然建議將IDOAS關(guān)于煙囪對(duì)稱設(shè)置,且兩臺(tái)DOAS與煙羽中心的連線之間成90°夾角.

5 外場(chǎng)實(shí)驗(yàn)

外場(chǎng)實(shí)驗(yàn)數(shù)據(jù)在淮南某電廠外獲得,電廠煙囪高210 m,在下風(fēng)約180 m處搭建由兩臺(tái)IDOAS組成的斷層掃描系統(tǒng),兩臺(tái)儀器距離約504 m.實(shí)驗(yàn)當(dāng)天天氣晴朗且能見(jiàn)度高,風(fēng)向?yàn)闁|北偏東風(fēng),風(fēng)速為 2 m/s.通過(guò) 307.5—318 nm波段反演SO2路徑積分濃度,參與反演的氣體包括SO2(293 K,vanDaele),NO2(294 K,vanDaele),O3(243 K,Voigt),O3(218 K,Brion)和ring光譜,選取煙囪上風(fēng)向的第一次測(cè)量譜作為參考譜.圖5以其中某一條光譜為反演實(shí)例,展示了SO2光譜數(shù)據(jù)的反演情況.

圖5 SO2柱濃度IDOAS擬合反演實(shí)例(a)SO2柱濃度;(b)擬合殘差Fig.5.An example of SO2 SCD from IDOAS retrieval:(a)SO2 SCD;(b)fitting residual.

圖5(a)中黑線(不規(guī)則曲線)為測(cè)量光譜,紅線(光滑曲線)為擬合光譜,SO2柱濃度為1.11×1017molecules·cm—2.圖5(b)所示的擬合殘差小于0.00322.選取其中一組SO2數(shù)據(jù)作為研究對(duì)象,由于氣體濃度不可能為負(fù),重建圖像時(shí),小于零的路徑積分濃度被設(shè)置為0[5].圖6給出了使用LTDTV法重建的污染氣體斷層圖像.

在圖6中,IDOAS#1位于X軸0點(diǎn),IDOAS#2位于X軸504 m處.煙囪中排出的SO2在空中分成兩股,符合現(xiàn)場(chǎng)目測(cè).能夠非常清楚地看到煙羽從兩個(gè)中心擴(kuò)散的情況,甚至可以觀察到風(fēng)向在略微向IDOAS#2的方向變化,SO2截面向該方向移動(dòng).重建圖中僅有右下角存在少量偽影,其他位置有零星偽影,不影響對(duì)重建圖像的觀察.IDOAS測(cè)得的柱濃度與重建圖像的柱濃度值之間的關(guān)系如圖7所示,直觀上看,柱濃度的測(cè)量值與重建值符合得很好.

當(dāng)氣體分布未知時(shí),一般使用一致性相關(guān)因子衡量重建效果[21].

圖6 煙囪煙羽SO2分布重建圖Fig.6.Reconstruction of SO2 of stack plume.

圖7 測(cè)量路徑積分濃度與重建路徑積分濃度對(duì)比(a)IDOAS #1;(b)IDOAS #2Fig.7.Comparison between measured path integrated concentration and reconstructed path integrated concentration:(a)IDOAS#1;(b)IDOAS #2.

式中ρ表示皮爾遜相關(guān)系數(shù),A是路徑積分濃度曲線發(fā)生位移時(shí)的校正因子,

其中S是IDOAS測(cè)得的路徑積分濃度,σS是測(cè)得的路徑積分濃度的標(biāo)準(zhǔn)差;是重建圖像的路徑積分濃度,是重建圖像路徑積分濃度的標(biāo)準(zhǔn)差,如果CCCF＞0.80 ,即可認(rèn)為重建圖像的路徑積分濃度和實(shí)驗(yàn)測(cè)量值擬合得很好,而圖7所示重建結(jié)果的一致性相關(guān)因子為0.9063.重建得到的路徑積分濃度比測(cè)量值稍低,符合數(shù)值模擬中出現(xiàn)的情況.

6 結(jié)果與討論

本文首次將CS理論引入氣體斷層重建領(lǐng)域,提出了一種新的氣體分布重建方法—LTD-TV算法.數(shù)值模擬表明,該方法將接近度由傳統(tǒng)LTD法的0.5—0.6降低到0.1—0.2,最大降低幅度達(dá)82.28%.針對(duì)加性噪聲,LTD-TV法表現(xiàn)出很強(qiáng)的抗噪聲能力,在路徑積分濃度中引入20%的加性噪聲的情況下,LTD-TV法仍然能重建出氣體擴(kuò)散的主要特征.在外場(chǎng)實(shí)驗(yàn)中使用IDOAS采集數(shù)據(jù),并重建了的煙羽分布圖形.該方法不但假設(shè)氣體的三階導(dǎo)數(shù)是稀疏的,而且隱含的假設(shè)氣體濃度值非負(fù),相當(dāng)于給投影方程組添加了更多的約束.因此與以往的煙羽重建方法相比,提高了氣體重建的時(shí)間分辨率和可信度,極大減少了偽影.該方法不但能用于豎直方向的煙羽分布重建,還能用于水平方向的區(qū)域氣體分布重建.不足之處在于: 使用二階多項(xiàng)式作為氣體擴(kuò)散模型,當(dāng)氣體近似按高斯函數(shù)分布時(shí),重建結(jié)果的峰值偏低,擴(kuò)散范圍偏大.