基于图像处理的颗粒羽流特性试验研究

江晨琦，刘 磊，2，李世舟

（1.上海交通大学a.海洋工程国家重点实验室；b.船舶海洋与建筑工程学院，上海 200240；2.上海交通大学三亚崖州湾深海科技研究院，海南三亚 572024）

0 引言

含颗粒的两相羽流和射流常见于各种自然及工程流动中，如深海采矿作业中的尾矿排放、工厂含颗粒废水排放、海洋疏浚作业、火山喷发等。两相流动时颗粒会在环境中扩散，对生态环境产生不利影响［1-2］。了解羽流中颗粒的扩散速度、浓度等特性，对于预测和控制羽流排放具有重要意义。

目前，对于单相羽流及射流的理论及试验研究已比较成熟［3］，但对于多相羽流尤其是复杂环境下（如横流、分层流等）的多相羽流特性研究尚不深入。国内外学者对两相羽流开展了试验研究。Jiang 等［4］使用三维粒子图像测试（PIV）技术研究了垂直含沙羽流的颗粒速度特性，但对于含较细颗粒（小于100 μm）的羽流，图像处理很难区分PIV粒子和羽流颗粒。Lee等［5］研究了水平含沙羽流特性，采用图像局部颗粒计数的方法确定羽流中颗粒浓度。Hall 等［6］采用接触式的光学探针测量垂直羽流颗粒的速度和浓度，初步研究颗粒速度和浓度特性、羽流扩散率受初始参数的影响以及颗粒动量通量的变化，但对于试验尺度羽流，此类接触式测量方法精度有限并会对流场造成影响。Mcconnochie 等［7］研究了染色流体中浮力方向对羽流扩散率的影响。James等［8］通过染色对羽流中的颗粒和流体分别进行可视化，基于图像结果研究了颗粒羽流在横流中的轨迹及颗粒与羽流流体的分离现象，并与单相羽流结果进行对比，但这类方法成本较高，并且难以获得羽流内部的颗粒速度特性。

设计了模拟固液两相羽流排放及扩散的试验装置，分析了基于高速图像采集和计算机图像处理的羽流扩散特性，并提出了颗粒浓度及速度的测试分析方法，用于静止及流动条件下羽流特性的准确提取。

1 试验装置和试验流程

1.1 试验装置

羽流排放及扩散模拟试验装置如图1 所示，主要由羽流发生装置、输送管路系统、标定池及试验水槽组成。羽流发生装置包括直径0.6 m、高0.4 m 的圆柱形玻璃桶和双桨叶搅拌装置，顶部设有注水口及颗粒投放漏斗。羽流通过管路系统被输送至试验水槽，管路设有流量计，用于实时测量排出流量。标定池用于确定羽流的初始浓度。

图1 试验装置

试验水槽长2.5 m、宽0.6 m、高1.0 m，为全透明玻璃材质。水槽上部设有拖曳系统，可拖曳羽流排放口做匀速直线运动，以模拟均匀横流中的羽流排放［8-9］。试验中，可根据试验目的改变排放深度、角度、拖曳速度以及排放口的直径、形状，研究不同排放条件下的羽流扩散特性。

试验中采用光学无接触式测量分析方法，图像数据使用SONYa7M3 4K 全画幅微单数码相机采集，分辨率为2 160 ×3 840；进一步地，为捕捉局部颗粒的运动特征，采用Photron FASTCAM 系列高速摄像机进行图像采集，配备105 mm 微距镜头，分辨率为1 024 ×1 024。

试验颗粒可根据试验需求进行选取。本试验采用密度约为2 500 kg／m3的白色玻璃微珠，粒径范围为100～150 目（目指每英寸筛网上的孔眼数目），如图2所示。

图2 模拟羽流颗粒

1.2 试验流程

试验流程如图3 所示。向试验水槽中注水至水面没过排放口5 cm。在发生装置中加入适量颗粒和水，开启搅拌装置至稳定转速，持续搅拌至混合物均匀、稳定。打开管路阀门，启动输送泵，待管路中羽流流量和浓度稳定后，将混合物引入标定池，标定排放羽流浓度和流量，并与流量计测量结果对比校验。标定完成后，关闭标定池阀门，打开羽流出口阀门，将混合物排入试验水槽，开展试验。排放达到稳定状态后，采用数码相机或高速相机采集羽流图像数据，流量计记录流量数据。对于每个工况，进行数次重复试验以验证重复性。待试验结束后，基于图像处理方法对采集的图像数据进行处理，得到羽流中颗粒的扩散及运动特性。

图3 试验流程

2 测试分析方法

基于非接触式的光学图像处理方法，以静止环境中向下排放的垂直颗粒羽流为例，阐述羽流扩散特性、颗粒速度和颗粒浓度的图像分析方法。图像处理算法均通过Matlab软件实现。

2.1 羽流轮廓

羽流轮廓是基于时均图像灰度阈值的方法进行的。以典型试验工况（羽流出口管径D＝8 mm、平均初始流速u0＝1.92 m／s、初始质量浓度c0＝15 g／L）为例，首先对摄像机拍摄得到的连续图像进行灰度值时间平均，然后将计算获得的图像去除背景灰度，最后通过二值化和边缘检测算法得到羽流的视觉边缘及由视觉边缘包络的羽流视觉半宽b（即羽流在某一高度处的半径），如图4 所示。

图4 颗粒羽流形态示意图

在羽流分析中，较常使用的是速度高斯半宽bp和浓度高斯半宽bc。在某一羽流高度z上，bp定义为流场速度uf＝（1／e）ufm上的点与中心轴线的距离，bc定义为颗粒浓度c＝（1／e）cm上的点与中心轴线的距离，其中ufm和cm分别表示羽流中心轴线上的流速和颗粒浓度。速度高斯半宽bp与视觉半宽b的关系可以表示为：2 种高斯半宽的比值λ ＝bc／bp通常取1.2［10］。理论上，在静止环境下羽流视觉半宽随高度线性增加［10］。因此，根据图像各高度处羽流边缘的像素点坐标进行线性拟合，得到羽流的视觉半宽扩散率αb＝db／dz，进一步按照比例关系得到速度高斯半宽扩散率αp＝dbp／dz。

图像处理过程及结果如图5 所示，其中x为羽流视觉边界与中心轴线的距离。

图5 羽流轮廓提取图像处理

2.2 羽流浓度

为确定羽流初始浓度，将均匀、稳定的混合物注入标定池内进行标定。标定池出口下方设置长方体水箱以及孔径远小于试验用颗粒的薄滤网。将混合物注入标定池数秒后关闭，记录水箱中混合物总体积；取出滤网中的颗粒进行清洗、烘干、称重，获得颗粒总质量，进而获得混合物的初始质量浓度。

利用我院自制的调查问卷对两组孕妇进行调查，主要包括：（1）孕妇的基本信息，如：年龄、民族、身高和体重，同时还包括文化背景、职业、孕产次和收入水平等。（2）孕前半年至孕前3个月的日常生活、工作、居住环境，同时还包括是否发生负性生活事件和感冒[3]。孕前期是否发生以下感染性疾病，如：慢性疾病和生殖系统感染，还需对孕妇的药物服用史进行和饮食状况进行掌握。营养摄入为：乳制品、蔬菜、肉、水果，在此基础上还包括叶酸和维生素。（3）对围产儿的具体状况进行了解，如：胎儿状况、分娩情况和诊断结果等[4]。

羽流浓度基于时均图像进行测试分析。Mingotti等［11］在对含颗粒的固液两相流浓度进行图像处理时指出，在相同视场宽度下，颗粒总数应与光强（图像灰度）线性相关。对于三维圆形羽流，在不同高度和不同水平位置羽流颗粒所占的厚度均不相同。由于采集的图像仅呈现二维信息，因此可将图像中二维信息依据高斯分布假设进行转换，获得该视场中的羽流中心浓度或平均浓度，如图6 所示。以羽流边缘延长线与中心轴线的交点为原点O，z轴与羽流出口方向平行，向下为正；x、y为水平坐标，x、y轴分别与相机视场方向垂直和平行。cm0表示x＝x1、y＝0 处的颗粒浓度。

图6 浓度图像处理信息转换

根据三维多相羽流的高斯浓度分布假设［12］：

式中，r为该点距离羽流中心的位置。由于本试验采用的是120～250 μm的细颗粒，在质量浓度小于20 g／L条件下重叠不明显，因此可以认为时均图像上某点的灰度值与该处视场中的颗粒总数正相关。根据高斯分布函数的性质，在任意高度z上，有

根据式（2），对于图像中颗粒相范围内的任意一点，该点所在高度的羽流视觉宽度2b范围内的视场平均浓度与该位置中心最大浓度cm之比为一定值A，故将图像各点的灰度值除以该高度处的羽流视觉宽度作为该处的相对粒子浓度强度，从而得到全局相对浓度分布，如图7 所示。根据标定得到的羽流出口浓度c0，可以得到全局的颗粒浓度分布。

图7 羽流浓度图像处理结果

浓度处理结果通过Lai等提出的颗粒羽流微分方程模型［13］并结合抽吸法［14］进行验证。抽吸法是一种精度较高的接触式浓度测量方法，如图8 所示。将内径为5 mm、外径为7 mm的透明抽吸管伸入待取样位置，在羽流试验进行的过程中同步抽吸到烧杯中，通过后续质量分析得到浓度。Lee［15］的研究表明，当抽吸管角度与流速平行、速度比（即抽吸速度us与流速uc的比值）q＝us／uc＝1 时效果最佳，且当0.8 ＜q＜2.0时，测量误差均在10%以下。垂直颗粒羽流验证试验中，抽吸管与羽流排放方向平行布置，蠕动泵抽吸速度的选择根据羽流在该处的流场时均理论速度值确定，使得q接近于1，以保证抽吸效率。

图8 抽吸法测颗粒浓度

以表1 中的工况3 为例，羽流充分发展段图像处理得到的相对浓度分布结果如图9 所示。可以发现，在羽流中心轴线上，图像处理得到的浓度分布与微分方程模型结果、抽吸法单点测量浓度（z／D＝10，20）均吻合良好；在各水平高度位置，图像处理结果均与高斯分布拟合结果吻合良好，同时拟合得到的高斯半宽值与基于视觉半宽确定的高斯半宽值具有良好的一致性，如图10 所示。

表1 试验工况参数

图9 羽流浓度图像处理结果及验证

图10 2种方法所得速度高斯半宽值对比

2.3 颗粒速度

羽流中颗粒的运动测量是基于高速相机采集和改进的粒子追踪测速（PTV）图像处理算法。由于羽流中颗粒速度的三维高斯分布及紊动特性，通过视场中特定位置处的颗粒速度并不稳定，而且整体向量场较为稀疏，因此应收集一定区域内的所有速度值来确定颗粒速度。将高速相机采集的结果逐帧提取，按8 帧为一组进行后续处理，处理流程和结果如图11 所示。具体过程为：①将灰度图像使用灰度阈值去除景深外的虚化颗粒，并通过滤波减小光照不均及颗粒反射对颗粒识别的影响；②进行自适应二值化；③较低浓度下对二值图像直接进行连通图提取，较高浓度下使用基于灰度图像的分水岭算法对重叠颗粒进行分离；④根据连通图所占像素点数，按一定阈值筛选符合条件的连通图作为颗粒，并根据连通图像素点平均坐标计算颗粒中心坐标；⑤对相邻帧使用相关系数法进行颗粒匹配；⑥若某一颗粒在8 帧内能够连续捕捉，则利用标定像素长度与位移、时间关系，计算颗粒通过该区域的速度值；⑦记录一定时间内通过该区域的所有颗粒的速度值。

图11 图像速度处理过程

对羽流势核段外（充分发展段）的局部区域进行拍摄，监测范围为3.35 cm×3.35 cm区域。部分结果如图12 所示，其中f为颗粒速度的柱状频数，uz为颗粒垂向速度，颗粒的平均沉降速度为0.41 m／s。

图12 羽流充分发展段颗粒速度图像处理结果

进一步地，在羽流的边缘处，通过对单个颗粒的持续追踪，可以确定颗粒在一定时间内的运动轨迹和速度，如图13 所示。据此，可对羽流的卷吸现象［3］开展进一步研究。

图13 羽流边缘颗粒运动检测（相邻帧时间Δt ＝0.8 ms）

基于上述颗粒速度的测试分析方法，在试验水槽内对200～250 mm 颗粒低浓度自由沉降过程进行分析，稳定段颗粒速度分布如图14 所示。可以发现，97%颗粒的速度范围为32～42 mm／s，均值为36 mm／s。作为对比，使用传统沉降法测得颗粒速度范围在30～45 mm／s，均值为39 mm／s；使用Stokes公式估算相应颗粒范围内的沉降速度区间为31～48 mm／s。可以认为，所提出的图像处理算法可以较好地分析颗粒的运动特征。

图14 颗粒自由沉降速度图像处理结果

需要说明的是，该速度图像处理方法不适用于颗粒浓度较大的羽流（质量浓度＞20 g／L），此类羽流颗粒重叠情况过多，难以使用算法对颗粒进行识别、分离、匹配等操作。

2.4 测量误差分析

试验中产生的误差主要来源于试验装置本身以及图像处理算法。

（1）羽流发生装置输送的混合物由于管道淤积等原因造成羽流初始浓度和流速存在小幅波动，影响羽流初始参数的精确确定。按照流程进行试验，先将混合物排放入标定池一段时间，再排入试验水槽中。采用质量分析方法，每隔5 s 记录一次标定池管口或试验水槽管口的质量浓度值，如图15 所示。结果表明，混合物经过管道输送需要一段时间（约10 s）才能达到稳定状态，同时标定池与试验水槽出口质量浓度存在一定偏差，重复测试结果表明该偏差值均在6%以内。因此，在试验开始前应尽量清除管内的残余颗粒，并且需要羽流排放10 s以上达到稳定状态后再进行数据采集工作。该发生装置通过调整输送泵的功率（0～150 W）控制流量，随着试验的进行混合物体积流量均呈现先升高再缓慢降低的趋势，如图16 所示。一般输送泵运行120 s 后，体积流量和颗粒质量浓度比初始稳定值降低7%以上，因此单次试验持续时间不宜过长。

图15 出口颗粒质量浓度变化

图16 不同泵功率下混合物体积流量曲线

（2）对于排放速度较大的工况，羽流触底后会扬起底部的部分颗粒，影响颗粒速度和浓度的图像分析结果。因此，每次试验前尽可能去除底部沉积的颗粒，并保证试验水槽中无其他悬浮杂质。

（3）浓度处理过程中，光照不均匀会显著影响浓度结果，因此试验过程中应尽量减少外界光源的影响，并采用均一无频闪强光源照射试验测试区域。

（4）速度处理过程中，算法对于较高浓度下重叠颗粒的处理存在误差，因此对于每一个工况应保证持续时间足够长，从而获得足够多的正确颗粒速度样本。

（5）在算法的二值化以及连通图提取步骤中，部分阈值需要根据颗粒大小以及光照强度确定。Muste等［16］在PTV 算法中也是根据经验选择二值化阈值。阈值选择不当会导致颗粒识别出现遗漏或将噪点识别为颗粒，因此试验参数及环境发生变化后，需要调整阈值以减少在该条件下试验产生的误差。

3 试验结果分析

3.1 试验工况

一般地，羽流数Np和弗劳德数Fr是影响颗粒羽流中颗粒相扩散率的重要参数［12］，分别表示为：

式中：Ws为颗粒的特征沉降速度；φ0为颗粒的体积分数；g′＝g（ρs－ρf）／ρf，其中ρs和ρf分别为颗粒和水的密度。由式（3）、（4）可知，u0同时控制Np和Fr这2 个无因次数。因此，u0对颗粒羽流特性影响更为显著。控制D不变，通过改变u0及c0来调整羽流初始状态，主要试验工况如表1 所示。

3.2 试验结果

Lai等［13］的试验结果表明，Np越大，羽流扩散率越大。针对工况1～6 开展试验，在其他条件不变的情况下通过改变初始流速改变Np，得到各工况下速度高斯半宽扩散率αp，如图17 所示。

图17 速度高斯半宽扩散率αp与羽流数Np关系

通过时均结果和瞬时结果对比可以发现，羽流存在局部湍流扩散。对于时均图像，颗粒相边界近似拟合为一条直线。瞬时图像与该拟合边界相比存在差异，如图18 所示。这是由羽流的卷吸现象导致的。当羽流初始流速较低时，瞬时图像与时均边界基本吻合；当羽流初始流速较大时，瞬时图像与时均边界存在明显的突出和凹陷，而且瞬时图像颗粒分布不均匀且不对称，羽流中心处的高浓度区域在边缘环境流体卷吸的作用下向一侧偏移。可以认为，初始流速高的羽流，卷吸或湍流扩散效应更加显著。

图18 瞬时图像与时均边界对比

针对工况1、3、5，将羽流中心z＝6D位置处的浓度作为标准质量浓度c0，则羽流中心相对浓度c／c0随高度z／D的变化规律如图19 所示。随着与羽流出口距离的增加，羽流中心相对浓度不断降低，不同工况下变化趋势基本相似。对于初始流速较大的羽流，卷吸效应更明显、颗粒相扩散率更大，因此颗粒的扩散浓度下降速度相对较快。

图19 不同初始流速下羽流中心相对浓度变化

羽流中颗粒在流体阻力的作用下动量逐渐减小，最终速度接近稳定自由沉降速度［12］。选择工况7、8、9，对低初始质量浓度（c0＝5 g／L）、不同初始流速（u0＝0.99～2.45 m／s）的羽流进行试验，分析距离羽流出口位置z／D＝17.5～22.5 范围内中心轴线颗粒平均速度um，测试分析区域为0.78 cm×0.20 cm的矩形区域，结果如图20 所示。在相同高度处，初始流速较高的羽流颗粒速度衰减率相对更高。

图20 不同初始流速下羽流中心颗粒相对速度变化

4 结语

提出了模拟多相羽流排放和扩散的试验装置与测试分析方法。基于时均灰度图像，采用边缘识别算法捕捉羽流边界和轮廓，获取颗粒相的扩散率；针对颗粒浓度，基于时均灰度图像并应用高斯分布假设，将二维平面图像中灰度值所表征的视场平均浓度转化为该位置羽流中心点浓度；针对颗粒速度，使用改进的PTV算法，提取非高浓度区域中颗粒速度，实现羽流边缘等较低浓度区域内颗粒速度和运动轨迹的连续捕捉。根据典型工况下试验，验证了该方法的准确性和可靠性。进一步研究了不同初始流速下羽流的扩散特征。结果表明，羽流初始流速越大，边缘卷吸效应越显著，颗粒扩散率越大，颗粒浓度与速度随高度的衰减越快。