基于CA-Markov和InVEST模型的土地利用格局變化對生境的影響研究——以北京淺山區為例

鄒天嬌

倪 暢

鄭 曦*

經濟發展和城市擴張伴隨著大規模土地利用方式、強度及格局的改變[1]。這會影響生境斑塊間的物質循環、能量流動和信息傳遞，使生境分布、面積和功能等都發生變化[2]。評價變化多以生境質量作為主要指標，即生態系統可以給物種的個體種群、群落及人類持續生存和繁殖所需條件的能力[3]。

近年來針對生境變化的研究，有學者選取過小城鎮[4]、河流[5]、自然保護區[6]等進行評估。隨著3S技術和生態模型發展，方法趨于空間化、定量化和精細化，使用較多的包括生境適宜性模型HIS[7]、IDRISI軟件中生物多樣性評價模塊[8]、SolVES模型[9]、InVEST模型[10]、ARIES[11]和MIMES[12]等。

其中InVEST模型被認為是目前最成熟、應用最多的生態系統服務評估模型[13]，它是自然資本項目(the Nature Capital Project)的成果，使用成本低、評估準確度高、空間分析能力強[14-15]。其中生境質量模型(Habitat Quality Model)更是廣泛應用于研究基于生境脅迫而評估生境質量[16]。Terrado等發現InVEST模型中，生境質量模塊與生物多樣性顯著相關[17]；Leh等用InVEST模型評估了非洲西部地區的生境質量[18]；包玉斌等分析了黃河濕地自然保護區的土地利用變化對生境的影響[19]；鞏杰等研究過甘肅白龍江流域生境質量時空分異[20]。

針對生境適應性管理，需要對未來的土地利用變化進行預測。目前最具有地理學代表意義的是以元胞自動機(Cellular Automata，CA)預測土地利用變化的模型[21-23]，GIS的發展更推動了其進步[24]。在國外，Batty等提出城市發展動態模型，并模擬了美國紐約州西部城市布法羅的城市擴張過程[25]；Nourqolipour等以馬來西亞油棕種植園為研究對象，利用模型模擬該地區2020年土地利用演變趨勢[26]等。我國關于CA模型研究熱潮始于20世紀90年代末期，后有孫賢斌等應用CA-Markov模型，研究撓力河流域不同時段土地利用對濕地的干擾[27]；趙建軍等運用3期土地利用數據，較好地預測出了向海保護區2010、2018年的土地覆被類型等[28]。

InVEST模型和CA-Markov模型雖然各自都有較久的使用歷史，但將二者結合的研究相對較少。本文擬結合InVEST模型和CA-Markov模型，先對1999和2017年北京淺山區生境質量、生境退化度等進行分析，并采用Logistic回歸模型進行景觀格局變化驅動力研究，模擬預測出淺山區2035年的生境情況，以期為相關部門進行區域生態規劃和出臺土地管理政策提供科學依據。

1 研究區概況與數據來源

1.1 研究區概況

北 京 位 于 華 北 平 原 西 北 部(115.7～117.4°E、39.4～41.6°N)，面積16 410.54km2。地勢西北高聳，東南低緩。本次所研究的淺山區位于平原城區與西北山脈的過渡地帶(圖1)，面積約為3 917.84km2。研究區主要以海拔和地形作為遴選依據，區內海拔高度100～300m，以丘陵和臺地地貌為主。氣候上，淺山地區沿西北山麓形成長條形暖溫帶，沿太行山、燕山山脈的迎風坡形成一條弧形多雨帶；風向受地形影響變化多樣。因位于生態交錯帶內，生態系統多樣化，物種種類和數量都勝于鄰近區域，并且還是地表水源涵養、地下水補充、農作物生產和防風固沙的重要區域。

1.2 數據來源

1.2.1 土地覆蓋類型分類

對1999、2017年的遙感影像進行處理，所用影像來源于地理空間數據云平臺(http://www.gscloud.cn/)，選擇云量少、月份接近的夏季數據，1999年采用Landsat7 ETM SLC-on影像，2017年采用Landsat8 OLI_TIRS影像，空間分辨率均為30m。

圖2 研究思路框架

在ENVI5.0中對這2期遙感影像進行預處理與影像拼接，使用研究范圍邊界裁剪，得到研究區遙感影像數據。根據土地資源用途、利用方式和覆蓋特征，參考《土地利用現狀分類》(GB/T 21010—2017)和相關研究文獻[29]，結合研究區實際情況，將土地分成5種類型，分別為水域、建設用地、林地、耕地和未利用地。采用監督分類與人工目視解譯結合的方法，對研究區1999和2017年2期景觀類別和格局信息進行解譯提取。

1.2.2 驅動因子

結合實際情況選取與景觀格局變化關聯較大的指標。驅動力因子選取高程、坡度、距主要道路距離、距居民點距離和距河流水庫距離5個因子。DEM數據來源于地理空間數據云平臺(http://www.gscloud.cn/)，經投影轉換、裁切得到研究區高程數據，并采用ArcGIS的Slope功能提取獲得坡度。交通道路、居民點和河流水庫分布圖均來自地圖，地圖數據是從國家基礎地理信息系統網站(http://nfgis.nsdi.gov.cn/)下載的1:400萬中國地圖和中國行政區劃圖，進行對應北京市環境圖層的提取。采用GIS中Buffer進行緩沖分析和賦值，獲得交通道路、居民點和河流水庫的區位緩沖圖。

1.2.3 威脅因子

參考前人對北京土地利用格局變化和對區域尺度生境質量的評估研究[29]，將城鎮用地、農村居民點、主要交通干道、耕地和礦區定義為生境的威脅源(表1)。城鎮用地、耕地的空間分布圖來源于對土地覆被類型分布圖重分類，交通道路和居民點的空間分布圖均來源于地圖，并確定了解譯出的5種土地利用類型對不同威脅源的相對敏感程度(表2)。

2 研究模型與研究方法

2.1 研究框架與技術路線

提取研究區1999、2017年的景觀類型信息，采用CA-Markov模型模擬淺山區2035年自然增長情景下的景觀格局情況；結合InVEST模型，評價2期景觀格局分布下的生境質量和變化，以及2035年模擬情景下的生境分布情況，總框架如圖2所示。

表1 威脅源的影響范圍及權重

表2 生境適宜度及其對不同威脅源的相對敏感程度

2.2 研究用到的模型

2.2.1 馬爾柯夫模型

馬爾柯夫(Markov)預測法來源于數學家馬爾柯夫的隨機過程理論，是一種特殊的隨機運動過程[30]。該模型通過計算各種不同狀態的初始概率和狀態見轉移概率關系，確定時間變化過程中的景觀格局演變趨勢，最終預測出景觀格局的變化。

景觀類型之間轉化無向性和雙向性強，類型和數量不斷隨機改變，難以用數學方式精確描述，符合了馬爾柯夫研究前提，關鍵點是確定好每種景觀之間轉化的概率。景觀格局變化過程中，每種狀態間轉化會形成一個轉移概率矩陣，公式如下：

式中，n為景觀類型數目；P x y為景觀x類型變為y類型的概率。Pxy需要同時滿足0≤Pxy≤1，且所有Pxy求和為1。

2.2.2 元胞自動機模型

元胞自動機(Cellular Automata，CA)模型，屬于不連續的時空動力學模型[31-33]。每個位于元胞空間的細胞，其特定狀態都是有限個的，會按所定義的局部規則同步更新。相互作用之下就產生了動態演化系統。CA模型的表達式如下：

式中，S為元胞有限而離散的狀態集合；t、t+1表示不同的時刻；N為元胞鄰域；f表示局部空間中元胞轉化規則。

2.2.3 CA模型與Markov模型的耦合

Markov模型和CA模型都有其局限性：Markov模型在預測土地利用變化時，只能預測到數量變化，而不能預測空間分布變化。CA模型主攻復雜系統空間演變，但沒有Markov模型長期預測優勢。二者結合則可以綜合考慮到各種因素影響和土地利用的歷史趨勢，從而更精確地模擬土地利用變化的時空演變過程。

在柵格圖中可以把一個個像元看作一個個元胞，元胞狀態表示景觀類型。在GIS中用轉換面積矩陣和條件概率圖像運算，探尋元胞狀態的轉移，據此模擬景觀格局變化[34]，步驟如下。

1)以疊置分析獲取景觀類型轉移概率矩陣、轉移面積矩陣和其他各種條件概率圖像(這些圖像來源于轉移概率矩陣，表示每個像元在下一時間點上被某一景觀類型覆蓋的概率)。

2)構造CA濾波器。分析鄰居離元胞距離的遠近，建立有顯著空間意義的權重因子，使其作用于元胞，確定元胞狀態改變情況。本研究采用5×5濾波器，即表示認定一個元胞周圍的5×5個元胞組成的矩形空間會對該元胞狀態改變產生顯著影響。

3)確定起始時刻和CA循環次數。以研究區域1999—2017年景觀類型轉移矩陣為基礎，將2017年作為預測的起始時刻，每年迭代一次，CA模型迭代次數定為18，在自然增長情景下模擬北京市淺山區2035年的景觀格局。

為保證模擬結果的準確，需要驗證模型，本研究中用Kappa指數檢驗了模型在模擬格局變化中的精度，公式如下：

式中，P0代表模擬正確的比例；Pc代表隨機情景下模型模擬正確的比例；Pp代表理想分類情景下模擬正確的比例。

再運用IDRISI軟件中的CROSSTAB模塊，輸入研究區2017年實際景觀類型圖與預測景觀類型圖，以對比檢驗模型精度。結果2017年Kappa系數結果為0.798 5，代表模擬效果良好，可使用驗證后的CA-Markov預測2035年的景觀類型。

2.2.4 InVEST模型

在InVEST模型中，Biodiversity模塊用棲息地質量的優劣來表征生物多樣性的彈性、恢復性、廣度和深度的變化。其計算結果生境質量和生境稀缺性能用來反映生物多樣性。

運行InVEST模型，按模型自帶參數在Biodiversity模塊里設置對應數據，模型會將景觀類型信息與生物多樣性威脅因子相結合計算，最終輸出生境質量指數和生境退化度圖層。

Habitatindex代表生境質量指數，公式如下：

式中，Qxj代表土地利用與土地覆蓋j中柵格x的生境質量；Hj代表土地利用與土地覆蓋j的生境適合性；k為半飽和常數，當k=0.5時，k值與D值相等；Dxj為景觀類型j柵格x的生境脅迫水平，公式如下：

柵格y中脅迫因子r對柵格x中生境的脅迫作用為irxy：

式中，dxy是柵格x與y之間的線性距離；drmax是威脅因子r的最大作用距離，是威脅因子的權重；βx是柵格x的可達性水平，1表示極容易達到；Sjr是景觀類型j對威脅因子r的敏感性，該值越接近1表示越敏感。

InVEST模型中假設認為，在一個生態系統中，地類對于威脅因子的敏感性程度越高，則該威脅因子對地類退化程度影響也就越大。生境退化程度由如下公式計算：

生境退化指數＝敏感性分布圖層×威脅強度分布圖層×權重值

3 結果與分析

3.1 景觀格局的時空變化

由圖3看出，1999—2017年，北京市淺山區的景觀空間格局主要變化表現在以下幾方面。1)在土地類型分布方面，農田和建設用地向淺山區內部擴張，在“鳳凰嶺-陽臺山-鷲峰-小西山”“天龍潭-白虎澗”和“泉水河-牤牛河”地區最為明顯，在密云水庫南部、平谷順義兩區山前地帶也多有碎片態建設用地出現。而多處林地退縮突出，整體而言景觀破碎化趨勢加強。2)在各類土地數量方面，建設用地和未利用地比例增長明顯，分別增長1.04%(40.55km2)和4.12%(161.61km2)，反映出淺山地帶不合理開發所導致的土地浪費現象較為嚴重。在5種研究的土地利用類型中，總量以耕地下降為最多，下降了140.65km2。消失的耕地大多轉變為了建設用地，也有少量轉為林地(表3、4)。

3.2 生境質量變化

生境質量高低代表的是生物多樣性高低，為便于判斷2個時間點的生境質量變化，根據輸出的生境質量取值，并用ArcMap中的自然斷裂法(Natural breaks)劃分等級。

由圖4看出，淺山區內生境質量較低的是農田和建設用地區域，生境質量較高的是林地和水域區域。生境質量變化與土地利用變化類型趨勢基本吻合，2017年較1999年，淺山區生境質量分布更碎片化，生境質量較差區域向質量較好區域更深楔入。還有許多整塊生境良好區域被低質量區域脈狀沖散，這些區域大多為溝域地帶。密云水庫北岸區域質量變差最為明顯。

3.3 生境退化度變化

生境退化程度分值的高低反映了該地類在當前保護程度下，受到人為威脅因子影響程度的大小，也就表示了其潛在的生境破壞與生境質量下降的可能性大小。平原與山區交界處的生境退化度十分突出，幾乎連成一線。生境退化度分值最高的是農田及建設用地區域附近(圖5)，表明其周邊生境具有較高潛在退化的可能性，1999年主要是密云、門頭溝區域高程在75～100m的地帶退化度較高。2017年較1999年，密云水庫、平谷區、房山區附近區域生境退化度明顯增強。

3.4 景觀生態格局預測與生境質量模擬

預測自然情景下淺山區2035年景觀格局(圖6)，發現2017—2035年的土地利用覆被變化基本會延續之前趨勢，農田和建設用地進一步向淺山腹地延展，侵占部分有林地和疏林地，但速度較之前有明顯減緩。之后預測分析了2035年淺山區的生境質量和生境退化度分布(圖7)。發現生境質量和退化度中等過渡地區減少，有兩極分化趨勢，且分布更趨于碎片化，將不利于地區完整生態過程的實現和可持續發展。

3.5 生境稀缺性分析

分析了1999—2017、2017—2035兩個時間段內的生境稀缺性(圖8)。生境稀缺性得分是通過對其空間位置和影響程度的計算，找出破碎化程度較高，且土地利用類型變化較為頻繁的地類斑塊或部分柵格，得分值的大小表示應受重視和保護的不同程度。淺山區密云水庫附近生境稀缺度較高，說明土地利用變化較頻繁且明顯，需要加強保護。其他得分值較高的區域大多分布在淺山近平原處的過渡地帶，并向淺山內部延伸，表明人為干擾在未來一段時間內仍將是影響淺山地帶生態環境的最主要因素，需要防微杜漸，分別因地制宜分析原因，制定保護性發展策略，實現土地精明利用和生態高效。

4 討論

4.1 研究方法局限性

本研究方法的探究實踐，尚存在一些局限性。

1)遙感數據主要來源于Landsat8 OLI夏季數據，遙感圖像時相差別及30m空間分辨率對提取結果有影響，導致景觀類型解譯存在一定誤差；解譯選取的數據準確度和數量有待提高。

2)基于CA-Markov模型預測景觀格局空間變化時，因為采用多地類模擬分析，各地類轉換概率的確定需要以幾期變化信息為主導，相關參數設定要參考地形、降雨和地區經濟發展政策等。故已知數據年份選取宜增多，使模擬參數設定更精確。

3)InVEST模型更多傾向于植被多樣性，認為生境質量與生物多樣性有很強的正相關關系，但現實中卻并非總是如此[35]。且InVEST模型中所有威脅因子均為疊加的，有時多種威脅的集體影響可能遠大于個體之和的威脅水平，但模型未能考慮。

4)數據獲取限制使得本研究的驅動因子指標體系有待完善。

4.2 研究結論與管理策略

4.2.1 研究結論

北京城市不斷擴張會對淺山景觀生態系統帶來一定的破壞，在追求社會和經濟發展的同時不能忽視保護生態環境。應限制盲目泛濫的建設，加強自然地區保育修復。研究發現1999—2017年，淺山區景觀空間趨于碎片化，一定程度上影響了生態涵養功能的實現。農田和建設用地向淺山內部擴張，逐漸與林地交織分布。當前山后區整體生境質量優于山前過渡地帶。預測2017—2035年仍將沿襲這一趨勢，但變化速度有所減緩。

就目前來看，最值得提起注意的是密云水庫北部廣大地區，在過去有大量林地被人工侵入，生境質量和退化度不容樂觀，而且稀缺性很高，即未來惡化可能性較大。還有一些線型的溝域地帶中，多處建設用地沿溝谷向山內延展，使原來完整的生境“龜裂”，主要是因為人類開發區楔入原自然背景，過去不合理的經營模式可能造成了人地不和諧因素的深化。但另一方面，研究也發現有一些區域有生境質量好轉的現象，推測應該屬于生態型保育規劃、退耕退工還林還草策略的成果。

4.2.2 對應管理策略

1)整體上：動態格局，持續精明發展。

從土地利用變化中提取規律，維護淺山生態安全和景觀連續完整性，可更多采用在山水基質上“嵌入斑塊”式發展，而非大面積地“沖撞、侵略”原自然生境。通過自然情景下的土地利用變化情況預測及生態效益評估，明確土地不同的規劃目標和開發強度閾值，對于開發土地集約高效利用。

2)當前生境質量不佳區域：分類究因，修復防護并進。

表3 1999—2017年北京市淺山區土地利用類型面積及比例

表4 1999—2017年土地利用轉移矩陣

對于當前生境質量表現為較差的區域，需對其深入調研，探究政策、開發和保護活動造成的影響，針對生物和非生物環境修正管理策略，控制不合理開發行為和不得當保護措施。探究其過去和未來的生境質量，如在過去或未來背景下有較好的情況，宜著重進行保育和修復，主要措施包括構建游憩廊道、綠道連接、對裸地增植補植、改造低質低效的林地為混合林等。對于過去和未來都情況不佳的地區，工作重點在于如何避免其不良態勢的擴大，主要措施以建設防護隔離帶等為主。

3)退化或稀缺度較高區域：規避風險，保障核心安全。

區內存在一些當前生境質量尚可，但有著較高退化度或稀缺度的區域，說明這些地區受到人為威脅因子影響程度大或土地變化頻繁，未來惡化風險也較大，應該嚴加防范，盡力避免“侵蝕”擴大化。具體措施包括杜絕威脅區域安全的行為、設生物圍欄，遠期升級為核心區，嚴格保護原有生境和物種等。

注：文中圖片均由作者繪制。

圖8 1999—2017年(8-1)、2017—2035年(8-2)生境稀缺度分布

[20] 鞏杰，馬學成，張玲玲，等.基于InVEST模型的甘肅白龍江流域生境質量時空分異[J].水土保持研究，2018，25(3)：191-196.

[21] Deadman P, Brown R D, Gimblett H R. Modelling Rural Residential Settlement Patterns with Cellular Automata[J].Journal of Environmental Management, 1993, 37(2): 147-160.

[22] Batty M, Xie Y. From cells to cities[J].Environment and Planning B: Planning and Design, 1994, 21(7): S31-S48.

[23] White R, Engelen G. Cellular automata as the basis of integrated dynamic regional modelling[J].Environment and Planning B: Planning and Design, 1997, 24(2): 235-246.

[24] Lett C, Silber C, Dube P, et al. Forest dynamics: A spatial gap model simulated on a cellular automata network[J].Canadian Journal of Remote Sensing, 1999, 25(4): 403-411.

[25] Batty M, Xie Y, Sun Z. Modeling urban dynamics through GIS-based cellular automata[J].Computers Environment & Urban Systems, 1999, 23: 205-233.

[26] N o u r q o l i p o u r R, A b d u l R a s h i d B M S, Balasundram S K, et al. A GIS-based model to analyze the spatial and temporal development of oil palm land use in Kuala Langat District, Malaysia[J].Environmental Earth Sciences, 2015, 73: 1687-1700.

[27] 孫賢斌，劉紅玉，李玉鳳，等.基于CA-Markov模型土地利用對景觀格局影響辨識[J].生態與農村環境學報，2009，25(1)：1-7；31.

[28] 趙建軍，張洪巖，喬志和，等.基于CA-Markov模型的向海濕地土地覆被變化動態模擬研究[J].自然資源學報，2009，24(12)：2178-2186.

[29] 陳妍，喬飛，江磊.基于InVEST模型的土地利用格局變化對區域尺度生境質量的影響研究：以北京為例[J].北京大學學報：自然科學版，2016，52(3)：553-562.

[30] Hulst R V. On the Dynamics of Vegetation: Markov Chains as Models of Succession[J].Vegetatio, 1979, 40(1): 3-14.

[31] 侯西勇，常斌，于信芳.基于CA-Markov的河西走廊土地利用變化研究[J].農業工程學報，2004(5)：286-291.

[32] 韓玲玲，何政偉，唐菊興，等.基于CA的城市增長與土地增值動態模擬方法探討[J].地理與地理信息科學，2003(2)：32-35.

[33] 藺卿，羅格平，陳曦.LUCC驅動力模型研究綜述[J].地理科學進展，2005，24(5)：79-87.

[34] 楊國清，吳志峰，祝國瑞.廣州地區土地利用景觀格局變化研究[J].農業工程學報，2006，22(5)：218-221.

[35] 吳健生，曹祺文，石淑芹，等.基于土地利用變化的京津冀生境質量時空演變[J].應用生態學報，2015，26(11)：3457-3466.