基于GIS-CA-Markov的白馬河流域土地利用變化模擬分析

荊延德1,李婷婷

(1.山東省高校南四湖濕地生態與環境保護重點實驗室,山東 濟寧 273165;2.曲阜師范大學 地理與旅游學院,山東 日照 276826)

1 引言

土地利用/植被變化(LUCC)對多種生態環境因素產生直接或間接的影響,進一步影響到區域和全球的可持續發展[1,2]。隨著社會經濟的快速發展,土地利用的變化速度和規模發生了巨大轉變,探討土地利用變化的時空變化規律,預測其未來變化速度與方向成為研究的熱點問題[3]。

目前進行土地利用變化研究的模型主要有遺傳神經網絡模型(GNN)[4]、CLUE模型[5]、CLUES-S模型[6]、系統動力學模型(SD)[7]、馬爾科夫鏈(Markov Process)[8]、元胞自動機模型(CA)[9,10]、CA-Markov模型[11,12]等。其中,元胞自動機—馬爾科夫(CA-Markov)模型利用Markov模型對土地利用變化進行預測并得出狀態轉移圖集,并利用元胞自動機模型對土地利用變化進行空間上的處理,這樣既結合了馬爾科夫鏈長期預測模擬土地類型數量結構的優勢,又具有元胞自動機模型對土地利用景觀格局的時空動態模擬功能,能很好地預測土地利用變化的發展趨勢[13]。目前學者們主要研究模型轉換規則的獲得[14,15],對時空尺度敏感性的研究很少考慮,這對土地利用模擬預測結果的準確性造成了很大的影響。研究表明,CA模型在預測過程中具有時空尺度敏感性,而只有在一定的時間尺度和空間尺度范圍內進行預測模擬才能得到與實際情況較符合的結果[16-20]。本文以已有的預測模擬研究中常用的像元大小和時間步長為候選條件,對其進行模擬精度的對比,選擇較合適的條件進行模擬預測,以增加結果的準確性。

白馬河發源于山東省鄒城市北部,是“南水北調”東線工程的重要水資源承載地之一,也是作為南四湖流域最重要的子流域之一,因此該區域的土地利用變化對南四湖的水資源調控和生態管理起著舉足輕重的作用。作為山東省區域經濟發展中的重要部分,該區域景觀格局變化對流經地經濟發展具有很大的影響[21,22]。本文以濟寧市白馬河流域為研究區域,在CA-Markov模型運行過程中通過改變研究區域的元胞空間尺度和時間尺度進行模擬預測,對比模擬精度,試圖探究元胞尺度、時間間隔等參數對CA-Markov模型模擬結果的影響機制,為白馬河流域的土地利用規劃和相關政策的制定提供理論依據。

2 研究區域與數據

2.1 研究區域

白馬河流域位于山東省西南部,是華北地區較大的淡水湖。白馬河全長60km,流經曲阜、兗州、鄒城、微山等縣(市)地區,于微山縣魯橋鎮九孔橋村入獨山湖。為了更好地研究土地利用變化機理,本文以白馬河及其所流經縣市區為研究區域進行研究，見圖1。

圖1 白馬河流域的地理位置

本研究區域地理位置位于116°41′—117°11′E、35°07′—35°33′N之間,氣候屬于暖溫帶大陸性季風氣候,四季分明、雨熱同期,降水多集中在夏秋兩季,年平均降水量約691.19mm,年平均氣溫約14.7℃;地勢較低平,大致由東向西傾斜,大支流有大沙河、望云河、石強河、十里溝和七里溝等。白馬河流域作為南四湖最重要的子流域之一,其土地利用/覆被變化研究可為流域周圍土地合理利用提供理論基礎,同時也對南四湖流域的水資源保護和生態管理具有重要意義。

2.2 數據來源

本研究需要的數據為分辨率30m的Landsat8 OLI-TIRS TM遙感影像和行政區劃圖。運用ENVI5.1和ArcGIS10.2軟件對圖像數據進行處理,再采用ENVI5.1軟件對遙感影像進行解譯,從而得到分辨率為30m×30m、100m×100m的1990年、2000年、2005年、2010年、2015年的五期土地利用類型圖。根據我國土地利用分類系統和白馬河流域實際土地利用情況對土地利用類型進行劃分,共分為耕地、林地、草地、水域、建設用地、未利用土地六大類。

3 研究方法

3.1 土地利用動態度

土地利用類型變化速率是反映區域土地利用變化劇烈程度的重要指標,一般可分為單一土地利用動態度和綜合土地利用動態度。動態度表示該土地利用類型的穩定情況[23],單一土地利用動態度模型的公式為:

式中,k表示研究時段內某一土地利用類型的動態度;At1、At2分別表示評測初期和末期各種土地利用類型的面積;t1、t2分別為變化初期和末期。

3.2 CA-Markov模型

CA-Markov模型既包括了元胞自動機模型模擬復雜空間變化的能力,也綜合了馬爾科夫模型長時間的預測優勢,可更精確地對土地利用類型變化進行模擬預測[30]。本研究中,像元可視為元胞自動機中的元胞,土地利用類型可視為元胞狀態,像元大小即為元胞的尺度大小,并利用濾波器對元胞鄰域進行定義[31]。通過Markov模型得到的兩期土地利用狀態轉移概率矩陣、轉移面積矩陣和狀態轉移圖集作為轉換規則進行預測下一時刻元胞狀態的變化趨勢,并對研究區域的土地利用/覆被變化進行模擬預測,具體步驟為:①數據格式轉化。本文在利用前期獲得的1990年、2000年、2005年、2010年、2015年五期土地利用分類圖的基礎上,為了方便將其導入IDRISI軟件中進行操作分析,轉換為軟件可用的數據格式(.rst的格式),包括矢量數據轉為柵格數據，柵格數據轉為ASCII文件，ASCII文件轉為IDRISI軟件識別的數據格式。②由于CA-Markov模型具有時空尺度敏感性,同時為了兼顧研究意義與效率,本文選取像元為30m和100m的兩種規格研究CA-Markov模型的空間尺度效應。CA濾波器是用來創建具有顯著空間意義的權重因子,作用于元胞以此來確定元胞下一時期的改變狀況,本文選用CA標準的5×5的濾波器。即確定中心元胞周圍為5×5的元胞矩形空間,迭代次數分別為5次和10次。③獲取轉移面積矩陣。將兩期土地利用類型圖導入IDRISI軟件中,利用Markov模塊得到該期間內土地利用轉移面積矩陣、轉移概率矩陣和狀態適宜圖集,其中轉移面積矩陣作為轉換規則參與下一步模擬與預測,將模擬得到的土地利用圖與實際圖進行一致性檢驗;以模擬精度較高時的元胞規格和迭代次數為條件對2020年的土地利用類型進行模擬預測。④模型預測的精度檢驗。為比較兩年土地利用圖的相似程度,利用Kappa系數[32]進行一致性檢驗,Kappa指數計算公式為:

式中,Pa為正確模擬的比例;Pc為隨機情況下期望的正確模擬比例;Pp為理想分類情況下的正確模擬比例,為100%。

本文將模型運行的迭代次數與其兩年間的時間步長設為一致。根據Kappa系數檢驗標準[33,34],當Kappa系數為0—0.20時,表明一致性極低;當Kappa系數為0.21—0.40時,表明一致性一般;當Kappa系數為0.41—0.60時,表明一致性中等;當Kappa系數為0.61—0.80時,表明具有高度的一致性;當Kappa系數為0.80—1.0時,表明兩圖像幾乎完全一樣。

4 模型檢驗與模擬分析

4.1 土地利用現狀分析

由ArcGIS10.2和ENVI5.1軟件分析所得到的五期土地利用類型圖,根據我國土地利用分類系統和白馬河流域實際土地利用情況進行土地利用類型的劃分,共分為六類土地利用類型,處理后的最終結果見圖2。從圖2可知,耕地面積廣布研究區域,草地和林地主要分布在研究區域的南部和東南部,東北部很少,而西部和北部幾乎沒有分布;建設用地面積主要占用其周圍的耕地而快速增加;未利用地在逐漸減少,水域面積變化較小。

圖2 白馬河流域五期土地利用類型

4.2 模擬精度檢驗

為了更加精確地預測研究區白馬河流域2020年的土地利用,首先針對該流域進行初步模擬精度檢驗分析,選擇一個較合適的時空尺度條件進行預測。本文對不同的時間步長分類進行了土地利用預測模擬。當時間步長即迭代次數設置為10次時,將1990年與2000年30m×30m的土地利用類型圖導入IDRISI軟件中,獲得1990—2000年期間內土地利用轉移面積矩陣、轉移概率矩陣和狀態適宜圖集,將轉移面積矩陣作為轉換規則參與2010年土地利用的模擬與預測,得到元胞為30m×30m的2010年土地利用預測圖。同理,將1990年與2000年100m×100m的土地利用類型圖導入IDRISI軟件進行相同操作,得到元胞為100m×100m的2010年土地利用預測圖。當時間步長即迭代次數設置為5次時,同上述步驟,依次將2005年與2010年元胞為30m×30m、100m×100m的土地利用類型圖導入IDRISI軟件,利用轉移面積矩陣分別得到2015年元胞為30m×30m、100m×100m的土地利用模擬預測圖(圖3)。首先，將2010年的土地利用模擬結果與實際2010年的土地利用數據在Database Query-CROSSTAB模塊中進行一致性檢驗,其中元胞為30m×30m、100m×100m所得到的Kappa系數值分別0.8876、0.8877;其次將2015年的模擬結果與實際2015年的土地利用數據進行一致性檢驗,其中元胞為30m×30m、100m×100m所得到的Kappa系數值分別為0.9121、0.9068。由此可以得出在元胞為30m×30m、迭代次數為5次的條件下對2015年進行模擬預測時得到的Kappa系數最高為0.9121。由Kappa系數檢驗標準可知,在此條件下CA-Markov模型預測模擬的2015年土地利用圖和實際幾乎完全一致,即此時空尺度為該研究區域預測條件中較適宜的模擬條件。

圖3 2015年、2010年白馬河流域土地利用預測

4.3 土地利用特征分析

由精度檢驗分析可知,分辨率為30m×30m、迭代次數為5次時預測精度達到最大,因此本文選擇在此條件下進行后續模擬分析,只對與該條件相關的幾期(2005年、2010年、2015年三期)土地利用進行特征分析。

土地利用面積變化:2005年、2010年、2015年三期的土地利用/覆被變化圖經過ArcGIS10.2統計分析,得到流域三個時期的土地利用類型面積與變化動態度,變化動態度即可計算各土地利用類型面積在某段時期內的變化速度(表1)。由表1可知,在2005—2015年研究期間,研究區域2005年、2010年、2015年不同土地利用類型的面積表現為耕地>建設用地>草地>水域>林地>未利用地,耕地是白馬河流域的主要土地利用類型(2005年、2010年、2015年分別占66.82%、66.87%、65.63%)。土地利用格局在該期間發生了巨大變化。整體看,耕地、林地、草地和未利用地均減少。其中,草地、林地減少較明顯。草地從2005年的446.57km2減少到2015年的328.38km2,減少幅度達26.48%;林地從2005年的114.83km2減少到2015年的92.31km2,減少幅度達19.67%。其次為未利用地和耕地,減少幅度分別為9.95%和1.93%。建設用地和水域逐漸擴張,建設用地呈現出十分強勁的上升趨勢,變化幅度在2005—2010年達28.2%,2010—2015年增加了6.91%;水域擴張不明顯,十年間僅增加了0.24%。

表1 三期土地利用面積

從整體來看(表1),各地類的土地變化速度相差較大。耕地在研究期間雖有增有減,但總體變化動態度較小;林地在整個研究期間呈減少趨勢,其中2005—2010年減少了19.52%,2010—2015年減少較緩慢;草地主要在2005—2010年減少且變化動態度較大,2010—2015年減少相對緩慢;水域整體增長較緩慢;建設用地急劇增長,在2005—2010年增加最快,動態度高達28.20%,2010—2015年增加相對緩慢;未利用地一直在減少,但變化動態度相對較慢。

土地利用結構的空間轉換分析:本文主要利用IDRISI軟件中的Markov模塊分別對2005—2010年和2010—2015年的土地利用類型圖進行轉移矩陣分析,得到轉移面積矩陣見表2、表3。由表2、表3可知,研究區域總體呈現出“兩增四減”的變化規律。即建設用地和水域呈現增加趨勢,耕地、林地、草地和未利用地呈現減少趨勢。水域變化不明顯,主要轉入地類為耕地、草地和建設用地,耕地轉入占優,主要轉出地類為耕地和建設用地。2005—2015年建設用地一直在增加,總體凈轉入576.34km2,轉入類型主要為耕地,占轉入的98%以上;轉出地類同樣主要為耕地,但轉出的耕地面積遠遠小于耕地的轉入面積。2005—2015年耕地一直在減少,轉出地類主要為建設用地和水域。其中,2005—2010年和2010—2015年耕地轉為建設用地,分別占總面積的25.23%和17.80%;轉入地類主要為草地和建設用地。林地轉換不平衡,2005—2010年凈轉出27.50km2,2010—2015年凈轉入0.07km2,轉出地類主要為耕地、草地和建設用地,其中主要為耕地和林地;轉入地類同樣主要為耕地和草地,尤其在2010—2015年草地轉林地面積為11.04km2。草地轉出地類為主要耕地,2005—2010年草地轉為耕地的面積占其總轉出面積的90%以上;草地轉入地類同樣主要為耕地。未利用地總體變化趨勢并不太明顯,但轉入轉出地類仍主要為耕地。

表2 2005—2010年土地利用轉移面積矩陣(km2)

表3 2010—2015年土地利用轉移面積矩陣(km2)

4.4 2020年土地利用模擬結果分析

本文選擇在元胞為30m×30m、迭代次數為5次的條件下，以2015年的土地利用數據為基礎對2020年的土地利用變化進行預測模擬。利用IDRISI軟件中的Markov模塊得到2010年和2015年兩期土地利用類型圖的轉移概率矩陣、轉移面積矩陣和狀態適宜性圖集。以2015年為基礎,以狀態適宜性圖集為轉換規則,迭代次數設置為5次(利用2010年和2015年預測2020年,時間間隔為五年),利用其轉移概率矩陣對2020年土地利用進行預測模擬。2015—2020年的土地利用轉化與2020年的土地利用類型見表4、表5和圖4。

表4 2015—2020年土地利用轉移面積矩陣(km2)

表5 2015—2020年土地利用類型變化與面積比

圖4 2020年土地利用變化狀況

由表4和表5可知,2015—2020年土地利用變化呈現出“三增三減”的變化規律,即耕地、林地、草地減少,水域、建設用地、未利用地面積增加。耕地變化幅度高達6.675%,面積比由65.63%減少到58.95%,主要轉為建設用地和水域,其中轉為建設用地的面積為488.95km2,占耕地總轉出的95.18%;轉入地類主要為草地和耕地周圍的林地。水域變化幅度較小,轉出地類主要為耕地、林地、建設用地和未利用地,其中林地占優;轉入地類主要為耕地和建設用地。草地面積比由9.18%減少到8.60%,轉出主要為耕地、林地和建設用地,其中轉為耕地的面積為48.58km2,占草地面積總轉出的77.59%;轉入地類主要為水域,變化幅度較小。建設用地變化幅度較大,面積比由18.26%增加到25.18%,轉出地類主要為水域;轉入地類主要為耕地、草地和水域,其中耕地轉入的面積占建設用地總轉出的97.96%,變化幅度和動態度都較大。林地和未利用地沒有明顯變化,變化幅度較小。

5 結論與討論

本文利用白馬河流域不同時期的遙感影像,以ArcGIS10.2、ENVI5.1和IDRISI17.00為計算平臺,結合CA-Markov模型,對研究區域2020年的地類進行了模擬預測。由于CA-Markov模型具有時空尺度敏感性,因此首先將該模型進行模擬所需的元胞規格和時間步長設為不同,分別對2010年和2015年地類進行預測;其次，將預測結果分別與實際年份的土地利用類型進行對比,得出模擬精度較高時的元胞規格與時間間隔;第三,在次較適宜的時空尺度基礎上,預測研究區域2020年的地類。研究分析表明:①以Markov模型中獲得的兩期土地利用狀態適宜圖集作為模型轉換規則具有可行性,CA-Markov模型在預測白馬河研究區域的土地利用類型時具有很好的適用性。②本文對比分析了像元大小為30m×30m,時間尺度5年;100m×100m,時間尺度5年;30m×30m,時間尺度10年;100m×100m,時間尺度10年共四種情況,分別在CA-Markov模型中進行模擬預測,并將預測結果與實際年份土地利用圖進行一致性檢驗。對2015年的土地利用/覆被變化進行模擬預測時得到的Kappa系數最高為91.21%,由此得出在本文研究中像元大小為30m×30m、時間尺度為5年為較適宜的研究模擬條件,模擬精度最高。③在研究期間內,耕地、建設用地和草地為研究區域的主要用地類型,三者合計占白馬河流域面積的90%以上?？傮w看,耕地、林地、草地和未利用地均呈減少趨勢,其中草地減少幅度最大,減少了31.06%;建設用地和水域呈增加趨勢,建設用地增加最多為47.10%。2005—2020年,白馬河流域建設用地面積呈現出急劇增加的趨勢,由2005年的476.66km2增加到2020年的901.09km2,主要是隨著經濟的發展城市化速度加快,城市人口增加,導致建設用地大幅度擴張,建設用地主要是占用周圍耕地面積而增加,耕地轉為建設用地占建設用地總轉入的94.30%。耕地面積雖然在2005—2010年有輕微增加,但總體減少了11.78%,除了被建設用地大量占用外,耕地同樣轉出為水域,占耕地總轉出的7.57%,耕地為了彌補巨大的減少面積而占用建設用地和周圍的草地。草地有減少趨勢,減少了31.06%,由2005年的446.57km2減少到2020年的307.88km2,草地主要被耕地占用。林地、水域、未利用地并無大幅度變化。

土地利用/覆被變化是一個極其復雜的變化過程,為了使預測結果更加準確,應綜合考慮自然因素和人為因素對地理空間系統的影響。在今后的研究中可嘗試建立不同的影響因子和不同決策者共同作用的影響因素,加入權重要素建立不同類型的土地利用/覆被變化模型。CA-Markov模型具有時空尺度敏感性,本文僅對比了四種不同時空尺度,選擇較適宜的條件對2020年的土地利用進行預測分析,對比分析的時空尺度條件較少。今后研究中可嘗試對時空尺度條件進行更細致的對比分析,以得出可靠性、準確性更高的研究結果。CA-Markov模型綜合了元胞自動機和馬爾科夫鏈各自的優勢,提高了模擬的精度,該模型的預測結果較好地反映了白馬河流域的土地利用變化趨勢,可在一定程度上對該流域的土地利用/覆被變化動態監測和可持續發展提供參考。