水文预报学习——01

水文预报学习——01

参考书籍：水文预报 第五版 河海大学 包为民

第四章 河道流量演算与洪水预报

一、用试算法求马斯京根法中的 K 、x

1、例子数据来源水文预报 P107

2、计算方法原理

来自于三峡大学水文预报课件

3、Matlab代码实现

clc;clear
%% 数据准备
delta_T=18;%上例图中时间间隔deltaT
%上例图中上游实测流量数据Q上
I_datas=[19900,24300,38800,50000,53800,50800,43400,35100,26900,22400,19600]';
I_datas=I_datas(2:end);
%上例图中下游实测流量数据Q下
Q_datas=[23100,25400,36600,47500,51400,49200,42600,35200,29000,23900,21100]';
Q_datas=Q_datas(1:end-1);
%% 计算过程
%1、河段蓄水量S的计算
delta_S=zeros(length(I_datas),1);
for i=1:length(I_datas)-1
    delta_S(i+1)=1/2*(I_datas(i)+I_datas(i+1))*delta_T-1/2*(Q_datas(i)+Q_datas(i+1))*delta_T;
end
S=cumsum(delta_S);
%2、流量Q的计算
best_cor=0;%存储拟合最好的相关系数
best_x=0;%存储拟合最好的x
for x=0:0.01:0.5
    Q=x*I_datas+(1-x)*Q_datas;
    temp=corrcoef(S(2:end),Q(2:end));%因为S(1)的数据未知，所以从第二项开始拟合
    temp=temp(1,2);
    if temp>best_cor
        best_cor=temp;
        best_x=x;
        best_Q=Q;
    end
end
%% 画出蓄水曲线和用S=KQ求出K值
p=polyfit(best_Q(2:end),S(2:end),1);
K=p(1);
plot(S(2:end),best_Q(2:end))
title('槽蓄曲线')
xlabel('河段蓄水量S')
ylabel('流量Q')
result=[K,best_x]

4、结果展示

K=19.67 x=0.15