礦區(qū)土地利用生態(tài)風(fēng)險的時空異質(zhì)性

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

1 山東農(nóng)業(yè)大學(xué)信息科學(xué)與工程學(xué)院,泰安 271018 2 山東農(nóng)業(yè)大學(xué)資源與環(huán)境學(xué)院,泰安 271018 3 山東省濟(jì)寧市國土資源局,濟(jì)寧 272001

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

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

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

1 研究區(qū)概況

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

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

2 數(shù)據(jù)來源與研究方法

2.1 數(shù)據(jù)來源與處理

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

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

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

2.2 采樣方法

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

2.3 土地利用生態(tài)風(fēng)險指數(shù)的構(gòu)建

2.3.1景觀干擾度指數(shù)Ei

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

Ei=aCi+bSi+cDOi

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

Ci=ni/A

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

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

Si=Di/Pi

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

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

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

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

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

2.3.2景觀脆弱度指數(shù)Fi

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

2.3.3土地利用生態(tài)風(fēng)險指數(shù)ERI

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

式中,ERIk為第k個采樣區(qū)內(nèi)土地利用生態(tài)風(fēng)險指數(shù)值,Aki為第k個采樣區(qū)內(nèi)土地利用類型i面積,Ak為第k個采樣區(qū)面積,Ei為土地利用類型i的干擾度指數(shù),Fi為相應(yīng)的脆弱度指數(shù)。

2.4 空間統(tǒng)計分析

2.4.1空間自相關(guān)分析

空間統(tǒng)計學(xué)是以具有地理空間信息特性的事物或現(xiàn)象的空間相互作用及變化規(guī)律為研究對象,以區(qū)域化變量為基礎(chǔ),研究既具有隨機(jī)性又具有結(jié)構(gòu)性,或空間相關(guān)性和依賴性的一門科學(xué)[30-31],其核心是空間自相關(guān)性、空間依賴性和空間異質(zhì)性。

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

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

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

2.4.2地統(tǒng)計分析

(1)半變異函數(shù)

本研究借助地統(tǒng)計學(xué)中的半變異函數(shù)深入探討景觀結(jié)構(gòu)變化等造成的土地利用生態(tài)風(fēng)險的時空變化情況和空間變異規(guī)律。半變異函數(shù)的數(shù)學(xué)公式[4,32-36]為：

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

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

(2)克里格插值

如果變異函數(shù)和相關(guān)分析的結(jié)果表明區(qū)域化變量存在空間相關(guān)性,就可以利用區(qū)域化變量的原始數(shù)據(jù)和變異函數(shù)的結(jié)構(gòu)特點,對未采樣點的區(qū)域化變量的取值進(jìn)行線性無偏、最優(yōu)估計[35]。

3 結(jié)果與分析

3.1 土地利用類型結(jié)構(gòu)變化特征

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

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

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

3.2 空間自相關(guān)分析

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

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

3.3 土地利用生態(tài)風(fēng)險的尺度效應(yīng)分析

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

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

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

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

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

綜上,在進(jìn)行微山縣11個鄉(xiāng)鎮(zhèn)的土地利用生態(tài)風(fēng)險時空異質(zhì)性研究時,1 km是最佳研究尺度。

3.4 適宜尺度下土地利用生態(tài)風(fēng)險的半方差分析

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

表5 2000年—2016年土地利用生態(tài)風(fēng)險半變異函數(shù)模型擬合參數(shù)

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

圖4 2000—2016年研究區(qū)土地利用生態(tài)風(fēng)險半變異函數(shù)模型擬合曲線圖Fig.4 The fitting curve of semivariogram model of land use ecological risk from 2000 to 2016 in the research area

3.5 土地利用生態(tài)風(fēng)險的克里格插值及精度評定

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

表6 插值預(yù)測誤差統(tǒng)計表

3.6 土地利用生態(tài)風(fēng)險的時空異質(zhì)性分析

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

圖5 2000—2016年研究區(qū)土地利用生態(tài)風(fēng)險空間分布圖Fig.5 Spatial distribution map of land use ecological risk from 2000 to 2016 in the research area

(1)時間變化

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

表7 土地利用生態(tài)風(fēng)險各年面積及比例

2005年以中低生態(tài)風(fēng)險區(qū)和中生態(tài)風(fēng)險區(qū)為主,占區(qū)域總面積的60.66%,其空間分布與2000年的大致一致,面積比例與2000年相比,中生態(tài)風(fēng)險區(qū)面積比減少了8.90%,中低生態(tài)風(fēng)險區(qū)面積比增加了11.98%。中高生態(tài)風(fēng)險區(qū)和高生態(tài)風(fēng)險區(qū)面積相比2000年都有所減少,分別減少了8.36%和3.50%,低生態(tài)風(fēng)險區(qū)面積19476.69 hm2,占區(qū)域總面積的16.51%,相比2000年面積比增加了8.77%。總體上2000—2005年,土地利用生態(tài)風(fēng)險由中、中高、高生態(tài)風(fēng)險向中低、低生態(tài)風(fēng)險轉(zhuǎn)變。

2010年以中低生態(tài)風(fēng)險區(qū)和中生態(tài)風(fēng)險區(qū)為主,占區(qū)域總面積的74.71%,和2005年相比中低生態(tài)風(fēng)險區(qū)和中生態(tài)風(fēng)險區(qū)面積分別增加了9.06%和4.99%,增幅較??；低生態(tài)風(fēng)險區(qū)面積和2005年相比變化不大,基本保持穩(wěn)定；中高生態(tài)風(fēng)險區(qū)和高生態(tài)風(fēng)險區(qū)面積相比2005年有所下降,特別是高生態(tài)風(fēng)險區(qū)面積降幅較大,下降了10.60%。2005—2010年土地利用生態(tài)風(fēng)險由中高、高生態(tài)風(fēng)險向中、中低生態(tài)風(fēng)險轉(zhuǎn)變。

表8 土地利用生態(tài)風(fēng)險各年面積變化及比例

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

總體上2000—2016年,土地利用生態(tài)風(fēng)險呈現(xiàn)出由中高、高生態(tài)風(fēng)險向中、中低生態(tài)風(fēng)險轉(zhuǎn)變的特征,其中高生態(tài)風(fēng)險區(qū)面積降幅最大,下降了10.56%,中低生態(tài)風(fēng)險區(qū)面積增幅最大,增加了7.05%,其次是中生態(tài)風(fēng)險區(qū),面積增加了5.89%,低生態(tài)風(fēng)險區(qū)和中高生態(tài)風(fēng)險區(qū)面積變化幅度不大。

(2)空間變化

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

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

圖6 2000—2016年研究區(qū)土地利用生態(tài)風(fēng)險空間變化分布圖Fig.6 Spatial distribution map of land use ecological risk in the research area from 2000 to 2016

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

表9 土地利用生態(tài)風(fēng)險時空面積變化及比例

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

4 討論與建議

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

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

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

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