基于混合像元分解的香格里拉市高山松空间分布变化研究

刘蓉姣，张加龙，陈培高

(西南林业大学 林学院，云南 昆明 650224)

香格里拉市地处世界自然遗产“三江并流”的核心地带，境内地形地貌复杂,森林资源丰富，但生态系统较为脆弱[1]。目前对其森林资源、生态环境和多树种分布等方面已有研究，但对该区域某单一树种分布变化的研究较少[2]。本研究围绕香格里拉市优势树种高山松空间分布变化展开。高山松(Pinusdensata)是我国的特有树种，源于第三纪云南松(P.yunnanensis)与油松(P.tabuliformis)的天然杂交种，耐寒、旱、贫瘠，广泛分布于西藏东南部，川西部，云南西北部和青海南部2 700～4 200 m的高山地带，多成纯林分布[3-4]有较高的经济价值和生态价值。

由于受传感器分辨率的限制和地面地形气候多样性等因素的影响，遥感影像中，像元通常都是几种地物的混合体，很少由单个均匀的地表覆盖类组成[5]。使用高分辨率影像进行目视解译所花费的代价太大，无法实现大范围的分类[6]。而对于中低分辨率影像来说，混合像元的存在很大程度上降低了分类精度。如果每个混合像元都可以被分解并且其覆盖类型组分(通常称为端元组分)占像元的百分比(丰度)也能够被求出，分类将会更精确。多年来，国内外学者利用混合像元分解技术提高遥感分类的准确性，并取得了有效进展。杨超[7]利用混合像元分解技术进行土地利用/覆被类型细分类。王铭[8]采用GOCI遥感影像提取时令水体及变化检测。H.Luy[9]利用混合像元分解对卫星时间序列数据准确地估算了泛太平洋25个城市的植被。M.Dawelbaitetal[10]利用混合像元分解技术，监测North Kurdufan State (Sudan)1987-2008年的沙漠化过程，并在时间序列下分析了研究区内沙漠与植被覆盖的空间特征及时间变化特征。陈虹兵等[11]利用像元二分模型法和线性光谱混合模型法对提取植被覆盖度，对比分析2种模型的适用性，最后利用线性光谱混合模型生成了延庆区植被覆盖分布图,为县域尺度应用Landsat 8 OLI影像估测植被覆盖度提供参考。博文杰等[12]利用线性光谱混合模型提取莆田市常太镇枇杷树信息，绘出枇杷树的(Eriobotryajaponica)丰度图，据此估计出枇杷的种植面积。

本研究运用混合像元分解结合决策树分类，提取香格里拉市高山松空间分布情况，分析出高山松面积的动态变化和高山松林地的发展趋势，对香格里拉市森林资源的及时保护和进行宏观决策提供参考依据。同时，建立的高山松分类指数对高山松分类研究具有一定参考价值。

1 材料与方法

1.1 数据源

香格里拉市是迪庆藏族自治州下辖市之一，位于云南省西北部的滇、川、藏大三角区域，地理位置99°20′-100°19′E、26°52′-28°52′N[13]。香格里拉市地形呈现西北高、东南低，最高点海拔5 545 m，最低点海拔1 503 m，海拔高差4 042 m，平均海拔3 459 m，地貌形态集山地、高原、盆地、河谷为一体[14]，市境内的土地总面积为1 141 739 hm2，林业用地面积950 911.7 hm2，占83.3%，森林覆盖率为74.99%，据2016年全国森林资源二类调查统计，云冷杉、高山松、高山栎、云南松为主要优势树种，占全市森林面积的90.8%[15](图1)。

图1 香格里拉市位置Fig.1 The location of Shangri-La

收集的数据有遥感数据和外业数据。遥感数据分别为香格里拉市1997、2002、2007、2012年的Landsat 5 TM数据和2017年Landsat 8 OLI影像。其中，2012年Landsat TM退役，采用时间相近的2011年的影像代替，形成每5 a间隔的时间序列影像集(表1)，研究区所在影像的轨道号为132/040、132/041、131/041，三景影像覆盖了整个研究区，按照研究区边界进行了裁剪，不存在拼图问题[16]。同时收集了2016年的全国森林资源二类调查数据，香格里拉市DEM数据。外业数据收集于2018年4月，对香格里拉市42块30 m×30 m的样地(168个样地点)进行坐标定位、周围地物情况记录、4个方位拍照和样地内树种进行树高胸径测量。

表1 研究区Landsat影像Table 1 Landsat images in the research area

1.2 研究方法

以Landsat遥感影像为数据源，首先，对5期影像进行辐射定标、大气校正、几何校正和地形校正的数据预处理；其次，通过最小噪声分离、纯净像元指数计算和N维可视化，获取植被、裸地和阴影的纯净像元；然后利用线性波谱分离和最小能量变化进行丰度提取；接下来，利用分量构建指数，获取植被指数及DEM数据；最后，利用获取的分类因子构建多元数据进行决策树分类并验证分类结果(图2)。

图2 技术路线Fig.2 Flowchart of methodology

1.2.1 数据预处理 利用Envi对Landsat遥感影像进行了辐射定标，大气校正，几何校正，地形校和镶嵌裁剪的预处理，最终获取研究所需影像。其中，几何校正利用一幅精纠正过的香格里拉市SPOT 5影像为基准影像，校正中每一景影像选择了30个控制点，平均RMS误差控制在30 m内。香格里拉市地形和地势落差较大，采用地形校正的方法来补偿由于不规则的地形起伏而造成的地物亮度的变化[17]，校正过程采用J.Nichol et al[18]的坡度匹配模型。一副研究区影像涉及3景数据，在完成各景的地形校正后，利用研究区的矢量边界进行切割，将裁剪后的3景图像进行镶嵌后形成研究区整体图像。外业数据处理，根据外业调查数据所记录的东西南北情况及所拍摄照片，运用ArcGIS获取实地高山松光谱曲线。

1.2.2 基于混合像元分解的决策树分类

1.2.2.1 高山松、土壤、阴影分量提取 样地调查中高山松多成纯林分布(图3)，森林中郁闭度不是很高的地方存在一些裸露的土壤。在日光斜照射下，具有一定高度的地物都会形成阴影,树木也不例外[19]，而且香格里拉市地势落差较大，阴影在影像中的存在就更明显了。所以在只包含某一个树种的影像中一个像元内就可能同时包含植被、土壤、阴影多个地物的光谱信息。

图3 外业调查Fig.3 Field survey

首先，对遥感影像进行最小噪声分离(minimum noise fraction)判断图像数据维数，分离数据中的噪声[20]，降噪后前5个波段已经集中了数据的85%以上的信息，对前5个波段进行PPI计算获取纯净像元。其次，通过N维可视化结合图像端元和参考端元选取高山松、裸地、阴影的端元。对照影像通过不断调整获取裸地和阴影的纯净端元，利用外业数据绘制的实地高山松波谱曲线和Landsat影像不断调整，获取高山松纯净端元。最后，在N维可视化中勾选的高山松绘制波谱曲线并与外业调查数据获取的高山松进行对比，得到高山松与实地外业数据获取的高山松波谱曲线对比情况(图4)。利用IDL编写的完全约束最小二乘算法[21]混合像元分解扩展工具对地物端元进行混合像元分解，避免了提取的丰度中出现负值的现象。

图4 高山松光谱曲线对比Fig.4 Comparison of spectral curves of Pinus densata from 1997 to 2017

1.2.2.2 决策树提取高山松空间分布 基于该特征运用丰度中3个分量构建植被-土壤指数、植被-阴影指数、植被-土壤-阴影指数和归一化多分量指数4种新的指数[22-23]，新指数用于决策树中来区别植被里的高山松、云南松、云冷杉和其他植被的关键数据。利用ENVI 5.3进行波段计算构建指数(表2)。

表2 新指数计算公式Table 2 New index calculation formulae

运用4种指数分离高山松与其他树种，决策最优指数进行分类。4种指数总体精度都不差，但相比较归一化多分量指数精度最高，归一化多分量指数总体精度都高于87%，其他3种指数精度在84%～87%，结合归一化多分量指数经过多次试验通过制定规则见图5，运用该规则进行分类，将影像分为高山松、云南松、云冷杉、其他林、非林地5类。香格里拉市云南松和云冷杉是占比较大的优势树种，占香格里拉市林地总面积的31.85%[24]，该规则通过海拔和其他树种的归一化指数先分类出了香格里拉市两大优势树种，这能更准确得到高山松空间分布情况。

图5 决策树分类示意Fig.5 Knowledge expert decision tree

1.2.2.3 精度评价 研究采用混淆矩阵[25]进行分类结果精度验证。通过统计树种的面积与香格里拉市面积并计算其比例[26]，运用公式(1)计算出验证点个数为507个。随机筛选出507个验证点，其中高山松的验证点76个，云冷杉105个，云南松42个，其他树种208个,非植被76个。对本研究分类结果进行了验证。

(1)

式中，n为验证点个数，B为卡方分布临界值，由(1-α)/k决定，其中α为置信度，k为分类个数。πi为第i类地物的比例，bi为每个像元的分类绝对误差(一般取5%)。

(2)

(3)

2 结果与分析

对分类结果进行了验证，精度评价结果见表3。各年分类结果总体精度分别为69.42%、76.73%、81.07%、78.90%和76.53%。其中高山松的分类精度较好，其次是云冷杉和云南松，非林地与其他林错分的情况相对多一些。

表3 精度评价结果Table 3 Accuracy evaluation results

2.1 高山松空间分布变化

运用基于混合像元分解的决策树分类，得到分类结果(图6)，红色区域为其他树种空间分布情况，其他树种包括了如车桑(Dodonaeaviscosa)、杜鹃(Rhododendronsimsii)、华山松(Pinusarmandii)、箭竹(Fargesiaspathacea)其他等多种树种；绿色区域为高山松空间分布情况；蓝色区域为云南松空间分布情况；电气石绿色区域为云冷杉空间分布情况。其他树种与云冷杉空间分布情况有重叠部分，而且两者的空间分布面积最大。

图6 分类结果Fig.6 Classification results

对高山松总面积进行了统计(图7)。1997年高山松的总体覆盖面积为1 608.67 km2，占香格里拉市植被面积的17.08%。2002年高山松总面积为1 393.07km2，占香格里拉市植被面积的14.79%。2007年高山松总面积为1 358.66 km2，占香格里拉市植被面积的14.42%。2012年高山松总面积为1 480.46 km2，占香格里拉市植被面积的15.72%。2017年高山松总面积为1 541 km2，占香格里拉市植被面积的16.36%。2002年高山松覆盖面积比1997年减少了13.40%，2002年比2007年减少了2.47%，2012年高山松覆盖面积比2007年增加了8.96%，2017年比2012年增加了4.06%。不同年份香格里拉市高山松分布状况不同，不同年份受降水、气温、人为行为、自然更替等不同因素的影响都可能影响香格里拉市高山松覆盖面积。

图7 各年香格里拉市高山松面积统计Fig.7 P.densata areas of the study years in Shangri-La city

应用ArcGIS中叠加分析的交集取反分析(symmetrical difference)[27]，将相邻年份高山松空间分布进行交集取反分析，获得相邻年份高山松增量或减量的空间分布情况，通过分析图8可知，1997-2002年5 a间高山松急剧减少，2002-2007年高山松总面积减少量低于1997-2002年。2012-2017年高山松面积增量小于2007-2012年的增量。

图8 高山松分布变化Fig.8 Distribution change maps of P.densata

2.2 香格里拉市高山松随高程的分布变化

结合DEM数据得知，每年高山松空间分布的海拔范围情况分别为：1997年高山松分布的高程范围是1 827.3～4 516.5 m，2002年高山松分布海拔范围1 995.2～4 483.8 m，2007年高山松分布海拔范围2 094.3～4 371.1 m，2012年高山松分布海拔范围2 100.4～4 360.5 m，2017年高山松分布海拔范围2 092.6～4 250 m。对同一海拔区间的不同年份高山松覆盖面积占当年高山松总面积的百分比进行了简单统计对比。海拔在1 800～2 300 m的不同年份高山松覆盖率分别为：0.27%、0.11%、0.02%、0.04% 和0.02%。2 300～2 800区间每年高山松森林覆盖情况分别是：4.28%、3.31%、2.38%、2.55%和2.59%。2 800～3 300 m区间高山松面积占当年高山松总面积的比率分别为：30.93%、29.95%、30.02%、35.61%和36.81%。3 300～3 800 m海拔区间内的高山松覆盖率是最高的，该区间高山松面积占当年高山松总面积的比率分别为：56.14%、58.00%、57.60%、58.45%和58.36%。3 800～4 300 m海拔区间的高山松面积与当年高山松总体面积之间的比值分别为8.30%、8.54%、3.23%、3.36%和2.22%。4 300～4 800 m海拔区间的高山松面积与当年高山松总体面积之间的比值分别为0.08%、0.09%、0.004%、0.002%和 0.003%。该区间是每年中高山松分布最少的区间。

每年高山松分布区间最低海拔为1 827.3 m，最高海拔为4 483.8 m。利用ArcGIS软件的空间分析结合重分类将高山松分布区域海拔高度划分6个高程区，分别为1 800～2 300、2 300～2 800、2 800～3 300、3 300～3 800、3 800～4 300、4 300～4 800 m。由图9可知，不同海拔区间对比高山松分布情况为：3 300～3 800 m>2 800～3 300 m>3 800～4 300 m>1 800～2 300 m>4 300～4 800 m。

图9 各年各海拔区间高山松空间分布Fig.9 Spatial distribution of P.densata at different altitudes in each year

3 结论与讨论

以香格里拉市高山松为研究对象，利用1997年至2017年的Landsat遥感影像，通过混合像元分解结合决策树分类，获取了1997年至2017年近20 a高山松空间分布，分析其长时间序列变化情况。1997、2002、2007、2012年和2017年分类总体精度分别为：69.42%、76.73%、81.07%、78.90%和76.53%。

通过多次试验将分类的阈值从一年推广到多年进行分类。其中，NDVI植被指数区分林地非林地阈值为0.45，高山松归一化多分量指数阈值为0.333，云南松归一化多分量指数阈值为0.208，区分云冷杉的归一化多分量指数值为0.362，高程为3 700～4 350 m。

香格里拉市高山松覆盖率在1997-2002年急剧减少，2002年高山松覆盖面积比1997年减少了13.40%，2007年高山松覆盖面积比2002年减少了2.47%，从2007-2012年和2017年香格里拉市高山松覆盖率都在逐年增加，增幅没有很大，2012年高山松覆盖面积比2007年增加了8.96%，2017年比2012年高山松覆盖面积增加了4.06%。其中1997-2002年变化最显著。

结合高程分析了高山松的空间分布情况。1997、2002、2007、2012年和2017年香格里拉市高山松主要分布在海拔高度为2 800～3 800 m的地区，高山松生长分布最低海拔1 800 m，最高海拔为4 400 m。不同海拔区间对比高山松分布情况为：3 300～3 800 m>2 800～3 300 m>3 800～4 300>1 800～2 300 m>4 300～4 800 m。其中高山松最适宜生存海拔为3 300～3 800 m。总之，不同年份香格里拉市高山松分布状况不同，不同年份的降水、气温、人为行为、自然更替等不同因素都会影响香格里拉市高山松覆盖面积。

随着遥感技术的不断发展，对树种分类的探究也越来越被需要。研究树种分布空间变化首先要能准确进行树种分布的提取，这就取决于分类结果，然而影响分类精度的关键因素之一是分类方法的选取。由于本研究区面积过大、地形差异大、植被覆盖率高、生物多样性丰富以及遥感影像分辨率低等因素的影响，运用光谱信息、坡度、坡向、纹理等特征难以区分出不同的树种，因此采用混合像元分解结合决策树进行分类，在解决多光谱数据由于大气散射和空间分辨率低等问题造成的光谱混合现象后，结合前人经验运用决策树进行分类，获取了香格里拉市高山松空间分布情况，为后序分布变化提高可靠数据。今后希望能用分辨率较高的数据对林分区别较大的区域进行混合像元结合决策树准确分类更多树种。