延河流域产沙过程对“源”“汇”景观格局的响应

侯 琨，王秀茹†，王计平，袁普金，程 复

(1.北京林业大学水土保持学院，100083，北京；2.中国林业科学研究院国家林业和草原局盐碱地研究中心，100091，北京；3.水利部水土保持监测中心，100053，北京)

由景观要素的空间分布和组合特征构成的景观格局与生态过程的相互联系是开展景观生态学研究的重要切入点[1]。格局决定过程，过程也会推进格局的演进。作为一种典型的面状物质能量流动过程，水土流失的发生和发展与景观格局的演变密不可分[2-3]。在人类干扰较强的景观格局中，剧烈波动变化的土地利用/覆被情况是景观格局演进的直接表现[4]，其引起的景观结构和功能变化，可对水土流失等生态过程产生重大影响[5]。

从景观类型角度，土地利用作用于水土流失过程主要体现在用地类型改变、用地方式差异、产沙源相对入河口的距离等方面[6]。研究表明，要达到切实有效控制土壤侵蚀的目的，除了改变景观类型，还要实现景观类型在空间上的合理配置[7-8]。异质景观的空间分布格局作用于水文结构和产沙系统，会造成景观破碎化程度、水文连通性等指标的改变，增强或削弱景观斑块产生或拦截水土流失的能力，最终直观表现在河流断面泥沙输出数量的变化[9]。

当前，关于景观格局与水土流失过程关系的研究大多数是基于小区或坡面试验[10]，集中于单一或少数几种土地利用类型、植被分布等因素对侵蚀产沙过程的独立或综合影响等[11]，并在此基础上推演更大尺度上的水沙过程，而从流域、区域尺度考虑景观空间配置因素开展的研究稍显薄弱。随着尺度上升至流域层面，流域内部地表径流和泥沙的产生机制及运移过程将会发生明显改变，在大尺度上应用有限小尺度模拟预测得出的模型可能会造成信息遗漏和丢失[12]；因此，在流域尺度上引入“源”“汇”理论，探究景观格局的时空分布特征变化与水沙过程的定量关系，可为揭示有利于流域可持续发展的合理景观格局配置提供新的角度和理论依据。

笔者以位于黄河多沙粗沙国家级水土流失重点治理区的延河流域作为研究对象，基于遥感与地理信息系统等技术，识别流域“源”“汇”景观类型，采取水文统计分析、景观指数提取、建立多元回归模型等方法，构建景观格局与平均输沙率的量化模型，以期为延河流域未来的生态恢复和流域科学治理工作探索新的方向。

1 研究区概况

延河是黄河中游的一级支流，地理坐标范围大致为E 108°37′48″～110°28′48″、N 36°21′～37°19′12″，全长约286.9 km，流域面积约7 725 km2[13]，95%为黄土丘陵沟壑地貌，土地利用方向主要为草地、林地和耕地。多年平均降水量514 mm，降水多集中在6—9月。土壤类型以黄绵土为主；植被类型主要为森林草原植被，分带明显。

为反映该流域的真实情况，严格控制变量，笔者针对甘谷驿水文站控制断面以上区域(以下简称“研究区”)开展工作；该水文站位于延安市宝塔区，为延河流域的把口水文站，建站时间早，长时间序列资料较完整，控制流域面积约5 891 km2。

2 数据与方法

2.1 数据来源

土地利用等地理信息数据来源于中国黄土高原资源与环境遥感系列图编委会以及中国科学院生态环境研究中心；甘谷驿站水文泥沙数据来源于水利部水文局(1969—1996年)和黄河水利委员会中游水文水资源局(1997—2019年)存档的逐日实测资料；降水量数据通过向国家气象科学数据中心申请获得，具体站点为位于研究区中部区域的安塞气象站(编号53841，地理位置E 109°19′12″、N 36°52′48″)。

2.2 研究方法

2.2.1 泥沙过程时变趋势 以年为单位对1969—2019年的逐日输沙实测资料进行整理统计，采用Mann-Kendall趋势检验法判断其变化趋势[14]，采用单位降水距平累积曲线法判断降水因素对研究区泥沙过程的影响程度[15]。

2.2.2 “源”“汇”景观类型分类 以GB/T 21010—2017《土地利用现状分类》为参考，结合研究区实际情况及本研究的需要，将土地利用类型划分为耕地、林地、草地、建设用地、水域和其他用地。根据“源”和“汇”理论，将耕地、建设用地和其他用地归类为“源”景观类型，将林地、草地和水域归类为“汇”景观类型(受遥感影像分辨率所限，将园地纳入林地一并考虑)。

2.2.3 景观指数筛选及提取 景观指数是体现景观几何特性和配置结构的重要指标，本次分别从类型水平和景观水平选取部分代表性较强的常见指数(表1)。景观指数基于景观结构数量化软件Fragstats 4.2进行提取，指数计算公式及其生态学意义详见文献[16]。

表1 选取的景观指数及尺度Tab.1 Selected landscape indexes and scales

2.2.4 线性相关分析 对于固定面积的流域，“源”“汇”景观在理论上应呈现出相反的变化趋势，即在类型景观水平上，“源”“汇”景观指数与流域产沙量若分别表现为反向相关，可认为是有效指标。在类型水平上，对同属于“源”或“汇”景观类型的各类景观要素，若景观指数与平均输沙率均表现为同向相关，则可认为是有效指标[17]，纳入多元回归模型的解释变量集合。

景观格局指数与产沙过程的相关性借助皮尔逊简单相关系数R进行判断

(1)

式中：X、Y为变量；Cov(X,Y)为X与Y的协方差；Var[X]为X的方差；Var[Y]为Y的方差。

2.2.5 多元回归模型 通过构建多元回归模型，确定景观格局指数与平均输沙率的定量相关关系。模型分析基于SPSS Statistics 21软件包进行，一般如式(2)所示：

y=β0+β1x1+β2x2+…+βpxp+ε。

(2)

式中：y为响应变量；x1，x2，…，xp为解释变量；β0为回归常数；β1，β2，…，βp为偏相关系数；ε为随机误差。

3 结果与分析

3.1 泥沙过程时变特征

3.1.1 泥沙变化趋势性 1969—2019年间，研究区年及汛期(6—9月)平均输沙率分别为996.58和2 930.92 kg/s(图1)，整体波动幅度较大，最大值(1977年)与最小值(2015年)相差约276倍，汛期输沙量占全年比例达到97.47%，年内及年际差异均极其显著。

图1 1969—2019年及汛期平均输沙率变化趋势Fig.1 Variation trend of annual and flood season average sediment transport rate from 1969 to 2019

年及汛期平均输沙率的Mann-Kendall趋势检验结果表明，统计量Z值分别为-4.76和-9.03，说明研究区年及汛期平均输沙率均呈现下降趋势，而|Z|>Z1-α/2=2.58，表明其通过了给定显著性水平α为0.01的检验，下降趋势极其显著。

3.1.2 泥沙变化阶段性 为判断降水因素对泥沙过程的影响程度，分别绘制输沙率(图2)和单位降水输沙率的距平累积曲线(图3)。通过对比可以看出，2条曲线均存在一定程度的波动，但波动规律和趋势几近一致，说明降水因素并不是造成研究区产沙输沙量波动变化的主要因素。

图2 年平均输沙率距平累积曲线Fig.2 Cumulative departure curve of annual sediment transport rate

图3 单位降水输沙率距平累积曲线Fig.3 Cumulative departure curve of annual sediment transport rate by unit precipitation

3.2 “源”“汇”景观的时变趋势

3.2.1 “源”“汇”景观类型的时变趋势 1985—2015年，研究区“源”“汇”景观各类型以2000年为界呈现较为明显的分时段变化趋势(图4)：1985—2000年，景观结构基本维持稳定，耕地、草地是优势景观要素；2000年之后，“源”“汇”景观类型结构呈现显著变化，草地、林地面积大幅增加，逐渐成为研究区内的优势景观要素，耕地面积急剧减少。

图4 1985—2015年各类“源”“汇”景观比例Fig.4 Percentage of the various “source” and “sink” landscape from 1985 to 2015

3.2.2 “源”“汇”景观整体时变趋势 整体上看，研究时段内研究区“汇”景观面积比例始终高于“源”景观(图5和图6)，其中，1985—2000年，“源”景观与“汇”景观分布呈现相对均势且较为稳定，“汇”景观的面积比例在55%左右浮动，略高于“源”景观；2000年之后，“源”景观与“汇”景观的均势分布状态开始发生倾斜，“汇”景观面积大幅增加，至2015年，“汇”景观面积比例达到91.0%。

图5 1985—2015年 “源”“汇”景观比例Fig.5 Percentage of the “source” and “sink” landscape from 1985 to 2015

图6 1985—2015年研究区“源”“汇”景观结构Fig.6 “Source” and “sink” landscape structure of the study region from 1985 to 2015

3.3 “源”“汇”景观格局变化对泥沙过程的影响

3.3.1 最优响应时段的选取 鉴于遥感影像解译获取的景观结构信息仅能代表研究区一段时间内相对稳定且分布平均的景观格局特征，而景观格局结构发生变化之后需要经历一段时间才能在水土资源循环系统中发挥作用，因此，研究以生态建设工程为主要驱动力的“源”“汇”景观格局变化与泥沙过程的关系，需要寻求最优响应时段。

如表2所示，不同年份耕地面积和林草地面积与若干不同时间段内的年均输沙率均有较强的相关性，其中，对应年份以来3年的平均输沙率与耕地面积、林草地面积相关性最为显著，由此确定本研究的最优响应时段为获取土地利用信息对应年份以来3年。

表2 不同年份耕地面积(X)、林草地面积(Y)与对应年份以来不同时段平均输沙率(Zi)决定系数Tab.2 Determination coefficient between cultivated land area (X), forest and grassland area (Y) in different years and average sediment transport rate (Zi) in different periods since corresponding year

3.3.2 泥沙过程对类型水平景观指数变化的响应 以1985、1995、2000、2008和2015年各“源”“汇”景观斑块类型以及合并后的“源”“汇”景观的类型水平景观指数为解释变量，以最优响应时段3年的平均输沙率为响应变量，分别计算皮尔逊简单相关系数，判断景观指数指示流域泥沙过程的有效性。

对于将各类景观类型合并后的“源”、“汇”景观，除了景观形状指数，其余参与统计的8个类型水平景观指数均呈现出反向相关的态势(图7)，表明“源”“汇”景观整体类型水平景观指数的波动变化作用于泥沙过程的机制较为复杂，需进一步分析各种“源”“汇”景观类型水平景观指数的指示能力。

图7 “源”“汇”景观类型水平景观指数与平均输沙率相关性Fig.7 Correlation between the class level index of the “source” and “sink” landscape and average sediment transport rate

对于各类“源”景观斑块类型，斑块密度、平均斑块面积、相似临近比例、聚合度与年均输沙率表现出同向相关性(图8)，这4个指标可有效指示“源”景观斑块变化对流域产沙量的影响。对于各类“汇”景观斑块类型，斑块密度、平均临近指数与年均输沙率表现出同向相关性(图9)，均为负相关，这4个指标可有效指示“汇”景观斑块变化对流域产沙量的影响。

图8 各类“源”景观类型水平景观指数与平均输沙率相关性Fig.8 Correlation between the class level index of the various “source” landscape and average sediment transport rate

图9 各类“汇”景观类型水平景观指数与平均输沙率相关性Fig.9 Correlation between the class level index of the various “sink” landscape and average sediment transport rate

3.3.3 泥沙过程对景观水平景观指数变化的响应 将最优响应时段的平均输沙率与研究区景观水平景观指数进行相关分析(图10)，结果显示，个别指数表现出较强的相关性，但相关系数整体较类型水平景观指数偏低，指示能力稍逊。

图10 景观水平景观指数与平均输沙率相关性Fig.10 Correlation between the class level landscape index and average sediment transport rate

此外，景观水平景观指数的衡量对象是流域整体景观，“源”景观和“汇”景观的类型组成和结构特征变化均可引起景观水平景观指数的波动，因而难以识别对流域产沙量变化起决定作用的景观指数。

3.4 “源”“汇”景观格局与泥沙过程的量化关系

分别以有效的“源”“汇”景观类型水平景观指数为解释变量，以最优响应时段的平均输沙率为响应变量，建立多元线性回归模型，定量分析景观指数对泥沙过程的指示程度。

对于“源”景观，提出假设模型为

y=β0+β1x1+β2x2+β3x3+β4x4+ε。

(3)

式中：y为最优响应时段的平均输沙率，kg/s；x为“源”景观类型水平景观指数，量纲为1，其中，x1为PD，x2为AREA_MN，x3为PLADJ，x4为AI。

在SPSS Statistics 21中输入数据进行多元回归分析，初步回归结果显示：x3变量被排除；x1与x4的方差膨胀因子(variance inflation factor，VIF)远大于10，存在显著的多重共线性，剔除x1，x4参与进一步回归。将保留的x2和x4与y再次进行回归，结果为

y=2 584.923+7.164x2-32.207x4+ε。

(4)

(2 878.665) (1.696) (35.008)

n=5R2=0.939F=15.304

回归方程复相关系数为0.969，拟合优度达到0.939，F检验在0.1的水平上显著，由此判断，回归方程整体较显著，在“源”景观类型中，平均斑块面积均与年均输沙率正相关，聚合度与年均输沙率负相关。研究时段内，“源”景观类型的平均斑块面积大幅减小，聚合度略有波动减小，这与多元回归的结果并不完全一致，表明2个因素并非独立作用于泥沙过程，而是共同影响年均输沙率的波动变化，平均斑块面积的影响程度相较更高。

对于“汇”景观，采取相同的分析步骤，提出假设模型为：

y=β0+β1x1+β2x2+ε。

(5)

式中：y为最优响应时段的平均输沙率，kg/s；x为“汇”景观类型水平景观指数，量纲为1，其中，x1为PD；x2为PROX_MN。回归结果如下：

y=371.439+5 765.753x1-0.001x2+ε。

(6)

(291.113) (2 130.384) (0.000)

n=5R2=0.945F=17.029

回归方程复相关系数为0.972，拟合优度达到0.945，F检验在0.1的水平上显著，由此判断，回归方程整体较显著，在“汇”景观类型中，斑块密度与年均输沙率正相关，平均临近指数与年均输沙率负相关。研究时段内，“汇”景观类型的斑块密度略有波动减小，平均临近指数成倍大幅增加，表明“汇”景观斑块破碎化程度降低，斑块间邻近度逐渐升高，景观连接性趋好，拦蓄泥沙的能力增强，利于削弱泥沙的产生和汇集过程，减少入河泥沙量。

4 结论

1)1969—2019年间，研究区汛期输沙量占全年比例超过97%；年及汛期平均输沙率通过了给定显著性水平α为0.01的Mann-Kendall趋势检验，下降趋势极其显著。对比输沙率和单位降水输沙率的距平累积曲线，排除了降水因素主导研究区产沙输沙量波动变化的可能性。

2)1985—2019年间，研究区“汇”景观逐渐占据主导格局，引起这种波动变化的主因是退耕还林还草等生态建设工程导致耕地面积大幅减少，林、草地面积大幅增加。在提出最优响应时段的基础上，通过计算皮尔逊简单相关系数，选取可有效指示研究区泥沙过程的景观指数，建立多元回归模型。结果表明：“源”景观类型的平均斑块面积和聚合度共同影响年均输沙率的波动，“汇”景观类型的斑块密度和平均临近指数共同影响年均输沙率的波动，模型拟合优度分别达到0.939和0.945，影响作用较为显著。

对于景观水平景观指数，相关系数整体较类型水平景观指数偏低，且无法区分引起指数波动的主导要素是“源”景观类型还是“汇”景观类型，因而难以识别对流域产沙量变化起决定作用的景观指数。

3)多元回归的结果表明，影响泥沙过程的景观格局要素主要为景观的破碎化程度和景观的聚合连接性，当“源”景观的破碎化程度降低，“汇”景观的邻近度升高、连接性增强时，形成的景观格局结构利于流域尺度水文循环系统的泥沙拦蓄和调控。