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 皮尔逊Ⅲ型频率曲线的离均系数 值表（摘录）

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 Ë®ÎÄÆàÂÛÇúÏßÏßDíÑD¾¿

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;

转载请注明：笑凌子 | » 水文频率曲线及MATLAB绘制

前面的基本概念主要拷贝于：http://jpkt.whu.edu.cn/jpkc2008/gcswx/swx/cc/hydrology04/4_4.htm （特别感谢）