一種基于灰色預測理論和抗差自適應Kalman濾波的滑坡監測算法

楊 旭，楊 旭，李 佳，王建國

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

0 引言

山體滑坡災害具有隱蔽性強、突發性強、破壞力大的特點。傳統的山體滑坡監測方法主要有宏觀地質觀測法、簡易觀測法、設站觀測法、儀表觀測法等，它們普遍存在的問題是需要人工定期到現場進行數據的采集，缺少實時性，監測精度較低，并且對滑坡的預測預警能力差[1-3]。

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

實時動態(Real-Time Kinematic，RTK)技術可實現實時定位，并且在靜態情況下的定位精度可達mm級。與傳統GNSS形變監測技術相比，RTK定位技術具有十分顯著的優勢。RTK技術被廣泛應用于安全監測、土地勘探、城市布局、水利等領域[7-11]。在我國的金川二礦區地裂縫[9]、黃土坡臨江I號崩滑體[10]等地質災害的變形監測中，均采用了高精度的全球定位系統(Global Positioning System，GPS)靜態相對定位技術，監測精度可以達到mm級；并逐漸開始嘗試采用GPS定位技術監測滑坡體的垂直形變，發現在基線較短、觀測時間較長和觀測條件較好的情況下，其精度也可以達到mm級[11]。

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

影響山體滑坡的因素復雜多變，有已知因素和未知因素，具有典型的灰色特征。灰色預測理論在解決信息不完全的系統問題時具有十分顯著的優勢。對于具有這樣特征的系統，可使用灰色預測理論對其進行深入的分析[16-19]。李曉鴿采用灰色預測理論對魯地拉水電站對外交通公路滑坡監測數據進行了擬合和預測研究，得到了很好的預測效果[16]；王旭昭等使用改進的灰色預測理論對新灘滑坡和黃茨滑坡進行預測分析，取得了良好預測效果[17]；此外，K. Yin[18]、王東岳[19]等均使用灰色預測理論進行了山體滑坡預測，并取得良好預測效果。

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

1 抗差自適應Kalman濾波

1.1 抗差自適應Kalman濾波[12-15]

在目標系統中，設Xi為由一組隨機變量組成的向量，稱為該系統在i時刻的狀態向量；設Zi為由一組觀測變量組成的向量，稱為該系統在i時刻的觀測向量。當目標系統為離散時間系統時，系統的狀態向量隨時間的變化規律可以通過狀態方程描述，目標系統的狀態向量和觀測向量的關系可以通過觀測方程描述

(1)

其中，Xi-1、Xi分別是目標系統在i-1和i時刻的狀態向量；Zi是目標系統在i時刻的觀測向量；U稱為控制向量；W稱為系統動態噪聲向量；V稱為系統觀測噪聲向量；A稱為系統從i-1時刻到i時刻的狀態轉移矩陣；B稱為控制向量系統矩陣；H稱為觀測值系統矩陣。

抗差自適應Kalman濾波技術將抗差估計和自適應估計相結合，以減小觀測誤差和系統建模誤差對濾波結果的影響，其技術原理可以通過如下方程加以描述：

時間更新方程

(2)

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

(3)

狀態更新方程

(4)

(5)

(6)

1.2 觀測噪聲等價協方差矩陣的構造

典型的觀測噪聲等價協方差矩陣D的構造方法[14]有：Huber構造法、丹麥法、IGG法、抗差貝葉斯估計法等。本文選用IGG1構造法，具體內容如下

(7)

其中，一般取k0=2.5，k1=3.5。Vi為觀測殘差，表達式如下

(8)

1.3 自適應因子的構造

本文采用兩段函數模型[14]進行自適應因子的構造

(9)

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

2 灰色預測理論

2.1 灰色預測模型(GM模型)

2.1.1 基本GM(1,1)模型

設X(0)為含有n個元素的序列，其中x(0)(k)≥0,k=1,2,…,n。上角標處括號內的數字為序列模型的階數，第二個括號內的數字為序列中元素的序號。對序列X(0)進行一次累加生成得到它的1-AGO序列X(1)。對序列X(1)進行一次緊鄰均值生成得到序列Z(1)。稱方程x(0)(k)+ax(1)(k)=b為GM(1,1)模型。其中a為發展系數，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 精度檢驗

(16)

(17)

在實際應用中，根據實際情況設定閾值A0，當A

(18)

在實際應用中，根據實際情況設定閾值P0，當P>P0時，則認為該模型為小概率誤差合格模型。關于均值方差比閾值A0和小概率誤差閾值P0的設定準則與模型精度等級的關系如表1所示。

表1 模型精度等級參考[19]

3 蠕變切線角判據

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

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

本文基于上述分析建立蠕變切線角預警判據[20]。通過一定時間的觀測得到穩定變化時期監測點的形變速度v0，將穩定變化時期的山體形變運動看作是勻速運動，可知v0是定值，定義

(19)

其中，S(k)為從起始監測時刻到當前監測時刻監測點處的累計位移量；T(k)為具有時間量綱的縱坐標值。由T(k)和t(k)繪制T-t曲線，如圖2所示，通過T-t曲線可以得到切線角的表達式

圖2 經過處理后得到的T-t曲線

(20)

其中，αk為切線角；t(k)為第k個監測時刻；Δt為采樣時間間隔。

根據行業內的長期觀測研究得到山體滑坡過程的預警判據表，如表2所示。該預警判據表綜合考慮了山體滑坡監測點的形變速率V和切線角α包含的信息內容，對山體滑坡預警等級進行了5層劃分，分別代表不同的山體滑坡可能性和危險性，從而實現對山體滑坡的等級預警功能。

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

表2 山體滑坡預警判據表[20]

4 實驗測試

4.1 滑坡監測系統

本文構建以RTK為核心的空天地一體化立體監測網絡，實現對滑坡災害的監測評估。

滑坡監測系統如圖3所示，使用RTK進行地質沉降監測，在監測區域布置土工帶傳感器，采集監測區域的位移信息并上傳到云端以待處理，利用無人機攝影測量技術提供空域增強信息，輔助RTK進行滑坡形變監測。北斗短報文通信技術可以實現全天候、全覆蓋域、高可靠性的通信，其通信時延約為0.5s，點對點通信時延約為1～5s，同時該技術的設備成本低，適用于廣泛的環境監測活動，所以RTK形變監測信息通過北斗短報文進行傳遞。

圖3 山體滑坡監測系統

設計后處理軌跡平滑算法剔除異常測量信息，使用數據融合[23]技術對RTK定位數據、土工帶傳感器數據、無人機三維成像數據進行綜合分析后，在滑坡監測點處計算得到位移、速度、加速度測量信息。本文選用抗差自適應Kalman濾波技術實現數據的融合分析，獲得精確的形變監測結果。在獲得經過融合處理的形變監測數據后，通過灰色預測理論對未來時刻的形變數據進行合理預測，得到未來時刻監測點處的形變信息，最后結合蠕變切線角判據實現對山體滑坡災害的預測預警功能。

設備布設的具體情況如圖4所示。在重點監測區域進行監測設備布設時，為合理兼顧設備成本和監測效果，主要將監測設備布設在可能滑坡體的頂部和底部，尤其是一些裂縫處。考慮到皖南山區植被覆蓋程度大、連續性強，植被覆蓋斷裂痕跡處也需重點布設監測設備。RTK接收機安置在網絡的中心位置。該土工帶傳感器與RTK接收機的排布網絡與無人機攝影測量得到的形變監測數據相結合，三種形變監測數據在空間上具備了數據融合的條件。

圖4 監測區域數據采集設備布設

4.2 滑坡監測模型

對山體滑坡進行監測分析時，監測點的形變可以使用運動模型進行分析。運動模型的描述方程如下

(21)

(22)

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

x(i)=dLtu(i)

(23)

(24)

加速度a(i)通過定期無人機攝影測量形變量和飛行間隔時間計算得到

(25)

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

4.3 數據采集

4.3.1 RTK定位數據

RTK定位數據包含RTK接收機的經度、緯度、天線高度、接收機的工作模式和觀測衛星數量等信息。通過RTK定位得到監測點的經度、緯度和天線高度，則監測點位移量dLRTK與豎直位移量dL1和水平位移量dL2滿足如下關系

(26)

dL1=H1-H2

(27)

(28)

其中θ1、φ1分別為在起始時刻監測點的經緯度；θ2、φ2為在觀測時刻監測點的經緯度；R為地球半徑，取平均值6371km。

已知目前RTK接收機動態定位精度可達cm級，平面精度最高可達1cm，高程精度可達2cm。觀測時間足夠長的情況下，RTK接收機的靜態定位精度可達mm級，平均誤差在1～2mm。考慮到山體滑坡在平穩形變時期形變速度極為緩慢，因此對于RTK接收機而言可以視為靜態定位，形變監測精度可達mm級，使mm級的數據融合具有可行性。

4.3.2 土工帶傳感器數據

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

(29)

(30)

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

因為土工帶排布網絡的相對尺寸較小，所以忽略山體凹凸性，在對一套監測設備進行分析時，可以看作將土工帶排布網絡安置在監測點處的切平面上。土工帶傳感器排布方式如圖5所示。10條土工帶分成縱橫兩組，各有5條土工帶。每一條土工帶上有10個土工帶傳感器。10條土工帶按照網格狀排列，空心點是土工帶傳感器，實心點是2個土工帶傳感器重合點，RTK接收機安置在網絡的中心位置。

圖5 土工帶傳感器與RTK接收機的排布網絡

在山體上布設土工帶排布網絡時，使方形網絡的一邊沿山體等高線，該方向定義為橫向；2個臨邊垂直于山體等高線，該方向定義為縱向。對土工帶數據進行中值濾波和平滑處理后得到監測點的橫向位移x和縱向位移y，則由土工帶傳感器測得的監測點的形變量dLtu為

(31)

土工帶傳感器測試設備如圖6所示，測試結果如表3所示。測試結果顯示，土工帶傳感器的形變監測數據可以達到mm級的監測精度，為高精度地質形變監測提供了可能，使mm級的數據融合具有可行性。

圖6 土工帶傳感器地質采集設備

表3 土工帶傳感器設備測試

4.3.3 無人機攝影測量數據

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

圖7 數字高程地圖

圖8 數字水平地圖

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

圖9 無人機攝影測量精度分析

5 結果分析

本文按照已知的數據內容和格式，綜合考慮典型案例中山體滑坡發生前后的形變規律，對數據進行模擬設置。基于模擬設置的警報級別滑坡監測數據對本文的數據處理算法進行仿真分析。最后通過對比模擬設置數據與算法處理結果，檢驗山體滑坡形變監測數據處理算法的功能。詳細分析內容如下，主要包括警報級別數據模擬設置、抗差自適應Kalman濾波融合結果分析、基于灰色預測理論和蠕變切線角判據的預測結果分析、抗差自適應Kalman濾波算法對預測結果的影響分析4個方面的內容。

5.1 實驗數據模擬設置

5.1.1 Kalman濾波效果展示數據模擬設置

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

表4 抗差自適應Kalman濾波數據模擬設置說明

5.1.2 警報級別數據模擬設置

為了檢驗滑坡監測算法在預警等級為警報等級時的監測效果，按照處理算法的結構并綜合考慮典型案例中山體滑坡發生前后的形變規律模擬出12期的數據進行仿真分析，前8期的數據為測量模擬數據，實現對第9、10、11、12期數據的預測和危險等級預警功能。警報數據模擬設置說明如表5所示，數據模擬結果如表6所示。

表6 警報級別數據模擬設置結果

5.2 抗差自適應Kalman濾波融合結果

上述內容已經說明每隔2d采集500條土工帶傳感器和RTK接收機的監測數據。針對每一期的大量監測數據，本文結合無人機近景攝影測量計算得到的監測點加速度信息，使用抗差自適應Kalman濾波進行數據融合，最終得到監測點處精確的形變量，為實現后續的山體滑坡預測預警功能做數據準備。抗差自適應Kalman濾波算法的數據融合效果主要通過某一期的數據融合結果進行展示。

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

圖10 抗差自適應Kalman濾波前后數據對比

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

5.3 基于灰色預測理論和蠕變切線角判據的預測結果

灰色預測理論形變預測結果與原模擬數據對比如圖11所示。對比結果顯示警戒級別數據經過Kalman濾波融合并經過灰色預測理論得到的警戒數據預測形變量曲線與警戒級別數據的形變量曲線基本重合，說明灰色預測理論預測的形變量與模擬設置數據的形變量接近，預測效果良好。

圖11 灰色預測理論的形變預測效果圖

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

表7 灰色預測理論預測等級與原數據等級對比

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

表8 警報級別數據實驗結果

5.4 抗差自適應Kalman濾波算法對預測預警的影響

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

圖12 灰色預測理論的形變預測效果圖(未經Kalman濾波融合)

表9 灰色預測理論預測等級與原數據等級對比(未經Kalman濾波融合)

6 總結

由上述分析可知，本文設計了一種基于抗差自適應Kalman濾波和灰色預測理論的山體滑坡監測技術。該技術通過多種監測方式對易發生滑坡的監測區域進行形變監測，通過抗差自適應Kalman濾波技術對包括RTK定位數據、無人機攝影測量數據、土工帶傳感器數據在內的多源數據進行數據融合，最后通過灰色預測理論實現對監測區域滑坡形變程度和危險等級的預測。

與以往研究工作相比，本文所設計的滑坡監測技術的價值主要體現在以下幾個方面：

1)本文所設計的滑坡監測技術具有很強的實時性。RTK定位技術、土工帶傳感器的使用，使得山體滑坡監測功能受天氣狀況和植被覆蓋等因素的影響大幅減弱，能夠及時發現和準確預測到可能發生的滑坡災害，有效預防財產損失和人員傷亡。

2)本文所設計的滑坡監測技術自動化程度高。場外的設備一經布設，不需人工定期到現場進行勘探。場外的設備故障能夠及時反饋到云端監測平臺，只需有針對性地進行修復即可。與此同時，少量故障設備回傳的錯誤信息也能通過算法濾波自行消化，對預測結果造成的影響很小。

3)本文所設計的滑坡監測技術的監測精度高。由上述仿真分析可知，由于土工帶傳感器形變監測、RTK靜態形變監測和無人機攝影測量形變監測的精度均能達到mm級，外加抗差自適應Kalman濾波融合處理，該技術對于監測點處的形變監測精度可達mm級。

4)本文所設計的滑坡監測技術的設備成本低。土工帶傳感器輔助RTK形變監測和無人機攝影測量形變監測，大大降低監測成本，為大范圍的山體滑坡監測提供了可行性。

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