南海鈣質砂微觀構造三維圖像分析與表征

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

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

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

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

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

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

1 鈣質砂CT掃描與圖像處理

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

1.1 “原位”側限壓縮試驗

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

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

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

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

1.2 CT圖像處理

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

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

2 鈣質砂三維形態表征

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

2.1 顆粒形態參數確定

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

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

圖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個代表性顆粒模型及其凸度計算結果。從圖5可以看出：凸度越大，顆粒表面越平整，凹凸程度越小。

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

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

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

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

2.2 結果與討論

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

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

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

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

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

圖9 鈣質砂與石英砂顆粒形態參數分布對比[24]Fig.9 Cumulative distribution curves of various morphological parameters of calcareous sand and quartz sand[24]

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

圖10 顆粒形態參數與粒徑關系Fig.10 Relationship between size and shape parameters

3 顆粒間接觸演變

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

圖11 典型鈣質砂配位數和接觸關系Fig.11 Typical grains with different coordination number and contact network

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

圖12 顆粒配位數與粒徑關系演化Fig.12 Evolution of relationship between grain coordination number and size

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

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

4 結論

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

2)圖像分析結果表明鈣質砂顆粒形狀與顆粒尺寸呈現顯著的相關性：較大尺寸顆粒外形更為復雜且棱角更明顯，壓縮過程使顆粒形態趨于規則演化。

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

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