生态系统为人类提供自然资源和赖以生存和发展的生态环境, 是社会繁荣和经济发展的基础。近几十年来, 人类对全球生态环境的破坏日趋严重, 生态系统健康受到极大的威胁, 生物多样性衰减、土地沙漠化、盐碱化、石漠化、酸雨、温室效应、水体污染等生态环境恶化问题日益突出, 这些问题不仅威胁到人类自身的生存, 而且会危害到子孙后代的生存环境^[1]。生态健康研究是环境保护和可持续发展领域的热点问题^[2], 专家学者对不同类型生态系统^[3-10]及不同区域尺度^[11-14]开展了生态健康评价研究, 取得了一定的进展。生态健康评价方法主要有指示物种法和指标体系法, 由于区域或流域是多种生态系统的空间载体, 往往很难找到恰当的指示物种(种群)对其进行监测评价；指标体系法综合考虑生态系统的特征及其提供的服务功能, 通过构建特定的指标体系进行定量评价, 评价结果具有综合性, 因而被国内外广泛应用^[15-22]。生态健康评价的指标分类方法很多, 众多国内外学者从自身的研究角度出发, 构建了不同的评价指标体系, 其中以压力-状态-响应(Press-State-Response, PSR)及其扩展修正模型指标体系最为常用。在常用的生态指标体系中, 往往人为忽略或简化上层指标和下层指标之间的重叠和交叉联系, 呈现“树杈”形状；生态等级网络(Ecological Hierarchy Network, EHN)框架继承了“树杈状”生态指标体系的层次结构, 并通过补充上下层指标间的重叠和交叉联系, 可以更好地反映生态系统的完整性并合理的分配指标权重^[23-24], 被越来越广泛地应用。如何自适应筛选出能反映区域特性的关键指标当前生态健康评价面临的重要问题之一。

常用的指标筛选方法有灰色粗糙集法^[25]、层次分析法^[26]等, 或者通过组合多种方法进行指标筛选^{[16, 27-30]}。近年来, 出现了定量化的指标筛选模型, 如Lin等^{[24, 31]}和黄洁等^[32]提出一种使用EHN框架进行指标筛选, 将网络分析与量化指标选择过程相结合, 通过矩阵化分析和定量选取模型, 为指标的筛选工作提供了一种新的思路。但是这些研究对指标采用等权重分配, 且对Dale和Beyeler^[33]总结的指标选取需满足的8个标准的判断具有较强的主观性, 在一定程度上影响了结果的客观性。

新疆开都河—孔雀河流域(简称开孔河流域)地处我国西北部极端干旱区, 生态环境十分脆弱。加之近年来, 人类对自然资源不合理地开发利用, 使得开孔河流域面临着水资源的极度匮乏及严重浪费, 出现森林、草地退化, 湿地萎缩, 土地沙漠化、盐碱化等问题, 流域生态健康面临着巨大的压力^[34]。近年关于开孔河流域的生态问题, 主要从气候变化^[35]、人居环境^[36]、水资源优化配置^[37]和水资源承载力评价^[38-39]等方面开展研究, 针对开孔河流域生态健康状况开展全面评价的研究较少。

本文以开孔河流域为研究区, 广泛收集能反映流域生态健康状况的指标, 融合PSR框架和EHN框架, 构建网状生态健康评价指标体系及指标定量化筛选模型, 并对新疆开孔河流域进行生态健康综合评价, 为区域生态可持续发展策略的制定提供参考依据。

1 研究区与数据预处理 1.1 研究区概况

开孔河流域位于新疆中部, 介于83°50′—90°45′E, 40°19—43°18′N之间(图 1), 主要由开都河流域、博斯腾湖和孔雀河流域组成。流域内海拔高度在637—4812 m之间, 地势西北高、东南低, 山地、盆地、河谷、平原、荒漠和戈壁交错分布, 地形较为复杂。开孔河流域地处欧亚大陆腹地, 因远离海洋和高山阻隔, 具有明显的干旱大陆性气候特征。整个流域干燥少雨、蒸发强烈；多沙尘暴、浮尘天气, 气候条件极为恶劣, 生态环境较为脆弱。近年来, 由于人口的不断增长、经济的快速发展、水土资源过度和不合理的开发利用, 使得人类活动的干扰已经超过该地区自然资源承载力, 出现了高山草甸和山地草原退化、湿地萎缩、灌区土地次生盐渍化、土地沙漠化等严重的生态环境问题, 区域生态健康面临的威胁日益增大^[34]。

1.2 数据及预处理

为开展开孔河流域生态健康评价, 收集了遥感数据、气象数据和其他数据(表 1)。

遥感数据以MODIS产品为主, 在NASA官方网站(https://ladsweb.modaps.eosdis.nasa.gov/)下载。所有产品先利用MRT软件转换为UTM-WGS84投影；使用LDOPE工具提取质量数据, 利用TIMESAT软件中的S-G(Savitzky-Golay)滤波算法去除时序数据中的噪声；然后进行数据级转换, 获得真实的物理量；最后进行研究区域裁剪。MOD17A2H数据产品包括8天总初级生产量(Gross Primary Production, GPP)和净初级生产量(Net Primary Production, NPP)的信息, MOD17A3H产品提供年度NPP的信息, 来自给定年份的8天NPP产品(MOD17A2H)的总和。由于只获取到2001—2015年的年度NPP数据, 因此, 将2016和2017年的MOD17A2H产品中NPP数据求和获得年度NPP值。

气象数据来源于中国气象数据网(http://data.cma.cn/)。其中中国地面气候资料年值、月值数据集选取的气象要素包括年均气温、平均最高气温、年降水量、最大日降水量等, 基于样条插值理论的Anusplin气象数据曲面拟合专业软件进行气象数据插值；中国气象背景数据集空间分辨率为500 m, 包括≥0℃积温、≥10℃积温、湿润指数、干燥度等参数。

地形数据使用的是美国地质调查局网站(https://gdex.cr.usgs.gov/)提供的空间分辨率为30m的ASTER GDEM V2全球数字高程数据, 用于小流域提取和地形因子指标计算。土地利用类型数据、陆地生态系统类型空间分布数据、土壤质地数据及1km格网的人口密度数据均来源于中国科学院资源环境数据云平台(http://www.resdc.cn/)。土地利用现状遥感监测数据包括耕地、林地、草地、水域、居民地和未利用土地6个一级土地利用类型以及25个二级土地利用类型类型, 结合小流域边界用于计算景观生态学指数、土地垦殖率等指标；陆地生态系统类型空间分布数据主要用于后续的结果分析。土壤质地数据通过百分比来反映砂土、粉砂土和黏土三大类不同质地颗粒的含量, 借鉴土壤可侵蚀性因子的计算方法^[40], 作为土壤质地指数来表征研究区内土地质地的分布差异特征。

2 研究方法

本文以新疆开孔河流域为研究区, 开展流域生态健康评价研究。技术路线从层次上可分为三大部分：一是网状指标体系构建, 收集研究区数据并计算相关指标, 基于PSR框架和EHN框架, 构建生态健康评价网状候选指标体系；二是定量化指标筛选模型构建, 包括候选指标权重计算和基于最优化理论的筛选模型；三是研究区生态健康评价与分析。

2.1 网状指标体系构建

借鉴相关研究成果, 融合PSR框架和EHN框架, 构建目标层-准则层-要素层-指标层的4层指标体系, 其中网状结构仅体现于要素层和指标层之间(图 2)。以生态健康评价作为目标层, 压力、状态、响应为准则层。压力从土地和人口两个角度选取评价指标, 状态从活力、组织力、恢复力、地形地貌、水热状况5个角度选取指标, 响应从沙漠化敏感性、盐碱化敏感性、土壤侵蚀敏感性3个角度选取指标。候选指标选取兼顾科学性、数据易获取性和可量化性等原则, 分析要素层指标与指标层指标之间的重叠和交叉联系, 在网状指标体系中, 指标层指标与要素层指标的联系包括一对一和一对多的关系。本案例中指标层有31个候选指标。

由于评价指标的量纲及其所代表的物理意义存在较大的区别, 在进一步处理前, 必须对评价指标数据进行标准化处理, 所有指标值都归一化到[0, 1], 栅格大小统一为500m×500m。

2.2 候选指标权重的客观计算

将主成分变换(Principal Component Analysis, PCA)和熵权法相结合进行候选指标权重的客观计算。其中要素层指标通过PCA法确定权重, 指标层指标通过熵权法确定相对权重, 然后按照网状指标体系结构, 通过要素层对指标层指标权重进行约束, 获取初始权重, 最后通过权重二次分配的方式获取候选指标层指标的组合权重。

要素层指标均为抽象且无法直接量化的指标, 因此, 将每个要素层指标对应的多个指标层指标通过主成分变换方法集成, 获取要素层指标的具体数值, 然后将集成后所获取的要素层指标再次进行主成分分析, 根据评价指标对各主成分的贡献度来自动、客观地确定指标权重, 可以避免因人的主观判断所确定的指标权重对评价结果造成的偏差。

指标层的指标均有具体的含义和数值, 其相对权重采用熵权法^[16]进行计算, 基本思路是根据指标相对变化程度的大小来确定客观权重, 指标相对变化程度越大, 指标所分配到的权重也应越大。

2.3 评价指标自动筛选模型构建

根据流域生态健康评价网状指标体系的等级结构, 构建要素层和指标层的关系矩阵R, 其将用于约束定量化指标筛选模型。通过线性代数的矩阵来表示指标体系中要素层与指标层的相互关联性, 形成关系矩阵R。

式中, 行与要素层指标对应, K表示要素层指标的个数；列与指标层指标对应, J表示指标层指标个数；k为要素层的指标序号, j为指标层的指标序号, r_kj表示要素层指标k与指标层指标j是否有关联关系。如果二者没有关联, r_kj=0；如果二者有关联则, 则r_kj=1。

基于最优化理论, 根据要素层指标和指标层指标之间的关系矩阵的分析结果, 并在保证指标体系完备性的同时, 减少指标冗余, 需使最终入选指标数(目标函数Z)达到最小值。

指标筛选模型的目标函数Z设为：

同时, 需要满足如下3个约束条件：

式中, Z代表指标层(j)指标要选择的数目；x_4j代表备选评价指标d_j是否被选；X₄代表方程解集向量, W₄表示候选指标的权重向量, ′表示转置；目标方程表示指标筛选满足选取的指标层指标数量最少；约束条件(1)表示基于保障生态系统完整性的要求, 每个要素层指标(c_k)(k=1, 2, …K)所衍生的指标层指标(x_4j)(即d_j, j=1, 2, …J)中至少有一个对应被选取, 保证下层指标与上层指标具有连续性；约束条件(2)表示, 为了达到指标体系完备性的要求, 指标体系的权重之和要达到一定的精度要求, 本文参照主成分分析法的要求, 被选指标权重占原始指标权重的比例不低于85%；约束条件(3)表示备选指标x_4j(j=1, 2, …J)只能被选择(=1)或者被排除(=0)。

该最优模型可以利用遗传算法求解, 获得中选指标。由于去除了指标层中的部分指标, 中选指标权重和小于1, 因此在保持要素层指标权重不变条件下, 对中选指标的权重进行再分配, 使权重和为1。

2.4 生态健康评价指数计算

采用综合指标法, 即基于中选指标的标准化值及相应的重分配后的权重, 通过加权求和的方法计算整个研究区范围内的各评价单元的综合指标值构建生态健康综合指数(Ecological Health Comprehensive Index, EHCI)。

式中, EHCI为生态系统健康综合指数, 其反映的是流域生态健康程度, 取值范围为[0, 1], n为中选的评价指标个数, E_i为第i中选指标的标准化值, w_i为第i个中选评价指标的重分配权重, 且。

3 结果与讨论 3.1 定量化指标筛选结果与中选指标合理性分析

计算指标层候选指标的多年平均值, 通过主成分分析法和熵权法分别计算要素层和指标层指标权重, 然后获得候选指标的组合权重(表 2), 获得候选指标的权重向量W₄(即“组合权重”一列)。

将网状指标体系中候选指标权重向量W₄和关系矩阵R, 带入指标筛选模型中, 通过遗传算法求解模型。在本文的实验中, 交叉率为0.8, 变异率为0.05, 种群规模为200, 遗传的最大代数为1000。

利用遗传算法获取最优解的数目minZ=17, 具体的解集为X=[1 1 1 0 1 1 0 1 1 0 0 1 1 1 1 0 0 1 1 0 0 0 1 0 0 0 0 1 1 0 1], 其中1表示该指标中选。对筛选后的指标进行权重再分配, 使其和为1。中选指标及再分配后的权重如表 2最后一列。

从表中可知, 要素层指标所衍生的指标层指标中至少有一个对应选取, 保证了指标体系的完整性。在本实例中, 将生态健康评价指标从31个优化到17个, 用54.8%的指标表达了85.98%的信息, 能够较全面地反映生态健康评价的各个重要方面, 并且排除了人为主观影响。

通过文献调研, 共收集了30篇关于干旱区、半干旱区生态健康评价相关的文献, 进行评价指标使用频次统计。17个中选指标在30篇文献中使用的频次统计如图 3所示。可以发现所有的17个指标在30篇文献中的使用频次百分比都在20%及以上。归一化植被指数(Normalized Difference Vegetation Index, NDVI)、年降水量和植被覆盖度3个指标的使用频次百分比均超过了50%。中选指标在干旱区、半干旱区生态健康评价与分析的研究中使用频率较高, 说明本文的定量化指标筛选模型具有较高的合理性。

本文的定量化指标筛选方法是基于各指标多年平均值来进行指标筛选, 后续的生态健康评价可以直接利用中选指标及再分配后的权重进行计算。因此, EHCI计算只需要收集和计算中选指标不同年份相关数据, 不需要计算所有的候选指标。而常用的PCA降维方法, 变换后各主成分分量的物理意义不够明确；PCA各分量不同年份的指标权重值不同, 有可能在一定程度上影响结果的可比性；所有进行评价的年份都需要收集和计算所有的候选指标。因此, 与常用的PCA降维方法相比, 本文提出的定量化指标筛选方法具有更明确的物理意义和可推广性。

3.2 开孔河流域生态健康空间差异

利用中选指标及其再分配后的权重, 计算开孔河流域2001—2017年各年份的流域生态健康评价综合指数EHCI, 其中2001年、2010年、2017年和17年的平均生态健康指数如图 4所示。

开孔河流域生态健康17年来的空间格局变化不大, 整个流域生态健康空间异质性显著, 总体分布特征为西北高、东南低。生态健康状况差的区域主要分布在流域的东南部、开都河绿洲和孔雀河绿洲两大绿洲外围, 这些区域主要为荒漠、戈壁区, 植被极度稀少, 17年来生态健康差的区域变化不大；生态健康一般的区域主要分布在流域中部两大绿洲与荒漠、戈壁区过渡带以及西北部的高山区, 这些区域植被覆盖比较稀少；生态健康良好的区域主要分布在开孔河流域西北部的河谷地带、中部两大绿洲区, 这些区域有大面积的草本、沼泽、绿洲和农田, 整体上植被覆盖比较好。17年来, 生态健康良好的区域有所增加, 其中改善最明显的区域主要分布在孔雀河绿洲区。

3.3 开孔河流域生态健康的时间变化趋势

利用Sen+Mann-Kendall法^[41]对开孔河流域2001—2017年的EHCI指数进行趋势分析。将Sen趋势度β绝对值小于0.005的位置判定为无明显变化, 采用检验统计量Z来检验趋势的显著性, 本文取显著水平α=α=0.05, Z_1-α/2 =Z_0.975=Z_0.975=1.96。当β>0.005, |Z|>1.96时, 序列呈显著上升趋势；当β>0.005, |Z|≤1.96时, 序列呈上升但不显著趋势；当β < -0.005, |Z|>1.96时, 序列呈显著下降趋势；当β < -0.005, |Z|≤1.96时序列呈下降但不显著趋势。开孔河流域生态健康评价结果的趋势分析结果如图 5所示。

图 5清晰地表明了开孔河流域EHCI变化趋势的空间分布格局：生态环境明显改善的区域集中在两大绿洲区、和硕县西北部的平原和丘陵区以及和静县的部分湿地河谷区域, 并以孔雀河绿洲区生态环境改善最为明显。这一现象说明2000—2013年期间, 塔河流域管理局开展的国家级生态综合治理工程^[42]效果显著, 该工程的具体措施包括河道治理、灌区节水改造、生态输水、生态建设等^[43], 受生态治理工程的影响, 流域内植被改善趋势十分明显。由于区域发展农业的需要, 人们开垦新的农田使得平原绿洲灌区面积增大, 此外, 农田防护林网增密, 也在一定程度上使得绿洲灌区植被生长季NDVI持续增加。生态健康严重退化的区域面积较少, 主要集中在尉犁县靠近库尔勒市的绿洲和荒漠之间的一部分过渡区, 该区域生态环境脆弱；库尔勒市市区和零散分布在开都河绿洲区以及和静县高山区, 对库尔勒市区而言, 则与近些年该地区社会经济快速发展、建设用地扩展有密切联系^[44]。轻微退化主要集中在和静县高山区, 该区域退化的原因是气候变化、草原过度放牧和虫鼠灾害造成的高山草甸、山地草原退化^[34]。此外, 开都河绿洲和孔雀河绿洲的边缘区域也存在轻微退化。

阴俊齐等^[34]对开孔河流域生态环境环境10年(2000—2010年)变化评估的成果表明, 生态系统质量发生改善的区域集中分布于库鲁克塔格山, 孔雀河中上游的绿洲区、流域西北部的部分湿地区。陈曦等^[43]利用NDVI表征塔里木河流域2000—2013年的生态健康状况, 其研究成果表明, 开孔河流域的平原耕地增加速率在塔里木河流域“四源一干”区域中最大, 植被覆盖度呈现上升的趋势, 流域西北部山区的植被出现退化的趋势。两者的研究成果与本文的研究成果在空间分布上较为吻合, 在一定程度上验证了本文生态健康评价方法的合理性和评价结果的可靠性。

对EHCI趋势分析结果进行统计(表 3)。表 3表明开孔河流域生态健康整体上呈改善趋势, 但是部分地区也存在着严重退化现象。42.23%的地区生态环境保持基本不变；生态环境改善区域面积占41.56%, 其中明显改善的区域的比例达到了10.26%；生态退化区域占16.21%, 其中严重退化面积为1.61%。生态环境改善区域面积远大于退化面积, 约是2.5倍, 其中明显改善面积约是严重退化面积的6.4倍。说明受塔河综合治理工程的影响, 区域生态环境得到较大改善。综上所述, 开孔河流域生态环境变化趋势整体上是在稳定的基础上, 向着上升的趋势发展。

3.4 不同生态系统生态健康指数分析

研究区有7大生态系统类型, 包括农田、森林、草地、水体与湿地、荒漠、聚落和其他生态系统(裸土地、裸岩石砾地)。计算7大生态系统在2001—2017年EHCI的平均值和标准差(图 6), 其中平均值在一定程度上反映了生态系统总体的生态健康状况, 标准差可以反映各生态系统内部生态健康的差异。从图中生态系统均值箱线图可以发现, 农田生态系统的总体生态健康状况最好, 其次是水体与湿地, 最差的为荒漠生态系统。从箱体的长短来看, 农田和聚落的箱体较长, 说明这2个生态系统的EHCI均值在这段时间内的波动较大, 也表明这两个生态系统健康状况发生了较大的变化。各生态系统17年间的EHCI标准差的差异较大, 森林、草地、水体与湿地和聚落标准差值均大于0.10, 说明这4种生态系统内部健康状况差异较大, 其中森林生态系统内部健康状况差异最大, 其次是草地；荒漠和其他生态系统的标准差较小, 基本上小于0.06, 说明这两种生态系统内部生态健康差异很小；农田生态系统的标准差值居中, 说明其内部生态健康状况的差异大于荒漠和其他生态系统。从标准差箱体的长短来看, 农田的箱体较长, 说明农田生态系统的EHCI标准差在这段时间内的波动较大。

统计17年农田、森林、草地和水体与湿地生态系统的EHCI分布直方图(图 7)。可以发现, 农田生态系统的EHCI值主要分布在0.45—0.6之间, 水体与湿地生态系统主要集中在0.3—0.4的范围内；森林和草地生态系统的EHCI分布较为分散, 森林在0.2—0.3低值区和0.5—0.6高值区占比均较高, 草地在0.2—0.3和0.5—0.6范围内分别有一个波峰, 在0.3—0.4范围内有一个波谷。上述结果也验证了森林和草地的内部生态健康状况差异较大, 而农田和水体与湿地生态系统内部差异较小, 特别是农田生态系统, EHCI值集中分布在高值区, 因此该流域总体生态健康状况最好的为农田生态系统。

对不同生态系统的趋势分析结果按变化类型进行统计(图 8)。可以看出, 农田生态系统健康状况明显改善的比例最多, 超过了50%, 而呈现退化区域的比例仅占10.8%, 这与该区域大规模进行绿洲开发和保护, 与农田防护林网密度不断增加的客观情况相吻合。聚落生态系统健康状况明显退化的比例最多, 达到了8.98%, 表明该地区聚落生态系统退化最严重, 其次为森林和草地生态系统, 这主要是因为该区域尤其是库尔勒地区快速城市化发展进程和人口的增长使得区域生态环境遭到破坏, 同时由于区域经济和农业的发展使得部分森林和草地被转化为建设用地和耕地。荒漠生态系统生态健康状况基本稳定的比例最高, 达到了67.98%, 也是明显改善的比例最低的生态系统, 其他生态系统健康状况基本稳定的比例也接近50%, 这与二者受人类干扰程度较低, 生态系统稳定性较高的客观情况吻合。水体与湿地生态系统健康状况的改善比例超过六成, 明显退化的比例仅为1.32%, 这与该区域相对于其他生态系统, 水资源相对充足, 生态环境较好的客观情况吻合。

4 结论

在文献调研与统计分析获取的用于流域生态健康评价的常用指标和通用指标的基础上, 构建了基于PSR和EHN框架的流域生态健康评价的网状指标体系, 能更准确地表达复杂生态系统各因素的相互关联。提出了一种集成主成分分析法和熵权法的客观确权法, 并通过权重二次分配的方式获取了指标层指标权重的方法, 避免了人为主观影响, 指标权重分配合理, 能体现不同指标之间重要性的差异。在关系矩阵和指标权重的约束下, 基于最优化理论构建了入选指标数目最少的指标筛选模型, 定量化地实现了区域生态健康评价指标的筛选, 提高了指标筛选过程的科学性和客观性。后续研究可以进一步考虑简化指标体系结构, 如不设置要素层, 构建更为合理的约束条件等, 从而降低定量化指标筛选模型的应用难度, 提高模型可推广性。

在开孔河流域生态健康评价中, 将生态健康评价指标从31个优化到17个, 用54.8%的指标表达了85.98%的信息, 有效解决了传统指标设置主观、冗余的问题。开孔河流域生态健康总体分布特征为西北高、东南低, 生态健康状况差的区域主要分布在开孔河流域东南部的大范围的荒漠和戈壁区和中部绿洲区外围等自然生态环境恶劣区域, 生态健康良好的区域主要分布在流域西北部的河谷湿润地带和中部的两大绿洲区。生态健康状况变化趋势整体上是朝着改善的方向发展, 生态环境明显改善的比例达到了10.26%, 仅有1.61%的地区出现明显退化, 其中孔雀河绿洲区生态环境改善最为明显, 尉犁县靠近库尔勒市的绿洲和荒漠之间的一部分过渡区及库尔勒市市区生态健康退化最为明显。生态健康状况最好的是农田生态系统, 最差的是荒漠生态系统。17年来, 农田生态系统生态健康状况明显改善的比例最高, 聚落生态系统明显退化的比例最高, 达到了8.98%。开孔河流域生态健康的总体好转趋势主要原因是塔河综合治理项目实施以后, 塔里木河流域生态环境得到了有效改善。