CN116736287B - 基于多探测模式的地波多径目标跟踪方法 - Google Patents
基于多探测模式的地波多径目标跟踪方法 Download PDFInfo
- Publication number
- CN116736287B CN116736287B CN202310725543.9A CN202310725543A CN116736287B CN 116736287 B CN116736287 B CN 116736287B CN 202310725543 A CN202310725543 A CN 202310725543A CN 116736287 B CN116736287 B CN 116736287B
- Authority
- CN
- China
- Prior art keywords
- target
- measurement
- state vector
- predicted value
- kth
- 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 57
- 238000001514 detection method Methods 0.000 title claims abstract description 33
- 238000005259 measurement Methods 0.000 claims abstract description 116
- 238000005070 sampling Methods 0.000 claims abstract description 75
- 239000013598 vector Substances 0.000 claims abstract description 69
- 239000011159 matrix material Substances 0.000 claims abstract description 28
- 238000006243 chemical reaction Methods 0.000 claims description 9
- 238000004364 calculation method Methods 0.000 claims description 2
- 230000001186 cumulative effect Effects 0.000 claims description 2
- 238000012545 processing Methods 0.000 abstract description 9
- 239000000523 sample Substances 0.000 description 6
- 239000005433 ionosphere Substances 0.000 description 5
- 238000004422 calculation algorithm Methods 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 238000004088 simulation Methods 0.000 description 4
- 230000000007 visual effect Effects 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 238000000342 Monte Carlo simulation Methods 0.000 description 2
- 238000009825 accumulation Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 230000004927 fusion Effects 0.000 description 2
- 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 description 1
- 230000001133 acceleration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 230000000977 initiatory effect Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000000465 moulding Methods 0.000 description 1
- 230000001629 suppression Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/66—Radar-tracking systems; Analogous systems
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Computer Networks & Wireless Communication (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
一种基于多探测模式的地波多径目标跟踪方法，属于雷达数据处理技术领域。本发明针对地波超视距雷达存在的多径数据无法准确关联的问题。包括：对采样状态向量求解权值，结合采样状态向量预测值得到目标状态向量预测值；将目标状态向量预测值转化到雷达坐标系下得到量测一步预测值，并保留处于地波超视距雷达最大威力范围内的传播模式的量测一步预测值；建立波门并得到总波门；对落入总波门内的目标量测值建立量测组，构造相关事件，得到量测组源于目标的条件概率的表达式，并计算得到量测组源于目标的条件概率；更新得到笛卡尔坐标系下的目标状态向量和误差协方差矩阵，实现地波多径目标跟踪。本发明用于地波多径目标跟踪。
Description
技术领域
本发明涉及基于多探测模式的地波多径目标跟踪方法，属于雷达数据处理技术领域。
背景技术
在军事和民用方面，雷达数据处理及目标跟踪都有着广泛的应用，数据关联的正确性和跟踪的准确性影响到整个数据处理过程和目标发现。目前对于电离层杂波的研究大多集中在杂波抑制方面，然而部分电离层散射的电磁波被目标反射，沿天波途径或地波途径返回接收机，经信号处理后可形成携带目标信息的类目标量测，进而导致一个目标生成多个量测的现象，即多径问题。
对于多径问题，天波超视距雷达方面已经具有了较为深入的研究，其多径数据关联算法分为以下四种：
(1)将单个目标的多径现象在雷达坐标系下形成多条航迹，进而转换到地理坐标系下进行融合，即基于坐标配准进行多径数据关联。这种方法成型较早，但受天波超视距雷达“三低一高”问题影响较为严重，因此存在较大的局限。
(2)通过EM迭代的方法针对多批数据进行循环处理，通过求解似然概率的方法求取未知参量。这类方法由于需要多批数据积累，因此常用于滑窗处理或离线处理。
(3)在贝叶斯框架下，将传统的数据关联算法扩展到多径条件下，进而实现对于多径数据的关联。这种方法在贝叶斯框架下根据概率进行加权融合，可以较好的实现实时处理。
(4)在随机有限集理论的基础上，概率假设密度算法(PHD)被提出并进行了一系列研究，这类方法目前还处于理论研究期间，具有较深的理论研究价值。
与天波超视距雷达的多径数据关联方法相比，目前对于地波超视距雷达的多径问题研究较少。对于地波多径问题，尚不能较好的进行数据关联，进而对目标跟踪的精度较低。
发明内容
针对地波超视距雷达存在的多径数据无法准确关联的问题，本发明提供一种基于多探测模式的地波多径目标跟踪方法。
本发明的一种基于多探测模式的地波多径目标跟踪方法，包括，
步骤一：对第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时刻的误差协方差矩阵，实现地波多径目标跟踪。
根据本发明的基于多探测模式的地波多径目标跟踪方法，步骤一中第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个采样状态点对应的权值为：
式中Wi m代表第i个采样状态点的均值权值，Wi c代表第i个采样状态点的协方差权值；β为预设常数。
根据本发明的基于多探测模式的地波多径目标跟踪方法，步骤一中，第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)：
式中ΔXi(k|k-1)为第i个采样状态点的预测值χi(k|k-1)与目标状态向量预测值X(k|k-1)的差值：
ΔXi(k|k-1)＝χi(k|k-1)-X(k|k-1)， (6)
Q(k)为第k时刻过程噪声矩阵。
根据本发明的基于多探测模式的地波多径目标跟踪方法，将第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为雷达坐标系下的径向速度；
式中RGG为地波途径射距，RE-SG为E层天地波途径射距，RF-SG为F层天地波途径射距，RE-SS为E层天波途径射距，REF-SS为EF天波途径射距，hE为E层高度，hF为F层高度，DGG为地波途径径向速度，DE-SG为E层天地波途径径向速度，DF-SG为F层天地波途径径向速度，DE-SS为E层天波途径径向速度，DEF-SS为EF天波途径径向速度；
步骤二中按照上述方法由第k时刻目标状态向量预测值X(k|k-1)得到第q种传播模式下的第k时刻量测一步预测值传播模式最多包括五种，分别是地波途径传播模式、E层天地波途径传播模式、F层天地波途径传播模式、E层天波途径和EF天波途径传播模式；
表示预测值/>中的射距：
式中RGM为地波超视距雷达最大威力范围，为第q种传播模式下射距的测量误差均方差。
根据本发明的基于多探测模式的地波多径目标跟踪方法，步骤三中，第k时刻第q种传播模式波门Gq(k)为：
式中z(k)为第k时刻的目标量测值，Sq(k)为第k时刻第q种传播模式的新息协方差；γ为波门范围；
则第k时刻总波门G(k)：
式中qmax表示第k时刻扫描周期存在的传播模式数量。
本发明的有益效果：本发明考虑了地波多径量测测量精度与检测概率不同的问题，在所建立的地波超视距雷达多径量测模型的基础上，将天波超视距雷达多径数据关联中的多探测模式扩展至地波多径条件下，推导得到了基于UKF的多探测模式状态估计方法，并利用地波范围对多径路径进行删减，从而实现了基于多探测模式的地波多径数据关联。
本发明方法较传统PDA方法相比在地波超视距雷达的多径现象下减少了虚假航迹的起始，并有效的提升了目标状态估计的精度。
附图说明
图1是本发明所述基于多探测模式的地波多径目标跟踪方法的单步执行流程图；
图2是笛卡尔坐标系与雷达坐标系进行转换的地波多径模型；图中F-layer表示F层，E-layer表示E层；
图3是假设存在三个波门和四个量测值的波门示意图；
图4是仿真实验中目标运动轨迹图；
图5是采用PDA方法生成航迹的示意图；
图6是采用本发明的MD-SGM方法生成航迹的示意图；
图7是目标1的X轴位置均方根误差；
图8是目标1的Y轴位置均方根误差；
图9是目标2的X轴位置均方根误差；
图10是目标2的Y轴位置均方根误差。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其它实施例，都属于本发明保护的范围。
需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。
下面结合附图和具体实施例对本发明作进一步说明，但不作为本发明的限定。
具体实施方式一、结合图1至图3所示，本发明提供了一种基于多探测模式的地波多径目标跟踪方法，包括，
步骤一：对第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时刻的误差协方差矩阵，实现地波多径目标跟踪。
进一步，步骤一中第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)的维数；表示为(ρ+κ)P(k-1|k-1)均方根矩阵的第i行或列；
则2ρ+1个采样状态点对应的权值为：
式中Wi m代表第i个采样状态点的均值权值，Wi c代表第i个采样状态点的协方差权值；β为预设常数。
步骤一中，第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)：
式中ΔXi(k|k-1)为第i个采样状态点的预测值χi(k|k-1)与目标状态向量预测值X(k|k-1)的差值：
ΔXi(k|k-1)＝χi(k|k-1)-X(k|k-1)， (6)
Q(k)为第k时刻过程噪声矩阵。
本实施方式中，结合图2所示，雷达站位于O点，目标位于T点，将第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为雷达坐标系下的径向速度；假设电离层没有明显倾斜效应，发射和回波信号在同一垂直平面内，所有途径的方位角转换公式可写为：
各途径射距的转换公式为：
由径向速度为射距的导数得：
式中RGG为地波途径射距，RE-SG为E层天地波途径射距，RF-SG为F层天地波途径射距，RE-SS为E层天波途径射距，REF-SS为EF天波途径射距，hE为电离层E层高度，hF为电离层F层高度，DGG为地波途径径向速度，DE-SG为E层天地波途径径向速度，DF-SG为F层天地波途径径向速度，DE-SS为E层天波途径径向速度，DEF-SS为EF天波途径径向速度；
在多探测模式中，需要针对每一个传播途径进行预测，以此构造量测组的一步预测值，可根据地波超视距雷达探测范围对多探测模式概率数据关联方法进行修正；步骤二中按照上述方法由第k时刻目标状态向量预测值X(k|k-1)得到第q种传播模式下的第k时刻量测一步预测值传播模式最多包括五种，分别是地波途径传播模式、E层天地波途径传播模式、F层天地波途径传播模式、E层天波途径和EF天波途径传播模式；
表示预测值/>中的射距：
式中RGM为地波超视距雷达最大威力范围，为第q种传播模式下射距的测量误差均方差。进行如下判断：如预测值大于/>则删减第q种传播模式预测值，否则保留。
再进一步，步骤三中，在多探测模式下，根据所存在的传播模式建立多个传播模式波门。其中第k时刻第q种传播模式波门Gq(k)为：
式中z(k)为第k时刻的目标量测值，Sq(k)为第k时刻第q种传播模式的新息协方差；γ为波门范围；
此时所有传播模式的波门的并集即为第k时刻总波门G(k)：
式中qmax表示第k时刻扫描周期存在的传播模式数量。
本实施方式的步骤四中，假设第k时刻存在三个传播模式，qmax＝3，形成三个波门G1，G2，G3；假设三个波门形成的总波门范围内落入m(k)个量测值，m(k)＝4，为[z1(k),z2(k),z3(k),z4(k)]，如图3所示；以表示选择的量测值个数，选择的最多量测值个数/>以/>表示在m(k)个量测值中选择/>个量测值的组合数量，并以量测组/>表示选择/>个量测值时的第/>个量测值组合；当/>时，量测组/>代表[z1(k)]，此时仅有z1(k)与目标相关联。
构造第k时刻相关事件
则得到量测组源于目标的条件概率/>的表达式：
式中p表示概率函数，Zk表示到第k时刻被确认量测值的累积集合；
根据贝叶斯公式，将公式(13)变形为：
式中c为归一化系数，为概率密度函数，为相关事件的条件概率函数。
再进一步，步骤五中，利用χi(k|k-1)和公式(7)-(9)对相关事件计算地波多径传播模式下的每个量测组第k时刻第i个采样点的采样状态一步预测值/>为：
式中表示第k时刻选择的/>个量测值中第/>个量测值对应的传播模式的量测转换公式；
则相关事件下量测组的第k时刻采样状态一步预测值为：
对应的量测组新息协方差矩阵为
式中为采样点的采样状态一步预测值/>与量测组的第k时刻采样状态一步预测值/>的差值：
为量测组噪声协方差矩阵：
表示选择的/>个量测值中第/>个量测值对应的传播模式的量测噪声协方差矩阵。步骤五中，量测组的卡尔曼增益/>为：
式中更新得到量测组的第k时刻目标状态向量/>
同时更新得到量测组的第k时刻误差协方差矩阵
步骤五中，将量测组源于目标的条件概率分为概率密度函数和相关事件的条件概率函数/>进行计算：/>
式中PG表示波门概率，N(·)表示似然函数，为量测组新息；
/>
式中u(·)为虚假量测值的概率质量函数，为目标在/>个传播模式下都能够被探测到的概率，PG为波门数，PD为目标被检测到的总概率：
最后，步骤六中，对所有相关事件中的量测组根据条件概率对目标状态和误差协方差进行更新，更新得到的第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)
下面通过仿真实验来说明本发明方法的效果：
数据关联算法仿真。在笛卡尔坐标系下产生两个目标，其中目标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所示。
假设目标已通过地波途径进行起始，若连续五个时刻没有点迹落入相关波门则对航迹进行删除，对未被使用的点迹重新进行起始，首先采用地波超视距雷达中常用的PDA进行跟踪，只采用地波模型时目标结果如图5所示，产生的多径量测生成了部分虚假航迹。
采用本发明方法的基于多探测模式的地波多径数据关联方法(Multi-detectiondata association methods for ground wave radar multipath data,MD-SRM)对目标跟踪结果如图6所示，本发明方法考虑了多条路径，利用所有目标量测，可以抑制虚假航迹的生成。
单次仿真具有偶然性，进行200次蒙特卡洛仿真对比本发明的MD-SGM方法与PDA方法的估计精度，采取均方根误差(RMSE)作为衡量精度性能的指标。航迹跟踪精度为目标航迹经过滤波后得到的状态信息与实际状态信息各维度之间的均方误差。对于第m次仿真，第t个目标的第v维的滤波误差为：
式中Xt(k)为第t个目标在第k时刻笛卡尔坐标系下的真实状态向量，Xt(k|k)为第t个目标在第k时刻笛卡尔坐标系下的目标状态向量；
则对于M次蒙特卡洛仿真可得第t个目标第v维的航迹跟踪精度(RMSE)为：
式中TraAcct,v为第t个目标第v维的航迹跟踪精度。
目标1与目标2采用PDA方法与MD-GM方法的RMSE结果如图7至图10所示。经对比，本发明方法与PDA方法相比，在X轴、Y轴距离上均有提升，滤波收敛后精度提升在10％以上。提升了目标航迹的估计精度。
具体实施例：
1)建立地波多径模型如图2所示，根据图2中的几何关系得到目标状态由笛卡尔坐标系转换为各途径雷达坐标系的转换公式；
2)根据多径数据位置关系对多径路径进行删除；
3)将多探测模式扩展到UKF方法中，对目标状态采样并求解权值，根据公式(1)可以确定共得到2ρ+1个采样点迹，对应的权值如公式(2)所示；
4)计算目标状态采样点的一步预测X(k|k-1)，计算误差协方差的一步预测P(k|k-1)；
5)基于k时刻目标状态X(k|k)根据式(2)至(6)对第i种传播模式进行预测，得到预测值
6)在多探测模式下，根据不同传播模式建立多个传播模式波门。假设在第k个扫描周期存在qmax种传播模式。所有传播模式的波门的并集即为总波门，假设存在三个传播途径形成三个波门，总波门范围内落入4个量测[z1(k),z2(k),z3(k),z4(k)]。如图4所示，qmax＝3，以表示选择的量测个数，最多选择量测个数/> 表示在m(k)个量测点迹中选择/>个的组合数量，在多探测模式下可分为以下四种情况：
时，此时量测组个数/>即不存在量测与目标关联；
时，此时量测组个数/>为[q1＝z1(k)]，[q1＝z2(k)]，[q2＝z2(k)]，[q2＝z3(k)]，[q3＝z4(k)]；
时，此时量测组个数/>为[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)]；
时，此时量测组个数/>为[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)]。当假设qmax＝3时，q1表示三种传播模式中的第一种，q2表示三种传播模式中的第二种，q3表示三种传播模式中的第三种。
7)进行如下判断：如预测值大于RGM+3σR，则删减第种传播模式预测值，否则保留。
8)对所有未删减的传播模式以预测值为中心，建立波门如式(10)，并对所有的波门求并集如，得到总波门如式(11)，并将总波门范围内的量测作为备选量测构造多探测量测组。
9)根据式(7)-(9)可得地波多径传播模式下的量测组状态采样点的一步预测由/> 可算得新息协方差/>
10)计算测量和状态向量的交互协方差和计算卡尔曼增益更新目标状态为更新误差协方差为/>
11)根据多探测量测组构造相关事件此时可得量测组/>源于目标的条件概率为/>此时目标被检测到的总概率为/>
12)对每个相关事件中的量测组进行状态更新与误差协方差更新，如式(22)和(23)所示，得到与/>
13)对所有相关事件中的量测组根据条件概率对目标状态和误差协方差进行更新，得到X(k|k)和P(k|k)。
虽然在本文中参照了特定的实施方式来描述本发明，但是应该理解的是，这些实施例仅仅是本发明的原理和应用的示例。因此应该理解的是，可以对示例性的实施例进行许多修改，并且可以设计出其他的布置，只要不偏离所附权利要求所限定的本发明的精神和范围。应该理解的是，可以通过不同于原始权利要求所描述的方式来结合不同的从属权利要求和本文中所述的特征。还可以理解的是，结合单独实施例所描述的特征可以使用在其它所述实施例中。
Claims (1)
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时刻的误差协方差矩阵，实现地波多径目标跟踪；
步骤一中第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个采样状态点对应的权值为：
式中Wi m代表第i个采样状态点的均值权值，Wi c代表第i个采样状态点的协方差权值；β为预设常数；
步骤一中，第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)：
式中ΔXi(k|k-1)为第i个采样状态点的预测值χi(k|k-1)与目标状态向量预测值X(k|k-1)的差值：
ΔXi(k|k-1)＝χi(k|k-1)-X(k|k-1) (6)，
Q(k)为第k时刻过程噪声矩阵；
将第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为雷达坐标系下的径向速度；
式中RGG为地波途径射距，RE-SG为E层天地波途径射距，RF-SG为F层天地波途径射距，RE-SS为E层天波途径射距，REF-SS为EF天波途径射距，hE为E层高度，hF为F层高度，DGG为地波途径径向速度，DE-SG为E层天地波途径径向速度，DF-SG为F层天地波途径径向速度，DE-SS为E层天波途径径向速度，DEF-SS为EF天波途径径向速度；
步骤二中按照上述方法由第k时刻目标状态向量预测值X(k|k-1)得到第q种传播模式下的第k时刻量测一步预测值传播模式最多包括五种，分别是地波途径传播模式、E层天地波途径传播模式、F层天地波途径传播模式、E层天波途径和EF天波途径传播模式；
表示预测值/>中的射距：
式中RGM为地波超视距雷达最大威力范围，为第q种传播模式下射距的测量误差均方差；
步骤三中，第k时刻第q种传播模式波门Gq(k)为：
式中z(k)为第k时刻的目标量测值，Sq(k)为第k时刻第q种传播模式的新息协方差；γ为波门范围；
则第k时刻总波门G(k)：
式中qmax表示第k时刻扫描周期存在的传播模式数量；
步骤四中，假设第k时刻存在三个传播模式，qmax＝3，形成三个波门G1，G2，G3；假设三个波门形成的总波门范围内落入m(k)个量测值，m(k)＝4；以表示选择的量测值个数，选择的最多量测值个数/>以/>表示在m(k)个量测值中选择/>个量测值的组合数量，并以量测组/>表示选择/>个量测值时的第/>个量测值组合；
构造第k时刻相关事件
则得到量测组源于目标的条件概率/>的表达式：
式中p表示概率函数，Zk表示到第k时刻被确认量测值的累积集合；
根据贝叶斯公式，将公式(13)变形为：
式中c为归一化系数，为概率密度函数，为相关事件的条件概率函数；
步骤五中，每个量测组第k时刻第i个采样点的采样状态一步预测值为：
式中表示第k时刻选择的/>个量测值中第/>个量测值对应的传播模式的量测转换公式；
则相关事件下量测组的第k时刻采样状态一步预测值为：
对应的量测组新息协方差矩阵为
式中为采样点的采样状态一步预测值/>与量测组的第k时刻采样状态一步预测值/>的差值：
为量测组噪声协方差矩阵：
表示选择的/>个量测值中第/>个量测值对应的传播模式的量测噪声协方差矩阵；步骤五中，量测组的卡尔曼增益/>为：
式中
更新得到量测组的第k时刻目标状态向量
同时更新得到量测组的第k时刻误差协方差矩阵
步骤五中，将量测组源于目标的条件概率分为概率密度函数和相关事件的条件概率函数/>进行计算：/>
式中PG表示波门概率，N(·)表示似然函数，为量测组新息；
式中u(·)为虚假量测值的概率质量函数，为目标在/>个传播模式下都能够被探测到的概率，PG为波门数，PD为目标被检测到的总概率：
步骤六中，更新得到的第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)。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310725543.9A CN116736287B (zh) | 2023-06-19 | 2023-06-19 | 基于多探测模式的地波多径目标跟踪方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310725543.9A CN116736287B (zh) | 2023-06-19 | 2023-06-19 | 基于多探测模式的地波多径目标跟踪方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116736287A CN116736287A (zh) | 2023-09-12 |
| CN116736287B true CN116736287B (zh) | 2024-04-05 |
Family
ID=87918213
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310725543.9A Active CN116736287B (zh) | 2023-06-19 | 2023-06-19 | 基于多探测模式的地波多径目标跟踪方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116736287B (zh) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107064901A (zh) * | 2017-04-27 | 2017-08-18 | 哈尔滨工业大学 | 一种舰载高频地波超视距雷达的目标方位估计方法 |
| CN107561503A (zh) * | 2017-08-28 | 2018-01-09 | 哈尔滨工业大学 | 一种基于多重渐消因子的自适应目标跟踪滤波方法 |
| CN108919253A (zh) * | 2018-05-17 | 2018-11-30 | 哈尔滨工业大学 | 一种基于帝国竞争算法的天波雷达多路径数据关联方法 |
| CN114942428A (zh) * | 2022-05-27 | 2022-08-26 | 哈尔滨工业大学 | 多目标跟踪的混合体制雷达数据关联算法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105607053B (zh) * | 2015-09-09 | 2018-05-04 | 湖北中南鹏力海洋探测系统工程有限公司 | 一种浮标式高频地波雷达系统 |
-
2023
- 2023-06-19 CN CN202310725543.9A patent/CN116736287B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107064901A (zh) * | 2017-04-27 | 2017-08-18 | 哈尔滨工业大学 | 一种舰载高频地波超视距雷达的目标方位估计方法 |
| CN107561503A (zh) * | 2017-08-28 | 2018-01-09 | 哈尔滨工业大学 | 一种基于多重渐消因子的自适应目标跟踪滤波方法 |
| CN108919253A (zh) * | 2018-05-17 | 2018-11-30 | 哈尔滨工业大学 | 一种基于帝国竞争算法的天波雷达多路径数据关联方法 |
| CN114942428A (zh) * | 2022-05-27 | 2022-08-26 | 哈尔滨工业大学 | 多目标跟踪的混合体制雷达数据关联算法 |
Non-Patent Citations (6)
| Title |
|---|
| A Method of Track Matching Based on Multipath Echoes in High-Frequency Surface Wave Radar;Jiazhi Zhang et al.;《IEEE Antennas and Wireless Propagation Letters》;20181031;第17卷(第10期);1852-1855 * |
| Wei Yinsheng et al..Experimental Analysis of a HF Hybrid Sky-Surface Wave Radar.《IEEE AEROSPACE AND ELECTRONIC SYSTEM MAGAZINE》.第33卷(第3期),32-40. * |
| 一种天波超视距雷达分级Hough变换航迹起始方法;金术玲;梁彦;潘泉;程咏梅;;电子与信息学报;20080815(第08期);194-198 * |
| 基于天波超视距雷达的多径数据关联算法研究与实现;王伟;《中国优秀硕士学位论文全文数据库 信息科技辑》;20190115;19-29 * |
| 天发舰收高频雷达数据关联算法研究;徐龙源;《万方数据知识服务平台》;20221114;52-53 * |
| 高频地波雷达机动多目标的信号处理方法;位寅生;边威;许荣庆;;哈尔滨工业大学学报;20070515(第05期);79-82 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116736287A (zh) | 2023-09-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109901153B (zh) | 基于信息熵权和最近邻域数据关联的目标航迹优化方法 | |
| CN108303692B (zh) | 一种解多普勒模糊的多目标跟踪方法 | |
| CN113064155B (zh) | 一种空中雷达多目标跟踪下航迹关联的优化方法 | |
| CN108734725B (zh) | 基于高斯过程的概率数据关联滤波扩展目标跟踪方法 | |
| CN112613532B (zh) | 基于雷达与循环神经网络补全红外融合的动目标跟踪方法 | |
| CN110146873B (zh) | 分布式非相参雷达的目标位置和速度的估计方法 | |
| KR101628154B1 (ko) | 수신 신호 세기를 이용한 다중 표적 추적 방법 | |
| CN114942428B (zh) | 多目标跟踪的混合体制雷达数据关联算法 | |
| CN112114308B (zh) | 一种扇扫雷达空时联合目标跟踪方法 | |
| CN110738275A (zh) | 基于ut-phd的多传感器序贯融合跟踪方法 | |
| CN116736287B (zh) | 基于多探测模式的地波多径目标跟踪方法 | |
| CN116679296A (zh) | 一种基于主辅复合滤波的远程高速目标雷达跟踪方法 | |
| CN111679270B (zh) | 一种反射点不确定场景下多路径融合目标检测算法 | |
| CN115204212A (zh) | 一种基于stm-pmbm滤波算法的多目标跟踪方法 | |
| Wang et al. | Research on multi-maneuvering target tracking JPDA algorithm | |
| CN113866754A (zh) | 基于高斯分布波门的运动目标航迹关联方法 | |
| CN114488104B (zh) | 基于交互一致性的天波超视距雷达目标跟踪方法 | |
| Yang et al. | Altitude estimation using multipath with a two-dimensional radar over spherical earth | |
| CN102426357B (zh) | 一种具有图像确认的多目标跟踪方法 | |
| CN112379366A (zh) | 一种基于粒子群滤波的雷达目标跟踪算法 | |
| CN116794646B (zh) | 基于变分贝叶斯推理的混合体制高频雷达目标跟踪方法 | |
| CN112285697A (zh) | 一种多传感器多目标空时偏差校准与融合方法 | |
| CN113514809B (zh) | 矢量脱靶量参数的测量方法、装置、电子设备和存储介质 | |
| CN114114170B (zh) | 一种雷达信号处理机目标测量的数学建模方法 | |
| CN114488116B (zh) | 一种基于两部两坐标外源雷达系统的3d目标跟踪方法 |
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 |