瀾滄江下游景觀破碎化時空動態及成因分析

劉世梁，安南南，侯笑云，董世魁，趙爽，許經緯

北京師范大學環境學院水環境模擬國家重點實驗室，北京 100875

運用遙感影像獲取地表植被覆蓋現狀及其變化信息，對于揭示地表空間分異規律、探討景觀格局動態具有重要的現實意義。目前，在景觀格局研究中應用最為廣泛的兩種方法：一是景觀格局指數方法，二是空間自相關分析法（張松林等，2007a，2007b；馮異星等，2010；吳莉等，2013；周天墨等，2013）。景觀格局主要是指構成景觀生態系統或土地利用/覆被類型的形狀、比例和空間配置（胡巍巍等，2008），通常利用景觀格局指數進行描述。目前，通過計算多時相不同景觀類型的景觀指數，分析時間序列景觀格局動態的研究較多（王小東，2011；崔世杰，2013）。而利用移動窗口法可以將景觀指數進行空間化展示，使景觀的空間格局更加清楚，因此成為研究景觀空間連續分布的重要方法。移動窗口法是通過設置窗口邊長的大小，來計算窗口內所選的景觀指數。目前移動窗口法已廣泛應用在流域、城市、干旱河谷和梯級水電站等空間梯度格局研究中（張琳琳等，2010；劉琦等，2012；陶宇等，2013；張玲玲等，2014）。

在區域尺度上，通過遙感影像反演植被指數來表征地表覆蓋信息的研究也越來越廣泛，歸一化差植被指數（Normalized Different Vegetation Index，NDVI）（崔林麗等，2010；宋富強等，2011；王曉利等，2013）是較為常用的指數。通過運用空間自相關法分析NDVI的空間聚集性及關聯性，來理解景觀異質性的研究較為廣泛。空間自相關法用于檢驗具有空間位置的某一要素的觀測值是否顯著地與其相鄰空間點上的觀測值之間表現出相對位置的依賴關系（邱炳文等，2007；張松林等，2007a，2007b）。張雪艷等（2009）、馬燕飛等（2010）、王曉利等（2013）采用空間自相關指數來分析區域NDVI的空間分布格局及聚集特征。

但是，運用兩種方法對比分析景觀異質性的研究相對較少。Chao et al.（2014）利用移動窗口法和空間自相關法分析城市景觀破碎化程度，認為兩種方法在反映景觀格局和空間異質性方面有著良好的可比性，但未利用NDVI進行指示性分析，同時基于多源數據的空間分析仍需要加強。本研究基于遙感數據和景觀指數，對比分析兩種方法在提供景觀結構特征信息方面的可靠性，以此來分析景觀的空間異質性。

另外，NDVI的變化與土地利用/土地覆蓋之間關系及其驅動因子分析也受到廣泛關注（毛德華等，2011；梁堯欽等，2012）。大量研究表明，植被分布隨地形梯度的變化可以揭示地形因子對景觀演變的影響程度（張秋玲等，2007；陳利頂等，2008）。因此，本研究通過綜合選擇多個地形因子，分析NDVI變化與外部環境的相互關系，為了解景觀格局的空間分布規律奠定基礎。

2015年，景洪市林業用地面積占全市土地總面積的 68.3%（李江龍，2015），隨著工業化、城鎮化等人類活動的加劇，該市林地生態系統監控與管理問題更加突出。本研究結合多時段遙感影像資料，通過提取連續型數據（NDVI）、離散型數據（景觀類型）和地形因子等信息，對比移動窗口法和空間自相關法定性分析區域景觀分布格局和 NDVI的空間關聯特征，了解植被景觀異質性，探尋可穩定反映景觀異質性的方法，以期從案例研究角度獲得林業資源管理的科學依據。

1 數據來源及研究方法

1.1 研究區概況

景洪市位于云南省西南部，地處橫斷山系縱谷區南端，瀾滄江大斷裂帶兩側，最高海拔約為2400 m，瀾滄江和流沙河為流經該區的主要河流。該地區屬于熱帶、亞熱帶季風氣候（劉世梁等，2006），長夏無冬，年均氣溫為 18.6～21.9 ℃，年均降水量為1200～1700 mm，季節分配極不均勻，一年之內有明顯的旱季和雨季。植被類型主要包括熱帶雨林、熱帶季雨林、南亞熱帶常綠闊葉林、南亞熱帶針闊葉混交林、竹木混交林、灌木林和草叢等。該地區是中國進入東南亞各國的主要通道和樞紐，交通建設發展迅速，人口密度增長快速。本研究區域為60 km×60 km的正方形，主要包括景洪市的中部和勐海縣的東部地區，地跨 21°40′～22°20′N，100°30′～101°10′E。以重點旅游干線和重點集鎮為骨架（耿紅生，2011）。2000—2010年，研究區內林地資源被大量占用，建設用地快速擴張，景觀格局發生了明顯變化。

1.2 數據源

采用多時相Landsat遙感數據分類，獲得所需要的景觀類型信息，根據2000年、2005年和2010年3個時段的TM影像（影像分辨率為30 m），進行人工目視判讀與監督分類，并結合2010年實地調研驗證，獲取3個時期研究區景觀類型圖，其影像的分類精度達 85%。本文參照 IGBP土地利用/土地覆被分類標準（張景華等，2011），將該區景觀類型劃分為常綠針葉林、常綠闊葉林、人工林、城市綠地、草地、水域、耕地、建筑用地和未利用地（圖1）。

圖1 研究區2000年、2005年和2010年景觀類型圖Fig. 1 The landscape type maps of study area in 2000, 2005 and 2010

利用2000—2010年的SPOT VEGETATION逐旬NDVI數據（空間分辨率為1 km），綜合考慮研究區域數據源及研究數據的有效性，在研究區域進行重采樣，生成間隔為1 km的點圖層，疊加3期NDVI柵格圖層提取每個點上的NDVI值。同時，基于DEM圖層，提取海拔、坡度、坡向、地形起伏度及地表粗糙度等5種地形因子。

1.3 研究方法

1.3.1 移動窗口法

利用 Fragstats 4.2選取反映景觀水平的景觀指數，選擇邊長為300 m矩形窗口進行計算，通過對窗口內選中的景觀特征進行計算，輸出數據為對應所選景觀指數的柵格圖。由于景觀指數類型繁雜，大部分指數在指示景觀格局特征時，存在靈敏度不高、指示含義重疊等問題（何鵬等，2009）。根據研究的實際需要，選取斑塊密度（PD）、邊緣密度（ED）、最大斑塊指數（LPI）、景觀形狀指數（LSI）、香農多樣性指數（SHDI）、蔓延度指數（CONTAG），擬從破碎度、形狀和多樣性角度探討研究區景觀格局的空間變化特征（何鵬等，

2009；張立平等，2014）。

1.3.2 空間自相關法

空間自相關法通常采用全局和局部指標進行度量。由于全局指標在一定程度上會掩蓋局部狀態的不穩定性，本文采用局部指標來分析研究區域內NDVI在不同空間單元之間的自相關程度。用以表示空間自相關的指數和方法很多，其中 Moran’s I指數是在分析自相關性中應用最廣泛的空間自相關判斷指標。局部Moran’s I指數計算公式如下：

式中，xi和xj分別是變量x在配對空間單元i和j上的取值；x是變量x的平均值；n是空間單元總數；wij是相鄰權重。本文空間權重矩陣采用的是K-nearest-neighbors，鄰域數目為4。局部Moran’s I指數的取值范圍為-1～1，正值表示該空間單元與鄰近單元的屬性值相似（“高—高”或“低—低”），負值表示該空間單元與鄰近單元的屬性值不相似（“高—低”或“低—高”）（李慧等，2011），該值為零時，該空間單元呈獨立隨機分布。對于空間是否存在自相關性，常采用z檢驗（Southworth et al.，2004）進行標準化。檢驗公式如下：

式中，E(I)為期望值；Var(I)為方差。

本文采用多時相SPOT-VEGETATION NDVI數據，通過對多個時段的NDVI與地形因子的空間自相關性進行比較，進一步解釋NDVI變化與地形因子之間的生態意義，這使得NDVI空間自相關研究更有價值。局部 Moran’s I的計算通常在 Geoda 0.9.5-i軟件中進行。

2 結果與分析

2.1 移動窗口法分析景觀格局的空間變化

2000—2010年，各景觀指數在變化不明顯（圖2）。其中，CONTAG和LPI變化明顯的區域主要分布在研究區域的中南部，該地區分布著廣泛的常綠闊葉林，在2000—2010年間，CONTAG和LPI值不斷增大，說明該地區的各類型的景觀斑塊離散程度較低，連接度較好。

以 2010年 CONTAG指數為例，研究區中心CONTAG值較低，LPI值也較低，ED和LSI值都較高，說明研究區中心的景觀破碎化程度較高，異質性高。主要是因為景洪市和景洪水電站距離研究區中心較近，工業化、城鎮化速度加快，人類干擾活動較多，土地利用隨之變化，斑塊類型更為豐富。在研究區中心的四周和北側中部區域，CONTAG值較高的區域，LPI值也較高，ED、LSI、SHDI及PD值則較低，說明該地區的斑塊較大且完整，連通性較好，景觀邊緣較規則且景觀多樣性較高。其中，研究區中心的四周出現這種現象的原因主要是該地區主要的景觀類型是人工林，在規范的管理模式下形成了較規則的斑塊。而在研究區域北側中部出現該現象主要是由于該地區景觀類型大多以常綠闊葉林、常綠針葉林為主，這些林地主要是天然林，林地形成歷史較長久，加上受到的人類活動干擾較少，因此景觀破碎化程度較低。在研究區的西側中部地區，CONTAG值較低，LPI值也較低，ED和LSI值都較高，說明該地區的景觀斑塊分散，連通性較差，其主要原因是該地區景觀類型較為多樣，主要包括灌叢、常綠闊葉林和常綠針葉林等，斑塊類型較為豐富。

2.2 研究區不同景觀類型的NDVI年際變化

利用 2000年、2005年和 2010年的年 NDVI最大值的平均值繪制時間變化柱狀圖（圖3）。整體上，研究區不同景觀類型歷年NDVI值均在0.6以上。2000—2010年間，各景觀類型的植被覆蓋趨于上升態勢。其中，未利用地的植被覆蓋呈明顯上升趨勢。

2.3 研究區時間序列NDVI整體空間自相關性

圖2 基于移動窗口法計算的2000年、2005年和2010年景觀指數柵格圖Fig. 2 Raster grid mapsof landscape indices based on moving window method in 2000, 2005 and 2010

分析研究區3期年最大NDVI的空間自相關指數（圖 4）。2000年、2005年和 2010年的年最大NDVI的局部 Moran’s I指數分別為 0.781、0.791和0.856，說明NDVI均呈現出顯著的空間正自相關，即該區域植被覆蓋都呈現集聚狀態，在空間上具有較強的關聯性。研究期間，NDVI的局部Moran’s I指數呈逐年升高趨勢，說明該區域NDVI的空間集聚格局隨著時間推進變得更加顯著。

進一步分析NDVI的局部空間自相關特征，以更好地探索NDVI的局部空間聚集模式和規律（圖5）。其中，以2010年為例，結合景觀類型圖，NDVI呈“高－高”自相關的地區大多為常綠闊葉林分布區，這些地區植被覆蓋總體較好，NDVI高值區集聚模式較明顯；NDVI呈“低—低”自相關的地區大多為建筑用地、耕地分布廣泛的地區，這些地區是人類活動聚集的地區，景洪市及景洪水電站均位于建筑用地區域，景洪市快速擴張及景洪水電站對周邊植被的破壞，都導致這些地區生態脆弱而敏感，植被覆蓋較少，而耕地呈現出“低—低”自相關模式的原因主要為該地區人類活動劇烈，農作物多為“一年兩熟”型，導致該地區的植被覆蓋不穩定，空間集聚性較差。NDVI呈“低—高”和“高—低”聚集模式的區域在研究區內分布較少。

圖3 2000—2010年間不同景觀類型NDVI變化Fig. 3 NDVI changes of different landscape types from 2000 to 2010

圖4 NDVI的局部空間自相關Moran’s I散點圖Fig. 4 Moran’s I scatterplot of local spatial autocorrelation of NDVI

圖5 不同時期NDVI空間關聯局部指標（LISA）分布圖Fig. 5 LISA distribution map of NDVI in different years

2.4 移動窗口法和空間自相關對比

分別統計各景觀類型的CONTAG和LSI值，以進一步對比分析局部Moran’s I指數和景觀指數在反映景觀破碎化時的可靠性。景觀破碎化程度較高時，表現為低CONTAG值，高LSI值。由于草地在2005年后幾乎不存在，本文忽略不計。以2010年為例（表 1，表 2），常綠闊葉林的局部Moran’s I指數平均值為0.53，表明常綠闊葉林空間分布較集聚，空間關聯性較好。同時，常綠闊葉林的蔓延度指數CONTAG為51.96，景觀形狀指數 LSI為1.35，這也表明常綠闊葉林的斑塊較少且完整，破碎化不明顯。裸地和建筑用地的局部 Moran’s I指數平均值分別為-3.44和-3.54，說明未利用地斑塊分布較為分散，空間聚集性差。同時，未利用地的蔓延度指數 CONTAG分別為24.15和39.24，景觀形狀指數LSI分別為1.62和1.41，在一定程度上表明未利用地的景觀破碎化較為嚴重。

表1 2010年不同景觀類型的局部Moran’s I指數描述性統計Table 1 Descriptive statistics of local Moran’s I index by landscape types in 2010

2.5 景觀空間分布與地形因子之間的關系

地形在一定程度上決定著景觀類型的空間分布格局，不同景觀類型的空間聚集和異常現象都與地形密切相關。地形因子與土地利用聚集或異常特征的出現關系密切，進而影響景觀類型空間自相關格局（谷建立等，2012）。提取了海拔、坡度、坡向、地形起伏度和地面粗糙度等地形因子（表3），分析了NDVI與地形因子之間的關系。海拔、坡度、坡向、地形起伏度和地面粗糙度與三期景觀類型的 NDVI的 Moran’s I指數分別為0.262、0.271、-0.018、0.302和 0.186。該研究區的 NDVI與海拔、坡度、地形起伏度和地面粗糙度的Moran’s I指數達到了0.18以上。此外，NDVI與這4種地形因子呈正相關關系。由此可以說明，本研究區植被覆蓋度一定程度上由海拔、坡度、形起伏度和地面粗糙度決定。

表3 不同時期地形因子與NDVI的Moran’s I統計值Table 3 Statistics of Moran’s I valuesbetween NDVI and topographical factors in different years

3 討論

在景觀類型和NDVI數據基礎上，結合遙感和GIS等技術，利用移動窗口法和空間自相關法，研究時間序列景觀格局的變化規律，并比較了兩種方法在指示景觀破碎化方面的可比性，為進一步探索指示景觀破碎化的合適方法提供了依據。同時，通過對景觀格局與地形因子的相關性進行分析進一步解釋了環境影響景觀格局的生態意義，使得景觀異質性研究更有價值。

分析多個景觀指數柵格圖，發現常綠闊葉林和未利用地等景觀類型，其局部自相關指數高的區域，CONTAG值也高，而LSI值較低；反之亦然。這種現象表明，移動窗口法計算的景觀指數和空間自相關指數在指示景觀破碎化時具有可比性。2000—2010年間，研究區域中心偏東南方向的人工林聚集區域CONTAG指數有逐年變大的趨勢，而LSI、ED等指數則有逐年變小的趨勢，這表明，在這期間該區域的植被景觀受到良好的保護，這和當地政府采取的“退耕還林”的政策密切相關，退耕還林后該地區植被類型的組成和結構發生了變化。

表2 2010年不同景觀類型的景觀指數描述統計Table 2 Descriptive statistics of landscape indices based on landscape types in 2010

4 結論

（1）研究區面積和移動窗口的大小，將會影響景觀破碎化的計算結果。在采用移動窗口法計算時分別嘗試150、200、300、500、1000和1500 m的邊長，通過對比分析發現，300 m的窗口邊長能保留梯度特征又不至于使景觀指數出現較大波動，可以通過景觀指標的變動特征真實反映空間格局的變化。

（2）統計對比景觀指數和空間自相關指數后發現，兩種方法在指示景觀破碎化上均具有良好的適用性。并且發現，區域內植被的分布一定程度上由海拔、坡度、地形起伏度和地面粗糙度所決定。但影響植被景觀變化的因素比較復雜，在今后的工作中可以結合土壤、水文、生物等其他自然地理要素和人為因素，分析植被景觀的多尺度變化特征，從而進行合理的景觀布局。

（3）景觀生態過程是動態的、發展延續的過程，將靜態的格局分析賦予動態變化的屬性，是把格局分析與生態過程研究結合起來的重要分析方法。今后，將結合時間序列的景觀數據和社會經濟數據等，建立過程影響因子動態變化圖譜，從而構建新的景觀動態變化格局，以分析其變化特征并實行空間格局的優化。