基于數(shù)字散斑相關(guān)的土體剪切帶自動識別方法與應(yīng)用

徐曉華 ，李元海 ，朱鴻鵠，王迎超

(1. 中國礦業(yè)大學 深部巖土力學與地下工程國家重點實驗室，江蘇 徐州，221116；2. 中國礦業(yè)大學 力學與土木工程學院，江蘇 徐州，221116；3. 南京大學 地球科學與工程學院，江蘇 南京，210023)

剪切帶的萌生、發(fā)展和貫通是巖土體災(zāi)變過程中應(yīng)變局部化的宏觀表現(xiàn)，是多數(shù)巖土體變形破壞的征兆和前提，如邊坡土體的滑動破壞，室內(nèi)土工試驗土樣的剪切破壞等，自20世紀80年代以來，一直受到固體力學、材料科學、巖土力學等領(lǐng)域科技人員的密切關(guān)注[1-3]。

目前，在土體剪切帶的形成機理、識別和變形分析等方面，研究者展開了許多細致的研究，并取得了豐富的成果，但在模型試驗中，土體剪切帶的自動識別方法與應(yīng)用方面還需要進一步研究。在巖土模型試驗中，土體剪切帶的分析研究主要依靠內(nèi)部觀測和表面觀測2種方法。內(nèi)部觀測有X 射線、CT 掃描以及光纖傳感等技術(shù)。X 射線技術(shù)可以用于分析土樣在剪切過程中的微觀力學行為，但該方法所使用的測試設(shè)備昂貴，試驗操作復(fù)雜，剪切帶的形狀識別也比較粗糙[4]。CT 掃描技術(shù)可以通過掃描得到3D圖像來定量分析剪切帶的擴展演化過程，但設(shè)備昂貴[5-6]；利用埋入式的應(yīng)變傳感光纖，可以推算土體剪切變形量，但分析結(jié)果受到纖土耦合性的影響，需要借助于較為繁瑣的標定工作來提高精度[7-8]。

土體剪切帶的表面觀測方法主要有描畫網(wǎng)格法和數(shù)字照相法。ALSHIBLI 等[9-10]采用描畫網(wǎng)格法分析剪切帶形成的過程，通過在薄膜上描繪網(wǎng)格，比較攝像機采集的試驗圖像，測量剪切帶取向角和厚度，但該方法對操作要求高，且識別精度較低。與之相比，數(shù)字照相技術(shù)克服了接觸式變形量測的不足，操作簡單且成本低。李元海等[11]提出了一種基于DSCM 的土體剪切帶識別與量測方法，并應(yīng)用于一個大型砂土剪切試驗，其結(jié)果表明了該方法的可行性。金愛兵等[12]采用數(shù)字圖像相關(guān)(digital image correlation, DIC)技術(shù)開展了礦巖散體直剪試驗，研究了剪切作用下礦巖散體顆粒剪切面顆粒運動、局部應(yīng)變和剪切帶的演變規(guī)律。王學濱等[13]提出了改善剪切帶邊界附近位移測量精度的環(huán)繞測點子區(qū)分割一階DIC方法，通過對虛擬剪切帶位移進行測量以及與基于粒子群優(yōu)化和Newton-Raphson迭代的粗-細方法結(jié)果進行了對比，驗證了該方法的有效性。上述基于數(shù)字照相的土體剪切帶測量都是在得到高清圖像的前提下，利用相關(guān)圖像分析軟件進行識別，其過程需要輔助大量的人工處理，尚未實現(xiàn)土體剪切帶的自動、準確識別。

李元海等[11,14-15]自主研發(fā)的數(shù)字散斑相關(guān)分析軟件PhotoInfor可以有效分析模型試驗圖像的剪切帶的演變過程，但其分析過程中仍需人工輔助處理。為此，本文作者對其進行二次開發(fā)，提出并實現(xiàn)了土體剪切帶的自動識別，以減少人工識別的誤差，對土體剪切帶形狀準確識別和變形精細分析具有重要意義。

1 剪切帶數(shù)字散斑相關(guān)分析方法

1.1 數(shù)字散斑相關(guān)方法

數(shù)字散斑相關(guān)方法的基本原理是：對于利用數(shù)字圖像設(shè)備連續(xù)采集的含位移的數(shù)字散斑圖像，在首張圖像中設(shè)置量測像素點(這里稱之為“測點”)，以測點為中心，構(gòu)建一個固定大小的像素區(qū)域(又稱“子區(qū)”[16])，通過圖像相關(guān)性搜索的方法計算與相鄰圖像上指定搜索區(qū)域內(nèi)所有測點子區(qū)像素點灰度的相關(guān)系數(shù)，并將最大相關(guān)性系數(shù)對應(yīng)的測點作為搜索目標點，以此實現(xiàn)序列圖像上設(shè)定測點位置的連續(xù)追蹤，然后，根據(jù)圖像間對應(yīng)測點的坐標變化進行微小位移計算，通過連續(xù)累加獲得測點的最終位移，從而實現(xiàn)觀測目標數(shù)字散斑圖像上測點的位移測量。

本文采用的PhotoInfor軟件是一款自行研制的高精度巖土變形量測軟件系統(tǒng)[17-18]，相對于國內(nèi)外同類軟件，PhotoInfor 針對巖土材料非連續(xù)與局部化兩大變形特點進行了專門處理。軟件采用穩(wěn)定的相關(guān)性計算公式進行圖像分析和位移計算，應(yīng)用有限元四邊形等參單元變換的方法實現(xiàn)圖像校準和應(yīng)變計算。同時，為適應(yīng)砂土模型復(fù)雜的變形特點，采用平移和旋轉(zhuǎn)搜索方法進行位移追蹤，可滿足土體剪切帶識別、分析的精度要求。

1.2 土體剪切帶分析流程

利用PhotoInfor軟件對試驗中采集的圖像進行剪切帶分析，主要包括剪切帶形狀和邊界點的自動識別、變形分析和厚度計算，具體分析流程見圖1。其中，圖像采集分辨率越高、測線布置越密集都可以有效提高剪切帶的識別精度，但會大大增加計算量，因此，在實際操作過程中，應(yīng)根據(jù)實際需求進行調(diào)整。由于本文研究的主要內(nèi)容在土體剪切帶形狀、邊界點自動識別以及剪切帶的變形分析和厚度計算，有關(guān)流程的詳細介紹可參見文獻[14]。

圖1 剪切帶分析流程Fig. 1 Shear band analysis process

2 砂土剪切帶自動識別與變形計算方法

2.1 剪切帶形狀識別

由于顆粒大小和顏色各不相同，砂土圖像呈現(xiàn)為明顯的散斑狀，無需手動繪制散斑場。因此，本文基于砂土直剪試驗采集的連續(xù)剪切變形圖像，采用PhotoInfor 進行剪切帶形狀分析，其過程如下。

1) 剪切帶范圍圈定?？紤]到多測點圖像分析計算量大且耗時長，并且在剪切帶范圍不明的條件下，布置較少的測點不利于剪切帶的識別，布置較多的測點又會增加圖像分析的時間和計算量，可以通過稀疏網(wǎng)格測點進行初步粗略分析，找到剪切帶的分布區(qū)域，縮小剪切帶測線的布置范圍，如圖2所示。

2) 剪切帶正式分析。在已知剪切帶分布范圍的基礎(chǔ)上，布置垂直穿越剪切帶的多條測線，每條測線由連續(xù)排列的測點組成，各測點之間相距1像素。通過測點的連續(xù)追蹤，計算測線上各測點的位移，獲得多組測線的“彎折”形態(tài)來實現(xiàn)剪切帶形狀的自動識別。方法原理如圖3 所示，其中，圖3(b)中的褐色填充區(qū)域為分析得到的剪切帶區(qū)域，將相鄰4個邊界點首尾相連構(gòu)成的封閉四邊形稱為測點單元。

剪切帶分析結(jié)束后，需要找到剪切帶邊界點來構(gòu)造測點單元以便分析變形和計算厚度，目前在邊界點識別方面主要依賴人眼識別和人工選取，效率低且受主觀性影響大，并且由于邊界點附近測點較為密集，難以依靠人眼選出邊界點。因此，本文對PhotoInfor軟件中的剪切帶識別功能進行完善，以實現(xiàn)土體剪切帶邊界點的自動識別功能。

2.2 剪切帶邊界點自動識別方法

針對上述砂土剪切帶分析中測線會在邊界點點處產(chǎn)生明顯“彎折”的特征，基于測點的坐標，提出一種基于最小二乘法擬合與曲率最大值判斷的剪切帶邊界點自動識別方法，具體算法如下。

1) 亞像素線性插值。依據(jù)單條測線的坐標可繪制其對應(yīng)的測點分布圖，如圖4所示。從圖4可以看出，剪切帶兩側(cè)的測點分布較為密集，而剪切帶內(nèi)部測點較為稀疏，其主要原因是剪切帶內(nèi)部位移變化較大，使得測點之間的間距也變大，測點數(shù)目變少。

圖4 單條測線變形后測點分布Fig. 4 Distribution of measuring points after deformation of a single measuring line

為了保證擬合曲線盡可能貼合測點分布，防止曲線有過大起伏，需要在距離過大的測點間進行插值?？紤]到本文測線繪制時測點的布置間隔設(shè)置為1像素，可以近似認為任意兩測點間的亞像素點變形是相同的。因此，本文采取線性的方式對測點間距高于一定閾值的部分進行插值處理(圖4)。

2) 擬合曲線。針對離散測點，采用多項式擬合方程和最小二乘法擬合數(shù)據(jù)。在實際擬合中發(fā)現(xiàn)，當多項式最高項次數(shù)過小時，擬合效果較差，無法達到預(yù)期要求；而當多項式次數(shù)過高時，又會造成過擬合，曲線彎折點較多。因此，最終將多項式最高項次數(shù)確定為8次。擬合線方程為

式中：i為冪指數(shù)，最小取值為0，最大取值為8；ai為第i項系數(shù)。測點擬合曲線見圖5。

圖5 測點擬合曲線Fig. 5 Schematic diagram of measuring point fitting curve

3) 邊界點識別。剪切帶邊界點通常為土體由應(yīng)變緩慢增加到急劇增大的轉(zhuǎn)折點，在擬合曲線上則表現(xiàn)為彎曲程度最大的2 個點。依據(jù)這種思想，本文采用曲率作為彎曲程度的衡量指標，通過式(2)計算曲線上各點的曲率K(x)，找到其中最大的2個極值點即為最大彎曲點。

式中：F′(x)和F″(x)分別為擬合曲線一階與二階導(dǎo)數(shù)，

圖6所示為各點曲率的分布圖。從圖6可以看出，有2個明顯較大的極值點，對應(yīng)曲線的最大彎折點；此外，還有其他較小的極值點，這主要是因為剪切帶外的部分測點分布不均勻，擬合時形成了較小的彎折，對識別效果影響較小。

圖6 擬合曲線的曲率分布圖Fig. 6 Curvature variation of fitting curve

實際上，曲率最大點并不是剪切帶測線上的邊界點，而是離邊界點最近的擬合曲線上的點。因此，需要進一步計算各測點到曲率最大點的距離，距離最小的測點為實際邊界點。圖7所示為某條測線穿越的剪切帶邊界點最終識別結(jié)果。

圖7 剪切帶邊界點識別結(jié)果Fig. 7 Identification result of shear band boundary point

2.3 不同傾角下剪切帶邊界點的識別

在實際剪切試驗中，剪切帶可能會形成不同的傾角，其測點分布也會呈現(xiàn)不同的形狀，影響曲線的擬合效果。為此，本文通過測點旋轉(zhuǎn)的方法模擬不同傾角下剪切帶測點分布，開展曲線擬合研究，并引入擬合優(yōu)度R2對擬合效果進行定量評價。R2的取值范圍為[-∞， 1)，其越接近1，則擬合效果越好。式(3)和式(4)分別為測點逆時針旋轉(zhuǎn)的坐標變換公式和R2計算公式。

式中：(X′i，Y′i)為旋轉(zhuǎn)后測點坐標；(Xi，Yi)為初始測點坐標；(XS，YS)為基準點坐標；θ為逆時針旋轉(zhuǎn)的角度；yi為擬合值；Yi測點坐標，-Y為測點平均坐標，n為測點個數(shù)。

圖8所示為通過旋轉(zhuǎn)前后不同傾角下的測點的擬合效果。從圖8可以看出，對于旋轉(zhuǎn)前測點分布圖，剪切帶傾角較大，擬合曲線的R2偏差值較大，測點擬合效果極差，無法繼續(xù)搜索邊界點；對于旋轉(zhuǎn)后測點分布圖，剪切帶傾角較小，擬合曲線的R2接近于1，曲線擬合效果較好。因此，亟需研究一種適用于不同傾角的邊界點識別算法。

圖8 旋轉(zhuǎn)前后測點的擬合效果Fig. 8 Fitting effect of measuring points before and after rotation

針對上述問題，本文提出一種最優(yōu)擬合線搜索的邊界點識別算法，即通過測點旋轉(zhuǎn)與擬合效果評定來確定擬合線，繼而求得邊界點坐標。具體算法如下。

1) 旋轉(zhuǎn)范圍的確定。為了保證搜索范圍能夠包含所有的傾角，本文選擇-90°～90°作為旋轉(zhuǎn)范圍，考慮到旋轉(zhuǎn)角度的間隔相差較小時對擬合曲線影響不大，且會增加計算時間，因此，將旋轉(zhuǎn)間隔確定為0.5°。

2) 測點旋轉(zhuǎn)與擬合。以第一個測點為基準點，依次在-90°～90°范圍內(nèi)每隔0.5°對測點分布圖執(zhí)行旋轉(zhuǎn)操作，然后按照上述所示的方法對旋轉(zhuǎn)后的測點進行插值與曲線擬合并計算相應(yīng)的R2。

3) 邊界點識別。找出R2最大時的旋轉(zhuǎn)角度β，將該角度下的擬合曲線確定為最佳擬合線，并對此曲線執(zhí)行邊界點識別。圖9所示為R2最大時的邊界識別圖，可見：邊界識別曲線擬合效果與識別結(jié)果都符合預(yù)期。最終將識別到的邊界點繞基準點順時針旋轉(zhuǎn)β即可得到實際邊界點。

圖9 最優(yōu)擬合線下的邊界點識別結(jié)果Fig. 9 Recognition results of boundary points under the optimal fitting line

2.4 剪切帶厚度計算方法

目前剪切帶的厚度測量主要有理論推導(dǎo)和試驗測量2種方法。本文在自動求得邊界點坐標的基礎(chǔ)上，提出一種求解剪切帶厚度的簡單方法。首先，將相鄰4 個邊界點相連構(gòu)成1 個四邊形單元，并將量測線上邊界點(下邊界點)相連的直線作為該等參單元的上邊界線(下邊界線)；然后，根據(jù)邊界點的坐標，計算得到每個邊界點到非所在剪切帶邊界線的距離，最終求取這4個距離平均值作為剪切帶的厚度。圖10 所示為剪切帶厚度計算原理的示意圖。圖10 中，A、B、C、D為相鄰兩條測線上的4個邊界點，分別求出A、B點到線段CD的距離和C、D點到線段AB的距離，以D點到線段AB距離D4的求解為例，假設(shè)A、B、D點的坐標分別為(xA，yA)、(xB，yB)、(xD，yD)，通過式(5)求得線段AB、BD和AD的長度LAB、LBD和LAD，則D4可由式(6)計算得到。

圖10 剪切帶厚度計算原理Fig. 10 Principle of shear band thickness calculation

同理，可求得A點和B點到線段CD的距離D1和D2以及C點到線段AB的距離D3，最終剪切帶的厚度可表示為

2.5 剪切帶變形計算方法

測點的位移可以通過不同圖像上對應(yīng)點坐標的變化直接求得。借用有限元中常用的四邊形等參單元的概念與方法，可以實現(xiàn)測點應(yīng)變的計算[17]。在自動識別得到剪切帶邊界點坐標后，以4個邊界點為角點的矩形可看作1個四邊形單元，在4個角點的坐標和位移已知的情況下，通過圖11所示的坐標變換可進行中心點的位移計算，計算公式如式(8)所示。

圖11 基于四邊形等參單元的坐標變換Fig. 11 Coordinate transformation based on quadrilateral isoparametric element

式中，(si，ti)為局部坐標系下第i個頂點的坐標；i=1，2，3，4；(s，t)為點(x，y)對應(yīng)的局部坐標；(ui，vi)為第i個角點在x和y方向上的位移分量。

由式(9)可求出四邊形中任意一點(x，y)的正應(yīng)變εx、εy以及剪應(yīng)變γxy和體積應(yīng)變εv。

式中：(X，Y)為變換后的坐標。

3 應(yīng)用試驗

3.1 試驗概況

采用自行設(shè)計的剪切試驗裝置，試驗?zāi)Ｐ烷L×寬×高分別為300 mm×200 mm×40 mm，采用日本標準豊浦砂，應(yīng)用空中下落法制作而成，砂漏出口在剪切盒上方約1 m處，以剪切盒豎向中心線為準，砂漏左右移動范圍各為30 cm，移動速度為30 cm/s，移動6～7次將砂土填滿剪切盒。具體試驗過程見文獻[11]。制作的模型相對密度可以達到85%以上，平均密度為1.60 g/cm3。水平剪切前后土體變化如圖12 所示，圖中，F(xiàn)V為垂直荷載，F(xiàn)T為水平荷載。經(jīng)過剪切后，土體上下部分發(fā)生錯動。

圖12 砂土剪切試驗示意圖Fig. 12 Schematic diagrams of sand shear test

3.2 剪切帶識別

圖13(a)和(b)分別為試驗圖像在剪切帶分析前后的測線布置情況。由圖13(a)和(b)可以看出：砂土中間部分出現(xiàn)明顯的應(yīng)變局部化現(xiàn)象，即剪切帶形成與發(fā)展區(qū)域。由自動分析得到的邊界點位置可見剪切帶在測線3 出現(xiàn)了明顯下降。圖13(c)所示為剪切帶邊界點自動識別效果和繪制的測點單元，可見，在剪切帶分析完成后，通過二次開發(fā)的自動識別功能，可自動計算并將邊界點位置與序號標在原圖上，將這些邊界點按順序依次相連，得到的帶狀區(qū)域即為剪切帶的形狀。

圖13 剪切帶分析與四邊形等參單元標識Fig. 13 Shear band analysis and drawing of quadrilateral isoparametric element

總體來看，識別程序可以較好地識別出剪切帶的邊界點和形狀。將相鄰的4 個邊界點依次相連，可以得到圖13(c)所示的5 個測點單元，分析這些測點單元可得到剪切帶不同位置的位移變化和剪切帶厚度。

3.3 邊界點識別結(jié)果對比

為了驗證邊界點自動識別的準確性，分別對邊界點采用人工選取和自動識別2種方法進行對比驗證，對比結(jié)果如表1所示，其中，ID號與測點一一對應(yīng)，人工選取結(jié)果與自動檢測結(jié)果越接近表示識別效果越好。

表1 邊界點識別結(jié)果(測點ID號)對比Table 1 Comparison between results of two boundary point identification methods

由表1 可得，人工識別與自動識別的測點ID號基本一致，僅在上邊界點2和上邊界點5處略有差異，其主要原因是這兩處測點較為密集，判斷難度大。自動識別能夠提高識別效率，且能夠有效避免人工識別的誤差。

3.4 剪切帶變形分析

識別出剪切帶邊界點后，每相鄰4個邊界點可構(gòu)成1 個四邊形單元。除了具有識別邊界點功能外，PhotoInfor 還有變形自動計算功能，通過分析剪切帶得到了正應(yīng)變、最大剪應(yīng)變γmax和體積應(yīng)變隨剪切位移的變化規(guī)律，如圖14所示。

圖14 剪切帶不同位置處應(yīng)變-剪切位移關(guān)系曲線Fig. 14 Strain and shear displacement relation curve at different positions of shear band

從圖14(a)和(b)可以看出：在砂土剪切過程中，x和y方向上正應(yīng)變的變化過程相似，當剪切位移s在臨界點前(s=55～60 mm)，正應(yīng)變保持為零或緩慢增加，此階段剪切帶變形相對緩慢；當s超過臨界值時，正應(yīng)變急速增長，此階段剪切帶變形急劇增加。需要注意的是，單元5的y方向正應(yīng)變在剪切位移較小時(s=4 mm)突然增加，然后保持穩(wěn)定，又在s=55 mm時開始激增。其主要原因是：圖13 中測線6 邊界點上方部分土體較早時產(chǎn)生了局部滑移，其邊界點位置也相應(yīng)地發(fā)生了變化。從圖14(c)和(d)可以看出：最大剪應(yīng)變在臨界點前時接近于0，而在臨界點后最大剪應(yīng)變開始急速上升，與正應(yīng)變變化過程相似；體積應(yīng)變也在臨界點前增長緩慢，而在臨界點后激增。由于剪切帶變形先后經(jīng)歷了緩慢增長和急速增長2個階段，并且各個單元的體積應(yīng)變均大于0，剪切帶內(nèi)土體都具有為體積膨脹的特點。

綜上分析可以看出：剪切帶內(nèi)有2 個變形特征，一個是剪切帶變形可分為緩慢增長和急速增長2個階段；另一個是在剪切過程中，剪切帶內(nèi)的土體始終呈現(xiàn)剪脹的特點。

3.5 剪切帶厚度分析

在識別出剪切帶邊界點后，PhotoInfor 軟件可根據(jù)各等參單元自動生成剪切帶不同位置處的厚度。表2所示為等參單元1～5中心點處的剪切帶最終厚度測量結(jié)果。由表2可以看出，剪切帶不同部位處的厚度有所差異，除單元2外，其余的剪切帶厚度都在3 mm以上，且差別較??；但在單元2處，剪切帶厚度為1.27 mm，與其他部分相比差異較大，主要原因是在圖13 中測線3 處剪切帶明顯下降，使得等參單元2較為扁平。最終得到除去單元2厚度后的剪切帶平均厚度為3.50 mm。

表2 剪切帶厚度測量結(jié)果Table 2 Results of shear band thicknessmm

4 結(jié)論

1) 基于數(shù)字散斑相關(guān)原理，提出了一種基于最小二乘法曲線擬合與曲率判斷的土體剪切帶自動識別新方法。針對不同傾角下剪切帶邊界點識別問題，設(shè)計了一種最佳擬合線搜索算法，并將其集成到自行研制的軟件系統(tǒng)PhotoInfor中。

2) 基于土體剪切帶邊界點坐標識別結(jié)果，進一步給出了一種土體剪切帶厚度的簡明計算方法，并分析了砂土在大型直剪試驗中剪切帶的厚度及其空間分布特點。

3) 在土體剪切帶自動識別基礎(chǔ)上，利用PhotoInfor 分析獲得了砂土在大型直剪中剪切帶內(nèi)應(yīng)變的演化特征，發(fā)現(xiàn)剪切帶變形可分為緩慢增長和急速增加2個典型階段，剪切帶內(nèi)土體變形始終呈現(xiàn)剪脹特征。

4) 人工手動所得結(jié)果和程序自動識別的結(jié)果基本一致，但程序自動識別的效率顯然更高，且能夠減少人工識別的人為誤差。