说明书一种动态级联耦合机电系统的可靠性建模与评估方法
技术领域
本发明属于可靠性技术领域，具体涉及一种动态级联耦合机电系统的可靠 性建模与评估方法。
背景技术
复杂机电系统大多在复杂服役环境下运行，系统行为与物理背景密切相关， 系统服役过程中性能逐渐劣化，零部件在环境载荷持续作用下损伤不断发展， 呈现时变和多态特点；系统内部层次复杂，各层之间存在上下级联、层内各单 元和各因素之间存在相互耦合的特点；系统故障呈现出多故障机理、多种故障 模式相互作用和竞争，且相互作用及竞争关系呈现动态变化的特点。复杂机电 系统的动态、级联、耦合、时变特点决定了传统可靠性理论的故障独立性、大 样本统计和静态假设条件不再成立。自然科学基金委员会2010年出版的《机械 工程学科发展战略报告》也明确指出：传统的系统可靠性方法不适用。
针对复杂机电系统动态可靠性建模方法，目前主要有：神经网络建模法， 马尔科夫分析法、GO-Flow法、Petri网、贝叶斯网等；其中，Petri网和贝叶斯 网是主流的两类方法。Petri网模型主要描述和分析离散事件动态系统的图形工 具和信息流的网络，能够描述系统动态特性，目前已在复杂机电系统可靠性方 面得以应用。贝叶斯网络方法(Bayesian networks，BN)能很好地表示变量的随机 不确定性和多种故障模式相关性，使系统和元件的关系以及状态表达更加直观、 清晰，已在可靠性分析评估、故障诊断、数据挖掘和经济等领域得到了应用。 然而，上述方法均基于事件发生概率统计思想，确定性信息利用不充分，与物 理背景结合不紧密。目前复杂机电系统的可靠性建模分析方法研究尚缺乏成熟 的理论体系支撑，更为严重的是复杂机电系统的某些功能突变、性能劣化、故 障演变与涌现等问题，难以在可靠性建模与分析过程中进行准确的分析表达。
发明内容
有鉴于此，本发明提供了一种动态级联耦合机电系统的可靠性建模与评估 方法，能够对动态级联耦合的复杂机电系统进行可靠性建模并对可靠度进行评 估。
本发明的一种动态级联耦合机电系统的可靠性建模与评估方法，包括如下 步骤：
步骤1、将机电系统自上而下分解成多个层次，每个层次包括若干不同单元， 由此建立机电系统的结构层次图；根据所述结构层次图，将机电系统的功能进 行划分，进而得到机电系统的功能层次图，其中，所述结构层次图中的一个单 元对应功能层次图中的一个功能；
步骤2、根据步骤1建立的机电系统结构层次图和功能层次图，参照 GJB1391-2006对机电系统进行故障模式及影响分析，根据该分析结果：
首先得到所述每个单元对应的故障模式以及每个故障模式对应的故障类 别，其中，故障模式由单元的性能响应和对应的性能阈值之间的关系表征；故 障类别包括退化型故障和突发型故障；
然后从所述结构层次图中由下至上逐层确定每个单元中的性能响应对本层 次各单元和其它层次各单元的性能响应的影响，即得到机电系统的故障传播结 果；
步骤3、根据步骤2的故障模式及影响分析结果，确定每个所述单元故障模 式对应的故障原因，进而确定每个单元的可靠性影响因素及其作用机理；
步骤4、建立各层次单元的性能响应模型，结合步骤3确定的可靠性影响因 素确定模型中的不确定性参数及其概率分布特征，具体为：
S4.1：建立各层次单元的性能响应模型
根据步骤1的所述结构层次图，以及步骤3所得的每个单元的可靠性影响 因素及其作用机理，并结合动力学方程，由下至上逐层建立各单元的性能响应 模型：
M q · · + M q · + Kq + φ q T λ = F 1 + F 2 - - - ( 1 ) ]]>
式中，M、C和K分别代表单元的质量矩阵、阻尼矩阵和刚度矩阵；和 q代表单元的运动状态参数；代表运动约束条件，其中上角标T表示矩阵的转 置；λ为拉格朗日系数；F1代表上、下层或本层的其他单元对本单元施加的外 力；F2代表本单元内静止部件对活动部件施加的作用力；
S4.2：建立各层次单元的响应状态模型
由S4.1中的公式(1)解得各单元的性能响应g(X)，其中，X＝(X1,X2,…，Xn)表 示不同层次单元的性能响应模型中的相关参数，n表示所述相关参数的个数；则 各层次单元的响应状态模型为：
Z＝b-g(X)＝b-g(X1,X2,…,Xn) (2)
其中，b为性能响应对应的性能阈值；
S4.3：将每个单元性能响应模型中的相关参数X与机电系统常用的不确定性 参数列表进行比对，确定每个所述单元的性能响应中的不确定性参数 XU＝(x1,x2,…,xk)以及不确定性参数的概率分布类型；其中，k为各性能响应中不 确定参数个数；
步骤5、建立机电系统可靠性网络模型，具体为：
S5.1：将S4.3中的不确定性参数作为底层节点，将所有单元的性能响应作 为上层节点，建立机电系统可靠性网络模型结构，其中，各性能响应在网络模 型结构中的层次位置与其所属的单元所在所述结构层次图中的位置对应；
S5.2：根据S4.2和S4.3确定的不确定参数与性能响应的对应关系，以及步 骤2得到的机电系统故障传播结果，在机电系统可靠性网络模型结构中由下至 上采用有向弧连接各层节点，建立机电系统可靠性网络模型；有向弧首、尾连 接的分别为父节点和子节点，其中，父节点为施加影响一方，子节点为接受影 响一方；
步骤6、对不确定信息从机电系统可靠性网络模型的底层到顶层的耦合与传 播进行递推，最终得到机电系统的可靠度，具体为：
S6.1：根据S4.3中得到的不确定参数及其概率分布类型，得到每个不确定 参数的概率密度函数
S6.2：首先，根据底层不确定参数的概率密度函数，采用公式(4)得到机 电系统可靠性网络模型底层的上一层中各性能响应的概率密度函数；
然后，利用低层次单元的性能响应概率密度函数，采用公式(4)对高层次 单元的性能响应概率密度函数进行逐层递推，直至得到顶层各性能响应的概率 密度函数；
f ( g c ) = Π j = 1 K f ( Y j c / parents ( Y j c ) ) - - - ( 4 ) ]]>
其中，f(gc)表示节点gc的概率密度函数，Π(·)表示求积符号，Ycj为与节点 gc通过有向弧连接的节点，j＝1,2,...,K，K为与节点gc通过有向弧连接的节点总 数，parents(Yjc)为与节点Ycj通过有向弧连接的节点，f(Yjc/parents(Yjc))表示节点 Ycj在其他节点为parents(Yjc)时的条件概率密度函数；
S6.3：根据下式得到网络模型顶层各性能响应的故障发生概率：
P i = P ( b i - g i ≤ 0 ) = 1 - P ( g i < b i ) = 1 - ∫ - ∞ ∞ f g i ( g i ) ∫ b i + ∞ f b i ( b i ) d g i d b i - - - ( 5 ) ]]>
其中，表示机电系统可靠性网络模型顶层中第i个性能响应的概率密 度函数，bi表示机电系统可靠性网络模型顶层中第i个性能响应对应的性能阈值， 表示网络模型顶层中第i个性能响应对应的性能阈值概率密度函数；
则网络模型顶层中第i个性能响应的可靠度为Ri＝1-Pi；
机电系统的可靠度为其中，h为网络模型顶层中性能响应的个 数；
利用机电系统的可靠度对机电系统的可靠性进行评估。
进一步的，利用不断获得的机电系统的试验数据，根据贝叶斯定理对所述 网络模型中的除底层节点之外的性能响应的概率密度函数进行更新：
f θ ( θ | Y * ) = f θ ( θ ) f θ ( Y * | θ ) ∫ f θ ( θ ) f θ ( Y * | θ ) dθ - - - ( 6 ) ]]>
其中，θ为需要更新的性能响应，Y*为试验数据中性能响应θ的值；fθ(θ)为 性能响θ的概率密度函数，fθ(Y*|θ)为Y*在θ条件下的概率密度函数，fθ(θ|Y*)为θ 在Y*条件下的条件概率密度函数；进而根据公式(5)得到各性能响应的故障发 生概率，最终得到更新后的机电系统的可靠度。
进一步的，所述步骤6的S6.3中，当与退化型故障相关的性能响应对应的 性能阈值为随时间退化的性能退化阈值bt时，网络模型的顶层中与退化型故障相 关的性能响应在t时刻的退化型故障发生概率为：
表示机电系统可靠性网络模型顶层中第i个性能响应在t时刻的概率 密度函数，表示网络模型顶层中第i个性能响应对应的性能阈值在t时刻 的概率密度函数；
当退化型故障中性能响应对突发型故障产生影响时，令机电系统可靠性网 络模型顶层中第s个性能响应t时刻的突发型故障的阈值为ls，其概率密度函数 为，与突发型故障相关的性能响应为yst，则突发型故障的发生概率为 退化型故障未发生、突发型故障发生的条件概率，即
则机电系统的可靠度为：
其中，max{·} 表示取最大值符号；I表示网络顶层中与退化型故障有关的性能响应个数，S表 示网络顶层中与突发型故障有关的性能响应个数。
进一步的，当机电系统先后执行多个任务时，机电系统的可靠度为： 其中，Ri′(t)为执行第i′个任务时的机电系统可靠度，根据执行该 任务时验数据并采用所述步骤1至6获得；I′表示任务总数。
本发明具有如下有益效果：
1)、本发明提供了一种动态级联耦合机电系统的可靠性建模与评估方法， 基于复杂机电系统的物理背景，以动力学响应模型为基础进行机电系统可靠性 的研究，所建立的机电系统可靠性网络模型更能准确描述复杂机电系统的动态、 级联和耦合的特点，使得本发明基于该网络模型获得的机电系统可靠度更加准 确可信；同时，利用动力学分析和不确定性分析结合的方法，描述不确定信息 从网络模型底层到顶层的传递，解决了由于实验费用大、周期长、技术难度大 及产品造价昂贵而导致航天机构系统响应不确定性特征、系统可靠性信息无法 通过实验手段得到的难题。
2)、本发明通过不断更新的机电系统试验数据对网络模型中的性能响应的 概率密度函数进行更新，从而得到机电系统各性能响应更加准确的概率密度函 数。
3)、本发明通过考虑性能阈值的时间退化特性,建立广义响应-阈值干涉模 型，得到实时的退化型故障的发生概率，由此可反应机电系统可靠性的时变特 性；同时考虑退化型故障对突发故障的影响，使得本发明对故障类型的考虑更 全面，更加符合机电系统的实际工作情况，最终得到更准确的系统可靠度。
附图说明
图1为本发明中复杂机电系统结构层次图；
图2为本发明中机电系统的结构和功能关系图；
图3(a)为本发明中退化型故障性能-时间曲线图，图3(b)为突发型故障 功能-时间曲线图；
图4(a)为本发明中网络模型的串连节点子拓扑结构，图4(b)为本发明 中网络模型的散连节点子拓扑结构，图4(c)为本发明中网络模型的会聚节点 子拓扑结构；
图5为本发明中网络模型节点的条件概率示意图；
图6为本发明中性能响应—阈值干涉模型；
图7为本发明中系统多任务可靠性模型；
图8为本发明实施例中双轴驱动机构层次结构划分模型；
图9为本发明实施例中双轴驱动机构功能分析；
图10为本发明实施例中输入轴输出转速概率密度函数；
图11为本发明实施例中柔轮输出转速概率密度函数；
图12为本发明实施例中齿轮啮合力概率密度函数；
图13为本发明实施例中横轴输出转速概率密度函数；
图14为本发明实施例中纵轴输出转速概率密度函数；
图15为本发明实施例中双轴驱动机构输出转速概率密度函数；
图16为本发明实施例中间隙随时间的变化趋势；
图17为本发明实施例中双轴驱动机构可靠度随时间的变化趋势；
图18为本发明实施例中间隙随时间变化趋势；
图19为本发明实施例中双轴驱动机构可靠度随时间变化趋势。
具体实施方式
下面结合附图并举实施例，对本发明进行详细描述。
本发明提供了一种动态级联耦合机电系统的可靠性建模与评估方法，包括 如下步骤：
步骤1：建立机电系统结构层次图，并进行系统功能分析
步骤1.1：根据机电系统的工作原理、任务剖面，设计手册、可靠性要求等 相关资料，将机电系统自顶向下分解成多个层次，建立机电系统结构层次图， 例如可分成系统层、子系统层和部件层三个层次，每个层次由不同单元组成， 如图1所示，子系统层中各单元为子系统1/2/3，部件层单元为各部件。
步骤1.2：根据机电系统结构层次图，进行机电系统分层次功能分析，使得 结构层次图中的一个单元对应功能层次图中的一个功能得到系统功能层次图。 如图2所示，图中各功能分别对应于图1中各层单元，如子功能1对应图1中 子系统层单元子系统1。
步骤2：进行机电系统故障行为分析
步骤2.1：进行故障模式及影响分析(FMEA)：
根据步骤1给出的机电系统结构层次图和功能层次图，参照GJB1391-2006 进行系统FMEA分析。首先得到每个单元对应的故障模式以及每个故障模式对 应的故障类别，其中，故障模式由单元的性能响应和对应的性能阈值之间的关 系表征；通过FMEA分析得到的故障模式将系统的故障进行分类——突发型故 障和退化型故障。其中，与系统性能指标阈值相关的故障为退化型故障，与系 统功能相关的故障为突发型故障(例如由于外部冲击或过应力引起的系统功能 故障)。突发和退化两类故障的示意图如图3所示。
步骤2.2：根据FMEA输出结果，得到每个单元对应的故障模式以及每个故 障模式对应的故障类别，其中，故障模式由单元的性能响应和对应的性能阈值 之间的关系表征。然后进行系统故障传播的定性分析，从底向上逐层分析当前 层次各单元对上层次单元的影响、对本层次单元的影响以及对系统顶层的最终 影响。故障传播方式可分为与退化型故障判据相关的性能退化故障传播，以及 与突发型故障判据相关的突发故障传播方式，并且这两种故障传播方式对于系 统整体状态的影响，还存在着竞争效应。最后，填写故障传播定性分析表，如 表1所示。
表1系统故障传播定性分析表
机电系统故障行为分析结果，可为后续进系统可靠性模型建立奠定基础。
步骤3：根据系统故障模式和故障原因，分析系统的可靠性影响因素及其作 用机理：
根据步骤2得到的故障模式及故障原因，通过定性或试验手段分析影响系 统可靠性的影响因素及其作用机理。例如在典型的复杂机电产品中，完成运动 功能和任务的系统各组成部件之间存在设计的原始间隙，以及制造和装配过程 中产生的误差，系统运行期间发生磨损等因素，使得部件间的间隙逐渐变大， 影响机构系统的运动功能和运动精度，导致系统可靠度下降，最终影响任务完 成。则设计、制造、装备和使用中的部件之间的间隙，即为复杂机电系统的可 靠性影响因素；而由于磨损变形造成的间隙逐渐变大进而影响系统可靠性，即 为可靠性影响因素的作用机理。可靠性影响因素及其作用机理将被用到后续建 立性能响应模型过程中。
步骤4：建立各层次单元的性能响应模型，结合步骤3确定的可靠性影响因 素确定模型中的不确定性参数及其概率分布特征：
S4.1：建立各层次单元的性能响应模型
复杂机电单元的动力学响应是影响系统性能和可靠性的主要因素。因此， 本发明将利用动力学方程建立性能响应模型来描述系统性能响应变化。根据步 骤1的系统结构层次图，以及步骤3所得的可靠性影响因素及其作用机理，从 底层到顶层，结合动力学方程，逐层建立各单元的性能响应模型，如式(1)所 示:
M q · · + M q · + Kq + φ q T λ = F 1 + F 2 - - - ( 1 ) ]]>
式中，M、C、K分别代表单元的质量矩阵、阻尼矩阵和刚度矩阵，和 q分别代表单元的运动状态参数；代表运动约束条件，其中T表示矩阵的转置； λ为拉格朗日系数；F1代表上、下层或本层的其他单元对本单元施加的外力；F2代表本单元内静止部件对活动部件施加的作用力；如步骤3磨损作用机理引起 的轴承对轴接触碰撞力。通过求解方程(1)可得出系统不同层次单元的性能响 应，如位移、速度等(即和q)，为步骤4.2提供性能响应关系。
S4.2：建立各层次单元的响应状态模型
由S4.1中的公式(1)解得各单元的性能响应g(X)(如位移、角速度和扭 矩等)，其中，X＝(X1,X2,…，Xn)表示不同层次单元的性能响应模型中的相关参数， 如结构尺寸、密度、材料属性、阻尼系数、接触力等,则各层次单元的响应状态 模型为：
Z＝b-g(X)＝b-g(X1,X2,…,Xn) (2)
其中，b为性能响应对应的性能阈值；当Z＞0时，单元为正常状态，当Z≤0 时，单元为故障状态。
S4.3：将每个单元性能响应模型中的相关参数X与机电系统常用的不确定性 参数列表进行比对，确定每个所述单元的性能响应中的不确定性参数 XU＝(x1,x2,…,xk),其中，k为各性能响应中不确定参数个数；并根据试验数据、历 史经验和相似单元等信息获得不确定性参数的概率分布类型。常用的不确定性 参数的概率分布类型如表2所示。
表2常用的不确定性参数的概率分布类型
步骤5：建立机电系统可靠性网络模型
本步骤集成系统结构层次图、功能层次图、FMEA分析结果、性能响应模 型，结合贝叶斯网理论，建立系统可靠性网络模型，描述系统状态演化过程。
步骤5.1：选取和定义系统可靠性网络模型节点；
S5.1：将步骤4.3中的不确定性参数作为机电系统可靠性网络模型底层节点， 将所有单元的性能响应作为机电系统可靠性网络模型的上层节点，建立机电系 统可靠性网络模型结构，其中，各性能响应在网络模型中的层次位置与其所属 的单元所在结构层次图中的位置对应；
S5.2：根据S4.2和S4.3确定的不确定参数与性能响应的对应关系，以及步 骤2得到的机电系统故障传播结果，在机电系统可靠性网络模型中由下到上采 用有向弧连接各层节点，建立一个有向无环图；有向弧首尾连接的分别为父节 点和子节点，表明了父节点和子节点之间的关系，即故障影响定性关系，其中， 父节点为施加影响一方，子节点为接受影响一方。如图4所示，所建的网络模 型共包含三种子拓扑结构，分别为串连节点、散连节点和会聚节点子拓扑结构。
S5.3：构造每个节点的条件概率：
用条件概率表来描述子系统节点与父节点间的逻辑关系。如图5所示，父 节点A的性能响应概率为P(A)，其子节点B和C的性能响应概率为条件概率 P(B,C|A)，即在A发生的情况下B和C发生的概率。
步骤6：对不确定信息从网络模型的底层到顶层的耦合与传播进行递推，最 终得到机电系统的可靠度，具体为：
S6.1：根据S4.3中得到不确定参数及其概率分布类型，结合试验数据、经 验数据和相似单元等信息，对机电系统可靠性网络模型底层节点进行不确定性 特征分析，得到每个不确定参数的概率密度函数
S6.2：首先，根据底层不确定参数的概率密度函数，采用公式(4)得到底 层的上一层次中各性能响应的概率密度函数；本实施例中即为子系统层中各性 能响应的概率密度函数；
然后，利用低层次单元的性能响应概率密度函数、性能响应模型，以及仿 真数据、试验数据，采用公式(4)对高层次单元的性能响应概率密度函数进行 逐层递推，直至得到顶层各性能响应的概率密度函数；
f ( g c ) = Π j = 1 K f ( Y j c / parents ( Y j c ) ) - - - ( 4 ) ]]>
其中，公式(4)基于S5.3中节点的条件概率得到，式中，f(gi)表示节点gi的概率密度函数，Π(·)表示求积符号，Yij为与节点gi通过有向弧连接的节点， j＝1,2,...,K，K为与节点gi通过有向弧连接的节点的总数，parents(Yij)为与节点Yij通过有向弧连接的节点，f(Yij/parents(Yij))表示节点Yij在其他节点为parents(Yij)时 的条件概率密度函数，通过核密度估计方法得到。
S6.3：根据下式得到网络模型顶层各性能响应的故障发生概率：
P i = P ( b i - g i ≤ 0 ) = 1 - P ( g i < b i ) = 1 - ∫ - ∞ ∞ f g i ( g i ) ∫ b i + ∞ f b i ( b i ) d g i d b i - - - ( 5 ) ]]>
其中，表示网络模型顶层中第i个性能响应概率密度函数，表示 网络模型顶层中第i个性能响应对应的阈值概率密度函数；
则网络模型顶层中第i个性能响应的可靠度为Ri＝1-Pi；机电系统的可靠度 为其中，h为网络模型顶层中性能响应的个数。利用机电系统的可 靠度对机电系统的可靠性进行评估。
进一步的，在步骤6中，利用不断获得的机电系统的试验数据，根据贝叶 斯定理(式(6))采用马尔科夫链蒙特卡罗方法(Markov Chain Monte Carlo， MCMC)对网络模型中的除底层节点之外的性能响应的概率密度函数进行更新：
f θ ( θ | Y ) = f θ ( θ ) f θ ( Y | θ ) ∫ f θ ( θ ) f θ ( Y | θ ) dθ - - - ( 6 ) ]]>
其中，θ为需要更新的性能响应，Y为试验数据中性能响应θ的值；fθ(θ)为 性能响θ的概率密度函数，fθ(Y|θ)为Y在θ条件下的概率密度函数，fθ(θ|Y)为θ在 Y条件下的条件概率密度函数；进而根据公式(5)得到更加准确的各性能响应 的故障发生概率，最终得到更新后的机电系统的可靠度。
由上述方案可知，机电系统的故障通常是由系统本身的作用机理和外部工 作环境共同作用的结果，该过程是一个极其复杂的过程。根据步骤2可知，主 要是退化故障和突发故障两类故障，故可以认为复杂机电系统的故障是这两类 故障模式不断竞争的结果，且故障是由最先发生的故障模式所导致。由于故障 过程复杂，为了避免常规分析方法对可靠性分析造成的误差，基于上述可靠性 网络模型提出更为有效的可靠性分析方法。下面根据假设给出竞争失效的一般 模型：
在步骤6的S6.3中，当与退化型故障相关的性能响应对应的性能阈值为随 时间退化的性能退化阈值bt时，结合时刻t的性能阈值或功能阈值及公式(2)， 建立广义响应—阈值干涉模型(如图6所示)，得到网络模型的顶层中与退化型 故障相关的性能响应在t时刻的退化型故障发生概率为：
表示网络模型顶层中第i个性能响应在t时刻的概率密度函数，根据 仿真数据和试验数据得到；表示网络模型顶层中第i个性能响应对应的性 能阈值在t时刻的概率密度函数，通过统计方法得到。
当退化型故障中性能响应对突发型故障产生影响时，令t时刻突发故障的阈 值为l，其概率密度函数为fl(l)(通过统计方法给出)，与突发型故障相关的性能 响应为yt，则突发型故障的发生概率Pt突发为退化型故障未发生、突发故障发生 的条件概率，即
则机电系统的可靠度为：
其中，I表 示网络顶层中与退化型故障有关的性能响应个数，S表示网络顶层中与突发型故 障有关的性能响应个数。
进一步的，当机电系统执行多个任务时，建立系统多任务可靠性模型，为 串联逻辑模型，如图7所示；根据上述步骤1至6的计算过程，结合各任务下 的试验数据，分别得到单个任务下的机电系统可靠度Ri′(t)，按照式(9)计算多 任务下的系统可靠度。
R m ( t ) = Π i ′ = 1 I ′ R i ′ ( t ) - - - ( 9 ) ]]>
实施例：
任务描述：双轴驱动机构动态级联耦合可靠性问题
卫星天线双轴驱动机构是空间机构中典型的机电系统，其可靠性水平直接 关系到卫星天线的安全和任务完成的可靠度。双轴驱动机构是天线双轴运动的 核心部件，包括横轴驱动组件、纵轴驱动组件，两个驱动组件完全一致。
双轴驱动机构“动态级联耦合可靠性”问题中：各部件、子系统和系统动力学 特性(力矩和转速)、运动副间隙随时间发生变化；系统响应(天线的位移和角 度)是由部组件(输入轴、刚柔轮、输出轴)、不同层次的子系统(纵轴和横轴) 的动力学响应构成；确定性和不确定性信息在不同层次间存在传递现象；考虑 机构运动副间隙影响，刚柔轮、轴与轴承存在动力学刚柔耦合问题。下面利用 本发明针对其动态级联耦合可靠性问题，分析时变间隙下，双轴驱动机构可靠 性动态变化规律。
具体步骤如下：
步骤1：建立双轴驱动机构结构层次图，并进行系统功能分析
步骤1.1：根据双轴驱动机构结构组成，将系统分为系统层、子系统层、和 部件层三层，如图8所示。
系统层单元为双轴驱动机构产品系统。
子系统层单元分别为横轴驱动组件、纵轴驱动组件和天线组件。
部件层单元分别为输入轴系1、刚柔轮运动副1、输出轴系1、输入轴系2、 刚柔轮运动副2以及输出轴系2。
步骤1.2：根据系统结构层次图，进行系统分层次功能分析，得到系统功能 层次图，如图9所示。
系统层功能为：实现跟踪和指向任务，从而完成卫星之间通信链路的建立 和连续通信。
子系统层各单元功能依次为：横轴组件输出力矩和转速带动纵轴组件级天 线组件转动；纵轴组件输出力矩和转速带动天线组件转动；天线组件在横轴组 件和纵轴组件的耦合运动下实现自身的平移和转动。
部件层各单元功能依次为：提供力矩和转速，传递力矩和转速，输出力矩 和转速，提供力矩和转速，传递力矩和转速以及输出力矩和转速。
步骤2：进行双轴驱动机构故障行为分析
步骤2.1：进行故障模式及影响分析(FMEA)进行。
通过FMEA分析，可知双轴驱动机构的关键件为：驱动组件中谐波减速器 柔轮、输入轴轴承、输出轴轴承。关键件的故障模式主要为磨损变形。
步骤2.2：根据FMEA输出结果，进行系统故障传播的定性分析，从底向上 逐层分析当前层单元对上层次单元的影响、对本层次单元的影响以及对系统顶 层的最终影响。对本层次单元的影响为横向传播，对上层次单元的影响为纵向 传播。纵向的传播关系又可分为与退化型故障判据相关的性能退化故障传播， 以及与突发型故障判据相关的突发故障传播方式，并且这两种故障传播方式对 于系统整体状态的影响，还存在着竞争效应。最后，填写故障传播定性分析表， 如表3所示。
表3系统故障传播定性分析表
为后续进系统可靠性模型建立奠定基础。
步骤3：根据系统故障模式和故障原因，分析系统的可靠性影响因素及其作 用机理
根据步骤2得到的故障模式及故障原因，确定可靠性主要影响参数为：电 机转速、柔轮弹性模量、柔轮齿轮厚度、刚轮刚度系数、齿轮直径、输入轴半 径。输入轴输出力矩和转速；柔轮输出力矩和转速；齿轮接触力。横轴输出力 矩、转速；纵轴输出力矩、转速。
可靠性影响因素的作用机理为由于磨损变形造成间隙变大导致运动副间产 生碰撞力，从而影响系统响应的稳定性，导致机构末端天线运动存在动态误差 及抖动。本案例中的故障机理为磨粒磨损，选用Archard磨损模型作为作用机理 模型。Archard磨损模型是工程中广泛接受的磨损形式之一，其形式如下所示：
V = kPL H - - - ( 10 ) ]]>
式中，V是磨损体积；P为接触载荷；L为滑动距离；H为材料硬度；k为 磨损系数。本案例接触载荷P利用动力学仿真手段得到。
在这个模型中，当配对材料和工况条件确定以后，磨损系数是决定磨损量 的关键。但磨损系数无理论值，只能通过实验或经验确定。对于金属滑动摩擦 副发生粘着磨损时，磨损系数可在10-2–10-10之间变化。Archard模型计算的磨 损量可以认为是一次相对运动产生，那么可以得到累积磨损量，如下所示：
V ( t ) = kPL H t T - - - ( 11 ) ]]>
式中，T为运动周期，t为部件运行时间。由于制造误差、使用环境等因素 影响，除T和t之外，其他参数均可看作随机变量。
步骤4：逐层建立各层级单元的性能响应模型，确定模型中的不确定性参数 及其分布特征
步骤4.1：建立各层次单元的性能响应模型
从底到上，考虑影响动力学响应的各个因素，包括构件间隙、材料属性、 装配精度、预紧力等，结合步骤3给出的可靠性关键影响因素，逐层建立双轴 驱动部件、子系统和系统的动力学响应模型，并确定组件与子系统，子系统与 系统之间的动力学耦合关系。
1)部件层：
(1)输入轴性能响应模型
M—轴的质量；J—轴的转动惯量矩阵；C—轴的阻尼系数矩阵；K—轴的 刚度系数矩阵；x、y、z—轴的平动位移；ψ、θ、—轴的转角位移；FD—驱动载 荷；FP—负载载荷；FBI—输入轴轴承载荷；FBF—柔性轴承载荷；TD—驱动转矩； TP—负载力矩；Tf—轴承的滑动摩擦力矩；Ty—轴承的碰撞力等效在y方向的转 矩；Tz—轴承的碰撞力等效在z方向的转矩。
J = J xx J xy J xz J yx J yy J yz J zx J zy J zz . ]]>[FPx FPy FPz TPx TPy TPz]′为耦合负载载荷。考虑本实施 例中的动力学方程高非线性无法直接求解，故利用动力学仿真模型进行仿真分 析得到，分析得到输入轴的转速和力矩响应结果。
(2)刚柔轮性能响应模型为
T＝(θgRg-(θg+θr)Rr)k (13)
其中，凸轮转动的角度为θg，扫过刚轮的弧长为θgRg。设柔轮相对刚轮的转 角为θr，则凸轮扫过柔轮的弧长为(θg+θr)Rr。Rg—刚轮分度圆半径；Rr—柔轮 分度圆半径；k—啮合刚度；考虑本实施例中的动力学方程高非线性无法直接求 解，故利用动力学仿真模型进行仿真分析得到，分析得到柔轮的转速和力矩响 应结果。
(3)输出轴性能响应模型
M—轴的质量；J—轴的转动惯量矩阵；C—轴的阻尼系数矩阵；K—轴的 刚度系数矩阵；x、y、z—轴的平动位移；ψ、θ、—轴的转角位移；FD—驱动载 荷；FP—负载载荷；FBI—输入轴轴承载荷；FBF—柔性轴承载荷；TD—驱动转矩； TP—负载力矩；Tf—轴承的滑动摩擦力矩；Ty—轴承的碰撞力等效在y方向的转 矩；Tz—轴承的碰撞力等效在z方向的转矩。
J = J xx J xy J xz J yx J yy J yz J zx J zy J zz ; C = C xx C xy C xz C yx C yy C yz C zx C zy C zz ; K = K xx K xy K xz K yx K yy K yz K zx K zy K zz ; ]]>
［FDx FDy FDz TDx TDy TDz］′为耦合驱动载荷；［FPx FPy FPz TPx TPy TPz］′为耦合 负载载荷。考虑本实施例中的动力学方程高非线性无法直接求解，故利用动力 学仿真模型进行仿真分析得到，分析得到输出轴的转速和力矩响应结果。
2)子系统层
横轴/纵轴性能响应模型
由于横轴/纵轴均由输入轴、刚柔轮和输出轴组成，故其性能响应模型由公 式(11)和公式(13)组合而成，其中两个公式中相同参数作为横轴/纵轴性能 响应模型的耦合参数。考虑本实施例中的子系统层耦合关系下动力学方程高非 线性无法直接求解，故利用动力学仿真模型进行仿真分析得到，分析得到纵轴/ 横轴的转速和力矩响应结果。
3)系统层
双轴驱动机构性能响应模型写成矩阵形式如下：
H b H bm H bm T H m x · · b θ · · + c b c m = F b τ m + J b T J m T F h - - - ( 15 ) ]]>
式中，基体惯量矩阵Hb
H b ∈ R 6 × 6 ≡ H v H vw H vw T H w - - - ( 16 ) ]]>
Hv∈R3×3≡wE3 (17)
H vw ∈ R 3 × 3 ≡ w r ~ 0 g T - - - ( 18 ) ]]>
H w ∈ R 3 × 3 ≡ Σ i = 1 n ( I i + m i r ~ 0 i T r ~ 0 i ) + I 0 - - - ( 19 ) ]]>
w = Σ i = 0 n ( m i ) - - - ( 20 ) ]]>
基体与双轴系统耦合惯量矩阵Hbm
H bm ∈ R 6 × n ≡ H vθ H wθ - - - ( 21 ) ]]>
H vθ ∈ R 3 × n ≡ Σ i = 1 n m i J Ti - - - ( 22 ) ]]>
H wθ ∈ R 3 × n ≡ Σ i = 1 n ( I i J Ri + m i r ~ 0 i J Ti ) - - - ( 23 ) ]]>
双轴系统惯量矩阵Hm
H m ∈ R n × n ≡ Σ i = 1 n ( J Ri R I i J Ri + m i J Ti T J Ti ) - - - ( 24 ) ]]>
雅克比矩阵
JTi∈R3×n≡[k1×(ri-p1),k2×(ri-p2),…，ki×(ri-pi),0,…,0] (25)
JRi∈R3×n≡[k1,k2,…,ki,0,…,0] (26)
J b ∈ R 6 × 6 ≡ E 3 p ~ 0 e T 0 E 3 - - - ( 27 ) ]]>
J m ∈ R 6 × n ≡ J Te J Re - - - ( 28 ) ]]>
JTe∈R3×n≡[k1×(pe-p1),…,kn×(pe-pn)] (29)
JRe∈R3×n≡[k1,k2,…,kn] (30)
r0g∈R3≡rg-r0 (31)
r0i∈R3≡ri-r0 (32)
p0e∈R3≡pe-r0 (33)
其中，Ii为轴i相对质心的惯量阵；ri为杆件i质心的位置矢量；ωi为角速度， vi为杆i质心的速度。xb∈R6为基体位置和方向， x · b = v 0 T ω 0 T T ; ]]>θ∈Rn＝[θ1 θ2 … θn]为关节角度变量；E3为3×3的单位矩阵；为矢量r的叉乘 矩阵，即满足JTi和JRi分别为平动和转动的雅克比矩阵。pi为 关节i的位置矢量；ri为杆件i质心的位置矢量；ki为关节i的转动轴方向单位向量。 Fb为基体所受外力及外力矩，τm为驱动关节力矩， F h = F e T τ e T T ]]>为执行末端所 受外力/外力矩。 c b = c v T c ω T T ]]>和cm分别为基体和双轴系统依赖速度的非线性项。 考虑本实施例中的系统层动力学方程高非线性无法直接求解，故利用动力学仿 真模型进行仿真分析得到，得到双轴驱动机构的转速和力矩响应结果。
将双轴驱动机构中的相关参数与表2所示的机电系统常用不确定参数列表 进行比对，发现所有相关参数均为不确定性参数，同时获得各不确定参数的概 率分布类型。
步骤5：建立系统可靠性网络模型
步骤5.1：结合步骤3，选取和定义系统可靠性网络模型节点；本案例可靠 性网络模型中，将节点进一步细分：圆圈节点具有随机特性，矩形节点表示试 验数据，各节点定义如下：
横轴/纵轴柔轮参数：θ1&2f＝[E,H]，其中E弹性模量和H齿轮厚度；
横轴/纵轴刚轮参数：θ1&2r＝[K,d]，K—刚度系数和d齿轮直径；
横轴/纵轴输入轴参数：θ1&2w＝[w,r1&2]，转速w及输入轴半径r1&2；
横轴/纵轴刚轮柔轮响应输出转速：由基本参数K,d，E,H作为 输入的响应函数模型。由有限元和动力学仿真得到。
横轴/纵轴输入轴响应输出转速：由输入轴参数w,r1&2组成的转 速的响应函数。由有限元和动力学仿真得到。
Xs1&2-横轴/纵轴输出轴的转速和力矩响应函数；
XA-双轴驱动机构输出的转速和力矩函数；
-试验数据：横轴/纵轴刚柔轮的性能试验数据—转速和力矩；
-试验数据：横轴/纵轴输出轴的性能试验数据—转速和力矩；
εf,εr,εw-通过仿真计算和试验测试结果得到的误差；
步骤5.2：定义系统可靠性网络模型拓扑结构；
根据步骤2的故障传播分析结果(见表3)，首先输入轴转速的波动会影响 柔轮的转速波动；其次，柔轮的转速波动及刚轮的不确定因素共同作用会影响 输出轴的转速波动；最终，横轴和纵轴输出轴的转速波动在耦合关系下影响系 统的转速波动从下到上用有向弧连接各层节点，建立一个有向无环图，如图10 所示。
步骤5.3：结合步骤5.2的可靠性网络模型拓扑结构，构造系统各节点的条 件概率关系如图11所示：
步骤6：进行不确定性推理，分析不确定信息从模型底层到顶层的耦合与传 播
步骤6.1：首先给出步骤5.1中各基本参数的概率分布如下表所示，不失一 般性，本发明假设各参数均服从正态分布，其均值和变异系数如下表4所示：
表4基本参数的均值和方差
参数/指标 均值 变异系数
电机转速(0/s) 200 0.01
柔轮弹性模量(GPA) 206 0.01
柔轮外半径(mm) 31.565 0.01
刚轮刚度系数(mm) 2.1 0.01
刚轮齿轮直径(mm) 65 0.01
输入轴半径(mm) 17.5 0.01
步骤6.2：结合步骤6.1不确定参数的概率密度函数，抽取10000个样本带 入动力学仿真得到输入轴转速、柔轮转速以及输出轴转速统计样本，利用核密 度估计方法得到各自的条件概率密度函数，其中g部件i为部 件层第i个单元的转速。
步骤6.3：根据Bayesian Network理论和公式(4)可知，本实施例中双轴驱 动机构输出转速的概率密度函数可以表示为
f ( g ) = f ( E ) f ( H ) f ( ρ ) f ( d ) f ( F ) f ( ϵ f ) f ( ϵ r ) f ( ϵ w ) × f ( Y M f | X M f , ϵ f ) f ( Y M r | ϵ r ) f ( Y v w | X v w , ϵ w ) - - - ( 34 ) ]]>
其中该等式第一行为各参数的概率密度函数，第二行为条件概率密度函数。
步骤6.4：结合部件层、子系统层的性能试验数据，利用马尔科夫蒙特卡洛 抽样方法和贝叶斯网络模型，利用式(5)和(34)，对可靠性网络模型中所有 节点的联合概率密度函数进行更新，从而得到各部件层、子系统层和系统层响 应参数的概率密度函数：输入轴输出转速的概率密度函数如图12所示，柔轮输 出转速概率密度如图13所示，齿轮接触力的概率密度函数如图14所示，横轴 和纵轴输出转速的概率密度函数如图15、16所示。以及双轴驱动机构输出转速 的概率密度函数如图17所示。
步骤7：利用系统各层次性能或功能响应模型，建立各层节点的广义响应— 阈值干涉模型，分析系统的可靠性。
考虑双轴驱动机构工作以1年为一个周期，前6个月保持匀速转动，第7 个月开始减速到静止状态，随后保持静止6个月。一个工作周期内主要由两个 任务组成：工作6个月基座轴(输出轴1)连续转动时输出转速(不得超过2.09°/s) (任务1)；第7个月起始时刻，基座轴(输出轴1)转动60°、载荷轴(输出轴 2)转动8°(即转动时间不得超过10s)(任务2)。可以看出，两个任务中均不 存在突发型故障，考虑由磨损导致的退化型故障，本案例采用时间离散方法， 将1年作为1个时刻。
步骤7.1：根据阈值要求和公式(6)，分别建立不同任务不同时刻对应的应 力强度干涉模型。结合MC抽样，可求出双轴驱动机构任务1初始时刻的可靠 度R1＝0.9996，任务2初始时刻的可靠度为R2＝0.9974。
步骤7.2：两个任务中均不存在突发型故障，故Pt突发≡0
步骤7.3：根据步骤7.2可知，两个任务的综合可靠度分别为R1＝0.9996， R2＝0.9974
步骤8：建立系统多任务可靠性模型，分析多任务下的系统可靠性
步骤8.1：根据步骤7可知，任务1和2为串联逻辑模型
步骤8.2：计算多任务下的系统可靠度
利用步骤7所计算的单个任务下的系统可靠度Ri(i＝1,2)，按照式(9)计算 多任务下的系统可靠度。
R = Π i = 1 2 R i = 0.9970 - - - ( 35 ) ]]>
考虑间隙变化如图18所示，考虑不同时刻的不同间隙，最终给出12年双 轴驱动机构转速精度可靠度变化规律，如图19所示。
综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保 护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等， 均应包含在本发明的保护范围之内。