生态系统健康是生态学研究领域的热点问题之一^{[1, 2]}，十八大以来，党和国家十分重视生态文明建设^[3]，而开展生态文明建设的基础就是生态系统的健康^[4]，因此，针对典型区域开展生态系统健康评价，及时发现社会发展中出现的生态环境问题，不仅对地区生态环境治理与恢复有重要意义，同时也对地区资源配置、土地利用、生态文明建设和社会经济可持续发展有着非常重要的参考和指导意义。

近年来国内外学者对生态系统健康评价主要采用指标体系法，即通过结合多学科多方面筛选出合适的指标，建立评价指标体系进行全面评价。Costanza^[5]早期提出的“活力-组织-弹性”(Vigor-Organization-Resilience，VOR)模型是该方法的原型，该模型从土地利用数据空间格局分布中选取相应指标构建评价指标体系，如Pan等^[6]、Chen等^[7]、张渊^[8]运用VOR模型分别对长江中游经济带、贵州省花江示范区和滇池流域等典型区域进行生态系统健康评价。VOR模型往往从土地利用数据空间格局分布中选取评价指标对生态系统健康进行评价，只关注土地利用方面的生态环境问题，评价结果与土地覆盖类型空间分布状态相关程度较高; 生态系统服务和权衡的综合评估(Integrated Valuation of Ecosystem Services and Trade-offs，InVEST)模型则是通过计算生境质量、碳储量等来评估生态系统服务价值，更关注于土地覆盖类型与自然环境之间的关系，可以较为准确地反映生态环境的实际承载力，评价结果更加侧重于生态服务价值方面，模型的准确性和有效性也已得到进一步证实。如李若玮等^[9]运用InVEST模型对青藏高原的碳储量进行估算，得到青藏高原东南部碳密度比西北部高的结论; 许宝荣等^[10]利用InVEST模型对兰州地区生境质量进行评价，得到兰州地区生境质量总体呈稳定状态，草地的生境质量变化较为频繁的结论; 尚俊等^[11]利用InVEST模型对鄱阳湖区1995-2015年生境质量进行时空变化特征分析，研究表明鄱阳湖区生境质量呈由高向低演化趋势。VOR模型和InVEST模型在生态健康评价方面都存在一定的片面性，因此，为了进一步克服它们在生态系统健康计算方面的局限性，通过合理的手段将这两种模型有效耦合，充分发挥它们的优势，将自然景观与生态系统服务价值相互统一、有机结合，可以更为全面、准确地反映所研究区域生态系统环境的实际状况。

本研究以广西北部湾经济区(以下简称“研究区”)2001-2020年的遥感数据、气象站数据为数据源，通过耦合VOR和InVEST模型，从自然环境角度和生态系统服务价值入手，选取生态系统健康评价指标构建“活力-组织力-恢复力-生态系统服务”多指标评价体系，利用最小二乘法优化组合赋权模型计算指标权重，对研究区整体和6个城市的生态系统健康状况进行连续长时间序列的定量评价，并提出促进生态系统健康发展的对策，以期为研究区可持续健康发展提供参考依据。

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

研究区于2006年3月成立，是我国面向东盟开放合作的重要区域和西南地区出海的主通道，以及对外开放的窗口。研究区范围包括南宁、北海、钦州、防城港4市所辖行政区和玉林、崇左2市物流区，逐渐形成“4+2”城市群，位于中国沿海西南端，地处北纬20°26′-24°02′，东经106°33′-110°53′，属亚热带季风气候区，气候温暖，阳光充足，多年年均气温21.6-23.0 ℃，雨水充沛，多年平均年降水量1 147.3-2 689.5 mm。研究区地理位置得天独厚，背靠大西南，面向东南亚，南部濒临北部湾，拥有广西全部的海岸线，是整个广西的海陆过渡带，是我国唯一临海、临边、临江的经济区，处在中国-东盟经济圈、华南经济圈、西南经济圈等多区域合作的交汇区域，是广西经济发展的龙头地区和对外开放的重点区域，对于国家实施区域发展总体战略和互利共赢的开放战略具有重要意义。

1.2 数据来源与处理

本研究数据来源主要由2001-2020年MODIS数据、Landsat数据和气象站点数据组成。MODIS数据主要包括MCD12Q1(土地利用数据)、MOD13A3[归一化植被指数(NDVI)数据]、MOD11A2[地表温度(LST)数据]、MOD16A3GF[潜在蒸散发(PET)数据]，来源于NASA(National Aeronautics and Space Administration)的数据共享网站(https://ladsweb.modaps.eosdis.nasa.gov/)。利用MODIS重投影工具对这几类数据进行镶嵌、投影和数据格式转换等处理，利用ArcGIS软件将这几类数据通过最大值合成法和重采样等操作，统一时间分辨率为1年，空间分辨率为500 m，其中MCD12Q1土地利用数据通过重分类方法分为7种土地利用类型^[12]，并利用Fragstats 4.2软件计算景观格局指数。Landsat数据通过谷歌遥感计算云平台(Google Earth Engine，GEE，https://code.earthengine.google.com)提取LANDSAT/LT05/C01/T1_SR、LANDSAT/LE07/C01/T1_SR和LANDSAT/LC08/C01/T1_SR 3个数据集，进行去云、裁剪、拼接等步骤，采用像元二分法对植被覆盖度(FVC)进行估算，最终重采样至分辨率为500 m，与MODIS数据相一致。气象站点数据来源于中国气象数据网(http://cdc.cma.gov.cn)，主要利用广西地区2001-2020年气象站降水数据，对缺失的数据采用线性插补方法进行插补并合成年数据，采用克里金插值方法对站点数据进行插值形成栅格数据，并利用研究区行政区划矢量数据进行裁剪，得到研究区2001-2020年降水栅格数据。

1.3 评价指标体系构建 1.3.1 评价指标选取

VOR模型从活力、组织力、恢复力3个方面对生态系统健康进行评价，活力是代表生态系统新陈代谢和生产力的重要指标，组织力是代表生态系统物种间构成复杂性的重要指标，而恢复力是代表生态系统抗干扰能力以及其稳定性的重要指标。InVEST模型主要有淡水、海洋、陆地3个板块的生态系统评估系统。本研究综合考虑研究区自然生态状况和生态服务价值两个方面，依据生态系统健康评价指标选取原则和相关研究^[13-16]，结合VOR模型中的活力、组织力、恢复力以及InVEST模型各个模块的特点对评价指标进行选取，从活力、组织力、恢复力和生态系统服务价值4个方面建立研究区生态系统健康评价体系。

评价体系的目标层为生态系统健康，准则层选择活力、组织力、恢复力和生态系统服务4个类别，目标层初步选取FVC、NDVI、LST、斑块密度(PD)、最大斑块占景观面积比例(LPI)、蔓延度指数(CONTAG)、香农均匀度指数(SHEI)、平均斑块面积(AREA_MN)、散布与并列指数(IJI)、生态弹性度、水源涵养、生境质量、碳循环等13个评价指标。

1.3.2 相关性分析和指标计算

在初步确定组成评价体系的评价指标后，为了满足评价指标选取原则，防止评价指标之间具有重复性，对选取的指标进行相关性分析，运用SPSS软件25.0版本中的相关性分析(Pearson相关系数)工具对目标层的各个指标进行相关性分析^[17]，结合研究区实际情况，去除LPI和CONTAG，最终优选11个评价指标，得到最终的生态系统健康评价体系并根据各自数据特点赋予生态系统健康意义，其中(+)代表指标具有正向驱动作用，(-)代表指标具有负向驱动作用。

各评价指标的含义和计算公式介绍如下。

$ \text { ① FVC }(+): \mathrm{FVC}=\frac{N-N_{\text {soil }}}{N_{\mathrm{veg}}-N_{\text {soil }}} \text {, } $

(1)

式(1)中，N为像元植被指数; N_soil为裸土或无植被覆盖像元NDVI值，N_veg为植被覆盖像元NDVI值。

$ \text { ② NDVI }(+): \mathrm{NDVI}=(\mathrm{NIR}-R) /(\mathrm{NIR}+R), $

(2)

式(2)中，NIR和R分别为近红外和红波段的反射率。

$ \text { ③ } \operatorname{LST}(-): T_{\mathrm{S}}=A_0+A_1 T_{31}-A_2 T_{32} \text {, } $

(3)

式(3)中，T_S为地表温度，T₃₁和T₃₂分别为MODIS第31和32波段的亮度温差，A₀、A₁和A₂是分裂窗算法的参数。

$ \text { ④ } \mathrm{PD}(+): \mathrm{PD}=\frac{M}{S} \text {, } $

(4)

式(4)中，M为景观中斑块数量，S为研究区总面积。

⑤ SHEI (+)：表示景观中优势斑块类型以及景观分布状况，其值在0到1之间。

$ {\rm{SHEI}} = \frac{{{\rm{SHDI}} \times {\rm{ }}TA}}{{CA}}, $

(5)

式(5)中，SHDI为香农多样性指数，TA为景观总面积，CA为某一斑块类型中所有斑块的面积之和。

$ \text {⑥ AREA_MN }(+): \text { AREA_MN }=\frac{S}{M} \text {, } $

(6)

式(6)中，S为研究区总面积，M为景观中斑块数量。

⑦ IJI (+)：用来描述不同类型景观斑块聚集程度，是描述景观空间格局的重要指标。

$ {\rm{IJI}} = \frac{{ - \sum {_{k = 1}^m\left[ {\frac{{{e_{ik}}}}{{\sum {_{k = 1}^m{e_{ik}}} }}\ln \frac{{{e_{ik}}}}{{\sum {_{k = 1}^m{e_{ik}}} }}} \right]} }}{{\ln (m - 1)}} \times 100, $

(7)

式(7)中，e_ik为i类型与k类型斑块相邻边长，m为景观斑块类型数目。

⑧ 生态弹性度(+)：用来表征生态系统受到干扰后维持或恢复其结构和功能的能力。

$ E=\sum\limits_{i=1}^n\left(A_i R_i\right){, } $

(8)

式(8)中，E为生态弹性度，n为景观数量，A_i为景观类型i所占面积比例，R_i为景观类型i的生态弹性度分值。

⑨ 水源涵养(+)：综合考虑年均降水量和实际蒸散发，结合气候、地形、根系深度和土壤深度等数据运用水量平衡法计算得到的指标。

$ Y_{i x}=\left(1-\frac{\mathrm{AET}_{i x}}{P_x}\right) \times P_x \text {, } $

(9)

式(9)中，Y_ix为第i种土地利用类型栅格x的年产水量，AET_ix为第i种土地利用类型栅格x的年潜在蒸散发; P_x为栅格x的年降水量。

⑩ 生境质量(+)：用于表征生境与土地利用类型的接近程度和土地利用强度。

$ H_{i x}=Q_i\left[1-\left(\frac{D_{i x}^2}{D_{i x}^2+K^2}\right)\right], $

(10)

式(10)中，H_ix为土地利用类型生境质量，Q_i为土地利用类型生境适应度，D_ix为土地利用类型生境胁迫水平，K为半饱和系数。

⑪ 碳循环(+)：用于表征总碳储量，它主要受各种土地利用植被碳密度值的影响。

$ C_{\text {all }}=\sum\limits_{i=1}^n C_i \times S_i, $

(11)

式(11)中，C_all表示总碳储量，S_i代表第i种土地利用类型的总面积，C_i表示第i种土地利用类型的碳密度。

1.4 最小二乘法优化组合赋权模型

通常在权重计算方面采用两种方法：层次分析法和熵权法。前者计算结果相对主观，后者计算结果过于客观。为了弥补单一赋权法的不足，达到主观与客观的统一，使计算的指标权重更加科学、合理、准确，因此建立最小二乘法优化组合赋权模型^[18]，以层次分析法和熵权法计算结果为基础，经过最小二乘法优化，增加权重计算结果的准确性。

假设评价指标综合权重为

$ R=\left(r_1, r_2, \cdots, r_n\right)^{\mathrm{T}} \text {, } $

(12)

建立最小二乘法优化组合赋权模型：

$ {\left\{ {\begin{array}{*{20}{c}} {\min T(R) = \sum {_{i = 1}^m\sum {_{j = 1}^n\left\{ {\left[ {\left( {{W_{j1}} - } \right.} \right.} \right.} } }\\ {\left. {{{\left. {\left. {{r_j}} \right){p_{ij}}} \right]}^2} + {{\left[ {\left( {{W_{j2}} - {r_j}} \right){p_{ij}}} \right]}^2}} \right\}}\\ {{\rm{ s}}{\rm{.t}}{\rm{. }}\sum {_{j = 1}^n} {r_j} = 1}\\ {{r_j}0,j = 1,2, \cdots ,n} \end{array}} \right.}\]{, } $

(13)

式(13)中，W_j1为层次分析法计算后的各指标权重，W_j2为熵权法计算后的各指标权重，r_j为各指标综合权重，p_ij为标准化后的评价矩阵。

1.5 多指标评价模型及等级划分

为避免指标数据来源和类型等不同造成各个指标间不能进行综合计算与相互比较的问题，本研究采用极差法^[19]对各个指标进行标准化处理，计算公式如下。

① 正向健康意义指标

$ Y=\left(X-X_{\min }\right) /\left(X_{\max }-X_{\min }\right), $

(14)

② 负向健康意义指标

$ Y=\left(X_{\max }-X\right) /\left(X_{\max }-X_{\min }\right), $

(15)

式(14)和式(15)中，Y为指标标准化值，X为指标实际值，X_max为指标观测数据中的最大值，X_min为指标观测数据中的最小值。

根据新构建的指标体系及权重，采用加权求和的方式构建生态系统健康多指标评价模型，计算公式如下：

$ \mathrm{EHI}=\sum\limits_{i=1}^n Y_i P_i, $

(16)

式(16)中，EHI为生态系统健康指数，Y_i为i项指标标准化后的值，P_i为i项指标权重，n为评价指标数量。

将研究区生态系统健康水平按照等间断法划分为5个级别^[20]，如表 1所示。

根据健康等级划分, 得到各等级面积占比和空间分布, 绘制研究区各等级面积占比图以及生态系统健康等级空间分布图, 对研究区整体生态系统健康状况从时间变化、空间分布以及空间变化角度进行分析。

2 结果与分析 2.1 评价指标相关性分析结果

评价指标相关性分析结果见表 2至表 4，表中数字1-13(不包括10)分别代表从目标层选取的FVC、NDVI、LST、PD、LPI、CONTAG、SHEI、AREA_MN、IJI、水源涵养、生境质量、碳循环等12个评价指标。其中1-3属于活力层内部指标，4-9属于组织力层内部指标，11-13属于生态系统服务层内部指标。10属于恢复力层指标，单指标无需计算相关性。

由表 2至表 4可以看出，在组织力层次LPI与CONTAG、CONTAG与SHEI、LPI与SHEI在0.01级别相关性分别达到1.000、0.994和0.996，均具有高度相关性; 而LPI与AREA_MN、CONTAG与AREA_MN、SHEI与AREA_MN在0.05级别相关性分别达到0.899、0.912、0.883，同样具有较高相关性。因此根据相关性分析结果，结合研究区实际情况，得到最终的生态系统健康评价体系并根据各自数据特点赋予生态系统健康意义。

2.2 最小二乘法计算指标权重结果

以层次分析法和熵权法计算结果为基础，经过最小二乘法优化计算出各指标的权重，结果如表 5所示。由表 5可知，在本研究所构建的生态健康评价体系中，组织力模块指标权重占比较大，为0.434 0;活力和恢复力模块指标权重占比大致持平，分别为0.150 8和0.155 8;通过耦合新增加的生态系统服务指标模块权重占比为0.259 4。

2.3 生态系统健康时间变化特征

图 1显示了研究区20年间各健康等级的面积占比。从图 1可以看出，研究区2001-2020年生态系统健康等级主要由二级、三级和四级构成，处于三级以上的区域面积占研究区总面积的85%左右，因此研究区生态系统健康整体处于一般健康水平以上。生态系统健康等级为二级和三级的区域面积占比变化起伏较大，而等级为四级的区域面积占比变化较为稳定。

研究区20年间生态系统健康等级为二级的区域(以下简称“二级健康区域”，其余等级区域的表述以此类推)面积占比变化整体呈现“M”型波动状态。研究区生态系统健康等级变化主要分为2001-2004年、2005-2008年和2009-2020年3个时间段，其中2001-2004年二级健康区域面积占比从68.51%上升到80.03%，年均上升2.88百分点; 但在2004-2005年，二级健康区域面积占比出现断崖式下跌，下降幅度为47.99百分点; 2005-2008年，二级健康区域面积占比相对上升，增长幅度为11.29百分点，说明该时间段研究区生态系统健康状况有向好趋势; 2008-2009年，二级健康区域面积占比又陡然升高，上升幅度为32.13百分点，虽然该时间段五级健康区域占比出现小幅度增长，但生态系统健康呈现由一般健康水平向亚健康水平转化的趋势; 2009-2020年，二级健康区域面积占比整体呈现阶梯式下降趋势，从2009年的75.45%下降到2020年的59.60%，年均下降约1.32百分点，三级和四级健康区域面积占比在这一时间段的变化起伏不定，均在2019年达到最大值，分别为25.33%和16.43%，因此2019年研究区生态系统健康处于最差状态。值得注意的是，2003年、2004年以及2009年的二级健康区域面积占比较高，这3年却均有一级健康区域，且面积占比分别为0.0003%、0.0003%以及0.0014%，虽然占比较小，但还是可以说明研究区生态系统健康在这3年处于较高水平。其中2009年生态系统健康处于一般健康水平以上的区域面积占比达到最大，为89.06%，该年为研究区近20年来生态系统健康最好状态年。

2.4 生态系统健康空间分布特征

图 2显示了研究区20年间生态系统各健康等级的空间分布。从图 2可以看出，研究区生态系统健康整体处于一般健康水平以上, 但2001-2020年生态系统健康各等级状态区域面积占比起伏变化均较为明显，研究区生态系统健康处于不稳定的状态，其按空间分布特征大体可以分为2001-2004年、2005-2008年、2009-2020年3个时间段。

2001-2004年研究区生态系统健康等级主要为二级，且二级健康区域面积占比呈现较为稳定的上升趋势，主要体现在崇左市和南宁市的三级健康区域向二级健康区域转变; 而这一时间段四级健康区域主要在玉林市的东南部，面积占比出现先下降后上升再下降的起伏变化趋势。2005-2009年研究区生态系统健康状况出现一定程度的好转，二级健康区域面积占比大幅度上升，变化最明显的区域位于南宁、崇左和钦州等市，主要是从三级健康区域向二级健康区域转变; 而健康等级为四级且面积占比变化最明显的区域主要位于玉林市东南部，从四级健康区域向三级健康区域转变。2009-2020年研究区二级健康区域面积占比呈现先阶梯式下降后上升的趋势，其中在2014-2015年和2018-2019年这2个时间段二级健康区域面积占比分别出现1个下降的“阶点”。2009-2014年研究区生态系统健康等级的变化主要发生在南宁市范围内，二级健康区域面积减少而三级增多; 在2014-2015年的“阶点”年份，南宁市西北部和东南部、玉林市西北部和西南部的二级健康区域面积占比下降，但四级健康区域面积占比上升; 从2015-2018年可以明显看到，南宁市北部、南宁市和钦州市的交界处三级健康区域面积占比增加，而四级健康区域面积占比减少; 在2018-2019年的“阶点”年份，南宁市北部、南宁市和钦州市的交界处以及玉林市中部地区，研究区生态系统健康水平出现小幅度下降，主要变化是二级健康区域面积占比下降和四级健康区域面积占比上升。

2.5 生态系统健康空间变化特征

为进一步研究2001-2020年间研究区生态系统健康的空间变化特征，对其生态系统健康等级空间变化情况较为突出的年份时间段进行分析，按照2001-2004年、2004-2005年、2005-2008年、2008-2009年、2009-2018年、2018-2019年和2019-2020年这7个时间段，将生态系统健康等级空间变化状态分为改善三级(五级→二级)、改善二级(五级→三级和四级→二级)、改善一级(五级→四级、四级→三级、三级→二级和二级→一级)、稳定、退化一级(一级→二级、二级→三级、三级→四级和四级→五级)、退化二级(二级→四级和三级→五级)和退化三级(二级→五级)7种类型，通过计算生成等级转移矩阵，并统计等级转移面积占比(表 6)，绘制生态系统健康等级变化空间分布图(图 3)。

由表 6和图 3可以看出，2001-2020年研究区生态系统健康等级转移面积除2004-2005年时间段外，其他时间段均呈现以稳定变化为主的状态。按改善型(改善一级、改善二级和改善三级)和退化型(退化一级、退化二级和退化三级)所占面积比例大小来分类，研究区生态系统健康等级转移可以分为改善型面积大于退化型面积类型(2001-2004年、2005-2008年、2008-2009年和2018-2019年)和退化型面积大于改善型面积类型(2004-2005年、2009-2018年、2019-2020年)。

在这20年期间，2004-2005年和2008-2009年研究区生态系统健康等级转移出现了较大的变动。其中2004-2005年研究区呈现退化型面积远远大于改善型面积和稳定型面积的空间变化特征。稳定型面积仅占研究区总面积的36.17%，退化型面积比改善型面积占比高47.93百分点，研究区生态系统健康状况剧烈恶化，生态系统健康等级退化的区域主要分布于崇左市、防城港市、钦州市的大部分地区和南宁市、北海市的部分地区。而2008-2009年研究区呈现改善型面积远远大于退化型面积的空间变化特征，前者比后者占比高34.03百分点，生态系统健康等级改善的区域主要分布于崇左、钦州和防城港等3市，研究区生态系统健康等级空间变化出现改善三级状态和退化三级状态，稳定型面积占研究区总面积的47.12%，但研究区整体生态系统健康状况呈现向好趋势。2001-2020年研究区总体呈现退化型面积略大于改善型面积的空间变化特征，而稳定型面积占研究区总面积的63.24%。生态系统健康等级改善的区域主要有南宁市边缘区域、北海市南部沿海区域、玉林市东南部地区、南宁市和钦州市的交界处。生态系统健康等级退化型区域面积占比为23.62%，分散分布于研究区各个市内。整体上，近20年来研究区生态系统健康等级转移主要以稳定型为主，生态系统健康呈现轻微恶化趋势。

2.6 研究区6市生态系统健康对比评价

为进一步探究研究区各市生态系统健康状况以及变化趋势，以研究区整体生态系统健康状况为基础，结合研究区行政区划矢量图，裁剪得到南宁市、崇左市、玉林市、防城港市、钦州市与北海市6个城市市区的生态系统健康等级空间分布图，将6市生态系统健康空间变化按2001年、2006年、2009年、2020年4个时间节点分为崇左市、防城港市、钦州市和南宁市、玉林市、北海市两组进行对比分析，如图 4、图 5所示。

2001-2020年南宁、玉林和北海3市与崇左、防城港、钦州3市的生态系统健康等级空间分布变化较为相同，生态系统健康水平均呈现先下降后上升再下降的趋势，但变化较为平缓。其中在2001-2006年，南宁市西北部地区三级健康区域面积增多，四级健康区域面积少量减少，南部地区三级健康区域面积明显增多，二级和四级健康区域面积明显减少，南宁市生态系统健康等级主要从二级和四级向三级转化; 玉林市中北部地区三级健康区域出现大范围扩张，二级和四级健康区域面积均出现一定程度的减少; 北海市的生态系统健康等级变化较为明显，北部、东部以及西部边缘地区三级健康区域面积大幅度增加，全市内四级健康区域处于较稳定状态，北海市的生态系统健康等级主要从二级向三级转化; 该时间段内3市的生态系统健康水平均出现一定程度的下降。在2006-2009年，南宁、玉林和北海3市三级健康区域面积均减少，二级健康区域面积增多，其中南宁市南部地区以及北海市东部地区四级健康区域面积出现小幅度增加，玉林市南部地区的生态系统健康等级从四级转化为三级，该时间段内3市生态系统健康水平均出现大幅度提高。在2009-2020年，南宁市和玉林市三级和四级健康区域面积均出现一定程度的增加，北海市中部地区四级健康区域呈破碎状扩张，在西部和南部地区出现小范围的五级健康区域，该时间段内3市生态系统健康状况均呈现恶化趋势。

综上所述，2001-2020年研究区6市生态系统健康等级空间分布呈现较为相同的变化趋势：在2001-2006年生态系统健康水平呈现下降趋势，在2006-2009年生态系统健康水平呈现改善趋势，而在2009-2020年生态系统健康水平又出现下降趋势; 崇左、防城港、钦州3市生态系统健康变化较为明显，而南宁、玉林、北海3市生态系统健康变化相对更平缓。

3 讨论

研究区的区位优势明显，区内具有丰富的生物资源和较高的植被覆盖率，开发潜力大，其生态环境受到学术界及政府的高度关注。本研究中，研究区近20年来生态系统健康整体处于一般健康水平以上，这与孙龙启^[21]利用“压力-状态-响应”(PSR)模型最终得到广西近海生态系统健康评价属于“良”的范围的结论基本一致; 研究区近20年来生态系统健康水平整体呈现“M”型波动变化，2005-2008年为二级健康区域面积占比史上最低时期，从2009年开始健康水平有较大程度的回升，这与研究区发展初期以牺牲环境为代价而追求经济高速增长的发展时间段相吻合。在2008年1月16日国家批准实施《广西北部湾经济区发展规划》(以下简称《规划》)，经济区被重新定位和贯彻可持续发展战略要求后，生态环境得到重视，总体生态系统健康得以恢复和进一步稳定。朱鹏飞^[22]研究分析得出广西沿海地区在2005-2010年建设用地和未利用地变化较为明显，总体生态环境状况呈下降趋势的结论，与本研究结果大致相符。刘少坤等^[23]测算出2008-2018年研究区在城市群一体化发展下，通过优化国土开发格局，各市的土地集约利用水平整体提升。相对应地，这期间生态系统健康变化较为平缓。然而不可否认的是，尽管有《规划》政策和具体的优化实施策略，但从长期来看，经济发展仍会对生态健康造成一些轻微损害，所以研究区在20年间生态系统健康呈现出退化型面积略大于改善型面积的时空变化特征，总体呈现轻微恶化趋势。本研究通过耦合VOR和InVEST模型，虽然增加了水源涵养、生境质量、碳循环等生态系统服务价值方面的评价，但是如上所述，经济发展对生态系统健康的影响不可忽略，因此，下一步计划扩展和丰富社会经济方面的相关指标，以期进一步完善评价指标体系，使评价结果更趋合理。

环境建设对研究区整体生态健康有着至关重要的影响，尤其是近年来，在“两山”理论和新发展理念指引下，生态环境保护与社会经济发展并重，因此相关部门应积极推动生态环境建设的有效落实，加大生态环境保护力度，合理调控和利用各类资源，以进一步改善研究区的生态环境质量。

4 结论

① 从时间变化来看，研究区近20年来生态系统健康整体处于一般健康水平以上，二级健康区域面积占比整体呈现“M”型波动变化，生态系统健康由亚健康水平向一般健康水平转变。

② 从空间分布来看，研究区生态系统健康总体上呈现北部、西部和南部地区优于中部和东部地区的空间分布差异。

③ 从空间变化来看，2001-2020年研究区生态系统健康表现出退化型面积略大于改善型面积的时空变化特征，总体呈现轻微恶化趋势。

④ 对研究区6市2001-2020年生态系统健康从时间尺度和空间尺度上对比分析发现，崇左、防城港、钦州3市生态系统健康状况变化较为明显，而南宁、玉林、北海3市相对较为平缓。