山丘区中小河流滩地时空演化与扰动因子分析

王金平，夏继红，汪颖俊，窦传彬，周之悦，叶继兵，刘 瀚

(1.河海大学水利水电学院，南京 210098；2.龙游县水利局，浙江 龙游 324400)

0 引 言

河流滩地作为河道水域与陆域的过渡区域，具有调蓄洪水、稳固河道等功能，对河流生态健康保护具有重要意义[1-3]。近年来，在对河流资源开发利用时，出现河流滩地面积减小、破碎化程度增加等现象引起不少学者的关注。例如：Crosato[4]基于河流形态二阶线性物理模型，以河道中滩地发育的数量，来预测蜿蜒型河道和辫状型河道的相互转换；Defina[5]基于无黏性河道二维形态动力演变有限元模型，模拟自由沙洲的形成和演变过程；李志威等[6]基于长江、汉江等7条河流不同河段的沙洲形态，将沙洲形态概化为椭圆形、竹叶形和镰刀形3种典型沙洲，并分析典型沙洲的平面形态与河道的关系；高进[7]对沙洲形状的概化，提出沙洲的发育长度公式；孙昭华等[8]通过顺直和分汊两种河型交界位置的洲滩和浅滩动态特征分析滩地形态演化；伊紫函等[9]通过形态指数定量划分滩地类型，分析不同类型滩地的时空演化等。然而，现有研究大多针对大型骨干河流单一滩地类型和形态特征，而对中小河流多滩地时空演化及扰动分析的定量研究则相对较少。本文以浙江省龙游县灵山港滩地为研究对象，应用景观格局指数分析分析滩地的时空演化，并在计算扰动度的基础上，进一步探讨主要扰动因子，以期为河道治理规划提供参考。

1 研究区概况

灵山港位于浙江省西部的仙霞岭山系东北缘，介于 E119°07′～119°12′，N28°49′～29°03′之间，是钱塘江上游衢江右岸的一条支流，发源于海拔1 265 m的遂昌县高坪乡和尚岭，流域总面积726. 9 km2。灵山港流域属亚热带季风气候区，雨量丰沛，四季明显。灵山港在龙游县境内的流域面积为367.6 km2，主流长43.79 km，河道平均比降为2.45‰。河道自南往北流经沐尘、溪口、灵山、官潭、寺后、上圩头、龙游县城，在驿前村附近注入衢江。河道上游沐尘村处建有大(二) 型水库——沐尘水库，流域内年平均降水量为1 815 mm，年平均径流总量为4.49 亿m3，多年平均径流深1 037.3 mm，径流系数为0.62。该河道滩地资源较为丰富，为维护河流健康发挥了重要作用。但是近年来，由于河流资源过度的开发利用，尤其是河流滩地采砂，滩地的结构稳定性遭到破坏，在水流冲刷的作用下滩地总体面积减小，破碎化程度增高，空间分布发生明显变化，河流生态功能退化严重。

2 材料与方法

2.1 数据获取方法

本文分别选取1982年、2003年龙游县地形图(1∶10 000)以及2010、2013年Google Earth影像图(1∶10 000)均用AutoCAD处理图像，并读取灵山港不同位置的河宽、主槽宽、河长等基本参数，读取河道沿线各滩地的面积、周长、纵径、横径、间距等几何参数以及各滩地的位置信息。通过对龙游县灵山港的现场勘测及相关资料整理，灵山港上游沐尘村建有沐尘水库，沿线建有16座堰且两岸有不同形式的护坡。在1982-2013年间，灵山港河岸线和河道走势未发生明显变化，河道滩地变化明显，此次选取的图像资料均是在10-11月，为枯水期；确定不同典型年份区间的扰动因子：水力坡度、弯曲度、主槽宽/河宽、平均堰高、水库距离和采砂强度。

2.2 数据处理与分析方法

将灵山港研究区域按水力坡降划分为上游、中游和下游，分别计算流域和区段的滩地景观格局指数，借鉴土地利用动态度模型，引入滩地扰动度。采用灰色关联分析法计算区段间在1982-2003年、2003-2010年、2010-2013年扰动度与各扰动因子之间的关联度，进而确定不同时间段内，扰动因子对滩地演化的贡献度。

2.2.1 景观格局指数计算方法

利用景观格局指数可以定量描述景观格局分布，结合景观格局与景观过程在时间维上的联系，来更好地阐述景观空间格局动态演化过程[10]。龙游县灵山港滩地属于典型的景观斑块，本文主要计算滩地密度PD、平均形态分形维数FD、平均邻近度指数MPI等滩地景观格局指数[10,11]。

PD=n/A

式中：A、n分别为研究区滩地总面积、数量；Pi和Ai分别为滩地i的周长、面积；h为滩地i与其相邻滩地的最小邻近距离。

2.2.2 扰动度计算方法

为寻求能反映滩地的面积、数量、边界形状及空间关系的变化的指数，本文借鉴土地利用动态度模型[12]，引入河流滩地扰动度(DI)。该指数以PD、FD、MPI三个景观格局指数相对变化值乘以一个权重系数，最后进行求和即可，扰动度大小反映滩地的面积、数量、边界形状及空间关系的变化程度，计算公式如下。该指数能反映出河流滩地所受扰动的程度大小。

(1)

式中：DI为扰动度指数；PDa和PDb分别为研究开始时期及结束时期滩地斑块密度；FDa和FDb分别为研究开始时期及结束时期斑块分形维数；MPIa和MPIb分别为研究开始时期及结束时期邻近度指数；T为研究时间范围；α、β和γ为权重系数。

2.2.3 影响因子分析方法

本文采用灰色关联分析法测1982-2003，2003-2010，2010-2013年滩地扰动度变量与扰动因子变量间的关联度，并进行综合比较和排序，确定各扰动因子在滩地时空演化的贡献度。灰色关联分析法是通过计算灰色关联度来确定系统因素间的影响程度或因子对系统主要行为的贡献测度的一种方法，其基本思想是根据序列曲线几何形状的相似程度来判断其联系是否紧密[13,14]。它对样本量的多少和样本有无规律都同样适用，不会出现量化结果与定性分析不符合的情况。灰色关联分析方法计算公式如下：

X0=[x0(1),x0(2),…,x0(n)]

X1=[x1(1),x1(2),…,x1(n)]

Xm=[xm(1),xm(2),…,xm(n)]

(2)

(3)

式中：X0和Xi分别为系统行为特征序列和相关因素序列；ξi为相关因素序列Xi的第k个元素与系统行为特征序列X0的第k个元素间的关联系数；ri为关联度；ρ为分辨率，用来减小Δ(max)过大而使关联系数失真的影响，人为引入此系数来提高关联系数之间的显著性；参数min和max的取值通过对各个比较序列不同时刻的绝对差值的比较来确定。

由于灰色关联分析系统中各数列数据量纲不同，不便比较或难以得到正确的结论，为保证结果的可靠性，在进行灰色关联分析时，需先采用均值变换对数据进行无量纲处理。

3 结果与分析

3.1 滩地演化分析

3.1.1 总体演化特征

在1982-2013年间，灵山港滩地数量明显增加，面积大幅度减少。灵山港滩地数量从1982年的78个增加到2013年的112个，滩地数量增加43.59%。滩地总面积呈现不断下降的趋势，从1982年的208.29 万m3下降到2013年的120.88 万m3，滩地总面积减少了41.97%。这说明从1982-2013的31年间，灵山港滩地资源严重萎缩，滩地对河道蓄洪、净化水质、生物栖息等生态功能在减弱。

3.1.2 格局演化特征

斑块密度反映景观破碎化程度，斑块密度越大，景观破碎化程度越高，反之则越低。从流域尺度分析，灵山港滩地斑块密度从1982年的0.374一直增加至2013年的0.960，说明1982-2013年，龙游县灵山港滩地呈现破碎化演变，空间分布不均匀性在逐渐增加。从区域尺度分析，由图1可得各区段滩地演化各有不同：上游段滩地斑块密度增加主要在2003-2010年间，滩地斑块密度由0.548(2003年)增加至0.936(2010年)，中游段滩地斑块密度一直在增加，滩地斑块密度由0.384(1982年)增加至1.236(2013年)，下游段滩地斑块密度从0.252(1982年)增加至0.799(2010年)，在2010-2013年间无明显增加。从统计的基础数据分析，上游段滩地破碎化表现为滩地数量增加，滩地的面积有轻微的减少；下游段滩地破碎化主要表现在滩地面积的急剧减少，中游段滩地破碎化在滩地数量的增加同时伴随着滩地面积减少。

图1 斑块密度演化特征Fig.1 Patch density evolution characteristics

形态分维数反映了景观斑块的形状复杂程度，分形维数越大，滩地形态变得越复杂、越不规则，在空间上呈扩张态势，反之则相反。从流域尺度分析，灵山港滩地形态分维数由1982年的1.118增加至2010年的1.151、再下降到2013年的1.142，总体呈现上升趋势，这表明滩地外围形态呈现复杂、不规则化发展，在空间上呈窄长化。从区段尺度分析，由图2可得上游段滩地形态分形维数在1982-2003年间减小，在2003-2013年间增加；中、下游段滩地形态分维数在1982-2010年间增加，在2010-2013年间减小，中游区段的形态分维数变化更显著。各区段滩地形态分维数在时间维度上虽然有减少的情况出现，但是各区段滩地的形态分形维数总体呈现增加的趋势。灵山港滩地外围形态呈现不规则、窄长化发展，尤其是2003-2010年间更加明显。2003-2010年间，对灵山港开发、建设与利用过程中，如河道滩地采砂、疏浚河道等，增加了滩地外围形态的不规则化，同时破坏滩地整体尤其是滩地边坡的结构稳定性，滩地外围形态在河道水流冲刷下呈现不规则、窄长化发展。

图2 形态分形维数演化特征Fig.2 fractal dimension evolution characteristics

平均邻近度指数既可以度量斑块的分离度也可以度量给定空间范围内斑块的破碎度，斑块面积越小且稀疏分布时，平均邻近度指数的值越小，反之斑块面积越大且聚集分布时，平均邻近度指数的值越大。从流域尺度分析，灵山港滩地平均邻近度指数由1982年的197.68增加到2013年的672.46，又滩地数量的增加和滩地面积减少，这说明滩地平均邻近度指数增加并非是滩地的分离度减小而是由于滩地破碎所致，滩地破碎后滩地的面积虽然在减小，但是滩地最小邻近距离减小使得平均邻近度指数在增加。从区段尺度分析，由图3可知上游段滩地平均邻近度指数先增加、后减小再增加的现象，中游段滩地平均邻近度指数一直在增加，下游段滩地平均邻近度指数在1982-2003年间减小，在2003-2013年间增大。在1982-2003年间，灵山港上游段滩地的面积和数量有少量的减少，但是在姜家至梅村的右岸由于泥沙淤积形成一个面积为0.135 km2，长轴长为1.77 km的滩地。下游区段平均邻近度指数在1982-2003年减小，主要是滩地面积减少使得在空间上表现稀疏分布。其余各区段滩地平均邻近度指数随时间呈现增加趋势，主要是由滩地破碎增加所致，而非滩地聚集分布趋势。

图3 平均邻近度指数演化特征Fig.3 Mean proximity evolution characteristics

通过分析1982-2013年灵山港滩地的斑块密度、形态分维数等景观格局指数的演化特征，灵山港滩地格局变化明显且在2003-2010年变化更显著，从上游到下游，滩地呈现窄长化的变化趋势，滩地分布均匀性减弱，破碎程度呈增大趋势。

3.2 影响因子分析

3.2.1 扰动度变化分析

灵山港滩地区段扰动度随时间上均呈现增加趋势，在2003-2010年时各区段扰动度均最大，1982-2003年间最小。由图4可知，在1982-2003年间，由于对灵山港的河流资源的开发在其自然承载力范围内，此时扰动度DI介于0.16～0.59；2003- 2010年间灵山港新建沐尘水库、桐子堰、橡胶坝等水利设施，全线河道在疏浚过程中且存在河流滩地大范围的采砂，此时扰动度DI介于1.53～3.25；2010年以后灵山港全面禁止采砂，但是由于之前的对滩地采砂，使得原有的滩地结构稳定性受到破化，在2010-2013年间的滩地的扰动度依旧较大，此时扰动度DI介于0.65～1.26。从1982-2013年，灵山港滩地数量与面积变化的主要原因是水利设施建设以及对河流滩地采砂。结合灵山港的实际，河流滩地扰动度能较好的反映灵山港滩地所受扰动程度。

图4 扰动度指数变化Fig.4 Disturbance index changes

3.2.2 影响因子分析

关联系数反映了两比较序列在某一时刻的紧密程度，关联度则反映了两比较序列计算时段内的平均紧密程度。在滩地演化过程中，关联度越大，表明该扰动因子对滩地格局演化的影响越大，反之则越小。由表1可得：研究区在1982-2003年间滩地扰动度关联系数由大到小排序，弯曲度(2.002)>主槽宽/河宽(1.861)>平均堰高(1.725)>水力坡度(1.546)；在2003-2010年间滩地扰动度关联系数由大到小排序，采砂强度(2.846)>平均堰高(2.811)>主槽宽/河宽(2.633)>弯曲度(2.464)>水力坡度(2.288)>水库距离(1.413)；在2010-2013年间滩地扰动度关联系数由大到小排序，主槽宽/河宽(2.787)>弯曲度(2.699)>水力坡度(2.548)>平均堰高(2.506)>水库距离(1.479)。

表1 典型年份影响因子的关联系数Tab.1 Correlation coefficient of typical year period

通过滩地扰动度关联系数的大小可以得出，在不同时间段水库距离对滩地格局演化贡献度都很低，其他扰动因子则相差不大。因而选取弯曲度、主槽宽/河宽及水力坡降的扰动度关联系数均值和采砂强度、平均堰高的扰动度关联系数均值分别来反映自然和人为扰动因子对河流滩地演化的贡献度。在1982-2003年和2010-2013年这两个时间段扰动度DI≤1.26，自然扰动因子对河流滩地演化的贡献度大；在2003-2010年时扰动度DI≥1.53，人为扰动因子对河流滩地演化的贡献度大；通过DI为1.26和1.53的自然和人为扰动因子对滩地演化的贡献度采用线性拟合求出扰动度DI临界值为1.35。

4 结 语

(1)1982-2013年间，灵山港滩地数量增加43.59%，滩地总面积减少41.97%，表明灵山港滩地发生了较为严重的萎缩和退化。

(2)1982-2013年间，斑块密度、形态分维数、平均邻近度指数等景观格局指数均增加，表面灵山港滩地格局变化明显，且在2003-2010年变化更显著，从上游到下游，滩地呈现窄长化发展趋势，滩地分布均匀性减弱，破碎程度呈增大趋势。

(3)灵山港滩地各区段扰动度随时间上均呈现先增加后减少，各个区段扰动度均表现在2003-2010年间最大，1982-2003年间最小。扰动度DI<1.35时，自然扰动因子如弯曲度、主槽宽/河宽和水力坡度对滩地的格局演化贡献度比较大，人为扰动因子的贡献度则相对较小；扰动度DI≥1.35时，人为扰动因子如采砂强度、平均堰高对滩地格局演化贡献度比较大，自然扰动因子的贡献度则相对较小。

□