引言

由于载荷、环境以及材料内部因素的作用,结构的性能一般随时间而逐渐退化^[1]. 为了便于评估结构服役期间的状态,常采用随机变量模型来描述结构性能的退化规律^[2,3]. 在工程结构分析及设计中,需考虑随机不确定和认知不确定共同发生影响的情形^[4].

构建初始物理模型可以利用同类型结构的数据信息,通过识别并提取有效信息得到模型参数的先验分布. 然而,仅利用一些先验数据难以得到准确的物理模型^[5,6,7]. 利用结构健康检测系统和无损检测技术可获得服役结构的一些检测数据,利用贝叶斯方法整合这些数据对先验分布进行更新,得到模型参数的后验分布,可有效提高物理模型的准确性.

在确定模型先验分布时,Guan等^[8]采用最大似然估计确定每条数据集的模型参数,然后通过数学统计确定模型参数的均值和协方差矩阵. He等^[9]采用证据理论(Dempster-shafer theory,DST)分配数据集的可信度来确定模型参数. 后验分布的计算会涉及到高维积分,解决这个问题通常采用解析近似和数值模拟方法^[10]. Luke等^[11]提出了拉普拉斯方法(Laplace)对目标分布做解析近似,以方便多维积分的计算. Guan等^[8]基于Laplace和逆一阶可靠性方法更新可靠性和系统响应. 然而Laplace解析近似方法虽缩短了计算时间,但仍然需要更稳健和更精确的数值模拟作为基准^[12]. Cappé等^[13]利用马尔科夫链蒙特卡洛法(MCMC)计算出积分结果. 粒子滤波算法(particle filter, PF)利用动态空间模型可以提供后验计算结果,尤其是可以使得似然函数具有无偏估计性^[14]. Fernández等^[15]首次将MCMC框架引入PF算法中,并应用于宏观经济模型的后验分布估计. Andrieu等^[16]指出,引入MCMC的PF算法是计算后验分布较通用的方法.

本文提出一种确定模型参数的混合PF方法,记为PF-DREAM. 该方法采用DST理论分析确定先验模型参数;在计算后验分布时,将差分进化自适应Metropolis算法(differential evolution adaptive Metropolis, DREAM)引入PF采样算法中,一方面解决了高维积分问题,另一方面可保证采样过程中粒子样本的多样性,使得参数确定及物理模型更加合理. 本文方法可为系统的可靠性评估提供有效工具.

1 贝叶斯原理

考虑用物理模型$M({\pmb y};{\pmb x})$来描述对应的观测事件 d , x 为模型不确定参数, y 为模型独立变量. 如果模型绝对完美,则认为 M ( y ; x ) = d .但是事实上,由于模型构建时的简化、模型参数 x 的统计识别、观测事件 d 中的测量噪声等不确定性,使得完美的模型几乎难以实现.

假设模型参数 x 的先验概率分布为 p ( x | M ) ,似然函数或者条件概率分布为 p ( d | x , M ) ,则 x 的后验概率分布 p ( x | d , M ) 可以运用贝叶斯理论表示如下

$ p({\pmb x}\vert d,M) = \dfrac{1}{Z}p({\pmb x}\vert M)p(d\vert {\pmb x},M) \propto \\ p({\pmb x}\vert M)p(d\vert {\pmb x},M) $(1)

其中 Z = ∫ x p x M p d x , M d x 是一个正则化常数.

若 m 表示模型的不确定性, e 表示测量的不确定性. 假设两者相互独立且为具有零均值,标准差分别为 σ m 和 σ e 的正态分布随机变量,则条件概率分布就可以采用如下的高斯分布

$ p(d\vert {\pmb x},M) = \\ \dfrac{1}{\sqrt {2\pi \sigma _M^2 } }exp \Bigg(\dfrac{ - [d - M({\pmb y};{\pmb x})]^2}{2\sigma _M^2} \Bigg )(2)$

其中 σ M 2 = σ m 2 + σ e 2 .

将式(2)代入式(1)就可以得到整合了观测事件 d 的关于模型参数 x 的后验概率分布. 如前所述,后验概率分布的计算需要考虑两个问题：(1)初始模型参数 x 的先验概率分布 p ( x | M ) 会极大地影响模型前期的预测精度;(2)对于高维问题,式(1)中的正则化常数 Z 是一个多维积分,一般很难通过解析的方法计算得出.

2 确定先验分布的DST方法

为合理描述结构的真实物理响应,对于数据集较少的情况,可采用DST的混合联合准则来确定模型参数的先验概率分布^[9]. DST理论有以下假设：如果两个个体的证据一致,那么它们应该具有更高的信任值,在联合单个证据中被赋予更大的权重.

假设描述一批同类型结构性能退化的数据集有 n 个. 以一维模型参数为例,确定初始化模型参数的DST分析步骤如下:

(a)根据选定的模型 M ,利用MATLAB软件中的高斯--牛顿算法的曲线拟合工具来拟合第 i 条数据集,模型参数 x 的估计值为 x i ,95%置信区间为 B i ( i = 1,2 , ⋯ , n ).

(b)将 x i 的置信区间 B i 按顺序进行比较,计算与其相关的信任函数 Bel ( x i ) .

假设所有的数据集具有相同的信任度,则

$m(x_i ) = \dfrac{1}{n} (3)$

其中, m ( x i ) 表示第 i 条数据集的基本概率分配,即估计值 x i 的基本概率分配. 参数估计值 x i 的信任函数 Bel ( x i ) 定义为

Bel ( x i ) = ∑ a ll B j ⊆ B i m ( x j ) ( 4 )

(c)利用信任函数更新基本概率分配 m ( x i ) ,即

m * ( x i ) = ∑ a ll B j ⊆ B i ( - 1 ) B i - B j Bel ( x j ) ( 5 )

其中 B i - B j 表示两个集合的基数差.

由于集合是采用置信区间来定义的,因此基数一般是不可数的,需要用到勒贝格测度. 如果区间值较小,区间的长度也会较小,则集合 B i 与 B j 基数差可以取零. 因此,在特殊情况下,基本概率分配的赋值可以采用如下近似公式

m * ( x i ) = ∑ a ll B j ⊆ B Bel ( x j ) ( 6 )

将重新赋值的基本概率分配进行正则化处理,即

$\hat {m}(x_i ) = \dfrac{m^\ast (x_i )}{\sum _{j = 1}^n {m^\ast (x_j )} } $(7)

(d)模型参数先验分布的均值可以用加权算术平均值得到

u 0 = ∑ i = 1 n m ̂ ( x i ) ⋅ x i ( 8 )

模型参数先验分布的方差 Σ 0 根据估计值 x i ( i = 1,2 , ⋯ , n )统计得到.

以上是利用DST理论确定一维模型参数先验分布的步骤. 如果物理模型中具有 N 个模型参数,则对每个模型参数重复上述步骤.

3 基于PF算法的贝叶斯更新

在结构使用期间,通常会对结构进行健康监测及性能检测,因而获得一系列检测数据,通过贝叶斯原理对这些信息进行分析,可得模型参数的后验分布,从而使结构的物理模型更加准确. 在后验分布的计算中,正则化积分 Z 通常是一个多维复杂积分,很难直接求解得出. 利用粒子滤波方法,即用一系列带权值的粒子来逼近后验概率密度函数,可避免对 Z 的直接计算^[17].

在PF算法中,通过生成一系列的粒子实现贝叶斯更新,这些粒子带有关于模型参数的概率信息. 当有新的检测数据时,将前一步 的后验分布作为当前步的先验分布,结合似然函数实现参数的逐次更 新^[18]. 实际上,该方法是使用蒙特卡洛法仿真来完成一个递推的贝叶斯滤波,因此PF算法常被称为序列蒙特卡洛贝叶斯滤波法^[19].

在贝叶斯滤波中,目标是估计参数的后验概率分布 p ( X k | d 1 : k ) ,因此,可以通过预测和更新2个步骤递推完成^[20],PF采样过程如图1所示.

(i)预测：假设 d 1 : k = [ d 1 , d 2 , ⋯ , d k ] 表示直到 k 时刻的响应观测值, p ( X k - 1 | d 1 : : k - 1 ) 表示在 k - 1 时刻的后验概率分布. 利用Chapman-Kolmogorov方程计算参数 X k 在 k 时刻的先验概率分布

p ( X k | d 1 : k - 1 ) = ∫ p ( X k | X k - 1 ) ⋅ p ( X k - 1 | d 1 : k - 1 ) d X k - 1 ( 9 )

式中, p ( X k | X k - 1 ) 表示状态转移 X k - 1 → X k 的概率.

(ii)更新：在 k 时刻,当有新的观测值时,使用贝叶斯方法获得 X k 的后验分布

$p({\pmb X}_k \vert {\pmb d}_{1:k} ) = \dfrac{p({\pmb X}_k \vert {\pmb d}_{1:k - 1} ) \cdot p(d_k \vert {\pmb X}_k )}{p(d_k \vert {\pmb d}_{1:k - 1} )} $(10)

式中, p ( d k | X k ) 表示似然函数,正则化常数 p ( d k | d 1 : k - 1 ) = ∫ p ( X k | d 1 : k - 1 ) ⋅ p ( d k | X k ) d X k .

利用递推求解式(9)和式(10)就可以得到贝叶斯最优估计. 为了避免计算积分,PF算法利用一系列带权值的粒子来逼近概率密度函数,即

$p({\pmb X}_k \vert {\pmb d}_{1:k} ) \approx \sum _{i = 1}^N {w_k^i \cdot \delta ({\pmb X}_k - {\pmb X}_k^i)} $(11)

式中, ${\pmb X}_k^i$, $i = 1,2,\cdots,N$为取自$p({\pmb X}_k \vert {\pmb d}_{1:k} )$的独立随机样本;$\delta (\cdot )$为狄拉克函数;$w_k^i $是与每个样本${\pmb X}_k^i $相关的贝叶斯重要权值.

实际上, $p({\pmb X}_k \vert {\pmb d}_{1:k} )$的分布是未知的, 因此, ${\pmb X}_k^i$常常依赖于重要性采样, 从任意分布$\pi ({\pmb X}_k^i \vert {\pmb d}_{1:k} )$中采集样本点${\pmb X}_k^i$, $i =1,2,\cdot,N$. 而${\pmb w}_k $可以用下式进行计算

${\pmb w}_k \propto \dfrac{p({\pmb X}_k \vert {\pmb d}_{1:k} )}{\pi ({\pmb X}_k \vert {\pmb d}_{1:k} )}$(12)

式中, $\pi ({\pmb X}_k \vert {\pmb d}_{1:k} )$叫做重要性函数.

将权重因子正则化

$\dfrac{w_k^i }{\sum _{j = 1}^N {w_k^j } } \Rightarrow w_k^i , \ \ \ \ i =1,2,\cdots,N (13)$

权值的更新常采用序列重点采样方法(sequential important sampling, SIS),更新公式如下

$w_k^i = w_{k-1}^i \cdot \dfrac{p({ d}_k \vert {\pmb X}_k^i ) \cdot p({\pmb X}_k^i {\vert }{\pmb X}_{k - 1}^i )}{\pi ({\pmb X}_k^i \vert {\pmb X}_{k - 1 }^i,{\pmb d}_{1:k} )} $(14)

如果选择的重要性函数为$\pi ({\pmb X}_k^i \vert {\pmb X}_{k - 1}^i ,{\pmb d}_{1:k} ) = p({\pmb X}_k^i {\vert}{\pmb X}_{k- 1}^i )$, 那么权值的更新规则变为

$w_k^i = w_{k-1}^i \cdot p({\pmb d}_k \vert{\pmb X}_k^i ) $(15)

随着更新过程的不断迭代,只有少数粒子的权重起显著作用,其余粒子的权重几乎降低到零,这些粒子对后验概率分布没有贡 献,即发 生粒子的退化现象.

为了解决退化问题,Gordon等^[21]提出重采样的方法,其基本思想是权重较大的粒子被复制,而权重较小的粒子被剔 除^[22]. 本文采用逆CDF重采样算 法^[18,23],基于似然函数构建累积分布函数(CDF),通过生成服从均匀分布 u ( 0,1 ) 的随机数作为CDF值,并反过来选取CDF相对应的粒子,实现粒子重采样. 经过重采样之后,所有的粒子具有相同的权重值,即$w_k^i = \dfrac{1}{N}$ ($i = 1,2,\cdots,N)$. 因此后验概率密度函数可以表示如下

$p({\pmb X}_k \vert {\pmb d}_{1:k} ) \approx \dfrac{1}{N}\sum _{i = 1}^N {\delta ({\pmb X}_k - {\pmb X}_k^i)} $(16)

基于以上粒子滤波方法计算模型参数的后验概率密度函数,随着迭代的进行,大部分粒子将分布于高概率密度的区域,粒子样本 多样性降低,进而出现粒子匮乏现象^[24]. 本文,粒子匮乏判断准则如下：若群体中相异粒子数目占群体总数目的比例低于某值(如 r =0.4)时,则认为粒子出现匮乏.

4 混合PF-DREAM算法

针对PF算法在采样过程中粒子样本多样性减少的问题,本文引入差分进化自适应Metropolis算法(differential evolution adaptive Metropolis, DREAM)移动步,用转移核矩阵 κ ( ⋅ ) 在空间中传递粒子,即在重采样步骤后引入基于DREAM方法的移动步来增加粒子群的多样性,从而提高计算精度. 如果马尔科夫链的稳态分布为 p ( X k | d 1 : k ) ,那么传递得到的新粒子仍然服从后验分布,但是新的粒子会移动到空间中的其他区域,而不仅限于先验采样空间的区域, 从而增加粒子样本的多样性.

DREAM算法是一种多条马尔科夫链并行采样的方法,各条链之间通过差分进化遗传算法自适应调整建议分布,从而实现每条链 上的采样^[25,26]. 这种采样算法不仅弥补了MCMC算法中无法判断何时收敛的问题,而且还解决了提议分布极大依赖于选取函数的问题^[27]. DREAM满足细致平稳条件且具有遍历性,在高维、非线性、多峰等问题上都具有极好的采样性能^[28]. 由于DREAM采用多链并行计算,可以及时判断马尔科夫链何时收敛;另外,在自适应链中引入遗传算法,可使粒子在空间做状态转移时, 有更大的可能移动到空间的其他区域上,从而能够极大地增加样本多样性^[29,30,31].

细致平稳性条件定理：对于非周期马尔科夫链转移矩阵 κ ( ⋅ ) 和分布 π ( x ) ,若对于所有的 i , j 满足 π ( i ) κ ij = π ( j ) κ ji ,那么 π ( x ) 是马尔科夫链的平稳分布.

已知目标分布为$p(x) = \pi (x)$, 需要产生服从分布$p(x)$的样本点. 由细致平稳条件可知, 实现上述结果的核心是构造转移矩阵${\pmb\kappa} ( \cdot )$, 使得$p(i) \cdot \kappa _{ij} = p(j) \cdot \kappa _{ji} $.假设已知一个转移矩阵${\pmb Q}$, 马氏链$q(i,j)$表示从状态$i$转移到状态$j$的概率, 通常情况下, $p(i) \cdot q(i,j) \ne p(j) \cdot q(j,i)$, 引入$\alpha (i,j) = p(j) \cdot q(j,i)$和$\alpha (j,i) = p(i) \cdot q(i,j)$, 可得

$p(i) \cdot q(i,j) \cdot \alpha (i,j) = p(j) \cdot q(j,i) \cdot \alpha(j,i) $(17)

令$\kappa _{ij} = q(i,j) \cdot \alpha (i,j)$, $\kappa _{ji} = q(j,i) \cdot \alpha (j,i)$, 则转移矩阵${\pmb\kappa} ( \cdot )$的平稳分布为$p(x)$.

混合PF-DREAM算法采样步骤如下：

(1)初始化马尔科夫链初始状态. 假设通过PF重采样后得到的后验目标分布为 p ( X ) , X 为 L 维模型参数. 初始产生 N 条链的起始点 X 0 = { X 0 1 , X 0 2 , ⋯ , X 0 N } ,假设每条链的长度为 T .

(2) 对 t = { 2,3 , ⋯ , T } 循环以下采样步骤.

对第 i 条链( i = 1,2 , ⋯ , N )实行以下循环.

(i)在初始空间中任意生成一个 L * 维的子集空间 A , R L * ⊆ R L ,链的移动步长 d X i 可以利用前一步的链合集 X t - 1 = X t - 1 1 , X t - 1 2 , ⋯ , X t - 1 N 的差分进化计算得到

$\left.\begin{array}{ll} d{\pmb X}_A^i = {\pmb \varsigma }_{L^\ast } + ({\pmb I }_{L^\ast } + {\pmb \lambda }_{L^\ast } ) \cdot \gamma (\delta ,L^\ast ) \cdot \sum _{j = 1}^\delta {({\pmb X}_A^{a_j } - {\pmb X}_A^{b_j } )} \\ d{\pmb X}_{ \ne A}^i = {\bf 0} \end{array} \right\} $(18)

其中, d X A i ≠ 0 ,跳跃率$\gamma = {2.38}/{\sqrt {2\delta L^\ast } }$, δ 表示用于产生跳跃的链的对数, a , b 是从 { 1,2 , ⋯ , i - 1 , i + 1 , ⋯ , N } 中组成链对 δ 的向量. λ 和 ς 为相互独立的随机变量, λ ~ U L * ( - c , c ) 和 ς ~ N L * ( 0 , c * ) ; I L * 表示 L * 维单位矩阵.

(ii)第 i 条链上的候选点为

X p i = X t - 1 i + d X i ( 19 )

(iii)采用Metropolis准则判断是否接受候选点

p a cc ( X t - 1 i → X p i ) = min [ 1 , p ( X p i ) / p ( X t - 1 i ) ] ( 20 )

实现准则如下

$\begin{array}{ll} {\rm if} \quad p_{\rm acc} ({\pmb X}_{t - 1}^i \to {\pmb X}_p^i ) \geqslant U(0,1) \\ {\rm then}\quad {\pmb X}_t^i = {\pmb X}_p^i \\ {\rm else}\quad {\pmb X}_t^i = {\pmb X}_{t - 1}^i \\ {\rm end} \end{array} $

通过上述步骤获得自适应马尔科夫链,当马尔科夫链 收敛后,采集的样本点服从目标分布^[18], 从而实现粒子的状态转 移,并且新粒子的多样性远优于旧粒子. 当出现下一个检测点用于模型更新时,后验分布结果和物理模型会更加精确合理.

5 PF-DREAM实现步骤及流程图

基于以上理论,PF-DREAM模型更新方法的流程如图2所示,程序实现的详细步骤如下:

(1) 根据选定的物理模型 M ,利用MATLAB中的高斯--牛顿算法的曲线拟合工具拟合先验数据集,得到每条曲线关于模型参数 的估计值及95%的置信区间;

(2) 利用DST联合准则确定模型参数的分布,作为先验分布 p ( x 0 ) ;

(3) 初始更新步 k = 1 ,从先验分布中生成等权值的 N 个粒子： { x i , 0 } i = 1 N ~ p ( x 0 ) { w i , 0 } i = 1 N = 1 / N ;

(4) 如果出现检测点,利用PF算法更新并执行重采样步骤得到粒子群： x i , k P i = 1 N , w i , k P i = 1 N = 1 / N ;若粒子出现匮乏现象,执行步骤(5);否则,执行步骤(6);

(5) 执行DREAM移动步,得到 N 个全新的粒子： x i , k * i = 1 N , w i , k * i = 1 N = 1 / N ;

(6) 计算参数后验估计值： x ̂ k = ∑ i = 1 N w i , k * ⋅ x i , k * ;

(7) 如果 k > K ( K 为检测数据点的个数),则停止;否则,令 k = k + 1 并返回步骤(4).

6 混合粒子滤波方法应用算例

6.1 电池退化性能预测

本算例,针对4组电池的容量退化,将其中3组数据作为先验数据信息,预测第4组电池容量退化轨迹. 数据来源说明如下： 对象为石墨阳极锂氧化钴阴极额定容量为1.1 Ah的商业锂离子电池,采用美国Arbin BT2000的锂电池实验系统在室温条件下进行循环充放电老化实验,放电电流为1.1 A,其中4组电池B1、B2、B3、B4的容量退化数据如图3中绿色十字符所示,数据来自于马里兰大学先进寿命周期工程中心.

电池的退化模型选择双指数模型^[9],即

$Q = a \cdot \exp (b \cdot N) + c \cdot \exp (d \cdot N) $(21)

式中, a , c 代表内阻抗, b , d 代表退化速率.

失效阈值取为额定电容的0.8倍. 为了方便起见,假设联合分布 x = [ a , b , c , d ] T 是一个多元正态分布. 基于B1、B2、B3数据集,采用DST理论确定模型参数的先验分布,似然函 数中的标准差采用最大似然估计(MLE)计算,得到的 [ a , b , c , d ] 的联合分布如下

$\begin{array}{ll} p_0 ({\pmb x}) = p_0 (a,b,c,d) =\dfrac{1}{\sqrt {(2\pi )^n\left| {{\pmb \varSigma }_0 } \right|} }\cdot \\ \exp\Big ( - \dfrac{1}{2} \left[ {{\pmb x} - {\pmb \mu }_0 } \right]^{\rm T}{\pmb \varSigma }_0 ^{ - 1}\left[ {{\pmb x} - {\pmb \mu }_0 } \right] \Big) \end{array} $(22)

式中

${\pmb \mu}_0 = \left[ \!\!\begin{array}{c} { - 2.643\times 10^{ - 3}} \\ {7.101\times 10^{ - 1}}\\ {1.098} \\ { - 2.397\times 10^{ - 2}} \end{array} \!\! \right] $

${\pmb\varSigma}_0 = \left[ \!\!\begin{array}{cccc} {6.033\!\times\!10^{ - 2}} & {8.714\!\times\!10^{ - 5}} & {7.355\!\times\!10^{ - 6}} & {7.171\!\times\!10^{ - 6}} \\ {8.714\!\times\!10^{ - 5}} & {3.661\!\times\!10^{ - 1}} & {1.135\!\times\!10^{ - 4}} & {1.107\!\times\!10^{ - 4}} \\ {7.355\!\times\!10^{ - 6}} & {1.135\!\times\!10^{ - 4}} & {1.559\!\times\!10^{ - 1}} & {9.345\!\times\!10^{ - 6}} \\ {7.171\!\times\!10^{ - 6}} & {1.107\!\times\!10^{ - 4}} & {9.345\!\times\!10^{ - 6}} & {6.253\!\times\!10^{ - 2}} \end{array}\!\! \right] $

在B4退化曲线中,每10个循环周期采集一个数据点,共取30个退化数据 Q i , N i i = 1 n = 30 ,这30个点作为实际观测值用来更新模型参数, [ a , b , c , d ] 的后验分布如下

$ p(a,b,c,d) \propto p_0 (a,b,c,d) \cdot \\ \exp \Bigg \{ - \dfrac{1}{2}\sum _{i = 1}^n {\left[ {\dfrac{Q_i - M(N_i;a,b,c,d)}{\sigma _M }} \right]^2}\Bigg \} $(23)

式中, M 为电池退化模型,似然函数中的标准差为 σ M = 0.01 .基于观测值 Q i , N i i = 1 n , n = 30 ,分别采用PF算法和混合PF-DREAM算法 对模型参数 [ a , b , c , d ] 进行更新,其均值和标准差,见表1和表2所示(每10次更新结果).

由以上计算结果可知,模型参数 [ a , b , c , d ] 的标准差随着检测点的增加而逐渐减小.表1结果表明,随着更新点数的增加, 模型参数的均值变化均不大,而标准差减小. 这是由于PF算法在每次更新的过程中会删除低权值的粒子,复制高权值的粒子,随着迭代更新次数的增加,样本集的多样 性大大降低,之后虽然样本集的方差减小,但样本群已经不具有一般性;表2结果表明,随着更新点数的增加,各参数均 值的变化较为明显,并且标准差随之减小. 由此说明混合PF-DREAM在计算后验分布时,能够显 著提高样本群的多样性,使后验分布计算结果更加准确.

为了形象地描述更新后的模型参数对响应预测结果的准确性,将基于 Q i , N i i = 1 n , n = 0,10,20,30 数据点更新后的物理模型预测轨迹均值绘制于图4、图5.图4为采用PF方法更新后的响应预测,图5为采用本文 方法更新后的响应预测.

图4和图5中的0点更新预测曲线(实线)是基于先验分布而绘制的曲线,初始阶段较为平缓. 可以看出,在初始阶段具有检测点后,更新后的预测曲线在初始阶段变得陡峭. 这是因为利用DST理论整合了所有的先验数据信息,得到的结果具有一般性;而B4组电池在初始阶段的容量退化趋势较快,因此更新后的模型在初始阶段更加接近实际退化轨迹.

图4中,利用 n = 10 检测点更新后的预测轨迹(点划线)与 n = 20 检测点更新后的预测轨迹(虚点线)基本上重合,且二者与 n = 30 检测点更新后预测的轨迹(虚线)变化也不大,但3条预测轨迹与实际退化轨迹误差较大.图5中, n = 10 检测点更新后预测轨迹(点划线)在充放电循环周期的前期阶段与实际退化轨迹较为接近,但在预测后期误差较大. 随着更新点数的增加(即 n =20虚点线, n =30虚线),预测轨迹与实际退化轨迹之间误差减小,表明混合PF-DREAM更新的预测结果比PF算法的预测结果更加接近实际,说明PF-DREAM更新方法能够提高物理模型的准确性.

6.2 疲劳裂纹扩展预测

本算例,采用12条裂纹扩展退化数据曲线,以其中10条退化数据作为先验信息,取curve 1 和curve 2前期部分数据作为实际观测数据,分别预测curve 1 和curve 2的退化轨迹. 相关数据如图6所示.

数据来源说明如下^[32,33]：选用的试验数据采用Wu和 Ni 所做的疲劳裂纹扩展试验,疲劳裂纹试验件为 2024-T351 铝合金板件. 试件的宽度和厚度分别为 b = 50 mm 和 t = 12 mm,预处理裂纹使其长度扩展至 18 mm. 在裂纹预 处理和疲劳裂纹扩展试验过程中采用正弦曲线信号,应力幅 Δ σ = 56.33 MPa.

裂纹的退化模型选择Paris模型

$\dfrac{da}{dN} = C \cdot \left( {\Delta K} \right)^m (24)$

式中, $ \Delta K = \Delta \sigma \cdot \sqrt { \uppi \cdot a \cdot {\rm sec}( \uppi \cdot a / b)}$, a / b < 0.7 , C , m 为材料参数. 一般为了简化计算,取 l n c 来替代^[13-14],并且假设联合分布 [ l n c , m ] 服从二元正态分 布^[13].

取其中的10条曲线作为先验信息集,如图6中的实线,另取两条曲线作为实际退化曲线,用曲线1,2表示. 为了实现这两条退化曲线的预测,分别在两条曲线上各取3个点 ( a i , N i ) , i = 1,2 , 3 作为实际的检测数据. 曲线1的检测数据分别为(18.912, 10 000)、(19.895, 20 000)和(21.088, 30 000);曲线2的检测数据分别为：(18.737, 10 000),(19.239, 20 000)和(20.186, 30 000).

基于先验数据集,采用DST理论确定模型参数的先验分布,似然函数中的标准差采用最大似然估计(MLE)计算,得到的 [ l n c , m ] 的联合分布如下

$\begin{array}{ll} p_0 ({\rm ln} c,m) = p_0 ({\pmb x}) = \\ \dfrac{1}{\sqrt {(2\pi )^n\left| {\pmb \varSigma } \right|} } \exp \Big( - \dfrac{1}{2}[{\pmb x} - {\pmb u}_0 ]^{\rm T} {\pmb\varSigma}^{ - 1}[{\pmb x} - {\pmb u}_0 ] \Big) \end{array} $(25)

其中

${\pmb\varSigma}= \left[ \!\!\begin{array}{cc} {1.409} & { - 0.1631} \\ { - 0.1631} & {3.132\times 10^{ - 4}}\end{array} \!\! \right] , \ \ {\pmb \mu }_{0} = \left[\!\!\begin{array}{c} {-46.39} \\ {5.998} \end{array}\!\!\right] $

利用检测数据分别更新模型参数, [ ln c , m ] 的后验分布如下

$ P (\ln c,m) \propto p_0 (\ln c,m) \cdot \\ \exp \left\{ - \dfrac{1}{2}\sum _{i = 1}^n \left[\dfrac{a_i - M(N_i ;\ln c,m)} {\sigma _M } \right]^2 \right\} $(26)

式中,似然函数中的标准差取为统计平均值 σ M = 0.304 .

分别基于曲线1,曲线2的观测值 ( a i , N i ) , i = 1,2 , 3 ,利用PF及混合PF-DREAM算法对参数 [ ln c , m ] 进行更新,结果如表3和表4所示(模拟样本点个数 N = 100 000 ).

从以上结果可以看出,随着更新点的增加,模型参数后验分布的统计标准差逐渐减小,相比于PF,PF-DREAM方法得到的结果, 其分散性更小.

若板件的临界裂纹尺寸 a c = 30 mm,则极限状态方程为

g = a c - M ( N ; ln c , m ) ( 27 )

对应于 N = 50 000 时,先验失效概率为 P f = 0.257 4 ;如基于曲线1数据进行3次更新,则PF¹和PF-DREAM²计 算得到的失效概率分别为 P f 1 = 0.494 1 , P f 2 = 0.524 8 .

为表示更新后的模型参数对裂纹扩展预测结果的影响,图7 ~ 图10为采用PF算法和PF-DREAM算法预测的中值曲线.

基于先验分布的裂纹扩展预测轨迹(如图中0点更新预测曲线),大体位于两边缘曲线之间,具有一般性. 从图中可以看出, 在前期只有一个检测点用于更新时,后验分布预测结果与先验分布(0个检测点)预测结果相差不大,都离实际裂纹扩展轨迹较远. 当继续增加检测点时,后验分布的预测结果越来越靠近实际轨迹. 在相同的更新点数(2个或3个检测点)的情况下,PF-DREAM算法预测的结果比传统的PF的结果更加接近实际退化轨迹. 原因在于PF算法在采样过程中保留并复制高权值粒子,删除低权值粒子,使得样本的多样性减少;而PF-DREAM算法通过引入 DREAM使粒子的多样性大大增加,因此得到的更新结果会更加准确. 通过此算例表明,PF-DREAM方法能够有效提高退化结构性能预测的准确性.

7 结 论

本文提出了一种混合粒子滤波方法来更新模型参数,使得描述结构退化的物理模型更加准确. 在确定模型参数的先验分布时,采用DST理论以充分利用数据集中的信息. 针对后验分布计算中的高维积分问题,本文提出的混合PF-DREAM方法相比于传统的PF算法,能够得到更加合理的物理模型. 此外,混合PF-DREAM可以有效提高样本粒子的多样性,解决重采样算法中粒子多样性匮乏的问题.

电池退化以及裂纹扩展的算例结果表明,PF-DREAM方法能够有效提高预测精度. 基于以上结果,可以得出以下几条结论：

(1) 结构未投入使用且先验的信息集较少的情况下,可以采用DST理论的混合联合法则分配权重,再利用统计方法得到模型参数的初始化先验分布.

(2) 结构投入使用后,获得的检测信息可以用来修正模型参数. 电池退化及裂纹扩展的算例中表明,随着用于更新的检测数据点的增加,预测结果也更加接近真实的退化轨迹.

(3) PF-DREAM算法解决了粒子多样性匮乏的问题,使更新结果更加准确合理.

此外,需要指出的是,当参数数目增加或需要处理隐式函数问题时,本文方法的应用步骤不发生变化,但数值运算次数随之增加,导致计算成本增加.

The authors have declared that no competing interests exist.

---

参考文献

文献选项

• 原文顺序
• 文献年度倒序
• 文中引用次数倒序
• 被引期刊影响因子

[1]	Tsai CC . Mis-specification analyses of gamma and Wiener degradation processes . Journal of Statistical Planning& Inference, 2011 , 141 ( 12 ): 3725 - 3735 [本文引用: 1]
[2]	Dan MF , Kallen MJ , Noortwijk JMV . Probabilistic models for life-cycle performance of deteriorating structures: Review and future directions . Progress in Structural Engineering& Materials, 2004 , 6 ( 4 ): 197 - 212 [本文引用: 1]
[3]	刘海波 , 姜潮 , 郑静 等 . 含概率与区间混合不确定性的系统可靠性分析方法 . 力学学报, 2017 , 49 ( 2 ): 456 - 466 [本文引用: 1] ( Liu Haibo , Jiang Chao , Zheng Jing , et al . A system reliability analysis method for structures with probability and interval mixed uncertainty . Chinese Journal of Theoretical and Applied Mechanics, 2017 , 49 ( 2 ): 456 - 466 (in Chinese)) [本文引用: 1]
[4]	Ge R , Chen J , Wei J . Reliability-based design of composites under the mixed uncertainties and the optimization algorithm . Acta Mechanica Solida Sinica, 2008 , 21 ( 1 ): 19 - 27 [本文引用: 1]
[5]	李帅兵 , 杨睿 , 罗喜胜 等 . 气流作用下同轴带电射流的不稳定性研究 . 力学学报, 2017 , 49 ( 5 ): 997 - 1007 [本文引用: 1] ( Li Shuaibing , Yang Rui , Luo Xisheng , et al . Instability study of an electrified coaxial jet in a coflowing gas stream . Chinese Journal of Theoretical and Applied Mechanics, 2017 , 49 ( 5 ): 997 - 1007 (in Chinese)) [本文引用: 1]
[6]	周春晓 , 汪锐琼 , 聂肇坤 等 . 基于最大熵方法的水下航行体结构动力响应概率建模 . 力学学报, 2018 , 50 ( 1 ): 114 - 123 [本文引用: 1] ( Zhou Chunxiao , Wang Ruiqiong , Nie Zhaokun , et al . Probabilistic modelling of dynamic response of underwater vehicle structure via maximum entropy method . Chinese Journal of Theoretical and Applied Mechanics, 2018 , 50 ( 1 ): 114 - 123 (in Chinese)) [本文引用: 1]
[7]	Huang X , Chen J . Time-dependent reliability model of deteriorating structures based on stochastic processes and Bayesian inference methods . Journal of Engineering Mechanics, 2015 , 141 ( 3 ): 04014123 [本文引用: 1]
[8]	Guan X , He J , Jha R , et al . An efficient analytical Bayesian method for reliability and system response updating based on Laplace and inverse first-order reliability computations . Reliability Engineering& System Safety, 2012 , 97 ( 1 ): 1 - 13 [本文引用: 2]
[9]	He W , Williard N , Osterman M , et al . Prognostics of lithium-ion batteries based on Dempster-Shafer theory and the Bayesian Monte Carlo method . Journal of Power Sources, 2011 , 196 ( 23 ): 10314 - 10321 [本文引用: 3]
[10]	Rebba R , Mahadevan S . Computational methods for model reliability assessment . Reliability Engineering& System Safety, 2008 , 93 ( 8 ): 1197 - 1207 [本文引用: 1]
[11]	Luke T , Joseph BK . Accurate Approximations for posterior moments and marginal densities . Journal of the American Statistical Association, 1986 , 81 ( 393 ): 82 - 86 [本文引用: 1]
[12]	王宇 . 贝叶斯参数更新在可靠性分析中的应用. [硕士论文] . 南京:南京航空航天大学, 2014 [本文引用: 1] ( Wang Yu . The application of Bayesian updating in reliability analysis. [Master Thesis] . Nanjing: Nanjing University of Aeronautics and Astronautics, 2014 (in Chinese)) [本文引用: 1]
[13]	Cappé O , Moulines E , Rydén T . Inference in Hidden Markov Models . Springer, 2005 : 574 - 575 [本文引用: 3]
[14]	Pitt MK , Silva RDS , Giordani P , et al . On some properties of Markov chain Monte Carlo simulation methods based on the particle filter . Journal of Econometrics, 2012 , 171 ( 2 ): 134 - 151 [本文引用: 2]
[15]	Fernándezvillaverde J , Rubioramrez JF . Estimating macroeconomic models: A likelihood approach . Review of Economic Studies, 2010 , 74 ( 4 ): 1059 - 1087 [本文引用: 1]
[16]	Christophe A , Arnaud D , Roman H . Particle Markov chain Monte Carlo methods . Journal of the Royal Statistical Society, 2010 , 72 ( 3 ): 269 - 342 [本文引用: 1]
[17]	Doucet A , Godsill S , Andrieu C . On sequential Monte Carlo sampling methods for Bayesian filtering . Statistics& Computing, 2000 , 10 ( 3 ): 197 - 208 [本文引用: 1]
[18]	An D , Choi JH , Kim NH . Prognostics 101: A tutorial for particle filter-based prognostics algorithm using Matlab . Reliability Engineering& System Safety, 2013 , 115 ( 1 ): 161 - 169 [本文引用: 3]
[19]	王鑫 , 胡昌华 , 暴飞虎 . 基于贝叶斯原理的粒子滤波算法 . 弹箭与制导学报, 2006 , 26 ( s5 ): 269 - 271 [本文引用: 1] ( Wang Xin , Hu Changhua , Bao Feihu . Particle filtering method based on Bayesian theorem . Journal of Projectiles, Roktets, Missiles and Guidance, 2006 , 26 ( s5 ): 269 - 271 (in Chinese)) [本文引用: 1]
[20]	Chen J , Yuan S , Qiu L , et al . On-line prognosis of fatigue crack propagation based on Gaussian weight-mixture proposal particle filter . Ultrasonics, 2017 , 82 : 134 [本文引用: 1]
[21]	Gordon NJ , Salmond DJ , Smith AFM . Novel approach to nonlinear/non-Gaussian Bayesian state estimation . IEE Proceedings F-Radar and Signal Processing, 2002 , 140 ( 2 ): 107 - 113 [本文引用: 1]
[22]	Djurić PM , Kotecha JH , Zhang J , et al . Particle filtering . Signal Processing Magazine IEEE, 2003 , 20 ( 5 ): 19 - 38 [本文引用: 1]
[23]	Zio E , Peloni G . Particle filtering prognostic estimation of the remaining useful life of nonlinear components . Reliability Engineering& System Safety, 2011 , 96 ( 3 ): 403 - 409 [本文引用: 1]
[24]	Bi H , Ma J , Wang F . An improved particle filter algorithm based on Ensemble Kalman filter and Markov Chain Monte Carlo method . IEEE Journal of Selected Topics in Applied Earth Observations& Remote Sensing, 2015 , 8 ( 2 ): 447 - 459 [本文引用: 1]
[25]	Vrugt JA , Braak CJFT , Clark MP , et al . Treatment of input uncertainty in hydrologic modeling: Doing hydrology backward with Markov chain Monte Carlo simulation . Water Resources Research, 2008 , 44 ( 12 ): 5121 - 5127 [本文引用: 1]
[26]	Vrugt JA , Braak CJFT , Diks CGH , et al . Accelerating Markov chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling . International Journal of Nonlinear Sciences& Numerical Simulation, 2009 , 10 ( 3 ): 273 - 290 [本文引用: 1]
[27]	Laloy E , Vrugt JA . High-dimensional posterior exploration of hydrologic models using multiple‐try DREAM(ZS) and high-performance computing . Water Resources Research, 2012 , 50 ( 3 ): 182 - 205 [本文引用: 1]
[28]	Vrugt J A . Markov chain Monte Carlo simulation using the DREAM software package: Theory, concepts, and MATLAB implementation . Environmental Modelling& Software, 2016 , 75 : 273 - 316 [本文引用: 1]
[29]	Wentworth MT , Smith RC , Williams B . Bayesian model calibration and uncertainty quantification for an HIV model using adaptive Metropolis algorithms . Inverse Problems in Science& Engineering, 2017 , 26 ( 2 ): 233 - 256 [本文引用: 1]
[30]	Post H , Vrugt JA , Fox A , et al . Estimation of community land model parameters for an improved assessment of net carbon fluxes at European sites . Journal of Geophysical Research Biogeosciences, 2017 , 122 ( 3 ): 661 - 689 [本文引用: 1]
[31]	Liang S , Jia H , Xu C , et al . A Bayesian approach for evaluation of the effect of water quality model parameter uncertainty on TMDLs: A case study of Miyun reservoir . Science of the Total Environment, 2016 , s560-561 : 44 - 54 [本文引用: 1]
[32]	Wu WF , Ni CC . A study of stochastic fatigue crack growth modeling through experimental data . Probabilistic Engineering Mechanics, 2003 , 18 ( 2 ): 107 - 118 [本文引用: 1]
[33]	Wu WF , Ni CC . Statistical aspects of some fatigue crack growth data . Engineering Fracture Mechanics, 2007 , 74 ( 18 ): 2952 - 2963 [本文引用: 1]