应用经典极值理论对风机极端载荷的预报

周 帅，王迎光，李昕雪

(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）从仿真时长及分块容量两方面检验了外推数据的稳定性，证明了本文所述方法不仅具有良好的外推结果还具有较高的可靠性。