摘要:烟羽断层重建一般使用两台光谱仪采集数据, 属于典型的不完全角度重建. 为了提高重建结果的稳定性和接近度, 将压缩感知理论引入气体分布重建领域. 提出了一种新的计算机层析算法——低三阶导数全变分法, 用于重建电厂烟囱排放的SO₂截面的二维分布. 使用低三阶导数模型模拟气体扩散, 认为气体浓度对位置的三阶导数是稀疏的. 将重建图像的全变分作为目标函数, 并通过数值最优化方法求得气体浓度分布的最优解. 数值模拟的结果表明, 与传统的低三阶导数法相比, 低三阶导数全变分法将接近度提高了80%以上. 外场实验表明, 重建图像的一致性相关因子达0.9023. 低三阶导数全变分法能有效消除测量误差对图像重建的影响, 提高重建图像的质量.
关键词: 大气吸收光谱/
烟羽重建/
低三阶导数全变分法/
成像差分吸收光谱仪

English Abstract

全文HTML

--> --> -->

燃煤发电在我国能源结构中占比高达65%, 为经济社会发展和人民生活提供坚实的能源保障的同时, 排放SO₂等也给环境带来了巨大的压力. 最新的研究表明, 大气中的有机物与硫酸在大气新粒子的形成和增长中具有重要作用^[1]. 而燃煤电厂排放的硫超过总量的40%, 因此研究电厂SO₂的排放情况及空间分布具有重要意义.
差分吸收光谱仪(differential optical absorption spectroscopy, DOAS)仅能测得气体浓度沿光传输路径的积分, 而断层图像重建使这种大气参数的间接测量成为实用技术^[2], 在大气遥感领域得以广泛应用. 用于断层重建的光学遥感设备通常仅有2台且固定在地面, 气体断层重建是典型的不完全角度重建. 早期的研究通过直接解方程来重建气体分布, 重建图像质量不高^[3-5]. 一种改进方法是对气体建模作为先验信息, 为方程组增加约束. 常用的模型有Drescher 等^[6]提出的高斯模型、Price等^[7]提出的低三阶导数(low third derivative method, LTD)模型和Olaguer等^[8-10]提出的基于流体力学的欧拉方程模型. LTD模型假设气体浓度对位置的三阶导数值全部为零, 相当于隐含地假设了气体浓度严格满足空间位置的二阶多项式. Johansson等^[11]分别将该方法用于重建火山和烟囱烟羽分布. Kazahaya等^[12]将LTD项乘以权重系数后和投影方程相加, 并使用最小二乘法求解, 但重建的结果仍然不够理想. Casaballe等^[13]利用LTD模型, 对投影方程组进行Tikhonov正则化, 并使用cvx优化工具箱进行求解, 取得了很好的数值模拟效果, 但该方法实际抗误差能力弱, 无法用于外场实验数据. 总的来说, 重建算法还很不成熟, 重建图像存在大量伪影. 压缩感知(compressed sensing, CS)理论为气体分布重建提供了新的思路: 如果能找到气体分布在某种变换下具有稀疏性, 就可以在采样数很少的情况下精确重建气体分布^[14-17].
本文使用成像差分吸收光谱仪(imaging differential optical absorption spectroscopy, IDOAS)采集数据, 提出了一种基于低三阶模型的全变分(low third derivative total variation, LTD-TV)法重建烟羽分布. 该方法基于LTD模型和压缩感知理论对烟羽分布进行重建, 首先使用代数重建算法(algebraic reconstruction technique, ART)对重建数据进行初始化, 用对数障碍函数法^[17]确定目标函数, 全变分法确定下降梯度, 最后使用Barzilai-Borwein(BB)算法^[18]确定步长后对重建结果进行优化. 本文对该方法进行了数值模拟, 并进行了外场实验, 重建了烟囱烟羽的断层2维分布, 并对重建结果进行了分析.

2

-->

2.1.测量原理

DOAS技术基于Lambert-Beer定律

$I(\lambda ) = {I_0}(\lambda )\exp \left[ { - \sum\limits_{m = 1}^n {{\sigma _m}\left( \lambda \right) \cdot {S_{m{\kern 1pt} }}} } \right],$

其中${I_0}\left( \lambda \right)$为大气层外的太阳光强, $I\left( \lambda \right)$为经过大气层后, DOAS接收到的太阳光强, ${\sigma _m}\left( \lambda \right)$为第m种气体的吸收截面,

${S_m} = \int\nolimits_0^L {{c_m}\left( r \right)} {\rm{d}}r,$

其中${S_m}$表示该种气体的斜柱浓度, 等于气体的浓度${c_m}$对光程r积分, 积分距离为L. 通过最小二乘法求解(1)式, 可以得到污染气体的斜柱浓度^[19,20].
2 -->

2.2.断层扫描系统

如图1所示, 扫描烟羽的过程中, 假设风沿Z轴方向, 虚拟的扫描平面垂直于大地, 烟羽被该虚拟平面截取了一个平面, 两台IDOAS放置在扫描平面与大地的交线上. 将IDOAS#1指向IDOAS#2的方向设为X轴的正方向, 竖直向上设为Y轴正方向. 扫描区域离散化如图1中的虚线所示. 两台IDOAS的视场角为30°, 相邻扫描线间隔0.625°, 在2 s内采样并存储48个点的柱浓度数据. 而采用以往常用的多轴差分吸收光谱仪(multi-axis differential optical absorption spectroscopy, MAX-DOAS), 采集相同数量的数据需要5 min以上. IDOAS采集数据的时间分辨率比MAX-DOAS提高了160多倍. 在正式扫描之前先进行预扫描, 调整IDOAS的仰角, 使烟羽位于IDOAS的视场角内. 将扫描区域划分为多个网格, 在这些网格上使用重建算法重建烟羽的浓度分布.
图 1 扫描区域离散化
Figure1. Scanning region discretization.

在离散情况下, (2)式可表示为

$S\left( i \right) = \sum\limits_j {H\left( {i,j} \right)C\left( j \right)} , $

$S\left( i \right)$为第i条射线的路径积分浓度, $C\left( j \right)$表示扫描区域第j个像素中的气体平均浓度, $H\left( {i,j} \right)$表示第i条射线穿越第j个像素的长度. 把整个系统的投影系数和浓度写成矩阵与向量相乘的形式

${S} = {HC},$

其中H表示投影矩阵, C表示图1中重建区域中像素中排列成的列向量, 通常C中像素的数目远大于S中射线的数目, 所以(4)式是一个欠定方程组.
2

-->

3.1.传统LTD算法

为了解方程组(4), LTD模型假设气体浓度相对位置的三阶导数等于0^[7], 于是有

$\begin{split}\frac{{{{\rm{d}}^3}{c}\left( {k,l} \right)}}{{{\rm{d}}{k^3}}} =\, & {c}\left( {k + 2,l} \right) - 3{c}\left( {k + 1,l} \right) \\ & + 3{c}\left( {k,l} \right) - {c}\left( {k - 1,l} \right) = 0, \end{split}$

$\begin{split}\frac{{{{\rm{d}}^3}{c}\left( {k,l} \right)}}{{{\rm{d}}{l^3}}} =\, & {c}\left( {k,l + 2} \right) - 3{c}\left( {k,l + 1} \right) \\ &+ 3{c}\left( {k,l} \right) - {c}\left( {k,l - 1} \right) = 0,\end{split}$

式中${c}\left( {k,l} \right)$是图1所示重建区域中像素矩阵c的第k行l列像素. 把(5)式和(6)式改写为矩阵的形式

${\bf{0}} = {LC},$

将(4)式和(7)式联立得到一个过定方程组, 求解就得到了浓度C的近似值.
2 -->

3.2.LTD-TV算法

在传统的LTD方法中, 隐含地假设气体的浓度是位置的二阶多项式. 然而, 气体的浓度不可能严格地按二阶多项式分布, (7)式左端直接设为零显然会导致图像像素间约束过强, 图像边缘出现大量伪峰, 严重影响重建气体的实际分布. 而本文提出的LTD-TV算法, 仅要求气体的浓度大体按照二阶多项式分布, 也就是假设气体浓度值在三阶导数下是稀疏的, 该假设显然比(7)式合理得多.
LTD-TV法首先使用ART算法对重建图像进行初始化, 然后使用基于全变分的优化算法优化目标函数. 目标函数为

${{C}^ * } = {\rm{argmin}}{\left\| {C} \right\|_{{\rm{TV}}}},\;\;{\rm s.t.}\,,\frac{1}{2}\left\| {{HC} - {S}} \right\|_2^2 \!\leqslant \!\varepsilon ,\;\; {C} \!\geqslant \! 0, $

式中${{C}^ * }$即为所求的气体浓度分布, ${\left\| {C} \right\|_{{\rm{TV}}}}$表示气体浓度分布的全变分的模, $\varepsilon $和镜头接收的光子数目的泊松分布有关, 这里简单设置为10^–12, 不作进一步的深入讨论. 使用对数障碍函数法^[17], 方程(8)中的约束条件可以写入优化问题中

$\begin{split}{{C}^ * } \!=\, &{\rm{argmin}}\left[\!{{{\left\| {C} \right\|}_{{\rm{TV}}}} \!+\! \frac{1}{t}\left(\!{ \!- \log \left(\!{\varepsilon \!-\! \frac{1}{2}\left\| {{HC} \!-\! {S}} \right\|_2^2}\!\right)}\!\!\right)}\!\right],\\ & {\rm{s}}{\rm{.t}}{\rm{.}}{\kern 1pt} \,{C} \geqslant 0.\end{split}$

如果使用低三阶导数模型, 根据(5)式和(6)式, (9)式中气体分布的全变分${\left\| C \right\|_{{\rm{TV}}}}$表达式为

$\begin{split}& \| {C} \|_{\rm TV} = \sum\limits_{k,l} {u}(k,l) \\ =\; & \sum\limits_{k,l} \Big[\big( {c}\left( {k + 2,l} \right) - 3{c}\left( {k + 1,l} \right) + 3{c}\left( {k,l} \right) \\ & - {c}\left( {k - 1,l} \right) \big)^2 + \big( {c}\left( {k,l + 2} \right) - 3{c}\left( {k,l + 1} \right) \\& + 3{c}\left( {k,l} \right) - {c}\left( {k,l - 1} \right) \big)^2\Big]^{1/2}.\end{split} $

对(9)式右边括号中的部分求导, 得到(9)式的梯度为

${g} = \nabla {\left\| {C} \right\|_{{\rm{TV}}}} + \frac{{{{H}^{\rm{T}}}\left( {{HC} - {S}} \right)}}{{\max \left[ {\beta ,t\left( {\varepsilon - \dfrac{1}{2}\left\| {{HC} - {S}} \right\|_2^2} \right)} \right]}}, $

(11)式中, 全变分的梯度的计算公式为

$\begin{split}& \nabla {\left\| {C} \right\|_{{\rm{TV}}}}= \\ & \frac{{3\left( {{c}\left( {k + 2,l} \right) - 3{c}\left( {k + 1,l} \right) + 3{c}\left( {k,l} \right) - {c}\left( {k - 1,l} \right)} \right) + 3\left( {{c}\left( {k,l + 2} \right) - 3{c}\left( {k,l + 1} \right) + 3{c}\left( {k,l} \right) - {c}\left( {k,l - 1} \right)} \right)}}{{{u}\left( {k,l} \right)}} \\ & + \frac{{{c}\left( {k,l} \right) - 3{c}\left( {k - 1,l} \right) + 3{c}\left( {k - 2,l} \right) - {c}\left( {k - 3,l} \right)}}{{{u}\left( {k - 2,l} \right)}} - \frac{{3\left( {{c}\left( {k + 1,l} \right) - 3{c}\left( {k,l} \right) + 3{c}\left( {k - 1,l} \right) - {c}\left( {k - 2,l} \right)} \right)}}{{{u}\left( {k - 1,l} \right)}} \\ & - \frac{{{c}\left( {k + 3,l} \right) - 3{c}\left( {k + 2,l} \right) + 3{c}\left( {k + 1,l} \right) - {c}\left( {k,l} \right)}}{{{u}\left( {k + 1,l} \right)}} + \frac{{{c}\left( {k,l} \right) - 3{c}\left( {k,l - 1} \right) + 3{c}\left( {k,l - 2} \right) - {c}\left( {k,l - 3} \right)}}{{{u}\left( {k,l - 2} \right)}} \\ & - \frac{{3\left( {{c}\left( {k,l + 1} \right) - 3{c}\left( {k,l} \right) + 3{c}\left( {k,l - 1} \right) - {c}\left( {k,l - 2} \right)} \right)}}{{{u}\left( {k,l - 1} \right)}} - \frac{{{c}\left( {k,l + 3} \right) - 3{c}\left( {k,l + 2} \right) + 3{c}\left( {k,l + 1} \right) - {c}\left( {k,l} \right)}}{{{u}\left( {k,l + 1} \right)}}.\end{split}$

(11)式中第二项的分母取$\beta $和$t\left( {\varepsilon - \dfrac{1}{2}\left\| {{HC} - {S}} \right\|_2^2} \right)$的最大值, 是为了防止$t\left( {\varepsilon - \dfrac{1}{2}\left\| {{HC} - {S}} \right\|_2^2} \right)$等于零时分母为零, 设$\beta $为10^–8. 利用(12)式所表示的梯度, 使用梯度下降法求解(9)式. 假设第n步的气体浓度分布用${{C}_n}$表示, 则第$n + 1$步的气体浓度可表示为${{C}_{n + 1}}$为

${{C}_{n + 1}} = {{C}_n} + {\alpha _n}{{p}_n}, $

使用优化算法求解(13)式, ${\alpha _n}$是优化算法的步长, ${{p}_n}$是优化算法的下降方向.

${{p}_n}\left( j \right)=\left\{ \begin{aligned}& {{g}_n}\left( j \right),\,\quad {{g}_n}\left( j \right) \leqslant 0,\text{或}{{C}_n}\left( j \right) > 0\\& 0,\,\quad \quad \quad \text{其他}\end{aligned} \right., $

式中的${{g}_n}\left( j \right)$由(11)式确定, j和(3)式中j的含义相同, 表示第j个像素, (14)式隐含地给气体浓度添加了非负约束. (13)式中${\alpha _n}$用BB算法得到

${\alpha _n} = \frac{{{s}_{n - 1}^{\rm{T}}{{y}_{n - 1}}}}{{{y}_{n - 1}^{\rm{T}}{{y}_{n - 1}}}}, $

式中${{s}_{n - 1}}{\rm{ = }}{{C}_n} - {{C}_{n - 1}},{{y}_{n - 1}} = {{p}_n} - {{p}_{n - 1}}$, 分别表示前后两次迭代的气体浓度分布差, 和前后两次迭代的梯度差. 当(13)式中两次迭代间的图像差别很小时, 说明算法收敛.

${\sigma _{{\rm{residual}}}} = \frac{1}{N}\sum\limits_{j = 1}^N {{{s}_{n - 1}}(j)} , $

式中N表示图像中像素的总数. 当${\sigma _{{\rm{residual}}}}$小于某个阈值, 或者迭代次数大于某个数时迭代停止.

假设真实的气体浓度分布服从高斯模型, 扫描系统如图1所示, 两台IDOAS分别位于X轴的左右两端, Y轴代表气体所在的高度, 扫描线间隔为0.625°, 重建图像解析度为20 × 20. 为了比较方便, 将高斯函数的浓度归一化为1. 在知道气体真实分布的情况下, 使用接近度${\sigma _{{\rm{nearness}}}}$作为重建效果的指标

${\sigma _{{\rm{nearness}}}} = \sqrt {\dfrac{{\displaystyle\sum\nolimits_j {{{\left( {{{C}^ * }\left( j \right) - {C}\left( j \right)} \right)}^2}} }}{{\displaystyle\sum\nolimits_j {{{\left( {{{C}^ * }\left( j \right) - {C}_{{\rm{avg}}}^*} \right)}^2}} }}} , $

其中${{C}^ * }\left( j \right)$是测试图像第j个像素中的气体浓度, ${C}\left( j \right)$是重建图像第j个像素的气体浓度, ${C}_{{\rm{avg}}}^*$是所有像素气体浓度的平均值. 接近度越小, 重建效果越好.
2

-->

4.1.LTD-TV算法与传统LTD算法重建结果对比

图2给出几幅用传统的LTD算法和LTD-TV算法重建气体分布的数值模拟等高线图. 图2(a)—图2(c)分别是单高斯、双高斯、三高斯气体扩散模型.
图 2 传统LTD法与LTD-TV法比较　(a), (b), (c)测试图形; (d), (e), (f)传统LTD法重建图形; (g), (h), (i) LTD-TV法20000次迭代重建图形
Figure2. Comparison between traditional LTD algorithm and LTD-TV algorithm: (a), (b), (c) Test distribution; (d), (e), (f) reconstruction of distribution using traditional LTD algorithm; (g), (h), (i) reconstruction of distribution using LTD-TV algorithm with 20000 iterations.

图2(d)—图2(f)使用传统LTD法重建, 接近度分别为0.6363, 0.5930, 0.5778; 图2(g)—图2(i)是在循环20000次的情况下用LTD-TV法重建, 接近度分别为0.1127, 0.1052, 0.1995, 分别比传统的LTD算法接近度减小82.28%, 82.25%, 65.47%. 从图2中可以看出, LTD-TV法重建图形相比传统的LTD法有很大的改善, 大大减小了重建的伪影. 该算法不但能重建出烟羽扩散的大致情况, 还能重建出丰富的细节. 然而, LTD法假设气体浓度是位置的二阶多项式, 重建过程其实是用二阶多项式去拟合高斯函数, 导致其峰值浓度比真实浓度低10%—15%, 扩散范围比测试图形稍大, 总体上有变圆的趋势.
2 -->

4.2.测量误差对重建结果的影响

外场实验中, 温度、湿度、气溶胶等都会影响DOAS的精确性. SO₂的测量值往往具有2%—20%的测量误差, 不好的重建算法可能会使误差在整个图形扩散, 造成重建结果偏离真实值, 甚至完全淹没在重建噪声中. 因此, 重建算法的抗误差能力是考察重建算法的一个重要方面. 由于误差来源的复杂性和不可重复性^[22,23], 一般采用叠加随机数的方法模拟误差对测量的影响^[13]. 给图2中测试图形的路径积分浓度${{S}_k}$加入加性噪声$\Delta {{S}_k}$

${{S}_k} = {{S}_k} + \Delta {{S}_k},\Delta {{S}_k}{\rm{ = }}f \cdot {{S}_k} \cdot {{R}_{{\rm{rand}}}},$

式中f表示误差系数, ${{R}_{{\rm{rand}}}}$是方差为1的随机数. 图3给出了f从1%—20%变化的情况下, 用LTD-TV重建图2中的测试图像得到的${\sigma _{{\rm{nearness}}}}$曲线.
图 3 接近度随误差系数变化曲线
Figure3. Nearness as the functions of error factors.

从图3可以看出, 随着误差系数f变大, LTD-TV算法的重建接近度也逐渐变大. 当误差系数达到某个值时, 接近度的上升变慢, 说明误差越大, LTD-TV算法的抗误差效果越明显. 从图3还可看出高斯数为2时重建误差最小, 对比图2(a)—图2(c)可以看出, 高斯数为2的测试图形的非零元素最少, 图像本身的稀疏性最强, 因此LTD-TV算法的性能也发挥得更加充分. 而在误差系数f变化的情况下, 传统LTD算法重建接近度仍然为0.6363, 0.5930, 0.5778, 比使用LTD-TV法重建的接近度更大. 总的来说, LTD-TV算法比传统的LTD算法的重建质量好得多.
2 -->

4.3.大气流动对重建结果的影响

已有的研究表明, 当烟羽位于重建图像的中心, 且两台DOAS与烟羽中心的连线之间成90°夹角时, 重建效果最好^[11,13]. 由于风向的不确定性, 该条件往往得不到满足. 图4中给出了图2(a)—图2(c)中高斯烟羽模型偏离初始位置时, 接近度随偏离距离变化的曲线, 偏移距离小于0时, 表示烟羽位置左移, 偏移距离大于0时, 表示烟羽位置右移.
图 4 接近度偏离距离变化曲线
Figure4. Nearness as the functions of offset distance.

从图4中可以看出, 当烟羽偏离图像中心时, 重建接近度和气体分布的稀疏性仍然相关. 当单高斯烟羽和双高斯烟羽偏离初始位置时, 重建图像的接近度发生振荡, 这种振荡是重建图像离散化造成的, 与以往的算法相类似^[11]. 除此之外, 接近度曲线未表现出明显的变化规律, 说明烟羽稀疏度较高的情况下, LTD-TV算法不受烟羽与IDOAS相对位置的影响. 而当三高斯烟羽模型偏离图像中心时, 重建接近度发生变化. 这是因为烟羽重建是一种极端的不完全角度重建, 当烟羽中气体浓度分布较为复杂时, 不同的烟羽位置会对IDOAS采集数据造成较大的影响, 从而造成重建效果发生起伏. 此外, 外场实验时, 特别是在湍流的影响下, 要表示烟羽中气体浓度的复杂分布, 需要更多的像素数目. 从信息论的角度来看, 在两台IDOAS信息采集能力有限的情况下, 低三阶导数模型增加的先验知识会将IDOAS采集的数据淹没, 导致重建结果偏离真实分布. 因此, 如果需要进一步提高LTD-TV算法的重建图像的分辨率, 只能进一步增加仪器的数量^[5,24].
基于烟羽位置与形状的复杂性, 使用LTD-TV法时仍然建议将IDOAS关于烟囱对称设置, 且两台DOAS与烟羽中心的连线之间成90°夹角.

外场实验数据在淮南某电厂外获得, 电厂烟囱高210 m, 在下风约180 m处搭建由两台IDOAS组成的断层扫描系统, 两台仪器距离约504 m. 实验当天天气晴朗且能见度高, 风向为东北偏东风, 风速为2 m/s. 通过307.5—318 nm波段反演SO₂路径积分浓度, 参与反演的气体包括SO₂(293 K, vanDaele), NO₂(294 K, vanDaele), O₃(243 K, Voigt), O₃(218 K, Brion)和ring光谱, 选取烟囱上风向的第一次测量谱作为参考谱. 图5以其中某一条光谱为反演实例, 展示了SO₂光谱数据的反演情况.
图 5 SO₂柱浓度IDOAS拟合反演实例　(a) SO₂柱浓度; (b)拟合残差
Figure5. An example of SO₂ SCD from IDOAS retrieval: (a) SO₂ SCD; (b) fitting residual.

图5(a)中黑线(不规则曲线)为测量光谱, 红线(光滑曲线)为拟合光谱, SO₂柱浓度为1.11 × 10¹⁷ molecules·cm^–2. 图5(b)所示的拟合残差小于0.00322. 选取其中一组SO₂数据作为研究对象, 由于气体浓度不可能为负, 重建图像时, 小于零的路径积分浓度被设置为0^[5]. 图6给出了使用LTD-TV法重建的污染气体断层图像.
图 6 烟囱烟羽SO₂分布重建图
Figure6. Reconstruction of SO₂ of stack plume.

在图6中, IDOAS#1位于X轴0点, IDOAS#2位于X轴504 m处. 烟囱中排出的SO₂在空中分成两股, 符合现场目测. 能够非常清楚地看到烟羽从两个中心扩散的情况, 甚至可以观察到风向在略微向IDOAS#2的方向变化, SO₂截面向该方向移动. 重建图中仅有右下角存在少量伪影, 其他位置有零星伪影, 不影响对重建图像的观察. IDOAS测得的柱浓度与重建图像的柱浓度值之间的关系如图7所示, 直观上看, 柱浓度的测量值与重建值符合得很好.
图 7 测量路径积分浓度与重建路径积分浓度对比　(a) IDOAS #1; (b) IDOAS #2
Figure7. Comparison between measured path integrated concentration and reconstructed path integrated concentration: (a) IDOAS #1; (b) IDOAS #2.

当气体分布未知时, 一般使用一致性相关因子衡量重建效果^[21].

${C_{{\rm{CCF}}}} = \rho A, $

式中$\rho $表示皮尔逊相关系数, A是路径积分浓度曲线发生位移时的校正因子,

$A = {\left[ {\frac{1}{2}\left( {\frac{{{\sigma _{{\rm{\hat S}}}}}}{{{\sigma _{\rm{S}}}}} + \frac{{{\sigma _{\rm{S}}}}}{{{\sigma _{{\rm{\hat S}}}}}} + {{\left( {\frac{{\hat S - S}}{{\sqrt {{\sigma _{{\rm{\hat S}}}}{\sigma _{\rm{S}}}} }}} \right)}^2}} \right)} \right]^{ - 1}}, $

其中S是IDOAS测得的路径积分浓度, ${\sigma _{\rm{S}}}$是测得的路径积分浓度的标准差; $\hat S$是重建图像的路径积分浓度, ${\sigma _{{\rm{\hat S}}}}$是重建图像路径积分浓度的标准差, 如果${C_{{\rm{CCF}}}} > 0.80$, 即可认为重建图像的路径积分浓度和实验测量值拟合得很好, 而图7所示重建结果的一致性相关因子为0.9063. 重建得到的路径积分浓度比测量值稍低, 符合数值模拟中出现的情况.

本文首次将CS理论引入气体断层重建领域, 提出了一种新的气体分布重建方法—LTD-TV算法. 数值模拟表明, 该方法将接近度由传统LTD法的0.5—0.6降低到0.1—0.2, 最大降低幅度达82.28%. 针对加性噪声, LTD-TV法表现出很强的抗噪声能力, 在路径积分浓度中引入20%的加性噪声的情况下, LTD-TV法仍然能重建出气体扩散的主要特征. 在外场实验中使用IDOAS采集数据, 并重建了的烟羽分布图形. 该方法不但假设气体的三阶导数是稀疏的, 而且隐含的假设气体浓度值非负, 相当于给投影方程组添加了更多的约束. 因此与以往的烟羽重建方法相比, 提高了气体重建的时间分辨率和可信度, 极大减少了伪影. 该方法不但能用于竖直方向的烟羽分布重建, 还能用于水平方向的区域气体分布重建. 不足之处在于: 使用二阶多项式作为气体扩散模型, 当气体近似按高斯函数分布时, 重建结果的峰值偏低, 扩散范围偏大.