python 马斯京根模型 马斯京根法参数

《水文预报》程序

• 写在前边
• 1.相应流量（水位）法河道洪水预报

• 1.1 要求
• 1.2 程序设计（原创）
• 1.3 输入和输出结果

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

• 2.1 要求和思路
• 2.2 程序设计（原创）
• 2.3 结果

• 3.降水保证率计算

• 3.1 要求和内容
• 3.2 程序设计（师给程序，微调）
• 3.3 结果

• 4.蓄满产流模型实验

• 4.1 要求
• 4.2 程序设计（师给程序，注释，微调画图）
• 4.3 结果

• 5.流域断面出口流量过程线计算以及反推单位线

• 5.1 要求及思路
• 5.2 程序设计（Python）（原创）
• 5.3 结果

• 5.超渗产流模型实验

• 5.1 内容
• 5.2 程序设计（师给程序，注释，微调画图）
• 5.3 运行结果

• 6.流域水文模型（蒸发模型实验&降雨径流计算）

• 6.1 要求
• 6.2 程序设计（师给程序，注释，微调画图）
• 6.3 结果

• 写在后边

写在前边

将专业课用目前所知&所学的工具结合在一起，也是件很有意义的事情。
在这里写下代码算是对自己本科最后一个科目的纪念吧❤
以下未标注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(3,1,2)
plot(up_t,up,'b-')
hold on
plot(up_t(locu_max),up_maxx,'r.')
ylabel('up stream','color','b');
m=[];
for ii=1:length(locu_max)
m=[m;ii];
text(up_t(locu_max(ii))+0.1,up_maxx(ii)+0.1,num2str(ii));
end
 
%% 在下游水位（流量)过程图标出极值点
subplot(3,1,3)
plot(down_t,down,'g-')
hold on
plot(down_t(locd_max),down_maxx,'k+')
n=[];
for ii=1:length(locd_max)
text(down_t(locd_max(ii))+0.1,down_maxx(ii)+0.1,num2str(ii));
n=[n;ii];
end
ylabel('down stream','color','g');

%% 输出数据
xlswrite('f',[m,up_t(locu_max),up_maxx],3,'a2')
xlswrite('f',[n,down_t(locd_max),down_maxx],3,'D2')
 
%% 设置寻找极值点的函数
function [maxx,loc_max]=peak(s,mph)
[maxx,loc_max]=findpeaks(s,'minpeakheight',mph);
end

2.3 结果

输入的上游和下游文件格式：

输出结果：

3.降水保证率计算

3.1 要求和内容

1. 降雨资料处理
计算水文统计参数： 、Cv、Cs；
2. 水文频率计算软件应用
① 学习软件运行原理及运行步骤；
② 学习参数调整原则和方法（初值用上一的结果）；
③ 适线调整，并推求不同保证率（25%、50%、75%和95%）下的水文设计变量；
3. 计算结果的合理性检验，并提交计算结果。

3.2 程序设计（师给程序，微调）

clear
clc
w=sort([281.1	463.7	210.1	398.4	419.4	509.1	542.2	424.1	416.6	270.1	622.9	562.1	338.4	485.6	466.3	336.9	375.5	235.3	585.5	295.8	269.6],'descend');%w为样本
%或者w=input('样本系列：')
ave=mean(w);                             %计算均值
cv=std((w)/mean(w))                     %计算cv值
cs=cv*2                                 %初步计算cs值
disp('Ctrl+c结束循环');                   %提示性
n=input('输入样本数：');                 %n为w（）的样本容量
while (1)                                %无限循环
    cv=input('输入Cv值：');              %输入cv cs值进行配线调整
    k=input('输入Cs/Cv值：');
    cs=k*cv;
    
%% 计算理论经验频率
for a=1:18
p=[ 0.01 0.1 0.2  1 2  5 10  20  30 40 50 60 70  80  90 95  99  99.9 (1:n)/(n+1)*100]/100;   %前面十八个为理论频率用,后面为经验频率用
Fp=(cs/2)*gaminv(1-p(a),4/cs*cs,1)-2/cs;                                                     %计算fai值
xp(a)=(Fp*cv+1)*ave;                                                                         %计算Xp值（设计值）   
end
yo=[1.5707962288 0.0370698 -8.364353589*10^(-4) -2.250947174*10^(-4)...                       %yo[]特定数组，进行坐标变换用？  
    6.841218299*10^(-6) 5.824238515*10^(-6) -1.045274970*10^(-6)...
    8.360937017*10^(-8) -3.231081277*10^(-9) 3.657763036*10^(-11)... 
    6.936233982*10^(-13)];

%% 将坐标进行海森几率格纸转换
for b=1:(n+18)                             %整个for循环式p数组所有元素循环              
u=-log(4*p(b)*(1-p(b)));                   %某计算数组u，计算坐标变换用？
for c=1:11;                                
    m(c)=(yo(c)*u^(c-1));                  %理论频率用的p>=50%,进行某计算，变换坐标系？
end
    xo(b)=sqrt(u*sum(m));                  %所有n+18个数参与计算，变换坐标系？  
    if (b<11)
        xo(b)=-xo(b);                       %p<50的理论用频率取负值
    elseif 18<b&b<(18+n/2+1)
         xo(b)=-xo(b);                      %p<50的经验用频率取负值
    end
    L(b)=xo(b)+3.71901645;                  %所有数的某种计算，变换坐标系？
end

%% 用来绘图
plot(L(1:18),xp)                            %输出经验频率曲线
xlabel('P(%)')                              %X轴名  P(%)
ylabel('Q(m^{3}/s)')                        %Y轴名  降水mm)
%set(gca,'ytick',[0:10:1400])
set(gca,'yMinorTick','on')
set(gca,'XTick',L(1:18) )                   %X轴按L()算法伸缩            
set(gca,'XTickLabel',p(1:18)*100)           %按p（1：18）标识伸缩后的X轴
grid on                                     %加网格线   
hold on                                     %保持图形     
title(strcat('Cv=',num2str(cv),';Cs/Cv=',num2str(k)))
plot(L(19:(18+n)),w,'*k')                   %绘制经验点 黑色的*
hold off                                    %取消图形保持  

%% 进行循环判断
mm=input('是否满意？（1\\2）');
if mm==1    
    break
else
    disp('好嘞，继续运行... ...')
    continue
end
end                                         %大循环结束，重新根据cv cs 绘制曲线

%% 以下代码原来求取保证率下对应的设计值
nn=input('是否要求设计值？(1\\2)');
if nn==1
    while(1)
        p=input('请输入所求保证率（小数样式）:');              %输入cv cs值进行配线调整
        Fpp=(cs/2)*gaminv(1-p,4/cs*cs,1)-2/cs;                                                     %计算fai值
        Xp=(Fpp*cv+1)*ave; 
        disp("Xp=")
        disp(Xp)
        disp('Ctrl+c结束循环');  
    end
else
    disp("end")
end

3.3 结果

4.蓄满产流模型实验

4.1 要求

① 掌握蓄满产流模型的运行过程、输入文件的整理、输出文件的格式设定；
② 掌握蓄满模型参数率定的方法；
③ 实用水利程序的开发训练。

4.2 程序设计（师给程序，注释，微调画图）

function RYUSSM(TNum,B,WM,FC,W0,PEfile,Rfile,RGfile,RSfile)
%函数功能：
% 已知流域初始蓄水容量WO，不均匀程度系数B，稳定下渗数值FC，最大蓄水容量WM，TNum个时段降水资料储存在PEfile中
%要求：TNum内各时段产生的径流总量、地下径流量、地表径流量，并且分别输出到Rfile,RGfile,RSfile中
echo off     
clc 

%% Part1：蓄满流模型  Runoff Yield Under Saturated Storage Model
%TNum表示计算时段长度
PE  =zeros(TNum,1);       W  =zeros(TNum,1);                                %净雨量PE和流域蓄水量W
A   =zeros(TNum,1);       R  =zeros(TNum,1);                                %前期影响雨量A和径流量R
RG  =zeros(TNum,1);       RS =zeros(TNum,1);                                %地表径流量RS和地下径流量RG

filename1=PEfile;
fp=fopen(filename1,'rt');                                                   %打开
if fp==-1
   error('Error when opening training set file !');
end
for i=1:TNum
   PE(i)=fscanf(fp,'%f',1);                                                 %打开文件，识别为浮点型数据
end
fclose(fp); 

WMM=WM*(1+B);                                                               %B表示流域蓄水容量分布的不均匀程度
W(1)=W0;


for i=1:TNum
   A(i)=WMM*(1-(1-W(i)/WM)^(1/(1+B)));                                      %将i时刻的蓄水容量转为i时刻的前期影响土湿度A
    
    if PE(i)<0.0                                                            %没有净雨产生            
        R(i)=0.0;
    end
    if PE(i)>=0.0                                                           %产生净雨
        if PE(i)+A(i)<WMM;
            R(i)=PE(i)-WM+W(i)+WM*(1-(PE(i)+A(i))/WMM)^(1+B);               %根据降雨径流关系求得产流量
        else
            R(i)=PE(i)-(WM-W(i));                                           %全流域蓄满时对应的产流量
        end
    end
        
    if PE(i)>FC;                                                            %FC表示流域的实际稳定下渗值
        RG(i)=FC*(R(i)/PE(i));                                              %地下径流量=稳定下渗率*产流面积（即R/PE)
        RS(i)=R(i)-RG(i);                                                   %超过稳定下渗值时产生地面径流=总径流量-地下径流量
    else
        RG(i)=R(i);                                                         %没有产生地面径流
        RS(i)=0.0;
    end
    
    W(i+1)=WMM/(1+B)*(1-(1-(A(i)+PE(i))/WMM)^(1+B));                        %计算流域的蓄水量
end
[R RG RS ]                                                                   %输出值
W

%% Part2 将数据保留在文件里：R——r.txt   RG--rg.txt     RS--rs.txt
filename3=Rfile;
fid=fopen(filename3,'wt');
if fid==-1
   error('Error when opening training set file !');                         
end
fprintf(fid,'%f\n',R);                                      
fclose(fid);

filename4=RGfile;
fid=fopen(filename4,'wt');
if fid==-1
   error('Error when opening training set file !');
end
fprintf(fid,'%f\n',RG);
fclose(fid);

filename5=RSfile;
fid=fopen(filename5,'wt');
if fid==-1
   error('Error when opening training set file !');
end
fprintf(fid,'%f\n',RS);
fclose(fid);

%% 画蓄水容量曲线图
WMM=(1+B)*WM;
Wm=1:(WMM-1)/1000:WMM;
a=1-(1-Wm./WMM).^B;
a_A=1-(1-A./WMM).^B;
mm=size(A);
m=7;
subplot(2,4,1)
plot(a,Wm)
xlabel("a")
ylabel("WM'")
legend("包气带蓄水容量曲线")
for i=1:m
    subplot(2,4,i+1)
    Aa=0:A(i)/100:A(i);
    aa=0:a_A(i)/100:a_A(i);
    for j=1:100
        plot(a_A(i),Aa(j),'r.')                                                             %标记此次降水各时段的前期土湿
        plot(a,Wm)
        hold on
        plot(aa(j),A(i),'r.')
    end
    text(a_A(i),A(i)+5,strcat("第",num2str(i),"个时刻"));
    xlabel("a")
    ylabel("WM'")
    legend("包气带蓄水容量曲线","前期土湿",'Location','southeast')
end

4.3 结果

5.流域断面出口流量过程线计算以及反推单位线

5.1 要求及思路

① 掌握地面径流单位线汇流计算方法和掌握单位线推求方法；
② 实用水利程序的开发训练。
思路
（1）水文分析法
计算原理：叠加假定和倍比假定
计算公式：qi=(Q_(d,i)-∑_(j=2)^(k_2)▒〖r_(d,j) q_(i-j+1) 〗)/r_(d,1)
式中：j表示时段单位线，qi表示单位线上第i个时段对应的径流值，Q_(d,i)表示第i时段在出口断面形成的径流量，r_(d,1)表示第1时段的所形成的净雨量，k2表示在i时段内所对应的净雨场数。
（2）试错法
原理：
将本流域中由一个时段净雨量分析得到的得先作为第一次假定的单位线q（t）
将假定的单位线带入，求解出除了最大试算净雨量外的各时段净雨量产生的部分地面径流过程
将其余时段求得的地面径流过程相加，与实际总径流过程求差值，得到的即为最大净雨量对应的地面径流量
因此利用公式 ▲Q*10/max( r)，得到单位线，
将4中的单位线与假定的对比，若相近，则假定合理，反之则重复1-4，再进行判断

5.2 程序设计（Python）（原创）

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

data=pd.read_excel(r"D:\1-3大四上\2-水文预报\sj\上机实验10\单位线.xlsx", encoding='gbk',index_col="时段")
t=3
r=list(data.iloc[:t,1])                         #r表示净雨过程
n=np.size(data,0)
qq=list(data.iloc[:,0])                         #qq表示题目单位线          

#%%------------------------------定义采用倍比法和叠加法求流量的函数--------------
def qiu_liuliang(r,qq):                         #r表示净雨  qq表示单位线  qtt表示出口断面流量
    t=len(r)                                    #确定输入净雨时段
    n=len(qq)
    qtt=np.zeros(n)
    for i in range(n):
        for j in range(t):
            if i-j>=0 and i-j+1<=n-t+1:
                qt=r[j]*qq [i-j]/10                  
            else:
                qt=0
            qtt[i]=qt+qtt[i]
    return qtt

#%%-----------------------------已知单位线r求出口断面流量过程Q_cal---------------------
Q_cal=qiu_liuliang(r,qq)
#plt.plot(Q_cal)
plt.xlabel('t')
plt.ylabel('Q/m3/s')
plt.show()

#%%-----------------------------水文分析法求流量---------------------------------
Q=list(data.iloc[:,5])
q=np.zeros(n)
q[0]=Q[0]/r[0]*10
QQ=np.zeros(n)
for i in range(n):
    for j in range(t):
        if i-j>=0 and i-j+1<=n-t+1:
            m=r[j]*q[i-j]/10
        else:
            m=0
        QQ[i]=m+QQ[i]                           #求总和 可以优化写
    q[i]=10*(Q[i]-QQ[i])/r[0]

#%%-----------------------------试算法求单位线----------------------------------
u=max(r)                                        #找到最大的净雨值
r[r.index(max(r))]=0                            #把最大值处净雨值替换为0    
qqq=np.linspace(0,3,n)                          #假定的单位线
Q_trial=qiu_liuliang(r,qqq)
qq_trial=(Q-Q_trial)/u*10        
h=max(qq_trial-qqq)
while h>0.0000001:				#判断假定和计算的单位线是否相似
    qqq=0.5*(qq_trial+qqq)
    Q_trial=qiu_liuliang(r,qqq)
    qq_trial=(Q-Q_trial)/u*10
    h=max(qq_trial-qqq)
a=Q_trial+u*qq_trial/10                         #试算求得的出口断面流量

#%%---------------------------------画图------------------------------------------
plt.plot(qq_trial,'b')
plt.plot(a,'b')
plt.plot(q,'r')
plt.plot(Q,'r+')
plt.legend(['q_trial-and-error method','Q_trial-and-error method','q_hydroanalytic methods','Q_hydroanalytic methods'])

5.3 结果

5.超渗产流模型实验

5.1 内容

超渗产流模型计算
①读超渗产流模型程序，标注核心程序段的功能；
②根据输入格式要求，编写输入文件；
③进行超渗模型参数选择和率定；
④整理输出计算结果。

5.2 程序设计（师给程序，注释，微调画图）

function RYUEIM(W0,F0,FC,K,TNum,ER,Rainfile,Rfile) 
%W0表示初始蓄水量，F0表示降雨开始时的下渗率，FC表示稳定下渗速率，K表示霍顿中的参数
%TNum表示降水时段数，ER表示容错率（error），Rainfile表示包含各时段的降水数据的文件，Rfile表示各时段形成的径流文件
%程序功能：
%在已知W0，F0，K，Thum的情况下，通过输入降水数据，求解各时段的径流量（忽略蒸发情况）

echo off     
clc 
%% Step 1：设定相关参数
P   =zeros(TNum,1);       W   =zeros(TNum,1);                               %P表示降水强度；W表示各时段的蓄水量
RR  =zeros(TNum,1);       FF  =zeros(TNum,1);                               %RR表示各时段的径流；FF表示各时段的下渗速率
%绘图相关参数：
PP  =zeros(2*TNum,1);     tt  =zeros(TNum,1);                               %PP用来绘制降水过程；tt表示时间T
Ft  =zeros(TNum,1);                                                         %Ft表示tt对应的下渗曲线上的点
filename1=Rainfile;                                                         
fp=fopen(filename1,'rt');
if fp==-1
   error('Error when opening training set file !');
end
for i=1:TNum
    P(i)=fscanf(fp,'%f',1);
end
fclose(fp);  

F=F0;

%% Step 2:求解降水过程对应的超渗产流过程
for i=1:TNum
%% Step 2.1:迭代求解T(i)和F（i）,其中WB表示下渗曲线上对应的点，（WB，T）不断更新，W0和F0在每个TNum中不变；
% 思路：
% 1.知道上一时刻的下渗量和上一时刻的初始蓄水容量，就知道这一时刻初的蓄水量（也就是知道下渗曲线围成的面积）
% 2.通过迭代求解得到这一时刻初的下渗率对应的T值
% 3.在曲线上再通过T读取这一时刻初的下渗率，并且作为这一时刻的下渗量。
    T=W0/F;                                                                 %给T一个初值，假设T=W0/F，也就是假设曲线围成的面积是矩形
    %最大迭代次数是10，防止运行过久
    for j=1:10
        WB=T*FC+(1-exp(-K*T))*(F0-FC)/K;                                    %霍顿下渗曲线公式；WB表示累积下渗量，求解假设的T对应的WB                                  
        X=abs(WB-W0);                                                       %表示假设的WB和W0的差值
        if X>ER                                                             %若误差大于0.01
            F=F0-K*(WB-FC*T);                                               %求解此时T对应的F
            T=T+(W0-WB)/F;                                                  %更新T
        end
    end
    %求解下渗曲线上的F（i），表示此时刻的下渗率(2种方法）
    F2=FC+(F0-FC)*exp(-K*T);                                                %方法1：直接利用T在下渗曲线上定位点F2
    F=F0-K*(WB-FC*T);                                                       %方法2：老师给的程序，求解T（i）对应的下渗率F（i）【其实化简后就是方法1】                                                                 
    tt(i)=T;
    Ft(i)=F;
%% Step 2.2:判断是否超渗产流
     %Step 2.2.1:求解产流量R（i)，并且确定实际过程中的F（i）
    if P(i)>F                                                               %超渗产流，如果考虑时段的蒸发量E，那么就是 P（i）-E（i）>F(i)
        R=P(i)-F;                                                           %表示时段的产流量，如果考虑时段的蒸发量E，那么就是 P（i）-E（i）-F(i)
        Fs=F;
    else
        Fs=P(i);
        R=0;
    end
      %Step 2.2.2:求解蓄水容量W（i)
    if Fs+W0>50
        W0=50;                                                              %流域最大蓄水容量 50
    else 
        W0=W0+Fs;                                                            %上一时刻的蓄水容量和时段下渗率之和                                                   
    end
    
    RR(i)=R;  W(i)=W0;  FF(i)=Fs;                        
end

[tt RR W FF];

%% Step 3:导出数据
filename3=Rfile;
fid=fopen(filename3,'wt');
if fid==-1
   error('Error when opening training set file !');
end

fprintf(fid,'%f  %f  %f\n',[RR';W';FF']);
fclose(fid);
ttt=[];
pp=[];
for i=1:TNum

    if  i<TNum
        ttt=[ttt,tt(i),tt(i)+0.00000001];
        pp=[pp,P(i),P(i+1)];
    else
        ttt=[ttt,tt(i)];
        pp=[pp,P(i)];
    end
end
%% Step 4:绘制图象
figure(1)
t=1:0.5:25;
f=FC+(F0-FC)*exp(-K*t);
plot(t,f);
hold on
title('下渗曲线');
xlabel('T');
ylabel('ft');
plot(tt,Ft,'r*')                                                           
plot(tt,FF,'ko')
plot(ttt,pp)
legend("下渗曲线","时刻对应的下渗率","实际下渗率（与雨强比大小之后)","降水过程");

5.3 运行结果

6.流域水文模型（蒸发模型实验&降雨径流计算）

6.1 要求

（1）蒸发模型实验
①读蒸发模型程序，标注核心程序段的功能；
②根据输入格式要求，编写输入文件；
③计算并、整理输出文件并打印。
（2）降雨径流计算
①编写具有超渗产流+单位线汇流两个环节的降雨径流计算程序；
②计算某一降雨过程的相应汇流过程；
③提交源程序和计算结果。

6.2 程序设计（师给程序，注释，微调画图）

function RYUEIM(W0,F0,FC,K,TNum,ER,Rainfile,Rfile,Unit_Rfile,TNum_R)
%输入降水P（Rainfile），输入单位线文件Unit_Rfile.txt,输出结果为出口断面流量Rr，存入Rfile中
%采用conv进行求解
echo off     

P   =zeros(TNum,1);       W   =zeros(TNum,1);         
RR  =zeros(TNum,1);       FF  =zeros(TNum,1);
Unit_R = zeros(TNum_R,1);

filename1=Rainfile;
fp=fopen(filename1,'rt');
if fp==-1
   error('Error when opening training set file !');
end
for i=1:TNum
    P(i)=fscanf(fp,'%f',1);
end
fclose(fp);  

filename1=Unit_Rfile;
fp=fopen(filename1,'rt');
if fp==-1
   error('Error when opening training set file !');
end
for i=1:TNum_R
    Unit_R(i)=fscanf(fp,'%f',1);
end
fclose(fp);  

%% 净雨过程计算
F=F0;

for i=1:TNum
    T=W0/F;
    
    for j=1:10
        WB=T*FC+(1-exp(-K*T))*(F0-FC)/K;
        X=abs(WB-W0);
        if X>ER
            F=F0-K*(WB-FC*T);
            T=T+(W0-WB)/F;
        end
    end
    
    F=F0-K*(WB-FC*T);
    if P(i)>F
        R=P(i)-F;
    else
        F=P(i);
        R=0;
    end
    
    if F+W0>50
        W0=50;
    else 
        W0=W0+F;
    end
    
    RR(i)=R;  W(i)=W0;  FF(i)=F;                                            %RR表示净雨
end
%% 单位线汇流
Rr=conv(Unit_R,RR)
RR
%% 输出数据
filename3=Rfile;
fid=fopen(filename3,'wt');
if fid==-1
   error('Error when opening training set file !');
end
fprintf(fid,'%f\n',[Rr']);
fclose(fid);

6.3 结果

写在后边

本科说长不长，说短不短
2022年，难以接受的数字
还无法相信新的一年确实到来了
而我也已经22岁了

回首过去一年，不甘与无奈算是缠绕我最最负面的情绪了吧
忘不了老师的鼓励和帮助，在我emo时候的安慰hhh

或许取舍很重要
没有天赋异禀
只能做时间管理大师

什么都想要
不如逐个去攻破

2022年，请善待我和我的家人
感恩2021
继续前进吧！！！