基于模式識別技術的鐵路道砟模型庫建立及離散元模擬

賀 勇，張宗堂，樊寶杰，柳光磊，張巨峰

(1.長沙市規劃勘測設計研究院，湖南 長沙 410007；2.湖南科技大學 巖土工程穩定控制與健康監測湖南省重點實驗室，湖南 湘潭 411201；3.湖南科技大學 資源環境與安全工程學院，湖南 湘潭 411201；4.湖南省地質礦產勘查開發局四一四地質隊，湖南 益陽 413000)

1 概述

有砟鐵路鋪設了大量的道砟顆粒[1-3]。傳統的顆粒形狀參數測量法通過手動測量或與標準圖表視覺比較來確定參數[4]。這些傳統方法存在很多局限性，例如易受噪聲，振動的影響[5]，手動測量和視覺比較法測量非常主觀且緩慢[6]。因此，需要自動、快速、準確的鐵路道砟顆粒測量方法來代替傳統的測量法。

當顆粒彼此不接觸的情況下，獲取圖像并進行圖像分析相對容易。RASCHKE、BANTA、ALTUHAFI等[7-9]都使用了這種分離方法。但在實際情況下，鐵路道砟顆粒都是彼此接觸的，因此難以對其進行顆粒分離。

在接觸顆粒的圖像中，計算鐵路道砟顆粒形狀更加困難，首先需要在圖像中識別出顆粒，并準確定義其邊界。此外，考慮三維鐵路道砟顆粒，道砟顆粒間不僅彼此接觸，而且會有顆粒重疊，在相機視圖中存在遮擋。因此，只有小部分道砟顆?？梢酝暾仫@示二維投影。只有圖像中完整地顯示出連續輪廓的顆粒才可用于形狀測量。因此，在ZHENG等[10]提出的一種半自動方法中，仍然需要科研人員在道砟顆粒的圖像中挑選出完整輪廓顆粒。

本文基于模式識別，提出了一種自動識別完整輪廓顆粒的算法，省去了人力工作，對道砟顆粒形狀評估結果進行統計分析。建立了顆粒形狀庫，用于在數值模型中隨機重建具有所需形狀值的鐵路道砟顆粒?；谒岢龅牡理念w粒形狀庫和離散元方法，對真實形狀的道砟顆粒進行了雙軸壓縮試驗，研究了道砟顆粒細長度對其力學特性的影響。

2 道砟顆粒的自動獲取

為節省時間和成本，并考慮到鐵路道砟樣本數量多、精度要求低的特點，本文采用了低精度的攝影設備(例如手機和數碼相機)進行拍攝，獲得了多個鐵路道砟顆粒圖像，如圖1所示。本文提取的目標是道砟顆粒的不規則形狀。

圖1 道砟圖像示例

由于算法的簡單性、高精度性和高計算效率，在本文中使用經過改進的Viola-Jones算法[11]來自動識別道砟顆粒。

整個過程主要分為4個主要步驟：訓練集準備，訓練集圖像的縮放和旋轉，Adaboost和Cascade訓練分類器，基于滑動窗口方法識別顆粒。

2.1 訓練集準備

首先，考慮到道砟顆粒的最常見形狀，其長寬比基本為1：1。拍攝1 206張道砟照片，對其中的一些道砟圖像進行裁剪、縮放并作為模型的訓練數據集。每個裁剪后的圖像都是原圖像中的一個小方塊。其中，包含道砟顆粒完整輪廓的圖像被標記為正圖像，而其他圖像被標記為負圖像。所有的裁剪圖像的大小均為32×32像素。將這些圖像作為輸入數據，計算機可以進行模型訓練和學習來獲得模式識別中的Haar-Like特征[12]。并將特征結果存儲在分類器中。試樣圖像如圖2所示。

(a)正圖像

2.2 訓練集圖像的縮放和旋轉

大多數鐵路道砟顆粒的長寬比為1:1至2:1。對于不處于該比例的道砟，很難通過拍攝的方式來大量收集圖像數據，但這些高長寬比的道砟顆粒是真實存在的，所以分類器也需要檢測到這些道砟。因此，需要一種新的獲取大量圖像的方法。為了快速獲得大量可用于訓練的細長道砟顆粒圖像，本文將道砟長寬比為1:1的現有圖像進行數字拉伸[13]。通過調整之前的1:1訓練圖像(32×32像素)的大小，生成1.5:1訓練圖像(48×32像素)，2:1訓練圖像(64×32像素)，3:1訓練圖像(96×32)像素)和4:1訓練圖像(128×32像素)，如圖3所示。

(a)1.5:1圖像

為了在任何方向上都能識別細長道砟顆粒，原始圖像以15°的增量逆時針旋轉，如圖4所示。對于每個旋轉的圖像，滑動窗依次掃描整個圖像。由于某些道砟顆粒傾斜放置，滑動窗法可能無法檢測到這些道砟。但是如果旋轉圖像，這些道砟將轉換為豎直方向上，易于識別。將圖像旋轉15°、30°、45°、60°、75°、90°、105°、120°、135°、150°、165°后，可以獲得13張旋轉的圖像。圖4給出了0°、30°、60°的簡單示例。

(a)0°

合并這些圖像的所有結果。在每個圖像中，假如成功識別出道砟，則添加一個道砟顆粒的邊界框。然后，在合并時所有邊界框會被疊加，但由于旋轉，同一道砟會生成一些不同的邊界框。此時，將圖像旋轉到原點0°位置，這相當于只旋轉了滑動窗。接著選擇面積最小的邊界框作為此道砟的最終邊界框，如圖5所示。

(a)所有旋轉合并結果

2.3 Adaboost和Cascade訓練分類器

為了比較數據集中的正負圖像，將所有標記的圖像用來訓練模型，得到正負圖像之間的特征差異。首先，使用Adaboost和Cascade方法[14]。然后通過滑動窗口方法[15]，即可成功識別道砟顆粒。

Adaboost是一種分類算法，可以將一些弱分類器整合為一個強分類器。而且該強分類器可用于Cascade方法。對于m個正像(yi=1)和n個負像(yi=0)，具有m+n個實例(xi，yi)的訓練數據，其中xi表示圖像，y表示圖像的正負。對于某些給定的值T，該算法迭代t=1，2，…，T次。對于數據實例的權重wi，將每個實例初始化為wi=1。在每次迭代中，將使用相同的分類方法(如決策樹)，但數據集的權重會在每次迭代后會發生變化，增加所有錯誤分類實例的權重，因此，在下一次迭代中會更注意這些錯誤分類的點。在第t次迭代中，使用當前權重wi構建函數ft并在訓練數據集上計算誤差errt，如果errt=0(即該方法在訓練數據集上非常準確)或errt≥0.5，算法終止；否則，對于每個正確分類的實例，本文保持其權重不變，而對于每個錯誤分類的實例，將權重乘一個因子(1-errt)/errt。最后，返回分類器的線性組合，其中第t個分類器的輸出權重為log(1-errt/errt)。將分類器用于測試一個實例，每一個分類器ft將輸出2類的概率分布。將這些概率按相關乘數的比例進行組合。詳細算法如下所示。

圖6 Adaboost算法

Cascade方法的訓練過程分為多個階段。在每個階段，通過Adaboost算法訓練一個強分類器。此外，需要預先給定3個參數：階段個數s，檢測率d和假陽性率f。其中d表示正確檢測的最小百分比，f表示錯誤識別為正的負圖像的最大百分比。每個階段中，需在Adaboost算法中搜索最小的T，以使由Adaboost訓練的強分類器C(x)的檢測率大于d且假陽性率小于f。

在如圖7所示的假定示例中，將100個正圖像和100個負圖像示例作為輸入數據。每個階段的原始參數為s=3，d=0.99，f=0.25。然后，在每個階段中，通過Adaboost訓練一個具有適當T的強分類器，并分別假設T為3、10、35。對于第一階段，模型有200個輸入圖像，且分類器滿足檢測率大于d且假陽性率小于f的要求。由于參數d=0.99，因此至少100×0.99=99個圖像被正確分類為正，而f=0.25，因此至多100×0.99=99個圖像被正確分類為正，而f=0.25，因此至多100×0.25=25被錯誤分類。因此，分類器將輸出124個正圖像，其中99個是真實顆粒，將124個圖像作為階段二的輸入。對于輸出為負的76個圖像，則直接拒絕。然后對第二階段進行類似的計算，由于只有99張正圖像，因此檢測率為98/99=0.99，假陽性率為6/25=0.24，可以滿足要求。然后，在階段三中，可以得到98個圖像的輸出，其中包含97個顆粒和1個非顆粒，檢測率和假陽性率分別為97/98=0.99和1/6=0.17。因此，在s=3個階段之后，可得到總檢測率D=d3=0.97和假陽性率F=f3=0.015。

圖7 Cascade算法的一個假定示例

2.4 基于滑動窗口方法識別顆粒

通過之前得到的分類器，可以獲得包含真實顆粒的圖像。接下來需要識別完整輪廓道砟顆粒。給定一個鐵路道砟顆粒新圖像，使用滑動窗口技術，在每個掃描位置，將滑動窗口內的區域輸入到檢測器中。檢測器將在窗口區域內找到并提取特征?；谶@些特征，分類器最終確定該圖像是正圖像還是負圖像。如果為正圖像，則在掃描位置添加一個邊界框，從而將其標識為包含完整輪廓顆粒的區域。具體說明如下。

由于顆粒的最小像素約為50×50，因此，首先使用50×50的滑動窗口掃描整個圖像。在每個窗口中，將裁剪圖輸入分類器。由于先前的分類器都是在32×32圖像的基礎上訓練的，因此，大于32×32的裁剪圖像應首先縮小為32×32，并利用Cascade算法中的C1(x)分類器進行測試。如果輸出為正，則進入Cascade方法的階段2。否則，將拒絕負圖像，且滑動窗口將以1個像素移動到下一個位置。如果裁剪圖像成功通過所有s個分類器，則它將有很大概率是包含完整輪廓道砟顆粒的正圖像。然后，將在此位置的滑動窗口標記為該顆粒的邊界框。接下來，繼續移動滑動窗口直到掃描完整個圖像?；瑒哟翱谠谒椒较蚝痛怪狈较蛏系囊苿釉隽烤鶠?個像素，并且移動路徑如圖8(a)所示。

當50×50滑動窗口成功掃描整個圖像時，在橫縱兩個方向上將滑動窗口的大小加5個像素，如圖8(b)所示。并重復之前的滑動過程，直到滑動窗口擴展到最大尺寸300×300像素為止。由于大多數裁剪后的圖像都會被分類為負圖像并被拒絕，因此只有極少數圖像可以成功通過所有s個分類器，故該算法計算效率較高。

(a)最小滑動窗口

當滑動窗口足夠大時，滑動窗口的幾次移動都可能包含同一個顆粒，從而產生多個邊界框。在這種情況下，取平均值作為一個新的單個窗口，如圖9所示。

(a)顆粒的多個邊界框

3 形狀分析和顆粒庫的建立

在分類器識別出完整輪廓顆粒后，本文通過一系列計算幾何算法來計算道砟顆粒形狀參數。主要通過3個幾何形狀評價指標來描述顆粒的形狀：細長度、棱角圓度和粗糙度，如圖10所示。

圖10 道砟顆粒形狀示例

3.1 幾何形狀的評價指標

圖10展示了道砟顆粒的幾何形狀評價指標。對于理想的圓形或等邊多邊形道砟顆粒，其長寬比等于1。對于細長度(EI)，其式如下：

EI=S/L

(1)

式中：S為最小主軸方向上的顆粒直徑；L為最大主軸方向上的顆粒直徑。

此外，棱角圓度(Rd)描述了顆粒表面棱角處接近圓的程度。

(2)

式中：Rinsc是最大內切圓的半徑；Rcirc是位于質心的外接圓的半徑。

粗糙度(Rg)反映了疊加在顆粒表面拐角處的形狀不規則性。粗糙度計算如下。首先，需要分解顆粒輪廓。根據極性系統中的角度θ和極性直徑d，將二維顆粒擴展為與θ和d相對應的連續點。然后，使用“局部加權回歸平滑”(LOESS)方法將這些離散點轉換為平滑曲線/曲面，得到以下式(3)：

(3)

式中：n是二維顆粒輪廓上的點數。yi是點i的y軸坐標。yi-LOESS是平滑后i的y軸坐標。

3.2 統計形狀分析

首先，分析在顆粒庫中獲得的二維顆粒的形狀指標信息。每個形狀指標的統計分布如圖11所示。

由圖11可見，棱角圓度和粗糙度的分布相對狹窄，近似于偏態分布；而細長度的分布較廣，近似于正態分布。

(a)棱角圓度的分布頻率

3.3 建立顆粒庫

為了便于檢索和制樣，將道砟顆粒的形狀指標與坐標一起存儲，存儲格式為[數字，輪廓，細長度，棱角圓度，粗糙度]。

從顆粒庫進行制樣的步驟如下：①輸入所需形狀參數的范圍；②然后，根據這些值的范圍，可以自動搜索出符合要求的二維輪廓；③選擇可視化的二維輪廓并將其輸出到“dxf”文件中。

顆粒庫建立完成后，選擇滿足特定要求的道砟顆粒二維輪廓將非常方便，以方便進行鐵路道砟的離散元建模，從而為后續研究道砟顆粒的幾何信息與力學性能之間的關系奠定基礎。

4 雙軸剪切行為仿真

4.1 雙軸剪切試驗的數值模擬

首先，構建5個具有不同細長度的試樣(EI=0.3、0.45、0.6、0.75、0.9)。試樣的初始大小為9 m×18 m。每個試樣包含大約5 000個顆粒。本文利用線彈性模型表示顆粒之間的接觸法則。通過移動試樣的4個剛性邊界(見圖12)并確保邊界上的力恒定在100 kPa(數值伺服控制的應力邊界條件)來壓縮試樣。為了使試樣在剪切過程中處于準靜態(為了使試樣的加載過程為準靜態加載)，剪切應變率應該非常小。因此，引入慣性參量如下：

(a)初始試樣

(4)

4.2 雙軸剪切試驗結果

4.2.1應力比隨軸向應變的變化

為了研究試樣的剪切強度，本文計算了試樣的平均應力p和偏應力q。在二維雙軸剪切試驗中，有效平均應力p'和偏應力q定義如下：

p'=(σ1+σ2)/2

(5)

q=σ1-σ2

(6)

圖 13(a)展示了在具有不同EI值的試樣中，應力比q/p'隨著軸向應變ε1的變化。初始階段隨著軸向應變的增加，所有試樣的應力比迅速增加，然后達到峰值階段，之后應力比逐漸減小并趨于穩定。當剪切達到穩定階段時，試樣的應力比與EI呈負相關關系。

4.2.2體積應變隨軸向應變的變化

此外，本文研究了剪切過程中不同試樣的體積變化?？紤]到本次模擬的邊界是剛性墻，可以根據剛性邊界的位移近似計算軸向應變和體積應變。

ε1=(h0-h)/h0

εv=(v0-v)/v0

(7)

式中：h0和h為試樣的初始高度和當前時刻的高度；v0和v為試樣的初始體積和當前時刻的體積；負體積應變表示體積膨脹。為了研究穩定后試樣的臨界狀態特性，將所有試樣剪切到軸向應變ε1=30%。在這種形變下，基本滿足臨界狀態的典型條件(即應力和體積隨應變保持恒定)。

如圖 13(b)所示，在初始階段，體積應變隨著軸向應變的增加而減小，而當軸向應變達到大約1%時，體積應變隨著軸向應變的增加而增加。當軸向應變在20%到30%的范圍內達到臨界狀態時，試樣的體積變形保持恒定，這表明所有試樣均表現出由剪切誘發的剪脹變形并伴隨應變軟化。此外，當試樣達到穩定狀態時，EI值越小(EI=0.3、0.45、0.6)，體積應變越大，而EI值越大(EI=0.75，0.9)，體積應變越小。因此，EI值較小的試樣的剪脹性越強。

4.2.3配位數隨軸向應變的變化

在模擬的每個時間步長中，計算每個顆粒的平均配位數(MCN)。圖 13(c)展示了具有不同EI值的試樣的平均配位數隨軸向應變的變化。在剪切初期，所有試樣的平均配位數隨軸向應變增加。當軸向應變大于2%時，平均配位數達到峰值并開始降低。當軸向應變大于10%時，平均配位數基本穩定。每個試樣的平均配位數達到穩定值后，隨著EI值從0.3增加到0.9，平均配位數從8.5減少到5.2，可以發現其與EI值呈負相關關系。

4.2.4滑動接觸的百分比隨軸向應變的變化

滑動接觸的百分比遵循Mohr-Coulomb定律?；瑒酉禂刀x為：

(8)

SP=NSC/NC×100%

(9)

式中：NSC表示顆粒系統中滑動接觸的數量，NC是接觸的總數。

圖 13(d)展示了每個試樣的SP隨軸向應變的變化。所有試樣的SP隨軸向應變的變化基本相同。隨著ε1的增加，SP迅速增加到峰值，然后逐漸下降到穩定值。

圖 13(d)還表明，接觸滑動率隨EI值的增加而降低。當EI值從0.3增加到0.9時，SP的峰值從0.42減小到0.2左右，SP的穩定性值從18%減小到5%左右。

(a)偏應力隨軸向應力的變化

4.2.5EI的摩擦角變化

本節計算了峰值狀態和臨界狀態內摩擦角Φ。非黏性顆粒材料中的Φ定義為：

(10)

式中：σ1和σ2分別為主要和次要主應力。

如圖 14所示，不同EI值的試樣的剪切強度有明顯的變化趨勢。在臨界狀態下，摩擦角與EI值負相關。當EI值從0.3增加到0.9時，臨界摩擦角從0.54減小到0.35左右。此外，峰值摩擦角首先隨EI值的增大而增大，然后減小。

圖14 摩擦角隨細長度的變化

5 結論

本文使用模式識別方法獲得了40 644個鐵路顆粒道砟的二維圖像，計算并存儲了各顆粒的3個形狀指標(細長度，棱角圓度和粗糙度)，建立了二維鐵路道砟顆粒模型庫，進行了形狀指標的統計分析。然后，基于所建立的模型庫，重建了二維道砟顆粒的虛擬試樣。最后，研究了二維道砟顆粒的雙軸剪切特性。結果表明，細長度對道砟顆粒的力學行為具有重要影響。

本文的主要創新如下：

a.提出了模式識別方法，全面獲取道砟顆粒的二維投影的信息。與傳統的方法相比，基于模式識別的方法具有簡單、經濟等優勢。

b.建立了二維道砟顆粒的形態指標模型庫，并通過統計方法進行了統計分析。

c.基于離散元算法研究二維道砟顆粒的力學性能。分析了細長度對剪切性能的影響，為進一步研究形狀指標與道砟顆粒的力學性能之間的定性和定量關系奠定了基礎。