庾 露,黄 艳 霞

(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部分。

本文提出基于异轨立体测量技术估算滑坡物质体积的方法,充分利用已有卫星遥感资料,以较低成本、较便捷的途径,实现对滑坡地形的高精度遥感测量,对常规山地地形测绘具有参考价值;通过更精细地描绘出滑坡构造,对进一步认识滑坡机理具有科学意义。

致 谢

本次研究在遥感数据获取和处理工作中得到了广西壮族自治区自然资源遥感院的大力支持和协助,在此表示感谢。