一种基于多数据源的机载卫星导航欺骗式干扰检测方法

1.本发明涉及一种卫星导航干扰检测方法。特别是涉及一种基于多数据源的机载卫星导航欺骗式干扰检测方法。


背景技术：

2.国际民航组织(icao)将全球导航卫星系统(gnss)视为新航行系统的主用导航系统，在民航中使用的卫星导航民用信号可为航空器提供连续的定位、测速和授时服务。但是，由于卫星导航民用信号码字信息和信号格式公开，极易受欺骗式干扰影响，这种能引导航空器偏离期望航迹的欺骗式干扰对航空安全构成严重威胁。
3.现有的卫星欺骗式干扰检测技术可分为基于信号处理的欺骗式干扰检测技术和基于信息解算的欺骗式干扰检测技术两大类。基于信号处理的欺骗式干扰检测技术需改变卫星导航信号体制，或者加装额外硬件等。而基于信息解算的欺骗式干扰检测技术是在卫星导航接收端利用检测算法对解算的导航信息进行可靠性判断，无需改变卫星导航信号体制和安装额外设备，是目前欺骗式干扰检测技术的研究重点。
4.常用的基于信息解算的机载欺骗式干扰检测技术由完好性监视技术发展而来，如接收机自主完好性监视技术(raim)和飞机自主完好性监视技术(aaim)。raim算法仅利用卫星导航信息，通过对卫星伪距的一致性检验实现对卫星异常的快速检测和识别，但是只能处理可见卫星6颗以上，同时异常卫星颗数较少的情形。aaim算法基于机载多导航数据源，如惯性导航系统(ins)、测距机(dme)接收机、gnss接收机等，通过多源数据融合实现对卫星异常的检测，可用于可见卫星颗数少的情形。然而，现有基于信息解算的机载欺骗式干扰检测技术较少同时区分卫星故障和卫星受欺骗这两种异常情况，并且对小幅值欺骗式干扰检测的漏警率较高。


技术实现要素：

5.本发明所要解决的技术问题是，提供一种能够用于在可见卫星颗数少、且大部分可见卫星异常情形下，区分卫星故障和欺骗式干扰两种异常情形的基于多数据源的机载卫星导航欺骗式干扰检测方法。
6.本发明所采用的技术方案是：一种基于多数据源的机载卫星导航欺骗式干扰检测方法，包括如下步骤：
7.1)利用载体位置信息计算当前时刻所有可见卫星伪距估计值：
[0008][0009]
其中，为第j颗可见卫星伪距估计值，n为可见卫星颗数，xi,yi,zi为来自惯性导航系统推算的载体位置，x
sat,j
,y
sat,j
,z
sat,j
为来自全球导航卫星系统接收机的第j颗可见卫星的位置；
[0010]
2)利用载体位置信息计算测距机导航台斜距估计值：
[0011][0012]
其中，为第i个测距机导航台的等效伪距，xi,yi,zi为来自惯性导航系统推算的载体位置，x
dme,i
,y
dme,i
,z
dme,i
为来自导航数据库的第i个测距机导航台的位置。
[0013]
3)计算所有可见卫星的异常检测量，具体包括：
[0014]
(a)建立所有可见卫星伪距差量测方程；
[0015]
(b)建立测距机斜距差量测方程；
[0016]
(c)结合卫星伪距差量测方程式(5)和测距机斜距差量测方程式(8)构建量测方程；
[0017]
(d)利用加权最小二乘法计算所有可见卫星的异常检测量
[0018]
(e)将所有可见卫星的异常检测量进行保存；
[0019]
4)计算第j颗可见卫星异常检测统计量λj*；
[0020]
5)对第j颗可见卫星进行异常检测判断；
[0021]
6)计算第j颗可见卫星异常识别统计量；
[0022]
7)对第j颗可见卫星欺骗故障进行识别；
[0023]
8)重复步骤4)到7)，对所有可见卫星进行序贯检测。
[0024]
本发明的一种基于多数据源的机载卫星导航欺骗式干扰检测方法，能够用于在可见卫星颗数少、且大部分可见卫星异常情形下，区分卫星故障和欺骗式干扰两种异常情形。本发明利用ins、dme和gnss提供的机载多源导航数据，应用故障检测领域中的贝叶斯序贯概率比检测算法，利用加权最小二乘法融合dme斜距量测值计算异常检测量并作为先验信息，用于计算异常检测统计量的自适应补偿值，在不增加虚警率的条件下缩短检测时间，有效检测小幅值欺骗式干扰。该方法采用dme辅助检测，能在可见卫星颗数少或大部分可见卫星异常情形下时保证检测有效进行。同时构建异常识别统计量识别卫星异常，用于区分卫星故障和欺骗。
附图说明
[0025]
图1是本发明一种基于多数据源的机载卫星导航欺骗式干扰检测方法的流程图；
[0026]
图2是可见卫星颗数为4颗时，不同伪距偏差下的卫星异常检测概率；
[0027]
图3是可见卫星颗数为5颗时，不同伪距偏差下的卫星异常检测概率；
[0028]
图4是可见卫星颗数为4颗时，不同多普勒频率偏差下的卫星导航欺骗式干扰识别概率；
[0029]
图5是可见卫星颗数为5颗时，不同多普勒频率偏差下的卫星导航欺骗式干扰识别概率。
具体实施方式
[0030]
下面结合实施例和附图对本发明的一种基于多数据源的机载卫星导航欺骗式干扰检测方法做出详细说明。
[0031]
如图1所示，本发明的一种基于多数据源的机载卫星导航欺骗式干扰检测方法，包括如下步骤：
[0032]
1)利用载体位置信息计算当前时刻所有可见卫星伪距估计值：
[0033][0034]
其中，为第j颗可见卫星伪距估计值，n为可见卫星颗数，xi,yi,zi为来自惯性导航系统(ins)推算的载体位置，x
sat,j
,y
sat,j
,z
sat,j
为来自全球导航卫星系统(gnss)接收机的第j颗可见卫星的位置；
[0035]
2)利用载体位置信息计算测距机(dme)导航台斜距估计值：
[0036][0037]
其中，为第i个测距机导航台的等效伪距，xi,yi,zi为来自惯性导航系统(ins)推算的载体位置，x
dme,i
,y
dme,i
,z
dme,i
为来自导航数据库的第i个测距机导航台的位置。
[0038]
3)计算所有可见卫星的异常检测量：
[0039]
所述卫星的异常检测量代表卫星异常引起的伪距差的最小二乘估计。以测距机(dme)斜距差为辅助量，联立卫星伪距差量测方程与测距机(dme)斜距差量测方程，建立量测信息与状态变量之间关系的量测方程，利用加权最小二乘法融合测距机(dme)斜距差量测值计算得到在卫星异常情况下更准确的伪距差估计值。
[0040]
所述可见卫星的异常检测量计算方法包括：
[0041]
(a)建立所有可见卫星伪距差量测方程：
[0042]
在载体真实位置(x,y,z)处对可见卫星伪距估计计算公式(1)右侧进行泰勒级数展开，保留一次项误差，得到：
[0043][0044]
式中，为第j颗可见卫星伪距估计值，r
gnss,j
代表载体到第j颗可见卫星的真实距离，l
gnss,j
、m
gnss,j
、n
gnss,j
为载体与第j颗可见卫星几何连线的三维方向余弦，(δx,δy,δz)＝(x
i-x,y
i-y,z
i-z)代表惯性导航系统(ins)推算的载体位置与载体真实位置之间的偏差，xi,yi,zi为来自惯性导航系统(ins)推算的载体位置，x
sat,j
,y
sat,j
,z
sat,j
为来自全球导航卫星系统(gnss)接收机的第j颗可见卫星的位置；
[0045]
来自全球导航卫星系统(gnss)接收机的伪距测量值ρ
gnss,j
表示为：
[0046]
ρ
gnss,j
＝r
gnss,j
+cδt
gnss-υ
ρ
ꢀꢀꢀ
(4)
[0047]
其中，cδt
gnss
为全球导航卫星系统(gnss)接收机钟差导致的位置偏差，δt
gnss
表示全球导航卫星系统(gnss)接收机钟差，c为光速；υ
ρ
为全球导航卫星系统(gnss)伪距测量噪声；
[0048]
式(3)减去式(4)，得到卫星伪距差量测方程：
[0049][0050]
其中，δρ
gnss,j
为伪距差量测值。
[0051]
(b)建立测距机(dme)斜距差量测方程；
[0052]
在载体真实位置x,y,z处对测距机(dme)导航台斜距估计值计算公式(2)右侧进行泰勒级数展开，保留一次项误差，得到：
[0053][0054]
式中，为第i个测距机(dme)导航台的等效伪距，r
dme,i
代表载体到第i个测距机(dme)导航台的真实距离，l
dme,i
、m
dme,i
、n
dme,i
为载体与第i个测距机(dme)导航台几何连线的三维方向余弦，(δx,δy,δz)＝(x
i-x,y
i-y,z
i-z)代表惯性导航系统(ins)推算的载体位置与载体真实位置之间的偏差，xi,yi,zi为来自惯性导航系统(ins)推算的载体位置，x
dme,i
,y
dme,i
,z
dme,i
为来自导航数据库的第i个测距机导航台的位置。
[0055]
来自测距机(dme)接收机的载体与第i个测距机(dme)导航台斜距测量值d
dme,i
表示为：
[0056]ddme,i
＝r
dme,i
+cδt
dme-υ
d,i
ꢀꢀꢀ
(7)
[0057]
其中，cδt
dme
为测距机(dme)接收机时间误差导致的位置偏差，δt
dme
表示测距机(dme)接收机时间误差，c为光速；υ
d,i
是第i个测距机(dme)斜距测量噪声；
[0058]
式(6)减去式(7)，可得测距机(dme)斜距差量测方程：
[0059][0060]
其中，δρ
dme,i
为斜距差量测值。
[0061]
(c)结合卫星伪距差量测方程式(5)和测距机(dme)斜距差量测方程式(8)构建量测方程：
[0062]
y＝gx+ε
ꢀꢀꢀ
(9)
[0063]
式(9)中，y为(n+2)维列向量量测信息，y＝[δρ
gnss,1
ꢀ…ꢀ
δρ
gnss,n δρ
dme,1 δρ
dme,2
]
t
，上标t表示转置，δρ
gnss,j
为gnss伪距差量测值，j＝1,2,
…
,n，δρ
dme,i
为斜距差量测值，i＝1,2；x为5维列向量状态变量，x＝[δx δy δz
ꢀ‑
cδt
gnss
ꢀ‑
cδt
dme
]
t
，(δx,δy,δz)＝(x
i-x,y
i-y,z
i-z)代表惯性导航系统(ins)推算的载体位置与载体真实位置之间的偏差；ε为(n+2)维测量噪声向量，ε＝[υ
ρ
ꢀ…ꢀ
υ
ρ υ
d,1 υ
d,2
]
t
，全球导航卫星系统(gnss)伪距测量噪声υ
ρ
服从零均值的高斯分布，σ
gnss
为全球导航卫星系统(gnss)伪距测量噪声标准差；两个测距机(dme)斜距测量噪声υ
d,1
、υ
d,2
同样服从零均值的高斯分布，σ
dme,i
表示第i个测距机(dme)斜距测量噪声标准差，i＝1,2；g为(n+2)
×
5维系统矩阵，表示为：
[0064][0065]
式(10)中，l
gnss,j
、m
gnss,j
、n
gnss,j
，j＝1,2,
…
,n，为载体与第j颗可见卫星几何连线的三维方向余弦，表示为：
[0066]
[0067]
其中，
[0068]
式(10)中，l
dme,i
、m
dme,i
、n
dme,i
，i＝1,2，为载体与第i个测距机(dme)导航台几何连线的三维方向余弦，表示为：
[0069][0070]
其中，
[0071]
计算式(11)和式(12)中的三维方向余弦时，式(11)和式(12)中的载体真实位置x,y,z使用惯性导航系统(ins)推算的载体位置xi,yi,zi值。
[0072]
为减小测距机(dme)和全球导航卫星系统(gnss)的测量噪声标准差不同而导致算法对卫星异常引起的伪距差敏感度下降，对所建立的量测方程式(9)进行归一化处理：
[0073]
wy＝wgx+wε
ꢀꢀꢀ
(13)
[0074]
其中，w为(n+2)
×
(n+2)维加权矩阵，
[0075]
由icao附件10规定，全球导航卫星系统(gnss)伪距测量噪声标准差σ
gnss
设定为30米。采用第i个测距机(dme)等效伪距计算第i个测距机(dme)斜距测量噪声标准差
[0076]
(d)利用加权最小二乘法计算所有可见卫星的异常检测量
[0077][0078]
其中，为(n+2)维列向量量测信息y的最小二乘估计的第j行；g为(n+2)
×
5维系统矩阵；w为(n+2)
×
(n+2)维加权矩阵。在可见星正常场景下，异常检测量的统计特性与伪距测量噪声一致，服从零均值的高斯分布；在可见星异常场景下，异常检测量均值不再为零，通过加权最小二乘法融合冗余的测距机(dme)导航数据辅助计算异常检测量，更准确地估计卫星异常引起的伪距偏差值。
[0079]
(e)将异常检测量进行保存。
[0080]
4)计算第j颗可见卫星异常检测统计量
[0081][0082]
其中，λj为第j颗可见卫星原始异常检测统计量，ζj为第j颗可见卫星自适应补偿值；
[0083]
(a)所述的第j颗可见卫星原始异常检测统计量λj的计算方法如下：
[0084]
公式(15)得到的第j颗卫星异常检测量为当前k时刻样本f
k,j
，与历史信息储存
模块中的历史异常检测量构成样本集{f
t,j
|t＝1,2,
…
k}，计算样本集在当前k时刻的均值和方差
[0085][0086]
基于贝叶斯序贯概率比检测算法，当前k时刻第j颗可见卫星原始异常检测统计量λj为：
[0087][0088]
其中，基于二元假设检验，h0为原假设，表示第j颗可见卫星正常；h1为备择假设，表示第j颗可见卫星异常；z
t,j
为利用当前k时刻样本集{f
t,j
|t＝1,2,
…
k}计算的对数概率比。p(f
t,j
|h0)为在原假设h0条件下的概率密度函数，p(f
t,j
|h1)为在备择假设h1条件下的概率密度函数，分别表示为：
[0089][0090]
式(18)的样本均值和方差由下式递推计算得到：
[0091][0092]
(b)第j颗可见卫星自适应补偿值ζj基于贝叶斯(bayes)参数估计理论得到：
[0093][0094]
其中，z
t,j
为利用当前k时刻样本集{f
t,j
|t＝1,2,
…
k}计算的对数概率比，此处作为先验信息；r为每次采样所付出的风险，k1为虚警所产生的风险，虚警产生的风险比每次采样所付出的风险要大。
[0095]
设k1与r有2倍的关系，即r/k1＝0.5。为在备择假设h1条件下利用当前k时刻样本f
k,j
计算的对数概率比的期望值，的计算公式为：
[0096][0097]
5)对第j颗可见卫星进行异常检测判断：
[0098]
将公式(16)得到的第j颗可见卫星异常检测统计量与异常检测门限t
ρ
进行比较，并根据比较结果判断第j颗可见卫星是否正常：
[0099]
若则判断第j颗可见卫星出现异常，即受到欺骗或出现故障，执行步骤6)；
[0100]
若则判断第j颗可见卫星正常，执行步骤8)。
[0101]
所述的异常检测门限t
ρ
由虚警率α和漏警率β确定：
[0102][0103]
按照icao附件10关于机载全球导航卫星系统(gnss)接收机巡航阶段完好性要求，虚警率α＝1.0
×
10-5
，漏警率β＝1.0
×
10-3
，所有可见卫星的异常检测门限设为t
ρ
＝11.5119。
[0104]
6)计算第j颗可见卫星异常识别统计量：
[0105]
在第j颗可见卫星未受欺骗的假设下，第j颗可见卫星异常识别统计量λ
d,j
满足自由度为1的中心卡方分布；在第j颗可见卫星受欺骗的假设下，第j颗可见卫星异常识别统计量λ
d,j
满足自由度为1的非中心卡方分布。根据卡方检测算法，计算第j颗可见卫星异常识别统计量λ
d,j
：
[0106][0107]
其中，δf
d,j
为第j颗可见卫星信号由于欺骗干扰源与接收机的相对运动产生的额外多普勒频移，计算公式为：
[0108][0109]
其中，c为光速；f
t
为卫星信号载波频率；为伪距率量测值：
[0110][0111]
其中，为来自全球导航卫星系统(gnss)接收机的第j颗可见卫星伪距率测量值；为第j颗可见卫星伪距率估计值：
[0112][0113]
其中，是来自惯性导航系统(ins)推算的载体速度，为来自全球导航卫星系统(gnss)接收机的第j颗可见卫星速度。
[0114]
7)对第j颗可见卫星欺骗故障进行识别：
[0115]
将公式(23)得到的第j颗可见卫星异常识别统计量λ
d,j
与异常识别门限td进行比
较，并根据比较结果区分故障和欺骗式干扰：
[0116]
若λ
d,j
＞td，则判断第j颗可见卫星受到了欺骗式干扰，并进行干扰告警；
[0117]
若λ
d,j
≤td，则判断第j颗可见卫星故障，并进行故障告警；
[0118]
所述的异常识别门限td由虚警率α确定：
[0119][0120]
按照icao附件10关于机载全球导航卫星系统(gnss)接收机巡航阶段完好性要求，虚警率α＝1.0
×
10-5
，由卡方分布分位数计算方法，所有可见卫星的异常识别门限为td＝19.5。
[0121]
8)重复步骤4)到7)，对所有可见卫星进行序贯检测。
[0122]
下面是对本发明一种基于多数据源的机载卫星导航欺骗式干扰检测方法的仿真结果及分析：
[0123]
按照icao附件10巡航阶段完好性告警时间的需求，设置30s的告警时间，将受到欺骗后的30秒内异常检测统计量超过异常检测门限作为一次成功检测。卫星受欺骗或卫星故障造成的卫星异常等效为伪距偏差值，伪距偏差值从10米递增到55米，递增步长为5米。
[0124]
图2是在可见卫星颗数少(4颗)的场景下，经过2000次蒙特卡洛实验得到的卫星异常检测概率随伪距偏差值的变化曲线。从图2可以看出，单星异常伪距偏差为10米时，异常检测概率为43％；单星异常伪距偏差大于25米时，异常检测漏警概率低于10-3
。双星异常伪距偏差为10米时，异常检测概率为24％；双星异常伪距偏差大于25米时，异常检测漏警概率低于10-3
。三颗星异常伪距偏差为10米时，异常检测概率趋近于0；三颗星异常伪距偏差为55米时，异常检测漏警概率低于10-3
。四颗星异常伪距偏差为10米时，异常检测概率趋近于0；四颗星异常伪距偏差大于55米时，异常检测漏警概率低于10-3
。不论哪种情况，只要伪距偏差大于55米，异常检测概率均接近100％，而漏警率低于10-3
，满足巡航阶段漏警率β＝1.0
×
10-3
的要求，可对单星或多星异常进行有效检测。
[0125]
图3是在可见卫星颗数少(5颗)的场景下，经过2000次蒙特卡洛实验得到的卫星异常检测概率随伪距偏差值的变化曲线。从图3可以看出，单星异常伪距偏差为10米时，异常检测概率为51.5％；单星异常伪距偏差大于25米时，异常检测漏警概率低于10-3
。双星异常伪距偏差为10米时，异常检测概率为29.5％；双星异常伪距偏差大于25米时，异常检测漏警概率低于10-3
。三颗星异常伪距偏差为10米时，异常检测概率为6％；三颗星异常伪距偏差为25米时，异常检测漏警概率低于10-3
。四颗星异常伪距偏差为10米时，异常检测概率趋近于0；四颗星异常伪距偏差大于50米时，异常检测漏警概率低于10-3
。五颗星异常伪距偏差为10米时，异常检测概率趋近于0；五颗星异常伪距偏差大于55米时，异常检测漏警概率低于10-3
。不论哪种情况，只要伪距偏差大于55米，异常检测概率均接近100％，而漏警率低于10-3
，满足巡航阶段漏警率β＝1.0
×
10-3
的要求，可对单星或多星异常进行有效检测。
[0126]
本发明提供的基于多数据源的机载卫星导航欺骗式干扰检测算法在对卫星异常进行检测后基于异常识别统计量对卫星欺骗和故障进行识别。在伪距偏差为45米条件下，异常识别统计量对卫星欺骗故障识别概率随多普勒频率偏差的变化曲线由图4和图5所示。
[0127]
图4是在可见卫星颗数少(4颗)的场景下，经过2000次蒙特卡洛实验得到的欺骗故障识别概率随多普勒频率偏差的变化曲线。从图4可以看出，单颗星、两颗星受欺骗多普勒
频率偏差为30hz时，欺骗故障识别概率接近100％；三颗星受欺骗多普勒频率偏差为60hz时，欺骗故障识别概率接近100％；四颗星受欺骗多普勒频率偏差为100hz时，欺骗故障识别概率接近100％；
[0128]
图5是在可见卫星颗数少(5颗)的场景下，经过2000次蒙特卡洛实验得到的欺骗故障识别概率随多普勒频率偏差的变化曲线。从图5可以看出，单颗星、两颗星受欺骗多普勒频率偏差为30hz时，欺骗故障识别概率接近100％；三颗星受欺骗多普勒频率偏差为40hz时，欺骗故障识别概率接近100％；四颗星受欺骗多普勒频率偏差为80hz时，欺骗故障识别概率接近100％；五颗星受欺骗多普勒频率偏差为130hz时，欺骗故障识别概率接近100％；
[0129]
表1说明不同欺骗和故障情形下欺骗检测概率的变化情况。以伪距偏差45米，多普勒频率偏差120hz为条件，在不同情形下欺骗检测概率等于对应异常检测概率和欺骗故障识别概率之积。
[0130]
表1
[0131][0132]
在可见卫星颗数少(4颗)的场景下实验4种情形：情形1为4颗星异常，其中4颗星受到欺骗；情形2为4颗星异常，其中3颗星受到欺骗，1颗星故障；情形3为3颗星异常，其中3颗星受到欺骗；情形4为3颗星异常，其中2颗星受到欺骗，1颗星故障。所有情形的欺骗检测概率大于97％，可有效识别欺骗和故障。
[0133]
在可见卫星颗数少(5颗)的场景下实验5种情形：情形5为5颗星异常，其中5颗星受到欺骗；情形6为5颗星异常，其中4颗星受到欺骗，1颗星故障；情形7为4颗星异常，其中4颗星受到欺骗；情形8为4颗星异常，其中3颗星受到欺骗，1颗星故障；情形9为3颗星异常，其中2颗星受到欺骗,1颗星故障。除了情形5和情形6所有可见卫星均异常时欺骗检测概率稍低外，其他情形的欺骗检测概率均大于99％，可有效识别欺骗和故障。

技术特征：
1.一种基于多数据源的机载卫星导航欺骗式干扰检测方法，其特征在于，包括如下步骤：1)利用载体位置信息计算当前时刻所有可见卫星伪距估计值：其中，为第j颗可见卫星伪距估计值，n为可见卫星颗数，x
i
,y
i
,z
i
为来自惯性导航系统推算的载体位置，x
sat,j
,y
sat,j
,z
sat,j
为来自全球导航卫星系统接收机的第j颗可见卫星的位置；2)利用载体位置信息计算测距机导航台斜距估计值：其中，为第i个测距机导航台的等效伪距，x
i
,y
i
,z
i
为来自惯性导航系统推算的载体位置，x
dme,i
,y
dme,i
,z
dme,i
为来自导航数据库的第i个测距机导航台的位置。3)计算所有可见卫星的异常检测量，具体包括：(a)建立所有可见卫星伪距差量测方程；(b)建立测距机斜距差量测方程；(c)结合卫星伪距差量测方程式(5)和测距机斜距差量测方程式(8)构建量测方程；(d)利用加权最小二乘法计算所有可见卫星的异常检测量(e)将所有可见卫星的异常检测量进行保存；4)计算第j颗可见卫星异常检测统计量5)对第j颗可见卫星进行异常检测判断；6)计算第j颗可见卫星异常识别统计量；7)对第j颗可见卫星欺骗故障进行识别；8)重复步骤4)到7)，对所有可见卫星进行序贯检测。2.根据权利要求1所述的一种基于多数据源的机载卫星导航欺骗式干扰检测方法，其特征在于，步骤3)中所述的建立所有可见卫星伪距差量测方程，具体：在载体真实位置(x,y,z)处对可见卫星伪距估计计算公式(1)右侧进行泰勒级数展开，保留一次项误差，得到：式中，为第j颗可见卫星伪距估计值，r
gnss,j
代表载体到第j颗可见卫星的真实距离，l
gnss,j
、m
gnss,j
、n
gnss,j
为载体与第j颗可见卫星几何连线的三维方向余弦，(δx,δy,δz)＝(x
i-x,y
i-y,z
i-z)代表惯性导航系统推算的载体位置与载体真实位置之间的偏差，x
i
,y
i
,z
i
为来自惯性导航系统推算的载体位置，x
sat,j
,y
sat,j
,z
sat,j
为来自全球导航卫星系统接收机的第j颗可见卫星的位置；来自全球导航卫星系统接收机的伪距测量值ρ
gnss,j
表示为：ρ
gnss,j
＝r
gnss,j
+cδt
gnss-υ
ρ
ꢀꢀꢀꢀꢀ
(4)
其中，cδt
gnss
为全球导航卫星系统接收机钟差导致的位置偏差，δt
gnss
表示自全球导航卫星系统接收机钟差，c为光速；υ
ρ
为全球导航卫星系统伪距测量噪声；式(3)减去式(4)，得到卫星伪距差量测方程：其中，δρ
gnss,j
为伪距差量测值。3.根据权利要求1所述的一种基于多数据源的机载卫星导航欺骗式干扰检测方法，其特征在于，步骤3)中所述的建立测距机斜距差量测方程，如下：在载体真实位置x,y,z处对测距机导航台斜距估计值计算公式(2)右侧进行泰勒级数展开，保留一次项误差，得到：式中，为第i个测距机导航台的等效伪距，r
dme,i
代表载体到第i个测距机导航台的真实距离，l
dme,i
、m
dme,i
、n
dme,i
为载体与第i个测距机导航台几何连线的三维方向余弦，(δx,δy,δz)＝(x
i-x,y
i-y,z
i-z)代表惯性导航系统(ins)推算的载体位置与载体真实位置之间的偏差，x
i
,y
i
,z
i
为来自惯性导航系统推算的载体位置，x
dme,i
,y
dme,i
,z
dme,i
为来自导航数据库的第i个测距机导航台的位置；来自测距机接收机的载体与第i个测距机导航台斜距测量值d
dme,i
表示为：d
dme,i
＝r
dme,i
+cδt
dme-υ
d,i
ꢀꢀꢀꢀꢀꢀ
(7)其中，cδt
dme
为测距机接收机时间误差导致的位置偏差，δt
dme
表示测距机接收机时间误差，c为光速；υ
d,i
是第i个测距机斜距测量噪声。式(6)减去式(7)，得到测距机斜距差量测方程：其中，δρ
dme,i
为斜距差量测值。4.根据权利要求1所述的一种基于多数据源的机载卫星导航欺骗式干扰检测方法，其特征在于，步骤3)中所述的结合卫星伪距差量测方程式(5)和测距机斜距差量测方程式(8)构建量测方程，如下：y＝gx+ε
ꢀꢀꢀꢀ
(9)式(9)中，y为(n+2)维列向量量测信息，y＝[δρ
gnss,1
ꢀ…ꢀ
δρ
gnss,n δρ
dme,1 δρ
dme,2
]
t
，上标t表示转置，δρ
gnss,j
为gnss伪距差量测值，j＝1,2,
…
,n，δρ
dme,i
为斜距差量测值，i＝1,2；x为5维列向量状态变量，x＝[δx δy δz
ꢀ‑
cδt
gnss
ꢀ‑
cδt
dme
]
t
，(δx,δy,δz)＝(x
i-x,y
i-y,z
i-z)代表惯性导航系统推算的载体位置与载体真实位置之间的偏差；ε为(n+2)维测量噪声向量，ε＝[υ
ρ
ꢀ…ꢀ
υ
ρ υ
d,1 υ
d,2
]
t
，全球导航卫星系统(gnss)伪距测量噪声υ
ρ
服从零均值的高斯分布，σ
gnss
为全球导航卫星系统(gnss)伪距测量噪声标准差；两个测距机(dme)斜距测量噪声υ
d,1
、υ
d,2
同样服从零均值的高斯分布，σ
dmei
表示第i个测距机(dme)斜距测量噪声标准差，i＝1,2；g为(n+2)
×
5维系统矩阵，表示为：
式(10)中，l
gnss,j
、m
gnss,j
、n
gnss,j
，j＝1,2,
…
,n，为载体与第j颗可见卫星几何连线的三维方向余弦，表示为：其中，式(10)中，l
dme,i
、m
dme,i
、n
dme,i
，i＝1,2，为载体与第i个测距机(dme)导航台几何连线的三维方向余弦，表示为：其中，计算式(11)和式(12)中的三维方向余弦时，式(11)和式(12)中的载体真实位置x,y,z使用惯性导航系统推算的载体位置x
i
,y
i
,z
i
值；为减小测距机和全球导航卫星系统的测量噪声标准差不同而导致算法对卫星异常引起的伪距差敏感度下降，对所建立的量测方程式(9)进行归一化处理：wy＝wgx+wε
ꢀꢀꢀꢀ
(13)其中，w为(n+2)
×
(n+2)维加权矩阵，全球导航卫星系统伪距测量噪声标准差σ
gnss
设定为30米，采用第i个测距机等效伪距计算第i个测距机斜距测量噪声标准差σ
dme,i
：(米）。5.根据权利要求1所述的一种基于多数据源的机载卫星导航欺骗式干扰检测方法，其特征在于，步骤3)中所述的利用加权最小二乘法计算所有可见卫星的异常检测量如下：其中，为(n+2)维列向量量测信息y的最小二乘估计的第j行；g为(n+2)
×
5维系统矩阵；w为(n+2)
×
(n+2)维加权矩阵。6.根据权利要求1所述的一种基于多数据源的机载卫星导航欺骗式干扰检测方法，其
特征在于，步骤4)中所述的计算第j颗可见卫星异常检测统计量公式如下：其中，λ
j
为第j颗可见卫星原始异常检测统计量，ζ
j
为第j颗可见卫星自适应补偿值；(a)所述的第j颗可见卫星原始异常检测统计量λ
j
的计算方法如下：公式(15)得到的第j颗卫星异常检测量为当前k时刻样本f
k,j
，与历史信息储存模块中的历史异常检测量构成样本集{f
t,j
|t＝1,2,
…
k}，计算样本集在当前k时刻的均值和方差方差基于贝叶斯序贯概率比检测算法，当前k时刻第j颗可见卫星原始异常检测统计量λ
j
为：其中，基于二元假设检验，h0为原假设，表示第j颗可见卫星正常；h1为备择假设，表示第j颗可见卫星异常；z
t,j
为利用当前k时刻样本集{f
t,j
|t＝1,2,
…
k}计算的对数概率比；p(f
t,j
|h0)为在原假设h0条件下的概率密度函数，p(f
t,j
|h1)为在备择假设h1条件下的概率密度函数，分别表示为：式(18)的样本均值和方差由下式递推计算得到：(b)第j颗可见卫星自适应补偿值ζ
j
基于贝叶斯参数估计理论得到：其中，z
t,j
为利用当前k时刻样本集{f
t,j
|t＝1,2,
…
k}计算的对数概率比，此处作为先验信息；r为每次采样所付出的风险，k1为虚警所产生的风险，虚警产生的风险比每次采样所付出的风险要大；
设k1与r有2倍的关系，即r/k1＝0.5；为在备择假设h1条件下利用当前k时刻样本f
k,j
计算的对数概率比的期望值，的计算公式为：7.根据权利要求1所述的一种基于多数据源的机载卫星导航欺骗式干扰检测方法，其特征在于，步骤5)所述的对第j颗可见卫星进行异常检测判断，包括：将公式(16)得到的第j颗可见卫星异常检测统计量与异常检测门限t
ρ
进行比较，并根据比较结果判断第j颗可见卫星是否正常：若则判断第j颗可见卫星出现异常，即受到欺骗或出现故障，执行步骤6)；若则判断第j颗可见卫星正常，执行步骤8)；所述的异常检测门限t
ρ
由虚警率α和漏警率β确定：所有可见卫星的异常检测门限设为t
ρ
＝11.5119。8.根据权利要求1所述的一种基于多数据源的机载卫星导航欺骗式干扰检测方法，其特征在于，步骤6)所述的计算第j颗可见卫星异常识别统计量，包括：在第j颗可见卫星未受欺骗的假设下，第j颗可见卫星异常识别统计量λ
d,j
满足自由度为1的中心卡方分布；在第j颗可见卫星受欺骗的假设下，第j颗可见卫星异常识别统计量λ
d,j
满足自由度为1的非中心卡方分布。根据卡方检测算法，计算第j颗可见卫星异常识别统计量λ
d,j
：其中，δf
d,j
为第j颗可见卫星信号由于欺骗干扰源与接收机的相对运动产生的额外多普勒频移，计算公式为：其中，c为光速；f
t
为卫星信号载波频率；为伪距率量测值：其中，为来自全球导航卫星系统(gnss)接收机的第j颗可见卫星伪距率测量值；为第j颗可见卫星伪距率估计值：其中，是来自惯性导航系统推算的载体速度，为来自全球导航卫星系统接收机的第j颗可见卫星速度。
9.根据权利要求1所述的一种基于多数据源的机载卫星导航欺骗式干扰检测方法，其特征在于，步骤7)所述的对第j颗可见卫星欺骗故障进行识别，包括：将公式(23)得到的第j颗可见卫星异常识别统计量λ
d,j
与异常识别门限t
d
进行比较，并根据比较结果区分故障和欺骗式干扰：若λ
d,j
＞t
d
，则判断第j颗可见卫星受到了欺骗式干扰，并进行干扰告警；若λ
d,j
≤t
d
，则判断第j颗可见卫星故障，并进行故障告警；所述的异常识别门限t
d
由虚警率α确定：所有可见卫星的异常识别门限为t
d
＝19.5。

技术总结
一种基于多数据源的机载卫星导航欺骗式干扰检测方法，利用惯性导航系统、测距机和全球导航卫星系统提供的机载多源导航数据，应用故障检测领域中的贝叶斯序贯概率比检测算法(Bayes SPRT)，利用加权最小二乘法融合DME斜距量测值计算异常检测量并作为先验信息，用于计算异常检测统计量的自适应补偿值，在不增加虚警率的条件下缩短检测时间，能够检测小幅值欺骗式干扰。本发明的方法采用测距机辅助检测，能在可见卫星颗数少或大部分可见卫星异常情形下时保证检测有效进行。同时构建异常识别统计量识别卫星异常，用于区分卫星故障和欺骗。本发明能够用于在可见卫星颗数少、且大部分可见卫星异常情形下，区分卫星故障和欺骗式干扰两种异常情形。干扰两种异常情形。干扰两种异常情形。


技术研发人员：钟伦珑 刘炅坡 刘永玉 张卓轩
受保护的技术使用者：中国民航大学
技术研发日：2021.11.28
技术公布日：2022/3/8