張恒敢　顧克軍　張斯梅

摘要構(gòu)建太湖水域較長(zhǎng)時(shí)間范圍的植被指數(shù)全時(shí)空數(shù)據(jù)有助于全面了解太湖藍(lán)藻水華的時(shí)空變化特征。以MODIS數(shù)據(jù)產(chǎn)品MOD13Q1為數(shù)據(jù)源，構(gòu)建了太湖水域MODIS-EVI全時(shí)空數(shù)據(jù)，并做了探索性空間數(shù)據(jù)分析。結(jié)果表明，該研究的時(shí)空數(shù)據(jù)獲取和分析方法是可行的，對(duì)該數(shù)據(jù)進(jìn)行經(jīng)驗(yàn)正交函數(shù)分解，能夠獲得太湖水面EVI典型的空間分布模式。

關(guān)鍵詞水華； 增強(qiáng)植被指數(shù)； 探索性空間數(shù)據(jù)分析

中圖分類(lèi)號(hào)S127文獻(xiàn)標(biāo)識(shí)碼

A文章編號(hào)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

水華是由于浮游植物的生物量顯著高于一般水體中的平均值，導(dǎo)致在水體表面大量聚集而成的藻類(lèi)聚積體[1]。太湖每年都會(huì)有或多或少的水華發(fā)生，大面積的水華暴發(fā)會(huì)造成水體的嚴(yán)重污染，對(duì)水生動(dòng)植物造成危害，甚至危及人體健康。

EOS-MODIS遙感影像數(shù)據(jù)因時(shí)間分辨率高、免費(fèi)且容易獲取常被用于水華監(jiān)測(cè)。黃君等[2]用2010—2013年MODIS影像資料，對(duì)太湖水華的時(shí)空分布規(guī)律進(jìn)行了研究，結(jié)果表明，每年水華的發(fā)生高峰在8—10月，西部沿岸是水華首次暴發(fā)最頻繁的區(qū)域。姜晟等[3]利用MODIS數(shù)據(jù)研究了太湖水華發(fā)生與水溫的關(guān)系。金焰等[4]應(yīng)用MODIS影像研究了太湖藍(lán)藻水華的時(shí)空分布規(guī)律。翁建中等[5]選擇2007和2008 年200幅MODIS太湖藍(lán)藻監(jiān)測(cè)遙感影像， 統(tǒng)計(jì)分析了梅梁灣、竺山灣宜興段、貢湖灣、東太湖胥口灣和湖州方向湖體藍(lán)藻水華暴發(fā)的空間和時(shí)間分布規(guī)律。

以上研究的共同特征是選取特定的時(shí)間點(diǎn)或時(shí)間片段的遙感影像作為分析的基礎(chǔ)數(shù)據(jù)，但都未對(duì)MODIS獲取的全部時(shí)間段數(shù)據(jù)作連續(xù)的時(shí)空變化分析。要做到這一點(diǎn)，首先需要構(gòu)建一個(gè)基于EOS-MODIS數(shù)據(jù)的完全時(shí)空數(shù)據(jù)集（Full Spatio-Temporal Data Set）。MOD13Q1是MODIS陸地?cái)?shù)據(jù)產(chǎn)品，空間分辨率為250 m，時(shí)間分辨率16 d，每年23個(gè)觀測(cè)值。太湖水華的聚集在時(shí)間上變化很快，采用16 d聚集數(shù)據(jù)相當(dāng)于部分平滑，使得水華變化的趨勢(shì)更加穩(wěn)定。MOD13Q1數(shù)據(jù)已經(jīng)過(guò)云處理，省去了數(shù)據(jù)處理的麻煩。美國(guó)航空航天局（NASA）的橡樹(shù)嶺數(shù)據(jù)中心提供該數(shù)據(jù)產(chǎn)品的subset web服務(wù)，大大提高了服務(wù)器的工作效率和客戶(hù)端下載數(shù)據(jù)的針對(duì)性，這使得獲取太湖水面完全時(shí)空數(shù)據(jù)成為可能。因此，筆者采用MOD13Q1數(shù)據(jù)產(chǎn)品的單像素增強(qiáng)植被指數(shù)（EVI）值構(gòu)建完全時(shí)空數(shù)據(jù)集，并進(jìn)行探索性空間數(shù)據(jù)分析（ESDA）[6-7]，以期發(fā)現(xiàn)時(shí)空數(shù)據(jù)中可能存在的規(guī)律，作為進(jìn)一步分析和研究的基礎(chǔ)。

1數(shù)據(jù)采集和預(yù)處理

1.1采集方法和流程

數(shù)據(jù)來(lái)源于美國(guó)國(guó)家橡樹(shù)嶺數(shù)據(jù)中心提供的web服務(wù)[8]，作為初步試驗(yàn)，筆者只下載MOD13Q1的EVI數(shù)據(jù)。下載過(guò)程分為3步：

第1步：空間采樣。對(duì)包含太湖水面的外接矩形區(qū)域進(jìn)行網(wǎng)格點(diǎn)采樣，其中落入太湖水面區(qū)域的點(diǎn)作為觀測(cè)樣點(diǎn)。該研究采用30×30的網(wǎng)格，900個(gè)樣點(diǎn)中落入太湖水面的共計(jì)440個(gè)數(shù)據(jù)點(diǎn)，構(gòu)成觀測(cè)點(diǎn)集。

第2步：下載數(shù)據(jù)。

以這440個(gè)樣點(diǎn)作為觀測(cè)點(diǎn)，下載每個(gè)樣點(diǎn)的單像元EVI數(shù)據(jù)。由于數(shù)據(jù)量較大，從2000年2月28日—2016年4月22日，共計(jì)373個(gè)觀測(cè)時(shí)間，440個(gè)位置，共計(jì)164 120個(gè)數(shù)據(jù)，DAAC數(shù)據(jù)中心為了防止服務(wù)器負(fù)擔(dān)過(guò)重，限制每次請(qǐng)求最多只能下載10條數(shù)據(jù)，但實(shí)際下載過(guò)程中發(fā)現(xiàn)，每次下載10條還是經(jīng)常出錯(cuò)，甚至無(wú)法下載。為解決這個(gè)問(wèn)題，采用每次請(qǐng)求一個(gè)數(shù)據(jù)多線程的下載方法，實(shí)踐證明可以正常下載，但為了減輕服務(wù)器負(fù)擔(dān)，同時(shí)為了能夠穩(wěn)定下載，下載線程數(shù)控制在10個(gè)以?xún)?nèi)。

第3步：查漏補(bǔ)缺。

下載的數(shù)據(jù)為科學(xué)數(shù)據(jù)格式，經(jīng)過(guò)解析處理后存入SQL Server數(shù)據(jù)庫(kù)。由于網(wǎng)絡(luò)連接不暢、數(shù)據(jù)提供端繁忙等原因，有些請(qǐng)求無(wú)法獲得結(jié)果。對(duì)下載失敗的請(qǐng)求必須暫時(shí)保存起來(lái)，過(guò)一段時(shí)間重新下載。經(jīng)過(guò)多次查漏補(bǔ)缺，完成MODIS-EVI全時(shí)空數(shù)據(jù)采集過(guò)程。

1.2缺失值處理

由于天氣等原因，某些位置或時(shí)間的數(shù)據(jù)缺失。如2002年9月30日18個(gè)觀測(cè)點(diǎn)缺失數(shù)據(jù)。2013年7月28日和2015年3月6日全部觀測(cè)點(diǎn)數(shù)據(jù)缺失。應(yīng)用R語(yǔ)言進(jìn)行統(tǒng)計(jì)分析時(shí)，多數(shù)函數(shù)可以處理缺值，而有些函數(shù)的輸入數(shù)據(jù)結(jié)構(gòu)需要完整的數(shù)據(jù)才能構(gòu)建。比如，進(jìn)行時(shí)空統(tǒng)計(jì)分析時(shí)也有些分析方法如spacetime包的經(jīng)驗(yàn)正交函數(shù)分解等需要完整的數(shù)據(jù)，否則無(wú)法進(jìn)行。因此需要首先對(duì)這些缺值進(jìn)行插值填充。

應(yīng)用R軟件的Zoo包提供的時(shí)間序列樣條插值進(jìn)行缺值估計(jì)。由于EVI基本不會(huì)受到前后1年數(shù)據(jù)的影響，為減少運(yùn)算量，提高插值準(zhǔn)確性，以缺值所在年份1年的數(shù)據(jù)作為插值環(huán)境數(shù)據(jù)進(jìn)行插值估計(jì)，對(duì)每個(gè)觀測(cè)點(diǎn)根據(jù)時(shí)間序列進(jìn)行樣條插值得到缺失值。

1.3使用的軟件工具

采樣點(diǎn)格數(shù)據(jù)生成使用Croplab軟件；數(shù)據(jù)下載與預(yù)處理使用ModisDataCollector程序；描述性數(shù)據(jù)分析和可視化件使用R和ggplot2包；空間數(shù)據(jù)分析使用spacetime，gstat包；時(shí)間序列分析使用xts包；經(jīng)驗(yàn)正交函數(shù)分解使用sinkr包；掃描統(tǒng)計(jì)量分析使用SaTScan獨(dú)立程序進(jìn)行，rSatScan包用于將SaTScan輸入輸出結(jié)果導(dǎo)出和導(dǎo)入R環(huán)境。

2探索性分析

先進(jìn)行描述性分析、時(shí)空數(shù)據(jù)可視化、時(shí)間序列分析，了解數(shù)據(jù)的基本特征，發(fā)現(xiàn)可能存在的規(guī)律，初步判斷太湖EVI具有一定的空間分布型和時(shí)空聚集特征，然后通過(guò)EOF分析提取出時(shí)空分布型，通過(guò)掃描統(tǒng)計(jì)量方法提取其時(shí)空聚集特征。

2.1描述性統(tǒng)計(jì)該試驗(yàn)采集了太湖水面共440個(gè)點(diǎn)單像素EVI數(shù)據(jù)，時(shí)間跨度從2000年2月28日—2016年4月22日共371個(gè)時(shí)間點(diǎn)。在該尺度下，2013年7月28日和2015年3月6日全部觀測(cè)點(diǎn)數(shù)據(jù)缺失。其余時(shí)間點(diǎn)缺失數(shù)據(jù)18條，共缺失記錄898條，有效記錄數(shù)163 222條。描述性數(shù)據(jù)分析以差值之前的原始數(shù)據(jù)分析，結(jié)果見(jiàn)表1。

安徽農(nóng)業(yè)科學(xué)2017年

由于MOD13Q1數(shù)據(jù)產(chǎn)品已經(jīng)對(duì)缺失值進(jìn)行了估計(jì)和處理，從網(wǎng)站下載的數(shù)據(jù)缺值較少，在該研究的時(shí)間段和分辨率條件下，只有2個(gè)觀測(cè)時(shí)間的EVI值全部缺失，占總數(shù)據(jù)總量的2/373。因此，進(jìn)行探索性數(shù)據(jù)分析時(shí)，可以忽略缺值的影響。

2.2數(shù)據(jù)可視化

時(shí)空數(shù)據(jù)可視化是時(shí)空探索性分析的基本方法。把各時(shí)間點(diǎn)的監(jiān)測(cè)數(shù)據(jù)繪成地圖，顯示2000—2016年371幀EVI空間分布圖。從這些空間分布圖可以大致看出，EVI高值分布的區(qū)域及年度、季節(jié)間的變化。從圖1可以看出，1—3月EVI處于低值，從3月末開(kāi)始，太湖西部沿岸EVI值開(kāi)始增加，7—8月達(dá)到高值，9月以后開(kāi)始下降。從單年度來(lái)看，各月份之間會(huì)有不同，從多年來(lái)看，這一趨勢(shì)是普遍存在的。2000—2016年的EVI空間分布組成一個(gè)較長(zhǎng)時(shí)間的序列，初步觀察其空間分布和時(shí)間變化，可以獲得一些感性的認(rèn)識(shí)：①太湖MODIS-EVI變異有明顯的季節(jié)變化，冬春季EVI值小，夏秋季EVI值大；②從整個(gè)觀測(cè)時(shí)間看，2006年開(kāi)始太湖EVI有較為明顯的上升，2010年以后又有較為明顯的降低；③從空間位置看，EVI高值區(qū)域總體分布在太湖西海岸、竺山灣、梅梁灣，有些年份會(huì)向湖心擴(kuò)散。

2.3時(shí)間序列

MOD13Q1數(shù)據(jù)產(chǎn)品每隔16 d 1條記錄，每年有23條記錄。圖2為全太湖440個(gè)樣點(diǎn)的EVI平均值隨著時(shí)間的變化曲線，每條曲線代表1年的變化。從圖2可以看出，EVI平均值存在明顯的季節(jié)變化，總體趨勢(shì)是從春季開(kāi)始EVI值不斷上升，到夏秋季達(dá)到高峰值，入冬以后迅速下降。

年度之間變化較大，每年內(nèi)的不同的觀測(cè)時(shí)間也有很大差異。1年內(nèi)各觀測(cè)時(shí)間點(diǎn)并未表現(xiàn)出穩(wěn)定的漸進(jìn)變化趨勢(shì)，而呈現(xiàn)出較強(qiáng)的隨機(jī)性。這可能是由云層覆蓋的不規(guī)則變化等造成。

對(duì)某一觀測(cè)時(shí)間全太湖所有觀測(cè)點(diǎn)的EVI進(jìn)行平均，再用趨勢(shì)分解（圖3）得到EVI變化的趨勢(shì)項(xiàng)、季節(jié)項(xiàng)和誤差項(xiàng)，結(jié)果見(jiàn)圖3。從圖3可以看出，觀測(cè)值顯得雜亂無(wú)章，但經(jīng)過(guò)分解后，趨勢(shì)項(xiàng)能夠看出明顯的單峰，從2006年開(kāi)始明顯升高，到2007年達(dá)到高峰，這和太湖2007年大規(guī)模水華暴發(fā)是一致的。EVI高值持續(xù)到2010年開(kāi)始下降，從2011年至今處于較低狀態(tài)，這與實(shí)際觀測(cè)也是一致的。藍(lán)藻大暴發(fā)后，我國(guó)各級(jí)政府十分重視，采用了相關(guān)治理措施，限制了太湖流域的污染排放，取得好的效果，自2010年至今，太湖未出現(xiàn)大規(guī)模水華暴發(fā)。

2.4經(jīng)驗(yàn)正交函數(shù)分解

經(jīng)驗(yàn)正交函數(shù)（Empirical Orthogonal Function， EOF）分解是地學(xué)統(tǒng)計(jì)中常用的分析方法。其實(shí)質(zhì)是基于對(duì)時(shí)空數(shù)據(jù)的方差-協(xié)方差矩陣或相關(guān)矩陣的主成分分析（Principal Component Analysis，PCA），把空間上的每個(gè)點(diǎn)作為1個(gè)變量來(lái)看，根據(jù)時(shí)間維的變異，經(jīng)過(guò)正交變換，建立若干新的變量線性組合，這些組合變量之間彼此正交。每個(gè)組合的空間分布即為1個(gè)空間分布型（一些文獻(xiàn)稱(chēng)作空間模態(tài)）。特定空間分布型在時(shí)間維變異的貢獻(xiàn)又叫做時(shí)間系數(shù)。因此，經(jīng)驗(yàn)正交函數(shù)分解也稱(chēng)為時(shí)空分解。這些線性無(wú)關(guān)的空間分布型可以理解為由彼此互不相關(guān)的因素獨(dú)立作用的結(jié)果。

通過(guò)經(jīng)驗(yàn)正交函數(shù)分解方法能夠降低數(shù)據(jù)繁多造成的干擾，更清楚地看出EVI變化的組成和時(shí)間分配。太湖MODIS-EVI數(shù)據(jù)EOF分解的結(jié)果見(jiàn)表2。

前8個(gè)主成分方差之和占總方差的41.62%（均顯著，置信概率0.05）。除第1主成分解釋方差為21.61%，其余各主成分解釋的方差相對(duì)平均，說(shuō)明太湖水域EVI特征的影響因素較多，作用模式復(fù)雜，具有決定作用的因子較少，因此隨機(jī)性較大。圖4顯示了EOF分解獲得的前8個(gè)空間分布型。結(jié)合表2可知，各空間分布型的差異較大。EOF1占總方差的21.61%，描述了太湖EVI變異的總體趨勢(shì)，其物理意義尚不明確。EOF2占總方差的5.34%，其聚集區(qū)位于太湖的湖心區(qū)，其空間分布模式與太湖的水深分布相近，推測(cè)其變異來(lái)源于水位深度。EOF3占總方差的5.00%，其聚集區(qū)域位于竺山灣和梅湖灣附近，推測(cè)其變異來(lái)源于北部的太湖區(qū)水系。EOF4占總方差的2.79%，其聚集區(qū)靠近太湖西岸的河流入湖口，推測(cè)其變異來(lái)源于湖西區(qū)水系。這種相關(guān)性在提高采樣分辨率（60×70網(wǎng)格）時(shí)更加明顯。當(dāng)然，僅憑分布的圖像模式和空間分布的距離難以確定這些變異來(lái)源于這些物理因素，還需要進(jìn)一步的空間相關(guān)分析和模式匹配才能進(jìn)一步確定。

這些空間分布型揭示了隱藏在總方差背后各個(gè)自主性的變異，雖然不一定代表某一個(gè)具體的物理因素效應(yīng)，但它們是1個(gè)或幾個(gè)協(xié)同的獨(dú)立變化因素的外在表現(xiàn)。作為探索性分析，這些分布模式為今后研究太湖水華的影響因素提供了一個(gè)新的視角。

2.5時(shí)空聚集

對(duì)整個(gè)太湖水面的MODIS-EVI平均值進(jìn)行時(shí)間序列分析能夠從整體上對(duì)EVI的變化做一個(gè)全局的了解，但是太湖EVI的變化是分區(qū)的，一些部分變化較大，如竺山灣、西部沿岸，而湖心區(qū)和東太湖變化較小，峰值也較底，表現(xiàn)出明顯的區(qū)域性、聚集性特征。因此，需要對(duì)湖面進(jìn)行分區(qū)，從空間上劃分成若干個(gè)聚集區(qū)作為獨(dú)立的單元進(jìn)行分析。因此，首先需要對(duì)太湖MODIS-EVI數(shù)據(jù)進(jìn)行聚集分析。

采用kuldoff時(shí)空掃描統(tǒng)計(jì)量方法進(jìn)行時(shí)空聚集分析，結(jié)果顯示有2個(gè)時(shí)空聚集區(qū)。第1個(gè)是西部聚集區(qū)，從2006年5月30日至2010年10月31日。時(shí)間上和時(shí)間序列分析的結(jié)果一致，空間上則確定了時(shí)空聚集的地理范圍在太湖西部沿岸、竺山灣、梅梁灣組成的區(qū)域。第2個(gè)是東部聚集區(qū)，位于東太湖，是水生植物生長(zhǎng)的地方（圖5）。

以上只是初步的掃描結(jié)果大致能夠看出EVI的時(shí)空分布，這與實(shí)際情況相符。據(jù)監(jiān)測(cè)數(shù)據(jù)，2010—2013 年太湖總氮濃度的逐年下降抑制了水華藍(lán)藻的生長(zhǎng)、繁殖，這與2010—2013年4—10月太湖平均藍(lán)藻密度的下降也是相吻合的[4]。這也說(shuō)明用時(shí)空掃描統(tǒng)計(jì)量方法可以用于太湖流域EVI的時(shí)空分析。但以上結(jié)果只是初步分析結(jié)果，尚未對(duì)掃描參數(shù)進(jìn)行篩選，數(shù)據(jù)點(diǎn)密度仍較低，有待于進(jìn)一步研究。

3結(jié)論與討論

該研究證實(shí)了該數(shù)據(jù)獲取方法的有效性，并對(duì)所得數(shù)據(jù)進(jìn)行初步分析，得出以下結(jié)論：

（1）太湖MODIS-EVI完全時(shí)空數(shù)據(jù)獲取和分析的方法是可行的，所得太湖EVI時(shí)空數(shù)據(jù)蘊(yùn)含大量有價(jià)值的信息，通過(guò)時(shí)空數(shù)據(jù)分析能獲得關(guān)于太湖水面EVI時(shí)空變化宏觀的認(rèn)識(shí)，而這些信息是片段的時(shí)空數(shù)據(jù)所不能提供的。

（2）要從這些數(shù)據(jù)中獲得可靠、準(zhǔn)確、有價(jià)值的信息，需要有效、恰當(dāng)?shù)臄?shù)據(jù)分析挖掘方法。經(jīng)驗(yàn)正交函數(shù)方法對(duì)于太湖EVI數(shù)據(jù)有很好的時(shí)空分解效果，能夠確定水華的

典型時(shí)空分布型，作為太湖的基本度量特征，如果與湖區(qū)氣象、水文、植被甚至流

域種植結(jié)構(gòu)等因素作空間相關(guān)分析，則可能更加精細(xì)地揭示太湖藍(lán)藻水華的形成機(jī)制。時(shí)間序列分析可使研究者對(duì)太湖EVI的多年來(lái)變化趨勢(shì)有一個(gè)量化的了解。時(shí)空聚集分析也是有效的，能夠找到水華在時(shí)間和空間上的聚集特征，但具體的參數(shù)設(shè)置則需要進(jìn)一步探索優(yōu)化。

（3）該研究的數(shù)據(jù)基于30×30網(wǎng)格的空間分辨率偏小，在進(jìn)行進(jìn)一步的主成分、因子分析、空間自相關(guān)分析、掃描統(tǒng)計(jì)量分析時(shí)需要更加精細(xì)的時(shí)空數(shù)據(jù)，增加空間分辨率。另外，低空間分辨率數(shù)據(jù)提取的水華會(huì)產(chǎn)生尺度誤差[9]，因此還需要考慮尺度變化的影響。

參考文獻(xiàn)

[1]

徐京萍，張柏，李方，等.基于MODIS數(shù)據(jù)的太湖藻華水體識(shí)別模式[J].湖泊科學(xué)，2008，20（2）：191-195.

[2] 黃君，宋挺，莊嚴(yán)，等.基于MODIS 數(shù)據(jù)的太湖藍(lán)藻水華時(shí)空分布規(guī)律研究[J].四川環(huán)境，2014，33（5）：27-33.

[3] 姜晟，張?jiān)?，蔣建軍，等.基于MODIS 數(shù)據(jù)的太湖藍(lán)藻變化與水溫關(guān)系研究[J].環(huán)境科技，2009，22（6）：28-31.

[4] 金焰，張?jiān)?，姜?EOS/MODIS數(shù)據(jù)在太湖藍(lán)藻水華時(shí)空分布規(guī)律提取中的應(yīng)用研究[J].環(huán)境科技，2009，22（S2）：9-11，14.

[5] 翁建中，李繼影，梁柱，等.太湖藍(lán)藻水華時(shí)空分布與預(yù)警監(jiān)測(cè)響應(yīng)的分析[J].環(huán)境監(jiān)控與預(yù)警，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] 羅伯特·海寧.空間數(shù)據(jù)分析理論與實(shí)踐[M].李建松，秦昆，譯.武漢：武漢大學(xué)出版社，2009.

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

ORNL DAAC，2008.

[9] 尚琳琳，馬榮華，段洪濤，等.利用MODIS 影像提取太湖藍(lán)藻水華的尺度差異性分析[J].湖泊科學(xué)，2011，23（6）：847-854.