基于異軌DSM提取技術的滑坡體積估計

庾 露,黃 艷 霞

(1.南寧師范大學 地理科學與規劃學院,廣西 南寧 530100; 2.廣西國土資源規劃設計集團有限公司,廣西 南寧 530029)

0 引 言

滑坡是常見的自然災害之一,極易誘發次生災害,對人類的生存環境造成極大的威脅,一直是災害防治工作中的重中之重[1]。據自然資源部統計,2020年全國共發生地質災害7 840起,其中滑坡4 810起,占總數的61%,滑坡災害呈常態化趨勢[2]。獲取滑坡的體積是研究滑坡災害的一個重要內容,滑坡體積對滑坡的破壞力、展布范圍及其可能引發的二次災害如堰塞湖和洪水具有指示意義[3-5],精確的滑坡測量能降低滑坡相關的危險度和風險性分析中的不確定因素并且有助于進一步的防災減災[6]。

數字表面模型(Digital Surface Model,DSM)差值計算是測量滑坡體積方法之一,通常使用遙感測繪技術獲取滑坡前后的高程分布,利用高程差和像元面積的乘積獲取滑坡體積,具有覆蓋范圍較大、時間和人力成本較低的優點。新一代遙感衛星和無人機平臺的成熟并大規模投入使用,以及LiDAR、InSAR等測量技術的發展,為地質災害監測提供了多樣的數據來源和分析手段[7]。如MA[8],彭大雷[9],CHANG[10]等使用低空無人機測量計算出滑坡體積,此方法基本不受云層遮擋,但需要攜帶設備到達現場開展工作。同時,受飛行平臺穩定性和載荷量限制,無人機遙感系統所獲取的圖像往往像幅較小、數量較多、畸變較大,圖像處理難度較高[11]。衛星遙感具有覆蓋范圍大、成本低、精度高的特點,也可作為滑坡體積估算的手段[12],但由于滑坡發生后,局部通常伴隨長時間降雨,會對衛星觀測造成諸多不利影響。傳統的光學立體測繪衛星受固有重訪周期限制、拍攝調度安排和降雨天氣云層的遮擋,往往在單顆衛星過境的短暫間隙內無法獲取有效影像。星載重復軌道InSAR測量在成像過程中,大氣水汽含量的時空變化會引起雷達信號傳播延遲,給干涉相位帶來較大的附加相位變化,導致干涉測量結果出現誤差[13]。因此,有必要根據研究區所收集的影像數據和地形環境情況,進一步拓展滑坡地形測量的手段。

基于光學衛星的異軌立體測量技術,是對傳統光學測繪衛星同軌測量方法的有效補充。該技術可將未設計用于立體測繪的衛星數據充分利用,借助單顆或多顆衛星重復過境縮短重訪周期的優勢,多期覆蓋研究區,將可能存在的有效影像篩選出來并組成異軌立體像對,生成DSM產品。已有學者使用該技術提取了地物高程,如胡芬等[14]對比分析了異軌和同軌立體影像生成的DSM,發現二者精度趨近。倪文儉等[15]使用高分2號衛星(GF-2)數據形成異軌立體像對,提取森林高度,與LiDAR結果的均方根誤差為3.6 m。

上述研究成果驗證了該技術的可行性,但目前的研究仍主要集中于單一成像時間點上的DSM提取,尚缺乏考慮滑坡事件發生前后高程劇烈變化產生的數據一致性問題,以及在地形起伏大、無實地控制點約束下的精度控制情況。因此,有必要就如何將異軌立體測量技術應用于滑坡體積提取場景做更深入的研究。本文以2020年7月11日湖北省恩施市馬者村沙子壩滑坡為例,重點分析使用高分2號衛星(GF-2)通過異軌立體成像立體測量技術在無實地控制點的場景下,提取滑坡后DSM的可行性和精度。同時,結合資源3號02星(ZY-302)影像提取的滑坡前DSM,估算滑坡體積,并對滑坡的構成進行分析。

1 研究區概況

受梅雨期連續降雨影響,2020年7月20～21 日湖北省恩施市清江上游馬者村沙子壩(東經109°18′,北緯30°22′)滑坡垮塌,發生縱向長1 200～1 500 m,橫向寬320～580 m的特大型滑坡,導致清江干流河道形成堰塞湖,水位的急速上升導致存在潰壩風險[16]。研究區地處巫山余脈和武陵山北上的余支交會部,地勢整體為西低東高,最大落差達1 380 m,地形崎嶇復雜,多為海拔 1 000 m 以上的高山。在氣候上,研究區屬亞熱帶大陸性季風氣候,雨量充沛,易發生暴雨洪澇,造成滑坡。具體的影像范圍與研究區范圍如圖1所示。

圖1 研究區域位置和影像范圍Fig.1 Location and image coverage of the study area

2 數據獲取與研究方法

2.1 數據來源

滑坡前的DSM數據由2020年3月19日的資源三號02(ZY-302)衛星影像生成,這是距離發生滑坡時間較近、成像質量較高的衛星影像數據。ZY-3和后續發射的ZY-302星均是中國民用高分辨率立體測繪衛星,通過星載三線陣光學相機,以同軌方式實現地形圖繪制。滑坡后覆蓋研究區域的可用影像較少,以ZY-3和ZY-302為例,在滑坡發生后180 d內僅有3景影像覆蓋,且滑坡區受到云層遮擋。而在此期間,高分二號(GF-2)衛星以不同角度拍攝了多景覆蓋滑坡區的影像,雖然該衛星并不是為立體測繪任務而設計的,但具有最高0.8 m的高空間分辨率和最大35°的側擺能力,在不同觀測角度和較短重復觀測周期的前提下可構成異軌立體觀測像對,從而實現DSM提取[15]。經過反復試驗,本文選取了距離滑坡發生后時間較接近的2020年11月08日和2020年12月22日的兩幅GF-2影像構成立體像對,用以提取滑坡后DSM。

本文使用ESRI公司提供的全球衛星影像網絡服務作為平面位置參考,使用GLO30 DSM數據作為高程參考。前者在研究區范圍內提供了分辨率為1 m,經過地形校正和拼接處理的衛星影像底圖,其分辨率與GF-2全色波段影像的分辨率接近,并優于ZY-302影像的分辨率,可用于生成控制點和驗證點的平面坐標;后者為歐空局哥白尼地球觀測計劃的數字表面高程產品,數據來源于TanDEM-X任務期間獲得的SAR衛星資料,GLO30的像元尺寸為30 m,絕對垂直精度<4 m,絕對水平精度<6 m[17],研究區最新數據更新至2014年可用于生成控制點和驗證點的高程值。

2.2 研究方法

2.2.1技術流程

對于滑坡區DSM提取問題,受其所處地理位置、地形地貌和自身不穩定地質環境等多方面因素制約,可能不利于開展實測工作。為此,本文提出完全基于衛星遙感數據源進行異軌DSM提取和精度校正的新思路,并在后文的論述中討論該方法的可行性。本次研究的技術流程如圖2所示,使用1景ZY-3 02衛星影像提取滑坡前DSM,使用2景異軌的GF-2衛星影像提取滑坡后DSM,使用ESRI衛星全球衛星影像作為平面位置參考,GLO30作為高程參考對提取的DSM進行位置和高程校正和精度評價,最后計算滑坡體積。

圖2 技術路線Fig.2 Technical route

2.2.2異軌立體測量基本原理

異軌立體測量是基于光學衛星所具備的側擺功能,由單臺相機在不同軌道觀測獲取立體影像的方法。如圖3所示,在兩個不同軌道上分別對目標成像,產生類似于雙目立體視覺,通過獲取兩個軌道的基線、焦距以及視差等信息可計算出地物的三維點云分布。異軌立體測量的精度包括平面精度和高程精度。其中平面精度和同軌立體測量方式接近[18-20],可通過衛星數據供應商提供的RPC(Rational Polynomial Cofficient)參數,運用成像幾何模型,并配合一定數量的地面控制點建立圖像坐標和地理坐標的幾何關系。而高程精度則主要受飛行平臺的軌道精度、影像空間分辨率、成像紋理清晰度、基高比等影響。對于異軌立體的測量,通常前3個因素基本相同或接近,而基高比可根據所收集的影像進行較靈活的調整。理論上基高比越大,交會條件越好,高程測量的精度越高。根據公式(1)基高比B/H由兩軌影像的側視角α和β決定:

圖3 基高比示意Fig.3 Baseline-to-height ratio schematic

(1)

式中:B為軌道間基線長度,H為攝站高度,α和β分別為側視角。經計算,GF-2異軌像對的基高比約為0.50,ZY-302衛星立體像對的基高比約為0.83。

2.2.3利用GF-2提取DSM

由于本次滑坡體積提取實驗未到達現場采集地面控制點,因此用于高精度幾何定位的控制點和用于DSM結果精度評估的驗證點均通過已公開的免費遙感數據產品采集獲得。本文使用ArcGIS Desktop載入全球衛星影像網絡服務和GLO30 DSM數據,采集地形起伏較小、地表紋理清晰穩定、易于辨認的人造地物(如道路交叉口),作為控制點和驗證點。使用Trimble Inpho軟件,導入控制點,對GF-2衛星異軌立體影像采用區域網平差的方式進行高精度的幾何定位,其流程如圖4所示。 首先,通過特征匹配方法自動提取大量均勻分布的連接點,實現立體影像密集匹配,提高匹配的效率和可靠性,降低誤匹配率;通過自由網平差定向,對觀測像點的殘差進行統計分析,進而采用選權迭代法實現誤匹配連接點的自動識別與剔除。然后,進行控制點半自動量測,通過區域網平差定向和定向殘差分析,迭代計算得到區域網平差定向結果。最后,在自由網平差或區域網平差的基礎上,利用檢查點對立體影像無控制或有控制條件下的幾何定位誤差進行統計分析。

圖4 DSM提取流程Fig.4 DSM extraction process

通過上述方法,在GF-2異軌像對重疊區共選取了13個控制點、78個連接點,生成分辨率為4.88 m 的滑坡后DSM。將相同方法應用在ZY-302的同軌像對上,共選取15個控制點、495個連接點,生成分辨率為14.19 m的滑坡前DSM。考慮到后續體積估算疊加分析,將滑坡前后的DSM統一重采樣為15 m,結果如圖5所示。滑坡后DSM影像中東北角和清江南岸存在一系列噪聲點,這是由于2020年12月22日的GF-2影像存在云層產生的測量誤差,不在滑坡區域內,對后續實驗并無影響。

圖5 滑坡前后DSMFig.5 DSM before and after the landslide

2.2.4DSM校正與精度評價

目前的研究成果通常只涉及單一時間點的DSM提取,且研究區多以低緩丘陵為主,適合開展實地精度控制工作。如胡芬等[14]利用異軌立體像對提取DSM時,使用了地面控制點半自動量測,并通過增加控制點數量提高精度。相似的,倪文儉等[15]在使用高分2號提取DSM的過程中,利用無人機搭載高精度激光雷達提取了研究區的DSM作為驗證數據,該區域80%為坡度15°的緩坡。而本文需要對滑坡前后兩個時間點提取的DSM進行差值計算,由于不同的影像數據源所生成的DSM結果存在高程基準差異,有必要將其統一校正到參考DSM的基準上。為此,本文提出在區域網平差的基礎上對滑坡前后的DSM均進行二次高程校正,從而進一步降低在復雜山地地形測量上出現的誤差。首先,采集一定數量的校正點,以GLO30數據作為統一的參考DSM,使用最小二乘擬合二次多項式分別對滑坡前后的DSM進行校正;然后再另采集一定數量的驗證點對校正后的DSM做精度評價。所有校正點和驗證點均以ESRI全球衛星影像網絡服務為底圖參考,通過目視解譯,在地物變化較少、無植被覆蓋且遠離滑坡崩塌區的硬質水泥路面上選取,以減少因各種因素造成的地物變形而產生的高程偏差。

如圖6所示,本文共選擇了44個校正點和24個驗證點,主要均勻分布在研究區四周的道路,高程較為穩定的區域。式(2)和式(3)分別為滑坡前后DSM與參考GLO30 數據的擬合方程:

圖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原始值,上述兩式中的一次項系數均十分接近于1,說明滑坡前后的DSM與GLO30的高程差異基本為一個常數平面。

經過高程校正后,滑坡前后的DSM高程中誤差分別為4.64 m和5.59 m,相較于校正前中誤差分別降低了5.64 m和2.45 m。考慮到GLO30 數據在非極地地區的平均誤差為1.92 m[17],滑坡前后的區域DSM誤差基本<10 m,滿足《基礎地理信息數字成果1∶5 000、1∶10 000、1∶25 000、1∶50 000、1∶100 000數字高程模型》[21]規定的1∶10 000 DEM高山地的測圖要求,可用于滑坡體積的估算。

2.2.5滑坡物質體積估算

滑坡發生后會形成滑坡凹地和滑坡堆積區。通過對滑坡后與滑坡前的高程差分,不但可以清晰地標識出二者的位置和范圍,還能通過式(4)估算出滑坡物質由凹地輸出或在堆積區匯集的體積。

(4)

式中:V為區域內滑坡物質體積總量,A,N和ΔDi分別為區域內單個像元面積、像元總數和第i個像元的高程差。經計算,沙子壩滑坡的凹地體積為237.01萬m3,堆積區體積為137.24萬m3,滑坡物質體積變化分布情況如圖7所示,其中負值代表滑坡凹地,正值代表滑坡堆積區。

圖7 滑坡體積變化Fig.7 Landslide volume change

理論上,滑坡凹地輸出的物質總量應等于滑坡堆積區匯集的物質總量,這種等量關系反映了“物質平衡原理”[22],而從本文計算結果可看出,堆積區體積遠小于凹地區體積,分析其原因:①該滑坡是典型的河岸滑坡,滑坡堆積區有很大一部分物質已經傾瀉入清江河道內,并逐漸累積形成滑坡壩,通過光學立體測量手段只能獲取滑坡壩露出水面的部分體積,因此出現堆積區體積偏低的情況;②GF-2影像拍攝時間距滑坡發生已間隔5個月以上,期間堆積物的清理工作已持續展開,大量堆積物已被人工清理,或被清江河水自然沖蝕搬運,致使測量獲得的堆積體積進一步減少。

綜合上述分析,要獲得較準確的滑坡物質體積,應從地面變化較小且較容易測量的滑坡凹地體積入手。在原凹地體積的基礎上,還應考慮夾帶和刮擦效應,即滑動物質在路徑運動過程中會夾帶和液化飽和土壤,增大實際的滑坡物質體積[8]。參考前人的研究成果[23-25],假設這一增量為25%,最終估算得到的滑坡物質體積為296.28萬m3。

3 結果分析與討論

根據圖7的體積變化分布結果,對沙子壩滑坡構造特征進行了劃分,其中:A為滑源區,是滑坡的啟動區,也是主要的物質輸出區,投影面積為7.29萬m2,平均高程差為-21.78 m,輸出物質體積占總量的65.05%。B為滑動帶,是滑坡物質移動的主要通道,后緣連接滑源區,前緣位于清江北岸的高陡坡,前后緣高程差為300 m,長約1 200 m,投影面積為24.88萬m2,整體體積變化輕微,呈中部減少兩側增加的分布。其原因是:滑坡體在通過該區域的過程中發生了鏟刮和攜帶效應,裹挾位于滑動帶中部的物質,造成中部體積損失;而部分位于滑坡體邊緣的物質由于速度和能量的衰減,逐漸停滯并堆積,這些堆積物沿滑坡東西側壁分布,造成了該處體積增加。C為滑坡壩,也是本次滑坡的主要堆積區,沿清江河道展布,長約1 000 m。受衛星攝影測量方法的局限性,圖7結果只反映了滑坡壩露出水面的體積,約為137.24萬m3。D為陡坡,位于清江北岸,連接滑動帶和滑坡壩。在已公開的對本次滑坡研究的文獻中,對滑坡體積的估算主要采用以下兩種方法:①平均厚度乘以面積法,估算得到的滑坡體積約為1 000萬m3[26];②現場估算法得到滑坡體積約為280萬m3[27]。方法①實施最簡單,但通常存在極大誤差,方法②則更具可信度,因此由方法②估算獲得的滑坡體積也被后續研究所采納。本文估算獲得的結果與方法②接近,偏差僅為5%,表明本文提出的滑坡體積估算方法不但具有較高的可行性和可靠的精度,且能直觀地反映滑坡體積變化的空間分布情況。

在方法設計上,本文所提出的數據處理思路可在數據缺乏、無實測控制點的情況下,充分挖掘已有數據資料,通過構建適合的異軌影像對,不依賴地面實測提取滑坡體積,可用于高山地等復雜地形的滑坡體積提取,也可進一步推廣應用于常規山地的高程測量。與此同時,我們仍應注意基高比對于異軌高程測量精度的影響。在數據選擇上應優先選擇大基高比的影像對。但由于衛星的運行軌道和載荷參數并未針對異軌測量而設計等多方面限制,大基高比的像對所覆蓋的公共區域通常較為狹窄,可能無法提取出足量的高質量地面控制點。因此,還可采取適當放寬基高比的策略,進而依靠更高分辨率的影像,并結合本文提出的二次高程校正方法彌補基高比不足帶來的精度損失。

4 結 論

本文以如何充分利用衛星遙感數據資源,并準確提取滑坡物質的體積為出發點,以2020年7月發生在湖北省恩施市馬者村沙子壩的大型滑坡為研究對象,將異軌DSM提取方法應用于傳統高分辨率光學衛星GF-2影像,并利用其生成的滑坡后DSM與ZY-302衛星生成的滑坡前同軌DSM進行差分計算獲取了滑坡體積,結果表明:

(1) 經過外部高程數據的二次校正精度后,異軌與同軌立體測量技術所生成的DSM精度相當,均滿足1∶10 000 DEM高山地的測圖要求。滑坡前后DSM與GLO30參考DSM均存在顯著的線性相關性,均滿足滑坡體積量算的精度要求。值得注意的是,GF-2全色影像的高分辨率,彌補了基高比偏低的不足,其精度接近于ZY-302同軌測量。

(2) 提取的滑坡體積為296.28萬m3,與已公開的文獻報道相吻合,具有較高的準確性。在生成的滑坡體積變化分布圖中,清晰地呈現此次滑坡的構造特征,可識別出滑源區、滑動帶、陡坡和滑坡壩4部分。

本文提出基于異軌立體測量技術估算滑坡物質體積的方法,充分利用已有衛星遙感資料,以較低成本、較便捷的途徑,實現對滑坡地形的高精度遙感測量,對常規山地地形測繪具有參考價值;通過更精細地描繪出滑坡構造,對進一步認識滑坡機理具有科學意義。

致 謝

本次研究在遙感數據獲取和處理工作中得到了廣西壯族自治區自然資源遙感院的大力支持和協助,在此表示感謝。