基于改进遥感生态指数的宁夏沿黄平原区生态环境质量评价

董春媛,乔荣荣,杨智程,罗立辉,常学礼,*

1 鲁东大学资源与环境工程学院, 烟台 264025 2 南京大学生命科学学院, 南京 210008 3 中国科学院西北生态环境资源研究院, 兰州 730000

生态环境质量是指一定时间范围和区域内对影响其社会发展与人类活动的部分或全部生态因素进行定性分析与评判的区域性生态环境优劣程度[1]。区域生态环境质量评价是国民经济建设与可持续性发展规划的基础,为此国家环境保护部在2006年以行业标准的形式颁布了《生态环境状况评价技术规范》[2],其后在全国一些省、县行政单元进行了生态环境评价[3—6]。在应用中由于规范中的核心指标生态指数(Ecological Index,EI)不具备空间表达且对主观因素(权重)依赖较大,存在无法精准落实到应用规划图中的短板[7]。针对这一问题徐涵秋在2013年提出了一种基于遥感影像信息提取的遥感生态指数(Remote Sensing Based Ecological Index,RSEI),该指数可以很好的弥补EI指数不能进行空间表达的缺憾,并借助遥感技术数据获取容易、计算程序化等优点迅速在诸多地区生态环境质量变化与评价中得到应用[7—10]。RSEI是以遥感数据计算的绿度(用归一化植被指数(Normalized Difference Vegetation Index, NDVI)表示)、湿度(Wet Index,WI)、干度(用归一化裸土指数和建筑指数(Normalized Difference Building-Soil Index,NDBSI)平均值表示)和热度(用地表温度(Land Surface Temperature, LST)表示)四个指标为输入量,经过主成分分析(Principal Component Analysis,PCA)获得分量,然后用各因子在分量轴上载荷向量值判断其代表性和对轴的影响力[11]。与EI相比,RSEI具有客观、易获得和二维空间表达的优势,对后续区域发展规划和生态环境质量评价具有重要支撑意义[12—14]。但是,RSEI评价因子NDBSI与LST之间存在的高度正相关引起导入因子(图层)生态学意义重复带来的偏差没有得到重视;同时也没有考虑基于栅格计算的RSEI如何概括表达相邻栅格性质差异所带来的影响,如某一栅格只与另一类栅格相邻或与多种类型栅格相邻所产生的RSEI不同无法表达。那么如何修正上述RSEI计算中的漏洞?本研究拟采用景观多样性指数(Landscape Diversity Index,LDI)替换LST对普遍采用的RSEI指数进行改进,称其为改进遥感生态指数(Modified Remote Sensing Ecological Index,MRSEI),主要依据包括:(1)与其它三个指标相比LST反映的是短时间尺度地表温度状态,而LDI因子与其他因子属性相似都是相对较长时间尺度地表属性累积状态表达;(2)LDI计算考虑到了基本分析单元周边的组成类型数和其面积比例信息量。

在MRSEI应用中,若用LDI替代LST首先需要解决的问题是如何确定适合的LDI分析尺度,即:确定LDI尺度依赖特征。已有的景观尺度效应研究表明,在水域分布为特点的巢湖和南四湖分别采用3000 m×3000 m和1000 m×1000 m分析网格(基本分析单元)进行研究获得了较为合理的计算结果[15—16];而在山西平朔矿山复垦生态恢复区的研究则以500 m×500 m网格获得了较为可信的结果[17]。由此可以看出不论是在景观相似或相异区域其尺度依赖规律都存在差异[15—18]。因此识别研究区LDI尺度依赖特征不仅是MRSEI计算中首先要确定的问题,同时也是景观变化研究中恰当的分析尺度应用的基础。

宁夏沿黄平原处于干旱、半干旱交错带属于典型人工绿洲。近几十年高速的非均质化的城镇和农田扩张导致域内生态环境质量空间异质性发生了极大变化[19—20]。为了深入解析宁夏沿黄平原区生态环境质量空间分布特征,本文拟采用景观多样性指数LDI替代LST的RSEI计算方法对研究区生态环境质量空间分布格局进行分析。拟解决的科学问题是:(1)在干旱半干旱人工绿洲生态环境质量评价中MRSEI计算的合理性与生态学意义解释;(2)宁夏沿黄平原生态环境质量空间异质性与总体生态环境现状。

1 研究区概况

宁夏沿黄平原位于宁夏回族自治区北部,范围介于北纬 37°20′—39°20′,东经 105°0′—107°0′,面积约10831.3 km2(图1)。该区北起石嘴山,南止黄土高原,东到鄂尔多斯台地,西接贺兰山。宁夏沿黄平原核心区是由黄河串联起来的卫宁灌域和银川灌域组成的人工绿洲区域,绿洲所占比例60%以上[20]。研究区内湖泊和湿地密布,人工灌渠纵横,主要有东干渠、西干渠、汉延渠、唐徕渠等。该区土壤以隐域性灌淤土和草甸土为主,天然植被以沿黄河分布的沙枣林和零散分布的灌丛湿地植被为主。主要植物种有沙枣(ElaeagnusangustifoliaLinn.)、枸杞(LyciumchinenseMiller)、柽柳(TamarixchinensisLour.)和芦苇(Phragmitesaustralis(Cav.) Trin. ex Steud.)等为主。该区是国家级沿黄经济区的核心地区,也是国家生态功能区划中的重点区域。

图1 研究区位置Fig.1 Location of the study area 研究区涉及到各级行政边界来源于全国地理信息资源目录服务系统

2 数据与方法

2.1 数据

为了获得准确的景观多样性尺度依赖特征,分析过程采用研究区多年土地利用/覆盖数据平均值进行。其中2000年、2010年和2020年数据来源于中国向全世界免费分享的GLOBELAND30数据(http://www.globallandcover.com/)。1975年和1987年由课题组下载的Landsat MSS和TM遥感数据,在遵循上述产品分类原则并参考中国科学院资源环境科学数据中心上世纪80年代土地利用/覆盖矢量数据,通过野外和历史访问调查在宁夏林科院专家指导下解译完成。研究区土地利用/覆盖类型为农田、乔木林地、灌木林地、草地、湿地、水体、人工地表和裸地等8类。在Google Earth Engine(GEE)平台上对2020年7月19日和28日的Landsat OLI数据进行计算(网址:https://developers.google.com/earth-engine/datasets/catalog/Landsat_LC08_C01_T1_SR),研究区涉及到的轨道号分别为:129—33、129—34和130—34。

2.2 方法

2.2.1景观多样性指数与阈值

在已有的研究中RSEI在计算流程和评价因子选择上基本固定[21—22],在本文研究中为了避免存在评价因子生态学含义重复表达的问题(NDBSI和LST),采用分析确定后最佳分析尺度的LDI替换LST。本文的NDVI、WI和NDBSI计算与已有文献一致[7],LDI和RSEI的计算如下:

景观多样性指数(LDI):借用植物多样性中Shannon-Weiner index进行计算,计算过程在ArcGIS的邻域分析中完成。

(1)

式中,Pi为景观中第i类土地利用类型(或类型分级)占分析单元面积比例。

由于本文采用基础土地利用数据分辨率为30 m,而且涉及到8类土地利用类型,故文中不同尺度LDI计算从3×3个栅格(90 m×90 m)起始,确保最小分析尺度也有包含8种类型的能力。同时,为了减少计算量,分析单位按照栅格倍数逐级增大,从90 m×90 m起,其后依次为300 m×300 m、600 m×600 m、900 m×900 m、1200 m×1200 m、1500 m×1500 m、3000 m×3000 m、4500 m×4500 m、6000 m×6000 m 共9个梯度。最终LDI尺度依赖性确定采用1975年、1987年、2000年、2010年和2020年共5年的土地利用/覆盖数据计算结果平均值判定。

2.2.2评价因子归一化与MRSEI计算

针对上述4个输入因子在赋值单位和数值变化幅度上存在差异,需要进行标准化处理使指标值统一到0—1之间。经极差归一化处理后获得输入因子图层(图2)。

图2 PCA分析中的评价层Fig.2 Evaluation layer in PCAPCA: 主成分分析 Principal component analysis

(2)

式中, MRSEI为改进遥感生态指数,MRSEI值越大,生态环境越好;反之亦然。n表示主成分特征根累积达到90%以上的分量数,在MRSEI分析中因为输入因子为4,故1≤n≤4(i为整数)。EVi和ETi分别为某个评价因子的特征值和特征向量。其它因子定义同前。全域可视化表达是建立在采用归一化差异水体指数剔除研究区内水体信息后的基础上。此外,文中MRSEI空间异质性分析是建立在分级基础上,分级原则采用均值标准差法[21],将分级间距分为4级(表1中第一和第二列)。

表1 MRSEI分级

最后,整个研究区MRSEI水平(评价)可通过公式(3)计算结果与表1中分级阈值比较确定研究区整体生态环境质量级别。

(3)

式中,MRSEI为时间t的遥感生态指数,Mi为评价区域MRSEI第i类级别中值,PAi为评价区域MRSEI第i类级别相对面积。文中统计分析在SPSS和Excel中完成,回归方程显著性检查采用F检查(P<0.001为极显著),关联系数显著性检查采用R显著性检查。阈值分别为R4,0.05=0.811,R4,0.01=0.917和R4,0.001=0.991。需要补充说明的是本文中的MRSEI计算考虑到了PC轴数目累计贡献率要超过90%,很显然参与计算PC轴数大于等于1,与其它研究在方法上存在一些不同[7—10]。

3 结果

3.1 LDI的尺度依赖性

从研究采用的5个时段以土地利用/覆盖类型为基础计算的LDI平均值的尺度依赖特征来看,存在显著一元二次方程变化规律(P<0.001,R=0.988)。在分析矩形边长6000 m范围内可以捕捉到LDI变化拐点(在3000 m处,图3中的三角点)。当分析尺度小于拐点时,LDI增加非常陡峭,大于拐点时增加趋于平缓。因此文中LDI指数采用3000 m×3000 m的基本单元计算,计算完成后用ArcMap数据管理中的重采样功能,用近邻法重采样至30 m×30 m与其它三个因子图层在分辨率上保持一致。

3.2 MRSEI主成分分析特征

从PCA输出评价因子与MRSEI的关联分析结果来看(表2),在所有评价因子中NDVI与MRSEI关联程度最大为0.7878,但仍未通过最低统计学显著性检查(R4,0.05=0.811),说明采用任何一个因子都无法获得通过统计学显著性检查的MRSEI状况。同时,从各评价因子与MRSEI的关联关系来看特点非常明显,WI和NDBSI与MRSEI呈负相关,NDVI和LDI与MRSEI呈正相关。此外,PCA结果显示PC1轴占特征值信息的68.98%,PC2轴占28.76%,二轴合计承载了特征值97.74%。可以确定,采用此2轴依据公式(2)可获得所有参与评价因子绝大部分信息。从各评价因子对PC1轴信息集成的贡献大小(特征向量大小和正负关系)来看(表3),PC1主要贡献者是NDVI,特征向量为0.8901;其它3个评价因子特征向量为负值,变化在-0.4146—-0.0317之间。PC2主要贡献者是LDI,特征向量为0.9100,其次是NDVI特征向量为0.4056;WI和NDBSI特征向量为负且最大不超过-0.0200。

表2 评价因子与MRSEI 和RSEI的相关系数矩阵

表3 不同评价因子PCA分析特征向量与特征值

从LDI替换LST对PCA结果的影响来看(表2),在MRSEI中LDI与MRSEI正相关与RSEI中LST与RSEI关系相反,说明LDI引入避免了RSEI中LST与NDBSI对结果的影响高度一致重复表达的不足(二者关联系数分别为-0.9617和-0.9330)。同时,从统计学显著性检查结果来看,在RSEI中4个评价因子关联系数显著性都超过了0.05水平,说明在人工绿洲区采用任一因子单独评价区域生态环境质量也可获得具有统计学意义的结果,这明显偏离了PCA分析应用的主旨。而反观MRSEI与评价因子的关联分析结果有效的避免了这种现象发生(表2)。此外,从PCA向量特征来看RSEI分析中的PC1累计贡献率就在91%以上,从多维投影向量合成角度来看,RSEI分析采用的评价因子在结果贡献上具有相似作用(聚集于PC1);而在MRSEI中PC1贡献率不足70%,引入的LDI因子在PC2中占优势并使特征向量贡献率达到了28%以上(表3)。这一结果从细节上辅证了LDI替代LST在人工绿洲区遥感生态环境质量评价中的合理性。

从研究区总体环境质量评价角度来看,在MRSEI评价中由于PC1和PC2轴累积贡献率超过90%以上,故依据表3中第一、二列数据采用公式(2)在ArcMap栅格计算器中对PCA输出的第一、二图层分别计算获得了图4中的研究区2020年MRSEI分布现状图,其中白色区域是采用归一化差异水体指数剔除的研究区内的水体信息。其后通过对MRSEI分布现状图的属性表进行分析获得研究区MRSEI平均值和标准差并依据表1的原则进行分类获得宁夏沿黄平原区MRSEI生态环境质量分级图。最后逐级进行MRSEI中值计算(表1第4列)。

图4 宁夏沿黄平原区2020年MRSEI分布现状与分级Fig.4 Distribution status and classification of MRSEI of Ningxia Plain along the Yellow River in 2020

3.3 MRSEI空间异质性特征

从图4 的MRSEI分布现状来看,MRSEI低值区主要分布在研究区周边,即贺兰山东麓向宁夏沿黄平原过渡区和黄河东岸鄂尔多斯台地与黄土高原向黄河阶地过渡区。从图4 的MRSEI生态环境质量分级结果来看,“差”和“较差”级别在研究区周边呈连续分布格局,其中西侧MRSEI以“差”级别为主,东侧以“较差”级别为主。从不同MRSEI级别分布格局的数量化特征来看(表4),在分布面积方面,以“较差”级别面积最大为3528.7 km2,占研究区的31.9%;其次为“较好”级别,面积为3370.7 km2,占31.9%;其后依次为级别“好”和“差”,面积分别为1855.4 km2和1804.8 km2,所占比例分别为17.6%和17.1%。在斑块密度方面,从“差”到“好”级别梯度上,斑块密度为减少趋势,由“差”级别的8.3个/km2减少到5.9个/km2。说明随MRSEI增加其分布趋向集中分布,这一点从平均斑块面积变化特点中可以得到印证,在同一梯度方向平均斑块面积由0.120 km2增加到0.169 km2。在同类斑块空间距离方面,“差”和“好”级别斑块间相隔相对较远,分别为253.4 m和236.2 m,相互之间的连通性相对较低;“较差”和“较好”级别分别为210.3 m和206.4 m。相互之间连通性相对较高。

表4 不同MRSEI级别空间格局数量特征

最后,研究区总体MRSEI水平可依据表4中第一列数据与研究区总面积比和表1中分级计算获取的中值采用公式(3)完成计算。研究结果表明研究区整体MRSEI为0.0117,查表1中分级阈值可知,研究区生态环境质量刚好达到较好级别(研究区MRSEI值为0.0117大于较好级别最低阈限值0.0116)。

4 讨论

从2013年RSEI在中国首次提出并应用于福建省长汀县生态环境质量评价以来[7],RSEI的计算方法被许多研究者应用来分析区域环境质量[23—28]。这些研究都延续采用NDVI、WI、NDBSI和LST四个指标,这不仅忽视了遥感手段获取的LST是表征的短时间尺度地表属性量[29],而其它三个指标是相对较长时间地表属性累积的反映,也忽视了在许多区域类型中NDBSI和LST在生态学意义解释上高度相似的性质。因为大量遥感指数反演研究表明,地表温度与人工地表和自然裸地有极高的正关联[30—31]。此外,在RSEI分析中把研究基本单元(栅格)设定成相互独立状态,忽视了栅格与同类或异类或多种异类相邻在生态环境质量评价中的意义,因此用能够表达相邻栅格异质性的生态环境评价指标替代LST成为提高RSEI在生态环境评价应用中的可靠性必须面对的问题。

LDI是景观生态学研究中使用频率最高的一个指数,它是基本分析单元(或研究区)中组成类型多样性的综合表达,其大小与组成类型数和面积比例信息均匀程度相关。LDI在不同研究区都表现出随分析单元大小而变化的尺度依赖特征[32]。从宁夏沿黄平原LDI尺度依赖特征结果来看,当分析尺度(正方形边长)达到3000 m时LDI变化趋于平缓(图3),在尺度分别增加1.5倍(4500 m)和2倍(6000 m)时,LDI仅分别增了0.069和0.092;在尺度减小0.5倍(1500 m)时,LDI减少了0.123。因此可以推断,在干旱半干旱人工绿洲区进行景观多样性研究中最佳阈值在3000 m。此外,需要指出的是宁夏沿黄平原面积为10831.3 km2,考虑到尺度效应不仅受分析单元的影响,而且还受研究区范围大小的影响[33],所以其它类似地区应用时要考虑到研究区规模大小,可选择3000 m×3000 m左右若干梯度做简约分析获得适宜分析尺度。

从MRSEI评价过程中的PCA结果来看,用LDI替代LST计算获得各因子与MRSEI的关联系数都未通过统计学最低显著性检查(R4, 0.05=0.811),说明用LDI替代LST后保持了各评价因子对MRSEI的相互独立状态(表2),满足PCA分析n维向量投影条件。从PCA结果来看,PC1轴占特征值信息的68.98%,其中NDVI具有最大向量投射(0.8901)是决定因子,是影响PC1最大的因子,LDI在PC1轴上向量投射居第二(-0.4146)。其与NDVI比值为-0.4770,作用方向相反,是影响PC1轴性质的次要因子。因此,PC1轴在MRSEI中的生态学解释是NDVI和LDI是影响生态环境质量的主要特征分量,其值越大生态环境质量越好,这种作用在PC1特征值中居主导作用;相反LDI越大生态环境质量越差,但这种作用在PC1特征值中居次要作用。对现状的解释为在NDVI高值区,LDI往往处于低水平。其生态学意义为大面积高NDVI农田和乔木林地分布区对应地是LDI低值区,这样的区域在MRSEI计算中分值较高。PC2轴占特征值信息的28.76%,其中LDI具有最大向量投射(0.9100)是决定因子,LDI是影响PC2最大的因子,NDVI在PC2轴上向量投射居第二(0.4056)与LDI比值为0.4457,作用与LDI方向一致是影响PC2轴性质的次要因子。因此,PC2轴在MRSEI中的生态学解释就是NDVI和LDI是影响生态环境质量的次要特征分量,NDVI和LDI越大生态环境质量越好,反之亦然。对现状的解释为在LDI高值区,NDVI往往处于低水平,这样的区域在MRSEI计算中分值较低。其生态学意义为若任一分析单元相邻有多种土地利用/覆盖类型,其植被生产力对MRSEI贡献被弱化。结合研究区整体MRSEI水平为0.0117略超过“较好”级别低限(0.0116),说明在研究区中NDVI和LDI高值区重叠现象发生概率较低。从MRSEI空间异质性特点来看,一是以“差”和“较差”级别多分布在贺兰山东麓—宁夏平原过渡带、鄂尔多斯台地—宁夏平原过渡带和黄土高原—宁夏平原过渡带三个区域,呈环状围绕整个研究区。二是格局指数斑块密度和平均斑块面积与MRSEI变化具有较好的规律性,主要表现为在从“差”到“好”梯度上,斑块密度为减少趋势,而平均斑块面积呈增加趋势(表4)。总的来看,LDI替代LST进行的MRSEI分析首先保证了与原有三个因子(NDVI、WI和NDBSI)的低关联度,其中与NDVI最高仅为-0.3221(表2),没有发生对MRSEI造成重复贡献表达的结果。同时,基于MRSEI分级中值阈值评价宁夏沿黄平原区生态环境质量总体上处于较好水平。从具体应用角度来看,在ArcMap环境中对研究区内任一指定范围(或行政区)的总体生态环境质量水平评估可通过表2中值阈值和选定范围内各级别相对面积用公式(3)完成计算。

在MRSEI分析中研究区自然特征必须予以高度重视。在大多数情况下采用LDI替代LST在评价因子属性代表性和生态学解释上是可行的,因为其有效地避免了RSEI分析中评价因子NDBSI和LST之间存在的生态学意义重复表达和多因子向量投影中的高度聚集。但是在特殊地理单元如青藏高原等常年冻土分布区和不同沙丘(固定、半固定等)类型为主的地貌类型区域,LST作为极其重要的指标是否适用需要慎重考虑。但是LST所反映的短时间尺度属性必须通过恰当方法(多数据统计等)改进,增强其与其它评价因子在时间尺度代表性上相对匹配。

5 结论

在宁夏沿黄平原区基于土地利用/覆盖分类的景观多样性指数具有显著的尺度依赖特征(P<0.001),阈值出现在3000 m×3000 m。PCA解释了研究区生态环境质量主要受到NDVI和LDI影响,其中NDVI是PC1(特征值贡献率68.98%)的决定因子,特征向量为0.8901;LDI为次要决定因子,特征向量为-0.4146,该类型区在MRSEI计算中分值较高。LDI是PC2(特征值贡献率28.76%)的决定因子,特征向量为0.9100;NDVI为次要决定因子,特征向量为0.4056。研究区空间异质性主要以“差”和“较差”级别分布在不同土地利用/覆盖类型交错区且呈环绕研究区为主要特点。在MRSEI由“差”到“好”梯度上,数量特征主要表现为斑块密度为减少趋势由8.3 个/km2减少到5.9 个/km2,而平均斑块面积呈增加趋势由0.120 km2增加到0.169 km2。综合来看宁夏沿黄平原区生态环境质量总体MRSEI值为0.0117处在较好水平。