CN112560278B - 卫星在轨故障传播和波及效应建模与预测方法及系统 - Google Patents

卫星在轨故障传播和波及效应建模与预测方法及系统 Download PDF

Info

Publication number
CN112560278B
CN112560278B CN202011537503.4A CN202011537503A CN112560278B CN 112560278 B CN112560278 B CN 112560278B CN 202011537503 A CN202011537503 A CN 202011537503A CN 112560278 B CN112560278 B CN 112560278B
Authority
CN
China
Prior art keywords
parameters
formula
correlation
satellite
test
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202011537503.4A
Other languages
English (en)
Other versions
CN112560278A (zh
Inventor
金光
陈思雅
尤杨
马心宇
孙鹏
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN202011537503.4A priority Critical patent/CN112560278B/zh
Publication of CN112560278A publication Critical patent/CN112560278A/zh
Application granted granted Critical
Publication of CN112560278B publication Critical patent/CN112560278B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00Details relating to CAD techniques
    • G06F2111/04Constraint-based CAD

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Monitoring And Testing Of Transmission In General (AREA)

Abstract

本发明涉及卫星状态监测技术领域,公开一种卫星在轨故障传播和波及效应建模与预测方法及系统,以分析出卫星的故障传播路径并确保预测的精度。方法包括:收集卫星遥测参数历史数据;使用格兰杰因果关系模型判断卫星遥测参数之间是否具有格兰杰因果关系;基于互相关函数判断参数之间的格兰杰因果关系的相关强度、相关方向、滞后时间,剔除相关强度低于阈值的相关关系;根据剔除之后剩余的遥测参数之间的关系构造邻接矩阵,画出参数之间的因果关系图;根据遥测参数之间的邻接矩阵,应用解释结构模型,建立卫星遥测参数故障传播层级图;根据格兰杰因果关系模型,定量构造故障传播模型,根据故障传播模型来判断卫星故障发生的波及效应。

Description

卫星在轨故障传播和波及效应建模与预测方法及系统
技术领域
本发明涉及卫星状态监测技术领域,尤其涉及一种卫星在轨故障传播和波及效应建模与预测方法及系统。
背景技术
目前关于故障传播和波及效应的建模与预测方法,主要有两种。一种是直接从系统架构出发,直接构建系统部件传播关系,系统分析方法用树和图来表示系统的结构,需要熟悉系统原理,较为复杂,需要很多的先验知识,难以建立。另一种是从数据的角度出发,从数据来分析出系统的故障传播关系,比如说贝叶斯网络,这是目前用的比较多的方法。但是这些方法在故障传播路径分析中并不总是成功的,而且难以建立起正确的故障传播路径。
由于卫星在轨普遍存在传播特性和波及效应。无论是部件还是分系统故障,若不加控制任其发展,都可能会降低卫星性能、功能或寿命。比如,某卫星发生陀螺的缓变故障后,导致卫星姿态逐渐偏离,控制系统失去自主故障判断机会;随着陀螺持续饱和,陀螺测量输出持续超差,卫星进入全姿态捕获模型,导致其无法执行部分功能。然而,关于故障传播和波及效应的先验数据是缺乏的,需要充分利用遥测参数之间的关系,分析在轨故障的传播特性和波及效应,预测故障及其处置措施对系统或整星性能或功能的影响。实践证明,卫星在轨故障是可以预测的,特别是一些缓变故障,可以通过遥测数据进行早期检测和预防。对于突发故障,虽然难以预测其发生的时机,但是可以通过遥测数据的分析预测其传播过程和波及效应,以便采取有效的故障初值措施降低故障影响。
发明内容
本发明目的在于公开一种卫星在轨故障传播和波及效应建模与预测方法及系统,以分析出卫星的故障传播路径,在卫星发生故障时采取措施降低故障影响,并确保预测的精度。
为达上述目的,本发明公开一种卫星在轨故障传播和波及效应建模与预测方法,包括以下步骤:
步骤S1、收集卫星遥测参数历史数据,使用前值填充的方法对缺失值进行填充;
步骤S2、使用格兰杰因果关系模型判断卫星遥测参数之间是否具有格兰杰因果关系;
步骤S3、基于互相关函数判断参数之间的格兰杰因果关系的相关强度、相关方向、滞后时间,剔除相关强度低于阈值的相关关系;
步骤S4、根据剔除之后剩余的遥测参数之间的关系构造邻接矩阵,画出参数之间的因果关系图;
步骤S5、根据遥测参数之间的邻接矩阵,应用解释结构模型,建立卫星遥测参数故障传播层级图;
步骤S6、根据格兰杰因果关系模型,定量的构造步骤S3剔除之后剩余的卫星遥测参数之间的故障传播模型,根据故障传播模型来判断卫星故障发生的波及效应。
优选地,所述步骤S2中,遥测参数之间的格兰杰因果关系的确定包括以下步骤:
步骤S21、对遥测数据进行去均值处理;
步骤S22、选用ADF对遥测数据进行单位根检验,进行ADF检验的步骤如下:
假设时间序列{yt}服从AR(P)过程:
yt=φ1yt-12yt-2+…+φpyt-pt 公式(1)
式中εt为白噪声,φ1、φ2、…、φp分别为回归系数,p代回归项数量;加入滞后算子B改写上式:
Ψ(B)yt=yt1yt-12yt-2-…-φpyt-p=(1-φ1B-φ2B2-…-φpBp)yt=εt 公式(2)
令:
ρ=φ12+…+φp 公式(3)
ζj=-(φj+1+…+φp);j=1,2,…,p-1 公式(4)
将滞后多项式Ψ(B)分解为:
Ψ(B)=(1-φ1B-φ2B2-…-φpBp)=(1-ρB)-(ζ1B+ζ2B2+…+ζp-1Bp-1)(1-B) 公式(5)
转化为:
ψ(B)yt={(1-ρB)-(ζ1B+ζ2B2+…+ζp-1Bp-1)(1-B)}yt=εt 公式(6)
整理可得:
yt=ρyt-11Δyt-12Δyt-2+…+ζp-1Δyt-p+1t 公式(7)
若服从AR(P)过程的时间序列有且只有一个单位根,则其特征方程为:
1-φ1z-φ2z2-…-φpzp=0 公式(8)
其中,z代表特征根,有且只有一个根则:
Ψ(1)=1-φ12-…-φp=1-ρ=0 公式(9)
等价于:ρ=1;
故对于服从AR(P)过程的时间序列单位根检验等价于检验公式(7)中是否有ρ=1;使得ADF检验根据以下四种情况进行单位根检验:
情况一:随机游走过程判定,时间序列不含常数项不含趋势项,由公式(7)检验存在单位根ρ=1;
情况二:单位根过程判定,时间序列含常数项不含趋势项,在公式(10)估计模型中检验ρ=1,其中,常数项为α;
yt=α+ρyt-11Δyt-12Δyt-2+…+ζp-1Δyt-p+1t 公式(10)
情况三:时间序列含趋势项不含常数项,在公式(10)估计模型中检验ρ=1;
情况四:趋势非平稳过程判定,时间序列含常数项含趋势项,在公式(11)估计模型中检验ρ=1,Lt表示趋势项;
yt=α+Lt+ρyt-11Δyt-12Δyt-2+…+ζp-1Δyt-p+1t 公式(11)
对于情况二,在ρ=1成立时,对公式(10)进行最小二乘估计,可得ρ的超极限估计
Figure BDA0002853956150000031
且有:
Figure BDA0002853956150000032
其中,N()代表正态分布函数;W()代表标准维纳过程;r代表0到1之间的参数;以最小二乘估计
Figure BDA0002853956150000033
代替ζj,得修正统计量ZADF
Figure BDA0002853956150000034
对于检验ρ=1的T统计量,有如下极限分布:
Figure BDA0002853956150000041
其中,
Figure BDA0002853956150000042
代表
Figure BDA0002853956150000043
的标准差,滞后项Δyt的系数估计量
Figure BDA0002853956150000044
极限分布是正态的,故对参数
Figure BDA0002853956150000045
的假设检验可采用T统计量和F统计量进行检验;
对于情况一和情况四,Dickey和Fuller同样证明检验ρ=1的统计量ZADF和T统计量;极限分布均非常规,且均与DF检验中对应极限分布完全一致,故直接使用DF检验对应临界值表;
对于情况三,T统计量的极限分布为常规的T分布,因此可用常规T检验,查T分布表;在实际应用中,用如下三种回归模型进行ADF检验:
Figure BDA0002853956150000046
Figure BDA0002853956150000047
Figure BDA0002853956150000048
故检验单位根的假设ρ=1在上述三种回归模型中变为回归系数
Figure BDA0002853956150000049
如果数据没有通过单位根检验,则需要进行差分处理,直到数据平稳为止,所述差分处理为取时间序列内连续相邻两项之差;
步骤S23、根据BIC准则,确定格兰杰因果关系检验中的最大延迟,BIC可由下式求得:
BIC=kln(n)-2ln(L) 公式(18)
式中,k为模型参数个数,n为样本数量,L为似然函数,取BIC最大时的个数作为格兰杰因果关系检验的最大延迟;
确定最大迟延后,采用多项式拟合的方法构建回归方程,多项式拟合的步骤如下;
多项式定义为:
Figure BDA00028539561500000410
其中,M是多项式的阶数,w0,…,wM是多项式的系数,记做W,使用下式来进行误差评估:
Figure BDA0002853956150000051
其中,E()代表误差,xn代表自变量,tn代表真实值,M代表多项式的数量;根据误差公式来进行数据拟合;
步骤S24、进行杜滨—沃特森检验,检测回归后的残差是否服从正态分布,如果没有通过检验,则参数之间不存在格兰杰因果关系;
所述杜滨—沃特森检验的步骤如下:
通过构建统计量d,依据d值所落的区间范围判定残差间是否存在自相关;d统计量可由下式求得:
Figure BDA0002853956150000052
式中ei为残差,检验结果存在依次由节点0、de、du、4-du、4-de及4所界定的正相关、无法判定、不存在正相关、无法判定和存在负相关共5个区间,de、du为DW中查到的值;
为了改进无法判定的情况,提出新统计量d*,可由下式求得:
d*=a+bdu 公式(22)
式中a、b满足:
Figure BDA0002853956150000053
式中,E(du)、Var(du)可查表得;且:
Figure BDA0002853956150000054
Figure BDA0002853956150000055
P=2(n-1)-tr[(X′AX)(X′X)-1] 公式(26)
Q=2(3n-4)-2tr[X′A2X(X′X)-1]+tr{[X′AX(X′X)-1]2} 公式(27)
Figure BDA0002853956150000061
其中,n为样本容量,k为参数个数,tr为求矩阵的迹,X为样本的值,X′为X的转置;判定:
若:de<d<du,如果d≤d*,则残差间存在正自相关;反之无自相关;
若:4-du<d<4-de,如果d≥d*,则残差间存在负相关,反之无自相关;
步骤S25、在进行格兰杰因果关系检验时,给定两个时间序列{Xs,s=1,2,3,…,t},{Ys,s=1,2,3,…,t},t是时间序列的长度,ε1s是白噪声,p和q分别为时间序列Y和X的最大滞后期数;具体步骤如下:
首先验证X是否是Y变化的格兰杰原因,该检验要求估计以下的无约束的回归模型:
Figure BDA0002853956150000062
以及有约束的回归模型:
Figure BDA0002853956150000063
其中,αi、βi为回归系数;提出检验的零假设,若:β1=β2=…=βq=0,则X不是引起Y变化的格兰杰原因;
分别计算所述无约束的回归模型和所述有约束的回归模型的残差平方和为:RSSu和RSSr
构造F统计量,n为样本容量:
Figure BDA0002853956150000064
检验零假设,如果F≥Fα,那么β12,…,βq显著不为零,则拒绝零假设,即X对Y有影响;否则,则接受原假设,即X对Y没有影响;
交换X与Y的位置,以同样的方式检验Y是否是引起X变化的格兰杰原因;
对一对遥测参数进行格兰杰因果关系检验后,得到它们的相关关系,从而得到节点之间边的关系。
优选地,所述遥测参数之间格兰杰因果关系强弱的确定包括以下步骤:
对遥测时间序列数据{xs,s=1,2,3,…,t}做z-score归一化处理
Figure BDA0002853956150000071
Figure BDA0002853956150000072
其中,μ为均值,δ为标准差;对于两条归一化的曲线G={gs,=1,2,3,…,t},H={hs,=1,2,3,…,t},令:
Figure BDA0002853956150000073
上式中,0的个数是|i|个,其中,-s<i<s;定义Gi与H的内积为:R(Gi,H)=Gi·H,·指的是向量之间的内积,定义相关强度CC(Gi,H):
Figure BDA0002853956150000074
于是:
Figure BDA0002853956150000075
Figure BDA0002853956150000076
则最小值或者最大值的指标分别是:
sc1=argmin-s<i<sCC(Gi,H) 公式(38)
sc2=argmax-s<i<sCC(Gi,H) 公式(39)
令:
Figure BDA0002853956150000077
FCC(G,H)为最终得到的衡量遥测参数之间相关性的元组,FCC(G,H)[0]∈[-1,1],FCC(G,H)[0]的绝对值越接近于1表示参数之间的相关强度越强,正值的FCC(G,H)[0]表示参数之间是正相关,负值的FCC(G,H)[0]表示参数之间是负相关,,FCC(G,H)[1]代表参数之间影响的滞后时间,当|FCC(G,H)[0]|<ε时,ε为阈值,代表参数之间的相关关系不成立,因此它们之间关系被剔除,这样就得到筛选之后的参数之间的相关关系;其中,FCC(G,H)[0]表示FCC(G,H)数组的第一个元素,FCC(G,H)[1]表示FCC(G,H)数组的第二个元素。
优选地,所述卫星遥测参数故障传播层级图的确定如下:
根根筛选之后的卫星参数之间的格兰杰因果关系,然后将其转化成矩阵的形式,构建要素之间的直接关系矩阵,最后利用Warshall算法计算可达矩阵,并进一步计算求解,得到系统的多级递阶层次结构模型。
优选地,所述故障传播模型的构造如下:
根据FCC(G,H)[1]得到的参数间的滞后时间,利用格兰杰因果关系模型,估计如下的回归:
Figure BDA0002853956150000081
Figure BDA0002853956150000082
公式中的e和f为FCC(G,H)[1]得到的卫星遥测参数之间的滞后时间,αj、βj、λj、δj为相应的回归系数,u1i、u2i为噪声,定量构建参数之间的故障传播模型,得出卫星发生故障时的波及效应。
为达上述目的,本发明还公开一种基于改进的格兰杰因果关系模型的卫星在轨故障传播和波及效应建模与预测系统,包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序,所述处理器执行所述计算机程序时实现上述方法的步骤。
本发明具有以下有益效果:
本发明基于互相关函数判断参数之间的格兰杰因果关系的相关强度、相关方向、滞后时间,剔除相关强度低于阈值的相关关系,解决了传统的因果关系建模无法判断参数之间相关关系的强弱问题,并基于剔除相关强度低于阈值的相关关系之后剩余的遥测参数得到故障传播层级图,分析遥测参数之间的层级关系,解决了有关故障传播图无法得出参数的层次结构的问题,并确保了根据故障传播模型来判断卫星故障发生的波及效应的精度。
下面将参照附图,对本发明作进一步详细的说明。
附图说明
构成本申请的一部分的附图用来提供对本发明的进一步理解,本发明的示意性实施例及其说明用于解释本发明,并不构成对本发明的不当限定。在附图中:
图1是本发明实施例的卫星在轨故障传播和波及效应建模与预测方法流程示意图。
图2是本发明实施例DW检验统计量d取值区间检验结果示意图。
图3是本发明实施例改进DW检验统计量d取值区间检验结果示意图。
图4是本发明实施例初始卫星参数之间的格兰杰因果关系图。
图5是本发明实施例筛选之后的卫星参数之间的故障传播图。
具体实施方式
以下结合附图对本发明的实施例进行详细说明,但是本发明可以由权利要求限定和覆盖的多种不同方式实施。
实施例1
本实施例公开一种卫星在轨故障传播和波及效应建模与预测方法,如图1所示,包括以下步骤:
步骤S1、收集卫星遥测参数历史数据,使用前值填充的方法对缺失值进行填充。
步骤S2、使用格兰杰因果关系模型判断卫星遥测参数之间是否具有格兰杰因果关系。
该步骤中,遥测参数之间的格兰杰因果关系的确定包括以下步骤:
步骤S21、对遥测数据进行去均值处理。以此保证卫星遥测数据满足格兰杰因果关系检验中对数据围绕y轴波动的要求。
步骤S22、选用ADF对遥测数据进行单位根检验。为确定卫星遥测数据是否满足格兰杰因果关系检验中待检测数据为平稳序列的要求,需对遥测数据进行单位根检验。进行ADF检验的步骤如下:
假设时间序列{yt}服从AR(P)过程:
yt=φ1yt-12yt-2+…+φpyt-pt 公式(1)
式中εt为白噪声,φ1、φ2、…、φp分别为回归系数,p代回归项数量;加入滞后算子B改写上式:
Ψ(B)yt=yt1yt-12yt-2-…-φpyt-p=(1-φ1B-φ2B2-…-φpBp)yt=εt 公式(2)
令:
ρ=φ12+…+φp 公式(3)
ζj=-(φj+1+…+φp);j=1,2,…,p-1 公式(4)
将滞后多项式Ψ(B)分解为:
Ψ(B)=(1-φ1B-φ2B2-…-φpBp)=(1-ρB)-(ζ1B+ζ2B2+…+ζp-1Bp-1)(1-B) 公式(5)
转化为:
ψ(B)yt={(1-ρB)-(ζ1B+ζ2B2+…+ζp-1Bp-1)(1-B)}yt=εt 公式(6)
整理可得:
yt=ρyt-11Δyt-12Δyt-2+…+ζp-1Δyt-p+1t 公式(7)
若服从AR(P)过程的时间序列有且只有一个单位根,则其特征方程为:
1-φ1z-φ2z2-…-φpzp=0 公式(8)
其中,z代表特征根,有且只有一个根则:
Ψ(1)=1-φ12-…-φp=1-ρ=0 公式(9)
等价于:ρ=1;
故对于服从AR(P)过程的时间序列单位根检验等价于检验公式(7)中是否有ρ=1;使得ADF检验根据以下四种情况进行单位根检验:
情况一:随机游走过程判定,时间序列不含常数项不含趋势项,由公式(7)检验存在单位根ρ=1;
情况二:单位根过程判定,时间序列含常数项不含趋势项,在公式(10)估计模型中检验ρ=1,其中,常数项为α;
yt=α+ρyt-11Δyt-12Δyt-2+…+ζp-1Δyt-p+1t 公式(10)
情况三:时间序列含趋势项不含常数项,在公式(10)估计模型中检验ρ=1;
情况四:趋势非平稳过程判定,时间序列含常数项含趋势项,在公式(11)估计模型中检验ρ=1,Lt表示趋势项;
yt=α+Lt+ρyt-11Δyt-12Δyt-2+…+ζp-1Δyt-p+1t 公式(11)
对于情况二,在ρ=1成立时,对公式(10)进行最小二乘估计,可得ρ的超极限估计
Figure BDA0002853956150000101
且有:
Figure BDA0002853956150000102
其中,N()代表正态分布函数;W()代表标准维纳过程;r代表0到1之间的参数;以最小二乘估计
Figure BDA0002853956150000103
代替ζj,得修正统计量ZADF
Figure BDA0002853956150000104
对于检验ρ=1的T统计量,有如下极限分布:
Figure BDA0002853956150000111
其中,
Figure BDA0002853956150000112
代表
Figure BDA0002853956150000113
的标准差,滞后项Δyt的系数估计量
Figure BDA0002853956150000114
极限分布是正态的,故对参数
Figure BDA0002853956150000115
的假设检验可采用T统计量和F统计量进行检验;
对于情况一和情况四,Dickey和Fuller同样证明检验ρ=1的统计量ZADF和T统计量;极限分布均非常规,且均与DF检验中对应极限分布完全一致,故直接使用DF检验对应临界值表;
对于情况三,T统计量的极限分布为常规的T分布,因此可用常规T检验,查T分布表;在实际应用中,用如下三种回归模型进行ADF检验:
Figure BDA0002853956150000116
Figure BDA0002853956150000117
Figure BDA0002853956150000118
故检验单位根的假设ρ=1在上述三种回归模型中变为回归系数
Figure BDA0002853956150000119
如果数据没有通过单位根检验,则需要进行差分处理,直到数据平稳为止,所述差分处理为取时间序列内连续相邻两项之差;
步骤S23、根据BIC准则,确定格兰杰因果关系检验中的最大延迟。BIC可由下式求得:
BIC=kln(n)-2ln(L) 公式(18)
式中,k为模型参数个数,n为样本数量,L为似然函数,取BIC最大时的个数作为格兰杰因果关系检验的最大延迟;
确定最大迟延后,采用多项式拟合的方法构建回归方程,多项式拟合的步骤如下;
多项式定义为:
Figure BDA00028539561500001110
其中,M是多项式的阶数,w0,…,wM是多项式的系数,记做W,使用下式来进行误差评估:
Figure BDA0002853956150000121
其中,E()代表误差,xn代表自变量,tn代表真实值,M代表多项式的数量;根据误差公式来进行数据拟合;
步骤S24、进行杜滨—沃特森检验,检测回归后的残差是否服从正态分布,如果没有通过检验,则参数之间不存在格兰杰因果关系。
所述杜滨—沃特森检验的步骤如下:
通过构建统计量d,依据d值所落的区间范围判定残差间是否存在自相关;d统计量可由下式求得:
Figure BDA0002853956150000122
式中ei为残差,如图2所示,检验结果存在依次由节点0、de、du、4-du、4-de及4所界定的正相关、无法判定、不存在正相关、无法判定和存在负相关共5个区间,de、du为DW中查到的值;
为了改进无法判定的情况,提出新统计量d*,可由下式求得:
d*=a+bdu 公式(22)
式中a、b满足:
Figure BDA0002853956150000123
式中,E(du)、Var(du)可查表得;且:
Figure BDA0002853956150000124
Figure BDA0002853956150000125
P=2(n-1)-tr[(X′AX)(X′X)-1] 公式(26)
Q=2(3n-4)-2tr[X′A2X(X′X)-1]+tr{[X′AX(X′X)-1]2} 公式(27)
Figure BDA0002853956150000131
其中,n为样本容量,k为参数个数,tr为求矩阵的迹,X为样本的值,X′为X的转置。可判定:
若:de<d<du,如果d≤d*,则残差间存在正自相关;反之无自相关。
若:4-du<d<4-de,如果d≥d*,则残差间存在负相关,反之无自相关。
具体判断逻辑如图3所示,其中,
Figure BDA0002853956150000132
步骤S25、在进行格兰杰因果关系检验时,给定两个时间序列{Xs,s=1,2,3,…,t},{Ys,s=1,2,3,…,t},t是时间序列的长度,ε1s是白噪声,p和q分别为时间序列Y和X的最大滞后期数;具体步骤如下:
首先验证X是否是Y变化的格兰杰原因,该检验要求估计以下的无约束的回归模型:
Figure BDA0002853956150000133
以及有约束的回归模型:
Figure BDA0002853956150000134
其中,αi、βi为回归系数;提出检验的零假设,若:β1=β2=…=βq=0,则X不是引起Y变化的格兰杰原因;
分别计算所述无约束的回归模型和所述有约束的回归模型的残差平方和为:RSSu和RSSr
构造F统计量,n为样本容量:
Figure BDA0002853956150000135
检验零假设,如果F≥Fα,那么β12,…,βq显著不为零,则拒绝零假设,即X对Y有影响;否则,则接受原假设,即X对Y没有影响;
交换X与Y的位置,以同样的方式检验Y是否是引起X变化的格兰杰原因;
对一对遥测参数进行格兰杰因果关系检验后,得到它们的相关关系,从而得到节点之间边的关系。
步骤S3、基于互相关函数判断参数之间的格兰杰因果关系的相关强度、相关方向、滞后时间,剔除相关强度低于阈值的相关关系。
该步骤中,所述遥测参数之间格兰杰因果关系强弱的确定包括以下步骤:
对遥测时间序列数据{xs,s=1,2,3,…,t}做z-score归一化处理
Figure BDA0002853956150000141
Figure BDA0002853956150000142
其中,μ为均值,δ为标准差;对于两条归一化的曲线G={gs,=1,2,3,…,t},H={hs,=1,2,3,…,t},令:
Figure BDA0002853956150000143
上式中,0的个数是|i|个,其中,-s<i<s;定义Gi与H的内积为:R(Gi,H)=Gi·H,·指的是向量之间的内积,定义相关强度CC(Gi,H):
Figure BDA0002853956150000144
于是:
Figure BDA0002853956150000145
Figure BDA0002853956150000146
则最小值或者最大值的指标分别是:
sc1=argmin-s<i<sCC(Gi,H) 公式(38)
sc2=argmax-s<i<sCC(Gi,H) 公式(39)
令:
Figure BDA0002853956150000147
FCC(G,H)为最终得到的衡量遥测参数之间相关性的元组,FCC(G,H)[0]∈[-1,1],FCC(G,H)[0]的绝对值越接近于1表示参数之间的相关强度越强,正值的FCC(G,H)[0]表示参数之间是正相关,负值的FCC(G,H)[0]表示参数之间是负相关,,FCC(G,H)[1]代表参数之间影响的滞后时间,当|FCC(G,H)[0]|<ε时,ε为阈值,代表参数之间的相关关系不成立,因此它们之间关系被剔除,这样就得到筛选之后的参数之间的相关关系;其中,FCC(G,H)[0]表示FCC(G,H)数组的第一个元素,FCC(G,H)[1]表示FCC(G,H)数组的第二个元素。
步骤S4、根据剔除之后剩余的遥测参数之间的关系构造邻接矩阵,画出参数之间的因果关系图。
步骤S5、根据遥测参数之间的邻接矩阵,应用解释结构模型,建立卫星遥测参数故障传播层级图。
该步骤中,所述卫星遥测参数故障传播层级图的确定如下:
根根筛选之后的卫星参数之间的格兰杰因果关系,然后将其转化成矩阵的形式,构建要素之间的直接关系矩阵,最后利用Warshall算法计算可达矩阵,并进一步计算求解,得到系统的多级递阶层次结构模型。
步骤S6、根据格兰杰因果关系模型,定量的构造步骤S3剔除之后剩余的卫星遥测参数之间的故障传播模型,根据故障传播模型来判断卫星故障发生的波及效应。
该步骤中,所述故障传播模型的构造如下:
根据FCC(G,H)[1]得到的参数间的滞后时间,利用格兰杰因果关系模型,估计如下的回归:
Figure BDA0002853956150000151
Figure BDA0002853956150000152
公式中的e和f为FCC(G,H)[1]得到的卫星遥测参数之间的滞后时间,αj、βj、λj、δj为相应的回归系数,u1i、u2i为噪声,定量构建参数之间的故障传播模型,得出卫星发生故障时的波及效应。
实施例2
依托上述实施例中方法,本实施例公开一种卫星在轨故障传播和波及效应建模与预测方法的具体算例,包括以下步骤:
S11:收集卫星遥测参数历史数据;使用前值填充的方法对缺失值进行填充。
参数数据包括温度、电流、电压等遥测参数,因为卫星在轨运行过程中,需要对其部件的多种物理量进行测量,这些物理量就是遥测参数,所得到的数据就是遥测数据,遥测参数以表示形式为时间序列,具体表现为
Figure BDA0002853956150000153
t为时间序列长度,n为参数个数。
S12:使用格兰杰因果关系模型判断卫星遥测参数之间是否具有格兰杰因果关系。遥测参数间的格兰杰因果关系进行两两判断,得出关于所有参数相关关系的邻接矩阵,画出其初步的格兰杰因果关系图,参见图4。
S13:基于互相关函数判断参数之间的格兰杰因果关系的相关强度、相关方向、滞后时间,剔除相关强度较弱的相关关系。
对已经判断出来具有相关关系的参数,使用FCC(G,H)对其进行相关强度、相关方向以及滞后时间的分析,参见表1,其中v,e代表遥测参数的名称。
表1:
v e 相关方向 相关强度 滞后时间
IN3 IN7 正相关 0.98610956 1
Psan IN7 正相关 0.93438408 1
Psas IN7 正相关 0.92383108 1
Psan VN2 正相关 0.81356213 2
Psas VN10 正相关 0.74096304 2
TK25 TG1 正相关 0.7356002 1
VN2 IN7 正相关 0.72157377 3
VN10 IN7 正相关 0.60260963 3
TG3 TK8 正相关 0.48262888 8
TK9 TG1 正相关 0.06938432 51
TG1 TK9 正相关 0.0687796 51
TK8 TK25 正相关 0.01806083 4
TK25 TG3 负相关 -0.02391981 27
TG3 TK25 负相关 -0.02446167 51
TK9 TK25 负相关 -0.0251751 56
TG3 TG1 负相关 -0.03613095 54
TK8 TG1 负相关 -0.05790689 57
TK8 IN7 负相关 -0.11074622 1
VN2 TG1 负相关 -0.16522871 1
VN10 TK25 负相关 -0.16780034 4
TG1 VN10 负相关 -0.20824498 1
TK9 VN10 负相关 -0.32048882 1
VN2 TK8 负相关 -0.37585549 1
TK8 VN10 负相关 -0.41694838 1
TK8 TK9 负相关 -0.46036104 4
S4:根据剔除之后的遥测参数之间的关系构造邻接矩阵,画出参数之间的因果关系图。
把相关强度过小的参数之间的相关关系剔除,得到新的邻接矩阵,画出筛选之后的格兰杰因果关系图如图5所示。
S5:根据遥测参数之间的邻接矩阵,应用解释结构模型,建立卫星遥测参数故障传播层级图。
对筛选之后的邻接矩阵,构建参数之间的直接关系矩阵,最后利用Warshall算法计算可达矩阵,并进一步计算求解,得到系统的多级递阶层次结构模型,得到的层次结构模型如表2所示。
表2:
Figure BDA0002853956150000161
Figure BDA0002853956150000171
S6:根据格兰杰因果关系模型,定量的构造剔除之后剩余的卫星遥测参数之间的故障传播模型,根据故障传播模型来判断卫星故障发生的波及效应,如表3所示。
表3:
Figure BDA0002853956150000172
实施例3
本实施例公开一种基于改进的格兰杰因果关系模型的卫星在轨故障传播和波及效应建模与预测系统,包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序,所述处理器执行所述计算机程序时实现上述两实施例中相对应方法的步骤。
综上,本发明上述各实施例所分别公开的卫星在轨故障传播和波及效应建模与预测方法及系统,至少具有以下有益效果:
本发明基于互相关函数判断参数之间的格兰杰因果关系的相关强度、相关方向、滞后时间,剔除相关强度低于阈值的相关关系,解决了传统的因果关系建模无法判断参数之间相关关系的强弱问题,并基于剔除相关强度低于阈值的相关关系之后剩余的遥测参数得到故障传播层级图,分析遥测参数之间的层级关系,解决了有关故障传播图无法得出参数的层次结构的问题,并确保了根据故障传播模型来判断卫星故障发生的波及效应的精度。
以上所述仅为本发明的优选实施例而已,并不用于限制本发明,对于本领域的技术人员来说,本发明可以有各种更改和变化。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。

Claims (6)

1.一种卫星在轨故障传播和波及效应建模与预测方法,其特征在于,包括以下步骤:
步骤S1、收集卫星遥测参数历史数据,使用前值填充的方法对缺失值进行填充;
步骤S2、使用格兰杰因果关系模型判断卫星遥测参数之间是否具有格兰杰因果关系;
步骤S3、基于互相关函数判断参数之间的格兰杰因果关系的相关强度、相关方向、滞后时间,剔除相关强度低于阈值的相关关系;
步骤S4、根据剔除之后剩余的遥测参数之间的关系构造邻接矩阵,画出参数之间的因果关系图;
步骤S5、根据遥测参数之间的邻接矩阵,应用解释结构模型,建立卫星遥测参数故障传播层级图;
步骤S6、根据格兰杰因果关系模型,定量的构造步骤S3剔除之后剩余的卫星遥测参数之间的故障传播模型,根据故障传播模型来判断卫星故障发生的波及效应。
2.根据权利要求1所述的一种卫星在轨故障传播和波及效应建模与预测方法,其特征在于,所述步骤S2中,遥测参数之间的格兰杰因果关系的确定包括以下步骤:
步骤S21、对遥测数据进行去均值处理;
步骤S22、选用ADF对遥测数据进行单位根检验,进行ADF检验的步骤如下:
假设时间序列{yt}服从AR(P)过程:
yt=φ1yt-12yt-2+…+φpyt-pt 公式(1)
式中εt为白噪声,φ1、φ2、…、φp分别为回归系数,p代回归项数量;加入滞后算子B改写上式:
Ψ(B)yt=yt1yt-12yt-2-…-φpyt-p=(1-φ1B-φ2B2-…-φpBp)yt=εt 公式(2)
令:
ρ=φ12+…+φp 公式(3)
ζj=-(φj+1+…+φp);j=1,2,…,p-1 公式(4)
将滞后多项式Ψ(B)分解为:
Ψ(B)=(1-φ1B-φ2B2-…-φpBp)=(1-ρB)-(ζ1B+ζ2B2+…+ζp-1Bp-1)(1-B) 公式(5)
转化为:
ψ(B)yt={(1-ρB)-(ζ1B+ζ2B2+…+ζp-1Bp-1)(1-B)}yt=εt 公式(6)
整理得:
yt=ρyt-11Δyt-12Δyt-2+…+ζp-1Δyt-p+1t 公式(7)
若服从AR(P)过程的时间序列有且只有一个单位根,则其特征方程为:
1-φ1z-φ2z2-…-φpzp=0 公式(8)
其中,z代表特征根,有且只有一个根则:
Ψ(1)=1-φ12-…-φp=1-ρ=0 公式(9)
等价于:ρ=1;
故对于服从AR(P)过程的时间序列单位根检验等价于检验公式(7)中是否有ρ=1;使得ADF检验根据以下四种情况进行单位根检验:
情况一:随机游走过程判定,时间序列不含常数项不含趋势项,由公式(7)检验存在单位根ρ=1;
情况二:单位根过程判定,时间序列含常数项不含趋势项,在公式(10)估计模型中检验ρ=1,其中,常数项为α;
yt=α+ρyt-11Δyt-12Δyt-2+…+ζp-1Δyt-p+1t 公式(10)
情况三:时间序列含趋势项不含常数项,在公式(10)估计模型中检验ρ=1;
情况四:趋势非平稳过程判定,时间序列含常数项含趋势项,在公式(11)估计模型中检验ρ=1,Lt表示趋势项;
yt=α+Lt+ρyt-11Δyt-12Δyt-2+…+ζp-1Δyt-p+1t 公式(11)
对于情况二,在ρ=1成立时,对公式(10)进行最小二乘估计,得ρ的超极限估计
Figure FDA0002853956140000023
且有:
Figure FDA0002853956140000021
其中,N()代表正态分布函数;W()代表标准维纳过程;r代表0到1之间的参数;以最小二乘估计
Figure FDA0002853956140000022
代替ζj,得修正统计量ZADF
Figure FDA0002853956140000031
对于检验ρ=1的T统计量,有如下极限分布:
Figure FDA0002853956140000032
其中,
Figure FDA0002853956140000033
代表
Figure FDA0002853956140000034
的标准差,滞后项Δyt的系数估计量
Figure FDA0002853956140000035
极限分布是正态的,故对参数
Figure FDA0002853956140000036
的假设检验采用T统计量和F统计量进行检验;
对于情况一和情况四,Dickey和Fuller同样证明检验ρ=1的统计量ZADF和T统计量;极限分布均非常规,且均与DF检验中对应极限分布完全一致,故直接使用DF检验对应临界值表;
对于情况三,T统计量的极限分布为常规的T分布,因此用常规T检验,查T分布表;
在实际应用中,用如下三种回归模型进行ADF检验:
Figure FDA0002853956140000037
Figure FDA0002853956140000038
Figure FDA0002853956140000039
故检验单位根的假设ρ=1在上述三种回归模型中变为回归系数
Figure FDA00028539561400000310
如果数据没有通过单位根检验,则需要进行差分处理,直到数据平稳为止,所述差分处理为取时间序列内连续相邻两项之差;
步骤S23、根据BIC准则,确定格兰杰因果关系检验中的最大延迟,BIC由下式求得:
BIC=kln(n)-2ln(L) 公式(18)
式中,k为模型参数个数,n为样本数量,L为似然函数,取BIC最大时的个数作为格兰杰因果关系检验的最大延迟;
确定最大迟延后,采用多项式拟合的方法构建回归方程,多项式拟合的步骤如下;
多项式定义为:
Figure FDA0002853956140000041
其中,M是多项式的阶数,w0,...,wM是多项式的系数,记做W,使用下式来进行误差评估:
Figure FDA0002853956140000042
其中,E()代表误差,xn代表自变量,tn代表真实值,M代表多项式的数量;根据误差公式来进行数据拟合;
步骤S24、进行杜滨—沃特森检验,检测回归后的残差是否服从正态分布,如果没有通过检验,则参数之间不存在格兰杰因果关系;
所述杜滨—沃特森检验的步骤如下:
通过构建统计量d,依据d值所落的区间范围判定残差间是否存在自相关;d统计量由下式求得:
Figure FDA0002853956140000043
式中ei为残差,检验结果存在依次由节点0、de、du、4-du、4-de及4所界定的正相关、无法判定、不存在正相关、无法判定和存在负相关共5个区间,de、du为DW中查到的值;
为了改进无法判定的情况,提出新统计量d*,由下式求得:
d*=a+bdu 公式(22)
式中a、b满足:
Figure FDA0002853956140000044
式中,E(du)、Var(du)查表得;且:
Figure FDA0002853956140000045
Figure FDA0002853956140000046
P=2(n-1)-tr[(X′AX)(X′X)-1] 公式(26)
Q=2(3n-4)-2tr[X′A2X(X′X)-1]+tr{[X′AX(X′X)-1]2} 公式(27)
Figure FDA0002853956140000051
其中,n为样本容量,k为参数个数,tr为求矩阵的迹,X为样本的值,X′为X的转置;判定:
若:de<d<du,如果d≤d*,则残差间存在正自相关;反之无自相关;
若:4-du<d<4-de,如果d≥d*,则残差间存在负相关,反之无自相关;
步骤S25、在进行格兰杰因果关系检验时,给定两个时间序列{Xs,s=1,2,3,…,t},{Ys,s=1,2,3,…,t},t是时间序列的长度,ε1s是白噪声,p和q分别为时间序列Y和X的最大滞后期数;具体步骤如下:
首先验证X是否是Y变化的格兰杰原因,该检验要求估计以下的无约束的回归模型:
Figure FDA0002853956140000052
以及有约束的回归模型:
Figure FDA0002853956140000053
其中,αi、βi为回归系数;提出检验的零假设,若:β1=β2=…=βq=0,则X不是引起Y变化的格兰杰原因;
分别计算所述无约束的回归模型和所述有约束的回归模型的残差平方和为:RSSu和RSSr
构造F统计量,n为样本容量:
Figure FDA0002853956140000054
检验零假设,如果F≥Fα,那么β12,…,βq显著不为零,则拒绝零假设,即X对Y有影响;否则,则接受原假设,即X对Y没有影响;
交换X与Y的位置,以同样的方式检验Y是否是引起X变化的格兰杰原因;
对一对遥测参数进行格兰杰因果关系检验后,得到它们的相关关系,从而得到节点之间边的关系。
3.根据权利要求1或2所述的一种卫星在轨故障传播和波及效应建模与预测方法,其特征在于,所述遥测参数之间格兰杰因果关系强弱的确定包括以下步骤:
对遥测时间序列数据{xs,s=1,2,3,...,t}做z-score归一化处理
Figure FDA0002853956140000061
Figure FDA0002853956140000062
其中,μ为均值,δ为标准差;对于两条归一化的曲线G={gs,=1,2,3,...,t},H={hs,=1,2,3,...,t},令:
Figure FDA0002853956140000063
上式中,0的个数是|i|个,其中,-s<i<s;定义Gi与H的内积为:R(Gi,H)=Gi·H,·指的是向量之间的内积,定义相关强度CC(Gi,H):
Figure FDA0002853956140000064
于是:
Figure FDA0002853956140000065
Figure FDA0002853956140000066
则最小值或者最大值的指标分别是:
sc1=argmin-s<i<sCC(Gi,H) 公式(38)
sc2=argmax-s<i<sCC(Gi,H) 公式(39)
令:
Figure FDA0002853956140000067
FCC(G,H)为最终得到的衡量遥测参数之间相关性的元组,FCC(G,H)[0]∈[-1,1],FCC(G,H)[0]的绝对值越接近于1表示参数之间的相关强度越强,正值的FCC(G,H)[0]表示参数之间是正相关,负值的FCC(G,H)[0]表示参数之间是负相关,,FCC(G,H)[1]代表参数之间影响的滞后时间,当|FCC(G,H)[0]|<ε时,ε为阈值,代表参数之间的相关关系不成立,因此它们之间关系被剔除,这样就得到筛选之后的参数之间的相关关系;其中,FCC(G,H)[0]表示FCC(G,H)数组的第一个元素,FCC(G,H)[1]表示FCC(G,H)数组的第二个元素。
4.根据权利要求3所述的一种卫星在轨故障传播和波及效应建模与预测方法,其特征在于,所述卫星遥测参数故障传播层级图的确定如下:
根根筛选之后的卫星参数之间的格兰杰因果关系,然后将其转化成矩阵的形式,构建要素之间的直接关系矩阵,最后利用Warshall算法计算可达矩阵,并进一步计算求解,得到系统的多级递阶层次结构模型。
5.根据权利要求4所述的一种卫星在轨故障传播和波及效应建模与预测方法,其特征在于,所述故障传播模型的构造如下:
根据FCC(G,H)[1]得到的参数间的滞后时间,利用格兰杰因果关系模型,估计如下的回归:
Figure FDA0002853956140000071
Figure FDA0002853956140000072
公式中的e和f为FCC(G,H)[1]得到的卫星遥测参数之间的滞后时间,αj、βj、λj、δj为相应的回归系数,u1i、u2i为噪声,定量构建参数之间的故障传播模型,得出卫星发生故障时的波及效应。
6.一种基于改进的格兰杰因果关系模型的卫星在轨故障传播和波及效应建模与预测系统,包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序,其特征在于,所述处理器执行所述计算机程序时实现上述权利要求1至5任一所述方法的步骤。
CN202011537503.4A 2020-12-23 2020-12-23 卫星在轨故障传播和波及效应建模与预测方法及系统 Active CN112560278B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011537503.4A CN112560278B (zh) 2020-12-23 2020-12-23 卫星在轨故障传播和波及效应建模与预测方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011537503.4A CN112560278B (zh) 2020-12-23 2020-12-23 卫星在轨故障传播和波及效应建模与预测方法及系统

Publications (2)

Publication Number Publication Date
CN112560278A CN112560278A (zh) 2021-03-26
CN112560278B true CN112560278B (zh) 2022-05-03

Family

ID=75031565

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011537503.4A Active CN112560278B (zh) 2020-12-23 2020-12-23 卫星在轨故障传播和波及效应建模与预测方法及系统

Country Status (1)

Country Link
CN (1) CN112560278B (zh)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104483958A (zh) * 2014-10-31 2015-04-01 中国石油大学(北京) 一种复杂炼化过程自适应数据驱动故障诊断方法及装置
CN107004212A (zh) * 2014-12-15 2017-08-01 微软技术许可有限责任公司 根据社交媒体和其它数字轨迹对动作、结果和目标实现进行建模
AU2020101478A4 (en) * 2020-06-30 2020-09-10 North China Electric Power University Method for analyzing factors influencing energy business evolution based on ism
CN112052496A (zh) * 2020-08-27 2020-12-08 河海大学 一种基于var模型的高拱坝谷幅变形影响因素分析系统的操作方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104483958A (zh) * 2014-10-31 2015-04-01 中国石油大学(北京) 一种复杂炼化过程自适应数据驱动故障诊断方法及装置
CN107004212A (zh) * 2014-12-15 2017-08-01 微软技术许可有限责任公司 根据社交媒体和其它数字轨迹对动作、结果和目标实现进行建模
AU2020101478A4 (en) * 2020-06-30 2020-09-10 North China Electric Power University Method for analyzing factors influencing energy business evolution based on ism
CN112052496A (zh) * 2020-08-27 2020-12-08 河海大学 一种基于var模型的高拱坝谷幅变形影响因素分析系统的操作方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
《卫星姿控系统基于模糊基函数网络与自回归模型的故障预测》;张茂林;《控制理论与应用》;20110829;472-478 *

Also Published As

Publication number Publication date
CN112560278A (zh) 2021-03-26

Similar Documents

Publication Publication Date Title
Simoen et al. Dealing with uncertainty in model updating for damage assessment: A review
Combastel Merging Kalman filtering and zonotopic state bounding for robust fault detection under noisy environment
Puig et al. Passive robust fault detection of dynamic processes using interval models
EP2342603B1 (en) Method and apparatus for creating state estimation models in machine condition monitoring
Beenstock et al. Polynomial cointegration tests of anthropogenic impact on global warming
Bücher et al. Multiplier bootstrap of tail copulas with applications
CN108171010A (zh) 基于半监督网络嵌入模型的蛋白质复合体检测方法与装置
CN111680398A (zh) 一种基于Holt-Winters模型的单机性能退化预测方法
Wang et al. Noise-dependent ranking of prognostics algorithms based on discrepancy without true damage information
CN108763250B (zh) 一种光伏电站监测数据修复方法
CN112560278B (zh) 卫星在轨故障传播和波及效应建模与预测方法及系统
CN111144650A (zh) 电力负荷预测方法、装置、计算机可读存储介质及设备
Chabane et al. Sensor fault detection and diagnosis using zonotopic set-membership estimation
Shutin et al. Distributed Superresolution Gas Source Localization Based on Poisson Equation
JP2016015586A (ja) Lrlsフィルタ
Simoen et al. Bayesian parameter estimation
Manss et al. Decentralized multi-agent entropy-driven exploration under sparsity constraints
Zhang et al. Sparse learning method with feature selection for sensor placement and response prediction
CN114417912B (zh) 一种野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法
Carratù et al. Cross-Correlation Estimation in Artificial Neural Network for Uncertainty Assessment
Fan et al. Parameter estimation for Hammerstein nonlinear controlled auto-regression models
Kim Local influence on a test of linear hypothesis in multiple regression model
Venugopal et al. Comparison of Methods for Granger Causality Network Estimation
Zhang et al. Adaptive unscented Kalman filter methods for identifying time‐variant parameters via state covariance re‐updating
Cai et al. Control variates with a dimension reduced Bayesian Monte Carlo sampler

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant