大区域TM影像归一化植被指数季相归一化处理方法

陈贤干，汪小钦，陈芸芝

（1.福州大学空间数据挖掘与信息共享教育部重点实验室，福州350002；2.福州大学空间信息工程研究中心，福州350002）

大区域TM影像归一化植被指数季相归一化处理方法

陈贤干1，2，汪小钦1，2，陈芸芝1，2

（1.福州大学空间数据挖掘与信息共享教育部重点实验室，福州350002；2.福州大学空间信息工程研究中心，福州350002）

针对高分辨率遥感监测植被时大范围获取影像的季相往往不一致，影响了归一化植被指数的连续性与可比性问题，借助MODIS归一化植被指数产品，开展TM归一化植被指数季相归一化方法研究。以福建省为例，分析不同季相MODIS归一化植被指数之间的一次项、二次项最小二乘关系及差值关系，对比分析转换结果。结果表明，二次项关系优于其他两种效果，较好地消除归一化植被指数时相不一致问题，获得了较为满意的效果。利用MODIS归一化植被指数构建二次项，完成福建省TM归一化植被指数季相归一化，消除了不同景影像间的季相不一致，使在大区域范围内归一化植被指数具有空间上的可比性和一致性。

MODIS NDVI；归一化方法；TM NDVI；季相差异；最小二乘关系

0 引 言

植被指数主要通过红光波段和近红外波段反射率的组合计算来设计并反映植被状态，常被用来表征地表植被覆盖、评价植被的生长及分布特征，并作为一种有效的遥感指标应用于植被监测、农情监测、土地覆盖变化等领域。目前遥感领域中提出的植被指数有100多种［1］，如归一化植被指数（NDVI）、比值植被指数（RVI）、差值植被指数（DVI）、垂直植被指数（PVI）、土壤修正植被指数（SAVI）等［2－5］，其中NDVI是最常用的。

高时间分辨率的遥感数据能以较短的重访周期进行大范围的重复访问，能够较好地反映植被的季相变化，但空间分辨率较低。如MODIS的NDVI产品能提供一致的、空间分布的、多时相的全球植被状态参数，在区域与全球植被动态监测中得到广泛应用［6－14］。然而其250m～1km的空间分辨率却不足以反映地表植被状况空间分布中高分辨率的差异，无法满足中高分辨率的植被监测应用需求。在应用中高分辨率遥感进行植被监测时，中高空间分辨率的遥感数据通常只覆盖较小的空间范围，在多云多雨的广大南方地区，同一区域或者较大区域所获取影像季相往往不一致，影响了植被指数的连续性与可比性，如何消除较大区域范围内植被指数的季相影响，是中高分辨率遥感应用的难题之一。

如何综合利用中高空间分辨率数据和高重访周期数据，生成同时具有高时间和高空间分辨率的遥感数据对于遥感的深入应用意义重大。大多数研究人员采用时空融合的方法，如Lorenzo等［15］将多时相低分辨率NDVI与基于高分辨率遥感数据的植被分类图相结合，进行亚像元分解，得到了高时空分辨率的NDVI数据；Gao等［16］提出时空适应性反射率融合模型（STARFM），获取MODIS数据相应时间的模拟 Landsat数据；蒙继华等［17］改进STARFM算法，发展出不同时空分辨率植被指数融合模型STAVFM，将MODIS NDVI时间变化信息与ETM NDVI空间差异信息的有机结合，从而构建高时空分辨率NDVI数据。但是，这些融合方法大多关注于算法的实现，应用区域较小，其对大范围的应用效果和效率有待于验证。针对不同传感器的植被指数，研究主要侧重于分析不同传感器由于光谱响应函数差异引起的植被指数间存在差异及相互转换［18－20］，没有考虑季相的差异。

在中高分辨率遥感数据中，TM数据应用最广，且时间序列长。本文以福建省为例，针对不同季相的TM数据，借助MODIS的NDVI产品，开展TM影像NDVI季相归一化方法研究。

1 研究区域与方法

1.1 研究区概况与数据预处理

（1）研究区概况

福建省地处我国东南沿海，介于23°33′N～28° 20′N、115°50′E～120°40′E之间，东与我国台湾省隔海相望，东北与浙江毗邻，西北横贯武夷山脉与江西相交，西南与广东相连，如图1所示。位于亚热带地区，属典型的亚热带季风气候，夏季高温多雨，冬季温和少雨，平均气温15.3℃～21.9℃，平均降雨量930mm～1843mm。全省土地面积121400km2，其中山地、丘陵占全省总面积的80%以上，地势总体上西北高、东南低。研究区主要的土地覆盖／利用类型为林地、耕地、裸地、居民地和水体，以林地为主，约占全省土地面积的2／3。

（2）Landsat5TM数据预处理

完整覆盖福建省TM影像需要13景，选择2010年左右无云和少云的影像，影像的获取时间如下：轨道号118／41、118／42的2009年4月28日TM影像；轨道号119／41、119／42、119／43的2009年6月6日TM影像；轨道号120／41、120／42、120／44的2009年10月3日TM影像；轨道号120／40、120／42、120／43的2009年12月6日TM影像；轨道号121／41、121／42、121／43的2010年1月14日TM影像。其中轨道号120／42的影像有2景，2009年12月6日的影像作为参与季相归一化的数据，2009年10月3日影像作为方法验证比较的参考数据。

利用遥感图像处理软件，对Landsat5TM遥感影像进行几何校正、大气辐射校正、影像裁剪等预处理，计算NDVI（文中用TM NDVI表示）。由于影像获取的季相各不相同，无法进行对比分析，拟将不同时相TM NDVI数据归一化到2009年10月中上旬。

（3）MODIS NDVI数据预处理

MODIS NDVI数据来源于美国国家航空航天局NASA免费提供的覆盖福建省的MOD13Q1级植被指数产品，时间范围为2009年～2010年，空间分辨率250m，时间分辨率为16d。采用MODIS网站上提供的MRT（MODIS Reprojection Tool）重投影工具对研究区数据进行格式转换和重投影，将影像的原始投影转换为与TM数据一致的投影，同时完成图像裁剪和重采样。

图1 福建省地理位置图

1.2 方法

首先对TM数据和MODIS NDVI数据进行预处理获得相应的NDVI数据，分析TM NDVI与MODIS NDVI时间变化趋势相似性，根据不同季相MODIS NDVI之间不同关系，将T0时刻TM NDVI分别模拟为Tt时刻的TM NDVI数据，并与参考数据进行对比分析，选择适宜的模拟方法进行TM NDVI数据季相归一化。所采用的总体技术路线如图2所示。

（1）TM影像植被指数计算

归一化植被指数［21］（Normalized Difference Vegetation Index，NDVI）是植被生长状态及植被覆盖度的最佳指示因子，在植被遥感中应用最为广泛。NDVI经比值处理可以部分消除与太阳高度角、卫星观测角、地形、云、阴影和大气条件有关的辐照度条件变化（大气程辐射）等的影响。其计算公式为：

（2）归一化方法

植被指数季相归一化，是基于同一时相的植被长势在空间地域上和同一空间地域的植被在不同时相上存在差异，为消除植被物候引起的季相差异而提出来的归一化方法，使其具备时间和空间上可比性［20］。季相归一化基本原理的前提是同一时相的MODIS NDVI与TM NDVI满足空间分布和动态变化的一致性［22－24］，根据选定 T0、Tt时刻的MODIS NDVI数据，利用它们之间的相互关系来模拟回归方程，将T0时刻的TM NDVI数据归一化到Tt时刻的TM NDVI数据。

本文主要基于最小二乘多项式（一次项和二次项）、差值关系开展季相的归一化。

由于多项式拟合方法简单，计算速度快，且能拟合所有像元，考虑到研究区范围和NDVI变化特征，因此采用最小二乘多项式回归方程来完成影像归一化：

其中，x是T0时相的影像；y是Tt时相影像；c0、c1、c2、c3…cn为拟合系数。利用T0、Tt时相的MODIS NDVI建立最小二乘多项式关系，求出相关的拟合系数，然后应用于TM影像。本文采用一次项式（3）、二次项式（4）方法。

式（3）、式（4）中，a、b、c为拟合系数。

差值关系是指不同时相数据之间存在差值增量，根据不同时相参考数据之间的差值关系来模拟实验数据，实现季相的归一化：

4.耦合协调评价表明，2001—2014年京津冀城市群土地综合承载力与区域经济发展的耦合协调度均较低，处于从中度失调向轻度失调、再向勉强协调、然后向中度协调过渡状态；京津冀城市群土地综合承载力与区域经济发展之间还存在着互动不足、协调不力等问题，有必要探索土地综合承载力约束下京津冀城市群经济发展的持续增效策略。

其中，x为已知T0时相的TM NDVI，y为待归一化的Tt时相的影像；Δm为Tt、T0时刻的MODIS NDVI数据之间的差值。

2 结果与分析

2.1 TM NDVI与MODIS NDVI时间变化趋势分析

TM NDVI是否与MODIS NDVI具有相似的时间变化趋势，是MODIS NDVI的变化方程能否应用于TM NDVI归一化的关键。利用轨道号120／42获取的2009年12月6日（T0）和10月3日（Tt）TM数据及相应时期的MODIS NDVI开展两者时间变化趋势分析。

分别将T0、Tt时刻的TM NDVI和MODIS NDVI的像元值投影到二维特征空间上，进行回归分析，求出转换方程（图3），统计它们的直方图（图4），并且利用它们之间的转换方程分别对T0时刻TM NDVI实验数据进行回归模拟，然后分别用各模拟数据减去参考数据，获得各模拟的差值直方图（图5）。

图3中，从转换方程来看，T0、Tt时刻TM NDVI一次项和二次项转换方程及MODIS NDVI的转换方程决定系数在0.60～0.71之间，相关性较高。TM NDVI一次项转换关系与MODIS NDVI一次项的差异小于0.03；在T0时刻TM NDVI的主要分布区域（0.15，0.55），二次项的差异也小于绝对值0.03。因此TM NDVI与MODIS NDVI的转换相关差异较小，两者具有较为一致的时间变化趋势。

图3 T0时刻NDVI（X轴）与Tt时刻NDVI（Y轴）散点图及转换方程

图4 不同时相NDVI直方图分布

图5 差值直方图

从图4直方图得出，12月份NDVI值小于10月份，MODIS NDVI总体上大于TM NDVI，在中高植被覆盖区域，T0和Tt时刻TM NDVI的变化趋势与MODIS NDVI的变化趋势非常相似。从差值直方图（图5）看出，T0、Tt时刻TM NDVI之间模拟结果与MODIS NDVI之间模拟结果的差值误差均在x＝0两侧分布，且变化趋势一致，说明满足其转换精度要求。综上分析，可知TM NDVI与MODIS NDVI随时间变化趋势具有很高的相似性，利用MODIS NDVI的变化趋势来进行TM NDVI的季相归一化是可行的。

2.2 归一化方法比较分析

利用真实数据，从误差分析、模拟数据统计分析等方面对不同归一化方法进行比较分析。

（1）误差分析

根据归一化方法原理，将 T0、Tt时刻的MODIS NDVI的像元值投影到二维特征空间上，然后对二者的关系进行回归分析，求出二者NDVI之间的转换方程（图3（b））。利用二者一次项、二次项以及差值关系的转换方程对T0时刻TM NDVI实验数据进行回归模拟，利用图3（b）关系方程以及公式（5）方法将实验数据归一化为Tt时刻TM NDVI模拟数据。

本文利用模拟数据与参考数据进行比较，从而对转换关系方程进行验证。在影像上随机选取67026个像元作为验证样本，然后分别用各模拟数据减去参考数据，获得各模拟的差值数据，并将差值数据与参考对应数据投影到二维特征空间上（图6），差值数据的值越趋近于零就越接近真实值。再将模拟数据与参考数据的验证样本像元值投影到二维特征空间（图7），通过计算均方根差（RMSE）来评价转换方程的模拟精度［25］。RMSE计算公式如下：

其中，Z′k为模拟数据的NDVI值，Zk为参考数据的NDVI值，N为验证样本个数。

图6中，一次项、二次项的误差数据基本上分布在y＝0两侧，点更加的聚集；而差值关系的误差值大部分分布在y＝0下方。二次项的斜率较一次项平缓，更接近于y＝0直线，表明二次项拟合的效果更接近真实值。图7中，一次项、二次项模拟数据的RMSE很小，只有0.0673和0.0671，表明拟合效果好，二者中又以二次项模拟的RMSE稍小些，而且二次项模拟值与参考值的散点基本沿y＝x线两侧分布，较一次项模拟、差值模拟精度更高。特别是在NDVI大于0.3的区域，模拟值与参考值更为接近，对于季相变化敏感的植被覆盖区域有较好的适宜性。

图6 各模拟的差值数据与参考数据散点图

图7 参考数据与模拟数据NDVI散点图

（2）TM NDVI统计分析

统计实验数据、参考数据与模拟数据，数据结果如表1所示。

从表1可见，一次项、二次项模拟数据的统计特征值与参考数据对应的统计特征值接近，差值关系的统计特征值与实验数据对应的统计特征值接近，这说明一次项、二次项模拟数据的结果比差值模拟结果理想，主要原因是差值关系的运算考虑到不同时相间隔时间要短，以及不同空间分辨率尺度转换。

表1 轨道号120／42 TM NDVI数据统计特征比较

从一次项、二次项模拟数据统计特征值得出，获取时间为12月实验数据其NDVI小于10月获取的值是符合季相变化特征的。经过归一化处理后，3种归一化方法的均值、最大值均有不同程度的提高，一次项和二次项的结果统计特征比较接近，其中二次项的模拟结果与参考影像最为相似。利用二次项关系进行季相归一化处理，基本可以满足中高分辨率大范围NDVI数据一致性的应用要求。

（3）典型地物NDVI比较

一次项与二次项模拟结果统计上比较接近，但他们对不同典型地物的适应程度如何还需要进一步分析。在模拟和参考影像上选取林地、耕地、裸地、不透水面、水体等不同地物的NDVI值作为样本，计算均值和RMSE，如表2所示。从表中可以看出：

林地：一次项、二次项归一化均值很接近参考均值0.598，RMSE小于0.05，二者林地模拟整体精度很高。

耕地：一次项、二次项的模拟均值较参考均值相差0.1左右，且两者RMSE相似，由于不同季相耕地的植被覆盖受季节种植作物变化影响很大，且空间尺度上也未能准确地反映其细节，所以在耕地模拟上未能达到较理想的效果。

非植被区域：裸地二次项比一次项的模拟均值接近参考均值，且RMSE也较小。而在不透水面和水体区域，从模拟的均值来看，二者归一化的结果都并不理想，但是二次项RMSE小于一次项，RMSE分别为0.1433、0.2138小于0.1679，0.2465，二次项模拟精度更高。

理论上，裸地、不透水面、水体在同一区域不同季相上的NDVI值变化不大，从这三类地物的模拟结果来看，二次项模拟结果更好；在同一区域林地占了80%以上植被覆盖类型，其模拟精度的高低决定了季相归一化效果的好坏；而耕地植被覆盖变化情况在同一区域随季节变化影响大，给季相归一化带来很大的难度。综上分析，总体上二次项模拟效果更好，本文利用二次项模拟方法对福建省进行大区域季相归一化处理。

表2 轨道号120／42不同地物TM NDVI均值及其RMSE比较

2.3 福建省TM NDVI季相归一化

从以上分析可知，二次项的拟合效果总体上优于一次项和差值方法。利用各轨道TM影像获取时相和2009年10月上旬的MODIS NDVI，计算二次拟合方程；利用这些二次拟合方程，完成福建省全省TM NDVI的季相归一化（图8（c））。图8（a）为福建省原始TM NDVI，图8（b）为福建省镶嵌影像生成的TM NDVI，图8（d）为福建省2009年10月中上旬的MODIS NDVI，主要用于对比分析，其中图8（a）和图8（c）的拼接方法相同，基本不改变影像原始的NDVI值，图8（b）为镶嵌影像生成的NDVI，为了使同一区域目视效果好，使其色调较一致，其镶嵌过程方法不同与图8（a）和图8（c）的拼接方法，从而改变影像原始的NDVI值。图9为图8中6个局部区域的放大图像，其中第1列为TM影像、第2列为原始TM NDVI、第3列为镶嵌影像生成的TM NDVI、第4列为季相归一化后TM NDVI、第5列为2009年10月中上旬的MODIS NDVI。

对比图8（a）、图8（c）和图9可以看出，当原始TM影像的季相差异较大时，归一化前的NDVI色调差异明显，相邻区域的数据出现明显的分块现象，影响了大范围区域数据的可比性；从图8（b）、图8（c）和图9得出，整体上，镶嵌影像生成的NDVI较图8（a）消除了一定的色调差异，与归一化后NDVI相比，色调过渡稍差，尤其是图8（b）中沿海地区NDVI值明显偏高（局部图9（f）可看出），这是由于影像镶嵌时改变了影像原始的NDVI值所造成；从图8（c）、图8（d）和图9中得出，归一化后NDVI色调过渡自然，与MODIS NDVI具有较为一致的空间变化特征，且归一化后TM NDVI更能体现原始影像的细节差异和植被纹理特征，空间分辨率较MODIS NDVI有较大的提高，能够满足较大比例尺的植被监测应用。

综上所述，借助MODIS NDVI数据，基于二次项最小二乘关系能够快速实现TM NDVI的季相归一化，减少不同时相NDVI的差异性，使得大区域范围内NDVI具有空间上的可比性和连续性。

图8 福建省NDVI对比分析

图9 福建省NDVI局部区域对比

4 结束语

TM NDVI的季相变化与MODIS NDVI具有较高的一致性，借助MODIS NDVI不同时相的相关方程，可以应用于TM影像NDVI的季相归一化。一次项、二次项模拟结果好于简单的差值关系；从不同方法归一化后全局和典型地类统计分析可以看出，二次项模拟效果最优。基于二次项的福建省TM NDVI影像季相归一化处理，很好地消除了不同季相NDVI的空间差异性，且较MODIS NDVI较好地反映了原始影像细节和纹理特征，使得大区域范围内NDVI不仅具有较高空间分辨率，而且有更好的空间可比性、连续性和一致性。如何改进归一化方法，对归一化结果进行更为客观的评价与分析需要进一步的深入研究。

［1］傅银贞，汪小钦.基于北京一号CCD数据的植被指数特性分析［J］.测绘科学，2010，35（6）：35－38.

［2］ROUSE J W，HAAS R H，SCHELL J A，et al.Monitoring vegetation systems in the Great Plains with ERTS［R］.Third Erts Symposium，NASA SP－351，1973，1：309－317.

［3］JORDAN C F.Derivation of leaf area index from quality of light on the forest floor［J］.Ecology，1969，50：663－666.

［4］RICHARDSON A J，WIEGAND C L.Distinguishing vegetation from soil background information［J］.Photogrammetric Engineering and Remote Sensing，1977，43：1541－1552.

［5］HUETE A R.A soil adjusted vegetation index（SAVI）［J］.Remote Sensing of Environment，1988，25：295－309.

［6］JUSTICE D H，SALOMONSON V，PRIVETTE J，et al.The Moderate Resolution Imaging Spectroradiometer（MODIS）：Land remote sensing for global change research［J］.IEEE Transactions on Geoscience and Remote Sensing，1998，36（4）：1228－1249.

［7］RUNNING S W，JUSTICE C，SALOMONSON V，et al.Terrestrial remote sensing science and algorithms planned for EOS／MODIS［J］.International Journal of Remote Sensing，1994，15（17）：3587－3620.

［8］陈燕丽，龙步菊，潘学标，等.MODIS NDVI和AVHRR NDVI对草原植被变化监测差异［J］.遥感学报，2011，15（4）：831－845.

［9］张峰，吴炳方，刘成林，等.利用时序植被指数监测作物物候的方法研究［J］.农业工程学报，2004，20（1）：155－159.

［10］王红说，黄敬峰.基于MODIS NDVI时间序列的植被覆盖变化特征研究［J］.浙江大学学报，2009，35（1）：105－110.

［11］何月，樊高峰，张小伟，等.浙江省植被物候变化及其对气候变化的响应［J］.自然资源学报，2013，28（2）：220－233.

［12］朱明明，侯西勇，毋亭.基于MODIS NDVI的中国东部植被覆盖时空变化特征［J］.遥感技术与应用，2013，28（6）：1027－1032.

［13］刘宪锋，任志远，林志慧，等.2000—2011年三江源区植被覆盖时空变化特征［J］.地理学报，2013，68（7）：897－908.

［14］穆少杰，李建龙，陈奕兆，等.2001—2010年内蒙古植被覆盖度时空变化特征［J］.地理学报，2012，67（9）：1255－1268.

［15］BUSETTO L，MERONI M，COLOMBO R.Combining medium and coarse spatial resolution satellite data to improve the estimation of subpixel NDVI time series［J］.Remote Sensing of Environment，2008，112（1）：118－131.

［16］GAO F，MASEK J，SCHWALLER M，et al.On the blending of the Landsat and MODIS surface reflectance：predicting daily Landsat surface reflectance［J］.IEEE Transactions on Geosciences and Remote Sensing，2006，44（8）：2207－2218.

［17］蒙继华，吴炳方，杜鑫，等.高时空分辨率NDVI数据集构建方法［J］.遥感学报，2011，15（1）：44－59.

［18］郭丽峰，高小红，亢健，等.伪不变特征法在遥感影像归一化处理中的应用［J］.遥感技术与应用，2009，24（5）：588－595.

［19］甘文霞，沈焕锋，张良培，等.采用6S模型的多时相MODIS植被指数NDVI归一化方法［J］.武汉大学学报（信息科学版），2014，39（3）：300－304.

［20］裴志远，杨邦杰.多时相归一化植被指数NDVI的时空特征提取与作物长势模型设计［J］.农业工程学报，2000，16（5）：20－22.

［21］郭铌.植被指数及其研究进展［J］.干旱气象，2003，21（4）：71－75.

［22］陈拉，黄敬峰，王秀珍.不同传感器的模拟植被指数对水稻叶面积指数的估测精度和敏感性分析［J］.遥感学报，2008，12（1）：143－151.

［23］徐涵秋，张铁军.ASTER与Landsat ETM＋植被指数的交互比较［J］.光谱学与光谱分析，2011，31（7）：1902－1907.

［24］杜加强，舒俭民，王跃辉，等.青藏高原MODIS NDVI与GIMMS NDVI的对比［J］.应用生态学报，2014，25（2）：533－544.

［25］赵凯，徐剑波，赵之重，等.HJ－1A／B CCD与Landsat TM／ETM＋植被指数的交互比较［J］.遥感技术与应用，2013，28（4）：674－680.

Seasonal Normalization Processing Methods of TM NDVI Images in Large Areas

CHEN Xian－gan1，2，WANG Xiao－qin1，2，CHEN Yun－zhi1，2
（1.Key Laboratory of Spatial Data Mining ＆Information Sharing of Ministry of Education，Fuzhou350002；2.Spatial Information Research Center of Fujian Province，Fuzhou University，Fuzhou350002）

Moderate resolution remote sensing data is widely used in vegetation monitoring.However，seasonal aspects of large areas is easily inconsistent in acquired images，which impair the continuity and comparability of normalized difference of vegetation index.MODIS data has the characteristics of high temporal resolution and wide coverage.The correlation analysis of TM NDVI and MODIS NDVI in different seasons shows that the variation of TM NDVI is quite similar to that of MODIS NDVI in the vegetation covered area.Therefore，with the MODIS NDVI products data，this paper studies TM NDVI on the seasonal normalization methods.By making a demonstration based on Fujian province，we analyze the relationship of quadratic least squares，the relationship of linear least squares and differential relationship between MODIS NDVI in different season phases.Subsequently，we analysis the converted results.The results show that the relationship of quadratic least squares method is superior to the other two methods，which can better eliminate the inconsistency of NDVI seasonal aspects and obtain satisfactory results.By using MODIS NDVI products，we build a quadratic term least squares relationship and fulfill the seasonal normalized of TM NDVI in Fujian province and eliminate the inconsistency of NDVI seasonal aspects at the same time，making the NDVI data in large areas has comparability and consistency on the space.

MODIS NDVI；normalization method；TM NDVI；seasonal difference；least squares

10.3969／j.issn.1000－3177.2015.06.016

TP751

A

1000－3177（2015）142－0081－08

2014－11－13

2015－01－08

国家科技支撑课题（2013BAC08B01）；福建省自然科学基金项目（2012J01166、2013J01158）。

陈贤干（1989—），男，硕士，主要从事遥感应用研究。

E－mail：549352703＠qq.com

汪小钦（1972—），女，博士，研究员，主要从事资源环境遥感应用研究。

E－mail：wangxq＠fzu.edu.cn