楊 旭，楊 旭，李 佳，王建國(guó)

(1.北京航空航天大學(xué)電子信息工程學(xué)院，北京 100089；2.北京航空航天大學(xué)合肥創(chuàng)新研究院，合肥 230000)

0 引言

山體滑坡災(zāi)害具有隱蔽性強(qiáng)、突發(fā)性強(qiáng)、破壞力大的特點(diǎn)。傳統(tǒng)的山體滑坡監(jiān)測(cè)方法主要有宏觀地質(zhì)觀測(cè)法、簡(jiǎn)易觀測(cè)法、設(shè)站觀測(cè)法、儀表觀測(cè)法等，它們普遍存在的問(wèn)題是需要人工定期到現(xiàn)場(chǎng)進(jìn)行數(shù)據(jù)的采集，缺少實(shí)時(shí)性，監(jiān)測(cè)精度較低，并且對(duì)滑坡的預(yù)測(cè)預(yù)警能力差[1-3]。

近年來(lái)，在地質(zhì)形變監(jiān)測(cè)技術(shù)不斷發(fā)展和完善的過(guò)程中，大地精密測(cè)量法、全球衛(wèi)星導(dǎo)航系統(tǒng)(Global Navigation Satellite System，GNSS)高精度地質(zhì)形變監(jiān)測(cè)法、近景攝影測(cè)量法、時(shí)間域反射測(cè)試技術(shù)、遙感監(jiān)測(cè)技術(shù)等方法常被用于地質(zhì)形變監(jiān)測(cè)中。目前主流的山體滑坡監(jiān)測(cè)方法是通過(guò)合成孔徑雷達(dá)(Synthetic Aperture Radar，SAR)技術(shù)進(jìn)行監(jiān)測(cè)，但是該方法受天氣狀況和植被覆蓋的影響大，實(shí)時(shí)性差，難以實(shí)現(xiàn)長(zhǎng)期監(jiān)測(cè)[4-6]。

實(shí)時(shí)動(dòng)態(tài)(Real-Time Kinematic，RTK)技術(shù)可實(shí)現(xiàn)實(shí)時(shí)定位，并且在靜態(tài)情況下的定位精度可達(dá)mm級(jí)。與傳統(tǒng)GNSS形變監(jiān)測(cè)技術(shù)相比，RTK定位技術(shù)具有十分顯著的優(yōu)勢(shì)。RTK技術(shù)被廣泛應(yīng)用于安全監(jiān)測(cè)、土地勘探、城市布局、水利等領(lǐng)域[7-11]。在我國(guó)的金川二礦區(qū)地裂縫[9]、黃土坡臨江I號(hào)崩滑體[10]等地質(zhì)災(zāi)害的變形監(jiān)測(cè)中，均采用了高精度的全球定位系統(tǒng)(Global Positioning System，GPS)靜態(tài)相對(duì)定位技術(shù)，監(jiān)測(cè)精度可以達(dá)到mm級(jí)；并逐漸開(kāi)始嘗試采用GPS定位技術(shù)監(jiān)測(cè)滑坡體的垂直形變，發(fā)現(xiàn)在基線較短、觀測(cè)時(shí)間較長(zhǎng)和觀測(cè)條件較好的情況下，其精度也可以達(dá)到mm級(jí)[11]。

作為形變監(jiān)測(cè)中的重要數(shù)據(jù)處理算法，Kalman濾波算法廣泛地應(yīng)用于建筑、橋梁、地質(zhì)災(zāi)害的形變監(jiān)測(cè)中，并取得了非常不錯(cuò)的監(jiān)測(cè)效果[12-15]。C.D.Ince等通過(guò)GPS定位技術(shù)和Kalman濾波方法對(duì)土耳其的伊斯坦布爾和馬爾馬拉進(jìn)行形變監(jiān)測(cè)，并取得預(yù)期的監(jiān)測(cè)效果[12]；朱建軍等使用變形監(jiān)測(cè)數(shù)據(jù)通過(guò)Kalman濾波模型動(dòng)態(tài)分析鏈子崖危巖體位移變形規(guī)律[13]；高雅萍、張勤使用自適應(yīng)Kalman濾波在地理信息系統(tǒng)(Geographic Information System，GIS)滑坡監(jiān)測(cè)中也取得了預(yù)期監(jiān)測(cè)效果[14]。

影響山體滑坡的因素復(fù)雜多變，有已知因素和未知因素，具有典型的灰色特征?；疑A(yù)測(cè)理論在解決信息不完全的系統(tǒng)問(wèn)題時(shí)具有十分顯著的優(yōu)勢(shì)。對(duì)于具有這樣特征的系統(tǒng)，可使用灰色預(yù)測(cè)理論對(duì)其進(jìn)行深入的分析[16-19]。李曉鴿采用灰色預(yù)測(cè)理論對(duì)魯?shù)乩娬緦?duì)外交通公路滑坡監(jiān)測(cè)數(shù)據(jù)進(jìn)行了擬合和預(yù)測(cè)研究，得到了很好的預(yù)測(cè)效果[16]；王旭昭等使用改進(jìn)的灰色預(yù)測(cè)理論對(duì)新灘滑坡和黃茨滑坡進(jìn)行預(yù)測(cè)分析，取得了良好預(yù)測(cè)效果[17]；此外，K. Yin[18]、王東岳[19]等均使用灰色預(yù)測(cè)理論進(jìn)行了山體滑坡預(yù)測(cè)，并取得良好預(yù)測(cè)效果。

針對(duì)當(dāng)前的山體滑坡監(jiān)測(cè)技術(shù)精度低、實(shí)時(shí)性差、缺少完善的預(yù)測(cè)預(yù)警機(jī)制的問(wèn)題，本文提出了一種基于抗差自適應(yīng)Kalman濾波和灰色預(yù)測(cè)理論的山體滑坡監(jiān)測(cè)技術(shù)。該技術(shù)通過(guò)抗差自適應(yīng)Kalman濾波進(jìn)行包括RTK定位技術(shù)、土工帶傳感器技術(shù)、無(wú)人機(jī)近景攝影測(cè)量技術(shù)在內(nèi)的多源數(shù)據(jù)融合，然后通過(guò)灰色預(yù)測(cè)理論與蠕變切線角判據(jù)實(shí)現(xiàn)對(duì)山體滑坡的預(yù)測(cè)預(yù)警功能。

1 抗差自適應(yīng)Kalman濾波

1.1 抗差自適應(yīng)Kalman濾波[12-15]

在目標(biāo)系統(tǒng)中，設(shè)Xi為由一組隨機(jī)變量組成的向量，稱為該系統(tǒng)在i時(shí)刻的狀態(tài)向量；設(shè)Zi為由一組觀測(cè)變量組成的向量，稱為該系統(tǒng)在i時(shí)刻的觀測(cè)向量。當(dāng)目標(biāo)系統(tǒng)為離散時(shí)間系統(tǒng)時(shí)，系統(tǒng)的狀態(tài)向量隨時(shí)間的變化規(guī)律可以通過(guò)狀態(tài)方程描述，目標(biāo)系統(tǒng)的狀態(tài)向量和觀測(cè)向量的關(guān)系可以通過(guò)觀測(cè)方程描述

(1)

其中，Xi-1、Xi分別是目標(biāo)系統(tǒng)在i-1和i時(shí)刻的狀態(tài)向量；Zi是目標(biāo)系統(tǒng)在i時(shí)刻的觀測(cè)向量；U稱為控制向量；W稱為系統(tǒng)動(dòng)態(tài)噪聲向量；V稱為系統(tǒng)觀測(cè)噪聲向量；A稱為系統(tǒng)從i-1時(shí)刻到i時(shí)刻的狀態(tài)轉(zhuǎn)移矩陣；B稱為控制向量系統(tǒng)矩陣；H稱為觀測(cè)值系統(tǒng)矩陣。

抗差自適應(yīng)Kalman濾波技術(shù)將抗差估計(jì)和自適應(yīng)估計(jì)相結(jié)合，以減小觀測(cè)誤差和系統(tǒng)建模誤差對(duì)濾波結(jié)果的影響，其技術(shù)原理可以通過(guò)如下方程加以描述：

時(shí)間更新方程

(2)

Pi/i-1=APi-1/i-1AT+Q

(3)

狀態(tài)更新方程

(4)

(5)

(6)

1.2 觀測(cè)噪聲等價(jià)協(xié)方差矩陣的構(gòu)造

典型的觀測(cè)噪聲等價(jià)協(xié)方差矩陣D的構(gòu)造方法[14]有：Huber構(gòu)造法、丹麥法、IGG法、抗差貝葉斯估計(jì)法等。本文選用IGG1構(gòu)造法，具體內(nèi)容如下

(7)

其中，一般取k0=2.5，k1=3.5。Vi為觀測(cè)殘差，表達(dá)式如下

(8)

1.3 自適應(yīng)因子的構(gòu)造

本文采用兩段函數(shù)模型[14]進(jìn)行自適應(yīng)因子的構(gòu)造

(9)

式中，c為常量，一般取c=1～2.5，0

2 灰色預(yù)測(cè)理論

2.1 灰色預(yù)測(cè)模型(GM模型)

2.1.1 基本GM(1,1)模型

設(shè)X(0)為含有n個(gè)元素的序列，其中x(0)(k)≥0,k=1,2,…,n。上角標(biāo)處括號(hào)內(nèi)的數(shù)字為序列模型的階數(shù)，第二個(gè)括號(hào)內(nèi)的數(shù)字為序列中元素的序號(hào)。對(duì)序列X(0)進(jìn)行一次累加生成得到它的1-AGO序列X(1)。對(duì)序列X(1)進(jìn)行一次緊鄰均值生成得到序列Z(1)。稱方程x(0)(k)+ax(1)(k)=b為GM(1,1)模型。其中a為發(fā)展系數(shù)，b為灰色作用量。令

(10)

Y=[x(0)(2),x(0)(3),…,x(0)(n)]T

(11)

(12)

則有

(13)

方程x(0)(k)+ax(1)(k)=b的解

(14)

取x(1)(0)=x(0)(1)，則

(15)

2.1.2 新陳代謝GM(1,1)模型

2.2 精度檢驗(yàn)

(16)

(17)

在實(shí)際應(yīng)用中，根據(jù)實(shí)際情況設(shè)定閾值A(chǔ)0，當(dāng)A

(18)

在實(shí)際應(yīng)用中，根據(jù)實(shí)際情況設(shè)定閾值P0，當(dāng)P>P0時(shí)，則認(rèn)為該模型為小概率誤差合格模型。關(guān)于均值方差比閾值A(chǔ)0和小概率誤差閾值P0的設(shè)定準(zhǔn)則與模型精度等級(jí)的關(guān)系如表1所示。

表1 模型精度等級(jí)參考[19]

3 蠕變切線角判據(jù)

通過(guò)大量的山體滑坡實(shí)例可知，山體滑坡的形變過(guò)程可分為3個(gè)時(shí)期：初始變化時(shí)期、穩(wěn)定變化時(shí)期和劇烈變化時(shí)期。劇烈變化時(shí)期進(jìn)一步細(xì)分為：初加速時(shí)期、中加速時(shí)期和臨破壞時(shí)期[19,21-22]。如圖1所示。

圖1 山體滑坡演化過(guò)程中的形變曲線

本文基于上述分析建立蠕變切線角預(yù)警判據(jù)[20]。通過(guò)一定時(shí)間的觀測(cè)得到穩(wěn)定變化時(shí)期監(jiān)測(cè)點(diǎn)的形變速度v0，將穩(wěn)定變化時(shí)期的山體形變運(yùn)動(dòng)看作是勻速運(yùn)動(dòng)，可知v0是定值，定義

(19)

其中，S(k)為從起始監(jiān)測(cè)時(shí)刻到當(dāng)前監(jiān)測(cè)時(shí)刻監(jiān)測(cè)點(diǎn)處的累計(jì)位移量；T(k)為具有時(shí)間量綱的縱坐標(biāo)值。由T(k)和t(k)繪制T-t曲線，如圖2所示，通過(guò)T-t曲線可以得到切線角的表達(dá)式

圖2 經(jīng)過(guò)處理后得到的T-t曲線

(20)

其中，αk為切線角；t(k)為第k個(gè)監(jiān)測(cè)時(shí)刻；Δt為采樣時(shí)間間隔。

根據(jù)行業(yè)內(nèi)的長(zhǎng)期觀測(cè)研究得到山體滑坡過(guò)程的預(yù)警判據(jù)表，如表2所示。該預(yù)警判據(jù)表綜合考慮了山體滑坡監(jiān)測(cè)點(diǎn)的形變速率V和切線角α包含的信息內(nèi)容，對(duì)山體滑坡預(yù)警等級(jí)進(jìn)行了5層劃分，分別代表不同的山體滑坡可能性和危險(xiǎn)性，從而實(shí)現(xiàn)對(duì)山體滑坡的等級(jí)預(yù)警功能。

課文大意：《綜合英語(yǔ)》第四冊(cè)中第七單元the Monster是一篇典型的描寫(xiě)文，栩栩如生地刻畫(huà)了十九世紀(jì)德國(guó)劇作家瓦格納的人物形象。作為劇作家，瓦格納是一個(gè)天才，畢生創(chuàng)作十三部歌劇和音樂(lè)劇，其中十一部仍繼續(xù)上演，八部被列入世界音樂(lè)名劇之列。但現(xiàn)實(shí)中，瓦格納可以稱得上是一個(gè)怪物，他病態(tài)，神經(jīng)質(zhì)，喜歡妄自尊大，絲毫不具責(zé)任感，喜新厭舊，甚至偷人妻子，賴賬不還。(張伊娜,216)這就是本文的主人公，一個(gè)集所有缺點(diǎn)于一身的怪物，同時(shí)他又是音樂(lè)界的天才，具備超凡的稟賦。作者這種先抑后揚(yáng)的寫(xiě)法極好地刻畫(huà)了這個(gè)杰出的人物。

表2 山體滑坡預(yù)警判據(jù)表[20]

4 實(shí)驗(yàn)測(cè)試

4.1 滑坡監(jiān)測(cè)系統(tǒng)

本文構(gòu)建以RTK為核心的空天地一體化立體監(jiān)測(cè)網(wǎng)絡(luò)，實(shí)現(xiàn)對(duì)滑坡災(zāi)害的監(jiān)測(cè)評(píng)估。

滑坡監(jiān)測(cè)系統(tǒng)如圖3所示，使用RTK進(jìn)行地質(zhì)沉降監(jiān)測(cè)，在監(jiān)測(cè)區(qū)域布置土工帶傳感器，采集監(jiān)測(cè)區(qū)域的位移信息并上傳到云端以待處理，利用無(wú)人機(jī)攝影測(cè)量技術(shù)提供空域增強(qiáng)信息，輔助RTK進(jìn)行滑坡形變監(jiān)測(cè)。北斗短報(bào)文通信技術(shù)可以實(shí)現(xiàn)全天候、全覆蓋域、高可靠性的通信，其通信時(shí)延約為0.5s，點(diǎn)對(duì)點(diǎn)通信時(shí)延約為1～5s，同時(shí)該技術(shù)的設(shè)備成本低，適用于廣泛的環(huán)境監(jiān)測(cè)活動(dòng)，所以RTK形變監(jiān)測(cè)信息通過(guò)北斗短報(bào)文進(jìn)行傳遞。

圖3 山體滑坡監(jiān)測(cè)系統(tǒng)

設(shè)計(jì)后處理軌跡平滑算法剔除異常測(cè)量信息，使用數(shù)據(jù)融合[23]技術(shù)對(duì)RTK定位數(shù)據(jù)、土工帶傳感器數(shù)據(jù)、無(wú)人機(jī)三維成像數(shù)據(jù)進(jìn)行綜合分析后，在滑坡監(jiān)測(cè)點(diǎn)處計(jì)算得到位移、速度、加速度測(cè)量信息。本文選用抗差自適應(yīng)Kalman濾波技術(shù)實(shí)現(xiàn)數(shù)據(jù)的融合分析，獲得精確的形變監(jiān)測(cè)結(jié)果。在獲得經(jīng)過(guò)融合處理的形變監(jiān)測(cè)數(shù)據(jù)后，通過(guò)灰色預(yù)測(cè)理論對(duì)未來(lái)時(shí)刻的形變數(shù)據(jù)進(jìn)行合理預(yù)測(cè)，得到未來(lái)時(shí)刻監(jiān)測(cè)點(diǎn)處的形變信息，最后結(jié)合蠕變切線角判據(jù)實(shí)現(xiàn)對(duì)山體滑坡災(zāi)害的預(yù)測(cè)預(yù)警功能。

設(shè)備布設(shè)的具體情況如圖4所示。在重點(diǎn)監(jiān)測(cè)區(qū)域進(jìn)行監(jiān)測(cè)設(shè)備布設(shè)時(shí)，為合理兼顧設(shè)備成本和監(jiān)測(cè)效果，主要將監(jiān)測(cè)設(shè)備布設(shè)在可能滑坡體的頂部和底部，尤其是一些裂縫處?？紤]到皖南山區(qū)植被覆蓋程度大、連續(xù)性強(qiáng)，植被覆蓋斷裂痕跡處也需重點(diǎn)布設(shè)監(jiān)測(cè)設(shè)備。RTK接收機(jī)安置在網(wǎng)絡(luò)的中心位置。該土工帶傳感器與RTK接收機(jī)的排布網(wǎng)絡(luò)與無(wú)人機(jī)攝影測(cè)量得到的形變監(jiān)測(cè)數(shù)據(jù)相結(jié)合，三種形變監(jiān)測(cè)數(shù)據(jù)在空間上具備了數(shù)據(jù)融合的條件。

圖4 監(jiān)測(cè)區(qū)域數(shù)據(jù)采集設(shè)備布設(shè)

4.2 滑坡監(jiān)測(cè)模型

對(duì)山體滑坡進(jìn)行監(jiān)測(cè)分析時(shí)，監(jiān)測(cè)點(diǎn)的形變可以使用運(yùn)動(dòng)模型進(jìn)行分析。運(yùn)動(dòng)模型的描述方程如下

(21)

(22)

位移x(i)通過(guò)土工帶傳感器直接測(cè)得

x(i)=dLtu(i)

(23)

(24)

加速度a(i)通過(guò)定期無(wú)人機(jī)攝影測(cè)量形變量和飛行間隔時(shí)間計(jì)算得到

(25)

為保證監(jiān)測(cè)系統(tǒng)的實(shí)時(shí)性，在平穩(wěn)形變階段，監(jiān)測(cè)區(qū)域加速度信息基本穩(wěn)定，無(wú)人機(jī)攝影測(cè)量周期TUAV可設(shè)置較大，TUAV≈30d≈15T。在加速形變易發(fā)生滑坡階段，應(yīng)按照實(shí)際情況進(jìn)行調(diào)整，對(duì)于監(jiān)測(cè)價(jià)值極高的區(qū)域，無(wú)人機(jī)攝影測(cè)量周期TUAV需大幅縮短；對(duì)于監(jiān)測(cè)價(jià)值一般的區(qū)域，可屏蔽無(wú)人機(jī)攝影測(cè)量數(shù)據(jù)，僅僅使用土工帶傳感器測(cè)量數(shù)據(jù)和RTK接收機(jī)測(cè)量數(shù)據(jù)進(jìn)行重點(diǎn)監(jiān)測(cè)，此時(shí)將加速度信息作為一種隨機(jī)噪聲進(jìn)行處理。在損失一定監(jiān)測(cè)精度的情況下降低監(jiān)測(cè)成本。

4.3 數(shù)據(jù)采集

4.3.1 RTK定位數(shù)據(jù)

RTK定位數(shù)據(jù)包含RTK接收機(jī)的經(jīng)度、緯度、天線高度、接收機(jī)的工作模式和觀測(cè)衛(wèi)星數(shù)量等信息。通過(guò)RTK定位得到監(jiān)測(cè)點(diǎn)的經(jīng)度、緯度和天線高度，則監(jiān)測(cè)點(diǎn)位移量dLRTK與豎直位移量dL1和水平位移量dL2滿足如下關(guān)系

(26)

dL1=H1-H2

(27)

(28)

其中θ1、φ1分別為在起始時(shí)刻監(jiān)測(cè)點(diǎn)的經(jīng)緯度；θ2、φ2為在觀測(cè)時(shí)刻監(jiān)測(cè)點(diǎn)的經(jīng)緯度；R為地球半徑，取平均值6371km。

已知目前RTK接收機(jī)動(dòng)態(tài)定位精度可達(dá)cm級(jí)，平面精度最高可達(dá)1cm，高程精度可達(dá)2cm。觀測(cè)時(shí)間足夠長(zhǎng)的情況下，RTK接收機(jī)的靜態(tài)定位精度可達(dá)mm級(jí)，平均誤差在1～2mm?？紤]到山體滑坡在平穩(wěn)形變時(shí)期形變速度極為緩慢，因此對(duì)于RTK接收機(jī)而言可以視為靜態(tài)定位，形變監(jiān)測(cè)精度可達(dá)mm級(jí)，使mm級(jí)的數(shù)據(jù)融合具有可行性。

4.3.2 土工帶傳感器數(shù)據(jù)

基于土工帶電阻值計(jì)算其形變量dL的計(jì)算公式如下

(29)

(30)

式中，R(mΩ)是土工帶在觀測(cè)時(shí)刻的電阻測(cè)量值；R0(mΩ)是土工帶在布置完成時(shí)的初始電阻值；L0(mm)是土工帶在布置完成時(shí)的初始長(zhǎng)度；k是比例系數(shù)；T(℃)為土壤溫度。

因?yàn)橥凉挪季W(wǎng)絡(luò)的相對(duì)尺寸較小，所以忽略山體凹凸性，在對(duì)一套監(jiān)測(cè)設(shè)備進(jìn)行分析時(shí)，可以看作將土工帶排布網(wǎng)絡(luò)安置在監(jiān)測(cè)點(diǎn)處的切平面上。土工帶傳感器排布方式如圖5所示。10條土工帶分成縱橫兩組，各有5條土工帶。每一條土工帶上有10個(gè)土工帶傳感器。10條土工帶按照網(wǎng)格狀排列，空心點(diǎn)是土工帶傳感器，實(shí)心點(diǎn)是2個(gè)土工帶傳感器重合點(diǎn)，RTK接收機(jī)安置在網(wǎng)絡(luò)的中心位置。

圖5 土工帶傳感器與RTK接收機(jī)的排布網(wǎng)絡(luò)

在山體上布設(shè)土工帶排布網(wǎng)絡(luò)時(shí)，使方形網(wǎng)絡(luò)的一邊沿山體等高線，該方向定義為橫向；2個(gè)臨邊垂直于山體等高線，該方向定義為縱向。對(duì)土工帶數(shù)據(jù)進(jìn)行中值濾波和平滑處理后得到監(jiān)測(cè)點(diǎn)的橫向位移x和縱向位移y，則由土工帶傳感器測(cè)得的監(jiān)測(cè)點(diǎn)的形變量dLtu為

(31)

土工帶傳感器測(cè)試設(shè)備如圖6所示，測(cè)試結(jié)果如表3所示。測(cè)試結(jié)果顯示，土工帶傳感器的形變監(jiān)測(cè)數(shù)據(jù)可以達(dá)到mm級(jí)的監(jiān)測(cè)精度，為高精度地質(zhì)形變監(jiān)測(cè)提供了可能，使mm級(jí)的數(shù)據(jù)融合具有可行性。

圖6 土工帶傳感器地質(zhì)采集設(shè)備

表3 土工帶傳感器設(shè)備測(cè)試

4.3.3 無(wú)人機(jī)攝影測(cè)量數(shù)據(jù)

通過(guò)無(wú)人機(jī)近景攝影測(cè)量技術(shù)[24-25]獲得監(jiān)測(cè)區(qū)域的三維點(diǎn)云數(shù)據(jù)和傾斜攝影處理結(jié)果。使用ArcMap對(duì)多期測(cè)量數(shù)據(jù)中的影像控制點(diǎn)進(jìn)行匹配后，在某一期的三維模型中加入另一期的相對(duì)三維模型，計(jì)算出前后兩期攝影測(cè)量中滑坡監(jiān)測(cè)點(diǎn)處的形變量dLUAV，并獲得形變區(qū)域的位置和大小。目標(biāo)區(qū)域的數(shù)字高程地圖(Digital Elevation Map，DEM)和數(shù)字水平地圖分別如圖7和圖8所示。無(wú)人機(jī)攝影測(cè)量技術(shù)的引入，在提供形變監(jiān)測(cè)信息的同時(shí)，還提供了地形地貌、建筑分布、交通情況等諸多信息要素，直觀地提供滑坡監(jiān)測(cè)點(diǎn)的位置信息，有利于災(zāi)害防治工作。

圖7 數(shù)字高程地圖

圖8 數(shù)字水平地圖

以大疆精靈Phantom無(wú)人機(jī)攝影測(cè)量為例，像元大小為2.41μm，相機(jī)焦距為8.8mm，飛行高度約為36.5m時(shí)，可以滿足客戶1cm的監(jiān)測(cè)精度要求，如圖9所示。考慮到在平穩(wěn)形變階段TUAV≈30d≈15T，所以在進(jìn)行數(shù)據(jù)融合計(jì)算時(shí),由無(wú)人機(jī)攝影測(cè)量得到的形變信息dLUAV的精度在1mm內(nèi)，在尺度與分辨率方面與土工帶傳感器形變量dLtu和RTK測(cè)量形變量dLRTKmm級(jí)精度保持一致。

圖9 無(wú)人機(jī)攝影測(cè)量精度分析

5 結(jié)果分析

本文按照已知的數(shù)據(jù)內(nèi)容和格式，綜合考慮典型案例中山體滑坡發(fā)生前后的形變規(guī)律，對(duì)數(shù)據(jù)進(jìn)行模擬設(shè)置?；谀M設(shè)置的警報(bào)級(jí)別滑坡監(jiān)測(cè)數(shù)據(jù)對(duì)本文的數(shù)據(jù)處理算法進(jìn)行仿真分析。最后通過(guò)對(duì)比模擬設(shè)置數(shù)據(jù)與算法處理結(jié)果，檢驗(yàn)山體滑坡形變監(jiān)測(cè)數(shù)據(jù)處理算法的功能。詳細(xì)分析內(nèi)容如下，主要包括警報(bào)級(jí)別數(shù)據(jù)模擬設(shè)置、抗差自適應(yīng)Kalman濾波融合結(jié)果分析、基于灰色預(yù)測(cè)理論和蠕變切線角判據(jù)的預(yù)測(cè)結(jié)果分析、抗差自適應(yīng)Kalman濾波算法對(duì)預(yù)測(cè)結(jié)果的影響分析4個(gè)方面的內(nèi)容。

5.1 實(shí)驗(yàn)數(shù)據(jù)模擬設(shè)置

5.1.1 Kalman濾波效果展示數(shù)據(jù)模擬設(shè)置

為直觀清晰地展示通過(guò)抗差自適應(yīng)Kalman濾波算法進(jìn)行數(shù)據(jù)融合的效果，本文模擬設(shè)置某期采集到的500條土工帶傳感器形變數(shù)據(jù)和500條RTK接收機(jī)形變數(shù)據(jù)，模擬設(shè)置無(wú)人機(jī)攝影測(cè)量到的滑坡監(jiān)測(cè)點(diǎn)的加速度信息，將模擬數(shù)據(jù)代入抗差自適應(yīng)Kalman濾波算法中進(jìn)行數(shù)據(jù)融合。數(shù)據(jù)模擬設(shè)置如表4所示。

表4 抗差自適應(yīng)Kalman濾波數(shù)據(jù)模擬設(shè)置說(shuō)明

5.1.2 警報(bào)級(jí)別數(shù)據(jù)模擬設(shè)置

為了檢驗(yàn)滑坡監(jiān)測(cè)算法在預(yù)警等級(jí)為警報(bào)等級(jí)時(shí)的監(jiān)測(cè)效果，按照處理算法的結(jié)構(gòu)并綜合考慮典型案例中山體滑坡發(fā)生前后的形變規(guī)律模擬出12期的數(shù)據(jù)進(jìn)行仿真分析，前8期的數(shù)據(jù)為測(cè)量模擬數(shù)據(jù)，實(shí)現(xiàn)對(duì)第9、10、11、12期數(shù)據(jù)的預(yù)測(cè)和危險(xiǎn)等級(jí)預(yù)警功能。警報(bào)數(shù)據(jù)模擬設(shè)置說(shuō)明如表5所示，數(shù)據(jù)模擬結(jié)果如表6所示。

表6 警報(bào)級(jí)別數(shù)據(jù)模擬設(shè)置結(jié)果

5.2 抗差自適應(yīng)Kalman濾波融合結(jié)果

上述內(nèi)容已經(jīng)說(shuō)明每隔2d采集500條土工帶傳感器和RTK接收機(jī)的監(jiān)測(cè)數(shù)據(jù)。針對(duì)每一期的大量監(jiān)測(cè)數(shù)據(jù)，本文結(jié)合無(wú)人機(jī)近景攝影測(cè)量計(jì)算得到的監(jiān)測(cè)點(diǎn)加速度信息，使用抗差自適應(yīng)Kalman濾波進(jìn)行數(shù)據(jù)融合，最終得到監(jiān)測(cè)點(diǎn)處精確的形變量，為實(shí)現(xiàn)后續(xù)的山體滑坡預(yù)測(cè)預(yù)警功能做數(shù)據(jù)準(zhǔn)備。抗差自適應(yīng)Kalman濾波算法的數(shù)據(jù)融合效果主要通過(guò)某一期的數(shù)據(jù)融合結(jié)果進(jìn)行展示。

抗差自適應(yīng) Kalman濾波前后的數(shù)據(jù)對(duì)比如圖10所示。圖10中,藍(lán)色實(shí)線為經(jīng)抗差自適應(yīng)Kalman濾波之后的滑坡監(jiān)測(cè)點(diǎn)的形變量，紅色點(diǎn)代表RTK定位的形變監(jiān)測(cè)數(shù)據(jù)，綠色點(diǎn)代表土工帶傳感器的形變監(jiān)測(cè)數(shù)據(jù)，黑色實(shí)線表示模擬設(shè)置的形變真值。調(diào)整無(wú)人機(jī)攝影測(cè)量的形變加速度信息，使之與土工帶數(shù)據(jù)和RTK數(shù)據(jù)相匹配，可以看出濾波融合效果良好，在真值附近收斂，本次數(shù)據(jù)融合效果展示中無(wú)人機(jī)近景攝影測(cè)量的形變加速度a=0.0015(mm·d-2)。

圖10 抗差自適應(yīng)Kalman濾波前后數(shù)據(jù)對(duì)比

從圖10可以發(fā)現(xiàn)，濾波收斂并穩(wěn)定之后數(shù)據(jù)的最大值為5.984mm，最小值為5.570mm，與最初設(shè)置的真值5.8mm之間的誤差不超過(guò)0.23mm，滿足mm級(jí)監(jiān)測(cè)精度要求。需要說(shuō)明的是，由于所選取的用于展示抗差自適應(yīng)Kalman濾波效果的這一期數(shù)據(jù)具有普適性，其監(jiān)測(cè)精度即可反映普遍數(shù)據(jù)經(jīng)抗差自適應(yīng)Kalman濾波之后的監(jiān)測(cè)精度，滿足本文數(shù)據(jù)處理精度要求，為之后的山體滑坡預(yù)測(cè)預(yù)警算法提供了精確的數(shù)據(jù)支撐。

5.3 基于灰色預(yù)測(cè)理論和蠕變切線角判據(jù)的預(yù)測(cè)結(jié)果

灰色預(yù)測(cè)理論形變預(yù)測(cè)結(jié)果與原模擬數(shù)據(jù)對(duì)比如圖11所示。對(duì)比結(jié)果顯示警戒級(jí)別數(shù)據(jù)經(jīng)過(guò)Kalman濾波融合并經(jīng)過(guò)灰色預(yù)測(cè)理論得到的警戒數(shù)據(jù)預(yù)測(cè)形變量曲線與警戒級(jí)別數(shù)據(jù)的形變量曲線基本重合，說(shuō)明灰色預(yù)測(cè)理論預(yù)測(cè)的形變量與模擬設(shè)置數(shù)據(jù)的形變量接近，預(yù)測(cè)效果良好。

圖11 灰色預(yù)測(cè)理論的形變預(yù)測(cè)效果圖

灰色預(yù)測(cè)理論預(yù)測(cè)等級(jí)與原數(shù)據(jù)等級(jí)對(duì)比如表7所示。對(duì)比結(jié)果顯示由灰色預(yù)測(cè)理論形變預(yù)測(cè)結(jié)果得到的預(yù)警等級(jí)與模擬設(shè)置數(shù)據(jù)的預(yù)警等級(jí)基本符合，僅在第4期數(shù)據(jù)上存在些許偏差，滑坡監(jiān)測(cè)算法具備有效的滑坡預(yù)測(cè)預(yù)警功能。

表7 灰色預(yù)測(cè)理論預(yù)測(cè)等級(jí)與原數(shù)據(jù)等級(jí)對(duì)比

警報(bào)級(jí)別數(shù)據(jù)實(shí)驗(yàn)結(jié)果如表8所示，均值方差比C1為0.0296，明顯小于0.35，即均值方差比精度等級(jí)為優(yōu)；小概率誤差P1為1，明顯大于0.95，即小概率誤差精度等級(jí)為優(yōu)，所以該灰色預(yù)測(cè)模型的精度等級(jí)為優(yōu)，形變預(yù)測(cè)結(jié)果精確。經(jīng)過(guò)抗差自適應(yīng)Kalman濾波的監(jiān)測(cè)點(diǎn)的精確形變數(shù)據(jù)通過(guò)灰色預(yù)測(cè)理論預(yù)測(cè)之后得到的多期數(shù)據(jù)dL_Kalman_Grey_predict呈現(xiàn)上升趨勢(shì)，即每隔2d滑坡監(jiān)測(cè)點(diǎn)的形變量越來(lái)越大，預(yù)計(jì)第12期(8d后)采集到的形變量能達(dá)到0.43714m；形變速度呈現(xiàn)明顯上升趨勢(shì)，即滑坡監(jiān)測(cè)點(diǎn)做加速形變，預(yù)計(jì)第12期(8天后)的形變速率能達(dá)到0.22388m/d；蠕變切線角angle呈現(xiàn)上升趨勢(shì)，預(yù)計(jì)第12期(8天后)能夠達(dá)到88.53°；預(yù)警等級(jí)warning_results也呈現(xiàn)上升趨勢(shì)，第8期(最近數(shù)據(jù)采集時(shí)刻)的預(yù)警等級(jí)為3級(jí)，預(yù)計(jì)第9期的預(yù)警等級(jí)為3級(jí)，第10、11、12期(未來(lái)8天)的預(yù)警等級(jí)為4級(jí)，即為警報(bào)等級(jí)，未來(lái)一周內(nèi)有非常大的概率出現(xiàn)山體滑坡災(zāi)害。

表8 警報(bào)級(jí)別數(shù)據(jù)實(shí)驗(yàn)結(jié)果

5.4 抗差自適應(yīng)Kalman濾波算法對(duì)預(yù)測(cè)預(yù)警的影響

為進(jìn)一步說(shuō)明抗差自適應(yīng)Kalman濾波算法對(duì)預(yù)測(cè)結(jié)果的影響，本文對(duì)未經(jīng)過(guò)抗差自適應(yīng)Kalman濾波融合的RTK形變監(jiān)測(cè)數(shù)據(jù)單獨(dú)通過(guò)灰色預(yù)測(cè)理論模型進(jìn)行形變預(yù)測(cè)和預(yù)警等級(jí)分析，將預(yù)測(cè)預(yù)警結(jié)果與模擬設(shè)置數(shù)據(jù)的形變量和預(yù)警結(jié)果進(jìn)行對(duì)比，同時(shí)與前文所述的經(jīng)過(guò)抗差自適應(yīng)Kalman濾波融合之后的滑坡預(yù)測(cè)預(yù)警效果進(jìn)行對(duì)比分析。如圖12所示，未經(jīng)過(guò)抗差自適應(yīng)Kalman濾波融合的RTK形變監(jiān)測(cè)數(shù)據(jù)單獨(dú)通過(guò)灰色預(yù)測(cè)理論進(jìn)行形變預(yù)測(cè)得到的結(jié)果與模擬設(shè)置數(shù)據(jù)的偏差較大，形變預(yù)測(cè)效果變差。如表9所示，未經(jīng)過(guò)抗差自適應(yīng)Kalman濾波融合的RTK形變監(jiān)測(cè)數(shù)據(jù)單獨(dú)通過(guò)灰色預(yù)測(cè)理論和蠕變切線角判據(jù)得到的預(yù)警等級(jí)與模擬設(shè)置數(shù)據(jù)的預(yù)警等級(jí)存在的偏差明顯增大，預(yù)警錯(cuò)誤率增加，監(jiān)測(cè)算法的預(yù)警效果變差。由表7和表9可知，使用抗差自適應(yīng)Kalman濾波融合之后的數(shù)據(jù),進(jìn)行形變預(yù)測(cè)和預(yù)警分析的效果明顯優(yōu)于使用單一形變監(jiān)測(cè)數(shù)據(jù)進(jìn)行預(yù)測(cè)預(yù)警分析。可見(jiàn)抗差自適應(yīng)Kalman濾波融合得到精確的形變監(jiān)測(cè)數(shù)據(jù),對(duì)于形變預(yù)測(cè)和預(yù)警分析是非常重要的。

圖12 灰色預(yù)測(cè)理論的形變預(yù)測(cè)效果圖(未經(jīng)Kalman濾波融合)

表9 灰色預(yù)測(cè)理論預(yù)測(cè)等級(jí)與原數(shù)據(jù)等級(jí)對(duì)比(未經(jīng)Kalman濾波融合)

6 總結(jié)

由上述分析可知，本文設(shè)計(jì)了一種基于抗差自適應(yīng)Kalman濾波和灰色預(yù)測(cè)理論的山體滑坡監(jiān)測(cè)技術(shù)。該技術(shù)通過(guò)多種監(jiān)測(cè)方式對(duì)易發(fā)生滑坡的監(jiān)測(cè)區(qū)域進(jìn)行形變監(jiān)測(cè)，通過(guò)抗差自適應(yīng)Kalman濾波技術(shù)對(duì)包括RTK定位數(shù)據(jù)、無(wú)人機(jī)攝影測(cè)量數(shù)據(jù)、土工帶傳感器數(shù)據(jù)在內(nèi)的多源數(shù)據(jù)進(jìn)行數(shù)據(jù)融合，最后通過(guò)灰色預(yù)測(cè)理論實(shí)現(xiàn)對(duì)監(jiān)測(cè)區(qū)域滑坡形變程度和危險(xiǎn)等級(jí)的預(yù)測(cè)。

與以往研究工作相比，本文所設(shè)計(jì)的滑坡監(jiān)測(cè)技術(shù)的價(jià)值主要體現(xiàn)在以下幾個(gè)方面：

1)本文所設(shè)計(jì)的滑坡監(jiān)測(cè)技術(shù)具有很強(qiáng)的實(shí)時(shí)性。RTK定位技術(shù)、土工帶傳感器的使用，使得山體滑坡監(jiān)測(cè)功能受天氣狀況和植被覆蓋等因素的影響大幅減弱，能夠及時(shí)發(fā)現(xiàn)和準(zhǔn)確預(yù)測(cè)到可能發(fā)生的滑坡災(zāi)害，有效預(yù)防財(cái)產(chǎn)損失和人員傷亡。

2)本文所設(shè)計(jì)的滑坡監(jiān)測(cè)技術(shù)自動(dòng)化程度高。場(chǎng)外的設(shè)備一經(jīng)布設(shè)，不需人工定期到現(xiàn)場(chǎng)進(jìn)行勘探。場(chǎng)外的設(shè)備故障能夠及時(shí)反饋到云端監(jiān)測(cè)平臺(tái)，只需有針對(duì)性地進(jìn)行修復(fù)即可。與此同時(shí)，少量故障設(shè)備回傳的錯(cuò)誤信息也能通過(guò)算法濾波自行消化，對(duì)預(yù)測(cè)結(jié)果造成的影響很小。

3)本文所設(shè)計(jì)的滑坡監(jiān)測(cè)技術(shù)的監(jiān)測(cè)精度高。由上述仿真分析可知，由于土工帶傳感器形變監(jiān)測(cè)、RTK靜態(tài)形變監(jiān)測(cè)和無(wú)人機(jī)攝影測(cè)量形變監(jiān)測(cè)的精度均能達(dá)到mm級(jí)，外加抗差自適應(yīng)Kalman濾波融合處理，該技術(shù)對(duì)于監(jiān)測(cè)點(diǎn)處的形變監(jiān)測(cè)精度可達(dá)mm級(jí)。

4)本文所設(shè)計(jì)的滑坡監(jiān)測(cè)技術(shù)的設(shè)備成本低。土工帶傳感器輔助RTK形變監(jiān)測(cè)和無(wú)人機(jī)攝影測(cè)量形變監(jiān)測(cè)，大大降低監(jiān)測(cè)成本，為大范圍的山體滑坡監(jiān)測(cè)提供了可行性。

本文目前尚存在一些不足之處，主要體現(xiàn)在受疫情影響和時(shí)間限制等多方面因素制約，缺乏實(shí)測(cè)數(shù)據(jù)進(jìn)行實(shí)驗(yàn)驗(yàn)證，算法的分析仍停留在仿真層面，該問(wèn)題待采得實(shí)測(cè)數(shù)據(jù)后即可解決。此外，本文所設(shè)計(jì)的滑坡監(jiān)測(cè)技術(shù)尚未考慮地理水文、人類(lèi)活動(dòng)等相關(guān)因素，后續(xù)工作計(jì)劃引入隨機(jī)森林模型對(duì)影響滑坡的諸多因素進(jìn)行分析。