黃秋霞,廖南燕,楊元征,3,4,俎佳星,3,4,王嘉麗,3,4,蔡文華,3,4,楊 健,3,4,蘇宏新,3,4,楊 婷
(1. 南寧師范大學地理科學與規劃學院,南 寧530001;2. 廣西壯族自治區防城金花茶國家級自然保護區管理處,廣西 防城港 538021;3. 北部灣環境演變與資源利用教育部重點實驗室,南寧 530001;4. 廣西地表過程與智能模擬重點實驗室,南寧 530001)
自然保護區的建立是保護珍稀動植物及其生境的重要舉措(Maiorano et al., 2008; Lung et al.,2010; Radeloff et al., 2010),其在緩解土地利用變化、保護生態多樣性方面起積極作用(Stein et al.,2008;Xu et al.,2017)。自1956年中國建立第一個自然保護區(鼎湖山國家級自然保護區,面積為11.33 km2)以來,自然保護區數量和面積不斷增加,建設體系日益成熟(張鐿鋰等,2015),了解自然保護區是否有效實現保護目的逐漸成為熱點話題,但其建設成效評估體系仍不完善(祝萍等,2018)。目前,評價保護區建設成效的方法很多,至今尚未形成統一標準。土地利用變化能夠直觀表現自然保護區保護成效(郭子良等,2021),但簡單分析保護區成立以來的土地利用類型面積變化并不能達到有效評價的目的。從長時間序列深入了解保護區及其鄰近區土地利用類型動態演變,可以綜合評估保護區建設成效,為土地資源規劃提供決策支撐,對維持生態系統平衡和資源環境的可持續發展具有重要意義(Lambert et al.,2017)。
土地利用是鏈接人類社會經濟活動與自然生態過程的紐帶(Liu et al.,2014),其利用方式轉變對生物多樣性(García et al., 2007)、生態系統服務(田義超等,2015;劉曉娟等,2019;Chen et al.,2020; Cao et al., 2021)、資源的可持續利用(姚靜韜等,2020)以及氣候變化都有深遠影響(Bu et al.,2014;Huang et al.,2015)。近年來,國內外學者對自然保護區土地利用變化做了大量研究,主要利用地理信息系統技術,結合數理統計與景觀生態學方法分析保護區建立后土地利用和景觀格局動態變化(Zhang et al.,2017;張宏靜等,2021),探討造成該動態變化的驅動因素(Su et al., 2015;滿蘇爾·沙比提等,2016,2017),并預測未來土地利用變化趨勢及其對生態系統服務的影響(井云清等,2016;Zhang et al.,2019;Li et al.,2021)。但目前研究較多集中在保護區內部,較少針對保護區內部與外部土地利用動態差異。
廣西防城金花茶國家級自然保護區(簡稱“金花茶保護區”)位于廣西壯族自治區防城港市境內,于1986年經廣西壯族自治區人民政府批準,建立廣西防城上岳金花茶自治區級自然保護區,1994年晉升為國家級自然保護區,主要保護對象為珍稀瀕危金花茶組植物及其生境,是中國乃至世界金花茶植物的分布中心。目前對金花茶保護區的研究主要集中在植物多樣性(龐潔等,2008)、植物資源保護及開發(謝釗等,2019)、昆蟲區系分析(莫璐瑕等,2014)等方面,對其土地利用變化的研究比較少。保護區的建立使人類活動逐漸減少,對金花茶生境質量提升起一定的積極作用,但在當地經濟發展的需求下,部分人類活動對金花茶保護區土地利用格局仍然造成影響(黃曉園等,2019)。保護區及鄰近區域土地利用變化會威脅保護區生態安全與生態系統服務功能。因而,有必要探索保護區的建立對區域土地利用產生的動態影響。鑒于此,以1986—2019 年7 期Landsat5 TM/Landsat8 OLI 遙感影像為數據源,采用決策樹分類方法進行土地利用分類,探討金花茶保護區非核心區和保護區外部2 km鄰近區的土地利用類型動態演變規律。以期為保護區的建設及監管、資源開發及生態安全保護提供一定的數據支撐。
金花茶保護區位于廣西壯族自治區防城港市防城區境內(李廣超等,2021),總面積為9 098.62 hm2,其核心區面積為1 479.10 hm2,非核心區面積為7 619.52 hm2,非核心區包括保護區內實驗區和緩沖區兩部分,其中緩沖區面積為3 459.21 hm2,實驗區面積為4 160.31 hm2,跨越那梭、華石、大菉和扶隆4個鄉鎮,地理位置為21°44′-21°48′N、108°03′-108°12′E(圖1),保護區內人口約1 900人(實地走訪調查所得),糧食作物以玉米和水稻為主,經濟作物主要有紅薯、蔬菜、八角、肉桂。保護區地處十萬大山山脈,屬熱帶季風氣候,全年氣候溫和濕潤,冬短夏長(李廣超等,2021),年均溫約21.8℃,年均降雨量高達2 900 mm。植被類型以闊葉混交林、亞熱帶常綠闊葉林為主。為更好地評估保護區建設與發展對當地生態環境的保護力度,將研究區劃分為保護區內的核心區(簡稱“核心區”)、保護區非核心區(簡稱“非核心區”),保護區外部鄰近區(簡稱“鄰近區”)三部分(見圖1)。其中,鄰近區為在保護區外部建立2 km 的緩沖區,外部緩沖區面積與保護區面積大致相同。通過對比保護區外部鄰近區與非核心區土地利用變化趨勢,可為綜合評價保護區建設成效提供數據支撐。

圖1 研究區示意Fig.1 Spatial distribution of the study area
Landsat 系列數據具有時間跨度長、重訪周期短、覆蓋范圍廣等優點,被廣泛用于土地利用時空格局演變研究中。基于Google Earth Engine(GEE)云平臺①https://earthengine.google.com/獲取金花茶保護區及其鄰近區1986、1994、1999、2004、2009、2014、2019 年7 期Landsat5 TM/Landsat8 OLI 地表反射率產品,影像行列號為125/45,空間分辨率為30 m。選取當年云量最小影像,消除云雨等因素的影響。該數據產品已進行大氣校正與正射校正處理。
因核心區人為干擾較少,土地利用狀況穩定,主要對比非核心區和鄰近區的土地利用變化。采用決策樹分類方法進行土地利用遙感解譯,并基于分類結果對比非核心區和鄰近區的土地利用演變規律。為定量評價土地利用時空演變規律,選取單一土地動態度、綜合土地利用動態度和土地利用程度綜合指數表征土地利用變化動態度,利用數量轉移軌跡與質心遷移軌跡定量展示土地利用時空演變特征。
1.3.1 土地利用類型解譯 基于金花茶保護區土地利用實際情況,并結合Landsat5 TM/Landsat8 OLI遙感數據的特點,參考中國科學院土地利用/土地覆蓋分類系統,利用ENVI 5.3 軟件進行決策樹分類,將研究區土地利用劃分為林地、草地、耕地、水域、建設用地和未利用地6類(未利用地包括裸地、裸巖等),并依據野外調查數據和Google Earth高分辨率遙感影像進行精度驗證。
1.3.2 土地利用類型動態演變 不同土地利用類型的變化速度可以通過土地利用動態度定量描述(鄔亞娟等,2020),單一土地動態度K可用于分析某一用地類型在一定時間范圍內的數量變化情況(張浚茂等,2019),綜合動態度LC可以表征研究區不同時間節點全部用地類型的整體變化速率(張浚茂等,2019),計算公式為:

式中:Uai、Ubi分別表示研究初期和研究末期某單一土地利用類型i的面積;T為研究時長;LUi為研究期初期第i類土地利用類型的面積;ΔLUi-j為研究期內第i類土地利用類型轉為其余土地利用類型面積的絕對值(張浚茂等,2019)。
土地利用程度可以反映人類活動對土地利用變化的影響程度(樊嘉琦等,2018),計算公式為:

式中:La為研究區土地利用程度綜合指數,數值越大,表明該區域人類活動越強烈;Ri為地類i的土地利用程度的分級指數,參考自文獻(鄔亞娟等,2020)(表1);Ci為地類i占研究區總面積的百分比;n為土地利用類型的數量。

表1 土地利用程度分級賦值Table 1 Grades of land use degree classification
1.3.3 土地利用數量變化軌跡 為認識非核心區和鄰近區不同時段內各土地利用類型之間轉移的方向和數量,計算不同時段的土地利用變化轉移矩陣(張佰發等,2020)。轉移矩陣的公式為:

式中:Yij表示某時期內第i種土地利用類型轉化為第j種土地類型的面積百分比;n為土地利用類型的數量。
通過R語言中的NetworkD3 包,采用桑基圖(Sankey diagram)將土地利用轉移矩陣進行可視化,構建土地利用數量變化軌跡模型。桑基圖是一種顯示物質或者能量流動的可視化方法(姜婷婷等,2016),在材料分析、金融、生態學等多個領域被廣泛運用。
1.3.4 土地利用遷移軌跡 為有效反映土地利用類型的時空遷移過程,采用標準差橢圓法(The Standard Deviational Ellipse,SDE)計算各地類的質心坐標。其質心坐標計算公式為(劉濤等,2013):

式中:SDE(x,y)表示某一地類的算數平均中心坐標,即重心;(xi,yi)表示某地類第i個圖斑的空間區位坐標;n為土地利用類型的數量。
利用ArcGIS 10.5 疊加各時段不同地類的重心坐標,得到鄰近區和非核心區土地利用重心遷移軌跡,運用歐式距離方法計算得到不同土地利用類型重心遷移距離。
采用決策樹分類法進行遙感解譯得到金花茶保護區及鄰近區1986、1994、1999、2004、2009、2014和2019年7期土地利用數據(圖2),根據野外調查數據和Google Earth 高分辨率遙感影像對土地利用分類結果進行精度驗證,7 期遙感圖像的分類結果總體精度均>85%,Kappa 系數均>0.81,分類精度能滿足研究需要。

圖2 1986—2019年研究區土地利用Fig.2 Land use pattern in the study area during 1986-2019
由7個時期土地利用時空分布情況可知,保護區(核心區和非核心區)土地利用類型主要為林地;未利用地和草地呈片狀分布于中部、中偏西部海拔較高區域,呈現一定的集聚性;耕地呈點帶狀環保護區邊緣較平緩地勢分布(見圖2)。鄰近區土地利用類型主要為林地、耕地;林地成片狀集中分布,東部、西部分布較多,南部分布較為破碎;南部有一條明顯的耕地帶,呈西南—東北走向,建設用地成片狀分布,與耕地的分布具有一定的空間協同性。整體而言,金花茶保護區及其鄰近區土地利用類型主要為林地、耕地,二者面積之和占總面積的93%以上(圖3),呈現集聚性分布;林地面積分布最廣,占比為76.73%~79.46%,主要分布在研究區中部及西南部;其次是耕地,面積比例為15.60%~18.36%,呈連片狀分布,南部分布較多;水域、未利用地、建設用地、草地占總面積比例較少,面積比例均小3%,其中建設用地面積分布最少。近35 年間1986-2019 年土地利用類型發生較大變化,建設用地面積占比增長幅度最大,且呈現逐年上升趨勢,面積占比從1986年的0.19%增長至2019 年的0.93%,2004—2009 年增長了0.29%,較為迅速;林地呈現波動性上升趨勢,從1986 年的78.09%上升到了2019 年的78.82%;草地面積呈現略微增加趨勢,面積占比從1986年的1.30%增長至2019年的1.68%;未利用地面積縮減幅度最大,面積占比從1986 年的2.02%減少至2019 年的1.07%,2004-2009 年面積占比減少了1.87%;耕地和水域呈現先增加后減少趨勢,總體上面積占比分別減少了0.76%和0.13%。

圖3 金花茶保護區及鄰近區土地利用面積變化趨勢Fig.3 Land use changes at various stages in Golden Camellia Reserve and its adjacent areas
根據非核心區和鄰近區單一土地利用類型動態度可以看出(表2),非核心區建設用地單一土地利用類型動態度居于第一,其中,建設用地1986-2019年平均動態度為161.60%,總體呈現顯著擴張趨勢,1994—1999年動態度最高為56%,2009年之后動態度為負值且呈現下降趨勢;草地的動態度相對于其余用地類型一直處于較高水平,年平均動態度為1.84%,總體位居第二且呈現擴張趨勢;林地呈現微增趨勢;未利用地1986—2019 年平均動態度為-1.58%,呈現波動減少趨勢;水域和耕地的動態度波動較平穩,呈微小降低趨勢。鄰近區建設用地的動態度一直處于較大值(見表2),1986-2019年平均動態度為11.42%,呈現擴張趨勢;其余5 種土地利用類型的動態度波動較為平緩,林地呈現微小增加趨勢,草地、未利用地、水域和耕地呈微小縮小趨勢。綜上,非核心區和鄰近區建設用地的動態度最高,但非核心區的建設用地動態度更大,2個地區近35年的單一土地利用類型動態變化除草地外其他類型趨于一致,非核心區草地為正值且排名第二,鄰近區草地則為負值。

表2 鄰近區和非核心區單一土地利用動態度Table 2 The rate of land use change at various stages in adjacent areas and non-core areas%
1986—2019年,鄰近區綜合土地利用動態度曲線波動較大(圖4-a),呈現先急劇減少再上升到最高值,之后波動下降的趨勢。1994—1999年處于較低值,1999—2004 年上升到最高值0.63%,表明該時期人類活動對鄰近區土地利用變化造成較大影響;2004年以后呈波動下降趨勢,土地利用變化趨于平緩。非核心區1986—2019年綜合土地利用動態度呈現先增加后遞減的趨勢,2004—2009年處于最高值0.55%,該時期土地利用變化劇烈。整體而言,鄰近區綜合土地利用動態度比非核心區變化更劇烈,表明保護區內部非核心區土地利用格局更為穩定,受人為活動影響較小,鄰近區人地矛盾更為劇烈,2 個區域2009—2019 年綜合動態度都趨于下降,表示土地利用變化都趨于平緩。

圖4 鄰近區和非核心區1986—2019年綜合土地利用動態度(a)和土地利用程度綜合指數(b)Fig.4 Comprehensive land use dynamic degree(a)and comprehensive index of land use degree(b)in the adjacent and non-core areas from 1986 to 2019
鄰近區和非核心區土地利用程度綜合指數變化趨勢相似,均呈現先上升再下降,之后較為平緩波動趨勢(圖4-b)。鄰近區在1999 年達到最大值229.58后略微下降,2004年以后平緩上升,這與該時期土地利用方式轉換相一致,表示土地利用程度高的建設用地在擴張,未利用地等土地利用程度較低的地類面積占比減少,土地利用強度顯著提高。非核心區土地利用程度綜合指數在1999年達到最大值,隨后呈波動下降趨勢,該區域土地利用程度綜合指數始終比鄰近區低,表明非核心的土地利用受人為干擾較少,建設用地、耕地等高土地利用強度所占比例得到控制并逐漸減少。
從金花茶保護區土地利用數量轉換軌跡看(圖5) (空間上的土地利用類型轉移情況見文末附圖1),鄰近區1986—2019 年土地利用類型轉換頻繁(圖5-a)。其中,林地轉換量最為顯著,主要流向耕地、草地,1986—1994 年流向耕地面積高達1 423.33 hm2,1994—1999 年變化相對較少,其他時期轉移面積波動減緩;耕地以流出為主,主要流向為林地、建設用地、草地,1986—1994年流向林地面積最大,達到1 149.14 hm2;建設用地流入行為明顯,主要由耕地、未利用地、林地流入,由耕地流入面積逐年增加。就土地利用類型轉換的軌跡線而言,1986—1994 年不同地類之間轉化最為活躍,轉換軌跡多樣化;1994—1999年不同地類之間轉換相對較少;1999—2019年4個時段不同土地利用類型間的轉變幅度較為一致。整體上,1986—2019年鄰近區地類轉換軌跡線比較稠密,各地類相互轉換比較多樣,不同類型間變化較活躍。

附圖1 金花茶國家級自然保護區鄰近區、非核心區1986—2019年土地利用變化情況Accompanying Fig.1 Land use change in adjacent and non-core areas of the Fangcheng Camellia National Nature Reserve from 1986 to 2019
非核心區1986—2019年土地利用類型之間轉換較少(圖5-b)。林地流入行為最為明顯,主要由耕地、草地、未利用地流入,耕地轉入面積最大,1986—1994 年流入面積達最大值212.75 hm2,其余年間轉換面積相對較少;草地以流入為主,主要由林地、未利用地、耕地轉入,2004—2009年由未利用地轉入52.73 hm2;未利用地流出明顯,主要轉化為草地、林地,2004—2009 年轉換為草地104.22 hm2;耕地、水域和建設用地與其他土地利用類型之間的轉移變化相對較少。與鄰近區對比,非核心區軌跡線較為稀疏,說明該區域土地利用類型變化較少,且轉變類型更簡單。綜上,保護區的建設與管理顯著效果,當地生態環境得到有效保護。

圖5 1986—2019年土地利用數量變化軌跡Fig.5 Land use change trajectory from 1986 to 2019
土地利用重心遷移軌跡分析結果表明(圖6),鄰近區林地重心呈現三角形狀遷移特征,總體向東北方向遷移,1999—2004年遷移幅度最大;耕地重心呈現“V”字型遷移特征,總體向西南方向遷移,1999—2004年遷移幅度最大;建設用地重心向西北方向遷移,1986—1994年遷移幅度最大;草地重心向東北方向遷移,2004—2009年遷移幅度最大;未利用地重心向東北方向遷移,1986—1994年遷移幅度最大;水域重心遷移幅度較小,向東北方向遷移。耕地與建設耕地總體遷移方向都為向西遷移。

圖6 鄰近區、非核心區土地利用重心遷移軌跡Fig.6 Migration trajectory of land use gravity center in adjacent areas and non-core areas
非核心區6種土地利用類型重心整體上都是向東北方向遷移,建設用地重心遷移幅度最大,在1986—1994年遷移距離最遠;林地在1999—2004年重心遷移幅度最大;耕地重心呈現“W”字型遷移特征,2004—2009年遷移幅度最大;草地在1986—1994年重心遷移幅度最大;未利用地在2014—2019年重心遷移幅度最大;水域重心呈現“Z”字型遷移特征。
綜上,保護區周邊鄰近區土地利用變化劇烈,土地利用類型轉換頻繁,非核心區土地利用變化較平緩,土地利用類型相互轉換較少,表明保護區建立后管護政策的貫徹實施,管護力度不斷加強,減少了人為活動的干擾,保護區的生態環境質量得到改善。1986—2019年非核心區和鄰近區建設用地動態度變化最大,都為擴張趨勢,主要由耕地流入,保護區為加強開發管理對自然保護區內的公路進行擴建,這是非核心區用地動態度變化的主要原因。鄰近區土地利用變化主要受那梭農場和沖皇溝生態旅游區的發展建設影響,西南方向原有居民地向西北方向不斷擴張,說明西部人類活動逐漸增多,社會經濟活動活躍,需對保護區西部加強管理;退耕還林、沼氣池建設工程等政策的實施,使非核心區和鄰近區林地面積都有所增加,非核心區的林地因得到保護增幅更為明顯,2 個區域林地都往東北方向遷移,與森林資源消耗的減少、東北方向耕地退耕有一定關聯;兩個區域的耕地均呈減少趨勢,與建設用地、林地的動態演變具有一定協同性;非核心區未利用地主要分布在海拔較高山頂區域;隨著保護區生態環境的改善,西南方向較低海拔的未利用地轉換為林地和草地,鄰近區的未利用地主要轉換為林地與耕地。
從時間梯度看,鄰近區綜合動態度在1994—1999 年處于較低值,土地類型轉換較少,這與90年代以后外出務工人員大量增加,土地利用受人為干擾下降有一定關系,當地居民通過務工收入水平顯著提高,擴大住宅面積,土地利用程度高的建設用地在擴張,使得土地利用程度上升;1999—2004年,那灣生態園的建設、那梭農場特色產業的開發,建設用地在西部顯著擴張,人類活動劇烈,土地利用類型之間轉換劇烈,以致土地綜合動態度較高,2002年退耕還林政策的實施,使得耕地面積減少,土地利用程度下降;2003年國家重點公益林補助試行之后,當地居民建設沼氣池的積極性提高,減少了森林資源消耗,退耕還林政策的繼續實施,有效控制了耕地的面積,但建設用地不斷擴張,土地利用動態度趨于波動增長,土地利用程度平緩提升。非核心區受道路建設、退耕還林等政策的影響,土地利用變化較強烈,土地利用類型之間轉換較頻繁,1986—2009年土地利用動態度一直處于增長趨勢,2009年之后呈現下降趨勢,表明人類活動逐漸較少,土地利用格局趨于穩定。
整體而言,保護區的設立,退耕還林還草等工程的實施,管護工作力度的不斷加強,使得保護區內部土地利用受人為干擾下降,土地利用格局更為穩定,建設用地、耕地等高土地利用強度所占比例得到控制并逐漸降低。同時,保護區的建立也促進了周邊生態旅游、礦產資源、水資源的開發以及經濟作物產業的發展,促使周邊人類活動增加,土地利用變化劇烈,土地利用動態度一直處于較高水平。總體上,保護區的管理取得較為顯著的成效,內部生態環境得到改善,但內部建設用地的緩慢增加與資源的不合理生產開發,仍然對保護區的生境質量產生負面影響。保護區在未來的管理建設中,應合理控制區內人類活動,強化土地集約化發展,有效實施生態移民搬遷、異地安置等政策。
金花茶保護區位于海拔較高區域,其地形特點在一定程度上也驅動著土地利用類型及其動態演變(崔步禮等,2011)。為更精準地評估保護區的建設成效,提取1986—2019年非核心區坡度<6°的土地利用數據,并進行土地利用動態度的計算,1986-2019 年7 個年份土地利用強度分別為224.41、221.54、223.96、223.02、224.15、222.01、221.20,均小于鄰近區的226.99、229.58、229.59、226.51、228.72、227.68、228.86,證明保護區的建設對生態環境的提高有促進作用。
利用1986—2019年遙感影像對金花茶保護區及其鄰近區進行土地利用遙感解譯,并基于GIS空間分析與R語言統計分析方法,探討了研究區近35年來土地利用動態演變規律,得到的主要結論有:
1)研究區土地利用類型整體分布具有一定的集聚性,保護區土地利用類型以林地為主,外部鄰近區以林地和耕地為主,近35年間土地利用類型發生較大變化,建設用地面積占比增長幅度最大,未利用地面積縮減幅度最大,其余地類呈現波動增減趨勢。
2)自保護區建立以來,鄰近區土地利用變化強烈,非核心區土地利用格局更為穩定;非核心區和鄰近區建設用地的動態度一直處于較大值,其余5 種土地利用類型的動態度波動較為平緩;鄰近區綜合土地利用動態度比非核心區變化更劇烈;鄰近區和非核心區土地利用程度綜合指數變化趨勢相似,均呈現先上升再下降,之后較為平緩波動趨勢。
3)鄰近區近35年間土地利用類型之間轉換頻繁,相對而言,非核心區土地利用類型相互轉換較少,2 個區域均是耕地與林地之間的轉移交換最為劇烈。鄰近區轉換軌跡線比較稠密,各地類相互轉換呈多樣化,非核心區軌跡線較為稀疏。
4)土地利用重心遷移具有一定的方向性,非核心區6 種土地利用類型總體上都是向東北遷移,鄰近區建設用地與耕地重心向西遷移,其余地類重心向東北遷移,表明建設用地、耕地在已有基礎上不斷擴張。
5)保護區的管理取得較為顯著的成效,保護區內部土地利用受人為干擾下降,土地利用格局更為穩定,生態環境得到改善,但內部建設用地的緩慢增加與資源的不合理生產開發,仍然對保護區的生境質量產生負面影響。
本研究定量對比分析了金花茶保護區非核心區和外部鄰近區在保護區建設與發展歷程中土地利用動態演變,探究保護區的建立對區域土地利用動態的影響,可為保護區的建設及監管、資源開發及生態安全保護提供一定的數據支撐。