基于成像差分吸收光谱技术和压缩感知理论的烟囱烟羽断层重建研究

引用本文

钟鸣宇, 周海金, 司福祺, 王煜, 窦科, 苏静明. 基于成像差分吸收光谱技术和压缩感知理论的烟囱烟羽断层重建研究[J]. 光谱学与光谱分析, 2021,41(6): 1708-1712.
ZHONG Ming-yu, ZHOU Hai-jin, SI Fu-qi, WANG Yu, DOU Ke, SU Jing-ming. Reconstruction of Stack Plume Based on Imaging Differential Absorption Spectroscopy and Compressed Sensing[J]. Spectroscopy and Spectral Analysis, 2021,41(6): 1708-1712.
Doi:10.3964/j.issn.1000-0593(2021)06-1708-05 复制到剪切板

Permissions

《光谱学与光谱分析》期刊社所有

基于成像差分吸收光谱技术和压缩感知理论的烟囱烟羽断层重建研究

钟鸣宇^1,^2,³, 周海金², 司福祺^2,^*, 王煜², 窦科², 苏静明^1,^2,³

1.安徽理工大学电气与信息工程学院, 安徽淮南 232001

2.中国科学院安徽光学精密机械研究所环境光学与技术重点实验室, 安徽合肥 230031

3.中国科学技术大学研究生院科学岛分院, 安徽合肥 230031

*通讯作者 e-mail: sifuqi@aiofm.ac.cn

作者简介: 钟鸣宇, 1982年生, 安徽理工大学电气与信息工程学院讲师 e-mail: ap0005238.c.b@163.com

收稿日期: 2020-05-08 修回日期: 2020-09-23

基金: 国家重点研发计划项目(2016YFC0200400, 2017YFB0503901)资助

摘要

烟羽断层重建质量受两方面条件限制: 其中一个限制条件是遥感设备的时间分辨率。以往的研究多使用多轴差分吸收光谱仪(MAX-DOAS)进行CT重建, 受采集数据速度的限制, 重建图像的时间分辨率较低。另一个限制条件是, 采集到的数据量有限, 是典型的不完全角度重建。过去多使用代数迭代重建算法或统计迭代重建算法, 重建图像受测量误差的影响比较大, 分辨率较低且伪影较多。构造了基于成像差分吸收光谱技术(IDOAS)的光谱数据采集系统, 与多轴差分吸收光谱仪构造的系统相比, 数据采集的时间分辨率提高了160多倍, 基本解决了时间分辨率的问题。提出了一种基于压缩感知理论和低三阶导数模型的烟羽断层重建算法——投影凸函数集低三阶导数法, 简称为POCS-LTD。在投影的过程中, 使用代数重建算法使重建图像符合投影方程; 在全变分迭代的过程中使用了优化算法, 将低三阶导数模型的全变分归一化值作为优化算法的迭代方向, 前次迭代运算结果与本次投影运算的差值的模作为迭代步长。对重建算法进行了数值模拟, 并以重建图像的接近度和一致性相关因子为指标, 对重建结果进行了分析。数值模拟表明, 算法具有良好的抗误差能力, 与传统的低三阶导数法相比, 本文提出的算法将重建接近度减小了80%以上。使用烟羽数据采集系统进行了外场实验, 用POCS-LTD算法对外场实验的数据进行了烟羽重建, 重建图像显示烟羽图像清晰, 伪影得到了较好的抑制。介绍的烟羽断层数据采集系统和烟羽断层重建算法, 提高了烟羽断层重建图像的时间分辨率, 减少了重建图像的伪影, 扩大了光谱测量技术的应用范围。

关键词: 成像差分吸收光谱仪; 大气吸收光谱; 计算机断层重建; 压缩感知

中图分类号:O433.1 文献标志码:A

Reconstruction of Stack Plume Based on Imaging Differential Absorption Spectroscopy and Compressed Sensing

ZHONG Ming-yu^1,^2,³, ZHOU Hai-jin², SI Fu-qi^2,^*, WANG Yu², DOU Ke², SU Jing-ming^1,^2,³

1. College of Electrical and Information Engineering, Anhui University of Science and Technology, Huainan 232001, China

2. Anhui Institute of Optics and Fine Mechanics, Chinese Academy of Sciences, Hefei 230031, China

3. University of Science and Technology of China, Hefei 230031, China

*Corresponding author

Abstract

The quality of computed tomography image of stack plume has two limits. Oneof the limitations isthe temporal resolution of data acquisition of remote sensing instrument. The conventional remote sensing instrument used for tomography is an multi-axis differential optical absorption spectrometer. Limited by its’ speed of data acquisition, the temporal resolution of the reconstructed image is low. The other limitation is the insufficient data acquired from remote sensing instruments. Algebraic reconstruction algorithms and iterative statistical algorithms are usually used in the reconstruction, but the reconstructed images have a low resolution and plenty of artifacts. To overcome the first limitation, data acquisition system in this paper is composed of imaging differential optical absorption spectrum technique. Compared with the system composed of the multi-axis differential optical absorption spectrum, the system’s temporal resolution increases above 160 times. An algorithm based on compressed sensing and a low third derivative model is introduced to overcome the second limitation. The algorithm is called projection on convex sets-low the third derivative method, which is POCS-LTD in short. The proposed algorithm belongs to gradient projection for sparse reconstruction algorithms, which is divide into two steps: projection and total variation iteration. In the process of projection, the algebraic reconstruction algorithm is used to make the reconstructed image conform to the projection equation, and the optimization algorithm is used in the process of total variation iteration. In the process of the total variation, the normalized value of the low third order derivative model is used as the iterative direction of the optimization algorithm, and the module of the difference between the result of the previous iteration and the present projection operation is used as the iterative step. According to nearness, the reconstructed images are evaluated, the relative difference of maximum and concordance correlation factor by numerical simulation. The numerical simulation shows that the proposed algorithm has good error resist,and reduces nearness by 80% compared with the traditional low third derivative method. The method is used to reconstruct the plume by the data gets from the field campaign. A clear plume and suppression of artifacts can be seen from the reconstructed image. The data acquisition system and algorithm introduced in this paper promote temporal resolution and reduce artifacts of the reconstructed images, and improve the technique’s practicability.

Keyword: Imaging differential optical absorption spectrometer; Spectral absorption by atmospheric gases; Computed tomography; Compressed sensing

文章图片

引言

大气中的二氧化硫对人民的身体健康造成了严重危害, 火电站化石能源的使用极大增加了空气中的二氧化硫的含量, 成为国内外的重点研究对象^{[1, 2, 3, 4]}。差分吸收光谱技术(DOAS)具有探测范围大、非接触、高灵敏度和可测量气体种类多等优点。将测量得到的光谱数据进行解释, 是DOAS技术的重要内容, 然而光谱反演得到的是路径积分浓度。为了得到浓度的空间分布, 一种方法是将DOAS技术与计算机断层技术(CT)相结合^[5], 称之为tom-DOAS技术。在该技术中, 常用两台被动DOAS采集光谱, 该条件下的CT重建属于极端的不完全角度重建。低三阶导数法(LTD)^{[6, 7]}给出了气体扩散的模型作为先验信息, 在某些条件下重建效果更好, 被广泛用于研究工业气体泄露、电厂烟囱烟羽等气体空间分布。

为了进一步提高重建图像的抗误差能力, 本文创新性的使用安光所的IDOAS组成实验系统采集光谱数据, 并提出了一种结合压缩感知理论和低三阶导数模型的烟羽断层重建算法— — 投影凸函数集低三阶导数法。根据压缩感知理论, 即使图像在严重的稀疏角度采样和噪声影响的条件下, 图像采样减少90%时, 仍然能重建出高质量的图像^{[8, 9]}, 本文通过数值模拟, 证实了算法在稀疏采样和噪声影响下对重建质量的改善。进行了外场实验, 并重建了气体的空间分布。本文介绍的方法提高了数据采集的时间分辨率, 改善了光谱数据CT重建的质量, 扩展了成像差分吸收光谱技术的应用范围。

1 DOAS反演原理及数据断层采集系统

1.1 测量原理

被动DOAS的测量服从朗伯比尔定律

$\begin{array}{l} I (λ) = I_{0} (λ) \exp {- \overset{n}{\sum_{m = 1}} [σ_{m} (λ) + σ'_{m} (λ) + \\ ε_{R} (λ) + ε_{M} (λ)] S_{m}} + B (λ) \end{array}$ (1)

式(1)中, I₀(λ )为光源发出的光强, 被动DOAS中常使用天顶光来代替。 I(λ )是IDOAS接收到的光强。 σ _m(λ )为第m种气体的窄带吸收截面; σ '_m(λ )为第m种气体的宽带吸收截面; ε _R(λ )为瑞利散射, ε _M(λ )为米散射, B(λ )表示各种噪声。 S_m表示该种气体的斜柱浓度, 等于气体的浓度c_m对光程r积分, 积分距离为L。

$S_{m} = \overset{L}{\int_{0}} c_{m} (r) d r$ (2)

将IDOAS接收到的光谱中的宽带吸收特征通过多项式拟合去除之后, 即得到表示窄带吸收特征的差分吸收光谱, 利用气体分子的差分吸收截面拟合差分吸收光谱, 即可得到 ${S_{m}}^{[1]}$ 。

1.2 断层数据采集系统

基于IDOAS的烟羽断层扫描系统如图1所示。假设风向垂直于纸面。以其中一台IDOAS为原点, 水平方向为x轴, 垂直于地面作为y轴建立坐标系。在IDOAS的视场角范围内, 仅需要2 s即可完成48条光谱的采集和存储, MAX-DOAS采集相同数量的光谱则需要约5 min。基于IDOAS的数据采集系统大大提高了采集烟羽的时间分辨率。

	Figure Option View Download New Window
	图1 地基IDOAS断层扫描系统Fig.1 Ground based IDOAS tomography system

2 重建算法

将式(2)按照图1离散化并写成矩阵的形式为

$S = HC (3)$

式(3)中, S∈ R^M表示反演光谱数据得到的气体路径积分浓度。 H∈ R^M^×^N的元素表示光路穿过图1中的像素的长度。 C∈ R^N表示图1中重建区域的离散区域的像素, 其数量为N, 通常N≽M。

2.1 传统LTD法

气体浓度的三阶导数为^[6]

$\frac{d^{3} c (k, l)}{d k^{3}} = c (k + 2, l) - 3 c (k + 1, l) + 3 c (k, l) - c (k - 1, l)$ (4)

$\frac{d^{3} c (k, l)}{d l^{3}} = c (k, l + 2) - 3 c (k, l + 1) + 3 c (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)中

$S_{α} = [\begin{array}{l} S \\ 0 \end{array}], H_{α} = [\begin{array}{l} H \\ αL \end{array}]$ (8)

式(8)中的α 表示相邻像素符合LTD模型的强度, 可根据具体情况调整。后面的数值模拟中, 取α =0.1。最小二乘解为

$C = \{\begin{array}{l} ({H^{T}}_{α} H_{α})^{- 1} {H^{T}}_{α} S, & C > 0 \\ 0, & C < 0 \end{array}$ (9)

2.2 POCS-LTD算法

压缩感知理论认为, 当信号在某种变换下具有稀疏性时, 只需要对原始信号进行少量的随机采样就可以精确恢复出原始信号。根据传统的LTD法, 认为气体浓度的三阶导数值不是全部为0, 而是大多数为0, 即气体浓度的三阶导数是稀疏的。基于以上分析, 设计POCS-LTD算法重建气体的浓度分布。

首先使用ART算法进行数据的初始化, 初始化后的浓度写作C_ART。根据式(4)和式(5), 全变分‖ C‖ _TV为

$\sum_{k, l} \sqrt{\left[\frac{\mathrm{d}^{3} \boldsymbol{c}(k, l)}{\mathrm{d} k^{3}}\right]^{2}+\left[\frac{\mathrm{d}^{3} \boldsymbol{c}(k, l)}{\mathrm{d} l^{3}}\right]^{2}}$(10)

‖ C‖ _TV的梯度∇‖ C‖ _TV为

$g = \nabla ‖ C ‖_{TV} = \frac{1}{u (k - 2, l)} \frac{d^{3} c (k - 2, l)}{d k^{3}} - \frac{3}{u (k - 1, l)} \frac{d^{3} c (k - 1, l)}{d k^{3}} + \frac{3}{u (k, l)} \frac{d^{3} c (k, l)}{d k^{3}} - \frac{1}{u (k + 1, l)} \frac{d^{3} c (k + 1, l)}{d k^{3}} + \frac{1}{u (k, l - 2)} \frac{d^{3} c (k, l - 2)}{d l^{3}} - \frac{3}{u (k, l - 1)} \frac{d^{3} c (k, l - 1)}{d l^{3}} + \frac{3}{u (k, l)} \frac{d^{3} c (k, l)}{d l^{3}} - \frac{1}{u (k + 1, l)} \frac{d^{3} c (k, l + 1)}{d l^{3}}$ (11)

假设第n-1步重建结果用Cⁿ^-1表示, 第n步用Cⁿ表示, 第n+1步用Cⁿ⁺¹表示, 则气体浓度分布的迭代公式为

$C^{n + 1} = C^{n} - γ d^{n} p^{n}$ (12)

式(12)中, γ 表示松弛因子。数值模拟表明γ 取0.2附近的值时取得最佳的收敛效果。 dⁿ表示步长

$\boldsymbol{d}^{n}=\sqrt{\left(\boldsymbol{C}^{n}-\boldsymbol{C}_{\mathrm{ART}}^{n+1}\right)^{\mathrm{T}}\left(\boldsymbol{C}^{n}-\boldsymbol{C}_{\mathrm{ART}}^{n+1}\right)}$(13)

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

$p^{n} = \frac{g^{n}}{| g^{n} |}$ (14)

式(14)中的g通过式(11)计算。当相邻两次重建均方差小于阈值σ _residual时迭代停止

$σ_{residual} = \frac{1}{N} \overset{N}{\sum_{j = 1}} (C_{j}^{n} - C_{j}^{n - 1})^{2}$ (15)

σ _residual设置为10^-10。有时收敛不到设置的阈值, 因此设置最大迭代次数作为另一个迭代停止条件。数值模拟表明, 迭代次数超过400次之后, 算法收敛速度缓慢, 所以设置算法迭代超过400次之后自动停止。

3 POCS-LTD算法数值模拟

使用高斯模型描述气体的浓度分布, 并将峰值浓度归一化。在气体分布重建中, 使用接近度σ _nearness痕量重建质量

$\sigma_{\text {nearness }}=\sqrt{\frac{\sum_{j}\left(\boldsymbol{C}^{*}(j)-\boldsymbol{C}(j)\right)^{2}}{\sum_{j}\left(\boldsymbol{C}^{*}(j)-\boldsymbol{C}_{\mathrm{avg}}^{*}\right)^{2}}}$(16)

式(16)中, C^*表示气体的真实分布, C表示气体的重建分布。 $C_{avg}^{*}$ 表示气体浓度的平均值。测试图像和重建图像如图2所示。

	Figure Option View Download New Window
	图2 气体浓度分布图像 (a), (b), (c): 测试图像; (d), (e), (f): LTD法重建图像; (g), (h), (i): POCS-LTD法重建图像Fig.2 Gas distribution images (a), (b), (c): Test distribution; (d), (e), (f): Reconstruction of distribution using LTD algorithm; (g), (h), (i): Reconstruction of distribution using POCS-LTD algorithm

图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%以内, 通过给S_i叠加随机误差Δ S_i的方式模拟测量误差

$S'_{i} = S_{i} + Δ S_{i}, i = 1, 2, \dots, M$ (17)

式(17)中, Δ S_i=fS_iR_rand 。 f控制误差的大小, R_rand是方差为1的随机数序列, R_rand∈ R^M。图3给出了重建接近度与误差系数f之间关系的曲线, T-LTD表示用LTD法重建, P-LTD表示用POCS-LTD重建。

	Figure Option View Download New Window
	图3 重建接近度随误差系数f变化曲线Fig.3 Nearness as the function of error factor f

可以看到, POCS-LTD法的接近度比传统的LTD法低得多。并且误差系数f越大优势越明显。这是由于压缩感知理论中, 约束等距条件把噪声对重建结果影响限制在一定范围内造成的。

4 外场实验

在淮南某电厂外取得实验数据。当烟羽位于重建图像几何中心, 且与两台IDOAS夹角为90° 时重建结果最好。因此IDOAS间的距离最好为烟羽高度的2倍。由于烟囱高210 m, 考虑到烟羽抬升, IDOAS距离为450~550 m时最符合要求。在烟囱下风处布置扫描系统, 烟羽近似垂直于重建平面。外场实验在天气晴朗的情况下进行, 烟羽中的水汽对测量值的影响可忽略。

SO₂路径积分浓度的波段为307.5~318 nm, 参与反演的气体分子包括SO₂(293 K), NO₂(294 K), O₃(243 K), O₃(218 K)和ring光谱。图4展示了其中一条光谱的反演情况。

	Figure Option View Download New Window
	图4 SO₂柱浓度IDOAS拟合反演实例 (a): SO₂柱浓度; (b): 拟合残差Fig.4 An example of SO₂ SCD from IDOAS retrieval (a): SO₂ SCD; (b): Fitting residual

图4(a)为测量光谱拟合实例, 浓度为1.11× 10¹⁷ molecules· cm^-2。图4(b)所示的拟合残差小于0.003 22。图5给出了POCS-LTD法重建的污染气体断层图像。

	Figure Option View Download New Window
	图5 POCS-LTD法重建SO₂气体浓度分布Fig.5 Reconstruction of SO₂ of stack plume by using POCS-LTD

在图5中仅有少量伪影, 不影响对烟羽的观测。重建图像的投影值和测量值如图6所示。根据图3所示的数值模拟结果, 重建接近度约为0.3。

	Figure Option View Download New Window
	图6 路径积分浓度的测量值与重建图像的投影的对比 (a): IDOAS #1; (b): IDOAS #2Fig.6 Comparison between measured path integrated concentration and reconstructed path integrated concentration (a): IDOAS #1; (b): IDOAS #2

可以看到图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气体扩散模型本质上是对高斯模型的近似, 重建结果的峰值浓度偏低, 重建结果有变圆的趋势。

参考文献

文献列表

[1]	LIU Jin, SI Fu-qi, ZHOU Hai-jin, et al(刘进, 司福祺, 周海金, 等). Acta Optica Sinica(光学学报), 2015, 35(6): 0630003. [本文引用:1]
[2]	Frins E, Bobrowski N, Osorio M, et al. Atmospheric Environment, 2014, 98: 347. [本文引用:1]
[3]	Lee H, Ryu J, Jeong U, et al. Bull. Korean Chem. Soc. , 2014, 35(12): 3427. [本文引用:1]
[4]	Yao L, Garmash O, Bianchi F, et al. Science, 2018, 361(6399): 278-281. [本文引用:1]
[5]	WEI Min-hong, TONG Min-ming, LI Su-wen, et al(韦民红, 童敏明, 李素文, 等). Spectroscopy and Spectral Analysis(光谱学与光谱分析), 2015, 35(8): 2252. [本文引用:1]
[6]	Price P N, Fischer M L, Gadgil A J, et al. Atmospheric Environment, 2001, 35(16): 2827. [本文引用:2]
[7]	Casaballe N, Osorio M, Martino M D, et al. Earth & Space Science, 2017, 4: 723. [本文引用:1]
[8]	Donoho D L. IEEE Transaction on Information Theory, 2006, 52(4): 1289. [本文引用:1]
[9]	Niu T, Zhu L. Medical Physics, 2012, 39(7): 4588. [本文引用:1]

2015

0.0

... 引言大气中的二氧化硫对人民的身体健康造成了严重危害, 火电站化石能源的使用极大增加了空气中的二氧化硫的含量, 成为国内外的重点研究对象^[1,2,3,4] ...

2014

0.0

2014

0.0

2018

0.0

2015

0.0

... 为了得到浓度的空间分布, 一种方法是将DOAS技术与计算机断层技术(CT)相结合^[5], 称之为tom-DOAS技术 ...

2001

0.0

... 低三阶导数法(LTD)^[6,7]给出了气体扩散的模型作为先验信息, 在某些条件下重建效果更好, 被广泛用于研究工业气体泄露、电厂烟囱烟羽等气体空间分布 ...

... 1 传统LTD法气体浓度的三阶导数为^[6] ...

2017

0.0

2006

0.0

... 根据压缩感知理论, 即使图像在严重的稀疏角度采样和噪声影响的条件下, 图像采样减少90%时, 仍然能重建出高质量的图像^[8,9], 本文通过数值模拟, 证实了算法在稀疏采样和噪声影响下对重建质量的改善 ...

2012

0.0