森林植被恢复光谱特征分析——以幕府山矿区为例

邱洁，张亚丽，李明诗，2

(1.南京林业大学 林学院，南京 210037；2.南京林业大学 南方林业协同创新中心，南京 210037)

0 引言

森林是地球上陆地生态系统的主体，也是陆地上最复杂、最庞大、物种最多的多功能与多效益统一的生态系统[1]。作为一种可再生资源，在自然和人为因素的共同作用下，森林经历着发生、生长、干扰和死亡的过程，也即森林生态系统无时不处于消长交替的动态过程之中[2]。因此，研究森林植被的变化，对于掌握森林植被干扰和恢复规律、分析森林植被变化原因、预测森林变化趋势等方面有重要意义。

在森林变化研究方面，传统基于野外连续时间观测方式会耗费大量的人力和时间。现今，众多学者关注于基于遥感的森林覆盖变化模式特征提取，来推导森林植被恢复态势。如吴雪琼等[3]对森林覆盖变化遥感检测方法进行了研究，并指出制图尺度多元化、动态变化检测定量化和监测方法工具化3个方面的森林覆盖遥感检测的发展趋势。尽管归一化燃烧指数(normalized burn ratio，NBR)最初被用于监测森林火灾，很多学者仍借助NBR指数进行森林干扰制图等相关研究。如White等[4]就利用NBR指数对加拿大地区野生大火所引起的干扰和恢复进行了分析。之后，NBR指数被指出对干扰比较敏感[5]。李晶等[6]利用各种植被指数指示森林变化状况，运用NBR指数和综合森林指数(integrated forest z-score，IFZ)对矿区土地利用和植被覆盖类型进行了分类。沈文娟等[7]利用NBR指数对南方人工林地区进行了森林干扰和恢复特征分析。吕莹莹等[8]对南京市老山和紫金山地区运用基于IFZ指数的VCT算法进行了森林干扰和恢复的研究。李洛晞等[9]对基于MODIS数据的森林扰动指数进行了比较。但是，这些研究对基于Landsat数据发展的NBR和IFZ综合森林指数2种指数在刻画森林植被恢复性能差异方面少有涉及，特别是对于这2种指数对森林恢复趋势和状态刻画性能差异方面的研究仍然有待推进。

而Pickell等[10]也是运用了VCT算法得出干扰像元计算光谱恢复时间，提出了干扰像素恢复的平均年限的定义，并据此对北美寒带森林进行了4种指数的植被光谱恢复时间的研究。研究发现，NBR指数10 a间的平均光谱恢复时间(5.6 a)最大，TCG(缨帽变换绿度)指数得到的平均光谱恢复时间则最小(1.7 a)。森林光谱恢复年限的准确定义为研究森林植被恢复时间提供了定量的指标，能够从时间角度指示森林恢复状况。但是在中纬度地区的矿区植被恢复时间测定中，这些指数的性能还需要进一步验证。

从20世纪50年代起，南京幕府山地区出现矿石开采，导致森林植被大面积消失，森林生态系统服务功能逐渐减弱。幕府山地区森林破坏的问题引起了南京市政府的高度关注。1999—2008年，南京市政府对幕府山地区进行了10期的植被恢复工程，并于2003年完全关闭了白云石矿。文献[11-12]对幕府山地区森林植被恢复的生态效益以及植物多样性等方面进行了评价。蒋美琛等[13]利用像元二分模型对北京市重点矿山开采区进行了植被覆盖度的反演。至今，幕府山地区的森林植被恢复工程已经进行了近20 a，其具体恢复状况尚未有空间意义明确的评价结果。本文的目的是利用长时间序列Landsat数据，从植被恢复光谱特征分析的角度，验证NBR指数和IFZ指数在矿区植被恢复监测方面的有效性，并给后续的实践应用推荐恰当的方法。

1 数据与方法

1.1 研究区概况

本文以南京市幕府山地区为研究区(图1)。幕府山(30°54′N～ 32° 12′ N，116° 22′E～ 121° 54′E)是一座位于南京市鼓楼区和栖霞区的丘陵山脉，全长约5.5 km，宽约2 km，最高海拔199.3 m。幕府山地区属于亚热带季风气候，温暖湿润，降水充沛，年平均气温 14.6～16.4 ℃，最热月平均温度 28.1 ℃，最冷月平均温度为 -2.9 ℃；年平均降雨量为 800～1 000 mm；夏季多雨，冬春干旱，无霜期 237 d，每年 6月下旬到 7月中旬为梅雨季节[14]。自20世纪50年代起，幕府山地区最多时分布有 9个采矿场，总采矿区面积为 0.6 km2。1998年起，幕府山采矿场陆续关闭，矿区的植被恢复工作随即开展。幕府山地区自然植被以落叶常绿阔叶林为主，且现有的植被都是次生群落。主要树种有麻栎(quercusacutissimacarruth)、构树(broussonetiapapyrifera)、刺槐(robiniapseudoacacia)等。幕府山周边有幕府山街道、燕子矶街道等，周围城市和乡镇人口稠密，对幕府山有较大的影响。

图1 研究区位置示意图

1.2 数据来源

本文使用的遥感数据包括1987—2017年30期Landsat TM/ETM+/OLI 时间序列影像(具体描述信息见表1)。目标地的轨道号为120/038。数据在美国地质调查局(United States Geological Survey，USGS)地球资源观测和科学中心(Earth Resources Observation and Sscience Center，EROS)网站免费下载(https：//glovis.usgs.gov/)。影像获取日期的选择为5—9月(少量数据扩大到4月底或10月初)，此时为植物生长旺季，有利于研究森林植被变化，同时尽量选择云量少的数据，以避免分析过程中引入额外的不确定性。此外，本研究还从天地图上获取到部分年份的幕府山地区影像数据，主要用来进行目视解译并验证Landsat数据分析结果的准确性。

表1 本文使用的 Landsat TM/ETM+/OLI影像

1)预处理。本文采用陆地卫星系统干扰自适应处理系统(LEDAPS)和陆地卫星表面反射率代码(LaSRC)对表 1 中的 Landsat 图像进行标准化处理，完成辐射定标和大气校正，从而构建陆地卫星图像时间序列堆栈(LTSS)。LEDAPS 程序使用MODTRAN 太阳能输出模型，将定标影像转换为表观(top-of -atmosphere，TOA)反射率影像。此模型转换方法是校正定标影像的太阳方位、TM或ETM+带通、日地距离以及太阳辐照度。

同时，对于 TM 图像，辐射处理采用了最新发表的参数，计算其大气上界辐射亮度和反射率；对于 ETM+图像，大气上界辐射亮度和反射率是利用其头文件中提供的增益，以及偏差参数进行计算。为了弥补大气散射和吸收对于大气上界辐射亮度和反射率的影响，采用“6S”模型进行地表反射率计算，并依据标准的 MODIS反射率产品进行验证[15]。利用浓密植被(DDV)算法对表观反射率影像插值生成的气溶胶光学厚度(AOT)，以及通过相关资料获得的大气压、臭氧(O3)浓度和水汽值等数据，使用 “6S” 辐射传输模型，生成地表反射率产品。而 Landsat-8 地表反射率数据来自 Landsat LaSRC，它使用了MODIS的辅助气候数据和一种独特的辐射传输模型，利用沿海气溶胶带进行气溶胶反演测试。

2)VCT算法。在分析植被恢复特征前，使用植被变化追踪(VCT)算法提取幕府山地区的森林干扰斑块[16]。

3)植被恢复光谱指数。NBR已被成功用于评估火灾燃烧的严重程度、检测分类森林干扰和表征森林属性上[17]，其计算公式如式(1)所示。

(1)

式中：ρNIR代表近红外波段反射率；ρSWIR代表短波红外波段反射率(2.08～2.35 μm)。健康森林的NBR值较大，且更高的NBR值通常与增加的森林结构和冠层覆盖有关[18]。

文献[17，19]均利用干扰5 a后的森林恢复期研究森林恢复状况，所以下列光谱恢复指标均在干扰发生5 a后进行森林恢复状况的研究。

IFZ[20]常被用来进行森林干扰和恢复检测制图工作，其计算公式如式(3)所示。计算IFZ指数前需计算FZi指数(式(2))。

(2)

(3)

式中：bpi为影像上任一像元的亮度值；bi和SDi分别为森林像元的波段光谱值的均值和标准差；N为使用波段的数量。IFZ指数主要使用的波段是Landsat TM/ETM+的3、5、7波段，故先计算这3个波段的FZ值(FZ3、FZ5、FZ7)。当像元IFZ指数值越低并接近0.2，成为森林像元的可能性就越高。确定森林类别主要方法是黑体目标方法[21]。

4)光谱恢复时间。光谱恢复时间参照Pickell等[10]提出的干扰像素恢复的平均年限。对于VCT方法得出的干扰区域，计算像元扰动后光谱值恢复至扰动前2 a平均值的80%所需要的年限，即为干扰区域像元恢复年限。本文计算并分析了NBR和IFZ的恢复年限。

5)Mann-Kendall趋势检验法。Mann-Kendall趋势检验法[22]是一种非参数统计检验方法，由于不要求数据正态分布、不易受异常值的干扰、计算方便等优点被广泛应用。这种趋势检验法对于样本数为n的时间序列Xt(t=1，2，3，…，n)，统计量Z的计算如式(4)所示。

(4)

式中：

(5)

式中：xi和xj为第i年和第j年的数据值。对于给定的置信水平∝，若|Z|≥Z∝，则表示时间序列变化趋势显著，反之不显著。本文选择∝为0.05，所对应的审阅临界值Z∝为1.96。

2 结果与分析

2.1 光谱指标趋势

1)NBR指数光谱趋势。由于篇幅限制，本文以5 a为时间间隔，列出部分NBR图像(图2)，以表示30 a间NBR指数的变化趋势。

图2 1987—2017年幕府山地区NBR值图像

从图2可见，1987—2002年间，幕府山地区西南部NBR值较低，说明矿场的存在给幕府山地区森林植被带来了严重破坏。1998年幕府山地区开始进行生态治理，2007年后幕府山地区NBR值有明显增加，说明此地区特别是矿区的森林植被恢复效果显著。此外，从2002年开始，幕府山边缘地区的NBR显示出较低水平，可能与周围城市扩张有关。

2)IFZ指数光谱趋势。Huang等[16]的研究指出，不同于其他土地覆盖类型(水除外)，落叶和针叶树森林IFZ值一般低于3。由于本文研究区幕府山地区并没有水体的分布，所以不考虑水体覆盖类型的影响。除此之外，本文对IFZ值进行了归一化处理，故以大于0.3作为判定为非森林的依据。IFZ值大于0.3时可以判定为非森林；小于0.3时则为森林。图3列出了以5 a为时间间隔，部分年份幕府山地区IFZ指数分类图像。总体上，1987—1997年间，森林植被覆盖呈现下降趋势，同时西南部在2002年之前的IFZ图像上一直显示为非森林。1997年后，森林植被逐渐恢复，这与前面其他指标导出的结果相一致。但是在IFZ图像中，幕府山边缘地区均为非森林部分，无法显示出城市化进程的影响，这可能与IFZ指数和NBR指数对于森林识别的敏感度不同有关。

图3 1987—2012年幕府山地区IFZ指数图像

2.2 光谱恢复时间

通过VCT算法得到了幕府山地区各年份像元干扰图，并且按照像元计算每年干扰区域干扰发生后5 a内光谱恢复时间，得到了不同光谱指数各年份5 a内光谱恢复率(图4)以及各年份5 a内恢复像元的光谱恢复时间(图5)。干扰发生后5 a内恢复像元的平均恢复率NBR指数为84.84%，低于IFZ指数的98.53%。NBR指数的5 a内像元恢复率均小于同年IFZ指数的像元恢复率，且IFZ指数像元恢复率全都超过90%。由此也佐证了IFZ指数能够监测更多森林恢复。同时各年份5 a内恢复像元平均恢复年限NBR指数为1.54 a，高于IFZ指数的1.04 a。1999年NBR指数恢复时间达到了最大(3.09 a)，其余年份像元平均恢复年限则与均值相差不大。这可能是由于南京市于1998年开始实施恢复政策，此时的植被干扰状况最为严重。

图4 干扰区域像元5 a内像元恢复率

图5 每年干扰区域恢复像元平均恢复年限

2.3 Mann-Kendall趋势检验

为了更好地监测幕府山地区30 a间的森林植被变化趋势，对每年的NBR指数和IFZ指数构成的时间序列逐个像元进行Mann-Kendall趋势分析和显著性检验，得到不同光谱指数的30 a趋势图(图6)和各部分趋势所占面积比例表(表2)。在NBR指数的趋势分析图中，幕府山矿区集中区域呈现显著上升趋势，幕府山外部地区(靠近城市)则呈现出不同程度的下降趋势，且边缘地区显著下降趋势较集中。IFZ指数得出的结果无论是面积比重还是分布上都与NBR具有高度相似性，不同的是IFZ指数检测得到趋势结果对于不显著的趋势更加敏感，不显著上升(26.85%)和不显著下降(19.86%)均大于NBR指数的比例，同时对于显著下降的地区，IFZ指数(19.47%)也表现得更加明显。

图6 NBR和IFZ时间序列Mann-Kendall趋势分析图

表2 NBR指数和IFZ指数Mann-Kendall趋势检验各部分趋势面积比 %

值得注意的是，Mann-Kendall趋势检验包含了所有的像元，包括持续森林或者持续非森林部分。对于这些地物持续不变像元，发现大多数像元均具有不显著的趋势，包括上升和下降趋势，这样的趋势特性与持续像元的特性相符合。

3 讨论

3.1 长时间序列影像

本研究利用长时间序列影像指示森林变化，存在以下优点：充分利用了30 a间的Landsat影像数据，避免了使用少量影像对比存在的影像代表性不足等问题，更能体现森林变化中的微小变化。沈文娟等[7]也指出，仅凭2幅或3幅遥感影像难以对南方人工林地区频繁的采伐更新活动进行很好的描述。

3.2 光谱恢复年限

在Pickell等[10]的研究中，NBR指数干扰发生后5 a内平均恢复了78%的干扰像元，NBR指数10 a间的平均光谱恢复时间最大(5.6 a)，TCG(缨帽变换绿度)指数得到的平均光谱恢复时间则最小(1.7 a)，在寒冷湿润的地区NBR指数光谱恢复年限与寒冷干燥地区有差别。发生干扰像元的比率不同可能与森林类别、土壤环境不同有关。本文对亚热带气候的矿区植被进行了光谱恢复时间的研究，得到的结果与上述结果相比光谱恢复时间较短，主要可能与气候、立地条件、树种等因素的影响有关。此外，也可能与pickell研究范围是干扰发生10 a间的光谱恢复有关。Brandt等[23]也在研究中指出，气候因子和干扰因子是2个塑造加拿大北方森林景观的主要因素。

3.3 可靠性分析

本研究通过NBR指数和IFZ指数这2组独立的指数分析，得出了高度相似的结果，证明了本研究中光谱恢复指标的可靠性。李林等[21]在对幕府山地区植物区系的研究中指出，目前幕府山地区的植被群落结构正走向复杂，此地区的森林植被恢复状况良好，这与本研究中幕府山地区森林植被恢复状况良好的结论一致。此外，潘叶等[11]对南京幕府山地区的生态修复工程进行了评价，与本研究的几个恢复指标指示的森林恢复状况整体趋势相同。

通过Mann-Kendall趋势检验可以明显看出，幕府山边缘地区森林植被显著下降，这与周围城市扩张有很大程度的关系。吕莹莹等[24]通过对1992—2011年间南京市幕府山地区DI指标(disturbance index)的变化趋势进行分析，发现幕府山东南部矿区集中区干扰逐渐减弱，而幕府山周边地区干扰程度逐渐增加，这也佐证了本研究结果的可靠性。

4 结束语

幕府山地区自1949年以来一直进行大面积采矿，最直接的表现就是在幕府山西部面积为79.19 hm2的废弃白云石矿矿坑，使得森林植被面积大面积减少[25]。自1998年开始，南京市市政府开始对幕府山地区进行植被恢复，实施了整改关停等措施。

郑玲等[12]在文章中指出，由于幕府山地区长期的生态恢复，植被群落处于重建和恢复阶段，初期恢复措施已经不适用于现有林分。因此，对此地区的生态恢复进行监测，从而调整恢复措施显得尤为重要。徐志超等[26]指出，由于矿区采矿活动，水土流失严重，而且土壤厚度和土壤肥力各不相同，加大了矿区生态植被恢复的难度。因此，针对此地区不同立地条件森林恢复状况进行分析，可为因地制宜地实施合理措施提供有效意见。王军等[27]对幕府山地区的群落物种多样性进行了研究，发现相对于灌木和乔本植物，草本植物较为丰富，正处于群落演替的初级阶段。本文研究了幕府山地区的光谱恢复时间，针对不同地力状况，选择生长速率合适的树种，从而做到适地适树，提高恢复效率。

本研究使用NBR和IFZ 2类恢复指数对幕府山地区的森林植被恢复进行了不同维度的分析。在对幕府山地区实地调查资料进行分析后发现，IFZ指数在幕府山地区能够监测到更多的森林恢复，光谱恢复时间也相应较短。从监测幕府山地区植被恢复状况的要求来说，对森林恢复更为敏感的IFZ指数相较于NBR指数更适合于现今的幕府山森林植被恢复状况的监测。同时，更短的光谱恢复时间与幕府山现处的群落演替初级阶段的状况更加相符，也有助于迅速调节生态恢复状况和实施更有效的生态修复手段。但是NBR指数对幕府山边缘地区森林植被的干扰比IFZ指数更为强烈，因此对监测边缘地区扩张导致森林干扰方面，NBR指数表现更优。本文对不同立地条件的森林植被恢复生态效益进行了评价，同时也为此类地区植被恢复选择合适的评价光谱指数提供了参考意见，进一步健全了矿区生态恢复效果监测方法体系。