实验一：径流的随机分析法（古尔德法）.doc

水资源学课程设计 04级水文水资源1班 王静 200431580185 实验一：径流的随机分析法（古尔德法）（一）、古尔德法的原理古尔德法又称直接求总库容的概率演算法，它以实测年、月径流资料为基础，假设入库年径流过程为独立随机序列，以年初蓄水状态为条件，对各年进行年内各月的调节演算求年末的蓄水状态，然后再根据实际演算结果，统计出蓄水状态的年转移概率，并沿用马氏过程的原理分析水库蓄水状态的年转移概率矩阵，并推求出水库稳定的蓄水概率分布。古尔德法的特点是它不仅考虑了水库对径流的年际调节，而且同时考虑了径流的年内调节。（二）、古尔德法步骤1、将水文年径流资料转为水利年径流资料：结果如表“径流资料” 2、状态离散： 将水库的兴利库容Vc从库空（0）到库满（Vn），其间等分为（m2）份，（m由精度等而定，本次取为10），每份为DVVn/（m2），各状态序号i所对应的状态值A（i）及状态的数值域x（i）按下式计算： il时 A(i)0 x(1)0 1im时 A（i）=DVi(3/2) (i-2)DVx(i)(i-1)DV im时 A（i）Vn x(m)Vn 3、推求转移概率矩阵和条件破坏概率以一个指定的年初蓄水状态A（i）为出发点，根据某年实测径流资料和供水资料，采用简化运行策略进行年内调节计算，求出水库各月的蓄水过程及其年末的蓄水状态。记K表示年份序号，L表示月份的序号，S（K，L）表示水库的第K年L月份的月初蓄水量，X（K，L）表示第K年L月份的月径流，d（L）表示第L月份的调节流量，则初始条件为S（K,1）A（i）逐月径流调节计算时，实际上是执行下式的水量平衡演算，即 S（K，L+1）=S（K，L）+X（K，L）d（L）演算中规定：水库的蓄水量应满足下式的约束，即 0S（K，L）Vn 当S（K，L+1）Vn时，应取S（K，L+1）=Vn ，其多余水量为S（K，L） x（ K，L）d（L）Vn，作为弃水泄放。当S（K，L+1）时，应取S（K，L+1）=，水库只能按S（K，L） x（ K，L）供水，由于S（K，L） x（ K，L）d（L），因此这个月将出现正常供水的破坏，规定若一年中任何一个月出现正常供水的不足则该年计为不能保证正常供水的年份。即这一年供水破坏。年末蓄水状态S（K，13）视其所在的状态数值域判断其状态序号j。 按上述演算方法，可对实测各水利年每年都以年初S（K，L）=A（i）为初始蓄水状态进行各年的演算，求出各年的水库蓄水过程及其相应的年末蓄水状态S（K，13）=A（j）。 然后统计出各年初蓄水状态序号i下各年末蓄水状态序号j的出现年数，并据此计算出蓄水状态的年转移概率qij，以及以年初蓄水状态为条件的供水破坏概率FY、FM。蓄水状态的年转移概率可按下式计算。 条件供水破坏概率按下式计算： 对于年初蓄水状态X（i）,il，2，3，10，逐一进行相同的分析，最后可求得蓄水状态的年转移概率矩阵。 Qqi,j i,jl，2，10及各状态点的条件供水破坏概率： FY=fy(1)fy(10) FMfm(1)fm(10) 4、推求水库稳定的蓄水概率分布（用迭代法进行演算）假定第一年初的蓄水状态的概率为： 1=1 0 0 0 0 0 0 0 0 0利用蓄水概率分布演算式t1=t Q进行逐年迭代演算，求出第一年末的蓄水概率分布2lQ，以第一年的年末的蓄水状态的概率分布2作为第二年的初蓄，求第二年末的蓄水概率分布3，直到此概率分布收效于一个稳定的蓄水概率分布，此时t1t。 5、水库正常供水保证率计算 前面在推求水库蓄水状态的年转移概率的同时，也求出了条件的洪水破坏概率FY、FM，这些条件的供水破坏概率只与年初蓄水状态有关。而水库的稳定的蓄水状态的概率分布代表水库正常运行条件下的年初蓄水概率分布，根据概率论原理可知，年初蓄水出现状态X（i）与该年出现正常供水破坏这二个事件同时出现的概率等于二者各自的出现概率的乘积，故可按全概率公式求水库长期运行的正常供水破坏概率。按相对破坏年数表示，供水破坏概率的计算式为： 按相对破坏月数表示，供水破坏概率： 分别求出相对年数表示的正常供水保证率PY PY=1HY 和相对月数表示的正常供水保证率PM PM=1HM（三）、简化运行策略水库的运行策略是指水库在已知来水的情况下，如何放水的问题，而简化运行策略则是运行策略的一种概化策略，它是一种将水库的入流量作为一个确定性的径流过程的简化的计算方法，适用于水资源系统规划阶段的水库模拟运行（或长系列操作），它与径流调节中的等流量调节时历法相似。图3.1为水库简化运行策略的示意图。图中abcd为水库供水指示线，Vn为有效库容，M为时段（月）额定供水量。MQpT，Qp为调节流量。S(K)+X(K)dVnbW(K)Mca 图3.l 水库简化运行策略示意图模拟实验方法 将时段入流作为已知，死库容取0。 设时段K入库水量为X（K），时段初水库蓄水量为S（K），时段末水库蓄水量为S（Kl），水库的实际时段供水量为W（K）。 以如下规则按简化运行策略进行操作演算。 第一条规则：S（K）X（K）M，W（K）S（K）X（K），S（K+1）0； 第二条规则：M S（K）X（K） MVn，W（K）= M，S（K1）S（K）+ X（K）M； 第三条规则：S（K）十X（K）M+Vn，W（K）S（K）X（K）Vn，S（K1）Vn； 第一条规则与指示线ab段相对应，即当时段初蓄水量S（K）与该时段入库水量X( K）之和小于该时段要求的额定供水量M时，水库只能按时段初蓄水量与时段入库水量之和进行供水，控制时段末蓄水量不致于消落至死水位（死库容）以下，因此ab段反映水库限制供水执行，执行结果是该时段不能满足正常供水要求（即正常供水遭到破坏）。 第二规则与指示线bc段相对应，即当时段初蓄水量S(K)与该时段入库水量X(K)之和大于该时段要求的额定供水量M时，水库在满足正常供水量的同时，尚有多余水量可充蓄水库。 第三条规则与指示线cd段相对应，即当时段初蓄水量S(K）与该时段入库水量X(K)之和除满足正常供水量M和充蓄有效库容Vn之外，仍有余水，可加大供水，但控制水库供水量不超过水库引水设备的最大过水能力，超过则通过泄洪建筑物泄放（弃水），使水库水位维持在正常蓄水位上。 (四)、程序/* x(i,j)-水库入流（年，月径流量） A(i)-蓄水状态变量 S(i,j)-月初蓄水量 S1(i)-年末蓄水量 Q(i,j)-年转移概率矩阵 IND(j)-年末出现小于等于某一状态的总次数 PI(i)-水库稳定蓄水概率分布 FY(i)-条件供水破坏概率（相对破坏年数） FM(i)-条件供水破坏概率（相对破坏月数） QP-调节流量(已转化为水量） Vn-有效库容 M-状态总数 K1-计算总年教 PY-年供水保证率 PM-月供水保证率*/#include#include#define M 60#define K1 49#define m 18void main()/*状态离散*/int i,j,k,l,e;float ZM,VM,Am;float Z0=162.0,Zn=207.5;float V0,Vn,DV;ifstream ZV(H:水位库容关系.txt);for(i=0;iZiVi; ZV.close(); for(i=0;i=Zi&Z0=Zi&ZnZi+1)Vn=Vi+(Vi+1-Vi)*(Zn-Zi)/(Zi+1-Zi);DV=Vn/(m-2);A0=0;Am-1=Vn;for(i=1;im-1;i+)Ai=DV*i-DV/2;/*年末蓄水状态及推求条件破坏概率*/float QP=7.8;float FYm,FMm,INDm;float SK113,xK112,S1K1,Qmm;int IY,IM,INF; ifstream Qin(H:径流资料.txt);for(i=0;iK1;i+)for(j=0;jxij;Qin.close();for(i=0;im;i+)IY=0;IM=0;for(k=0;kK1;k+)Sk0=Ai; INF=0;for(l=0;l12;l+)Skl+1=Skl+xkl*0.0263-QP; if(Skl+1=Vn) Skl+1=Vn; if(INF0) IY=IY+1; IM=IM+INF;S1k=Sk12;FYi=IY/float(K1);FMi=IM/float(K1*12);/*推求水库蓄水状态的年转移概率*/for(j=0;jm;j+)INDj=0;for(k=0;kK1;k+)if(S1k=0.0)IND0+=1;else if(S1k=Vn)INDm-1+=1;elsee=int(S1k/DV)+1; INDe+=1;for(j=0;jm;j+)Qij=INDj/float(K1);/*推求水库稳定的蓄水概率分布*/float PAIm,PIm;for(j=0;jm;j+)if(j=0)PAIj=1.0;elsePAIj=0.0;for(j=0;jm;j+)PIj=0.0;for(i=0;im;i+)PIj=PIj+PAIi*Qij; for(e=0;e=0.001)PAIj=PIj;/*计算水库正常供水保证率*/float HY=0.0,HM=0.0;float PY,PM;for(i=0;im;i+)HY+=PIi*FYi;HM+=PIi*FMi;PY=1-HY;PM=1-HM;/*输出文件*/ofstream file1(H:状态转移矩阵.txt);for(i=0;im;i+)for(j=0;jm;j+)file1Qijt; file1n;file1.close();ofstream file2(H:条件破坏概率.txt);for(i=0;im;i+)file2FYitFMiendl;file2.close();ofstream file3(H:水库稳定蓄水概率分布.txt);for(i=0;im;i+)file3PIiendl;file3.close();ofstream file4(H:状态点各值.txt); for(i=0;im;i+) file4Aiendl;file4.close();ofstream file5(H:正常供水保证率.txt); file5PYtPMendl;file5.close();输入数据：水位库容关系：水位库容水位库容水位库容水位库容1551.361703.31856.720011.71561.461713.471866.9920112.11571.561723.651877.2920212.521581.661733.831887.5920312.961591.761744.021897.920413.421601.871754.221908.2220513.91611.991764.431918.5420614.381622.121774.651928.8620714.861632.251784.881939.1920815.381642.391795.111949.5220915.931652.531805.351959.8621016.51662.671815.619610.2121117.071672.821825.8619710.5721217.651682.971836.1319810.9421318.241693.131846.4119911.3221418.84径流资料年份四月五月六月七月八月九月十月十一月十二月一月二月三月1951-195271761526515002151621911361071041203121952-1953545154041087315101110621288165851304501953-19544835925909582361827973691151901271331954-19557831440160028601480142195911391001593411955-195626295723104466523361061125256961461956-1957584134039047964122816175441402262051957-1958717844652101014508585203322148951531958-19594171600631936113050779112877652443681959-19607467087892761541221573982181101172401960-196143052512601210130187902706546607661961-19625683714423892271384217282571671142721962-19635401250150050654536142550431610688851963-19649531470813163012401194646192473001794791964-196512901150209074754358462331112093941391965-1966636490626310511630690151169861211231966-19676842828723951627048413872601863691967-19684932010165011005954345774242141181246251968-19698097813571230532950411410166131991951969-19705428376732000103071827432587741231901970-1971461119089412102227431901662062072471861971-1972592900127035747936239721486751813881972-197363712707211981004397806541171473302761973-19741110102014907335081180186805548701291974-197558112309581140565259573987877771931975-197668013701420485244240356539188911523251976-1977579110014901020153212441406113711425251977-1978145011701240153065418344253610322576991978-197921660512603994153761596901337782661979-1980104498130770521639105545279851591980-1981481799167118692459284749180151891191701981-198280867377919725811412723885521572881982-19835167591406560109710353976482552091981831983-1984680918114722655296945503001087557691984-198528361411059476343376811722831171583891985-19865907037597412351881572741651401641891986-19875705171163152627125614015211410877811987-1988411778111414051040401628249122110821531988-198910258659428447611431769268681744341989-199011785359213862117145726591801432965181990-1991572821965654144891763791251702753721991-199240098779320885574451408176551073901992-199340185795350918314911069661612643181993-199442664665415901582509280486193106971221994-19954943568443001372074452543532673381931995-19964858701856105073218051217699103861611996-1997256566157622545535604517111411182843241997-19987345034867983081831991572032011924481998-1999735117385819821712325287125828163521999-200058266021781740414392209468137129177262程序输出结果离散状态点各值00.4731.4182.3633.3084.2535.1986.1437.0888.0338.9789.92310.86811.81312.75813.70314.64815.120状态转移矩阵：0.5100.0410.0410.1020.0200.0410.0610.0410.0410000.02000.0410.02000.0200.5100.0410.0410.1020.0200.0410.0610.0410.0410000.02000.0410.02000.0200.5100.0410.0410.1020.0200.0410.0610.0410.0410000.02000.0410.02000.0200.5100.0410.0410.1020.0200.0410.0610.0410.0410000.02000.0410.02000.0200.5100.0410.0410.1020.0200.0410.0610.0410.0410000.02000.0410.02000.0200.5100.0410.0410.1020.0200.0410.0610.0410.0410000.02000.0410.02000.0200.5100.0410.0410.1020.0200.0410.0610.0410.0410000.02000.0410.02000.0200.5100.0410.0410.1020.0200.0410.0610.0410.0410000.02000.0410.02000.0200.5100.0410.0410.1020.0200.0410.0610.0410.0410000.02000.0410.02000.0200.5100.0410.0410.1020.0200.0410.0610.0410.0410000.02000.0410.02000.0200.5100.0410.0410.1020.0200.0410.0610.