基于土地利用變化的黃河流域景觀格局及生態風險評估

劉希朝，李效順，蔣冬梅

（1. 中國礦業大學中國資源型城市轉型發展與鄉村振興研究中心，徐州 221116；2. 自然資源部海岸帶開發與保護重點實驗室，南京 210095；3. 江蘇自然資源智庫中國礦業大學研究基地，徐州 221116）

0 引 言

流域是一種結構復雜的地理區域，隸屬于生態、經濟及社會多個系統，具有生產、生活、生態及文化等多項相互聯系的功能[1-2]，過度利用其中一種或幾種功能勢必對流域整體功能產生負面影響[3-4]。在城鎮化迅速發展的背景下，對流域的開發利用主要聚焦在其經濟功能的發展，而對生態功能關注較少。研究流域生態風險能夠為流域管理和風險防控提供決策依據，進而促進流域可持續高質量發展[5]。

生態風險評價是為生態風險管理服務的，美國環境保護署于1992年定義生態風險為評估由于一種或多種外界因素導致可能發生或正在發生的不利生態影響的過程[6-7]。早期的生態風險評價主要針對環境污染展開研究，評價規模多數是單一風險源和單一風險受體。近年來，有關生態風險評價的研究較多關注生態系統的整體影響以及生態風險的空間相關性，評價規模也擴展到區域尺度，如流域、城市群、經濟帶等[8-9]。隨著對土地利用和生態風險研究的深入，從土地利用視角評估生態風險逐漸成為研究的主流，主要評價方法有景觀指數法和相對風險模型（Relative Risk Model，RRM）[6]。景觀指數法是選擇與生態風險相關的景觀格局指數構建評價模型[10]，如張玉娟等[3,7,11-13]選取不同類型的景觀格局指數構建生態風險指數（Ecological Risk Index，ERI），分別評價省級區域、流域、海岸帶等區域的生態風險。奚世軍等[14-15]在計算ERI的基礎上，分析ERI的空間相關性，并研究其驅動因素。相對風險模型是一種區域復合壓力風險評價模型，應用于水域、陸地等多個生態系統的生態風險評價。如張天華等[4,16-18]運用RRM分別評價水系工程、流域和縣級區域的生態風險，并提出相應的風險防控對策。王麗萍等[19]構建了基于RRM思想的多級模糊綜合評價模型，得到生態風險模糊評價分值，為改善區域生態環境提供參考。綜上，國內外學者主要通過建立評價模型來分析生態風險，評價方法和體系比較成熟，但從多尺度分析區域生態風險時空特征的研究較少。

黃河流域生態環境脆弱，地形復雜，水資源短缺，水土流失嚴重[20-22]。雖然近年來實施退耕還林、坡改梯等修復工程，但其生態保護工作仍面臨挑戰。在2019年召開的黃河流域生態保護和高質量發展座談會上，習近平總書記提出既要合理開發利用資源，又要兼顧資源的可持續性，提升資源的循環利用[23]，自此，黃河流域生態保護正式上升為國家戰略。分析長時間序列黃河流域景觀格局并從多尺度診斷其生態風險，對流域高質量發展具有重要意義。基于此，本文采用黃河流域2000、2010和2018年的土地利用數據，運用Fragstats4.2計算景觀格局指數，構建生態風險評價模型，從網格和縣域不同尺度診斷生態風險，揭示生態風險的時空分異規律、轉移特征和地類分布，并分析生態風險的空間關聯性，通過劃分生態風險管控區，提出黃河流域生態保護分區管控建議，以期為黃河流域風險防控和生態保護提供定量參考，并為其高質量發展提供決策支持。

1 材料與方法

1.1 研究區概況

黃河是中國第二長河，全長約5 464 km，流經青海、四川、甘肅、寧夏、內蒙古、陜西、山西、河南及山東9個省（自治區），約390個縣（區）[20]，如圖1所示。流域面積約79萬km2，占國土面積的8.3%，橫跨青藏高原、內蒙古高原、黃土高原和黃淮海平原4個地貌單元，地勢西高東低。2018年黃河流域（九省區）GDP為2.39×105億元，占全國GDP總量的26.5%。

1.2 數據來源與預處理

為深入診斷黃河流域景觀及生態風險時序特征，本文將研究時段定為2000—2018年。數據均來源于中國科學院資源環境科學數據中心（http://www.resdc.cn），包括2000、2010、2018年全國土地利用遙感監測數據（柵格數據，分辨率為1 km）和中國行政區劃數據（矢量數據）。按照研究區域的行政邊界裁剪土地利用數據，得到黃河流域2000、2010、2018年的土地利用數據，運用ENVI5.3對其進行精度檢驗，三期影像的綜合精度達91%以上。根據本文的研究需要，在ArcGIS中對其進行重分類等處理，將土地利用類型劃分為耕地、林地、草地、水域、建設用地和未利用地6大類（如圖2）。

1.3 研究方法

1.3.1 景觀格局指數

景觀格局包括景觀單元的類型、數量和空間分布，用景觀格局指數表達。為避免指標選取過多帶來的重復性研究[24-28]，根據黃河流域土地利用情況，綜合考慮各景觀格局指數的生態學意義，本文從景觀水平和類型水平兩方面選取7個景觀格局指數（表1）分析黃河流域景觀要素的數量、形狀和空間分布等特征。景觀水平指數可以反映研究區整體的景觀格局特征，類型水平可以反映研究區內不同斑塊類型的結構特征。運用Fragstats4.2軟件處理經預處理后的土地利用數據（柵格數據，1 km×1 km），可得到各景觀指數的數值。

1.3.2 生態風險評價模型

從景觀生態學角度出發，為建立景觀結構和生態風險的關系[7,11]，選取能夠反映區域生態損失的景觀格局指數[25-27,29-30]：景觀干擾度指數和景觀脆弱度指數，構建生態風險評價模型，具體如下：

1）景觀干擾度指數（Di）

Di反映不同景觀生態系統受到外界干擾的程度，由景觀破碎度指數（Fi）、景觀分離度指數（Si）、景觀分維數（FDi）計算得出[7,11]。Fi是指景觀的空間分裂破碎程度，可以反映自然或人為干擾對景觀的影響，其值越大，表明景觀的穩定性越差，受到干擾的程度越大。Si是指景觀中不同斑塊的分離程度，其值越大，表明景觀越分散，景觀分布越復雜。FDi是用來測定斑塊形狀對內部斑塊生態過程影響的指標，其值越大，表明斑塊形狀越復雜。指數計算公式如下：

表1 景觀格局指數應用尺度、生態學意義及計算方法Table 1 Application scale, ecological significance and calculation method of landscape pattern index

式中a、b、c為權重，且a+b+c=1，根據其作用程度[25-27,29-30]，分別賦值為0.5、0.3、0.2；ni為第i種景觀的斑塊數；Ai為第i種景觀的面積，km2；A為所有景觀的總面積，km2；Pi為第i種景觀的周長，km。

2）景觀脆弱度指數（Vi）

Vi反映不同景觀生態系統對外部干擾的敏感性，Vi值越大，生態系統越不穩定，越容易受到損害[25-27,29-30]。結合研究區的特點，對6類景觀賦值（耕地-4，林地-2，草地-3，水域-5，建設用地-1，未利用地-6）[25-27]，并做歸一化處理。

3）生態風險評價模型（ERI）

根據景觀干擾度指數和景觀脆弱度指數構建ERI，將景觀空間結構轉化為空間生態風險。ERI用于描述評價單元內的生態損失程度，其值越大，生態風險越高[7,11,29-30]。計算公式如下：

式中ERIk為評價單元k的生態風險指數；n為景觀類型的數量；Aki為景觀生態風險評價單元k中第i類景觀的面積，km2；Ai為景觀生態風險評價單元k的面積，km2。

為充分體現生態風險的空間分布情況，本文從網格和縣域兩個尺度展開研究，兩尺度的研究結果可以進行對比，也可以互相印證，使生態風險評價結果更具有科學性。在網格尺度上，以黃河流域景觀斑塊平均面積的5倍（30 km×30 km）創建漁網[13-14]，將其劃分成1 078個風險單元；在縣域尺度上，黃河流域共包含391個縣區，每個縣區作為1個風險單元，如圖3所示。

1.3.3 空間分析法

采用Moran’s I指數（I）和LISA指數（Ii）分析黃河流域生態風險的空間相關性。Moran’s I統計量是一種應用非常廣泛的空間自相關統計量，可以反映空間相鄰或相近單元屬性值的相似程度，公式如下：

式中Yi、Yj為變量在相鄰配對空間單元的取值；ωij為空間權重矩陣；ˉ為屬性值的平均值。I的取值在[-1，1]之間，當I＞0時，表明研究單元的觀測值趨于空間聚集，空間正相關；當I＜0時，表明空間呈離散分布狀態，空間負相關；當I=0時，表明空間不相關。

LISA指數又稱局部Moran’s I指數，可以反映某區域與相鄰區域間的差異程度及顯著性，公式如下：

式中n′為樣本數量，即研究單元的數量，S2為統計量的方差。當Ii＞0時，表示一個觀測值高（低）的區域被一個高（低）觀測值的區域包圍，即“高-高”（“低-低”）聚集；當Ii＜0時，表示一個觀測值高（低）的區域被一個低（高）觀測值的區域包圍，即“高-低”（“低-高”）聚集；Ii=0時，表示觀測區域與相鄰區域無關聯，即不顯著。

2 結果與分析

2.1 黃河流域景觀格局分析

2.1.1 土地利用變化

根據數據預處理的結果，統計2000、2010和2018年黃河流域各地類的面積及變化情況。如表2所示，黃河流域土地利用類型以草地和耕地為主，草地占比最大。從時序特征來看，2000—2018年耕地和未利用地面積分別減少14 243 km2和9 410.3 km2，林地、草地、水域面積在2000—2010年間維持相對穩定，2010—2018年間有小幅度增加，這與國家自2003年開始實行的退耕還林（草）政策相關。退耕還林（草）在一定程度上調整了土地利用結構，改善了農業生產條件。2000—2018年建設用地面積大幅度增加，增加面積為12 179.8 km2，尤其在2010—2018年，面積增長了68.36%，這是城鎮化建設和經濟發展的共同需求。現階段，黃河流域仍處在快速城鎮化發展的時期，對建設用地的剛性需求強勁。

表2 2000—2018年黃河流域地類面積及變化情況Table 2 Area and changes of land use types in the Yellow River Basin from 2000 to 2018

2.1.2 整體景觀特征

由表3可以看出，2000—2018年，黃河流域整體斑塊數（NP）呈上升趨勢，增長幅度為5.83%，整體斑塊密度（PD）呈上升趨勢，增長幅度為6.12%。由此可知，黃河流域景觀破碎化日益嚴重，破碎化進程加快，景觀分布分散。最大斑塊指數（LPI）呈上升趨勢，增長幅度為7.14%，說明黃河流域的優勢斑塊類型控制景觀的作用增強。景觀形狀指數（LSI）呈上升趨勢，說明景觀格局的形狀變得更為復雜和不規則，形狀種類多樣化。香農多樣性指數（SHDI）呈上升趨勢，增長幅度為1.49%，說明景觀類型豐富并趨于多樣化。蔓延度指數（CONTAG）呈下降趨勢，下降幅度為3.80%，說明流域景觀聚集程度降低，景觀破碎化加劇，景觀的抗干擾能力下降。整體來看，黃河流域受人類活動影響嚴重，生態功能有所損害，景觀格局趨于復雜化、破碎化和分散化。

表3 黃河流域景觀水平的景觀格局指數計算結果Table 3 Calculation results of landscape pattern index of landscape level of the Yellow River Basin

2.1.3 分類型景觀特點

由表4可以看出，2000—2018年黃河流域各類型的斑塊數量（NP）和斑塊密度（PD）均有不同程度的增減，其中，耕地、林地、水域和建設用地的斑塊數量和斑塊密度持續增加，草地和未利用地減少，表明各地類的破碎化程度及變化趨勢不同。其中，林地、耕地破碎化程度較高，主要是因為建設用地的不斷擴張和交通的發展，大片的耕地、林地被占用、分割。各類型的景觀形狀指數（LSI）與斑塊數量（NP）的變化趨勢一致，但變化程度較小，其中，草地和耕地的形狀最復雜。耕地、林地和未利用地的分離度指數（SPLIT）呈上升趨勢，草地、水域和建設用地呈下降趨勢，且建設用地的下降幅度最大，下降速度最快，說明建設用地斑塊不斷趨于聚集，但水域和建設用地的指數值仍較大，說明其斑塊間的分離程度與其他類型相比仍較高。

2.2 黃河流域生態風險分異

2.2.1 生態風險時空演變

根據1.3.2中的生態風險評價模型，測算得到每個風險單元的生態風險值，在ArcGIS中將其賦值給各單元的中心點，運用半方差變異函數優化數據，通過克里金插值，得到2000—2018年黃河流域在網格尺度和縣域尺度的生態風險空間分布圖（圖4），可以看出黃河流域生態風險的時空分異特征十分明顯。

基于自然斷點法，將生態風險分為5級，對應5個等級的風險區：低風險區（ERI≤0.46）、較低風險區（0.46＜ERI≤0.50）、中等風險區（0.50＜ERI≤0.55）、較高風險區（0.55＜ERI≤0.60）、高風險區（ERI ＞0.60），如圖4所示。時間序列上，2000—2018年黃河流域生態風險值的最高值呈上升趨勢，高風險區所占比例逐漸增加，網格尺度生態風險最高值由0.62增長至0.74；縣域尺度生態風險最高值由0.62增長至0.73。空間分布上，2000年高風險區主要集中在下游的山東段、河南段，上游的四川段（網格尺度）；2010年，中游的陜西段南部（西安、咸陽、渭南、商洛），山西段南部（晉城、運城、長治）等地區也呈現高風險；2018年網格尺度上，內蒙古的鄂爾多斯市等地區也呈現高風險。這些地區生態用地分散，建設用地面積不斷擴張，造成生態環境壓力，導致生態風險持續升高。

統計得到各風險等級的面積占比，如表5所示。網格尺度上，2000年以低風險和較低風險區為主，2010年和2018年以中等風險和高風險區為主，整個研究期間，中等風險、較高風險和高風險區面積不斷增加，其余等級區面積減少，其中高風險區面積占比提高了22.87%，低風險區面積占比下降了32.03%。縣域尺度上，2000和2010年以較低風險區和中等風險為主，2018年以中等風險和較高風險區為主，整個研究期間，較高風險和高風險面積不斷增加，其余等級區面積減少，其中高風險區面積占比提高了22.22%，低風險區面積占比下降了24.55%。

表5 生態風險各等級面積占比Table 5 Area proportion of each ecological risk level %

兩種尺度對比可以看出，網格尺度各等級面積的變化幅度較大，且高風險區和低風險區面積占比均大于縣域尺度，說明網格尺度下的生態風險變化敏感，較易產生兩端值。

2.2.2 生態風險轉移特征

通過疊加2000年和2018年生態風險分布圖，得到2000—2018年黃河流域各生態風險等級轉移情況。由表 6可以看出，在網格尺度上，面積轉移前3名依次為：低風險轉為中等風險（轉移面積119 437.00 km2），低風險轉為較高風險（轉移面積98 816.10 km2），較低風險轉為中等風險（轉移面積77 808.70 km2）。在縣域尺度上，面積轉移前3名依次為：中等風險轉為較高風險（轉移面積123 095.00 km2），較低風險轉為中等風險（轉移面積108 254.00 km2），中等風險轉為高風險（轉移面積101 085.00 km2）。由此可見，黃河流域的生態風險大都是由低等級向高等級轉移，生態風險加劇，進而導致生態問題嚴峻。

2.2.3 生態風險地類分布

運用ArcGIS統計工具得到2000—2018年黃河流域不同地類在各等級生態風險中的面積分布情況，如表7所示。2000年耕地主要分布在低風險區，占耕地總面積的47.97%，2010年和2018年主要分布在高風險區，占比分別為33.70%、37.79%，說明耕地所在區域的生態風險加劇，主要是人類活動干擾的結果，由于開發建設的需求，城鎮和農村周邊建設用地增加，導致耕地趨于破碎化，生態風險值增加。2000年林地主要分布在低風險區，占林地總面積的43.24%，2010年和2018年主要分布在高風險區，占比分別為40.30%、31.07%，說明林地所在區域的生態風險先增加后降低，生態風險的降低主要與國家實施的退耕還林及林地保護等政策相關。2000年草地主要分布在低風險區，占比43.24%，2010年主要分布在較低風險區，占比35.63%，2018年主要分布在中等風險區，占比37.00%，整體來看草地所在區域的生態風險增加，但分布在高風險等級區的面積比例較小。2000年和2010年水域主要分布在中等風險區，2018年分布在高風險區，占比46.18%，說明水域在2010—2018年生態風險劇增，相關水利設施的建設以及建設用地的開發對水域生態造成影響。2000—2018年，建設用地和未利用地的分布區域逐漸以高風險為主，這與建設用地規模增加，未利用地開發利用和存量土地更新相關。

表6 2000—2018年黃河流域生態風險轉移情況Table 6 Ecological risks transfer of the Yellow River Basin from 2000 to 2018 km2

表7 2000—2018年黃河流域各地類的生態風險面積Table 7 Ecological risk area of various types in the Yellow River Basin from 2000 to 2018

2.2.4 生態風險空間關聯

運用GeoDa軟件，根據黃河流域2000、2010、2018年的生態風險情況，分析其全局自相關性，結果如表8所示。在整個研究期間，黃河流域生態風險的Moran’s I大于0，呈現空間正相關，即生態風險在空間上相互影響，具有空間相似性。網格尺度上，三期的Moran’s I先上升后下降，縣域尺度上，三期的指數值持續下降，且網格尺度的Moran’s I大于縣域尺度，表明網格尺度生態風險的空間正相關性更強。總體來看，黃河流域生態風險的Moran’s I值在近20年間呈下降趨勢，這表明隨著社會經濟的發展和土地利用的變化，生態風險的空間聚集程度和空間分異性減弱。

同理，運用GeoDa軟件分析生態風險的局部自相關性。局部自相關可以反映生態風險的聚集類型和空間位置，如圖5 所示。網格尺度的“高-高”聚集單元（“熱點”區域）主要集中在下游的河南和山東，中游的山西南部和陜西南部，上游的四川等區域，2000—2010年“低-低”聚集單元（“冷點”區域）主要集中在上游的甘肅中部（白銀、蘭州、慶陽、定西）、中游的陜西北部、山西北部（朔州、沂州）、內蒙的鄂爾多斯、呼和浩特以及青海的玉樹、果洛等區域。2018年網格尺度的“冷點”區域減少，主要集中在陜西北部、山西北部以及內蒙的呼和浩特。縣域尺度2000年的“熱點”區域主要集中在下游的河南和山東以及上游的四川，2010年“熱點”區域集中在中游的山西、陜西南部以及下游的大部分區域，2018年“熱點”區域僅剩21個，集中在河南的洛陽、三門峽以及山東東營區域。縣域尺度的“冷點”區域與網格尺度基本一致。

表8 黃河流域生態風險的Moran’s I值Table 8 Moran’s I of ecological risk in the Yellow River Basin

整體來看，“熱點”區域和“冷點”區域的數量先增加后減少，總體上呈減少趨勢。空間聚集不顯著單元數量增加，說明生態風險局部空間聚集程度降低。

2.3 黃河流域生態風險管控區劃分

黃河流域生態風險的空間差異較大，本文綜合網格尺度和縣域尺度生態風險診斷結果（圖4），根據生態風險防控的需要，將兩尺度空間上風險等級不一致區域以較高的風險等級為主進行調整，且為方便區域管理，按照不跨越市級行政區的原則，將黃河流域劃分為生態風險重點管控區、嚴格管控區和一般管控區，如表9所示。

表9 黃河流域生態風險空間管控區Table 9 Ecological risk spatial control area of the Yellow River Basin

1）生態風險重點管控區：該區域生態用地分散，建設用地所占比例較大，且面積不斷擴張，造成生態環境壓力，生態風險也持續升高，應減少對建設用地的增量開發，轉向存量更新和集約高效利用，并全面保護林地、草地、水域等生態用地和高質量耕地。

2）生態風險嚴格管控區：該區域各類用地分布較為均衡，應發揮國土空間規劃的引領作用，建設用地的增量開發與存量更新同行，在促進經濟發展的同時，盡量減少對生態的負面影響，維持當前生態風險水平。

3）生態風險一般管控區：該區域主要以農牧業為主，生態風險較低。應合理開發布局建設用地，加強基礎設施建設，同時維持農業發展優勢，推進農業現代化。

3 討 論

2000、2010和2018年黃河流域（九省區）GDP分別為2.41×104、1.18×105和2.39×105億元，年均增長率14.43%，城鎮化率分別為32.11%、45.16%、56.20%，雖然城鎮化率不斷提升，但均低于全國平均水平（36.22%、49.68%、59.58%），說明黃河流域現階段仍處在城鎮化快速發展時期。與此同時，由于發展方式粗放、產業結構不合理以及建設用地無序擴張等問題，導致黃河流域的整體景觀格局破碎、生態風險加劇，生態環境壓力加重[31]。在黃河流域內部，區域發展不平衡現象突出，山東、河南等下游城市的經濟發展和城鎮化水平高于大部分中游、上游的城市，土地利用變化較為劇烈[20-21,31]，這也導致了生態風險的空間分布差異：高風險區主要集中在下游的山東、河南區域，低風險區主要集中在上游的寧夏、甘肅區域。總體而言，黃河流域的生態環境保護水平滯后于城鎮化發展，且快速城鎮化加劇了區域生態風險。

劃分生態風險評價單元是診斷區域生態風險的基礎。基于網格的劃分方式有利于生態風險空間異質性的表達，是現階段應用較廣的劃分方法。多數學者在研究中、小流域尺度或市、縣域尺度時（如喀斯特山區小流域、南四湖流域），都會通過劃分網格單元研究生態風險[3,14-15]。也有部分學者通過設置不同大小的網格單元，研究生態風險變化的尺度特征。基于行政區的劃分方式可輔助決策者制定更具有適宜性的風險防控及管理政策[32]，其應用較少且多應用于大尺度的研究。本文從網格尺度和縣域尺度分別診斷黃河流域生態風險的時空分異特征，且由上述分析結果可知，兩尺度生態風險的空間分布情況和變化趨勢相近，可以互相印證，提升生態風險評價結果的科學性。綜合兩尺度的結果，可以為流域高質量發展提供更為合理的建議。但兩尺度的結果又有不同，網格尺度更微觀，對生態風險變化的反映更敏感，可為縣域尺度或更高一級的行政區內部微觀生態風險管理提供參考，這也是論文需要進一步探索的方向。

4 結 論

1）土地利用結構分析表明，黃河流域的主要土地利用類型為草地和耕地。2000—2018年，耕地面積減少14 243 km2，未利用地面積減少9 410.3 km2，建設用地面積增加12 179.8 km2。現階段黃河流域仍處在快速城鎮化建設時期，對建設用地的剛性需求強勁。

2）景觀格局計算結果顯示，2000—2018年，黃河流域整體斑塊數、斑塊密度、最大斑塊指數和香濃多樣性指數的增長幅度分別為5.83%、6.12%、7.14%和1.49%，蔓延度指數的下降幅度為3.80%。黃河流域整體景觀格局趨于復雜化、破碎化和分散化。其中，林地、耕地破碎化程度最高，草地、耕地的景觀形狀和結構最復雜。

3）生態風險評價結果顯示，2000—2018年，黃河流域網格尺度和縣城尺度生態風險的最高值由0.62分別增長至0.74和0.73，網格和縣域尺度高風險占比分別提高了22.87%、22.22%。高風險主要集中在下游的山東段、河南段，上游的四川段，中游的陜西段南部和山西段南部。且生態風險大多由低等級向高等級轉移。

4）流域空間關聯分析看出，2000—2018年，黃河流域生態風險的Moran’s I均大于0，呈現空間正相關，但Moran’s I呈下降趨勢，說明生態風險的空間聚集程度和空間分異性減弱。網格尺度的Moran’s I大于縣域尺度，表明網格尺度生態風險的空間正相關性更強。