基于土地利用變化的河北省壩上地區(qū)景觀生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)

劉孟竹，王彥芳，裴宏偉

(1.河北建筑工程學(xué)院 市政與環(huán)境工程系，河北 張家口 075000;2.河北地質(zhì)大學(xué)，河北 石家莊 050031)

土地生態(tài)系統(tǒng)作為地區(qū)生態(tài)環(huán)境的載體，其長(zhǎng)期的積累性變化是全球環(huán)境問(wèn)題的一大誘因[1]，其中土地生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)作為生態(tài)安全問(wèn)題里的重點(diǎn)內(nèi)容[2]，國(guó)內(nèi)外學(xué)者就其已展開(kāi)了廣泛研究。目前，針對(duì)已有的研究成果分析，相關(guān)研究在空間尺度選擇上多為市、區(qū)、縣域甚至小型湖泊流域，時(shí)間跨度從以往幾年到幾百年不等。例如，高賓[3]等基于景觀尺度分析了錦州灣沿海經(jīng)濟(jì)開(kāi)發(fā)區(qū)15 a來(lái)的生態(tài)風(fēng)險(xiǎn)水平，張學(xué)斌[4]等以遙感影像為基礎(chǔ)，構(gòu)建了石羊河流域的景觀生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)機(jī)制。從國(guó)內(nèi)學(xué)者廣泛地應(yīng)用該評(píng)價(jià)機(jī)制來(lái)看，通過(guò)景觀指數(shù)構(gòu)建的區(qū)域生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)體系提供了一個(gè)區(qū)域生態(tài)評(píng)價(jià)的新視角，是切實(shí)可行并具有參考意義的，但同時(shí)該方法也無(wú)可避免地缺乏對(duì)具體生態(tài)過(guò)程的考量。此外，區(qū)域生態(tài)評(píng)價(jià)體系還可對(duì)區(qū)域內(nèi)社會(huì)經(jīng)濟(jì)、水文氣候、土地性質(zhì)等多角度對(duì)生態(tài)脅迫因素進(jìn)行定性選取[5]，通過(guò)主成分法、熵值法、層次分析法等確定各因素權(quán)重從而構(gòu)建一套完整的評(píng)價(jià)模型[6]，應(yīng)用較為廣泛的包括PSR，DPSIR，SSD等模型。針對(duì)區(qū)域景觀生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)的研究，普遍是針對(duì)過(guò)去的生態(tài)風(fēng)險(xiǎn)變化進(jìn)行分析，了解其變化趨勢(shì)以及驅(qū)動(dòng)機(jī)理仍不足以對(duì)風(fēng)險(xiǎn)進(jìn)行預(yù)警和有效管控，尤其在生態(tài)敏感區(qū)，無(wú)論是壩上地區(qū)還是其所屬的農(nóng)牧交錯(cuò)帶，這類(lèi)地區(qū)往往存在地廣人稀的人文特點(diǎn)，因而對(duì)全局的生態(tài)風(fēng)險(xiǎn)把控難以面面俱到?？焖佟⒕珳?zhǔn)地識(shí)別生態(tài)風(fēng)險(xiǎn)區(qū)域從而進(jìn)行管理，可以通過(guò)基于土地變化視角的景觀生態(tài)風(fēng)險(xiǎn)預(yù)警得以實(shí)現(xiàn)，CLUE-S，CA-Markov等模型均為土地預(yù)測(cè)模型并且得到了國(guó)內(nèi)外學(xué)者廣泛的應(yīng)用。因此，構(gòu)建風(fēng)險(xiǎn)預(yù)測(cè)模型對(duì)研究區(qū)以及北方農(nóng)牧交錯(cuò)帶生態(tài)環(huán)境類(lèi)似的區(qū)域，在生態(tài)建設(shè)、生態(tài)風(fēng)險(xiǎn)預(yù)警等均具有切實(shí)的指導(dǎo)意義。

河北省壩上地區(qū)作為北方農(nóng)牧交錯(cuò)帶典型的干旱半干旱過(guò)渡區(qū)，其生態(tài)環(huán)境極為脆弱，一旦被破壞其生態(tài)恢復(fù)與重建將難以逆轉(zhuǎn)。多年來(lái)壩上地區(qū)土地沙漠化嚴(yán)重，土壤風(fēng)蝕、水土流失、水域萎縮等生境問(wèn)題一再惡化[7]，同時(shí)該區(qū)極端天氣頻發(fā)，短歷時(shí)暴雨頻現(xiàn)。河北省壩上地區(qū)惡劣的氣候環(huán)境條件，也破壞了該區(qū)土地生態(tài)環(huán)境。作為首都生態(tài)屏障，壩上地區(qū)肩負(fù)的涵養(yǎng)水源、防風(fēng)固沙等生態(tài)功能與其現(xiàn)有的惡劣生境條件之間的矛盾非常尖銳[8]。而目前關(guān)于壩上地區(qū)生態(tài)的研究較多關(guān)注其生態(tài)環(huán)境問(wèn)題及其成因，時(shí)間也相對(duì)滯后，特定地針對(duì)當(dāng)下壩上地區(qū)的生態(tài)風(fēng)險(xiǎn)評(píng)估以及預(yù)測(cè)的研究卻不多見(jiàn)?；诖?，本文以河北省壩上地區(qū)為研究對(duì)象，結(jié)合土地利用及景觀生態(tài)學(xué)視角，對(duì)其近40 a來(lái)的土地利用動(dòng)態(tài)變化和生態(tài)風(fēng)險(xiǎn)進(jìn)行分析評(píng)價(jià)并對(duì)未來(lái)趨勢(shì)作出預(yù)測(cè)，以期為壩上地區(qū)生態(tài)建設(shè)和治理、可持續(xù)發(fā)展提供一定科學(xué)依據(jù)，同時(shí)壩上地區(qū)處于北方農(nóng)牧交錯(cuò)帶中段，其氣候條件以及農(nóng)、牧業(yè)交錯(cuò)的土地格局，可以作為農(nóng)牧交錯(cuò)帶最典型的代表區(qū)，該地區(qū)的研究可適用于整個(gè)北方農(nóng)牧交錯(cuò)帶地區(qū)。

1 材料與方法

1.1 研究區(qū)概況

壩上地區(qū)位于河北省西北部地帶，接壤內(nèi)蒙古高原南緣，地理位置介于114°35′—116°45′E和41°00′—42°20′N(xiāo)之間。其涵蓋的行政區(qū)域包含沽源縣、張北縣、康?？h的全區(qū)以及尚義縣、圍場(chǎng)縣、豐寧縣的部分區(qū)域，總面積約1.90×104km2。壩上地區(qū)屬大陸季風(fēng)高原氣候，終年以干冷、多風(fēng)氣候?yàn)橹?，年均氣溫?～2 ℃，多年平均降雨量在400 mm左右，70%的雨量集中在6—9月，年均蒸發(fā)量高達(dá)1 800 mm。該地區(qū)地貌以丘陵、平原居多研究區(qū)，整體海拔為834～2 229 m。壩上地區(qū)因其顯著的干旱多風(fēng)等氣候特點(diǎn)，加上全年少雨，地下水位下降、土地荒漠化、土壤風(fēng)蝕沙化逐年加劇。該區(qū)年均蒸發(fā)量顯著大于年降水量，地表水域現(xiàn)已大多萎縮甚至干涸[9]。近20 a壩上地區(qū)在市場(chǎng)經(jīng)濟(jì)的驅(qū)動(dòng)下農(nóng)業(yè)轉(zhuǎn)型開(kāi)始大面積種植經(jīng)濟(jì)效益較高的錯(cuò)季蔬菜，對(duì)地下水消耗量巨大，加上外界過(guò)量地下水開(kāi)采以及過(guò)度農(nóng)業(yè)灌溉，導(dǎo)致其地下水位近40 a來(lái)下降了3～5 m[10]，生態(tài)承載功能不堪重負(fù)，曾經(jīng)的華北第一大高原內(nèi)陸湖安固里淖早已干涸。壩上處于生態(tài)敏感區(qū)，與其所處北方農(nóng)牧交錯(cuò)帶面臨的生態(tài)問(wèn)題類(lèi)似，如今均水土流失、氣候惡劣、土地沙漠化趨于劇烈，給該區(qū)的生態(tài)環(huán)境恢復(fù)與重建帶來(lái)巨大挑戰(zhàn)。

1.2 數(shù)據(jù)源及處理

本研究所用的土地?cái)?shù)據(jù)來(lái)源于中國(guó)科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心，其均由Landsat衛(wèi)星影像目視解譯制成，分辨率為30 m。時(shí)間選取了1980，1990，2000，2010和2018年共計(jì)5期。該土地?cái)?shù)據(jù)由專(zhuān)業(yè)人員通過(guò)隨機(jī)選點(diǎn)后實(shí)地調(diào)查，數(shù)據(jù)精度在85%以上，能滿足本文研究需要。土地利用分類(lèi)根據(jù)《土地利用現(xiàn)狀分類(lèi)》標(biāo)準(zhǔn)并參考?jí)紊系貐^(qū)的實(shí)際情況將其分為6類(lèi)：耕地、林地、草地、水域、建設(shè)用地、未利用地。2026年預(yù)測(cè)數(shù)據(jù)以上述數(shù)據(jù)為基礎(chǔ)參數(shù)在IDRISI軟件中模擬所得，所有數(shù)據(jù)均統(tǒng)一為WGS 1984坐標(biāo)系，分辨率與原始數(shù)據(jù)保持一致。

1.3 方 法

1.3.1 CA-Markov模型預(yù)測(cè)模擬 元胞自助機(jī)(cellular automaton，CA)模型具有在時(shí)間、空間、狀態(tài)呈離散分布的特性，是一種在空間上相互作用、時(shí)間上又具有因果關(guān)系的網(wǎng)格動(dòng)力學(xué)模型，具有處理復(fù)雜空間系統(tǒng)的能力[11]。Markov模型基于Markov鏈的空間概率模型，因其較好的穩(wěn)定性和無(wú)后效性，能夠預(yù)測(cè)土地利用變化中各個(gè)時(shí)刻的變動(dòng)過(guò)程[12]。CA-Markov模型綜合了Markov模型的時(shí)序預(yù)測(cè)和CA模型的空間分布模擬，被廣泛應(yīng)用于土地變化的預(yù)測(cè)模擬中[13]。對(duì)于CA-Markov模型多數(shù)研究常結(jié)合二元Logistic回歸模型和MCE模型于構(gòu)成耦合模型用于土地利用預(yù)測(cè)中。以分布適宜性概率和多準(zhǔn)則約束作為模型預(yù)測(cè)模擬的轉(zhuǎn)換規(guī)則，其模擬結(jié)果也較理想。然而在實(shí)際中其驅(qū)動(dòng)因子數(shù)據(jù)來(lái)源不一，時(shí)效性差，空間量化不合理以及人為主觀對(duì)各因子權(quán)重設(shè)置的標(biāo)準(zhǔn)不統(tǒng)一往往給預(yù)測(cè)結(jié)果帶來(lái)更多不確定性。考慮到本文重點(diǎn)不在于預(yù)測(cè)，因此僅考慮未來(lái)的土地利用變化按照以往變化趨勢(shì)發(fā)展的自然情景。具體方法是在IDRISI軟件中通過(guò)Markov模型得到2000—2010年土地面積轉(zhuǎn)移矩陣和轉(zhuǎn)移概率矩陣作為CA-Markov模型模擬的初始參數(shù)進(jìn)行壩上地區(qū)2018年土地利用的預(yù)測(cè)模擬，模擬結(jié)果與2018實(shí)際分類(lèi)圖進(jìn)行疊置分析，利用像元正確預(yù)測(cè)個(gè)數(shù)占比與Kappa系數(shù)等指標(biāo)驗(yàn)證模型可信度。在模型可信度良好的情況下進(jìn)行壩上地區(qū)2026年土地預(yù)測(cè)。

1.3.2 土地利用時(shí)空變化分析 本研究采用綜合單一土地利用動(dòng)態(tài)度、土地利用轉(zhuǎn)移矩陣以及重心模型等指標(biāo)分別從土地利用變化的速率、數(shù)量和空間方向3個(gè)角度對(duì)壩上地區(qū)土地格局的變化進(jìn)行整體分析，以揭示其多年來(lái)土地結(jié)構(gòu)變化的特征以及未來(lái)變化的趨勢(shì)。

1.3.3 土地利用動(dòng)態(tài)度 土地利用動(dòng)態(tài)度能定量反映出土地變化的速度，其中單一土地利用動(dòng)態(tài)度側(cè)重研究期內(nèi)某種土地利用類(lèi)型數(shù)量的年變化率，綜合土地利用動(dòng)態(tài)度則是對(duì)研究區(qū)內(nèi)土地利用整體變化情況的刻畫(huà)。為深入細(xì)化分析多年來(lái)壩上地區(qū)土地利用的變化情況，本文選擇土地利用動(dòng)態(tài)度這一指標(biāo)，其計(jì)算公式[14]為：

(1)

式中：K為研究時(shí)段內(nèi)某一土地利用類(lèi)型動(dòng)態(tài)度；Ua，Ub分別為研究初期和研究末期某種土地利用類(lèi)型的數(shù)量；T為研究時(shí)長(zhǎng)。

1.3.4 土地利用轉(zhuǎn)移矩陣 土地轉(zhuǎn)移矩陣反映了某一區(qū)域在研究初期和研究末期各類(lèi)土地類(lèi)型面積互相轉(zhuǎn)化的動(dòng)態(tài)信息，不僅可以定量地表明不同土地利用類(lèi)型之間的轉(zhuǎn)化情況，還可以揭示不同土地利用類(lèi)型間的轉(zhuǎn)移速率。轉(zhuǎn)移矩陣表達(dá)式[15]為：

(2)

式中：S為面積；n為轉(zhuǎn)移前后土地利用類(lèi)型數(shù)；i，j(i，j=1，2，3，…，n)為轉(zhuǎn)移前后的土地類(lèi)型；Sij為轉(zhuǎn)移前的i類(lèi)土地轉(zhuǎn)換成轉(zhuǎn)移后j類(lèi)土地類(lèi)型的面積。

1.3.5 重心模型 土地格局的空間變化情況可以用各土地利用類(lèi)型的重心轉(zhuǎn)移變化來(lái)進(jìn)行表達(dá)。各土地類(lèi)型的重心坐標(biāo)一般用影像或者地圖的地理坐標(biāo)來(lái)表示，計(jì)算公式[16]為：

(3)

(4)

式中：Xt，Yt為第t年某種土地利用類(lèi)型分布重心的經(jīng)緯度坐標(biāo)；Cti為t年第i個(gè)區(qū)域土地利用類(lèi)型面積；Xi，Yi為第i個(gè)區(qū)域重心地理坐標(biāo)。

1.4 區(qū)域生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)

本研究采用景觀生態(tài)學(xué)角度構(gòu)建壩上地區(qū)域生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)模型。利用FRAGSTATS軟件計(jì)算出1980—2026年6期土地利用數(shù)據(jù)的各景觀總面積(CA)、各景觀斑塊數(shù)(NP)、斑塊密度(PD)以及景觀分維數(shù)(FRACT)4個(gè)指標(biāo)作為初始參數(shù)。在ArcMap軟件中對(duì)整個(gè)研究區(qū)進(jìn)行風(fēng)險(xiǎn)小區(qū)劃分，單元格大小根據(jù)平均斑塊面積的2～5倍原則進(jìn)行采樣[17]，最終確定采樣面積大小為4 km×4 km，整個(gè)壩上地區(qū)被劃分為1 399個(gè)風(fēng)險(xiǎn)小區(qū)。利用初始參數(shù)構(gòu)建的生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)模型對(duì)每一個(gè)風(fēng)險(xiǎn)小區(qū)進(jìn)行賦值計(jì)算，得到風(fēng)險(xiǎn)指數(shù)作為每一個(gè)單元格內(nèi)的參數(shù)值再經(jīng)過(guò)克里金(Kriging)插值法實(shí)現(xiàn)整個(gè)研究區(qū)生態(tài)風(fēng)險(xiǎn)空間可視化。最后采用空間自相關(guān)分析和進(jìn)行生態(tài)風(fēng)險(xiǎn)指數(shù)進(jìn)行空間分析。

1.4.1 區(qū)域生態(tài)風(fēng)險(xiǎn)評(píng)價(jià) 景觀生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)體系中，景觀指數(shù)法常以多尺度的土地利用數(shù)據(jù)為基礎(chǔ)進(jìn)行構(gòu)建，因其依托遙感數(shù)據(jù)而具有了與之同步的高度時(shí)空性，在生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)領(lǐng)域被廣泛應(yīng)用。其具體評(píng)價(jià)范式是將景觀生態(tài)風(fēng)險(xiǎn)定義為風(fēng)險(xiǎn)概率與景觀損失度的累乘[18]。本文基于已有研究中理論與方法較為成熟的生態(tài)損失度(Rs)模型構(gòu)建生態(tài)風(fēng)險(xiǎn)指數(shù)(ERIm)。其具體表達(dá)式[19]為：

(5)

Rs=D×Qs

(6)

Qs=aMs+bNs+cKs

(7)

式中：ERIm為第m個(gè)風(fēng)險(xiǎn)小區(qū)生態(tài)風(fēng)險(xiǎn)指數(shù)值；Am為第m個(gè)風(fēng)險(xiǎn)小區(qū)面積；Ams為m小區(qū)第s類(lèi)景觀面積；Rs為第s類(lèi)景觀損失度；D為景觀脆弱度，由專(zhuān)家賦值并歸一化處理所得；Qs為第s類(lèi)景觀干擾度；a，b，c為權(quán)重(0.5，0.3，0.2)；Ms，Ns，Ks分別為第s類(lèi)土地利用類(lèi)型景觀破碎度、景觀分離度、景觀分維數(shù)。

上述景觀指數(shù)的具體內(nèi)容和意義以及因子權(quán)重的分配在相關(guān)研究中已被準(zhǔn)確地定義并廣泛使用，此處不作贅述，具體可參考文獻(xiàn)[20]。

1.4.2 景觀生態(tài)風(fēng)險(xiǎn)空間分析 區(qū)域景觀生態(tài)風(fēng)險(xiǎn)指數(shù)作為一種空間變量常用地統(tǒng)計(jì)法對(duì)其空間分異特征進(jìn)行分析評(píng)估[21]。本文采用空間自相關(guān)分析法對(duì)生態(tài)風(fēng)險(xiǎn)指數(shù)在分布上的集聚性進(jìn)行分析并確定其與鄰域空間變量的依賴性。全局空間自相關(guān)與局部空間自相關(guān)常用的指數(shù)分別為Moran’sI指數(shù)與LISA指數(shù)。Moran’sI指數(shù)主要反映全部數(shù)據(jù)的某項(xiàng)屬性值在空間中的相關(guān)性大小，其值范圍為-1～1，其絕對(duì)值越接近1代表全局空間自相關(guān)性越強(qiáng)，等于0時(shí)不相關(guān)[21]。局部空間自相關(guān)LISA指數(shù)可以很好描述空間變量極值的局部空間聚集程度，其具體表現(xiàn)為：LL型(低低聚集)、HH型(高高聚集)、HL型(低值包高值聚集)、LH型(高值包低值聚集)4種類(lèi)型[22]。其對(duì)于局部空間異常特征能有效識(shí)別。上述指標(biāo)表達(dá)式[23]為：

(8)

(9)

2 結(jié)果與分析

2.1 壩上地區(qū)2026年土地模擬與驗(yàn)證

在IDRISI軟件中通過(guò)CA-Markov模型對(duì)壩上地區(qū)土地利用進(jìn)行預(yù)測(cè)模擬，結(jié)果如附圖4(見(jiàn)封3)所示。其中，對(duì)2018年土地利用預(yù)測(cè)模擬精度驗(yàn)證的結(jié)果詳見(jiàn)表1。經(jīng)混淆矩陣統(tǒng)計(jì)，2018年兩幅圖像中，像元正確預(yù)測(cè)個(gè)數(shù)占比88.74%，kappa系數(shù)為0.83，說(shuō)明模擬精度良好，該期預(yù)測(cè)結(jié)果和實(shí)際狀況一致性較高，模型具有可信度，可用于下一時(shí)期的模擬。壩上地區(qū)2026年土地利用預(yù)測(cè)結(jié)果如附圖4(見(jiàn)封3)所示，經(jīng)過(guò)上述精度驗(yàn)證，該期土地利用預(yù)測(cè)數(shù)據(jù)可作為本研究數(shù)據(jù)以作參考。

表1 壩上地區(qū)2018年土地利用預(yù)測(cè)模擬混淆矩陣 km2

2.2 土地利用時(shí)空數(shù)量變化分析

經(jīng)過(guò)ArcMap軟件對(duì)壩上地區(qū)1980—2026年土地利用分類(lèi)圖統(tǒng)計(jì)分析得到6期各土地類(lèi)型面積數(shù)據(jù)如表2所示。結(jié)合圖1分析，整體上看，1980—2018年壩上地區(qū)主要以耕地為主，其面積占比常年穩(wěn)居在50%左右，其在2000—2010年和2010—2018年有快速上升后下降的趨勢(shì)。草地面積居其次，近40 a來(lái)浮動(dòng)在22.61%～26.18%，在2000—2010年有一個(gè)較為顯著的下降，隨后8 a間處于較小幅度上升。林地在壩上地區(qū)覆蓋度不如草地高，其變化主要表現(xiàn)在2000—2018年期間較急劇的擴(kuò)張，在2018年覆蓋度達(dá)到最高水平，占全區(qū)面積的近19%。壩上地區(qū)水域在2000—2010年顯著萎縮，究其原因是該期間內(nèi)，位于張北縣的華北第一內(nèi)陸湖淖安固里淖全完干涸，造成了該區(qū)水域面積的斷崖式下降。壩上6縣居民住地、工礦用地、交通用地等面積在每一階段穩(wěn)步上升，2000—2018年為顯著增加階段。由于改革開(kāi)放以來(lái)人口數(shù)量的激增以及城鎮(zhèn)化水平的不斷提高，壩上地區(qū)總建設(shè)用地面積由1980年初374.84 km2持續(xù)上升到2018年的540.90 km2，近40 a增加了44.30%。多年來(lái)土地荒漠化、土地沙化等一直侵?jǐn)_著壩上地區(qū)土地生態(tài)安全，荒草地、鹽堿地、沙地等未利用地在該區(qū)面積多達(dá)900～1 000 km2，但是情況好在多年來(lái)呈現(xiàn)逐步下降趨勢(shì)。根據(jù)模型結(jié)果，預(yù)計(jì)在2026年，壩上地區(qū)耕地、建設(shè)用地、未利用地均有較小上升趨勢(shì)，未來(lái)時(shí)期除林地將有所下降外，其余地類(lèi)變化不會(huì)太明顯。

圖1 壩上地區(qū)1980-2026年土地利用動(dòng)態(tài)度

表2 壩上地區(qū)1980-2026年土地利用面積變化

1980—2018年期間，壩上地區(qū)耕地主要由草地轉(zhuǎn)入而來(lái)，多達(dá)350.55 km2，這與該區(qū)多年來(lái)毀草種糧的農(nóng)業(yè)耕作方式密不可分。此外分別約有100 km2的建設(shè)用地和未利用地被用于開(kāi)墾種植。同時(shí)期內(nèi)，分別約有448.32和204.17 km2的耕地轉(zhuǎn)出林草地和建設(shè)用地，說(shuō)明該區(qū)近40 a來(lái)退耕還林還草政策取得實(shí)質(zhì)性的成果同時(shí)也以犧牲耕地的代價(jià)加速著城鎮(zhèn)化水平的進(jìn)程。多年來(lái)林地的擴(kuò)張主要由草地和耕地轉(zhuǎn)入而來(lái)，在2018年林地覆蓋度接近19%，說(shuō)明對(duì)該地區(qū)實(shí)行的綠化工程以及三北防護(hù)林工程已初見(jiàn)成效。壩上地區(qū)地表水域逐年萎縮，約62.48 km2的水域在近40 a間干涸繼而演變成荒廢地，其主要原因如上述所言與張北縣安固里淖的干涸有直接關(guān)聯(lián)。除水域外，草地的轉(zhuǎn)入也是未利用地增加的一大因素，壩上地區(qū)惡劣的氣候條件加速了土地荒漠化，水土流失、土壤風(fēng)蝕、地表水資源急劇下降，加上人為地對(duì)草地不合理墾殖后又廢棄等種種原因造成該區(qū)大量可開(kāi)發(fā)的土地演變成難以被利用的局面。預(yù)測(cè)在2018—2026年期間，耕地和草地的互相轉(zhuǎn)化將會(huì)保持活躍，轉(zhuǎn)移面積在200～250 km2不等。林地轉(zhuǎn)入為耕地和草地的面積均超過(guò)100 km2，建設(shè)用地的增長(zhǎng)預(yù)計(jì)主要由耕地轉(zhuǎn)入。由此觀之，該區(qū)的未來(lái)的土地利用變化主要發(fā)在在農(nóng)業(yè)用地之間互相的流轉(zhuǎn)，除此以外的其他地類(lèi)變化程度不會(huì)太大。總體來(lái)看，該區(qū)作為干旱半干旱區(qū)區(qū)，水域的萎縮表明該區(qū)生態(tài)的水源涵養(yǎng)功能在下降，干旱的氣候條件進(jìn)一步加劇，對(duì)該區(qū)生態(tài)環(huán)境造成一定影響。同時(shí)耕地和草地的互相轉(zhuǎn)換較為明顯，但草地面積減少嚴(yán)重，追求地區(qū)經(jīng)濟(jì)增長(zhǎng)而進(jìn)行草地墾殖會(huì)顯著減弱該區(qū)的水土保持、固碳等生態(tài)功能，更不利于該區(qū)荒漠化治理。同時(shí)，耕地大面積流入為林地，“退耕還林”等政策抑制該區(qū)生態(tài)環(huán)境惡化，將沙化嚴(yán)重或者生產(chǎn)能力低的耕地因地制宜來(lái)植樹(shù)造林，有效防患水土流失、洪澇干旱以及風(fēng)沙等自然災(zāi)害。從草地、水域明顯萎縮分析，盡管林地有所增加，但該區(qū)生境質(zhì)量總體是下降的，生態(tài)風(fēng)險(xiǎn)進(jìn)而趨于不穩(wěn)定(見(jiàn)表3)。

表3 壩上地區(qū)1980-2026年土地利用類(lèi)型面積轉(zhuǎn)移矩陣 km2

2.3 土地利用時(shí)空空間變化分析

通過(guò)ArcMap軟件的重心工具對(duì)1980—2026年6個(gè)時(shí)期的壩上地區(qū)6種土地利用類(lèi)型進(jìn)行重心分析得到結(jié)果如圖2所示。由圖2可以觀察到，6種地類(lèi)的重心在2000—2010年期間發(fā)生了顯著轉(zhuǎn)移變化。以2000年為界，耕地、草地、水域、建設(shè)用地4種類(lèi)型土地重心在該年以前轉(zhuǎn)移程度均較小。在2010年之后重心保持微小變動(dòng)。耕地、未利用地在每一時(shí)期變動(dòng)都較為明顯。從重心轉(zhuǎn)移方向上看，耕地和未利用地的重心轉(zhuǎn)移軌跡為復(fù)雜無(wú)序，但耕地最終是向東北偏移，預(yù)計(jì)未來(lái)會(huì)向東南方向發(fā)展；未利用地的重心最終落腳點(diǎn)轉(zhuǎn)向西南，接下來(lái)的8 a可能沿東北方向遷移。草地和水域的重心變化趨勢(shì)非常一致，在2000年以前基本保持相對(duì)穩(wěn)定，在2000—2010年10 a間兩者重心均向東北方向變遷后，隨之保持近似不變，預(yù)測(cè)結(jié)果顯示到2026年重心位置變化程度將不太明顯。林地重心的變化主要表現(xiàn)在2000—2018年向著西南方向轉(zhuǎn)變，兩階段內(nèi)轉(zhuǎn)變方向基本相同，前一時(shí)期變動(dòng)程度大于后一時(shí)期。建設(shè)用地的重心呈現(xiàn)兩極分化，2000年以前主要穩(wěn)定在西北方小幅度變動(dòng)，2000年以后遷向東南方，未來(lái)有向東北方向發(fā)展的趨勢(shì)。

圖2 壩上地區(qū)1980-2026年各土地利用類(lèi)型重心轉(zhuǎn)移

2.4 生態(tài)風(fēng)險(xiǎn)指數(shù)空間自相關(guān)分析

2.4.1 全局空間自相關(guān) 利用GeoDa軟件對(duì)1980—2026年6期壩上地區(qū)1 399個(gè)風(fēng)險(xiǎn)小區(qū)土地生態(tài)風(fēng)險(xiǎn)值的全局Moran’sI指數(shù)進(jìn)行計(jì)算，結(jié)果顯示，1980—2026年6個(gè)時(shí)期全局Moran’sI指數(shù)值分別為0.546，0.544，0.527，0.495，0.495，0.505。且均通過(guò)了p=0.05水平的顯著性檢驗(yàn)。6個(gè)時(shí)期Moran’sI指數(shù)值均保持在0.500左右，表明壩上地區(qū)風(fēng)險(xiǎn)指數(shù)值在空間上保持較強(qiáng)的正相關(guān)，即高風(fēng)險(xiǎn)區(qū)域其領(lǐng)域風(fēng)險(xiǎn)值也高，低風(fēng)險(xiǎn)區(qū)域亦然，呈現(xiàn)高度的相似性。在1980—2018年，全局Moran’sI指數(shù)持續(xù)下降，表明壩上地區(qū)生態(tài)風(fēng)險(xiǎn)值在空間上的趨同集聚性一直在減弱，預(yù)計(jì)到2026年，該趨勢(shì)會(huì)發(fā)生逆轉(zhuǎn)。

2.4.2 局部空間自相關(guān) 壩上地區(qū)1980—2026年6期的局部空間自相關(guān)LISA指數(shù)如圖3所示。由圖3可知，1980—2018年壩上地區(qū)的局部的生態(tài)風(fēng)險(xiǎn)空間集聚性階段性增強(qiáng)，“熱點(diǎn)”(高值集聚區(qū))和“冷點(diǎn)”(低值集聚區(qū))由初期相對(duì)分散的集中狀態(tài)逐步過(guò)渡到末期區(qū)域顯著集中狀態(tài)，“高—低”和“低—高”等局部“極值點(diǎn)”個(gè)數(shù)明顯少于“高—高”區(qū)(H—H)和“低—低”區(qū)(L—L)，分布較為分散，規(guī)律不明顯。生態(tài)風(fēng)險(xiǎn)指數(shù)“高—高”區(qū)主要分布在張家口壩上4縣城鎮(zhèn)區(qū)域和未干涸的安固里淖以及豐寧縣北部未利用地等區(qū)域。由于城鎮(zhèn)區(qū)域相比郊區(qū)人口較為集中，建設(shè)用地占比較高，經(jīng)濟(jì)水平更發(fā)達(dá)，對(duì)自然生態(tài)環(huán)境的破壞的風(fēng)險(xiǎn)性高；在2010年前后，位于張北縣的安固里淖隨著其在2000—2010年期間的完全干涸該湖淖地區(qū)由高風(fēng)險(xiǎn)區(qū)轉(zhuǎn)變?yōu)椴伙@著區(qū)，說(shuō)明水域具有較高生態(tài)風(fēng)險(xiǎn)性。生態(tài)風(fēng)險(xiǎn)“低—低”區(qū)較為顯著地分布在尚義縣南部、豐寧縣東南部和圍場(chǎng)縣東部，該范圍區(qū)域林草地覆蓋度高，建設(shè)用地相對(duì)較少，人類(lèi)活動(dòng)干擾作用小，因此生態(tài)風(fēng)險(xiǎn)程度低。“熱點(diǎn)”與“冷點(diǎn)”的分布與土地結(jié)構(gòu)的土地格局分布具有很強(qiáng)的關(guān)聯(lián)性。從風(fēng)險(xiǎn)區(qū)個(gè)數(shù)比例來(lái)看，“高—高”區(qū)和“低—低”區(qū)個(gè)數(shù)均有減小，且逐漸分散，表明壩上地區(qū)局部生態(tài)風(fēng)險(xiǎn)值兩極分化程度在減弱，局部風(fēng)險(xiǎn)急劇變化區(qū)域變小，表明該地區(qū)較研究初期風(fēng)險(xiǎn)地域性更加鮮明。由圖3可知，上述趨勢(shì)預(yù)計(jì)在未來(lái)2026年會(huì)小幅度地持續(xù)下去。

圖3 壩上地區(qū)1980-2026年生態(tài)風(fēng)險(xiǎn)值LISA指數(shù)分布

2.5 生態(tài)風(fēng)險(xiǎn)空間指數(shù)空間分布特征

近40 a來(lái)，隨著壩上地區(qū)人口增多、城鎮(zhèn)化進(jìn)程加快、氣候條件惡化，在綜合因素影響下，該區(qū)生態(tài)風(fēng)險(xiǎn)值最高值由1980年的0.34升至2018年的0.48，增加了41.18%，多年來(lái)區(qū)域生態(tài)風(fēng)險(xiǎn)等級(jí)有明顯上升，到2026年該值有可能上升至0.55。在ArcMap軟件中采用普通克里金(Ordinary Kriging)工具對(duì)壩上地區(qū)生態(tài)風(fēng)險(xiǎn)值進(jìn)行空間插值，其中克里格半變異函數(shù)采用的球面模型，步長(zhǎng)為4 km，得到1980—2026年3期該區(qū)生態(tài)風(fēng)險(xiǎn)圖，根據(jù)生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)結(jié)果采用等間隔分割法將本區(qū)內(nèi)生態(tài)風(fēng)險(xiǎn)等級(jí)分為5級(jí)，即：低風(fēng)險(xiǎn)區(qū)(ERI≤0.10)、較低風(fēng)險(xiǎn)區(qū)(0.100.40)，劃分結(jié)果如封3附圖5所示，經(jīng)過(guò)對(duì)各等級(jí)面積統(tǒng)計(jì)分析得到結(jié)果如表4。由封3附圖5觀察可知，生態(tài)風(fēng)險(xiǎn)值空間分布有著明顯的地域特征，主要由張北、尚義、康保、沽源、豐寧、圍場(chǎng)6縣城鎮(zhèn)區(qū)域?yàn)橹行娘L(fēng)險(xiǎn)程度逐漸降低。1980—2018年，壩上地區(qū)在研究初期無(wú)生態(tài)高風(fēng)險(xiǎn)區(qū)，研究末期出現(xiàn)的高風(fēng)險(xiǎn)區(qū)分布在張北縣中心城區(qū)。近40 a來(lái)中等風(fēng)險(xiǎn)區(qū)主要分布在壩上地區(qū)的城鎮(zhèn)、水域、未利用地等區(qū)域，具有明顯擴(kuò)張趨勢(shì)。豐寧縣和圍場(chǎng)縣以及尚義縣南部等植被覆蓋度較高的區(qū)域處于低風(fēng)險(xiǎn)或較低風(fēng)險(xiǎn)區(qū)，到2018年該區(qū)域臨近部分區(qū)域風(fēng)險(xiǎn)等級(jí)升高為中等風(fēng)險(xiǎn)區(qū)。預(yù)計(jì)在2018年，壩上地區(qū)較低風(fēng)險(xiǎn)和低風(fēng)險(xiǎn)區(qū)出現(xiàn)較明顯的萎縮趨勢(shì)，中等及中等以上風(fēng)險(xiǎn)區(qū)域呈持續(xù)擴(kuò)張趨勢(shì)，豐寧、圍場(chǎng)兩縣將分別出現(xiàn)不同程度面積的高風(fēng)險(xiǎn)區(qū)和較高風(fēng)險(xiǎn)區(qū)。

由表4可見(jiàn)，在1980—2018年，壩上地區(qū)低風(fēng)險(xiǎn)、較低風(fēng)險(xiǎn)區(qū)面積分別減少了720.52和637.77 km2，占地比例分別減少了3.75%和3.32%。中等風(fēng)險(xiǎn)區(qū)面積增加值最為顯著，由研究初期的906.15 km2至末期的1 834.49 km2，在原有基礎(chǔ)上增加了1.02倍多。較高風(fēng)險(xiǎn)區(qū)面積近40 a間擴(kuò)張了初期規(guī)模的近14.75倍，高風(fēng)險(xiǎn)區(qū)面積在研究期間由0增至100 km2多。預(yù)計(jì)在2026年，除低風(fēng)險(xiǎn)區(qū)面積有明顯減少外，其余風(fēng)險(xiǎn)等級(jí)土地均有不同程度的增加，主要表現(xiàn)在中等風(fēng)險(xiǎn)區(qū)擴(kuò)張。整體分析，壩上地區(qū)與其所屬的農(nóng)牧交錯(cuò)帶地區(qū)的土地結(jié)構(gòu)相似，均以耕地、林地、草地為主，其三者的變化決定了該區(qū)的景觀生態(tài)風(fēng)險(xiǎn)水平。該區(qū)域內(nèi)水土流失嚴(yán)重，干旱、荒漠化趨勢(shì)加劇，草地退化、墾殖造成了草地大面積銳減，從生態(tài)的角度加大了該區(qū)生態(tài)風(fēng)險(xiǎn)，使得區(qū)域生境質(zhì)量降低；從景觀指數(shù)的角度來(lái)看，壩上地區(qū)土地利用變化使得該區(qū)景觀破碎度加大，草地的脆弱度指數(shù)遠(yuǎn)低于其他大部分地類(lèi)，而其大量的轉(zhuǎn)換為脆弱度較高的未利用地和耕地，使得景觀風(fēng)險(xiǎn)值在綜合加權(quán)計(jì)算上升高，從而使得區(qū)域內(nèi)的景觀風(fēng)險(xiǎn)增高。

表4 壩上地區(qū)1980-2026年生態(tài)風(fēng)險(xiǎn)等級(jí)面積及比例

3 討論與結(jié)論

3.1 討 論

對(duì)比北方農(nóng)牧交錯(cuò)帶壩上地區(qū)的早期已有的研究分析，其中宋素青等[25]通過(guò)對(duì)張家口壩上4縣的景觀格局分析得到其研究區(qū)在退耕退林還草政策實(shí)施后耕地減少、林草地增加的結(jié)論，在本研究中2010—2018年時(shí)段的趨勢(shì)對(duì)其得到了很好的印證，但實(shí)際情況中該趨勢(shì)應(yīng)在2010年以前已發(fā)生；文獻(xiàn)[26]中所述壩上3縣2000—2015土地利用變化情況與本文2000—2010年時(shí)段趨勢(shì)基本一致。此外，本研究重點(diǎn)關(guān)注的是近40 a來(lái)整個(gè)時(shí)段初末期的壩上地區(qū)基于景觀生態(tài)角度的土地風(fēng)險(xiǎn)時(shí)空變化及預(yù)測(cè)，故未對(duì)更細(xì)化的時(shí)段進(jìn)行分析。

北方農(nóng)牧交錯(cuò)帶壩上地區(qū)作為生態(tài)脆弱區(qū)，多年來(lái)土地格局發(fā)生了較大改變，處于相對(duì)不穩(wěn)定態(tài)勢(shì)，對(duì)其未來(lái)情況進(jìn)行預(yù)測(cè)具有一定的必要性。一般而言，短期內(nèi)社會(huì)政策因素往往對(duì)區(qū)域內(nèi)土地的流轉(zhuǎn)起著關(guān)鍵作用。本研究選用的CA-Markov模型對(duì)土地的預(yù)測(cè)是基于以前一時(shí)段的土地轉(zhuǎn)移參數(shù)作為預(yù)測(cè)時(shí)期土地各地類(lèi)對(duì)應(yīng)像元的適宜性分布，因此僅考慮自然情景下的未來(lái)的土地生態(tài)風(fēng)險(xiǎn)變化的預(yù)估。未來(lái)會(huì)向著多情景模式下的預(yù)估進(jìn)行更綜合的評(píng)價(jià)。

區(qū)域的景觀生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)方法是基于土地利用變化的基礎(chǔ)上建立的，土地的變化繼而影響了景觀指數(shù)的變化，通過(guò)景觀指數(shù)構(gòu)建的景觀生態(tài)風(fēng)險(xiǎn)值也隨之改變。對(duì)于生態(tài)評(píng)價(jià)而言，該方法提供了一種便捷而高效的途徑，從高賓等[3]、張學(xué)斌等[4]等國(guó)內(nèi)學(xué)者廣泛地使用該評(píng)價(jià)模型可以認(rèn)為基于土地利用變化視角的景觀生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)體系是合理的、可適用的。但一般而言，生態(tài)評(píng)價(jià)中更多地應(yīng)考慮生態(tài)過(guò)程，從影響區(qū)域生態(tài)環(huán)境的各角度進(jìn)行綜合評(píng)判，本研究該不足之處在未來(lái)研究中會(huì)加以改進(jìn)?；谕恋乩米兓木坝^生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)對(duì)壩上地區(qū)甚至整個(gè)北方農(nóng)牧交錯(cuò)帶在生態(tài)評(píng)價(jià)中提供了的一種手段，對(duì)于本區(qū)內(nèi)的的整體生態(tài)建設(shè)、特定區(qū)域生態(tài)風(fēng)險(xiǎn)管理具有較為實(shí)用的參考價(jià)值。

3.2 結(jié) 論

(1) 1980—2018年期間，河北省壩上地區(qū)土地面積大小按其利用類(lèi)型排列為：耕地>草地>林地>未利用地>建設(shè)用地>水域，其中，耕地占據(jù)研究區(qū)面積近1/2，林草地面積共計(jì)占比41%左右。未利用地面積基本控制在在900～1 000 km2，其余水域兩者面積占比均為超過(guò)3%，預(yù)計(jì)在2026年，各類(lèi)土地利用變化仍然保持以往趨勢(shì)變化。

(2) 1980—2018年，壩上地區(qū)在最初20 a內(nèi)土地利用變化不夠明顯，但2000年以后，土地格局發(fā)生顯著變化，尤其表現(xiàn)在2000—2010年期間，土地更迭速度、各土地類(lèi)型重心轉(zhuǎn)移程度均處于最劇烈態(tài)勢(shì)。2000—2018年，耕地在前后階段內(nèi)發(fā)生明顯上升后下降的趨勢(shì)，林地持續(xù)增加，草地表現(xiàn)為明顯下降后不變的趨勢(shì)，說(shuō)明在該地區(qū)2000年后實(shí)施的退耕還林還草、三北防護(hù)林等相關(guān)政策取得了一定效果。水域在2000—2010年銳減的原因主要是華北地區(qū)第一大高原湖泊安固里淖該期間徹底干涸。未利用地多年來(lái)一直處于被開(kāi)發(fā)中，其面積處于平穩(wěn)下降趨勢(shì)。預(yù)計(jì)在2018—2026年，壩上地區(qū)土地之間的流轉(zhuǎn)與近40 a來(lái)的趨勢(shì)基本相同，主要發(fā)生在耕地、林地和草地之間互相的轉(zhuǎn)入轉(zhuǎn)出。

(3) 1980—2026年6個(gè)時(shí)期壩上地區(qū)生態(tài)風(fēng)險(xiǎn)值Moran’sI指數(shù)均在0.500左右，空間單元生態(tài)風(fēng)險(xiǎn)指數(shù)值與鄰域之間具有較強(qiáng)相關(guān)性?！盁狳c(diǎn)”與“冷點(diǎn)”分布與土地利用格局具有較高關(guān)聯(lián)性，其集聚性逐時(shí)期減弱，呈現(xiàn)“中心—分散”的趨勢(shì)。近40 a來(lái)壩上地區(qū)生態(tài)風(fēng)險(xiǎn)等級(jí)升高，局部生態(tài)風(fēng)險(xiǎn)區(qū)集中化，高風(fēng)險(xiǎn)、較高風(fēng)險(xiǎn)區(qū)域具有較強(qiáng)地域性，主要分布在張北、康保、尚義、沽源以及豐寧縣的城鎮(zhèn)、水域，未利用地區(qū)域。預(yù)計(jì)未來(lái)8 a壩上地區(qū)較高風(fēng)險(xiǎn)區(qū)域繼續(xù)保持?jǐn)U張趨勢(shì)，豐寧縣和圍場(chǎng)縣將分別出現(xiàn)小范圍的高風(fēng)險(xiǎn)區(qū)和較高風(fēng)險(xiǎn)區(qū)，研究區(qū)無(wú)論在目前階段還是未來(lái)時(shí)期，生態(tài)風(fēng)險(xiǎn)性均較大?？梢哉J(rèn)為，北農(nóng)農(nóng)牧交錯(cuò)帶壩上地區(qū)生態(tài)風(fēng)險(xiǎn)值分布與其土地格局具有很強(qiáng)的相關(guān)性，土地利用變化與土地結(jié)構(gòu)組成各異均會(huì)顯著影響區(qū)域內(nèi)的景觀生態(tài)風(fēng)險(xiǎn)水平。