基于回归分析的表面流场粒子匹配算法

陈 红，周国梁，闫 静，嵇 阳，晏成明

(1. 河海大学水利水电学院，江苏 南京 210098； 2. 广东水利电力职业技术学院，广东 广州 510635)

我国河流上兴建了大量水利、水运、交通等涉河工程，受河流复杂水沙运动作用，大中型涉河工程关键设计参数多通过物理模型试验研究确定。其中，流速是河流物理模型试验中最基本的参量[1-2]，试验工况确定、水流结构分析、模型验证或工程方案比选均依赖流速数据分析，流速测量准确性直接关系到试验成果水平。目前常用流速测量仪器有毕托管[1]、微型旋桨流速仪[3]、超声波多普勒流速仪(acoustic Doppler velocimetry，ADV)[3]、大尺度流场图像测速仪( large-scale particle image velocimetry，LSPIV)[3]、热线热膜流速仪(hot wire/hot film anemometry, HWFA)[4]、激光多普勒流速仪 (laser Doppler velocimetry/anemometry, LDV/LDA)[5]、粒子图像测速仪(particle image velocimetry，PIV)[6-10]等。物理模型试验水体含杂质较多，易阻塞毕托管或损坏HWFA敏感元件，因此毕托管和HWFA不适用于物理模型试验。PIV和LDV测量范围有限，布置烦琐，易受杂质干扰，也无法适用于物理模型试验[11-13]。微型旋桨流速仪和ADV常用于物理模型流速测量，微型旋桨流速仪操作简便、价格低廉、测量精度较高，然而，对流态干扰较大，存在测量盲区，流速较小时测量误差大。ADV测量精度高，可同时测量三维点流速和流向，但超声波在空气和水体中传播速度不同，测量时应将发射接收探头置于水面下才能进行准确测量，对流态有一定干扰，而且其价格昂贵，难以大规模应用。LSPIV通过测量水流表面示踪粒子计算对应水流质点流速，具有测量范围大、测量效率高、不干扰流场形态等显著优点，目前已广泛应用于物理模型试验。

LSPIV系统主要包括示踪粒子、图像采集和图像处理等[14]，王兴奎等[15]基于图像采集卡研发了LSPIV系统，建立了最邻近算法实现了粒子匹配，并测量了三峡工程物理模型表面流场，分析了工程前后流场变化。唐洪武等[16]改进了粒子匹配算法，采用Hopfield神经网络对流场数据进行校核和修正，并应用到快速自航船模船尾流场及航道断面流场的实时测量。陈诚[17]对示踪粒子跟随性和图像处理技术进行了深入研究，并应用于黄河模型试验中。LSPIV最核心技术是示踪粒子匹配，常规采用最邻近算法，但当粒子密度较大时，粒子匹配误差会大幅增加，降低流场测量精度，影响工程前后流速变化分析。

本文针对粒子密度大时粒子匹配误差大等问题，改进了最邻近匹配算法，提出多帧回归粒子匹配算法，在此基础上开发了改进型大尺度模型表面流场图像测试系统(ILSPIV)，并将该系统应用于八卦洲河工模型，同时与旋桨流速仪测量结果进行了对比。

1 流场图像特征

在模型水流表面播散跟随性能好的示踪粒子跟随水流运动，利用CCD摄像头记录一定时间间隔的流场图像(如图1所示)，经图像处理提取示踪粒子像，通过粒子匹配，获取t1时刻示踪粒子经过Δt时间后在t2时刻的图像，得到Δt时间内粒子位移，从而求得示踪粒子的速度，即对应水流点的速度。

图1 流场粒子图

图1流场图像中示踪粒子为白色，RGB灰度值较大。在图像中选取5个粒子像，依次读取中心点、边缘点及相邻背景点灰度值，RGB灰度值如表1所示。灰度值表明不同位置示踪粒子RGB灰度值存在一定差异，粒子中心点及边缘点处灰度值大于背景点灰度值，且随着背景点灰度值增大，粒子中心点及边缘点处灰度值也相应增大，粒子中心点、边缘点、背景点灰度值变化趋势表现一致，示踪粒子中心点灰度值与对应背景点灰度值存在一定相关关系。

表1 示踪粒子RGB灰度值

2 粒子像提取

粒子像提取是将示踪粒子从图像中分离出来，获取示踪粒子中心坐标。传统粒子提取算法主要采用阈值分割，分割算法如下式:

(1)

式中：gi(x,y)为第i帧图像在坐标(x,y)处灰度值；ki(x,y)为阈值分割后图像灰度值；T为阈值。

图2 阈值分割结果对比

阈值分割处理前景和背景差异较小的图像时分割效果较差。对图1直接进行阈值分割，处理结果如图2(a)所示。帧差法的基本原理是首先对相邻两帧图像做差分运算，再通过阈值变换消除背景、分离运动目标信息，该方法能提高背景和前景差异，有利于阈值分割，适用于运动目标和背景差异较小的图像检测，具有算法简单、运算速度快等优点。帧差法计算方法如下:

pi(x,y)=gi+1(x,y)-gi(x,y)

(2)

式中：pi(x,y)为帧差法处理后图像灰度值。

帧差法能消除背景影响，即使背景存在因外部条件导致的灰度差异，只要该差异小于背景与示踪粒子间的差异，仍可采用帧差法和阈值分割组合(帧差法阈值分割)提取示踪粒子。帧差法阈值分割先采用式(2)计算出pi(x,y)，然后将gi(x,y)=pi(x,y)代入式(1)得出ki(x,y)。对图1进行帧差法阈值分割，结果如图2(b)所示，帧差法阈值分割能有效消除模糊背景，示踪粒子提取准确度更高。

选取帧差法阈值分割处理图像，对第i帧和i+1帧图像进行帧差法阈值分割得到ki(x,y)，对第i+1帧和i+2帧图像进行帧差法阈值分割得到ki+1(x,y),对ki(x,y)和ki+1(x,y)进行和运算，得到hi(x,y)：

(3)

hi(x,y)是第i帧、i+1帧、i+2帧粒子像全部叠加在一起得到的粒子图像，背景灰度值为0，粒子灰度值为255。沿x、y方向依次读取坐标灰度值，当灰度值为255时，即可判别该坐标为粒子像坐标，记为(xm,yn),m=0,1,…,l1;n=0,1,…,l2;l1、l2分别为提取出的x、y方向上粒子像的最大数。

3 粒子匹配

示踪粒子像包含了第i帧、i+1帧、i+2帧示踪粒子像，粒子浓度较大时，无法通过最邻近匹配算法进行粒子同源匹配。根据流体运动连续性，多帧同源示踪粒子像遵循以下原则：①同源示踪粒子运动位移量较小；②同源示踪粒子运动轨迹连续光滑，不存在突变，表现为多帧同源示踪粒子像运动轨迹近似直线，回归分析时相关系数最大。本文改进了最邻近匹配算法，提出了多帧回归粒子匹配算法，算法的步骤如下：

步骤1搜索第i+1时刻同源疑似点。原则①表明示踪粒子运动范围是以初始时刻坐标为圆心，最大运动位移量d0(与运动速度相关，可由断面平均流速估计)为半径的圆，根据该原则初步确定同源点，称为第i+1时刻同源疑似点。从示踪粒子像坐标提取起始点(xm,yn)，为提高图像处理速度，设立正方形查询窗口，边长2r，在该正方形范围内提取点(xm+k,yn+h)，其中k=-r, -r+1,…,r，h=-r, -r+1,…,r，两个点的距离d1为

(4)

如果d1

步骤2搜索第i+2时刻同源疑似点。以点(xm+k,yn+h)为中心，在边长为2r的正方形查询窗口内提取点(xm+k+g,yn+h+f)，其中g=-r, -r+1,…,r，f=-r, -r+1,…,r。计算这两个点的距离d2，如果d2

(5)

总平方和为

(6)

假设上述三点拟合直线方程为

y=ax+b

(7)

(8)

(9)

相关系数Rm为

(10)

图3 同源疑似示踪粒子像

依次对其他同源疑似点与起始点(xm,yn)进行回归分析，相关系数最大的作为同源点，完成全流场粒子匹配。

图3为根据原则①处理获取的同源示踪粒子像，黑色点坐标为(x0,y0)，其同源疑似点有点(x0+1,y0-1)、(x0+1,y0)和(x0,y0+1)；点(x0+1,y0-1)的同源疑似点有(x0+2,y0-2)、(x0+2,y0-1)和(x0+1,y0);点(x0+1,y0)的同源疑似点有(x0+1,y0-1)、(x0,y0+1)和(x0+2,y0-1);点(x0,y0+1)的同源疑似点有(x0-1,y0+2)和(x0+1,y0)，则对应相关系数最高的点为点(x0,y0)、(x0+1,y0-1)和(x0+2,y0-2)，其拟合直线为

y=x

(11)

拟合直线相关系数为1，依次类推可得，黑色点与红色点属于同源粒子。

基于回归分析的粒子匹配算法开发了改进型大尺度模型表面流场图像测试系统(ILSPIV)，对图1表面流场粒子图像进行处理，得到流速矢量如图4所示，总体上矢量图错误较少，矢量与粒子轨迹重合度高。

图4 表面流场矢量图

4 实际应用

长江南京河段八卦洲左汊河道沿程阻力远大于右汊，导致左汊逐渐淤积、萎缩，既影响左汊沿岸生产生活，又影响南京长江河段整体河势稳定。因此，需要采取有效工程措施稳定八卦洲汊道河势，提高左汊分流比。为对比研究工程效果，采用物理模型研究了工程前后洲头流速分布等变化情况。根据试验要求，确定八卦洲汊道整治物理模型模拟范围，上边界为大胜关，下边界至九乡河，模拟河段总长度约43 km。选定模型平面比尺λL= 480，垂直比尺λH= 120，流速比尺λv= 10.95。

试验工况流量为15 290 m3/s，选定测流断面位于分流口，分别采用旋桨流速仪和ILSPIV对流场进行对比测试。旋桨流速仪测量数据为断面垂线平均流速，ILSPIV测量流速为表面流速，根据有关理论两者流速关系为

vc=0.85vb

(12)

式中：vc为断面垂线平均流速；vb为表面流速。

依据式(12)将ILSPIV测量的表面流速数据换算成断面垂线平均流速，图5给出了两种测试技术获得的工程前后断面垂线平均流速。总体上，两种测试技术测量的断面垂线平均流速分布形态相似，局部流速点存在一定偏差，相对误差小于2.7%。ILSPIV测量的工程前后流速变化范围为-0.25～0.22 m/s，旋桨流速仪测量的工程前后流速变化范围为-0.18～0.22 m/s，表明工程前后两种技术测量结果相对变化量基本一致，对工程影响评估无差异。

图5 两种测试技术获得的工程前后断面垂线平均流速对比

5 结 语

粒子匹配是影响LSPIV测量精度最重要的因素之一，为提高粒子匹配精度，改进了最邻近匹配算法，提出了多帧回归粒子匹配算法，从动态图像中提取连续帧，分组叠加再做减运算，采用一元线性回归，以相关系数为评估参数，确保了高粒子浓度条件下粒子的精确匹配。运用多帧回归粒子匹配算法开发了ILSPIV，将ILSPIV应用到八卦洲河工模型试验，ILSPIV和旋桨流速仪测量的断面垂线平均流速分布形态一致，相对偏差小于2.7%，工程前后两种技术测量结果相对变化量偏差小，表明多帧回归粒子匹配算法能提升高浓度粒子条件下流速测量精度，更有利于工程前后流速变化分析和评估。