一种库区区间来沙量的计算方法与流程

1.本发明涉及一种库区区间来沙量的计算方法，用于开展库区区间来沙量计算研究，主要应用于水力学及河流动力学领域。


背景技术：

2.库区淤积量是一个关系到水库调度和水库使用寿命的重大问题，以往对水库淤积进行计算时，没有充分考虑库区区间来沙。然而，近年来由于干流大型水利工程的建设以及水土保持工作推进，世界上多数河流的输沙量呈现减小趋势，库区区间来沙的权重相对增大，对水库排沙比影响研究不可忽略。在实际水库中，如河道型水库的代表三峡水库，库区回水长度达到500～660km，区间支流众多。根据估算，2003～2016年三峡库区区间来沙约占总入库沙量的10.3％；2013年后由于长江上游向家坝、溪洛渡水电站相继蓄水运用，库区区间来沙占比在2013～2016年间增大至26.9％。
3.然而，库区区间支流目前仍存在水文资料缺测，或其控制站控制面积无法覆盖完整支流流域的情况，且流域干支流往往还建有梯级水库，梯级水库群拦沙又使得库区区间来沙量的计算无法准确并快速的确定。如三峡水库，库区区间重要支流的水文控制站多在支流中上游，无法通过各支流测站的输沙量相加来估算区间输沙量，又由于三峡水库蓄水后河道淤积，也无法用干流出口站减去进口站的方法求得区间输沙量。而且三峡区间受人类活动影响频繁且强度很大，存在支流梯级水库建设拦沙，且下垫面条件变化剧烈，不同时期其径流量～输沙量关系存在较大差异，也难以直接通过这一关系来求区间的输沙量。要较为准确地估算区间来沙量很困难。故至今三峡水库区间鲜见有定量的月、年来沙量估算成果。
4.因此，如何快速、准确计算完整库区区间支流来沙量，是一个需要解决的问题。


技术实现要素：

5.为了克服现有技术的问题，本发明提出了一种库区区间来沙量的计算方法。所述的方法基于已知支流水文测站的输沙量以及雨量站的降雨量数据，计算水文站控制流域的还原输沙模数以及未控区域的输沙模数，从而合理定量计算完整库区区间支流来沙量。在库区区间支流缺乏可覆盖完整支流汇流区域实测数据情况下，实现库区区间来沙量的合理计算。
6.本发明的目的是这样实现的：
7.一种库区区间来沙量的计算方法，包括以下步骤：
8.步骤1，库区区间分区及数据收集：
9.收集待计算水库区间内各支流控制站的实测水沙数据，以及各支流控制站附近雨量站的同期降雨数据，筛选获得资料基础好的若干控制站作为代表，结合待计算水库区间内流域、地质地貌条件情况，将待计算水库区间分为若干片区。
10.步骤2，计算各片区代表站控制区域还原输沙模数：
11.假定各水文控制站控制流域内输沙模数相等，各支流上若建有若干水库，这些水库的拦沙作用不可忽略。因此，在根据控制站输沙数据计算输沙模数时，还需考虑控制站上游水库的拦沙作用，还原该控制站的输沙量，得到更为可靠的输沙模数计算结果。
12.故按照以下公式，分别计算步骤1筛选获得的各代表站控制区域的还原输沙模数：
[0013][0014]
式中，sr是代表控制区域的还原输沙模数；gs为水文控制站输沙量；a为水文站控制流域面积；a’为水库控制流域面积；η为水库拦沙率。
[0015]
步骤3，建立降雨与输沙模数的经验关系：
[0016]
降雨侵蚀力对流域输沙模数的作用通常不是线性效应，故采用二阶多项式回归模型，拟合步骤2计算得到的各代表站控制区域还原输沙模数与其对应雨量站的降雨侵蚀力关系，所述二阶多项式回归模型为：
[0017]
s＝a0+a1·er
+a2·er2
+ε
[0018]
式中，s为输沙模数；er为降雨侵蚀力；a0、a1、a2为模型参数；ε为模型误差。
[0019]
步骤4，计算水文站未控区域输沙模数：
[0020]
按照以下公式，分别计算待计算水库区间内各未控区域的输沙模数：
[0021]
si＝s
·k[0022]
式中，si为待求未控区域输沙模数；s为步骤3计算获得的代表站控制区域输沙模数；k为修正系数，当未控区域与其所在区域的代表站控制区域下垫面条件有明显差异时，采用k值进行修正。
[0023]
上述二阶多项式回归模型是基于实测数据回归拟合得到，反映了在相似地貌条件、土地利用及植被条件下输沙模数与降雨侵蚀力的关系，因此可用于未控区域的输沙模数计算。
[0024]
步骤5，计算支流来沙量：
[0025]
支流来沙量计算以支流流域各区域的输沙模数为基础，并考虑沿程水库的拦沙情况，按照以下公式，对待计算水库区间内各片区的来沙量进行计算：
[0026][0027]
式中，ηj为水文控制站下游至入库河段上的第j个水库拦沙率；aj为第j个水库的流域控制面积减去第j+1个水库的流域控制面积；sj为对应流域面积上的输沙模数。
[0028]
进一步的，步骤2中，所述水库拦沙率η的确定方法包括：
[0029]
(1)根据水库水文监测资料计算获得；
[0030]
(2)对于水库拦沙量及拦沙率缺乏实测资料的小型水库，参照同流域有实测资料的水库淤积情况，确定其水库拦沙率。
[0031]
进一步的，流域产沙与降雨侵蚀力密切相关，而降雨侵蚀力同时受降雨量、降雨强度等因素影响，为表征降雨侵蚀力，步骤3中，以降雨天数、平均雨强为变量的所述降雨侵蚀力er的经验表达式为：
[0032]er
＝d
·
ra[0033]
式中，d为统计时段内出现降雨天数的概率；r为雨强，即统计时段内总降雨量与降雨天数的比值；a为待定参数，通过率定确定。
[0034]
本发明的优点和有益效果是：
[0035]
本发明利用已知支流水文测站的输沙量以及雨量站的降雨量数据，计算水文站控制流域的还原输沙模数以及未控区域的输沙模数，从而合理定量计算完整库区区间支流来沙量。在库区区间支流缺乏可覆盖完整支流汇流区域实测数据情况下，实现了库区区间来沙量的合理计算。
附图说明
[0036]
下面结合附图和实施例对本发明作进一步说明。
[0037]
图1是本发明实施例所述三峡水库区间输沙量计算片区划分图；
[0038]
图2是本发明实施例所述输沙模数与降雨侵蚀力关系结果图。
具体实施方式
[0039]
实施例：
[0040]
本实施例以三峡库区为具体算例，提供了其库区区间来沙量的计算方法，对本发明的技术方案进行详细的说明，同时通过分析，说明本发明的结果合理性。所述方法包括以下步骤：
[0041]
步骤1，库区区间分区及数据收集：
[0042]
通过收集三峡库区间内各支流控制站的实测水沙数据，确定在三峡库区区间支流中，龙河石柱站、小江温泉站、磨刀溪长滩站、大宁河巫溪站和香溪河兴山站有较好的资料基础，以这5个控制站作为代表。再结合三峡库区内流域、地质地貌条件情况，如图1、表1所示，本实施例以龙河、小江、香溪河、大宁河以及磨刀溪为代表将三峡库区分为5大片。
[0043]
表1三峡区间分区情况表
[0044][0045]
根据上述5个代表站的水文数据统计，2003～2016年上述5个流域年均来流量占入库总流量的约4％，而相关研究表明三峡库区区间来流约占入库总流量10％左右，因此，上述5个流域来流约占三峡区间入库流量的40％，可一定程度代表三峡库区区间的来水来沙情况。同时，收集到了4个雨量站的逐日降雨信息，分别为达川站、奉节站、巴东站和万州站。本实施例收集到的水文站流量及输沙数据均为2010～2016年实测数据，其中，水文站数据为月均数据，雨量站数据为日均数据。
[0046]
步骤2，计算各片区代表站控制区域还原输沙模数：
[0047]
按照以下公式(1)，分别计算步骤1筛选获得的各代表站控制区域的还原输沙模
数：
[0048][0049]
式中，sr是代表控制区域的还原输沙模数(t/(km2.a)；gs为水文控制站输沙量(t/a)；a为水文站控制流域面积(km2)；a’为水库控制流域面积(km2)；η为水库拦沙率。
[0050]
各片区的具体计算结果为：
[0051]
龙河片区：龙河片区的代表性水文站石柱站控制面积为898km2，其以上区域属于龙河流域上游，占流域总面积的32.3％，地貌条件，土地利用及植被条件相差不大，可以粗略地认为各地土壤侵蚀模数相同。龙河上游小型水库拦沙能力有限，水库拦沙可以只考虑滕子沟水库(2005年开始运行)。由于藤子沟水库建成后，水库拦沙量及拦沙率尚无实测资料，可以参照长江流域其它水库的淤积情况。金沙江毛家村水库运用方式与藤子沟水库相似，毛家村水库拦沙率为80％(金沙江流域水沙变化原因调查与变化趋势分析，2006)。参照毛家村水库的拦沙情况，假设藤子沟水库滕拦沙率为80％，假设石柱上游各区域侵蚀模数相同，则可将石柱站的输沙量按流域面积比例分解为滕子沟上游流域(600km2)和其它区域(298km2)两个部分，则可由式(1)求得龙河片区代表站控制区域2010-2016年间逐月还原输沙模数。
[0052]
磨刀溪片区：磨刀溪是三峡库区长江一级支流，流域面积3049km2。磨刀溪流域梯级水库已经基本建成，其中鱼背山水库具有较大的拦沙库容，在计算拦沙作用时仅考虑鱼背山水库(1998年8月开始拦沙)，拦沙率按0.8计，鱼背山坝址控制流域面积为1389km2。鱼背山水库下游有长滩水文站，控制面积为2034km2。假设长滩水文站上游各区域侵蚀模数相同，则可将长滩站的输沙量按流域面积比例分解为鱼背山上游流域(1389km2)和其它区域(645km2)两个部分，则可由式(1)求得磨刀溪片区代表站控制区域2010-2016年间逐月还原输沙模数。
[0053]
小江片区：小江是长江的一级支流，流域面积5205km2。温泉水文站位于小江上游，控制面积为1158km2，控制面积占小江流域比例为22.24％。温泉水文站上游已建成的大中型水库主要有鲤鱼塘水库，鲤鱼塘坝址控制面积为235.8km2，2008年建成蓄水，水库拦沙率均按0.8计。假设温泉水文站上游各区域侵蚀模数相同，则可将温泉站的输沙量按流域面积比例分解为鲤鱼塘上游流域(235.8km2)和其它区域(922.2km2)两个部分，则可由式(1)求得小江片区代表站控制区域2010-2016年间逐月还原输沙模数。
[0054]
大宁河片区：大宁河又名巫溪水，是长江北岸一级支流，流域面积4407km2。巫溪水文站位于大宁河干流，控制面积为2001km2，占流域总面积的45.43％。位于巫溪水文站上游有两座电站具有较大库区，分别是位于支流后溪河的孔梁水库，控制面积58.5km2，2005年开始拦沙，以及位于支流西溪河的中梁水库，控制面积525km2，2011年开始拦沙。孔梁水库和中梁水库流域地貌条件，土地利用及植被条件相差不大，可以粗略地认为各地土壤侵蚀模数相同。参照毛家村水库的拦沙情况，假设各水库拦沙率为80％，则可将巫溪站输沙量按流域面积比例分解为水库上游流域(2005～2010年为58.5km2，2011～2015年为58.5+525km2)和其它区域两个部分，则可由式(1)求得大宁河片区代表站控制区域2010-2016年间逐月还原输沙模数。
[0055]
香溪河片区：香溪河是长江左岸西陵峡北面一级支流，由香溪河及西河、东河、高
岚河三条支流汇合而成，干流全长101km，流域面积3214km2。兴山站为香溪控制性水文站，控制流域面积1900km2，占流域总面积的59.2％。兴山站上游需要考虑古洞口梯级(古洞口一级电站和古洞口二级电站)的拦沙作用，古洞口梯级2005开始投入运用，总控制流域面积为1161km2。古洞口水库建成后，水库拦沙量及拦沙率尚无实测资料，参照毛家村水库的拦沙情况，假设各水库拦沙率为80％。可将兴山站输沙量按流域面积比例分解为古洞口水库上游、和兴山上游其它区域两个部分，则可由式(1)求得香溪河片区代表站控制区域2010-2016年间逐月还原输沙模数。
[0056]
步骤3，建立降雨与输沙模数的经验关系：
[0057]
采用二阶多项式回归模型，拟合步骤2计算得到的各代表站控制区域还原输沙模数与其对应雨量站的降雨侵蚀力关系，所述二阶多项式回归模型为：
[0058]
s＝a0+a1·er
+a2·er2
+ε
ꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0059]
式中，s为输沙模数；er为降雨侵蚀力；a0、a1、a2为模型参数；ε为模型误差。
[0060]
以降雨天数、平均雨强为变量的所述降雨侵蚀力er的经验表达式为：
[0061]er
＝d
·
raꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0062]
根据4个雨量站数据，由式(3)计算2010-2016年各站逐月降雨侵蚀力(er)。然后，采用式(2)二阶多项式回归降雨侵蚀力与各雨量站对应片区的代表水文站还原输沙模数关系，见图2。
[0063]
由图2可见，输沙模数与降雨侵蚀力之间相关关系较好，此时，a取1.47。
[0064]
步骤4，计算水文站未控区域输沙模数：
[0065]
以图2回归模型计算式，以及下述公式(4)，分别计算待计算水库区间内各未控区域的输沙模数：
[0066]
si＝s
·kꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0067]
式中，si为待求未控区域输沙模数；s为步骤3计算获得的代表站控制区域输沙模数；k为修正系数，当未控区域与其所在区域的代表站控制区域下垫面条件有明显差异时，采用k值进行修正。k值需要根据各片区具体情况来确定，对于三峡库区区间流域，土地利用及植被条件相差不多，主要需考虑地形坡度对产沙输沙的影响。
[0068]
各片区的具体计算结果为：
[0069]
龙河片区：龙河流域与磨刀溪流域相邻，地质地貌条件相近，因磨刀溪流域的长滩和龙角有对比输沙量观测资料，上游部分输沙模数是下游部分的一半，因此参考磨刀溪上下游输沙模数关系，石柱下游输沙模数取石柱上游还原输沙模数的2倍，即k＝2。龙河片区其它区域的输沙模数的地区分布也参考磨刀溪流域的分布，上游地区(与石柱站控制流域比例相同)输沙模数与石柱上游还原输沙模数相同(k＝1)，下游地区以2倍输沙模数计(k＝2)。
[0070]
磨刀溪片区：磨刀溪流域的长滩和龙角有对比输沙量观测资料，上游部分输沙模数是下游部分的一半，长滩站下游输沙模数取长滩上游的2倍(k＝2)。磨刀溪片区其它流域的输沙模数的地区分布也参考磨刀溪的分布，上游地区(与长滩站控制流域比例相同)输沙模数与长滩站上游还原输沙模数相同(k＝1)，下游地区以2倍输沙模数计(k＝2)。
[0071]
小江片区：小江片区水库数量众多，水库控制流域面积比例大，地形高差相对较小，坡度也较小，侵蚀强度较弱，因此，该片区的上下游输沙模数可认为一致，均取k＝1。
[0072]
大宁河片区和香溪河片区：大宁河片区和香溪河片区的输沙模数分布同样可参考磨刀溪流域，上游地区(与巫溪站控制流域比例相同)输沙模数与巫溪上游还原输沙模数相同(k＝1)，下游地区以2倍输沙模数计(k＝2)。
[0073]
步骤5，计算支流来沙量：
[0074]
按照以下公式(5)，对待计算水库区间内各片区的来沙量进行计算：
[0075][0076]
式中，ηj为水文控制站下游至入库河段上的第j个水库拦沙率；aj为第j个水库的流域控制面积减去第j+1个水库的流域控制面积；sj为对应流域面积上的输沙模数。
[0077]
各片区的来沙进行计算，主要需确定拦沙水库的影响，各片区的具体计算结果为：
[0078]
龙河片区：龙河流域在石柱水文站下游还有牛栏口、石板水和鱼剑口三座较大的水库，其中，位于龙河干流中段的石板水水库起主要拦沙作用，因此龙河中下游的输沙仅考虑石板水水库的拦沙作用。石板水坝址以上流域面积1640km2，无水文监测资料，按照长江上游水库的拦沙情况，拦沙率可按80％计。龙河片区其它区域的水库拦沙只考虑大溪河中游的鱼跳水库，水库拦沙率仍按0.8计。大溪河流域面积2065km2，鱼跳水库坝址以上控制集水面积1380km2。
[0079]
磨刀溪片区：磨刀溪流域在长滩水文站下游没有需要考虑拦沙的水库。该片区其他流域的水库拦沙只考虑长滩河的清水湖水库(2013年开始拦沙)，水库拦沙率按0.8计。长滩河流域面积1526km2，清水湖水库坝址以上控制集水面积1077km2。
[0080]
小江片区：小江流域温泉水文站下游需要考虑开县水位调节坝水库的拦沙作用，开县调节坝控制面积为3198.6km2，2012年5月开始蓄水，水库拦沙率均按0.8计。小江片区其它流域水库密布，很多流域都完成了梯级开发，水库控制流域面积比例很大。进行拦沙计算时只考虑狮子滩、大洪河、白石、河流水上游梯级及白水河梯级。狮子滩位于龙溪河，控制流域面积3020km2，1956年10月建成发电。大洪河水库位于大洪河，控制流域面积1317km2，水库于1960年建成投产，此外，大洪河流域合流水及四海上游有五个梯级，四海水文站控制面积1153km2，合流水控制面积930km2，白水河292km2，上游有三个梯级。黄金河有六个梯级水库，白石水库是唯一调节水库，坝址以上控制集雨面积173km2，白石水库扩建工程于2011年1月试通水成功。水库拦沙率均按0.8计。
[0081]
大宁河片区：大宁河流域巫溪水文站下游需考虑刘家沟水库的拦沙作用，刘家沟水库位于支流小溪河，控制面积67.3km2，2013年开始拦沙，拦沙率以0.8计。大宁河片区其他流域主要有梅峰水库、团滩河(水库)、草坪河水库、青莲溪水库及渡口坝电站水库。由于其它水库控制流域面积较小，水库拦沙作用只考虑渡口坝和团滩水库，其中梅溪河渡口坝水库控制流域面积764.9km2，库容9254万m3，2011年蓄水拦沙。汤溪河团滩水库控制流域面积226.1km2，库容1687万m3，2011年蓄水拦沙。水库拦沙率仍按0.8计。
[0082]
香溪河片区：香溪河流域兴山站下游没有需要考虑拦沙的水库。香溪河片区需要考虑拦沙的有清港河流域的升平和观音堂两座中型水库。升坪一级电站为可调节混合式电站，总库容2740万m3，坝址以上承雨面积117km2，2002年建成。观音堂水电站水库总库容4618万m3，控制流域面积509km2，电站2009年底全面开工，2011年4月试运行发电。水库拦沙率仍按0.8计。
[0083]
经过上述步骤后，最后可计算得到各片区2010-2016年逐月入库沙量，表2展示了
2010-2016年逐年入库沙量结果：
[0084]
表2三峡库区区间来沙量(万t)
[0085]
年份香溪河片区小江片区磨刀溪片区大宁河片龙河片区区间来沙合计20108634314724124410612011242836363420509237020123055872541463951687201326411853331346382553201441314873174096013226201524230117212129411302016 4294194046991950平均2597382862684832034
[0086]
计算结果合理性分析：
[0087]
长江科学院基于遥感影像(三峡水库区间坡面侵蚀产沙及调控研究，2013)，采用土壤侵蚀量法估算的2010年三峡水库区间入库沙量为1009万t，与采用本实施例所述方法计算的2010年三峡区间入库沙量1061万t十分接近。
[0088]
长江水利委员会水文局(三峡水库区间来沙量分析研究报告，2016)对三峡水库寸滩至宜昌区间的来沙量进行了估算，其2003～2015年年均来沙量约为2000万t，与本实施例所述方法计算的2010～2016年年均来沙量2034万t接近。
[0089]
综上，说明本发明提出的库区区间来沙量计算方法是合理的。本发明的计算方法成本低，计算便捷，且可根据雨量、水文站数据精度提供时间精度更高的数据结果，如逐月区间入库沙量数据，同时计算结果能够与多家采用不同方法的计算值相吻合，能够准确计算库区区间来沙量大小。
[0090]
最后应说明的是，以上仅用以说明本发明的技术方案而非限制，尽管参照较佳布置方案对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围。