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

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

1. 安徽理工大學電氣與信息工程學院， 安徽 淮南 232001 2. 中國科學院安徽光學精密機械研究所環境光學與技術重點實驗室， 安徽 合肥 230031 3. 中國科學技術大學研究生院科學島分院， 安徽 合肥 230031

引 言

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

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

1 DOAS反演原理及數據斷層采集系統

1.1 測量原理

被動DOAS的測量服從朗伯比爾定律

(1)

(2)

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

1.2 斷層數據采集系統

基于IDOAS的煙羽斷層掃描系統如圖1所示。 假設風向垂直于紙面。 以其中一臺IDOAS為原點， 水平方向為x軸， 垂直于地面作為y軸建立坐標系。 在IDOAS的視場角范圍內， 僅需要2 s即可完成48條光譜的采集和存儲， MAX-DOAS采集相同數量的光譜則需要約5 min。 基于IDOAS的數據采集系統大大提高了采集煙羽的時間分辨率。

圖1 地基IDOAS斷層掃描系統

2 重建算法

將式(2)按照圖1離散化并寫成矩陣的形式為

S=HC

(3)

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

2.1 傳統LTD法

氣體濃度的三階導數為[6]

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

(4)

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

(5)

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

0=LC

(6)

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

Sα=HαC

(7)

式(7)中

(8)

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

(9)

2.2 POCS-LTD算法

壓縮感知理論認為， 當信號在某種變換下具有稀疏性時， 只需要對原始信號進行少量的隨機采樣就可以精確恢復出原始信號。 根據傳統的LTD法， 認為氣體濃度的三階導數值不是全部為0， 而是大多數為0， 即氣體濃度的三階導數是稀疏的。 基于以上分析， 設計POCS-LTD算法重建氣體的濃度分布。

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

(10)

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

(11)

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

Cn+1=Cn-γdnpn

(12)

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

(13)

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

(14)

式(14)中的g通過式(11)計算。 當相鄰兩次重建均方差小于閾值σresidual時迭代停止

(15)

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

3 POCS-LTD算法數值模擬

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

(16)

圖2 氣體濃度分布圖像

圖2中(a)， (b)， (c)代表氣體的真實分布。 (d)， (e)， (f)為LTD法重建的圖像， 接近度為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測量過程中， 總的誤差控制在20%以內， 通過給Si疊加隨機誤差ΔSi的方式模擬測量誤差

(17)

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

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

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

4 外場實驗

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

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

圖4 SO2柱濃度IDOAS擬合反演實例

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

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

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

圖6 路徑積分濃度的測量值與重建圖像的投影的對比

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

σCCF=ρA

(18)

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

5 結 論

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