新安江模型程序C 代码

合集下载

新安江模型程序C代码

新安江模型程序C代码新安江模型程序C++代码以下是类的声明：class XinanjiangModel{private:// FORCINGdouble *m_pP; // 降水数据double *m_pEm; // 水面蒸发数据//long m_nSteps; // 模型要运行的步长(一共m_nSteps步) long steps;// OUTPUTdouble *m_pR; // 流域内每一步长的产流量(径流深度) double *m_pRs; // 每一步长的地表径流深(毫米) double *m_pRi; // 每一步长的壤中流深（毫米）double *m_pRg; // 每一步长的地下径流深(毫米) double *m_pE; // 每一步长的蒸发(毫米)double *m_pQrs; // 流域出口地表径流量double *m_pQri; // 流域出口壤中流径流流量double *m_pQrg; // 流域出口地下径流量double *m_pQ; // 流域出口的总流量double m_U; // for 24h. U=A(km^2)/3.6/delta_t// SOILdouble *m_pW; // 流域内土壤湿度double *m_pWu; // 流域内上层土壤湿度double *m_pWl; // 流域内下层土壤适度double *m_pWd; // 流域内深层土壤湿度double m_Wum; // 流域内上层土壤蓄水容量double m_Wlm; // 流域内下层土壤蓄水容量double m_Wdm; // 流域内深层土壤蓄水容量，WDM=WM-WUM-WLM // EVAPORATIONdouble *m_pEu; // 上层土壤蒸发量（毫米）double *m_pEl; // 下层土壤蒸发量（毫米）double *m_pEd; // 深层土壤蒸发量（毫米）//runoffdouble *RF;// PARAMETERdouble m_Kc; // 流域蒸散发能力与实测蒸散发值的比double m_IM; // 不透水面积占全流域面积之比double m_B; // 蓄水容量曲线的方次，小流域（几平方公里）B0.1左右// 中等面积（平方公里以内）.2~0.3，较大面积.3~0.4 double m_WM; // 流域平均蓄水容量（毫米）(WM=WUM+WLM+WDM) double m_C; // 流域内深层土壤蒸发系数，江南湿润地区：0.15-0.2，//华北半湿润地区：.09-0.12double m_SM; //自由水蓄水容量double m_EX; //自由水蓄水容量～面积分布曲线指数double m_KG; //地下水日出流系数double m_KI; //壤中流日出流系数double m_CG; //地下水消退系数double m_CI; //壤中流消退系数double *m_UH; // 单元流域上地面径流的单位线double m_WMM; // 流域内最大蓄水容量double m_Area; // 流域面积int m_DeltaT; // 每一步长的小时数int m_PD; // 给定数据，用以判断是否时行河道汇流计算public:XinanjiangModel(void);~XinanjiangModel(void);// 初始化模型void InitModel(long nSteps, double Area,int DeltaT, int PD, char *ForcingFile);// 设置模型参数void SetParameters(double *Params);// 运行新安江模型void RunModel(void);// 保存模拟结果到文件void SaveResults(char *FileName);// 记录出流数据，用以作图分析void Runoff(char *runoff);private:// 进行汇流计算，将径流深度转换为流域出口的流量void Routing(void);};以下是类的定义#include"stdafx.h"#include"xinanjiangmodel.h"#include#include#includeusing namespace std;#include"math.h"#include"stdio.h"#include"conio.h"XinanjiangModel::XinanjiangModel(void){this->m_pP = NULL;this->m_pEm = NULL;this->m_pE = NULL;this->m_pEd = NULL;this->m_pEl = NULL;this->m_pEu = NULL;this->m_pW = NULL;this->m_pWd = NULL;this->m_pWl = NULL;this->m_pWu = NULL;this->m_pR = NULL;this->m_pRg = NULL;this->m_pRi = NULL;this->m_pRs = NULL;this->m_pQ = NULL;this->m_pQrg = NULL;this->m_pQri = NULL;this->m_pQrs = NULL;}XinanjiangModel::~XinanjiangModel(void) { delete[] this->m_pP;delete[] this->m_pEm;delete[] this->m_pE;delete[] this->m_pEd;delete[] this->m_pEl;delete[] this->m_pEu;delete[] this->m_pW;delete[] this->m_pWd;delete[] this->m_pWl;delete[] this->m_pWu;delete[] this->m_pR;delete[] this->m_pRg;delete[] this->m_pRi;delete[] this->m_pRs;delete[] this->m_pQ;delete[] this->m_pQrg;delete[] this->m_pQrs;delete[] this->m_pQri;}// 初始化模型void XinanjiangModel::InitModel(long nSteps, double Area, int DeltaT,int PD, char* ForcingFile){FILE * fp;int i;this->m_nSteps = nSteps;this->steps = this->m_nSteps + 18;// 驱动数据this->m_pP = new double[this->steps];this->m_pEm = new double[this->steps];// 模型输出，蒸散发项this->m_pE = new double[this->steps];this->m_pEd = new double[this->steps];this->m_pEl = new double[this->steps];this->m_pEu = new double[this->steps];// 模型输出，出流项，经过汇流的产流this->m_pQrg = new double[this->steps];this->m_pQrs = new double[this->steps];this->m_pQri = new double[this->steps];this->m_pQ = new double[this->steps];// 模型输出，产流项this->m_pR = new double[this->steps];this->m_pRg= new double[this->steps];this->m_pRi= new double[this->steps];this->m_pRs = new double[this->steps];// 模型状态量，土壤湿度this->m_pW = new double[this->steps];this->m_pWd = new double[this->steps];this->m_pWl = new double[this->steps];this->m_pWu = new double[this->steps];//runoff值this->RF = new double[this->steps];for(i = 0;isteps;i++ ){// 驱动数据this->m_pP [i] = 0.00;this->m_pEm [i] = 0.00;// 模型输出，蒸散发项this->m_pE [i] = 0.00;this->m_pEd [i] = 0.00;this->m_pEl [i] = 0.00;this->m_pEu [i] = 0.00;// 模型输出，出流项，经过汇流的产流this->m_pQrg[i] = 0.00; this->m_pQrs[i] = 0.00;this->m_pQri[i] = 0.00;this->m_pQ[i] = 0.00;// 模型输出，产流项this->m_pR [i] = 0.00;this->m_pRg [i] = 0.00;this->m_pRi [i] = 0.00;this->m_pRs [i] = 0.00;// 模型状态量，土壤湿度this->m_pW [i] = 0.00;this->m_pWd[i] = 0.00;this->m_pWl[i] = 0.00;this->m_pWu[i] = 0.00;}this->m_Area = Area;this->m_DeltaT = DeltaT;this->m_PD = PD;this->m_U = this->m_Area/(3.6 * this->m_DeltaT);// Forcing文件格式：第一列：降水（单位毫米）空格第二列水面蒸发（毫米）if((fp = fopen(ForcingFile,"r")) == NULL) {printf("Can not open forcing file!\n");return; }for(i = 0;im_nSteps;i++ ){ fscanf(fp,"%lf%lf",&(this->m_pP[i]),&(this->m_pEm[i])); } fclose(fp);}// 设置模型参数void XinanjiangModel::SetParameters(double* Params){this->m_Kc = Params[0]; // (1) 流域蒸散发能力与实测水面蒸发之比this->m_IM = Params[1]; // (2) 流域不透水面积占全流域面积之比this->m_B = Params[2]; // (3) 蓄水容量曲线的方次this->m_Wum = Params[3]; // (4) 上层蓄水容量this->m_Wlm = Params[4]; // (5) 下层蓄水容量this->m_Wdm = Params[5]; // (6) 深层蓄水容量this->m_C = Params[6]; // (7) 深层蒸散发系数this->m_SM = Params[7]; // (8)自由水蓄水容量this->m_EX = Params[8]; // (9)自由水蓄水容量～面积分布曲线指数this->m_KG = Params[9]; // (10)地下水日出流系数this->m_KI = Params[10]; // (11)壤中流日出流系数this->m_CG = Params[11]; // (12)地下水消退系数this->m_CI = Params[12]; // (13)壤中流消退系数this->m_WM = this->m_Wum + this->m_Wlm + this->m_Wdm;this->m_WMM = this->m_WM * (1.0 + this->m_B)/(1.0 - this->m_IM); }// 运行新安江模型void XinanjiangModel::RunModel(void){long i;// 模型的状态变量double PE; // > 0 时为净雨量;< 0 为蒸发不足量（mm）double Ep; //m_Kc * m_pEm[i]double P;double R; // 产流深度，包括地表径流、壤中流和地下径流（mm）double RB; // 不透水面上产生的径流深度（mm）double RG; // 地下径流深度（mm）double RI; // 壤中流深度（mm）double RS; // 地表径流深（mm）double A; //土壤湿度为W时土壤含水量折算成的径流深度（mm）double E = 0.0; // 蒸散发(mm)double EU = 0.0; // 上层土壤蒸散发量（mm）double EL = 0.0; // 下层土壤蒸散发量（mm）double ED =0.0; // 深层土壤蒸散发量（mm）double S;double FRo;double FR;double MS;double AU;double WU = 5.0; // 流域内上层土壤湿度double WL = 55.0; // 流域内下层土壤适度double WD = 40.0; // 流域内深层土壤湿度double W = 100.0;double So = 5.0;MS = m_SM * (1 + m_EX);FRo = 1 - pow((1 - So/MS),m_EX);for(i = 0;im_nSteps;i++ ){// ——————蒸散发计算————————————//RB = m_pP[i] * m_IM; // RB是降在不透水面的降雨量P = m_pP[i] * (1 - m_IM);Ep = m_Kc * m_pEm[i];if ((WU + P)>= Ep){EU = Ep; EL = 0; ED = 0; }else if((WU + P)<ep)< bdsfid="342" p=""></ep)<>{EU = WU + P;if(WL>= (m_C * m_Wlm)){ EL = (Ep - EU) * WL/m_Wlm; ED = 0; }else if ((m_C * (Ep - EU))<= WL&&WL<(m_C * m_Wlm)) { EL = m_C * (Ep - EU); ED = 0; }else if (WL<="" p="">{ EL = WL; ED = m_C * (Ep - EU) - EL; }}E = EU + EL + ED;PE = P - E;/* ———————蒸散发计算结束——————————— *///——————子流域产流量计算————————————// if(PE<= 0){ R = 0.00; W = W + PE; }else{A = m_WMM * (1 - pow( (1.0 - W/m_WM), 1.0/(1 + m_B) ) );// 土壤湿度折算净雨量+降水后蒸发剩余雨量<流域内最大含水容量if((A + PE)m_WMM){// 流域内的产流深度计算R = PE /* 降水蒸发后的剩余量*/+ W /* 流域内土壤湿度*/+ m_WM * pow((1 - (PE + A)/m_WMM),(1 + m_B))- m_WM /* 减去流域平均蓄水容量（m_WM:参数）*/+ RB; /* 不透水面上产生的径流*/}// 土壤湿度折算净雨量+降水后蒸发剩余雨量<流域内最大含水容量else{// 流域内的产流深度计算R = PE /* 降水蒸发后的剩余量+ W /* 流域内土壤湿度*/- m_WM /* 减去流域平均蓄水容量*/+ RB; /* 不透水面上产生的径流*/}}//三层蓄水量的计算: WU, WL, WDif(WU + P - EU - R <= m_Wum){ WU = WU + P - EU - R; WL = WL - EL; WD = WD – ED; }else{WU = m_Wum;if(WL - EL + ( WU + P - EU - R - m_Wum ) <= m_Wlm ){WL = WL – EL + ( WU + P - EU - R - m_Wum );WD = WD - ED;}else{WL = m_Wlm;if(WD - ED + WL - EL + ( WU + P - EU - R - m_Wum )- m_Wlm <= m_Wdm )WD = WD - ED + WL - EL+ (WU + P - EU - R - m_Wum ) - m_Wlm;elseWD = m_Wdm;}}W = WU + WL + WD;////三水源划分汇流计算if(PE>0){FR = (R - RB) / PE;AU = MS * (1 - pow((1 - So * FRo/FR/m_SM),1/(1 + m_EX)));if(PE + AU<ms)< bdsfid="408" p=""></ms)<>RS = FR * (PE + So * FRo/FR - m_SM + m_SM * pow((1 - (PE + AU) / MS),m_EX + 1));else if(PE + AU>= MS)RS = FR * ( PE + So * Fro / FR - m_SM );S = So * Fro / FR + ( R – RS ) / FR;RI = m_KI * S * FR;RG = m_KG * S * FR;RS += RB;R = RS + RI + RG;So = S * ( 1 - m_KI - m_KG );FRo = FR;}else{S = So;FR = 1 - pow((1 – S / MS) , m_EX );RI = 0.00;RG = 0.00;So = S * ( 1 - m_KI - m_KG );RS = RB;R = RS + RI + RG;FRo = FR;}////三水源划分计算结束/* 以下部分是状态量：总蒸发量、上、下和深层土壤的蒸发的保存*//* 1 */this->m_pE[i] = E; // 当前步长的蒸发（模型重要输出）/* 2 */this->m_pEu[i] = EU; // 当前步长上层土壤蒸发/* 3 */this->m_pEl[i] = EL; // 当前步长下层土壤蒸发/* 4 */this->m_pEd[i] = ED; // 当前步长深层土壤蒸发/* 5 */this->m_pW[i] = W; // 当前步长流域平均土壤含水量/* 6 */this->m_pWu[i] = WU; // 当前步长流域上层土壤含水量/* 7 */this->m_pWl[i] = WL; // 当前步长流域下层土壤含水量/* 8 */this->m_pWd[i] = WD; // 当前步长流域深层土壤含水量/* 9 */this->m_pRg[i] = RG; // 当前步长流域地下径流深度/* 10 */this->m_pRi[i] = RI; // 当前步长流域壤中流深度/* 11 */this->m_pRs[i] = RS; // 当前步长流域地表径流径流深度/* 12 */this->m_pR[i] = R; // 当前步长的总产流径流深度}this->Routing();}// 保存模拟结果到文件void XinanjiangModel::SaveResults(char* FileName){int i;FILE * fp;if((fp = fopen(FileName,"w")) == NULL){printf("Can not create output file!\n");return;}fprintf(fp," - -- -- ---- - - -- -- ---- - ------ --- -- ------ - - -- ---- ----- -- - - ------- -- -- -\n");fprintf(fp," E(mm) EU(mm) EL(mm) ED(mm) W(mm) WU(mm) WL(mm) WD(mm) R(mm) RS(mm) RI(mm) RG(mm) Q(m3/d) QS(m3/d) QI(m3/d) QG(m3/d)\n");fprintf(fp," - -- -- -- - - --- -- -- - - -- ----- -- - - -- -- -- - - -- -- --------- - - -- --------- ---\n");for(i = 0;isteps;i++ ){ fprintf(fp,"%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3l f%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf\n",this->m_pE[i],this->m_pEu[i],this->m_pEl[i],this->m_pEd[i], this->m_pW[i],this->m_pWu[i],this->m_pWl[i],this->m_pWd[i],this->m_pR[i],this->m_pRs[i],this->m_pRi[i],this->m_pRg[i], this->m_pQ[i],this->m_pQrs[i],this->m_pQri[i],this->m_pQr g[i]);}fclose(fp);}// 进行汇流计算，将径流深度转换为流域出口的流量void XinanjiangModel::Routing(void){///////////// 地面径流汇流计算：单位线法/////////////////////// int i,j;double B[10000] = {0.00};if (this->m_PD == 1){double UH[] ={3.71,12.99,38.96,94.63,131.74,154.00,166.99,176.27,178.12, 172.55,146.58, 90.91,53.80, 31.54,18.55, 9.27, 3.71,0.00};for(i = 0;im_nSteps;i++ ){for(j = 0;j<18;j++ ){ B[i + j] += this->m_pRs[i] * UH[j]/10.0; }}}else{double UH[] ={ 7.18,23.38,63.20,143.10,221.75,365.18,447.40,491.29, 506.93,504.82,468.46,388.56,309.91,166.49,84.26,40.37,17.56 ,3.46};for(i = 0;im_nSteps;i++ ){for(j = 0;j<18;j++ ){ B[i + j] += this->m_pRs[i] * UH[j]/10.0; }}}for(i = 0;isteps;i++ )this->m_pQrs[i] = B[i];///// 壤中流汇流计算:线性水库for(i = 1;isteps;i++ ) {this->m_pQri[i] = this->m_CI * this->m_pQri[i - 1] + (1.0 - this->m_CI) * this->m_pRi[i] * this->m_U; }///// 地下径流汇流计算：线性水库for(i = 1;isteps;i++ ) {this->m_pQrg[i] = this->m_pQrg[i - 1] * this->m_CG + this->m_pRg[i] * (1.0 - this->m_CG) * this->m_U; } //////单元面积总入流计算for(i = 0;isteps;i++ ){ this->m_pQ[i] = this->m_pQrs[i] + this->m_pQri[i] + this->m_pQrg[i]; }}void XinanjiangModel::Runoff(char * runoff){int i;ofstream outfile;outfile.open(runoff);if(outfile.is_open()){for(i = 0;isteps;i++ ){ outfile<<setprecision(3)<<setiosflags(ios::fixed)<m_pQ[i] <<=""></setprecision(3)<<setiosflags(ios::fixed)<outfile.close();}以下是main()函数语句int _tmain(int argc, _TCHAR * argv[]){ long nSteps = 942;int DeltaT = 24;double Area1 = 1603;XinanjiangModel Model1;Model1.InitModel(nSteps, Area1,DeltaT,1,"LFForcingfile.txt");//模型参数/*Kc,IM,B,m_Wum,Wlm,Wdm,C,SM,EX,KG,KI,CG,CI double Params1[] = {0.50,0.01,0.30,10,60,40,0.18,32,1.2,0.075,0.072,0.94,0.7};Model1.SetParameters(Params1);Model1.RunModel();Model1.SaveResults("来凤站日模型计算结果.txt");Model1.Runoff("LF_Q.txt");Model1.~XinanjiangModel();double Area2 = 2991;XinanjiangModel Model2;Model2.InitModel(nSteps, Area2,DeltaT,0,"YCForcingfile.txt");//模型参数/*Kc,IM,B,m_Wum,Wlm,Wdm,C,SM,EX,KG,KI,CG,CI double Params2[] = {0.75,0.01,0.32,10,60,40,0.18,27,1.2,0.065,0.067,0.96,0.8};Model2.SetParameters(Params2);Model2.RunModel();Model2.SaveResults("file.txt");Model2.Runoff("YC_Q.txt");Model2.~XinanjiangModel();FILE *fp1,*fp2;double Q1[1000],Q2[1000],Q[1000]={0.00};ofstream outfile;outfile.open("Q.txt");if((fp1=fopen("LF_Q.txt","r"))==NULL){ printf("Can not open the file!\n"); return 0;}if((fp2=fopen("YC_Q.txt","r"))==NULL){ printf("Can not open the file!\n"); return 0; }if(outfile.is_open()){for(int i=0;i<960;i++){fscanf(fp1,"%lf",&Q1[i]);fscanf(fp2,"%lf",&Q2[i]);Q[i]=Q1[i]+Q2[i];outfile<<setprecision(3)<<setiosflags(ios::fixed)<<q[i]<<en dl;< bdsfid="568" p=""></setprecision(3)<<setiosflags(ios::fixed)<<q[i]<<endl;<> }fclose(fp1);fclose(fp2);}outfile.close() ;return 0;}。

新安江模型程序核心源代码

%%%新安江模型程序核心源代码function Qr=XAJ_JUN(DAREA,DT,EM,WwU,WwL,WwD,P,S0, FR0, Qrs0, Qrss0, Qrg0,parameter,Qm) % XAJ是新安江的运行程序,用于单纯形和遗传算法调用,也用于新安江模型的预报Imp1=parameter.IMP ; %流域不透水面积比：次洪Kc= parameter.Kc ; %流域蒸散发折算系数：多年总径流量决定WMU=parameter.WMU ; %流域上层蓄水容量WML=parameter.WML ; %流域中层蓄水容量WMD = parameter.WMD ; %流域下层蓄水容量B = parameter.B ; %流域蓄水容量分布曲线指数C = parameter.C ; %流域深层蒸发系数Ex = parameter.Ex; %流域自由水分布曲线指数SM = parameter.SM ; %流域自由水平均蓄水容量Ki = parameter.Ki ; %自由水箱壤中流出流系数Kg = parameter.Kg ; %自由水箱地下水出流系数Cs = parameter.Cs ; %地面水线性水库汇流系数Ci = parameter.Ci ; %壤中流线性水库汇流系数Cg = parameter.Cg ; %地下水线性水库汇流系数Ke = parameter.Ke ; %马斯京根法河段传播时间Xe = parameter.Xe ; %马斯京根法流量比重系数L = parameter.L ; %滞后演算法参数%次洪决定：WM，B，Imp%WwU(0)=WwU;WwL(0)=WwL;WwD(0)=WwD;%由于日模型与次洪模型的计算时段长不同，参数值不能全部通用，但K、WM、WUM、WLM、B、IMP、EX、C与时段长无关，可以直接引用，%Kc SM、KG、KSS、CS、CI、CG与时段长相关，不能直接引用，需要另外率定%junjunzhu-XAJ-MODELU=DAREA/(DT*3.6); %单位转换D=24/DT;KSSD = (1 - (1 - (Ki + Kg)) ^ (1 / D)) / (1 + Kg / Ki); % 'KSSD，ki出流系数KGD消退系数KGD = KSSD * Kg / Ki;%A_WM=A_WUM+A_WLM+A_WDM;%WMM=(1+B).*WM/(1-IMP);Epp=Kc*EM;% PE=P-K.*EM;for T=1:size(EM,1) %% T以时段为单位计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%三层蒸散发计算if (WwU + P(T)) >= Epp(T)EU(T) = Epp(T); %上层蒸发%Epp为EMEL(T) = 0; %中层ED(T) = 0; %下层elseEU(T) = WwU + P(T) ; %'Ww(1) + P为EUEL(T) = (Epp(T) - EU(T)) * WwL / WML; %要求计算的下层蒸发量与剩余蒸散发能力之比不小于深层蒸散发系数cED(T) = 0;if WwL <= (C * WML) %第二层水量小于蒸散发能力if WwL >= C * (Epp(T) - EU(T)) %'要求计算的下层蒸发量与剩余蒸散发能力之比小于深层蒸散发系数cEL(T) = C * (Epp(T) - EU(T)) ;ED(T) = 0;elseEL(T) = WwL;ED(T) = C * (Epp(T) - EU(T)) - EL(T) ;endendendPE(T) = P(T) - EU(T) - EL(T) - ED(T); %%%%%%%%%%%%%%%%%%%%%%%%==========================================%产流计算部分%%%%%%%%%%%%%%%%%%%%%%%%%%%%%===================================== =====Wm0 = WMU + WML + WMD; % '平均蓄水容量W0 = WwU + WwL + WwD; %'初始含水量R = 0;Rimp = 0;Wmm = (1 + B) * Wm0 / (1 - Imp1) ; % 'Imp1不透水面积比，Wmm为蓄水容量极值if PE(T) >0 %Then GoTo 1000 '降雨小于蒸发，B为蓄水容量曲线的指数if abs(Wm0 - W0) <= 0.0001 % 'Wmm为蓄水容量极值A = Wmm;elseA = Wmm * (1 - (1 - W0 / Wm0) ^ (1 / (1 + B))); %'A为与W0对应的在蓄水容量曲线的纵坐标endif (PE(T) + A) < Wmm % '部分产流R = PE(T) - Wm0 + W0 + Wm0 * ((1 - (PE(T) + A) / Wmm) ^ (1 + B));elseR = PE(T) - (Wm0 - W0) ; % '全部产流endif abs(R - PE(T)) <= 0.0001R = PE(T);Rimp = PE(T) * Imp1 ; % '直接径流endWwU = WwU + P(T) - R - EU(T); %% '第一层蓄水变化WwL = WwL - EL(T) ; % '第二层蓄水变化WwD = WwD - ED (T); %'第三层蓄水变化elseWwU = WwU + P(T) - EU(T); %% '第一层蓄水变化WwL = WwL - EL(T) ; % '第二层蓄水变化WwD = WwD - ED(T) ; %'第三层蓄水变化endif WwU > WMU % '由Ww(1) = Ww(1) + P - R－E(1)：E(1)两断Epp和Ww1WwL = WwL + WwU - WMU; % '由Ww(2) = Ww(2) + Ww(1) - WM(1)检查是否超标WwU = WMU; % '纠正if WwL > WMLWwD = WwD + WwL - WML;WwL = WML;endendif WwU < 0WwU = 0;end%'＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝%'汇流计算部分%'＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝%'水源划分X = FR0 ; % 'FR0产流面积if PE(T) <= 0 %'认为单是地下自由水在产流面积上的深为Rs(T) = 0;Rss(T) = S0 * KSSD * FR0 ; %'KSSD，ki，KGD（KG地下水出流）出流系数Rg(T) = S0 * KGD * FR0;S0 = S0 - (Rss(T) + Rg(T)) / FR0 ; % 's表示自由水在产流面积上的平均蓄水深elseFR0 = R / PE(T); % '用流量除以单位面积上的净雨（可以理解为产流深）即得产流面积S0 = X * S0 / FR0 ; % '产流面积变化的影响SS = S0;Q = R / FR0 ; % '为产流面积上的平均值NN = fix(Q / 5) + 1 ; % '每次入流按5毫米分成并取整数NN为了消除前向差分误差Q = Q / NN; % '一天分为CSng(NN)个时段Kssdd = (1 - (1 - (KGD + KSSD)) ^ (1 / NN)) / (1 + KGD / KSSD);Kgdd = Kssdd * KGD / KSSD;Rs(T) = 0;Rss(T) = 0;Rg(T) = 0;% ' SM流域的平均自由水容量Smm = (1 + Ex) * SM ; % ' Smm全流域最大的自由水蓄水容量if Ex < 0.001 ThenSmmf = Smm ; % ' Smmf表示产流面积最大一点的自由蓄水容量elseSmmf = Smm * (1 - (1 - FR0) ^ (1 / Ex)); % ' Ex表示流域自有水容水容量曲线的指数endSmf = Smmf / (1 + Ex); %' Smf表示产流面积上一点的自有水平均蓄水容量深for j = 1:NNif S0 > Smf %'s 表示自由水在产流面积上的平均蓄水深S0 = Smf;endAU = Smmf * (1 - (1 - S0 / Smf) ^ (1 / (1 + Ex)));if Q + AU <= 0Rsd(T) = 0 ; %' 当径流与此时刻的平均蓄水深之和小于0时不产流Rssd(T) = 0;Rgd(T) = 0;S0 = 0;else if Q + AU >= Smmf % ' 当径流与此时刻的平均蓄水深之和大于最大平均蓄水深全面产壤中流Rsd(T) = (Q + S0 - Smf) * FR0 ; % ' Rsd中d为分段的地面流Rssd(T) = Smf * Kssdd * FR0 ; % ' Rsd中d为分段的壤中流Rgd(T) = Smf * Kgdd * FR0 ; % ' Rsd中d为分段的地下径流S0 = Smf - (Rssd(T) + Rgd(T)) / FR0; % ' s表示自由水在产流面积上的平均蓄水深else if Q + AU < Smmf % ' 当径流与此时刻的平均蓄水深之和大于最大平均蓄水深部分产壤中流Rsd(T) = (Q - Smf + S0 + Smf * (1 - (Q + AU) / Smmf) ^ (1 + Ex)) * FR0;Rssd(T) = (S0+ Q - Rsd(T) / FR0) * Kssdd * FR0;Rgd(T) = (S0 + Q - Rsd(T) / FR0) * Kgdd * FR0;S0 = S0 + Q - (Rsd(T) + Rssd(T) + Rgd(T)) / FR0;endendendRs(T) = Rs(T) + Rsd(T) ; % '累计三流Rss(T) = Rss(T) + Rssd(T) ; % '累计Rg(T) = Rg(T) + Rgd(T);clear Rsd Rssd RgdendendOUT=[Rs;Rss;Rg];%Rs=OUT(:,1); Rss=OUT(:,2);Rg=OUT(:,3);Rs(T) = Rs(T) * (1 - Imp1) ; % '扣除不透水面积Rss(T) = Rss(T) * (1 - Imp1);Rg(T) = Rg(T) * (1 - Imp1);%'Qrs = (Rs + Rimp) * U%'Qrss = Rss * U * (1 - Ci) + Qrss0 * Ci%'Qrg = Rg * U * (1 - Cg) + Qrg0 * Cg%'＝＝＝＝＝＝＝＝＝＝◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎××%' '坡面汇流－－－－－－－－－－－汇流%'＝＝＝＝！＝＝＝＝＝＝￥￥￥￥￥￥￥￥￥￥￥￥￥￥￥￥￥￥￥￥￥￥￥￥￥￥￥￥＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＝＝＝＝＝＝＝＝Qrs(T) = (Rs(T) + Rimp) * U * (1 - Cs) + Qrs0 * Cs; % '地面水线性水库汇流系数CS Qrss(T) = Rss(T) * U * (1 - Ci) + Qrss0 * Ci ; % '壤中流线性水库汇流系数CI Qrg(T) = Rg(T) * U * (1 - Cg) + Qrg0 * Cg ; % '地下水线性水库汇流系数Cg Qtr(T) = Qrs(T) + Qrss(T) + Qrg(T);QsN(T) = (Rs(T) + Rimp) * U ; %'地面径流总和QIIGG(T) = Qrss(T) + Qrg(T) ; % '地下和壤中总和Qm(T) = Qtr(T); %马斯金根Qfm = Qtr(T); %非马斯金根Qrs0 = Qrs(T);Qrss0 = Qrss(T);Qrg0 = Qrg(T);Rs0 = Rs(T);Qr=Qtr' ;clear Qrs Qrss Qrg Rs R Rimp Rs Rss Rgend。

新安江模型介绍

新安江模型介绍：三水源新安江模型蒸散发计算采用三层模型；产流计算采用蓄满产流模型；用自由水蓄水库结构将总径流划分为地表径流、壤中流和地下径流三种；流域汇流计算采用线性水库。

模型结构：模型计算：在新安江模型中，流域蒸散发计算没有考虑流域内土壤含水量在面上分布的不均匀性，而是按土壤垂向分布的不均匀性将土层分为三层，用三层蒸散发模型计算蒸散发量。

参数有流域平均张力水容量WM（mm），上层张力水容量UM(mm)，下层张力水容量LM(mm)，深层张力水容量DM(mm)，蒸散发折算系数KC和深层蒸散发扩散系数C。

具体计算为若P+WU>=EP,则EU=EP,EL=0,ED=0;若P+WU<EP,则EU=P+WU;若WL>C*LM,则WL=(EP-EU)WL/LM,ED=0;若WL<C*LM且WL>=C*(EP-EU),则EL=C*(EP-EU),ED=0;若WL<C*LM且WL<C*(EP-EU),则EL=WL,ED=C*(EP-EU)-WL;水源划分中，本小组采用的是三水源划分。

三水源用自由水蓄水库结构解决水源划分问题，自由水蓄水库结构考虑了包气带的垂向调蓄作用，按蓄满产流模型计算出总径流量R，先进入自由水蓄水库调蓄，再划分水源。

模型参数调整：1蒸散发能力折算系数KCKC是影响产流量计算最为重要和敏感的参数，产流计算中KC控制着水量平衡，因此，对水量计算是最重要的。

KC主要反映流域平均高程与蒸发站高程之间差别的影响和蒸发皿蒸散发于陆面蒸散发间差别的影响。

在实际模拟计算中KC值往往变化很大，最后需经模型调试并验证后确定。

2流域平均张力水容量WM流域平均张力水容量WM表示流域干旱程度，分为UM,LM,DM。

根据经验，南方湿润地区WM约为120~150mm，半湿润地区WM约为150~200mm。

3流域蓄水容量—面积分布曲线指数BB值反映划分单元流域张力水蓄水分布的不均匀程度。

matlap新安江三水源模型程序

.新安江模型程序核心源代码function [fit,dc,result]=XAJ(XX)% XAJ是新安江的运行程序,用于单纯形和遗传算法调用,也用于新安江模型的预报% XX是调用的优化参数% fit 返回目标函数的适值% dc返回有效性系数.% result是一个数组,返回格式为[时间,雨量,实测流量,计算流量];.%% $Date: 2005/5/25 $%email:******************.cn% 输入起始值W,WU,WL,WD,QGWU=20;WL=50;WD=10;FR=0.89; S=2; AREA=7547;U=AREA/3.6;W=WU+WL+WD;%输入雨量E,蒸散发能力P,实测流量QSglobal DA TA;TIME=DA TA(:,1);P=DA TA(:,2);EM=DATA(:,3);QS=DATA(:,4);TRSS0=0.3.*QS(1);TRG0=0.4.*QS(1);% 参数处理[num,numvars]=size(XX);% 优化参数A_K=XX(:,1);A_SM=XX(:,2);A_KG=XX(:,3);A_KSS=XX(:,4);A_KKG=XX(:,5);A_KKSS=XX(:,6);A_CS=XX(:,7);A_WUM=XX(:,8);A_WLM=XX(:,9);A_WDM=XX(:,10);A_IMP=XX(:,11);A_B=XX(:,12);A_C=XX(:,13);A_EX=XX(:,14);A_L=XX(:,15);A_WM=A_WUM+A_WLM+A_WDM;for I=1:num %%%% %%% 对每组数计算K=A_K(I);SM=A_SM(I);KG=A_KG(I);KSS=A_KSS(I);KKG=A_KKG(I);KKSS=A_KKSS(I);CS=A_CS(I);WUM=A_WUM(I);WLM=A_WLM(I);WDM=A_WDM(I);WM=WUM+WLM+WDM;IMP=A_IMP(I);B=A_B(I);C=A_C(I);EX=A_EX(I);L=A_L(I);L=round(L);WMM=(1+B).*WM/(1-IMP);M=size(P,1);PE=P-K.*EM;for T=1:M %% T以时段为单位计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %以下为产流计算if PE(T)<0R=0;elseif W>=WMA=WMM;elseA=WMM*(1-(1-W/WM).^(1/(1+B)));endif A+PE(T)>0if A+PE(T)<WMMR=PE(T)-WM+W+WM.*(1-(PE(T)+A)./WMM).^(1+B);elseR=PE(T)+W-WM;endelseR=0;endend %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5 % 以下为蒸发计算zhengfaif PE(T)<0if WU+PE(T)>0EU=K*EM(T);ED=0;EL=0;WU=WU+PE(T);elseEU=WU+P(T);WU=0;if WL>C*WLMEL=(K.*EM(T)-EU).*WL/WLM;WL=WL-EL;ED=0;elseif WL>C.*(K.*EM(T)-EU)EL=C.*(K.*EM(T)-EU);WL=WL-EL;ED=0;elseEL=WL;WL=0;ED=C.*(K*EM(T)-EU)-EL;WD=WD-ED;endendendelseEU=K.*EM(T);ED=0;EL=0;if WU+PE(T)-R<WUMWU=WU+PE(T)-R;elseif WU+WL+PE(T)-WUM>WLMWU=WUM;WL=WLM;WD=W+PE(T)-R-WU-WL;elseWU=WUM;WL=WU+WL+PE(T)-R-WUM;endendendE=EU+EL+ED;W=WU+WL+WD;% 以下为分水计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SMM=(1+EX).*SM;if (PE(T)<=0)|(R<=0)RS=0;RG=S.*KG.*FR;RSS=RG.*KSS./KG;elseX=FR;FR=(R-PE(T).*IMP)./PE(T);S=X.*S./FR;SS=S;Q=R./FR;G=fix(Q./5)+1;Q=Q./G;%KSSD=KSS.^(1/G);KGD=KSSD.*KG./KSS;RS=0;RG=0;RSS=0;for J=1:Gif S>=SMAU=SMM;elseAU=SMM.*(1-(1-S./SM).^(1./(1+EX)));endif AU+Q<SMMRS=(Q-SM+S+SM.*(1-(Q+AU)./SMM).^(1+EX)).*FR+RS;elseRS=(Q+S-SM).*FR+RS;endS=J.*Q-RS./FR+S;RG=S.*KGD.*FR+RG;RSS=S.*KSSD.*FR+RSS;S=J.*Q+SS-(RS+RSS+RG)./FR;endendOUT(T,:)=[RS,RSS,RG];end % 一次数据演算完%以下为汇流计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% RS=OUT(:,1); RSS=OUT(:,2);RG=OUT(:,3);TRS(1)=RS(1).*U;TRSS(1)=TRSS0 ;TRG(1)=TRG0 ;TR(1)=TRS(1)+TRSS(1)+TRG(1);for T=2:MTRS(T)=RS(T).*U;TRSS(T)=TRSS(T-1).*KKSS+RSS(T).*(1-KKSS).*U;TRG(T)=TRG(T-1).*KKG+RG(T).*(1-KKG).*U;TR(T)=TRS(T)+TRSS(T)+TRG(T);endQJ=TR;if L<0 L=0;endfor T=L+2:MQJ(T)=CS.*QJ(T-1)+(1-CS).*TR(T-L);end%以下为目标函数计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% alf=0.6;y1=0;y2=0;n1=1;n2=1;for T=1:Mif QJ(T)>800y1=(QJ(T)-QS(T)).^2+y1;n1=n1+1;elsey2=(QJ(T)-QS(T)).^2+y2;n2=n2+1;endendq0=mean(QS);q1=mean(QJ);y=(y1*alf/n1+y2*(1-alf)/n2)*(1+abs(q0-q1)/q0);fit(I)=y;%以下为(有效性系数)确定性系数计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%f1=sum( (QS-QJ').^2);f2=sum((QS-mean(QS).*ones(M,1)).^2);dq=1-f1/f2;dc(I)=dq;result =[TIME,P,QS,QJ'];end %一组参数计算结束Ifit=-fit'; 遗传算法为了求最大值,在此加负号.dc=dc';。

WinBTOPMC模型和新安江模型的比较

The Comparison Between Winbtopmc And Xin’AnjiangModelXin PengleiGraduate student，College of water resources and environment，Hohai University，Nanjing（210098)AbstractThe Xin’anjiang model is used widely and perfected in China during the last thirty years, and WinBTOMC model is refers to the Blockwise application of TOPMODEL, which needs more high resolution hydrologic data and seems have a better performance. This study compared these two models on a same semi-arid basin, and the result is: the Xin’anjaing model performed much better than the WinBTOPMC model on yearly data while the WinBTOPMC model performs well on monthly metrological data.Key w ords：Xin’anjiang Model，WinBTOPMC model，distributed1. IntroductionXin’anjiang model[1] is used widely and perfectly in china for the last thirty years, however, as a conceptual model, it seems can’t catch up with the world now. Nowadays, there are more and more high resolution digital data available and more and more distributed hydrologic models are flourishing. May be it is high time we modify the Xin’anjiang mdoel to a distributed one.The WinBTOPMC model[2] is refers to the Blockwise application of TOPMODEL, using Muskingham-Cunge flow routing, and because of its sub-division block, it overcome a commonly accepted maximum basin area of 1500 km2 for TOPMODEL applications.This study tried to find out the merit of this WinBTOPMC comparing with the Xin’anjiang model on a semi-arid watershed－the Dongwan basin which located in the south of the Yellow river in China. It is a middle larged watershed with an area about 2700 quare kilometers and the terrain of the region is moderately sloping with soils. The annual average precipitation of the region is 500-800mm which varies greatly among different years, the greatest annual precipitation would be 3 times than that of others. There are 11 precipitation gauging stations and 3 discharge gauging stations.However, the result is not perfect as expected, The Xin’anjiang model performed better on the yearly data; otherwise, the WinBTOPMC model performed well on monthly data.2. Xin’Anjiang ModelThe main feature of Xin’anjiang Model is the concept of runoff formation on repletion of storage, which means that runoff is not produced until the soil moisture content of the aeration zone reaches field capacity, and then the runoff equals the rainfall excess without further loss.2.1The Model structureThe basin is divided into several blocks. The outflow from each block is first simulated and then routed down the channels to the mail basin outlet. Based on the concept of runoff formation on the repletion of storage, the simulation of outflow from each block is consisted of four major parts:(a)the evapotranspiration which generates the deficit of soil storage which is divided into three layers: upper, lower and deep;(b)the runoff production which produces the runoff according to the rainfall and soil storage deficit;(c)the runoff separation which divides the above so determined into three components: surface, interflow and groundwater;(d) the flow routing which transfers the local runoff to the outlet of each block forming the outflow of the block.2.2The Model performanceXin’anjiang Model performs well in this basin,which can be illustrated by the following figures. The blue line shows the observed flood hydrograph, the yellow line shows the discharge from the upper stream, the black line shows the calculated discharge.Fig.1 The flood hydrograph of 1995Fig.2 The flood hydrograph of 1996Fig.3The flood hydrograph of 1997Fig.4 The flood hydrograph of 1998Fig.5 The flood hydrograph of 1999Fig.6 The flood hydrograph of 2000From the figures above, we can see that the simulation results are good, all of the Nash coefficients are greater than 80%.3. WinbtopmodelWinBTOPMC refers to the Blockwise application of TOPMODEL, usingMuskingham-Cunge flow routing. The two primary extensions made in WinBTOPMC are: (a) Overcoming a commonly accepted maximum basin area of 1500 km 2 for TOPMODEL applications by representing large basins as a collection of "blocks" or sub-basins.(b) Inclusion of Muskingham-Cunge flow routing, since the timing of discharge delivery to the outlet becomes increasingly important for large basin sizes.The Model structure:WinBTOPMC retains the assumption of spatial lumping of the water table over a modelled area. However, on the assumption that this lumping is an important factor in limiting the maximum application area for TOPMODEL, the water table is only lumped at the block (sub-basin) scale, rather than the scale of the entire basin. Therefore, for n blocks there will be n equations describing the relationship between local values of the saturation deficit and topographic index and block-average values of saturation deficit,topographic index and m, of the form:where i is the block number of the location (x,y).Fig.7 The vertical profile of each grid The vertical profile of each grid cell is conceptualised as a three zone system: a root zone of fixed capacity that directly receives precipitation and is subject to evapotranspiration, an unsaturated zone that receives all "overflow" from the root zone, and a saturated zone (water table) that receives water from the unsaturatedzone at a rate given by:where K 0 is the vertical unsaturated conductivity,conventionally assumed to be equal to the saturated transmissivity T 0, and S uz is the unsaturated zone storage. Overland flow q of occurs if the saturation deficit minus the unsaturated zone storage is less than zero.In addition, WinBTOPMC permits grid-by-grid spatial variability in soil and land cover properties, as well as climatic forcing (i.e. precipitation and potential evapotranspiration).4. Data PreparationData required include land cover classification; Digital Elevation Model; soil map and other meteorogical datas.4.1 Digital elevation modelThe DEM dataset are obtained from /mgg/topo/globe.html with a horizontal grid spacing of 30 arc seconds. See figure 8.4.2 Land cover map (Matched with IGBP land cover classification scheme)Downloaded the GLOBAL LAND COVER CHARACTERIZATION from:/landdaac/glcc/glcc.ht mlSee figure 9.4.3 soil map(figure 10.)4.4 NDVI (Pathfinder AVHRR Land Data)From the web: ftp:///data/avhrrFig.8 DEMFig.9 Land cover mapFig.10 Soil map4.5 Monthly metrological dataMonthly metrological data are mainly used for calculating the evaporation on the basin. Include Mean daily temperature, Mean diurnal temperature range, vapor pressure, cloud cover, wind speed. : Source form CRU Climatic Datasets (UEA/CRU) and Extra terrestrial Radiation and Possible duration of sunshine : calculate by equation.5. ParameterizationThe specific parameters used byWinBTOPMC are as follows:(a) Saturated transmissivity, T 0 , which describesthe potential rate of lateral flow for a completely saturated soil profile for a given hydraulicgradient.This parameter is obtained through calibration. That is reclassifying the soil classification in order to reduce the numbers of classes and thencalibrate this parameter.(b) A decay coefficient, m , which describes how actual transmissivity decreases when the soil is not completely saturated. And this parameter is obtained by try and error.(c) Manning's roughness parameter, n , which is required by the Muskingham-Cunge flow routing algorithm. And this parameter is obtained by try and error.(d) Rooting depth, S r,max , which is used to represent the vegetation rooting depth, but also integrates the interception capacity of the canopy. And here also calibrate the S r,max for each land cover classes.6. Results6.1 Parameter resoultsT Able1 Basin charactersBasinID m N0 alpha SDbar 0 0.03 0.03 0.01 0.2 1 0.02 0.005 0.01 0.15 2 0.01 0.005 0.01 0.13 3 0.01 0.01 0.01 0.2T Able2 Maximum Storage in the Root zoneID Landcover Area (%)Sramx4 DBF 33.075 M F 5.04 0.56 CS 14.06 8 WS 11.49 9 S 0.19 10 G 0.01 0.1 12 C 21.55 14 C/N 14.590.1T Able3 Soil TransmisivitySoil_Type T0Clay 110 Sand 200 Silt 1506.2 Flood hydrographsFig.11 The flood hydrograph of 1997Fig.12 The flood hydrograph of 1998Fig.13 The flood hydrograph of 1999Fig.14 The flood hydrograph of 2000Fig.15 The flood hydrograph of 1995Fig.16 The flood hydrograph of 1996The first 4 floods are used for calibration, the last 3 for validation. From the figures above, we can see that the simulation results are good but not perfect as Xin’anjiang model.However, the fig.17 and 18 provide us a good detail performance in calibration months.And the validation month are also very good! (Figure 19-20)Fig.17 the flood hydrolograph of AUG. 1998Fig.18 the flood hydrolograph of jul. 2000Fig.19 the flood hydrolograph of aug. 1995Fig.20 the flood hydrolograph of sep. 1996 AcknowledgementIt is my great pleasure to have been participated in the 21st Century COE Virtual Academy 2006(VA) of University of Yamanashi, in the past 4 months, I learned a lot from the program, and put them into practice. This work would not have been possible without the help of the BTOPMC-TEAM; I’d like to express my great thanks to them, for their impatience and great help in many aspects, thanks!References[1]: ZHAO R J. The Xinanjiang Model applied in China[J].Journal of Hydrology, 1992,135(4):371-381. [2]:http://coeinav.cec.yamanashi.ac.jp/inavi/scripts/inredir.dll.Author introduction ：Xin Penglei (1982-)，female ，born in Shan Dong province ，Graduate student ，Study on floodforecasting.。

新安江流域水文模型

第二章新安江流域水文模型60年代初，河海大学(原华东水利学院)水文系赵人授等开始研究蓄满产流模型，配合一定的汇流计算，将模型应用于水文预报和水文设计。

1973年，他们在对新安江水库做人库流量预报的工作中，把他们的经验归纳成一个完整的降雨径流流域模型——新安江模型。

模型可用于湿润地区和半湿润地区的湿润季节径流模拟和计算。

最初的新安江模型为两水源模型，只能模拟地表径流和地下径流。

80年代初期，模型研制者将萨克拉门托模型与水箱模型中，用线性水库函数划分水源的概念引入新安江模型，提出了三水源新安江模型，模型可以模拟地面径流、壤中流、地下径流。

1984至1986年，又提出了四水源新安江模型，可以模拟地面径流、壤中流、快速地下径流和慢速地下径流。

三水源新安江模型一般应用效果较好，但模拟地下水丰富地区的日径流过程精度不够理想。

在新安江三模型中增加慢速地下水结构就成为四水源新安江模型。

当流域面积较小时，新安江模型采用集总模型，当面积较大时，采用分块模型。

分块模型把流域分成许多块单元流域，对每个单元流域做产、汇计算，得到单元流域的出口流量过程。

再进行出口以下的河道洪水演算，求得流域出口的流量过程。

把每个单元流域的出流过程相加，就求得了流域出口的总出流过程。

划分单元流域的主要目的是处理降雨分布的不均匀性，因此单元流域应当大小适当，使得每块面积上的降雨分布比较均匀．并有一定数目的雨量站。

其次尽可能使单元流域与自然流域相一致，以便于分析与处理问题，并便于利用已有的小流域水文资料。

如果流域内有大中型水库，则水库以上的集水面积即应作为一个单元流域。

因为各单元流域的产汇、流计算方法基本相同，以下只讨论一个单元流域的情况。

2.1新安江两水源模型1．模型结构和参数新安江两水源模型的产流子模型采用蓄满产流模型，蒸发计算采用三层蒸发计算模型。

利用稳定下渗率FC将径流划分为地面径流和地下径流两种水源。

地面径流采用单位线汇流，地下径流采用一次线性水库汇流。

洪水调节设计(试算法和半图解法)模板 - 带试算C语言程序.

《洪水调节课程设计》任务书一、设计目的1.洪水调节目的：定量地找出入库洪水、下泄洪水、拦蓄洪水的库容、水库水位的变化、泄洪建筑物型式和尺寸间的关系，为确定水库的有关参数和泄洪建筑型式选择、尺寸确定提供依据；2.掌握列表试算法和半图解法的基本原理、方法、步骤及各自的特点；3.了解工程设计所需洪水调节计算要解决的课题；培养学生分析问题、解决问题的能力。

二、设计基本资料1.某水利枢纽工程以发电为主，兼有防洪、供水、养殖等综合效益，电站装机为5000KW，年发电量1372×104 kw·h，水库库容0.55亿m3。

挡水建筑物为混凝土面板坝，最大坝高84.80m。

溢洪道堰顶高程519.00m，采用2孔8m×6m（宽×高）的弧形门控制。

水库正常蓄水位525.00m。

电站发电引用流量为10 m3/s。

2.本工程采用2孔溢洪道泄洪。

在洪水期间洪水来临时，先用闸门控制下泄流量q并使其等于洪水来水量Q，使水库水位保持在防洪限制水位不变；当洪水来水量Q继续增大时，闸门逐渐打开；当闸门达到全开后，就不再用闸门控制，下泄流量q随水库水位z的升高而增大，流态为自由流态，情况与无闸门控制一样。

3.上游防洪限制水位524.8m（注：X=524.5+学号最后1位/10，即524.5m-525.4m），下游无防汛要求。

三、设计任务及步骤分别对设计洪水标准、校核洪水标准，按照上述拟定的泄洪建筑物的类型、尺寸和水库运用方式，分别采用列表试算法和半图解法推求水库下泄流量过程，以及相应的库容、水位变化过程。

具体步骤：1.根据工程规模和建筑物的等级，确定相应的洪水标准；2.用列表试算法进行调洪演算：①根据已知水库水位容积关系曲线V～Z和泄洪建筑物方案，用水力学公式求出下泄流量与库容关系曲线q～Z，并将V～Z，q～Z绘制在图上；②决定开始计算时刻和此时的q1、V1，然后列表试算，试算过程中，对每一时段的q2、V2进行试算；③ 将计算结果绘成曲线：Q ～t 、q ～t 在一张图上，Z ～t 曲线绘制在下方。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

新安江模型程序C++代码以下是类的声明：class XinanjiangModel{private:// FORCINGdouble *m_pP; // 降水数据double *m_pEm; // 水面蒸发数据//long m_nSteps; // 模型要运行的步长(一共m_nSteps步)long steps;// OUTPUTdouble *m_pR; // 流域内每一步长的产流量(径流深度)double *m_pRs; // 每一步长的地表径流深(毫米)double *m_pRi; // 每一步长的壤中流深（毫米）double *m_pRg; // 每一步长的地下径流深(毫米)double *m_pE; // 每一步长的蒸发(毫米)double *m_pQrs; // 流域出口地表径流量double *m_pQri; // 流域出口壤中流径流流量double *m_pQrg; // 流域出口地下径流量double *m_pQ; // 流域出口的总流量double m_U; // for 24h. U=A(km^2)/3.6/delta_t// SOILdouble *m_pW; // 流域内土壤湿度double *m_pWu; // 流域内上层土壤湿度double *m_pWl; // 流域内下层土壤适度double *m_pWd; // 流域内深层土壤湿度double m_Wum; // 流域内上层土壤蓄水容量double m_Wlm; // 流域内下层土壤蓄水容量double m_Wdm; // 流域内深层土壤蓄水容量，WDM=WM-WUM-WLM // EVAPORATIONdouble *m_pEu; // 上层土壤蒸发量（毫米）double *m_pEl; // 下层土壤蒸发量（毫米）double *m_pEd; // 深层土壤蒸发量（毫米）//runoffdouble *RF;// PARAMETERdouble m_Kc; // 流域蒸散发能力与实测蒸散发值的比double m_IM; // 不透水面积占全流域面积之比double m_B; // 蓄水容量曲线的方次，小流域（几平方公里）B0.1左右// 中等面积（平方公里以内）.2~0.3，较大面积.3~0.4 double m_WM; // 流域平均蓄水容量（毫米）(WM=WUM+WLM+WDM) double m_C; // 流域内深层土壤蒸发系数，江南湿润地区：0.15-0.2，//华北半湿润地区：.09-0.12double m_SM; //自由水蓄水容量double m_EX; //自由水蓄水容量～面积分布曲线指数double m_KG; //地下水日出流系数double m_KI; //壤中流日出流系数double m_CG; //地下水消退系数double m_CI; //壤中流消退系数double *m_UH; // 单元流域上地面径流的单位线double m_WMM; // 流域内最大蓄水容量double m_Area; // 流域面积int m_DeltaT; // 每一步长的小时数int m_PD; // 给定数据，用以判断是否时行河道汇流计算public:XinanjiangModel(void);~XinanjiangModel(void);// 初始化模型void InitModel(long nSteps, double Area,int DeltaT, int PD, char *ForcingFile);// 设置模型参数void SetParameters(double *Params);// 运行新安江模型void RunModel(void);// 保存模拟结果到文件void SaveResults(char *FileName);// 记录出流数据，用以作图分析void Runoff(char *runoff);private:// 进行汇流计算，将径流深度转换为流域出口的流量void Routing(void);};以下是类的定义#include"stdafx.h"#include"xinanjiangmodel.h"#include<iostream>#include<fstream>#include<iomanip>using namespace std;#include"math.h"#include"stdio.h"#include"conio.h"XinanjiangModel::XinanjiangModel(void){this->m_pP = NULL;this->m_pEm = NULL;this->m_pE = NULL;this->m_pEd = NULL;this->m_pEl = NULL;this->m_pEu = NULL;this->m_pW = NULL;this->m_pWd = NULL;this->m_pWl = NULL;this->m_pWu = NULL;this->m_pR = NULL;this->m_pRg = NULL;this->m_pRi = NULL;this->m_pRs = NULL;this->m_pQ = NULL;this->m_pQrg = NULL;this->m_pQri = NULL;this->m_pQrs = NULL;}XinanjiangModel::~XinanjiangModel(void) {delete[] this->m_pP;delete[] this->m_pEm;delete[] this->m_pE;delete[] this->m_pEd;delete[] this->m_pEl;delete[] this->m_pEu;delete[] this->m_pW;delete[] this->m_pWd;delete[] this->m_pWl;delete[] this->m_pWu;delete[] this->m_pR;delete[] this->m_pRg;delete[] this->m_pRi;delete[] this->m_pRs;delete[] this->m_pQ;delete[] this->m_pQrg;delete[] this->m_pQrs;delete[] this->m_pQri;}// 初始化模型void XinanjiangModel::InitModel(long nSteps, double Area, int DeltaT,int PD, char* ForcingFile){FILE * fp;int i;this->m_nSteps = nSteps;this->steps = this->m_nSteps + 18;// 驱动数据this->m_pP = new double[this->steps];this->m_pEm = new double[this->steps];// 模型输出，蒸散发项this->m_pE = new double[this->steps];this->m_pEd = new double[this->steps];this->m_pEl = new double[this->steps];this->m_pEu = new double[this->steps];// 模型输出，出流项，经过汇流的产流this->m_pQrg = new double[this->steps];this->m_pQrs = new double[this->steps];this->m_pQri = new double[this->steps];this->m_pQ = new double[this->steps];// 模型输出，产流项this->m_pR = new double[this->steps];this->m_pRg= new double[this->steps];this->m_pRi= new double[this->steps];this->m_pRs = new double[this->steps];// 模型状态量，土壤湿度this->m_pW = new double[this->steps];this->m_pWd = new double[this->steps];this->m_pWl = new double[this->steps];this->m_pWu = new double[this->steps];//runoff值this->RF = new double[this->steps];for(i = 0;i<this->steps;i++ ){// 驱动数据this->m_pP [i] = 0.00;this->m_pEm [i] = 0.00;// 模型输出，蒸散发项this->m_pE [i] = 0.00;this->m_pEd [i] = 0.00;this->m_pEl [i] = 0.00;this->m_pEu [i] = 0.00;// 模型输出，出流项，经过汇流的产流this->m_pQrg[i] = 0.00;this->m_pQrs[i] = 0.00;this->m_pQri[i] = 0.00;this->m_pQ[i] = 0.00;// 模型输出，产流项this->m_pR [i] = 0.00;this->m_pRg [i] = 0.00;this->m_pRi [i] = 0.00;this->m_pRs [i] = 0.00;// 模型状态量，土壤湿度this->m_pW [i] = 0.00;this->m_pWd[i] = 0.00;this->m_pWl[i] = 0.00;this->m_pWu[i] = 0.00;}this->m_Area = Area;this->m_DeltaT = DeltaT;this->m_PD = PD;this->m_U = this->m_Area/(3.6 * this->m_DeltaT);// Forcing文件格式：第一列：降水（单位毫米）空格第二列水面蒸发（毫米）if((fp = fopen(ForcingFile,"r")) == NULL){printf("Can not open forcing file!\n");return; }for(i = 0;i<this->m_nSteps;i++ ){ fscanf(fp,"%lf%lf",&(this->m_pP[i]),&(this->m_pEm[i])); }fclose(fp);}// 设置模型参数void XinanjiangModel::SetParameters(double* Params){this->m_Kc = Params[0]; // (1) 流域蒸散发能力与实测水面蒸发之比this->m_IM = Params[1]; // (2) 流域不透水面积占全流域面积之比this->m_B = Params[2]; // (3) 蓄水容量曲线的方次this->m_Wum = Params[3]; // (4) 上层蓄水容量this->m_Wlm = Params[4]; // (5) 下层蓄水容量this->m_Wdm = Params[5]; // (6) 深层蓄水容量this->m_C = Params[6]; // (7) 深层蒸散发系数this->m_SM = Params[7]; // (8)自由水蓄水容量this->m_EX = Params[8]; // (9)自由水蓄水容量～面积分布曲线指数this->m_KG = Params[9]; // (10)地下水日出流系数this->m_KI = Params[10]; // (11)壤中流日出流系数this->m_CG = Params[11]; // (12)地下水消退系数this->m_CI = Params[12]; // (13)壤中流消退系数this->m_WM = this->m_Wum + this->m_Wlm + this->m_Wdm;this->m_WMM = this->m_WM * (1.0 + this->m_B)/(1.0 - this->m_IM); }// 运行新安江模型void XinanjiangModel::RunModel(void){long i;// 模型的状态变量double PE; // > 0 时为净雨量;< 0 为蒸发不足量（mm）double Ep; //m_Kc * m_pEm[i]double P;double R; // 产流深度，包括地表径流、壤中流和地下径流（mm）double RB; // 不透水面上产生的径流深度（mm）double RG; // 地下径流深度（mm）double RI; // 壤中流深度（mm）double RS; // 地表径流深（mm）double A; //土壤湿度为W时土壤含水量折算成的径流深度（mm）double E = 0.0; // 蒸散发(mm)double EU = 0.0; // 上层土壤蒸散发量（mm）double EL = 0.0; // 下层土壤蒸散发量（mm）double ED =0.0; // 深层土壤蒸散发量（mm）double S;double FRo;double FR;double MS;double AU;double WU = 5.0; // 流域内上层土壤湿度double WL = 55.0; // 流域内下层土壤适度double WD = 40.0; // 流域内深层土壤湿度double W = 100.0;double So = 5.0;MS = m_SM * (1 + m_EX);FRo = 1 - pow((1 - So/MS),m_EX);for(i = 0;i<this->m_nSteps;i++ ){// ——————蒸散发计算————————————//RB = m_pP[i] * m_IM; // RB是降在不透水面的降雨量P = m_pP[i] * (1 - m_IM);Ep = m_Kc * m_pEm[i];if ((WU + P)>= Ep){EU = Ep; EL = 0; ED = 0; }else if((WU + P)<Ep){EU = WU + P;if(WL>= (m_C * m_Wlm)){ EL = (Ep - EU) * WL/m_Wlm; ED = 0; }else if ((m_C * (Ep - EU))<= WL&&WL<(m_C * m_Wlm)){ EL = m_C * (Ep - EU); ED = 0; }else if (WL<m_C * (Ep - EU)){ EL = WL; ED = m_C * (Ep - EU) - EL; }}E = EU + EL + ED;PE = P - E;/* ———————蒸散发计算结束——————————— */ //——————子流域产流量计算————————————// if(PE<= 0){ R = 0.00; W = W + PE; }else{A = m_WMM * (1 - pow( (1.0 - W/m_WM), 1.0/(1 + m_B) ) );// 土壤湿度折算净雨量+降水后蒸发剩余雨量<流域内最大含水容量if((A + PE)<this->m_WMM){// 流域内的产流深度计算R = PE /* 降水蒸发后的剩余量*/+ W /* 流域内土壤湿度*/+ m_WM * pow((1 - (PE + A)/m_WMM),(1 + m_B))- m_WM /* 减去流域平均蓄水容量（m_WM:参数）*/+ RB; /* 不透水面上产生的径流*/}// 土壤湿度折算净雨量+降水后蒸发剩余雨量<流域内最大含水容量else{// 流域内的产流深度计算R = PE /* 降水蒸发后的剩余量+ W /* 流域内土壤湿度*/- m_WM /* 减去流域平均蓄水容量*/+ RB; /* 不透水面上产生的径流*/}}//三层蓄水量的计算: WU, WL, WDif(WU + P - EU - R <= m_Wum){ WU = WU + P - EU - R; WL = WL - EL; WD = WD – ED; }else{WU = m_Wum;if(WL - EL + ( WU + P - EU - R - m_Wum ) <= m_Wlm ){WL = WL – EL + ( WU + P - EU - R - m_Wum );WD = WD - ED;}else{WL = m_Wlm;if(WD - ED + WL - EL + ( WU + P - EU - R - m_Wum )- m_Wlm <= m_Wdm )WD = WD - ED + WL - EL+ (WU + P - EU - R - m_Wum ) - m_Wlm;elseWD = m_Wdm;}}W = WU + WL + WD;////三水源划分汇流计算if(PE>0){FR = (R - RB) / PE;AU = MS * (1 - pow((1 - So * FRo/FR/m_SM),1/(1 + m_EX)));if(PE + AU<MS)RS = FR * (PE + So * FRo/FR - m_SM + m_SM * pow((1 - (PE + AU) / MS),m_EX + 1));else if(PE + AU>= MS)RS = FR * ( PE + So * Fro / FR - m_SM );S = So * Fro / FR + ( R – RS ) / FR;RI = m_KI * S * FR;RG = m_KG * S * FR;RS += RB;R = RS + RI + RG;So = S * ( 1 - m_KI - m_KG );FRo = FR;}else{S = So;FR = 1 - pow((1 – S / MS) , m_EX );RI = 0.00;RG = 0.00;So = S * ( 1 - m_KI - m_KG );RS = RB;R = RS + RI + RG;FRo = FR;}////三水源划分计算结束/* 以下部分是状态量：总蒸发量、上、下和深层土壤的蒸发的保存*//* 1 */this->m_pE[i] = E; // 当前步长的蒸发（模型重要输出）/* 2 */this->m_pEu[i] = EU; // 当前步长上层土壤蒸发/* 3 */this->m_pEl[i] = EL; // 当前步长下层土壤蒸发/* 4 */this->m_pEd[i] = ED; // 当前步长深层土壤蒸发/* 5 */this->m_pW[i] = W; // 当前步长流域平均土壤含水量/* 6 */this->m_pWu[i] = WU; // 当前步长流域上层土壤含水量/* 7 */this->m_pWl[i] = WL; // 当前步长流域下层土壤含水量/* 8 */this->m_pWd[i] = WD; // 当前步长流域深层土壤含水量/* 9 */this->m_pRg[i] = RG; // 当前步长流域地下径流深度/* 10 */this->m_pRi[i] = RI; // 当前步长流域壤中流深度/* 11 */this->m_pRs[i] = RS; // 当前步长流域地表径流径流深度/* 12 */this->m_pR[i] = R; // 当前步长的总产流径流深度}this->Routing();}// 保存模拟结果到文件void XinanjiangModel::SaveResults(char* FileName){int i;FILE * fp;if((fp = fopen(FileName,"w")) == NULL){printf("Can not create output file!\n");return;}fprintf(fp," - -- -- ---- - - -- -- ---- - ------ --- -- ------ - - -- ---- ----- -- - - ------- -- -- -\n");fprintf(fp," E(mm) EU(mm) EL(mm) ED(mm) W(mm) WU(mm) WL(mm) WD(mm) R(mm) RS(mm) RI(mm) RG(mm) Q(m3/d) QS(m3/d) QI(m3/d) QG(m3/d)\n");fprintf(fp," - -- -- -- - - --- -- -- - - -- ----- -- - - -- -- -- - - -- -- --------- - - -- --------- ---\n");for(i = 0;i<this->steps;i++ ){ fprintf(fp,"%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf\n",this->m_pE[i],this->m_pEu[i],this->m_pEl[i],this->m_pEd[i],this->m_pW[i],this->m_pWu[i],this->m_pWl[i],this->m_pWd[i],this->m_pR[i],this->m_pRs[i],this->m_pRi[i],this->m_pRg[i],this->m_pQ[i],this->m_pQrs[i],this->m_pQri[i],this->m_pQrg[i]);}fclose(fp);}// 进行汇流计算，将径流深度转换为流域出口的流量void XinanjiangModel::Routing(void){///////////// 地面径流汇流计算：单位线法///////////////////////int i,j;double B[10000] = {0.00};if (this->m_PD == 1){double UH[] ={3.71,12.99,38.96,94.63,131.74,154.00,166.99,176.27,178.12,172.55,146.58, 90.91,53.80, 31.54,18.55, 9.27, 3.71,0.00};for(i = 0;i<this->m_nSteps;i++ ){for(j = 0;j<18;j++ ){ B[i + j] += this->m_pRs[i] * UH[j]/10.0; }}}else{double UH[] ={ 7.18,23.38,63.20,143.10,221.75,365.18,447.40,491.29,506.93,504.82,468.46,388.56,309.91,166.49,84.26,40.37,17.56,3.46};for(i = 0;i<this->m_nSteps;i++ ){for(j = 0;j<18;j++ ){ B[i + j] += this->m_pRs[i] * UH[j]/10.0; }}}for(i = 0;i<this->steps;i++ )this->m_pQrs[i] = B[i];///// 壤中流汇流计算:线性水库for(i = 1;i<this->steps;i++ ) {this->m_pQri[i] = this->m_CI * this->m_pQri[i - 1] + (1.0 - this->m_CI) * this->m_pRi[i] * this->m_U; }///// 地下径流汇流计算：线性水库for(i = 1;i<this->steps;i++ ) {this->m_pQrg[i] = this->m_pQrg[i - 1] * this->m_CG + this->m_pRg[i] * (1.0 - this->m_CG) * this->m_U; }//////单元面积总入流计算for(i = 0;i<this->steps;i++ ){ this->m_pQ[i] = this->m_pQrs[i] + this->m_pQri[i] + this->m_pQrg[i]; }}void XinanjiangModel::Runoff(char * runoff){int i;ofstream outfile;outfile.open(runoff);if(outfile.is_open()){for(i = 0;i<this->steps;i++ ){ outfile<<setprecision(3)<<setiosflags(ios::fixed)<<this->m_pQ[i]<<endl; } }outfile.close();}以下是main()函数语句int _tmain(int argc, _TCHAR * argv[]){ long nSteps = 942;int DeltaT = 24;double Area1 = 1603;XinanjiangModel Model1;Model1.InitModel(nSteps, Area1,DeltaT,1,"LFForcingfile.txt");//模型参数/*Kc,IM,B,m_Wum,Wlm,Wdm,C,SM,EX,KG,KI,CG,CIdouble Params1[] = {0.50,0.01,0.30,10,60,40,0.18,32,1.2,0.075,0.072,0.94,0.7};Model1.SetParameters(Params1);Model1.RunModel();Model1.SaveResults("来凤站日模型计算结果.txt");Model1.Runoff("LF_Q.txt");Model1.~XinanjiangModel();double Area2 = 2991;XinanjiangModel Model2;Model2.InitModel(nSteps, Area2,DeltaT,0,"YCForcingfile.txt");//模型参数/*Kc,IM,B,m_Wum,Wlm,Wdm,C,SM,EX,KG,KI,CG,CIdouble Params2[] = {0.75,0.01,0.32,10,60,40,0.18,27,1.2,0.065,0.067,0.96,0.8}; Model2.SetParameters(Params2);Model2.RunModel();Model2.SaveResults("file.txt");Model2.Runoff("YC_Q.txt");Model2.~XinanjiangModel();FILE *fp1,*fp2;double Q1[1000],Q2[1000],Q[1000]={0.00};ofstream outfile;outfile.open("Q.txt");if((fp1=fopen("LF_Q.txt","r"))==NULL){ printf("Can not open the file!\n"); return 0;}if((fp2=fopen("YC_Q.txt","r"))==NULL){ printf("Can not open the file!\n"); return 0; }if(outfile.is_open()){for(int i=0;i<960;i++){fscanf(fp1,"%lf",&Q1[i]);fscanf(fp2,"%lf",&Q2[i]);Q[i]=Q1[i]+Q2[i];outfile<<setprecision(3)<<setiosflags(ios::fixed)<<Q[i]<<endl;}fclose(fp1);fclose(fp2);}outfile.close() ;return 0;}。

新安江模型程序C 代码