0   引　言

近年来，随着研究者们在闭环控制、流场降阶等实时问题^[1-4]上的深入研究，对PIV（Particle Image Velocimetry）图像处理算法的快速性与鲁棒性产生了更高需求。PIV图像处理算法主要分为互相关算法以及光流算法两类，互相关算法是在连续2幅时间间隔极短的图像上，假设查询窗口内的粒子运动状态一致，通过计算图像内相应位置的相关度（用互相关函数表示），获取查询窗口内粒子的平均位移，从而获得所需的速度场。Ye等^[5]利用合成流场分析讨论了互相关算法误差的影响因素，将基于现有设备改进观察方法应用于基于LIF（Laser-induced fluorescence）图像的上升气泡尾流可视化。Koroteeva等^[6]利用互相关软件进行液体流动红外热成像实验，对红外吸收液体的近壁湍流进行有效的定量分析和真正的非接触式测速。尽管互相关算法精度较高，但计算精度以及成本受到查询窗口尺寸大小、窗口之间的重叠率以及粒子位移大小的影响，实时性仍有提升空间^[7-9]。

Horn和Schunck^[10]在上世纪80年代开创性提出经典的变分光流估计模型后，光流技术被逐步应用于基于粒子的测速。光流算法利用图像匹配技术等方法，通过优化问题求解速度场，部分算法计算效率高，具有可输出单像素密集速度场等优点^[11-12]。Lu等^[13]提出了一种基于速度场分割的变分光流算法，证明该光流算法在解算非均匀流特性方面具有可靠性。Nicolas等^[14]在壁面边界的湍流流动测量中，使用光流算法获得了基于粒子图像的湍流信息，并且使用 PIV 算法对光流法进行了验证。Giannopoulos等^[15]对比分析基于标准互相关算法和基于迭代多金字塔光流算法发现，尽管光流算法准确度比互相关低5%，但更强的数值扩散性改善了子窗口未分辨梯度的计算，实现湍流中粘性子范围的更直接、稳定的起始测量。

光流算法中的LK（Lucas-Kanade）算法^[16]因其具有的天然高并行结构，可在GPU（Graphic Processing Unit）中进行计算，与在CPU（Central Processing Unit）中计算相比，能达到至少50倍的加速效果^[17]，而LS（Liu-Shen）算法^[18]可描述流体运动，具有更高的精度及鲁棒性，还可以与其他光流算法耦合^[19]。李成美等^[20]将LK光流算法与三帧差分法结合，并且在图像处理中引入了运动补偿，通过实际应用发现该方法能够有效地估计出目标物体的运动信息。Lin等^[21]研究了LK、PIV和PTV（Particle Tracking Velocimetry）三种方法对于使用热像仪和热示踪剂的组合进行表面速度测量的性能，发现与经典的互相关方法及其变体相比，LK算法考虑了连续图像面块（如剪切和旋转）的仿射变形，使得结果更加准确和稳健。Luís P. N. Mendes等^[22-23]提出了流体力学应用的光流方法基准研究，评估Lucas-Kanade、Liu-Shen、Farnebäck/Liu-Shen三种光流方法的性能，LK与LS 组合的性能在所有流量类型、图像条件和图像位深度下都是最佳的。在Luís P. N. Mendes等^[24]的最新研究中提出了基于相关性的粒子图像测速（PIV）和LK和LS光流方法的混合的新型测速算法，发现混合PIV技术合成图像的PIV密度增加了3倍，表明混合PIV 技术尤其是集成 Lucas-Kanade方法，可以扩大PIV 可实现的湍流测量范围，具有广阔的应用潜力。

综上所述，大多数研究针对光流算法的性能进行了系统测试，并着重注意了算法的实时性，但大多测试算法类别单一，且较少在钝体绕流实验流场中进行验证。本研究以钝体绕流流场作为应用背景，采用GPU加速的LK光流算法和LK&LS光流算法进行测试，算法在标准粒子图像中进行验证后，利用拖曳水槽装置搭建实验流动环境，测试算法性能，并与互相关算法进行比较。

1   光流算法及有效性验证

1.1   基于GPU的光流算法及实现

1.1.1   LK算法

LK算法本质一种基于窗口匹配的算法，通过查询窗口的图像变形，使得2张图像（模板图像和目标图像）匹配，该算法需要满足亮度守恒假设、空间一致性假设和时间连续性3个基本假设条件。所有的LK算法均采用叠加式（Additive Methods）的方法，即算法在迭代修正速度时，在原有速度上增加微小量作为新的速度。根据模板图像与目标图像建立的映射关系的方向不同，可简单将LK算法分为3类，表示如下：

1）正向图像变形技术，用符号LK-F表示，原始公式如下所示:

2）逆向图像变形技术，用符号LK-I表示，原始公式如下所示:

3）对称图像变形技术，用符号LK-S表示，原始公式如下所示:

式中，I⁰和I^t分别代表0和t时刻位于(x, $y $)处的图像亮度，W(x, y)为窗口权重函数（即高斯函数），查询窗口IW大小为(2r + 1) × (2r + 1)，r为半径。

第一次迭代计算时设定初始速度估计${\boldsymbol{u}}_{0} $设置为0，通过泰勒展开对原始公式变形，并进一步计算出修正因子，即速度场${\text{Δ}} {\boldsymbol{u}} $。LK算法中一般通过设定有限迭代次数或者根据多次迭代后的${\text{Δ}} {\boldsymbol{u}} $大小作为收敛条件。当达到最大迭代次数或者||${\text{Δ}} {\boldsymbol{u}} $|| ≤ ε（ε为设置的某个较小的值）时，迭代终止。为同时满足计算耗时和计算精度的要求，设定三次高斯-牛顿迭代后即可终止。计算得到速度场${\boldsymbol{u}} $后，修正速度估计${\boldsymbol{u}}_{0} $，使得${\boldsymbol{u}}_{0} $不断逼近${\boldsymbol{u}} $，当满足收敛条件后，迭代终止。

以正向变形为例，假设初始速度估计为${\boldsymbol{u}}_{0} $，且初始速度${\boldsymbol{u}}_{0} $是${\boldsymbol{u}} $的近似值，即${\boldsymbol{u}} $−${\boldsymbol{u}}_{0} $ ≈ 0，为了公式简洁，${\boldsymbol{u}}_{0} $−${\boldsymbol{u}}_{0} $用${\text{Δ}} {\boldsymbol{u}} $表示。使用泰勒展开对I^t(${\boldsymbol{x}} $ + $ {\boldsymbol{u}} $)在${\boldsymbol{x}} $ + ${\boldsymbol{u}}_{0} $处进行逼近，取一阶项并舍去高阶项，可以将方程组变形为如下所示：

其中，对图像光强的空间梯度展开为$\nabla {I}^{t}{(\boldsymbol{x + }{\boldsymbol{u}})} $ = [$\partial $I/$\partial $x,$\partial $I/$\partial $y]^T。${I}^{t}(\boldsymbol{x + }{\boldsymbol{u}}_{\text{0}} ) $为根据初始速度进行变形估计得到的中间辅助图像，利用最小二乘法，可变形为如下所示方程组：

式中，对角矩阵W = diag[W（${\boldsymbol{x}} $₁）,W（${\boldsymbol{x}} $₂）,…,W（${\boldsymbol{x}}_{\boldsymbol{n} }$）]。方程组（5）常常被简写为：

式中，${\boldsymbol{H}} $为2 × 2的伪海森矩阵，在编写程序时可以利用像素乘法和卷积对矩阵${\boldsymbol{H}} $和矩阵${\boldsymbol{c}} $的分量进行求解。因此，将原式用卷积可以表示为：

式中，*代表卷积符号。上式是通过对${\text{Δ}} {\boldsymbol{u}} $求解后，更新初始速度估计${{{\boldsymbol{u}}}_{0}} $ + ${\text{Δ}} {{\boldsymbol{u}}} $→${{{\boldsymbol{u}}}_{0}} $，多次修正速度估计${{{\boldsymbol{u}}}_{0}} $后，使得${{{\boldsymbol{u}}}_{0}} $不断逼近${\boldsymbol{u}} $，从而求得最优速度场${\boldsymbol{u}}$。

为了提高迭代计算的鲁棒性，使用${\boldsymbol{u}} $ − ${{{\boldsymbol{u}}}_{0}} $替换${\text{Δ}} {{\boldsymbol{u}}} $，即变形为：

在迭代计算时，利用式（1）获取速度场$\boldsymbol{u} $后，将其作为下次迭代中的修正速度${\boldsymbol{u}}_{\boldsymbol{0}} $，满足一定的收敛条件后迭代终止。

1.1.2   LS算法

Liu等^[18]通过方程建立了流体运动与三维对象在图像平面上的投影的联系，可用于各种流动可视化中，即LS光流算法。其推导出的基于物理的光流方程如下：

式中，方程左侧项中，$I $为图像亮度，$\boldsymbol{u} $ = （u, v）^T为相机拍摄物体在图像平面中的投影的速度，$\nabla $即求图像空间梯度，$\boldsymbol{x} $指图像坐标矢量。方程右侧项代表了由标量ψ（与可视化介质相关）引起的亮度变化。计算采用的LS算法使用HS约束，即变分光流，泛函表达形式如下：

泛函中，α为拉格朗日乘子。引入任意光滑函数，通过在图像边界上施加黎曼边界条件，并利用欧拉-拉格朗日方程可求解该问题，因此可得：

求解，可得：

式中，h代表空间步长，大小为1。该线性方程可用雅可比迭代进行求解。

1.1.3   图像金字塔结构

求解PIV图像时，同一速度下，相机采集频率越低，跨帧图像中PIV粒子位移越大。针对大位移运动问题，需要将原算法耦合至高斯图像金字塔结构中，以3层图像金字塔为例（图1）：

图  1  高斯图像金字塔结构

Fig.  1  Gaussian image pyramid structure

下载: 全尺寸图片 幻灯片

1）第一层放置原尺寸的粒子图像对，使用高斯核对图像卷积并降采样（本文为去除偶数行和偶数列），重复两次，获得长宽为原始尺寸一半和长宽为原始尺寸四分之一的小尺度图像，分别将其作为金字塔结构的第二层和第三层（即最顶层）。

2）在图像最顶层，对小尺寸的模板图像和目标图像进行变形操作，设置初始速度场为0，使用LK算法迭代计算速度场（u₃, v₃），该速度场主要代表了大位移的粒子运动。

3）对步骤2最顶层图像计算出的速度场（u₃, v₃）进行升采样（即将速度变为顶层速度场的两倍，并将分辨率扩展至第二层图像的尺寸大小），以升采样后的速度场作为初始速度，使用LK算法迭代计算得到速度场（u₂, v₂），该速度场含有更低频的信息。

4）与步骤3类似，对速度场（u₂, v₂）升采样后，将其作为初始速度，使用LK算法迭代计算得到最终的速度场（u₁, v₁），该速度场分辨率与原图像一致。

与其他光流算法不同，本算法将速度增量作为迭代结果，并根据速度增量算出最终速度，即为未后处理的高分辨速度场。经过中值滤波、高斯滤波等后处理手段，并耦合Liu-Shen算法后，可得到更高品质的速度场。

1.2   基于标准粒子图像的有效性验证

测试所使用粒子图像数据集来自Cai等^[25]圆柱绕流数据集，所有算法均在CPU型号为i5-10400F以及GPU型号为Tian V的工作站上运行测试。图2给出了圆柱绕流的合成粒子图像对和真实速度场矢量。测试采用的粒子图像大小为256像素 × 256像素，用于图像处理部分的大小为150像素 × 256像素。雷诺数（Re = ρuL/μ，L是特征长度，即圆柱直径）在150～300之间，图像总数量为900。LK算法设置中，窗口半径为9像素，金字塔层数为3层。

图  2  圆柱绕流测试结果

Fig.  2  Test results of cylindrical flow

下载: 全尺寸图片 幻灯片

根据圆柱绕流粒子图像数据集，验证LK算法和LK&LS算法的精度及耗时，并与互相关算法进行对比。互相关算法来自OPENPIV中的windef函数，查询窗口大小为64像素 × 64像素，32像素 × 32像素和16像素 × 16像素，窗口重叠率为50%，算法迭代中通过中值滤波去除异常值。采用的误差判定准则包括均方根误差（$ RMSE = \sqrt {\dfrac{1}{{mn}}\sum\limits_{x,y \in \Omega } {{{\left\| {{\boldsymbol{\,u}_{pre}}(\,x,y\,) - {\boldsymbol{u}_{ref}}(\,x,y\,)}\, \right\|}^2}} } $，Root Mean Squared Error）和角误差（$ A\,A\,E = \dfrac{1}{{mn}}{\,{\mathrm{c\,o\,s}} ^{ - 1}} \left( {\dfrac{{{\boldsymbol{u}_{pre}}(x,y){\boldsymbol{u}_{ref}}(x,y)}}{{\left\| {{\boldsymbol{u}_{pre}}(x,y)} \right\|\left\| {{\boldsymbol{u}_{ref}}(x,y)} \right\|}}} \right) $，Average Angular Error）。

如表1所示，根据RMSE以及AAE的结果，LK算法对于求解理想情况下低雷诺数的圆柱绕流要优于互相关算法，同时由于算法在GPU中实现，LK算法计算耗时明显更小。不同图像变形下的光流算法精度从低到高依次是LK-I、LK-I&LS、LK-F、LK-F/LS、LK-S及LK-S&LS算法，计算耗时由低到高依次是LK-I、LK-F、LK-S、LK-I&LS、LK-F/LS及LK-S&LS算法。LK算法在耦合LS后，RMSE及AAE均有小幅度的降低，且计算耗时进一步提高。

表  1  圆柱绕流合成粒子图像测试结果

Table  1  Results of composite particle image test for cylindrical flow

算法名称	均方根误差/像素	角误差/(°)	计算耗时/ms
LK-I	0.0829	1.806	15（GPU）
LK-I&LS	0.0813	1.784	88（GPU）
LK-F	0.0797	1.780	23（GPU）
LK-F/LS	0.0774	1.753	96（GPU）
LK-S	0.0760	1.670	37（GPU）
LK-S&LS	0.0730	1.561	107（GPU）
互相关算法	0.1277	3.131	1502（CPU）

下载: 导出CSV 

| 显示表格

图3分别给出了理论值、LK算法（以LK-S为例）、LK&LS算法（以LK-S&LS为例）和互相关算法的圆柱绕流的流向速度云图及横向速度云图结果。互相关算法的空间分辨率与查询窗口大小及窗口重叠率相关，求解得到的流场大小仅18像素 × 31像素；LK以及LK&LS算法均可以逐像素输出得到密集的速度场和涡量场，速度场大小为150像素 × 256像素，流动细节相比互相关算法更丰富，在流向速度云图中可以较为清晰的显示出回流区的摆动。但单一的LK算法的空间分辨率仍与查询窗口大小相关，LK&LS与单LK算法相比，流向速度场及横向速度场均更为光滑，异常值也较少，同时更接近真实解。此外，这几种算法在边界处均稍有失真，互相关算法的失真程度更大。

图  3  圆柱绕流速度云图

Fig.  3  Velocity cloud image of cylindrical flow

下载: 全尺寸图片 幻灯片

图4是为不同粒子浓度、不同粒子半径和不同高斯噪声下算法的均方根误差。综合来看，LK算法对粒子半径和粒子浓度敏感性最高，其次是互相关算法，在大部分合适的工况下，LK算法精度要优于互相关算法，而在测试的几乎所有工况下，LK-S&LS算法拥有最高的鲁棒性和精度。其中，LK-I算法仍是所有算法中最快的算法。

图  4  不同算法的均方根误差

Fig.  4  Root-mean-square error of different algorithms

下载: 全尺寸图片 幻灯片

2   基于拖曳水槽的实验测试

2.1   实验设置与方法

拖曳式水槽实验平台主要由电机、导轨、水槽、激光、相机等组成。如图5所示，水槽为有机玻璃制成，长宽高尺寸为2400 × 500 × 400 mm，厚度为20 mm，水深保持在300 mm，水温为18℃。2400 mm长的水槽能保证特征长度为10 mm的钝体在水槽中运动时（低雷诺数下），行程长度足够流场充分发展。水槽支撑结构采用6060L标准铝型材，具有强度高、不易氧化且耐腐蚀能力强等特点，实验中PIV流场测量基于该拖曳式水槽实验平台进行。

图  5  PIV测量系统

Fig.  5  PIV measurement system

下载: 全尺寸图片 幻灯片

实验中流场显示测量系统主要由以下几部分组成：激光器（Photontec Berlin公司MGL-N-532A DPSS Laser）、示踪粒子（Dantec Dynamics公司，S-HGS-10 Silver Coated Hollow Glass Spheres Description）、工业相机（维视数字图像技术有限公司，MV-HS130G）、后处理工作站。激光器工作模式为连续激光，波长大小为532 ± 1 nm，功率4 W，激光颜色为绿色，激光强度连续可调。通过线性激光产生器和30°透镜产生相应的约1 mm厚的激光平面光源。示踪粒子的平均直径大小为10 μm，密度与水密度接近，不易沉降，熔点为 740℃，材料为硼硅酸盐玻璃，通过表面镀银增加了整体的反光性。工业相机像素尺寸为4 × 4 μm，其最大拍摄像素为1024像素 × 1280像素，两帧图像时间间隔为0.05 s，单次采集时间为5 s。

钝体模型选为圆柱体，材料采用有机玻璃制作，表面喷涂哑光漆。圆柱的特征长度统一取L = 10 mm，均垂直于拖曳速度放置，如图5（c）所示，装置阻塞比（钝体特征长度与水槽宽度的比值）为0.02，该阻塞比下钝体几乎不受壁面影响。实验采集了在雷诺数Re = 100、150和200下圆柱尾迹处的 PIV 粒子图像，如图所示，由于激光平面厚度约为 1 mm，可将其近似为二维流动。

2.2   钝体绕流流场快速计算

为获取时均流场，以约3个涡旋脱落周期下的瞬时数据作为统计量（包括LK-S、LK-S&LS和LS算法），在获得沿y = 0处的平均速度曲线后，与其他文献中的结果进行对比。图6（a）为LK-S算法下不同雷诺数的圆柱绕流平均流向速度场，图6（b）为平均流向速度场中沿y = 0的速度曲线，3种不同雷诺数下，3种算法处理得到的沿y = 0处的轴向速度分布重合度均较高，其中Re = 150时，在x/L > 6处互相关算法的曲线与LK-S及LK-S&LS算法出现小幅度分离。3种算法与Qu^[26]、Williamson^[27]和Persillon^[28]等的三维仿真或者实验结果进行了对比，不同雷诺数下处理得到的速度与文献结果有较高的一致性，这进一步证实了所搭建的拖曳式水槽装置及配套的图像测速系统产生的误差对实验结果影响较小。

图  6  不同雷诺数下算法处理结果

Fig.  6  Algorithm processing results under different Reynolds numbers

下载: 全尺寸图片 幻灯片

图7中，从x/L = 2的无量纲速度分量u和v的曲线（Re = 150，T = 0时刻）来看，LK-S与LK-S&LS算法无明显差异。其中，和流向速度u相比，横向速度v的量级较小，因此互相关算法与光流算法在横向速度之间的差异较大，特别是在y/L = −0.46附近，受限于分辨率等原因，互相关算法无法准确的得出峰值速度。由于该位置在回流区最远点附近，因此最小流向速度u接近于零，且光流算法与互相关算法速度场分量u基本相同，但是LK-S以及LK-S&LS算法更光滑。从x/L = 4的u和v的曲线（Re = 150，T = 0时刻）来看，可得到相似的结论。从速度分量u的曲线来看，在速度梯度较小的区域（−2.5 ≤ y/L ≤ −1），LK-S和LK-S&LS算法比互相关算法的速度型更光滑，从速度分量v的曲线来看，在速度峰值处（y/L = 1.3），互相关算法出现异常点，而LK-S和LK-S&LS算法在峰值处速度曲线更为精细且光滑。

图  7  速度剖面（Re = 150，T = 0）

Fig.  7  Velocity profile （Re = 150, T = 0）

下载: 全尺寸图片 幻灯片

图8表示Re = 150时，单个涡脱落周期T内4个瞬时时刻的涡量等值线图（尾涡完全发展后呈周期性脱落）。互相关算法、LK-S算法和LK-S&LS算法得到的涡量等值线图中，均可以看到圆柱尾部由上下剪切层依次卷起和发展形成的旋涡，涡旋交替脱落且轮廓较为完整。经对比可以发现，互相关算法计算得到的涡量等值线图分辨率低，且存在较多噪点（异常涡量值），LK-S算法与LK-S&LS算法由于输出的密集速度场，得到的涡量场分辨率更高，噪点更少，LK-S&LS计算得到的涡量场中异常涡量值数据是3种算法内最少的，涡旋的形貌也更为清晰。

图  8  圆柱绕流瞬时涡量等值线图（Re = 150）

Fig.  8  Isogram of instantaneous vorticity around a cylinder （Re = 150）

下载: 全尺寸图片 幻灯片

3   结　论

本文采用图像变形技术和图像金字塔技术，并将LK算法与描述流体物理运动的LS算法耦合，开发了基于GPU的多类光流算法程序，实现了流场的高精度实时显示，研究结论如下：

1）开展了基于人工合成粒子图像的光流算法有效性验证。使用具有真实值的粒子图像，评估了在圆柱绕流中光流算法和互相关算法的计算精度及效率。结果表明，LK-I算法具有最高的计算效率，耗时为15 ms，LK-S&LS算法的RMSE和AAE为0.073和1.56，具备更高精度和鲁棒性，抗噪声能力最强。

2）开展了基于小型拖曳水槽的光流算法实验测试，利用LK光流法和LK&LS耦合算法快速获取近尾流区速度场及涡量场，流场光滑且涡旋清晰明显，其中LK&LS算法的流场噪点最少，LK光流算法其次。

3）开展了光流算法与互相关算法的性能对比测试，精度最高的LK-S&LS光流算法RMSE和AAE为0.073和1.56，精度最低的LK-I算法RMSE和AAE为0.083和1.80，互相关算法的RMSE和AAE为0.128和3.13，且在实验测试中互相关算法流场噪点明显更多。因此在钝体绕流流场下，光流算法优于互相关算法。

本文在标准粒子图像的光流算法有效性验证中，未考虑粒子浓度和粒子直径相互之间的影响，未来研究可考虑通过不同粒子直径和粒子浓度两两组合研究对算法精度的影响。另外，拖曳式水槽实验中流动显示与测量技术手段较为单一，后续可以使用纹影法、热线等其他测量手段，进行进一步的误差分析，并验证本文的光流算法应用在纹影图像上的可行性。