董春媛, 黄海涛, 乔荣荣, 罗立辉, 常学礼
(1.鲁东大学资源与环境工程学院,山东 烟台 264025;2.山东省烟台第五中学,山东 烟台 264025;3.中国科学院西北生态环境资源研究院,甘肃 兰州 730000)
景观过程是当前全球变化研究中一个主要视角,在经典(狭义)景观生态学中介于区域和生态系统之间,是受人类活动影响和自然条件限制的重要表达尺度[1-2]。景观过程包括了斑块类型组成、结构在时空尺度上的多维度变化[3-4]。这种变化在时间上依赖研究过程中时间序列长短和分析节点疏密[5],在空间上取决于研究单元空间构型(组成要素大小、数量、形状和空间距离)和景观要素表达的基本单元大小[6-7]。在度量景观变化中有许多数量化指数可用来描述景观特征[8-10],其中基于景观组成要素面积计算的指数(多样性和分裂度)采用频率最高,因为该类型指标是基于景观组成要素面积信息熵和面积与斑块数量集成来实现的[11]。进一步看,空间异质性(多样性)尺度依赖特征(阈值)是景观生态学研究必须首先了解的问题[5,12-13],因为在任何地区进行景观变化研究中,恰当的时空分析尺度判断是保障研究结果可靠的基础。在已有以水域分布为特点的巢湖和南四湖研究中,李莹莹等[14]和梁佳欣等[15]分别采用3 km×3 km 和1 km×1 km 的分析网格获得了较为合理的计算结果。杨庚等[16]在山西平朔矿区景观研究中,采用0.5 km×0.5 km 分析网格获得较好的结果。毛忠超等[17]在对荒漠绿洲过渡带景观识别中采用1.5 m×1.5 km 网格获得了较好的结果。而杨馗等[5]在重庆江津研究中发现了2 个最佳分析尺度,分别为90 m 和3 km。由此可以看出,无论是在景观相似(不同湖区)或相异区(湖区、矿区和荒漠绿洲区),景观特征的空间尺度依赖规律都存在差异。同时,从景观特征时序变化趋势研究现状来看,目前对时间尺度上景观结构变化异质性考虑较少,多通过现有数据时间间隔分段分析比较,没有说明分段研究的合理性[12-13]。
此外,尺度效应还受到基础数据粒度(分辨率)和研究范围(研究区面积)的影响[9,18],由于篇幅限制这部分内容将在后续研究中完成。总的来看,在已有景观格局研究中,建立在景观多样性尺度依赖特征基础上的时空格局变化研究相对较弱[12-18]。所以,在分析尺度阈值分析明晰基础上,深入分析景观多样性格局时空变化对更好理解景观多样性变化机制是一个重要补充。
宁夏沿黄绿洲处于干旱、半干旱交错带,属于典型人工绿洲。近几十年高速的城镇化和农田扩张过程导致域内各种景观要素转换强烈,景观结构在绿洲规模扩展、水资源消耗和生物多样性保护等方面扮演重要角色[19-20]。为了深入解析宁夏沿黄绿洲区景观变化机制,本文拟在分析阈值确认的前提下采用斑块类型面积(CA)、景观多样性指数(LDI)、最大斑块指数(LPI)和相对分裂度指数(RSI)进行景观变化分析。拟解决的科学问题是:(1)宁夏沿黄绿洲景观多样性的尺度依赖特征;(2)在分析阈值约束下景观多样性格局时空变化与流转特征。
宁 夏 沿 黄 绿 洲(37°20′~39°20′N,105°00′~107°00′E)位于宁夏回族自治区北部,北起石嘴山、南止黄土高原、东到鄂尔多斯台地、西接贺兰山,面积约10831.3 km2(图1)。宁夏沿黄绿洲核心区是由黄河串联起来的卫宁灌域和银川灌域组成的人工绿洲区域,绿洲面积占总面积比例60%以上[21]。研究区内湖泊和湿地密布,人工灌渠纵横,主要有东干渠、西干渠、汉延渠、唐徕渠等。该区土壤以隐域性灌淤土和草甸土为主,天然植被以沿黄河分布的沙枣林和零散分布的灌丛湿地植被为主。主要植物物种有沙枣(Elaeagnus angustifolia)、枸杞(Lycium chinense)、柽柳(Tamarix chinensis)和芦苇(Phragmites australis)等[22]。该区是国家级沿黄经济区的核心地区,也是国家生态功能区划中的重点区域。
图1 研究区位置Fig.1 Location of the study area
本研究分析年份选择1975、1987、2000、2010 年和2020 年。土地利用类型数据来源分别是中国科学院资源环境科学数据中心基于Landsat MSS、TM/ETM 和Landsat 8卫星遥感数据解译的全国1980s以来的土地利用数据集(分辨率30 m×30 m),选取2000、2010 年和2020 年数据进行类型整合,分别为农田、乔木林地、灌木林地、草地、湿地、水体、人工地表和裸地8 类[23],总体解译精度(Kappa 系数)为0.772[24]。1975 年和1987 年由课题组下载的Landsat MSS 和TM 遥感数据(https://earthexplorer.usgs.gov/),在遵循上述产品分类原则并参考上述数据集中1980 年和1990 年的数据,在宁夏农林科学院专家野外指导和回顾访问调查下解译完成。
2.2.1 研究区边界确定分别以卫宁灌域最高引水口(沙坡头引水枢纽)和银川灌域最高引水口(青铜峡引水枢纽)的高程为限定条件,在ArcGIS 环境中用“con”命令对数字高程数据进行初级判定,这一过程生成的结果由多个图斑组成。其后以多年土地利用类型数据中的农田、湿地、水域等为辅助条件(缓冲源),在ArcGIS 环境中用“buffer”命令进行正负缓冲分析,最终形成宁夏绿洲区边界。
2.2.2 景观指数本文选取斑块类型面积(CA)、最大斑块指数(LPI)、相对分裂度指数(RSI)和景观多样性指数(LDI)4个指标完成相关分析。
LPIj为j类型中最大斑块占该类型的百分比(%),计算公式如下:
式中:Amax为研究区或分析单元中第j类土地利用类型(或类型分级)中最大斑块的面积(km2);CAj为研究区或分析单元中第j类土地利用类型(或类型分级)总面积(km2)
RSIj为研究区中第j类土地利用类型(或类型分级)的修改分裂度,RSIj的提出是为了强调原公式中存在的不同斑块(土地利用或其他分类)类型面积大小可比性未知而明确化的表述,变化在0~100%之间,RSIj越大分裂度越大,反之亦然。计算公式如下:
式中:SIj为McGarigal 等[25]原始定义;SIjmax表示CAj由相互独立分布的栅格组成,拥有j类型最大的分裂度;CAj为研究区或分析单元中第j类土地利用类型(或类型分级)面积(km2);BA 为输入数据的基础栅格面积,在本文研究中为9 km2(3000 m×3000 m);i为j类型中的斑块序列号。
借用植物多样性指数中Shannon-Weiner指数计算LDI,计算公式如下:
式中:Pj为j类土地利用类型(或类型分级)占分析单元总面积的比例(%);i为j类型中的斑块序列号。
2.2.3 LDI分级 本文采用基础土地利用数据分辨率为30 m,而且涉及到8类土地利用类型,故文中空间尺度确定从3×3个栅格(90 m×90 m)起始,确保最小分析尺度也有包含8种类型的保障。其后为减少计算量分别用10 倍栅格递增梯度进行设计。具体为: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个梯度来确定。分析过程在ArcMap 环境下采用邻域分析中block 和focal 命令分别对上述梯度对应单元中各土地利用类型的相对面积(Pi)进行计算,然后导出属性表用式(4)完成最终计算。景观多样性格局分析是建立在分级基础上,分级原则采用均值标准差法[26],将分级间距分为4级(表1)。景观多样性典型时段变化采用不同年份景观多样级别相减来进行,若前一个年份减去后一个年份LDI 增加(差值为正),则为好转区,反之(差值为负)为退化区,不变(差值为0)为稳定区。文中数据统计分析在SPSS中完成,显著性检查采用F检验,P<0.01 为较显著,P<0.001为显著。
表1 LDI分级Tab.1 Classification of LDI
从1975、1987、2000、2010 年和2020 年研究区全域LDI 变化来看(图2),2000 年是LDI 变化转折点,1975—2000 年LDI 呈下降趋势,由1.528 减少到1.142;2000—2020 年LDI 变化呈增加趋势,由1.142增加到1.274。总体来看,研究区LDI 变化呈倒抛物线下降趋势但不显著(P=0.130,R2=0.590),在2000年以后虽有增加但仍处在平均线之下。从LDI变化的尺度依赖特征来看,在分析的5 个年份中都存在显著一元二次方变化规律且都达到显著水平(P<0.001)。采用最大矩形(边长6000 m)捕捉到了LDI变化拐点,5 个年份重复检查可以判定边长3000 m矩形处前后是研究区LDI 变化转折点(图3)。当分析尺度小于阈值时,LDI增加非常陡峭,大于阈值时增加趋于平缓。所以,2000 年是研究区景观多样性变化转折年,3000 m×3000 m 的分析单元是研究区景观多样性格局变化分析的最佳尺度。
图2 1975—2020年LDI变化Fig.2 Change of LDI from 1975 to 2020
图3 1975—2020年LDI变化趋势与尺度依赖特征Fig.3 Change trends and scale dependent characteristics of LDI from 1975 to 2020
在确定最佳分析尺度基础上,对1975 年以来5个年份的宁夏沿黄绿洲区景观多样性进行分析。从图4 可见,1975 年和1987 年研究区LDI 高值区优势明显,2000 年和2010 年低值区优势明显,2020 年高低值优势趋向均衡,但与1975 年和1987 年相比高值区明显持续偏小,而低值区明显偏大。从变化的数量化特点来看(图5),LDI 高值区CA 在1975 年和1987 年明显处于高位,分别为5562 km2和4905 km2(图5a),而同期对应的RSI 处于低位,分别为0.18%和0.22%(图5c)。从2000 年开始变化特点明显相异,LDI 高值区CA 在2000 年和2010 年明显处于低位,分别为2403 km2和2565 km2(图5a),而同期对应的RSI 则相反,处于高位,分别为4.71%和4.24%(图5c)。到2020 年,LDI 高值区的CA 略有增加,为3843 km2,而RSI略有减小,为0.82%。
图4 1975—2020年LDI分级特征Fig.4 Classification characteristics of LDI from 1975 to 2020
从LDI 低值区的CA 变化来看与高值区明显相反(图5a),在1975 年和1987 年明显处于低位,分别为2452 km2和2857 km2(图5a),而同期对应的RSI处于高位,分别为11.70%和12.70%(图5c)。与高值区相反,从2000 年开始LDI 低值区CA 增加,在2000 年和2010 年分别达到6349 km2和6250 km2(图5a),而同期对应的RSI 落入低位,分别为0.52%和0.57%(图5c)。到2020 年,LDI 低值区CA 减少到4576 km2,而对应RSI则略又增加至1.08%。
与LDI 高、低值区格局指数(CA 和RSI)时序变化相比,较高值区和较低值区变化明显不同(图5b、d),其中较高值区CA 在1975 年和1987 年相对较高同为2421 km2(图5b),同期的RSI 分别为15.24%和17.32%。2000 年和2010 年较高值区CA 处在低谷,分别为1741 km2和1705 km2,到2020 年时较高值区CA 又增加到2128 km2(图5b);而同期对应的RSI 也呈相反规律,在2000 年和2010 年相对较高,分别为22.53%和25.10%,到2020年RSI略又减小至14.86%(图5d)。较低值区CA 在1975 年和1987 年也处在相对较高位,分别为2074 km2和2335 km2(图5b),同期的RSI 分别为21.12%和19.79%。2000—2020 年较低值区CA 呈平缓下降,分别为2025 km2、1989 km2和1962 km2;而同时期的RSI 变化相对复杂,分别为23.08%、20.25%和23.66%。
研究区LDI 时序变化分析表明,2000 年是景观多样性变化转折年(图2),为了减少赘述,以下仅选取1975—2000 年和2000—2020 年2 个时段进行分析。
结果表明,在1975—2000 年LDI 以退化为主;2000—2020年以稳定占优(图6)。从变化细节来看(表2),在1975—2000 年,退化区面积最大为6840 km2,占比54.3%;稳定区面积居中为4419 km2,占比35.1%;好转区最小为1332 km2,占总面积10.6%。2000—2020 年,稳定区面积最大为7848 km2,占比62.5%;好转区面积居中3915 km2,占比31.2%;退化区最小为792 km2,占总面积6.3%。从不同时期3种变化分区的流向特点来看(图7),1975—2000 年好转区有796 km2(60.5%)转换为后一阶段的稳定区;稳定区有3533 km2(80.1%)维持性质不变;退化区有3519 km2(51.5%)和3036 km2(44.4%)分别流向稳定区和好转区。其他类型间转化比例不高,最多不超过553 km2(稳定区流向好转区,占比12.5%)。
图6 不同转换期LDI分级变化特征Fig.6 Hierarchical change characteristics of LDI in different transformation periods
图7 不同转换期景观多样性分级转换流向Fig.7 Transformation direction of different landscape diversity levels in different transformation periods
表2 景观多样性分级格局变化Tab.2 Change of landscape diversity classification pattern
从分级后组成的景观特点来看(表2),在1975—2000 年退化区分布面积最大,其格局特点以斑块数量少(22 个)、RSI 小(0.18%)和LPI 高(46.5%)为标志。这一特点同样表现在2000—2020 年分布面积最大的稳定区,斑块数量仅为11 个、RSI 为0.12%、LPI 高达61.3%。而分别处在前后2 个时期分布面积最小的好转区和退化区,分别具有最小的LPI(2.2%和0.6%)和最大的RSI(8.6%和35.48%)。
基于宁夏沿黄绿洲土地利用数据(30 m 分辨率)的景观多样性时序变化分析表明,在1975—2020 年期间LDI 变化具有阶段性,在1975—2000 年呈下降趋势LDI减少了0.386,在2000—2020年呈增加趋势LDI 增加了0.132。这一过程主要受到了农田面积扩张、城镇化过程加剧以及生态环境建设(植树造林与子植被恢复治理)等因素的影响。这种现象在中国干旱半干旱地区人工绿洲研究中都可得到证实[13,19]。同时,从基础数据的关键指标农田、林地和人工地表面积变化也可得到辅证[27],在1975年到2020年农田由4597.6 km2增加到6342.5 km2,林地由335.2 km2增加到585.3 km2,人工地表由505.8 km2增加到1177.0 km2。需要指出的是研究区农田面积在2000 年达到最大为6700.7 km2,而同期LDI 处在研究期最低,说明绿洲农田扩展是导致景观多样性变化的最主要因素之一,也表明在区域景观多样性分析中时间尺度上的转折点判定对分时段研究具有重要参考价值。
基于1975—2020 年5 个年份尺度效应分析指出,在干旱半干旱地区发育于河流沿岸的人工绿洲,在景观类型分类保持在8类的条件下,当分析单元(矩形)为3000 m×3000 m 时LDI 变化趋于平缓,对图3 各年份LDI 尺度依赖特征平均水平进一步分析表明,在尺度分别增加1.5 倍(4500 m×4500 m)和2倍(6000 m×6000 m)时,LDI仅分别增加了0.069和0.092;在尺度减小0.5倍(1500 m×1500 m)时,LDI减少了0.123。因此可以推断,在干旱半干旱沿河绿洲区进行景观分析中最适合分析单元阈值在3000 m。此外,需要指出的是宁夏沿黄绿洲面积为10831.3 km2,考虑到尺度效应不仅受分析单元的影响,而且还受研究区基础数据分辨率和范围大小的影响[14-18]。所以,在其他干旱半干旱区人工绿洲进行景观规划设计和研究时可依据数据特点和绿洲规模大小,选择3000 m×3000 m 左右若干梯度做简约分析即可获得适宜分析尺度。
近10 a 来随着中高分辨率数据(2.5 m×2.5 m~30 m×30 m)土地利用与覆盖数据普遍使用,由此带来的景观多格局复杂性成为判断研究结果的主要干扰因素[28]。如何合理分级归类成为简化景观多样性图形表达的主要途径,也是等级理论认知的体现[29]。从采用均值标准差法分类结果的研究可以看出宁夏沿黄绿洲区LDI变化具有以下几个主要特点。首先,从不同时期各级别CA 和RSI 关系来看,在高值区两者之间关联系数为-0.927(P=0.023),在较高值区两者关联系数为-0.696(P=0.192),在较低值区两者关联系数为-0.886(P=0.045),在低值区关联系数为-0.931(P=0.021)。说明宁夏绿洲不同LDI级别CA 越大,其RSI 越小,这种现象在景观多样性分级格局变化分析中同样存在(表2)。其次,从LDI空间格局转换阶段性特点来看,1975—2000 年稳定区面积是2000—2020 年的56.3%,说明研究区在1975—2000 年景观多样性变化迅速且以退化为主,因为该周期LDI 好转区面积仅为退化区的19.3%。同时,2000—2020 年好转区面积由前一周期的1322 km2增加到3915 km2并不能说明该周期LDI 状况好于1975—2000 年,因为这一结果是基于2000年总体LDI 在研究时段最低而导致(图2)。上述分析表明,时序变化分析中在数据量不足以进行统计学趋势分析时,采用分级变化仍可实现LDI 可视化和数量化表达,这一特点对区域景观规划设计、生态系统健康评价和稳定性分析具有额外支撑作用。
最后需要补充说明的是,在采用文中所涉及的SI(Splitting index)指数进行景观表述时,其计算过程要慎重应用fragstats 软件,因为其结果出现异常概率时有发生与其定义相悖(1<SI≤景观或某斑块类型的斑块数量)[25],这种现象在最新研究中仍有出现[30]。
在区域景观变化研究中,确定分析指标在时空尺度上的变化转折点既是保障研究结果具有借鉴和共享的必要条件,也是区域景观多样性可视化表达与分析的基础。
(1)宁夏沿黄绿洲5 个年份重复分析表明,LDI在空间上具有明显尺度依赖特征,转折拐点在3000 m 左右,在其他相似地区应用时可针对绿洲规模大小在此尺度左右做简化分析进行判断。
(2)在时间尺度上,1975—2000 年LDI 呈降低趋势,其变化格局表现为退化区CA 最大、好转区最小;2000—2020 年LDI 呈增加趋势,表现为稳定区CA 最大、退化区最小。由于后期初始LDI(2000 年)在整体上处在最低,其好转程度没有达到前期(1975年)初始水平。
(3)景观多样性分级面积转换特征以前期好转区向后期稳定区流转(占好转区60.8%)和退化区分别向稳定区、好转区流转(分别占退化区51.5%、44.4%)为主要特征。不同景观多样性变化区格局以景观CA与RSI呈负关联为特征且具有普遍性。