基于改进稳态遗传算法的桥梁有限元模型修正*

秦世强，张亚州，康俊涛

（武汉理工大学土木工程与建筑学院，湖北武汉430070）

由于结构设计参数的不确定性、建模过程简化及误差、结构损伤和材料退化等因素，按照设计图纸建立的结构有限元模型难以代表实际结构，需要根据试验结果对有限元模型进行修正［1-2］。有限元模型修正是通过修正模型中的设计参数，如材料弹性模量、密度、截面惯性矩、边界条件等，不断缩小结构响应的有限元计算值与实测值之间的误差。常用于模型修正的结构响应包括位移、应变、模态频率、振型和频响函数等［3-6］。因此，模型修正可看作成优化问题，其目标函数是由结构响应计算值和实测值之间的残差函数构成。现有的模型修正算法大多是在修正参数的设计空间内，寻找一个使目标函数最小的解（即全局最优解）；然而，结构动力学反问题存在多解属性；考虑到建模、测试误差以及修正参数的物理意义，有限模型修正的局部最优解可能比全局最优解更符合实际情况，即局部最优解也有可能是模型修正的合理解［7-8］。已有文献中得到模型修正单一解的，作者也认为可能存在更合理的解［9-10］。因此，需要同时找到目标函数的全局最优和局部最优解，交由分析者根据实际情况和工程经验，选择一个解作为模型修正最合理的解，这类问题被定义为有限元模型的多解问题。

有限元模型修正多解问题已经得到关注并取得了一些进展；Teughels 等［11］提出一种耦合局部最小化算法用于模型修正，并应用于一座钢筋混凝土简支梁的模型修正；该算法能够同时识别损伤工况设置的两个局部最优解；Zárate 和Caice⁃do［7］将Hook Skip and Jump（HSJ）算法引入模型修正，实现对Bill Emerson 大桥的模型修正；该算法可以同时获取目标函数在设计空间内的多个解，作者选定一个局部最优解为最终修正结果，因为相比全局最优解，该局部最优解更具物理意义和符合实际情况；Caicedo 和Yun［8］提出了一种稳态遗传算法（Steady State Genetic Algorithm， SSGA）用于寻找模型修正问题的多个解，该算法通过一种单一角度算子来搜索可能的局部最优解；Shab⁃bir 和Omenzetter［12-13］利用群体智能算法结合连续小生镜技术获取了模型修正的多个解，并成功应用于一座人行桥的模型修正，作者同样选取了一个局部最优解作为最终的结果。然而，现有的模型修正多解算法往往会出现局部最优解遗漏现象。本文基于文献［8］提出的SSGA 算法，提出一种改进的稳态遗传算法（Improved SSGA， ISSGA），用以同时寻找模型修正中的多个解。它相比SS⁃GA，ISSGA 对目标函数的峰值在设计空间内的分布适应性更好，在保证近似的计算效率前提下，有效地避免了局部最优解的遗漏。

1 模型修正理论

有限元模型修正是通过修正结构设计参数，缩小响应计算值与试验值之间的误差，其目标函数J(x)为：

式中，ri(x)表示第i（i=1，2，3…，m，m 表示模型修正中考虑的结构响应类型数量）类结构响应计算值与试验值之间的残差函数，是待修正参数向量x 的函数。αi表示各类响应残差函数的权重系数。一般可以用响应的计算值和试验值之间的相对误差平方和来构建残差函数，如频率的残差函数可以表示为：

式中，fa，i和ft，i分别表示第i阶频率的计算值和试验值；nf表示模型修正中考虑的频率阶次。考虑到修正参数的物理意义，设定修正过程中的修正参数的上下限，即模型修正可定义为约束优化问题。对这类约束优化问题，由于结构响应与设计参数之间的高度非线性关系，模型修正的目标函数通常在解空间内存在多个局部极值［14］；另一方面，由于结构建模误差和实测误差，局部最优解可能比全局最优解能更好地代表实际结构。因此，同时识别目标函数全局最优解和局部最优解，并交由分析者根据工程经验或先验知识来判定解的合理性，是获得桥梁模型修正合理结果的关键。

2 ISSGA算法

SSGA 用于同时寻找模型修正多个解是由Ca⁃icedo 和Yun［8］首先提出，这种算法的最大优点是计算效率高，这是由于在迭代过程中，SSGA 通过一种可行解更新算子，选出一部分种群及其伴侣解用于遗传，这样在每次迭代过程中，并不是所有种群都需要用来计算目标函数值，从而提高了计算效率。然而，SSGA 算法中可行解更新算子是基于单一角度来控制的，即每次迭代过程中计算空间两点对角点所形成的夹角余弦值来判定这两点是否满足可行解的要求。这种单一角度控制方法的缺点是：当目标函数有多个局部最优解近似的位于同一直线上时，单一角度控制方法容易遗漏可能的可行解。如图1所示，目标函数的等高线在搜索空间内有两个位于同一直线上的可行解，如果只从角点A 出发，计算两个可行解与角点A 的夹角余弦值，则会由于角度过小导致可行解在迭代过程中遗漏。为此，本文提出一种双角度算子来更新可行解，即在B 点处增加一个角度控制算子，来避免可行解在迭代过程中的遗漏。图1所示的情况，尽管可以利用在B 的单一角度完成搜索，但是实际情况中分析者对目标函数的局部最优解分布并无先验知识，局部最优解可能位于搜索空间的对角线上，单一角度方法必然会导致可行解的遗漏。基于双角度算子的ISSGA算法介绍如下。

图1 单一角度及双角度算子示意图Fig.1 The sketch map of single angle and double angle operator

2.1 双角度算子

在每个迭代子步中，首先将遗传算法中的个体按照目标函数进行排序，并在搜索空间内任意取一点作为第一个可行解；然后，在搜索空间内设定两个角点A、B，每个个体及当前可行解与角点A分别形成两个向量，计算这两个向量的夹角余弦值；同理，计算角点B处对应的两个向量的夹角余弦值，只要二者余弦值中有一个小于设定的阈值α，即认为该个体是一个候选的可行解。举例说明如下：如图2 所示，假定已经找到l 个可行解，通过双角度算子寻找第l+1 个可行解的两个判定条件如下：

图2 寻找第l+1个可行解的示意图Fig.2 The sketch map of searching the （l+1）th feasible solution

对于获得的候选可行解，需进一步判定其目标函数值的大小，并与当前的最优可行解（当前子步中目标函数值最小值）进行比较，即

式中，J( ⋅)表示目标函数；θl表示当前最优可行解；θl+1表示候选可行解，是从可行解位置来定义目标函数，β 为候选可行解的最优性阈值，其取值与目标函数有关；图3 示意了阈值β 的物理意义及取值方法；假定一个目标函数在修正参数的设计范围内有三个极值，其中全局最优解对应的目标函数值为a，两个局部最优解对应的目标函数值分别为b 和c，若b/a < β < c/a，则根据式（5），全局最优及局部最优1被定义为可行解；若β > c/a，则全局最优、局部最优1 和2 均被定义为可行解。因此，可行解最优性判定阈值β实际上是通过目标函数值来控制解的数量，当多个极值的目标函数值较为接近时，可以通过β 将他们定义为一个可行解，从而确保选择的可行解具有不同的物理意义。

图3 可行解最优性判定阈值β的物理意义Fig.3 The physical meaning of the threshold β for optimality judgement of the alternative solutions

2.2 可行解的伴侣解

为了获得更优位置的可行解，需要通过交叉、变异操作生成新的子代，即可行解需要定义伴侣解，从而作为父代生成新的子代。在获得可行解集合后，需要找到可行解的伴侣解，每个可行解及其伴侣解将作为父代，并通过交叉、变异等遗传算法基本过程完成整个优化问题。每个可行解均对应一个伴侣解，其选择的依据如下：

2.3 算法流程

得到可行解及其伴侣解之后，实施交叉、变异等不断迭代更新可行解，并评估每次可行解更新后的目标函数值，直至整个迭代过程触发终止准则结束。ISSGA的算法流程如下：

1）设定参数，初始化种群。待设定的参数包括可行解角度阈值α、最优性阈值β、种群数量N；最大迭代次数M；设定好参数后，给每个个体随机取值初始化种群；并根据目标函数给初始化的种群进行排序；

2）更新可行解。在搜索空间内确定双角度算子的两个角点，一般可选搜索空间的边界点。假定第一个个体为第一个可行解；并按照式（3）、（4）、（5）进行判断，只要其满足式（3）、（4）中的一个，并且满足式（5），则认为其为一个可行解；按照类似的方法，不断更新可行解；

3）寻找可行解的伴侣解。按照式（6）、（7）给每个可行解选择其伴侣解；

4）交叉、变异。将可行解及其伴侣解作为父代，实施交叉和变异以生成新的子代；交叉算法可选取均匀交叉算法；变异之后，用新生成的子代替换当前可行解集合中最差的2个可行解（即其对应的目标函数值最大）；

5）评估目标函数值。由于可行解及其伴侣解对应的目标函数值均已经评估过，因此这里只需要评估新替换的可行解；

6）终止准则判定。当触发最大迭代次数M，则迭代终止。

3 测试函数

3.1 测试函数形式

采用Ackley′s Path 函数和Rastrigin′s Path 函数，对比ISGGA 算法和SSGA 算法在解决多解问题中的效果。这两个函数是Matlab 遗传算法及进化算法工具箱GEATbx［15］中的测试函数，广泛用于各种随机搜索算法的搜索精度、稳定性的测试和验证。Ackley′s Path 函数、Rastrigin′s Path 函数的函数形式分别为：

式中，x 为自变量，其定义域分别为［-1.5，1.5］，［-1，1］；e 为自然指数；n 为函数的维度，为了使搜索过程能够图形化显示，取函数的维度n为2。分别利用SSGA 和ISSGA 同时寻找两个测试函数的全局和局部最优解，可以发现：Ackley′s Path 函数共有9 个极小值，其中位于（0，0）处的极小值为全局最优解，对应的函数值为15.0，其余8 个极小值为局部最优解；Rastrigin′s Path 函数有4 个极小值，且对称分布，4 个极值均为全局最优值，对应的函数值为-40.502 5。

3.2 算法参数选取

为了保证SSGA 算法和ISSGA 算法的优化结果具有可对比性，两种算法中的参数保持一致。表1列出了数值算例中的算法参数的取值情况。

种群数量，建议的取值范围为500 ～1 000，取的过高，会导致计算时间增加；取的过低，会导致可行解可能无法找到伴侣解；最大迭代次数建议取值范围为100 ～200，相比标准遗传算法，稳态遗传算法的最大特点是通过较少的迭代即可找到结果；交叉概率和变异概率分别取0.9 及0.1，这是标准遗传算法中的默认值；可行解判定阈值α取0.98，已经十分接近1.0，对应的角度约为11.47°，避免了因角度过大导致的可行解遗漏；可行解最优性判定阈值β 取为1.5，取值依据如下。根据图3 可知，阈值β 的定义是为了避免函数值过于接近的解被认定为不同的局部最优解。对Ack⁃ley′s Path 函数，8 个局部最优解最大的函数值为15.903 2，全局最优解对应的函数值为15.000 0，两者比值为1.06，即只要β 大于1.06，该函数的9个极值均被定义为可行解；对Rastrigin′s Path 函数，4 个极值对称分布，均为全局最优，函数值为-40.502 5，比值为1.0；即只要β 大于1.0，该函数的4个极值均被定义为可行解；因此，在两个数值算例中，阈值β均取为1.5。

表1 数值算例中的算法参数取值Table 1 The values of algorithm parameters in numerical simulation

3.3 结果分析

首先，通过图形直观地显示两种算法的收敛结果。图4 给出了两种算法收敛时Ackley′s Path 函数的等高线图；图中，“o”表示当前迭代步的种群分布情况，“☆”表示已找到的可行解。图4 中，两种算法经过50 次左右的迭代，即可收敛；SSGA识别了5 个极值，ISSGA 识别了全部的9 个解；这表明双角度算子的引入，有效地避免了单一角度下解的遗漏。需要指出的是，SSGA 算法遗漏的解均为局部最优解，这表明SSGA 算法在获取全局最优解方面较为成熟。Rastrigin′s Path 函数的搜索结果见表2，可得到与上述类似的结论。

作为一种随机搜索算法，ISSGA 的搜索精度和稳定性值得关注。将两种算法分别独立运行200次，计算200 次搜索结果的均值及标准差，如表2所示。为便于对比，表2也列出了两个测试函数极值的理论解。可以看出，与理论解相比，两种算法的搜索精度较高；标准差方面，除个别极值点外，ISSGA 算法的的搜索结果标准差更小，表明其稳定性要略高于SSGA 算法；总体来说，两种算法的搜索精度及稳定性均较好。然而，在识别的数量上，ISSGA优势明显。

最后，考查两种算法的计算效率；从表2数据可以看出，ISSGA 所需时间略长，这是由于ISGA采用双角度算子评估可行解，对每个可行解需要多评估一次角度值，因而计算时间略有增加。对于实际工程问题，计算时间大多消耗在计算复杂的目标函数值上；相比之下，消耗在算法本身的时间较少［16］。因此，ISSGA 算法略微增加的时间并不会显著影响实际工程问题的计算效率。

图4 Ackley’s Path函数两种算法的收敛结果Fig.4 The convergence result of the two algorithms for ackley′s path function

表2 两种算法独立运行200次的统计结果Table 2 The statistical results of 200 independent runs of two algorithms

4 工程案例

4.1 混凝土小箱梁桥

混凝土小箱梁是我国市政桥梁最常用的桥型之一。如图5所示，试验研究的预应力混凝土连续小箱梁桥的总长为90 m，共计3跨，每跨30 m；桥面宽12.75 m，布置3 个行车道。由于该类桥型的试验资料较多，限于篇幅，仅介绍与模型修正相关的试验工况。选定第一跨为试验跨，测试截面为第一跨跨中截面。图6给出了桥梁标准横断面及测试截面车辆布置情况，在该测试截面的偏载工况下，共布置6 辆加载汽车，横向3 辆，纵向两排；加载汽车的前轴轴重为60 kN，中后轴轴重均为120 kN，试验过程中，记录4片箱梁梁底的位移（d1～d4）；另外，通过环境振动试验测试桥梁的加速度，并结合随机子空间识别桥梁前4阶频率。表3 给出了位移、频率的试验值和初始模型的计算值。可以看到，未经修正的初始模型计算值与试验值相对误差较大，频率最大相对误差达到13.2%，位移最大相对误差为-28.2%；由于静力位移和频率测试技术成熟，一般具有较高的精度；因此，初始模型修正需要进一步修正，减小建模过程简化和模拟误差，从而更好的代表实际结构。

图5 混凝土箱梁立面图Fig.5 The elevation of the concrete box girder bridge

图6 桥梁横截面及测试截面车辆布置Fig.6 The cross section and the layout of test trucks

4.2 修正参数及目标函数

选择如下7个设计参数进行灵敏度分析：箱梁的弹性模量（E1）、横隔板的弹性模量（E2）、箱梁密度（DEN）、四个支座的竖向刚度（K1-K4）。图7 给出了7 个参数的频率和位移灵敏度。可以看出，频率对E1、DEN 两个参数最为敏感，四个支座刚度的频率灵敏度均较低。相对来说，频率对K1、K2更为敏感；位移对E1、K1、K2较敏感，对其余参数均不敏感。因此，选定E1、DEN、K1、K2 四个参数为修正参数；初始模型中修正参数的取值为3.45×1010Pa、2.5×103kg/m3、8×105N/mm、8×105N/mm。在修正过程中，各参数的上下限为其初始值的±30%。利用前4 阶频率和4 个测点的位移，构建如下的目标函数：

表3 频率、位移的初始模型计算值和试验值对比Table 3 The comparison between analytical and experimental values of the frequencies and displacements

式中，f、d表示频率和位移，上标a、t分别表示计算值和试验值；i、j表示频率和位移的阶次或测点位置。为了提高计算效率，目标函数中的频率和位移计算值由代理模型预测；本文采用多项式响应面作为有限元模型的代理模型，构件响应面的过程主要包括：通过中心复合设计和有限元计算获得修正参数和结构响应的样本、利用F检验测试多项式的交叉项和二次项的显著性、最后利用二阶多项式拟合并进行精度测试。由于多项式响应面用于桥梁模型修正已经有较多工程案例［16］，其理论基础较为成熟，此处不再介绍多项式响应的具体建立过程；仅以一阶频率为例，简要介绍其多项式响应面函数和精度评价指标。构建的一阶频率响应面函数表达式为：

进一步比较式（11）预测值和有限元模型计算值，可以计算出一阶频率响应面精度评价指标R2为0.991；其余频率和位移响应面具有相同精度水平，这表明构建响应面模型具有较高精度，可以替代有限元模型预测频率和位移。

图7 修正参数灵敏度分析Fig.7 Sensitivity Analysis of Updating Parameters

4.3 算法参数选择

分别利用SSGA 和ISSGA 两种算法寻找式（10）所示的目标函数在设计空间内的极小值。两种算法的主要参数确定过程介绍如下：种群数量取1 000；最大迭代次数取200；自变量个数为修正参数的数量，取为4；遗传算法交叉和变异概率分别取0.9 和0.1；上述各参数的取值方法与数值函数相似。

可行解判定阈值α 取为0.98，这样可行解判定的角度值为11.47°，能够避免因角度过大导致的可行解遗漏。可行解的最优性判定阈值β 取为1.8；根据图3 所示的β 的物理意义可知，β 为分析者希望获取的目标函数最大值与目标函数最小值的比值。对于实际工程的目标函数而言，分析者并不知道其在设计空间的极值分布情况，一般可采用标准遗传算法先独立运行50 次，分析50 次的搜索结果对应的目标函数值，进而确定β；对于式（10）的目标函数，通过标准遗传算法50次的独立搜索可知，其最小值在0.14 附近，搜索到的最大值在0.25 附近，两者之比为1.786，因此取β为1.8。

搜索空间角点选择上，SSGA 算法需确定一个角点，取四个参数的下限值，即（2.415，1.75，5.6，5.6），四个参数的单位分别为1010Pa、103kg/m3、105N/mm、105N/mm；ISSGA 算法需确定两个角点，除了上述角点外，另一角点将四个参数分成两组，前两个取下限值，后两个取上限值，即另一个角点为（2.415，1.75，10.4，10.4）。与二维的数值函数类比可知，对于多维问题，将修正参数分成两组，即可按照同样的思路确定搜索空间的角点；这种分组的方式在模型修正多解问题中常常用到，如文献［12］在利用ASCE benchmark 作为算例时，同样将多个修正参数分成两组，从而使多维问题转换为二维问题。

4.4 修正结果分析

采用上述参数，分别利用SSGA 和ISSGA 独立运行100 次，寻找目标函数在设计空间的极小值。 表4 列出了两种算法优化结果的统计分析；考虑到随机搜索算法每次搜索并不能收敛到完全相同的解，在进行统计分析时，若两个极值点对应的4 个修正参数之间相对差异不超过5%，并且目标函数值相对差异不超过5%，则认为两个极值点是同一个解，对其进行合并处理。分析表4 的数据可知： （1） 通过两种算法迭代计算，目标函数值由修正前的0. 479 6 （根据表3 和公式（10） 计算得到）降低至最小0. 137 1，使得修正后的模型计算值和试验值相对误差降低，其中ISS⁃GA 识别的解1 即为全局最优解；（2） 在相同的参数情况下，SSGA 共寻找到3 个极小值，ISSGA 寻找到5 个极小值；对比两者的结果可知，SSGA 遗漏了目标函数在0. 14～0. 17 之间的解，这表明ISSGA 算法中引入的双角度算子能够避免解的遗漏；（3） 相比支座刚度（K1、K2），主梁弹性模量和密度的标准差要更小，表明其稳定性更好；这主要是由于支座刚度的频率和位移灵敏度更小，其变化对频率和位移的影响更小；（4） 在平均单次运行时间上，ISSGA 算法由于引入了双角度算子，运行时间略长，这主要是由于用于算法的时间远小于目标函数的计算时间，即ISSGA 并未显著降低计算效率。综上所述，相比SSGA 算法，ISSGA 算法能保证相似的计算效率，并为分析者提供更多有意义的解，这提高了获得模型修正合理解的可能性。

进一步考查表4中的修正参数修正前后的变化情况，箱梁弹性模量（E1）相对其初始值有比较大幅度的提升，变化范围在18.9% ～30% 之间，表明初始模型对实际结构的刚度模拟过低，这主要是由于初始模型中未考虑桥面板对结构刚度的贡献。此外，靠近支座附近箱梁截面增大的部分也未在初始模型中考虑；箱梁混凝土的密度既有提升的也有降低的，各个解的变化范围在-6.8%～10.7%之间，对于修正后密度提升的解，主要是由于初始模型中未考虑桥面铺装的质量；而对于密度降低的解，主要是由于初始模型结构刚度过低，导致频率偏低，优化过程设置了弹性模量的上限值，因此只有通过降低密度来提高频率值，因此这些解对刚度的贡献约为16.7%；ISSGA 解3 中，箱梁弹性模量在修正后提升为18.9%，与16.7%最为接近，因此，ISSGA的解3被选定为最终的解。

表5列出了频率和位移在修正前后的值和初始模型值，其中修正后的计算值由ISSGA 解3 对应的模型计算得到；经过修正，频率相对误差控制在5%以内，位移相对误差控制在10%以内，表明修正后的模型相比初始模型能够更好的代表实际结构。由于频率识别算法较为成熟，因此取5%为可接受的相对误差上限；考虑到位移现场测试环境，取10%为可接受的相对误差上限；上述关于频率和位移的相对误差范围在文献［18-19］亦有类似表述。

表4 两种算法100次独立运行搜索到的解的统计结果Table 4 The statistical results obtained from 100 independent runs of the two algorithms

表5 模型修正前后频率和位移值的对比Table 5 The comparison of frequency and displacement before and after model updating

针对上述结果分析，需要说明的是：（1）各个修正参数的修正变化幅度反映的是结构整体刚度、质量分布的等效变化情况，并非指各参数的真实变化情况，例如箱梁弹性模量修正后的值提升了18.9%～30%，仅仅表明初始模型的刚度被低估，修正后的参数可以理解为等效弹性模量，而并非指混凝土的弹性模量发生了变化；（2）论文通过改进SSGA 算法，能够同时获得模型修正多个具有不同意义的解，避免局部最优解的遗漏，以供分析者根据实际情况和工程经验进行选择。最终选择哪一个解作为模型修正的最终解，一般可以按照如下原则评估多解：（1）优先选择全局最优解作为最终解；（2）在有检测数据支撑时，评估局部最优解修正参数对应的物理意义，结合实际工程选择局部最优解作为最终解；（3）根据后续分析目的选择最终解。例如利用修正后的模型进行动力分析时，可以选择频率更为接近的局部最优解作为最终解。

5 结 论

论文提出了一种改进的稳态遗传算法用于结构有限元模型修正多解识别，采用双角度算子来寻找可行解，有效地避免了可行解的遗漏现象。通过测试函数及混凝土小箱梁桥的模型修正，论述了算法的原理、参数的确定原则，并验证了ISS⁃GA算法相比SSGA算法的优势。结论如下：

（1） ISSGA 算法利用双角度算子来判定可行解，可同时识别目标函数在设计空间内的多个解，避免了解的遗漏现象；

（2） 测试函数表明：相比于SSGA，ISSGA 识别解的精度、稳定性和计算效率与SSGA 相近，但ISSGA 可成功识别全部9个解，而SSGA 只能识别5个解，体现了双角度算子的优势；

（3） 混凝土箱梁的模型修正表明：对实际工程，可以通过将修正参数分组的方式形成搜索空间，角点确定的原则以分组后的修正参数上下界为依据；ISSGA 算法各参数取值方式与数值模拟近似；所提算法能同时获取模型修正的多个解，为最终确定一个合理的模型修正结果提供了可能。