拉萨市(91° 06'E, 29° 36'N)位于青藏高原中南部, 地势整体由东向西倾斜, 地形以高原山地为主, 平均海拔约3 650 m。 行政区域内现辖5县3区, 东西长约277 km, 南北宽约202 km。 该地区属高原温带半干旱季风气候, 农业主要种植青稞、 小麦等喜凉作物, 植被覆盖垂直地带性显著。 受喜马拉雅山脉对当地大气环流的影响主下沉气流, 全年多晴朗天气, 被称为“ 日光城” , 作为重要生态区域青藏高原的核心发展城市, 人口聚集程度高, 社会基础设施完善, 在政教文卫各方面起着带头作用^{[14, 15]}。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 拉萨市DEM及行政区域分布Fig.1 Digital elevation model (DEM) and administrative division of Lhasa City

Sentinel-2(S2)卫星搭载多光谱成像仪(MSI), 覆盖4个10 m、 6个20 m及3个60 m波段, 两颗卫星互补重访周期为5 d。 GEE云平台提供的地表反射率产品Sentinel-2 MSI已经过正射校正和大气校正, 可节约数据获取及处理时间成本, 保证数据质量。

拉萨市地区夏季植被生长旺盛, 闲荒的农田、 裸地等易与人造地表混淆的覆盖类型较少^[16]。 选取的拉萨市数据的时间范围为2022年7月1日至8月31日, 使用QA60波段对其进行去云处理, 并将质量较差的像元掩膜, 根据需求选择Green, Red, Blue, NIR, SWIR1, SWIR2共六个波段。

基于S2结合Google Earth Pro高分辨率遥感影像进行目视解译, 采用随机抽样的方法, 在拉萨市的遥感影像上进行训练数据的选取。 通常将目视解译得到的全部样本点通过分层随机抽样的方式, 按7：3的比例分为训练及验证数据^{[5, 11]}, 但这种方式易使模型过拟合, 从而产生过高的验证精度, 影响研究结果的判断。 因此, 本研究的测试样本点在研究区内随机生成, 再利用遥感影像目视解译对其进行类型判别赋予属性, 同时考虑到各地表覆盖类型样点数据的比例, 在人造地表聚集区重点生成。 本研究共生成了2 028个样本点, 其中训练样本1 420个, 验证样本608个, 样本点具体分布状况见图2。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 采样点分布状况Fig.2 Distribution of sampling points

大量研究结果表明, 分类特征构建会产生巨大影响, 使用全面的特征组合是克服遥感影像“ 异物同谱、 同谱异物” 现象的重要手段^[17]。 拉萨市居民点主要分布在地形平坦、 水热条件优越的河谷地区及交通主干道沿线, 可知该地城镇化形态以地形为主要驱动因素。 地形和纹理特征可以从不同方面描述数据的表面形态信息, 添加如高程、 坡度、 坡向的形态特征作为辅助数据能够提高分类精度^[11]。 在GEE云平台调用研究区30 m分辨率的SRTM数据, 重采样至10 m后通过该DEM数据计算提取高程、 坡度、 坡向因子。

空间纹理特征能够在中分辨率影像的分类上提供光谱信号所不具备的信息, 通常被认为是提取低密度人造地表覆盖的重要辅助特征^[13]。 在基于像元及面向对象的分类方法中对比加入和不加入纹理特征的结果精度, 分析纹理特征在10 m分辨率数据源下进行城市尺度人造地表提取中的作用力。 使用基于统计的灰度共生矩阵(gray-level co-occurrence matrix, GLCM)算法提取纹理特征, 基于S2影像信息通过式(1)加权线性组合得到灰度图像。

Gray=(0.3×NIR)+(0.59×Red)+(0.11×Green)(1)

式(1)中, NIR为近红外波段反射率, Red为红波段反射率, Green为绿波段反射率。

对计算得到的纹理信息采用主成分分析(principal component analysis, PCA)压缩, 其中第一主成分(principal component 1, PC1)被认为包含绝大多数的结构信息^{[11, 18]}, 因此使用PC1作为本研究的纹理特征。

夜间灯光指数是拉萨市人口分布的主导因子^[19], 人类聚集性活动不断推动人造地表覆盖的发展, 因此认为夜间灯光(nighttime light, NTL)遥感数据可直接反映城镇化的分布概况, 是人造地表覆盖提取的重要辅助数据^{[20, 21]}。 EVI调整夜间灯光指数(EVI adjusted nighttime light index, EANTLI)是EVI参与调整后的NTL指数, 该指数能够有效克服裸土效应, 突显城市内部差异^[22], 对DMSP-OLS和NPP-VIIRS两种夜间灯光数据集的处理都产生了很好的校正效果^[20]。 VIIRS数据时空分辨率均高于DMSP数据, 且很大程度上降低了DMSP数据中的过饱和问题^[21]。 因此, 通过GEE提供的2022年EVI数据调整NPP-VIIRS数据, 重采样至10m计算EANTLI指数参与分类。

多种光谱指数根据不同的研究需求不断发展, 在增加可用信息的基础上改善最终分类结果^[11]。 本研究以光谱特征为主, 结合地形、 纹理共23个特征(见表1), 为探明在特征集合中加入纹理特征对人造地表提取的精度影响, 另一组特征剔除纹理特征。

除地形因子、 夜间灯光指数及纹理特征外, 还使用了在城镇用地信息相关提取中被广泛应用的12种光谱指数。 使用夏季植被生长旺盛时期数据提取NDVI、 EVI有助于提升植被区分度; SAVI、 MSAVI、 BSI、 IBI、 EIBI、 NDBI、 NDBBI可减弱大气及地形的影响, 将易被混淆为人造地表的土地类型如土壤及裸地等与建筑区差异增强; NDWI、 MNDWI能够提升水体识别能力, 避免大片建筑物阴影误判, 这些特征组合可增强对研究区内大片水体、 植被、 裸地、 农田等覆盖类型的区分能力, 优化人造地表覆盖提取结果^{[23, 24]}。 将上述使用的所有相关数据重采样至10 m分辨率参与分类提取过程。

表1

Table 1

表1(Table 1)

表1 人造地表提取特征选择 Table 1 Features selection for artificial surface extraction 类型 特征 描述 波段反射率 Blue(B2) 蓝波段反射率 Green(B3) 绿波段反射率 Red(B4) 红波段反射率 NIR(B8) 近红外波段反射率 SWIR1(B11) 短波红外1波段反射率 SWIR2(B12) 短波红外2波段反射率 地形因子 高程(Elevation) 坡度(Slope) 坡向(Aspect) 纹理特征 第一主成分(PC1) 光谱指数 归一化植被指数 (normalized difference vegetation index, NDVI) NIR-RedNIR+Red 增强型植被指数 (enhanced vegetation index, EVI) 2.5×(NIR-Red)NIR+6×Red-7.5×Blue+1 归一化差值建筑用地指数 (normalized difference built-up index, NDBI) SWIR1-NIRSWIR1+NIR 归一化差值裸地与建筑用地指数 (normalized difference bareness and built-up index, NDBBI) [1.5×SWIR2-(NIR+Green)/2][1.5×SWIR2+(NIR+Green)/2] 归一化差值水体指数 (normalized difference water index, NDWI) Green-NIRGreen+NIR 改良归一化水体指数 (modified normalized difference water index, MNDWI) Green-SWIR1Green+SWIR1 土壤调整指数 (soil adjusted vegetation index, SAVI) (NIR-Red)(1+0.5)NIR+Red+0.5 改良土壤调整指数 (modified soil adjusted vegetation index, MSAVI) (2×NIR+1)-[(2×NIR+1)2-8×(NIR-Red)]2 裸地指数 (bare soil index, BSI) [(SWIR1+Red)-(NIR+Blue)][(SWIR1+Red)+(NIR+Blue)] 建筑用地指数 (index-based built-up index, IBI) IBIA= 2×SWIR1NIR+SWIR1 IBIB= NIRNIR+Red+ GreenGreen+SWIR1 IBI= IBIA-IBIBIBIA+IBIB 增强指数型建筑用地指数 (enhanced index-based built-up index, EIBI) NDBBI-4×BSI-MNDWIBSI+MNDWI+SAVI+MNDWI/6NDBBI+4×BSI-MNDWIBSI+MNDWI+SAVI+MNDWI/6 增强型裸土指数 (enhanced bare soil index, EBSI) BSI-MNDWI)(BSI+MNDWI) EVI调整夜间灯光指数 (EVI adjusted nighttime light index, EANTLI) (1-EVIm+NTLn)(1+EVIm-NTLn)× NTL Note: EVIm represants the Annual Mean MODIS EVI and NTLn represents normalizaed nighttime light value 注: EVIm为MODIS年度植被指数EVI的平均值, NTLn为夜间灯光归一化值

表1 人造地表提取特征选择 Table 1 Features selection for artificial surface extraction

1.5.1 研究方法

1.5.1.1 研究方法技术流程

使用2022年7月— 8月S2遥感数据, 基于光谱特征对拉萨市人造地表进行提取。 对比输入特征组合是否包含纹理特征, 使用基于像元还是面向对象分类方法, 基于SVM或RF分类器时的分类结果, 以期获取10 m分辨率在生态城市尺度进行人造地表覆盖提取时的最优方法组合。 共形成了八种分类方法组合, 包括基于像元方法中(1)SVM分类器(PB_SVM), (2)输入纹理特征的SVM分类器(PB_GLCM_SVM), (3)RF分类器(PB_RF), (4)输入纹理特征的RF分类器(PB_GLCM_RF); 面向对象方法中通过简单非迭代聚类(simple non-iterative clustering, SNIC)算法分割影像后(5)SVM分类器(OO_SVM), (6)输入纹理特征的RF分类器(OO_GLCM_SVM), (7)RF分类器(OO_RF), (8)输入纹理特征的RF分类器(OO_GLCM_RF)。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 研究方法流程图Fig.3 Research methodology flowchart

在过去机器学习的研究过程中, 地表覆盖提取方法被发展为基于像元及面向对象两大部分, 两者区别主要表现为处理数据的基本单元不同。

1.5.1.2 基于像元的分类方法

SVM分类算法在上世纪九十年代被开发^[25], 作为一个二值分类器, 该方法在高维空间中利用输入向量寻找最优超平面, 将两类训练向量分离到两个子空间中以获得分类结果, 遵循结构风险最小化原则, 能够很好的平衡风险与分类器性能。 SVM以其良好的表现在城镇用地相关研究中被广泛使用, GlobeLand30、 MSMT_IS30-2015等数据集产品的处理过程都使用了该方法^[17]。

RF分类算法是由Breiman^[26]2001年在决策树的基础上发展而来的集成分类器, 通过控制内部决策树数量, 能够有效避免模型出现过拟合并提高泛化性及预测能力^[27]。 在实际应用过程中发现, 当决策树数量达到一定程度后, 棵数增加对结果精度只会带来微弱波动, 但会增加数据量造成计算存储的负担。 RF以其灵活性、 高效率、 准确性等各方面的优势, 被应用于多源遥感及地理数据的处理, 在遥感、 生态等领域得到了广泛应用^[28]。

1.5.1.3 面向对象的分类方法

以对象为单位的分类方法能够优化基于像元结果中的“ 椒盐” 现象, 在地表覆盖提取、 生态遥感、 动态监测等方面都得到了广泛应用。 影像分割尺度参数是面向对象方法中最重要的部分, 分割原则是保证影像分割对象内部异质性低, 同时不同对象之间的差异性显著, 以提升不同覆盖类型之间的区分度^[29]。

为了更好的发掘面向对象方法的优势, 获得更有意义的分割对象, 多种图像分割算法被开发, 主要包括基于图、 梯度上升两类方法, 如NC05、 GS04、 SL08、 WS91等^[30]。 与这些传统方法相比, Achanta等^[30]提出基于K-means聚类生成的超像素方法(SLIC), 该方法在处理速度、 分割质量、 参数调优等方面具有绝对性的优势, 更适用于海量数据处理。 本研究使用的SNIC算法是SLIC的进阶版本, 兼具性能稳定与高效率内存计算的优势。

尺度分割的参数对分类结果精度影响极大。 综合考虑研究区地表覆盖分布特性、 计算内存压力及分类结果精度, 衡量数据计算存储压力与分割结果, 设置SNIC分割方法参数为: 种子距离(Size): 10, 分割紧密度(Compactness): 0, 连通性(Connectivity): 8, 邻域大小(Neighborhoodsize): 256^[31]。

1.5.2 精度评定

使用混淆矩阵, 从不同角度进行分类结果的精度评定, 包括: 总体精度(overall accuracy, OA)、 Kappa系数(Kappa coefficient, KC)、 生产者精度(producer accuracy, PA)、 用户精度(user accuracy, UA)等指标。 其中OA代表总体数据的精度, 可直接判断分类效果的优劣。 而Kappa系数兼顾正确、 错分及漏分数据, 更能够体现分类与实地类型的适配程度, 对分类结果精度更具指示性意义, 高于0.8则认为分类结果基本正确。 PA代表分类结果中该类别得到的正确数量与该类别在测试集中实际具有的数量的比值, UA代表在分类结果中该类别正确分类的数量与影像中被分类判定为该类别的总数的比值。

这些精度评定指标的计算公式如式(2)— 式(5)

OA=∑i=1naii∑i=1n∑j=1naij(2)

Kappa=N∑i=1naii-∑i=1n(ai+·a+i)N2-∑i=1n(ai+·a+i)(3)

PA=aij∑a+j(4)

UA=aij∑ai+(5)

式(2)— 式(5)中, ∑i=1naii代表正确分类样本数量和, ∑i=1nai+和 ∑i=1na+i分别代表第i行、第i列的样本数量和。

不同分割尺度会产生不同的结果, 影响到最终提取结果的精度。 由表2不同种子距离尺度分割影像提取精度可知, 总体精度及Kappa系数随种子像素距离增大而逐渐提升, 在种子像素距离为10像素时达到峰值, 随后逐渐减小至最小值附近波动。 因此, 选取研究区的最优种子距离为10像素。

表2

Table 2

表2(Table 2)

表2 不同种子距离尺度分割影像提取精度对比 Table 2 The comparison of extraction accuracy for image segmentation at different seed distance scale 种子距离 /像素 OO_SVM OO_RF(150) OA KC OA KC 2 0.940 8 0.855 1 0.973 7 0.936 0 4 0.940 8 0.854 2 0.972 0 0.931 9 6 0.942 4 0.858 0 0.977 0 0.943 9 8 0.947 4 0.870 4 0.973 7 0.936 0 10 0.955 6 0.891 9 0.980 3 0.952 0 20 0.942 4 0.858 0 0.965 5 0.915 9 30 0.935 9 0.842 8 0.955 6 0.891 5 40 0.939 2 0.850 4 0.962 2 0.907 9 50 0.929 3 0.825 5 0.945 7 0.866 6 60 0.926 0 0.818 6 0.940 8 0.854 7 70 0.916 1 0.788 9 0.945 7 0.865 7 80 0.916 1 0.790 3 0.940 8 0.854 7

表2 不同种子距离尺度分割影像提取精度对比 Table 2 The comparison of extraction accuracy for image segmentation at different seed distance scale

共生成8组提取组合方法, 各分类方法获得的OA、 Kappa系数、 UA、 PA见表3。 对比结果可发现, 使用不同分类器时面向对象方法结果精度基本都高于基于像元方法, 不使用纹理特征参与RF分类时, 面向对象OA、 Kappa分别提升了0.005 7、 0.012 2; 不使用纹理特征参与SVM分类时, 面向对象OA、 Kappa分别提升了0.019 9、 0.049 1; 使用纹理特征参与SVM分类的方法, 面向对象OA、 Kappa分别提升了0.011 4、 0.025 2。 仅在使用纹理特征参与RF分类的方法时有差异, 面向对象OA、 Kappa分别降低了0.004 4、 0.011 9。

表3

Table 3

表3(Table 3)

表3 拉萨市人造地表覆盖提取不同方法结果精度 Table 3 Accuracies of different methods for extracting artificial surface coverage in Lhasa City 方法 种子距离 /像素 无纹理特征 有纹理特征 OA KC UA PA OA KC UA PA PB_SVM 0.935 7 0.842 8 0.833 3 0.949 4 0.942 5 0.862 5 0.883 3 0.924 4 PB_RF(150) 0.974 6 0.939 8 0.950 0 0.966 1 0.979 7 0.951 8 0.955 6 0.977 3 OO_SVM 10 0.955 6 0.891 9 0.900 0 0.947 4 0.953 9 0.887 7 0.894 4 0.947 1 OO_RF(150) 10 0.980 3 0.952 0 0.944 4 0.988 4 0.975 3 0.939 9 0.933 3 0.982 5

表3 拉萨市人造地表覆盖提取不同方法结果精度 Table 3 Accuracies of different methods for extracting artificial surface coverage in Lhasa City

分类组合方法中的最优方法及参数设置为: 在面向对象方法中使用SNIC算法对影像进行分割, 设置分割参数为Size: 10, Compactness: 0, Connectivity: 8, Neighborhoodsize: 256, 使用150棵树的RF分类器进行分类, 得到总体精度为98.03%, Kappa系数为0.952 0, 人造地表覆盖提取准确。

各组合方法通过参数设置获得的最优人造地表覆盖提取结果图见图4及图5, 其中图4为基于像元方法, 图5为面向对象方法。 对比图4(a)、 图4(b)与图4(c)、 图4(d)及图5(e)、 图5(f)与图5(g)、 图5(h)可发现, 对相同尺度的遥感数据, SVM比RF分类器提取结果范围更广, 结果高估于RF算法。 对于城市内部流域纵穿区域, RF对其分隔能力更强, 而SVM将更多的流域内干涸硬土地识别为人造地表覆盖类型, 识别能力较弱。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 基于像元方法拉萨市人造地表覆盖提取结果图 (a): PB_SVM; (b): PB_GLCM_SVM; (c): PB_RF; (d): PB_GLCM_RFFig.4 Artificial surface extraction in Lhasa City based on pixel-based method (a): PB_SVM; (b): PB_GLCM_SVM; (c): PB_RF; (d): PB_GLCM_RF

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 面向对象方法拉萨市人造地表覆盖提取结果图 (e): OO_SVM; (f): OO_GLCM_SVM; (g): OO_RF; (h): OO_GLCM_RFFig.5 Artificial surface extraction in Lhasa City based on object-oriented method (e): OO_SVM; (f): OO_GLCM_SVM; (g): OO_RF; (h): OO_GLCM_RF

为了分析在输入分类器的特征中是否加入纹理信息, 共产生了包括基于像元/面向对象、 使用SVM/RF分类器的四组结果。 由结果可知, 纹理特征的作用效果不大, 甚至可能会对分类结果的精度产生负面影响。 分析其原因, 地形及纹理特征都是用于描述遥感数据表面特征信息, 两者之间相关性较高, 相互作用使得数据冗余, 在一定程度上减弱了纹理信息的作用能力。 使用更高分辨率的数据, 能获取更为精细的信息, 纹理特征可能会有更好的表现。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 纹理特征输入与否分类结果对比Fig.6 Comparison of accuracy between method with and without texture features input

对比使用不同分类器的结果精度可知, 不管使用基于像元还是面向对象方法, 是否输入纹理特征参与分类, 使用RF分类器都比SVM分类器能够获得更好的结果。 随着数据时空分辨率的发展, 大多使用长时序、 多源数据结合进行遥感数据分析处理。 RF分类器在处理高维多源数据时具有稳定性, 能够在较短的计算时间内获得较高的精度, 具有高处理效率。 通过目视解译手段获取训练数据仍被广泛使用, 这种方式获得的数据有概率存在主观误差, 使用具有鲁棒性的RF分类器, 可抵抗该误差对结果的影响^[10]。 因此, 在人造地表相关信息提取的研究中, 选择不同机器学习算法获得结果的差异明显, 使用RF能够获得更好的结果。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 SVM/RF分类结果精度对比Fig.7 Comparison of accuracy between SVM and RF classification results

保持基于像元和面向对象方法使用的数据源、 输入特征及分类器相同, 对比了两种分类方法提取人造地表的效果。 仅基于整体分类精度判别时, 两种分类结果精度不存在较大的差异, 但面向对象比基于像元方法占用更多运算时间^[32]。 对使用相同分类器的分类结果进行目视解译判别, 从结果图细节对比可发现, 面向对象方法比基于像元的结果图噪点更少, 图像效果更好。 使用基于像元方法对人造地表提取时, 考虑到城市内部异质性较强, 难以保证单一像元内部均为纯净, 在一定程度上导致提取结果的误差; 而基于对象的方法以分割后的均质对象为单位, 充分考虑不同对象之间的邻域拓扑关系, 能够减少分类结果中的“ 椒盐” 现象^[30], 因此使用面向对象方法获得的图像的细节更为准确。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 结果细节噪点情况对比图Fig.8 Comparison of noise details

不透水面由人类活动建造, 是城镇化的主要的地表覆盖类型, 对人造地表的提取结果具有重要参考意义。 使用认可度高、 更新及时的高分辨率全球不透水地表相关数据集作为提取结果的参照数据。 各数据集在拉萨市行政区域内的分布状况见图9, 包括30 m全球土地覆盖数据集(GlobeLand30, 2020)^[11]; 全球精细分辨率土地覆盖观测与监测(FROM_GLC, 2017)^[33]; 全球30米不透水表面动态数据集(GISD30, 2020)^[34]。

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 高分辨率全球不透水地表相关数据集合对比Fig.9 Comparison of high-resolution global impervious surface datasets

通过图像对比可发现, 结果与各数据集城市范围分布一致性强, 仅在交通轨道沿线及城郊、 村镇等人口密集程度低的区域差异较大, 证明本研究结果具有可靠性, 是进行城市乃至更大尺度的人造地表覆盖提取的有效手段。 而离散居民点处的区别较大的原因, 主要在于原始数据影像的空间分辨率, 及不可避免存在混合像元使得地物识别能力受限; 加上各研究使用的采样点数量、 质量、 分布等状况不一致, 导致模型训练的效果有差异。 在更大范围的研究中, 考虑到目视解译遥感影像采样的方式耗费较多人工, 且易产生主观判别错误, 可使用统一时间范围的各相关数据集叠加, 筛选后进行采样数据的选取。