基于數字圖像灰度相關性的類巖石材料損傷分形特征研究

鄒 飛，李海波，周青春，朱小明，莫振澤，何恩光，趙 羽,

(1. 中國科學院武漢巖土力學研究所 巖土力學與工程國家重點實驗室，武漢 430071；2. 貴州省質安交通工程監控檢測中心有限責任公司，貴陽 550000；3. 北方重工盾構機公司，沈陽 110041；4. 浙江大學 現代制造工程研究所，杭州 310027)

1 引 言

巖石的損傷和破壞是一個漸進的過程，往往表現為微裂紋的萌生、擴展和貫通，通過微觀測量來定量刻畫巖石損傷演化所導致的力學性質的劣化是十分必要的。數字散斑相關方法作為一種全場非接觸式測量手段，其關鍵在于所選擇基準圖像與變形后圖像的相關性分析。該方法也被眾多研究者[1－7]用于巖石類材料損傷、破壞機制的研究中。馬少鵬等[1]采用數字圖像灰度相關性分布的標準差為統計指標表征巖石變形和破裂的發展，并通過試驗證明灰度性的分布演化與結構的變形局部化演化在空間和時間上都存在對應關系，該方法介于定性描述與定量分析之間。趙永紅等[4]采用數字圖像相關技術研究了含微裂紋巖石的變形，用相關技術處理灰度分布圖獲得位移場，對加載和卸載試件表面的灰度分布圖進行相關計算，獲得位移分布與裂紋分布的關系。朱珍德等[5]通過對紅砂巖試樣應變響應和表面裂紋圖像的同步觀測試驗以及試樣表面細觀裂紋萌生、擴展圖像信息的采集和量化，提出了用數字圖像相關技術處理不同含水狀態下紅砂巖試樣動態變形過程中的灰度分布圖，探討了紅砂巖動態損傷狀態與其宏觀力學之間的關聯。在對細觀結構劣化的描述中，很多學者采用簡單分形維和多重分形譜來進行刻畫[8－13]。彭瑞東[8]通過灰巖拉伸過程中表面細觀形貌的SEM數字圖像，采用分形維數描述損傷的發展過程。曹茂森[9,13]采用多重分形譜理論研究了混凝土結構表面裂紋的分形特征，通過分形特征量來衡量結構損傷的程度。

本文以石膏試件的單軸壓縮試驗為基礎，通過試件表面數字圖像灰度相關系數的變化來表征試件表面損傷狀態的演化過程。對不同時刻的灰度相關系數圖進行圖像處理，采用分形維數和多重分形譜定量描述材料表面細觀損傷狀態在時間和空間上的演化過程。

2 數字圖像原理

2.1 數字圖像灰度的相關性

圖像可以定義為一個二維函數f(x, y)，其中x和y是平面坐標，f是圖像在該點的振幅，也稱為該點處的亮度。將圖像轉化為數字形式的結果是獲得一個實數矩陣。對一幅圖像f(x, y)進行取樣后，得到m行n列的數組，這幅圖像大小即為m×n，數組的每一個元素稱為像素?；叶葎t是指黑白圖像中點的顏色深度，范圍為0～255，白色為255，黑色為0。

數字圖像的灰度分布 f(x, y)和g(x′, y′)，分別為初始點F(x, y)和移動點G(x′, y′)處的灰度值，在兩幅圖上的兩點間的水平、垂直位移分量u(x, y)和v(x, y)可以由下式[4]給出：

在荷載作用下，變形前后相關點的灰度強度具有如下映射關系：

在基準圖像中取一個包含點F(x, y)的具有m×n像素的子集f(m, n)。在變形圖像中同樣劃分m×n的小區域g(m, n)。這兩個圖像子集的灰度相關性按下式[6]進行計算：

2.2 灰度相關性和損傷的關系

由式（4）可知，相關系數 C是變形前后圖像小區域的灰度分布函數 f(m, n)和g(m, n)的函數。在荷載作用下，由于材料的不均一性、原始微裂紋的擴展，局部應力集中產生的顆粒剝離等導致試件拍攝面的圖像灰度在時間和空間上發生很大的變化，導致g改變增大，C值減小。文獻[1]通過相關性與變形的關系模擬試驗證明，大變形與低相關度是對應的，因此，灰度相關性在一定程度上反映了試件表面的變形程度。0

圖1 散斑圖及其相關系數等值線圖Fig.1 Speckle images and contours of correlation coefficient

按相關性公式（4），在基準圖像（見圖 1(a)）中間區域選擇一個20×20像素點的方格區域，在圖1(a)、圖 1(b)中進行匹配。在對兩幅散斑圖進行灰度相關性計算后，由圖 1(c)可知，在圖像中心區域有很好的單峰性，即在該區域有最佳匹配點，其灰度相關系數峰值為 1，說明采用的相關系數計算公式有很高的精度。圖1(b)是在圖1(a)基礎上添加隨機噪聲獲得的。采用該方法來模擬荷載作用下材料的劣化引起的灰度變化，由圖1(d)可知，峰值明顯降低，但單峰性依然明顯，其峰值為0.7?？梢娛剑?）對荷載作用下由于材料的劣化所引起灰度變化在空間上仍然有很好的匹配定位性。從而使得基于灰度相關系數峰值改變量和描述巖石空間損傷的損傷狀態之間建立起對應關系成為可能。

3 試驗結果及分析

3.1 試件、試驗裝置

為了獲得材料力學性質相對均一的試樣，采用高強度石膏在鋼模中澆注試件（保證試件表面較高的平整度），試件尺寸為5 cm×5 cm×10 cm的立方柱。為了使整個拍攝面的灰度反差較大有利于空間匹配，在試件表面用黑色玻璃微珠漆進行噴涂，形成隨機人工散斑場。采用RMT-150型巖石力學試驗系統進行單軸壓縮加載，加載速率為0.005 mm/s。在整個試件加載過程中用 CCD相機實時記錄試件表面散斑場的變化，采集分析軟件為DSCM觀測系統[1－3]，該系統的核心是圖像灰度相關性分析，對整個圖像先進行像素搜索，后進行亞像素搜索，獲得最佳匹配點，并獲得拍攝面的位移場和應變場。圖像采集系統的采集速率可控，最高為17幀/s，整個試驗過程共記錄約450幅圖像。

3.2 相關性分布演化與荷載對應關系

單軸壓縮應力-應變曲線如圖2所示，整個曲線階段性明顯，將曲線分為以下5個階段：

圖2 加載過程中應力-應變曲線及St演化曲線Fig.2 Evolution curve of Stand stress-strain curve during loading process

（1）OA：彈性階段，試件變形符合彈性變形，此段近似直線，符合線彈性本構關系；

（2）AB：硬化階段，A點為屈服極限，從該點開始試件產生塑性變形；

（3）BC：快速軟化階段，應變迅速增加，而應力不斷減??；

（4）CD：近似理想塑性階段，其特點是應力無明顯變化，應變不斷增加，該階段應力-應變呈理想塑性本構關系；

（5）DE：破壞階段，應力迅速跌落，試件壓碎破壞，僅保持一定的殘余強度。

在采集圖像中截取感興趣的代表性區域（ROI）進行計算，截取過程中盡量靠近試件的中間部位，又必須包含盡可能多的物面信息，減小加載過程中端部效應的影響。

采用文獻[1]的方法，通過不同時刻圖像與基準圖像對比計算，以獲得的灰度相關性系數的標?準差St為統計指標，其定義為

式中：Ckl為t時刻變形圖像中第k個小區域的相關系數值為t時刻變形圖像中m×n個小區域的相關系數值的均值。

從圖2中可知，在加載初期由于試件原始裂隙的壓實導致試件表面灰度調整，使St值有輕微擾動。在 OA、AB段標準差值變化不是很明顯；BC段中期由于微裂紋的出現，St值有明顯的陡增；CD段由于微裂紋的大量產生和擴展，St值迅速增加；DE段，試件變形急劇增大，表面宏觀裂紋產生和貫通，導致C值大幅度降低，St值以更大的速率增加。可見該統計指標能有效反映試件損傷變形的各個階段的特征。

3.3 損傷閾值的選取

假定OA階段為完全彈性階段，以AB段起點時刻A點（彈性應力上限）的試件表面灰度圖相關系數圖為閾值選取的基準圖，如圖3(a)所示（以O點時刻的灰度圖為基準圖像計算獲得）。為了獲得該時刻灰度相關系數的概率分布密度，采用核密度估計法。

一元連續的總體樣本在任意點x處的總體密度函數f(x)的核函數密度估計定義為

式中：K()稱為核函數；h稱為窗寬。

具體方法為：將通過DSCM系統計算獲得的灰度相關系數圖由二維矩陣轉換為一維矩陣。將該一維矩陣作為總體樣本，采用 Gaussian核函數，h=0.0011，求得核密度估計曲線如圖3(b)所示。核密度估計曲線與灰度相關系數的頻率直方圖符合的很好，為了進行對比，另外采用正態分布估計進行擬合，發現灰度相關系數的頻率直方圖與均值為0.9793，標準差為0.0101的正態分布曲線也非常的接近。

圖3 相關系數等值線圖及頻率直方圖、密度分布圖Fig.3 Contour, frequency histogram, density distribution of correlation coefficient

為了建立損傷和灰度相關系數改變量之間的關系，依據灰度相關系數的密度分布圖（見圖3(b)），并認為占總數0.05的最小值部分為初始損傷和采集過程中噪聲污染之和，該部分對應的最大 C值為0.9605。以0.9605為二值化轉換閾值對整個灰度相關系數圖進行圖像處理，C＜0.9605則認為該點發生大的變形，處于損傷或破壞狀態，在圖像矩陣中將其賦值為1；C＞0.9605則認為未損傷，將其賦值為 0。建立以時間為序列的損傷演化二值圖（二值圖中1代表白色像素，0代表黑色像素）。從AB、BC、CD、DE段中間時刻各選取一幅具有代表性的圖，如圖4所示。

圖4 損傷演化二值圖Fig.4 Binary images of damage evolution

這里將發生損傷或破壞的微元數Nd（白色像素總數）與微元總數N（像素總數）的比值定義為損傷變量 Dd，其范圍為 0～1。則試件表面的損傷演化如圖5所示，在加載初期，處于彈性階段，試件無明顯損傷，隨著應變的增加，損傷也隨之明顯增加，在破壞時有明顯的急速上升。

圖5 加載過程中應力-應變曲線及Dd演化曲線Fig.5 Evolution curve of Ddand stress-strain curve during loading process

4 損傷分形特征

分形幾何是研究自然界不規則現象及其內在規律的學科。巖石材料損傷和破壞過程是其內部缺陷不斷萌生發育、擴展、聚集和貫通的最終結果，這個從細觀損傷發展到宏觀破碎的過程具有分形性質[8]。但整個發展過程并不嚴格滿足自相似性的分形，往往只具有統計自相似性。對于上節損傷演化二值圖中損傷微元的分布，本文從統計自相似性出發，用分維數定量描述試件表面的損傷程度。

4.1 分形測定方法及特征

分維數的測定采用基于圖像處理的盒維數法[9]，具體方法如下：將相關系數分布圖進行閾值處理以后，獲得二值圖（本文中圖像像素為177×399），用邊長為l的方格（在圖像中為l個像素點），計算圖像被邊長為l的方格覆蓋占有的方格數N(l)，不斷改變方格邊長l的大小，并計算其相應的N(l)，可得到一組N(l)和l數據。根據盒維數的基本定義，損傷分維數D可由下式計算：

lgN(l)-lg(1/l)的曲線如圖 6所示，兩者呈線性關系，相關系數 R=0.9865，表明損傷微元分布具有自相似性，可以用分形進行表述。

為獲得圖2中各個階段的損傷與應變的對應關系，從塑性變形開始往后各階段的中間部分（避免轉折點處的突變）以5 s為時間間隔，獲得各個階段損傷演化二值圖。AB段選擇7幅圖，BC階段9幅圖，CD階段9幅圖，DE階段6幅圖。分別計算各階段在不同應變條件下的分形維數。

在對不同階段下相應損傷演化二值圖進行分形分析后，提取各個所選二值圖的分形維數。其維數的變化如圖7所示，隨著應變增加，各階段分形維數的范圍為 1.011～1.055、1.092～1.208、1.210～1.362、1.365～1.469，最小值為 1.011，最大值為1.469?？梢婋S著應變的增大，材料表面細觀結構的損傷程度增加，其對應的損傷演化二值圖上的損傷微元越來越多，分形維數也隨之增大；損傷微元由整體隨機分布，向局部集中發展。宏觀上表現為微裂紋的萌生、擴展，宏觀裂紋的出現。盒維數的大小反應了損傷的集中程度。不同階段分形維數和應變的一元線性回歸擬合系數分別為 0.0305、0.0567、0.0738、0.0930，由此可知，在不同階段損傷的發展是一個逐漸加快的過程，也從另外一個方面反應了試件加載作用下，力學性質加速劣化直到破壞的過程。

圖7 損傷分形維數與應變關系圖Fig.7 Relationships between fractal dimension and strain

4.2 多重分形測定方法

簡單的分形維數對試件表面損傷的演化只能作整體性、平均性的描述與表征。多重分形是定義在分形結構上的由無窮多個標度指數所組成的集合，是通過一個譜函數來描述分形結構上不同的局域條件或分形結構在演化過程中不同層次所導致的特殊結構行為與特征，是從系統的局部出發來研究其整體的特征，并借助統計物理學的方法來討論特征參量概率測度的分布規律[10]。

本文采用Ghhabra和Jensen[11]于1989年提出的計算多重分形譜的方法，首先用邊長為l的方格對圖像進行分網，損傷分布的總像素點為N，第i個方格所包括的損傷區域（白色像素點）為Ni，從而得到每個尺度為l的方格中分形圖像的概率測度：

Pi(l)的歸一化概率測度為

式中：q為權重因子，不同的q表示不同大小的概率測度 Pi(l)在配分函數pi(l)q中所具有的比重。

a(q)為描述概率測度局部奇異性強度的參數，表示多重分形集整體奇異性的平均值，其表達式為

f(a)為用奇異性標度指數標識的分形子集的Hausdorff維數，其表達式為

奇異譜 f[ a (q)]為 q∈(-∞,+∞)的單參數曲線,由于 f[ a (q)]與a(q)有惟一對應性，不同的a(q)及其所對應的 f[ a (q)]便構成多重分形的維數譜。

多重分形譜繪制的具體步驟為：首先程序讀入損傷二值圖，統計損傷微元總數N（圖像中白色像素點總數），再對二值圖像進行網格劃分（2n，n=1～8），程序對圖像中劃分的網格進行逐個識別統計，獲得每一方格損傷微元數目Ni，建立概率測度Pi=Ni/N，讀入q值，通過式（9）建立 ui(q, l)，計算獲得 ∑ui(q, l) lg[ui(q, l)]和 ∑ui(q, l) lg[pi(l)]的數值，并分別與lgl建立數據對，對上述兩組數據對采用最小二乘法進行一元線性回歸擬合，擬合出a(q)和 f[ a (q)]值，即可繪制多重分形譜圖?；谝陨纤悸罚肕atlab編制了基于數字圖像處理的多重分形譜計算程序。

4.3 損傷演化多重分形特征

圖8 ε=0.35%時多重分形譜Fig.8 Mutifractal spectrum with ε=0.35%

表1中僅列出應變為0.35%、0.70%、0.88%時的多重分形譜值。實際計算中 q取值范圍為-10～10，步長為 0.1。多重分形譜中大的a(q)反映的是小概率測度區域的性質，小的a(q)反映的是大概率測度區域的性質。圖9為圖4中4個不同應變條件下損傷演化二值圖的多重分形譜圖，其曲線由開始的對稱單峰形式，逐漸發展為左側相對集中，右側稀疏的鉤狀。 a(q)-f[ a(q)]曲線上的數據點在左側逐漸相對集中，表明大概率分布逐漸占主導地位，損傷演化從整體隨機均勻分布，到局部集中的過程。

圖9 多重分形譜Fig.9 Multifractal spectrum

多重分形譜比簡單分形包含了更豐富的結構信息，可以從分形譜中提取低維特征量來表征損傷的內在分形機制。當q=0、1、2時，D0為容量維，D1為信息維，D2為關聯維，其數值見表 2，三者的關系為D0>D1>D2。

表2 不同應變下D0、D1和D2的值Table 2 Values of D0, D1and D2at different strains

多重分形譜的寬度Δa=a (q)max-a(q)min的大小則反映了整個分形結構上概率測度分布的不均勻性的程度。Δa越大，表示損傷分布越不均勻，Δa越小，則表示損傷分布相對集中。ε=0.35%時，譜寬為0.381， a(q)-f(q)曲線呈對稱分布；ε=1.00%時，譜寬為 0.345，a(q)-f[ a (q)]曲線相對集中，其多重分形性相對減弱，損傷演化的奇異性也減?。ㄒ姳?）。

表3 不同應變下Δα的值Table 3 Values of Δα at different strains

5 結 論

（1）荷載作用下，試件表面的數字圖像灰度相關性和試件表面的損傷演化在時間和空間上存在對應關系，并可以用灰度相關系數C值的變化來對損傷狀態進行表述。

（2）通過簡單分形維數對由灰度相關系數建立的損傷演化二值圖進行刻畫，發現在各個不同階段，損傷分維數與應變呈線性關系。不同階段的增長速率不同，越接近破壞，增長速率越大。反映了材料在荷載作用下的加速劣化、破壞過程。

（3）通過多重分形譜反映了試件表面損傷演化在分形結構上不均勻分布的性質。其曲線形狀從開始的對稱光滑單峰形式，逐漸向左側集中發展，其譜寬逐漸減小，大概率分布逐漸占主導地位，反映了損傷演化從整體隨機分布，到局部集中的過程。

本文采用概率統計的方法，建立圖像灰度相關系數的改變和試件表面損傷狀態之間的關系，并對材料表面的細觀結構破壞進行了表述。文中僅僅對損傷狀態進行表述，更深層次的損傷機制及程度則需要更進一步的研究。

[1] 馬少鵬, 劉善軍, 趙永紅. 數字圖像灰度相關性用以描述巖石試件損傷演化的研究[J]. 巖石力學與工程學報,2006, 25(3): 590－595.MA Shao-peng, LIU Shan-jun, ZHAO Yong-hong. Gray correlation of digital images from loaded rock specimen surface to evaluate its damage evolution[J]. Chinese Journal of Rock Mechanics and Engineering, 2006,25(3): 590－595.

[2] 馬少鵬, 周輝. 巖石破壞過程中試件表面應變場演化特征研究[J]. 巖石力學與工程學報, 2008, 27(8): 1667－1673.MA Shao-peng, ZHOU Hui. Surface strain field evolution of rock specimen during failure process[J]. Chinese Journal of Rock Mechanics and Engineering, 2008,27(8): 1667－1673.

[3] MA S P, WANG L G, JIN G C. Damage evolution inspection of rock using digital speckle correlation method[J]. Key Engineering Materials, 2006, 326－328:1117－1120.

[4] 趙永紅, 梁海華, 熊春陽, 等. 用數字圖像相關技術進行巖石損傷的變形分析[J]. 巖石力學與工程學報, 2002,21(1): 73－76.ZHAO Yong-hong, LIANG Hai-hua, XIONG Chun-yang,et al. Deformation measurement of rock damage by digital image correlation method[J]. Chinese Journal of Rock Mechanics and Engineering, 2002, 21(1):73－76.

[5] 朱珍德, 張勇, 李術才, 等. 用數字圖像相關技術進行紅砂巖細觀裂紋損傷特性研究[J]. 巖石力學與工程學報, 2005, 24(7): 1123－1128.ZHU Zhen-de, ZHANG Yong, LI Shu-cai, et al. Studies of microcosmic crack damage properties of a red sandstone with digital image technique[J]. Chinese Journal of Rock Mechanics and Engineering, 2005,24(7): 1123－1128.

[6] CHU T C, RANSON W F, SUTTON M A. Applications of digital-image-correlation techniques to experimental mechanics[J]. Experimental Mechanics, 1985, 25(3):232－244.

[7] MA S P, JIN G C. Digital speckle correlation method improved by genetic algorithm[J]. Acta Mechanica Solida Sinica, 2003, 16(4): 366－373.

[8] 彭瑞東, 鞠楊, 謝和平. 灰巖拉伸過程中細觀結構演化的分形特征[J]. 巖土力學, 2007, 28(12): 2579－2582.PENG Rui-dong, JU Yang, XIE He-ping. Fractal characterization of meso-structural evolution during tension of limestone[J]. Rock and Soil Mechanics, 2007,28(12): 2579－2582.

[9] 曹茂森, 任青文, 翟愛良, 等. 混凝土結構損傷的分形特征實驗分析[J]. 巖土力學. 2005, 26(10): 1570－1574.CAO Mao-sen, REN Qing-wen, ZHAI Ai-liang, et al.Experimental study of fractal characterization in damages of concrete structures[J]. Rock and Soil Mechanics,2005, 26(10): 1570－1574.

[10] 李平, 胡可樂, 汪秉宏. 多重分形譜在材料分析中的應用研究[J]. 南京航空航天大學學報, 2004, 36(1): 77－81.LI Ping, HU Ke-le, WANG Bing-hong. Design and application about computing program of material multifractal spectrum[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2004, 36(1): 77－81.

[11] CHBABRA A, JENSEN R V. Direct determination on the f(α)singularity spectrum[J]. Physical Review Letters,1989, 62(12): 1327－1330.

[12] 王金安, 謝和平. 巖石斷裂面的各向異性分形和多重分形研究[J]. 巖土工程學報, 1998, 20(6): 16－21.WANG Jin-an, XIE He-ping. On anistropic fractal and multifractal properties of rock fracture surface[J].Chinese Journal of Geotechnical Engineering, 1998,20(6): 16－21.

[13] 曹茂森, 任青文. 鋼筋混凝土結構損傷檢測的分行特征因子[J]. 土木工程學報, 2005, 38(12): 59－64.CAO Mao-sen, REN Qing-wen. Damage detection of reinforced concrete structures based on fractal characteristic factor[J]. China Civil Engineering Journal, 2005, 38(12): 59－64.