P3分布推求设计频率年降水量的python程序

最新推荐文章于 2024-02-21 22:34:25 发布

eco_hydro

最新推荐文章于 2024-02-21 22:34:25 发布

阅读量5.4k

点赞数

本文链接：https://blog.csdn.net/qq_33015657/article/details/108630737

版权

水文上常用皮尔逊三型分布（P3分布）推算设计频率的水文变量，比如降水量、流量、径流量，等等。本文给出了一个利用P3分布推算设计频率下的年降水量的python程序。

P3分布的概率密度函数如下：

$f(x)=\tfrac{\beta ^{\alpha }}{\Gamma (\alpha )}(x-x_{0})^{\alpha -1}e^{-(x-x_{0})\beta }$

其中：

$\alpha =\tfrac{4}{Cs^{2}}$

$\beta =\tfrac{2}{E(x)CsCv}$

$x_{0}=E(x)(1-\tfrac{2Cv}{Cs})$

Cv是变差系数， $\large Cv=\tfrac{\sigma }{\bar{x}}$ ，Cs是偏态系数。

以下是程序全部代码，代码运行所需文件已经上传。

import pandas as pd
import math as mt
import numpy as np


ws=pd.read_excel('年降水量系列.xlsx')
pdata=ws.iloc[:,1]#导入降水系列数据

#计算样本统计参数
mean=np.mean(pdata) #降水系列的均值
std=np.std(pdata) #降水系列的标准偏差
Cs=pdata.skew() #降水系列的偏态系数
Cv=std/mean #降水系列的变差系数
max=max(pdata) #最大值
min=min(pdata) #最小值

#计算P3分布统计参数
alpha=4/Cs**2 #参数alpha
beta=2/(mean*Cs*Cv) #参数beta
x0=mean*(1-2*Cv/Cs) #参数x0

pmp=max+3*std #构造极端最大降水量
pseq=[i for i in range(mt.ceil(pmp/10)*10,mt.floor(x0/10)*10,-10)] #构造降水等差序列
seq=[i for i in range(1,len(pseq)+1)] #将降水等差序列排序

quantiles=[]#计算p3分布各分位点的概率密度函数值
for i in range(len(pseq)):
    quantile=beta**alpha/mt.gamma(alpha)*(pseq[i]-x0)**(alpha-1)*mt.exp(-(pseq[i]-x0)*beta)
    quantiles.append(quantile)
    
probabilities=[]#计算各分位点之间的区间概率
probabilities.append(0)
for i in range(1,len(quantiles)):
    probability=(quantiles[i-1]+quantiles[i])*(pseq[i-1]-pseq[i])/2
    probabilities.append(probability)
    
accum_probblts=[]#计算累计概率
for i in range(len(probabilities)):
    accum_probblts.append(sum(probabilities[:i+1]))

def Linterpo(x,y,xi):
    if xi<x[0]:
        yi=y[0]+(xi-x[0])/(x[0]-x[1])*(y[0]-y[1])
    if xi>x[-1]:
        yi=y[-1]+(xi-x[-1])/(x[-1]-x[-2])*(y[-1]-y[-2])
    else:
        for j in range(len(x)):
            if xi>=x[j] and xi<x[j+1]:
                yi=y[j]+(xi-x[j])/(x[j+1]-x[j])*(y[j+1]-y[j])
    return yi

p=[1/10000,1/1000,1/100,1/50,1/30,1/20,1/10,1/5] #设计频率
for i in range(len(p)):
    pp=round(Linterpo(accum_probblts,pseq,p[i]),1)
    print(str(round(p[i]*100,1))+'%设计频率的年降水量是'+str(pp)+'mm')

程序运行结果如下：