基于成像差分吸收光譜技術(shù)和壓縮感知理論的煙囪煙羽斷層重建研究

鐘鳴宇， 周海金， 司福祺*， 王 煜， 竇 科， 蘇靜明

1. 安徽理工大學(xué)電氣與信息工程學(xué)院， 安徽 淮南 232001 2. 中國(guó)科學(xué)院安徽光學(xué)精密機(jī)械研究所環(huán)境光學(xué)與技術(shù)重點(diǎn)實(shí)驗(yàn)室， 安徽 合肥 230031 3. 中國(guó)科學(xué)技術(shù)大學(xué)研究生院科學(xué)島分院， 安徽 合肥 230031

引 言

大氣中的二氧化硫?qū)θ嗣竦纳眢w健康造成了嚴(yán)重危害， 火電站化石能源的使用極大增加了空氣中的二氧化硫的含量， 成為國(guó)內(nèi)外的重點(diǎn)研究對(duì)象[1-4]。 差分吸收光譜技術(shù)(DOAS)具有探測(cè)范圍大、 非接觸、 高靈敏度和可測(cè)量氣體種類多等優(yōu)點(diǎn)。 將測(cè)量得到的光譜數(shù)據(jù)進(jìn)行解釋， 是DOAS技術(shù)的重要內(nèi)容， 然而光譜反演得到的是路徑積分濃度。 為了得到濃度的空間分布， 一種方法是將DOAS技術(shù)與計(jì)算機(jī)斷層技術(shù)(CT)相結(jié)合[5]， 稱之為tom-DOAS技術(shù)。 在該技術(shù)中， 常用兩臺(tái)被動(dòng)DOAS采集光譜， 該條件下的CT重建屬于極端的不完全角度重建。 低三階導(dǎo)數(shù)法(LTD)[6-7]給出了氣體擴(kuò)散的模型作為先驗(yàn)信息， 在某些條件下重建效果更好， 被廣泛用于研究工業(yè)氣體泄露、 電廠煙囪煙羽等氣體空間分布。

為了進(jìn)一步提高重建圖像的抗誤差能力， 本文創(chuàng)新性的使用安光所的IDOAS組成實(shí)驗(yàn)系統(tǒng)采集光譜數(shù)據(jù)， 并提出了一種結(jié)合壓縮感知理論和低三階導(dǎo)數(shù)模型的煙羽斷層重建算法——投影凸函數(shù)集低三階導(dǎo)數(shù)法。 根據(jù)壓縮感知理論， 即使圖像在嚴(yán)重的稀疏角度采樣和噪聲影響的條件下， 圖像采樣減少90%時(shí)， 仍然能重建出高質(zhì)量的圖像[8-9]， 本文通過(guò)數(shù)值模擬， 證實(shí)了算法在稀疏采樣和噪聲影響下對(duì)重建質(zhì)量的改善。 進(jìn)行了外場(chǎng)實(shí)驗(yàn)， 并重建了氣體的空間分布。 本文介紹的方法提高了數(shù)據(jù)采集的時(shí)間分辨率， 改善了光譜數(shù)據(jù)CT重建的質(zhì)量， 擴(kuò)展了成像差分吸收光譜技術(shù)的應(yīng)用范圍。

1 DOAS反演原理及數(shù)據(jù)斷層采集系統(tǒng)

1.1 測(cè)量原理

被動(dòng)DOAS的測(cè)量服從朗伯比爾定律

(1)

(2)

將IDOAS接收到的光譜中的寬帶吸收特征通過(guò)多項(xiàng)式擬合去除之后， 即得到表示窄帶吸收特征的差分吸收光譜， 利用氣體分子的差分吸收截面擬合差分吸收光譜， 即可得到Sm[1]。

1.2 斷層數(shù)據(jù)采集系統(tǒng)

基于IDOAS的煙羽斷層掃描系統(tǒng)如圖1所示。 假設(shè)風(fēng)向垂直于紙面。 以其中一臺(tái)IDOAS為原點(diǎn)， 水平方向?yàn)閤軸， 垂直于地面作為y軸建立坐標(biāo)系。 在IDOAS的視場(chǎng)角范圍內(nèi)， 僅需要2 s即可完成48條光譜的采集和存儲(chǔ)， MAX-DOAS采集相同數(shù)量的光譜則需要約5 min。 基于IDOAS的數(shù)據(jù)采集系統(tǒng)大大提高了采集煙羽的時(shí)間分辨率。

圖1 地基IDOAS斷層掃描系統(tǒng)

2 重建算法

將式(2)按照?qǐng)D1離散化并寫成矩陣的形式為

S=HC

(3)

式(3)中，S∈RM表示反演光譜數(shù)據(jù)得到的氣體路徑積分濃度。H∈RM×N的元素表示光路穿過(guò)圖1中的像素的長(zhǎng)度。C∈RN表示圖1中重建區(qū)域的離散區(qū)域的像素， 其數(shù)量為N， 通常NM。

2.1 傳統(tǒng)LTD法

氣體濃度的三階導(dǎo)數(shù)為[6]

3c(k,l)-c(k-1,l)

(4)

3c(k,l)-c(k,l-1)

(5)

式(4)和式(5)當(dāng)中的c表示氣體濃度矩陣， 式(3)中的C是通過(guò)矩陣c重新排列而成的向量。 傳統(tǒng)的LTD法認(rèn)為式(4)和式(5)的值為0

0=LC

(6)

相當(dāng)于假設(shè)氣體在整個(gè)空間中嚴(yán)格的按照二階多項(xiàng)式分布， 會(huì)引起圖像邊緣的劇烈振蕩而產(chǎn)生大量偽影。 為了削弱式(6)產(chǎn)生的偽峰， 將式(6)與式(3)用以下方式聯(lián)擬得到

Sα=HαC

(7)

式(7)中

(8)

式(8)中的α表示相鄰像素符合LTD模型的強(qiáng)度， 可根據(jù)具體情況調(diào)整。 后面的數(shù)值模擬中， 取α=0.1。 最小二乘解為

(9)

2.2 POCS-LTD算法

壓縮感知理論認(rèn)為， 當(dāng)信號(hào)在某種變換下具有稀疏性時(shí)， 只需要對(duì)原始信號(hào)進(jìn)行少量的隨機(jī)采樣就可以精確恢復(fù)出原始信號(hào)。 根據(jù)傳統(tǒng)的LTD法， 認(rèn)為氣體濃度的三階導(dǎo)數(shù)值不是全部為0， 而是大多數(shù)為0， 即氣體濃度的三階導(dǎo)數(shù)是稀疏的。 基于以上分析， 設(shè)計(jì)POCS-LTD算法重建氣體的濃度分布。

首先使用ART算法進(jìn)行數(shù)據(jù)的初始化， 初始化后的濃度寫作CART。 根據(jù)式(4)和式(5)， 全變分‖C‖TV為

(10)

‖C‖TV的梯度‖C‖TV為

(11)

假設(shè)第n-1步重建結(jié)果用Cn-1表示， 第n步用Cn表示， 第n+1步用Cn+1表示， 則氣體濃度分布的迭代公式為

Cn+1=Cn-γdnpn

(12)

式(12)中，γ表示松弛因子。 數(shù)值模擬表明γ取0.2附近的值時(shí)取得最佳的收斂效果。dn表示步長(zhǎng)

(13)

pn表示式(12)中第n步迭代的下降方向， 使用最速下降法

(14)

式(14)中的g通過(guò)式(11)計(jì)算。 當(dāng)相鄰兩次重建均方差小于閾值σresidual時(shí)迭代停止

(15)

σresidual設(shè)置為10-10。 有時(shí)收斂不到設(shè)置的閾值， 因此設(shè)置最大迭代次數(shù)作為另一個(gè)迭代停止條件。 數(shù)值模擬表明， 迭代次數(shù)超過(guò)400次之后， 算法收斂速度緩慢， 所以設(shè)置算法迭代超過(guò)400次之后自動(dòng)停止。

3 POCS-LTD算法數(shù)值模擬

使用高斯模型描述氣體的濃度分布， 并將峰值濃度歸一化。 在氣體分布重建中， 使用接近度σnearness痕量重建質(zhì)量

(16)

圖2 氣體濃度分布圖像

圖2中(a)， (b)， (c)代表氣體的真實(shí)分布。 (d)， (e)， (f)為L(zhǎng)TD法重建的圖像， 接近度為0.598 3， 0.530 8和0.590 9。 (g)， (h)， (i)用POCS-LTD法重建的接近度分別為0.103 3， 0.147 9和0.297 4。 在最佳情況下POCS-LTD算法將接近度減小了83%以上， 在最差情況下也減小了近50%。 IDOAS測(cè)量過(guò)程中， 總的誤差控制在20%以內(nèi)， 通過(guò)給Si疊加隨機(jī)誤差ΔSi的方式模擬測(cè)量誤差

(17)

式(17)中， ΔSi=fSiRrand。f控制誤差的大小，Rrand是方差為1的隨機(jī)數(shù)序列，Rrand∈RM。 圖3給出了重建接近度與誤差系數(shù)f之間關(guān)系的曲線， T-LTD表示用LTD法重建， P-LTD表示用POCS-LTD重建。

圖3 重建接近度隨誤差系數(shù)f變化曲線

可以看到， POCS-LTD法的接近度比傳統(tǒng)的LTD法低得多。 并且誤差系數(shù)f越大優(yōu)勢(shì)越明顯。 這是由于壓縮感知理論中， 約束等距條件把噪聲對(duì)重建結(jié)果影響限制在一定范圍內(nèi)造成的。

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

在淮南某電廠外取得實(shí)驗(yàn)數(shù)據(jù)。 當(dāng)煙羽位于重建圖像幾何中心， 且與兩臺(tái)IDOAS夾角為90°時(shí)重建結(jié)果最好。 因此IDOAS間的距離最好為煙羽高度的2倍。 由于煙囪高210 m， 考慮到煙羽抬升， IDOAS距離為450～550 m時(shí)最符合要求。 在煙囪下風(fēng)處布置掃描系統(tǒng)， 煙羽近似垂直于重建平面。 外場(chǎng)實(shí)驗(yàn)在天氣晴朗的情況下進(jìn)行， 煙羽中的水汽對(duì)測(cè)量值的影響可忽略。

SO2路徑積分濃度的波段為307.5～318 nm， 參與反演的氣體分子包括SO2(293 K)， NO2(294 K)， O3(243 K)， O3(218 K)和ring光譜。 圖4展示了其中一條光譜的反演情況。

圖4 SO2柱濃度IDOAS擬合反演實(shí)例

圖4(a)為測(cè)量光譜擬合實(shí)例， 濃度為1.11×1017molecules·cm-2。 圖4(b)所示的擬合殘差小于0.003 22。 圖5給出了POCS-LTD法重建的污染氣體斷層圖像。

圖5 POCS-LTD法重建SO2氣體濃度分布

在圖5中僅有少量偽影， 不影響對(duì)煙羽的觀測(cè)。 重建圖像的投影值和測(cè)量值如圖6所示。 根據(jù)圖3所示的數(shù)值模擬結(jié)果， 重建接近度約為0.3。

圖6 路徑積分濃度的測(cè)量值與重建圖像的投影的對(duì)比

可以看到圖6(b)中， 重建圖像的投影值和測(cè)量值并不完全符合， 這是由于POCS-LTD算法在投影方程和LTD模型之間取得平衡， 而測(cè)量誤差只能影響投影方程， 卻不能影響LTD模型， 算法利用LTD模型對(duì)誤差進(jìn)行了修正。 用一致性相關(guān)因子來(lái)衡量未知分布?xì)怏w的重建結(jié)果

σCCF=ρA

(18)

式(18)中的ρ為相關(guān)系數(shù)，A是曲線位移校正因子。 圖6中σCCF=0.916 5， 大于0.8， 說(shuō)明重建結(jié)果較好， 不需要重新計(jì)算。

5 結(jié) 論

利用IDOAS搭建了光譜采集系統(tǒng)， 提出了用全變分改進(jìn)LTD的氣體濃度分布模型， 設(shè)計(jì)了POCS-LTD算法， 說(shuō)明了將壓縮感知理論引入氣體重建領(lǐng)域的可行性。 在文中設(shè)置的參數(shù)條件下， 測(cè)量誤差越大， POCS-LTD抑制測(cè)量誤差的優(yōu)勢(shì)越明顯。 外場(chǎng)實(shí)驗(yàn)表明， 該方法能夠重建出清晰的氣體斷層圖像， 重建圖像的投影值與光譜測(cè)量值反演結(jié)果之間的一致性相關(guān)因子大于0.9。 POCS-LTD法能夠修正部分測(cè)量誤差。 本文介紹的方法擴(kuò)展了IDOAS技術(shù)的應(yīng)用范圍。 然而POCS-LTD算法的運(yùn)算復(fù)雜度較大。 LTD氣體擴(kuò)散模型本質(zhì)上是對(duì)高斯模型的近似， 重建結(jié)果的峰值濃度偏低， 重建結(jié)果有變圓的趨勢(shì)。