基于SEM的土體微觀結構三維分析與分維計算方法

張豫川,高旭龍,劉東發,黃鴻偉

(1.蘭州大學 土木工程與力學學院,蘭州 730000; 2.蘭州大學 西部災害與環境力學教育部重點實驗室,蘭州 730000)

0 引 言

土體在長期的歷史沉積過程中,形成了多孔的微觀結構形態。微觀結構直接決定土體宏觀工程力學性質[1],也決定著土體外部環境改變時新結構的變化趨勢[2]。土體微觀結構研究一般包括顆粒形態與分布[3]、顆粒之間的接觸與排列、孔隙的形態與分布、孔隙間的孔喉連接[4]與迂曲程度[5-6]等內容。土體微觀結構的研究,對認識滲流[7]、變形、物質運移與水土保持等宏觀現象規律,揭示土體強度的維持與喪失等物理力學特性機理有極其重要的意義。

掃描電子顯微鏡(Scanning Electron Microscope,SEM)能無損、有效地獲取微觀信息,已廣泛應用于土體微觀結構的分析中[1,8-9]。一般將SEM圖像二值化,用黑白兩色分隔孔隙與顆粒后進行量化研究[10-11],但是閾值的選取對分析結果有直接影響,得到的視孔隙率、顆?；蛘呖紫兜钠矫婷娣e與周長等參數不具備直接分析的意義[11],將原本灰度值豐富的SEM圖像轉為只有0與255兩個灰度值的黑白圖也舍棄了大量可供處理與分析的信息。為了避免閾值選取,盡可能保留SEM圖像信息,王寶軍等[12]利用圖像處理技術結合地理信息處理系統(Geographic Information System,GIS),從黏性土微觀SEM中得到了三維視孔隙率等結構參數;徐日慶等[13]通過分析軟土SEM圖像,建立三維空間計算模型,開展了SEM分析對軟土孔隙率計算的影響研究;Li等[14]利用GIS技術分析固化污泥的SEM圖像,研究了固化污泥在生物降解過程中的微觀結構變化。

復雜的土體微觀結構難以用傳統歐式幾何準確量化,而分形理論在描述土體等結構復雜且具有自相似特征的多孔介質結構方面有著獨特優勢。分形理論由Mandelbort[15]提出后,最早由Tyler和Wheatcraft[3]引入土體學的研究中,分形維數可以定量表征土體結構,已在土體理化性質描述[16]、滲透、空間變異統計[17-18]等諸多方面有了成功的應用。土壤學研究中很早便建立了三維空間的分形維數模型及其分形維數的計算方法[3],但所需的計算參數很難由常規土工試驗數據得出,導致其在土體研究中的應用受阻。Tyler等[3]和楊培嶺等[19]假設同一土體顆粒具有相同密度,從土體顆粒質量分布計算分形維數,這一假設與實際情況不符[20],應用效果也差強人意[21]。土體孔隙分布的分形維數與顆粒不同,準確測量需要使用分子吸附法[22]。

本文對陜西定邊與甘肅蘭州的原狀土樣進行SEM試驗,基于SEM圖像三維化處理技術以及三維曲面表面積、體積計算方法,計算了土體視孔隙率、比面等結構參數,提出了三維化處理后微觀結構分形維數計算的實現方法。本文提出的方法可以僅由SEM圖像得出土體顆粒與孔隙的分形維數,無需額外的測定試驗,是土體微觀結構SEM圖像處理與分形維數計算的新路徑。

1 SEM圖像三維化處理方法

SEM圖像是灰度圖像,用純白(灰度值為255)、純黑(灰度值為0)以及二者之間一系列過渡色記錄觀測結果。掃描中發射出的電子束在樣品觀測面激發二次電子,被探測器收集,表現為或亮(白)或暗(黑)的差別,這個掃描成像過程與自然光照射在物體表面后反射類似,距離探測器較近的區域(凸起的顆粒)亮度較高,在圖像中更趨近于白色,距離較遠處(下凹的孔隙)更趨近于黑色。

二值化處理是預先選定灰度閾值,低于此閾值的像素點灰度值歸為0,記錄為孔隙范圍,高于閾值者歸為255,記錄為顆粒范圍,由此將灰度圖像轉為黑(孔隙)-白(顆粒)二值圖。

SEM圖像三維化處理無需使用閾值劃分,保留了灰度信息,其原理是在平面圖像x-y兩個維度基礎上,增添由灰度值轉換而成表征觀測面起伏的z維度坐標,形成三維圖像。如圖1所示,取3 pixel×3 pixel的平面像素空間為例,在圖像平面內建立x-y坐標系以描述各像素點位置,將每個像素點的灰度值轉為z坐標值,即

圖1 二值化與三維化處理SEM圖像原理對比Fig.1 Comparison between binary and three-dimensional SEM image processing principles

z(x,y)=g(x,y)。

(1)

式中:z(x,y)為三維化后(x,y)像素點的z坐標;g(x,y)為SEM圖像中(x,y)像素點的灰度值(無量綱)。

2 分形維數計算理論基礎及實現原理

2.1 分形維數計算理論基礎

分形維數是基于分形理論研究土體微觀結構的重要參數,具有明確的統計意義,土體的微觀結構可以用孔隙、顆粒等方面的分形描述。下面給出三維化處理后土體微觀分形模型的定義與計算方法。

歐式維數為E(E∈[1,3] )的幾何體如果是分形體,那么它的量度可用式(2)表示[23],即

G(ε)=G0εE-D。

(2)

式中:ε為尺度;G為ε尺度下的量度;G0為分形系數(常數);D為分形維數。

式(2)適用于分形曲線、面積與體積的量度。當E=1,G與ε對應于長度,此時D∈[1,2);當E=2,G與ε對應于面積,D∈[2,3);當E=3,G與ε對應于體積,D∈[3,4)。以E=1為例,分形維數D在幾何上描述了曲線的不規則程度,D越小,曲線越光滑,一維歐氏空間中的任意光滑(規則)曲線的分形維數D=1;D越大,曲線越曲折,越不規則,D→2的極限情況下,曲線極不規則,趨近于填充整個二維空間。

對于E=2的情況,分形體的量度G變成了面積S,此時分形維數DS表征了表面積的不規則程度,即“表面起伏分形維數”

S(ε)=S0ε2-DS。

(3)

兩邊取對數即可得到分形維數DS的計算式為

(4)

可改寫為

(5)

式中S0為對應于表面積S的分形系數。

假設土體微觀顆粒(孔隙)表面具有分形性質,SEM圖像三維化處理后,利用三角形不規則網格原理求取表面積S(ε),計算尺度ε取SEM圖像最小像素尺寸的整數倍,則相鄰像素點水平間距j=ε,j或者ε取值越大,表面積S(ε)計算精度越低,反之精度越高,可得到不同ε對應的S(ε)計算值。由式(5),作S(ε)/ε2-ε雙對數圖,若可線性擬合,則證實土體微觀顆粒(孔隙)表面積具有分形性質,分形維數DS即為擬合直線的負斜率。

上述DS的定義與用正方形格子覆蓋海岸線計算分形維數的方法[23-24]類似,故也被稱為“盒維數法”或“覆蓋法”,得到的分形維數DS也被稱為土體結構的表面起伏分形維數。

2.2 分形維數計算的實現原理

2.2.1 土體孔隙與顆粒的識別

三維化處理后灰度的明暗轉為z方向坐標值的高低,形成了試樣觀測面的直接映像,此三維化的起伏表面即是顆粒與孔隙的分界面。如圖2所示,在三維化圖像的最低點作一與x-y面平行的平面,此最低面與分界面之間的空間認定為顆粒,類似地,最高面與分界面之間的空間認定為孔隙。

圖2 土體孔隙與顆粒的識別示意Fig.2 Identification of soil pores and particles

2.2.2 顆粒(孔隙)的表面積與體積計算

孔隙與顆粒的表面積、體積等幾何參數是計算分形維數的基礎,可利用三角形不規則網格原理[25-26]計算。如圖3所示,三維化起伏面上取4點A、B、C和D,將坐標點相互連接可得到兩個空間三角形,它們的面積之和可作為這4點所圍區域的表面積近似值,4點所圍區域與任意平面之間的空間體積,可用空間三角形與該平面所夾體積近似,縮小4點之間的間距,即可逐漸逼近精確值。

圖3 三角形不規則網格計算原理示意Fig.3 Calculation principle of triangular irregular mesh

用海倫公式計算圖3中空間三角形面積之和,即

(6)

其中:

式中:Si為第i個不規則網格空間面積之和;|AB|代表圖3中A、B兩點距離,以此類推。各點之間距離用空間坐標計算。

令Δx=xi+j-xi,Δy=yi+j-yi,假設SEM圖像由m×m個像素點構成,則三維起伏面的表面積S可由式(7)得到

(7)

將圖3所示2個五面體沿頂面的最低點分割為三棱柱與三棱錐后計算空間三角形與x-y平面所夾空間體積,即

(xi+j-xi)(yi+j-yi) 。

(8)

由m×m個像素點構成的SEM圖像,三維起伏面與x-y平面所夾空間體積V為

(9)

分形維數DS定義的一個特點是分形維數由測量尺度ε與計算值S(ε)、V(ε)之間的關系得到。SEM圖像經過三維化處理后,聯合三角形不規則網格原理計算不同精度的S(ε)、V(ε)計算值,可以實現DS的計算。

3 三維化處理與分形維數計算實例

3.1 試驗材料與儀器

試驗用原狀土樣分別取自陜西省榆林市定邊縣與甘肅省蘭州市城關區青白石。定邊土樣編號為DB,蘭州土樣編號為LZ?；就亮W指標與顆粒的粒度分布如表1所示。

表1 試驗用土樣基本土力學指標與粒度分布Table 1 Basic mechanical indexes and particle size distribution of soil samples for the test

已有的試驗研究結果[7]表明,兩地土體都為濕陷性黃土,但滲透、應力響應等宏觀性質有較大差異,這正是二者微觀結構不同的體現。故而選擇定邊與蘭州兩地土樣進行微觀結構觀測,對比二者的微觀結構量化參數,并在此過程中說明本文提出的三維分析與計算分形維數等結構參數的方法。

原狀樣自然風干,避免迅速失水引起的土體結構變化[27],之后切削成直徑7 mm左右的柱體,使試樣可以放置在基座上。在柱體試樣中間用小刀輕劃一圈劃痕,沿劃痕掰斷試樣,得到保持結構原始狀態的新鮮斷面。噴金處理后移入JSM-6510型掃描電子顯微鏡,觀察試樣微觀結構特征。

3.2 SEM圖像三維化處理

DB、LZ這2個樣品放大到500倍(10 pixel/μm)SEM圖像如圖4所示。圖像尺寸均為2 560 pixel×1 920 pixel。

圖4 樣品的SEM圖像Fig.4 SEM images of samples

ArcGIS技術可將灰度值轉換為數字高程模型(Digital Elevation Model,DEM)中的高程,即z維度坐標值,實現第1節所述SEM圖像三維化處理。具體操作步驟與DEM的詳細內容,可參見參考文獻[28] ,三維化處理后的圖像如圖5所示,三維圖像根據z坐標值映射了色彩,用白色小箭頭標示觀察方向。

圖5 樣品SEM圖像三維化處理Fig.5 Three-dimensional processing of SEM images

3.3 分形維數計算結果

利用三角形不規則網格原理求取顆粒(孔隙)表面積S(ε)與體積V(ε),設定尺度ε=1,2,5,10,20,40,80,相應地改變計算網格大小可得到不同計算尺度下的S(ε)與V(ε)計算值。尺度為1時的計算值即為此SEM圖像所能計算的最精確值,放大尺度,計算量減小,但精度降低。不同計算尺度的SEM圖像示意如圖6所示。

圖6 部分計算尺度的SEM圖像Fig.6 SEM images of samples in partial calculation scale

圖7(a)是顆粒(孔隙)表面積計算值與量測尺度雙對數關系曲線,表面積計算值隨尺度不斷增大而大幅降低,DB樣在尺度取值范圍內可線性擬合,LZ樣在尺度取值1～20范圍內線性擬合效果較好,尺度>20后擬合效果逐步變差。圖7(b)是孔隙體積計算值與量測尺度關系曲線,體積計算值精度變化不如表面積計算值敏感,且存在一臨界值,臨界值內計算精度不會有明顯降低,臨界值與觀測土體微結構本身、觀測比例以及SEM圖像的像素總數等因素有關,本次試驗的臨界值為20。

圖7 計算值與量測尺度關系Fig.7 Relationship between calculated value and scale

土體孔隙的S(ε)/ε2-ε雙對數圖如圖8所示。

圖8 分形維數計算的雙對數圖Fig.8 Double logarithmic curves for fractal dimension calculation

圖8中DB樣與LZ樣均可線性擬合,且擬合效果較好,表明二者顆粒(孔隙)的表面起伏具有分形特征。分形維數計算結果匯總如表2所示。

表2 分形維數計算結果Table 2 Results of fractal dimension calculation

表面起伏分形維數DS描述孔隙表面起伏的不規則程度,DB樣的DS為2.729 9,大于LZ樣的值2.594 8,表明DB樣的孔隙表面起伏不規則程度大于LZ樣,DB樣孔隙表面較LZ更加“粗糙”。

3.4 分形維數的驗證與應用舉例

土水特征曲線是非飽和土本構關系的重要組成,以分形維數在土水特征曲線預測中的應用為例,驗證本文提出分形維數實現方法的計算結果。

用分形理論描述土體孔隙結構特性,可以得到非飽和土的土水特征曲線預測公式[21],即

(10)

式中:Se為有效飽和度, 由殘余含水量θr與飽和含水量θs標準化得到, 若認為θr為0, 則有效飽和度Se與飽和度Sr相同;ψ為基質吸力;ψa為進氣值, 指土體脫濕過程中氣體開始涌入孔隙時的基質吸力值;D為孔隙表面分形維數, 用本文中的DS表征。

試驗值由SOILMOISTURE公司生產的2100f小型張力計測得,分形模型預測結果與試驗值比較如圖9所示,其中進氣值ψa由最小二乘法對試驗點的最佳擬合得到。分形預測曲線與試驗值較為一致,印證了本文提出的三維化分形維數計算實現方法的可靠性。

圖9 分形模型預測結果與試驗值比較Fig.9 Comparison between fractal model calculation results and experimental results

4 視孔隙率與比面的計算

SEM圖像經三維化處理不僅可以得到分形維數,在實現土體結構表面的顆粒(孔隙)表面積S(ε)與體積V(ε)計算后,取ε=1時(最精確)的計算值還可以得到土體的視孔隙率與比面這2個觀測面的重要微觀結構參數。

4.1 視孔隙率與比面

孔隙率是土體最重要的微觀結構參數,SEM圖像只能得到觀測面的微觀結構信息,所得孔隙率可以稱之為“視孔隙率”,計算式為

(13)

式中:n為試樣觀測面的視孔隙率;Vv為孔隙(三維化圖像的最高面與分界面之間的空間)體積;Vs為顆粒(最低面與分界面之間)體積。

比面一般定義為土體固體(顆粒)骨架的總表面積與土體總體積之比,是土體等多孔介質傳熱、滲流、增脫濕等過程中十分重要的結構參數。三維化處理后可得到觀測面的顆粒(孔隙)表面積以及顆粒與孔隙的體積,由此得到比面,即

(14)

式中:Ω為比面,量綱為長度的倒數L-1;S為顆粒(孔隙)的表面積。土體等多孔介質的比面可以理解為單位體積內的孔隙表面積,比面越大,顆粒分布越分散,結構越復雜。

4.2 視孔隙率與比面計算結果

SEM目前常見的二值化處理也可得到視孔隙率,其值為黑色與白色區域的面積之比。三維化、二值化以及由常規土工試驗所得孔隙率對比見表3。

表3 孔隙率計算結果Table 3 Calculation results of apparent porosity

本文提出的三維化處理結果與試驗孔隙率更為接近,這是因為二值化處理通過設定閾值分割顆粒與孔隙,計算模型精度較低;三維化處理無需此種操作,計算精度主要取決于原始SEM圖像的像素點多少,以及設定的計算網格大小,因而相同觀測比例,三維化視孔隙率與試驗孔隙率更為接近。

需要說明的是,視孔隙率是觀測面孔隙狀態的反映,而試驗孔隙率表征的是試驗土樣體單元的孔隙狀態,二者概念不同,無法完全等同。從試驗結果可知,本文的三維化結果較二值化更加接近試驗孔隙率,表明三維化視孔隙率一定程度上能夠表征真實孔隙狀態,有望成為土體微觀結構定量分析的工具之一。

三維化處理后的比面計算值如表4所示。DB樣的比面Ω為2.809 6 μm-1,LZ樣為22.318 1 μm-1,DB樣Ω僅為LZ樣的12.56%。LZ樣比面大于DB樣,表明LZ樣的顆粒更小,顆粒分布更為分散,這一結論從二者的粒度分布試驗結果中可以得到驗證。

表4 比面計算結果Table 4 Calculation results of specific surface area

5 結 論

本文在已有的分形維數定義基礎上,提出了針對土體微觀SEM圖像的三維化處理以及分形維數計算實現方法,通過定邊與蘭州兩地原狀土樣的微觀結構觀測實例驗證了該方法的有效性。主要得出以下結論:

(1)三維化處理無需選擇閾值,保留了SEM圖像的灰度信息,可以實現微觀結構表面的三維重建和可視化,計算得到觀測范圍內土體顆粒(孔隙)的表面積與體積,繼而獲取視孔隙率、比面等土體微觀結構參數。

(2)分形維數DS的計算均可由測量尺度與計算值之間關系得到,本文提出的三維化處理與三角形網格原理聯合方法,可以計算土體顆粒(孔隙)不同測量尺度ε對應的S(ε)與V(ε)計算值,由此實現分形維數的計算。

(3)試驗所用定邊與蘭州土樣的微觀顆粒(孔隙)表面均具有分形特征。定邊土樣表面起伏分形維數DS為2.729 9,蘭州土樣DS為2.594 8,定邊土樣的孔隙表面起伏不規則程度大于蘭州土樣。

(4)計算得到的分形維數在土水特征曲線預測方面應用效果良好,印證了本文提出的三維化分形維數計算實現方法的有效性與可靠性。