付杰 ,宋伦,2 ,于旭光,2 ,雷利元

(1.辽宁省海洋水产科学研究院 大连 116023;2.辽宁省海洋环境监测总站 大连 116020)

0 引言

《中华人民共和国海岛保护法》中的海岛是指“四面环海水并在高潮时高于水面的自然形成的陆地区域,包括有居民海岛和无居民海岛”。海岛作为重要的国土资源,在维护国家安全与海洋权益、保障海洋生态环境安全等方面具有不可替代的作用。近年来,随着我国沿海经济的快速发展,近岸海岛成为海洋经济发展的前沿阵地,然而海岛资源有限,生态承载力低、生态系统脆弱,如何实施海岛及其周边海域可持续利用,是当前重大而又紧迫的任务。当今,遥感监测技术发展主要呈现“高空间分辨率”“高光谱分辨率”与“多时相/高时相”等优势,广泛应用于海岛土地利用监测领域。隋玉正等[1]利用SPOT-5 影像结合实地调查,分析了2004年、2008年 和2010 年浙江洞头县43个海岛填海造地时空变化。李利红等[2]以西门岛SPOT-5影像为数据源,通过不同单尺度纹理指数结合光谱信息得出不同精度的土地利用分类,筛选出最佳波段和最佳纹理指数。杨曦光等[3]探讨利用决策树分类器对镆铘岛土地利用进行分类,相较其他分类器得到较好的分类结果。王忠芳等[4]利用南威岛QuickBird遥感影像,通过加权均值方差法和最大面积法确定了岛礁上各类典型地物的最优分割尺度。上述方法代表遥感图像分类的4个阶段:人工解译、基于像素的光谱自动分类、基于像素的专家决策树分类和面向对象图像分类。人工解译需要一种先验知识,在遥感图像寻找对应关系。基于像素的图像分类以像元为研究尺度,主要利用光谱特征,使用统计学方法、支持向量机、最大似然法等进行分类。但是这些方法存在明显局限,首先单纯地采用光谱特征未考虑其纹理特征和几何特征,造成分类精度低下;其次分类器需要处理大量的影像数据,效率低下;最后未考虑影像的空间结构信息和相邻像元之间的影响关系,造成“椒盐”及错分现象,直接影响分类精度。面向对象分类对原有影像按照一定的规则进行分割,从而提取同质对象,将最小的工作单元由像元变成同质对象。用同质对象间的关系反映地面实体间的联系,充分利用影像的光谱、纹理、形状等特征,能够十分有效地保持地类的完整性,是目前主流的分类方法。

面向对象分割尺度参数的选择直接决定着同质对象斑块的大小。参数设置上,如果选择的分割尺度过大,分割得到的对象数量就越少,易将不同类型的地物分割成为一个同质对象,降低分类精度,地物类别处于“欠拟合”状态。反过来,选择的分割尺度越小,分割得到的对象数量就越多,大大降低了运算效率,地物处于“过拟合”状态。如何智能高效地选择最优分割尺度来提高海岛土地利用的分类精度是个值得研究的问题。

分类器的选择对于分类结果也是至关重要的,针对简单地物类别,拟合能力强泛化能力弱的分类器就可以通过较少的训练样本得到较好分类结果。反之,如果是复杂地物,那么分类器学习就需要大量的训练样本和泛化能力强的学习算法。好的分类器能够根据分类对象复杂度和训练样本的数量自动地调整拟合能力和泛化能力之间的平衡。

本研究提出一种局部方差变化率算法计算海岛地类所对应的多级最优分割尺度,借助随机森林算法分类器对海岛土地利用进行智能化分类,目的是提高分类效率和精度,为开展海岛土地利用适宜性、开发强度和生态承载力评价提供前期技术参考。

1 区域概况与数据源

1.1 研究区概况

觉华岛,俗称“大海山”,又称“菊花岛”,位于辽东湾西部,连山湾与长山寺湾之间,辽宁省兴城市曹庄镇对面,距大陆最近(辽宁省兴城市曹庄镇五城子村)6.7 km,其地理坐标120°48'49.4″E,40°30'08.7″N。整体呈NE-SW 走向,形似长葫芦,两头宽阔,中间窄细。最长处约6 km,宽1～4 km,海岛岸线长245 km,海岛面积11.25 km2,是辽东湾内最大的岛屿。岛上最高点为大架山,海拔198.2 m。海岛周边海域分布有磨盘山岛、杨家山岛、张家山岛、双砬子礁岛、无名岛XCL4等岛礁。

1.2 数据源及预处理

本研究采用的遥感数据源是2018年6月28日获取的北京二号卫星0.8 m 全色波段和3.2 m 多光谱波段影像,覆盖面积62 km2。波段范围:蓝色,440～510 nm;绿色,510～590 nm;红色,600～670 nm;全色,450～650 nm;近红外,760～910 nm。多光谱与全色融合后的图像最高分辨率为0.8 m,换算成比例尺为1∶(0.8×3 779.52)=1∶3 023.616,所以成图比例尺为1∶5 000。本研究还搜集了1∶5 000 研究区内第三次土地利用现状调查矢量数据,主要用于该研究方法分类精度评价。

数据预处理首先对获取影像进行辐射定标、FLAASH 大气校正和几何配准,配准误差在一个像素之内。用搜集12.5 m 分辨率的DEM 进行正射校正,对研究区影像同时进行倾斜改正和投影差改正。为了提高影像的空间分辨率,将多光谱和与之对应的全色波段进行Gram-Schmidt算法的波段融合,该融合方法原始多光谱和辐射校正后的多光谱各自融合后的色彩效果相差极小,能保持融合前后影像波谱信息的一致性,是一种高保真的遥感影像融合方法。

2 研究方法

2.1 多尺度分割算法

图像分割是面向对象分类技术最为基础和重要的问题,尺度直接决定着分类后的精度。本研究采用多尺度分割(multiresolution segmentation)算法,能够达到给定尺度较好的地物斑块提取效果。多尺度分割是一种自下而上(bottom-up)的分割技术,其目的是实现分割后影像对象内部异质性最小化。在分割过程中,需要综合考虑光谱异质性和形状异质性,以实现分割对象形状紧凑以及边界光滑的效果。形状异质性又包括光滑度与紧致度两个因子[5]。

(1)影像对象的光谱异质性hspectrum:

式中:C为参与分割波段的总波段数;Wc为对应波段权重;σc为该波段的影像对象标准差,由构成一个影像对象的所有像素值计算得到。标准差代表了影像灰度分布的波动情况,可以用来衡量整体差异。

(2)影像对象的形状异质性hshape:

式中:hsmooth代表光滑度异质性参数;hcompact代表紧致度参数;W代表各自的权重值,二者和为一。光滑度是对象的周长l与最小外包矩形周长b的比值,用来表示对象的平滑程度,是衡量对象规则与否的一种指标。紧致度是对象的周长l与对象大小(对象的像素数n)的平方根的比值,衡量对象接近圆的程度。若光滑度权重越高,分割后的图像斑块边界就越光滑,反之,紧致度权重越高,分割后的图像斑块边界就越接近矩形。

(3)影像对象合并前后的异质性:

式中:nobj1、nobj2为两个相邻的子对象像素数,nmerge为合并新生成的对象像素数;为两个相邻的子对象标准差为合并新生成的对象标准差。得出结果为对象合并后得到的光谱异质性值。

(4)影像对象的总异质性f:

f=Wspectrum×hspectrum+Wshape×hshape

式中:对象的总体异质性f由光谱异质性和形状异质性的加权平均值所构成。Wspectrum和Wshape代表光谱和形状的权重值,二者之和为一。

2.2 最优分割尺度选择

多尺度分割仅仅是高分辨率影像解译的第一步,针对海岛每一种的土地利用类型,理论意义上都会存在一个最佳尺度与之相吻合,尚不存在单一的普适尺度满足所有地类,智能快速高效获取最佳尺度仍是面向对象影像分割的研究热点。过往研究中,有不断反复“试错法”,通过目视解译定性判断不同地类的最佳分割尺度;于欢等[6]提出了矢量距离指数法来解决面向对象遥感影像分类中的最优分割尺度选择问题,取得较好效果;后续国外学者Drǎguţ等[7]引入局部方差变化率(ROC-LV)并开发了尺度参数估算2.0(ESP2)模块,是对最早采用局部方差(LV)法的改进,用于评价不同尺度分割结果整体最大异质性。该模块能直观显示局部方差变化率随分割尺度变化曲线,曲线对应的峰值指向了可能的最优分割尺度。由于影像的复杂性,最优分割尺度往往存在很多个,这是针对影像所包含不同地物类别得出的。本研究利用ESP2确定多个最佳分割尺度,目的使得分割对象内同质性最大且对象间异质性最大。LV和ROC-LV的计算公式为:

式中:CL为单个影像对象在第L波段的亮度值;为影像所有对象在第L波段的亮度均值;m为影像对象总个数;LVL-1表示将L目标层当作基准的下一层局部方差。

2.3 随机森林算法

最佳尺度分割得到影像对象作为面向对象影像分类的最小单元,后期必须经过不同地物类别样本选取、最优特征筛选、使用最优特征子集训练样本、选择合适分类器对多尺度影像对象进行分类。分类器的选择也是本研究的关键。目前主流的监督分类器有:贝叶斯、K最近邻、支持向量机、决策树、随机森林和深度学习等。针对特定分类问题,选择最适合分类器的最佳途径是测试算法,即选择在交叉验证中表现最好的,首要考虑的是分类精度、训练和分类速度、易用性和原始数据集大小等因素。支持向量机、随机森林和K最近邻精度高,效率不是很高,不利于解释;决策树精度和效率适中,容易发生过拟合,可解释性好,易用性好;贝叶斯能对大规模数据进行分类,效率高,但需要计算先验概率,决策分类存在错误率;深度学习分类精度高,可移植性好且学习能力强,但模型过于复杂、计算量大和对软硬件要求较高,目前在AI领域很有前景,属于比较“奢侈”的分类器。本研究采用随机森林(RF)分类器对基于最佳尺度的海岛土地利用类别进行分类。RF 分类核心思想[8]是以指定决策树数量(Ntree)为基础的集成分类器。分类过程是对所有决策树的决策结果使用“少数服从多数”原则进行打分,得到最终的分类结果。袋外数据(OOB)误差决定分类精度。

随机森林方法主要包括训练与分类两个过程。生成RF的步骤如下[9]:①训练过程通过Bootstrap自助抽样方法从训练样本中有放回地随机抽取K个样本子集,构建K棵树,没有抽中的样本作为袋外数据,共计K个袋外数据。②每一棵树的所有节点有放回的抽取m个(m小于对应样本子集中总特征数量M)特征,通过计算每一个特征所包含的信息量,在m个特征中选择一个最具分类能力的特征进行节点分裂。③每棵树最大限度生长,不做修剪。④由生成的每棵树组建随机森林,用随机森林对没有参与采样数据进行分类,按照如下算式产生分类结果:

式中:H(x)为RF 最终分类结果;hi(x)为单棵树分类结果;Y为输出变量;I为示性函数。分类后,结合验证样本得到每一地类分类的Kappa系数和总体精度。采用多尺度可连续性技术,保证相邻地类边界的一致性,最终所有土地利用分类结果呈现在一幅图上。采用聚类技术,对分类图像进一步完善,技术流程如图1所示。RF 具备复杂地物分类的能力,对于噪声和存在缺损值的数据具有良好的鲁棒性,同时具有较快的学习速度,其对多维特征数据重要性程度进行度量,依据特征贡献率进行最优特征筛选,从而达到对高维特征空间进行降维的目的。相较当前流行的分类器具有较高的准确性和稳健性。因此,本研究使用RF进行分类。

图1 技术流程

3 结果与分析

3.1 最佳分割尺度获取

执行多尺度分割需要首先设置形状因子和紧致度因子,针对森林覆盖度较高区域,过往研究得出在因子0.5 附近[10],分割结果更为符合研究区覆盖特征。本研究在50～545尺度参数范围内开展以单位步长为5 进行递增的多尺度分割测试,将形状因子和紧致度因子参数设置为3种情况来计算最优分割尺度:①Wshape=0.5,Wcompact=0.4;②Wshape=0.5,Wcompact=0.6;③Wshape=0.4,Wcompact=0.6。局部方差及变化率随分割尺度变化曲线见图2。由图2 可以得出:Wshape=0.5,Wcompact=0.4,对应最佳尺度为75、135、230、415;Wshape=0.5,Wcompact=0.6,对应最佳尺度为65、90、180、520;Wshape=0.4,Wcompact=0.6,对应最佳尺度为80、160、355、465。

图2 多分割尺度LV、ROC变化曲线

经目视解译,在影像上选取7处位置依次分别对应典型7类地物类型:淤泥质海岸、农用地、林草地、建设用地、针叶林地、水体和基岩海岸。利用3组形状因子和紧致度因子对应最佳尺度对研究区影像进行分割,为保证同一组内不同尺度分割结果得到对象边界的连续性,首先需要实现在最小尺度下分割出若干子对象,然后不断向上(Create above)重新聚类实现较大尺度下分割出中间对象,最后在最大尺度下生成一个父对象(图3至图5,表1)。

表1 最佳尺度分割层次

经过实地踏勘,结合图3、图4、图5和表1可以得到:当分割尺度小于160 时,影像产生对象均在15 000以上,在7类典型地物类别中,真实地物斑块要比分割后的影像对象大,地物类别处于过拟合状态,代表地类有淤泥质海岸、水体和基岩海岸;尺度大于180时,其他4类典型地类分割数量骤减,真实地物斑块要比分割后的影像对象小,地物类别处于欠分割状态。因此,图3至图5中农用地应继续细分为有茬(长有农作物)耕地和无茬耕地。林草地应细分为落叶阔叶林地和草地。建设用地应细分为建筑物、道路和裸地。上述确定海岛地类的核心思想是:选择经ESP2模块计算的最优分割尺度,再将分割结果对象与影像进行叠加分析,通过实地判

图3 Wshape=0.4,Wcompact=0.6

图4 Wshape=0.5,Wcompact=0.4

图5 Wshape=0.5,Wcompact=0.6

读研究分割对象与地类斑块两者的边界形状、面积大小以及分割个数的吻合程度,从而判定该尺度针对该地类分割结果的好坏。如果一种地类,随着分割尺度减小,生成对象数量急剧上升,表明该地类“不纯”,需要细分。为便于对比研究不同地类对应最佳尺度,本研究将觉华岛及周边岛屿划分11个土地利用类别,分别为建筑物、道路、裸地、常绿针叶林地、落叶阔叶林地、灌草地、有茬耕地、无茬耕地、粉砂淤泥质海岸、基岩海岸和水域。

3.2 样本对象最优特征筛选

本研究共提取3类特征,分别是光谱、纹理和几何特征。

(1)光谱特征包含所有波段像元亮度均值(brightness)、标准差(standard deviation)、最大差分(Max.diff.)、每个波段像元亮度均值(mean)、贡献率(ratio)、归一化植被指数(NDVI)、归一化水体指数(NDWI)、土壤指数(NDSI)和建筑物面积指数(BAI),共4个波段累计18个特征。

(2)纹理特征使用灰度共生矩阵(GLCM)来描述,选取特征有均质性(homogeneity)、对比度(contrast)、相异性(dissimilarity)、熵(entropy)和角二阶矩(Ang.2nd moment),共4个方向累计20个特征。

(3)几何特征是基于构成影像对象的像素的空间分布统计之上的。使用协方差矩阵作为统计处理工具。选取特征有紧致度(compactness)、密度(density)、面积(area)、矩形适合性(rectangular Fit)和长宽比(length/width),共5个特征。

上述选择了43个针对海岛地物类别的典型特征。对于分类器训练样本来讲,并不是特征越多分类精度越高,而是特征越多模型越容易过拟合,分类效率和精度越低。因此需要从过多特征中筛选出对地物类别识别作用最大的特征子集。如何高效智能地对高维特征空间精准降维是当前研究的热点,国内外学者认为[10-11]分离阈值法(SEa TH)是目前比较有代表的特征优化方法,该方法最大优点是:①能够自动选择特征,自动确定分离阈值;②J-M 距离适用于专题信息提取,以及类别数目较少的情况(两两类别区分)。当然此法仍需学者不断改进[12-13],这里不再赘述。本研究依照SEa TH,结合电脑硬件配置综合考虑,筛选出13个特征作为最优特征组合,按照特征贡献度依次是:ratio(Blue)、length/width、Max.diff.、density、GLCM entropy(0)、ratio(Red)、mean(Nir)、compactness、mean(Blue)、NDVI、brightness、ratio(Nir)、mean(Green),能有效提高分类精度及减少分类时间。

3.3 分类结果与精度评价

前期最佳分割结果产生的对象是地物分类的基础。经实地踏勘,确定11个土地利用类别后,需建立遥感信息解译标志。研究区内海岛植被占比较大,主要植被类型是农作物群落、常绿针叶林、落叶阔叶林和灌草丛。农作物种为玉米和果蔬,分布在海岛居民区附近的平地;常绿针叶林主要是油松林、油松刺槐林,以油松纯林为主;落叶阔叶林为刺槐林,分布在平缓低地;灌草丛包括崖椒灌丛和荆条灌丛,以荆条灌丛为主,在林缘坡地和海岸坡地处分布。针叶林在北京二号卫星R4G3B2 假彩色影像上呈深暗红色及暗红色,纹理呈不规则的条带或面状,纹理较粗糙,有立体感。阔叶林在R3G2B1真彩色影像上呈亮绿色,与针叶林比较,色调更亮一些。农作物群落划分在有茬耕地类中,在影像中形状规则,呈格状、弧形或条带状,纹理粗糙,假彩色中呈深灰、浅红色或鲜红。灌草地是磨盘山岛、杨家山岛和张家山岛的主要植被类型,连片分布,由于覆盖度不同,在假彩色中分别呈暗红、红色、浅红和浅黄绿色。觉华岛海岸带主要有基岩和粉砂淤泥质海岸,基岩坡度较大,宽度较窄,干燥的部分亮度较高,常被海水冲刷部分呈暗灰色。淤泥质形状呈带状分布,宽度不等,从数百米到数千米,坡度平缓,有潮沟发育,影像上多呈树枝状或蛇形状,假彩色上呈灰黑色条带,色调较暗,无立体感。

随机森林重要参数的确定:参考刘敏等[14]利用加州大学欧文分校(University of CaliforniaIrvine,UCI)数据集对随机森林中决策树的数量与数据集的关系进行实验分析,实验结果表明当树的数量为100时,就可以使分类精度达到要求。顾海燕等[9]实验得出在特征数量保持常量为5时,随着树的数量的增加,误差变小;当树的数量为200时,误差达到最小;当大于200时,误差差别很小,趋于稳定,计算时间随着树的数量增加而增加。考虑到运算效率与分类精度平衡,本研究实验设计如下:决策树设置范围为(150,200,250),特征数分别为(11,12,13)时,参数配对9组,通过样本精度验证来寻找参数最佳组合。最终确定使用筛选好的13个特征对520、465、415、230、180、160这6种尺度选取一定数量样本进行训练,考虑到运算效率与分类后精度可对比性,将随机森林中决策树的数量Ntree选定为200。本研究通过不同尺度分类结果的混淆矩阵提取4类指标来进行总体和某一类地物的精度评价,分别为总体精度(overall accuracy,OA)、Kappa系数、生产者精度和用户精度。6组分割尺度对应整体精度评价如表2 所示,11 种地物分类精度与Kappa系数如表3所示。

表2 6组分割尺度对应整体精度评价

表3 6组分割尺度对应11种地物类别精度评价汇总

研究区内建筑物对应最佳分割参数为(0.5,0.6,180),单类Kappa系数为0.735 7;道路对应最佳分割参数为(0.5,0.4,230),单类Kappa系数为0.752 8。本研究建筑物提取精度不高主要原因是高分辨率影像形状细节特征冗余,阴影遮挡,屋顶材质不同导致光谱信息多变从而影响建筑物提取精度。目前高分辨率影像建筑物精准获取是研究热点问题,国内一些学者进行了深入研究[15-18]。不同分割尺度下道路提取的生产者精度较低而用户精度均高于83.87%,道路验证样本被错分到建筑物、裸地和有茬耕地的比例较高,其他验证样本错分到道路的比例较低。研究区内主次干道的光谱特征和几何形状与局部建筑物屋顶和裸地表面呈现极大相似性,支路由于宽度较窄,与有茬耕地在160以上的分割尺度难以区分。水域包含养殖的坑塘水面和岛屿周边海域[19],对应的最佳分割参数为(0.5,0.4,230),83个验证样本全部得到正确分类,与海域邻接的13个海岸带地类被错误划分到水域中,导致用户精度降低至86.46%。裸地在假彩色影像上色调呈青灰或亮白色,形状为块状、条带状,纹理上有凸凹不平的立体感,对应最佳分割参数为(0.4,0.6,465),提取用户精度低是由于43个建筑物、道路和基岩样本被错误划分到此类中,类别特征空间距离较近。其他作为海岛特色的地类,采用本研究分类器,提取精度较高:①海岛植被,其中常绿针叶林地对应最佳分割参数为(0.5,0.4,415),单类Kappa系数为0.923 2;落叶阔叶林地对应最佳分割参数为(0.4,0.6,465),单类Kappa 系数为0.895 0;灌草地对应最佳分割参数为(0.5,0.4,230),单类Kappa系数为0.889 0。海岛植被提取的生产者精度均在90.40%以上。②海岸带,其中基岩海岸对应最佳分割参数为(0.5,0.6,520),单类Kappa系数为0.909 2;粉砂淤泥质海岸对应最佳分割参数为(0.4,0.6,465),单类Kappa系数为0.979 9。海岸带提取的生产者精度均在92.07%以上。③耕地,其中有茬耕地和无茬耕地对应最佳分割参数为(0.5,0.6,180),单类Kappa系数分别是0.858 7和0.915 3。耕地提取的生产者精度均在87.70%以上。

4 结论

使用高分辨率遥感影像对海岛土地利用类型进行遥感分类研究必然涉及最优分割尺度和分类器选择问题,不同地类要求使用与之对应的绘制尺度,一种尺度不可能满足所有地物类别划分。本研究通过对候选最佳分割尺度范围内的样本参与随机森林分类所取得的分类精度进行评价分析,得到以下结论。

(1)基于高分辨率的北京二号卫星融合后0.8 m分辨率数据,采用随机森林分类进行觉华岛及周边岛屿土地利用类型分类,当形状因子、紧致度因子和分割尺度3参数组合为(0.5,0.6,180)时,整幅图分类精度表现最佳,总分类精度达到81.73%,总Kappa系数为0.798 0。

(2)本研究针对海岛特色地类,选择最佳分割尺度,筛选特征贡献率最高的13个特征,随机森林中决策树的数量Ntree为200,可以得到单个地类分类精度达到效果最好。海岛植被提取的单类Kappa系数为0.889 0以上,生产者精度均在90.40%以上;海岸带提取单类Kappa系数为0.909 2以上,生产者精度均在90.40%以上;耕地提取单类Kappa系数为0.858 7以上,生产者精度均在87.70%以上。

(3)本研究针对高分辨率建筑物类别提取精度不高。主要原因是高分辨率影像形状细节特征冗余,阴影遮挡,屋顶材质不同导致光谱信息多变。