近年来，我国桥梁建设发展迅速．施工不规范或桥梁老化，会导致桥梁结构出现一定程度的损伤，进而带来安全隐患，因此包含桥梁在内的结构损伤识别越来越受到研究者的重视[1-2]．当前，大多数损伤识别方法通常使用动力测量数据．动力损伤识别操作简单、获得信息量较多，但是不仅受结构刚度的影响，还受质量和阻尼比的影响．静力损伤识别只须考虑结构刚度且静力测试成本低、测量精度高．根据静力测量数据的不同，研究者提出了基于实测静位移、静应变及结合静位移和静应变的损伤识别方法．Lu等[3]利用有限的静力位移来定位梁结构的单元损伤．杨秋伟等[4]结合静力测试位移数据和结构有限元模型的刚度矩阵生成静力残余力向量来确定损伤的位置．Liu等[5]介绍了一种结合准静力应变影响线的时域分析技术来定位梁的损伤．王艺霖等[6]探讨了弯曲正应变和剪切应变在桥梁损伤识别中的应用．Maity等[7]结合应变和位移，使用静力信息作为反向传播神经网络损伤识别方法的可能备选方案，以检测简单悬臂梁的现有损伤．为了获得更多的可用数据，一些研究者还采用了基于静动力结合的识别算法．Wang等[8]提出了一种利用结构固有频率和静态测量位移的变化来识别结构损伤的两阶段识别算法．Lu等[9]试图通过结合噪声固有频率和位移测量来识别损伤．李雪艳等[10]通过应变脉冲在时间区间上积分，提出了应变脉冲响应协方差参数来表征结构状态．当考虑初始模型误差和测量误差的不确定性时，上述方法不再适用．基于概率的结构损伤识别方法通过求解随机逆问题，有望实现工程结构损伤的统计识别．目前，已经发展了许多随机静力损伤识别方法，包括蒙特卡罗模拟方法[11-12]和摄动方法[13-14]等．通常蒙特卡罗模拟用来研究识别结果的随机特性和验证其他随机识别方法的准确性．单德山等[11]通过将模态应变能与蒙特卡罗法相结合提出了一种新的随机损伤识别方法．Hu等[12]利用蒙特卡罗模拟研究测量误差的影响，将测量误差建模为零均值高斯分布．然而，将蒙特卡罗模拟应用于大规模复杂问题时对计算能力和工作量要求都很高．低阶摄动法是解决随机损伤识别的另一种常用而有效的方法．He等[13]使用曲率模态形状和频率的一阶和二阶摄动方程来识别损伤．黄斌等[14]考虑模型误差和测量误差不确定性的影响，利用递推随机有限元方法获取结构单元损伤系数的统计特性．摄动法处理小损伤时效率很高，但处理较大损伤时，往往不能保证损伤识别的准确性．文献[15]提出了基于低阶摄动法的随机梁式结构静力损伤识别方法，将结构可靠度中定义的失效概率用来评估结构损伤概率，但对小损伤识别会出现误判的现象．文献[16]在文献[15]基础上引入了L1正则化方法，减小了对结构小损伤的误判，并通过试验证明了在大损伤情况下其识别效果优于最小二乘法，但精度难以控制．本研究在随机有限元模型框架下，采用了同伦分析方法解决损伤不确定性较大的静力随机识别问题，推导了同时考虑初始模型误差和静力测量误差不确定性的随机静力损伤识别方程，引入同伦分析方法来求解该方程．利用静力凝聚技术和L1正则化解决因测量信息不完整和测量的不确定性引起的不适定性问题．1 基于同伦分析的随机损伤识别1.1　随机静力损伤识别方程对于一个具有N个自由度的梁结构体系，初始状态和损伤状态下的静力荷载是相同的，可以得到以下公式KsXs=KdXd=Fs，(1)式中：Ks和Kd分别为初始状态和损伤状态下的N×N维刚度矩阵；Xs和Xd分别为结构在初始状态和损伤状态下的位移向量；Fs为静荷载向量．梁结构损伤引起结构刚度变化，损伤梁结构的刚度矩阵Kd和位移向量Xd可以分别表示为Kd=Ks+∑i=1nαiKi；(2)Xd=Xs-ΔX，(3)式中：n为梁结构单元的总数；αi为第i个单元的损伤系数，表示第i个单元的刚度的变化比例；Ki为梁结构第i个单元刚度矩阵的N×N维扩阶矩阵；ΔX为梁结构初始状态和损伤状态下位移的差值向量．将式(2)和(3)代入式(1)得到KsXs=(Ks+∑i=1nαiKi)Xd．(4)为识别梁结构的损伤，须建立结构的初始模型．然而，与真实结构相比，初始模型会出现误差(结构参数的不确定性)．对于所提出方法，结构参数可以是任意分布，如对数正态分布．但实际的结构参数是有界的，无论是材料参数还是几何参数，都可以假定是服从有界分布的随机变量ξj(j=1,2,⋯,m)且是完全相关的，m为与初始模型误差相对应的随机变量的数量．在上述假设下，Ks和Xs可以分别表示为随机变量的函数，Ks=Κsa+∑j=1mξjKsb；(5)Xs=Xsa+∑j=1mξjXsb，(6)式中：Κsa和Xsa分别为Ks和Xs的均值；∑j=1mξjKsb和∑j=1mξjXsb分别为Ks和Xs的不确定部分；Ksb和Xsb为随机变量ξj的伴随矩阵．假定随机变量ξ(ξ1,ξ2,...,ξm)完全相关，则有ξj=ξ(j=1,2,...,m)，即：Ks=Κsa+ξKsb；(7)Xs=Xsa+ξXsb．(8)在实际应用中，实测静力测量数据也会存在不确定性．静力测量误差和初始模型误差是相互独立的．静力数据的测量误差可用另一个与ξ相互独立的随机变量η来描述，且η同样满足有界分布如贝塔分布，则Xd可以表示为Xd=Xda+ηXdb，(9)式中：Xda为Xd的确定性部分；ηXdb为Xd的随机部分；Xdb为随机变量η的伴随矩阵．令Ksb=[Ks1，Ks2,⋯,Ksm]，Xsb= [Xs1,Xs2,⋯,Xsm]，考虑到αi同样是关于ξ和η的联合函数，将式(9)代入式(4)，可得    (Κsa+ξKsb)(Xsa+ξXsb)=(Κsa+ξKsb)+∑i=1nαiKi(Xda+ηXdb). (10)式(10)即为同时考虑初始模型误差和静力测量误差不确定性的静力随机损伤识别方程．1.2　随机损伤识别方程的同伦解许多随机有限元方法[17-19]，如广义谱随机有限元方法[17]和随机缩聚方法[18]，都可以用来求解方程(10)．在实际测量过程中，往往会存在一些测量误差的变异性比较大以及结构参数与位移之间往往存在非线性关系的情况，因此使用同伦分析方法来确定方程(10)的解．文献[19-20]中的同伦分析方法可用于建立确定性损伤系数与随机损伤系数之间的同伦关系．为了便于同伦分析，式(10)可以改写为      (R0+ηR1)α=ΚsaXsa-M0+ξ(M1+M2-M3)+ξ2M4-ξηΜ5-ηM6, (11)式中：α=[α1,α2,⋯,αn]；R0=[K1Xda，K2Xda，⋯，KnXda]；R1=[K1Xdb,K2Xdb,⋯,KnXdb]；Kn=Kii=n；M0=KsaXda；M1=KsaXsb；M2=KsbXsa；M3=KsbXda；M4=KsbXsb；M5=KsbXdb；M6=KsaXdb．现在可以重构随机变量ξ，η和参数h，p的新函数Γ(ξ,η,h,p)．当p=0时，Γ(ξ,η,0)=α0；当p=1时，Γ(ξ,η,h,1)=α(ξ,η,h)．α0为确定性损伤系数向量，α0=R0-1ΚsaXsa-M0．然后，构造α0与α(ξ,η,h)之间的关系     (1-p)(R0Γ(ξ,η,h,p)-R0α0)=ph[M0+ηM1Γξ,η,h,p-KsaXsa-M0+T0], (12)式中，h为辅助参数，且h≠0；T0=ξ(M1+M2-M3)+ξ2M4-ξηΜ5-ηM6．方程(12)是零阶变形方程，表示α0和α(ξ,η,h)之间的同伦关系．Γ(ξ,η,h,1)关于参数p的麦克劳林级数，展开式可以表示为     Γ(ξ,η,h,1)=Γ(ξ,η,h,0)+∑k∞Γkξ,η,h,0k!pk, (13)式中：k为同伦级数展开的阶数；Γk(ξ,η,h,0)为Γ(ξ,η,h,1)关于p在p=0处的k阶偏导数，可通过取零阶变形方程的k阶偏导数，再令p=0得到，零阶变形方程对p求一次偏导可得         -(R0Γ-R0α0)+(1-p)R0Γ1(KsaXsa-M0+T0)+ph(R0+ηR1)Γ1, (14)其中Γ为Γ(ξ,η,h,p)的缩写．令p=0，式(14)可以写成     R0Γ1p=0=h[(R0+ηR1)α0-(KsaXsa-M0+T0)]. (15)考虑R0α0-(KsaXsa-M0)=0，则Γ1p=0=-h(ηα1+T1)，式中：α1=-R0-1R1α0；T1=R0-1T0．然后，可以得到Γ1(ξ,η,h,p)．一阶变形方程(14)对p求一次偏导可得     -2R0Γ2+(1-p)R0Γ2=2h(R0+ηR1)Γ1+ph(R0+ηR1)Γ2. (17)式(17)是二阶变形方程，令p=0，整理可得      Γ2p=0=2(-h)(1+h)(ηα1+T1)+2(-h)2(η2α2+ηT2), (18)式中：α2=-R0-1R1α1；T2=R0-1R1T1．同理，可以分别确定函数Γ(ξ,η,h,p)关于p的三阶导数．当p=0时，     Γ3p=0=6(-h)(1+h)2(ηα1+T1)+12(-h)21+h(η2α2+ηT2)+6(-h)3(η3α3+η2T3), (19)式中：α3=-R0-1R1α2；T3=R0-1R1T2．其他高阶变形方程Γkp=0可以通过取零级变形方程相对于参数p在p=0处的k次偏导数以相同的方式获得．将Γkp=0代入式(19)，并设p=1，函数Γ(ξ,η,h,p)的麦克劳林级数展开式可表示为     α(ξ,η,h)=Γ(ξ,η,h,1)=α0+βk,1(h)γ1+βk,2(h)γ2+⋯+βk,j(h)γj+⋯, (19)式中：γj=ηjαj+ηj-1Tj (j≥1)；βk,j(h) (j=1,2，⋯,k)为趋近函数，表示为βk,j(h)=0     (qk);(-h)q∑j=0k-qj+q-1j(1+h)j    (1≤q≤k);1     (q0). (21)趋近函数βk,j(h)的取值可以控制同伦级数的收敛域和收敛速度．在损伤识别过程中，采用有限阶项来提高计算效率．αl(ξ,η,h)可用来表示同伦级数的l阶近似解．对于第i个单元的损伤系数，同伦级数解αi(ξ,η,h)与l阶近似解αl,i(ξ,η,h)之间存在不可避免的近似误差，它们之间的相对误差可定义为ψi(ξ,η,h)=αi(ξ,η,h)-αl,i(ξ,η,h)αi(ξ,η,h)．(22)假设随机变量{ξ,η}的均值和均方差分别为{μ1,μ2}和{σ1,σ2}，则ξφ=μ1+φσ1，ηφ=μ2+φσ2  (φ=1,2,3)，它们是随机变量ξ和η的样本空间中的三个样本．通过求解方程(10)，可以得到第i个单元相对于这三个样本的确定性损伤系数αdi(ξφ,ηφ)．这样求解得到的损伤系数就是精确解．然后，将样本点坐标带入到同伦级数表达式进行计算，可得到第i个单元损伤系数的l阶同伦近似解αl,i(ξφ,ηφ,h) ，该近似解仅依赖于参数h．为了获得h的最优解，将目标函数Wφh定义为Wφh=∑i=1nαdi(ξφ,ηφ)-αl,i(ξφ,ηφ,h)2，(23)可以通过优化目标函数Wφ(h)找到h的适当值，使得误差函数ψi(ξ,η,h)的取值尽可能小．1.3　损伤识别中的静力凝聚在实际结构的损伤识别中，静力测量位移的数据是不完备的．为了消除有限的测量信息对损伤识别的效果产生影响，使用静力凝聚法[21-22]来消除未测量的位移．将位移向量Xs分为竖向位移向量Xsv和转角位移向量Xsθ，将式(1)矩阵中的竖向自由度和转角自由度分离，可表示为KsvvKsvθKsθvKsθθXsvXsθ=FsvFsθ, (24)式中：v和θ分别为测量和未测量位移，对于欧拉-伯努利梁，分别表示竖向位移自由度和转角自由度；Ksvv和Ksvθ为Ks中对应于竖向位移自由度的子矩阵；Ksθv和Ksθθ为Ks对应于转角自由度的子矩阵；Fsv和Fsθ分别为对应于竖向位移自由度和转角自由度的力向量．因实际结构中，通常会在梁结构竖直方向进行加载，故力的转动分量Fsθ=0，代入式(24)得到Xsθ=-Ksθθ-1KsθvXsv，(25)将式(25)代入式(24)可得(Ksvv-KsvθKsθθ-1Ksθv)Xsv=Fsv，(26)为方便表达，上式可简写为KsvXsv=Fsv，(27)式中Ksv=Ksvv-KsvθKsθθ-1Ksθv．同理，可以得到损伤后结构经过静力凝聚的刚度矩阵Kdv，Kdv=Kdvv-KdvθKdθθ-1Kdθv．(28)将式(2)代入(28)，可得      Kdv=Ksvv+∑i=1nαiKivv-Ksvθ+∑i=1nαiKivθ∙Ksθθ+∑i=1nαiKiθθ-1Ksθv+∑i=1nαiKiθv. (29)将式(29)对αi在αi=0处求偏导可得      ∂Kdv∂αiαi=0=Kivv-KiθvKsθθ-1Ksθv-KsθvKsθθ-1Kiθv+KsvθKsθθ-1KiθθKsθθ-1Ksθv, (30)对Kdv在αi=0处进行一阶泰勒展开，可得Kdv=Kdvαi=0+∑i=1nαi∂Kdv∂αiαi=0，(31)将式(31)代入式(30)，可得KsvXsv=Kdvαi=0+∑i=1nαi∂Kdv∂αiαi=0Xdv，(32)式中αi=0表示结构损伤前后的刚度不发生改变，即Kdvαi=0=Ksv．令Kiv=∂Kdv∂αiαi=0，式(32)可以表示为KsvXsv=Ksv+∑i=1nαiKivXdv．(33)式(33)为经过静力凝聚的确定性静力损伤识别方程．考虑初始模型误差和静力测量误差的不确定性，按照前文所示步骤，式(10)可以改写为       ∑i=1nαi(ξ,η)Kiv(Xdav+ηXdbv)=(Κsav+ξKsbv)[(Xsav+ξXsbv)-(Xdav+ηXdbv)], (34)式中：Xsav,Κsav和Xdav分别为静力凝聚后的初始位移矢量Xsv、初始刚度矩阵Κsv和实测位移矢量Xsv的均值部分；Xsbv和Κsbv为ξ的伴随矩阵；Xdbv为η的伴随矩阵．然后，按照1.2节中的步骤，可以得到随机损伤系数α．1.4　L1正则化在实际的损伤识别过程中，测量的不确定性及不完备的测量信息会导致随机损伤识别方程中的损伤系数与方程的数量不匹配，产生随机损伤识别方程中的不适定性问题．引入文献[23]中的L1正则化算法来求解静力随机损伤识别方程．式(34)可以简化为Rα-M的形式，其中，R和M都是已知矩阵．α可以通过最小化目标函数Rα-M来求解．采用L1正则化技术，目标函数被重新定义为H(α)=Rα-M22+λα1，式中：Rα-M22为残差范数；α11为解范数；λ为正则化参数，可以影响正则化算法的约束效果．正则化参数可以由L曲线准则确定．1.5　基于概率的损伤识别为了描述结构的损伤概率，将结构可靠度中的失效概率定义引入到损伤概率的评估中[24]．单元级别的结构损伤概率可定义为：对于第i个单元，初始刚度参数损伤状态下的刚度参数kdi大于ksi的概率或者随机损伤系数αi0的概率．因此，结构中第i个单元的损伤概率Pdi可以写为Pdi=P(ksikdi)=P(αi0)．(36)为提高损伤识别的灵敏度，提出了损伤概率指标Φdi，Φdi=(Pdi-0.5)/0.5，Φdi取值在0~1之间，当Φdi值接近1时，表示单元i发生损伤．当Φdi趋于0时，说明单元i发生损伤的可能性很小．2 数值算例2.1　简支梁静力损伤识别算例假定结构损伤是由于梁单元抗弯刚度的降低引起的，同时考虑模型误差和测量误差．一个矩形变截面的简支梁力学模型示意图如图1所示．每跨梁梁长L=1.9 m，左跨和右跨梁段的截面分别为0.15 m×0.25 m和0.15 m×0.35 m，材料假定为混凝土，强度为C25，弹性模量为28 GPa．简支梁分为12个欧拉-伯努利梁单元和13个节点．每个节点包含一个竖向位移和一个转角．假定在初始状态下，梁单元的抗弯刚度为服从于贝塔分布的随机变量，分别在节点4和节点10施加两个竖向集中荷载F=100 kN，只测量每个节点的竖向位移，利用静力凝聚技术去除节点的转动自由度．10.13245/j.hust.239242.F001图1变截面简支梁力学模型示意图给出一种损伤工况：单元1~12的抗弯刚度折减量分别为5%，10%，20%，25%，20%，5%，5%，10%，20%，25%，20%，5%，用抗弯刚度的折减率来描述单元的损伤程度．每个单元的抗弯刚度折减率定义为实际单元的抗弯刚度与初始完好状态下相同单元的抗弯刚度的相对误差．随机变量均假定为贝塔分布．用随机损伤识别(homotopy damage identification，HDI)、基于一阶摄动法的随机损伤识别(first-order perturbation damage identification，FPDI)和蒙特卡罗(Monte Carlo，MC)方法进行损伤识别，以1×105个样本的MC方法识别的结果为基准．利用Matlab软件建立有限元模型，根据1.3节中静力凝聚方法去除节点的转向自由度．2.2　测量误差的影响为了研究测量误差不确定性对简支梁结构损伤识别结果的影响．假定位移测量误差的变异系数(CV)分别为0.15和0.20，测量误差和建模误差均为贝塔分布，初始模型误差的变异系数为0.1．图2为单元损伤系数均值α¯(CVC=0.15)．可以看出：HDI方法计算的每个单元的损伤系数的均值与MC方法的一致．随着测量误差变异系数的增大，各单元损伤系数的均值或波动都会增大，但相比FPDI方法，HDI方法与MC方法的结果更加符合．10.13245/j.hust.239242.F002图2单元损伤系数均值(CVC=0.15)图3为单元10损伤系数(α10)的概率密度函数(f)．当CVC=0.15时，采用HDI方法得到的损伤系数的概率密度函数曲线与MC方法非常符合，然而HPDI方法结果有不小偏差．当CVC=0.20时，可发现同样现象．10.13245/j.hust.239242.F003图3单元10损伤系数的概率密度函数HDI和MC方法的每个单元的损伤概率结果相差不大，当测量误差的CVC=0.15~0.20时，无法确认第1，2，6，7，8，11和12个单元中假定的小损伤．测量误差变异系数增大，FPDI方法计算得到每个单元的损伤概率也与MC方法相差很远．2.3　模型误差的影响为了说明较高的初始建模误差对损伤识别的影响，假定简支梁初始模型误差的变异系数CVM=0.20，测量误差的变异系数CVC=0.15，使用HDI，FPDI和MC方法的识别结果如图4所示．可以看出：与FPDI方法相比，HDI方法得到的各单元损伤系数的均值和单元10损伤系数的概率密度函数更接近MC方法；与图2相比，当模型误差较大时，HDI方法有着更好的计算结果。图4HDI，FPDI和MC方法的识别结果10.13245/j.hust.239242.F4a110.13245/j.hust.239242.F4a22.4　三种方法的计算效率比较HDI，FPDI和MC方法的计算时间为534，443和4 158 s．FPDI方法的时间消耗约为MC方法的1/10．尽管HDI方法的计算时间比FPDI方法的计算时间慢91 s，但精度明显高于FPDI方法．HDI方法的计算效率远高于MC方法．大型复杂结构的单元数量可能数以万计，HDI方法将会比MC方法的计算效率更高．3 试验验证为验证所提出的识别方法的有效性，进行了混凝土连续梁的静力损伤识别试验．试验符合文献[25]中的测试标准．混凝土连续梁的静载试验中混凝土梁为两跨连续梁，主要材料参数如下：试验梁每跨的长度为2 m，计算长度为1.9 m，钢筋混凝土梁截面为矩形截面，横截面为0.15 m×0.25 m，混凝土强度为C25，混凝土弹性模量为2.8×1010 Pa．梁的力学模型如图5所示，连续梁被分成12单元．假定混凝土梁结构的初始模型不确定性是由于抗弯刚度的不确定性引起的，并且随机性服从贝塔分布．试验中用结构单元抗弯刚度的减少来表示结构单元的损伤．10.13245/j.hust.239242.F005图5混凝土连续梁的力学模型在加载试验中，只测量了梁的竖向挠度．当F=30 kN时，使用千分表多次测量11个节点的竖向挠度，见表1．这里假定所有测量的挠度都是完全相关的，并为贝塔分布的随机变量．根据实测的混凝土弹性模量，假定初始梁单元抗弯刚度的变异系数为0.1．使用HDI和MC方法计算了每个单元的损伤系数的均值和均方差(Eα)，结果如图6所示．可以看出：识别出的损伤单元是第3，4，6，7，9和10单元，它们位于连续梁的跨中区域和中间支座．对于同伦分析方法，6个单元的损伤系数均值为0.188~0.239，均方差为0.10~0.17．第1，2，5，8，11，12单元的损伤系数均值为负值或接近于零，因此下面只讨论第3，4，6，7，9，10单元的损伤情况．此外，HDI方法的统计结果接近MC方法．第3，4，6，7，9和10单元的损伤系数的概率密度函数曲线中，HDI方法的结果始终与MC模拟的结果相符合．6个单元的破坏概率都大于0.5．10.13245/j.hust.239242.T001表1连续梁在30 kN下竖向挠度测量数据节点号竖向挠度/mm变异系数单元1单元2单元3单元4单元5均值2-0.365-0.381-0.416-0.442-0.479-0.4160.1103-0.627-0.673-0.711-0.770-0.852-0.7270.1204-0.661-0.718-0.767-0.821-0.871-0.7680.1085-0.443-0.470-0.509-0.542-0.590-0.5110.1146-0.168-0.185-0.204-0.231-0.256-0.2090.1698-0.100-0.113-0.129-0.149-0.167-0.1320.2059-0.333-0.359-0.389-0.467-0.502-0.4100.17610-0.496-0.533-0.617-0.675-0.744-0.6130.16511-0.424-0.472-0.509-0.568-0.612-0.5170.14512-0.262-0.274-0.308-0.334-0.363-0.3080.13510.13245/j.hust.239242.F006图6两种不同方法中单元损伤系数的均值和均方差梁在不同载荷下的裂缝状态如图7所示．当F=30 kN时，蓝色部分出现裂纹．梁高方向的裂缝长度差异不大，即第3，4，6，7，9和10对应的梁单元损伤程度基本相同．这与所提出的同伦分析方法的识别结果是一致的，说明了同伦分析方法的有效性．10.13245/j.hust.239242.F007图7混凝土连续梁裂缝状态4 结论a．本文方法适用于不同损伤程度的结构．但是，若测量误差和建模误差的不确定度较高，则会导致无法识别小损伤．b．当测量误差和建模误差较大时，本文方法比FPDI方法具有更高的损伤识别精度．c．在保证损伤识别精度的前提下，HDI方法的计算效率高于MC方法，约为MC所耗时长的13%．d．本文方法对连续梁的识别结果与静力试验中观察到的现象一致．