基于遥感影像的水域违章建筑识别研究

陈建国 ，王晨辉 ，徐绪堪 ，嵇庆才

（1. 宁夏水利电力工程学校，宁夏 银川 750006；2. 河海大学商学院，江苏 常州 213022；3. 常州工业大数据挖掘与知识管理重点实验室，江苏 常州 213022；4. 镇江新区城乡建设局，江苏 镇江 212132）

0 引言

随着城乡建设的快速发展，城市现代化进程不断加快，大量河湖违章建筑影响城乡水域公共空间，河湖流域违法的圈圩和建设严重影响城市防洪排涝，给城市社会经济稳定发展带来较大潜在危害。城市河湖违建整治是城市全面推进河（湖）长制、建设人水和谐生态的重要举措之一，根据江苏省委办公厅要求，全省各地对辖区内河湖流域逐一排查，确保各类违法圈圩和建设不遗漏，2020 年底将全面整治到位。

在过去的河湖“清四乱”专项活动中，水域违章建筑检测主要依赖人工对 0.3 m 分辨率遥感影像进行比对排查，为解决人力资源有限、人工比对耗时长的问题，提出一种基于遥感影像变化技术的水域违章建筑识别模型，并通用实例验证识别的效果，为水域违章建筑检测提供科学快捷的方法。

1 相关研究

众多学者在违章建筑检测方面提出了很多方案，詹金瑞[1]提出了以历年卫星影像为基础，结合实地测绘和规划信息的违章建筑检测总体架构；朱建伟等[2]提出利用无人机低空环拍影像构建三维模型进而识别违章建筑的系统；黄蓉等[3]以青岛市为例，设计了基于高分辨率遥感影像变化检测技术的违章建筑检测技术路线，解决了实地测绘的低效率和无人机环拍的高技术门槛问题，成为了目前主流的违建检测方法。

遥感影像变化检测技术是指通过图像处理等技术对不同时期同一区域的遥感影像数据进行分析以确定 2 个时期影像间的变化[4]。依据处理单元粒度大小不同，可以将现有的变化检测算法分为两大类：1）基于像素点的变化检测。以像素为最小单位进行处理，试图将所有像素归为发生和未发生变化2 类的方法，例如：Malila 最早提出变化矢量分析法（CVA）[5]，利用多波段的遥感影像数据将每个像素描述为 1 个一维列向量，进而计算前后时期同位置像素之间的向量差值，代表变化的强度；陈晋等[6]将 CVA 算法用于国内土地利用/覆盖变化检测，并提出一种双窗口变步长阈值搜寻的方法，在变化检测的自动程度上更进一步；黄维等[7]结合主成分分析（PCA）和 CVA 算法，对多波段影像数据提取第一主成分后进行差值运算和阈值划分，降低了图像噪声的影响；赵秋菊[8]进一步改进 PCA 算法，以差值图像的一倍标准差作为变化阈值，自动确定出变化区域。这类算法一般对输入影像向量表示并求差值后，人工设定或自动搜索阈值，确定每个像素是否发生变化。2）面向对象的变化监测。常常采用先分类后检测变化的方法，首先对某一类型像素集合看作对象，研究改对象范围内是否发生变化。例如：李亮等[9]通过影响分割获取像斑，结合 EM 算法和贝叶斯最小错误理论确定每个像斑的变化/非变化类别归属；张亚亚等[10]提出一种基于对象的语义网遥感图像知识分类框架通过制定规则解决对象分割问题。近年来，随着计算机运算性能的提高，基于深度学习的监督分类方法往往在实际中取得更好的效果，张晓东等[11]将目标检测领域的主流网络Faster R-CNN 应用到高分辨率遥感图像变化检测中取得了理想的效果。

综上所述，国内外学者为违章建筑检测贡献了有效的系统架构和技术方法，然而传统的基于像素的变化检测算法容易出现“椒盐”现象，难以满足违章建筑粗筛选的精度需求，而基于深度学习的运算时间高度依赖计算机的性能，难以适用实时的检测系统。因此，本研究提出基于改进主成分分析和k-Means 聚类的遥感影像快速变化检测算法以实现水域违章建筑粗筛选。

2 遥感影像快速变化检测算法设计

2.1 问题定义

给定当前时期某区域遥感影像，欲筛选出违建区域，首先获取过去某时期同区域遥感影像；接着将 2 个时期影像数据代入算法，获得变化区域图斑，并在当前时期遥感影像中标记出变化位置；最后将标注有变化区域的遥感影像推送至核查人员处，审查变化区域是否发生违章建筑。其中，获取变化区域是核心算法，X1，X2是同一地区不同时期的 2 张遥感影像图片，满足以下条件：

式中：H和W分别表示图像的高度和宽度；x1(i，j)和x2(i，j) 表示像素点 (i，j) 的灰度值。

目标是生成一张变化区域图片Xc= {xc(i，j) | 1 ≤i≤H，1 ≤j≤W}，满足：

2.2 算法流程

算法流程主要包括计算差分图像、构造图像变化特征矩阵和像素点二分类 3 步。1）差分图像通过2 个时期遥感影像灰度值差值计算；2）图像变化特征通过 2 次划分像素块，将差分图像的每个像素点变化信息表示成向量，并利用 PCA 算法降维；3）像素点二分类是利用改进的k-Means 聚类将图像像素点分为发生和未生变化 2 类。算法流程如图 1 所示。

图 1 算法流程图

2.2.1 计算差分图像

在变化检测模型中有 2 种主流方法：合并 2 个不同时间的图像；用后一时期影像减去前一时期影像得到差分图像。

本研究提出的快速检测的算法是基于差分图像的，差分图像Xd表示图像X1和X2每个像素点灰度值的绝对差值。

2.2.2 构造图像变化特征矩阵

将差分图像Xd划分成大小为h×h(h≥2) 的像素块，则

式中：xd(m，n) 表示位于m行n列的像素块，xd(m，n) 简化写作将h×h大小的像素块依据行列顺序展开为h2×1 的向量表示第k个像素块。

式中：K为划分后像素块的总个数每个向量距离平均向量Ψ的距离

PCA 算法试图寻找一组N(N≤h2) 个线性无关的综合指标代替原来的h2个指标，对集合Δk使用PCA 算法，首先计算协方差矩阵其中i和j的协方差cij= Cov(Δi，Δj)，i，j=1，2，…，h2。然后计算协方差矩阵C的特征向量与特征值，由于矩阵C= (cij)h2×h2，所以共有h2个特征向量es与特征值λs，按照特征值递减顺序排列，选取前N(N≤h2) 个特征向量组成特征向量空间EVS= [e1，e2，…，eN]T，N≤h2。

最后重新将差分图像Xd划分成大小为h×h(h≥2) 的像素块，这次划分采用重叠的方式，即同一行或同一列相邻的 2 个像素块之间的间隔为定长，用stride表示，设置每个stride为 1 个像素，划分后Xd={xd(m，n) |m+h=W，n+h=H}，每个像素块用来表示其中心位置像素点的变化特征，同样按行列展开为h2×1 的向量，记作(H-h)，按式 (6) 和 (7) 将映射到EVS中去：

最终通过 PCA 算法，将差分图像大部分区域像素的变化特征转化为维度间线性无关的N维向量，称为图像变化特征矩阵。

2.2.3 像素点二分类

使用聚类算法依据上一步得到的图像变化特征矩阵对差分图像像素点进行二分类，相较于传统的k-Means 聚类（k= 2）[12]，做出如下改进：

1）比较使用多种k值，由于不同时期的卫星图片在整体上区别不大，而改变总在细微处，所以可以认为包含像素点最少的一类，是发生改变的一类。

2）为选择初始聚类中心，优化k-Means 聚类方法，获得更快的迭代速度和更好的聚类效果，采用如下算法：a. 输入数据集X，聚类中心个数k，在数据集中随机选取 1 个样本作为初始聚类中心C1。b. 计算每个样本到与之最近的 1 个聚类中心的距离D(x)，x∈X；计算每个样本被选为下一聚类中心的概率选择概率最大的样本作为下一聚类中心。c. 重复步骤 b 直到选择出k个聚类中心。

设定“1”表示某像素点发生改变，“0”表示该像素点没有发生改变，因此分类结果图中，黑色区域表示没有变化，白色区域表示发生变化。

3 实验与分析

3.1 数据来源与预处理

用一组拍摄时间分别为 2015 年 6 月和 2018 年6 月的江苏省常州市某区域遥感影像数据进行实验，该区域遥感影像图片如图 2 所示。

图 2 常州市某区域遥感影像图片

图片大小为 8 000 像素×7 000 像素，图片格式为 TIFF；影像信息格式为 TFW。图 2 所包含影像信息内容及含义如表 1 所示，其中：像素分辨率表示单位像素宽度代表的实际距离，单位为 m/像素，在 0.3 m/像素的分辨率尺度下，每张图片的实际地理长宽分别是H= 0.3 ×8 000 = 2 400 m，W= 0.3 ×7 000 = 2 100 m。

数据预处理主要包括：

表 1 某遥感影像图片信息

1）成对图片位置校正。由于获取的遥感影像可能存在图片大小误差、地理位置偏差等异常情况，误差示意图如图 3 所示。

图 3 误差示意图

若A，B点的像素坐标分别为 (x1，x2)，(y1，y2)，地理坐标分别为 (x1'，x2')，(y1'，y2')，2 个图像素分辨率都是α，那么A，B像素坐标及地理信息满足下式：

分别对 2 个时期遥感影像图片剪裁出图 3 中黄色区域，代替原始数据，完成图片位置校正。

2）切割图像。由于原始图像过大，容易造成内存溢出，所以将原始图像切割成长宽都为 1 000 像素的矩形。

3）转化成灰度图。将红、绿、蓝三通道图像转换为单通道灰度图。

预处理完成后，每次输入算法的都是 2 张1 000 像素×1 000 像素的单通道同区域不同时期遥感影像图片。

3.2 对比实验

在本研究提出的方法中，主要超参数包括像素块大小h和聚类中心个数k，选用如表 2 所示的 6 组超参数进行实验。

在实验中，前 3 个主成分的累积贡献率达到了99.7%，故取主成分个数N= 3，对比实验结果如图4 所示。

对比方案 1，2，3 或 4，5，6，发现像素块大小h对实验结果的影响不显著；对比方案 1，2 或3，4 或 5，6，发现聚类中心个数k决定了变化区域的阈值，k值越小，结果图像的椒盐化现象就越显著，误分类率也随之提高。图 5 所示是图像上部一块区域的实际图片与实验结果，从实际图像中可以看出，该区域变化主要是下方农田改建为道路的施工用地，在k= 3 个时，算法将上方一块池塘误归类为变化区域，这是由于前后 2 个时期遥感影像光线不同导致的图像颜色深浅不一致；而k= 5 个的方案很好地解决了这一误分类。

表 2 实验超参数选择

图 4 各参数方案提取的变化图像

图 5 某区域原始图片与实验结果

在定性分析的基础上，保持像素块大小h= 5 像素恒定，选用 2～6 这 5 种不同k值，使用样点检验的方式进行定量分析。依据特定原则选取 500 个样本点比较各算法的精度[4]。人工判断这 300 个样本点的集合归属，其中变化样本 166 个，未变化样本334 个，计算各个方法的召回率R、准确率P和总体精度结合传统 CVA 方法比较结果如表 3所示。

表 3 算法精度比较

从表 3 可以看出：随着k值增加，准确率提高至收敛，召回率却不断降低；同时，本研究算法普遍优于传统的 CVA 算法，在k= 4 个取得最优的总体精度，通过上述分析，实验的最优参数选择为h=5 像素，k= 4 个。

3.3 结果分析

从图 5 中可以看出：常州市某区 2015—2018 年的主要变化体现在区域上部新修一条主干道，其他的一些典型变化列举如图 6 所示。

上述用地变化是否涉及违章建筑需要人工进一步比对审查，但在总体上，遥感影像快速变化检测算法提供了可能发生违建的大致位置，节省了人力的比对。

4 结论

本研究为解决违章建筑的粗筛选问题，以常州市某区域遥感影像为例，改进了基于 PCA 和k-Means 聚类的变化检测算法，通过选择不同的分割像素块边长和聚类中心个数的超参数组合进行比对实验。结果表明：在精度水平上该算法优于传统CVA 算法，并且在像素块边长为 5 像素、聚类中心个数为 4 个时，既能有效抑制噪声，减小伪变化的出现，又能保证一定的召回率，使得整体精度达到理想效果，并在测试用例中筛选出池塘新修小道、池塘变为农田等变化区域。

图 6 典型变化示例