基于快速滑面重構的滑坡體積計算方法研究

摘 要：在滑坡等應急測繪中，快速計算滑坡體的體積對災害評估和應對、資源調配和救援指導有重要的意義。隨著無人機測繪技術的軟、硬件發(fā)展，現(xiàn)場數(shù)據(jù)的獲取變得容易、高效。本文旨在以無人機測繪實景三維數(shù)據(jù)為基礎，通過利用滑坡后壁暴露的部分幾何信息重構滑動面，結合Voronoi圖和Delaunay三角剖分進行插值計算，模擬下滑面的三維形狀，通過實景三維表面和下滑面數(shù)據(jù)計算滑坡體積。縮短了計算滑坡體體積的時間，提高了滑坡體積計算的準確性和可靠性，具有較好的實用價值和應用前景。

關鍵詞：應急測繪；滑坡；滑面重構

中圖分類號：P 642 文獻標志碼：A

滑坡是自然災害中最具破壞性的一種，由巖石、土壤或人造物質在重力作用下不穩(wěn)定運動引發(fā)，對人類生活和環(huán)境構成威脅，全球每年造成數(shù)千人死亡，數(shù)十億美元經(jīng)濟損失。中國2019年共發(fā)生地質災害6181次，滑坡占68.27%，2022年共發(fā)生滑坡等地質災害5659起，主要發(fā)生在中南、華南、西南，特別是四川、甘肅和江浙一帶。

無人機測繪技術的普及和發(fā)展，使獲取DEM、DOM和實景三維模型的技術手段更成熟、速度更快、精度更高，尤其是基于RTK技術在無人機設備上的應用，使無人機外業(yè)工作中獲取POS數(shù)據(jù)的精度已優(yōu)于厘米級，可以在重建工作中免相控進行測繪[1-3]，進一步提高了測繪效率。然而，滑坡體積分析至少需要兩期地表數(shù)據(jù)，即原地面和滑面數(shù)據(jù)，但滑坡前地表數(shù)據(jù)往往未被有效測量，衛(wèi)星遙感數(shù)據(jù)高程精度不高，在無勘察配合下難以快速獲取坡滑面數(shù)據(jù)。

本文通過無人機航測獲得地表實景三維數(shù)據(jù)，利用原地表和滑面連續(xù)、光滑的自然特性，快速模擬地表或滑面數(shù)據(jù)，構建兩期三角網(wǎng)，使用三角網(wǎng)法或方格網(wǎng)法快速計算滑坡體積。后期使用地質災害勘察成果中的鉆孔確定滑面特征點后，可以較準確地模擬滑面并準確計算滑坡體積。

1 滑面構建的原理

在滑面的構建中，獲得滑面的幾何特征點比較重要。然而，在許多情況下，無法獲取詳細的滑動面信息。一種可能的解決方案是利用滑坡后壁暴露的部分幾何信息來重構滑動面。這是因為滑動面作為一個連續(xù)的整體，其在滑坡后壁上的部分幾何信息可以為滑動面在地下的延續(xù)形態(tài)提供一定的指示。因此，如果能夠知道滑動面的幾何形態(tài)，并將后壁形態(tài)擬合到這個特定的幾何面上，那么得到的擬合結果可以作為獲取計算滑坡體積所需的地下信息。

本方法的數(shù)據(jù)主要利用無人機測繪實景三維的成果為基礎：滑坡的邊界一般可以從實景三維中判斷得到，能方便地讀取滑坡邊界的三維坐標。獲得滑面上少量的特征點即可使用邊界加特征點構建滑坡滑面。特征點的深度Z可以通過實景三維數(shù)據(jù)中滑面既有的形狀來推斷。

通過輸入的邊界點坐標集以及下滑區(qū)滑面特征點坐標，構建Voronoi圖初始的三角形網(wǎng)格并進行Delaunay三角剖分，刪除多余的三角形，在被剖分的三角形內根據(jù)等分或者等距內插點，計算點的三維坐標，形成模擬下滑面。

對離散分布的輸入特征數(shù)據(jù)來說，根據(jù)計算者給定的內插間距或等分數(shù)量進行剖分計算，本文構建線性grid函數(shù)進行插值計算。首先，對二維空間進行線性插值。將其中任意三角形頂點坐標分別隨機命名為（x1，y1）、（x2，y2，）、（x3，y3，），選擇任意1條邊作為基底，給定點的間距為d，將三角形的高除以d來確定在高方向上需要生成的點的數(shù)量，沿著高方向上點的坐標如公式（1）所示。

[x3+i×（xbase-x3）/n，y3+i×（ybase+x3）/n] （1）

式中：（xbase，ybase）為基底上的點，通過線性插值得到；n為需要等分的數(shù)量；i為正整數(shù)，i ∈ [0，n]。

沿著高方向上的每一個點，在其與基底的連線上等間距地生成點。假設需要在點[x3+i×（xbase-x3）/n，y3+i×（ybase+y3）/n]與基底的連線上生成m個點，坐標如公式（2）所示。

[x3+i×（xbase-x3）/n+j×（xbase-x1）/m，y3+i×（ybase-y3）/n+j×

（ybase-y1）/m] （2）

式中：m為在高方向上等分點的數(shù)量；j為正整數(shù)，j∈[0，m]。

其中，內插生成平面新點的高程值根據(jù)周圍臨近輸入的高程加入權重計算，權重的確定根據(jù)離內插點最近的特征點的距離確定。權重關系在后期數(shù)據(jù)驗證階段可以根據(jù)實際需要加入其他算法，供使用者選擇。如果點p處有k個臨近點，（i=1，2，3，…，n），那么點p的高程值如公式（3）所示。

（3）

式中：Zp為p點的高程值；Hi為p點臨近的第i個點的高程值；Li為p點臨近的第i個點距p點的距離。

如果有勘察成果，通過勘察鉆孔中對實際地層分布的揭露，就可以精確確定下滑面相對于地表的位置，通過從實景三維模型中獲得勘察點的高程值（即坐標Z值）減去勘察中確定的滑面深，得到較準確的滑面特征點三維坐標，勘察點數(shù)量越多越能準確地模擬下滑面的三維形狀。配合航測獲得的三維影像，獲取準確的地表特征，使用構建的滑面和地表特征并采用方格網(wǎng)法或者三角網(wǎng)法計算下滑體體積。

根據(jù)構建原理，使用Python編制了計算程序，程序界面如下，支持簡單的四點輸入坐標模擬插值和多點文件輸入插值；并可以通過等分和等距2種方式控制插值點間隔。通過繪圖窗口展示插值后點的分布情況（如圖1所示），如果出現(xiàn)與預期不同的點位分布，可以通過控制特征點的數(shù)量、分布、高程等調整輸入數(shù)據(jù)，重新進行插值計算。直到與專家系統(tǒng)對滑坡體下畫面預期一致，即可進行下滑量的計算。圖1中深色表示輸入的已知特征點，淺色表示由已知點通過程序構建的面特征點，點的密度通過點間距離或者等分間隔數(shù)來動態(tài)控制。

本文方法可以通過Python方便、高效地實現(xiàn)程序功能，通過numpy和pandas庫實現(xiàn)數(shù)據(jù)的輸入和輸出功能，并根據(jù)需要調整輸入、輸出數(shù)據(jù)的格式。使用數(shù)學工具包SciPy，SciPy包括的模塊有最優(yōu)化、線性代數(shù)、積分、插值、特殊函數(shù)、快速傅里葉變換、信號處理和圖像處理、常微分方程求解等科學計算功能。

為了提高程序的可用性，設計程序使用主界面，如圖2所示。SciPy中spatial模塊提供了對空間數(shù)據(jù)進行操作的工具。包括計算點與點的距離、構造空間數(shù)據(jù)結構以及進行空間操作的功能。程序中使用spatial模塊提供的功能構建Voronoi圖初始的三角形網(wǎng)格、進行Delaunay三角剖分，刪除多余的三角形等計算，開源庫穩(wěn)定性和通用性比自編算法穩(wěn)定、可靠，節(jié)省開發(fā)時間。使用matplotlib庫實現(xiàn)后期的繪圖功能，有助于計算完成后實時判斷構建數(shù)據(jù)的質量，調整輸入數(shù)據(jù)。使用標準圖形用戶界面（GUI）庫tkinter制作簡單的程序界面，實現(xiàn)輸入輸出文件路徑的選擇、命名、選擇模擬數(shù)據(jù)時內插的方式等，可以提高程序的可用性、便捷性。

2 案例計算

本計算案例以甘肅省定西市境內某古滑坡體為例，首先需要繪制滑坡體邊界，不需要刻意準確找到滑坡體的邊界，只需要比滑坡體大一點。當計算滑坡體體積時，兩期數(shù)據(jù)疊加，會根據(jù)高程之間的差值變化關系，自動找到填挖分界線，即為滑坡體的邊界。

滑坡體滑面上的特征點取值對滑坡體滑面的影響較大。本滑坡體選取山頭，特征點的選取可根據(jù)專家的經(jīng)驗以及滑面的情況綜合決定，由于模擬過程已經(jīng)編制好了程序，因此過程耗時極短。可以經(jīng)過調整后試算，通過試算結果調整特征點數(shù)據(jù)，也可以通過不同特征點求得的滑坡體積，根據(jù)專家判斷得到最優(yōu)解。

本文計算案例中特征點選擇示意圖如圖3所示，當然不同測區(qū)需要根據(jù)實際地形來選擇特征點，也可以試選多種方案，根據(jù)專家意見擇優(yōu)使用成果。基于滑坡體的走向，可以通過滑坡體上下緣特征點的三維坐標并通過拋物面獲得特征點相對地面的距離，從而得到特征點在滑面上的高程值。在有滑坡勘察資料的情況下，直接采用勘察鉆孔點為特征點，在mesh模型上獲取勘察鉆孔點的地面高程，扣除勘察資料顯示的滑面深度即為下滑面特征點的高程值。特征點的高程獲取方法一般有以下2種。1）可以在應急救援過程中快速模擬下滑面，獲得下滑體體積，從而為應急救援提供數(shù)據(jù)支持。2）可以較精確地獲得滑坡體下滑面，即可精確計算下滑體積。

此案例中使用第一種方法獲得特征點1、2的高程值。特征點1、2的選擇是構建滑面數(shù)據(jù)能否反應真實滑面的關鍵，本文根據(jù)筆者的經(jīng)驗選取了2個特征點。在實際案例計算中，可以根據(jù)現(xiàn)場專家對滑面的判斷，選擇任意數(shù)量的特征點并根據(jù)下滑體前后緣的情況給定特征點的Z值。根據(jù)計算結果圖判斷構建滑面與實際滑面之間的符合程度。根據(jù)可能滑面的形態(tài)模擬若干滑面成果。現(xiàn)場借鑒不同專家對各構建滑面的認可程度選取一個滑面進行下滑量計算。本案例獲取案例區(qū)域投影面積4.2萬m2，插值點間距5m，共插值出1769個高程點。通過構建滑面的1769個高程點在三維空間中采用Delaunay三角剖分法構建三角網(wǎng)。在航測成果實景三維模型中根據(jù)計算范圍給定高程點間距提取三維坐標，使用相同的方法構建三角網(wǎng)。提取高程點的間距可以是任意間距，筆者測試后建議使用與構建滑面相同的間距。根據(jù)兩期三角網(wǎng)可以計算下滑量，如圖4所示。再使用兩期表面數(shù)據(jù)獲取方格網(wǎng)角點及兩期表面間的高程差，利用方格網(wǎng)面積乘以方格網(wǎng)4個角點的平均高差獲得棱柱體積，對區(qū)域內所有棱柱體體積進行求和，即可得到測區(qū)土方，其中2個方格網(wǎng)數(shù)據(jù)見表1。

3 結語

本文外業(yè)使用精靈4RTK無人機作業(yè)獲得570張照片。三維重建獲得高精度實景三維模。繪制大致滑坡區(qū)域，提取滑坡邊界三維坐標并選擇特征點，計算特征點在畫面上的高程數(shù)據(jù)，模擬滑面高程點數(shù)據(jù)。使用兩期高程點計算滑體體積。

本方法最顯著的特點在于應急測繪中可以快速使用無人機測繪的成果估算下滑體體積。每個滑體的計算時間小于20min，可以為應急救援工作節(jié)省寶貴的時間；使用本文的方法可以根據(jù)不同專家對滑體的判斷，修改特征點數(shù)據(jù)，得到不同專家對滑體滑面情況的模擬結果，再經(jīng)過綜合判斷后選取認可度較高的特征點作為滑面模擬的輸入數(shù)據(jù)，得到最優(yōu)的滑體體積估算值；可以根據(jù)后期勘察的鉆孔數(shù)據(jù)作為特征點模擬滑面，可以較精確地估算滑體體積。精度取決于勘察鉆孔的分布和點位的數(shù)量，可以對應急階段估算值的精度進行后評估。本文方法對邊界的選擇不敏感，只要計算范圍比滑體范圍大即可。

該方法支持大多數(shù)只能獲得一期地表數(shù)據(jù)的體積計算情況，根據(jù)不同場景選擇合適的特征點，模擬任意形狀面，從而根據(jù)實際面和模擬構建面的兩期數(shù)據(jù)，計算包裹在兩期面之間的體積。

參考文獻

[1]陳富強.基于傾斜攝影與BIM技術的土方量算方法[J].地理空間信息，2024，22（1）：56-59.

[2]趙永康.消費級無人機傾斜攝影測量技術在土方量計算中的應用研究[J].測繪與空間地理信息，2024，47（1）：183-185.

[3]王山，何承錦，胡炫宇，等.基于BIM和GIS的無人機傾斜攝影技術在大體積土石方工程中的應用[J].綠色建造與智能建筑，2023（11）：66-70.