從NDVI序列檢測(cè)植被“綠化”或“褐化”變化趨勢(shì)的新方法

唐亮, 趙忠明, 唐娉

(1.海南熱帶海洋學(xué)院,三亞 572022; 2.中國(guó)科學(xué)院遙感與數(shù)字地球研究所,北京 100101)

0 引言

衛(wèi)星數(shù)據(jù)是當(dāng)今陸地生態(tài)系統(tǒng)監(jiān)測(cè)的主要信息源。幾十a(chǎn)的衛(wèi)星時(shí)間序列為定量描述全球植被活動(dòng)及其變化趨勢(shì)提供了數(shù)據(jù)基礎(chǔ)。植被作為陸地生物圈的主要組成部分，是氣候變化中的關(guān)鍵因素[1]，而且人們認(rèn)為氣候變暖將強(qiáng)烈影響陸地生物圈[2]。植被狀態(tài)常用于評(píng)估自然和農(nóng)用土地的生產(chǎn)率和退化[3-5]，其變化趨勢(shì)被認(rèn)為與土地退化關(guān)聯(lián)[6-8]。陸地生態(tài)系統(tǒng)的時(shí)間動(dòng)態(tài)包括漸變、突然變化和沒(méi)有變化。沒(méi)有變化是很少見(jiàn)的。在幾十a(chǎn)的時(shí)間框架內(nèi)，一般認(rèn)為突然的植被減少，可能由一些短期的過(guò)程導(dǎo)致，如火災(zāi)、作物收割或發(fā)生了其他災(zāi)害； 而降雨事件或雪蓋的減少也會(huì)引起突然的植被變綠。植被漸變會(huì)發(fā)生在更長(zhǎng)的時(shí)間內(nèi)，體現(xiàn)植被對(duì)全球變化的適應(yīng)過(guò)程，包括大洋振蕩，或持續(xù)的氣候變化，如年際降水減少或大氣中二氧化碳濃度增加。整個(gè)序列的漸變是最不復(fù)雜的變化類型，這種變化也指單調(diào)的“綠化(greening)”或“褐化(browning)”。歸一化植被指數(shù)(normalized difference vegetation index，NDVI)與植被生產(chǎn)率強(qiáng)相關(guān)[9-10]，NDVI的趨勢(shì)可近似表達(dá)植被“綠化”或“褐化”的趨勢(shì)[11-12]。對(duì)于長(zhǎng)時(shí)間序列，植被活動(dòng)及其趨勢(shì)分析方法已經(jīng)有了一些研究[13-14]，其中最常使用的長(zhǎng)時(shí)間序列趨勢(shì)分析方法包括參數(shù)化和非參數(shù)化方法，前者有最小二乘的一元回歸變化斜率法； 后者包括Sen趨勢(shì)度估計(jì)、Mann-Kendall(MK)方法及季節(jié)性MK方法等。Beurs等[13]指出，基于最小二乘回歸得到的NDVI曲線趨勢(shì)的有效性會(huì)受到影響，因?yàn)樯镂锢韰?shù)的時(shí)間序列在時(shí)間上有相關(guān)性，會(huì)導(dǎo)致最小二乘回歸方法應(yīng)用的假設(shè)條件遭到破壞，這些條件是： ①所有的函數(shù)值應(yīng)當(dāng)是彼此獨(dú)立的； ②殘差有隨機(jī)性； ③殘差是零均值的； ④殘差的方差應(yīng)當(dāng)對(duì)所有時(shí)間點(diǎn)是一樣的。

另一方面，Belle等[15]指出Sen趨勢(shì)度估計(jì)、MK方法及季節(jié)性MK方法具有一個(gè)隱含的假設(shè)，即季節(jié)間趨勢(shì)必須是一致的，否則有實(shí)例表明每個(gè)季節(jié)有明顯的趨勢(shì)，但整體統(tǒng)計(jì)卻沒(méi)有趨勢(shì)。因此在趨勢(shì)性分析時(shí)需同時(shí)對(duì)時(shí)間序列的周期性或季節(jié)性建模，以去除整體統(tǒng)計(jì)對(duì)季節(jié)變化趨勢(shì)的掩蓋。故而進(jìn)行趨勢(shì)性分析是當(dāng)前時(shí)間序列趨勢(shì)分析的主要內(nèi)容之一。因地球繞太陽(yáng)公轉(zhuǎn)的周期性影響引起的氣候、地下水等多種周期性影響及地形的影響，NDVI序列具有強(qiáng)烈的周期性或多種周期疊加效應(yīng)，周期性影響是NDVI趨勢(shì)性分析必須考慮的因素。時(shí)頻分析是周期性分析的基本手段，其中Fourier分析和小波分析是最常使用的時(shí)頻分析工具，但對(duì)時(shí)間序列進(jìn)行Fourier分解和小波分解均尚不能直接得到時(shí)間序列的趨勢(shì)分量。Verbesselt等[14]在Fourier分解基礎(chǔ)上提出了季節(jié)—趨勢(shì)分解模型，將具有周期性特征的時(shí)間序列曲線擬合分解為趨勢(shì)項(xiàng)(一階線性函數(shù))、季節(jié)項(xiàng)(周期函數(shù)組合)及殘差的疊加，從而將季節(jié)項(xiàng)和趨勢(shì)項(xiàng)分離出來(lái)，直接得到分段表示的NDVI變化趨勢(shì)，但因?yàn)槭菙M合方法，其趨勢(shì)項(xiàng)的估計(jì)仍然采用了最小二乘估計(jì)。不同于小波分析，同樣適合于非平穩(wěn)信號(hào)的新的時(shí)頻分析工具經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition，EMD)，可將時(shí)間序列中不同頻率和振幅的信息很好地分解出來(lái)，隨著分解的分量中頻率逐漸降低而振幅增大，可分解得到表征原序列整體趨勢(shì)的分量，該分量是近似單調(diào)變化的序列，即通常所說(shuō)的線性趨勢(shì)項(xiàng)，具有非常明確的物理意義。

本文以2006—2015年間中國(guó)局部地區(qū)的NDVI數(shù)據(jù)為數(shù)據(jù)源，探求EMD方法結(jié)合MK顯著性檢驗(yàn)進(jìn)行植被單調(diào)性變化分析的方法。實(shí)驗(yàn)證明該方法是一種有效的時(shí)間序列變化趨勢(shì)分析方法。

1 NDVI變化趨勢(shì)分析方法

1.1 EMD方法

EMD方法是Huang等[16]于1998年提出的一種時(shí)頻信號(hào)分析方法，屬于自適應(yīng)局部時(shí)頻分析方法，它能根據(jù)數(shù)據(jù)信號(hào)本身的特性將其分解為有限個(gè)本征模態(tài)函數(shù)(intrinsic mode function，IMF)的疊加，所分解出來(lái)的各IMF分量包含了原信號(hào)不同時(shí)間尺度的局部特征信號(hào)，特別適合非平穩(wěn)信號(hào)的分析[17]。目前，EMD方法已經(jīng)廣泛應(yīng)用于天氣、地震和醫(yī)療等信號(hào)處理的各個(gè)領(lǐng)域。其基本思想是將一個(gè)不規(guī)則信號(hào)表示成若干IMF與單調(diào)殘差函數(shù)相加的形式，一維信號(hào)x(t)的分解可表示為

(1)

式中：imfi(t)為第i個(gè)IMF；n為IMF的個(gè)數(shù)；rn(t)為殘差函數(shù)。

EMD方法的本質(zhì)是通過(guò)數(shù)據(jù)的時(shí)間尺度特征來(lái)獲得本征波動(dòng)模式，進(jìn)行數(shù)據(jù)分解。可以形象地稱這種分解過(guò)程為“篩選(sifting)”過(guò)程，用x(t)代表原始數(shù)據(jù)信號(hào)，其分解步驟如下：

1)求取原始數(shù)據(jù)信號(hào)的極大值點(diǎn)和極小值點(diǎn)。

2)分別由極大值點(diǎn)和極小值點(diǎn)利用三次樣條插值函數(shù)擬合獲得上包絡(luò)線s1和下包絡(luò)線s2。

3)計(jì)算上、下包絡(luò)線的均值m1，m1=(s1+s2)/2，如圖1所示。

圖1 原始NDVI時(shí)間序列和上、下包絡(luò)線及其均值

4)將原始數(shù)據(jù)序列x(t)減去包絡(luò)線均值m1，得到一個(gè)新的數(shù)據(jù)序列h0，h0=x(t)-m1。

5)將新的序列重復(fù)步驟1)—4)進(jìn)行迭代，直到滿足迭代停止準(zhǔn)則，獲得第一個(gè)IMF，imf1=hk。迭代停止準(zhǔn)則計(jì)算公式為

(2)

當(dāng)sdT超過(guò)給定的門限閾值時(shí)，迭代停止，此時(shí)hk代表第k次迭代后的新序列。

圖2為第1次迭代后的結(jié)果。

圖2 第1次迭代后的上、下包絡(luò)線及其均值

從圖2中可以看出，包絡(luò)線起伏較大，但總體均值趨于0值。圖3為第4次迭代后的結(jié)果。

圖3 第4次迭代后的上、下包絡(luò)線及其均值

從圖3中可以看出，包絡(luò)線均值起伏不大，已基本上趨于0值。

6)從原始數(shù)據(jù)序列中減去獲取的第一個(gè)IMF，得到殘差r1，r1=x(t)-imf1。

7)如果r1不是單調(diào)函數(shù)，重復(fù)步驟1)—6)，獲取imf2，imf3，…，直到r1為單調(diào)函數(shù)結(jié)束。

EMD方法分解結(jié)果如圖4所示。從圖4可以看出，通過(guò)EMD方法首先得到的是原始數(shù)據(jù)信號(hào)中的高頻分量imf1，后續(xù)的imf2—imf6分量頻率逐漸降低，殘差為單調(diào)函數(shù)，表示平均趨勢(shì)。由于數(shù)據(jù)長(zhǎng)度的有限性，殘差中的極值點(diǎn)數(shù)量逐漸減少，使得原始數(shù)據(jù)信號(hào)經(jīng)過(guò)幾次提取后可以得到有限個(gè)IMF分量和最后的殘差分量RES。

圖4 EMD方法分解結(jié)果

一般趨勢(shì)可以用殘差表示，也可以用殘差疊加若干IMF得到。圖5所示的曲線是根據(jù)中國(guó)植被區(qū)劃(1980)圖[18]選取的2條曲線及其EMD方法分解的殘差曲線RES和imf6+RES的結(jié)果，一條曲線位置在廣東省南部的熱帶季雨林、雨林區(qū)域； 另一條在黃土高原的暖溫帶草原區(qū)域。

(a) 熱帶季雨林、雨林區(qū)域 (b) 黃土高原暖溫帶草原區(qū)域

從圖5可以看出，如果只是把殘差RES看做平均趨勢(shì)項(xiàng)，2006—2015年間NDVI相對(duì)呈現(xiàn)單調(diào)變化趨勢(shì)。如果把imf6+RES看做趨勢(shì)項(xiàng)，2006—2008年間NDVI會(huì)呈現(xiàn)分段變化，圖5中2區(qū)域都呈先升后降的現(xiàn)象，圖5(a)NDVI在2006—2012年間先升后降，在2012—2015年間緩慢上升； 圖5(b)NDVI在2006—2014年間一直呈現(xiàn)上升趨勢(shì)，在2014—2015年間呈下降態(tài)勢(shì)，而平均趨勢(shì)一個(gè)下降一個(gè)則是上升。

圖5中的NDVI上升或下降趨勢(shì)是目視判斷的，是否呈現(xiàn)顯著的單調(diào)增、單調(diào)降趨勢(shì)還是沒(méi)有顯著變化，需要采用顯著性檢驗(yàn)方法進(jìn)行檢驗(yàn)。

1.2 MK檢驗(yàn)方法

由于本文主要針對(duì)單調(diào)性的變化趨勢(shì)進(jìn)行研究，因此主要針對(duì)殘余量，即平均趨勢(shì)進(jìn)行趨勢(shì)檢驗(yàn)，檢驗(yàn)是顯著單調(diào)增、單調(diào)降還是沒(méi)有顯著變化。檢驗(yàn)方法采用MK趨勢(shì)檢驗(yàn)方法。MK檢驗(yàn)方法最初由Mann[19]和Kendall[20]提出，用來(lái)檢測(cè)水域降水的長(zhǎng)期變化趨勢(shì)和突變情況。在時(shí)間序列趨勢(shì)分析中，許多學(xué)者應(yīng)用MK檢驗(yàn)方法分析氣溫、水質(zhì)、降水和徑流等要素時(shí)間序列趨勢(shì)變化[21-24]。MK檢驗(yàn)方法不受少數(shù)異常值的干擾，也不需要樣本遵循一定的分布，適用于非正態(tài)分布的數(shù)據(jù)。在MK檢驗(yàn)中，設(shè)一時(shí)間序列數(shù)據(jù)(x1，x2,...,xn)是n個(gè)獨(dú)立的、隨機(jī)變量同分布的樣本； 檢驗(yàn)的統(tǒng)計(jì)量S計(jì)算公式為

(3)

(4)

當(dāng)n>10時(shí)，S近似正態(tài)分布，其均值為0，方差var(S)=n(n-1)(2n+5)/18，標(biāo)準(zhǔn)的正態(tài)系統(tǒng)變量ZMK計(jì)算公式為

(5)

這樣，在雙邊的趨勢(shì)檢驗(yàn)中，在α顯著水平上，如果|ZMK|>Z1-α/2，則原假設(shè)是不可接受的，即在α顯著水平上，時(shí)間序列數(shù)據(jù)存在明顯上升或下降趨勢(shì)，否則接受原始時(shí)間序列無(wú)變化趨勢(shì)的假設(shè)。對(duì)于統(tǒng)計(jì)量ZMK，當(dāng)ZMK>0時(shí)為上升趨勢(shì)； 當(dāng)ZMK<0時(shí)為下降趨勢(shì)。當(dāng)顯著性水平為α?xí)r，置信度為(1-α)×100%，本文要求95%的置信度，則ZMK=1.96。

2 實(shí)驗(yàn)與分析

2.1 實(shí)驗(yàn)區(qū)域及實(shí)驗(yàn)數(shù)據(jù)

為了檢驗(yàn)本文方法的魯棒性，根據(jù)氣候和緯度的不同，在中國(guó)北部、西部和南部的區(qū)域任意固定大小的3個(gè)樣區(qū)，通過(guò)實(shí)驗(yàn)結(jié)果與樣區(qū)實(shí)際情況的對(duì)比，分析評(píng)價(jià)該方法的適用性。樣區(qū)A位于榆林附近，是黃土高原與內(nèi)蒙古高原的過(guò)渡區(qū)。北部是毛烏素沙漠南緣風(fēng)沙草灘區(qū)，南部是黃土高原腹地，溝壑縱橫，丘陵峁梁交錯(cuò)，海拔為1 000～1 200 m。植被分區(qū)屬于溫帶草原區(qū)域和溫帶荒漠區(qū)域。氣溫四季變化明顯。經(jīng)過(guò)長(zhǎng)期退耕還林還草的治理，該區(qū)域如今沙灘地植被茂盛； 海子(湖泊)星羅棋布。水土侵蝕逐步得到治理，水土流失得到初步控制，生態(tài)環(huán)境有較大改善。樣區(qū)B位于川、青、藏三省區(qū)的結(jié)合部，大陸性季風(fēng)高原氣候，氣溫低，年均氣溫在0 ℃以下，日照時(shí)間長(zhǎng)，晝夜溫差大，年均降水量為400～700 mm，海拔在4 000 m以上，植被分區(qū)屬于青藏高原高寒植被區(qū)域，包括山地寒溫性針葉林區(qū)域及高寒灌叢草甸地帶。樣區(qū)C位于廣州附近，屬于東亞季風(fēng)區(qū)，從北向南分別為中亞熱帶、南亞熱帶和熱帶氣候，是中國(guó)光、熱和水資源最豐富的地區(qū)之一。該地區(qū)年平均降水量在1 300～2 500 mm之間，降水充沛，空間分布基本上也呈南高北低的趨勢(shì)。改革開(kāi)放后經(jīng)濟(jì)發(fā)展迅速。

3個(gè)樣區(qū)實(shí)驗(yàn)數(shù)據(jù)選取的是2006—2015年MODIS 16 d合成的NDVI數(shù)據(jù)。從http: //ladsweb.nascom.nasa.gov/data/search.html網(wǎng)站下載MODIS產(chǎn)品數(shù)據(jù)MOD13A2數(shù)據(jù)，即全球1 km空間分辨率植被指數(shù)16 d合成數(shù)據(jù)。MOD13A2數(shù)據(jù)16 d的組合，1 a有23個(gè)時(shí)相，共10 a的數(shù)據(jù)。

實(shí)驗(yàn)數(shù)據(jù)的預(yù)處理包括重投影、數(shù)據(jù)重排和數(shù)據(jù)重構(gòu)。重投影即將正弦投影轉(zhuǎn)換為經(jīng)緯度投影； 數(shù)據(jù)重排即將按空間順序組織的數(shù)據(jù)按照時(shí)間順序進(jìn)行了重排； 數(shù)據(jù)重構(gòu)采用Chen等[25]的方法,核心是通過(guò)Savitzky-Golay濾波器的迭代濾波使NDVI不斷逼近上包絡(luò)來(lái)減弱云噪聲和大氣變化的影響。

2.2 結(jié)果與分析

每個(gè)樣區(qū)的衛(wèi)星影像及經(jīng)過(guò)趨勢(shì)提取后的NDVI中值和NDVI趨勢(shì)如圖6。圖6中(c)，(f)和(i)為3樣區(qū)NDVI趨勢(shì)圖，深灰表示10 a間植被沒(méi)有顯著變化，淺灰表示NDVI值小于0.1，即植被很少，橘色表示10 a間植被有變褐的趨勢(shì)，綠色表示10 a間植被有變綠的趨勢(shì)。由圖6(a)—(c)可知，樣區(qū)A植被從北到南NDVI從低到高分布，即該地區(qū)植被的覆蓋程度從南向北變差，這與該地區(qū)的地貌地形和氣候是相吻合的，但從變化趨勢(shì)看，從南到北，植被總體上是呈現(xiàn)變綠的趨勢(shì)。這種現(xiàn)象與這10 a間退耕還林、水土侵蝕得到治理、水土流失得到控制、生態(tài)環(huán)境改善相吻合。由圖6(d)—(f)可知，樣區(qū)B植被的分布情況相對(duì)復(fù)雜，因地形復(fù)雜，有常年積雪區(qū)域存在，相應(yīng)植被變化的趨勢(shì)也比較復(fù)雜。在NDVI值比較低的周圍既有呈現(xiàn)“褐化”的趨勢(shì)，也有呈現(xiàn)“綠化”的趨勢(shì)，產(chǎn)生這種現(xiàn)象的原因需要后續(xù)進(jìn)行更深入細(xì)致的分析。由圖6(g)—(i)可知，樣區(qū)C除城市區(qū)域外，植被覆蓋度很高。這與該地區(qū)的氣候吻合。從植被變化的趨勢(shì)看，存在“褐化”的區(qū)域也存在“綠化”的區(qū)域?！昂只钡膮^(qū)域大多與城市擴(kuò)展、經(jīng)濟(jì)開(kāi)發(fā)相關(guān)聯(lián)。

(d) 樣區(qū)B 衛(wèi)星影像 (e) 樣區(qū)B NDVI中值 (f) 樣區(qū)B NDVI趨勢(shì)圖

(g) 樣區(qū)C 衛(wèi)星影像 (h) 樣區(qū)C NDVI中值 (i) 樣區(qū)C NDVI趨勢(shì)圖

3 結(jié)論

本文提出了一種利用EMD方法結(jié)合 MK檢驗(yàn)方法從NDVI序列檢測(cè)植被“綠化”或“褐化”變化趨勢(shì)的新方法，該方法在檢測(cè)變化趨勢(shì)時(shí)無(wú)需去除序列的周期性或季節(jié)性直接進(jìn)行趨勢(shì)分析，對(duì)序列不需要作任何隱含假設(shè)，且趨勢(shì)檢測(cè)過(guò)程中不需要使用最小二乘進(jìn)行一元回歸。本文通過(guò)對(duì)中國(guó)北部、西部和南部3個(gè)實(shí)驗(yàn)區(qū)域2006—2015年10 a間植被變化趨勢(shì)的分析，表明本文方法是一種有效的時(shí)間序列變化趨勢(shì)的分析方法，對(duì)生態(tài)長(zhǎng)期變化研究有一定參考價(jià)值。

本文只研究了10 a間NDVI變化的總體平均趨勢(shì)，事實(shí)上10 a間植被的變化不一定是單調(diào)變化的，可能會(huì)先上升后下降，也可能先下降后上升，也可能分段上升或分段下降，有些變化可以直接關(guān)聯(lián)到變化的驅(qū)動(dòng)力，有些變化還不能直接關(guān)聯(lián)到驅(qū)動(dòng)力，而僅憑NDVI單要素的時(shí)間序列變化趨勢(shì)分析生態(tài)驅(qū)動(dòng)機(jī)理尚有不確定性。后續(xù)的研究一方面將擴(kuò)展到對(duì)全國(guó)的植被“綠化”、“褐化”的趨勢(shì)進(jìn)行研究； 另一方面將結(jié)合數(shù)據(jù)的時(shí)間尺度特征通過(guò)EMD方法分解的本征模態(tài)的疊加，分析NDVI在10 a間的波動(dòng)模式，期望對(duì)生態(tài)變化的時(shí)間尺度特征有更清晰的了解，從而評(píng)價(jià)中國(guó)10 a間退耕還林、退耕還草所取得的實(shí)效，為陸地生態(tài)系統(tǒng)變化研究提供支撐。