基于多探测模式的地波多径目标跟踪方法

1.本发明涉及基于多探测模式的地波多径目标跟踪方法，属于雷达数据处理技术领域。


背景技术：

2.在军事和民用方面，雷达数据处理及目标跟踪都有着广泛的应用，数据关联的正确性和跟踪的准确性影响到整个数据处理过程和目标发现。目前对于电离层杂波的研究大多集中在杂波抑制方面，然而部分电离层散射的电磁波被目标反射，沿天波途径或地波途径返回接收机，经信号处理后可形成携带目标信息的类目标量测，进而导致一个目标生成多个量测的现象，即多径问题。
3.对于多径问题，天波超视距雷达方面已经具有了较为深入的研究，其多径数据关联算法分为以下四种：
4.(1)将单个目标的多径现象在雷达坐标系下形成多条航迹，进而转换到地理坐标系下进行融合，即基于坐标配准进行多径数据关联。这种方法成型较早，但受天波超视距雷达“三低一高”问题影响较为严重，因此存在较大的局限。
5.(2)通过em迭代的方法针对多批数据进行循环处理，通过求解似然概率的方法求取未知参量。这类方法由于需要多批数据积累，因此常用于滑窗处理或离线处理。
6.(3)在贝叶斯框架下，将传统的数据关联算法扩展到多径条件下，进而实现对于多径数据的关联。这种方法在贝叶斯框架下根据概率进行加权融合，可以较好的实现实时处理。
7.(4)在随机有限集理论的基础上，概率假设密度算法(phd)被提出并进行了一系列研究，这类方法目前还处于理论研究期间，具有较深的理论研究价值。
8.与天波超视距雷达的多径数据关联方法相比，目前对于地波超视距雷达的多径问题研究较少。对于地波多径问题，尚不能较好的进行数据关联，进而对目标跟踪的精度较低。


技术实现要素：

9.针对地波超视距雷达存在的多径数据无法准确关联的问题，本发明提供一种基于多探测模式的地波多径目标跟踪方法。
10.本发明的一种基于多探测模式的地波多径目标跟踪方法，包括，
11.步骤一：对第k-1时刻笛卡尔坐标系下的目标状态向量x(k-1|k-1)采样，得到第k-1时刻的采样状态向量，对采样状态向量求解权值；采用目标运动模型基于第k-1时刻的采样状态向量进行一步预测得到第k时刻的采样状态向量预测值；再基于权值和第k时刻的采样状态向量预测值进行目标状态的一步预测，得到第k时刻目标状态向量预测值x(k|k-1)；
12.步骤二：将第k时刻目标状态向量预测值x(k|k-1)转化为雷达坐标系下第q种传播模式下的第k时刻量测一步预测值判断第k时刻量测一步预测值是
否处于地波超视距雷达最大威力范围，保留处于地波超视距雷达最大威力范围内的第q种传播模式下的第k时刻量测一步预测值
13.步骤三：在多探测模式下，根据保留的传播模式建立对应个数的第k时刻传播模式波门，并基于所有传播模式波门建立第k时刻总波门；
14.步骤四：对落入总波门内的所有第k时刻的目标量测值根据传播模式波门的个数建立量测组，根据量测组构造相关事件，得到量测组源于目标的条件概率的表达式；
15.步骤五：基于相关事件计算地波多径传播模式下每个量测组第k时刻第i个采样点的采样状态一步预测值，并计算得到所述相关事件下量测组的第k时刻采样状态一步预测值；再计算得到对应的量测组新息协方差矩阵；
16.基于量测组新息协方差矩阵计算量测组的卡尔曼增益，并进行量测组的第k时刻目标状态向量更新；
17.进一步计算得到量测组源于目标的条件概率；
18.步骤六：更新得到第k时刻笛卡尔坐标系下的目标状态向量x(k|k)和第k时刻的误差协方差矩阵，实现地波多径目标跟踪。
19.根据本发明的基于多探测模式的地波多径目标跟踪方法，步骤一中第k-1时刻的采样状态向量包括2ρ+1个采样状态点：
[0020][0021]
式中χi(k-1|k-1)表示第k-1时刻采样状态向量的第i个采样状态点，κ为尺度参数，ρ+κ≠0并且κ＝ρ(α
2-1)，α为常数，取值范围为0.00001至1之间；p(k-1|k-1)为k-1时刻的误差协方差矩阵；ρ为第k-1时刻目标状态向量x(k-1|k-1)的维数；
[0022]
则2ρ+1个采样状态点对应的权值为：
[0023][0024]
式中w
im
代表第i个采样状态点的均值权值，w
ic
代表第i个采样状态点的协方差权值；β为预设常数。
[0025]
根据本发明的基于多探测模式的地波多径目标跟踪方法，步骤一中，第k时刻采样状态向量第i个采样状态点的预测值χi(k|k-1)为：
[0026]
χi(k|k-1)＝f(k,χi(k-1|k-1))， (3)
[0027]
式中f为目标运动模型；
[0028]
对第k时刻采样状态向量所有采样状态点的预测值进行融合，得到第k时刻目标状
态向量预测值x(k|k-1)：
[0029][0030]
再计算得到第k时刻误差协方差矩阵预测值p(k|k-1)：
[0031][0032]
式中δxi(k|k-1)为第i个采样状态点的预测值χi(k|k-1)与目标状态向量预测值x(k|k-1)的差值：
[0033]
δxi(k|k-1)＝χi(k|k-1)-x(k|k-1)， (6)
[0034]
q(k)为第k时刻过程噪声矩阵。
[0035]
根据本发明的基于多探测模式的地波多径目标跟踪方法，将第k-1时刻笛卡尔坐标系下的目标状态向量x(k-1|k-1)转换为第k-1时刻雷达坐标系下量测值z(k-1|k-1)通过使用量测转换公式(7)至公式(9)实现：
[0036][0037]
式中x为笛卡尔坐标系的x轴位置坐标，y为笛卡尔坐标系的y轴位置坐标；
[0038]
z(k-1|k-1)＝[r,a,d]，
[0039]
式中r为雷达坐标系下的射距，a为雷达坐标系下的方位角，d为雷达坐标系下的径向速度；
[0040][0041]
[0042][0043]
式中r
gg
为地波途径射距，r
e-sg
为e层天地波途径射距，r
f-sg
为f层天地波途径射距，r
e-ss
为e层天波途径射距，r
ef-ss
为ef天波途径射距，he为e层高度，hf为f层高度，d
gg
为地波途径径向速度，d
e-sg
为e层天地波途径径向速度，d
f-sg
为f层天地波途径径向速度，d
e-ss
为e层天波途径径向速度，d
ef-ss
为ef天波途径径向速度；
[0044]
步骤二中按照上述方法由第k时刻目标状态向量预测值x(k|k-1)得到第q种传播模式下的第k时刻量测一步预测值传播模式最多包括五种，分别是地波途径传播模式、e层天地波途径传播模式、f层天地波途径传播模式、e层天波途径和ef天波途径传播模式；
[0045]
表示预测值中的射距：
[0046][0047]
式中r
gm
为地波超视距雷达最大威力范围，为第q种传播模式下射距的测量误差均方差。
[0048]
根据本发明的基于多探测模式的地波多径目标跟踪方法，步骤三中，第k时刻第q种传播模式波门gq(k)为：
[0049][0050]
式中z(k)为第k时刻的目标量测值，sq(k)为第k时刻第q种传播模式的新息协方差；γ为波门范围；
[0051]
则第k时刻总波门g(k)：
[0052][0053]
式中q
max
表示第k时刻扫描周期存在的传播模式数量。
[0054]
本发明的有益效果：本发明考虑了地波多径量测测量精度与检测概率不同的问题，在所建立的地波超视距雷达多径量测模型的基础上，将天波超视距雷达多径数据关联中的多探测模式扩展至地波多径条件下，推导得到了基于ukf的多探测模式状态估计方法，
并利用地波范围对多径路径进行删减，从而实现了基于多探测模式的地波多径数据关联。
[0055]
本发明方法较传统pda方法相比在地波超视距雷达的多径现象下减少了虚假航迹的起始，并有效的提升了目标状态估计的精度。
附图说明
[0056]
图1是本发明所述基于多探测模式的地波多径目标跟踪方法的单步执行流程图；
[0057]
图2是笛卡尔坐标系与雷达坐标系进行转换的地波多径模型；图中f-layer表示f层，e-layer表示e层；
[0058]
图3是假设存在三个波门和四个量测值的波门示意图；
[0059]
图4是仿真实验中目标运动轨迹图；
[0060]
图5是采用pda方法生成航迹的示意图；
[0061]
图6是采用本发明的md-sgm方法生成航迹的示意图；
[0062]
图7是目标1的x轴位置均方根误差；
[0063]
图8是目标1的y轴位置均方根误差；
[0064]
图9是目标2的x轴位置均方根误差；
[0065]
图10是目标2的y轴位置均方根误差。
具体实施方式
[0066]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其它实施例，都属于本发明保护的范围。
[0067]
需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。
[0068]
下面结合附图和具体实施例对本发明作进一步说明，但不作为本发明的限定。
[0069]
具体实施方式一、结合图1至图3所示，本发明提供了一种基于多探测模式的地波多径目标跟踪方法，包括，
[0070]
步骤一：对第k-1时刻笛卡尔坐标系下的目标状态向量x(k-1|k-1)采样，得到第k-1时刻的采样状态向量，对采样状态向量求解权值；采用目标运动模型基于第k-1时刻的采样状态向量进行一步预测得到第k时刻的采样状态向量预测值；再基于权值和第k时刻的采样状态向量预测值进行目标状态的一步预测，得到第k时刻目标状态向量预测值x(k|k-1)；
[0071]
步骤二：将第k时刻目标状态向量预测值x(k|k-1)转化为雷达坐标系下第q种传播模式下的第k时刻量测一步预测值判断第k时刻量测一步预测值是否处于地波超视距雷达最大威力范围，保留处于地波超视距雷达最大威力范围内的第q种传播模式下的第k时刻量测一步预测值
[0072]
步骤三：在多探测模式下，根据保留的传播模式建立对应个数的第k时刻传播模式波门，并基于所有传播模式波门建立第k时刻总波门；
[0073]
步骤四：对落入总波门内的所有第k时刻的目标量测值根据传播模式波门的个数
建立量测组，根据量测组构造相关事件，得到量测组源于目标的条件概率的表达式；
[0074]
步骤五：基于相关事件计算地波多径传播模式下每个量测组第k时刻第i个采样点的采样状态一步预测值，并计算得到所述相关事件下量测组的第k时刻采样状态一步预测值；再计算得到对应的量测组新息协方差矩阵；
[0075]
基于量测组新息协方差矩阵计算量测组的卡尔曼增益，并进行量测组的第k时刻目标状态向量更新；
[0076]
进一步计算得到量测组源于目标的条件概率；
[0077]
步骤六：更新得到第k时刻笛卡尔坐标系下的目标状态向量x(k|k)和第k时刻的误差协方差矩阵，实现地波多径目标跟踪。
[0078]
进一步，步骤一中第k-1时刻的采样状态向量包括2ρ+1个采样状态点：
[0079][0080]
式中χi(k-1|k-1)表示第k-1时刻采样状态向量的第i个采样状态点，κ为尺度参数，ρ+κ≠0并且κ＝ρ(α
2-1)，α为常数，取值范围为0.00001至1之间；p(k-1|k-1)为k-1时刻的误差协方差矩阵；ρ为第k-1时刻目标状态向量x(k-1|k-1)的维数；表示为(ρ+κ)p(k-1|k-1)均方根矩阵的第i行或列；
[0081]
则2ρ+1个采样状态点对应的权值为：
[0082][0083]
式中w
im
代表第i个采样状态点的均值权值，w
ic
代表第i个采样状态点的协方差权值；β为预设常数。
[0084]
步骤一中，第k时刻采样状态向量第i个采样状态点的预测值χi(k|k-1)为：
[0085]
χi(k|k-1)＝f(k,χi(k-1|k-1))， (3)
[0086]
式中f为目标运动模型；
[0087]
对第k时刻采样状态向量所有采样状态点的预测值进行融合，得到第k时刻目标状态向量预测值x(k|k-1)：
[0088][0089]
再计算得到第k时刻误差协方差矩阵预测值p(k|k-1)：
[0090][0091]
式中δxi(k|k-1)为第i个采样状态点的预测值χi(k|k-1)与目标状态向量预测值x(k|k-1)的差值：
[0092]
δxi(k|k-1)＝χi(k|k-1)-x(k|k-1)， (6)
[0093]
q(k)为第k时刻过程噪声矩阵。
[0094]
本实施方式中，结合图2所示，雷达站位于o点，目标位于t点，将第k-1时刻笛卡尔坐标系下的目标状态向量x(k-1|k-1)转换为第k-1时刻雷达坐标系下量测值z(k-1|k-1)通过使用量测转换公式(7)至公式(9)实现：
[0095][0096]
式中x为笛卡尔坐标系的x轴位置坐标，y为笛卡尔坐标系的y轴位置坐标；
[0097]
z(k-1|k-1)＝[r,a,d]，
[0098]
式中r为雷达坐标系下的射距，a为雷达坐标系下的方位角，d为雷达坐标系下的径向速度；假设电离层没有明显倾斜效应，发射和回波信号在同一垂直平面内，所有途径的方位角转换公式可写为：
[0099][0100]
各途径射距的转换公式为：
[0101][0102]
由径向速度为射距的导数得：
[0103][0104]
式中r
gg
为地波途径射距，r
e-sg
为e层天地波途径射距，r
f-sg
为f层天地波途径射距，r
e-ss
为e层天波途径射距，r
ef-ss
为ef天波途径射距，he为电离层e层高度，hf为电离层f层高度，d
gg
为地波途径径向速度，d
e-sg
为e层天地波途径径向速度，d
f-sg
为f层天地波途径径向速度，d
e-ss
为e层天波途径径向速度，d
ef-ss
为ef天波途径径向速度；
[0105]
在多探测模式中，需要针对每一个传播途径进行预测，以此构造量测组的一步预测值，可根据地波超视距雷达探测范围对多探测模式概率数据关联方法进行修正；步骤二中按照上述方法由第k时刻目标状态向量预测值x(k|k-1)得到第q种传播模式下的第k时刻量测一步预测值传播模式最多包括五种，分别是地波途径传播模式、e层天地波途径传播模式、f层天地波途径传播模式、e层天波途径和ef天波途径传播模式；
[0106]
表示预测值中的射距：
[0107][0108]
式中r
gm
为地波超视距雷达最大威力范围，为第q种传播模式下射距的测量误差均方差。进行如下判断：如预测值大于则删减第q种传播模式预测值，否则保留。
[0109]
再进一步，步骤三中，在多探测模式下，根据所存在的传播模式建立多个传播模式波门。其中第k时刻第q种传播模式波门gq(k)为：
[0110][0111]
式中z(k)为第k时刻的目标量测值，sq(k)为第k时刻第q种传播模式的新息协方差；γ为波门范围；
[0112]
此时所有传播模式的波门的并集即为第k时刻总波门g(k)：
[0113][0114]
式中q
max
表示第k时刻扫描周期存在的传播模式数量。
[0115]
本实施方式的步骤四中，假设第k时刻存在三个传播模式，q
max
＝3，形成三个波门g1，g2，g3；假设三个波门形成的总波门范围内落入m(k)个量测值，m(k)＝4，为[z1(k),z2(k),
z3(k),z4(k)]，如图3所示；以表示选择的量测值个数，选择的最多量测值个数以表示在m(k)个量测值中选择个量测值的组合数量，并以量测组表示选择个量测值时的第个量测值组合；当时，量测组代表[z1(k)]，此时仅有z1(k)与目标相关联。
[0116]
构造第k时刻相关事件
[0117][0118]
则得到量测组源于目标的条件概率的表达式：
[0119][0120]
式中p表示概率函数，zk表示到第k时刻被确认量测值的累积集合；
[0121]
根据贝叶斯公式，将公式(13)变形为：
[0122][0123]
式中c为归一化系数，为概率密度函数，为相关事件的条件概率函数。
[0124]
再进一步，步骤五中，利用χi(k|k-1)和公式(7)-(9)对相关事件计算地波多径传播模式下的每个量测组第k时刻第i个采样点的采样状态一步预测值为：
[0125][0126]
式中表示第k时刻选择的个量测值中第个量测值对应的传播模式的量测转换公式；
[0127]
则相关事件下量测组的第k时刻采样状态一步预测值为：
[0128][0129]
对应的量测组新息协方差矩阵为
[0130]
[0131]
式中为采样点的采样状态一步预测值与量测组的第k时刻采样状态一步预测值的差值：
[0132][0133]
为量测组噪声协方差矩阵：
[0134][0135]
表示选择的个量测值中第个量测值对应的传播模式的量测噪声协方差矩阵。步骤五中，量测组的卡尔曼增益为：
[0136][0137]
式中更新得到量测组的第k时刻目标状态向量
[0138][0139]
同时更新得到量测组的第k时刻误差协方差矩阵
[0140][0141]
步骤五中，将量测组源于目标的条件概率分为概率密度函数和相关事件的条件概率函数进行计算：
[0142]
式中pg表示波门概率，n(
·
)表示似然函数，为量测组新息；
[0143]
[0144]
式中u(
·
)为虚假量测值的概率质量函数，为目标在个传播模式下都能够被探测到的概率，pg为波门数，pd为目标被检测到的总概率：
[0145][0146]
最后，步骤六中，对所有相关事件中的量测组根据条件概率对目标状态和误差协方差进行更新，更新得到的第k时刻笛卡尔坐标系下的目标状态向量x(k|k)为：
[0147][0148]
第k时刻的误差协方差矩阵p(k|k)为：
[0149]
p(k|k)＝p1+p2+p3， (27)
[0150]
式中p1、p2和p3均为中间变量：
[0151]
p1＝β
0,1
(k)p(k|k-1)， (28)
[0152][0153]
p3＝-x(k|k)x(k|k)'。 (30)
[0154]
下面通过仿真实验来说明本发明方法的效果：
[0155]
数据关联算法仿真。在笛卡尔坐标系下产生两个目标，其中目标1初始状态为[50km,200m/s,100km,100m/s]，沿近似径向的路径进行运动，目标2初始状态为[300km,-150m/s,50km,100m/s]，沿近似切向的路径进行运动，积累时间为9s，假设e层电离层高度为100km，f层电离层高度为220km。地波模式测量误差为[1km,0.52m/s,0.3
°
]，检测概率为0.9，天地波测量误差为[1.5km,0.52m/s,0.5
°
]，检测概率为0.7，天波模式测量误差为[2km,0.52m/s,0.8
°
]，检测概率为0.5，产生100个周期测量数据，杂波在区域内均匀分布，数量服从泊松分布，密度为0.0001/km2，过程噪声以加速度的形式产生，均方差为0.01m/s，此时笛卡尔坐标系下，目标运动轨迹如图4所示。
[0156]
假设目标已通过地波途径进行起始，若连续五个时刻没有点迹落入相关波门则对航迹进行删除，对未被使用的点迹重新进行起始，首先采用地波超视距雷达中常用的pda进行跟踪，只采用地波模型时目标结果如图5所示，产生的多径量测生成了部分虚假航迹。
[0157]
采用本发明方法的基于多探测模式的地波多径数据关联方法(multi-detection data association methods for ground wave radar multipath data,md-srm)对目标跟踪结果如图6所示，本发明方法考虑了多条路径，利用所有目标量测，可以抑制虚假航迹的生成。
[0158]
单次仿真具有偶然性，进行200次蒙特卡洛仿真对比本发明的md-sgm方法与pda方法的估计精度，采取均方根误差(rmse)作为衡量精度性能的指标。航迹跟踪精度为目标航迹经过滤波后得到的状态信息与实际状态信息各维度之间的均方误差。对于第m次仿真，第t个目标的第v维的滤波误差为：
[0159][0160]
式中x
t
(k)为第t个目标在第k时刻笛卡尔坐标系下的真实状态向量，x
t
(k|k)为第t个目标在第k时刻笛卡尔坐标系下的目标状态向量；
[0161]
则对于m次蒙特卡洛仿真可得第t个目标第v维的航迹跟踪精度(rmse)为：
[0162][0163]
式中traacc
t,v
为第t个目标第v维的航迹跟踪精度。
[0164]
目标1与目标2采用pda方法与md-gm方法的rmse结果如图7至图10所示。经对比，本发明方法与pda方法相比，在x轴、y轴距离上均有提升，滤波收敛后精度提升在10％以上。提升了目标航迹的估计精度。
[0165]
具体实施例：
[0166]
1)建立地波多径模型如图2所示，根据图2中的几何关系得到目标状态由笛卡尔坐标系转换为各途径雷达坐标系的转换公式；
[0167]
2)根据多径数据位置关系对多径路径进行删除；
[0168]
3)将多探测模式扩展到ukf方法中，对目标状态采样并求解权值，根据公式(1)可以确定共得到2ρ+1个采样点迹，对应的权值如公式(2)所示；
[0169]
4)计算目标状态采样点的一步预测x(k|k-1)，计算误差协方差的一步预测p(k|k-1)；
[0170]
5)基于k时刻目标状态x(k|k)根据式(2)至(6)对第i种传播模式进行预测，得到预测值
[0171]
6)在多探测模式下，根据不同传播模式建立多个传播模式波门。假设在第k个扫描周期存在q
max
种传播模式。所有传播模式的波门的并集即为总波门，假设存在三个传播途径形成三个波门，总波门范围内落入4个量测[z1(k),z2(k),z3(k),z4(k)]。如图4所示，q
max
＝3，以表示选择的量测个数，最多选择量测个数表示选择的量测个数，最多选择量测个数表示在m(k)个量测点迹中选择个的组合数量，在多探测模式下可分为以下四种情况：
[0172]
时，此时量测组个数即不存在量测与目标关联；
[0173]
时，此时量测组个数为[q1＝z1(k)]，[q1＝z2(k)]，[q2＝z2(k)]，[q2＝z3(k)]，[q3＝z4(k)]；
[0174]
时，此时量测组个数为[q1＝z1(k),q2＝z2(k)]，[q1＝z1(k),q2＝z3(k)]，[q1＝z1(k),q3＝z4(k)]，[q1＝z2(k),q2＝z3(k)]，[q1＝z2(k),q3＝z4(k)]，[q2＝z2(k),q3＝z4(k)]，[q2＝z3(k),q3＝z4(k)]；
[0175]
时，此时量测组个数为[q1＝z1(k),q2＝z3(k),q3＝z4(k)]，[q1＝z2(k),q2＝z3(k),q3＝z4(k)]，[q1＝z1(k),q2＝z2(k),q3＝z4(k)]。当假设q
max
＝3时，q1表示三种传播模式中的第一种，q2表示三种传播模式中的第二种，q3表示三种传播模式中的第三种。
[0176]
7)进行如下判断：如预测值大于r
gm
+3σr，则删减第种传播模式预测值，否则保留。
[0177]
8)对所有未删减的传播模式以预测值为中心，建立波门如式(10)，并对所有的波门求并集如，得到总波门如式(11)，并将总波门范围内的量测作为备选量测构造多探测量测组。
[0178]
9)根据式(7)-(9)可得地波多径传播模式下的量测组状态采样点的一步预测
由由可算得新息协方差
[0179]
10)计算测量和状态向量的交互协方差和计算卡尔曼增益更新目标状态为更新误差协方差为
[0180]
11)根据多探测量测组构造相关事件此时可得量测组源于目标的条件概率为此时目标被检测到的总概率为
[0181]
12)对每个相关事件中的量测组进行状态更新与误差协方差更新，如式(22)和(23)所示，得到与
[0182]
13)对所有相关事件中的量测组根据条件概率对目标状态和误差协方差进行更新，得到x(k|k)和p(k|k)。
[0183]
虽然在本文中参照了特定的实施方式来描述本发明，但是应该理解的是，这些实施例仅仅是本发明的原理和应用的示例。因此应该理解的是，可以对示例性的实施例进行许多修改，并且可以设计出其他的布置，只要不偏离所附权利要求所限定的本发明的精神和范围。应该理解的是，可以通过不同于原始权利要求所描述的方式来结合不同的从属权利要求和本文中所述的特征。还可以理解的是，结合单独实施例所描述的特征可以使用在其它所述实施例中。

技术特征：
1.一种基于多探测模式的地波多径目标跟踪方法，其特征在于包括，步骤一：对第k-1时刻笛卡尔坐标系下的目标状态向量x(k-1|k-1)采样，得到第k-1时刻的采样状态向量，对采样状态向量求解权值；采用目标运动模型基于第k-1时刻的采样状态向量进行一步预测得到第k时刻的采样状态向量预测值；再基于权值和第k时刻的采样状态向量预测值进行目标状态的一步预测，得到第k时刻目标状态向量预测值x(k|k-1)；步骤二：将第k时刻目标状态向量预测值x(k|k-1)转化为雷达坐标系下第q种传播模式下的第k时刻量测一步预测值判断第k时刻量测一步预测值是否处于地波超视距雷达最大威力范围，保留处于地波超视距雷达最大威力范围内的第q种传播模式下的第k时刻量测一步预测值步骤三：在多探测模式下，根据保留的传播模式建立对应个数的第k时刻传播模式波门，并基于所有传播模式波门建立第k时刻总波门；步骤四：对落入总波门内的所有第k时刻的目标量测值根据传播模式波门的个数建立量测组，根据量测组构造相关事件，得到量测组源于目标的条件概率的表达式；步骤五：基于相关事件计算地波多径传播模式下每个量测组第k时刻第i个采样点的采样状态一步预测值，并计算得到所述相关事件下量测组的第k时刻采样状态一步预测值；再计算得到对应的量测组新息协方差矩阵；基于量测组新息协方差矩阵计算量测组的卡尔曼增益，并进行量测组的第k时刻目标状态向量更新；进一步计算得到量测组源于目标的条件概率；步骤六：更新得到第k时刻笛卡尔坐标系下的目标状态向量x(k|k)和第k时刻的误差协方差矩阵，实现地波多径目标跟踪。2.根据权利要求1所述的基于多探测模式的地波多径目标跟踪方法，其特征在于，步骤一中第k-1时刻的采样状态向量包括2ρ+1个采样状态点：式中χ
i
(k-1|k-1)表示第k-1时刻采样状态向量的第i个采样状态点，κ为尺度参数，ρ+κ≠0并且κ＝ρ(α
2-1)，α为常数，取值范围为0.00001至1之间；p(k-1|k-1)为k-1时刻的误差协方差矩阵；ρ为第k-1时刻目标状态向量x(k-1|k-1)的维数；则2ρ+1个采样状态点对应的权值为：
式中w
im
代表第i个采样状态点的均值权值，w
ic
代表第i个采样状态点的协方差权值；β为预设常数。3.根据权利要求2所述的基于多探测模式的地波多径目标跟踪方法，其特征在于，步骤一中，第k时刻采样状态向量第i个采样状态点的预测值χ
i
(k|k-1)为：χ
i
(k|k-1)＝f(k,χ
i
(k-1|k-1))，(3)式中f为目标运动模型；对第k时刻采样状态向量所有采样状态点的预测值进行融合，得到第k时刻目标状态向量预测值x(k|k-1)：再计算得到第k时刻误差协方差矩阵预测值p(k|k-1)：式中δx
i
(k|k-1)为第i个采样状态点的预测值χ
i
(k|k-1)与目标状态向量预测值x(k|k-1)的差值：δx
i
(k|k-1)＝χ
i
(k|k-1)-x(k|k-1)，(6)q(k)为第k时刻过程噪声矩阵。4.根据权利要求3所述的基于多探测模式的地波多径目标跟踪方法，其特征在于，将第k-1时刻笛卡尔坐标系下的目标状态向量x(k-1|k-1)转换为第k-1时刻雷达坐标系下量测值z(k-1|k-1)通过使用量测转换公式(7)至公式(9)实现：式中x为笛卡尔坐标系的x轴位置坐标，y为笛卡尔坐标系的y轴位置坐标；z(k-1|k-1)＝[r,a,d]，式中r为雷达坐标系下的射距，a为雷达坐标系下的方位角，d为雷达坐标系下的径向速度；
式中r
gg
为地波途径射距，r
e-sg
为e层天地波途径射距，r
f-sg
为f层天地波途径射距，r
e-ss
为e层天波途径射距，r
ef-ss
为ef天波途径射距，h
e
为e层高度，h
f
为f层高度，d
gg
为地波途径径向速度，d
e-sg
为e层天地波途径径向速度，d
f-sg
为f层天地波途径径向速度，d
e-ss
为e层天波途径径向速度，d
ef-ss
为ef天波途径径向速度；步骤二中按照上述方法由第k时刻目标状态向量预测值x(k|k-1)得到第q种传播模式下的第k时刻量测一步预测值传播模式最多包括五种，分别是地波途径传播模式、e层天地波途径传播模式、f层天地波途径传播模式、e层天波途径和ef天波途径传播模式；表示预测值中的射距：式中r
gm
为地波超视距雷达最大威力范围，为第q种传播模式下射距的测量误差均方差。5.根据权利要求4所述的基于多探测模式的地波多径目标跟踪方法，其特征在于，步骤三中，第k时刻第q种传播模式波门g
q
(k)为：
式中z(k)为第k时刻的目标量测值，s
q
(k)为第k时刻第q种传播模式的新息协方差；γ为波门范围；则第k时刻总波门g(k)：式中q
max
表示第k时刻扫描周期存在的传播模式数量。6.根据权利要求5所述的基于多探测模式的地波多径目标跟踪方法，其特征在于，步骤四中，假设第k时刻存在三个传播模式，q
max
＝3，形成三个波门g1，g2，g3；假设三个波门形成的总波门范围内落入m(k)个量测值，m(k)＝4；以表示选择的量测值个数，选择的最多量测值个数以表示在m(k)个量测值中选择个量测值的组合数量，并以量测组表示选择个量测值时的第个量测值组合；构造第k时刻相关事件构造第k时刻相关事件则得到量测组源于目标的条件概率的表达式：式中p表示概率函数，z
k
表示到第k时刻被确认量测值的累积集合；根据贝叶斯公式，将公式(13)变形为：式中c为归一化系数，为概率密度函数，为相关事件的条件概率函数。7.根据权利要求6所述的基于多探测模式的地波多径目标跟踪方法，其特征在于，步骤五中，每个量测组第k时刻第i个采样点的采样状态一步预测值为：式中表示第k时刻选择的个量测值中第个量测值对应的传播模式的量测转换公式；则相关事件下量测组的第k时刻采样状态一步预测值为：
对应的量测组新息协方差矩阵为对应的量测组新息协方差矩阵为式中为采样点的采样状态一步预测值与量测组的第k时刻采样状态一步预测值的差值：的差值：为量测组噪声协方差矩阵：为量测组噪声协方差矩阵：表示选择的个量测值中第个量测值对应的传播模式的量测噪声协方差矩阵。8.根据权利要求7所述的基于多探测模式的地波多径目标跟踪方法，其特征在于，步骤五中，量测组的卡尔曼增益为：式中更新得到量测组的第k时刻目标状态向量更新得到量测组的第k时刻目标状态向量同时更新得到量测组的第k时刻误差协方差矩阵同时更新得到量测组的第k时刻误差协方差矩阵9.根据权利要求8所述的基于多探测模式的地波多径目标跟踪方法，其特征在于，步骤五中，将量测组源于目标的条件概率分为概率密度函数和相关事件的条件概率函数进行计算：
式中p
g
表示波门概率，n(
·
)表示似然函数，为量测组新息；式中u(
·
)为虚假量测值的概率质量函数，为目标在个传播模式下都能够被探测到的概率，p
g
为波门数，p
d
为目标被检测到的总概率：10.根据权利要求9所述的基于多探测模式的地波多径目标跟踪方法，其特征在于，步骤六中，更新得到的第k时刻笛卡尔坐标系下的目标状态向量x(k|k)为：第k时刻的误差协方差矩阵p(k|k)为：p(k|k)＝p1+p2+p3，(27)式中p1、p2和p3均为中间变量：p1＝β
0,1
(k)p(k|k-1)，(28)p3＝-x(k|k)x(k|k)'(30)。

技术总结
一种基于多探测模式的地波多径目标跟踪方法，属于雷达数据处理技术领域。本发明针对地波超视距雷达存在的多径数据无法准确关联的问题。包括：对采样状态向量求解权值，结合采样状态向量预测值得到目标状态向量预测值；将目标状态向量预测值转化到雷达坐标系下得到量测一步预测值，并保留处于地波超视距雷达最大威力范围内的传播模式的量测一步预测值；建立波门并得到总波门；对落入总波门内的目标量测值建立量测组，构造相关事件，得到量测组源于目标的条件概率的表达式，并计算得到量测组源于目标的条件概率；更新得到笛卡尔坐标系下的目标状态向量和误差协方差矩阵，实现地波多径目标跟踪。本发明用于地波多径目标跟踪。本发明用于地波多径目标跟踪。本发明用于地波多径目标跟踪。


技术研发人员：徐龙源 位寅生 童鹏 于雷 徐丝雨
受保护的技术使用者：哈尔滨工业大学
技术研发日：2023.06.19
技术公布日：2023/9/14