矿山恢复治理区植被物候与健康状况遥感监测

帅 爽，张 志，吕新彪，陈 思，马梓程，谢翠容

（1. 中国地质大学（武汉）地质调查研究院，武汉 430074；2. 中国地质大学（武汉）地球物理与空间信息学院，武汉 430074；3. 湖北省国土测绘院，武汉 430010）

0 引 言

矿产资源的开发为国民经济发展做出了巨大贡献，但粗放式开发导致的耕地损毁、土壤侵蚀、粉尘污染、噪音污染、水污染以及对生物多样性的影响逐步显现。矿山恢复治理工程对降低矿产资源开发活动对生态环境的影响至关重要[1]。随着国家生态文明建设战略的实施，“使矿山生态保护、修复工作贯穿整个矿山生命周期”是矿山地质环境修复与治理工作的新理念。

对因采矿活动而退化的土地进行植被重建和恢复工作具有挑战性，特别是矿山恢复治理区域土壤养分缺乏，不适当的护理或长期的搁置可能会永久性损害植被的健康[2]。恢复治理区植被健康状况是评价恢复治理工作效果、研究矿山退化土地环境变化的重要指标之一[3]。

传统的植被覆盖变化调查和植被健康调绘主要采用野外调查和实验室研究等常规方法，通常耗费大量人力物力，而且这些方法仅适合小区域的调查和评估[1]。前人利用多期次多光谱遥感数据结合各类型植被指数开展了大量矿山恢复治理区的监测和评估，取得了良好效果。Zhang等[4]验证了矿山恢复治理区植被覆盖度与植被指数的关系。Zenkov等[5]利用1984年至2016年多期次遥感数据研究了保加利亚TROJANOVO露天煤矿区的土地复垦动态变化过程。Xie等[6]提出了一种用于评估矿山损毁和恢复治理状况的新型MRAIs （Termed Mining and Restoration Assessment Indicators）指数，相较传统植被指数，MRAIs对矿山开发、修复引起的地类细微变化有更好显示。Yang等[7]基于NDVI时间序列，根据矿产资源开发和恢复治理特征，定义了8种NDVI时间序列变化模板，再利用一系列阈值和动态匹配方法，对研究区像元进行归类，获取矿山开发扰动时间和恢复治理状况。然而，这些研究仅针对矿区内土地利用类型变化和植被覆盖度变化进行遥感监测，忽略了恢复治理区植被健康状况的评估，另外在遥感图像上区分恢复治理区植被与正常植被也是急需解决的问题。

矿山恢复治理区土壤质地、结构、含水量、孔隙度和重金属含量等土壤质量指标的差异会引起植株生长稀薄、叶片短小、严重黄化[8-9]。由于反应植被健康程度的生理要素（如叶绿素含量、细胞结构、水分或氮含量）的变化可以通过反射光谱来监测[10]，部分研究者利用遥感技术评估恢复治理区植被健康程度和生长状况。Lu等[11]通过对TM3与TM4波段构成的二维特征空间分析，获得土壤线，将TM4波段坐标轴向土壤线旋转，消除土壤影响，进而评价特定时间植被的生长状态。胡秀娟等[12]利用WorldView-2卫星影像，基于归一化山地植被指数、氮素反射指数和黄光波段反射率3个因子，提出了一种新型植被健康指数（VHI），以快速、大面积地监测与评价水土流失区的植被健康状况。Karan等[1]利用2000年和2015年Landsat遥感数据和RVI、EVI和NDVI等植被指数监测露天煤矿恢复治理区植被变化，并利用遥感方法评估植被水分与植被健康程度的关系。然而，这些研究大多利用单一时相或多年相近月份的反射率光谱数据进行植被健康程度的评估，导致研究结果只能反映植被一个或几个生长阶段对土壤质量的响应，无法获得植被整个生长阶段对土壤质量的响应特征，可能导致结论的随机性。

除此之外，矿山恢复治理区土壤质量退化还会导致作物物候特征的改变[13]。作物生长期长度和生长速率等物候参数是作物生长模拟的重要指标[13-15]，也能反映土壤的质量状况。遥感植被指数的时间序列常被用于作物物候参数的提取和监测，该技术已经成功应用于作物类型识别、产量估算和生长状态监测[16]，成为物候学与其他学科联系的桥梁。但目前还鲜有通过研究矿山恢复治理区作物物候特征，评估恢复治理区植被健康状况及土壤质量，进而评价矿山恢复治理工程效果的研究。另外，对于单个矿山区域监测，MODIS数据难以满足作物识别的精度要求[17]，Landsat数据又难以满足植被物候信息提取所需的时间分辨率。而哨兵2号遥感数据具备小区域植被物候特征研究所需的时间分辨率和空间分辨率，同时其在植被探测特征的可见光至近红外波段具有光谱分辨率优势[18]，与Landsat、ASTER和其他多光谱卫星数据相比，哨兵2号数据在植被红边区域包含更多波段，使得最初为高光谱数据开发的植被红边位置（REP）计算技术可用于哨兵2号数据。

植被指数与作物产量密切相关，同时能反应作物的养分和生理信息差异。植被指数作为单产估算模型因子与作物产量建立统计相关模型是一种常用的作物产量估算方法[19]。GVI、NDVI、EVI、NDVI累积值[20]作为模型输入变量，建立玉米、小麦等作物产量统计相关模型，获得了较高的精度。当作物受到胁迫作用时，相应的氮、色素、酶等发生变化，可应用各种植被指数监测这些生理指标变化[21]。MCARI、NDVI、NDMI等植被数据常被用于作物色素、氮含量和冠层叶片含水量的检测[22]。通过分析恢复治理区和正常耕作区作物各类型植被指数值的差异，有利于获取不同区域作物养分和生理信息差异。

本研究使用2020-4-12-10-31日采集的13期哨兵2号遥感数据开展玄武岩矿山恢复治理区植被健康状况遥感诊断指标研究。主要内容包括：1）利用哨兵2号NDVI时间序列，基于NDVI时序重构、动态阈值、曲线曲率等方法，分别提取并对比恢复治理区和正常耕作区作物关键物候指标，分析恢复治理区作物与正常作物生长期内生长状况的差异；2）提取研究区多类型植被指数、植被覆盖度、植被红边位置等遥感指标参数，对比恢复治理区作物与正常作物各遥感指标参数差异，评价各遥感指标区分矿山恢复治理区作物与正常作物健康状况的能力。为评价矿山恢复治理工程效果提供依据。

1 研究区及数据源

1.1 研究区概况

研究区位于黑龙江省七台河市桃山区一处玄武岩采石场，地貌上为玄武岩丘陵，采石场于2012年开始开采，占地类型包括开采面、中转场地、固体废弃物和矿石堆等。玄武岩的开采活动造成了采石场及其周边耕地的挖损和压占。采石场北侧两处固体废弃物区域（图 1中A，B区域）于2019年进行了恢复治理，恢复治理类型为土地复垦，在对A、B区域固体废弃物堆进行削坡和表面平整的基础上，使用附近耕地土壤进行了堆垫，表层堆垫土壤厚度约为60 cm，堆垫土有机质、氮磷钾等土壤化学特征与周围耕地近似。但由于堆垫过程中，重型机械的碾压，导致表层客土被压实而土壤容重升高，进而影响植被扎根和水分与肥料的供给。A、B区域在土地复垦的基础上，种植了玉米，玉米的播种、施肥与周边农田同步。采石场周边土地利用类型主要为玉米旱地和田间道路。该区域属于温带季风气候，年平均气温2.4～3.9 ℃，年均降雨量约490 mm。

1.2 遥感数据及处理

哨兵2号遥感数据由欧空局（https://scihub.copernicus.eu/）免费提供。由于黑龙江省玉米作物种植期主要为4月至10月，收集了研究区2020年4月至10月云覆盖量小于5%的共13期（4月12日—10月31日）哨兵2号L1C级原始数据。利用Sen2cor软件对L1C级哨兵2号数据进行辐射定标和大气校正，获取10 m空间分辨率的地表反射率数据，用于研究区植被物候参数、植被指数、红边位置、植被覆盖度等信息的提取。不同时相的高空间分辨率遥感数据（高分2号、Google Earth数据）作为补充数据被用于对比分析，判断恢复治理时间，勾绘恢复治理区域范围。

1.3 野外调查

2020年8月24日开展野外调查工作，主要调查恢复治理区作物与正常作物的生长状态，调查中发现矿山恢复治理区与周边正常耕地玉米作物长势存在明显差异。矿山恢复治理区玉米植株生长稀薄，植株高度明显低于周边正常玉米作物，叶片短小且黄化现象明显，生长状态明显劣于周边正常耕作区域玉米作物（图2）。

2 研究方法

研究方法包含两部分，首先是建立研究区作物NDVI时序重构曲线并提取物候参数，在对哨兵2号影像进行大气校正获取反射率数据的基础上，计算研究区NDVI时间序列，对比、选取适当滤波方法，重构NDVI时序曲线，提取关键物候参数，对比分析整个生长期内恢复治理区作物与正常作物生长状况差异。其次是分别提取并对比恢复治理区作物与正常作物植被指数、植被覆盖度、红边位置等遥感指标值差异，评估各遥感指标区分恢复治理区作物与正常作物的能力，和诊断作物健康程度的可行性和敏感性。具体研究方法和流程如图3。

2.1 NDVI时序重构

在对研究区哨兵2号数据集进行大气校正获取反射率数据的基础上，计算获得NDVI时间序列。由于云层和阴影的影响，NDVI时序数据仍然存在较大噪声，因此在应用之前必须对其进行重构[23]。以往研究中对NDVI时序重构方法的优劣没有统一观点[10,23]，本研究中，对比常用滤波方法（Savitzky-Golay（S-G）滤波、Double Logistic（D-L）拟合与Asymmetric Gaussian（A-G）滤波）的去除噪声、保留植被生长特征的效果，对NDVI时序重构效果进行定性评价，同时使用统计分析对各滤波方法NDVI时序重构效果进行定量评价，选取最适合的去噪方法。S-G滤波器是一种基于滑动窗口的加权平均算法，它通过n阶拟合多项式来计算某一点附近固定数量的点的平滑值，达到降噪效果。在S-G滤波器中，多项式次数和窗口大小的设置至关重要。D-L拟合和A-G滤波方法是采用特定函数对时间序列中最大值和最小值附近的数据进行拟合，并将局部拟合函数组合成整体拟合函数，增强拟合函数的灵活性，使其符合NDVI时间序列数据的复杂行为。

玉米营养生长期生长速率表现出“慢—快—慢”的规律，其中拔节期生长速率最快[24]。从图4可以看出，这3种拟合方法均去除了大部分噪声，在NDVI达到峰值前，3种拟合曲线均呈“S”型，符合玉米的生长趋势，其中S-G拟合曲线与原始NDVI数据更接近，尤其在生长起始阶段更能反映原始NDVI曲线的变化趋势，而D-L和A-G的拟合曲线更平滑，这可能是因为S-G滤波器强调细节，而D-L拟合和A-G滤波方法重建了整体变化趋势[10]。3种滤波方法NDVI时序重构效果的统计结果如表1所示，D-L拟合方法的对原始数据的增大比其他两种方法更为明显。标准差、协方差和相关系数反映了重构曲线与原始曲线的接近程度（即保真能力），标准差、协方差越小，相关系数越大，曲线的保真能力越好[10]。3种方法中，S-G滤波结果的标准差和均方根误差均最小、相关系数最大，分别为0.191、0.024 2和0.993。与其他方法相比，S-G滤波器不仅能去除噪声，而且能更好地拟合NDVI原始数据。

表1 三种滤波方法NDVI时序重构效果的统计参数Table 1 Statistical parameters of NDVI time series reconstruction effect of three filtering methods

2.2 物候参数提取

玉米的生长大致分为营养生长（种子发芽到抽雄期结束）和生殖生长（抽穗开始到成熟）两个阶段[24]。在对NDVI时序重构的基础上，利用曲线曲率法和动态阈值法[25]分别提取恢复治理区和正常耕作区玉米的关键物候期（出苗期、拔节期、抽雄期、成熟期和物候期间隔时间）。

2.2.1 动态阈值法

利用Jönsson等[25]提出的动态阈值法，提取研究区玉米作物的出苗期和抽雄期。将作物生长曲线上升阶段，距离NDVI最小值为NDVI最大值与最小值差值10%的位置定义成作物营养生长的开始期，即出苗期；在作物生长曲线下降阶段，将距离NDVI最大值为NDVI最大值与最小值差值10%的位置定义为生殖生长的开始期，即抽雄期，如图5a所示。

2.2.2 曲线曲率法

利用李艳等[24]、Jönsson等[25]使用的曲线曲率法，通过计算NDVI时序重构曲线曲率极值确定玉米的拔节期和成熟期，NDVI时序重构曲线的曲率通过计算重构曲线的一阶导数获得，曲率绝对值最大值对应时刻为农作物生长（衰老）速率变化最大的时期，重构曲线上升曲线曲率最大值点为玉米营养生长的拔节期；重构曲线下降曲线曲率最小值点为玉米生殖生长的成熟期，如图5b。

2.2.3 物候期间隔时间

此次研究中，还选取了不同物候期之间的间隔天数作为物候指标，分析整个生长期内恢复治理区作物与正常作物生长期整体生长状况差异。具体指标计算如下

式中DOY为关键物候期的年积日天数，L1、L2和L3分别指4个关键物候期（出苗期、拔节期、抽雄期和成熟期）之间的间隔天数。

2.3 植被指数

野外调查中发现，相比正常玉米作物，恢复治理区玉米作物明显黄化。为对比各植被指数区分恢复治理区作物与正常作物、诊断矿山恢复治理区植被健康状况的能力，选取6种常用于绿色植被色素检测的经典植被指数，对比其对于区分恢复治区作物与正常耕作区作物和反应作物生长状况的能力。6种植被指数分别为归一化植被指数（NDVI）[26]，绿色归一化植被指数（GNDVI）[27]，用以估算叶绿素含量的多时相中分辨率成像光谱仪（MERIS）的陆地叶绿素指数（MTCI）[28]，与叶绿素a、叶绿素b和类胡萝卜素的单位面积浓度具有较强相关性的特征色素简单比值指数（PSSRA）[29]，用以估算植被冠层叶绿素含量的倒红边叶绿素指数（IRECI）[30]和反应叶绿素含量变化的改进的叶绿素吸收指数（MCARI）[31]，如表2。

表2 植被指数及其计算公式Table 2 Vegetable indices and their formulas

2.4 植被红边

绿色植被在680至780 nm之间反射光谱中最大斜率对应的波长位置被定义为植被的红边位置（REP，Red-edge Position）[32]。REP对于植被叶片与冠层绿叶素含量差异、环境条件差异对叶绿素含量估算的影响不敏感，可用于评估植被的叶绿素含量、物候期生长状况、健康状况以及重金属污染状况等[32]。通过计算REP可以判断植被冠层的生长期、病虫害状况以及重金属含量等。为了从不同类型的光谱数据中提取REP参数，前人提出了多种REP提取方法，其中最常用的方法包括最大一阶导数法、多项式拟合技术、逆高斯方法和线性外推技术等。

根据哨兵2号数据在红边区域的波段设置，本次研究中使用四点线性插补法[32]提取目标像元红边位置。该方法假设红边的反射率ρ曲线可以用一条位于780 nm和670 nm反射率中点附近的直线来近似，定义红边位置为

对于哨兵2号遥感数据，红边位置则定义为

2.5 植被覆盖度

前人研究证实了植被覆盖度和NDVI之间的显著线性相关关系[4]，所以通过建立二者之间的转换关系可以提取植覆盖度信息。利用像元二分法模型对植被覆盖度进行估算

式中FVC为植被覆盖度，NDVIsoil为全裸土（无植被覆盖）区域的归一化植被指数值，NDVIveg为全植被覆盖区域的归一化植被指数值。研究中通常选取累积像元数占比5%和95%处的NDVI值作为纯裸土（NDVIsoil）和纯植被（NDVIveg）的NDVI值[33]。

3 结果与分析

3.1 不同区域作物物候特征提取与对比

利用曲线曲率法和动态阈值法在恢复治理区（A区和B区）和正常耕作区（C区）内各像元重构的NDVI曲线上分别提取了出苗期、拔节期、抽雄期和成熟期等4个玉米作物关键物候指标。A、B、C三区域各像元出苗期、拔节期、抽雄期和成熟期均值如图6。可以看出，A区和B区玉米作物各关键物候期均晚于C区，C区平均出苗期为第147天，A区和B区平均出苗期分别为第159天和第152天；A区和B区平均拔节期分别为第195天和第192天，晚于C区的平均第183天；A、B、C三区域平均抽雄期比较接近，分别为第249天、第249天和第248天；A区和B区平均成熟期分别为第296天和第297天，晚于C区的平均第275天。显示出恢复治理区与正常耕作区玉米作物生长状况的差异。从3个区域各物候指标的统计参数（表3）来看，正常玉米种植区内部作物各物候指标标准差较小，物候指标较稳定，而恢复治理区内部作物各物候指标标准差较大，各像元间各物候指标变化明显，显示出复垦土地区域内作物生长状况存在差异。另外，三区域各物候指标值的振幅区间存在较大重叠，不利于区分恢复治理区和正常耕作区作物。

表3 训练区像元关键物候期的统计参数Table 3 The statistical parameters of the key phenological periods o the training areas

同时对比了训练区像元各关键物候期间间隔的时间长度，如图7。L1、L2和L3的统计参数如表4。整体而言，C区域拔节期与抽雄期间时间长度（L2）比A、B区域长，而A、B区域出苗期与拔节期间时间长度（L1）和抽雄期与成熟期间时间长度（L3）大于C区域。其中A、B区出苗期与拔节期间时间长度（L1）与C区相差较小，A、B区域平均出苗期与拔节期间时间长度（L1）分别比C区域分别长约1～2 d和3～4 d，这可能是恢复治理区土壤质地的差异对玉米作物生长前期影响较弱。随着作物逐渐生长，土壤质地的差异也逐步显现，表现为A、B区域平均L2比C区域分别短约10～11 d和7～8 d，抽雄期与成熟期间时间长度（L3）分别比C区域长20～21 d和21～22 d。相比之下，C区域L3参数与A、B区域相差更大。同时，统计参数（表4）也反映，3个区域内部像元L3标准差最小，且其L3取值振幅重叠度也最小，利于区分恢复治理区 和正常耕作区植被及其生长状况。

表4 训练区像元关键物候期间隔时间长度的统计参数Table 4 The standard deviation, maximum and minimum of the four phenological periods from the three training areas

3.2 不同区域作物植被指数对比

利用A、B、C三区域NDVI均值差异最大的第235天（2020年8月22日）哨兵2号数据计算A、B、C三区各像元NDVI、GNDVI、MTCI、PSSRA、IRECI、MCARI等植被指数，并对比恢复治理区（A区和B区）和正常耕作区（C区）像元各植被指数值的分布情况，如图8。

结果显示正常作物区像元NDVI、GNDVI、PSSRA和IRECI值振幅区间与恢复治理区无重叠，能较好区分恢复治理区与正常耕作区作物及其生长状态差异，而恢复治理区与正常耕作区像元MTCI值和MCARI值分布区间有较大重叠，难以区分三区域作物及其生长状态的差异。C区像元NDVI、GNDVI、PSSRA和IRECI值整体大于A区和B区，其中C区像元PSSRA值振幅区间与A区和B区差异最大，C区PSSRA值振幅下限（累积像元数占比5%处值）比A区和B区PSSRA值振幅上线（累积像元数占比95%处值）分别高12.61%和44.74%，C区NDVI、GNDVI和IRECI值振幅下限比A区振幅上线分别高1.94%，3.82%和5.81%，比B区振幅上线分别高7.20%、5.73%和42.59%。

3.3 不同区域作物红边位置对比

同样利用2020年8月22日哨兵2号数据计算A、B、C三区各像元植被红边位置指数（REIP），A、B、C三区各像元REIP值分布如图9。可见，相较正常健康作物（C区），恢复治理区作物（A、B）红边明显向短波方向移动，作物光谱的“红边蓝移”现象反映了恢复治理区作物与正常耕作区作物健康状况的差异。从REIP值的振幅区间上看，正常耕作区（C区）像元REIP值较集中，分布在728.58 nm至729.93 nm之间，而恢复治理区像元REIP值较分散，A区REIP值分布于714.98 nm至726.22 nm之间，B区REIP值分布于716.90 nm至728.38 nm之间，同时，恢复治理区像元REIP值振幅区间与正常作物区无重叠，表明植被红边位置指数（REIP）能较好区分恢复治理区、正常耕作区作物及其生长状态。C区REIP值振幅下限比A区和B区振幅上线分别高0.03%和0.32%。

3.4 不同区域作物植被覆盖度对比

通过对研究区NDVI图像像元值进行统计，确定累积像元数为5%和95%处的NDVI值，即纯裸土（NDVIsoil）和纯植被（NDVIveg）的NDVI值，分别为0.221和0.879，并依据公式6计算了A、B、C三区域的植被覆盖度（FVC），三区域像元FVC值分布如图10。结果显示正常耕作区（C区）像元FVC分布十分集中（0.937至0.968），而恢复治理区（A区、B区）像元FVC分布比较分散，A区像元FVC值分布在0.376至0.914，B区像元FVC值分布在0.537至0.851。C区像元FVC值整体大于A区和B区，同时C区FVC值振幅区间与A区和B区没有重叠，表明FVC能较好区分恢复治理区、正常耕作区作物及其生长状态。C区FVC值振幅下限比A区、B区FVC值振幅上限分别高2.53%和10.17%。

4 讨 论

4.1 不同区域作物物候参数差异分析

本研究利用哨兵2号影像时间序列和S-G滤波方法重构了研究区NDVI时间序列，并基于动态阈值、曲线曲率方法等方法提取并对比了矿山恢复治理区和正常耕作区玉米作物关键物候期和各物候期间时间间隔等关键物候参数，验证了矿山恢复治理区作物与正常耕作区作物物候特征的差异。相较正常作物，恢复治理区作物出苗期、拔节期、抽雄期、成熟期不同程度推迟，这可能是矿山恢复治理区土壤质地、含水量差异，甚至重金属胁迫减缓了玉米的生长速率，抑制了胚根的生长[34]。另外，恢复治理区与正常耕作区作物生长状态的差异是逐步显现的，表现为A、B、C区域像元各关键物候期间的时间间隔差异逐渐增大，这与前人对于重金属胁迫下水稻物候特征变化的研究结论类似[10,15]。各物候参数中，抽雄期与成熟期间时间长度（L3）对于诊断恢复治理区与正常耕作区作物生长状态效果最好，训练区内部像元L3取值较稳定（标准差最小，表4），且A区、B区和C区的L3取值振幅重叠度也最小。同时，此次研究验证了哨兵2号数据用于提取、分析作物物候信息的可行性，哨兵2号数据在可见光、近红外波段空间分辨率为10 m和20 m，重返周期5 d，一定程度上能够弥补Landsat数据（重返周期16 d）用于植被物候信息研究中时间分辨率不足以捕捉到快速发生的变化（如农田绿化）的问题[18]，以及MODIS、AVHRR和MERIS等数据用于非均质区域植被物候特征研究时的混合像元问题[17]。研究区周边佳木斯、集贤农业气象监测站点监测数据显示，该区域春玉米平均拔节期与成熟期分别为第180天至190天和第265至275天[35]，与此次提取的正常耕作区玉米作物拔节期（平均第183天）和成熟期（平均第275天）数据大致相符，验证了本文数据源和方法应用于玉米物候信息提取的可靠性。仍需要注意的是，由于部分时间段内（如第250天至275天，第275天至300天）影像云量过大无法使用，导致此次研究中影像时间序列间隔不一致，可能会影响物候信息提取的精度。下一步工作中，将尝试采用时空融合算法，利用MODIS等数据弥补部分期次影像云量过大的问题。

4.2 遥感植被指标敏感性分析

研究结果表明，REIP、NDVI、GNDVI、PSSRA、IRECI和FVC等遥感指标能较好区分矿山恢复治理区与正常耕作区作物及其生长状态差异。整体表现为，正常耕作区像元遥感指标值振幅区间相对收敛，而恢复治理区作物像元各遥感指标值振幅区间则比较发散（图8-10；表5），表明正常耕作区作物生长状态比较统一，而恢复治理区作物生长状态有一定差异，反应了恢复治理区内作物养分生理要素状况的差异。同时，正常耕作区像元REIP、NDVI、GNDVI、PSSRA、IRECI和FVC值整体大于恢复治理区。研究结果显示NDVI能较好地识别恢复治理区植被及其生长状况，与前人得出的NDVI能可靠评估植被健康的结论一致[36]。利用正常耕作区像元各遥感指标值振幅区间下限（累积像元数占比5%处值）高于恢复治理区像元各遥感指标值振幅区间上限（累积像元数占比95%处值）的百分比来定量对比各遥感指标对于区分恢复治理区作物及其生长状况的敏感程度（如图11），同时计算研究区像元各遥感指标值的标准差（表5）。结果显示C区像元PSSRA值振幅区间下限高于A区和B区PSSRA值振幅区间上限的百分比最高，分别为12.61%和44.72%，表明PSSRA指数对于区分恢复治理区和正常耕作区作物及其生长状况差异最敏感。各遥感指标中，恢复治理区像元的GNDVI值标准差最小（表5），表明恢复治理区像元的GNDVI值范围更集中，对于反应恢复治理区作物的整体生长状况更稳定；而正常耕作区像元的NDVI值标准差最小（表 5），表明NDVI对于反应正常耕作区作物的整体生长状况更稳定。

表5 训练区像元各遥感指标统计Table 5 Statistics of the remote sensing indexes from the three training areas

4.3 不同区域作物养分生理信息差异分析

不同类型区作物遥感物候参数、植被指数等指标差异反应了作物光谱特征的差异，而这种差异实质上是作物养分生理要素的差异造成的。研究结果表明，正常耕作区像元REIP、NDVI、GNDVI、PSSRA、IRECI和FVC值整体大于恢复治理区。前人研究显示NDVI、GNDVI、PSSRA等遥感指标与植被健康状况、叶片叶绿素含量呈显著正相关[27,29,36]，表明恢复治理区作物叶绿素含量小于正常耕作区作物，这与野外验证的情况相符（恢复治理区玉米作物叶片明显黄化）。

另外，研究表明GNDVI对于反应恢复治理区作物的整体生长状况更稳定，而PSSRA指数对于区分恢复治理区和正常耕作区作物及其生长状态差异最敏感。前人证实了PSSRA、GNDVI与玉米作物绿色叶面积指数的强相关性[37-38]，Nguy-Robertson等通过对多个玉米种植区实验研究，发现当玉米作物绿色叶面积指数大于2m2/m2时，PSSRA与玉米作物绿色叶面积指数的相关性强于DNVI、EVI、SAVI（Soil-Adjusted Vegetation Index）等11种植被指数[38]。王来刚等利用无人机多光谱数据进行玉米叶面积指数反演，发现GNDVI与玉米作物绿色叶面积指数的相关性强于DNVI、EVI、SAVI、OSAVI（Optimized Soil-Adjusted Vegetation Index）4种植被指数[37]。所以恢复治理区与正常耕作区玉米作物光谱特征、PSSRA、GNDVI值的差异可能与两区域玉米作物叶面积指数的差异有关。野外调查也证实了这一点，恢复治理区玉米作物株高和叶片大小明显小于正常耕作区，且叶片黄化明显，导致恢复治理区玉米叶面积指数与绿色叶面积指数小于正常耕作区玉米。

同时，恢复治理区表层堆垫土土壤容重的升高会影响植被扎根和水分的吸收。NDMI（Normalized Difference Moisture Index）与冠层的水分高度相关，与NDVI相比，NDMI可以更紧密地跟踪植物生物量和水分胁迫的变化，多项研究证实了NDMI在非均质环境中描绘植被叶片含水量的能力[22,39]。为验证恢复治理区与正常耕作区玉米叶片含水量的差异，对比了A，B，C区像元NDMI值的分布情况（图12），同时分析了各色素相关植被指数与NDMI的相关性（表 6，图 13）。图12显示，正常耕作区玉米作物NDMI值明显高于恢复治理区玉米作物，表明正常耕作区玉米作物叶片含水量大于恢复治理区玉米作物，这是由于恢复治理区表层堆垫土土壤容重的升高影响了玉米作物水分的吸收。另外表6与图13显示，各与色素相关的植被指数中，PSSRA与NDMI相关性最强（R2=0.778），表明相较其他植被指数，PSSRA能更好的反应作物叶片含水量的差异。

表6 植被指数与NDMI统计关系Table 6 Statistical correlation between vegetation indices and NDMI

综上，恢复治理区与正常耕作区玉米作物光谱特征、各遥感指标值的差异可能受两区域玉米作物叶绿素含量、叶面积指数和叶片含水量等生理信息差异的综合影响。另外矿山开采引起的水土污染和重金属富集也可能导致植被的光谱胁迫特征，下一步工作中进一步补充获取土壤质量、重金属元素数据，分析其对玉米作物光谱的影响。本次研究方法是通过检测因作物养分生理信息差异引起的光谱特征差异，来区分恢复治理区作物与健康作物，并评估作物的健康状况，可用于检测其他矿山恢复治理区域内植被生理信息差异引起的光谱特征差异。但由于其他矿山恢复治理区影响植被养分生理信息差异的具体类型（叶绿素类、氮素、水分、杂草、重金属等）可能不同，需要对各遥感诊断指标的敏感性进行再分析。

5 结 论

基于哨兵2号遥感数据，分别提取了恢复治理区和正常耕作区作物遥感物候指标、各类型植被指数、红边位置、植被覆盖度等遥感指标，并评估了各遥感指标对于区分正常耕作区与矿山恢复治理区作物，诊断作物生长状态和健康程度的可行性，主要取得了以下结论：

1）研究区矿山恢复治理区与正常耕作区作物生长状况的差异明显。恢复治理区玉米作物出苗期、拔节期、抽雄期、成熟期等关键物候期不同程度推迟。同时，恢复治理区与正常耕作区作物生长状态的差异随着作物生长阶段的推移逐步显现，表现为恢复治理区与正常耕作区作物像元各关键物候期间的间隔天数差异逐渐增大。

2）对比评估了红边位置指数（REIP）、植被覆盖度（FVC）、各类型植被指数等遥感指标对于区分正常耕作区与矿山恢复治理区作物，诊断作物健康程度的能力。REIP、NDVI、绿色归一化植被指数（GNDVI）、特征色素简单比值指数（PSSRA）、倒红边叶绿素指数（IRECI）和FVC等遥感指标能较好区分矿山恢复治理区与正常耕作区作物及其健康状况差异，其中恢复治理区像元的GNDVI值标准差最小，对于反应恢复治理区作物的整体生长状况更稳定。恢复治理区与正常耕作区像元的PSSRA值振幅区间差异最大，PSSRA指数对于区分恢复治理区和正常耕作区作物及其健康状况差异最敏感。

3）通过野外调查结合遥感分析，发现恢复治理区与正常耕作区玉米作物光谱特征、各遥感指标值的差异可能受两区域玉米作物叶绿素含量、叶面积指数和叶片含水量等生理信息差异的综合影响。

本次研究验证、评估了各遥感指标对于诊断矿山恢复治理区植被健康状况的可行性，为利用遥感技术进行矿山恢复治理效果评估提供了一种技术思路。下一步尝试将本文方法应用于区域性矿山恢复治理区植被健康状况的评估，进一步验证遥感诊断指标对其他矿山恢复治理区域植被健康状况评估的适用性。