王佃來(lái)， 宿愛(ài)霞， 劉文萍

1.首鋼工學(xué)院信息工程系, 北京100144

2.中國(guó)軟件評(píng)測(cè)中心, 北京100048

3.北京林業(yè)大學(xué)信息學(xué)院, 北京100083

植被是陸地生態(tài)系統(tǒng)的重要組成部分，是對(duì)生態(tài)環(huán)境影響的敏感指標(biāo)[1-8]，因此基于長(zhǎng)時(shí)間序列遙感數(shù)據(jù)的植被變化監(jiān)測(cè)與評(píng)估一直是生態(tài)學(xué)和全球變化的重要研究領(lǐng)域.相關(guān)系數(shù)法被認(rèn)為是分析長(zhǎng)期植被趨勢(shì)變化的最優(yōu)方法[4,8]，通常包括Pearson 相關(guān)系數(shù)法和Spearman 等級(jí)相關(guān)系數(shù)法.部分學(xué)者以Pearson 相關(guān)系數(shù)法分析長(zhǎng)時(shí)間序列植被趨勢(shì)[4-5,8]，該方法屬于參數(shù)檢驗(yàn)的范疇[3,5]，但存在以下局限性：1）僅適用于兩個(gè)隨機(jī)變量服從二元正態(tài)分布的數(shù)據(jù)；2）對(duì)異常值敏感且只能發(fā)現(xiàn)變量之間的線性關(guān)系.Spearman 等級(jí)相關(guān)系數(shù)屬于非參數(shù)檢范疇，對(duì)變量數(shù)據(jù)無(wú)正態(tài)分布假設(shè)要求，既可以發(fā)現(xiàn)線性關(guān)系，又能發(fā)現(xiàn)單調(diào)的非線性關(guān)系，且對(duì)異常值不敏感，因此適用范圍更廣.長(zhǎng)時(shí)間序列的植被變化監(jiān)測(cè)與評(píng)估受以下因素制約：1）遙感數(shù)據(jù)長(zhǎng)度較短，通常不超過(guò)30；2）數(shù)據(jù)較難滿足正態(tài)分布假設(shè)，3）遙感數(shù)據(jù)普通存在較多噪聲，因此Spearman 等級(jí)相關(guān)系數(shù)更適合長(zhǎng)時(shí)序遙感植被監(jiān)測(cè).通過(guò)文獻(xiàn)檢索發(fā)現(xiàn)以Spearman 等級(jí)相關(guān)系數(shù)進(jìn)行長(zhǎng)時(shí)序植被變化趨勢(shì)分析的報(bào)道很少，于是本文基于1998—2013年SPOT vegetaiton 的植被指數(shù)(normalized vegetation index, NDVI)數(shù)據(jù)，使用Spearman 等級(jí)相關(guān)系數(shù)法分析了內(nèi)蒙古自治區(qū)近16年的植被變化趨勢(shì)，并將分析結(jié)果與Pearson 相關(guān)系數(shù)法和Mann-Kendall 假設(shè)檢驗(yàn)法進(jìn)行對(duì)比，驗(yàn)證了該方法在長(zhǎng)時(shí)序植被監(jiān)測(cè)中的可行性和適用性.

1 研究數(shù)據(jù)與研究方法

1.1 研究數(shù)據(jù)

本文所用的遙感數(shù)據(jù)是1998—2013年SPOT-4 vegetation 空間分辨率為1 km 的S10 NDVI 數(shù)據(jù)，該數(shù)據(jù)每年12 個(gè)月，每月3 組數(shù)據(jù)(上旬、中旬和下旬的數(shù)據(jù))，即每年36 組數(shù)據(jù).該數(shù)據(jù)下載網(wǎng)址為：http://www.vgt.vito.be/, 數(shù)據(jù)所在區(qū)域?yàn)闁|南亞(SE-Asia).

為了避免其他區(qū)域數(shù)據(jù)干擾而突出研究區(qū)域和數(shù)據(jù)配準(zhǔn)，對(duì)下載數(shù)據(jù)進(jìn)行以下預(yù)處理：一是數(shù)據(jù)裁剪，使用VGT Extract 軟件按內(nèi)蒙古自治區(qū)的經(jīng)緯度邊界范圍進(jìn)行數(shù)據(jù)裁剪，其中經(jīng)度范圍為東經(jīng)97.160?～126.090?，緯度范圍為北緯37.380?～153.400?.二是數(shù)據(jù)掩模，根據(jù)內(nèi)蒙古自治區(qū)的行政區(qū)劃矢量圖使用ENVI 軟件的掩模功能裁剪出內(nèi)蒙古自治區(qū)的精確遙感數(shù)據(jù).數(shù)據(jù)的裁剪是按照同一經(jīng)緯度范圍對(duì)1998—2013年的16 幅數(shù)據(jù)進(jìn)行裁剪.

研究植被覆蓋變化趨勢(shì)時(shí)，為了取得更好的效果，本文數(shù)據(jù)選自一年中植被最茂盛時(shí)期的NDVI 值，即在每年36 組NDVI 數(shù)據(jù)中選取最大值，最后生成一幅圖像用于植被變化趨勢(shì)分析.

1.2 研究方法

1.2.1 Spearman 等級(jí)相關(guān)系數(shù)法介紹

Spearman 等級(jí)相關(guān)系數(shù)法是一種非參數(shù)檢驗(yàn)方法，可以度量變量之間的強(qiáng)弱關(guān)系，通常用rs表示.將數(shù)據(jù)xi、yi從小到大排序，記為xi、yi在排序后的位置，稱(chēng)為xi、yi的秩次，則秩次差Spearman 等級(jí)相關(guān)系數(shù)為

式中，n 為時(shí)間序列的長(zhǎng)度.為便于理解，給出如表1 所示的Spearman 等級(jí)相關(guān)系數(shù)的算例.

表1 Spearman 等級(jí)相關(guān)系數(shù)算例Table 1 Example for Spearman rank correlation coefficient

rs=1 ?((6×29)/(73?7))=0.482 .在秩次計(jì)算過(guò)程中，如果序列中存在相同的值，則秩次為所有相等數(shù)據(jù)秩次的平均值.例如，Data2 中的數(shù)據(jù)199、199，其秩次原本為2、3，因?yàn)橐笃渲迪嗤?，所以取其秩次的平均?.5，類(lèi)似的數(shù)據(jù)還有200、200.

檢驗(yàn)rs的顯著性可以分為以下兩種情況：

1）當(dāng)樣本數(shù)n>50 時(shí)，以t 檢驗(yàn)rs的顯著性，其計(jì)算公式為

式中，n 為樣本數(shù)量，rs為Spearman 等級(jí)系數(shù).在顯著水平α 條件下，比較式(2)的|t|值與t值表中查得的自由度為n ?2 的臨界值t'，決定rs是否顯著.

1.2.2 算法的抗噪性分析

因?yàn)檫b感數(shù)據(jù)中普遍存在噪聲，所以將抗噪性作為算法的一個(gè)重要指標(biāo).從理論上來(lái)看，Spearman 等級(jí)相關(guān)系數(shù)是通過(guò)數(shù)據(jù)序列排序后計(jì)算數(shù)據(jù)的秩次來(lái)實(shí)現(xiàn)的，因此較小的噪聲對(duì)數(shù)據(jù)的排序結(jié)果產(chǎn)生的影響較小，可見(jiàn)該方法具有一定的抗噪性.

為了進(jìn)一步驗(yàn)證這一結(jié)論，本文隨機(jī)選取研究區(qū)30 組數(shù)據(jù)作為原始數(shù)據(jù)，以每組數(shù)據(jù)極差的10%為隨機(jī)噪聲添加到原始數(shù)據(jù)上生成對(duì)比數(shù)據(jù)，實(shí)驗(yàn)結(jié)果如圖1 所示.

圖1 Spearman等級(jí)相關(guān)系數(shù)法的原始數(shù)據(jù)與添加噪聲數(shù)據(jù)的對(duì)比圖Figure 1 Comparison between original data and added noise data based on Spearman rank correlation coefficient method

從圖1 中可以看出，原始數(shù)據(jù)與添加噪聲數(shù)據(jù)的計(jì)算結(jié)果吻合度很高，只是在個(gè)別點(diǎn)上有較小的偏移，表明Spearman 等級(jí)相關(guān)系數(shù)法有較好的抗噪性.

1.2.3 算法流程

基于Spearman 等級(jí)相關(guān)系數(shù)法的植被變化趨勢(shì)分析算法流程如下：

步驟1讀取多年經(jīng)過(guò)幾何處理的n 幅NDVI 圖像數(shù)據(jù)，并存儲(chǔ)該數(shù)據(jù).

步驟2在每幅圖像對(duì)應(yīng)的相同位置取NDVI 值，形成NDVI 序列Y = {y1,y2,··· ,yn}；按時(shí)間形成時(shí)間序列X ={x1,x2,··· ,xn}.

步驟3將時(shí)間序列X 排序并計(jì)算其秩次如果按年份從小到大排列，則的秩次分別為1,2,··· ,n.

步驟4將NDVI 序列Y 排序并計(jì)算其秩次

步驟5基于式(1)計(jì)算Spearman 等級(jí)相關(guān)系數(shù)rs，其中d=X'?Y '.

步驟6從Spearman 秩相關(guān)系數(shù)界值表查得臨界值t(n50)，也可以根據(jù)式(2)計(jì)算t(n>50).如果rs>0 且rs>t，則表明植被明顯改善；如果rs>0 且rst，則表明植被輕微改善；如果rs= 0，則表明植被無(wú)變化；如果rs<0 且rs> t，則表明植被嚴(yán)重退化；如果rs<0 且rst，則表明植被輕微退化.

步驟7重復(fù)步驟2～6，直至計(jì)算NDVI 圖像的所有序列，即可得到多年某區(qū)域的植被變化趨勢(shì).

2 實(shí)驗(yàn)結(jié)果與討論

2.1 實(shí)驗(yàn)環(huán)境與實(shí)驗(yàn)工具

實(shí)驗(yàn)的開(kāi)發(fā)環(huán)境如下：操作系統(tǒng)為Window8.1，開(kāi)發(fā)工具為Eclipse4.5.1 和JDK1.7，遙感圖像裁剪軟件VGTExtract，遙感圖像掩模軟件為ENVI4.5，開(kāi)發(fā)語(yǔ)言為Java.

2.2 實(shí)驗(yàn)結(jié)果

為了驗(yàn)證Spearman 等級(jí)相關(guān)系數(shù)法在植被變化趨勢(shì)分析中的適用性，本文引入植被變化趨勢(shì)分析常用的Pearson 相關(guān)系數(shù)法和Mann-Kendall 趨勢(shì)分析算法對(duì)1998—2013年SPOT-4 vegetation 空間分辨率為1 km 的S10 NDVI 數(shù)據(jù)進(jìn)行分析，得出了相應(yīng)的實(shí)驗(yàn)結(jié)果.

2.2.1 基于Pearson 相關(guān)系數(shù)的植被變化趨勢(shì)分析結(jié)果

利用Pearson 相關(guān)系數(shù)對(duì)內(nèi)蒙古自治區(qū)進(jìn)行長(zhǎng)時(shí)序植被覆蓋變化分析，所得結(jié)果如表2所示，其中rp是對(duì)內(nèi)蒙古自治區(qū)1998—2013年際每個(gè)配準(zhǔn)的NDVI 序列值與年份序列值進(jìn)行相關(guān)分析得到的Pearson 相關(guān)系數(shù).在顯著水平α = 0.05 條件下，從相關(guān)系數(shù)顯著性檢驗(yàn)表查得相關(guān)系數(shù)值為0.497，其中自由度為n?2=16?2=14.

表2 Pearson 相關(guān)系數(shù)法、Mann-Kendall 檢驗(yàn)法和Spearman 等級(jí)相關(guān)系數(shù)法的趨勢(shì)分析結(jié)果Table 2 Results of vegetation cover changes based on Pearson correlation coefficient,Mann-Kendall test and Spearman rank correlation coefficient

由表2 數(shù)據(jù)可以看出：在1998—2013年期間，內(nèi)蒙古自治區(qū)的植被覆蓋變化以改善趨勢(shì)為主，植被明顯改善和輕微改善的區(qū)域達(dá)到內(nèi)蒙古自治區(qū)總面積的83.30%，植被退化的區(qū)域占總面積的16.63%，植被嚴(yán)重退化的區(qū)域僅為0.38%.

圖2 為基于Pearson 相關(guān)系數(shù)法生成的內(nèi)蒙古自治區(qū)植被覆蓋變化趨勢(shì)圖，圖2 中深綠色代表植被明顯改善，淺綠色代表植被輕微改善，淺褐色代表植被輕微退化，深褐色代表植被嚴(yán)重退化，黑色代表植被無(wú)變化趨勢(shì)，可以看出植被明顯改善的區(qū)域集中在阿拉善盟的阿拉善右旗的大部分地區(qū)、鄂爾多斯市的大部分地區(qū)、赤峰市的南部、通遼市的西南部、興安盟的部分地區(qū)、呼倫貝爾市的北部以及南部的扎蘭屯市.植被輕微退化的地區(qū)分布在阿拉善盟的額濟(jì)納旗的西北部、錫林郭勒盟的大部分地區(qū)以及呼倫貝爾市西部地區(qū).植被嚴(yán)重退化的區(qū)域主要分布在人類(lèi)生活密集的市區(qū)，如通遼市區(qū)、赤峰市區(qū)、烏蘭察布市區(qū)、呼和浩特市區(qū)、包頭市區(qū)、巴彥淖爾市區(qū)、鄂爾多期市區(qū)、烏海市區(qū).植被輕微改善的區(qū)遍布內(nèi)蒙古自治區(qū)的大部分地區(qū)，占內(nèi)蒙古自治區(qū)總面積的60.03%，如表2 所示.

圖2 Pearson 相關(guān)系數(shù)法所得內(nèi)蒙古自治區(qū)植被變化趨勢(shì)圖Figure 2 Trend of vegetation cover changes in the Inner Mongolia Autonomous Region based on Pearson correlation coefficient

2.2.2 基于Mann-Kendall 檢驗(yàn)法的植被變化趨勢(shì)分析結(jié)果

以Mann-Kendall 檢驗(yàn)法分析內(nèi)蒙古自治區(qū)的植被變化趨勢(shì)，所得結(jié)果如表2 所示，其中植被明顯改善的區(qū)域?yàn)?9.27%，植被輕微改善的區(qū)域?yàn)?3.39%，植被嚴(yán)重退化的區(qū)域?yàn)?.30%，植被輕微退化的區(qū)域?yàn)?5.36%.從表2 中可以看出，內(nèi)蒙古植被變化以改善趨勢(shì)為主，改善區(qū)域占總面積的82.66%，退化區(qū)域占總面積的15.66%.

利用Mann-Kendall 檢驗(yàn)法生成的內(nèi)蒙古自治區(qū)植被覆蓋變化結(jié)果如圖3 所示，圖中顏色標(biāo)注與圖2 相同，可以看出整體的植被變化分布與圖2 基本一致.

圖3 Mann-Kendall 檢驗(yàn)法所得內(nèi)蒙古自治區(qū)植被變化趨勢(shì)圖Figure 3 Trend of vegetation cover changes in the Inner Mongolia Autonomous Region based on Mann-Kendall test

2.2.3 基于Spearman 等級(jí)相關(guān)系數(shù)法的植被變化趨勢(shì)分析結(jié)果

利用基于Spearman 等級(jí)相關(guān)系數(shù)對(duì)內(nèi)蒙古自治區(qū)進(jìn)行植被變化趨勢(shì)分析的結(jié)果如表2所示，其中rs為Spearman 等級(jí)相關(guān)系數(shù).在顯著水平α=0.05 條件下，由Spearman 秩相關(guān)系數(shù)界值表查得臨界值為0.538，其中自由度為n ?2=16 ?2=14.

從表2 中可以看出內(nèi)蒙古自治區(qū)植被變化活動(dòng)在增強(qiáng)，其中明顯改善和輕微改善的區(qū)域占總面積的84.58%，植被退化區(qū)域占總面積的15.33%.

圖4 為基于Spearman 等級(jí)相關(guān)系數(shù)法生成的內(nèi)蒙古自治區(qū)植被覆蓋變化趨勢(shì)分布圖，與圖2 和3 在整體上有較高的一致性，植被明顯改善、輕微退化和嚴(yán)重退化的區(qū)域分布大致相同.

圖4 Spearman 等級(jí)相關(guān)系數(shù)法所得內(nèi)蒙古自治區(qū)植被變化趨勢(shì)圖Figure 4 Trend of vegetation cover changes in the Inner Mongolia Autonomous Region based on Spearman rank correlation coefficient

2.2.4 3 種方法在研究區(qū)域的時(shí)間復(fù)雜度對(duì)比

為直觀了解3 種植被趨勢(shì)分析方法的時(shí)間復(fù)雜度，本文比較了研究區(qū)域3 241×1 795 組數(shù)據(jù)的運(yùn)行時(shí)間，結(jié)果見(jiàn)表3.實(shí)驗(yàn)結(jié)果顯示：從運(yùn)行時(shí)間來(lái)看，Pearson 相關(guān)系數(shù)法最長(zhǎng)，Spearman 等級(jí)相關(guān)系數(shù)法次之，Mann-Kendall 檢驗(yàn)法最短.從時(shí)間復(fù)雜度層面來(lái)分析，Spearman 等級(jí)相關(guān)系數(shù)法需要對(duì)數(shù)據(jù)進(jìn)行二次排序，其時(shí)間復(fù)雜度在理論上高于Pearson 相關(guān)系數(shù)法，但實(shí)驗(yàn)結(jié)果恰恰相反.究其原因如下：Spearman 等級(jí)相關(guān)系數(shù)法雖然進(jìn)行二次排序，但數(shù)據(jù)組的長(zhǎng)度較短，計(jì)算機(jī)排序速度較快，且該算法主要執(zhí)行加減運(yùn)算；Pearson 相關(guān)系數(shù)法要進(jìn)行浮點(diǎn)型運(yùn)算，因此運(yùn)行時(shí)間較長(zhǎng).

2.3 實(shí)驗(yàn)結(jié)果對(duì)比與討論

從表2 中可以看出，基于Spearman 等級(jí)相關(guān)系數(shù)法的植被覆蓋變化趨勢(shì)分析結(jié)果與Pearson 相關(guān)系數(shù)法的分析結(jié)果基本一致.若采用Spearman 等級(jí)相關(guān)系數(shù)法，則植被明顯改善和輕微改善區(qū)域所占比例的和為84.58%；若采用Pearson 相關(guān)系數(shù)法，則植被明顯改善和輕微改善區(qū)域所占比例的和為83.31%.兩種方法在植被改善區(qū)域的差值僅為1.27%；同樣在植被退化區(qū)域也只有細(xì)微差別，僅為1.30%；在無(wú)趨勢(shì)區(qū)域幾乎無(wú)差別，僅為0.03%.

表3 3種植被趨勢(shì)分析方法的運(yùn)行時(shí)間對(duì)比表Table 3 Comparison of runtime to three trend analysis methods

分析表2 可以發(fā)現(xiàn)，Spearman 等級(jí)相關(guān)系數(shù)法與Mann-Kendall 法的植被覆蓋變化趨勢(shì)分析結(jié)果基本一致，兩者在植被改善區(qū)域的差為1.44%，在植被退化區(qū)域的差異更少，僅為0.33%.3 種方法的統(tǒng)計(jì)量都是在顯著水平α = 0.05的條件下計(jì)算的，且3 種方法結(jié)果差異均小于5%，因此有理由認(rèn)為這3 種方法對(duì)內(nèi)蒙古自治區(qū)的植被覆蓋變化趨勢(shì)分析結(jié)果的差異不明顯，說(shuō)明Spearman 等級(jí)相關(guān)系法是有效的.

觀察圖2～4 可以看出，3 種方法的植被覆蓋變化趨勢(shì)分析結(jié)果在空間分布上有較好的一致性，植被明顯改善的區(qū)域集中在阿拉善右旗的大部分地區(qū)、鄂爾多斯市除環(huán)城的大部分地區(qū)、赤峰市的南部、通遼市的西南部、興安盟的部分地區(qū)、呼倫貝爾市的北部和南部地區(qū)；植被嚴(yán)重退化的區(qū)域主要集中在人類(lèi)活動(dòng)密集的城區(qū).植被輕微退化的區(qū)域在圖2～4 中的分布也基本一致，因此Spearman 等級(jí)相關(guān)系數(shù)法有較高的正確性和較好的適用性.

3 結(jié) 語(yǔ)

本文基于SPOT VEGATATION 歸一化植被指數(shù)數(shù)據(jù)研究了Spearman 等級(jí)相關(guān)系數(shù)法在長(zhǎng)時(shí)間序列植被覆蓋變化趨勢(shì)分析中的可行性，并分析了算法的抗噪性，同時(shí)將趨勢(shì)分析結(jié)果與Pearson 相關(guān)系數(shù)法和Mann-Kendall 檢驗(yàn)實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，得出以下3條主要結(jié)論.

1) Spearman 等級(jí)相關(guān)系數(shù)法在30 組10% 的隨機(jī)噪聲數(shù)據(jù)模擬條件下具有良好的抗噪性.

2) Spearman 等級(jí)相關(guān)系數(shù)法、Pearson 相關(guān)系數(shù)法、Mann-Kendall 檢驗(yàn)法在內(nèi)蒙古自治區(qū)的植被覆蓋變化趨勢(shì)分析中均表現(xiàn)出良好的一致性：首先，植被覆蓋變化空間分布與對(duì)比算法有較高的一致性；其次，植被改善、植被不變和退化區(qū)域數(shù)據(jù)與對(duì)比算法的最大差異不超過(guò)2%.Spearman 等級(jí)相關(guān)系數(shù)法的特點(diǎn)一是對(duì)數(shù)據(jù)分布無(wú)要求，二是能發(fā)現(xiàn)線性及非線性關(guān)系，三是對(duì)異常值不敏感，因此該方法是一種很有潛力的長(zhǎng)時(shí)序遙感植被覆蓋變化趨勢(shì)分析方法.

3) 3 種算法的實(shí)驗(yàn)結(jié)果表明：在1998—2013年期間，內(nèi)蒙古自治區(qū)的植被覆蓋變化趨勢(shì)整體表現(xiàn)為增強(qiáng).植被明顯改善和輕微改善的區(qū)域占內(nèi)蒙古自治區(qū)總面積的83%以上，植被明顯改善的區(qū)域包括阿拉善右旗的大部分地區(qū)、鄂爾多斯市除環(huán)城的大部分地區(qū)、赤峰市的南部、通遼市的西南部、興安盟的部分地區(qū)、呼倫貝爾市的北部和南部地區(qū)，植被嚴(yán)重退化的區(qū)域集中在人類(lèi)活動(dòng)頻繁的城市區(qū)域.