土地利用生态风险评价是针对土地利用过程中发生不利生态影响可能性的评价，对认知土地利用生态风险产生机制，降低生态风险概率以及优化土地利用具有重要意义^[1-2]，随着城市化进程的快速发展，导致环境管理目标的转变，多风险源、多受体及区域景观水平阶段已经成为当前土地利用生态风险评价研究的重点^[3]。而基于人类开发导致的景观格局转变、土地利用类型变化的生态风险评价更加关注生态风险的时空异质性，其重点为评价人类活动在生态系统中产生的不良影响^[4]。近年来，由于大、中型城市建设以及“城乡融合”等战略的实施，土地开发活动愈加频繁，导致各种生态环境问题的出现，生态风险评价逐渐成为研究和解决环境问题的重要手段。国内大量学者从不同角度对区域生态风险进行研究，提出了不同的生态风险评价模型：包括基于P-S-R模型设计土地利用生态安全(LUES)分级标准构建生态风险模型^[5-6]，利用土地类型面积比重构建生态安全评价模型^[7-8]，从景观结构角度利用景观扰动指数与景观脆弱指数构建模型对区域生态风险进行评价等^[9-11]，其中从景观角度出发进行生态风险的研究越来越多，研究尺度大到城市群^[12]、市域^[13-15]、小到县域^[4]、湿地公园等^[11]，研究区域包括流域^[16-17]、干旱区^[10]、海岸带^[18]、丘陵区^[19]以及农林交错带等^[9]，尽管前人研究尺度以及研究区域包含甚广，但对黄河沿岸土地利用生态风险的研究还有所缺乏，黄河流域作为中华民族的发源地，从西到东横跨9个省区以及青藏高原、内蒙古高原、黄土高原和黄淮海平原4个地貌单元，由于河水冲积水土流失严重，土地类型变化频繁。2019年9月习近平总书记在河南郑州主持召开黄河流域生态保护和高质量发展座谈会并发表重要讲话，明确指出黄河流域在我国经济社会发展和生态安全方面具有十分重要的地位，深刻阐明黄河流域生态保护和高质量发展的重大意义，因此研究黄河流域沿岸土地利用生态风险变化特征，对加强黄河治理保护、降低黄河沿岸土地利用生态风险具有重大意义。为此，本文以河南省黄河沿岸7个城市为研究对象，以2005，2010，2015以及2018年4期土地利用数据为基础，根据景观扰动指数与景观脆弱指数构建生态风险评价模型，借助土地转移矩阵与空间自分析方法对土地利用生态风险进行分析与评价，研究结果可为正确认识黄河流域沿岸城市土地利用生态风险，对土地景观资源保护、生态系统优化与生态环境保护、管控措施制定等提供一定的理论依据。

1 数据和方法 1.1 研究区概况

黄河流域途径河南省7个市区，包括三门峡市、洛阳市、焦作市、新乡市、郑州市、开封市以及濮阳市，属北亚热带向暖温带过渡的大陆性季风气候，南北过渡性明显。研究区域大多为平原、丘陵地貌，耕地是主要的景观基质，林地和草地主要位于研究区的西南方向，以郑州为中心的7个市区，近几年来发展迅速，建设用地面积快速增长。

1.2 数据来源

本文所采用的2005，2010，2015和2018年土地利用类型栅格数据均来源于中国科学院资源环境科学数据中心。所用数据的分辨率为30 m，在此基础上参考《土地利用现状分类(GB/T21010-2017)》并结合研究区实际情况将研究区土地类型分为耕地、林地、草地、水域、建设用地以及未利用地(见封3，附图 4)。

1.3 研究方法 1.3.1 构建生态风险指数

(1) 景观扰动指数(E_i)。同一地区的生物多样性、生态环境整体结构和功能、生态系统的更新和演替受到不同景观类型构成的影响，表现出不同的特征，同时生态系统抗干扰能力的大小也与该地区景观类型的构成密切相关^[4]。景观扰动指数(E_i)用来反映外界干扰(主要是人类活动)对生态环境(由不同的景观类型组成)的影响程度，可以通过对景观破碎度指数(C_i)、景观分离指数(N_i)和景观优势指数(D_i)三者赋予权重累加获得，以反映出不同景观构成的干扰程度，景观扰动指数的大小主要受人类活动的影响，受人类开发活动影响越大，扰动指数越大，随之生态风险越大，反之则越小^{[2, 4, 9-10]}，景观扰动指数(E_i)可以用下面的公式表示：

$ {E_i} = a{C_i} + b{N_i} + c{D_i} $

(1)

式中：a, b, c为各景观指数权重, 且其和为1, 根据前人研究结果^[9-11]并结合研究区特点, 认为景观破碎度的分布特征比较重要, 分离度重要程度, 结合专家打分法将景观破碎度权重赋值为0.5;分离度赋值为0.3;优势度赋值为0.2。

景观破碎度指数(C_i)可表示为：

$ {C_i} = \frac{{{n_i}}}{{{A_i}}} $

(2)

景观分离指数(N_i)可表示为：

$ {N_i} = \frac{A}{{2{A_i}}}\sqrt {\frac{{{n_i}}}{A}} $

(3)

景观优势指数(D_i)可表示为：

$ {D_i} = \frac{{{Q_i} + {M_i} + 2{L_i}}}{4} $

(4)

式中：A_i为景观类型i的面积；A为景观总面积; n_i为景观类型i的斑块数; Q_i为斑块i出现的风险小区数/总风险小区数; M_i为斑块i的数目/斑块总数; L_i为斑块i的面积/风险小区的总面积^[2]。

(2) 景观脆弱指数(F_i)。景观脆弱指数用来表示不同生态系统的易损性，不仅反映了土地本身的自然属性，更是人为因素和自然因素综合作用的结果，景观脆弱指数越大，表明对区域生态系统的干扰越大，区域生态风险越高。结合前人研究以及研究区实际情况认为建设用地最为稳定，不易产生土地利用类型的变化，故将建设用地赋值为1；未利用地是指农用地和建设用地以外的荒草地、沙地、裸岩、裸土地等，受环境变化、人类活动影响较大，土地利用类型改变相对容易，其景观脆弱性最强，最为敏感，故将未利用地赋值为6；水域、耕地、林地、草地则根据前人研究结果和各自脆弱性特点分别赋值为5, 4, 3, 2^[2]。进行归一化处理后得到研究区的各土地类型景观脆弱指数。

(3) 生态风险指数(ERI)。在对上述景观扰动指数和景观脆弱指数研究分析的基础上，为了描述样地内生态损失的大小，通过采样方法将空间格局转变为生态风险变量，以此构建土地利用生态风险指数ERI。

$ {\text{ERI}} = \sum\limits_{i = 1}^n {\frac{{{S_{ki}}}}{{{S_k}}}} \sqrt {{E_i}{F_i}} $

(5)

式中：n为景观类型数量; S_ki为第k个风险小区第i类景观组分的面积; S_k为第k个风险小区总面积。

1.3.2 空间自相关分析

(1) 全局Moran’s I指数。空间自相关分析的目的是度量某一变量是否在空间上相关，并描述其相关程度如何。本研究采用Moran’s I指数图对土地利用生态风险全局空间自相关进行分析描述。具体计算公式^[20]为：

$ \begin{gathered} {\text{Global}}\;{\text{Moran's}}\;I = \hfill \\ \frac{n}{{\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{W_{ij}}} } }} \times \frac{{\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{W_{ij}}} } \left( {{y_i} - \bar y} \right)\left( {{y_j} - \bar y} \right)}}{{\sum\limits_{i = 1}^n {{{\left( {{y_i} - \bar y} \right)}^2}} }}(i \ne j) \hfill \\ \end{gathered} $

(6)

式中：n为空间单元总个数；y_i，y_j是变量y在相邻配对空间单元(或栅格细胞)的取值; 是所有单元的平均值；W_ij是单元i, j的空间权重矩阵；Moran’s I指数的取值范围在-1~1之间。

(2) 局部Moran’s I指数。利用局部Moran’s I指标局部空间自相关性，对空间集聚和空间异质进行识别，并生成LISA聚类图。局部Moran’s I指数计算公式为^[20]：

$ \begin{gathered} {\text{Local Moran's}}\;I = (\frac{{{y_i} - \bar y}}{m})\sum\limits_{j = 1}^n {{W_{ij}}} {\left( {{y_i} - \bar y} \right)^2} \hfill \\ (i \ne j) \hfill \\ \end{gathered} $

(7)

$ m = \frac{{\sum\limits_{i = 1}^n {{{\left( {{y_i} - \bar y} \right)}^2}} }}{n} $

(8)

式中：y_i，y_j是标准化后的空间单元观测值；其余变量与全局Moran’s I指数变量含义相同。

1.3.3 基于格网的生态风险小区划分与生态风险空间分析

本文根据研究区特点，在参考前人研究以及考虑研究区空间异质性和斑块大小的基础上^{[4, 9, 12]}，借鉴传统景观生态学中关于景观取样大小的方法，即网格划分大小应当为斑块平均面积的2~5倍，为减少网格对大斑块的切割，选取耕地平均斑块面积的5倍，利用GIS软件将研究区划分6 km×6 km的采样网格，采样方式为等间距系统采样法，共计1 737个生态风险小区。在利用空间自相关分析证明了生态风险评价值有空间自相关性的基础上，在GIS中计算每个生态风险小区的生态风险值，然后将每个网格的值赋予该网格的中心点，最后采用克里金插值法进行空间插值得到研究区生态风险的空间分布值。

2 结果与分析 2.1 土地利用变化分析

通过2005-2018年土地利用转移矩阵表(表 1)分析可知，从总体上看，13 a间，研究区域内建设用地面积增加了38.05%，增长趋势较快，这与城市化的快速发展、以及城市的扩张是分不开的，其中草地面积减少了10.70%、林地面积减少了2.41%，一部分是土地开垦导致的面积减少，还有一部分是建设用地的扩张占用了城市周围的草地和林地；水域面积减少了4.81%，这可能是由于河流水位的逐渐下降，开垦河畔导致水域面积减少；未利用地面积减少了34.75%，大部分是由于开垦耕地导致的。从土地类型转变上看，13 a间有277 495.47hm²耕地转化为建设用地，这是建设用地增加的主要来源，表明早期的城市扩张大多是以牺牲耕地为主；在研究时期内，有96 759.81 hm²的建设用地转化为耕地，其中主要是农村居民点用地转化为耕地，这是由于在“城乡融合”、“美丽乡村建设”和“乡村振兴”等战略的实施下，河南省开展大量村庄规划、整治工作，在土地整治、迁村并点等具体工程中使得原本凌乱分布、用地浪费的农村居民点变得更加集中，同时将废的宅基地复垦为耕地，使得耕地面积增加；而草地、林地和水域分别有44 580.87，35 392.32，31 014.36 hm²转化为耕地，是补充建设用地占用耕地的主要来源。

2.2 景观格局指数变化分析

由表 2景观指数统计结果可知，2005-2018年，耕地是研究区域内最主要的景观，班块面积最大，分布最为广泛；其次是建设用地。从破碎度和分离度来看：在研究期限内耕地的破碎度先减小后保持不变，分离度先减小后增加，减小的拐点都出现在2010年，这是由于在2010年耕地斑块数量出现较大幅度减少；建设用地面积和斑块数量不断增加，破碎度和分离度相对2005年都有所降低，这是由于建设用地的扩张占用了周围其他用地类型，使的建设用地集聚性不断增强；林地面积和斑块数量在2005-2010年期间减少较快，2010-2018年期间变化幅度较小，所以其破碎度和分离度先减小后基本保持不变；草地面积整体上在减小，而斑块数量变化幅度不大，因此其破碎度和分离度整体上在增大；水域面积先增后减，而斑块数量一直在增多，这可能是水位下降导致水域连通性减低，导致其破碎度和分离度在不断增加。

2.3 生态风险分析 2.3.1 土地利用生态风险空间自相关分析

根据每个网格计算的生态风险值，利用GeoDa软件计算出研究区域内生态风险指数的2005，2010，2015以及2018年4期全局Moran’s I值，分别为0.561 9，0.573 7，0.568 2以及0.566 3，且通过显著性检验，结果表明土地利用生态风险值在空间上呈现出显著的正相关性，存在空间集聚特征，同时为下一步的插值分析奠定了理论前提。与全局自相关相似，利用GeoDa软件对4个时期土地利用生态风险指数的局部空间聚集情况进行分析，得到LISA聚类图(图 1)，由图 1可知，研究区域内以“高-高”、“低-低”聚集模式为主，说明土地利用生态风险急剧变化的区域相对较小，“高-高”区域主要分布在黄河沿岸以及农村居民点散乱分布的区域，“低-低”区域主要分布在研究区域边界处，主要是大面积的林地以及草地分布区。

2.3.2 土地利用生态风险时空变化分析

根据每个网格计算的生态风险值，将值赋予网格中心点，采用克里金插值法进行空间插值得到4个时期土地利用生态风险值的空间分布图(图 2)。利用自然断点法将研究区划分为5个生态风险等级：低风险区(0.047≤ERI＜0.087)，较低风险区(0.087≤ERI＜0.105)，中风险区(0.105≤ERI＜0.119)，较高风险区(0.119≤ERI＜0.131)，高风险区(0.131≤ERI＜0.151)，结果见图 2。从生态风险值大小来看：2005-2018年，4个时期生态风险最小值依次为0.047 2，0.047 5，0.047 4以及0.047 6，最大值依次为0.147 9，0.150 5，0.146 2以及0.150 7，最小值和最大值大致呈N形变化，并且变动幅度较小。从空间分布来看，在4个时期内，高风险区分布区域大致相同，呈条状和块状分布，条状大多分布在黄河沿岸，由于黄河水位的变化以及河水对两岸土壤的冲击，土地利用类型不断变化，大部分在耕地和水域之间转变，使得人类干扰活动增强，导致土地利用生态风险较高；块状高风险区大多分布在农村居民点零散分部的平原区，由于居民点斑块分布过于零散，破碎度以及分离度较高，使得干扰度较高，土地利用生态风险值也高。低和较低风险区大多分布在研究区域四周以及西南大部分区域，其中西南、西北区域大部分为大面积的林地和草地，受人为干扰较小，生态风险值低；其他低和较低风险区所处的区域耕地和建设用地是主要的景观基质，并且农村居民点面积都较大，耕地和建设用地破碎度和分离度较小，生态风险值较低；在较高风险区包裹的研究区域中部出现部分较低风险区，大多位于各城市建成区，由于大面积都是建设用地，受人类干扰后的损失度低，生态风险值相对较低。

由表 3各年份土地利用生态风险等级面积分析可知，在研究时限内，较低风险区和较高风险区所占面积最大，其次是低风险区和中等风险区，面积最小的是高风险区；从面积变化来看，低风险区和中等风险区面积在逐年增加，分别增加了335 223.45和149 024.16 hm²，较低风险区面积先减少后增加，总体上较2005年减少了143 642.72 hm²；高风险区和较高风险区面积在逐年减少，分别减少了179 419.86和161 185.05 hm²，风险区面积变化说明，研究区内土地利用生态风险整体上在不断改善；由表 4各类风险区面积转移矩阵分析可知，有342 231.48 hm²的较低风险区转化为低风险区，是低风险区面积增加的唯一来源；中等风险区主要向较低风险区转化，占中等风险区转出面积的98.4%；较高风险区转出对象涉及到3个风险区，最主要是中等风险区，面积达到325 801.08 hm²；50%左右的高风险区全部转化为较高风险区，其中块状高风险区面积减少较多，主要是因为城镇居民点斑块在扩张，农村居民点小斑块数目减少，大班块数目增加，耕地和建设用地的连通性都在增强，破碎度和分离度降低，生态风险值随之降低，其他块状高风险区面积减少的原因也与之相同。从总体转移方向来看，各风险区主要向相邻风险区转化，仅有少量地区出现跨区转移，说明研究区土地利用生态风险变化相对稳定，生态风险急剧变化的区域较小。

3 讨论与结论

(1) 2005-2018年，从地类变化来看：研究区域内建设用地面积增加，且增长趋势较快；其他几类用地面积均有不同程度的减少；从地类转移来看，耕地是建设用地增加的主要来源，林地、草地和水域主要转化为耕地，不同程度上补充了建设用地对耕地的占用。

(2) 在研究时段内，研究区域内四期生态风险指数全局Moran’s I值分别为0.561 9，0.573 7，0.568 2以及0.566 3，土地利用生态风险在空间上呈现出显著的正相关性，存在空间集聚特征；且LISA聚类图显示研究区域内“高-高”、“低-低”是生态风险主要的空间聚集模式，说明土地利用生态风险急剧变化的区域相对较小。

(3) 2005-2018年，生态风险最小值和最大值呈N形变化，且变化幅度较小；4个时期内，高风险区分布区域大致相同，呈条状和块状分布，条状大多分布在黄河沿岸，块状高风险区大多分布在农村居民点零散分部的平原区；低和较低风险区大多分布在研究区域四周以及西南大部分区域，这些区域大部分为大面积的林地和草地。在研究时限内，低风险区和中等风险区面积在逐年增加，其中较低风险区是低风险区面积增加的唯一来源，中等风险区转入的主要来源是较高风险区；较低风险区面积先减少后增加，总体呈减小趋势；高风险区和较高风险区面积在逐年减少；各风险区主要向相邻风险区转化，仅有少量地区出现跨区转移，说明研究区土地利用生态风险变化相对稳定，生态风险急剧变化的区域较小。

研究区地势西高东低，除三门峡和洛阳两市有大面积集中分布的林地和草地外，其余各市林地和草地面积相对较少，因此研究区内低和较低风险区主要集中分布在三门峡和洛阳两个地区；黄河沿岸高风险区主要位于焦作、新乡、郑州和开封4个地区，主要是由于这4个地区黄河沿岸多为耕地，河水冲积容易造成水土流失，并且水位下降使得沿岸土地类型在耕地和水域之间转变，受人为干扰活动较大，生态风险较高，并且在13 a间面积变化不大，非黄河沿岸的高风险区多处于农村居民点斑块破碎，且斑块面积较小的地区，主要分布在洛阳、新乡和开封，建设用地和耕地连通性弱，破碎度和分离度大，生态风险值较高，在研究时限内高风险区面积减少部分多来源于块状高风险区，尤其是焦作和新乡两市，自2010年开始高风险区几乎全部转化为较高风险区，这部分地区耕地和农村居民点用地是主要土地类型；在河南省开展迁村并点、土地整治等活动下，建设用地和耕地斑块总数量在减少，但大斑块数量在增多，因此连通性增强，生态风险值在降低，因此科学的村庄规划、整治工作在改善农村人居环境的同时也可以降低土地利用生态风险；洛阳、三门峡两个地区黄河沿岸多为林地和草地，可以降低水土流失风险，减小了土地类型转变的频率以及人为干扰活动的影响，生态风险值较低。因此就研究区而言，洛阳、三门峡两个地区应注重林地和草地的保护，稳定该地区黄河沿岸土地利用低风险现状，同时开展科学合理的村庄规划，可以有效降低平原地区土地利用生态风险值；焦作、新乡、郑州和开封应注重对黄河沿岸水土流失的治理，稳定黄河沿岸土地类型，降低人为活动干扰，同时对散乱分布的农村居民点进行科学合理的规划、整治，提高建设用地连通性，降低耕地破碎度，以达到降低土地利用生态风险的目的。