引言

我国大陆发生的M＞6.5地震一般会在地表产生数千米至数十千米的地表破裂带（铁瑞等，2016），并会对沿线的房屋建筑、建设工程等产生直接破坏，这种地震灾害主要是由于活动断层突发的快速错动而形成的强烈地面运动和地表破裂变形引起，同时还会在地貌上形成线性陡坎、断错水系、扭动冲沟和地表裂隙等地震遗迹，这些永久保留下来的、类型丰富的断错地貌记录着断层活动的重要信息，如地震期次、活动类型及位移量等，结合年代学和震级经验公式，可进一步讨论地震震级、计算大地震复发间隔等（邓起东，闻学泽，2008），这对于掌握震后灾害情况、预测未来地震复发模式、评估地震危险性和提高城市防震减灾能力等具有重要意义。因此震后快速获取地震地表破裂带的数字高程模型（digital elevation model，缩写为DEM）和地貌数据（正射影像），提取发震断层的定量参数，是分析地震灾情、评估地震断层危害性等工作的重要依据。

影像数据一般通过传统的卫星遥感平台和航空平台（李云等，2011）来获取，但这样获取的影像数据只能从中提取水平方向的形变参数，不能获得垂直方向的形变参数（Zhou et al，2015）。随着机载激光雷达（light detection and ranging，缩写为LiDAR）技术的出现，获取地震形变参数有了新的技术手段。Oskin等（2012）利用LiDAR技术首次获得了El Mayor-Cucapah 地震地表破裂在垂直方向的形变参数。这种遥感技术不仅能够主动实现大范围的影像数据采集，而且可以剔除植被覆盖的影响，因此适用于断错地貌的精细化和定量化研究中（刘静等，2013；Lin et al，2013），但由于该方法成本高、数据处理流程复杂，广泛应用受限。近年来在计算机视觉理论及特征匹配算法的发展中，出现了一种新型数字摄影测量技术SfM （Snavely et al，2008；Harwin，Lucieer，2012；James，Robson，2012；Fonstad et al，2013）。该技术通过自动解算相机的位置和姿态可以快速地获取高精度的三维地形地貌数据，具有作业成本低、获取效率高和操作方式便捷等特点，已被广泛用于获取DEM数据（Klinger et al，2011），同时还被应用于多个地震研究领域，如滑坡监测（Lucieer et al，2014）、古地震分析（Ren et al，2015；高伟等，2017）、微构造地貌解译（Johnson et al，2014；孙稳等，2019）和活动构造（Bemis et al，2014；Angster et al，2016；毕海芸等，2017；艾明等，2018）。前人利用遥感影像获取地震应急灾情信息的大量研究主要集中在灾区建筑物、道路和滑坡方面（李金香等，2019；荆帅军等，2019），鲜有基于无人机摄影测量技术提取地震地表破裂带的定量参数研究。

为此，本文针对大震后如何快速获取地表破裂带的影像数据以及如何生成高精度、高分辨率的DEM数据这一实际问题，拟利用大疆无人机采集的影像数据，以集成SfM算法的PhotoScan软件为平台，通过对比不同质量密集点云生成DEM数据的优缺点来确定软件处理的最优参数配置，然后以中卫天景山活动断裂带上1709年中卫南M7½地震的主体地表破裂带为例，提取其上地震断层的定量参数，并结合前人研究结果讨论所获取参数的可靠性及该方法的可行性，以期揭示无人机摄影测量技术在地震应急方面的应用前景。

1.   影像数据的获取与处理

1.1   影像数据的获取

本次使用的无人机为精灵Phantom 4 Pro V2.0，该无人机搭载了大疆制造的FC6310S型号的相机，配备了1英寸2 000万像素的影像传感器和8.8 mm焦距的广角镜头，确保获取到分辨率高、色彩还原度好及物体成像变形小的影像数据；其质量（含电池和桨）为1 375 g，质量轻、体积小使得在地震应急中便于搬运、装卸及携带；最大飞行海拔高度为6 000 m，最大可承受风速为10 m/s，工作环境温度为0℃—40℃，再加上FlightAutonomy安全系统，该无人机能够满足高原、山地上的正常作业；集成了全球定位系统(Global Positioning System，缩写为GPS）和全球卫星导航系统（Global Navigation Satellite System，缩写为GLONASS）卫星定位双模块，GPS定位正常时水平误差为±1.5 m，垂直误差为±0.5 m，在震后的恶劣环境下，无需满足对控制点精度（魏占玉等，2015）的要求，便可以快速获取地表破裂带高精度、高分辨率的影像数据。

无人机摄影图像的质量与环境条件、无人机的航迹、高度、重叠度等有关（陆博迪等，2011）。此次采集的地表破裂带地处山地，选择晴朗弱风天气的中午时分进行航拍，以保证地貌上的光线最佳，无人机飞行的稳定性好。由于山地地形起伏较大且地面如电线杆等障碍物较多，利用自动飞行软件设置飞行高度为120 m、航向和旁向的重叠度均为80%，这样既能保障飞行的安全性，又能确保影像的分辨率。受限于电池的电量，无人机一次飞行的覆盖区域不超过500 m×500 m，耗时一般在15 min左右，研究区域可分三四个架次完成，相邻架次采集的覆盖区域必须有一定的重叠，这样便于后期的影像拼接。

1.2   影像数据的处理

本文选取中卫市西天景山断裂带某段的41张影像照片作为试验数据，其覆盖面积为1.28×10⁵ m²，无人机的平均飞行高度为146 m。利用集成SfM算法的PhotoScan软件，在8 G内存、英特尔酷睿i7双核处理器的图形工作站处理数据，经过对齐照片、生成密集点云、生成网格、建立DEM和获取正射影像等五个步骤完成拼接，其中前三个步骤的参数设置直接决定了所生成的DEM和正射影像的分辨率，详见表1。各步骤分述如下：

表  1  对齐照片、生成密集点云和生成网格的参数设置

Table  1.  Parameter setting for aligning photo，generating dense point clouds，and generating grids

对齐照片			生成密集点云			生成网格
精度	成对预选		质量	深度过滤		表面类型	源数据	面数
高	参考		超高，高，中	进取		高度场	密集点云	中

下载: 导出CSV 

| 显示表格

1） 对齐照片。当精度设置为“高”时，直接处理的是原始影像。由于影像是具有测量坐标的照片，故成对预选为“参考”模式。

2） 生成密集点云。这一步是获取点高程的关键，点云密度越大，生成的DEM分辨率越高，软件处理所需时间越长。本研究通过设置“超高”、“高”和“中”三种点云密度，在保持其它参数不变的情况下对比生成的DEM在水平位置和高程上的差异，以挑选出最优的参数配置。超高质量表示处理原始影像，高质量处理的影像尺寸为超高质量的四分之一，中等质量处理的影像尺寸则为超高质量的十六分之一；深度过滤是为了挑选出影像对焦不好的异常点，以滤掉地表植被覆盖产生的异常点，这里选择适用于航空摄影数据处理的“进取”过滤模式。

3） 生成网格。处理航空摄影测量数据时，选择“高度场”表面类型模式，其对内存的要求更低；数据源选择上一步的密集点云，这样处理结果会更好，但也意味着处理时间长。最终生成网格的最大面数是基于上一步生成的密集点云中点的数量计算出来，“中”的比率为十五分之一最为合适，因为面值太小生成的网格过于粗糙，而面值太大即超过1千万会影响模型的可视化。

基于上述各步骤设定的参数，获取了相应的DEM和正射影像，将分辨率和各类参数进行对比，结果列于表2。可以看出：密集点云的质量仅对DEM的点云密度和分辨率产生影响，而对正射影像的分辨率并无影响；如果不考虑处理时长，密集点云的质量显然越高越好，但为了震后快速获取DEM的同时又能保证厘米级的分辨率，处理时长不应太长。从超高质量到高质量密集点云生成的DEM，虽然其分辨率降低了50%，但获取时间缩短了80%，故从处理时间和分辨率方面综合考虑高质量密集点云为最佳选择。

表  2  不同质量密集点云生成的DEM和正射影像的分辨率及其它相关参数

Table  2.  Resolution of DEM and orthophoto map generated by dense point clouds of different masses and the other parameters

密集点云 的质量	处理时间/h	DEM		正射影像分辨率 /（cm·pix−1）
		点云密度/（点·m−2）	分辨率/（cm·pix−1）
超高	15	876	3.38	3.38
高	3	219	6.76	3.38
中	0.6	55	13.50	3.38

下载: 导出CSV 

| 显示表格

1.3   不同质量密集点云生成的DEM水平位置和高程的差异对比

在缺少地面控制点的情况下，首先通过图像特征匹配与跟踪解算出相机的空间位置，然后进一步计算出三维密集点的相对坐标，该坐标为图像的空间坐标系并非现实世界的空间坐标系，故解算出的相机位置决定了密集点云的坐标系。1.2节中步骤一“对齐照片”用于解算相机位置，因此三种质量的密集点云均是基于相同的相机位置而生成，其相机位置的水平（xy）误差和高程（z）误差均相同（表3），但密集点的差异使得生成DEM时进行插值计算所得的坐标位置和高程略有不同，由表4可见三种质量密集点云生成的DEM中点的投影坐标和高程范围的差异均小于0.100 m。

表  3  相机位置误差估计值

Table  3.  The error estimates of camera locations

x误差/m	y误差/m	xy误差/m	z误差/m
0.540	0.338	0.637	3.127

下载: 导出CSV 

| 显示表格

表  4  DEM中点的投影坐标和高程范围

Table  4.  The projected coordinates and elevation ranges of midpoint in DEM

不同质量密集点云	x/m	y/m	最低高程/m	最高高程/m	高程差/m
超高	495639.212	4142958.128	1281.045	1338.719	57.674
超高与高的差异	0.000	0.007	0.057	0.079	0.022
高	495639.212	4142958.135	1281.102	1338.798	57.696
高与中等的差异	0.089	0.041	0.008	0.046	0.038
中等	495639.123	4142958.176	1281.094	1338.752	57.658

下载: 导出CSV 

| 显示表格

综合以上数据分析，在利用PhotoScan软件处理时，基于不同质量密集点云生成的DEM具有差异性。虽然超高质量密集点云可以处理原始影像，能够将软件引起的误差降至最低，但在实际工作中耗时太长无法满足地震应急的需求；中等质量密集点云虽然处理影像的时间短，但获取的DEM分辨率无法达到厘米级；高质量较超高质量密集点云的处理时间大大缩短，且获取的DEM分辨率能够达到厘米级，水平位置和高程与另外两者的差值均小于0.100 m，是快速获取高精度、高分辨率DEM和正射影像的最优参数配置。

2.   1709年中卫南M7½地震地表破裂带定量参数的提取

以1709年宁夏中卫市南发生的M7½ 地震为例来提取地表破裂带的定量参数。中卫地震产生的整个地表破裂带长约60 km，其主体地表破裂带为位于天景山活动断裂带的中段粉石沟—双井子一带（图1），该破裂带长约30 km （周俊喜，刘百箎，1987；张维歧等，1988）。前人对天景山活动断裂带地震事件的研究结果表明该断裂带上唯一的历史大地震即为1709年中卫南M7½地震（汪一鹏等，1990；闵伟等，2001），因此近千年以来形成的小冲沟位错仅为1709年地震造成，这为本研究提供了合适的场点。

图  1  天景山活动断裂带及1709年中卫南M7½地震的宏观震中

Figure  1.  Tianjingshan active fault zone and macroseismic epicenter of 1709 M7½ southern Zhongwei earthquake

下载: 全尺寸图片 幻灯片

在1709年中卫南地震主体地表破裂带上的形变遗迹有地震陡坎及水系错动，特别是青蛇崖西至红谷梁一带（图1）发育最好。因此在红谷梁西的洪积扇面上，进行无人机摄影测量获取影像数据，之后利用PhotoScan软件，基于高质量密集点云生成DEM和正射影像，按照1.2节的方法进行参数设置，最终生成分辨率为6.33 cm/pix的DEM，其点云密度为249 点/ m²，有效重叠约7.6次，影像覆盖度在9张以上的面积占区域覆盖总面积的75% （图2a），正射影像分辨率为3.16 cm/pix （图2b），通过厘米级、高分辨率的DEM和正射影像可以清晰地显示冲沟和陡坎地貌，在此基础上能够准确地解译地震地表破裂并提取相关参数。

图  2  影像数据处理生成的成果图

（a） 相机位置和摄影照片覆盖度；（b） 正射影像图

Figure  2.  Results of processing image data

（a） Camera location and image overlapping；（b） Digital orthophoto map

下载: 全尺寸图片 幻灯片

2.1   垂直位错量

从获取的高精度DEM和正射影像上来看，在山前的洪积扇上，有一条明显的东西向延伸线，其展布不受地形影响并伴随一系列冲沟的左旋扭动，汪一鹏等（1990）证实该线性地貌是跨过1709年地震地表破裂带的一条断层陡坎。在该陡坎附近保存完好、未被植被破坏的同一时期沉积扇面上，提取垂直于断层方向的地形剖面数据，以确定该次地震的垂直断错量，剖面位置见图3a。结合地形剖面的高度和坡度，不仅能够快速地分离断层的上、下盘，而且地形剖面的坡度能够反映断层陡坎的坡降是否连续。由图3b和图3c可见l₁剖面的坡度曲线仅有一个拐点，而l₂剖面的坡度曲线有两个拐点，这说明l₁剖面仅仅经历和记录了1709年地震的信息，而l₂剖面除记录最新一次地震事件外，可能还记录了之前的地震或滑坡等事件（孙稳等，2019），但是相距仅15 m的两剖面所经历的事件不应有较大差别，那么多出来的拐点应该是剖面中的异常点，实际上为后期人为和植被对地形的破坏。经对比可知l₁剖面上原始地形保留得较好，所以将其上所测量的垂直断距作为地震的垂直位错量更合理。于是，将l₁剖面图上已划分出的上、下盘地形线进行最小二乘法拟合，分别用拟合公式（图3b，c）计算两点的高程之差h₁和h₂，然后利用断层陡坎高度＝（h₁＋h₂）/2、断距测量误差＝（h₁－h₂）/2确定垂直位错量为（0.704±0.293）m。

图  3  由地形剖面确定的垂直断错量

图（b）和（c）中，黑色虚线为上、下盘拟合线，蓝色虚线是根据坡度绘制的上、下盘分界线（a） 地形剖面位置；（b） 由l₁剖面确定的垂直断错量；（c） 由 l₂剖面确定的垂直断错量

Figure  3.  The vertical displacement of the fault extracted from the terrain profiles

In Figs. （b） and （c），the black dashed lines are the fitting lines for the hanging wall and footwall，and the blue vertical lines are the boundary lines between hanging wall and footwall determined by slope. （a） The position of the terrain profiles l₁ and l₂；（b） The vertical displacement of the fault from the l₁；（c） The vertical displacement of the fault from the l₂

下载: 全尺寸图片 幻灯片

2.2   水平位移量

水系对于断层新活动反应十分灵敏，而且沿断层分布也最为广泛，因此在实际测量时以冲沟为标志。结合在高精度DEM基础上生成的坡度图和1 m等高线图，确定断层的位置和走向，然后以沟心线为参照测量该断层的水平走滑位移。大地震发生后，通过确定宏观震中的最大水平位移量，可以结合其它相关参数进一步求取大地震的平均重复间隔。

在高精度DEM的基础上，生成坡度图和1 m等高线图（图4a，4b），确定断裂带的位置和走向之后，测量两条相距30 m的冲沟的水平位移，结果为10.4 m和5.1 m。虽然两条冲沟的长度基本相当，但东侧冲沟上游逐渐变小、消失、未切到基岩，西侧冲沟的深度和宽度较大、沟头已下切入石炭系基岩，所以西侧冲沟的形成年代显然要早，且位移量为东侧冲沟的两倍，故推测其可能经历了包括1709年地震在内的两次地震事件。从横向对比来看，海原断裂带1920年M8.6地震震中水平位移普遍处于6—8 m之间，而1709年震中位移多介于4.0—4.5 m之间，其最大水平位移不会超过1920年地震，所以10.4 m的位移可能为两次错动的叠加，而东侧5.1 m的位移才是1709年地震的水平位移量。

图  4  基于坡度图和等高线图确定的水平位移量及野外实测水平位移值

（a） 坡度图；（b） 坡度图与1 m间隔等高线的叠置图；（c） 水平位错量；（d） 野外实测水平位移

Figure  4.  The horizontal displacement of the fault extracted from the superposition of slope map and 1 m interval contour map as well as that measured in the field

（a） The slope map；（b） Overlapping of slope map and 1 m interval contour map；（c） The horizontal displacement of the fault；（d） The horizontal displacement measured in the field

下载: 全尺寸图片 幻灯片

2.3   地表破裂带定量参数的探讨与无人机摄影测量技术的应用前景

张维歧等（2015）认为1709年大地震的宏观震中在主体地表破裂段的阴洞沟至红谷梁之间（图1），断错地貌上的垂直断距在0.6—0.8 m之间，水平位移在3.4—5.6 m之间，最大水平位移为5.6 m，位于红谷梁西，平均水平位移为4.05 m。本文基于无人机摄影测量技术提取红谷梁西的垂直位错为（0.704±0.293） m，水平位移为5.1 m，均与前人研究结果相吻合，同时与东侧冲沟野外实测水平位移值亦相同（图4d），因此该数据真实可靠，可以代表1709年地震的同震位移。

该方法利用无人机摄影测量技术快速获取高精度、高分辨率的DEM数据和正射影像，能够精细地判识亚米级断错地貌，准确地定位规律性冲沟扭动，清晰地再现地震断层陡坎；在大地震发生后的恶劣环境下，当人员无法到达宏观震中时，利用无人机第一时间获取发震断层的影像数据，不仅可为震后获取灾害情况提供有价值的图件和数据，更是为今后的科研保留一份珍贵的资料，且兼具安全性、及时性、重要性和科学性。利用无人机提取的定量参数，结合年代学确定的断裂在对应地质时期以来的平均滑动速率，通过提取大量的位移值确定地震的最大位移和平均位移，再根据地质方法（华莱士方法）确定的计算公式可以估算大地震的平均重复间隔；结合进一步的槽探结果可以讨论地震事件的震级，这些均展现出该方法获取的地震地表破裂带定量参数在未来地震危险性评估及地震地质等方面的研究中具有良好的应用前景。

3.   讨论与结论

本文简述了使用无人机摄影测量技术获取大地震后地表破裂带影像数据的流程，详细介绍了在集成SfM算法的PhotoScan软件平台上影像数据处理时各步骤的参数设置，深入讨论了不同质量密集点云生成的DEM数据所需处理时间、分辨率、水平位置和高程的差异，最终确定了快速获取高精度、高分辨率DEM和正射影像的最优参数配置。之后以1709年中卫南M7½大地震的主体地表破裂带为例，提取其上地震断层的垂直位错量和水平位移量，为大地震发生后基于无人机摄影测量技术快速提取发震断层的定量参数提供了可行、高效和科学的技术方法，并探讨了获取的定量参数在进一步研究和应用中的广阔前景。本文研究结果如下：

1） 无人机摄影测量技术获取影像数据具有成本低、安全性高、操作简易且获取速度快的特点。利用PhotoScan软件进行数据处理时，高质量密集点云能够生成厘米级分辨率的DEM，处理时间不需太长，水平位置和高程与另外两种质量密集点云生成的DEM差异均小于0.100 m，能够满足地震应急中快速获取地震地表破裂带影像数据及处理的需求。

2） 利用无人机摄影测量技术获取1709年中卫南M7½地震地表破裂带的影像数据，基于SfM算法生成DEM和正射影像，所提取地震断层的垂直位错量为（0.704±0.293） m、水平位移量为5.1 m，与张维歧等（2015）的研究结果相吻合，因此该数值能够代表1709年地震的同震位移，证实了该方法的可靠性和科学性。

使用无人机摄影测量技术在震后快速获取地震地表破裂带的影像数据、提取地震断层的定量参数，是将该技术用于活动构造定量化、精细化研究的延伸。利用生成的DEM能够提取地形剖面和坡度信息，在地震应急中获取传统遥感影像难以直接获得的破裂带垂直位移信息，显示出该技术巨大的应用潜力，但目前还存在一定的局限性，一些重要的因素如天气、飞行高度、测量面积、续航时间和地面控制点等，在大地震后恶劣的环境条件下，还会对影像的获取和精度的提高产生较大影响，尚需要进一步研究。未来工作可从飞行高度和地面控制点方面入手，探讨以适宜的飞行高度来获取高分辨率的影像数据以及选取适量的地面控制点来提取高精度的DEM，使无人机摄影测量技术在地震应急中发挥出更大的作用。

两位审稿专家提出了宝贵的修改意见和建议，作者在此表示衷心的感谢！