基于數字圖像處理的顆粒流厚度動態提取方法研究

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

（1.西南交通大學地質工程系，四川 成都 610031；2.西南交通大學高速鐵路運營安全空間信息技術國家地方聯合工程實驗室，四川 成都 611756；3.高鐵線路工程教育部重點實驗室，四川 成都610031；4.中鐵二院工程集團有限責任公司，四川 成都 610031）

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

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

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

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

1 斜槽試驗方案

1.1 試驗裝置及監測

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

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

圖1 試驗裝置Fig.1 Diagram of the experimental apparatus

1.2 試驗工況設計

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

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

表1 試驗工況設置及顆粒物理參數Table 1 Test condition setting and physical properties of the granular material

2 數字圖像處理方法

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

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

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

圖4 數字圖像處理流程與相對應處理結果影像Fig.4 Flow chart of digital image processing and the corresponding processing result image

2.1 自適應中值濾波

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

圖5 自適應中值濾波流程圖Fig.5 Flow chart of adaptive median filtering

其中進程A包含兩個算法，分別為：

進程B包含兩個算法，分別為：

式中：Zxy—坐標為（x,y）處的灰度值；

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

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

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

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

2.2 圖像二值化

圖像二值化可以通過適當的閾值反映出感興趣的圖像輪廓，使圖像輪廓特征更加明顯，且二值化后圖像灰度值只有0或255，可以起到簡化圖形的作用。圖像二值化是先選定一個灰度值作為閾值Zt，通過對比Zxy與Zt，若ZxyZt，則Zxy值替換為255。

通過圖像二值化使圖像中顆粒流與背景對比更明顯，同時二值化之前圖像灰度值分布在0~255之間。通過適當的閾值將圖像的像素點的灰度值設置為0或255，閾值的大小需要根據光照強度進行調整，若光照強度一定，可以設定一個相同的閾值。圖像二值化可以使圖像數據量大大減少，提高計算機圖像處理效率，結果如圖4（c）所示。

2.3 圖像腐蝕

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

式中：z—所有點的集合；

DΘE—E對D的腐蝕；

Dc—D的補集；

?—空集。

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

2.4 種子填充

種子填充又稱邊界填充算法，是數字圖像處理方法中常用的填充算法。基本思想是從封閉區間內部一點開始，由內逐步向外在各個方向進行搜索，對搜索區域進行判斷，不是邊界則給定填充顏色，若是邊界則停止填充。

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

3 結果與討論

3.1 試驗結果分析

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

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

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

3.2 誤差分析

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

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

圖8 各試驗工況下誤差比例Fig.8 Error ratio for the tests

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

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

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

4 結論與展望

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

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

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

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