應用經典極值理論對風機極端載荷的預報

周 帥，王迎光，李昕雪

(1. 上海交通大學 船舶海洋與建筑工程學院，上海 200240；2. 海洋工程國家重點實驗室，上海 200240；3. 高新船舶與深海開發裝備協同創新中心，上海 200240)

0 引 言

為了保證風機結構的完整性，國際電工委員會發布的陸上風機標準IEC 61400-1[1]和海上風機標準IEC 61400-3[2]，對風機的多個部分均提出了極限強度分析的要求，需要計算構件在一年或多年重現周期下可能承受的最大載荷。

夏一青等[3]利用區組模型和超越門檻值法聯合Gumbel分布求解了某5 MW單樁型和浮式型風機葉片根部的極端面外彎矩；張友文等[4]以分塊法聯合Gumbel分布求解了某5 MW陸上風機和海上風機的極端塔筒基底側向彎矩等載荷；李昕雪等[5]比較了廣義極值分布和廣義帕累托分布在風機極端載荷預報中的差別。但是上述研究工作主要集中在短期載荷極值的求解上，對載荷的長期超越概率分布考慮較少。

P.J. Moriarty等[6]使用超越門檻值法選取WP_Baseline 1.5 MW陸上風機葉片根部的短期載荷極值，然后運用多種類型的分布進行擬合，求解了具有1年和50年重現周期的極端載荷。為了對比驗證，進行了時長1年的直接仿真，獲取了相應載荷的“真實數據”。但是對比之后，發現外推結果與“真實數據”相差較大。P.J. Moriarty等分析其主要原因是短期載荷分布對數據尾部的擬合效果不佳，建議采用擬合效果更好的分布來減少估計高分位數時的誤差。

本文承接P.J. Moriarty的建議，以經典極值理論為基礎，采用分塊法選取樣本點，以廣義極值分布進行擬合，由線性矩法估計分布參數，得到了與“真實數據”貼近的外推值，并且對外推的穩定性進行檢驗，從而證實了P.J. Moriarty想法的可行性。此外，Korn Saranyasoontorn等[7]另辟蹊徑，利用環境等值線法直接求取長期載荷，避開了擬合分布的問題，同樣得到了較好的效果。

1 基本原理

1.1 經典極值理論

設 X1， X2，…， Xn是獨立同分布的隨機變量序列，如果存在常數序列 {an＞0}和{ bn}，使得

成立，其中 H(x)是非退化的分布函數，那么 H 必屬于下列3種類型之一。

Ⅰ型分布：

Ⅱ型分布：

Ⅲ型分布：

其中： an和 bn為規范化常數，為極值指標。Ⅰ型分布稱為Gumbel分布，Ⅱ型分布稱為Fréchet分布，Ⅲ型分布稱為Weibull分布，這3種分布統稱為極值分布。由累積分布函數容易求得各分布的概率密度函數，分別記為 h1(x)， h2(x)， h3(x)。當=3.6時，這3種分布的密度函數如圖1所示。注意極值Ⅱ型分布和極值Ⅲ型分布的密度函數分別具有有限的下端點和上端點。

圖 1 3 種極值分布的概率密度函數Fig. 1 Probability density functions of 3 extreme value distributions

則可以用統一的形式表示上述3種不同類型的極值分布

其中 μ ,ξ∈ R ， σ ＞0。 H 稱為廣義極值分布，簡記為GEV分布， ξ稱為形狀參數。當 ξ ＞0時，取 α =1/ξ，當 ξ ＜0時，取 α =?1/ξ ，當ξ趨近于0時，有因此極值分布的類型完全由形狀參數確定，與位置參數和尺度參數無關。

設獨立同分布序列 X 滿足的分布函數為 F (x)，稱F(x)的廣義反函數

為它的分位數函數。 xp=F?1(p)稱為F的p分位數。在實際應用中，F(x)在其支撐上一般總是單調連續的，因(此)其廣義反函數 F?1(p)是普通的反函數，因此有Fxp=p。

假設 X1， X2,···,XT為年觀測最大值序列，現對某個閾值 u 研究超閾值事件 {Xi＞u}。定義 u 為 T 年重現水平，則在 T 年的時間內，超閾值事件 {Xi＞u}發生的平均次數為1，即

顯然，由

即可得

為式（1）–式（5）的解。由此可知一方面 u(T)是F的(1?1/T)分位數，另一方面年觀測極值超過 T 年重現水平 u(T)的概率為1 / T 。通常用某一小概率值 p代替分位數 u (T)，稱 T=1/(1?p)為 T 年重現周期。

在風機極值載荷的求解中，仿真時長通常設為10 min。現考慮重現周期 T 為1年，則在1年的時長內，總共有 N =1×365×24×6=52 560個 10 min，由

得

稱1/N 為重現周期 T 對應的超越概率，則1年重現周期對應的超越概率為1.902 6E–5，同理有20年重現周期對應的超越概率為9.512 9E–7。

1.2 直接積分法極限載荷外推原理

在直接積分法中，要求的重現周期 T 確定了相應的超越概率 PT，然后以 PT為目標求取相應的極限載荷 lT。

其中： X為多維的環境變量，包含風、浪、流等因素；fX(x)為環境變量的聯合概率密度函數。對于陸上風機，一般僅以風參數作為影響葉片載荷的主要因素。若單次隨機仿真的有效時長為10 min，則上式可記為

其中： L10min為從10 min仿真數據中選取的極大值；Vin， Vout分別為切入和切出風速。在實際求解中，將區間 Vin至 Vout劃分成多個子區間，對每個子區間進行多次隨機仿真，以獲取足夠多的極值點來擬合該風速下短期載荷極值的分布 FL10min|V(l|v)，以此求解式（1）～式（9）中短期極值載荷的超越概率Pr{L10min＞lT|V=v}。將所有子區間的超越概率累加起來可得總體的超越概率，選取不同 lT的值使等式（1）～式（9）左右兩邊相等，此時的 lT即為要求的具有年重現周期的極端載荷。

1.3 線性矩法參數估計

在求解 FL10min|V(l|v)的過程中，需要估計分布參數。通常的方法有矩法、極大似然法、最小二乘法，線性回歸法等[8–10]。羅純[11]在外推求解某水文站最高水位時，曾用概率加權矩法解決了Gumbel分布參數估計的問題，效果良好，因此本文采用與概率加權矩法類似的線性矩法進行參數估計。

線性矩法由Hosking[12]于1990年定義，目前在水文、氣象等領域已廣泛應用。陳元芳等[13–14]利用線性矩法分別研究了在P-Ⅲ型分布和GEV分布下可考慮歷史洪水信息的參數估計問題。梁玉音等[15]、吳俊梅等[16]利用線性矩法研究了太湖流域暴雨頻率的問題。上述研究工作均通過線性矩法取得了不錯的成效。

設 X 為實數域上的隨機變量，其分布函數為F(x)。X1:n≤ X2:n≤ ···≤ Xn:n為遞增的次序統計量，定義總體的線性矩如下：

樣本的線性矩可表示為：

則前4階樣本線性矩計算公式為：

為了更好地描述統計曲線的分布特征，類似常規矩法中的變異系數、偏態系數、峰度系數，定義線性矩系數如下：

廣義極值分布的線性矩及其系數與分布參數之間的關系如下：

尺度參數 σ 和位置參數 μ的估計公式如下：

2 仿真基本參數

本文的仿真工具以美國國家可再生能源實驗室開發的湍流風仿真軟件TurbSim和風機仿真軟件FAST為基礎。風機模型為WP_Baseline 1.5MW陸上風機，基本參數見表1。

風機等級IECⅠ-A。計算工況為IEC 61400-1中正常發電工況類的DLC 1.1。假定平均風速服從瑞利分布。

其中 μV為平均風速，對于IECⅠ-A型風機 μV取10 m/s。只考慮風機在正常發電時失效的情況，因此對5 m/s以下和25 m/s以上的風速進行截斷，將風速的累積分布函數 FV(v)改寫為

表 1 風機基本參數Tab. 1 Basic parameters of wind turbine

將工作風速區間按照2 m/s的分辨率劃分成11個子區間，各子區間的平均風速為 5 m/s，7 m/s，…，25 m/s。對每一個子區間進行8次隨機仿真，單次仿真時長630 s，其中前30 s是為了消除風機啟動時的影響，在選取極值點時會被忽略。湍流譜為IEC Kaimal譜，冪指數風切。

3 極限載荷求解

3.1 短期載荷極值

采用分塊法選取短期載荷極值，即將短期載荷數據分為若干塊，然后從每一塊中選取各自的最大值。設定輸出步長0.05 s，則單次仿真輸出12 000個有效數據。以風速11 m/s，塊容量400為例，選取的短期極值點如圖2所示。

圖 2 短期載荷極值點Fig. 2 Short-term load extremes

設分塊數為m，則在累積分布函數上，分塊極值點與全局極值點的關系如下：

選取了短期載荷極值點后，即可根據前述原理對樣本點進行擬合。風速11 m/s，塊容量400時，擬合分布的概率密度函數和累積分布函數如圖3和圖4所示。

圖 3 擬合分布概率密度函數Fig. 3 Probability density function of fitted distribution

圖 4 擬合分布累積分布函數Fig. 4 Cumulative distribution function of fitted distribution

在統計分析中，Q-Q圖是檢驗模型分布最常用也是最直觀的方法。同樣以風速11 m/s，塊容量400為例，經驗分布函數和擬合分布函數所形成的Q-Q圖如圖5所示，其中經驗頻率公式為中值公式。

圖 5 短期載荷極值 Q-Q 圖Fig. 5 Quantile-Quantile plot of short-term load extremes

由圖可知，擬合分布與經驗分布總體較接近，曲線右側高分位數之間的誤差較小，這對提高外推精度無疑是有利的。因為總體的超越概率是由多個子區間的超越概率累加起來的，因此每個子區間誤差的細微減小，也能在總體上較大地減少誤差。

3.2 長期載荷分布

利用1.3節原理對選取的短期載荷極值進行統計分析，得到不同風速下的形狀參數，位置參數和尺度參數 σ 。注意樣本服從廣義極值分布需要滿足1+ξ(x?μ)/σ＞ 0的要求，當有樣本點不滿足該條件時，本文采取的做法是在Q-Q圖中尋找與擬合曲線差值最大的樣本點，將其刪除，然后重新計算分布參數，如果仍不滿足條件則繼續迭代上述步驟，直至所有樣本點均滿足要求。從實際計算來看，僅在極少數情況下需要刪除不合理的樣本點。以塊容量400為例，參數估計結果如表2所示。

表 2 擬合分布參數估計Tab. 2 Estimated parameters of fitted distribution

由表可知，各平均風速下的形狀參數均為負數，因此所有的短期載荷分布均為廣義極值Ⅲ型分布。

將載荷短期分布和平均風速分布聯合起來，即可求得葉片根部面外彎矩的長期分布。選取不同的值，得到相應的超越概率曲線，與“真實數據”和其他同類研究結果的對比如圖6和表2所示。

圖6中實線部分為“真實數據”，短虛線為本文外推結果，長虛線為P.J. Moriarty外推結果，方形點為Korn Saranyasoontorn外推結果。由圖6可知，本文計算出的超越概率曲線與“真實數據”整體較接近，兩者在500～2 500 kN·m的區間內，都具有相似的形態。P.J. Moriarty 計算出的超越概率曲線在 1 750 kN·m之前與“真實數據”差別不大，但是在該點之后曲線斜率出現了較大的偏差，以致在1年和20年重現周期時較遠的偏離了“真實數據”點。Korn Saranyasoontorn的計算結果沒有超越概率曲線，是通過環境等值線法外直接外推1年和20年重現周期時的極端載荷，由圖可知，其結果與“真實數據”較吻合。同時由表3可知，在20年重現周期時，本文的外推結果與環境等值線法也大致相當。

圖 6 外推結果與“真實數據”對比Fig. 6 Comparison of extrapolation results and real data

表 3 1年和20年重現周期外推結果對比Tab. 3 Comparison of extrapolation in 1 year and 20 years return periods

3.3 數據穩定性

為了檢驗外推結果的穩定性，將仿真時長分別增加至20 min和40 min，同時保持其他設置參數不變。多個分塊容量下不同仿真時長在1年重現周期時的結果對比如表4所示。

表 4 多個分塊容量及仿真時長下的外推結果Tab. 4 Extrapolations of multiple simulation times and block capacities

當仿真時長增加后，外推結果仍然具有較好的穩定性。在多個塊容量下，外推值都集中在2 300 kN·m左右，波動幅度5%以下。同時，對于每個仿真時長，塊容量的變化對外推值的影響也較小。

4 結 語

以經典極值理論和積分法為基礎，采用分塊法選取短期載荷極值，廣義極值分布進行樣本擬合，線性矩法估計分布參數，外推求解了WP_Baseline 1.5 MW陸上風機葉片根部面外彎矩的長期分布，得出以下結論：

1）采用廣義極值分布對短期載荷極值進行擬合，解決了P.J. Moriarty使用其他分布時遇到的擬合問題，減少了估計高分位數時的誤差，最終在極端載荷長期分布的外推上獲得了更好的效果。

2）利用線性矩法估計分布參數，得到了較好的參數估計值，并且經過了分位數對比圖的檢驗。

3）從仿真時長及分塊容量兩方面檢驗了外推數據的穩定性，證明了本文所述方法不僅具有良好的外推結果還具有較高的可靠性。