周小荃，余宏亮，魏玉杰，胡 節(jié)，蔡崇法 ※

（1.華中農(nóng)業(yè)大學(xué)資源與環(huán)境學(xué)院，武漢 430070；2.農(nóng)業(yè)部長(zhǎng)江中下游耕地保育重點(diǎn)實(shí)驗(yàn)室，武漢 430070；3.湖北省電力勘測(cè)設(shè)計(jì)院有限公司，武漢 430040）

0 引 言

崩崗是指在水力和重力的綜合作用下，山坡土石體受破壞而崩塌和受沖刷的侵蝕現(xiàn)象[1]。崩崗主要由崩壁，崩積堆和洪（沖）積扇組成[2]，崩崗侵蝕復(fù)雜多變，包含重力崩塌、雨滴擊濺、片流沖刷、徑流侵蝕等多種形式[3]，所以對(duì)崩崗侵蝕量的精確監(jiān)測(cè)對(duì)研究崩崗的發(fā)生與發(fā)育機(jī)理有著重要意義。然而崩崗流域地貌具有深溝陡坡，形態(tài)復(fù)雜多變，且存在其特有的掏蝕現(xiàn)象[4]，使得監(jiān)測(cè)崩崗的侵蝕量變化十分困難。因此，目前對(duì)崩崗侵蝕量監(jiān)測(cè)的研究非常少。1

在侵蝕地貌的監(jiān)測(cè)中，研究人員嘗試將測(cè)繪儀器和設(shè)備應(yīng)用于量化侵蝕過(guò)程，通過(guò)獲取同一區(qū)域不同時(shí)期重復(fù)的地形測(cè)量數(shù)據(jù)，生成多時(shí)相的數(shù)字高程模型（DEM），并建立相應(yīng)的算法用以模擬侵蝕過(guò)程、估算侵蝕量和預(yù)測(cè)地形環(huán)境變化[5-6]。根據(jù)獲取地形數(shù)據(jù)方式的不同，可將現(xiàn)有的監(jiān)測(cè)質(zhì)量運(yùn)移的方法分為3 種：基于全球定位系統(tǒng)（GPS）[7]地形測(cè)量法、數(shù)字航空攝影測(cè)量法(DAP)和三維激光掃描法[8]。在崩崗侵蝕監(jiān)測(cè)中，單點(diǎn)GPS 測(cè)量可用于獲得精確的參考數(shù)據(jù)或者有限區(qū)域的地形信息，這種方法可以保證單點(diǎn)精度，但形成DEM 的分辨率很低[7]。三維激光掃描法可以獲得精度極高的點(diǎn)云數(shù)據(jù)，現(xiàn)有的技術(shù)已經(jīng)可以達(dá)到毫米級(jí)的位移監(jiān)測(cè)[9]，但設(shè)備價(jià)格昂貴且儀器龐大沉重，需要固定且平穩(wěn)的觀測(cè)點(diǎn)，在野外工作中十分不便于使用，同時(shí)點(diǎn)云數(shù)據(jù)會(huì)受地形和植被等遮擋影響，形成掃描死角從而產(chǎn)生點(diǎn)云空洞，使用時(shí)多站點(diǎn)大量的點(diǎn)云數(shù)據(jù)拼接難度較大、耗時(shí)費(fèi)力[10]。相反，航空攝影可以快速、連續(xù)地覆蓋大尺度的區(qū)域，獲取高分辨率的正射影像[11]，但是樣本密度不足，無(wú)法到達(dá)精度，且容易被植被和地形因素影響[12]。

近年來(lái)，無(wú)人機(jī)(unmanned aerial vehicle, UAV)航空攝影測(cè)量[13]作為一項(xiàng)測(cè)繪新技術(shù)，使得地形數(shù)據(jù)的獲取快捷方便，監(jiān)測(cè)范圍更大效率更高，現(xiàn)在已經(jīng)被廣泛地運(yùn)用在了很多需要實(shí)景三維建模的領(lǐng)域[14-15]，在崩崗侵蝕方面也開(kāi)始萌芽[16]，但是在崩崗崩壁崩塌侵蝕量監(jiān)測(cè)方面運(yùn)用并不成熟。運(yùn)用在崩崗侵蝕監(jiān)測(cè)中的優(yōu)勢(shì)是：1）無(wú)人機(jī)可以貼近崩壁、溝道等人為無(wú)法靠近的區(qū)域飛行，從而獲取全方位、多角度的影像信息，避免視野盲區(qū)，解決崩壁崩塌侵蝕數(shù)據(jù)無(wú)法測(cè)算的問(wèn)題，獲取完整的地形信息甚至細(xì)節(jié)紋理[17]；2）影像經(jīng)過(guò)同名點(diǎn)匹配生成的點(diǎn)云數(shù)據(jù)量遠(yuǎn)小于三維激光掃描法獲取的數(shù)據(jù)量，因而數(shù)據(jù)處理、分析更加易于操作；3）相比三維激光掃描儀，無(wú)人機(jī)價(jià)格便宜、易于攜帶，更適合野外測(cè)量。

本文對(duì)研究區(qū)崩崗進(jìn)行低空傾斜攝影獲取多角度航空影像，通過(guò)空三加密處理生成目標(biāo)崩崗崩壁和崩積堆的三維點(diǎn)云模型。為了達(dá)到最優(yōu)效果，得到更多有用的描述侵蝕過(guò)程的數(shù)據(jù)，使用了兩種方法配合使用來(lái)測(cè)算崩崗侵蝕量和提取崩崗的空間數(shù)據(jù)，分析在監(jiān)測(cè)周期內(nèi)崩崗侵蝕的變化規(guī)律；最后使用三維激光掃描儀與無(wú)人機(jī)同期實(shí)測(cè)目標(biāo)崩崗作為參考數(shù)據(jù)驗(yàn)證精度。克服航空及傳統(tǒng)測(cè)量無(wú)法對(duì)垂直崖壁進(jìn)行定量描述的限制，為崩崗侵蝕量定量研究提供比三維激光掃描更簡(jiǎn)單快捷的方法，使深入研究崩崗侵蝕成因機(jī)理成為可能[18]。

1 材料與方法

1.1 研究區(qū)概況

試驗(yàn)監(jiān)測(cè)崩崗位于湖北省咸寧市通城縣五里鎮(zhèn)程鳳村（113°43'26"E，29°12'18"N），屬于北亞熱帶季風(fēng)氣候區(qū)，年平均氣溫15.5～16.7 ℃，雨季集中在4—8 月，年平均降水量1 450～1 600 mm。該崩崗為活動(dòng)型崩崗，崩壁向分水嶺侵蝕后退仍在進(jìn)行中，具有典型瓢形崩崗的各種特征，土壤類(lèi)型為棕紅壤，崩崗內(nèi)部植被稀少。

1.2 數(shù)據(jù)獲取

1.2.1 控制點(diǎn)的布設(shè)與選取

在目標(biāo)崩崗內(nèi)采用大地坐標(biāo)系布設(shè)控制點(diǎn)，首次測(cè)量前使用Trimble R4 GNSS 差分GPS 布設(shè)控制網(wǎng)，為避免崩崗快速侵蝕對(duì)控制點(diǎn)坐標(biāo)的影響，共布設(shè)10 個(gè)控制點(diǎn)（理論上每平方千米3 個(gè)控制點(diǎn)即可實(shí)現(xiàn)坐標(biāo)定位），選擇了崩崗附近監(jiān)測(cè)室的4 個(gè)角點(diǎn)和崩壁、崩積堆打入六根1 m 長(zhǎng)的鐵桿并留出0.3 m 的高度，作為固定點(diǎn)位的控制點(diǎn)。

1.2.2 無(wú)人機(jī)數(shù)據(jù)的獲取

數(shù)據(jù)采集的時(shí)間為2018 年6 月到9 月雨季期間，在每月選擇天氣晴朗、風(fēng)力較小的時(shí)間段進(jìn)行。無(wú)人機(jī)型號(hào)為大疆精靈4pro 準(zhǔn)專(zhuān)業(yè)無(wú)人機(jī)，該機(jī)型自帶的1 英寸2 000 萬(wàn)像素CMOS 影像傳感器與云臺(tái)集成，配有主、副 2 個(gè)遙控器云臺(tái)，可控轉(zhuǎn)動(dòng)范圍為俯仰-90°至+30°可實(shí)現(xiàn)全方位的崩崗監(jiān)測(cè)，無(wú)人機(jī)及相機(jī)的基本參數(shù)如表 1 所示。

表1 無(wú)人機(jī)及相機(jī)的基本參數(shù) Table 1 Basic parameters of UAV and camera

在目標(biāo)崩崗操控?zé)o人機(jī)如圖1 所示，從距離集水區(qū)頂部20 m 處開(kāi)始操控云臺(tái)，在垂直方向上以-10°、-30°、-50°、-70°、-90°，這5 個(gè)固定角度，水平方向上360°進(jìn)行環(huán)形全景拍攝，在一個(gè)高度拍攝完畢之后降低5 m 再次進(jìn)行同樣操作；當(dāng)無(wú)人機(jī)降低到崩崗集水區(qū)高度以下時(shí)，操控云臺(tái)在垂直方向改為+10°、-10°、-30°、-50°、-70°、-90°這6 個(gè)固定角度。從而在每一個(gè)高度都能獲得崩壁的完整影像數(shù)據(jù)，以確保獲取目標(biāo)崩崗不同位置和不同角度的影像，根據(jù)具體地形的不同可以變化拍攝方法，本次實(shí)驗(yàn)最高高度為高于集水區(qū)20 m 處，最低高度為距離溝底5 m 處。若崩崗有很深的侵蝕溝，則操控?zé)o人機(jī)時(shí)需緩慢向溝谷內(nèi)移動(dòng)并重復(fù)以上操作，移動(dòng)過(guò)程中也要進(jìn)行拍攝。在操作過(guò)程中確保在同一高度采集的相鄰影像重疊率超過(guò)90%，不同高度采集的影像重疊率超過(guò) 60%，在保證安全距離的情況下按照以上操作進(jìn)行近距攝影，從而獲取崩崗全部影像信息。

圖1 無(wú)人機(jī)傾斜攝影示意圖 Fig.1 Diagram of UAV oblique photography

依照此方案，對(duì)目標(biāo)崩崗在2018 年6 月15 日，2018年7 月4 日，2018 年8 月1 日，2018 年8 月13 日，2018 年9 月30 日進(jìn)行了5 次無(wú)人機(jī)攝影。

1.2.3 參考數(shù)據(jù)的獲取

依據(jù)劉希林等[19]提出的方法， 利用 Z+F IMAGER5006h 三維激光掃描儀與無(wú)人機(jī)同期使用如圖2所示，同時(shí)獲取崩崗高精度點(diǎn)云數(shù)據(jù)作為參考。測(cè)量的關(guān)鍵步驟如下：

圖2 三維激光掃描儀測(cè)量目標(biāo)崩崗 Fig.2 Laser scanner for measuring Benggang of target

1） 控制點(diǎn)布設(shè)。使用與采集無(wú)人機(jī)數(shù)據(jù)一致的控制點(diǎn)。

2）數(shù)據(jù)拼接與配準(zhǔn)。將目標(biāo)崩崗分為2 站測(cè)量，再利用海聯(lián)達(dá)三維軟件進(jìn)行多站點(diǎn)拼接，并使用控制點(diǎn)和標(biāo)靶的大地坐標(biāo)對(duì)儀器內(nèi)置坐標(biāo)進(jìn)行轉(zhuǎn)換，將所有點(diǎn)云數(shù)據(jù)統(tǒng)一在真實(shí)大地坐標(biāo)系下，以實(shí)現(xiàn)不同站點(diǎn)點(diǎn)云數(shù)據(jù)的拼接和不同時(shí)期點(diǎn)云數(shù)據(jù)的配準(zhǔn)。

3）體積的計(jì)算。為了達(dá)到能夠作為參考數(shù)據(jù)的毫米級(jí)精度，不抽希點(diǎn)云，通過(guò)Arcgis10.0 的3D 分析模塊，采樣密度大，使用LAS 數(shù)據(jù)集轉(zhuǎn)TIN（不規(guī)則三角網(wǎng)）形成表面。使用ArcGIS 三維分析模塊中的 Area and Volume 工具能夠計(jì)算出TIN 表面到任意高程水平面垂直投影下的不規(guī)則立柱體的體積。本文以程鳳崩崗最低點(diǎn)主溝溝口高程為底層，從崩崗表面向下作垂直投影形成不規(guī)則立柱體，計(jì)算前后兩次立柱體的體積之差，即為崩崗侵蝕或沉積量。

1.3 研究方法

如圖3 所示，本文對(duì)無(wú)人機(jī)獲取的影像數(shù)據(jù)運(yùn)用Photoscan 進(jìn)行空三加密處理，生成目標(biāo)崩崗的三維點(diǎn)云模型；運(yùn)用SOR 去噪法和人工手動(dòng)去除植被影響，得到崩壁和崩積堆的完整點(diǎn)云模型；利用點(diǎn)云模型生成數(shù)字地形模型（DTM），提取崩崗地形信息，運(yùn)用多時(shí)相連續(xù)DEM 相減獲得DOD（Difference of DEMs），從而獲取監(jiān)測(cè)周期內(nèi)崩崗連續(xù)的侵蝕沉積動(dòng)態(tài)變化并找到侵蝕嚴(yán)重的部位，再使用2.5D 體積測(cè)算方法來(lái)細(xì)化侵蝕嚴(yán)重的崩壁和溝頭部位的侵蝕量作為補(bǔ)充，最終獲得監(jiān)測(cè)期內(nèi)的總侵蝕/沉積量。并以三維激光掃描儀得到崩崗侵蝕體積數(shù)據(jù)作為參考值，對(duì)崩崗的侵蝕與沉積量測(cè)算結(jié)果進(jìn)行精度檢驗(yàn)，分析在監(jiān)測(cè)周期內(nèi)短期崩崗侵蝕過(guò)程。

圖3 無(wú)人機(jī)數(shù)據(jù)監(jiān)測(cè)崩崗侵蝕過(guò)程技術(shù)路線圖 Fig.3 Technical roadmap of UAV data monitoring collapse erosion process

1.3.1 無(wú)人機(jī)數(shù)據(jù)處理

以Photoscan 1.2.5 攝影測(cè)量軟件為工具，分別對(duì)同一目標(biāo)崩崗5 期影像數(shù)據(jù)進(jìn)行處理。以2018 年6 月15 日的崩崗拍攝影像為例對(duì)主要處理步驟進(jìn)行如下說(shuō)明：

1） 空三加密

導(dǎo)入影像數(shù)據(jù)和控制點(diǎn)數(shù)據(jù)，將參考系設(shè)置為WGS-84 坐標(biāo)系，再導(dǎo)入相機(jī)校驗(yàn)參數(shù)文件完成屬性設(shè)置。通過(guò)Photoscan 軟件對(duì)所有照片進(jìn)行空三加密處理，軟件自動(dòng)恢復(fù)每每張影像片的方位姿態(tài)[20]。

2） 點(diǎn)云生成

經(jīng)過(guò)多次試驗(yàn)嘗試，為了避免數(shù)據(jù)量過(guò)大噪點(diǎn)過(guò)多而導(dǎo)致實(shí)驗(yàn)數(shù)據(jù)處理困難，最終確定將點(diǎn)云生成的“圖像比例”參數(shù)設(shè)置為“1/2”，“匹配最低數(shù)值”設(shè)置為“3”，“點(diǎn)密度”設(shè)置為“中等”。

3） 點(diǎn)云去噪與輸出

利用觀察法人工手動(dòng)框選的方式將生成的點(diǎn)云數(shù)據(jù)中植被、天空、建筑等多余噪聲點(diǎn)刪除進(jìn)行初步去噪，再使用SOR（statistical outlier removal）去噪法精細(xì)去除相機(jī)拍攝產(chǎn)生的噪點(diǎn)，最后輸出目標(biāo)崩崗的三維點(diǎn)云模型圖4a，去噪后的正視圖、俯視圖、斜視圖如圖4b-4d所示。以LAS1.2 格式對(duì)包含坐標(biāo)和顏色信息的點(diǎn)云成果進(jìn)行輸出，LAS 格式便于使用 CloudCompare 和Arcgis10.4.1 進(jìn)行多時(shí)相點(diǎn)云模型對(duì)齊和地形數(shù)據(jù)提取分析。2018 年6 月15 日的拍攝影像建立的崩崗三維點(diǎn)云模型成果的點(diǎn)云數(shù)為14 419 534 個(gè)，去噪后只留下崩崗區(qū)的點(diǎn)云數(shù)為7 419 534 個(gè)。

圖4 崩崗三維點(diǎn)云模型 Fig.4 3D point cloud model of Benggang

4） 不同時(shí)期點(diǎn)云對(duì)齊

由于測(cè)算侵蝕量需要不同時(shí)期數(shù)據(jù)的對(duì)比，將導(dǎo)出的不同時(shí)期點(diǎn)云模型的LAS 文件導(dǎo)入CloudCompare 軟件中，利用CloudCompare 軟件的對(duì)齊拼接工具手動(dòng)選擇已經(jīng)標(biāo)定的控制點(diǎn)進(jìn)行點(diǎn)云配準(zhǔn)，使用點(diǎn)與點(diǎn)的距離工具來(lái)判定兩期點(diǎn)云的復(fù)現(xiàn)性，得到平均均方根誤差為0.000 37 m，可以達(dá)到監(jiān)測(cè)崩崗侵蝕要求。

1.3.2 侵蝕量測(cè)算方法

計(jì)算崩崗侵蝕量的思路是將后一時(shí)期的點(diǎn)云數(shù)據(jù)或柵格數(shù)據(jù)與前一時(shí)期的進(jìn)行對(duì)比，計(jì)算點(diǎn)云模型差值的體積。沉積是侵蝕土體在監(jiān)測(cè)期間沒(méi)有被運(yùn)移出崩崗地貌而形成，所以最終得到的正值為沉積的體積量而負(fù)值為侵蝕的體積量，兩者相加得到在監(jiān)測(cè)周期內(nèi)目標(biāo)崩崗范圍中總共的侵蝕或沉積量。

1） 基于多時(shí)相DEM 差分方法

依據(jù)文獻(xiàn)[21]曾用方法來(lái)計(jì)算。多期地形數(shù)據(jù)中的變化通過(guò)創(chuàng)建一個(gè)DEM 差值（DOD）來(lái)計(jì)算。以高程值作為每一個(gè)柵格的計(jì)算標(biāo)準(zhǔn)，用每一個(gè)柵格高程的差值ZΔ 與柵格大小的乘積表示體積量的變化。總的體積變化 0VΔ 計(jì)算為[22]

式中AV為參與計(jì)算的總面積（m2）；為在柵格（i，j）處的高程變化（m）；n 為在總面積AV中的柵格總數(shù)。

2）基于點(diǎn)云2.5D 體積測(cè)量法

運(yùn)用三維點(diǎn)云數(shù)據(jù)計(jì)算體積是很難實(shí)現(xiàn)的，因?yàn)槿绾味x物體的邊界是很難解決的問(wèn)題，一種方法是提取三維點(diǎn)云的凸殼（包絡(luò)體），然后再進(jìn)行計(jì)算。但是在點(diǎn)云存在孔、洞的情況下提取和計(jì)算就變得十分復(fù)雜，所以將點(diǎn)云降維成2.5D 即平面的函數(shù)值計(jì)算體積。

本文將引入一種基于開(kāi)源點(diǎn)云處理軟件cloudcompare 的2.5D 體積算法解決崩壁陡坡上崩塌侵蝕和坡腳掏蝕無(wú)法計(jì)算的問(wèn)題，將2.5D 點(diǎn)云底面劃分成離散的網(wǎng)格如圖5 所示，計(jì)算每個(gè)網(wǎng)格對(duì)應(yīng)頂面點(diǎn)云單元的體積并相加求和得到 VΔ 計(jì)算方程為

式中S 為網(wǎng)格劃分步長(zhǎng)step（m2）；ΔH 為兩期點(diǎn)云中從定義的底面點(diǎn)云到頂面點(diǎn)云在設(shè)定方向上的距離(m)。

圖5 點(diǎn)云2.5D 體積計(jì)算模型 Fig.5 Point cloud 2.5D volume calculation Model

式中V 為1 個(gè)監(jiān)測(cè)周期內(nèi)的總侵蝕/沉積量。

3）不確定度計(jì)算

數(shù)字高程模型（DEM）的不確定性因素取決于點(diǎn)云、點(diǎn)密度、表面復(fù)雜性和插值方法的準(zhǔn)確性，這些不確定因素都會(huì)對(duì)DEM 的精度產(chǎn)生影響，而DOD 是計(jì)算兩組DEM 數(shù)據(jù)的差值，所以需要考慮DEM 的不確定性并進(jìn)行誤差分析才能得到可靠地結(jié)果。

通常情況下通過(guò)平差得到的控制點(diǎn)的總均方根誤差用來(lái)評(píng)估轉(zhuǎn)換誤差數(shù)值[23]。由兩組DEM 將不確定度傳播到DOD 中的共同誤差估算公式如下[24]：

在地貌分析領(lǐng)域關(guān)于DOD 不確定性閾值估算的研究已有較為成熟的研究，根據(jù)本試驗(yàn)的特點(diǎn)選用最低等級(jí)的檢測(cè)方法（LODmin），即假定DOD 不能準(zhǔn)確的計(jì)算高程范圍在- LODmin～+ LODmin之間的土體變化，該范圍內(nèi)有許多不正確的高程變異點(diǎn)，在最終的計(jì)算結(jié)果中刪去。在本次試驗(yàn)中設(shè)定[25]。

4）泥沙量的計(jì)算

將得到的侵蝕體積根據(jù)不同部位不同的容重進(jìn)行產(chǎn)沙量的計(jì)算，公式如下[26]：

式中SP 為產(chǎn)沙量（kg），ED 為高程差之和(m)，GR 為柵格分辨率(m)，Bd=沉積體或侵蝕體的容重(kg/m3)，V為侵蝕或沉積的總量（m3）。

2 結(jié)果與分析

5 次飛行任務(wù)后采集的影像結(jié)果如下表2 所示。每次采集的影像在300 張左右，圖像校準(zhǔn)成功率接近100%。地面分辨率最優(yōu)可以達(dá)到0.016 5 m 到達(dá)厘米級(jí)。地面控制點(diǎn)在所有圖像中都至少有2 個(gè)可以標(biāo)記。5 組地面控制點(diǎn)的平均重投影誤差為0.00 679 m，總均方根誤差的平均值為29.46 mm。根據(jù)式（3）將總均方根誤差的平均值帶入每期DEM 中計(jì)算得出LODmin為42 mm。

表2 拍攝影像的基本信息 Table 2 Basic information on photographing images

2.1 崩崗數(shù)字地形模型的構(gòu)建與地形信息的提取

地形因素是影響崩崗侵蝕的關(guān)鍵因素，崩崗侵蝕發(fā)生的方式與速率都與之相關(guān)，也是崩崗發(fā)育過(guò)程的重要判斷依據(jù)[27]。本文將已經(jīng)對(duì)齊的不同時(shí)期的崩崗點(diǎn)云LAS 文件導(dǎo)入Arcgis，使用3D Analyst 模塊下的LAS 轉(zhuǎn)多點(diǎn)要素，再使用反距離權(quán)重法構(gòu)建分辨率為1 cm 的DTM 柵格，由于獲得的點(diǎn)云數(shù)量足夠多，該方法可以快速便捷的得出高精度柵格。以2018 年6 月15 日測(cè)量數(shù)據(jù)為例，將DTM 輸出為DEM并生成等高線（圖6a），通過(guò)Spatial Analyst 工具提取坡度，坡向和匯流信息（圖6b～6d)），實(shí)現(xiàn)對(duì)崩崗地形的量化分析。計(jì)算得出，崩崗高程在246.782～261.604 m 相對(duì)高差14.822 m，崩崗壁高13.55 m，平均深度20 m，崩口寬1.8 m，崩壁沿線內(nèi)部面積236.52 m2屬小型崩崗。5 條切溝呈爪狀分布侵蝕主溝1 條，溝道長(zhǎng)12.45 m，溝道平均寬度1.14 m。由于崩崗較小且陡峭，坡度大于60°近乎直立的崩壁面積占總面積的33.35%，坡度大于25°的陡坡面積占85.7%。

圖6 地形信息圖 Fig.6 Topographic information map

2.2 崩崗侵蝕量的測(cè)算與分析

1）多時(shí)相DEM 差值

運(yùn)用Arcgis 柵格計(jì)算器將前一期的表面DEM 從后一期的表面DEM 中減去建立DOD 圖，利用柵格計(jì)算器的編程工具相應(yīng)的編寫(xiě)代碼使DOD 計(jì)算時(shí)忽略最小閾值±LODmin，從圖7 中可直觀的看出崩崗侵蝕強(qiáng)度的空間分布，綠色正值表示沉積區(qū)域，而紅色負(fù)值表示侵蝕區(qū)域。運(yùn)用Arcgis 的Spatial Analyst 工具的填挖方功能如式（1）計(jì)算崩崗整體侵蝕量 0VΔ 見(jiàn)表3。

由圖7 可知，研究區(qū)崩崗的主要侵蝕形式為溝道下切，最強(qiáng)烈的侵蝕區(qū)發(fā)生在溝道旁崩壁的坡腳處,可以很清楚地看出2018-07-04 至2018-08-01 崩壁土體崩落散落在溝道內(nèi)形成沉積，2018-08-13 至2018-09-30 由于溝道深切溝道兩側(cè)的土體失穩(wěn)向溝道內(nèi)部垮塌。

圖7 DODLOD 侵蝕量計(jì)算圖 Fig.7 Calculation of erosion amount of DODLOD

但是在崩崗這種特殊的侵蝕地貌中崩壁陡峭而且由于淘蝕作用向內(nèi)凹[4]，如圖8 所示在地形坡度大于45°的區(qū)域，由于X、Y 方向的微小的水平誤差，會(huì)造成X/S（其中S 為坡度）出現(xiàn)超過(guò)1 m 的垂直誤差[28]，本研究中受崩壁上層懸掛土體凸出影響，計(jì)算會(huì)出現(xiàn)超過(guò)4 m 的誤差（圖7a～圖7d 中左下角的紅綠斑點(diǎn)處），而在未發(fā)生崩塌的崩壁上這種規(guī)模的實(shí)際變化是不可能的，并且在實(shí)際發(fā)生崩塌的區(qū)域由于坡腳土體受水蝕影響向壁內(nèi)凹陷，如果用高程差來(lái)計(jì)算體積的話就會(huì)使計(jì)算的侵蝕量偏大，所以為了準(zhǔn)確計(jì)算需要在DEM中將這些部位刪去。

所以本文使用CloudCompare 軟件根據(jù)圖7 將主要發(fā)生侵蝕的溝頭處和左側(cè)壁分別分割出來(lái)，并且使用2.5D體積計(jì)算功能來(lái)測(cè)算壁侵蝕的體積為DOD 數(shù)據(jù)做重要的補(bǔ)充，如圖8 所示分割后的點(diǎn)云模型包括全部的崩積堆和部分崩壁，使用DOD 法制作侵蝕量分布圖重新計(jì)算侵蝕量，記錄為ΔV1見(jiàn)表3。

如圖8b～圖8e 所示分割后的DOD 圖能更加清晰和精確的監(jiān)測(cè)崩崗侵蝕過(guò)程的變化，2018-07-04—2018-08-01 崩積堆發(fā)生了坡面侵蝕也可以清晰的看出，2018-06—15—2018-07-04 時(shí)期開(kāi)始溝道不斷的下切侵蝕直到2018-08-13—2018-09-30 監(jiān)測(cè)期間溝道側(cè)壁發(fā)生了崩塌現(xiàn)象，崩落的土體在溝道內(nèi)堆積。

2）基于點(diǎn)云2.5D 體積測(cè)量崩壁侵蝕體積

將崩壁方向定義為底面進(jìn)行點(diǎn)云2.5d 體積測(cè)量，運(yùn)用 CloudCompare 軟件在體積計(jì)算界面中定義地面（Ground）及頂面（Ceil）：這里后一期點(diǎn)云選為常量，值為0；頂面選為前一期點(diǎn)云，空單元（不包含數(shù)據(jù)點(diǎn)）不參與體積計(jì)算，這里將其選為leave empty，本文選擇0.01step 即1cm 的單元。本試驗(yàn)側(cè)壁崩塌現(xiàn)象發(fā)生在2018-07-04—2018-08-01 之間（圖9a～9a b），溝頭溯源侵蝕在2018-0813 到2018-09-30 之間比較明顯（圖9c～圖d），最終可以達(dá)到很好的可視化效果，紅色正值表示侵蝕區(qū)域，藍(lán)色負(fù)值表示沉積區(qū)域（在計(jì)算時(shí)設(shè)侵蝕為負(fù)沉積為正），體積測(cè)算方式如公式（2），如圖9b 只分割出了發(fā)生崩塌的側(cè)壁侵蝕量記為ΔV2，圖9d 為溝頭處坡度大于80°的陡坡侵蝕量記為ΔV3 結(jié)果見(jiàn)表3。

圖8 分割后的點(diǎn)云模型與分割后的侵蝕分布圖 Fig.8 Segmented point cloud model and DOD

表3 侵蝕量計(jì)算結(jié)果 Table 3 Erosion calculation results

圖9 2.5D 體積侵蝕量圖 Fig 9 2.5D volume erosion chart

3）總侵蝕和泥沙量的計(jì)算

根據(jù)以往的研究[29]得知目標(biāo)崩崗崩壁不同層次土體的容重和崩積堆的土體容重得到平均干容重為1.37 g/cm3。將ΔV1、ΔV2、ΔV3相加得到最終的總侵蝕量V，再將V 代入式（4）進(jìn)行計(jì)算得到每個(gè)監(jiān)測(cè)周期侵蝕/沉積的總泥沙量見(jiàn)表3，最終得到該崩崗在監(jiān)測(cè)期2018-06-15—2018-09-30期間近4 個(gè)月流失的總泥沙量為5356.289 kg。

4）監(jiān)測(cè)期內(nèi)崩崗侵蝕過(guò)程

降雨是崩崗侵蝕的主要驅(qū)動(dòng)力[30]，為進(jìn)一步分析崩崗侵蝕過(guò)程與降雨量的關(guān)系，本次試驗(yàn)獲取了武漢立方在通城程鳳村監(jiān)測(cè)的降雨量數(shù)據(jù)如表4。

由表4 可知崩塌侵蝕量相對(duì)較小，所以侵蝕量的主要來(lái)源是崩積堆，而崩積堆土體松散所以受降雨的影響很大，但是降雨量最多的監(jiān)測(cè)時(shí)段卻不是侵蝕量最大的 時(shí)段，造成這種現(xiàn)象的原因是連續(xù)的降雨發(fā)生在6 月30日到7 月8 日之間，前一時(shí)期的大量降雨使土體含水率增大更易侵蝕[31]，在連續(xù)降雨的條件下土壤抗侵蝕能力有先增大后減小的趨勢(shì)[32]，所以在圖8 中可以明顯看出后一監(jiān)測(cè)期發(fā)生了更嚴(yán)重的侵蝕現(xiàn)象。在2018-06-15—2018-07-04 監(jiān)測(cè)期內(nèi)的降雨量是最大的，發(fā)生了一次降雨量大于80 mm 的大暴雨，由圖8b-圖8d 所示監(jiān)測(cè)期內(nèi)主要的侵蝕形式是水力侵蝕為主的溝道下切加深，主要發(fā)生區(qū)域?yàn)闇系纼?nèi)部，其他部位距離溝道越遠(yuǎn)侵蝕強(qiáng)度逐漸減弱，崩積堆發(fā)生坡面片狀侵蝕，崩壁只發(fā)生了少量的剝落；而在2018-07-04—2018-08-01 期間崩壁發(fā)生大量崩塌，崩落的土體散落在溝道上形成沉積，崩塌由于溝道受徑流影響側(cè)蝕加寬侵蝕坡腳使崩壁達(dá)到崩塌的條件最終形成，說(shuō)明崩壁崩塌事件與降雨量沒(méi)有固定的關(guān)系；2018-08-01—2018-08-13 的總侵蝕量值為沉積，因?yàn)檫@個(gè)監(jiān)測(cè)周期相對(duì)較短且期間總降雨量為22 mm 所以幾乎沒(méi)有侵蝕發(fā)生，但是在之前崩落的土體在崩崗系統(tǒng)內(nèi)移動(dòng)造成了一些體積變化，這說(shuō)明短期監(jiān)測(cè)很難發(fā)現(xiàn)規(guī)律性，只能作為長(zhǎng)期監(jiān)測(cè)的有效補(bǔ)充[28]；由圖9d 可知2018-08-13—2018-09-30 期間溝頭后退發(fā)生在崩壁的跌水口處和跛腳處，這是由于集水區(qū)徑流從崩壁跌水口沖出發(fā)生溯源侵蝕，沿崩壁向下沖刷坡腳形成淘蝕現(xiàn)象，且由圖8e 可以看出坡腳處溝道的側(cè)壁發(fā)生了崩塌現(xiàn)象。

表4 降雨量與侵蝕量數(shù)據(jù)對(duì)照表 Table 4 Rainfall and erosion data comparison

所以在監(jiān)測(cè)周期2018-06-15—2018-09-30 內(nèi)，目標(biāo)崩崗先以水力侵蝕為主，溝道和崩積堆為主要侵蝕部位；到達(dá)崩塌臨界條件后坡腳處發(fā)生了重力侵蝕崩塌；之后又以水力侵蝕為主，發(fā)生侵蝕的部位都在溝道附近，溝頭處的溯源侵蝕量隨著崩崗侵蝕的發(fā)育持續(xù)增加，溝道側(cè)的坡腳再次發(fā)生崩塌，可見(jiàn)崩崗侵蝕是一個(gè)水力侵蝕與重力侵蝕交替發(fā)生的侵蝕地貌，這種循環(huán)導(dǎo)致了崩崗侵蝕的發(fā)育。

2.3 測(cè)算結(jié)果精度檢驗(yàn)

以三維激光掃描儀實(shí)測(cè)獲取的試驗(yàn)崩崗的三維點(diǎn)云為基礎(chǔ)，根據(jù)文獻(xiàn)[10]的方法如1.22 所述將精度調(diào)整至0.01m的TIN 柵格獲得TIN 表面，分別求得 5 次監(jiān)測(cè)日期崩崗體積的實(shí)測(cè)值分別為929.875 2，928.922 5，925.985，926.960 6，925.706 4 m3，將前后兩期體積差作為監(jiān)測(cè)期內(nèi)的總侵蝕或沉積量體積，以此為參考值，對(duì)本文方法測(cè)算崩崗侵蝕/沉積體積的結(jié)果進(jìn)行精度檢驗(yàn)，結(jié)果如表5 所示。

表5 崩崗侵蝕量計(jì)算值精度 Table 5 Accuracy of calculation value of erosion amount of Benggang

由表5 知，用三維激光掃描儀測(cè)量后使用TIN 構(gòu)建投影面體積法測(cè)算的侵蝕量參考值與本試驗(yàn)的測(cè)算值相差很小，相對(duì)誤差在5.58%～14.71%的范圍內(nèi)，絕對(duì)誤差在0.053 2～0.3033 m3之間在年侵蝕量巨大的崩崗地貌中具有很好的準(zhǔn)確度，說(shuō)明無(wú)人機(jī)低空傾斜航空攝影監(jiān)測(cè)能夠準(zhǔn)確計(jì)算出崩崗侵蝕量。

3 討 論

以往在對(duì)崩崗侵蝕的監(jiān)測(cè)中，排樁法和在流域出口修建谷坊監(jiān)測(cè)泥沙的方法都無(wú)法對(duì)崩壁地形進(jìn)行局部監(jiān)測(cè)，且精度很低；地面近景攝影測(cè)量需要永久固定攝影位置而在侵蝕劇烈的崩崗地貌中很難完成且只能固定局部范圍進(jìn)行監(jiān)測(cè)；航空遙感監(jiān)測(cè)周期長(zhǎng)適用的尺度很大，而崩崗侵蝕的面積小而且分散，難以滿足精細(xì)的研究要求；GPS 測(cè)量技術(shù)需要人工到達(dá)很多危險(xiǎn)的地形，且在崩崗中測(cè)量衛(wèi)星信號(hào)難以保證精度。三維激光掃描作為近年來(lái)的創(chuàng)新測(cè)量技術(shù)，具有高精度，穿透性，掃描速度快等特點(diǎn)能很準(zhǔn)確的獲得點(diǎn)云數(shù)據(jù)、但是儀器價(jià)格高昂不方便攜帶，在崩崗中有很多深溝使其會(huì)有很多掃描死角，且多站點(diǎn)的視角在拼接和配準(zhǔn)時(shí)不易操作。

本文使用的無(wú)人機(jī)傾斜航空攝影測(cè)量方法，不僅價(jià)格低廉而且空中視角可以全方位的獲取崩崗地形點(diǎn)云數(shù)據(jù)，獲取數(shù)據(jù)速度更快，在崩崗這種小尺度的侵蝕地貌監(jiān)測(cè)中會(huì)有廣闊的應(yīng)用前景和實(shí)際的使用價(jià)值。但是，在操作中對(duì)控制點(diǎn)的要求很高否則很難生成高精度的點(diǎn)云數(shù)據(jù)，且飛行時(shí)需要手動(dòng)操控拍照容易出現(xiàn)影像數(shù)據(jù)獲取質(zhì)量不高的情況，目前只以單體面積較小的崩崗作為研究對(duì)象。在點(diǎn)云的處理方法上，本文使用的方法較為復(fù)雜但是能夠很精細(xì)的測(cè)算出崩崗局部的侵蝕情況，然而如果使用在大型崩崗上時(shí)就會(huì)顯得繁瑣復(fù)雜很難定義侵蝕與沉積，所以如何實(shí)踐對(duì)結(jié)構(gòu)更為復(fù)雜的大型崩崗的三維點(diǎn)云建模并提高侵蝕量測(cè)算精細(xì)度，有待進(jìn)一步研究。

4 結(jié) 論

本文將基于無(wú)人機(jī)傾斜航空影像建立三維點(diǎn)云模型測(cè)算侵蝕量的方法運(yùn)用在了崩崗侵蝕方面并確定了其適用性和點(diǎn)云數(shù)據(jù)處理方法。該方法以消費(fèi)級(jí)多旋翼無(wú)人機(jī)搭載普通數(shù)碼相機(jī)獲取目標(biāo)崩崗的多角度傾斜航空影像為基礎(chǔ)，根據(jù)傾斜攝影測(cè)量原理對(duì)影像進(jìn)行空三加密處理，生成點(diǎn)云模型的平均總均方根誤差僅為0.042 m 達(dá)到厘米級(jí)測(cè)量，再使用Arcgis 提取崩崗地形信息，基于點(diǎn)云模型運(yùn)用多時(shí)相DEM 差分法和點(diǎn)云2.5D 體積算法兩種方法配合測(cè)算崩崗侵蝕量，監(jiān)測(cè)崩崗不同部位的侵蝕量動(dòng)態(tài)變化，并對(duì)結(jié)果進(jìn)行精度檢驗(yàn)。4 個(gè)周期的崩崗侵蝕量計(jì)算的平均相對(duì)誤差為9.69%，但在崩崗這個(gè)尺度的侵蝕地貌中絕對(duì)誤差最大值僅為0.3033 m3。利用無(wú)人機(jī)傾斜航空攝影監(jiān)測(cè)崩崗侵蝕沉積量的動(dòng)態(tài)變化是有效且便捷的方法。