基于CA－Markov模型的海灣土地利用動態預測

吳季秋，俞花美，肖明，葛成軍，唐文浩*

1.熱帶作物種質資源保護與開發利用教育部重點實驗室，海南大學環境與植物保護學院，海南 海口 570228

2.中國科學院生態系統網絡觀測與模擬重點實驗室，北京 100101

3.海南大學經濟與管理學院，海南 海口 570228

土地作為人類生存和發展的基本物質條件，其資源配置的合理性已成為地區經濟的主要制約因素，關乎著土地資源的可持續發展與國家和民族的命運前途。2009年全國土地利用變更調查結果顯示，我國人均耕地面積只有1.38畝，還不到世界人均水平的1/3，其中九個省區人均耕地面積低于1畝，三個省區人均耕地面積低于0.5畝［1］。我國的土地資源本身內部稟賦不足，而土地退化、耕地銳減、生態環境質量惡化、用地結構不夠合理等現象進一步加劇了經濟發展訴求與土地資源稀缺間的矛盾。面臨嚴峻的現狀，開展土地利用趨勢變化分析，探求與生態環境協調發展相適宜的土地利用模式，一方面可為區域土地利用的可持續發展提供理論依據，另一方面為土地管理和資源優化配置提供支持和借鑒。運用GIS等先進手段進行分析，既是GIS與土地管理決策結合的理論嘗試，又是落實《國土資源部關于加快推進土地利用規劃數據庫建設的通知》(國土資發〔2011〕3號)的指示精神的具體體現。

元胞自動機-馬爾科夫(CA-Markov)模型綜合了元胞自動機(cellular automata，CA)和馬爾科夫模型各自的優點，即元胞自動機模型強大的模擬復雜系統空間演變功能和馬爾科夫模型長期預測的優勢，將二者有機結合在一起，有利于精確模擬土地利用變化的時空演變過程，是目前區域土地利用空間格局動態模擬預測研究中應用最多且較為有效的研究方法［2-3］。筆者以海南省文昌市八門灣地區為研究對象，采用馬爾科夫鏈分析研究區域1998年、2009年土地利用面積轉移矩陣及土地利用轉移概率矩陣，并根據馬爾科夫的結論，依據土地適宜性原則和環境友好型土地利用模式［4］，運用CA-Markov模型動態模擬預測該區域2020年的土地利用變化趨勢，進一步揭示不同的土地利用約束條件對于土地利用/土地覆蓋變化的影響，由此分析和探討土地利用與生態環境協調發展的理論途徑。

1 研究區域概況

研究區域八門灣位于海南省文昌市東南部(19°28'15″N ～ 19°43'44″N，110°42'35″E ～110°57'14″E)。八門灣因東邊的文教河與西邊文昌河在海灣出海口交匯，形成“八”字而得名，所轄灣內海域面積約680多hm2。流經研究區域的河流眾多，主要以文教河、文昌河及其各支流所形成的水系為主。研究范圍包括清瀾鎮、文城鎮、頭苑鎮、東閣鎮、文教鎮和東郊鎮，總面積約35377.73 hm2。

2 數據及預處理

2.1 數據來源

選取1988年、1998年、2009年三期TM影像，目視解譯標志采用海南省生態監測2009年、2010年兩期全國生態環境監測野外定點數據(海南站)建立目視解譯數據庫。數據處理按照《海洋生態調查指南》、《土地利用現狀調查技術規范》要求制定。相關參考資料有《文昌市土地利用總體規劃(2006—2020年)》、《文昌市產業布局圖》、《海南省環境公報(1998—2009)》、《海南統計年鑒》。

2.2 數據預處理

采用ERDAS IMAGINE 2011的raster模塊，在其中完成影像的波段組合、影像鑲嵌、大氣校正、去云等操作。利用數字化海南省1∶10萬地形圖對影像進行幾何校正，校正后的像元誤差控制在0.5個像元以內進行重采樣。土地覆蓋類型分類參照《土地利用現狀調查技術規范》和全國生態環境監測評價體系要求分為19個類型(表1)。

表1 研究區域土地覆蓋分類體系Table 1 Land covers classification system in study region

3 分析方法

CA-Markov的模擬過程大致分為兩個階段:先用Markov模型對于土地利用的時空推移概率進行計算，然后建立土地利用適宜性規則集，綜合前一階段的概率集在CA-Markov模型中推導出最終結果(圖1)。

圖1 研究路線Fig.1 Research route

3.1 Markov模型描述

目前，在眾多的土地利用分析模型、土地利用/土地覆蓋變化模擬和城市景觀格局變化預測研究中，Markov模型因其長期預測的明顯優勢而成為應用較為廣泛的一種土地利用變化建模方法［5-6］。Markov過程是基于馬爾科夫鏈的具有無后效性特征的一種狀態和時間都離散的隨機運動過程，即一個運動系統在未來某個t+1時刻的狀態僅與當前已知的t時刻所處的狀態有關，而與t時刻之前的狀態無關。馬爾科夫鏈分析的原理是通過一個運動系統的不同初始狀態及各種狀態的轉移概率來研究各種狀態之間的轉換情況，進而實現對未來狀態變化趨勢的預測。在馬爾科夫鏈中，系統狀態的轉移可用概率矩陣來描述，其數學表達式為［7］:

式中，P為某一土地利用類型從初始狀態到下一狀態的轉移概率矩陣;m和n為土地利用類型個數;Pij為土地利用類型i轉變為土地類型j的轉移概率。

自任何一個狀態開始，經過一次轉移，即會出現該系統所有狀態中的一種，該過程可表示為:

式中，Rt為初始狀態;Rt+1為下一狀態。因此，一次馬爾科夫模型可由已知t時刻的狀態推知t+1時刻的狀態。

在Idrisi的 Modeling菜單下點擊 Ecological environmental/Simulation Models，選擇 Markov模塊，選擇1998年、2009年八門灣區域的土地利用/土地覆蓋圖像，將這兩期圖像的間隔時間和預測時間周期均設置為11年，比例誤差設置為0.15，運算得出該區域1998年、2009年的土地利用面積轉移矩陣及土地利用轉移概率矩陣。

3.2 CA-Markov模型原理

元胞自動機(CA)，是時間、空間、狀態都離散，每個變量都只有有限多個狀態，且狀態改變的規則在空間的相互作用及時間上的關系皆表現為局部特征的網格動力學模型［8-12］。元胞自動機實際上是一個包含元胞(cell)及狀態信息(state)、元胞空間(lattice)、鄰居(neighbor)、規則(rule)及與時間(time)相關的轉換機制的函數［13］。

CA 模型可表示為［14］:

式中，S為元胞有限、離散的狀態集合;N為元胞的鄰域;f為局部空間的元胞狀態的轉化規則。

在Idrisi的Modeling菜單中調用CA-Markov模型。以2009年土地利用類型圖為基礎，預測時間間隔為11年，模擬2020年八門灣的土地利用類型轉化狀況。上述元胞自動機構成的六部分參數設置如下。

元胞:是元胞自動機的最小運算單元，按照一定的空間大小劃分元胞單位。文中元胞自動機的元胞即為圖像柵格，元胞大小設置為30 m×30 m。

狀態:每個元胞的狀態可以定義為多個變量。在空間分析中，其狀態是指元胞所具有的空間位置、大小、土地利用類型、周長等空間屬性。文中的元胞狀態即為其土地利用類型數字化屬性。

元胞空間:是指元胞在空間位置分布上網點的集合。在地理空間模擬的應用實踐中，元胞空間是具有空間位置屬性的概念，與空間柵格數緊密相關。

鄰居:元胞下一時刻的狀態是由元胞及其鄰居這一時刻的狀態所決定的。濾波器的工作原理是依據鄰居距離元胞的遠近位置關系來創建具有明顯空間意義的權重因子，并作用于元胞，進而促使元胞狀態的轉變，文中選用的是5×5濾波器。

規則:元胞自動機的函數運算準則是依據其設定的規則，將形成的約束函數用以指導元胞轉變的方向。在地理空間分析中，規則可由一系列的土地利用控制條件構成，按照經濟適用性、距離、高程和環境保護等適宜性實際需要，制定運行規則。Idrisi中依據多目標決策模塊(multi-criteria evaluation，MCE)確定的適宜性圖像集來確定［15］。以距離、坡度、降雨量為約束條件，將道路、水庫、居民點、坡度、降雨七個因子二值標準化為布爾圖像，即將適宜開發的區域賦值為1，反之賦值為0。

距離約束條件，從降低開發費用的角度出發，設置距離道路小于500 m區域適宜開發，賦值為1;大于或等于500 m的區域不適宜開發，賦值為0［16］。考慮到生產、生活產生的垃圾會導致地下水污染和周邊居民點的非點源污染，將距離水庫100 m緩沖區范圍設為非適宜開發，其余區域設置為適宜開發。

坡度約束條件，我國《水土保持工作條例》規定，25°以上的陡坡地，禁止開荒種植農作物。據此，采用數字化高程DEM數據，利用ArcGIS中的Slope工具得到坡度圖，再重分類為25°以上和以下兩個區域劃分禁耕區和適宜種植區。

降雨約束條件，采集了分布于區域內文昌河、文教河、清瀾紅樹林保護區的五個氣象監測站的降雨數據，統計了2003—2009年的降雨記錄，結果顯示，年均降雨量為913～1300 mm的西部農業生產總值明顯低于年均降雨量為1220～1935 mm的東部。對五個監測站點的降雨量進行克里格插值，并重新分類后得到降雨分布圖。年均降雨量在1200 mm以下的劃定為不適宜農業種植區域，年均降雨量大于1200 mm的劃定為農業生產種植適宜區。

時間:元胞自動機是一個在時間變化上離散的動力系統，即時間t是具有連續等間距特性的整數值。元胞在t+1時刻的狀態只取決于t時刻該元胞及其鄰近元胞的狀態。設置模擬預測圖像的時間間隔為11年。

4 結果與分析

4.1 Markov模型模擬結果分析

分析表2、表3可知，在1998年、2009年的土地利用變化中，養殖水面面積增加較大，其主要是由水田、果園和旱作耕地轉化而來，其中水田的轉移面積最大，占696.60 hm2，轉移概率為0.100，果園和旱作耕地轉移面積次之，分別為332.91和298.62 hm2，其轉移概率分別為0.099和0.093;水田面積的增加主要是由果園和養殖水面改變用途轉變而來，其轉移概率分別為0.340和0.264，總面積為1721.07 hm2，占該土地利用類型總轉化面積的35.54%。分析2009年土地利用類型面積減少的狀況，紅樹林減少的面積主要變成了水田和養殖水面，其轉移概率分別為0.262和0.269，共計267.66 hm2，占該土地類型轉化總面積的57.88%;減少的旱作耕地主要轉變為果園，沿海附近的耕地轉變用途變成了養殖水塘，轉移概率分別為0.156和0.093，共計799.47 hm2，占該土地利用類型轉移面積的28.98%。

通過對 Markov模型模擬結果的分析，對比1998年和2009年土地利用轉移矩陣，探討引起八門灣區域土地利用轉移特征的主要原因:1)八門灣海灣地區的養殖水面迅速增長，表明該區域沿海附近為追求經濟效益，大力發展高位池塘養殖，從2009年的遙感衛星影像上可以清晰看到沿文教河入海口至八門灣出海口沿岸一線布設建造了大量的圍欄養殖塘。這種局面雖然短期內帶來了一定的經濟回報，但卻帶來了一系列較為嚴重的生態問題，未經處理或處理不充分的養殖廢水的成倍增加給海洋生態環境造成了較嚴重的威脅和破壞，危及該地區的生態安全。2)八門灣的紅樹林地帶是海南島著名的紅樹林景觀之一，但是，近年來由于高位池養殖的方興未艾，導致大片的紅樹林生態系統轉變為海產品養殖池塘，同時受到經濟利益的驅使，不斷圍海造田，把沿岸紅樹林灘涂濕地改造成為農田，進而導致紅樹林面積的減少。

4.2 CA-Markov模型結果分析

通過對于馬爾科夫土地利用轉移矩陣的特征分析發現，當前的土地利用方式仍然存在許多不合理之處，1998—2009年的11年間，養殖水面和水田面積呈穩定增長的態勢，而紅樹林和旱作耕地面積卻明顯減少。這對當前八門灣生態環境的保護以及社會經濟環境的可持續發展構成了一定的威脅。基于此，在運用CA-Markov模型模擬預測該地區2020年的土地利用狀況時，相應地設置了一系列土地開發利用適宜性規則，以距離、坡度和降雨三方面的約束條件為出發點來制定元胞運行規則，其計算結果如圖2和表4所示。

圖2 CA-Markov模型預測2020年八門灣土地利用類型Fig.2 Land use type in 2020 in Bamen bay

預測結果表明:在土地利用適宜性規則形成的函數約束條件下，2020年防風林的面積增長迅速，較2009年增加了307.91 hm2，土地利用類型變化單一動態度指數達0.1585，天然林和紅樹林的面積也有所增加，動態度指數分別為0.0243和0.0184;伴隨防風林、天然林和紅樹林面積的增長，養殖水面、果園和水田的面積呈減少趨勢，養殖水面減少了1102.75 hm2，動態度指數為0.0322，果園和水田的面積分別減少了1088.18和2057.43 hm2，動態度指數分別達到0.0174和0.0172。由圖2可知，沿文昌河、文教河至八門灣海灣沿岸一線，紅樹林和防風林的覆蓋面積有所恢復和增長，相應的養殖水面有所縮減，陸域部分果園周邊的防風林和紅樹林植被覆蓋也有所增加。這一預測結果說明人為制定的一系列適宜土地開發利用和保護的經濟實用性及生態環境保護規則，不僅有利于保護該區域生態系統的平衡，維護生物多樣性及生態安全，而且可以保障該區域經濟社會與生態環境的協調可持續發展。

?

?

表4 2009年和2020年八門灣土地利用類型空間變化Table 4 Spatial distribution of land use type change between 2009 and 2020 in Bamen bay hm2

4.3 預測精度檢驗

CA-Markov模擬結果的精度驗證在目前還未形成較為統一的方法，通常使用的主要有整體逐個檢測與隨機抽樣檢測兩種。隨機抽樣檢測適用于驗證對象的分布較為均勻，代表性較強的情況。雖然該方法因為對象數據量少而處理快速，但其主觀隨意性大，容易因為樣本選取的典型性不夠而產生偏差。整體逐個檢測可以檢驗到參與模擬的每一個單元，可信度較高，但是計算冗繁，需要較長的處理時間。考慮文中研究區域數據量較小，故采用整體逐個檢測法。將模擬出的2009年預測結果〔圖3(b)〕與實際解譯的2009年土地利用類型結果〔圖3(a)〕進行求差柵格運算，前后一致的區域求差后土地利用類型屬性編碼值為0〔圖3(c)〕。最后以提取出土地利用類型編碼值為0的區域總柵格數，除以研究區域總柵格數，即得模擬精度。由表5可見，2009年預測的模擬精度為77.30%，反映出該模擬結果可信度較好，CA-Markov預測的土地利用類型具有較好的可靠性和適用性。

圖3 八門灣土地利用類型預測精度驗證Fig.3 Verification of forecast accuracy for land cover types in Bamen bay

表5 2009年預測精度對比Table 5 Comparison of forecast simulation accuracy in 2009

5 結論

(1)馬爾科夫鏈分析得出，研究區域養殖水面和水田的面積呈增長趨勢，相應地使得紅樹林和旱作耕地等土地利用類型的面積減少。

(2)基于以上的發展變化趨勢，從可持續發展的角度出發，根據土地利用適宜性原則，調整驅動因子空間分布規則，由此得到2020年的預測結果，防風林、天然林和紅樹林的面積將有所增加，動態度指數分別達0.1585、0.0243和0.0184;而養殖水面、果園和水田面積將有所減少，動態度指數分別為0.0322、0.0174 和0.0172。

(3)按照有利于生態環境改善，不影響當前經濟發展，降低經濟開發成本的原則，將距離、坡度和降雨三個因素對土地利用適宜性約束條件進行了調整，模擬結果顯示，研究區域的生態環境在該土地適宜性原則的約束下，天然植被有所恢復，人為的生產性土地擴張得到相應的控制。

6 不足與展望

對于CA-Markov的適宜性原則主要來自大量文獻的比較研究和土地利用政策綜合得出，由于缺少八門灣其他生態影響因子的實際研究成果和依據，因此只從距離、坡度和降雨三個方面的自然條件，同時僅結合政策或經濟約束條件進行模擬，未考慮更多其他可能的潛在影響因素，如人口數量的快速增長與產業布局的不斷調整都會使得區域的土地利用發生變化，從而對整個區域的生態格局產生各種影響。與此同時，隨著八門灣地區部分吸引旅游開發投資的產業政策的出臺，也會產生與之對應的土地利用空間格局的變化。因此，在今后的研究中，不僅要豐富其他潛在影響因子，如放牧、捕魚等人為影響因子的影響半徑測定，而且需要進一步探討按照不同的自然資源條件和經濟發展水平、政策導向影響下的土地利用空間變化及其與土地承載力和適宜性的評價機制，以完善土地利用變化趨勢研究的理論和手段。

［1］顏玉華:對我國18億畝耕地“紅線”的憂思［EB/OL］.中國改革論壇.(2012-01-21)［2012-03-01］.http://www.chinareform.org.cn/Economy/Agriculture/Forward/201201/t20120122_132973.htm.

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

［3］趙冠偉，陳穎彪，陳健飛，等.CA-Markov模型的空間尺度敏感性研究［J］.地理科學，2011，31(8):897.

［4］張秋玲，李東敏，邵亞杰.基于環境友好型模式的土地利用與生態環境協調研究:以河南省信陽市為例［J］.安徽農業科學，2009，37(23):11130-11132.

［5］ALEJANDRO F S，MIGUEL M R，OMAR R M.Assessing implications of land-use and land-cover change dynamics for conservation of a highly diverse tropicalrain forest［J］.Biological Conservation，2007，138(1/2):131-145.

［6］楊國清，劉耀林，吳志峰.基于CA-Markov模型的土地利用格局變化研究［J］.武漢大學學報:信息科學版，2007，32(5):414-418.

［7］唐啟義，馮明光.實用統計分析及其DPS數據處理系統［M］.北京:科學出版社，2000:95-105.

［8］朱海涌，李新琪，仲嘉亮.基于CA-Markov模型的艾比湖流域平原區景觀格局動態模擬預測［J］.干旱環境監測，2008，22(3):134-139.

［9］李書娟，曾輝，夏潔，等.景觀空間動態模型研究現狀和應重點解決的問題［J］.應用生態學報，2004，15(4):701-706.

［10］LI X.Neural-network-based cellular automata for simulating multiple land use changes using GIS［J］.International Journal of Geographical Information Science，2002，16(4):323-343.

［11］LOPEZ E，BOCCO G，MENDOZA M，et al.Predicting land-cover and land-use change in the urban fringe:a case in Morelia City，Mexico［J］.Landscape and Urban Planning，2001，55(4):271-285.

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

［13］HAN J，HAYASHIA Y，CAO X，etal.Application of an integrated system dynamics and cellular automata model for urban growth assessment:a case study of Shanghai，China［J］.Landscape and Urban Planning，2009，91(3):133-141.

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

［15］何丹，金鳳君，周璟.基于Logistic-CA-Markov的土地利用景觀格局變化:以京津冀都市圈為例［J］.地理科學，2011，31(8):903-909.

［16］毛蔣興，閆小培.城市交通干道對土地利用的廊道效應研究:以廣州大道為例［J］.地理與地理信息科學，2004，20(5):58-61.?