庾 露,黃 艷 霞
(1.南寧師范大學 地理科學與規(guī)劃學院,廣西 南寧 530100; 2.廣西國土資源規(guī)劃設計集團有限公司,廣西 南寧 530029)
滑坡是常見的自然災害之一,極易誘發(fā)次生災害,對人類的生存環(huán)境造成極大的威脅,一直是災害防治工作中的重中之重[1]。據(jù)自然資源部統(tǒng)計,2020年全國共發(fā)生地質(zhì)災害7 840起,其中滑坡4 810起,占總數(shù)的61%,滑坡災害呈常態(tài)化趨勢[2]。獲取滑坡的體積是研究滑坡災害的一個重要內(nèi)容,滑坡體積對滑坡的破壞力、展布范圍及其可能引發(fā)的二次災害如堰塞湖和洪水具有指示意義[3-5],精確的滑坡測量能降低滑坡相關的危險度和風險性分析中的不確定因素并且有助于進一步的防災減災[6]。
數(shù)字表面模型(Digital Surface Model,DSM)差值計算是測量滑坡體積方法之一,通常使用遙感測繪技術獲取滑坡前后的高程分布,利用高程差和像元面積的乘積獲取滑坡體積,具有覆蓋范圍較大、時間和人力成本較低的優(yōu)點。新一代遙感衛(wèi)星和無人機平臺的成熟并大規(guī)模投入使用,以及LiDAR、InSAR等測量技術的發(fā)展,為地質(zhì)災害監(jiān)測提供了多樣的數(shù)據(jù)來源和分析手段[7]。如MA[8],彭大雷[9],CHANG[10]等使用低空無人機測量計算出滑坡體積,此方法基本不受云層遮擋,但需要攜帶設備到達現(xiàn)場開展工作。同時,受飛行平臺穩(wěn)定性和載荷量限制,無人機遙感系統(tǒng)所獲取的圖像往往像幅較小、數(shù)量較多、畸變較大,圖像處理難度較高[11]。衛(wèi)星遙感具有覆蓋范圍大、成本低、精度高的特點,也可作為滑坡體積估算的手段[12],但由于滑坡發(fā)生后,局部通常伴隨長時間降雨,會對衛(wèi)星觀測造成諸多不利影響。傳統(tǒng)的光學立體測繪衛(wèi)星受固有重訪周期限制、拍攝調(diào)度安排和降雨天氣云層的遮擋,往往在單顆衛(wèi)星過境的短暫間隙內(nèi)無法獲取有效影像。星載重復軌道InSAR測量在成像過程中,大氣水汽含量的時空變化會引起雷達信號傳播延遲,給干涉相位帶來較大的附加相位變化,導致干涉測量結果出現(xiàn)誤差[13]。因此,有必要根據(jù)研究區(qū)所收集的影像數(shù)據(jù)和地形環(huán)境情況,進一步拓展滑坡地形測量的手段。
基于光學衛(wèi)星的異軌立體測量技術,是對傳統(tǒng)光學測繪衛(wèi)星同軌測量方法的有效補充。該技術可將未設計用于立體測繪的衛(wèi)星數(shù)據(jù)充分利用,借助單顆或多顆衛(wèi)星重復過境縮短重訪周期的優(yōu)勢,多期覆蓋研究區(qū),將可能存在的有效影像篩選出來并組成異軌立體像對,生成DSM產(chǎn)品。已有學者使用該技術提取了地物高程,如胡芬等[14]對比分析了異軌和同軌立體影像生成的DSM,發(fā)現(xiàn)二者精度趨近。倪文儉等[15]使用高分2號衛(wèi)星(GF-2)數(shù)據(jù)形成異軌立體像對,提取森林高度,與LiDAR結果的均方根誤差為3.6 m。
上述研究成果驗證了該技術的可行性,但目前的研究仍主要集中于單一成像時間點上的DSM提取,尚缺乏考慮滑坡事件發(fā)生前后高程劇烈變化產(chǎn)生的數(shù)據(jù)一致性問題,以及在地形起伏大、無實地控制點約束下的精度控制情況。因此,有必要就如何將異軌立體測量技術應用于滑坡體積提取場景做更深入的研究。本文以2020年7月11日湖北省恩施市馬者村沙子壩滑坡為例,重點分析使用高分2號衛(wèi)星(GF-2)通過異軌立體成像立體測量技術在無實地控制點的場景下,提取滑坡后DSM的可行性和精度。同時,結合資源3號02星(ZY-302)影像提取的滑坡前DSM,估算滑坡體積,并對滑坡的構成進行分析。
受梅雨期連續(xù)降雨影響,2020年7月20~21 日湖北省恩施市清江上游馬者村沙子壩(東經(jīng)109°18′,北緯30°22′)滑坡垮塌,發(fā)生縱向長1 200~1 500 m,橫向?qū)?20~580 m的特大型滑坡,導致清江干流河道形成堰塞湖,水位的急速上升導致存在潰壩風險[16]。研究區(qū)地處巫山余脈和武陵山北上的余支交會部,地勢整體為西低東高,最大落差達1 380 m,地形崎嶇復雜,多為海拔 1 000 m 以上的高山。在氣候上,研究區(qū)屬亞熱帶大陸性季風氣候,雨量充沛,易發(fā)生暴雨洪澇,造成滑坡。具體的影像范圍與研究區(qū)范圍如圖1所示。
圖1 研究區(qū)域位置和影像范圍Fig.1 Location and image coverage of the study area
滑坡前的DSM數(shù)據(jù)由2020年3月19日的資源三號02(ZY-302)衛(wèi)星影像生成,這是距離發(fā)生滑坡時間較近、成像質(zhì)量較高的衛(wèi)星影像數(shù)據(jù)。ZY-3和后續(xù)發(fā)射的ZY-302星均是中國民用高分辨率立體測繪衛(wèi)星,通過星載三線陣光學相機,以同軌方式實現(xiàn)地形圖繪制?;潞蟾采w研究區(qū)域的可用影像較少,以ZY-3和ZY-302為例,在滑坡發(fā)生后180 d內(nèi)僅有3景影像覆蓋,且滑坡區(qū)受到云層遮擋。而在此期間,高分二號(GF-2)衛(wèi)星以不同角度拍攝了多景覆蓋滑坡區(qū)的影像,雖然該衛(wèi)星并不是為立體測繪任務而設計的,但具有最高0.8 m的高空間分辨率和最大35°的側(cè)擺能力,在不同觀測角度和較短重復觀測周期的前提下可構成異軌立體觀測像對,從而實現(xiàn)DSM提取[15]。經(jīng)過反復試驗,本文選取了距離滑坡發(fā)生后時間較接近的2020年11月08日和2020年12月22日的兩幅GF-2影像構成立體像對,用以提取滑坡后DSM。
本文使用ESRI公司提供的全球衛(wèi)星影像網(wǎng)絡服務作為平面位置參考,使用GLO30 DSM數(shù)據(jù)作為高程參考。前者在研究區(qū)范圍內(nèi)提供了分辨率為1 m,經(jīng)過地形校正和拼接處理的衛(wèi)星影像底圖,其分辨率與GF-2全色波段影像的分辨率接近,并優(yōu)于ZY-302影像的分辨率,可用于生成控制點和驗證點的平面坐標;后者為歐空局哥白尼地球觀測計劃的數(shù)字表面高程產(chǎn)品,數(shù)據(jù)來源于TanDEM-X任務期間獲得的SAR衛(wèi)星資料,GLO30的像元尺寸為30 m,絕對垂直精度<4 m,絕對水平精度<6 m[17],研究區(qū)最新數(shù)據(jù)更新至2014年可用于生成控制點和驗證點的高程值。
2.2.1技術流程
對于滑坡區(qū)DSM提取問題,受其所處地理位置、地形地貌和自身不穩(wěn)定地質(zhì)環(huán)境等多方面因素制約,可能不利于開展實測工作。為此,本文提出完全基于衛(wèi)星遙感數(shù)據(jù)源進行異軌DSM提取和精度校正的新思路,并在后文的論述中討論該方法的可行性。本次研究的技術流程如圖2所示,使用1景ZY-3 02衛(wèi)星影像提取滑坡前DSM,使用2景異軌的GF-2衛(wèi)星影像提取滑坡后DSM,使用ESRI衛(wèi)星全球衛(wèi)星影像作為平面位置參考,GLO30作為高程參考對提取的DSM進行位置和高程校正和精度評價,最后計算滑坡體積。
圖2 技術路線Fig.2 Technical route
2.2.2異軌立體測量基本原理
異軌立體測量是基于光學衛(wèi)星所具備的側(cè)擺功能,由單臺相機在不同軌道觀測獲取立體影像的方法。如圖3所示,在兩個不同軌道上分別對目標成像,產(chǎn)生類似于雙目立體視覺,通過獲取兩個軌道的基線、焦距以及視差等信息可計算出地物的三維點云分布。異軌立體測量的精度包括平面精度和高程精度。其中平面精度和同軌立體測量方式接近[18-20],可通過衛(wèi)星數(shù)據(jù)供應商提供的RPC(Rational Polynomial Cofficient)參數(shù),運用成像幾何模型,并配合一定數(shù)量的地面控制點建立圖像坐標和地理坐標的幾何關系。而高程精度則主要受飛行平臺的軌道精度、影像空間分辨率、成像紋理清晰度、基高比等影響。對于異軌立體的測量,通常前3個因素基本相同或接近,而基高比可根據(jù)所收集的影像進行較靈活的調(diào)整。理論上基高比越大,交會條件越好,高程測量的精度越高。根據(jù)公式(1)基高比B/H由兩軌影像的側(cè)視角α和β決定:
圖3 基高比示意Fig.3 Baseline-to-height ratio schematic
(1)
式中:B為軌道間基線長度,H為攝站高度,α和β分別為側(cè)視角。經(jīng)計算,GF-2異軌像對的基高比約為0.50,ZY-302衛(wèi)星立體像對的基高比約為0.83。
2.2.3利用GF-2提取DSM
由于本次滑坡體積提取實驗未到達現(xiàn)場采集地面控制點,因此用于高精度幾何定位的控制點和用于DSM結果精度評估的驗證點均通過已公開的免費遙感數(shù)據(jù)產(chǎn)品采集獲得。本文使用ArcGIS Desktop載入全球衛(wèi)星影像網(wǎng)絡服務和GLO30 DSM數(shù)據(jù),采集地形起伏較小、地表紋理清晰穩(wěn)定、易于辨認的人造地物(如道路交叉口),作為控制點和驗證點。使用Trimble Inpho軟件,導入控制點,對GF-2衛(wèi)星異軌立體影像采用區(qū)域網(wǎng)平差的方式進行高精度的幾何定位,其流程如圖4所示。 首先,通過特征匹配方法自動提取大量均勻分布的連接點,實現(xiàn)立體影像密集匹配,提高匹配的效率和可靠性,降低誤匹配率;通過自由網(wǎng)平差定向,對觀測像點的殘差進行統(tǒng)計分析,進而采用選權迭代法實現(xiàn)誤匹配連接點的自動識別與剔除。然后,進行控制點半自動量測,通過區(qū)域網(wǎng)平差定向和定向殘差分析,迭代計算得到區(qū)域網(wǎng)平差定向結果。最后,在自由網(wǎng)平差或區(qū)域網(wǎng)平差的基礎上,利用檢查點對立體影像無控制或有控制條件下的幾何定位誤差進行統(tǒng)計分析。
圖4 DSM提取流程Fig.4 DSM extraction process
通過上述方法,在GF-2異軌像對重疊區(qū)共選取了13個控制點、78個連接點,生成分辨率為4.88 m 的滑坡后DSM。將相同方法應用在ZY-302的同軌像對上,共選取15個控制點、495個連接點,生成分辨率為14.19 m的滑坡前DSM??紤]到后續(xù)體積估算疊加分析,將滑坡前后的DSM統(tǒng)一重采樣為15 m,結果如圖5所示?;潞驞SM影像中東北角和清江南岸存在一系列噪聲點,這是由于2020年12月22日的GF-2影像存在云層產(chǎn)生的測量誤差,不在滑坡區(qū)域內(nèi),對后續(xù)實驗并無影響。
圖5 滑坡前后DSMFig.5 DSM before and after the landslide
2.2.4DSM校正與精度評價
目前的研究成果通常只涉及單一時間點的DSM提取,且研究區(qū)多以低緩丘陵為主,適合開展實地精度控制工作。如胡芬等[14]利用異軌立體像對提取DSM時,使用了地面控制點半自動量測,并通過增加控制點數(shù)量提高精度。相似的,倪文儉等[15]在使用高分2號提取DSM的過程中,利用無人機搭載高精度激光雷達提取了研究區(qū)的DSM作為驗證數(shù)據(jù),該區(qū)域80%為坡度15°的緩坡。而本文需要對滑坡前后兩個時間點提取的DSM進行差值計算,由于不同的影像數(shù)據(jù)源所生成的DSM結果存在高程基準差異,有必要將其統(tǒng)一校正到參考DSM的基準上。為此,本文提出在區(qū)域網(wǎng)平差的基礎上對滑坡前后的DSM均進行二次高程校正,從而進一步降低在復雜山地地形測量上出現(xiàn)的誤差。首先,采集一定數(shù)量的校正點,以GLO30數(shù)據(jù)作為統(tǒng)一的參考DSM,使用最小二乘擬合二次多項式分別對滑坡前后的DSM進行校正;然后再另采集一定數(shù)量的驗證點對校正后的DSM做精度評價。所有校正點和驗證點均以ESRI全球衛(wèi)星影像網(wǎng)絡服務為底圖參考,通過目視解譯,在地物變化較少、無植被覆蓋且遠離滑坡崩塌區(qū)的硬質(zhì)水泥路面上選取,以減少因各種因素造成的地物變形而產(chǎn)生的高程偏差。
如圖6所示,本文共選擇了44個校正點和24個驗證點,主要均勻分布在研究區(qū)四周的道路,高程較為穩(wěn)定的區(qū)域。式(2)和式(3)分別為滑坡前后DSM與參考GLO30 數(shù)據(jù)的擬合方程:
圖6 校正點和驗證點分布Fig.6 Layout of correction points and verification points
Y=0.9837X1+6.7636
(2)
Y=0.9816X2+11.2100
(3)
式中:Y為GLO30高程值,X1和X2分別為滑坡前和滑坡后DSM原始值,上述兩式中的一次項系數(shù)均十分接近于1,說明滑坡前后的DSM與GLO30的高程差異基本為一個常數(shù)平面。
經(jīng)過高程校正后,滑坡前后的DSM高程中誤差分別為4.64 m和5.59 m,相較于校正前中誤差分別降低了5.64 m和2.45 m??紤]到GLO30 數(shù)據(jù)在非極地地區(qū)的平均誤差為1.92 m[17],滑坡前后的區(qū)域DSM誤差基本<10 m,滿足《基礎地理信息數(shù)字成果1∶5 000、1∶10 000、1∶25 000、1∶50 000、1∶100 000數(shù)字高程模型》[21]規(guī)定的1∶10 000 DEM高山地的測圖要求,可用于滑坡體積的估算。
2.2.5滑坡物質(zhì)體積估算
滑坡發(fā)生后會形成滑坡凹地和滑坡堆積區(qū)。通過對滑坡后與滑坡前的高程差分,不但可以清晰地標識出二者的位置和范圍,還能通過式(4)估算出滑坡物質(zhì)由凹地輸出或在堆積區(qū)匯集的體積。
(4)
式中:V為區(qū)域內(nèi)滑坡物質(zhì)體積總量,A,N和ΔDi分別為區(qū)域內(nèi)單個像元面積、像元總數(shù)和第i個像元的高程差。經(jīng)計算,沙子壩滑坡的凹地體積為237.01萬m3,堆積區(qū)體積為137.24萬m3,滑坡物質(zhì)體積變化分布情況如圖7所示,其中負值代表滑坡凹地,正值代表滑坡堆積區(qū)。
圖7 滑坡體積變化Fig.7 Landslide volume change
理論上,滑坡凹地輸出的物質(zhì)總量應等于滑坡堆積區(qū)匯集的物質(zhì)總量,這種等量關系反映了“物質(zhì)平衡原理”[22],而從本文計算結果可看出,堆積區(qū)體積遠小于凹地區(qū)體積,分析其原因:①該滑坡是典型的河岸滑坡,滑坡堆積區(qū)有很大一部分物質(zhì)已經(jīng)傾瀉入清江河道內(nèi),并逐漸累積形成滑坡壩,通過光學立體測量手段只能獲取滑坡壩露出水面的部分體積,因此出現(xiàn)堆積區(qū)體積偏低的情況;②GF-2影像拍攝時間距滑坡發(fā)生已間隔5個月以上,期間堆積物的清理工作已持續(xù)展開,大量堆積物已被人工清理,或被清江河水自然沖蝕搬運,致使測量獲得的堆積體積進一步減少。
綜合上述分析,要獲得較準確的滑坡物質(zhì)體積,應從地面變化較小且較容易測量的滑坡凹地體積入手。在原凹地體積的基礎上,還應考慮夾帶和刮擦效應,即滑動物質(zhì)在路徑運動過程中會夾帶和液化飽和土壤,增大實際的滑坡物質(zhì)體積[8]。參考前人的研究成果[23-25],假設這一增量為25%,最終估算得到的滑坡物質(zhì)體積為296.28萬m3。
根據(jù)圖7的體積變化分布結果,對沙子壩滑坡構造特征進行了劃分,其中:A為滑源區(qū),是滑坡的啟動區(qū),也是主要的物質(zhì)輸出區(qū),投影面積為7.29萬m2,平均高程差為-21.78 m,輸出物質(zhì)體積占總量的65.05%。B為滑動帶,是滑坡物質(zhì)移動的主要通道,后緣連接滑源區(qū),前緣位于清江北岸的高陡坡,前后緣高程差為300 m,長約1 200 m,投影面積為24.88萬m2,整體體積變化輕微,呈中部減少兩側(cè)增加的分布。其原因是:滑坡體在通過該區(qū)域的過程中發(fā)生了鏟刮和攜帶效應,裹挾位于滑動帶中部的物質(zhì),造成中部體積損失;而部分位于滑坡體邊緣的物質(zhì)由于速度和能量的衰減,逐漸停滯并堆積,這些堆積物沿滑坡東西側(cè)壁分布,造成了該處體積增加。C為滑坡壩,也是本次滑坡的主要堆積區(qū),沿清江河道展布,長約1 000 m。受衛(wèi)星攝影測量方法的局限性,圖7結果只反映了滑坡壩露出水面的體積,約為137.24萬m3。D為陡坡,位于清江北岸,連接滑動帶和滑坡壩。在已公開的對本次滑坡研究的文獻中,對滑坡體積的估算主要采用以下兩種方法:①平均厚度乘以面積法,估算得到的滑坡體積約為1 000萬m3[26];②現(xiàn)場估算法得到滑坡體積約為280萬m3[27]。方法①實施最簡單,但通常存在極大誤差,方法②則更具可信度,因此由方法②估算獲得的滑坡體積也被后續(xù)研究所采納。本文估算獲得的結果與方法②接近,偏差僅為5%,表明本文提出的滑坡體積估算方法不但具有較高的可行性和可靠的精度,且能直觀地反映滑坡體積變化的空間分布情況。
在方法設計上,本文所提出的數(shù)據(jù)處理思路可在數(shù)據(jù)缺乏、無實測控制點的情況下,充分挖掘已有數(shù)據(jù)資料,通過構建適合的異軌影像對,不依賴地面實測提取滑坡體積,可用于高山地等復雜地形的滑坡體積提取,也可進一步推廣應用于常規(guī)山地的高程測量。與此同時,我們?nèi)詰⒁饣弑葘τ诋愜壐叱虦y量精度的影響。在數(shù)據(jù)選擇上應優(yōu)先選擇大基高比的影像對。但由于衛(wèi)星的運行軌道和載荷參數(shù)并未針對異軌測量而設計等多方面限制,大基高比的像對所覆蓋的公共區(qū)域通常較為狹窄,可能無法提取出足量的高質(zhì)量地面控制點。因此,還可采取適當放寬基高比的策略,進而依靠更高分辨率的影像,并結合本文提出的二次高程校正方法彌補基高比不足帶來的精度損失。
本文以如何充分利用衛(wèi)星遙感數(shù)據(jù)資源,并準確提取滑坡物質(zhì)的體積為出發(fā)點,以2020年7月發(fā)生在湖北省恩施市馬者村沙子壩的大型滑坡為研究對象,將異軌DSM提取方法應用于傳統(tǒng)高分辨率光學衛(wèi)星GF-2影像,并利用其生成的滑坡后DSM與ZY-302衛(wèi)星生成的滑坡前同軌DSM進行差分計算獲取了滑坡體積,結果表明:
(1) 經(jīng)過外部高程數(shù)據(jù)的二次校正精度后,異軌與同軌立體測量技術所生成的DSM精度相當,均滿足1∶10 000 DEM高山地的測圖要求。滑坡前后DSM與GLO30參考DSM均存在顯著的線性相關性,均滿足滑坡體積量算的精度要求。值得注意的是,GF-2全色影像的高分辨率,彌補了基高比偏低的不足,其精度接近于ZY-302同軌測量。
(2) 提取的滑坡體積為296.28萬m3,與已公開的文獻報道相吻合,具有較高的準確性。在生成的滑坡體積變化分布圖中,清晰地呈現(xiàn)此次滑坡的構造特征,可識別出滑源區(qū)、滑動帶、陡坡和滑坡壩4部分。
本文提出基于異軌立體測量技術估算滑坡物質(zhì)體積的方法,充分利用已有衛(wèi)星遙感資料,以較低成本、較便捷的途徑,實現(xiàn)對滑坡地形的高精度遙感測量,對常規(guī)山地地形測繪具有參考價值;通過更精細地描繪出滑坡構造,對進一步認識滑坡機理具有科學意義。
致 謝
本次研究在遙感數(shù)據(jù)獲取和處理工作中得到了廣西壯族自治區(qū)自然資源遙感院的大力支持和協(xié)助,在此表示感謝。