研究区位于河北省张家口市怀来县东花园镇, 东西宽约200 m, 南北长约250 m。 根据研究区大小设计无人机飞行轨迹如图1所示。 该区地处中温带半干旱地区, 四季分明, 光照充足, 雨热同期, 属于温带大陆性季风气候, 年均日照时数3 027 h, 平均气温9.1 ℃, 年均降水396 mm, 昼夜温差大。 实地考察表明该区域包含城镇不透水面的基本类型, 具有典型性。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 研究区地理位置与航迹规划Fig.1 Geographic location of the study area and route planning

采用M600无人机搭载pikaL高光谱传感器采集研究区影像数据, 光谱分辨率为2.7 nm, 空间分辨率为0.1 m, 像元尺寸为5.86 μ m。 为减少阴影对分类结果的影响, 选择光线良好的时段进行航拍。 对影像数据进行波段选择, 剔除受水汽和二氧化碳影响严重、 噪声多、 质量差的波段, 最终保留150个有效波段(394~1 032 nm)。 之后进行大气校正和影像拼接, 并去除畸变严重区域, 得到影像如图2所示。 结合影像进行目视解译, 发现不透水表面主要包括三类材质: 混凝土(包括道路、 停车场、 高层建筑)、 金属面(蓝色屋顶)、 瓦片(红色屋顶), 透水面主要包括两类裸土、 植被(作物、 灌木和树木)。 相关地物类型如表1。 典型不透水地面和透水面的光谱特征如下: 混凝土表面反射率较高, 光谱曲线在可见光波段较为平坦, 但在近红外波段有明显反射率变化; 金属屋顶在可见光波段有特定的吸收特征, 反射率较低, 但在近红外波段反射率较高; 瓦片屋顶在红光波段有较高的反射率, 而在其他可见光波段反射率较低; 植被在近红外波段反射率显著高于可见光波段, 具有典型的“ 红边” 特征; 土壤光谱特征较为复杂, 但通常在可见光和近红外波段有较高的反射率, 与植被有明显区别, 这些光谱特征的差异为后续的分类提供了依据。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 研究区高光谱影像处理结果Fig.2 Preprocessed hyperspectral image of the study area

表1

Table 1

表1(Table 1)

表1 城镇典型地物类型 Table 1 Typical land cover types of towns

表1 城镇典型地物类型 Table 1 Typical land cover types of towns

为了降低模型处理负担并消除因波段间高度相关性引起的数据冗余, 需要对无人机高光谱数据中的150个波段进行筛选, 选出最适合提取城镇不透水面的特征波段。 波段选择应符合以下要求: (1) 所选波段具有较大的信息量; (2) 所选波段之间的相关性应尽可能低; (3) 波段对不同地物类别有明显的光谱差异。 基于这些原则, 采用逐步判别分析法进行波段选择, 同时使用波段标准差分析、 波段间相关性和主成分分析综合分析判断。

2.1.1 逐步判别分析法

逐步判别分析是一种动态调节变量的方法, 用于选择高光谱影像中最具判别能力的特征波段。 该方法初始状态下不包含任何变量, 逐步引入具有最强判别能力的变量, 同时剔除那些因引入新变量而判别能力下降的旧变量^[8]。 通过对判别式的统计检验, 当没有符合标准的变量可引入时, 判别过程终止, 最终的判别函数仅包含少量但判别能力强的变量。 逐步判别分析使用Wilk's Lambda来度量引入或剔除变量的标准, 用Λ 表示。 Wilk's Lambda的取值范围在0到1之间, 数值越小代表变量对判别模型的贡献越大。 因此, 在每一步中, 选择Λ 值最小的变量加入模型。 逐步判别分析的整体过程如图3所示。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 逐步判别分析流程图Fig.3 Flow chart of the step-by-step discriminant analysis

在使用逐步判别分析法选择波段前, 首先需要确定不同地物类别在所有波段中的反射率。 选定了城镇遥感影像中的五种常见地物类型, 包括不透水面的混凝土、 金属材质屋顶、 瓦片屋顶, 以及透水面的植被和裸土。 根据实地调研信息和地理坐标, 在遥感影像中为五类地物各选择50个样本区域, 记录这些样本点在高光谱影像150个波段中的反射率值, 作为标准分类组。 然后, 使用SPSS软件中的逐步判别分析模型, 对150个波段分别进行分析, 最终选择了14个满足条件的特征波段, 分别为394、 407、 443、 509、 517、 533、 592、 625、 680、 689、 710、 749、 848和911 nm。

2.1.2 波段信息量分析

使用波段标准差来评价其信息量, 标准差越大, 包含的地物信息量越丰富。 图4显示了东花园镇图像各波段的标准偏差变化曲线。 总体来看, 从可见光到近红外波长范围内, 各波段的标准差随波段增加先上升后下降, 在可见光绿色波段达到最大。 对150个波段的标准差进行排序, 前十的波段分别为554、 538、 550、 546、 533、 562、 558、 529、 542和592 nm, 表明可见光绿色波段信息量大。 同时在近红外波段848 nm和可见红光波段749nm有明显的峰值, 这些波段皆为该类波段中信息量大的波长。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 试验区高光谱数据波段标准差变化Fig.4 Standard deviation variation of the hyperspectral data in the test area

2.1.3 相关性分析

波段相关性通常通过计算相关系数进行评估。 图5显示了高光谱图像150个波段之间的相关系数。 总体来看, 150个波段间的相关性较高, 但有一定的块状特征。 A区块为B1到B30波段区间(394~513 nm), B区块为B31到B90波段区间(517~766 nm), C区块为B91到B130区间(770~941 nm), D区块为B131到B150区间(946~1 032 nm)。 这些区块内部波段的相关系数非常高, 超过0.95, 但不同区块之间的相关系数较低, 如AD区块和AC区块的平均相关系数分别为0.78和0.82。 为确保所选波段能表达独立且全面的特征信息, 最佳特征波段应从不同区块中挑选。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 各波段间相关系数矩阵图Fig.5 Matrix diagram of the correlation coefficient between the various bands

2.1.4 主成分分析

对高光谱影像的150个波段进行主成分分析(PCA), 以解决波段过多和信息重叠的问题。 图6展示了各个主成分的特征值分布。 结果表明, 前5个PCA分量占原始影像光谱信息的99.17%, 其中前两个主成分PC1和PC2占比之和为97.06%。 这表明经过PCA运算后, 150个波段的信息几乎全部集中在前5个PCA分量中, 且以PC1和PC2为主导。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 影像PCA各分量所对应的特征值Fig.6 The characteristic values corresponding to each component of the image PCA

表2统计了各波段对PC1— PC5的贡献。 表3展示了每个PCA成分的贡献率排序的前10个波段。 从表中可知, 对PC1贡献最大的波段主要集中在近红外波段。 对PC2贡献最大波段主要集中在红光波段。 由于PC1和PC2占原始光谱信息的97.06%, 并包括大部分原始图像信息, 因此这些主导波长在区分不透水表面和透水表面时具有重要意义。

表2

Table 2

表2(Table 2)

表2 影像PCA分量特征值所占比例 Table 2 Proportion of image PCA components 主成分分量 特征值所占比例/% PC1 92.20 PC2 4.86 PC3 1.61 PC4 0.37 PC5 0.13 总计 99.17

表2 影像PCA分量特征值所占比例 Table 2 Proportion of image PCA components

表3

Table 3

表3(Table 3)

表3 与前五个主成分相关的波段 Table 3 The spectral regions well explained by each principal component PCA波段 波长/nm 主导波段 PC1 394、 937、 398、 933、 946、 951、 910、 955、 407、 960 近红外 PC2 753、 748、 745、 775、 779、 771、 784、 758、 788、 740 红光 PC3 480、 476、 472、 484、 467、 463、 459、 455、 492、 488 蓝光 PC4 680、 459、 685、 455、 676、 463、 451、 672、 467、 668 蓝光、 红光 PC5 749、 753、 745、 758、 740、 736、 771、 775、 779、 732 红光

表3 与前五个主成分相关的波段 Table 3 The spectral regions well explained by each principal component

2.1.5 特征波段选择结果分析

为了对逐步判别分析得到14个特征波段进行更全面的考量, 对比了不同方法的结果。 分析显示, 这些波段涵盖了标准差排名靠前, 信息量丰富的绿色波段(如533和592 nm)、 近红外波段(848 nm)及红光波段(749 nm)。 根据相关性分析结果表明, 属于不同区间的波段的相关性更小。 选择的14个波段中394、 407、 443和509 nm在A区, 517、 533、 592、 625、 680、 689、 710和749 nm属于B区, 848和911 nm属于D区。 这三个区块两两组合相关性较小。 主成分分析揭示, 主导波段为近红外与红光波段的PC1和PC2占原始光谱信息的97.06%, 且所选波段中半数为这两种类型, 其中入选波段394与407 nm在PC1中贡献显著。 图7展示了各类地物在150个波段下的光谱特征, 所选波段为灰色虚线所在位置。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 各类地物的光谱特征曲线Fig.7 Spectral characterisation curves for various land covers features

卷积神经网络近年来在众多研究领域得到应用, 其在图像分类、 目标检测等诸多任务上都表现优异。 本研究主要使用改进的AlexNet网络架构进行城镇不透水面的提取工作。 AlexNet网络结构相对简单直观, 易于理解, 相比目前一些深度和宽度更大的模型, 其参数量较少, 训练和推理速度会更快, 在面对中小规模的数据集或资源受限的环境时, 具有一定的优势。 AlexNet原网络共有8个层次, 其中包括5个卷积层和3个全连接层, 此外网络还有三个最大池化层。 由于池化层的降维作用会使训练样本在高层次的尺寸过小^[28], 引起信息丢失和限制模型的表征能力。 为了提高对阴影区域的分类精度, 本工作设计的卷积网络仅保留最后一层池化层, 同时设计4个卷积层, 参数设置分别为32@7× 7、 64@7× 7、 64@5× 5和1 28@5× 5, 卷积核大小分别为1、 1、 3和1, 步长统一设置为1。 在池化层部分选择使用2× 2大小、 步长为2的最大池化方法, 最后设计两个节点数分别为256和128的全连接层, 整体模型结构如图8所示。 需要说明的是, 由于实验数据具有14个波段, 因此需要使用GDAL将遥感影像转换为Numpy格式, 并将其转换为深度学习框架的张量, 然后输入模型进行处理。 整个网络的学习率初始化为0.001, 权重衰减因子为0.001, 全连接层前后使用的Dropout分别设置为0.25和0.5, 采用了Adam优化算法。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 卷积神经网络模型结构图Fig.8 Structure diagram of the convolutional neural network model

本研究选取了混凝土、 金属材质、 瓦片、 植被和裸土五类地物, 各类选取50个样本区域, 共计16, 160个像素级训练样本。 由于样本数量差异明显, 使用重采样来平衡各类别的样本数量, 如表4所示。 每类地物随机选取70%的样本用于训练, 剩余30%的样本用于验证。

表4

Table 4

表4(Table 4)

表4 各类型的像素级样本数量 Table 4 Number of pixel-level samples of each type 类型 混凝 土地面 金属 屋顶 瓦片 屋顶 植被 裸土 原始样本数 28 191 18 263 25 386 15 241 15 105 重采样后样本数 15 105 18 463 15 105 15 241 15 105

表4 各类型的像素级样本数量 Table 4 Number of pixel-level samples of each type

为了验证本方法的有效性, 设计两步验证思路: (1) 对比使用高光谱原始影像和特征波段影像的实验结果, 评估高光谱特征波段选择方法对城镇不透水面提取的效果。 (2) 将提出的卷积神经网络模型的提取结果与AlexNet网络及三种常用机器学习方法: 朴素贝叶斯(NB)、 支持向量机(SVM)和随机森林(RF)的结果进行对比, 证明所设计网络模型的优越性。 论文精度评估指标包括: 总体精度(OA), 用户精度(UA), 生产者精度(PA), Kappa系数和Macro F1(F1)。

OA=TPTP+FP+FN×100%(1)

Kappa=p0-pe1-pe(2)

F1=2×UA×PAUA+PA(3)

其中, 真阳性(TP)表示正确分类的特定类型像素数量, 假阳性(FP)表示错误分类为特定类型的像素数量, 假阴性(FN)表示特定类型像素被错误分类的数量。 P₀表示观察到的一致比例, P_e表示偶然期望的一致比例。