劉 利, 王小青, 陳鵬真, 種勁松

(1. 中國科學(xué)院電子學(xué)研究所微波成像技術(shù)重點實驗室, 北京 100190; 2. 中國科學(xué)院大學(xué), 北京 100190)

0 引 言

大量實驗和理論研究表明，水下地形、內(nèi)波和漩渦等海洋現(xiàn)象都是通過其產(chǎn)生的表面弱流場與表面波的相互作用而被微波遙感手段所觀測到[1]。而對風(fēng)浪環(huán)境下的弱流場進行測量是研究波流相互作用機理的核心步驟。因此，研究風(fēng)浪環(huán)境下表面弱流場測量的方法和裝置對于研究波流相互作用機理具有重要意義。

基于圖像處理技術(shù)的傳統(tǒng)流場測量方法，需要在水面散布示蹤粒子，根據(jù)示蹤粒子的移動來測量流場[2-3]，這類方法稱為PIV(Particle Image Velocimetry，粒子圖像測速)方法。風(fēng)浪環(huán)境下，風(fēng)生波會對示蹤粒子的運動產(chǎn)生嚴重影響，用PIV方法測得的是波浪振動速度與表面弱流場的疊加結(jié)果，難以將二者區(qū)分。而需要測量的表面弱流場是與風(fēng)生波無關(guān)的，因此傳統(tǒng)的流場測量方法無法工作在風(fēng)浪環(huán)境下，難以對流場與風(fēng)生波相互作用展開研究。

本文提出的基于線陣CCD(Charge-coupled Device，電荷耦合組件)的流場測量方法作為一種非接觸式的光學(xué)測量方法，在風(fēng)浪環(huán)境下可以高精度地獲取表面弱流場信息而不受波浪振動的影響。本方法是在基于線陣CCD的波浪測量方法上發(fā)展而來?；诰€陣CCD的波浪測量方法是近年來在國際上出現(xiàn)的一種高精度、高速的空時二維測量方法[4-7]，其基本原理是在特定光照條件下，水面的光強與水面的傾角有一一對應(yīng)的映射關(guān)系，通過獲取水面光強可以反演出水面斜率，進而反演出水面波高。之前的研究主要集中于使用基于線陣CCD或面陣CCD的光學(xué)方法對波浪斜率、波高分布以及波譜、方向譜和能量譜等進行測量[4-11]。目前尚無利用線陣CCD測量水面流場的相關(guān)報告。

線陣CCD由于其時間、空間分辨率高，可以對振動頻率很高的微尺度波進行測量，而且可以獲得波浪相位的連續(xù)變化。本文所提出的測量方法就是根據(jù)表面波浪的相位變化獲得表面波的相速度，其值是波的本征相速度與外界流場的疊加，外界流場包括風(fēng)生流場及所要測量的表面弱流場。根據(jù)經(jīng)典的波浪理論，波的本征相速度與波浪的振幅高低無關(guān)；風(fēng)速穩(wěn)定后風(fēng)生流場也是不變的，因此可以通過表面波相速度的變化計算表面弱流場速度。最后用此方法對水槽實驗獲得的實際CCD數(shù)據(jù)進行處理，驗證了本方法的正確性。

1 基于線陣CCD的風(fēng)浪環(huán)境下水面流速測量方法

風(fēng)浪環(huán)境下，表面波的相速度可以表示為[12]

c=c0+ud+uc

(1)

式中：c0為本征相速度，ud為風(fēng)引起的風(fēng)生流場速度，uc即為所要測量的表面弱流場速度。根據(jù)經(jīng)典波浪理論，波的本征相速度為：

(2)

式中：g為重力加速度，T為表面張力系數(shù)，ρ為水的密度，k為表面波波數(shù)。由公式(2)可知，波的本征相速度只與波長有關(guān)，與波的振幅無關(guān)。而且風(fēng)速穩(wěn)定后風(fēng)生流場的速度也是穩(wěn)定不變的，因此可以通過表面波相速度的變化來測量表面弱流場速度。

表面波的相速度c為波長λ與周期T之比

c=λf

(3)

(4)

式中：波數(shù)k=2π/λ，角頻率ω=2πf=2π/T。

因此，為了計算某一波長表面波的相速度，可以有以下兩種方法：

(1) 波頻率法

此方法是計算某一波長的波頻率f，通過公式(3)得到相速度。在計算頻率f時，需要對斜率數(shù)據(jù)進行時間傅里葉變換，頻率分辨率與所選時間軸長度有關(guān)。由于頻率分辨率Δf與時間軸長度ts關(guān)系為：Δf=1/ts，因此頻率分辨率若要達到1Hz，時間軸長度須選為1s；頻率分辨率若要達到0.1Hz，時間軸長度須選為10s。因此采用此方法要達到高的頻率分辨就必然會損失時間精度，但此方法的抗噪性能較好。

(2) 波相位法

此方法是計算相位差Δφ，通過公式(4)計算相速度。本方法的優(yōu)勢是時間分辨率很高，但抗噪性能較差。

由于本研究需要提取較高時間分辨率的流場信息，因此采用第二種計算方法，數(shù)據(jù)處理中為了減小噪聲對測量結(jié)果的影響，對數(shù)據(jù)進行自適應(yīng)的二維頻域濾波，對相位差分時間間隔進行了優(yōu)化，并對不同波長的表面波相速度測量結(jié)果進行統(tǒng)計平均得到最終較好的測量結(jié)果。

所提出的測量方法主要步驟如圖1所示：

圖1 算法總體流程圖

下面詳細描述圖1中的4個步驟。

(1) 幾何校正及光強斜率轉(zhuǎn)換

在特定光照條件下，水面的光強與水面的傾角有一一對應(yīng)的映射關(guān)系。在標定映射關(guān)系后，就可以通過線陣CCD相機獲取水面的光強變化來反演出水面斜率。由于CCD相機的位置及光學(xué)相機拍攝物體的“近大遠小”特征，在標定映射關(guān)系前，需要進行幾何校正，將CCD相機拍攝圖像中的像素位置準確投影到空間位置。

(2) 斜率噪聲抑制

由于實驗設(shè)備以及測量環(huán)境的影響，獲得的數(shù)據(jù)中不可避免地會混雜有噪聲。噪聲通常較為均勻地分布在時-空二維頻域；而水面波浪由于其特殊的頻散關(guān)系，波浪斜率信號在時-空二維頻域分布較為集中。因此利用此原理可以通過二維頻域自適應(yīng)濾波的方法大幅度提高波浪斜率數(shù)據(jù)的信噪比。

假設(shè)線陣CCD掃描儀測得的斜率表示為s(x,t)，其中x表示空間位置，t表示時間，則斜率的自相關(guān)函數(shù)可以表示為：

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

(5)

斜率譜S(kx,ω)為

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

(6)

式中：kx為表面波的波數(shù)矢量，ω為表面波的角頻率。根據(jù)經(jīng)典的波浪理論，斜率數(shù)據(jù)在頻域中主要集中于頻散關(guān)系曲線附近，

(7)

式中：ρ為水體密度，T為表面張力系數(shù)，g為重力加速度?？捎妙l域二維濾波器對斜率二維頻譜濾波，濾波范圍可自適應(yīng)地選擇為其譜寬。

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

假設(shè)在時刻t0，沿距離向的斜率分布為s(x,t0)，其空間傅里葉變換為S(k,t0)。傅里葉變換之后可以對不同波數(shù)的表面波相位進行計算。對于波長為λ0(波數(shù)為k0=2π/λ0)的表面波，假設(shè)其傅里葉變換形式為：

S(k0,t0)=SR+jSI

(8)

其相位為：

(9)

由于線陣CCD相機自身存儲空間有限，在其離線采集完一幀數(shù)據(jù)后需要將數(shù)據(jù)傳輸回計算機才可以開始下一幀數(shù)據(jù)的采集。因此相對于幀數(shù)據(jù)內(nèi)的時間間隔，幀數(shù)據(jù)與幀數(shù)據(jù)間的時間間隔并不穩(wěn)定。由此導(dǎo)致的結(jié)果是相位數(shù)據(jù)在幀數(shù)據(jù)間會有小幅度的跳變，這種跳變使得相位差在幀數(shù)據(jù)切換時刻變化比較劇烈。使用未經(jīng)修正的相位數(shù)據(jù)計算出的相速度會出現(xiàn)明顯異常，因此需要對原始的相位數(shù)據(jù)進行修正。

具體修正方法為：由相位數(shù)據(jù)得到相位差后，幀數(shù)據(jù)切換時刻的相位差修正為采用其相鄰數(shù)據(jù)經(jīng)過插值后得到的新的相位差。

(4) 表面波相速度計算

計算得到相位數(shù)據(jù)后，根據(jù)公式(4)就可以計算表面波相速度。對于波長為λ0(波數(shù)為k0=2π/λ0)的表面波，假設(shè)在時刻t0其相位為φ(k0,t0)，在時刻t1=t0+Δt其相位為φ(k0,t1)，則其角頻率為

(10)

在根據(jù)公式(4)計算時，由于存在相位噪聲的影響，選擇不同的時間間隔Δt計算得出的相速度精度是不同的。時間間隔過小，相機內(nèi)部的熱噪聲影響較大，使得最終計算結(jié)果中誤差占比過大；加大時間間隔可以減小熱噪聲的影響，但加大時間間隔后由于波浪的形變會導(dǎo)致表面波的相干性降低，相位噪聲增加；此外加大時間間隔還會引起時間分辨率的下降。因此在時間間隔的選擇上需要折中考慮。

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

(11)

其中，角頻率精度為

(12)

為了選擇合適的時間間隔，2.2節(jié)中描述了實驗中所選時間間隔與測量精度的關(guān)系。

在本文第2節(jié)中，將以實驗數(shù)據(jù)為例對圖1的算法流程進行詳細的說明。

2 實驗驗證

2.1實驗裝置

實驗于2013年3月在解放軍理工大學(xué)分層流實驗水槽中進行。該水槽的長、寬和高分別為12，1.5和1.2m。系統(tǒng)在實驗過程中會產(chǎn)生內(nèi)波，內(nèi)波是指發(fā)生在密度穩(wěn)定層化的流體內(nèi)部的一種波動。本實驗中使用鹽水和淡水產(chǎn)生密度分層，如圖2所示。內(nèi)波的傳播會在表面產(chǎn)生流場。實驗中風(fēng)速可變，內(nèi)波振幅為10cm，傳播方向與水槽長度方向平行，內(nèi)波傳播過程中遇到水槽壁會發(fā)生反射現(xiàn)象。圖2為實驗裝置示意圖。

圖中L表示有效掃描區(qū)域。通過特殊的光源以及散光板產(chǎn)生特殊的光照條件，使得水面光強與表面波斜率有近似的一一映射關(guān)系。光線經(jīng)過水面反射后進入圖像采集器件——線陣CCD相機。實驗中CCD相機距離水面的高度為44cm，觀測角度(與垂直方向夾角)為34°。

圖2 水槽內(nèi)波觀測實驗光學(xué)裝置示意圖

光學(xué)系統(tǒng)的主要參數(shù)如下：有效掃描長度：L=23cm；分辨率：0.38mm(空間分辨率)，1/300s(時間分辨率)；掃描頻率：300Hz；掃描點數(shù)：2048。

實驗中采用了6臺CCD相機進行并行數(shù)據(jù)采集,圖3是實際實驗裝置圖。

圖3 實驗裝置圖

2.2實驗結(jié)果及分析

以下給出第1節(jié)中所述的4個步驟的具體處理過程及部分結(jié)果。

(1) 幾何校正及光強斜率轉(zhuǎn)換

實驗中得到的時間-空間二維圖像數(shù)據(jù)如圖4所示。橫坐標x表示距離(單位：cm)，縱坐標t表示時間(單位：s)，灰度值大小表示經(jīng)過水面反射后進入相機的光強。每幅圖像的尺寸為300pixel×2048pixel，掃描頻率為300Hz，因此一幅圖代表1s時間內(nèi)掃描區(qū)域內(nèi)波動變化。表面波的運動在圖像中表現(xiàn)為傾斜的條紋，傾斜條紋的斜率即表示表面波的傳播速度。

CCD相機通過記錄掃描區(qū)域的光強度從而獲得掃描區(qū)域的實時波斜率，如圖5所示。CCD在空間上的掃描范圍約為23cm，每個采樣點的間隔約為0.23m/600=0.38mm,那么CCD可以記錄的表面波的波數(shù)范圍理論上為30～8200rad/m。

比較圖4與5，圖4中的傾斜條紋存在一定的幾何畸變(條紋彎曲)，經(jīng)過幾何校正后這種畸變消失(圖5)；而且圖4中光強分布沿距離方向是不均勻的，但經(jīng)過光強斜率轉(zhuǎn)換后的斜率分布比較均勻。

圖4 實驗獲得的CCD圖像數(shù)據(jù)

圖5 由圖4計算得出的波斜率分布

(2) 斜率噪聲抑制

實驗中測得的斜率s(x,t)和斜率譜S(kx,ω)分別如圖5和6所示。圖6中的黑色曲線為根據(jù)頻散關(guān)系(公式(7))得出的頻率與波數(shù)關(guān)系曲線。圖6中波頻集中區(qū)域與頻散關(guān)系曲線吻合很好，這一方面驗證了實驗所獲取數(shù)據(jù)的有效性，另一方面也表明了頻率成分主要分布于頻散關(guān)系曲線附近，因此可用圖7所示的自適應(yīng)二維頻域濾波器進行濾波處理，濾波后的斜率分布如圖8所示，與圖5相比斜率分布更加清晰。

某一距離點上波斜率隨時間的變化如圖9所示。

圖6 斜率譜

圖7 頻域二維濾波器

圖8 濾波后的斜率分布

圖9 1s內(nèi)波斜率隨時間變化

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

以波長λ=4.57cm為例，在觀測時間第28s內(nèi)相位變化如圖10所示。由于直接計算所得相位值位于[-π,π]區(qū)間，在計算相位差時需要對相位進行解卷繞。

圖10 第28s內(nèi)相位變化

圖11為某波長表面波的相位在一段時間內(nèi)的變化，可以看出在幀數(shù)據(jù)與幀數(shù)據(jù)切換時相位會有小幅度跳變，這個跳變使得相位差變化非常劇烈。

由第1節(jié)介紹的方法修正后的相位及相位差如圖12所示，與圖11相比較有明顯的改善。

(4) 表面波相速度計算

在計算得到相位數(shù)據(jù)后，可由公式(4)計算不同波長的表面波相速度。計算過程中對6臺CCD相機的結(jié)果進行了平均處理。

如第1節(jié)所述，表面波相速度的測量精度主要與角頻率的精度有關(guān)，而角頻率的精度與所選時間間隔有關(guān)。圖13為波長分別為2.08和1.52cm的表面波相速度測量精度與所選時間間隔的關(guān)系。

由圖13可以明顯看出，在選擇時間間隔Δt≥2.5s時，表面波相速度測量精度優(yōu)于1cm/s。

(a)

(b)

(a)

(b)

圖14為對波長為4.57cm的表面波采用不同的時間間隔時計算得出的相速度結(jié)果。

由圖14可以看出，間隔為2.5s的計算結(jié)果要明顯優(yōu)于間隔為0.5s的計算結(jié)果。

計算時采用時間間隔2.5s，經(jīng)過平滑濾波后得到的波長分別為2.08cm、1.52cm和1.14cm的表面波的相速度如圖15(a)、(b)和(c)所示。

由圖15及圖14(b)可以看出，速度有4次明顯的變化(圖14(b)中A、B、C、D四點)。這是因為在實驗過程中在水槽中有內(nèi)孤立波經(jīng)過探測區(qū)域，內(nèi)波的傳播引起了表面波相速度上的變化。圖15及圖14(b)中虛線表示在內(nèi)波經(jīng)過之前(前80s)表面波相速度均值。這4次表面波相速度的變化與實驗中觀察到的現(xiàn)象十分吻合，即內(nèi)波首先逆風(fēng)向經(jīng)過探測區(qū)域(100s左右)，而后經(jīng)過3次反射。逆風(fēng)向經(jīng)過探測區(qū)域時內(nèi)波傳播使得表面波相速度下降(圖14(b)中A、C)，順風(fēng)向經(jīng)過探測區(qū)域時內(nèi)波傳播使得表面波相速度增大(圖14(b)中B、D)。而且由圖15可以看出，雖然不同波長下相速度大小不同，但在內(nèi)波經(jīng)過時引起的相速度變化幅值是基本一致的，約為3～4cm/s，整體幅值變化約為6～7cm/s，這種表面波相速度上的變化就是內(nèi)波傳播引起的表面弱流場導(dǎo)致的。

(a) 波長為2.08cm

(b) 波長為1.52cm

(a) Δt=0.5s

(b) Δt=2.5s

(a) 2.08cm

(b) 1.52cm

(c) 1.14cm

圖16為不同風(fēng)速下內(nèi)波經(jīng)過探測區(qū)域前(前80s)不同波長表面波相速度對比。

可以看出，對于波長相同的表面波，風(fēng)速越大測得的相速度越大。這是由于風(fēng)速越大，風(fēng)生流越大，根據(jù)公式(1)，表面波相速度越大。同時，可以看出隨著風(fēng)速增大，相速度曲線在圖中近似于“平移”，這也印證了波的本征相速度與風(fēng)速無關(guān)，亦即與波高無關(guān)。

得到如圖15所示的表面波相速度后，將此相速度與前80s(內(nèi)波經(jīng)過前)的相速度均值相減，便可得到內(nèi)波傳播所引起的表面弱流場。由于可以計算不同波長的表面波相速度，將多個波長表面波計算得出的表面弱流場進行平均處理，得到結(jié)果如圖17所示。

圖16 不同風(fēng)速下不同波長表面波相速度(Vwind為風(fēng)速)

圖17 內(nèi)波激發(fā)的表面弱流場速度(正值表示內(nèi)波與表面波運動方向相同)

由圖17可以看出，在前80s左右，內(nèi)波尚未經(jīng)過觀測區(qū)域時，流場速度在零值上下波動且變化幅度很小。對前80s的數(shù)據(jù)進行統(tǒng)計，其精度達到0.3cm/s，足以應(yīng)用到對波-流相互作用機理的研究中。

3 結(jié) 論

提出了一種基于線陣CCD的水面流場光學(xué)測量方法，其最大的優(yōu)點是可以在風(fēng)浪環(huán)境下高精度地獲取表面弱流場信息而不受波浪振動的影響，克服了傳統(tǒng)PIV流場測量方法不能工作于風(fēng)浪環(huán)境下的缺點。應(yīng)用此方法對水槽內(nèi)波觀測實驗數(shù)據(jù)的處理結(jié)果表明，表面弱流場的測量精度達到0.3cm/s。已有學(xué)者利用線陣CCD測量表面波浪譜等信息[4-7]，而本文在此基礎(chǔ)上進一步同時測量表面彈流場速度，并且保持實驗中流場方向基本與水槽長度方向平行，因此可以同時獲得風(fēng)浪環(huán)境下的表面弱流場與波浪信息，這在波-流相互作用機理的研究中將會發(fā)揮重要作用。

所采用的實驗設(shè)備及方法也存在一定的局限性。實驗中對光照條件要求較高，這在實驗室環(huán)境中容易實現(xiàn)，但實際測量時測量條件較為惡劣，這對該方法的實用性會造成不利影響。因此今后的研究內(nèi)容主要是改進該方法使其適用于實際測量以及討論該方法應(yīng)用于復(fù)雜的隨機波浪表面測量可能面臨的問題和對測量精度造成的影響。

參考文獻:

[1]何宜軍. 成像雷達海浪成像機制[J]. 中國科學(xué): 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]王浩, 曾理江. 二維及三維流場的光學(xué)測量方法[J]. 光學(xué)技術(shù), 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.

作者簡介:

劉利(1989-)，男，山東泰安人，中國科學(xué)院電子學(xué)研究所碩士研究生。研究方向：信號與信息處理。通訊地址：北京市海淀區(qū)北四環(huán)西路19號中國科學(xué)院電子學(xué)研究所一室(100190)。E-mail：kelly07_11@163.com