基于高速攝像技術的推移質床面沙條帶結構試驗研究

王 浩，范潔鈴，鐘 強，陳 槐，陳 銘，彭國平

（1.福州大學 土木工程學院，福建 福州 350116；2．中國農(nóng)業(yè)大學 水利與土木工程學院，北京 100083；3.南京水利水電科學研究院，江蘇 南京 210029）

1 研究背景

推移質運動是水沙相互作用的結果，泥沙顆粒沖刷輸移過程中會造成河床沖淤演變、河岸變形、海岸后退等實際工程問題[1]。中低強度水流條件下，大尺度的紊流相干結構顯著，將對泥沙顆粒運動產(chǎn)生很大影響，如圖1所示的床面沙帶狀結構，同時泥沙顆粒運動過程所形成的床面形態(tài)也將對水流結構產(chǎn)生反饋作用。推移質運動與紊流結構之間的相互作用規(guī)律也一直是國內外學者重點關注和研究的對象[2-3]。

圖1 明渠試驗順直沙條帶結構

很多學者通過試驗和野外觀測[4-7]均發(fā)現(xiàn)河床床面受紊流相干結構影響，會形成沿水流方向順直的沙條帶結構，并開始對紊流相干結構和推移質運動及其相互作用關系開展研究。Kline 等[8]首先通過氫氣泡示蹤法，開創(chuàng)性的發(fā)現(xiàn)近壁區(qū)黏性底層中存在相互間隔的高低速條帶結構和猝發(fā)現(xiàn)象，拉開了人們對紊流相干結構研究的序幕。水沙兩相運動過程中沙顆粒受水流相干結構的猝發(fā)事件影響，能夠形成穩(wěn)定的沙粒條帶結構[9]?；谖闪飨喔山Y構，Gyr 等[10]指出在壁面光滑和粗糙條件下，水流中均存在發(fā)夾渦及發(fā)夾渦群對泥沙顆粒運動影響。Colombini 等[11]觀測到沙槽和沙脊越靠近水槽中部越明顯，越靠近水槽邊壁越不明顯。王殿常[12]利用圖像處理方法開展明渠紊流近壁區(qū)稀疏沙帶狀結構性質的研究，指出沙帶狀結構的形成與紊流大尺度結構有關。但Nezu 等[13]認為紊流猝發(fā)尺度較小，不足以引起大尺度的相干結構和沙帶狀結構，邊壁效應二次流是形成大尺度渦結構的主要誘導因素。而Zhong 等[14]則認為推移質運動在床面上形成的沙壟與沙脊與紊流的Q2/Q4事件和大尺度流向渦結構有關，認為發(fā)夾渦群是紊流猝發(fā)的本質，是維持湍動能的基礎。上述可見由于紊流結構的復雜性和對沙條帶結構認識不夠深入，前人研究總體偏向于現(xiàn)象的觀察和機理的猜測，缺乏在顆粒尺度下定量研究泥沙顆粒運動規(guī)律及其與床面形態(tài)的關系。隨著高速攝像技術、數(shù)字圖像處理、運動識別及三維地形重構等人工智能技術的蓬勃發(fā)展，為研究推移質運動規(guī)律及床面形態(tài)提供了技術支持。目前，如何將先進的測量技術應用到水沙運動研究中，仍需對相關過程中的處理技術和分析方法進行不斷完善和改進。

本文擬在前人研究的基礎上，開展推移質平衡輸沙試驗，借助高速攝像技術，采用數(shù)字圖像處理技術獲取運動泥沙顆粒精確質心，分析顆粒運動速度，并從顆粒運動軌跡和區(qū)域對推移質沙顆粒運動特征進行研究。最后結合SFM 技術獲取床面三維地形結構形態(tài)，將推移質運動與床面形態(tài)進行耦合，建立紊流相干結構與推移質運動相互作用關系機理預測，為相關研究提供參考。

2 實驗設計

2.1 試驗水槽試驗在高精度明渠水槽中開展，水槽系統(tǒng)如圖2。水槽長11 m、寬25 cm、高0.4 m，變坡范圍0～1%。邊壁和底板均為透明玻璃，水槽整體安裝誤差約為±0.2 mm。水槽的啟動與關閉由電腦控制，試驗過程中可通過變頻器控制水泵轉速調節(jié)流量，流量大小可通過電磁流量計測得，電磁流量計安裝于水槽下方管道處，誤差在0.4%以內。水槽左端進水口處設3個矩形蜂窩狀整流器，使水流平穩(wěn)地流入水槽，出口設置活頁尾門，通過尾門的開閉程度控制水深。沿水流流向在水槽中心線安裝5個超聲水位計測量瞬時水深。測量段距水槽入口7 m，以保證來流流態(tài)達到均勻流，距出口4 m可消除尾門對測量段水流的擾動。定義x軸沿水流方向，y軸為垂直水流方向，z軸沿水槽橫向方向。

圖2 試驗水槽

試驗進行約30 min 后，推移質運動將形成穩(wěn)定的水流條件和床面形態(tài)并保持平衡輸沙。此時采用高速CMOS 相機（分辨率為2560×1920pixels2）對測量段進行拍攝，相機安裝于水槽中心上方，豎直向下與水槽流向方向中心線垂直。本文試驗中，高速相機拍攝過程中畫面區(qū)域覆蓋整個水槽寬度的25 cm，流向（x方向）為18.75 cm，進行滿畫幅拍攝，成像分辨率為10.2 像素/mm。連續(xù)兩幀圖像的時間間隔Δt是圖像識別技術的一個重要參數(shù)，若Δt太小，連續(xù)拍攝圖像中泥沙顆粒位移過小，識別累計誤差將增大；若Δt太大，泥沙顆粒運動幅度大，不利于泥沙顆粒質心的識別與跟蹤。當Δt＜Tm時，連續(xù)兩幀圖像將較準確追蹤床面推移質的運動[15]，其中Tm=15D/u*，為中間時間尺度[16]，指泥沙顆粒連續(xù)兩次處于靜止狀態(tài)的時間間隔。因此高速相機最小幀頻應當滿足：f＞1/Tm。依據(jù)相機最小幀頻，試驗過程連續(xù)采樣5000幀顆粒運動圖像，拍攝頻率為200 Hz。

2.2 試驗方案使用孔徑為2.5 mm和3 mm 的篩網(wǎng)篩選出平均粒徑D=2.75 mm 的天然沙，其密度r=2650 kg/m3。在水槽中鋪設3 cm 等厚沙顆粒。距水流入口處1 m 處布設可變頻加沙機，通過改變加沙機電機頻率和出沙口寬度調節(jié)加沙速率，以確保試驗過程中床面高程基本不變，實現(xiàn)平衡輸沙。距水流進口1.5 m 處沿床面橫向豎直向下放置一薄鋼尺，其高度與沙厚度一致，使沙粒均勻流向下游。加沙裝置距測量段約5.5 m，對測量段數(shù)據(jù)的獲取無影響。本文共開展10種雷諾數(shù)的恒定均勻流平衡輸沙推移質試驗，各試驗組次的水流條件見表1。

表1 推移質試驗水流條件

3 顆粒運動分析方法

3.1 運動顆粒質心識別由于床面顆粒的相似性、顆粒運動的復雜性（如顆粒偏離水流方向與其他顆粒交疊）以及隱蔽顆?；叶容^難識別等問題，全面識別床面顆粒質心用以分析、判定顆粒運動較難實現(xiàn)。經(jīng)嘗試，先識別運動顆粒，再尋找運動顆粒質心具有很好效果。圖3（a）為原始圖像，由于圖像在拍攝過程中，存在光線照亮不均勻，對原始圖像進行頂帽變換，以校正不均勻光照的影響，同時用直方圖均衡化增強圖像的局部對比度，結果如圖3（b）所示。

圖3 運動顆粒識別

推移質顆粒在床面滑動、翻滾過程中，其不斷變化的不規(guī)則感光面上存在光線的反射差異。當高速相機連續(xù)拍攝兩幀圖像的時間間隔為Δt時，兩張圖像灰度相減所得的灰度變化，可認定為在Δt時間內推移質顆粒的運動區(qū)域[17-18]。如圖3（c）所示，灰度相減所得運動區(qū)域包含兩幀圖像中顆粒運動的信息，即包含運動前和運動后顆粒區(qū)域之和。使用高斯濾波對圖3（c）中的噪聲進行平滑處理，獲取精確的運動顆粒區(qū)域圖像如圖3（d）所示。

獲得兩幀顆粒運動區(qū)域圖像后，對區(qū)域進行質心提取，將此灰度相減區(qū)域的質心作為兩幀圖像所共有的粗略質心，如圖4（a）中的藍色圓點所示。本文在粗略質心2倍粒徑范圍內，尋找距離粗略質心最近的泥沙顆粒質心，以識別兩幀圖像中運動顆粒的精確質心，如圖4（b）中的紅色圓點所示。圖4（c）為精確質心與粗略質心的對比圖，從圖中看出紅色圓點比藍色圓點更精準的定位泥沙顆粒中心位置，表明這一改進方法能精確識別運動顆粒的質心。

圖4 運動顆粒精確質心識別識別

3.2 顆粒運動速度顆粒運動速度是推移質運動規(guī)律的重要內容之一，本文采用PTV（Particle Track?ing Velocimetry）方法對連續(xù)兩幀圖像進行相關運算，分析顆粒運動速度。本試驗中兩幀圖像拍攝時間間隔為0.005 s，假定顆粒從前一幀圖像所在位置運動到后一幀圖像所在位置做直線運動，其線性距離與時間的比值為顆粒運動速度。具體步驟如下：（1）獲得連續(xù)兩幀圖像中運動顆粒的質心，如圖5（a）（b）所示。（2）采用匹配幾率法對連續(xù)兩幀圖像的運動顆粒進行匹配計算。（3）進行剔錯和插值，求出泥沙顆粒運動速度，速度矢量圖如圖5（c）所示。

圖5 顆粒運動速度計算

3.3 運動顆粒跟蹤本文采用卡爾曼濾波方法分析推移質顆粒在拍攝區(qū)域內的運動軌跡?？柭鼮V波基本不受脈沖信號影響，它能從一系列帶有不確定性的數(shù)據(jù)中算出物體在下一時刻的最優(yōu)位置值，適用于對運動物體進行跟蹤。在運用卡爾曼濾波之前，需確定預估值和測量值，以顆粒當前位置為中心，到下一幀圖像中流向上3倍粒徑，展向上1倍粒徑的范圍內作互相關運算尋找顆粒，互相關系數(shù)最大值所對應的顆粒位置即為顆粒的測量值。以上述所求的運動顆粒精確質心為中心，加上PTV 方法求出的顆粒速度，運用運動學公式求出顆粒在下一幀的位移作為預估值。最后融合預估值和測量值，得到最優(yōu)值。預估值需用卡爾曼濾波狀態(tài)預測方程求出，公式如下：

預估值算出后，需用測量值對預估值進行更新，得到泥沙顆粒當前時刻位移的最優(yōu)值，狀態(tài)更新方程為：

式中：Zk為當前時刻的測量值；為根據(jù)預估值和測量值計算出的當前時刻位置的最優(yōu)值；Hk為狀態(tài)變量到測量變量的轉換矩陣；Kk為卡爾曼增益，用于確定測量值和預估值的權重，其方程如下式：

式中：Rk為測量噪聲的協(xié)方差矩陣；Qk-1為系統(tǒng)噪聲的協(xié)方差矩陣；Pk-1為上一時刻的誤差協(xié)方差矩陣。

為使卡爾曼濾波不斷運行，狀態(tài)更新后還需更新誤差協(xié)方差矩陣，公式如下：

式中：Pk表示根據(jù)預測誤差協(xié)方差矩陣P -

k計算出的當前時刻誤差協(xié)方差矩陣。

由于卡爾曼濾波需確定初始的顆粒位移值和誤差協(xié)方差矩陣，以第一幀圖像運動顆粒的精確質心為初始的顆粒位移值，同時根據(jù)圖像的顆粒信息算出初始協(xié)方差矩陣，定義預測誤差根據(jù)方差的定義，誤差協(xié)方差矩陣其中X1為第一幀圖像中運動顆粒的真實位置，本文初始誤差協(xié)方差矩陣為

從第二幀圖像開始運用卡爾曼濾波預測顆粒位置。根據(jù)泥沙運動特性，本文的顆粒運動狀態(tài)預測方程為：

PTV 方法計算所得的泥沙顆粒的流向速度和展向速度分別為u、v。xk、yk是泥沙顆粒在流向和展向上的位移。根據(jù)上述求出的運動顆粒精確質心，利用式（6）、式（7），可算出顆粒在下一幀圖像中的位移預估值。本文Q、R的取值均為2[19]。由于僅預測軌跡，不涉及其他參數(shù)，故轉換矩陣H取為對運動顆粒圖像重復上述步驟，最終獲得所有運動顆粒不同時刻的位置，實現(xiàn)對運動顆粒的拉格朗日概念的連續(xù)跟蹤。

4 結果與分析

圖6 運動顆粒鏈軌跡

4.1 顆粒運動軌跡當水流強度高于臨界水流強度時，泥沙顆粒將在河床床面上輸移運動，其運動軌跡是本文研究的重要方面。運用卡爾曼濾波方法獲得推移質顆粒運動軌跡后，依次列出由不同幀數(shù)圖像計算得到的運動顆粒軌跡累積圖，如圖6所示。其中同一顏色表征泥沙運動在相機拍攝范圍內的連續(xù)運動軌跡。100幀圖像的顆粒運動軌跡較為零散，隨著圖像幀數(shù)的增加，沿水流方向顆粒運動軌跡不斷疊加形成明顯的推移質運動軌跡條帶狀結構。如200幀圖像時顆粒的運動軌跡稀疏與密集相間分布已經(jīng)較為明顯，形狀與沿水流方向高低速條帶結構類似[20]，250幀顆粒運動軌跡更加稠密。

4.2 推移質運動時均特征為進一步研究時間尺度上推移質泥沙顆粒運動特征，在Δt時間內，以床面運動顆粒精確質心為圓心，泥沙粒徑為直徑的圓形區(qū)域，可認定為運動泥沙顆粒所占區(qū)域。對床面運動顆粒所占區(qū)域賦值為1，則靜止顆粒所占區(qū)域賦值為0。將5000幀圖像計算所得的顆粒鏈中運動顆粒所占區(qū)域進行累加，累加值可看作在連續(xù)拍攝時間內，顆粒運動在時間尺度上的累加。圖7即為顆粒鏈累加圖，顯示了四種不同水流條件下床面推移質運動強度等值線圖，其中b代表水槽中部到兩邊的距離，B為水槽寬度，L為水流方向長度。用各水流條件下的5000幀圖像中顆粒運動次數(shù)最大值歸一化床面泥沙顆粒運動次數(shù)，表征無量綱的推移質運動強度。從圖7來看，在水槽邊壁推移質運動強度較弱或幾乎沒有。當水流強度較弱時（Θ=0.04），推移質運動規(guī)律性不明顯；隨著水流強度（Θ=0.048，0.06）增大，在水槽中部位置，顆粒運動強度沿橫向呈有規(guī)律性的高低交替出現(xiàn)；當水流強度繼續(xù)增大時（Θ=0.080），推移質運動強度的條帶結構特征變得模糊。

圖7 顆粒鏈累加

本文認為，紊流相干結構在較低水流強度條件下，其大尺度相干結構較弱，不足以影響泥沙顆粒形成規(guī)則的形態(tài)。隨著水流強度增大，紊流中大尺度相干結構發(fā)展充分，其尺度為水深尺度，能夠影響泥沙顆粒運動形成特定規(guī)律。而當水流強度進一步增大時，紊動劇烈導致相干結構不明顯，規(guī)則的泥沙運動也隨之消失。

4.3 推移質運動與床面地形關系試驗沖刷后，床面沖刷地形同樣是分析推移質運動規(guī)律的重要方面。使用基于SFM（Structure From Motion）方法的三維重構來獲取試驗中C1345（Θ=0.060）的三維床面地形[21]。圖8（a）顯示了C1345 沖刷后重構的水槽床面三維地形，可以看出本試驗條件下推移質床面不是平坦的床面地形。這反映出紊流相干結構對推移質運動的影響，最終在床面上呈現(xiàn)出沿水槽橫向高低起伏的沙壟和沙脊。圖8（b）顯示高速相機拍攝推移質顆粒運動區(qū)域的床面三維地形。與推移質顆粒運動強度規(guī)律一致，沙脊和沙壟在水槽中部較明顯，從水槽中部向兩側逐漸模糊。

圖8 三維地形云圖

可見，推移質運動強度和床面結構均出現(xiàn)沿水槽橫向高低間隔交替的波動規(guī)律。為深入研究兩者規(guī)律，圖9（a）將C1345 同一床面位置上，推移質運動強度分布圖與三維床面地形圖進行重疊比較，圖中可以看出推移質運動概率較大的區(qū)域對應著床面的凹槽，推移質運動概率較小的區(qū)域對應于床面的凸槽。沿展向方向歸一化推移質運動強度和床面地形高程，同樣推移質運動強度較大時對應于床面地形低洼處，床面地形較高處對應于推移質運動強度較小處，如圖9（b）所示。

5 討論

以往對于床面條帶結構的形成有兩種解釋，分別是二次流模型和流向渦模型。二次流發(fā)源于水槽邊壁，在邊壁具有最大尺度，但試驗中推移質運動多在水槽中部附近運動，二次流無法解釋這一現(xiàn)象。而由發(fā)夾渦群產(chǎn)生的Q2/Q4 猝發(fā)事件所形成的大尺度流向渦模型能夠更合理地解釋上述現(xiàn)象。根據(jù)象限分析法[22]，猝發(fā)現(xiàn)象包含Q2/Q4事件，低速帶舉升的脈動流速為u＜0，v＞0，位于uv坐標系的第二象限，稱為Q2 事件；高速流體掃蕩的脈動流速u＞0，v＜0，位于uv 坐標系的第四象限，稱為Q4事件。Imamoto 等[23]在1986年提出流向渦模型，高低速條帶結構和一對運動方向相反的相鄰流向渦的存在，很好地解釋了猝發(fā)現(xiàn)象中Q2/Q4事件的維持機理以及推移質床面凹凸形態(tài)間隔排列的現(xiàn)象。大尺度流向渦向上和向下旋轉分別對應著Q2/Q4事件，Q2/Q4事件使得水流發(fā)生變化，在水面上出現(xiàn)相間排列的高低速水流條帶結構。高速帶區(qū)域的流體速度大，向下清掃過程中帶走泥沙顆粒，形成凹形沙結構，低速帶區(qū)域流體速度小，泥沙顆粒運動較少，形成凸形沙結構。如圖10所示。

圖9 床面地形和推移質運動強度對比圖

圖10 條帶結構與流向渦的概念模型

推移質運動特征和最終床面結構形態(tài)均受紊流大尺度結構的影響。而大尺度流向渦是一種瞬時結構，其尺度、位置均存在紊動隨機性，時間平均后瞬時的相干結構也將消失。但本試驗條件下，床面形成了穩(wěn)定的時間平均尺度下的沙脊沙壟帶狀結構，這一固定結構說明紊流相干結構與河床形態(tài)存在密切的相互作用、影響和制約的關系。紊流相干結構促使推移質運動形成沙脊沙壟結構，形成的沙條帶結構將反作用于紊流相干結構，使原本在紊流中紊動變化的流向渦被沙脊沙壟結構所固定，兩者相互耦合。而對于本文工況中的流場特性，仍需進一步研究。

6 結論

本文在明渠水槽中開展中低水流強度推移質平衡輸沙實驗，借助高速攝像和分析方法獲取運動顆粒質心，引入卡爾曼濾波方法實現(xiàn)對泥沙顆粒運動的軌跡跟蹤。在此基礎上，建立顆粒運動強度與床面地形結構的耦合關系，從顆粒尺度探討分析泥沙顆粒運動與紊流相干結構的相互作用機理。具體結論如下：（1）本文方法較好地實現(xiàn)在顆粒尺度下對推移質顆粒運動的識別，為定量研究推移質運動規(guī)律提供參考。（2）顆粒鏈累加所得推移質運動強度分布特征表明，中低水流強度條件下推移質運動強度在床面上是非均勻化的，顆粒運動強度在床面沿水槽橫向形成高低間隔的帶狀結構。從總體趨勢上看，水槽中部推移質運動強度最大，至水槽兩側強度逐漸減弱，在水槽邊壁顆粒運動較少，這是二次流大尺度模型無法解釋的現(xiàn)象。（3）床面沖刷后沿水流方向形成交替間隔的沙脊和沙壟結構，將推移質運動強度與床面地形作比較發(fā)現(xiàn)，在沙壟區(qū)域推移質運動強度較大，在沙脊區(qū)域運動強度較小。（4）流向渦模型很好地解釋推移質運動強度條帶結構和沙脊沙壟結構的形成，流向渦使得原本平整的床面變?yōu)樯臣股硥砰g隔的結構，而床面結構也將限制流向渦的紊動性質，兩者相互影響最終形成耦合結構。