结构动力特性对结构动力设计和振动分析有重要意义，它反映了结构系统刚度、质量分布及边界条件等信息. 由于构件制造误差、老化损伤及自身固有的随机性等因素影响，结构参数具有不确定性，从而导致结构动力特性具有随机性^[1]. 为了提供更准确的结构动力信息^[2-3]，须定量计算结构参数传递到结构动力的不确定性. 蒙特卡罗模拟法（Monte Carlo simulation, MCS）是结构动力不确定性量化常用的方法，其通过大量的有限元模型计算得到动力响应的统计特性^[4]，不足之处在于计算成本过高. 为了克服MCS方法的不足，代理模型方法构建用于替代结构有限元模型的近似数学模型. 不确定性量化的代理模型包括响应面^[5-6]、多项式混沌展开^[7]、高斯过程模型（Gaussian process model, GPM）^[8]等. 高斯过程模型是一种非参数概率模型，广泛用于不确定性量化^[9-11]，具有模拟灵活性和预测不确定性定量等优点.

代理模型的预测精度受试验设计方法的影响. Santner等^[12-13]研究不同试验设计方法，比较不同的试验设计和代理模型的拟合效果. 对于代理模型建模，空间填充设计是较好的选择，包括拉丁超立方设计^[14-15]、Sobol序列抽样^[16]、Hammersley序列抽样^[17]等，这些能够较好地均匀填充参数变化空间. 采用传统试验设计构建GPM，是事先设定样本的数量，一次生成所有训练样本，会因为过度采样导致计算资源浪费. 序贯设计^[18]是把样本选择和代理模型结合起来，基于填充准则依次选择最佳样本点，自适应地更新当前代理模型，以较少样本建立一个高精度的代理模型.

本研究提出基于序贯设计和高斯过程模型的结构动力不确定性量化方法. 基于GPM的预测方差（mean squared error, MSE）信息，通过最大化预测方差（maximizing mean squared error, MMSE）^[19]依次选择最佳设计点，逐步自适应地更新GPM. 基于建立的自适应GPM，动力特性统计矩的复杂高维积分可以转化为简单一维积分，进而可以解析计算出结构动力特性统计矩.

高斯过程模型是基于贝叶斯理论的一种非参数概率模型，利用高斯过程先验，得到待预测点的后验概率分布. GPM完全由其均值函数m(x)和平方指数协方差函数C(x, x')决定^[9-10]，通常均值函数采用常数形式μ ，平方指数协方差函数计算式为

(1)$ C({\boldsymbol{x}},{\boldsymbol{x}}') = {\eta ^2}\exp \left[ { - \frac{1}{2}\sum\limits_{k = 1}^d {{{\left( {\frac{{{{{x}}_k} - {{x}}_k{'}}}{{{l_k}}}} \right)}^2}} } \right]. $

式中：η²为协方差函数的变化尺度；d为输入参数的维度；x_k为参数x的第k个元素；l_k为协方差函数的变化速率；协方差函数的参数$\theta $为超参数， $\theta = \left\{ {{l_1},{l_2},\cdots,{l_d},{\eta ^2}} \right\}$.

假设有n个观测值的样本点集$ {{D}} = \left\{{{\boldsymbol{X}},{\boldsymbol{Y}}} \right\} $，其中$ {\boldsymbol{X}} = {\left[ {{\boldsymbol{x}}_1^{\rm{T}},{\boldsymbol{x}}_2^{\rm{T}},\cdots,{\boldsymbol{x}}_n^{\rm{T}}} \right]^{\rm{T}}} $、$ {\boldsymbol{Y}} = {\left[ {{{{y}}_1},{{{y}}_2},\cdots,{{{y}}_n}} \right]^{\rm{T}}} $，根据高斯过程先验，模型输出服从高斯分布：

(2)$ P({\boldsymbol{Y}}) = N\left( {{\boldsymbol{\mu}} ,C({\boldsymbol{X}},{\boldsymbol{X}})} \right). $

待预测点${{\boldsymbol{x}}_*}$处的预测值${{{y}}_{{*}}}$与Y也服从高斯分布，计算式为

(3)$ P({\boldsymbol{Y}},{{{y}}_*}) = N\left( {\left[ {\begin{array}{*{20}{c}} {\boldsymbol{\mu}} \\ {{{{\mu}} _*}} \end{array}} \right],\left[ {\begin{array}{*{20}{c}} {C({\boldsymbol{X}},{\boldsymbol{X}})}&{C({\boldsymbol{X}},{{\boldsymbol{x}}_*})} \\ {C({{\boldsymbol{x}}_*},{\boldsymbol{X}})}&{C({{\boldsymbol{x}}_*},{{\boldsymbol{x}}_*})} \end{array}} \right]} \right). $

根据贝叶斯原理，预测值${y_{{*}}}$的后验分布为

(4)$ P({y_*}) = {{P({\boldsymbol{Y}},{{{y}}_*})}}/{{P({\boldsymbol{Y}})}}. $

将式(2)、(3)代入式(4)，得到

(5)$ P({{{y}}_*}) = N({\hat {{y}}_*},{v_{{y_*}}}). $

式中：$ {\hat y_*} $为预测均值，$ {v_{{y_*}}} $为预测方差.

(6)$ {\hat {{y}}_*} = {\mu _*}+{{\boldsymbol{\alpha}} ^{\rm{T}}}{{\boldsymbol{C}}_*}, $

(7)$ {v_{{y_*}}} = {\eta ^2} - {\boldsymbol{C}}_*^{{\rm{T}}}{{\boldsymbol{C}}^{ - 1}}{{\boldsymbol{C}}_*}. $

式中：$ \;{\mu _{{*}}} $为待预测点的均值，${\mu _{{*}}}{\text{ = }}{\left( {{{\boldsymbol{e}}^{\rm{T}}}{{\boldsymbol{C}}^{ - 1}}{\boldsymbol{e}}} \right)^{ - 1}}{{\boldsymbol{e}}^{\rm{T}}}{{\boldsymbol{C}}^{ - 1}}{\boldsymbol{Y}}$，e为长度为n的单位列向量；${\boldsymbol{\alpha}} = {{\boldsymbol{C}}^{ - 1}}({\boldsymbol{Y}} - {\boldsymbol{e}}{\mu _*})$；C_*为待预测点与训练样本点之间的协方差矩阵，表示待预测点与训练样本点的相关性，${{\boldsymbol{C}}_*}{\text{ = }} \left[ {C\left( {{{\boldsymbol{x}}_1},{{\boldsymbol{x}}_*}} \right),}\right. {C}\left.{\left( {{{\boldsymbol{x}}_2},{{\boldsymbol{x}}_*}} \right),\cdots,C\left( {{{\boldsymbol{x}}_n},{{\boldsymbol{x}}_*}} \right)} \right]^{\rm{T}}$；C为训练样本点之间的协方差矩阵，表示训练样本点自身的相关性，${\boldsymbol{C}} = C\left( {{\boldsymbol{X}},{\boldsymbol{X}}} \right)$.

GPM的超参数$\theta = \left\{ {{l_1},{l_2},\cdots,l_d,{\eta ^2}} \right\}$可通过最大化边缘似然函数求得，即最小化负对数边缘似然函数${L}$^[9]：

(8)$ \left. \begin{gathered} {{\hat \varTheta }} = \mathop {\arg \min } \limits_{{{\mathbf{\varTheta }}}} \; {L}({\varTheta }), \\ {L}\left( {{\varTheta }} \right) = \frac{{{n}}}{2}\ln \;\left( {2\text{π} } \right)+\frac{1}{2}\ln \;\left( {\left| {{{\boldsymbol{C}}}} \right|} \right)+ \\ \qquad \frac{1}{2}{\left( {{{\boldsymbol{Y}}} - {\boldsymbol{e}}^{}{\mu }} \right)^{\rm{T}}}{\boldsymbol{C}}^{ - 1}\left( {{{\boldsymbol{Y}}} - {\boldsymbol{e}}^{}{\mu}} \right). \\ \end{gathered} \right\} $

序贯设计是基于样本填充准则，依次选择最有用的样本点填充至初始样本集中，不断地通过新样本点更新代理模型. 序贯设计的关键是如何根据样本点和当前代理模型选择新的样本点，即如何定义样本填充准则. GPM能够提供预测方差信息（mean squared error, MSE），MSE反映GPM与原始物理模型之间的差异. MSE越大，表明GPM的预测误差越大. 在预测方差最大处增加样本点，能够减小GPM的预测误差. 因此提出基于最大化预测方差的样本填充准则，计算式为

(9)$ {{\boldsymbol{x}}_{{\rm{new}}}} = \mathop {\arg \max }\limits_{\boldsymbol{x}}\; v_{{y_*}}^{}({\boldsymbol{x}}). $

式中：$ v_{{y_*}}^{}({\boldsymbol{x}}) $为GPM的预测方差MSE.

样本迭代采用的停止准则^[20]为

(10)$ \max \; v_{{y_*}}^{}({\boldsymbol{x}}) \leqslant r\Delta {\boldsymbol{Y}}. $

式中：r为比例系数，取值为0.01%~1.00%；$ \Delta {\boldsymbol{Y}} $为输出响应中最大值与最小值的差值.

通过MMSE填充准则的动态序贯设计，构建自适应GPM，由于包含最优设计点，可以有效地提高GPM的预测精度. 建立自适应GPM的具体步骤如下.

1) 获取初始训练样本点X；

2) 以X作为输入值计算有限元模型，得到相应的响应输出值Y；

3) 得到初始样本集$ D = \left\{ {{\boldsymbol{X}},{\boldsymbol{Y}}} \right\} $ ，建立初始GPM；

4) $ {{\boldsymbol{x}}_{\rm{new}}} \leftarrow {{\boldsymbol{x}}_{\rm{new}}} = \mathop {\arg \max }\limits_{\boldsymbol{x}} \; v_{{y_*}}^{}({\boldsymbol{x}}) $；

5) 以$ {{\boldsymbol{x}}_{\rm{new}}} $作为输入值计算相应的响应输出值$ {y_{\rm{new}}} $；

6) 更新初始样本集${D'} = \left\{{{\boldsymbol{X}'}},{{\boldsymbol{Y}}'}\right\} \leftarrow {{\boldsymbol{X}}'} = \left\{{\boldsymbol{X}},{{\boldsymbol{x}}}_{\rm{new}}^{\rm{T}}\right\} ， {{\boldsymbol{Y}}'}=\left\{{\boldsymbol{Y}},{y}_{\rm{new}}\right\}$；

7) 利用步骤6) 中的样本集$ {D'} = \left\{ {{{\boldsymbol{X}}'},{{\boldsymbol{Y}}'}} \right\} $更新当前GPM；

8) 判断是否满足停止准则$\max \;{v_{{y_*}}}({\boldsymbol{x}}) \leqslant r\Delta {\boldsymbol{Y}}$：若满足，则停止迭代；否则重复步骤4) ~7) ，直至满足停止准则.

将结构动力特性的输入参数x与输出响应y的关系用自适应GPM来表达，其中输入参数x的概率分布函数为$ p({\boldsymbol{x}}) $. 根据统计原理，动力特性的均值和方差表达式为

(11)$ \begin{split} E(y) =& \int {yp(y)} {\rm{d}}y = \\ &\int {\left[ {yp(y|{\boldsymbol{x}},{{D}}){\rm{d}}y} \right]} p({\boldsymbol{x}}){\rm{d}}{\boldsymbol{x}} = \int {\hat y} p({\boldsymbol{x}}){\rm{d}}{\boldsymbol{x}},\end{split} $

(12)$ \begin{split} V(y) =& \int {{y^2}p(y)} {\rm{d}}y - {E^2}(y)= \\ & \int {{y^2}[p(y|{\boldsymbol{x}},{{D}})p({\boldsymbol{x}}){\rm{d}}{\boldsymbol{x}}]} {\rm{d}}y - {E^2}(y)= \\ & \int {[{y^2}p(y|{\boldsymbol{x}},{{D}}){\rm{d}}y]} p({\boldsymbol{x}}){\rm{d}}{\boldsymbol{x}} - {E^2}(y)= \\ & \int {({v_y}+{{\hat y}^2}} )p({\boldsymbol{x}}){\rm{d}}{\boldsymbol{x}} - {E^2}(y). \end{split} $

利用平方指数协方差函数的分离特性，式（6）、（7）可以转换为

(13)$ {\hat y_*} = {\mu _*}+a\sum\limits_{i = 1}^n {\alpha _i^{}\mathop \Pi \limits_{k = 1}^d N_{{x_k}}^{}\left( {x_k^i,l_k^2} \right)} , $

(14)$ {v_{y_*}} = {\eta ^2} - {a^2}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{A_{i,j}}} } \mathop \Pi \limits_{k = 1}^d N_{{x_k}}^{}\left( {x_k^i,l_k^2} \right)N_{{x_k}}^{}\left( {x_k^j,l_k^2} \right). $

将式（13）、（14）代入式（11）、（12），得到

(15)$ \begin{split} E(y) & =\displaystyle\int\left(\mu+a \sum_{i=1}^n \alpha_i \prod_{k=1}^d N_{x_k}\left(x_k^i, l_k^2\right)\right) p(\boldsymbol{x}) \mathrm{d} \boldsymbol{x}= \\& \mu+a \sum_{i=1}^n \alpha_i \int \prod_{k=1}^d N_{x_k}\left(x_k^i, l_k^2\right) p\left(x_k\right) \mathrm{d} x_k= \\& \mu+a \sum_{i=1}^n \alpha_i \prod_{k=1}^d I_k^i,\end{split} $

(16)$ \begin{split} V(y) =& \int \{ {\eta ^2} - {a^2}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{A_{i,j}}} } \mathop \Pi \limits_{k = 1}^d N_{{x_k}}^{}({{x}}_k^i,2l_k^2)N_{{{{x}}_k}}^{}({{x}}_k^j,2l_k^2)+ {[\mu +a\sum\limits_{i = 1}^n {\alpha _i^{}\mathop \Pi \limits_{k = 1}^d N_{{x_k}}^{}} ({{x}}_k^i,2l_k^2)]^2}\} p({{{\boldsymbol{x}}}}){\rm{d}}{{{\boldsymbol{x}}}} - {E^2}(y)= \\ &{\eta ^2} + {\mu ^2} + 2\mu a \sum\limits_{i = 1}^n {{\alpha _i}\mathop \Pi \limits_{k = 1}^d I_k^i} + {a^2}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{\alpha _i}} } {\alpha _j}\mathop \Pi \limits_{k = 1}^d {N_{x_k^i}}({{x}}_k^j,2l_k^2) \int {{N_{x_k}} \left( {\frac{{{{x}}_k^i + {{x}}_k^j}}{2},\frac{{l_k^2}}{2}} \right)} p({{{x}}_k}){\rm{d}}{{{x}}_k} - {a^2} \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{A_{i,j}}} } \mathop \Pi \limits_{k = 1}^d {N_{x_k^j}}({{x}}_k^i,2l_k^2) \times\\ &\int {{N_{x_k^{}}}\left( {\frac{{{{x}}_k^i+{{x}}_k^j}}{2},\frac{{l_k^2}}{2}} \right)} p({{{x}}_k}){\rm{d}}{{{x}}_k} - {E^2}(y) = \eta _{}^2+{\mu ^2}+2\mu a\sum\limits_{i = 1}^n {{\alpha _i}\mathop \Pi \limits_{k = 1}^d I_k^i} +{a^2}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{\alpha _i}} } {\alpha _j}\mathop \Pi \limits_{k = 1}^d {N_{x_k^i}}({{x}}_k^j,2l_k^2)I_k^{i,j} - \\ & {a^2}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{A_{j,i}}\mathop \Pi \limits_{k = 1}^d } } {N_{x_k^i}}({{x}}_k^j,2l_k^2)I_k^{i,j} - {E^2}(y). \\[-18pt] \end{split} $

式中：$I_k^i$和$I_k^{i,j}$为一维积分；${a_{1(2)}} = \eta _{1(2)}^2{(2\text{π} )^{d/2}}\Pi _{k = 1}^dl_{1(2),k}^{}$；${{{\alpha}} _i}$为${\boldsymbol{\alpha}} $的第i个元素；${A_{i,j}}$为$ {\boldsymbol{C}}_{}^{ - 1} $第i行第j列的元素；${N_{{x_k}}}({{x}}_k^{i},l_k^2) = \dfrac{1}{{{l_k}\sqrt {2\text{π} } }}\exp\; \left[ { - \dfrac{1}{{2l_k^2}}{{\left( {{{{x}}_k} - {{x}}_k^{i}} \right)}^2}} \right]$；${{{x}}_{*,k}}$为待预测点${{\boldsymbol{x}}_*}$的第k个元素；$x_k^{i}$为${\boldsymbol{X}}$的第i 行第k列的元素.

当参数${{{x}}_k}$服从正态分布时，一维积分$I_k^i$和$I_k^{i,j}$计算式为

(17)$ \begin{split} I_k^i =& \int {N_{{x_k}}^ip({x_k}){\rm{d}}{x_k}} = \int {N_{{x_k}}^{}\left( {x_k^i,l_k^2} \right){N_{{x_k}}}(\xi ,{\theta ^2}){\rm{d}}{x_k}}= \\ & {N_{x_k^i}}\left( {\xi ,l_k^2+{\theta ^2}} \right) ,\end{split} $

(18)$ \begin{split} I_k^{i,j} =& \int {{N_{x_k}}\left( {\frac{{x_k^i+x_k^j}}{2},\frac{{l_k^2}}{2}} \right)p({x_k}){\rm{d}}{x_k}} =\\ & \int {N_{{x_k}}\left( {\frac{{x_k^i+x_k^j}}{2},\frac{{l_k^2}}{2}} \right){{N}_{{x_k}}}(\xi ,{\theta ^2}){\rm{d}}{x_k}} = \\ &{N_{{\tfrac{{x_k^i+x_k^j}}{2}}}}\left( {\xi ,\frac{{l_k^2}}{2}+{\theta ^2}} \right). \end{split} $

当参数${x_k}$服从其他分布时，可根据概率相等的原则将其转化为服从正态分布的参数$x_k^{'}$，采用上述解析结果. 不同概率分布转化的表达式为

(19)$ {x'_k} = {\varPhi ^{ - 1}}\left( {F_{{x_k}}^{}({{{x}}_k})} \right). $

式中：$F_{{x_k}}^{}$为参数${{{x}}_k}$的累积分布函数；$\varPhi$为标准正态累积分布函数.

采用一维函数$ {f_1}(x) $来说明本研究所提自适应GPM的拟合过程，表达式为

(20)$ \begin{split} &{f_1}(x) = - (1.4 - 3x)\sin \;(18x);\quad x \in [0,1.2].\end{split} $

初始样本设为4个，根据MMSE填充准则，在每次迭代过程中选择最优的设计点填充至初始样本集中，并更新当前的GPM，停止准则的比例系数r=0.02%. 整个迭代建模的部分动态过程如图1所示，可以看出，随着迭代次数的增加，GPM的预测均值与真实函数曲线的拟合程度越来越高；GPM的预测方差逐渐减小，即图中95%的置信区间面积不断减小，表明所建立的GPM精度越来越高.

二维函数$ {f_2}(x) $用来进一步展示自适应GPM方法的迭代拟合过程，表达式为

(21)$ {{f}_{2}}(x)=\sin\; (0.83\text{π} {{x}_{1}})\cos\; (1.25\text{π} {{x}_{2}});\;\text{ } \text{ }{{x}_{1}},{{x}_{2}}\text{ } \text{ }\in [-1.0,1.0]. $

初始样本设为10个，停止准则的比例系数r=0.02%. 在每次迭代过程中，通过MMSE选择最佳设计点，迭代过程如图2所示. 图中，实线为真实函数，虚线为自适应GPM预测值，圆点为初始设计点，菱形为序贯设计点，正三角形为最优设计点. 对于二维函数，以等高线图呈现拟合程度来提高可视化. 与一维函数类似，随着迭代次数增加，GPM的预测曲线与真实函数曲线之间的拟合优度提高.

为了更清晰地说明基于MMSE的样本填充机制，迭代过程的部分预测方差如图3所示. 可以看出，随着迭代次数的增加，GPM的预测方差明显减小，进一步表明基于MMSE填充准则建立的自适应GPM可以有效地提高模型预测精度.

将基于序贯设计的自适应高斯过程模型方法用于一个双层柱面网壳（见图4）的动力特性不确定性量化. 该网壳跨度为30 m，矢高为6.8 m，长度为50 m，采用两边支承，整个网壳均由截面直径d=80 mm、薄壁厚度为2 mm的圆形薄壁钢管组成. 在ANSYS环境下建立该网壳的有限元模型，采用LINK8杆单元模拟全部杆件，所建有限元模型共326个节点，1 200个杆单元. 该网壳的有限元模型及前4阶振型如图5所示.

假定钢管半径、钢材密度和弹性模量为不确定性参数（见表1），以网壳的前4阶固有频率为分析对象，计算不确定性参数下该结构固有频率的统计特性. 基于MMSE建立自适应GPM用于结构固有频率的统计矩计算.与传统GPM方法对比，大量次数（$2 \times {10^4}$）采样的MCS方法用于近似固有频率统计矩真值. 计算平台：电脑品牌为DELL，处理器为Intel(R) Core(TM) i5-10500CPU@3.10 GHz，有限元软件为ANSYS @19.2. 自适应GPM法、GPM法和MCS法的计算结果如表2~5所示，均采用Sobol序列抽样获取初始样本点，初始样本数为12. 自适应GPM的样本填充停止准则的比例系数r=0.05%. 由表2~5可知，在增加了15个样本点后，自适应GPM已经达到较高的计算精度，均值最大误差为0.003 7%，方差最大误差为0.561 7%，而传统GPM方法在任意增加32个样本点之后才达到相似的精度，表明所提基于序贯设计的自适应高斯过程模型方法能够在保证高精度的同时显著降低计算成本.

(1) 所提自适应GPM方法可以通过较少样本数量达到较高的预测精度，均值和方差的相对误差均较低，而传统GPM则需要较多样本点才能达到相当的计算精度，表明所提方法能够明显降低计算成本.

(2) 迭代过程中统计矩的相对误差的变化表明，所提自适应GPM法的相对误差在几次迭代后明显降低，而传统GPM的相对误差在增加大量样本后仍然较高，表明所提自适应GPM法相对于传统GPM法具有高效率的优势.

(3) 所提自适应高斯过程模型方法只涉及一种样本填充准则，可进一步研究其他序贯抽样方法，如最大熵准则、整体均方误差准则，比较几种抽样方法的计算结果，得到不同抽样准则的使用范围。