新安江模型程序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 曲线绘制在下方。

6 新安江模型汇总

2018/11/1
9
（2）超蓄产流模型的结构 a）点模型以含气层缺水量为控制条件，就流域中某点而言：
蓄满前: P E WW2 WW1 蓄满后: P E R 式中: P : 时段降雨量 E : 时段蒸散发量 R : 时段产流量
(6 - 1)
WW1 , WW2 : 时段初末的土壤含水量
21.11
30.51 24.27 84.62
0.0427
0.0632 0.0495 0.2091
0.50
0 7.46 17.61
0.50
0 2.22 5.84
0
0 5.23 11.76
29.57
23.68 76.57 113.20
23
24
20.27
-2.79
138.85
156.00
0.4844
1.0000
c）流域产流计算 P－E>0时，产流，否则不产流，产流时：
P E A WWMM 时 : R P E ( WM W ) P E A WWMM 时 : R P - E - (WM - W) - WM1 - (P - E A)/WWMM

1 B

13
2018/11/1
10
式(6 1)中 , R RS RG , 即 RG FC RS R - RG P - E - FC FC : 时段稳定下渗量
b）流域蓄水容量曲线（超蓄产流模型的核心）
(6 - 2)
WWM：流域蓄水容量 WWMM:流域最大蓄水容量
WM：流域平均蓄水容量
I
E
2018/11/1
18
2、用试算法求fc f RS i R i i f c t i F

新安江模型参数的分析

一、模型的结构与参数三水源新安江模型的流程图如图1所示。

图1 三水源新安江模型流程图图1 中输入为实测雨量P ，实测水面蒸发EM ；输出为流域出口流量Q ，流域蒸散发E 。

方框内是状态变量，方框外是参数变量。

模型结构及计算方法可分为以下四大部分。

1．蒸散发计算用三个土层的模型，其参数为上层张力水容量UM ，下层张力水容量LM ，深层蒸散发系数C ，蒸散发折算系数K ，所用公式如下：当上层张力水蓄量足够时，上层蒸散发EU 为EM E EU ⨯=当上层已干，而下层蓄量足够时，下层蒸散发EL 为LM WL EM K EL /⨯⨯=当下层蓄量亦不足，要触及深层时，蒸散发ED 为EM K C ED ⨯⨯=2．产流量计算据蓄满产流概念，参数为包气带张力水容量WM ，张力水蓄水容量曲线的方次B ，不透水面积的比值IM ，所用公式为)1/()1(IM B WM WM -+⨯=))/1(1()1/(1B WM W MM A +--=当0≤⨯-EM K P ，则R=0不然，则当MM A EM K P <+⨯-，B MM A EM K P WM W WM EM K P R ++⨯--⨯++-⨯-=1)/)(1(不然，则W WM EM K P R +-⨯-=式中 R ——产流量；MM ——流域最大点蓄水容量。

3．分水源计算分三种水源，即地面径流RS 、地下径流RG 和壤中流RI 。

参数为表层土自由水蓄水容量SM ，表层自由水蓄水容量曲线的方次EX ，表层自由水蓄量对地下水的出流系数KG 及对壤中流的出流系数KI ，所用公式为SM EX MS ⨯+=)1())/1(1()1/(1EX SM S MS AU +--⨯=)/())((EM K P EM K P IM R FR ⨯-⨯-⨯-=FR KG S RG ⨯⨯=FR KI S RI ⨯⨯=当 0,0=≤⨯-RS EM K P不然，当MS AU EM K P <+⨯-，则FR MS AU EM K P SM S SM EM K P RS EX ⨯+⨯--⨯++-⨯-=+))/)(1((1 当MS AU EM K P ≥+⨯-，则FR SM S EM K P RS ⨯-+⨯-=)(4．汇流计算地下径流用线性水库模拟，其消退系数为CG ，出流进入河网。

新安江流域水文模型

－2002.12－
世界上第一个流域水文模型－Stanford模型出现在20世纪60年代，目前全世界已提出数以百计的流域水文模型。主要包括由美国天气局V. T. Sitten提出的API模型、N. H. Crawford和R. K. Linsley提出的斯坦福模型以及R. J. C. Bernash等提出的萨克拉门托模型，日本国立防灾科学研究中心菅原正已教授提出的水箱模型，丹麦技术大学提出的 NAM模型，以及原华东水利学院赵人俊教授提出的新安江模型。这些概念性水文模型对流域的降雨径流过程进行了较为细致的模拟。由于这些模型具有较好的结构形式和良好的模拟预报精度，因此在洪水实时预报中得到广泛地应用。本文主要介绍国内应用最为广泛的新安江三水源模型。
－2002.12－
新安江模型的结构
蒸散发计算原理
各层蒸散发的计算原则是，上层按蒸散发能力蒸发，上层含水量蒸发量不够蒸发时，剩余蒸散发能力从下层蒸发，下层蒸发与蒸散发能力及下层含水量成正比，与下层蓄水容量成反比。要求计算的下层蒸发量与剩余蒸散发能力之比不小于深层蒸散发系数。否则，不足部分由下层含水量补给,当下层水量不够补给时，用深层含水量补充。(上层以蒸散发能力蒸发,直到上层水分耗尽, 才蒸发下层;下层土壤蒸散发量与剩余蒸散发能力(流域蒸散发能力与上层蒸散发量之差)及下层土壤实际含水量成正比。)
－2002.12－
流域蓄水容量曲线
新安江模型的结构－产流量计算
－2002.12－
新安江模型的结构－产流量计算
据此可求得流域平均蓄水容量为
WM

Wmm (1
0
f F
)dWm

Wmm B 1
与流域初始平均蓄水量W0 相应的纵坐标A(
)为 A

基于SCE-UA算法的英那河水库流域新安江模型参数优化及应用

科技论坛
・３・３
基于ＳＥＵＡ算法的英那河水库流域新安江模型Ｃ — 参数优化及应用
Байду номын сангаас、
李博
辛江
王志刚
（宁省水文水资源勘测局大连分局，宁大连１６２）辽辽１０３
摘要：以大连市英那河水库为例，用ＳＥｕＡ算法对新安江模型参数进行优化，应Ｃ— 确定出适合英那河水库的场次洪水的模型参数，实现新安江模型在英那河水库洪水预报中的准确应用。关键词：新安江模型；Ｃ — ＳＥｕＡ算法；参数优化概念性模型的参数是有一定的物理意义的，因此在原则上应该按其物理意义来定量，或者实测，或者实验，或者参考类似的情况。但实际上并不可能做得到，因为许多量没有观测值，者观测值代或表性不好。因此，常用的办法是先以实测值或类似经验定好的参数作为初值，然后用模型计算出产流过程，与实测过程进行比较，再作优化调试，以误差最小为原则，确定参数的最优值。然而模型计算最重要的部分就是确定参数。对于参数的优选直接决定着模型计算的序号成败。参数优选是水文模型研究和应用的一个重要方面，在水文模型发展初期一般都采用手动率定参数，由于模型参数较多，但手动率定不但费时费力，而且还可能得不到理想的结果ｌｌｌ。近年来，随着计算机技术的发展，自动优化方法被广泛应用到模型参数率定中。就新安江模型的参数自动优化方法来说，目前国内使用地比较普遍的自动优化方法主要有粒子群算法，遗传算法以及ＳＥｕＣ — Ａ算法日。本研究采用ＳＥＵＣ — Ａ算法对新安江模型参数进行优选。表１三水源新安江模型产流参数初值范围

c语言程序设计水利水电

C语言程序设计水利水电概述水利水电工程是指利用水资源进行电力生产以及水利枢纽工程和渠道建设等水利工程。

C语言是一种通用的高级编程语言，具有简洁、高效、可移植等特点。

本文档将介绍如何运用C语言进行水利水电工程的程序设计，以实现相关功能。

目录1.水利水电工程简介2.C语言程序设计基础知识1.数据类型2.控制流程3.函数4.数组和指针5.结构体3.水文数据处理程序设计1.数据输入与输出2.数据处理与计算3.数据存储与读取4.水电站模拟程序设计1.水库水位与流量计算2.水电站发电量预测3.系统运行状态检测与控制5.水利工程可视化设计1.图形界面设计2.数据可视化展示3.用户交互与操作6.水利水电工程实例分析1.坝体稳定性分析2.溢洪道流量计算3.水电站效益评估7.总结与展望1.水利水电工程简介水利水电工程是人类利用水资源进行发电和水资源利用的重要项目。

包括水库、大坝、水电站、引水渠等建设和管理。

C语言程序设计作为一种强大的编程语言，可以应用于水利水电工程的自动化控制和数据处理等方面。

2. C语言程序设计基础知识在进行水利水电工程的程序设计之前，我们需要掌握一些C语言的基础知识。

这些知识包括数据类型、控制流程、函数、数组和指针以及结构体等。

2.1数据类型C语言中有多种数据类型，包括整型、浮点型、字符型等。

我们需要了解这些数据类型的用途、取值范围以及相应的格式符号。

2.2控制流程控制流程是程序执行的路线控制，包括顺序结构、选择结构和循环结构。

我们需要掌握如何使用i f语句、swi t ch语句和f or循环、wh il e 循环等进行流程控制。

2.3函数函数是C语言程序的基本组成单元，可以方便地对代码进行模块化设计。

我们需要学会如何定义函数、调用函数以及传递参数和返回值。

2.4数组和指针数组和指针是C语言中常用的数据结构，可以用于存储和操作一组相关的数据。

我们需要了解如何定义和使用数组以及指针的作用和使用方法。

新安江模型介绍

水文第七组成员：刘俊、冯远、曹胜荣、杨春智
数据输入：刘俊、冯远、曹胜荣、杨春智
程序编写：杨春智
总结报告：杨春智、曹胜荣、冯远、刘俊
小组成员：杨春智、曹胜荣、冯远、刘俊
2010年11月08日
{
WU[i+1]=0;
WL[i+1]=WL[i]-EL[i];
WD[i+1]=WD[i];
}
else if(WD[i]-ED[i]>0)
{
WU[i+1]=0;
WL[i+1]=0;
WD[i+1]=WD[i]-ED[i];
}
}
else
{
WU[i+1]=WU[i]+PE[i]-R[i];
WL[i+1]=WL[i];
7自由水蓄水容量—面积分布曲线指数EX
EX反映流域自由水蓄水分布的不均匀程度，一般EX=1.0~1.5
8自由水蓄水库对地下水和壤中流的日出流系数KG+KI
KG的大小反映基岩和深层土壤的渗透性，KI的大小反映表层土的渗透性。
9地下水消退系数CG
CG可根据枯季地下径流退水规律来推求。若以日作为计算时段长，CG=0.950~0.998
a[i]=WMM*(1.0-pow((1.0-W[i]/WM),(1.0/(1.0+B))));
//时段产流量的计算
if(PE[i]>0)
{if(a[i]+PE[i]<=WMM)
R[i]=PE[i]+W[i]-WM+WM*pow((1.0-(PE[i]+a[i])/WMM),(B+1.0));
else
4不透水面积占全流域面积的比例IM

全国河流代码表之欧阳引擎创编

河流代码表
黑龙江流域海拉尔河库都尔河特尼河
伊敏河辉河鄂依那河
威特很高勒河
锡尼很高勒河
宁高勒河
免渡河
乌尼日河
克鲁伦河
呼仑湖
额尔古纳河
莫尔格勒高勒河
乌尔逊河
贝尔湖
哈拉哈河
根河
图里河
额根河
得耳布干河
莫尔道嘎河
古纳河
激流河
安格林河
敖鲁古雅河
阿龙山河
金河
阿巴河
恩和哈达河
黑龙江
北极村河
额木尔河(克波河、大林河)
古莲河
老槽河
二龙河
大林河
盘古河
西尔根气河
小西尔根气河
呼玛河
阿吉羊河
卡马兰河
亚里河
呼玛尔河
欧阳引擎创编 2021.01.01
欧阳引擎创编 2021.01.01。

新安江水文模型简介

《流域水文模拟》结课报告新安江模型的原理、结构及应用、发展历程The principle, structure, application anddevelopment process of Xin anjiang Model作者姓名：孔旭学科、专业：水文学及水资源学号：指导教师：王国利完成日期： 2016年8月30日大连理工大学Dalian University of Technology摘要新安江模型是河海大学提出的一个概念性降雨径流模型，具有原创性，是我国为数不多的被国际上广泛认可的水文模型。

新安江水文模型在我国湿润与半湿润地区广为应用，取得了良好的效果。

经过近50年的发展，新安江模型已经从最初的专门从事水库入库洪水预报的单一功能模型发展为适合用于水文预报、水资源管理、水土资源评价、面源污染预测、气候变化和人类活动影响研究的多功能的水文模型；其部分参数已从靠经验率定发展为可以进行物理推求。

总之，新安江模型是一个不断发展的模型体系。

本文主要由三部分构成。

第一部分为新安江模型简介，回顾了新安江模型产生的历史背景和发展历程，介绍了新安江模型的基本原理和结构体系；第二部分讲述了新安江模型参数的物理意义及其率定；第三部分为新安江水文模型在英那河流域防洪规划编制当中的应用。

关键词：水文模型；新安江模型；洪水预报The principle, structure, application and development processof Xin anjiang ModelAbstractXin anjiang Model originally proposed by Hehai University is a conceptual rainfall runoff model and is also one of the few widely recognized international hydrological model in China. Xin anjiang hydrological model was widely used in our humid and semi-humid areas, and achieved good results.After nearly 50 years study, Xin anjiang model has been developed from the single-function of reservoir flood forecasting into multi-purpose model including hydrological forecasting, water resources management, water and soil resources evaluation, non-point source pollution prediction, climate change and human activities versatile hydrological model studies. And part of its parameters can be acquired through physical calculation instead of experience. In short, Xin anjiang model is an evolving model system.This paper consists of three parts. The first part is about the brief introduction of Xin anjiang model, which recalls the historical background and the development, as well as introduces the basic principles and architecture; the second part describes the physical meaning of Xin anjiang model parameters and calibration; the third part is about the application of Xin anjiang model in Ying Na River Basin flood control planning.Key Words: hydrological model; Xin anjiang model; Flood forecasting目录摘要.............................................. 错误!未定义书签。

马斯京根法及新安江模型

斯京根连续演算
x－流量比重因子，与河道调蓄能力呈反比流量比重因子，
马斯京根法一般无预见期， =0时马斯京根法一般无预见期，仅当C0=0时（即 Δt=2Kx）,有1个时段预见期
新安江模型
新安江三层蒸发三水源模型简介
流域蒸发：流域蒸发：三层模式
土壤含水量：土壤含水量：递推公式流域产流：流域产流：分水源：分水源：坡面汇流：坡面汇流：河网汇流：河网汇流：河道汇流：河道汇流：蓄满产流模式自由水水库线性水库无因次时段单位线马斯京根法
WUt +1 =
(WUt + P − EUt − Rt ) t WUt +1 ≤ WUm
WUm
WUm
WUt +1 > WUm
WLm
WLt+1= Wt+1 – WUt+1
蓄满产流公式
P + A 1+b P +W0 −Wm +Wm(1− e ' ) e Wmm R= R = P +W −Wm e 0
TRss
河网汇流无因次时段单位线
u
t
河道汇流马斯京根法
Q下， 2 ＝Ｃ０ Q上， 2 ＋Ｃ１ Q上， 1 ＋Ｃ２ Q下， 1
C0 =
C1 =
0 . 5 ∆ t − Kx K − Kx + 0 . 5 ∆ t
0 . 5 ∆ t + Kx K − Kx + 0 . 5 ∆ t
C2 =
E水
P
S0 、Fr0
新安江三水源模型运行框图
β
Ep
Wm、WUm、WLm、C

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 代码