李琰 王立海 邢艳秋
(东北林业大学,哈尔滨,150040)
责任编辑:王广建。
森林类型识别是林业资源保护和林业资源动态监测的基础。高光谱遥感具有较高的光谱分辨率,被广泛应用于各个领域。高光谱遥感具有数据量大、光谱分辨率高、光谱范围广、波段窄、对地物可以连续成像的特点,在林业调查中具有重要的作用[1-2]。由于不同地物属性的先验知识不易获得,通常采用非监督分类方法,研究数据的本来结构及自然点群的分布情况,以降低处理数据的工作量[3-5]。谭炳香等[6-7]用EO-1 Hyperion高光谱数据,对吉林省汪清林业局境内森林类型进行识别,发现高光谱分类比多光谱分类效果好,但未研究聚类对其的影响;仲青青等[8]利用粒子群算法进行多光谱聚类研究,只考虑了对参数选取的改进。猫群算法(CSO)是一种模拟猫的日常行为而建立的新兴的群智能算法[9]。猫群算法收敛速度快,比传统进化算法、粒子群算法更加简洁。因此,将猫群仿生算法应用到基于高光谱遥感图像的森林类型识别上,对提高森林分类精度、减少分类时间、降低专业理解难度、减少人工干预等具有重要作用。
研究区为吉林省汪清林业局经营区,地理坐标为129°56'~131°04'E,43°05'~43°40'N。汪清林区属综合性山川地貌,地处长白山支脉老爷岭和哈尔巴岭的山麓,四周群山环绕、中间为河谷平地、群山与河谷平地之间为丘陵。海拔360~1 477 m,平均海拔806 m,坡度变化范围为0~45°。平均降水量547 mm,其中5—9月的降水量为438 mm,占全年总降水量的80%。年平均气温3.9℃。土壤种类中,山地土壤主要为暗棕色土壤,谷地土壤主要为草甸土。经营区总面积3.04×105hm2,林地面积为1.87×105hm2,属典型的天然次生林。该区资源丰富,森林类型以针阔混交林为主,带状分布。针叶树主要有红松(Pinus koraiensis)、云杉(Picea)、臭松(Symplocarpus salisb)、落叶松(Larix kaempferi);阔叶树为水曲柳(Fraxinus maudschurica)、胡桃楸(Juglans maudshurica)、蒙古栎(Quercus monglica)、椴树(Tilia)、色木(Acermono)、榆树(Ulmus pumila)、白桦(Betula platyphylla)、杨树(Populus)和枫桦(Betuladavuric)等[9]。主要纯林有蒙古栎、白桦和落叶松,占纯林树种的90%以上。
环境与灾害监测预报小卫星(HJ1A/B星)搭载的超光谱成像仪(HSI)和多光谱CCD成像仪,HJ-1A图像幅宽大于50 km,地面像元分辨率为100 m,共有115个波段。光谱范围为0.45~0.9μm,平均光谱分辨率为4.32 nm,且具备±30°侧视能力及星上定标功能。重访周期96 h,可实现对研究区的快速重复观测[10]。CCD图像4个波段,地面像元分辨率30 m。本文处理的高光谱数据为2010年6月9日的影像,HSI与CCD光谱景序列号分别为321 381、318 468。数据为2级产品,HSI中1~20以及113、114、115波段是坏数据,所以实际参与分析的波段为92个(21~88波段为可见光部分,89~112为近红外部分)。数据经过对数残差大气纠正、去条带、通过地面控制点进行几何精校正等预处理,利用参考文献[11-12]的重构和评价方法进行插值,将HSI高光谱数据影像与同时成相的CCD多光谱影像融合,形成115个波段空间分辨率为30 m的新高光谱数据。经过相对反射率的分析,典型地物HSI波谱与CCD波谱具有相似的走势,一致性较高,符合地面光谱仪的测量。
分别在2006年9月、2007年9月和2010年9月进行了野外地面调查。按照森林类型、坡度等情况随机布设79个水平投影面积为500 m2的圆形样地,样地分布情况见图1。对样地内林木进行每木测量,测量参数包括胸径、树高、树种和郁闭度。将野外调查样地的类型分为混交林、阔叶林、针叶林3个类型组,混交林53块样地,阔叶林11块样地,针叶林15块样地。再将野外调查数据与汪清林业局的森林二类调查资料相结合。在覆盖研究区域内较均匀地选取阔叶林、针叶林、混交林、水体、其他各40个验证样本[13],二类调查资料选取的验证样本为每个900 m2。
图1 研究区样地分布图
猫群优化算法是一种模仿日常生活中个体猫和猫群的行为模式建立起来的群智能优化算法。猫群优化算法主要分跟踪和搜寻两种模式[14]。在搜索模式下,猫对自身位置的每一个副本应用变异算子产生新的位置(邻域),并将新产生的位置放在记忆池中,并进行适应度值计算[15]。在记忆池中,选择适应度最高的候选点,作为猫所要移动到的下一个位置点,从而利用全局最优的位置更新猫的当前位置。在迭代过程中,根据适应度函数计算并保留当前最优解。
搜索模式的参数。SMP表示记忆池,定义了每只猫所观察的范围,即搜寻记忆的大小,它用来存放猫所搜寻的位置点,猫将依据适应度值的大小从记忆池中选择一个最好的位置点;SRD表示所选维数的变化域,即个体上每个基因的改变范围;CDC表示维度的变化比率,区间为(0,1)。
当猫处于跟踪模式时,按自身速度进行移动,用速度——位移模型来移动每一位基因的值[14]。整个猫群经历过的最好位置,为目前搜索到的最优解xbest,即利用全局最优位置来更新猫的当前速度,再通过更新后的速度值改变猫的当前位置,不断趋向最优解方向[15]。具体活动行为描述如下:
式中:vk,d为改变每维的速度;xk,d(t)为t时刻第k只猫在第d维的位置。d为某只猫的基因位;k表示某一只猫;C为常量;r为0~1之间的随机数。Vi=表示每只猫都有一个速度,L为第i只猫的基因总长度数。
跟踪模式下,用适应度值最高的猫的第d位基因位,替换了当前猫的同一位置的基因位的值,从而更新了当前猫的速度,继而由公式(2)得出第k只猫的位置。
式中:xk,d(t)为t时刻第k只猫在第d维的位置;xk,d(t+1)为更新后第k只猫的第d个位置分量;vk,d(t+1)是更新后第k只猫的第d位基因的速度值。
2.4.1 结构设计
假设已知聚类数目,把一景高光谱数据分成非林地、针叶林、阔叶林、混交林等几种类别,每只猫都表示为一种可能的聚类结果。通过考察每个像元的光谱特征、空间特征或图像的形状、颜色特征等(如:DN值、方差、灰度梯度、反射率、饱和度、光谱吸收指数、熵、能量、分维、精度等),通过猫群优化算法搜索,将像元自动归类。设这些待分类像元集为X={Xα,α=1,2,…M},其中Xα代表第α个像元的总特征。Xα为N维模式,Xα的每个分量代表所需的某一种特征。假设Xα的第一个特征为灰度特征,用高光谱各个波段的灰度值表示,则当有n个波段时,前n个分量分别为该像元各个波段的灰度值。第n+1个特征位放入第2个特征,依次类推,原理如表1所示:
表1 待分类样本集结构
本文主要依据高光谱较为丰富的光谱特征来进行分类研究,采用灰度特征从待分类像元(样品集)X中找到一个划分C={C1…CK},使得总类内的离散度和达到最小。每只猫作为一个可行解,这些可行解组成猫群,也即解集。该解集可以通过聚类中心的集合作为每只猫所代表的解,即每只猫的位置是由C={C1…CK}中的K个聚类中心组成的。高光谱聚类中每个猫包含位置、速度和适应度。以某只猫的第1位、第j位、第k位为例,其编码结构如表2所示。
表2 个体猫的基本编码
Cj为第j个聚类的中心;每只猫的速度Vj为第j个聚类中心的速度值,它和Cj一样和Xi保持一致,都为一个N维向量。编程实现时,猫的位置编码结构表示为Cat(i).location[]=[C1,…,Cj…,CK];Cat(i).velocity[]=[V1,…,Vj,…,VK],即Vij=(Vi1,Vi2…,Vik),j=1,2,…,k;每只猫的速度V为每个聚类中心的速度共同作用形成。位置和速度之间的更新和交换只限于对应的基因位之间的交换,当每维的速度超过该维的最大速度,则该维的速度被限制为最大速度(Vmaxj)。
2.4.2 适应度值选取
猫的适应度值f为一个实数。用公式表示为:f=g/JC。式中:JC为总类内的离散度之和;g为一个较小的整数。因为目标是在待分类像元集X中找到一个划分C={C1…CK},使得总类内的离散度之和达到最小,即公式中的JC最小,f越大,该猫所代表的聚类划分的总类间离散度越小。
2.4.3 基于猫群算法的聚类实现
步骤1:设置参数。最大迭代次数为M,待分类别数K。依据经验,SMP=5;SRD=0.2(即每类的聚类中心的值每次的变化范围在0.2以内)表示只在每只猫的自身邻域内搜索;维数变化率CDC=1,在计算过程中,每一类的聚类中心都是动态变化的。
步骤2:初始化猫群。种群大小为m,每只猫为一个D维向量,聚类类型数目为K时,D=K。随机初始化每一维位置(xi,d),随机选取聚类中心,对于第i只猫,将每一个像元随机指定为某一类,作为最初的聚类划分。该聚类中心集合即作为猫i的位置编码;并为每一位置初始化速度vi,d为0;同时初始化一个随机分布矩阵,记录待分类像元的类别分布情况。
步骤3:计算每只猫的适应度值,反复进行,生成m只猫,即形成了整个猫群的搜索空间。
步骤4:引进分组率,将猫群分为搜寻模式和跟踪模式,该值可随机设定。猫经常处于懒散状态,因此分组率MR=0.02。
步骤5:为每只猫设置一个标志位。根据标志位信息确定属于哪种模式,如果是搜寻模式,标志位为0,进入第6步。否则,进入第7步。
步骤6:进入搜寻模式。复制每只猫自身SMP份,放入记忆池中。根据每只个体猫的突变率(CDC)和变化域(SRD)。根据程序Cat(i).location(j)=Cat(i).location(j)+Cat(i).location(j)(SDR*(r*2-1))更新其位置,Cat(i).location(j)表示第i只猫在第j个基因位的位置,r为随机数。对每一只猫,通过CDC的值,确定该猫的哪几个维度位置需要改变,可随机或设定。根据公式f=g/JC计算记忆池内所有个体猫副本的适应度(f)。通过排序或选择算子,从复制的SMP只猫(副本)中,选择适应度值最大的副本来代替当前猫的位置。此时,这个副本猫的每一个分量位置所代表的聚类中心对应代替原来猫的各个聚类中心。
步骤7进入跟踪模式。计算适应度的值,记录猫群的全局最优位置(C_gd.location[])。每只猫会根据更新速度更新自己的位置,找出全局最优个体,达到最优位置。程序为:
Cat(i).velocity(j)为第i只猫在第j个基因位的速度。计算t+1时刻第i只猫的位置,得到每只猫基因位(分量),得到每个类的聚类中心的解,不断迭代,猫逐步趋近最优。聚类中心的集合内的值是要求出的不同类别的聚类中心。
步骤8:对于每个待分类像元,根据猫的聚类中心编码,可以使用最邻近法则来确定待分像元的聚类划分,即对待分像元Xa。某待分像元到聚类中心的距离为d(Xa,Cj)=min d(Xa,Cl),l=1,2,3…K。若第j类的聚类中心Cj满足公式d(Xa,Cj)=min d(Xa,Cl),则Xa属于类j。比如:对第1只猫,第1个像元X1,初始时并不知道属于哪一类,将其随机指定为第3类阔叶林,即Cj中的j=3,其他像元也随机指定。对每个待分像元计算与每只猫里第j类聚类中心的距离,比如第1个像元实际和第1只猫的第2类针叶林的距离最小,那就说明它在第2次迭代中应被归为第2类针叶林。这时重新计算聚类中心,由于个体猫和整个猫群的位置和适应度值都在不断变化,每只猫的聚类中心集合C也在不断变化。每个待分像元不断计算和某只猫的每个聚类中心的距离,调整自己所属的类别。
步骤9:计算所有猫的适应度值,记录和保留适应度值最优的猫,即全局最优解。
步骤10:判断是否满足终止条件,如设置的最大迭代次数,或者得到的解符合某种阈值,则本次循环结束。否则,继续执行步骤4。
除上述步骤中的参数选取,其余参数设置如表3所示。
表3 基于猫群优化算法的高光谱聚类算法参数
本文利用剔除坏波段的92个波段,空间分辨率30 m的高光谱影像,通过MATLAB编程与ENVI实现仿真实验。根据针叶林、阔叶林、混交林等地面采集数据,分析森林类型的光谱特征,提取出几种森林类型的光谱曲线(见图2)。
图2 地物光谱曲线
利用猫群优化算法迭代20次,其迭代次数与适应度收敛情况如图3所示。
图3 迭代20次与适应度关系
选取迭代第2次与第20次的不同分类结果进行比较,当猫群算法经过第2次迭代后,结果如图4所示。
图4 第2次迭代后森林类型聚类结果
其混淆矩阵见表4。
表4 第2次迭代分类精度
从表4中可看出,经过2次迭代,200个样本中138个样本被正确聚类。混交林40个验证样本中,24个落在混交林分类区,7个样本落在针叶林分类区,5个被错分成阔叶林,3个落在水体分类区,1个被错分为非林地他,制图精度为60%。针叶林的40个样本中,有29个样本落在针叶林的分类区,6个落在阔叶林分类区,其余落在水体和非林地,制图精度为72.5%;阔叶林的40个样本中,有22个样本落在阔叶林的分类区,7个落在针叶林分类区,5个落在混交林的分类区,1个落在水体分类区,1个落在非林地,制图精度55%;水体的40个验证样本中,有33个分类正确,2个落在混交林分类区,2个落在针叶林分类区,制图精度为82.5%;非林地30个落在非林地分类区,4个落在混交林分类区,2个落在针叶林分类区,4个落在阔叶林分类区,制图精度为75%。总体精度=69%,Kappa系数=0.612 5。分类总体精度较低,主要因为猫群聚类算法尚未趋向收敛,聚类中心还接近于随机状态,高光谱各个像元到聚类中心距离尚未寻找到最小。
当猫群算法经过第20次迭代后结果如图5所示。其混淆矩阵如表5所示。
图5 第20次迭代后森林类型聚类结果
表5 第20次迭代分类精度
从表3中可看出,经过20次迭代,200个样本中167个样本被正确聚类,总体分类精度为83.5%,Kappa系数=0.793。混交林验证样本中,28个落在混交林分类区,4个样本落在针叶林分类区,5个被错分成阔叶林,2个落在水体分类区,1个被错分成非林地,制图精度为70%;针叶林的40个样本中,有33个样本落在针叶林的分类区,1个落在混交林的分类区,1个落在阔叶林分类区,其余落在水体和非林地区域,制图精度为82.5%;阔叶林的40个样本中有30个样本落在阔叶林的分类区,2个落在混交林的分类区,5个落在水体分类区,制图精度为75%;水体的40个验证样本中有39个分类正确,1个落在阔叶林分类区,制图精度为97.5%;非林地37个落在非林地分类区,3个落在混交林分类区,制图精度为92.5%。
混交林、针叶林、阔叶林、水体、非林地的制图精度在迭代20次后比迭代2次时,分别提高了10%、10%,20%、15%和17.5%,混交林和针叶林提高相对较小。总的来说,精度有了大幅提高。主要是因为算法经20次时收敛,适应度值为0.345。但也存在部分错分、漏分现象,原因是研究区部分山地地形复杂,且算法以像元光谱特征为主要分类依据,植被光谱存在同谱异物现象,较为突出的在波长0.5~0.525、0.65~0.7μm。混交林与阔叶林光谱特征相近,干扰了算法对类内距离的判断,导致混交林与阔叶林有混分现象。在波长0.6~0.63、0.7~0.74μm,混交林、阔叶林与针叶林光谱有重叠,较难区分。因此,对聚类总体精度造成了一定影响。但在21~24、28~31、39~47、55~60、62~77、85~112波段的67个波段内,森林类型的光谱相对有明显区别,所以混淆程度较低,能较好的区分开。虽然在可见光区域及近红外区域的波长0.7~0.92μm,针叶林的灰度值最低,但在波长0.48~0.7μm处光谱特征置于阔叶林和混交林之间,使得他们的灰度值差别不大,从而引起混交林的验证样本中有几个样本落在针叶林的分类区,降低了针叶林的分类精度。
为了减少由于光谱特征重叠造成的影响,得到更好的精度,本研究对猫群聚类算法中的部分参数在基于一般经验基础上,不断调整,进行实验,最终选择了表3的参数设置。首先,当种群m=5时,虽然收敛速度快,但分类精度降低,这是由于算法搜索空间变小,容易陷入局部僵局,在并未得到最优聚类划分时便已经收敛。当m=20时,发现所用时间比m=10时大幅增加,但相应的精度却未见明显增加。这说明在m=10时,算法在20次迭代时会收敛,可以满足高光谱森林类型聚类的问题复杂度,过多增大搜索空间,增加搜索次数,反而影响搜索的时间和效率。种群参数m对分类结果的影响如图6所示。
图6 参数m改变的精度对比
由图7可知,当SRD=0.4时,总体精度为80.5%,精度略有下降,且在不同迭代次数上,每次得到的精度中大部分比SRD=0.2时低,且算法收敛速度变慢,这主要是因为当SRD=0.4时搜寻模式下邻域范围扩大,在前期搜索较好,但波动较大,后期收敛速度变慢;当SRD=0.1时,分类精度为79%,搜寻模式邻域过于狭窄,容易陷入局部僵局。因此,根据比较验证了经验值0.2是较为合理的。
图7 参数S RD改变的分类精度对比
本论文在新型的群智能算法猫群优化算法的基本思想的基础上,对其加以改进,并设计出了编码规则及个体结构,将其应用到林业遥感的聚类方法中,将高光谱影像的聚类过程映射为猫群算法的寻优过程。利用不同地物的光谱特征,从HJ/1A遥感影像的115个波段中挑选出类内距离较小的集合,经过验证样本验证,对阔叶林、针叶林、混交林、非林地、水体的分类精度达到83.5%,比传统K-means方法和ISODATA方法分别高7%和6.5%。结果表明基于猫群算法的高光谱森林分类能较好地区分阔叶林、针叶林和混交林等样地类型。猫群优化算法拥有自适应、自组织、自学习性强的特点,虽然聚类不如监督分类精度更高,但使用猫群算法进行自学习聚类为森林类型识别提供了一种新思路。
[1]童庆禧,张兵,郑兰芬,等.高光谱遥感:原理、技术与应用[M].北京:高等教育出版社,2006.
[2]Fazakas Z,Nilsson M,Olsson H.Regional forest biomass and wood volume estimation using satellite data and ancillary data[J].Agricultural and forest Meteorology,1999,98/99:417-425.
[3]舒清态,唐守正.国际森林资源监测的现状与发展趋势[J].世界林业研究,2005,18(3):33-37.
[4]董连英.Hyperion影像森林植被分类方法与应用研究[D].长春:吉林大学,2013.
[5]吴玉明.森林对人类生活环境的影响和作用[J].现代农业科学,2008,15(12):66-67.
[6]谭炳香.高光谱遥感森林应用研究探讨[J].世界林业研究,2003,16(2):33-37.
[7]谭炳香,李增元,陈尔学,等.高光谱与多光谱遥感数据的森林类型识别[J].东北林业大学学报,2005,33(S):61-64.
[8]仲青青.基于粒子群优化的遥感图像聚类研究[D].杭州:浙江大学,2011.
[9]李俊明,邢艳秋,杨超,等.基于环境与灾害监测预报小卫星的树种识别[J].东北林业大学学报,2013,41(11):41-45,50
[10]韦玉春,汤国安,杨昕.遥感数字图像处理教程[M].北京:科学出版社,2007.
[11]高晓岚,汪小钦.多源遥感数据在植被识别和提取中的应用[J].资源科学,2008,30(1):153-158.
[12]熊文成,魏斌,孙中平,等.一种针对环境一号卫星A星高光谱与CCD数据融合的方法[J].遥感信息,2011(6):79-82.
[13]李俊明,邢艳秋,杨超.基于森林类型光谱特征的最佳波段选择研究:以HJ/1A高光谱影像为例[J].森林工程,2013,29(4):42-46.
[14]Shu-chuan Chu,Pei-Wei Tsai.Computational intelligence based on the behavior of cats[J].International Journal of Innovative Computing,Information and Control,2007,3(1):163-173.
[15]C.Chu-Chuan,T.PeiWei,P.Jeng-shyang.Cat Swarm Optimization[C]//China:In Proceedings of the 9th Pacific Rim International Conference on Artificial Intelligence,2006:854-858.