基于空間統計的珠三角地區土地利用生態風險時空特征

王 彬, 梅志雄, 魏俊超, 馬 杰

(華南師范大學 地理科學學院, 廣州 510631)

當前,我國的城市化進程加快推進,加速了對土地的索取,土地資源利用程度持續增加,對自然和生態環境造成了極大壓力,人地矛盾日漸突出,其中以土地利用結構不合理、土地生態風險加劇等問題最為明顯[1-3]。研究表明:土地利用是人類活動與自然環境交互過程中的關鍵部分,與很多生態系統問題密切相關,甚至對區域生態安全起到首要作用[4-5]。不同土地利用方式和過程對區域內生態系統產生累積性的影響,其影響直觀地反映在生態系統的結構和組成變化當中,并涉及到許多生態過程,由此產生了生態風險問題[6-7]。在此背景下,為探究區域土地利用變化對生態環境引起的不利影響,尋找降低其影響的方法和途徑,開展區域土地利用生態風險分析并提出相應的管控措施,有利于促進土地資源可持續利用和保障生態安全。

土地利用生態風險是指土地利用及其變化造成的生態系統結構破壞及功能降低的可能性[8],近年來受到諸多學者的關注,且日益注重表現其時空分異特性[9-10]。土地利用生態風險在不同尺度上表現出的空間格局差異較大,具有尺度效應[11]。但現有相關研究通常基于一種空間尺度,如以格網作為評價單元(樣方尺度)[12-13]或以行政區(省、區縣、鄉鎮)為評價單元[14],或小流域單元尺度等[15]。這些研究一定程度上推動了土地利用生態風險研究的進展,但仍存在如下缺陷:僅限于單一尺度的土地利用生態風險時空變化分析,普遍缺乏多種空間尺度下的土地利用生態風險規律的對比分析而未真正揭示其尺度效應,且主要依賴經驗及主觀因素選取研究尺度,缺乏客觀判斷依據[16]。目前土地利用生態風險評價方法主要有基于風險源匯和基于景觀格局兩類[4,8]。基于風險源匯的方法需要識別風險受體和風險源,涉及變量多且變量間關系復雜,數據獲取困難,且對生態風險的空間表達能力不足、未考慮空間異質性;基于景觀格局的土地利用生態風險評價方法是以土地利用變化為誘因,直接從景觀格局角度定量識別和評估生態本底安全風險[4,17-18]。這類方法所需數據易獲取,便于通過空間采樣將景觀空間結構轉化為空間化的生態風險變量,與GIS結合可深入表達生態風險的空間特征、空間異質性及其時序變化[4,7-8,19-20],逐漸成為近年土地利用生態風險評價的主流方法。目前國內外相關研究熱點區域主要集中在城市[18]、流域[7,21]、沿海地帶[22]以及重點自然保護區[23]等,而針對珠三角地區開展的相關研究較少。僅見葉長盛等[24]評價了1990—2006年珠三角地區土地利用生態風險,但缺乏2006年后的分析,且只基于5 km×5 km格網尺度。

為此,本文再次從土地利用視角,選取城鎮化快速發展、土地利用結構劇烈變化導致生態風險較為突出的珠三角地區為研究區,以2000年、2010年、2020年的土地利用數據為基礎構建生態風險指數,通過半變異函數工具對比分析樣方尺度、行政單元尺度和小流域尺度下土地利用生態風險的空間異質性,并借助空間自相關和GIS空間分析方法揭示其時空變化特征。研究將力圖揭示不同尺度下區域土地利用生態風險的空間結構模式、時空規律及其差異表現,探測其尺度效應,彌補現有相關研究的不足,并識別和探尋出需要關注的重點區域和潛在生態風險區,進而根據生態風險等級制定不同土地規劃管理措施,為區域土地資源合理利用和規劃管理及促進區域生態安全提供決策參考。

1 研究區及數據

珠三角(21°40′—24°30′N,111°20′—115°30′E)位于廣東省中南部,毗鄰南海,與東南亞隔海相望,總面積約5.48萬km2,包括廣州、深圳、佛山、惠州、東莞、珠海、中山、肇慶和江門9個地級市。改革開放尤其是新世紀以來,伴隨著經濟和城鎮化的快速發展,該區大量農業用地和生態用地轉為建設用地[25]、土地利用結構變化劇烈,嚴重擾動了生態環境,區域內生態風險問題備受關注[26]。因此,從土地利用視角分析新世紀以來該區土地利用變化所引起的生態風險時空特征,對生態風險防控具有重要理論及實踐指導意義。

本文所用的數據主要包括數字高程模型(DEM)和土地利用數據。從地理空間數據云平臺(http:∥www.gscloud.cn/)獲得2009年30 m分辨率的DEM數據,用于劃分分水嶺小流域單元。通過全球地理信息公共產品平臺(http:∥www.globallandcover.com/)獲取2000年、2010年、2020年30 m分辨率的地表覆蓋柵格數據,根據研究區實際和本文需要,將該地表覆蓋數據重新歸類得到100 m分辨率下耕地、林地、草地、建設用地、水體和未利用地6種地類。

2 研究方法

2.1 評價單元劃分

選取格網單元、行政區劃單元和小流域單元,分別對應樣方、行政單元和小流域3種空間尺度,來實現多尺度土地利用生態風險值的空間表達。具體劃分情況為:(1) 根據景觀生態學原理,通常將研究區內斑塊平均面積的2～5倍來確定格網單元的劃分尺度[19]。根據研究區土地景觀斑塊大小實際,選取0.5 km為步長,從0.5 km×0.5 km到5 km×5 km創建10種格網單元,分別得到216 060,53 940,23 998,13 488,8 629,6 025,4 542,3 383,3 061,2 160個格網評價單元;(2) 以研究區下轄569個鄉鎮為行政區評價單元;(3) 使用ArcGIS的水文分析模塊,從DEM提取出3 415個小分水嶺作為小流域評價單元。

2.2 土地利用生態風險指數

為綜合表達不同景觀類型在外界影響下產生的潛在危害性,參考前人研究[6,8,27],選擇土地利用類型面積比例、景觀干擾度指數和景觀脆弱度指數來構建土地利用生態風險指數,計算公式為:

(1)

式中:ERIk為單元k的土地利用生態風險指數;Ak和Aki分別為單元k的面積(km2)和單元k內地類i的面積(km2);Ei和Fi分別為i地類的景觀干擾度指數和景觀脆弱度指數。

2.2.1 景觀干擾度指數 景觀干擾度指數Ei表示景觀要素受到外部環境干擾的程度[28-29],可表達為:

Ei=aCi+bSi+cUi

(2)

式中:Ci,Si,Ui分別為i地類的景觀破碎度指數、分離度指數和優勢度指數;a,b,c為權重,參考相關研究[28]并結合專家咨詢分別賦值為0.5,0.3,0.2。

(1) 景觀破碎度指數(Ci)指某一地類在外界干擾下從均勻的整體狀態轉變為異質的破碎狀態,其值越大,單元內部越不穩定,可由i地類的斑塊數與單元面積的比值表示。

(2) 景觀分離度指數(Si)指某地類的斑塊單元在空間上的分離程度,其值越高,表示該地類的空間散布程度越大,容易受到外界干擾。其計算公式為:

(3)

(4)

(5)

式中:Si為i地類的景觀分離度指數;Di為i地類的距離指標;Pi為i地類的面積比值;A和Ai分別為單元面積(km2)和單元內i地類的面積(km2);ni為單元內i地類的斑塊數。

(3) 景觀優勢度指數(Ui)指某一地類對景觀整體的影響程度,其值越大,表明該地類對景觀格局的影響程度越大,可表達為:

(6)

式中:Ui為i地類的景觀優勢度;n為單元內總斑塊數;ni,Pi含義同上。

2.2.2 景觀脆弱度指數 景觀脆弱度指數Fi表示景觀要素結構受到外部因素擾動時的易損性,一般通過專家咨詢法并結合實際對不同地類賦值并歸一化獲得。參考相關研究[20,28]和專家意見,并進行歸一化得到6種地類的脆弱度指數如下:耕地0.190 5,林地0.095 2,草地0.142 9,建設用地0.047 2,水體0.238 1,未利用地0.285 7。

2.3 半變異函數

半變異函數是地統計學中研究區域化變量空間異質性的重要方法,能夠表述變量的結構性和隨機性變化特征[30]。其計算公式為:

(7)

式中:γ(h)為半變異函數;h為樣本點間隔距離即步長;N(h)為步長為h時的樣本點對數量;z(xi)和z(xi+h)分別為位置xi和(xi+h)處的土地利用生態風險值。

半變異函數有3個重要擬合參數:(1) 塊金值C0,表示區域化變量的不確定性;(2) 基臺值C+C0(包括塊金值C0和結構方差C),為半變異函數保持平穩時對應的數值,描述變量整體變異情況;(3) 變程A,表示變量保持空間自相關性的距離。常用的半變異函數模型有指數模型、球狀模型和高斯模型等,可通過模型擬合后的決定系數R2和殘差RSS大小對比來確定最優擬合模型[30]。

2.4 空間自相關方法

空間自相關用來探測空間單元與鄰域單元在相同屬性值上的關聯程度,并揭示其空間異質和空間集聚,以檢驗區域化變量的空間結構關系,包括全局和局部空間自相關,常用Moran′sI指數表示[31]。

全局Moran′sI反映區域單元觀測值在整個研究區的空間自相關特征[32],其表達式為:

(8)

式中:Xi為i單元的觀測值;Wij為行標準化的空間權重矩陣;I值范圍落在[-1,1],在定義的顯著性水平下,當I顯著為正時,單元觀測值呈空間集聚分布;當I顯著為負時,單元觀測值呈空間散布狀態;當I值為0時,表明單元觀測值為獨立隨機分布。

全局空間自相關雖可探測區域中具有同一屬性單元的總體集聚程度,但無法表示局域單元的空間集聚特征。而局部空間自相關不僅可探測空間異質,也可準確反映出集聚區的具體位置和范圍[31]。局部Moran′sI表達式如下:

(9)

式中:S為各單元觀測值的標準差。

局部Moran′sI計算結果可劃分為高—高集聚、低—低集聚、高—低集聚和低—高集聚4種類型。高—高(低—低)型即某單元與鄰域單元間表現為正空間自相關,高(低)值在局域空間呈現集聚現象;高—低(低—高)型即某單元與鄰域單元間表現為負空間自相關,高(低)值單元被鄰域的低(高)值包圍,是空間異質區域[31]。

3 結果與分析

3.1 多尺度比較及適宜評價尺度的確定

根據公式(1)—(6),借助Fragstats軟件計算得到選擇的不同空間尺度下2000年、2010年、2020年各評價單元的土地利用生態風險值,將其作為樣本點數據進行后續研究。以2000年土地利用生態風險指數計算結果為例,使用GS+軟件計算得到其半變異函數擬合結果(表1)并對比不同尺度的參數變化特征,以選取最適宜研究尺度。

表1 2000年不同尺度下土地利用生態風險指數半變異函數模型擬合結果Table 1 Land use ecological risk index semi-variance function model fitting results at different scales in 2000

由表1可知,0.5 km樣方尺度下,土地利用生態風險指數的C0和C+C0均最小,說明該尺度下土地利用生態風險的隨機性較小;隨著樣方尺度評價單元從0.5 km到5 km變化,C0和C+C0逐漸增加,說明土地利用生態風險的隨機性逐漸增大,人類活動或自然災害等外部干擾因素即偶然因素造成的空間異質性逐漸增強。行政單元尺度下的C0和C+C0均最大,說明該尺度下土地利用生態風險的隨機性特征比其他尺度都明顯。0.5～3 km樣方尺度下,其結構方差/基臺值即C/(C+C0)呈小幅波動,但均維持在0.82以上。根據地統計學原理,這說明結構性因素引起的空間異質性在樣本間總變異中起主導作用,也反映了土地利用生態風險的空間自相關性較強[16]。而3 km～5 km樣方尺度下,C/(C+C0)總體呈快速下降態勢,說明結構性因素引起的空間異質性逐漸減弱,隨機因素的影響程度加深,土地利用生態風險的隨機性變異增強。小流域尺度的C/(C+C0)達到0.839,行政單元尺度的C/(C+C0)也較高。從不同尺度擬合結果的決定系數R2來看,0.5～2.5 km樣方尺度下的R2總體上比3.5～5 km樣方尺度下的R2小,而3 km樣方尺度下的R2為0.829,為所有尺度中的最大值;小流域尺度下的R2為0.617,擬合效果尚可;行政單元尺度下的R2為0.453,擬合效果較差。因此,3 km樣方尺度是研究區2000年土地利用生態風險分析的最佳尺度。

以同樣的方法分別對2010年、2020年的土地利用生態風險進行適宜評價尺度分析,均得到3 km樣方尺度為最佳評價尺度。但考慮到小流域尺度的半方差函數擬合效果也較好以及為了保持研究區地理表面特征,最終選擇3 km樣方尺度和小流域尺度作為本文的分析尺度。

3.2 適宜尺度下土地利用生態風險的空間變異分析

由表2可知,3 km樣方尺度和小流域尺度下2000年、2010年、2020年半變異函數最優擬合模型均為指數模型。這兩個尺度下各年擬合結果的決定系數R2均較高,且殘差RSS均較小,表明各年擬合效果均較好。觀察各年最佳擬合參數中的塊金值C0,不論是3 km樣方尺度還是小流域尺度,均為2010年C0值最小、2020年C0值最大,表明2020年土地利用生態風險的空間異質性受到隨機因素的影響程度較強。兩個尺度下基臺值在2000—2020年均逐漸降低,表明土地利用生態風險的空間變異程度有所減弱,C+C0變化不大說明結構趨于穩定。兩個尺度下結構方差和基臺值的比值C/(C+C0)均為2010年最大,且2000年和2010年的C/(C+C0)均達80%以上,表明研究區土地利用生態風險的空間異質性受隨機部分因素影響較小,以結構性因素影響為主,也反映該區域土地利用生態風險的空間自相關性較強;2010—2020年,兩個尺度的C/(C+C0)均明顯下降,說明隨機因素對土地利用生態風險的影響逐漸增強。兩個尺度下2000年變程A均最大,表明該年研究區土地利用生態風險的空間相關性的作用范圍較大。

表2 2000-2020年土地利用生態風險指數半變異函數模型擬合結果Table 2 Land use ecological risk index semi-variance function model fitting results for 2000-2020

3.3 土地利用生態風險空間自相關分析

3.3.1 全局空間自相關分析 根據公式(8)計算各年兩個尺度下土地利用生態風險的全局Moran′sI值和相關指標(表3)。由表3可知,兩個尺度下各年的全局Moran′sI值均顯著為正,說明2000—2020年珠三角土地利用生態風險整體呈現正空間集聚模式,即高—高值集聚、低—低值集聚。3 km樣方尺度下全局空間正相關性整體比小流域尺度下更強;相較于2000年、2010年兩個尺度下的全局Moran′sI值小幅上升,說明2000—2010年兩個尺度下的這種正空間集聚程度均略微增強。但2020年兩個尺度下的全局Moran′sI值又略有下降,說明2010—2020年兩個尺度下的這種正空間集聚程度趨于減弱。

表3 3 km樣方尺度和小流域尺度下土地利用生態風險的全局Moran′s I和相關指標Table 3 Global Moran′s I and related indicators of land use ecological risk at the 3 km quadratic scale and small watershed scale

3.3.2 局部空間自相關分析 利用公式(9)計算各年兩個尺度下各單元土地利用生態風險值的局部Moran′sI值,然后借助ArcGIS軟件空間化得到局部空間集聚格局演化圖(圖1)。研究期內,兩個尺度下研究區大部分評價單元土地利用生態風險的局部空間結構變化較小,高—高集聚和低—低集聚單元相對較多且表現相對突出。高—高集聚區主要集中在較高風險和高風險區域,其主要地類為耕地和水域等,但其面積隨著時間推移有減少趨勢,說明土地利用生態風險有一定程度的減弱趨向;低—低集聚區主要集中分布在研究區中部,其土地利用類型較單一、多為建設用地,同時該集聚區有向外擴散趨向;低—高集聚區零星分布于高—高集聚區周圍,主要地類為耕地和草地等;高—低集聚區零星分布于低—低集聚區周圍,主要地類為建設用地和耕地,建設用地的景觀脆弱度水平較低、內部結構穩定,對減輕土地利用生態風險值有一定作用。低—高和高—低集聚區零星分布且無多大變化,說明研究期內區域土地利用生態風險的空間異質區域并未發生劇烈的變動,區域內部土地利用導致的生態風險不均衡性雖然總體表現不強烈,但依然存在。

圖1 3 km樣方尺度和小流域尺度下土地利用生態風險的局部空間集聚格局演化Fig. 1 Evolution of localized spatial agglomeration patterns of land use ecological risks at the 3 km quadratic scale and small watershed scale

3.4 土地利用生態風險的時空變化分析

前文已證明3 km樣方尺度和小流域尺度下3個年份的土地利用生態風險指數均具有顯著的空間自相關性,因此可利用克里金法進行空間化插值。由于兩個尺度下劃分的評價單元大小不同,生態風險指數值計算結果相差較大,但兩個尺度的計算結果是相對的而不影響生態風險等級對比分析。為方便對比各年土地利用生態風險空間分布狀況,采用自然斷裂點法將兩個尺度下的插值結果均劃分為5個等級,其中3 km樣方尺度下插值結果劃分為低風險區(ERI<0.065)、較低風險區(0.065≤ERI<0.086)、中風險區(0.086≤ERI<0.112)、較高風險區(0.112≤ERI<0.150)、高風險區(0.150≤ERI<0.230),小流域尺度下的等級劃分為低風險區(ERI<0.006 2)、較低風險區(0.006 2≤ERI<0.008 3)、中風險區(0.008 3≤ERI<0.010 6)、較高風險區(0.010 6≤ERI<0.014)、高風險區(0.014≤ERI<0.024)。最終得到3個年份土地利用生態風險等級空間格局圖(圖2)。

圖2 2000-2020年3 km樣方尺度和小流域尺度下的研究區土地利用生態風險空間格局Fig. 2 Spatial patterns of land use ecological risks in the study area at the 3 km quadratic scale and small watershed scale in 2000-2020

從圖2可以看出,3 km樣方尺度和小流域尺度下的土地利用生態風險值的空間分布狀況比較相似,都隨著時間的推移,總體風險程度有所降低。研究區大多都以低風險和較低風險等級為主,主要分布在珠三角的中部、西部和北部地區,包括廣州市中心城區、東莞、深圳、肇慶和佛山等地區;低風險區的土地利用類型主要是景觀脆弱度較低的建設用地,且有向其周邊擴張的趨勢;較低風險區主要分布在研究區邊緣的山地和丘陵地區,且集中在有大片林地的區域;中風險區散布于整個研究區,其范圍從2000—2020年逐漸縮小,以耕地、林地和草地類型為主;較高風險區和高風險區分布范圍相對較少,且兩個尺度下均有收縮趨勢,主要分布在順德區、南沙區、從化區、博羅縣、中山市北部以及珠海市南部的斗門區和金灣區等地。研究區中部的較低和低風險區域范圍向外擴張,面積增加,與建設用地集中分布區域有所吻合,而建設用地從碎片化逐漸演變為有序規模模式,能夠降低損失程度,對于降低生態風險有較大作用。研究區內較高和高風險區的主要土地利用類型以耕地、草地、水域和碎片化的建設用地等為主,由于耕地和水域的景觀脆弱度較大,且耕地、草地與零碎小斑塊的建設用地交錯分布在一起會造成景觀破碎度增加,同時這些區域的人類活動較為頻繁,其原有景觀的完整性遭到人為破壞,受到的外界干擾程度較大,導致生態風險值相對較高。

由表4可知,3 km樣方尺度下2000年和2010年以較低風險區和中風險區為主,面積占比分別為70.67%,69.79%;2020年較低風險區和低風險區占主導,面積占比達到75.88%。2000—2020年,土地利用生態風險值變化總體趨于降低,表現在較低和低風險區范圍變大,其中低風險區的面積增加幅度達到了8.91%;而中等風險以上區域范圍均呈縮小態勢,中風險區面積減少10.06%,較高風險區面積減少4.51%。小流域尺度下,2000年、2010年、2020年研究區主導風險區均為較低風險和中風險區,面積占比均達67%以上,最高達到72.98%。2000—2020年,土地利用生態風險值變化同樣趨于降低,最主要表現為低風險區面積增幅達到了13.63%,中風險區面積減少8.08%,較高風險區面積減少4.97%;2000—2010年的土地利用生態風險變化幅度比2010—2020年變化幅度較小,尤其是低風險區面積增加了11.11%,整體表現出風險降低的趨勢。研究區整體的土地利用生態風險值有所降低,原因在于珠三角地區加快產業升級轉換,以節約集約用地為原則,大力開展“三舊”改造工程,以提高土地利用集約程度,同時落實耕地保護與生態建設有機結合,優化土地利用結構和布局,促進區域經濟生態協調發展。

表4 2000-2020年3 km樣方尺度和小流域尺度下的土地利用生態風險等級面積和比例Table 4 Area and proportion of land-use ecological risk classes at the 3 km quadratic scale and small watershed scale in 2000-2020

總體上,兩個尺度下研究期內土地利用生態風險均逐漸減弱,且生態風險等級在2000—2010年變化幅度較小而2010—2020年變化幅度較大,其中低風險區的面積增加幅度最大,而中風險區面積減少幅度最大。

4 討論與結論

(1) 研究期內,珠三角地區土地利用生態風險分析的最佳尺度是3 km樣方尺度,小流域尺度也較為合適。這兩個尺度下,土地利用生態風險的空間異質性均以結構性因素影響為主且存在較強的空間自相關性,而在2020年隨機因素引起的空間變異程度略有加強。

(2) 總體上研究區土地利用生態風險具正空間集聚特征,兩個尺度下2010年的土地利用生態風險全局正空間集聚程度均比其他兩年略強;而局部空間集聚格局變化不大,高—高集聚和低—低集聚單元較多且表現相對突出,分別集中在高風險區和低風險區位置;低—高集聚和高—低集聚區零星分布且變化較小。

(3) 研究期內研究區土地利用生態風險以低風險和較低風險等級為主,整體上逐漸減弱,并且2010—2020年生態風險減少幅度較大。兩個尺度下生態風險空間分布格局較為相似,低風險和較低風險區范圍均逐漸擴大,以集中分布的建設用地和林地兩種地類為主;中風險、較高風險和高風險區范圍均逐漸縮小,以耕地、草地和水域為主。由于珠三角地區持續推出了一系列優化土地利用結構政策及生態保護措施,該地生態風險得到一定的控制。

基于以上結論,本文的政策啟示表現在:對于低等級土地利用生態風險區,應持續全面推進節約集約建設用地政策,并加強實施重要生態功能保護區和生態安全屏障功能區等生態修復工作;對于中、高等級土地利用生態風險區,應適當調整土地利用管理措施,提高土地集約利用水平,以增強生態系統結構穩定性。同時應合理配置生產、生活和生態空間區域,在加強耕地保護、守好耕地紅線的前提下鞏固退耕還林等政策成果,并劃定生態保護區和生態緩沖區等措施來緩解土地利用生態風險壓力,優化國土空間布局以促進區域可持續發展。

本文利用多個景觀指數構建土地利用生態風險評價模型,借助空間統計方法和GIS空間分析技術,揭示了珠三角地區土地利用生態風險的時空變化特征,并得到需要關注的重點區域和潛在的風險區域,研究結果可為區域國土空間優化利用和生態風險管控指明方向。但是,本文仍存在一些后續需進一步完善的地方:首先,本文主要從土地利用視角探討生態風險,而區域生態風險演替過程較為復雜,其影響因素眾多,未來可考慮將風險源等相關生態因素納入研究中,使生態風險研究更豐富和全面;其次,在劃分不同土地利用生態風險等級時主要參考前人的方法,目前還沒有統一的標準來確定生態風險等級,對于不同區域的研究結果可比性不強,這需要進一步深入研究。