太湖MODIS—EVI完全时空数据获取与探索性分析

张恒敢　顾克军　张斯梅

摘要构建太湖水域较长时間范围的植被指数全时空数据有助于全面了解太湖蓝藻水华的时空变化特征。以MODIS数据产品MOD13Q1为数据源，构建了太湖水域MODIS-EVI全时空数据，并做了探索性空间数据分析。结果表明，该研究的时空数据获取和分析方法是可行的，对该数据进行经验正交函数分解，能够获得太湖水面EVI典型的空间分布模式。

关键词水华； 增强植被指数； 探索性空间数据分析

中图分类号S127文献标识码

A文章编号0517-6611（2017）19-226-05

Data Harvesting and Exploratory Analysis on MODISEVI Full Spatiotemperal Data Set for Taihu Lake

ZHANG Henggan，GU Kejun，ZHANG Simei（Institute of Agricultural Resources and Environment， Jiangsu Academy of Agricultural Sciences， Nanjing， Jiangsu 210014）

AbstractLong time spatiotemperal data set was conductive to a comprehensive understanding of the spatiotemperal characteristics of cyanobacterria bloom in Taihu Lake. Taking the MODIS data product MOD13Q1 as data source， an long time， full grid spatiotemporal MODISEVI data set were constructed and analyzed. The results showed that the method described in this paper works well and can be used to obtain the typical spatial patterns of EVI distribution in Taihu Lake.

Key wordsBloom；EVI；ESDA

水华是由于浮游植物的生物量显著高于一般水体中的平均值，导致在水体表面大量聚集而成的藻类聚积体[1]。太湖每年都会有或多或少的水华发生，大面积的水华暴发会造成水体的严重污染，对水生动植物造成危害，甚至危及人体健康。

EOS-MODIS遥感影像数据因时间分辨率高、免费且容易获取常被用于水华监测。黄君等[2]用2010—2013年MODIS影像资料，对太湖水华的时空分布规律进行了研究，结果表明，每年水华的发生高峰在8—10月，西部沿岸是水华首次暴发最频繁的区域。姜晟等[3]利用MODIS数据研究了太湖水华发生与水温的关系。金焰等[4]应用MODIS影像研究了太湖蓝藻水华的时空分布规律。翁建中等[5]选择2007和2008 年200幅MODIS太湖蓝藻监测遥感影像， 统计分析了梅梁湾、竺山湾宜兴段、贡湖湾、东太湖胥口湾和湖州方向湖体蓝藻水华暴发的空间和时间分布规律。

以上研究的共同特征是选取特定的时间点或时间片段的遥感影像作为分析的基础数据，但都未对MODIS获取的全部时间段数据作连续的时空变化分析。要做到这一点，首先需要构建一个基于EOS-MODIS数据的完全时空数据集（Full Spatio-Temporal Data Set）。MOD13Q1是MODIS陆地数据产品，空间分辨率为250 m，时间分辨率16 d，每年23个观测值。太湖水华的聚集在时间上变化很快，采用16 d聚集数据相当于部分平滑，使得水华变化的趋势更加稳定。MOD13Q1数据已经过云处理，省去了数据处理的麻烦。美国航空航天局（NASA）的橡树岭数据中心提供该数据产品的subset web服务，大大提高了服务器的工作效率和客户端下载数据的针对性，这使得获取太湖水面完全时空数据成为可能。因此，笔者采用MOD13Q1数据产品的单像素增强植被指数（EVI）值构建完全时空数据集，并进行探索性空间数据分析（ESDA）[6-7]，以期发现时空数据中可能存在的规律，作为进一步分析和研究的基础。

1数据采集和预处理

1.1采集方法和流程

数据来源于美国国家橡树岭数据中心提供的web服务[8]，作为初步试验，笔者只下载MOD13Q1的EVI数据。下载过程分为3步：

第1步：空间采样。对包含太湖水面的外接矩形区域进行网格点采样，其中落入太湖水面区域的点作为观测样点。该研究采用30×30的网格，900个样点中落入太湖水面的共计440个数据点，构成观测点集。

第2步：下载数据。

以这440个样点作为观测点，下载每个样点的单像元EVI数据。由于数据量较大，从2000年2月28日—2016年4月22日，共计373个观测时间，440个位置，共计164 120个数据，DAAC数据中心为了防止服务器负担过重，限制每次请求最多只能下载10条数据，但实际下载过程中发现，每次下载10条还是经常出错，甚至无法下载。为解决这个问题，采用每次请求一个数据多线程的下载方法，实践证明可以正常下载，但为了减轻服务器负担，同时为了能够稳定下载，下载线程数控制在10个以内。

第3步：查漏补缺。

下载的数据为科学数据格式，经过解析处理后存入SQL Server数据库。由于网络连接不畅、数据提供端繁忙等原因，有些请求无法获得结果。对下载失败的请求必须暂时保存起来，过一段时间重新下载。经过多次查漏补缺，完成MODIS-EVI全时空数据采集过程。

1.2缺失值处理

由于天气等原因，某些位置或时间的数据缺失。如2002年9月30日18个观测点缺失数据。2013年7月28日和2015年3月6日全部观测点数据缺失。应用R语言进行统计分析时，多数函数可以处理缺值，而有些函数的输入数据结构需要完整的数据才能构建。比如，进行时空统计分析时也有些分析方法如spacetime包的经验正交函数分解等需要完整的数据，否则无法进行。因此需要首先对这些缺值进行插值填充。

应用R软件的Zoo包提供的时间序列样条插值进行缺值估计。由于EVI基本不会受到前后1年数据的影响，为减少运算量，提高插值准确性，以缺值所在年份1年的数据作为插值环境数据进行插值估计，对每个观测点根据时间序列进行样条插值得到缺失值。

1.3使用的软件工具

采样点格数据生成使用Croplab软件；数据下载与预处理使用ModisDataCollector程序；描述性数据分析和可视化件使用R和ggplot2包；空间数据分析使用spacetime，gstat包；时间序列分析使用xts包；经验正交函数分解使用sinkr包；扫描统计量分析使用SaTScan独立程序进行，rSatScan包用于将SaTScan输入输出结果导出和导入R环境。

2探索性分析

先进行描述性分析、时空数据可视化、时间序列分析，了解数据的基本特征，发现可能存在的规律，初步判断太湖EVI具有一定的空间分布型和时空聚集特征，然后通过EOF分析提取出时空分布型，通过扫描统计量方法提取其时空聚集特征。

2.1描述性统计该试验采集了太湖水面共440个点单像素EVI数据，时间跨度从2000年2月28日—2016年4月22日共371个时间点。在该尺度下，2013年7月28日和2015年3月6日全部观测点数据缺失。其余时间点缺失数据18条，共缺失记录898条，有效记录数163 222条。描述性数据分析以差值之前的原始数据分析，结果见表1。

安徽农业科学2017年

由于MOD13Q1数据产品已经对缺失值进行了估计和处理，从网站下载的数据缺值较少，在该研究的时间段和分辨率条件下，只有2个观测时间的EVI值全部缺失，占总数据总量的2/373。因此，进行探索性数据分析时，可以忽略缺值的影响。

2.2数据可视化

时空数据可视化是时空探索性分析的基本方法。把各时间点的监测数据绘成地图，显示2000—2016年371帧EVI空间分布图。从这些空间分布图可以大致看出，EVI高值分布的区域及年度、季节间的变化。从图1可以看出，1—3月EVI处于低值，从3月末开始，太湖西部沿岸EVI值开始增加，7—8月达到高值，9月以后开始下降。从单年度来看，各月份之间会有不同，从多年来看，这一趋势是普遍存在的。2000—2016年的EVI空间分布组成一个较长时间的序列，初步观察其空间分布和时间变化，可以获得一些感性的认识：①太湖MODIS-EVI变异有明显的季节变化，冬春季EVI值小，夏秋季EVI值大；②从整个观测时间看，2006年开始太湖EVI有较为明显的上升，2010年以后又有较为明显的降低；③从空间位置看，EVI高值区域总体分布在太湖西海岸、竺山湾、梅梁湾，有些年份会向湖心扩散。

2.3时间序列

MOD13Q1数据产品每隔16 d 1条记录，每年有23条记录。图2为全太湖440个样点的EVI平均值随着时间的变化曲线，每条曲线代表1年的变化。从图2可以看出，EVI平均值存在明显的季节变化，总体趋势是从春季开始EVI值不断上升，到夏秋季达到高峰值，入冬以后迅速下降。

年度之间变化较大，每年内的不同的观测时间也有很大差异。1年内各观测时间点并未表现出稳定的渐进变化趋势，而呈现出较强的随机性。这可能是由云层覆盖的不规则变化等造成。

对某一观测时间全太湖所有观测点的EVI进行平均，再用趋势分解（图3）得到EVI变化的趋势项、季节项和误差项，结果见图3。从图3可以看出，观测值显得杂乱无章，但经过分解后，趋势项能够看出明显的单峰，从2006年开始明显升高，到2007年达到高峰，这和太湖2007年大规模水华暴发是一致的。EVI高值持续到2010年开始下降，从2011年至今处于较低状态，这与实际观测也是一致的。蓝藻大暴发后，我国各级政府十分重视，采用了相关治理措施，限制了太湖流域的污染排放，取得好的效果，自2010年至今，太湖未出现大规模水华暴发。

2.4经验正交函数分解

经验正交函数（Empirical Orthogonal Function， EOF）分解是地学统计中常用的分析方法。其实质是基于对时空数据的方差-协方差矩阵或相关矩阵的主成分分析（Principal Component Analysis，PCA），把空間上的每个点作为1个变量来看，根据时间维的变异，经过正交变换，建立若干新的变量线性组合，这些组合变量之间彼此正交。每个组合的空间分布即为1个空间分布型（一些文献称作空间模态）。特定空间分布型在时间维变异的贡献又叫做时间系数。因此，经验正交函数分解也称为时空分解。这些线性无关的空间分布型可以理解为由彼此互不相关的因素独立作用的结果。

通过经验正交函数分解方法能够降低数据繁多造成的干扰，更清楚地看出EVI变化的组成和时间分配。太湖MODIS-EVI数据EOF分解的结果见表2。

前8个主成分方差之和占总方差的41.62%（均显著，置信概率0.05）。除第1主成分解释方差为21.61%，其余各主成分解释的方差相对平均，说明太湖水域EVI特征的影响因素较多，作用模式复杂，具有决定作用的因子较少，因此随机性较大。图4显示了EOF分解获得的前8个空间分布型。结合表2可知，各空间分布型的差异较大。EOF1占总方差的21.61%，描述了太湖EVI变异的总体趋势，其物理意义尚不明确。EOF2占总方差的5.34%，其聚集区位于太湖的湖心区，其空间分布模式与太湖的水深分布相近，推测其变异来源于水位深度。EOF3占总方差的5.00%，其聚集区域位于竺山湾和梅湖湾附近，推测其变异来源于北部的太湖区水系。EOF4占总方差的2.79%，其聚集区靠近太湖西岸的河流入湖口，推测其变异来源于湖西区水系。这种相关性在提高采样分辨率（60×70网格）时更加明显。当然，仅凭分布的图像模式和空间分布的距离难以确定这些变异来源于这些物理因素，还需要进一步的空间相关分析和模式匹配才能进一步确定。

这些空间分布型揭示了隐藏在总方差背后各个自主性的变异，虽然不一定代表某一个具体的物理因素效应，但它们是1个或几个协同的独立变化因素的外在表现。作为探索性分析，这些分布模式为今后研究太湖水华的影响因素提供了一个新的视角。

2.5时空聚集

对整个太湖水面的MODIS-EVI平均值进行时间序列分析能够从整体上对EVI的变化做一个全局的了解，但是太湖EVI的变化是分区的，一些部分变化较大，如竺山湾、西部沿岸，而湖心区和东太湖变化较小，峰值也较底，表现出明显的区域性、聚集性特征。因此，需要对湖面进行分区，从空间上划分成若干个聚集区作为独立的单元进行分析。因此，首先需要对太湖MODIS-EVI数据进行聚集分析。

采用kuldoff时空扫描统计量方法进行时空聚集分析，结果显示有2个时空聚集区。第1个是西部聚集区，从2006年5月30日至2010年10月31日。时间上和时间序列分析的结果一致，空间上则确定了时空聚集的地理范围在太湖西部沿岸、竺山湾、梅梁湾组成的区域。第2个是东部聚集区，位于东太湖，是水生植物生长的地方（图5）。

以上只是初步的扫描结果大致能够看出EVI的时空分布，这与实际情况相符。据监测数据，2010—2013 年太湖总氮浓度的逐年下降抑制了水华蓝藻的生长、繁殖，这与2010—2013年4—10月太湖平均蓝藻密度的下降也是相吻合的[4]。这也说明用时空扫描统计量方法可以用于太湖流域EVI的时空分析。但以上结果只是初步分析结果，尚未对扫描参数进行筛选，数据点密度仍较低，有待于进一步研究。

3结论与讨论

该研究证实了该数据获取方法的有效性，并对所得数据进行初步分析，得出以下结论：

（1）太湖MODIS-EVI完全时空数据获取和分析的方法是可行的，所得太湖EVI時空数据蕴含大量有价值的信息，通过时空数据分析能获得关于太湖水面EVI时空变化宏观的认识，而这些信息是片段的时空数据所不能提供的。

（2）要从这些数据中获得可靠、准确、有价值的信息，需要有效、恰当的数据分析挖掘方法。经验正交函数方法对于太湖EVI数据有很好的时空分解效果，能够确定水华的

典型时空分布型，作为太湖的基本度量特征，如果与湖区气象、水文、植被甚至流

域种植结构等因素作空间相关分析，则可能更加精细地揭示太湖蓝藻水华的形成机制。时间序列分析可使研究者对太湖EVI的多年来变化趋势有一个量化的了解。时空聚集分析也是有效的，能够找到水华在时间和空间上的聚集特征，但具体的参数设置则需要进一步探索优化。

（3）该研究的数据基于30×30网格的空间分辨率偏小，在进行进一步的主成分、因子分析、空间自相关分析、扫描统计量分析时需要更加精细的时空数据，增加空间分辨率。另外，低空间分辨率数据提取的水华会产生尺度误差[9]，因此还需要考虑尺度变化的影响。

参考文献

[1]

徐京萍，张柏，李方，等.基于MODIS数据的太湖藻华水体识别模式[J].湖泊科学，2008，20（2）：191-195.

[2] 黄君，宋挺，庄严，等.基于MODIS 数据的太湖蓝藻水华时空分布规律研究[J].四川环境，2014，33（5）：27-33.

[3] 姜晟，张咏，蒋建军，等.基于MODIS 数据的太湖蓝藻变化与水温关系研究[J].环境科技，2009，22（6）：28-31.

[4] 金焰，张咏，姜晟.EOS/MODIS数据在太湖蓝藻水华时空分布规律提取中的应用研究[J].环境科技，2009，22（S2）：9-11，14.

[5] 翁建中，李继影，梁柱，等.太湖蓝藻水华时空分布与预警监测响应的分析[J].环境监控与预警，2010，2（3）：1-4.

[6] SOLANO R，DIDAN K，JACOBSON A，et al.MODIS vegetation indices （MOD13） C5 users guide[M].Tucson，AZ，USA：University of Arizona，2010.

[7] 罗伯特·海宁.空间数据分析理论与实践[M].李建松，秦昆，译.武汉：武汉大学出版社，2009.

[8] ORNL DAAC.MODIS collection 5 land products global subsetting and visualization tool[Z].Oak Ridge， Tennessee， USA：

ORNL DAAC，2008.

[9] 尚琳琳，马荣华，段洪涛，等.利用MODIS 影像提取太湖蓝藻水华的尺度差异性分析[J].湖泊科学，2011，23（6）：847-854.