Xu an's Web
0%

水文预报课设

发表于 更新于 分类于 Valine：
接续奋斗

任务

本学期第二个课设为水文预报，根据流域选取标准，选择旬河流域作为产汇流计算流域。

设计内容

  • 洪水场次资料的整理
  • 产流方案设计
  • 汇流方案
  • 精度评价

最终需要给出预报过程图，包括降雨、实测流量、模拟流量。

设计所需数据

由打包文件.DAT给出，其中退水曲线蓄泄系数取 K1=48h ，流域集水面积6448km²。

洪水场次资料的整理

  • 暴雨洪水场次的划分：利用预留的1987——1990划分场次洪水资料，根据已给出的洪水场次计算各场洪水的实测洪峰、实测峰现时间和实测洪水总量。

具体时间如下表： |洪号|流量起止时间| | ----- | --------- | |870513 |05.10-05.20| |870614 |06.11-06.20| |870804 |08.02-08.12| |870903 |09.01-09.09| |880407 |04.01-04.20| |880505 |05.05-05.19| |880705 |07.02-07.16| |880819 |08.17-08.24| |890429 |04.26-05.09| |890711 |07.05-07.22| |890820 |08.14-08.25| |890928 |09.23-10.04| |900501 |04.28-05.13| |900701 |06.29-07.11| |900816 |08.12-08.28|

  • 为了研究暴雨与洪水之间的关系，必须对流量过程线加以分割（次洪划分），又由于不同水源运动规律不同，要把洪水径流分割为地上径流和地下径流（径流成分划分）。
计算各场洪水的洪水总量和次洪径流深

由退水指数方程

\[ Q_t = Q_0e^{\frac{t}{K}} \]

计算各次洪水的退水曲线。

以870513为例，利用蓄泄关系法对洪水径流深进行计算，计算公式如下：

\[ R_0 = 3.6\Delta t(\sum_{i=2}^{n-1}Q_i+\frac{Q_1+Q_n}{2})/A+R_{e末}-R_{e初} \]

其中，\(\Delta t\)为24h，\(R_{e末}-R_{e初}=3.6K(Q_末-Q_初)/A\)。

IMP：不透水面积比例

计算顺序调整:建立产流计算模型

重新研读任务书后，现将设计顺序改为先建立产流模型，明确计算过程，再梳理需要优化率定的参数。

  • 首先根据三层蒸发模型计算出该时段的EU、EL、ED，并计算总蒸散发量E。三层蒸发模型具体计算情况分为以下四种：
  1. 当\(WU+P \ge E_P\)时

\[ E_U = E_P , E_L = 0 , E_D = 0 \]

  1. 当 \(WU+P < E_P , WL \ge C \cdot WLM\)时

\[ E_U = WU+P , E_L = (E_P-E_U)WL/WLM , E_D = 0\]

  1. 当\(WU+P<E_P , C(E_P-E_U) \ge C \cdot WLM\)时

\[E_U = WU+P , E_L = C(E_P-E_U) , E_D = 0\]

  1. 当\(WU+P<E_P , WL < C(E_P-E_U)\)时

\[E_U = WU+P , E_L = WL , E_D = C(E_P - E_U) - E_L\]

式中，\(E_P\)为流域蒸发能力，WL为下层土壤含水量，WLM为下层土壤含水容量，C为蒸发扩散系数，WU为上层土壤含水量，P为降雨量。

  • 算得蒸散发量之后可以得到净雨量PE，通过PE推求产流量。

首先需要计算初始土湿分布以计算流域产流比例，确定后续的产流量。

\[ a = WMM[1-(1-W/WM)^{\frac{1}{1+b}}]\]

在初始土湿为W的条件下，可以建立降雨径流的关系，分为以下两种情况：

  1. $a + PE WMM $

\[R = PE + W - WM + WM(1 - \frac{PE+a}{WMM})^{1+b}\]

  1. \(a + PE > WMM\)

\[R = PE - (WM - W)\]

  • 根据研究，流域坡地上的降雨产流量因产流过程的条件和运动路径不同，受流域的调蓄作用不同，各径流成分在流量过程线上的反映是不一样的。在实际工作中，因为需要按各种径流成分分别进行计算或模拟，因而要对产流量进行水源划分。

本次课程设计需要将其划分为地面径流、壤中流和地下径流以备后续计算。

经过分析概化河槽一侧的土壤剖面,可以得出径流特性可用水箱概念模型来描述和分水源。与蓄满产流类似，由于下垫面的不均匀性，自由水蓄量也存在空间分布不均匀性。因此，应考虑产流面积和自由水蓄量空间分布不均匀的影响。2

根据上述理论，可以列出流域平均自由水蓄积容量的关系：

\[S_m = \frac{S_{mm}}{1+EX}\]

设时段初始自由蓄水量为\(S_1\),其相应纵坐标为AU，则考虑上时段和本时段产流面积有不同的转换有：

\[AU = S_{mm}[1-(1-\frac{S_1 \cdot FR_1/FR}{S_m})^{\frac{1}{1+EX}}]\]

有了上列计算式，即可划分水源。设扣除雨期蒸发后的降雨量为PE，可分为以下两种情况： 1. \(PE+AU<S_{mm}\)

\[RS = FR[PE+\frac{S_1 \cdot FR_1}{FR} - S_m + S_m(1-\frac{PE+AU}{S_{mm}})^{1+EX} ]\]

  1. \(PE+AU \ge S_{mm}\)

\[RS = FR(PE+\frac{S_1 \cdot FR_1}{FR} - S_m)\]

本时段的自由蓄水量为

\[S = \frac{S_1 \cdot FR_1}{FR} + \frac{R-RS}{FR}\]

相应的壤中流和地下径流为

\[ RSS = KSS \cdot S \cdot FR\] \[ RG = KG \cdot S \cdot FR \]

本时段末即下一时段初的自由水蓄量变为

\[S_1 = S(1-KSS-KG)\]

对十五场降水蒸发情况进行迭代计算，得到最终的水源划分结果。 产流计算模型具体代码如下：

产流计算模型 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
%%产流模型
%P时段降雨
%ET时段潜在蒸散发
%%产流参数
WM = 126;
WUM = 63;
WLM = 13;
WDM = 50;
KC = 0.71;%流域蒸散发折算系数
C = 0.17;
B = 2;
IMP = 0.001;
FE = 0.8;%初始土壤水容量折算系数
WMM = WM*(1+B);

W(1,:) = WU(1,:)+WL(1,:)+WD(1,:);
SM = 36;
EX = 0.46;
KG = 0.046;
KKG = 0.995;
KSS = 0.83;
KKSS = 0.06;
SMM = SM*(1+EX);
S1(1,:) = FE*SM;

k = 1;
WU = zeros(20,k);
WL = zeros(20,k);
WD = zeros(20,k);
R = zeros(20,k);
WU(1,:) = FE*WUM;
WL(1,:) = FE*WLM;
WD(1,:) = FE*WDM;
n = 0;
for i = 1:20
    if E0(i,k) > 0
        n = n+1;
    end
end
%%三层蒸发模型
for i = 1:n
    Ep(i,k) = KC*E0(i,k);
    %第一层
    if WU(i,k) + P(i,k) >= Ep(i,k)
        Eu(i,k) = Ep(i,k);
        El(i,k) = 0;
        Ed(i,k) = 0;
    %第二层
    elseif (WU(i,k) + P(i,k) <Ep(i,k)) && (WL(i,k) >= C*WLM)
        Eu(i,k) = WU(i,k) + P(i,k);
        El(i,k) = (Ep(i,k)-Eu(i,k))*WL(i,k)/WLM;
        Ed(i,k) = 0;
    elseif (WU(i,k) + P(i,k) < Ep(i,k)) && (WL(i,k) > C*(Ep(i,k)-Eu(i,k))) && (WL(i,k) < C*WLM)
        Eu(i,k) = WU(i,k) + P(i,k);
        El(i,k) = C*(Ep(i,k) - Eu(i,k));
        Ed(i,k) = 0;
    %第三层
    elseif (WU(i,k) + P(i,k) < Ep(i,k)) && (WL(i,k) < C*(Ep(i,k)-Eu(i,k)))
        Eu(i,k) = WU(i,k) + P(i,k);
        El(i,k) = WL(i,k);
        Ed(i,k) = C*(Ep(i,k) - Eu(i,k)) - El(i,k);
    end
    E(i,k) = Eu(i,k)+El(i,k)+Ed(i,k);
    PE(i,k) = P(i,k) - E(i,k);
    
    a(i,k) = WMM*(1-(1-W(i,k)/WM)^(1/(1+B)));
    %计算i时刻RS，RSS，RG,R
    %S为本时段的自由蓄水量
    if PE(i,k) > 0%产流
        %透水面积产流计算
        if a(i,k)+PE(i,k)<=WMM
            R(i,k) = PE(i,k)+W(i,k)-WM+WM*(1-(PE(i,k)+a(i,k))/WMM)^(B+1);
        else
            R(i,k) = PE(i,k)-(WM-W(i,k));
        end
        %不透水面积产流计算
        RIM(i,k) = P(i,k) * IMP;
        FR(i,k) = R(i,k)/PE(i,k);
        if i == 1
            FR1 = 1-(1-W(i,k)/WMM)^B;
            AU(i,k) = SMM*(1-(1-S1(i,k)*FR1*FR(i,k)/SM)^(1/(1+EX)));%S1(1,k)为时段初始自由水蓄量
        else
             AU(i,k) = SMM*(1-(1-S1(i,k)*FR(i,k)*FR(i,k)/SM)^(1/(1+EX)));%S1(1,k)为时段初始自由水蓄量
        end
        if PE(i,k)+AU(i,k)<SMM
            RS(i,k) = FR(i,k)*(PE(i,k)+S1(i,k)-SM+SM*(1-(PE(i,k)+AU(i,k))/SMM)^(EX+1)) + RIM(i,k);
            RSS(i,k) = (SM - SM*(1-(PE(i,k)+AU(i,k))/SMM)^(EX+1))*KSS*FR(i,k);
            RG(i,k) = (SM - SM*(1-(PE(i,k)+AU(i,k))/SMM)^(EX+1))*KG*FR(i,k);
        else
            RS(i,k) = FR(i,k)*(PE(i,k)+S1(i,k)*FR(i-1,k)/FR(i,k)-SM) + RIM(i,k);
            RSS(i,k) = (PE-SM+S1(i,k))*FR(i,k);
            RG(i,k) = SM*KG*FR(i,k);
            S1(i,k) = (1-KSS-KG)*SM;
        end
    else%不产流
        R(i,k) = 0;
        FR(i,k) = R(i,k)/PE(i,k);
        FRt(i,k) = 1 - (1-W(i,k)/WM)^(B/(1+B));
        RS(i,k) = 0;%地面径流
        RSS(i,k) = S1(i,k)*KSS*FR(i,k);%壤中流
        RG = S1(i,k)*KG*FRt(i,k);%地下径流
        S1(i+1,k) = (1-KSS-KG)*S1(i,k);
    end

    %计算土壤含水量的变化
    WU(i+1,k)=PE(i,k)+WU(i,k)-R(i,k);
    if WU(i+1,k)>=WUM
        WL(i+1,k)=WU(i+1,k)-WUM+WL(i,k);
        if WL(i+1,k)>=WLM
            WD(i+1,k)=WL(i+1,k)-WLM+WD(i,k);
            WL(i,k)=WLM;
            if WD(i+1,k)>WDM
                WD(i+1,k)=WDM;
            end
        end
        WU(i+1,k)=WUM;
    end
    W(i+1,k)=W(i,k)+PE(i,k)-R(i,k);
end

%水量平衡检验
x1 = sum(E(:,k));
x2 = sum(P(:,k));
x3 = sum(PE(:,k));
x = x3-(x2-x1);

if x < 0.001
    Y = ['误差为',num2str(x),',符合水量平衡'];
    disp(Y);
else
    disp('错误！不符合水量平衡 。');
end

进行汇流计算

由于不同水源流速等因素对汇流的影响不同，故需要使用不同的汇流方法计算结果。通常，由地面净雨直接推得的地面径流使用单位线进行计算，地下径流汇流使用线性水库方法计算。

直接应用优选的汇流参数和已给出的单位线对场次暴雨洪水进行单位线计算。壤中流和地下径流用线性水库模型进行求解。计算汇流具体代码如下：

汇流计算 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
%%单位线汇流计算
A = 6448;
U = A/(3.6*24);
k = k;
q = length(RS(:,k));
for j = 1:q
    for i = 1:n
        danweixian(i+j-1,j+2) = danweixian(i,1) * RS(j,k)/10;
    end
end
for i = 1:n
    danweixian(i,j+3) = sum(danweixian(i,3:j+2));
end

m = length(danweixian(:,j+3));
Qend = zeros(m,k);
Qg = zeros(m,k);
Qss = zeros(m,k);
%%线性水库汇流
Qg(1,k) = RG(1,k);
for i = 1:n
    %Qg(1,k) = RO(1,k);
    Qg(i+1,k) = RG(i,k)*(1-KKG)*U+Qg(i,k)*KKG;
end
Qss(1,k) = RSS(1,k);
for i = 1:n
    Qss(i+1,k) = RSS(i,k)*(1-KSS)*U+Qss(i,k)*KSS;
end
for p = 1:m
    Qend(p,k) = danweixian(p,j+3)+Qg(p,k)+Qss(p,k);
end

按时段进行计算，得出每时段预报结果图

已经完成各个时段预报结果的初步计算，现在存在的问题有，时段8和时段13洪峰量很底且峰现时间出现了错位，需要进一步研究出现该问题的原因。

现将检验所用代码插入如下：

预报结果检验 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
%%退水曲线检验
A = 6448;
K = 48;
t = 1;
k = 8;
m = 8;
Qtjy(1,k) = Qend(m,k);
R = 0;
if Qend(1,k) > Qend(2,k) || Qend(1,k) == 0
        d = 2;
else
        d = 1;
end

while Qtjy(t,k) > Qend(d,k)
    t = t+1;
    Qtjy(t,k) = Qtjy(t-1,k) * exp(-1/K);
end

l = length(Qtjy(:,k));

for i = 1:m
    Qtsjy(i,k) = Qend(i,k);
end
for i = m:m+t-1
    Qtsjy(i,k) = Qtjy(i-m+1,k);
end

%R(1,k) = sum(Qend(2:m,k))+sum(Qtsjy(m+1:m+t-1,k));
R0jy(1,k) = 3.6*24*sum(Qtsjy(1:m+t-1,k))/A;
%R0jy(1,k) = 3.6*24*sum(Qend(1:m,k))/A;
%R0jy(1,k) = 3.6*24*(R(1,k)+(Qend(1,k)+Qtsjy(l,k))/2)/A + K*3.6*(Qend(m,k)-Qend(d,k))/A;

%RO(1,k) = R0jy(1,k) - 3.6*24*(Qts(1,k)+Qts(t-1,k))*(m+t-1-2)/(2*A);

%%作图
plot(Qend(1:m,k));
xlim([1 m+t]);
hold on;
plot(m:m-1+t,Qtjy(1:t,k));
title("场次一退水曲线");
ylabel("Q/(m³/s)");
xlabel("t");
line([d d],[0 Qtsjy(d,k)],'LineStyle','--');
line([m-d+t m-2+t],[0 Qtsjy(m-d+t,k)],'LineStyle','--');
line([2 m-d+t],[Qtsjy(d,k) Qtsjy(m-d+t,k)],'LineStyle','--');

在最终绘制对比图代码如下：

作图 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
%%画图
%%绘制实测流量过程、预测流量过程和流域平均降雨量过程示意图
k = 8;
m = 8;
%RS(m,k) = 0;
yyaxis left;%激活左边的轴
bar(1:m,RS(1:m,k));
set(gca,'YDir','reverse'); %反转y轴
xlabel("日期");
ylabel("left 降雨量(mm)");
ylim([0 50]);

yyaxis right;%激活右边的轴
plot(1:m,Q(1:m,k));
hold on
plot(1:m,Qtsjy(1:m,k));%Qtsjy
ylabel("right 流量(m³/s)");
xlim([1,m]);
ylim([0 3000]);
xticks(1:1:m);
title("场次二预报结果")
%{
xticklabels(["5-10","5-11","5-12","5-13","5-14","5-15","5-16","5-17","5-18","5-19","5-20"]);
%}

在分析结果后，发现整个预报的参数还有可以率定的余地。在不改变所有计算初值的情况下，现在需要研究各个参数对预报结果的影响，以期改善预报结果。

尝试人工试错优选后，发现参数的改变会出现不好的预报变好了但是好的预报又变差了，所以现在考虑编程，使用Rockenbrock函数法进行参数优选。

[Rockenbrock函数]http://t.csdnimg.cn/3ETMb

7.4正式完成课设。

总结

本次课设的重点和难点并不在于如何编程和设置数据结构，因为一切使用的公式书中都有给出。重点在于如何通过改变参数使结果能达到一个较好的预期，真正的任务量在于调参，算是一种全新的设计体验。

  1. K为一确定常数，由\(W_t=KQ_t\)确定。该式表明，当泄流流量恒定为\(Q_t\)，K是泄完蓄水量\(W_t\)所需的时间。由于蓄量分布在流域上，距出口断面的距离远近不同。汇集时间大小不等，其平均汇集时间应等于K。从这个意义上讲，K又可以解释为流域水流平均汇集时间。↩︎

  2. 《水文预报（第五版）》，主编：包为民，2017↩︎