太湖MODIS—EVI完全時空數(shù)據(jù)獲取與探索性分析

張恒敢　顧克軍　張斯梅

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

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

中圖分類號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]。太湖每年都會有或多或少的水華發(fā)生，大面積的水華暴發(fā)會造成水體的嚴重污染，對水生動植物造成危害，甚至危及人體健康。

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

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

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

1.1采集方法和流程

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

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

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

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

第3步：查漏補缺。

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

1.2缺失值處理

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

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

1.3使用的軟件工具

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

2探索性分析

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

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

安徽農業(yè)科學2017年

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

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

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

2.3時間序列

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

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

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

2.4經驗正交函數(shù)分解

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

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

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

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

2.5時空聚集

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

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

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

3結論與討論

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

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

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

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

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

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

參考文獻

[1]

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

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

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

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

[5] 翁建中，李繼影，梁柱，等.太湖藍藻水華時空分布與預警監(jiān)測響應的分析[J].環(huán)境監(jiān)控與預警，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ù)分析理論與實踐[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.