基于EnKF-MCRP模型的生態(tài)用地?cái)U(kuò)張模擬研究

于　強(qiáng)　岳德鵬　Di Yang　張啟斌　馬　歡　李宇彤

(1.北京林業(yè)大學(xué)精準(zhǔn)林業(yè)北京市重點(diǎn)實(shí)驗(yàn)室， 北京 100083； 2.佛羅里達(dá)大學(xué)地理系， 蓋恩斯維爾 32611；3.國土資源部土地整治中心， 北京 100035)

?

基于EnKF-MCRP模型的生態(tài)用地?cái)U(kuò)張模擬研究

于強(qiáng)1岳德鵬1Di Yang2張啟斌1馬歡1李宇彤3

(1.北京林業(yè)大學(xué)精準(zhǔn)林業(yè)北京市重點(diǎn)實(shí)驗(yàn)室， 北京 100083； 2.佛羅里達(dá)大學(xué)地理系， 蓋恩斯維爾 32611；3.國土資源部土地整治中心， 北京 100035)

以生態(tài)脆弱區(qū)典型縣域磴口縣為研究區(qū)，基于2002、2007、2012、2015年4期遙感影像解譯數(shù)據(jù)，將數(shù)據(jù)同化技術(shù)引入生態(tài)源地的變化模擬中，考慮生態(tài)障礙和生態(tài)阻力構(gòu)建EnKF-MCRP模型，進(jìn)行磴口縣生態(tài)用地的變化模擬。結(jié)果表明，引入數(shù)據(jù)同化技術(shù)的EnKF-CA/Markov模型的模擬總精度達(dá)到82.4%，數(shù)據(jù)同化能夠減少誤差的積累。根據(jù)擴(kuò)展能力磴口縣生態(tài)源地共分為5個(gè)等級(jí)，其中3、4、5級(jí)生態(tài)源地的空間布局形成東北-西南、西北-西南的沙漠化防護(hù)格局。引入生態(tài)源地變化的EnKF-MCRP模型的生態(tài)用地?cái)U(kuò)張模擬精度最高，生態(tài)用地面積與空間布局最接近實(shí)際情況，方差達(dá)到0.4。此研究可為當(dāng)前以及未來的生態(tài)用地規(guī)劃和管理提供科學(xué)根據(jù)。

生態(tài)用地； 數(shù)據(jù)同化； EnKF-MCRP模型； 磴口縣

引言

生態(tài)用地是具有重要的生態(tài)系統(tǒng)服務(wù)功能或生態(tài)敏感性較高、生態(tài)環(huán)境較為脆弱、對(duì)維護(hù)區(qū)域關(guān)鍵生態(tài)過程發(fā)揮重要作用的土地單元[1]。在西北荒漠綠洲過渡區(qū)，覆蓋區(qū)域的生態(tài)系統(tǒng)極其脆弱，生態(tài)用地作為干旱區(qū)防治和減緩?fù)恋鼗哪铀贁U(kuò)展的“緩沖劑”[2]，還具有維護(hù)綠洲生態(tài)系統(tǒng)健康和持續(xù)發(fā)展、提供生態(tài)服務(wù)等多重功能，因此荒漠綠洲生態(tài)脆弱區(qū)生態(tài)用地變化模擬已成為全世界學(xué)者們的研究熱點(diǎn)。生態(tài)用地的重建、模擬及預(yù)測可以為荒漠綠洲區(qū)生態(tài)用地資源的合理利用、科學(xué)規(guī)劃以及改善區(qū)域生態(tài)環(huán)境提供有效示范和支撐。

最小累計(jì)阻力模型(Minimal cumulative resistance，MCR)被廣泛應(yīng)用于物種保護(hù)和景觀格局分析等生態(tài)領(lǐng)域[3-4]，眾多學(xué)者利用MCR模型進(jìn)行了旅游用地規(guī)劃、城市用地?cái)U(kuò)張模擬、生態(tài)保護(hù)用地?cái)U(kuò)張趨勢分析、城市土地生態(tài)適宜性評(píng)價(jià)、山區(qū)城市增長邊界劃定等研究[5-6]，在構(gòu)建生態(tài)阻力模型時(shí)均考慮了源、距離和阻力3個(gè)因素，目前考慮源地等級(jí)研究較少，考慮源地的變化則基本沒有[7]。作為源地的生態(tài)用地意義重大，其隨著時(shí)間及外界條件不斷發(fā)生變化，數(shù)據(jù)同化能夠同化不同來源的觀測數(shù)據(jù)以修正動(dòng)態(tài)過程模型，目前眾多學(xué)者應(yīng)用不同的數(shù)據(jù)同化算法進(jìn)行了城市擴(kuò)張模擬[8]、土地利用模擬[9]、陸面系統(tǒng)模擬[10]等研究，研究表明數(shù)據(jù)同化能更新(修正)模型并提高模型的模擬精度[8]。本文將目前的研究熱點(diǎn)數(shù)據(jù)同化技術(shù)引進(jìn)到生態(tài)源地的變化模擬中[11]，利用同化值修正源地變化的模擬結(jié)果，并將源地的變化引入到生態(tài)阻力面模型中，另外考慮生態(tài)源地?cái)U(kuò)張能力因子修正MCR模型，構(gòu)建適合荒漠綠洲區(qū)的生態(tài)阻力模型。

本文選擇荒漠綠洲區(qū)典型城市磴口縣為研究區(qū)，采用數(shù)據(jù)同化技術(shù)中的集合卡爾曼濾波算法模擬修正生態(tài)源地的變化，構(gòu)建考慮源地變化、源地等級(jí)、距離、生態(tài)阻力/障礙4方面的生態(tài)阻力模型，對(duì)磴口縣的生態(tài)用地?cái)U(kuò)張進(jìn)行模擬。

1　材料與方法

1.1研究區(qū)概況

磴口縣地處中國西北部(東經(jīng)107°05′，北緯40°13′)，位于內(nèi)蒙古河套平原源頭，黃河中上游，背靠狼山山脈，西鄰烏蘭布和沙漠。磴口縣氣候干旱少雨，水資源較為短缺，土地沙漠化嚴(yán)重，土地退化嚴(yán)重，區(qū)域蒸發(fā)量大，導(dǎo)致土地鹽漬化程度深，境內(nèi)海拔高度1 030～2 046 m，整個(gè)地形除山區(qū)外，呈東南高西北低，逐漸傾斜。屬中溫帶大陸性季風(fēng)氣候，歷年平均風(fēng)速3.0 m/s，瞬間最大風(fēng)速28 m/s，多年平均降水143.9 mm，多年平均蒸發(fā)量2 327 mm，多年平均氣溫7.6℃，無霜期136 d。全縣有6個(gè)土類，10個(gè)亞類，31個(gè)土屬，258個(gè)土種。黃河流經(jīng)磴口縣52 km，年徑流量3.1×1010m3，共有水域面積約24.07 km2。河套地區(qū)地下水埋深0.5～3 m，沙區(qū)地下水埋深3～10 m，山前洪積扇地下水埋深3～30 m，相對(duì)豐富的地表水與地下水對(duì)磴口縣發(fā)展以及生態(tài)用地?cái)U(kuò)張?zhí)峁┝擞辛ΡＵ稀?/p>

1.2數(shù)據(jù)來源與處理

選取磴口縣夏季且少云的TM影像(2002、2007、2012年)和OLI影像(2015年)為研究對(duì)象，同時(shí)利用研究區(qū)2012年1∶5 萬植被分布圖、2015年1∶5 萬土地利用圖、空間分辨率30 m的DEM等作為輔助數(shù)據(jù)。利用ENVI 5.1軟件對(duì)影像進(jìn)行波段合成、圖像增強(qiáng)和幾何校正處理，選擇最大似然監(jiān)督分類法對(duì)4期遙感影像進(jìn)行目視解譯[12]，提取磴口縣的生態(tài)用地類型信息，所提取的磴口縣生態(tài)用地為基礎(chǔ)性生態(tài)用地，其在改善環(huán)境、維持區(qū)域生態(tài)平衡方面具有不可替代的重要作用，包括湖泊、灘地、水庫坑塘等水域和自然荒漠灌林地、人工綠地等林草地。使用ArcMap 10.2進(jìn)行細(xì)碎斑塊處理，運(yùn)用疊加分析工具進(jìn)行空間數(shù)據(jù)分析，最終在ArcInfo 10.2環(huán)境下完成拓?fù)浜透腻e(cuò)處理。

1.3技術(shù)路線

利用所構(gòu)建的EnKF-MCRP模型進(jìn)行磴口縣生態(tài)用地?cái)U(kuò)張模擬，總體技術(shù)路線如圖1所示。

圖1　總體技術(shù)路線Fig.1　Technical roadmap

1.4基于EnKF的源地變化模擬

磴口縣的生態(tài)用地類型包括人工型生態(tài)用地和自然型生態(tài)用地，其中能作為生態(tài)源地的生態(tài)用地類型為濕地型生態(tài)用地和荒漠植被型生態(tài)用地[2]。生態(tài)源地是隨著時(shí)間不斷變化的，且生態(tài)源地的變化會(huì)影響區(qū)域生態(tài)用地格局。通過引入集合卡爾曼濾波的數(shù)據(jù)同化方法，使用CA/Markov模型為數(shù)據(jù)同化模型算子模擬源地的變化，并將源地變化引入到生態(tài)用地?cái)U(kuò)張模擬中[13]。

1.4.1集合卡爾曼濾波

集合卡爾曼濾波(Ensemble Kalman filter, EnKF)將模型狀態(tài)預(yù)報(bào)看成近似隨機(jī)動(dòng)態(tài)預(yù)報(bào)，用一個(gè)狀態(tài)總體(設(shè)數(shù)目為N，即集合數(shù))代表隨機(jī)動(dòng)態(tài)預(yù)報(bào)中的概率密度函數(shù)[14]，通過向前積分計(jì)算下一時(shí)刻狀態(tài)總體的概率密度函數(shù)，并得到該時(shí)刻的統(tǒng)計(jì)特性(如均值與協(xié)方差)[15-16]。

集合卡爾曼濾波算法計(jì)算步驟如下：

(1)模型狀態(tài)集合初始化為

(1)

式中n——模型狀態(tài)變量的個(gè)數(shù)

k——時(shí)間刻度

(2)計(jì)算每個(gè)狀態(tài)變量在第k+1時(shí)刻的預(yù)報(bào)值

(2)

ωk——期望為0、方差為Wk的高斯白噪聲

Wk——模型的誤差方差矩陣

M()為模型算子，本文指源的演變CA/Markov模型，在源地變化研究中，將源地變化過程看作Markov過程[17]，如利用2002年和2007年的源地?cái)?shù)據(jù)預(yù)測2012年的源地空間分布，并與2012年源地實(shí)際數(shù)據(jù)進(jìn)行對(duì)比，即2007年的源地斑塊對(duì)應(yīng)于Markov過程中的可能狀態(tài)，它只與2002年的源地斑塊狀態(tài)相關(guān)，不同等級(jí)源地相互轉(zhuǎn)換的面積數(shù)量或比例即為狀態(tài)轉(zhuǎn)移概率[18-19]。

(3)計(jì)算k+1時(shí)刻的卡爾曼增益矩陣Kk+1為

(3)

其中

H——觀測算子

R——預(yù)測誤差協(xié)方差矩陣

(νk+1～N(0,R))

(4)

(5)

其中

(6)

式中Yk+1——k+1時(shí)刻的狀態(tài)變量的觀測值

νk+1——期望為0、方差為R的高斯白噪聲

(5)判斷是否到結(jié)束時(shí)刻，即目標(biāo)年份。如果未到結(jié)束時(shí)刻，則返回步驟(2)，否則結(jié)束。

1.4.2EnKF與CA/Markov模型實(shí)現(xiàn)

利用集合卡爾曼濾波數(shù)據(jù)同化算法融合高分辨率遙感影像源地的觀測數(shù)據(jù)得出同化結(jié)果，其輸入數(shù)據(jù)為非布爾型，而CA/Markov模型模擬元胞的結(jié)果狀態(tài)只有兩種，即源地與非源地[13-20]。故將整個(gè)研究區(qū)劃分為若干個(gè)正方格區(qū)域，每個(gè)方格里有若干個(gè)像元(元胞)，得出每一個(gè)方格內(nèi)的源地元胞擴(kuò)張強(qiáng)度。

(7)

式中con(sij=source)——區(qū)域內(nèi)源地元胞的數(shù)量

利用ArcGIS軟件的Fishnet模塊將全區(qū)分為513個(gè)規(guī)則方塊，從中選出66個(gè)觀測點(diǎn)和56個(gè)檢測點(diǎn)(圖2)。

圖2　觀測點(diǎn)與檢測點(diǎn)空間分布圖 Fig.2　Spatial distribution of observation points and detection points

誤差統(tǒng)計(jì)方式采用方差(VSE)，其計(jì)算式為

(8)

式中y——真實(shí)值

x——同化值或預(yù)測值，y和x維數(shù)一致

1.5MCRP模型

生態(tài)用地?cái)U(kuò)展可以看作是生態(tài)綠色空間對(duì)其他用地空間的競爭性控制過程，且這種控制和覆蓋必須通過克服阻力來實(shí)現(xiàn)[5，21-22]，這樣生態(tài)用地?cái)U(kuò)展就可以模擬為從源到匯克服阻力做功的水平過程[23]?；旧鷳B(tài)阻力模型，即MCR模型的基本公式為

(9)

式中VMCR——最小累積阻力面值

fmin——某土地單元對(duì)不同的生態(tài)用地取累積阻力最小值Dij——從生態(tài)用地j到土地單元i的空間距離

Ri——用地單元i對(duì)運(yùn)動(dòng)過程的阻力系數(shù)

考慮不同生態(tài)源地有不同的生態(tài)擴(kuò)張能力，即不同生態(tài)源地的影響力是不同的，引入不同等級(jí)生態(tài)源地的擴(kuò)張能力因子Pj，修正MCR模型，為簡化計(jì)算，按照生態(tài)源地斑塊的面積對(duì)生態(tài)源地進(jìn)行分級(jí)，并賦予擴(kuò)張能力值。

利用EnKF-CA/Markov模型預(yù)測2017年的生態(tài)源地，利用EnKF算法將預(yù)測的2017年的生態(tài)源地與2012年的現(xiàn)狀生態(tài)源地進(jìn)行數(shù)據(jù)同化，得到同化后的生態(tài)源地空間分布值，該值既包含2012年的生態(tài)源地信息，又包含生態(tài)源地的未來變化趨勢，即將源地的變化引入到生態(tài)阻力模型中。

修正后建立源地?cái)U(kuò)張能力約束下的生態(tài)阻力模型——MCRP模型(Minimal cumulative resistance power)考慮了生態(tài)用地?cái)U(kuò)張的4方面因素，即生態(tài)源地變化、生態(tài)源地等級(jí)、距離和基面阻力特征，修正后公式為

(10)

式中VMCRP——生態(tài)用地?cái)U(kuò)展最小生態(tài)累積阻力面值

Pj——生態(tài)用地j所屬等級(jí)的擴(kuò)張能力因子

生態(tài)用地的等級(jí)越高，擴(kuò)張能力越強(qiáng)。

2　結(jié)果與分析

2.1源地變化模擬與分級(jí)

根據(jù)2002、2007、2012年影像解譯數(shù)據(jù)結(jié)合磴口縣實(shí)際情況提取生態(tài)源地，在荒漠綠洲區(qū)，濕地的生態(tài)作用巨大，故將所有水域提取出來均作為源地，在磴口縣荒漠灌林地的防風(fēng)固沙作用顯著，生態(tài)作用十分重要，故將面積大于0.1 km2的荒漠灌林地提取出來作為源地。

利用EnKF以CA/Markov模型為數(shù)據(jù)同化算子模擬源地的時(shí)空演變，集合大小設(shè)為30，數(shù)據(jù)同化算子參數(shù)主要包括距離因子、地形因子、水文因子，如圖3所示。

圖3　各類適宜性因子評(píng)價(jià)圖Fig.3　Evaluation of various types of suitability factors

分別利用CA/Markov模型和EnKF-CA/Markov模型，基于2002、2007年的源地?cái)?shù)據(jù)，預(yù)測2012年的源地空間分布(圖4b、4c)，與2012年的生態(tài)用地實(shí)際情況進(jìn)行對(duì)比驗(yàn)證(圖4a)。

圖4　2012年生態(tài)源地實(shí)際與模型模擬空間分布圖Fig.4　Actual and model simulated spatial distribution map of ecological source in 2012

利用IDRISI軟件進(jìn)行逐像元對(duì)比分析得出EnKF-CA/Markov模型的模擬總精度為82.4%，模擬結(jié)果方差為0.6，Kappa系數(shù)為0.513。傳統(tǒng)CA/Markov模型的模擬總精度為65.4%，模擬結(jié)果方差為1.3，Kappa系數(shù)為0.342。在模擬過程中，利用數(shù)據(jù)同化技術(shù)修正的模擬結(jié)果在精度與方差上都有改善，本研究通過引入集合卡爾曼濾波，模型能夠減少誤差的積累，提高模擬精度。

圖5　生態(tài)源地及等級(jí)劃分空間分布圖Fig.5　Distribution map of ecological sources and classification of space

利用EnKF-CA/Markov模型，基于2007、2012年的源地預(yù)測2017年的源地空間分布情況，得到2017年源地空間分布如圖5a所示。將2017年的生態(tài)源地作為EnKF模型的觀測值，與2012年的實(shí)際生態(tài)源地進(jìn)行數(shù)據(jù)同化，得到2017、2012年的生態(tài)源地同化數(shù)據(jù)(圖5b)，其中既包括2012年的生態(tài)源地信息，也包含生態(tài)用地變化至2017年的變化趨勢。

進(jìn)而根據(jù)源地?cái)U(kuò)展能力，對(duì)源地進(jìn)行等級(jí)劃分，共分為5個(gè)等級(jí)(圖5c)，級(jí)別越高表明源地的擴(kuò)張能力因子越大，5級(jí)源地主要分布在磴口縣包勒浩特嘎查、巴音溫都爾嘎查、溫都爾毛道嘎查、烏蘭布和農(nóng)場耕地外圍，自磴口縣東北貫穿至磴口縣西南，與奈倫湖相連。4級(jí)生態(tài)源地主要分布在巴音博日那嘎查、哈騰套海農(nóng)場、納林套海農(nóng)場的耕地外圍。3級(jí)生態(tài)源地主要分布在磴口縣防沙林場、奈倫湖、包爾蓋農(nóng)場附近。3、4、5級(jí)生態(tài)源地的空間布局形成東北-西南、西北-西南的沙漠化防護(hù)格局，對(duì)于防治烏蘭布和沙漠的擴(kuò)張具有重要意義。

2.2生態(tài)阻力面構(gòu)建

由于基面特性差異，生態(tài)用地在擴(kuò)張過程中所受的阻力是不同的，阻力系數(shù)反映了生態(tài)用地控制和覆蓋其他用地類型難易程度。將基面阻力分為2類：生態(tài)障礙和生態(tài)阻力。生態(tài)障礙是生態(tài)用地?zé)o法擴(kuò)展到的地塊，如建設(shè)用地、基本農(nóng)田等，構(gòu)成生態(tài)用地?cái)U(kuò)張的剛性約束，阻力系數(shù)無窮大。生態(tài)阻力是對(duì)生態(tài)用地?cái)U(kuò)張有不同等級(jí)阻礙作用的地塊，考慮地塊表面的固有屬性和附加屬性構(gòu)建了生態(tài)阻力評(píng)價(jià)體系，將生態(tài)阻力分為5個(gè)等級(jí)，1級(jí)表示阻力最小，5級(jí)表示阻力最大。固有屬性包括用地類型和地質(zhì)地貌特征。附加屬性包括生態(tài)系統(tǒng)服務(wù)價(jià)值和NDVI-水體距離。所構(gòu)建的阻力評(píng)價(jià)體系如表1、2所示。

首先確定生態(tài)用地空間擴(kuò)張障礙，第①類是基本農(nóng)田，是保證磴口縣區(qū)域經(jīng)濟(jì)、社會(huì)穩(wěn)定的保障，面積為796.28 km2。第②類是建設(shè)用地，是生態(tài)用地?cái)U(kuò)張的剛性限制。生態(tài)障礙分布圖如圖6a所示。

進(jìn)而確定生態(tài)阻力面，根據(jù)生態(tài)阻力評(píng)價(jià)體系，使用MCRP模型算法，利用ArcGIS軟件中的cost-distance模塊分別生成5個(gè)等級(jí)源地生態(tài)阻力面，進(jìn)行疊加計(jì)算生成生態(tài)用地?cái)U(kuò)張生態(tài)阻力面，如圖6b所示，累計(jì)阻力值最小為0，位于生態(tài)用地內(nèi)部，累計(jì)阻力值最大為77 780.8，位于烏蘭布和沙漠內(nèi)部。由累計(jì)阻力面可以看出，生態(tài)用地外圍阻力均較大，形成明顯的累計(jì)阻力“山脊線”，這是由于在荒漠灌林地外圍是大面積的沙漠，累計(jì)阻力低的區(qū)域較為破碎，整體性差，嚴(yán)重影響景觀流的流動(dòng)。干旱區(qū)由于生態(tài)環(huán)境脆弱，荒漠化嚴(yán)重使得景觀生態(tài)流流動(dòng)不暢，固有屬性導(dǎo)致其景觀流阻力較大，而城市化的加劇導(dǎo)致了人工不透水表面的增加，進(jìn)一步加劇了景觀破碎，增加了景觀流動(dòng)阻力。

表1　生態(tài)阻力評(píng)價(jià)體系Tab.1　Evaluation system of ecological resistance

表2　NDVI與水體距離生態(tài)阻力評(píng)價(jià)體系Tab.2　Evaluation system of ecological resistance between NDVI and water distance

圖6　生態(tài)障礙、阻力面及2015年生態(tài)用地空間分布圖Fig.6　Ecological barrier, resistance surface and spatial distribution of ecological land in 2015

2.3生態(tài)用地?cái)U(kuò)張模擬

利用ArcGIS 10.2軟件，使用python腳本語言編寫程序?qū)崿F(xiàn)模型的模擬。以2012年為基期數(shù)據(jù)，設(shè)置時(shí)間閾值為3 a，分別使用2012年實(shí)際生態(tài)源地、經(jīng)過數(shù)據(jù)同化后的生態(tài)源地?cái)?shù)據(jù)，利用生成的阻力面使用MCR模型和MCRP模型分別進(jìn)行生態(tài)用地?cái)U(kuò)張模擬，即模擬得出至2015年時(shí)的生態(tài)用地情況(圖6c)，為便于驗(yàn)證，以2015年的生態(tài)用地面積為檢驗(yàn)標(biāo)準(zhǔn)，磴口縣2015年的濕地型生態(tài)用地共137.92 km2。荒漠植被型、天然綠洲型、城鎮(zhèn)綠洲型、夾荒地型生態(tài)用地總計(jì)為1 042.96 km2。

圖7　生態(tài)用地模型模擬空間分布圖Fig.7　Spatial distribution map of ecological land use through model simulation

經(jīng)過迭代計(jì)算，利用EnKF-MCRP模型模擬的2015年生態(tài)用地面積共1 260.95 km2(圖7a)， MCRP模型模擬的生態(tài)用地面積為983.57 km2(圖7b)，EnKF-MCR模型模擬的生態(tài)用地面積為1 631.35 km2(圖7c)。與2015年的實(shí)際生態(tài)用地面積相對(duì)比，利用數(shù)據(jù)同化技術(shù)加入多時(shí)期生態(tài)源地變化的MCRP模型精度最高，逐像元對(duì)比方差為0.4，且整體的生態(tài)用地空間布局接近實(shí)際情況，利用MCRP模型的模擬精度其次，方差為1.7?？紤]生態(tài)源地的變化后，MCRP模型模擬結(jié)果比2015年生態(tài)用地真實(shí)面積略大，不使用數(shù)據(jù)同化技術(shù)的MCRP模型模擬結(jié)果要比實(shí)際偏小，說明數(shù)據(jù)同化技術(shù)能夠提高模型模擬精度，將生態(tài)源地的變化作為模型因素十分必要。與EnKF-MCRP模型相比，EnKF-MCR模型的模擬結(jié)果精度較低，逐像元對(duì)比方差為3.1，模擬面積比真實(shí)值高450.37 km2，表明在構(gòu)建生態(tài)阻力面時(shí)一方面要考慮不同等級(jí)源地的生態(tài)阻力，另一方面要考慮生態(tài)障礙的剛性限制，由模擬結(jié)果圖可知，未考慮生態(tài)障礙的生態(tài)用地?cái)U(kuò)張中明顯侵占了基本農(nóng)田和建設(shè)用地。

對(duì)56個(gè)檢測點(diǎn)的用地類型與實(shí)際用地類型進(jìn)行對(duì)比驗(yàn)證，其中相一致的賦值為1，不一致的賦值為2，56個(gè)檢測點(diǎn)的用地類型一致性對(duì)比發(fā)現(xiàn)，EnKF-MCRP模型模擬結(jié)果與實(shí)際情況相符程度最高，MCRP模擬結(jié)果有11個(gè)點(diǎn)與實(shí)際不一致，EnKF-MCR模擬結(jié)果有19個(gè)點(diǎn)與實(shí)際不一致。

在56個(gè)檢測點(diǎn)中選出16個(gè)具有代表性的檢測點(diǎn)并對(duì)其模擬結(jié)果進(jìn)行重點(diǎn)分析(表3)，其中檢測點(diǎn)1、3、4、5、6、11、14、15、16處EnKF-MCRP、MCRP、EnKF-MCR模型的模擬結(jié)果均與實(shí)際相符合，經(jīng)過同化作用后，同化的結(jié)果較明顯，特別是檢測點(diǎn)2，2015年檢測點(diǎn)2處為生態(tài)用地，EnKF-MCRP、EnKF-MCR模型的模擬結(jié)果與實(shí)際相符，但是MCRP模型模擬結(jié)果與實(shí)際不符合，表明數(shù)據(jù)同化對(duì)于模擬精度的提高作用明顯，檢測點(diǎn)10、12位于生態(tài)用地的邊緣，結(jié)合數(shù)據(jù)同化技術(shù)的模型模擬結(jié)果均與實(shí)際相符合，未利用數(shù)據(jù)同化技術(shù)的MCRP模型模擬結(jié)果則精度較低，若迭代年份改變可能會(huì)導(dǎo)致模擬結(jié)果錯(cuò)誤，總體來看，整個(gè)同化過程中各點(diǎn)均能從觀測點(diǎn)中獲取正確的增益信息，因而能夠得出一個(gè)較好的誤差相關(guān)矩陣，作用到增益矩陣中，從而得出較好的同化值。所以數(shù)據(jù)同化能夠根據(jù)觀測點(diǎn)的信息得出比模擬值更接近真實(shí)值的同化值。故引入集合卡爾曼濾波的數(shù)據(jù)同化方法到模型中，能取得比較好的模擬結(jié)果。檢測點(diǎn)7、8、9處于生態(tài)障礙中，檢測點(diǎn)13位于生態(tài)障礙邊緣，EnKF-MCRP、MCRP模型的模擬結(jié)果均與實(shí)際相符，EnKF-MCR模型的模擬結(jié)果為生態(tài)用地，表明生態(tài)用地發(fā)生了過度擴(kuò)張，在構(gòu)建生態(tài)阻力面模型時(shí)考慮不同等級(jí)生態(tài)源地的生態(tài)障礙十分必要。

表3　部分檢測點(diǎn)用地類型Tab.3　Land type for some detection points

注：1為生態(tài)用地，2為生態(tài)障礙，3為沙漠，4為非生態(tài)用地。

3　結(jié)論

(1) 磴口縣生態(tài)源地包括濕地型生態(tài)用地和荒漠植被型生態(tài)用地，生態(tài)源地對(duì)于維持區(qū)域生態(tài)用地穩(wěn)定發(fā)展有重要作用。本文將數(shù)據(jù)同化技術(shù)應(yīng)用于生態(tài)源地的變化模擬中，將生態(tài)源地的變化趨勢量化，EnKF-CA/Markov模型的模擬總精度達(dá)到82.4%，通過引入集合卡爾曼濾波，能夠減少誤差的積累，提高模擬的精度。

(2) 根據(jù)生態(tài)源地的擴(kuò)展能力，將磴口縣生態(tài)源地劃分為5個(gè)等級(jí)，研究表明：磴口縣的3、4、5級(jí)生態(tài)源地的空間布局形成東北-西南、西北-西南的沙漠化防護(hù)格局，對(duì)于防治烏蘭布和沙漠的擴(kuò)張具有重要意義，對(duì)于維持區(qū)域生態(tài)環(huán)境穩(wěn)定具有重大作用。

(3) 在構(gòu)建生態(tài)阻力面時(shí)考慮生態(tài)障礙和生態(tài)阻力兩個(gè)方面，將生態(tài)障礙引入阻力面構(gòu)建中能夠防止生態(tài)用地的過度擴(kuò)張，考慮生態(tài)源地變化、生態(tài)源地等級(jí)、距離、基面阻力4方面構(gòu)建MCRP模型，累計(jì)阻力值最小為0，最大為77 780.8，生態(tài)障礙阻力為無窮大。結(jié)果表明引入數(shù)據(jù)同化技術(shù)的EnKF-MCRP模型的生態(tài)用地?cái)U(kuò)張模擬精度最高，生態(tài)用地面積與空間布局最接近實(shí)際情況，數(shù)據(jù)同化對(duì)于模擬精度的提高作用明顯。

1周銳，王新軍，蘇海龍，等.平頂山新區(qū)生態(tài)用地的識(shí)別與安全格局構(gòu)建[J].生態(tài)學(xué)報(bào)，2015，35(6)：2003-2012.

ZHOU Rui, WANG Xinjun, SU Hailong, et al. Identification and security pattern of ecological land in Pingdingshan newly developed area[J]. Acta Ecologica Sinica,2015,35(6):2003-2012.(in Chinese)

2張紅旗，王立新，賈寶全.西北干旱區(qū)生態(tài)用地概念及其功能分類研究[J].中國生態(tài)農(nóng)業(yè)學(xué)報(bào)，2004，12(2)：5-8.

ZHANG Hongqi, WANG Lixin, JIA Baoquan. A conception of ecological land use and its function classification in arid area in Northwest China[J]. Chinese Journal of Eco-Agriculture,2004,12(2):5-8.(in Chinese)

3岳德鵬，王計(jì)平，劉永兵，等.京郊西北地區(qū)農(nóng)地利用與景觀格局時(shí)空特征[J].農(nóng)業(yè)工程學(xué)報(bào)，2008，24(6)：89-95.

YUE Depeng, WANG Jiping, LIU Yongbing, et al. Spatial and temporal characteristics of agricultural land-use and landscape pattern changes in northwest of Beijing[J]. Transactions of the CSAE,2008,24(6):89-95.(in Chinese)

4梁曼，孫璇，楊振.土地景觀生態(tài)流空間通達(dá)性評(píng)價(jià)——以荊州市為例[J].河南科學(xué)，2014，32(8)：1654-1658.

LIANG Man, SUN Xuan, YANG Zhen. Evaluation on spatial accessibility of ecological circulation of landscape—a case of Jingzhou City[J]. Henan Science,2014,32(8):1654-1658.(in Chinese)

5葉玉瑤，蘇泳嫻，張虹鷗，等.生態(tài)阻力面模型構(gòu)建及其在城市擴(kuò)展模擬中的應(yīng)用[J].地理學(xué)報(bào)，2014，69(4)：485-496.

YE Yuyao, SU Yongxian, ZHANG Hong’ou, et al. Ecological resistance surface model and its application in urban expansion simulations[J]. Acta Geographica Sinica,2014,69(4):485-496.(in Chinese)

6YU Q, YUE D P, WANG J P, et al. The optimization of urban ecological infrastructure network based on the changes of county landscape patterns: a typical case study of ecological fragile zone located at Deng Kou (Inner Mongolia)[J]. Journal of Cleaner Production, doi: 10.1016/j.jclepro.2016.05.014.

7LI F, Ye Y P, SONG B W. Evaluation of urban suitable ecological land based on the minimum cumulative resistance model: a case study from Changzhou, China[J]. Ecological Modelling,2015,318:194-203.

8張亦漢，黎夏，劉小平，等. 基于數(shù)據(jù)同化的元胞自動(dòng)機(jī)[J].遙感學(xué)報(bào)，2011，15(3)：475-491.

ZHANG Yihan, LI Xia, LIU Xiaoping, et al. The CA model based on data assimilation[J]. Journal of Remote Sensing,2011,15(3):475-491.(in Chinese)

9桂預(yù)風(fēng)，紀(jì)元，楊慧. 基于粒子濾波的隨機(jī)元胞自動(dòng)機(jī)城市用地?cái)U(kuò)張模型研究[J].數(shù)學(xué)的實(shí)踐與認(rèn)識(shí)，2015，45(17)：222-227.

GUI Yufeng, JI Yuan, YANG Hui. The research of urban sprawl with the stochastic cellular automata model based on particle filter[J]. Mathematics in Practice and Theory,2015,45(17):222-227.(in Chinese)

10張生雷，陳良富，蘇林，等.一種基于數(shù)據(jù)同化優(yōu)化陸面過程模式參數(shù)化方案的方法[J].中國科學(xué)：地球科學(xué)，2015，45(12)：1907-1922.

ZHANG Shenglei, CHEN Liangfu, SU Lin, et al. A data assimilation-based method for optimizing parameterization schemes in a land surface process model[J]. Science China: Earth Sciences,2015,45(12):1907-1922.(in Chinese)

11黃健熙，李昕璐，劉帝佑，等. 順序同化不同時(shí)空分辨率LAI的冬小麥估產(chǎn)對(duì)比研究[J].農(nóng)業(yè)機(jī)械學(xué)報(bào)，2015，46(1)：240-248.

HUANG Jianxi, LI Xinlu, LIU Diyou,et al. Comparison of winter wheat yield estimation by sequential assimilation of three different spatial-temporal resolution remotely sensed LAI datasets[J]. Transactions of the Chinese Society for Agricultural Machinery,2015,46(1):240-248.(in Chinese)

12李保杰.礦區(qū)土地景觀格局演變及其生態(tài)效應(yīng)研究——以徐州市賈汪礦區(qū)為例[D].徐州：中國礦業(yè)大學(xué)，2014.

LI Baojie. Study on landscape pattern changes and ecological effects in mining area—A case study in Jiawang mining area [D]. Xuzhou: China University of Mining and Technology,2014.(in Chinese)

13張亦漢，喬紀(jì)綱，艾彬.基于集合卡爾曼濾波動(dòng)態(tài)優(yōu)化CA模型參數(shù)的方法[J].測繪學(xué)報(bào)，2013，42(1)：123-130.

ZHANG Yihan, QIAO Jigang, AI Bin. Parameter optimization for CA model using ensemble Kalman filter[J]. Acta Geodaetica et Cartographica Sinica,2013,42(1):123-130.(in Chinese)

14王文，劉永偉，寇小華，等.基于集合卡爾曼濾波和HYDRUS-1D模型的土壤剖面含水量同化試驗(yàn)[J].水利學(xué)報(bào)，2012，43(11)：1302-1311.

WANG Wen, LIU Yongwei, KOU Xiaohua,et al. EnKF and HYDRUS-1D based data assimilation experiments for improving soil moisture profile prediction[J]. Journal of Hydraulic Engineering,2012,43(11):1302-1311.(in Chinese)

15李淵，基于數(shù)據(jù)同化的太湖葉綠素濃度遙感估算[D].南京：南京師范大學(xué)，2014.

LI Yuan. Chlorophyll a concentration remote evaluation in Lake Taihu based on data assimilation[D]. Nanjing： Nanjing Normal University,2014.(in Chinese)

16高斌斌.集合Kalman濾波在T106L19中期數(shù)值預(yù)報(bào)譜模式中的應(yīng)用研究[D].南京：南京信息工程大學(xué)，2008.

GAO Binbin. Application research on the ensemble Kalman filter (EnKF) with a medium-range numerical weather prediction (NWP) spectral model at a T106L19 resolution[D]. Nanjing： Nanjing University of Information Science & Technology, 2008.(in Chinese)

17HALMY M W, GESSLER P E, HICKE J A, et al. Land use/land cover change detection and prediction in the north-western coastal desert of Egypt using Markov-CA[J]. Applied Geography,2015,63:101-112.

18張瀅，丁建麗.綠洲土地利用變化未來趨勢預(yù)測及其調(diào)控研究[J].干旱區(qū)資源與環(huán)境，2006，20(6)：29-35.

ZHANG Ying, DING Jianli. Study on forecasting and regulation of future trend of oasis land use[J]. Journal of Arid Land Resources and Environment,2006,20(6):29-35.(in Chinese)

19王友生，余新曉，賀康寧，等.基于CA-Markov模型的藉河流域土地利用變化動(dòng)態(tài)模擬[J].農(nóng)業(yè)工程學(xué)報(bào)，2011，27(12)：330-336.

WANG Yousheng, YU Xinxiao, HE Kangning, et al. Dynamic simulation of land use change in Ji He watershed based on CA-Markov model[J]. Transactions of the CSAE,2011,27(12):330-336.(in Chinese)

20趙冬玲，杜萌，楊建宇，等. 基于CA-Markov模型的土地利用演化模擬預(yù)測研究[J].農(nóng)業(yè)機(jī)械學(xué)報(bào)，2016，47(3)：278-285.

ZHAO Dongling, DU Meng, YANG Jianyu, et al. Simulation and forecast study of land use change based on CA-Markov model[J]. Transactions of the Chinese Society for Agricultural Machinery,2016,47(3):278-285.(in Chinese)

21尹發(fā)能.江漢平原四湖流域景觀生態(tài)規(guī)劃與流域生態(tài)管理研究[D].上海：華東師范大學(xué)，2008.

YIN Faneng. Study of landscape ecological Planning and drainage basin ecological management in the Four-Lake basin[D].Shanghai:East China Normal University,2008.(in Chinese)

22RAY N, LEHMANN A, JOLY P. Modeling spatial distribution of amphibian populations: a GIS approach based on habitat matrix permeability[J]. Biodiversity & Conservation,2002,11(12):2143-2165.

23HERZOG F, LAUSCH A, MULLER E, et al. Landscape metrics for assessment of landscape destruction and rehabilitation[J]. Environmental Management,2001,27(1):91-107.

Simulation on Ecological Land Use Expansion Based on EnKF-MCRP Model

Yu Qiang1Yue Depeng1Di Yang2Zhang Qibin1Ma Huan1Li Yutong3

(1.BeijingKeyLaboratoryofPrecisionForestry,BeijingForestryUniversity,Beijing100083,China2.DepartmentofGeographyUniversityofFlorida，GainesvilleFL32611,USA3.CenterforLandConsolidation,MinistryofLandandResources,Beijing100035,China)

In the ecologically vulnerable area which locates in northwest arid and semiarid areas, ecological land is the important guarantee to maintain the security and stability of regional ecological environment. Studying the simulation of ecological land use change is of great significance. Accordingly, based on the typical ecologically vulnerable area—Dengkou County, this paper simulates ecological land use evolution of Dengkou County, using the four remote sensing image interpretation data of 2002, 2007, 2012 and 2015. Taking ecological eco-resistance barriers and EnKF-MCRP model into account, an ecological resistance surface was built. The EnKF-MCRP model was constructed to simulate the evolution of ecological sources considering the ecological sources change, ecological source level, distance and base surface resistance characteristics. The results showed that the combination of data assimilation and EnKF-CA/Markov model made a total accuracy of 82.4%, by using ensemble Kalman filter, the model can accumulate less errors and improve the accuracy of simulation, i.e., data assimilation can reduce the accumulation of errors. According to the expandability, ecological sources of Dengkou County were divided into five grades, of which the spatial layout of 3, 4, 5 grades formed the northeast-southwest and northwest-southwest pattern of desertification prevention. The building of EnKF-MCRP which takes the ecological source evolution into consideration made the highest precision. And the area of ecological sources and spatial distribution were the closest to the reality, of which the variance met 0.4. Different levels of ecological sources and ecological barriers were used to modify the model, which can improve the accuracy of the simulation results. This study could provide a scientific basis for the current and future ecological land planning and management.

ecological land; data assimilation; EnKF-MCRP model; Dengkou County

10.6041/j.issn.1000-1298.2016.09.039

2016-06-16

2016-06-30

國家自然科學(xué)基金項(xiàng)目(41371189)和“十二五”國家科技支撐計(jì)劃項(xiàng)目(2012BAD16B00)

于強(qiáng)(1987—)，男，博士生，主要從事3S技術(shù)在生態(tài)環(huán)境中的應(yīng)用研究，E-mail： yuqiang@bjfu.edu.cn

岳德鵬(1963—)，男，教授，博士生導(dǎo)師，主要從事景觀生態(tài)學(xué)和土地評(píng)價(jià)研究，E-mail： yuedepeng@126.com

K903

A

1000-1298(2016)09-0285-09