CN106969756A - 河流悬移质泥沙输沙量改正计算方法 - Google Patents

Abstract

本发明提供一种河流悬移质泥沙输沙量改正计算方法，其特征在于，包括以下步骤：步骤1.计算断面平均含沙量，得到两种断面含沙量的平均值D_CS_I和D_CS_II；步骤2.确定边缘分布，计算联合分布函数参数：根据D_CS_I和D_CS_II资料序列，确定边缘分布函数F₁(CS_I)和F₂(CS_II)以及Copula函数C_θ(u,v)的参数θ；步骤3.区间离散：将取值区间均匀的分为n等分，每等分h＝(D_CS_IImax‑D_CS_IImin)/n，并记子区间[D_CS_II(k)，D_CS_II(k+1)]的中点为x_k+1/2＝x_k+h/2；步骤4.建立基于Copula函数的联合分布模型：根据给定的断面平均含沙量x_CS_I和x_CS_II，边缘分布F₁(CS_I)和F₂(CS_II)，计算得到和以及y的反函数和Copula函数c(u,v)；步骤5.基于条件期望改正计算含沙量。本方法具有较强的适用性，能够有效提高含沙量改正计算的精度，并减小水文测验的工作强度。

Description

河流悬移质泥沙输沙量改正计算方法

技术领域

本发明属于水文测验领域，具体涉及河流悬移质泥沙输沙量改正计算方法。

技术背景

河流中的泥沙，可以淤积河道，造成河流的泛滥和游荡。也可以改良土壤，使盐碱沙荒变为良田。对一个流域或一个地区，为了达到兴利除害的目的，就要了解泥沙的特性、来源、数量及其时空变化，为流域的开发和国民经济建设，提供可靠的依据。为此，必须开展泥沙测验工作，系统地搜集泥沙资料。

《河流悬移质泥沙测验规范》(GB50159-92)中规定，对于悬移质输沙率测验，测验内容包括断面输沙率测验和单位水样含沙量(简称单沙)测验。断面输沙率是指单位时间内通过河渠某一断面的悬移质沙量。单位水样含沙量是指断面上有代表性的垂线或测点的含沙量。

断面输沙率的测验是为了准确推求断面平均含沙量，测验时根据泥沙在横向分布变化情况，布设若干条垂线(见图1)。取样方法有：在每条垂线的不同测点上，逐点取样，称积点法；各点按一定容积比例取样，并予混合称定比混合法；各点按其流速比例确定取样容积，并予混合，称流速比混合法；用瓶式或抽气式采样器在垂线上以均匀速度提放，采取整个垂线上的水样，称积深法等。可根据水情、水深和测验设备条件合理选用。

如图1所示，P1、P2、P3、P4、P5为代表测验垂线上不同水深处的测点，不同的测验断面数量也可能不同，较为常见的主要有以下几种类型，即一点法：在水面下0.6相对深处测速取沙；二点法：在水面下0.2和0.8相对水深处，测速取沙；三点法：在水面下0.2、0.6、0.8相对深处测速取沙；五点法：在水面、水面下0.2、0.6、0.8相对水深和河底附近测速取沙；七点法：水面，水面下0.2、0.6、0.8、0.9、河底以上0.5米、河底以上0.1～0.2米附近测速取沙。

断面输沙率测验须与流量测验同时进行，需要进行颗粒分析的测次，同时加测水温。由于断面输沙率测验工作量大，费时较多，不可能把断面输沙率变化的每一个转折点都实地测到，更不能在泥沙变化大时逐时实测。因此，运用实测断面输沙率与测定单位水样含沙量两者相结合的方法，即在测得的断面输沙率资料中，选取1条或2～3条垂线的平均含沙量同断面平均含沙量建立稳定对应关系；这样，只要在此选定的1条或2～3条垂线的位置上测取水样，求得此单位水样含沙量后，通过上述稳定的对应关系，即可求得断面平均含沙量，并与相对应的时段平均流量相乘，即得该时段的平均输沙率，然后乘上所经的历时并累积相加，即得各种时段如日、月、年等的输沙量。

长期以来，在进行河道冲淤量计算时，输沙量法计算的冲淤量与同期体积法河道泥沙冲淤量相比不太匹配。其中无法准确获得临底层悬移质含沙量(简称临底悬沙)是问题的关键。

目前悬移质泥沙测验测量范围多是在距水面0.8倍水深以上，较精密者也就采用五点法。

但是限于仪器构造，距河床0.5米以下浓度较高或粒径较粗的沙重，却被漏测了。例如，图1中的临底悬沙区，通常距河床0.5米以下，该区域内悬移质泥沙的浓度最高，但限于仪器的构造，通常会被漏测，从而影响到最终输沙量的计算成果(图1中仅给出了3条测线，实际工作中根据测验断面的宽度、断面形态以及水文特性的不同，代表性测验垂线的条数会有不同)。

这样就导致实际测算的输沙率与真实值不符，往往相差较大，从而使年输沙量的成果精度受到影响，必须通过实测资料的试验与分析计算，改正实测悬移质输沙率，以便得到比较符合实际的数值。

现行输沙量改正计算方法存在的问题是：

(1)适用性问题：现行的输沙量改正计算方法多以含沙量垂线分布概化曲线公式为核心，但由于河流特性和测站特性的不同，天然河道中断面几何形状不规则，流速、含沙量横向分布不均匀等的影响因素，概化垂线公式的适用性不强。

(2)参数选取问题：流速和含沙量分布公式对改正计算成果的影响较大，流速分布公式多采用对数公式，公式参数包括水深、比降和粗糙度等；含沙量分布公式多采用Rouse公式，需要在对数格纸上点绘直线求得水面流量和分布指数等参数。故参数的选取过程比较繁琐，且往往具有较大的不确定性。

发明内容

本发明是为了解决上述问题而进行的，目的在于提供一种河流悬移质泥沙输沙量改正计算方法，具有较强的适用性，能够提高含沙量改正计算的精度，并减小水文测验的工作强度。

本发明为了实现上述目的，采用了以下方案：

本发明提供一种河流悬移质泥沙输沙量改正计算方法，其特征在于，包括以下步骤：

步骤1.计算断面平均含沙量

将各实测垂线i同一相对水深η的测点流速V_η用面积Ω加权：

将含沙量CS_η用面积Ω与流速V_η二者加权：

再根据两组水深权值，分别计算出两种断面含沙量的平均值D_CS_I和D_CS_II，

式中，η为垂线上测点的相对水深值；Ω_i为以垂线i间距的中点为分界线的部分测验断面面积；Ω为全测验断面面积；V_η-i为第i条垂线，在相对水深η处的测点流速；V_η为各垂线在同一相对水深η处的平均流速；CS_η为各垂线在同一相对水深η处的平均测点含沙量；K_i'和K_i”均为面积比系数；I对应第一组水深权值，II对应第二组水深权值；

步骤2.确定边缘分布，计算联合分布函数参数

根据计算得到的D_CS_I和D_CS_II资料序列，通过拟合优度评价指标，确定其边缘分布函数F₁(CS_I)和F₂(CS_II)，并确定联合分布Copula函数C_θ(u,v)的参数θ；

步骤3.区间离散

根据D_CS_II资料序列的最大值D_CS_IImax和最小值D_CS_IImin，将取值区间[D_CS_IImax，D_CS_IImin]均匀的分为n等分，每等分h＝(D_CS_IImax-D_CS_IImin)/n，并记子区间[D_CS_II(k)，D_CS_II(k+1)]的中点为x_k+1/2＝x_k+h/2；

步骤4.建立基于Copula函数的联合分布模型

根据给定的断面平均含沙量x_CS_I和x_CS_II，边缘分布F₁(CS_I)和F₂(CS_II)，计算得到和以及y的反函数和Copula函数c(u,v)；

步骤5.基于条件期望改正计算含沙量

对于第i个区间，令利用复化辛普森公式，通过数值积分的方法计算得到，给定断面平均含沙量条件x_CS_I下，断面平均含沙量x_CS_II的期望值E(y|x)。

本发明提供的河流悬移质泥沙输沙量改正计算方法，还可以具有以下特征：在步骤2中：离差平方和准则为式中，P_i和P_ei分别为理论频率与经验频率；n为序列个数；Copula函数f(cs_I,cs_II)为f(cs_I,cs_II)＝C_θ(F₁(cs_I),F₂(cs_II))＝C_θ(u,v)；选用Archimedean Copula函数族中的G-H Copula函数作为联合分布函数，选用RMSE准则评价Copula函数的拟合情况；G-H Copula函数的数学表达式为：C_θ(u,v)＝exp{-[(-ln u)^θ+(-ln v)^θ]^1/θ}，θ＞1；参数θ通过Kendall秩相关系数τ求得：θ＝1/(1-τ)。

本发明提供的河流悬移质泥沙输沙量改正计算方法，还可以具有以下特征：在步骤1中所描述的两组水深权值为用三点法、五点法、七点法测量得到的三组深权值中的任意两组。

本发明提供的河流悬移质泥沙输沙量改正计算方法，还可以具有以下特征：第一组水深权值为用五点法测量得到的一组深权值，所述第二组水深权值为用七点法测量得到的一组深权值。

本发明提供的河流悬移质泥沙输沙量改正计算方法，还可以具有以下特征：期望值E(y|x)的计算公式为：

发明的作用与效果

与现有技术相比，本发明的有益效果在于：

(1)本发明从各测验断面实测历史样本特性出发，利用统计学方法优选模型和参数，较传统的含沙量垂线分布经验公式更具有适用性；

(2)本发明的核心计算方法基于联合分布模型及条件期望公式，具有较强的理论基础，从而显著提高含沙量改正计算的精度；

(3)本发明通过统计学及数学公式的推导，可以减少泥沙测验代表垂线上的测点样本数量，从而减小水文测验的工作强度。

综上，本发明提供的河流悬移质泥沙输沙量改正计算方法，可广泛应用于水文测验技术领域内的输沙量改正计算，不仅更为简洁高效，而且能够提高含沙量改正计算的精度，并减小水文测验的工作强度。

附图说明

图1为传统泥沙测验方法测点分布图；

图2为本发明实施例中河流悬移质泥沙输沙量改正计算方法的流程图。

具体实施方式

以下结合附图对本发明涉及的河流悬移质泥沙输沙量改正计算方法的具体实施方案进行详细地说明。

<实施例>

如图2所示，本实施例所提供的河流悬移质泥沙输沙量改正计算方法包括以下步骤：

步骤1.计算断面平均含沙量

将各实测垂线i同一相对水深η的测点流速V_η用面积Ω加权：

将含沙量CS_η用面积Ω与流速V_η二者加权：

以上各式中，η为垂线上测点的相对水深值(由河底为零起称)；Ω_i为以垂线i间距的中点为分界线的部分测验断面面积；Ω为全测验断面面积；V_η-i为第i条垂线，在相对水深η处的测点流速；V_η为各垂线在同一相对水深η处的平均流速；CS_η为各垂线在同一相对水深η处的平均测点含沙量；K_i'和K_i”均为面积比系数，

然后，根据两组水深权值，分别计算出两种断面含沙量的平均值D_CS_I和D_CS_II；这里，两组水深权值为用三点法、五点法、七点法测量得到的三组深权值中的任意两组；本实施例中，第一组水深权值为用五点法测量得到的一组深权值，第二组水深权值为用七点法测量得到的一组深权值：

D_CS_I＝D_CS₅＝(CS₀+3CS_0.2+3CS_0.6+2CS_0.8+CS₁)/10 (3)

D_CS_II＝D_CS₇＝(CS₀+CS_0.2+2CS_0.6+2CS_0.8+2CS_0.9+CS_h-0.5+CS_h-0.2)/10 (4)

步骤2.优选边缘分布，计算联合分布函数参数

根据计算得到的D_CS₅和D_CS₇资料序列，通过拟合优度评价指标，确定其边缘分布函数F₁(CS₅)和F₂(CS₇)，并确定联合分布Copula函数C_θ(u,v)的参数θ，具体计算公式如下：

1)离差平方和准则(RMSE准则)为

式中，P_i和P_ei分别为理论频率与经验频率；n为序列个数；

2)Copula函数

Copula函数是构建联合分布的一种有效方法，假定五点法和七点法断面平均含沙量的边缘分布函数分别为F₁(CS₅)以及F₂(CS₇)，由Sklar定理，存在一个二维的Copula函数f(cs₅,cs₇)，使得：

f(cs₅,cs₇)＝C_θ(F₁(cs₅),F₂(cs₇))＝C_θ(u,v) (6)

由于五点法和七点法断面平均含沙量存在正相关关系，本实施例中选用Archimedean Copula函数族中的Gumbel-Hougaard Copula函数(G-H Copula函数)作为联合分布函数，选用RMSE准则评价Copula函数的拟合情况；G-H Copula函数的数学表达式为：

C_θ(u,v)＝exp{-[(-ln u)^θ+(-ln v)^θ]^1/θ}，θ＞1 (7)

式中，参数θ可通过Kendall秩相关系数τ求得：θ＝1/(1-τ)。

表1不同类型分布函数统计表

步骤3.区间离散

根据D_CS₇资料序列的最大值D_CS_7max和最小值D_CS_7min，将取值区间[D_CS_7max，D_CS_7min]均匀的分为n等分，每等分h＝(D_CS_7max-D_CS_7min)/n，并记子区间[D_CS_7(k)，D_CS_7(k+1)]的中点为x_k+1/2＝x_k+h/2；

步骤4.建立基于Copula函数的联合分布模型

根据给定的五点法和七点法断面平均含沙量x_CS₅和x_CS₇，边缘分布F₁(CS₅)和F₂(CS₇)，计算得到和以及y的反函数和Copula函数c(u,v)；

这里，边缘分布函数及反函数根据函数的不同(见表1)，具体计算表达式也不同，这里就不逐一列出；

步骤5.基于条件期望改正计算含沙量

对于第i个区间，令利用复化辛普森公式，通过数值积分的方法计算得到，给定五点法断面平均含沙量条件x_CS₅下，断面平均含沙量x_CS₇的期望值E(y|x)：

下面，选用宜昌水文站1973～1977年间的65次临底悬沙测验成果，作为分析的样本数据。采用前文所述的K-S检验，对分布函数进行验证，并通过RMSE指标，检验边缘分布的拟合情况，将统计结果列于表2。

从表2可以看出，无论是五点法还是七点法断面平均含沙量样本序列，给定置信水平α＝0.05时，12种分布函数的统计量D_KS均小于临界值16.87％，即在统计意义上，这12种分布函数均可用于描述样本序列的概率分布情况。为了能在所有满足条件的分布函数中选择拟合最好的，本文采用了前述的RMSE准则，即RMSE值越小，该型分布函数对样本的拟合度就越高。12中分布函数中，GNO、GPA和PIII分布的RMSE值均较小，考虑到P-III分布(PIII)的RMSE值最小，同时由前文分析可知宜昌站的水沙一致性较好，而P-III分布又是水文频率分析计算中应用最为广泛的分布形式，故本文选用P-III分布作为五点法和七点法断面平均含沙量样本序列的边缘分布。

表2五点法和七点法断面平均含沙量边缘分布拟合统计表

注：n＝65，α＝0.05时，K-S检验的临界值为

确定七点法和五点法断面平均含沙量样本的边缘分布并进行参数估计后，便可构建基于Copula函数的联合分布。本文选用G-H Copula函数，通过秩次相关系数τ＝0.95，推算copula函数的参数θ＝1/(1-τ)＝20，则构建的Copula函数为：

C_θ＝20(u,v)＝exp{-[(-ln u)²⁰+(-ln v)²⁰]^1/20} (9)

其中，u和v值分别通过PIII分布的累积密度函数f_PIII(g)求得(该密度函数具有标准的函数表达式，文中不再详述)，即：

将七点法断面平均含沙量作为真值，利用前述的基于Copula函数的含沙量改正计算，对五点法断面平均含沙量进行改正，将改正计算的结果列于表3。

从表3可以看出，改正计算前，五点法和七点法断面平均含沙量样本均值的绝对误差为0.092kg/m³，相对误差为5.948％；改正计算后，二者样本均值的绝对误差减小为0.006kg/m³，相对误差减小为0.38％。从样本分布特征值方面看，各种含沙量量级条件下，改正计算后的测验误差均有不同程度的减小，降幅最大的是P＝10％时，相对误差由原来的7.92％，减小到0.92％；降幅最小的是P＝90％时，相对误差由原来的4.061％，减小到1.024％，降幅也超过了3％。

表3基于Copula函数的含沙量改正计算成果表

以上实施例仅仅是对本发明技术方案所做的举例说明。本发明所涉及的河流悬移质泥沙输沙量改正计算方法并不仅仅限定于在以上实施例中所描述的内容，而是以权利要求所限定的范围为准。本发明所属领域技术人员在该实施例的基础上所做的任何修改或补充或等效替换，都在本发明的权利要求所要求保护的范围内。

在上述实施例中，是以五点法和七点法为例进行说明的，本发明所提供的改正计算方法对于其他任意两种测点方式同样适用。并且，采用五点法和三点法来进行测点，会更将简单快速。

Claims (5)

1.一种河流悬移质泥沙输沙量改正计算方法，其特征在于，包括以下步骤：

步骤1.计算断面平均含沙量

将各实测垂线i同一相对水深η的测点流速V_η用面积Ω加权：

$V_{\mathrm{\&eta;}}=\underset{i}{\&Sigma;}{K_i}^{\&prime;}{V}_{\mathrm{\&eta;}-i}=\underset{i}{\&Sigma;}\frac{{\mathrm{\&Omega;}}_i}{\mathrm{\&Omega;}}{V}_{\mathrm{\&eta;}-i}$

将含沙量CS_η用面积Ω与流速V_η二者加权：

${\mathrm{CS}}_{\mathrm{\&eta;}}=\underset{i}{\&Sigma;}{K_i}^{\&prime;\&prime;}{\mathrm{CS}}_{\mathrm{\&eta;}-i}=\underset{i}{\&Sigma;}\frac{{\mathrm{\&Omega;}}_i}{\mathrm{\&Omega;}}\frac{V_{\mathrm{\&eta;}-i}}{V_{\mathrm{\&eta;}}}{\mathrm{CS}}_{\mathrm{\&eta;}-i}$

再根据两组水深权值，分别计算出两种断面含沙量的平均值D_CS_I和D_CS_II，

式中，η为垂线上测点的相对水深值；Ω_i为以垂线i间距的中点为分界线的部分测验断面面积；Ω为全测验断面面积；V_η-i为第i条垂线，在相对水深η处的测点流速；V_η为各垂线在同一相对水深η处的平均流速；CS_η为各垂线在同一相对水深η处的平均测点含沙量；K_i'和K_i”均为面积比系数；I对应第一组水深权值，II对应第二组水深权值；

步骤2.确定边缘分布，计算联合分布函数参数

根据计算得到的D_CS_I和D_CS_II资料序列，通过拟合优度评价指标，确定其边缘分布函数F₁(CS_I)和F₂(CS_II)，并确定联合分布Copula函数C_θ(u,v)的参数θ；

步骤3.区间离散

根据D_CS_II资料序列的最大值D_CS_IImax和最小值D_CS_IImin，将取值区间[D_CS_IImax，D_CS_IImin]均匀的分为n等分，每等分h＝(D_CS_{II max}-D_CS_{II min})/n，并记子区间[D_CS_II(k)，D_CS_II(k+1)]的中点为x_k+1/2＝x_k+h/2；

步骤4.建立基于Copula函数的联合分布模型

根据给定的断面平均含沙量x_CS_I和x_CS_II，边缘分布F₁(CS_I)和F₂(CS_II)，计算得到和以及y的反函数和Copula函数c(u,v)；

步骤5.基于条件期望改正计算含沙量

对于第i个区间，令利用复化辛普森公式，通过数值积分的方法计算得到，给定断面平均含沙量条件x_CS_I下，断面平均含沙量x_CS_II的期望值E(y|x)。

2.根据权利要求1所述的河流悬移质泥沙输沙量改正计算方法，其特征在于：

其中，所述步骤2中：

离差平方和准则为式中，P_i和P_ei分别为理论频率与经验频率；n为序列个数；

Copula函数f(cs_I,cs_II)为f(cs_I,cs_II)＝C_θ(F₁(cs_I),F₂(cs_II))＝C_θ(u,v)；

选用Archimedean Copula函数族中的G-H Copula函数作为联合分布函数，选用RMSE准则评价Copula函数的拟合情况；

G-H Copula函数的数学表达式为：C_θ(u,v)＝exp{-[(-lnu)^θ+(-lnv)^θ]^1/θ}，θ＞1；

参数θ通过Kendall秩相关系数τ求得：θ＝1/(1-τ)。

3.根据权利要求1所述的河流悬移质泥沙输沙量改正计算方法，其特征在于：

其中，所述两组水深权值为用三点法、五点法、七点法测量得到的三组深权值中的任意两组。

4.根据权利要求3所述的河流悬移质泥沙输沙量改正计算方法，其特征在于：

其中，所述第一组水深权值为用五点法测量得到的一组深权值，

所述第二组水深权值为用七点法测量得到的一组深权值。

5.根据权利要求4所述的河流悬移质泥沙输沙量改正计算方法，其特征在于：

其中，所述期望值E(y|x)的计算公式为：

$E\left(y\right|x)=\underset{0}{\overset{1}{\&Integral;}}c(u,v)F_Y^{-1}\left(v\right)dv=\frac{1}{6}\&lsqb;f\left(v_0\right)+4\underset{i=0}{\overset{n-1}{\mathrm{\&Sigma;}}}f\left(v_{i+1/2}\right)+2\underset{i=0}{\overset{n-1}{\mathrm{\&Sigma;}}}f\left(v_i\right)+f\left(v_n\right)\&rsqb;.$