CN112114308B - 一种扇扫雷达空时联合目标跟踪方法 - Google Patents

一种扇扫雷达空时联合目标跟踪方法 Download PDF

Info

Publication number
CN112114308B
CN112114308B CN201910535351.5A CN201910535351A CN112114308B CN 112114308 B CN112114308 B CN 112114308B CN 201910535351 A CN201910535351 A CN 201910535351A CN 112114308 B CN112114308 B CN 112114308B
Authority
CN
China
Prior art keywords
target
track
state
radar
time
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
CN201910535351.5A
Other languages
English (en)
Other versions
CN112114308A (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.)
Harbin Institute of Technology
Original Assignee
Harbin Institute of 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 Harbin Institute of Technology filed Critical Harbin Institute of Technology
Priority to CN201910535351.5A priority Critical patent/CN112114308B/zh
Publication of CN112114308A publication Critical patent/CN112114308A/zh
Application granted granted Critical
Publication of CN112114308B publication Critical patent/CN112114308B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Systems 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/66Radar-tracking systems; Analogous systems
    • G01S13/72Radar-tracking systems; Analogous systems for two-dimensional tracking, e.g. combination of angle and range tracking, track-while-scan radar
    • G01S13/723Radar-tracking systems; Analogous systems for two-dimensional tracking, e.g. combination of angle and range tracking, track-while-scan radar by using numerical data
    • G01S13/726Multiple target tracking

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.不同方位上的目标预测时间(得到某个目标点迹后,预测同一目标下次出现的时间)不同,同一目标相邻点迹之间采样时间间隔不同,使得滤波与预测不能在固定重访间隔的条件下进行。
3.目标与雷达本身之间存在相对运动,同时存在噪声等因素影响,因此目标下一次被照射的时间存在不确定性。
传统多目标跟踪方法中,天线对目标的重访间隔往往是确定且已知的,这对扇扫雷达显然是不适用的。目前已有一些研究成果对于扇扫雷达目标重访间隔进行粗略计算,但均是建立在雷达天线方位角速度大小不变的前提下,由于扇扫雷达天线扫描方向会发生改变,显然这种前提是不合理的,存在较大的误差。
发明内容
(一)要解决的技术问题
本发明要解决的技术问题是解决传统目标跟踪方法对于扇扫雷达存在模型不匹配、检测偏差大的问题。
(二)技术方案
为了解决上述技术问题,本发明提供了一种扇扫雷达空时联合目标跟踪方法,对每个扇扫周期依次执行如下步骤:
S1、从观测雷达处获得当前周期的量测信息;
S2、依据上一周期更新后的跟踪结果判断当前周期是否存在暂时航迹和/或真实航迹;若存在,则顺序执行步骤S3,若不存在,则跳转执行步骤S7;
S3、对当前存在的所有航迹进行一步预测,根据其上一周期的航迹存在性概率得到当前周期的航迹存在性概率预测值,并根据扇扫雷达空时联合系统模型得到航迹目标状态的一步预测值和量测预测值;
S4、对每一条航迹,依据所述量测预测值建立相关波门,从步骤 S1中得到的量测信息中获取与航迹相关的点迹;
S5、对每一条航迹,利用与该航迹相关的点迹采用空时联合不敏卡尔曼滤波器进行滤波,获得多个状态估计值;
S6、对于每一条航迹,将其对应的多个状态估计值及状态的一步预测值进行概率互联,获得目标最终状态估计值及航迹存在性概率,更新为当前周期跟踪结果,依据航迹存在性概率判断航迹状态,输出确认航迹并删除终结航迹;
S7、对所有未与任一航迹相关的点迹进行滤波初始化,获得航迹目标状态值并设置初始航迹存在性概率,更新为当前周期跟踪结果。
优选地,所述步骤S7中初始化后,各个航迹初始均为暂时航迹,所述步骤S6中依据航迹存在性概率判断航迹状态时,若一暂时航迹的航迹存在性概率大于预先设置的确认门限tc,则该航迹被确认为真实航迹,并且保持确认状态直到被终结;若一暂时航迹或真实航迹的航迹存在性概率跌至小于预先设置的航迹终结门限tt,则该航迹被终结。
优选地,所述步骤S3中的扇扫雷达空时联合系统模型为线性加减速扇扫模式的雷达跟踪系统模型或三角函数加减速扇扫模式的雷达跟踪系统模型。
优选地,所述线性加减速扇扫模式的雷达跟踪系统模型中,天线的扇扫角度范围为[-β,β],匀速阶段天线扫描角速度为α,加速阶段天线加速度为
Figure GDA0003685640310000031
天线加速阶段扫描的角度为
Figure GDA0003685640310000032
天线扫描一周时间为Tradar;设当前周期为第k周期,目标第k-1次扫描到目标的时间为tk-1,第k周期内目标被扫到时的目标状态向量为
Figure GDA0003685640310000033
方位角为A,目标方位角与目标状态向量之间关系为
Figure GDA0003685640310000034
目标状态方程组的表达式为:
Figure GDA0003685640310000035
其中,F(k)是状态转移矩阵;v(k)是过程噪声向量;Γ(k)是噪声分布矩阵;g(tk-1,x(k))表示目标重访间隔与目标状态的非线性关系:
Figure GDA0003685640310000036
其中,Tk为目标重访间隔;f与目标扫描方向相关,若第k周期目标顺时针扫描则f=-1,若目标逆时针扫描则f=1。
优选地,所述三角函数加减速扇扫模式的雷达跟踪系统模型中,天线扫描角速度加减速以三角函数形式进行,天线加速阶段或减速阶段扫描的角度范围为α/2,加速阶段天线扫描角速度从0增加到α,加速阶段步长为t=1,速度曲线采用α(t)=α/2(1-cos(πt)),天线扫描一周时间为Tradar;设当前周期为第k周期,目标第k-1次扫描到目标的时间为tk-1,第k周期内目标被扫到时的目标状态向量为
Figure GDA0003685640310000041
方位角为A,目标方位角与目标状态向量之间关系为
Figure GDA0003685640310000042
当A|<β-Δa时,目标状态方程组为:
Figure GDA0003685640310000043
当A|>β-Δa时,目标状态方程组为:
Figure GDA0003685640310000044
其中,F(k)是状态转移矩阵;v(k)是过程噪声向量;Γ(k)是噪声分布矩阵;目标重访间隔Tk=tk-tk-1,Δt为天线在第k周期内天线扫过目标到距离目标较近的扫描边界之间的角度所用时间:
Figure GDA0003685640310000045
优选地,所述线性加减速扇扫模式的雷达跟踪系统模型或三角函数加减速扇扫模式的雷达跟踪系统模型中,量测方程为:
Figure GDA0003685640310000046
其中,z(k)为目标的量测向量,包括目标相对所述观测雷达坐标系原点的距离量测rk、方位角量测θk;hk(x(k))表示目标状态与量测之间的数学关系;W(k)为零均值的高斯白噪声序列,表示k时刻量测噪声,包括
Figure GDA0003685640310000047
两个分量,分别为距离、方位角量测噪声;
量测噪声协方差矩阵为:
Figure GDA0003685640310000048
其中,Rk,rr、Rk,θθ分别表示k时刻各量测噪声分量的自协方差,其值分别为
Figure GDA0003685640310000049
Rk,rθ分别表示k时刻各量测噪声分量的互协方差,假设各量测之间是互不相关的,因此各分量互协方差均为0。
优选地,所述步骤S3中根据其上一周期的航迹存在性概率得到当前周期的航迹存在性概率预测值时,通过马尔科夫链1阶模型表述航迹存在或不存在两种状态的转移概率:
Figure GDA0003685640310000051
其中,P11表示上一周期航迹存在,下一周期航迹依然存在的概率; P12表示上一周期航迹存在,下一周期航迹不存在的概率;P21表示上一周期航迹不存在,下一周期航迹存在的概率;P22表示上一周期航迹不存在,下一周期航迹依然不存在的概率;
已知k-1周期航迹τ存在的概率为
Figure GDA0003685640310000052
表示k-1周期航迹τ存在;定义与航迹存在性相关的两种状态:
Figure GDA0003685640310000053
表示目标在第k周期航迹τ存在,
Figure GDA0003685640310000054
表示目标在第k周期航迹τ不存在;
对航迹存在性概率进行预测,得到k周期该航迹存在性概率的预测值为:
Figure GDA0003685640310000055
Figure GDA0003685640310000056
其中,Zτ,k-1表示直到k-1时刻落入航迹τ波门的所有量测的集合,
Figure GDA0003685640310000057
Zτ(k)表示k时刻落入航迹τ相关波门内的
Figure GDA0003685640310000058
个量测的集合,
Figure GDA0003685640310000059
优选地,所述步骤S3中根据扇扫雷达空时联合系统模型得到状态的一步预测值时,选择均值周围的采样点δ作为非线性变换的输入,对输出结果求取统计特性,得到目标状态的一步预测:
设第k-1周期中目标τ的状态估计值为
Figure GDA00036856403100000510
协方差为 Pτ(k-1|k-1);依据该状态和协方差,产生一组长度为2L+1的采样点δ,各采样点δ均匀分布在第k-1周期状态估计值
Figure GDA00036856403100000512
的附近,其中L 是状态向量
Figure GDA00036856403100000511
的维数:
Figure GDA0003685640310000061
其中,
Figure GDA0003685640310000062
是(L+λ)Pτ(k-1|k-1)的矩阵平方根的第i 列;
将上述2L+1个采样点分别作为扇扫雷达空时联合系统模型中目标状态方程组的输入,采用迭代法求解得到2L+1组解,获得2L+1个重访间隔及状态预测值;
目标的重访间隔是2L+1个采样点重访间隔的加权总和:
Figure GDA0003685640310000063
对2L+1个状态预测值进行加权,获得目标最终状态预测值及状态预测协方差:
Figure GDA0003685640310000064
Figure GDA0003685640310000065
其中,
Figure GDA0003685640310000066
Γi,k是依据
Figure GDA0003685640310000067
得到的噪声分布矩阵,状态和协方差的加权权重分别是:
Figure GDA0003685640310000068
α和κ控制采样点δ的传播;β与x的分布有关;
根据扇扫雷达空时联合系统模型中的量测方程,预测测量的采样点是:
Figure GDA0003685640310000071
量测的预测值和相应的协方差分别为:
Figure GDA0003685640310000072
Figure GDA0003685640310000073
其中,
Figure GDA0003685640310000074
状态预测值及量测预测值的互协方差矩阵是:
Figure GDA0003685640310000075
优选地,所述步骤S3中,将2L+1个采样点分别作为扇扫雷达空时联合系统模型中目标状态方程组的输入,采用迭代法求解得到2L+1 组解时,包括如下步骤:
以第i个采样点
Figure GDA0003685640310000076
为输入,若在第k-1周期里目标状态更新时间为tk-1,则天线扫过目标到达扫描边界时间为T1=Tradar·(k-1)-tk-1,迭代法求解状态方程组中目标重访间隔的初始值为:
T′=2·T1=2·(Tradar·(k-1)-tk-1);
以T′为目标重访间隔,依据目标状态方程组中的第二个方程获得迭代初始的状态预测值x′k|k-1为:
Figure GDA0003685640310000077
其中,状态转移矩阵F′(k)依据T′得到:
Figure GDA0003685640310000078
计算第k周期目标被天线扫描到时的目标方位角:
Figure GDA0003685640310000081
依据目标方位角,计算天线在第k周期里从天线扫描边界到扫过目标的时间间隔T2
存在时间差为ΔT=T′-(T1+T2),设置时间差的最大阈值为σ;若 |ΔT|<σ,则认为T′与目标真实的重访间隔差距可忽略不计,重访间隔
Figure GDA0003685640310000082
目标状态预测值为
Figure GDA0003685640310000083
否则对T′修正为T′=T′-0.5·ΔT,并返回至以T′为目标重访间隔,依据目标状态方程组中的第二个方程获得迭代初始的状态预测值x′k|k-1步骤,直到|ΔT|<σ;最终获得以第i个采样点为输入获得的目标状态的一步预测值
Figure GDA0003685640310000084
及重访间隔
Figure GDA0003685640310000085
其中,计算天线在第k周期里从天线扫描边界到扫过目标的时间间隔T2时,T2的计算方式与目标方位角和天线扫描方向有关,对于所述线性加减速扇扫模式的雷达跟踪系统模型,T2的表达式为:
Figure GDA0003685640310000086
对于三角函数加减速扇扫模式的雷达跟踪系统模型,若目标位于匀速区即|A|<β-α/2,则T2=1+(β-Δa+f·A)/α;
若目标位于加速或减速区,则通过求解
Figure GDA0003685640310000087
求取T2,若目标位于加速区则T2=Δt,目标位于减速区则T2=Tradar-Δt。
优选地,所述步骤S5中利用与该航迹相关的点迹采用空时联合不敏卡尔曼滤波器进行滤波时,设一航迹τ有mk个相关的点迹,分别利用mk个点迹应用空时联合不敏卡尔曼滤波器进行滤波,获得该航迹τ相关的mk个状态估计值;
其中,滤波包括获得航迹目标状态的一步预测值和量测预测值后,计算卡尔曼增益:
Figure GDA0003685640310000088
依据一点迹i得到的状态更新值是预测状态与被卡尔曼增益加权的新息之和,得到:
xτ,i(k|k)=xτ(k|k-1)+K(k)[zi(k)-zτ(k|k-1)];
协方差更新值为:
Figure GDA0003685640310000091
(三)有益效果
本发明的上述技术方案具有如下优点:本发明提供了一种扇扫雷达空时联合目标跟踪方法,该方法采用递推循环的方式实现多目标跟踪,采用增加了描述时间的方程对状态方程进行约束的扇扫雷达空时联合系统模型对扇扫雷达跟踪系统进行描述,并提出了相应的滤波方法,有效解决了传统目标跟踪方法中系统模型与扇扫雷达存在模型不匹配的问题,提高了目标跟踪在准确性以及真实航迹确认速度。
附图说明
图1示出了本发明实施例中扇扫雷达空时联合目标跟踪方法 (STJ-LMIPDA方法)步骤示意图;
图2示出了本发明实施例中迭代法求解目标状态方程组解的步骤示意图;
图3示出了仿真实验中构造的笛卡尔坐标系下三个目标沿直线匀速运动的目标轨迹;
图4示出了线性加减速扇扫机制下,对目标2运用STJ-LMIPDA 方法与区分奇偶周期以2倍天线扫描周期为重访间隔运用传统 LMIPDA方法进行比较,得到的位置分量估计的均方根误差结果;
图5示出线性加减速扇扫机制下两种方法对目标2速度估计均方根误差对比结果;
图6示出了三角函数加减速扇扫机制下两种方法对目标2位置分量估计的均方根误差结果;
图7示出了三角函数加减速扇扫机制下两种方法对目标2速度分量估计的均方根误差结果;
图8显示了线性加减速扇扫机制下两种方法真实航迹随扫描次数增加的关系图;
图9显示了线性加减速扇扫机制下两种方法虚假航迹随扫描次数增加的关系图;
图10显示了三角函数加减速扇扫机制下两种方法真实航迹随扫描次数增加的关系图;
图11显示了三角函数加减速扇扫机制下两种方法虚假航迹随扫描次数增加的关系图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明的一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例,都属于本发明保护的范围。
如图1所示,本发明实施例提供的一种扇扫雷达空时联合目标跟踪方法,是递推循环的过程,设当前周期为第k周期,对每个扇扫周期依次执行如下步骤:
S1、从观测雷达处获得当前周期的量测信息;
S2、依据上一周期更新后的跟踪结果判断当前周期是否存在暂时航迹和/或真实航迹;若存在,则顺序执行步骤S3,若不存在,则跳转执行步骤S7;
S3、对当前存在的所有航迹,即所有暂时航迹和/或真实航迹,进行一步预测,根据其上一周期的航迹存在性概率得到当前周期的航迹存在性概率预测值,并根据扇扫雷达空时联合系统模型得到航迹目标状态的一步预测值和量测预测值;
S4、对每一条航迹,依据所述量测预测值建立相关波门,从步骤 S1中得到的量测信息中获取与航迹相关的点迹;
S5、对每一条航迹,利用与该航迹相关的点迹采用空时联合不敏卡尔曼滤波器(STJ-UKF)进行滤波,获得多个状态估计值;
S6、对于每一条航迹,将其对应的多个状态估计值及状态的一步预测值进行概率互联,获得目标最终状态估计值及航迹存在性概率,更新为当前周期跟踪结果,依据航迹存在性概率判断航迹状态,输出确认航迹并从内存里删除终结航迹;
S7、对所有未与任一航迹相关的点迹进行滤波初始化,获得航迹目标状态值并设置初始航迹存在性概率,更新为当前周期跟踪结果。
由于航迹包括暂时航迹和真实航迹,优选地,在(上一周期)所述步骤S7中进行初始化后,各个航迹初始均为暂时航迹,所述步骤S6 中依据航迹存在性概率判断航迹状态时,预先设置门限,若一暂时航迹的航迹存在性概率大于预先设置的确认门限tc,则该航迹被确认为真实航迹,并且保持确认状态直到被终结;若一暂时航迹或真实航迹的航迹存在性概率跌至小于预先设置的航迹终结门限tt,则该航迹被终结,对判断终结的航迹从内存中删除,减少占用空间。
优选地,所述步骤S3中根据其上一周期的航迹存在性概率得到当前周期的航迹存在性概率预测值时,通过马尔科夫链1阶模型表述航迹存在或航迹不存在这两种状态的转移概率:
Figure GDA0003685640310000111
其中,P11表示上一周期航迹存在,下一周期航迹依然存在的概率; P12表示上一周期航迹存在,下一周期航迹不存在的概率;P21表示上一周期航迹不存在,下一周期航迹存在的概率;P22表示上一周期航迹不存在,下一周期航迹依然不存在的概率;P11、P12、P21和P22的具体数值可由工程师根据实际情况进行设置。
已知k-1周期航迹τ存在的概率为
Figure GDA0003685640310000112
表示k-1周期航迹τ存在;定义与航迹存在性相关的两种状态:
Figure GDA0003685640310000121
表示目标在第k周期航迹τ存在,
Figure GDA0003685640310000122
表示目标在第k周期航迹τ不存在;
对航迹存在性概率进行预测,得到k周期该航迹存在性概率的预测值为:
Figure GDA0003685640310000123
Figure GDA0003685640310000124
其中,Zτ,k-1表示直到k-1时刻落入航迹τ波门的所有量测的集合,
Figure GDA0003685640310000125
Zτ(k)表示k时刻落入航迹τ相关波门内的
Figure GDA0003685640310000126
个量测的集合,
Figure GDA0003685640310000127
考虑到扇扫雷达的通用性,本发明针对线性加减速扇扫雷达及三角函数加减速扇扫雷达进行了建模。优选地,所述步骤S3中的扇扫雷达空时联合系统模型为线性加减速扇扫模式的雷达跟踪系统模型或三角函数加减速扇扫模式的雷达跟踪系统模型。
优选地,所述线性加减速扇扫模式的雷达系统模型中,天线的扇扫角度范围为[-β,β],匀速阶段天线扫描角速度为α,加速阶段天线加速度为
Figure GDA0003685640310000128
天线加速阶段扫描的角度为
Figure GDA0003685640310000129
天线扫描一周时间为Tradar;设当前周期为第k周期,目标第k-1次扫描到目标的时间为tk-1,第k周期内目标被扫到时的目标状态向量为
Figure GDA00036856403100001210
方位角为A,目标方位角与目标状态向量之间关系为
Figure GDA00036856403100001211
目标状态方程组的表达式为:
Figure GDA00036856403100001212
其中,F(k)是状态转移矩阵;v(k)是过程噪声向量;Γ(k)是噪声分布矩阵;g(tk-1,x(k))表示目标重访间隔与目标状态的非线性关系:
Figure GDA0003685640310000131
其中,Tk为目标重访间隔;
Figure GDA0003685640310000132
f与目标扫描方向相关,若第k周期目标顺时针扫描则f=-1,若目标逆时针扫描则f=1。
优选地,所述三角函数加减速扇扫模式的雷达跟踪系统模型中,天线扫描角速度加减速以三角函数形式进行,天线加速阶段或减速阶段扫描的角度范围为α/2,加速阶段天线扫描角速度从0增加到α,加速阶段步长为t=1,速度曲线采用α(t)=α/2(1-cos(πt)),天线扫描一周时间为Tradar;设当前周期为第k周期,目标第k-1次扫描到目标的时间为tk-1,第k周期内目标被扫到时的目标状态向量为
Figure GDA0003685640310000133
方位角为A,目标方位角与目标状态向量之间关系为
Figure GDA0003685640310000134
当A|<β-Δa时,目标状态方程组为:
Figure GDA0003685640310000135
当A|>β-Δa时,目标状态方程组为:
Figure GDA0003685640310000136
其中,F(k)是状态转移矩阵;v(k)是过程噪声向量;Γ(k)是噪声分布矩阵;目标重访间隔Tk=tk-tk-1,Δt为天线在第k周期内天线扫过目标到距离目标较近的扫描边界之间的角度所用时间:
Figure GDA0003685640310000137
优选地,线性加减速扇扫模式的雷达跟踪系统模型或三角函数加减速扇扫模式的雷达跟踪系统模型中,两种扇扫方式的量测方程是一致的,量测方程的表达式为:
Figure GDA0003685640310000141
其中,z(k)为目标的量测向量,包括目标相对所述观测雷达坐标系原点的距离量测rk、方位角量测θk;hk(x(k))表示目标状态与量测之间的数学关系;W(k)为零均值的高斯白噪声序列,表示k时刻量测噪声,包括
Figure GDA0003685640310000142
两个分量,分别为距离、方位角量测噪声;
量测噪声协方差矩阵为:
Figure GDA0003685640310000143
其中,Rk,rr、Rk,θθ分别表示k时刻各量测噪声分量的自协方差,其值分别为
Figure GDA0003685640310000144
Rk,rθ分别表示k时刻各量测噪声分量的互协方差,假设各量测之间是互不相关的,因此各分量互协方差均为0。
优选地,所述步骤S3中根据扇扫雷达空时联合系统模型得到状态的一步预测值时,依据UT变换,选择均值周围的采样点δ作为非线性变换的输入,对输出结果求取统计特性,得到目标状态的一步预测。
设第k-1周期中目标τ的状态估计值为
Figure GDA0003685640310000145
协方差为 Pτ(k-1|k-1);依据该状态和协方差,产生一组长度为2L+1的采样点δ,各采样点δ均匀分布在第k-1周期,即上一周期,更新后的状态估计值
Figure GDA0003685640310000146
的附近,其中L是状态向量
Figure GDA0003685640310000147
的维数;
Figure GDA0003685640310000148
其中,
Figure GDA0003685640310000149
是(L+λ)Pτ(k-1|k-1)的矩阵平方根的第i 列;
将上述2L+1个采样点分别作为扇扫雷达空时联合系统模型中目标状态方程组的输入,采用迭代法求解得到2L+1组解,获得2L+1个重访间隔及状态预测值;
目标的重访间隔是2L+1个采样点重访间隔的加权总和:
Figure GDA0003685640310000151
对2L+1个状态预测值进行加权,获得目标最终状态预测值及状态预测协方差:
Figure GDA0003685640310000152
Figure GDA0003685640310000153
其中,
Figure GDA0003685640310000154
Γi,k是依据
Figure GDA0003685640310000155
得到的噪声分布矩阵,状态和协方差的加权权重分别是:
Figure GDA0003685640310000156
α和κ控制采样点δ的传播;β与x的分布有关;
根据扇扫雷达空时联合系统模型中的量测方程,预测测量的采样点是:
Figure GDA0003685640310000157
量测的预测值和相应的协方差分别为:
Figure GDA0003685640310000158
Figure GDA0003685640310000161
其中,
Figure GDA0003685640310000162
状态预测值及量测预测值的互协方差矩阵是:
Figure GDA0003685640310000163
优选地,如图2所示,所述步骤S3中,将2L+1个采样点分别作为扇扫雷达空时联合系统模型中目标状态方程组的输入,采用迭代法求解得到2L+1组解时,对于线性加减速扇扫模式的雷达跟踪系统模型,包括如下步骤:
1)以第i个采样点
Figure GDA0003685640310000164
为输入,若在第k-1周期里目标状态更新时间为tk-1,则天线扫过目标到达扫描边界时间为T1=Tradar·(k-1)-tk-1,迭代法求解状态方程组中目标重访间隔的初始值为:
T′=2·T1=2·(Tradar·(k-1)-tk-1);
2)以T′为目标重访间隔,依据目标状态方程组中的第二个方程获得迭代初始的状态预测值x′k|k-1为:
Figure GDA0003685640310000165
其中,状态转移矩阵F′(k)依据T′得到:
Figure GDA0003685640310000166
3)计算第k周期目标被天线扫描到时的目标方位角:
Figure GDA0003685640310000167
其中,x′k|k-1(1)、x′k|k-1(3)分别表示4x1维向量x′k|k-1的第1个值和第3 个值;
4)依据目标方位角,计算天线在第k周期里从天线扫描边界到扫过目标的时间间隔T2,T2的计算方式与目标方位角和天线扫描方向有关,表达式为:
Figure GDA0003685640310000171
5)存在时间差为ΔT=T′-(T1+T2),设置时间差的最大阈值为σ,若 |ΔT|<σ,则认为T′与目标真实的重访间隔差距非常小可以忽略不计,重访间隔
Figure GDA0003685640310000172
目标状态预测值为
Figure GDA0003685640310000173
否则对T′修正为T′=T′-0.5·ΔT,并返回至以T′为目标重访间隔,依据目标状态方程组中的第二个方程获得迭代初始的状态预测值x′k|k-1步骤,即在修正T′后重新进行步骤2)到5),直到|ΔT|<σ;最终获得以第 i个采样点为输入获得的目标状态的一步预测值
Figure GDA0003685640310000174
及重访间隔
Figure GDA0003685640310000175
优选地,所述步骤S3中,将2L+1个采样点分别作为扇扫雷达空时联合系统模型中目标状态方程组的输入,采用迭代法求解得到2L+1 组解时,对于三角函数加减速扇扫模式的雷达跟踪系统模型,包括如下步骤:
1)以第i个采样点
Figure GDA0003685640310000176
为输入,若在第k-1周期里目标状态更新时间为tk-1,则天线扫过目标到达扫描边界时间为T1=Tradar·(k-1)-tk-1,目标重访间隔的初始值为:
T′=2·T1=2·(Tradar·(k-1)-tk-1);
2)以T′为目标重访间隔,依据状态方程组中的第二个方程获得迭代初始的状态预测值x′k|k-1,即为:
Figure GDA0003685640310000177
其中,状态转移矩阵F′(k)依据T′得到,即:
Figure GDA0003685640310000181
3)计算第k周期目标被天线扫描到时目标方位角:
Figure GDA0003685640310000182
4)判断目标所在的区域,依据目标方位角,计算在天线k周期内,从天线扫描边缘扫到目标的时间T2
若目标位于匀速区即|A|<β-α/2,则T2=1+(β-Δa+f·A)/α;
若目标位于加速或减速区,求取T2需要求解
Figure GDA0003685640310000183
鉴于该方程是非线性的,同样可利用迭代法求解该方程的解,求解过程如下:
1、天线在加速或减速区域扫描的时间为1s,令迭代初始值T2′=1,由于T2′的值存在误差,存在角度误差
Figure GDA0003685640310000184
2、比较角度误差的阈值是σ1。它可能是一个非常小的值,可按需自行设定。如果角度误差小于阈值即|ε|<σ1,那么认为角度误差可以忽略,即Δt=Δt′。否则需要对Δt′进行调整,令
Figure GDA0003685640310000185
重新执行第1 步骤,计算角度误差并对时间Δt′进行调整,直至角度误差满足|ε|<σ1,若目标位于加速区那么T2=Δt,目标位于减速区T2=Tradar-Δt。
5)依据目标方位角A获得时间T2,求取时间差为ΔT=T′-(T1+T2),设置时间差的最大阈值为σ2。若|ΔT|<σ2则认为T′与目标真实的重访间隔差距非常小可以忽略不计,重访间隔
Figure GDA0003685640310000186
目标状态预测值为
Figure GDA0003685640310000187
否则对T′修正为T′=T′-0.5·ΔT,并返回至以T′为目标重访间隔,依据状态方程组中的第二个方程获得迭代初始的状态预测值 x′k|k-1步骤,即在修正T′后重新进行步骤2)到5),直到|ΔT|<σ2。至此我们得到第i个采样点为输入获得的目标状态的一步预测值
Figure GDA0003685640310000188
及重访间隔
Figure GDA0003685640310000191
优选地,步骤S4中对每一条航迹,依据所述量测预测值建立相关波门,从步骤S1中得到的量测信息中获取与航迹相关的观测点时,对 k时刻雷达接收所有的点迹进行波门相关,相关波门依据被跟踪目标在当前时刻所在位置的预测值为中心,确立一个量测向量可能出现的范围的决策门限,所有落入该门限的量测值均考虑为候选的量测,而在门限值之外的量测则认为是杂波。相关波门通常用来滤除杂波,找到所有可能来自于目标的量测值。常用的相关波门有矩形波门、环形波门、椭圆(球)波门、极坐标系下的扇形波门等,这里我们将使用椭圆波门。
对于航迹对于传感器在第k周期内收到的任何点迹
Figure GDA0003685640310000192
若存在:
Figure GDA0003685640310000193
则认为该点迹落在目标τ的波门,那么存在一个假设即,点迹i为目标τ的量测值,并且
Figure GDA0003685640310000194
将用于更新目标τ的状态估计值。给定目标的波门概率γ。
一条航迹在k时刻落入波门内的点迹可能有多个,也可能没有。
优选地,若k时刻有mk个点迹落入航迹τ的相关波门,对于其中的第i个点迹zi(k)(这里i与前文中的i意义不同,表征落入相关波门内的某一点迹,前文中指代UT变换过程中某一采样点),利用点迹zi(k) 对航迹τ进行状态更新,得到该点迹对航迹τ的状态更新值
Figure GDA0003685640310000195
同时得到该点迹来自航迹τ所跟踪目标的似然函数
Figure GDA0003685640310000196
进一步优选地,所述步骤S5中对每一条航迹,利用与该航迹相关的点迹进行滤波过程,获得多个状态估计值时,设一航迹τ有
Figure GDA0003685640310000197
个相关的点迹,分别利用
Figure GDA0003685640310000198
个点迹应用空时联合不敏卡尔曼滤波器 (STJ-UKF)进行滤波,获得该航迹τ相关的
Figure GDA0003685640310000199
个状态估计值。STJ-UKF 方法前半部分内容即状态一步预测过程,这里不再重复描述。
滤波还包括获得航迹目标状态的一步预测值和量测预测值后,计算卡尔曼增益:
Figure GDA0003685640310000201
依据某一点迹i,即zi(k),得到的状态更新值是预测状态与被卡尔曼增益加权的新息之和,得到:
xτ,i(k|k)=xτ(k|k-1)+K(k)[zi(k)-zτ(k|k-1)];
相应的协方差更新值为:
Figure GDA0003685640310000202
优选地,步骤S6中对于每一条航迹,将其对应的多个状态估计值及状态的一步预测值进行概率互联时,如果点迹i落在目标τ的波门内,则表示点迹i与目标τ相关,点迹i可能是航迹τ的观测值,定义事件
Figure GDA0003685640310000203
为波门内的第i个点迹来自目标τ,其他点迹对该航迹为杂波。
点迹i来自航迹τ的先验概率为:
Figure GDA0003685640310000204
其中
Figure GDA0003685640310000205
分别表示目标被检测的概率及波门概率,可由工程师根据需要进行设定。
Figure GDA0003685640310000206
为航迹存在性概率的预测值,前文已经介绍。
Figure GDA0003685640310000207
N表示目标数量,ρk,i是目标的杂波密度,
Figure GDA0003685640310000208
是在zk(i)是目标τ的量测的假设下的预测测量概率密度函数,其由下式定义:
Figure GDA0003685640310000209
在多目标跟踪系统中,目标τ波门中的点迹i可能来自另一个目标。在线性化方法下,通过排除目标τ的所有可能来源来修改量测i所在位置的杂波密度。跟踪τ的验证门中的测量i的先验散射体测量密度,由Ωτ(i)表示,可以表示为:
Figure GDA0003685640310000211
结合综合概率数据互联算法(IPDA),目标存在的后验概率被修改为:
Figure GDA0003685640310000212
Figure GDA0003685640310000213
其中,
Figure GDA0003685640310000214
Figure GDA0003685640310000215
表示落入目标门的测量总数,量测i是目标τ(i=0表示该目标没有量测)的量测的后验概率由下式给出:
Figure GDA0003685640310000216
Figure GDA0003685640310000217
其中,
Figure GDA0003685640310000218
至此获得落入目标τ的波门内的各个量测属于该目标的后验概率,以及航迹τ存在的后验概率值。
假设有
Figure GDA0003685640310000219
个点迹落在目标τ的波门内,则存在
Figure GDA00036856403100002110
种假设:目标没有量测被探测到
Figure GDA00036856403100002111
个点迹均为杂波;这
Figure GDA00036856403100002112
个点迹分别是目标τ的量测值,其他
Figure GDA00036856403100002113
个点迹为杂波。后验跟踪状态估计是:
Figure GDA00036856403100002114
其中
Figure GDA00036856403100002115
是以点迹
Figure GDA00036856403100002116
作为目标的量测值对目标进行状态估计得到的估计值,i=0表示目标没有量测被探测,目标状态的预测值,即
Figure GDA0003685640310000221
状态估计协方差矩阵由下式给出:
Figure GDA0003685640310000222
目标τ的更新时间是
Figure GDA0003685640310000223
至此,目标τ的状态更新已经完成。
对所有航迹均进行以上航迹更新过程,依据航迹存在性概率对航迹进行管理。某些点迹可能未落入任何一条航迹的相关波门内,则需要进行航迹初始化。
考虑到未被用来更新航迹的点迹可能代表着出现新航迹,对其执行航迹初始化,可获得新航迹目标初始状态。优选地,步骤S7中对所有未与任一航迹相关的点迹,依据观测雷达当前周期和上一周期所获得的未与航迹相关的点迹进行滤波初始化时,利用未与航迹相关的点进行航迹初始化获得新航迹,以第一周期和第二周期的点迹为例,第一次扫描雷达得到n1个点迹,第二次扫描雷达得到n2个点迹,利用第一次扫描时刻的第i(1≤i≤n1)个点迹和第二次扫描时刻的第j(1≤j≤n2)个点迹获得第τ(1≤τ≤n1×n2)条航迹在笛卡尔坐标系下第二次扫描时刻的状态向量。
Figure GDA0003685640310000224
对应的状态协方差为:
Figure GDA0003685640310000231
其中下标1、2分别表示第一周期和第二周期对应的值。T表示雷达扫描周期。
Figure GDA0003685640310000232
Figure GDA0003685640310000233
是笛卡尔坐标下的量测沿x,y方向的位置量测信息,是通过无偏量测转换方法将雷达所获得量测(径向距离、方位角) 转换到笛卡尔坐标系下得到的转换位置量测,转换公式为:
Figure GDA0003685640310000234
其中rk,θk是从雷达获取的距离,方位角量测;
Figure GDA0003685640310000235
是转换后得到的沿x,y方向的笛卡尔坐标量测,
Figure GDA0003685640310000236
是转换后的量测向量;μθ是去偏系数,可通过方位角测量噪声方差
Figure GDA0003685640310000237
求得:
Figure GDA0003685640310000238
对应的协方差矩阵为:
Figure GDA0003685640310000239
其中各分量的自协方差及互协方差:
Figure GDA00036856403100002310
上标“c”代表与转换量测相关的向量、矩阵和函数。初始化的状态误差协方差为:
Figure GDA0003685640310000241
T可以通过两个量测的时间差得到。
航迹存在性概率初始值可由实际需要进行设定。
综上,本申请提供了一种扇扫雷达空时联合目标跟踪方法,该方法包括如下优势:
(1)针对扇扫雷达重访间隔不确定导致系统模型不匹配的问题,通过增加时间约束方程,分别针对线性加减速扇扫雷达与三角函数扇扫机制雷达建立了适用于扇扫雷达的状态方程组,准确描述两种扫描机制下目标状态随时间变化的关系。
(2)针对新建立的系统模型,研究了对应的空时联合状态估计滤波方法STJ-UKF,通过迭代法求取非线性状态方程组的解,同时得到目标状态预测值及重访间隔,该过程可以看作一个广义的非线性函数,因此可应用不敏变换计算其输出变量的统计特性,最后获得目标状态估计结果。
(3)将滤波方法应用于线性多目标联合概率密度互联(LMIPDA) 框架,研究了适应于的扇扫雷达的多目标跟踪方法STJ-LMIPDA方法
为验证本发明的效果,利用仿真数据进行蒙特卡洛实验。仿真实验监视区域内设置三个目标,仿真试验中的目标沿直线匀速运动,运动轨迹如图3所示。三个目标初始状态如表1所示。
表1 三 个目标的初始状态
Figure GDA0003685640310000242
进行100次蒙脱卡洛仿真,每次模拟运行有100次扫描。监视区域中的杂波密度为6.0×10-6/scan/m2且杂波满足均匀的泊松分布。目标初始位置和速度总结在表1 中给出。过程噪声设置为ds=0.01m/s2,范围和方位角测量标准偏差分别为σr=50m,σθ=1deg。天线匀速扫描阶段扫描角速度为α=π/6 rad/s,天线扫描范围为[-π/3,π/3],迭代过程时间差的阈值为σ2=0.001s。
线性加减速扇扫方式加速阶段天线角度加速度为
Figure GDA0003685640310000251
三角函数加减速扇扫方式加速时间为1s,角度误差的阈值为σ1=10-8rad。
将区分奇偶扫描周期,以两倍天线扫描周期作为重访间隔的传统 LMIPDA和本文所提方法进行比较用于评估STJ-LMIPDA的性能。
图4示出了线性加减速扇扫方式两种方法对目标2位置滤波结果对应的估计均方根误差对比,图5示出了线性加减速扇扫方式两种方法对目标2速度滤波结果对应的估计均方根误差对比。图6示出了三角函数加减速扇扫方式两种方法对目标2位置滤波结果对应的估计均方根误差对比,图7示出了三角函数加减速扇扫方式两种方法对目标2 速度滤波结果对应的估计均方根误差对比。从图4到图6中可以看到不管哪种扫描方式STJ-LMIPDA方法对目标估计精度有显著提高。这是由于通过对扇扫雷达系统模型的准确建立,通过迭代法获得目标的状态预测值及重访间隔较为精确,综合考虑了目标运动及目标方位角对重访间隔的影响,使得数据互联较为准确。在精准建模的基础上,状态滤波结果也较为准确。但本发明中示出的方法在提供状态估计结果的同时,参照图8及图10可以看出本发明所提方法还能够显著提高真实航迹确认速度;参照图9及图11可以看出除开始几个周期本文所提方法虚警稍高外,后续周期显著降低了目标跟踪过程的虚警率。
最后应说明的是:以上实施例仅用以说明本发明的技术方案,而非对其限制;尽管参照前述实施例对本发明进行了详细的说明,本领域的普通技术人员应当理解:其依然可以对前述各实施例所记载的技术方案进行修改,或者对其中部分技术特征进行等同替换;而这些修改或者替换,并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (7)

1.一种扇扫雷达空时联合目标跟踪方法,其特征在于,对每个扇扫周期依次执行如下步骤:
S1、从观测雷达处获得当前周期的量测信息;
S2、依据上一周期更新后的跟踪结果判断当前周期是否存在暂时航迹和/或真实航迹;若存在,则顺序执行步骤S3,若不存在,则跳转执行步骤S7;
S3、对当前存在的所有航迹进行一步预测,根据其上一周期的航迹存在性概率得到当前周期的航迹存在性概率预测值,并根据扇扫雷达空时联合系统模型得到航迹目标状态的一步预测值和量测预测值;
S4、对每一条航迹,依据所述量测预测值建立相关波门,从步骤S1中得到的量测信息中获取与航迹相关的点迹;
S5、对每一条航迹,利用与该航迹相关的点迹采用空时联合不敏卡尔曼滤波器进行滤波,获得多个状态估计值;
S6、对于每一条航迹,将其对应的多个状态估计值及状态的一步预测值进行概率互联,获得目标最终状态估计值及航迹存在性概率,更新为当前周期跟踪结果,依据航迹存在性概率判断航迹状态,输出确认航迹并删除终结航迹;
S7、对所有未与任一航迹相关的点迹进行滤波初始化,获得航迹目标状态值并设置初始航迹存在性概率,更新为当前周期跟踪结果;
所述步骤S3中的扇扫雷达空时联合系统模型为线性加减速扇扫模式的雷达跟踪系统模型或三角函数加减速扇扫模式的雷达跟踪系统模型;
所述线性加减速扇扫模式的雷达跟踪系统模型中,天线的扇扫角度范围为[-β,β],匀速阶段天线扫描角速度为α,加速阶段天线加速度为
Figure FDA0003685640300000012
天线加速阶段扫描的角度为
Figure FDA0003685640300000011
天线扫描一周时间为Tradar;设当前周期为第k周期,目标第k-1次扫描到目标的时间为tk-1,第k周期内目标被扫到时的目标状态向量为
Figure FDA0003685640300000027
方位角为A,目标方位角与目标状态向量之间关系为
Figure FDA0003685640300000021
目标状态方程组的表达式为:
Figure FDA0003685640300000022
其中,F(k)是状态转移矩阵;v(k)是过程噪声向量;Γ(k)是噪声分布矩阵;g(tk-1,x(k))表示目标重访间隔与目标状态的非线性关系:
Figure FDA0003685640300000023
其中,Tk为目标重访间隔;f与目标扫描方向相关,若第k周期目标顺时针扫描则f=-1,若目标逆时针扫描则f=1;
所述三角函数加减速扇扫模式的雷达跟踪系统模型中,天线扫描角速度加减速以三角函数形式进行,天线加速阶段或减速阶段扫描的角度范围为α/2,加速阶段天线扫描角速度从0增加到α,加速阶段步长为t=1,速度曲线采用α(t)=α/2(1-cos(πt)),天线扫描一周时间为Tradar;设当前周期为第k周期,目标第k-1次扫描到目标的时间为tk-1,第k周期内目标被扫到时的目标状态向量为
Figure FDA0003685640300000024
方位角为A,目标方位角与目标状态向量之间关系为
Figure FDA0003685640300000025
当| A|<β-Δa时,目标状态方程组为:
Figure FDA0003685640300000026
当| A|>β-Δa时,目标状态方程组为:
Figure FDA0003685640300000031
其中,F(k)是状态转移矩阵;v(k)是过程噪声向量;Γ(k)是噪声分布矩阵;目标重访间隔Tk=tk-tk-1,Δt为天线在第k周期内天线扫过目标到距离目标较近的扫描边界之间的角度所用时间:
Figure FDA0003685640300000032
2.根据权利要求1所述的扇扫雷达空时联合目标跟踪方法,其特征在于:所述步骤S7中初始化后,各个航迹初始均为暂时航迹,所述步骤S6中依据航迹存在性概率判断航迹状态时,若一暂时航迹的航迹存在性概率大于预先设置的确认门限tc,则该航迹被确认为真实航迹,并且保持确认状态直到被终结;若一暂时航迹或真实航迹的航迹存在性概率跌至小于预先设置的航迹终结门限tt,则该航迹被终结。
3.根据权利要求1所述的扇扫雷达空时联合目标跟踪方法,其特征在于:
所述线性加减速扇扫模式的雷达跟踪系统模型或三角函数加减速扇扫模式的雷达跟踪系统模型中,量测方程为:
Figure FDA0003685640300000033
其中,z(k)为目标的量测向量,包括目标相对观测雷达坐标系原点的距离量测rk、方位角量测θk;hk(x(k))表示目标状态与量测之间的数学关系;W(k)为零均值的高斯白噪声序列,表示k时刻量测噪声,包括
Figure FDA0003685640300000034
两个分量,分别为距离、方位角量测噪声;
量测噪声协方差矩阵为:
Figure FDA0003685640300000041
其中,Rk,rr、Rk,θθ分别表示k时刻各量测噪声分量的自协方差,其值分别为
Figure FDA0003685640300000042
Rk,rθ分别表示k时刻各量测噪声分量的互协方差,假设各量测之间是互不相关的,因此各分量互协方差均为0。
4.根据权利要求3所述的扇扫雷达空时联合目标跟踪方法,其特征在于:所述步骤S3中根据其上一周期的航迹存在性概率得到当前周期的航迹存在性概率预测值时,通过马尔科夫链1阶模型表述航迹存在或不存在两种状态的转移概率:
Figure FDA0003685640300000043
其中,P11表示上一周期航迹存在,下一周期航迹依然存在的概率;P12表示上一周期航迹存在,下一周期航迹不存在的概率;P21表示上一周期航迹不存在,下一周期航迹存在的概率;P22表示上一周期航迹不存在,下一周期航迹依然不存在的概率;
已知k-1周期航迹τ存在的概率为
Figure FDA0003685640300000044
Figure FDA0003685640300000045
表示k-1周期航迹τ存在;定义与航迹存在性相关的两种状态:
Figure FDA0003685640300000046
表示目标在第k周期航迹τ存在,
Figure FDA0003685640300000047
表示目标在第k周期航迹τ不存在;
对航迹存在性概率进行预测,得到k周期该航迹存在性概率的预测值为:
Figure FDA0003685640300000048
Figure FDA0003685640300000049
其中,Zτ,k-1表示直到k-1时刻落入航迹τ波门的所有量测的集合,
Figure FDA00036856403000000410
Zτ(k)表示k时刻落入航迹τ相关波门内的
Figure FDA00036856403000000411
个量测的集合,
Figure FDA00036856403000000412
5.根据权利要求4所述的扇扫雷达空时联合目标跟踪方法,其特征在于,所述步骤S3中根据扇扫雷达空时联合系统模型得到状态的一步预测值时,依据UT变换,选择均值周围的采样点δ作为非线性变换的输入,对输出结果求取统计特性,得到目标状态的一步预测:
设第k-1周期中目标τ的状态估计值为
Figure FDA0003685640300000051
协方差为Pτ(k-1|k-1);依据该状态和协方差,产生一组长度为2L+1的采样点δ,各采样点δ均匀分布在第k-1周期状态估计值
Figure FDA0003685640300000052
的附近,其中L是状态向量
Figure FDA0003685640300000053
的维数:
Figure FDA0003685640300000054
Figure FDA0003685640300000055
Figure FDA0003685640300000056
其中,
Figure FDA0003685640300000057
是(L+λ)Pτ(k-1|k-1)的矩阵平方根的第i列;
将上述2L+1个采样点分别作为扇扫雷达空时联合系统模型中目标状态方程组的输入,采用迭代法求解得到2L+1组解,获得2L+1个重访间隔及状态预测值;
目标的重访间隔是2L+1个采样点重访间隔的加权总和:
Figure FDA0003685640300000058
对2L+1个状态预测值进行加权,获得目标最终状态预测值及状态预测协方差:
Figure FDA0003685640300000059
Figure FDA00036856403000000510
其中,
Figure FDA00036856403000000511
Γi,k是依据
Figure FDA00036856403000000512
得到的噪声分布矩阵,状态和协方差的加权权重分别是:
Figure FDA0003685640300000061
α和κ控制采样点δ的传播;β与x的分布有关;
根据扇扫雷达空时联合系统模型中的量测方程,预测测量的采样点是:
Figure FDA0003685640300000062
量测的预测值和相应的协方差分别为:
Figure FDA0003685640300000063
Figure FDA0003685640300000064
其中,
Figure FDA0003685640300000067
状态预测值及量测预测值的互协方差矩阵是:
Figure FDA0003685640300000065
6.根据权利要求5所述的扇扫雷达空时联合目标跟踪方法,其特征在于,所述步骤S3中,将2L+1个采样点分别作为扇扫雷达空时联合系统模型中目标状态方程组的输入,采用迭代法求解得到2L+1组解时,包括如下步骤:
以第i个采样点
Figure FDA0003685640300000066
为输入,若在第k-1周期里目标状态更新时间为tk-1,则天线扫过目标到达扫描边界时间为T1=Tradar·(k-1)-tk-1,迭代法求解状态方程组中目标重访间隔的初始值为:
T′=2·T1=2·(Tradar·(k-1)-tk-1);
以T′为目标重访间隔,依据目标状态方程组中的第二个方程获得迭代初始的状态预测值x′k|k-1为:
Figure FDA0003685640300000071
其中,状态转移矩阵F′(k)依据T′得到:
Figure FDA0003685640300000072
计算第k周期目标被天线扫描到时的目标方位角:
Figure FDA0003685640300000073
依据目标方位角,计算天线在第k周期里从天线扫描边界到扫过目标的时间间隔T2
存在时间差为ΔT=T′-(T1+T2),设置时间差的最大阈值为σ;若|ΔT|<σ,则认为T′与目标真实的重访间隔差距可忽略不计,重访间隔
Figure FDA0003685640300000074
目标状态预测值为
Figure FDA0003685640300000075
否则对T′修正为T′=T′-0.5·ΔT,并返回至以T′为目标重访间隔,依据目标状态方程组中的第二个方程获得迭代初始的状态预测值x′k|k-1步骤,直到|ΔT|<σ;最终获得以第i个采样点为输入获得的目标状态的一步预测值
Figure FDA0003685640300000076
及重访间隔
Figure FDA0003685640300000077
其中,计算天线在第k周期里从天线扫描边界到扫过目标的时间间隔T2时,T2的计算方式与目标方位角和天线扫描方向有关,对于所述线性加减速扇扫模式的雷达跟踪系统模型,T2的表达式为:
Figure FDA0003685640300000078
对于三角函数加减速扇扫模式的雷达跟踪系统模型,若目标位于匀速区即|A|<β-α/2,则T2=1+(β-Δa+f·A)/α;
若目标位于加速或减速区,则通过求解
Figure FDA0003685640300000079
求取T2,若目标位于加速区则T2=Δt,目标位于减速区则T2=Tradar-Δt。
7.根据权利要求6所述的扇扫雷达空时联合目标跟踪方法,其特征在于:所述步骤S5中利用与该航迹相关的点迹采用空时联合不敏卡尔曼滤波器进行滤波时,设一航迹τ有mk个相关的点迹,分别利用mk个点迹应用空时联合不敏卡尔曼滤波器进行滤波,获得该航迹τ相关的mk个状态估计值;
其中,滤波包括获得航迹目标状态的一步预测值和量测预测值后,计算卡尔曼增益:
Figure FDA0003685640300000081
依据一点迹i得到的状态更新值是预测状态与被卡尔曼增益加权的新息之和,得到:
xτ,i(k|k)=xτ(k|k-1)+K(k)[zi(k)-zτ(k|k-1)];
zi(k)表示点迹i;
协方差更新值为:
Figure FDA0003685640300000082
CN201910535351.5A 2019-06-20 2019-06-20 一种扇扫雷达空时联合目标跟踪方法 Active CN112114308B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910535351.5A CN112114308B (zh) 2019-06-20 2019-06-20 一种扇扫雷达空时联合目标跟踪方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910535351.5A CN112114308B (zh) 2019-06-20 2019-06-20 一种扇扫雷达空时联合目标跟踪方法

Publications (2)

Publication Number Publication Date
CN112114308A CN112114308A (zh) 2020-12-22
CN112114308B true CN112114308B (zh) 2022-08-02

Family

ID=73795759

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910535351.5A Active CN112114308B (zh) 2019-06-20 2019-06-20 一种扇扫雷达空时联合目标跟踪方法

Country Status (1)

Country Link
CN (1) CN112114308B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112285696B (zh) * 2020-12-29 2021-05-07 北京海兰信数据科技股份有限公司 一种雷达目标跟踪方法及系统
CN113406617B (zh) * 2021-06-22 2022-09-30 哈尔滨工业大学 一种扇扫雷达多目标连续跟踪方法
CN113406618B (zh) * 2021-06-22 2022-10-25 哈尔滨工业大学 一种tws雷达多目标连续跟踪方法
CN114839626B (zh) * 2021-11-01 2024-06-14 北京遥测技术研究所 一种基于线性多目标方法的航迹关联方法
CN114792379A (zh) * 2022-03-29 2022-07-26 华中光电技术研究所(中国船舶重工集团公司第七一七研究所) 一种基于lmb的红外全景图像多目标航迹处理方法及系统
CN115184896B (zh) * 2022-09-09 2022-11-29 北京海兰信数据科技股份有限公司 一种基于雷达目标幅度加权均值滤波跟踪处理方法及系统

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS58204379A (ja) * 1982-05-25 1983-11-29 Mitsubishi Electric Corp レ−ダ装置
CN103926583A (zh) * 2014-05-04 2014-07-16 中国人民解放军空军预警学院监控系统工程研究所 一种岸基近程雷达自动跟踪处理方法及计算机
CN106980114A (zh) * 2017-03-31 2017-07-25 电子科技大学 无源雷达目标跟踪方法
CN108303692A (zh) * 2018-01-30 2018-07-20 哈尔滨工业大学 一种解多普勒模糊的多目标跟踪方法
CN109557532A (zh) * 2018-10-18 2019-04-02 西安电子科技大学 基于三维霍夫变换的检测前跟踪方法、雷达目标检测系统
CN109633589A (zh) * 2019-01-08 2019-04-16 沈阳理工大学 目标跟踪中基于多模型优化多假设的多目标数据关联方法
CN109655822A (zh) * 2018-11-09 2019-04-19 上海无线电设备研究所 一种改进的航迹起始方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS58204379A (ja) * 1982-05-25 1983-11-29 Mitsubishi Electric Corp レ−ダ装置
CN103926583A (zh) * 2014-05-04 2014-07-16 中国人民解放军空军预警学院监控系统工程研究所 一种岸基近程雷达自动跟踪处理方法及计算机
CN106980114A (zh) * 2017-03-31 2017-07-25 电子科技大学 无源雷达目标跟踪方法
CN108303692A (zh) * 2018-01-30 2018-07-20 哈尔滨工业大学 一种解多普勒模糊的多目标跟踪方法
CN109557532A (zh) * 2018-10-18 2019-04-02 西安电子科技大学 基于三维霍夫变换的检测前跟踪方法、雷达目标检测系统
CN109655822A (zh) * 2018-11-09 2019-04-19 上海无线电设备研究所 一种改进的航迹起始方法
CN109633589A (zh) * 2019-01-08 2019-04-16 沈阳理工大学 目标跟踪中基于多模型优化多假设的多目标数据关联方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
"Multiple hypothesis tracking for multiple target tracking";Blackman 等;《 IEEE Transactions on Aerospace and Electronic Systems》;20041230;正文第5-16页 *
扇扫雷达多目标航迹跟踪算法研究;于文才;《电子世界》;20161230(第08期);正文第183-184页 *

Also Published As

Publication number Publication date
CN112114308A (zh) 2020-12-22

Similar Documents

Publication Publication Date Title
CN112114308B (zh) 一种扇扫雷达空时联合目标跟踪方法
CN108303692B (zh) 一种解多普勒模糊的多目标跟踪方法
CN109901153B (zh) 基于信息熵权和最近邻域数据关联的目标航迹优化方法
CN107045125B (zh) 一种基于预测值量测转换的交互多模型雷达目标跟踪方法
CN109633589A (zh) 目标跟踪中基于多模型优化多假设的多目标数据关联方法
CN104297748B (zh) 一种基于轨迹增强的雷达目标检测前跟踪方法
CN113064155B (zh) 一种空中雷达多目标跟踪下航迹关联的优化方法
CN112613532B (zh) 基于雷达与循环神经网络补全红外融合的动目标跟踪方法
CN107526070A (zh) 天波超视距雷达的多路径融合多目标跟踪算法
CN114662535B (zh) 一种基于变分贝叶斯学习的井下传感器网络目标跟踪方法
CN111274529B (zh) 一种鲁棒的高斯逆威沙特phd多扩展目标跟踪算法
CN110233608B (zh) 一种基于权值自适应的粒子滤波方法和雷达系统
US7928898B2 (en) Method for determining the kinematic state of an object, by evaluating sensor measured values
CN104777465B (zh) 基于b样条函数任意扩展目标形状及状态估计方法
CN108490429B (zh) Tws雷达多目标跟踪方法及系统
CN113280821B (zh) 基于斜率约束和回溯搜索的水下多目标跟踪方法
CN113466817A (zh) 基于b样条曲线拟合的雷达目标航迹自动起始方法
CN111262556B (zh) 一种同时估计未知高斯测量噪声统计量的多目标跟踪方法
US20230118390A1 (en) Method and device for providing tracking data for recognizing the movement of persons and hands for controlling at least one function of a technical system, and sensor system
CN115421153B (zh) 一种基于扩展卡尔曼滤波的激光雷达和uwb组合定位方法及系统
CN113866754B (zh) 基于高斯分布波门的运动目标航迹关联方法
CN108572362B (zh) 一种tws雷达空时联合关联跟踪方法及装置
CN115544425A (zh) 一种基于目标信噪比特征估计的鲁棒多目标跟踪方法
CN110031797B (zh) 用于被动传感系统对具有非连续特性目标的检测跟踪方法
CN109035301B (zh) 一种基于斥力模型修正随机矩阵算法的群组目标跟踪方法

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