CN112114308B - 一种扇扫雷达空时联合目标跟踪方法 - Google Patents
一种扇扫雷达空时联合目标跟踪方法 Download PDFInfo
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 77
- 238000005259 measurement Methods 0.000 claims abstract description 82
- 238000001914 filtration Methods 0.000 claims abstract description 23
- 230000001133 acceleration Effects 0.000 claims description 69
- 238000005070 sampling Methods 0.000 claims description 32
- 239000011159 matrix material Substances 0.000 claims description 30
- 239000013598 vector Substances 0.000 claims description 29
- 230000000875 corresponding effect Effects 0.000 claims description 17
- 230000008569 process Effects 0.000 claims description 15
- 230000007704 transition Effects 0.000 claims description 13
- 230000009466 transformation Effects 0.000 claims description 6
- 230000000737 periodic effect Effects 0.000 claims description 5
- 238000010408 sweeping Methods 0.000 claims description 5
- 238000012790 confirmation Methods 0.000 claims description 4
- 230000001419 dependent effect Effects 0.000 claims description 3
- 230000009191 jumping Effects 0.000 claims description 3
- 230000002596 correlated effect Effects 0.000 claims description 2
- 230000006870 function Effects 0.000 description 23
- 230000007246 mechanism Effects 0.000 description 10
- 238000004088 simulation Methods 0.000 description 6
- 238000006243 chemical reaction Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000000342 Monte Carlo simulation Methods 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 230000028161 membrane depolarization Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Images
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
- G01S13/72—Radar-tracking systems; Analogous systems for two-dimensional tracking, e.g. combination of angle and range tracking, track-while-scan radar
- G01S13/723—Radar-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/726—Multiple 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中的扇扫雷达空时联合系统模型为线性加减速扇扫模式的雷达跟踪系统模型或三角函数加减速扇扫模式的雷达跟踪系统模型。
优选地,所述线性加减速扇扫模式的雷达跟踪系统模型中,天线的扇扫角度范围为[-β,β],匀速阶段天线扫描角速度为α,加速阶段天线加速度为天线加速阶段扫描的角度为天线扫描一周时间为Tradar;设当前周期为第k周期,目标第k-1次扫描到目标的时间为tk-1,第k周期内目标被扫到时的目标状态向量为方位角为A,目标方位角与目标状态向量之间关系为
目标状态方程组的表达式为:
其中,F(k)是状态转移矩阵;v(k)是过程噪声向量;Γ(k)是噪声分布矩阵;g(tk-1,x(k))表示目标重访间隔与目标状态的非线性关系:
其中,Tk为目标重访间隔;f与目标扫描方向相关,若第k周期目标顺时针扫描则f=-1,若目标逆时针扫描则f=1。
优选地,所述三角函数加减速扇扫模式的雷达跟踪系统模型中,天线扫描角速度加减速以三角函数形式进行,天线加速阶段或减速阶段扫描的角度范围为α/2,加速阶段天线扫描角速度从0增加到α,加速阶段步长为t=1,速度曲线采用α(t)=α/2(1-cos(πt)),天线扫描一周时间为Tradar;设当前周期为第k周期,目标第k-1次扫描到目标的时间为tk-1,第k周期内目标被扫到时的目标状态向量为方位角为A,目标方位角与目标状态向量之间关系为
当A|<β-Δa时,目标状态方程组为:
当A|>β-Δa时,目标状态方程组为:
其中,F(k)是状态转移矩阵;v(k)是过程噪声向量;Γ(k)是噪声分布矩阵;目标重访间隔Tk=tk-tk-1,Δt为天线在第k周期内天线扫过目标到距离目标较近的扫描边界之间的角度所用时间:
优选地,所述线性加减速扇扫模式的雷达跟踪系统模型或三角函数加减速扇扫模式的雷达跟踪系统模型中,量测方程为:
其中,z(k)为目标的量测向量,包括目标相对所述观测雷达坐标系原点的距离量测rk、方位角量测θk;hk(x(k))表示目标状态与量测之间的数学关系;W(k)为零均值的高斯白噪声序列,表示k时刻量测噪声,包括两个分量,分别为距离、方位角量测噪声;
量测噪声协方差矩阵为:
优选地,所述步骤S3中根据其上一周期的航迹存在性概率得到当前周期的航迹存在性概率预测值时,通过马尔科夫链1阶模型表述航迹存在或不存在两种状态的转移概率:
其中,P11表示上一周期航迹存在,下一周期航迹依然存在的概率; P12表示上一周期航迹存在,下一周期航迹不存在的概率;P21表示上一周期航迹不存在,下一周期航迹存在的概率;P22表示上一周期航迹不存在,下一周期航迹依然不存在的概率;
对航迹存在性概率进行预测,得到k周期该航迹存在性概率的预测值为:
优选地,所述步骤S3中根据扇扫雷达空时联合系统模型得到状态的一步预测值时,选择均值周围的采样点δ作为非线性变换的输入,对输出结果求取统计特性,得到目标状态的一步预测:
设第k-1周期中目标τ的状态估计值为协方差为 Pτ(k-1|k-1);依据该状态和协方差,产生一组长度为2L+1的采样点δ,各采样点δ均匀分布在第k-1周期状态估计值的附近,其中L 是状态向量的维数:
将上述2L+1个采样点分别作为扇扫雷达空时联合系统模型中目标状态方程组的输入,采用迭代法求解得到2L+1组解,获得2L+1个重访间隔及状态预测值;
目标的重访间隔是2L+1个采样点重访间隔的加权总和:
对2L+1个状态预测值进行加权,获得目标最终状态预测值及状态预测协方差:
α和κ控制采样点δ的传播;β与x的分布有关;
根据扇扫雷达空时联合系统模型中的量测方程,预测测量的采样点是:
量测的预测值和相应的协方差分别为:
其中,
状态预测值及量测预测值的互协方差矩阵是:
优选地,所述步骤S3中,将2L+1个采样点分别作为扇扫雷达空时联合系统模型中目标状态方程组的输入,采用迭代法求解得到2L+1 组解时,包括如下步骤:
T′=2·T1=2·(Tradar·(k-1)-tk-1);
以T′为目标重访间隔,依据目标状态方程组中的第二个方程获得迭代初始的状态预测值x′k|k-1为:
其中,状态转移矩阵F′(k)依据T′得到:
计算第k周期目标被天线扫描到时的目标方位角:
依据目标方位角,计算天线在第k周期里从天线扫描边界到扫过目标的时间间隔T2;
存在时间差为ΔT=T′-(T1+T2),设置时间差的最大阈值为σ;若 |ΔT|<σ,则认为T′与目标真实的重访间隔差距可忽略不计,重访间隔目标状态预测值为否则对T′修正为T′=T′-0.5·ΔT,并返回至以T′为目标重访间隔,依据目标状态方程组中的第二个方程获得迭代初始的状态预测值x′k|k-1步骤,直到|ΔT|<σ;最终获得以第i个采样点为输入获得的目标状态的一步预测值及重访间隔
其中,计算天线在第k周期里从天线扫描边界到扫过目标的时间间隔T2时,T2的计算方式与目标方位角和天线扫描方向有关,对于所述线性加减速扇扫模式的雷达跟踪系统模型,T2的表达式为:
对于三角函数加减速扇扫模式的雷达跟踪系统模型,若目标位于匀速区即|A|<β-α/2,则T2=1+(β-Δa+f·A)/α;
优选地,所述步骤S5中利用与该航迹相关的点迹采用空时联合不敏卡尔曼滤波器进行滤波时,设一航迹τ有mk个相关的点迹,分别利用mk个点迹应用空时联合不敏卡尔曼滤波器进行滤波,获得该航迹τ相关的mk个状态估计值;
其中,滤波包括获得航迹目标状态的一步预测值和量测预测值后,计算卡尔曼增益:
依据一点迹i得到的状态更新值是预测状态与被卡尔曼增益加权的新息之和,得到:
xτ,i(k|k)=xτ(k|k-1)+K(k)[zi(k)-zτ(k|k-1)];
协方差更新值为:
(三)有益效果
本发明的上述技术方案具有如下优点:本发明提供了一种扇扫雷达空时联合目标跟踪方法,该方法采用递推循环的方式实现多目标跟踪,采用增加了描述时间的方程对状态方程进行约束的扇扫雷达空时联合系统模型对扇扫雷达跟踪系统进行描述,并提出了相应的滤波方法,有效解决了传统目标跟踪方法中系统模型与扇扫雷达存在模型不匹配的问题,提高了目标跟踪在准确性以及真实航迹确认速度。
附图说明
图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阶模型表述航迹存在或航迹不存在这两种状态的转移概率:
其中,P11表示上一周期航迹存在,下一周期航迹依然存在的概率; P12表示上一周期航迹存在,下一周期航迹不存在的概率;P21表示上一周期航迹不存在,下一周期航迹存在的概率;P22表示上一周期航迹不存在,下一周期航迹依然不存在的概率;P11、P12、P21和P22的具体数值可由工程师根据实际情况进行设置。
对航迹存在性概率进行预测,得到k周期该航迹存在性概率的预测值为:
考虑到扇扫雷达的通用性,本发明针对线性加减速扇扫雷达及三角函数加减速扇扫雷达进行了建模。优选地,所述步骤S3中的扇扫雷达空时联合系统模型为线性加减速扇扫模式的雷达跟踪系统模型或三角函数加减速扇扫模式的雷达跟踪系统模型。
优选地,所述线性加减速扇扫模式的雷达系统模型中,天线的扇扫角度范围为[-β,β],匀速阶段天线扫描角速度为α,加速阶段天线加速度为天线加速阶段扫描的角度为天线扫描一周时间为Tradar;设当前周期为第k周期,目标第k-1次扫描到目标的时间为tk-1,第k周期内目标被扫到时的目标状态向量为方位角为A,目标方位角与目标状态向量之间关系为
目标状态方程组的表达式为:
其中,F(k)是状态转移矩阵;v(k)是过程噪声向量;Γ(k)是噪声分布矩阵;g(tk-1,x(k))表示目标重访间隔与目标状态的非线性关系:
优选地,所述三角函数加减速扇扫模式的雷达跟踪系统模型中,天线扫描角速度加减速以三角函数形式进行,天线加速阶段或减速阶段扫描的角度范围为α/2,加速阶段天线扫描角速度从0增加到α,加速阶段步长为t=1,速度曲线采用α(t)=α/2(1-cos(πt)),天线扫描一周时间为Tradar;设当前周期为第k周期,目标第k-1次扫描到目标的时间为tk-1,第k周期内目标被扫到时的目标状态向量为方位角为A,目标方位角与目标状态向量之间关系为
当A|<β-Δa时,目标状态方程组为:
当A|>β-Δa时,目标状态方程组为:
其中,F(k)是状态转移矩阵;v(k)是过程噪声向量;Γ(k)是噪声分布矩阵;目标重访间隔Tk=tk-tk-1,Δt为天线在第k周期内天线扫过目标到距离目标较近的扫描边界之间的角度所用时间:
优选地,线性加减速扇扫模式的雷达跟踪系统模型或三角函数加减速扇扫模式的雷达跟踪系统模型中,两种扇扫方式的量测方程是一致的,量测方程的表达式为:
其中,z(k)为目标的量测向量,包括目标相对所述观测雷达坐标系原点的距离量测rk、方位角量测θk;hk(x(k))表示目标状态与量测之间的数学关系;W(k)为零均值的高斯白噪声序列,表示k时刻量测噪声,包括两个分量,分别为距离、方位角量测噪声;
量测噪声协方差矩阵为:
优选地,所述步骤S3中根据扇扫雷达空时联合系统模型得到状态的一步预测值时,依据UT变换,选择均值周围的采样点δ作为非线性变换的输入,对输出结果求取统计特性,得到目标状态的一步预测。
设第k-1周期中目标τ的状态估计值为协方差为 Pτ(k-1|k-1);依据该状态和协方差,产生一组长度为2L+1的采样点δ,各采样点δ均匀分布在第k-1周期,即上一周期,更新后的状态估计值的附近,其中L是状态向量的维数;
将上述2L+1个采样点分别作为扇扫雷达空时联合系统模型中目标状态方程组的输入,采用迭代法求解得到2L+1组解,获得2L+1个重访间隔及状态预测值;
目标的重访间隔是2L+1个采样点重访间隔的加权总和:
对2L+1个状态预测值进行加权,获得目标最终状态预测值及状态预测协方差:
其中,
α和κ控制采样点δ的传播;β与x的分布有关;
根据扇扫雷达空时联合系统模型中的量测方程,预测测量的采样点是:
量测的预测值和相应的协方差分别为:
其中,
状态预测值及量测预测值的互协方差矩阵是:
优选地,如图2所示,所述步骤S3中,将2L+1个采样点分别作为扇扫雷达空时联合系统模型中目标状态方程组的输入,采用迭代法求解得到2L+1组解时,对于线性加减速扇扫模式的雷达跟踪系统模型,包括如下步骤:
T′=2·T1=2·(Tradar·(k-1)-tk-1);
2)以T′为目标重访间隔,依据目标状态方程组中的第二个方程获得迭代初始的状态预测值x′k|k-1为:
其中,状态转移矩阵F′(k)依据T′得到:
3)计算第k周期目标被天线扫描到时的目标方位角:
其中,x′k|k-1(1)、x′k|k-1(3)分别表示4x1维向量x′k|k-1的第1个值和第3 个值;
4)依据目标方位角,计算天线在第k周期里从天线扫描边界到扫过目标的时间间隔T2,T2的计算方式与目标方位角和天线扫描方向有关,表达式为:
否则对T′修正为T′=T′-0.5·ΔT,并返回至以T′为目标重访间隔,依据目标状态方程组中的第二个方程获得迭代初始的状态预测值x′k|k-1步骤,即在修正T′后重新进行步骤2)到5),直到|ΔT|<σ;最终获得以第 i个采样点为输入获得的目标状态的一步预测值及重访间隔
优选地,所述步骤S3中,将2L+1个采样点分别作为扇扫雷达空时联合系统模型中目标状态方程组的输入,采用迭代法求解得到2L+1 组解时,对于三角函数加减速扇扫模式的雷达跟踪系统模型,包括如下步骤:
T′=2·T1=2·(Tradar·(k-1)-tk-1);
2)以T′为目标重访间隔,依据状态方程组中的第二个方程获得迭代初始的状态预测值x′k|k-1,即为:
其中,状态转移矩阵F′(k)依据T′得到,即:
3)计算第k周期目标被天线扫描到时目标方位角:
4)判断目标所在的区域,依据目标方位角,计算在天线k周期内,从天线扫描边缘扫到目标的时间T2:
若目标位于匀速区即|A|<β-α/2,则T2=1+(β-Δa+f·A)/α;
2、比较角度误差的阈值是σ1。它可能是一个非常小的值,可按需自行设定。如果角度误差小于阈值即|ε|<σ1,那么认为角度误差可以忽略,即Δt=Δt′。否则需要对Δt′进行调整,令重新执行第1 步骤,计算角度误差并对时间Δt′进行调整,直至角度误差满足|ε|<σ1,若目标位于加速区那么T2=Δt,目标位于减速区T2=Tradar-Δt。
5)依据目标方位角A获得时间T2,求取时间差为ΔT=T′-(T1+T2),设置时间差的最大阈值为σ2。若|ΔT|<σ2则认为T′与目标真实的重访间隔差距非常小可以忽略不计,重访间隔目标状态预测值为否则对T′修正为T′=T′-0.5·ΔT,并返回至以T′为目标重访间隔,依据状态方程组中的第二个方程获得迭代初始的状态预测值 x′k|k-1步骤,即在修正T′后重新进行步骤2)到5),直到|ΔT|<σ2。至此我们得到第i个采样点为输入获得的目标状态的一步预测值及重访间隔
优选地,步骤S4中对每一条航迹,依据所述量测预测值建立相关波门,从步骤S1中得到的量测信息中获取与航迹相关的观测点时,对 k时刻雷达接收所有的点迹进行波门相关,相关波门依据被跟踪目标在当前时刻所在位置的预测值为中心,确立一个量测向量可能出现的范围的决策门限,所有落入该门限的量测值均考虑为候选的量测,而在门限值之外的量测则认为是杂波。相关波门通常用来滤除杂波,找到所有可能来自于目标的量测值。常用的相关波门有矩形波门、环形波门、椭圆(球)波门、极坐标系下的扇形波门等,这里我们将使用椭圆波门。
一条航迹在k时刻落入波门内的点迹可能有多个,也可能没有。
优选地,若k时刻有mk个点迹落入航迹τ的相关波门,对于其中的第i个点迹zi(k)(这里i与前文中的i意义不同,表征落入相关波门内的某一点迹,前文中指代UT变换过程中某一采样点),利用点迹zi(k) 对航迹τ进行状态更新,得到该点迹对航迹τ的状态更新值同时得到该点迹来自航迹τ所跟踪目标的似然函数
进一步优选地,所述步骤S5中对每一条航迹,利用与该航迹相关的点迹进行滤波过程,获得多个状态估计值时,设一航迹τ有个相关的点迹,分别利用个点迹应用空时联合不敏卡尔曼滤波器 (STJ-UKF)进行滤波,获得该航迹τ相关的个状态估计值。STJ-UKF 方法前半部分内容即状态一步预测过程,这里不再重复描述。
滤波还包括获得航迹目标状态的一步预测值和量测预测值后,计算卡尔曼增益:
依据某一点迹i,即zi(k),得到的状态更新值是预测状态与被卡尔曼增益加权的新息之和,得到:
xτ,i(k|k)=xτ(k|k-1)+K(k)[zi(k)-zτ(k|k-1)];
相应的协方差更新值为:
优选地,步骤S6中对于每一条航迹,将其对应的多个状态估计值及状态的一步预测值进行概率互联时,如果点迹i落在目标τ的波门内,则表示点迹i与目标τ相关,点迹i可能是航迹τ的观测值,定义事件为波门内的第i个点迹来自目标τ,其他点迹对该航迹为杂波。
点迹i来自航迹τ的先验概率为:
在多目标跟踪系统中,目标τ波门中的点迹i可能来自另一个目标。在线性化方法下,通过排除目标τ的所有可能来源来修改量测i所在位置的杂波密度。跟踪τ的验证门中的测量i的先验散射体测量密度,由Ωτ(i)表示,可以表示为:
结合综合概率数据互联算法(IPDA),目标存在的后验概率被修改为:
其中,
其中,
至此获得落入目标τ的波门内的各个量测属于该目标的后验概率,以及航迹τ存在的后验概率值。
对所有航迹均进行以上航迹更新过程,依据航迹存在性概率对航迹进行管理。某些点迹可能未落入任何一条航迹的相关波门内,则需要进行航迹初始化。
考虑到未被用来更新航迹的点迹可能代表着出现新航迹,对其执行航迹初始化,可获得新航迹目标初始状态。优选地,步骤S7中对所有未与任一航迹相关的点迹,依据观测雷达当前周期和上一周期所获得的未与航迹相关的点迹进行滤波初始化时,利用未与航迹相关的点进行航迹初始化获得新航迹,以第一周期和第二周期的点迹为例,第一次扫描雷达得到n1个点迹,第二次扫描雷达得到n2个点迹,利用第一次扫描时刻的第i(1≤i≤n1)个点迹和第二次扫描时刻的第j(1≤j≤n2)个点迹获得第τ(1≤τ≤n1×n2)条航迹在笛卡尔坐标系下第二次扫描时刻的状态向量。
对应的状态协方差为:
其中下标1、2分别表示第一周期和第二周期对应的值。T表示雷达扫描周期。和是笛卡尔坐标下的量测沿x,y方向的位置量测信息,是通过无偏量测转换方法将雷达所获得量测(径向距离、方位角) 转换到笛卡尔坐标系下得到的转换位置量测,转换公式为:
对应的协方差矩阵为:
其中各分量的自协方差及互协方差:
上标“c”代表与转换量测相关的向量、矩阵和函数。初始化的状态误差协方差为:
T可以通过两个量测的时间差得到。
航迹存在性概率初始值可由实际需要进行设定。
综上,本申请提供了一种扇扫雷达空时联合目标跟踪方法,该方法包括如下优势:
(1)针对扇扫雷达重访间隔不确定导致系统模型不匹配的问题,通过增加时间约束方程,分别针对线性加减速扇扫雷达与三角函数扇扫机制雷达建立了适用于扇扫雷达的状态方程组,准确描述两种扫描机制下目标状态随时间变化的关系。
(2)针对新建立的系统模型,研究了对应的空时联合状态估计滤波方法STJ-UKF,通过迭代法求取非线性状态方程组的解,同时得到目标状态预测值及重访间隔,该过程可以看作一个广义的非线性函数,因此可应用不敏变换计算其输出变量的统计特性,最后获得目标状态估计结果。
(3)将滤波方法应用于线性多目标联合概率密度互联(LMIPDA) 框架,研究了适应于的扇扫雷达的多目标跟踪方法STJ-LMIPDA方法
为验证本发明的效果,利用仿真数据进行蒙特卡洛实验。仿真实验监视区域内设置三个目标,仿真试验中的目标沿直线匀速运动,运动轨迹如图3所示。三个目标初始状态如表1所示。
表1 三 个目标的初始状态
进行100次蒙脱卡洛仿真,每次模拟运行有100次扫描。监视区域中的杂波密度为6.0×10-6/scan/m2且杂波满足均匀的泊松分布。目标初始位置和速度总结在表1 中给出。过程噪声设置为ds=0.01m/s2,范围和方位角测量标准偏差分别为σr=50m,σθ=1deg。天线匀速扫描阶段扫描角速度为α=π/6 rad/s,天线扫描范围为[-π/3,π/3],迭代过程时间差的阈值为σ2=0.001s。
将区分奇偶扫描周期,以两倍天线扫描周期作为重访间隔的传统 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中的扇扫雷达空时联合系统模型为线性加减速扇扫模式的雷达跟踪系统模型或三角函数加减速扇扫模式的雷达跟踪系统模型;
所述线性加减速扇扫模式的雷达跟踪系统模型中,天线的扇扫角度范围为[-β,β],匀速阶段天线扫描角速度为α,加速阶段天线加速度为天线加速阶段扫描的角度为天线扫描一周时间为Tradar;设当前周期为第k周期,目标第k-1次扫描到目标的时间为tk-1,第k周期内目标被扫到时的目标状态向量为方位角为A,目标方位角与目标状态向量之间关系为
目标状态方程组的表达式为:
其中,F(k)是状态转移矩阵;v(k)是过程噪声向量;Γ(k)是噪声分布矩阵;g(tk-1,x(k))表示目标重访间隔与目标状态的非线性关系:
其中,Tk为目标重访间隔;f与目标扫描方向相关,若第k周期目标顺时针扫描则f=-1,若目标逆时针扫描则f=1;
所述三角函数加减速扇扫模式的雷达跟踪系统模型中,天线扫描角速度加减速以三角函数形式进行,天线加速阶段或减速阶段扫描的角度范围为α/2,加速阶段天线扫描角速度从0增加到α,加速阶段步长为t=1,速度曲线采用α(t)=α/2(1-cos(πt)),天线扫描一周时间为Tradar;设当前周期为第k周期,目标第k-1次扫描到目标的时间为tk-1,第k周期内目标被扫到时的目标状态向量为方位角为A,目标方位角与目标状态向量之间关系为
当| A|<β-Δa时,目标状态方程组为:
当| A|>β-Δa时,目标状态方程组为:
其中,F(k)是状态转移矩阵;v(k)是过程噪声向量;Γ(k)是噪声分布矩阵;目标重访间隔Tk=tk-tk-1,Δt为天线在第k周期内天线扫过目标到距离目标较近的扫描边界之间的角度所用时间:
2.根据权利要求1所述的扇扫雷达空时联合目标跟踪方法,其特征在于:所述步骤S7中初始化后,各个航迹初始均为暂时航迹,所述步骤S6中依据航迹存在性概率判断航迹状态时,若一暂时航迹的航迹存在性概率大于预先设置的确认门限tc,则该航迹被确认为真实航迹,并且保持确认状态直到被终结;若一暂时航迹或真实航迹的航迹存在性概率跌至小于预先设置的航迹终结门限tt,则该航迹被终结。
3.根据权利要求1所述的扇扫雷达空时联合目标跟踪方法,其特征在于:
所述线性加减速扇扫模式的雷达跟踪系统模型或三角函数加减速扇扫模式的雷达跟踪系统模型中,量测方程为:
其中,z(k)为目标的量测向量,包括目标相对观测雷达坐标系原点的距离量测rk、方位角量测θk;hk(x(k))表示目标状态与量测之间的数学关系;W(k)为零均值的高斯白噪声序列,表示k时刻量测噪声,包括两个分量,分别为距离、方位角量测噪声;
量测噪声协方差矩阵为:
4.根据权利要求3所述的扇扫雷达空时联合目标跟踪方法,其特征在于:所述步骤S3中根据其上一周期的航迹存在性概率得到当前周期的航迹存在性概率预测值时,通过马尔科夫链1阶模型表述航迹存在或不存在两种状态的转移概率:
其中,P11表示上一周期航迹存在,下一周期航迹依然存在的概率;P12表示上一周期航迹存在,下一周期航迹不存在的概率;P21表示上一周期航迹不存在,下一周期航迹存在的概率;P22表示上一周期航迹不存在,下一周期航迹依然不存在的概率;
对航迹存在性概率进行预测,得到k周期该航迹存在性概率的预测值为:
5.根据权利要求4所述的扇扫雷达空时联合目标跟踪方法,其特征在于,所述步骤S3中根据扇扫雷达空时联合系统模型得到状态的一步预测值时,依据UT变换,选择均值周围的采样点δ作为非线性变换的输入,对输出结果求取统计特性,得到目标状态的一步预测:
将上述2L+1个采样点分别作为扇扫雷达空时联合系统模型中目标状态方程组的输入,采用迭代法求解得到2L+1组解,获得2L+1个重访间隔及状态预测值;
目标的重访间隔是2L+1个采样点重访间隔的加权总和:
对2L+1个状态预测值进行加权,获得目标最终状态预测值及状态预测协方差:
α和κ控制采样点δ的传播;β与x的分布有关;
根据扇扫雷达空时联合系统模型中的量测方程,预测测量的采样点是:
量测的预测值和相应的协方差分别为:
其中,
状态预测值及量测预测值的互协方差矩阵是:
6.根据权利要求5所述的扇扫雷达空时联合目标跟踪方法,其特征在于,所述步骤S3中,将2L+1个采样点分别作为扇扫雷达空时联合系统模型中目标状态方程组的输入,采用迭代法求解得到2L+1组解时,包括如下步骤:
T′=2·T1=2·(Tradar·(k-1)-tk-1);
以T′为目标重访间隔,依据目标状态方程组中的第二个方程获得迭代初始的状态预测值x′k|k-1为:
其中,状态转移矩阵F′(k)依据T′得到:
计算第k周期目标被天线扫描到时的目标方位角:
依据目标方位角,计算天线在第k周期里从天线扫描边界到扫过目标的时间间隔T2;
存在时间差为ΔT=T′-(T1+T2),设置时间差的最大阈值为σ;若|ΔT|<σ,则认为T′与目标真实的重访间隔差距可忽略不计,重访间隔目标状态预测值为否则对T′修正为T′=T′-0.5·ΔT,并返回至以T′为目标重访间隔,依据目标状态方程组中的第二个方程获得迭代初始的状态预测值x′k|k-1步骤,直到|ΔT|<σ;最终获得以第i个采样点为输入获得的目标状态的一步预测值及重访间隔
其中,计算天线在第k周期里从天线扫描边界到扫过目标的时间间隔T2时,T2的计算方式与目标方位角和天线扫描方向有关,对于所述线性加减速扇扫模式的雷达跟踪系统模型,T2的表达式为:
对于三角函数加减速扇扫模式的雷达跟踪系统模型,若目标位于匀速区即|A|<β-α/2,则T2=1+(β-Δa+f·A)/α;
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)
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)
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 | 上海无线电设备研究所 | 一种改进的航迹起始方法 |
-
2019
- 2019-06-20 CN CN201910535351.5A patent/CN112114308B/zh active Active
Patent Citations (7)
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)
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 |