基于三维图像分析的浮力扩散火焰几何与辐射特性研究

余 辉，雷 佼，邓文扬，李元洲

(中国科学技术大学火灾科学国家重点实验室，合肥，230026)

0 引言

火蔓延预测预报是森林火灾防治工作的重要内容。准确计算火前锋的辐射热流分布可显著提升火蔓延物理模型的预测精度。目前，经典火焰辐射预测模型包括点源辐射模型和固体火焰模型[1-6]。Orloff和De Ris[7]在计算燃料表面的辐射热反馈和向火焰目标的辐射传热时，提出了轴对称火焰辐射优化模型。该模型利用二维图像数据估算火焰三维特征，即将二维火焰图像的每一层视为单独的圆柱体，获取了火焰的平均体积和表面积，从而简化计算火焰的近场辐射，并发现火焰体积热释放速率为定值。Hu等[8]将上述方法用于计算射流火焰的平均体积和体积热释放速率。但是，这些传统火焰形态和辐射模型计算结果针对的是形态与辐射的平均值，忽略了火焰的脉动和不稳定性。在真实火场燃烧行为中，火焰的非轴对称性(径向方向沿着轴线非对称)和强脉动性使得经典辐射模型适用性较差。

近年来，随着实验测量技术的不断发展(激光技术、断层技术)，研究者们开始研究火焰三维几何特性的变化规律。Bheemul等[9]利用多台同步相机重构了轴对称火焰的三维轮廓，获得了火焰体积、表面积、方向、长度和圆度。Mason和Rogers[10]利用两台相机拍摄沙发火和丙烷火，利用所获取的火焰三维概率分布云图计算了火焰周边辐射场。Doi和Sato[11]使用滤光片获得了特定波长下的火焰图像，并利用层析断层重建和比色方法获得了三维火焰结构和温度场分布。Grauer等[12]基于背景纹影技术(BOS)方法测量了旋流火焰区域的三维折射率场。Rossi等[13]最早使用双目视觉技术估算了森林火灾火前锋的表面积和体积，并利用已知立方体估计误差。虽然前人已经研究了火焰的三维结构，但多注重三维重构算法和火焰几何特性研究，较少关注火焰几何与辐射特性的关系。

本文使用7台CCD工业相机同步采集了小尺度圆形燃烧器上的二维瞬时丙烷火焰图像，利用Visual Hull方法重构了三维瞬时火焰图像，得到了火焰高度、表面积、体积随时间和热释放速率的变化规律。在此基础上，计算了三维瞬态火焰表面面元对外部目标的视角系数和辐射热流分布，提出了具有普遍适用性的三维火焰辐射预测模型，并验证了该模型的可靠性。

1 三维重构内容概括

1.1 三维重构原理

在测量火焰时，为确定空间物体表面某点的三维几何位置与其在图像中对应点之间的相互关系，必须建立相机成像的几何模型(图1)。

图1 相机成像模型Fig.1 Camera imaging model

图1包含4个坐标系。其中，Ow-XYZ、Oc-XYZ分别是世界坐标系和相机坐标系，o-xy和o-uv分别是图像坐标系和像素坐标系。P是世界坐标系中的任意点，在空间中代表火焰区域。p是点P的图像投影点坐标(x,y)。利用Zhang[14]的标定算法，通过对标定板进行网格成像，确定每一台相机的内部参数矩阵、畸变系数和旋转平移向量。利用相机内部参数的矩阵变换得到P在相机上的投影位置p。

三维重建过程将火焰视为唯一光源，图像中每一个亮点代表火焰在图像上的投影，该区域代表火焰在平行于相机光轴方向上的外形轮廓信息，对应于空间的一个椎体。利用交集方法可以获得火焰的外包体，其过程如图2所示。当有足够多的投影图像时，便可得到火焰的三维形状。通过Cubes[15]算法得到的等值面，可用于表征火焰的三维表面。该表面由多个三角形面片构成(图3)。火焰的高度、体积根据剩余体素计算得到，根据每个三角形面片的顶点P0，P1和P2的三维坐标，利用累加法计算火焰表面积和视角系数[10]。

图2 火焰图像和反投影所相交的轮廓Fig.2 Intersected contours of flame images and back projections

图3 Marching Cube重构的三角形表面Fig. 3 Reconstructed triangular surface by Marching Cube algorithm

(1)

其中r是三角形重心到目标点的距离，θS是三角形法线与三角形重心到目标的直线所指向的矢量之间的角度，θT为目标法线与目标到三角形重心的直线所指向的向量之间的角度，Si为三角形的面积。因此，对整个火焰表面进行积分，得到总体的视角系数：

(2)

1.2 三维重构过程

如图4所示，本研究中的火焰三维重构过程包括：

1)相机标定：根据不同方向相机的标定图片，利用软件进行标定，获取相机的内、外参数，进而通过投影算法获取火焰轮廓。本文使用Zhang[14]提出的单平面棋盘格标定算法，对标定图像进行处理。

图4 火焰重构步骤Fig.4 Steps of flame reconstruction

2)图像采集：获取不同方向CCD相机标定图片和瞬态火焰图像。

3)构建网格区域：根据实验工况对火焰存在的区域进行网格划分。

4)图像分割：将原始彩色图像转为灰度图像，通过灰度像素值与阈值的比较，再将灰度图像转换为二值图像。

5)反投影计算：根据网格点在各个方向的投影，判断该区域是否在火焰区域内。

6)统计目标点个数：根据对网格的遍历，统计在火焰区域内的坐标点个数，计算火焰体积。

7)等值面计算：建立三维离散坐标系，利用Marching Cube算法计算火焰等值面，并进行火焰三维表面重构。

8)数据后处理：通过三角形面元积分方法，计算火焰的表面积和角系数等。

由于火焰三维重构将用于计算火焰体积和表面积，图像处理过程需要较高的精度，因此必须精细处理火焰底部较暗的区域(无法识别的火焰区)。本研究使用差异化阈值对火焰进行分割，如图5所示，根据不同区域(图像中水平和竖直像素坐标区域)火焰亮度特点使用差异化阈值进行图像二值化。

图5 分区差异化域值的效果图Fig. 5 Differentiated threshold effect in different areas

1.3 误差分析

本研究结果可能包含火焰形态和视角系数的计算误差。由于无法获知火焰真实的三维轮廓，本研究通过重构已知规格的圆柱体，计算重构物体高度、表面积、体积误差。该圆柱体直径为7.5 cm，高度分别为30.0 cm、40.0 cm和50.0 cm。如表1所示，体积和表面积的预测误差均在5%以内，最佳精度为2.5%。Hankinson[16]模拟了地面上的直圆柱体(D=4 m,H=10 m)和面向圆柱中心的垂直目标，发现视角系数计算精度取决于圆柱表面划分元素的数量和目标与圆柱体的距离。当三角形数量在9 524～27 106范围时，视角系数计算误差会随着目标与圆柱距离的减小而增大。不过，即使距离减小到130 mm，计算误差仍保持在1%以内。本研究中火焰的三角形面片数量均大于9 524，辐射热流计与火源距离满足相似变换下精度为1%的要求。

表1 七个投影面重构圆柱的误差Table 1 Errors of reconstructed cylinders by seven projection planes

2 实验设置

图6 实验布置与实验平台图Fig. 6 Schematic diagram of the experimental setup

在距离火源中心1.7 m～2.0 m的周向上等间距(26°)布置了7台CCD工业相机(分辨率：1 920 ×1 200 pixels，图6(a))，用于记录7个方向的瞬时火焰形态。使用高精度同步器实现CCD相机和辐射热流采集模块的同步。利用水冷式辐射热流计(响应时间：25 ms)测量火焰外部竖直方向(热流计测量面竖直，编号q1～q4)和水平方向(热流计测量面水平，编号q5)的辐射热流。两组辐射热流计布置在两个正交的半径方向上。辐射热流计距离火源中心为0.3 m和0.6 m，分别用于表征火源近场辐射[17](L<5D，其中L是辐射热流计与火源中心的距离)和火源远场辐射(L>5D)。为避免热流计对火焰辐射的遮挡，L=0.6 m处的辐射热流计略高于(0.1 m)L=0.3 m处的辐射热流计。在燃烧器中心轴线0.05 m～0.80 m高度上均匀布置16根结点直径0.374 mm的K型热电偶，用于获取中心轴线火焰温度变化规律。实验中浮力火焰脉动频率约为5 Hz，火焰图像和辐射热流数据采集频率取为40 Hz，完全满足奈奎斯特采样定律(采集频率>2×样本频率)。

表2 实验工况Table 2 Experimental conditions

3 结果与讨论

3.1 火焰形态特性

图7 火焰原图(a)、二值图(b)和三维重构图(c, d)Fig. 7 Flame images (a), binary images (b) and 3D reconstructed images (c, d)

图8 火焰的高度、表面积和体积随时间变化Fig. 8 Variations of flame height, surface area and volume with time kW)

图9 火焰几何参数与参数极值比随变化Fig. 9 Flame geometry parameters and ratios between maximum and minimum values vs

图10 火焰表面积热释放速率与体积热释放速率随热释放速率变化规律Fig. 10 Flame area and volume heat release rate as a function of heat release rate

表3 火焰特征参数拟合结果Table 3 Fitting results of flame characteristic parameters

3.2 火焰辐射特性

3.2.1 火焰几何参数与辐射热流

图11给出了平均火焰高度、表面积和体积随平均辐射热流的变化趋势。可以看到，平均火焰高度、表面积和体积与平均辐射热流之间存在较好的幂函数关系，且拟合指数依次增大，即火焰辐射热流与对应火焰形态的敏感性依次降低。随着与火源距离的增大拟合指数减小，即拟合指数随距离不断变化，因此很难用这种与火焰平均几何形态直接关联的方法来预测外部辐射分布，本研究将使用面元积分方法直接精确计算火焰辐射热流分布。

图11 火焰几何参数平均值随辐射热流变化情况Fig. 11 Mean flame geometric parameters vs radiative heat flux

3.2.2 面元积分的辐射热流计算

利用重构的瞬态火焰三维面元精细结构，基于辐射理论可以直接计算火焰周围的瞬时辐射热流分布。固体火焰模型假设辐射表面和被辐射表面均为漫反射表面，火焰表面发射功率均匀，并根据下式计算火焰辐射热流分布：

(3)

其中τ是空气透射率，一般取为1。FS→T是火焰与目标物的视角系数[10]。Eav=εσT4是火焰发射功率，其中ε是火焰表面发射率，取决于燃料性质和火焰厚度，斯蒂芬-玻尔兹曼常数σ=5.67×10-8W/m2·K4，T是火焰表面温度。

目前文献中尚无小尺度丙烷火焰发射率的精细实验测量值。Wan等[5]使用15 cm×15 cm油盘的热流数据和火焰辐射模型估算得到了丙烷火焰发射率。其中，热释放速率10 kW～32 kW范围内的火焰发射率均值为0.17。本文将其作为丙烷火焰的发射率。图12给出了不同热释放速率条件下火焰中心轴线过余温度随无量纲火焰高度的变化规律。使用结点直径0.374 mm 的K型热电偶测得连续火焰区的平均火焰温度为1 103 K。利用Cox和Dupuy等[21,22]提出的热电偶辐射误差修正方法，得到火焰特征温度为1 229 K。据此，结合火焰三维重构获得的视角系数就可以计算瞬时辐射热流。

图12 火焰中心轴线过余温度随无量纲火焰高度的变化Fig. 12 Centerline excess temperature versus the normalized height

图13(a)和图13(b)分别给出了15 kW热释放速率下丙烷火焰外围瞬时辐射热流的计算值与测量值。显然，不同位置的热流测量值具有较好的同步性。基于面元积分得到的瞬时辐射热流计算值很好地反映了瞬时脉动火焰辐射的波动性，且不依赖于辐射热流计的径向位置和朝向。图14给出了基于面元积分的平均辐射热流计算值与经典的点源和圆柱辐射模型计算值对比。其中，点源辐射模型和圆柱辐射模型分别根据公式(4)和公式(5)计算[2]：

(4)

(5)

图13 瞬时辐射热流计算值与测量值Fig. 13 Calculated radiation values and measured values at different heat flow meters

相比于经典的点源和圆柱辐射模型，火焰面元积分方法能够更好地预测火焰的平均辐射热流分布，火焰辐射预测值均位于±15%误差线以内。圆柱辐射模型的预测效果次之，其水平辐射热流预测值(q5)出现了明显误差。当距离火焰较远时(L/D=6.6>5.0)，点源辐射模型计算准确度较高。当距离火焰很近时(L/D=3.33<5.0)，点源辐射模型对竖直方向辐射热流的预测值偏小，而对水平辐射热流的预测值偏大。

图14 三种模型的平均辐射热流计算值与测量值对比Fig. 14 Comparison of the calculated and measured radiative heat fluxes of the three models

4 结论

准确预测火焰辐射热流分布可提升火蔓延预测精度，且为确定人员和建筑的安全距离提供基础支撑。因此，本文采用三维火焰重构方法获取了不同热释放速率下圆形火源上浮力扩散火焰几何特性和辐射热流的变化规律，得到以下结论：

1)通过三维火焰重构方法能够获得浮力扩散火焰三维几何参数的瞬态变化规律。

2)平均火焰高度、平均火焰表面积和平均火焰体积均可较好地拟合为热释放速率的幂函数。火焰表面积与热释放速率具有较好的线性关系，即火焰表面积热释放速率趋近于常数(77.6 kW/m2)。

3)平均火焰高度、表面积和体积与平均辐射热流之间有较好的幂函数关系，且拟合指数随着与火源距离的增大而减小。

4)通过面元积分的方法可以预测火焰对周边的瞬时与平均辐射热流，预测效果相比于经典辐射模型更佳，且目标物朝向不影响预测精度，该方法具有较好的适用性。