随着人类活动进一步加剧及社会经济快速发展，全球生态环境问题日益突出，对人们生活产生严重影响。在气候变化背景下，生态脆弱性评价研究日益引起相关学者关注^[1-2]，是可持续发展研究核心问题之一。为实现人与自然和谐共处，世界各国都将生态脆弱性研究作为研究重点，大力推进生态文明建设，并将脆弱性、韧性和社区减灾纳入减灾框架的关键词^[3]。对生态脆弱性的动态评价及其变化规律分析，可为生态恢复治理和环境保护提供重要的科学支持。国内外学者基于不同方法对不同区域的生态脆弱性展开了相关研究，研究区域主要涵盖行政区划^[4-5]、河流流域^[6-7]、山地高原等^[8]，但针对东北地区典型流域的研究较少；生态脆弱性评价方法众多，常见有主成分分析法^[9]、层次分析法^[10]、模糊综合评价法^[11]、熵权法等^[12]。各评价方法都有自身的不足，如层次分析法受主观思想影响较多，熵权法在数据量较大时，计算量大且权重差较小。近年来，随着生态脆弱性模型相关研究的不断深入，对应的评价方法由单一方法转向多种方法集成。比如主成分分析法和层次分析法的结合^[7]、层次分析法与熵权法的结合等^[13-14]，但这些方法的集成均侧重于二维点状数据间的优化，缺乏对空间数据的定量分析，而空间主成分分析能直观地将空间数据信息映射到二维空间^[15]。因此，本文借鉴已有研究成果，力图将层次分析法与空间主成分分析法相结合^[16](AHP-SPCA)，对生态脆弱性进行评价，实现多方法的优势互补。然而，在研究视角上，已有研究多集中于生态脆弱性评价、时空格局演变以及影响因素分析^[9-13]，对生态脆弱性未来发展的模拟预测较少，CA-Markov模型最早广泛应用于土地利用变化的预测中，但目前在其他领域的应用也取得了较好的效果^[17]。

松花江流域自然条件复杂，洪涝灾害频发。近年来，人类过度追求经济发展和人居环境改善，导致生态系统负担增重，水污染、湿地破坏、水土流失等生态问题加剧，为避免区域出现“返贫”现象，研究本区域的生态环境可持续性具有重要意义。鉴于此，本文以松花江流域为研究区，利用地学信息图谱的思维特征^[18]，综合评价2005，2010，2015和2020年4期松花江流域的生态脆弱性，揭示其时空变化规律；并结合CA-Markov模型，对2025年松花江流域生态脆弱性进行模拟预测，从而实现该地区生态脆弱状况的变化监测^[19-20]，以期为研究区内的生态保护和修复提供科学依据。

1 研究区概况及数据来源 1.1 研究区概况

松花江流域位于中国东北北部，位于东经119°52′—132°31′，北纬41°42′—51°38′，是中国主要的江河流域之一。流域面积广阔，约为5.57×10⁵ km²，覆盖黑、吉两省大部、内蒙古东部以及辽宁少部，共包括25个地级市。整个流域地形向东北方向倾斜，三面以群山环绕，河道可划分为三段。流域地处温带季风气候^[21]，四季较分明，温度变化较明显，年均降水量不超过500 mm，夏季降水以及春季积雪融水形成的春汛为河流主要来源。流域内水资源丰富，包含众多河流分支。

1.2 数据来源及处理

研究数据主要包括地形、地表、气象、植被和社会经济数据，年份为2005，2010，2015和2020年共4期数据。

(1) 气象数据、植被数据以及土地利用栅格数据均来源于中国科学院资源环境科学与数据中心，利用ArcGIS 10.2软件中的掩膜提取方法得到研究区所需数据；

(2) DEM来源于地理空间数据云(http://www.gscloud.cn/)的30 m分辨率的数字高程数据，坡度和地形起伏度通过ArcGIS 10.2软件提取；

(3) 景观多样性数据是在土地利用栅格数据基础上利用Fragstas 4.2软件进行多样性指数提取；

(4) 社会经济数据以市级为单位，均来源于《吉林省统计年鉴》《辽宁省统计年鉴》《黑龙江省统计年鉴》《内蒙古自治区统计年鉴》(2005，2010，2015和2020年)；利用反距离权重插值法(IDW)对社会经济数据进行空间化。为确保各指标在空间位置上有较好的重合性与协调一致性，所有指标均采用Krasovsky_1940_Albers坐标系统。

2 研究方法

通过分析松花江流域生态环境特征，构建因地制宜的生态脆弱性评价指标，选择层次分析法(AHP)、空间主成分分析法(SPCA)和最小相对信息熵对指标进行权重计算，使松花江流域生态脆弱性更具科学性。结合生态脆弱性图谱，进一步分析松花江流域生态脆弱性的时空特征，在此基础上，利用CA-Markov模型对松花江流域2025年生态脆弱性进行预测。

2.1 指标体系构建与标准化处理 2.1.1 指标体系构建

结合相关研究成果^[22-23]并充分考虑数据可获取性、典型性，本文从生态敏感性、生态恢复力和生态压力度3个维度构建生态脆弱性评价指标体系(表 1)。其中地形、地表条件是造成生态脆弱性的基础条件；气象、植被因子对生态具有重要的指示作用；社会因子表示人类活动的干扰强度。在此基础上，选取高程、坡度、地形起伏、年降水等11个正向指标和3个负向指标。

2.1.2 评价指标标准化

为进一步提高各指标数据收敛性，根据指标正负相关性完成数据标准化处理^[24]。正向指标包括高程、坡度、地形起伏、年降水量、人均GDP值、人口密度、建成区面积、人均拥有道路面积、农林牧渔业生产总值、在校学生数、医院床位数、社会因子主要体现在面临环境灾害时所施加的外部压力，指标值越高生态脆弱性压力越大；负向指标包括景观多样性、年均气温、归一化植被指数，其中气温越高越有利于植被的生长，景观多样性越大，归一化植被指数越高，生态系统越稳定^[25]。

2.2 生态脆弱性综合评价模型

在参考有关文献和考虑所选指标数据的基础上，采用空间主成分分析法(SPCA)和层次分析法(AHP)法分别计算主客观权重，利用最小相对信息熵确定综合指标权重。

(1) AHP是一种常用的主观评价方法^[26]，将各种元素划分为多个层次，通过对指标层的两两对比，构造判断矩阵以及一致性的检验最后得出合理的计算结果。

(2) SPCA是在主成分分析的基础上改变原有统计数据，利用空间数据、空间分析技术将多变量转换成几个不相关的综合指标，本文利用ArcGIS以及Matlab软件计算得到〔公式(1)—(3)〕，具体公式为：

$ A=\sqrt{\alpha_1 \beta_1}, \sqrt{\alpha_2 \beta_2} \cdots \sqrt{\alpha_n \beta_n} $

(1)

式中：α为各主成分特征值；β为主成分的正交化特征变量。

$ H_j=\sum\limits_{j=1}^m \lambda_{j k}^2 \quad(j=1, 2, 3 \cdots m) $

(2)

式中：H_j是第j个指标的公因子方差；λ_jk是第k个指标在主成分j上的特征值；m为主成分数(本研究m=4)。

$ W_{2 j}=H_j / \sum\limits_{j=1}^m H_j \quad(j=1, 2, 3 \cdots 14) $

(3)

式中：W_2j为各项指标的权重系数，最终计算出研究区2005—2020年生态脆弱性评价权重系数。

(3) 最小相对信息熵。运用最小相对信息熵得出优化后的权重，计算公式^[27]为。

$ W_j=\frac{\sqrt{W_{1 j} W_{2 j}}}{\sum\limits_\limits{j=1}^n \sqrt{W_{1 j} W_{2 j}}} $

(4)

式中：W_1j，W_2j分别为利用AHP和SPCA计算的各指标权重，优化后的权重为W_j。

(4) 生态脆弱性指数。利用ArcGIS中的栅格计算器计算，计算公式为：

$ \mathrm{EVI}=\sum\limits_{j=1}^m\left(K_i w_i\right) \quad(j=1, 2 \cdots m) $

(5)

式中：EVI为生态脆弱性指数；K_i为指标标准化结果；W_i为各指标权重；EVI值越大生态脆弱程度越高；反之生态越稳定。

2.3 脆弱性分等定级及综合指数

参照国内外研究成果^[28]，将松花江流域生态脆弱性指数划分为5级：潜在脆弱、微度脆弱、轻度脆弱、中度脆弱、重度脆弱，采用自然间断法选取同等级阈值的均值为最终分类标准^[29](表 2)。计算生态脆弱性综合指数(EVSI)，比照区域不同年份生态脆弱性整体差异^[30]，计算公式为：

$ \mathrm{EVSI}=\sum\limits_{i=1}^n P_i \frac{A_i}{S} $

(6)

式中：EVSI为生态脆弱性综合指数；P_i为脆弱性等级；A_i为第i个等级的栅格数量；n为等级总数；S为栅格总数。

2.4 生态脆弱性图谱构建

生态脆弱性图谱分为涨势图谱和落势图谱，利用ArcGIS 10.2对两期脆弱性评价结果进行叠加处理，如十位数大于个位数为转出即落势，十位数小于个位数为转入即涨势，从而进一步掌握流域整体变化趋势。本文以2005—2020年为例进行分析，公式^[31]如下：

$ Z=10 X+Y $

(7)

式中：Z为生态脆弱性图谱代码值；X, Y为两期生态脆弱性等级代码值(1—5), 因此, 在运算后Z表示2005年生态脆弱性分级类型转化为2020年生态脆弱性分级类型(如21指2005年的微度脆弱类型转变为2020年的潜在脆弱类型)。

2.5 CA-Markov模型

元胞自动机(Cellular automat，CA)加入了空间特征分析，具有较强空间模拟演变能力，在时间、空间状态都离散，模型为：

$ S_{i j}^{t+1}=f_q\left(S_{i j}^t\right) $

(8)

式中：S为第ij元胞的状态; t，t+1为时刻; f为转换函数; q为邻域。

马尔科夫模型是对一个系统离散状态的转移概率来模拟未来状况，与过去状态无关^[32]，通过两个不同时间的生态脆弱性转换概率矩阵来创建。

$ S_e=S_0 P_e $

(9)

式中：S_e为e个周期后的状态; S₀为初始时刻的状态; e为周期个数; P_e为转换概率矩阵。

本研究利用IDRISI软件，以2010和2015年数据为基础数据，预测2020年生态脆弱性。将2020年预测结果与生态脆弱性指数计算后的真实结果进行精度检验，验证CA-Markov模型的可用性，最后完成研究区2025年生态脆弱性的模拟预测。

3 结果与分析 3.1 生态脆弱性时空分布

结合空间分布特征(图 1)分析可知，流域生态脆弱性整体呈现东部向西部逐渐变低的空间变化，西部脆弱性比东部脆弱性相对较低。流域西部主要位于内蒙古地区，该地区人口较少，干扰程度较低且植被覆盖度较高，下渗量增加，地表径流量减少，生态压力较少；东部为黑龙江、吉林地区，相比之下经济、人类活动均高于内蒙古地区，但植被覆盖度较低。

松花江流域内主要以微度脆弱、轻度脆弱、中度脆弱为主，约占整个区域的75%。2005年松花江流域重度脆弱主要分布在西部和东部，2010年分布在西部和中部，其中西部重度脆弱面积减少且东部大面积减少，总体上有一定程度的改善；2010—2015年，生态环境继续呈好转的趋势，主要是西部中度和重度脆弱程度大面积减少；2015—2020年，西北部呼伦贝尔市和大兴安岭市脆弱等级急剧上升，生态环境不断恶化，主要是由于随着经济的发展建成区面积和人均拥有道路逐渐增多等原因。潜在脆弱区与中度脆弱性区都约占整个区域的10%，潜在脆弱区大多集中于白城市、松原市、通辽市附近，其地形起伏较低、植被覆盖较高，重度脆弱区大多集中在伊春市、抚顺市、白山市及通化市，区域内植被覆盖较低人口密度较大。

利用公式(3)计算松花江流域2005—2020年EVSI值，分别为3.105 2，2.898 9，2.840 8，3.182 3。从时间尺度上看出，2005—2015年EVSI整体呈现下降趋势，表明2005—2015年松花江流域生态环境整体呈现好转的发展趋势，潜在脆弱性比例由0.84%升到17.84%，其中年均温由4.53 ℃升高到4.81 ℃，年降水量由608.3 mm增长到616.7 mm，这两项指标的升高对生态脆弱性变化具有促进作用。随着近几年各市县生产总值以及农林牧渔业生产总值的增高，生态保护意识不断加强，合理利用资源和生态文明建设仍是可持续发展的两个重要手段。而2020年EVSI整体呈上升趋势，向着恶化的方向发展，重度比例由10.36%增长到13.83%，2015—2020年人口密度由119.9人/km²增长到130.5人/km²，建成区面积由112.6 km²升高到129.4 km²，人为干扰严重，人类对自然的不合理开发与利用正是生态恶化的主导因素。

3.2 生态脆弱性变化趋势分析

通过对EVI均值定级，对比不同等级脆弱区松花江流域2005—2020年各时期的栅格数量及比例，从而进一步了解该区域生态脆弱性状况的结果差异和变化特征。

从表 3可知，各等级栅格比例在结构分布上其比例大小从大到小依次约为轻度、微度、中度、潜在、重度。2005年轻度及以下脆弱性程度总栅格数占全部栅格数的66.3%，中度和重度所占比例分别为23.78%和10%，不足总栅格数的2/5，表明该区域生态脆弱性整体处于中等水平。2010和2015年轻度及以下脆弱区百分比均高于65%，只有2020年稍低为60.43%，但2005—2020年重度脆弱区百分比均在10%左右，所占比例较少，进一步证明松花江流域仍处于中等水平。

为掌握流域脆弱性整体变化趋势，利用变化图谱进一步了解其等级升高和降低情况。从表 4可知，2005—2020年生态脆弱性变化图谱共生成25类图谱，其中变化的图谱有20类。变化栅格占总栅格数的61.15%，等级上升和下降的栅格分别为31 379，28 869个，占变化总栅格的51.80%，48.20%。

根据图 2可知，2005—2020年，松花江流域生态脆弱性轻度脆弱涨势最强，转入栅格数为48 257个，占变化栅格的20.28%，主要位于黑河市和呼伦贝尔市。潜在脆弱转入栅格数最少为1 396个仅占总栅格的0.59%，主要位于大庆市、齐齐哈尔市、松原市和绥化市，呈零星式分布。2005—2020年中度落势最强，转出栅格数为33 683个，占变化栅格总数的14.15%，主要集中于东部的长春市、哈尔滨市、牡丹江市以及延边朝鲜自治州附近。微度、轻度落势较弱，分别占变化栅格总数的9.33%和9.54%，位于白城市、通辽市和大庆市附近。通过对比各等级栅格变化情况，涨势较强大部分位于西部，落势较强大部分位于东部。

3.3 生态脆弱性预测结果与分析

2020年生态脆弱性指数计算后的真实结果与模拟成果的kappa系数检验值为0.81，系数大于75%，模型适用于该预测研究。对该地区2025年进行模拟预测(图 3)，从而对该区域未来生态脆弱性状况进行评价。

2025年预测结果表明，潜在、微度、轻度、中度、重度各类别的栅格比例分别为1.57%，26.17%，32.17%，24.27%和15.82%，微度和轻度仍比例最大，潜在和重度比例最小，与往期数据趋势吻合。相比2020年潜在轻度和中度比例呈减少的趋势，微度和重度呈增加的趋势。脆弱性较高的区域主要位于佳木斯市、双鸭山市、兴安盟、白城市和松原市；脆弱性较低的区域为黑河市、呼伦贝尔市、大庆市和通辽市。呈现这一趋势与自然环境有一定的关联，但人类活动及此区域政策的实施同样有着较大影响。在生态脆弱性图谱中，2025年轻度脆弱区涨势最强，占转入栅格总数的27.32%，潜在脆弱区涨势最弱，仅占9.17%，2020—2025年涨势图谱的变化情况与2005—2020年变化情况大致相同；微度脆弱区落势最强，占转出栅格总面积的29.8%，但中度和重度涨势所占比例仍高于其落势所占比例。同时，经预测该区域2025年EVSI值为3.265 9，仍呈上升的趋势(图 4)。从而表明松花江流域生态脆弱性仍向着恶化的趋势发展，预示着该区域政府应大力加强植树造林、退耕还林等项目，高效及时改善生态环境。

4 讨论与结论 4.1 讨论

本文以松花江流域为研究区，利用地学信息图谱分析研究区内2005—2020年4期生态脆弱性时空演变规律，结合CA-Markov模型对2025年生态脆弱性进行模拟预测，为资源有效利用和生态系统可持续发展提供理论参考。研究结果显示，在时间上，生态脆弱性综合指数呈现上升趋势；在空间上，研究区西部生态脆弱性等级涨势较强，东部落势较强。在政策上^[33]，2020年为我国脱贫攻坚决胜年，我国西部贫困县数量最多，内蒙古自治区正位于我国西北部。因此，该区域应注重对生态系统的积极响应，在加大扶贫力度的同时还需因地制宜地开展生态保护工作，实现生态环境可持续发展。相较于已有研究^[16-17]，本文采用AHP-SPCA熵权模型、地学信息图谱以及CA-Markov模型相结合的方法，能够更加客观地反映区域生态脆弱性时空演变特征以及生态脆弱性的模拟预测。空间主成分分析(SPCA)运用栅格插值数据计算其正交化的特征向量以及公因子方差得到各指标权重，不仅提高了计算效率也提高了准确性，同时结合层次分析法(AHP)并利用最小相对信息熵^[34]对权重进行综合优化，达到主客观结合。生态脆弱性图谱^[35]不仅从数值上分析其涨落趋势，更能从空间上对比不同区域的差异及同一区域变化情况。CA与Markov模型的结合，弥补了两种方法的局限性，最终既能对空间格局的演变进行分析，又能对未来的发展状况进行模拟预测。虽然本文借鉴国内外研究基础^{[7, 10]}，并考虑实际情况构建了市域生态脆弱性评价体系，但由于系统复杂性和生态脆弱性受多重因素干扰，指标体系的构建还需进一步优化。

4.2 结论

(1) 结合RS遥感技术和GIS空间分析技术，从主观和客观双重角度出发，利用层次分析法、空间主成分分析法和最小相对信息熵对松花江流域生态脆弱性进行评价，以避免主观权重的盲目性和客观权重的片面性，从而提高结果的准确性。

(2) 在空间尺度上，松花江流域西部大部分区域的生态脆弱性等级低于东部的地区，但研究期间西部生态脆弱性等级出现上升、脆弱性加重的变化趋势，主要原因可能为部分地区自然环境复杂，某些特殊区域人工无法干预实施。以2015—2020年为例，在行政区尺度上，大兴安岭市、呼伦贝尔市生态脆弱性变化最大，各政府应加大力度并采取有效的措施。在时间尺度上，松花江流域2005—2020年EVSI值均在2.5~3.5之间，整体均处于中等脆弱水平。主要以微度和轻度脆弱为主，所占比例超过研究区的3/5。

(3) 生态脆弱性涨落势图谱中，轻度和中度转出数比例最大，微度和轻度转入数比例最大，中度和重度转入数比例均较高，均超出20%，总体而言生态脆弱性图谱能较好体现空间上的变化。利用CA-Markov模型对2025年进行预测可知微度和轻度仍比例最大，潜在和重度比例最小，EVSI值仍呈增长的趋势，整体生态脆弱性朝着恶化的方向发展。