近期忙，忙于生计，三千出头的工资早已力不从心~公众号自上次更新以来已久久没有动静。今天决定再更新一篇：关于水文P-III曲线的一些内容。P-III曲线这里只讲一下绘图的问题。最近在做了几个频率曲线适线的事情，比较了手里好几个P-III曲线绘图软件：里面的几种都不失为优秀的P-III曲线分析和绘图系统，除了可以直接有系列值计算频率及绘图外，还可以实现均值、Cv、Cs的实时调整与同步绘图，有的工具甚至有一些可以用到的参数优化方法。但有时候还是有些不方便：如对经验点的定制化(定制化频率)、对微调的步长大小限制，当然其中一些工具excel版的可以通过二开实现。因此就涉及到定制化的要求。结合在最近的适线过程中个人认为：在某些特殊情况下应手动计算各经验点，然后再由绘图系统去单纯的绘图适线。因此可以计算出经验点后由仅工具二、三、四的配图系统去适线。说到定制化可以想到编程：结合vba的excel、MATLAB、Python。下面给出MATLAB、Python的P-III曲线绘图系统代码。代码来源Github，在此向原作者感谢。MATLAB版

function HessianProbabilityGrid(dataArray, estCv, estCsDivCv)    % 控制Y轴坐标上限    MAX_Y = max(dataArray);    MAX_Y = ceil(MAX_Y / 10^floor(log10(MAX_Y))) * 10^floor(log10(MAX_Y));    % 定义坐标函数    axisProb = [0.001, 0.01:0.01:0.1, 0.2:0.1:1, 1.2:0.2:2, 3:1:20, 22:2:80, 81:1:98, 98.2:0.2:99, 99.1:0.1:99.9]; %横坐标的值    axisX = norminv(axisProb / 100, 0, 1);    axisX = axisX - axisX(1);    % 纵向划分线    xLabel = [0.001 0.01 0.1 0.2 0.333 0.5 1 2 3 5 10 20 25 30 40 50 60 70 75 80 85 90 95 97 99 99.9 100];    xTick = norminv(xLabel / 100, 0, 1);    xTick = xTick - xTick(1);    for i = 1:1:size(xTick, 2)        line([xTick(i), xTick(i)], [0, MAX_Y]);    end    % 横向划分线    for i = MAX_Y / 10:MAX_Y / 10:MAX_Y        line([0, axisX(end)], [i, i]);    end    % 绘制坐标格纸    h = findobj('type', 'axes');    % 去掉x,y轴的刻度，并设置x,y轴的上下界    set(h, 'xtick', [], 'ytick', [], 'xlim', [0, axisX(end)], 'ylim', [0, MAX_Y]);    % 横坐标的值    xLabel = [0.001 0.01 0.1 0.2 0.333 0.5 1 2 3 5 10 20 25 30 40 50 60 70 75 80 85 90 95 97 99 99.9 100];    xTick = norminv(xLabel / 100, 0, 1);    xTick = xTick - xTick(1);    % 重新设置坐标轴的刻度    set(h, 'Xtick', xTick);    set(h, 'XticklabeL', {0.001 0.01 0.1 0.2 0.3 0.5 1 2 3 5 10 20 25 30 40 50 60 70 75 80 85 90 95 97 99 99.9 100});    set(h, 'Ytick', 0:MAX_Y / 10:MAX_Y);    set(h, 'YticklabeL', {0:MAX_Y / 10:MAX_Y});    % 对x,y轴进行标注    xlabel('P(％)');    ylabel('Qm(m^3/s)');    hold on;    % 原始数据处理    sortedData = sort(dataArray, 'descend'); % 排序后的数据    EX = mean(dataArray);    expProb = [1:length(sortedData)] / (length(sortedData) + 1);    % 实际数据横坐标的值    dataX = norminv(expProb, 0, 1);    dataX = dataX - norminv(0.00001, 0, 1);    % 绘制实际数据散点图    scatter(dataX, sortedData, 'k');    hold on;    % 绘制理论曲线    estCs = estCsDivCv * estCv;    % 估计的P-III型曲线参数    P3_alpha = 4 / estCs^2;    P3_beta = 2 / (EX * estCv * estCs);    P3_a0 = EX * (1 - 2 * estCv / estCs);    % 估计的gamma函数参数    gamma_alpha = P3_alpha;    gamma_beta = 1 / P3_beta;    theoryY = gaminv((100 - xLabel) / 100, gamma_alpha, gamma_beta) + P3_a0; % 注意超过制概率和不及制概率    plot(xTick, theoryY, 'k');    hold on;    % grid on;end

Python版

# -*- coding: utf-8 -*-import matplotlib.pyplot as pltimport numpy as npimport probscaleimport mathfrom scipy.stats import pearson3#plt.rcParams['font.sans-serif'] = ['Sarasa Gothic CL']plt.rcParams['font.sans-serif'] = 'Microsoft YaHei'def p3ProbScatter(data, yLogScale=False, probLeftLim=0.1, probRightLim=99.9, show=False):    """    该函数用于绘制水文数据的散点图    data                输入的原始水文数据，建议使用np.array类    yLogScale=False     y轴是否对数化显示，默认关闭    probLeftLim=0.1     x轴，概率起始点，默认为0.1%    probRightLim=99.9   x轴，概率结束点，默认为99.9%    show=False          是否立即显示，默认不显示    fig， ax            返回图形和坐标轴对象    """    # 创建单幅    fig, ax = plt.subplots()    # 自适应Y轴上限，最高位向上取整    MAX_Y = max(data)    MAX_Y = math.ceil(MAX_Y / 10 ** math.floor(math.log10(MAX_Y))                      ) * 10 ** math.floor(math.log10(MAX_Y))    # 设置y轴的范围，以及是否对数化显示    ax.set_ylim(0, MAX_Y)    if yLogScale == True:        ax.set_yscale('log')    # 设置x轴的范围，参数应该直接为百分数    # 比如50%，直接传入50    ax.set_xlim(probLeftLim, probRightLim)    ax.set_xscale('prob')    # 对数据进行降序排序，并计算经验频率    sortedData = np.sort(data)[::-1]    expProb = np.arange(1, len(sortedData) + 1) / (len(sortedData)+1) * 100    # 绘制散点图    plt.scatter(expProb, sortedData, label='经验概率点')    if show == True:        plt.show()    # 返回图像，轴    return fig, axdef p3Stat(data):    """    该函数用于统计数据的数学期望Ex，变差系数Cv，偏态系数Cs    """    # 数学期望    Ex = np.average(data)    # 模比系数    K = data / Ex    # 变差系数    Cv = np.sqrt(np.sum((K-1)**2) / (len(data) - 1))    # 偏态系数    Cs = np.sum((K - 1)**3)/((len(data) - 3) * Cv**3)    return Ex, Cv, Csdef p3ProbPlot(Ex, Cv, Cs, show=False):    """    该函数用于绘制P3曲线    Ex          曲线的数学期望    Cv          曲线的变差系数    Cs          曲线的偏态系数    show=False  是否显示图像，默认不显示    """    # 获取理论曲线的控制点    x = np.linspace(0.1, 99.9, 1000)    theoryY = (pearson3.ppf(1 - x / 100, Cs) * Cv + 1) * Ex    # 绘制理论曲线    plt.plot(x, theoryY, 'r-', lw=2, label='理论频率曲线')    plt.legend()    if show == True:        plt.show()# 测试例子data = np.array([778, 752, 709, 697, 650, 649, 644, 642, 623, 616, 612, 605, 575,                 573, 558, 551, 541, 526, 526, 521, 519, 510, 508, 506, 494, 486,                 480, 477, 472, 461, 442, 420, 415, 410, 398, 395, 382, 379, 376,                 370, 368, 359, 350, 349, 341, 337, 336, 315, 315, 292])statResult = p3Stat(data)print(statResult)p3ProbScatter(data)p3ProbPlot(statResult[0], statResult[1], 2 * statResult[1])plt.grid(True)plt.show()

将借鉴此版MATLAB此版P-III曲线绘图函数及自己以前的部分成果，最近计划利用MATLAB App Designer设计基于MATLAB App Designer的P-III曲线绘图系统，此版：仅为绘图系统，不包含经验点的计算；包含Ex、Cv、Cs的自由步长设置及微调；暂时不包含参数优化内容。目前已经取得一些进展。余下的工作将会在以后几个月的业余时间完成。

---