基于改進切線角實時算法的滑坡監測預警

賈 丁，陳小紅，周永貴

（1.航天宏圖信息技術股份有限公司，北京 100089；2.自然資源部地質災害智能監測與風險預警工程技術創新中心，北京 100081）

0 引言

我國地質災害頻發，滑坡作為我國最主要的地質災害類型之一，嚴重威脅人民生命財產安全。我國滑坡數量龐大、分布廣泛、發生頻率高、突發性強且危害嚴重，滑坡災害形勢復雜嚴峻［1-4］。監測預警是主動防范地質災害的重要手段，滑坡以巖土體變形失控后產生運動的方式造成破壞，變形作為滑坡過程中的顯著表征，是滑坡預警的主要依據和關鍵的監測參數，變形監測是分析滑坡危險程度和演變規律的重要環節［5-6］。全球衛星導航系統（Global Navigation Satellite System，GNSS）作為監測變形和位移的重要方法在滑坡變形監測預警中得以廣泛應用［7-10］。

GNSS 的滑坡監測結果主要通過位移—時間數據呈現，對該數據的分析處理是監測預警中的重要流程［11-14］。對于一個完整的滑坡生命周期而言，理想的滑坡監測預警曲線分為初始變形階段、等速變形階段和加速變形階段，如圖1（a）所示。切線角指位移—時間曲線各點斜率和時間軸的夾角，是變形速度的體現，位移—時間曲線切線角的變化監測已經成為《崩塌、滑坡、泥石流監測規范DZ/T0221-2006》中推薦方法的重要技術參數。

Fig.1 Ideal landslide displacement-time curve and real complicated situation圖1 理想的滑坡位移—時間曲線和實際復雜情況

許強等［15］指出傳統切線角定義和計算的不足，提出了改進切線角，并在已有的滑坡應用中進行了檢驗。王珣等［16］通過建立等速變形速率和臨滑切線角的聯系，并結合未破壞滑坡變形過程的最大切線角建立了預警判據，得到了較好的應用效果。劉富有等［17］將NGM（1，1，k，c）模型與改進切線角結合對危巖變形趨勢進行分析預測。李洋［18］建立了以臨界位移量和T-t 曲線切線角“雙指標”為閾值進行滑坡預測預報的方法。Tan 等［19］使用切線角計算的思路對微變形雷達監測數據進行處理，并與滑坡預測中的速度和位移雙閾值相結合，進行礦區邊坡監測預警。

然而，在前人研究中，對切線角計算過程和方法的研究較少。由于計算切線角的求導過程會放大曲線本身的噪聲和波動，計算結果往往呈現零散分布狀態，且常常有負值出現。同時，改進切線角計算需要在勻速變形階段對速度進行量綱處理，而勻速變形階段需要有較長時間的累積以進行趨勢判斷。但實際情況往往復雜多變，如圖1（b）所示，監測曲線中可能遇到多個勻速階段或者短暫加速及減速階段，有的勻速速率過大，已經不能用于切線角計算。現有的計算方法中，勻速階段難以判定，隨機性強。上述因素對切線角計算的準確性和時效性均有影響。本文在對位移—時間數據進行移動平均、數據重采樣的基礎上，采用拉格朗日插值微分法實現了計算結果降噪；同時，提出將變形階段的含量加以實時統計，以滿足切線角計算中勻速階段的速度更新，提升滑坡監測預警可靠性。

1 改進切線角計算原理

由于環境噪聲和各種不確定因素影響，GNSS 得到的原始位移—時間數據具有較大波動性，因此在計算切線角計算前需要進行平滑。本文采用指數移動平均和重采樣的方法對數據進行降噪處理。

1.1 指數移動平均

滑坡的GNSS 位移—時間曲線是典型的時間序列。移動平均法在時間序列分析中應用非常廣泛［20］，該方法可以對時間序列起到較好的平滑作用，黃智偉［11］驗證了該方法在滑坡監測數據中有具有較好的應用效果。本文采用指數移動平均方法對以小時為單位的位移—時間曲線進行移動平滑。

指數移動平均法是對過去所有數據加權平均，權由平滑因子α決定，而且隨著時間的推移加權指數遞減，故稱為指數滑動平均。其公式如式（1）所示。

其中，t為窗口大小，0 ＜α≤1 為平滑因子，α可根據窗口大小計算，如2/(t+1)，也可自定義。(1-α)i為呈指數增加的權重，離預測時刻越近權重越大。圖2 為不同大小的平滑窗口得到的平滑結果，可以看出，越大的平滑窗口平滑力度越大，對原始曲線的降噪作用也越好。但過大的平滑窗口也會使曲線喪失一些關鍵的波動特征，需要適當選取。需要注意的是，指數移動平均結果會損失前面一個平滑窗口數量的數據，在對目標時間段數據進行處理時，需要往前冗余選取一個時間窗口數量的數據。

Fig.2 Smoothing effect of exponential moving average on displacement-time curve圖2 指數移動平均對位移—時間曲線的平滑效果

1.2 數據重采樣

數據重采樣是信號處理領域的常用方法。在時間序列中，重采樣可以將時間序列從一個頻率轉化為另一個頻率進行處理。將高頻率數據轉化為低頻率數據為降采樣，將低頻率轉化為高頻率為升采樣。本文采用平均采樣方法將每日以小時為單位的數據降采樣為以天為單位的位移—時間數據。圖3 為重采樣前后的結果，重采樣可以減少對曲線局部變化的關注，突出曲線整體變化趨勢。

Fig.3 Results of data resampling by day圖3 按天進行數據重采樣結果

1.3 切線角計算

改進切線角反映的是每一變形階段的變形速度與勻速變形階段速度的相對大小。切線角的改進可以統一不同時間單位，如周、天、小時下的計算結果。其計算原理為：通過用位移除以勻速變形階段的速度v的做法將曲線S-t的縱坐標變換為與橫坐標相同的時間量綱，即定義如式（2）所示。

其中，ΔS(i)表示某一單位時間段（一般采用一個監測周期，如1 天、1 周等）內的位移變化量；v表示等速變形階段的位移速率；T(i)表示變換后與時間相同量綱的縱坐標值。這樣，T替代S后的曲線稱為T-t曲線。根據T-t曲線，可以得到改進的切線角αi的表達式，如式（3）所示。

其中，αi表示改進切線角；ti表示某一監測時刻；Δt表示對應ΔS的單位時間段；ΔT表示單位時間段內T(i)的變化量。顯然，根據上述定義有：當αi＜45°時，滑坡處于初始變形階段；當αi≈45°時，滑坡處于等速變形階段；當αi＞45°時，滑坡處于加速變形階段。

2 改進切線角算法優化及動態計算方法

2.1 拉格朗日插值微分算法

上述切線角計算算法中的重要步驟是斜率計算，現有的斜率計算公式如式（4）所示。

通常采用式（4）這種用差分來代替微分的方法，求導本身會放大曲線波動和誤差，因此該計算方式得到的切線角結果波動較大。拉格朗日插值微分是一種典型的離散點求導方法，能夠同時考慮被計算點左右兩側點的情況，提高計算結果的精度。為了降低切線角計算結果的波動性，采用拉格朗日插值微分五點法中的中心插值微分公式進行計算，如式（5）所示。

在實際應用中，對于被計算數據左側早期數據點而言，可以選擇其左側數據點進行補充，因而可以始終維持f′(x2)計算方法。如圖4 所示，當步長為2 時，計算x2位置點的導數需要用到x0，x1，x3，x44 個點的值。對于右側點而言，當其右側沒有足夠的點維持f′(x2)的計算方法時，經過驗證可以通過其與往前2h點的差分作為其求導結果。本文結合前人工作，將式（5）修改為式（6）。

Fig.4 Schematic diagram of Lagrange interpolation differential calculation（h=2）圖4 拉格朗日插值微分計算示意圖（h=2）

其中，h為拉格朗日插值微分的計算步長，n為數據點總量。

2.2 階段判識與階段變形平均速度計算

如上所述，通過變形速度對位移—時間曲線的變形階段劃分為減速變形階段、勻速變形階段和加速變形階段。由于切線角計算要使用勻速階段速度，因此變形階段劃分判識至關重要。對于較長時間尺度和時間序列的數據而言，曲線的凸凹性可以很好地用來判斷曲線速度整體變化趨勢。王一帆等［21］提出了通過位移—時間曲線的凸凹性進行滑坡變形階段判識的方法：定義時間從t1～t3的數據對應的總位移為s1～s3，選取t2為t1～t3時刻的中點，其對應的位移為s2，定義R1=0.5(s3-s1)，R2=s2-s1。定義無量綱參數γ=R2/R1，當0.9 ≤γ≤1.1 時，認為曲線處于勻速變形階段；當γ＞1.1 時認為曲線處于減速變形階段；當γ＜0.9 時，認為曲線處于加速變形階段（見圖5［21］）。取t1～t2和t2～t3兩段速度平均值作為該階段的速度，即v=。

Fig.5 Deformation stage division according to displacement-time curve of convexity and concavity圖5 根據凸凹性的位移—時間曲線變形階段劃分

2.3 切線角計算勻速速度更新算法

通常，切線角計算方法中勻速階段需要長期積累，且計算結果隨機性強，對于長期運行、實時監測的監測預警系統而言，其時效性很難達到實用要求。本研究在變形階段劃分基礎上，通過一種變形階段實時統計方法對勻速階段的選取進行實時判斷。

對于某一天ti位移數據，選取其之前m天數據，即ti-m～ti位移—時間數據，通過上述方法可以判定得到該范圍位移—時間曲線的變形階段類型ci。向前統計窗口大小為d，即ci-d～ci的所有階段類型的含量，記減速、勻速、加速變形階段的含量為p1、p2、p3，同時計算ci-d～ci每一階段切線角大小的平均速度。

在上述3 種階段類型含量和切線角平均值大小的基礎上，對用于切線角計算的勻速速率是否更新作為條件進行判斷。這種判斷可以排除曲線進入加速階段、減速階段以及新的勻速階段變形速率過大等情況。

圖6 展示了當m=30、d=30 和d=10 時典型滑坡位移—時間曲線p1、p2、p3的變化趨勢。除完全進入加速階段后加速階段百分含量p3為1 外，其余時刻由于曲線的波動，3 種變形階段都會出現。即當曲線整體處于勻速變形階段時，局部的階段判識由于曲線波動影響可能會得到減速階段或者加速階段的結果。但整體趨勢為減速階段逐漸減少，在進入加速階段后逐漸上升，進入預警階段時為加速階段。

Fig.6 Description of stage content percentage and uniform speed algorithm parameters圖6 階段含量百分比統計和勻速速度算法參數

因此，可以通過限制減速變形階段百分含量（p1）和加速變形階段百分含量（p3）的上界、勻速變形階段百分含量（p2）的下界作為范圍條件。同時，為了避免在即將進入加速階段時勻速速度過大，即如果具有過大的切線角平均值，即使判定為勻速階段也不對勻速速度進行更新。設置該時段切線角平均值范圍大小。當p1、p2、p3和滿足上述條件時，對切線角計算所使用的勻速速度進行更新，以此為條件便可以形成勻速速度更新判斷方法，在此基礎上對時間序列進行迭代即可實現切線角計算結果的動態計算和更新。

在Ubuntu20.04 LTS 系統，Python3.10 環境下完成了優化后的切線角計算和勻速速度更新的算法編程實現。算法如下：

3 實驗結果與分析

本文對河南省鶴壁市滑坡災害監測數據進行實驗分析，選取位于靈山地質災害隱患點的監測數據，時間從2022 年7 月1 日至2023 年6 月30 日。在為期一年的GNSS位移—時間曲線數據中，選取其中一個月的曲線數據進行實驗分析，該時間階段屬于滑坡勻速變形階段。

3.1 計算參數對切線角結果的影響

在拉格朗日插值微分中，步長h插值微分的重要參數對計算結果有很大影響。同時，為了驗證指數移動平均法和數據重采樣在該算法中所起的作用，采用拉格朗日插值微分方法進行計算，比較不同計算步長h條件下和是否采取數據平滑與重采樣的條件下勻速變形階段的切線角大小。不同計算步長切線角計算結果如圖7 所示，移動平均對切線角計算的影響如圖8所示。

Fig.7 Calculation results of tangent angles with different calculation steps圖7 不同計算步長切線角計算結果

Fig.8 Influence of moving average method on calculation of tangent angle圖8 移動平均對切線角計算的影響

重采樣前后計算結果比較如圖9 所示。以小時為單位時間尺度進行比較，以72 個點的間距為步長進行計算，按每天的平均值進行重采樣后，72 小時的步長變為3 天。結果表明，重采樣的計算結果能夠包含更多數據信息，得到更為精確的結果。

Fig.9 Influence of data resample on calculation of tangent angle圖9 數據重采樣對切線角計算結果的影響

3.2 拉格朗日插值微分評價

對經過指數移動平均和重采樣的數據進行計算比較，圖10（a）和圖10（b）分別為采用式（4）和式（6）對一段勻速變形階段的位移—時間曲線的切線角計算結果，步長均設置為5。根據前述切線角計算原理，其理論切線角應為45°左右。式（4）的計算結果波動性大，比較分散，且有負值。而式（6）的計算結果則相對較集中，沿著45°分布。由此可見，該計算方法對切線角計算結果起到了很好的降噪作用，切線角計算結果波動性顯著下降。

Fig.10 Comparison of calculation results between general difference method and Lagrangian interpolation five-point method圖10 差分計算方法和拉格朗日插值五點法求導計算結果比較

為了量化改進后的效果，計算了兩組數據相較于勻速階段應有切線角45°的均方根誤差。其計算公式如式（7）所示。

其中，f(xi)為計算值或測量值；yi為真實值，在此處為45°。利用式（7）計算得到算法改進前切線角結果的均方根誤差為32.24°，改進后均方根誤差為7.31°，精度提高了77.33%。

為了驗證該算法在滑坡勻速和加速連續變形階段的適宜性，在已有一年期的GNSS 位移—時間曲線數據基礎上，根據李忠君等［22］提出的非線性蠕變模型，生成了滑坡勻速和加速連續變形的測試數據。依據上述切線角算法，對連續變形階段進行計算，結果顯示在勻速階段切線角值在45°附近上下波動，進入加速階段后波動減小并迅速增大至預警值，與實際情況較為吻合，如圖11所示。

Fig.11 Calculation result of tangent angle of landslide simulation curve圖11 滑坡曲線的切線角計算結果

4 結語

改進切線角作為滑坡監測預警中的重要判據指標，在各種監測預警模型中得以廣泛應用。本文對改進的切線角計算方法進行了優化研究，采用格朗日插值微分求導算法進行改進切線角的計算，該方法可以盡可能多地利用監測數據，從而提高了計算準確度。實驗表明，改進的拉格朗日插值微分算法相對于傳統的差分方法計算準確度提高了77.33%。在此基礎上分析了計算步長、指數移動平均法、數據重采樣對切線角計算結果的影響，最后通過加入階段判識和用勻速速度更新算法，實現了切線角動態計算。結果表明，指數移動平均法和數據重采樣均可以濾除位移—監測曲線數據中的原始噪聲，減少計算結果的隨機性。通過滑坡階段判識和階段含量統計，輔以切線角勻速速度更新的條件判斷，實現了切線角動態實時計算，為滑坡監測預警提供實時、動態判據，提高了滑坡監測預警實用性。

改進的算法需在后續滑坡監測預警中針對加速變形階段進一步加以驗證，優化切線角計算中的步長等參數，以適應復雜地質條件的滑坡監測預警。