景观格局分析是景观格局与生态过程相互作用研究中的重要方法, 也是景观生态学研究的重要内容^[1]。其目的是将景观结构与其对生态过程的作用联系起来^[2]。景观指数是景观格局分析常用的工具, 因其简单、方便的优点而在景观生态学研究中得到广泛应用。尤其是在Fragstats软件出现后, 景观指数的应用更为广泛, 形成了“Fragstats指数家族”。但这些景观指数的生态相关性常常未经证实或者存在疑问, 许多是没有内在生态意义的数学构建, 其生态相关性更多时候是假设的而非真正建立的, 缺少经验性的证据^[2]。在土壤侵蚀研究中, 这类景观指数也在小区尺度^[3]到坡面^[4]、流域尺度^[5-6]得到大量应用。由于景观指数本身的性质、景观数据属性和土壤侵蚀过程自身的特点, 用景观指数描述土壤侵蚀变化存在固有的缺陷^[7]。认识到这些局限性, 在深入理解土壤侵蚀与景观格局的关系的基础上发展出多尺度土壤侵蚀评价指数^[8-9]、方向性渗透指数^[10]等耦合土壤侵蚀过程和景观格局的指数, 发展了基于土壤侵蚀过程和景观功能耦合景观格局的方法。尽管如此, “Fragstats指数家族”仍然是描述景观格局的简捷工具。对这些景观指数的有效性进行评估, 寻找正确应用的一些原则、使用前提和适用范围, 遴选出适用于指示景观格局变化与土壤侵蚀响应的景观指数, 有助于推进景观格局与土壤侵蚀相互关系研究。景观指数最根本的局限与联系空间格局和生态过程来建立景观指数的生态相关性有关^[2]。本文以碾庄沟流域为研究区, 通过景观指数与用侵蚀模型模拟的土壤侵蚀量和输沙量的关系分析, 阐释景观指数对土壤侵蚀过程的表征能力和适用性, 探讨有效景观指数的筛选原则。

1 研究区概况

碾庄沟位于黄土高原北部(36°37′－36°45′N, 109°26′－109°37′E), 为延河中游一级支沟, 在延安市宝塔区李渠镇汇入延河。碾庄沟流域总面积54.2 km²。多年平均降水量527 mm, 地处温带半湿润半干旱地带, 大陆性季风气候。海拔926－1278m, 地形破碎, 属典型的黄土丘陵沟壑区, 梁峁、沟谷占总面积的90%以上。沟壑密度达到2.74 km/km^{2 [11]}。主要为可侵蚀性高的黄绵土覆盖。流域土壤平均颗粒组成为砂粒20%, 粉粒55%, 粘粒25%^[12]。碾庄沟是小流域治理试点流域, 淤地坝建设和植被恢复并重^[13-14]。近年来, 随着经济的发展, 居民点、工矿交通等建设用地的扩张成为流域内土壤侵蚀的重要影响因素。

2 研究方法 2.1 景观格局指数选择

对土壤侵蚀而言, 斑块水平景观指数指示性不强。因此, 本文主要从斑块类型水平和景观水平两个水平遴选景观指数。表征景观格局同一特征属性的不同景观指数相互冗余^[15], 部分指数即可代表全部指数。与土壤侵蚀相关的景观特征主要是景观破碎度, 因此, 选取的景观指数包括4个类型：(1) 斑块、边界密度指数；(2) 形状类指数；(3) 描述聚集与分散的指数；(4) 斑块类型多样性指数。选取的各类景观指数及其含义和分析尺度如表 1所示。

2.2 流域侵蚀模拟

流域土壤侵蚀模拟采用WATEM/SEDEM模型 。WATEM/SEDEM模型是基于栅格的土壤侵蚀模拟模型, 由侵蚀量模拟、泥沙传输能力估算和泥沙在流域内传输模拟三部分组成, 在年尺度上模拟产沙和泥沙传输。土壤总侵蚀量为每个栅格年产生的泥沙总量, 采用修正通用土壤流失方程(RUSLE)^[19]计算。

(1)

式中, E是栅格年土壤侵蚀强度(kg m^-2 a^-1)；R是降雨侵蚀力(MJ mm m^-2 h^-1 a^-1), 为排除降雨变化对模拟结果的影响以更好地反映土地利用格局改变的侵蚀效应, 在模拟中都取1992年的降雨侵蚀力; K是土壤抗蚀性因子(kg h MJ^-1 mm^-1)；L为坡长因子, S为坡度因子, 二者合称地形因子；C是植被覆盖管理因子；P是土壤保持措施。C、P因子参考^[20]和刘宝元等^[21]在黄土高原地区的研究结果。

WATEM/SEDEM模型中年泥沙传输能力(T_C)(kg m^-1 a^-1)指地表径流向下坡搬运泥沙的最大潜能, 定义为单位长度等高线上能流出的最大泥沙量^[17]。对每个栅格来说, T_C就是该栅格处水流能搬运的最大泥沙量与栅格下坡方向边长的比值, 单位为kg m^-1。T_C计算公式为：

(2)

式中, K_tc为传输能力系数(m), 理论上为均一坡度和径流条件下产生足够的泥沙以达到栅格传输能力所需的上坡距离^[17]；R、K、LS与RUSLE中相同, 分别为降雨侵蚀力因子、土壤可蚀性因子和地形因子。

在确定年侵蚀量后和泥沙传输能力后, 基于侵蚀量和传输能力确定泥沙在传输路径上的再分配。泥沙传输的路径是从泥沙源区沿地形梯度向河道的水流路径。WATEM/SEDEM模型采用八邻单一流向算法获取路径。当泥沙供给(上游输入的泥沙与当地产生的泥沙之和)大于泥沙传输能力时(传输制约), 泥沙向下游传输的量达到最大(等于T_C), 同时沉积发生, 沉积量为输入量与传输量之差；当获得的泥沙供给量小于泥沙传输能力时(侵蚀制约), 没有沉积发生, 向下游的传输量为获得的全部泥沙供给量。模型最终输出结果包括每个栅格的上坡泥沙输入量、输出量、侵蚀量、沉积量和到达河道或流域出口的泥沙总量。

通常利用流域泥沙观测数据进行模型校验来寻求最佳的参数组合来获得K_tc值^[22-25]。将土地利用/覆盖划分为具有K_tc高值的低植被覆盖类型(耕地、裸地、工矿、交通建设用地)和具有K_tc低值的高植被覆盖区(林地、草地、灌丛)来获取适合区域的K_tc值是常用的方法, 对应的K_tc分别为K_{tc_high}和K_{tc_low}.根据Liu和Fu^[11]的研究, 碾庄沟流域的 K_{tc_high}和K_{tc_low}分别为20和15。

2.3 景观格局数据获取和景观指数计算

土地利用/覆盖通过Landsat卫星TM影像解译获取, 类型包括坝地、耕地(包括沟坡、峁坡的梯田和坡耕地)、林地、灌丛、水域、居民地(包括工矿等建设用地)、草地、裸土。共提取1992、1995、1999、2002、2006、2009年六期土地利用/覆盖数据。DEM为1 : 10000地形图数字化生成。在流域尺度上, 淤地坝对地形的改变最为显著。利用遥感解译的淤地坝图层, 对地形进行修正以获取当年地形, 每座淤地坝库内的高程取坝库内最高点的高程值。降雨采用邻近沟口的延安气象站观测数据。景观类型在土壤侵蚀过程中发挥着抑制土壤侵蚀发生发展或促进侵蚀过程发展的作用, 前者为“汇”类型, 后者为“源”类型。因此, 土地利用/覆被类型可划分为“源”“汇”两大类型。属于“源”类型的土地利用/覆盖类型包括林地、灌丛、草地、坝地和水域。“汇”类型则包括裸土、居民地、耕地(不包括坝地中的耕地)。景观指数利用Fragstat 4.2计算, 分别计算按“源”、“汇”类型归并前后景观尺度和斑块类型尺度的景观指数。

3 结果分析 3.1 流域景观格局变化

从1992年到2009年, 流域景观类型结构的发生了较大的改变。主要表现为：林地、草地、灌丛等植被覆盖区的面积比例大幅度增加和耕地、裸土持续减少, 主要转变为林地、草地；沟谷内居民地和坝地逐渐扩张。2002年以前, 草地的扩张幅度最大, 之后林地增幅最大, 部分草地也转变为林地或灌丛。1992年林地、草地、耕地和裸土占碾庄沟流域面积的绝大部分。到2009年, 林地成为流域内比例最大的景观类型, 占流域总面积的50%以上(图 1)。

这一时期, 流域景观格局发生了显著的变化。在景观水平上, 斑块密度(PD)、边界密度(ED)和景观形状指数(LSI)持续降低, 斑块结合度(COHESION)、相似邻近百分比(PLADJ)、分割指数(SPLIT)以及斑块类型多样性指数指数(SHDI)都呈现下降趋势(图 2)。这表明流域景观破碎度降低, 聚合度在增强。

在类型水平上, 除居民地外, 土壤侵蚀的各类“源”景观类型破碎度在增加, 聚合度在降低(图 3), 而各类“汇”景观类型破碎度都呈下降趋势, 聚合度持续上升(图 4)。对所有的景观类型, 最大斑块指数(LPI)是一个极不稳定的指数, 没有稳定的变化趋势。将“源”、“汇”类型分别合并后计算得到的景观指数变化趋势表明：“汇”类型PLAND上升, PD、ED、LSI、SPLIT呈下降趋势；“源”类型PLAND下降, PD、ED、LSI、SPLIT呈上升趋势(图 5)。“源”、“汇”两个类型的COHESION值都处于较高的水平, 各年全都高于90%, 表明它们都具有较高的聚合度。对于“源”、“汇”两种类型, 最大斑块指数(LPI)也是一个极不稳定、没有稳定的变化趋势的景观指数。

3.2 流域土壤侵蚀变化

WATEM/SEDEM模型模拟的流域总侵蚀量为每个栅格产生的泥沙量之和。整个流域总侵蚀量逐渐下降, 流域平均侵蚀模数从1992年的10.05 t hm^-2 a^-1下降到2009年的7.41 t hm^-2 a^-1(图 6)。WATEM/SEDEM模拟的流域泥沙输出量等于到达与沟口连通的河道的泥沙总量。受分布于沟谷中的淤地坝控制, 流域内大部分区域并不是流域出口泥沙的来源。直接与沟道连接的泥沙产生区(有效集水区)对流域泥沙输出的贡献量是最主要的。由于流域内产生的泥沙大量被淤地坝拦截, 加上坡面植被对土壤侵蚀的抑制作用, 流域内与出口畅通连接的有效泥沙源区面积远小于流域总的泥沙源面积。因此, 流域泥沙输出到河道的总量远小于侵蚀总量。由于近河道区域是对侵蚀贡献较大的居民地和沟坡裸土的主要分布区, 且居民地持续扩张, 尽管流域侵蚀总量大幅下降, 但流域泥沙输出总量却呈现波动变化趋势, 2009年泥沙输出最高, 达到3906.25 t/a(图 6)。

4 讨论 4.1 景观指数指示流域土壤侵蚀变化的有效性 4.1.1 类型水平景观指数

景观指数与土壤侵蚀具有显著的相关性是其适用于表征景观格局与土壤侵蚀关系、指示景观格局变化的土壤侵蚀效应的前提。根据景观类型在土壤侵蚀过程中的功能, 可将其分为“源”、“汇”两种类型^[26]。一般而言, 土壤侵蚀中起“源”作用的景观类型比例的减少和破碎化有助于降低侵蚀, 反之则增加侵蚀；而对于“汇”功能类型, 其比例增加和破碎化会增强侵蚀, 反之则减弱侵蚀^[27]。林地、灌丛、草地、坝地和水域在流域土壤侵蚀过程中主要发挥“汇”功能, 而耕地、居民地和其他裸露地表都是土壤侵蚀过程中主要的径流和泥沙源区^[26], 是土壤侵蚀过程中的“源”景观类型。对“源”景观类型：PLAND、LPI、PLADJ、COHESION指数值的增加意味着其面积比例增加, 破碎化降低, 聚合度增加, 对应着土壤侵蚀的增强；PD、ED、LSI、SPLIT指数值增大则意味着破碎化增大, 对应着侵蚀的减弱。对“汇”景观类型, PLAND、LPI、PLADJ、COHESION指数值的增加意味着面积比例增加, 破碎化降低, 聚合度增加, 对应着土壤侵蚀的减弱；PD、ED、LSI、SPLIT指数值增大则意味着破碎化增大, 对应着侵蚀的增强。理论上, 对“源”、“汇”两类景观类型, 有效表征景观格局变化的土壤侵蚀效应的景观指数与土壤侵蚀表征变量间的相关关系应具有一致性, 即对同为“汇”或“源”景观类型的多个景观类型, 景观指数与土壤侵蚀变量之间都为正相关或负相关, 相关系数具有相同的符号。

如图 7所示, 类型水平景观指数与流域土壤侵蚀模数相关系数的符号一致性存在差异。对多个“源”景观类型的, 除了PD、LSI外, 同一景观指数与侵蚀模数的相关系数符号都不一致(图 7)。对多个“汇”景观类型, 各景观指数与流域平均侵蚀模数相关系数的符号都不一致。即使是对林地、灌丛、草地这三类起土壤侵蚀过程中的“汇”景观类型, 除LSI与土壤侵蚀模数都呈正相关外, 其他指数都与侵蚀模数相关系数符号都不一致(图 7)。在多个“源”景观类型格局与流域输沙量的相关关系表达上, 除SPLIT指数外, 都表现为一致的正相关或负相关, 且正确反映了“源”景观类型格局与输沙量的关系(图 8)。在表达多个“汇”景观类型格局与流域输沙量的相关关系时, PD、ED和LSI都具有符号一致性, 都为正相关或负相关, 而其余指数则不具有这种一致性(图 8)。

正确响应对土壤侵蚀起正、负效应的景观类型的消长是景观指数能耦合景观格局变化与土壤侵蚀效应的前提。土壤侵蚀中起“源”功能的景观类型的破碎化对应侵蚀的减弱, 反之则侵蚀加强；“汇”功能景观类型的破碎化意味着侵蚀的加强, 反之则侵蚀减弱。由此, 在类型水平上, 对“源”、“汇”两类景观类型, 表征破碎化的景观指数与侵蚀产沙的相关系数的符号应该相反。将林地、灌丛、草地、水域、坝地、耕地、坝地、居民地和裸土按“源”、“汇”功能归并为两类后计算景观指数。将这些指数与侵蚀模数做相关分析, 得到图 9所示的相关系数。对“源”、“汇”两类景观类型, ED、LPI、LSI和SPLIT与侵蚀模数的相关系数符号相同, 不满足正确指示景观变化的土壤侵蚀效应的条件。PD、PLADJ、PLAND和COHESION与土壤侵蚀模数的相关系数满足符号相反的条件, 可作为指示土壤侵蚀产沙的指数。如图 9所示, 对按“源”、“汇”功能归并后的两类景观类型, LPI、PD、PLAND、SPLIT、COHESON与流域输沙量的相关系数符号相反, 能指示景观格局变化对流域输沙量的影响；而对“源”、“汇”两类景观类型, ED、LSI、PLADJ与流域输沙量的相关系数符号相同, 不符合指示景观格局变化对输沙量的影响。

4.1.2 景观水平景观指数

在景观水平上, 景观指数与土壤侵蚀的相关系数整体较高(图 10)。从指数的计算原理上看, 景观水平指数的变化并不能表明侵蚀到底是增强还是减弱。植被斑块的破碎化是侵蚀增强的重要原因^[28-30]。如果是由于对土壤侵蚀起抑制作用的“汇”斑块类型的扩张、“源”斑块类型趋于破碎化引起景观破碎度的降低, 则侵蚀减弱；反之, 如果破碎度降低是由“源”斑块类型的连结成片而导致的, 则侵蚀加强。例如斑块密度PD与侵蚀产沙模数和流域输沙量都呈较好的相关性, 但是PD的变化既可能是“源”类型斑块破碎化导致的, 也可能是“汇”斑块破碎化引起, 或者是二者的共同变化形成的。同样, 在景观水平上, ED、PLADJ、COHESION、SPLIT、SHDI指数都存在分不清是“源”变化还是“汇”变化引起指数值变化的问题。在景观水平上, 尽管景观指数与土壤侵蚀表征变量之间可存在很好的相关关系, 但难以正确解释侵蚀的变化。

4.2 影响景观指数有效性的原因

景观指数在表达景观格局与土壤侵蚀之间的关系时的有效性取决于指数的性质和侵蚀的过程特征^[7]。在类型水平上景观指数指示景观格局与土壤侵蚀的关系的有效性要好于景观水平。从指数的性质来看, 斑块、边界密度和形状类的指数都是以斑块边界为基础建立的。在划分为“源”“汇”两个功能相反的斑块类型的景观中, “源”的边界同时也是“汇”的边界, 因此二者的变化方向一致, 无法正确反映景观格局变化对平均土壤侵蚀模数的影响。其中的LPI是描述最大斑块的指标, 反映的是景观中的极端斑块值, 在指示景观格局变化方面具有很大的不稳定性。植被斑块的空间集聚和分散直接决定径流和泥沙“源”区之间的连通性^[30-31], 因此描述景观类型聚集与分散的指数比其他类型的指数在指示景观格局的侵蚀产沙效应时更有效。

降雨侵蚀的产沙过程和输沙过程都是沿水流路径方向进行的, 具有方向性。上游景观类型单元的组成和格局影响径流量、流速等与径流侵蚀力和挟沙能力密切相关的指标, 因而是立地侵蚀产沙和输沙的决定性因素之一。在坝地等地形洼地和植被斑块、植被带等景观单元的拦截、阻碍下, 流域内产生的泥沙能够到达出口的比例很小。只有与河道之间存在畅通泥沙输送通道的区域(有效集水区)内的泥沙才有可能对流域出口泥沙量有贡献, 大部分区域并不是泥沙输出的有效贡献区域^[32]。由于泥沙的输移是沿汇流路径方向的, 从分水岭经坡面、沟道到流域出口是一个汇聚的过程, 淤地坝、植被带等一些位于关键位置的少数几个“汇”景观类型单元往往控制了流域输沙量。而基于斑块-廊道-基底模式的景观指数是各向同性的, 不能表达景观类型单元沿水流路径方向的空间排列, 因而在反映景观内泥沙“源”“汇”单元的空间分布对泥沙输出的控制作用时有效性不足。鉴于此, 在土壤侵蚀研究中将景观指数作为比较和分类的工具而非预测模型更为合适^[33]。近年来, 针对土壤侵蚀发展具有指示、预测功能的景观指数成为一个新的研究热点。在国内, 傅伯杰等^[8]提出了从坡面到流域的多尺度土壤侵蚀评价指数。研究^[34]表明, 多尺度土壤侵蚀评价指数比RUSLE中的C因子更能有效识别需要进行格局优化的流域。Ludwig等^[10]在水土保持功能二元区分的基础上, 构建了适用于规则径流小区的方向性渗透指数(DLI)。基于DLI的框架, 又进一步发展出植被覆盖度、地形为参数, 适用于从小区、坡面到流域和区域的渗透指数(LI)^[35-36]。Mayor等^[37]在将坡面划分为水土流失的“汇”和“源”的基础上, 提出了“源”到“汇”之间的水流路径长度来表征景观水土保持能力。而Borselli等^[38]则更进一步, 将上坡集水区面积、植被和地形因素确定的侵蚀产沙可能性与泥沙传输路径的植被、地形等特征确定的泥沙输出可能性相结合, 直接将景观格局与侵蚀产沙和输沙风险结合起来。这些指数以显式或隐含的方式将土壤侵蚀过程、景观类型功能和景观类型的空间分布耦合起来, 耦合了景观格局和土壤侵蚀过程, 具有明确的意义和较高的有效性。

4.3 结果的不确定性

本文论证了“Fragstats家族”景观指数在土壤侵蚀应用中的有效性及其正确遴选和使用原则, 但由于本文采用的侵蚀模型主要模拟面蚀和细沟侵蚀^[17], 不具模拟沟道侵蚀的能力。而在研究流域所在的黄土丘陵沟壑区, 冲沟、切沟等沟道侵蚀是主要的泥沙来源^[39]。沟道网络的连通性等指标可能更适合耦合沟道侵蚀与景观格局。此外, 模型不确定性可能会对本文的结果有一定的影响, 主要表现在泥沙传输能力参数K_tc的率定上, 校验采用的是流域出口处监测的输沙数据, 监测数据本身就存在一定的误差。且对景观类型按水土保持功能进行了类型的合并, 会造成一定程度的信息遗漏, 导致一定的泥沙输移计算误差。但模拟中将降雨、地形、土壤属性都作为不变的参数处理, 因此这些信息遗漏不会导致结论的错误。在研究区, 坡面是流域面积比重最大的部分, 是主要的径流贡献区。而坡面景观类型配置是径流的重要影响因素。因此, 有效描述流域景观格局与坡面侵蚀相互作用的指数在一定程度上也反映了沟道侵蚀的强度。因此, 本研究揭示的景观指数的有效性程度可能会随流域沟壑密度的增加而有所改变, 但由此得出的遴选和使用原则仍具有普适性。

5 结论

用Fragstats家族景观指数解释景观格局变化的土壤侵蚀效应时, 在类型水平上有一定的有效性, 但在景观水平上景观指数变化不能正确解释景观格局变化对土壤侵蚀的影响。景观类型聚集与分散的指数比其它指数在指示景观格局的侵蚀产沙效应时更有效。在研究景观格局变化的土壤侵蚀效应时, 有效的景观指数需要同时满足两个条件：(1) 对“源”、“汇”两类景观类型, 景观指数与土壤侵蚀状况表征变量的相关系数符号相反；(2) 对同为“源”或“汇”景观类型的多个景观类型, 景观指数与土壤侵蚀表征变量的相关系数应具有符号一致性。尽管景观指数在类型水平上具有一定的有效性, 但用其预测景观格局变化的侵蚀效应时存在很大的不确定性。因此, 通过景观格局分析将景观格局与土壤侵蚀过程耦合起来需要发展基于侵蚀过程机制、耦合格局和过程的景观指数。