王玮，郭铌，沙莎，胡蝶，王小平，李耀辉

(中国气象局兰州干旱气象研究所，甘肃省干旱气候变化与减灾重点实验室，中国气象局干旱气候变化与减灾重点实验室，甘肃 兰州 730020)



我国西北地区东部时间序列NDVI数据集重建方法比较研究

王玮，郭铌*，沙莎，胡蝶，王小平，李耀辉

(中国气象局兰州干旱气象研究所，甘肃省干旱气候变化与减灾重点实验室，中国气象局干旱气候变化与减灾重点实验室，甘肃 兰州 730020)

高质量、长时序归一化植被指数(NDVI)数据集不仅是连续监测陆地表面特征的基础，也是研究气候与陆地生态系统变化的重要参数。本研究以生态环境较为脆弱的西北地区东部为例，借助多种时间序列重建方法对LTDR NDVI数据集中的噪声进行拟合重建，并结合农业气象资料和高质量NDVI数据，对不同重建方法的拟合结果开展适用性评价分析，结果表明，1)下垫面类型是影响重建方法拟合效果的重要因素。根据不同植被类型或作物生长特点，每种重建方法对其噪声消除能力有所不同；2)在年均NDVI较高(NDVI≥0.3)，且NDVI曲线具有明显季节变化的草地、林地以及牧草等作物种植区域内，经过D-L拟合重建的NDVI具有较高的保真能力和适应性；3)在年均NDVI较低(NDVI<0.3)，且植被季节生长变化不明显或NDVI曲线不呈季节对称性变化的稀疏植被区，以及以冬小麦为典型作物种植的区域内，经过S-G滤波重建的NDVI数据表现出相对较好的保真能力和适应性。

NDVI；时间序列重建；植被遥感；AVHRR

植被在地球系统中扮演着重要的角色，是气候和人文因素对环境影响的敏感指标[1-2]。遥感技术的快速发展为人类准确获取陆地表面大范围、多尺度的植被信息提供了强有力的手段。利用可见光/近红外波段形成的归一化植被指数(NDVI)已被广泛地应用于区域和全球环境变化、农作物估产、干旱监测以及模拟新的植被指数等方面，是全球应用最广的一种植被指数[3-6]。

获取高质量、长时序NDVI数据不仅是长期监测陆地表面特征的基础，也是研究气候与陆地生态系统变化的前提。由于卫星传感器成像本身具有瞬时性和周期性(卫星轨道周期和卫星生命周期)的特点，同时还会受到潜在噪声的影响(如云、气溶胶等)，因此在使用卫星传感器来获取NDVI时间序列数据时，数据的质量和连续性将会受到极大的挑战[7]。特别是在大空间尺度的时间序列监测研究中，这一问题尤为严重。研究表明，大气中的水汽、臭氧、气溶胶和瑞利散射等会影响红光波段与近红外波段的反射率，使得卫星传感器在接收到来自地面目标信号的同时，还会收到部分大气噪声信号，导致NDVI数值与实际情况出现较大的偏差[8-10]。尽管研究对象是植被，提取的是NDVI数据，但在植被生长的不同阶段，传感器实际接收到的信息也包括除植被以外的土壤背景。而土壤背景影响的本质是土壤本身对红光和近红外波段具有不同的反射率，包括潮湿地面、雪、枯叶、粗糙度、土壤类型等因素[11]。同时，在每次观测时“太阳-地物-传感器”三者之间的关系变化都可能直接影响植被指数计算带来的不确定性[12-13]。因此，建立高质量、长时序的NDVI数据集是进一步开展植被定量遥感监测、气候变化及陆地生态系统响应等相关研究的关键。

虽然目前常用的时间序列NDVI数据集在发布前都经过了较为严格的预处理来控制数据的质量，但由于受云、气溶胶、土壤背景、太阳/观测角度、数据传输错误、植被空间结构以及轨道和传感器信号衰退等因素的影响，使这些数据集中仍然残留较多噪声，导致时间序列NDVI曲线出现陡升或者陡降的情况，从而影响了数据集的质量与连续性，造成在宏观生态系统研究中基础观测数据存在不可靠的问题[14]。尤其是AVHRR传感器，在设计之初并不是以监测地表特征为目的，导致其生成的NDVI时间序列数据更容易受到潜在噪声的影响[15-17]。作为AVHRR NDVI的完善和升级，自2000年美国航空航天局开始免费分发全球MODIS NDVI产品后，使其快速成为应用研究的热点[18-19]。虽然MODIS NDVI具有较窄的波段设置和较高的精度，但时间序列长度有限，在年代际尺度的变化趋势研究中，仍然需要结合AVHRR数据对其进行时间序列的扩展和完善[20]。然而由于AVHRR数据已经积累了35年的时间长度，是目前时间序列跨度最长的光学遥感数据，对研究全球变化具有重要作用。因此，只有不断订正AVHRR NDVI数据集中的噪声，插补缺测资料，发挥该数据长时序的优势，为相关领域研究提供高质量的NDVI数据集将具有十分重要的现实意义。

时间序列重建是利用多种统计和数值分析方法，实现数据的插补和去噪，并提高数据质量的一种手段，目前已被成功地应用于NDVI数据时序重建工作中[21-22]。迄今为止，国内外学者在NDVI时间序列数据集重建方面开展了大量研究，目前时序重建方法已超过20种，依据计算方式可将其归并为：阈值法、滤波法、函数拟合法以及综合方法等[8,11]。虽然重建方法众多，但每种重建方法对不同性质噪声点去噪、平滑和保真能力不同。Hird等[23]通过6种重建算法对MODIS 16d NDVI产品进行重建比较分析发现，6种重建方法均有较好的去噪能力，但随植被类型的变化重建效果具有明显的差异性。黎治华[24]利用多种重建算法相结合的方案对MODIS的时间序列产品进行重构后发现，不同重建算法各有优劣，相关研究工作还存在很大的不确定性。由此可见，多种重建算法的提出，虽然为长时间序列NDVI数据集的构建提供了良好的技术支持，但对这些重建算法进行系统分析和比较的研究工作还不够深入。对不同植被覆盖或作物类型而言，研究者很少关注已有重建方法的适用性，在重建效果定量分析和精度验证方面还需进一步开展相关研究工作。因此，系统比较和分析常用重建方法的拟合效果，开展以业务应用为目的重建算法适用性评价研究是很有必要的。

2014年6月美国戈达德太空中心发布了V4版本的LTDR(Land Long Term Data Record)产品，其中包括1981-2010年的地表反射率数据，这为建立高质量、长时序的NDVI数据集奠定了较好的基础。因此，本研究以生态环境相对脆弱的西北地区东部为例，借助时间序列重建软件(TIMESAT)[25]，通过结合农业气象台站资料，在系统比较与分析常用重建算法拟合效果的基础上，对不同重建方法展开适用性评价，确定研究区内最佳拟合重建方法，以期为建立高质量、长时序的NDVI数据集提供参考。

1　材料与方法

1.1研究区概况

本研究以地处干旱区和半干旱区过渡带的我国西北地区东部省份作为研究区。该地区介于 32°44′11″ N-44°07′24″ N、 88°47′17″ E-112°43′00″ E 之间，面积为1.48×106km2，约占我国陆地总面积的1/6。我国西北地区东部幅员辽阔，自然条件恶劣、生态环境较为脆弱，以山区降水和冰川融水补给为主，年降水量空间分布极不均匀，从南部450 mm以上过渡到北部沙漠地带的50 mm以下，耕地灌溉主要依靠自然降水。特殊的地理位置和水热气候条件，造成研究区自然植被空间分布差异较大，形成了由西南草地到东北荒漠逐渐倾斜过渡的自然植被空间分布特征(图1)。尤其是自20世纪70年代中期全球变暖加剧以来，该地区气候条件复杂多样，极端水文事件增加，地表植被覆盖对环境变化敏感。因此，建立高质量、长时序的NDVI数据集，为准确监测该地区的植被生长状况，了解气候变化对植被的影响，维护干旱半干旱区生态系统的稳定具有重要指导意义。

图1　西北地区东部植被分类图Fig.1　The classification of vegetation in the east of Northwest China　SⅠ：牧草 Pasture；SⅡ：春小麦+玉米Spring wheat+Corn；SⅢ：冬小麦+玉米Winter wheat+Corn；SⅣ：冬小麦Winter wheat.

1.2研究数据与预处理

1.2.1LTDR地表反射率数据通过美国航空航天局戈达德太空中心(http://ltdr.nascom.nasa.gov/cgi-bin/ltdr/ltdrPage.cgi)，下载了1981-2010年的LTDR 每日地表反射率产品(AVH09 Surface Reflectance Product)。该数据是以NOAA07(1981-1985年)、NOAA09(1985-1988年)、NOAA11(1988-1994年)、NOAA14(1995-1999年)、NOAA16(2000-2006年)和NOAA18(2005-2010年)上的AVHRR观测资料为基础，通过轨道筛选、辐射定标、云检测、大气校正、卫星漂移校正及双向反射率分布函数(bidirectional reflectance distribution function, BRDF)处理后生成的全球逐日格网数据集。数据分发采用与MODIS相类似的业务流程，以HDF分层格式存储，包括1～3波段地表反射率、4～5波段大气上界亮温、太阳高度角、传感器高度角、相对方位角以及质量控制文件。其中地表反射率数据是以空间分辨率为0.05°(约5 km)的格网进行存储，质量控制文件对每个格网制定1～10级的质量标记。为了便于计算，本研究利用ENVI 5.1软件将1～2波段地表反射数据处理为Geo-TIFF格式，投影方式转为Albers，空间分辨率定义为5 km。

1.2.2植被类型数据植被类型数据来源于2001年中国科学院中国植被图编辑委员会发布的中国植被数据集。该数据集是由我国著名植被生态学家侯学煜院士主编出版的《1∶1000000 中国植被图集》数字矢量化而来，目前在国家自然科学基金委员会资助的“中国西部环境与生态科学数据中心”(http://westdc.westgis.ac.cn)提供下载共享，数据格式为Shapefile格式，投影方式为Albers。中国植被图是开展自然资源和自然地理特征研究的基本图件，它全面反映了我国11个植被类型组、54个植被型的796个群系和亚群系植被。在我国西北东部地区形成了以草地或草甸(以下简称草地植被)、稀疏植被、林地为主的植被分布格局(图1)。

1.2.3气象数据通过中国气象科学数据共享服务网(http://cdc.nmic.cn)，下载了研究区101个农业气象站在1992-2010年的逐旬资料。资料主要包括种植作物名称、发育期名称、发育期时间、发育程度等观测指标。经初步统计分析表明，研究区主要种植作物分为4种类型，分别是牧草、冬小麦(Triticumaestivum)、冬小麦+其他种植作物[玉米(Zeamays)、马铃薯(Solanumtuberosum)等]、春小麦(Triticumaestivum)+其他种植作物(玉米、一季稻等)。

1.2.4典型作物样区数据根据研究区作物种植特点，分别选取资料记录较为完整，且具有典型代表性的青海河南SⅠ(牧草)、甘肃张掖SⅡ(春小麦+玉米)、甘肃西峰SⅢ(冬小麦+玉米)和陕西蒲城SⅣ(冬小麦)4个农业气象台站资料作为本研究结果的验证数据。与此同时，为了尽量避免混合像元等不确定性因素对验证结果的影响，本研究利用地面实测资料和Landsat8 OLI高分辨率资料，选取了台站附近地势较为平坦、耕地分布均匀，种植作物类型与台站记录相同、且种植面积大于5 km×5 km的下垫面作为典型样区(图1)。

1.3研究方法

虽然AVHRR传感器均可提供时间分辨率为1 d的数据，但逐日数据易受云、气溶胶及坏线等因素的影响。因此，本研究首先提取1981-2010年LTDR地表反射率产品中的1～2波段资料生成每日NDVI数据，其次结合数据的质量控制文件，根据最大合成法(max value composite, MVC)[26]，将逐日NDVI数据集中的明显噪声进行初步去除，并生成逐旬初始NDVI数据集，然后利用TIMESAT软件对初始LTDR NDVI数据进行时序拟合重建，最后结合典型样区资料和同时期高质量NDVI数据，对不同重建方法开展定量分析和适用性评价，确定研究区不同植被或作物类型上的最佳重建拟合方法。具体流程如图2。

图2　技术路线图Fig.2　The flow chart of research method

1.3.1初始NDVI数据合成本研究利用LTDR 1～2波段地表反射率资料，根据NDVI计算公式生成1981-2010年每日的NDVI数据(O-N)，公式如下：

NDVI=(ρNir-ρRed)/(ρNir+ρRed)

(1)

式中，ρNir和ρRed分别为红光波段和近红外波段的反射率，LTDR资料相对应为第1波段和第2波段的地表反射率值。

结合数据自带的质量控制文件(quality assessment, QA)，根据MVC合成法[26]，将每日LTDR NDVI合成初始逐旬NDVI数据，并生成相应时期的QA文件。QA文件是在像元尺度上能够综合评价和反映该时段遥感数据质量信息的资料。

1.3.2基于TIMESAT软件的时间序列重建方法TIMESAT软件是由Jönsson等[25]共同开发，并用于植被指数时间序列数据集重建及提取植被生长物候信息的程序包。软件集成了目前已被国内外学者广泛应用的非对称高斯函数拟合(asymmetric Gaussian model function)、双逻辑斯蒂拟合(double Logistic function)以及Savizky-Golay 滤波3种重建算法，并取得较好的效果[27-28]。该软件有专门网页介绍相关信息(http://www.nateko.lu.se/timesat/timesat.asp)，可同时处理时间序列数据和影像数据(二维空间阵列)。最初的版为FORTRAN语言编译，之后由FORTRAN和MATLAB两种语言编译成升级兼容版本，在功能上有了较大的提升。主要包括数据预处理、数据处理以及后处理3个大模块。可结合数据的质量控制文件实现不同遥感时序数据的重建。在数据预处理模块中可以选择其中一种或几种重建算法进行重建，并可在同一窗口中显示不同重建方法的预览效果。TIMESAT软件最大的特点是可根据不同的植被覆盖类型设置不同的重建算法及参数，提取不同植被特征下的物候参数。

1)非对称高斯函数拟合法(A-G)。A-G拟合法是一个从半局部拟合到整体拟合的过程，该算法可分为区间提取、局部拟合和整体连接三部分[22]。首先根据时间序列点相对应的值来提取区间，然后使用局部拟合函数即高斯函数对每个区间的时间序列数据进行拟合，最后通过整体拟合实现时间序列的重建，具体步骤如下。

第1步：通过局部拟合公式进行拟合，公式如下:

f(t;c1,c2,a1,……，a5)=c1+c2g(t;a1,……，a5)

(2)

式中，f(t;c1,c2,a1,……，a5)为拟合结果；c1和c2控制曲线的基准和幅度；g(t;a1,……，a5)为高斯函数，计算公式为：

(3)

式中：a1决定最大值和最小值的位置；a2、a3、a4、a5控制曲线左半部分和右半部分的宽度和陡峭度。

第2步：参数确定，通过极小化局部拟合公式进行拟合计算，公式如下：

(4)

式中：σi为第i个数据点的测量不确定性，可认为已知(通过质量控制文件获取)，Ii为第i点的观测值。

第3步：调整拟合曲线至植被生长的上包络线：若拟合曲线的NDVI值低于原始观测值，则通过调整低值拟合点的σi值进行再次拟合，从而使拟合曲线接近植被生长曲线的上包络线。

第4步：通过整体拟合实现时间序列重建，将局部拟合的左边最小值，中间最大值和右边最小值分别表示为fL(t)，fC(t)和fR(t)，得到一个VI*(t)函数来纠正局部拟合函数，拟合公式为：

(5)

式中：α(t)和β(t)为裁切函数。

2)双逻辑斯蒂拟合法(D-L)。D-L算法与A-G拟合法类似，也是一种半局部拟合方法。首先将整个时间序列中时间点对应的值按极大或极小值分成多个区间，分别对各区间进行双逻辑斯蒂函数局部拟合，其局部拟合方式与非对称高斯拟合方法类似，公式如下：

f(t;c1,c2,a1,……，a4)=c1+c2g(t;a1,……，a4)

式中：g(t;a1,……，a4)根据双逻辑斯蒂函数拟合得到，计算公式为：

(6)

式中：a1和a3决定曲线左边和右边变化的位置，a2和a4决定曲线左边和右边的变化率。以上参数的获取通过极小化x2获得，如下式：

(7)

最后利用公式(5)进行整体拟合，并实现长时间序列的重建。

3)Savizky-Golay 滤波法(S-G)。S-G滤波是一种权重滑动平均滤波方法，主要是基于最小二乘卷积拟合法来实现噪声消除的。其权重取决于一个滤波窗口内做多项式最小二乘拟合的次数，如下式：

(8)

式中：Y为原始的NDVI值；Y*为平滑后的NDVI值；Ci为平滑窗口中第i个NDVI值的权重系数；N为卷积数，等于平滑窗口的大小(2m+1)；指数j是原始坐标数据的滑动指数；m为滑动窗口的一半宽度。Ci可根据Madden等[29]提供的等式计算得到。具体步骤如下：

第1步：噪声NDVI值的线性内插，通过质量控制文件识别标记受残留噪声影响的数据，使用线性内插对这些数据进行修正，得到新的时间序列数据Ni0。

第2步：使用S-G滤波拟合得到长时间变化的趋势曲线Nitr；通过选择适当的参数值m和j拟合NDVI序列长期变化的趋势线。

第3步：确定NDVI时序每一个点的权重：根据以上两步得到的新时间序列数据Ni0和序列的趋势曲线Nitr来确定新时间序列上每个点的权重Wi，计算公式如下：

(9)

式中：Ni0和Nitr分别为新序列和趋势线上第i点的NDVI值，di=|Ni0-Nitr|，dmax为di的最大值。

第4步：根据Ni0和Nitr生成新的时间序列Ni1，用于修正Ni0中存在的某些明显的噪声点，生成的新序列更真实地反映NDVI的变化趋势，接近原始NDVI数据的上包络线。计算公式为：

(10)

第5步：使用S-G滤波对Ni1曲线进行拟合，得到新的拟合曲线Ni2。

第6步：计算拟合效果指数Fk，用于估计拟合过程中拟合的NDVI值接近较高权重值的程度，拟合公式为：

(11)

本研究在利用TIMESAT软件对初始逐旬LTDR NDVI数据进行A-G、D-L和S-G 3种拟合重建处理时，还需要对一些相关处理参数进行调试，主要包括：NDVI有效值域(range，0～10000)、噪声去除阈值(spike，2)、滑动窗口大小(w-size，5～6～7)、拟合峰值参数(altitude，5)、迭代次数(fitting steps，3)等。

1.3.33种重建方法拟合效果比较分析为了定量分析以上3种重建方法的拟合效果和保真性能，本研究分别对典型耕地样区和不同植被类型区的时序重建NDVI数据进行分析与比较。

1)典型耕地样区拟合前后物候期比较分析。已有研究表明NDVI对植被的物候期非常敏感[30]，当作物从一个物候阶段转向另一个阶段时，NDVI曲线曲率的变化能够被用来确定物候期的转变[30-32]。Zhang 等[33]将植被返青日期定义为NDVI的曲率导数首次达到局部最大值的日期。在作物生长周期内，通过计算NDVI拟合曲线的曲率，并求其极值，由此能够确定植被生长阶段转变的时间起点有6个，分别对应返青、生长速率加快、成熟、衰老、衰老速率加快和枯黄死亡期(图3)。曲率计算公式如下：

(12)

式中：da是沿时间曲线移动单位弧长时切线转过的角度；ds为单位弧长；z=exp(a+bt)，a和b为时间t的拟合参数；c为NDVI最大值。

因此，根据以上植被物候期特征遥感识别方法[33]，本研究在典型样区分别提取A-G、D-L、S-G以及初始NDVI时间序列曲线物候特征信息的基础上，利用蒲城、海宴、河南、西峰和张掖4个样区农业气象台站实际观测的作物生育期资料，分析比较每个典型区内不同NDVI时间序列重建曲线表征的物候期(如返青期、抽穗和枯黄期)时间节点与农气台站实测作物生育期之间的误差。为了使误差定量可比，本研究将农气台站记录的作物返青、成熟和枯黄的日期作为实测值，将通过NDVI重建曲线提取到的各相应生育期的时间节点作为模拟值，最后利用均方根误差(RMSE)和平均绝对误差(Pd)等相关统计量，定量分析模拟值与实测值之间的误差精度。

其中，RMSE可以描述重建前后NDVI数值之间的平均差异程度，其值越小，拟合值的代表性越强；Pd能够深入反映各算法彼此间保真能力的差异。RMSE和Pd公式如下：

(13)

(14)

式中：xi、yi分别表示实测样本数据及其对应的重建样本数据；Zai为第i个样点的实际观测值，Zei为重建拟合值；n为用于验证的样本数目。

2)不同土地类型区误差定量分析。参照植被覆盖类型和数字高程模型，在研究区下垫面较为均一、地势较为平坦的草地、耕地、林地以及稀疏植被区分别选取初始NDVI时序数据中高质量的样本数据(即QA=1的样本数据)，并利用这些样本数据定量分析重建前后NDVI数值之间的相关系数和均方根误差等相关统计量。

2　结果与分析

2.1典型作物样区拟合重建结果特征比较

通过分析1981-2010年青海河南(SⅠ)、甘肃张掖(SⅡ)、甘肃西峰(SⅢ)和陕西蒲城(SⅣ)4个典型作物样区内的时间序列NDVI拟合重建数据，结果表明，3种拟合重建算法在对初始NDVI像元中的残留噪声进行消除的同时，也将典型像元中的时序NDVI进行了不同程度的拟合调整。其中经过D-L和A-G方法拟合曲线的结果比较接近，仅在NDVI峰值处两者之间的拟合结果略有不同，相对误差为±0.7%(图4)。S-G滤波法对噪声比较敏感，其重建的NDVI值与前两种拟合重建结果有较大的差异，尤其是在生长季开始和结束时段，S-G拟合曲线相对于初始(O-N)、D-L和A-G拟合曲线的峰值存在较大程度的“左偏”和“右移”的情况(图4)。

植被从返青期到枯黄期阶段，经过3种重建算法拟合后的NDVI整体较初始NDVI值提高了11.6%～23.7%，且重建序列曲线“拐点”对应的时段也不一致。在生长季阶段(4-10月)，青海河南和甘肃张掖典型样区的初始NDVI曲线呈较为明显的年内对称性，且3种拟合重建后的结果也较为接近。但以冬小麦为主要作物种植区的陕西蒲城和甘肃西峰(图4)，在生长季阶段的NDVI曲线具有明显的年内不对称性，波峰波谷更替较为频繁，且初始NDVI曲线与3种拟合重建后的曲线具有较大差异。由此可见，不同作物的种植结构对拟合重建效果具有较大影响。

2.2典型作物样区拟合重建结果误差分析

通过比较分析1992-2010年各典型区初始NDVI曲线(O-N)，以及3种重建NDVI曲线曲率所表征的物候期(返青、成熟和枯黄期)节点与相应台站记录作物生育期之间的差异，可以看出：由于受到噪声的影响，初始NDVI曲线反映的各物候期节点与台站记录的相应生育期之间存在较大误差，且随着作物不断生长，即由返青期到枯黄期阶段，平均绝对误差Pd和均方根误差RMSE也随之升高，其中Pd和RMSE分别介于2.36～3.68和1.73～5.95之间。

与初始NDVI曲线表征的各物候期相比，经过拟合重建后的3种NDVI曲线所反映的物候期节点与台站记录的相应生育期之间的各项误差均有所降低，其中在青海河南(SⅠ)和甘肃张掖(SⅡ)典型区域内，通过D-L重建后的NDVI拟合曲线反映的各生育期节点与台站记录的相应时期最为接近，Pd和RMSE达到最小，分别在1.89和4.21以下；在陕西蒲城(SⅣ)以及甘肃西峰(SⅢ)典型区域内，经过S-G滤波重建的NDVI曲线所反映的各生育期节点与实测值之间具有较小Pd和RMSE，分别介于2.35～2.47和1.57～4.69之间(表1)。

由此可见，3种拟合重建算法均可以对初始NDVI数据中的残留噪声进行消除，但根据不同作物生长特点或作物种植结构，每种重建方法对其噪声消除能力是不同的。其中在生长季阶段(4-10月)，作物生长曲线具有明显年内对称性的SⅠ和SⅡ典型区域内，经过D-L拟合重建结果在相应生育期内的各项误差较小，其适应性较好；而在以冬小麦为典型种植区的SⅢ和SⅣ区域内，经过S-G滤波拟合重建的NDVI曲线对各生育期节点的反映与实测值更为接近，具有较好的适应性。

2.3不同植被或作物类型NDVI拟合重建结果误差分析

根据NDVI质量控制文件，在获取1981-2010年各植被类型区高质量(QA=1)NDVI样本数据，以及经过上述3种拟合重建后同期NDVI样本数据的基础上，通过比较分析重建前后NDVI数据之间的相关系数和RMSE(图5)，结果表明，在草地(Ⅰ)和林地(Ⅲ)为主的自然植被区域内，以及牧草(Ⅳ)和春小麦+其他作物(Ⅴ)的耕地作物区域中，经过D-L拟合重建的NDVI数据与重建前同期高质量NDVI数据之间的相关系数达到最高(≥0.680)，且RMSE较低(≤0.026)。经过A-G拟合后的结果与D-L较为接近，而S-G拟合结果与同期高质量NDVI数据之间具有较大差别，其相关系数和RMSE分别介于0.620～0.780和0.019～0.028之间。与此同时，在稀疏植被区(Ⅱ)、冬小麦+其他作物(Ⅵ)和冬小麦(Ⅶ)区域中，经过S-G滤波拟合重建后的NDVI数据与同期重建前高质量NDVI数据之间具有较高的相关系数(≥0.56)和较低的RMSE(≤0.027)。

图4　1981-2010年典型样区逐旬NDVI时间序列重建前后对比Fig.4　Comparison of original and reconstructed NDVI of ten-days period from 1981 to 2010 in sample areas

表1　典型作物样区NDVI曲线物候期误差分析

图5　不同植被或作物类型区内高质量NDVI重建前后保真能力分析Fig.5　Comparison of the ability on keeping the main characters by NDVI in different vegetation coverage areasⅠ～Ⅶ分别表示草地、稀疏植被、林地、牧草、春小麦+其他作物、冬小麦+其他作物和冬小麦区Ⅰ to Ⅶ are respectively grassland, sparse vegetation, woodland, pasture, spring wheat+other crops, winter wheat+other crops, and winter wheat.

由此可见，随着下垫面植被或作物类型的变化，3种拟合重建结果具有较大的差异性。其中，在年均NDVI较高(NDVI≥0.3)的草地、林地区、牧草和春小麦+其他作物区域中，经过D-L拟合重建的NDVI数据具有较高的保真能力和适应性，而在年均NDVI较低(NDVI<0.3)的稀疏植被和冬小麦+其他作物区域中，S-G滤波重建方法表现出相对较好的保真能力和适应性。

3　讨论

综上所述，由于目前已发布的时间序列 LTDR NDVI数据中仍然残留较多噪声，会造成NDVI曲线出现陡升或者陡降的情况，由此导致植被生长曲线表征的各物候期时间节点与台站实际记录的生育期之间存在较大偏差。因此，在实际应用过程中需要对该数据集中的残留噪声进行消除。利用时间序列重建算法对NDVI数据集进行拟合重建，不仅能够消除数据中的残留噪声，同时还可以插补缺测数据，是提高数据集质量，获取连续时间序列的一种有效手段。

根据本研究可以看出，当植被生长具有明显的季节变化特征时(如草地或部分作物区域)，利用D-L或A-G拟合方法重建的NDVI时间序列曲线，对返青期、成熟期以及枯黄期等时间节点的反映较为准确，并与重建前高质量NDVI数据具有较小的误差和较高的保真性。这主要是由于A-G与D-L拟合方法均是一个从半局部拟合到整体拟合的过程，两者拟合时序数据的出发点是相同的，只是后期分别使用了高斯函数和逻辑斯蒂函数进行拟合。因此，经过以上两种方法得到的NDVI曲线结果较为接近，仅仅只是对NDVI峰值的拟合存在差异。而S-G滤波重建主要是利用移动加权滑动平均滤波对整个NDVI时序数据的长期变化趋势进行拟合，并将NDVI值分为“真值”点和“伪真值”点两类。最后通过局部迭代拟合的方式将“伪真值”点用滤波值进行代替，从而得到新的拟合曲线。与A-G与D-L拟合方法相比较，S-G滤波方法对曲线局部细节的波动更为敏感，因此，在植被季节生长变化特征不明显的荒漠和部分耕地作物地区，利用S-G滤波方法重建的NDVI曲线具有较好的拟合效果。Hird等[23]在加拿大落基山脉地区，利用MODIS 16d NDVI产品进行NDVI重建比较分析发现，A-G和D-L拟合结果优于S-G。曹云锋等[27]在长白山自然保护区进行NDVI重建比较研究也得出了以上相似结果。而在本研究中可以看出，只有在植被生长具有明显季节变化特征的草地、林地、牧草和部分作物类型上，经过A-G和D-L拟合方法才具有较好的拟合效果，在其他植被类型条件下，S-G滤波法均优于A-G和D-L的拟合效果。宋春桥等[28]对我国西藏北部地区的MODIS NDVI数据进行重建分析发现，A-G和D-L方法在草原和灌丛类型上具有较好的重建效果，相对而言，S-G滤波对噪声比较敏感，在处理局部较强的噪声时，可能使其“过度”拟合而引入新的噪声。本研究一方面通过重建前高质量NDVI数据和GIS空间分析等技术手段，从整个空间尺度上对上述观点进行了定量化的分析与比较，另一方面利用相应时期的农业气象资料，反映出植被季节生长曲线变化特征对不同重建方法的拟合重建效果具有重要影响。由此可见，在系统比较和分析多种重建方法的基础上，融合多种重建方法的优势，构建一种综合重建方法是下一步需要完善的工作。

由于本研究使用的农气资料记录时间有限，而无法对1981-1991年的NDVI数据进行定量分析，同时在自然植被类型区，由于缺乏长期的地面观测资料，造成精度评价并不完整。因此，收集更长时间序列的相关实测资料，从各个方面定量分析评价不同拟合重建方法的结果，对建立高质量、长时序的NDVI重建数据集具有重要的现实意义。

4　结论

本研究以生态环境相对脆弱的西北东部地区为例，结合农业气象台站资料和高质量NDVI数据，在系统比较与分析A-G、D-L和S-G 3种重建算法拟合效果的基础上，对不同重建方法展开适用性评价，得出如下结论：

1)下垫面类型是影响不同重建方法效果的一个重要因素。经过3种拟合重建算法获取的NDVI数据，均可以对初始NDVI数据中的残留噪声进行消除，并提高数据集的质量，但根据不同植被类型或作物生长特点，每种重建方法对其噪声消除能力是不同的。

2)在年均NDVI较高(NDVI≥0.3)，且具有明显季节变化特征的草地、林地、牧草以及春小麦+其他作物区域中，经过D-L拟合重建的NDVI曲线，在植被或作物的返青期、成熟期以及枯黄期等时间节点上与台站实测生育期较为接近，且与重建前高质量NDVI数据的误差较小，具有较高的保真能力和适应性。A-G拟合重建结果与D-L结果较为接近，只是根据D-L和A-G拟合算法重建的NDVI峰值结果略有差异。

3)在年均NDVI较低(NDVI<0.3)，且在生长季内植被季节生长变化特征不明显或曲线不具对称性的稀疏植被区，或以冬小麦为典型种植作物的区域中，经过S-G滤波重建的NDVI数据与台站实测资料较为接近，且与重建前高质量NDVI数据的误差较小，表现出相对较好的保真能力和适应性。

References：

[1]Sudipta S, Menas K. Interannual variability of vegetation over the Indian sub-continent and its relation to the different meteorological parameters. Remote Sensing of Environment, 2004, 90: 268-280.

[2]Li F, Zeng Y, Li X S,etal. Remote sensing based monitoring of interannual variations in vegetation activity in China from 1982 to 2009. Science China: Earth Sciences, 2014, 57: 1800-1806.

[3]Hou M T, Zhao H Y, Wang Z,etal. Vegetation responses to climate change by using the satellite-derived normalized difference vegetation index: A Review. Climatic and Environmental Research, 2013, 18(3): 353-364.

[4]Arnon K, Nurit A, Rachel T P,etal. Use of NDVI and land surface temperature for drought assessment: merits and limitations. Journal of Climate, 2010, 23: 618-633.

[5]Barbosa H A, Huete A R, Baethgen W E. A 20-year study of NDVI variability over the Northeast region of Brazil. Journal of Arid Environments, 2006, 67: 288-307.

[6]Guo N, Wang X P. Advances and developing opportunities in remote sensing of drought. Journal of Arid Meteorology, 2015, 33(1): 1-18.

[7]Piao S L, Wang X H, Ciais P. Changes in satellite-derived vegetation growth trend in temperate and boreal Eurasia from 1982 to 2006. Global Change Biology, 2011, 17(10): 3228-3239.

[8]Gu J, Li X, Huang C L. Research on the reconstructing of Time-series NDVI Data. Remote Sensing Technology and Application, 2006, 21(4): 391-395.

[9]Han P, Yao J, Li T H. Comparison of 3 NDVI datasets and the application at Yanhe basin, China. Journal of Basic Science and Engineering, 2014, 4: 661-674.

[10]Wang Z X, Liu C, Huete A. From AVHRR-NDVI to MODIS-EVI: Advances in vegetation index research. Acta Ecologica Sinica, 2003, 23(5): 979-987.

[11]Geng L Y, Ma M G. Advance in method comparison of reconstructing remote sensing time series data sets. Remote Sensing Technology and Application, 2014, 29(2): 362-368.

[12]Ma M G, Song Y, Wang X F,etal. Development status and application research of the time series remote sensing data products based on AVHRR、VEGETATION and MODIS. Remote Sensing Technology and Application, 2012, 27(5): 663-670.

[13]Du J Q, Shu J M, Wang Y H,etal. Comparison of GIMMS and MODIS normalized vegetation index composite data for Qinghai-Tibet Plateau. Chinese Journal of Applied Ecology, 2014, 25(2): 533-544.

[14]Michishita R, Zhenyu J, Jin C,etal. Empirical comparison of noise reduction techniques for NDVI time-series based on a new measure. ISPRS Journal of Photogrammetry and Remote Sensing, 2014, 91: 17-28.

[15]Beck P S, Goetz S J. Satellite observations of high northern latitude vegetation productivity changes between 1982 and 2008: Ecological variability and regional differences. Environmental Research Letters, 2011, 6: 5501-5511.

[16]Jeremy L W, David S G, Julia E A,etal. Long-term vegetation monitoring with NDVI in a diverse semi-arid setting, central New Mexico, USA. Journal of Arid Environments, 2004, 58(2): 249-272.

[17]Steven M D, Timothy J M, Frédéric B,etal. Intercalibration of vegetation indices from different sensor systems. Remote Sensing of Environment, 2003, 88(4): 412-422.

[18]Chen Y L, Long B J, Pan X B,etal. Differences between MODIS NDVI and AVHRR NDVI in monitoring grasslands change. Journal of Remote Sensing, 2011, 15(4): 831-845.

[19]Fensholt R, Proud S R. Evaluation of earth observation based global long term vegetation trends—Comparing GIMMS and MODIS global NDVI time series. Remote Sensing of Environment, 2012, 119(16): 131-147.

[20]Mao D H, Wang Z M, Luo L,etal． Integrating AVHRR and MODIS data to monitor NDVI changes and their relationships with climatic parameters in Northeast China. International Journal of Applied Earth Observation and Geoinformation, 2012, 18: 528-536.

[21]Geng L, Ma M, Wang X,etal. Comparison of eight techniques for reconstructing multi-satellite sensor time-series NDVI data sets in the heihe river basin, China. Remote Sensing, 2014, 6: 2024-2049.

[22]Li R, Zhang X, Liu B,etal. Review on methods of remote sensing time-series data reconstruction. Journal of Remote Sensing, 2009, 13(2): 335-341.

[23]Hird J N, McDermid G J. Noise reduction of NDVI time series: An empirical comparison of selected techniques. Remote Sensing of Environment, 2009, 113: 248-258.

[24]Li Z H. A Study on the Eco-environment Evolution of Yangtze River Delta Region Based on the Retrieval & Reconstruction of MODIS Time Series Datasets[D].Shanghai: East China Normal University, 2011.

[25]Jönsson P, Eklundh L. TIMESAT-a program for analyzing time-series of satellite sensor data. Computers & Geosciences, 2004, 30(8): 833-845.

[26]Holben B N. Characterististics of maximun-value composite images from temporal AVHRR data. International Journal of Remote Sensing, 1986, 7(11): 1417-1434.

[27]Cao Y F, Wang Z X, Deng F P. Fidelity performance of three filters for high quality NDVI time-series analysis. Remote Sensing Technology and Application, 2010, 1(1): 118-125.

[28]Song C Q, You S C, Ke L H,etal. Analysis on three NDVI time-series reconstruction methods and their applications in North Tibet. Journal of Geo-Information Science, 2011, 1(1): 133-143.

[29]Madden H H, Chem A. Comments on the Savitzky-Golay convolution method for least-squares-fit smoothing and differentiation of digital data. Analytical Chemistry, 1978, (9): 1383-1386.

[30]Moulin S, Kergoat L, Viovy N,etal. Global-scale assessment of vegetation phenology using NOAA/AVHRR satellite measurements. Journal of Climate, 2010, 10(6): 1154-1170.

[31]Liu H. Spring Phenology Model of Grassland Based on Soil Moisture and Air Temperature and Vegetation Reactions to Drought[D]. Beijing: Tsinghua University, 2012.

[32]Wang L, Ding J L. Vegetation index feature change and its influencing factors and spatial-temporal process analysis of desert grassland in the Ebinur Lake Nature Reserve, Xinjiang. Acta Prataculturae Sinica, 2015, 24(5): 4-11.

[33]Zhang X, Friedl M A, Schaaf C B,etal. Monitoring vegetation phenology using MODIS. Remote Sensing of Environment, 2003, 84(3): 471-475.

[3]侯美亭, 赵海燕, 王筝, 等. 基于卫星遥感的植被NDVI对气候变化响应的研究进展. 气候与环境研究, 2013, 18(3): 353-364.

[6]郭铌, 王小平. 遥感干旱应用技术进展及面临的技术问题与发展机遇. 干旱气象,2015, 33(1): 1-18.

[8]顾娟, 李新, 黄春林. NDVI 时间序列数据集重建方法述评. 遥感技术与应用, 2006, 21(4): 391-395.

[9]韩鹏, 姚娟, 李天宏. 3种不同数据源NDVI的比较分析及其在延河流域的应用研究. 应用基础与工程科学学报, 2014, 4: 661-674.

[10]王正兴, 刘闯, Huete A. 植被指数研究进展: 从AVHRR-NDVI到MODIS-EVI. 生态学报, 2003, 23(5): 979-987.

[11]耿丽英, 马明国．长时间序列NDVI数据重建方法比较研究进展．遥感技术与应用, 2014, 29(2): 362-368.

[12]马明国, 宋怡, 王旭峰, 等. AVHRR、VEGETATION和MODIS时间序列遥感数据产品现状与应用研究进展. 遥感技术与应用, 2012,27(5): 663-670.

[13]杜加强, 舒俭民, 王跃辉, 等. 青藏高原MODIS NDVI与GIMMS NDVI的对比. 应用生态学报, 2014, 25(2): 533-544.

[22]李儒, 张霞, 刘波, 等. 遥感时间序列数据滤波重建算法发展综述.遥感学报, 2009, 13(2): 335-341.

[24]黎治华. 基于MODIS反演重构时间序列数据的长江三角洲地区生态环境演变研究[D]. 上海： 华东师范大学, 2011.

[27]曹云锋, 王正兴, 邓芳萍. 3种滤波算法对NDVI高质量数据保真性研究. 遥感技术与应用, 2010, 1(1): 118-125.

[28]宋春桥, 游松财, 柯灵红, 等. 藏北地区三种时序NDVI重建方法与应用分析. 地球信息科学学报, 2011, 1(1): 133-143.

[31]刘慧. 基于土壤水分和气温的草地返青模型及植被干旱研究[D]. 北京： 清华大学, 2012.

[32]王璐, 丁建丽. 艾比湖保护区荒漠植被时空过程变化及其植被指数影响因素分析. 草业学报, 2015, 24(5): 4-11.

*Comparative studies of reconstruction methods for the long term NDVI dataset in the east of Northwest China

WANG Wei, GUO Ni*, SHA Sha, HU Die, WANG Xiao-Ping, LI Yao-Hui

KeyLaboratoryofAridClimaticChangeandReducingDisasterofGansuProvince,KeyOpenLaboratoryofAridChangeandDisasterReductionofCMA,InstituteofAridMeteorology,ChinaMeteorologicalAdministration,Lanzhou730020,China

A high-level time-series NDVI dataset is not only the basis for continuous monitoring of the land surface, but also an important tool for studying change related to climate and land use factors in terrestrial eco-systems. We reconstructed the noise component of the LTDR NDVI data for the east of Northwestern China where the ecosystem is fragile, using various time-series reconstruction methods. This paper use agrometeorological data and high-level NDVI data to evaluate the accuracy of different reconstruction methods. The results show that: 1) The vegetation or crop land cover is an important factor affecting fitted results of the various reconstruction methods. Each reconstruction method has a different noise reduction ability depending on differences in vegetation or crop growth characters; 2) The D-L reconstruction method has a better noise reduction ability and applicability in those areas of grassland, and woodland for which the annual average NDVI data is higher (NDVI≥0.3) and the NDVI curve has obvious seasonal changes; 3) The S-G reconstruction method has better fidelity ability and applicability in some areas of crop land in winter wheat and in areas of sparse vegetation for which annual average NDVI data are lower (NDVI<0.3) and where the NDVI curve have no obvious seasonal changes.

NDVI; time-series data reconstruction; vegetation remote sensing; AVHRR

10.11686/cyxb2015489http://cyxb.lzu.edu.cn

王玮, 郭铌, 沙莎, 胡蝶, 王小平, 李耀辉. 我国西北地区东部时间序列NDVI数据集重建方法比较研究. 草业学报, 2016, 25(8): 1-13.

WANG Wei, GUO Ni, SHA Sha, HU Die, WANG Xiao-Ping, LI Yao-Hui. Comparative studies of reconstruction methods for the long term NDVI dataset in the east of Northwest China. Acta Prataculturae Sinica, 2016, 25(8): 1-13.

2015-10-22；改回日期：2016-01-04

甘肃省气象局气象科研项目(GSMAMs2016-10)，中国博士后科学基金项目(2015M582734)和公益性行业(气象)科研专项(重大专项)(GYHY201506001-5)资助。

王玮(1985-)，男，甘肃兰州人，助理研究员。 E-mail: wangwei9969@163.com

Corresponding author. E-mail: guoni0531@126.com