——以三峽庫區(qū)張家灣滑坡為例"/>
杜磊, 陳潔,2, 李敏敏, 鄭雄偉, 李京, 高子弘
(1.中國自然資源航空物探遙感中心,北京 100083; 2.中國科學(xué)院遙感與數(shù)字地球研究所,北京 100101)
滑坡作為一種全球范圍內(nèi)普遍存在而且頻繁發(fā)生的地質(zhì)災(zāi)害現(xiàn)象,對其規(guī)律性進(jìn)行研究和預(yù)警一直是地學(xué)科研熱點之一。目前滑坡調(diào)查的手段主要有野外實地勘查和遙感調(diào)查2種,前者作為傳統(tǒng)的工作方法,可近距離現(xiàn)場研究滑坡形態(tài),圈定滑坡要素,但其工作耗時久、效率低,抵達(dá)調(diào)查區(qū)域也受限較多; 衛(wèi)星遙感技術(shù)支持下的滑坡調(diào)查能快速獲取大尺度、多時相的數(shù)字遙感數(shù)據(jù),現(xiàn)已形成一套較為完善的“數(shù)字滑坡”識別技術(shù)[1]。但該技術(shù)的不足之處在于受天氣影響較大、且對被植被覆蓋的滑坡識別能力有限,源數(shù)據(jù)生成的數(shù)字高程模型(digital elevation model,DEM)精度較低,難以刻畫地表微小變形的精細(xì)特征。
機載激光雷達(dá)(light detection and ranging,LiDAR)作為一種主動式航空遙感技術(shù),理論上可不受天氣限制,其發(fā)射的激光點云能夠穿透一定厚度的植被并獲取地表信息。與傳統(tǒng)航空攝影相比,機載LiDAR數(shù)據(jù)采集方式能在很大程度上減弱地形切割陰影的影響,其構(gòu)建的DEM垂直精度可達(dá)10 cm,可以精確表達(dá)微地貌地形特征[2],具備提升滑坡體識別精度的條件。近年來,美國地質(zhì)調(diào)查局及歐洲太空局等機構(gòu)均開展了基于機載LiDAR衍生DEM的滑坡識別研究,取得了一定的進(jìn)展[3-4],包括: 基于高精度DEM計算滑坡坡度、坡向、地表粗糙度、半方差和分形等特征; 結(jié)合傅里葉變換、數(shù)據(jù)小波分析等多種算法開展地表地形有關(guān)參數(shù)精細(xì)特征研究[5],并探索建立滑坡等地質(zhì)災(zāi)害分布與上述有關(guān)因子的聯(lián)系。以上研究充分顯示了機載LiDAR技術(shù)在滑坡調(diào)查中的應(yīng)用潛力。
本文基于機載LiDAR技術(shù)獲取了研究區(qū)點云數(shù)據(jù),通過對點云數(shù)據(jù)先后開展檢校、濾波和分類等數(shù)據(jù)處理,構(gòu)建研究區(qū)精細(xì)的DEM產(chǎn)品,再利用這些DEM產(chǎn)品對研究區(qū)的滑坡體進(jìn)行識別以及重要參數(shù)提取,取得了較好的應(yīng)用效果。
三峽庫區(qū)一直以來是我國滑坡等地質(zhì)災(zāi)害的頻發(fā)區(qū)和重災(zāi)區(qū)。張家灣滑坡是庫區(qū)較為典型的一處滑坡群。它由6個傾向長江的單體滑坡組成,位于童莊河右岸,隸屬湖北省宜昌市秭歸縣郭家壩鎮(zhèn)頭道河村一組,與入江口相距約1.5 km?;麦w總體前緣高程約為120 m,后緣高程約290 m,水位以上總面積約29×104m2,對應(yīng)總體積約550×104m3。滑坡區(qū)地勢整體東高西低,后側(cè)山頂高程約500 m,滑坡前緣為三峽庫區(qū)水面,中部及前緣地帶坡面相對平緩,坡度為15°~25°,后緣及后側(cè)山坡較陡峭,坡度為25°~35°?;聟^(qū)內(nèi)出露的地層主要為侏羅系砂頁巖及三疊系中統(tǒng)巴東組(T2b)砂巖、泥灰?guī)r等半堅硬巖類組成。
選擇2009年由中國自然資源航空物探遙感中心獲取的巴東—秭歸地質(zhì)災(zāi)害易發(fā)多發(fā)重點區(qū)1 000 km2LiDAR數(shù)據(jù)作為遙感數(shù)據(jù)源,其真實記錄了三峽庫區(qū)自移民遷建、175 m水位蓄水和2008年汶川震后的壩體微變化,具有開展庫區(qū)庫岸穩(wěn)定性、地質(zhì)災(zāi)害和生態(tài)地質(zhì)環(huán)境變化調(diào)查等研究的潛力。
針對現(xiàn)有滑坡調(diào)查技術(shù)在復(fù)雜地質(zhì)條件下難以識別滑坡,以及傳統(tǒng)基于DEM數(shù)據(jù)識別滑坡方法存在的問題,本文提出了一種從機載LiDAR構(gòu)建的精細(xì)DEM中表征滑坡地表特征參數(shù)的技術(shù)方法,實現(xiàn)了滑坡邊界的自動識別和相應(yīng)參數(shù)提取。新方法步驟如下:
1)針對目標(biāo)區(qū)域制定作業(yè)計劃,使用定位定向系統(tǒng)(position orientation system,POS)輔助機載LiDAR系統(tǒng)對目標(biāo)區(qū)域開展航空遙感數(shù)據(jù)采集。
2)對飛行獲取的機載LiDAR點云數(shù)據(jù)先后進(jìn)行檢校、濾波和分類等數(shù)據(jù)處理,得到消除植被覆蓋信息的純地表精細(xì)DEM。
3)從DEM中分析、提取出地貌特征并建立紋理特征參數(shù),形成特征參數(shù)文件。
4)利用計算機智能分類算法計算目標(biāo)區(qū)域中已知滑坡像元以及非滑坡像元特征參數(shù),結(jié)合上一步獲取的特征參數(shù)文件,確定最優(yōu)特征參數(shù)。
5)選取目標(biāo)區(qū)域內(nèi)的部分已知滑坡像元和非滑坡像元作為算法訓(xùn)練集,不斷優(yōu)化訓(xùn)練集各元素,綜合確定最優(yōu)特征參數(shù)組合。利用智能分類算法預(yù)測訓(xùn)練集中的像元屬性同時進(jìn)行精度評價,獲得符合預(yù)設(shè)精度條件的平衡系數(shù)。
6)利用獲取平衡系數(shù)的訓(xùn)練集以及最優(yōu)特征參數(shù)組合,訓(xùn)練分類模型,預(yù)測已知滑坡像元和非滑坡像元的數(shù)據(jù)集,并分別計算數(shù)據(jù)集的平均用戶精度、平均生產(chǎn)者精度和總體精度。
7)若所計算的精度滿足設(shè)定閾值要求,則使用滿足該平衡系數(shù)和最優(yōu)特征參數(shù)組合的智能分類模型預(yù)測整個研究區(qū),再采用邊緣檢測算子計算滑坡邊界,最終實現(xiàn)滑坡識別。
8)對識別出的滑坡像元集,通過對應(yīng)的精細(xì)DEM提取出滑坡坡度和地表粗糙度等滑坡參數(shù),并完成定量分析。
2009年4—6月,使用運5-B型飛機作為航空遙感平臺,搭載Leica Geosystem公司制造的ALS50-II機載LiDAR系統(tǒng),獲取了包含張家灣滑坡在內(nèi)約1 000 km2研究區(qū)綜合遙感數(shù)據(jù)資料,包括機載LiDAR點云數(shù)據(jù)、數(shù)字航空影像、機載GPS數(shù)據(jù)、機載IMU數(shù)據(jù)和地面GPS參考站定位數(shù)據(jù)等。已獲取的機載LiDAR點云數(shù)據(jù)包含激光距離、角度和多次回波強度等重要信息。通過對獲取的原始數(shù)據(jù)進(jìn)行數(shù)據(jù)解算、系統(tǒng)誤差檢校等預(yù)處理,得到了含地物三維坐標(biāo)信息的地表激光點云數(shù)據(jù)。
此次ALS50-II機載LiDAR系統(tǒng)飛行相對航線高度約為1 000 m,航線旁向重疊度為30%~39%,飛行地速為150~180 km/h,機載激光空中掃描角度為45°,其最大掃描頻率設(shè)置為45 Hz,機載LiDAR平均點云密度為3.2個/m2。
受系統(tǒng)誤差和偶然誤差影響,機載LiDAR系統(tǒng)獲取的原始LiDAR點云坐標(biāo)存在系統(tǒng)偏差,需先通過數(shù)據(jù)檢校消除誤差。機載LiDAR點云數(shù)據(jù)檢校一般指檢校側(cè)滾角、俯仰角和航偏角3個安置角和測距誤差。
3個安置角的檢校方法如下: 首先,對檢校場的數(shù)據(jù)開展人工初檢校,分別獲取較為準(zhǔn)確的3個安置角檢校參數(shù); 然后,對初檢校的機載LiDAR數(shù)據(jù)開展進(jìn)一步精檢校,獲取檢校參數(shù)的改正值; 最后,將上述2次檢校參數(shù)相加得到安置角最終的檢校參數(shù)。
完成安置角檢校后,利用實測GPS點對測距誤差進(jìn)行檢校,獲取測距檢校參數(shù)。利用上述安置角和測距檢校參數(shù),對研究區(qū)全部機載LiDAR點云數(shù)據(jù)進(jìn)行檢校,以消除研究區(qū)點云數(shù)據(jù)系統(tǒng)誤差。經(jīng)驗證,研究區(qū)域各航帶點云匹配效果較好,高程方向相對高差優(yōu)于10 cm,平面位置相對差優(yōu)于25 cm,結(jié)果穩(wěn)定可靠。
對于數(shù)據(jù)檢校解算依然無法較好消除的殘余航帶系統(tǒng)誤差,本文采用目前已相對較為成熟的Burman航帶平差方法進(jìn)行改善,通過TerraSolid商業(yè)軟件系列中的TerraMatch模塊實現(xiàn)。
生成DEM前,先對機載LiDAR點云數(shù)據(jù)進(jìn)行濾波分類處理,以區(qū)分地面點和非地面點,再利用地面點數(shù)據(jù)制作所需的高精度DEM。LiDAR點云分類處理采用了美國Bentley公司TerraSolid軟件中TerraScan和TerraModeler等模塊。其中,TerraScan模塊地面點分類是基于Axelsson提出的不規(guī)則三角格網(wǎng)(triangulated irregular network,TIN)加密改進(jìn)算法實現(xiàn)[6]。TerraScan模塊提供了地面點分類的Ground工具命令,該方法利用迭代算法建立TIN表面模型來區(qū)分地面點和非地面點。
開展數(shù)據(jù)迭代運算前需首先選定若干初始地面點,通過Max building size參數(shù)調(diào)節(jié)初始點的選擇。實驗表明: 當(dāng)Max building size取值30 m時,不僅可以很好地去除植被和房屋并保留微地形特征,而且還可以減少手工編輯的工作量。
通過分類前后對比(圖1)可以看出,經(jīng)濾波分類后的DEM中大部分建筑物和植被信息被自動濾除。少量建筑物和植被殘余信息則用人工編輯的方法予以剔除。
(a) 分類前(b) 分類后
圖1分類前后DEM渲染對比
Fig.1ComparisonofDEMshadingmapsbeforeandaftercloud-pointclassification
先利用機載LiDAR系統(tǒng)獲取的海量LiDAR點云數(shù)據(jù)快速開展目標(biāo)區(qū)域地面三維數(shù)字地形模型(digital terrain model,DTM)采集,再對原始點云數(shù)據(jù)進(jìn)行系統(tǒng)誤差糾正和點云數(shù)據(jù)分類,有效濾除了點云中的非地面點,進(jìn)而獲得目標(biāo)區(qū)域高精度DEM。為了評價DEM產(chǎn)品的精度,在研究區(qū)內(nèi)均勻布設(shè)30個地面檢查點,通過DEM讀取記錄地面檢查點的三維坐標(biāo)后,再和實測同名點三維坐標(biāo)進(jìn)行對比計算,完成研究區(qū)DEM的精度質(zhì)量分析(表1)。結(jié)果表明,生成的DEM高程中誤差滿足1: 1 000比例尺規(guī)范精度要求[7],可應(yīng)用于滑坡識別調(diào)查。
表1 DEM成果精度評價Tab.1 DEM precision evaluation (m)
滑坡的解譯標(biāo)志主要為“箕”狀色調(diào)帶和紋理異常帶?;麦w一般位于地質(zhì)較穩(wěn)定的自然斜坡凸突的負(fù)地形中,滑坡后壁與滑體交接處往往形成洼地,中部則有多級垂直滑動方向的臺坎。圖2是張家灣滑坡的數(shù)字正射影像圖(digital orthophoto map, DOM),它是通過機載LiDAR同步獲取的航空遙感影像制作而成。由圖2可以看出,張家灣滑坡是一傾向長江的單斜順向滑坡群,影像中滑體一般呈淺色調(diào),通常具有一定的擠壓、擾動或松脫等特征,巖(土)體破碎,結(jié)構(gòu)松散,耕地、園地及林地被擾動或者破壞比較嚴(yán)重,滑體紋理與背景呈突變關(guān)系。
圖2 張家灣滑坡區(qū)域數(shù)字正射影像圖Fig.2 DOM of Zhangjiawan Village landslide area
本文的野外工作包括前期野外踏勘和解譯成果野外驗證2部分。前期野外踏勘主要是對研究內(nèi)容進(jìn)行地面初步感性認(rèn)識,結(jié)合相應(yīng)遙感影像特征,建立遙感解譯標(biāo)志; 解譯成果野外驗證則是對室內(nèi)目視解譯結(jié)果在實地開展查驗,并對解譯標(biāo)志做進(jìn)一步修正,從而達(dá)到理性認(rèn)識,使遙感解譯成果與實際情況更加吻合。
通過機載LiDAR點云數(shù)據(jù)制作的張家灣滑坡區(qū)域DEM可精細(xì)表達(dá)該地區(qū)微地形地貌信息,本文通過對其進(jìn)行定性及定量分析完成了滑坡調(diào)查有關(guān)應(yīng)用研究。
定性分析是通過張家灣DEM制作不同視角山體陰影圖完成。山體陰影圖是通過模擬光線計算DEM每個柵格單元的照明值。通過對比分析不同角度山體陰影圖,不僅可較好地掌握張家灣滑坡區(qū)域的立體形態(tài),還可以有效地提取出某些地形遮蔽信息,如滑坡、斷層崖等線性特征。制作張家灣滑坡區(qū)域山體陰影圖需確定3個重要參數(shù),分別是: 太陽方位角、太陽高度角和表面灰度值。考慮到研究區(qū)內(nèi)干流及支流走向,選擇以相同的太陽高度角(45°),用180°,135°和90°這3個不同的方位角,制作山體陰影圖,表面灰度值設(shè)定為0~255。
定量方法則是指通過目標(biāo)區(qū)域的DEM數(shù)據(jù)提取精細(xì)微地形地貌參數(shù),分析滑坡坡度、地表粗糙度、半方差和分維等滑坡要素信息。其中,地表粗糙度是指在某確定范圍內(nèi)地表表面積與其在水平方向上的投影面積之比,是反映地形起伏和被侵蝕程度的一個重要指標(biāo),可有效反映地表形態(tài)的宏觀變化。半方差是指滑坡體相對于周圍地形表現(xiàn)在DEM數(shù)據(jù)值上的突變,它反映地表性質(zhì)的不同距離觀測值之間的變化。當(dāng)DEM設(shè)置為半方差變量時,半方差函數(shù)則能夠衡量地形的空間自相關(guān)性。相同尺度下,半方差值越大,對應(yīng)空間地形變化越大,滑坡體與非滑坡體識別則越容易。分形是指部分與整體以某種形式相似的形; 分維是指分形的維數(shù),它不僅可以衡量地形粗糙度,還可以體現(xiàn)地形的復(fù)雜度和空間自相關(guān)性。相同條件下,目標(biāo)區(qū)域分維數(shù)值變化越大,其對應(yīng)地形越復(fù)雜,地表則越粗糙[3]。
通過張家灣地區(qū)高精度DEM制作不同視角下的山體陰影圖,結(jié)合張家灣滑坡的坡度圖和地表粗糙度圖完成滑坡識別; 通過計算張家灣滑坡體半方差函數(shù)和分維完成對滑坡的定量分析。圖3(a)是張家灣滑坡野外實地勘查照片,紅線是其中2處滑坡的邊界線; 圖3(b)是張家灣滑坡DEM彩色暈渲圖,通過對其進(jìn)行數(shù)值分析可了解滑坡地貌形態(tài)特征。
(a) 野外實地勘查照片(b) 滑坡DEM暈渲圖
圖3張家灣滑坡野外實地勘查照片及其滑坡DEM暈渲圖
Fig.3LandslidefieldphotoanditsDEMshadingmapofZhangjiawanVillage
圖4是張家灣滑坡2009年山體陰影系列圖,圖4(a),(b),(c)分別是太陽方位角設(shè)置為90°,135°和180°時制作出的山體陰影圖(太陽高度角均設(shè)置為45°,表面灰度值設(shè)定為0~255)。對比分析3幅山體陰影圖可知: ①與DEM數(shù)據(jù)柵格圖像相比,3幅山體陰影圖均可較好地反映滑坡地形的立體形態(tài),對沖溝等線性地物特征的表現(xiàn)效果尤為突出; ②對滑坡體局部地形地貌特征的表現(xiàn),3幅山體陰影圖則差異明顯,比如: 滑坡沖溝等線性地物在太陽方位角為90°和180°的山體陰影圖上明顯比135°的山體陰影圖要清晰可見; ③太陽方位角為90°和180°的山體陰影圖反映出的局部地貌特征相似但又不完全相同,但2幅影像信息呈互補關(guān)系。綜上分析,通過高精度DEM數(shù)據(jù)制作的不同太陽方位角山體陰影系列圖,不僅能有效地識別出張家灣6個單體滑坡的滑動范圍,還可較為準(zhǔn)確地確定其滑坡后緣、側(cè)緣和滑舌等滑坡要素信息。
(a) 方位角為90°(b) 方位角為135°(c) 方位角為180°
圖4張家灣滑坡山體陰影系列圖
Fig.4SlidemountainshadowseriesmapsofZhangjiawanVillage
圖5是張家灣滑坡的坡度圖和地表粗糙度圖。由張家灣滑坡坡度圖分析可知: ①滑坡側(cè)緣或滑坡兩側(cè)坡度跳躍變化明顯,且跳躍變化幅度大體一致,可在坡度圖圈定出滑坡側(cè)緣邊界(圖5(a)中紅線所示); ②滑坡壁與滑舌的坡度值較為集中,兩者色彩較均勻,且在形態(tài)上表現(xiàn)為面狀分布,亦可較好地將6個不同的滑坡體區(qū)分開; ③在滑坡體與滑舌、滑坡壁與滑坡體的分界處,滑坡側(cè)緣及其附近地形地貌處,坡度變化明顯并出現(xiàn)跳躍式變大或縮小,且跳躍變化幅度大致相同。綜上,通過滑坡坡度圖也可以有效圈定出張家灣滑坡的3大要素。此外,結(jié)合圖4和圖5可知,在張家灣滑坡群上方,地形輪廓清晰且呈封閉狀,“雙溝同源”現(xiàn)象比較明顯,表明此處具有誘發(fā)滑坡地質(zhì)災(zāi)害的地形特征。結(jié)合山體陰影系列圖,可圈出滑坡隱患區(qū)域(圖5中藍(lán)線所示區(qū)域)。
(a) 滑坡坡度圖(b) 地表粗糙度圖
圖5張家灣滑坡坡度圖和地表粗糙度圖
Fig.5LandslideslopemapandlandslidesurfaceroughnessmapofZhangjiawanVillage
對張家灣地區(qū)地表粗糙度圖進(jìn)行解譯,結(jié)果顯示: ①滑坡側(cè)緣地表粗糙度與兩側(cè)地表粗糙度存在明顯跳躍,地表粗糙度跳躍幅度也基本相同,據(jù)此從地表粗糙度圖上可圈定出張家灣地區(qū)6個不同滑坡體滑坡側(cè)緣的邊界(圖5(b)中紅線所示); ②對于張家灣地區(qū)滑坡圖像灰階分布,其滑坡體后部與滑坡壁的地表粗糙度值較集中,色彩也較勻滑,形態(tài)上表現(xiàn)為面狀分布,據(jù)此亦可以將它們與滑坡體其他部分區(qū)別開; ③在滑坡壁與滑坡體、滑坡側(cè)緣與其附近地貌處,地表粗糙度明顯出現(xiàn)跳躍式變大或縮小,且跳躍變化幅度大致相同。綜上,通過地表粗糙度圖亦可有效圈定出張家灣滑坡的3大要素。
分別選擇滑坡后緣、滑坡體和滑舌3部分作為半方差函數(shù)計算范圍。張家灣地區(qū)滑坡后緣、滑坡體和滑舌的半方差曲線見圖6。
圖6 張家灣滑坡要素半方差曲線Fig.6 Landslide elements semivariogram curve of Zhangjiawan Village
由圖6可知: ①滑坡后緣、滑坡體和滑舌半方差曲線斜率存在明顯區(qū)別,這是滑坡要素識別的數(shù)學(xué)理論基礎(chǔ); ②滑坡后緣、滑坡體和滑舌的半方差函數(shù)曲線隨變程增加均變大,表明張家灣地區(qū)各滑坡要素的地形地貌相關(guān)性較強; ③滑坡后緣半方差值最大,表明滑坡后緣地形空間變化最大; 滑坡體半方差值次之,表明滑坡體地形空間變化幅度介于滑坡后緣與滑舌之間; 滑舌半方差值最小,表明滑舌部分對應(yīng)的空間地形變化也最小。雖然3種滑坡要素的半方差函數(shù)曲線均隨著變程的增大而增大,但也存在明顯區(qū)別: 在相同變程條件下,滑坡后緣值最大,滑坡體次之,而滑舌最小。
張家灣3種滑坡要素分維函數(shù)的計算范圍與半方差相同,其分維曲線與同區(qū)域2006年分維曲線對比可知: 滑坡后緣的分維值變化最大,表明滑坡后緣對應(yīng)地形空間的變化最大; 滑坡體的分維值變化次之,表明滑坡體對應(yīng)地形空間的變化幅度介于滑坡后緣與滑舌之間; 滑舌的分維值變化最小,說明滑舌空間地形變化最小(如圖7所示)。滑坡要素分維值的計算結(jié)果與半方差計算結(jié)果具有較好的一致性。
(a) 滑坡后緣分維變化(b) 滑坡體分維變化(c) 滑舌分維變化
圖7張家灣滑坡要素分維曲線對比
Fig.7ComparisonoflandslidefractaldimensioncalculationresultsofZhangjiawanVillage
機載激光雷達(dá)技術(shù)作為一種新的航空遙感調(diào)查與監(jiān)測手段,在地質(zhì)災(zāi)害調(diào)查與監(jiān)測領(lǐng)域已獲得較廣泛關(guān)注和應(yīng)用。與衛(wèi)星和航空影像相比,機載激光雷達(dá)具有一定的植被穿透能力,多次回波可以用來精確測量和表達(dá)微地貌特征信息,在霧氣較大或陰雨天氣時也能快速獲取真實地形信息,具備提升滑坡體識別精度的條件。具體表現(xiàn)在2方面: ①與邊坡的整體變化相關(guān)的特征,如存在可以被機載激光雷達(dá)點云數(shù)據(jù)派生的DEM識別的斜坡凹陷和突變; ②結(jié)合機載激光雷達(dá)點云數(shù)據(jù)派生的精細(xì)DEM,結(jié)合地表粗糙度變化,可以有效識別某些特殊地表特征,如地表內(nèi)部結(jié)構(gòu)異常、地裂隙或滑痕等。
本文基于機載激光雷達(dá)系統(tǒng)獲取的航空遙感數(shù)據(jù)構(gòu)建了精細(xì)的DEM,通過對DEM產(chǎn)品進(jìn)行深入挖掘,進(jìn)行綜合對比分析,參考有關(guān)地質(zhì)資料,結(jié)合一些必要的野外驗證,有效地獲得了三峽庫區(qū)張家灣滑坡群的滑坡信息,應(yīng)用效果良好,表明采用機載激光雷達(dá)技術(shù)對地質(zhì)災(zāi)害進(jìn)行調(diào)查與監(jiān)測是一種直觀、快速而經(jīng)濟(jì)的方法,尤其對于災(zāi)害應(yīng)急反應(yīng)和災(zāi)后救援具有重要意義。