南海鈣質(zhì)砂微觀構(gòu)造三維圖像分析與表征

萬佳怡，孔德瓊，蘇思楊，彭宇，朱斌,3

(1.浙江大學(xué) 超重力研究中心，浙江 杭州，310058；2.浙江大學(xué) 建筑工程學(xué)院，浙江 杭州，310058；3.浙江大學(xué) 巖土工程研究所，浙江 杭州，310058)

隨著國家“一帶一路”倡議的深入開展，在珊瑚礁區(qū)開展的海洋油氣、漁業(yè)資源開發(fā)和國防建設(shè)過程規(guī)模不斷增大。鈣質(zhì)砂作為南海島礁的主要沉積物，主要礦物成分為碳酸鈣，其沉積過程大多未經(jīng)長途搬運因而保留有豐富的原生生物骨架細(xì)小孔隙[1]，表現(xiàn)出形狀復(fù)雜、內(nèi)孔隙豐富、易破碎、易膠結(jié)[2]的特征，與一般陸相、海相沉積物存在顯著差異。

國內(nèi)外學(xué)者針對鈣質(zhì)砂力學(xué)性質(zhì)進行了大量試驗研究。WU等[3]通過室內(nèi)三軸試驗定量分析了圍壓、顆粒級配、初始密度、排水條件和含水量等因素對顆粒破碎總量的影響。JAⅤDANIAN等[4-5]對鈣質(zhì)砂顆粒進行動三軸試驗，分析了剪切模量、剪切剛度和阻尼比等動態(tài)力學(xué)性質(zhì)。張季如等[6-7]開展室內(nèi)壓縮試驗分析了鈣質(zhì)砂的壓縮性能和破壞特征。WANG等[8]基于滲透試驗結(jié)果，對顆粒級配、孔隙比和滲透性之間的關(guān)系進行了評估。丁選明等[9]結(jié)合振動臺模型試驗,對孔隙水壓力、加速度、位移和動應(yīng)變等動力響應(yīng)進行了測試和分析。基于室內(nèi)試驗，相關(guān)學(xué)者在鈣質(zhì)砂本構(gòu)模型的發(fā)展方面也取得了進展[10-11]，但由于目前在試驗過程中缺少實時監(jiān)測顆粒級配變化的方法，大都只關(guān)注了試樣被加載至破壞或臨界狀態(tài)時的顆粒破碎量，而較少探討具體應(yīng)力路徑下顆粒破碎的中間發(fā)展過程，因此目前所建立的顆粒破碎模型對于具體應(yīng)力路徑的適用性仍有待驗證[11]。此外，當(dāng)前對鈣質(zhì)砂力學(xué)特性的研究主要集中于宏觀層面的室內(nèi)試驗與原位測試，在顆粒尺度對其運動、接觸、破碎等微觀構(gòu)造演變行為方面的研究仍較為缺乏[12]，對顆粒級配分布、孔隙分布、顆粒形狀、孔隙形狀等影響滲蝕發(fā)生的內(nèi)部因素研究較少[13]。

近年來，高分辨率成像設(shè)備和圖像處理技術(shù)的迅速發(fā)展為巖土材料內(nèi)部結(jié)構(gòu)的高精度無損檢測和荷載作用下顆粒尺度力學(xué)行為評估提供了新的可能。當(dāng)前，高精度成像技術(shù)在顆粒材料的尺寸、方向和形態(tài)表征[14-15]以及接觸關(guān)系[16]、顆粒運動[17]、局部應(yīng)變[18]等顆粒尺度行為評估方面受到了廣泛應(yīng)用。相關(guān)學(xué)者[13,19-21]結(jié)合CT、SEM 掃描試驗結(jié)果和圖像分析算法對南海鈣質(zhì)砂顆粒三維圖像形貌和內(nèi)部孔隙結(jié)構(gòu)分布規(guī)律進行了研究，而在鈣質(zhì)砂顆粒形狀與尺寸參數(shù)相關(guān)性分析和顆粒間接觸關(guān)系的研究等問題仍有待進一步探索。

本文作者將側(cè)限壓縮試驗與CT斷層掃描技術(shù)相結(jié)合，獲取特定應(yīng)力狀態(tài)下樣本中鈣質(zhì)砂顆粒三維圖像。然后，運用自適應(yīng)圖像分割方法和主成分分析方法提取鈣質(zhì)砂顆粒尺寸、主軸方向等基本特征，進而量化鈣質(zhì)砂方向比、球度、凸度、棱角度等表征參數(shù)，并結(jié)合壓縮過程探究受力狀態(tài)下顆粒形態(tài)演化規(guī)律。同時利用膨脹算法獲取顆粒間接觸網(wǎng)絡(luò)，并實現(xiàn)受力狀態(tài)下接觸對數(shù)量、方向演化規(guī)律量化和可視化的統(tǒng)一，以期為后續(xù)鈣質(zhì)砂本構(gòu)模型開發(fā)與離散元數(shù)值模型驗證與標(biāo)定提供參考。

1 鈣質(zhì)砂CT掃描與圖像處理

本研究選取南海鈣質(zhì)砂顆粒為研究對象。試驗前，取適量砂樣經(jīng)清水洗掉鹽分后烘干備用。對原始試樣進行初篩后從(0.85，1.27]mm、(1.27，2.36]mm 和(2.36，3.35]mm 三個粒徑級別中取等量砂樣均勻混合，配置級配均勻的砂樣進行試驗。

1.1 “原位”側(cè)限壓縮試驗

本次試驗在浙江大學(xué)Nikon XTH 225/320 LC CT 掃描儀進行，并配合自主設(shè)計的微型側(cè)限壓縮裝置對鈣質(zhì)砂顆粒進行初始表征與原位壓縮試驗。側(cè)限壓縮裝置示意圖如圖1所示。裝置制作材料為聚甲醛，該材料密度較低，有效減少了因高吸收率造成的X 射線穿透過程中的能量耗散，同時和其他非金屬材料相比在抗壓強度、耐磨損性能上更具優(yōu)勢。加載過程通過旋轉(zhuǎn)頂部驅(qū)動螺桿，繼而推動活塞向下運動控制，螺桿與裝置外殼通過螺紋連接，螺距設(shè)置為0.5 mm，每旋轉(zhuǎn)一圈螺桿向下移動0.5 mm。螺桿底部設(shè)計為光滑的半圓球體，從而減少因摩擦而引起的下方活塞旋轉(zhuǎn)。裝置頂部、底部均設(shè)有力傳感器，量程為2 000 N。裝置內(nèi)樣品的直徑d=25 mm，高度H=25 mm，共有鈣質(zhì)砂顆粒2 000余個。

圖1 CT“原位”側(cè)限壓縮裝置Fig.1 In situ confining compression device

試驗過程中，壓縮裝置固定于X 射線光源與探測器間的旋轉(zhuǎn)平臺上。本試驗所選用的X 射線加速電壓為160 keⅤ，在樣品旋轉(zhuǎn)180°范圍內(nèi)均勻采集2 500 幀投影圖，其單幀曝光時間為500 ms。掃描結(jié)果分辨率設(shè)置為17.1 μm(單位體素長度)，組成顆粒的最小體素數(shù)量超過3 000，保證了砂粒真實形態(tài)表征的準(zhǔn)確性。

試驗步驟為：通過裝置頂部的驅(qū)動旋鈕對鈣質(zhì)砂顆粒壓縮至目標(biāo)荷載?？紤]到本試驗為位移控制加載，且在掃描過程中試樣反力可能隨時間發(fā)生松弛現(xiàn)象，故等待20 min 直至傳感器示數(shù)基本穩(wěn)定后啟動CT進行拍攝，記錄此目標(biāo)荷載下樣品三維影像。試驗共設(shè)置有100，1 500 和3 000 kPa 三級荷載。

1.2 CT圖像處理

為獲取試樣中單個顆粒的形態(tài)信息，本文采用自適應(yīng)分水嶺法[13]對掃描獲得的樣品三維CT圖像進行分割。X 射線穿過樣本后，受樣本物質(zhì)成分、密度差異影響其強度發(fā)生衰減變化，根據(jù)衰減系數(shù)生成CT掃描灰度矩陣(圖2(a))；然后基于雙重閾值算法確定合適的閾值對圖像進行二進制化(圖2(b))，從而實現(xiàn)實體材料與背景部分的分離。此時顆粒材料因相互接觸而無法被識別為獨立顆粒，通常采用基于分水嶺法的圖像處理技術(shù)進行分割。圖像分割過程通常存在過度分割(oversegmentation)和欠分割(under-segmentation)的問題，需對二進制圖像計算所得集水盆地深度進行調(diào)整[22-23]。因鈣質(zhì)砂形狀極不規(guī)則且尺寸分布較廣，設(shè)定全局控制參數(shù)調(diào)整集水盆地深度往往不適用。本文采用自適應(yīng)分水嶺分割法迭代進行圖像分割，在每一次迭代過程中，針對新分割出的區(qū)域計算相應(yīng)的局部調(diào)整參數(shù)，且在迭代過程中保留欠分割顆粒至下一次迭代繼續(xù)處理，從而實現(xiàn)圖像的精細(xì)分割。圖2(c)所示為最終分割圖像的二維切片，其中不同顏色區(qū)域即對應(yīng)不同顆粒?？梢钥吹皆摲椒ㄔ诜指钼}質(zhì)砂這一類不規(guī)則顆粒土樣的圖像中具有較好的效果。將掃描結(jié)果與石英砂掃描結(jié)果對比[24]，可以看出鈣質(zhì)砂顆粒形狀復(fù)雜，內(nèi)部保留有豐富的原生生物骨架細(xì)小孔隙。以上過程均在MATLAB 中進行，其他具體細(xì)節(jié)參考文獻(xiàn)[13]。

圖2 CT圖像分割(僅展示二維切片)Fig.2 Illustration of image segmentation(showing a 2D slice)

2 鈣質(zhì)砂三維形態(tài)表征

本文基于圖像分割結(jié)果對裝置內(nèi)全部鈣質(zhì)砂顆粒進行逐個分析，進一步提取鈣質(zhì)砂顆粒形心坐標(biāo)、體積和三維尺寸等基本信息，進而計算出鈣質(zhì)砂顆粒的主軸方向、平坦度、伸長度、球度和凸度等典型形態(tài)參數(shù)，實現(xiàn)鈣質(zhì)砂顆粒三維形態(tài)的定量表征。

2.1 顆粒形態(tài)參數(shù)確定

CT 掃描重構(gòu)可獲得代表鈣質(zhì)砂顆粒的一系列相同標(biāo)號三維頂點坐標(biāo)。由于鈣質(zhì)砂顆粒具有形狀復(fù)雜、不規(guī)則的特點，采用主成分分析(principal component analysis,PCA)方法對樣本中每個單獨顆粒的頂點坐標(biāo)進行正交變換，從而將原始坐標(biāo)轉(zhuǎn)化為互不相關(guān)的新變量長軸、中軸和短軸，進而實現(xiàn)不規(guī)則顆粒尺寸、方向等基本信息的標(biāo)準(zhǔn)化測量。根據(jù)PCA 計算結(jié)果，將其從小到大依次定義為顆粒短軸模量c、中軸模量b和長軸模量a。參考KONG等[14]對鈣質(zhì)砂顆粒軸長和篩分試驗對比結(jié)果，采用長軸模量a作為鈣質(zhì)砂顆粒等效粒徑D。

基于長軸模量a、中軸模量b和短軸模量c計算顆粒形狀比參數(shù)，從而對顆粒形狀進行度量。形狀比參數(shù)包括伸長度E和平坦度F，通常采用E=b/a和F=c/b分別進行計算[25]。圖3所示為2個代表性顆粒模型及其伸長度、平坦度計算結(jié)果?？梢钥闯錾扉L度E、平坦度F接近于1 時，顆粒各方向軸長模量近似相同。采用球度對顆粒球狀性進行定義，通常采用(其中，S0為顆粒球度，V為顆粒體積，SA為顆粒等效球表面積)進行計算[24]，但當(dāng)顆粒表面較為粗糙時使用該方法會使得計算結(jié)果產(chǎn)生較大誤差。因此，由于鈣質(zhì)砂具有多棱角、形狀復(fù)雜的特征，本文定義球度指標(biāo)S=Vfill/VS(其中，Vfill為填充內(nèi)孔隙后的顆粒體積，VS表示顆粒外接球體積)。圖4所示為2個代表性顆粒模型及其球度計算結(jié)果。

圖3 典型顆粒平坦度與伸長度Fig.3 Typical grains with different flatness and elongation values

圖4 典型顆粒球度Fig.4 Typical grains with different sphericity values

采用凸度對顆粒表面凹凸程度進行定義，采用公式C=V/VCH(其中，C為顆粒凸度，VCH為包裹住顆粒所需的最小凸包體積)計算。圖5所示為2個代表性顆粒模型及其凸度計算結(jié)果。從圖5可以看出：凸度越大，顆粒表面越平整，凹凸程度越小。

圖5 典型顆粒凸度Fig.5 Typical grains with different convexity values

為了對顆粒整體棱角特征進行評估，通常采用棱角度對顆粒表面形態(tài)進行量化。首先需對顆粒表面局部棱角特征進行定義，通過MATLAB 內(nèi)置函數(shù)vertexNormal 和vertexAttachments 計算某一頂點與其相鄰頂點矢量，獲取該頂點最大曲率k1和最小曲率k2。之后，本文通過顆粒表面頂點平均曲率來判斷該頂點是否為一個棱角，若該頂點處平均曲率km,j對應(yīng)的平均曲率半徑1/km,j小于顆粒內(nèi)切球半徑Rin，則判定頂點為一棱角。圖6所示為2個典型顆粒表面平均曲率計算情況。從圖6可以看出：棱角處對應(yīng)的平均曲率較大，平緩處平均曲率較小，表明平均曲率能較好表征鈣質(zhì)砂顆粒局部棱角特征。在此基礎(chǔ)上，定義棱角度A計算方法如下：

式中：kin為顆粒內(nèi)接球半徑的曲率；km,j為該顆粒第j個面元的平均曲率；Aj為第j個面元的面積。圖6中典型顆粒棱角度計算結(jié)果顯示，表面棱角豐富的鈣質(zhì)砂顆粒棱角度較大，而表面棱角較為平緩的顆粒棱角度較小，表明棱角度能夠較好地表征鈣質(zhì)砂顆粒整體棱角情況。

圖6 典型顆粒棱角度Fig.6 Typical grains with different angularity values

2.2 結(jié)果與討論

圖7所示為試驗過程中3次掃描獲得的鈣質(zhì)砂顆粒等效粒徑的累積分布曲線演化情況。從圖7可以看出：在加載過程中，鈣質(zhì)砂顆粒等效粒徑中位數(shù)D50僅發(fā)生微弱減小，即在壓力P為3 MPa 條件下，級配均勻的鈣質(zhì)砂樣本并未表現(xiàn)出明顯的顆粒破碎現(xiàn)象，與一般情況下大顆粒鈣質(zhì)砂在該壓力下出現(xiàn)較大的顆粒破碎不符，這可能是良好的級配顯著降低了鈣質(zhì)砂的破碎性所導(dǎo)致[26]。

圖7 壓縮過程顆粒等效粒徑演化Fig.7 Grain size evolution during compression process

FONSECA 等[16]建議采用顆粒主軸a方向作為具體顆粒的方向向量。將每個顆粒形心空間坐標(biāo)作為向量起點，主軸方向作為顆粒方向。為了探究壓縮過程中鈣質(zhì)砂顆粒方向演化關(guān)系，本文分別求出顆粒主軸方向與X軸、Z軸的夾角，并繪制出壓縮過程中顆粒方向的演化情況，如圖8所示。從圖8可以看出：在加載過程中，顆粒與X軸夾角逐漸減小，與Z軸夾角趨于垂直，即加載過程中樣品中顆粒方向趨于向垂直于受力方向的平面演化。

圖8 壓縮過程顆粒方向演化Fig.8 Grain orientation evolution in compression process

圖9所示為南海鈣質(zhì)砂與英國近海石英砂[24]三維形態(tài)參數(shù)的累積分布曲線，包括凸度和伸長度。從圖9可以看出：南海鈣質(zhì)砂的凸度小于英國近海石英砂的凸度，而伸長度更大，這定量地說明了相較于石英砂顆粒，南海鈣質(zhì)砂外形更加復(fù)雜，棱角更加明顯。此外，可以注意到相較于石英砂顆粒，南海鈣質(zhì)砂顆粒的所有三維形態(tài)參數(shù)的累積分布曲線的分布范圍更大。上述研究結(jié)果表明，南海鈣質(zhì)砂顆粒受風(fēng)化和剝蝕作用影響，形態(tài)更加復(fù)雜，棱角顯著，形態(tài)特征變異性較大；而石英砂顆粒在地質(zhì)搬運過程中趨于規(guī)則演化，形態(tài)特征變異性較小[27]。

圖9 鈣質(zhì)砂與石英砂顆粒形態(tài)參數(shù)分布對比[24]Fig.9 Cumulative distribution curves of various morphological parameters of calcareous sand and quartz sand[24]

圖10所示為南海鈣質(zhì)砂不同形態(tài)參數(shù)與顆粒尺寸的分布關(guān)系，包括球度和棱角度。從圖10 可以看出：隨著顆粒尺寸的減小，顆粒球度不斷增加，而棱角度出現(xiàn)明顯減小。該結(jié)果表明，細(xì)顆粒砂由于受地質(zhì)搬運過程歷時更長，顆粒形態(tài)趨于規(guī)則演化，棱角更加平緩。同時，受頂端壓縮影響，鈣質(zhì)砂顆粒球度分布范圍略有所提高，而棱角度分布范圍則出現(xiàn)整體下移，表明壓縮過程使鈣質(zhì)砂顆粒形態(tài)發(fā)生改變，且變化方向趨于規(guī)則。此外，分析結(jié)果表明伸長度等其他形態(tài)參數(shù)與顆粒尺寸間不存在明顯的相關(guān)關(guān)系，因此本文僅給出球度、棱角度與顆粒尺寸分布關(guān)系。

圖10 顆粒形態(tài)參數(shù)與粒徑關(guān)系Fig.10 Relationship between size and shape parameters

3 顆粒間接觸演變

顆粒間接觸關(guān)系是影響砂土力學(xué)性質(zhì)的重要因素之一。本研究采用基于MATLAB 膨脹算法的顆粒接觸關(guān)系捕捉方法，將目標(biāo)顆粒進行膨脹運算后與圖像分割后的周邊顆粒疊加，標(biāo)記出現(xiàn)接觸的部分從而確定配位數(shù)(coordination number)。此外，本方法僅對連接在同一分水線上的顆粒執(zhí)行接觸捕捉算法，有效避免了將非接觸顆粒識別為接觸對的錯誤發(fā)生，且與文獻(xiàn)[16]采用的算法相比極大地提高了運算效率。圖11所示為4 組典型顆粒接觸對識別結(jié)果，圖中藍(lán)色箭頭表示接觸對方向向量，其中向量起點表示顆粒間接觸點，向量方向指向兩顆粒接觸面法線方向。

圖11 典型鈣質(zhì)砂配位數(shù)和接觸關(guān)系Fig.11 Typical grains with different coordination number and contact network

圖12(a)所示為加載前顆粒配位數(shù)與粒徑分布情況。從圖12(a)可見：顆粒配位數(shù)與顆粒尺寸呈明顯的正相關(guān)關(guān)系，且顆粒最大可能配位數(shù)與等效粒徑呈明顯的線性關(guān)系。為研究壓縮過程中配位數(shù)與粒徑分布關(guān)系演化情況，分別計算特定粒徑范圍內(nèi)配位數(shù)均值，得到平均粒徑與平均配位數(shù)的關(guān)系(圖12(b))。由于粒徑在6～8 mm范圍內(nèi)的顆粒數(shù)量極少，為減小計算誤差僅對粒徑在1～6 mm范圍顆粒進行分析?？梢姡殡S壓縮過程土體密實程度增加，大顆粒鈣質(zhì)砂配位數(shù)出現(xiàn)明顯增加。

圖12 顆粒配位數(shù)與粒徑關(guān)系演化Fig.12 Evolution of relationship between grain coordination number and size

將樣品中所有顆粒的接觸向量識別結(jié)果繪制于同一坐標(biāo)系下，得到整個樣品接觸網(wǎng)絡(luò)如圖13(a)所示。本文將接觸面法向量方向定義為接觸方向，計算樣品中各接觸方向與X軸夾角，并繪制出壓縮過程接觸方向的演化情況，如圖13(b)所示。可見：樣品中顆粒間接觸方向明顯垂直于加載平面，且隨著加載過程進行，接觸方向繼續(xù)向垂直于加載平面方向演化，這與文獻(xiàn)[16，28]中認(rèn)為接觸方向趨向于平行于主應(yīng)力方向從而形成穩(wěn)定力鏈的結(jié)論一致。

圖13 樣品接觸網(wǎng)絡(luò)及壓縮過程接觸方向演化Fig.13 Contact network and evolution of contact vector in compression process

4 結(jié)論

1)在CT掃描室內(nèi)開展“原位”加載試驗并結(jié)合有效的圖像分析技術(shù)，能夠有效獲取鈣質(zhì)砂樣本顆粒尺寸、形狀、接觸等信息，實現(xiàn)對鈣質(zhì)砂等“問題性”土的微觀結(jié)構(gòu)表征。

2)圖像分析結(jié)果表明鈣質(zhì)砂顆粒形狀與顆粒尺寸呈現(xiàn)顯著的相關(guān)性：較大尺寸顆粒外形更為復(fù)雜且棱角更明顯，壓縮過程使顆粒形態(tài)趨于規(guī)則演化。

3)級配良好的南海鈣質(zhì)砂在3 MPa 側(cè)限壓縮壓力下未出現(xiàn)明顯顆粒破碎，可能由于壓縮過程中土樣微觀結(jié)構(gòu)重組使小粒徑顆粒填充于大顆粒間孔隙所致。

4)顆粒配位數(shù)與尺寸成明顯的正相關(guān)關(guān)系，顆粒最大可能配位數(shù)與等效粒徑呈明顯的線性關(guān)系；在土樣壓密過程中顆粒配位數(shù)明顯增加，但與等效粒徑相關(guān)性未發(fā)生明顯改變。