郭明珠，梁洲婕，王天成，劉 晃

(北京工業(yè)大學(xué) 城市建設(shè)學(xué)部，北京 100124)

0 引言

由于特殊的氣候、水文、地質(zhì)等自然條件，西藏地區(qū)地質(zhì)災(zāi)害頻發(fā)，具有災(zāi)害類型多、影響范圍廣、破壞程度大等特點(diǎn)，其中滑坡、泥石流、崩塌產(chǎn)生的直接危害最為嚴(yán)重。滑坡的發(fā)生常會(huì)導(dǎo)致房屋倒塌、道路中斷、河流堵塞，威脅下游水利工程建設(shè)，阻礙土地、礦產(chǎn)、水利等資源的開(kāi)發(fā)利用(張梁等，1998；曹叔尤等，2013；郭明珠等，2021)?；聻?zāi)害危及群眾生命財(cái)產(chǎn)安全，嚴(yán)重影響山區(qū)經(jīng)濟(jì)社會(huì)發(fā)展。災(zāi)害發(fā)生后，當(dāng)務(wù)之急是進(jìn)行災(zāi)后重建工作及災(zāi)害評(píng)估工作。準(zhǔn)確測(cè)量滑坡堆積體體積，有利于災(zāi)后重建工作的開(kāi)展，地質(zhì)災(zāi)害評(píng)價(jià)工作的進(jìn)行以及后期土地、礦產(chǎn)等資源的開(kāi)發(fā)。

現(xiàn)有計(jì)算滑坡堆積體體積的幾何方法主要有以下幾種：一是將滑坡視作一個(gè)半橢球，代入半橢球體積公式計(jì)算滑坡體積(Cruden，Varnes，1996)；二是基于積分近似計(jì)算原理，利用不規(guī)則形體體積的計(jì)算方法，手工繪制網(wǎng)格計(jì)算滑坡堆積體體積(熊道錕，胡濟(jì)珍，1990)；三是以統(tǒng)計(jì)關(guān)系計(jì)算為主，通過(guò)建立滑坡體體積和投影面積、滑坡體積與地震震級(jí)、滑坡體積與前期降水量、滑坡體積對(duì)數(shù)lg與滑坡最大垂直、水平距離、等效摩擦系數(shù)等經(jīng)驗(yàn)公式計(jì)算滑坡堆積體體積(Malamuda，2004；Guzzetti，2008；李軍，周成虎，2002；尉德新等，2019)；四是用平行斷面法計(jì)算滑坡堆積體體積，利用勘探剖面將地質(zhì)體分為不同塊段，分別計(jì)算并累加得到滑坡堆積體體積(王湘桂等，2012)；五是對(duì)比滑坡前后數(shù)字高程模型，計(jì)算出滑坡堆積體的體積(殷躍平等，2011；張小詠等，2013；韓健楠，2018；Chen，2005)；六是采用三維激光掃描儀或雙目立體視覺(jué)技術(shù)獲取滑坡堆積體的三維圖像信息構(gòu)造幾何模型，以堆積體頂部頂點(diǎn)為三棱錐頂點(diǎn)、堆積體底面為三棱錐底面，自上而下將堆積體剖分成若干，計(jì)算滑坡堆積體體積(陳展鵬等，2013；于歡，2018)。

現(xiàn)有的幾何計(jì)算堆積體體積主要是將堆積體擬合成多面體或拆分成三棱柱、三棱錐進(jìn)行計(jì)算，計(jì)算方法較為繁瑣。本文嘗試將石根華(1997)提出的單純形積分理論與Dlaunay三角剖分理論相結(jié)合，提出了一種新的計(jì)算堆積體體積的方法，以西藏昌都地區(qū)芒康縣索多西鄉(xiāng)的雪隆囊滑坡為例，通過(guò)Matlab的Dlaunay三角剖分算法將堆積體剖分成多個(gè)三棱錐代入單純形積分理論進(jìn)行計(jì)算。

1 研究區(qū)及滑坡堆積體概況

1.1 研究區(qū)概況

雪隆囊滑坡堆積體位于青藏高原東南緣金沙江上游的西藏昌都地區(qū)芒康縣索多西鄉(xiāng)雪隆囊村?；聟^(qū)位于金沙江深切割區(qū)的峽谷之中，氣候垂直變化顯著，屬于典型的高山峽谷地貌，山體龐大密集，山峰陡峻，兩岸地勢(shì)高差懸殊、地形起伏很大，相對(duì)高程在2 000 m以上，切割深度在800 m上下，兩岸山體高程可達(dá)2 000～4 000 m。河谷干熱、高山冷濕，屬于典型的亞熱帶氣候區(qū)，物理風(fēng)化作用極其強(qiáng)烈。滑坡所在區(qū)域降雨稀少，多年平均值為497～620 mm，主要集中在5—9月，占全年的85%以上?；聟^(qū)內(nèi)裸地、荒山草坡居多，部分區(qū)域基巖出露，整體植被覆蓋率低，不利于水土保持及斜坡穩(wěn)定性，故該流域?yàn)榛乱装l(fā)區(qū)。滑坡所在位置的區(qū)域構(gòu)造活動(dòng)強(qiáng)烈，地形地貌復(fù)雜，屬于我國(guó)一級(jí)階梯與二級(jí)階梯的過(guò)渡地帶。滑坡所處的金沙江段，江水總體流向由N20°W轉(zhuǎn)為N50°W，河谷深切，深度達(dá)700～1 000 m。同時(shí)，河谷狹長(zhǎng)，水流湍急，坡降不明顯。河流兩岸沖刷痕跡明顯，有河相沉積。

1.2 滑坡堆積體概況

雪隆囊滑坡堆積體位于蘇洼龍壩址上游5.88 km崗達(dá)村對(duì)岸，前緣臨江。由于堆積體前緣的阻隔，金沙江在此形成繞堆積體前緣的弧形河道，堆積體后緣及兩側(cè)均以干溝與山體相隔離，在地形上形似圈椅。

堆積體表面地形平緩，坡度一般為5°～20°，前緣直達(dá)江邊，與河床覆蓋層直接相接。堆積體中部為一低洼地帶，將堆積體分成兩級(jí)臺(tái)階狀，臨江部分頂面高程2 440～2 500 m，為堆積體I區(qū)；靠山側(cè)頂面高程2 500～2 620 m，為堆積體II區(qū)。根據(jù)中國(guó)電建集團(tuán)北京勘測(cè)設(shè)計(jì)研究院有限公司的地質(zhì)勘察報(bào)告(以下簡(jiǎn)稱“地質(zhì)勘察報(bào)告”)，繪制了雪隆囊滑坡剖面和平面圖(圖1、2)。

圖1 雪隆囊滑坡剖面圖(據(jù)地質(zhì)勘察報(bào)告修改)Fig.1 Section of the Xuelongnang landslide(modified according to Geological Survey Report)

圖2 雪隆囊滑坡平面圖(據(jù)地質(zhì)勘察報(bào)告修改)Fig.2 Plane of the Xuelongnang landslide(modified according to Geological Survey Report)

2 計(jì)算理論

2.1 單純形積分理論

單純形積分理論由美國(guó)數(shù)學(xué)力學(xué)家石根華(1997)提出，該方法可用邊界頂角的坐標(biāo)來(lái)計(jì)算塊體體積和重心，其基本原理如下：

一般的維單純形有+1個(gè)有序頂點(diǎn)，，，…，，表示為：

：(，，，…，0)

：(，，，…，1)

：(，，，…，2)

……

：(1，2，3，…，)

維坐標(biāo)單純形有+1個(gè)頂點(diǎn)，，，…，，表示為：

：(0，0，0，…，0)

：(1，0，0，…，0)

：(0，1，0，…，0)

……

：(0，0，0，…，1)

雅可比行列式為：

(1)

一個(gè)維單純形積分為：

(2)

式中：表示的次數(shù)。

標(biāo)準(zhǔn)的三維積分式如下：

(3)

(4)

(5)

(6)

(7)

式中：表示點(diǎn)(=1，2，3，…，)的坐標(biāo)；表示點(diǎn)(=1，2，3，…，)的坐標(biāo)；表示點(diǎn)(=1，2，3，…，)的坐標(biāo)。

2.2 Delauany三角剖分

令={，，…，}為平面域()上個(gè)離散點(diǎn)的集合，俄國(guó)數(shù)學(xué)家Dlaunay在1934年的相關(guān)研究證明：必定存在且僅存在一種三角剖分(一般稱之為Delauany三角剖分)算法，使得所有三角形的最小內(nèi)角之和最大。

一般情況下，一個(gè)完整的Voronio圖應(yīng)該由多個(gè)Voronio多邊形組成，第個(gè)Voronio多邊形的數(shù)學(xué)表達(dá)形式如下：

={∈：=-=≤=-==1，2，…，；≠，=1，2，…，}

(8)

式中：=-=表示平面域上點(diǎn)和節(jié)點(diǎn)之間的歐氏距離。

從式(8)可知，Voronio多邊形內(nèi)任意點(diǎn)到節(jié)點(diǎn)的距離比到點(diǎn)集中任何其他節(jié)點(diǎn)的距離更近。因此由節(jié)點(diǎn)和每個(gè)相鄰節(jié)點(diǎn)的垂直平分線所形成的開(kāi)式半平面的交集組成，所以，此時(shí)必為凸多邊形。一般情況下，Voronio圖的一個(gè)頂點(diǎn)同時(shí)屬于3個(gè)Voronio多邊形，每個(gè)Voronio多邊形內(nèi)有且僅有1個(gè)節(jié)點(diǎn)。連接3個(gè)共點(diǎn)Voronio多邊形分別對(duì)應(yīng)的3個(gè)節(jié)點(diǎn)則形成一個(gè)Delauany三角形，所有這樣的三角形的集合就是著名的Delauany三角剖分(丁永祥等，1994)。

2.3 單純形積分求重心及體積

根據(jù)數(shù)值流形方法的單純形積分以及重心的定義，可得到體積、重心(武艷強(qiáng)等，2011)的表達(dá)式為：

(9)

(10)

(11)

(12)

3 計(jì)算方法和結(jié)果

根據(jù)現(xiàn)場(chǎng)勘測(cè)獲得滑坡后壁、滑坡周界、堆積體表面等關(guān)鍵位置的數(shù)據(jù)(如經(jīng)緯度及高程)，確定多個(gè)關(guān)鍵點(diǎn)(表1)。由于數(shù)據(jù)量過(guò)多，篇幅有限，表1僅給出部分關(guān)鍵點(diǎn)數(shù)據(jù)，其中、值由關(guān)鍵點(diǎn)的WGS84坐標(biāo)轉(zhuǎn)換西安80坐標(biāo)得到。

表1 關(guān)鍵點(diǎn)數(shù)據(jù)Tab.1 The data of some critical points

將關(guān)鍵點(diǎn)經(jīng)緯度轉(zhuǎn)換為大地坐標(biāo)，使用Matlab的曲線擬合工具箱對(duì)關(guān)鍵點(diǎn)坐標(biāo)進(jìn)行擬合(圖3)，得出二元一次函數(shù)、二元二次函數(shù)和二元三次函數(shù)曲線擬合結(jié)果(表2)。由表2可知：①根據(jù)相關(guān)系數(shù)和均方根誤差可知，二元三次曲面擬合效果最優(yōu)，二元二次曲面次之，二元一次曲面較差；②滑動(dòng)面擬合精度的提高依賴于勘察報(bào)告以及關(guān)鍵點(diǎn)的個(gè)數(shù)。

圖3 滑動(dòng)面模擬二元一次(a)、二次(b)、三次(c)曲面圖及殘差Fig.3 The sliding surface simulating binary primary(a)，quadratic(b)，and cubic(c)surface diagram and residuals

表2 滑動(dòng)面曲線擬合結(jié)果Tab.2 Curve fitting of the sliding surface

結(jié)合地質(zhì)勘察報(bào)告結(jié)果對(duì)擬合結(jié)果進(jìn)行評(píng)價(jià)，選出最佳擬合滑動(dòng)面曲線。根據(jù)擬合滑動(dòng)面曲線算出滑動(dòng)面計(jì)算點(diǎn)坐標(biāo)，與堆積體表面數(shù)據(jù)結(jié)合，整理得出堆積體形狀及其三維坐標(biāo)。利用Matlab的三維Delauany三角剖分算法，將堆積體剖分成三棱錐(圖4)，并將三棱錐相對(duì)應(yīng)的點(diǎn)坐標(biāo)依次帶入公式(1)、(3)進(jìn)行計(jì)算，得到每個(gè)三棱錐體積，累加得到滑坡體體積。擬合計(jì)算程序如圖5所示。

圖4 三角剖分示意圖Fig.4 Schematic diagram of triangulation

圖5 滑坡堆積體體積計(jì)算程序圖Fig.5 The flow chart of the landslide volume calculation

利用單純形積分以及平行斷面法、經(jīng)驗(yàn)公式法分別計(jì)算堆積體體積，計(jì)算結(jié)果及具體計(jì)算公式見(jiàn)表3。由表3可知：

表3 不同方法得出的堆積體體積Tab.3 The volume of accumulation body obtained by different methods

(1)單純形積分計(jì)算結(jié)果與地質(zhì)勘察報(bào)告值對(duì)比，兩者相差5.63%，吻合程度較好，說(shuō)明該方法計(jì)算結(jié)果可信度較高。

(2)單純形積分計(jì)算結(jié)果與傳統(tǒng)幾何方法計(jì)算結(jié)果相差27.4%。傳統(tǒng)幾何方法計(jì)算滑坡體積時(shí)將滑坡的滑動(dòng)面當(dāng)作橢球面，未考慮滑動(dòng)面的實(shí)際情況；同時(shí)將滑坡體看作橢球形，對(duì)滑坡形態(tài)過(guò)度簡(jiǎn)化，未考慮不同滑坡之間的差異性。

(3)單純形積分計(jì)算結(jié)果與經(jīng)驗(yàn)公式法相比，分別相差5.2%、17.9%、75.3%。由于經(jīng)驗(yàn)公式法選取的地域、樣本類型、樣本個(gè)數(shù)等都不相同，不同經(jīng)驗(yàn)公式得出的結(jié)果相差較大。經(jīng)驗(yàn)公式法具有明顯區(qū)域性，不適用于所有滑坡堆積體體積計(jì)算。

(4)其他學(xué)者(陳松，2016；陳劍，崔之久，2015)對(duì)此滑坡體進(jìn)行計(jì)算，得到滑坡體體積為4.1×10m，本文方法計(jì)算結(jié)果與其相差15.9%。

(5)單純形積分計(jì)算方法精確度與關(guān)鍵點(diǎn)的個(gè)數(shù)、野外實(shí)地勘查的詳細(xì)程度以及構(gòu)建幾何模型的精確程度相關(guān)。關(guān)鍵點(diǎn)個(gè)數(shù)越多，構(gòu)建的滑坡體幾何模型越貼近實(shí)際，計(jì)算結(jié)果越準(zhǔn)確。

4 結(jié)論

本文利用Matlab軟件對(duì)滑動(dòng)面曲線進(jìn)行擬合得出完整堆積體，將滑坡堆積體剖分成三棱錐，利用單純形積分計(jì)算堆積體體積。得出結(jié)論如下：

(1)相較于其它方法，單純形積分計(jì)算方法準(zhǔn)確度較高，可用于滑坡堆積體體積的計(jì)算。

(2)結(jié)合Delauany三角剖分算法，將任意多面體剖分成三棱錐并代入單純形積分進(jìn)行計(jì)算的方法可用于雙目立體視覺(jué)技術(shù)或三維激光掃描技術(shù)建立滑坡堆積體幾何模型后的滑坡體體積計(jì)算，同時(shí)還可用于任意不規(guī)則形體體積計(jì)算。

(3)基于Matlab現(xiàn)有的視覺(jué)圖像識(shí)別處理技術(shù)，凹陷塊體的外部邊界無(wú)法準(zhǔn)確識(shí)別，計(jì)算凹陷塊體體積時(shí)需考慮缺失部分，計(jì)算后進(jìn)行相應(yīng)刪減。

感謝中國(guó)電建集團(tuán)北京勘測(cè)設(shè)計(jì)研究院有限公司提供雪隆囊滑坡堆積體(崗達(dá)堆積體)的勘探資料和體積估算。