CN112560278B - 卫星在轨故障传播和波及效应建模与预测方法及系统 - Google Patents
卫星在轨故障传播和波及效应建模与预测方法及系统 Download PDFInfo
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 55
- 230000000694 effects Effects 0.000 title claims abstract description 36
- 230000001364 causal effect Effects 0.000 claims abstract description 66
- 239000011159 matrix material Substances 0.000 claims abstract description 27
- 238000010586 diagram Methods 0.000 claims abstract description 13
- 238000005314 correlation function Methods 0.000 claims abstract description 7
- 230000008030 elimination Effects 0.000 claims abstract description 5
- 238000003379 elimination reaction Methods 0.000 claims abstract description 5
- 238000012360 testing method Methods 0.000 claims description 54
- 238000007689 inspection Methods 0.000 claims description 18
- 230000008569 process Effects 0.000 claims description 16
- 230000000875 corresponding effect Effects 0.000 claims description 10
- 238000012545 processing Methods 0.000 claims description 9
- 238000001744 unit root test Methods 0.000 claims description 7
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 6
- 238000004590 computer program Methods 0.000 claims description 6
- 230000002596 correlated effect Effects 0.000 claims description 6
- 238000012216 screening Methods 0.000 claims description 6
- 230000006870 function Effects 0.000 claims description 5
- 230000005653 Brownian motion process Effects 0.000 claims description 3
- 230000008859 change Effects 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 238000012937 correction Methods 0.000 claims description 3
- 238000005315 distribution function Methods 0.000 claims description 3
- 238000005295 random walk Methods 0.000 claims description 3
- 230000005654 stationary process Effects 0.000 claims description 3
- 239000013598 vector Substances 0.000 claims description 3
- 238000002759 z-score normalization Methods 0.000 claims description 3
- 230000017105 transposition Effects 0.000 claims description 2
- 238000012544 monitoring process Methods 0.000 abstract description 2
- 102100023600 Fibroblast growth factor receptor 2 Human genes 0.000 description 6
- 101000827688 Homo sapiens Fibroblast growth factor receptor 2 Proteins 0.000 description 6
- 238000004458 analytical method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000006872 improvement Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/04—Constraint-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-1+φ2yt-2+…+φpyt-p+εt 公式(1)
式中εt为白噪声,φ1、φ2、…、φp分别为回归系数,p代回归项数量;加入滞后算子B改写上式:
Ψ(B)yt=yt-φ1yt-1-φ2yt-2-…-φpyt-p=(1-φ1B-φ2B2-…-φpBp)yt=εt 公式(2)
令:
ρ=φ1+φ2+…+φ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-1+ζ1Δyt-1+ζ2Δyt-2+…+ζp-1Δyt-p+1+εt 公式(7)
若服从AR(P)过程的时间序列有且只有一个单位根,则其特征方程为:
1-φ1z-φ2z2-…-φpzp=0 公式(8)
其中,z代表特征根,有且只有一个根则:
Ψ(1)=1-φ1-φ2-…-φp=1-ρ=0 公式(9)
等价于:ρ=1;
故对于服从AR(P)过程的时间序列单位根检验等价于检验公式(7)中是否有ρ=1;使得ADF检验根据以下四种情况进行单位根检验:
情况一:随机游走过程判定,时间序列不含常数项不含趋势项,由公式(7)检验存在单位根ρ=1;
情况二:单位根过程判定,时间序列含常数项不含趋势项,在公式(10)估计模型中检验ρ=1,其中,常数项为α;
yt=α+ρyt-1+ζ1Δyt-1+ζ2Δyt-2+…+ζp-1Δyt-p+1+εt 公式(10)
情况三:时间序列含趋势项不含常数项,在公式(10)估计模型中检验ρ=1;
情况四:趋势非平稳过程判定,时间序列含常数项含趋势项,在公式(11)估计模型中检验ρ=1,Lt表示趋势项;
yt=α+Lt+ρyt-1+ζ1Δyt-1+ζ2Δyt-2+…+ζp-1Δyt-p+1+εt 公式(11)
对于检验ρ=1的T统计量,有如下极限分布:
对于情况一和情况四,Dickey和Fuller同样证明检验ρ=1的统计量ZADF和T统计量;极限分布均非常规,且均与DF检验中对应极限分布完全一致,故直接使用DF检验对应临界值表;
对于情况三,T统计量的极限分布为常规的T分布,因此可用常规T检验,查T分布表;在实际应用中,用如下三种回归模型进行ADF检验:
如果数据没有通过单位根检验,则需要进行差分处理,直到数据平稳为止,所述差分处理为取时间序列内连续相邻两项之差;
步骤S23、根据BIC准则,确定格兰杰因果关系检验中的最大延迟,BIC可由下式求得:
BIC=kln(n)-2ln(L) 公式(18)
式中,k为模型参数个数,n为样本数量,L为似然函数,取BIC最大时的个数作为格兰杰因果关系检验的最大延迟;
确定最大迟延后,采用多项式拟合的方法构建回归方程,多项式拟合的步骤如下;
多项式定义为:
其中,M是多项式的阶数,w0,…,wM是多项式的系数,记做W,使用下式来进行误差评估:
其中,E()代表误差,xn代表自变量,tn代表真实值,M代表多项式的数量;根据误差公式来进行数据拟合;
步骤S24、进行杜滨—沃特森检验,检测回归后的残差是否服从正态分布,如果没有通过检验,则参数之间不存在格兰杰因果关系;
所述杜滨—沃特森检验的步骤如下:
通过构建统计量d,依据d值所落的区间范围判定残差间是否存在自相关;d统计量可由下式求得:
式中ei为残差,检验结果存在依次由节点0、de、du、4-du、4-de及4所界定的正相关、无法判定、不存在正相关、无法判定和存在负相关共5个区间,de、du为DW中查到的值;
为了改进无法判定的情况,提出新统计量d*,可由下式求得:
d*=a+bdu 公式(22)
式中a、b满足:
式中,E(du)、Var(du)可查表得;且:
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)
其中,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变化的格兰杰原因,该检验要求估计以下的无约束的回归模型:
以及有约束的回归模型:
其中,αi、βi为回归系数;提出检验的零假设,若:β1=β2=…=βq=0,则X不是引起Y变化的格兰杰原因;
分别计算所述无约束的回归模型和所述有约束的回归模型的残差平方和为:RSSu和RSSr。
构造F统计量,n为样本容量:
检验零假设,如果F≥Fα,那么β1,β2,…,βq显著不为零,则拒绝零假设,即X对Y有影响;否则,则接受原假设,即X对Y没有影响;
交换X与Y的位置,以同样的方式检验Y是否是引起X变化的格兰杰原因;
对一对遥测参数进行格兰杰因果关系检验后,得到它们的相关关系,从而得到节点之间边的关系。
优选地,所述遥测参数之间格兰杰因果关系强弱的确定包括以下步骤:
对遥测时间序列数据{xs,s=1,2,3,…,t}做z-score归一化处理
其中,μ为均值,δ为标准差;对于两条归一化的曲线G={gs,=1,2,3,…,t},H={hs,=1,2,3,…,t},令:
上式中,0的个数是|i|个,其中,-s<i<s;定义Gi与H的内积为:R(Gi,H)=Gi·H,·指的是向量之间的内积,定义相关强度CC(Gi,H):
于是:
则最小值或者最大值的指标分别是:
sc1=argmin-s<i<sCC(Gi,H) 公式(38)
sc2=argmax-s<i<sCC(Gi,H) 公式(39)
令:
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]得到的参数间的滞后时间,利用格兰杰因果关系模型,估计如下的回归:
公式中的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-1+φ2yt-2+…+φpyt-p+εt 公式(1)
式中εt为白噪声,φ1、φ2、…、φp分别为回归系数,p代回归项数量;加入滞后算子B改写上式:
Ψ(B)yt=yt-φ1yt-1-φ2yt-2-…-φpyt-p=(1-φ1B-φ2B2-…-φpBp)yt=εt 公式(2)
令:
ρ=φ1+φ2+…+φ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-1+ζ1Δyt-1+ζ2Δyt-2+…+ζp-1Δyt-p+1+εt 公式(7)
若服从AR(P)过程的时间序列有且只有一个单位根,则其特征方程为:
1-φ1z-φ2z2-…-φpzp=0 公式(8)
其中,z代表特征根,有且只有一个根则:
Ψ(1)=1-φ1-φ2-…-φp=1-ρ=0 公式(9)
等价于:ρ=1;
故对于服从AR(P)过程的时间序列单位根检验等价于检验公式(7)中是否有ρ=1;使得ADF检验根据以下四种情况进行单位根检验:
情况一:随机游走过程判定,时间序列不含常数项不含趋势项,由公式(7)检验存在单位根ρ=1;
情况二:单位根过程判定,时间序列含常数项不含趋势项,在公式(10)估计模型中检验ρ=1,其中,常数项为α;
yt=α+ρyt-1+ζ1Δyt-1+ζ2Δyt-2+…+ζp-1Δyt-p+1+εt 公式(10)
情况三:时间序列含趋势项不含常数项,在公式(10)估计模型中检验ρ=1;
情况四:趋势非平稳过程判定,时间序列含常数项含趋势项,在公式(11)估计模型中检验ρ=1,Lt表示趋势项;
yt=α+Lt+ρyt-1+ζ1Δyt-1+ζ2Δyt-2+…+ζp-1Δyt-p+1+εt 公式(11)
对于检验ρ=1的T统计量,有如下极限分布:
对于情况一和情况四,Dickey和Fuller同样证明检验ρ=1的统计量ZADF和T统计量;极限分布均非常规,且均与DF检验中对应极限分布完全一致,故直接使用DF检验对应临界值表;
对于情况三,T统计量的极限分布为常规的T分布,因此可用常规T检验,查T分布表;在实际应用中,用如下三种回归模型进行ADF检验:
如果数据没有通过单位根检验,则需要进行差分处理,直到数据平稳为止,所述差分处理为取时间序列内连续相邻两项之差;
步骤S23、根据BIC准则,确定格兰杰因果关系检验中的最大延迟。BIC可由下式求得:
BIC=kln(n)-2ln(L) 公式(18)
式中,k为模型参数个数,n为样本数量,L为似然函数,取BIC最大时的个数作为格兰杰因果关系检验的最大延迟;
确定最大迟延后,采用多项式拟合的方法构建回归方程,多项式拟合的步骤如下;
多项式定义为:
其中,M是多项式的阶数,w0,…,wM是多项式的系数,记做W,使用下式来进行误差评估:
其中,E()代表误差,xn代表自变量,tn代表真实值,M代表多项式的数量;根据误差公式来进行数据拟合;
步骤S24、进行杜滨—沃特森检验,检测回归后的残差是否服从正态分布,如果没有通过检验,则参数之间不存在格兰杰因果关系。
所述杜滨—沃特森检验的步骤如下:
通过构建统计量d,依据d值所落的区间范围判定残差间是否存在自相关;d统计量可由下式求得:
式中ei为残差,如图2所示,检验结果存在依次由节点0、de、du、4-du、4-de及4所界定的正相关、无法判定、不存在正相关、无法判定和存在负相关共5个区间,de、du为DW中查到的值;
为了改进无法判定的情况,提出新统计量d*,可由下式求得:
d*=a+bdu 公式(22)
式中a、b满足:
式中,E(du)、Var(du)可查表得;且:
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)
其中,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变化的格兰杰原因,该检验要求估计以下的无约束的回归模型:
以及有约束的回归模型:
其中,αi、βi为回归系数;提出检验的零假设,若:β1=β2=…=βq=0,则X不是引起Y变化的格兰杰原因;
分别计算所述无约束的回归模型和所述有约束的回归模型的残差平方和为:RSSu和RSSr。
构造F统计量,n为样本容量:
检验零假设,如果F≥Fα,那么β1,β2,…,βq显著不为零,则拒绝零假设,即X对Y有影响;否则,则接受原假设,即X对Y没有影响;
交换X与Y的位置,以同样的方式检验Y是否是引起X变化的格兰杰原因;
对一对遥测参数进行格兰杰因果关系检验后,得到它们的相关关系,从而得到节点之间边的关系。
步骤S3、基于互相关函数判断参数之间的格兰杰因果关系的相关强度、相关方向、滞后时间,剔除相关强度低于阈值的相关关系。
该步骤中,所述遥测参数之间格兰杰因果关系强弱的确定包括以下步骤:
对遥测时间序列数据{xs,s=1,2,3,…,t}做z-score归一化处理
其中,μ为均值,δ为标准差;对于两条归一化的曲线G={gs,=1,2,3,…,t},H={hs,=1,2,3,…,t},令:
上式中,0的个数是|i|个,其中,-s<i<s;定义Gi与H的内积为:R(Gi,H)=Gi·H,·指的是向量之间的内积,定义相关强度CC(Gi,H):
于是:
则最小值或者最大值的指标分别是:
sc1=argmin-s<i<sCC(Gi,H) 公式(38)
sc2=argmax-s<i<sCC(Gi,H) 公式(39)
令:
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]得到的参数间的滞后时间,利用格兰杰因果关系模型,估计如下的回归:
公式中的e和f为FCC(G,H)[1]得到的卫星遥测参数之间的滞后时间,αj、βj、λj、δj为相应的回归系数,u1i、u2i为噪声,定量构建参数之间的故障传播模型,得出卫星发生故障时的波及效应。
实施例2
依托上述实施例中方法,本实施例公开一种卫星在轨故障传播和波及效应建模与预测方法的具体算例,包括以下步骤:
S11:收集卫星遥测参数历史数据;使用前值填充的方法对缺失值进行填充。
参数数据包括温度、电流、电压等遥测参数,因为卫星在轨运行过程中,需要对其部件的多种物理量进行测量,这些物理量就是遥测参数,所得到的数据就是遥测数据,遥测参数以表示形式为时间序列,具体表现为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:
S6:根据格兰杰因果关系模型,定量的构造剔除之后剩余的卫星遥测参数之间的故障传播模型,根据故障传播模型来判断卫星故障发生的波及效应,如表3所示。
表3:
实施例3
本实施例公开一种基于改进的格兰杰因果关系模型的卫星在轨故障传播和波及效应建模与预测系统,包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序,所述处理器执行所述计算机程序时实现上述两实施例中相对应方法的步骤。
综上,本发明上述各实施例所分别公开的卫星在轨故障传播和波及效应建模与预测方法及系统,至少具有以下有益效果:
本发明基于互相关函数判断参数之间的格兰杰因果关系的相关强度、相关方向、滞后时间,剔除相关强度低于阈值的相关关系,解决了传统的因果关系建模无法判断参数之间相关关系的强弱问题,并基于剔除相关强度低于阈值的相关关系之后剩余的遥测参数得到故障传播层级图,分析遥测参数之间的层级关系,解决了有关故障传播图无法得出参数的层次结构的问题,并确保了根据故障传播模型来判断卫星故障发生的波及效应的精度。
以上所述仅为本发明的优选实施例而已,并不用于限制本发明,对于本领域的技术人员来说,本发明可以有各种更改和变化。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (6)
1.一种卫星在轨故障传播和波及效应建模与预测方法,其特征在于,包括以下步骤:
步骤S1、收集卫星遥测参数历史数据,使用前值填充的方法对缺失值进行填充;
步骤S2、使用格兰杰因果关系模型判断卫星遥测参数之间是否具有格兰杰因果关系;
步骤S3、基于互相关函数判断参数之间的格兰杰因果关系的相关强度、相关方向、滞后时间,剔除相关强度低于阈值的相关关系;
步骤S4、根据剔除之后剩余的遥测参数之间的关系构造邻接矩阵,画出参数之间的因果关系图;
步骤S5、根据遥测参数之间的邻接矩阵,应用解释结构模型,建立卫星遥测参数故障传播层级图;
步骤S6、根据格兰杰因果关系模型,定量的构造步骤S3剔除之后剩余的卫星遥测参数之间的故障传播模型,根据故障传播模型来判断卫星故障发生的波及效应。
2.根据权利要求1所述的一种卫星在轨故障传播和波及效应建模与预测方法,其特征在于,所述步骤S2中,遥测参数之间的格兰杰因果关系的确定包括以下步骤:
步骤S21、对遥测数据进行去均值处理;
步骤S22、选用ADF对遥测数据进行单位根检验,进行ADF检验的步骤如下:
假设时间序列{yt}服从AR(P)过程:
yt=φ1yt-1+φ2yt-2+…+φpyt-p+εt 公式(1)
式中εt为白噪声,φ1、φ2、…、φp分别为回归系数,p代回归项数量;加入滞后算子B改写上式:
Ψ(B)yt=yt-φ1yt-1-φ2yt-2-…-φpyt-p=(1-φ1B-φ2B2-…-φpBp)yt=εt 公式(2)
令:
ρ=φ1+φ2+…+φ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-1+ζ1Δyt-1+ζ2Δyt-2+…+ζp-1Δyt-p+1+εt 公式(7)
若服从AR(P)过程的时间序列有且只有一个单位根,则其特征方程为:
1-φ1z-φ2z2-…-φpzp=0 公式(8)
其中,z代表特征根,有且只有一个根则:
Ψ(1)=1-φ1-φ2-…-φp=1-ρ=0 公式(9)
等价于:ρ=1;
故对于服从AR(P)过程的时间序列单位根检验等价于检验公式(7)中是否有ρ=1;使得ADF检验根据以下四种情况进行单位根检验:
情况一:随机游走过程判定,时间序列不含常数项不含趋势项,由公式(7)检验存在单位根ρ=1;
情况二:单位根过程判定,时间序列含常数项不含趋势项,在公式(10)估计模型中检验ρ=1,其中,常数项为α;
yt=α+ρyt-1+ζ1Δyt-1+ζ2Δyt-2+…+ζp-1Δyt-p+1+εt 公式(10)
情况三:时间序列含趋势项不含常数项,在公式(10)估计模型中检验ρ=1;
情况四:趋势非平稳过程判定,时间序列含常数项含趋势项,在公式(11)估计模型中检验ρ=1,Lt表示趋势项;
yt=α+Lt+ρyt-1+ζ1Δyt-1+ζ2Δyt-2+…+ζp-1Δyt-p+1+εt 公式(11)
对于检验ρ=1的T统计量,有如下极限分布:
对于情况一和情况四,Dickey和Fuller同样证明检验ρ=1的统计量ZADF和T统计量;极限分布均非常规,且均与DF检验中对应极限分布完全一致,故直接使用DF检验对应临界值表;
对于情况三,T统计量的极限分布为常规的T分布,因此用常规T检验,查T分布表;
在实际应用中,用如下三种回归模型进行ADF检验:
如果数据没有通过单位根检验,则需要进行差分处理,直到数据平稳为止,所述差分处理为取时间序列内连续相邻两项之差;
步骤S23、根据BIC准则,确定格兰杰因果关系检验中的最大延迟,BIC由下式求得:
BIC=kln(n)-2ln(L) 公式(18)
式中,k为模型参数个数,n为样本数量,L为似然函数,取BIC最大时的个数作为格兰杰因果关系检验的最大延迟;
确定最大迟延后,采用多项式拟合的方法构建回归方程,多项式拟合的步骤如下;
多项式定义为:
其中,M是多项式的阶数,w0,...,wM是多项式的系数,记做W,使用下式来进行误差评估:
其中,E()代表误差,xn代表自变量,tn代表真实值,M代表多项式的数量;根据误差公式来进行数据拟合;
步骤S24、进行杜滨—沃特森检验,检测回归后的残差是否服从正态分布,如果没有通过检验,则参数之间不存在格兰杰因果关系;
所述杜滨—沃特森检验的步骤如下:
通过构建统计量d,依据d值所落的区间范围判定残差间是否存在自相关;d统计量由下式求得:
式中ei为残差,检验结果存在依次由节点0、de、du、4-du、4-de及4所界定的正相关、无法判定、不存在正相关、无法判定和存在负相关共5个区间,de、du为DW中查到的值;
为了改进无法判定的情况,提出新统计量d*,由下式求得:
d*=a+bdu 公式(22)
式中a、b满足:
式中,E(du)、Var(du)查表得;且:
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)
其中,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变化的格兰杰原因,该检验要求估计以下的无约束的回归模型:
以及有约束的回归模型:
其中,αi、βi为回归系数;提出检验的零假设,若:β1=β2=…=βq=0,则X不是引起Y变化的格兰杰原因;
分别计算所述无约束的回归模型和所述有约束的回归模型的残差平方和为:RSSu和RSSr;
构造F统计量,n为样本容量:
检验零假设,如果F≥Fα,那么β1,β2,…,βq显著不为零,则拒绝零假设,即X对Y有影响;否则,则接受原假设,即X对Y没有影响;
交换X与Y的位置,以同样的方式检验Y是否是引起X变化的格兰杰原因;
对一对遥测参数进行格兰杰因果关系检验后,得到它们的相关关系,从而得到节点之间边的关系。
3.根据权利要求1或2所述的一种卫星在轨故障传播和波及效应建模与预测方法,其特征在于,所述遥测参数之间格兰杰因果关系强弱的确定包括以下步骤:
对遥测时间序列数据{xs,s=1,2,3,...,t}做z-score归一化处理
其中,μ为均值,δ为标准差;对于两条归一化的曲线G={gs,=1,2,3,...,t},H={hs,=1,2,3,...,t},令:
上式中,0的个数是|i|个,其中,-s<i<s;定义Gi与H的内积为:R(Gi,H)=Gi·H,·指的是向量之间的内积,定义相关强度CC(Gi,H):
于是:
则最小值或者最大值的指标分别是:
sc1=argmin-s<i<sCC(Gi,H) 公式(38)
sc2=argmax-s<i<sCC(Gi,H) 公式(39)
令:
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算法计算可达矩阵,并进一步计算求解,得到系统的多级递阶层次结构模型。
6.一种基于改进的格兰杰因果关系模型的卫星在轨故障传播和波及效应建模与预测系统,包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序,其特征在于,所述处理器执行所述计算机程序时实现上述权利要求1至5任一所述方法的步骤。
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)
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模型的高拱坝谷幅变形影响因素分析系统的操作方法 |
-
2020
- 2020-12-23 CN CN202011537503.4A patent/CN112560278B/zh active Active
Patent Citations (4)
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)
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 |