基于深度学习与多源遥感数据的新增建设用地自动检测*

张泽瑞，刘小平，张鸿辉，罗伟玲

1. 中山大学地理科学与规划学院，广东 广州 510275

2. 广东国地规划科技股份有限公司，广东 广州 510650

随着我国经济发展和人民生活水平的提高，城镇快速扩张，常住人口城镇化率从1949 年的10.6%增长到2018 年的59.58%［1］，建设用地从2001 年的3 641.2 万hm2到2016 年末的3 906.82 万hm2［2］。建设用地的快速扩张下伴随着大量的违法占地现象，如非法占用土地、破坏耕地和破坏森林资源，其中新增建设用地是违法用地中最常见的类型。仅2016 年间，全国各级国土资源主管部门查处了涉及2.69 万hm2的共7 万余违法用地案件［3］。2020年1月，自然资源部发布《自然资源调查监测体系构建总体方案》，方案中指出以遥感监测为主要手段的技术体系为重点，查清我国各类自然资源状况和变化情况。因此，自动检测新增建设用地将为违法用地的查处提供一种新的思路，可以成为保护我国自然资源的重要技术手段。

早期的违法用地查处与监督主要依靠人工实地检验，不仅费时、费力、效率低，而且难以做到全面的查验、及时的发现和准确的划定。随着遥感技术的成熟与广泛的应用，已有学者利用遥感影像实现了违法用地的检测［3-6］。例如，曹端广等利用无人机倾斜摄影技术完成了对违法用地监测系统的研究［7］，余永欣等基于地理国情数据进行了违法用地检测［8］，但依旧离不开大量的人工判别的方式。因此，陈鹏等基于高分辨率遥感影像提出了用于违法用地检测的变化检测方法［9］，有效地提高了土地监察执法的效率；孙笑古等在土地督查地理信息服务平台中提出了图像分割、信息提取算法以提高违法用地的检查效率［10］。上述研究虽利用了遥感技术来获取违法用地的信息，但未能发挥出多源遥感影像的优势，实现大范围、快速地检测。

人工智能中的深度学习技术是通过深层次的神经网络，从给定的目标任务中获得目标数据源由低层到高层的特征，提升了语音识别、视觉对象识别、目标检测和其他众多领域的技术水平［11］。在图像应用中，主要的应用方向是目标检测、人脸识别、图像分割、姿态估计、图像分类［12-13］，其中深度学习与遥感影像相结合，已经应用于土地利用与土地覆盖分类、目标检测、场景识别、影像分割、变化检测、影像融合、影像配准等各个方面。实现了对特定地物类别、场景类别的目标识别、语义分割、分类等多项应用［14］。例如，深度学习在建筑物提取中克服了传统方法的其他地物的误分类与误提取问题，提取精度得到了提升［15］；在土地利用与土地覆盖分类任务中，融合深度学习网络提取的特征在识别农田、灌木地、草地、不透水面、湿地等地类中比以往的机器学习模型得到了更优分类精度［16］。但在违法用地查处与管理中，深度学习尚未得到很好的应用。因此，发挥深度学习模型的优势将为新增建设用地自动检测提供一种新的方法。但至今仍缺少将多源遥感数据与深度学习技术相结合，应用于新增建设用地自动检测中的研究。

针对以上问题，本文提出了在新增建设用地变化规则的基础上对数据进行标记，对深度学习模型进行训练、调优，并结合Sentinel-2A 遥感影像的变化区域对新增建设用地进行大范围的快速检测，及时发现用地违法现象，提高建设用地管理效率，以利于自然资源调查和国土资源维护。

1 研究区概况及数据

1.1 研究区概况

本文研究区位于广西省贵港市，位于北纬22°39 ′～24°02′、东经109°11′～110°40′之间。总面积10 602 km2，市中心城区建成区面积73 km2，全市总人口约555 万。在2018 年贵港市政府工作报告中，2017 年贵港市共获得新增建设用地指标约2 000 hm2，同比增长117%。同年，贵港市有关部门开展违法用地与整治违法建筑的行动，查处违法占地面积42.69 hm2，违法建筑面积34.86 hm2。本文选择广西贵港市作为研究区域，并将研究区域划分为训练区与测试区（图1）。

图1 贵港市研究区划分Fig.1 The study area division of Guigang city

1.2 数据来源与预处理

本文使用的原始遥感影像数据包括：1）2018年1 月与2019 年1 月贵港市高分辨率可见光波段遥感影像，空间分辨率为2 m。高分辨率遥感影像的预处理以时相一影像为基础对时相二期影像进行处理，根据两期影像中出现的问题，具体预处理内容包括：地理配准，直方图匹配，拉普拉斯算子滤波操作，使得两幅影像在纹理特征与光谱特征上一致（如图2）；2） 2017 年12 月与2018 年底贵港市Sentinel-2A 星遥感影像，Sentinel-2A 星影像从可见光到短波红外共有12 个，空间分辨率10～60 m 之间，此次研究选择10 m 分辨率的可见光波段和近红外波段。Sentinel-2A 星遥感影像下载来源于美国地质调查局（https：//earthexplorer. usgs.gov/）。该网站上所提供的是经过几何精校正的大气表观反射率产品Sentinel-2A L1C，通过欧空局提供的Sen2Cor 插件可对L1C 产品进行大气校正，生成L2A地表反射率产品。

图2 训练区高分辨率影像预处理Fig.2 High resolution image preprocessing in training area

2 研究方法

首先对训练区中出现的新增建设用地在高分辨率影像上进行标记，对标记的训练数据进行数据增广和差值计算操作，选择深度学习中的Deep-Labv3+模型进行训练并调优；为达到快速、大范围的新增建设用自动检测的目的，利用迭代加权多元变化检测（IRMAD，iteratively reweighted multivariate alteration detection），对两期Sentinel-2A 影像进行变化检测，再用高斯混合模型（GMM，Gaussian mixture model）对变化检测结果进行阈值划分，得到变化区域结果；最后，在变化区域中选择高分辨率遥感影像在训练后的模型上实现新增建设用地的自动检测并评定精度，技术路线如图3。

图3 技术路线Fig.3 Technical route

2.1 训练数据处理

2.1.1 训练样本标记新增建设用地在不同时刻的遥感影像上表现出来的图斑变化可以分为3 种：①植被转变为裸地、空地；②植被转变为建设用地；③空地、裸地转变为建设用地。基于上述3种新增建设用地的变化规则，对两期高分辨率遥感影像利用人工目视判别进行标记。此外，为区分农业生产活动导致的周期性植被变化，需要重点考虑变化规则一，本文采取下述方法对变化规则一进行筛选：①若裸地变化明显且裸地周围或内部存在建设用地，则勾选该裸地图斑；②若裸地变化不明显，但裸地内部存在建设用地，仅勾选内部建设用地图斑；③若裸地变化明显，但周围不存在建设用地，则不勾选该裸地图斑。图4展示了3种在训练区中图斑变化规则的新增建设用地标记。

图4 训练区新增建设用地标记（红色标记）Fig.4 Labels of newly increased construction land in training area(marked in red)

2.1.2 数据增广与处理对两期高分辨率影像以特定窗口大小和重叠率进行裁剪得到原始训练数据集，考虑到数据大小有限，而深度学习中需要依靠大量的数据来训练以提高模型的拟合能力，需通过数据增广操作对数据进行扩充，数据增广操作的目的是提高模型的泛化能力减少过拟合［17］，传统的数据增广方法是利用现有数据对空间域或色彩域进行变换等操作来生产新样本，如旋转、平移、变形、缩放、裁剪、色彩抖动、对比变换、噪声等一系列操作，其他的数据增广方法还包括生成式对抗网络GAN （generative adversarial networks）生成样本、幻觉图像生成样本、基于特征空间的数据增广等［18］。此次研究中，仅对裁剪后的影像进行镜像翻转、旋转、缩放3种数据增广操作。

对于在两期高分辨率遥感影像上的进行新增建设用地的检测，本质上是变化检测的一种，传统变化检测中对于双时相数据的处理主要有：基于构造差异的方法（差值法、比值法、差异主成分法、CVA）、基于变换特征间差异的方法（K-T变换法、主成分差异法、MAD）、基于分类后比较的方法［19-22］。本次研究仅利用双时相影像进行差值处理，通过以下方法对增广操作后的数据进行简单的差值处理，

其中X，Y分别为时相一和时相二的可见光波段高分辨率影像，差值图D是两期影像各位置像元相减再利用根据像元可能的取值范围归一化至［0，1］。

2.2 DeepLabv3+模型

深度学习中图像分割任务的目标，是根据图像的纹理、颜色、形状、领域关系等特点将具有相同特征的像素点进行聚类，达到像素级别的分类效果。经典的图像分割模型有：FCN［23］、UNet［24］、SegNet［25］、DeepLab 系列［26-30］、RefineNet［31］、PSPNet［32］。DeepLab 系列如今共有v1［26］、v2［27］、v3［28］、v3+［29］、auto［30］等5 种版本，其中v3+继承了前3 种版本的空洞卷积、ASPP 模块、批量正则化层的优点，同时引入了SegNet 中的Decoder-Encoder 结构，分割结果有了更好的表现。图5 为本文使用的DeepLabv3+模型结构。

图5 DeepLabv3+模型结构Fig.5 DeepLabv3+model structure

DeepLabv3+模型以公式（1）构建的差值图作为输入，首先用主干特征提取网络（以图中的ResNet 为例）提取出下采样4 倍的底层特征(64×64×2 048)与下采样16 倍高层特征(16×16×2 048)，对高层特征分别进行不同步长的ASPP（Atrous spatial pyramid pooling）和全局池化操作后进行连接，并通过卷积上采样4 倍到与底层特征相同的尺寸。将卷积上采样后的特征与底层特征进行拼接，再次经过卷积与上采样4倍恢复原始尺寸并得到预测结果。DeepLabv3+模型框架的优点在于将模型划分为Encoder-Decoder 结构，Encoder 结构中对高层特征进行ASPP操作，ASPP通过不同大小空洞卷积与全局池化得到不同尺度下的语义信息；Decoder结构中，结合高层特征提供的语义信息与底层特征逐步恢复空间信息得到预测结果。

2.3 变化区域提取算法

2.3.1 IRMAD变化检测为实现快速、大范围的新增建设用地检测而忽略两期影像间变化不明显的区域，本文采用利用IRMAD 变化检测算法提取Sentinel-2A 影像上变化区域。IRMAD［33］算法是在MAD［34］基础上提出的。最初的MAD 算法是对影像波段进行线性变换得到两个特定变量，变量间相关性反映了原始影像波段之间的相关关系的方法，具体计算公式为

由公式（5）可知，U-V之间的方差越大，则它们的相关系数越小。但对原始影像进行变化后所得的U-V的相关系数是递减的，所以对所求结果按照相关系数从小到大排列，即按照方差从大到小排列得到MAD算法的各个分量

其中不同的MADi表示影像中不同特征间的关系。IRMAD 算法是在MAD 算法基础上引入迭代加权思想，通过每次迭代计算像元的权重用于更新下一轮中平均值与方差的计算，提高变化检测精度。每个像元的初始权重为1，假设对每个像元j的MADij进行标准化到单位方差后的平方和是服从自由度为p的卡方分布χ2(p)，

公式（7）描述了通过像元j的各个MAD 分量与对应的方差来计算像元j在卡方分布中的位置Tj。通过计算得到Tj后，再根据公式（8）确定像元j在卡方位置上百分位值，查表可得像元j在下一轮的计算中的权重

当迭代轮次完成或者当前轮次与上轮次相关系数的绝对差值中的最大值小于阈值后停止迭代，本文所用迭代阈值为10-6。

2.3.2 GMM阈值划分变化检测中需确定最后的变化阈值对结果进行划分，确定变化与非变化区域。本文利用GMM［35］对IRMAD 算法的变化检测结果进行自动划分阈值，减少人工选取阈值的干扰。GMM 模型是将数据假设由几个不同的正态分布的随机变量组合而成，考虑不同随机变量之间的协方差以及潜在的正态分布信息，

其中w1+w2= 1，x表示的是变化检测后的结果，w1和w2分别是变化类与非变化类的权重，g(x|u1，Σ1)和g(x|u2，Σ2)分别表示变化类与非变化类服从高斯分布的概率密度函数，通过最大似然函数对参数：λ={w1，w2，u1，u2，Σ1Σ2}进行估计得到各参数值。

3 研究结果

3.1 模型训练

本文实验的硬件环境为Intel(R) Core(TM)i7-8700K 3.70GHz 的CPU 与NVIDIA GeForce GTX 1070 Ti 8GB 的GPU； 软件环境为Python3.6、PyTorch 1.6.0。实验中采用256×256大小的窗口对两期高分辨率影像进行裁剪，对裁剪后的数据进行数据增广操作和差值图处理后训练模型。此外，本文中对不同的主干特征提取网络进行实验，分别选择VGG16［36］、ResNet152［37］、DenseNet161［38］、GoogLeNet［39］、MobileNetV2［40］等5 种经典图像分类网络，以上网络均使用在ImageNet 上预训练的权重以加速模型收敛。精度评定选择IoU、总体精度和F1 分数等3 种指标，训练参数均为0.000 1 的初始学习率、100 轮训练轮次、损失函数为Binary Cross Entropy、Adam 优化器、学习率衰减策略为每10轮次下降0.5、采取早停法进行训练。不同主干特征提取网络下，在测试集上得到的精度结果见表1。

表1 不同主干特征提取网络精度Table 1 The accuracy of different backbone feature extraction networks

从3种指标精度上表现来看，各类主干特征网络在总体精度上相差不大，都在97%左右；F1 分数也较高，分布范围在85.14%～90.57%之间；对于IoU 指标，最高精度与最低精度则相差8.1%。相对于其他网络，DenseNet161 作为特征提取网络不仅模型参数量适中而且能在3种指标上取得较好的结果，所以选择DenseNet161 作为DeepLabv3+模型的主干特征提取网络作为后续研究中新增建设用地自动检测的模型。

3.2 变化区域提取

变化区域提取结合了IRMAD 算法和GMM 算法，在两期影像上可自动得到变化与非变化区域，并用形态学操作与NDVI 掩模进行后处理优化提取结果。首先对两期Sentinel-2A 影像经过预处理后，利用IRMAD 算法对可见光波段与近红外波段进行变化检测。由公式可知，IRMAD 算法能得到4 个MAD 分量的结果，图6 展示了由IRMAD 算法得到4 个分量结果。其中每个分量分别对应不同特征成分的变化，变化程度则代表了该特征下地类在两时相内发生变化的大小。通过各个分量与影像上变化的地类进行对比，分量二与所期望提取的变化成分相似度最高。相对比其他分量，分量二中包含了较为全面的地类变化的信息，但所含的不仅是从植被变化到裸地、建筑用地等地类变化，也包含了裸地到植被的转变，这种情况将为检测带来的非必要的计算量，本文还利用时相二影像的NDVI 进行掩模操作，对分量二中结果NDVI 大于0.5的像元将它重新归为未变化像元，最后使用形态学操作中的开操作消除细小的像元得到变化区域，得到最终的二值变化检测结果区域。

图6 变化检测分量结果Fig.6 The results of change detection component

Sentinel-2A 影像提取的变化区域用于判断该区域是否存在变化像元，即当前待检测区域中若变化区域的像元数量占比大于该区域总像素数量的1%，则对该区域进行检测，否则不进行检测。因此本文以像素为单位，对不同大小的矩形结构元素进行实验，在测试区中对参考新增建设用地图斑进行漏检率和模型计算量减少率进行分析，结果如表2。

表2 不同结构元素对计算量和漏检率的影响Table 2 The influence of different structural elements on calculation amount and omission rate

结构元素越大，消除更多的细小像元的同时能减少更多的计算量，但也会排除一些变化较小但仍有新增建设用地的区域。综合考虑漏检率与计算量减少率，选择大小为3的矩形结构元素对变化区域进行开操作处理，图7展示了对变化检测结果进行未处理的变化区域、NDVI 掩膜处理后的变化区域、开操作处理后的变化区域。

图7 变化检测结果处理Fig.7 The results of change detection before and after processing

3.3 新增建设用地检测

本文结合变化区域与训练调优后的模型，在测试区中实现新增建设用地自动检测，并对结果进行精度评定。将测试区的影像按照模型训练窗口大小（256×256）与一定的裁剪重叠率进行划分，若划分的影像中变化区域的像素大于区域总像素的1%，则对该划分区域的两期影像计算按照公式（1）计算差值并进行检测。为了使得检测结果与实际影像相符合，采用了条件随机场与形态学闭操作对检测结果进行边界优化、连接相近图斑和去除图斑内部空洞等后处理方法。

新增建设用地自动检测精度评定过程以建设用地图斑为基本单位进行精度评价。由于参考图斑中最小面积为1 004 m2，所以仅将新增建设用地检测结果中面积大于1 000 m2的图斑与测试区中的参考图斑进行检验，若检测图斑与参考图斑相交则标记为正确检测。表3 统计了经过后处理方法后，不同裁剪重叠率下新增建设用地检测结果中划分区域数量、图斑正确检测率、漏检率等。

通过表3 可以看出，随着重叠率的逐渐增加，新增建设用地检测图斑的正确率上升，但是划分区域数量与错分率也随之增加。在70%和80%与80%和90%重叠率之间，划分区域数量有明显的跳跃式增长，虽然80%与90%重叠率的检测正确率上能达到很好的精度，但划分区域数量和错分率明显高于其他裁剪重叠率。相对于其他裁剪重叠率，70%的裁剪重叠率在计算量、检测正确率和错分率上能得到很好的平衡，所以选择70%裁剪重叠率的检测结果用于后续实验。图8中展示了测试区中新增建设用地的图斑检测结果。

表3 不同裁剪重叠率的新增建设用地检测结果Table 3 The detection results of newly increased construction land with different cropping overlap rate

70%裁剪重叠率下的新增建设用地检测结果中，图斑检测正确率为85.16%，图斑漏检率为16.73%。从图8 中可以看出，变化明显的新增建设用地检测结果较好，面积较大的用地在城区周围分布多，面积较小的用地有少数分布在郊区。图斑误判数量较高，错分率为36.57%，大部分错分图斑分布在城市周围或远离城市的区域，错分图斑虽然能检测出被破坏的植被，但此类植被被破坏并不是因为建设用地的建造而是农业生产活动导致的。对图8 的测试区中的检验结果进行统计，90%的检测图斑IoU 精度大于0.2，其中有62.68%的检测图斑IoU大于50%。在检测的图斑中，IoU最大的为95.77%，最小的为1.57%，图斑平均IoU为57.23%，参考图斑总面积795.06 hm2，检测图斑相交面积592.53 hm2，检测面积率74.52%。

图8 测试区中的检测结果Fig.8 The detection results in test area

4 结 论

本文对现有的违法用地检测方法进行了梳理，针对当前对违法用地管理与研究不足的现状，提出了一种基于深度学习和多源遥感影像的新增建设用地自动检测方法。以广西贵港市为研究区，在高分辨率遥感影像上对新增建设用地图斑进行标记、处理、训练深度学习模型，对模型的不同特征提取网络进行分析，选择DenseNet161 作为DeepLabv3+的主干特征提取网络；将高分辨率影像与Sentinel-2A 影像的变化区域相结合，实现快速定位变化区域、减少计算量的目的；此外，研究中还分析了开操作中不同参数量对图斑漏检率和模型计算量的影响，最后还针对不同的裁剪重叠率对图斑检测精度的影响进行分析。结果表明，70%的裁剪重叠率能在计算量与图斑检测精度、错分率之间能取得较好的平衡，在测试区中图斑检测正确率85.16%，错分率36.57%，图斑平均IoU 为57.23%，62.68%图斑的IoU 在50%以上，总体面积检测率为74.52%。本文的研究利用多源遥感影像的优势与深度学习的特点，在同类研究中能达到快速、大范围对违法用地进行自动检测的目的，在一定程度上减少了人工目视判断与边界划定的工作量，为城市规划、土地管理提供一种新的工作手段，在一定程度上提高违法用地的监测准确度及工作效率。

本文的研究内容是辅助违章用地执法监察的一次成功实践，虽然此类方法能检测出大部分正确的图斑，但仍有较高的错分率，整体IoU 过低。导致错分图斑产生的主要原因是错分图斑的变化属于农业生产活动导致的植被周期性变化，虽然在样本标记环节中对符合此变化规则的图斑进一步筛选，但在测试区中仍未能很好地解决此问题。此外，在后续实验中可以增加对其他类型向建设用地变化的检测，如水体变裸地或水体变建筑等。由于实验中时相二高分辨率影像由多幅影像拼接而成，在部分区域中与时相二影像在光谱、纹理方面差异较大，处理两期影像的方式也是较为简单的差值相减的方法。在后续研究中，使用质量相似的两期高分辨率影像能进一步减少错分率、提高图斑IoU，并可结合孪生神经网络结构对新增建设用地自动检测进行进一步研究。