海南島東北部土地利用與生態系統服務價值空間自相關格局分析

雷金睿,陳宗鑄,吳庭天,李苑菱,楊 琦,陳小花

海南省林業科學研究所, 海口 571100

生態系統服務是指通過生態系統的結構、過程和功能直接或間接得到的生命支持產品和服務[1],對人類健康與生存以及區域和全球生態安全至關重要[2- 4],是生態學和地理學研究的前沿和熱點[5]。土地利用/覆被變化(Land Use/Cover Change,LUCC)作為人類活動行為與地球陸表自然生態系統最直接的表征形式[6-7],對生態系統服務價值(Ecosystem Service Values,ESV)起決定性作用,其評估結果是高效、合理配置競爭性環境資源的基礎[8],也是制定生態補償[9]、可持續發展規劃政策[10- 12]以及生態系統服務權衡[13- 15]的重要前提。

目前,國內外對于ESV定量評估的研究日趨成熟,Costanza 等[1]于1997年率先在全球尺度上對生態系統服務功能進行劃分與價值評估,成為后來研究ESV的重要基礎。謝高地等[16-17]在Costanza評價模型基礎上改進并制定了“中國陸地生態系統單位面積生態服務價值當量表”,被廣泛應用在我國的森林[18]、草原[8,16]、湖泊[19-20]、城市圈[11,21-22]等區域ESV評估中。白楊等[18]利用市場價值法、影子工程法和生產成本法等定量評價了海河流域森林生態系統服務功能的經濟價值；李晉昌等[8]以修正后的價值系數為基礎,對若爾蓋高原ESV進行了評估；姚小薇等[21]采用價值當量法核算了武漢城市圈ESV,并運用空間計量方法探討不同城鎮化水平對ESV空間分異的影響；Zhang等[23]采用雙變量Moran等指數分析了武漢市ESV和城市化之間的空間相關性,結果顯示兩者存在空間負相關關系和空間溢出效應。相關研究表明,社會經濟活動引起的LUCC是導致生態系統服務供給變化的主要因素,不同土地利用格局會產生相應的生態過程,從而對生態系統服務造成影響[5,24]。近年來隨著研究的逐步深入以及3S 技術的發展,越來越多的學者也開始注重土地利用對生態服務價值時空分布格局及其影響因素的探討[5,11,20],但依然缺乏基于空間統計與分析來探索土地利用和ESV在空間上的聚集規律、關聯模式等分布特征的定量研究[20,25]。

空間統計分析可定量揭示地理變量的空間關聯性問題,是研究區域土地利用變化與ESV空間格局特征的重要工具,也是地理學領域研究的重要方向[19,26-28]。海南島作為一個獨立的地理單元,生態系統豐富多樣,但隨著國際旅游島建設的不斷推進,土地利用變化劇烈,與生態系統之間的矛盾也逐步凸顯。位于海南島東北部的海口、文昌2市是省會經濟圈建設和國家航天發展戰略所在地,承擔著區域經濟和社會可持續發展重任。本文以海南島東北部為研究對象,運用空間自相關理論分析土地利用與ESV空間格局的分布特征,旨在解決：(1)采用空間統計方法探索研究區土地利用空間自相關關系；(2)分析海南島東北部地區ESV的空間分布變化；(3)通過雙變量Moran′sI探討土地利用與ESV之間的空間相關性,分析ESV對土地利用的空間依賴性。研究結果可為區域土地利用規劃管理、生態系統服務價值的保育與恢復提供實際指導。

1 研究區概況

研究區位于海南島東北部(110°07′—111°03′E,19°20′—20°10′N),包含海口、文昌2市全境,是“海澄文”一體化經濟圈的重要組成部分,下轄41個鎮、農場,面積4748.27 km2,約占全省土地總面積的14%(圖1)。區域北瀕瓊州海峽,東面南海,西南接澄邁、定安和瓊海3市縣；地形西南高東北低,以鍋蓋嶺、馬鞍嶺為隆起高點,大部以濱海平原地貌為主,海拔在0—388 m；年均氣溫24℃,年降水量1900 mm,平均相對濕度84%,屬于熱帶海洋性季風氣候。區域內有海南島最長的河流南渡江從海口市中部穿過而入海,水資源豐富,且有海南東寨港、清瀾港兩大紅樹林自然保護區。2016年,區域常住人口278萬,地區生產總值1332億元,分別占海南省的31%和36%,是海南省政治、經濟、交通、科技、文化中心和我國文昌航天發射場所在地。

圖1 研究區域示意圖Fig.1 Schematic map of the study area

2 數據來源與研究方法

2.1 數據來源及預處理

研究采用的土地利用/覆被數據來源于2016年Landsat 8 OLI遙感影像(中國科學院計算機網絡信息中心地理空間數據云平臺(http://www.gscloud.cn)下載)解譯數據,并結合海南省森林資源二類調查矢量數據進行修正,數據精度達到95%以上。依據全國土地利用調查分類體系和研究區土地利用整體特征,將土地利用類型分為7類,分別為林地、耕地、草地、建設用地、濕地、水體及未利用地。借助ArcGIS 10.3 Create Fishnet工具將研究區劃分為1 km×1 km正方形網格單元,共計4967個,計算每個網格單元中各土地利用類型的面積。

以研究區單元格各土地利用類型的面積百分比作為空間探索變量,采用經驗貝葉斯(Empirical Bayes)平滑方法[29-30],對面積百分比數據進行標準化處理,并對標準化后的數據進行空間自相關分析。

2.2 研究方法

2.2.1土地利用程度指數

土地利用程度在一定程度上反映著人為活動的干擾程度[25]。根據莊大方和劉紀遠[31]提出的數量化土地利用程度分析方法,將土地利用程度按照土地自然綜合體在社會因素影響下自然平衡保持狀態分為4級,并分級賦予指數,計算公式為：

(1)

式中,L為研究區域土地利用程度綜合指數；n為土地利用類型的數量,這里取7；Ai為第i類土地利用類型的面積；AT為研究區域總面積；Pi為不同類型的土地利用程度參數,其中：未利用地賦值1,林地、草地、濕地和水體賦值2,耕地賦值3,建設用地賦值4[11,31]。

2.2.2生態系統服務價值評估

生態系統服務價值包含供給服務、調節服務、支持服務和文化服務4大類型所提供的價值[2]。目前,ESV評估方法大致可分為“直接市場價值法”和“條件價值法(假想市場法)”兩種,但國內相關研究仍以直接市場價值法為主[20,32],特別適用于區域及全球大尺度ESV評估[33- 35]；而后者在具體理論和方法上尚未形成突破[32]。因此,本研究依據謝高地等[16-17]的價值當量因子換算方法,并結合研究區實際情況進行修正[26]。根據《海南統計年鑒》,以2000—2015年海口、文昌兩市的平均實際糧食產量和糧食平均價格,計算出研究區ESV的當量為1307.34元,從而確定了不同土地利用類型單位面積的生態系統服務價值系數(表1),其中建設用地的價值系數為0。

表1 研究區土地利用類型生態系統服務價值系數/(元hm-2a-1)

對每個網格分別計算生態系統服務功能價值和價值強度(按面積消除行政區邊緣不完整網格的價值高低差異),并進行空間特征分析[1,36]。

每個網格的生態系統服務功能價值：

(2)

每個網格的生態系統服務功能價值強度：

(3)

研究區生態系統服務功能總價值：

(4)

式中,Aij為第j個網格i種土地利用類型的分布面積(hm2)；Cvi為第i種地類生態價值系數(元hm-2a-1)；S為每個網格的面積；n為土地利用類型；m為網格個數。

2.2.3空間自相關分析

空間自相關分析是用于衡量空間變量的分布是否具有集聚性,包含全局空間自相關和局部空間自相關兩方面[21]。為了揭示多個變量之間的空間相關性,Anselin[37]在此基礎上提出雙變量空間自相關,揭示空間單元屬性值與鄰近空間上其他屬性值的相關性[22]。本文空間自相關分析采用GeoDa 1.6.7軟件完成。

全局空間自相關能夠反應整個研究區域內各個地域單元與鄰近地域單元之間的相似性。全局Moran′sI是被廣泛應用的全局自相關統計量,計算公式為：

(5)

局部空間自相關指標(Local Indicators of Spatial Association,LISA)常采用局部Moran′sI統計量進行度量,用以準確地把握局部空間要素的聚集性和分異特征,在z檢驗的基礎上(P<0.05)繪制LISA 分布圖[30],計算公式為：

(6)

2.2.4半變異函數分析

半變異函數分析是反映區域化變量的空間相關性,進行空間變量結構分析和最優模擬的重要工具,計算公式及詳細論述參看文獻[38]。變異函數理論模型的最優選擇用決定系數R2來決定,并綜合考慮RSS、塊金值C0和變程A0。空間變量經log變換后,其分布符合正態分布規律,滿足半變異函數分析的前提條件,借助地統計學軟件GS+ 9.0應用最小二乘法以正北方向呈0°、45°、90°和135°4個角度為典型方向,分別對實測半方差進行球面模型、指數模型、線性模型和高斯模型4種不同模型的變異函數理論模型擬合,估計不同距離的半變異值,進而采用ArcGIS對變量進行普通克里格法(Ordinary Kriging)插值分析[36,39]。

3 結果與分析

3.1 土地利用格局

3.1.1全局空間自相關

從表2檢驗結果顯示,海南島東北部各土地利用類型全局Moran′sI值均大于0,P值均小于0.001,說明研究區內各土地利用類型整體上呈顯著的正向空間自相關關系,具有非常明顯的聚集性,在空間分布上并非完全隨機(表2)。其中建設用地和林地表現最為突出,空間聚集性最強,濕地和未利用地的空間自相關性相對較弱。

表2 各土地利用類型全局空間自相關顯著性檢驗表

zscore是標準差的倍數,Pvalue表示概率,z與P相關聯

3.1.2局部空間自相關

土地利用LISA分布圖(P<0.05)直觀反映了7種土地利用類型在空間聚集與分異的位置分布特征(圖2)。其中,林地主要聚集在研究區西南部鍋蓋嶺和馬鞍嶺等較高海拔地區(HH聚集,高-高聚集),植被覆蓋程度高,面積百分比均在50%以上；而林地LL聚集(低-低聚集)主要分布在海口市建成區,由于城市建設擴張造成林地占比很低。耕地HH聚集主要出現在文昌市北部的鋪前、羅豆、錦山3鎮,在城市建成區、海岸帶等不適宜耕作或農業發展的區域表現為LL聚集。草地HH型“組團”出現在城鄉結合部或庫塘周邊,這些區域很適合熱帶雜草生長或傳播,而在西南部則因喬木林地的高度覆蓋造成純草地分布很少。建設用地HH聚集明顯分布于海口、文昌城區以及重要鄉鎮。濕地和水體則“組團”聚集分布在東寨港、清瀾港、海岸帶、南渡江及庫塘等區域。而未利用地主要是城市待開發建設區域所形成的人工堆掘地或裸露地表,HH聚集集中分布在海口城區周邊、文昌北部木蘭灣及東部月亮灣等區域。

3.1.3土地利用程度空間自相關

以土地利用程度綜合指數為觀測變量,從圖3可以看出,Moran散點主要分布在第一象限(HH)和第三象限(LL),第二象限(LH)和第四象限(HL)散點分布相對較少,全局Moran′sI指數為0.5687,這說明海南島東北部土地利用程度指數具有很強的空間正相關性。

從LISA分布圖來看(圖4),土地利用程度HH聚集主要分布在海口、文昌城區及交通便利的重要鄉鎮,而在文昌市北部的鋪前、羅豆、錦山三地因有大面積集約經營的耕地,因此土地利用程度也較高。LL聚集區主要出現在未開發沿海地帶以及東寨港、清瀾港和庫塘區域,可能與自然保護區限制、利用困難或交通不便等因素有關,造成土地利用程度較低。而HL聚集和LH聚集類型則在研究區內呈零星分布。

圖3 土地利用程度Moran散點圖Fig.3 Moran scatter-plot of land use degree

圖4 土地利用程度LISA分布圖Fig.4 LISA distribution map of land use degree

3.2 生態系統服務價值評估

3.2.1各土地利用類型生態系統服務價值分析

從表3可以看出,2016年海南島東北部生態系統服務總價值為132.09億元,占2016年研究區GDP的9.92%。其中林地價值最高,占67.09%,是研究區ESV的貢獻主體；除建設用地外,未利用地價值最低,僅占0.13%。從ESV強度來看,濕地和水體最高。每個網格單元(1 km2)的平均ESV為265.44萬元,最小值為0.02萬元,最大值達716.03萬元,各網格單元價值差異較大。

表3 海南島東北部土地利用類型生態系統服務價值統計

3.2.2生態系統服務價值空間自相關分析

為避免研究區邊緣網格單元面積大小不同造成ESV相差較大,采用ESV強度進行空間自相關分析。結果顯示,海南島東北部ESV強度的全局空間自相關Moran′sI指數為0.5041,同樣表現出很強的空間正向自相關。從LISA分布圖看(圖5),ESV強度較高的區域主要分布于東北部海岸帶、東寨港和清瀾港紅樹林、以及大型水庫等,在這些濕地和水體的主要分布地帶,表現為顯著的HH聚集,也反映出濕地和水體的單位ESV很高。LL聚集“組團”出現在海口城區和南渡江近入海口兩側,人口分布密集,建設用地面積比例很高,導致ESV很弱。其次為文昌北部的耕地集中分布區和重要鄉鎮分布點,同樣因城鎮用地和耕地集中分布也造成ESV較低。而在林地集中分布的研究區西南部區域也出現有零星的HH聚集,但因差異不大,表現并不十分顯著。

3.2.3半變異函數分析

圖5 生態系統服務價值LISA分布圖Fig.5 LISA distribution map of ecosystem services value

圖6 經log變換的VES的半變異函數曲線 Fig.6 Semi-variable function sketch of ecosystem services value after log transformation

圖7所示,海南島東北部ESV強度呈現以海口城區、文昌北部農耕區為主要低值核心,東寨港、清瀾港及海岸帶為高值分布地帶的總體分布格局,模型擬合的結果與ESV強度LISA分布格局(圖5)高度一致。高值區域明顯聚集在濕地、海岸帶等重要生態保護區,生態系統服務功能很強；而低值區則容易出現在交通干線連接起來的人類活動密集區,這些區域往往是城市建成區或耕地分布的重點區域,在帶來較好經濟與社會效益的同時,也容易導致ESV損失。由此可見,ESV強度在空間分布上與土地利用格局、區位因素和經濟發展水平有很大關聯。

圖7 海南島東北部生態系統服務價值克里格插值空間分布圖Fig.7 Distribution map of Kriging interpolation for ecosystem services value in northeast Hainan island

3.3 雙變量空間自相關關系

為進一步探索海南島東北部土地利用程度與ESV之間的空間關系,采用雙變量空間自相關檢驗結果顯示,Moran′sI指數為-0.3937,表明土地利用程度與ESV之間存在顯著的空間負相關關系,即隨著土地利用程度的增強,ESV總體上呈現出下降趨勢。從土地利用程度與ESV雙變量LISA空間分布圖來看(圖8),HL型主要分布在海口城區、文昌北部以及重要鄉鎮等,LH型分布在東寨港和清瀾港紅樹林自然保護區、海岸帶、以及大型水庫等區域。從雙變量LISA顯著性水平上看(圖9),在研究區城鎮建設區、東寨港和清瀾港紅樹林分布區、庫塘、海岸帶等土地利用強度和ESV呈HL型或LH型分布的區域表現為極顯著相關(P<0.01),而在其周邊部分區域表現為顯著相關(P<0.05)。

4 討論

4.1 土地利用空間自相關格局

LUCC作為全球及區域氣候與生態環境變化的主要根源,受到自然、社會、經濟等諸多因素影響[30,40]。Verburg等[41]研究認為,自然生態環境因子在土地利用變化的空間分布上起主導作用,而社會經濟因子在決定土地利用變化的數量特征方面為主。本文研究結果表明,在研究區西南部的自然生態環境更適合林木生長的山地地區,林地呈明顯聚集狀態,水體則聚集在地勢平緩東北部區域。受區域發展和交通等因素的影響,在文昌北部的鋪前、羅豆、錦山一帶,耕地呈“組團”聚集分布；而城鎮建設用地則明顯聚集在交通干線便利的人口密集區,反映了區域土地利用的空間自相關格局和分布特征。事實上,LUCC是“人類-環境”綜合作用的動態復雜系統[42],驅動因子和驅動機制多樣、復雜[43],也造成土地利用呈明顯的區域差異性或聚集性,且具有很強的空間相關性。城鎮用地擴張、人口聚集使得土地利用程度的高值在海口城區及周邊地區的空間自相關性不斷加強,耕地的“組團”分布也形成文昌北部區域的土地利用程度較高；同時,得益于近年來海南省海岸帶生態修復和濕地保護政策,使得東部海岸帶及重要濕地周邊的土地利用程度呈顯著LL聚集。

圖8 雙變量LISA分布圖Fig.8 Bivariate LISA distribution map

圖9 雙變量LISA顯著性水平Fig.9 Bivariate LISA significance level

4.2 生態系統服務價值空間分布特征

本文采用基于修正系數價值當量因子換算的直接市場方法進行ESV評估,研究結果表明,2016年海南島東北部生態系統服務功能總價值為132.09億元,其中林地構成ESV的主體。與其他研究相比,喻露露等[26]對2012年海口市海岸帶ESV的評估結果為39.24億元,從面積換算來看兩者的評估結果比較接近；高玲等[44]采用市場價格法等方法評估2008年海口市ESV達106.02億元,評估結果有較大偏差。這是由于不同評估方法在評估模型、參數設定等方面的不同,造成結果差異較大[8,25],但這并不影響區域內ESV空間上的比較,研究結論依舊具有可靠性。在空間分布上,ESV強度表現出明顯的自相關分布特征,高值區(HH聚集)位于濕地和水體的主要分布地帶,這可能跟生態系統內豐富生物多樣性和環境組分之間物質、能量、信息相互作用有關；低值區(LL聚集)位于海口城區、文昌北部耕地分布區以及重要鄉鎮,該區域受城市土地利用開發等人為活動影響,打亂了原本系統內平衡的生態過程,在帶來較好經濟與社會效益的同時,也容易導致ESV減弱甚至退化。

4.3 生態系統服務價值對土地利用變化的響應

地理要素的空間交互作用中普遍存在空間溢出效應,它是指由基于位置的鄰近性所產生的空間外部性,即一個單位從其鄰居帶來收益或成本[23,45]。大量研究表明,土地利用程度對ESV呈負相關[20,25],且存在明顯的空間自相關性和空間溢出效應。總體而言,土地利用程度對ESV有負面的空間溢出效應,也就是說,一個地區的土地利用程度的提高可能會導致周邊地區的ESV的退化。土地利用程度較高的如建設用地、耕地等是ESV的低值分布區,而土地利用程度較低的如濕地、水域則為ESV的高值分布區。因此,在高ESV附近的區域,應避免耕地活動和城市建設開發產生的溢出效應對自然生態系統造成的負面影響。相反,城市自然邊緣應該致力于生態友好的土地用途(例如郊野公園、綠地系統)[46]。根據陳利頂等[47]提出的“源-匯”景觀理論,本文研究認為,控制建設用地、耕地等“匯”(負服務)景觀類型,保護與恢復濕地、水域、林地等“源”(正服務)景觀類型,預防“匯”景觀類型向“源”景觀類型擴展或滲透,搭建生態系統內部良性循環,積極引導土地利用向ESV保值和增值方向發展,是優化與提高區域ESV的有效途徑。

4.4 LUCC與ESV空間關系對區域景觀規劃的啟示

LUCC與ESV是一個相輔相成的綜合系統,開發建設還是生態保護？這是區域生態規劃者和管理者必須面對的問題,如何在不破壞生態可持續性的前提下分配建設活動,應當是生態系統服務權衡的重要任務[23,48- 49]。

根據土地利用程度與ESV雙變量空間關系的四種聚類模式認為,研究區內具有低土地利用程度和高ESV的區域(LH型)是生態優勢區域,可以被認為是區域性的“生態核心區”,比如東寨港和清瀾港紅樹林自然保護區、海岸帶以及大型水庫等,應該采取最嚴格的保護措施,禁止或限制城市擴張[23,26,50]。在高土地利用程度和低ESV的區域(HL型),比如海口城區和重要鄉鎮,應當營造綠色空間或減少建筑強度來緩解人類對生態系統的壓力,引導土地利用向高ESV發展。在高土地利用程度和高ESV地區(HH型),應鼓勵生態友好的土地利用類型(例如城市森林、草地等),而限制高污染土地利用類型。而在低土地利用程度和低ESV地區(LL型),則應考慮更多的生態修復工程,也可考慮在這些地區進行一些建設項目,增強土地利用效率。此外,還應規劃和實施有效的綜合生態修復工程,根據生態系統退化程度的優先次序有計劃地實施。

5 結論與展望

(1)海南島東北部各土地利用類型的全局Moran′sI值均大于0.4,P值也均小于0.001,具有顯著的空間自相關性,其中建設用地和林地的空間聚集性最強；受人為干擾和環境作用的綜合影響,土地利用程度有顯著差異與空間聚集特征,不同土地利用類型所表現的空間聚集或異常的區域及范圍也明顯不同。

(2)研究區ESV有顯著的空間正相關性,空間集聚程度較高；高值區聚集于海岸帶、東寨港和清瀾港紅樹林、以及大型水庫等,低值區集中在海口城區、文昌北部農耕區域。

(3)研究區土地利用程度與ESV呈空間負相關關系,即土地利用程度越高,ESV越低,且具有空間溢出效應。根據空間聚類模式提出區域保護措施,為城市景觀規劃與修復提供實際指導。

本文運用空間自相關理論探討了土地利用與ESV空間格局的分布特征與自相關關系,根據研究結論建議繼續實施海岸帶生態修復與保護工程、紅樹林保護管理等政策,優化建成區土地利用結構,提升生態核心的服務價值,對實現區域可持續發展和維護生態安全具有重要的實踐意義。本研究也存在不足之處。首先,生態學和地理學的特征與現象往往伴隨著復雜的尺度效應,土地利用與ESV的空間自相關也表現出明顯的尺度依賴特征,如何更科學的確定空間自相關的尺度還有待于深入研究。其次,生態系統服務具有內部關聯性,但社會經濟因子與生態系統之間聯系不應被忽視,尤其是不同土地利用變化驅動情景下的生態系統過程與服務的關系,如何揭示其背后的生態過程和響應機理應當是今后重點研究的方向之一。

致謝：感謝海南大學宋希強、王華鋒和清華大學楊軍老師對本文寫作的幫助。