生态安全是国家安全的重要组成部分，是区域与国家经济安全与社会安定的生态环境基础和支撑^[1]。近年来，随着生态文明建设思想的不断深入，生态安全已成为区域可持续发展的重要关注点和研究热点^[2]。随着城市化进程的加快，区域生态环境和景观格局正经历着前所未有的变化，引发了诸如生物多样性减少、景观破碎化、土地荒漠化等一系列生态环境问题，深刻影响着区域生态安全^{[3, 4]}。

漓江流域作为桂北重要生态屏障、南岭地区生物多样性保护重地、国家重要生态功能区以及喀斯特世界自然遗产保护地，其生态战略地位十分突出^[5]。然而，漓江流域是一个相对封闭且生态环境极为敏感的喀斯特区域。近年来，快速城市化建设以及高强度旅游开发等因素对原本脆弱的喀斯特景观生态系统造成了巨大压力。胡金龙等^[6]研究表明，漓江流域高生态风险区域正在迅速扩张；此外，由于研究的相对滞后和认知上的偏差，漓江流域旅游业开发与保护未能有效统一^[7]，这些都对漓江流域生态安全造成了威胁。因此，深入探讨漓江流域生态安全对维持其区域可持续发展具有重要意义。

当前，生态安全的评价方法主要包括生态足迹法^[8]、压力-状态-响应(PSR)模型^[9]、最小累积阻力模型^[10]和景观指数模型^[11]。其中，PSR模型发展最早，理论基础最为成熟，但在生态安全空间格局的评价上存在一定的不足^[12]。例如，采用PSR模型评估某一区域的生态安全时，由于PSR模型采用静态面板数据进行评价^[13]，无法充分体现区域内部的空间差异，导致评价结果的准确性下降。最小累积阻力模型可以较好地弥补这一不足，然而其主要基于高程、坡度等自然因子，对人类活动的影响考虑不够充分^[14]。景观作为自然过程与人类活动相互作用的重要载体，其格局变化直接影响区域生态系统的演变过程和功能^[15]。景观指数模型充分利用景观指数能够有效反映景观的功能、动态和空间异质性的特点，从景观角度研究区域生态安全(这种方法评价的生态安全被称为景观生态安全)，其既能体现人与环境的相互作用，又能进行良好的空间展示^[16]，已成为当前生态安全研究的一种主流方法。

生态安全在时间和空间上具有差异性。通过分析其时空变化，可了解生态安全格局的形成过程和演变规律^[17]。因此，区域生态安全时空演化特征分析已成为当前研究的热点^{[18, 19]}，并取得丰硕成果，如张湘等^[20]对长江经济带景观生态安全的研究以及施馨雨等^[21]对云南省景观生态安全的研究，等等。然而，上述研究都是从时空演变的角度探讨景观生态安全的变化过程及变异特征，而且基于现有数据的评价只能揭示景观生态安全格局的过去和现在，无法有效预测和展示未来的状况，制约了相关研究成果的进一步应用。因此，通过将景观指数模型与景观模拟技术结合^[22]，实现对未来景观生态安全格局的评价，对区域生态系统管理与建设具有重大指导意义。

此外，景观生态安全格局是在多种影响因素的相互作用下形成的，这些影响因素在时间和空间上存在差异^{[23, 24]}。漓江流域是生态战略地位极其重要而生态环境又极其敏感的区域，但在人们对漓江流域生态安全评价的过程中，却极少关注其变化影响因素^[25]，这限制了人们对漓江流域景观生态安全的全面认识。与其他流域相比，漓江流域喀斯特地貌高度发育，旅游业是其最重要的经济支柱，区域内农业生产较为发达，这些特殊性使得漓江流域生态安全的影响因素与其他流域有所区别，需进一步研究。鉴于此，本研究运用景观指数模型对漓江流域过去和未来的景观生态安全格局变化进行分析，并采用随机森林(RF)模型明晰漓江流域景观生态安全格局的影响因素，以期为漓江流域及相似区域的生态环境保护和可持续发展规划等提供科学依据。

1 材料与方法 1.1 研究区概况

漓江流域位于广西东北部的桂林市境内，其地理位置为110°5′-110°44′E、24°38′-25°56′N^[26]。漓江发源于兴安县华江乡猫儿山东北面老山界南侧，自北向南流经灵川县、桂林市辖区及阳朔县，在平乐县平乐镇北与恭城河汇合，属于珠江水系西江干流桂江支流的上游河段。流域北部为猫儿山，地势较高；东部连接海洋山脉，主要由碎屑岩构成的中低山地貌组成；中部地势相对平坦，以低山、丘陵及河谷平原为主，主要由碳酸盐岩构成，发育有峰林平原、峰丛洼地、峰丛谷地等多种地貌类型；西部与架桥岭山脉相接，为碎屑岩构成的中低山地貌；南部地势逐渐降低，河流纵横，以低山、丘陵和河谷平原为主。流域内多年平均气温为16.5-20 ℃，多年平均降水量为1 367.5-1 932.9 mm^[27]。受季风气候影响，降雨主要集中在春夏两季，占全年降水量的70%以上。漓江流域植被资源丰富，类型多样且分布广泛。

1.2 数据来源

本研究应用到景观类型、地形、植被、土壤、气候、水文等多源遥感数据，数据概况如表 1所示。景观类型数据由李何英等^[28]提供，包含2000年、2005年、2010年、2015年以及2020年这5年的历史景观类型数据，数据精度为92%-97%；另外，还包含2030年的景观模拟数据，数据精度为83.57%。景观模拟数据考虑了两种情景：自然发展情景和生态保护优先情景。在自然发展情景下，假设漓江流域的景观格局不受新的政策因素影响，依据历史景观格局变化特征，按自然过程进行发展。而生态保护优先情景则以可持续发展为导向，严格控制建设用地的无序扩展，减缓林地和耕地面积的下降趋势。具体而言，生态保护优先情景要求林地、草地向建设用地的转化率减少约40%，耕地向建设用地的转化率减少约30%^[28]。

地形地貌[数字高程模型(DEM)、坡度(Slope)、坡向(Aspect)]数据采集自日本对地观测卫星(ALOS)。气候数据中的气温(Temperature)和降雨(Precipitation)数据来源于第5代ECMWF大气再分析全球气候数据(ERA5)。饱和水汽压亏缺(VPD)数据来源于美国加州大学计算的地表干旱气象数据集(GRIDMET)。蒸散发(ET)数据来源于中分辨率成像光谱仪(MODIS)卫星的MOD16A2产品。植被数据中叶面积指数(LAI)来源于MODIS卫星的MOD15A2H产品。总初级生产力(GPP)数据来源于MODIS卫星的MOD17A2H产品。植被覆盖度(FVC)基于Landsat5 TM和Landsat8 OLI/TIRS卫星影像光谱反射率数据，并采用像元二分法计算得到^[29]。日光诱导的叶绿素荧光(SIF)数据来源于全球太阳诱导叶绿素荧光产品(GOSIF)^[30]。土壤数据中的土壤类型(ST)和土壤有机质(SOM)指标分别来自中国科学院资源环境科学与数据中心和OpenLandMap土壤有机碳数据集。土壤湿度(SM)、土壤全氮(TN)和土壤全磷(TP)数据来自国家青藏高原科学数据中心。地表径流(Runoff)数据来源于ERA5。距离[到湖泊河流距离(D-river)和到人工表面距离(D-artificial)]是根据景观类型数据提取水域或建设用地(人工表面)景观斑块，并通过欧氏距离算法计算得出。旅游[旅游点密度(T-density)和宾馆酒店密度(H-density)]数据来源于去哪儿网(hwz ttps://hotel.qunar.com/)，并通过核密度计算得出。人口密度数据(P-density)来源于南安普顿大学发布的全球人口空间分布网格数据(WorldPop)。收入[农民人均纯收入(NetIncome)和人均财政收入(FiscalRevenue)]、农业发展[单位面积粮食产量(G-production)和单位面积肉类产量(M-production)]数据来自桂林经济社会统计年鉴(2000年、2010年和2020年)。自相关系数(LISA-I)通过对景观类型数据进行局部自相关计算获得。发电厂密度(PP-density)来源于世界资源研究所(WRI)的电厂分布数据，经核密度计算得到。石漠化敏感性(Karst)通过结合坡度、降雨、植被覆盖度和岩石性质等要素并参考张景华等^[31]、肖荣波等^[32]的研究成果计算得出。由于各指标的来源不同，空间分辨率存在较大差异，为了便于计算，统一重采样为1 000 m×1 000 m。

1.3 方法 1.3.1 景观生态安全评价模型

生态安全的状态，取决于生态系统所面临的外部干扰强度与其内部抵抗能力之间的平衡。通常情况下，外部干扰越强，生态系统自身的承受能力越弱，生态系统的安全性则越低^[33]。本研究选取景观干扰度指数(I_i)、景观脆弱度指数(V_i)、景观损失度指数(L_i)和景观生态安全程度(ES_j)来构建景观生态安全模型。

① I_i

I_i反映生态系统受人类活动因素干扰的程度，由景观破碎度指数(F_i)、景观分离度指数(S_i)以及景观优势度(D_i)三者加权求和得到，表达式为

$ I_i=w_1 F_i+w_2 S_i+w_3 D_i, $

(1)

$ F_i=\frac{N_i}{A_i}, $

(2)

$ S_i=\frac{1}{2} \sqrt{\frac{N_i}{A}} \frac{A}{A_i}, $

(3)

$ D_i=\frac{1}{4}\left(\frac{N_i}{N}+\frac{M_i}{M}\right)+\frac{1}{2}\left(\frac{A_i}{A}\right), $

(4)

式中，w₁、w₂和w₃分别为F_i、S_i以及D_i的权重，参照Liu等^[9]、杨伶等^[33]的研究成果，本研究将上述权重分别设置为0.5、0.3、0.2。N_i为i类景观的斑块个数，N为景观斑块总数；A_i为i类景观的面积，A为景观总面积；M_i为i类景观斑块出现的样方数，M为研究区总样方数。

② V_i

V_i用于衡量生态系统抵御外界干扰的能力，其值越大，表示生态系统的抗干扰能力越弱。通常，水域更易受到污染，而建设用地和林地的抗干扰能力相对较强。参考刘希朝等^[34]的研究，将水域赋值为5、耕地赋值为4、草地赋值为3、林地赋值为2、建设用地赋值为1，并进行归一化处理，归一化后的值定义为V_i。

③ L_i

L_i是对景观受到干扰后，生态系统受损失程度的描述，用I_i和V_i的乘积来表示:

$ L_i=I_i V_i。$

(5)

④ ES_j

ES_j是对区域景观生态安全状况的全面反映，可以通过公式(6)来表示:

$ E S_j=\sum\limits_{j=1}^n \frac{A_{j i}}{A_j}\left(1-10 L_i\right), $

(6)

式中，ES_j为第j个样区的景观生态安全，A_ji为第j个样区的第i类景观的面积，A_j为第j个样区的总面积。

由于景观在空间上的分布并不均匀，因此需要将研究区划分为多个小样区进行评价，这样才能更客观地反映不同区域的生态安全水平。根据漓江流域的范围和景观特征(平均斑块面积约0.28 km²，最大斑块面积约4 km²)，研究区被划分为1 588个2 km×2 km的样区。对每个样区进行生态安全评估后，使用普通克里金法进行空间插值，最终得到漓江流域的景观生态安全格局。本研究利用2000-2020年漓江流域的景观类型数据，采用景观指数法对2000年、2005年、2010年、2015年和2020年的漓江流域景观生态安全进行评价。以2015年的数据为基础，使用自然间断法将评价结果划分为5个等级^[35]：不安全(Ⅰ)、较不安全(Ⅱ)、中等(Ⅲ)、较安全(Ⅳ)和安全(Ⅴ)。

1.3.2 变化趋势

变化趋势可用Sen+Mann-Kendall方法进行分析。Sen+Mann-Kendall方法由Sen′s slope趋势分析和Mann-Kendal检验两部分组成。Sen′s slope趋势分析可表示为

$ \beta=\operatorname{median}\left(\frac{x_j-x_i}{j-i}\right), \forall j>i, $

(7)

式中，x_i和x_j为时间序列数据；β表示斜率，β>0时，表示时间序列数值呈上升的趋势，反之呈下降的趋势。Mann-Kendall检验是一种非参数统计方法，用于检测时间序列数据中的趋势，并且不需要数据遵循特定的分布，对异常值不敏感。通过结合Sen′s slope趋势分析和Mann-Kendall检验方法，可以对时间序列数据的趋势及其显著性进行综合评估。

1.3.3 RF模型

RF模型由Breiman^[36]于2001年提出，是由多棵决策树组成的集成模型，每棵决策树都有投票权，通过投票方式选择最优分类结果。使用RF的这种投票机制可以评估特征变量(因子)的重要性。该模型不仅具备处理高维数据、避免过拟合、训练效率高等优点，而且还能评估特征的重要性，并可分析不同特征之间的相互影响。即使在特征缺失较为严重的情况下，RF模型仍能保持较高准确度，因此被广泛应用于多个领域。

2 结果与分析 2.1 景观生态安全时空动态

漓江流域景观生态安全等级为Ⅰ等级的区域主要集中在流域的中西部地区，2000-2020年这些区域的面积呈现出先减少后增加的趋势(图 1)。景观生态安全等级为Ⅴ等级的区域主要位于流域的北部和东部，其面积在研究期间持续减少。从变化趋势来看，景观生态安全等级提升的区域主要位于漓江流域上游区域，包括青狮潭水库附近的九屋镇、潭下镇以及兴安县北部华江瑶族乡等区域(P < 0.05, 图 1)。而景观生态安全等级下降的区域分布较广，并以漓江下游的白沙镇、阳朔镇和高田镇等地最为显著(P < 0.05)。综合来看，2000-2010年漓江流域景观生态安全总体水平上呈现好转趋势；然而，自2010年起，景观生态安全总体水平开始下降，尤其是2015-2020年，Ⅰ和Ⅱ等级的区域面积迅速扩展，而Ⅳ等级的区域面积则大幅度缩减。

2.2 不同情景下的景观生态安全

基于2000-2020年的景观格局，对2030年漓江流域的景观生态安全进行情景预测，结果如图 2所示，各安全等级的面积分布详见表 2。通过对比图 1和图 2可以发现，在自然发展情景下，漓江流域的景观生态安全水平将进一步下降。尽管Ⅴ等级区域的分布面积与2020年相比变化不大，但Ⅰ等级区域的面积将进一步扩展，预测到2030年，Ⅰ等级的总面积将达到1 206.60 km²(表 2)，相较于2020年的1 009.17 km²，增加了197.43 km²。与之相比，在生态保护优先情景下，预测到2030年，Ⅰ等级区域的面积为1 074.76 km²，与2020年相比基本不变。

在预测的空间分布方面，生态保护优先情景与自然发展情景大体相似。高等级生态安全区域(Ⅳ和Ⅴ)主要集中在流域北部和东部，而低等级区域(Ⅰ和Ⅱ)主要分布在西部和南部。值得注意的是，在自然发展情景下，Ⅰ等级区域的扩展主要朝西和南方向进行，特别是桂林市辖区的秀峰往临桂方向区域以及阳朔县的白沙镇、高田镇，这些区域可能大面积扩展。然而，在生态保护优先情景下，这些区域可以实现平稳发展。与自然发展情景相比，生态保护优先情景将有效保障漓江流域的景观生态安全。

2.3 景观生态安全影响因素

从图 3可以看出，所选择的29个因子对漓江流域景观生态安全格局的解释力非常好，各年份的R²值均大于或等于0.97。在影响因素中，LISA-I是最重要的，对景观生态安全格局的影响远高于其他因子。除了LISA-I外，VPD、Precipitation、T-density、D-river、P-density、LAI、SIF以及DEM等因素也对漓江流域的景观生态安全格局产生较大影响。

漓江流域景观生态安全格局在不同时间段的影响因素存在差异。在2000年，除了LISA-I外，主要受DEM、T-density、D-river和PP-density等因素的影响(图 3)。到2010年，主要影响因素变为LISA-I、P-density、VPD、SM、Precipitation和SIF，地形、水文及旅游因素的相对重要性有所下降。到了2020年，VPD和LAI的影响力显著上升，分别跃居至第2和第4位。这些变化表明，漓江流域景观生态安全格局从2000年主要受地形、旅游、水文和农业生产等因素的影响，到2022年逐步转为受气候、人口密度和植被等因素的影响。

3 讨论

漓江流域是中国南方生态安全屏障的重要组成部分。本研究发现，2000-2010年，漓江流域的景观生态安全水平逐渐提高，但在2010年后呈现出下降趋势。李宁等^[37]、何毅等^[38]对过去20年漓江流域生态服务价值的研究也验证了这一趋势。2002年起，随着“退耕还林”政策的全面实施以及石漠化治理、天然林保护政策的逐步落实，漓江流域植被得到快速恢复，生态安全水平明显提升。然而，林明明等^[39]研究表明，近些年来桂林市生态环境治理和保护力度逐渐落后于生态环境压力的增长，这可能是研究后期漓江流域生态安全水平下降的原因。

在空间分布方面，本研究发现，景观生态安全水平较高的区域主要集中在漓江流域北部和东部，而景观生态安全水平较低的区域则主要分布在漓江流域西部和南部。景观生态安全的变化与景观类型的转移密不可分^[40]，不同景观类型下生态系统的组成和结构不同，从而使景观生态安全水平呈现不同的状态。漓江流域北部和东部大部分区域是山地，植被覆盖良好。其中，北部广西猫儿山国家级自然保护区和东部广西海洋山自治区级自然保护区长期以来得到较好的保护^[41]，成为漓江流域重要的水源涵养地，因此，北部和东部景观生态安全一直保持较高水平。而漓江中下游区域是整个流域人口最密集区域，同时该区域喀斯特地貌高度发育，农业生产和旅游开发活动频繁^[42]，导致景观生态安全等级相对较低。

对景观生态安全进行情景预测是区域生态环境管理与决策的重要依据^[43]。漓江流域在自然发展情景下，未来生态安全水平将会进一步恶化，且以桂林市城区以及流域中下游地区更为突出(图 2)。这主要是因为这些区域建设用地的迅速扩张大量侵占生态用地从而降低该区域的生态安全水平^[37]。相较之下，在生态保护优先情景中，通过严格限制建设用地的增长，并有效保护耕地和林地，可有效改善漓江流域景观生态安全水平。刘小琼等^[44]对长江经济带生态安全格局的多情景预测中也有类似的发现，生态保护优先情景下侵占的生态用地最少，生态安全水平优于其他情景。因此，实施生态优先政策是维护漓江流域生态安全格局稳定的可行策略。

景观生态安全的时空演变规律和驱动机制受多种因素的制约^[45]。本研究采用RF模型对漓江流域景观生态安全时空演变的影响因素进行剖析，发现人类活动(如旅游点密度、人口密度等)、水文、气候以及植被等因素对漓江流域景观生态安全的影响较为显著。这一结果与黄河流域主要受人均水资源量以及自然保护区面积占比等自然因素影响^[46]有所不同，而与三峡库区主要受人口密度、年平均气温和年降水量影响^[47]有相似之处。与黄河流域以及三峡库区不同，桂林市长期以来以旅游业为主导产业，旅游业的发展深刻影响漓江流域的各个方面，因此会对漓江流域景观生态安全产生较大影响。然而，随着桂林市产业结构的升级和转型^[48]，旅游开发对环境造成的压力逐渐减少，这一点可以从RF模型中旅游因素的相对重要性逐渐降低得到证实(图 3)。除旅游因素外，人口密度对漓江流域生态安全的影响也较为显著。近年来，桂林市社会经济快速发展，人类活动对区域景观的改造作用也不断加强^[49]，这导致人口密度因素对景观生态安全产生较大影响。此外，漓江流域还受到水文、气候和空间自相关等因素的影响，这与其特殊的生态环境密切相关。例如，水文因素影响漓江流域的旅游、养殖和种植业等多个产业的发展，气候因素则影响着漓江流域的植被生长状况，从而影响区域生态安全。在未来漓江流域生态安全的建设过程中，需继续关注旅游发展和其他人类活动对生态安全的影响，同时应高度重视区域气候变化和土壤状况的作用。

虽然本研究评估了漓江流域生态安全的时空变化，并综合分析了多种因素对漓江流域景观生态安全格局的影响，但是仍有诸多不足之处。如，在情景假设方面，只考虑了自然发展和生态保护优先两种情景，未来可进一步考虑耕地保护优先等情景。此外，虽然本研究综合了29种可能对景观生态安全格局造成影响的因素进行分析，但是未对这些因素之间可能的多重共线性进行评估，给研究的结果带来一定的不确定性。

4 结论

本研究基于景观指数模型对2000-2030年漓江流域景观生态安全格局进行分析，同时采用RF模型对影响漓江流域景观生态安全时空演变规律的因素进行分析，得出的主要结论如下：

① 2000-2020年，漓江流域景观生态安全的总体水平呈现出先上升后下降的趋势。漓江上游区域青狮潭水库附近的九屋镇、潭下镇以及兴安县北部华江瑶族乡等区域景观生态安全水平呈显著上升趋势(P < 0.05)，而漓江下游的白沙镇、阳朔镇及高田镇等地的景观生态安全则显著下降(P < 0.05)。

② 在自然发展情景下，未来漓江流域景观生态安全水平将会进一步下降，到2030年，景观生态安全处于Ⅰ等级的总面积将达到1 206.60 km²，相较于2020年增加197.43 km²；而在生态保护优先情景下，这一下降趋势将得到遏制，低等级区域(Ⅰ)无明显增加。

③ 漓江流域的景观生态安全格局受多种因素共同影响，其中人类活动(如旅游点密度、人口密度等)、水文、气候以及植被等因素对漓江流域景观生态安全的影响较为显著。此外，2000-2020年，影响漓江流域的景观生态安全格局的因素逐步由主要受地形、旅游、水文和农业生产等因素影响向主要受气候、人口密度和植被等因素影响转变。