榕江流域多年水沙特性分析

敖亮挺

(广东省水利电力勘测设计研究院有限公司，广东 广州 510635)

榕江俗称南河，又称揭阳江，位于广东省东部，濒临南海，南海水系外河流，其东北和西南方向分布有韩江、练江。榕江是粤东地区的第二大河流，仅次于韩江，也是广东省著名深水河，仅次于珠江，可进出3 000～5 000 t级货轮，直航香港和广州、上海、湛江等地，历史上有“黄金水道”和“状元港”的美誉。榕江上游由南、北两河组成，于砲台镇双溪咀汇合形成榕江干流，经汕头港内的牛田洋注入南海。榕江全长175 km，多年平均径流量61.4亿m3，河床平均坡降0.49‰。其流域面积4 408 km2，山区面积占47.8%，丘陵占16.2%，平原占36%。

南河为榕江主流，发源于汕尾市陆河县凤凰山南麓，发展到揭阳市三洲拦河闸以下进入感潮区，每天出现两次高潮和两次低潮，相邻两次高潮或低潮的潮位不等，涨落潮历时也不等，属不规则半日潮。北河是榕江最大的支流，发源于梅州市丰顺县桐子洋，发展到北河桥闸以下为感潮河段。榕江南河上游河道坡降较陡,洪水传播迅速，而中游河槽平衍淤浅，沙洲碍流，河道弯曲，易受潮汐沙顶托，洪水宣泄不畅，下游则河面宽，吹程大，风暴潮威胁严重。20世纪90年代以前河槽淤浅严重,沙洲不断增多。榕江北河河道弯曲狭窄，坡降平缓，水土流失严重。1957—1959年，由于新岭矿场采砂冲洗泥沙大量入河，河槽淤积较为严重，1960年以后浚河固堤，河道有所疏通。1970年大洪水，新西河水库泄洪道二、三级跌水被冲垮，大量泥沙流入北河，致使赤坎河段淤高1 m左右。榕江河段洪水期受暴雨影响来水较丰，枯水期基本没有径流汇入，由于受到潮流影响，泥沙运动特点较为复杂[1]。榕江流域历年来水土流失比较严重，河道冲淤变化明显，并且水流受到潮汐的影响较大，涨落急时刻流速可达1.5 m/s[2]。在现代经济快速发展的过程中，人类不合理的开发、利用自然资源和工业无限发展对生态环境造成的破坏，导致了河道污染、水质下降、水资源短缺、水土流失等系列系统性问题。地表植被的水源涵养力下降，气候和季节的突变，生态重构对水域变异等问题都迫切地等待人们去不断研究和探讨。

1 研究资料与方法

本文研究区域见图1，南河至东桥园水文站，北河至赤坎水文站，两河于双溪咀汇合，其中南河为榕江主流，为主要研究区域。径流量分析资料取自榕江南河东桥园站1953—2017年的流量数据，榕江北河赤坎站1967—2017年的流量数据；泥沙资料取自榕江南河东桥园站1955—2017年的输沙量与含沙量数据，北河流量资料缺少而与南河流量资料不同步，因此仅可作为辅助分析。因北河与南河水土流失情况近似，故将历年泥沙演变与南河作同等程度对待。

图1 研究区域地理位置

趋势性、突变性和周期性是水文时间序列的重要特征[3-4]。本文采用5 a滑动平均法、Mann-Kendall检验等方法对水沙序列的趋势性与突变性进行分析，利用Morlet小波分析进行周期性分析。下文对这些方法的主要特点和适用性作简要的介绍。

1.1 滑动平均法

滑动平均法是以递推的形式，不断滑动取一定量的相邻数据进行算术平均[5]，以削弱序列中的短周期波动，能有效抑制频繁波动起伏的随机误差，最后得到的滑动平均值随时间的变化可反映序列的变化趋势[6]。滑动平均法的主要特点是简单便捷，但存在一定的主观性和随意性。

1.2 Mann-Kendall检验法

Mann-Kendall检验法是一种非参数检验方法，该方法适用范围广、避免人为干扰、高度定量化，是水文要素评估常会用到的方法。该方法可以检测分析时间序列的变化趋势，并且可以分析出序列发生突变的位置[7]。M-K方法采用原时间序列得到秩序列，由秩序列计算得到随时间变化的统计量过程线UF与UB，其交点通常为突变点。当UF与UB过程线超过临界线时，表明变化趋势显著，并根据UF与UB的正负情况判断序列的升降[6]；当其交点位于临界线外时，表明时间序列显著变异。M-K检验法的显著水平一般取为α=0.05，此时临界值为±1.96。

1.3 有序聚类分析法

有序聚类分析法是通过递推的方式设分割点，以此点将序列分为两部分，求得突变点前后序列的离差平方和，再求出总的离差平方和S，S最小时的分割点为最优分割点，即为序列的突变点[8]。有序聚类分析法简易直观，计算结果较精确。有序聚类分析法识别突变点简单有效，但在水文序列的突变点识别过程中可能产失漏和偏差而导致计算结果不准确。

1.4 小波分析法

小波分析方法源于对信号分析的伸缩与平移[9]，在水文中主要用于具有多时间尺度和振荡特征的时间序列，识别其时间-频率尺度，能够清楚的显示出隐含在时间序列中随时间变化的尺度周期，也能反映时间序列的变化趋势，并对未来的演变趋势作出定性估计[10]。另外，由小波系数得到的小波方差图也可反映出该序列演变过程的主周期。小波分析法通过小波变换系数，能快速获得任意长度的水沙序列的变化周期，且计算结果精确全面，适用性好。

2 研究结果

2.1 榕江流域水沙分布特性

根据东桥园水文站流量资料(1953—2017年)统计，榕江南河多年平均径流量27.62亿m3，多年平均流量87.44 m3/s；1961年对应年径流量和年均流量最大，分别为48.60亿m3和154.12 m3/s；2004年对应年径流量和年均流量最小，分别为12.04亿m3和38.17 m3/s。

根据赤坎水文站流量资料(1968—2017年)统计，榕江北河多年平均径流量为8.17亿m3，多年平均流量26.10 m3/s；1983年对应年径流量和年均流量最大，分别为12.99亿m3和41 m3/s；2004年对应年径流量和年均流量最小，分别为4.83亿m3和15.40 m3/s。

根据东桥园水文站泥沙资料(1955—2017年)统计，榕江南河多年平均含沙量为0.157 kg/m3，最大年平均含沙量为0.321 kg/m3(1959年)，最小年平均含沙量为0.034 kg/m3(2014、2017年)。多年平均输沙率为14.543 kg/s，最大年均输沙率为37.8 kg/s(1973年)，最小年均输沙率为1.288 kg/s(2004年)。因北河与南河水土流失情况近似，故可作同等程度对待。

2.2 径流泥沙的趋势性分析

趋势性分析主要是用于分析时间序列总的走势，结论分别仅有增加、减少或不变3种[11]。趋势分析对预测未来流域演变的趋势有着重要意义。榕江流域多年径流量过程线见图2a、2b，采用5 a滑动平均辅助分析，由图可见，近60年来的南、北河径流量整体上均呈现为年际波动的交替曲线，5 a滑动平均值在线性趋势附近上下波动，表明径流量无明显增加或减少的趋势。

a)南河径流量

多年来沙统计过程线见图2c、2d，由图可见，近60年来的输沙量及含沙量均呈年际波动，且随时间推移均呈整体下降趋势；2005—2013年期间，两者的年际波动现象较为剧烈，2013年后较为稳定。对来沙资料进行线性拟合，得到输沙量以0.914万t/a的趋势减小，含沙量以0.003 16 kg/(m3·a)的趋势减小。

采用Mann-Kendall趋势检验法进一步对水沙资料进行趋势性分析，计算结果见表1，由表可见，东桥园站1953—2017年径流量的统计量为0.47，赤坎站1967—2017年径流量的统计量为0.655，两站统计量|Z|均小于显著性水平α=0.05时的临界值1.96，即表明南河与北河的多年径流量均有增加的趋势，但趋势不显著；输沙量与含沙量的统计量|Z|均大于显著性水平，表明输沙量与含沙量均有显著减小的趋势。

表1 榕江流域径流量、输沙量与含沙量逐年变化趋势M-K检验结果

2.3 径流泥沙的突变性分析

2.3.1泥沙序列突变分析

水文要素的突变性是指：突变位置的前后序列的变化趋势不同[12]。即时间序列从一个均值状态或者趋势突然跳跃到另一个均值状态或趋势[13]。据此对造成突变的自然因素与人为因素进行分析，有利于对将来流域的管理治理。

由图2及表1可知，南河与北河多年径流量无明显趋势性变化，时间序列不存在显著突变，而泥沙序列存在明显的减小趋势，因此，下文仅对输沙量与含沙量进行显著性突变分析。在分析榕江流域近60年来水沙趋势的基础上，通过Mann-Kendall突变检验法确定榕江流域历史输沙量与含沙量的突变时间。

图3a为输沙量的M-K突变检验统计量的曲线。图中可看出在0.05显著水平内，输沙量的统计曲线UF与UB曲线在1961、1975—1980年出现交点，而在显著区间外无交点，表明输沙量的突变点可能在1961年或1975—1980年，但均不显著。含沙量的统计曲线见图3b，可看出在1960、1968—1971、1980—1987、2013—2015年UF与UB曲线均出现交点，而1980—1987年、2013—2015年的交点在显著区间外，为显著性突变可能发生的时间段。

a)输沙量

以上仅仅确定了泥沙序列突变可能发生的时间段，为进一步确定输沙量与含沙量的突变点，采用有序聚类分析法对输沙量与含沙量序列进行的突变检验。检验结果见图4，输沙量与含沙量离差平方和的最小值均在1987年，表明榕江流域的输沙量与含沙量可能在1987年发生了突变。

a)输沙量

为进一步直观展示输沙量与含沙量的突变情况，采用累计距平法对泥沙序列的变化规律进行分析。图5为输沙量与含沙量的累计距平曲线，结合图中信息可知：输沙量与含沙量的累计距平曲线趋势大致相同，均呈先升后降，其由升转为降的年份出现在1987年，表明在1987年前，榕江流域整体处于多沙期，1987—2017年间整个阶段处于少沙期。

图5 输沙量与含沙量累计距平曲线

因此，在结合多种突变分析方法的基础上，得到榕江流域的输沙量在1987年变化不大，而含沙量在1987年变化显著。

2.3.2沙量突变成因分析

a)持续枯水年的作用。根据东桥园站径流资料作差积曲线分析，见图6，东桥园站与赤坎站呈现丰枯水年交替变化，而20世纪80年代末开始，接近20年几乎为连续枯水时间段，河段的输沙量与含沙量因此受到一定的影响。

图6 东桥园站年径流差积曲线

b)建筑工程的大量采砂作用。20世纪90年代前，榕江南河河槽淤浅严重,沙洲不断增多；榕江北河河道弯曲狭窄，坡降平缓，水土流失严重。20世纪80年代后期，榕江流域的整治疏浚开始受到重视，20世纪90年代后，由于建筑工程的大量采砂，河槽开始逐年下切，沙洲面积减小，上游来沙落淤至深槽内，使得流域的输沙量与含沙量均减小。

c)航道疏浚的作用。榕江流域属于平原地带，防洪要求较高，航道整治通常采用挖槽进行[14]。20世纪90年代初，为满足运移量不断增长的需求，榕江流域进行了长距离的航道疏浚，榕江河道水流的挟沙能力随着河道水深的增大和流速的减小而减小，泥沙受水流挟带而运动，水流的挟沙能力与流速成正比。

2.4 径流泥沙周期性分析

水文序列的多年变化周期十分复杂，提取出来的一般为近似周期。周期性跟随研究尺度的不同而发生相应的变化，这种变化一般表现为大尺度的变化周期嵌套着小时间尺度的变化周期[15-16]。得到周期性规律能够更好地预测序列未来的变化规律与未来年份的水文特征，而且可以为流域的防洪与灌溉等工程的设计提供理论帮助[17]。

本文采用Morlet复小波变换，以榕江流域近60 a的水沙序列为基本资料，来探讨榕江流域的水沙周期性变化特征。图7为各序列距平值小波变换的实部等值线图，图中为实线且颜色较浅的部分为正，即为丰水期(多沙期)，虚线且颜色较深的部分为负，即为枯水期(少沙期)。等值线为0(实线与虚线相互过渡的区域)的地方表示丰枯交替(多沙-少沙交替)的转折点。

图7a可看出：南河径流量在1953—2017年的演变过程中可分为3个时间尺度周期，分别为5 a以下、5～8 a与10～15 a。其中10～15 a在全局振荡最强烈，并从1953—2017年经历了枯-丰交替的11个阶段，且2017年刚好出现正值等值线，表明2017年后将进入丰水期；5～8 a则是随机出现。图7b为北河径流量在1967—2017年的小波变换实部等值线，可看出北河径流量在1967—2017年的演变过程中尺度周期为5a以下、5～8 a与10～15 a，其中10～15 a从1967—2017年经历了枯-丰交替的9个阶段，同样2017年后迎来丰水期；80年代末，10～15 a的尺度周期振荡减弱，5～8 a开始活跃。由图7a、7b可得出：榕江流域径流的演变过程中，10～15 a的尺度周期在整个时间段内的发展最稳定，为径流的主周期。

图7c、7d为泥沙序列小波变换的实部等值线，可见图中左侧等值线密度均大于右侧，表明随时间推移榕江流域的输沙量与含沙量均减小，与前文整体趋势性分析的结论一致，而进入21世纪后输沙量与含沙量的等值线较密，表明输沙量与含沙量有所回升，这与前文泥沙序列逐年变化过曲线的状态相符。图中可看出输沙量与含沙量的振荡周期主要为5 a以下、5～8 a与10～15 a，主振荡周期为10～15 a，从1956—2017年经历了“少沙-多沙”交替的11个阶段，与径流演变相对应，而2017年等值线为负值，表明2017年后的一段时间仍为少沙期。由图7c、7d可得出：输沙量与含沙量的演变过程具有与径流演变相似的特点，即10～15 a的尺度周期在整个时间段内的发展最稳定，为演变过程的主周期。但进入21世纪后10～15 a的振荡程度减弱，5～8 a的尺度周期逐渐活跃。

a)南河径流量

为更加直观表现水沙时间序列演变过程的主周期，以小波方差对小波变换进行分析，小波方差的大小反应各时段时间尺度周期振荡的强弱，即方差越大，该时间对应的尺度周期振荡越强[18]。小波方差图的峰值的对应周期尺度即为序列的主周期。从图8可看出，所有序列的方差峰值均出现在时间尺度周期为10～15 a内，约为12 a，表明该尺度周期的振荡最强，为序列的第一主周期，与上述小波实部等值线分析结论一致，图中其他极值点也表明各序列还存在小周期。

a)南河径流量

结合小波系数实部等值线与小波方差分析，取水沙序列的第一主周期12 a附近的小波系数进行趋势分析，见图9，纵坐标为第一主周期的小波系数，其为正值时表示径流量、沙量在此时段偏多，对应为丰水期(多沙期)，为负值时表示径流量、沙量在此时段偏少，对应为枯水期(少沙期)。图中可直观看出南河、北河径流量分别经历了11次与9次丰枯交替，且2017年后小波系数出现正值，表明2017年后将进入丰水期，这与上述径流量小波系数等值线分析结论相互印证；输沙量与含沙量均经历了11次交替，但在2017年的小波系数为负值，表明2017年后一段时间仍为少沙期，与上述泥沙序列的小波系数等值线分析结论相互印证。

a)南河径流量

d)含沙量

2.5 降雨量和径流量分析

降雨量是榕江流域地表径流的唯一来源，地表径流量的变化与降雨量基本一致(表2)[19]。榕江上游地势陡峻，山高谷深，溪润密布，河流落差大，多急流、瀑布，降水强度大，洪水汇流快；揭西县河婆镇以下，河谷逐渐开阔，比降较为平缓，地势平坦；揭阳市三洲拦河闸以下为潮感区，河道更为平缓，两岸农田高程多在3 m以下，不同程度地受台风暴雨、风暴潮影响，历史上洪(潮)涝灾较为严重。

表2 榕江流域代表站多年平均降水量年内分配统计

2.6 水文潮汐分析

流入海洋的河流称为外流河，流入内陆湖泊或消失于沙漠之中的这类瞎尾河称为内流河。榕江由汕头市牛田洋入海，是南海热带气旋常影响和登陆的地区之一，每年受其影响，少则1～2次，多则5～6次，热带气旋(包括热带低压、热带风暴、强热带风暴、台风)带来狂风暴雨暴潮，这是榕江流域的一个重要气候特征。

榕江多年平均水位3.32 m(珠江基面，下同)，历史最高洪水位9.92 m(1970年9月14日)，最低水位2.29 m(1955年3月22日)。据东桥园水文站1954—1979年统计，榕江年平均流量87.4 m3/s，历史上最大洪峰流量4 830 m3/s(1970年9月14日)，最枯流量为0(1954年3月31日断流)。年平均含沙量0.224 kg/m3，年均输沙量63.9万t。北河与南河水土流失情况相似，但北河虽苦旱数月，也从没断流。

榕江受潮汐影响较大，涨潮时北河回水线(感潮区)至琅山上边的北良。南河涨潮时回水线至三洲，枫江涨潮的回水线至潮州市浮岗。

3 结论分析

榕江上游由南、北两河组成，于砲台镇双溪咀汇合形成榕江干流，经汕头港内的牛田洋注入南海。对榕江流域水沙过程的进行趋势性、突变性以及周期性的分析研究，以期对该流域的水资源管理与河道治理方面提供理论参考。研究主要得到以下几点认识和结论。

a)采用Mann-Kendall趋势检验法对榕江流域的水沙资料进行趋势性分析，得出东桥园站与赤坎站两站统计量均为正数，但绝对值小于显著性水平α=0.05，表明南河与北河的多年径流量均有增加的趋势，但并不显著，径流量整体上围绕一定值呈年际波动；输沙量与含沙量的统计量均为负数，但绝对值大于显著性水平，表明输沙量与含沙量均有显著减小的趋势。

b)通过组合Mann-Kendall突变检验法、有序聚类分析法与累计距平法对有显著变化的输沙量与含沙量序列进行突变性分析，经检验分析，输沙量与含沙量的突变点均为1987年，进一步确定榕江流域的输沙量在1987年发生不显著突变，而含沙量在1987年发生显著突变；以1987年为转折点，榕江来沙过程整体上经历了多沙、少沙2个过程。

c)采用Morlet复小波变换分析榕江流域的水沙周期性变化特征，以小波系数方差图与其相互印证，得到径流量、输沙量与含沙量的第一主周期约为12 a。以12 a为尺度周期，径流量从1953—2017年经历了枯-丰交替的11个阶段，并在2017年后将进入丰水期；来沙过程从1955—2017年经历了少沙-多沙交替的11个阶段，在2017年后一段时间仍为少沙期。

综上几种分析方法和数据可以看出，全球变暖，气候变化，频繁的人类活动等各类生产经济的建设，会直接对江河流域的趋势性、突变性以及周期性带来不可恢复的变化。在一定的空间尺度上，这种变化具有一定的重复性，也就是周期性。Mann-Kendall趋势检验法，Morlet 小波分析法，在某周期内能较为准确分析计算出流域的趋势和流域的周期突变，并开展流域水沙趋势的预测。