桥梁极值应力的贝叶斯动态耦合线性预测

摘要

基于含有趋势性和随机性的监测极值应力时间序列数据，进行了桥梁极值应力的贝叶斯动态耦合线性预测。将历史监测耦合极值应力数据解耦，建立解耦极值应力的动态耦合线性模型（DCLM)，并将其分解为有限个动态线性模型(DLMs)。基于动态监测的解耦极值应力数据，结合贝叶斯方法，进行动态耦合线性模型的概率递推，进而实现桥梁解耦极值应力的动态预测，并将预测的解耦极值应力进行相加，即可实现桥梁极值应力的动态预测。最后，通过在役桥梁的含有随机性和趋势性的监测极值应力数据，利用ARIMA（auto-regressive integrated moving average）模型，验证了所提方法的合理性与有效性。

关键词

桥梁结构; 监测数据; 动态耦合线性模型(DCLM); 动态线性模型(DLM); 贝叶斯方法; 极值应力预测

桥梁健康监测(BHM)系统在长期运营过程中积累了大量的、全方位的监测信息，找寻隐藏在这些信息中的特定规律，并根据所得规律来预测信息的变化趋势，进而预估桥梁结构的状态，及早规划桥梁的养护策略，可以有效避免事故的发生。发展至今，桥梁健康监测研究主要经历了数据采集和监测数据的合理应用两个阶段^[1]。数据采集阶段主要是传感器和结构健康监测系统的研制，从而获得动态监测数据，数据采集阶段的大量研究主要集中在数据传送系统^[2]、数据压缩系统^[3]、数据恢复系统^[4]、数据获得技术^[5]、测点布置和系统组装技术^[6]等，侧重于数据采集方法、数据实时性和准确性等方面，目前已处于相对成熟阶段。监测数据的合理应用阶段，目前大量研究主要集中在车辆荷载建模和荷载效应分离^[7]、结构模态参数识别^[8]、结构模型修正^[9]、结构损伤识别^[10]、结构可靠性评估及动态预测^[11]等领域，但仍难以结合监测数据对在役桥梁极值应力进行合理动态预估。由于桥梁可靠性的度量主要依赖于结构抗力与荷载效应，因此荷载效应极值的动态预测可以为桥梁可靠性预测提供理论基础。考虑到监测数据的随机性和趋势性等特性，研究在役桥梁极值应力的动态预测就成为桥梁健康监测的关键问题之一。

对于在役桥梁结构，鉴于有限元模型精细化建模复杂、模型动态更新困难以及智能算法耗时等特点，从动态监测数据的有模型极值应力预测方法到实际的高精度应用还需要大量的研究，尤其是如何采用含有随机性和趋势性等特性的监测数据直接预测桥梁极值应力。因此，有模型方法的特点使得无模型方法逐渐成为桥梁极值应力预测的发展趋势。

采用无有限元模型的分析方法，基于已有监测信息建立合理的数学模型，给出模型参数与监测信息之间的动态变化关系，进而结合动态监测数据，可实现桥梁信息的动态预测。国内外学者已经给出了一些无模型分析方法，如基于数据驱动的随机子空间算法^[12]、基于主元分析的故障诊断算法^[13]等，这些方法主要应用于桥梁结构模态参数识别与损伤识别，无法进行动态预测。鉴于上述算法的不足之处，基于贝叶斯动态模型的分析方法具有独特的先进性，能够实现桥梁极值应力的动态预测。目前，桥梁极值应力预测已取得一些研究成果。对于健康监测条件下基于混合高斯粒子滤波器的桥梁极值应力贝叶斯动态线性预测^[1]，分析过程中假定状态差值为定值，不考虑动态性和随机性；对于基于均匀监测数据的桥梁极值应力贝叶斯动态非线性预测方法^[14]，分析过程中所建立的状态方程依赖于历史数据的回归方程，没有普遍性；对于基于非均匀监测数据的桥梁动力响应贝叶斯动态线性离线预测方法^[15]，分析过程中未考虑状态的多峰分布形式；对于动力响应的解耦及其预测方法^[16]，分析过程中耦合信息未实现高精度解耦。鉴于上述已有桥梁极值应力预测方法存在的缺点，考虑监测数据的随机性、动态性、趋势性以及状态多样性等特点，建立适合于普遍应用的动态耦合线性模型(DCLM)，并结合贝叶斯方法实现桥梁极值应力的动态预测，已成为桥梁健康监测预防性维修决策的关键所在。

对于基于多种荷载效应耦合的含有趋势性和随机性的监测极值应力时间序列数据，首先采用简单移动平均法进行解耦，分别对各种荷载效应下的极值应力数据建立动态耦合线性模型，并将其分解为有限个动态线性模型(DLM)。然后，对于基于动态监测的极值应力数据，结合贝叶斯方法，进行动态耦合线性模型的概率递推，进而实现桥梁解耦极值应力的动态预测，将预测的解耦极值应力进行相加，即可得到桥梁极值应力的预测值。最后，通过在役桥梁的监测极值应力数据对所提模型和方法的合理性与有效性进行验证。

1 桥梁极值应力的动态耦合线性模型及对应动态线性模型

1.1　桥梁极值应力的动态耦合线性模型

桥梁健康监测系统在长期运营过程中积累了大量监测应力信息，应力是在温度荷载、车辆荷载以及桥梁自重恒载的耦合作用下产生的，具有随机性、动态性以及趋势性等特点。采用简单移动平均法来进行应力监测数据时间序列的趋势项提取，将原始时程曲线与趋势项曲线的差值作为车辆荷载引起的极值应力；将趋势项的均值作为由自重恒载作用引起的极值应力，并认为由自重恒载引起的极值应力在监测的时间长度内不随时间变化；将剔除均值后趋势项的剩余部分作为温度引起的极值应力。通过上述过程实现了车辆荷载效应、温度荷载效应和恒载效应的解耦。

针对解耦的车辆荷载效应与温度荷载效应数据，经过五点三次平滑处理之后的数据为极值应力状态数据^[1]，并定义 $\left\{X_t,t=\mathrm{1,2},\cdots \right\}$ 为极值应力状态时间序列， $X_t$ 为t时刻的极值应力状态。所建动态耦合线性模型基于以下两点假设：

(1) $\left\{X_t,t=\mathrm{1,2},\cdots \right\}$ 具有马尔科夫性。

(2) 监测应力 $\left\{Z_t,t=\mathrm{1,2},\cdots \right\}$ 之间相互独立，并且 $Z_t$ 仅和状态变量 $X_t$ 相关， $Z_t$ 与 $X_t$ 成线性关系。

桥梁监测解耦极值应力数据存在一定程度的非线性特征，采用非线性动态模型会得到更好的预测效果，但是非线性模型存在过程复杂、耗时长等问题。因此，对于弱非线性数据，可以将其转化为线性问题，从而能够使用线性预测模型处理^[1]。若监测极值应力的变化趋势曲线能够用一个简单的二次曲线 $h(t)$ 表示，并将其作为结构变量的状态变化曲线，则可以产生以下推导：

$$X_t=h(t)$$，$$X_{t+\mathrm{1}}=h(t+\mathrm{1})$$

采用泰勒级数展开将 $h(t+\mathrm{1})$ 在 $t$ 处展开，可得

$$X_{t+\mathrm{1}}\approx X_t+h\text{'}(t)((t+\mathrm{1})-t)=X_t+h\text{'}(t)$$

式中： $h\text{'}(t)$ 为 $h(t)$ 在 $t$ 时刻的一阶差分。

对于所得的极值应力状态数据，由Augmented Dickey‒Fuller方法^[17]检验可知：状态时间序列为非平稳时间序列，而状态时间序列的一阶差分时间序列为平稳时间序列，得到的状态方程所用随机变量包括状态变量和状态的一阶差分变量，进而可得桥梁解耦极值应力的动态模型。动态模型的状态方程如下所示：

$$\begin{array}{c}{X}_t=X_{t-\mathrm{1}}+{\beta }_{t-\mathrm{1}}+w_{t,\mathrm{1}},{\beta }_t={\beta }_{t-\mathrm{1}}+w_{t,\mathrm{2}},{\mathit{w}}_t~N(\mathrm{0},{\mathit{W}}_t)\\ \Rightarrow {\mathit{\theta }}_t=\mathit{G}{\mathit{\theta }}_{t-\mathrm{1}}+{\mathit{w}}_t,{\mathit{w}}_t~N(\mathrm{0},{\mathit{W}}_t)\\ {\mathit{\theta }}_t=（X_t,{\beta }_t{）}^{\mathrm{T}},{\mathit{w}}_t=（w_{t,\mathrm{1}},w_{t,\mathrm{2}}{）}^{\mathrm{T}},\mathit{G}=\left[\begin{array}{cc}\mathrm{1}& \mathrm{1}\\ \mathrm{0}& \mathrm{1}\end{array}\right],{\mathit{W}}_t=\left[\begin{array}{cc}{W}_{t,\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& W_{t,\mathrm{2}}\end{array}\right]\end{array}$$

(1)

动态模型的监测方程为

$$Z_t=X_t+v_t,v_t~N(\mathrm{0},V_t)$$

(2)

动态模型的初始状态信息为

$$\begin{array}{c}p(X_{t-\mathrm{1}}\left|{\mathit{D}}_{t-\mathrm{1,1}}\right.)\approx {\alpha }_{i_{\mathrm{1}},t-\mathrm{1}}N(X_{t-\mathrm{1}};{\mu }_{i_{\mathrm{1}},t-\mathrm{1}},{\varphi }_{i_{\mathrm{1}},t-\mathrm{1}})\\ p({\beta }_{t-\mathrm{1}}\left|{\mathit{D}}_{t-\mathrm{1,2}}\right.)\approx {\alpha }_{i_{\mathrm{2}},t-\mathrm{1}}N({\beta }_{t-\mathrm{1}};{\mu }_{i_{\mathrm{2}},t-\mathrm{1}},{\varphi }_{i_{\mathrm{2}},t-\mathrm{1}})\end{array}$$

(3)

式(1)~(3)中： $v_t$ 为监测误差； $V_t$ 为监测误差方差； $w_{t,\mathrm{1}}$ 为状态误差； $w_{t,\mathrm{2}}$ 为状态一阶差分误差； ${\mathit{w}}_t$ 为包括 $w_{t,\mathrm{1}}$ 和 $w_{t,\mathrm{2}}$ 的向量； ${\mathit{W}}_t$ 为 ${\mathit{w}}_t$ 的协方差矩阵； $W_{t,\mathrm{1}}$ 和 $W_{t,\mathrm{2}}$ 为状态误差方差； $p(X_{t-\mathrm{1}}\left|{\mathit{D}}_{t-\mathrm{1,1}}\right.)$ 、 $p({\beta }_{t-\mathrm{1}}\left|{\mathit{D}}_{t-\mathrm{1,2}}\right.)$ 为概率密度函数，可根据极值应力状态时间序列样本及其一阶差分样本的概率统计得到； $N(·)$ 为正态概率密度函数(PDF)； ${\alpha }_{i_{\mathrm{1}},t-\mathrm{1}}$ 和 ${\alpha }_{i_{\mathrm{2}},t-\mathrm{1}}$ 分别为对应正态概率密度函数的权重； ${\mu }_{i_{\mathrm{1}},t-\mathrm{1}}$ 和 ${\mu }_{i_{\mathrm{2}},t-\mathrm{1}}$ 分别为对应正态概率密度函数的均值； ${\varphi }_{i_{\mathrm{1}},t-\mathrm{1}}$ 和 ${\varphi }_{i_{\mathrm{2}},t-\mathrm{1}}$ 分别为对应正态概率密度函数的方差； ${\mathit{D}}_{t-\mathrm{1,1}}$ 为 $(t-\mathrm{1})$ 时刻之前状态信息的向量； ${\mathit{D}}_{t-\mathrm{1,2}}$ 为 $(t-\mathrm{1})$ 时刻之前状态一阶差分信息的向量; ${\mathit{\theta }}_t$ 为状态变量； ${\beta }_t$ 为状态变化量。

1.2　动态耦合线性模型对应的动态线性模型

基于式(1)~(3)，将动态耦合线性模型转化为g ₁×g ₂ 个动态线性模型，每个动态线性模型的状态方程和监测方程分别与式(1)和式(2)相一致，动态线性模型与动态耦合线性模型的不同之处在于初始状态信息的不同。动态线性模型的初始状态信息如下所示：

$$\begin{array}{c}p(X_{t-\mathrm{1}}\left|{\mathit{D}}_{t-\mathrm{1,1}}\right.)\approx N(X_{t-\mathrm{1}};{\mu }_{i_{\mathrm{1}},t-\mathrm{1}},{\varphi }_{i_{\mathrm{1}},t-\mathrm{1}}),i_{\mathrm{1}}\in \left\{\mathrm{1,2},\cdots ,g_{\mathrm{1}}\right\}\\ p({\beta }_{t-\mathrm{1}}\left|{\mathit{D}}_{t-\mathrm{1,2}}\right.)\approx N({\beta }_{t-\mathrm{1}};{\mu }_{i_{\mathrm{2}},t-\mathrm{1}},{\varphi }_{i_{\mathrm{2}},t-\mathrm{1}}),i_{\mathrm{2}}\in \left\{\mathrm{1,2},\cdots ,g_{\mathrm{2}}\right\}\\ \Rightarrow (\left.{\mathit{\theta }}_{t-\mathrm{1},J}\right|{\mathit{D}}_{t-\mathrm{1},J})~N({\mathit{M}}_{t-\mathrm{1},J},{\mathit{C}}_{t-\mathrm{1},J}),J\in \left\{\mathrm{1,2},\cdots ,g_{\mathrm{1}}\times g_{\mathrm{2}}\right\}\\ {\mathit{\theta }}_{t-\mathrm{1},J}=(X_{t-\mathrm{1}},{\beta }_{t-\mathrm{1}}{)}_J^{\mathrm{T}}，{\mathit{D}}_{t-\mathrm{1},J}=({\mathit{D}}_{t-\mathrm{1,1}},{\mathit{D}}_{t-\mathrm{1,2}}{)}_J^{\mathrm{T}}，{\mathit{M}}_{t-\mathrm{1},J}=({\mu }_{i_{\mathrm{1}},t-\mathrm{1}},{\mu }_{i_{\mathrm{2}},t-\mathrm{1}}{)}_J^{\mathrm{T}}，{\mathit{C}}_{t-\mathrm{1},J}={\left[\begin{array}{cc}{\varphi }_{i_{\mathrm{1}},t-\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& {\varphi }_{i_{\mathrm{2}},t-\mathrm{1}}\end{array}\right]}_J\end{array}$$

(4)

2 基于贝叶斯方法的动态耦合线性模型概率递推

首先基于式(1)、(2)和式(4)，采用贝叶斯方法^[18]对每个动态线性模型进行动态概率递推，然后结合动态线性模型的动态组合权重，实现动态耦合线性模型的概率递推。

2.1　动态线性模型的贝叶斯动态概率递推

(1) 利用式 (1) 和 式(4)，可得到状态 ${\mathit{\theta }}_{t,J}$ 的先验概率密度函数为

$$({\mathit{\theta }}_{t,J}\left|{\mathit{D}}_{t-\mathrm{1},J}\right.)~N({\mathit{a}}_{t,J},{\mathit{R}}_{t,J}),J\in \left\{\mathrm{1,2},\cdots ,g_{\mathrm{1}}\times g_{\mathrm{2}}\right\}$$

(5)

$${\mathit{a}}_{t,J}=\mathit{G}{\mathit{M}}_{t-\mathrm{1},J}$$，$${\mathit{R}}_{t,J}=\mathit{G}{\mathit{C}}_{t-\mathrm{1},J}{\mathit{G}}^{\mathrm{T}}+{\mathit{W}}_{t,J}$$

(2) 利用式 (2) 和式 (5)，可得监测变量的一步预测概率密度函数为

$$(\left.y_{t,J}\right|{\mathit{D}}_{t-\mathrm{1},J})~N(f_{t,J},Q_{t,J}),J\in \left\{\mathrm{1,2},\cdots ,g_{\mathrm{1}}\times g_{\mathrm{2}}\right\}$$

(6)

$$f_{t,J}=\mathit{F}{\mathit{a}}_{t,J},Q_{t,J}=\mathit{F}{\mathit{R}}_{t,J}{\mathit{F}}^{\mathrm{T}}+V_{t,J},V_{t,J}=V_t$$

式中： F 为系数矩阵; $y_{t,J}$ 为t时刻的监测数据。因此，t时刻第J个动态线性模型的含有95%置信率的预测区间^[18]为 $\left[f_{t,J}-\mathrm{1.645}\sqrt[]{Q_{t,J}},f_{t,J}+\mathrm{1.645}\sqrt[]{Q_{t,J}}\right]$ 。

(3) 综合利用贝叶斯方法^[18]和式(5)、(6)，基于y_t,J =y_t ，状态 ${\mathit{\theta }}_{t,J}$ 的后验概率密度函数为

$$({\mathit{\theta }}_{t,J}\left|{\mathit{D}}_{t,J}\right.)~N({\mathit{M}}_{t,J},{\mathit{C}}_{t,J}),J\in \left\{\mathrm{1,2},\cdots ,g_{\mathrm{1}}\times g_{\mathrm{2}}\right\}$$

(7)

$${\mathit{M}}_{t,J}={\mathit{a}}_{t,J}+{\mathit{A}}_{t,J}{e}_{t,J},{\mathit{C}}_{t,J}={\mathit{R}}_{t,J}-{\mathit{A}}_{t,J}{Q}_{t,J}{\mathit{A}}_{t,J}^{\mathrm{T}},{\mathit{A}}_{t,J}={\mathit{R}}_{t,J}{\mathit{F}}^{\mathrm{T}}{Q}_{t,J}^{-\mathrm{1}},e_{t,J}=y_{t,J}-f_{t,J}$$

式中： $e_{t,J}$ 为一步预测误差； $Q_{t,J}^{-\mathrm{1}}$ 为第J个动态线性模型的预测精度。

基于式(1)和式(2)以及式(4)~(7)，即可以得到每个动态线性模型的循环动态概率修正递推过程。

2.2　动态耦合线性模型的贝叶斯动态概率递推

利用式(1)~(7)，并结合动态监测的解耦极值应力数据，可得含有时变权重的动态耦合线性模型的贝叶斯动态概率修正过程。

(1) 状态 ${\mathit{\theta }}_{t-\mathrm{1}}$ 的后验概率密度函数为

$$({\mathit{\theta }}_{t-\mathrm{1}}\left|{\mathit{D}}_{t-\mathrm{1}}\right.)~\stackrel{g_{\mathrm{1}}\times g_{\mathrm{2}}}{\sum _{J=\mathrm{1}}}(p_{t-\mathrm{1},J}N({\mathit{M}}_{t-\mathrm{1},J},{\mathit{C}}_{t-\mathrm{1},J}))$$

(8)

式中： $p_{t-\mathrm{1},J}$ 由式(3)、(4)来确定。

(2) 状态 ${\mathit{\theta }}_t$ 的先验概率密度函数为

$$({\mathit{\theta }}_t\left|{\mathit{D}}_{t-\mathrm{1}}\right.)~\stackrel{g_{\mathrm{1}}\times g_{\mathrm{2}}}{\sum _{J=\mathrm{1}}}(p_{t-\mathrm{1},J}N({\mathit{a}}_{t,J},{\mathit{R}}_{t,J}))$$

(9)

(3) t时刻的一步向前预测概率密度函数为

$$(\left.y_t\right|{\mathit{D}}_{t-\mathrm{1}})~\stackrel{g_{\mathrm{1}}\times g_{\mathrm{2}}}{\sum _{J=\mathrm{1}}}(p_{t-\mathrm{1},J}N(f_{t,J},Q_{t,J}))$$

(10)

可得预测值、预测精度分别为

$$f_t=\stackrel{g_{\mathrm{1}}\times g_{\mathrm{2}}}{\sum _{J=\mathrm{1}}}(p_{t-\mathrm{1},J}{f}_{t,J})$$

$$Q_{t,J}={\left((y_t-f_t{)}^{\mathrm{2}}\stackrel{g_{\mathrm{1}}\times g_{\mathrm{2}}}{\sum _{J=\mathrm{1}}}(p_{t-\mathrm{1},J}N(f_{t,J},Q_{t,J}))\mathrm{d}{y}_t\right)}^{-\mathrm{1}}$$

含有95%置信率的预测区间为

$$\left[\stackrel{g_{\mathrm{1}}\times g_{\mathrm{2}}}{\sum _{J=\mathrm{1}}}(p_{t-\mathrm{1},J}{f}_{t,J})-\mathrm{1.645}\sqrt[]{\stackrel{+\mathrm{\infty }}{\underset{-\mathrm{\infty }}{\int }}(y_t-f_t{)}^{\mathrm{2}}\stackrel{g_{\mathrm{1}}\times g_{\mathrm{2}}}{\sum _{J=\mathrm{1}}}(p_{t-\mathrm{1},J}N(f_{t,J},Q_{t,J}))\mathrm{d}{y}_t}，\stackrel{g_{\mathrm{1}}\times g_{\mathrm{2}}}{\sum _{J=\mathrm{1}}}（p_{t-\mathrm{1},J}{f}_{t,J}）+\mathrm{1.645}\sqrt[]{\stackrel{+\mathrm{\infty }}{\underset{-\mathrm{\infty }}{\int }}(y_t-f_t{)}^{\mathrm{2}}\stackrel{g_{\mathrm{1}}\times g_{\mathrm{2}}}{\sum _{J=\mathrm{1}}}(p_{t-\mathrm{1},J}N(f_{t,J},Q_{t,J}))\mathrm{d}{y}_t}\right]$$

(4) 综合利用贝叶斯方法和式(5)、(6)，基于监测值y_t ，状态 ${\mathit{\theta }}_t$ 的后验概率密度函数为

$$({\mathit{\theta }}_t\left|{\mathit{D}}_t\right.)~\stackrel{g_{\mathrm{1}}\times g_{\mathrm{2}}}{\sum _{J=\mathrm{1}}}{p}_{t,J}N({\mathit{M}}_{t,J},{\mathit{C}}_{t,J})$$

(11)

$$p_{t,J}=\frac{p_{t-\mathrm{1},J}p(\left.y_{t,J}\right|{\mathit{D}}_{t-\mathrm{1},J})}{\stackrel{g_{\mathrm{1}}\times g_{\mathrm{2}}}{\sum _{J=\mathrm{1}}}(p_{t-\mathrm{1},J}p(\left.y_{t,J}\right|{\mathit{D}}_{t-\mathrm{1},J}))}=\frac{p_{t-\mathrm{1},J}N(f_{t,J},{\mathit{Q}}_{t,J})}{\stackrel{g_{\mathrm{1}}\times g_{\mathrm{2}}}{\sum _{J=\mathrm{1}}}(p_{t-\mathrm{1},J}N(f_{t,J},{\mathit{Q}}_{t,J}))}$$

3 动态耦合线性模型参数的确定方法

动态耦合线性模型需要确定的参数主要有 $V_t$ 、 $W_{t,\mathrm{1}}$ 、 $W_{t,\mathrm{2}}$ 、 $({\mu }_{i_{\mathrm{1}},t-\mathrm{1}},{\varphi }_{i_{\mathrm{1}},t-\mathrm{1}})$ 、 $({\mu }_{i_{\mathrm{2}},t-\mathrm{1}},{\varphi }_{i_{\mathrm{2}},t-\mathrm{1}})$ 、 ${\alpha }_{i_{\mathrm{1}},t-\mathrm{1}}$ 以及 ${\alpha }_{i_{\mathrm{2}},t-\mathrm{1}}$ 。

监测误差方差V_t 可通过监测数据五点三次平滑处理后的数据与监测数据之间差值序列的方差估计得到；由于状态变量的不可监测性，状态误差方差 $W_{t,\mathrm{1}}$ 、 $W_{t,\mathrm{2}}$ 由初始状态信息的方差结合折扣因子来近似确定，如下所示：

$${\mathit{W}}_{t,J}=-\mathit{G}{\mathit{C}}_{t-\mathrm{1},J}{\mathit{G}}^{\mathrm{T}}+\frac{{\mathit{C}}_{t-\mathrm{1},J}}{\delta }$$

(12)

式中： ${\mathit{W}}_{t,J}$ 为第J个动态线性模型的状态误差和状态一阶差分误差所估计的方差矩阵； $\delta$ 为折扣因子，一般取0.48~0.98^[11]，此处取值为0.98。

对于初始状态信息和状态一阶差值信息的均值 ${\mu }_{i_{\mathrm{1}},t-\mathrm{1}}\mathrm{、}{\mu }_{i_{\mathrm{2}},t-\mathrm{1}}$ 与方差 ${\varphi }_{i_{\mathrm{1}},t-\mathrm{1}}\mathrm{、}{\varphi }_{i_{\mathrm{2}},t-\mathrm{1}}$ ，可以对平滑后的数据及其一阶差分序列进行处理，即可得到初始状态信息和状态一阶差分信息的概率分布。

(1) 由核密度估计方法^[19]近似得到平滑后数据（状态信息）实际概率密度函数的点估计为 $f(X_{t-\mathrm{1}})$ ，可以利用高斯函数^[20]进行g ₁阶拟合并依据最小二乘原理确定参数，如下所示：

$$f(X_{t-\mathrm{1}})=a_{i_{\mathrm{1}}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\left(\frac{X_{t-\mathrm{1}}-b_{i_{\mathrm{1}}}}{c_{i_{\mathrm{1}}}}\right)}^{\mathrm{2}}\right)$$

(13)

式中： $a_{i_{\mathrm{1}}}$ 为对应单峰分布的峰值； $b_{i_{\mathrm{1}}}$ 为对应单峰分布的均值； $c_{i_{\mathrm{1}}}$ 是与峰宽相关的值。

(2) 假定 $c_{i_{\mathrm{1}}}=\sqrt[]{\mathrm{2}}{\sigma }_{i_{\mathrm{1}}}$ ，对式(13)的右项进行变换，变换后计算式如下所示： $\stackrel{g_{\mathrm{1}}}{\sum _{i_{\mathrm{1}}=\mathrm{1}}}{a}_{i_{\mathrm{1}}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\left(\frac{X_{t-\mathrm{1}}-b_{i_{\mathrm{1}}}}{c_{i_{\mathrm{1}}}}\right)}^{\mathrm{2}}\right)=\stackrel{g_{\mathrm{1}}}{\sum _{i_{\mathrm{1}}=\mathrm{1}}}\sqrt[]{\mathrm{2}\mathrm{\pi }}{\sigma }_{i_{\mathrm{1}}}{a}_{i_{\mathrm{1}}}\frac{\mathrm{1}}{\sqrt[]{\mathrm{2}\mathrm{\pi }}{\sigma }_{i_{\mathrm{1}}}}exp\left(-{\left(\frac{X_{t-\mathrm{1}}-b_{i_{\mathrm{1}}}}{\sqrt[]{\mathrm{2}}{\sigma }_{i_{\mathrm{1}}}}\right)}^{\mathrm{2}}\right) =$

$$\stackrel{g_{\mathrm{1}}}{\sum _{i_{\mathrm{1}}=\mathrm{1}}}\sqrt[]{\mathrm{2}\mathrm{\pi }}{\sigma }_{i_{\mathrm{1}}}{a}_{i_{\mathrm{1}}}\varphi \left(\frac{X_{t-\mathrm{1}}-b_{i_{\mathrm{1}}}}{{\sigma }_{i_{\mathrm{1}}}}\right)$$

(14)

式中： ${\sigma }_{i_{\mathrm{1}}}$ 为标准差， ${\sigma }_{i_{\mathrm{1}}}=c_{i_{\mathrm{1}}}/\sqrt[]{\mathrm{2}}$ ； $\varphi (·)$ 为标准正态概率密度函数。由式(14)即可得到n个正态概率密度函数的均值 $b_{i_{\mathrm{1}}}$ 与标准差 ${\sigma }_{i_{\mathrm{1}}}$ 。

(3) 对高斯函数式中的各项分别进行积分，可将各项积分值近似作为各项分布的权重，如下所示：

$$s_{i_{\mathrm{1}}}={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\sqrt[]{\mathrm{2}\mathrm{\pi }}{\sigma }_{i_{\mathrm{1}}}{a}_{i_{\mathrm{1}}}\varphi \left(\frac{X_{t-\mathrm{1}}-b_{i_{\mathrm{1}}}}{{\sigma }_{i_{\mathrm{1}}}}\right)\mathrm{d}x=\sqrt[]{\mathrm{2}\mathrm{\pi }}{\sigma }_{i_{\mathrm{1}}}{a}_{i_{\mathrm{1}}}$$

(15)

式中： $s_{i_{\mathrm{1}}}$ 为对应项的积分值。

(4) 考虑到拟合误差以及计算误差，归一化权重，即可得 ${\alpha }_{i_{\mathrm{1}},t-\mathrm{1}}$ ，如下所示：

$${\alpha }_{i_{\mathrm{1}},t-\mathrm{1}}=\frac{s_{i_{\mathrm{1}}}}{\stackrel{g_{\mathrm{1}}}{\sum _{i=\mathrm{1}}}{s}_{i_{\mathrm{1}}}}$$

(16)

(5) 将g ₁个正态概率密度函数加权组合，即可得到初始状态信息的概率分布,如下所示：

$$（X_{t-\mathrm{1}}）=\stackrel{g_{\mathrm{1}}}{\sum _{i_{\mathrm{1}}=\mathrm{1}}}{\alpha }_{i_{\mathrm{1}},t-\mathrm{1}}N(b_{i_{\mathrm{1}}},{\sigma }_{i_{\mathrm{1}}}^{\mathrm{2}})$$

(17)

基于状态一阶差分信息， $f(X_{t-\mathrm{1}})$ 可以参考式(13)~(17)近似得到。

结合第1~3节，将解耦极值应力的动态预测值进行相加，即可实现桥梁极值应力的动态预测。

4 算例分析

采用天津富民桥主梁危险截面（D截面）监测数据对所提算法进行验证。天津富民桥结构示意图、D截面位置示意图以及D截面传感器布置示意图分别如图1,2,3所示。图2中，A~F分别表示主梁截面的编号。天津富民桥总长340.3 m，桥宽40.0 m。桥梁动力响应十分复杂，主梁受温度荷载和车辆随机荷载影响较大，同时也存在自重恒载，即极值应力主要是由温度荷载和车辆荷载、自重恒载耦合产生的，满足研究动态模型的数据要求。

图1 天津富民桥

Fig.1 Tianjin Fumin Bridge

图2 主梁的监测截面

Fig.2 Monitoring sections for main girder

图3 D截面的传感器布置

Fig.3 Sensor layout of section D

在D截面横梁截面部分安装3个FBG应变传感器(FBG01015、FBG01005、FBG01012)，对应变进行动态监测，然后结合钢材的弹性模量，即可得到对应3个监测点的应力，进而可得极值应力数据。根据材料力学的基本原理可知,FBG01012监测的部位为最为危险的位置，因而本算例主要基于D截面FBG01012的监测极值应力数据对极值应力进行动态预测分析。

极值应力数据为每小时极值应力所形成的时间序列数据，总时长为2 905 h。对1~2 000 h监测极值应力数据（见图4）进行解耦，然后分别建立动态耦合线性模型，进而基于2 000~2 904 h的解耦极值应力数据动态预测2 001~2 905 h的极值应力，最后进行耦合。结合对应差分整合移动平均自回归模型，采用均方误差(MSE)进行准确度比较。

图4 1~2 000 h监测极值应力

Fig.4 Monitoring extreme stresses in 1~2 000 h

4.1　基于动态耦合线性模型的极值应力预测分析

对测点处1~2 000 h的极值应力数据利用简单移动平均法取5个时点进行平滑，得到解耦之后的车辆荷载效应和温度荷载效应时程曲线，如图5,6所示。监测点的自重恒载效应为43.08 MPa。

图5 车辆荷载效应时程曲线

Fig.5 Time history curve of vehicle load effects

图6 温度荷载效应时程曲线

Fig.6 Time history curve of temperature load effects

五点三次平滑处理后，即可近似得到初始状态信息，并将平滑后的数据进行一阶差分处理，可得状态一阶差分信息，如图7,8所示。

图7 车辆荷载效应监测数据、状态数据以及状态一阶差分值

Fig.7 Moritoring data, state data and the first-order differential data of vehicle load effects

图8 温度荷载效应监测数据、状态数据以及状态一阶差分值

Fig.8 Moritoring data, state data and the first-order differential data of temperature load effects

对初始状态信息及其一阶差分信息（即趋势信息）进行核密度估计^[19]，基于式(13)~(17)即可近似得到状态变量和状态一阶差分值的概率密度分布函数，如式(18)~(21)和图9,10,11,12所示。

图9 车辆荷载效应状态变量概率密度

Fig.9 Probability density of state variables of vehicle load effects

图10 车辆荷载效应状态一阶差分值的概率密度

Fig.10 Probability density of the first-order differencial data for the state of vehicle load effects

图11 温度荷载效应状态变量概率密度

Fig.11 Probability density of state variables of temperature load effects

图12 温度荷载效应状态一阶差分值的概率密度

Fig.12 Probability density of the first-order differencial data for the state of temperature load effects

（1）车辆荷载效应

$$X_{t-\mathrm{1}}|{\mathit{D}}_{t-\mathrm{1,1}}=\mathrm{0.42}N(X_{t-\mathrm{1}},-\mathrm{0.020}\mathrm{9},\mathrm{0.139}{\mathrm{7}}^{\mathrm{2}}/\mathrm{2})+\mathrm{0.58}N(X_{t-\mathrm{1}},\mathrm{0.013}\mathrm{0},\mathrm{0.221}{\mathrm{3}}^{\mathrm{2}}/\mathrm{2})$$

(18)

$${\beta }_{t-\mathrm{1}}|{\mathit{D}}_{t-\mathrm{1,2}}=\mathrm{0.46}N({\beta }_{t-\mathrm{1}},\mathrm{0.002}\mathrm{6},\mathrm{0.126}{\mathrm{8}}^{\mathrm{2}}/\mathrm{2})+\mathrm{0.54}N({\beta }_{t-\mathrm{1}},-\mathrm{0.002}\mathrm{8},\mathrm{0.187}{\mathrm{2}}^{\mathrm{2}}/\mathrm{2})$$

(19)

（2）温度荷载效应

$$X_{t-\mathrm{1}}|{\mathit{D}}_{t-\mathrm{1,1}}=\mathrm{0.64}N(X_{t-\mathrm{1}},-\mathrm{0.705}\mathrm{5},\mathrm{0.906}{\mathrm{2}}^{\mathrm{2}}/\mathrm{2})+\mathrm{0.36}N(X_{t-\mathrm{1}},\mathrm{1.087}\mathrm{0},\mathrm{1.046}{\mathrm{0}}^{\mathrm{2}}/\mathrm{2})$$

(20)

$${\beta }_{t-\mathrm{1}}|{\mathit{D}}_{t-\mathrm{1,2}}=\mathrm{0.38}N({\beta }_{t-\mathrm{1}},-\mathrm{0.008}\mathrm{9,0}.\mathrm{052}{\mathrm{4}}^{\mathrm{2}}/\mathrm{2})+\mathrm{0.62}N({\beta }_{t-\mathrm{1}},\mathrm{0.004}\mathrm{5},\mathrm{0.081}{\mathrm{9}}^{\mathrm{2}}/\mathrm{2})$$

(21)

结合式(1)~(3)，可建立动态耦合线性模型。

(1) 状态方程

$${\mathit{\theta }}_t=\mathit{G}{\mathit{\theta }}_{t-\mathrm{1}}+{\mathit{w}}_t,{\mathit{w}}_t~N(\mathrm{0},{\mathit{W}}_t)$$

(22)

$$\begin{array}{c}{\mathit{\theta }}_t=(X_t,{\beta }_t{)}^{\mathrm{T}},{\mathit{w}}_t=(w_{t,\mathrm{1}},w_{t,\mathrm{2}}{)}^{\mathrm{T}},\mathit{G}=\\ \left[\begin{array}{cc}\mathrm{1}& \mathrm{1}\\ \mathrm{0}& \mathrm{1}\end{array}\right],{\mathit{W}}_t=\left[\begin{array}{cc}{W}_{t,\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& W_{t,\mathrm{2}}\end{array}\right]\end{array}$$

(2) 监测方程

车辆荷载效应为

$$Z_t=X_t+v_t,v_t~N(\mathrm{0,0.128}\mathrm{7})$$

(23)

温度荷载效应为

$$Z_t=X_t+v_t,v_t~N(\mathrm{0,1.412}\mathrm{6})$$

(24)

初始状态信息如式(18)~(21)所示。进一步，结合式(1)、（2)和式(4)以及式(18)~(21)，动态耦合线性模型可转化为4个动态线性模型，并且4个动态线性模型的状态方程和监测方程均相同。

(1) 车辆荷载效应

第1~4个动态线性模型的初始信息分别为

$$X_{t-\mathrm{1,1}}|{\mathit{D}}_{t-\mathrm{1,1}}=N(X_{t-\mathrm{1}},\mathrm{0.020}\mathrm{9,0}.\mathrm{139}{\mathrm{7}}^{\mathrm{2}}/\mathrm{2})$$

$${\beta }_{t-\mathrm{1,1}}|{\mathit{D}}_{t-\mathrm{1,2}}=N({\beta }_{t-\mathrm{1}},\mathrm{0.002}\mathrm{6},\mathrm{0.126}{\mathrm{8}}^{\mathrm{2}}/\mathrm{2})$$

(25)

$$\begin{array}{c}{X}_{t-\mathrm{1,2}}|{\mathit{D}}_{t-\mathrm{1,1}}=N(X_{t-\mathrm{1}},-\mathrm{0.020}\mathrm{9},\mathrm{0.139}{\mathrm{7}}^{\mathrm{2}}/\mathrm{2})\\ {\beta }_{t-\mathrm{1,2}}|{\mathit{D}}_{t-\mathrm{1,2}}=N({\beta }_{t-\mathrm{1}},-\mathrm{0.002}\mathrm{8},\mathrm{0.187}{\mathrm{2}}^{\mathrm{2}}/\mathrm{2})\end{array}$$

(26)

$$\begin{array}{c}{X}_{t-\mathrm{1,3}}|{\mathit{D}}_{t-\mathrm{1,1}}=N(X_{t-\mathrm{1}},\mathrm{0.012}\mathrm{6},\mathrm{0.221}{\mathrm{3}}^{\mathrm{2}}/\mathrm{2})\\ {\beta }_{t-\mathrm{1,3}}|{\mathit{D}}_{t-\mathrm{1,2}}=N({\beta }_{t-\mathrm{1}},\mathrm{0.002}\mathrm{6},\mathrm{0.126}{\mathrm{8}}^{\mathrm{2}}/\mathrm{2})\end{array}$$

(27)

$$\begin{array}{c}{X}_{t-\mathrm{1,4}}|{\mathit{D}}_{t-\mathrm{1,1}}=N(X_{t-\mathrm{1}},\mathrm{0.012}\mathrm{6},\mathrm{0.221}{\mathrm{3}}^{\mathrm{2}}/\mathrm{2})\\ {\beta }_{t-\mathrm{1,4}}|{\mathit{D}}_{t-\mathrm{1,2}}=N({\beta }_{t-\mathrm{1}},-\mathrm{0.002}\mathrm{8},\mathrm{0.187}{\mathrm{2}}^{\mathrm{2}}/\mathrm{2})\end{array}$$

(28)

(2) 温度荷载效应

第1~4个动态线性模型的初始信息分别为

$$X_{t-\mathrm{1,1}}|{\mathit{D}}_{t-\mathrm{1,1}}=N(X_{t-\mathrm{1}},-\mathrm{0.705}\mathrm{5},\mathrm{0.906}{\mathrm{2}}^{\mathrm{2}}/\mathrm{2})$$ $${\beta }_{t-\mathrm{1,1}}|{\mathit{D}}_{t-\mathrm{1,2}}=N({\beta }_{t-\mathrm{1}},-\mathrm{0.008}\mathrm{9,0}.\mathrm{052}{\mathrm{4}}^{\mathrm{2}}/\mathrm{2})$$

(29)

$$\begin{array}{c}{X}_{t-\mathrm{1,2}}|{\mathit{D}}_{t-\mathrm{1,1}}=N(X_{t-\mathrm{1}},-\mathrm{0.705}\mathrm{5},\mathrm{0.906}{\mathrm{2}}^{\mathrm{2}}/\mathrm{2})\\ {\beta }_{t-\mathrm{1,2}}|{\mathit{D}}_{t-\mathrm{1,2}}=N({\beta }_{t-\mathrm{1}},\mathrm{0.004}\mathrm{5},\mathrm{0.081}{\mathrm{9}}^{\mathrm{2}}/\mathrm{2})\end{array}$$

(30)

$$\begin{array}{c}{X}_{t-\mathrm{1,3}}|{\mathit{D}}_{t-\mathrm{1,1}}=N(X_{t-\mathrm{1}},\mathrm{1.087}\mathrm{0},\mathrm{1.046}{\mathrm{0}}^{\mathrm{2}}/\mathrm{2})\\ {\beta }_{t-\mathrm{1,3}}|{\mathit{D}}_{t-\mathrm{1,2}}=N({\beta }_{t-\mathrm{1}},-\mathrm{0.008}\mathrm{9,0}.\mathrm{052}{\mathrm{4}}^{\mathrm{2}}/\mathrm{2})\end{array}$$

(31)

$$\begin{array}{c}{X}_{t-\mathrm{1,4}}|{\mathit{D}}_{t-\mathrm{1,1}}=N(X_{t-\mathrm{1}},\mathrm{1.087}\mathrm{0},\mathrm{1.046}{\mathrm{0}}^{\mathrm{2}}/\mathrm{2})\\ {\beta }_{t-\mathrm{1,4}}|{\mathit{D}}_{t-\mathrm{1,2}}=N({\beta }_{t-\mathrm{1}},\mathrm{0.004}\mathrm{5},\mathrm{0.081}{\mathrm{9}}^{\mathrm{2}}/\mathrm{2})\end{array}$$

(32)

对于第J(J=1,2,3,4)个动态线性模型，结合式(4)，对应的状态误差方差可由下式近似计算：

$${\mathit{W}}_{t,J}=-\left[\begin{array}{cc}\mathrm{1}& \mathrm{1}\\ \mathrm{0}& \mathrm{1}\end{array}\right]{\mathit{C}}_{t-\mathrm{1},J}{\left[\begin{array}{cc}\mathrm{1}& \mathrm{1}\\ \mathrm{0}& \mathrm{1}\end{array}\right]}^{\mathrm{T}}+\frac{{\mathit{C}}_{t-\mathrm{1},J}}{\mathrm{0.98}}$$

(33)

采用式(25)~(33)，参考式(1)~(11)，对

2 000~2 904 h的监测极值应力数据采用简单移动平均法逐步进行解耦，动态预测2 001~2 905 h的极值应力，得到如图13,14所示的一步预测结果。可以看出，预测值与实时监测值相接近，预测效果较好。将车辆荷载效应、温度荷载效应和自重恒载效应的预测结果进行相加，得到最终的耦合预测结果，如图15所示。可以看出，耦合后的预测值与实时监测值相接近，预测效果较好。

图13 车辆荷载效应监测极值应力与一步预测值

Fig.13 Monitoring and predicted extreme stress of vehicle load effects

图14 温度荷载效应监测极值应力与一步预测值

Fig.14 Monitoring and predicted extreme stress of temperature load effects

图15 监测极值应力与3种荷载效应加和的一步预测值

Fig.15 Monitoring and predicted extreme stress of three load effects

4.2　极值应力预测精度比较分析

为验证所建动态耦合线性模型的合理性和优越性，采用ARIMA（auto-regressive integrated moving average）模型，对预测精度进行比较分析。

利用ARIMA模型进行预测时，为了与动态耦合线性模型的预测结果进行比较，采用先差分再预测后还原的办法进行一步动态预测，两组数据预测结果如图16所示。可以看出，ARIMA模型预测值距离实时监测值较远，预测效果不如动态耦合线性模型。

图16 监测极值应力、ARIMA(10,1,15)与动态耦合线性

Fig.16 Monitoring extreme stresses and predicted extreme stresses based on ARIMA(10,1,15) and dynamic coupled linear model

注：

模型预测值

为进一步验证动态耦合线性模型的预测准确度，采用均方误差(MSE)进行比较。动态耦合线性模型的均方误差为0.117 2，ARIMA模型的均方误差为0.297 9。由预测结果可知，动态耦合线性模型具有较好的预测准确度。

5 结语

建立了解耦极值应力的动态耦合线性模型，并将其转化为有限个动态线性模型，然后采用贝叶斯方法对动态线性模型及其相关权重进行动态概率递推，进而实现动态耦合线性模型的动态概率递推。将动态预测的解耦极值应力进行相加，即可得到动态预测的桥梁极值应力，并利用实际桥梁监测数据进行验证分析。结果表明，动态耦合线性模型能够对由温度荷载和车辆荷载耦合产生的桥梁极值应力进行合理预测，预测值和监测值的变化趋势一致，大小近似相等。因此，动态耦合线性模型具有较好的预测准确度，可为桥梁预防性养护维修决策提供一定的理论基础。