MAX-DOAS重建NO2竖直平面分布的方法研究
常振1, 钟鸣宇2,*, 苏静明1,2, 司福祺1, 王煜3, 周海金1, 窦科1, 张泉1
1.中国科学院合肥物质科学研究院安徽光学精密机械研究所, 环境光学与技术重点实验室, 安徽 合肥 230031
2.安徽理工大学电气与信息工程学院, 安徽 淮南 232001
3.安徽大学物质科学与信息技术研究院, 信息材料与智能感知安徽省实验室, 安徽 合肥 230601
*通讯作者 e-mail: ap0005238.c.b@163.com

作者简介: 常 振, 1988年生, 中国科学院安徽光学精密机械研究所副研究员 e-mail: zhchang@aiofm.ac.cn

摘要

多轴差分吸收光谱仪(MAX-DOAS)结合计算机断层重建算法可获取目标痕量气体的空间分布情况。 为研究在具有背景浓度的条件下, 如城市背景下某个竖直截面上重建NO2空间分布的可行性, 设计了气体浓度可控条件下的验证性实验; 证明了利用MAX-DOAS在竖直平面重建NO2气体分布的可行性。 将充入标准气体的JGS1石英玻璃样品池作为研究对象, 使用两台MAX-DOAS采集光谱数据。 将气体浓度的梯度作为先验信息, 利用经典的ABOCS算法和Barzilai-Borwein算法重建了竖直平面内的NO2气体分布, 验证了利用MAX-DOAS在竖直平面内重建NO2气体空间分布的可行性, 同时确定了背景浓度对重建结果的影响。 研究结果表明, 以天空为背景的光谱作为参考谱和以空样品池为背景作为参考谱, 反演得到的NO2浓度非常接近, 因此研究对象中的样品池容器在NO2竖直平面分布重建方法中对实验结果的影响可以忽略。 实验中以市区为背景的MAX-DOAS具有较高的背景浓度, 特别是在仰角较低的情况下NO2背景浓度几乎达到6×1016 molec·cm-2, 以城市郊区没有明显的污染源为背景的MAX-DOAS, 背景浓度较低可以忽略。 重建结果显示, 当仰角为28°时, 气体沿光路的平均分子数密度为3.932 7×1015 molec·cm-2, 且在样品池内下部密度大, 上部密度小; 重建得到的SCD和测量得到的SCD符合比较好, 计算结果显示重建得到的气体分子数密度的峰值为5.77×1015 molec·cm-2, 与以城市郊区为背景的MAX-DOAS反演结果较为接近, 而以市区为背景时, 特别是仰角较小时, NO2背景浓度特别明显, 重建结果比测量结果的值小很多。 结果表明, 背景浓度在重建图像中表现为伪影, 影响对气体分布的观察, 而如果在重建算法时加入利用样品池内外气体存在浓度突变这一先验信息, 能够减轻背景浓度对重建结果造成的影响。

关键词: 差分吸收光谱; 浓度重建; 多轴差分吸收光谱; 迭代算法; 数据拟合
中图分类号:O433.4 文献标志码:A
Study on the Reconstructing the NO2 Gas Distribution in a Vertical Plane Using MAX-DOAS
CHANG Zhen1, ZHONG Ming-yu2,*, SU Jing-ming1,2, SI Fu-qi1, WANG Yu3, ZHOU Hai-jin1, DOU Ke1, ZHANG Quan1
1. Key Laboratory of Environmental Optical and Technology, Anhui Institute of Optics and Fine Mechanics, Hefei Institutes of Physical Science, Chinese Academy of Sciences, Hefei 230031, China
2. College of Electrical and Information Engineering, Anhui University of Science and Technology, Huainan 232001, China
3. Information Materials and Intelligent Sensing Laboratory of Anhui Province, Institutes of Physical Science and Information Technology, Anhui University, Hefei 230601, China
*Corresponding author
Abstract

The multi-axis differential absorption spectrometer (MAX-DOAS) combined with computed tomography (CT) reconstruction algorithm can be used to obtain the spatial distribution of the target trace gases. In order to study the feasibility of reconstructing the spatial distribution of NO2 on a vertical cross-section under the condition of background concentration (such as the urban background), a confirmatory experiment was designed under the condition of controllable gas density. The feasibility of using MAX-DOAS to reconstruct the distribution of NO2 in the vertical plane is proved. The JGS1 quartz glass sample cell filled with standard gas was used as the research object, and two MAX-DOAS were used to collect spectral data. Taking the gradient of gas density as a priori information, and using the classical ABOCS algorithm and the Barzilai-Borwein algorithm, the NO2 distribution in the vertical plane is reconstructed. The feasibility of using MAX-DOAS to reconstruct the spatial distribution of NO2 in the vertical plane was verified, and the influence of background density on the reconstruction results was determined. The results show that the NO2 density obtained by the retrieving is very close to that by using the sky as the reference spectrum and the empty sample pool as the reference spectrum. Therefore, the influence of the sample pool container on the experimental results can be neglected in the reconstruction method of NO2 vertical plane distribution. The background density of MAX-DOAS with the urban background was high, especially when the observation angle was low, the background density of NO2 was almost 6×1016 molec·cm-2, the MAX-DOAS, whose background density has no obvious pollution source in the city suburb, is low enough to be ignored. The reconstruction results show that the average molecular number density along the optical path is 3.932 7×1015 molec·cm-2 when the observation angle is 28°, and the density in the lower region of the sample pool is higher than that in the upper region. The reconstructed SCD is in good agreement with the measured SCD. The calculated results show that the peak value of molecular number density is 5.77×1015 molec·cm-2, which is close to the MAX-DOAS retrieving results with the background of the urban suburbs. However, it is not close to the MAX-DOAS retrieving results with the background of the urban areas, especially when the elevation angle is small, whose background density of NO2 is especially obvious, and its reconstructed result is much lower than the measured one. To sum up, background density is an artifact in the reconstructed image, which affects the observation of gas distribution. If the prior information of gas density mutation in and out of the sample pool is added to the reconstruction algorithm, the effect of background density on the reconstruction results can be reduced.

Keyword: Differential optical absorption spectrometer; Gas reconstructing; MAX-DOAS; Iterative algorithms; Data fitting
引言

多轴差分吸收光谱(multi-axis differential optical absorption spectroscopy, MAX-DOAS)基于差分吸收光谱技术通过获取目标区域的光谱信息反演痕量气体浓度, 应用广泛, 如测量NO2对流层垂直分布和垂直柱浓度[1]、 基于车载测量城市NO2分布和排放评估、 测量气溶胶时空变化[2]等。 MAX-DOAS的测量结果是沿光路的斜柱浓度(slant column density, SCD)[1, 3], 为了研究其空间分布情况, 可利用计算机断层重建算法(computed tomography, CT)[4, 5, 6]将SCD重建为气体的空间分布, 被称为tom-DOAS技术。 Mettendorf等以石英样品池为研究对象, 利用LP-DOAS(long-path DOAS)在室内实验研究了水平方向上的气体分布情况, 证实了LP-DOAS被用于tom-DOAS技术的可行性[7]

由于SO2和有机物等物质在大气中的含量比较小, 背景浓度对光谱测量及重建结果的影响不大, 因此以往基于MAX-DOAS研究的主要对象为烟囱当中的S O2[8, 9]。 但是NO2作为重要的大气污染物, 对其空间分布的研究同样意义重大。 对于NO2气体的空间分布, 国内外进行了广泛的研究。 Laepple等利用LP-DOAS研究了公路上空NO2的竖直空间分布[10], 结果表明当研究对象远离地面时, LP-DOAS技术不再适用。 另外研究NO2的空间分布时, 由于背景浓度的影响难以消除, 只能在背景浓度较低的条件下进行研究[11]。 然而当空间分布的研究目标为城市上空的NO2时, 背景浓度较高, 其对重建结果的影响无法避免。

为了研究城市背景下某个竖直截面上重建NO2空间分布的可行性, 设计了气体浓度可控条件下的验证性实验。 用JGS1石英玻璃制作样品池, 在样品池中充入标准气体, 通过对比光谱反演的结果, 排除了样品池对实验结果的影响, 并估计了样品池中的气体浓度范围。 利用样品池内外气体存在浓度突变这一先验信息, 作为额外的约束用于重建气体的空间分布, 重建结果说明了背景浓度在气体重建结果中表现为伪影, 先验信息对背景浓度造成的伪影有较强的抑制作用。

1 实验部分
1.1 数据采集系统

实验数据采集设备由两台多轴差分吸收光谱仪(MAX-DOAS)组成; 两台MAX-DOAS分别放置在1 m高的金属支架上, 之间的距离为1.8 m, 镜头在同一直线上。 在两台MAX-DOAS的中间, 放置一个棱长为45 cm的石英玻璃样品池, 样品池被放置在高度为1.4 m的铝合金支架上。 数据采集系统如图1所示。

图1 实验数据采集系统
(a): 立体示意图; (b): 俯视方位图
Fig.1 The sketch of experimental configuration
(a): Dimensional schematic diagram; (b): System top view diagram

图1中, MAX-DOAS-1指向的方向是合肥市的市区, 有一定的NO2背景浓度。 MAX-DOAS-2指向的是合肥市的郊区, 几乎没有污染。 实验现场照片如图2所示。

图2 实验现场Fig.2 Photo of experimental scen

由于实验在室外进行, 石英玻璃对光的吸收和反射影响光谱的采集。 为了检验样品池外壁对测量结果的影响, 实验时首先以天空为背景采集了一组数据, 然后以没有充入气体的空样品池为背景再采集一组数据, 将这两组数据分别作为参考光谱。 然后在样品池中充入标称体积分数为1 508 ppm的样气, 大约5分钟之后, 估计样品池中的气体浓度已经基本恒定时, 开始采集数据。

1.2 测量区域离散化及光路几何

以MAX-DOAS-1的镜头为原点O, 以该原点指向MAX-DOAS-2为横轴X, 以竖直向上为纵轴Y, 以深度方向为Z轴建立三维坐标。 MAX-DOAS的镜头扫描范围为0° ~90° 。 将两台MAX-DOAS之间区域及样品池在XOY面的投影, 进行区域离散化, 分成20× 20个像素, 所得到的区域离散化示意图如图3所示。

图3 区域离散化示意图Fig.3 Discretization of reconstructed region

图3中的红线表示测量光路, 穿过该区域时与相应的像素相交, 该像素到O点的长度用mh, k表示, h, k分别表示第h条测量光路和其中的第k个像素, 不在该测量光路内的像素用0表示, 则所有测量光路中的像素到O的距离可用矩阵H表示

H=m1, 10000m2, 1m2, 2m2, 2000mh, 1mh, 2mh, k00mp, 1mp, 2mp, 2000(1)

矩阵H为中的每行元素对应一条扫描路径, mh, k表示第h条扫描路径下第k个像素的测量光路长度, 总共有p条扫描路径; 最长的扫描路径为第h行, 路径中总共有k个像素, 其他行不足数量的元素位置, 用0表示。 ci表示第i个元素内的气体平均浓度, 所有像素组成列向量c, 则测量得到的所有路径上的斜柱浓度SCD可用S表示为

S=Hc(2)

式(2)中, c中的元素个数远大于扫描路径数, 因此需要使用重建算法进行计算。

2 重建算法

经典的迭代算法需要根据SCD进行迭代, 重建图像的投影值和测量值几乎完全拟合, 因此从理论上来说, 背景浓度所造成的伪影完全无法消除。 然而由于NO2装在样品池中, 和样品池外NO2的分子数密度存在梯度, 该信息作为先验知识用式(3)描述为

CTV=j, ku(k, l)=j, k[c(k, l)-c(k-1, l)]2-[c(k, l)-c(k-1, l-1)]2(3)

式(3)中, C为矩阵c排列成的向量, ‖ CTVC的全变分, c(k, l)表示图3中重建区域内第k行、 第l列元素。 利用经典的ABOCS算法[12]重建气体的浓度分布为

C* =argminCTVs.t. 12HC-S22ε,  C(4)

ε 在此处表示SCD的测量值和投影值之间的差, 式(4)表示最终得到的浓度分布是投影方程和气体分布的先验知识共同决定的。 公式可转化为式(5)

C* =argminCTV+1t-logε-12HC-S22s.t. C0(5)

在经过初始化之后, 利用优化算法, 根据气体浓度梯度全变分的梯度对重建图像进行进一步优化。 全变分的梯度为

g=CTV+HT(HC-S)maxβ, tε-12HC-S22(6)

重建算法的迭代公式为

Cn+1=Cn+αnpn(7)

式(7)中, Cn表示迭代算法第n步的迭代结果, Cn+1表示第n+1步的迭代结果。 α n表示步长, 由Barzilai-Borwein算法确定。 而

pn=gn, gn0Cn> 00, 其他(8)

当迭代超过一定次数, 或者CnCn+1之间的差值小于某个阈值时迭代停止。

2 结果与讨论
2.1 光谱反演

反演NO2所使用的波段范围为338~370 nm, 吸收分子如表1所示。

表1 参与反演的气体分子及吸收截面和Ring光谱 Table 1 List of trace gas references used for the DOAS analysis

其中一条光谱的反演结果如图4所示。

图4 光谱反演实例
(a): 拟合结果; (b): 拟合残差
Fig.4 Example of spectral inversion
(a): Fitting result; (b): Fitting residuals

图4(a)是光谱拟合的结果, 黑色的谱线表示测量光谱, 红色的谱线表示拟合光谱, 反演得NO2的SCD为1.54× 1017 molec· cm-2。 图4(b)表示反演之后的残差, 大小为1.68× 10-3。 MAX-DOAS-1和MAX-DOAS-2测量得到的SCD如图5所示。

图5 反演得到的SCD
(a): MAX-DOAS-1的测量值; (b): MAX-DOAS-2的测量值
Fig.5 SCD results of inversion of sepctra taken
(a): MAX-DOAS-1; (b): MAX-DOAS-2

从图5中可以看出, 不管是将以天空为背景的光谱作为参考谱, 还是以空样品池为背景作为参考谱, 反演得到的NO2浓度非常接近, 因此认为石英玻璃样品池对实验结果没有影响。 由于MAX-DOAS-1正对的是合肥市区, 特别是在仰角较低的情况下背景浓度几乎达到6× 1016 molec· cm-2。 而MAX-DOAS-2正对合肥郊区, 没有明显的污染源, 观测不到NO2存在, 这种污染分布情况在合肥秋冬季中等污染情况下具有一定的代表性。 另外需要注意图5(a)中仰角较低时展现出的背景浓度对重建结果的影响。

2.2 样品池中NO2浓度估计

在估计样品池中NO2浓度的过程中, 由于MAX-DOAS-1的测量结果受到背景浓度的影响, 此处主要将MAX-DOAS-2的结果作为主要研究对象。 样气的标称值1 508 mol(mol)-1, 实验当日天气晴朗, 光照条件比较好, 有微风。 气温为18 ℃, 气压为P=101.18 kPa。 假如样品池中的气体体积分数为1 508 mol(mol)-1, 计算SCD的理论值如式(9)

PV=nRT(9)

式(9)中, R=8.31 J· (mol· k)-1为普适气体常数, n=1, 得到1 mol气体的体积为23.900 1 L。 标准气体的标注为: C=1 508× 10-6 mol(mol)-1, NA=6.02× 1023, 为阿伏伽德罗常数, 经换算得到气体的分子数密度为

c0=C×NAV=3.7984×1016molec·cm-3(10)

当仰角为28° 时, 经过样品池的光路长度为Ld=42.21 cm, 气体路径积分浓度为

CSCD=Ldc0=1.6033×1018molec·cm-2(11)

当仰角为44° 时, 光路长度为Ld=47.45 cm, 同理可得路径积分浓度为: 1.802 3× 1018 molec· cm-2

然而如图5(b)所示, 当仰角为28° 时反演得到的SCD为1.66× 1017 molec· cm-2, 只有理论值的1/10, 显示样品池中NO2的体积分数并没有达到理论值的大小。 当仰角为44° 时反演得到的SCD为4.94× 1016 molec· cm-2, 只有理论值的1/36。 由此可见样品池底部的气体体积分数较大, 样品池顶部的体积分数较小。 根据式(11), 当仰角为28° 时, 气体沿光路的平均分子数密度为3.932 7× 1015 molec· cm-3, 且分布为下部分子数密度大, 上部分子数密度小。

2.3 重建结果

样品池中NO2的分布情况如图6所示, 白色方框表示样品池的位置。 由于像素离散化和测量误差, 在图6的左上侧和右上侧有少量伪影。 在样品池的左边有部分的气体分布在样品池之外, 这是画等高线时插值造成的。 样品池中气体表现为下部的分子数密度高, 上部的分子数密度低。 分子数密度的最大值约为5.77× 1015 molec· cm-3, 和前文的估计基本相符。

图6 NO2气体分布重建图Fig.6 Reconstruction of NO2 distribution

在样品池右边显示有较多的气体分布在样品池之外, 在图5中仰角较低时, MAX-DOAS-1的测量值与样品池中的气体浓度无关, 完全取决于背景浓度的影响。 SCD的测量值和重建值如图7所示。

图7 SCD的测量值与重建值
(a): MAX-DOAS-1; (b): MAX-DOAS-2
Fig.7 The measured and reconstructed SCD
(a): MAX-DOAS-1; (b): MAX-DOAS-2

图7显示, 重建得到的SCD和测量得到的SCD符合得比较好, 特别是图7(b)对应的是几乎没有污染的郊区。 而图7(a)的测量结果对应的是污染比较大的市区, 特别是仰角较小时, 污染特别明显。 然而仰角较小时, 重建结果比测量结果的值小很多, 说明引入先验信息之后, 重建算法能部分纠正背景浓度对重建结果的影响。 否则图6中, 样品池右侧的伪影强度会更高。

为了进一步证明MAX-DOAS在竖直平面重建NO2空间分布的可行性, 使用510 mol(mol)-1的标准气体进行了实验。 实验当日天气晴朗, 光照条件比较好, 有微风, 气温为17 ℃, 气压与前述实验相同。 当天空气质量良好, 在低空几乎没有NO2分布。 得到样品池中NO2的空间分布如图8所示。 可见, 在没有低空背景浓度影响的条件下, NO2所在的位置与样品池位置符合得更好。 利用3.2中的方法理论分析得到的仰角28° 时的平均分子数密度为9.828 5× 1014 molec· cm-3, 利用CT算法重建NO2空间分布得到的分子数密度的最大值为3.036 5× 1015 molec· cm-3

图8 NO2气体分布重建图Fig.8 Reconstruction of NO2 distribution

图9给出了SCD的测量值和重建值对比。

图9 SCD的测量值与重建值
(a): MAX-DOAS-1; (b): MAX-DOAS-2
Fig.9 The measured and reconstructed SCD
(a): MAX-DOAS-1; (b): MAX-DOAS-2

3 结论

研究了利用MAX-DOAS在竖直平面重建NO2空间分布的可行性。 利用JGS1石英玻璃样品池装入了NO2气体作为研究对象, 反演结果显示样品池本身不会对测量结果产生影响。 当仰角为28° 时, 气体沿光路的平均分子数密度为3.932 7× 1015 molec· cm-2, 且在样品池内下部密度大, 上部密度小。 重建得到的气体分子数密度的峰值为5.77× 1015 molec· cm-2, 与预先的估计比较接近。 研究表明, 背景浓度在重建图像中表现为伪影, 影响对气体分布的观察。 如果在重建算法时加入先验信息, 能够减轻背景浓度对重建结果造成的影响。

参考文献
[1] WANG Yang, LI Ang, XIE Pin-hua, et al(王杨, 李昂, 谢品华, ). Acta Physica Sinica(物理学报), 2013, 62(20): 200705. [本文引用:2]
[2] LIU Hao-ran, HU Qi-hou, TAN Wei, et al(刘浩然, 胡启后, 谈伟, ). Spectroscopy and Spectral Analysis(光谱学与光谱分析), 2021, 41(1): 11. [本文引用:1]
[3] Cheng X H, Ma J Z, Jin J L, et al. Atmospheric Chemistry & Physics, 2020, 20: 10757. [本文引用:1]
[4] Hartl A, Song B C, Pundt I. Atmospheric Chemistry and Physics, 2006, 6(3): 847. [本文引用:1]
[5] Frins E, Bobrowski N, Platt U, et al. Appl. Opt. , 2006, 45(24): 6227. [本文引用:1]
[6] Pundt I, Kai U M. Appl. Opt. , 2005, 44(23): 4985. [本文引用:1]
[7] Mettendorf K U, Hartl A, Pundt I. Journal of Environmental Monitoring, 2006, 8: 279. [本文引用:1]
[8] ZHONG Ming-yu, XI Liang, SI Fu-qi, et al(钟鸣宇, 奚亮, 司福祺, ). Acta Physica Sinica(物理学报), 2019, 68(16): 164205. [本文引用:1]
[9] Olaguer E P, Stutz J, Erickson M H, et al. Atmospheric Environment, 2017, 150(2): 220. [本文引用:1]
[10] Laepple T, Knab V, Mettendorf K U, et al. Atmospheric Chemistry and Physics, 2004, 4(5): 1323. [本文引用:1]
[11] Frins E, Bobrowski N, Osorio M, et al. Atmospheric Environment, 2014, 98: 347. [本文引用:1]
[12] Niu Tianye, Zhu Lei. Accelerated Barrier Optimization Compressed Sensing (ABOCS) Reconstruction: Performance Evaluation for Cone-beam CT, 2012 Nuclear Science Symposium & Medical Imaging Conference. IEEE, 2013: 4588(doi: 10.1109/NSSMIC.2012.6551535). [本文引用:1]