0 引言
CT成像技术原理是使用多台相机在不同角度下同步拍摄火焰形成投影图像组,基于计算机断层层析算法进行反演实现三维燃烧场参数的重建[9,10,11]。CT成像技术可以与基于不同原理的成熟光学测量方法组合使用,如CT与粒子图像测速仪（particle image velocimetry, PIV）[12]结合,可以测量三维速度场;与火焰化学发光[13,14]结合,可测量三维温度场;与体积激光诱导荧光（volumetric laser-induced fluorescence, VLIF）[15,16]结合,可获得自由基团OH、CH等的空间分布;与体积激光诱导白炽光（volumetric laser-induced incandescence, VLII）[17]结合,可获得碳烟浓度、粒径三维分布。通常来说,在CT成像技术应用中,图像拍摄角度越多,三维重建结果越精确,但所配置的相机与计算机硬件成本越高。而角度下限受到火焰与拍摄条件影响和算法效率的制约[8,9]。因此,优选拍摄机位与优化重建算法成为CT成像燃烧诊断技术的关键。
根据实验条件的差异,在燃烧场CT成像诊断的策略方面,大致可分为单相机旋转-对象固定[9]与多相机-对象固定[10,14]两种方式。单相机旋转方法主要针对稳态火焰的测量,其优势在于高性价比,以及能提供无限多个冗余拍摄角度,可用于算法验证与优化。而多相机方法可以测量非稳态燃烧场,但受限于相机空间布置与高品质相机的价格,通常可获得的不同角度投影图像有限,对层析重建算法要求更高。在层析重建算法方面,研究者们建立了最大似然估计法（maximum likelihood estimation method, MLEM）[18]、反投影算法（back-projection, BP）[19,20]、迭代重建算法（algebraic reconstruction technique, ART）[9,21-22]、体素扩散模型（voxel spread function, VSF）[23]等方法,并进行针对性的优化[24,25,26,27],使所需的拍摄角度数从40个以上降低至8 ~ 17个。通过对不同算法模型的比较,研究者指出其中的VSF模型精度最高,但计算成本也最高,当燃烧场的规模不大时,截线长模型搭配ART算法具有良好的性价比[23]。此外,通过反射镜、中继棱镜、光纤传像束等方式[28,29,30,31],可将所需的相机数量降低至1 ~ 3台,能极大降低硬件成本。
在基于双色法碳烟火焰温度场的二维测量方面,国内外学者开展了大量的研究,双波长选择范围从可见光到红外光都被验证[32,33,34,35,36,37],其中的550 nm和650 nm的组合被广泛采用。而在碳烟火焰三维温度场测量方面,相关文献报道较少。基于此,本研究提出一种火焰旋转-单相机固定的多角度图像采集实验系统,与单相机旋转方法相比,该火焰旋转方法具有所需空间小、定位精度高的优势。采用单一增强型电荷耦合装置（intensified charge-coupled device, ICCD）相机对非旋转轴对称的稳定乙烯扩散碳烟火焰分别从28个不同角度依次进行自发光全波段和两组窄波段滤波后的拍摄,能够在较高的分辨率下采集数据。并通过建立与之相应的反演重建算法,再结合碳烟火焰修正双色法[38],对非规则稳定碳烟火焰三维温度场进行测量,以期实现更高性价比层析成像技术在碳烟火焰燃烧诊断中的应用。
1 研究方法
1.1 实验系统与火焰构造
实验系统由火焰旋转装置和相机固定的多角度图像采集系统组成,实物装置如图1所示。具体包括防震台、旋转台、同轴射流燃烧器及支架、高性能ICCD相机及两组窄带滤波片。防震台配有两组升降台和水平移轴,旋转台侧面标有角刻度,可以精确控制燃烧器的旋转角度和相机的拍摄位置。
Fig. 1 Photograph of the experimental system图1 实验系统实物图 |
Fig. 2 Co-flow jet burner with rotary table (a) and the C2H4 non-premixed soot flame (b)图2 可旋转同轴射流燃烧器（a）及其形成的乙烯扩散碳烟火焰（b） |
为保证火焰旋转后仍然维持稳定状态,燃烧器周围布置有防风围挡。火焰的旋转会导致周边流场改变,从而干扰燃烧场,因此需要以较慢的速度旋转燃烧器。在实验中,火焰旋转是通过用手拨动旋转台实现,旋转的角速度约为2 (°)/s,即火焰外缘的线速度约为0.03 m/s,燃烧器外缘的线速度约为0.15 m/s,对周围空气扰动较小。此外,每次旋转后静置时间约为20 ~ 30 s,然后以10 Hz连续拍摄100张照片,可以认为多角度下拍摄的火焰对象为同一稳定燃烧场。
1.2 相机标定与多角度投影图像获取
拍摄所使用的相机为单台ICCD相机（德国,LaVision,Pro X 4M）,分辨率为1 024 × 1 024,搭配窄带滤波片[（550 ± 10） nm与（650 ± 10）nm]。由于常温常压下乙烯扩散碳烟火焰温度低于2 500 K,火焰中碳烟热辐射的信号强度在可见光波段随波长增加而增加,因此,采用碳烟火焰双色法测温常用的550 nm和650 nm波段组合可以保证比色计算的可靠性[32,38]。在拍摄前,确定了最合适的快门速度（5 000 ns）和信号增益（80%）的组合。同时,在相应相机参数设置下,采用标准卤钨灯光源（美国,Thorlabs,SLS201L/M）对不同波段下光学探测系统的量子效率进行标定。
实验中,在 [0, 90) 角度区间内,每隔5° 进行一次拍摄,在 [90, 180] 角度区间内,每隔10° 进行一次拍摄,共获得28 × 3组结果,单一工况结果如图3所示。拍摄100张图像数据,先用LaVision软件进行时均化的预处理,然后导出用MATLAB软件进行后处理。在拍摄中,不同波段、不同角度均保持相机参数设置不变。
Fig. 3 Photos of the soot flame from different angles in full band, (550 ± 10) nm, and (650 ± 10) nm narrowband, respectively图3 全波段与（550 ± 10）nm、（650 ± 10）nm波段下不同角度拍摄的碳烟火焰图像 |
1.3 三维重建方法
建立一套三维反演重构的算法组,包括图像数据预处理（模块 ①）、光学坐标系定位（模块 ②）、重构计算算法（模块 ③）与重构结果展示（模块 ④）四个模块,其流程图如图4所示。首先是图像数据预处理模块,对输入图像进行自动背景去除与平滑滤波。与其并行的是光学坐标系定位模块,在输入所有拍摄相关参数后,通过射线与体网格的快速求交线方法,实现待求火焰空间与相机投影几何关系的系数匹配计算,从而实现投影矩阵、测量矢量与图像矢量的对应。其中,射线与网格的快速求交线部分采用了自主编写并优化的蛇形寻路算法,与枚举法和简单判别法相比,可以分别节约90%和70%的计算耗时[22]。算法组的核心部分为重构计算算法模块,包含迭代重建的ART算法、改进双色测温计算方法等,求解计算出火焰空间中每一点的信号值和温度值。最后是重构结果展示模块,对所得的燃烧场三维模型进行滤波,再采用自主编写的三维半透视散点剖窗成像算法等方法,进行全方位可视化的展示。
Fig. 4 Flow chart of 3D inversion and reconstruction algorithm图4 三维反演重构算法流程图 |
该反演重构算法组的核心为ART算法,其原理为：将待测火焰及周边区域,即感兴趣区域（region of interest, ROI）进行有限元化,分解为n × n × h的三维体素（voxel）阵列,作为图像数组。而拍摄所得的二维图片的每个像素（pixel）为测量数组,是由相机镜头在该角度下的每一条采样射线穿过ROI阵列后的层累结果,该结果由每一体素的信号值和射线在体素内的截线长共同决定,为一个n × n × h × m × h的五元投影数组,其中m为有效采样射线的数量。通过对光路的平行近似,可以将图像数组和测量数组的Z方向展开,成为图像矩阵和测量矢量,并且此时的图像矩阵为n阶方阵,可以进一步展开为1 × n2的图像矢量,投影数组也随之变为投影矩阵,是一个m × n2的满秩矩阵。此时的投影矩阵加权因子计算原理示意图如图5所示。
Fig. 5 Schematic diagram of approximate calculation of weighted factors for the projection matrix图5 投影矩阵加权因子近似计算原理图 |
图中由射线与矩阵求交算法（即蛇形寻路算法）数值计算得到LFG长度值,所得到的每个体素长度值可以组成投影矩阵A。此时存在线性关系,测量矢量p = 投影矩阵A·图像矢量x + 误差e。由于层析成像技术通常面临欠定问题,该线性关系若用最小二乘法求解,会因为秩亏（rank deficiency, RD）而无解。因此需要引入迭代计算的方法,即ART算法。其公式如式（1）所示：
$\overrightarrow{\boldsymbol{x}}^{(k+1, i)}=\overrightarrow{\boldsymbol{x}}^{(k, i)}-\boldsymbol{\varepsilon}_{\mathrm{ART}} \frac{\boldsymbol{A}_{i} \vec{x}^{(k, i)}-\overrightarrow{\boldsymbol{p}}_{i}}{\left\|\boldsymbol{A}_{i}\right\|_{2}^{2}} A_{i}^{\mathrm{T}}$ (1)
$T=\frac{C_{2}\left(\frac{1}{\lambda_{n+1}}-\frac{1}{\lambda_{n}}\right)}{\ln \frac{E_{n}}{\eta E_{n+1}}+(5+\alpha) \ln \frac{\lambda_{n}}{\lambda_{n+1}}}$ (2)
2 结果与讨论
2.1 三维火焰轮廓
构造的样例火焰的尺寸约为20 mm × 15 mm × 25 mm,因此反演重建选择的ROI尺寸设为30 mm × 30 mm × 30 mm。ROI的体素分辨率n × n × h在两次试计算阶段,分别设为10 × 10 × 10和30 × 30 × 30。后经过对算法的迭代优化,最终将体素分辨率选定为50 × 50 × 50,即每一体素为棱长0.6 mm的正立方体,在运算条件许可的情况下保证了较高的结果精度。输出二维图像时,采用3级插值加密,图像分辨率为393 × 393;输出三维图像时,火焰轮廓采用1级插值加密,分辨率为99 × 99 × 99,温度场采用2级插值加密,分辨率为197 × 197 × 197。
采用28张全波段图像反演重建出的火焰三维轮廓形状如图6所示,已对亮度值进行归一化处理,分不同显示阈值和视角进行展示。当取亮度显示阈值较低时,两个火焰趋近于连为一体;取亮度显示阈值较高时,两个火焰中间较暗的部分不再相连,火焰轮廓与实际火焰外形相符,证明本文所采用的重构方法可行。
Fig. 6 3D scatter images of flame profile under different thresholds and perspectives: (a) lower brightness threshold; (b) higher brightness threshold图6 不同显示阈值与视角下火焰轮廓的三维散点图：（a）较低亮度显示阈值;（b）较高亮度显示阈值 |
2.2 三维温度场
Fig. 7 3D semi-perspective sectional images of temperature field distributions: (a) temperature distributions in different cross-sections; (b) temperature distributions in different longitudinal-sections图7 温度场的三维半透视截面：（a）不同横截面温度分布;（b）不同纵截面温度分布 |
Fig. 8 2D sectional images of temperature field distributions: (a) temperature distributions in different cross-sections; (b) temperature distributions in different longitudinal-sections图8 温度场的二维截面：（a）不同横截面温度分布;（b）不同纵截面温度分布 |
Fig. 9 3D semi-perspective scattered windowed images of temperature field under different window situation and perspectives:(a) windowed in the second octant; (b) windowed in the forth octant 图9 不同剖窗位置与视角下的三维温度场半透视散点剖窗成像图：（a）剖窗在第二卦限;（b）剖窗在第四卦限 |
此外,将反演重建所得的三维光强信号场,向某一个角度累加所获得的二维层累图像,与该角度下直接拍摄的层累图像进行对比,同样可以用于评估重建结果的准确性。图10是投影角度为90° 时的对比图,在采用红蓝影叠加的基础上,分别描绘了以25%信号为阈值的轮廓线。可见重建结果与拍摄结果在核心部分吻合度很高,边缘部分存在一定的不吻合现象。对重合程度进行量化计算,重合面积占比为88.9%,加权重合面积占比为97.6%。火焰最顶部和尾部的尖锐部分,重建结果存在一定的缺失。这一现象可能是由于其本身信号值偏弱、拍摄过程中存在细微抖动和体素分辨率还不够高等多个因素共同导致。
Fig. 10 Comparison between the photo shot directly and the image summed by 3D reconstruction field from 90°图10 投影角度为90° 时,直接拍摄图像与重建信号场层累后图像的对比 |
2.3 冗余验证
在实验过程中,以下三个原因可能会导致重构结果产生误差。其一是火焰不完全稳定,不同时序所拍摄的不同角度图像之间有细微差异而造成的误差,可以通过控制燃料量和火焰高度以提高稳定性、每次旋转燃烧器后须静置待其稳定再拍摄、增加拍摄时的均值化张数等方法来减小该误差。其二是燃烧器调平欠佳时会随着旋转而偏移所造成的误差,可以通过分层精确调平来规避误差,或者在每一拍摄角度都分别设置标尺从而进行后期校准来规避误差。其三是算法本身在迭代残差取舍时产生的误差,以及对前两项误差可能产生的局部积累,可以用模型滤波算法来减小该误差。
在ROI体素分辨率分别取10 × 10 × 10和30 × 30 × 30的两组预实验中,离散病态性比值分别为1.1%和2.3%,在体素分辨率取50 × 50 × 50的正式实验中,离散病态性比值有所提高,为4.6%,依然处于较低的水平。考虑到火焰对象集中分布于ROI中30 × 25 × 40的中央区域,外围的噪点均可以舍去,则中央离散病态性比值为2.0%,证明该实验结果的可靠性较高。
在实验设计方案下,选取的拍摄角度数目存在冗余,可用于进一步研究离散病态性与拍摄角度数目的关系。将拍摄角度数目由28顺序递减至8的过程中,中央区域离散病态点占比的变化如图11所示。此外,在拍摄角度为27和26个时,通过比较减少不同角度的图像后离散病态性比值的变化,发现有27个角度的比值变化在0.2%以内,这一结果验证了不同角度所拍摄火焰图像的一致性。
Fig. 11 Relationship of discrete ill-posedness voxel proportion in central zone varied with number of shoot angles图11 中心区域DIP体素占比随拍摄角度数目变化关系图 |
3 结论
搭建了一套火焰三维计算层析成像测试系统,并建立了与之匹配的反演重构算法及程序。在同轴射流燃烧器上用导流片构造出非规则的乙烯-空气扩散火焰,采用单一ICCD相机结合旋转台从28个不同角度依次进行拍摄,先根据自发光的全波段拍摄信号,反演重构得到碳烟火焰的三维轮廓,再拍摄550 nm和650 nm窄波段的自发光信号,结合改进双色法计算得到碳烟火焰的三维温度场,并验证了重建结果的准确性。
本方法当前可用于任意非规则形状的稳态火焰自发光三维燃烧场参量测量,具有高效率、高性价比的优势。进一步结合多角度图像同步采集,所建立的反演重构方法后续将在湍流火焰三维测试中应用与验证,考察其对非稳态火焰测量的适用性。
