降水/径流等水文要素年内分配指标

1 不均匀系数nuniformity coefficient (Cn)

1.1 原理

1.2 MATLAB代码

MATLAB调用函数代码如下：

function Cn = GetCn( X )
% GetCn函数用于计算得到不均匀系数Cn
% 输入参数
% X  数据系列，可为降雨或径流
% 输出参数
% Cn  不均匀系数n = length(X)/12;
Xmean = mean( X );
XX = reshape(X,[],12);
sigma = zeros(n,1);
for iN=1:nsigma(iN,1) = ( sum( (XX(iN,:)-Xmean).^2 )/12 )^0.5;
end
Cn = sigma./Xmean;
end

1.3 案例

MATLAB主函数代码如下：

clc
close all
clear
%% 导入数据：某站点1954-2011年共58年月径流数据
load('X.mat');
Xlength = length( X(:,1) );
yearStart =1954;
yearEnd = 2011;
nYear = yearEnd-yearStart+1;%% 计算气象水文要素年内分配指数1-不均匀系数Cn
% 指数计算------------------------------------
Cn = GetCn( X );% 线性拟合
% [P,S] = POLYFIT(X,Y,N)   polyfit中S变量
[P , S] = polyfit(1:nYear,Cn' ,1);
Cnfit = polyval(P, 1:nYear );% 绘图-----------------------------------------
figure(1)
hold on;box on;
h(1) = plot(1:nYear,Cn,'k.-','linewidth',1.5,'markersize',10);
h(2) = plot(1:nYear,Cnfit,'r--','linewidth',1.5);
xlabel("Year");
ylabel("Cn");
set(gca, 'XTick', [2:5 : 57],'XTickLabel',{'1955','1960','1965','1970','1975','1980','1985','1990','1995','2000','2005','2010'}) ;
set(gca,'Layer','top','FontSize',12,'Fontname', 'Times New Roman');

图形如下：

2 完全调节系数complete accommodation coefficient (Cc)

2.1 原理

2.2 MATLAB代码

MATLAB调用函数代码如下：

function Cc = GetCc( X )
% GetCc函数用于计算得到完全调节系数Cc
% 输入参数
% X  数据系列，可为降雨或径流
% 输出参数
% Cc  完全调节系数n = length(X)/12;
Xmean = mean( X );
for in=1:length(X)if X(in,1)>=Xmeanj(in) = 1;elsej(in) = 0;end
endXX = reshape(X,[],12);
jj =reshape(j,[],12);
for in=1:nCc(in,1) = sum( jj(in,:).*(XX(in,:)-Xmean ) )/sum( XX(in,:) );
endend

2.3 案例

MATLAB主函数代码如下：

clc
close all
clear
%% 导入数据：某站点1954-2011年共58年月径流数据
load('X.mat');
Xlength = length( X(:,1) );
yearStart =1954;
yearEnd = 2011;
nYear = yearEnd-yearStart+1;%% 计算气象水文要素年内分配指数2-完全调节系数Cc
% 指数计算------------------------------------
Cc = GetCc( X );% 线性拟合
% [P,S] = POLYFIT(X,Y,N)   polyfit中S变量
[P , S] = polyfit(1:nYear,Cc' ,1);
Ccfit = polyval(P, 1:nYear );% 绘图-----------------------------------------
figure(2)
hold on;box on;
h(1) = plot(1:nYear,Cc,'k.-','linewidth',1.5,'markersize',10);
h(2) = plot(1:nYear,Ccfit,'r--','linewidth',1.5);
xlabel("Year");
ylabel("Cc");
set(gca, 'XTick', [2:5 : 57],'XTickLabel',{'1955','1960','1965','1970','1975','1980','1985','1990','1995','2000','2005','2010'}) ;
set(gca,'Layer','top','FontSize',12,'Fontname', 'Times New Roman');

图形如下：

3 集中度concentration degree (Cd)

3.1 原理

3.2 MATLAB代码

MATLAB调用函数代码如下：

function Cd = GetCd( X )
% GetCd函数用于计算得到集中度
% 输入参数
% X  数据系列，可为降雨或径流
% 输出参数
% Cd   集中度 concentration degreen = length(X)/12;
XX = reshape(X,[],12);
for in=1:nfor imonth=1:12Xxtemp(imonth,1) = XX(in,imonth)*cos(imonth*2*pi/12);Xytemp(imonth,1) = XX(in,imonth)*sin(imonth*2*pi/12);endXx = sum( Xxtemp );Xy = sum( Xytemp );Xxy = ( Xx^2+Xy^2 )^0.5;Cd(in,1) = Xxy/sum( XX(in,:) );
endend

3.3 案例

MATLAB主函数代码如下：

clc
close all
clear
%% 导入数据：某站点1954-2011年共58年月径流数据
load('X.mat');
Xlength = length( X(:,1) );
yearStart =1954;
yearEnd = 2011;
nYear = yearEnd-yearStart+1;
%% 计算气象水文要素年内分配指数3-集中度Cd
% 指数计算------------------------------------
Cd = GetCd( X );% 线性拟合
% [P,S] = POLYFIT(X,Y,N)   polyfit中S变量
[P , S] = polyfit(1:nYear,Cd' ,1);
Cdfit = polyval(P, 1:nYear );% 绘图-----------------------------------------
figure(3)
hold on;box on;
h(1) = plot(1:nYear,Cd,'k.-','linewidth',1.5,'markersize',10);
h(2) = plot(1:nYear,Cdfit,'r--','linewidth',1.5);
xlabel("Year");
ylabel("Cd");
set(gca, 'XTick', [2:5 : 57],'XTickLabel',{'1955','1960','1965','1970','1975','1980','1985','1990','1995','2000','2005','2010'}) ;
set(gca,'Layer','top','FontSize',12,'Fontname', 'Times New Roman');

图形如下：

4 集中期concentration period (Cp)

4.1 原理

4.2 MATLAB代码

MATLAB调用函数代码如下：

function Cp = GetCp( X )
% GetCp函数用于计算得到集中期
% 输入参数
% X  数据系列，可为降雨或径流
% 输出参数
% Cp   集中期 concentration periodn = length(X)/12;
XX = reshape(X,[],12);
for in=1:nfor imonth=1:12Xxtemp(imonth,1) = XX(in,imonth)*cos(imonth*2*pi/12);Xytemp(imonth,1) = XX(in,imonth)*sin(imonth*2*pi/12);endXx = sum( Xxtemp );Xy = sum( Xytemp );Cp(in,1) = atan( Xy/Xx);
end%求集中期，求得的是（-pi/2,pi/2)的值，要根据Xx，Xy来判定象限，再求侯数（度数）
D =(180*Cp/pi+180)';end

4.3 案例

MATLAB主函数代码如下：

clc
close all
clear
%% 导入数据：某站点1954-2011年共58年月径流数据
load('X.mat');
Xlength = length( X(:,1) );
yearStart =1954;
yearEnd = 2011;
nYear = yearEnd-yearStart+1;%% 计算气象水文要素年内分配指数4-集中期Cp
% 指数计算------------------------------------
Cp = GetCp( X );% 线性拟合
% [P,S] = POLYFIT(X,Y,N)   polyfit中S变量
[P , S] = polyfit(1:nYear,Cp' ,1);
Cpfit = polyval(P, 1:nYear );% 绘图-----------------------------------------
figure(4)
hold on;box on;
h(1) = plot(1:nYear,Cp,'k.-','linewidth',1.5,'markersize',10);
h(2) = plot(1:nYear,Cpfit,'r--','linewidth',1.5);
xlabel("Year");
ylabel("Cp");
set(gca, 'XTick', [2:5 : 57],'XTickLabel',{'1955','1960','1965','1970','1975','1980','1985','1990','1995','2000','2005','2010'}) ;
set(gca,'Layer','top','FontSize',12,'Fontname', 'Times New Roman');

图形如下：

5 相对变化幅度relative variation range (Cr)

5.1 原理

5.2 MATLAB代码

5.3 案例

6 绝对变化幅度absolute variation range (Ca)

6.1 原理

6.2 MATLAB代码

6.3 案例

7 洛伦兹不对称系数

7.1 原理

7.2 MATLAB代码

MATLAB调用函数代码：

7.3 案例

8 基尼系数Cg

8.1 原理

8.2 MATLAB代码

MATLAB调用函数代码：

function Cg = GetCg( X )
% GetCg函数用于计算得到基尼系数
% 输入参数
% X  数据系列，可为降雨或径流n = length(X)/12;
XX = reshape(X,[],12);
Xmean =mean( XX, 1 );
makeplot = false;  %  不绘图
V=[];
for i=1:nval = X(12*i-11:12*i);V = [V val];Cg(i,1) = gini(Xmean,val,makeplot);
endend

Gini函数调用代码：

% GINI computes the Gini coefficient and the Lorentz curve.
%
% Usage:
%   g = gini(pop,val)
%   [g,l] = gini(pop,val)
%   [g,l,a] = gini(pop,val)
%   ... = gini(pop,val,makeplot)
%
% Input and Output:
%   pop     A vector of population sizes of the different classes.
%   val     A vector of the measurement variable (e.g. income per capita)
%           in the diffrerent classes.
%   g       Gini coefficient.
%   l       Lorentz curve: This is a two-column array, with the left
%           column representing cumulative population shares of the
%           different classes, sorted according to val, and the right
%           column representing the cumulative value share that belongs to
%           the population up to the given class. The Lorentz curve is a
%           scatter plot of the left vs the right column.
%   a       Same as l, except that the components are not normalized to
%           range in the unit interval. Thus, the left column of a is the
%           absolute cumulative population sizes of the classes, and the
%           right colun is the absolute cumulative value of all classes up
%           to the given one.
%   makeplot  is a boolean, indicating whether a figure of the Lorentz
%           curve should be produced or not. Default is false.
%
% Example:
%   x = rand(100,1);
%   y = rand(100,1);
%   gini(x,y,true);             % random populations with random incomes
%   figure;
%   gini(x,ones(100,1),true);   % perfect equality
%
% Explanation:
%
%   The vectors pop and val must be equally long and must contain only
%   positive values (zeros are also acceptable). A typical application
%   would be that pop represents population sizes of some subgroups (e.g.
%   different countries or states), and val represents the income per
%   capita in this different subgroups. The Gini coefficient is a measure
%   of how unequally income is distributed between these classes. A
%   coefficient of zero means that all subgroups have exactly the same
%   income per capital, so there is no dispesion of income; A very large
%   coefficient would result if all the income accrues only to one subgroup
%   and all the remaining groups have zero income. In the limit, when the
%   total population size approaches infinity, but all the income accrues
%   only to one individual, the Gini coefficient approaches unity.
%
%   The Lorenz curve is a graphical representation of the distribution. If
%   (x,y) is a point on the Lorenz curve, then the poorest x-share of the
%   population has the y-share of total income. By definition, (0,0) and
%   (1,1) are points on the Lorentz curve (the poorest 0% have 0% of total
%   income, and the poorest 100% [ie, everyone] have 100% of total income).
%   Equal distribution implies that the Lorentz curve is the 45 degree
%   line. Any inequality manifests itself as deviation of the Lorentz curve
%   from the  45 degree line. By construction, the Lorenz curve is weakly
%   convex and increasing.
%
%   The two concepts are related as follows: The Gini coefficient is twice
%   the area between the 45 degree line and the Lorentz curve.
%
% Author : Yvan Lengwiler
% Release: $1.0$
% Date   : $2010-06-27$
function [g,l,a] = gini(pop,val,makeplot)% check argumentsassert(nargin >= 2, 'gini expects at least two arguments.')if nargin < 3makeplot = false;endassert(numel(pop) == numel(val), ...'gini expects two equally long vectors (%d ~= %d).', ...size(pop,1),size(val,1))pop = [0;pop(:)]; val = [0;val(:)];     % pre-append a zeroisok = all(~isnan([pop,val]'))';        % filter out NaNsif sum(isok) < 2warning('gini:lacking_data','not enough data');g = NaN; l = NaN(1,4);return;endpop = pop(isok); val = val(isok);assert(all(pop>=0) && all(val>=0), ...'gini expects nonnegative vectors (neg elements in pop = %d, in val = %d).', ...sum(pop<0),sum(val<0))% process inputz = val .* pop;[~,ord] = sort(val);pop    = pop(ord);     z    = z(ord);pop    = cumsum(pop);  z    = cumsum(z);relpop = pop/pop(end); relz = z/z(end);% Gini coefficient% We compute the area below the Lorentz curve. We do this by% computing the average of the left and right Riemann-like sums.% (I say Riemann-'like' because we evaluate not on a uniform grid, but% on the points given by the pop data).%% These are the two Rieman-like sums:%    leftsum  = sum(relz(1:end-1) .* diff(relpop));%    rightsum = sum(relz(2:end)   .* diff(relpop));% The Gini coefficient is one minus twice the average of leftsum and% rightsum. We can put all of this into one line.g = 1 - sum((relz(1:end-1)+relz(2:end)) .* diff(relpop));% Lorentz curvel = [relpop,relz];a = [pop,z];if makeplot   % ... plot it?area(relpop,relz,'FaceColor',[0.5,0.5,1.0]);    % the Lorentz curvehold onplot([0,1],[0,1],'k');                        % 45 degree lineaxis tight      % ranges of abscissa and ordinate are by definition exactly [0,1]axis square     % both axes should be equally longset(gca,'XTick',get(gca,'YTick'))   % ensure equal tickingset(gca,'Layer','top');             % grid above the shaded area% grid on;title(['\bfGini coefficient = ',num2str(g)]);xlabel('累积时间/%');ylabel('累计降雨量/%');end
end

8.3 案例

MATLAB主函数代码如下：

clc
close all
clear
%% 导入数据：某站点1954-2011年共58年月径流数据
load('X.mat');
Xlength = length( X(:,1) );
yearStart =1954;
yearEnd = 2011;
nYear = yearEnd-yearStart+1;%% 计算气象水文要素年内分配指数8-基尼系数Cg
Cg = GetCg( X );% 线性拟合
% [P,S] = POLYFIT(X,Y,N)   polyfit中S变量
[P , S] = polyfit(1:nYear,Cg' ,1);
Cgfit = polyval(P, 1:nYear );% 绘图-----------------------------------------
figure(8)
hold on;box on;
h(1) = plot(1:nYear,Cg,'k.-','linewidth',1.5,'markersize',10);
h(2) = plot(1:nYear,Cgfit,'r--','linewidth',1.5);
xlabel("Year");
ylabel("Cg");
set(gca, 'XTick', [2:5 : 57],'XTickLabel',{'1955','1960','1965','1970','1975','1980','1985','1990','1995','2000','2005','2010'}) ;
set(gca,'Layer','top','FontSize',12,'Fontname', 'Times New Roman');

图形如下：

9 差积曲线

10 丰平枯统计分析

参考资料

1.论文-J2017-Inter- and intra- annual environmental flow alteration and its implication in the Pearl River Delta, South China
2.论文-J1989-水文时间序列不均匀系数的分析与计算
3.代码参考-Matlab计算气象水文要素年内分配指数

MATLAB计算气象水文要素年内分配指数相关推荐

1. 应用matlab计算线性定常系统的矩阵指数

   syms s t A = [5 0 3;0 5 2;9 1 3]; Fs = inv(seye(3)-A); %求预解矩阵FS=[(sI-A)]^(-1),eye(3)为33矩阵 eAt = ilap ...

2. Matlab处理气象数据（四）观测数据的预处理和计算

   观测数据为中国气象数据网上获得的中国地面温度月值0.5°×0.5°格点数据集,时间范围是1961年1月至2013年12月.这套数据为txt格式,包含头文件.头文件信息为: NCOLS 128 NROW ...

3. 【物理应用】基于matlab GUI气象参数计算综合指标和IAQI【含Matlab源码 2116期】

   ⛄一.获取代码方式 获取代码方式1: 完整代码已上传我的资源:[物理应用]基于matlab GUI气象参数计算综合指标和IAQI[含Matlab源码 2116期] 点击上面蓝色字体,直接付费下载,即可 ...

4. Matlab处理气象数据（十一）数据的异常值计算

   %平均温度的异常值计算 load('Tem1.mat');%导入NCEP数据的面积加权年平均 load('Tem2.mat');%导入观测数据的面积加权年平均 m1=mean(Tem1); %求Tem ...

5. matlab画复变函数,科学网—复数复变函数的Matlab计算与绘图 - 周铁戈的博文

   复数复变函数的Matlab计算与绘图 周铁戈 复数的表示 存在两种表示方法,一种是代数式,一种是指数式,在Matlab中的方式如下: >> z=1+2i            #代数式,1 ...

6. P2 Matlab计算基础-《Matlab/Simulink与控制系统仿真》程序指令总结

   上一篇 回到目录 下一篇 <Matlab/Simulink与控制系统仿真>程序指令总结 Matlab_Simulink_BookExample 2. Matlab 计算基础 表2.1 Ma ...

7. matlab计算复活节概率,复活节日期的计算方法

   复活节(主复活日)是一个西方的重要节日,在每年春分月圆之后第一个星期日.基督徒认为,复活节象征着重生与希望,为纪念耶稣基督于公元30到33年之间被钉死在十字架之后第三天复活的日子. 算法 复活节是西方 ...

8. matlab——计算VPD（vapor pressure defict）

   需求: 计算VPD(vapor pressure defict). 介绍: 饱和水汽压差(简称VPD) 是指在一定温度下,饱和水汽压与空气中的实际水汽压之间的差值(百度百科). 因此,温室中VPD的理 ...

9. matlab计算macd_[转载]4.K线图以及常用技术指标的Matlab实现-基于Matlab的量化投资...

   本次主要讲解用Matlab来实现K线图以及常用的技术指标. 相关前导帖子浏览: 系列帖子目录 Matlab的Bar图中Bar颜色灵活设置的一点总结[by faruto] [原创]坐标轴刻度标签旋转升级 ...

最新文章

1. 电路中滤波电容和退耦电容_电子电路中电容的作用，滤波消抖，充放电，耦合，退耦...
2. STM32F407VG uCOS-II2.91 IAR工程 以及uCOS使用库编译的方法
3. window.open 和showModalDialog的返回值
4. RTP格式图 NNEXB格式和RTP格式
5. Docker 遇到swapon failed Operation not permitted
6. Hibernate面试总结
7. 10.28-11.1-广州软件所-实习工作日记
8. 修改strut默认的action后缀
9. 中国GPS开发工具市场现状研究分析与发展前景预测报告（2022）
10. 四参数拟合曲线_Origin进行体外释药规律的拟合
11. datavideo切换台说明书_GoPro结合洋铭切换台现场节目制作
12. 计算机应用专业，报软考应该选什么？
13. 2021.11.20【读书笔记】|差异可变剪接事件及DTU分析
14. 马云在大学学计算机,IT大佬高考成绩单：李彦宏是状元 马云数学仅1分
15. 如何浏览与下载全球免费的地图高分辨率(亚米级)的遥感影像？
16. C++运算符重载 ++，--，+，-，+=，-=，输出输入运算符
17. 道里云公司网络虚拟化架构NVI技术开放源代码--序言
18. Chrome开发工具Network没有显示完整的http request和response对话
19. p5.js创意绘图(2)自画像
20. 从零入门云计算（1）：云计算究竟是个啥？

热门文章

1. VPO微珀宣布完成千万人民币Pre-A融资，专注电子烟场景化消费市场...
2. 如何统一Windows和Macbook的使用体验，无缝切换，避免精分？
3. OGNL表达式的基本语法和用法
4. 【Python】获取roc、auc时候报错：raise ValueError({0} format is not supported.format(y_type))
5. 狂神说 SpringBoot笔记
6. 前端加载动画-三点加载
7. 5分钟玩转Axure之中继器（表格篇）
8. Razer Fintech 委任林祥源先生为顾问委员会成员
9. 常用xshell5基本命令
10. 第四十五天 百度地图定位SDK