流速测量技术是开展流体相关实验研究和工程应用的重要技术，对于研究基本物理现象和深入解释流动机理具有重要意义。在过去几十年中，流速测量技术不断发展，经历了从介入式到非介入式、从单点到平面再到全场、从低速测量到超声速测量的过程[1]。
近年来，基于图像处理技术的流速测量方法得到了快速发展，特别是粒子图像测速技术(particle image velocimetry, PIV)[7]。在前计算机时代，由于缺乏图像处理能力，只能通过测量匹配点和流体结构在不同时间步长的位移来获得总体对流速度[8]。随着技术的发展，纹影测速的研究和应用不断深入。纹影测速主要形成了两类算法：纹影PIV互相关算法和光流(optical flow，OF)算法。纹影PIV(schlieren PIV, SPIV)技术利用流场中自然发生的折射湍流和涡流作为虚拟示踪粒子进行速度测量，并通过互相关算法以与PIV技术类似的方式估计速度场。然而，SPIV通过窗口粒子互相关算法得到速度的空间分辨率有限，以往研究大多展示多帧内的平均结果，仅能反映主流运动速度[9]。Kegerise等[10]使用SPIV测量轴对称热羽流时，认为纹影技术应用在三维流场时存在的积分效应影响了测量结果。
Fu等[15]是最早将OF算法应用于纹影图像测速的一批研究者。他们提出了方向性平滑约束以替代传统的空间项约束，提升了流场的平滑度，但该方法在计算精度方面仍有改进空间。该方法依赖于亮度恒定假设，并且缺乏物理约束，因此在处理复杂的运动模式时，尤其是对于纹影图像，往往难以准确捕捉流场的真实运动。
为了解决这一问题，本文引入了一种具有明确物理意义的光流模型，依据流体特性简化Navier-Stokes (N-S)方程，推导出一般物理意义上动量守恒项。由适用于流体的约束项和空间平滑项组建成正则项，并运用在基本光流模型之中。利用该算法，通过获取不同高度的火焰羽流纹影图像，可以计算流速并基于线性拟合结果预测流速。结果表明，与基本光流算法(以下简称光流算法或OF算法)相比，在火焰羽流纹影图像上，本文提出的优化光流测速算法(以下简称优化算法)能够展现更加平滑的光流场，对微弱流型火焰区域的检测更加灵敏，所计算的涡量场分布也更加均匀；在火焰羽流测速实验中的统计结果更稳定、更加接近Heskestad预测模型[16]，是一种稳定和准确的火焰羽流纹影测速算法。
1 数学模型
1.1 基本模型
基本OF算法基于连续两帧来估计二维光流场。最常见的算法是在Horn和Schunck[11]方法的基础上形成的。该算法通过估计增量位移最小化的全局能量函数，获得二维光流场，包含数据守恒假设和空间平滑性假设。
其中：$ \boldsymbol{I}_x=\mathrm{d} \boldsymbol{I} / \mathrm{d} x, \boldsymbol{I}_y=\mathrm{d} \boldsymbol{I} / \mathrm{d} y, \boldsymbol{I}_t=\mathrm{d} \boldsymbol{I} / \mathrm{d} t$ I是图像点的亮度; x和y是图像点的坐标; t是时间; u和v是图像点在x和y方向上的运动分量。数据守恒项HD可以写为
其中惩罚函数通常为二次函数。在正则项中，假设速度在空间分布上具有连续性，即在小邻域中的点及其周围的点具有相似的速度，速度场在空间上变化平缓。该约束假设的表达方式是最小化速度梯度的平方，可以表示为将二次函数应用于速度一阶导数。
全局的能量函数可以写成
在大多数情况下，数据守恒假设(HD)和空间平滑性假设(HS)难以同时满足，因为这两种假设存在矛盾，导致全局能量函数最小化效果会折中，所以引入权重参数α控制每个假设的相对重要性，即α是数据项与正则项的权重参数。
基本光流模型能够很好地满足刚体物体的位移场求解，然而在面对扭曲、旋转的流体流场时，数据项在严格意义上并不是守恒的，这是由于连续的亮度在图像上展现出来是离散的[17]。
1.2 优化模型
对于正则项部分，本文保留了空间平滑项，并且考虑加入新的物理约束。
二维流动的动量守恒方程可以用N-S方程表示，对于在极短时间内的流动，本文作了一些简化。二维流动的N-S方程一般可以用如下形式表示：
在极短时间内，流动可能不会导致明显的密度和压力变化; 其次，流动是不可压缩的，流动处于高Reynolds数下，惯性力比黏性力更为显著。因此，可以忽略黏性项。此外，极短时间内的流动可以不考虑流体的加速度。式(5)可以简化为：
基于动量守恒，本文引入新的约束项标量场，
组建成新的能量方程为
其中β是空间平滑项和动量守恒项的权重系数。然后，构造Euler-Lagrange方程，使用梯度下降的方法更新光流估计值。光流增量的迭代公式如下：
其中：f′为鲁棒函数的一阶导函数，V(k)=(u(k), v(k))表示迭代次数为k的速度场，$\nabla$ 为Laplace算子，$\nabla$ 2u、$\nabla$ 2v分别为速度分量的二阶梯度。V(0)=(u(0), v(0))的初始值设置为 0，每一次迭代产生的增量都会作为下一次迭代的初始值，每次迭代产生的光流增量都会累加，直到完成迭代次数，得到当前层的光流场。最后，以原始分辨率的图像计算出整个图像尺寸下的光流场。本方法参考文[18]使用了多分辨率金字塔算法更精细地计算位移。
2 实验设计
纹影技术的原理是将偏转的光线与未受干扰的光线分开，通常使用的方法是阻挡偏转的光线，使它们不会到达测试部分成像的记录屏幕上。如果这种偏转发生在测试部分的一点上，并且偏转的光线随后被适当的刀口所拦截，图像平面上的对应点就显得更暗。在单色纹影系统中，空气密度梯度反映了火焰亮度的变化。
密度梯度具有矢量特征，其完整描述需要测量密度大小和梯度方向。典型的纹影系统只能检测该矢量的一部分，即在预定方向上的梯度(x和y方向的梯度)或者其分量，这是因为在记录的图像上，亮度是以二维网格分布的离散状态。实验装置如图 1所示。
光学部件包含：10 W的LED光源、2个相同抛物面镜(直径12英寸, 即30.48 cm)、刀口、高速相机。高速相机为XG1600MG-C-662工业相机，帧率设为500帧/s，曝光时间设为1.1 μs，图像的分辨率为960×1 024像素。
实验采用乙醇作为燃料，这是因为乙醇燃烧效率高，火焰羽流中固体烟颗粒较少，羽流几乎是透明的。
3 实验结果
3.1 点火过程
图 2给出了油盘直径30 mm下乙醇点火过程的纹影序列图。图像采集时间间隔为50 ms。
图 7所示的是优化前后的光流算法的收敛过程，图中5条曲线对应于多分辨率金字塔算法中不同的金字塔层。其中：Level1代表多分辨率图像金字塔中原始分辨率的图像，Level5则代表分辨率最低的图像。每个金字塔层内的能量残差的定义为
其中：M和N分别是该层图像水平像素值和垂直像素值，Vi, j是迭代所产生的光流增量，n为迭代次数。能量残差的物理意义是连续迭代的光流增量的平均方差，是量纲为1的量，用于评估光流估计的误差，即光流场的准确性。相比于优化前的算法，优化后的算法收敛速度很快，第1次迭代完成后能量残差就可以小于10-2。
从图 7中可以看出，两种算法在每个层级的光流场都拥有很高的准确性，优化前的算法在最终迭代出的光流场下能量残差维持10-4以下(Level1)，而优化后的算法则维持在10-5以下(Level1)，表明优化后的算法应用于纹影图像上具有更高的精度。
3.2 火焰羽流测速
其中：u为火焰中心线位置火羽流的速度，g为重力加速度，cp为空气定压比热容，ρ∞为空气密度，Qc为对流热释放速率，z为火焰正上方高度，z0为虚拟点源高度。火焰中心线并非羽流中心线，对于实际的弱羽流，羽流速度最大处应该是在火焰上方中心线附近呈正态分布；在自由羽流情况下，速度最大处出现在火焰上方中心线附近的概率最高[20]。根据纹影拍摄获得的相邻两帧图像以及间隔时间，通过优化后的算法可以计算出其光流场，即火焰羽流的位移场，结合像素距离与实际距离的关系(本研究中每个像素对应于0.366 mm)，再结合采样时间间隔，即可获得相邻两帧图像的速度场。
为了避免光流过耦合的情况，需要降低图像噪声，本文对原始速度场进行窗口平滑，平滑尺寸为1/16纹影直径大小。在平滑后的速度场选取中心区域(约4.1 cm2)计算速度，在同一工况下每隔50 ms取样计算一次，最终得到平均速度。
表 1 实验工况 |
| 参数 | 值 |
| 燃料类型 | 乙醇 |
| 室内温度/℃ | 18 |
| 空气比热容cp/ (kJ·kg-1·K-1) | 1.004 |
| 油盘直径/mm | 25, 40, 55, 65, 75, 80, 90, 100, 105, 110 |
| 热释放速率/kW | 0.19, 0.52, 0.91, 1.29, 1.69, 2.08, 2.63, 3.18, 3.92, 4.18 |
实验结果如图 8所示。本文对结果进行了量纲分析，两种方法的拟合结果都与经典Heskestad模型(式(12))较为一致。在光流算法的结果中，k=4.1，tanθ=-0.31，在优化光流算法的结果中，k=3.9，tanθ=-0.32。式(12)中，k=3.4，tanθ=-1/3。然而，优化算法的结果采用Heskestad模型拟合的拟合优度(R2=0.98), 要明显高于光流算法(R2=0.90)。
图 9展示了两种算法的计算结果与式(12)模型预测结果的对比。可以看出，光流算法在羽流区域内的结果与Heskestad模型预测结果的误差在-20%~30%之间，优化算法大大减少了误差，将误差范围缩小到-5%~20%左右，其计算结果与Heskestad模型预测结果更加接近。从平均绝对误差来看，光流算法的结果比Heskestad模型预测结果大13%左右，而优化光流算法的结果比Heskestad模型预测结果大8%左右，优化算法得到的结果更紧密地集中在Heskestad模型预测结果的周围，表明算法的准确性有所提高。这可能是由于基本光流算法在计算部分羽流纹影图像时的结果是发散的，导致统计结果偏离Heskestad模型，而优化算法则具有更好的鲁棒性和稳定性。
4 结论
本文利用纹影技术对火焰羽流进行可视化，并通过构建新的约束方程，提出了一种优化光流测速算法。该算法在火焰羽流测速方面比基本光流算法具有更好的表现：
1) 优化算法迭代得到的光流场更加平滑，对于弱流型火焰区域更加敏感，能够检测到弱小涡旋结构，可以得到更加均匀平滑的涡度场。
2) 优化算法的能量残差收敛速度很快，第1次迭代完成的能量残差就可以小于10-2，最终迭代出的光流场能量残差低于10-5。
3) 优化算法计算结果与Heskestad模型预测结果的误差范围缩小到-5%~20%左右，绝对平均误差减小到8%左右，优化算法的准确性和稳定性比基本光流算法有所提高。
尽管优化算法在羽流测速中展现了优异的性能，但其应用范围仍有进一步扩展的空间。下一步将在更多的火源类型、更大的油盘直径、更复杂的环境下进行火焰羽流速度场测量，为更准确地分析火灾动力学行为提供参考。
