基于真實細觀尺度的瀝青混合料三維重構算法

萬成，張肖寧，賀玲鳳，王端宜，張吉慶

(華南理工大學 土木與交通學院，廣東 廣州，510640)

瀝青混合料具有損傷、非均勻、非連續的特性，其性能既受瀝青、集料和空隙的體積含量影響，也受這些因素的空間分布影響[1]。在大多數的情況下，瀝青混凝土被假設為均勻的，忽略其微細觀結構的影響。造成這種現象的一個主要原因是在非均勻性或微細觀結構的定量觀測和評價方面所遇到的困難。例如，在很多研究中，采用統計數學描述材料非均勻性的方法(如Weibull分布等)[2]，但是，基于統計學的非均勻性和微細觀結構研究其缺點在于其不能準確地反映實際材料的非均勻性和微細觀結構。認識到基于瀝青混合料細觀尺度進行研究的重要性以來，計算機層析攝影技術或計算機層析識別技術(X-ray CT技術)在土木工程領域中的應用和發展，使得研究瀝青混合料內部微細觀結構成為可能。Masad等[3]在1999年采用X-ray CT獲取圓柱形瀝青混合料試件的內部結構組成。從此，瀝青混合料細觀領域的研究成為一個新的發展方向。此后，Wang等[4-6]利用X-ray CT考慮瀝青混合料的實際分布形態，進行混合料內部二維結構的定量或定性分析，特別是 You等[7-8]利用高精度掃描儀采集瀝青混合料試件的內部結構信息后，建立二維數值模型，并采用有限元方法建立二維黏彈力學數值模型，在細觀力學數值模擬方面取得豐富的成果。陳俊[9]依據概率理論，推導瀝青混合料二維截面上粗集料數量級配的計算式，運用計算機隨機投放技術，并綜合考慮集料不規則形狀和級配特征，實現混合料二維數字試件的生成。但目前對混合料內部微細觀研究多為基于二維的識別與分析。由于尚不能完全考慮集料、瀝青及空隙的三維的體積特征，因而其研究結果在多單元體積元統計意義上具有代表性，而在考慮某個單元體時仍然具有一定的隨機性與變異性。瀝青混合料內部結構空間分布狀況的三維重構將是今后研究的重點[10-11]。瀝青混合料三維重構可用于分析混合料的三維體積組成，更可進行細觀力學模擬，模擬各種荷載形式作用下混合料內部微細觀結構的變化與發展。在此過程中，基于空間實際分布形態的真三維重構試樣顯得尤為重要。本文利用X-ray CT掃描實際瀝青混合料試件，直接獲取試樣內部結構，根據瀝青混合料三維重構原理，把握從三維可視化模型向三維數值模型轉化的本質特點，基于X-ray CT技術對瀝青混合料進行三維數值重構算法進行重點研究。

1 三維重構理論

圖像三維重建就是要從二維的圖像序列中提取三維對象的信息，使得用戶能夠直觀地看到三維對象的組織結構，并且加強圖像中原有的各種細節。三維重建的任務是從采樣數據中恢復物體的三維結構，即建立物體的立體模型。重建的方法分為基于外表面的三維重建和直接由體素三維重構2大類。本文為基于體素的三維重構。

體數據的最基本的單元是體素，概念與像素類似。體數據集是定義在三維空間網格上的標量或向量數據，這些網絡通常是正交網格。數據定義在網格節點上，相鄰的8個網格節點構成1個立方體。如圖1所示。假設在空間有一個長為(L-1)×VX，寬為(M-1)×VY，高為(N-1)×VZ的長方體，此長方體的面平行于坐標面。對于節點空間坐標(i, j, k)，i=0, 1, 2, …, L-1；j=0, 1, 2, …, M-1；k=0, 1, 2, …, N-1，式中 VX，VY，VZ為長度、寬度和高度方向的單位長度；i，j，k為編號。體數據與三維數組相對應，大小為L×M×N，定義在此三維空間網格的網格節點上，在節點(i, j, k)處，有體數據值f(i, j, k)[12]。

圖1 體素定義示意圖Fig.1 Sketch map of pixel definition

定義在此三維空間網格的各立方體上的單一體數據，稱之為體素，即體素(i, j, k)有體數據值f(i, j, k)。針對斷層圖像序列而言，即每幅圖像的1個像素被擴展為三維空間的1個勻質的體素，該體素的體數據值f(i, j, k)即為其所對應像素的灰度、顏色索引值或顏色值。可以根據1個算子，從體數據值f(i, j, k)得出體素(i, j, k)所對應的顏色分量C(i, j, k)及非透明度O(i, j,k)。根據上述原理，可知，每個體素由8個相鄰節點組成，實現重建的關鍵技術是求出每個節點的三維坐標(i, j, k)，以及每個體素與自身相鄰8個節點形成的對應關系。

2 三維重構算法

圖2 CT系列連續圖像Fig.2 Series of continuous CT images

利用CT掃描混合料試件后可以得到N幅一系列連續CT圖像，如圖2所示。對圖像進行疊加，利用圖像處理軟件如VGStudio，IDL和Image Pro Plus等軟件實現試樣的三維可視化重構，如圖3所示。必須指出的是三維可視化模型不能用來進行力學數值模擬，開發基于 CT掃描的三維圖像技術來獲取瀝青混合料試件的內部結構，除了用于分析試樣的內部體積組成，更重要的是為進行力學數值模擬，即虛擬力學試驗，這就需要將三維可視化試樣轉換成三維數值試樣。由于沒有專用于瀝青混合料的圖像處理軟件，并且處理過程相對復雜，實際上現有工業CT所附帶的圖像處理軟件主要用于工業生產領域，如工業產品缺陷檢查等。本文重點研究考慮空間真實細觀結構的瀝青混合料三維數值模型的重構算法，根據提出的算法開發瀝青混合料試樣三維重構程序，建立反映空間實際分布形態的瀝青混合料三維數值模型，為虛擬力學試驗奠定基礎。

CT圖像數量可以通過參數設定控制，理論上，只要獲得的CT圖像數量越多，相鄰間距越小，越接近真實混合料試樣。例如，對于一個標準馬歇爾試件，其標準高度為63.5 mm，若獲得1 001幅CT圖像，則相鄰2幅CT圖像之間的間距為0.063 5 mm，這樣小的圖像間距對于建立的三維結構誤差可以忽略不計。CT圖像可以轉換成256級的灰度圖像，即8位黑白圖像。對于每幅CT圖像的灰度像素矩陣，假設其矩陣大小為mn，即像素個數為mn，因此1個像素可以代表1個長方體單元，每個長方體單元具有8個節點，多幅連續CT圖像轉換成的灰度像素矩陣，疊加形成的灰度單元體網格，如圖4所示。取其中任意一幅CT圖像的灰度像素矩陣，如圖 5所示，其 1個像素由1個字節描述。0表示黑色，255為白色。從圖5的像素矩陣示意圖可見：每1個像素點是1個微小的正方形區域，可以用來表征材料的二維平面網格。根據上述原理和重構思想，可以得出結論：實現重建的關鍵技術是求出每個節點的三維空間坐標(i, j, k)，以及每個單元體與自身相鄰8個節點的匹配邏輯關系。

圖3 連續CT圖像實現三維可視化試樣重構Fig.3 3D reconstruction of visualization model with continuous CT images

圖4 多幅CT圖像的灰度單元體網格Fig.4 Gray unit net of multiple CT images

圖5 單幅CT圖像像素矩陣示意圖Fig.5 Sketch map of pixel matrix of one image

2.1 像素單元體的節點編號

取圖4中第1幅CT灰度圖像的像素矩陣，其單元和節點編號示意圖如圖6所示。單元體編號I=1, 2,3, …, mn，節點編號J=1, 2, 3, …, (m+1)(n+1)。單元和節點編號規則為：從矩陣左下角第1個單元、節點起，從左至右、從下往上編號，每一行的末尾編號與下一行的起始單元、節點編號相連，從圖 6可見：1幅CT灰度圖像的像素矩陣的單元數目為mn，節點數目為(m+1)(n+1)。

取圖4中第2幅CT灰度圖像的像素矩陣，其單元和節點編號示意圖如圖7所示，單元、節點編號續接第1幅CT圖像像素矩陣最末1個單元、節點編號，所以該圖的第 1個單元、節點編號分別為 mn+1，(m+1)(n+1)+1，其余編號規則與第1幅相同。

比較第1幅與第2幅示意圖的單元、節點編號可知：相鄰2幅圖片對應位置的單元編號相差mn，節點編號相差(m+1)(n+1)，所以對于圖4中第k幅圖片，其單元和節點編號如圖8所示，為方便示意圖中標記，記 P=(k-1)mn，Q=(k-1)(m+1)(n+1)。

2.2 單元體與節點對應

提取相鄰兩層圖片中相鄰的 8 個節點，構成 1個單元體。根據2.1節中單元和節點編號規則，可方便地用數學描述出每個單元體與自身8個節點之間的對應邏輯關系。

圖6 第1幅CT圖像像素矩陣單元和節點編號示意圖Fig.6 Sketch map of elements and nodes label for the first CT image pixel matrix

圖7 第2幅CT圖像像素矩陣單元和節點編號示意圖Fig.7 Sketch map of elements and nodes label for kth CT image pixel matrix

圖8 第k幅CT圖像像素矩陣單元和節點編號示意圖Fig.8 Sketch map of elements and nodes label for kth CT image pixel matrix

2.3 節點的空間三維坐標

單元體與節點的匹配關系確定以后，如果能確定節點的空間三維坐標，相應的單元體位置也就得到確定。由2.1中節點編號可知，對于任意第k幅圖片，其節點編號在(k-1)(m+1)(n+1)+1到k(m+1)(n+1)之間，設

則對于任意第k幅圖片的節點三維坐標(i, j, k)可以表示為

至此，單元和節點的編號、對應關系已經確立，節點的三維空間坐標也已經確定。混合料的三維數值模型初步建立起來，如圖9所示。必須要指出的是，此時建立的模型包含所有背景像素(CT圖像中的黑色區域)和混合料像素所代表的單元和節點，所以呈現的是立方體形狀，而實際混合料試件是圓柱體形狀。

圖9 三維數值初始模型Fig.9 3D initial digital model

2.4 混合料區域內單元的識別和提取

初始三維數值模型建立以后，需要對屬于混合料區域內的單元進行識別和提取。背景區域內像素灰度值都為0，混合料區域內像素灰度在0～255之間，如果單純以閾值分割的方法進行單元提取，如(0, T1)屬于空隙，[T1, T2)屬于砂膠，[T2, 255]屬于集料，則會把混合料區域內屬于空隙的單元(灰度為 0)遺漏或者把背景區域的像素包括進去。為此，采用指針掃描和閾值分割相結合的辦法，對屬于混合料區域內的單元進行提取。

所謂指針掃描即采用一定半徑的指針，以試樣中心為圓心點，不斷旋轉掃描，根據指定的判別標準，對區域內單元像素進行識別，區分出屬于混合料區域的單元和背景區域的單元，并對區分出來的背景區域的單元像素賦予大于255的灰度，例如1 000，這樣2個區域內的灰度分布完全區分開來，如圖10所示。

2.5 閾值分割

混合料區域內的單元提取出來以后，需要對單元進行識別，區分出集料、砂膠和空隙。能否準確區分3者，閾值T1和T2的選取顯得尤為重要。在有限元法構造實體模型中，最重要的是確定不同材料的幾何邊界，采用的圖象處理技術包括圖像的三值化和不同組分的邊界提取技術。

CT圖像的三值化首先要確定的是區分空隙與瀝青砂膠的閥值T1及區分瀝青砂膠和集料的閥值T2。閥值 T1和 T2的確定是數字轉換過程中的重點問題，較為常用的方法有雙峰法[13]導數法[14]。雙峰法認為：灰度分布圖中，2個峰值之間的應為同一介質；而導數法認為：灰度一階偏導數分布圖中，2個峰值之間的應為同一介質。李曉軍等[15]分別利用雙峰法和導數法進行了試算，結果發現對于CT圖像，這2種方法差別不大，本文采用雙峰法進行了計算。

還可以采用物質密度對比的方法進行閥值分割。分別對集料顆粒與配制的瀝青砂膠進行CT掃描，獲得的圖像中確定2種物質的灰階范圍。

圖10 指針掃描Fig.10 Pointer scanning

初步確定 T1和 T2后，采用反復試驗方法微調閥值T1和T2，圖11為三值化圖CT圖像。圖像中紅色部分分別為集料、砂膠和空隙。

2.6 三維數值試樣的建立

利用開發的 MATLAB程序實現上述過程，最后利用MATLAB中的DLWRITE命令將數據文件寫入構造有限元模型的INP文件，該INP文件包含了連續斷層圖像的細觀結構信息，最終建立三維細觀模型，如圖12所示，區別于傳統的均質體有限元數值模型，該模型真實地反映實際試樣中集料、砂膠和空隙的三維細觀分布，即考慮3者的空間真實分布形態，只需對3者分別賦予不同的材料屬性，為以后的瀝青混合料虛擬力學試驗奠定基礎。

圖11 3組分的分割Fig.11 Image segmentation of three constituents

圖12 三維數值模型Fig.12 3D digital specimen

3 結論

(1) 對基于X-ray CT技術的瀝青混合料三維重構方法進行研究。分析瀝青混合料三維重構原理，提出了對瀝青混合料進行三維數值重構的基本算法，開發瀝青混合料試樣三維重構程序，該數值試樣區別于傳統的均質體有限元模型，而是基于集料、瀝青砂膠和空隙的真實三維空間分布，因而更加貼近實際情況。

(2) 基于瀝青混合料真實細觀尺度建立的三維數值模型，是整個瀝青混合料虛擬力學試驗的前提與基礎，有助于以后從三維細觀尺度研究瀝青混合料的力學行為。

(3) 限于計算機運算速度限制，對 CT 圖像只能采取降低分辨率的方法，從而影響到了集料、砂膠和空隙的分割精度，隨著計算機技術的快速發展以及算法的優化，可以預見將獲得更為高精度的三維數值模型，但本研究從方法上解決基于真實三維細觀分布的瀝青混合料重構問題。

[1] TAN Yi-qiu, XU Hui-ning, LI Xiao-min. Is normal distribution the most appropriate statistical distribution for volumetric properties in asphalt mixtures?[J]. Journal of Testing and Evaluation, 2009, 27(5): 458-478.

[2] 邱延峻, 閆常赫, 艾長發. 非均質瀝青混合料劈裂試驗全過程數值模擬[J]. 交通運輸工程學報, 2009, 9(2): 12-16.QIU Yan-jun, YAN Chang-he, AI Chang-fa. Numerical simulation of split test process for asphalt mixture under heterogeneous state[J]. Journal of Traffic and Transportation Engineering, 2009, 9(2): 12-16.

[3] Masad E, Muhunthan B, Shashidhar N, et al. Quantifying laboratory compaction effects on the internal structure of asphalt concrete[C]//Proceeding of Transportation Research Record(TRB) 78th Annual Meeting. Washington D C: Transportation Research Board, 1999: 179-185.

[4] Wang L B, Wang Y, Mohammad L, et al. Voids distribution and performance of asphalt concrete[J]. International Journal of Pavements, 2002, 1(3): 22-33.

[5] Tashman L, Masad E, Angelo J, et al. X-ray tomography to characterize air void distribution in superpave gyratory compacted specimens[J]. International Journal of Pavement Engineering, 2002, 3(1): 19-28.

[6] Kim H, Wagoner M P, Buttlar W G. Numerical fracture analysis on the specimen size dependency of asphalt concrete using a cohesive softening model[J]. Construction and Building Materials, 2009, 23(5): 2112-2120.

[7] You Z P, Dai Q L. Dynamic complex modulus predictions of hot-mix asphalt using a micromechanical-based finite element model[J]. Canadian Journal of Civil Engineering, 2007, 34(12):1519-1528.

[8] Dai Q L. Prediction of dynamic modulus and phase angle of stone-based composites using a micromechanical finite-element approach[J]. Journal of Materials in Civil Engineering, 2010,22(6): 618-627.

[9] 陳俊. 基于離散元方法的瀝青混合料虛擬疲勞試驗研究[D].南京: 東南大學交通學院, 2010: 101-120.CHENG Jun. Virtual fatigue tests of asphalt mixture based on discrete element method[D]. Nanjing: Southeast University.School of Transportation, 2010: 101-120.

[10] 汪海年, 郝培文. 瀝青混合料微細觀結構的研究進展[J]. 長安大學學報: 自然科學版, 2008, 28(3): 11-14.WANG Hai-nian, HAO Pei-wen. Advances in microstructure study on asphalt mixture[J]. Journal of Changan University:Natural Science Edition, 2008, 28(3): 11-14.

[11] You Z P, Dai Q L. Review of advances in micromechanical modeling of aggregate–aggregate interactions in asphalt mixtures[J]. Canadian Journal of Civil Engineering, 2007, 34(2):239-252.

[12] 王旭, 趙歆波, 桂業英, 等. 醫學圖像序列的直接體繪制算法研究[J]. 生物醫學工程學雜志, 1999, 16(3): 324-328.WANG Xu, ZhAO Xin-bo, GUI Ye-ying, et al. A new algorithm for direct volume rendering of medic image series[J]. Journal of Biomedical Engineering, 1999, 16(3): 324-328.

[13] Kranz R L. Crack-crack and crack-pore interactions in stressed granite[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1979, 16(1): 37-47.

[14] Wu F J, Thomsen L. Microfracturing and deformation of westerly granite under creep condition[J]. International Journal of Rock Mechanics and Mining Sciences, 1975, 12(5): 167-173.

[15] 李曉軍, 張肖寧, 武建民. 瀝青混合料單軸重復加卸載破損CT識別[J]. 哈爾濱工業大學學報, 2005, 37(9): 1228-1230.LI Xiao-jun, ZHANG Xiao-ning, WU Jian-min. Internal structure discrimination of asphalt mix with X-ray computerized tomography under repetitive uniaxial loadings[J]. Journal of Harbin Institute of Technology, 2005, 37(9): 1228-1230.