聞 巖，黃泊霖，袁 林，劉浩然，鞠艷旭

(1.燕山大學(xué)機(jī)械工程學(xué)院，秦皇島 066004;2.燕山大學(xué)信息科學(xué)與工程學(xué)院，秦皇島 066004)

0 引 言

顆粒物質(zhì)廣泛存在于制藥、水泥、化工等工業(yè)生產(chǎn)中。但是涉及顆粒運動行為時，許多基本問題還沒有弄清楚。顆粒聚集體流動時表現(xiàn)出類似流體的運動行為，但是使用流體力學(xué)理論對其描述并不合適，究其原因在于顆粒體系復(fù)雜的力鏈網(wǎng)絡(luò)。力鏈網(wǎng)絡(luò)不僅支撐著顆粒體系的部分重量還能承擔(dān)一定的外載荷，其在顆粒接觸時開始出現(xiàn)，并隨著顆粒流動不斷演變[1]。研究顆粒動力學(xué)問題首先需要建立靜力學(xué)模型，對基于數(shù)值模擬的顆粒流動研究而言，建立顆粒體系初始堆積狀態(tài)的離散元模型十分必要。為了建立能夠反映顆粒接觸關(guān)系的聚集體模型，國內(nèi)外學(xué)者在建模方法上展開了一系列的研究。在單顆粒方面，Nolan等[2]提出了利用數(shù)個球體的包絡(luò)面逼近特定顆粒輪廓的單顆粒建模方法。對于特定顆粒的輪廓特征，多數(shù)研究者先通過掃描技術(shù)獲取顆粒的外形特征，再利用較小粒度的顆粒球?qū)ζ涮畛鋄3-5]。以上研究都著眼于對單個顆粒建立精準(zhǔn)模型。但是有關(guān)復(fù)雜顆粒體系建模的研究還比較少。在體系內(nèi)部通常包含大量顆粒物質(zhì)，粒度分布存在一定范圍，并且各個顆粒間又存在著相互作用。為了對特定顆粒群展開研究，所建立的離散元模型應(yīng)考慮體系內(nèi)顆粒的真實接觸情況，應(yīng)盡可能與真實的顆粒模型等價[6-8]。

水泥熟料是比較典型的顆粒物質(zhì)，熟料堆積體在生產(chǎn)中因受推動棒作用而往復(fù)運動，為了利用離散元軟件研究水泥熟料的運動行為，就需要先建立堆積狀態(tài)下的離散元模型，但是關(guān)于水泥熟料堆積體建模的報道較少。其中，張文明等[9]先利用雙目視覺技術(shù)重建熟料層，再通過網(wǎng)格劃分技術(shù)還原熟料幾何外形，實現(xiàn)了水泥熟料顆粒體系的建模。但是并非以力鏈形成為線索。

本文提出一種基于神經(jīng)網(wǎng)絡(luò)的顆粒體系建模方法。首先通過試驗獲取顆粒的粒度與質(zhì)心三維坐標(biāo)信息，使用神經(jīng)網(wǎng)絡(luò)建立一個預(yù)測模型，以顆粒粒度分布及質(zhì)心三維坐標(biāo)作為模型的輸入，以顆粒質(zhì)心三維坐標(biāo)作為輸出，訓(xùn)練神經(jīng)網(wǎng)絡(luò)并利用自相關(guān)函數(shù)和互相關(guān)函數(shù)驗證網(wǎng)絡(luò)。預(yù)測模型只需輸入顆粒集合的粒度分布，網(wǎng)絡(luò)就可以預(yù)測出各個顆粒的質(zhì)心三維坐標(biāo)，進(jìn)而實現(xiàn)對顆粒堆積體的建模。這種建模方法可以推廣到其他顆粒物質(zhì)。

1 模型重建方法

高溫環(huán)境中的水泥熟料在工業(yè)生產(chǎn)時易受到風(fēng)室中冷空氣的影響而突然冷卻，變?yōu)榘肴廴跔顟B(tài)，隨后在凝結(jié)破裂后不斷堆積，最終聚集成一定數(shù)量的顆粒集合，即顆粒堆積體[10]。這個過程中不規(guī)則的熟料因邊角鈍化，近似成球體。本文在進(jìn)行處理時亦將顆粒等效成球體。熟料堆積體內(nèi)部含有大量形狀各異的孔隙結(jié)構(gòu)，常稱其為多孔介質(zhì)。此類結(jié)構(gòu)可以采用切片法、立體成像法和圖像法進(jìn)行三維重建[11]。

本文擬采用切片法對熟料堆積體進(jìn)行三維重建，以期獲取顆粒集合中各個顆粒的粒度值和質(zhì)心三維坐標(biāo)信息。研究方法如圖1所示，以工業(yè)生產(chǎn)的水泥熟料為試驗樣本構(gòu)建堆積體，制作切片并采集圖像。通過圖像及數(shù)據(jù)處理過程，獲取包含顆粒三維坐標(biāo)及粒度值的空間信息。將空間信息進(jìn)行歸一化處理，利用動態(tài)神經(jīng)網(wǎng)絡(luò)建立預(yù)測模型并對網(wǎng)絡(luò)進(jìn)行訓(xùn)練，訓(xùn)練完成后對網(wǎng)絡(luò)進(jìn)行驗證。

圖1 研究方法Fig.1 Research method

2 試驗方法及數(shù)據(jù)處理

為了獲取相關(guān)數(shù)據(jù)，試驗時以工業(yè)生產(chǎn)的水泥熟料為試驗材料，模擬熟料在篦冷機(jī)內(nèi)的堆積過程。

首先選取樣本顆粒，考慮到在堆積過程中會出現(xiàn)邊壁效應(yīng)，即顆粒與硬邊界直接接觸而導(dǎo)致的對堆積成型效果的影響，試驗中在盒體內(nèi)壁附著一層厚達(dá)20 mm的保溫棉，承載熟料-石膏混合體容器軟邊壁的邊長為150 mm,料盒硬邊界總體尺寸達(dá)到190 mm。使用樣本顆粒密實填充總體尺寸為190 mm的無蓋立方體料盒，使用石膏浸漬劑固化顆粒堆積體，固化后的堆積體如圖2所示。因為顆粒粒度分布絕大部分在20～60 mm之間，在試驗操作中，首先沿固化熟料體Z軸方向，使用切割機(jī)按20 mm等距切割堆積體，然后用相機(jī)采集切片圖像并編號，切片圖像如圖3所示。其中乳白色部分為浸漬劑，黑色部分為水泥熟料。

在數(shù)據(jù)處理之前，先對所采集的圖像進(jìn)行預(yù)處理。利用ImageJ軟件對圖像尺寸進(jìn)行標(biāo)定，使計算機(jī)中圖像的像素所代表的尺寸與試驗標(biāo)尺測定的切片真實長度相等，因采集到圖像為RGB模式需進(jìn)行二值化處理，處理后如圖4所示。

在二值化處理結(jié)束后，使用ImageJ軟件進(jìn)行圖像的定量分析，得到切片中某一顆粒截面面積。隨后按面積等效原則求出相應(yīng)的等效圓直徑。將切片內(nèi)顆粒編號，然后進(jìn)行切片內(nèi)所屬各顆粒的粒度判斷。在任一所切截面區(qū)域，顆粒的切片數(shù)量及其最大等效圓直徑與標(biāo)準(zhǔn)粒度在Z方向上的關(guān)系如表1所示。按面積等效原則求出的最大等效圓直徑應(yīng)該接近標(biāo)準(zhǔn)粒度，當(dāng)切片數(shù)相同時按最大等效直徑確定標(biāo)準(zhǔn)粒度。通過某一顆粒所占切片數(shù)及最大等效直徑確定顆粒粒度范圍，并給出其所在區(qū)間的粒度初值。對于只占一個切面的顆?；蜉^小空白處，可按10 mm顆粒給定初值。

圖2 堆積體等效圖Fig.2 Plane equivalent image ofaccumulation body

圖3 切片原始圖Fig.3 Original slice image

圖4 切片二值化圖Fig.4 Slice binarization image

表1 顆粒的切片數(shù)量及其最大等效圓直徑與標(biāo)準(zhǔn)粒度的關(guān)系Table 1 Relationship between the slice number and maximum equivalent diameter of particles and the normal diameter

圖5 切片質(zhì)心圖Fig.5 Slice centroid image

圖6 顆粒切片中的幾何關(guān)系Fig.6 Geometric relationship of particle slices

在按上述方法處理之后，一些顆粒會出現(xiàn)過接觸或懸浮狀，需要對模型內(nèi)顆粒進(jìn)行接觸化后處理。采用基于離散單元法(Discrete Element Method，DEM)開發(fā)的離散元軟件EDEM支持的顆粒工廠編譯方式，將模型內(nèi)置到顆粒工廠中以EDEM依賴的顆粒接觸理論重構(gòu)堆積體。

試驗共制作了10組熟料堆積體，每個固化熟料獲得6～8個有代表性的切片，共獲得65張有效圖片，共統(tǒng)計到1 363個數(shù)據(jù)點。

3 預(yù)測模型建立

熟料堆積過程與神經(jīng)網(wǎng)絡(luò)中的時間序列問題較為符合，按這類問題的模式建立熟料堆積預(yù)測模型。堆積體一般由各個顆粒在一段時間內(nèi)逐漸聚合而成，在顆粒堆積完成后，每一顆粒的終態(tài)都對應(yīng)堆積過程中某個時間點，各個顆粒在時間上表現(xiàn)出明顯的先后順序，集合中顆粒Z坐標(biāo)的逐漸增加就是一個有力支撐。因此可按照神經(jīng)網(wǎng)絡(luò)中的時間序列模式解決熟料堆積問題。為了獲得足夠長的序列，按照組別在Z方向上隨機(jī)排列，組內(nèi)順序不變原則組合數(shù)據(jù)集。構(gòu)建完成后顆粒集合體的輸入向量是時間的函數(shù)，堆積體的半徑與坐標(biāo)數(shù)據(jù)如圖7所示。

圖7 堆積體的半徑與坐標(biāo)數(shù)據(jù)Fig.7 Radius and coordinate data of the accumulation body

按以上分析，可選擇非線性自回歸(Nonlinear Autoregressive with Exogenous Input,NARX)神經(jīng)網(wǎng)絡(luò)來構(gòu)建預(yù)測模型。將顆粒的粒度值作為輸入，質(zhì)心三維坐標(biāo)作為輸出，顆粒質(zhì)心三維坐標(biāo)的下一個值可由顆粒坐標(biāo)的先前值和粒度的先前值回歸得到。將采集到的數(shù)據(jù)進(jìn)行分割，然后利用貝葉斯正則化方法訓(xùn)練網(wǎng)絡(luò)，以均方差為性能指標(biāo)，訓(xùn)練完成后使用自相關(guān)函數(shù)和互相關(guān)函數(shù)進(jìn)行驗證。

3.1 網(wǎng)絡(luò)結(jié)構(gòu)

數(shù)據(jù)集收集完成后需進(jìn)行歸一化處理，網(wǎng)絡(luò)結(jié)構(gòu)選定后其中的輸入和輸出需要定義。

nm=2(n-nmin)/(nmax-nmin)-1

(1)

式中：nm為歸一化后的輸入向量；n為輸入向量元素；nmin為輸入向量每一維度上最小值的向量；nmax為輸入向量每一維度上最大值的向量。首先將數(shù)據(jù)集按式(1)進(jìn)行歸一化。歸一化后的數(shù)據(jù)通常在[-1,1]標(biāo)準(zhǔn)的范圍內(nèi)，便于網(wǎng)絡(luò)的訓(xùn)練，歸一化后的半徑與坐標(biāo)數(shù)據(jù)如圖8所示。

NARX神經(jīng)網(wǎng)絡(luò)是一個帶有外部輸入的非線性自回歸網(wǎng)絡(luò)，廣泛用于預(yù)測問題和動態(tài)建模[12]。選擇其作為堆積體建模的網(wǎng)絡(luò)結(jié)構(gòu)。NARX網(wǎng)絡(luò)結(jié)構(gòu)如圖9所示。

NARX模型的定義方程如下：

y(t)=f(y(t-1),y(t-2),…,y(t-ny),u(t-1),u(t-2),…,u(t-nu))

(2)

式中：y為因變量輸入信號；t為時間；u為自變量輸入信號；ny為跟因變量y相關(guān)的自然數(shù)；nu為跟自變量u相關(guān)的自然數(shù)。公式(2)可以這樣理解，t時刻的輸入信號與t-1到t-n時刻的自變量和因變量有關(guān)。

圖8 歸一化后的半徑與坐標(biāo)數(shù)據(jù)Fig.8 Normalized radius and coordinate data

圖9 NARX網(wǎng)絡(luò)結(jié)構(gòu)Fig.9 NARX network structure

網(wǎng)絡(luò)訓(xùn)練時將顆粒半徑R作為輸入，將顆粒的質(zhì)心三維坐標(biāo)q作為輸出。構(gòu)造一個元胞數(shù)組p作為輸入值，該元胞數(shù)組包含先前的輸入和輸出元素。任一顆粒的當(dāng)前質(zhì)心三維坐標(biāo)既受其粒度值又受相接觸顆粒的粒度值和質(zhì)心三維坐標(biāo)影響，在網(wǎng)絡(luò)中表現(xiàn)為抽頭延長線(Tapped Delay Line,TDL)[13]的長度，最終抽頭延長線的長度需要根據(jù)訓(xùn)練效果而確定。使用串并聯(lián)結(jié)構(gòu)配置網(wǎng)絡(luò)。

輸入值p:

(3)

目標(biāo)為下一個輸出值t：

(4)

3.2 網(wǎng)絡(luò)訓(xùn)練

訓(xùn)練時采用貝葉斯正則化方法，選擇均方誤差(MSE)作為性能指標(biāo)。

訓(xùn)練之前需要將數(shù)據(jù)集進(jìn)行分割，切片試驗中共采集到1 363個有效數(shù)據(jù)，按式(5)進(jìn)行數(shù)據(jù)集劃分[14],其中ropt為最優(yōu)劃分率，m為網(wǎng)絡(luò)結(jié)構(gòu)參數(shù)[15]，即網(wǎng)絡(luò)結(jié)構(gòu)中權(quán)值和偏置值的總和。因為訓(xùn)練時使用了貝葉斯正則化訓(xùn)練算法，所以不需要單獨劃分驗證集，由于將熟料堆積問題處理為時間序列問題，選擇最后連續(xù)的集合作為測試集。使用Widrow/Nguyen方法初始化權(quán)值[16]。

圖10 均方誤差隨迭代次數(shù)的關(guān)系Fig.10 Relationship between mean square error and iteration number

使用貝葉斯正則化訓(xùn)練方法，該方法包含貝葉斯分析和正則化法兩個概念，Mackay[17]對其進(jìn)行過細(xì)致的研究，并提出了描述性改進(jìn)。貝葉斯分析最早由Bayes提出[18]，廣泛應(yīng)用于統(tǒng)計學(xué)，正則化法常用于提升網(wǎng)絡(luò)的泛化能力。貝葉斯正則化方法在訓(xùn)練時會自動選擇正則化參數(shù)[19]，并將神經(jīng)網(wǎng)絡(luò)置于貝葉斯統(tǒng)計框架中。正確的正則化參數(shù)可使網(wǎng)絡(luò)具備良好的泛化能力，正則化概念最早由Tikhonov提出[20]。另外，貝葉斯正則化方法在訓(xùn)練時會計算一個有效參數(shù)數(shù)量[21]。當(dāng)網(wǎng)絡(luò)結(jié)構(gòu)參數(shù)總數(shù)接近有效參數(shù)數(shù)量時，說明網(wǎng)絡(luò)規(guī)模太小，應(yīng)該增大神經(jīng)元數(shù)或者抽頭延長線數(shù)；當(dāng)網(wǎng)絡(luò)結(jié)構(gòu)參數(shù)總數(shù)遠(yuǎn)大于有效參數(shù)數(shù)量時，說明網(wǎng)絡(luò)規(guī)模太大，可以適當(dāng)減小網(wǎng)絡(luò)規(guī)模。通過對兩者進(jìn)行比較可以幫助尋找到最簡網(wǎng)絡(luò)，最簡網(wǎng)絡(luò)可以避免數(shù)據(jù)過擬合并減少計算時間。

對網(wǎng)絡(luò)進(jìn)行多次訓(xùn)練發(fā)現(xiàn)性能函數(shù)MSE曲線相似，說明網(wǎng)絡(luò)沒有收斂到局部極小值。均方誤差隨迭代次數(shù)關(guān)系如圖10所示。經(jīng)驗證網(wǎng)絡(luò)中抽頭延長線設(shè)為20，神經(jīng)元設(shè)為10,此時網(wǎng)絡(luò)最小，訓(xùn)練效果較好。這里網(wǎng)絡(luò)中有843個結(jié)構(gòu)參數(shù)，按式(5)計算，訓(xùn)練集取總數(shù)據(jù)的98%，測試集取2%。

(5)

3.3 網(wǎng)絡(luò)驗證

使用自相關(guān)函數(shù)和互相關(guān)函數(shù)驗證網(wǎng)絡(luò)。訓(xùn)練好的網(wǎng)絡(luò)應(yīng)該具備兩個特性，預(yù)測誤差應(yīng)當(dāng)在時間上不相關(guān)，預(yù)測的誤差應(yīng)當(dāng)與輸入序列不相關(guān)[22]。

自相關(guān)函數(shù)Re(τ)：

(6)

式中：e(t)表示t時刻的誤差；e(t+τ)表示t+τ時刻的誤差；τ為時間增量；Q為總時間。

定義Re(τ)的置信區(qū)間：

(7)

網(wǎng)絡(luò)在標(biāo)準(zhǔn)化后訓(xùn)練的預(yù)測誤差的自相關(guān)函數(shù)如圖11所示，虛線為置信區(qū)間范圍。從圖中可以看出自相關(guān)函數(shù)在0時是一個脈沖，其他情況下，函數(shù)值為0，所以預(yù)測誤差在時間上不相關(guān)。

互相關(guān)函數(shù)Rre(τ)：

(8)

式中：r(t)表示t時刻的輸入；e(t+τ)表示t+τ時刻的誤差；τ為時間增量；Q為總時間。

定義Rre(τ)的置信區(qū)間：

(9)

式中：Rr(0)表示輸入r(t)的自相關(guān)函數(shù)。

圖11 自相關(guān)誤差函數(shù)圖Fig.11 Autocorrelation error function diagram

圖12 互相關(guān)誤差函數(shù)圖Fig.12 Cross correlation error function diagram

網(wǎng)絡(luò)在標(biāo)準(zhǔn)化后訓(xùn)練的預(yù)測誤差的互相關(guān)函數(shù)如圖12所示，虛線為置信區(qū)間范圍。從圖中可以看出互相關(guān)函數(shù)在置信區(qū)間內(nèi)，所以預(yù)測誤差與輸入序列不相關(guān)。

通過以上分析表明已經(jīng)建立完成一個精確的預(yù)測模型。最終預(yù)測模型的誤差如圖13所示。

圖13 預(yù)測誤差Fig.13 Predictive error

從圖13可以看出歸一化后的訓(xùn)練集和測試集誤差接近，說明網(wǎng)絡(luò)沒有進(jìn)行外推，總體誤差較小，在-0.2～0.2之間，說明網(wǎng)絡(luò)能力足夠。網(wǎng)絡(luò)模型在一定程度上可以滿足堆積體建模的要求。

4 結(jié) 論

(1)水泥熟料堆積體的模型比較復(fù)雜，神經(jīng)網(wǎng)絡(luò)功能強(qiáng)大且善于解決工業(yè)問題。熟料堆積過程與神經(jīng)網(wǎng)絡(luò)中的時間序列問題類似，按照時間序列模式處理熟料堆積過程。選擇NARX網(wǎng)絡(luò)結(jié)構(gòu)，構(gòu)建一個包含輸入和輸出元素的元胞數(shù)組，使用貝葉斯正則化方法訓(xùn)練網(wǎng)絡(luò)，以均方誤差作為性能指標(biāo)。

(2)訓(xùn)練神經(jīng)網(wǎng)絡(luò)時需要先獲取原始數(shù)據(jù)。為了獲取熟料堆積體的結(jié)構(gòu)信息，以工業(yè)生產(chǎn)的水泥熟料為材料進(jìn)行了堆積試驗，并制作了一系列熟料切片，通過對切片圖像進(jìn)行處理和統(tǒng)計獲得了堆積體中顆粒集合的粒度和質(zhì)心三維坐標(biāo)。

(3)對訓(xùn)練后的神經(jīng)網(wǎng)絡(luò)進(jìn)行驗證分析，預(yù)測誤差的自相關(guān)函數(shù)和互相關(guān)函數(shù)都位于置信區(qū)間內(nèi)，說明誤差與誤差，誤差與輸入之間沒有相關(guān)性，模型可以正確運行。歸一化后預(yù)測模型的總誤差在-0.2～0.2之間，網(wǎng)絡(luò)模型在一定程度上可以滿足堆積體建模的要求。