李珍?劉曉星?安百鋼

摘要：本論文工作針對移動床內顆粒物料卸料特性的CCD攝像圖片，開發(fā)了基于OpenCV的圖像識別處理程序和顆粒追蹤測速程序，以刻畫移動床卸料過程中的流型演化和顆粒運動信息。為檢驗程序的準確性，以離散單元法數(shù)值模擬得到的顆粒圖像為對象，對比了所開發(fā)圖像后處理程序分析結果與離散單元法數(shù)值模擬數(shù)據(jù)，發(fā)現(xiàn)兩者能定量吻合。

關鍵詞：OpenCV；圖像識別；信息處理；顆粒流

一、前言

本研究工作的研究目的之一是基于CCD高速攝像技術來研究移動床卸料過程中床層內顆粒物料的流動特征，進而得到流動區(qū)輪廓以及顆粒速度隨時間的變化。針對這一研究目的，開發(fā)了圖片后處理程序來對圖片進行處理，以得到相關顆粒運動信息。程序開發(fā)引用了OpenCV視覺庫里的API，基于粒子圖像測速法（PIV）的思想，監(jiān)測流場內所有顆粒的瞬時運動狀態(tài)，使用最近鄰法算法來匹配批量圖像上的同一個顆粒，實現(xiàn)對移動床內顆粒物料的宏觀流動區(qū)的演化以及微觀單個顆粒的運動狀態(tài)的監(jiān)測。為此編寫了兩套程序，一為流動區(qū)識別程序，二為顆粒識別追蹤測速程序。

二、實驗裝置

本論文工作所開展實驗的主要目的是探究移動床厚度對床層內顆粒卸料行為包括流動區(qū)輪廓、顆粒速度分布等參數(shù)的影響。參考文獻中所報道的移動床實驗裝置，結合工業(yè)項目需求，搭建了移動床，安裝了CCD高速攝像儀、LED燈、電腦等圖像拍攝和記錄裝備，用于研究顆粒物料卸料過程。實驗中通過CCD高速攝像儀拍攝一系列不同卸料時刻床層物料照片，然后用自主開發(fā)的流動區(qū)識別程序以及顆粒追蹤測速程序來處理這些照片，監(jiān)測提取卸料過程中床層物料內流動區(qū)輪廓的變化以及單個顆粒的運動狀態(tài)等信息，以探究床層厚度對顆粒物料卸料特性的影響。

設計搭建的移動床卸料實驗平臺如圖1所示。移動床裝置被固定在鋁合金支架上，移動床壁面的材質是亞克力板，鋁合金支架下方預留自由空間以接取卸料后的顆粒。高速攝像儀鏡頭用的是尼康AF Nikkor 50mm f/1.4D，最高幀頻可達每秒120000幀，燈源為Godox的LED308，電腦是聯(lián)想品牌電腦，操作系統(tǒng)使用的是Windows10操作系統(tǒng)。

三、程序開發(fā)

（一）流動區(qū)識別程序

所編寫的流動區(qū)識別程序實現(xiàn)邏輯如下：用CCD高速攝像儀拍攝的移動床卸料的原始圖片為灰度圖像，顏色值為單通道，在0到255之間。imread讀取相鄰的兩張圖片，灰度值做差之后result_img得到新的灰度圖像[1]。如果前后兩張圖片上同一位置的灰度值相同，代表兩張圖片中該位置處顆粒未發(fā)生運動或者改變。由于兩張圖片同一位置處灰度值無差別，在生成的新圖片上，這個位置的灰度值相減變?yōu)?，用黑色表示，代表靜止區(qū)。如果前后兩張圖片上同一位置的灰度值不同，則代表該位置處顆粒發(fā)生改變或者產(chǎn)生運動。新生成圖片的灰度值為兩張圖片的灰度值之差，用非黑色，即不同顏色的灰色（1到255）表示，代表流動區(qū)。

因為預處理圖片較為模糊，且有噪聲，所以采用threshold調節(jié)閾值，使圖片變?yōu)楹诎讏D像。實際操作中選定一個閾值，灰度值大于這一閾值時，顏色值設置為1，是白色，反之為0，即黑色。經(jīng)過測試，閾值為33時，流動區(qū)整體輪廓最為清晰。選取前后兩張顏色值做差的圖片的時間間隔也很重要。如果時間間隔取得比較密，例如0.01s，最終流動區(qū)與后文分析的速度場相匹配，可以看到最明顯的速度波現(xiàn)象；如果時間間隔取得稍微大一些，例如0.05s，則看不到速度波現(xiàn)象，而是整體輪廓最為清晰的流動區(qū)。獲得流動區(qū)之后，為進一步監(jiān)測統(tǒng)計流動區(qū)的具體演化，需要對流動區(qū)的輪廓進行提取。

圖2是提取顆粒流動區(qū)整體輪廓的示意圖。獲得流動區(qū)輪廓坐標的程序實現(xiàn)方法如下：首先調用OpenCV的庫函數(shù)morphologyEx，對圖2中左圖的流動區(qū)做開運算處理，使得流動區(qū)內部的空隙以及孔洞完全被填充，成為一個整體，進而得到如圖2中間圖片所示的白色區(qū)域。形態(tài)學運算中有腐蝕，膨脹，開運算和閉運算幾種常見的方法，其中腐蝕方法可以抹除輪廓邊界像素點，消除毛邊，毛刺，以及一些無意義的擾動點。對輪廓邊界使用腐蝕，會使得輪廓邊界減小，向內收縮。膨脹方法是將輪廓所臨近相接的所有像素點吸收到該輪廓里，使得輪廓邊界向外延展。它可以使得輪廓邊界內部的孔洞斑點得以填充。將這兩種方法結合，可以先使用腐蝕后使用膨脹，稱為開運算，它的功能是去掉離散的小斑點、分開接近的兩個斑點、消除大輪廓邊界的毛糙點并且盡可能保證輪廓面積不變。先使用膨脹后使用腐蝕，稱為閉運算，它的功能是去掉大輪廓內部的小斑點、連接靠近的幾個斑點、消除大輪廓邊界的毛糙點并且盡可能保證輪廓面積不變[2]。本論文工作采用開運算策略。得到圖2中的中間圖片后，調用OpenCV的庫函數(shù)findContours，對該圖片中的白色區(qū)域進行輪廓識別，可以得到大致的輪廓坐標信息，如圖2中右圖的紅色閉合曲線。

得到流動區(qū)輪廓坐標后，可以進一步對流動區(qū)的一些幾何信息進行提取，來具體描述流動區(qū)的演化。圖3是提取流動區(qū)寬度、高度等信息的示意圖。獲得自由表面高度h、流動區(qū)高度H、流動區(qū)寬度W、特征區(qū)寬度W（1/2h）的程序實現(xiàn)思路如下：遍歷輪廓縱坐標，當輪廓縱坐標Y1最小時，得到此時的橫坐標X，即為卸料口的橫坐標，接著尋找橫坐標為X時，縱坐標的最大值Y2，Y2減去Y1得到H。因為流動區(qū)從卸料口處發(fā)展至顆粒初始堆積高度（自由表面處）時流動區(qū)高度會有一個先增高的變化，之后隨著自由表面高度的降低而降低，所以自由表面高度h的變化是流動區(qū)高度H發(fā)展到顆粒初始堆積高度后開始降低的部分。得到流動區(qū)高度后，從卸料口位置往上檢索縱坐標相同的輪廓坐標，它們的橫坐標相減即可得到流動區(qū)的寬度W，當輪廓縱坐標的值為h的一半時，可得到特征區(qū)寬度W（1/2h）。

流動區(qū)輪廓信息提取的關鍵代碼如下：

for （int j = 0; j < contours[row].size（）; j++）

{

if （contours[row][j].y > maxY）

{

maxY = contours[row][j].y;

X = contours[row][j].x;

}

}

cout << "x = " << X? << endl;

for （int j = 0; j < contours[row].size（）; j++）

{

if （X == contours[row][j].x && contours[row][j].y < minY）

minY = contours[row][j].y;

}

h = maxY - minY;

W = 0.5*h + minY;

獲得流動區(qū)寬度、高度等的信息后，有一個拐點值得注意，就是靜止區(qū)與移動床邊壁交界的最高點，即有效轉捩點（ETP）。ETP點的高度會隨著卸料的進行發(fā)生改變。根據(jù)流動區(qū)輪廓坐標提取獲得ETP點信息的邏輯：確定流動區(qū)輪廓的最大寬度widthMAX，遍歷不同縱坐標下的流動區(qū)寬度W，當W等于widthMAX時，并且此時的縱坐標最小，縱坐標的高度即為ETP點的位置，ETP點是靜止區(qū)與床層交界處的最高點的位置。

從而得到ETP點的位置隨自由表面高度h的變化曲線，曲線的變化趨勢分為三個階段：（1）當自由表面的高度隨著卸料的進行從最高處開始降低，有效轉捩點縱坐標位置y處于相對平穩(wěn)期，在某一值附近波動，對坐標y進行累計平均得到ETP點；（2）一段時間后卸料流型由半整體流型轉捩為漏斗流，有效轉捩點縱坐標y開始急劇增大，此處即為流型轉捩點，此時靜止區(qū)的最高點會到達自由表面；（3）ETP點到達自由表面之后，此時它的高度與流動區(qū)整體高度一致，它的縱坐標y隨著自由表面高度h的降低而降低，直至卸料完成。至此，流動區(qū)識別程序可以得到宏觀流動區(qū)演化過程中的大部分運動信息。

（二）顆粒識別程序和顆粒追蹤測速程序

通過流動區(qū)識別程序監(jiān)測得到宏觀流動區(qū)的演化信息后，為更清晰地了解卸料過程中顆粒的運動特性，進一步對微觀單個顆粒進行識別、追蹤匹配、測速等。要追蹤單個顆粒的軌跡，首先要對單個顆粒輪廓進行識別。

顆粒質心識別程序的實現(xiàn)操作如下[3]：首先基于for循環(huán)，imread讀入批量的數(shù)字圖片，再將每張圖片進行cvtColor灰度化處理以減少運算量，然后通過GaussianBlur高斯模糊減少噪聲，再進行morphologyEx形態(tài)學處理[4]，通過腐蝕膨脹減少毛刺，填充球形內部的黑洞，使球體進一步接近球形；之后進行threshold自適應閾值化操作，選取合適的閾值進一步突出顆粒的外邊界輪廓；最后對顆粒的輪廓和邊界進行識別，使用findContours，得到顆粒的坐標信息[5]。

顆粒質心信息提取的關鍵代碼如下：

for （int t = 0; t < contours.size（）; t++）

{

mu[t] = moments（contours[t]， false）;

if （mu[t].m00 == 0）

{

mc[t] = Point2f（0， 0）;

}

else {

mc[t] = Point2f（mu[t].m10 / mu[t].m00， mu[t].m01 / mu[t].m00）;

}

cout << "x： " << mc[t].x << "? " << "y： " << mc[t].y << endl;

drawContours（src， contours， t， Scalar（0， 0， 255）， 2， 8）;

}

顆粒質心識別程序識別顆粒整體輪廓的效果如圖4所示，除了極個別的顆粒無法識別，99.9%以上的顆粒都可以被識別，表明程序整體識別性能優(yōu)良，可為后續(xù)具體的移動床卸料實驗的圖像處理提供優(yōu)秀的技術支持。

同時，當CCD高速攝像儀的攝像頭需要拍攝的顆粒數(shù)目較多，分辨率不夠，拍攝得到的顆粒輪廓較為模糊時，可以通過識別單個顆粒表面反射的光斑來確定顆粒的坐標[6]。識別顆粒光斑確定顆粒坐標，統(tǒng)計得到程序準確率高達99.9%，這兩種識別方法可根據(jù)實驗不同需求進行選擇。

識別單個顆粒，得到顆粒的坐標信息之后，對批量圖片上處于不同位置處的同一個顆粒進行匹配、測速，得到顆粒的運動軌跡和速度信息[7]。顆粒追蹤測速程序的算法邏輯如下：輸入之前程序處理得到的顆粒質心坐標，選取第tn時刻的第i個顆粒的坐標與第tn+1時刻該位置附近的顆粒坐標進行比較，如果兩者之間的距離小于顆粒的半徑，則說明是同一個顆粒。這是因為選取的兩個時刻的間隔非常短，顆粒的位移不超過顆粒的半徑，這也是最近鄰法的思路。CCD高速攝像儀可以輔助高速拍攝，滿足實驗的需求。完成第tn時刻的所有顆粒的匹配，輸出結果，將已匹配顆粒之間的移動距離除以兩張圖片的時間間隔可獲得該顆粒的移動速度。自此，實驗中比較重要的顆粒運動信息都可以通過程序對圖片進行后處理提取得到[8]。

四、程序處理結果精度驗證

為驗證編寫的顆粒識別程序和顆粒追蹤測速程序的準確性，使用編寫的程序對離散單元法數(shù)值模擬得到的顆?？梢暬繄D像進行信息處理。之所以采用該策略是因為離散單元法數(shù)值模擬中不同時刻每個顆粒的信息都是已知的，因此，用自主編寫的程序識別批量模擬圖片得到的顆粒信息，與數(shù)值模擬監(jiān)測得到的顆粒信息進行對比，如若二者的誤差較小，即可證明編寫的圖像后處理程序精確度較高，可用于顆粒流的卸料實驗研究。

首先對比分析模擬得到的顆粒坐標和程序處理相同實驗條件下模擬照片得到的顆粒坐標信息，發(fā)現(xiàn)二者的坐標位置重合度較高，顆粒坐標的識別準確度高達100%，但是二者的顆粒坐標位置有細微差別，出現(xiàn)的誤差在可接受范圍之內[9]，證明顆粒識別程序準確度比較高，可為后續(xù)的顆粒追蹤測速程序進行高質量服務。

識別顆粒坐標之后，需要對單個顆粒進行匹配測速，得到不同高度位置處的顆粒速度分布曲線，從而監(jiān)測顆粒的運動軌跡和運動信息的連續(xù)變化[10]。圖5為顆粒追蹤測速程序處理批量模擬圖片得到的速度分布和離散單元法數(shù)值模擬得到的速度分布，兩條速度分布曲線基本重合，統(tǒng)計得到，顆粒追蹤測速程序的準確度大于99.9%，表明顆粒追蹤測速程序準確性良好，可以為后續(xù)更加深入的實驗研究提供較好的數(shù)值統(tǒng)計分析服務。

五、結語

本文詳述了對CCD攝像得到的大批量圖像進行數(shù)值處理的算法思想和程序實現(xiàn)策略，基于OpenCV開發(fā)了圖像識別程序和顆粒追蹤測速程序，并對所開發(fā)后處理程序顆粒識別功能的準確性和識別精度進行了檢驗。驗證結果表明所開發(fā)的圖像后處理程序能夠精確地對移動床卸料過程中運動狀態(tài)隨時間變化的顆粒進行識別和追蹤，從而為接下來具體移動床卸料實驗的數(shù)值分析奠定了扎實的基礎。本文的工作對顆粒流的數(shù)值處理技術具有一定的指導意義，并在顆粒物料相關的工業(yè)應用中具有廣闊前景。

參考文獻

[1]曹建海，路長厚.基于灰度圖像和矢量的圓心定位[J].光電子·激光，2004（06）：714-718.

[2]李冬冬，胡浩然，肖明，王小毛.基于PFC-FLAC離散-連續(xù)耦合的隧洞開挖與錨桿支護顆粒流模型[J/OL].武漢大學學報（工學版）：1-9[2023-03-27].http：//kns.cnki.net/kcms/detail/42.1675.T.20230215.1743.009.html

[3]鐘華勇，葉建生，何高清.基于OpenCV的圓心坐標定位的優(yōu)化設計[J].制造技術與機床，2021（05）：110-115.

[4]肖劍雄峰，趙漢理，吳承文.局部對比度增強的彩色圖像灰度化參數(shù)化算法[J].計算機時代，2016（03）：45-49.

[5]趙飛翔，遲世春.基于離散元的碎片尺寸隨機的顆粒破碎模擬方法[J].東北大學學報（自然科學版），2023，44（03）：408-414.

[6]Na Zhao，F(xiàn)ang Zhao，Wen-Zhen Zhong. EDEM Simulation Study of Lignite Pressing Molding with Different Particle Size[P]. Proceedings of the 2015 4th International Conference on Computer， Mechatronics， Control and Electronic Engineering，2015.

[7]Feng Wenyu，Han Yanlong，Chen Peiyu，Li Anqi，Wang Yinglong，Zhang Jincheng，F(xiàn)ei Jiaming，Hao Xianzhi，Shen Shaohang，Jia Fuguo. Effect of eccentricity on the particle flow characteristics in a double-orifice silo[J]. Powder Technology，2023，421.

[8]Wei Wei，Ding Lieyun，Luo Hanbin，Li Chen，Li Guowei. Automated bughole detection and quality performance assessment of concrete using image processing and deep convolutional neural networks[J]. Construction and Building Materials，2021，281.

[9]Sakamaki Shouta，Krengel Dominik，Mueller Jan，Chen Jian，Matuttis Hans Georg. Implementation of numerically exact Coulomb friction for discrete element simulations in three dimensions[J]. EPJ Web of Conferences，2021，249.

[10]Laurent Babout，Krzysztof Grudzien，Eric Maire，Philip J. Withers. Influence of wall roughness and packing density on stagnant zone formation during funnel flow discharge from a silo： An X-ray imaging study[J]. Chemical Engineering Science，2013，97.

作者單位：李珍，遼寧科技大學化學工程學院、中國科學院過程工程研究所多相復雜系統(tǒng)國家重點實驗室；劉曉星，中國科學院過程工程研究所多相復雜系統(tǒng)國家重點實驗室、中國科學院大學化學工程學院；安百鋼，遼寧科技大學化學工程學院