基于數(shù)字圖像處理的顆粒流厚度動態(tài)提取方法研究

吳 越，李 坤，程謙恭,2,3，王玉峰，龍艷梅，姜潤昱，宋 章，劉 毅

（1.西南交通大學(xué)地質(zhì)工程系，四川 成都 610031；2.西南交通大學(xué)高速鐵路運(yùn)營安全空間信息技術(shù)國家地方聯(lián)合工程實(shí)驗(yàn)室，四川 成都 611756；3.高鐵線路工程教育部重點(diǎn)實(shí)驗(yàn)室，四川 成都610031；4.中鐵二院工程集團(tuán)有限責(zé)任公司，四川 成都 610031）

流動型滑坡[1?2]是高山峽谷區(qū)普遍存在的一類地質(zhì)災(zāi)害，常見的流動型滑坡包括高速遠(yuǎn)程滑坡、泥石流、火山碎屑流等。長期以來，流動型滑坡的運(yùn)動機(jī)理都是地質(zhì)災(zāi)害領(lǐng)域研究的熱點(diǎn)問題[3]。由于這一類災(zāi)害的發(fā)生往往具有突發(fā)性、瞬時性，對于其運(yùn)動過程難以進(jìn)行直接觀測，通常采用物理模型試驗(yàn)的方法對其啟動-運(yùn)動-停積的全過程進(jìn)行再現(xiàn)[4?6]。隨著監(jiān)測方法不斷完善，得以從試驗(yàn)中提取更多、更全面的運(yùn)動學(xué)參數(shù)。其中，顆粒流厚度及其演化是碎屑流物理模型試驗(yàn)中重點(diǎn)關(guān)注的要素之一。顆粒流厚度的變化可以直觀反映顆粒流形態(tài)演化過程[7]，同時，顆粒流厚度也可用于計算顆粒流的其他運(yùn)動學(xué)參數(shù)[8?10]。物理模型試驗(yàn)中，顆粒流厚度定義為底板到底板以上觀察到的大量顆粒的最高點(diǎn)之間的垂直距離，不包括表面不連續(xù)、松散的顆粒[11]。

根據(jù)流體運(yùn)動特征的觀測方法，可將顆粒流厚度監(jiān)測分為拉格朗日法和歐拉法兩種。目前顆粒流試驗(yàn)中厚度的動態(tài)提取主要以歐拉法為基礎(chǔ)。Ancey[12]采用超聲波傳感器獲取一定面積內(nèi)的顆粒流厚度的平均值；Forterre等[7]采用光吸收法，通過試驗(yàn)中獲取的光照強(qiáng)度與標(biāo)定的光照強(qiáng)度曲線進(jìn)行對比間接獲取顆粒流厚度，但由于標(biāo)定過程與運(yùn)動過程中顆粒密實(shí)度不同，最終獲取的顆粒流厚度值比真實(shí)厚度值偏小；Russell等[13]和Saingier等[14]通過激光掃描方法獲取了連續(xù)高精度的顆粒流厚度，但該方法不能排除離散顆粒；Zhou等[8]、Rognon等[10]和Takagi等[15]利用激光測距傳感器測量顆粒流厚度，但該方法無法排除離散顆粒；Sanjadp等[16]和Augenstein等[17]通過點(diǎn)探針的方法獲取連續(xù)的顆粒流厚度，其厚度不受離散顆粒影響，但該方法在測量過程中探頭會對顆粒流表面形態(tài)產(chǎn)生破壞，使測量厚度偏??；Bryant等[11]通過對高速相機(jī)影像進(jìn)行手動測量獲取顆粒流厚度，結(jié)果準(zhǔn)確，但工作量大，采樣頻率較低。

隨著圖像采集設(shè)備與計算機(jī)水平的不斷提升，學(xué)科跨領(lǐng)域的實(shí)踐，數(shù)字圖像處理技術(shù)被逐漸應(yīng)用于工程地質(zhì)領(lǐng)域并展現(xiàn)出廣闊的應(yīng)用前景[18?20]。苗得雨等[21]提出了基于Matlab環(huán)境下土體SEM圖像處理方法；吳凱等[22]使用Image-Pro Plus圖像處理軟件對壓實(shí)黃土孔隙微觀結(jié)構(gòu)進(jìn)行了定量評價；張嘎等[23]在觸面試驗(yàn)中通過數(shù)字圖像技術(shù)，獲得了土顆粒的運(yùn)動過程；彭雙麒等[24]運(yùn)用PCAS圖像處理軟件獲取白格滑坡-顆粒流堆積體形態(tài)系數(shù)；White等[25]和Stanier等[26?27]基于MATLAB編制了GeoPIV程序，以此獲取顆粒流的速度場與位移場。PIV分析中顆粒流速度場的范圍變化反映了顆粒流厚度的變化，但獲取的厚度只是實(shí)際厚度的近似值。因此，可通過數(shù)字圖像處理方法進(jìn)一步處理高速相機(jī)圖像序列，以獲取精確的顆粒流厚度。通過數(shù)字圖像處理的圖像，相對于原圖像可以更直接地向讀者解釋圖像中的信息[28]，更便于機(jī)器的自動理解。

本文通過開展顆粒流斜槽試驗(yàn)，采用高速相機(jī)對顆粒流的運(yùn)動過程進(jìn)行監(jiān)測并獲取了顆粒流通過監(jiān)測斷面時的圖像序列，在此基礎(chǔ)上，通過自適應(yīng)中值濾波、二值化、圖像腐蝕和種子填充等[29]圖像處理方法對采集的圖像進(jìn)行處理分析，并通過編制計算機(jī)程序?qū)崿F(xiàn)了顆粒流厚度的動態(tài)提取，形成一套完整的基于數(shù)字圖像處理獲取顆粒流厚度的方法。該方法減少監(jiān)測傳感器，節(jié)約人力，使圖像處理技術(shù)在顆粒流物理模型試驗(yàn)中得到進(jìn)一步應(yīng)用。通過改變監(jiān)測位置，該方法可同時應(yīng)用于顆粒流運(yùn)動過程中其他形態(tài)參數(shù)的動態(tài)提取，因此具有一定的應(yīng)用前景和實(shí)踐意義。對于實(shí)際滑坡，該方法對于滑坡演化過程研究具有一定價值。

1 斜槽試驗(yàn)方案

1.1 試驗(yàn)裝置及監(jiān)測

本試驗(yàn)采用常規(guī)的顆粒流物理模型試驗(yàn)裝置（圖1）。該裝置主要分為物料區(qū)、斜槽段及堆積槽段三個區(qū)段，分別用以對應(yīng)滑坡原型的源區(qū)、運(yùn)動區(qū)及堆積區(qū)。其中，斜槽段和堆積槽段通過同軸裝置連接，可實(shí)現(xiàn)斜槽傾角調(diào)節(jié)（本試驗(yàn)中傾角為45°）；物料供給區(qū)與斜槽段通過擋板分隔，試驗(yàn)開始前，關(guān)閉擋板并將物料置于擋板之后，通過快速開啟擋板實(shí)現(xiàn)碎屑顆粒的快速釋放。模型各部分的尺寸如圖1（a）所示，其中，斜槽寬度為20 cm，堆積槽寬度為22.2 cm。模型槽兩側(cè)壁均采用可透視的有機(jī)玻璃板，以從側(cè)向監(jiān)測顆粒流運(yùn)動及堆積過程。

為實(shí)現(xiàn)顆粒流厚度的動態(tài)測量，通過側(cè)向局部監(jiān)測的方法進(jìn)行顆粒流運(yùn)動過程圖像采集。圖像序列的獲取采用AcutEye高速攝影系統(tǒng)（16M-148CXP）（幀速率1 000 fps，圖像像素1 280×800 px）。高速相機(jī)布設(shè)于擋板以下132.5 cm處，在該位置處，顆粒流接近完全展開的狀態(tài)，且不受顆粒堆積過程的影響，監(jiān)測區(qū)范圍為14 cm×10 cm，如圖1（a）所示。

圖1 試驗(yàn)裝置Fig.1 Diagram of the experimental apparatus

1.2 試驗(yàn)工況設(shè)計

本試驗(yàn)采用棱角狀石英砂顆粒作為試驗(yàn)材料，根據(jù)土工試驗(yàn)標(biāo)準(zhǔn)篩篩分孔徑，共設(shè)定5種不同的粒徑條件（GS1—GS5），如圖2所示。各粒徑范圍及不同粒徑條件下測定所得物理力學(xué)參數(shù)如表1所示；另一方面，為對比不同體積條件下顆粒流厚度的變化特征，在每種粒徑條件下設(shè)定3種體積：V1表示2 L，V2表示5 L，V3表示8 L，共開展15組（GS1V1、GS1V2、······、GS5V3）工況的顆粒流試驗(yàn)。

圖2 各粒徑影像Fig.2 Images of particle sizes

表1 試驗(yàn)工況設(shè)置及顆粒物理參數(shù)Table 1 Test condition setting and physical properties of the granular material

2 數(shù)字圖像處理方法

圖3為同一拍攝范圍內(nèi)不同時刻、不同工況顆粒流的高速相機(jī)影像，其中，圖3（a）展示了GS3V2條件下顆粒流不同位置通過監(jiān)測區(qū)的圖像。從圖中可以看出，同一顆粒流不同位置其厚度、顆粒運(yùn)動狀態(tài)均不相同，自顆粒流前緣向尾部，顆粒流厚度先增大后減小，伴隨厚度的變化，顆粒流態(tài)也呈現(xiàn)前緣和尾部較為離散、主體部分較為密集的狀態(tài)。圖3（b）—（d）為9種不同工況條件下顆粒流主體部分通過監(jiān)測區(qū)域時的高速相機(jī)圖像。從圖中可以看出，對于同一體積、不同粒徑的顆粒流，其厚度近似相等，隨著粒徑的增大，顆粒流表面顆粒逐漸呈現(xiàn)離散的狀態(tài)；另一方面，當(dāng)顆粒粒徑相同時，隨著顆粒流體積的增加，顆粒流厚度顯著增大，且表面離散顆粒逐漸減少。由此可見，顆粒流厚度的變化是影響顆粒流態(tài)及顆粒間相互作用的關(guān)鍵因素之一。

圖3 顆粒流高速相機(jī)影像Fig.3 Granular flow behavior recorded by the high-speed camera

由于拍攝范圍較小，且拍攝方向垂直于斜槽側(cè)壁，因此，通過適當(dāng)調(diào)節(jié)拍攝角度，可以使側(cè)向拍攝的顆粒流厚度反映顆粒流在y方向的最大厚度值。同時，通過在圖像中某一位置選取一個垂直于顆粒流運(yùn)動方向的縱剖面（采樣剖面），并測量顆粒流不同時刻流經(jīng)該剖面的高度值，即可得到顆粒流厚度的整體變化趨勢。由于單次試驗(yàn)中所獲取的圖像序列數(shù)量龐大，且圖像中顆粒流與背景灰度差異顯著，因此可通過編寫計算機(jī)程序來實(shí)現(xiàn)不同時刻顆粒流厚度的提取。為更加精確地獲取顆粒流主體的厚度變化情況，需采用一系列圖像處理方法對原始圖像進(jìn)行處理。常用的數(shù)字圖像方法包括圖像變換、圖像編碼壓縮、圖像增強(qiáng)和復(fù)原、圖像分割、圖像描述和圖像識別等。本文采用的具體數(shù)字圖像方法包括自適應(yīng)中值濾波、圖像二值化、圖像腐蝕和種子填充。首先通過中值濾波消除顆粒流內(nèi)部的椒鹽噪聲，再依次通過圖像二值化、圖像腐蝕及種子填充的方法生成由顆粒流邊界包圍的閉合區(qū)域，該區(qū)域內(nèi)，顆粒流上邊界定義為與下部顆粒仍保持接觸的最上部顆粒的最高點(diǎn)，因此不包括表面不連續(xù)、松散的顆粒，下邊界則定義為斜槽面與顆粒的接觸面位置處，進(jìn)而通過在采樣剖面測量上下邊界之間的像素個數(shù)來獲取顆粒流流經(jīng)該采樣剖面的厚度變化。數(shù)字圖像處理的流程如圖4所示。根據(jù)高速相機(jī)拍攝的幀速率，本試驗(yàn)顆粒流厚度的采樣頻率為1 000 Hz。顆粒流厚度的初始值為顆粒流前緣通過采樣剖面時的厚度值。

圖4 數(shù)字圖像處理流程與相對應(yīng)處理結(jié)果影像Fig.4 Flow chart of digital image processing and the corresponding processing result image

2.1 自適應(yīng)中值濾波

自適應(yīng)中值濾波可用于去除椒鹽（脈沖）噪聲，平滑非脈沖噪聲，并減少圖像邊界失真[29]。自適應(yīng)中值濾波算法是先選取一個a×b的矩形窗口Sab進(jìn)行噪聲監(jiān)測，算法中包含兩個進(jìn)程A和B，其中，Smax為矩形窗口Sab允許的最大尺寸，算法流程如圖5所示。

圖5 自適應(yīng)中值濾波流程圖Fig.5 Flow chart of adaptive median filtering

其中進(jìn)程A包含兩個算法，分別為：

進(jìn)程B包含兩個算法，分別為：

式中：Zxy—坐標(biāo)為（x,y）處的灰度值；

Zmin—矩形窗口中最小灰度值；

Zmed—矩形窗口中灰度值的中值；

Zmax—矩形窗口中最大灰度值。

自適應(yīng)中值濾波是在圖像二值化以前對原圖像進(jìn)行的預(yù)處理，以此來增強(qiáng)圖像與背景的對比度，去除由顆?？紫杜c雜色顆粒造成的噪聲干擾，從而得到可直接用于分析顆粒運(yùn)動特征量的圖像。從圖4（a）可以看出，由于顆粒之間的孔隙灰度值與背景接近，染色顆?；叶戎递^低，在圖像中顆粒流內(nèi)部會形成與背景灰度接近的暗斑，這些暗斑在圖像二值化時會在顆粒流中產(chǎn)生空洞，從而對厚度的提取造成影響，這種暗斑即可視為圖像中的椒鹽噪聲。通過自適應(yīng)中值濾波可以有效地去除椒鹽噪聲造成的干擾，并能較好地保留顆粒流的邊界特性，提高剖面高度提取算法的精確度，通過自適應(yīng)中值濾波處理后的圖像如圖4（b）所示。

2.2 圖像二值化

圖像二值化可以通過適當(dāng)?shù)拈撝捣从吵龈信d趣的圖像輪廓，使圖像輪廓特征更加明顯，且二值化后圖像灰度值只有0或255，可以起到簡化圖形的作用。圖像二值化是先選定一個灰度值作為閾值Zt，通過對比Zxy與Zt，若ZxyZt，則Zxy值替換為255。

通過圖像二值化使圖像中顆粒流與背景對比更明顯，同時二值化之前圖像灰度值分布在0~255之間。通過適當(dāng)?shù)拈撝祵D像的像素點(diǎn)的灰度值設(shè)置為0或255，閾值的大小需要根據(jù)光照強(qiáng)度進(jìn)行調(diào)整，若光照強(qiáng)度一定，可以設(shè)定一個相同的閾值。圖像二值化可以使圖像數(shù)據(jù)量大大減少，提高計算機(jī)圖像處理效率，結(jié)果如圖4（c）所示。

2.3 圖像腐蝕

圖像腐蝕是數(shù)學(xué)形態(tài)學(xué)中的一種最基本的形態(tài)學(xué)算子[29]，可以快速將圖形邊界灰度值處理為背景灰度值。圖像腐蝕技術(shù)可以通過式（5）表示。

式中：z—所有點(diǎn)的集合；

DΘE—E對D的腐蝕；

Dc—D的補(bǔ)集；

?—空集。

通過圖像腐蝕技術(shù)可以快速提取圖像中顆粒流的輪廓，具體操作步驟為：選定適當(dāng)大小的結(jié)構(gòu)元素，對二值化圖像進(jìn)行圖像腐蝕處理，得到失去邊界的圖像，通過原二值化圖像P1與經(jīng)過圖像腐蝕圖像P2做差得到圖像P3。圖像腐蝕掉的邊界圖像P3即為顆粒流輪廓圖像，如圖4（d）所示。

2.4 種子填充

種子填充又稱邊界填充算法，是數(shù)字圖像處理方法中常用的填充算法?；舅枷胧菑姆忾]區(qū)間內(nèi)部一點(diǎn)開始，由內(nèi)逐步向外在各個方向進(jìn)行搜索，對搜索區(qū)域進(jìn)行判斷，不是邊界則給定填充顏色，若是邊界則停止填充。

基于上一步圖像腐蝕獲取的圖像中顆粒流閉合輪廓，選定顆粒流取樣剖面位于底板上方輪廓內(nèi)一點(diǎn)作為種子點(diǎn)，給定填充顏色灰度值，通過種子填充算法對各像素點(diǎn)逐個填充，直到抵達(dá)顆粒流輪廓邊界時停止填充顏色，從而實(shí)現(xiàn)對顆粒流主體部分填充的效果，使顆粒流主體部分與背景分離，種子填充后所得圖像如圖4（e）所示。

3 結(jié)果與討論

3.1 試驗(yàn)結(jié)果分析

根據(jù)以上方法，基于MATLAB平臺編寫了顆粒流厚度動態(tài)提取程序，進(jìn)行不同試驗(yàn)條件下顆粒流厚度的自動提取，得到的厚度（h）隨時間（t）的演化曲線如圖6所示。從圖中可以看出，從顆粒流前緣開始h迅速增大，然后達(dá)到一個厚度變化較小的相對穩(wěn)定狀態(tài)，之后h逐漸減小直至到達(dá)一個臨界值，隨后h開始第二次上升之后逐漸減小到零，表明顆粒流已完全通過監(jiān)測區(qū)域。另一方面，從圖6可以看出，顆粒流厚度減小到臨界值以后，相鄰時間間隔的厚度值波動幅度顯著增大，說明顆粒流表面開始出現(xiàn)顯著的顆粒離散化，對應(yīng)于圖4中的顆粒流尾部區(qū)段。而在臨界值以前，顆粒流厚度變化較為連續(xù)，顆粒流主要呈密集狀態(tài)，因此該區(qū)段可視為顆粒流的主體部分。

通過對比圖6中不同體積條件下的厚度變化曲線，可以看出：粒徑相同時，顆粒流體積越大，整體厚度越大，顆粒流持續(xù)時間越長，主要表現(xiàn)為主體區(qū)段運(yùn)動時間的延長，而離散段持續(xù)時間近似相等。顆粒體積相同時，不同粒徑顆粒流主體部分持續(xù)時間較為接近，且厚度演化趨勢近于一致，而對于粒徑較大的顆粒，離散段持續(xù)時間顯著減小。此外，隨著粒徑的增大和體積的減小，顆粒流主體部分相鄰時間間隔的厚度波動幅度逐漸增大，表明顆粒流表面更加離散，該變化趨勢與圖3（b）—（d）中所反映的不同工況條件下的顆粒流流態(tài)變化趨勢一致。

圖6 各工況顆粒流厚度曲線Fig.6 Variations in granular flow thickness with time for the tests

3.2 誤差分析

為分析本文中所提出的數(shù)字圖像處理方法的合理性與準(zhǔn)確性，通過將上述結(jié)果對比人工測讀的顆粒流厚度值對該方法所產(chǎn)生的誤差進(jìn)行分析。針對每一組試驗(yàn)所采集的圖像序列，每隔10張選取一張圖像進(jìn)行厚度的人工讀取。人工讀取的位置選取與圖像處理中的同一列像素，以圖像中斜槽表面與顆粒接觸面為底邊界，自下而上等間距（10 px）添加輔助點(diǎn)（圖7），通過輔助點(diǎn)目測讀取像素個數(shù)，并將其換算為實(shí)際厚度。通過對比圖像處理獲取的顆粒流厚度與目測讀取獲取的顆粒流厚度，可得到各組試驗(yàn)誤差比例圖（圖8）（采用目測讀取厚度減圖像處理厚度）。從圖8中可以看出，圖像處理獲取的顆粒流厚度能夠較為準(zhǔn)確地與手動測量的顆粒流厚度相匹配，但也存在個別誤差較大的數(shù)據(jù)。還可以看出，同一體積工況下，誤差比例相近；同一粒徑工況下，體積增大，誤差比例有一定增加，但相對趨勢相同。由此可知粒徑、體積等因素對顆粒流厚度誤差產(chǎn)生的影響較小。

圖7 顆粒流厚度目測讀取Fig.7 Granular flow thickness which is manually measured

圖8 各試驗(yàn)工況下誤差比例Fig.8 Error ratio for the tests

為查明產(chǎn)生較大誤差原因，將每組試驗(yàn)圖像識別的顆粒流厚度與目測讀取的顆粒流厚度進(jìn)行對比，如圖9所示（對比差異等于目測讀取厚度減去圖像處理厚度）?？梢钥闯?，圖像處理獲取的顆粒流厚度在顆粒流主體部分與實(shí)際厚度吻合度較高，而在顆粒流尾部二者的吻合程度相對較低，在某些工況下呈現(xiàn)出較大的差異。說明該方法對于密集流態(tài)的顆粒流具有良好的適用性。為分析顆粒流尾部差異較大的原因，對顆粒流尾部的原始圖像與經(jīng)過一系列處理后的圖像進(jìn)行了對比分析，結(jié)果如圖10所示。可以看出，在對顆粒流進(jìn)行側(cè)向拍攝時，離散段的顆粒流顆粒在圖像中只能呈現(xiàn)二維形態(tài)。在圖像處理過程中，由于圖像處理處于二維平面，而實(shí)際顆粒流斜槽試驗(yàn)是在三維空間中進(jìn)行，三維空間降維到二維平面，本處于不同空間位置的部分離散顆粒降維之后在圖像上呈現(xiàn)出重疊的現(xiàn)象，造成顆粒連續(xù)的視覺假象，圖像處理由此將不連續(xù)的顆粒處理為連續(xù)的顆粒，從而導(dǎo)致種子填充區(qū)域遠(yuǎn)大于實(shí)際顆粒區(qū)域（圖10），使獲取的顆粒流厚度大于實(shí)際高度。另外，結(jié)合圖8中可以發(fā)現(xiàn)，存在較大計算誤差的均為目測讀取小于圖像處理，由此可推斷較大誤差的產(chǎn)生均為圖像處理時識別的顆粒流邊界過大，從而導(dǎo)致圖像處理獲取的顆粒流厚度偏大。

圖9 數(shù)字圖像處理厚度曲線與目測讀取厚度曲線對比Fig.9 Comparison between the digital image processing thickness curve and the manual measurement thickness curve

圖10 圖像處理誤差圖像Fig.10 Image processing error image

4 結(jié)論與展望

（1）本文提出了一種基于高速相機(jī)獲取顆粒流剖面圖像，并通過自適應(yīng)中值濾波、圖像二值化、圖像腐蝕及種子填充等一系列圖像處理方式獲取動態(tài)顆粒流厚度的方法。

（2）通過顆粒流的厚度剖面，可以清晰地看出厚度剖面整體趨勢為先增高后降低，再升高之后降低，能夠反映出顆粒流的流態(tài)化形態(tài)，且后一次增高時對比圖片發(fā)現(xiàn)飛濺顆粒明顯，由此可以推斷后一次增高是由于顆粒流離散段進(jìn)入了拍攝區(qū)域。

（3）利用圖像處理獲取的顆粒流厚度，能夠準(zhǔn)確地反映顆粒流密實(shí)段的厚度，能夠較準(zhǔn)確地反映顆粒流離散段的厚度，實(shí)現(xiàn)了顆粒流厚度獲取的自動化，提高了效率，降低了勞動強(qiáng)度，精度可靠，為快速獲取大量顆粒流厚度提供了新的方法。

基于數(shù)字圖像處理的方式可以獲取動態(tài)顆粒流厚度，由此類推，可進(jìn)一步將該方法推廣用于顆粒流試驗(yàn)形態(tài)學(xué)參數(shù)獲取，如堆積體的堆積輪廓、堆積面積及堆積體擴(kuò)散速率的計算等。同時，在實(shí)際滑坡中，該方法可以用于滑坡演化過程影像的處理，通過處理結(jié)果獲取滑坡前緣擴(kuò)散速率、覆蓋范圍等運(yùn)動特征參數(shù)。