水文预报程序

写在前边

将专业课用目前所知&所学的工具结合在一起，也是件很有意义的事情。
在这里写下代码算是对自己本科最后一个科目的纪念吧❤
以下未标注python程序的都是用matlab写的，部分是原创，部分是在老师给的代码基础上修改

1.相应流量（水位）法河道洪水预报

1.1 要求

① 采用EXCEL确定马斯京根法参数x；
② 根据参数x；求出马斯京根法方程；
③ 打印试算表和马斯京根预报结果。

1.2 程序设计（原创）

%% 输入数据
filename = "D:\1-3大四上\2-水文预报\sj\上机实验2\马斯京根法.xls";
data=xlsread(filename,'A7:d24');
t=1; %单位时段为12h
Qu=data(:,2);
n=length(Qu);
Qd=data(:,3);
q=data(:,4);
I=Qu;%input
O=Qd-q; %output
bianhua=I-Qd;
 
%% 选取带入的x个数为k^2个，进行k^2次运算
k=3;
kk=k^2;
xx=0.2:0.2/kk:0.4;
poa=[];% 储存试算的预报流量结果
 
for m=1:kk
    delt_W=zeros(1,n);
    W=zeros(n,1);
    x=xx(m);
    subplot(k,k,m)
    for i=2:n
        delt_W(i)=t*(bianhua(i)+bianhua(i-1))/2;
        W(i)=delt_W(i)+W(i-1);
    end
    Q_e=x*I+(1-x)*O; % Q‘
    r=polyfit(Q_e,W,1);
    K=r(1);
    x=xx(m);
    
  %直线拟合R2
    W=polyval(r,QQ);
    R2=1-sum((QQ-Q_e).^2)/sum((Q_e-mean(Q_e)).^2);
    
%开始预报、检验
C0=(0.5*t-K*x)/(0.5*t+K-K*x);
    C1=(0.5*t+K*x)/(0.5*t+K-K*x);
    C2=(-0.5*t+K-K*x)/(0.5*t+K-K*x);
    delt_Q=zeros(1,n);
    Q2=zeros(n,1);
    Q2(1)=O(1);
    for i=2:n
        I1=I(i-1);
        I2=I(i);
        Q1=Q2(i-1)-q(i-1);
        Q2(i)=C0*I1+C1*I2+C2*Q1; %预报值  
    end
    
% 绘制绳套曲线 
    ColorMaker=1:1:n;
    plot(Q_e,W);
    hold on
    scatter(Q_e,W,10,ColorMaker,'o','filled');
    xlabel('Q/(m^3/s)');
    ylabel('W(12h.m^3/s)');
 

% 绘制拟合直线
plot(QQ,W,':');

%计算辅助判断的参数
text(1000,1000,strcat('R2=',num2str(R2)));%R^2
    poa=[poa,Q2];
    delt_Q=Q2-O;  
    delt_Q_max=max(delt_Q)/1000;
    delt_Q_sum=sum(abs(delt_Q))/1000;
    title(strcat('x=',num2str(x),'K=',num2str(round(r(在这里插入代码片1),3)),'误差max',num2str(round(delt_Q_max,3))));
end

1.3 输入和输出结果

输入

输出

2.相应流量（水位）法河道洪水预报

2.1 要求和思路

①处理数据；
②EXCEL绘图；
③添加趋势线及拟合公式；
④检验预报误差，并进行误差统计分析；
⑤打印计算结果。

分析思路
1 相应水位（流量）法：
2 利用程序绘制水位（流量）过程线，找到极值点（在此处寻找波峰），并且提出对应时间、水位和流量序列
3 利用波峰来判断上下游对应同一个传播质点，并且提出数据
4 绘制上下游水位（流量）相关图
5 绘制上游水位（流量）-传播时间关系图

2.2 程序设计（原创）

%% 输入数据
filename = "D:\1-3´óËÄÉÏ\2-Ë®ÎÄÔ¤±¨\sj\ÉÏ»úʵÑé3\ÏÂÓÎ-¼ÐºÓ̲1980.xls";%下游流量文件
down_data=xlsread(filename,'b:F'); %表示文件名
filename2 = "D:\1-3´óËÄÉÏ\2-Ë®ÎÄÔ¤±¨\sj\ÉÏ»úʵÑé3\ÉÏÓÎ-»¨Ô°¿Ú1980.xls";%表示文件名，上游流量文件
up_data=xlsread(filename2,'b:F');
up_t=up_data(:,1);
up_s=up_data(:,2);
up_q=up_data(:,3);
down_t=down_data(:,1);
down_s=down_data(:,2);
down_q=down_data(:,3);
%% 
fig='discharge';% 相应水位（流量）法
up=up_q;% 相应水位法(up_s)或者是相应流量法（up_q)
down=down_q;
mph=2;
[up_maxx,locu_max]=peak(up,mph);
[down_maxx,locd_max]=peak(down,mph);
  
%% 绘制水位（流量)过程图
subplot(3,1,1)
[AX,H1,H2]=plotyy(up_t,up,down_t,down,'plot');
set(get(AX(1),'Ylabel'),'String','Up stream','color','b') ;
set(get(AX(1),'Yaxis'),'color','b') ;
set(get(AX(2),'Ylabel'),'String','Down stream','color','g'); 
set(get(AX(2),'Yaxis'),'color','g') ;
xlabel('Time') ;
title(fig) ;
set(H1,'LineStyle','-','color','b');
set(H2,'LineStyle','-','color','g');
 
%% 在上游水位（流量)过程图标出极值点
subplot