0 引言

由于实际材料的变异性及建造过程中的随机性，有限元模型与实际桥梁结构之间存在较大差异，具体表现在结构静动力响应上均与实际结构存在差别。合理的有限元模型是进行结构分析和结构评估的基础，因此，有限元模型修正是桥梁评估的基础工作^[1-2]。

有限元模型修正经过了从矩阵到参数、从低维到高维、从动力响应到静力或联合静动力响应、从局部求解到全局优化、从确定性到不确定性等发展历程^[3-7]，各种新的模型修正方法层出不穷。在所有的模型修正方法中，基于响应面法的代理模型修正方法由于其精度好、效率高的优点被研究人员广泛采用^[8-16]，应用最多的是基于多项式建立的响应面模型^{[9, 14]}。但目前在多项式阶次选取、模型精度等方面的研究较少，因此在实际运用中难以准确把握多项式的具体建立规则。

基于响应面的有限元模型修正以结构动力或静动力试验结果为优化目标^{[1, 4, 13, 17]}，通过求解建立的响应面模型，得到有限元模型参数的变化量，从而达到修正有限元模型的目的。多数模型修正过程以桥梁模态参数作为修正目标，桥梁模态参数一般通过分析桥梁振动加速度数据获得，由于环境干扰等因素，实际识别得到的桥梁结构各阶模态参数往往与有限元结果存在不均衡差异，有的阶次差异较大，有的阶次差异较小^[18-20]。现有模型修正方法往往关注修正过程本身，忽略优化目标与有限元结果之间的差异。采用通用方法处理修正过程中的迭代优化问题很难控制局部误差，很容易将差异较大阶次的误差传递到其他差异较小的阶次，对其他阶次修正结果产生影响，从而影响整体修正效果。

本研究针对现有基于响应面的有限元模型修正过程存在的多项式阶次选取和指标误差不能局部控制的求解问题，依托某座斜拉桥，采用基于响应面法的有限元模型修正过程，分析多项式阶次与响应面精度、计算量之间的关系。将稳健估计法引入响应面模型的优化求解过程来解决优化指标波动较大时的误差传递问题，在优化目标波动较大时，仍能得到较好的修正结果，保证模型修正的修正精度。

1 基于响应面的模型修正 1.1 响应面法

响应面法(Response Surface Methodology，RSM) 是利用合理的试验设计方法并通过试验得到一定数据，采用多元二次回归方程来拟合因素与响应值之间的函数关系，通过对回归方程的分析来寻求最优工艺参数，解决多变量问题的一种统计方法。响应面方法被广泛用于影响因素与响应没有显示关系的优化问题中。

1.2 基于响应面的模型修正

为了避免模型修正过程中反复调用有限元模型，提高计算效率，并将修正参数和修正结果的关系可视化，响应面模型被用于桥梁结构模型修正。通过建立结构参数与结构频率响应之间的响应面模型。基于响应面的模型修正包含参数灵敏度分析和显著性检验、响应面拟合、响应面精度检验和优化求解4步。

(1) 参数敏感性分析和显著性检验。参数敏感性分析的目的是选取待修正参数，通过参数敏感性分析挑选出对修正结果影响较大的参数作为待修正参数，并将对修正结果影响较小的参数剔除。在筛选出待修正参数并完成试验设计后，采用显著性检验法确定参与多项式响应面拟合的参数及参数组合。通过参数敏感性分析和显著性检验能够保证在不影响修正结果的情况下减小响应面模型的待求解参数，大大减小响应面模型求解的计算量。

敏感性分析一般采用灵敏度分析法，采用如下指标。

(1)

式中, S_ij^f为第j个响应特征对参数x_i的灵敏度；x_i⁰, x_i¹分别为2次计算的参数取值；f_j⁰, f_j¹分别为2次计算的第j个响应值；n为参数总个数，m为响应特征总个数，由实际采集的数据确定。

显著性检验一般采用F检验法进行，首先计算各个参数的显著性水平值，随后通过检验假设，筛选出对结构响应影响显著的设计参数。显著性水平指标F_m的计算式为:

(2)

式中, SSE(x₁, x₂, x₃, …, x_m-1)为m-1个参数x₁, x₂, x₃, …, x_m-1的回归模型的误差平方和; SSE(x₁, x₂, x₃, …, x_m-1, x_m)为含有m个参数x₁, x₂, x₃, …, x_m-1, x_m的回归模型的误差平方和; n为参数总个数。

一般给定临界显著性水平0.05，如果样本参数概率P≤0.05, 则为特征量影响显著，应该参与响应面拟合；否则，认为不显著，不将该参数组合用于响应面拟合。

(2) 响应面拟合。响应面拟合的目的是建立待修正参数和结构响应之间的显式关系。响应面模型可选的基函数有多项式函数、径向基函数和非线性函数等，一般常用多项式函数建立响应面模型。常用的二阶多项式形式为：

(3)

式中β₀, β_i, β_ii, β_ij为各个待定系数。

在响应面拟合过程中，合理选择响应面模型的阶次是基于响应面的模型修正的关键，多项式模型阶次较低可能引起响应面模型精度不足。阶次较高时，响应面模型需要求解的系数参量大幅增加，计算量显著增加。因此，选择合适的多项式阶次是平衡求解精度和求解效率的关键。

(3) 响应面精度检验。该检验的目的是确定响应面模型的拟合精度是否满足要求，是否能代表实际结构的真实响应。具体过程为在整个设计空间内随机选取有限样本点，并尽量避免选择的样本点是初始试验设计的样本点，将选出的样本点分别代入有限元模型和拟合的模型，得到各自模型的响应特征值，通过一定的统计评价指标，判断响应面的拟合效果。响应面精度的评估方法有复相关系数R²检验和均方根误差RMSE检验等，其中R²检验及RMSE检验是比较常用的检验方法，具体方法为：

(4)

(5)

式中，N为检验的样本点个数；y_RS(j)为第j个样本点的响应面模型响应值；y(j)为第j个样本点的有限元模型响应值；y为样本点所有值的平均值。检验的计算结果在0与1之间，复相关系数R²值若接近1，相对均方根误差RMSE值若接近0，则证明响应面模型的效果较好，可以较好地反映原有的模型特性；若计算检验的结果不好，则需从头开始重新进行响应面拟合。

(4) 响应面优化。响应面优化指在参数设定范围内，通过优化算法对响应面方程组进行有约束的优化求解过程，可归结为以下有约束的优化问题：

(6)

式中，F(x)为实测响应与响应面模型的差值函数；f_exp为实测响应；f_RSM为响应面模型得到的结构响应；lb, ub分别为x的下限和上限。求解上述问题主要采用最小二乘法。

2 稳健估计法

基于响应面的模型修正最终转化为多元方程优化求解问题，传统求解主要采用最小二乘法。基于动力的模型修正以桥梁结构实际模态参数为优化目标，由于桥梁实际运营过程中的监测数据存在多种不可量化的干扰因素，基于监测数据识别得到的桥梁模态参数很可能具有较强的离散性，各阶修正目标与有限元结果之间的偏差分布不均匀，现有修正方法在修正过程中将各阶目标函数统一对待，在这样的修正过程下，不可避免地会将偏差较大的阶次的影响传递到其他阶次，从而影响整体修正结果。为控制修正目标的误差影响程度，尽可能使模型修正结果不受偏差较大的模态参数指标影响，本研究提出基于稳健估计的响应面模型优化方法。

稳健估计的提出是为了解决观测值中存在的异常值引起优化结果的整体偏离问题。通过稳健估计，当观测值与目标值的偏差较小时，能得到精确的估计值；当观测值与理想值的偏差较大时，参数估计值的性能也不会太差。稳健估计中常用的是M估计，通过定义一个ρ函数对不同观测值进行赋权，进而求取参数的最佳估计，其过程如下。

设观测值为L，待求解的未知参数向量为，观测方程为：

(7)

式中，V为n维测值残差向量；A为n×m阶设计矩阵；a_i为A的第i行；L为对应响应面模型的优化目标，即桥梁模态参数指标；A为对应响应面模型的系数矩阵；为对应响应面模型中的待求解的未知参数。

先验权对角矩阵P：

(8)

M估计采用极大似然估计的思想，一组独立小误差联合出现的概率密度最大，即：

(9)

上式等价为：

(10)

选用函数ρ代替-ln f，方程求解的目标是残差最小，残差可以表示为：

(11)

上式对x求导，并令为0。同时计ϕ(ν_j)=∂ρ/∂ν_i，则有：

(12)

令ϕ(ν_i)/ν_i=w_i(权因子)，p_ii=p_iw_i为权元素，则式(12)可写成：

(13)

因此，式(6)的迭代法方程可表示为：

(14)

式中每次平差的权矩阵为：

(15)

式中，P(v)为权矩阵；j为第j个观测值的权函数；i为第i次平差后第j个值的残差。如果权矩阵P(v)是常数阵，则该M估计过程得到的参数估计值为最小二乘估计。

由此解得参数的抗差M估计为：

(16)

式(15)的求解过程与LS求解方法相同，只是将LS求解方法中的权矩阵换成了不断更新的权矩阵。上式一般通过迭代的方法求解，k+1步的迭代解为：

(17)

式中P^k为不断更新的权矩阵。

基于稳健估计的有限元模型修正的整体流程见图 1。

3 斜拉桥有限元模型修正 3.1 桥梁概况

背景桥梁为(90+143+295+143+90) m双塔双索面钢箱梁斜拉桥，桥梁建成时进行了动力荷载试验，测试了桥梁结构的动力性能，桥梁概况见图 2。

建立了该桥有限元模型，根据实际情况，采用半漂浮体系，主梁及主墩单元采用Beam4单元，斜拉索采用Link8单元建立该桥的空间杆系有限元模型。全桥模型共有单元338个，其中主梁136个，主塔54个，主墩28个，斜拉索120个。采用空间鱼刺模型，所有截面几何特性及质量集中于主梁上(横梁质量采用mass21模拟)，刚臂单元材料特性为10 000倍钢箱梁的弹模。有限元模型的频率和实桥前5阶频率见表 1。为建立合理的有限元模型便于桥梁后期评估服务，以识别的前5阶频率为修正目标，对该有限元模型进行动力模型修正。

3.2 修正参数选取

选取9个影响结构的模态参数作为待选参数并进行灵敏度分析。采用Ansys建模，并结合Matlab混合编程技术，通过编写循环调用命令，自动计算9个参数在多种参数变化工况条件下的分析结果。图 4为9个参数的灵敏度。

图 4表明，主塔弹模、主塔密度、主塔面积、拉索面积对频率的影响整体较小，因此不考虑将它们作为待修正参数。主梁密度与主梁面积对各阶频率的影响程度相近，除第4阶频率外，对其余4阶频率均比较敏感，主梁弹模对第3阶和第5阶频率较为敏感，拉索弹模对第1阶频率最为敏感，拉索密度对第4阶和第5阶频率最为敏感。因此选择以下5个参数作为待修正参数进行中心复合设计，待修正参数设定上下浮动30%，待修正参数取值范围见表 2。

3.3 试验设计及响应面拟合

采用中心复合设计，将5个待修正参数在设定范围内进行组合生成待分析参数组合点，调用编写的自动分析程序，共产生100个样本点。采用F检验法对样本数据进行显著性分析，通过P值判断参数和参数组合的显著性，若P≤0.05，则该项对响应有显著影响，以此决定参与多项式拟合的变量和变量组合。

为探究多项式阶次对响应面模型的精度和计算量的影响，分别计算了参与二阶、三阶和四阶多项式的参数和参数组合的显著性。各回归方程的参数组合P值均小于0.05，表明参数组合均具有较强显著性，各参数组合都应该参与多项式拟合。

通过中心复合设计和显著性分析，依据得到的样本点进行响应面多项式拟合。因篇幅限制，本研究只给出二阶响应面模型的第5阶频率的个别参数的响应面模型，见图 5。

采用R²和RMSE指标对响应面模型精度进行检验。图 6分别为二阶、三阶和四阶响应面模型在不同阶频率的响应面拟合精度。

图 6表明，不同阶次多项式下R²指标基本都接近1，RMSE指标均小于0.008，精度检验结果表明各阶响应面模型精度都较高。但是，不同阶次频率的响应面模型的精度并没有随着多项式的阶次的提高而显著提高，说明提高多项式阶次的方法并不能有效提高响应面的拟合精度。但随着多项式阶次的提升，多项式需求解的未知系数个数大幅增加，这一方面增加了响应面的拟合难度，另一方面增加了响应面优化求解时的矩阵维度，显著增加了响应面优化求解的难度和计算量。

上述分析结果表明，提高响应面多项式的阶次，并不一定能全面提升响应面精度，但会大幅增加计算量。因此，在选用低阶多项式能得到较高精度时，应尽量采用低阶响应面模型。

3.4 响应面函数优化求解

实际桥梁的模态参数与有限元模型存在差异，不同的桥梁结构及监测数据识别得到的模态参数与有限元结果的差异不同，有的阶次差异较大，有的阶次差异较小。不同响应面的求解过程会对这种差异进行不同处理。

对比分析引入的稳健估计法和最小二乘法的优化结果，稳健估计法中的权函数采用Hubber权函数。表 3为在背景桥梁模态参数误差基本均衡的情况下，响应面模型基于最小二乘法和稳健估计法分别得到的优化结果。

表 3优化结果说明，在识别结果无误的条件下，基于最小二乘法和稳健估计法得到相同的修正结果，这是因为稳健估计法在进行优化求解时，通过计算误差大小，不断更新权矩阵，给不同阶次的残差赋予不同的权值。当所有参数误差在正常范围内时，权值函数赋予不同阶次的残差相同的权值，此时，稳健估计法退化为最小二乘法。

为了表示稳健估计法在响应面模型优化过程中的优越性，假设识别得到的二阶桥梁频率为0.512 Hz，与有限元结果的偏差较大。以该频率为目标频率进行响应面模型优化求解，表 4给出最小二乘法和稳健估计法分别得到的修正结果。

基于稳健估计的方程求解时，通过引入不断更新的权矩阵，不断降低误差较大的项的权值，将误差较大项的误差控制在本身阶次内，因此，系统方程求解时没有将该阶次的误差分摊到其他阶次，保证了其他阶次结果的准确性。最小二乘法求解时，将数据异常阶次的误差分配到了其他阶次，使其他阶次误差上升。

稳健估计法能保证当识别结果的差异在合理范围内时与最小二乘法识别得到相同的结果，当识别结果与有限元结果在某阶相差较大时，求解结果保证误差不会传递到其他阶，提高了方程求解的鲁棒性。该过程说明，本研究提出的基于稳健估计的响应面优化求解方法能减小优化目标误差较大时误差对修正结果的干扰风险，能控制误差传播，最大程度保证优化结果的可靠性，可用于实际有限元模型修正过程。

4 结论

(1) 响应面模型的精度并不一定随着多项式阶次的升高而升高，但随着多项式阶次的增高，待求解参数大量增加，增加计算量。二阶响应面模型已经具有较高精度，实际应用中在响应面已具有较高精度时，应避免建立高阶响应面模型。

(2) 采用稳健估计法优化响应面模型，能有效控制目标指标误差引起的优化结果错误。基于稳健估计法的优化方法在目标指标偏差不均衡的情况下能将误差控制在指标误差较大的阶次，能有效控制误差传播，避免其他阶次修正结果受到干扰，保证其他阶修正结果的可靠性。

(3) 当目标指标与有限元结果偏差比较均衡时，稳健估计法退化为最小二乘法，得到与最小二乘法相同的优化结果。背景桥梁修正结果表明，稳健估计法能用于有限元模型修正的过程，提高修正结果的可靠性。