基于規劃求解的暴雨強度公式推求方法研究

李 莉， 冷 藝， 黃 俊， 張 賽， 王子奇， 鄒浩東， 洪毅怡暉， 符榮松

(1.重慶大學 三峽庫區生態環境教育部重點實驗室, 重慶 400045； 2.中國交通建設股份有限公司， 北京 100088； 3.廣東海洋大學 海洋與氣象學院， 廣東 湛江 524088； 4.重慶工業職業技術學院， 重慶 401120)

1 研究背景

近5年來，我國300多個城市均發生了不同程度的內澇災害，而暴雨是形成洪澇災害的主要原因[1]。 研究反映降雨規律的暴雨強度公式和參數對于城鎮防洪、海綿城市建設[2]有著重要的意義。目前我國暴雨強度公式的推求基本按照《室外排水設計規范》(GB50014-2006)推薦的基本步驟和方法進行，另外還有不少文獻也對暴雨強度公式編制中的一些技術問題進行了研究[2]，使得暴雨強度公式的編制逐步規范化。但目前對編制過程中基礎數據異常值的剔除、Cv、Cs值對擬合曲線的影響、有約束的規劃求解法在暴雨強度公式推求的應用等問題的研究較少。

現有的皮爾遜-III分布適線研究[6]和公式推求方法[8-11]很少涉及異常值的剔除，隨著數據的更新、資料年限的增加，降雨資料中出現異常值的概率也隨之增加，這直接影響暴雨強度公式的準確度。Cv、Cs值的大小不僅反映實測數據的分布規律，還影響暴雨強度-降雨歷時-重現期(q-t-P)的關系，只有在明確把握Cv、Cs的意義和影響的基礎上，才能根據樣本分布情況正確估計Cv、Cs的大小關系，得到準確可靠的(q-t-P)表。此外，皮爾遜-III型曲線適線過程中推薦采用的目估適線法具有主觀性，工作量大，精度不夠高，通過約束條件采用規劃求解不僅能解決研究者發現參數出現負值的問題[12]，還能提高擬合速度。因此本文以重慶市酉陽縣1993-2013年共21 a的年最大降雨資料為例，研究了異常值的剔除方法、Cv和Cs值對擬合曲線的影響規律以及目標規劃求解的應用，實現對公式的快速精確的推求。

2 降雨數據資料

本文收集了重慶市酉陽縣1993-2013年共計21 a的逐分鐘降雨資料，按年最大值法取樣后將不同年份5～120 min共9個降雨歷時下的最大雨量統計如下(表1)。

表1 降雨數據資料

3 暴雨強度公式擬合與分析

3.1 異常值的剔除

我國制定的城市暴雨強度公式編制和設計暴雨雨型確定技術導則(以下簡稱《導則》)[13]中提到了數據篩選，但是并未提出數據篩選的方法。異常值可能是由于雪融水、颶風、人為等非正常因素造成的，將它們納入計算范圍會導致擬合曲線與普遍降雨規律偏離較遠，因此需要根據氣象歷史資料剔除高異常值和低異常值。美國關于洪水控制的技術導則[14](以下簡稱Bulletin 17B)中提出，選取樣本時會出現遇到特大值的情況，需要進行專門分析確認。一般特大值是指比相應歷史資料序列的平均值(計入特大值后的均值)大2倍以上的稀遇暴雨值，合理取舍特大值有利于提升理論頻率曲線擬合的準確度。通過分析研究，當數據服從對數皮爾遜分布且偏態值在-3～3之間時，采用如下的方法剔除異常值較為合理：

(1)

(2)

小于qL的所有樣本數據應被視為異常值而被剔除，而大于qH的數據只能在有氣象資料支撐的基礎上剔除最大值。當樣本偏態系數大于0.4時，Bulletin 17B建議先剔除高異常值，偏態系數小于0.4時，先剔除低異常值，當偏態系數介于-0.4～0.4時，對高、低異常值的剔除順序不作要求。

鑒于本文中酉陽的樣本量僅有21 a，異常值出現的概率較低，同時修訂公式時樣本量不應小于20 a，對于每一個降雨歷時在本次擬合中僅需檢驗是否存在一個異常值。同時，酉陽城市防澇僅需達到20年一遇的等級，故采用0.1顯著水平下的KN值。在擬合酉陽暴雨強度公式時，按照上述方法首先分析數據是否存在異常值，然后根據計算分別剔除了原始資料中45、60、90、120 min的最大暴雨強度值。

3.2 Cv、Cs值對擬合曲線的影響規律

皮爾遜-III型分布的概率密度函數為：

(3)

由于累積分布函數表示的是小于等于某一暴雨強度q的概率，而理論頻率是指大于等于某一暴雨強度q的概率,故按下列積分式計算理論頻率：

Ptheoretical=P(q≥qtheoetical)=

(4)

對于公式(4)的求解方法可以參考羅雅文[17]的方法，令t=β(q-b),公式(4)變為：

Ptheoretical=P(t>ttheoretical)

(5)

通過以上分析，對Cv、Cs賦不同的數值計算皮爾遜-III型概率密度曲線及其積分后的累積概率曲線，曲線作圖詳見圖1所示。

由圖1(a)可知，當Cs不變時，Cv越大，概率密度曲線矮且寬，累積分布曲線越來越緩，由此推導Cv越大，皮爾遜-III型曲線可擬合的范圍也越大，此時曲線能擬合出離散程度高的數據。當實際降雨數據計算的概率集中在均值附近時，在進行適線時可適當減小Cv，而當數據離散程度較大時，可以適當增大Cv。由圖1(b)可知，當Cv一定時，Cs越大，概率密度曲線越高且窄，左偏程度越高,累積分布曲線中部越向左偏，上段越陡，下段越平緩,由此可推導，當數據正偏態程度越大，即小于平均值的數據占總數據的比例越大時，應該用更大的Cs/Cv曲線去擬合數據。

綜上所述，在皮爾遜-III型曲線適線過程中，可通過計算樣本均值，分析數據的離散與偏態程度，結合圖1中數據離散程度與Cv、Cs變化的相關關系，來調整Cv/Cs比值，提高適線擬合速度。

3.3 基于規劃求解的分布曲線擬合及適線

為得到指定重現期及指定降雨歷時下的暴雨強度，需按照理論頻率分布擬合曲線，將分散的原始數據擬合得出內在趨勢，為暴雨強度公式的外延計算創造條件。理論頻率分布曲線包括：皮爾遜-III型分布曲線(P-Ⅲ型曲線)、耿貝爾分布曲線(Gumbel曲線)和指數分布曲線。由于指數分布曲線適用于非年最大值法取樣的數據[18]，耿貝爾曲線是皮爾遜-III型曲線Cs=1.140的一個特例，其彈性不足，對樣本適應性較差[19]，因此本研究采用皮爾遜-III型曲線來擬合暴雨強度-降雨歷時-重現期(q-t-P)的關系。

將通過規劃求解優化擬合的皮爾遜-III型分布曲線與經驗頻率點共同體現在海森頻率格紙[20]上，得到如圖2所示的暴雨強度-降雨歷時-重現期(q-t-P)關系圖。

根據以上方法按照皮爾遜-III型分布擬合曲線，擬合的絕對均方誤差為0.0557 mm/min，相對均方誤差為6.27%。同時對未剔除異常值時誤差進行了對比分析，發現未剔除異常值時理論頻率曲線的絕對均方誤差為0.0704 mm/min，相對均方誤差為7.03%；剔除異常值后絕對均方誤差值減小0.0147 mm/min，相對均方誤差值減小0.76%。

3.4 暴雨強度公式推求

通過計算qtheoretical或者由圖2查詢得到一組包含5個重現期(酉陽縣管網按20年一遇防澇等級設計，故只需20、10、5、3、2年的重現期)、9個降雨歷時(5、10、15、20、30、45、60、90、120 min)的數據，分別采用最小二乘法和高斯牛頓法推導暴雨強度公式。

依據《室外排水設計規范》(GB50014-2006)[21]，暴雨強度公式定義為：

(6)

式中：q為暴雨強度，L/(s·hm2)；P為重現期，a；t為降雨歷時，min；A1、b、C、n為與地方暴雨特性有關且需求解的參數，其中A1為雨力參數，即重現期為1 a時的1 min設計雨量，mm ；C為雨力變動參數；b為降雨歷時修正參數，即對暴雨強度公式兩側求對數后能使曲線轉為直線所加的時間參數，min；n為暴雨衰減指數，與重現期有關。

3.4.1 最小二乘法推求結果 為求得暴雨強度公式，對于公式(6)兩側取對數得：

lnq=ln 167A1+ln(1+ClgP)-nln(t+b)

(7)

令y=lnq，b0=ln 167A1，x1=ln(1+ClgP)，b2=-n,x2=ln (t+b)，公式(7)變為：

y=β0+x1+β2x2

(8)

(9)

(10)

應用最小二乘法求解暴雨強度公式曲線與皮爾遜-III頻率分布設計值的比較如圖2所示。最小二乘法求解暴雨強度公式的絕對均方誤差為0.0654 mm/min，相對均方誤差為3.96%。最小二乘法是以殘差平方和最小為目標，對估量值求偏導從而求解參數，這種方法對異常值特別敏感，在計算中收斂速度不是很高，此外對于非線性函數需要進行線性轉化，故最小二乘法在暴雨強度公式的擬合上具有一定局限性。

圖1 不同Cv、Cs對概率密度函數和累積分布函數的影響

圖2 暴雨強度-降雨歷時-重現期(q-t-P)關系的經驗頻率點和理論頻率曲線

(11)

高斯牛頓法采用Matlab程序求解，得暴雨強度公式為：

(12)

應用高斯牛頓法求解暴雨強度公式曲線與皮爾遜-III頻率分布設計值的比較如圖4所示。高斯牛頓法求解暴雨強度公式的絕對均方誤差為0.0502 mm/min，相對均方誤差為5.38%。

圖3最小二乘法求解暴雨強度公式推求值 圖4高斯牛頓法求解暴雨強度公式推求值

與P-III擬合值的比較 與P-III擬合值的比較

4 結 論

針對暴雨強度公式推求中存在的問題，以實際暴雨強度數據結合理論知識對推求過程的幾個問題進行了研究，得到以下幾個結論：

(1)通過研究分析提出了異常值的剔除方法，結合分析Cv、Cs的意義總結其對適線的影響規律。當數據服從對數皮爾遜分布且偏態值在-3～3之間時，采用對數形式的高低限值方法剔除異常值，合理取舍異常值有利于提升理論頻率曲線擬合的準確度；當數據集中分布在均值左右時適當減小Cv，當數據正偏程度較大時，適當增大Cs/Cv，為方便準確地擬合出暴雨強度-降雨歷時-重現期(q-t-P)曲線提供了支撐。

(2)通過約束條件，設置目標函數采用Excel規劃求解Cv、Cs，與目估適線法相比，明顯提高了擬合速度和精度。

(3)基于規劃求解的最小二乘法推導的暴雨強度公式絕對均方誤差為0.0654 mm/min，相對均方誤差為3.96%；高斯牛頓法推導的暴雨強度公式絕對均方誤差為0.0502 mm/min，相對均方誤差為5.38%，通過高斯牛頓法推求的酉陽暴雨強度公式較合理。高斯牛頓法可采用最小二乘法的計算值作為初始值進行迭代計算，相比將而言高斯牛頓法具有更快的收斂速度和更高的擬合精度。