基于SNIC-CNN-SVM模型的京津风沙源二期工程区土地利用/土地覆盖遥感识别研究

李长龙,李增元,高志海,孙 斌,*,王丝丝

1 广州商学院信息技术与工程学院, 广州 511363 2 中国林业科学研究院资源信息研究所, 北京 100091 3 国家林业和草原局林业遥感与信息技术重点实验室, 北京 100091 4 国家遥感中心, 北京 100084

土地覆被信息不仅是单一地表类型的表达,而且体现一定区域内生态环境和人类活动等自然和非自然属性的综合特征[1]。因此,不同的土地覆盖类型代表着地球表面丰富多样的生态系统和自然人文景观,它是生态系统内部长期动态演变的积累,也是量变到质变的外部反映[2—3]。通过对土地覆被长时间序列变化研究能够分析不同生态系统内部生物和非生物的综合生态效应,从而掌握生态系统的长期演变规律,为全球气候变化及其人类活动影响等提供重要信息支撑[4—6]。

随着遥感技术的飞速发展,特别是深度学习算法的不断优化,依靠遥感技术自动半自动地提取土地覆被信息已然成为当前的重要手段[7—9]。近年来,在国家和区域尺度,部分学者研发了多套地表覆盖遥感产品,例如刘良云团队GLC_FCS30[10]、宫鹏团队FROM_GLC30[11]和FROM_GLC10[12]、国家自然资源部GlobeLand30[13]、欧空局全球陆地覆盖数据(ESA GlobCover)[14]等,这些都极大地提高了气候变化和生态环境评估等综合性、大尺度深度研究的效率,为地理国情监测等提供了重要的基础信息。

稀疏植被区主要是指植被覆盖度在10%—50%的地表区域,包括草原、苔原、荒漠草原、荒漠、岩溶地区和部分高海拔植被区等。其中,零星的植被由草本和(或)木质和半木质植物组成,其余区域为自然裸地[15—16]。京津风沙源治理二期工程区从东南向西北降水逐步减少,稀疏植被分布广泛。由于稀疏植被区内的自然裸地多为反射率较高的沙砾或岩石,而植被又零散分布在自然裸地内部,这就造成稀疏植被的光谱信息易被淹没,提取难度较大[17—18],而针对不同地类的有效性特征是提高地表覆盖制图精度的重要因素[19]。目前,全球或全国的地表覆盖产品尺度较大,针对稀疏植被区的特异性较差,类型识别精度较低,而且缺乏对不同沙化土地类型的精细划分[20—21]。这些都极大地限制了稀疏植被区生态效应评价、系统管理和科学决策等方面的应用需求,因此迫切需要建立基于大尺度的稀疏植被区土地覆被信息的遥感提取方法。

根据遥感分类的影像基础可以将土地覆盖分类方法细化为基于像元、基于超像素和面向对象三类[22]。基于像元的方法较为传统,主要依靠影像的光谱信息对不同地物进行监督或非监督的聚类,但是缺乏对几何和周围影像信息考虑[23—25]。面向对象的方法主要有分水岭变换算法[26]、均值漂移算法[27]、分形网络进化算法[28]和马尔可夫随机场分割算法[29]等,分类的基础是对遥感图像信息提取后分割成的对象,图像分割优劣直接影响分类精度,不同地物类型对应不同分割尺度,计算量较大,效率较低,最优分割尺度的选择多以人工或人机交互方式进行,不利于建立自动化提取模型[30—31]。超像素分割则弥补两者缺点,通过对像素进行相似性分组,从而形成相对规则和均匀大小对象,以类似过分割的方式对图像信息进行解析[32],去掉了最优尺度选择和对象形状特征提取等工作量,分析作用更强,更有利于沙化土地类型中植被和沙地边界模糊或渐变性边界的提取和划分,目前常用的算法包括基于图论的算法和基于梯度下降的算法等[33]。以卷积神经网络(Convolutional Neural Networks,CNN)为代表的深度学习方法解决了最邻近、极大似然等传统遥感分类方法在遥感图像中隐藏的深层次语义信息缺乏的问题,能够从本身的深层结构中学习高度抽象的图像语义区别特征[34—35],对于稀疏植被区土地覆盖分类具有较强的适用性,能够一定程度解决该区域类型复杂,变化多样,植被分布很不均一,渐变型或混杂型地类较多等难题。

综上,本研究以京津风沙源治理二期工程区为研究区,基于多源遥感数据,通过建立惩罚性机制优化超像素分割算法,研究构建用于稀疏植被区典型土地利用/土地覆盖分类识别的SCS模型,以期提高稀疏植被的遥感信息提取能力,获取科学准确的土地覆被类型信息,为研究区生态工程效应评价奠定基础。

1 研究区概况

本研究以京津风沙源治理二期工程区为研究区,该区域东西距离约1300km,南北距离约1000km,总面积约为70.98万km2,包括陕西省、内蒙古自治区、北京市、天津市、山西省和河北省的共计138个县(旗、市、区),地理坐标为105°8′—120°59′E,36°48′—46°47′N,具体位置如图1所示。

研究区地貌由平原、山地、高原三大部分组成,地形起伏较为复杂,最低处仅有几十米,最高处超过了3000m。区域内的沙化土地聚集分布,形成了“三沙地一沙漠”的格局(科尔沁沙地、浑善达克沙地、毛乌素沙地和库布齐沙漠)。研究区受到降水、气温和太阳直射时间等自然条件的影响,从东南向西北的气候分异明显,依次包括有亚湿润区、干旱亚湿润区、半干旱区和干旱区,植被则依次分布有森林草原、典型草原、荒漠草原、草原化荒漠和荒漠。由于降水和地形等自然因素的影响,研究区中部和西北部呈现明显植被和荒漠交错的破碎化自然景观,属于典型的稀疏植被区,其中植被沿低洼处或河流周围聚集性分布,种类多样且林灌草多混杂生长；裸露地表的差异性也较大,依次分布有栗钙土、砂质土、粗沙、砾石等。另外,再叠加人类活动和林农牧交错区影响,多种地表类型在一定区域内共存或多次更迭,光谱信息均一性较差,遥感信息提取难度较大[36—37]。

2 研究方法

2.1 数据获取与处理

本研究主要获取的遥感数据为Landsat影像,其中,2000年和2010年采用的是Landsat 5影像,总景数分别为274景和255景；2020年采用的是Landsat 8影像,总景数为322景,数量分布如图2所示。每一年影像均选取研究区6—9月生长季影像,然后基于像元的最小云量进行影像合成,从而最大程度减少云对影像分类的干扰。影像的波段选择为蓝、绿、红、近红外、短波红外和热红外,均采用经过大气校正和几何校正等预处理后的数据。

为了能够更好地识别沙化土地,需要去除表层植被的干扰和影响,因此本研究采用SAR数据穿透表层植被,强化沙化土地的识别能力。根据实验分析,极化方式为VV+VH组合的双极化影像对沙化土地的敏感性更强,因此数据选择Sentinel- 1A中的GRD一级产品,极化方式为VV+VH组合的双极化,成像时间为2020年7月15日。数字高程模型(DEM)数据来自NASA的地球观测系统数据信息系统,空间分辨率30m,成像时间为2010年,包含高程、坡度、坡向三类数据。为了与Landsat数据进行匹配,SAR数据经正射校正和重采样,DEM数据经过尺度转换,均统一到空间分辨率为30m。真值数据的来源主要分为两部分：野外调查和高分辨率遥感影像判读,各数据的数量和年份分布如表1所示,空间分布如图1所示。每年份的野外调查数据采样时间均为植被的生长季(7月或8月),影像判读数据则是基于Google Earth上更高分辨率遥感影像进行。

表1 真值数据来源和年份分布

2.2 研究方法

本研究在原有超像素分割算法——简单非迭代聚类(Simple Non-Iterative Clustering,SNIC)的基础上,通过建立惩罚性机制,将更多的植被光谱信息引入影像分割中,增强稀疏植被遥感信息提取能力,然后利用SAR数据本身的穿透能力建立不同纹理特征对沙化土地进行识别,最后基于CNN机器学习方法和SVM分类器对整个研究区进行不同类别的确定,这样就形成了SCS稀疏植被区遥感分类模型,系统地将(稀疏)植被、沙化土地和其他类别进行精准划分。具体主要分为五步：分类体系的建立,特征选择,超像素分割,影像分类和精度评价。

第一,分类体系的建立。参考自然资源部《土地利用现状分类》(GB/T 21010—2017)和《第三次全国国土调查技术规程》(TD/T 1055—2019)分类标准,综合Landsat遥感影像土地覆盖信息提取能力,研究区分类体系定为8大类,分别是：乔木林、灌木林、草地、沙地、耕地、水体、戈壁和其他,含义如表2所示：

表2 土地覆盖分类体系及含义

第二,特征选择。主要分为光谱特征、指数特征和纹理特征,其中光谱特征为影像数据的波段信息；指数特征除了归一化水体指数(NDWI)和归一化植被指数(NDVI)外,增加了土壤调节植被指数(SAVI),其目的是去除稀疏植被区(特别是戈壁)背景土壤亮度较高的影响,三者的计算公式如下：

(1)

(2)

(3)

式中,NIR、Red和Green分别为近红外、红光和绿光波段的光谱值；L为植被密度参数。纹理特征是基于SAR数据(Sentinel- 1A)形成的七种,分别为差异性、对比度、方差、角二阶矩、逆差矩、熵和相关性,示意图如图3所示。

图3 SAR影像的7种纹理特征示意图Fig.3 The seven texture features of SAR imageSAR：合成孔径雷达 Synthetic Aperture Radar

第三,超像素分割。SNIC算法是将影像转化为CIELAB颜色空间的五维特征向量,表征空间距离,即位置信息X、Y的两维数据；表征颜色距离,即色彩信息的三维数据,最后基于这五维特征向量建立距离度量标准,进而对图像进行局部聚类的过程。算法的实现过程主要有四步：生成种子点；种子点优化；最优距离的求解；优先队列算法像元聚类。其中最优距离的求解如公式4—6所示,优先队列算法则是基于K个种子点建立K个元素ei={(xi,yi),(li,ai,bi),k,Di,k},这K个元素建立最原始的优先队列Q,然后Q每次返回队列中到第k个种子点的距离Di,k最小的元素ei,这样可以将多次迭代优化成一次迭代,极大地提高运行效率。

(4)

(5)

(6)

式中,Ds、Dc和D分别代表颜色距离、空间距离和最优距离,s和m为颜色距离和空间距离的归一化参数。

原有SNIC算法将空间距离和颜色距离作为同等重要参数进行计算,更易形成较为规则的分割图像,这样在稀疏植被区土壤和植被的交界处,高亮的背景土壤会淹没稀疏植被本身的光谱信息,造成影像分割误差。因此,本研究在原有算法后增加第五步,惩罚性机制和连通性检查。即对于某个属于种子点a的像元jk,假设它的空间距离参数属于种子点a,而颜色距离参数属于种子点b,且满足下面公式：

Da≥ε1,Db≤ε2

(7)

式中,Da和Db分别为jk属于种子点a和种子点b的类内距离,ε1和ε2为最大类内距离。经连通性检查,如果像元jk与种子点b的类内像元相连接,那么该像元jk将被标签为种子点b的类内元素。

第四,影像分类。运用CNN算法挖掘深层特征,支持向量机(Support Vector Machine,SVM)进行精确分类,分类过程如图4所示。CNN算法通过两次卷积,两次采样和一次全连接算法,获得卷积层权重参数,将一维特征影像转换成256维。假设在CNN的卷积层中,输入的一维特征向量表示为xk-1,则其经过卷积计算加上偏置参数的非线性变换,输出特征向量为xk,计算过程可以表示为：

xk=f(wk×xk-1+bk)

(8)

式中,wk和bk表示卷积层k的卷积核参数(5×5)和偏置参数；*表示卷积计算符号,f()表示激活函数,本研究采用ReLU函数f(x)=max(0,x)进行非线性计算。

图4 CNN-SVM算法分类过程示意图Fig.4 The classification process of CNN+SVM algorithmCNN：卷积神经网络 convolutional neural networks；SVM：支持向量机 support vector machine

第五,精度评价。采用混淆矩阵(Confusion matrix),包括生产者精度(Producer accuracy)、用户者精度(User accuracy)、总体精度(Overall accuracy)和Kappa系数(Kappa coefficient)。

3 结果与分析

3.1 超像素分割结果

为了呈现超像素分割过程优化效果,选取浑善达克沙地腹地部分区域,对比了原始影像、SNIC算法优化前后的分割结果(图5)。结果显示,优化的SNIC算法在地类内部分割效果与未优化前相比,变化不大。但在不同地类间,特别是在稀疏植被和沙化土地类别的边界(图5的红圈部分)区域,原有算法区分度不够,在多个位置的地类界限具有明显的混杂现象,而优化算法有效解决了这一问题,使分割效果更加清晰,能够准确细致地对稀疏植被和沙化土地的边界进行刻画,有利于稀疏植被和沙化土地信息的区分。

整体而言,对于整个工程区中诸如沙地和草地、沙地和戈壁、草地和灌木林地、草地和森林等易混分的土地覆盖类型,优化算法的遥感信息提取能力更强,分割效果更好,特别是各易混类别的边界区分度更高,这样更有利于下一步影像分类,进而提升最终的信息提取精度。

图5 SNIC算法优化前后的影像分割对比Fig.5 Comparison of image segmentation before and after SNIC algorithm optimizationSNIC：简单非迭代聚类 Simple non-iterative clustering

3.2 分类精度评价

优化前后SCS分类结果的混淆矩阵分别如表3和表4所示。结果表明,基于原有SCS模型的工程区典型土地利用/土地覆盖的分类识别总体精度为78.24%,Kappa系数为0.7452,而经过优化,模型的总体分类精度显著提高,达到89.41%,Kappa系数为0.8760,总体提高了11.17%,Kappa系数提高了0.1308。其中对于工程区面积最大的两种土地利用/土地覆盖类型——草地和沙地,两者的用户者精度分别提高了14%和9%。综合整体分类精度和典型地类单一识别精度结果发现,本研究提出的SCS模型优化方案能够有效提高稀疏植被区域的植被与沙地信息提取精度,能够满足大区域尺度对稀疏植被区生态效应、环境变化和人类活动等深入研究的需要。

表3 2020年原有模型土地覆盖分类结果混淆矩阵

表4 2020年优化模型土地覆盖分类结果混淆矩阵

3.3 2020年工程区典型土地利用/土地覆盖空间分布

采用优化的SCS模型,获得了2020年工程区典型土地利用/土地覆盖的空间分布情况(图6)。结果表明,总体上,研究区多以草地、林地和沙地为主,其中草地面积最大,占研究区总面积的51.52%,广布于研究区的大部,东北部的乌珠穆沁草原区域是最重要的分布范围；林地范围次之,乔木林和灌木林分别占研究区总面积的19.55%和7.23%,分布范围上来说主要位于河北省北部和内蒙古锡林郭勒盟南部以及零星分布于沙地和草地周边；沙地则主要集中于一个沙漠三个沙地区域(库布齐沙漠、毛乌素沙地、浑善达克沙地和科尔沁沙地),总面积占研究区总面积的11.96%；耕地主要分布于黄河沿岸和北京市、天津市以及河北山西的平原上,约占研究区总面积的1.54%；戈壁则在研究区的西北部乌拉特后旗,约占研究区总面积的5.20%。稀疏植被覆盖(草地、沙地、戈壁)区域面积占比约为68.68%,研究区大部分区域属于脆弱生态区。

图6 研究区2020年土地利用/土地覆盖类型分布Fig.6 The distribution of land use/land cover types in 2020

3.4 近20年工程区典型土地利用/土地覆盖类型空间变化格局

自2000年以来,研究区相继实施了退耕还林还草工程、野生动植物保护及自然保护区建设工程、天然林保护工程、京津风沙源治理工程等生态保护重大工程[38],在森林面积增加,沙化土地减少,区域经济发展等方面取得了丰硕的成果,充分发挥了生态修复和改善,维护国家生态安全等方面的重要作用。由于京津风沙源治理工程(一期)实施时间为(2001—2010年),二期工程实施时间为(2013—2022年),因此本研究选择2000年和2010年作为对比时间节点,深入分析近20年京津风沙源治理二期工程区土地利用/土地覆盖类型变化情况。图7所示为2000年和2010年土地利用/土地覆盖类型分布,图8所示为近20年工程区典型土地利用/土地覆盖类型的变化情况。结果表明,整体上,近20年来乔木林、灌木林和草地的面积增加最多,其中,乔木林占研究区的比例由2000年的10.17%增加到2010年的19.22%,再到2020年的19.55%,增加了9.38%；灌木林的面积则从1.43%(2000年)增加到6.65%(2010年)再到2020年的7.23%,增加了5.79%；草地的面积由46.05%(2000年)下降到44.26%(2010年)再增加到51.52%(2020年),增加了5.47%。与之相应的是,耕地和沙地的面积明显减少,两者占工程区总面积的比例分别由2000年的13.33%和16.54%,减少到2020年2.43%和11.96%,分别减少了10.89%和4.58%。

图7 研究区2000年和2010年土地利用/土地覆盖类型分布Fig.7 The distribution of land use/land cover types in 2000 and 2010

图8 研究区2000—2020年土地利用/土地覆盖类型变化统计Fig.8 The statistics of land use/land cover changes in the study area from 2000 to 2020

在自然条件下,土地覆盖类型的变化是量变到质变的过程。植被逐步退化最终变为荒地或沙地,而沙化土地随着植被的逐步改善,最终变为草原或灌木林[39]。因此,本研究进一步从乔灌草植被覆盖区和沙化土地两方面的相互转化情况来反映近20年来工程区的实施效果。变化情况如图9所示。从2000年到2020年,位于浑善达克沙地、毛乌素沙地和科尔沁沙地的腹地的沙化土地相对稳定,而沙地边缘植被恢复明显,特别是杭锦旗、达尔罕茂明安联合旗等沙地边缘地带,沙地转换区域中的86.45%转换为乔灌草植被覆盖区；由其他土地利用/土地覆盖类型转为沙地的区域主要集中在锡林郭勒盟的苏尼特左旗和苏尼特右旗等地,转换为沙地的区域中76.14%来自于原有乔灌草植被类型覆盖区。植被方面,乔灌草土地覆盖类型不变区域主要位于研究区北部和东部大部分区域；乔灌草土地覆盖类型减少区域主要位于浑善达克沙地和科尔沁沙地边缘地带,减少区域中的33.41%转换成沙化土地；乔灌草土地覆盖类型增加区域主要位于毛乌素沙地、库布齐沙漠区域,增加区域的25.62%来源于沙化土地。总体来说,京津风沙源治理工程实施的十多年时间内,沙化土地面积大幅度降低,乔灌草总面积大幅度上升,土地退化的趋势得到遏制,草地质量和生态环境改善明显,表明研究区植被状况不断改善,实施的各项生态工程作用显著,能够更有效地服务于多维度生态功能的发挥和评价。

图9 2000—2020年研究区植被类型(乔灌草) 和沙化土地变化分布Fig.9 The change distribution of vegetation types (arbor, shrub, grass) and desertification land in the study area from 2000 to 2020

4 讨论

前人的研究结果[40—41]表明,京津风沙源治理工程区土地利用/土地覆盖变化在2000—2010年主要表现为沙地、耕地和林地的增加,草地的减少,2010—2018年则表现为耕地、沙地的减少,草地和林地的增加,恢复区域中以鄂尔多斯高原、浑善达克沙地和科尔沁沙地为主,这些结论都基本与本研究获得的结果一致。只是部分认为2018年研究区内草地的分布占比为55.25%,略高于本研究51.52%(2020年)而沙地的占比略低。可能的原因是本研究利用SAR数据(Sentinel- 1A)的穿透能力,将研究区内表层具有稀疏植被,但土壤为沙砾的土地覆盖类型划分为沙化土地,一定程度上减少了草地占比而增加了沙地占比,但本研究的提取结果能够更准确地表征沙化土地,具有更强的科学性和实际意义。

全球尺度方面,根据与刘良云、宫鹏和陈军等人研究的结果[10—13]对比来看,在各类别的区域分布上与本研究的结果基本一致,但因三人的研究成果均基于全球尺度,在区域尺度上并未细致研究各类别的特异性特征,因此也缺乏对京津风沙源区域复杂的植被土壤环境的深入研究,易出现植被和沙化土地错分的情况,特别是在稀疏植被覆盖的沙化土地区域,而本研究在该方面有较为明显的改善。

由于大尺度、长时间序列工程措施生态效应评价的需要,本研究采用的超像素分割算法一方面可以增加像元分类中缺少的对象信息,另外一方面又减少面向对象分割中最优尺度等人工参与,实现稀疏植被区沙化土地类型识别和土地覆盖分类的自动化、全流程的快速提取,以SAR数据表征沙化土地为基础的机器学习方法能够更加细致刻画稀疏植被区中沙化土地和稀疏植被的边界。但由于验证样本不足,本研究在部分区域的精度还需进一步进行论证,特别是本研究主要聚焦于稀疏植被和沙化土地等自然地表类型,对于建筑物、道路等人工地表类型的识别能力还需深入研究。此外,可以拓展研究该区域外非工程区土地覆盖类型的变化,从而与本研究的结果进行对比分析,进一步论证生态工程措施的积极作用。下一步,通过融合更多时序信息、专题信息等对研究区的土地利用/土地覆盖类型进行更深入的分析,进一步提高该区域的遥感识别精度。

5 结论

本研究以典型大尺度稀疏植被区——京津风沙源治理二期工程区为研究区,通过构建惩罚性机制优化超像素分割算法,研建了SCS模型方案,实现了研究区典型林草沙要素信息提取和主要土地利用/土地覆盖识别。研究结果表明：

1)引入惩罚性机制优化后的SNIC分割算法,有效提升了稀疏植被区与沙地区的边界区分度,有助于分类精度的提升；

2)基于改进SNIC-CNN-SVM模型方案的研究区总体分类精度达89.41%,较优化前提高了11.17%,特别是针对草地和沙地这两种工程区核心类型的识别精度分别提高了14%和9%,表明该模型具有较好的沙化土地和稀疏植被识别能力,该优化方案在以研究区为代表的稀疏植被区域分类中具有较好的应用效果和推广价值；

3)2020年工程区中以草地面积最大,占研究区总面积的51.52%；沙地主要集中于一个沙漠三个沙地区域(库布齐沙漠、毛乌素沙地、浑善达克沙地和科尔沁沙地),总面积占工程区总面积的11.96%；稀疏植被覆盖(草地、沙地、戈壁)区域面积占比约为68.68%,表明工程区处在林地-稀疏植被-沙地的过渡地带,生态环境保护压力与防沙治沙形势依然严峻。

4)通过对比近20年来的典型土地利用/土地覆盖类型空间分布格局发现,乔灌草等植被增加面积约占工程区20.64%,主要由沙化土地转化,沙化土地减少面积约占研究区的4.58%,沙地减少的区域中86.45%转换为乔灌草等植被类型。植被覆盖范围增加,质量提升,沙化土地面积降低,特别是在毛乌素沙地和库布齐沙漠等区域改善效果突出,表明经过20余年的治理,研究区土地退化趋势得到遏制,草地质量和生态环境改善明显,实施的各项生态工程作用显著,能够更有效地服务于多维度生态服务功能的发挥。该研究以期能够为京津风沙源二期工程区的生态系统演变规律研究及生态工程评价等工作提供重要科学支撑。