CN108595736B - 一种机构可靠性建模方法 - Google Patents

Abstract

本发明涉及一种机构可靠性建模方法，主要通过四个步骤实现：第一步，通过对样本数据使用随机过程拟合得到退化模型；第二步，对退化特征量进行时变相关性分析，使用D‑VineCopula函数进行多元相关性建模，时变相关性使用Kendall系数时变模型描述；第三步对相关特征量进行抽样；第四步，建立失效模式的极限状态方程，将所有样本带入功能函数进行分析，并计算得到可靠度。本发明基于所提出的多元时变相关性退化建模方法，解决了现有建模方法由于将多特征量独立假设导致机构可靠度结果误差偏大的缺陷。同时本发明可以灵活构造多个退化特征量之间的相关关系，对相关关系进行时变建模，得到的可靠性评估结果更为准确。

Description

一种机构可靠性建模方法

技术领域

本发明涉及机械产品可靠性设计与分析领域，尤其涉及机构功能可靠性设计中一种机构可靠性建模方法。

背景技术

随着精密机械工程和航空航天工程的不断发展，机械机构系统向着高精度、高效率、高可靠性和长寿命的目标迈进。机构是机械系统的重要组成部分，机构中的零部件通过铰链相互连接，共同传递运动和载荷，以实现机构的功能，因此机构的功能可靠性成为研究的重点。随着机构的运行，机构的退化和损伤成为影响机构性能的主要因素之一，如机构中运动副的磨损。同时一个机构的多个部件在相同的环境下工作，受到相同的影响因素，各个部件之间的损伤与退化相互影响，它们之间是相关的，随着机构运动周期的增加，这种影响愈加强烈，对机构性能的影响也就越大。

现有技术中，王华伟等提出了多失效模式下多模型集成航空发动机运行可靠性评估方法(CN 103778295)，对性能退化失效、结构强度失效和突发失效进行了分析，并对三种失效选择最优模型，建立了多模型集成航空发动机运行可靠性评估模型，采用马尔科夫链蒙特卡洛仿真求解可靠度；张志鹏等提出了一种面向产品多维相关性退化失效的建模方法(CN 105117550)，利用D-VineCopula函数建立了多维相关退化模型，使用Wiener过程进行了退化建模，有效处理了具有复杂相关性的多维退化问题；李绍军等基于C-VineCopula函数进行了多模态的相关性建模(CN 104914775)，采用马尔科夫蒙特卡洛方法对联合密度函数进行抽样，利用密度分位数法等实现故障实时监控；谷东伟等提出了一种Copula函数可靠性动态影响度分析方法(CN 105700474)，考虑了数控机床子系统之间的相关关系，采用人工鱼群算法进行了相关系数优化，提出了采用可靠性动态核心影响度衡量子系统的改进潜力。

然而，上述现有技术的技术方案均是基于变量间的相关性不随时间变化开展的，随着对相关性研究的深入，变量间的相关性并不是不变的，而是随着时间动态变化的。机构退化过程中，工作环境和退化速率并不是恒定不变的，退化特征量之间的相关性也可能随时间不断变化。因此，在建立多元退化相关模型时，有必要考虑这些退化特征量之间相关关系随时间变化的时变特征。

发明内容

针对上述问题，本发明提供了一套能够考虑机构多部位退化之间的相互影响，同时考虑这些相关影响的时变特征的机构可靠性分析方法。

具体地，本文提出的方法主要通过四个步骤实现：第一步，通过对样本数据使用随机过程拟合得到退化模型；第二步，对退化特征量进行时变相关性分析，使用D-VineCopula函数进行多元相关性建模，时变相关性使用Kendall系数时变模型描述；第三步对相关特征量进行抽样；第四步，建立失效模式的极限状态方程，将所有样本带入功能函数进行分析，并计算得到可靠度。

本发明提出的机构可靠性分析方法包括以下详细步骤：

步骤1：确定机构的n个退化特征量，通过退化试验收集试验数据y＝(y₁,y₂,...,y_n)；

步骤2：利用步骤1得到的特征量的退化数据，在不考虑各个特征量之间的相关性情况下，利用随机过程，如Wiener过程或Gamma过程，建立每一个特征量的退化模型；

步骤3：利用马尔科夫链蒙特卡洛方法(MCMC)估计步骤2中退化模型的参数；

步骤4：建立多元时变相关模型，分两步进行：

步骤4(a)：建立多元相关模型：利用D-VineCopula建立n个特征量的相关性模型，将n元分布分解为多个二元相关的形式；

步骤4(b)：建立时变相关模型：使用基于Kendall系数τ的时变相关性模型描述特征量间的时变相关性；

步骤5：对步骤4的时变相关模型进行参数估计，具体包括以下五个步骤：

步骤5(a)：选取工程中常用的二元Copula函数：Normal Copula、ArchimedeanCopula,作为Copula函数备选集；

步骤5(b)：确定多元相关Copula函数类型：将步骤3得到的特征量退化参数带入步骤4(a)多元相关模型，用最大似然估计法得到每个备选Copula函数中的参数值，利用AIC准则从备选函数中选取最优Copula，作为第一个二元Copula；

步骤5(c)：重复步骤5(a)和5(b)，直到确定步骤4中模型的所有二元Copula函数类型；

步骤5(d)：使用最小二乘法求解步骤4(b)模型的参数；

步骤5(e)：根据Kendall系数τ与Copula参数θ之间的转换关系，根据步骤5(c)得到的Copula函数的类型，将τ(t)转换为θ(t)，由此得到每一个Copula函数的时变相关系数值；

步骤6：根据步骤5得到的各Copula类型和时变相关系数，基于Rosenblatt变换方法，生成随机向量的m组样本，具体包括以下五个步骤：

步骤6(a)：令u₁,u₂,...,u_n分别表示n维连续随机变量X＝(X₁,X₂,...,X_n)的边缘分布函数，即u₁＝F₁(x₁),...,u_n＝F_n(x_n)，则x₁＝F₁ ^-1(u₁),...,x_n＝F_n ^-1(u_n)。假设n维随机变量X的一个样本为(x₁,x₂,...,x_n)，r₁,r₂,...,r_n分别为[0,1]之间均匀分布的随机变量的一个样本，根据Rosenblatt变换，令u₁＝r₁，可得：x₁＝F₁ ^-1(u₁)；

步骤6(b)：由

得

则有

步骤6(c)：由

得

则有

步骤6(d)：由

可得

则有

步骤6(e)：以此类推得到

由此得到随机向量X的一组样本(x₁,x₂,...,x_n)；

上述步骤6中，F(·)为随机变量的边缘分布函数，F^-1(·)为边缘分布函数的逆函数，F_i|j(x_i|x_j)为随机变量x_j条件下随机变量x_i的条件分布函数；h_ij(i,j)为条件Copula函数。

步骤7：重复步骤6，对变量u₁,u₂,...,u_n进行m组随机抽样，共得到m组随机向量X，即X_m×n；

步骤8：建立失效模型的极限状态方程，对隐式的极限状态方程，使用响应面拟合各失效模式的极限状态方程，使用最小二乘法进行参数估计；

步骤9：将步骤7得到的样本带入步骤8的极限状态方程，使用蒙特卡洛法求解可靠度，落入失效域F内的样本点的个数N_um与总样本点的个数m之比即为失效概率的估计值P_f，即

步骤4(a)中的D-VineCopula模型的解析表达式如下：

其中，f(x₁,x₂,...,x_n)为n维特征量的联合概率密度函数，f_k(x_k)为第k个特征量的边缘概率密度函数，c_{i,i+j|i+1,...,i+j-1}为相应的二元Copula函数，θ_{i,i+j|i+1,...,i+j-1}为相关系数。

图2为基于D-VineCopula函数的一个四元可靠性退化特征量的树形表示，此时的联合概率密度函数为：

f(x₁,x₂,x₃,x₄)＝f₁(x₁)·f₂(x₂)·f₃(x₃)·f₄(x₄)·c₁₂·c₃₄·c_1,3|2·c_2,4|3·c_1,4|23 (2)

步骤4(b)的基于Kendall系数的时变相关性模型为：

τ(t)＝L·[sin(2π·(1+γ·e^-ξt))]^k+ψ (3)

其中，L,γ,ξ,k,ψ均为未知参数，L代表饱和度，γ为微调因子，ξ为速率系数，k是形状参数，ψ为稳定值，所有未知参数通过最小二乘回归方法获得，t→∞,τ(t)→ψ。

步骤6的Rosenblatt变换公式为：

步骤8中对于基本变量为X的系统，第k个失效模式的响应面函数为：

其中，a₀,a_ij,a_ii均为未知参数，通过最小二乘回归进行参数估计。

本发明基于所提出的多元时变相关性退化建模方法，解决了现有建模方法由于将多特征量独立假设导致机构可靠度结果误差偏大的缺陷。同时本发明可以灵活构造多个退化特征量之间的相关关系，对相关关系进行时变建模，得到的可靠性评估结果更为准确。

附图说明

图1为本发明蒙特卡洛求解系统可靠度的流程图

图2为基于D-VineCopula函数的一个四元可靠性退化特征量的树形表示法示意图

图3是舱门锁试验机构的平面示意图

图4各旋转副累积磨损量图

图5试验机构的D-VineCopula模型

图6试验机构可靠度

附图标记：1-摇臂；2-连杆；3-导轨；4-驱动滑块；5-内部加载弹簧；6-连杆AB；7-连杆BC；8-连杆BD；9-位移滑块；10-外部加载弹簧。

具体实施方式

下面结合附图，对实施例作详细说明。

以一种飞机舱门锁试验机构为例，图3是从舱门上位锁试验机构的平面示意图，利用位移滑块9的运动模拟锁钩运动形式，采用伺服电机和曲柄滑块机构驱动滑块4在导轨上做竖直方向的运动。通过连杆组件传递到位移滑块9上，使其沿导轨作水平方向的运动。连杆组件是一套平面连杆机构，用于将驱动滑块4上的驱动力传递到位移滑块9上，并保证两个滑块之间运动关系相协调。用外部加载弹簧10模拟作用在位移滑块9上的外部载荷。

步骤1：提取退化特征量，收集试验数据；

将摇臂1旋转一周作为一个耗损周期，试验平台中，共考虑4处旋转副的磨损，由于点B处是一个复合旋转副，该旋转副包括两个旋转副。将驱动滑块4与连杆6的旋转副定义为A，将连杆6与连杆7之间的旋转副定义为旋转副B1，连杆6与旋转副8之间的旋转副定义为旋转副B2，将连杆7与位移滑块9之间的旋转副定义为旋转副C，各随机变量初始时刻的尺寸如表1所示。

表1各运动副初始时刻尺寸(mm)

对各个随机参量进行参数化处理，建立试验机构的参数化多体动力学仿真模型，通过仿真试验得到四个运动副的磨损数据，结果如图4所示。

步骤2：分别建立每一个运动副的退化模型；

使用Gamma过程分别拟合四个运动副的磨损过程，其模型如下式：

退化增量表示为：

Δx_ik(t_j)的概率密度函数为：

步骤3：特征量k(1≤k≤4)对应的退化模型的参数为(μ_k,γ_k,λ_k)，采用MCMC方法进行退化参数估计，结果如表2所示。

表2各运动副磨损退化参数

步骤4：建立四元时变相关性模型：

步骤4(a)：随着运动周期的增加，各运动副磨损增加，出现损伤累积，由于各运动副受相同的连杆作用，承受相同的载荷，相同的工作环境，所以各运动副的损伤累积会互相影响，进而影响到机构的性能。使用D-VineCopula进行该试验机构的旋转副磨损量的相关建模，其模型如图5所示，四个运动副的联合概率密度函数为：

f(x_A,x_B1,x_B2,x_C)＝f_A(x_A)·f_B1(x_B1)·f_B2(x_B2)·f_C(x_C)·c_A,B1·c_B2,C·c_A,B2|B1·c_B1,C|B2·c_A,C|B1,B2 (8)

步骤4(b)采用Kendall系数的时变稳定模型

拟合运动副之间的时变相关性。

步骤5：确定步骤4中所有的二元Copula函数的类型和时变相关系数。

步骤5(a)：将Normal Copula、Gumbel、Frank、ClaytonCopula，作为Copula函数备选集；

步骤5(b)：将步骤3得到的特征量退化参数(μ_k,γ_k,λ_k)带入各边缘概率密度函数，四个旋转副的对数似然函数为：

通过

得到常相关系数估计值，从备选函数中选取AIC值最小的作为最优Copula，作为第一个二元Copula；

步骤5(c)：重复步骤5(a)和5(b)，直到确定步骤4中模型的所有二元Copula函数类型，结果如表3所示，六个copula函数的类型分别为：Normal、Clayton、Clayton、Clayton、Frank和Gumbel Copula。

表3基于D-VineCopula函数参数估计结果

步骤5(d)：采用最小二乘法估计时变模型中的参数，其结果如表4所示。

表4时变相关模型参数值

由Kendall系数τ与Copula参数θ之间的转换关系，如表5所示，将τ(t)转换为θ(t)，由此得到每一个Copula函数的时变相关系数值。

表5各类Copula函数的参数

其中，

步骤6：由表3和表4中Copula函数的类型及时变相关系数的模型，根据Rosenblatt变换，得到具有相关关系的样本。

步骤7：重复步骤6，获得m组相关样本。

步骤8：该试验机构主要的失效模式为位移滑块9向左运动的最大位移l达不到指标要求。最大位移受到各连杆长度和运动副间隙的影响，运动副磨损致使运动副间隙增大，导致最大位于逐渐减小，当l小于规定指标时，机构发生性能失效。

机构发生位移失效的极限状态方程为：g(t)＝l(t)-l₀。式中，l₀为位移额定值，l₀＝200mm。

最大位移l使用二次响应面法拟合得到，具体形式如式(14)，其结果如表6所示。

表6二次响应面拟合结果

步骤9：将抽样的样本带入极限状态方程，计算舱门锁试验机构的可靠度，其可靠度曲线如图6所示。

从可靠度曲线来看，试验机构的使用寿命主要集中在9000-9600次之间。由于舱门锁在收放次数达到9000次之后失效件数量急剧增加，其可靠性也开始急剧下降。

上述实施例仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (1)

1.一种机构可靠性分析方法，其特征在于，所述方法包括以下步骤：

步骤1：确定机构的n个退化特征量，通过退化试验收集试验数据y＝(y₁,y₂,...,y_n)；

步骤2：在不考虑各个特征量之间的相关性情况下，利用随机过程建立每一个特征量的退化模型；

步骤3：利用MCMC方法估计所述步骤2中退化模型的参数；

步骤4：建立多元相关模型和时变相关模型；

步骤5：对所述时变相关模型进行参数估计，选取二元Copula函数作为Copula函数备选集，确定多元相关Copula函数类型，得到每一个Copula函数的时变相关系数；

步骤6：根据步骤5得到的各Copula函数类型和时变相关系数，基于Rosenblatt变换方法，生成随机向量的m组样本；

步骤7：重复步骤6，对变量u₁,u₂,...,u_n进行m组随机抽样，得到m组随机向量X_m×n；

步骤8：建立失效模型的极限状态方程；

步骤9：将步骤7得到的样本带入步骤8的极限状态方程，使用蒙特卡洛法求解可靠度，落入失效域F内的样本点的个数N_um与总样本点的个数m之比即为失效概率的估计值P_f；

所述步骤2中所述随机过程为Wiener过程或Gamma过程；

所述步骤4具体包括：

步骤4(a)：建立多元相关模型：利用D-VineCopula模型建立n个特征量的相关性模型，将n元分布分解为多个二元相关的形式；

步骤4(b)：建立时变相关模型：使用基于Kendall系数τ的时变相关性模型描述特征量间的时变相关性；

所述步骤5具体包括：

步骤5(a)：选取二元Copula函数Normal Copula和/或Archimedean Copula,作为Copula函数备选集；

步骤5(b)：确定多元相关Copula函数类型：将步骤3得到的特征量退化参数带入步骤4(a)多元相关模型，用最大似然估计法得到每个备选Copula函数中的参数值，利用AIC准则从备选函数中选取最优Copula，作为第一个二元Copula；

步骤5(c)：重复步骤5(a)和5(b)，直到确定步骤4中模型的所有二元Copula函数类型；

步骤5(d)：使用最小二乘法求解步骤4(b)模型的参数；

步骤5(e)：根据Kendall系数τ与Copula参数θ之间的转换关系，根据步骤5(c)得到的Copula函数的类型，将τ(t)转换为θ(t)，由此得到每一个Copula函数的时变相关系数值；

所述步骤6包括：

步骤6(a)：令u₁,u₂,...,u_n分别表示n维连续随机变量X＝(X₁,X₂,...,X_n)的边缘分布函数，即u₁＝F₁(x₁),...,u_n＝F_n(x_n)，则x₁＝F₁ ^-1(u₁),...,x_n＝F_n ^-1(u_n)，设n维随机变量X的一个样本为(x₁,x₂,...,x_n)，r₁,r₂,...,r_n分别为[0,1]之间均匀分布的随机变量的一个样本，根据Rosenblatt变换，令u₁＝r₁，可得：x₁＝F₁ ^-1(u₁)；

步骤6(b)：由

得

则有

步骤6(c):由

得

则有x₃＝F₃ ^-1(u₃)；

步骤6(d)：由＝h_4|23,1|23(h_4|3,2|3(h₄₃(u₄,u₃),h₂₃(u₂,u₃)),h_1|2,3|2(h₁₂(u₁,u₂),h₃₂(u₃,u₂)))＝r₃，可得

则有

步骤6(e)：以此类推得到

由此得到随机向量X的一组样本(x₁,x₂,...,x_n)；

其中，F(·)为随机变量的边缘分布函数，F^-1(·)为边缘分布函数的逆函数，F_i|j(x_i|x_j)为随机变量x_j条件下随机变量x_i的条件分布函数；h_ij(i,j)为条件Copula函数。

Images