孟永東,袁昌緯,田 斌,蔡征龍,張偉杰

(1. 三峽大學(xué)湖北長江三峽滑坡國家野外科學(xué)觀測研究站,湖北 宜昌 443002; 2. 三峽大學(xué)水利與環(huán)境學(xué)院,湖北 宜昌 443002)

滑坡作為一種多發(fā)性的地質(zhì)災(zāi)害,時刻威脅著人民的生命和財產(chǎn)安全。對滑坡易發(fā)區(qū)域進(jìn)行動態(tài)實時監(jiān)測,是及時發(fā)現(xiàn)并預(yù)報滑坡災(zāi)害發(fā)生、避免人員生命和財產(chǎn)安全受到嚴(yán)重侵害的重要手段之一。隨著無人機(jī)傾斜攝影技術(shù)的不斷進(jìn)步,以運動恢復(fù)結(jié)構(gòu)算法(structure from motion, SfM)為代表的一系列三維重建技術(shù)逐漸成熟,使以小型無人機(jī)為監(jiān)測工具的地表監(jiān)測系統(tǒng)開始承擔(dān)各種復(fù)雜地形條件下的滑坡監(jiān)測任務(wù)[1-2]。

目前,對滑坡表面位移的監(jiān)測主要有單點式和整體式兩種方法。其中,滑坡表面位移單點式監(jiān)測方法對滑坡的固定監(jiān)測點進(jìn)行位移監(jiān)測,進(jìn)而達(dá)到對滑坡實施動態(tài)監(jiān)測的目的[3-5],該方法雖然可對滑坡實施高精度動態(tài)監(jiān)測,但無法做到對滑坡的變動區(qū)域進(jìn)行整體性把握。因此,國內(nèi)外學(xué)者開展了關(guān)于滑坡表面位移整體式監(jiān)測的調(diào)查研究[6]。目前常用的方法主要有基于DEM的差分比較法(DEM of difference,DoD)、兩云直接比較法(direct cloud-to-cloud comparison with closest point technique,C2C)、云到網(wǎng)格和云到模型距離比較法(cloud-to-mesh distance or cloud-to-model distance,C2M)。DoD方法是在像素的基礎(chǔ)上將兩期DEM進(jìn)行比對區(qū)分,從而展示坡體的位移變化[7];但該方法無法用于三維模型的比較,對于復(fù)雜地形下的變化展示存在一定不足[8-9]。C2C方法具有簡單快捷的特點,無需對數(shù)據(jù)進(jìn)行匹配和網(wǎng)格化,對表面法線依賴度低,可用于三維邊坡的位移監(jiān)測;但該方法測得的位移也極易受點云精細(xì)度、誤差點及點間距的影響,常用于三維密集點云的快速變化檢測,不能實現(xiàn)對點云變化的精細(xì)測量。C2M方法通過點云并參考三維網(wǎng)格或理論模型之間的距離計算表面位移[10];但該方法只在地形平坦條件下表現(xiàn)良好,且需要人工創(chuàng)建網(wǎng)格,對丟失的數(shù)據(jù)進(jìn)行插值也會產(chǎn)生監(jiān)測位移的誤差[11-12]。鑒于此,文獻(xiàn)[13]提出多尺度模型對模型點云比對(multiscale model-to-model cloud comparison,M3C2)算法。相較于其他點云比對算法而言,該算法能直接在點云上檢測復(fù)雜地形變化,無需生成網(wǎng)格;且進(jìn)行變化計算時受空間點密度、表面粗糙度及不同采樣位置影響較小,目前已被成功應(yīng)用于湖濱暗溝侵蝕、支護(hù)結(jié)構(gòu)位移監(jiān)測及近海危巖監(jiān)測等領(lǐng)域[14-16]。

針對現(xiàn)有滑坡監(jiān)測方法所存在的缺陷,本文借助無人機(jī)傾斜攝影技術(shù),結(jié)合SfM算法與M3C2算法,實現(xiàn)復(fù)雜地形下的滑坡監(jiān)測需求。通過無人機(jī)采集的影像進(jìn)行運動恢復(fù)結(jié)構(gòu)分析,還原滑坡三維點云模型,并采用點云比對算法處理兩期航測點云數(shù)據(jù),以色值深淺表征滑坡區(qū)域點云變動距離,進(jìn)而識別坡體表面位移變化。

1 M3C2算法原理

為對比分析不同時期三維點云模型,采用M3C2點云比對算法,步驟包括選取核心點云、計算三維曲面法線、點云距離計算、確定空間變量置信區(qū)間等,如圖1所示。該算法可以直接在點云上檢測復(fù)雜地形的變化,無需網(wǎng)格劃分;且進(jìn)行變化計算時,受空間點密度、表面粗糙度及不同采樣位置的影響較小[17]。

由于選取合適的核心點云可顯著提高計算效率,因此計算初始階段需根據(jù)設(shè)定的距離將原始數(shù)據(jù)下采樣,得到密度較低且分布均勻的核心點云,并將其作為變化識別的基礎(chǔ)(如圖1(a)所示)。該方法能夠顯著提升計算效率,優(yōu)化時間復(fù)雜度,從而大幅縮短計算所需時間。

在選取合適的核心點云后,對于任何給定的核心點Pcore,在半徑為D/2的范圍內(nèi),均可與鄰域其他點云數(shù)據(jù)擬合出一個平面,并由此得出兩期點云的局部法向量N。如圖1(b)所示,記錄核心點Pcore半徑D/2范圍內(nèi)所有點到最佳擬合平面距離,并使用粗糙度σ(D) 表征標(biāo)準(zhǔn)差大小[18],即

(1)

沿法線方向從Pcore所在的擬合平面出發(fā),以d/2為投影半徑,存在一個通過Pcore、以法向量N為軸線且與兩期點云相交的圓柱體。搜索柱面內(nèi)包含的所有點云n1、n2,沿法向量分別計算兩期點云柱內(nèi)的平均位置,此時兩期點云在柱內(nèi)的平均位置分別為M1、M2,兩平均位置的差值即為間距LM3C2,也即坡體在Pcore點變化的距離。對整個坡體進(jìn)行迭代運算直至遍歷所有點云,可得出整個目標(biāo)區(qū)域的點云變化情況,即坡體表面位移,進(jìn)而識別滑坡等地質(zhì)災(zāi)害的發(fā)生,如圖1(c)所示。在計算過程中,能否設(shè)置合適的算法參數(shù),直接影響后續(xù)的計算效果。其中,投影半徑d/2、法向量半徑D/2,以及最大計算深度H為影響計算精度與效率的3個主要參數(shù)。

如圖1(d)所示,計算點云距離后,為估計局部距離變化量測精度,避免各種誤差導(dǎo)致坡體災(zāi)害識別誤判,需要進(jìn)一步確定空間置信區(qū)間,降低地質(zhì)災(zāi)害誤判的概率。在多次測量的誤差遵從獨立高斯分布的前提條件下,n1、n2≥ 30時采用z-雙尾差異顯著性檢驗公式計算置信水平0.95以上的置信區(qū)間,即

(2)

式中,n1、n2為d/2投影半徑下兩期點云的核心點點數(shù);REG代表兩期點云的配準(zhǔn)誤差;LOD95%(d)為投影半徑d/2下置信水平0.95以上置信區(qū)間的最小變化距離。其中,兩期點云的配準(zhǔn)誤差REG的計算公式為

(3)

式中,RMSE1為參照點云的均方根誤差;RMSE2為對比點云的均方根誤差。

當(dāng)4

(4)

式中,當(dāng)n1、n2小于4時,無需置信區(qū)間。

2 監(jiān)測試驗與數(shù)據(jù)處理

2.1 監(jiān)測試驗方案

為研究基于M3C2算法的無人機(jī)滑坡表面位移整體式監(jiān)測方法的可行性與識別精度,以及該算法相較于其他常用典型點云比對算法的優(yōu)劣性,開展室外試驗分析驗證。首先布設(shè)監(jiān)測點;其次規(guī)劃航線并采集影像;然后在此基礎(chǔ)上借助SfM算法還原目標(biāo)區(qū)域三維點云模型;最后用點云比對算法識別點云模型變化,與實際位移數(shù)據(jù)進(jìn)行對比分析,驗證該方法用于滑坡表面位移監(jiān)測的可行性,并探究其識別精度問題。

采用大疆Phantom 4 RTK無人機(jī)配套監(jiān)測系統(tǒng)采集影像數(shù)據(jù)(如圖2所示)。其機(jī)身搭載云臺相機(jī)型號為DJI-FC6310R,并具備三軸云臺增穩(wěn)系統(tǒng);采用1英寸大幅面CMOS傳感器,影像分辨率為5472×3648像素,搭配8.8 mm定焦鏡頭,能提供精準(zhǔn)的航測影像輸出;厘米級定位精度的GNSS定位系統(tǒng)搭配內(nèi)置RTK模塊,可精準(zhǔn)獲取每張圖像的POS數(shù)據(jù),即外方位元素X、Y、Z、ω、κ、γ,并保存在xmp字段中;操作端帶屏遙控器可規(guī)劃航測路徑,同時實時反饋無人機(jī)飛行狀況。

圖2 無人機(jī)監(jiān)測系統(tǒng)

2.2 野外試驗數(shù)據(jù)獲取

以三峽庫區(qū)某邊坡為試驗對象,在坡體表面布設(shè)4個監(jiān)測點,監(jiān)測點設(shè)計為柱狀結(jié)構(gòu),頂面、底面圓直徑均為115 cm。為精準(zhǔn)控制監(jiān)測點的移動距離,監(jiān)測點底部用木板保持水平并鋪設(shè)坐標(biāo)紙。4個監(jiān)測點位隨機(jī)分布,相互獨立,現(xiàn)場布置如圖3所示。

為保證無人機(jī)飛行安全,試驗前獲取該區(qū)域的DEM模型數(shù)據(jù),并設(shè)置仿地飛行模式,飛行航線以“井”字形規(guī)劃,航行高度為50 m。為保證建模的精細(xì)程度,飛行航向、旁向重疊率均設(shè)置為80%,相機(jī)傾斜角度為-60°。

在保證底板水平的前提下,同時水平移動4個監(jiān)測點,模擬滑坡發(fā)生時的坡表水平位移;或同時替換4個高度更低的監(jiān)測點,模擬滑坡發(fā)生時的坡表豎向位移。在移動或替換4個監(jiān)測點前后分別采用無人機(jī)監(jiān)測系統(tǒng)對目標(biāo)坡體進(jìn)行傾斜攝影,并獲取影像數(shù)據(jù)。試驗中實際控制的水平移動量和豎直變高都為10 mm。

2.3 建模流程及數(shù)據(jù)處理

建模及數(shù)據(jù)處理流程如圖4所示,大致可分為兩步:①三維點云模型獲取,基于尺度不變特征變換(scale-invariant feature transform,SIFT)算法、SfM算法及空中三角測量原理,對坡體進(jìn)行三維重建,得到坡體三維點云數(shù)據(jù),從而恢復(fù)坡體三維結(jié)構(gòu)。其主要包括圖像采集與處理、獲取內(nèi)外方位元素(POS數(shù)據(jù))、空中三角測量、三維點云模型構(gòu)建。②點云比對,將獲取的兩期三維點云模型進(jìn)行裁剪,采用M3C2算法比對兩期點云數(shù)據(jù),識別點云變動區(qū)域,進(jìn)而實現(xiàn)坡體表面位移的整體監(jiān)測。關(guān)鍵環(huán)節(jié)如下。

(1)多視圖影像采集環(huán)節(jié),采用傾斜攝影方式能有效解決航拍影像的變形及遮擋問題。但受相機(jī)鏡頭制造精度影響,拍攝的像片會出現(xiàn)不同程度的畸變,因此試驗關(guān)閉“畸變校正”后飛行,以獲取圖像畸變參數(shù)。

(2)參數(shù)獲取環(huán)節(jié),使用高斯函數(shù)將采集的影像序列卷積下采樣,以獲取高斯金字塔,并采用SIFT算法捕捉多視圖影像的特征點并進(jìn)行匹配。該過程提取的特征點又稱SIFT特征,具有高度穩(wěn)定性,不隨圖像拉伸、平移、明暗變化而變化,且能在一定程度上減小噪聲與視角變換帶來的負(fù)面影響,能有效改善無人機(jī)監(jiān)測系統(tǒng)在航測過程中因光照、坡度、轉(zhuǎn)向等各種因素變化而導(dǎo)致的圖像尺度差異、旋轉(zhuǎn)等問題。由當(dāng)前圖像與不同尺度核參數(shù)σ進(jìn)行卷積運算產(chǎn)生的一系列圖像即尺度空間圖像[19]。公式分別如下

L(x,y,σ)=I(x,y)*G(x,y,σ)

(5)

(6)

D(x,y,σ)=[G(x,y,kσ)-G(x,y,σ)]·I(x,y)

(7)

D(x,y,σ)=L(x,y,kσ)-L(x,y,σ)

(8)

式中,L(x,y,σ)為卷積函數(shù);(x,y)為特征點像素坐標(biāo);*為卷積運算符;G(x,y,σ)表示高斯函數(shù);σ為尺度空間參數(shù),決定圖像相應(yīng)的尺度及被平滑的程度,而尺度參數(shù)又能影響圖像的特征概貌與特征細(xì)節(jié)。

(3)增量式三維重建環(huán)節(jié),增量式三維重建方法通過迭代圖像配準(zhǔn)—三角化—采用光束法區(qū)域網(wǎng)平差解算地面坐標(biāo),并濾除誤差較大的點,遍歷所有像片后建立稀疏點云模型。在此基礎(chǔ)上,利用多視圖密集匹配(patch-based multi-view stereo,PMVS)算法,通過迭代延伸、擴(kuò)散稀疏點云構(gòu)建面元,最終形成密集點云模型。

(4)點云比對環(huán)節(jié),得到密集點云模型后,為節(jié)約算力資源,裁剪不同時期的點云,保留重點關(guān)注區(qū)域,對兩期點云進(jìn)行最近點配準(zhǔn)(iterative closest point,ICP)并迭代,將其轉(zhuǎn)換到同一坐標(biāo)系下并對齊。將對齊后的兩期點云圖采用M3C2算法進(jìn)行點云比對,以色值深淺變化表征位移距離大小。為了能夠在短時間內(nèi)對大型點云的兩次掃描之間的變化進(jìn)行估計,降低時間復(fù)雜度,M3C2算法使用特定的八叉樹(Octree)結(jié)構(gòu),以確保內(nèi)存使用和點鄰域提取速度之間存在良好比例。

2.4 坡體監(jiān)測試驗結(jié)果分析

為實現(xiàn)變化區(qū)域的可視化表達(dá),對邊坡點云數(shù)據(jù)進(jìn)行處理后,分別采用同為三維點云比對算法的C2C算法與M3C2算法對變動區(qū)域進(jìn)行識別,并比較識別效果。

2.4.1 C2C算法

在對坡體進(jìn)行三維重建的基礎(chǔ)上,將兩期點云裁剪、對齊,得到如圖5(a)、(d)所示的區(qū)域點云。采用 C2C 點云算法對位移前后的兩期點云數(shù)據(jù)進(jìn)行變化識別。設(shè)置合適參數(shù)后運行計算,得到兩期點云水平、豎直方向數(shù)據(jù)變化監(jiān)測結(jié)果(如圖5(b)、(e)所示),變化區(qū)域顯示設(shè)置為綠色。

圖5 C2C點云變化監(jiān)測(單位為m)

結(jié)果顯示,C2C算法能夠捕捉到監(jiān)測點水平方向和豎直方向上1 cm的位移變化。為進(jìn)一步探究其識別效果,控制變化顯示范圍為0.5～1.5 cm(如圖5(c)、(f)所示)。結(jié)果表明,雖然C2C三維點云比對算法能夠捕捉到監(jiān)測點在水平及豎直方向上厘米級的位移變化,但變化分布較為離散,并不能精準(zhǔn)提供其位移距離并采用色值大小體現(xiàn)。

2.4.2 M3C2算法

同理,對試驗區(qū)域完成坡體三維重建后,將兩期點云裁剪、對齊,得到如圖6(a)所示的區(qū)域點云。采用 M3C2 點云比對算法識別兩期點云數(shù)據(jù)變化。參數(shù)設(shè)置及步驟為:對第1期點云數(shù)據(jù)以 1.0 cm為間距進(jìn)行下采樣,得到稀疏、均質(zhì)的核心點云;設(shè)置步長為1.0 cm,從初始值12.8 cm開始,到終值51.6 cm結(jié)束,選取最合適的法向量擬合半徑,參考實際情況將最大計算深度設(shè)置為15 cm。兩期點云數(shù)據(jù)變化監(jiān)測結(jié)果如圖6(b)所示,可以發(fā)現(xiàn)位移變化集中在0～2 cm,顯示顏色為綠色;而灰色區(qū)域?qū)儆跓o變動區(qū)域。

圖6 M3C2點云變化監(jiān)測(單位為m)

為更加直觀地顯示監(jiān)測點的變動狀況,控制顯示范圍為0.5～1.5 cm。如圖6(c)所示,可以明顯觀測到0.5～1.5 cm的變動區(qū)域只存在于監(jiān)測點所處位置范圍內(nèi),即綠色區(qū)域?；诖?有理由認(rèn)為該方法可以準(zhǔn)確顯示水平方向的厘米級位移變化。同理,由圖6(d)—(f)可知,該方法也能夠反映垂直方向的厘米級位移變化,且變動區(qū)域顯示為粉色。雖然有極少部分噪聲存在,但總體而言顯示效果良好,能清晰地反映出變動區(qū)域和非變動區(qū)域,且能采用色值深淺體現(xiàn)其位移距離大小。

因此,基于M3C2算法的無人機(jī)滑坡表面位移整體式監(jiān)測方法可以準(zhǔn)確觀測并捕捉到厘米級的地物位移變化,相較于C2C點云比較算法而言,雖然花費時間更長,具有更高的時間復(fù)雜度,但噪點更少精度更優(yōu)。

2.5 滑坡現(xiàn)場試驗

野外試驗雖然證明了基于M3C2點云比對算法識別地物位移的有效性,但該方法對于滑坡表面位移監(jiān)測的有效性仍不夠明確,同時為比較M3C2算法和C2C算法分別運用于實地滑坡監(jiān)測的效果,以湖北省宜昌市龍泉鎮(zhèn)某地質(zhì)災(zāi)害點為現(xiàn)場試驗對象,對上述問題進(jìn)行驗證。

采用基于無人機(jī)航測點云比對的滑坡表面位移監(jiān)測方法對該區(qū)域坡體進(jìn)行多次監(jiān)測,分別于2021年12月18日、2022年1月8日,以及強(qiáng)降雨后的2022年3月29日,對目標(biāo)區(qū)域進(jìn)行圖像數(shù)據(jù)采集。綜合考慮坡高、樹木遮擋及飛行安全因素,3次試驗飛行均設(shè)置航高為60 m(無人機(jī)到坡底的距離);航線規(guī)劃的航向及旁向重疊率均設(shè)置為80%,相機(jī)傾斜角度為-60°。所獲取的變形區(qū)現(xiàn)場照片如圖7所示,獲得的監(jiān)測對象重點監(jiān)測區(qū)域的點云數(shù)據(jù)如圖8所示。

圖7 監(jiān)測區(qū)域現(xiàn)場照片

圖8 監(jiān)測區(qū)域點云數(shù)據(jù)

對3期點云進(jìn)行裁剪、配準(zhǔn),分別對3次采集的點云數(shù)據(jù)采用M3C2和C2C點云算法進(jìn)行量化識別,以判別兩種算法的應(yīng)用效果,如圖9所示。

圖9 變化區(qū)識別(單位為m)

M3C2算法以1.0 cm為間隔進(jìn)行下采樣;設(shè)置步長為1 cm,以初始值19 cm開始,到終值78 cm結(jié)束,計算M3C2距離,最大計算深度為5 m,以色值深淺表征位移大小。比對3期點云數(shù)據(jù)發(fā)現(xiàn),地質(zhì)災(zāi)害點坡體表面位移主要發(fā)生在3月強(qiáng)降雨后,比對結(jié)果如圖9(c)所示;圖中存在白色、藍(lán)色、紅色3種顏色區(qū)域。參考顏色標(biāo)尺,白色區(qū)域?qū)?yīng)的變動值為0,即該區(qū)域未發(fā)生位移變化,可認(rèn)為屬于坡體穩(wěn)定區(qū);相對白色而言的其他顏色區(qū)域表明有位移發(fā)生,可認(rèn)為屬于坡體不穩(wěn)定區(qū)域。觀察災(zāi)害識別圖可以看出,目標(biāo)區(qū)域自上而下經(jīng)歷了藍(lán)色、淺白色、紅色區(qū)域。其中,藍(lán)色區(qū)域為坍落區(qū),紅色區(qū)域為滑坡堆積區(qū),最大變化值大于1.5 m。基于此可以判斷,由于強(qiáng)降雨導(dǎo)致藍(lán)色區(qū)域坍落區(qū)土體發(fā)生滑動,滑動土體最后于坡腳處堆積,堆積區(qū)標(biāo)識為紅色。此外有其他極小范圍的稀疏紅色區(qū)域和藍(lán)色區(qū)域存在,造成該現(xiàn)象的原因為植被自然生長及人為搬運物體干擾造成的點云變化。

C2C算法的識別結(jié)果如圖9(d)所示?？梢钥闯?雖然C2C算法也能夠識別出變動區(qū)域,但相對于M3C2算法而言,無法通過顏色區(qū)分坍落區(qū)與堆積區(qū),同時具有更多的噪聲,識別效果較差。而M3C2點云比對算法的識別精度更高,顯示的變動區(qū)域更準(zhǔn)確,更適用于滑坡發(fā)育過程中的早期變形判識。

3 結(jié) 論

本文基于無人機(jī)航測點云比對的滑坡表面位移監(jiān)測方法,開展室外試驗,驗證了該方法的識別精度,并有效應(yīng)用于邊坡地質(zhì)災(zāi)害的實地調(diào)查與識別。結(jié)論如下:

(1)借助無人機(jī)傾斜攝影技術(shù),可以獲取地質(zhì)災(zāi)害隱患點的多視圖影像序列,借助SfM和MVS算法能夠獲取目標(biāo)區(qū)域高密度三維點云,并以此構(gòu)建精細(xì)化三維模型。

(2)采用M3C2多尺度模型與模型點云比對算法,可以捕捉到1 cm的水平、垂直位移變化,且以色值差異表征滑坡面域變形狀況,成功識別出地質(zhì)災(zāi)害隱患點位的滑坡變動區(qū)。相對于C2C點云比對算法而言,雖然時間復(fù)雜度更高,但噪點更少,精度更高,性能更優(yōu)。

(3)基于無人機(jī)航測點云比對的滑坡表面位移監(jiān)測方法,可以建立不同時期、不同狀態(tài)下的滑坡三維模型,不僅能整體性把握滑坡的位移變化過程,且能及時對滑坡災(zāi)害的發(fā)生作出預(yù)警,可實現(xiàn)復(fù)雜地形條件下滑坡的自動化監(jiān)測和地質(zhì)災(zāi)害識別。