基于Spearman 等級系數(shù)的植被變化趨勢分析

王佃來， 宿愛霞， 劉文萍

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

2.中國軟件評測中心, 北京100048

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

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

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

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

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

為了避免其他區(qū)域數(shù)據(jù)干擾而突出研究區(qū)域和數(shù)據(jù)配準(zhǔn)，對下載數(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)緯度范圍對1998—2013年的16 幅數(shù)據(jù)進(jìn)行裁剪.

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

1.2 研究方法

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

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

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

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

rs=1 ?((6×29)/(73?7))=0.482 .在秩次計算過程中，如果序列中存在相同的值，則秩次為所有相等數(shù)據(jù)秩次的平均值.例如，Data2 中的數(shù)據(jù)199、199，其秩次原本為2、3，因為要求其值相同，所以取其秩次的平均值2.5，類似的數(shù)據(jù)還有200、200.

檢驗rs的顯著性可以分為以下兩種情況：

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

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

1.2.2 算法的抗噪性分析

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

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

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

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

1.2.3 算法流程

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

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

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

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

步驟4將NDVI 序列Y 排序并計算其秩次

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

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

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

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

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

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

2.2 實(shí)驗結(jié)果

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

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

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

表2 Pearson 相關(guān)系數(shù)法、Mann-Kendall 檢驗法和Spearman 等級相關(guān)系數(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ū)的植被覆蓋變化以改善趨勢為主，植被明顯改善和輕微改善的區(qū)域達(dá)到內(nèi)蒙古自治區(qū)總面積的83.30%，植被退化的區(qū)域占總面積的16.63%，植被嚴(yán)重退化的區(qū)域僅為0.38%.

圖2 為基于Pearson 相關(guān)系數(shù)法生成的內(nèi)蒙古自治區(qū)植被覆蓋變化趨勢圖，圖2 中深綠色代表植被明顯改善，淺綠色代表植被輕微改善，淺褐色代表植被輕微退化，深褐色代表植被嚴(yán)重退化，黑色代表植被無變化趨勢，可以看出植被明顯改善的區(qū)域集中在阿拉善盟的阿拉善右旗的大部分地區(qū)、鄂爾多斯市的大部分地區(qū)、赤峰市的南部、通遼市的西南部、興安盟的部分地區(qū)、呼倫貝爾市的北部以及南部的扎蘭屯市.植被輕微退化的地區(qū)分布在阿拉善盟的額濟(jì)納旗的西北部、錫林郭勒盟的大部分地區(qū)以及呼倫貝爾市西部地區(qū).植被嚴(yán)重退化的區(qū)域主要分布在人類生活密集的市區(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ū)植被變化趨勢圖Figure 2 Trend of vegetation cover changes in the Inner Mongolia Autonomous Region based on Pearson correlation coefficient

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

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

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

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

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

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

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

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

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

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

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

2.3 實(shí)驗結(jié)果對比與討論

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

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

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

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

3 結(jié) 語

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

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

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

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