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(%)Cs0.115205080959999.9
0.03.092.331.640.840.00-0.84-1.64-2.33-3.09
0.13.231.672.00.84-0.02-0.85-1.62-2.25-2.95
0.23.382.471.700.83-0.03-0.85-1.59-2.18-2.81
0.33.522.541.730.82-0.05-0.85-1.55-2.10-2.67
0.43.672.621.750.82-0.07-0.85-1.52-2.03-2.54
0.53.812.681.770.81-0.08-0.85-1.40-1.96-2.40
0.63.962.751.800.80-0.10-0.85-1.45-1.88-2.27
0.74.102.821.820.79-0.12-0.85-1.42-1.81-2.14
0.84.242.891.840.78-0.13-0.85-1.38-1.74-2.02
0.94.392.961.860.77-0.15-0.85-1.35-1.66-1.90
1.04.533.021.880.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对应的三个参数进行计算。
参考：程银才等，一种新的水文频率计算方法，水文，，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的优化工具进行曲线拟合，本人主要用了三种方法：matlab自带的非线性最小二乘拟合函数lsqcurvefit
这种方法，是根据实际的观测值(排序过后)计算出的理论频率值与经验频率值进行拟合的方法。
%计算经验频率曲线
real_data=sort(max_data,‘descend’);
fori=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) ;
matlab带有的全局优化方法Multistart
%% Global optimization
ms = MultiStart;
C = run(ms, problem, 100);
上述两种方法，拟合结果都一样。
matlab 带有的带条件优化方法fmincon
为什么需要这种方法，看下面这张通过上面两种方法拟合的图：不难发现，在曲线的末端会出现负值，与水文现象的物理意义出现矛盾，这种拟合结果是因为Cs<2Cv,同时当Cs>2时，密度曲线呈乙型，因此参数合理范围，2Cv
因此，正确的拟合方式，应该是在上面的条件下进行拟合，对于这种带有条件的非线性拟合，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);
这样拟合的结果，就不会出现上述情况。当然如果在实际的应用中，没有出现上面的结果，应该同时采用不同的方法，然后比较不同方法产生的均方根误差，来最后决定最优的方法、
fori=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-
% Created by Ling Zhang, zhanglingky@, CAREERI, /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-
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;
转载本文请联系原作者获取授权，同时请注明本文来自张凌科学网博客。
链接地址：/blog-922140-884754.html
上一篇：高斯非线性最小二乘法(全局优化)拟合
下一篇：Hydrologic Vs Hydraulic Modelling
