崩塌條件下顆粒重排區域的堆積結構圖像測量方法

中圖分類號：TB126文獻標志碼：A

Abstract： The triggering and propagation mechanisms of particle material avalanche constitute an exceedingly complex issue， with the crucial factor being the variation in the packing structure of the internal static accumulation area. To delve further into the packing structure of the rearrangement region within the avalanche system， this study employed an image-based approach to identify the spatial coordinates of three-dimensional cylindrical particles on the side and utilized Voronoi diagrams to characterize their packing structure. The results indicate that the Voronoi polygon area of the side-stacked particles follows a gamma distribution. Linear structures are more prone to participating in particle rearrangement， and they exhibit greater fragility in the direction perpendicular to the radius of the cylinder. During the avalanche interval， the changes in the scale of the passive layer particle packing structure are positively linearly correlated with the magnitude of the next avalanche and negatively linearly correlated with the time interval until the next avalanche occurs.

Keywords： granular flow; packing structure; Voronoi diagram; avalanche; rearrangement

引言

顆粒物質是由大量固體顆粒堆積而成的離散系統，廣泛存在于自然界與工業生產中[1-5]。顆粒物質的崩塌是其中一種常見的運動模式，常見于雪崩、泥石流等自然災害中[]。顆粒崩塌可以劃分為表面崩塌顆粒流與內部靜態堆積顆粒。內部顆?？傮w上可以認為是相對靜正的，但其中會有少數顆粒發生局部、微小的非仿射位移，導致局部堆積結構的變化，這種局部的結構上的變化也被稱為顆粒重排[8-9]。內部顆粒的局部堆積結構變化很大程度上影響著整個系統的動力學行為，研究其在崩塌前后的分布、變化特征可以更好地對崩塌行為進行預測。

目前研究人員在已知顆粒具體位置的模擬系統中已經提出了許多描述堆積結構的方法，例如傳統方法的Voronoi圖用來對堆積顆粒的體積進行劃分，新提出的Cell網格方法將堆積顆粒按照接觸顆粒的中心點沿著逆時針方向連接起來圍成的區域進行劃分，保留了接觸顆粒的連通性信息，在顆粒統計力學中發揮著重要作用[10-12]，這些方法的有效性已被廣泛證明。但在實驗系統中，如何獲得每個顆粒的位置依然是一個難題。圖像法在二維圓盤的堆積結構測量中已被廣泛使用，華云松等[13]通過連通域特征提取二維圓盤顆粒的空間坐標，并通過Cell量化了其在二維筒倉中的幾何堆積結構。Seguin等[14]使用顆粒追蹤以及光彈性技術考察了二維壓縮系統中接觸力網絡的變化。但對三維緊密堆積的顆粒系統來說，獲得其堆積結構等顆粒微觀特性的實驗手段是有限的。Brodu等[15使用激光片光照射浸沒在折射率匹配的浴液中的浸透熒光染料的水凝膠顆粒，采集顆粒系統橫截面圖像，通過一組橫截面圖像疊加重建出顆粒的表面。但該方法對顆粒材質以及顆粒系統限制較大且時間分辨率較小，難以捕捉到堆積結構的變化。Bom等[1使用散射方法，用太赫茲電磁輻射在顆粒介質中的傳輸特性量化顆粒堆積結構，但該方法不夠直觀且還處于理論階段，距離實際使用還有一段距離。

在三維顆粒系統成像困難的情況下，可以考慮用表面顆粒的二維信息表征三維系統的特征。如在顆粒流測速領域，研究人員使用粒子圖像測速（particle imagevelocimetry，PIV）方法測量表面顆粒的二維流場[17]，從而一定程度上反映了三維顆粒流的運動特征。三維滾筒中的顆粒系統沿滾筒軸向大致是均勻分布的，滾筒側面顆粒的堆積結構一定程度上也可以反映整個系統的堆積結構。但之前在PIV中對顆粒表面圖像的處理不涉及具體顆粒的識別，如何準確識別表面單層顆粒，濾除堆積孔隙中露出的內部顆粒的干擾，還需要進一步優化。

本文使用Voronoi圖研究了三維滾筒中的顆粒體系在崩塌周期中側面顆粒的堆積結構。首先使用圖像法對滾筒側面表層顆粒進行識別，使用具有顆粒反光形狀特征的特征模板去除第二層顆粒的影響；再使用Voronoi圖來表征滾筒側面表層顆粒的結構特征，以此來反映三維滾筒內部顆粒的結構；然后研究崩塌周期中重排區域堆積結構的變化，并通過圖像差分的方法估算出崩塌規模，分析其與堆積結構之間的相關性。

1實驗

旋轉滾筒內的顆粒崩塌相對容易控制與測量，并且可以持續產生周期性的崩塌流，其已成為研究顆粒崩塌流動力學的常用模型系統[18-19]。實驗中使用低速旋轉的三維滾筒為后續研究提供周期性的崩塌流，如圖1（a）所示。滾筒中顆粒崩塌系統可以劃分為表面崩塌流（主動層）以及內部靜態堆積的區域（被動層），如圖1（b）所示。顆粒速度梯度為0的區域為主被動層交界處，根據之前對該實驗系統的測量，距離自由表面垂直深度大于 10mm 的區域為被動層[20]，此處自由表面的角度取上、下休止角的均值。

本實驗中使用的面陣高速相機為千眼狼5F04M型黑白攝像機，該相機基于CMOS圖像傳感器，最大內存為 16GB 。在實驗過程中選取了兩種幀率。采集整個滾筒側面圖像時，目標主要為表面崩塌顆粒流，其運動速度較大，不需要太高的幀率。同時為了消除電機振動、轉速不穩定帶來的誤差，相機以5000幀/s的幀率采樣，每幀圖像分辨率為 1984×1680 像素，每10幀圖像作平均，最終得到的圖像序列為500幀/s。在采集內部顆粒重排事件圖像時，顆粒重排的持續時間大概在 1ms 數量級，故為進一步捕捉顆粒重排過程，相機采用10000幀/s的幀率進行采樣，圖像分辨率為 256×1396 像素。

這兩種幀率下采集得到的單幀圖像并沒有明顯差異，區別主要體現在相鄰兩幀圖像之間的差異上：高幀率下相鄰兩幀圖像差異小，便于觀察重排過程中顆粒詳細的運動軌跡；低幀率下相鄰兩幀圖像之間差異大，直接展示出重排顆粒的起始位置與終止位置，省略了重排的具體過程，便于重排事件的統計，并且低幀率下可以采集更多的連續的崩塌周期，便于統計具有相同歷史堆積結構的重排事件。在后文中，2.1節對重排細節的討論中使用的是高幀率下采集的圖像，2.2節對重排與崩塌之間相關性的討論中使用的是低幀率下采集的圖像。

1.1 實驗裝置

實驗裝置由高速相機、LED光源、三維滾筒與直流電機組成，如圖1（a）所示。滾筒由透明有機玻璃制成，直徑 140mm ，長 20mm 。兩根金屬支撐桿水平固定且相互平行，均可以旋轉，兩者之間通過一根傳動皮帶連接。旋轉電機采用直流電源供電，控制其中一根支撐桿旋轉，另一根支撐桿在傳動帶作用下以相同的角速度旋轉。支撐桿上固定4個圓柱形橡膠塊，滾筒放置于橡膠塊上，軸線與支撐桿平行，通過橡膠塊的摩擦力帶動旋轉。改變直流電源輸出電壓可以控制滾筒轉速，滾筒轉速設置為 0.19r/min 。滾筒內填充粒徑為 0.95～1.05mm 的干燥、無黏性電鍍玻璃珠，平均粒徑為 1mm 。滾筒裝置放置在水平實驗平臺上，采用高速相機沿著滾筒軸向采集滾筒側面圖像，盡可能保證相機與滾筒中心處于一條水平線上。

同時在采集過程中用LED燈照射滾筒側面進行補光，以提高采集圖像的對比度。實驗開始時，滾筒內顆粒自由堆積，自由表面處于水平狀態；在滾筒勻速旋轉 5min 后認為系統進入穩定的崩塌周期，開始圖像采集。圖1（b）給出了其中一幀典型的實驗圖像。

1.2 顆粒堆積圖像處理方法

圖2給出了5000frame/s幀率下獲得的圖像序列平均方法示意圖。原始圖像每10幀分為一組，針對一組內的10張圖像的每個像素點的灰度值取均值，得到500frame/s的圖像序列。平均后的圖像序列與10000frame/s幀率下采集的圖像序列一起進行圖像畸變校正。

圖3（a）給出了校正前的圖像。由于相機或實驗裝置擺放位置有偏差，圖像存在一定程度的畸變。選取比例尺上、下沿與滾筒圓周的交點

A、B、C、D4個點作為關鍵點進行透視變換。已知滾筒直徑為 140mm ，根據 ?_AB 、 CD 的弦長即可計算出4點的真實坐標。透視變換的結果如圖3（b）所示。根據圖中比例尺可知，一個像素對應的實際長度約為 0.667mm 。

圖像為 I ，特征模板為 B ，則互相關系數可以表示為

式中： M 為特征模板的邊長； k_ν 1為模板在原始圖像上的位置?；ハ嚓P計算結果如圖4（c）所示。

互相關計算結果依然保持了圓形特征，故采用霍夫變換對互相關計算結果的灰度圖像進行識別?；舴蜃儞Q是將 x-y 平面上的點映射到參數空間的變換，參數空間是由待識別目標的形狀決定的。對于圓形目標而言，其在 x-y 平面上的方程為

r²=（x-a）²+（y-b）²

滾筒側面顆粒的堆積結構一定程度上可以反映整個系統的結構特征。測量其堆積結構首先需要對每個顆粒的位置進行識別。在滾筒側面圖像中選取被動層邊長為 17mm 的正方形區域為例，顆粒識別過程如圖4所示。

顆粒表面在LED燈光照射下的反光具有豐富的灰度信息，表層顆粒與第二層顆粒的反光具有顯著差異。取10個典型表層顆粒的反光圖像的平均值作為顆粒反光的特征模板，如圖4（b）所示。用該模板遍歷原始圖像上的每一個位置，計算原始圖像和特征模板的互相關系數。記原始

可以看出圓有3個參數 r ， Ψ_a 和 b ，其中a 和 b 分別為 x 和 y 方向上的圓心位置， r 為圓的半徑。圓的參數表示為

由式（3）可知，圓的參數空間為 R³ ，即 x. y 平面同一個圓形上所有的點可以經過霍夫變換映射到參數空間中的一個點 （a，b，r） 。霍夫變換檢測圓的基本算法步驟為

1）對原始圖像進行Canny算子邊緣檢測，得到邊緣檢測后的二值圖像；

2）創建一個大小與參數空間一致的三維累加器矩陣 A[a，b，r]

3）對于圖像中的每個邊緣點 （x_i，y_i） ，畫出所有通過該點并且 r 在指定范圍內的圓，并映射到參數空間即可得到一個曲面，曲面上所有點在累加器矩陣中對應位置的值加一；

4）在累加器矩陣 中設置一個閾值搜索局部最大值，獲得圖像中圓的中心位置及半徑。

霍夫檢測結果如圖4（d）所示，紅色圓形即為檢測到的圓形亮點。

最終根據圓心位置得到的所有顆粒中心的坐標如圖4（e）所示。該方法一方面去除了第二層顆粒的反光對表層顆粒識別的干擾[圖4（d）中的紅色方框]，另一方面與直接閾值分割的方法相比，修正了反光點與顆粒中心位置的偏移，有效識別出表層顆粒的位置。

完成每個顆粒位置的識別后采用Voronoi細胞劃分方法對其堆積結構進行表征。給定一組離散的顆粒，每個顆粒的Voronoi細胞是比任何其他顆粒更接近它的空間區域。對于顆粒 作其連線的中垂線將平面劃分為兩個半空間，其中一個半空間中所有的點距離 s_i 更近。根據Voronoi細胞的定義， s_i 的Voronoi細胞由所有最靠近它的點組成，所有距離 s_i 比 s_j 更近的點組成的半空間 H_ij 可以表示為

式中， d 表示空間中點 x 到顆粒 s_i/s_j 的距離。將 s_i 的Voronoi單元定義為所有這些半空間的交集，表示為

V（s_i）=∩_jH_ij

式中， j 表示系統中所有其他顆粒。而遠離中心顆粒的顆粒不影響其Voronoi細胞，Voronoi細胞的形狀通常僅由附近的顆粒決定。由式（4）和式（5）提供的定義迭代地使用相交半空間切割圖4中顆粒 s_i 周圍的區域即可得到整個系統的Voronoi圖，如圖5所示。

圖6為被動層區域Voronoi多邊形面積與邊數的分布，其中S表示單個Voronoi細胞的面積， Ps 表示面積 s 的概率密度函數， n 表示單個Voronoi細胞多邊形的邊數， N 表示每種邊數Voronoi細胞的個數。Voronoi多邊形面積服從伽馬分布，圖中紅色曲線為式（6）預測的伽馬分布，該式是由Aste等[21]提出的三維隨機堆積顆粒系統中Voronoi體積的分布，其中 V 為每個顆粒Voronoi多面體的體積， V_min 為整個顆粒系統內Voronoi多面體的最小體積， 為平均體積。該式是變量 （V-V_min） 的伽馬分布，其特征為形狀參數 k 和尺度參數 ，此處用面積 s 代替式（6中的體積 V ，實驗測量結果與理論預測具有良好的一致性，驗證了側面單層顆粒的堆積結構信息可以反映三維系統的堆積結構特征。邊數服從正態分布，平均邊數為6，最大邊數為8，最小邊數為4。

2 結果與討論

2.1顆粒重排對堆積結構的影響

實驗結果顯示，在崩塌周期中，滾筒中被動層區域靜態堆積的顆粒并不是完全相對靜止的，局部顆粒會發生一些非仿射位移，導致局部結構的變化。這種局部堆積結構的變化也被稱為顆粒重排。重排顆粒的非仿射位移非常微小，從原始圖像中難以直接觀察到其位置變化，故以1.2節中識別出的顆粒坐標為基礎，認為在 0.1ms 內與其相鄰顆粒相對位移超過粒徑 20% 的顆粒參與了顆粒重排，相鄰的重排顆粒組成了一個單次重排事件。滾筒被動層區域的顆粒系統具有許多不同的堆積結構，這些堆積結構具有不同的穩定性，即有些結構在外力作用下相對更容易發生重排。為了考察局部重排區域的結構特征以及重排前后非仿射位移帶來的結構的變化，本節研究了不同重排區域的結構特征。首先對單次重排規模進行統計。圖7給出了每次重排過程中參與重排顆粒個數的統計結果。

圖7（a）給出了測量區域示意圖。Li等[20]的研究顯示顆粒重排事件的分布在崩塌系統上下游之間存在明顯區別，但其在被動層中的傳播沒有固定的方向，故測量區域的劃分參考地面坐標。選取崩塌開始前自由表面最高點的水平位置作為左邊界，如圖7（a）中 A 點所示；崩塌結束時自由表面最低點的水平位置作為右邊界，如圖7（b）中 B 點所示。等間隔將其劃分為七個區域，得到每個區域水平方向上寬度為 17mm 豎直方向上包含該寬度內所有顆粒。統計115個崩塌周期每個區域被動層中每次顆粒重排涉及的顆粒個數，持續時間約 480s 。由于區域六、七的被動層中顆粒在表面顆粒撞擊作用下堆積相對較為緊密，基本不發生顆粒重排，統計結果顯示115個周期中區域六只發生了2次重排，區域七沒有發生重排。可見，區域六、七的結果具有偶然性，故后續的討論中不考慮這兩個區域。區域二至區域四的被動層面積根據區域一的面積進行歸一化，得到歸一化后每個區域中參與重排顆粒個數分布，如圖7（c）所示。

由于只有 1～2 個顆粒參與的顆粒重排對其鄰域顆粒的堆積結構不會造成影響，可以認為其是在“顆粒牢籠”中相對自由的顆粒，不參與后續重排區域結構變化的討論，因此統計結果從3個顆粒參與的顆粒重排開始展示。結果顯示，重排涉及的范圍大部分為 3～6 個顆粒，最大范圍為10個顆粒，且出現頻率大致隨參與重排顆粒數量增加而減小；比較不同區域的結果可以看出，區域一內重排發生的次數與規模都是最大的，且區域一至區域五重排發生的次數與規模逐漸減小。

在參與重排顆粒個數分布結果的基礎上可以對二維顆粒系統中參與重排顆粒的堆積結構進行進一步考察。對重排區域內顆粒進行二維Voronoi劃分，根據顆粒的運動趨勢以及Voronoi圖結構可以將重排顆粒劃分為三種基本團簇：線形、三角形以及四邊形，認為具有相近方向位移的顆粒處于同一顆粒團簇。圖8為重排前后相鄰兩幀圖像中三種基本顆粒團簇對應的顆粒位置及Voronoi圖變化。藍色點為重排前顆粒位置，紅色點為重排后顆粒位置，團簇內顆粒具有相近的運動趨勢，可以認為在重排過程中團簇內顆粒作為一個整體一起運動。

由圖7（b）中給出的參與重排的顆粒個數分布可知，除了只包含單個基本團簇的小范圍重排（ 3～4 個顆粒）之外，還存在許多較大范圍的顆粒重排（ 5～10 個顆粒）。對于這些較大范圍的重排，可以根據不同的運動趨勢將其劃分為幾個基本團簇。圖9給出了一次由8個顆粒參與的顆粒重排過程。

如圖9所示，藍色點為重排前顆粒位置，紅色點為重排后顆粒位置。將該顆粒重排區域劃分為I、Ⅱ和Ⅲ3個基本團簇，其中I、Ⅱ為三角形，Ⅲ為四邊形。區域I作為一個整體向右下方運動，區域Ⅲ作為一個整體向左下方運動，區域Ⅱ作為區域I、Ⅲ之間的橋梁存在，將區域I指向右下方的壓力傳遞到區域Ⅲ中。類似于圖9中給出的劃分方法，對所有顆粒重排進行基本團簇的劃分，統計在115個崩塌周期被動層顆粒重排中3種基本團簇出現的次數，如表1所示。

根據表1給出的結果可知，3種基本團簇中線性團簇出現次數最多，三角形團簇出現次數最少。即線性團簇結構相對最不穩定，更容易出現在重排區域中；而三角形團簇結構相對較為穩定，較少參與顆粒重排。

對區域一中線形團簇指向的方向進行統計，由于該系統中線形團簇在重力作用下均向下運動，取線形團簇中最上端顆粒為起始點，最下端顆粒為終點，團簇的方向向量記為 ，如圖10中插圖所示。計算方向向量 與 z 軸負方向的夾角 θ ，夾角位于 z 軸右側記為正值，位于左側記為負值。 θ 的分布情況如圖10所示。由圖10可知，大部分 θ 分布在 0^°～60^° 之間，即大部分線形團簇的方向指向右下方向。這說明線形團簇在右下方向上具有更高的自由度，而在接近于指向圓心的方向不易發生結構變化。

2.2 崩塌與堆積結構的相關性分析

顆粒系統內部堆積結構的變化在崩塌的發展過程中起到十分重要的作用。本小節主要分析了崩塌間隔階段被動層堆積結構變化的劇烈程度與下一次崩塌之間的相關性。一個崩塌周期中參與崩塌顆粒的最大數量可以用來衡量該次崩塌的規模。圖11給出了崩塌顆粒數量的估計過程。圖11（a）、11（b）為崩塌過程中相鄰兩幀圖像，間隔時間為 2ms 。在滾筒圓周上選取5個點，根據其坐標使用最小二乘法擬合出圓心 o 的位置。已知滾筒以 0.19r/min 的速度逆時針旋轉，為消除滾筒旋轉帶來的顆粒位移的干擾，將圖11（b）以 o 為旋轉中心順時針旋轉 3.98× 10^-5 rad，旋轉前后的圖像中心點保持不變，截取與圖11（a）相同大小的部分，兩者作差，差分結果的偽彩色圖如圖11（c）所示。

統計了100張崩塌間隔期間相對靜止的顆粒系統的差分圖像灰度值分布，結果顯示灰度值最大為0.17。認為灰度值大于0.2的區域的顆粒堆積結構發生了變化。對差分圖像進行二值化處理，如圖12（a）所示；再以直徑為20像素的圓形結構元素對其進行閉運算，以填補二值化圖像中崩塌區域的空白，得到崩塌區域如圖12（b）中白色區域所示。該區域在崩塌周期中的最大深度約為 9.87mm ，與主動層最大深度基本一致，說明該閾值選擇比較合理。估算參與崩塌的顆粒數量為

式中： Ψ_c 為崩塌區域像素數； L 為滾筒長度；ν 為單個顆粒的體積； φ 為顆粒體系的體積分數，根據之前的數值模擬結果可知約為0.58。

被動層中顆粒堆積結構變化的劇烈程度可以用發生重排的顆粒數量來衡量。重排顆粒的統計與崩塌顆粒采用相同的閾值，即認為差分圖像中灰度值大于0.2的被動層顆粒發生了重排。統計22個連續的崩塌周期的最大崩塌顆粒數量以及每個崩塌階段之前一個崩塌間隔中堆積結構發生變化的顆粒的個數，如圖13（a）所示。其中，橫坐標 N_arr 為崩塌間隔期間發生重排的顆粒個數，縱坐標 N_ava 為每個崩塌間隔之后緊跟的一次崩塌中崩塌顆粒的數量?？梢钥闯鰞烧咧g具有較強的線性相關關系，計算得到相關系數 R= 0.8968，擬合出線性相關斜率 k=20.58 。不難看出，崩塌規模與該次崩塌前顆粒系統內重排顆粒個數呈正線性相關，即崩塌前靜止抬升階段重排顆粒數量越多，參與崩塌的顆粒數量越多，崩塌規模越大。崩塌顆粒個數與重排顆粒個數之間大致服從 N_ava=20.58×N_arr+1428 ，線性擬合結果如圖13（a）中直線所示。

認為差分圖像中剛出現灰度值大于0.2的區域時為崩塌開始時刻，沒有灰度值大于0.2的區域時為崩塌終止時刻。圖13（b）給出了崩塌終止后0.1s內發生重排的顆粒個數與該崩塌終止時刻到下個崩塌開始時刻之間的時間間隔之間的關系。縱坐標 t 為崩塌終止時刻到下一次崩塌開始時刻之間的時間間隔。可以看出，兩者具有較強的負相關性，計算得到相關系數 R=-0.8875 ，擬合出線性相關斜率 k=-1.982 。顆粒系統在靜止抬升階段發生重排的次數的增加會使得下一次崩塌的發生提前，兩次崩塌之間的時間間隔 t 與重排顆粒個數 N_arr 之間的線性負相關大致服從t=-1.982×N_arr+275.9 ，線性擬合結果如圖13（b）中直線所示。這說明顆粒重排帶來的微小擾動會使還未達到上休止角的顆粒系統提前發生崩塌，重排顆粒數量越多，擾動越大，崩塌的發生越提前。

由圖13的結果可知，崩塌間隔期間被動層區域顆粒堆積結構變化越劇烈，接下來的一次崩塌規模會越大，并且下一次崩塌的發生會提前。據此，可以通過崩塌前顆粒系統內部堆積結構變化的信息對崩塌行為進行一定程度的預測，為雪崩、泥石流等自然災害的防制提供了更充足的參考依據。

3結論

研究了低速旋轉滾筒中周期性顆粒崩塌系統中被動層的堆積結構，使用顆粒反光的特征信息有效去除了內部顆粒反光對表層顆粒識別的干擾，準確識別出滾筒側面單層顆粒的位置；研究了滾筒側面被動層區域顆粒的堆積結構的分布特征及其在崩塌周期中的變化；使用圖像差分方法估算了崩塌規模，分析了堆積結構變化與崩塌之間的相關性。

1）被動層中局部的堆積結構變化可以劃分為線形、三角形、四邊形3種基本團簇，其中線形團簇最不穩定，更容易參與到局部重排中，而三角形團簇相對最穩定，更不容易參與到重排中。其中線形團簇在垂直于半徑的方向上具有更高的自由度，而在接近于指向圓心的方向不易發生結構變化。

2）實驗結果揭示了滾筒中顆粒被動層靠上游的位置結構相對更脆弱，更容易發生重排。這可能是由于顆粒在主動層下游區域進入被動層時具有較強的各向異性，可以支撐指向圓心方向的壓縮而難以支撐垂直于半徑方向的剪切作用。滾筒內顆粒在崩塌間隔期間重排次數越多意味著主力鏈受到破壞的程度越大，導致被動層堆積結構難以再為主動層顆粒提供足夠的支撐，使得下一個崩塌的規模擴大，崩塌開始時刻提前。這一結論為人們更好地預測與防治雪崩、泥石流、山體滑坡等自然災害提供了進一步的理論支撐。

參考文獻：

[1]CAMPBELL C S. Rapid granular flows[J]. Annual ReviewofFluidMechanics，1990，22：57-90.

[2]楊猛，胡志超，張延化，等.農業顆粒物料氣力清選裝 置研究現狀與展望[J].中國農機化學報，2020，41（3）： 121 -127.

[3]LIK，WANG YF，LINQW，et al.Experiments on granular flow behavior and deposit characteristics： implicationsforrockavalanchekinematics[J]. Landslides，2021，18（5）： 1779-1799.

[4]XUD，LIJ，LIUJ，et al.Advances in continuous flow aerobic granular sludge： a review[J]. Process Safety and Environmental Protection，2022，163： 27-35.

[5]CHAKSMITHANONTP，MILMANF，LEUNG C， etal.Scale-upof granularmaterial flowinanagitated filter dryer[J]. Powder Technology，2022， 407： 117684.

[6］孫其誠，金峰，王光謙.密集顆粒物質的多尺度結構 [J].力學與實踐，2010，32（1）：10－15.

[7]HUTTER K. Avalanche dynamics[M]//SINGHVP. Hydrology of disasters. Dordrecht： Springer，1996： 317-394.

[8]VENIERCM，MARQUEZ DAMIAN S，BERTONE S E，et al. Discrete and continuum approaches for modeling solids motion inside a rotating drumat different regimes[J].Applied Sciences，2021， 11（21）： 10090.

[9]LI K，WANG Y F，LIN Q W， et al. Experiments on granular flow behavior and deposit characteristics： Implicationsforrockavalanchekinematics[J]. Landslides，2021，18：1779-1799

[10]LAZAREA，LUJY，RYCROFTCH.Voronoicell analysis：the shapes of particle systems[J].American Journal of Physics，2022， 90（6）： 469-480.

[11].ZHANGCH，ZHAO SW，ZHAOJD， et al.Threedimensional Voronoi analysis of realistic grain packing： anXCT assisted set Voronoi tessellation framework[J]. Powder Technology，2021，379：251-264.

[12]PHUA A，SMITH J，DAVIES C H J，etal. Understanding the structure and dynamics of local powder packing density variations in metal additive manufacturing using set Voronoi analysis[J]. Powder Technology，2023，418：118272.

[13]華云松，戴紅陽，李宏澤，等.二維顆粒系統中指數分 布的幾何結構[J].光學儀器，2022，44（4）：67－72.

[14].SEGUIN A. Experimental study of some properties of thestrong and weak force networks in a jammed granularmedium[J].GranularMatter，2020，22（2）：48.

[15]BRODU N，DIJKSMANJ A，BEHRINGER R P. Spanning the scales of granular materials through microscopic force imaging[J]. Nature Communications， 2015，6：6361.

[16]BORN P，HOLLDACK K.Analysis ofgranular packing structure by scattering of THz radiation[J]. Review of Scientific Instruments，2017， 88（5）： 051802.

[17]ZENGQ，LIR，LIYM， etal.Recognition of a quasistatic region in a granular bed impacted witha sphere[J].PowderTechnology，2022，407：117612.

[18]王昆，凡鳳仙.滾筒裝備內顆?；旌吓c粉磨研究進展 [J].中國粉體技術，2020，26（4）：38-45.

[19]楊暉，張國華，王宇杰，等.密集顆粒體系的顆粒運動 及結構測量技術[J].力學進展，2018，48（1）：541- 590.

[20]LI R，YANG H， ZHENG G，et al.Double specklevisibility spectroscopy for the dynamics of a passive layerinarotatingdrum[J].Powder Technology，2016， 295：167-174.

[21]ASTE T，DI MATTEO T. Emergence of gamma distributionsin granularmaterials and packing models[J].PhysicalReviewE，2008，77（2）：021309.

（編輯.張磊）