基于CA-Markov模型的惠州市生态与经济协调度时空演变特征及模拟预测

李 龙,吴大放①,刘艳艳,龚建周,刘毅华，郑金育

(1.广州大学地理科学学院,广东 广州 510006；2.广东园梦生态环境工程规划设计有限公司，广东 惠州 516000)

随着城镇化快速发展,土地利用强度日益增大,过度开发土地给生态环境造成极大的压力[1]。快速城市化区域土地资源利用约束下的生态环境与经济效益协调发展关系也成为土地科学的焦点[2]。通过模拟分析生态系统服务价值与区域经济发展协调度的时空演变格局和空间异质性,合理优化土地利用布局和土地资源配置,实现生态和经济效益最大化,有利于快速城市化区域的自然生态系统和社会经济系统长期稳定发展[3]。国外对于生态服务价值定量化测算[4]、变化动态[5]、空间演变[6]等研究有较为成熟的成果,已有研究主要探讨生态工业园区[7]、生态村[8]或者生态城市[9]的生态经济协调度空间分布特征。国内学者对于生态服务价值定量化测算[10]进行了深入探讨,并针对生态与经济协调度的定量评价开展模型和方法研究,如生态足迹法[11]、综合评价指标法[12]以及引用年度经济环境协调度指数(YEEH)和经济环境协调度变化指数(EEHC)等模型[13],同时运用GIS手段进行不同地区生态系统服务价值和生态经济协调度的时空演变分析,主要探讨中国陆地[14]、自然保护区[15]、城市群[16]等区域,对快速城市化区域的研究相对较少。国内外学者较多基于宏观视角,以县域为单元分析生态经济协调度的时空演变[17],并以过去一定时段的数据进行时空分布演变分析[18],然而历史时期的生态经济协调度演变规律分析存在一定的滞后效应,计算结果与未来政策实施背景有一定差异。因此,需要基于未来土地利用预测数据,对区域生态经济协调度的聚集规律、关联模式等分布特征展开更为细致的定量研究和模拟预测,进一步预测未来生态经济协调度的分布规律,尽可能贴近未来规划背景,为区域生态系统保护和优化管理决策提供更好的科学支撑。目前基于CA-Markov模型对土地利用进行预测的研究较多,也有部分学者用此模型预测未来的生态系统服务价值,但运用此模型预测未来生态经济协调度的研究较少。

近20 a惠州市经济快速发展,生态用地因城市扩张被割裂、占用,城市景观破碎化严重,经济效益提升给城市生态安全带来巨大承载压力。因此研究该市的环境保护与经济效益的协同耦合发展是目前亟需解决的重要问题。笔者以惠州市为研究对象,基于耦合协调视角,通过CA-Markov预测惠州市2025年土地利用数据,并引入EEHC模型研究该区域生态与经济的耦合协调关系,基于重采样的网格单元,微观模拟和分析土地利用生态服务价值与区域经济发展的协调度、空间演变格局和聚集规律等,为促进土地安全、生态环境和经济效益协调发展研究提供理论依据,对保障自然、经济和社会系统稳定发展,实现粤港澳大湾区以及其他快速城市化地区的可持续发展具有重要指导意义。

1 研究方法与数据来源

1.1 研究区概况

惠州市位于珠三角东部(图1),地处北纬22°24′～23°57′,东经113°51′～115°28′之间,总面积约1.13万km2,地势地貌较为复杂,属亚热带海洋性季风气候区。该市GDP从2012年的2 407.0亿元上升至2017年的3 830.6亿元,第三产业产值从2012年的893.2亿元上升至1 559.6亿元,表明该市经济基础较好,产业结构正向第三产业主导转型。惠州市属于广东省新兴的工业化城市,处于快速城市化阶段,耕地和林地大量转化为建设用地,是快速城市化区域的典型代表。

1.2 数据来源与处理

1.2.1数据来源

研究选用2000年11月8日、2005年12月1日、2010年9月12日Landsat系列的TM影像和2015年11月6日的OLI影像作为主要土地利用数据源。社会经济数据主要以地区GDP以及粮食作物单产、面积、价格等为代表,来源于2018年《惠州市统计年鉴》及相关经济统计公报。

图1 惠州市区位图

1.2.2数据处理

惠州市秋冬季节云量较少,适合运用ENVI图像处理平台进行遥感影像解译分类,并按照全国土地利用分类标准划分地类[19],利用Google Earth高分辨率遥感影像、野外考察GPS样点及土地利用现状图, 对分类结果进行精度验证, 4期分类的Kappa系数均在85%以上, 符合分类要求。同时为了打破行政单元的限制,更好地体现空间上的差异性,细致分析该区域生态与经济协调度的空间分异特征,运用ArcGIS软件的Fishnet和叠加分析工具,将各时点惠州市土地利用数据属性赋值于1 km×1 km的格网单元内,统一投影为WGS_1984_UTM_Zone_49N坐标系。运用线性函数对2012—2017年的GDP数据进行线性回归处理,回归方程为y=268.69x-538 177(y为GDP,x为年份),R2=0.987 7,据此预测2025年GDP为5 920.25亿元。参考张元栋等[20]基于格网的GDP空间化方法,首先使用2000、2005、2010和2015年惠州市各区县GDP除以建设用地面积求取地均GDP,然后用每个格网建设用地面积乘以地均GDP得到各格网的GDP。

1.3 研究方法

1.3.1生态系统服务价值计算方法

目前普遍认为根据COSTANZA等[21]确定的基准单价进行生态系统服务价值计算较为权威。谢高地等[22]基于其成果,依据中国实际情况得出了中国生态系统单位面积的服务价值当量表。基于其标准当量表,同时依据惠州市经济发展水平,利用2000—2015年研究区主要粮食作物的平均年单产、收购价格以及种植面积,计算每1 hm2农田粮食作物的年经济价值,用于修正惠州市单位面积生态系统服务价值,使其更符合研究区现状。修正模型参照丁丽莲等[23]的研究成果。经过实际修正之后,惠州市单位面积的生态系统服务价值系数见表1。

表1 惠州市单位面积生态系统服务价值系数

基于生态系统服务价值系数,计算该区域生态系统服务价值(ESV,VES),公式为

(1)

式(1)中,Aj为j类生态系统面积，hm2;CV，ij为j类生态系统i类生态服务功能的价值系数，元·hm-2。

1.3.2生态服务价值与经济协调度模型

当前生态保护与经济效益的协调发展程度没有统一的衡量标准,普遍运用生态与经济协调度等类似的相对指标表征生态服务价值与经济发展水平之间的相互协调关系。以乔标等[24]和魏伟等[25]的研究为基础,构建经济环境协调变化指数(EEHC,CEEH),以反映土地资源约束下生态环境与经济发展过程中的相互影响和耦合互动程度。

(2)

(3)

(4)

式(2)～(4)中,RESV为格网单元的生态系统服务价值变化率，%;RESV,m和RESV,n分别为不同时期初末格网单元的生态系统服务价值，元;RGDP为格网单元的GDP变化率，%;RGDP,m和RGDP,n分别为不同时期初末格网单元的GDP值，元。

由于惠州市在2000—2025年处于经济持续增长时期,RGDP为正值,EEHC值存在2种情况:正值表明该区域生态与经济发展较协调,负值表明该区域生态与经济发展不协调。参考前人研究成果[26],根据EEHC值区分生态与经济的协调程度:≥1,表明该时段经济增长速度不高于ESV增长速度,处于高度协调状态;0～<1,表明经济增长速度高于ESV增长速度,但仍处于低度协调状态;-1～<0,表明ESV呈负增长,但速度小于经济增长速度,为低度失调状态;<-1,表明ESV下降速度大于经济增长速度,属于高度失调状态。

1.3.3CA-Markov模拟预测

土地利用转移矩阵可表征某一时期区域内各土地利用类型之间的转换关系。利用转移矩阵分析土地利用结构特征及其变化,并对各类土地相互转化进行追踪研究。转移矩阵的具体表达式及说明参照张晓娟等[27]的研究。

CA模型的优点是能够较好地模拟一定复杂程度系统的时空演变过程,缺点是模型转换规则构建成功与否决定了模拟精度的高低,但转换规则构建难度大且变化数量无法精确模拟[28]。Markov模型通过转移概率矩阵实现数量变化的模拟预测,然而缺乏相关空间变量,无法得知研究对象的空间变化程度[29]。因此,笔者将CA模型的空间拓展特性和Markov模型的数量预测2个优势结合起来,可以实现土地利用变化的空间和数量模拟。

元胞邻域如要发生理想状态的转化,需构建转换规则,规则表达式为

St+1=f(St,N)。

(5)

式(5)中,t和t+1为元胞所处的前后时刻;S为元胞状态聚集;f为转化规则;N为元胞邻域。

Markov模型公式为

St+1=PijSt,

(6)

(7)

(8)

式(6)～(8)中,St为随机事件在t时刻的状态;St+1为随机事件在t+1时刻的状态;P为转移矩阵,取值为0～1;n为土地利用类型。

CA-Markov具体参数设置:Markov模块,将影像时间间隔和周期均设置为5 a,比例误差设置为0.15;适宜性图集选择因子为坡度、离城镇建设用地距离、离主要道路距离和离水域距离,同时将水域设置为限制条件;滤波器为5×5。以2005和2010年遥感影像数据为基础,模拟2015年土地利用变化时,起始时刻为2010年,CA循环次数为5次;预测2025年土地利用变化时,起始时刻为2015年,CA循环次数为10次。

1.3.4空间自相关模型Getis-Ord Gi*

生态与经济协调度的空间分异特征可用空间自相关的局部指标进行分析,如果空间单元与相邻单元属性值较为接近,则为空间正相关,反之为空间负相关。该研究运用Getis-Ord Gi*指数计算某一空间单元和邻近单元关于EEHC属性的关联程度[30],分析相邻格网EEHC年际变化的高低空间聚类分布,用于识别惠州市生态环境与经济发展协调度的热点区和冷点区,在一定程度上反映生态与经济协调度的空间分布关系。

(9)

式(9)中,Xi为格网i的EEHC值;Wij(d)为格网i和j之间的空间权重矩阵,空间相邻为1,不相邻为0。若Gi*(d)为高正值,表明i格网单元和相邻单元的EEHC值高于周边,属热点区域,在空间上表现为高值聚集;反之则为冷点区域。

2 结果与分析

2.1 2000—2025年惠州市土地利用变化

首先,基于2005和2010年土地分类数据模拟预测2015年土地利用情况,对比2015年实际土地分类数据,检验Kappa系数达89%,预测精度符合要求。其次,以2010和2015年土地利用数据为基础,运用CA-Markov模型预测2025年惠州市土地利用状况(表2、图2)。

从图2可知,该市2000年建设用地较少,主要集中在惠阳区。随着城市化进程加快,2005年惠城区和惠阳区经济得到快速发展,耕地、林地因城市经济发展需要被大量侵占,同时建设用地沿着交通干线纵轴直线蔓延。2005—2015年是惠州市经济高度提升时期,建成区不断向北延伸,博罗县和龙门县大量耕地和林地等生态用地被占用为建设用地,景观破碎程度愈加严重。2005—2010年惠州市建设用地面积急剧上升,但2010年后扩张速度减缓,土地利用变化强度大为减弱。

纵观2000—2025年土地利用变化转移矩阵(表2),建设用地扩张最为突出,其中耕地和林地被侵占程度最高,面积分别为1 643.61和639.97 km2;耕地面积从2000年的2 804.69 km2下降至2025年的1 502.50 km2,减少46.42%;林地面积从2000年的6 406.37 km2下降至2025年的5 837.83 km2,减少8.87%。将2025年惠州市土地利用格局的预测结果与土地利用现状相比,可知2025年时该区域土地利用更为高效,生态破碎化现象得到有效缓解。

表2 2000—2025年惠州市土地利用转移矩阵

2.2 惠州市生态与经济协调度时空演变及模拟

基于2000—2025年惠州市生态服务价值和GDP动态变化数据,利用EEHC指数反映经济与生态发展变化过程的时空耦合协调度。从图3可知,2000—2010年,该区域的生态与经济协调度总体上处于较低水平,惠城区、惠阳区、博罗县大部分区域和龙门县东部处于高度失调区域,高度协调和低度协调区域面积较少,共占7.81%,惠东县大部分区域和龙门县东部区域处于低度失调区域。2010—2025年,该区域的生态与经济协调程度总体较为低下,大部分经济发达区域仍处于高度失调阶段,城区部分零散区域表现为高度协调和低度协调,同时惠东县和龙门县部分零散区域从低度失调区转化为低度协调区,虽然低度失调和高度失调区在往惠州市的生态屏障方向,即西北方向和东部扩张,但低度协调和高度协调区域面积也开始增加,零散分布在惠东县、龙门县以及博罗县部分区域。纵观整个研究时段,2000—2025年的EEHC指数整体处于较低水平,低度协调区和高度协调区分布较为散乱,高度失调区主要位于惠城区、惠阳区等经济发达地区,低度协调区主要分布在龙门县、惠东县以及博罗县地区。

从表3可知,2000—2010年高度失调、低度失调和低度协调区域面积分别为1 771.17、8 671.86和884.82 km2,面积占比分别为15.64%、76.55%和7.81%。该时期惠州市过于注重经济发展,区域生态系统服务价值功能发挥受到抑制,导致经济增长速度高于生态服务价值提升速度,生态与经济的协调程度不断下降。2010—2025年高度失调、低度失调、低度协调和高度协调区域面积分别为2 707.17、6 845.86、1 132.38和642.44 km2,面积占比分别为23.90%、60.43%、10%和5.67%。生态与经济的协调发展程度仍然较差,形势比较严峻。但局部出现高度协调区域,主要是受益于惠城区和惠阳区良好的经济发展状况,其生态用地较少受到侵占,所以部分低度失调区转化为低度协调和高度协调区。2000—2025年低度协调区和高度协调区主要分布在惠东、惠阳、博罗县等区域,少许散落于城区部分区域,面积分别为1 527.55和366.79 km2,面积占比分别为13.48%和3.24%;高度失调区主要位于惠城区和惠阳区等经济发达地区,面积为3 273.49 km2,面积占比为28.90%;低度失调区主要分布在龙门县和惠东县,面积为6 160.02 km2,面积占比为54.38%。结果表明,惠州市失调区域面积大于协调区域,且在空间分布上失衡,两极化现象较为严重,该市生态环境和经济发展的协调水平总体上不高,但有局部提升且向全市发展的趋势,说明政府采取的生态保护措施开始显现出生态和经济效益。

图2 惠州市2000—2015年土地利用格局及2025年土地利用预测

图3 2000—2025年惠州市经济环境协调变化指数(EEHC)时空演变及模拟

表3 2000—2025年惠州市环境与经济不同协调程度区域面积及占比

2.3 惠州市生态与经济协调度空间演变聚集和分异特征

通过Getis-Ord Gi*对2000—2025年惠州市EEHC时空变化差异进行比较,以1 km格网内EEHC作为研究对象,相应统计2000—2010年和2010—2025年这2个相邻时段EEHC的Gi*热度,基于置信区间将其划分为7类,其中Z值>2.58为热点区,Z值<-2.58为冷点区。从图4和表4可知,该时期总体呈现冷点聚集分布模式,冷点区域占20.19%,主要分布在惠城区、惠阳区、博罗县东西部和龙门县东部,几乎没有热点出现。2010—2025年,该市表现为冷点聚集主体分布、热点聚集散乱分布的总体特征,冷点仍然主要分布在主城区域,热点分布在冷点周围,占比为9.22%,EEHC波动起伏较为剧烈,热点集聚效应逐渐加强,聚集特征也更为明显。从2000—2025年整个时期来看,冷点区域占比为11.61%,热点区域占比为2.65%。惠州市的EEHC空间分布模式从2000年的冷点单一聚集发展为2025年的冷热点2种集聚模式共同分布,同时热点向冷点逐渐迁移,进一步说明惠州市经济发展和生态环境的协调度朝着两极化方向发展,虽然目前以生态经济低协调或者高度不协调为主,但未来趋势为从不协调向高度协调的良性转变。

图4 2000—2025年惠州市经济环境协调变化指数(EEHC)冷、热点区域分布

表4 2000—2025年惠州市冷点、热点区域面积与比例

3 讨论

3.1 生态系统服务价值预测结果探讨

目前对于生态系统服务价值以及生态与经济协调度的研究主要是对历史时期的土地利用数据进行分析,笔者引入CA-Markov模型预测未来土地利用数据,探索未来生态与经济的协调发展趋势,为目前快速城市化地区生态修复、改善生态系统服务功能发挥提供一定参考。但是基于土地利用数据预测值难以精准计算未来的生态系统服务价值,笔者参考其他研究,设置一定数量的限制因子和条件对未来土地利用状况进行预测[30],使其更符合未来发展趋势,预测结果较为合理,但较少考虑政策发展、规划执行等其他因素的影响。该研究结合粮食产量、价格和播种面积等因素修正生态系统服务价值,使其更符合区域实际情况,该方法被广泛运用于生态系统服务价值计算中[31-32],结果较为合理。生态系统服务价值的计算较为依赖当年土地利用数据,虽然解译多时段的土地利用数据能更细致地体现时空变化,但由于客观条件限制,精确数据获取难度较大,通过30 m分辨率的遥感影像难以将土地利用类型划分更为精细的二级地类,该研究仅划分为一级地类,导致结果有一定偏差,但生态系统服务价值总体偏差较小,加之指标RESV为相对值,因此该误差可忽略不计,但在后续研究中应进一步细化土地利用分类,提高研究结果的精准性。

3.2 GDP预测结果探讨

计算生态与经济协调度时,笔者在总结前人研究的基础上引用EEHC模型进行简单评价,其客观性需要进一步研究。笔者采用的GDP统计数据的准确性达到研究要求,但是GDP统计来源多样,还需增加其余经济指标以便更准确地反映经济水平,这是未来需要探索的重点。对2012—2017年GDP进行线性回归处理,虽然线性回归结果较好(R2=0.987 7),但基于此对2025年GDP进行预测,预测结果仍不够精确,如何引入其他模型预测并修正GDP,使空间化结果更为具体,也是未来关注的重点。

3.3 评价单元选择导致的尺度效应问题

目前国内关于生态与经济协调度的研究主要局限于全国、省域、长江经济带或者城市群等尺度,对于市域尺度的研究鲜有报道,笔者研究也是对此进行初步探索。以往研究尺度较为宏观,且考虑到经济数据的可获取性,以及评价单元需有权利调配资源,所以多数以县域为评价单元,但是研究结果难以体现县域内部生态经济协调度的空间差异。因此该研究以格网为评价单元,打破行政单元的限制,可以更好地体现空间上的差异性,更细致地分析生态与经济协调度的空间分异特征。参考已有文献,发现常用的格网单元有 500 m×500 m[33]、1 km×1 km[26]、2 km×2 km[34]、3 km×3 km[35]、5 km×5 km[36]等,考虑到图幅限制以及张元栋等[20]得出的1 km尺度格网GDP空间化结果误差最小的结论,笔者最终选择1 km格网作为评价单元。然而评价单元和研究尺度不同,研究结论也具有一定差异。如魏伟等[25]开展中国陆地经济与生态环境协调发展的空间演变研究,以1 km为格网单元,认为东部经济发展较好的区域协调度总体较差,中部和西部等经济发展速度相对缓慢的区域协调度好于东部,笔者的研究结论与其大致相符,该研究中经济发展较好的区域大部分为生态与经济失调区域。王振波等[26]以县域为研究单元,开展1991年以来长三角快速城市化地区生态经济系统协调度评价及其空间演化模式研究,认为研究区空间异质性较高,内部空间差异较大。该研究与其部分结论不一致,笔者认为惠州市两极分化严重,内部空间分异不大。由于目前评价单元没有统一标准,尺度效应问题仍未得到解决,生态与经济协调度的尺度效应值得进一步深入研究。

4 结论

土地利用、生态环境和经济社会系统3者密切相关,该研究基于2000—2015年惠州市遥感影像获取各年土地利用数据,运用CA-Markov模型模拟预测2025年土地利用状况,通过利用不同时段土地利用数据计算生态系统服务价值,同时运用EEHC等生态经济协调度模型,探讨在土地利用资源约束下快速城市化地区的生态环境与经济发展协调度的时空演变格局。

(1)纵观2000—2025年惠州市土地利用格局,建设用地增加2 240.76 km2,林地和耕地共减少1 870.73 km2,未利用地逐年减少,表明该区域土地利用效率较高,虽然景观破碎化现象突出,但相对于前期已得到有效缓解。这反映了快速城市化区域的土地利用正在转型,同时当地政府更加注重协同发展生态与经济,这有助于经济、社会、生态效益的统一提升。

(2)近20 a惠州市EEHC变化较大,失调区域面积大于协调区域,空间分布上失衡且两极化现象较为严重,生态与经济协调水平总体上不高,发达地区周边局部地区从不协调向协调转变,并有逐渐扩大趋势,表明快速城市化区域愈发重视生态与经济的协调发展,未来协调发展重心应放在城市用地之外的生态空间中。

(3)基于格网分析法和Getis-Ord Gi*,从微观角度分析2000—2025年惠州市生态与经济协调发展的空间格局及其集聚特征,EEHC空间分布模式从2000年的冷点单一聚集发展为2025年的冷热点2种集聚模式共同分布,同时热点向冷点逐渐迁移,表明虽然目前惠州市空间差异不均衡和两极化发展问题较为突出,但生态与经济协同发展趋势在向良性发展。发达地区在发展经济的同时需要加大生态修复和保护力度,实现生态与经济协同可持续发展。