风工程中，建筑结构不但要承受过去某一段时间的风速，还要保证在某一规定的时间期限内安全可靠地承受可能经受的风速.自然界中的风速具有随机性，不同时间有不同的规律，因此有必要根据数理统计的方法来求出建筑结构的设计风速，尤其是对一些重要的对风敏感的结构，如输电塔、桥梁、桅杆等^[1-3].大多数荷载规范只能较好地用于建筑结构的风荷载静力分析或是拟静力分析，所以，估算工程场地处重现期内的极值风速是工程抗风设计的首要任务.当以某种极值分布概型拟合风速母样的极值渐近分布时，对重现期内极值风速的估算结果往往与拟合概型和抽样数量有关，桥梁设计规范中规定极值风速分布服从极值-Ⅰ型^[4].

目前进行极值-Ⅰ型分布风速预测的方法主要有最大似然估计法、矩估计法和概率权矩法，这3种估计方法均属于经典统计范畴.在对极值-Ⅰ型分布风速预测时，采用矩估计法获得的极值风速偏保守，概率权矩法偏危险，最大似然估计法虽然较前两种方法的精度高，但是公式复杂^[5].此外，经典统计有3个共同的局限性：一是提高统计推断的精度，主要靠数据多少决定，这对于小样本，往往发生很大困难甚至无能为力；二是在对极值-Ⅰ型风速预测的过程中均假定位置参数和尺度参数是各自独立的参数，而在理论上的极值-Ⅰ型分布模型中，位置参数和尺度参数不是相互独立的；三是仅仅依靠样本信息对参数进行估计，而没有依靠模型的先验信息.因此，前述的最大似然估计法、矩估计法和概率权矩法在对极值-Ⅰ型分布风速进行预测时的精度就受到了限制.为了弥补现有极值-Ⅰ型风速预测方法的不足，本文采用Bayes统计理论^[6]建立了极值-Ⅰ型风速预测方法.该方法有以下特点：1) Bayes统计理论利用样本信息对先验信息进行修正而得到后验信息；2) Bayes统计由于利用了模型的先验信息，因而对于小样本一般也有较好的统计推断效果；3) Bayes统计对于极值-Ⅰ型模型中的位置参数和尺度参数是否相互独立均适用.最后，通过算例验证了该方法的准确性与有效性.

1 基于Bayes理论的极值-Ⅰ型风速预测方法 1.1 极值-Ⅰ型分布

极值-Ⅰ型分布^[7-8]的概率密度函数和累计分布函数为

$ f\left( t \right) = \frac{1}{\sigma }\exp \left[ { - \frac{{t - \mu }}{\sigma } - \exp \left( { - \frac{{t - \mu }}{\sigma }} \right)} \right], $

(1)

$ F\left( t \right) = \exp \left\{ { - \exp \left[ { - \left( {\frac{{t - \mu }}{\sigma }} \right)} \right]} \right\}. $

(2)

式中μ、σ分别为位置参数和尺度参数.

对 (2) 式两边取对数，得到重现期为T(保证率为1-$\frac{1}{T}$) 的风速预测值为

$ {t_{\frac{1}{T}}} = \mu - \sigma \ln \left[ { - \ln \left( {1 - \frac{1}{T}} \right)} \right]. $

(3)

从而，得到百年一遇的风速预测值t_0.01为

$ {t_{0.01}} = \mu - \sigma \ln \left[ { - \ln 0.99} \right] = \mu + 4.6\sigma . $

(4)

1.2 Bayes理论

采用Bayes理论将μ和σ作为随机变量来估计μ和σ的联合概率密度函数π(μ, σ)，下面采用Jeffreys无信息先验分布来对t_{$\frac{1}{T}$}进行估计.

Jeffreys用Fisher信息矩阵行列式的平方根作为 (μ, σ) 先验密度的核，用Jeffreys准则寻找无信息先验分布^{[6, 9]}的步骤如下.

步骤1 写出样本似然函数的对数：

$ \begin{array}{l} L = \ln \left[ {L\left( {t\left| {\mu ,\sigma } \right.} \right)} \right] = \ln \left[ {\prod\limits_{i = 1}^n {f\left( {{t_i}\left| {\mu ,\sigma } \right.} \right)} } \right] = \\ \;\;\;\;\;\;\;\sum\limits_{i = 1}^n {\ln f\left( {{t_i}\left| {\mu ,\sigma } \right.} \right)} . \end{array} $

(5)

步骤2 求Fisher信息矩阵：

$ \mathit{\boldsymbol{I}}\left( {\mu ,\sigma } \right) = E\left( { - \frac{{{\partial ^2}L}}{{\partial \mu \partial \sigma }}} \right). $

(6)

步骤3 求 (μ, σ) 的无信息先验密度函数：

$ \pi \left( {\mu ,\sigma } \right) = {\left[ {\det \mathit{\boldsymbol{I}}\left( {\mu ,\sigma } \right)} \right]^{\frac{1}{2}}}. $

(7)

对于极值-Ⅰ型分布，Fisher信息矩阵为

$ \mathit{\boldsymbol{I}}\left( {\mu ,\sigma } \right) = - E\left| {\begin{array}{*{20}{c}} {\frac{{{\partial ^2}\ln L}}{{\partial {\mu ^2}}}} & {\frac{{{\partial ^2}\ln L}}{{\partial \mu \partial \sigma }}}\\ {\frac{{{\partial ^2}\ln L}}{{\partial \mu \partial \sigma }}} & {\frac{{{\partial ^2}\ln L}}{{\partial {\sigma ^2}}}} \end{array}} \right|. $

(8)

由Jeffreys准则可得

$ \pi \left( {\mu ,\sigma } \right) = \frac{1}{{{\sigma ^2}}}. $

(9)

下面基于π(μ, σ)，采用Lindley近似^[9-11]方法推导t_{$\frac{1}{T}$}的解析表达式.令

$ I = \frac{{\int {u\left( \theta \right)v\left( \theta \right){{\rm{e}}^{L\left( \theta \right)}}{\rm{d}}\theta } }}{{\int {v\left( \theta \right){{\rm{e}}^{L\left( \theta \right)}}{\rm{d}}\theta } }}. $

(10)

其中θ=(θ₁, θ₂, …, θ_k) 为参数向量，L为似然函数的对数.需要注意的是，I为在给定先验分布v(θ) 的情况下u(θ) 的后验期望.

根据Lindley近似方法，由式 (10) 可得

$ \begin{array}{*{20}{c}} {\;{u_1} = \frac{{\partial u}}{{\partial {\theta _1}}},\;\;\;\;{u_2} = \frac{{\partial u}}{{\partial {\theta _2}}};}\\ {{u_{11}} = \frac{{{\partial ^2}u}}{{\partial \theta _1^2}},\;\;\;\;{u_{22}} = \frac{{{\partial ^2}u}}{{\partial \theta _2^2}};}\\ {p = \pi \left( {{\theta _1},{\theta _2}} \right);}\\ {\;{p_1} = \frac{{\partial p}}{{\partial {\theta _1}}},\;\;\;\;{p_2} = \frac{{\partial p}}{{\partial {\theta _2}}};}\\ {{L_{20}} = \frac{{{\partial ^2}L}}{{\partial \theta _1^2}},\;\;\;\;{L_{02}} = \frac{{{\partial ^2}L}}{{\partial \theta _2^2}};}\\ {{L_{30}} = \frac{{{\partial ^3}L}}{{\partial \theta _1^3}},\;\;\;\;{L_{03}} = \frac{{{\partial ^3}L}}{{\partial \theta _2^3}};}\\ {{\sigma _{11}} = {{\left( { - {L_{20}}} \right)}^{ - 1}},\;\;\;\;{\sigma _{22}} = {{\left( { - {L_{02}}} \right)}^{ - 1}}.} \end{array} $

进一步可得

$ \begin{array}{l} E\left( {u\left( \theta \right)\left| {\overrightarrow t } \right.} \right) = u\left( {\overline {{\theta _1}} ,\overline {{\theta _2}} } \right) + \frac{1}{2}\left( {{u_{11}}{\sigma _{11}} + {u_{22}}{\sigma _{22}}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{p_1}{u_1}{\sigma _{11}} + {p_2}{u_2}{\sigma _{22}} + \frac{1}{2}\left( {{L_{30}}{u_1}\sigma _{11}^2 + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{L_{03}}{u_2}\sigma _{22}^2 + {L_{21}}{u_2}{\sigma _{11}}{\sigma _{22}} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{L_{12}}{u_1}{\sigma _{22}}{\sigma _{11}}} \right). \end{array} $

(11)

其中$\overline {{\theta _1}} $、$\overline {{\theta _2}} $分别为θ₁和θ₂的最大似然估计值.

极值-Ⅰ型风速预测表达式为

$ {t_B} = \mu - b\sigma = u\left( {\mu ,\sigma } \right). $

(12)

式中μ=θ₁，σ=θ₂.由此可得u₁=1，u₂=－b，其中b=ln (－ln$\frac{1}{T}$)，u₁₁=0，u₂₂=0.因此

$ \begin{array}{*{20}{c}} {p\left( {{\theta _1},{\theta _2}} \right) = \pi \left( {\mu ,\sigma } \right) = \frac{1}{{{\sigma ^2}}},}\\ {{p_1} = 0,{p_2} = - \frac{2}{{{\sigma ^3}}}.} \end{array} $

令μ、σ分别为μ、σ的最大似然估计值，可得

$ \ln L = - n\ln \sigma - \sum\limits_{i = 1}^n {\frac{{{t_i} - \mu }}{\sigma }} - \sum\limits_{i = 1}^n {{{\rm{e}}^{ - \frac{{{t_i} - \mu }}{\sigma }}}} . $

(13)

所以

$ \begin{array}{l} \frac{{\partial \ln L}}{{\partial \mu }} = \frac{n}{\sigma } - \frac{1}{\sigma }\sum\limits_{i = 1}^n {{{\rm{e}}^{ - \frac{{{t_i} - \mu }}{\sigma }}}} ,\\ {L_{20}} = \frac{{{\partial ^2}\ln L}}{{\partial {\mu ^2}}} = - \frac{1}{{{\sigma ^2}}}\sum\limits_{i = 1}^n {{{\rm{e}}^{ - \frac{{{t_i} - \mu }}{\sigma }}}} ,\\ \frac{{\partial \ln L}}{{\partial \sigma }} = \frac{n}{\sigma } + \sum\limits_{i = 1}^n {\frac{{{t_i} - \mu }}{{{\sigma ^2}}}} - \frac{1}{{{\sigma ^2}}}\sum\limits_{i = 1}^n {{{\rm{e}}^{ - \frac{{{t_i} - \mu }}{\sigma }}}} \left( {{t_i} - \mu } \right),\\ {L_{02}} = \frac{{{\partial ^2}\ln L}}{{\partial {\sigma ^2}}} = \frac{n}{{{\sigma ^2}}} - \left\{ {\sum\limits_{i = 1}^n {2\left( {{t_i} - \mu } \right)\left[ {1 - {{\rm{e}}^{ - \frac{{{t_i} - \mu }}{\sigma }}}} \right]} } \right\}\frac{1}{{{\sigma ^3}}} + \\ \;\;\;\;\;\;\;\;\;\left[ {\sum\limits_{i = 1}^n {{{\rm{e}}^{ - \frac{{{t_i} - \mu }}{\sigma }}}} \left( {{t_i} - \mu } \right)2} \right]\frac{1}{{{\sigma ^4}}},\\ {L_{03}} = \frac{{{\partial ^3}\ln L}}{{\partial {\sigma ^3}}} = \frac{{ - 2n}}{{{\sigma ^3}}} + 6\left\{ {\sum\limits_{i = 1}^n {\left( {{t_i} - \mu } \right)\left[ {1 - {{\rm{e}}^{ - \frac{{{t_i} - \mu }}{\sigma }}}} \right]} } \right\}\frac{1}{{{\sigma ^4}}} + \\ \;\;\;\;\;\;\;\;6\left[ {\sum\limits_{i = 1}^n {{{\rm{e}}^{ - \frac{{{t_i} - \mu }}{\sigma }}}} \left( {{t_i} - \mu } \right)2} \right]\frac{1}{{{\sigma ^5}}} - \left[ {\sum\limits_{i = 1}^n {{{\rm{e}}^{ - \frac{{{t_i} - \mu }}{\sigma }}}} \left( {{t_i} - \mu } \right)3} \right]\frac{1}{{{\sigma ^6}}},\\ {L_{21}} = \frac{\partial }{{\partial \sigma }}\left( {\frac{{{\partial ^2}\ln L}}{{\partial {\mu ^2}}}} \right) = \left[ {\sum\limits_{i = 1}^n {{{\rm{e}}^{ - \frac{{{t_i} - \mu }}{\sigma }}}} } \right]\frac{2}{{{\sigma ^3}}} - \left[ {\sum\limits_{i = 1}^n {{{\rm{e}}^{ - \frac{{{t_i} - \mu }}{\sigma }}}} \left( {{t_i} - \mu } \right)} \right]\frac{1}{{{\sigma ^4}}},\\ {L_{12}} = \frac{\partial }{{\partial \mu }}\left( {\frac{{{\partial ^2}\ln L}}{{\partial {\sigma ^2}}}} \right) = \frac{2}{{{\sigma ^3}}}\left[ {n - \sum\limits_{i = 1}^n {{{\rm{e}}^{ - \frac{{{t_i} - \mu }}{\sigma }}}} } \right] + \\ \;\;\;\;\;\;\;\;\frac{1}{{{\sigma ^4}}}\left[ {\sum\limits_{i = 1}^n {{{\rm{e}}^{ - \frac{{{t_i} - \mu }}{\sigma }}}} \left( {{t_i} - \mu } \right)} \right] - \\ \;\;\;\;\;\;\;\;\frac{1}{{{\sigma ^5}}}\left[ {\sum\limits_{i = 1}^n {{{\left( {{t_i} - \mu } \right)}^2}{{\rm{e}}^{ - \frac{{{t_i} - \mu }}{\sigma }}}} } \right]. \end{array} $

因此，可得t_{$\frac{1}{T}$}的Bayes估计值为

$ \begin{array}{l} \overline {{t_B}} = \overline {{t_{\frac{1}{T}}}} + {p_2}{\mu _2}{\sigma _{22}} + \frac{1}{2}\left( {{L_{30}}\sigma _{11}^2 + {L_{03}}{u_2}\sigma _{22}^2 + } \right.\\ \;\;\;\;\;\;\left. {{L_{21}}{u_2}{\sigma _{11}}{\sigma _{22}} + {L_{12}}{\mu _1}{\sigma _{22}}{\sigma _{11}}} \right). \end{array} $

(14)

1.3 极值-Ⅰ型风速的预测

从上述Bayes估计理论可得极值-Ⅰ型风速预测值为

$ \begin{array}{l} \overline {{t_B}} = \bar \mu + 4.6\bar \sigma + {p_2}{\mu _2}{\sigma _{22}} + \frac{1}{2}\left( {{L_{30}}\sigma _{11}^2 + {L_{03}}{u_2}\sigma _{22}^2 + } \right.\\ \;\;\;\;\;\;\left. {{L_{21}}{u_2}{\sigma _{11}}{\sigma _{22}} + {L_{12}}{\mu _1}{\sigma _{22}}{\sigma _{11}}} \right). \end{array} $

(15)

式中μ、σ分别为μ和σ的最大似然估计值.

2 数值算例

本文调查和收集了安徽安庆宿松县、望江县两个气象站1971—2011年实测的风速资料 (共计744个风速样本)，以此提供伪风速母样概率分布模型中位置参数μ和尺度参数σ的合理取值.基于上述伪风速母样概率分布模型参数的合理取值，建立伪风速母样理论模型，然后将Bayes估计和最大似然估计^[12-16]的重现期为100 a极值-Ⅰ型风速预测值与理论模型值进行比较分析.

在伪风速母样理论模型的建立过程中，首先，假定伪风速母样概率分布模型中的位置参数μ和尺度参数σ为相互独立的随机变量，位置参数μ服从正态分布，尺度参数σ服从均匀分布^[7]；其次，考虑样本数量、位置参数μ的先验方差和尺度参数σ的变化对极值-Ⅰ型风速预测结果的影响；最后，基于位置参数μ和尺度参数σ，采用Monte Carlo法产生伪风速母样.本文采用的极值-Ⅰ型风速预测流程如图 1所示.

极值-Ⅰ型风速预测结果见表 1~3，其中，μ_B表示采用m个先验样本计算位置参数的Bayes估计值，而μ和σ分别表示位置参数和尺度参数的最大似然估计值，t_{$\frac{1}{T}$}表示伪风速母样理论值，t_{$\frac{1}{T}$}表示最大似然估计极值-Ⅰ型风速预测值，t_B表示Bayes估计极值-Ⅰ型风速预测值，$\bar e = \left| {{t_{\frac{1}{T}}} - \overline {{t_{\frac{1}{T}}}} } \right|/{t_{\frac{1}{T}}}$和${e_{\rm{B}}} = \left| {{t_{\frac{1}{T}}} - \overline {{t_{\rm{B}}}} } \right|/{t_{\frac{1}{T}}}$分别表示采用最大似然估计和Bayes估计的极值-Ⅰ型风速预测值的误差.

由表 1~3可以看出：1) 当位置参数先验样本数为50，极值-Ⅰ型风速Bayes估计值比最大似然估计值更接近伪风速母样理论值.2) 随着位置参数先验样本数和伪风速母样样本数的增加，尺度参数的增大，极值-Ⅰ型风速Bayes估计值与最大似然估计值之间的差异越来越小.3) 极值-Ⅰ型风速Bayes估计精度随着伪风速母样样本数的增加而提高.4) 位置参数先验样本数量的多少和先验方差的大小对Bayes估计精度没有影响.5) 在大多数情况下，极值-Ⅰ型风速Bayes估计比最大似然估计精度高.

3 工程场地极值风速预测

选择安徽安庆市宿松县和望江县气象站作为采样测站，调查和收集了两个气象站1971—2011年原始风速记录共2×372个，包含了1971年1月至2011年12月的全部372个月的月最大风速值.选取31个年最大风速值进行百年一遇极值风速预测，采用本文提出的Bayes估计方法预测的安徽安庆宿松县和望江县的百年一遇最大风速值分别为27.89 m/s和24.20 m/s，为安全起见，取27.89 m/s作为本文贝叶斯理论预测的安徽安庆市百年一遇风速值.该计算结果与《公路桥梁抗风设计规范》^[4]附表A规定的安徽安庆市百年一遇风速值27.1 m/s相比误差较小，这表明采用本文提出的贝叶斯方法进行实际工程场地极值风速预测是合理可行的.

4 结论

1) 基于Bayes理论提出了极值-Ⅰ型风速预测方法，采用Monte Carlo法产生伪风速母样，分别进行极值-Ⅰ型风速Bayes估计和最大似然估计，并将两者的估计结果与伪风速母样理论值进行比较.

2) 与最大似然估计相比，采用Bayes估计进行极值-Ⅰ型风速预测精度更高.随着极值-Ⅰ型伪风速母样样本数增加，Bayes估计极值-Ⅰ型风速的误差变小.极值-Ⅰ型分布中位置参数的先验样本数和先验方差均不影响Bayes估计极值-Ⅰ型风速预测精度.

3) 在大样本和大尺度参数下，采用Bayes估计极值-Ⅰ型风速预测值与最大似然估计值的差异较小.