水文频率曲线及MATLAB绘制

水文频率曲线及MATLAB绘制

 

  1. 水文频率曲线

 水文分析计算中使用的概率分布曲线俗称水文频率曲线,习惯上把由实测资料(样本)绘制的频率曲线称为经验频率曲线,而把由数学方程式所表示的频率曲线称为理论频率曲线。所谓水文频率分布线型是指所采用的理论频率曲线(频率函数)的型式(水文中常用线型为正态分布型、极值分布型、皮尔逊Ⅲ型分布型等),它的选择主要取决于与大多数水文资料的经验频率点据的配合情况。分布线型的选择与统计参数的估算,一起构成了频率计算的两大内容。

2. Pearson Type III 频率曲线

2.1 皮尔逊型曲线的概率密度函数
    皮尔逊Ⅲ型曲线是一条一端有限一端无限的不对称单峰、正偏曲线(见图4-4-3),数学上常称伽玛分布,其概率密度函数为:
        

式中:Γα―α的伽玛函数;
     
αβ、a0分别为皮尔逊Ⅲ型分布的形状尺度和位置未知参数,
     
α﹥0, β﹥0 。

显然,三个参数确定以后,该密度函数随之可以确定。可以推论,这三个参数与总体三个参数 、Cv、CS具有如下关系:
        
                          (4-4-9)
        


2.2 皮尔逊型频率曲线及其绘制

水文计算中,一般需要求出指定频率P所相应的随机变量取值xp,也就是通过对密度曲线进行积分,即:
            (4-4-10)
    求出等于及大于xp的累积频率P值。直接由式(4-4-10)计算P值,非常麻烦,实际做法是通过变量转换,变换成下面的积分形式 :
            (4-4-11)
    式(4-4-11)中被积函数只含有一个待定参数CS,其它两个参数 、Cv都包含在 中。 ,x是标准化变量,称为离均系数。 的均值为0,标准差为1。因此,只需要假定一个CS值,便可从式(4-4-11)通过积分求出  之间的关系。对于若干个给定的CS值, 的对应数值表,已先后由美国福斯特和前苏联雷布京制作出来,见附表1″皮尔逊Ⅲ型频率曲线的离均系数 值表”。由 就可以求出相应频率 的x值: 
            (4-4-12)

 

附表1 皮尔逊型频率曲线的离均系数 值表(摘录)

P(%)Cs 

0.1

1

5

20

50  

80

95

 99 

99.9

0.0

3.09

2.33

1.64

0.84

0.00

-0.84

-1.64

-2.33

-3.09

0.1

3.23

1.67

2.0

0.84

-0.02

-0.85

-1.62

-2.25

-2.95

0.2

3.38

2.47

1.70

0.83

-0.03

-0.85

-1.59

-2.18

-2.81

0.3

3.52

2.54

1.73

0.82

-0.05

-0.85

-1.55

-2.10

-2.67

0.4

3.67

2.62

1.75

0.82

-0.07

-0.85

-1.52

-2.03

-2.54

0.5

3.81

2.68

1.77

0.81

-0.08

-0.85

-1.40

-1.96

-2.40

0.6

3.96

2.75

1.80

0.80

-0.10

-0.85

-1.45

-1.88

-2.27

0.7

4.10

2.82

1.82

0.79

-0.12

-0.85

-1.42

-1.81

-2.14

0.8

4.24

2.89

1.84

0.78

-0.13

-0.85

-1.38

-1.74

-2.02

0.9

4.39

2.96

1.86

0.77

-0.15

-0.85

-1.35

-1.66

-1.90

1.0

4.53

3.02

1.88

0.76

-0.16

-0.85

-1.32

-1.59

-1.79

 

3. 经验频率曲线

上述各种频率曲线是用数学方程式来表示的, 属于理论频率曲线。在水文计算中还有一种经验频率曲线, 是由实测资料绘制而成的, 它是水文频率计算的基础, 具有一定的实用性。根据实测水文资料,按从大到小的顺序排列,如图4-4-4所示,然后用经验频率公式计算系列中各项的频率,称为经验频率。以水文变量x为纵坐标,以经验频率为横坐标,点绘经验频率点据,根据点群趋势绘出一条平滑的曲线,称为经验频率曲线。

 对经验频率的计算,目前我国水文计算上广泛采用的是数学期望公式 :
           
 (4-4-13)
    式中 p – 等于和大于xm的经验频率;
         m – xm的序号,即等于和大于xm的项数;
         n – 系列的总项数。

   经验频率曲线计算工作量小,绘制简单,查用方便,但受实测资料所限,往往难以满足设计上的需要。为此,提出用理论频率曲线来配合经验点据,这就是水文频率计算适线(配线)法。

4. 频率曲线参数估计

在概率分布函数中都含有一些表示分布特征的参数, 例如皮尔逊III型分布曲线中就包含有 、Cv、Cs三个参数。水文频率曲线线型选定之后, 为了具体确定出概率分布函数, 就得估计出这些参数。目前,由样本估计总体参数的方法主要有矩法、三点法、权函数法等。

 

    矩法是用样本矩估计总体矩,并通过矩和参数之间的关系,来估计频率曲线参数的一种方法。
    前述,一阶原点矩的计算公式就是均值 ,均方差σ的计算式为二阶中心矩开方,偏态系数CS计算式中的分子则为三阶中心矩。由此,得到计算参数的公式(4-3-6)、(4-3-9)、(4-3-10)、(4-3-13)。它们与相应的总体同名参数不一定相等。
    但是,我们希望由样本系列计算出来的统计参数与总体更接近些,因此,需要将上述公式加以修正,修正后的参数计算式为:
                         (4-5-1)

                 (4-5-2)
           
  (4-5-3)
                  (4-5-4)
    水文计算上习惯称上述公式为无偏估值公式,并用它们估算总体参数,作为配线法的参考数值(配线法将在下面介绍)。

 

5. 水文频率曲线的拟合

5.1 目估适线法(传统方法)

目估配线法又称目估适线法,是以经验频率点据为基础,给它们选配一条符合较好的理论频率曲线,并以此来估计水文要素总体的统计规律。具体步骤如下:
 —- 将实测资料由大到小排列,计算各项的经验频率,在频率格纸上点绘经验点据(纵坐标为变量的取值,横坐标为对应的经验频率)
 —- 选定水文频率分布线型(一般选用皮尔逊Ⅲ型)。
 —- 先采用矩法或其它方法估计出频率曲线参数的初估值 Cv,而Cs凭经验初选为Cv的倍数。
 —- 根据拟定的 Cv和Cs,查附表1或附表2,计算xp值。以xp 为纵坐标, p为横坐标,即可得到频率曲线。将此线画在绘有经验点据的图上,看与经验点据配合的情况 。若不理想,可通过调整 、cv和cs点绘频率曲线。
 —- 最后根据频率曲线与经验点据的配合情况,从中选出一条与经验点据配合较好的曲线作为采用曲线,相应于该曲线的参数便看作是总体参数的估值。
 —- 求指定频率的水文变量设计值。

5.2 优化适线法

优化适线法是在一定的适线准则(即目标函数)下,求解与经验点据拟合最优的频率曲线的统计参数的方法。优化适线法按不同的适线准则分为三种,即离差平方和最小准则(OLS)、离差绝对值和最小准则(ABS)、相对离差平方和最小准则(WLS),其中以离差平方和最小准则(OLS)最为常用。

6. 水文频率曲线的MATLAB拟合

6.1 利用MATLAB进行PIII型水文频率曲线拟合,主要包括两类:

一类是仅仅将MATALB作为一种绘图工具,离均系数 通过查表的方法获取,然后 求出相应频率 的x值: 
            (4-4-12)

这种方法和传统基于EXCEL的方法相比,并没有体现出MATLAB的优势

 

    另一类是通过MATLAB身所带有的库函数及优化工具,对PIII水文频率曲线进行拟合。这种方法更能发挥MATLAB的优势,而且离均系数 不需要认为的进行查表,能直接计算出来,这样大大缩减了绘制水文频率曲线的时间。

 

6.2 下面我就主要介绍第二类MATLAB的拟合方法。

6.2.1 首先,用MATLAB如何进行理论频率的计算。

 

参考文献,计算方法包括两类,其实是两种不同的解法,即通过不同的变量代换方法求解。

一种解法仍然是通过离均系数,然后更加4-4-12的方程,计算某个频率下的x值,或者反计算出x值所对应的理论p值。

参考:林莺等,水文频率曲线简洁计算和绘图技巧,2002(33)

一种方法,是不通过离均系数,即直接通过4-4-9对应的三个参数进行计算。

参考:程银才等,一种新的水文频率计算方法,水文,2008,1(28)

 

相对来说,第一种方法,比较简单,因此后面我的例子,选择第一方法。

这种方法,设t=β(x-a0)将P-III累计频率简化为标准Gamma分布(如下图所示),然后利用标准Gamma分布的返回数,gammain计算不同设计频率对应的tp值,然后根据公式(3)计算离均系数,最后通过公式(4)计算x值。

通过(2)式求解tp,实际上是gamma函数的逆函数求解(gammainv):tp=gammainv(1-p,a,1),a=4/Cs^2


 

截图于:郑伟,基于MATLAB GUI 技术的P-III型频率曲线实现,科技论坛

 

同时可以根据(3),(4)反解除tp, 利用(2)式计算出不同x对应的p值,即p=1-gammacdf(tp,a,1)

6.2.2 P-III频率曲线拟合

利用MATLAB的优化工具进行曲线拟合,本人主要用了三种方法:

  1. matlab自带的非线性最小二乘拟合函数lsqcurvefit

这种方法,是根据实际的观测值(排序过后)计算出的理论频率值与经验频率值进行拟合的方法。

% 计算经验频率曲线

real_data=sort(max_data,‘descend’);

for i=1:n;

p_real_data(i)=length(find(real_data>=real_data(i)))/(n+1); % units: 1

end

% 计算理论频率曲线

1-gamcdf((2*(real_data-p)/(p*C(1)*C(2))+4/C(2)^2),4/C(2)^2,1)

% 进行lsqcurvefit拟合

%% lsqcurvefit

objectFun=@(C,real_data) 1-gamcdf((2*(real_data-p)/(p*C(1)*C(2))+4/C(2)^2),4/C(2)^2,1); % Cv=C(1);Cs=C(2)

x0=[Cv0,Cs0];

xdata=real_data;

ydata=p_real_data’;

problem = createOptimProblem(‘lsqcurvefit’,


‘objective’, objectFun,


‘xdata’, xdata, ‘ydata’, ydata,


‘x0’,x0,


‘lb’, [0 0],


‘ub’, [inf inf]);

% local optimization

C=lsqcurvefit(problem) ;

 

  1. matlab带有的全局优化方法Multistart

%% Global optimization

ms = MultiStart;

C = run(ms, problem, 100);

上述两种方法,拟合结果都一样。

 

  1. matlab 带有的带条件优化方法fmincon

为什么需要这种方法,看下面这张通过上面两种方法拟合的图:不难发现,在曲线的末端会出现负值,与水文现象的物理意义出现矛盾,这种拟合结果是因为Cs<2Cv,同时当Cs>2时,密度曲线呈乙型,因此参数合理范围,2Cv<Cs<2 (参考:梁季阳等1986 水文频率曲线的线型研究)


因此,正确的拟合方式,应该是在上面的条件下进行拟合,对于这种带有条件的非线性拟合,maltab带有fmincon函数。这里的目标函数需要更改:

objectFun=@(C) sum(((1-gamcdf((2*(real_data-p)/(p*C(1)*C(2))+4/C(2)^2),4/C(2)^2,1))’-p_real_data).^2); % Cv=C(1);Cs=C(2)

x0=[Cv0,Cs0];

A_ineq=[2 -1;0 2]; % refered to Jiyang Liang, 1986 Ë®ÎÄÆÀÂÛÇúÏßÏßÐÍÑо¿

b=[0;2];

C=fmincon(objectFun,x0,A_ineq,b);

 

这样拟合的结果,就不会出现上述情况。当然如果在实际的应用中,没有出现上面的结果,应该同时采用不同的方法,然后比较不同方法产生的均方根误差,来最后决定最优的方法、

 

for i=1:length(x);

tp2=gaminv(1-x(i)/100,4/C(2)^2,1);

Z1(i)=((C(2)*tp2/2-2/C(2))*C(1)+1)*p; % Resuls optimizated

end

 

tp_real=2*(real_data-p)/(p*C(1)*C(2))+4/C(2)^2;

P_theory=1-gamcdf(tp_real,4/C(2)^2,1);

RMSE=sqrt(sum((P_theory’-p_real_data).^2)/n); % error evaluation

 

 

附件 1 : 完整的P-III曲线MATLAB拟合代码

%——————————————————————————————————————————

% This Program was used to plot the Peason Type III (P-III) distrigution of

% the annnal max and min streamflows from 1979-2009

% Created by Ling Zhang, zhanglingky@lzb.ac.cn, CAREERI, 2015/4/16

%——————————————————————————————————————————

% We have tried three methods, namely lsqcurvefit, fmincon and Multistart, to fit the P-III

% curve. The results are same for both methods.

%——————————————————————————

 

%% Annal max and min streamflow data extraction from 1979-2009

data_xls=xlsread(‘Yingluoxia_Daily_Flow_SEP.xls’);

max_data=[];

temp_locat=0;

for i=1:(length(data_xls(:,3))-1);

if data_xls(i,3)>data_xls(i+1,3)

max_data_tmp=max(data_xls((temp_locat+1):i,6));

max_data=[max_data;max_data_tmp];

temp_locat=i;

end

end

max_end=max(data_xls((temp_locat+1):end,6));

max_data=[max_data;max_end];

%% Inital parameters calculation

n=length(max_data);

p=mean(max_data);

m=std(max_data); %s^2=[(x1-x)^2+(x2-x)^2+……(xn-x)^2]/(n-1)

Cv0=m/p;

Cs0=sum((max_data-p).^3)/((n-3)*m^3); %Cs, Coefficient of Skewness

A=4/Cs0^2;

% 计算不同频率的设计值

x=[0.01 0.05 0.5 1 2 5 10 20 30 40 50 60 70 80 90 95 98 99.5 99.9 99.99]; %units: %

for i=1:length(x);

tp1=gaminv(1-x(i)/100,A,1);

Z0(i)=((Cs0*tp1/2-2/Cs0)*Cv0+1)*p; % Resuls not optimizated

end

%计算经验频率

real_data=sort(max_data,’descend’);

for i=1:n;

p_real_data(i)=length(find(real_data>=real_data(i)))/(n+1); % units: 1

end

 

%% Prameters (Cv and Cs) Optimization

%% fmincon

 

% tp_real=2*(real_data-p)/(p*Cv*Cs)+4/Cs^2;

% P_theory=1-gamcdf(tp,4/Cs^2,1);

% Object=sqrt(sum((P_theory’-p_real_data).^2)/n); % error evaluation

objectFun=@(C) sum(((1-gamcdf((2*(real_data-p)/(p*C(1)*C(2))+4/C(2)^2),4/C(2)^2,1))’-p_real_data).^2); % Cv=C(1);Cs=C(2)

x0=[Cv0,Cs0];

A_ineq=[2 -1;0 2]; % refered to Jiyang Liang, 1986 水文评论曲线线型研究

b=[0;2];

C=fmincon(objectFun,x0,A_ineq,b);

 

%% lsqcurvefit

objectFun=@(C,real_data) 1-gamcdf((2*(real_data-p)/(p*C(1)*C(2))+4/C(2)^2),4/C(2)^2,1); % Cv=C(1);Cs=C(2)

x0=[Cv0,Cs0];

xdata=real_data;

ydata=p_real_data’;

problem = createOptimProblem(‘lsqcurvefit’, …

‘objective’, objectFun, …

‘xdata’, xdata, ‘ydata’, ydata, …

‘x0’,x0,…

‘lb’, [0 0],…

‘ub’, [inf inf]);

% local optimization

C=lsqcurvefit(problem) ;

 

%% Global optimization

ms = MultiStart;

C = run(ms, problem, 100);

 

%% 计算不同频率的设计值 (Opitimizated results)

for i=1:length(x);

tp2=gaminv(1-x(i)/100,4/C(2)^2,1);

Z1(i)=((C(2)*tp2/2-2/C(2))*C(1)+1)*p; % Resuls optimizated

end

 

tp_real=2*(real_data-p)/(p*C(1)*C(2))+4/C(2)^2;

P_theory=1-gamcdf(tp_real,4/C(2)^2,1);

RMSE=sqrt(sum((P_theory’-p_real_data).^2)/n); % error evaluation

 

%% Another way to set normal frequency grids

y=norminv(x./100,0,1); %数据y所对应的正态分布的 X 值

% x./100 ia the requency

y=y-norminv(0.0001,0,1); % axis conversions (Theoritical frequency)

set(gca,’xtick’,[],’xlim’,[0 y(end)],’ylim’,[0 ,1900]);%1900, max y axis which should chagne according to the observed data

set(gca,’xtick’,y);

set(gca,’xticklabel’,x);

grid on;

hold on;

 

%% Plot results

yy=norminv(p_real_data,0,1); %数据p_real_data所对应的正态分布的 X 值

yy=yy-norminv(0.0001,0,1); % axis conversions (Emperical frequency)

plot(yy,real_data,’r*’);

hold on;

plot(y,Z0,’g’); % Resuls not optimizated

hold on;

plot(y,Z1,’r’); % Resuls optimizated

hold off;

Updated: 2015-04-21 — am10:33

1
Leave a Reply

1 Comment threads
0 Thread replies
0 Followers
 
Most reacted comment
Hottest comment thread
1 Comment authors
  Subscribe  
newest oldest most voted
Notify of

好高深……