基于圖像處理的道砟顆粒二維廓形追蹤提取算法

李 朋,劉鑒興,肖杰靈,王 平

(1.西南交通大學(xué)高速鐵路線路工程教育部重點實驗室,成都 610031； 2.西南交通大學(xué)土木工程學(xué)院,成都 610031)

隨著我國鐵路提速與高速化發(fā)展，有砟軌道道床病害日益突出。通過離散元法來模擬有砟道床細(xì)觀、宏觀力學(xué)行為與劣化機(jī)制，是發(fā)展和完善高速鐵路有砟軌道結(jié)構(gòu)設(shè)計與運營維護(hù)技術(shù)的一項重要研究手段[1-3]。其中，精確提取道砟顆粒廓形是保障仿真結(jié)果可靠性的重要基礎(chǔ)。除服務(wù)于離散元仿真分析外，準(zhǔn)確的道砟顆粒廓形還可以應(yīng)用于評估散體道床道砟堆積體的物理特性，如堆積密實度、道砟棱角系數(shù)、道床級配等[4]。

在提取顆粒的幾何形態(tài)時，通常采用數(shù)碼相機(jī)、激光掃描或CT設(shè)備等計算機(jī)視覺測量技術(shù)。其中，由于數(shù)碼相機(jī)的技術(shù)復(fù)雜度及成本低于激光掃描等方式，因此其被廣泛用于獲取顆粒的二維特征[5-6]。然而，在運用數(shù)碼相機(jī)拍攝的數(shù)字圖像來提取道砟顆粒二維廓形時，使用的傳統(tǒng)邊緣檢測方法，如Sobel算子、Roberts算子等雖然具有實現(xiàn)簡單，運算速度快等優(yōu)點，但其抗干擾能力差，邊緣寬度往往比實際大，且邊緣像素孤立或分小段連續(xù)[7-9]。不僅降低了道砟顆粒廓形的提取精度，而且孤立和分?jǐn)嗟臒o序點云也不利于后期道砟形態(tài)表征的評估。一些學(xué)者提出了不同的邊界檢測改良算法[10-12]，但這些算法大都針對特定的情況和物體，完全套用這些算法，亦不能獲得理想的道砟顆粒二維廓形。

為提高道砟顆粒二維廓形的連續(xù)性和準(zhǔn)確性，同時得到有序的廓形數(shù)據(jù)點列以彌補(bǔ)無序數(shù)據(jù)點云在分析道砟形態(tài)特征時的不足，本文針對性地提出一種基于圖像處理的道砟二維廓形追蹤提取算法。算法特點在于抗外界環(huán)境干擾能力高，能夠一次性處理多個道砟顆粒，提取的道砟顆粒二維廓形不僅連續(xù)性好，而且檢測出來的邊緣是單像素寬。通過“追蹤”的思想，使提取的廓形數(shù)據(jù)形成有序點列。通過將算法提取的廓形數(shù)據(jù)與游標(biāo)卡尺測量的實際道砟廓形數(shù)據(jù)對比分析，發(fā)現(xiàn)兩種粒徑的差值均在1 mm之內(nèi)?？紤]測量誤差，可以肯定本文算法的有效性與可靠性。

1 搭建消影拍照臺

隨意地將道砟放置于多變的自然環(huán)境中，由于光照條件不易控制，會有很多復(fù)雜的噪聲影響，降低圖像的拍攝質(zhì)量。為減小不良因素對圖像的影響，同時避免繁瑣的圖像預(yù)處理步驟，需要搭建一個特殊的消影拍照臺，該拍照臺由數(shù)碼相機(jī)、白色背景板、透明玻璃板和LED燈組等構(gòu)件組成，如圖1所示。

圖1 消影拍照臺

2 標(biāo)定與圖像預(yù)處理

在拍攝道砟顆粒的二維圖像之前，要進(jìn)行標(biāo)定以獲得相機(jī)的內(nèi)外參數(shù)，進(jìn)而獲得道砟顆粒的尺寸等信息。本文采用了經(jīng)典的棋盤格標(biāo)定法進(jìn)行標(biāo)定。首先制作了棋盤格標(biāo)定板，并拍攝了20張其不同角度的圖像。不同角度拍攝所得道砟顆粒的圖像有較大區(qū)別。針對二維廓形的研究，需考察道砟顆粒具有代表性角度下拍攝的圖像。使道砟顆粒從距地面0.5 m自由落下，待其在平坦地面上穩(wěn)定后，將其移至拍照臺進(jìn)行拍攝。此后運用MATLAB標(biāo)定工具箱獲得了相機(jī)的內(nèi)外參數(shù)，用以計算得出道砟顆粒的真實尺寸。

拍攝完標(biāo)定板后，將6個道砟顆粒均布于白色背景板上(后期實踐證明，選擇與道砟顆粒顏色差異大的背景色即可，如綠色、藍(lán)色等)，每批次擺放的道砟顆粒大致在同一位置。

拍攝的道砟顆粒圖像無需進(jìn)行圖像預(yù)處理，而且由于拍攝時使用的背景與道砟顆粒對比度很大，也不需要進(jìn)行灰度變換。如有特殊情況，本文推薦采用的主要圖像預(yù)處理方法有：圖像平滑、圖像去噪、圖像銳化等[13]。進(jìn)行預(yù)處理的目的是對不符合要求的圖像進(jìn)行某些適當(dāng)?shù)淖儞Q，從而使獲得的道砟顆粒圖像清晰并且滿足后期處理的要求。

3 道砟顆粒二維廓形追蹤提取算法

算法主要包括顏色空間轉(zhuǎn)換與二值化處理、獲取道砟邊界、初始化、追蹤循環(huán)等過程。其流程如圖2所示。

圖2 算法流程

3.1 顏色空間的轉(zhuǎn)換與二值化處理

為了方便分辨圖片中的像素是屬于背景色還是屬于道砟顆粒，注意到明度最能反映道砟顆粒與背景的差異，且RGB顏色空間不方便描述顏色之間的距離，因而考慮將RGB顏色空間轉(zhuǎn)化到Lab顏色空間[14-16]。兩者之間的轉(zhuǎn)化用XYZ顏色空間傳遞，即先將RGB顏色空間轉(zhuǎn)化到XYZ顏色空間，再將XYZ顏色空間轉(zhuǎn)化到Lab顏色空間。本文采用一種通用的空間轉(zhuǎn)化算法[17]，具體參數(shù)取值通過試驗獲得。

在將RGB顏色空間轉(zhuǎn)化到Lab顏色空間后，設(shè)置Lab顏色空間的明度(即L、a、b中的L矩陣)閾值，將提取的明度矩陣轉(zhuǎn)化為二值矩陣，用來表征道砟顆粒與背景色?；谏鲜霾僮鳎煽焖偻ㄟ^梯度算子求出二值矩陣的梯度，從而得到二維圖像中道砟顆粒的邊界像素點集。

3.2 追蹤提取二維道砟廓形

在一幅二維圖像中為每一個道砟選取一個位于道砟顆粒輪廓線以內(nèi)的點作為起點(該點應(yīng)大致位于道砟顆粒中心)。因前文要求每批次道砟顆粒擺放位置大致相同，所以后續(xù)不同批次的道砟顆粒都可以用第一次選點作為起點，這樣可提高數(shù)據(jù)處理的速度。

3.2.1 初始化

(1)尋找起始點

用50×50的像素矩陣(本文結(jié)合二維相片中道砟顆粒的邊界寬度，選取50×50的像素矩陣?？筛鶕?jù)實際情況調(diào)整，如設(shè)置40×40或60×60等不同形式的像素矩陣)從起點向正上方移動，找到道砟顆粒的邊界，即判斷像素矩陣框內(nèi)是否存在道砟顆粒邊界像素點。具體方法如下：

在Lab顏色空間中，邊界像素點的值為非0數(shù)(一般接近1)，其余空白像素點部分的值均為0。所以在50×50的像素矩陣范圍內(nèi)求和

(1)

其中，Ni為第i個像素點的值。當(dāng)C非零時，就可以說明在此范圍內(nèi)存在至少一個像素點，從而將矩陣中心移動至此，該像素點稱為0號中心像素點；當(dāng)C為零時，則說明在此范圍內(nèi)不存在像素點，則繼續(xù)移動像素框。

(2)初次擬合

此步驟為追蹤的關(guān)鍵，目的在于運用最小二乘法原理確定像素矩陣沿道砟廓形的移動方向，如圖3所示。

圖3 初次擬合示意

已知通過顏色空間轉(zhuǎn)換和梯度算子求得的道砟顆粒邊界像素點集為

(2)

將像素矩陣中的道砟顆粒邊界像素點集令為[A]，即

(3)

設(shè)直線擬合系數(shù)為(k，m)，直線表達(dá)式為y=kx+m，則有

(4)

(5)

(3)初次步進(jìn)

在求得直線后，以10個像素為步長(為保障提取精度及效率，10個像素為宜)將中心像素點按照直線任意方向移動(“任意”只適用于初次步進(jìn))，此時選取的方向線與水平正方向形成夾角，如圖4所示為正方向情況。

圖4 步進(jìn)示意

移動以后判斷像素矩陣框內(nèi)是否存在道砟顆粒邊界像素點，若存在，則將50×50的像素矩陣中心移至距離最近的道砟顆粒邊界像素點上。

3.2.2 追蹤

此步驟是本文算法的核心，與初始化步驟中的內(nèi)容大致相同，但需做一些改動，具體內(nèi)容如下。

(1)擬合

運用最小二乘法原理，求出移動后新像素矩陣內(nèi)廓形像素點的擬合直線，根據(jù)此直線確定下一次步進(jìn)方向。

(2)步進(jìn)

此時需保證步進(jìn)方向與初次步進(jìn)時相同(順時針或逆時針)，但確定步進(jìn)方向時的線為直線而非射線，會有兩個相反方向。觀察到像素框移動步長小時，道砟顆粒的二維輪廓梯度改變量不會太大，所以前一次產(chǎn)生的“α1”和本次擬合直線與水平正方向產(chǎn)生的夾角“α2”應(yīng)相差不會太大，即如圖5所示，其中

α3=α2+180°

(6)

圖5 α2和α3示意

基于上述條件，可設(shè)定夾角參數(shù)θ來判斷50×50的像素矩陣沿順時針(或逆時針)運動的方向。其中

θ=|α1-α2|<180°

(7)

(3)校準(zhǔn)

步進(jìn)后，若在50×50的像素矩陣范圍內(nèi)存在道砟顆粒邊界像素點，則將像素矩陣中心移至距離最近的道砟顆粒邊界像素點上，如圖6所示。移動到距離最近點的方法如下。

設(shè)50×50的像素矩陣中心所在位置坐標(biāo)為(a，b)，像素矩陣范圍內(nèi)存在的任意一個道砟顆粒邊界像素點的坐標(biāo)為(xj，yj)，兩點之間的距離可以表示為

(8)

利用計算機(jī)可以快速計算出minLj，則校準(zhǔn)后的點可表示為

(a′，b′)=(a，b)+minLj·δ

(9)

其中，δ表示單位方向向量，則定義移動的方向。

圖6 校準(zhǔn)示意

3.2.3 輸出有序二維廓形數(shù)據(jù)點列

進(jìn)行下一次移動，每次移動后判斷新的邊界位置是否與起點邊界位置重合或接近(小于5像素)，即是否回到0號中心像素點，如果回到0號中心像素點則停止移動并輸出所記錄的點列，將每一次的中心像素點坐標(biāo)連線，繪制擬合的道砟顆粒二維廓形圖；否則繼續(xù)追蹤。在第一顆道砟廓形尋找完成后，由于前期將各批道砟均擺放在相同位置，便可定義各顆道砟50×50的像素矩陣的初始點，由此循環(huán)上述“初始化—追蹤”步驟，即可一次性獲取6顆道砟的二維廓形。

3.3 二維廓形點列數(shù)據(jù)優(yōu)化

通過上述廓形提取算法獲取道砟顆粒二維廓形表示為

C={(xi，yi)|i=1，2，…，n}

(10)

n的取值根據(jù)道砟大小及相機(jī)精度不同而變化，本文中n在100至300內(nèi)波動。

如圖7所示，由于直接獲取的廓形數(shù)據(jù)由大量密集點列組成，而在一定范圍內(nèi)，大多數(shù)點對于廓形變化的表征作用并不明顯，造成了點列的冗雜。點列數(shù)據(jù)越多，對后期仿真效率的影響也就越大，因此針對一些不需要精細(xì)道砟顆粒二維廓形的情況，可對數(shù)據(jù)進(jìn)行壓縮優(yōu)化。優(yōu)化的目的是在不影響道砟顆粒整體形狀的情況下，對邊界節(jié)點個數(shù)進(jìn)行精簡，也即忽略部分紋理、變化率較小的表面波動，而將道砟顆粒的形狀表達(dá)出來。下面說明優(yōu)化方法。

圖7 數(shù)據(jù)優(yōu)化示意

設(shè)置一個角度變化閾值，每積累超過該閾值則引入一個新的點，在小于該閾值時則減少點數(shù)以省略不關(guān)注的細(xì)節(jié)。以此達(dá)到用最優(yōu)的方式記錄廓形邊界數(shù)據(jù)的目的。

步驟1：設(shè)定一個角度控制閾值，該閾值用于表征邊界廓形的變化劇烈程度。閾值越小，對廓形的描述越精確，但廓形數(shù)據(jù)的數(shù)量就相對較多；角度越大，對廓形的描述則越粗糙，廓形數(shù)據(jù)的數(shù)量就相對較少。

步驟2：定義角度記錄變量β，用于存儲角度累計值。

步驟3：定義角度增量βi，用于表示相鄰兩點之間的角度增量值；初始β=0。

步驟4：選取任意一點為迭代初始點，計算輸入的廓形點列中相鄰兩個點之間的角度βi

(11)

判斷β=β+βi是否超過角度控制閾值ρ

(12)

超過則記錄當(dāng)前節(jié)點坐標(biāo)，并重置β=0；否則累計角度記錄變量β，再進(jìn)行下一個點的計算。

步驟5：每次計算后判斷是否已經(jīng)遍歷了所有節(jié)點，即所計算的點與迭代初始點坐標(biāo)信息是否相同。相同則停止計算，輸出記錄點，形成新的多邊形廓形，否則繼續(xù)循環(huán)此計算。

步驟6：設(shè)求出的經(jīng)過壓縮的點為

P={(xj，yj)|k1，k2，k3,…，kN}

(13)

則在遍歷處理前的所有點后，所求得的解的相鄰兩個點之間應(yīng)該滿足下列條件

(14)

4 實例分析及可靠性驗證

4.1 實例分析

拍攝照片：本文選用了遮光紙箱、鋁合金架、數(shù)碼相機(jī)、A3大小的白紙、透明玻璃板以及LED光源等組成消影拍照臺。利用藍(lán)牙遙控器控制數(shù)碼相機(jī)拍攝道砟顆粒的照片。

顏色空間轉(zhuǎn)換：將照片信息導(dǎo)入Matlab軟件，進(jìn)行從RGB顏色空間到Lab顏色空間的轉(zhuǎn)化并求取廓形像素點集。

追蹤提取二維廓形：根據(jù)算法編寫程序?qū)D(zhuǎn)換后的圖片進(jìn)行道砟顆粒二維廓形追蹤提取工作。圖8中藍(lán)色的閉合曲線為提取的6個道砟顆粒二維廓形?？梢钥闯鎏崛〉睦闻c道砟顆粒擬合度良好。

二維廓形點列數(shù)據(jù)優(yōu)化：為方便敘述，現(xiàn)用提取的二維廓形圖中的第一行第一列道砟顆粒輪廓單獨說明。編寫程序，對提取的道砟顆粒二維廓形進(jìn)行數(shù)據(jù)優(yōu)化處理。圖中紅線為優(yōu)化廓形結(jié)果，紅色圓圈為最終保留的輪廓數(shù)據(jù)。

圖8 實例分析

4.2 算法可靠性驗證

除通過觀察提取的廓形與道砟顆粒外形擬合度來了解算法可靠性外，本文在像素坐標(biāo)系下獲得道砟顆粒廓形數(shù)據(jù)后，通過標(biāo)定得到的相關(guān)信息，反算出廓形在真實坐標(biāo)系下的尺寸，再利用游標(biāo)卡尺量取真實顆粒的最大粒徑與反算結(jié)果進(jìn)行比較，旁證本文算法的有效性和可靠性。

通過上述二維廓形追蹤提取算法得到圖8示例顆粒的輪廓數(shù)據(jù)點共186組，同時求得其中兩點之間最大距離為1 007.542 1像素距離。由標(biāo)定反算得到像素距離與真實距離之間的關(guān)系為

(15)

由此可以得到像素坐標(biāo)系下最大粒徑反算結(jié)果為

L=距離×S=1 007.542 1×0.094 0=94.709 0 mm游標(biāo)卡尺測得的道砟真實粒徑為95.26 mm。通過對比，可知獲得的道砟廓形結(jié)果與真實粒徑相差0.551 0 mm。

隨后用游標(biāo)卡尺測量了多個道砟的最大粒徑，兩種粒徑的差值均在1 mm之內(nèi)，考慮到測量誤差的存在，可以認(rèn)為提取的道砟顆粒廓形數(shù)據(jù)有效可靠。部分測量結(jié)果和誤差分析結(jié)果如表1所示。

此外，棱角指數(shù)[18]是定量描述顆粒幾何特征的最為關(guān)鍵的參數(shù)，本文采用Rao C等[19]提出的方法計算顆粒的棱角指數(shù)?；趫D8所提取出的實際廓形，將其輪廓簡化為n邊形。把區(qū)間[0°，180°]分為18個長度為10°的等值區(qū)間，計算n邊形相鄰頂點對應(yīng)的內(nèi)角變化值在這18個區(qū)間上的概率分布，按下式計算棱角度A

表1 部分顆粒測量結(jié)果及誤差分析 mm

(16)

式中，e為區(qū)間E=[e，e+10]的起點值；P(e)為多邊形相鄰內(nèi)角變化值α=|αi+1αi|位于區(qū)間E內(nèi)概率。

用式(17)求三個投影圖棱角度的加權(quán)平均值，即為棱角指數(shù)AI。

(17)

式中，k為投影圖編號；S(k)為投影圖k的面積。

表2列出了Masad E等[20]測定的典型碎石的棱角指數(shù)值和本文中所采取算法計算的棱角指數(shù)，可看出同種碎石的棱角指數(shù)相近，進(jìn)一步旁證本文算法的有效性和可靠性。

表2 典型碎石的棱角指數(shù)范圍

5 結(jié)論

針對道砟顆粒的二維廓形，提出了基于圖像識別的道砟顆粒二維廓形追蹤提取算法。通過實例，分析了算法的精度及耗時，得出以下結(jié)論。

(1)RGB顏色空間不方便描述顏色之間的距離，Lab顏色空間更加便于分辨圖片中的像素是屬于背景色還是屬于道砟顆粒，且明度最能反映道砟顆粒與背景的差異。

(2)改變像素矩陣和移動步長的大小等參數(shù)，可以對道砟顆粒二維廓形數(shù)據(jù)提取的精度和時間進(jìn)行控制。精度越高，則耗時越多，反正亦然。

(3)本算法提取的二維廓形數(shù)據(jù)為有序的點列，而非無序的點云數(shù)據(jù)。

(4)提取出來的廓形與實際廓形的最大粒徑在考慮測量誤差及精度的條件下二者數(shù)值基本相同。本算法適用于提取道砟顆粒的二維廓形。