松嫩高平原黑土区土地利用／覆被变化及其景观生态风险分析

孙丽娜，宋 戈

（东北农业大学 资源与环境学院，哈尔滨150030）

近年来，土地利用／覆被变化在全球环境变化和可持续发展研究中占有重要的地位，绝大部分自然生态系统的空间已被人类各种土地利用所改变，土地利用／覆被变化在很大程度上记录了人与自然相互作用的过程[1]。景观生态风险分析是在景观基础上，描述人类活动或自然灾害对区域内生态系统结构、功能等产生不利生态效应的可能性和危害程度的过程[2]。由于自然条件和人类活动干扰程度的差异，不同区域生态功能与土地利用变化过程是不一样的，土地利用变化必然导致景观生态风险的变化，进而影响区域土地生态结构和功能。因此，探索区域土地利用变化及其产生的生态风险效应对于实现区域生态环境保护和恢复具有重要意义。

目前，国内外学者相关的研究主要有：一是对土地利用／覆被变化的研究，代表学者有李秀彬、葛全胜、史培军、摆万奇等人，分别对不同区域的土地利用变化、对区域的影响及动态演变模拟进行了研究[3－6]；二是对生态风险的研究，代表学者有Preston B、陈鹏、高永年、卢宏玮等，分别对生态风险评价及生态功能分区等方面进行研究[7－12]；三是土地利用变化及其景观生态风险研究，张银辉、王娟、高永年等人的研究区多集中于流域地区[13－16]，肖杨、李景刚、谢花林等人的研究则是在分析土地利用变化的基础上，对景观生态风险进行空间分析[16－20]。本文突破研究区限制，以松嫩高平原黑土区巴彦县为研究区，该区域耕地面积大、质量好，是国家重要的商品粮生产基地，研究该区域的土地利用变化及其对景观生态风险的影响，对于保护耕地生态环境，确保粮食安全以及实现土地可持续利用意义重大。

1 研究区概况

松嫩高平原黑土区位于黑龙江省西部松嫩平原腹地，属于由小兴安岭余脉向松嫩平原过渡地带，此区域具有土壤肥沃、腐殖质深厚、有机质含量高、团粒结构好等特点，成为重要的农业土壤，是宝贵的土地资源，年均粮食商品率高达70%以上。巴彦县是松嫩高平原黑土区内具有很强代表性的区域，辖10镇、8乡、116个行政村，总人口为70.0万人，其中农业人口56万人，总土地面积31.33万hm2，耕地占土地总面积的72.5%，是国家重点商品粮生产县。全县地势东高西低、北岗南平、中部多丘陵，素有“三山一水六分田”之称，地理形状为北宽南窄的楔形。该区域属于中温带大陆性季风气候，夏季温暖多雨，冬季寒冷干燥，雨热同季，并以黑土和草甸土为主，是典型波状缓坡漫岗平原耕地利用的地貌类型。

2 数据与方法

2.1 数据获取

研究所需的资料为巴彦县1991年、2006年和2009年三个时期的土地利用数据，数据来源于对1991年、2006年和2009年3个年份的TM／ETM＋遥感影像的解译。土地利用分类为耕地、林地、园地、牧草地、水域、建设用地及其它用地7种景观斑块类型，土地利用数据用COVERAG格式存放。

2.2 评价单元的划分

根据研究区范围和采样的工作量，将研究区划分为4km×4km的正方形作为评价单元，采样方式为等间距，共有244个采样区，计算每一样区的景观损失指数，以此作为样地中心点的生态风险值。

2.3 研究方法

在GIS技术的支持下，以景观格局指数为评价指标分析松嫩高平原黑土区巴彦县土地利用变化及其对区域景观生态风险的影响。

2.3.1 土地利用变化分析方法 土地利用变化的实质是区域内不同用地之间的相互转换，本文在Arc－GIS平台下分别将1991年、2006年和2006年、2009年的土地利用矢量数据进行空间叠加分析，得出研究时段内各类用地之间的转换数量及空间分布。并引用各地类的转化率[11]（CRijk）来研究各地类之间的相互转换趋势和速率，计算公式为：

式中：CRijk——在时间t内，地类i转变为地类j占地类j变化面积的比例，即地类i对地类j面积变化的贡献率；Zit－Zjt——在时间t内地类i转变为地类j的面积；ΔZjt——在t时间段内景观类型j的变化面积（增加或减少）；i，j——初期、末期的地类，分别为耕地、林地、园地、牧草地、水域、建设用地和其它用地。

2.3.2 区域景观生态风险分析方法 为建立景观结构与区域生态风险之间的联系，本文引用生态风险指数[13]，描述样地内综合生态风险的相对大小，通过采样的方法将景观空间结构转化为空间化的生态风险变量，计算公式为：

式中：ERi——生态风险指数；n——景观分组类型的数量；Ai——区域内第i类景观组分的面积；A——景观的总面积；Ri——第i中景观组分所反映的景观损失指数。

式中：Ri——景观损失指数，遭遇干扰时各类型景观所受到的生态损失的差别，即其自然属性损失的程度，是某一景观类型的景观干扰指数Ui和景观敏感度指数Si的综合；Ui——景观干扰指数，在自然、人为干扰下，景观一般由单一、连续和均质的整体向复杂的、不连续和异质的斑块镶嵌体变化[13，21]；Si——景观敏感度指数，不同生态系统的易损性，与其在景观自然演替过程中所处的阶段有关，一般情况下，处于初级演替阶段、食物链结构简单、生物多样性指数小的生态系统较为脆弱，本文借鉴他人研究成果[8，14，21]，并结合研究区实际将该区7类景观类型敏感度指数排序由高到低依次为：1建设用地、2林地、3耕地、4牧草地、5园地、6水域、7其它用地，进行归一化处理后得到各自的脆弱度指数；Ci——景观破碎度；Fi——景观分离度；Di——景观优势度；Ni——景观的斑块数；Ai——斑块面积；A——景观总面积；m——景观中斑块类型的总数；Pi——各类景观的面积比；a，b，c——别为破碎度、分离度和优势度的权重，借鉴前人的相关研究成果并结合研究区的实际情况，a，b，c的权重值为0.5，0.3，0.2，其中其它用地的权重为0.3，0.2，0.5；i——耕地、林地、园地、牧草地、水域、建设用地、其它用地7种景观类型。

2.4 空间分析方法

地统计学是一系列检测、模拟和估计变量在空间上的相关关系和格局的统计方法。区域生态风险指数本身是一种空间变量，空间变化特征具有结构性和随机性，本文采用地统计学中的变异函数方法，借助半方差函数进行区域生态风险程度的空间分析，计算公式如下：

式中：r（h）——样本距为h的半方差；h——样本间隔距离；n（h）——抽样间距为h时的样点对总数；Z——某一系统属性的随机变量；x——空间位置；Z（xi），Z（xi＋h）——变量在xi和xi＋h点的取值。以半方差r（h）为纵坐标，样本间距h为横坐标作图即半方差图。运用ArcGIS的空间分析和地统计功能，通过求和、采样、普通Kriging空间插值，以及对所得的实际半方差图用球状模型拟合，从而得到生态风险指数空间分布图。

3 结果与分析

3.1 研究区土地利用变化

3.1.1 研究区土地利用结构及特点 通过对巴彦县1991年、2006年及2009年土地利用数据的统计及分析，由表1可以看出，研究区在这18a间，耕地、牧草地、水域和建设用地面积都有所增长，其中耕地面积增加比较明显，占土地比重由1991年的70.53%增加到2009年的74.83%，其次是水域、牧草地和建设用地；其它地类面积均不同程度地减少，其中林地面积减少后又出现增加，但总数2009年的较1991年的仍是减少，其它用地面积减少明显，占土地比重由1991年的5.58%减少到2009年的0.08%。上述变化表明，巴彦县作为农业大县，对耕地的保护力度较大，且其耕地数量的增加主要来源于林地和其它用地的转化。

表1 不同年份巴彦县各类用地比重 %

3.1.2 研究区各地类时空转化 运用ArcGIS的空间叠加分析功能，将1991年、2006年和2009年土地利用图相互叠加，得出各地类面积之间的相互转化数量关系，并运用转化率定量研究各地类之间的转化速度与方向，揭示各地类相互转化规律及程度。

1991—2006年，耕地、水域和建设用地面积增加，林地、园地、牧草地和其它用地的面积减少（表2）。由公式（1）可以算出各地类之间的转化率，其中，林地、建设用地向耕地的转入率较高，分别为3.42%，2.51%，耕地、牧草地和其它用地向水域的转入率较高，分别为9.13%，10.02%，11.58%，耕地向建设用地的转入率较高，为11.30%；林地向耕地的转出率为3.40%，园地向耕地、林地的转出率分别为35.64%和6.70%，牧草地向耕地和水域的转出率分别为7.34%和2.35%，其它用地向耕地、林地、水域和建设用地的转化率分别为11.87%，3.00%，2.93%，1.36%。

2006—2009年，耕地、林地、牧草地、水域面积增加，园地、建设用地和其它用地面积减少（表3）。同样可以算出各地类之间的转化率，其中，林地、牧草地和其它用地向耕地的转入率分别为1.28%，1.49%，3.01%，耕地和其它用地向林地的转入率分别为11.19%，2.72%，耕地和其它用地向牧草地的转入率分别为48.58%，25.07%，耕地、其他用地向水域的转入率分别为33.50%，17.64%，耕地向建设用地的转入率为13.71%；园地向耕地、林地的转出率分别为74.37%，14.04%，其它用地向耕地、林地和建设用地的转出率分别为46.50%，8.56%，6.77%。

表2 1991－2006年研究区土地利用转移矩阵 hm2

表3 2006－2009年研究区土地利用转移矩阵 hm2

通过对研究区18a年间的土地利用类型转化情况进行分析可以得知：（1）地类之间的转化比较复杂，既有其它地类向某一地类的转入，又有各地类向其它地类的转出；（2）在各地类相互转化关系中，其它地类向耕地的转化率较高，表明研究区对耕地保护的力度加大，并取得一定成效；（3）林地、其它用地面积的减少对耕地增加的贡献率较大，说明研究区产业结构调整向着发展粮食种植的方向进行，这与巴彦县作为我国产粮大县的职能定位相符合；（4）建设用地面积先增加后减少，18a间共增加了70.31hm2，且主要来源于耕地、林地和其它用地的转化，建设用地面积的增加是人类社会发展的必然结果。

3.2 研究区景观生态风险评价

3.2.1 研究区景观生态风险的时序变化 根据公式（2）—（7），对研究区1991年、2006年及2009年244个评价单元的生态风险指数进行计算，从结果（表4）可以看出各评价单元的生态风险指数变化趋势为先减后增、持续减少、持续增加和先增后减，这四种趋势所占的面积比例分别为20.08%，10.24%，37.30%，32.38%。18a间研究区有68.03%的评价单元的生态风险指数增大，生态安全降低。

为了便于比较研究区生态风险的大小，运用SPSS软件绘制生态风险指数分布图，根据曲线斜率的突变点（拐点）确定等级数目和临界点。巴彦县的生态风险指数划分为五个级别，即高生态风险区（Ⅰ：ER≥0.2730）、中生态风险区（Ⅱ：0.2417≤ER＜0.2730）、低生态风险区（Ⅲ：0.2000≤ER＜0.2417）、弱生态风险区（Ⅳ：0.1273≤ER＜0.2000）和无生态风险区（Ⅴ：ER＜0.1273）。本文生态风险级别的划分是针对该研究区的相对标准，便于比较研究区内的风险大小情况，并非绝对标准，不能与其它研究区相比较。

表4 1991年、2006年及2009年研究区生态风险评价结果

对研究区244个评价单元的生态风险指数各级所占面积进行统计（图1），可以看出，1991—2009年间，各级生态风险变化趋势不一，无生态风险区和弱生态风险区面积逐年减少，平均年下降率分别为2.99%和0.76%，低生态风险区面积增加，平均年增长率为0.34%，中生态风险区面积先减少后增加，平均年增长率为1.07%，高生态风险区面积增加后又出现略微减少，平均年增长率为2.34%。

图1 不同年份各级生态风险面积所占比例

3.2.2 研究区景观生态风险的空间分异 本文运用地统计学方法对研究区生态风险空间分异进行研究。首先对研究区生态风险指数进行常规性统计分析，得知生态风险指数的整体特征，然后对数据分布类型进行正态分布检验，得知数据的集中性和变异性情况，再运用地统计软件GS＋3.0对数据进行半方差分析，反映不同距离观测值之间的变化趋势，同时根据数据的不同类型的半方差模型进行拟合，根据拟合度和残差的大小来判断选取最佳理论模型，分析发现研究区的理论模型符合指数模型，其基台值为0.585，基台值明显说明生态风险指数具有明显的空间依赖性和空间结构，变程为1.089km，块金值与基台值的比值为0.502，表明变量的空间相关性中等，即生态风险的空间变异除了受内在结构性因子影响之外，还受外在因素的影响，如人为干扰、社会经济因素等。通过分析发现研究区生态风险指数在一定的空间范围内具有空间相关性，其理论函数拟合的决定系数经统计检验达到显著水平，半方差函数模型的选取符合要求。

应用ArcGIS 9.2的地统计分析集模块（Geostatistical Analyst）中的Kriging（克里格插值）方法对1991年、2006年和2009年采样数据进行空间最优内插，并用研究区区域边界进行分割，得到三个时期巴彦县生态风险指数的空间分布情况，结果如图2所示。1991—2009年巴彦县生态风险逐年增加，通过分析发现，研究区生态风险逐年加大，生态风险分布主要受地形和地类的影响较大。巴彦县东部地势较高，主要为林地，其生态风险较小；中部生态风险较高区域主要与建设用地分布有关；周边水域附近生态风险相对较小。具体分布情况如下：1991年生态风险整体较低，相对较高的区域主要分布在龙庙镇、巴彦镇、兴隆镇、德祥乡、天增镇和山后乡，且均处于中生态风险区，其它乡镇均属于低生态风险区、弱生态风险区，其中黑山镇部分区域属于无生态风险区；2006年出现高生态风险区，主要分布在德祥乡、天增镇和龙庙镇，黑山镇、富江乡、松花江乡、红光乡及西集镇的部分区域，属于弱生态风险区，其它乡镇的大部分区域均属于中生态风险区和低生态风险区；2009年高生态风险区面积继续增加，主要分布在德祥乡、天增镇、龙庙镇和巴彦镇，黑山镇、富江乡及松花江乡的部分区域属于弱生态风险区，且生态风险区在各乡镇所占的面积较少，其它乡镇的大部分区域均属于中生态风险区和低生态风险区。

图2 不同年份生态风险克吕格插值空间分布

4 结论与讨论

（1）18a间，巴彦县土地利用变化表现为耕地、牧草地、水域和建设用地面积有所增长，其它地类面积均不同程度减少。上述变化表明，巴彦县作为农业大县，对耕地的保护力度较大，且其耕地数量的增加主要来源于林地和其它用地。

（2）通过土地利用类型相互转化分析可以看出，18a间，巴彦县各地类之间的转化比较复杂，既有某一地类向其它地类的转入，又有各地类向其它地类的转出。其中，其它地类向耕地的转化率较高，林地、其它用地面积的减少对耕地增加的贡献率较大，建设用地面积先增加后减少，且主要来源于耕地、林地和其它用地的转化。导致以上结果的原因主要是研究区对耕地保护的力度加大，产业结构调整向着发展粮食种植的方向进行，建设用地面积的增加是人类社会发展的必然结果。

（3）1991年、2006年及2009年，研究区生态风险指数计算结果表明：各评价单元的生态风险指数变化表现为先减后增、持续减少、持续增加和先增后减，有68.03%的评价单元的生态风险指数增大，说明全县整体生态安全降低。通过生态风险分区可以看出，18 a间，无生态风险区和弱生态风险区面积逐年减少，低生态风险区面积增加，中生态风险区面积先减少后增加，高生态风险区面积增加后又出现略微减少。

（4）1991年、2006年及2009年巴彦县生态风险指数克吕格插值结果表明，研究区生态风险逐年加大，生态风险分布主要受地形和地类的影响较大。巴彦县东部生态风险较小，中部生态风险较高，周边水域附近生态风险相对较小。1991年，生态风险整体较低，相对较高的区域主要分布在龙庙镇、巴彦镇等6个乡镇，且均处于中生态风险区，其它乡镇均属于低生态风险区、弱生态风险区，其中黑山镇部分区域属于无生态风险区；2006年，出现高生态风险区，主要分布在德祥乡、天增镇等8个乡镇的部分区域，属于弱生态风险区，其它乡镇的大部分区域均属于中生态风险区和低生态风险区；2009年，高生态风险区面积继续增加，主要分布在德祥乡、天增镇等4个乡镇，黑山镇、富江乡及松花江乡的部分区域属于弱生态风险区，且生态风险区在各乡镇所占的面积较少，其它乡镇的大部分区域均属于中生态风险区和低生态风险区。

本文以景观格局指数作为巴彦县相对生态风险的评价指标，没有考虑社会、经济及环境等其它因素对生态风险的影响，但由于土地利用景观格局的改变必定会引起区域生态功能的变化，故通过研究不同土地利用景观类型之间的转移变化特征来识别区域生态环境的变化趋势及其内在因素是可行的。此外，本文的生态风险分级是为了便于比较研究区内生态风险大小的不同，是相对标准，与其它区域相比较不具有可比性。本文的研究结果可以为区域生态环境管理提供数量化的决策依据和理论支持，根据景观生态风险的高低程度，进行相应的生态建设和环境保护，从而提高研究区的土地生产功能和环境功能；同时，对于土地生态风险区域的生态建设也不能忽视，确保低生态风险区不向高生态风险区演化是实现区域生态环境良好的基本保障，最终实现巴彦县生态环境和社会经济建设的协调和可持续发展。

[1]吴次芳，宋戈.土地利用学[M].北京：科学出版社，2009：420－431.

[2]Barrell S M，Gardner R H，O′Neill R V.Ecological Risk Estimation[M].Boca Raton：Lew Publishers，1992.

[3]李秀彬.全球环境变化的核心领域：土地利用／覆被变化的国际研究动向[J].地理学报，1996，51（6）：553－557.

[4]葛全胜.20世纪中国土地利用变化研究[J].地理学报，2000，55（6）：698－706.

[5]史培军，袁艺，陈晋.深圳市土地利用变化对流域径流的影响[J].生态学报，2001，21（7）：1041－1050.

[6]摆万奇，张永民，阎建忠，等.大渡河上游地区土地利用动态模拟分析[J].地理研究，2005，24（2）：206－212.

[7]Preston B，Shackelford J.Multiple stressor effects on benthic biodiversity of Chesapeake Bay：implication for ecological risk assessment[J].Ecotoxicology，2002，11（2）：85－99.

[8]陈鹏，潘晓玲.干旱区内陆河流域区域景观生态风险分析[J].生态学杂志，2003，22（4）：116－120.

[9]高永年，高俊峰.太湖流域水生态功能分区[J].地理研究，2010，29（1）：111－117.

[10]卢宏玮，曾光明，谢更新，等.洞庭湖流域区域生态风险评价[J].生态学报，2003，23（12）：2520－2530.

[11]Wu L Y，HuangY X.Landscape ecological risk assessment of Dongshan Island[J].Journal of Oceanography in TaiwanStrait，2005，24（1）：35－42.

[12]Kapustka L A，Galbraith H，Luxon M，et al.Using landscape ecology to focus ecological risk assessment and guide risk management decision－making[J].Toxi－cology and Industrial Health，2001，17（5／10）：236－246.

[13]张银辉，罗毅，刘纪远，等.内蒙古河套灌区土地利用变化及其景观生态效应[J].资源科学，2005，27（2）：141－145.

[14]王娟，崔保山，刘杰，等.云南澜沧江流域土地利用及其变化对景观生态风险的影响[J].环境科学学报，2008，28（2）：269－277.

[15]高永年，高俊峰，许妍.太湖流域水生态功能区土地利用变化的景观生态风险效应[J].自然资源学报，2010，25（7）：1088－1096.

[16]Abell R，Thieme M L，Revenga C，et al.Freshwater eco－regions of the world：A new map of bio－geographic units for fresh－water biodiversity conservation[J].Bio－Science，2008，58（5）：403－414.

[17]于兴修，杨桂山，李恒鹏.典型流域土地利用／覆被变化及其景观生态效应：以浙江省苕溪流域为例[J].自然资源学报，2003，18（1）：13－19.

[18]肖杨，毛显强.区域景观生态风险空间分析[J].中国环境科学，2006，26（5）：623－626.

[19]李景刚，何春阳，李晓兵.快速城市化地区自然／半自然景观空间生态风险评价研究：以北京为例[J].自然资源学报，2008，23（1）：33－47.

[20]谢花林.基于景观结构和空间统计学的区域生态风险分析[J].生态学报，2008，28（10）：5020－5026.

[21]李晓燕，张树文.基于景观结构的吉林西部生态安全动态分析[J].干旱区研究，2005，22（1）：57－62.