基于位移觀測序列聚類的邊坡潛在滑移區(qū)早期辨識(shí)方法研究

白 宇，孫郡慶，劉明明，李衛(wèi)鵬

(國能新疆托克遜能源有限責(zé)任公司，新疆維吾爾自治區(qū)吐魯番市，830011)

0 引言

邊坡安全穩(wěn)定問題長期困擾著露天煤礦的高效開采，是露天煤礦安全生產(chǎn)不可回避的關(guān)鍵問題[1-2]。近年來，隨著GPS、邊坡雷達(dá)、柔性測斜儀等新型自動(dòng)化監(jiān)測設(shè)備的廣泛應(yīng)用，露天煤礦邊坡變形監(jiān)測技術(shù)得到長足發(fā)展，為主動(dòng)防控邊坡災(zāi)害提供越來越多的現(xiàn)場數(shù)據(jù)支持[3-6]。相比而言，露天煤礦邊坡風(fēng)險(xiǎn)判識(shí)技術(shù)明顯滯后，大量觀測數(shù)據(jù)沒有得到充分挖掘與利用[7-9]。目前，邊坡位移觀測數(shù)據(jù)主要用于臨滑險(xiǎn)情預(yù)警，難以在潛在風(fēng)險(xiǎn)區(qū)的早期識(shí)別方面發(fā)揮作用?；骂A(yù)測缺乏多源監(jiān)測數(shù)據(jù)融合分析手段，通常僅圍繞某個(gè)監(jiān)測點(diǎn)的觀測數(shù)據(jù)進(jìn)行分析，邊坡狀態(tài)判識(shí)存在主觀性和盲目性。

從系統(tǒng)辨識(shí)的角度看，邊坡立體監(jiān)測網(wǎng)不同空間點(diǎn)位處觀測到的監(jiān)測數(shù)據(jù)集蘊(yùn)含邊坡系統(tǒng)在外部應(yīng)力作用下，由無序穩(wěn)態(tài)向有序非穩(wěn)態(tài)轉(zhuǎn)換過程中時(shí)空全域演化的信息[10-12]。為挖掘隱藏在多源監(jiān)測數(shù)據(jù)背后的邊坡系統(tǒng)演化規(guī)律，一些學(xué)者引入數(shù)據(jù)挖掘領(lǐng)域中的聚類分析概念，對邊坡觀測數(shù)據(jù)進(jìn)行特征提取和知識(shí)挖掘。秦雨樵等[13]綜合考慮邊坡各點(diǎn)的位移以及點(diǎn)安全系數(shù)，提出基于K-means聚類算法的滑面搜索方法；雷從雨[14]將邊坡節(jié)點(diǎn)位移聚類分析與變點(diǎn)分析相融合，提出改進(jìn)的特征點(diǎn)位移突變判據(jù)；陳波等[15]采用劃分聚類和投影聚類對邊坡位移監(jiān)測資料進(jìn)行測點(diǎn)區(qū)域劃分和特征提取，以達(dá)到邊坡監(jiān)測數(shù)據(jù)簡約的目的。

作為按觀測時(shí)間先后順序排列的一組數(shù)據(jù)序列，邊坡位移觀測序列表現(xiàn)出十分明顯的非線性、非平穩(wěn)特征，不僅序列內(nèi)部相鄰觀測值存在關(guān)聯(lián)性，而且序列之間存在不等長、偏移性、時(shí)間軸扭曲等結(jié)構(gòu)差異。然而，現(xiàn)有研究工作普遍忽略了位移觀測序列結(jié)構(gòu)的復(fù)雜性，或者僅對觀測序列中特定時(shí)刻的觀測值進(jìn)行聚類[13-14]，或者雖然將整個(gè)觀測序列看作聚類對象，但采用歐氏距離進(jìn)行同步硬性度量[15]，聚類結(jié)果的準(zhǔn)確性和可靠性無法得到保證。為此，筆者提出一種將動(dòng)態(tài)時(shí)間彎曲(DTW)距離和K-means(K均值)聚類相結(jié)合的無監(jiān)督數(shù)據(jù)挖掘方法，從多源監(jiān)測數(shù)據(jù)集中聚類“觀測域上相近、空間域上連續(xù)、變形趨勢明顯”的監(jiān)測點(diǎn)子集，將其所在區(qū)域圈定為潛在滑移區(qū)。實(shí)驗(yàn)結(jié)果表明，上述聚類方法能夠在滑動(dòng)面尚未完全形成階段，提前辨識(shí)出邊坡潛在滑移區(qū)域。

1 邊坡動(dòng)力災(zāi)害演化規(guī)律

露天煤礦邊坡屬于復(fù)雜非線性系統(tǒng)范疇，其動(dòng)力演化遵循一般非線性系統(tǒng)演化所具有的普遍規(guī)律，即從無序到有序的演化歷程[10-11]。在邊坡災(zāi)害的早期階段，各子系統(tǒng)基本保持力學(xué)平衡，整個(gè)系統(tǒng)處于最無序的平衡態(tài)；隨著邊坡高度、暴露面積、暴露時(shí)間的增加，在大氣降水、開采擾動(dòng)等外部應(yīng)力的作用下，局部區(qū)域巖土體單元力學(xué)平衡狀態(tài)被打破，應(yīng)力重新分布，表現(xiàn)出短程相關(guān)和局部有序；在邊坡災(zāi)害的中后期階段，隨著微損傷從小尺度到大尺度的串級(jí)發(fā)展以及主控滑動(dòng)面的逐漸形成，各子系統(tǒng)趨向協(xié)同作用，表現(xiàn)為長程相關(guān)和全局有序，直至整個(gè)邊坡系統(tǒng)發(fā)生突變而形成滑坡。

當(dāng)邊坡子系統(tǒng)發(fā)展到有序階段時(shí)，所處區(qū)域上的監(jiān)測點(diǎn)位移觀測序列必然表現(xiàn)出相似的演化行為，即有序區(qū)域具有“觀測域上相近、空間域上連續(xù)、變形趨勢明顯”的特點(diǎn)。因此，通過邊坡立體監(jiān)測網(wǎng)各測點(diǎn)的時(shí)空域數(shù)據(jù)進(jìn)行聚類分析，可以對露天煤礦邊坡系統(tǒng)的有序程度進(jìn)行定量化表征，為邊坡災(zāi)害的早期辨識(shí)與晚期預(yù)警預(yù)報(bào)提供支撐。

2 位移觀測序列相似性度量方法

邊坡位移觀測序列是按觀測時(shí)間先后順序排列的數(shù)據(jù)序列，可以采用時(shí)間序列相似性度量方法進(jìn)行相似性評(píng)價(jià)。時(shí)間序列相似性度量方法主要分為歐氏距離和DTW距離2種。歐式距離和DTW距離數(shù)據(jù)匹配如圖1所示。

圖1 歐式距離和DTW距離數(shù)據(jù)匹配

圖1中藍(lán)色曲線與紅色曲線代表2個(gè)時(shí)間序列，灰色曲線代表2個(gè)時(shí)間序列元素之間的相似匹配關(guān)系，曲線兩端代表2個(gè)元素。歐式距離是一對一元素相似匹配(藍(lán)色曲線上的一個(gè)元素匹配紅色曲線上的一個(gè)元素)，DTW距離是一對多元素相似匹配(藍(lán)色曲線上的一個(gè)元素有可能匹配紅色曲線上的多個(gè)元素)。

2.1 歐氏距離

歐氏距離(ED)是一種對時(shí)間序列進(jìn)行同步硬性度量(“一對一”匹配)的方法。設(shè)A=(a1，a2，…，an)和B=(b1，b2，…，bn)分別代表2個(gè)具有相同長度的位移觀測序列，那么它們之間的歐氏距離定義為：

(1)

式中：d(A，B)——?dú)W氏距離；

n——位移觀測序列的長度，即觀測時(shí)間的總次數(shù)；

ai——位移觀測序列A中第i個(gè)觀測值；

bi——位移觀測序列B中第i個(gè)觀測值。

d值越小說明序列A和B相似度越高，當(dāng)d=0時(shí)，可認(rèn)為A=B。

2.2 DTW距離

DTW距離是一種對時(shí)間序列進(jìn)行異步相似形態(tài)度量(“一對多”匹配)的方法[16]。假定2個(gè)位移觀測序列A=(a1，a2，…，an)和B=(b1，b2，…，bn)，長度分別為n和m，則序列A和B的DTW距離計(jì)算過程如下。

2.2.1 距離矩陣

構(gòu)造大小為n×m的距離矩陣：

(2)

式中：d(ai，bj)——序列A中第i個(gè)觀測值ai到序列B中第j個(gè)觀測值bj的距離。

2.2.2 彎曲路徑

定義彎曲路徑：

(3)

式中：wr——序列A中第i個(gè)點(diǎn)與序列B中第j個(gè)點(diǎn)之間的匹配關(guān)系。

并且，彎曲路徑需要滿足3個(gè)約束條件：

(1)邊界約束。彎曲路徑起于距離矩陣D的左下角，止于距離矩陣D的右上角，即w1=(1，1)，wk=(n，m)。

(2)單調(diào)性約束。對于相鄰元素wr=(i，j)和wr+1=(i′，j′)，需要滿足i≤i′和j≤j′，即彎矩路徑只能沿著時(shí)間軸單向移動(dòng)。

(3)連續(xù)性約束。對于相鄰元素wr=(i，j)和wr+1=(i′，j′)，需要滿足i′≤i+1和j′≤j+1，即彎曲路徑上任意2個(gè)相鄰的元素在距離矩陣D中也相鄰。

2.2.3 DTW距離

定義序列A和B之間的DTW距離：

DTW(A，B)=min(∑(i，j)∈Wd(ai，bj))

(4)

由式(4)可知，DTW(A，B)對應(yīng)的彎曲路徑為最優(yōu)路徑，在此路徑上的累積距離最小?；趧?dòng)態(tài)規(guī)劃，求解式(2)等價(jià)于構(gòu)造累積距離矩陣Γ，其中的元素服從如下遞歸關(guān)系：

(5)

式中：γ(i，j)——累積距離。

當(dāng)完成累積距離矩陣Γ計(jì)算后，元素γ(n，m)便為序列A和B的DTW距離，即：

DTW(A，B)=γ(n，m)

(6)

進(jìn)一步，以Γ的右上角為起點(diǎn)，即wk=(n，m)，以累積距離最小為選取原則，逆序?qū)ふ襑中各元素，直至Γ的左下角，即w1=(1，1)，可以得到完整的最優(yōu)彎曲路徑。

DTW距離不僅能實(shí)現(xiàn)不等長時(shí)間序列的度量，而且對時(shí)間序列的偏移、形狀縮放、振幅變化、時(shí)間點(diǎn)異常數(shù)據(jù)具有較強(qiáng)的魯棒性。因此，筆者采用DTW距離對位移觀測序列的相似性進(jìn)行分析。

3 位移觀測序列相似性聚類方法

聚類分析一直是數(shù)據(jù)挖掘和機(jī)器學(xué)習(xí)研究的熱點(diǎn)之一，其目的是將數(shù)據(jù)集中的數(shù)據(jù)點(diǎn)按各自的相似度劃分為不同的簇類，并確保簇內(nèi)的數(shù)據(jù)點(diǎn)相似程度最大、簇間數(shù)據(jù)點(diǎn)相似程度最小，從而判識(shí)數(shù)據(jù)集合的結(jié)構(gòu)，提取需要的重要信息[17]。為深度挖掘位移觀測序列背后蘊(yùn)含的邊坡系統(tǒng)演化規(guī)律，采用基于DTW距離的K-means聚類算法，對位移觀測序列進(jìn)行相似性聚類分析。

3.1 算法原理

假定數(shù)據(jù)集X由n個(gè)位移觀測序列組成，Χ=(Χ1，Χ2，…，Χn)，每個(gè)序列含m個(gè)觀測值，Χi=(χi1，χi2，…，χim)。設(shè)定聚類數(shù)k，則聚類算法目標(biāo)是找到一組簇類集合{S1，S2，…，Sk}，使得簇間觀測序列盡可能不同，而簇內(nèi)觀測序列盡可能相似。其中，每個(gè)簇的聚類中心由cj=(cj1，cj2，…，cjm)表示，1≤j≤k。具體K-means算法框架如下。

(1)初始化階段。

步驟1：在數(shù)據(jù)集X中，隨機(jī)選取k個(gè)觀測序列作為初始聚類中心。

(2)聚簇階段。

步驟2：遍歷數(shù)據(jù)集X中所有的觀測序列，根據(jù)各觀測序列到每個(gè)聚類中心的DTW距離，將觀測序列分配到距離最近的聚類中心所代表的簇中，構(gòu)建k個(gè)子集{S1，S2，…，Sk}。

步驟3：計(jì)算每個(gè)子集的位移觀測均值，作為新的聚類中心。

重復(fù)步驟2～3，直到k個(gè)聚類中心不再變化，或達(dá)到指定的迭代次數(shù)。在此過程中，每個(gè)觀測序列有且僅屬于一個(gè)簇，并且與該簇聚類中心之間的距離最小。

3.2 聚類中心計(jì)算

每個(gè)簇的聚類中心由該簇內(nèi)觀測序列各個(gè)時(shí)間點(diǎn)觀測值的均值構(gòu)成，相應(yīng)計(jì)算公式如下：

(7)

式中：cjt——簇Sj聚類中心cj第t個(gè)元素；

χlt——簇Sj內(nèi)觀測序列Χl第t個(gè)時(shí)間點(diǎn)的觀測值；

|Sj|——簇Sj中觀測序列的個(gè)數(shù)。

4 邊坡潛在滑移區(qū)智能辨識(shí)

依據(jù)露天煤礦邊坡系統(tǒng)所表現(xiàn)出的從無序到有序的動(dòng)力演化特征，提出一種基于多源數(shù)據(jù)融合的邊坡潛在滑移區(qū)辨識(shí)方法，通過對監(jiān)測數(shù)據(jù)集的聚類分析，將“觀測域上相近、空間域上連續(xù)、變形趨勢明顯”的監(jiān)測點(diǎn)選定為有序測點(diǎn)，將有序測點(diǎn)所在區(qū)域圈定為潛在滑移區(qū)。具體辨識(shí)流程如下。

(8)

(2)對監(jiān)測數(shù)據(jù)集F進(jìn)行數(shù)據(jù)清洗，消除數(shù)據(jù)集中的缺失值、重復(fù)值、異常值。

(3)設(shè)定聚類數(shù)k，采用基于DTW距離的K-means聚類算法，對F觀測域上的位移序列進(jìn)行相似性度量及聚類分析，將F分解為k個(gè)簇，每個(gè)簇中測點(diǎn)位移觀測序列行為相近。計(jì)算同簇測點(diǎn)位移觀測序列的均值，獲得每個(gè)簇的聚類中心。

(4)對各簇中心進(jìn)行排序，選出[t0，tn]區(qū)間位移變化最大的簇中心，并對該簇測點(diǎn)的空間連續(xù)性進(jìn)行評(píng)估。如果存在連續(xù)性，則將該簇測點(diǎn)選定為有序測點(diǎn)，將該簇測點(diǎn)所在區(qū)域圈定為潛在滑移區(qū)?；禄瑒?dòng)面智能辨識(shí)的流程如圖2所示。

5 實(shí)驗(yàn)研究與結(jié)果分析

采用有限元數(shù)值建模方式仿真模擬典型邊坡漸進(jìn)失穩(wěn)全過程，輸出邊坡不同演化階段的位移場作為實(shí)驗(yàn)測試樣本，對提出的邊坡潛在滑移區(qū)辨識(shí)方法展開適用性研究。

5.1 數(shù)值建模與計(jì)算

采用文獻(xiàn)[18-19]中的邊坡算例：邊坡高度H=20 m，坡角β=26.57°，楊氏模量E=100 MPa，泊松比ν=0.3，容重?=20 kN/m3，黏聚力c=10 kPa，內(nèi)摩擦角φ=20°。土體材料為理想彈塑性材料，強(qiáng)度準(zhǔn)則采用平面應(yīng)變條件下莫爾-庫侖結(jié)合DP準(zhǔn)則(DP4準(zhǔn)則)。邊坡尺寸及網(wǎng)格劃分如圖3所示，邊坡左、右兩側(cè)邊界為法向約束，底邊為雙向固定約束，共劃分四邊形單元5 686個(gè)、節(jié)點(diǎn)5 878個(gè)。采用有限元強(qiáng)度折減方法，按照c/ω，tanφ/ω方式不斷折減土體材料強(qiáng)度參數(shù)，模擬邊坡漸進(jìn)破壞過程及其位移場演化信息，直至邊坡發(fā)生整體失穩(wěn)。其中，ω為折減系數(shù)，可看作廣義時(shí)間，其初始值ω0=1.0，步長△ω=0.002，第m步的折減系數(shù)ωm=ωm-1+△ω。

圖3 邊坡尺寸及網(wǎng)格劃分

折減系數(shù)分別為ω90=1.18、ω150=1.30、ω185=1.36和ω192=1.384時(shí)，邊坡等效塑性應(yīng)變云圖和局部化帶擴(kuò)展路徑如圖4所示。由圖4可以看出，ω90、ω150、ω185和ω192時(shí)刻分別對應(yīng)于邊坡漸進(jìn)演化的早期、中期、晚期和失穩(wěn)階段；當(dāng)邊坡處于前3個(gè)階段時(shí)，雖然滑帶自坡腳向坡體內(nèi)部逐漸擴(kuò)展，但整個(gè)滑帶尚未貫通，滑移區(qū)也沒有完全形成，無法通過塑性云圖和滑帶擴(kuò)展路徑提前辨識(shí)潛在滑移區(qū)；當(dāng)邊坡進(jìn)入失穩(wěn)階段、滑帶完全貫通后，通過塑性云圖和滑帶擴(kuò)展路徑可以準(zhǔn)確獲取滑移區(qū)范圍。

圖4 不同折減系數(shù)工況下邊坡內(nèi)部等效塑性應(yīng)變云圖及局部化帶擴(kuò)展路徑

5.2 監(jiān)測樣本構(gòu)建

將邊坡模型中的所有節(jié)點(diǎn)選做監(jiān)測點(diǎn)，以ωi0為觀測起始時(shí)間，以ωn為觀測截止時(shí)間，從數(shù)值模型計(jì)算結(jié)果中提取觀測數(shù)據(jù)，構(gòu)建邊坡監(jiān)測數(shù)據(jù)集。

(9)

5.3 辨識(shí)算法測試

以觀測起始時(shí)間為分類標(biāo)準(zhǔn)，將監(jiān)測數(shù)據(jù)集細(xì)分為自邊坡漸進(jìn)演化的早期、中期、晚期啟動(dòng)觀測的監(jiān)測數(shù)據(jù)集，按照前文所述方法構(gòu)建監(jiān)測數(shù)據(jù)集。設(shè)定聚類數(shù)k=3，開展監(jiān)測數(shù)據(jù)集觀測域聚類。聚類結(jié)果通過測點(diǎn)空間分布圖和簇中心曲線圖可視化展示。其中，同簇測點(diǎn)及其簇中心曲線使用相同顏色。

5.3.1 邊坡漸進(jìn)演化早期階段啟動(dòng)的觀測數(shù)據(jù)集(簡稱“早期監(jiān)測數(shù)據(jù)集”)

圖5 同簇測點(diǎn)空間分布與簇中心曲線(早期監(jiān)測數(shù)據(jù)集聚類結(jié)果)

5.3.2 邊坡漸進(jìn)演化中期階段啟動(dòng)的觀測數(shù)據(jù)集(簡稱“中期監(jiān)測數(shù)據(jù)集”)

圖6 同簇測點(diǎn)空間分布與簇中心曲線(中期監(jiān)測數(shù)據(jù)集聚類結(jié)果)

5.3.3 邊坡漸進(jìn)演化晚期階段啟動(dòng)的觀測數(shù)據(jù)集(簡稱“晚期監(jiān)測數(shù)據(jù)集”)

5.3.4 辨識(shí)結(jié)果分析

通過對早期、中期和晚期監(jiān)測數(shù)據(jù)集的聚類挖掘，可以得到如下結(jié)論：

(2)在監(jiān)測點(diǎn)數(shù)量足夠多的前提下，邊坡潛在滑移區(qū)辨識(shí)效果將主要取決于觀測的起始時(shí)間和觀測的累積時(shí)長。觀測起始時(shí)間越早，滑移區(qū)可辨識(shí)的時(shí)間也越早，但所需的觀測時(shí)長也越多。以ω90為觀測起始時(shí)間，于ω155時(shí)刻可辨識(shí)出潛在滑移區(qū)，需要65個(gè)時(shí)間步；以ω150為觀測起始時(shí)間，于ω175時(shí)刻可辨識(shí)出潛在滑移區(qū)，需要25個(gè)時(shí)間步；以ω185為觀測起始時(shí)間，滑移區(qū)辨識(shí)僅需要5個(gè)時(shí)間步，但辨識(shí)時(shí)間也最晚。

(3)現(xiàn)場位移監(jiān)測數(shù)據(jù)中的本底噪聲對辨識(shí)效果有很大影響，只有當(dāng)觀測累積位移值明顯大于本底噪聲時(shí)，才有可能對潛在滑移區(qū)進(jìn)行聚類辨識(shí)。因此，選用觀測精度高的監(jiān)測設(shè)備對于潛在滑移區(qū)的早期辨識(shí)十分重要。

6 結(jié)論

(1)邊坡立體監(jiān)測網(wǎng)不同空間點(diǎn)位處觀測到的監(jiān)測數(shù)據(jù)集蘊(yùn)含邊坡系統(tǒng)在外部應(yīng)力作用下，由無序穩(wěn)態(tài)向有序非穩(wěn)態(tài)轉(zhuǎn)換過程中時(shí)空全域演化信息。為深度挖掘隱藏在多源監(jiān)測數(shù)據(jù)背后的邊坡系統(tǒng)演化規(guī)律，提出一種將DTW距離和K-means聚類相結(jié)合的無監(jiān)督數(shù)據(jù)挖掘方法，從多源監(jiān)測數(shù)據(jù)集中聚類“觀測域上相近、空間域上連續(xù)、變形趨勢明顯”的監(jiān)測點(diǎn)子集，將其所在區(qū)域圈定為潛在滑移區(qū)。

(2)采用強(qiáng)度折減有限元數(shù)值模型仿真模擬典型邊坡漸進(jìn)失穩(wěn)全過程，輸出邊坡不同演化階段的位移場作為實(shí)驗(yàn)測試樣本，對潛在滑移區(qū)辨識(shí)方法展開適用性研究。實(shí)驗(yàn)結(jié)果表明，該聚類方法能夠捕捉邊坡系統(tǒng)從無序到有序的漸進(jìn)演化特征，可以在滑動(dòng)面尚未完全形成階段，提前辨識(shí)出邊坡潛在滑移區(qū)域。

(3)需要指出的是，采用數(shù)值模型構(gòu)建監(jiān)測樣本與邊坡現(xiàn)場觀測數(shù)據(jù)之間存在差異，并且將所有節(jié)點(diǎn)都作為監(jiān)測點(diǎn)是理想情況。因此，加快構(gòu)建露天煤礦邊坡地表-深部立體監(jiān)測網(wǎng)，開展基于現(xiàn)場有限測點(diǎn)數(shù)據(jù)的潛在滑移區(qū)早期識(shí)別，將成為未來研究的重點(diǎn)。