基于线阵CCD的风浪环境下水面弱流场速度测量方法

2014-03-29王小青陈鹏真种劲松

实验流体力学 2014年2期

刘利, 王小青, 陈鹏真, 种劲松

(1. 中国科学院电子学研究所微波成像技术重点实验室, 北京 100190; 2. 中国科学院大学, 北京 100190)

0 引言

大量实验和理论研究表明，水下地形、内波和漩涡等海洋现象都是通过其产生的表面弱流场与表面波的相互作用而被微波遥感手段所观测到[1]。而对风浪环境下的弱流场进行测量是研究波流相互作用机理的核心步骤。因此，研究风浪环境下表面弱流场测量的方法和装置对于研究波流相互作用机理具有重要意义。

基于图像处理技术的传统流场测量方法，需要在水面散布示踪粒子，根据示踪粒子的移动来测量流场[2-3]，这类方法称为PIV(Particle Image Velocimetry，粒子图像测速)方法。风浪环境下，风生波会对示踪粒子的运动产生严重影响，用PIV方法测得的是波浪振动速度与表面弱流场的叠加结果，难以将二者区分。而需要测量的表面弱流场是与风生波无关的，因此传统的流场测量方法无法工作在风浪环境下，难以对流场与风生波相互作用展开研究。

本文提出的基于线阵CCD(Charge-coupled Device，电荷耦合组件)的流场测量方法作为一种非接触式的光学测量方法，在风浪环境下可以高精度地获取表面弱流场信息而不受波浪振动的影响。本方法是在基于线阵CCD的波浪测量方法上发展而来。基于线阵CCD的波浪测量方法是近年来在国际上出现的一种高精度、高速的空时二维测量方法[4-7]，其基本原理是在特定光照条件下，水面的光强与水面的倾角有一一对应的映射关系，通过获取水面光强可以反演出水面斜率，进而反演出水面波高。之前的研究主要集中于使用基于线阵CCD或面阵CCD的光学方法对波浪斜率、波高分布以及波谱、方向谱和能量谱等进行测量[4-11]。目前尚无利用线阵CCD测量水面流场的相关报告。

线阵CCD由于其时间、空间分辨率高，可以对振动频率很高的微尺度波进行测量，而且可以获得波浪相位的连续变化。本文所提出的测量方法就是根据表面波浪的相位变化获得表面波的相速度，其值是波的本征相速度与外界流场的叠加，外界流场包括风生流场及所要测量的表面弱流场。根据经典的波浪理论，波的本征相速度与波浪的振幅高低无关；风速稳定后风生流场也是不变的，因此可以通过表面波相速度的变化计算表面弱流场速度。最后用此方法对水槽实验获得的实际CCD数据进行处理，验证了本方法的正确性。

1 基于线阵CCD的风浪环境下水面流速测量方法

风浪环境下，表面波的相速度可以表示为[12]

c=c0+ud+uc

(1)

式中：c0为本征相速度，ud为风引起的风生流场速度，uc即为所要测量的表面弱流场速度。根据经典波浪理论，波的本征相速度为：

(2)

式中：g为重力加速度，T为表面张力系数，ρ为水的密度，k为表面波波数。由公式(2)可知，波的本征相速度只与波长有关，与波的振幅无关。而且风速稳定后风生流场的速度也是稳定不变的，因此可以通过表面波相速度的变化来测量表面弱流场速度。

表面波的相速度c为波长λ与周期T之比

c=λf

(3)

(4)

式中：波数k=2π/λ，角频率ω=2πf=2π/T。

因此，为了计算某一波长表面波的相速度，可以有以下两种方法：

(1) 波频率法

此方法是计算某一波长的波频率f，通过公式(3)得到相速度。在计算频率f时，需要对斜率数据进行时间傅里叶变换，频率分辨率与所选时间轴长度有关。由于频率分辨率Δf与时间轴长度ts关系为：Δf=1/ts，因此频率分辨率若要达到1Hz，时间轴长度须选为1s；频率分辨率若要达到0.1Hz，时间轴长度须选为10s。因此采用此方法要达到高的频率分辨就必然会损失时间精度，但此方法的抗噪性能较好。

(2) 波相位法

此方法是计算相位差Δφ，通过公式(4)计算相速度。本方法的优势是时间分辨率很高，但抗噪性能较差。

由于本研究需要提取较高时间分辨率的流场信息，因此采用第二种计算方法，数据处理中为了减小噪声对测量结果的影响，对数据进行自适应的二维频域滤波，对相位差分时间间隔进行了优化，并对不同波长的表面波相速度测量结果进行统计平均得到最终较好的测量结果。

所提出的测量方法主要步骤如图1所示：

图1 算法总体流程图

下面详细描述图1中的4个步骤。

(1) 几何校正及光强斜率转换

在特定光照条件下，水面的光强与水面的倾角有一一对应的映射关系。在标定映射关系后，就可以通过线阵CCD相机获取水面的光强变化来反演出水面斜率。由于CCD相机的位置及光学相机拍摄物体的“近大远小”特征，在标定映射关系前，需要进行几何校正，将CCD相机拍摄图像中的像素位置准确投影到空间位置。

(2) 斜率噪声抑制

由于实验设备以及测量环境的影响，获得的数据中不可避免地会混杂有噪声。噪声通常较为均匀地分布在时-空二维频域；而水面波浪由于其特殊的频散关系，波浪斜率信号在时-空二维频域分布较为集中。因此利用此原理可以通过二维频域自适应滤波的方法大幅度提高波浪斜率数据的信噪比。

假设线阵CCD扫描仪测得的斜率表示为s(x,t)，其中x表示空间位置，t表示时间，则斜率的自相关函数可以表示为：

Rs(ξ,τ)=〈s(x,t)s(x+ξ,t+τ)〉

(5)

斜率谱S(kx,ω)为

∬Rs(ξ,τ)e-j(kxξ+ωτ)dξdτ

(6)

式中：kx为表面波的波数矢量，ω为表面波的角频率。根据经典的波浪理论，斜率数据在频域中主要集中于频散关系曲线附近，

(7)

式中：ρ为水体密度，T为表面张力系数，g为重力加速度。可用频域二维滤波器对斜率二维频谱滤波，滤波范围可自适应地选择为其谱宽。

(3) 相位信息提取及修正

假设在时刻t0，沿距离向的斜率分布为s(x,t0)，其空间傅里叶变换为S(k,t0)。傅里叶变换之后可以对不同波数的表面波相位进行计算。对于波长为λ0(波数为k0=2π/λ0)的表面波，假设其傅里叶变换形式为：

S(k0,t0)=SR+jSI

(8)

其相位为：

(9)

由于线阵CCD相机自身存储空间有限，在其离线采集完一帧数据后需要将数据传输回计算机才可以开始下一帧数据的采集。因此相对于帧数据内的时间间隔，帧数据与帧数据间的时间间隔并不稳定。由此导致的结果是相位数据在帧数据间会有小幅度的跳变，这种跳变使得相位差在帧数据切换时刻变化比较剧烈。使用未经修正的相位数据计算出的相速度会出现明显异常，因此需要对原始的相位数据进行修正。

具体修正方法为：由相位数据得到相位差后，帧数据切换时刻的相位差修正为采用其相邻数据经过插值后得到的新的相位差。

(4) 表面波相速度计算

计算得到相位数据后，根据公式(4)就可以计算表面波相速度。对于波长为λ0(波数为k0=2π/λ0)的表面波，假设在时刻t0其相位为φ(k0,t0)，在时刻t1=t0+Δt其相位为φ(k0,t1)，则其角频率为

(10)

在根据公式(4)计算时，由于存在相位噪声的影响，选择不同的时间间隔Δt计算得出的相速度精度是不同的。时间间隔过小，相机内部的热噪声影响较大，使得最终计算结果中误差占比过大；加大时间间隔可以减小热噪声的影响，但加大时间间隔后由于波浪的形变会导致表面波的相干性降低，相位噪声增加；此外加大时间间隔还会引起时间分辨率的下降。因此在时间间隔的选择上需要折中考虑。

以标准差表示精度，标准差越小，精度越高。由公式(4)可知，

(11)

其中，角频率精度为

(12)

为了选择合适的时间间隔，2.2节中描述了实验中所选时间间隔与测量精度的关系。

在本文第2节中，将以实验数据为例对图1的算法流程进行详细的说明。

2 实验验证

2.1实验装置

实验于2013年3月在解放军理工大学分层流实验水槽中进行。该水槽的长、宽和高分别为12，1.5和1.2m。系统在实验过程中会产生内波，内波是指发生在密度稳定层化的流体内部的一种波动。本实验中使用盐水和淡水产生密度分层，如图2所示。内波的传播会在表面产生流场。实验中风速可变，内波振幅为10cm，传播方向与水槽长度方向平行，内波传播过程中遇到水槽壁会发生反射现象。图2为实验装置示意图。

图中L表示有效扫描区域。通过特殊的光源以及散光板产生特殊的光照条件，使得水面光强与表面波斜率有近似的一一映射关系。光线经过水面反射后进入图像采集器件——线阵CCD相机。实验中CCD相机距离水面的高度为44cm，观测角度(与垂直方向夹角)为34°。

图2 水槽内波观测实验光学装置示意图

光学系统的主要参数如下：有效扫描长度：L=23cm；分辨率：0.38mm(空间分辨率)，1/300s(时间分辨率)；扫描频率：300Hz；扫描点数：2048。

实验中采用了6台CCD相机进行并行数据采集,图3是实际实验装置图。

图3 实验装置图

2.2实验结果及分析

以下给出第1节中所述的4个步骤的具体处理过程及部分结果。

(1) 几何校正及光强斜率转换

实验中得到的时间-空间二维图像数据如图4所示。横坐标x表示距离(单位：cm)，纵坐标t表示时间(单位：s)，灰度值大小表示经过水面反射后进入相机的光强。每幅图像的尺寸为300pixel×2048pixel，扫描频率为300Hz，因此一幅图代表1s时间内扫描区域内波动变化。表面波的运动在图像中表现为倾斜的条纹，倾斜条纹的斜率即表示表面波的传播速度。

CCD相机通过记录扫描区域的光强度从而获得扫描区域的实时波斜率，如图5所示。CCD在空间上的扫描范围约为23cm，每个采样点的间隔约为0.23m/600=0.38mm,那么CCD可以记录的表面波的波数范围理论上为30～8200rad/m。

比较图4与5，图4中的倾斜条纹存在一定的几何畸变(条纹弯曲)，经过几何校正后这种畸变消失(图5)；而且图4中光强分布沿距离方向是不均匀的，但经过光强斜率转换后的斜率分布比较均匀。

图4 实验获得的CCD图像数据

图5 由图4计算得出的波斜率分布

(2) 斜率噪声抑制

实验中测得的斜率s(x,t)和斜率谱S(kx,ω)分别如图5和6所示。图6中的黑色曲线为根据频散关系(公式(7))得出的频率与波数关系曲线。图6中波频集中区域与频散关系曲线吻合很好，这一方面验证了实验所获取数据的有效性，另一方面也表明了频率成分主要分布于频散关系曲线附近，因此可用图7所示的自适应二维频域滤波器进行滤波处理，滤波后的斜率分布如图8所示，与图5相比斜率分布更加清晰。

某一距离点上波斜率随时间的变化如图9所示。

图6 斜率谱

图7 频域二维滤波器

图8 滤波后的斜率分布

图9 1s内波斜率随时间变化

(3) 相位信息提取及修正

以波长λ=4.57cm为例，在观测时间第28s内相位变化如图10所示。由于直接计算所得相位值位于[-π,π]区间，在计算相位差时需要对相位进行解卷绕。

图10 第28s内相位变化

图11为某波长表面波的相位在一段时间内的变化，可以看出在帧数据与帧数据切换时相位会有小幅度跳变，这个跳变使得相位差变化非常剧烈。

由第1节介绍的方法修正后的相位及相位差如图12所示，与图11相比较有明显的改善。

(4) 表面波相速度计算

在计算得到相位数据后，可由公式(4)计算不同波长的表面波相速度。计算过程中对6台CCD相机的结果进行了平均处理。

如第1节所述，表面波相速度的测量精度主要与角频率的精度有关，而角频率的精度与所选时间间隔有关。图13为波长分别为2.08和1.52cm的表面波相速度测量精度与所选时间间隔的关系。

由图13可以明显看出，在选择时间间隔Δt≥2.5s时，表面波相速度测量精度优于1cm/s。

(a)

(b)

(a)

(b)

图14为对波长为4.57cm的表面波采用不同的时间间隔时计算得出的相速度结果。

由图14可以看出，间隔为2.5s的计算结果要明显优于间隔为0.5s的计算结果。

计算时采用时间间隔2.5s，经过平滑滤波后得到的波长分别为2.08cm、1.52cm和1.14cm的表面波的相速度如图15(a)、(b)和(c)所示。

由图15及图14(b)可以看出，速度有4次明显的变化(图14(b)中A、B、C、D四点)。这是因为在实验过程中在水槽中有内孤立波经过探测区域，内波的传播引起了表面波相速度上的变化。图15及图14(b)中虚线表示在内波经过之前(前80s)表面波相速度均值。这4次表面波相速度的变化与实验中观察到的现象十分吻合，即内波首先逆风向经过探测区域(100s左右)，而后经过3次反射。逆风向经过探测区域时内波传播使得表面波相速度下降(图14(b)中A、C)，顺风向经过探测区域时内波传播使得表面波相速度增大(图14(b)中B、D)。而且由图15可以看出，虽然不同波长下相速度大小不同，但在内波经过时引起的相速度变化幅值是基本一致的，约为3～4cm/s，整体幅值变化约为6～7cm/s，这种表面波相速度上的变化就是内波传播引起的表面弱流场导致的。

(a) 波长为2.08cm

(b) 波长为1.52cm

(a) Δt=0.5s

(b) Δt=2.5s

(a) 2.08cm

(b) 1.52cm

图16为不同风速下内波经过探测区域前(前80s)不同波长表面波相速度对比。

可以看出，对于波长相同的表面波，风速越大测得的相速度越大。这是由于风速越大，风生流越大，根据公式(1)，表面波相速度越大。同时，可以看出随着风速增大，相速度曲线在图中近似于“平移”，这也印证了波的本征相速度与风速无关，亦即与波高无关。

得到如图15所示的表面波相速度后，将此相速度与前80s(内波经过前)的相速度均值相减，便可得到内波传播所引起的表面弱流场。由于可以计算不同波长的表面波相速度，将多个波长表面波计算得出的表面弱流场进行平均处理，得到结果如图17所示。

图16 不同风速下不同波长表面波相速度(Vwind为风速)

图17 内波激发的表面弱流场速度(正值表示内波与表面波运动方向相同)

由图17可以看出，在前80s左右，内波尚未经过观测区域时，流场速度在零值上下波动且变化幅度很小。对前80s的数据进行统计，其精度达到0.3cm/s，足以应用到对波-流相互作用机理的研究中。

3 结论

提出了一种基于线阵CCD的水面流场光学测量方法，其最大的优点是可以在风浪环境下高精度地获取表面弱流场信息而不受波浪振动的影响，克服了传统PIV流场测量方法不能工作于风浪环境下的缺点。应用此方法对水槽内波观测实验数据的处理结果表明，表面弱流场的测量精度达到0.3cm/s。已有学者利用线阵CCD测量表面波浪谱等信息[4-7]，而本文在此基础上进一步同时测量表面弹流场速度，并且保持实验中流场方向基本与水槽长度方向平行，因此可以同时获得风浪环境下的表面弱流场与波浪信息，这在波-流相互作用机理的研究中将会发挥重要作用。

所采用的实验设备及方法也存在一定的局限性。实验中对光照条件要求较高，这在实验室环境中容易实现，但实际测量时测量条件较为恶劣，这对该方法的实用性会造成不利影响。因此今后的研究内容主要是改进该方法使其适用于实际测量以及讨论该方法应用于复杂的随机波浪表面测量可能面临的问题和对测量精度造成的影响。

参考文献:

[1]何宜军. 成像雷达海浪成像机制[J]. 中国科学: D辑, 2000, 30(5): 554-560.

He Yijun. Imaging rader imaging mechanism of ocean waves[J]. Science in Chian, 2000, 30(5): 554-560.

[2]王浩, 曾理江. 二维及三维流场的光学测量方法[J]. 光学技术, 2001, 27(2): 139-142.

Wang Hao, Zeng Lijiang. Optical methods for measuring 2D and 3D flow field[J]. Optical Technology, 2001, 27(2): 139-142.

[3]Willert C E, Gharib M. Digital particle image velocimetry[J]. Experiments in Fluids, 1991, 10(4): 181-193.

[4]Titov V I, Bakhanova V V, Kemarskaja O N, et al. Investigation of sea roughness with complex of optical devices[C]//Proceeding SPIE 7473, Remote Sensing of the Ocean, Sea Ice and Large Water Regions, 2009, 74730T.

[5]Titov V I, Zuikova E M, Luchinin A G, et al. Investigation of surface roughness with optical methods[C]//Proceeding SPIE 7825, Remote swnsing of the Ocear, Sea Ice and Large Water Regions, 2010, 78250G.

[6]Titov V, Bakhanov V, Zuikova E, et al. Remote sensing of water basins using optical range-time images of water surface[C]//Proceeding SPIE 8175, Remote Sensing of the Sea Ice, Coastal Water and Large Water Regions, 2011, 81751D.

[7]Titov V I, Bakhanov V V, Luchinin A G, et al. Remote sensing of sea surface features by optical RTI images[C]//Proceeding SPIE 8888, Remote Sensing ot the Ocean, Sea Ice, Coastal Water and Large Water Regions, 2013, 88880J.

[8]Stockdon H F, Holman R A. Estimation of wave phase speed and nearshore bathymetry from video imagery[J]. Journal of Geophysical Research: Oceans (1978-2012), 2000, 105(C9): 22015-22033.

[9]Jähne B, Klinke J, WAAS S. Imaging of short ocean wind waves: a critical theoretical review[J]. JOSA A, 1994, 11(8): 2197-2209.

[10] Stilwelljr D, Pilon R O. Directional spectra of surface waves from photographs[J]. Journal of Geophysical Research, 1974, 79(9): 1277-1284.

[11] Young I R, Rosenthal W, Ziemer F. A three-dimensional analysis of marine radar images for the determination of ocean wave directionality and surface currents[J]. Journal of Geophysical Research: Oceans (1978-2012), 1985, 90(C1): 1049-1059.

[12] Rozenberg A D, Ritter M J, Melville W K, et al. Free and bound capillary waves as microwave scatterers: Laboratory studies[J]. IEEE Transactions on Geoscience and Remote Sensing, 1999, 37(2): 1052-1065.

作者简介: