基于图像重建方法获得的气体浓度二维分布图是研究室内空气环境的一个重要手段，气体浓度分布图揭示了环境中气体扩散的特征并用于定位气体扩散源^[1-3]。因此，众多学者提出了改进的图像重建算法^[4]，代数重建技术(algebraic reconstruction technique, ART)是一种在不完备投影数据重建图像时具有明显优势的图像重建算法^[5]，但ART算法的计算量大、重建速度慢等限制了它的实际应用。为加快重建速度和提高重建精度，学者们主要从迭代公式^[6-8]、松弛因子^[9-11]、投影系数、加入合理的正则化约束^[12]等方面改进ART算法。

投影系数计算的时间在整个重建过程中占90%以上^[13]，快速有效地计算投影系数是提高重建速度的关键，因此不断有学者提出新的投影系数快速计算方法。计算投影系数即计算光线在重建区域离散化后各网格内的光程。Siddon^[14]于1985年提出有效的投影系数计算方法，但由于该方法包含大量求交、排序操作，导致整体的计算时间成本大。针对遥测起点固定在重建区域外或边界上的情况，以下学者提出了快速的投影系数计算方法：张顺利等^[15]及Xue等^[16]在经典的Siddon算法基础上做了改进，通过增量计算确定光线穿过的网格编号及投影系数，该方法比Siddon算法快近6倍，但引入了大量的分支判断，不适合在GPU端使用；Li等^[17]根据线性关系使用2个数组求投影系数，该方法在GPU端也有较好的加速效果，但在投影系数计算方面仍存在较多的计算量；程志勇等^[18]在此基础上进一步减少了计算量；杨文良和魏东波^[19]通过简单的加减法及比较运算完成投影系数的快速递推计算，进一步提高了投影系数的计算速度；Zhang等^[20]提出一种用于基于线积分模型的迭代重建算法的快速计算正投影算法，并通过消除排序操作改进反投影计算方法，提高了迭代重建算法的重建速度；Wang等^[21]提出一种锥形光束的快速投影系数计算方法。

随着气体浓度二维分布探测技术的发展，Li等^[22]和王阳等^[3]利用位置固定的平行光束重建水平和竖直截面的有害气体浓度二维分布图，该方法响应速度快且检测精度高，可用于气体浓度的实时检测，但需要同时使用多台气体遥测收发装置。Li等^[1]利用固定点发出的扇形光束重建不同水平和竖直截面的有害气体浓度二维分布图，该方法仅利用1台气体遥测仪，但重建精度较低。本文基于移动激光气体遥测技术重建室内气体浓度二维分布图，仅需1台气体检测仪可以获得大量平行和扇形光束，重建精度高。但在实际中遥测路线不是严格沿着网格线的，遥测起点会随机分布在重建区域内的网格内或网格线上，故不能直接使用现有的快速投影系数计算方法。

基于上述问题，本文选择ART算法和改进的自适应代数重建算法(modified adaptive ART, MAART)^[10]作为图像重建算法，在现有的快速投影系数递推计算方法^[19]的基础上，提出一种无起点约束的快速投影系数递推计算方法。通过设计3个光程累积因子，在不降低原算法计算速度的情况下，统一了遥测起点在网格内和网格线上的递推过程，实现了基于移动激光气体遥测技术的快速投影系数计算。

1 ART算法基本原理

利用ART重建气体浓度二维分布图，首先需要离散化图像，即将整个气体浓度二维分布图像f(x, y)离散化为N个网格，网格内的气体浓度由吸收系数反演得到^[7]，其网格值f_i(1≤i≤N)表示第i个网格内的气体吸收系数α_i。假设一共有M条激光，第j条经过重建区域的光线投影值即为光线路径的积分吸光度，记为A_j(1≤j≤M)，由激光气体遥测技术测量得到。因此重建图像和投影值之间的线性方程组为

$\left\{\begin{aligned} A_1 & =\sum\nolimits_{i=1}^N \omega_{i 1} \alpha_i, \\ A_2 & =\sum\nolimits_{i=1}^N \omega_{i 2} \alpha_i, \\ & \vdots \\ A_M & =\sum\nolimits_{i=1}^N \omega_{i M} \alpha_i . \end{aligned}\right.$

(1)

其中：ω_ij为投影系数，它在数值上等于第j条光线在第i个网格内的长度。因此图像重建问题可转化为求解高维线性方程组，但实际投影数据中包含各种噪声，这导致线性方程组通常是无解的或有无穷解^[21]。因此线性方程组不能直接用解析法求解，需要用迭代算法逼近各网格的真实吸收系数。一般设置吸收向量α=(α₁, α₂, …, α_N)的初值向量为零向量，代入式(1)进行迭代计算，将α^(j-1)代入式(1)的第j个方程时，得到α^(j)的迭代公式^[23-24]为

$\alpha_i^{(j)}=\alpha_i^{(j-1)}+\lambda \frac{A_j-\sum\nolimits_{n=1}^N \omega_{n j} \alpha_n^{(j-1)}}{\sum\nolimits_{n=1}^N \omega_{n j}^2} \omega_{i j}.$

(2)

其中：λ(0 < λ < 2)是松弛因子，与收敛速度和重建质量有关；n(1≤n≤N)表示网格序号。将ART中的λ改成一个自动调整的松弛参数，称为MAART^[10]，式(2)的λ表示为

$\lambda_j^{(k)}=\beta \times\left(\omega_{i, j} \alpha_j^{(k)}\right) /\left(\sum\limits_{j=1}^J \omega_{i, j} \alpha_j^{(k)}\right) .$

(3)

其中：k代表MAART迭代过程中的迭代次数；β是一个常数，本文取值0.25^[10]。依照上述迭代，由式(1)第M个方程得到α^(M)，则认为完成了一次完整的迭代。每次迭代中判断吸收向量的变化是否大于阈值且迭代次数是否小于10 000。若式(2)收敛且收敛速度较快，则在10 000次迭代前因吸收向量的变化小于阈值而提前结束迭代；若式(2)收敛但收敛速度较慢或者不收敛，则会因超过10 000次迭代而结束迭代。

ART算法的核心过程是正投影和反投影，正投影是网格沿光线方向的吸收系数的总和，反投影是将测量投影值与正投影值之间的差异根据松弛因子更新网格的吸收系数^[25]。两者都反复用到投影系数，因此快速准确地计算投影系数ω是决定该算法能否快速重建的关键因素。

投影系数有多种数学计算模型，综合考虑计算速度和精度后本文选择长度加权模型。在这种模型下，多用递推方法提高投影系数的计算速度，本文将优化其递推过程以实现无起点约束的投影系数快速计算。

2 改进的投影系数计算方法 2.1 算法原理

为统一遥测起点在网格内和网格线上的递推过程，本文改进了现有的快速投影系数递推计算方法来获得投影系数。在改进前方法中，用当前网格的相邻网格线间的光程长度递推计算当前网格的投影系数和网格序号，当遥测起点在网格内，其在相邻网格线间的光程长度无法由改进前方法中的公式来计算。本文通过从遥测起点至当前网格的3个光程累积因子间大小关系判定光线下一个相交的网格序号，同时利用3个光程累积因子递推计算出当前网格的投影系数，因光程累积因子本身包含遥测起点位置信息，所以遥测起点在网格内或网格线上对递推过程无影响。

如图 1所示，根据ART算法需要将二维图像离散化为N=m×n个网格，设网格边长δ=1，网格序号按从左到右、从下到上的顺序编号，光线1的遥测起点坐标为(x_start, y_start)、倾斜角为θ∈(-π, π]。则光线贯穿相邻2个垂直网格线的光程为$v_0=\sqrt{k^2+1}$，光线贯穿相邻2个水平网格线的光程为$h_0=\sqrt{1+\frac{1}{k^2}}$，其中k=tanθ，如图 2所示。起点所在网格为初始网格，递推方向为从初始网格开始依次遍历光线经过的网格，直至到达重建区域的边界。

下面定义3个重要的光程因子：一是垂直光程累积因子v，表示从起点到当前网格即将穿出的竖直网格线之间的光程；二是水平光程累积因子h，表示从起点到当前网格即将穿出的水平网格线之间的光程；三是已遍历的光程累积因子d，表示从起点到当前网格之前的光程。如图 2所示，当光线1从起点开始递推到序号为q的网格时，垂直光程累积因子为v=v_start+3v₀，水平光程累积因子为h=h_start+h₀，已遍历的光程累积因子d=v_start+2v₀。其中v_start为初始网格的垂直光程累积因子，h_start为初始网格的水平光程累积因子，其值由起点坐标和倾斜角计算得到。

根据当前网格q的v和h的大小关系，可划分成光线从水平边界、垂直边界、顶点穿出的3种情况，全文以光线倾斜角$\theta \in\left(0, \frac{\mathsf{π}}{2}\right)$为例，同理可推广到$\theta \in\left\lceil\frac{\mathsf{π}}{2}, \mathsf{π}\right\rceil \cup(-\mathsf{π}, 0]$的情况。如图 3所示，分别用光线2~光线4、光线5~光线7、光线8~光线10表示对应的3种情况。接下来分情况讨论递推计算投影系数和网格序号的方法。

1) 若网格q的h < v，网格q的投影系数s=h-d，下一个穿过网格的序号为q+n。如图 3(a)所示, 其中光线2在网格q处的3个光程累积因子分别为v_{2, q}、h_{2, q}、d_{2, q}，此时h_{2, q} < v_{2, q}；网格q的投影系数s_{2, q}=h_{2, q}-d_{2, q}，下一个穿过网格的序号为q+n，网格q+n的v_{2, q+n}=v_{2, q}、h_{2, q+n}=h_{2, q}+h₀、d_{2, q+n}=h_{2, q}。同理可以计算光线3、4在网格q处的投影系数和网格序号。

2) 若网格q的h>v，网格q的投影系数s=v-d，下一个穿过网格的序号为q+1。如图 3(b)所示, 其中光线5在网格q处的3个光程累积因子分别为v_{5, q}、h_{5, q}、d_{5, q}，此时h_{5, q}>v_{5, q}；网格q的投影系数s_{5, q}=v_{5, q}-d_{5, q}，下一个穿过网格的序号为q+1，网格q+1的v_{5, q+1}=v_{5, q}+v₀、h_{5, q+1}=h_{5, q}、d_{5, q+1}=v_{5, q}。同理可以计算光线6、7在网格q处的投影系数和网格序号。

3) 若网格q的h=v，网格q的投影系数s=h-d=v-d，下一个穿过网格的序号为q+n+1。如图 3(c)所示, 其中光线8在网格q处的3个光程累积因子分别v_{8, q}、h_{8, q}、d_{8, q}，此时h_{8, q}=v_{8, q}；网格q的投影系数s_{8, q}=h_{8, q}-d_{8, q}=v_{8, q}-d_{8, q}，下一个穿过网格的序号为q+n+1，网格q+n+1的h_{8, q+n+1}=h_{8, q}+h₀、v_{8, q+n+1}=v_{8, q}+v₀、d_{8, q+n+1}=h_{8, q}=v_{8, q}。同理可以计算光线9、10在网格q处的投影系数和网格序号。

通过上述讨论可知，光线在当前网格q的投影系数s及下一个网格的序号q_next的递推公式如下

$s= \begin{cases}h-d, & h <v, \\ v-d, & h>v, \\ h-d=v-d, & h=v\end{cases}$

(4)

$q_{\text {next }}= \begin{cases}q+n, & h <v, \\ q+1, & h>v, \\ q+(n+1), & h=v .\end{cases}$

(5)

其中：h为光线在网格q处的水平光程累积因子；v为光线在网格q处的垂直光程累积因子；d为光线在网格q处的已遍历的光程累积因子；n为重建图像离散化后的网格列数。为便于计算下一个遍历到的网格q_next时的投影系数，需要在每次递推过程中更新光程累积因子，公式如下

$\left\{\begin{array}{llll} h_{\text {next }}=h+h_0, & v_{\text {next }}=v, & d_{\text {next }}=h, & h <v, \\ v_{\text {next }}=v+v_0, & h_{\text {next }}=h, & d_{\text {next }}=v, & h>v, \\ h_{\text {next }}=h+h_0, & v_{\text {next }}=v+v_0, & d_{\text {next }}=h=v, & h=v . \end{array}\right.$

(6)

其中: h_next为光线在网格q_next处的水平光程累积因子；v_next为光线在网格q_next处的垂直光程累积因子；d_next为光线在网格q_next处的已遍历的距离；h₀为光线贯穿相邻2个水平网格线的长度；v₀为光线贯穿相邻2个垂直网格线的长度。

2.2 算法步骤

本文计算投影系数的算法流程如图 4所示。本算法利用3个光程累积因子统一了遥测起点在网格内和网格线上的投影系数和网格序号的快速递推计算，每次循环体中只有少量的比较和加减运算，省去了耗时严重的求交、排序运算，因而相比传统的Siddon算法大大减少了计算时间。

2.3 算法复杂度分析

假设遥测光线在重建区域内分别与一系列水平和垂直网格线相交，n为光线经过的网格数，则本文算法的递推次数为n。设一次比较、加法、减法、乘法、除法运算和floor()函数调用开销分别为: t₀、t₁、t₂、t₃、t₄和t₅, 由式(4)~式(6)可知，计算光线经过网格的投影系数和网格序号至少需进行n次比较、2n次加法和n次减法，最多需进行2n次比较、3n次加法和n次减法，所以本文计算一条光线经过网格的投影系数和网格序号的开销为[nt₀+2nt₁+nt₂, 2nt₀+3nt₁+nt₂]。现有快速投影系数递推计算方法计算一条光线经过网格的投影系数和网格序号的开销为[nt₀+nt₁+2nt₂, 2nt₀+3nt₁+nt₂]^[19]。Siddon算法计算一条光线的开销为nt₀+5nt₂+3n(t₁+t₃+t₄)+2nt₅^[15]。

由此可见，本文算法和现有的快速投影系数递推计算方法仅使用了比较、加法和减法运算，Siddon算法中多次使用乘法、除法和取整运算，所以本文算法和现有的快速投影系数递推计算方法的开销远少于Siddon算法。相比现有的快速投影系数递推计算方法, 本文算法多1次加法运算，少1次减法运算。计算机执行减法操作时，将负数用其补码表示，从而将减法运算转换成加法运算^[26]，因此相比现有的快速投影系数递推计算方法, 本文算法没有增加算法复杂度。

3 实验与结果

实验采用的计算机配置为Intel Core i7-10710U (1.10 GHz)，软件开发环境为C++开发程序。

3.1 基于移动激光气体遥测技术重建室内气体分布的模拟实验

由气体扩散模型^[27-28]可知，当发生微小泄漏时，可燃气体无法瞬间充满整个密闭房间，气体的扩散是一个长时间变化的过程，一般监测气体泄漏设备相邻2次测量时间间隔约为10 min，而基于移动式气体遥测仪检测室内气体浓度的时间取决于移动平台的行驶速度。通过优化遥测路线可在1 min左右完成全部光强数据的采集工作，采集期间密闭房间内的气体浓度分布近似不变。

本文采用gazebo仿真软件模拟移动激光气体遥测仪在20 m×20 m密闭空房间内的遥测路线，以模拟随机分布的遥测起点。设计遥测路线近似为正方形，以获得近似正交分布的光束。遥测路线如图 5(a)所示，遥测路线的起点随机选取靠近房间墙壁的一点，用图中左上角的倒三角形表示，方向为顺时针，其遥测时的光线方向与行驶方向垂直，行驶时的平均速度约为2.5 m/s，在拐角处停止行驶，仅以约0.628 rad/s的角速度旋转。实验获得的光线分布如图 5(b)所示，图 5(b)选择较大的遥测间隔(0.4 s)，以更清楚地显示光线分布，其中圆点表示遥测起点，遥测起点随机分布在网格内和网格线上，大部分遥测起点分布在网格内，少部分遥测起点分布在网格线上。

设某次光强数据采集工作内，单个气体扩散源在密闭空间下的气体浓度分布符合Gauss分布，气体质量浓度从1.0×10⁵ mg/m³变化到2.0×10⁵ mg/m³。建立1、2、3 Gauss峰模型分别模拟1、2、3个气体扩散源在密闭空间下的气体浓度分布，如图 6(a)、6(c)、6(e)所示。为获得较好的重建结果，应选择较小的遥测间隔，当遥测间隔为0.025 s时，共获得1 751条遥测光线，利用本文算法结合代数重建技术得到的重建结果如图 6(b)、6(d)、6(f)所示。本实验的初值向量为α⁽⁰⁾=(0.1, 0.1, …, 0.1)。

用相对均方根误差(relative root mean square error, RRMSE)量化数值模拟图和重建结果图之间的一致性。RRMSE^[1]表示相对于平均模拟值的重建浓度与模拟值之间的误差大小，RRMSE是非负的且越低越好。其公式为

$\text { RRMSE } \left.=\sqrt{\frac{\sum\nolimits_{i=1}^n\left(C_{\text {ree }, i}-C_{\text {sim }, i}\right)^2}{n}} \right/\overline{C_{\text {sim }}} .$

(6)

其中：C_{rec, i}和C_{sim, i}分别表示第i个网格的重建浓度值和模拟浓度值；n为网格总数；$\overline{C_{\text {sim }}}$为平均模拟浓度值。投影系数的快速计算不会影响重建精度，但如果投影系数计算错误会导致重建误差增大。在1、2、3个气体扩散源模型下，对比本文算法和Siddon算法在ART算法和MAART算法下的RRMSE，如表 1所示。由表 1可知，在3个模型下，MAART算法的RRMSE均小于ART算法的RRMSE，MAART算法能减小重建误差；本文算法和Siddon的重建误差一致，验证了本文算法的有效性。本文算法能实现无起点约束的投影系数计算，可应用到基于移动激光气体遥测技术中，为研究室内气体环境提供帮助。

进一步讨论遥测光线的数量与重建结果的关系，设采样间隔分别为0.4、0.3、0.2、0.1、0.05、0.04、0.025 s，共采集7组(110、146、219、438、876、1 095、1 751条)遥测光线，在不同遥测光线数量下的基于ART算法的RRMSE值和重建时间如图 7所示。由图可知，随着采样间隔减小，遥测光线数量增加，重建结果的RRMSE逐渐减小，但重建时间逐渐增加。当光线数量大于1 095条，RRMSE减小的速度变慢；当遥测光线数量大于876条，重建时间增加的速度变慢。由此可以看出采样间隔减小到一定程度后，光线数据存在冗余对重建图像的修正作用减弱。

3.2 投影系数计算速度的比较实验

对比本文算法、Siddon算法和改进的投影系数计算方法的平均投影系数计算时间，以验证本文算法能实现投影系数的快速计算。其中改进的投影系数计算方法的遥测路线需要沿着网格线。因此，以网格大小为20×20为例，设计遥测路线和遥测方式如图 8所示：遥测路线为从点(1, 1)—(1, 19)—(19, 19)—(19, 1)—(1, 1)的折线；遥测方式为从起点开始每隔1个网格为1个遥测起点，1个遥测起点上有1组光线，光线方向为从行驶方向到行驶方向顺时针旋转180°(不含180°)，相邻2个光线间的夹角为5°，拐点上因行驶方向改变有2组光线，共有4×18×36条遥测光线。

在上述遥测路线和遥测方式下，计算不同网格大小(20×20、40×40、80×80、128×128、256×256)下一次迭代中的平均投影系数计算时间。如表 2所示，其第2列为对应网格大小下的遥测光线数量，第3~7列分别为Siddon算法^[14]、文献[18]和文献[15]中改进的算法、现有的快速投影系数递推计算方法^[19]和本文方法的1次迭代中的平均投影系数计算时间。由表可知，划分的网格数量越多本文算法的加速效果越好；在256×256网格下，本文算法的平均投影系数计算时间为1.434 s，传统Siddon算法的平均投影系数计算时间为19.95 s，本文算法的计算速度约是传统Siddon算法的13.9倍；本文算法的耗时均少于现有改进的快速投影系数递推计算方法以及文献[18]和文献[15]中改进的算法，证明了本算法是快速的且相较现有的快速算法没有增加计算量。

4 结论

室内气体二维分布图重建方法中，与固定激光气体遥测技术相比，移动激光气体遥测技术更适用于多场景和复杂场景中。由于移动激光气体遥测技术的遥测起点会随机分布在重建区域的网格内和网格线上，给投影系数的递推计算过程带来较大的挑战。现有的快速投影系数递推计算方法未充分考虑遥测起点在网格内和网格线上的不一致性，导致其不适用于基于移动激光气体遥测技术的图像重建过程。

针对上述问题，本文提出一种二维图像重建中无起点约束的快速投影系数递推计算方法，通过设计3个光程累积因子：垂直光程累积因子、水平光程累积因子和已遍历光程累积因子，统一了遥测起点在网格内和网格线上的投影系数及网格序号的快速递推过程，每次递推过程仅涉及简单的比较和加减运算，相比改进前的算法未增加每次递推的计算量。

基于室内气体浓度二维分布模拟数据和ART算法开展了重建实验，实验结果验证了本方法的有效性。在重建精度方面，本文重建结果中最小RRMSE约1.16%。在计算速度方面，本文算法优于传统Siddon算法和现有的快速投影系数递推计算方法，最快约是经典Siddon算法的13.9倍，且网格数量越多速度优势越突出。本文提出的计算方法实现了无起点约束投影系数的快速计算，并为基于移动激光气体遥测技术的二维气体浓度分布快速重建提供了理论依据。