荧光油膜速度场的自适应快速光流解法

蔡章博，张征宇，余皓，占书恒

1.西南科技大学 信息工程学院，绵阳 621000

2.中国空气动力研究与发展中心 高速空气动力研究所，绵阳 621000

降低摩阻意味着飞行器的油耗下降、航程增加［1-3］。基于荧光油膜的全局摩阻直接测量法，具有测量设备简单（仅需紫外激发光源、相机与镜头）、空间分辨率高、对模型物面无特殊要求等优点，已成为研究的热点和发展趋势［4-6］。该法是基于油流法［5］，在油膜中加入荧光分子，通过荧光油膜控制方程，描述荧光油膜（受紫外光激发的发光成像）灰度与其运动速度的关系，推导出基于荧光油膜的摩阻测量方程［6］。因其测量方程与光流方程［7］具有相似的形式，因此，采用光流法解算风洞试验中模型表面荧光油膜运动的时序图像，即可获得试验模型的近壁面流动结构，为掌握试验模型壁面发生流动分离的位置、分离方式与特点、漩涡形成机制等提供重要的研究数据。

自1981 年Horn和Schunck［7］提出了Horn-Schunck（HS）光流法后，Wildes等［8］引入流体力学原理率先将其用于云运动估计，2008 年Liu和Shen［9］建立了光流法和流体流动之间的数学联系。2015年，Zhong等［10］重建了机翼连接处的全局摩擦力拓扑结构；Liu等［11］研究表明，光流法较粒子图像测速法可获取更密集的流体运动速度场。2016、2017年，黄湛等［12］、王鹏等［13］先后采用变分迭代方法求解表面摩阻分布，验证了基于荧光油膜的全局表面摩阻测量技术的有效性。2019年，邹易峰等［4］通过添加人工网格线利用离散匹配，实现了模型振动与油膜运动的解耦，提高了测量精度。

但是上述光流法的数值计算中，都采用了类雅克比的矩阵分割迭代形式，存在计算复杂度高、收敛速度慢、耗时长等问题［14］，故现有光流法尚难在生产型风洞的荧光油膜试验现场实时测量并显示模型壁面流动过程，降低了工程应用的价值。

因此，高效的光流法对于基于荧光油膜的全局摩阻实时测量至关重要。尤其是随着高分辨率（可达2 500 万～6 500 万像素）的高速工业相机在荧光油膜风洞试验中的应用，使用现有光流法处理2 GB/s 的海量荧光油膜时序图像，需要在试验结束后花1～2 h 才能获得一次试验的精细测量结果。

为此，本文提出荧光油膜速度场的空间自适应快速光流解法，一方面在数值计算方式上，通过构造共轭迭代式，较现有光流法的雅可比矩阵分割迭代式有更高的计算效率；另一方面，借鉴金字塔思想，对荧光油膜图像进行空间的灰度梯度与速度场自适应降/升采样，有效降低计算的复杂度，进一步提升荧光油膜速度场的解算速度。

1 光流法的类雅可比迭代式

当前风 洞试验 中常采 用Liu等［15-17］基于Horn-Schunck 光流所改进的光流法，该方法追求一个全局能量泛函最小化的解，公式为

式中：Ix、Iy、It分别为图像灰度沿x、y 方向空间域及时域上的梯度信息；α 为全局平滑参数，该值越大代表追求更平滑的解，即颗粒度更精细的速度场；(u，v)为给定坐标点(x，y)的速度向量。通过欧拉-拉格朗日方程最小化求解可得

式中：Δ 为拉普拉斯算子，通常以有限差分形式近似Δu(x，y)=n((x，y)-u(x，y))；()为给定点坐标下n 邻域内(u，v)的加权值［7］，代入式（2）可得

后续为便于构造共轭迭代式，此处定义一个系数a ≈α 且a ≤α。当a=α时，经由式（3）构造的现有光流迭代式为

式（4）为一种类雅可比的矩阵分割迭代形式，收敛误差为

式中：eu、ev分别为(u，v)当前步和上一步的残差。

现有光流法在迭代过程中，会运用到迭代点局部邻域的信息，其中灰度梯度大的像素点在式（4）中的加权值高，即局部邻域施加的影响大。较共轭迭代式，该迭代方式所需存储量大、收敛速度慢得多［18］，其收敛速度与迭代矩阵的谱半径相关，尤其是随着高分辨率高速工业相机应用，在高分辨率条件下追求更为平滑的速度场时，相邻像素点图像灰度梯度信息变化趋于平缓，迭代矩阵谱半径增大，进而导致收敛耗时增加。

2 空间自适应快速光流解法

2.1 共轭迭代式推导

考虑到油膜速度场中速度梯度大的像素点其灰度梯度大，为此，约定(u，v)等于邻域内的速度最大值(um，vm)，有

鉴于共轭梯度法所需存储量小，具有步收敛性、稳定性高的优点［18］，为此，借鉴局部光流法中“邻域内速度具有相关性”的思想［19］有

代入式（6），有

式中：λ 的大小反映了邻域内的平滑程度，为接近0 的正值，通常取0.001～0.010。故式（3）可化为

化为Ax=b 形式，即

显然A 为一个对称正定矩阵，必然存在共轭向量ri满足

设rk+1=rk-αkApk，其中，pk为残差方向的单位向量；αk为步长。Schmidt 正交化可得

为了求解αk，构造如式（14）的二次型：

式中：x=(uk，vk)T。式（14）与式（10）同解，将f (xk+akpk)记作g(ak)，令g′(ak)=0 可得

至此，求解光流方程的共轭迭代式为

2.2 油膜图像的空间自适应采样

为提高光流法在大位移下的精度，有学者提供了网格匹配［21］、金字塔采样［22-24］的思路。本文借鉴金字塔思想，通过对荧光油膜图像进行空间的灰度梯度与速度场自适应降/升采样，以有效提高大运动估计的精度并降低计算的复杂度。为此，本文提出了基于八邻域的空间尺度的灰度梯度与速度场自适应降/升采样方法。

如图1 所示，沿(u，v)方向将给定图像的灰度梯度(Iu，Iv)矩阵按照3×3 大小进行分块，边界未能整除分块的部分以0 填充对齐，对于给定的3×3 图像块，通过自适应降采样后得到1 个像素，其灰度梯度为如此，即可解算得到空间降采样层的速度场(u′，v′)，再借助原八邻域梯度信息对(u′，v′)升采样恢复到原有分辨率下的(u，v)。

图1 自适应采样示意图Fig.1 Process of proposed adaptive sampling

2.2.1 基于八邻域的灰度梯度自适应降采样

对于给定的荧光油膜图像上3×3 图像块，为了保留八邻域内权重最大图像灰度梯度信息，可将其灰度梯度分为两类，即

式中：Iumax、Iumin分别为原始图 像(x，y)的n邻域内u 方向的梯度最大、最小值；v 方向的梯度计算同理。降采样示意图如图2 所示，在该八邻域内梯度高于均值的点更多，满足式（14）中取Iumax的条件，则取该区域内的梯度最高值12 作为采样值。

图2 灰度梯度自适应降采样示例Fig.2 Example of proposed gray-gradient down-sampling

较传统的双线性插值的降采样方法，本文方法得到灰度梯度矩阵会保留灰度梯度显著区域的信息。

2.2.2 基于八邻域的速度场自适应升采样

通过共轭光流计算，得到灰度梯度自适应降采样图像的速度场。在后续计算中需要将速度场恢复到原图像空间分辨率。鉴于传统的双线性插值方法难以在空间上升采样后保留具有突变特征的速度场，为此，提出八邻域的速度场自适应升采样法，通过原有梯度信息对应恢复邻域内速度向量，计算公式为

图3 速度场八邻域恢复示例Fig.3 Example of proposed velocity up-sampling

考虑到油膜速度场中速度梯度大的像素点其灰度梯度大，按照局部光流法理论用u、v 代替后，既减少了速度矩阵的存储量和硬件线程同步耗时，又能采用共轭迭代法提升解算效率，快速得到表达大流动结构的速度场初值。再将速度场初值按照其八邻域内点中灰度梯度分布插值，得到原图空间分辨率速度场初值，再进行全局速度场优化，即可准确捕捉到精细的流动结构，有效提高荧光油膜速度场的解算速度。

2.3 算法流程

本文算法全部流程如图4 所示，两帧时序荧光油膜图像分别对应图像Pt、Pt+1。

图4 本文方法流程Fig.4 Process of proposed method

3 仿真实验与分析

3.1 Ossen 涡对实验

参照文献［4，6］设计仿真实验，利用图5（a）、图5（b）所示的荧光油膜模拟图像表征给定的速度场为

图5 荧光油膜仿真图像及Ossen 涡对速度场Fig.5 Fluorescent oil film simulation images and Ossen vortex pair velocity field

式中：Γ 为涡核强度，本实验设定为5 000 pixel2/s；r0=30 pixel；u=10 pixel/s。在Ossen 涡对速 度场上再叠加一个匀速速度场，式（20）速度场叠加效果如图5（c）所示，对图5（a）利用双线性插值，在演化步长为0.000 1 s 的条件下演化1 000步，取间隔0.1 s 后演化至5（b）所示荧光油膜模拟图。

为验证本文算法准确性，对上述流场下的演变图像进行测试，限定迭代次数为300 进行对比实验，为了扣除边界条件对迭代计算的影响，忽略边界像素点（x<20 pixel，x>200 pixel）的结果，得到的仿真实验结果如图6 所示。如表1 所示，本文方法计算的误差更小，相比于文献［4］方法最大误差减小了2.3%，平均误差减小1.4%。从图6 也可看出本文方法较现有方法与理论曲线更为贴合，在一些速度变化较为平缓的区域本文方法效果更好。

表1 迭代300 次沿涡对分布的速度相对误差Table 1 Average errors distributing along vortex pair under 300 iterations

图6 不同方法沿Ossen 涡对分布的测量速度Fig.6 Estimated velocities distributing along vortex cores under various methods

3.2 壁面射流实验

参照文献［25］选取2 张壁面射流图像模拟流体油膜运动进行仿真实验，为方便观察，论文图像在实验图像基础上做了灰度拉伸（见图7）。为验证比较收敛速度，统计文献［25］方法收敛到给定误差限的迭代次数。壁面射流实验共设定了3 个误差 限：1.0×10-6、5.0×10-7、1.0×10-7，统计了不同方法收敛到误差限所需迭代次数及所用时间。统计结果如表2 所示。

表2 壁面射流试验在不同误差限下收敛的迭代次数Table 2 Iterations to convergence with various error limits in wall-jet experiment

图7 壁面射流图像Fig.7 Images of wall-jet

可以看出在误差限1.0×10-6、5.0×10-7的标准下，以耗时为例，本文方法效率分别提升了40.72%、26.66%。在1.0×10-7的误差限下效率提升了5.73%，这是由于后续全局优化的解算性质，随着迭代次数提高，整体的收敛速度放缓。如图8 所示，在相同的全局收敛误差限下，本文方法解得的流场特性会更加明显。

图8 同误差限下实验结果对比Fig.8 Comparison of results with same error limit

4 风洞试验

4.1 某机翼表面荧光油膜风洞试验

以某大展弦比静弹性试验机翼为对象，在2.4 m 跨声速式风洞中开展验证试验（马赫数Ma=0.82），相机分辨率为2 352 pixel×1 728 pixel，曝光时间为16 ms，镜头焦距为35 mm，采集帧率60 Hz，最终解算油膜速度场分辨率为2 000 pixel×1 200 pixel，采用本文算法处理该试验采集的时序图像，以便与文献［4］的方法进行对比。

模型迎角α=0°时，试验采集相邻两帧图像如图9（a）、图9（b）所示，本文求解到的油膜运动速度场，通过视频测量技术［26-27］可得到速度的标准单位（m/s），如图9（c）所示。

图9 机翼试验图像及解算结果Fig.9 Images and estimated result of wing test

由表3 可知，在误差限为1.0×10-6、5.0×10-7、1.0×10-7的标准下，以耗时为例，本文方法效率分别提升了103.25%、53.66%、26.17%。

表3 机翼表面荧光油膜试验在不同误差限下收敛的迭代次数（α=0°）Table 3 Iterations to convergence with various error limits in fluorescent oil film test on wing surface（α=0°）

模型迎角α=10°时，取试验采集相邻两帧图像（见图10（a）、图10（b）），误差限为1.0×10-7标准下求得的摩擦应力线图谱如10（c）所示，以图中矩形区域为例，着重与文献［4］对应结果进行比较，其油流分离现象一致，本文方法解算效率提升26.17%。相比于传统油流试验图像，本文方法得到了更丰富、清晰的流动细节。该试验表明本文方法适用于曲面试验模型。与文献［4］中算法效率比较如表4 所示，在误差限为1.0×10-6、5.0×10-7、1.0×10-7的标准下，以耗时为例，本文算法效率分别提升了9.76%、69.47%、19.89%。

表4 机翼表面荧光油膜试验在不同误差限下收敛的迭代次数（α=10°）Table 4 Iterations to convergence with various error limits in fluorescent oil film test on wing surface（α=10°）

图10 机翼试验图像及解算结果Fig.10 Images and estimated result of wing experiment

4.2 平板上圆柱绕流油流试验

在中国空气动力研究与发展中心（CARDC）的0.6 m 连续式风洞中开展了“平板上的圆柱绕流油流”试验（Ma=0.6）。图11（a）中所示为试验所用平板，红色区域为荧光油膜，图11（b）、图11（c）为相邻两帧平板上圆柱绕流引起的油膜运动图像，相机分辨率为5 120 pixel × 5 120 pixel，曝光时间为12.5 ms，镜头焦距为30 mm，采集频率80 Hz，解算油膜速度场分辨率为1 200 pixel×780 pixel。

图11 平板试验图像及解算结果Fig.11 Images and estimated result of flat plate experiment

如图11（d）所示，当平板上的来流遇到圆柱障碍时，平板上圆柱附近的壁面边界层速度开始下降、涡量开始增加，随着涡量的增加，边界层发生流动分离、不断卷起形成漩涡并试图脱离壁面时与来流相互作用，形成如图11（d）、图11（e）所示分离线，当漩涡前进一段距离再附时，将出现如图11（d）、图11（e）所示绕流附着线，其中图11（e）所示分离点S1、S2，附着点A1 对应图11（d）所标记点。以来流方向视角观察气流在S1 处发生分离时，边界层气流如图11（e）所示由S1 流动至附着点A1，经过A1 的壁面气流按图示方式向前流动并经过S2。试验结果与文献［16］试验测得流动现象相似，表明本文算法不仅能正确解算风洞试验中油流路径，而且提供了清晰的摩擦力线图谱。

5 结论

推导了构造共轭迭代式，创建迭代效率更高的共轭光流法，将解得速度场作为初值有效提高后续全局优化效率。对荧光油膜图像进行空间尺度的灰度梯度与速度场自适应降/升采样，进一步提升高了荧光油膜速度场解算速度。

1）Oseen 涡对的荧光油膜路径速度场仿真试验结果表明：本文方法在计算结果在相同迭代次数下相比现有方法收敛效果更好，限定迭代300 次的情况下，比现有方法最大误差减小了2.3%，平均误差减小1.4%。壁面射流仿真试验结果进一步表明，在同误差限下本文方法解算出来的流场特征更为清晰。

2）风洞试验中，本文方法不仅测得了正确流动现象，并且给出了清晰的摩擦力线图谱，在误差限较高标准下相较于文献［4］方法，计算效率提升了26.17%。试验结果证明本文方法相较于现有方法优势明显，工程应用价值大。