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

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

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

1 三維重構(gòu)理論

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

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

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

定義在此三維空間網(wǎng)格的各立方體上的單一體數(shù)據(jù)，稱之為體素，即體素(i, j, k)有體數(shù)據(jù)值f(i, j, k)。針對(duì)斷層圖像序列而言，即每幅圖像的1個(gè)像素被擴(kuò)展為三維空間的1個(gè)勻質(zhì)的體素，該體素的體數(shù)據(jù)值f(i, j, k)即為其所對(duì)應(yīng)像素的灰度、顏色索引值或顏色值?？梢愿鶕?jù)1個(gè)算子，從體數(shù)據(jù)值f(i, j, k)得出體素(i, j, k)所對(duì)應(yīng)的顏色分量C(i, j, k)及非透明度O(i, j,k)。根據(jù)上述原理，可知，每個(gè)體素由8個(gè)相鄰節(jié)點(diǎn)組成，實(shí)現(xiàn)重建的關(guān)鍵技術(shù)是求出每個(gè)節(jié)點(diǎn)的三維坐標(biāo)(i, j, k)，以及每個(gè)體素與自身相鄰8個(gè)節(jié)點(diǎn)形成的對(duì)應(yīng)關(guān)系。

2 三維重構(gòu)算法

圖2 CT系列連續(xù)圖像Fig.2 Series of continuous CT images

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

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

圖3 連續(xù)CT圖像實(shí)現(xiàn)三維可視化試樣重構(gòu)Fig.3 3D reconstruction of visualization model with continuous CT images

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

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

2.1 像素單元體的節(jié)點(diǎn)編號(hào)

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

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

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

2.2 單元體與節(jié)點(diǎn)對(duì)應(yīng)

提取相鄰兩層圖片中相鄰的 8 個(gè)節(jié)點(diǎn)，構(gòu)成 1個(gè)單元體。根據(jù)2.1節(jié)中單元和節(jié)點(diǎn)編號(hào)規(guī)則，可方便地用數(shù)學(xué)描述出每個(gè)單元體與自身8個(gè)節(jié)點(diǎn)之間的對(duì)應(yīng)邏輯關(guān)系。

圖6 第1幅CT圖像像素矩陣單元和節(jié)點(diǎn)編號(hào)示意圖Fig.6 Sketch map of elements and nodes label for the first CT image pixel matrix

圖7 第2幅CT圖像像素矩陣單元和節(jié)點(diǎn)編號(hào)示意圖Fig.7 Sketch map of elements and nodes label for kth CT image pixel matrix

圖8 第k幅CT圖像像素矩陣單元和節(jié)點(diǎn)編號(hào)示意圖Fig.8 Sketch map of elements and nodes label for kth CT image pixel matrix

2.3 節(jié)點(diǎn)的空間三維坐標(biāo)

單元體與節(jié)點(diǎn)的匹配關(guān)系確定以后，如果能確定節(jié)點(diǎn)的空間三維坐標(biāo)，相應(yīng)的單元體位置也就得到確定。由2.1中節(jié)點(diǎn)編號(hào)可知，對(duì)于任意第k幅圖片，其節(jié)點(diǎn)編號(hào)在(k-1)(m+1)(n+1)+1到k(m+1)(n+1)之間，設(shè)

則對(duì)于任意第k幅圖片的節(jié)點(diǎn)三維坐標(biāo)(i, j, k)可以表示為

至此，單元和節(jié)點(diǎn)的編號(hào)、對(duì)應(yīng)關(guān)系已經(jīng)確立，節(jié)點(diǎn)的三維空間坐標(biāo)也已經(jīng)確定。混合料的三維數(shù)值模型初步建立起來，如圖9所示。必須要指出的是，此時(shí)建立的模型包含所有背景像素(CT圖像中的黑色區(qū)域)和混合料像素所代表的單元和節(jié)點(diǎn)，所以呈現(xiàn)的是立方體形狀，而實(shí)際混合料試件是圓柱體形狀。

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

2.4 混合料區(qū)域內(nèi)單元的識(shí)別和提取

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

所謂指針掃描即采用一定半徑的指針，以試樣中心為圓心點(diǎn)，不斷旋轉(zhuǎn)掃描，根據(jù)指定的判別標(biāo)準(zhǔn)，對(duì)區(qū)域內(nèi)單元像素進(jìn)行識(shí)別，區(qū)分出屬于混合料區(qū)域的單元和背景區(qū)域的單元，并對(duì)區(qū)分出來的背景區(qū)域的單元像素賦予大于255的灰度，例如1 000，這樣2個(gè)區(qū)域內(nèi)的灰度分布完全區(qū)分開來，如圖10所示。

2.5 閾值分割

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

CT圖像的三值化首先要確定的是區(qū)分空隙與瀝青砂膠的閥值T1及區(qū)分瀝青砂膠和集料的閥值T2。閥值 T1和 T2的確定是數(shù)字轉(zhuǎn)換過程中的重點(diǎn)問題，較為常用的方法有雙峰法[13]導(dǎo)數(shù)法[14]。雙峰法認(rèn)為：灰度分布圖中，2個(gè)峰值之間的應(yīng)為同一介質(zhì)；而導(dǎo)數(shù)法認(rèn)為：灰度一階偏導(dǎo)數(shù)分布圖中，2個(gè)峰值之間的應(yīng)為同一介質(zhì)。李曉軍等[15]分別利用雙峰法和導(dǎo)數(shù)法進(jìn)行了試算，結(jié)果發(fā)現(xiàn)對(duì)于CT圖像，這2種方法差別不大，本文采用雙峰法進(jìn)行了計(jì)算。

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

圖10 指針掃描Fig.10 Pointer scanning

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

2.6 三維數(shù)值試樣的建立

利用開發(fā)的 MATLAB程序?qū)崿F(xiàn)上述過程，最后利用MATLAB中的DLWRITE命令將數(shù)據(jù)文件寫入構(gòu)造有限元模型的INP文件，該INP文件包含了連續(xù)斷層圖像的細(xì)觀結(jié)構(gòu)信息，最終建立三維細(xì)觀模型，如圖12所示，區(qū)別于傳統(tǒng)的均質(zhì)體有限元數(shù)值模型，該模型真實(shí)地反映實(shí)際試樣中集料、砂膠和空隙的三維細(xì)觀分布，即考慮3者的空間真實(shí)分布形態(tài)，只需對(duì)3者分別賦予不同的材料屬性，為以后的瀝青混合料虛擬力學(xué)試驗(yàn)奠定基礎(chǔ)。

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

圖12 三維數(shù)值模型Fig.12 3D digital specimen

3 結(jié)論

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

(2) 基于瀝青混合料真實(shí)細(xì)觀尺度建立的三維數(shù)值模型，是整個(gè)瀝青混合料虛擬力學(xué)試驗(yàn)的前提與基礎(chǔ)，有助于以后從三維細(xì)觀尺度研究瀝青混合料的力學(xué)行為。

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

[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] 邱延峻, 閆常赫, 艾長(zhǎng)發(fā). 非均質(zhì)瀝青混合料劈裂試驗(yàn)全過程數(shù)值模擬[J]. 交通運(yùn)輸工程學(xué)報(bào), 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] 陳俊. 基于離散元方法的瀝青混合料虛擬疲勞試驗(yàn)研究[D].南京: 東南大學(xué)交通學(xué)院, 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] 汪海年, 郝培文. 瀝青混合料微細(xì)觀結(jié)構(gòu)的研究進(jìn)展[J]. 長(zhǎng)安大學(xué)學(xué)報(bào): 自然科學(xué)版, 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] 王旭, 趙歆波, 桂業(yè)英, 等. 醫(yī)學(xué)圖像序列的直接體繪制算法研究[J]. 生物醫(yī)學(xué)工程學(xué)雜志, 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] 李曉軍, 張肖寧, 武建民. 瀝青混合料單軸重復(fù)加卸載破損CT識(shí)別[J]. 哈爾濱工業(yè)大學(xué)學(xué)報(bào), 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.