CN104297516B - 一种流体表面二维流速场测量方法 - Google Patents
一种流体表面二维流速场测量方法 Download PDFInfo
- Publication number
- CN104297516B CN104297516B CN201410617463.2A CN201410617463A CN104297516B CN 104297516 B CN104297516 B CN 104297516B CN 201410617463 A CN201410617463 A CN 201410617463A CN 104297516 B CN104297516 B CN 104297516B
- Authority
- CN
- China
- Prior art keywords
- image
- fluid
- visual
- velocity field
- rectangle
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Landscapes
- Indicating Or Recording The Presence, Absence, Or Direction Of Movement (AREA)
- Image Analysis (AREA)
Abstract
本发明公开了一种流体表面二维流速场测量方法。针对现有技术中测量运动流体表面二维流速场的方法存在不适用于野外流体表面场测量、任意角度拍摄影像测量、表面矢量流速场测量的缺陷,本发明提供了一种流体表面二维流速场测量方法。本方法利用光流算法实现在较大程度减少计算量的情况下获得像素级别的流动表面二维流速场图。在此基础上,本发明方法利用伪矩形图像透视投影变换方法,实现对从迎流体运动方向上任意角度捕捉流体运动影像进行图像分析与表面流速场测量。本发明还提供了一种运动流体图像透视投影变换计算方法。本发明方法原理可靠,计算过程科学简便,实施仪器简单,安装方便,测算结果精度高。尤其能够应用于泥石流表面二维流速测量。
Description
技术领域
本发明涉及一种流体运动特征测量方法,特别是涉及一种流体表面二维流速场的非接触式测量方法,属于水利工程、测量技术领域。
背景技术
流体的表面流速场是流体表面的速度分布特征,是描述流体运动、运行特征的重要物理量。流体表面二维流速场具有非定常特性,是一种复杂而又普遍地反映自然本质的物理过程。流体力学中许多疑难问题的突破都取决于流速场测试技术的进展。
泥石流是一种特殊流体。在其相关研究试验中,泥石流流速的分布和计算是泥石流运动力学研究的核心问题之一,泥石流流速场也是各类泥石流防治工程设计中的核心参数。作为一种复杂的多相非牛顿体,无论是在野外原型观测条件下还是室内实验条件下,泥石流表面二维流速场的测量,特别是精度测量,都存在较多困难。再加之泥石流运动时产生的强破坏力,使得泥石流表面二维流速场的测量愈加困难。目前,泥石流表面二维流速场研究主要采用数值模拟进行,尚无能对泥石流表面二维流速场进行精确测量的方法。
现有技术中用于测量流体表面流速的方法主要有三种,分别是粒子图像测速技术(PIV,particle image velocimetry)、基于网格对象跟踪的测量技术、基于多普勒效应的电磁波测量技术。
粒子图像测速技术是七十年代末发展起来的一种瞬态、多点、无接触式的激光流体力学测速方法。近几十年来得到了不断完善与发展,并且成为主流的流体表面流速场测量方法。该技术缺陷主要在于两点:一是必须在流体中掺入密度与流体相当并具有很好的跟随性的示踪粒子,通过追踪示踪粒子来测得流速场。这在许多实际场景中是不可能实现的,比如在沟道中实际运动的泥石流。二是PIV设备多限于实验室内部,流体规模很小,成像设备视角狭窄,无法胜任野外流体表面场的测量。三是PIV实验无法满足泥石流试验的需求,泥石流流体为浑浊的流体,泥浆包裹着示踪粒子,无法通过图像设备来追踪示踪粒子的运动轨迹。
基于网格对象跟踪的实时流体表面速度测量方法是由网络摄像头、路由器和计算机组成的监测系统中,采用图像输入和网格化,网格对像跟踪,实现实时流体速度的监测。该技术缺陷主要在于两点:一是跟踪对象是一个个的网格对象,无法对每个像素进行跟踪,从而无法得到像素级的流速场图。二是该方法需要相机从运行流体上方正对着下方才能得到比较精确的速度,如果相机从别的角度拍摄流体就无法实施其技术。而在实际场景中很难达到从流体正上方往下拍摄这样苛刻的条件。
基于多普勒效应的电磁波测量方法主要利用多普勒效应(Doppler Effect),即移动物体对所接收的电磁波有频移效应的原理,通过借由频率的改变数值计算出目标与雷达的相对速度,进而得出被测物体的运动速度。该方法可用于测量流体表面流速值,但同样存在两点技术缺陷:一是电磁波发射器或雷达发射器价格昂贵,不可能在野外条件下大量安装使用,因此方法不适用于推广使用。二是该方法仅能测量出一个点或者多个点的表面速度值,不能测量出流体表面的矢量流速场。
发明内容
本发明的目的就是针对现有技术的不足,提供一种测量流体表面二维流速场的方法。该方法通过对从迎流体运动方向上的任意角度捕捉的流体运动视频进行图像分析,可以实现对运动中流体表面流速场的精确测量。
为实现上述目的,本发明首先提供一种利用正投影拍摄图像测量运动流体表面流速场的方法,其技术方案如下:
一种流体表面二维流速场测量方法,其特征在于:依以下步骤实施:
步骤S1、流体运动影像获取及初始化处理
步骤S11、现场设备安装
在运动流体测量现场安装数字拍摄设备,使拍摄镜头自正上方垂直向下拍摄获得运动流体正投影方向上的流体运动影像;
步骤S12、获取流体运动影像
获取流体运动影像,将影像输入数字处理设备;
步骤S13、图像栅格化处理
拾取流体运动影像上任意静止图像,设置像素参数对图像进行栅格化处理,得到图像栅格点阵;
步骤S14、设置图像处理参数
设置图像处理速率、图像实际处理分辨率、像素参数与实际位移换算比例;
步骤S2、计算流体表面在视觉坐标系xyz中的二维流速场
步骤S21、建立流体运动影像三维空间视觉坐标系xyz
在流体运动影像上建立三维空间视觉坐标系xyz获得视觉平面xoy;
步骤S22、计算流体运动图像稠密光流
拾取流体运动影像播放中任意两帧图像,两帧图像间隔时间Δt,计算两帧图像上所有对应像素栅格点的稠密光流其中(Δi,Δj)是测量区内点(i,j)在时间Δt内在视觉坐标系xyz中经过的位移;
步骤S23、计算图像点阵在视觉坐标系xyz中的像素流速场
计算所有点阵在图像视觉坐标系xyz中的移动速度,得到式1所示视觉平面像素流速场:
式中,—点阵在视觉平面像素流速场;
步骤S24、计算点阵在视觉平面xoy中的二维流速场
根据像素参数与实际位移换算比例,将步骤S23所得点阵在视觉平面中的像素流速场换算得到点阵在视觉平面xoy中的视觉投影变换流速场;所得流体表面在视觉坐标系xyz中的视觉投影变换流速场即为运动流体表面二维流速场。
上述测量方法中,像素参数设置为1pix~100pix为优,图像处理速率设置为10ms/帧~100ms/帧为优。帧间隔参数设置根据成像设备(一般为相机为摄像机)的帧率和流体运动速度加以确定,成像设备帧率越高或流体运动速度越低,则帧间隔设置越高。帧间隔一般设置区间在1~10,且以1~5为优。如果数字处理设备(通常为计算机)性能较好应该尽可能按图像原分辨率处理视频。如果处理实时处理视频但是数字处理设备性能不足出现丢帧,则适当降低分辨率。
上述测量方法的基本技术原理在于:
原理一:利用数字摄像设备拍摄一组关于流体运动的连续图像序列,即流体运动影像。由于运动流体表面存在对比度差异的颜色显著区域,因此可以观察到流体流动。将连续流体充分分段后,流体的各段表面近似平面,此平面内流体表面速度可以计算。计算所有连续平面中的流体表面速度,则能计算出流体表面二维速度分布。
原理二:由于用于分析的流体运动图像是自运动流体正上方垂直向下拍摄而获得的正投影方向上的影像,所以计算获得的流体表面在视觉坐标系xyz中的视觉投影流速场即为测量得到的运动流体表面二维流速场。
原理三:对于单目相机拍摄得到的流体运行影像,如果流体在流动,那么某测量点(i,j)在连续的帧与帧图像之间的对应位置存在微小位移,即该测量点(i,j)成为移动的点。根据稠密光流法(DOF)理论,移动的点在帧时间间距内灰度保持不变。因此可以凭借其灰度值追踪该测量点在下一帧画面中(即下一时刻)的坐标位置。因此,在流体运动影像上建立三维空间视觉坐标系xyz,再计算图像中所有像素点的光流就可得到稠密光流。其中,和分别是任一测量点在x方向和y方向的位移变化率。
此处,对于某测量点(icurrent,jcurrent)定义的最小化残差ε,有:
由于测量点(i,j)的灰度保持不变,因而有:
ft(i,j)=ft+Δt(i+Δi,j+Δj) 式3
其中,由于图像在前后帧之间变化并不明显,因而可近似认为:
结合式3、式4有:
并最终得到:
对运动流体影像分析,可以得到大量观测值数据矩阵对其求矩阵伪逆得到从而测算出图像上所有点在x、y方向上的移动速度,得到所有像素点的速度指标并显示为像素级别的流速场图(即流体在视觉坐标系xyz中的像素流速场)。进一步利用步骤S1中设定的像素参数与实际位移换算比例,可以得到流体在图像视觉坐标系xyz中的实际流速场。
基于上原理,上述测量方法步骤S22中计算流体运动图像稠密光流可利用图像空间灰度梯度来有效地匹配前后图像中的对应点,具体可采用模版匹配算法实现。
本发明还提供上述测量方法的优化技术方案,具体是透视投影变换的计算的技术方案。
上述测量方法是将摄像机从平面x′oy′正上方垂直于该平面的角度往下拍摄,则视觉平面xoy与流体表面平面x′oy′重合,由此得到流体在z′方向速度为0的影像。在此条件下计算所得运动流体表面在视觉平面的流速场即为流体表面二维流速场。
然而在实际场景中,受限于现场硬件条件,很难从自然流体正上方向下拍摄获得正投影特定方向的流体运动影像,往往只能获得从其它任意一个角度(通常是迎流体运动方向上的任意一个角度)拍摄至流体运动影像。此种情况下,无论从哪一个角度拍摄,摄像设备拍摄时的成像平面都存在透视变形,得到的流体运动影像也都是经过透视投影的,则步骤S2中计算所得的视觉平面中的投影变换流速场并不是流体表面实际的二维流速场。为解决这一技术问题,使本发明测量方法具有更广的适用性与推广价值,本发明进一步提供一种可以从迎流体运动方向上的任意角度拍摄运动流体影像并通过图像透视投影变换计算运动流体表面二维流速场的方法。
为实现上述目的,本发明首先提供一种运动流体图像透视投影变换计算方法,该方法利用流体运动测量现场布置的矩形参照物实现图像透视投影变换,其技术方案如下:
运动流体图像透视投影变换计算方法:依如下步骤实施:
步骤S1a、设置矩形参照物A″B″C″D″
在运动流体测量现场安装数字拍摄设备,从迎流体运动方向上的任意角度获取流体运动影像;在测量流体附近能够被摄像设备完全拍摄且不影响流体运动的位置设置矩形参照物A″B″C″D″,确定矩形参照物A″B″C″D″边长A″B″、C″D″长度;
步骤S2a、确定伪矩形A′B′C′D′坐标
获取流体运动影像,确保影像中有矩形参照物A″B″C″D″的完整图像;拾取流体运动影像上任意静止图像,矩形参照物A″B″C″D″在静止图像上对应伪矩形A′B′C′D′;在流体运动图像上建立三维空间视觉坐标系xyz,计算机确定伪矩形A′B′C′D′四角点在视觉平面xoy中的图像坐标A′(i1′,j1′)、B′(i2′,j2′)、C′(i3′,j3′)、D′(i4′,j4′);
步骤S3a、设置投影矩形ABCD
在视觉平面xoy中建立投影矩形ABCD,投影矩形ABCD四角点在视觉坐标系xyz中坐标分别是A(i1,j1)、B(i2,j2)、C(i3,j3)、D(i4,j4);所述A(i1,j1)与A′(i1′,j1′)重合、B(i2,j2)与B′(i2′,j2′)重合,确定边长AB;根据式7确定边长BC及C(i3,j3)、D(i4,j4)坐标;
步骤S4a、计算透视投影变换矩阵PPM
根据A(i1,j1)、B(i2,j2)、C(i3,j3)、D(i4,j4)与A′(i1′,j1′)、B′(i2′,j2′)、C′(i3′,j3′)、D′(i4′,j4′)的映射关系建立投影矩形ABCD与伪矩形A′B′C′D′间的透视投影变换矩阵PPM:
步骤S5a、计算比例尺k
依式9计算比例尺k
步骤S6a、计算静止图像中点的实际坐标
将静止图像上任意点的图像坐标经透视投影变换矩阵PPM变换得到该点在视觉平面xoy中的投影坐标,再将投影坐标与比例尺k相乘得到静止图像上任意点在流体表面真实平面x′oy′中的实际坐标。
上述运动流体图像透视投影变换计算方法的技术原理在于:
几乎所有的摄像设备成像都是将对像经过透镜投影到相机成像平面,这一过程就是一个投射投影过程,导致图像原来的集合形状产生了变化。由于采用单目摄像设备采集的视频只有平面数据无深度数据,即只有视觉坐标系xyz中dx、dy速度,无dz速度。根据上文原理一中分析,当流体充分分段后,其各段表面近似平面。进一步地,如果在视觉坐标系xyz中已知一个点位于某平面ax+by+cz+d=0上(系数a、b、c、d已知),那么当x、y已知,则可求出z。由此则,在流体各段视觉平面xoy上,根据dx与dy可确定dz速度。进一步地,视觉坐标系xyz与流体表面真实坐标系x′y′z′间存在透视投影变换的换算关系,其转换关系用透视投影变换矩阵PPM(perspective projection matrix)表示。利用投影变换矩阵PPM可以把视觉平面xoy上的点坐标转换为其在流体表面平面x′oy′上的坐标。于是用模型描述这种变换:假设(i,j)是某点在视觉平面xoy上的位置坐标,(i′,j′)是该点在流体表面真实平面x′oy′上的物理坐标位置,则投影变换可以描述为:
对于上述变换,只需要4对观测点就可以完全求出系数矩阵:
因为i·w=p00·i'+p01·j'+p02以及w=p20·i′+p21·j′+1 式11、式12
所以有i=p00·i'+p01j'+p02-p20·i·i'-p21·i·j' 式13
可以得到等式:j=p10·i′+p11·j′+p12-p20·j·i′-p21·j·j′ 式14
然后通过4对观测点得到等式:
对式15右边的系数矩阵求逆,然后左乘等式两遍可求得系数,得到投影变换矩阵PPM即式8,进一步可完成坐标转换。
基于上述原理,本发明运动流体图像透视投影变换计算方法首先在测量现场布置矩形参照物A″B″C″D″,参照物实际边长A″B″、C″D″可以确定。矩形参照物A″B″C″D″在拍摄图像上变形成为伪矩形A′B′C′D′,其四角点图像坐标A′(i1′,j1′)、B′(i2′,j2′)、C′(i3′,j3′)、D′(i4′,j4′)可以通过计算机确定(计算机自动读取或手动标定)。视觉平面xoy中建立投影矩形ABCD,其四个角点坐标分别是A(i1,j1)、B(i2,j2)、C(i3,j3)、D(i4,j4)。设置投影矩形ABCD与伪矩形A′B′C′D′底边重合,即A(i1,j1)与A′(i1′,j1′)重合、B(i2,j2)与B′(i2′,j2′)重合,由此可以确定投影矩形ABCD底边AB边长。设置投影矩形ABCD与矩形参照物A″B″C″D″有相同的长宽比例,因此利用式7可以确定边长BC及C(i3,j3)、D(i4,j4)坐标。由此确定了伪矩形A′B′C′D′与投影矩形ABCD上4对观点(即四个角点)的对应坐标,因而可以建立透视投影变换矩阵PPM。
假设一个伪矩形A′B′C′D′中的点在前后两帧图像的视觉坐标系xyz中的坐标为(x0′,y0′)与(x1′,y1′),则其投影到投影矩形ABCD中的坐标为(x0,y0)与(x1,y1)。由于投影矩形ABCD的长宽比与矩形参照物A″B″C″D″相同,且投影矩形ABCD内的点坐标是根据透视投影变换矩阵PPM计算得出,因此投影矩形ABCD内的几何形状与真实的物体几何形状是相同的。也即可以把投影矩形ABCD看作是拍摄制备从矩形参照物A″B″C″D″正上方垂直往下拍摄得到的图像。在此基础上,只需要将投影矩形ABCD内点坐标乘以比例尺k就可以得到这些点在真实矩形A″B″C″D″内的坐标。由此,也能最终得到静止图像上任意点在流体表面真实平面x′oy′中的实际坐标。进一步地,利用PPM可以直接将伪矩形A′B′C′D′内的点的坐标(图像失真条件下的坐标值)转换到投影矩形ABCD。
本发明测量方法,如果流体各段不在一个平面上,比如快流下坡度较陡区域后又缓缓流过较平缓区域,下坡区域和平缓区域不在一个平面上,此时应该分成2段处理,则应该针对不同平面分别划分透视投影变换区域,建立独立的透视投影变换PPM进行后期计算。
本发明方法是一种非接触式精确测量运动流体表面二维流速场的方法。可以对实现对普通流体的测量,以及以泥石流为代表的复杂多相流体的表面二维流速场的测量。
与现有技术相比,本发明的有益效果是:(1)本发明提供了一种流体表面二维流速场测量方法,以正投影拍摄图像为基础,利用光流法测量运动流体表面流速场,可以在较大程度减少计算量的情况下获得像素级别的流动表面二维流速场图。(2)本发明进一步提供的流体表面二维流速场测量方法可以通过从迎流体运动方向上的任意角度捕捉的流体运动影像进行图像分析,通过伪矩形图像透视投影变换实现对运动中流体表面流速场的精确测量。(3)本发明还提供了一种运动流体图像透视投影变换计算方法。
附图说明
图1是测量自然条件下弯曲运动流体表面流速场装置示意图。
图2是测量模拟条件下直线运动流体表面流速场装置示意图。
图3是实施例一测量方法流程图(示循环计算)。
图4是栅格化处理及测量区点阵位置示意图。
图5视觉坐标系xyz建立示意图(箭头表示流体运动方向)。
图6是平面演示点二维跟踪示意图(图中,(i,j)是t时刻像素点位置及坐标,(i+^i,j+^j)是(t+^t)时刻像素点位置及坐标)。
图7是样本数据所在测量区像素栅格点图像速度计算结果。
图8是流体表面在视觉平面上像素流速场效果图(第111帧、第4.48s时刻图)。
图9真实坐标系x′y′z′建立示意图(箭头表示流体运动方向)。
图10是透视投影变换关系示意图。
图11是流体伪彩色表面二维流速场效果图。
附图中的数字标记分别是:
1实验水槽 2矩形参照物 3数字摄像设备
具体实施方式
下面结合附图,对本发明的优选实施例作进一步的描述。
实施例一
如图1~图11所示,用本发明方法测量运动流体表面二维流速场。
1、设备安装与连接
图1是测量自然条件下弯曲运动流体表面流速场装置示意图,图2是测量模拟条件下直线运动流体表面流速场装置示意图。
本实施例在实验室条件下进行,采用实验水槽模拟自然状态下运动流体,采用图2所示装置测量运动流体表面流速场。实验针对特殊多相流体进行。实验水槽中流体采用水、土、石块混合成而,模拟泥石流组成。本实施例中采用本发明提供的图像透视投影变换计算方法实现测量中的透视投影变换。
如图2所示,设备包含实验水槽1、矩形参照物2、数字摄像设备3、计算机、必要的数据连接线。实验水槽1是长矩形槽结构,坡度15°,底部与横截面均呈矩形。实验水槽长5m,横截面宽1m。数字摄像设备3安装在实验水槽1上方,其镜头对准实验水槽1且迎流体运动方向,保证镜头稳定不发生明显晃动。数字摄像设备3采集实验水槽1中运动流体影像,并通过数据连接线传输至计算机。计算机接收运动流体影像视频,对其进行实时的数字处理与运算,同时依需要保存输出分析结果。在实验水槽1与数字摄像设备3相对的一端出口处布置矩形参照物2(即矩形参照物ABCD),矩形参照物ABCD能够被摄像设备镜头完全拍摄到。矩形参照物AB边长1m,BC边长4m。
2、测量流体表面二维流速场
图3是本实施例测量方法流程图(示循环计算)。
步骤S1、流体运动影像初始化处理
计算机获取数字摄像设备拍摄的流体运动影像。
步骤S12、图像栅格化处理
拾取流体运动影像上任意静止图像,对图像进行栅格化处理,像素参数=10pix,得到图像栅格点阵。
如图4所示,本具体实施方法中,选取部分数据作为样本具体演示表面流速场测量过程。因面在图像栅格点阵中设定该数据样本对应位置点阵为测量区点阵,并用白色矩形框圈定出测量区点阵。
步骤S13、设置图像处理参数
设置图像处理速率=40ms/帧、帧间隔=1,像素参数与实际位移换算比例100pix=1m。
步骤S2、计算流体表面测量区在视觉坐标系xyz中的流速场
步骤S21、建立流体运动影像三维空间视觉坐标系xyz
在流体运动影像上建立三维空间视觉坐标系xyz获得影像视觉平面xoy。图5是视觉坐标系xyz建立示意图。
步骤S22、计算流体运动图像稠密光流
拾取流体运动影像播放中任意两帧图像,第一帧图像时刻t=360ms,两帧图像间隔时间Δt=40ms。某像素栅格点(演示点)在第一帧图像的坐标系xyz中(i,j)=(25.10,40.20),到第二帧图像时移动了(Δi,Δj)=(9.60,6.40)。由此计算出该演示点在两帧图像间的平面二维图像速度(Vx=0.024pix/s,Vy=0.016pix/s)。图6是平面演示点二维跟踪示意图。
同理,可以计算出两帧图像中所有像素栅格点的图像速度,得到流体稠密光流其中(Δi,Δj)是点(i,j)在时间Δt内经过的位移,分别是在x方向和y方向的变化率,即是所得光流。样本数据所在测量区像素栅格点图像速度计算结果如图7显示。
步骤S23、计算图像点阵在视觉坐标系xyz中的像素流速场
计算图像中所有点阵在图像视觉坐标系xyz中的移动速度,得到式1所示视觉平面像素流速场,能够测算出流体表面所有点在视觉平面上的二维移动速度,显示效果图为像素级别的流速场图,且每一帧图像的像素级别流速场图均可得到。图8为第111帧、第4.48s时刻流速场图。
步骤S24、计算图像点阵在图像视觉坐标系xyz中的视觉投影流速场
根据已设置的像素参数与实际位移换算比例(100pix=1m)以及图像处理速率40ms/帧,将像素视觉平面流速场换算为测量区点阵在图像视觉坐标系xyz中的实际速度。计算得到上述演示点的实际二维速度为(Vx=2.4m/s,Vy=1.6m/s)。由此类推,得到流体表面在图像视觉坐标系xyz中的视觉投影流速场。
步骤S3、计算流体表面二维流速场
步骤S31、计算投影变换矩阵PPM
在影像上建立流体表面三维空间真实坐标系x′y′z′,获得流体表面平面x′oy′。本实施方法中,y′轴为流体流向、z′轴与x′oy′平面正交(图8)。
如图10所示,矩形参照物A″B″C″D″在静止图像上成像为伪矩形A′B′C′D′,确定静止图像上伪矩形A′B′C′D′位置,通过计算机读取四角点在视觉平面xoy上的图像坐标A′(i1′,j1′)、B′(i2′,j2′)、C′(i3′,j3′)、D′(i4′,j4′)。在视觉平面xoy中建立投影矩形ABCD,投影矩形ABCD四角点在视觉坐标系xyz中坐标分别是A(i1,j1)、B(i2,j2)、C(i3,j3)、D(i4,j4)。
设置投影矩形ABCD与伪矩形A′B′C′D′底边重合,即(i1,j1)与(i1′,j1′)重合,(i2,j2)与(i2′,j2′)与重合,确定边长AB=0.8m。依式7计算确定边长BC=3.2m及C(i3,j3)、D(i4,j4)坐标。
根据A(i1,j1)、B(i2,j2)、C(i3,j3)、D(i4,j4)与A′(i1′,j1′)、B′(i2′,j2′)、C′(i3′,j3′)、D′(i4′,j4′)的映射关系建立投影矩形ABCD与伪矩形A′B′C′D′间的透视投影变换矩阵PPM。
步骤S32、计算流体表面二维流速场
依式9计算比例尺k=1/0.8=1.25。
根据投影变换矩阵PPM及比例尺k,将步骤S2所得点阵在图像视觉坐标系xyz中的视觉投影流速场换算为其在真实坐标系x′y′z′中的实际流速场,即得到流体表面二维流速场。
以步骤S22中演示点为例,将位于伪矩形区域内的演示点(i,j)坐标乘以透视投影矩阵PPM,则得到该演示点在流体表面真实平面x′oy′内的真实坐标。
演示点真实坐标(i0,j0,z0)=(i,j,z)×PPM=(25.10,40.20,0)×PPM=(24.80,41.00,0)。再根据演示点在在流体表面真实坐标系x′y′z′中的坐标变化计算该点的实际速度,演示点在流体表面的实际速度为步骤S24中已经测算出演示点在图像视觉坐标系xyz中的视觉投影速度(Vx=2.4m/s,Vy=1.6m/s,Vz=0m/s),根据PPM的矩阵反算演示点在流体表面真实坐标系x′y′z′的实际速度为:投影变换速度×PPM×k=(Vx=2.4m/s,Vy=1.6m/s,Vz=0m/s)×PPM×k=(Vx=2.8m/s,Vy=3.2m/s,Vz=0m/s)。
同理,可以计算出伪矩形范围内其他像素栅格点在流体表面真实坐标系x′y′z′中的二维实际速度,进而得到整个流体表面的二维流速场(图11)。
Claims (7)
1.一种流体表面二维流速场测量方法,其特征在于:依以下步骤实施:
步骤S1、流体运动影像获取及初始化处理
步骤S11、现场设备安装
在运动流体测量现场安装数字拍摄设备,使拍摄镜头自正上方垂直向下拍摄获得运动流体正投影方向上的流体运动影像;
步骤S12、获取流体运动影像
获取流体运动影像,将影像输入数字处理设备;
步骤S13、图像栅格化处理
拾取流体运动影像上任意静止图像,设置像素参数对图像进行栅格化处理,得到图像栅格点阵;
步骤S14、设置图像处理参数
设置图像处理速率、图像实际处理分辨率、像素参数与实际位移换算比例;
步骤S2、计算流体表面的二维流速场
步骤S21、建立流体运动影像三维空间视觉坐标系xyz
在流体运动影像上建立三维空间视觉坐标系xyz获得视觉平面xoy;
步骤S22、计算流体运动图像稠密光流
拾取流体运动影像播放中任意两帧图像,两帧图像间隔时间△t,计算两帧图像上所有对应像素栅格点的稠密光流其中(△i,△j)是测量区内点(i,j)在时间△t内在三维空间视觉坐标系xyz中经过的位移;
步骤S23、计算图像点阵在三维空间视觉坐标系xyz中的像素流速场
计算所有点阵在三维空间视觉坐标系xyz中的移动速度,得到式1所示视觉平面像素流速场:
式中,—点阵在视觉平面像素流速场,
—点(i,j)的灰度值ft(i,j)的偏导数;
步骤S24、计算点阵在视觉平面xoy中的视觉投影流速场
根据像素参数与实际位移换算比例,将步骤S23所得点阵在视觉平面中的像素流速场换算得到点阵在视觉平面xoy中的视觉投影流速场;
步骤S24所得视觉平面xoy中的视觉投影流速场即为测量得到的流体表面二维流速场。
2.一种流体表面二维流速场测量方法,其特征在于:依如下步骤实施:
步骤S1、流体运动影像获取及初始化处理
步骤S11、现场设备安装
在运动流体测量现场安装数字拍摄设备,从迎流体运动方向上的任意角度获取流体运动影像;在测量流体附近能够被数字拍摄设备完全拍摄且不影响流体运动的位置设置矩形参照物A″B″C″D″,确定矩形参照物A″B″C″D″边长A″B″、C″D″长度;
步骤S12、获取流体运动影像
获取流体运动影像,确保影像中有矩形参照物A″B″C″D″的完整图像,将影像输入数字处理设备;
步骤S13、图像栅格化处理
拾取流体运动影像上任意静止图像,设置像素参数对图像进行栅格化处理,得到图像栅格点阵;
步骤S14、设置图像处理参数
设置图像处理速率、图像实际处理分辨率、像素参数与实际位移换算比例;
步骤S2、计算视觉投影流速场
步骤S21、建立流体运动影像三维空间视觉坐标系xyz
在流体运动影像上建立三维空间视觉坐标系xyz获得视觉平面xoy;
步骤S22、计算流体运动图像稠密光流
拾取流体运动影像播放中任意两帧图像,两帧图像间隔时间△t,计算两帧图像上所有对应像素栅格点的稠密光流其中(△i,△j)是测量区内点(i,j)在时间△t内在三维空间视觉坐标系xyz中经过的位移;
步骤S23、计算图像点阵在三维空间视觉坐标系xyz中的像素流速场
计算所有点阵在三维空间视觉坐标系xyz中的移动速度,得到式1所示视觉平面像素流速场:
式中,—点阵在视觉平面像素流速场,
—点(i,j)的灰度值ft(i,j)的偏导数;
步骤S24、计算点阵在视觉平面xoy中的视觉投影流速场
根据像素参数与实际位移换算比例,将步骤S23所得点阵在视觉平面中的像素流速场换算得到点阵在视觉平面xoy中的视觉投影流速场;
步骤S3、运动流体图像透视投影变换
步骤S31、确定伪矩形A′B′C′D′坐标:
拾取流体运动影像上任意静止图像,矩形参照物A″B″C″D″在静止图像上对应伪矩形A′B′C′D′;在静止图像上建立三维空间视觉坐标系xyz,计算机确定伪矩形A′B′C′D′四角点在视觉平面xoy中的图像坐标A′(i′1,j′1)、B′(i′2,j′2)、C′(i′3,j′3)、D′(i′4,j′4);
步骤S32、设置投影矩形ABCD
在视觉平面xoy中建立投影矩形ABCD,投影矩形ABCD四角点在三维空间视觉坐标系xyz中坐标分别是A(i1,j1)、B(i2,j2)、C(i3,j3)、D(i4,j4);所述A(i1,j1)与A′(i′1,j′1)重合、B(i2,j2)与B′(i′2,j′2)重合,确定边长AB;根据式7确定边长BC及C(i3,j3)、D(i4,j4)坐标;
步骤S33、计算透视投影变换矩阵PPM
根据A(i1,j1)、B(i2,j2)、C(i3,j3)、D(i4,j4)与A′(i′1,j′1)、B′(i′2,j′2)、C′(i′3,j′3)、D′(i′4,j′4)的映射关系建立投影矩形ABCD与伪矩形A′B′C′D′间的透视投影变换矩阵PPM:
步骤S4、计算流体表面二维流速场
步骤S41、计算比例尺k
依式9计算比例尺k
步骤S42、计算流体表面二维流速场
根据投影变换矩阵PPM及比例尺k,将步骤S2所得点阵在三维空间视觉坐标系xyz视觉平面xoy中的视觉投影流速场换算为其在真实坐标系x'y'z'真实平面x′oy′中的流速场,即得到流体表面二维流速场。
3.根据权利要求1或2所述的测量方法,其特征在于:所述像素参数设置为1pix~100pix,图像处理速率设置为10ms/帧~100ms/帧,帧间隔1~10。
4.根据权利要求3所述的测量方法,其特征在于:所述帧间隔1~5。
5.根据权利要求1、2、4任一所述的测量方法,其特征在于:所述流体是泥石流。
6.一种运动流体图像透视投影变换计算方法:依如下步骤实施:
步骤S1a、设置矩形参照物A″B″C″D″
在运动流体测量现场安装数字拍摄设备,从迎流体运动方向上的任意角度获取流体运动影像;在测量流体附近能够被数字拍摄设备完全拍摄且不影响流体运动的位置设置矩形参照物A″B″C″D″,确定矩形参照物A″B″C″D″边长A″B″、C″D″长度;
步骤S2a、确定伪矩形A′B′C′D′坐标
获取流体运动影像,确保影像中有矩形参照物A″B″C″D″的完整图像;拾取流体运动影像上任意静止图像,矩形参照物A″B″C″D″在静止图像上对应伪矩形A′B′C′D′;在流体运动影像上建立三维空间视觉坐标系xyz获得视觉平面xoy,计算机确定伪矩形A′B′C′D′四角点在视觉平面xoy中的图像坐标A′(i′1,j′1)、B′(i′2,j′2)、C′(i′3,j′3)、D′(i′4,j′4);
步骤S3a、设置投影矩形ABCD
在视觉平面xoy中建立投影矩形ABCD,投影矩形ABCD四角点在三维空间视觉坐标系xyz中坐标分别是A(i1,j1)、B(i2,j2)、C(i3,j3)、D(i4,j4);所述A(i1,j1)与A′(i′1,j′1)重合、B(i2,j2)与B′(i′2,j′2)重合,确定边长AB;根据式7确定边长BC及C(i3,j3)、D(i4,j4)坐标;
步骤S4a、计算透视投影变换矩阵PPM
根据A(i1,j1)、B(i2,j2)、C(i3,j3)、D(i4,j4)与A′(i′1,j′1)、B′(i′2,j′2)、C′(i′3,j′3)、D′(i′4,j′4)的映射关系建立投影矩形ABCD与伪矩形A′B′C′D′间的透视投影变换矩阵PPM:
步骤S5a、计算比例尺k
依式9计算比例尺k
步骤S6a、计算静止图像中点的实际坐标
将静止图像上任意点的图像坐标经透视投影变换矩阵PPM变换得到该点在视觉平面xoy中的投影坐标,再将投影坐标与比例尺k相乘得到静止图像上任意点在流体表面真实平面x′oy′中的实际坐标。
7.根据权利要求6所述的运动流体图像透视投影变换计算方法,其特征在于:所述流体是泥石流。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410617463.2A CN104297516B (zh) | 2014-11-06 | 2014-11-06 | 一种流体表面二维流速场测量方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410617463.2A CN104297516B (zh) | 2014-11-06 | 2014-11-06 | 一种流体表面二维流速场测量方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104297516A CN104297516A (zh) | 2015-01-21 |
CN104297516B true CN104297516B (zh) | 2017-07-18 |
Family
ID=52317327
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201410617463.2A Expired - Fee Related CN104297516B (zh) | 2014-11-06 | 2014-11-06 | 一种流体表面二维流速场测量方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104297516B (zh) |
Families Citing this family (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105588546A (zh) * | 2016-03-15 | 2016-05-18 | 清华大学 | 一种测量浆体形状参数及其随时间变化的方法及系统 |
CN107102165B (zh) * | 2017-04-14 | 2020-03-20 | 重庆大学 | 一种基于粒子图像测速的表面流场测量方法 |
CN107084926B (zh) * | 2017-05-04 | 2019-10-11 | 北京清环智慧水务科技有限公司 | 液体清澈度检测方法、系统和装置 |
KR101996992B1 (ko) * | 2018-11-13 | 2019-07-08 | 주식회사 하이드로셈 | 옵티컬 플로우 영상 처리를 이용하는 하천 유속 측정 장치 및 방법 |
KR101978351B1 (ko) * | 2018-11-13 | 2019-05-15 | 주식회사 하이드로셈 | Cctv 영상 기반의 실시간 자동 유량계측 시스템 및 방법 |
KR101975476B1 (ko) * | 2018-11-13 | 2019-05-08 | 주식회사 하이드로셈 | 측위 자료 필터링을 통한 실시간 하천 수위 측정 장치 및 방법 |
CN111141927B (zh) * | 2019-12-31 | 2021-07-13 | 清华大学 | 泥石流示踪粒、内部流速实验系统、内部流场测算方法 |
CN113804917A (zh) * | 2021-09-17 | 2021-12-17 | 山东新一代信息产业技术研究院有限公司 | 一种基于点扩散估计的河流表面流速测量方法 |
CN114877954B (zh) * | 2022-07-12 | 2022-09-23 | 杭州春来科技有限公司 | 一种固定污染源测量方法及系统 |
CN116183957A (zh) * | 2023-04-23 | 2023-05-30 | 中国科学院合肥物质科学研究院 | 一种超导磁体间隙填充材料流速测量系统及方法 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102288778A (zh) * | 2011-05-16 | 2011-12-21 | 西南交通大学 | 基于网格对象跟踪的实时泥石流表面速度测量方法 |
CN103996171A (zh) * | 2014-05-05 | 2014-08-20 | 河海大学 | 基于时空图像的流体运动矢量估计方法 |
-
2014
- 2014-11-06 CN CN201410617463.2A patent/CN104297516B/zh not_active Expired - Fee Related
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102288778A (zh) * | 2011-05-16 | 2011-12-21 | 西南交通大学 | 基于网格对象跟踪的实时泥石流表面速度测量方法 |
CN103996171A (zh) * | 2014-05-05 | 2014-08-20 | 河海大学 | 基于时空图像的流体运动矢量估计方法 |
Non-Patent Citations (4)
Title |
---|
Large-scale particle image velocimetry for flow analysis in hydraulic engineering applications;ICHIRO FUJITA 等;《JOURNAL OF HYDRAULIC RESEARCH》;19981231;第36卷(第3期);第397-414页 * |
Quantification of swash flows using video-based particle image velocimetry;K.T. Holland 等;《Coastal Engineering》;20011231;第44卷(第2期);第65-77页 * |
Using Line Integral Convolution for Flow Visualization: Curvilinear Grids, Variable-Speed Animation, and Unsteady Flows;Lisa K. Forssell 等;《IEEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS》;19950630;第1卷(第2期);第133-141页 * |
基于图像处理的山区河道表面流场测算研究;李蔚 等;《人民长江》;20140831;第45卷(第15期);第89-92、104页 * |
Also Published As
Publication number | Publication date |
---|---|
CN104297516A (zh) | 2015-01-21 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104297516B (zh) | 一种流体表面二维流速场测量方法 | |
KR101996992B1 (ko) | 옵티컬 플로우 영상 처리를 이용하는 하천 유속 측정 장치 및 방법 | |
CN110118528B (zh) | 一种基于棋盘靶标的线结构光标定方法 | |
CN105354875B (zh) | 一种室内环境二维与三维联合模型的构建方法和系统 | |
CN106289106B (zh) | 一种线阵相机和面阵相机相结合的立体视觉传感器及标定方法 | |
CN105066909B (zh) | 一种手持式多激光条纹快速三维测量方法 | |
CN103439230B (zh) | 一种气泡参数测量方法及测量装置 | |
CN107525945B (zh) | 基于集成成像技术的3d-3c粒子图像测速系统及方法 | |
KR102523451B1 (ko) | 드론 영상 기반의 하천 유속 측정을 위한 장치 및 방법 | |
CN107102165B (zh) | 一种基于粒子图像测速的表面流场测量方法 | |
KR101305305B1 (ko) | 시공간 영상의 상관 분석을 이용한 표면 유속 측정 시스템 및 방법 | |
CN102721376A (zh) | 一种大视场三维视觉传感器的标定方法 | |
CN109559355A (zh) | 一种基于相机组的无公共视场的多相机全局标定装置及方法 | |
CN103903263B (zh) | 一种基于Ladybug全景相机图像的360度全方位测距方法 | |
CA2529044A1 (en) | Three-dimensional modeling from arbitrary three-dimensional curves | |
Yuan et al. | Development of a robust Stereo-PIV system for 3-D soil deformation measurement | |
CN110207666A (zh) | 一种气浮平台上模拟卫星的视觉位姿测量方法及装置 | |
CN112785654A (zh) | 轨道几何检测系统标定方法及装置 | |
CN109855559A (zh) | 一种全空间标定系统及方法 | |
CN116047104A (zh) | 一种基于转镜式高速相机的运动目标速度测量方法 | |
CN115393537A (zh) | 输电通道三维可视化建模的精准性评估系统及方法 | |
JP3512894B2 (ja) | 相対的移動量算出装置及び相対的移動量算出方法 | |
Alekseev et al. | Visual-inertial odometry algorithms on the base of thermal camera | |
TW201620698A (zh) | 獲取物體三維模型的方法及設備 | |
CN109932281A (zh) | 基于视觉的液体黏度在线测量方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20170718 Termination date: 20181106 |