基于地統計與空間自相關的山區草地資源空間分布格局及變化特征分析

孫曉莉，袁 磊，趙 然，和萬榮，田淑靜，郭朝元

(1.云南國土資源職業學院測繪地理信息學院,云南 昆明 652501; 2.云南師范大學信息學院, 云南 昆明 650500；3. 徐州工程學院管理學院，江蘇 徐州 221000)

【研究意義】草地資源是陸地生態系統的重要組成部分，在水土流失防治、生態環境保護等方面起著舉足輕重的作用[1-2]。草地資源不僅是發展國民經濟的物質基礎，還具有維護生態環境的天然屏障功能[3]。近年來，受經濟發展、人口增加等因素影響，中國草地資源退化嚴重[4]，草地資源面積減少。但隨著退耕還草、退林還草等草地保護政策的實施，部分區域草地資源面積出現了增加趨勢，針對草地資源的這一變化情況，國內外學者開展了大量研究工作。【前人研究進展】主要集中在草地資源類型退化時空變化特征及其原因分析[5-8]、微生物對草地資源數量、質量的影響[9-10]、基于遙感數據的生物量反演[11]、草地植被覆蓋變化[12]等方面，傳統的研究側重于草地資源類型、草地資源質量等方面的較多。【本研究切入點】基于地統計和空間自相關模型從縣級、鄉(鎮)級、村級區域單元3種不同尺度定量地對山區草地資源空間分布格局及變化特征研究較少。【擬解決的關鍵問題】以云南省會澤縣為研究區域，從縣級、鄉(鎮)級、村級區域單元不同尺度，基于OpenGeoda、ArcGIS等軟件，從草地資源動態變化度、草地資源集中度指數、核密度分析、空間自相關模型對2009-2016年草地資源空間分布格局及變化特征進行分析，為草地資源清查、生態環境保護以及全面深化草原生態文明體制改革提供數據支撐和決策依據。

1 材料與方法

1.1 研究區概況

會澤縣位于云南省東北部、金沙江東岸、曲靖市西北部，東經103°03′～103°55′、北緯25°48′～27°04′，地處滇東高原與黔西高原結合部，位于烏蒙山系主峰地段，地勢南高北低，呈階梯狀下降，境內最高海拔大海梁子牯牛寨4017.3 m，最低海拔小江四級電站695 m，相對高差3322.3 m。境內草地資源豐富，主要包括山地草甸類、灌木叢類等。會澤縣共有21個鄉(鎮)，2009年土地總面積為588 912.70 hm2，其中草地資源75 928.38 hm2，占研究區土地總面積的12.89 %。

1.2 數據來源與處理

數據來源：2009、2016年草地資源數據來源于遙感影像，圖像空間分辨率為0.5 m，主要采用目視解譯、監督分類和外業核查的方式獲取草地資源圖斑；行政區劃等數據來源于會澤縣第二次全國土地調查數據庫。文中主要基于ArcGIS10.0和OpenGeoda軟件對數據進行處理為下文研究提供數據基礎。

圖1 研究區位置示意圖Fig.1 The location of research area

1.3 研究方法

1.3.1 草地資源動態變化度 草地資源動態變化度用于定量描述研究區域草地資源面積數量和空間上的變化情況以及在一定的監測時間段內變化的劇烈程度。草地資源動態變化度計算模型如下[13]：

式中：K為草地資源動態變化度，Sb為某區域監測結束時期的草地資源面積，Sa為某區域監測開始時期草地資源面積，T為監測時間段。

1.3.2 草地資源集中度指數 草地資源集中度指數的概念來源于地理集中度指數[14]，是指草地資源在研究區域空間地域上集中程度的測算指標，主要用來分析草地資源在地域空間上的分布情況。草地資源集中度指數計算模型如下[15-16]：

式中：G為草地資源集中度指數，xi為空間區域草地資源分布面積，T為研究區草地資源總面積，n為研究區空間區域個數。

1.3.3 核密度分析 核密度分析屬于一種非參數密度估計的統計方法，其是熱點和冷區識別與分析的一種有益的探測性方法[17]。在草地資源空間分布格局及變化特征分析中，核密度分析主要是通過輸入研究區域的草地資源數據來測算整個研究區域草地資源圖斑的集聚情況。假定x1,x2,x3.....xn為獨立同分布f的n個樣本，設其概率密度函數為f，則核密度分析的計算模型如下[18]：

式中：fn為草地資源地塊核密度計算值，n為觀測草地資源地塊的數量，k為核密度函數，x～xi為測算草地資源地塊x到樣本草地資源地塊xi處的距離，h為核密度計算帶寬或平滑參數。

1.3.4 空間自相關 空間自相關是指一些變量在同一個分布區內的觀測數據之間潛在的相互依賴性。而Moran’sI是衡量空間自相關集聚程度的基本測度，其包括全局空間自相關和局部空間自相關，在進行空間自相關分析時應首先確定空間權重矩陣。

(1)全局空間自相關：全局空間自相關(Global Moran’sI)是全面測度研究區域地理要素屬性值聚集或離散的程度。Global Moran’sI計算模型如下[14]：

標準化Z值經常用于檢查全局空間自相關指數的顯著性水平，其計算模型如下：

式中：E(I)為自相關性的期望值，VAR(I)為方差，Zscore表示標準化統計量的閾值。

主要選取Global Moran’sI來測算各村級區域單元草地資源面積比例(即草地資源面積占所在村級區域單元面積的比例)的全局空間自相關程度。

(2)局部空間自相關：局部空間自相關(Local Moran’sI)指數可以更好的分析研究區域內各村級區域單元在空間位置上與各自鄰近單元屬性的相關性，其可以使用圖形的形式更直觀地表示研究區域各村級區域單元的空間集聚情況，并根據LISA聚集圖分析研究區域村級區域單元的聚集或離散的空間位置分布，LISA圖的結果常分為HH(高-高型)、LL(低-低型)、HL(高-低型)、LH(低-高型)以及NN(非顯著相關型)，其中HH和LL型均表現為空間正相關類型，HL和LH型均表現為空間負相關型，NN則表示無明顯的聚集或離散特征。局部空間自相關計算模型如下[19-20]：

主要選取局部 Moran’sI來測算各村級區域單元草地資源面積比例(即草地資源面積占所在村級區域單元面積的比例)的局部空間自相關程度。若Moran’sI>0,則說明草地分布存在空間集聚，差異較小；反之，若Moran’sI<0,則說明草地分布存在空間顯著性差異。

2 結果與分析

2.1 草地資源區域分布及變化特征

2.1.1 草地資源面積變化分析 會澤縣草地資源

圖2 研究區草地資源集中度指數雷達圖Fig.2 Radar map of the research area grassland resources concentration index

動態變化度在2009-2016年全縣草地資源數量呈現增加趨勢，草地資源動態變化度為3.98 %，共增加21 134.51 hm2。研究區域各個鄉(鎮)草地資源動態變化度最高是田壩鄉為20.01 %、其次是雨碌鄉(14.24 %)、新街鄉(13.41 %)、紙廠鄉(12.09 %)、大井鎮(11.31 %)。全縣21個鄉(鎮)中，2009-2016年間草地資源增加的鄉(鎮)有16個，減少的有5個，減少的鄉(鎮)中動態變化度較大的是馬路鄉(-2.93 %)、礦山鎮(-2.18 %)和老廠鄉(-1.85 %)。各個鄉(鎮)草地資源面積的變化有一定的差異。

2.1.2 草地資源空間分布變化特征 2009年草地資源分布較多的鄉(鎮)為中西部的大海鄉、待補鎮、金鐘鎮，北部的迤車鎮和火紅鄉，南部的田壩鄉，而中部的者海鎮、大井鎮、新街鄉則分布的較少，整體呈現出“南北兩端高、中間低”的空間分布特征。2016年草地資源分布較多的鄉(鎮)為大海鄉、待補鎮、田壩鄉、金鐘鎮等，而位于研究區域中上方的老廠鄉、大橋鄉、樂業鎮、礦山鎮則分布較少，與2009年相比，整體呈現“南北兩端高、中間低”的空間分布特征更為明顯。

如圖2所示，全縣草地資源集中度2009年為7.70，2016年為7.98，2009-2016年增加了0.28。2016與2009年相比大多數鄉(鎮)草地資源集中度指數變動幅度較小，其中增加較明顯的鄉(鎮)主要有田壩鄉、老廠鄉、礦山鎮和娜姑鎮等；減少較為明顯的鄉(鎮)主要有駕車鄉、者海鎮、馬路鄉和五星鄉等。從2009和2016年草地資源集中度指數空間分布來看，高值區主要集中在北部的紙廠鄉、西北部的五星鄉和大橋鄉、中部的新街鄉、東南部的魯納鄉以及東北部的礦山鎮區域，整體上呈現出“西北、中東部區域高，其他區域稍低”的空間分布特征。

2.2 草地資源核密度空間分布及變化特征

借助ArcGIS10.0軟件，首先將2009和2016年的草地資源圖斑通過“要素轉點”工具，將草地資源圖斑(面數據)轉換為點數據，然后使用“核密度分析”功能進行草地資源核密度分析，并將結果采用自然間距分類法劃分為不同的區間，如圖3所示。2009年草地資源地塊核密度計算結果，共劃分了0～2、2～4、4～6.52塊/km23個不同的等級區間。2016年草地資源地塊核密度計算結果，共分為了0～2、2～4、4～5.11塊/km23個不同的等級區間。2009、2016年研究區草地資源空間分布核密度測算最大值分別為6.52和5.11塊/km2。從核密度計算結果來看，部分區域在2009-2016年期間草地資源地塊空間分布數量呈現減少的趨勢。2009年核密度計算結果高值區域集中在大橋鄉的中西部、五星鄉的東部、雨碌鄉的西南部以及魯納鄉的西部，2016年高值區域則集中在五星鄉的東部。2009-2016年核密度計算高值區域發生明顯的減少；研究區域草地資源核密度分析中值區域(2～4塊/km2)：2009年主要分布在除研究區兩側區域外的其他區域，主要集中在西北部的大橋鄉、五星鄉、樂業鎮，中部區域的新街鄉、雨碌鄉和魯納鄉，南部區域的田壩鄉等區域；2016年主要集中在五星鄉的竹箐村、鉛廠村、大坪子村、紅石巖村、野豬沖村、干松林村、石龍村，少量分布在紙廠鄉的紙廠村、燈草塘村和渾水塘村，魯納鄉的窩坡村和哈克村，田壩鄉的白巖村、田壩村等區域。中值區域，除五星鄉和田壩鄉變化不大外，大橋鄉、樂業鎮、雨碌鄉、新街鄉、上村鄉等發生了不同程度的減少；低值區域(0～2塊/km2)，2009-2016年出現了由外圍向中心擴張的趨勢，但整體變化不太明顯。

圖3 研究區草地資源核密度計算空間分布Fig.3 Spatial distribution of grassland resources kernel density of reseatch area

2016與2009年相比草地資源地塊分布密度較為分散，除五星鄉外，無明顯集中區域，且研究區兩側草地資源核密度明顯低于其他區域，整體呈現出以五星鄉為中心的西北部區域較集中，其他區域較分散的情況。從草地資源圖斑分布情況來看，2009-2016年期間草地資源面積增加一部分來源于退耕還草和退林還草，空間分布相對較為分散，且增加的草地資源圖斑約45 %分布于坡度25°以上的區域；草地資源面積增加僅僅是單個圖斑面積增加，草地圖斑數量卻在減少，導致2009-2016年草地資源核密度測算值出現了減少的情況。

2.3 草地資源分布空間關聯性分析

2.3.1 全局空間自相關 借助ArcGIS10.0軟件，首先計算出各個村級區域單元的草地資源面積比例，然后再基于相同字段使用空間連接進行關聯，分別得出2009和2016年草地資源面積比例圖層數據；然后再使用OpenGeoda軟件對研究區2009和2016年各村級區域單元草地資源面積比例值進行全局Moran’sI指數，并對其進行顯著性檢驗，結果如表1所示。2009和2016年研究區各村級區域單元草地資源面積比例值全局Moran’sI指數均>0且對比標準化Zscore(2.58)與顯著水平的臨界值，2009和2016年的Zscore分別為13.63、12.88，均明顯大于檢驗臨界值，并通過了顯著性水平a=0.05的檢驗，充分說明研究區域各村級區域單元的草地資源面積比例值呈現出顯著的空間自相關性，表現出顯著的集聚分布情況，即研究區各村級區域單元草地資源面積比例值高的村級區域單元其周圍區域草地資源面積比例值也高，反之，各村級區域單元草地資源面積比例值低的村級區域單元其周圍區域草地資源面積比例值也低。

表1 研究區2009和2016年各村級區域單元草地資源面積比例值全局Moran’s I指數Table 1 Global Moran's I index of proportion of grassland resources area at village level in 2009 and 2016

圖4 2009和2016年各村級區域單元草地資源面積比例值局部Moran’s I散點圖Fig.4 Scatter-plots of local Moran’s I for the proportion of grassland resources area at village level in 2009 and 2016

2.3.2 局部空間自相關 研究區各村級區域單元草地資源面積比例值的全局空間自相關分析結果顯示2009和2016年草地資源雖然在空間分布上呈現出集聚狀態，但卻出現了減弱的趨勢。從草地資源動態變化度、草地資源集中度指數、草地資源核密度等空間分布情況來看，各個鄉(鎮)內部存在較大的差異性。因此，通過局部空間自相關分析來研究區域內516個村級區域單元草地資源面積比例的空間自相關水平的差異性顯得尤為重要。

首先利用OpenGeoda軟件獲取研究區2009和2016年各村級區域單元草地資源面積比例的局部Moran’sI指數，然后以其為橫軸，各村級區域單元草地資源面積比例的局部Moran’sI指數相對應的空間滯后向量值為縱軸，繪制出研究區2009和2016年各村級區域單元草地資源面積比例觀測值的局部Moran’sI指數散點圖，如圖4所示。

從圖4可以得出，如果橫坐標>0，則說明空間單元為草地資源面積比例觀測值較高的村級區域單元，反之則為較低的村級區域單元；如果縱坐標>0，則空間單元的鄰近單元為草地資源面積比例觀測值較高的村級區域單元，反之則為較低的村級區域單元。根據圖4的局部Moran’sI散點圖結果，將其劃分為四個象限。其中第I象限為正相關的“高-高”類型區、第Ⅱ象限為負相關的“低-高”類型區、第Ⅲ象限為正相關的“低-低”類型區、第Ⅳ象限為負相關的“高-低”類型區。

從研究區2009-2016年間草地資源面積比例局部Moran’sI指數散點圖中可以得出：①2009和2016年位于第I象限的“高-高”類型區、第Ⅲ象限“低-低”類型區的村級區域單元占比較重，其中2009和2016年的個數分別為370、357個，即位于這2個象限的草地資源面積比例較高(或較低)的區域呈現出較明顯的空間集聚分布特征，且位于“低-低”類型區的村級區域單元個數高于“高-高”類型區的個數；2009和2016年位于第Ⅱ象限的“低-高”類型區、第Ⅳ象限“高-低”類型區的村級區域單元占比較輕，其中2009和2016年的個數分別為146、159個，且位于“低-高”類型區的村級區域單元個數高于“高-低”類型區的個數。②2016與2009年相比，位于各個類型區的村級區域單元發生了不同程度的變化，其中“高-高”、“低-低”類型區的個數有所減少，說明呈現空間集聚分布的空間單元在減少，呈現空間分散特征的空間單元在增加。2016與2009年相比位于第I象限的“高-高”類型區單元減少了7個，但是從空間分布情況來看，除大海鄉附近空間集聚分布情況未發生較大變化外，其他區域的空間分布發生了較大的變化，由研究區域的東北部轉變為了中東部和南部區域；2016與2009年相比位于第Ⅱ象限的“低-高”類型區單元增加了8個，從空間分布情況來看除增加了研究區的南部區域外，其余區域未發生較大的變化；2016與2009年相比位于第Ⅲ象限“低-低”類型區單元減少了6個，空間分布未發生較大的變化；2016與2009年相比位于第Ⅳ象限“高-低”類型區增加了6個，空間分布增加了北部區域，其他區域未發生較大變化。③從整體來看，研究區域2016與2009年相比，盡管位于第I和Ⅲ象限的村級區域單元個數在減少，整體上呈現集聚分布的趨勢在減弱，但總體變化不太大。

為了更好、更加直觀的表示研究區域各村級區域單元草地資源面積比例自相關水平差異程度以及變化趨勢等，在a=0.05的顯著性水平下，使用OpenGeoda軟件，并借助于ArcGIS10.0軟件獲得研究區各村級區域單元草地資源面積比例LISA集聚圖，如圖5所示。以期更好的表示空間單元與相鄰單元同質和異質的局部特征變化及其地理分布情況。

從LISA集聚圖中可以看出“高-高”、“低-低”、“高-低”、“低-高”類型區所表現的特征主要有如下。

(1)“高-高”類型區：該類型區是指其自身屬于草地資源面積比例高的村級區域單元，周圍相鄰單元也是草地資源面積比例高的村級區域單元，屬于局部均質分布的高值集聚區。2009和2016年“高-高”類型區的村級區域單元數量分別為52、40個。2009年主要分布于大海鄉，娜姑鎮的東南部，金鐘鎮的西部，待補鎮的中東部、中西部，駕車鄉的西北部，上村鄉的中部，迤車鎮中東部，馬路鄉的西南部，火紅鄉東部，礦山鎮北部等區域；2016年主要分布于大海鄉、娜姑鎮東南部、待補鎮西南部、駕車鄉西北部、大井鎮的中西部和西南部、者海鎮南部、雨碌鄉中部、田壩鄉東南部等區域。2016與2009年相比，大海鄉周邊區域仍保持“高-高”集聚區不變，主要是由于大海草山旅游景區一定程度上促進了草地資源保護，該區域地勢平緩，草地集中連片，且草地資源十分豐富。而其他區域呈現出由研究區東北部(迤車鎮、礦山鎮等)向中東部(大井鎮、雨碌鄉等)以及南部(田壩鄉)偏移的趨勢。另外，受到退草還林以及退草還耕等因素的影響，也使得2009-2016年期間出現“高-高”集聚區的變化中既有增加也有減少(縮減)情況的出現。

(2)“低-低”類型區：該類型區是指其自身屬于草地資源面積比例低的村級區域單元，周圍相鄰單元也是草地資源面積比例低的村級區域單元，屬于局部性空間差異較小的低值集聚區。2009和2016年“低-低”類型區的村級區域單元數量分別為44、42個。2009和2016年“低-低”類型區主要分布于者海鎮、新街鄉和金鐘鎮3個鄉(鎮)的交界處、樂業鎮和礦山鎮的交界處以及大橋鄉的西北部等區域。2016與2009年相比，雖然“低-低”類型區總體位置分布差異不大，但較為明顯的出現了縮減和擴大，主要表現為者海鎮、新街鄉和金鐘鎮3個鄉(鎮)交界處出現了縮減情況，新街鄉不再有“低-低”類型區，者海鎮也出現了一定的縮減；擴大情況出現在以樂業鎮東側為邊界，向樂業鎮東北部以及礦山鎮的西部進行擴大，其他區域變化不大。研究區“低-低”類型區出現縮減和擴大的原因主要是：影像判讀時，當種植的樹木較小時，影像判讀為草地；隨著樹木的生長，茂密的林地遮擋了草地，便判讀為了林地，這樣便出現了部分區域草地資源面積的增加和減少，從而影響了研究區“低-低”類型區的縮減和擴大情況。

(3)“高-低”類型區：該類型區是指其自身屬于草地資源面積比例高的村級區域單元，周圍相鄰單元是草地資源面積比例低的村級區域單元，呈現空間負相關性，空間異質性較強。2009和2016年“高-低”類型區的村級區域單元數量分別為2、0個。該類型區分布較少，2009年主要分布在大井鎮的井田村和金鐘鎮的仙龍社區，2016年研究區沒有該類型區的分布。該區域出現縮減的原因主要是出現了村莊及建制鎮的占用情況導致該區域草地資源面積減少，而鄰近單元草地資源面積出現了增加，使得彼此草地資源面積比例對比不再明顯。

(4)“低-高”類型區：該類型區是指其自身屬于草地資源面積比例低的村級區域單元，周圍相鄰單元是草地資源面積比例高的村級區域單元，呈現空間負相關性，空間異質性較強。2009和2016年“低-高”類型區的村級區域單元數量分別為10、16個。“低-高”類型區在空間上分布較為分散，2016與2009年相比，“低-高”類型區出現了不同程度的縮減和擴大變化。

3 討 論

許多學者以縣域為研究尺度，從不同草地資源類型的植被特征來研究草地資源空間分布特征[6]、草地資源退化時空分布特征[21]等，對空間分布關聯性研究較少。但近年來隨著GIS技術的快速發展[22]，運用GIS和相關地統計模型對草地資源空間分布格局進行研究成為一種可能，從不同尺度[縣、鄉(鎮)、村級區域單元]、不同角度(宏觀和微觀)對研究區草地資源動態變化進行定量分析，突出了草地資源空間分布的動態變化及其空間分布的關聯性，不僅對開展草地資源空間分布及變化特征研究具有重要的借鑒意義，而且對傳統研究中缺乏草地資源動態變化及其空間分布關聯性進行了補充。

本文從數量和空間上對草地資源空間分布格局及變化特征進行研究，但草地資源動態變化是一個長期的過程，空間分布特征及其關聯性往往受研究區溫度、降雨量、地形條件、耕作半徑等因素影響。為了更加精細化、定量化、動態化的表達草地資源的動態變化特征，下一步將基于遙感影像建立實時的草地資源動態變化監測體系[23]，進一步探討影響草地資源退化/恢復的驅動力因子，不斷促進草地資源的可持續利用。

圖5 各村級區域單元草地資源面積比例值LISA集聚圖Fig.5 LISA plots for the proportion of grassland resources area at village level

4 結 論

從縣級、鄉(鎮)級、村級區域單元3種不同的尺度，借助于地統計與空間自相關模型對研究區2009-2016年草地資源在宏觀、中觀和微觀不同層次的空間集聚和差異性特征進行定量分析，揭示了研究區在2009-2016年期間草地資源數量和空間上的動態變化。結果表明草地資源動態變化度較大的為田壩鄉，草地資源相對較為穩定的區域為大海鄉、娜姑鎮及待補鎮，與當地通過草原生態獎補機制項目的貫徹落實加快該區域草原生態恢復、防止草原退化以及與多部門聯合執法杜絕破壞草原現象發生等舉措具有密切關系。此外，研究結果可作為當地草地資源保護區劃分以及為今后生態建設和可持續發展等提供重要依據。