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

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

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

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

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

1 搭建消影拍照臺(tái)

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

圖1 消影拍照臺(tái)

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ū)別。針對(duì)二維廓形的研究，需考察道砟顆粒具有代表性角度下拍攝的圖像。使道砟顆粒從距地面0.5 m自由落下，待其在平坦地面上穩(wěn)定后，將其移至拍照臺(tái)進(jìn)行拍攝。此后運(yùn)用MATLAB標(biāo)定工具箱獲得了相機(jī)的內(nèi)外參數(shù)，用以計(jì)算得出道砟顆粒的真實(shí)尺寸。

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

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

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

算法主要包括顏色空間轉(zhuǎn)換與二值化處理、獲取道砟邊界、初始化、追蹤循環(huán)等過(guò)程。其流程如圖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ù)取值通過(guò)試驗(yàn)獲得。

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

3.2 追蹤提取二維道砟廓形

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

3.2.1 初始化

(1)尋找起始點(diǎn)

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

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

(1)

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

(2)初次擬合

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

圖3 初次擬合示意

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

(2)

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

(3)

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

(4)

(5)

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

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

圖4 步進(jìn)示意

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

3.2.2 追蹤

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

(1)擬合

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

(2)步進(jìn)

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

α3=α2+180°

(6)

圖5 α2和α3示意

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

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

(7)

(3)校準(zhǔn)

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

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

(8)

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

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

(9)

其中，δ表示單位方向向量，則定義移動(dòng)的方向。

圖6 校準(zhǔn)示意

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

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

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

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

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

(10)

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

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

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

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

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

步驟2：定義角度記錄變量β，用于存儲(chǔ)角度累計(jì)值。

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

步驟4：選取任意一點(diǎn)為迭代初始點(diǎn)，計(jì)算輸入的廓形點(diǎn)列中相鄰兩個(gè)點(diǎn)之間的角度βi

(11)

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

(12)

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

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

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

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

(13)

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

(14)

4 實(shí)例分析及可靠性驗(yàn)證

4.1 實(shí)例分析

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

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

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

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

圖8 實(shí)例分析

4.2 算法可靠性驗(yàn)證

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

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

(15)

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

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

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

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

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

(16)

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

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

(17)

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

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

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

5 結(jié)論

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

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

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

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

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