礦區土地利用生態風險的時空異質性

常小燕,李新舉,*,李西燦,郭 鵬,高 峰

1 山東農業大學信息科學與工程學院,泰安 271018 2 山東農業大學資源與環境學院,泰安 271018 3 山東省濟寧市國土資源局,濟寧 272001

生態風險是指一個種群、生態系統或整個景觀的正常功能受外界脅迫,系統內部某些要素或其本身的健康、生產力、遺傳結構、經濟價值和美學價值發生退化的可能性[1]。1992年美國環境保護署(U.S. Environmental Protection Agency)將生態風險評價(Ecological Risk Assessment,ERA)定義為“評估暴露于一種或多種壓力因子后,可能出現或正在出現的負面生態效應的可能性過程”[2]。區域生態風險評價作為生態風險評價的一個分支,強調在區域尺度上描述和評估環境污染、人為活動或自然災害對生態系統及其組分產生不利作用的可能性及其大小的過程[3]。空間異質性是指某種生態學變量在空間分布上的不均勻性及復雜程度,即系統特征在時間和空間上的復雜性和變異性[4]。

我國區域生態風險評價起步較晚,主要是對一些生態敏感脆弱區[1,5-19]及城市[20-28],以“風險源-風險受體-暴露危害分析-生態終點”這一概念模型,構建綜合生態風險評價指標體系,并借助遙感(RS)和地理信息系統(GIS)技術完成研究區生態風險評價[1,3,5-16,20]；也有一些學者通過對遙感影像解譯獲取土地利用類型數據,運用景觀格局指數構建土地利用生態風險指數,借助空間統計學、地統計學等研究方法,評價土地利用變化所帶來的區域生態風險[17-19,21-29]。但不足之處主要在于基于概念模型進行土地利用生態風險評價時,對土地利用生態風險的空間異質性缺乏深層次的分析,僅限于生態風險的狀態分析；運用地統計學進行土地利用生態風險空間異質性分析時,大多數研究者人為地根據研究區大小及經驗選取研究尺度,缺乏科學依據。

從景觀生態學角度將土地利用格局與生態風險相結合以評價區域生態環境,結合空間統計學理論進行土地利用生態風險空間關聯性及時空變異性研究,有利于揭示區域土地利用生態風險的時空分布格局及空間變異規律,分析生態風險與土地利用類型、人類活動之間的耦合關系,為土地資源的優化配置及開發整治提供理論基礎,以促進土地資源的可持續利用。微山縣作為山東省的“南大門”,境內有豐富的礦產資源,同時山東省最大的淡水湖—南四湖也坐落其中；本文從微山縣自身生態風險的特殊性入手,在現有礦區生態風險評價研究成果的基礎上,選取適宜的景觀格局指數構建景觀意義上的土地利用生態風險指數,結合空間自相關分析及半變異函數理論,選取適宜的研究尺度,進行研究區土地利用生態風險半變異函數模型擬合,在此基礎上進行空間插值分析,以探究2000—2016年微山縣適宜尺度下土地利用生態風險的時空變化情況及空間異質性,實現礦區土地利用生態風險空間異質性的定量化和可視化表達,以期為后續的礦區生態恢復與重建工作提供理論依據。

1 研究區概況

微山縣位于山東省濟寧市南部,地處116°34′—117°24′E,34°27′—35°20′N,南北相距120 km,東西相距8—30 km,總面積1779.8 km2；其中湖面面積1266 km2,占全縣總面積的三分之二,由北到南依次為南陽湖、獨山湖、昭陽湖、微山湖,統稱南四湖,占全省淡水量的45%。區內轄10個鎮,2個鄉和1個縣經濟開發區(2014年行政區劃)。在選取研究范圍時,主要選擇煤礦較集中、地表塌陷較嚴重的地區,同時考慮研究區的連貫性,最終選取微山縣的11個鄉鎮作為研究區,總面積1176.86 km2。研究區地理位置及范圍見圖1。

圖1 研究區地理位置及范圍Fig.1 The location and scope of the research area

2 數據來源與研究方法

2.1 數據來源與處理

選取2000年8月21日、2005年9月4日、2010年9月18日的Landsat5 TM遙感影像和2016年9月2日的Landsat8 OLI遙感影像作為數據源,其云量均小于2%,空間分辨率為30 m；影像的時間間隔大致為5年左右,以保證礦區土地利用生態風險評估時間粒度的一致性；同時選用降雨量較充沛、植被覆蓋度較高的9月份左右的遙感影像,便于礦區土地利用信息的提取。其他輔助數據包括：微山縣土地利用現狀圖、微山縣礦區分布圖、微山縣2014—2016年國民經濟和社會發展統計公報。

參照《土地利用現狀分類和編碼》(GB-T21010—2007),考慮研究區土地的用途、利用方式、覆蓋特征和遙感影像的分辨率,將研究區景觀類型分為耕地、其他農用地、城鄉建設用地、水域、塌陷積水區、灘涂沼澤6類。利用ENVI 5.3,采用監督分類中的支持向量機分類、人工目視判讀與決策樹分類、NDWI、NDVI相結合的分層分類,得到2000、2005、2010、2016年4期土地利用類型分布圖(圖2)。在4期分類圖上分別隨機選取200個檢查點,通過實地調查并結合相近年份的土地利用現狀圖和天地圖網站上的高清衛星圖,獲取檢查點的實際土地利用類型,利用ENVI軟件計算混淆矩陣和Kappa系數,經驗證Kappa系數均在85%以上,能滿足本研究的需要。

圖2 2000—2016年研究區土地利用類型分布圖Fig.2 Distribution maps of land use type in the research area from 2000 to 2016

2.2 采樣方法

運用ArcGIS中的Fishnet(漁網)工具創建0.5 km×0.5 km、1 km×1 km、1.5 km×1.5 km、2 km×2 km的正方形規則格網(樣區),其樣區個數分別為：4710個、1172個、521個、296個,根據公式計算每一個樣區的土地利用生態風險指數值,并將其作為樣區中心點的值[18-19,23-29],從而實現生態風險值的空間化表達,以更好地對研究區土地利用生態風險的空間異質性進行研究分析。

2.3 土地利用生態風險指數的構建

2.3.1景觀干擾度指數Ei

不同景觀類型在保護物種、維護生物多樣性、保持生態系統結構和功能完整性、促進景觀結構自然演替等方面的作用是有差別的,同時會受到外界環境的干擾,而不同景觀類型對外界環境的抗干擾能力也是不同的。參考相關文獻[17-18,23-24,27-29],選取景觀破碎度指數、景觀分離度指數和景觀優勢度指數,通過這三個指數的疊加構建景觀干擾度指數Ei,以反映不同景觀類型所代表的生態系統受到干擾(主要是人類開發活動)的程度,其表達式為：

Ei=aCi+bSi+cDOi

(1)景觀破碎度指數Ci：景觀破碎化是由于自然或人為干擾所導致的景觀由單一、均質和連續的整體趨向于復雜、異質和不連續的斑塊鑲嵌體的過程,景觀破碎度是景觀異質性的重要組成部分,破碎度值越大,表明景觀單元內部穩定性越低。公式為：

Ci=ni/A

式中,Ci為景觀類型i的破碎度,ni為景觀類型i的斑塊數,A為景觀(一個樣區)的總面積。

(2)景觀分離度指數Si：是指某一景觀類型中不同元素或斑塊個體分布的分離程度,值越大,表明景觀在地域分布上越分散,景觀分布越復雜[17],受到的干擾程度越大。

Si=Di/Pi

式中,Si為景觀類型i的分離度,Di為景觀類型i的距離指數,Pi為景觀類型i的面積指數,Ai為一個樣區中景觀類型i的面積,ni、A含義同上。

(3)景觀優勢度指數DOi：景觀優勢度表征景觀結構中某一類型支配景觀的程度,反映了該景觀類型對景觀格局形成和變化影響的大小[17]。公式為：

DOi=(1+景觀類型i的密度+景觀類型i的面積指數)/3

其中,密度=(一個樣區中景觀類型i的斑塊數/此樣區中總的斑塊數),景觀類型i的面積指數同上。

根據以上公式計算出Ci、Si、DOi指標后,由于量綱不同,需進行歸一化處理。a、b、c為各景觀指數的權重,且三者相加為1。權重值的大小反映各景觀指數解釋景觀類型受干擾的程度,根據研究區實際情況,并結合前人研究成果[17-18,23-24,26-29],認為破碎度指數最重要,其次為分離度指數和優勢度指數,三個指數分別賦值為0.502、0.301、0.197。

2.3.2景觀脆弱度指數Fi

景觀脆弱度表示不同生態系統的易損性,與其在景觀自然演替過程中所處的階段有關。一般情況下,生態系統處于初級演替階段時,食物鏈結構簡單、生物多樣性指數小,其較為脆弱。在土地利用生態系統中,人類活動是主要的干擾因素之一,土地利用程度不僅反映土地本身的自然屬性,同時也反映人為因素與自然環境因素的綜合效應。參考相關文獻[11,17-18,23-26]并根據研究區實際情況,本研究區6種景觀類型,以塌陷積水區最為脆弱,其次是其他農用地,最穩定的是城鄉建設用地,6種景觀類型分別賦值為：塌陷積水區為6、其他農用地為5、灘涂沼澤為4、耕地為3、水域較穩定為2、城鄉建設用地為1。采用反正切函數歸一化方法處理[11,17,23-25]之后,得到各景觀類型的脆弱度指數Fi值分別為：0.8949、0.8743、0.8440、0.7952、0.7048、0.5。

2.3.3土地利用生態風險指數ERI

基于上述所建立的景觀干擾度指數和景觀脆弱度指數,結合各景觀類型的面積比重構建土地利用生態風險指數,用于描述一個評價單元內綜合生態損失度的相對大小,以便通過采樣方法將景觀空間結構轉化為空間化的生態風險變量。土地利用生態風險指數ERI計算公式如下：

式中,ERIk為第k個采樣區內土地利用生態風險指數值,Aki為第k個采樣區內土地利用類型i面積,Ak為第k個采樣區面積,Ei為土地利用類型i的干擾度指數,Fi為相應的脆弱度指數。

2.4 空間統計分析

2.4.1空間自相關分析

空間統計學是以具有地理空間信息特性的事物或現象的空間相互作用及變化規律為研究對象,以區域化變量為基礎,研究既具有隨機性又具有結構性,或空間相關性和依賴性的一門科學[30-31],其核心是空間自相關性、空間依賴性和空間異質性。

所謂空間自相關性是指空間上越靠近的事物或現象越相似[4]?？臻g自相關的基本度量是空間自相關系數,可用全局和局部兩種指標來度量。

全局空間自相關(Global Spatial Autocorrelation)用于描述區域單元某種屬性值的整體分布狀況,判斷該屬性值在空間上是否存在聚集性的特點。常用的指數是全局Moran′sI指數。

全局Moran′sI指數用來評估整個研究區域內所有空間對象是集聚分布、離散分布還是隨機分布。其計算公式為[30-35]：

2.4.2地統計分析

(1)半變異函數

本研究借助地統計學中的半變異函數深入探討景觀結構變化等造成的土地利用生態風險的時空變化情況和空間變異規律。半變異函數的數學公式[4,32-36]為：

式中,r(h)是變異函數,h是兩樣本配對抽樣的分隔距離,即步長,Z(xi)和Z(xi+h)分別是隨機變量Z在空間位置xi和xi+h上的取值,N(h)為分隔距離為h時的樣本點總對數。

以r(h)為縱坐標,h為橫坐標做圖,得到半方差圖(semivariogram)。半方差圖中包含3個重要參數：1)塊金值(nugget)C0,表示區域化變量在小于抽樣尺度時的非連續性變異[35],2)基臺值(sill)C+C0,即平穩值,描述變量在研究區域范圍中總的空間變異程度,基臺值包括塊金值(C0)和結構方差(C)兩部分。3)變程(range)A,表示極限距離,即某特征在空間上自相關的空間幅度。

(2)克里格插值

如果變異函數和相關分析的結果表明區域化變量存在空間相關性,就可以利用區域化變量的原始數據和變異函數的結構特點,對未采樣點的區域化變量的取值進行線性無偏、最優估計[35]。

3 結果與分析

3.1 土地利用類型結構變化特征

由表1及圖3可知,17年間研究區耕地面積呈波動式減少,2000—2005年面積減少4604.94 hm2,2005—2010面積又有所增加,至2016年耕地面積為28706.85 hm2；城鄉建設用地2000—2005年增幅最大,增加了2830.41 hm2,其后面積呈緩慢增加趨勢,2010—2016年間建設用地減少了804.96 hm2,17年間建設用地共增加4589.19 hm2；塌陷積水區2000—2005年面積增幅最大,之后有所緩解,2010—2016年面積減少了132.57 hm2,整個研究期內塌陷積水區面積共增加473.04 hm2；其他農用地和水域面積總體呈增加趨勢,灘涂沼澤呈減少趨勢。

表1 2000—2016年研究區土地利用類型面積變化

圖3 2000—2016年研究區土地利用類型面積變化統計圖 Fig.3 Statistic map of land use type area change in the research area from 2000 to 2016

3.2 空間自相關分析

根據空間概念化模型中的“共享邊或角”規則建立空間權重矩陣,并進行行標準化,計算四種尺度下土地利用生態風險的全局Moran′sI指數(表2),結合P值和z得分,進行自相關顯著性檢驗[37]。以1 km×1 km尺度下的土地利用生態風險指數值為例,首先進行零假設驗證,以正態分布95%置信區間雙側檢驗閾值1.96為臨界值,根據表2,研究期的P值均小于0.05,z得分均大于1.96,說明研究區土地利用生態風險的空間分布不是隨機模式,存在空間自相關性。1 km×1 km尺度下2000年、2005年、2010年、2016年的Moran′sI指數分別為0.4576、0.4950、0.4279、0.4677,說明土地利用生態風險的空間分布呈集聚分布模式,具有很強的空間正相關性。

表2 1 km×1 km尺度下全局Moran′s I 指數及驗證值

3.3 土地利用生態風險的尺度效應分析

空間異質性依賴于尺度(粒度和幅度),粒度和幅度對空間異質性的測量和理解有重要的影響[4]。為了使建立的半變異函數模型能準確地反映各種尺度上的變化特征,要確定最佳的采樣尺度[4,38]。根據土地利用生態風險指數的構成,以2005年土地利用類型分布圖為試驗數據,選取破碎度指數,并計算研究區Shannon多樣性指數值,借助GS+地統計分析軟件,進行不同尺度下的破碎度指數和Shannon多樣性指數半變異函數模型擬合(表3、表4),以確定土地利用生態風險空間變異的最佳研究尺度。

表3 不同尺度下破碎度指數半變異函數擬合模型參數

根據表3,0.5 km尺度下,破碎度指數的C0為0.0014,C+C0為0.0040,相比較1 km、1.5 km和2 km尺度下的C0和C+C0,其值最小,說明該幅度內尺度效應不明顯；從0.5 km到1.5 km,C0從0.0014上升到0.0056,說明人類活動等隨機因素引起的空間異質性在逐漸增強[4],1.5 km到2 kmC0從0.0056驟然下降到0.0004,尺度的增加導致破碎度發生了明顯變化。C+C0從0.5 km到1 km上升最快,1 km到2 km緩慢上升,尺度的增加會降低空間自相關性。0.5 km到1 kmC/(C+C0)上升最快,之后緩慢上升,也反映了空間相關性在逐漸減弱。綜上可知,0.5 km尺度過小,1.5 km到2 km空間自相關性明顯減弱,1 km是適宜的研究尺度。

表4 不同尺度下Shannon多樣性指數半變異函數擬合模型參數

分析不同尺度下多樣性指數的半變異函數模型擬合參數(表4),尺度由0.5 km上升到1.5 km時,C0從0.0228下降到0.0206,到2 km時C0降至最低值0.0183。根據半變異函數的理論,塊金值較大,說明較小尺度上的某種過程不可忽視,隨機因素引起的空間異質性起主要作用。0.5 km到1 kmC+C0增加至0.0773,1 km到2 kmC+C0由0.0773降低至0.0633,說明隨著尺度的增加,多樣性指數的變化逐漸減弱。0.5 km到1 kmC/(C+C0)由0.606上升至0.730,說明0.5 km到1 km,多樣性指數的空間自相關性逐漸增強；1 km到2 kmC/(C+C0),波動式下降,由結構性因素所引起的空間異質性在不斷減弱。因此,1km能較直觀的反映多樣性指數的空間變異情況。

綜上,在進行微山縣11個鄉鎮的土地利用生態風險時空異質性研究時,1 km是最佳研究尺度。

3.4 適宜尺度下土地利用生態風險的半方差分析

在GS+中,分別對4期1 km×1 km的樣點數據進行半變異函數值的計算,并進行土地利用生態風險半變異函數模型擬合(表5及圖4),其中步長設置為格網的間距1 km,有效滯后距離39 km為采樣點最大間距的1/2。

表5 2000年—2016年土地利用生態風險半變異函數模型擬合參數

由表5及圖4看出,研究區2000年和2016年土地利用生態風險的最優擬合模型為球狀模型,2005年和2010年的最優擬合模型為指數模型。4個時期的土地利用生態風險的空間變異特征主要表現在以下幾方面：1)從各個年份最優擬合模型中塊金值C0的變化來看,2016年C0最大,表明土地利用生態風險受人類活動等隨機因素引起的空間異質性較大,較小尺度上的某種過程不能忽視。2)基臺值C+C0,從2000年的0.01000上升至2005年的最大值0.01276,到2016年有所降低,說明土地利用生態風險度空間變異程度隨時間的推移在不斷增強,其中2005年基臺值相對突出,表明這個時段影響土地利用生態風險度的某些因子的空間變異性明顯增強,2005年土地利用生態風險總的空間變異程度較高。3)結構方差與基臺值的比值C/(C+C0),2000、2005、2010、2016年分別為57.5%、68.8%、67.9%和54.5%,均大于50%,說明氣溫、降水、地形地貌、土壤類型等空間結構性因素引起的空間異質性占主導地位,從2000年到2005年,其值先逐漸增大至最大值68.8%,說明研究區土地利用生態風險的空間分異從2000年到2005年受結構性因素的影響逐漸加大,生態風險指數在小尺度上的隨機變異逐漸被較大尺度上的空間結構性變異所取代。2005—2016年呈緩慢下降趨勢,表明隨著社會經濟的不斷發展,人類活動對土地利用干擾程度的不斷加深,導致2005年到2016年土地利用生態風險由隨機因素引起的空間變異程度在增強。4)變程A表示土地利用生態風險在空間上自相關的空間幅度,在大于該變程的空間范圍內,土地利用生態風險的空間自相關性消失,2000年、2005年、2010年、2016年土地利用生態風險的變程分別為13 km、23.16 km、12.75 km、15.19 km,由2005年變程值較大,也可以看出土地利用生態風險值的變化2005年較其他年份劇烈。

圖4 2000—2016年研究區土地利用生態風險半變異函數模型擬合曲線圖Fig.4 The fitting curve of semivariogram model of land use ecological risk from 2000 to 2016 in the research area

3.5 土地利用生態風險的克里格插值及精度評定

利用GS+完成2000—2016年土地利用生態風險半變異函數最優模型擬合,獲取各模型相關參數之后,結合ArcGIS地統計分析模塊,在步長1 km下,選擇普通克里格插值方法進行空間估值,通過交叉驗證進行插值精度評定(表6)。通過驗證分析,誤差平均值、誤差標準平均值最接近于0,均方根預測誤差最小,平均標準誤差最接近于均方根預測誤差,標準均方根預測誤差接近于1,預測誤差皆符合要求,預測精度較高。

表6 插值預測誤差統計表

3.6 土地利用生態風險的時空異質性分析

將插值結果轉化為柵格格式,并按照研究區范圍進行不規則裁剪。通過ArcGIS數據處理,發現4期土地利用生態風險預測值的范圍都介于0.2718—0.5814值域范圍內,為了便于比較4期土地利用生態風險的空間分布情況,利用ArcGIS軟件中的重分類工具,進行土地利用生態風險等級劃分,具體劃分5個等級：低生態風險區(ERI<0.32)、中低生態風險區(0.32≤ERI<0.43)、中生態風險區(0.43≤ERI<0.48)、中高生態風險區(0.48≤ERI<0.5)、高生態風險區(ERI≥0.5),最終得到2000、2005、2010和2016年四期土地利用生態風險等級圖(圖5)。

圖5 2000—2016年研究區土地利用生態風險空間分布圖Fig.5 Spatial distribution map of land use ecological risk from 2000 to 2016 in the research area

(1)時間變化

從圖5及表7、表8可知,2000年研究區以中低生態風險區和中生態風險區為主,占總面積的57.58%,主要分布在微山湖湖區東西兩側,呈集中連片分布,另外獨山湖、昭陽湖湖區及周邊也以中低生態風險和中生態風險為主,景觀類型主要為湖泊水面及其他農用地類型中的臺田魚塘；中高生態風險區和高生態風險區面積相對較少,分別占區域總面積的17.32%和17.36%,呈散列式分布,主要分布于留莊鎮東南部、歡城鎮北部和中部一小部分、傅村街道中東部、昭陽街道中東部、韓莊鎮東部、高樓鄉西北部及西南部的帶狀區域及西平鄉的北部,對照同時期的土地利用類型分布圖,高生態風險區主要分布在城鄉建設用地及工礦企業周邊、道路沿線,另外因煤礦塌陷造成零散分布的耕地其生態風險值也較高。低生態風險區面積9125.48 hm2,比例最少,主要分布在微山湖湖心區域。

表7 土地利用生態風險各年面積及比例

2005年以中低生態風險區和中生態風險區為主,占區域總面積的60.66%,其空間分布與2000年的大致一致,面積比例與2000年相比,中生態風險區面積比減少了8.90%,中低生態風險區面積比增加了11.98%。中高生態風險區和高生態風險區面積相比2000年都有所減少,分別減少了8.36%和3.50%,低生態風險區面積19476.69 hm2,占區域總面積的16.51%,相比2000年面積比增加了8.77%?？傮w上2000—2005年,土地利用生態風險由中、中高、高生態風險向中低、低生態風險轉變。

2010年以中低生態風險區和中生態風險區為主,占區域總面積的74.71%,和2005年相比中低生態風險區和中生態風險區面積分別增加了9.06%和4.99%,增幅較?。坏蜕鷳B風險區面積和2005年相比變化不大,基本保持穩定；中高生態風險區和高生態風險區面積相比2005年有所下降,特別是高生態風險區面積降幅較大,下降了10.60%。2005—2010年土地利用生態風險由中高、高生態風險向中、中低生態風險轉變。

表8 土地利用生態風險各年面積變化及比例

2016年研究區也以中低生態風險區和中生態風險區為主,占區域總面積的70.51%,和2010年相比,中低生態風險區面積有所減少,減少了14%,但中生態風險區面積增幅較大,增加了9.8%；低生態風險區面積和2010年相比減少了8.51%；中高生態風險區和高生態風險區面積相比2010年又有所升高,特別是中高生態風險區面積增幅較大,增加了9.17%。2010—2016年,土地利用生態風險由低、中低生態風險向中、中高和高生態風險轉變。

總體上2000—2016年,土地利用生態風險呈現出由中高、高生態風險向中、中低生態風險轉變的特征,其中高生態風險區面積降幅最大,下降了10.56%,中低生態風險區面積增幅最大,增加了7.05%,其次是中生態風險區,面積增加了5.89%,低生態風險區和中高生態風險區面積變化幅度不大。

(2)空間變化

為了更好地分析2000—2016 年間研究區土地利用生態風險的空間異質性演變規律,在ArcGIS中運用Spatial Analyst模塊中的Raster Calculator工具,將2000、2005、2010和2016年土地利用生態風險空間分布(圖5)進行空間運算,獲取2000—2005年、2005—2010年、2010—2016年、2000—2016年的土地利用生態風險空間變化圖(圖6),將差值運算結果等于0的區域定義為生態風險“穩定”區,等于-1的區域定義為生態風險“改善”區,等于-2、-3、-4的區域定義為生態風險“明顯改善”區,等于1的區域定義為生態風險“惡化”區,等于2、3、4的區域定義為生態風險“明顯惡化”區。

根據圖6及表9,2000—2005年研究區土地利用生態風險程度較低,生態系統較穩定,生態風險穩定區和改善區面積占區域總面積的78.92%,且分布較集中,主要分布在獨山湖、昭陽湖、微山湖湖區及周邊,景觀類型主要為湖泊水域,這部分區域景觀連通性較高,人為干擾性較低,生態系統較穩定；生態風險惡化區和明顯惡化區分布較分散,主要分布在留莊鎮主城區以南新安煤礦周邊,歡城鎮作為煤礦較集中的鄉鎮,其生態風險惡化及明顯惡化的區域較大,主要集中在蔣莊煤礦、岱莊生建煤礦,柴里煤礦等礦區及周邊,另外傅村街道東南部、昭陽街道東北部塌陷積水較嚴重的區域、趙廟鎮孔莊煤礦及周邊生態風險惡化程度也較高,生態風險惡化區和明顯惡化區主要位于礦區周邊及道路沿線,這些區域土地利用人為干擾因素較強且離湖區較近,地面塌陷較嚴重,景觀類型的破碎化程度較高。

圖6 2000—2016年研究區土地利用生態風險空間變化分布圖Fig.6 Spatial distribution map of land use ecological risk in the research area from 2000 to 2016

2005—2010年生態風險明顯惡化區比例降低,穩定區比例較2000—2005年增幅較大,增加了10.52%,且煤礦較集中的歡城鎮、傅村街道和昭陽街道生態風險有了大幅改善,主要原因是2005年以來濟寧市政府開展了大規模的塌陷地復墾治理工程,微山縣作為采煤塌陷較嚴重的區域,也開展了大規模的復墾治理工程,經過5年多土地復墾,塌陷地的生態環境狀況有了一定的改善,耕地經過整治后變得集中連片,塌陷積水區整治后變成規整的魚塘等水產養殖用地。

表9 土地利用生態風險時空面積變化及比例

2010—2016年生態風險惡化區面積比例大幅增加,主要分布于獨山湖、昭陽湖等湖周邊的灘涂沼澤,另外歡城鎮中東部、傅村街道東部、夏鎮街道東北部、韓莊鎮東北部等煤礦集中區生態風險惡化區也有成片分布。主要原因是2014年南四湖遭遇了罕見的旱災,水面和濕地大面積萎縮導致水質緩沖能力差,持續的旱情使獨山湖、昭陽湖、微山湖等湖區及周邊生態環境遭到嚴重破壞；同時隨著科學技術的發展及采煤機械化程度的提高,礦山的回采率逐步提高,煤炭資源的大規模開采,產生了越來越嚴重的地面塌陷、矸石堆積、水體污染等生態環境問題。

4 討論與建議

(1)在進行微山縣11個鄉鎮的土地利用生態風險時空異質性研究時,1 km是最佳研究尺度,在該尺度下,土地利用生態風險的空間分布呈集聚模式,且具有很強的空間正相關性。2000—2005年生態風險指數在小尺度上的隨機變異逐漸被較大尺度上的空間結構性變異所取代,2005年土地利用生態風險總的空間變異程度較高,2005—2016年,隨著社會經濟的不斷發展,人類活動對土地利用干擾程度的不斷加深,土地利用生態風險由隨機因素引起的空間變異程度在增強。

(2)2000—2016年,土地利用生態風險呈現由中高、高生態風險向中、中低生態風險轉變的特征。其中高生態風險區面積降幅最大,中低生態風險區面積增幅最大,其次增幅較大的是中生態風險區,低生態風險區和中高生態風險區面積變化幅度不大。17年間土地利用生態風險有了大幅降低,生態風險穩定區和改善區占區域總面積的79.21%,生態風險惡化區和明顯惡化區主要分布在東部礦區周邊及道路沿線等人類活動較頻繁、地面塌陷較嚴重的區域。

根據以上結論,結合《濟寧市采煤塌陷地治理總體規劃(2009—2020年)》,在后續的土地利用過程中,宜采取以下治理措施：1)對于東部煤礦集中區,采煤塌陷較嚴重,這些區域應根據地形、土壤、積水區深淺及地表塌陷的實際狀況因地制宜地采取復墾治理措施。對于塌陷較淺的區域,應充分發揮南四湖湖底淤泥養分充足的資源優勢,在保護河湖生態平衡的前提下,對沿湖(河)礦井采用引湖填充方法復墾塌陷地,并最大限度的恢復為耕地,增加農用地面積；對于塌陷較深,地面坡度較大的區域,應采取挖深墊淺的治理措施,整理成魚塘,用于水產養殖業。2)在研究區西南部微山湖沿岸東西兩側,耕地的生態風險惡化程度2010—2016年較高,對此,在保證耕地總量動態平衡的前提下,應積極實行退耕還湖、退耕還濕的治理措施,增加蘆葦沼澤濕地的面積,從而保護生物物種的多樣性,提高生態系統的穩定程度。

本研究基于景觀生態學構建土地利用生態風險指數,但在景觀干擾度指數的權重賦值時,沒有進一步分析社會、經濟因素與景觀干擾度指數值間的定量關系,只是概念性的進行綜合性度量,所以計算出的土地利用生態風險指數值也是相對的,并不具有絕對性；另外,在進行研究區土地利用生態風險半變異函數模型擬合時,雖然考慮了土地利用生態風險的尺度效應,但缺乏對生態風險的變化方向效應分析,未來可在這兩個方面開展進一步的研究。