简 卿, 张立强, 张蓬涛, 赵 丽

(河北农业大学 国土资源学院, 河北 保定 071001)

生态风险评价是一种具有多风险因子、多风险受体，注重对复杂生态系统特征进行分析的综合评价，能够明晰制约区域生态持续性的主要生态要素与过程，已成为生态学研究的热点[1-2]。景观作为具有高度空间异质性的区域生态系统的集合，其整体结构的改变和动态过程均比较缓慢，但组成景观的空间组分在受到干扰时可以不同速度和强度发生变化。这就使得生态风险评价可以通过景观这一慢速改变量来衡量风险受体范围和危害状况等关键要素，景观生态风险评价为生态风险评价的研究提供了新的视角[2-3]。国内外学者十分关注景观生态风险研究，从不同的评价对象、评价方法等方面进行了探讨，并对景观格局变化进行量化分析，根据评价结果来表征风险值变化。评价单元根据研究区域和评价目的的差异而进行划分，热点研究区域包括流域[4-6]、城市地域及行政单元[7-9]、道路沿线[10]、自然保护区[11-13]、矿区[14-15]、海岸带[16]等，也有针对青藏高原东部等特定地貌区[17]的研究。评价方法总体上可以概括为基于源汇理论评价法[8-9]和基于景观格局指数评价法[4,6]两种。对于生态功能区的研究较少，且评价角度较为单一，基于景观格局指数的评价模型构建往往局限于通过景观干扰度和景观脆弱度计算景观损失度进而得到的景观生态风险值，其全面性有待进一步完善。借鉴知识管理对知识的分类，知识可分为显性知识和隐性知识两方面，风险作为一种认知，其本身也应具有显性风险和隐性风险两个方面。从显隐性的角度看，目前大部分研究都只分析了显性景观生态风险，而对隐性景观生态风险评价不全，所以从景观生态风险的显隐性两个方面出发才能完整地分析景观生态风险。

生态功能区可为自身和周边地域提供生态服务功能，对维护区域生态稳定和生态平衡具有重要意义。青龙满族自治县于2016年9月被认定为国家重点生态功能区，对于保障京津冀地区的生态安全具有重要意义。本文从风险的本质属性出发，提出显性景观生态风险和隐性景观生态风险的概念，并分别从显、隐性景观生态风险角度构建评价指标体系，以青龙满族自治县为研究区域，评价和揭示区域景观生态风险等级和时空分异，旨在为研究区未来的城镇发展和土地利用规划生态布局提供参考。

1 研究区概况

青龙满族自治县(118°34′—119°36′E，40°04′—40°36′N)，位于河北省东北部，隶属秦皇岛市，县域总面积为350 557.93 hm2，辖11乡14镇，常住人口506 989人。县域东、西、北三面高，南面较低，呈簸箕状分布。地貌主要是中低山和丘陵，占全县面积的80%左右，河谷面积占10%，耕地面积占10%左右，素有“八山一水一分田”之称。该县为暖温带半湿润大陆性季风气候，四季分明，年均温为9.3℃，年均降水量为759.1 mm，一日最大降水量为121.5 mm。境内河流总长478 km，其中较大的河流共有6条，最大最主要的为青龙河。土壤类型主要为褐土和棕壤土，其中褐土所占全县总面积比例达79.6%，主要分布在低山、丘陵地区，棕壤占17.6%，主要分布在较高的山地地区。该县森林覆盖率较高，达56.24%，且树种比较齐全，有野生动植物达590多种。但随着人口的快速增加，耕地减少速度过快，导致区域土地供需矛盾突出；人类活动对景观的干扰加剧，导致生态环境敏感度升高。

图1 青龙满族自治县地理位置

2 数据与方法

2.1 数据来源与处理

Landsat TM遥感影像和DEM影像数据来源于中国科学院计算机网络信息中心的地理空间数据云平台(http:∥www.gscloud.cn)。其中2001年、2007年和2013年Landsat TM影像空间分辨率为30 m，成像时间为8月和9月，每期2景，条带号/行编号为121/32，122/32。DEM数据来源于中国科学院计算机网络信息中心的地理空间数据云平台的ASTER GDEM数据，空间分辨率为30 m，共2景。3期NDVI数据来源于美国NASA网站(http:∥reverb.echo.nasa.gov)下载的MOD13Q1数据，空间分辨率为250 m，每期3景，成像时间为7月28日—8月12日。社会经济统计数据来源于青龙满族自治县统计资料。

土地利用类型的划分依据全国《土地利用现状分类》(GB/T21010—2017)，结合研究区实际土地利用结构，分为耕地、园地、林地、建设用地、未利用地和水域6大类。以2012年土地利用现状图为基础，经过定义坐标系、几何校正、融合、镶嵌和裁剪等预处理后，采用人机交互解译，分别得到3期土地利用现状数据，并抽取2%的图斑经过精度验证，解译精度均大于85%，可以满足景观分析的需要。

2.2 研究方法

基于研究区生态系统的景观类型，从风险的本质属性出发，提出显性景观生态风险和隐性景观生态风险的概念，并从这两个方面选取因子计算景观生态风险指数，建立景观生态风险评价模型；基于格网GIS划分风险小区，保证风险小区的区间异质性和区内同质性，求得每个风险小区的生态风险指数，以生态风险值作为小区中心点的值，在ArcGIS中利用普通克里金插值法，获得研究区生态风险空间分布格局图；结合自然断点分级法，将研究区生态风险水平划分为低风险、较低风险、中等风险、较高风险和高风险5个等级，并对其时空演变和空间相关性进行分析。

2.2.1 景观生态风险显、隐性分析 景观是由多种要素构成的、具有高度异质性的集合体，区域生态系统的变化反映在景观上就是景观要素的空间结构变化和景观功能的变化。自然状态下的生态系统，区域内景观要素形状相对规则、景观结构相对完整、连接度良好，但在外界干扰的长时期作用下，外在上会表现出景观斑块的破碎、同质景观的分离以及景观斑块形状复杂化加剧等，本文将这种已表现出来的景观斑块外在变化定义为显性景观生态风险；而在外界干扰下，区域内在的生态系统敏感性、不稳定性也随之增加，尽管这种风险在一定时期内没有外在的表现形式，但长久以往必然导致景观格局的变化、生物物种多样性的减少和生态风险的增高，本文将这种暂无外在表现形式但对景观结构和功能造成影响的风险定义为隐性景观生态风险。为保证评价对象的一致性，显隐性景观生态风险均指景观本身的风险变化，而非讨论外界干扰的显隐性和潜在性。隐性景观生态风险在达到一定阈值后会转化为显性景观生态风险。本文暂不考虑潜在的自然灾害。

显性风险和隐性风险共同构成了完整的区域景观生态风险，通过对显、隐性两方面的风险评价，最终得到研究区内景观生态风险的空间分布及演变。显性景观生态风险和隐性景观生态风险的评价指标选取主要从景观结构和景观功能两个方面考虑。

(1) 显性景观生态风险。

① 景观破碎度指数(Ci)用以表述在受到自然或人为干扰后一个完整的生态系统或该生态系统的组成部分被分割的破碎程度[18]，与景观生态风险指数呈正相关。其公式可以表述为:

Ci=Ni/Ai

(1)

式中：Ni为第i类景观斑块数目；Ai为第i类景观斑块面积。

② 景观分离度指数(Di)用以表述不同类型景观斑块在空间上的离散程度和相互演变的速度[18]，与景观生态风险指数呈正相关。其公式可以表述为：

(2)

式中：A为所有景观类型的总面积。

③ 景观分维数指数(Fi)用以描述景观斑块受干扰后形态变化的复杂程度，值域[1,2]，景观斑块的形状越复杂，景观分维数指数越高[19]，且与景观生态风险指数呈正相关。其公式可以表述为：

Fi=2ln(Pi/4)/lnAi

(3)

式中：Pi为第i类景观斑块的周长。

④ 植被覆盖度指数(PL)可以定量地说明植被的生长状况，反映了生态系统的生物量和生态系统稳定程度[20]，生物量越大往往其恢复力稳定性越强，对外界干扰的抵抗能力也越强。植被覆盖度指数与景观生态风险值呈负相关，值域[0,1]，其公式可以表示为：

PL=(NDVI-NDVIsoil)/(NDVIveg-NDVIsoil)

(4)

式中：NDVI为归一化植被指数；NDVIsoil为裸土或无植被覆盖区域的NDVI值；NDVIveg为纯植被区域的NDVI值。

(2) 隐性景观生态风险指数。

① 景观脆弱度指数(Vi)用以表述不同生态系统受外界干扰的易损性和敏感性，与景观生态风险指数呈正比。不同景观类型对外界干扰的敏感程度不同，抵抗外界的干扰能力也是不同的。景观类型抵抗外界干扰的能力越弱，景观脆弱性越高。景观类型的脆弱度和景观结构及功能都有着紧密的联系，但要对景观类型的脆弱度有准确的定位往往不容易得到[21]。根据前人研究成果[19,22]，生物多样性比较单一、食物链结构简单以及处于群落演替初级阶段的景观生态系统的脆弱度往往较高。因此，参考相关研究成果[4,6,19]，结合研究区实际，将6种土地利用景观类型脆弱度从低到高依次赋值1～6，得到各类景观类型脆弱度赋值如表1所示。

表1 景观脆弱度赋值标准

② 结构稳定性指数(SIi)用以表述土地利用景观状态的稳定程度[23-25]，该指数与景观类型受到干扰后的稳定性呈正相关，与景观生态风险指数呈负相关。该值理论值域为[0,0.5]，其公式可以表述为：

SIi=|1.5-Di|

(5)

式中：D为二维欧式空间的面积周长分维数，其值介于[1,2]，且D值越大，结构越复杂；但D值越接近1.5时景观结构越不稳定，即其处于类似布朗运动的随机运动状态[23]。D值的计算通常采用在ArcGIS中运用最小二乘法模拟面积周长线性函数得到[24]，本文用上文分析得到的景观分维数指数Fi替代各类景观的D值进行估算。

2.2.2 景观生态风险评价模型

(1) 指标标准化。由于各指数的性质不同，不同指数存在不同的量纲和量纲单位，直接加权计算会影响到评价结果的准确性。为了消除不同量纲的影响，需要对指标进行标准化处理。本文采取极差标准化对景观破碎度指数、景观分离度指数、景观分维数指数、植被覆盖度指数、景观脆弱度指数和结构稳定性指数进行标准化处理。

对于正向指标，需采用：

(6)

对于负向指标，需采用：

(7)

(2) 风险小区的划分。景观格局在空间上分布不均，具有空间异质性。通过对风险小区进行划分，可以在一定程度上弥补这种不足，使每个风险小区既具有区内同质性又具有区间异质性[19]。结合相关文献及研究区面积，选取2 km×2 km的格网作为评价单元不仅可以保证工作量适中，又可以对区内同质性和区间异质性有很好的区分，且评价单元精度合理准确，因此按此粒度将研究区分为973个风险小区。

(3) 景观生态风险评价。景观生态风险是区域内显性生态风险和隐性生态风险的结合，每个风险小区内其显性景观生态风险和隐性景观生态风险的计算公式分别为：

(8)

式中：ERIk显为第k个风险小区的显性景观生态风险指数；Aki为第k个风险小区中第i类景观的面积；Ak为第k个风险小区的总面积；PLk为第k个风险小区的植被覆盖度；a，b，c，d分别为相应的权重，由层次分析法得到；PLk是根据植被覆盖度PL和风险小区2 km×2 km的格网，在ArcGIS中进行分区统计得到，取平均值为每个风险小区的植被覆盖度。

(9)

式中：ERIk隐为第k个风险小区的隐性景观生态风险指数；e，f分别为相应的权重值，由层次分析法得到。

那么，每个风险小区的景观生态风险指数可以表示为：

ERIk=w1ERIk显+w2ERIk隐

(10)

式中：w1和w2分别为层次分析法得到显性指数的权重和隐性指数的权重。

表2 景观生态风险评价指标体系

2.2.3 空间相关性分析 空间的自相关性是景观格局分布的最大特征之一，空间相关性分析可以检验在空间分布上具有一定规律性的空间变量在各自不同的空间位置上的离散或聚集程度[26]。常用的测度空间自相关的模型有Moran′s I，Geary′s C，Join count等，在此本文基于Geoda软件，选取全局空间自相关指数Moran′s I和局部空间自相关指数LISA进行空间相关性分析。

(1) 全局空间自相关分析。全局自相关分析是对空间变量的某种属性在整个区域内的空间分布特征的整体描述，其计算公式如下[27]：

(11)

(2) 局部空间自相关分析。局部空间自相关是对局部小范围内的空间变量的某一属性值与周边单元相同属性值相似性和相关性的描述，也可以用来识别“热点区域”和检验数据的“异质性”，其公式表述为[28]:

(12)

式中：LISAi为局部空间自相关指数。由此绘制LISA聚类图，得到局部高—高、低—低等空间聚类现象。

3 结果与分析

3.1 景观生态风险时空变化分析

由图2可知，从空间分布上看，高风险区域主要集中在县域的西南部、中部以及东北部地区，该区域景观斑块破碎、分散，植被覆盖度不高，且未利用地较多，因此景观生态风险最高；较高风险区面积最大，几乎遍布于整个县域；而景观生态风险较低的区域主要分布在县域的县界周边，尤其是西北部的都山林场和东南部的祖山林场，该区域植被覆盖度高，景观斑块连续且较为完整，地类以林地为主。总体看来，青龙满族自治县西部的景观生态风险值高于县域东部地区。

图2 青龙满族自治县景观生态风险空间分布

由图3可知，2001—2013年，青龙满族自治县景观生态风险等级发生变化的区域面积较大，遍布全县各乡镇，整体上景观生态风险呈降低趋势，且东部降低较为明显。风险等级下降区域面积共121 882.90 hm2，且以降低1等别为主。风险等级上升区域面积为11 118.31 hm2，主要分布在马圈子镇、土门子镇、青龙镇及龙王庙乡；其中马圈子镇、土门子镇及青龙镇风险升高主要是开发建设活动使区域景观斑块破碎度增大导致显性景观生态风险升高明显，而龙王庙乡风险上升主要是荒草地退化所致。风险等级未发生变化的区域面积为217 556.70 hm2，占县域面积的62.06%。对比2001—2007年和2007—2013年风险等级变化，前者风险等级变化面积明显高于后者。2001—2007年风险等级上升区域面积为17 118.31 hm2，主要分布县域北部和东部区域；其中北部风险升高区也是由建设活动导致的景观破碎化加剧引起，而东部风险升高主要是人类干扰严重的低山丘陵区的荒草地退化引起。2007—2013年风险上升区域面积为21 208.74 hm2，比2001—2007年增长4 090.43 hm2，且主要分布在县域西北部和南部地区；其中西北部马圈子镇风险等级依旧升高，青龙镇风险等级升高主要由于建设用地增长和未利用地的退化导致植被覆盖度减小所致，而南部风险升高区则是由于水域面积减少引起。总体来看，景观生态风险等级的变化主要发生在人类活动干扰密集区，而都山林场、祖山林场以及县界周边区域由于植被覆盖度较高、高程较大而人类活动干扰较小景观生态风险往往较低。

由表3可知，2001—2013年，除林地和水域外各景观类型斑块数均呈先增加后减少趋势，但相比2001年的斑块数量均大幅度增加(未利用地除外)。各类型景观面积也均有波动，其中耕地、水域和未利用地面积均有不同程度的减少，且未利用地减少幅度最大，由2001年的108 387.09 hm2减至2013年的87 923.35 hm2，减少了20 463.74 hm2，减少的未利用地主要转化为耕地和建设用地，其次面积减少的速率依次是水域和耕地；建设用地、林地和园地面积均有不同程度的增加，且以林地面积增加最快，由2001年的141 220.89 hm2增至2013年的156 374.24 hm2，增加了15 153.35 hm2，其次按增速排列依次是园地和建设用地。由此可见，青龙满族自治县十分重视林地的保护，积极响应国家退耕还林的号召，林地的增加导致植被覆盖率升高，间接增加了生态系统的稳定性，降低了生态风险；但耕地斑块数的大量增加导致耕地破碎化严重，尽管由于地形坡度使耕地大规模集中连片较为困难，然而减少耕地破碎化、加强高标准农田建设仍应是工作重点。

图3 青龙满族自治县景观生态风险等级变化

景观类型年份斑块数面积/hm2破碎度指数分离度指数景观分维数脆弱度指数干扰度指数2001416430861.990.13490.01241.45130.19050.3614耕地20071809630638.880.59060.00591.50890.19050.598920131476930404.990.48570.00661.50070.19050.54502001308314003.740.22020.03081.45090.04760.4095建设用地20071002314678.200.68280.01681.49640.04760.64572013576315992.790.36030.02061.47550.04760.4814200115831141220.890.11210.00141.48120.09520.3527林地200717543150345.090.11670.00121.48590.09520.3559201317141156374.240.10960.00121.48650.09520.3525200113708116.290.16880.08181.43600.23810.3961水域200716457003.170.23490.08461.43920.23810.4307201320246951.870.29110.07751.43740.23810.4563200121105108387.090.19470.00161.49750.28570.3973未利用地20072449895235.210.25720.00171.50980.28570.431120131965587923.350.22350.00211.49970.28570.41232001933447967.930.19460.00491.47540.14290.3939园地20071637652657.380.31100.00351.49750.14290.456120131619852910.680.30610.00371.49720.14290.4536

由图4可知，2001—2013年，低风险区域面积略有增加，但基本变化不大；较低风险区域面积占比逐年升高，但最高也仅占5.00%左右；中等风险区域面积比例逐年升高，且在2001—2007年增速明显，由2001年的9.90%增至2007年的19.89%，增加了近1倍，2007—2013年增加幅度不明显，增幅在5.00%以内；较高风险等级在2001—2013年间面积占比均最大，达到40%以上，虽然每年增幅不大，但均呈增加趋势，其中2013年较高景观生态风险面积占比为47.47%；高风险区域面积占研究区总面积的比例呈减小趋势，由2001年的40.81%下降至2013年的24.98%，尤其是2001—2007年减少幅度较大，减少了11.84%，2007—2013年减少了3.99%。总体来看，青龙满族自治县的景观生态风险在这13年间有所下降，但整体风险仍较高。

图4 青龙满族自治县各年份景观生态风险面积占比

3.2 景观生态风险格局变化空间相关性分析

基于GeoDa空间分析软件，利用青龙满族自治县973个风险小区3期景观生态风险指数数据进行空间自相关分析，得到3期全局自相关Moran′sI指数，以及3期局部空间自相关LISA结果。结果显示，在显著性水平p<0.05下，青龙满族自治县2001年、2007年、2013年3年的Moran′sI分别是0.533 7，0.571 5，0.562 7，呈先增加后减少的趋势，尽管后来有所下降，但在这13年间整体上升了0.029，可以看出，青龙满族自治县的景观生态风险指数在空间上存在着较强烈的空间自相关性，相邻样地在空间上呈较高度的聚集分布。

从景观生态风险的局部空间自相关聚集来看，青龙满族自治县景观生态风险指数空间分布以高—高聚类和低—低聚类为主，高—低聚类和低—高聚类呈零星分布。2001—2013年13年间，景观生态风险指数高值聚集区主要分布在研究区中部、中西部以及东北部地区；其中2001—2007年，高值聚集区向西部、南部和东部扩展，呈分散演化，但东北部的高值聚集区则向中部靠拢，其原因是东北部地区在此期间加大了林业建设，森林覆盖率不断提高，苹果等经济林的种植也是景观生态风险指数下降的主要原因，而高风险值扩散区主要是区域人口的快速增加导致建设用地扩张、景观斑块的破碎度明显加大，人类活动剧烈；2007—2013年南部和东部的高值风险区减少，整体向中部地区收缩，主要是因为在此期间青龙满族自治县大力开展了土地整治，建设用地集中整理、复垦使该景观类型逐渐连片，景观斑块数量迅速缩减，且中南部的桃林口水库区域更加注重生态建设，加大水土涵养防护林建设，景观生态风险值下降。景观生态风险低值聚集区主要分布在县域周围，且在研究期间变化不明显，一般位于地形坡度稍高的山区，植被覆盖度高，人类活动干扰不明显，景观斑块完整，明显聚集地为都山林场和祖山林场。高—低聚类“奇异点”零星分布于县域的边界地带，且以西南部分布稍集中，原因是桃林口水库下游的林地和未利用地相间分布，未利用地植被覆盖度低，生态系统不稳定。

4 讨论与结论

4.1 讨 论

景观生态风险评价是复杂的系统性的评价，本文通过景观生态风险的表现对其风险的显隐性进行定义，进而选取评价因子进行评价，明确风险变化。未来的景观生态风险评价可以加上对风险本身的显隐性和潜在性的分析，如自然灾害预警等。生态功能区对维护区域生态安全和稳定有着重要意义，只有全面分析其风险，才能防患于未然，以期长久有效地发挥生态功能区的功能。

4.2 结 论

(1) 本文从景观生态风险的本质属性出发，根据风险的显隐性表现，从不同的角度选取评价指标，将显性景观生态风险和隐性景观生态风险相结合，得到研究区的整体景观生态风险，使风险分析更为全面，为景观生态风险评价研究提供了新的思路。不难发现，景观斑块的破碎化等显性风险是风险升高的直接原因，而隐性风险例如长久以往的人类活动干扰而导致的生态系统结构稳定性的破坏无疑是诱因。

(2) 研究结果显示，2001—2013年，青龙满族自治县虽然整体景观生态风险下降，但生态风险等级仍处于较高状态，且呈西部高东部低的特征。从空间上看，高风险区域主要集中在县域的西南部、中部以及东北部地区，主要由于人口的快速增长导致了人类活动干扰加剧，不同地类之间演变速度加快，建设用地面积的增加导致了景观斑块的破碎和形变，景观生态风险指数也相应升高。景观生态风险指数在空间上存在着较强烈的空间自相关性，且整体上呈高—高、低—低聚集的空间分布。因此，在未来的土地利用规划中，应加快对该区域的产业调整，重视生态环境建设，加大土地整治力度，使景观斑块趋于连续、完整，确保生态经济协调发展。

(3) 通过对风险变化区的地类变化分析可知，林地对于区域景观生态风险值的变化起着至关重要的作用，林地面积的增加可有效降低区域景观生态风险值；未利用地整体生态脆弱度高，将其开发为林地等其他地类可以明显降低景观脆弱度指数，但未利用地极快速地演变为建设用地在一定程度上反而会因为破坏自身生态系统结构而加剧区域景观生态风险指数的上升。因此，未利用地开发前应充分考虑开发前后的景观生态风险指数变化，控制未利用地合理适时适度地开发，得到更大的生态效益。