陶子寅， 李 然， 陳 泉， 修文正， 華云松， 楊 暉

（上海理工大學(xué) 光電信息與計算機工程學(xué)院，上海 200093）

顆粒物質(zhì)在自然界、日常生活及生產(chǎn)技術(shù)中普遍存在，是地球上存在最多、最為人們所熟悉的物質(zhì)類型之一[1]。與流體運動相比，顆粒的運動非常復(fù)雜，兼具流體和固體的某些特征，但與它們又都不相同，顆粒流與流體重要的區(qū)別之一是粒料中存在與流體不同的混合與分離現(xiàn)象，特別是不同性質(zhì)顆粒的分離和分層[2]。由于不同顆粒物質(zhì)的混合與分聚對于社會生產(chǎn)中的諸如糧食干燥和儲存、藥品加工等生產(chǎn)工藝有著重要的影響，因此對于滾筒內(nèi)顆粒的混合與分選過程的研究具有很重要的理論和實際意義：合理地應(yīng)用顆粒的混合與分離作用能夠節(jié)省能源與加工費用，同時也能在一定程度上促進顆粒力學(xué)的發(fā)展[3]。

針對不同顆粒所組成的顆粒體系的混合與分聚現(xiàn)象，近年來國內(nèi)外學(xué)者作了許多研究。Kwapinska等[4]通過離散元模擬方法研究了旋轉(zhuǎn)速度、轉(zhuǎn)筒直徑、裝載量等參數(shù)對顆粒徑向混合的影響。趙永志等[5]利用三維三方程線性彈性–阻尼離散單元模型討論了滾筒轉(zhuǎn)速、顆粒裝載率等參數(shù)對薄滾筒內(nèi)二元S型顆粒體系分層的影響，當轉(zhuǎn)速較高時，滾筒內(nèi)形成大顆粒在外、小顆粒在內(nèi)、具有圓形界面的月亮模式，當轉(zhuǎn)速較低時形成具有波浪形界面的花瓣模式。Arntz等[6]通過玻爾茲曼方程提出了一種通過混合熵來量化顆粒體系混合特性的方式。歐陽鴻武等[7]系統(tǒng)地研究顆粒混合物體系在半填充轉(zhuǎn)鼓中的花瓣斑圖形成過程，探討規(guī)則花瓣斑圖的形成條件。

本文在上述研究內(nèi)容的基礎(chǔ)上，針對二維轉(zhuǎn)筒，創(chuàng)新性地將光流法和圖像法應(yīng)用于顆粒體系的表面流動層速度與傾斜角檢測上，同時結(jié)合混合熵這一量化顆粒體系分選程度的指標，重點探究在不同轉(zhuǎn)速下顆粒體系的混合熵、表面流動層平均速度和傾斜角三者之間的關(guān)系。

1 實驗方法

1.1 實驗裝置

圖1為本實驗中所采用的實驗裝置，由二維滾筒、可調(diào)光源、步進電機和高速面陣CCD相機組成。滾筒采用有機玻璃材料加工制成，內(nèi)徑R=140 mm，筒壁厚度為5 mm，兩個壁面之間的間隙d=10 mm。滾筒與步進電機相連，通過輸出電壓的方波頻率來驅(qū)動電機，進而控制滾筒轉(zhuǎn)速。實驗過程中滾筒的轉(zhuǎn)速范圍為0～0.12 rad/s。將轉(zhuǎn)筒與高速面陣CCD相機（3 F04 M型號）放置在同一水平面上，對轉(zhuǎn)筒側(cè)壁面進行拍照，獲得顆粒體系的正視圖。相機分辨率為 2320×1720，單位像元尺寸長度為7.4 μm，采樣頻率設(shè)置為330幀/s。轉(zhuǎn)筒內(nèi)填充的顆粒材料為兩種顏色的電鍍玻璃球形顆粒，綠色顆粒的均值粒徑d1=2.6 mm，紅色顆粒的均值粒徑d2=1.25 mm。將兩種顆粒以體積比1∶1進行填充，轉(zhuǎn)筒內(nèi)顆粒材料的填充度（顆粒體系體積與二維滾筒總體積的比值）為38%。

圖 1 二維滾筒實驗裝置圖Fig.1 Experimental device diagram for two-dimensional drum

圖 2 表面流動層的區(qū)域劃分示意圖Fig. 2 Schematic diagram of the regional division on the surface flow layer

圖2為二維滾筒內(nèi)的顆粒體系流動層的位置劃分示意圖。圖中劃分了10個區(qū)域，每個區(qū)域長度為170像素，寬度為32像素，相鄰兩個區(qū)域中心位置的間隔為200像素，并分別命名為Ⅰ～Ⅹ；以流動層的中點位置為原點o建立直角坐標系，平行于表面流動層的方向為坐標x軸，垂直于流動層的方向為坐標y軸。

1.2 流動層顆粒的速度測量算法

目前光流法大致有兩種計算方法：稠密光流和稀疏光流[7]。其區(qū)別在于稠密光流是將每個像素與速度關(guān)聯(lián)；稀疏光流則是先指定一組具有明顯特征的像素點進行匹配。因此稀疏光流相比于稠密光流不僅結(jié)果會更加穩(wěn)定，相應(yīng)的運算量也會減少許多。本文所采用的測速算法即是一種基于金字塔分層的Lucas-Kanade光流法，屬于稀疏光流這一類，相應(yīng)的特征像素點的選取采用Shi-Tomasi角點檢測算法。

角點為圖像中輪廓邊界的相交點[8]，也被定義為在任意方向的一個微小變動都會引起灰度很大變化的點?，F(xiàn)有的角點檢測算法主要有基于邊緣輪廓和基于灰度圖像兩類。其中由于基于灰度的角點檢測算法避免了在提取輪廓時存在的誤差，因而應(yīng)用更為廣泛[8]。本文所采用的是Shi-Tomasi角點檢測算法，這是一種針對Harris角點檢測[9]的改進算法，獲取Harris角點中的強角點。

Harris角點檢測算法是一種直接基于灰度圖像的角點檢測算法，其主要原理是對于強度信息為I(x,y)的圖像，在像素點 (x,y)處將一個局部的小窗口w(x,y)沿任意方向移動一段微小位移(u,v)，w(x,y)為滑動窗口函數(shù)，通過引入灰度變化函數(shù)E(x,y)來表征圖像的自相關(guān)性。Harris對于灰度變化函數(shù)E(x,y)的定義如下：

對位移后的I(x+u,y+v)進行泰勒展開，可得

O表示泰勒展開式中的高階項。忽略高階無窮小項，可將E(x,y)化簡為二次型：

對于實對稱矩陣M，通過對其對角化處理后可以得到它的兩個特征值 λ1和 λ2，由于這兩個特征值表征了圖像在x軸和y軸的表面曲率，因此可以通過這兩個特征值的大小來判別圖像中的角點、邊緣以及平坦區(qū)域。Harris在此基礎(chǔ)上給出了角點響應(yīng)函數(shù)（CRF）的定義：

式中， D et(M)表示矩陣M的行列式，Det(M)=λ1λ2； T r(M)表示矩陣M的跡， T r(M)=λ1+λ2。當目標點所求得的CRF值大于閾值時，Harris認為其為角點。Shi-Tomasi在Harris研究的基礎(chǔ)上進行了改進，只選取矩陣M的兩個特征值中較小的一個與閾值進行比較，如果大于閾值，即可以得到強角點。

圖3為通過Shi-Tomasi算法檢測到的區(qū)域角點結(jié)果。檢測到的角點將作為后續(xù)Lucas-Kanade稀疏光流法的特征點基礎(chǔ)。

圖 3 Shi-Tomasi角點檢測算法處理混合顆粒體系示意圖Fig.3 Schematic diagram of Shi-Tomasi corner detection algorithm for processing the mixed particle system

光流場的計算最初由Horn和Schunck提出[10]，其約束方程的建立基于兩個假設(shè)：相鄰圖像之間的灰度基本保持不變以及相鄰圖像之間的時間間隔很小，從而時間的變化不會引起位置的劇烈變化。假定某一時刻t圖像上某一像素點 (x,y)所對應(yīng)的灰度值為I(x,y,t)，經(jīng)過一個很小的時間間隔Δt后，運動后的像素點所對應(yīng)的灰度值變?yōu)镮(x+Δx,y+Δy,t+Δt)，根據(jù)灰度基本保持不變的假設(shè)可得

將運動后像素點處的灰度值通過泰勒展開，并忽略高階項后可得

式 中 ，Ix，Iy和It分 別 表 示 灰 度 信 號 在x，y和t上 的偏導(dǎo)。

僅有一個光流約束方程是無法確定u和v這兩個光流分量的，必須引入其他的附加約束條件。于是Lucas和Kanade在此基礎(chǔ)上引入了新的假設(shè)：假設(shè)在一個小的空間鄰域 Ω 上運動矢量保持恒定。通過空間一致性的假設(shè)，可以考慮借助領(lǐng)域內(nèi)的n個 像素點來建立n個光流約束方程用來估計光流分量。下面以一個 3 ×3的領(lǐng)域窗口為例，可以建立9個方程：

對于上式的超定方程A[u,v]T=-B，可以采用最小二乘法來求解 [u,v]T的最優(yōu)解：

由于Lucas-Kanade光流法在計算過程中假設(shè)了相鄰圖像之間的位置不會發(fā)生劇烈的變化，因而在物體運動速度較快時，該假設(shè)不成立，從而導(dǎo)致后續(xù)的計算產(chǎn)生較大偏差。為了解決上述弊端，本文采用由Jean-Yves Bouguet提出的一種基于金字塔分層的改進Lucas-Kanade光流法[11]。該方法的主要原理是將圖像的長寬每次縮放為原來的一半，通過對圖像的縮放分層來起到減小像素位移的目的，同時金字塔這種結(jié)構(gòu)也可以起到逐層分解、自上而下修正運動量的作用。

計算時先計算出圖像金字塔的頂層（即分辨率最低的圖片）的光流，將得到的結(jié)果代入相鄰兩層圖像的遞歸公式中計算出下一層圖像的猜測光流向量gL-1：

式中：dL表示像素點在第L層迭代運算中的剩余光流向量，通過標準的Lucas-Kanade光流算法迭代得到；gL表示第L層的猜測光流向量，并且假設(shè)頂層的g=[0,0]T。如此循環(huán)就能通過頂層的光流向量計算得到底層的猜測光流向量g0和剩余光流向量d0，從而得到最終的光流結(jié)果，即

式中，Lm表示金字塔的最大層數(shù)。

金字塔式的Lucas-Kanade光流法可以帶來很高的像素位移增益，從而可以在不使用大窗口的情況下處理運動幅度較大的圖片。

1.3 二元顆粒體系的傾斜角測量方法

當滾筒內(nèi)的二元顆粒體系混合達到穩(wěn)定狀態(tài)后，利用高速CCD相機抓拍二元滾筒內(nèi)的顆粒體系正視圖，并將抓拍到的圖像進行灰度化等預(yù)處理，隨后再利用邊緣輪廓提取以及曲線的擬合得到最終的表面曲線輪廓，從而進一步測量出此時堆積顆粒的傾斜角 θ[12-14]。圖4為在 ω =0.031 rad/s時通過圖像法所測得的顆粒體系傾斜角示意圖。

圖 4 滾筒顆粒體系傾斜角示意圖Fig.4 Schematic diagram of the angle of slope in the particle system

1.4 顆?；旌铣潭鹊撵?/h3>

針對顆粒物質(zhì)隨時間變化的混合或是分離程度如何量化這一問題，早期有學(xué)者們提出過一些不同的方法[15-24]。而在2010年時，Arntz等[6]提出了一種基于混合熵的方法，該方法有著統(tǒng)計學(xué)的理論基礎(chǔ)，相較于先前的其他方法也很容易應(yīng)用于不同類型的分離過程。于是借用此方法（即總體熵）來對不同轉(zhuǎn)速下的顆粒分選程度進行量化。首先在二維滾筒的正視圖上劃分1 2 ×12個網(wǎng)格單元，用黑色直線標注出來，如圖5所示。在定義網(wǎng)格之后，通過玻爾茲曼公式計算每個網(wǎng)格下的混合熵s(k)，計算公式如下：

式中，xa(k)和xb(k)分別指代在網(wǎng)格k中顆粒a所占的數(shù)量比以及顆粒b所占的數(shù)量比。在求解出每個網(wǎng)格下的混合熵后，通過加權(quán)求和的方式計算得出在某一時刻t下的總體熵S(t)，計算公式如下

式中：N表示整個二維滾筒內(nèi)的總顆粒數(shù)；n(k)表示單個網(wǎng)格k中的顆粒個數(shù)。

根據(jù)Arntz等[6]的研究結(jié)果：在同一轉(zhuǎn)速下的兩種顆粒所組成的顆粒體系，它的總體熵最終會在一個值附近波動并趨于穩(wěn)定。因此在本實驗中，通過改變轉(zhuǎn)速，選擇在不同轉(zhuǎn)速下顆粒體系達到穩(wěn)定狀態(tài)時的總體熵作為該轉(zhuǎn)速下的顆粒分選程度的量化指標。

圖 5 二維滾筒的正視圖Fig. 5 Front view of the two-dimensional drum

2 結(jié)果與討論

2.1 顆粒體系的總體熵

圖6是不同轉(zhuǎn)速下處于穩(wěn)定狀態(tài)的二維滾筒顆粒體系的總體熵的變化趨勢，以及分別選取轉(zhuǎn)速為0.023，0.065和0.102 rad/s下所拍攝的顆粒體系正視圖（分別對應(yīng)于圖 6中的（a），（b）和（c）），圖6（d）中橫軸表示二維滾筒的轉(zhuǎn)速，縱軸表示顆粒體系的總體熵。在0～0.051 rad/s的轉(zhuǎn)速范圍內(nèi)，顆粒的運動狀態(tài)表現(xiàn)為離散雪崩，顆粒體系的總體熵隨著轉(zhuǎn)速的增加逐漸減小。由圖6（a）可以看到，顆粒體系的分選狀態(tài)表現(xiàn)為花瓣模式（petals pattern），并且隨著轉(zhuǎn)速的增加逐漸轉(zhuǎn)向月亮模式（moon pattern）。0.051～0.079 rad/s的轉(zhuǎn)速范圍內(nèi)，顆粒處于由間歇性的雪崩向連續(xù)雪崩轉(zhuǎn)變的過渡狀態(tài)，顆粒體系的總體熵在定值0.614的附近波動。而在0.079～0.120 rad/s的轉(zhuǎn)速范圍內(nèi)，顆粒的運動狀態(tài)表現(xiàn)為連續(xù)雪崩，顆粒體系的總體熵隨著轉(zhuǎn)速的增加逐漸增大。值得注意的是，對于圖6的（b）和（c），滾筒內(nèi)的顆粒體系在這兩個轉(zhuǎn)速下所表現(xiàn)出的分選狀態(tài)從拍攝的圖像上很難看出差異，但從總體熵這個量化指標上看卻存在著較為明顯的差異，因此可以通過總體熵這個指標來進一步地描述顆粒分選過程中的細微變化。

2.2 表面流動層的速度特征

本實驗在測量出混合總體熵的變化趨勢后，進一步地對不同轉(zhuǎn)速下的顆粒體系的表面流動層的速度進行研究，并最終希望將顆粒速度與之前所研究的總體熵相對應(yīng)。

圖 6 二維滾筒內(nèi)顆粒體系總體熵的變化曲線Fig. 6 Variation curve of the total entropy of the particle system in the two-dimensional drum

實驗分別選取了 0.058，0.065 和 0.072 rad/s 3 個轉(zhuǎn)速，且在保證顆粒體系的分選過程達到穩(wěn)定狀態(tài)的情況下對本文所劃分的10個區(qū)域內(nèi)的速度分布進行測量研究，速度的測量方法遵循前文所介紹的Lucas-Kanade光流法，同時對于3個轉(zhuǎn)速下的區(qū)域速度的空間分布進行橫向?qū)Ρ?。圖7為在0.058，0.065和0.072 rad/s這3個轉(zhuǎn)速下，區(qū)域Ⅱ內(nèi)的速度分布圖。從圖7可以看出，當分選過程達到穩(wěn)定狀態(tài)時，顆粒體系的表面流動層上的速度呈現(xiàn)出正態(tài)分布，且標準差都在0.035左右。因而后續(xù)對速度的空間分布進行橫向?qū)Ρ葧r，選擇用每個區(qū)域內(nèi)的平均速度作為該區(qū)域位置下的流動層速度。

圖8為在3個不同的轉(zhuǎn)速下，顆粒速度沿著表面流動層的分布情況，圖中橫軸表示各區(qū)域內(nèi)的平均速度矢量距離流動層中心點的相對位置，縱軸表示表面流動層的平均速度。由圖8可以看出，顆粒速度的空間分布在不同轉(zhuǎn)速下呈現(xiàn)出近似相同的變化趨勢，靠近坐標原點的區(qū)域由于顆粒經(jīng)過加速的過程而導(dǎo)致速度較大，但總體的速度變化比較平穩(wěn)；位于表面流動層兩側(cè)靠近邊緣的區(qū)域由于顆粒在該處會發(fā)生碰撞等現(xiàn)象，顆粒之間的交換更多而導(dǎo)致速度較小，速度的變化較大；于此同時，3根曲線都在區(qū)域Ⅴ處產(chǎn)生了一個向下的突變，這與同種顆粒的速度空間分布[25]相比存在不同。通過對圖8的分析可以認為在區(qū)域Ⅴ處發(fā)生速度上的突變之前，流動層的速度變化是比較理想的，并且靠近坐標原點附近的區(qū)域相比于兩側(cè)的區(qū)域，其顆粒速度更能反映在當前轉(zhuǎn)速條件下的流動層顆粒速度的真實情況，因而在接下來的分析中選擇突變之前的區(qū)域Ⅵ作為主要的測量區(qū)域。

圖 7 區(qū)域Ⅱ內(nèi)的速度分布直方圖Fig. 7 Histogram of the velocity distribution in region II

圖 8 表面流動層的速度空間分布隨轉(zhuǎn)速的變化曲線Fig. 8 Spatial distribution of the velocity on the surface flow layer at different rotating speeds

2.3 顆粒速度、傾斜角與總體熵的關(guān)系

圖9是在不同轉(zhuǎn)速條件下，顆粒運動的平均速度與堆積顆粒傾斜角的變化曲線，圖中橫軸表示二維滾筒的轉(zhuǎn)速，左側(cè)縱軸表示區(qū)域Ⅵ下的流動層平均速度，右側(cè)縱軸表示顆粒體系的傾斜角。從圖中可以看出，隨著轉(zhuǎn)速的增加，顆粒運動速度與傾斜角的變化趨勢一致，呈現(xiàn)出一個先增大后平穩(wěn)再增大的變化過程。雪崩過程中的顆粒運動速度與傾斜角呈正比例關(guān)系，這是因為在轉(zhuǎn)筒內(nèi)顆粒流中，當堆積顆粒的傾斜角增大時，顆粒床的重心位置也抬高，使得雪崩前顆粒的重力勢能增加，導(dǎo)致崩塌過程中顆粒的運動速度變大。在0～0.051 rad/s范圍內(nèi)，轉(zhuǎn)筒內(nèi)顆粒的運動狀態(tài)為離散雪崩，顆粒的運動速度與傾斜角隨著轉(zhuǎn)速的增加呈線性增大。而在0.079～0.120 rad/s的范圍內(nèi)，轉(zhuǎn)筒內(nèi)顆粒的運動狀態(tài)為連續(xù)雪崩，此時顆粒的運動速度與傾斜角也隨著轉(zhuǎn)速的增加呈線性增大。值得注意的是，在0.051～0.079 rad/s的范圍內(nèi)，平均速度與傾斜角為定值，分別為0.085 m/s和32.85°。這可能是由于轉(zhuǎn)筒內(nèi)混合顆粒處于過渡運動狀態(tài)（過渡運動狀態(tài)是由間歇性的雪崩向連續(xù)雪崩轉(zhuǎn)變的中間狀態(tài)）導(dǎo)致的。

圖 9 顆粒運動平均速度與傾斜角隨轉(zhuǎn)速的變化曲線Fig. 9 Average speed and angle of slope of the particle motion at different rotating speeds

對比圖9與圖6發(fā)現(xiàn)，轉(zhuǎn)速大于0.051 rad/s時，顆粒的運動速度、傾斜角與顆粒體系的總體熵的變化趨勢一致，在0.051～0.079 rad/s的范圍內(nèi)運動速度、傾斜角與總體熵三者均為定值，而在0.079～0.120 rad/s的范圍內(nèi)三者均隨著轉(zhuǎn)速的增加而增大。這說明了顆粒的運動速度和傾斜角影響著顆粒體系的總體熵。另一方面，也可以根據(jù)顆粒的運動速度或堆積顆粒的傾斜角去判斷總體熵的變化，即顆粒體系的分選程度。

3 結(jié) 論

本文在滾筒轉(zhuǎn)速為0～0.120 rad/s，滾筒內(nèi)二維顆粒體系填充度為38%的條件下，采用Lucas-Kanade稀疏光流法與圖像法觀察均值粒徑分別為2.6 mm和1.25 mm的規(guī)則球形顆粒所組成的顆粒系統(tǒng)的分選過程。通過測量顆粒體系的總體熵、顆粒體系表面流動層的平均速度與傾斜角，得到如下結(jié)論：

a. 針對二維滾筒內(nèi)混合顆粒體系的分選程度，可以通過圖像法與總體熵相結(jié)合的方法進行判別，當圖像法觀察不出明顯的差異（如本實驗過程中轉(zhuǎn)速為 0.065 rad/s以及 0.102 rad/s時），可以通過測量顆粒體系的總體熵來區(qū)別分選過程中的細微變化。

b. 顆粒體系表面流動層的平均速度與傾斜角的正比例關(guān)系，說明了在崩塌過程中平均速度受傾斜角的影響。

c. 顆粒體系的總體熵、表面流動層平均速度和傾斜角與轉(zhuǎn)速之間的關(guān)系，在轉(zhuǎn)速大于0.051 rad/s時三者的變化趨勢相一致，說明在混合顆粒的運動狀態(tài)轉(zhuǎn)為連續(xù)雪崩的條件下，顆粒體系的總體熵，即分選程度與平均速度和傾斜角有關(guān)，可以借助平均速度和傾斜角來判斷總體熵的變化。