極值Ⅰ型分布設計安全修正值系數的研究

王曉惠，潘曉春，沈旭偉

（中國能源建設集團江蘇省電力設計院有限公司，江蘇 南京 211102）

根據現行《港口與航道水文規范》，海港工程極端高水位、極端低水位、設計波高及波周期等要素的分析計算可分別采用極值Ⅰ型、Pearson-Ⅲ型概率分布模型，也可以根據實測資料擬合最佳的原則，在兩個分布模型之間擇優選取[1]。此外，包括海港工程在內的結構抗風設計風速統計，也常采用兩個分布模型進行分析計算[2]。

我國設計洪水計算推薦采用Pearson-Ⅲ型概率分布模型，并以設計值的均方誤差作為安全修正值來估算抽樣誤差，通過設計值加上對應頻率安全修正值，以確保計算成果的可靠[3]。安全修正值由安全修正值系數B 值、樣本標準差及樣本容量綜合確定，同時與概率分布模型的參數估計方法密切相關[4-10]。目前關于Pearson-Ⅲ型概率分布模型的安全修正值系數B 值已有較多的研究[4-10]，有代表性的成果已被我國規范采用[3]，并應用到工程實踐中。然而，現行相關技術標準[1-2]在推薦采用極值Ⅰ型分布模型進行頻率分析計算時，僅進行不同頻率設計值的點估計，未要求對設計值進行區間估計，忽略了樣本抽樣誤差對設計值的影響，這對海洋水文實測資料短缺海區的海港工程設計極端高低水位和設計波要素可能帶來可觀的誤差。

目前對極值Ⅰ型分布模型的抽樣誤差的研究較少，未見關于該分布設計安全修正值及系數的相關研究，本文通過統計試驗（即Monte Carlo Method）推求基于Gumbel 參數估計方法的極值Ⅰ型分布模型的安全修正值系數B 值、B′值諾模圖，并與Pearson-Ⅲ型概率分布模型的B 值諾模圖進行比較驗證，供工程實踐尤其資料短缺的海港工程設計中參考使用。

1 理論與方法

1.1 設計值抽樣誤差的一般表達

潮位、波浪等海洋水文極值系列是未知的無限總體，而實測資料系列只能是有限樣本，用有限樣本估計總體參數必然存在不可避免的抽樣誤差，即由樣本的特征推算和估計總體的特征而引起的誤差。

用各種估計方法可對隨機變量系列X 估得不同頻率P 時的變數xP值，并繪出頻率曲線。根據抽樣原理，從頻率曲線上讀得的值是點估計值，在其兩側有一定的波動范圍。如圖1 所示，相應A 點對應的設計值xP隱含一個條件分布，其密度函數為f（xP|P）。此時，設在頻率曲線上讀得的估計值為xP，即A 點的縱坐標值，這是該條件分布的均值、中值、眾值或別的特征值。

圖1 頻率分析設計值xP 的條件密度分布圖[4]

1.2 統計試驗推求B 值諾模圖

極值Ⅰ型分布模型表達式如下。

式中，λP，n為極值Ⅰ型分布離均系數，其值可按式（12）計算或查JTS 145—2015《港口與航道水文規范》[1]附錄表D.0.2。可見，JTS 145—2015 與GB 50009—2012 中的極值Ⅰ型分布的參數估計式的表達形式不同，但本質一致。

1.2.2 統計試驗設計

現行技術標準推薦采用Gumbel 法進行極值Ⅰ型分布參數估計，為便于研究成果的實踐應用，下文基于該法估計參數進行統計試驗，研究樣本設計成果的抽樣誤差，進而計算B 值。統計試驗具體步驟如下。

步驟1：樣本序列的生成。為了保證計算結果的一致性和有效性，E（X）= 1，CV=（0.1，0.15，0.20），運用蒙特卡洛方法生成容量n=5、10、15、20、30、50、70 和100 的序列，每組生成N=500個隨機序列作為1 個計算小組。

步驟2：統計參數的確定。按照上述式（8）、式（9）計算每組隨機序列的尺度參數和位置參數。

周淑蘋女士為前任中西女塾皇后。民國十六年間畢業。各報爭刊其玉影。……是晚女士之裝束華麗冠全場。粉紅色的衣綴以一大紅花。衣外復披白狐黑大衣，色彩與式樣可稱時髦已極。女士之頭上浩浩乎平滑無疆，是又不得不歸功于司丹康矣。聞今日為其出閨佳期。記者深愿女士能始終得條件之保障。而享受攜手白頭之樂也。……聞夏璐敏女士為滬大西劇社社員，因家庭反對男女合演，故未加入云。 ”[12]5（圖 8、圖 9）《新聞報》1929年 3月6日刊載《值得登廣告》，詳細記載了該劇的本事：

步驟3：計算xP的均方誤差SxP。采用式（4）確定每組隨機序列的設計值xP，i（P=0.01%，0.1%，1%，2%，3.3%，5%，10%），對于n、E（X）、CV相同的500 組隨機序列，計算其理論設計值x0P、小組均方誤差每500組序列作為一個小組，分別計算10 000 組的SxP，求得10 000 組的平均值。

步驟4：計算B 值。對n = 5、10、15、20、30、50、70 和100 的序列，根據下式計算對應的B值序列。

步驟5：計算相對偏差。對n = 5、10、15、20、30、50、70 和100 的序列，根據下式計算對應的相對偏差δ。

1.3 統計試驗結果分析

1.3.1 CV對B 值的影響

樣本容量取10、20、30、50 和100，對應不同的CV值分別進行試驗，B 值成果如表1 所示。可見同一樣本容量下相同頻率，B 值相當接近，可以認為B 值與CV值無關。

表1 不同CV 值與B 值的關系表（樣本容量n=10、20、30、50 和100 共5 組）

1.3.2 樣本容量對B 值的影響

1.3.3 B 值諾模圖

根據上述分析，B 值與CV值無關，與樣本容量n、頻率P 關系密切。通過上述的15 組試驗方案，統計試驗10 000 次后，通過式（12）計算得到B 值表2 所示，繪制成B 值諾模圖如圖2。為使參數更加清晰，查圖使用更加方便，令從而式（1）變為式（15）。

圖2 極值Ⅰ型分布B 值諾模圖

表2 不同樣本容量與B 值的關系表（CV=0.10、0.15、0.20 共3 組，取均值）

B′值及諾模圖，分別如表3 和圖3 所示。

表3 不同樣本容量與B′值的關系表（CV=0.10、0.15、0.20 共3 組，取均值）

圖3 極值Ⅰ型分布B′值諾模圖

1.3.4 樣本容量與相對偏差的關系

樣本容量取5、 10、 15、 20、 30、 50、 70 和100， 對應不同的CV值（CV=0.1、 0.15、 0.2） 分別進行試驗， 采用式（14） 計算樣本容量引起估計值與真值的相對偏差， 如圖4所示。 在值相同的情況下， 樣本容量越小相對偏差越大； 樣本容量相同的情況下， CV值越大相對偏差越大。

圖4 不同樣本容量引起的估計值與真值的相對偏差

2 極值Ⅰ型分布B 值的合理性分析

Pearson-Ⅲ型概率分布模型也被常用于水文氣象要素極值的頻率分析，許多學者及相關的規范采用了優化適線法、絕對值離差和最小法、最小二乘法、線性矩法等進行參數估計，并通過蒙特卡洛試驗，得到了安全修正值系數B 值[3-10]。Pearson-Ⅲ型分布的B 值是CS與頻率P 的函數，且與參數的估算方法密切相關。為與極值Ⅰ型的B 值對比，本文將n=5、10、15、20、30、50、70 和100 統計得到的B 值取平均值，并與Pearson-Ⅲ型CS= 0、0.5、1.0、1.5 的B 值進行比較分析。

叢樹錚等[7]根據絕對值離差和最小（φ2-W）優化適線法，由統計試驗方法計算得到了B 值圖，如圖5（a）所示；金光炎等[4]用統計試驗方法，采用絕對值離差和最小法及最小二乘法適線法進行分布參數估計，繪制B 值諾模圖，如圖5（b）和圖5（c）所示；劉攀等[8]認為線性矩法估計頻率分布曲線的參數有較好的無偏性和有效性，并采用線性矩法估計Pearson-Ⅲ型的分布參數，繪制了用安全修正值系數B 值諾模圖，如圖5（d）所示。

圖5 基于Gumbel 法估計的極值Ⅰ型與不同參數估計方法的Pearson-Ⅲ型分布的B 值諾模圖

對于Pearson-Ⅲ型分布，在不同的參數估計方法中，當CS相同時，B 值由大至小依次為：最小二乘適線法、絕對值離差和最小法、線性矩法、絕對值離差和最小（φ2-W）優化適線法，其中絕對值離差和最小法與線性矩法參數估計相應的B 值非常接近。極值Ⅰ型分布的B 的均值介于Pearson-Ⅲ型分布CS= 0.5～1.5 的B 值之間。在頻率P 較小時，極值Ⅰ型分布的B 值介于Pearson-Ⅲ型分布CS=0.5～1.0 的B 值之間；在頻率P 較大時，極值Ⅰ型分布的B 值介于Pearson-Ⅲ型分布CS=1.0～1.5的B 值之間；頻率P ＞4%，極值Ⅰ型分布的B 值與Pearson-Ⅲ型分布CS=1.5 的B 值大體相當。

3 實例分析

以某潮位站1989—2018 年共30 年歷年最高潮位資料為例，進行設計高潮位實例分析。該潮位站1989—2018 年歷年最高潮位序列的均值為1.77 m（珠江基面，下同），標準差為0.389，具體系列值如圖6 所示。

圖6 某潮位站1989—2018 年歷年最高潮位

3.1 不同頻率分布模型的對比

采用極值Ⅰ型分布進行頻率分析，參數估計采用我國標準推薦的Gumbel 法估計，得到位置參數μ=1.563，尺度參數α=0.359。采用Pearson-Ⅲ型分布進行頻率分析，采用常用的φ2-W 適線參數估計方法得出CV= 0.26，CS/ CV= 4，對應的位置參數、尺度參數、形狀參數分別為0.885、3.70、4.18。兩個分布考慮安全修正值修正后得到的設計值及B 值匯總如表4 所示。

表4 極值Ⅰ型與Pearson-Ⅲ型分布修正后的設計值

極值I 型分布計算得到的設計潮位值xP較Pearson-Ⅲ型分布得到的成果略大，不同頻率下差異約-0.02～0.25 m；在考慮安全修正值SxP后，兩個分布得到的設計潮位值x*P相當，差異約0.01～0.13 m；頻率越小，兩個分布的差異越大，常用的P=1%、2%下，兩個分布模型得到的xP、x*P基本一致。

在采用極值I 型分布分析計算設計潮位時，安全修正值SxP約為xP的7.6%～14.6%，約為x*P的7.1%～12.7%，即在不考慮安全修正值時，極值I 型分布分析計算設計潮位較x*P小約7.1%～12.7%。在采用Pearson-Ⅲ型分布分析計算設計潮位時，安全修正值SxP約為xP的6.3%～18.0%，約為x*P的5.9%～15.3%，即在不考慮安全修正值時，Pearson-Ⅲ型分布分析計算設計潮位較x*P小約5.9%～15.3%。本例充分說明了為得到更加可靠的設計潮位值，考慮安全修正系數修正抽樣誤差是十分必要的。

3.2 不同樣本容量的對比

采用Gumbel 法對近30 年、20 年、15 年、10年、5 年的歷年最高潮位資料進行設計潮位的計算，考慮安全修正值修正后得到的設計值及B′值匯總如表5、圖7 和圖8 所示。

圖7 不同樣本容量的設計潮位xP 成果對比

圖8 不同樣本容量的設計潮位x*P 成果對比

表5 不同樣本容量對應的設計潮位成果匯總表

隨著樣本容量的增加，設計潮位值xP、 x*P大致呈現減小的趨勢；近10 年、15 年、20 年的設計潮位值xP大小相當，不同樣本容量的xP差異在0.11～0.14 m；以近30 年的xP為基礎，頻率越小，不同樣本長度的xP差距越大，近5 年資料得到的設計值大了約20.3%～52.5%，近10 年、15 年、20 年的設計值大了約2.15%～11.9%。近10 年、15 年、20 年的設計潮位值x*P大小相當，不同樣本容量的xP差異在-0.01～0.23 m；以近30 年的x*P為基礎，頻率越小，不同樣本長度的x*P差距越大，近5 年資料得到的設計值大了約40.8%～87.3%，近10 年、15年、20 年的設計值大了約8.6%～20.0%。可見不同樣本容量對設計值的影響較大，尤其是樣本容量不足10時，引起的偏差達到了50%以上。

在不同樣本容量時，不同頻率的安全修正值SxP約為xP的7.1%～28.9%，約為x*P的7.6%～40.7%，且樣本容量越小，SxP占xP、x*P的比例越大，具體如表6 所示；在使用近5 年的資料時，不同頻率的安全修正值SxP約為xP的20.6%～28.9%，約為x*P的26.0%～40.7%；在使用近10 年的資料時，不同頻率的安全修正值SxP約為xP的12.6%～20.8%，約為x*P的14.5%～26.2%；在使用近15 年的資料時，不同頻率的安全修正值SxP約為xP的10.5%～17.7%，約為x*P的11.7%～21.5%；在使用近20 年的資料時，不同頻率的安全修正值SxP約為xP的9.4%～15.7%，約為x*P的10.3%～18.6%。可見，樣本容量越短，考慮安全修正系數修正抽樣誤差越有必要。對于《港口與航道水文規范》規定的不少于連續20 年的資料要求，本例的誤差達到了10.3%～18.6%，可見考慮抽樣誤差十分必要。

表6 不同樣本容量設計潮位相對誤差匯總表

4 結 論

本文針對基于Gumbel 法進行參數估計的極值Ⅰ型分布，通過統計試驗，研究繪制了安全修正值系數B 值的諾模圖及考慮樣本容量的B′值的諾模圖，與多種方法參數估計的Pearson-Ⅲ型的B 值諾模圖對比分析，并以某潮位站設計最高潮位計算為例加以驗證，得到下列結論。

（1）海港工程設計與結構抗風設計中，常采用極值Ⅰ型分布推求稀遇頻率的極端高水位、極端低水位、設計波要素和設計最大風速。本文基于現行技術標準推薦的Gumbel 參數估計方法進行統計試驗，證明樣本容量過小引起設計值的抽樣誤差不容忽視。

（2）極值Ⅰ型分布的安全修正值系數B 值與離差系數無關，為頻率P、樣本容量n 的函數，可采用本文統計試驗得出的B 值或B′值諾模圖近似估計該分布設計值的抽樣誤差。

（3）根據叢樹錚等、劉攀等[7-8]的研究，對Pearson-Ⅲ型分布，在CS較小時，樣本容量n 對B值的影響較小，可忽略n 對B 值的影響。然而，根據本文統計試驗及實例分析，對于極值Ⅰ型分布，樣本容量n 對B 值進而對極端水位設計值的影響相對顯著，不宜忽略。

（4）本文提出的基于Gumbel 法參數估計的極值Ⅰ型分布的設計安全修正系數B 值介于Pearson-Ⅲ型分布CS=0.5～1.5 的B 值之間。

（5）本文實例驗證了采用Gumbel 法進行極值Ⅰ型分布參數估計，獲得的B 值、B′值諾模圖可為海港工程設計極端高水位、極端低水位、設計波要素，以及結構荷載計算中的最大風速分析計算的設計值抽樣誤差估計提供可靠的依據，并可供相關設計標準修訂時參考。