基于FY3/MERSI数据的辽宁省植被指数重建及时空变化分析

冯锐 纪瑞鹏 武晋雯 于文颖 刘丹 陈妮娜 王莹 张玉书

(1.中国气象局沈阳大气环境研究所，辽宁 沈阳 110166； 2.辽宁省农业气象灾害重点实验室，辽宁 沈阳 110166；3.黑龙江省气象科学研究所，黑龙江 哈尔滨 150030； 4.辽宁省生态气象和卫星遥感中心，辽宁 沈阳 110166)

引言

辽宁省地处世界三大黑土区之一的东北黑土区南端，横跨黄、渤两海，属温带湿润、半湿润大陆性季风气候，是我国重要的粮食生产基地，也是气候变化敏感地区之一[1]。对植被状况变化的研究从一定程度上可以反映区域气候变化的趋势，可为未来粮食安全生产提供基础分析数据[2]。归一化植被指数(Normalized Difference Vegetation Index，NDVI)是开展植被动态变化、生长及覆盖状况研究的最佳指数之一，是监测生态环境变化和区域尺度地表改变的有效指标[3]。目前，MOD13Q1(MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m SIN Grid)产品在大范围植被动态监测、地表植被分类和物候信息提取等方面广泛应用[4-6]。随着卫星遥感资料的积累，虽然同系列的卫星传感器波段范围、通道宽度及分辨率相同，但受到卫星传感器状态、云、气溶胶、太阳高度角和地物二向反射等因素影响，卫星数据集依然会存在着很多噪声，给数据进一步应用与分析带来限制[7]。国内外研究人员针对不同下垫面状况、不同研究区域开展了大量的数据重建研究，比较常见的方法包括Savitzky-Golay滤波法[8]、傅里叶变换法(Fourier Transform，FT)[9]、均值迭代滤波(Mean-value Iteration Filter，MVI)[10]、时间序列谐波分析法[11]、Whittaker平滑法[12]、非对称高斯函数拟合法[13]和双Logistic函数拟合法[14]等，由于研究区域、植被覆盖类型和噪声程度不同，无法确定某种方法是否为普适的最优数据重建方法[15]。2008年5月27日中国第二代极轨气象卫星风云三号A星升空[16]，到目前为止已经发射了5颗系列卫星，数据积累年份超过12 a，其携带的中分辨率光谱成像仪(MERSI)具有5个(Ⅰ型)或6个(Ⅱ型)250 m分辨率通道，可以更好地实现陆地多光谱连续观测[17]，在大范围遥感监测中具有较高的实用价值[18]，但目前针对此传感器的植被指数数据开展长序列重建研究成果还未见报道。因此，本文基于中国风云系列FY3/MERSI数据，选择滤波和函数拟合方法对辽宁区域进行NDVI数据重建，并定量分析数据重建效果，以大田作物的生长阶段为节点，分析辽宁省植被指数时空变化，揭示区域植被变化特征，为辽宁省乃至东北地区的环境监测、气候变化对农业生产的影响评估、精细化农业气候区划以及粮食安全提供精准、可信的国产气象卫星数据。

1 资料与方法

1.1 研究区概况

辽宁省位于中国东北地区南部，在38.72°—43.43°N、118.88°—125.77°E之间，总面积14.59×104km2，东部和西部为山地丘陵，中部为平原(图1)，四季分明，雨热同季，年降水量为433—1100 mm，降水分布不均[19]。

图1 研究区位置示意图(a)和7类地物提取点位置图(b)Fig.1 Locations of the research area (a) and extraction points in seven types of ground features (b)

1.2 主要数据源

归一化植被指数计算采用2009—2020年FY3A/MERSI、FY3B/MERSI和FY3D/MERSI数据，数据分辨率250 m，由于FY3C/MERSI数据只有不到2 a的数据，因此，数据选择主要以FY3B/MERSI为主，FY3A/MERSI和FY3D/MERSI为补充，2009—2010年选用FY3A/MERSI数据，2011—2018年选用FY3B/MERSI数据，2019—2020年选用FY3D/MERSI数据。FY3A为上午星，FY3B和FY3D为下午星，FY3A-C携带MERSI(Ⅰ型)传感器，FY3D携带MERSI(Ⅱ型)传感器，本文计算植被指数所用通道为3通道和4通道，Ⅰ型和Ⅱ型的各参数基本一致(表1)。

表1 FY3/MERSI部分通道参数Table 1 Parameters in some channels of FY3/MERSI

FY3/MERSI数据来源于风云卫星遥感数据服务网(http://satellite.nsmc.org.cn/PortalSite/Default.aspx)，数据处理采用卫星监测分析遥感应用系统(Smart2.0)和ENVI软件，采用等经纬度投影，进行逐日NDVI计算，通过最大值合成法建立旬值植被指数原始数据集。植被指数原始数据集建立时，利用A-D星的光谱响应函数和地面光谱观测数据进行NDVI的相关分析，相关系数达到了0.99以上。因此，本文的NDVI数据重建中直接应用各卫星的NDVI计算值。

行政区划数据来自于中国气象局下发的1∶250000基础地理信息。

1.3 归一化植被指数计算

(1)

式(1)中，NDVI为归一化植被指数；ρnir为近红外波段反射率；ρred为红光波段反射率。在FY3/MERSI中，分别对应4通道和3通道。

1.4 数据重建方法

1.4.1 Savitzky-Golay滤波法(SG)

SG方法是基于最小二乘法原理的卷积滤波方法，基本公式为[20]：

(2)

式(2)中，Y*是拟合值；Y是原始值；系数j是指原NDVI数组的系数；滑动数组包含有2m+1个点；Ci是滑动窗口内第i个点的系数。

1.4.2 时间序列谐波分析法(HANTS)

HANTS方法是以傅立叶变换为基础的函数拟合法，表达式为[21]：

(3)

式(3)中，y是拟合值；A0是谐波余项；Aj是谐波振幅；ωj是谐波频率；i是y的时间标记；N是时间序列长度；θj是谐波初相位；m是谐波个数；j是谐波标记。

1.4.3 非对称高斯函数拟合法(AG)

AG方法是基于最小二乘法原理的拟合算法，采用不对称高斯函数进行非线性拟合，表达式为[22]:

(4)

式(4)中，g(x，a1,…,a5)是高斯函数；a1是确定最大或最小值位置的参数，a2、a3和a4、a5分别是确定右、左两边的拟合曲线宽度和陡峭度的参数，一般情况下，为避免高斯函数的匹配点处出现尖点x=a1，a3和a5设定值应大于2。

1.4.4 双Logistic函数拟合法(DL)

DL方法利用冬季植被指数、最大植被指数、曲线上升和下降两个拐点和拐点处植被指数变化速率进行拟合，表达式为[7,23]：

(5)

式(5)中，g(x,a1,…a4)是双Logistic函数；a1、a2和a3、a4分别是确定拟合曲线左、右两边的拐点位置和变化速率的参数。

1.4.5 四种方法的参数设置

SG、DL和AG数据重建是在TIMESAT软件中实现，而HANTS数据重建是在hants软件中实现，其参数设置对植被指数滤波后效果有很大的影响，参数的设定一般通过经验或者反复试验确定[24]，本文通过查阅相关的研究成果和反复试验，确定了各方法的参数[7,25-27](部分参数设置见表2)。对于SG、

表2 四种方法的参数设置Table 2 Parameter settings of four fitting methods

DL和AG方法，本文尖峰方法选择1，即中值滤波方法，与此方法相对应的尖峰值设置为2，尖峰值主要决定尖峰去除程度，低值将消除更多尖峰；季节设置为1，一年一季；滑动窗口是SG滤波需要设置的参数，这个值过小容易产生冗余数据，过大会遗漏细节信息，选择5；迭代次数的设置让时间序列数据在迭代的过程中更接近上部包络线。对于HANTS方法，频率数决定拟合偏差的大小，数值越大拟合偏差越小，但也会带来更多噪声，经过对比试验选取5效果较好；相位偏移代表数据起始月份，本系列数据从1月份开始，因此取值0。

1.5 数据重建效果评价方法

分别从时间序列曲线分析和误差统计分析两方面对数据重建效果进行评价，误差统计分析选择相关系数(R)和均方根误差(RMSE)2个参数进行比较和评价[28]。

2 结果分析

2.1 典型地物NDVI数据重建时间序列曲线比较

为评价辽宁地区SG、DL、AG、HANTS等方法数据重建效果和适用性，选择林地、芦苇湿地、水稻、玉米、大豆、城市和水体等不同地物，进行2009—2020年重建数据与原始数据时间序列曲线对比分析。从图2可以看出，四种方法进行数据重建时均表现出相对较好的去噪能力。对于高植被区(图2a至图2e)，在每年7—8月和每年11月中旬到次年2月底均出现较多噪声，对于城市和水体(图2f和图2g)，每年11月中旬到次年2月底出现噪声较多，这是由于每年11月中旬到次年2月底是辽宁省的冬季，噪声的形成除了云的影响，还有冰雪的因素，而7—8月是植被生长茂盛阶段，也恰好是辽宁省的夏季，噪声主要来源于云的影响。

图2a中林地的原始NDVI曲线在高值区波动较大，四种方法对高值区和低值区的突降点均做了平滑处理，有所不同的是，AG和DL方法拟合后的曲线对高值区和低值区都做了较好的平滑，SG和HANTS方法受到噪声影响，更多的保留了突降点的变化特征，曲线波形变化剧烈，甚至在某些年份出现了双峰现象，尤其是在2011年、2013年、2015年和2017年7—8月出现的低值噪声滤波时表现不佳，以2017年8月上旬为例，林地原始植被指数为0.09，存在着严重偏低情况，滤波后SG方法NDVI为0.57，DL方法为0.71，AG方法为0.65，HANTS方法为0.58，结合上一旬和下一旬植被指数，分别为0.68和0.72，DL方法滤波后的结果更接近真实生长过程，SG和HANTS存在着低估现象。

图2 2009—2020年辽宁地区不同地物4种重建数据与林地(a)、湿地(b)、水稻(c)、玉米(d)、大豆(e)、城市(f)和水体(g)原始数据对比曲线Fig.2 Comparisons between reconstruction data and original data of woodland (a),wetland (b),rice (c),corn (d),soybean (e),urban (f),water body (g) and their reconstruction data from 2009 to 2020

图2b是湿地的NDVI变化曲线，辽宁省的芦苇湿地位于辽河三角洲南端，是亚洲第一大芦苇湿地。四种方法重建数据的湿地数据与林地数据比较相似，除了HANTS方法在高值区有了较好的平滑处理外，SG方法受到噪声影响，仍出现比较多的扰动，AG和DL方法对高、低值区都做了较好的平滑，DL方法的峰值更接近原始峰值。

图2f和图2g中，城市和水体的原始NDVI值都比较低，基本集中在-0.20～0.25，除HANTS方法在低值区受噪声影响，曲线波动比较大以外，SG、AG和DL三种方法均表现不错，做了较好的平滑处理，更接近实际的情况。

图2c、图2d和图2e分别为辽宁省主要的三大粮食作物水稻、玉米和大豆的NDVI变化曲线，辽宁省的农作物为一年一熟制，因此年际间的曲线变化均为单峰形式，AG和DL方法平滑效果较好，低值区表现优秀。在高值区，玉米和大豆两种方法均存在着低估的情况，水稻AG方法也存在着低估情况，但DL方法没有明显低估现象。SG和HANTS方法在农田高值区中有很好的表现，但在低值区受突变点影响，扰动明显。

2.2 典型地物数据重建时间序列数据定量分析

为定量评价2009—2020年辽宁地区SG、DL、AG和HANTS方法的适用性和保真能力，在林地、湿地、水稻、玉米、大豆、城市和水体等不同地物共选取采样点498个，分别计算相关系数和均方根误差(图3)。从图3可以看出，在林地、湿地、大豆、水稻和玉米这五类高植被覆盖区和季节性作物区，SG方法重建后的NDVI与初始值的相关系数均高于其他三种方法，达到了0.93以上。同时，均方根误差低于其他三种方法，尤其是大豆，RMSE小于0.05。在城市和水体低植被指数区，HANTS方法重建后的NDVI与初始值相关系数最高，为0.87，其余三种方法重建后的NDVI与初始值相关系数差别不大，为0.79—0.82，四种方法的均方根误差的差别不大，在0.06左右。因此，从数据保真性来看，SG方法效果更优。

图3 2009—2020年辽宁地区不同地物4种重建数据与原始数据的相关系数和RMSEFig.3 The correlation coefficients and RMSE between reconstruction data and original data based on the four methods in various ground features from 2009 to 2020

2.3 辽宁省植被指数空间分布

对比时间序列曲线分析和基于统计的定量分析结果可以发现，不同分析四种数据重建方法的适用性并不一致，曲线分析认为，不论在高植被覆盖区还是在低植被覆盖区DL方法均有较好的表现，统计定量分析认为，SG方法在高植被覆盖区更适合，HANTS方法在低植被覆盖区表现更为突出。造成这一结果的原因，与SG方法更多的保留了原始NDVI数据的细节，而DL方法对高值区和低值区都做了较好的平滑有关。综合考虑曲线和定量分析结果，并且数据应用时更多的关注高植被覆盖区和作物种植区，因此，选取数据保真性更高的SG方法对辽宁省植被指数数据集进行数据重建。

将2009—2020年重建后的逐旬植被指数进行年内最大值计算，得到2009—2020年逐年植被指数，再将计算结果平均，得到辽宁省植被指数空间分布(图4)，从图4可以看出，植被指数的高低分布与下垫面植被类型相符合，辽宁东部林地、中部水稻种植区、中北部旱田种植区及西部少量林地的植被指数在0.70以上，尤其是东部山区林地达到0.75以上；辽西灌木林区、盘锦芦苇湿地及全省大部分旱田种植区，植被指数均在0.60—0.70，仅有小部分辽宁西部和辽宁南部的旱田植被指数在0.50—0.60；水体、城市、滩涂等，植被指数均低于0.30。从植被指数分布来看，高值集中，0.55以上的像元占到了92%，占比最高的数值是0.65—0.70，为26%。

图4 2009—2020年辽宁省植被指数空间分布Fig.4 Spatial distributions of vegetation indices in Liaoning province from 2009 to 2020

2.4 辽宁省植被指数年际变化

基于上文2009—2020年逐年植被指数，将辽宁省范围内的植被指数和林地、城市、大豆、湿地、水稻、玉米及水体七类地物提取点植被指数平均，进行辽宁省植被指数年际分析(图5)，从图5可以看出，2009—2020年辽宁省NDVI年均值存在波动，最大值出现在2012年，为0.77，最低值出现在2014年，为0.59；12 a中，有5 a的相对变化百分率为负值，负增长最大值出现在2012—2013年，为-18.0%，正值年份有6 a，增长率最大为13.3%，出现在2015—2016年。

从图5b可以看出，2009—2020年，辽宁省林地植被指数为0.67—0.89，2009年出现最大值，为0.89，2019年出现最小值，为0.67；湿地和水稻植被指数数值接近，低于林地的植被指数，大部分值集中在0.65左右；旱田作物(玉米、大豆)的植被指数各年均比水稻植被指数值低，且两者的植被指数非常接近，受干旱年的影响植被指数变化稍大，在2014年和2017年出现低值[29-31]；水体和城市植被指数在12 a间变化相对较小，两者之间的数值差别也很小，一般情况下，在0.07以内，且远低于其他五类地物，为0.14—0.33。

图5 2009—2020年辽宁省平均植被指数(a)和不同地物均值(b)年际变化Fig.5 Inter-annual variations of the average vegetation indices in Liaoning Province (a) and mean values of various ground features (b) from 2009 to 2020

2.5 辽宁省大田作物植被指数年内变化

将2009—2020年辽宁省玉米、大豆和水稻植被指数进行平均，得到三种作物年内变化曲线(图6)。从图6可以看出，辽宁省一般在5—9月是作物生长季，植被的变化也主要在这个时间段体现。在生长季中，三种作物NDVI均呈单峰分布，玉米和大豆的生长曲线基本重合，植被指数均在5月中旬开始进入持续增长阶段，在8月上旬达到最大值，分别是0.57和0.53；水稻的生长曲线与玉米、大豆有很大不同，水稻在移栽期(5月下旬)之后NDVI快速增大，出现生长陡坡，在8月上旬达到最大值0.67。从作物植被指数曲线变化来看，利用中分辨率卫星数据进行旱田作物分类的可能性很低，但可利用水稻快速生长期的植被指数变率与移栽期植被指数相结合开展水稻面积提取。

图6 2009—2020年辽宁省玉米、大豆和水稻年内旬值植被指数变化Fig.6 Annual vegetation index variations in ten-day values of corn,soybean and rice from 2009 to 2020

3 结论与讨论

(1)SG、DL、AG和HANTS四种方法重建后的时间序列曲线与原始曲线相比，均表现出相对较好的去噪能力。SG方法对噪声比较敏感，会保留更多原始曲线的细节信息。HANTS方法在低值区受噪声影响大。AG和DL方法平滑效果较好，但存在峰值低估的现象。从定量角度看，在高植被覆盖区和季节性作物区，SG方法相关系数最高(>0.93)、均方根误差最低(<0.1)。综合考虑曲线和定量分析结果，选取SG方法进行辽宁省植被指数数据集数据重建，为风云三数据进行植被动态变化分析提供科学、准确的基础数据。

(2)重建后的辽宁省植被指数空间分布与下垫面植被类型相符合，SG方法在辽宁省植被重建中具有较好的效果。从植被指数年际变化来看，以2014年为时间节点，存在着先下降后上升的现象，旱田作物(玉米、大豆)的植被指数受干旱年的影响植被指数变化较大，重建后的植被指数数据集可为气候变化影响评估、灾害监测评估提供基础数据。

(3)水稻、玉米和大豆三种作物的植被指数年内均呈单峰分布，水稻存在着6月上中旬生长陡坡期等作物生长曲线变化，可为中分辨率卫星数据开展作物分类提供定量信息依据。

(4)目前，国内外学者已发展多种时间序列数据重建方法[15]，本研究选取了一种滤波方法(SG)和三种函数拟合方法(AG、DL和HANTS)，从拟合效果来看，四种方法均表现出相对较好的去噪能力，AG方法和DL方法在进行植被指数长时间序列数据重建时效果相近，这与众多研究结果一致[7,27,32-34]，这两种方法在植被的高低值区均有明显的平滑表现，这与AG方法和DL方法在进行数据拟合时，考虑植被生长变化的整体特征有关，尤其是在辽宁省，作物是一年一熟的单峰曲线，这两种方法均能对突然出现的噪声进行平滑，但也会对一些作物生长过程中的局部变化造成忽略[7]，比如水稻在返青—分蘖期的生长陡坡。从曲线分析来看，DL和AG方法在耕地的重建中均有较好的表现，与王乾坤等[35]和卫炜等[36]研究结果一致，但是从相关系数和均方根误差的数据来看，这两种方法的表现并不出色。HANTS方法在低植被覆盖区和水体中表现突出，相关系数远高于其他三种方法，这与邹明亮等[37]针对疏勒河流域的数据重建时结论一致。SG方法对噪声敏感，在数据波动明显时，重建结果会存在波形扰动，但同时对局部细节也更加敏感，数据的保真性更好[7,35,38]。每种数据重建方法均有其适用性，因此，应综合考虑植被类型、数据应用目的选择数据重建方法。

从辽宁省植被指数空间分布来看，数值高低分布与东部山区、西部丘陵和中部平原的地形相吻合。从时间变化来看，年均归一化植被指数存在着波动，且低值年(2014年、2017年和2019年)均与辽宁的干旱年相吻合，国志兴等[39]分析1982—2003年东北地区植被指数变化趋势认为森林、农田的植被指数均呈下降趋势，本研究中，植被指数的趋势线以2014年为分界点，呈现先下降后上升的趋势，但由于时间序列较短，数据的变化是否能支撑此论点还有待于进一步研究。辽宁主要大田作物的植被指数年内变化均呈单峰曲线[40]，在8月上旬达到峰值，这一结论与国志兴等[39]一致。