小水电梯级开发河段地貌单元识别与异质性分析

于子铖，张 晶，赵进勇，彭文启，丁 洋，王 琦

（中国水利水电科学研究院 水生态环境所，北京 100038）

1 研究背景

我国小水电的建设在不同历史时期，其定位、发展导向和实际效用也不尽相同。从建国后至改革开放前，我国小水电主要是为了解决无电缺电地区人民生活用电问题；从改革开放后到20世纪末，小水电主要解决贫困地区脱贫致富及农村地区电气化问题；进入21世纪后，小水电则主要是为缓解资源约束，改善生态环境，促进可再生能源利用。不同历史时期的小水电建设在促进经济社会发展、能源结构转型等方面发挥了重要作用［1］。但与此同时，部分流域、河流存在着无序、过度开发等问题，使河流生态系统遭到了严重的破坏。随着生态文明理念的深入，国家和相关部门高度重视小水电的相关工作，先后印发了 《水利部关于推进绿色小水电发展的指导意见》 《关于开展长江经济带小水电清理整改工作的意见》《关于进一步做好小水电分类整改工作的意见》等。我国小水电行业发展的历史跨度大、遗留问题多，部分地区在理解、执行偏差下将小水电整治变成不加区分的 “一刀切拆除”，给生态环境保护和经济社会发展均造成负面影响。此外，小水电对碳达峰、碳中和的贡献不可或缺，小水电整治亟须正确处理好生态环境保护、经济发展、社会稳定之间的关系［2］。2021年12月，7部委发布的 《关于进一步做好小水电分类整改工作的意见》中提出小水电的整改需深入细致地做好评估分类，科学开展综合评估，合理确定整改和修复目标，明确了小水电退出、整改、保留的条件。

小水电建设导致河流栖息地受到强烈扰动，破坏了河流地貌空间异质性。河流是复杂的系统，非生物和生物成分在不同的时空尺度上相互作用。不同空间尺度的生态系统之间形成嵌套层级结构，较小的时空层次嵌套在较大的时空尺度内［3］。董哲仁［4］将其分为流域、河流廊道、河段、地貌单元、微栖息地5个层级，分别对应不同的空间尺度。其中，地貌单元（100～102m）是水流和泥沙输送过程的物理表现形式［3］，是河段（102～104m）的组成部分，在这个尺度内能够反映河流生态系统的结构、功能和过程［4］。地貌单元由河流不同位置的不同过程组合形成（河源区以侵蚀过程为主，下游以沉积形式为主），在任一自然河段内，特定位置都有明显的特征，如弯曲河道内侧一般都有深潭等。地貌单元的自然多样性是河流栖息地适宜生存的重要决定因素，河段地形越复杂和多样化，地貌异质性越高［5］，提高地貌异质性是维持生物多样性的前提［6］。从地貌单元过程-形式关系出发，Fryirs等［5］将地貌单元分为河道中段的沉积、滨岸沉积单元等四类；Belletti等［3］将其分为主河槽与河漫滩单元两大类。对河流地貌单元进行分类，可以简化河流廊道中复杂的连续过程、减少描述河流地貌过程所需的时间与工作量［7-8］。

长江流域小水电的无序开发使部分河流出现减脱水现象［9-10］。小水电梯级开发更是加剧了栖息地的破碎化，破坏了自然生态系统的生物多样性［11］，对底栖生物、鱼类、爬行动物与植物群落等空间分布及组成形成一定的影响［12-13］。大宁河是长江的一级支流，同时又位于三峡库区，其生态环境效应一直是相关领域的研究热点，主要包括库区土地利用、生态环境因子、生态环境质量以及生态系统服务等方面的响应变化［14-15］。本文以大宁河巫溪段为例，依据实地调查观测与文献调研，基于地貌单元过程-形式关系，对其地貌单元进行分类与识别，并从地貌单元尺度，探讨分析小水电梯级开发对河流地貌异质性的影响，旨在从地貌角度为小水电拆除、整改、保留等分类策略以及生态修复具体措施提供科学依据。

2 研究区域与研究方法

2．1 研究区域 巫溪县位于长江上游地区、重庆东北部，在大巴山东段南麓。地处东经108°44′至109°58′，北纬31°14′至31°44′之间。巫溪县主要河流均属长江水系，包括大宁河、西溪河、柏杨河、后溪河、东溪河等。巫溪县地形以山地为主，山地占93%，最低海拔139．4 m，最高海拔2796．8 m。起伏变化的地貌使其在水电站建设上具有得天独厚的优势，因此巫溪县内的小水电数量众多。据统计，截止到2020年，巫溪县有130座小水电，1座中型水电站，其开发方式大都为引水式。

大宁河发源于大巴山南麓，一支为加担湾河，源头在巫溪县高楼乡，河长约39．6 km，控制流域面积292 km2。一支为汤家坝河，源头为巫溪县和平乡小龙洞，河长约31．1 km，控制流域面积194 km2。两源东流至中梁乡龙头嘴汇合后称西溪河，再东流至两河口，与东溪河汇合始称大宁河，转而南流至宁厂镇纳后溪河、至巫溪县城下首纳柏杨河，经龙溪、大昌等地至巫山县城注入长江。河流基本属峡谷型河道，河弯多，坡降大，平均坡降均在10%以上。西溪河干流上有2座梯级电站，分别为中梁一级电站、中梁三级电站，其中中梁二级电站未在干流上，其他还有白马穴、关庙、崇溪沟等3座小型水电站，其挡水坝均未在干流上。大宁河干流上有1座梯子口电站，其挡水坝未在干流上，具体情况见表1。

表1 西溪河、大宁河干流水电站统计

结合2018—2021年的实地调研数据，本次研究河段选取西溪河中梁一级大坝—大宁河巫溪境内段，总长为60．34 km。为分析大宁河地貌单元的布局与分布以及水电梯级开发对地貌异质性的影响，依据支流汇入、干流小水电空间布局（白马穴、关庙、崇溪沟、梯子口电站附近流量变化较小，中梁二级电站未在干流上，故均未作为分段节点）、河流坡降等，将其分为中梁一级大坝下游—中梁一级电站、中梁一级电站—中梁二级大坝、中梁二级大坝—中梁三级大坝、中梁三级大坝—中梁三级电站、中梁三级电站—两河口（东溪河汇入口）、两河口—后溪河汇入口、后溪河汇入口—柏杨河汇入口、柏杨河汇入口—大宁河巫溪段境内8个河段，分别用R1—R8描述，见图1，其中，R1—R4属于梯级开发直接影响河段，R5—R8属于较自然河段（流量全部下泄，干流已无挡水坝）。

图1 研究区域图

2．2 水文分析 基于巫溪水文站1972—2016年水文数据，利用累积距平法、滑动秩和检验法对水文数据进行突变检验，明确突变点，据此将水文序列分为变异前、变异后等不同阶段，利用IHA-RVA法分析变异前后各月份平均流量（对应IHA指标的前12个指标）改变程度。地貌单元由时空尺度定义与限制，水文改变程度越小，识别的地貌单元分布越具有代表性。巫溪水文站位置见图1。

式中：IAi为水文改变度；Ni为变异后水文系列统计的水文参数值实际落在流量管理目标范围内的年数；Ne变异后水文系列统计的水文参数值期望落在流量管理目标范围内的年数；r＝50%；NT为变异后水文序列的总年数。

2．3 实地调研与测量 依据Bigmap下载高分辨率遥感数据，利用大疆精灵4 PRO无人机、华测i70 RTK、D390测深仪、LS1260B流速仪等仪器进行空-地实地测量，获取研究河段的遥感影像、河床地形变化、水深流速分布等数据，利用标尺法、筛分法对底质状况进行调查，依据 《河流泥沙颗粒分析规程》（SL42-2010）进行底质分类。S1—S10为实地调研的典型河段，位于R1—R8河段减水、深浅相间、弯曲等典型区域，其航拍图见图2。

图2 典型河段无人机航拍图

2．4 河流地貌单元分类与识别

2．4．1 Fryirs和Brierley地貌单元分类体系 地貌单元是构成水系的基石，可依据形态、沉积物成分（底质）、边界表面（植被覆盖程度）和在河流中的空间位置来划分［5］。本文以Fryirs和Brierley地貌单元分类体系［5］作为依据，该体系基于河流地貌单元的过程-形式组合，将河流地貌单元归为4类：冲积与侵蚀的基岩和巨石单元、河道中段的沉积单元、滨岸的沉积单元、冲积与侵蚀的细颗粒单元。其中，冲积与侵蚀的基岩和巨石单元分为台阶跌水（Bedrock step）、跌水-深潭（Step-pool）、串联跌水（Cascade）、急滩（Rapid）、滑水（Run）等9种地貌单元，河道中段的沉积单元分为深潭（Pool）、浅滩（Riffle）、心滩（Medial bar）等13种地貌单元，滨岸的沉积单元包括侧向滩（Lateral bar）、点滩（Point bar）等15种地貌单元，冲积与侵蚀的细粒单元包括侧向侵蚀滩（Sculpted lateral bar）、冲刷潭（Scour pool）等5种地貌单元，分类体系见图3。

图3 Fryirs和Brierley地貌单元分类体系［5］

2．4．2 地貌单元类型确定及其空间分布 不同的地貌单元，其空间分布、几何形状与边界、形态特征、沉积物组成及其植被覆盖不一致，其过程-形式组合不同，这为识别地貌单元提供了基本原则［16］。如急滩是非常稳定、陡峭的山脊状序列，由不规则横肋中的巨石排列而成，在基岩受限的情况下，这些横肋部分或全部横跨河道。在高流量下，巨石在结构上重新排列，形成稳定的序列，与发散或汇聚的流动都不相关。通常，15%～50%的水流呈现超临界流动［16］。本文将野外观测、高分辨率遥感数据与每种类型地貌单元的过程-形式关联的理论知识联系起来，再结合沉积物的组成分析，对河流地貌单元进行识别。

2．5 基于地貌单元类型的地貌异质性分析方法 地貌单元作为河流嵌套层级结构中重要的一环，其类型、丰度的差异影响着河流的空间异质性，空间异质性是指系统特征在空间分布上的复杂性和变异性［17］。河流空间异质性越高，河流生态系统越稳定。基于前人的相关研究，结合研究区域实际概况，用地貌单元丰富度指数、地貌单元密度指数、地貌单元多样性指数3个指标对其进行量化评价［9，18］。

（1）地貌单元丰富度指数（R）。地貌单元丰富度指数（R）指在给定范围内观察到的地貌单元类型数与整体最大类型数相比的情况，表达式如下：

式中：NTGU为所研究区域内地貌单元类型总数（例如，在存在浅滩，深潭和急滩的情况下，NTGU＝3）；n为地貌单元类型的最大数量。

（2）地貌单元密度指数（D）。地貌单元密度指数（D）是指计算每单位长度调查范围内的地貌单元总数（与类型无关），表达式如下：

式中：NGU为沿调查的范围观察到的地貌单元总数（如，对于5个浅滩、5个深潭和5个边滩，NGU＝15）；L为所调查的范围（或子范围），km。

（3）地貌单元多样性指数（Hd）。Shannon多样性指数可用来反映景观异质性的程度，其值越大，景观异质性程度越好，景观内的生物越能稳定生存，景观异质性包含空间异质性。地貌单元多样性指数（Hd）是基于Shannon多样性指数提出的，表达式如下：

式中：m为研究区域地貌单元类型的总数；Pi为地貌单元类型i的面积（Ai）所占研究河段总面积（TAj）的比例。

3 结果分析

3．1 水文分析 基于巫溪水文站1972—2016年径流数据，分析径流量变化趋势。从图4可知，随着时间的变化，年径流量趋势变小，年径流量最大值出现在1983年，为40．82亿m3，年径流量最小值出现在在2006年，为10．83亿m3。

图4 径流量变化趋势

依据巫溪站年径流量数据，计算距平值，建立累积距平曲线，见图5。由图可见，年径流量距平曲线整体呈现一个先上升后下降的趋势，在1983年显著上升，在2005年之后表现为显著下降，1993年累计年径流量最大。据此，确定1983、1993、2005为可能变异点。经滑动秩和检验算法，确定1983年、2005年为可能变异点，见表2。

图5 累积距平曲线

表2 秩和检验计算结果表（1972—2016年）

基于上述方法，确定1983、2005年为巫溪水文站控制区域的水文变异点。因此，将水文序列分为1972—1982年（变异前）、1983—2004年（变异阶段Ⅰ）、2005—2016年（变异阶段Ⅱ）。利用IHA-RVA方法，分别计算变异阶段Ⅰ与变异前、变异阶段Ⅱ与变异前各月份平均流量改变程度，见图6。变异阶段Ⅰ与变异前相比，各月份水文改变度的平均值为15．53%，5、6月份水文改变度最小，为0，7月份改变程度最大，为43．75%；变异阶段Ⅱ与变异前相比，各月份水文改变度的平均值为24．42%，1、11月份改变程度最小，为8．33%，10月份改变程度最大，为47．62%，6月水文改变度为19．79%。

图6 水文改变度变化趋势

3．2 地貌单元空间分布

3．2．1 地貌单元初步识别 根据上述水文分析，选择水文改变程度相对较低的时段（6月份），在实测流量下进行地貌单元的识别，8个河段的实测流量分别为4．43、7．91、3．44、3．14、20、38．7、46．2、50．6 m3/s。基于Fryirs和Brierley的分类体系，对比不同地貌单元的过程—形式，结合实际调查与资料梳理结果，最终识别了串联跌水（Cascade）、急滩（Rapid）、滑水（Run）、平滑流动区（Sculpted run）、深潭（Pool）、浅滩（Riffle）、缓流区（Slack water）、边滩（Point bar、Lateral bar）、江心洲（Island）、心滩（Medial bar）等10个地貌单元。其中，平滑流动区与滑水的区别在于底质颗粒组成的粗细，细颗粒的为平滑流动区。S1—S3典型河段的分布如图7所示。

图7 S1—S3地貌单元空间分布

串联跌水具备在陡峭、基岩封闭的环境中观察到的非常稳定的、粗粒的（或基岩）特征，水流在大石块上层层叠叠，形成一系列高约一个碎石直径的短台阶，并由范围小于一个渠道宽度的水流区域分隔；急滩是非常稳定、陡峭的序列，由不规则巨砾排列而成，水流状态比较急；滑水是由基岩或鹅卵石、砾石组成的均匀且相对无特征的岩层［3，5］。

3．2．2 地貌单元空间分布 依据现场调研结果，利用Arcgis形成地貌单元空间分布图，见图8。统计各河段地貌单元面积、类型，见表3与图9。

图8 各河段地貌单元分布

深潭多出现在河流弯曲段，对岸相邻地貌单元往往是边滩，上下游多是跌水与浅滩单元，如R2段可以明显看出跌水—深潭、深潭—浅滩序列。急滩、滑水多在河流顺直段，上下游多是浅滩单元，对岸紧邻单元往往是边滩。由图8可以看出，河流弯曲段地貌单元远远丰富于河流顺直段。

由表3与图9可知，自上游至下游，R1河段边滩和心滩、江心洲最多，可占56．97%，缓流区与深潭单元次之；R2河段地貌单元分布相对均匀，以串联跌水与缓流区最多，分别占22．3%与21．88%，浅滩与深潭次之；R3、R4河段地貌单元缺少急滩、滑水单元，以边滩、心滩或江心洲、缓流区、浅滩单元为主；R5—R8河段地貌单元分布较为均匀，以深潭、急滩单元为主，缓流区与浅滩单元相对较少，分别为6．88%、3．90%、4．41%、2．60%。针对地貌单元类型，利用SPSS软件，对8河段地貌单元空间分布数据进行方差分析，得到P值为0．014（＜0．05），不同河段地貌单元类型之间存在显著差异。对8个河段地貌单元进行相关性分析，见表4。深潭与滑水、急滩与滑水、浅滩与缓流区存在正相关性，相关系数分别为0．889、0．823、0．779；急滩与缓流区、急滩与滑水存在负相关性，相关系数分别为-0．870、-0．882。这表明，不同地貌单元类型之间存在相关性，但地貌单元类型的具体关系需结合研究区域的实际情况进行量化分析。

表3 各河段地貌单元面积占比 （单位：%）

3．3 地貌异质性分析 基于地貌单元类型的地貌异质性分析方法，梳理实地调研与模型模拟数据，分别计算各河段地貌单元丰富度、地貌单元密度、地貌单元多样性，对其异质性进行量化。具体结果见图10，3个指标所对应范围均为各河段总长。

图10 地貌单元异质性分析（1～8：R1—R8）

从地貌单元丰富度指数来看，R3、R4段地貌异质性明显低于其他河段，分别为0．75、0．88，缺少急滩、滑水地貌单元，其余河段丰富度指数均为1；从地貌单元密度指数、地貌单元多样性指数来看，R1、R3、R4段明显低于其他5段，此3段均为挡水坝下游河段，流量较小，说明小水电的梯级开发对河流地貌异质性造成一定的影响。第R6段、R7段分别对应地貌单元密度指数与多样性指数最大，为32．76与2．66，大于流量最大的R8段，说明除水文影响地貌异质性以外，还有其他因素对其造成影响。地貌单元密度指数与地貌单元多样性指数趋势更加明显，对地貌单元异质性的量化更加准确。但结合3个指数的概念，地貌单元丰富度指数可以直观的反映地貌单元的类型。

4 讨论

4．1 引水式电站梯级开发对河流地貌异质性的影响 引水式电站的梯级开发改变了河流原有基本形态与水文状况，导致河流流量、水温、物质传输等方面发生变化［19］。地貌单元是水流和泥沙共同作用的表现，水文状况发生变化，塑造地貌单元的过程就会发生变化，从而使河流地貌单元异质性发生改变，影响水生生物的多样性［20］。针对梯级开发直接影响河段（R1—R4）和较自然河段（R5—R8）两种类型，分析并比较其地貌异质性，见表5。

表5 梯级开发河段地貌异质性分析

由表5知，R1—R4的地貌单元丰富度（R）、密度（D）、多样性指数（Hd）的平均值均小于R5—R8河段，分别相差9%、32．1%、17．2%。即受梯级开发的影响，地貌异质性明显偏低。R1—R4段的标准差、标准误差的平均值均大于R5—R8段，其中，R2段随着下泄流量的变大，相比R1段，地貌单元丰富度、密度、多样性指数明显变大。R2段为中梁一级电站下游、R5段为中梁三级电站下游（下泄流量增大），其地貌异质性均高于R1、R3、R4段，对于这种可以通过优化调度、增大下泄流量提高地貌异质性的小水电影响河段，在满足 《关于进一步做好小水电分类整改工作的意见》相关要求下，建议对小水电整改保留，在识别地貌单元空间分布基础上，结合生态保护目标进行生态修复。对于地貌异质性低、生态意义重要、对珍稀特有水生生物造成严重影响，且通过修复措施无法达到相关生态要求的小水电影响河段，建议对相关小水电做退出处理。

4．2 河流形态特征与地貌单元的关系 Rosgen［21］依据宽窄率、汊道数、蜿蜒度、宽深比、河流水面坡降和河床底质对河流进行分类，将河流分为7个主要类型、42个亚类，认为不同类型的自然河流都应有其典型的地貌单元类型，如C类型包含深潭-浅滩单元。Montgomery等［22］根据河流比降与河床底质，经过调研与归纳，绘制了理想化的山区冲积河流长剖面图，自上游至下游，随着河床底质组成颗粒的逐渐变小，将河流依次分为底质以基岩为主（粗颗粒）、连续的、高度湍流的串联跌水，水流急缓相间的跌水-深潭，突出单块或几块巨石、水流较为均匀的平面河床，深浅相间的深潭-浅滩序列，底质以沙为主的沙丘-波纹序列等5段。从图8可知，本研究在R2段可以明显看出跌水-深潭序列，但在靠下游的第R5段下游并未发现。位于大宁河干流的第R6、R7、R8段，存在多个深潭-浅滩序列，且多在弯曲段。

蜿蜒河流地貌具有复杂性，在不同的时间尺度上又具备动态性，其往往包含多种地貌单元，有较强的异质性［23］。董哲仁［4］将蜿蜒度大于1．3的河流归为蜿蜒性河流，蜿蜒度为河段两端点之间河流弯曲弧线长度与直线长度的比值。为分析蜿蜒度与地貌单元的关系，结合遥感影像，计算R6、R7、R8段的蜿蜒度，分别为1．22、1．18、1．12。结合图10，在河流水面坡降、底质组成相对一致的情况下，R6、R7、R8段的地貌异质性未随流量的增大而递增。拟合分析蜿蜒度与地貌单元丰富度（R）、密度（D）、多样性指数（Hd）的关系，见图11，蜿蜒度与上述3指数均呈正相关，随着蜿蜒度的增加，地貌单元密度指数明显变高。

图11 不同蜿蜒度下地貌单元异质性分析（R6—R8）

4．3 基于水动力阈值识别地貌单元 不同类型的地貌单元往往表现出离散的水力特征，相比野外直接观测，先确定地貌单元水动力阈值，再结合水动力模型识别的地貌单元分布的方法具有效率快、相对客观等优点［24］。基于Pasternack［25］与Fryirs等［26］的相关研究，河流地貌单元水动力阈值的确定应结合研究区域的实际背景，分区（空间差异性）、分时（时间差异性）、分类确定。（1）分区：不同区域泥沙侵蚀与沉积过程、坡度、底质组成等不同，存在空间差异性。在实地调研中发现R1—R8段同类地貌单元水动力阈值存在明显差异，随着坡度的下降与流量的增大，水动力阈值范围、最值都呈现不同的变化；（2）分时：对于自然、较为自然河流，不同时间如丰水期、平水期、枯水期，水位、水面线存在明显差异，地貌单元水动力阈值需进行动态调整与合理确定，地貌单元空间分布需言明在何种水文条件下。如padmore的分析中，在底质类型组成差异较小的自然河段，随着流量的变化，静水（deadwater）转化为滑水（run），不同流量下滑水与静水的空间分布存在明显差异［27］。（3）分类：不同类型河流，如平原、山区河流等，在地貌单元类型、水动力阈值确定存在明显区别。相对于山区河流，平原河流需兼顾考虑河漫滩地貌单元类型及其水动力阈值的确定。Wyrick等［28］基于水动力模型，以洪峰流量（21 100 ft3/s）作为输入条件确定河漫滩地貌单元类型及其水动力阈值，以全年非洪水期常出现的低流量（880 ft3/s）作为输入条件确定主河道地貌单元类型及其水动力阈值，并提出0．2～0．4倍平滩流量下的水深流速能更好反映地貌变化。

5 结论

本文从区域水文变异分析，河流地貌单元分类、识别、分析等4个方面进行展开，水文变异分析确定了突变点，明确了变异前后各月份水文改变程度；在河流地貌单元分类上，结合前人研究与实地调查数据，确定了以串联跌水、深潭、浅滩等10种单元作为分类依据；在识别上，基于现场观测与高分辨率遥感数据，结合地貌单元过程-形式关系，识别了R1—R8河段的地貌单元空间分布；在分析上，用3个指标对地貌单元异质性进行了定量化，并从地貌单元尺度，探讨了受梯级开发影响的R1、R3、R4段与其余各段地貌异质性差异，尤其是R3、R4段受引水式电站梯级开发影响，缺少急滩、滑水单元，地貌异质性较低。通过明晰河段的地貌单元类型、空间分布及其地貌异质性，可以更好地分析评估小水电对河流生态系统的影响、判别生物栖息地受损情况，结合生态与社会的双目标，从地貌角度为小水电退出、整改、保留的具体实施及其生态修复提供科学依据。