基于動態(tài)神經(jīng)網(wǎng)絡(luò)的水泥熟料堆積體建模方法

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

(1.燕山大學(xué)機械工程學(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]。以上研究都著眼于對單個顆粒建立精準模型。但是有關(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ì)心三維坐標信息，使用神經(jīng)網(wǎng)絡(luò)建立一個預(yù)測模型，以顆粒粒度分布及質(zhì)心三維坐標作為模型的輸入，以顆粒質(zhì)心三維坐標作為輸出，訓(xùn)練神經(jīng)網(wǎng)絡(luò)并利用自相關(guān)函數(shù)和互相關(guān)函數(shù)驗證網(wǎng)絡(luò)。預(yù)測模型只需輸入顆粒集合的粒度分布，網(wǎng)絡(luò)就可以預(yù)測出各個顆粒的質(zhì)心三維坐標，進而實現(xiàn)對顆粒堆積體的建模。這種建模方法可以推廣到其他顆粒物質(zhì)。

1 模型重建方法

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

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

圖1 研究方法Fig.1 Research method

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

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

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

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

在二值化處理結(jié)束后，使用ImageJ軟件進行圖像的定量分析，得到切片中某一顆粒截面面積。隨后按面積等效原則求出相應(yīng)的等效圓直徑。將切片內(nèi)顆粒編號，然后進行切片內(nèi)所屬各顆粒的粒度判斷。在任一所切截面區(qū)域，顆粒的切片數(shù)量及其最大等效圓直徑與標準粒度在Z方向上的關(guān)系如表1所示。按面積等效原則求出的最大等效圓直徑應(yīng)該接近標準粒度，當(dāng)切片數(shù)相同時按最大等效直徑確定標準粒度。通過某一顆粒所占切片數(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ù)量及其最大等效圓直徑與標準粒度的關(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)顆粒進行接觸化后處理。采用基于離散單元法(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坐標的逐漸增加就是一個有力支撐。因此可按照神經(jīng)網(wǎng)絡(luò)中的時間序列模式解決熟料堆積問題。為了獲得足夠長的序列，按照組別在Z方向上隨機排列，組內(nèi)順序不變原則組合數(shù)據(jù)集。構(gòu)建完成后顆粒集合體的輸入向量是時間的函數(shù)，堆積體的半徑與坐標數(shù)據(jù)如圖7所示。

圖7 堆積體的半徑與坐標數(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ì)心三維坐標作為輸出，顆粒質(zhì)心三維坐標的下一個值可由顆粒坐標的先前值和粒度的先前值回歸得到。將采集到的數(shù)據(jù)進行分割，然后利用貝葉斯正則化方法訓(xùn)練網(wǎng)絡(luò)，以均方差為性能指標，訓(xùn)練完成后使用自相關(guān)函數(shù)和互相關(guān)函數(shù)進行驗證。

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

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

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

(1)

式中：nm為歸一化后的輸入向量；n為輸入向量元素；nmin為輸入向量每一維度上最小值的向量；nmax為輸入向量每一維度上最大值的向量。首先將數(shù)據(jù)集按式(1)進行歸一化。歸一化后的數(shù)據(jù)通常在[-1,1]標準的范圍內(nèi)，便于網(wǎng)絡(luò)的訓(xùn)練，歸一化后的半徑與坐標數(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 歸一化后的半徑與坐標數(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ì)心三維坐標q作為輸出。構(gòu)造一個元胞數(shù)組p作為輸入值，該元胞數(shù)組包含先前的輸入和輸出元素。任一顆粒的當(dāng)前質(zhì)心三維坐標既受其粒度值又受相接觸顆粒的粒度值和質(zhì)心三維坐標影響，在網(wǎng)絡(luò)中表現(xiàn)為抽頭延長線(Tapped Delay Line,TDL)[13]的長度，最終抽頭延長線的長度需要根據(jù)訓(xùn)練效果而確定。使用串并聯(lián)結(jié)構(gòu)配置網(wǎng)絡(luò)。

輸入值p:

(3)

目標為下一個輸出值t：

(4)

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

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

訓(xùn)練之前需要將數(shù)據(jù)集進行分割，切片試驗中共采集到1 363個有效數(shù)據(jù)，按式(5)進行數(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]對其進行過細致的研究，并提出了描述性改進。貝葉斯分析最早由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ù)遠大于有效參數(shù)數(shù)量時，說明網(wǎng)絡(luò)規(guī)模太大，可以適當(dāng)減小網(wǎng)絡(luò)規(guī)模。通過對兩者進行比較可以幫助尋找到最簡網(wǎng)絡(luò)，最簡網(wǎng)絡(luò)可以避免數(shù)據(jù)過擬合并減少計算時間。

對網(wǎng)絡(luò)進行多次訓(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ò)在標準化后訓(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ò)在標準化后訓(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ò)沒有進行外推，總體誤差較小，在-0.2～0.2之間，說明網(wǎng)絡(luò)能力足夠。網(wǎng)絡(luò)模型在一定程度上可以滿足堆積體建模的要求。

4 結(jié) 論

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

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

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