1. 引言

近年来，植被恢复逐渐成为矿区生态重建以及矿区周围环境改善的一种重要手段[1]。传统的现场实地调查方法往往耗时耗力，并且很难对矿区周围生态环境进行实时动态监测。因此，利用遥感技术研究矿区周围植被覆盖的变化是矿区生态重建的常用方法[2]。

归一化植被指数又被称作标准化植被指数，它是由近红外和红外反射率的比率进行归一化得来的，反映为植被在电磁波谱中红光与近红外的吸收与反射特征，因此它是植物生长状态以及植被空间分布密度的最佳指示因子，与植被分布密度呈线性相关关系[3]。大量实验数据表明，NDVI数据对土壤背景的变化较为敏感，可以综合反映单位像元内的植被类型、覆盖形态、生长状况，其大小取决于植被覆盖度和叶绿素含量等生物物理参数，此外NDVI数据对植被覆盖度的检测幅度较宽，有较好的时间和空间适应性[4]。

由于在除去外界条件影响的前提下，植被覆盖随着时间的变化是连续的，因此NDVI时序变化曲线在理论上应当是一条连续且光滑的曲线，且能够反映由于季节变化或病虫害等影响因素引起的曲线波动[5]。然而，由于卫星传感器在获取地表信息时受到云覆盖、大气干扰、双向反射(太阳入射角、传感器观测角)等因素的影响，NDVI时间序列数据比一般数据更容易存在严重的噪声，很大程度上，限制了NDVI时间序列数据的进一步解译和更好的应用[6]。

因此在使用NDVI时间序列数据之前，必须要采取一定的处理方法对NDVI数据进行修复，以去除NDVI时间序列数据中的噪声。迄今为止，国内外学者提出了大量的遥感时间序列数据的去噪重建算法，这些方法可大致分为两类：阈值去除法，例如最佳指数斜率提取法(BISE) [7]以及基于滤波的去噪方法，例如Savitzky-Golay (S-G)滤波法[8]和小波低通滤波算法[9] [10]。后来的学者在这些算法的基础上又提出了联合改进算法，例如B-W算法[11]即最佳指数斜率提取法(BISE)和小波变换(WT)的结合。尽管这些方法对于噪声都能起到降低的效果，但每种算法还是有自身的局限性，在使用之前还需进一步讨论这些算法的适用情况，究竟哪种方法最优，国内外学者并未达成共识。

针对以上现有的去噪效果对比方法，本研究将通过目视判别法、定量分析法、和变化检测三种方法对目前三种常用NDVI时间序列数据的滤波去噪算法进行对比实验，并分析比较了不同算法的特点和去噪效果。

2. 研究区域概况及数据

宝日希勒露天矿区位于内蒙古自治区呼伦贝尔市陈巴尔虎旗境内，距离宝日希勒镇10.5 km，距离呼伦贝尔市海拉尔区20 km (见图1)。地势东北高西南低，属于“海拉尔内陆断陷盆地”(东经119˚41'56''~119˚49'23''；北纬49˚20'24''~49˚27'31'')，海拔高度为600.00~650.00 m，矿区内地形起伏较为平缓，多以丘陵为主。气候方面，宝日希勒矿区形成了独特的中温带半干旱大陆性草原气候：春季多大风而少雨，蒸发量大；夏季清凉且短促，降水集中但量较少；秋季降温快，霜冻早；冬季严寒漫长，无霜期短，平均为125天，地面积雪时间长。全年极端最高气温38.4℃，极端最低气温为−49.0℃。年平均气温为−2.6˚C，年平均降水总量290~350 mm，其中1月和7月的平均降雨量分别为3.79 mm和99.38 mm。植被方面，首先矿区处于呼伦贝尔草原的中部，是典型的草原生态薄弱区域。其次矿区内部地带性的植被多为草甸草原植被，包括羊草、冰草、披碱草等。植被生态适应性多以中旱生、旱生植物为主。

Figure 1. Map of the location of the study area

图1. 研究区的位置示意图

本文数据主要采用来自USGS美国地质调查局所发布的Landsat图像所提取的NDVI数据。

Landsat卫星旧称为地球资源技术卫星，是美国NASA陆地卫星计划中的系列卫星，其重复周期短、图像覆盖范围广。

对于每幅陆地卫星影像，NDVI都是使用红色和近红外波段按照公式(1)计算而来[12]：

$NDVI=\frac{{\rho }_{NIR}-{\rho }_R}{{\rho }_{NIR}+{\rho }_R}$ (1)

其中 ${\rho }_{NIR}$ 是植被在近红外波段的反射值， ${\rho }_R$ 是植被在红色波段的反射值，以上的研究区数据就是2002年至2019年间NDVI时间序列数据值的组合。

3. 研究方法

3.1. 滤波方法简介

本文使用了S-G滤波、BISE滤波和BISE-WT滤波3种NDVI时间序列数据去噪算法。下面从基本原理和数学模型两个角度，简单介绍这3种去噪算法及优缺点。

3.1.1. Savitzky-Golay滤波(S-G)

Savitzky和Golay [1]在1964年提出S-G滤波，该滤波又被称为最小二乘法或数据平滑多项式滤波。S-G滤波的核心思想是：通过随机选取某点周围的一部分点，应用最小二乘卷积算法拟合出一个n阶多项式，利用该多项式计算出这个点的平滑值。实际上简单来讲，就是对原有序列的一个加权平均，其中权重的大小由n阶多项式模型来决定。需注意一点，该算法只适用于云处理之后数据的去噪。NDVI时间序列数据的S-G滤波过程可以用公式(2)描述：

$Y_j^{*}=\frac{{\displaystyle {\sum }_{i=-m}^{i=m}{C}_i{Y}_{j+i}}}{N}$ . (2)

式中， $Y_j^{*}$ 代表拟合后的时间序列数据； $Y_{j+i}$ 代表原始时间序列数据； $C_i$ 指的是第i个时序数据经滤波时的系数；N为滑动窗口所包含的数据点。

在处理时间序列数据时，运用S-G滤波算法需要满足两个基本条件：1) 所处理的NDVI时间序列曲线是连续的，且服从植被的年际动态变化趋势。2) 云、大气和雨雾等外界因素会使NDVI值降低，因此绝大部分噪声会低于整体的NDVI序列数据的均值。

此外，在使用S-G算法时，需要事先人为确定两个参数，即滑动窗口的宽度以及多项式的拟合阶数，合理的参数能够保证NDVI时间序列数据拟合的精确性。倘若滑动窗口的宽度设定的偏小，则会出现大量的数据冗余，不利于获取数据集的长期变化趋势；若滑动窗口的宽度设定偏大，则会忽略很多细节信息。多项式的拟合阶数通常设置为2~4，若阶数偏低，则时序曲线会变得更加平滑，与此同时，会保存一些异常的数值；反之，能够去掉这些异常值，但也会出现过度拟合的现象，产生新的噪声。

3.1.2. 最佳指数斜率提取算法(BISE)

最佳指数斜率提取算法(BISE)是Viovy [9]等人于1992年提出的，用于去除NDVI时间序列中的噪声。这是一种基于阈值的方法，通常需要根据主观的经验策略，并取决于分析员的技能和经验，确定滑动期和滑动期内NDVI再生可接受百分比增加的阈值。Viovy等人认为NDVI时间序列曲线中的突变可能是由多种因素造成的。例如，有云到无云状况的转变，或视场角的切换。BISE方法的核心思想是：在一个滑动窗口内向前查询，若下一个点的值比起点的值高，则该点的值保持不变；若下一个点的值比起点的值低，且该低值点的1.2倍高于滑动窗口内的最大值，则该点的值不变；如果该低值点的1.2倍低于滑动窗口内的最大值，则需要用滑动窗口内的最大值代替该点的值。

BISE滤波法的关键是滑动窗口的大小，如果滑动窗口过小，则会包含大量的噪声；反之，如果滑动窗口过大，有些植被的生长速度较快，就会忽略掉一些重要的植被动态变化有效信息。VIovy等人经过大量的实验探究发现，选择30天作为滑动窗口是最理想的。

这里，我们再通过一个大小为3像素的滑动窗口进一步描述BISE算法，将第一个像素的值和第二个像素的值之间的减少值设为 $D_{i-1,j}$ ，将最后一个像素的值于第二个像素的值之间的减少值设为 $D_{i+1,j}$ ，分别按照等式(3)和等式(4)计算：

$D_{i-1,j}=\frac{x_{i-1}-x_i}{x_{i-1}}$ (3)

$D_{i+1,j}=\frac{x_{i+1}-x_i}{x_{i+1}}$ (4)

其中， $x_{i-1}$ 表示第一个像素的值， $x_i$ 和 $x_{i+1}$ 分别是第二个像素的值和最后一个像素的值。如果 $D_{i-1,j}$ 、 $D_{i+1，j}$ 均大于阈值， $x_i$ 将被重新分配如下：

$x_i=\frac{x_{i-1}+x_{i+1}}{2}$ (5)

随着滑动窗口向前移动，BISE将对序列进行去噪。

3.1.3. 小波变换法(Wavelet Transform)

Lu等人提出的小波变换法[10]是一种基于频域的滤波方法，其核心思想是：将连续的信号分解为不同频率分量的函数，每个分解后的量具有与其对应的尺度。经过这一分解，会将原时间序列信号分为低频和高频两部分，然后需要将低频部分进行再分解，再次分解为低频和高频部分，如此下去，最后通过低通滤波器滤除高频信号，得到平滑的NDVI时序曲线。

长度固定的且快速降低的振荡波形(母函数)在不同尺度上转换即得到小波函数(子函数)，而小波变换往往都是由小波函数进行表示的。小波的基本形式即称作母函数，通过对小波母函数 $\text{Ψ}\left(t\right)$ 。进行伸缩平移，函数可通过公式(6)得到：

${\text{Ψ}}_{a,\tau }\left(t\right)=\frac{1}{\sqrt{a}}\Psi \left(\frac{t-\tau }{a}\right)$ (6)

式中：a,τ∈R；a > 0，a为尺度因子， $\tau$ 为位移参数； $\text{Ψ}\left(t\right)$ 为母函数。

对于函数f(t)的连续小波变换其表达式为(7)：

${\Psi }_{a,\tau }\left(t\right)\ge \frac{1}{\sqrt{a}}$ $\underset{R}{\overset{ }{{\displaystyle \int }}}f\left(t\right){\Psi }^{*}\left(\frac{t-\tau }{a}\right)\text{d}t$ (7)

式中：*代表共轭函数。

3.2. 定性分析

为比较S-G滤波算法、BISE算法以及B-W算法的NDVI时序曲线平滑算法在去噪效果上的差异，本研究提取了Landsat NDVI最大合成数据中的5组样点数据进行对比分析。由于空间限制，本研究图中显示的为宝日希勒矿区2000~2021年的曲线数据。

宝日希勒矿区植被覆盖较为分散，这点可从图2中明显看出，原时间序列数据曲线上下波动较大，NDVI值在0.1~0.7之间浮动，2008年NDVI值猛降至0.3，2013年降至0.1。

从整体看来，图2中3种去噪重建算法针对原数据得到的NDVI曲线的变化趋势都较为平滑，说明S-G滤波、BISE滤波和B-W滤波均具有一定的去噪能力。其中较为显著的噪点如2008年和2013年，在图2中可以看到经过三种算法去噪后，NDVI值都有明显的抬升，说明在重建时均已被剔除。

从图2(a)中可以明显看到，BISE算法消除了2008年的显著噪声，并保留了因地表开采造成的NDVI的突变的情况，与S-G滤波和B-W滤波相比，它与原图的整体趋势基本相同。BISE得到的NDVI时序曲线较为粗糙，保留了连续噪声，与原始数据之间的差异最小，说明其去噪效果一般，仅实现了初步去噪；从图2(b)、图2(c)中可以看出S-G滤波在很大程度上去除了噪声，同时去除了2004、2008、2013年的连续噪声，得到的NDVI时序曲线比其余两个方法得到的曲线更加平滑，存在过度拟合的现象，这容易造成波峰部分高值低估、波谷部分低值高估的情况，去噪效果较差；从图2(c)中可以看出，BISE-WT去噪方法通过贝叶斯风险函数计算了小波阈值，将幅值差异较大的信号分割开来，然后将幅值较低的噪声信号剔除，同时将幅值较大的有效信号予以保留。即很好地综合了BISE滤波保留原曲线趋势较好以及WT保证曲线平滑度的特点，不仅去除了噪声，更保留了曲线原本的趋势，去噪效果最好。具体效果图如下，其中数字代表测试像素的顺序。

Figure 2. The NDVI time series denoising curves of BISE, S-G and B-W for test pixels

图2. BISE、S-G和B-W对测试像素的NDVI时间序列去噪曲线

3.3. 定量分析

为了定量评价3种方法对NDVI时序数据的去噪重建效果，本研究选取400个测试像素进行RMSE分析。表1显示了噪声等级分别为10%、40%、70%时，3种去噪方法的RMSE值。

对于Landsat-NDVI数据而言，不同平滑滤波的NDVI重建效果具有一定的差异。从表1中可以看出，与对照组相比，三种方法的RMSE值均低于对照组，证明三种去噪方法都可以抑制原NDVI时序数据的噪声；对于噪声为10%的NDVI时序数据，BISE的RMSE值相比其他两种方法更低，去噪效果最好；对于噪声为40%和70%时的NDVI时序数据，B-W去噪方法的RMSE值更低，该方法的去噪性能最好。同时，无论噪声百分比是多少，S-G滤波的RMSE值都是三者中最高的，证明其在三者中去噪效果最差。

Table 1. RMSE values of three denoising methods

表1. 3种去噪方法的RMSE值

3.4. 变化检测分析

在存有噪声影响时，会使得原本因开采导致NDVI值下降的年份提前或向后推迟，去除噪声后，则会恢复到原本应该产生变化的年份。基于此原理，进行变化监测分析3种方法的去噪效果。图3显示的为研究区域用于变化检测分析后的结果。其中0值代表的为2000年之前的采矿情况。我们通过可视化解释确定了所选样本点的实际扰动年数，并将上一年因非植被而变化的原始NDVI最大合成影像的检测日期定义为每个像素的干扰年。从图3中，我们可以注意到，3种去噪方法得到的结果与原始NDVI最大合成图像得到的结果非常不同。从(d)图中可以看出原图已经受到了噪声的影响，2000年之前开采的地方，有一些被识别成其他的年份。与(d)图相比，图(b)所显示的S-G滤波识别变化年份时出错较多，过度拟合现象较为严重、去噪效果较差；图(c)所显示的B-W算法，可以看出部分年份被识别成其他年份的情况较其他两种方法少，同时平滑了因2000年以前因开采导致的噪声。

为了进一步评估变化检测结果，我们从样本区域随机选择了200个样本点，计算各种方法的整体精度，具体结果如表2所示。不同的滤波方法均存在着过度拟合的现象，可以看到经三种方法去噪后，变化检测结果的整体精度反而小于去噪前的影像，其中S-G滤波的过度拟合现象最为严重，分类精度最低，B-W的整体分类精度最好。

4. 结论

本研究以宝日希勒矿区为研究区，从定性分析、定量分析、以及变化检测三个角度对比分析了Savitzky-Golay滤波(S-G)、最佳指数斜率算法(BISE)以及BISE-WT算法对NDVI最大合成时间序列的去噪重建效果，研究结果如下：

1) 3种平滑算法都有效地对原NDVI时序数据进行了去噪，且效果明显；2) 对于本文所选的NDVI最大合成数据而言，通过三方面的对比分析，3种去噪算法中B-W算法的去噪重建效果最好，定性分析中其保留了与特征变化相关的真实变化信息，同时去除了噪声；定量分析中其RMSE值最小与理想NDVI

图中(a)为BISE滤波的变化检测结果；(b)为S-G滤波的变化检测结果；(c)为B-W滤波的变化检测结果；(d)为NDVI最大合成影像的变化检测结果。

Figure 3. Three kinds of noise filtering and change detection results of the original image

图3. 3种去噪滤波以及原图的变化检测结果

Table 2. Change detection and analysis results of raw and de-noised NDVI images

表2. 原始和去噪后NDVI影像的变化检测分析结果

序列最接近；变化检测中，其识别变化年份的精度最为准确。其次去噪效果一般的为BISE算法，最后去噪效果最差为S-G算法。

参考文献