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

钟鸣宇， 周海金， 司福祺*， 王 煜， 窦 科， 苏静明

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