本研究以位于新疆维吾尔自治区哈密市天山南坡的堆积型戈壁为研究对象，通过自一条洪积扇的扇心到靠近扇缘分别选取3个典型的样区(图 1)，开展戈壁近地面无人机飞行实验。样区1(扇心)、样区2(扇中)、样区3(靠近扇缘)的平均海拔分别为1 735 m，1 436 m，1 209 m。

图  1  研究区地理位置及样区布设

Figure 1.  Location of study area and sample zones

2016年8月21日到8月29日，利用大疆无人机(DJI Inspire 1 RAW)在3个样区分别进行了戈壁表面图像采集，并利用高精度GNSS设备(南方极智X6)进行地面控制点测量。

每个样区宽度固定为30 m，长度为1 km左右。无人机飞行的相对高度设置为20 m，共4条航线，两条相邻航线的间距为10 m。为了尽可能保证无人机距离地面的相对高度恒定，即保证采集图像的地面分辨率一致，在每个飞行小区内，沿样线方向均匀地内插了4个飞行高度控制点，更加精确地约束了无人机的飞行轨迹。

无人机的飞行速度为4 m·s^-1，拍照间隔为2 s。无人机获取的单景航片大小为4 000×2 250像素，对应的地面大小为29.06 m×16.35 m。根据飞行设置及航片大小，计算得到无人机航拍的航向重叠率为72.47%，旁向重叠率为38.84%。

为了在后期无人机图像处理中提高正射图像和地形数据的精度，使用中心喷绘有十字的方形发泡板进行地面控制点布设。每个飞行任务地面布设5个控制点(上端2个，下端2个，中间1个)，并运用连接新疆CORS网络系统后精度可达厘米级的南方极智X6进行地面控制点的经度、纬度和高程信息采集。

本研究无人机正射图像的生成利用基于SfM技术的AgiSoft LLC公司的软件PhotoScan 1.2.6完成。图像处理流程如图 2所示，具体步骤如下：采用通用成对预选进行高精度的照片对齐；利用添加的地面控制点坐标对相机对齐进行优化；采用进取模式进行高精度密集点云构建并生成数字表面模型(Digital Surface Model, DSM)；最后以DSM模型作为表面模型进行正射校正，以生成正射图像。

图  2  无人机正射图像拼接流程图

Figure 2.  Flow chart of UAV image processing

按照“砾石形貌特征高效提取技术”^[15]的基本流程，进行了基于无人机图像的砾石特征参数提取。砾石特征提取主要包含训练数据采集、构建决策树分类模型、砾石分类、砾石分割和特征参数计算5个关键步骤。

利用ArcGIS 10.2(ESRI, 美国)的矢量工具分别进行了前景(砾石)和背景(非砾石)的训练数据勾绘，采集不同颜色、光照条件下的砾石像元及多种颜色和光照下的植被、土壤等非砾石像元。

通过色彩空间变换，将原来的数据由RGB空间扩展到HSV和CIE L^*a^*b^*空间，即每个像元对应了红(R)、绿(G)、蓝(B)、色调(H)、饱和度(L)、明度(V)、亮度(L)以及有关色彩的a和b共9种特征量。然后利用CART算法，进行决策树分类模型构建。

为了降低噪声的影响，减少错分像元，对分类后图像进行前景和背景的噪声去除。将前景中像元数量小于10的对象去除，背景中像元数量小于10的空洞填充。通过砾石像元与图像中非空像元的比值计算得到砾石覆盖度。

为了提取出单个砾石，对分类后的二值图像利用分水岭算法进行图像分割。分割完成后，单个对象即为单个砾石。

对每个对象进行特征参数计算，以进行筛选及分析。计算的参数包括面积(mm²)，紧实度，Feret最大和最小粒径(mm)^[16]4项参数。根据正射图像的分辨率(7 mm·像元^-1左右)，筛选出紧实度不小于0.8且不大于0.95、Feret最小粒径大于32 mm、砾石面积大于10个像元的砾石。最后，输出所有砾石的Feret最大粒径即为砾石粒径。

为了分析砾石的空间分布，利用PhotoScan在生成密集点云之后，对点云进行分类得到地面点点云(单元大小1.5 m，最大距离0.1 m，最大角度15°)。利用地面点点云生成高精度的DTM，并在此基础上利用ArcGIS 10.2生成坡度和坡向数据。由于DTM数据地面分辨率过高，而坡度和坡向数据都是根据目标像元为中心的9个像元值来计算的，对原始DTM数据以100×100像元为单元进行聚合，以去除砾石自身造成的地形变化。

经过照片对齐，3个样区对齐的照片数量分别为253张，346张，357张。利用地面控制点对相片对齐的结果进行了优化，得到地面控制点的误差结果如表 1所示。样区控制点总体误差最大的是样区1，误差为28 cm，整体误差最小的样区是样区2，误差为8.19 cm。

通过以DSM作为表面模型进行正射校正，生成了正射图像。3个样区正射图像的地面分辨率依次是6.98，7.61，9.46 mm·像元^-1。根据PhotoScan的结果报告显示，样区1到样区3的正射拼接图像的地面覆盖面积分别为0.022 8，0.034 7，0.038 6 km²。

如图 3内放大图所示，洪积扇表面主要包含砾石、土壤以及植被3种地物类型。大多数砾石颜色为深蓝色，包含少部分黄色、红色、白色砾石。因此在进行砾石训练数据采集时，对各个颜色的砾石都进行了采样，而把土壤和植被的训练数据都归入背景数据集。根据构建的决策树分类模型，分别得到了3个样区的砾石分类图像。如图 3中砾石分类结果图所示，图像中砾石呈白色显示，植被和土壤呈黑色显示，决策树分类模型较好的区分出了砾石和背景像元。

图  3  砾石(粒径≥7 mm)覆盖提取结果

Figure 3.  Extraction map of gravel coverage

分别对3个样区分类后的图像进行了图像分割，按照确定的砾石筛选规则，提取出了形状较为完整且粒径大于32 mm极粗砾。如图 4中放大图所示，提取出来的砾石与正射图像中砾石形状和大小较为一致。

图  4  单个砾石(粒径≥32 mm)提取结果

Figure 4.  Extraction map of individual gravel

应用PhotoScan依据地面点点云生成了样区DTM。3个样区DTM的地面分辨率分别为1.4，1.5，1.9 cm·像元^-1。如图 5所示，样区1的平均海拔为1 735.21 m，样区2的平均海拔为1 436.43 m，样区3的平均海拔为1 208.94 m。整体而言，3个样区的海拔变化均在25 m左右，海拔变化方向都是自东北方向向西南方向递减。

图  5  样区数字地面模型(从左到右分别为样区1、样区2、样区3)

Figure 5.  DTM of sample zones. Sample zone 1, sample zone 2, sample zone 3 is arranged from left to right.

应用ArcGIS 10.2以DTM数据为基础，生成了坡度和坡向数据。如图 6所示，3个样区的地形大体平坦，基本坡度都小于10°。样区1坡度以2°~10°为主，样区2坡度以2°~5°为主，样区3坡度以小于2°为主。如图 7所示，样区1和样区2的坡向都是以西南坡为主；而样区3中各个方向的坡向分布较为均匀。在样区1和样区2中坡向变化还较为连续，而在样区3中由于地形起伏较小，因此坡向变化杂乱。由此可知，从洪积扇扇心到扇缘，地形起伏越来越平缓，冲沟和地垄逐渐消失。

图  6  样区坡度图(从左到右分别为样区1、样区2、样区3)

Figure 6.  Slope distribution map of sample zone. Sample zone 1, sample zone 2, sample zone 3 is arranged from left to right.

图  7  样区坡向图(从左到右分别为样区1、样区2、样区3)

Figure 7.  Aspect distribution map of sample zone. Sample zone 1, sample zone 2, sample zone 3 is arranged from left to right.

经过计算各样区的总体砾石覆盖度(表 2)，样区1的砾石覆盖度为34.22%，样区2的砾石覆盖度为26.85%，样区3的砾石覆盖度为21.88%。通过样区间砾石覆盖度的变化可知，随着海拔下降，砾石覆盖度逐渐下降且下降幅度减小。

如表 2所示，通过统计3个样区的砾石数量，粒径均值、中位数、众数以及标准差，发现：样区1内大砾石较多而砾石数量最少，粒径均值130 mm，中值128 mm，以56 mm粒径砾石最多，标准差最大，说明砾石粒径差异最大；样区2砾石粒径居中而砾石数量最多，粒径均值95 mm，中值85 mm，以61 mm粒径的砾石最多；样区3砾石普遍较小而砾石数量居中，粒径均值78 mm，中值69 mm，以51 mm粒径砾石最多，标准差最小，说明砾石粒径差异最小。

整体而言，从洪积扇扇心到靠近扇缘，砾石覆盖度逐渐降低，砾石粒径均值也逐渐减小，而且变化幅度都逐渐减小，该趋势与海拔变化趋势基本一致，说明在洪积扇尺度上砾石覆盖度与粒径的变化主要受海拔影响。

以南北方向上约30 m为间隔将单个样区划分为样方，计算了每个样方内的砾石覆盖度，砾石Feret最大粒径的均值、中位数、众数、方差。

在单个样区内，如图 8所示，从北向南，样区1内砾石覆盖度先升后降，而样区2和样区3内砾石覆盖度在样区靠下位置上升明显。砾石覆盖度变化与海拔关系不明显。

图  8  单个飞行小区内砾石覆盖度变化

Figure 8.  Changes of gravel coverage in different sample regions

如图 9所示，在单个样区内，自北向南，样区1内砾石粒径均值、中位数和标准差整体呈上升趋势；样区2内，砾石粒径均值和中位数略有下降；样区3内，砾石粒径均值和中位数变化不大。在靠近洪积扇扇心区域砾石粒径相比其他区域变化较大。

图  9  单个飞行小区内砾石粒径变化

Figure 9.  Changes of gravel diameter in different sample regions

整体而言，在单个样区内砾石覆盖度和粒径的变化与海拔关系不大。

基于数字图像来获取戈壁表面的砾石形貌特征参数已被证明具有可行性和适用性^[15-16]，然而地面拍摄的图像往往覆盖面积过小(1 m×1 m)，不利于进行大面积的戈壁表面砾石分布研究。而利用遥感影像估算的戈壁表面砾石粒径受遥感影像空间分辨率的限制，只能估算像元内砾石的平均粒径。利用无人机近地面飞行，可以实时获取较大面积连续覆盖的戈壁表面高清图像，对于分析戈壁面砾石特征具有重要意义。然而受无人机平台本身以及图像的地面分辨率的影响，在获取和分析砾石特征时应当有所取舍。

首先，无人机平台虽然相对于人为调查可以进行更大面积的调查，但是由于砾石相对植被、建筑等地物类型尺寸要小得多，因此无人机的飞行高度不宜过高。受当前无人机电池的续航能力限制，其飞行时间有限，如本研究使用的大疆无人机DJI Inspire 1 RAW，理论载重飞行时间30 min左右，而受戈壁区天气影响，往往实际只能飞行20 min左右。在飞行高度较低且飞行速度较慢情况下，执行一次飞行任务覆盖的地面面积一般小于10 hm²，更适合进行典型样区研究。

其次，当前无人机平台搭载的相机分辨率越来越高，但是因为距离地面较远，获取的拼接正射图像的地面分辨率多在5 mm·像元^-1到10 mm·像元^-1左右(飞行高度距离地面20 m)，相对直接在地面拍摄的照片地面分辨率有所下降，因此不适合识别粒径较小的砾石，但是对于粗砾以上的砾石可以较好的识别。

戈壁表面砾石覆盖度和粒径的分布往往与局部地形间存在联系。Dietze等^[14]通过研究戈壁表面砾石的方位角分布和坡度、坡向的关系，认为坡面过程，如非汇集的坡面流和蠕变流对戈壁的形成和恢复有作用。除了地形的影响之外，戈壁表面植被的分布格局对砾石覆盖度的影响也不可忽略。Quade等^[12]对莫哈维戈壁的研究发现，植被是砾石减少的主要原因。随着植被增多，砾石覆盖度下降明显。

本研究分析了分别位于洪积扇扇心、扇中和靠近扇缘位置的3个样区间以及样区内的砾石覆盖度及粒径变化，发现自扇心到靠近扇缘，3个样区的砾石覆盖度和粒径均值都逐渐下降，但是下降趋势减缓，该趋势基本与海拔下降的趋势相一致。而在各样区内部，砾石覆盖度和粒径的变化与海拔关系不明显，推测可能受到植被分布和局部地形的影响。在不同的研究尺度上，决定砾石特征空间分布的主要因素具有差异。

本研究基于无人机高清图像，并运用砾石形貌特征高效提取技术，对戈壁表面砾石特征进行了研究。通过自洪积扇扇心到靠近扇缘选取3个典型样区，计算并分析了不同样区间以及样区内的砾石覆盖度及粒径变化。研究发现，在不同的研究尺度上，决定砾石特征空间分布的主要因素具有差异。在洪积扇尺度上，砾石覆盖度与粒径的变化整体与海拔变化相一致。而在样区尺度上，砾石覆盖度和粒径分布与海拔关系不大，可能受到植被分布和局部地形的影响。利用无人机和图像处理技术，可以高效且准确地评估大范围戈壁表面砾石特征及分布，为戈壁形成演化研究提供基础数据及技术支撑。