基于粒子滤波的隧道火灾烟气速度估计方法

黄琼，丁兆云

（1.中国消防救援学院 基础部，北京 102202；2.国防科技大学 系统工程学院，长沙 410073）

0 引言

随着公路、铁路和城市地铁隧道数量的剧增，隧道火灾事故发生频率大幅增加，成为危害最严重的隧道交通事故类型之一。由于隧道属于狭长密闭空间，火灾产生的热量和浓烟大量聚集，会对隧道内人员的生命安全造成严重威胁。目前关于隧道火灾的烟气运动和控制已得到研究人员的广泛关注，并对此进行了大量研究［1-3］。作为火灾烟气运动的重要参数，火灾烟气速度的准确估计具有重要意义。但由于火灾的复杂性，烟雾的传播存在诸多不确定性，高温烟气速度估计是一项困难且具有挑战性的任务。目前主要有两种方法：一是基于物理模型设计测速设备进行测量估计。比如利用皮托管根据压差完成测量［4］，但该方法不适合低流速测量，误差较大；采用热线式设备测量［5］的准确度较高，但探头易被烟气中的炭黑粒子和细小液滴堵塞损坏，测量成本较高。也有许多研究人员提出了低成本且简单可靠的测量方法［6-8］，但由于需要预先布设相关测量设备，在实际火灾应用中受到一定限制。二是基于火灾模型进行有效的火灾蔓延和烟气扩散模拟估计。目前的火灾模型有经验模型、区域模型及计算流体动力学（Computational Fluid Dynamics，CFD）模型。CFD 模型以火焰动力学模拟器（Fire Dynamics Simulator，FDS）［9］为代表，在模拟火灾时，精确度高，但对计算资源要求也高，往往需要耗费比真实火灾更多的时间实现模拟。区域模型，如综合火灾增长和烟雾传输的CFAST（Consolidated model of Fire And Smoke Transport）［10］，将火灾区域分为两区，且假设各区的物理量是均匀的，通过开口计算烟气分布及温度。在大多数情况下，现有的火灾模拟工具的准确性很大程度上依赖于用户预先设定的模型参数，如火源的热释放率、通风条件、火源位置等［11-14］和系统的初始状态。在实际火灾中，由于火灾的复杂性，影响火灾蔓延和烟气扩散的参数可能随时间变化，导致传统的火灾模型的模拟结果与实际情况有很大不同。许多学者试图结合模型与传感器数据［15-17］以改进模拟模型，但仍然难以保证实时性。近年来，数据同化技术（Data Assimilation，DA）如集合卡尔曼滤波（Ensemble Kalman Filter，EnKF）、粒子滤波器等［18-20］在许多领域广泛应用。大部分研究［21-23］基于EnKF 结合观测数据与模型，实时模拟火灾烟气的温度分布和烟层高度，能够获得较为精确的模拟估计性能，且对计算机资源的需求显著降低。本文选择粒子滤波算法模拟估计隧道火灾烟气速度，以解决目前测量成本较高、模拟精度较低与较难保证实时性等问题。实验结果表明，本文方法能够基本满足实时性要求，且具有较高的模拟精度，能为消防救援人员提供有用的关键信息，也可为排烟系统和消防规划策略提供理论依据。

1 粒子滤波原理

火灾烟气运动作为复杂的非线性动态系统，受多种不确定因素的影响。粒子滤波算法采用贯序的蒙特卡洛方法，能够对非线性的状态空间进行有效估计，适用于非线性非高斯系统。动态时变系统的状态方程和观测方程可以描述为：

其中：Xk是k时刻系统的状态量；Zk是k时刻的系统观测量；Uk和Rk分别是系统的状态噪声和观测噪声，表征系统状态估计及观测过程的不确定性程度；Fk和Hk分别表示状态转移函数和系统观测方程，用于描述系统在状态空间中的状态变化过程、系统状态与观测量间的映射关系。

1.1 火灾烟气运动的状态方程

当隧道发生火灾时，热烟气上升和卷吸的空气混合形成火羽流，受到顶棚撞击后，发生水平流动。纪杰等［24］将隧道火灾烟气的运动区域划分为：自由上升羽流阶段、顶棚径向蔓延阶段、过渡阶段和一维水平蔓延阶段。本文主要针对烟气一维水平蔓延过程中的速度估计展开。已知烟气速度在火羽流内几乎一致，但超出火羽流，急剧下降为0［25］。与火羽流中轴线距离r的气流速度可表示为：

其中：Q为热释放速率；H为隧道顶部到火源的高度。

基于离散时间序列，将k时刻的热释放速率Qk泰勒展开，如式（2）所示，同时记为常数。

结合式（2）～（4），建立烟气运动状态方程：

其中：dt为k时刻和k-1 时刻的时间间隔；系统状态Xk为通过添加噪声参数U=(U0，U1，U2)T使状态方程能够更好地反映烟气运动过程中的不确定性，提高烟气速度估计的可靠性和准确性。

1.2 基于温度的观测方程

在离开火羽流中轴线的任意距离r处，温度分布体现为：随着垂直顶棚的距离增加，温度升高，当达到最高值后，又逐渐降至环境温度T0。根据一系列全尺寸火灾实验结果，导出如下温度分布关系式［25］：

将式（1）代入式（6），得到系统观测方程，如式（7）所示，并记为常数。

为了更好地模拟隧道火灾环境下的观测过程，式（7）中引入测量噪声以表征传感器的不确定性观测误差分别为正态分布的均值和方差。监测离开火羽流中轴线任意距离处的温度值，根据式（7）建立烟气运动过程中某点温度与该点烟气速度之间的特征关系，以准确反映烟气特性间的变化规律。

基于上述建立的状态方程和观测方程，运用粒子滤波的隧道火灾烟气速度估计方法实现流程如下。

步骤6 若满足终止条件，结束流程；否则k=k+1，转步骤2，继续估计k+1 时刻的烟气速度。

2 实验与结果分析

2.1 隧道模型和测点布置

本文通过FDS 火灾模拟器获取实验数据。模拟隧道长50 m，宽7 m，高5 m，隧道两端口的属性设置为“OPEN”，即隧道端口直接连通开放空间。通过网格独立性分析，确定合适的网格尺寸，以保证模拟的可靠性和准确性。McGrattan等［26］建议取值为［4，16］，其中：δx是网格尺寸；D*为火源特征长度，因此网格尺寸设为0.5 m。隧道顶棚和侧壁的材料属性设置为“CONCRETE”。隧道内初始环境温度设为20 ℃。着火点布置在隧道纵向中心线处。选择方形作为模拟火源形状，与隧道地面齐平，燃烧面尺寸为1 m×1 m，燃料属性为正庚烷。考虑到公路隧道中小型汽车和载货卡车的火灾热释放速率（Heat Release Rate，HRR）的范围约为3～5 MW 和10～20 MW［27］，为使模拟实验更贴近真实火灾，令最大热释放速率分别为3 MW 和15 MW。隧道内设置了温度和气流速度的测点。由于竖直分布的温度最大值在距离顶棚0.02H（H为隧道顶部到火源的高度）处［24］，且系统状态方程和观测方程的建立条件为至少距离火羽流中轴线0.18H外的任意处，因此为保证算法的精确性和有效性，在隧道顶棚下方10 cm，气流速度测点布置在距离火源3 m 和6 m 处，温度测点则在火源一侧6 m 范围内每0.5 m 为间隔纵向排列。隧道模型和测点布置如图1 所示。在该模型的基础上进行4 个场景的实验，如表1 所示。

图1 隧道模型和测点布置Fig.1 Tunnel model and measuring point layout

表1 实验场景设置Tab.1 Experimental scene setting

2.2 结果分析

将隧道火灾模拟实验产生的数据，通过插值法进行预处理，并依据隧道模型参数建立系统状态方程和观测方程，再基于粒子滤波算法实现烟气速度估计验证与分析。

2.2.1 敏感性分析

在粒子滤波中，粒子数目越高，粒子多样性越好，能有效延缓退化，但必然增加算法的时间复杂度，在实时性要求较高的应用领域无法容忍。因此有必要对粒子数n进行讨论，以确定n对烟气速度估计的影响。在场景A1 下，粒子数分别为100、120，150、180、200 时，烟气速度估计的精度及计算时间，结果如图2、3 所示。

从图2 可看出，烟气速度估计曲线在n=150，180，200 时，基本趋于一致，接近于真实状态，而n=100，120 时，估计状态与真实状态偏离较远。这一点同样可以在图3 中清楚地反映出来，实验中利用均方根误差（Root Mean Square Error，RMSE）衡量估计值与真实值间的偏差。RMSE 越低，估计精度越高。可以看出n从100 增加到150，RMSE 急剧下降，而150 增加到200，RMSE 基本相当，均处于较低值。若总模拟时间为3 000 s，每10 s 获取一次观测数据，为满足实时性要求，每次获取观测数据后，将响应时间控制在10 ms 以内，即在本实验中的计算时间应小于3 s，是可接受的。由表2 可知，当n=180，200 时，计算时间超过3 s；当n=100，120，150时，计算时间少于3 s，可在每次获取观测数据后进行实时估计。综合考虑估计精度与计算时间，当n=150 时，实时估计隧道火灾烟气速度的效果最佳。

图2 粒子数对估计精度的影响Fig.2 Influence of particle number on estimation accuracy

图3 不同烟气速度的RMSE对比Fig.3 Comparison of RMSE for different smoke velocities

表2 粒子数对计算时间的影响Tab.2 Influence of particle number on calculation time

2.2.2 观测数据分析

本节从两个方面研究观测数据对烟气速度估计效果的影响如图4 所示。

图4 观测数据对估计精度的影响Fig.4 Influence of observation data on estimation accuracy

1）由于观测位置的偏差，必须对观测数据进行近似处理。如在场景A2 和A3 下，温度Ts1分别根据Ts4和Ts5、Ts3和Ts6通过线性插值方法获得，讨论近似误差对Vs1估计精度的影响，如图4（a）所示。

2）观测数据的不确定性对估计结果的影响至关重要，主要来源于传感器的测量误差。因此为解释传感器的内在特性对估计效果的影响，在场景A4 下，往观测数据中分别加入协方差为Rt1=10和Rt2=30 的噪声，这是影响粒子滤波效果的一个重要参数，如图4（b）所示。

图4（a）的结果表明，总体上估计精度能够满足要求，尤其在火灾稳定燃烧阶段；观测点位置距离估计点越近，估计值越接近真实状态值。主要原因在于观测数据的误差削弱了它在模拟估计过程中对系统状态的修正作用。下一步将研究更为有效的近似方法，以提高观测数据质量，保证较好的估计效果。从图4（b）可看出，观测不确定性程度越小，估计效果越好。由于观测器件的内在特性，远离火源的传感器的观测误差相较于火源附近的探测器的观测误差更小。因此以后的研究将对火源附近的估计点获取的观测数据进行有效合理的降噪处理，以降低它对估计精度的影响。

2.2.3 烟气速度估计效果测试

在A1 和A5 场景下实现烟气速度估计测试，如图5 所示。可以看出，在t=250 s和t=560 s 之前，估计误差相对较大，之后烟气速度的估计曲线与真实状态随时间的变化趋势基本一致，估计值在真实烟气速度的±20%误差范围内。使用平均绝对误差（Mean Absolute Error，MAE）进行评估，如表3 所示，各场景下烟气速度的MAE 均低于真实状态值的20%（ts0.2）。总体模拟估计效果较好，而且误差基本能够控制在真实值的20%以内。因为观测方程基于稳态方程建立，火灾燃烧的初始阶段为非稳态过程，在一定程度上影响了估计精度。为解决这一问题，可以通过调整观测方程中的噪声参数以减少对系统状态修正的影响。

表3 A1与A5场景下的MAE与真实状态值20%的平均值的对比 单位：m/sTab.3 Comparison of MAE and 20% of true state value in scenario A1 and A5 unit：m/s

图5 不同场景下的烟气速度估计结果Fig.5 Smoke velocity estimation results in different scenarios

3 结语

本文提出一种基于粒子滤波算法的隧道火灾烟气速度估计方法，实验结果表明，该方法可通过传感器的观测值修正实现烟气速度的精确估计。同时讨论了影响模拟估计精度的主要因素，主要结论如下：1）粒子数目反映了粒子多样性程度，粒子数越高，烟气速度估计的精确度越高，但带来的是计算时间的大幅增加，无法满足实时性要求；2）观测位置与观测的不确定性同样也会影响估计精度。观测位置越接近估计点，观测的不确定程度越低，估计性能越好，因此可对观测数据进行预处理，以提高模拟估计效果。3）在火灾稳定燃烧阶段的估计效果优于燃烧初期，但可通过调整方程中的噪声参数，使模拟估计效果达到较好水平。本文方法可用于一般隧道火灾场景，但由于实验环境建立在FDS 模拟基础上，为进一步验证估计方法在实际火灾中的模拟效果，需要再进行小尺寸或全尺寸的火灾实验，将模拟估计值与真实火灾状态值进行比较，以应用于实际工作，为消防救援和人员疏散提供有用的关键信息，同时也为排烟系统和消防规划策略提供理论依据。