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

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

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

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

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

1 引 言

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

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

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

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

2.1 測量原理

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

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

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

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

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

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

3 重建算法

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

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

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

3.1 傳統(tǒng)LTD算法

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

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

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

3.2 LTD-TV算法

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

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

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

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

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

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

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

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

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

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

4 數(shù)值模擬

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

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

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

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

圖2 傳統(tǒng)LTD法與LTD-TV法比較(a),(b),(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法有很大的改善,大大減小了重建的偽影.該算法不但能重建出煙羽擴散的大致情況,還能重建出豐富的細節(jié).然而,LTD法假設(shè)氣體濃度是位置的二階多項式,重建過程其實是用二階多項式去擬合高斯函數(shù),導(dǎo)致其峰值濃度比真實濃度低10%—15%,擴散范圍比測試圖形稍大,總體上有變圓的趨勢.

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

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

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

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

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

4.3 大氣流動對重建結(jié)果的影響

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

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

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

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

5 外場實驗

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

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

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

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

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

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

圖7 測量路徑積分濃度與重建路徑積分濃度對比(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ā)生位移時的校正因子,

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

6 結(jié)果與討論

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