估計(jì)、適線法與Python實(shí)現(xiàn))
簡(jiǎn)介這份資源聚焦皮爾遜III型曲線的擬合與繪圖面向水文學(xué)、統(tǒng)計(jì)學(xué)及水資源管理領(lǐng)域的學(xué)習(xí)者與工程人員幫助解決流量、降雨等連續(xù)隨機(jī)變量的分布建模與頻率計(jì)算問(wèn)題。壓縮包共3個(gè)文件均為m腳本文件體積約3KB分別承擔(dān)曲線擬合、參數(shù)計(jì)算與繪圖插值等核心環(huán)節(jié)結(jié)構(gòu)精簡(jiǎn)、便于直接運(yùn)行與二次修改。內(nèi)容圍繞形狀參數(shù)κ、尺度參數(shù)θ和位置參數(shù)μ的估計(jì)展開(kāi)涵蓋最大似然估計(jì)思路、概率密度函數(shù)繪制以及通過(guò)插值獲取不同累積概率對(duì)應(yīng)流量值的方法可用于風(fēng)險(xiǎn)評(píng)估與防洪規(guī)劃等場(chǎng)景。目前已有1303人學(xué)習(xí)下載適合希望快速掌握皮爾遜III型曲線參數(shù)求解、圖形展示與頻率插值實(shí)操的讀者參考。1. 皮爾遜III型曲線擬合與繪圖從一組水位數(shù)據(jù)說(shuō)起手上有一組 30 年的年最大洪峰流量或者一批降雨極值領(lǐng)導(dǎo)要你給出百年一遇的設(shè)計(jì)值還要一張能放進(jìn)報(bào)告的頻率曲線圖。這時(shí)候繞不開(kāi)的就是皮爾遜III型曲線P-III 型曲線。它本質(zhì)上是伽馬分布的一個(gè)變種靠三個(gè)參數(shù)——均值、變差系數(shù) Cv、偏態(tài)系數(shù) Cs——把經(jīng)驗(yàn)點(diǎn)據(jù)擬合成一條光滑的理論頻率曲線再外延到稀遇頻率。標(biāo)題里那串doingni2、thusnwe看著像隨手敲的亂碼但核心訴求很清楚p3曲線的參數(shù)怎么定、皮爾遜曲線怎么擬合、圖怎么畫(huà)出來(lái)。這篇筆記就按我實(shí)際干活的順序把參數(shù)估計(jì)、擬合、繪圖、避坑一條線講透適合水文水利、可靠性分析、極值統(tǒng)計(jì)方向、手上正拿著 Excel 或 Python 準(zhǔn)備動(dòng)手的人。2. 皮爾遜III型曲線的三個(gè)參數(shù)到底怎么估2.1 先搞清楚 P-III 型曲線的數(shù)學(xué)形態(tài)皮爾遜III型分布的概率密度函數(shù)長(zhǎng)這樣f(x) [β^α / Γ(α)] * (x - a0)^(α-1) * exp[-β(x - a0)]其中 α 是形狀參數(shù)β 是尺度參數(shù)a0 是位置參數(shù)。水文上不直接用這三個(gè)而是換成均值 x?、變差系數(shù) Cv、偏態(tài)系數(shù) Cs換算關(guān)系是α 4 / Cs2β 2 / (x? · Cv · Cs)a0 x? · (1 - 2Cv / Cs)這三個(gè)參數(shù)里均值最穩(wěn)Cv 次之Cs 最玄學(xué)。為什么說(shuō) Cs 玄學(xué)因?yàn)?Cs 反映的是分布的不對(duì)稱(chēng)程度而實(shí)測(cè)樣本里極端值就那么幾個(gè)去掉一個(gè)最大值Cs 可能從 2.0 掉到 1.2。這就是為什么規(guī)范里對(duì) Cs 有各種經(jīng)驗(yàn)約束而不是讓你純靠矩法硬算。2.2 矩法、權(quán)函數(shù)法、線性矩法選哪個(gè)參數(shù)估計(jì)主流三條路我按實(shí)操推薦度排方法核心思路優(yōu)點(diǎn)缺點(diǎn)適用場(chǎng)景矩法用樣本矩直接算簡(jiǎn)單、快Cs 偏差大需大量樣本樣本 50 年粗算權(quán)函數(shù)法對(duì)矩法加權(quán)修正Cs 估計(jì)更穩(wěn)實(shí)現(xiàn)稍復(fù)雜樣本 30~50 年推薦線性矩法用線性組合代替冪矩對(duì)極端值不敏感需排序、公式多樣本含特大值首選我一般會(huì)均值用矩法Cv 用矩法Cs 先用線性矩法算再對(duì)照規(guī)范做適線調(diào)整。線性矩法對(duì)異常值穩(wěn)健這是血淚經(jīng)驗(yàn)——早年直接用矩法算 Cs一個(gè)特大洪水把 Cs 拉到 3.5適線時(shí)怎么調(diào)都別扭。2.3 用 Python 算三個(gè)參數(shù)的最小代碼import numpy as np from scipy import stats # 假設(shè) data 是年最大洪峰流量序列單位 m3/s data np.array([320, 450, 510, 380, 620, 290, 540, 470, 390, 580, 410, 360, 500, 440, 330, 610, 480, 420, 350, 560]) # 1. 均值 mean_x np.mean(data) # 2. 變差系數(shù) Cv 標(biāo)準(zhǔn)差 / 均值 cv np.std(data, ddof1) / mean_x # 3. 偏態(tài)系數(shù) Cs矩法 n len(data) cs_moment n * np.sum((data - mean_x)**3) / ((n-1)*(n-2)*np.std(data, ddof1)**3) # 4. 線性矩法估 Cs更穩(wěn)健 sorted_data np.sort(data) n len(sorted_data) # 線性矩l1均值, l2尺度, l3偏度 b0 np.mean(sorted_data) b1 np.sum([(i)/(n-1) * sorted_data[i] for i in range(1, n)]) / n b2 np.sum([(i*(i-1))/((n-1)*(n-2)) * sorted_data[i] for i in range(2, n)]) / n l1 b0 l2 2*b1 - b0 l3 6*b2 - 6*b1 b0 t3 l3 / l2 # L-skewness # 線性矩與 Cs 的近似換算水文常用 cs_lmom t3 * 3.0 # 粗略換算精確需查表或迭代 print(f均值 {mean_x:.2f}) print(fCv {cv:.3f}) print(fCs(矩法) {cs_moment:.3f}) print(fCs(線性矩) {cs_lmom:.3f})邏輯說(shuō)明先算均值和 Cv這兩個(gè)基本沒(méi)爭(zhēng)議。Cs 給了兩套矩法公式里(n-1)*(n-2)是無(wú)偏修正樣本小于 30 時(shí)這個(gè)修正很關(guān)鍵。線性矩部分b0/b1/b2 是概率權(quán)重矩t3 是 L-偏度乘 3 是水文上常用的粗略換算精確值要用多項(xiàng)式或迭代反解。參數(shù)說(shuō)明ddof1表示樣本標(biāo)準(zhǔn)差除以 n-1別用默認(rèn)的 0。線性矩里i從 1 開(kāi)始sorted_data[i]是第 i1 個(gè)排序值索引別搞錯(cuò)這是翻車(chē)高發(fā)區(qū)。3. 適線法讓理論曲線穿過(guò)經(jīng)驗(yàn)點(diǎn)據(jù)3.1 經(jīng)驗(yàn)頻率公式怎么選適線的前提是算出每個(gè)實(shí)測(cè)值對(duì)應(yīng)的經(jīng)驗(yàn)頻率。常用公式數(shù)學(xué)期望公式P m / (n1)海森公式P (m - 0.5) / n切哥達(dá)耶夫公式P (m - 0.3) / (n 0.4)我一般用數(shù)學(xué)期望公式因?yàn)樗鼰o(wú)偏且規(guī)范推薦。m 是降序排列的序號(hào)n 是樣本數(shù)。注意洪水、降雨極值用降序m1 對(duì)應(yīng)最大值的頻率最小別排反了。3.2 適線法的操作步驟適線不是讓程序自動(dòng)擬合就完事而是人機(jī)交互調(diào)參。標(biāo)準(zhǔn)流程用矩法/線性矩法算出初始 x?、Cv、Cs固定 x?微調(diào) Cv一般 ±0.05 步長(zhǎng)再調(diào) Cs一般 ±0.1 步長(zhǎng)觀察曲線頭部和尾部目測(cè)理論曲線與經(jīng)驗(yàn)點(diǎn)據(jù)的貼合度重點(diǎn)看頭部幾個(gè)大值和尾部外延段確定最終參數(shù)計(jì)算設(shè)計(jì)頻率對(duì)應(yīng)的設(shè)計(jì)值為什么固定 x?因?yàn)榫凳菍?shí)測(cè)樣本的客觀統(tǒng)計(jì)量適線時(shí)動(dòng)它沒(méi)有物理依據(jù)。Cv 控制曲線的陡緩Cs 控制曲線的彎曲方向。3.3 用 scipy 做適線并計(jì)算設(shè)計(jì)值import numpy as np from scipy import stats import matplotlib.pyplot as plt data np.array([320, 450, 510, 380, 620, 290, 540, 470, 390, 580, 410, 360, 500, 440, 330, 610, 480, 420, 350, 560]) n len(data) mean_x np.mean(data) cv np.std(data, ddof1) / mean_x cs 2.0 # 先給個(gè)初值后面手動(dòng)調(diào) # 經(jīng)驗(yàn)頻率數(shù)學(xué)期望公式 sorted_desc np.sort(data)[::-1] m np.arange(1, n1) p_emp m / (n 1) # P-III 型理論頻率scipy 的 pearson3 用偏度參數(shù) # scipy.stats.pearson3 的 skew 參數(shù)就是 Cs # 換算scipy 的 scale 標(biāo)準(zhǔn)差loc 均值 std_x mean_x * cv p_theory np.linspace(0.001, 0.999, 500) # pearson3.ppf 輸入累積概率輸出分位數(shù) q_theory stats.pearson3.ppf(1 - p_theory, skewcs, locmean_x, scalestd_x) # 計(jì)算設(shè)計(jì)值P1% 百年一遇 q_1pct stats.pearson3.ppf(1 - 0.01, skewcs, locmean_x, scalestd_x) q_01pct stats.pearson3.ppf(1 - 0.001, skewcs, locmean_x, scalestd_x) print(fP1% 設(shè)計(jì)值 {q_1pct:.1f}) print(fP0.1% 設(shè)計(jì)值 {q_01pct:.1f})邏輯說(shuō)明pearson3.ppf輸入的是累積概率所以用1 - p把頻率轉(zhuǎn)成累積概率。skewcs直接對(duì)應(yīng)偏態(tài)系數(shù)scale是標(biāo)準(zhǔn)差不是 Cv這里容易搞混。設(shè)計(jì)值就是查 ppf。參數(shù)說(shuō)明cs初值給 2.0 是水文常見(jiàn)范圍1.5~3.0實(shí)際要調(diào)。p_theory從 0.001 到 0.999 是為了畫(huà)圖時(shí)曲線覆蓋全范圍。注意ppf在 p 接近 0 或 1 時(shí)可能返回 inf畫(huà)圖時(shí)要截?cái)唷?. 繪圖把頻率曲線畫(huà)成能進(jìn)報(bào)告的樣子4.1 海森概率格紙是標(biāo)配水文頻率曲線不能畫(huà)在普通坐標(biāo)紙上要用海森概率格紙——縱軸是流量線性橫軸是頻率正態(tài)分位刻度。這樣 P-III 型曲線在圖上接近直線方便目測(cè)適線效果。用 matplotlib 實(shí)現(xiàn)import numpy as np import matplotlib.pyplot as plt from scipy import stats data np.array([320, 450, 510, 380, 620, 290, 540, 470, 390, 580, 410, 360, 500, 440, 330, 610, 480, 420, 350, 560]) n len(data) mean_x np.mean(data) cv np.std(data, ddof1) / mean_x cs 2.0 std_x mean_x * cv # 經(jīng)驗(yàn)點(diǎn) sorted_desc np.sort(data)[::-1] m np.arange(1, n1) p_emp m / (n 1) # 理論曲線 p_theory np.linspace(0.001, 0.999, 1000) q_theory stats.pearson3.ppf(1 - p_theory, skewcs, locmean_x, scalestd_x) fig, ax plt.subplots(figsize(10, 6)) # 橫軸用正態(tài)分位刻度把頻率轉(zhuǎn)成標(biāo)準(zhǔn)正態(tài)分位數(shù) def freq_to_z(p): return stats.norm.ppf(1 - p) # 畫(huà)理論曲線 z_theory freq_to_z(p_theory) ax.plot(z_theory, q_theory, b-, linewidth2, labelfP-III (Cs{cs})) # 畫(huà)經(jīng)驗(yàn)點(diǎn) z_emp freq_to_z(p_emp) ax.scatter(z_emp, sorted_desc, colorred, zorder5, label經(jīng)驗(yàn)點(diǎn)據(jù)) # 設(shè)置橫軸刻度標(biāo)簽為頻率值 ticks_p [0.01, 0.1, 0.5, 1, 2, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 99] ticks_z [freq_to_z(p/100) for p in ticks_p] ax.set_xticks(ticks_z) ax.set_xticklabels([f{p} for p in ticks_p]) ax.set_xlabel(頻率 P (%)) ax.set_ylabel(流量 (m3/s)) ax.set_title(皮爾遜III型頻率曲線) ax.grid(True, linestyle--, alpha0.5) ax.legend() plt.tight_layout() plt.savefig(p3_curve.png, dpi300) plt.show()邏輯說(shuō)明核心技巧是freq_to_z把頻率轉(zhuǎn)成標(biāo)準(zhǔn)正態(tài)分位數(shù)這樣橫軸就變成了概率格紙的刻度。理論曲線用plot連線經(jīng)驗(yàn)點(diǎn)用scatter。ticks_p是你要顯示的頻率刻度轉(zhuǎn)成 z 值后設(shè)置xticklabels。參數(shù)說(shuō)明dpi300是報(bào)告插圖的最低要求低了打印發(fā)虛。zorder5讓經(jīng)驗(yàn)點(diǎn)蓋在曲線上。linspace(0.001, 0.999, 1000)保證曲線兩端平滑點(diǎn)數(shù)少了會(huì)折線感明顯。4.2 適線調(diào)整的目測(cè)要點(diǎn)圖畫(huà)出來(lái)后重點(diǎn)看三處頭部P 5%理論曲線是否穿過(guò)最大的幾個(gè)經(jīng)驗(yàn)點(diǎn)。如果曲線明顯低于最大點(diǎn)說(shuō)明 Cs 偏小或 Cv 偏小。中部P 20%~80%曲線應(yīng)穿過(guò)點(diǎn)群中心偏離說(shuō)明均值或 Cv 有問(wèn)題。尾部P 90%外延段是否合理不能出現(xiàn)負(fù)值或急劇上翹。我一般會(huì)先調(diào) Cv 讓曲線整體陡緩合適再調(diào) Cs 讓頭部貼合。如果怎么調(diào)都不順回頭檢查經(jīng)驗(yàn)頻率公式和數(shù)據(jù)里有沒(méi)有混入非極值樣本。5. 避坑皮爾遜III型擬合里最容易翻車(chē)的 5 個(gè)地方5.1 現(xiàn)象Cs 算出來(lái)是負(fù)數(shù)曲線反向彎曲原因樣本里小值偏多或存在系統(tǒng)性偏小導(dǎo)致偏度統(tǒng)計(jì)量為負(fù)。但水文極值分布理論上 Cs 應(yīng)大于 0。解決檢查數(shù)據(jù)是否混入了非年最大值比如把月最大值當(dāng)年最大值。如果數(shù)據(jù)沒(méi)問(wèn)題用線性矩法重算或按規(guī)范取 Cs 2~3 倍 Cv 作為約束。別硬用負(fù) Cs 去適線外延會(huì)出負(fù)值。5.2 現(xiàn)象設(shè)計(jì)值 P0.1% 算出來(lái)比 P1% 還小原因pearson3.ppf的累積概率方向搞反了或者頻率和累積概率換算錯(cuò)誤。解決記住P頻率越小設(shè)計(jì)值越大。ppf(1-p)里 p 是頻率1-p 是累積概率。P0.001 對(duì)應(yīng)累積概率 0.999分位數(shù)最大。檢查代碼里有沒(méi)有寫(xiě)成ppf(p)。5.3 現(xiàn)象經(jīng)驗(yàn)點(diǎn)畫(huà)在圖上位置全亂和理論曲線對(duì)不上原因排序方向錯(cuò)了。頻率最小的點(diǎn)應(yīng)該對(duì)應(yīng)最大值如果升序排列m1 對(duì)應(yīng)最小值頻率和值就反了。解決用np.sort(data)[::-1]降序排列或者np.sort(data)[::-1]后 m 從 1 開(kāi)始。畫(huà)圖前打印前 5 個(gè)(p_emp, sorted_desc)對(duì)確認(rèn) P 小值大。5.4 現(xiàn)象matplotlib 畫(huà)出來(lái)橫軸刻度擠在一起頻率標(biāo)簽重疊原因概率格紙的刻度是非線性的兩端密中間疏默認(rèn)刻度會(huì)重疊。解決手動(dòng)指定ticks_p只保留關(guān)鍵頻率點(diǎn)0.01, 0.1, 0.5, 1, 2, 5, 10, 20, 50, 80, 90, 95, 99用set_xticksset_xticklabels控制。必要時(shí)旋轉(zhuǎn)標(biāo)簽rotation45。5.5 現(xiàn)象換一組數(shù)據(jù)后 Cs 從 2.0 跳到 3.5設(shè)計(jì)值差一倍原因樣本容量太小 30Cs 對(duì)極端值極度敏感這是 P-III 型曲線的固有缺陷。解決小樣本時(shí)用地區(qū)綜合法——把鄰近站點(diǎn)的 Cs 綜合成一個(gè)區(qū)域值或者用線性矩法降低敏感度。如果必須單站適線在報(bào)告里注明 Cs 的不確定性范圍給出設(shè)計(jì)值的區(qū)間而不是單點(diǎn)。6. 進(jìn)階用線性矩法自動(dòng)適線并做不確定性評(píng)估手動(dòng)調(diào)參適合出報(bào)告但如果你要批量處理幾十個(gè)站點(diǎn)或者想量化設(shè)計(jì)值的不確定性得換自動(dòng)化思路。我現(xiàn)在的習(xí)慣是線性矩法估初值 優(yōu)化算法微調(diào) Bootstrap 重采樣給置信區(qū)間。線性矩法的優(yōu)勢(shì)前面提過(guò)對(duì)極端值穩(wěn)健。自動(dòng)適線的目標(biāo)函數(shù)用經(jīng)驗(yàn)點(diǎn)與理論曲線的加權(quán)殘差平方和權(quán)重取經(jīng)驗(yàn)頻率的倒數(shù)讓頭部大值權(quán)重更高import numpy as np from scipy import stats, optimize data np.array([320, 450, 510, 380, 620, 290, 540, 470, 390, 580, 410, 360, 500, 440, 330, 610, 480, 420, 350, 560]) n len(data) sorted_desc np.sort(data)[::-1] m np.arange(1, n1) p_emp m / (n 1) def objective(params): mean_x, cv, cs params if cv 0 or cs 0: return 1e10 std_x mean_x * cv q_theory stats.pearson3.ppf(1 - p_emp, skewcs, locmean_x, scalestd_x) # 加權(quán)殘差頭部權(quán)重高 weights 1.0 / p_emp return np.sum(weights * (sorted_desc - q_theory)**2) # 初值矩法 mean0 np.mean(data) cv0 np.std(data, ddof1) / mean0 cs0 2.0 result optimize.minimize(objective, [mean0, cv0, cs0], methodNelder-Mead, options{xatol: 1e-4, fatol: 1e-4}) mean_fit, cv_fit, cs_fit result.x print(f適線結(jié)果: 均值{mean_fit:.1f}, Cv{cv_fit:.3f}, Cs{cs_fit:.3f}) # Bootstrap 不確定性評(píng)估 n_boot 500 q_1pct_boot [] for _ in range(n_boot): sample np.random.choice(data, sizen, replaceTrue) s_desc np.sort(sample)[::-1] p_b np.arange(1, n1) / (n 1) def obj_b(params): mx, cv, cs params if cv 0 or cs 0: return 1e10 q_t stats.pearson3.ppf(1 - p_b, skewcs, locmx, scalemx*cv) return np.sum((1.0/p_b) * (s_desc - q_t)**2) res_b optimize.minimize(obj_b, [np.mean(sample), np.std(sample, ddof1)/np.mean(sample), 2.0], methodNelder-Mead) mx, cv, cs res_b.x q_1pct_boot.append(stats.pearson3.ppf(0.99, skewcs, locmx, scalemx*cv)) q_1pct_boot np.array(q_1pct_boot) print(fP1% 設(shè)計(jì)值 90% 置信區(qū)間: [{np.percentile(q_1pct_boot, 5):.1f}, {np.percentile(q_1pct_boot, 95):.1f}])邏輯說(shuō)明objective是加權(quán)殘差平方和權(quán)重1/p_emp讓 P 小的點(diǎn)大值在擬合中占更大比重這符合水文關(guān)注稀遇事件的邏輯。Nelder-Mead不需要梯度適合這種參數(shù)有約束的優(yōu)化。Bootstrap 重采樣 500 次每次重新適線得到設(shè)計(jì)值的分布取 5% 和 95% 分位數(shù)作為 90% 置信區(qū)間。參數(shù)說(shuō)明xatol和fatol控制收斂精度1e-4 對(duì)工程夠用。n_boot500是精度和耗時(shí)的折中要求高可以上 1000。Bootstrap 里每次都要重新優(yōu)化計(jì)算量是 500 倍樣本大時(shí)要有耐心。這套流程跑下來(lái)你不僅能給出設(shè)計(jì)值還能告訴領(lǐng)導(dǎo)「這個(gè)值有 90% 把握落在某個(gè)區(qū)間內(nèi)」比單點(diǎn)值有說(shuō)服力得多。我現(xiàn)在的習(xí)慣是任何要簽字的設(shè)計(jì)值都附一個(gè) Bootstrap 置信區(qū)間這是給自己留的后悔藥——萬(wàn)一后面數(shù)據(jù)更新至少能說(shuō)清楚當(dāng)時(shí)的不確定性范圍。皮爾遜III型曲線擬合不難難的是知道它的邊界在哪希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取 皮爾遜III型曲線擬合與繪圖：參數(shù)估計(jì)、適線法與Python實(shí)現(xiàn)
簡(jiǎn)介這份資源聚焦皮爾遜III型曲線的擬合與繪圖面向水文學(xué)、統(tǒng)計(jì)學(xué)及水資源管理領(lǐng)域的學(xué)習(xí)者與工程人員幫助解決流量、降雨等連續(xù)隨機(jī)變量的分布建模與頻率計(jì)算問(wèn)題。壓縮包共3個(gè)文件均為m腳本文件體積約3KB分別承擔(dān)曲線擬合、參數(shù)計(jì)算與繪圖插值等核心環(huán)節(jié)結(jié)構(gòu)精簡(jiǎn)、便于直接運(yùn)行與二次修改。內(nèi)容圍繞形狀參數(shù)κ、尺度參數(shù)θ和位置參數(shù)μ的估計(jì)展開(kāi)涵蓋最大似然估計(jì)思路、概率密度函數(shù)繪制以及通過(guò)插值獲取不同累積概率對(duì)應(yīng)流量值的方法可用于風(fēng)險(xiǎn)評(píng)估與防洪規(guī)劃等場(chǎng)景。目前已有1303人學(xué)習(xí)下載適合希望快速掌握皮爾遜III型曲線參數(shù)求解、圖形展示與頻率插值實(shí)操的讀者參考。1. 皮爾遜III型曲線擬合與繪圖從一組水位數(shù)據(jù)說(shuō)起手上有一組 30 年的年最大洪峰流量或者一批降雨極值領(lǐng)導(dǎo)要你給出百年一遇的設(shè)計(jì)值還要一張能放進(jìn)報(bào)告的頻率曲線圖。這時(shí)候繞不開(kāi)的就是皮爾遜III型曲線P-III 型曲線。它本質(zhì)上是伽馬分布的一個(gè)變種靠三個(gè)參數(shù)——均值、變差系數(shù) Cv、偏態(tài)系數(shù) Cs——把經(jīng)驗(yàn)點(diǎn)據(jù)擬合成一條光滑的理論頻率曲線再外延到稀遇頻率。標(biāo)題里那串doingni2、thusnwe看著像隨手敲的亂碼但核心訴求很清楚p3曲線的參數(shù)怎么定、皮爾遜曲線怎么擬合、圖怎么畫(huà)出來(lái)。這篇筆記就按我實(shí)際干活的順序把參數(shù)估計(jì)、擬合、繪圖、避坑一條線講透適合水文水利、可靠性分析、極值統(tǒng)計(jì)方向、手上正拿著 Excel 或 Python 準(zhǔn)備動(dòng)手的人。2. 皮爾遜III型曲線的三個(gè)參數(shù)到底怎么估2.1 先搞清楚 P-III 型曲線的數(shù)學(xué)形態(tài)皮爾遜III型分布的概率密度函數(shù)長(zhǎng)這樣f(x) [β^α / Γ(α)] * (x - a0)^(α-1) * exp[-β(x - a0)]其中 α 是形狀參數(shù)β 是尺度參數(shù)a0 是位置參數(shù)。水文上不直接用這三個(gè)而是換成均值 x?、變差系數(shù) Cv、偏態(tài)系數(shù) Cs換算關(guān)系是α 4 / Cs2β 2 / (x? · Cv · Cs)a0 x? · (1 - 2Cv / Cs)這三個(gè)參數(shù)里均值最穩(wěn)Cv 次之Cs 最玄學(xué)。為什么說(shuō) Cs 玄學(xué)因?yàn)?Cs 反映的是分布的不對(duì)稱(chēng)程度而實(shí)測(cè)樣本里極端值就那么幾個(gè)去掉一個(gè)最大值Cs 可能從 2.0 掉到 1.2。這就是為什么規(guī)范里對(duì) Cs 有各種經(jīng)驗(yàn)約束而不是讓你純靠矩法硬算。2.2 矩法、權(quán)函數(shù)法、線性矩法選哪個(gè)參數(shù)估計(jì)主流三條路我按實(shí)操推薦度排方法核心思路優(yōu)點(diǎn)缺點(diǎn)適用場(chǎng)景矩法用樣本矩直接算簡(jiǎn)單、快Cs 偏差大需大量樣本樣本 50 年粗算權(quán)函數(shù)法對(duì)矩法加權(quán)修正Cs 估計(jì)更穩(wěn)實(shí)現(xiàn)稍復(fù)雜樣本 30~50 年推薦線性矩法用線性組合代替冪矩對(duì)極端值不敏感需排序、公式多樣本含特大值首選我一般會(huì)均值用矩法Cv 用矩法Cs 先用線性矩法算再對(duì)照規(guī)范做適線調(diào)整。線性矩法對(duì)異常值穩(wěn)健這是血淚經(jīng)驗(yàn)——早年直接用矩法算 Cs一個(gè)特大洪水把 Cs 拉到 3.5適線時(shí)怎么調(diào)都別扭。2.3 用 Python 算三個(gè)參數(shù)的最小代碼import numpy as np from scipy import stats # 假設(shè) data 是年最大洪峰流量序列單位 m3/s data np.array([320, 450, 510, 380, 620, 290, 540, 470, 390, 580, 410, 360, 500, 440, 330, 610, 480, 420, 350, 560]) # 1. 均值 mean_x np.mean(data) # 2. 變差系數(shù) Cv 標(biāo)準(zhǔn)差 / 均值 cv np.std(data, ddof1) / mean_x # 3. 偏態(tài)系數(shù) Cs矩法 n len(data) cs_moment n * np.sum((data - mean_x)**3) / ((n-1)*(n-2)*np.std(data, ddof1)**3) # 4. 線性矩法估 Cs更穩(wěn)健 sorted_data np.sort(data) n len(sorted_data) # 線性矩l1均值, l2尺度, l3偏度 b0 np.mean(sorted_data) b1 np.sum([(i)/(n-1) * sorted_data[i] for i in range(1, n)]) / n b2 np.sum([(i*(i-1))/((n-1)*(n-2)) * sorted_data[i] for i in range(2, n)]) / n l1 b0 l2 2*b1 - b0 l3 6*b2 - 6*b1 b0 t3 l3 / l2 # L-skewness # 線性矩與 Cs 的近似換算水文常用 cs_lmom t3 * 3.0 # 粗略換算精確需查表或迭代 print(f均值 {mean_x:.2f}) print(fCv {cv:.3f}) print(fCs(矩法) {cs_moment:.3f}) print(fCs(線性矩) {cs_lmom:.3f})邏輯說(shuō)明先算均值和 Cv這兩個(gè)基本沒(méi)爭(zhēng)議。Cs 給了兩套矩法公式里(n-1)*(n-2)是無(wú)偏修正樣本小于 30 時(shí)這個(gè)修正很關(guān)鍵。線性矩部分b0/b1/b2 是概率權(quán)重矩t3 是 L-偏度乘 3 是水文上常用的粗略換算精確值要用多項(xiàng)式或迭代反解。參數(shù)說(shuō)明ddof1表示樣本標(biāo)準(zhǔn)差除以 n-1別用默認(rèn)的 0。線性矩里i從 1 開(kāi)始sorted_data[i]是第 i1 個(gè)排序值索引別搞錯(cuò)這是翻車(chē)高發(fā)區(qū)。3. 適線法讓理論曲線穿過(guò)經(jīng)驗(yàn)點(diǎn)據(jù)3.1 經(jīng)驗(yàn)頻率公式怎么選適線的前提是算出每個(gè)實(shí)測(cè)值對(duì)應(yīng)的經(jīng)驗(yàn)頻率。常用公式數(shù)學(xué)期望公式P m / (n1)海森公式P (m - 0.5) / n切哥達(dá)耶夫公式P (m - 0.3) / (n 0.4)我一般用數(shù)學(xué)期望公式因?yàn)樗鼰o(wú)偏且規(guī)范推薦。m 是降序排列的序號(hào)n 是樣本數(shù)。注意洪水、降雨極值用降序m1 對(duì)應(yīng)最大值的頻率最小別排反了。3.2 適線法的操作步驟適線不是讓程序自動(dòng)擬合就完事而是人機(jī)交互調(diào)參。標(biāo)準(zhǔn)流程用矩法/線性矩法算出初始 x?、Cv、Cs固定 x?微調(diào) Cv一般 ±0.05 步長(zhǎng)再調(diào) Cs一般 ±0.1 步長(zhǎng)觀察曲線頭部和尾部目測(cè)理論曲線與經(jīng)驗(yàn)點(diǎn)據(jù)的貼合度重點(diǎn)看頭部幾個(gè)大值和尾部外延段確定最終參數(shù)計(jì)算設(shè)計(jì)頻率對(duì)應(yīng)的設(shè)計(jì)值為什么固定 x?因?yàn)榫凳菍?shí)測(cè)樣本的客觀統(tǒng)計(jì)量適線時(shí)動(dòng)它沒(méi)有物理依據(jù)。Cv 控制曲線的陡緩Cs 控制曲線的彎曲方向。3.3 用 scipy 做適線并計(jì)算設(shè)計(jì)值import numpy as np from scipy import stats import matplotlib.pyplot as plt data np.array([320, 450, 510, 380, 620, 290, 540, 470, 390, 580, 410, 360, 500, 440, 330, 610, 480, 420, 350, 560]) n len(data) mean_x np.mean(data) cv np.std(data, ddof1) / mean_x cs 2.0 # 先給個(gè)初值后面手動(dòng)調(diào) # 經(jīng)驗(yàn)頻率數(shù)學(xué)期望公式 sorted_desc np.sort(data)[::-1] m np.arange(1, n1) p_emp m / (n 1) # P-III 型理論頻率scipy 的 pearson3 用偏度參數(shù) # scipy.stats.pearson3 的 skew 參數(shù)就是 Cs # 換算scipy 的 scale 標(biāo)準(zhǔn)差loc 均值 std_x mean_x * cv p_theory np.linspace(0.001, 0.999, 500) # pearson3.ppf 輸入累積概率輸出分位數(shù) q_theory stats.pearson3.ppf(1 - p_theory, skewcs, locmean_x, scalestd_x) # 計(jì)算設(shè)計(jì)值P1% 百年一遇 q_1pct stats.pearson3.ppf(1 - 0.01, skewcs, locmean_x, scalestd_x) q_01pct stats.pearson3.ppf(1 - 0.001, skewcs, locmean_x, scalestd_x) print(fP1% 設(shè)計(jì)值 {q_1pct:.1f}) print(fP0.1% 設(shè)計(jì)值 {q_01pct:.1f})邏輯說(shuō)明pearson3.ppf輸入的是累積概率所以用1 - p把頻率轉(zhuǎn)成累積概率。skewcs直接對(duì)應(yīng)偏態(tài)系數(shù)scale是標(biāo)準(zhǔn)差不是 Cv這里容易搞混。設(shè)計(jì)值就是查 ppf。參數(shù)說(shuō)明cs初值給 2.0 是水文常見(jiàn)范圍1.5~3.0實(shí)際要調(diào)。p_theory從 0.001 到 0.999 是為了畫(huà)圖時(shí)曲線覆蓋全范圍。注意ppf在 p 接近 0 或 1 時(shí)可能返回 inf畫(huà)圖時(shí)要截?cái)唷?. 繪圖把頻率曲線畫(huà)成能進(jìn)報(bào)告的樣子4.1 海森概率格紙是標(biāo)配水文頻率曲線不能畫(huà)在普通坐標(biāo)紙上要用海森概率格紙——縱軸是流量線性橫軸是頻率正態(tài)分位刻度。這樣 P-III 型曲線在圖上接近直線方便目測(cè)適線效果。用 matplotlib 實(shí)現(xiàn)import numpy as np import matplotlib.pyplot as plt from scipy import stats data np.array([320, 450, 510, 380, 620, 290, 540, 470, 390, 580, 410, 360, 500, 440, 330, 610, 480, 420, 350, 560]) n len(data) mean_x np.mean(data) cv np.std(data, ddof1) / mean_x cs 2.0 std_x mean_x * cv # 經(jīng)驗(yàn)點(diǎn) sorted_desc np.sort(data)[::-1] m np.arange(1, n1) p_emp m / (n 1) # 理論曲線 p_theory np.linspace(0.001, 0.999, 1000) q_theory stats.pearson3.ppf(1 - p_theory, skewcs, locmean_x, scalestd_x) fig, ax plt.subplots(figsize(10, 6)) # 橫軸用正態(tài)分位刻度把頻率轉(zhuǎn)成標(biāo)準(zhǔn)正態(tài)分位數(shù) def freq_to_z(p): return stats.norm.ppf(1 - p) # 畫(huà)理論曲線 z_theory freq_to_z(p_theory) ax.plot(z_theory, q_theory, b-, linewidth2, labelfP-III (Cs{cs})) # 畫(huà)經(jīng)驗(yàn)點(diǎn) z_emp freq_to_z(p_emp) ax.scatter(z_emp, sorted_desc, colorred, zorder5, label經(jīng)驗(yàn)點(diǎn)據(jù)) # 設(shè)置橫軸刻度標(biāo)簽為頻率值 ticks_p [0.01, 0.1, 0.5, 1, 2, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 99] ticks_z [freq_to_z(p/100) for p in ticks_p] ax.set_xticks(ticks_z) ax.set_xticklabels([f{p} for p in ticks_p]) ax.set_xlabel(頻率 P (%)) ax.set_ylabel(流量 (m3/s)) ax.set_title(皮爾遜III型頻率曲線) ax.grid(True, linestyle--, alpha0.5) ax.legend() plt.tight_layout() plt.savefig(p3_curve.png, dpi300) plt.show()邏輯說(shuō)明核心技巧是freq_to_z把頻率轉(zhuǎn)成標(biāo)準(zhǔn)正態(tài)分位數(shù)這樣橫軸就變成了概率格紙的刻度。理論曲線用plot連線經(jīng)驗(yàn)點(diǎn)用scatter。ticks_p是你要顯示的頻率刻度轉(zhuǎn)成 z 值后設(shè)置xticklabels。參數(shù)說(shuō)明dpi300是報(bào)告插圖的最低要求低了打印發(fā)虛。zorder5讓經(jīng)驗(yàn)點(diǎn)蓋在曲線上。linspace(0.001, 0.999, 1000)保證曲線兩端平滑點(diǎn)數(shù)少了會(huì)折線感明顯。4.2 適線調(diào)整的目測(cè)要點(diǎn)圖畫(huà)出來(lái)后重點(diǎn)看三處頭部P 5%理論曲線是否穿過(guò)最大的幾個(gè)經(jīng)驗(yàn)點(diǎn)。如果曲線明顯低于最大點(diǎn)說(shuō)明 Cs 偏小或 Cv 偏小。中部P 20%~80%曲線應(yīng)穿過(guò)點(diǎn)群中心偏離說(shuō)明均值或 Cv 有問(wèn)題。尾部P 90%外延段是否合理不能出現(xiàn)負(fù)值或急劇上翹。我一般會(huì)先調(diào) Cv 讓曲線整體陡緩合適再調(diào) Cs 讓頭部貼合。如果怎么調(diào)都不順回頭檢查經(jīng)驗(yàn)頻率公式和數(shù)據(jù)里有沒(méi)有混入非極值樣本。5. 避坑皮爾遜III型擬合里最容易翻車(chē)的 5 個(gè)地方5.1 現(xiàn)象Cs 算出來(lái)是負(fù)數(shù)曲線反向彎曲原因樣本里小值偏多或存在系統(tǒng)性偏小導(dǎo)致偏度統(tǒng)計(jì)量為負(fù)。但水文極值分布理論上 Cs 應(yīng)大于 0。解決檢查數(shù)據(jù)是否混入了非年最大值比如把月最大值當(dāng)年最大值。如果數(shù)據(jù)沒(méi)問(wèn)題用線性矩法重算或按規(guī)范取 Cs 2~3 倍 Cv 作為約束。別硬用負(fù) Cs 去適線外延會(huì)出負(fù)值。5.2 現(xiàn)象設(shè)計(jì)值 P0.1% 算出來(lái)比 P1% 還小原因pearson3.ppf的累積概率方向搞反了或者頻率和累積概率換算錯(cuò)誤。解決記住P頻率越小設(shè)計(jì)值越大。ppf(1-p)里 p 是頻率1-p 是累積概率。P0.001 對(duì)應(yīng)累積概率 0.999分位數(shù)最大。檢查代碼里有沒(méi)有寫(xiě)成ppf(p)。5.3 現(xiàn)象經(jīng)驗(yàn)點(diǎn)畫(huà)在圖上位置全亂和理論曲線對(duì)不上原因排序方向錯(cuò)了。頻率最小的點(diǎn)應(yīng)該對(duì)應(yīng)最大值如果升序排列m1 對(duì)應(yīng)最小值頻率和值就反了。解決用np.sort(data)[::-1]降序排列或者np.sort(data)[::-1]后 m 從 1 開(kāi)始。畫(huà)圖前打印前 5 個(gè)(p_emp, sorted_desc)對(duì)確認(rèn) P 小值大。5.4 現(xiàn)象matplotlib 畫(huà)出來(lái)橫軸刻度擠在一起頻率標(biāo)簽重疊原因概率格紙的刻度是非線性的兩端密中間疏默認(rèn)刻度會(huì)重疊。解決手動(dòng)指定ticks_p只保留關(guān)鍵頻率點(diǎn)0.01, 0.1, 0.5, 1, 2, 5, 10, 20, 50, 80, 90, 95, 99用set_xticksset_xticklabels控制。必要時(shí)旋轉(zhuǎn)標(biāo)簽rotation45。5.5 現(xiàn)象換一組數(shù)據(jù)后 Cs 從 2.0 跳到 3.5設(shè)計(jì)值差一倍原因樣本容量太小 30Cs 對(duì)極端值極度敏感這是 P-III 型曲線的固有缺陷。解決小樣本時(shí)用地區(qū)綜合法——把鄰近站點(diǎn)的 Cs 綜合成一個(gè)區(qū)域值或者用線性矩法降低敏感度。如果必須單站適線在報(bào)告里注明 Cs 的不確定性范圍給出設(shè)計(jì)值的區(qū)間而不是單點(diǎn)。6. 進(jìn)階用線性矩法自動(dòng)適線并做不確定性評(píng)估手動(dòng)調(diào)參適合出報(bào)告但如果你要批量處理幾十個(gè)站點(diǎn)或者想量化設(shè)計(jì)值的不確定性得換自動(dòng)化思路。我現(xiàn)在的習(xí)慣是線性矩法估初值 優(yōu)化算法微調(diào) Bootstrap 重采樣給置信區(qū)間。線性矩法的優(yōu)勢(shì)前面提過(guò)對(duì)極端值穩(wěn)健。自動(dòng)適線的目標(biāo)函數(shù)用經(jīng)驗(yàn)點(diǎn)與理論曲線的加權(quán)殘差平方和權(quán)重取經(jīng)驗(yàn)頻率的倒數(shù)讓頭部大值權(quán)重更高import numpy as np from scipy import stats, optimize data np.array([320, 450, 510, 380, 620, 290, 540, 470, 390, 580, 410, 360, 500, 440, 330, 610, 480, 420, 350, 560]) n len(data) sorted_desc np.sort(data)[::-1] m np.arange(1, n1) p_emp m / (n 1) def objective(params): mean_x, cv, cs params if cv 0 or cs 0: return 1e10 std_x mean_x * cv q_theory stats.pearson3.ppf(1 - p_emp, skewcs, locmean_x, scalestd_x) # 加權(quán)殘差頭部權(quán)重高 weights 1.0 / p_emp return np.sum(weights * (sorted_desc - q_theory)**2) # 初值矩法 mean0 np.mean(data) cv0 np.std(data, ddof1) / mean0 cs0 2.0 result optimize.minimize(objective, [mean0, cv0, cs0], methodNelder-Mead, options{xatol: 1e-4, fatol: 1e-4}) mean_fit, cv_fit, cs_fit result.x print(f適線結(jié)果: 均值{mean_fit:.1f}, Cv{cv_fit:.3f}, Cs{cs_fit:.3f}) # Bootstrap 不確定性評(píng)估 n_boot 500 q_1pct_boot [] for _ in range(n_boot): sample np.random.choice(data, sizen, replaceTrue) s_desc np.sort(sample)[::-1] p_b np.arange(1, n1) / (n 1) def obj_b(params): mx, cv, cs params if cv 0 or cs 0: return 1e10 q_t stats.pearson3.ppf(1 - p_b, skewcs, locmx, scalemx*cv) return np.sum((1.0/p_b) * (s_desc - q_t)**2) res_b optimize.minimize(obj_b, [np.mean(sample), np.std(sample, ddof1)/np.mean(sample), 2.0], methodNelder-Mead) mx, cv, cs res_b.x q_1pct_boot.append(stats.pearson3.ppf(0.99, skewcs, locmx, scalemx*cv)) q_1pct_boot np.array(q_1pct_boot) print(fP1% 設(shè)計(jì)值 90% 置信區(qū)間: [{np.percentile(q_1pct_boot, 5):.1f}, {np.percentile(q_1pct_boot, 95):.1f}])邏輯說(shuō)明objective是加權(quán)殘差平方和權(quán)重1/p_emp讓 P 小的點(diǎn)大值在擬合中占更大比重這符合水文關(guān)注稀遇事件的邏輯。Nelder-Mead不需要梯度適合這種參數(shù)有約束的優(yōu)化。Bootstrap 重采樣 500 次每次重新適線得到設(shè)計(jì)值的分布取 5% 和 95% 分位數(shù)作為 90% 置信區(qū)間。參數(shù)說(shuō)明xatol和fatol控制收斂精度1e-4 對(duì)工程夠用。n_boot500是精度和耗時(shí)的折中要求高可以上 1000。Bootstrap 里每次都要重新優(yōu)化計(jì)算量是 500 倍樣本大時(shí)要有耐心。這套流程跑下來(lái)你不僅能給出設(shè)計(jì)值還能告訴領(lǐng)導(dǎo)「這個(gè)值有 90% 把握落在某個(gè)區(qū)間內(nèi)」比單點(diǎn)值有說(shuō)服力得多。我現(xiàn)在的習(xí)慣是任何要簽字的設(shè)計(jì)值都附一個(gè) Bootstrap 置信區(qū)間這是給自己留的后悔藥——萬(wàn)一后面數(shù)據(jù)更新至少能說(shuō)清楚當(dāng)時(shí)的不確定性范圍。皮爾遜III型曲線擬合不難難的是知道它的邊界在哪希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取