應(yīng)用經(jīng)典極值理論對風(fēng)機(jī)極端載荷的預(yù)報(bào)

周 帥，王迎光，李昕雪

(1. 上海交通大學(xué) 船舶海洋與建筑工程學(xué)院，上海 200240；2. 海洋工程國家重點(diǎn)實(shí)驗(yàn)室，上海 200240；3. 高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240)

0 引 言

為了保證風(fēng)機(jī)結(jié)構(gòu)的完整性，國際電工委員會(huì)發(fā)布的陸上風(fēng)機(jī)標(biāo)準(zhǔn)IEC 61400-1[1]和海上風(fēng)機(jī)標(biāo)準(zhǔn)IEC 61400-3[2]，對風(fēng)機(jī)的多個(gè)部分均提出了極限強(qiáng)度分析的要求，需要計(jì)算構(gòu)件在一年或多年重現(xiàn)周期下可能承受的最大載荷。

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

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

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

1 基本原理

1.1 經(jīng)典極值理論

設(shè) X1， X2，…， Xn是獨(dú)立同分布的隨機(jī)變量序列，如果存在常數(shù)序列 {an＞0}和{ bn}，使得

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

Ⅰ型分布：

Ⅱ型分布：

Ⅲ型分布：

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

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

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

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

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

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

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

顯然，由

即可得

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

在風(fēng)機(jī)極值載荷的求解中，仿真時(shí)長通常設(shè)為10 min?，F(xiàn)考慮重現(xiàn)周期 T 為1年，則在1年的時(shí)長內(nèi)，總共有 N =1×365×24×6=52 560個(gè) 10 min，由

得

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

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

在直接積分法中，要求的重現(xiàn)周期 T 確定了相應(yīng)的超越概率 PT，然后以 PT為目標(biāo)求取相應(yīng)的極限載荷 lT。

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

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

1.3 線性矩法參數(shù)估計(jì)

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

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

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

樣本的線性矩可表示為：

則前4階樣本線性矩計(jì)算公式為：

為了更好地描述統(tǒng)計(jì)曲線的分布特征，類似常規(guī)矩法中的變異系數(shù)、偏態(tài)系數(shù)、峰度系數(shù)，定義線性矩系數(shù)如下：

廣義極值分布的線性矩及其系數(shù)與分布參數(shù)之間的關(guān)系如下：

尺度參數(shù) σ 和位置參數(shù) μ的估計(jì)公式如下：

2 仿真基本參數(shù)

本文的仿真工具以美國國家可再生能源實(shí)驗(yàn)室開發(fā)的湍流風(fēng)仿真軟件TurbSim和風(fēng)機(jī)仿真軟件FAST為基礎(chǔ)。風(fēng)機(jī)模型為WP_Baseline 1.5MW陸上風(fēng)機(jī)，基本參數(shù)見表1。

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

其中 μV為平均風(fēng)速，對于IECⅠ-A型風(fēng)機(jī) μV取10 m/s。只考慮風(fēng)機(jī)在正常發(fā)電時(shí)失效的情況，因此對5 m/s以下和25 m/s以上的風(fēng)速進(jìn)行截?cái)?，將風(fēng)速的累積分布函數(shù) FV(v)改寫為

表 1 風(fēng)機(jī)基本參數(shù)Tab. 1 Basic parameters of wind turbine

將工作風(fēng)速區(qū)間按照2 m/s的分辨率劃分成11個(gè)子區(qū)間，各子區(qū)間的平均風(fēng)速為 5 m/s，7 m/s，…，25 m/s。對每一個(gè)子區(qū)間進(jìn)行8次隨機(jī)仿真，單次仿真時(shí)長630 s，其中前30 s是為了消除風(fēng)機(jī)啟動(dòng)時(shí)的影響，在選取極值點(diǎn)時(shí)會(huì)被忽略。湍流譜為IEC Kaimal譜，冪指數(shù)風(fēng)切。

3 極限載荷求解

3.1 短期載荷極值

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

圖 2 短期載荷極值點(diǎn)Fig. 2 Short-term load extremes

設(shè)分塊數(shù)為m，則在累積分布函數(shù)上，分塊極值點(diǎn)與全局極值點(diǎn)的關(guān)系如下：

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

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

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

在統(tǒng)計(jì)分析中，Q-Q圖是檢驗(yàn)?zāi)Ｐ头植甲畛Ｓ靡彩亲钪庇^的方法。同樣以風(fēng)速11 m/s，塊容量400為例，經(jīng)驗(yàn)分布函數(shù)和擬合分布函數(shù)所形成的Q-Q圖如圖5所示，其中經(jīng)驗(yàn)頻率公式為中值公式。

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

由圖可知，擬合分布與經(jīng)驗(yàn)分布總體較接近，曲線右側(cè)高分位數(shù)之間的誤差較小，這對提高外推精度無疑是有利的。因?yàn)榭傮w的超越概率是由多個(gè)子區(qū)間的超越概率累加起來的，因此每個(gè)子區(qū)間誤差的細(xì)微減小，也能在總體上較大地減少誤差。

3.2 長期載荷分布

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

表 2 擬合分布參數(shù)估計(jì)Tab. 2 Estimated parameters of fitted distribution

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

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

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

圖 6 外推結(jié)果與“真實(shí)數(shù)據(jù)”對比Fig. 6 Comparison of extrapolation results and real data

表 3 1年和20年重現(xiàn)周期外推結(jié)果對比Tab. 3 Comparison of extrapolation in 1 year and 20 years return periods

3.3 數(shù)據(jù)穩(wěn)定性

為了檢驗(yàn)外推結(jié)果的穩(wěn)定性，將仿真時(shí)長分別增加至20 min和40 min，同時(shí)保持其他設(shè)置參數(shù)不變。多個(gè)分塊容量下不同仿真時(shí)長在1年重現(xiàn)周期時(shí)的結(jié)果對比如表4所示。

表 4 多個(gè)分塊容量及仿真時(shí)長下的外推結(jié)果Tab. 4 Extrapolations of multiple simulation times and block capacities

當(dāng)仿真時(shí)長增加后，外推結(jié)果仍然具有較好的穩(wěn)定性。在多個(gè)塊容量下，外推值都集中在2 300 kN·m左右，波動(dòng)幅度5%以下。同時(shí)，對于每個(gè)仿真時(shí)長，塊容量的變化對外推值的影響也較小。

4 結(jié) 語

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

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

2）利用線性矩法估計(jì)分布參數(shù)，得到了較好的參數(shù)估計(jì)值，并且經(jīng)過了分位數(shù)對比圖的檢驗(yàn)。

3）從仿真時(shí)長及分塊容量兩方面檢驗(yàn)了外推數(shù)據(jù)的穩(wěn)定性，證明了本文所述方法不僅具有良好的外推結(jié)果還具有較高的可靠性。