刘明潇,Michele Guala ,孙东坡

(1. 华北水利水电大学水利学院,河南 郑州 450046；2. St. Antony Fall Laboratory,University of Minnesota,Minneapolis 55414,USA)

河流床面泥沙颗粒(以下简称粒子)的随机状态受粒子与近底水流相互作用以及床面边界效应的影响,是研究推移质输移规律的基础。从微观视角揭示粒子行为特征与近底水沙相干机制是泥沙科学继续发展的一个重要理论突破点[1]。借助先进的PIV/PTV(激光测速/粒子跟踪)技术开展对床面粒子间歇运动精细追踪观测,获取准确可靠的粒子运动特征值和等待参量,剖析“粒子运动- 休息循环间歇过程”的物理机制等,这些基础研究工作都离不开可靠的试验数据支撑。

在微观层面进行粒子尺度的推移质实验研究方面,国内学者胡春宏[2]较早使用高速摄影技术获取粒子运动轨迹,探索粒子间歇位移的随机性,分析了粒子步长、跃高的统计规律,指出粒径、水流强度、边界条件是床面粒子随机运动的重要影响因素；白玉川等[3]在试验研究基础上考虑床面粒子各阶段受力情况,分析了跃移参数的确定性关系及其统计分布规律。国外学者Hassan等[4]通过在天然河流中放入大量示踪剂,观测洪水过后示踪粒子位移的概率分布,进而研究推移质的扩散特征；Roseberry等[5]和Fathel等[6]借助水槽实验研究发现,粒子速度顺流向分量和横向分量的概率分布均表现为类指数函数特征,分布特征参数受水流强度影响,粒子速度的指数分布特征也已有初步试验验证。但目前对试验研究结果尚存一些分歧,有研究者认为粒子运动停时应服从指数分布形式,也有研究者认为粒子停时应服从“双隆起分布”,由快、慢速输移机制所引发的时间尺度分离效应与分布特征有关[7]。Fan等[8]研究指出粒子速度、单步距离、单步运动时间(简称“单步时长”)、等待时间这些变量的概率密度函数(PDF)并不全是细尾分布；在微观尺度上,床面粒子各自的单步位移均不相同；在宏观空间- 时间尺度上,床面群体粒子输移呈现出扩散性。Martin等[9]通过水槽实验发现具有长尾分布特征的粒子随机停时会导致群体泥沙的奇异扩散。限于试验观测技术与图像分析难度,目前难以给出足量的可靠试验数据支撑建立理论模型。因此,改进推移质试验观测技术,开展床面粒子状态特征的系统试验研究,获取更精准的试验观测数据,对厘清床面粒子随机状态与群体输移规律十分必要。另外在粒子随机运动特征研究领域,涉及粒子尺度的运动要素(速度,加速度,单步时长,等待时间等)概率分布及其对水流强度的依赖性[10- 11]方面还存有认识的差异。例如多数研究认为纵向粒子速度的概率密度函数呈指数分布,但存在这一特征的机理以及这种分布特征在低水流强度时的适用性还不清楚；尽管研究表明粒子单步运动距离(Ls)和单步时长(Ts)的PDF呈细尾分布,但分布函数形式以及与水流强度的关系尚存在争议[12]；虽然学界都认为床面待动粒子的状态对其后续运动很重要,但目前对等待状态参量的概率分布特征还研究甚少,粒子动- 静的判别阈值定义也并不明确[13]。目前明槽流床面泥沙运动的精细观测多采用基于图像识别的PTV技术,俯视图像采集受水面波动影响效果欠佳；采用压波板改善也避免不了水珠、气泡等干扰,还是难以获得床面粒子高清图像。因此,有必要对这些问题深入开展研究。

本文拟利用开发的水下摄影技术(UP),借助水槽系列试验,精细观测床面泥沙运动状态;在数据滤波处理、误差分析基础上,利用试验数据进行床面粒子运动状态的分析,探寻粒子运动要素遵循的统计学规律。

1 水下摄像试验观测系统与工作流程

1.1 试验观测系统

(1) 试验水槽。床沙运动观测试验在可调坡精密水槽中进行。水槽为矩形断面,宽0.3 m,高0.5 m,长16 m；水槽前端采用多层卧管梳理稳定来流,如图1所示,其中,H为水深,δ为测点距底高度。水槽顺序设置上游粗砾过渡段L1,动床铺沙段L2,近底泥沙信息采集区L3、动床近底流速采集区L4,下游过渡段L5。水槽采用循环供水方式,由变频器调控试验水流流量。为观测推移质运动状态,试验床面为厚8～10 cm的泥沙,级配相对均匀,见图2。

图1 试验水槽及观测设施示意Fig.1 Experimental Flume and measure instruments

(2) 主要量测仪器。流量采用E- mag电磁流量计测量,水深测量采用超声波水位仪,流速测量使用ADV多普勒测速仪,分上下2个区域采集,仪器布置见图1。采用UP技术,对床面泥沙运动进行图像采集,装在自由旋转卡上的水下相机设置在近底泥沙信息采集区上方,推移质输移量采用电子天平动态记录。

1.2 试验水沙条件及图像采集

有研究认为[14- 15],适当控制床面粒子运动强度,有助于采集清晰的粒子运动图像,减小识别难度,能够获取粒子随机运动的较完整信息。故本试验按满足泥沙初动条件控制水流强度:近底流速u=0.3～0.4 m/s,H=0.15～0.20 m,弗劳德数Fr=0.25～0.32,相应床面切应力τ0=0.3～0.7 Pa,其中摩阻流速u*利用实测脉动流速确定,并考虑测速位置及床面剪切力分布特性影响[16]。

为避免粒径差异对粒子起动及等待的干扰,采用粗细相对均匀的床沙,密度ρs=2 600 kg/m3,中值粒径D50=1.15 mm,泥沙级配见图2。采用表1所示水流控制条件开展5种工况的床沙运动试验,粒子视频采集密度不低于30个/组次；利用UP技术采集的水下粒子状态见图3。

图2 床沙粒径级配Fig.2 Sediment gradation used in experiments

表1 试验控制水流条件

图3 试验床沙的水下图像Fig.3 Underwater image of the test bed sands

1.3 工作流程

1.3.1 基于UP技术的床面粒子运动图像采集

采用分辨率为640×480像素的超薄小型防水相机,将相机固定在360°全方位可调的相机卡架上；在水槽试验段布置LED补偿光源和卤素灯晕光装置,增加光漫射效果。图像采集频率采用120帧/s,每2帧图像间隔0.008 3 s；能够完整采集粒子动、静2种状态及交替转换过程。

为不影响近底流场同时又可获取高清图像,分析Musa等[17]、王浩等[18]的水槽试验研究成果,结合Liu等[19]在以往水槽试验中对粒子图像采集的分析,认为在水深较大、流速较低时,相机镜头距床面高度(hs)按hs/D50>90且H/hs<1.7控制为宜。本试验中一般hs>10 cm,相机对近底流场干扰的影响在可忽略的范围内；水下拍摄粒子图像没有自由面水纹噪点,虽然粒子很细但图像清晰度依然很高,见图3。采用UP方式可实现对无压明流床面泥沙运动的精细观测,同时又避免水面波动对采集图像质量的影响。

1.3.2 基于动态阈值的粒子运动跟踪识别流程

粒子追踪识别基本过程包括对采集视频图像进行读取、解析及运算等,具体工作流程为:导入床面粒子运动图像文件；利用Matlab工具,识别目标粒子及位置信息；粒子命名,多目标动态锁定、追踪；计算粒子运动状态参量；筛选提取有效长序列粒子踪迹数据；输出粒子Lagrange运动过程线；输出粒子各状态参量。为提高粒子识别精度,采用动态阈值和优化的Baek- Lee松弛算法,开发粒子追踪专用识别程序；利用数据滤波技术,建立有效粒子的时空信息数据集。

2 试验成果及分析

2.1 床面粒子状态分析

采用前述试验设备及控制条件进行多组次试验观测,基于UP/PTV技术获取了大量清晰的粒子运动图像信息；利用采集视频中连续帧上移动的粒子质心像素坐标(xp,yp),理论上可以获取粒子运动轨迹及相应的运动学指标,但实际上还必须通过数据筛选,才能得到可供分析的床面粒子运动时空信息及各状态特征值,包括粒子运动的空间轨迹、运动要素、等待时间及单步时间等。

2.1.1 床面粒子的活跃性与速度阈值的设定

在床面粒子随机运动过程中,运动与休息(等待再起动)是2种基本存在状态,但很难准确界定。很多情况下粒子处于动与不动之间的非稳态,如何科学定义和度量床面粒子的这2种状态非常重要,这直接影响粒子各运动要素的构成及随机分布特性。为精准反映床面粒子行为特点、合理筛选分类数据,根据获取的床面粒子图像特征,需要区分活性粒子与惰性粒子。活性粒子在床面表现为活跃运动状态,而惰性粒子在床面处于上静止或仅有小幅度震颤。

利用跟踪获取到床面上单颗泥沙的运动轨迹,这里给出床面上一个典型粒子运动和静止状态交错出现的图像轨迹与时空序列图,如图4(a)所示。从粒子运动轨迹的空间重构中,可以看到粒子间歇运动(步进、跳跃)和休息(等待)的时间特征,如图4(b),进而获取粒子单步运动时间及单步步长,up为粒子速度。

图4 运动- 休息状态循环交替的粒子间歇运动Fig.4 Particle trajectory in intermittent motion with the circulate transitions between motion and rest states

为区分粒子运动与休息2种基本存在状态,需要区分粒子的当地摇摆与运动的差异。通过水槽试验采集的水下摄影图像分析表明,不少床面粒子处于在原地摇摆而没有明显净位移的状态；为了剔除这些干扰因素,应设置界定床面粒子发生有效运动的临界指标,即粒子运动阈值(uc)(临界速度)。根据本试验数据分析,可以选择单次采样时间内纵向净位移达到0.25D50的粒子速度为运动阈值,这与Lajeunesse[20]等采用的速度阈值基本一致,一般约为实测运动粒子质心不确定所产生误差的2倍。

对于处于不断振荡摇摆状态的粒子,虽然它们似乎具有运动特征,但只在自身位置处震颤,净位移很小；故将这种(惰性粒子)状态定义为非运动的活跃等待状态。区分粒子处于运动或等待状态时,要考虑床面颗粒晃动状态,还须关注颗粒起动临界特征。根据大量试验数据分析认为,判别粒子运动应满足:①沿流向的局部累积净位移必须超过D50/4才被认为是运动；②仅当粒子速度高于设定的粒子速度阈值。

滤除惰性粒子的干扰数据,在视频连续帧上采集达到起动识别标准的床面活性粒子像素的平面坐标；识别每个活性粒子在各帧间的运动轨迹,获取粒子的帧间纵横向位移和相应流场中纵横向分速度与加速度。

2.1.2 粒子间歇运动状态特征及精度评估

根据试验统计分析[19],以粒子运动轨迹10D50为界,当床面活性粒子运动轨迹累积超过此界以后,通常粒子都会显示出运动与等待交替出现的行为过程特征。在对运动粒子的图像识别中,粒子图像在每个后续视场中的移动都会发生独具特色的变化,从而被识别与跟踪,见图5。

图5 颗粒在1个有效样本中的运动状态及定位特征Fig.5 Particle states in one reliable sample of particle trajectory and particle locating characteristics

图5(a)为追踪到的1颗粒子的运动轨迹,图中清晰显示了间歇运动粒子的床面位置(x,y)变化过程；图中标识的5处“停”表示在这些位置粒子进行了短暂休息,随后又继续运动。2个休息(停)之间为1个单步运动,即1个步长,其中的每1小步为1次步进。粒子运动跟踪精度取决于对粒子定位的准确性,这可以通过对沿运动轨迹采集到的等效粒径(当量粒径)标准偏差(D′)来估计,图5(b)显示了运动粒子当量粒径沿流向的变化情况。经估算分析,D′约为0.11Dave,Dave为沿流向运动粒子检测面的平均直径(本试验约为1.25 mm)；统计分析样本误差约为0.048Dave(0.06 mm),实验观测精度满足要求。为了减小定位误判,试验中要求单颗泥沙的直径在采集图像中不小于3个像素,试验泥沙粒径不小于0.25 mm。

2.1.3 粒子运动的拉格朗日时间变化特征

通过观察沿运动轨迹的活性粒子数量、粒子速度,可以分析粒子在拉格朗日时空体系中的运动特征。从起动开始分析活性粒子的单步运动,提取单步运动时间点上(单步运动中每1帧时间间隔内)的平均粒子速度,探寻粒子速度(up)与拉格朗日特征时间(τl)的变化规律。图6表明了运动粒子数量(Np)与粒子运动要素沿τl的变化。在图中只选取粒子处于运动状态的数据,不牵涉止动(等待)阶段的影响。图6(a)中up变化曲线体现了床面粒子典型运动特征:在τl<0.3s时段内,粒子运动速度先是急剧增加然后表现为缓慢衰减,随着拉格朗日时间增长,运动粒子数量显著减少；在τl>0.3 s后只有很少粒子还在维持运动状态,同时近底紊流的强脉动性导致残余粒子的运动速度呈现随机波动。粒子纵向加速度(as)在τl<0.3 s的变化特征与up相近,只是衰减更快,如图6(b)所示；在τl>0.3 s后(右侧)少量粒子速度急剧变化引起加速度更明显的随机波动。这样的粒子运动场景本质上是近底紊动猝发机制的产物:在拉格朗日体系中,开启粒子单步运动旅程的粒子位移,最初源于粒子- 粒子相互碰撞[21],或紊流强扫荡[22]的突然“激发”,瞬时紊流脉冲力使粒子突然加速进入床面群体运动层中；但在持续的床面摩阻力作用下,在1个步长内粒子速度逐渐减慢,同时加速度在由正变负过程中随机波动；这和粒子与床面的碰触摩擦有关,同时也受紊流强扫荡后的喷射影响,粒子加速度的阵发性特征也是对近底紊流脉动的响应。

图6 颗粒数量、纵向速度和加速度随Lagrange时间的变化Fig.6 Development of particles quantity,velocities and acceleration as a function of the Lagrange time

2.2 床面粒子运动特性的统计学规律分析

2.2.1 粒子纵向运动速度的概率分布

利用系列水槽试验结果,分析了不同水流条件下床面粒子运动速度的概率分布特征。分析表明,5种试验工况(见表1)下所得到的粒子速度的概率分布具有相似特征,主要体现在概率密度曲线都具有指数分布的尾部变化特征。图7仅给出试验工况D和工况E条件下粒子纵向速度的概率密度分布曲线,基于统计学理论分析,可用数学上的细尾Gamma函数来表示其概率密度分布特征:

(1)

式中:Γ(a)为Gamma函数;a和b是决定概率分布特性的2个参数。

图7中对应2种工况的特性参数(a,b)分别为:工况D(2.827,0.030),工况E(2.618,0.033)。如果考虑粒子原地颤动、前后摇摆对起动判别的影响,即使较低的粒子速度也遵循指数分布,这与Roseberry等近年的试验结果基本一致[4]。

图7 粒子运动速度纵向分量概率密度分布曲线Fig.7 Probability distribution of the stream wise particle velocity component

在水流强度逐渐增大的过程中,床面剪切速度(摩阻流速)也相应增加,从0.067 m/s逐渐增大到0.086 m/s(工况A—工况E)。对比各种摩阻流速条件下粒子速度概率分布密度与速度均值的变化,发现水流摩阻流速是粒子速度的密切相关因子。试验结果分析表明,粒子纵向运动速度与摩阻流速间存在基本稳定的比率关系:up/u*=3.4～3.5。这2个物理量是影响推移质输移率的关键因素。

2.2.2 粒子纵向加速度的概率分布

利用采集的粒子轨迹,分析了纵向加速度的统计特征,图8给出了工况B和工况C中as的PDF曲线。在5种试验工况中,粒子加速度的概率密度分布均呈现拉普拉斯分布特征,其函数表达如式(2)所示:

图8 颗粒加速度纵向分量的概率分布曲线Fig.8 Probability distribution of the particle acceleration component

(2)

式中:μ为位置参数,可以用中值来表征,试验数据的中值在0.23～0.33之间,已接近于预期值0,这对于变化范围很大(±40)的粒子加速度而言,样本数据体现了很好的统计特征；b0为尺度参数,反映随机量的主要分布范围,可用as与μ的标准差来表征。

利用数据较多的4个试验组次(工况B—工况E),得到加速度的统计特征值,如表2所示；其中asrms为加速度的标准差,asrmsTs/u*为量纲一加速度；这里使用摩阻流速与单步时长的水沙组合因子u*/Ts对加速度进行量纲一化处理。表2的统计数据分析表明,决定PDF曲线分布宽度的加速度标准差和量纲一加速度均呈现随水流强度增加(工况B—工况E)而增加的变化趋势。

表2 粒子加速度统计特征值

2.2.3 粒子单步运动时间的概率分布

单步运动时间(单步时长)既是体现单个粒子运动状态又是可以表征群体粒子输移强度的重要指标,在5种工况中选取B和C,图9给出了Ts的PDF分布特征,曲线一致呈现出指数分布的细尾特征,可用式(3)表示:

图9 颗粒单步运动时间的概率分布曲线Fig.9 Probability distribution of the particle step time

(3)

式中:参数λ为变量Ts的均值。由试验数据分析得到5种水流条件对应的Ts平均值分别为0.11 s、0.12 s、0.12 s、0.13 s、0.13 s,表明在试验输沙强度范围内,单步时长随水流强度增加略有增大；这意味着参与运动的粒子总数增加和对应粒子等待时间减少。

粒子单步步长和单步时长是相关的,试验数据分析表明步长的PDF也呈指数细尾分布,且二者的联合分布有较好的幂指数变化特征,如图10所示。分析4种水流条件(工况B—工况E)的数据拟合趋势线(约有22%的偏差),各趋势线斜率(即幂律指数)随水流强度增大呈减小趋势,变化范围在1.25～1.30之间。与Roseberry等研究发现的指数1.67相比[4],本试验幂指数略小。两者的差异可能与前者的试验Fr数较低有关,粒子步长对水流强度的响应使幂指数偏大；另外床沙异质性也会对粒子运动特性产生影响。本试验床沙粒径80%为0.7～1.3 mm,前者试验床沙多为0.45～0.55 mm。从图可以看出作者试验数据得到的联合分布点据丰富、集中且分布趋势明显,与Roseberry曲线(图10中SR趋势线)基本一致,说明试验数据具有较好的质量与精度。

图10 颗粒单步步长与单步运动时间的联合分布Fig.10 Joint scatter plot of particle step time and step length

3 结 论

为精细研究床面粒子状态特征,研发了一种基于水下摄影和粒子跟踪技术的床沙运动试验观测方法,探讨了水下图像采集的关键技术与数据处理方法,提出了活性粒子与考虑震颤效应的起动阈值,可以保证粒子图像识别精度；通过5种水流强度多组次动床水槽试验表明,采集试验数据质量比较高,能满足基于拉格朗日体系的粒子状态特性研究需求。研究表明:

(1) 在1个单步步长内,粒子速度先是急剧增加然后缓慢衰减；同时,随着拉格朗日时间(τl)增长,粒子运动数量也在显著衰减,在τl>0.3 s后只有很少量粒子的运动还在持续;同时粒子速度与加速度都呈现强烈的随机波动,这种特性与粒子间的碰撞和近底紊流的扫荡现象有关。

(2) 床面粒子运动速度的概率密度函数(PDF)曲线表现出很好的细尾Gamma函数分布特性,粒子加速度的PDF曲线具有拉普拉斯函数分布特征；粒子的当地震颤摇摆效应会影响其速度概率分布特性；粒子运动速度与摩阻流速之间关系密切,比率基本在3.4～3.5之间。

(3) 粒子的单步运动时间与单步步长的概率密度分布曲线都一致地呈现出指数分布的细尾特征,两者的联合分布则呈幂函数规律变化。试验拟合曲线的指数一般为1.25～1.30,研究表明幂律指数会随推移质输沙率增加而逐渐降低,床面粒子的异质性也对幂律指数的大小产生影响。