基于改進(jìn)切線角實(shí)時(shí)算法的滑坡監(jiān)測(cè)預(yù)警

賈 丁，陳小紅，周永貴

（1.航天宏圖信息技術(shù)股份有限公司，北京 100089；2.自然資源部地質(zhì)災(zāi)害智能監(jiān)測(cè)與風(fēng)險(xiǎn)預(yù)警工程技術(shù)創(chuàng)新中心，北京 100081）

0 引言

我國(guó)地質(zhì)災(zāi)害頻發(fā)，滑坡作為我國(guó)最主要的地質(zhì)災(zāi)害類型之一，嚴(yán)重威脅人民生命財(cái)產(chǎn)安全。我國(guó)滑坡數(shù)量龐大、分布廣泛、發(fā)生頻率高、突發(fā)性強(qiáng)且危害嚴(yán)重，滑坡災(zāi)害形勢(shì)復(fù)雜嚴(yán)峻［1-4］。監(jiān)測(cè)預(yù)警是主動(dòng)防范地質(zhì)災(zāi)害的重要手段，滑坡以巖土體變形失控后產(chǎn)生運(yùn)動(dòng)的方式造成破壞，變形作為滑坡過(guò)程中的顯著表征，是滑坡預(yù)警的主要依據(jù)和關(guān)鍵的監(jiān)測(cè)參數(shù)，變形監(jiān)測(cè)是分析滑坡危險(xiǎn)程度和演變規(guī)律的重要環(huán)節(jié)［5-6］。全球衛(wèi)星導(dǎo)航系統(tǒng)（Global Navigation Satellite System，GNSS）作為監(jiān)測(cè)變形和位移的重要方法在滑坡變形監(jiān)測(cè)預(yù)警中得以廣泛應(yīng)用［7-10］。

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

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

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

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

1 改進(jìn)切線角計(jì)算原理

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

1.1 指數(shù)移動(dòng)平均

滑坡的GNSS 位移—時(shí)間曲線是典型的時(shí)間序列。移動(dòng)平均法在時(shí)間序列分析中應(yīng)用非常廣泛［20］，該方法可以對(duì)時(shí)間序列起到較好的平滑作用，黃智偉［11］驗(yàn)證了該方法在滑坡監(jiān)測(cè)數(shù)據(jù)中有具有較好的應(yīng)用效果。本文采用指數(shù)移動(dòng)平均方法對(duì)以小時(shí)為單位的位移—時(shí)間曲線進(jìn)行移動(dòng)平滑。

指數(shù)移動(dòng)平均法是對(duì)過(guò)去所有數(shù)據(jù)加權(quán)平均，權(quán)由平滑因子α決定，而且隨著時(shí)間的推移加權(quán)指數(shù)遞減，故稱為指數(shù)滑動(dòng)平均。其公式如式（1）所示。

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

Fig.2 Smoothing effect of exponential moving average on displacement-time curve圖2 指數(shù)移動(dòng)平均對(duì)位移—時(shí)間曲線的平滑效果

1.2 數(shù)據(jù)重采樣

數(shù)據(jù)重采樣是信號(hào)處理領(lǐng)域的常用方法。在時(shí)間序列中，重采樣可以將時(shí)間序列從一個(gè)頻率轉(zhuǎn)化為另一個(gè)頻率進(jìn)行處理。將高頻率數(shù)據(jù)轉(zhuǎn)化為低頻率數(shù)據(jù)為降采樣，將低頻率轉(zhuǎn)化為高頻率為升采樣。本文采用平均采樣方法將每日以小時(shí)為單位的數(shù)據(jù)降采樣為以天為單位的位移—時(shí)間數(shù)據(jù)。圖3 為重采樣前后的結(jié)果，重采樣可以減少對(duì)曲線局部變化的關(guān)注，突出曲線整體變化趨勢(shì)。

Fig.3 Results of data resampling by day圖3 按天進(jìn)行數(shù)據(jù)重采樣結(jié)果

1.3 切線角計(jì)算

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

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

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

2 改進(jìn)切線角算法優(yōu)化及動(dòng)態(tài)計(jì)算方法

2.1 拉格朗日插值微分算法

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

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

在實(shí)際應(yīng)用中，對(duì)于被計(jì)算數(shù)據(jù)左側(cè)早期數(shù)據(jù)點(diǎn)而言，可以選擇其左側(cè)數(shù)據(jù)點(diǎn)進(jìn)行補(bǔ)充，因而可以始終維持f′(x2)計(jì)算方法。如圖4 所示，當(dāng)步長(zhǎng)為2 時(shí)，計(jì)算x2位置點(diǎn)的導(dǎo)數(shù)需要用到x0，x1，x3，x44 個(gè)點(diǎn)的值。對(duì)于右側(cè)點(diǎn)而言，當(dāng)其右側(cè)沒(méi)有足夠的點(diǎn)維持f′(x2)的計(jì)算方法時(shí)，經(jīng)過(guò)驗(yàn)證可以通過(guò)其與往前2h點(diǎn)的差分作為其求導(dǎo)結(jié)果。本文結(jié)合前人工作，將式（5）修改為式（6）。

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

其中，h為拉格朗日插值微分的計(jì)算步長(zhǎng)，n為數(shù)據(jù)點(diǎn)總量。

2.2 階段判識(shí)與階段變形平均速度計(jì)算

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

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

2.3 切線角計(jì)算勻速速度更新算法

通常，切線角計(jì)算方法中勻速階段需要長(zhǎng)期積累，且計(jì)算結(jié)果隨機(jī)性強(qiáng)，對(duì)于長(zhǎng)期運(yùn)行、實(shí)時(shí)監(jiān)測(cè)的監(jiān)測(cè)預(yù)警系統(tǒng)而言，其時(shí)效性很難達(dá)到實(shí)用要求。本研究在變形階段劃分基礎(chǔ)上，通過(guò)一種變形階段實(shí)時(shí)統(tǒng)計(jì)方法對(duì)勻速階段的選取進(jìn)行實(shí)時(shí)判斷。

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

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

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

Fig.6 Description of stage content percentage and uniform speed algorithm parameters圖6 階段含量百分比統(tǒng)計(jì)和勻速速度算法參數(shù)

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

在Ubuntu20.04 LTS 系統(tǒng)，Python3.10 環(huán)境下完成了優(yōu)化后的切線角計(jì)算和勻速速度更新的算法編程實(shí)現(xiàn)。算法如下：

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

本文對(duì)河南省鶴壁市滑坡災(zāi)害監(jiān)測(cè)數(shù)據(jù)進(jìn)行實(shí)驗(yàn)分析，選取位于靈山地質(zhì)災(zāi)害隱患點(diǎn)的監(jiān)測(cè)數(shù)據(jù)，時(shí)間從2022 年7 月1 日至2023 年6 月30 日。在為期一年的GNSS位移—時(shí)間曲線數(shù)據(jù)中，選取其中一個(gè)月的曲線數(shù)據(jù)進(jìn)行實(shí)驗(yàn)分析，該時(shí)間階段屬于滑坡勻速變形階段。

3.1 計(jì)算參數(shù)對(duì)切線角結(jié)果的影響

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

Fig.7 Calculation results of tangent angles with different calculation steps圖7 不同計(jì)算步長(zhǎng)切線角計(jì)算結(jié)果

Fig.8 Influence of moving average method on calculation of tangent angle圖8 移動(dòng)平均對(duì)切線角計(jì)算的影響

重采樣前后計(jì)算結(jié)果比較如圖9 所示。以小時(shí)為單位時(shí)間尺度進(jìn)行比較，以72 個(gè)點(diǎn)的間距為步長(zhǎng)進(jìn)行計(jì)算，按每天的平均值進(jìn)行重采樣后，72 小時(shí)的步長(zhǎng)變?yōu)? 天。結(jié)果表明，重采樣的計(jì)算結(jié)果能夠包含更多數(shù)據(jù)信息，得到更為精確的結(jié)果。

Fig.9 Influence of data resample on calculation of tangent angle圖9 數(shù)據(jù)重采樣對(duì)切線角計(jì)算結(jié)果的影響

3.2 拉格朗日插值微分評(píng)價(jià)

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

Fig.10 Comparison of calculation results between general difference method and Lagrangian interpolation five-point method圖10 差分計(jì)算方法和拉格朗日插值五點(diǎn)法求導(dǎo)計(jì)算結(jié)果比較

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

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

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

Fig.11 Calculation result of tangent angle of landslide simulation curve圖11 滑坡曲線的切線角計(jì)算結(jié)果

4 結(jié)語(yǔ)

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

改進(jìn)的算法需在后續(xù)滑坡監(jiān)測(cè)預(yù)警中針對(duì)加速變形階段進(jìn)一步加以驗(yàn)證，優(yōu)化切線角計(jì)算中的步長(zhǎng)等參數(shù)，以適應(yīng)復(fù)雜地質(zhì)條件的滑坡監(jiān)測(cè)預(yù)警。