陳 喬 張 闊 陳吉龍 張少杰 呂 君? 周儒夫 韋方強,3)
?(礦山地質災害成災機理與防控重點實驗室,西安 710054)
?(中國科學院重慶綠色智能技術研究院,重慶 400714)
??(四川省川建勘察設計院,成都 610041)
??(中國科學院水利部成都山地災害與環境研究所,成都 610041)
???(中國科學院聲學研究所,北京 100190)
???(荊江水文水資源勘測局,湖北荊州 434099)
滑坡是一種全球性的地質災害,它對世界上很多國家都造成了巨大的危害,促使這些國家的研究人員不斷探索滑坡災害監測和預警的技術手段.追溯滑坡監測的研究歷史,學界比較公認的時間起點為日本學者齋藤迪孝[1]在1963 年提出滑坡預報經驗公式時,之后的幾十年各種高新技術在滑坡監測領域中得到應用.在滑坡監測技術的發展過程中,有不少專家學者已經從不同角度以不同觀點進行了系統的總結[2-3],并對其發展趨勢做出了展望.而就目前我國的監測技術,預警預報的成功率仍然較低,如2016 年共發生地質災害9710 起,而全國共成功預報地質災害只有676 起[4],因此開展滑坡監測新技術的探索具有較高的參考價值.
次聲波是一種在自然界非常常見的聲波,其頻率在0~20 Hz 之間,低頻率和長波長的特點使得次聲波具有傳播距離遠、衰減小、抗干擾能力強及穿透能力強等優點,在傳播過程中大氣等介質對它的吸收系數很小,適用于遠距離觀測,被廣泛的用于火山噴發[5-8]、地震[9-11]、泥石流 [12-15] 等自然災害監測中.前期研究表明巖土體在變形破裂過程中因能量釋放必然會產生彈性波和在空氣中傳播的聲波,滑坡過程會產生明顯的次聲波信號[16],因此次聲監測可以作為判斷滑坡的一種技術手段.
近幾年,國內外學者圍繞基于次聲波的滑坡監測技術展開了大量研究.在次聲滑坡監測預警系統研究方面.Li 等[17]基于譜熵和支持向量機對次聲信號進行分類,可用于滑坡的次聲實時監測.任際周等[18-19]從軟硬件設計的角度研究了基于次聲波技術的滑坡監測預警系統,通過對近年來我國數十起重大滑坡災害臨滑預警中收集到的數據進行分析,建立了滑坡預警模型.由于該模型是利用信號分析和統計的方法建立起來的,缺乏對滑坡次聲波產生機理的考慮,因此預警模型的準確度依賴于大量現場數據的統計,而各個區域滑坡因素的多樣性決定了該模型的現場應用具有一定的局限性.
在滑坡破壞過程中次聲波監測機理研究方面.Chai 等[20]通過單軸壓縮試驗發現花崗巖破裂過程中也產生次聲信號.朱星等[21-22]以巖體穩定性監測為目的,對6 種典型的巖石試樣單軸全過程加載試驗過程中的次聲波信號進行了實時監測,并對數字次聲波信號進行了處理與分析.徐洪等[23]對巖石變形破壞次聲異常的能量特征進行分析,研究結果為巖石破壞前兆預警提供了重要依據.賈炳等[24-30]以次聲波技術應用于監測、預警礦井煤巖動力災害為目的,相繼開展了大量的研究,通過進行煤巖試樣的單軸加載試驗,對采集的次聲信號的相對能量和頻率進行分析,發現相對能量和頻率具有明顯的階段性.除了進行實驗室物理模型試驗次聲波的研究外,朱星等[31]還對大光包滑坡側緣拉裂壁巖崩進行次聲監測,并分析了巖崩次聲的基本特征.上述研究結果為巖質滑坡次聲波監測提供了理論參考,但對土質滑坡次聲波監測的研究鮮有報道.基于此,考慮到滑坡的主要力學破壞方式為剪切破壞,開展土剪破壞次聲監測試驗研究,分析了剪切過程中的微觀聲學響應[32-34],探索土質滑坡中次聲產生機理,為實現土質滑坡次聲波監測提供理論參考.
在實驗中,土體直剪實驗力?聲聯合監測系統主要由變尺寸直剪盒、力學數據采集模塊、聲學數據采集模塊以及位移采集模塊四部分組成,如圖1 所示.
傳統的應變控制式直剪儀在進行土體剪切實驗時,采用的土樣規格一般為面積30 cm2高2 cm 的圓柱,可以實時獲取剪切應力與剪切位移數據,但不能采集聲波信號,而在此尺寸下的土樣在剪切過程中所產生的聲波能量又是非常微小,因此本研究需要設計一個大尺寸且尺寸可變的剪切儀器,在本論文中稱之為變尺寸直剪盒.
如圖2 所示,整個變尺寸剪切盒包括上部剪切盒、上部擋土板、固定上部擋土板的螺紋桿、下部擋土板、導向輪(滾珠花鍵)、下部剪切盒、固定上下剪切盒的螺紋桿、固定上部擋土板的螺紋桿,其中的上部剪切盒、上部擋土板、下部剪切盒、下部擋土板是4 個主要結構體.在設計變尺寸直剪盒的各個部件和整體的協調組合時,根據研究條件重點考慮了基本的剪切運動、土樣尺寸可變、防銹蝕和防變形的性能.圖3 為加工廠制作后的實際成品.傳統直剪儀的固有結構不利于安裝導波管放大次聲信號,故本次研究所設計的變尺寸直剪儀更適合土體直剪實驗過程中的次聲?力?位移聯合監測.

圖2 變尺寸剪切盒整體設計圖Fig.2 Overall design of the variable size shear box

圖3 變尺寸剪切儀成品圖Fig.3 Finished product of variable size shear
在實驗中需要對土樣剪切過程中的推力進行實時采集,因此需要一套合適的力學數據采集系統.本研究中使用的力學采集儀器主要為平膜盒式力學傳感器,除此之外,還有配套的數字轉換器以及數據處理軟件等,圖4 為整個系統的示意圖.

圖4 力學數據采集系統示意圖Fig.4 Schematic diagram of mechanical data acquisition system
次聲數據采集模塊包括剛性導波管、次聲傳感器、次聲數據存儲器以及相應的數據處理軟件,其中最核心的部件當屬次聲傳感器.本次實驗所使用的次聲傳感器來自于中國科學院聲學研究所研制的IDS2016 型次聲傳感器,該傳感器具有體積小、靈敏度高、對機械振動不敏感、頻響好的優點.其中3 dB 平穩帶寬保持在0.5~200 Hz 范圍內,靈敏度為50 mV/Pa,在1~100 Hz 頻率范圍內靈敏度曲線平直,可以完全覆蓋剪切過程中產生次聲波的整個頻段.
針對傳統的基于傳感器的接觸式測量方法來測量剪切位移,存在儀器安裝困難的現狀,故結合本研究的實際情況,實驗采用自行研發的非接觸式光學測量方法結合數字圖像識別技術來進行位移數據的采集.
測量的原理是通過計算物體位移變化前后呈現的兩幅圖像之間的相關性來計算位移值.根據本實驗的實際情況,利用此測量方法時需要在上部剪切盒表面貼一張長方形白紙(白色與鋼板的深色在二值化后會形成鮮明反差)作為定位標記,然后使用數碼攝像機錄制整個剪切盒移動過程,得到一串數字圖像序列,再利用計算機程序對每一副圖像進行預處理,包括圖形矯正、二值化處理、選取特征像素行(特征像素行指,圖像序列中的所有圖像中,作為定位標記的像素數量都相等的那一行像素,之所以直接取行像素作為計算量,是因為在本次研究的實驗中,上部剪切盒只做水平直線運動)等步驟(流程如圖5所示).
圖像預處理的每一步需使用相同的變換參數.顯然,預處理后的圖像的特征像素行是一個一維二進制行向量(稱為特征行向量).如果將第一張圖像的特征像素行右移(這里的移位為循環移位,即移位后末位變首位),每次右移后得到新的特征行向量序列,將這一特征行向量序列依次與其他圖像的特征行向量進行相關性計算(見圖6),每次計算后都會得到一組相關性系數,在一組相關性系數中的每一個值都對應一個移位的像素數,而移位的像素數又對應一個位移距離,因此,每一組相關性系數序列中的最大值對應的位移距離就是這幅圖像中被測物體在當前時刻的位移值.
由于位移數據的采集最初來自數碼攝像機,所以本測量方法的精度主要依賴于攝像機的分辨率.本次研究采用的4K 攝像機橫向分辨率為3840 個像素,在白紙所在的焦距處整個屏幕所囊括的總距離在所有實驗中平均約為1000 mm,由此可以算得測量精度為每像素0.26 mm,在之后的研究中可以采用更高分辨率的攝像機來提高測量精度.整個流程可以通過Photoshop 圖像處理軟件和Matlab 編程實現.

圖5 圖像序列第1 幀和第n 幀的處理流程圖Fig.5 Processing flowchart of the first and nth frames of the image sequence

圖6 相關性系數計算示意圖Fig.6 Schematic diagram of correlation coefficient calculation
本次實驗選用的土樣均來自三峽庫區[35-36]重慶市奉節縣,共采集了5 種土類,包括棕壤(黏土)、黃壤(黏土)、黃棕壤(黏土)、粗骨土(含礫石黏土)、紫色土(含礫石黏土),部分土樣如圖7 所示.其中棕壤取自云霧鄉羅圈村,黃壤取自安坪鄉新浦滑坡后緣,黃棕壤取自長安鄉保安村,粗骨土取自草堂鎮廟灣,紫色土取自康樂鎮平皋村,取樣點位置見圖8,重塑土樣是原狀的黃壤土經烘干、碾碎而成,各類土的相關物理力學性質見表1.

圖7 土樣類別圖Fig.7 Types of soil samples
首先,將野外取回的所有土樣分批放入烘箱中,以105?C 的溫度持續24 h 鼓風干燥,并將烘干后的土樣研磨并過2 mm 篩備用; 接下來,按具體實驗目的,將變尺寸剪切儀的上部擋土板和下部擋土板固定在規定的位置,采用分層壓實的方法將處理好的土樣堆填于剪切盒中;然后,在土樣頂部插入剛性導波管,在導波管頂部安裝次聲傳感器探頭,將力學傳感器放置在千斤頂與上部剪切盒之間,將記錄位移用的攝像頭擺放在合適位置;最后,記錄實驗開始時間,并打開數據采集軟件,在記錄30 s 背景噪音的基礎上,開始按壓千斤頂為上部剪切盒加載直至土樣剪切破壞.

圖8 土樣取樣點位置圖Fig.8 Location of soil sampling points

表1 各類土的相關物理力學性質Table 1 Related physical and mechanical properties of various soils
3.2.1 背景噪聲分析與處理
在室內試驗過程中,會出現很多外在條件干擾的情況.例如開、關門和抽風機運轉時會導致空氣不穩定流通,實驗用水所發出的超低頻聲等.針對上述可控的低頻聲源,本試驗時間均選定在凌晨,關閉與試驗無直接關系的電器設備,以減少外界帶來的干擾.在預處理的基礎上,考慮到本次研究中關注的重點是聲波信號的次聲部分,所以需要對聲波信號進行了低通濾波處理.在文中利用巴特沃斯濾波器對信號進行處理,設置通帶截止頻率為20 Hz,阻帶截止頻率為25 Hz,通帶內最大衰減為0.5 dB,阻帶衰減為40 dB.同時,利用背景噪聲強度較低的現狀,設置聲壓強度閾值為0.000 5 Pa,最后,對閾值處理后小波信號進行小波逆變換.
3.2.2 次聲事件拾取
為了提取次聲事件,將整個時域按次聲事件分段,本文提出了一種按次聲事件由計算機程序自動劃分時域的方法.首先根據每個采樣點的信噪比初步判斷該采樣點是否位于一個次聲事件的時段內,再歸并比較密集的有效采樣點形成一個次聲事件,具體步驟如下.
第一步,初步判斷一個采樣點是否位于一個次聲事件的時段內.
根據離散聲波信號的能量公式

其中,E表示信號能量;s表示聲波通過的面積,可以認為是常量; ρ 表示聲波傳播介質的密度,這里即指空氣的密度;v表示音速;I表示采樣點總數;x(i) 表示第i個采樣點的聲壓值,可以得到單個采樣點所代表的能量

根據功率的定義(單位時間內的能量),本文定義離散聲波信號中第i個采樣點的功率為

本文對“均方功率”的定義由式(1)右側去掉s,ρ 和v再開方得到,因為式中的s,ρ 和v為常量,所以去掉這3 個量不影響信號點間的相對取值且避免了對它們的測量,再開方的目的是為了使最后得到的值與聲壓值量綱相同.

由式(3)定義一個判據式

其中,Pj表示第j個采樣點的功率,Pmax表示前30 s背景噪聲的最大功率,L表示用于計算信噪比的時窗長度.
式(4)左側的意義為以第i個采樣點為中心的L時長范圍內的信號能量,右側的意義為在長度同樣為L的時段內噪聲的最大能量,整個式子表示將局部的信號能量與相同長度噪聲能量進行比較.針對新信號中的每一個點代入式(5) 中判斷,結果為真時,把該點標記為1,否則標記為0,這樣就形成一個由0 和1 組成的序列,稱為第一輪標記序列,其長度等于原始信號.于是整個時域按信噪比被初步劃分出來.此時的劃分是不夠的,因為第一輪標記序列中的1 值點過于零碎,需要將間隔很短的信號段再一次歸并.
第二步,對第一輪標記序列進行歸并整合.這一步會用到一個所有值均為1 的行向量r作為參考向量,把這個參考向量作為時窗在第一輪標記序列的時域中平移,并與和時窗長度相同的局部按下式做計算

其中,ci表示判別該采樣點(第i個點)是否需要重新標記為1 或0 的值,J表示移動時窗的采樣點序列總數,sj表示第一輪標記序列中的第j個點,??表示向上取整,??表示向下取整.
由式(6)計算得到結果ci后,再判別是否需要重置采樣點的聲壓值,當ci> 0 時,說明第一輪標記序列中與時窗長度相同的局部向量與作為參考樣本的向量相互匹配,該點可以標記為1,否則當ci< 0 時應該標記為0.最后再次得到一個基于時域的由若干0 和1 分別連續分布的新序列,稱為次聲事件標記序列,其中連續的1 指明該時段內存在能量明顯大于背景噪音的次聲信號,而連續的0 指明該時段內不存在有效的次聲信號.
式(5)中的時窗長度L和式(6)中時窗長度J是兩個對最終劃分結果影響較大的計算因子.其中L應設置一個較小值,若設置值過大,一些能量較小的次聲事件將不會被識別出.而時窗長度J影響最終次聲事件劃分出來的數量,若J取值過小,次聲事件會劃分的過密,無疑會增加后期分析的工作量,而J若取值過大,對次聲事件的劃分又太過于簡單,一些間隔很遠的有效信號也被歸并為一個次聲事件.經過反復實驗,本文認為L=5,J=100 是合適的.最后根據式(5)和式(6),使用MATLAB 編制了對信號按次聲事件自動劃分時域的程序.
本次共完成了不同尺寸、不同土性和不同飽和度共13 組土剪切實驗(見表2).圖9 展示了所有實驗的“推力?時間”曲線和聲波信號時域圖土樣的剪切力與位移原始關系.從推力變化的整體趨勢來看,每組推力都隨時間先陡然升高,超過峰值后逐漸降低.通過對推力進行規格化處理,統一在0~7000 N 范圍內,可以進一步反映推力在時域上的變化趨勢.以具有代表性的sy3 實驗為例(如圖10),圖中實線為實測推力值,推力曲線呈現鋸齒狀是因為本次實驗加載方式采用的是千斤頂人工加載,所以導致加載不連續,反映出了該類土的應力松弛特點.
人工加載對剪切速度的影響可以認為是微乎其微的,在本實驗中推力與剪切帶受到的剪力幾乎相等,由牛頓第二定律

其中,f是土體受到的剪力,單位kN;F是千斤頂的推力,單位kN;m是土體上半部分質量,實測sy3 的土體總質量49.24 kg;a為位移加速度,單位為m/s2.經計算ma的值在?0.046 6 到0.049 9 之間,遠遠小于推力大小,因此可大致認為在千斤頂加載的過程中上部剪切盒保持勻速運動.

表2 實驗組別一覽Table 2 List of experimental groups

圖9 “推力?時間”曲線和聲波信號時域圖Fig.9 “Thrust-time”curve and time domain diagram of acoustic signal

圖10 sy3 推力?位移曲線Fig.10 Thrust-displacement curve of sy3
同時,推力的趨勢線隨位移變化呈現升高后降低的規律,當sy3 發生剪切破壞時,首先在某一點處超過抗剪強度而發生破裂,然后從此破裂點開始向整個剪切面蔓延擴展,隨著剪切面的蔓延擴展,整個土體的抗剪強度也在降低,應變局部化現象增強,直至整個剪切面貫通,整體表現出明顯的應變軟化特性.
根據漸進破壞理論和推力趨勢線的形態,本文將其劃分為了彈性變形階段、彈塑性變形階段以及峰值后的應變軟化階段3 個階段.為了探究與聲學信號的關系,將推力?位移曲線在時域上重新描述,繪制出的“推力?時間”曲線圖,以sy3 為例(見圖11),推力趨勢線在經過前30 s 的背景噪音錄制階段后隨時間逐漸增大,加載后先是呈直線狀態增大,處于彈性變形階段;在74 s 后曲線斜率開始減小,曲線本身呈向下彎曲狀,處于彈塑性變形階段; 在118 s 處推力趨勢線達到峰值6.1 kN 后開始減少,隨后隨著時間推移逐漸減小,對應著應變軟化階段.
sy3 在次聲事件拾取的基礎上獲取了26 個次聲事件(見圖11),將每個事件對應的次聲信號均方功率峰值用虛線連接起來形成均方功率的包絡線,隨著時間推移包絡線開始逐漸上升,在105 s 處均方功率達到0.047 9 Pa 后開始下降,之后偶有起伏但都沒有超過105 s 處的峰值,并且很明顯的,該包絡線的時域形態與推力趨勢線較為相似.

圖11 sy3 均方功率圖Fig.11 sy3 mean square power diagram
從圖11 來看,剪切試驗過程中的時域信號在105 s 處能量達到了峰值,同時應力曲線也達到峰值,實際上從現場結果來看試樣在達到應力峰值時已經破壞,但此時應力并沒有大幅度下降,表明土樣還有一定殘余強度.破壞過程中的次聲信號強度比破裂后剪切面滑移所產生的次聲信號強度大,說明土樣內部裂隙發展所產生的次聲信號與土樣顆粒摩擦所產生的次聲信號有一定的區別.
在此基礎上,將所有組實驗的均方峰值與推力峰值之間的時間差進行計算,從結果來看(見表3 和圖12),次聲信號的均方功率峰值都在推力峰值之前,提前的時間分布在2.30~47.11 s,平均值為22.95 s.該結論可為土質滑坡次聲預警模型研究提供理論參考,具有重大意義和工程應用價值.

表3 各組實驗次聲事件次數與時間差計算結果Table 1 Calculation results of the number and time difference of infrasound events in each group of experiments

圖12 所有實驗均方功率圖Fig.12 Mean square power plots for all experiments
在“推力?位移” 曲線的3 個階段中,彈性變形階段推力趨勢線呈近直線型,斜率幾乎不變,彈塑性變形階段趨勢線的增量逐漸降低,線形向下凹,應變軟化階段趨勢線經過峰值后逐漸降低,但較彈性變形階段下降趨勢更為平緩.根據應變軟化類黏性土在剪切破壞中滿足漸進性破壞理論,剪切面破壞面的形成呈漸進性,另外參考一些對大型直剪實驗剪切帶形態的相關研究[37],本文認為彈性變形階段對應于圖13 中的剪切前到初始剪切之間的狀態,彈塑性變形階段和應變軟化階段對應初始剪切到剪切完成之間的狀態.

圖13 剪切過程示意圖(陰影區域表示壓密區,1#、2#和3#三根粗直線表示鋁絲,用于表征土樣平行于剪切方向的截面的形變)Fig.13 Schematic diagram of the cutting process
在剪切過程初期,土體整體變形以上部試樣中壓密區土體的彈性變形為主,剪切帶中的剪切變形為輔,因此推力趨勢線呈近直線狀,在此階段,土顆粒主要發生彈性變形,顆粒間接觸面相互擠壓產生的聲波能量較小,但隨著剪切過程的發展,土體整體變形進入到以剪切帶塑性變形為主的階段,此階段中主要以黏土膠體團粒本身的變形為主,這種膠體團粒的變形比彈性變形所產生的聲波能量要大,從以土顆粒彈性變形為主到黏土團粒本身變形為主,后者產生的聲波信號成分也逐漸超過前者,所以在宏觀上表現為彈塑性階段次聲信號幅度值比彈性階段更大(見圖14),均方功率峰值一般也發生在此階段.在推力峰值時,土體剪切面開始漸進性破壞,并進入到應變軟化階段,在此階段剪切面從某一點發生破裂并逐步蔓延到整個剪切面,此時剪切帶的運動形式在微觀上發展成以黏土團粒拉張破裂后土顆粒間的摩擦為主,黏土團粒的破裂為輔,此時的信號能量也逐漸減弱,偶爾局部團粒破裂活動較為集中時,在時域信號上還會產生相應的脈沖,但包含此脈沖的次聲事件的幅值一般不會再超過彈性階段的幅值.

圖14 不同階段聲波幅值變化情況Fig.14 Changes in sound wave amplitude at different stages
(1)針對常規土體應變式直剪儀用于聲發射實驗時聲學信號微弱的缺陷,研發了專門用于探測次聲信號的變尺寸直剪儀,設計了非接觸式位移數據的自動采集模塊,搭建了一套大尺寸土體直剪實驗力?聲?位移聯合采集系統,為土質滑坡次聲監測技術室內實驗研究奠定了基礎.
(2) 將聲學信號結合力學和位移數據進行分析,發現次聲信號在時域上的脈沖形態與直剪實驗中的“推力?時間” 曲線具有明顯的相關性,并在應變軟化類黏性土直剪過程中剪切帶的漸進性破壞理論的基礎上,分析了剪切過程中彈性階段、彈塑性階段和應變軟化階段的微觀聲學響應,認為黏性土剪切破壞時發出的次聲信號主要來源于彈性階段中黏土顆粒的相互擠壓、彈塑性階段和應變軟化階段中由黏土顆粒構成的膠體團粒本身的變形以及黏土膠體團粒拉張破裂后包含其中的土顆粒間的相互摩擦三種形式.
(3)在土剪不同階段,次聲信號的幅值大小呈現出應變軟化階段< 彈性階段< 彈塑性階段的規律,同時,次聲信號的均方功率峰值都在推力峰值之前,平均提前為22.95 s,這說明將次聲技術作為土質滑坡預警手段的具有理論可行性.
本研究對于土質邊坡的動態監測和穩定性預測都具有重要的指導意義,可為進一步基于次聲的土質滑坡預警模型研究提供一定的理論參考.