基于结构随机场的桥梁移动荷载统计矩识别

饶勇平，张富博，雷 鹰; RAO Yong-ping; ZHANG Fu-bo; LEI Ying

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

摘要

由于结构、材料不确定等因素，桥梁结构往往具有随机性，因此，基于桥梁结构的响应识别得到的桥梁所受移动荷载也是不确定的，但相关研究还较少。为此，本文研究桥梁结构为空间随机场的情况下，有效识别移动荷载的统计矩。提出的识别方法基于桥梁结构随机场的Karhunen⁃Loeve（KL）展开、未知输入的卡尔曼滤波与改进的两点估计法的结合。桥梁结构参数为空间相关的随机场，通过KL展开将随机场转化为多随机变量的组合。在多随机变量下的不确定传播中，利用改进的两点估计法，将识别移动荷载统计矩问题转化为若干个确定性识别移动荷载识别的逆问题。采用未知输入的卡尔曼滤波进行确定性识别移动荷载，可有效估计识别移动荷载的统计矩。通过数值模拟算例验证了提出的估计方法的有效性。

关键词

移动荷载; Karhunen⁃Loeve展开; 改进两点估计法; 未知激励; 卡尔曼滤波

引言

移动荷载信息无论是在桥梁设计还是桥梁性能评估中都起着重要的作用^［

1⁃3］。然而，由于移动荷载的变化在时间和空间上同时发生，很难直接测量车辆和桥梁之间的相互作用力^{［参考文献 4

百度学术}4］。因此，发展根据桥梁结构的动态响应进行移动荷载的间接识别方法具有重要意义，近几十年也得到了广泛的研究。大多数研究将移动荷载和桥梁结构的相互作用视为确定性问题^{［参考文献 5⁃7}5⁃7］，例如Ding等^{［参考文献 8

百度学术}8］提出了一种基于平均加速度离散算法的离散力识别方法；Pan等^{［参考文献 9

百度学术}9］提出了一种改进的Tikhonov正则化方法来处理离散反问题；文献［10］提出了未知输入下的卡尔曼滤波方法（Kalman Filter with Unknown Input， KF⁃UI），可在部分观测结构响应的情况下，同时识别结构响应状态与未知输入。

在实际工程中，由于加工、制造及环境的影响，桥梁结构往往具有随机性，因此考虑结构随机性的识别方法更为适用。目前，随机性的考虑主要是将结构参数视作随机变量，如桥梁结构的质量线密度和阻尼系数为服从具体概率分布的随机变量。结构参数为随机变量的不确定性传播问题中，Liu等^［

11］研究了作用在薄壁圆柱壳结构上动态集中荷载边界的识别问题；方圣恩等^{［参考文献 12

百度学术}12］提出了一种结构参数识别的区间反演算法，通过构建两个方程来分别识别结构参数的中值和结构参数的半径；Wang等^{［参考文献 13

百度学术}13］提出了一种基于时域的分布式动态荷载识别方法。以上学者对结构参数为随机变量的研究中，大都采用基于一阶泰勒展开的区间分析方法，展开过程中舍弃了高阶项，展开项中用有限差分进行灵敏度的求解，都会导致不确定传播的计算结果存在误差。

实际工程中桥梁结构参数体现的随机性往往为具有空间相关性的随机场，随机场在其场域内的每个位置均为随机变量，即包含无限个随机变量，因此，基于桥梁结构参数为随机变量进行移动荷载的识别具有局限性。结构参数为随机场的不确定性问题中，Wu和Law^［

14⁃15］对移动荷载下的不确定性简支梁开展了研究，研究中将结构的弹性模量和质量线密度视为空间相关的随机场，通过KL展开将随机场转化为若干个随机变量的组合。孙燕伟等^{［参考文献 16

百度学术}16］研究了弹性模量和密度为随机场的悬臂梁模型，提出了在随机系统下集中荷载均值和方差的识别改进算法。上述研究方法中，需要大量的载荷识别，得到足够的识别荷载样本估计荷载的统计矩。因此，识别效率仍需提高。

在不确定性传播的方法中，点估计方法（Point Estimation Method， PEM）^［

17⁃20］是一种近似估计随机函数统计矩的方法。与上述识别移动荷载识统计矩的方法相比较，采用点估计方法，所需识别的样本少，估计效率高。在解决实际问题中，Hong等^{［参考文献 19

百度学术}19］提出的两点估计法（2 Point Estimation Method， 2PEM）和三点估计法（3 Point Estimation Method， 3PEM）得到广泛的使用^{［参考文献 21⁃22}21⁃22］。由于三点估计法（3PEM）在非对称概率密度的多随机变量情况下会出现虚数解的问题，两点估计（2PEM）又存在精度低的不足，车玉龙等^{［参考文献 23

百度学术}23］提出了增加统一概率的改进两点估计法（Improved 2 Point Estimation Method， I2PEM），在增加点估计方法精度的同时又不增加高阶矩的使用。

本文针对相关研究现状，尤其是随机结构情况下识别逆问题涉及不确定性传播和逆问题计算的双环过程，计算效率低的问题，研究桥梁结构参数为空间随机场情况下识别移动荷载统计矩的有效方法。提出的方法基于桥梁结构参数随机场的KL展开、KF⁃UI识别和I2PEM的结合。通过KL展开可以将桥梁结构参数的随机场转化为多个随机变量的组合。为了避免耗时的双环过程，不确定性传播采用I2PEM，可以将不确定性逆问题转化为若干个确定性逆问题。对移动荷载识别的逆问题，采用作者最近提出的基于数据融合的KF⁃UI方法^［

10］，可以在观测桥梁结构部分响应的情况下，识别未知移动荷载。最后，进行确定性识别移动荷载，可有效地估计移动荷载的统计矩。本文通过数值模拟识别算例，验证所提出方法的有效性。

1 提出的识别方法

1.1 基于KL展开考虑结构参数的随机场

实际工程中桥梁结构参数体现的随机性往往为具有空间相关性的随机场分布，随机场在其场域内的每个位置均为随机变量，即包含无限个随机变量。通过KL展开，可将随机场转化为若干个随机变量的组合。

假设 $d (x, θ)$ 是一个随机场，其中 $x$ 约束在场域 $D$ 内， $θ$ 是随机事件 $Ο$ 的空间。随机场 $d (x, θ)$ 在场域内的均值部分表示为 $\bar{d} (x)$ ，随机部分表示为 $\tilde{d} (x, θ)$ 。该随机场的协方差函数可以用谱分解表示为

C (x_{1}, x_{2}) = \sum_{n = 1}^{\infty} λ_{n} φ_{n} (x_{1}) φ_{n} (x_{2})

（1）

式中 $x_{1}$ 和 $x_{2}$ 分别表示场域内的两个位置； $λ_{n}$ 和 $φ_{n} (x)$ 分别表示协方差函数的特征值和特征函数，Ghanem和Spanos^［

24］证明了可以通过下面的积分方程来求解特征值和特征函数：

\int_{D} C (x_{1}, x_{2}) φ_{n} (x_{1}) d x_{1} = λ_{n} φ_{n} (x_{2})

（2）

协方差函数 $C (x_{1}, x_{2})$ 是对称的、正定的、有界的，其特征函数 $φ_{n} (x)$ 是正交的，并且所有的特征函数形成了用来表示协方差函数 $C (x_{1}, x_{2})$ 的完备集合，特征函数 $φ_{n} (x)$ 可以根据下式进行归一化

\int_{D} φ_{m} (x) φ_{n} (x) d x = δ_{m n}

（3）

式中 $δ_{m n}$ 是克罗内克函数，满足：

δ_{m n} = \{\begin{array}{l} 1, m = n \\ 0, m \neq n \end{array}

（4）

随机场 $d (x, θ)$ 可以写成^［

25］：

\begin{array}{l} d (x, θ) = \bar{d} (x) + \tilde{d} (x, θ) = \bar{d} (x) + \\ \sum_{n = 1}^{k_{d}} ξ_{n} (θ) \sqrt[]{λ_{n}} φ_{n} (x) \end{array}

（5）

式中 $k_{d}$ 是KL展开截断项的项数，由所截取的特征值之和占所有特征值总和的比例来确定^［

14⁃16］；

ξ_{n} (θ)

是不相关的随机变量，当

d (x, θ)

是一个高斯随机场时，

ξ_{n} (θ)

具有以下特性^{［参考文献 26⁃27}26⁃27］：

E (ξ_{n} (θ)) = 0, E (ξ_{m} (θ) ξ_{n} (θ)) = δ_{m n}

（6）

式中 $E (\cdot)$ 表示对 $\cdot$ 取期望。

1.2 考虑桥梁结构随机性的改进两点估计法

本文研究的桥梁模型为Bernoulli⁃Euler简支梁，移动荷载在结构上以匀速运动，Bernoulli⁃Euler简支梁的模型如图1所示。

图1 移动荷载下的Bernoulli-Euler简支梁

Fig.1 A simply supported Bernoulli-Euler beam under moving loads

在实际工程中，梁模型的结构参数往往是不确定的，在识别计算中，将会引起结构参数不确定到未知移动荷载不确定的传递。采用点估计法，可以较有效地估计由于结构参数不确定导致的未知移动荷载不确定性的统计矩。车玉龙等^［

23］针对传统三点估计法（3PEM）和两点估计法（2PEM）的不足，提出了改进两点估计法（I2PEM）。本文采用该方法，以结构参数的统计矩为输入信息，移动荷载的统计矩为输出信息。用于未知移动荷载识别的随机结构模型如下：

y = h (f^{u}, S)

（7）

式中 $f^{u}$ 为需要通过逆问题识别的作用在桥梁的未知移动荷载； $S$ 表示 $n$ 维随机变量， $S = (S_{1}, S_{2}, S_{3}, \dots, S_{n})$ ； $y$ 是观测的桥梁结构响应矢量， $h (\cdot)$ 表示观测方程。实际工况下，通过测量的结构响应的许多样本，统计得到响应的均值作为上述逆问题中观测的响应，进行未知移动荷载的识别。

从观测方程 $y = h (f^{u}, S)$ 出发，可以得到识别荷载的逆问题表达式：

f^{u} = \overset{⃖}{h} (y, S)

（8）

式中 $\overset{⃖}{h} (\cdot)$ 表示基于结构响应 $y$ 和随机变量 $S$ 对未知移动荷载 $f^{u}$ 进行识别的逆问题函数。

令 $μ_{k}$ 和 $σ_{k}$ 分别表示随机变量 $S_{k}$ 的均值和标准差， $M_{k, i}^{'} (S_{k})$ 表示随机变量 $S_{k}$ 的第 $i$ 阶中心矩， $i = 1,2, 3, \dots$ ； $k = 1,2, 3, \dots, n$ 。

M_{k, i}^{'} (S_{k}) = \int_{- \infty}^{+ \infty} {(S_{k} - μ_{k})}^{i} f_{S_{k}} (S_{k}) d S_{k}

（9）

式中 $f_{S_{k}} (S_{k})$ 为随机变量 $S_{k}$ 的概率密度函数。

令 $λ_{k, i}$ 为 $M_{k, i}^{'} (S_{k})$ 和 $σ_{k}^{i}$ 的比，即：

λ_{k, i} = M_{k, i}^{'} (S_{k}) / σ_{k}^{i}

（10）

式中 $λ_{k, 1} = 0$ ， $λ_{k, 2} = 1$ 。

值得注意的是，若随机变量 $S_{k}$ 服从高斯分布，第 $i$ 阶中心矩 $M_{k, i}^{'} (S_{k})$ 满足

M_{k, i}^{'} (S_{k}) = \{\begin{array}{l} 0, i 为奇 数 \\ σ_{k}^{i} (i - 1) (i - 3) \dots 3 \cdot 1 \end{array} \begin{matrix} , i 为偶 数 \end{matrix}

（11）

结合式（10）和（11），可以得出服从高斯分布的随机变量 $S_{k}$ 的偏斜系数 $λ_{k, 3} = 0$ 。

随机变量 $S_{k}$ 在点估计法中第 $i$ 个选点 $s_{k, i}$ 可以表述为： $s_{k, i} = μ_{k} + ξ_{k, i} σ_{k}$ ， $i = 1,2$ ； $k = 1,2, 3, \dots, n$ ，其中 $ξ_{k, i}$ 为随机变量 $S_{k}$ 在点估计法中的第 $i$ 个选点的系数。令 $p_{k, 1}$ 和 $p_{k, 2}$ 分别对应于 $s_{k, 1}$ 和 $s_{k, 2}$ 的两个待定的集中概率。此时，每个随机变量 $S_{k}$ 都可以建立如下四个方程：

\begin{array}{l} p_{k, 1} ξ_{k, 1}^{i} + p_{k, 2} ξ_{k, 2}^{i} = λ_{k, i}, \\ i = 0,1, 2,3; k = 1,2, 3, \dots, n \end{array}

（12）

同时，两个集中概率 $p_{k, 1}$ 和 $p_{k, 2}$ 满足：

\{\begin{array}{l} \sum_{k = 1}^{n} (p_{k, 1} + p_{k, 2}) = 1 \\ p_{k, 1} + p_{k, 2} = 1 / n, k = 1,2, 3, \dots, n \end{array}

（13）

联立式（12）和（13）可求解得到：

\begin{matrix} ξ_{k, i} = \frac{λ_{k, 3}}{2} + {(- 1)}^{3 - i} \sqrt[]{n + \frac{λ_{k, 3}^{2}}{4}}, \\ p_{k, i} = \frac{\frac{1}{n} {(- 1)}^{i} ξ_{k, 3 - i}}{2 \sqrt[]{n + λ_{k, 3}^{2} / 4}}, \\ i = 1,2; k = 1,2, 3, \dots, n \end{matrix}

（14）

多随机变量下 $f^{u}$ 的第 $j$ 阶矩即为：

E [{(f^{u})}^{j}] = \sum_{k = 1}^{n} \sum_{i = 1}^{2} p_{k, i} {(\overset{⃖}{h} (μ_{1}, μ_{2}, μ_{2}, \dots, s_{k, i}, \dots, μ_{n}))}^{j}

（15）

上述2PEM^［

19］又称

2 n

集合方案，2PEM需要用到随机变量的前三阶中心矩。理论上，点估计方法中用到的点数越多，待估值的准确度越高。但是，随着估计点数的增多，需要用到随机变量的更高阶中心矩，高阶中心矩会使随机变量标准位置和权重系数的求解出现非实值，标准位置和权重系数没有解析解。对此，车玉龙等^{［参考文献 23

百度学术}23］提出了改进的两点估计方法（I2PEM），增加估计点的点数以提高估计准确度，但不需要用到随机变量的高阶中心矩。I2PEM基于上面描述的2PEM，对每个随机变量

S_{k}

新增三个点，其中两个为同一概率的估计点，另一个为均值估计点。

新增的两个同一概率的估计点为：

\begin{array}{l} ξ_{k, i}^{'} = {(- 1)}^{3 - i} \sqrt[]{2 n} \\ p_{k, i}^{'} = \frac{1}{2 n}, i = 1,2; k = 1,2, 3, \dots, n \end{array}

（16）

式中随机变量 $S_{k}$ 新增的第 $i$ 个位置 $s_{k, i}^{'}$ 为 $s_{k, i}^{'} = μ_{k} + ξ_{k, i}^{'} σ_{k}$ ， $ξ_{k, i}^{'}$ 为随机变量 $S_{k}$ 新增的第 $i$ 个标准位置； $p_{k, 1}^{'}$ 和 $p_{k, 2}^{'}$ 分别对应于 $s_{k, 1}$ 和 $s_{k, 2}$ 的两个新增权重系数。

每个随机变量 $S_{k}$ 新增一个均值估计点，即重复进行 $n$ 次同样的识别运算（ $f^{u} = \overset{⃖}{h} (μ_{1}, μ_{2}, μ_{2}, \dots, μ_{n})$ ），故对待估值的计算只需增加一个权重系数 $p_{0}^{} = 1$ 的均值点。结合式（15）和（16），对多随机变量下 $f^{u}$ 第 $j$ 阶矩进行估计：

\begin{array}{l} E [{(f^{u})}^{j}] = \\ \frac{1}{3} [\sum_{k = 1}^{n} \sum_{i = 1}^{2} p_{k, i} {(\overset{⃖}{h} (μ_{1}, μ_{2}, μ_{2}, \dots, s_{k, i}, \dots, μ_{n}))}^{j} + \\ \sum_{k = 1}^{n} \sum_{i = 1}^{2} p_{k, i}^{'} {(\overset{⃖}{h} (μ_{1}, μ_{2}, μ_{2}, \dots, s_{k, i}^{'}, \dots, μ_{n}))}^{j} + \\ p_{0}^{} {(\overset{⃖}{h} (μ_{1}, μ_{2}, μ_{2}, \dots, μ_{n}))}^{j}] \end{array}

（17）

上式即为I2PEM的公式，又称为 $4 n + 1$ 集合方案。原有的两点用到了随机变量的前三阶中心矩，新增三点用到了随机变量的前两阶中心矩，可以看出 $4 n + 1$ 集合方案无需随机变量的高阶中心矩，即可增加估计点数。

1.3 基于桥梁部分响应的移动荷载识别

基于桥梁部分响应的识别方法中，采用文献［

10］提出的KF⁃UI，结构运动方程如下：

M \ddot{x} + C \dot{x} + K x = η^{u} f^{u}

（18）

式中 $x$ ， $\dot{x}$ 和 $\ddot{x}$ 分别代表结构的位移、速度和加速度； $K$ 和 $M$ 分别为简支梁结构的刚度矩阵和质量矩阵； $C$ 为结构的阻尼矩阵，研究中采用瑞利阻尼； $f^{u} = {[\begin{matrix} f_{1}^{u} & f_{2}^{u} \end{matrix}]}^{T}$ 为未知的两个移动荷载， $η^{u}$ 是对应于未知移动 $f^{u}$ 的定位矩阵。

通过设置状态向量 $Z = {[\begin{matrix} x^{T} & {\dot{x}}^{T} \end{matrix}]}^{T}$ ，系统的状态方程可以表示为：

\dot{Z} = G Z + E f^{u}

（19）

式中 $G = [\begin{matrix} 0 & I \\ - M^{- 1} K & - M^{- 1} C \end{matrix}]$ ， $E = [\begin{matrix} 0 \\ M^{- 1} η^{u} \end{matrix}]$ 。

考虑系统模型误差 $w_{k}$ ，其均值为 $0$ ，方差为 $Q$ ，离散式（19）可以得到下式：

Z_{k + 1} = A_{k} Z_{k} + B_{k} f_{k}^{u} + w_{k}

（20）

式中 $A_{k} = e^{G Δ t}$ ， $B_{k} = (A_{k} - I) {(G Δ t)}^{- 1} (E Δ t)$ 。

系统的观测方程可以表示为：

\begin{array}{l} y_{k + 1} = [\begin{matrix} L^{s} Θ x_{k + 1} \\ L^{a} {\ddot{x}}_{k + 1} \end{matrix}] + v_{k + 1} = H_{k + 1} Z_{k + 1} + \\ D_{k + 1} f_{k + 1}^{u} + v_{k + 1} \end{array}

（21）

式中 $y_{k + 1}$ 代表观测量，即观测结构应变响应和加速度响应； $L^{s}$ 和 $L^{a}$ 分别表示应变和加速度的位置观测矩阵； $v_{k + 1}$ 是观测误差，其均值为 $0$ ，方差为 $R$ ； $Θ$ 表示结构位移和应变之间的转换矩阵。

结构状态的时间预测和观测更新如下所示：

{\tilde{Z}}_{k + 1 | k}^{} = A_{k} {\hat{Z}}_{k | k} + B_{k} {\hat{f}}_{k | k}^{u}

（22）

{\hat{Z}}_{k + 1 | k + 1}^{} = {\tilde{Z}}_{k + 1 | k}^{} + K_{k + 1} (y_{k + 1} - H_{k + 1} {\tilde{Z}}_{k + 1 | k}^{} - D_{k + 1} f_{k + 1}^{u})

（23）

估计误差的协方差矩阵为：

{\tilde{P}}^{Z}_{k + 1 | k} = [\begin{matrix} A_{k} & B_{k} \end{matrix}] [\begin{matrix} {\hat{P}}^{Z}_{k | k} & {\hat{P}}^{Z f}_{k | k} \\ {\hat{P}}^{f Z}_{k | k} & {\hat{P}}^{f}_{k | k} \end{matrix}] [\begin{matrix} A_{k}^{T} \\ B_{k}^{T} \end{matrix}] + Q

（24）

\begin{matrix} {\hat{P}}^{Z}_{k + 1 | k + 1} = {\tilde{P}}^{Z}_{k + 1 | k} - K_{k + 1} H_{k + 1} {\tilde{P}}^{Z}_{k + 1 | k} - \\ {\tilde{P}}^{Z}_{k + 1 | k} H_{k + 1}^{T} K_{k + 1}^{T} + K_{k + 1} \cdot \\ (H_{k + 1} {\tilde{P}}^{Z}_{k + 1 | k} H_{k + 1}^{T} + R) K_{k + 1}^{T} \end{matrix}

（25）

通过最小化 ${\hat{P}}_{k + 1 | k + 1}^{Z}$ 的迹来估计卡尔曼增益矩阵 $K_{k + 1}$ ：

K_{k + 1} = {\tilde{P}}^{Z}_{k + 1 | k} H_{k + 1}^{T} {(H_{k + 1} {\tilde{P}}^{Z}_{k + 1 | k} H_{k + 1}^{T} + R)}^{- 1}

（26）

将其代入式（25）后得出以下结果：

\begin{matrix} {\hat{P}}^{Z}_{k + 1 | k + 1} = (I + K_{k + 1} D_{k + 1} S_{k + 1} D_{k + 1}^{T} R_{k + 1}^{- 1} H_{k + 1}) \cdot \\ (I - K_{k + 1} H_{k + 1}) {\tilde{P}}^{Z f}_{k + 1 | k} \end{matrix}

（27）

S_{k + 1} = {[D_{k + 1}^{T} R_{k + 1}^{- 1} (I - H_{k + 1} K_{k + 1}) D_{k + 1}]}^{- 1}

（28）

根据式（23），估计输出为：

{\hat{y}}_{k + 1 | k + 1} = H_{k + 1} {\hat{Z}}_{k + 1 | k + 1}^{} + D_{k + 1} f_{k + 1}^{u}

（29）

然后，将观测误差定义为：

\begin{array}{l} Δ_{k + 1} = y_{k + 1} - {\hat{y}}_{k + 1 | k + 1} = \\ (I - H_{k + 1} K_{k + 1}) (y_{k + 1} - H_{k + 1} {\overset{⌢}{Z}}_{k + 1 | k} - D_{k + 1} f_{k + 1}^{u}) \end{array}

（30）

其中， ${\overset{⌢}{Z}}_{k + 1 | k} = A_{k} {\hat{Z}}_{k | k} + B_{k} {\hat{f}}_{k | k}^{u}$ 。

基于最小二乘法估计未知移动荷载 $f_{k + 1}^{u}$ ：

{\hat{f}}_{k + 1 | k + 1}^{u} = S_{k + 1} D_{k + 1}^{T} R_{k + 1}^{- 1} (I - H_{k + 1} K_{k + 1}) \cdot (y_{k + 1} - H_{k + 1} {\overset{⌢}{Z}}_{k + 1 | k})

（31）

相应地，未知力的估计误差协方差由下式给出：

{\hat{P}}^{f}_{k + 1 | k + 1} = S_{k + 1}

（32）

交叉项的协方差为：

\begin{array}{l} {\hat{P}}^{Z f}_{k + 1 | k + 1} = - K_{k + 1} D_{k + 1} S_{k + 1}, \\ {\hat{P}}^{f Z}_{k + 1 | k + 1} = {({\hat{P}}^{Z f}_{k + 1 | k + 1})}^{T} \end{array}

（33）

上述KF⁃UI利用桥梁部分观测的应变响应和加速度响应，可实现对未知移动荷载的识别。

结合上一节的I2PEM，随机变量 $S_{k}$ 的第 $i$ 个位置 $s_{k, i}$ 为 $s_{k, i} = μ_{k} + ξ_{k, i} σ_{k}$ ， $i = 1,2$ ； $k = 1,2, 3, \dots, n$ ，每个位置 $s_{k, i}$ 的移动荷载识别如下：

f_{k, i}^{u} = \overset{⃖}{h} (μ_{1}, μ_{2}, μ_{2}, \dots, s_{k, i}, \dots, μ_{n})

（34）

同理，随机变量 $S_{k}$ 新增的第 $i$ 个位置 $s_{k, i}^{'}$ 为 $s_{k, i}^{'} = μ_{k} + ξ_{k, i}^{'} σ_{k}$ ， $i = 1,2$ ； $k = 1,2, 3, \dots, n$ ，每个位置 $s_{k, i}^{'}$ 的移动荷载识别如下：

f_{k, i}^{u'} = \overset{⃖}{h} (μ_{1}, μ_{2}, μ_{2}, \dots, s_{k, i}^{'}, \dots, μ_{n})

（35）

移动荷载在随机变量均值处的识别如下：

f_{μ}^{u} = \overset{⃖}{h} (μ_{1}, μ_{2}, μ_{2}, \dots, μ_{n})

（36）

多随机变量下移动荷载 $f^{u}$ 第 $j$ 阶矩为：

\begin{array}{l} E [{(f^{u})}^{j}] = \frac{1}{3} [\sum_{k = 1}^{n} \sum_{i = 1}^{2} p_{k, i} {(f_{k, i}^{u})}^{j} + \\ \sum_{k = 1}^{n} \sum_{i = 1}^{2} p_{k, i}^{'} {(f_{k, i}^{u'})}^{j} + p_{0}^{} {(f_{μ}^{u})}^{j}] \end{array}

（37）

1.4 桥梁结构参数为随机场的移动荷载统计矩的识别

综合1.1节至1.3节描述，桥梁结构参数为随机场的移动荷载统计矩的识别步骤如下：

（1 ）对服从随机场的桥梁结构参数进行KL展开，将随机场转化为若干随机变量的组合；

（2）每个随机变量 $S_{k}$ 的 $λ_{k}$ ， $k = 1,2, 3, \dots, n$ ；利用式（14）计算随机变量的两个标准位置、两个权重系数，然后计算其位置 $s_{k, i}^{} = μ_{k} + ξ_{k, i}^{} σ_{k}$ ， $i = 1,2$ ；利用式（16）计算随机变量的两个新增标准位置、两个新增权重系数，然后计算其新增位置 $s_{k, i}^{'} = μ_{k} + ξ_{k, i}^{'} σ_{k}$ ， $i = 1,2$ ；

（3）基于桥梁结构的部分响应，采用KF⁃UI对两个位置，两个新增位置，随机变量均值位置进行移动荷载识别计算，根据式（37）估计识别的移动荷载的统计矩；

（4）采用蒙特卡洛方法，在桥梁结构参数随机的情况下利用结构部分响应进行识别，模拟得到真实值。

2 数值模拟验证

本节采用的Bernoulli⁃Euler简支梁如图2所示，划分为10个等长梁单元，节点标号从左到右依次编号。简支梁总长度 $L = 15 m$ ，每个梁单元的长度为 $L_{e} = 1.5 m$ 。梁系统中共有20个自由度，包含9个垂直自由度和11个旋转自由度。结构截面惯性矩 $I_{e} = 2.304 m^{4}$ ，梁横截面高度 $h = 2.4 m$ 。算例中采用瑞利阻尼，前两阶阻尼比取值为0.02。简支梁的弹性模量 $E$ 和质量线密度 $m$ 为随机场，均服从高斯分布。弹性模量 $E$ 的均值为 $5 \times 10^{10} P a$ ，变异系数为5%。质量线密度 $m$ 的均值为 $1.2 \times 10^{5} k g / m$ ，变异系数也为5%。弹性模量 $E$ 和质量线密度 $m$ 具有空间相关性，两者的协方差函数均假定为指数型模型^［

16］：

C (x_{1}, x_{2}) = σ^{2} e x p (- |x_{1} - x_{2}| / a)

（38）

式中 $σ$ 为弹性模量 $E$ 和质量线密度的标准差； $a$ 为空间相关长度，本算例中取为结构单元的长度。

图2 移动荷载下的Bernoulli-Euler简支梁

Fig. 2 A simply supported Bernoulli-Euler beam under moving loads

以形成简支梁模型的刚度矩阵为例，弹性模量 $E$ 为随机场时，对其进行KL展开：

\begin{array}{l} E (x, θ) = \bar{E} (x) + \tilde{E} (x, θ) = \bar{E} (x) + \\ \sum_{i 1 = 1}^{k_{E}} ξ_{i 1} (θ) \sqrt[]{λ_{i 1}} φ_{i 1} (x) \end{array}

（39）

式中 $\bar{E} (x)$ 和 $\tilde{E} (x, θ)$ 分别表示随机场 $E$ 的均值部分和随机部分； $k_{E}$ 是KL展开截断项，算例中确定为20； $λ_{i 1}$ 和 $φ_{i 1} (x)$ 分别表示协方差函数的特征值和特征函数。

梁单元刚度矩阵如下：

{\bar{K}}_{e} = \int_{L_{e}} B^{T} \bar{E} (x) I B d x

（40）

{\tilde{K}}_{e} = \int_{L_{e}} B^{T} \tilde{E} (x, θ) I B d x

（41）

式中 ${\bar{K}}_{e}$ 和 ${\tilde{K}}_{e}$ 分别表示单元刚度矩阵的均值部分和随机部分； $B$ 为对应的形函数。

将式（39）代入式（41）可得出：

\begin{array}{l} {\tilde{K}}_{e} = \int_{L_{e}} B^{T} [\sum_{i 1 = 1}^{k_{E}} ξ_{i 1} (θ) \sqrt[]{λ_{i 1}} φ_{i 1} (x)] I B d x = \\ \sum_{i 1 = 1}^{k_{E}} ξ_{i 1} (θ) \int_{L_{e}} B^{T} \sqrt[]{λ_{i 1}} φ_{i 1} (x) I B d x = \sum_{i 1 = 1}^{k_{E}} ξ_{i 1} (θ) K_{i 1}^{} \end{array}

（42）

则单元刚度矩阵为：

K_{e} = {\bar{K}}_{e} + {\tilde{K}}_{e} = {\bar{K}}_{e} + \sum_{i 1 = 1}^{k_{E}} ξ_{i 1} (θ) K_{i 1}

（43）

同理可得到单元质量矩阵如下

M_{e} = {\bar{M}}_{e} + {\tilde{M}}_{e} = {\bar{M}}_{e} + \sum_{i 2 = 1}^{k_{m}} ξ_{i 2} (θ) M_{i 2}

（44）

分别对单元刚度矩阵和单元质量矩阵进行组装，可得到整体刚度矩阵和整体质量矩阵。阻尼矩阵采用瑞利阻尼。

移动荷载由车辆重量和简谐荷载叠加而成，其时程表达式为：

\begin{array}{l} f_{}^{u} (t) = [\begin{array}{l} f_{1}^{u} (t) \\ f_{2}^{u} (t) \end{array}] = \\ [\begin{array}{l} 2000 s i n (10 π t) + 1000 s i n (40 π t) + 20000 \\ 2000 s i n (10 π t) - 1000 s i n (40 π t) + 20000 \end{array}] \end{array}

（45）

式中 $f_{1}^{u} (t)$ 表示前轮移动荷载； $f_{2}^{u} (t)$ 表示后轮移动荷载。两个移动荷载的速度为 $v = 15 m / s$ ，前后轮距 $l_{a} = 3 m$ 。

未知移动荷载统计矩识别值与准确值的相对误差表示为：

R E \underset{̲}{} S^{f} = \frac{{‖S_{i d e n t i f i e d}^{f} - S_{e x a c t}^{f}‖}_{2}}{{‖S_{e x a c t}^{f}‖}_{2}} \times 100 %

（46）

式中 ${‖\cdot‖}_{2}$ 表示2⁃范数； $S^{f}$ 为移动荷载的统计矩；下标“identified”和“exact”分别表示识别值和真实值。采用蒙特卡洛方法，在结构参数随机的情况下利用结构部分响应模拟识别10000次得到真实值。

采样持续时间为0.8 s，传感器采样频率为1000 Hz。移动荷载后轮位于简支梁左端时采样开始，移动荷载前轮行驶至简支梁末端时采样结束。加速度传感器部署在结构第2，4，5，6，7，9，10节点，应变传感器部署在结构第2，5，6，10节点。且获取的传感器信号都受到2% RMS噪声。

在多重随机场（弹性模量、质量线密度）下，利用KF⁃UI和I2PEM，可得到未知移动荷载的前三阶矩如图3~8所示。

图3 前轮移动荷载一阶矩识别图

Fig.3 The 1st moment of identified front moving load

图4 前轮移动荷载二阶矩识别图

Fig.4 The 2nd moment of identified front moving load

图5 前轮移动荷载三阶矩识别图

Fig.5 The 3rd moment of identified front moving load

图6 后轮移动荷载一阶矩识别图

Fig.6 The 1st moment of identifed rear moving load

图7 后轮移动荷载二阶矩识别图

Fig.7 The 2nd moment of identified rear moving load

图8 后轮移动荷载三阶矩识别图

Fig.8 The 3rd moment of identified rear moving load

图3~5为前轮移动荷载统计矩（一阶矩、二阶矩和三阶矩）的识别情况，图6~8为后轮移动荷载统计矩的识别图。图中表明，利用I2PEM和KF⁃UI得到的识别值与准确值非常接近，说明了即便是存在双随机场（弹性模量和质量线密度），I2PEM结合KF⁃UI识别算法仍然具备有效性。

表1研究了结构双随机场下，简支梁在I2PEM方法下前后轮移动荷载统计矩识别的相对误差。横向研究对象包含识别的一阶矩、二阶矩和三阶矩，竖向研究对象分别为前后轮移动荷载的识别。移动荷载统计矩识别的相对误差随着阶数提高而增大，但是前轮移动荷载识别的最大相对误差不超过1.5%，后轮移动荷载识别的最大相对误差不超过0.9%，说明所用识别算法的有效性。

表1 移动荷载统计矩的相对误差/%

Tab. 1 The relative error of statistical moving loads/%

	一阶矩	二阶矩	三阶矩
前轮移动荷载	0.415	0.903	1.471
后轮移动荷载	0.295	0.567	0.825

3 结论

本文考虑了桥梁结构随机性识别的移动荷载统计矩有效估计。随机结构的逆问题研究通常涉及带有不确定性传播和逆问题计算的双环过程，提出的方法基于桥梁结构参数随机场的KL展开与KF⁃UI识别和I2PEM的结合，首先通过KL展开将服从空间分布的随机场（弹性模量和质量线密度）转化为若干个随机变量的组合。然后基于I2PEM对每个不确定随机变量的近似位置进行选择，将不确定性逆问题转化为多个确定性逆问题。最后基于桥梁结构的部分响应对每个随机变量选择的位置进行确定性的移动荷载识别，估计移动荷载的统计矩。利用移动荷载在梁桥模型下的数值模拟，验证了桥梁结构参数为随机场情况下识别的移动荷载统计矩的估计方法的有效性。

本文考虑桥梁结构的弹性模量和质量线密度为服从高斯分布的随机场，实际工程下桥梁结构参数服从其他类型分布时，如非高斯分布，可用混沌多项式展开，类似地进行移动荷载识别的统计矩的有效估计。本文中采用点估计方法进行桥梁移动荷载统计矩的估计，由于点估计方法在随机变量均值处采用泰勒展开，当随机变量的变异系数较大时，统计矩估计精度会受到影响。这需要进一步研究加以改进。另外，本文假定车辆重量为确定，也需研究拓展到考虑车辆荷载为随机变量的实际情况。

参考文献

Yi T H， Li H N， Gu M. Wavelet based multi-step filtering method for bridge health monitoring using GPS and accelerometer［J］. Smart Structures & Systems， 2013， 11（4）： 331-348. [百度学术]

Chen H P. Residual flexural capacity and performance assessment of corroded reinforced concrete beams［J］. Journal of Structural Engineering， 2018， 144（12）： 04018213. [百度学术]

Chen B， Ye Z N， Chen Z， et al. Bridge vehicle load model on different grades of roads in China based on Weigh-in-Motion （WIM） data［J］. Measurement， 2018， 122： 670-678. [百度学术]

Liu H， Ling Y. Sparse regularization for traffic load monitoring using bridge response measurements［J］. Measurement， 2019， 131： 173-182. [百度学术]

Law S S， Chan T， Zeng Q H. Moving force identification： a time domain method ［J］. Journal of Sound & Vibration， 1997， 201（1）： 1-22. [百度学术]

Yu Y， Zhao X， Shi Y， et al. Design of a real-time overload monitoring system for bridges and roads based on structural response［J］. Measurement， 2013， 46（1）： 345-352. [百度学术]

Pan C D， Yu L， Liu H L， et al. Moving force identification based on redundant concatenated dictionary and weighted l1-norm regularization［J］. Mechanical Systems and Signal Processing， 2018， 98： 32-49. [百度学术]

Ding Y， Law S S， Wu B， et al. Average acceleration discrete algorithm for force identification in state space［J］. Engineering Structures， 2013， 56： 1880-1892. [百度学术]

Pan C D， Yu L， Liu H L. Identification of moving vehicle forces on bridge structures via moving average Tikhonov regularization［J］. Smart Materials and Structures， 2017， 26（8）： 085041. [百度学术]

Liu L J， Zhu J J，Su Y， et al. Improved Kalman filter with unknown inputs based on data fusion of partial acceleration and displacement measurements［J］. Smart Structures & Systems， 2016， 17（6）： 903-915. [百度学术]

Liu J， Han X， Jiang C， et al. Dynamic load identification for uncertain structures based on interval analysis and regularization method［J］. International Journal of Computational Methods， 2011， 8（4）： 667-683. [百度学术]

方圣恩，张秋虎，林友勤.基于泰勒级数展开的区间反演方法［J］.计算力学学报， 2015， 32（6）： 796-802. [百度学术]

Fang Shengen， Zhang Qiuhu， Lin Youqin. An interval inverse solution method based on Taylor series expansion［J］. Chinese Journal of Computational Mechanics， 2015， 32（6）： 796-802. [百度学术]

Wang L， Liu Y， Liu Y. An inverse method for distributed dynamic load identification of structures with interval uncertainties［J］. Advances in Engineering Software， 2019， 131：77-89. [百度学术]

Wu S Q， Law S S. Moving force identification based on stochastic finite element model［J］. Engineering Structures， 2010， 32（4）： 1016-1027. [百度学术]

Wu S Q， Law S S. Dynamic analysis of bridge-vehicle system with uncertainties based on the finite element model［J］. Probabilistic Engineering Mechanics， 2010， 25（4）： 425-432. [百度学术]

孙燕伟，吴邵庆，李彦斌，等.不确定性动力学系统上随机载荷识别改进算法［J］.振动工程学报， 2019， 32（2）： 206-214. [百度学术]

Sun Yanwei， Wu Shaoqing， Li Yanbin， et al. An improved algorithm for stochastic load identification for random system［J］. Journal of Vibration Engineering， 2019， 32（2）： 206-214. [百度学术]

Rosenblueth E. Point estimation for probability moments［J］. Proceedings of the National Academy of Sciences of the United States of America， 1975， 72： 3812-3814. [百度学术]

Rosenblueth E. Two-point estimates in probabilities［J］. Applied Mathematical Modelling， 1981， 5（5）： 329-335. [百度学术]

Hong H P. An efficient point estimate method for probabilistic analysis［J］. Reliability Engineering & System Safety， 1998， 59（3）： 261-267. [百度学术]

刘祥，蒋丽忠，向平，等. 基于改进点估计法的随机车桥竖向振动分析［J］. 振动与冲击， 2020， 39（6）： 15-21. [百度学术]

Liu Xiang， Jiang Lizhong， Xiang Ping， et al. Analysis of train-bridge vertical random vibration based on a new point［J］. Journal of Vibration and Shock， 2020， 39（6）：15-21. [百度学术]

Su C L. Probabilistic load-flow computation using point estimate method［J］. IEEE Transactions on Power Systems， 2005， 20（4）： 1843-1851. [百度学术]

Morales J M， Pérez-Ruiz J， et al. Point estimate schemes to solve the probabilistic power flow［J］. IEEE Transactions on Power Systems， 2005， 22（4）： 1594-1601. [百度学术]

车玉龙，吕晓琴，王晓茹，等. 含非正态分布概率潮流计算的改进型两点估计法［J］. 电力自动化设备， 2019， 39（12）： 128-133. [百度学术]

Che Yulong， Lü Xiaoqin， Wang Xiaoru， et al. Improved two point estimation method for probabilistic power flow calculation with non-normal distribution［J］. Electric Power Automation Equipment， 2019， 39（12）： 128-133. [百度学术]

Ghanem R G， Spanos P D. Stochastic Finite Elements： A Spectral Approach［M］. Springer-Verlag， 1991. [百度学术]

Wu S Q， Law S S. Statistical moving load identification including uncertainty［J］. Probabilistic Engineering Mechanics， 2012， 29： 70-78. [百度学术]

Wu S Q， Law S S. Vehicle axle load identification on bridge deck with irregular road surface profile［J］. Engineering Structures， 2011， 33（2）： 591-601. [百度学术]

Wu S Q， Law S S. A reduced polynomial chaos expansion model for stochastic analysis of a moving load on beam system with non-Gaussian parameters［J］. Journal of Vibroengineering，2015， 17（3）： 1560-1577. [百度学术]

基于结构随机场的桥梁移动荷载统计矩识别 PDF

摘要

关键词

引 言