结合FCN和DenseCRF模型的无人机梯田识别方法研究

杨亚男，张宏鸣，3，李杭昊，杨江涛，全 凯

1.西北农林科技大学 信息工程学院，陕西 杨陵712100

2.西北农林科技大学 水利与建筑工程学院，陕西 杨陵712100

3.宁夏智慧农业产业技术协同创新中心，银川750004

梯田是在山丘坡地上沿等高线修筑的阶台式或波浪式断面农田，修建梯田不仅是最快捷有效的水土流失控制措施[1]，还能极大地提高坡耕地的农业生产能力，因此成为黄土高原最重要的水土保持工程[2]。根据水利部发布的《全国坡耕地水土流失综合治理“十三五”专项建设方案》，2016—2020 年专项建设490.71 万亩坡改梯即将完成，及时、准确地掌握区域梯田分布信息成为水土保持监测和评价的重要内容[3]。

传统的梯田信息统计主要通过实地调查实现，但人工统计费时费力。随着高分辨率影像处理技术的不断发展，利用高分辨率遥感影像提取梯田信息成为可能[4-5]。与高分辨率卫星遥感数据相比，利用无人机航摄系统获取高分辨率遥感影像具有高效、及时、精确、低成本的优势[6-7]，在梯田工程的监管和规划等方面具有广泛的应用前景[8]。

无人机高分辨率遥感梯田识别研究中，张宏鸣等[9]基于改进的Canny边缘检测方法通过粗轮廓提取、伪边缘剔除和区域生长算法的结合对梯田进行分割，但是该方法提取精度较低，分类效果一般。Eckert 等[10]采用面向对象分类的方法，结合光谱特征、纹理特征以及形状特征通过SVM 方法对梯田进行识别，该方法仅从区域划分的角度识别梯田建设区域，不能对梯田田面进行精准的识别。薛牡丹等[11]采用面向对象多尺度分割与机器学习结合的梯田提取方法，但该方法的应用效果高度依赖于人工选择的分割尺度，不具有普适性。以上研究方法均依赖人工进行梯田特征的提取，效率不高，并且仅能进行浅层学习[12]，难以适应海量高分辨率遥感影像梯田的识别需求。

自2006 年Hinton 等[13]提出采用深度学习模型实现数据的分类并取得成功后，深度学习因其自动从海量影像数据中学习最有效特征的优势在语义分割、图像分类领域快速发展[14-15]，近年来被广泛应用于遥感地物识别[16-17]和遥感场景分类[18-19]中。陈锋军等[20]以VGG-16模型为基础应用全卷积神经网络分割图像的算法，对无人机采集的云杉图像进行分割，取得了较好的识别效果。夏梦等[21]采用将卷积神经网络与条件随机场相结合的方法，对中国广西南宁市地区的遥感图像进行分类，通过模型的融合在预分类结果的基础上进一步抑制噪声，提高了分类结果的准确性。

目前，基于高分辨率无人机遥感影像梯田识别的方法研究仍然停留在面向对象技术层面，尚未有公开的研究成果将深度学习技术应用于遥感梯田的识别。然而应用于遥感地物识别和遥感场景分类的深度学习方法研究更侧重多种地物类别的区分效果[22-23]，不适用于以精准识别梯田田面为目标的遥感梯田识别研究。

针对现有无人机遥感梯田识别方法的不足之处，本文通过无人机航空摄影获取的高分辨率遥感梯田影像，制作一套梯田正射影像样本集并设计构建全卷积神经网络8s 模型（Fully Convolutional Networks-8s，FCN-8s）与全连接条件随机场模型（Dense Conditional Random Field，DenseCRF）结合的梯田识别方法，实现了梯田特征的自动学习，避免了传统分类方法依赖人工的特征提取过程，并取得了精确的识别效果。

1 数据材料

1.1 研究区域及数据获取

1.1.1 研究区域

中国甘肃省生态环境脆弱，水土流失严重，大量坡耕地既是农业生产的严重制约，也是水土流失的主要发源地。为改善人民生产生活条件，甘肃省榆中地区从二十世纪六七十年代开始兴修梯田，县区内梯田面积较大，范围清晰[24]。

本文的研究区域是中国甘肃省榆中地区龙泉乡黄土丘陵区典型旱梯田区，该区域地理坐标范围为东经104°10′58″～104°19′51″，北纬35°34′4″～35°40′56″，平均海拔2 160 m，总面积84 km2，地形以山沟梁峁为主。该研究区内梯田形态种类丰富，具有一定代表性。图1为该地区所在位置及区域正射影像展示。

图1 研究区域所在位置及区域正射影像展示

1.1.2 数据获取

研究区图像拍摄于2016 年3 月，拍摄期天气晴朗，能见度高，风力小于4 级。使用安翔动力AF1000 型无人机，该机翼展开27 m，质量14 kg，有效载荷2 kg，航高90 m 左右，搭载SONYA5100 相机，有效像素2 430 万，单幅影像面积约340 m×500 m。飞行采用自动起飞/规划航线飞行/自动降落模式，全程总计耗时约24 h。

获取原始数据后采用Agisoft photoscan软件，导入影像，控制点数据及POS数据，将整个区域分成25个测试区块进行处理，并最终拼接生成空间分辨率为0.5 m的正射影像。依据无人机影像制作的1∶500 地形图的精度满足《1∶500 1∶1 000 1∶2 000地形图航空摄影测量内业规范》[25]对1∶500平地、丘陵的成图要求。

1.2 数据预处理

本文从研究区域截取1 块样本区和3 块测试区，其正射影像如图2 所示。样本区（图2（a））为10 000×10 000 像素，范围最大，包含梯田种类全面，特征清晰，作为本实验正射影像样本集来源。测试区1（图2（b））、测试区2（图2（c））、测试区3（图2（d））作为本实验的测试集，均为1 000×1 000像素，包含了研究区内的主要梯田类型。其中测试区1的梯田田块形状较短宽，位于山脊周边，田面无覆盖物，样区内包含大量山丘，部分建筑物及道路，为典型山脊区梯田。测试区2的梯田田块形状舒长较窄，分布密集，田块边缘曲线光滑，田面有部分庄稼堆积物，样区内包含少量建筑物，为密集水平梯田。测试区3的梯田田块形状蜿蜒不规则，田面有大量积雪覆盖，样区内包含部分建筑物和少许道路，属于不规则梯田。

图2 样本区及各测试区正射影像图

本文实验所使用的正射影像数据空间分辨率高（0.5 m）、细节清晰、特征明显、类别丰富，因此基于目视解译使用图像标注工具labelme 软件，制作实验所需要的标签图像，并将梯田田面的像素值设置为1，用红色（rgb（255，0，0））标注，非梯田部分的像素值设置为0，用黑色（rgb（0，0，0））标注。为防止过拟合现象，采用翻转（水平、垂直和左右）和旋转（90°、180°和270°）的方式对数据集进行数据增强。经过数据增强后生成样本集4 000组，其中训练集3 000组，验证集1 000组。实验样本的正射影像（图3（a））和标签图像（图3（b））一一成组对应，如图3所示。

图3 实验样本集图像

2 方法设计

梯田在高分辨率遥感影像上表现为形状相似且排列连续的面状地物，梯田的田坎表现为台阶状等高线条纹理特征，其边界处常有路或沟壑[25]。不同于遥感地物识别与遥感场景分类的多目标类别检测，旱梯田区遥感影像的梯田部分地形特征复杂但有一定规律性，而非梯田部分均质化程度高[27]，因此梯田识别的重点在于通过学习更深层次的特征实现梯田边界更精细的分割和梯田田面更准确的提取。

针对遥感影像梯田识别的特殊性，本文方法首先制作正射影像样本集，构造并训练FCN-8s模型，然后将测试集输入FCN-8s 模型进行预识别，输出的梯田预识别结果作为DenseCRF 模型的一阶势函数，然后利用高斯核函数的线性组合定义DenseCRF 模型的二阶势函数，通过平均近似场方法求解模型，得到最终的梯田识别结果，最后进行精度评价，如图4。

图4 本文方法流程图

方法流程如下：

（1）数据预处理。利用无人机拍摄的高分遥感梯田样本区正射影像，生成图片大小为224×224像素的正射影像数据集，并标注每个实验样本的标签图像，实验样本的标签图像与原始正射影像数据一一成组对应。经数据增强处理后分为训练集和验证集。

（2）构建FCN-8s模型。根据全卷积神经网络原理，在VGG-19（Visual Geometry Group，VGG）网络模型的基础上构建FCN-8s模型。

（3）训练FCN-8s 模型。同时将训练集和验证集放入构建好的模型进行训练及验证，通过参数对比调优，完成模型的训练。

（4）FCN-8s 模型预识别。将测试集输入训练好的模型，经计算后输出FCN-8s梯田影像预识别结果。

（5）DenseCRF 模型。输入梯田预识别结果，构造DenseCRF 模型二阶势函数，计算后输出本文方法的最终识别结果。

（6）精度评价。对比各测试区的真值图，对本文方法的实验结果计算总体精度、F1分数和Kappa系数进行精度评价。

2.1 FCN-8s模型

本文方法中的FCN-8s 模型以VGG-19 作为特征提取器，采用FCN-8s架构进行构建。FCN-8s网络架构是FCN架构[28]中分割结果最细致最准确的网络架构，能够精确预测遥感影像中每个像素的分割结果[29]。VGG[30]是基于Alexnet网络[31]对深度神经网络在深度和宽度上做了改进后形成的具备更强表达能力的卷积网络模型，去掉了Alexnet的局部响应归一化方法（Local Response Normalization，LRN），并采用了更小的3×3 的卷积核、2×2的池化核以及步长为1，因此参数量更少。VGG-19是最深的VGG 网络[32]，由16 个卷积层和3 个全连接层构成，能够更精准有效地提取特征对图像进行分类并得到输入图像的类别概率值。

根据全卷积神经网络的原理，将VGG-19卷积层之后的3 个全连接层转化成3 个卷积层（其中前2 个卷积层卷积核大小为4 096×1×1，第3 个卷积核大小为1 000×1×1）并添加反卷积操作，对特征提取步骤的最后一个卷积层特征图像进行8倍上采样，实现无人机遥感梯田影像的端到端、像素到像素级别的分类，并得到与输入图像大小相同的语义分割结果图。

2.1.1 特征提取

深度学习网络模型的特征提取通过卷积层实现，卷积层利用卷积核计算降维提取图像不同频段的特征。卷积核即滤波器filter，带有一组固定权重的神经元，多个卷积核叠加形成卷积层。卷积层常用的激活函数有tanh、sigmoid 和Relu 函数。但是sigmoid 和tanh 函数都要对输入进行归一化，而Relu函数则不需要，并且Relu函数可以去除数据中的冗余[33]，最大限度地保留数据的特征，因此本文方法中的FCN-8s 模型特征提取器每一层卷积都选用Relu 函数作为激活函数。Relu 函数的数学表达式为：

该特征提取器采用逐层训练的方式，16个卷积层的卷积核均为3×3，每层特征提取的特征图示例如图5。

2.1.2 池化

图5 卷积层特征提取图

池化的目的是对输入的特征图进行压缩，同时改善结果防止过拟合现象的发生。下采样是池化层的主要操作，通过下采样的压缩简化了网络计算的复杂度并在压缩过程中再一次提取特征图的主要特征。池化操作一般有Avy Pooling 和Max Pooling，本文方法FCN-8s模型的池化层采用Max Pooling，池化核是2×2，步长为2，池化操作如图6所示。模型共有5个池化层，pool1在Conv1-1～1-2 后面，pool2 在Conv2-1～2-2 后面，pool3 在Conv3-1～3-4 后面，pool4 在Conv4-1～4-2 后面，pool5 在Conv5-1～5-4后面。

图6 池化图

2.1.3 反卷积

反卷积通过上采样操作从最后一个卷积层得到的特征图恢复到输入图像的原始尺寸，并通过卷积+补位的操作还原多出的空白像素点。本文方法中FCN-8s模型反卷积操作采用转置卷积的方式并添加crop 操作辅助反卷积，最后使用softmax 函数作为分类器对每个像素进行分类，softmax函数的数学表达式为：

其中，yi是softmax 输出的分类结果也称为预测结果，Zi是该网络层的第i个输出，wij是第i个神经元的第j个权重，bj是其bias值。

softmax的交叉熵损失函数表达式为：

本文方法中的FCN-8s 模型结构图如图7 所示。输入样本图像首先在Conv1-1～1-2 卷积后在pool1 层第一次下采样尺寸缩小为原图像的1/2，接着再经过Conv2-1～2-2、pool2 缩小为原图像的1/4，之后从第三组卷积层开始每组4 次卷积，Conv3-1～3-4、pool3，此时缩小为原图像的1/8，并储存pool3 的feature map，继续卷积Conv4-1～4-2、pool4，并储存pool4 的feature map，此时图像已缩小为原图像的1/16，最后进行Conv5-1～5-4 四层卷积，并将pool5 池化后缩小为原图像大小的1/32的图像改称为heat map。Heat map 继续进行conv6、conv7、conv8三层卷积，完成全卷积过程。

完成全卷积的heat map 经过soaftmax 像素级分类生成预测结果图，此时的heat map是原图像大小的1/32，其预测结果图也是原图像大小的1/32。对该结果图进行第一次反卷积2倍上采样恢复到原图像大小的1/16，并与pool4的feature map生成的预测结果叠加，将此结果进行第二次反卷积，经2倍上采样恢复到原图像大小的1/8，再与pool3 的feature map 预测结果叠加后进行最后一次反卷积8倍上采样，最终得到与输入梯田影像大小相同的梯田预识别结果图。

2.2 FCN-8s模型的训练

本文方法的FCN-8s 模型在Ubuntu 16.04.5LTS 系统，Python2.7 环境下，采用TensorFlow1.13 版本进行训练，将训练集和验证集输入网络，设置Batch Size 值为2，一次输入所有数据，设置学习率为0.000 1，训练次数为66 000次，训练过程的交叉熵损失函数值曲线图如图8 所示，其横轴为训练次数，纵轴为损失函数值。前5 000 次迭代过程中损失函数值快速减小，经过30 000次迭代已基本收敛，虽然仍有波动出现，但随着迭代次数的不断增加，波动的幅度和频率都逐渐下降，66 000次的交叉熵损失函数值稳定在6.236 51×10-4。

2.3 DenseCRF模型

图7 FCN-8s模型结构图

图8 交叉熵损失函数曲线

条件随机场（Conditional Random Field，CRF）[34]是图模型的一种特定类型，通过最小化用户定义的能量函数实现模型预测结果的后验分布估计，使图像分割结果更加精准。与由一阶势函数和相邻元素构成的势函数所组成的基本条件随机场（Basic Conditional Random Field，BasicCRF）模型不同，DenseCRF 模型的二阶势函数描述的是每一个像素与其他所有像素的关系，使用该模型在图像中的所有像素对上建立点对势能从而实现极大地细化和分割。

本文方法将FCN-8s模型的预识别结果作为一阶势函数构建DenseCRF 模型，针对全卷积神经网络不能充分考虑全局上下文信息，忽略了空间正则化（Spatial Regularization，SR）步骤，缺乏空间一致性的缺点，通过将每个像素点与其他所有像素点均构成一条边的稠密全连接模型，进一步提高梯田识别结果。

2.3.1 一阶势

DenseCRF 作为一种判别式模型，可以在给定观测场的条件下，对标记场的后验概率直接建模。一阶势函数对单像素点的预测信息和相应的标签之间的关系进行建模，计算每个像素点的类别概率。

随机场的Gibbs分布定义如下：

标签值x对应的Gibbs能量表示为：

这里，P(xi=li)为FCN-8s模型预识别结果中，像素i的标签。

2.3.2 二阶势

式中，μ(xi,xj)=[xi≠xj]为标签兼容函数，用来惩罚相邻相似的像素被标记成不同类型，ω(m)为每个高斯核对应的权重，fi、fj分别表示像素i和像素j具有的特征向量，k(m)为高斯核，公式为：

其中，pi，pj是像素点i，j的位置，θα，θγ是接近度和相似度。

2.3.3 平均近似场算法

本文的DenseCRF 模型使用平均场近似方法[35]求解，计算复杂度为O(N2)，该目标函数的表达式为：

式中，向量fi、fj是像素点i、j的特征向量，ω(m)是线性组合权重，μ是标签兼容性函数，

推断算法如下：

（2）迭代更新：

（3）当算法收敛时，停止迭代。

2.4 精度评价

精度评价是指比较实际数据与方法结果以确定该方法的准确性。本文通过目视解译，在arcmap10.5中绘制接近于实际梯田的真值图作为评价样本。为了做出有效的评估，实验使用总体精度（Overall Accuracy，OA）、F1分数和Kappa系数[36]对本文方法的实验结果进行量化评价[37]，三种精度的计算方法为：

式中，TP（True Positive）为分类正确的梯田像素点；TN（True Negative）为分类正确的非梯田像素点；FN（False Negative）为分类错误的非梯田像素点；FP（False Positive）为分类错误的梯田像素点。

3 实验结果与分析

为验证本文方法的有效性，分别选取山脊区梯田、密集水平梯田和不规则梯田正射影像作为测试样区，并在使用相同样本集的基础上将本文方法的最终实验结果分别与FCN-8s梯田影像预识别结果，SegNet、DeepLab方法结果进行对比。图9 为三个测试区梯田识别结果对比图。

图中第一行是测试区1 山脊区梯田的识别结果对比图，识别的误差主要出现在以红框1为例的不平整梯田田面和蓝框2、黄框3 为例的邻山脊处梯田。其中DeepLab方法在红框1处断裂点最多，识别效果最差，蓝框2、黄框3 处，由于紧挨山脊，纹理特征和颜色与其他部分田面有所不同，均错分为非梯田。SegNet方法在红框1处识别效果相对于DeepLab方法有所提升，但依然有明显断裂，并且在蓝框2 和黄框3 处也出现错分情况。以上两种方法对于梯田遥感影像的纹理颜色特征容错性较差，在山脊区梯田识别中表现不佳。FCN-8s预识别结果图中，红框1处不平整田面部分识别效果明显提升，蓝框2、黄框3 处略有提升，部分田块边缘像素点识别正确，但未能识别整块田面；在本文方法的最终实验结果图中，红框1处和FCN-8s识别效果相当，蓝框2、黄框3 处识别效果明显提升，可以正确识别整块田面，识别效果较好。

图中第二行是测试区2 密集水平梯田的识别结果对比图，识别的误差主要出现在以红框4为例的建筑物群和蓝框5、黄框6为例的土坎陡坡地形处。其中Deep-Lab 方法在红框4 处将大部分建筑物错分为梯田，在蓝框5 和黄框6 处将土坎陡坡也错分为梯田。而SegNet方法在红框1 处识别效果略有提升，错分情况有所好转，但是在蓝框5 土坎陡坡处也将土坎陡坡错分为梯田，并且将蓝框5 左侧的非梯田部分错分为梯田，在黄框6 处土坎陡坡错分为梯田的部分相对DeepLab 方法有所减少，但还是明显存在。FCN-8s 预识别结果图中红框4 处建筑物群基本都能正确识别为非梯田，蓝框5和黄框6 的土坎陡坡也正确识别为非梯田，但是蓝框5和黄框6之间的梯田区域都错分为了非梯田；本文方法的最终实验结果图中，红框4处建筑物群基本能正确识别为非梯田，蓝框5 和黄框6 处大部分土坎陡坡正确识别为非梯田，两处之间的梯田区域错分情况也有所减轻，相比于前三种方法，取得了最好的识别效果。

图中第三行是测试区3 不规则梯田的识别结果对比图，由于该样区内梯田田面有积雪覆盖，积雪覆盖区域纹理平整，颜色接近，识别效果误差较大，误差主要出现在以红框7 为例的断崖下的梯田，蓝框8 为例的积雪未完全消融的湿土梯田田面，黄框9为例的纹理多样的不规则小块梯田区。其中DeepLab 方法和SegNet 方法在红框7处、蓝框8处均出现错分情况，在黄框9处，Seg-Net 方法的识别效果比DeepLab 方法略好，但识别效果依然不够理想。FCN-8s 预识别结果图中红框7 处依然没有能正确识别断崖下梯田，蓝框8 处和黄框9 处田面的识别效果相比于DeepLab和SegNet方法有所提升，但是比较粗糙；本文方法的最终实验结果图中，红框7 处能够正确识别出断崖下梯田，但是断崖部分也一并识别为梯田，由于该影像来自于俯视航拍，近于垂直的断崖不能拍摄到完整断崖面特征，而拍到的断崖顶颜色纹理和断崖下梯田相似，从影像上难以区分。蓝框8处和黄框9 处基本能正确识别，其中蓝框8 积雪未完全消融的湿土田面识别效果较为精准，黄框9处多块不规则小块梯田由于纹理颜色各异，识别效果不够精准，有待进一步提升。

从实验结果图中可以看出，本文方法的识别结果最接近于人工标记的真值图像，Deeplab 方法和SegNet方法在测试区1、2、3的识别中错分和漏分现象明显，效果难以令人满意。FCN-8s 模型预识别效果在测试区1的识别中表现良好，而测试区2 和测试区3 的识别结果相对较为粗糙，结合DenseCRF 模型优化处理的最终结果在3个测试区的识别结果更细致准确，具有良好的鲁棒性。

图9 测试区梯田识别结果对比图

表1 测试区梯田识别结果精度评价%

此外，根据测试区正射影像图分析本文方法的误差，测试区1中紧挨山脊处的梯田田面颜色与山脊相近且与非邻近山脊的其他梯田差异较大，容易被错识别为非梯田；测试区2 中土坎陡坡处颜色、纹理和轮廓特征与邻近的梯田田面十分接近，只是因为坡度变化大且坡平面较窄不利于作物种植所以一般不在土坎陡坡处修建梯田，目视解译时凭借人工经验将此处划分为非梯田，如何进一步解决梯田识别中这种典型的同谱异物现象有待在今后的研究中改进。测试区3 由于受田面积雪覆盖影响，梯田田面之间颜色和纹理等特征相差较大，进一步提升田面有覆盖物情况下的梯田识别也要在之后的研究中进行探索。

表1为测试区梯田识别结果精度评价表。从表1可以看出，相比于本文方法，Deeplab 方法和SegNet 方法在总体精度、F1分数和Kappa系数三种精度的表现都比较差。

与FCN-8s 梯田影像预识别结果相比，本文方法将FCN-8s与DenseCRF模型结合后，三个测试区在总体精度、F1 分数和Kappa 系数表现上平均分别提高了5.41、7.29、10.07个百分点。

结果表明，本文方法在三项精度指标上表现最佳，平均总体精度达到86.85%，平均F1 分数达到87.28%，平均Kappa系数值达到80.41%，取得了较好的识别效果。

4 结束语

本文针对无人机高分辨率遥感梯田识别研究中无法进行梯田特征自动学习的问题，根据梯田识别的特点，制作一套梯田正射影像样本集，提出一种FCN-8s与DenseCRF模型结合的深度学习识别方法。

该方法首先在TensorFlow 深度学习框架上搭建FCN-8s网络架构，以VGG-19深度神经网络模型作为特征选择器对梯田影像进行预识别，通过模型中的卷积层在梯田正射影像样本集上完成梯田特征的自动学习，避免了传统提取方法中复杂的人工处理过程。该步骤将FCN全卷积网络中分割最细致的8s架构与VGG系列最深的网络结合进行特征提取与预识别，充分利用VGG-19深层学习特征的优势以进行梯田影像二分类的精准区分。然后以预识别的分类结果定义每个像素点一阶势函数，用两种高斯核函数的线性组合定义二阶势函数，利用全连接条件随机场考虑丰富的空间上下文信息的特点，抑制更多的分类噪声，进一步提高梯田识别精度。实验结果表明，该方法取得了较好的梯田识别效果，推动了深度学习技术应用于无人机遥感梯田识别方法研究领域。

本文仍有可以改进之处，后续的研究将从以下几个方面进行尝试：（1）能否改进模型以适应多源数据的输入和学习，通过将梯田的地形因子数据与影像数据结合学习提高识别精度。（2）能否将集成学习与深度学习方法结合，通过集成学习的投票技术进一步优化梯田识别结果。（3）在深度学习技术之外，能否采用半监督学习方式进行无人机遥感梯田的识别有待进一步探索。