CN111708015B - 一种多径效应下的低空目标跟踪滤波方法 - Google Patents

一种多径效应下的低空目标跟踪滤波方法 Download PDF

Info

Publication number
CN111708015B
CN111708015B CN202010671798.8A CN202010671798A CN111708015B CN 111708015 B CN111708015 B CN 111708015B CN 202010671798 A CN202010671798 A CN 202010671798A CN 111708015 B CN111708015 B CN 111708015B
Authority
CN
China
Prior art keywords
target
state
single pulse
ratio
steps
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
CN202010671798.8A
Other languages
English (en)
Other versions
CN111708015A (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 CN202010671798.8A priority Critical patent/CN111708015B/zh
Publication of CN111708015A publication Critical patent/CN111708015A/zh
Application granted granted Critical
Publication of CN111708015B publication Critical patent/CN111708015B/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
    • 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
    • G01S7/00Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/41Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
    • G01S7/418Theoretical aspects

Abstract

本发明涉及雷达目标跟踪滤波技术领域。本发明的一种多径效应下的低空目标跟踪滤波方法包括:建立基于单脉冲比观测的系统模型,所述系统模型中的观测方程为单脉冲比关于目标状态的函数;应用基于分段二分法的初始化方法进行滤波器初始化,得到目标可能的初始状态;对得到的每个可能的初始状态分别应用不敏卡尔曼滤波器进行滤波;将所述不敏卡尔曼滤波器的滤波结果应用静态多模型估计器得到目标高度估计,实现低空目标跟踪。本发明能够减小目标高度估计的误差,对目标进行有效地跟踪,并且上述初始化方法使运行时间得以缩短,同时也减少了可能遗漏的目标状态。

Description

一种多径效应下的低空目标跟踪滤波方法
技术领域
本发明涉及雷达目标跟踪滤波技术领域,尤其涉及一种多径效应下的低空目标跟踪滤波。
背景技术
雷达领域经典的测角方法是单脉冲和差波束比幅法,但是应用这个方法跟踪海面低空飞行目标时,来自海面的镜面反射和漫反射使得雷达跟踪目标性能大幅下降,这种多径效应特别是给目标俯仰角的测量带来很大偏差。而原因就在于直接回波和反射回波的干涉作用,当二者相位反相互抵消时,总回波包含的有用信息很少,这时单脉冲和差波束比幅系统的测角结果出现大的尖峰误差,可能因此丢失待跟踪的目标。
在多径效应下海面低空目标跟踪问题中,由于来自海面反射的回波也进入了天线主波束的半功率波瓣宽度内,因而使得常规单脉冲和差比幅雷达测角性能严重下降,目标高度估计中出现近似周期性的大的尖峰误差和近似固定的偏差;而应用窄波束技术时仍有一些场景存在反射回波无法从总回波中剔除的问题,比如雷达高度100米,目标距离海面60米,距离雷达水平距离50公里时,反射直射路径差约为0.24米,而C波段6GHz雷达带宽600MHz对应的距离分辨力约0.25米,此情况或者目标更低更远的情况下C波段6GHz雷达便很难从时间延迟上区分直射路径和反射路径,导致雷达系统对目标角度及高度测量性能严重下降。目前的信号处理技术是以雷达发射多个频率来解决多径场景中出现的高度模糊问题,而单频率单脉冲雷达则很难对高度解模糊。
因此,针对以上不足,需要提供一种能解决低站址单频率单脉冲雷达跟踪超低空远距离目标时性能下降的方法,实现对目标有效的跟踪。
发明内容
本发明要解决的技术问题在于应用低站址单频率单脉冲雷达跟踪超低空远距离目标时俯仰角及高度量测效果由于多径效应而严重下降,针对现有技术中的缺陷,提供一种多径效应下的低空目标跟踪滤波方法。
本发明提供了一种多径效应下的低空目标跟踪滤波方法,所述方法包括:
建立基于单脉冲比观测的系统模型,所述系统模型中的观测方程为单脉冲比关于目标状态的函数;
应用基于分段二分法的初始化方法进行滤波器初始化,得到目标可能的初始状态;
对得到的每个可能的初始状态分别应用不敏卡尔曼滤波器进行滤波;
将所述不敏卡尔曼滤波器的滤波结果应用静态多模型估计器得到目标高度估计,实现低空目标跟踪。
可选地,所述的建立基于单脉冲比观测的系统模型包括:
建立低空目标跟踪的状态方程:
Xk=Fk-1Xk-1+Vk-1
其中,Xk为k时刻目标状态向量,Fk-1为k-1时刻状态转移矩阵;Vk-1为k-1时刻过程噪声,所述噪声为零均值高斯白噪声,协方差矩阵为Qk-1
对于恒速度模型目标,状态转移矩阵为:
Figure BDA0002581600740000021
其中,T为雷达采样时间间隔;
将所述目标状态向量作为所述观测方程的输入,将单脉冲和差波束比幅雷达获得的单脉冲比量作为所述观测方程的输出,建立k时刻观测方程为:
Zk=hk(Xk)+Wk
其中,Zk为k时刻观测向量:
Zk=[ρk rk]T
其中,ρk是k时刻径向距离量测值,rk是k时刻实单脉冲比量测值;
Wk为k时刻观测噪声,所述观测噪声为零均值高斯白噪声,协方差矩阵为Rk
hk()满足下式:
hk(Xk)=[hρ(Xk) hratio(Xk)]k T
其中,
Figure BDA0002581600740000031
其中,hr为雷达高度,xk为k时刻目标水平距离,yk为k时刻目标垂直高度;
hratio()是表征单脉冲比量测量与状态向量之间关系的函数。
可选地,所述的hratio()的获得方法包括:通过以下表达式由目标状态向量得到单脉冲比以表征此函数关系:
Figure BDA0002581600740000032
其中,F和FΔ分别为雷达和波束与差波束的电压增益方向性函数,即天线方向图函数;θ为目标偏向于天线视轴的角度;
Figure BDA0002581600740000033
为镜像目标偏向于天线视轴的角度;θ0为天线视轴的偏置角度;θt为目标相对于海面水平线的真实俯仰角值;θm为镜像目标相对于海面水平线的俯仰角值;v为海面总反射系数;Δρ为反射路径与直射路径之间的路径长度差;φl是由反射路径与直射路径之间的路径长度差造成的相位差;λ为电磁波波长。
可选地,所述的应用基于分段二分法的初始化方法进行滤波器初始化,得到目标可能的初始状态包括:
确定目标俯仰角的搜索范围;
按照每个单脉冲比波峰波谷点对应的高度值将所述目标俯仰角的搜索范围分段,得到若干段单脉冲比的单调曲线;
在单脉冲比的每一段单调曲线上用二分法搜索目标可能的初始状态,即:对于所述的每一段单调曲线,计算目标状态对应的单脉冲比量与单脉冲比量测量的差值,将第一个差值符合单脉冲比误差范围的目标状态作为该段单调曲线上目标可能的初始状态。
可选地,所述的确定目标俯仰角的搜索范围包括:
选取天线主波束半功率波束宽度内雷达能照射到目标的区域作为所述目标俯仰角的搜索范围;
若天线主波束半功率波束宽度下边界照射到海面以下,则将所述下边界调整为海面水平线处的俯仰角。
可选地,所述的单脉冲比波峰波谷点对应的高度值的确定方法包括:
当反射路径与直射路径之间的路径长度差是半波长的偶数倍时,对应为单脉冲比波峰点;当路径长度差是半波长的奇数倍时,对应为单脉冲比波谷点;由反射路径长度ρm即可得到目标此刻高度值h。
可选地,所述的在单脉冲比的每一段单调曲线上用二分法进行搜索目标可能的初始状态的具体方法为:
步骤S231、计算目标状态对应的单脉冲比量
Figure BDA0002581600740000041
与单脉冲比量测量r的差值,将其记作f(h):
Figure BDA0002581600740000042
对于任意一个单脉冲比单调段,其两个端点对应的目标高度值分别记作hmax和hmin,且hmax>hmin
若f(hmin)*f(hmax)<0且min(|f(hmin)|,|f(hmax)|)>e,e为允许误差,则执行步骤S232;否则,执行步骤S235;
步骤S232、计算hmax、hmin的中点值hhalf
hhalf=(hmin+hmax)/2
若|f(hhalf)|<e,则执行步骤S233;否则执行步骤S234;
步骤S233、输出hhalf对应的目标状态作为一个可能的目标初始状态,完成此单调段的搜索;
步骤S234、判断f(hhalf)*f(hmin)<0是否成立:
若成立,则令hmax=hhalf,然后返回至步骤S232,直至找到符合误差条件的目标状态,完成此单调段的搜索;
否则,令hmin=hhalf,然后返回至步骤S232,直至找到符合误差条件的目标状态,完成此单调段的搜索;返回至计算hhalf的步骤,然后重新执行后续步骤直至找到符合误差条件的目标状态,完成此单调段的搜索;
步骤S235、判断|f(hmin)|<|f(hmax)|是否成立;
若成立,则输出hmin对应的目标状态作为一个可能的目标初始状态,完成此单调段的搜索;
否则,输出hmax对应的目标状态作为一个可能的目标初始状态,完成此单调段的搜索。
可选地,所述的静态多模型估计器满足如下公式:
Figure BDA0002581600740000051
其中,P为概率密度;k代表采样时间;ξk为k时刻工作模式,ξk∈M={m1,m2,…,mz};M为工作模式集合,mi为第i个工作模式;z为工作模式的个数。
可选地,所述的对得到的每个可能的初始状态分别应用不敏卡尔曼滤波器进行滤波的具体方法为:
将每个可能的初始状态作为一个模式;
对各模式分别进行不敏卡尔曼滤波
Figure BDA0002581600740000061
模式j的状态一步预测
Figure BDA0002581600740000062
和一步预测协方差矩阵
Figure BDA0002581600740000063
分别为:
Figure BDA0002581600740000064
Figure BDA0002581600740000065
其中,P为目标状态向量的协方差矩阵;
生成模式j的Sigma点
Figure BDA0002581600740000066
为:
Figure BDA0002581600740000067
其中,nx为状态向量维数;γ为尺度参数;
模式j的Sigma点的量测预测
Figure BDA0002581600740000068
为:
Figure BDA0002581600740000069
模式j的量测预测
Figure BDA00025816007400000610
为:
Figure BDA00025816007400000611
其中,i表示第i个Sigma点;
Figure BDA00025816007400000612
为第i个Sigma点的均值加权所用的权值;
模式j的量测的预测协方差矩阵
Figure BDA00025816007400000613
和状态及量测之间的协方差矩阵
Figure BDA00025816007400000614
为:
Figure BDA00025816007400000615
Figure BDA00025816007400000616
其中,
Figure BDA00025816007400000617
为第i个Sigma点的协方差加权所用权值;
模式j的滤波增益为:
Figure BDA00025816007400000618
模式j的状态更新和协方差矩阵更新分别为:
Figure BDA0002581600740000071
Figure BDA0002581600740000072
可选地,所述的应用静态多模型估计器得到目标高度估计的具体方法为:
模式j的概率更新为:
Figure BDA0002581600740000073
其中,
Figure BDA0002581600740000074
N表示正态分布;
Figure BDA0002581600740000075
表示以
Figure BDA0002581600740000076
为均值、
Figure BDA0002581600740000077
为协方差矩阵的正态分布;
最终的目标状态估计结果为:
Figure BDA0002581600740000078
目标的协方差矩阵估计结果为:
Figure BDA0002581600740000079
实施本发明的多径效应下的低空目标跟踪滤波方法,具有以下有益效果:充分利用了多径效应中包含的目标位置信息,通过对静态多模型估计器中各个不同初始状态模型的单脉冲比估计误差实时比较可逐步解决高度模糊问题,减小目标高度估计的误差,最终能够对目标进行有效地跟踪。此外,结合多径场景中单脉冲比的特点在初始化方案中引入了分段二分法,使得上述方法的运行时间得以缩短,同时也减少了可能遗漏的目标状态。
附图说明
图1是本发明实施例一的多径效应下的低空目标跟踪滤波方法的流程图;
图2是本发明实施例二的多径效应下的低空目标跟踪滤波方应用的场景示意图;
图3是本发明实施例二的目标在地平线反射区的电磁波路径示意图;
图4是本发明实施例三的基于分段二分法的初始化方案的应用原理示意图,其中(a)为单脉冲比随目标高度的变化关系,(b)为路径差随目标高度的变化关系;
图5是本发明实施例三的基于分段二分法的初始化方法的示意性流程图;
图6是本发明实施例三的静态多模型估计器的框架示意图;
图7是本发明实施例三的地平线反射区目标经本发明所述方法滤波后的运动轨迹图;
图8是本发明实施例三的地平线反射区目标经本发明所述方法滤波后在四个维度下的滤波误差,图中(a)、(b)、(c)、和(d)分别对应水平距离x、水平方向速度vx、垂直高度y、以及垂直方向速度vy四个维度;
图9是本发明实施例三的地平线反射区目标分别应用常规单脉冲法、多径误差条件交互式多模型方法、以及本发明所述方法滤波后的均方根误差;
图10是本发明实施例三的本发明所述方法在不同径向距离量测噪声标准差下的性能比较。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明的一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例,都属于本发明保护的范围。
实施例一
如图1所示,本发明实施例提供的多径效应下的低空目标跟踪滤波方法一般性地可包括:
步骤S1、建立基于单脉冲比观测的系统模型,所述系统模型中的观测方程为单脉冲比关于目标状态的函数;
步骤S2、应用基于分段二分法的初始化方法进行滤波器初始化,得到目标可能的初始状态;
步骤S3、对得到的每个可能的初始状态分别应用不敏卡尔曼滤波器进行滤波;
步骤S4、将所述不敏卡尔曼滤波器的滤波结果应用静态多模型估计器得到目标高度估计,实现低空目标跟踪。
实施例二
本实施例与上面的实施例基本相同,相同之处不再赘述,不同之处在于:
本实施例对场景加以具体描述,平面地球模型下岸基雷达对海面低空目标的俯仰角测量的简化场景在图2中示出。A为岸基雷达天线,T为海面低空目标,T’为目标的镜像,R为镜面反射的反射点。θt为目标相对于地面水平线的真实俯仰角值;θm为镜像目标相对于地面水平线的俯仰角值;单脉冲和差比幅测角法下两个偏置天线的中心轴线,以下称作天线视轴,θ0为天线视轴的偏置角度。实际上来自海面的反射有两种,镜面反射和漫反射。它们的强度和海面情况有关。平静的海面有一个光滑的表面,产生的镜面反射会强一些,而漫反射很弱。而粗糙的海面漫反射就会增强一些。假设场景为平静海面,以镜面反射为主,漫反射作为量测噪声处理。路径ATA为雷达电磁波直接传播路径,路径ATRA、路径ART、路径ARTRA为三条可能的镜面反射传播路径,电磁波回波的强度同回波进入雷达天线的方向有关。
目标偏向于天线视轴的角度可以写为
θ=-(θ0t)
镜像目标偏向于天线视轴的角度可以写为
Figure BDA0002581600740000101
对于偏向天线视轴的角度,规定图2中天线视轴以下的角度为负角度。
四条路径的单脉冲比值同两条路径的单脉冲比值相同,即单脉冲比值r可以写做如下形式
Figure BDA0002581600740000102
其中,F和FΔ分别为雷达和波束与差波束的电压增益方向性函数,即天线方向图函数;ξ表征目标信号具有的特性,包括目标的RCS起伏等,对于低空目标的场景,认为直射路径和反射路径的ξ是相同的;v为海面总反射系数,表征电磁波在海面反射的特性情况,是个复数;φl是反射路径和直射路径之间的距离差造成的相位差,φl可表示为下式
Figure BDA0002581600740000103
其中,λ为雷达电磁波的波长;l为目标到雷达的距离;lr为雷达到反射点的距离;lt为目标到反射点的距离。
本发明实施例以解决反射路径在地平线反射区的低空目标跟踪问题为例,地平线反射区即直射路径信号和反射路径信号都进入了天线主波束半功率波束宽度以内的高接收增益区,如图3所示,它是低空目标跟踪问题中误差最为严重的一种情况。
所述的步骤S1中的建立基于单脉冲比观测的系统模型包括:
步骤S11、建立低空目标跟踪的状态方程:
Xk=Fk-1Xk-1+Vk-1
其中,Xk为k时刻目标状态向量,以二维平面(无方位维度)下目标跟踪为例,则目标状态向量可表示为
Figure BDA0002581600740000111
其中,xk为k时刻目标水平距离,
Figure BDA0002581600740000112
为k时刻目标水平速度,yk为k时刻目标垂直高度,
Figure BDA0002581600740000113
为k时刻目标垂直速度。
Fk-1为k-1时刻状态转移矩阵;Vk-1为k-1时刻过程噪声,所述噪声为零均值高斯白噪声,协方差矩阵为Qk-1
对于恒速度模型目标,状态转移矩阵为:
Figure BDA0002581600740000114
其中,T为雷达采样时间间隔;
步骤S12、将所述目标状态向量作为所述观测方程的输入,将单脉冲和差波束比幅雷达获得的单脉冲比量作为所述观测方程的输出,建立k时刻观测方程为:
Zk=hk(Xk)+Wk
其中,Zk为k时刻观测向量:
Zk=[ρk rk]T
其中,ρk是k时刻径向距离量测值,rk是k时刻实单脉冲比量测值;
Wk为k时刻观测噪声,所述观测噪声为零均值高斯白噪声,协方差矩阵为Rk
hk()是表征观测向量与状态向量之间关系的函数,满足下式:
hk(Xk)=[hρ(Xk) hratio(Xk)]k T
其中,hρ()是表征径向距离与状态向量之间关系的函数:
Figure BDA0002581600740000115
其中,hr为雷达高度;
hratio()是表征单脉冲比与状态向量之间关系的函数;通过以下表达式由目标状态向量得到单脉冲比以表征hratio()函数关系:
Figure BDA0002581600740000121
其中,F和FΔ分别为雷达和波束与差波束的电压增益方向性函数,即天线方向图函数;θ为目标偏向于天线视轴的角度;
Figure BDA0002581600740000122
为镜像目标偏向于天线视轴的角度;θ0为天线视轴的偏置角度;θt为目标相对于海面水平线的真实俯仰角值;θm为镜像目标相对于海面水平线的俯仰角值;v为海面总反射系数,表征电磁波在海面反射的特性情况,是个复数;φl是由反射路径与直射路径之间的距离差造成的相位差;Δρ为反射路径与直射路径的路径差;λ为电磁波波长;以上物理量在图2和图3中对应示出。
例如设定雷达高度为230m,天线发射电磁波频率6GHz,天线视轴固定0度,天线半功率波束宽度1.1度,两单脉冲天线偏离天线视轴0.55度;海面总反射系数设定为1.0e-jπ;天线方向图函数采用如下的高斯函数
Figure BDA0002581600740000123
其导函数为
Figure BDA0002581600740000124
其中,x为目标距单天线轴线的偏角;θ3dB为天线的半功率波瓣宽度。
下面的公式给出了在如上参数下hratio()的解析表达式:
Figure BDA0002581600740000131
Figure BDA0002581600740000132
Figure BDA0002581600740000133
Figure BDA0002581600740000134
Figure BDA0002581600740000135
Figure BDA0002581600740000136
实施例三
本实施例与上面的实施例二基本相同,相同之处不再赘述,不同之处在于:
如图4所示,所述的步骤S2应用基于分段二分法的初始化方法进行滤波器初始化,得到目标可能的初始状态包括具体包括以下步骤:
步骤S21、确定目标俯仰角的搜索范围:
选取天线主波束半功率波束宽度内雷达能照射到目标的区域作为所述目标俯仰角的搜索范围;
若天线主波束半功率波束宽度下边界照射到海面以下时,下边界调整为海面水平线处的俯仰角;
综上,所述目标俯仰角的搜索范围为:
Figure BDA0002581600740000137
其中,θ3dB为天线的半功率波瓣宽度;
Figure BDA0002581600740000138
为时刻1目标到雷达的径向距离估计量。
步骤S22、按照每个单脉冲比波峰波谷点对应的高度值将所述目标俯仰角的搜索范围分段,得到若干段单脉冲比的单调曲线,具体方法为:
在海面总反射系数的相位为π的假设下,当反射路径与直射路径之间的路径长度差是半波长的偶数倍时,对应为单脉冲比波峰点;当路径长度差是半波长的奇数倍时,对应为单脉冲比波谷点;由反射路径长度ρm即可得到目标此刻高度值h;
即由以下两个公式得到单脉冲比波峰点对应的高度值:
ρm=k*0.5λ+ρ,k=0,2,4,…,2n,…,n∈N
Figure BDA0002581600740000141
再由以下两个公式得到单脉冲比波谷点对应的高度值:
ρm=k*0.5λ+ρ,k=1,3,5,…,2n+1,…,n∈N
Figure BDA0002581600740000142
这里N表示自然数。
图4示出了本发明中基于二分法的初始化方案的应用原理图,展示了相同径向距离50公里处单脉冲比、路径差随目标高度的变化。参数设置同实施例二中的参数设置,不同之处在于雷达高度为100米。图4的上下两幅图中加粗曲线段是相对应的,这一段曲线就是应用二分法搜索时的其中一个研究对象,因为曲线单调的特点,故可在这段曲线上应用二分法搜索。
步骤S23、在单脉冲比的每一段单调曲线上用二分法搜索目标可能的初始状态,即:对于所述的每一段单调曲线,计算目标状态对应的单脉冲比量与单脉冲比量测量的差值,将第一个差值符合单脉冲比误差范围的目标状态作为该段单调曲线上目标可能的初始状态。
如图5所示,所述步骤S23一般性地可以包括:
步骤S231、计算目标状态对应的单脉冲比量
Figure BDA0002581600740000143
与单脉冲比量测量r的差值,将其记作f(h):
Figure BDA0002581600740000144
对于任意一个单脉冲比单调段,其两个端点对应的目标高度值分别记作hmax和hmin,且hmax>hmin
若f(hmin)*f(hmax)<0且min(|f(hmin)|,|f(hmax)|)>e,e为允许误差,则执行步骤S232;否则,执行步骤S235;
步骤S232、计算hmax、hmin的中点值hhalf
hhalf=(hmin+hmax)/2
若f(hhalf)的绝对值小于允许误差e,即|f(hhalf)|<e,则执行步骤S233;否则执行步骤S234;
步骤S233、输出hhalf对应的目标状态作为一个可能的目标初始状态,完成此单调段的搜索;
步骤S234、判断f(hhalf)*f(hmin)<0是否成立:
若成立,则令hmax=hhalf,然后返回至步骤S232,直至找到符合误差条件的目标状态,完成此单调段的搜索;
否则,令hmin=hhalf,然后返回至步骤S232,直至找到符合误差条件的目标状态,完成此单调段的搜索;返回至计算hhalf的步骤,然后重新执行后续步骤直至找到符合误差条件的目标状态,完成此单调段的搜索;
步骤S235、判断|f(hmin)|<|f(hmax)|是否成立;
若成立,则输出hmin对应的目标状态作为一个可能的目标初始状态,完成此单调段的搜索;
否则,输出hmax对应的目标状态作为一个可能的目标初始状态,完成此单调段的搜索。
所述的步骤S4中的静态多模型估计器,其特点是把每一个可能的目标初始状态向量对应的目标运动轨迹看作是系统的一种工作模式,系统虽然可能有多种工作模式,但系统的实际工作模式只会是众多模式中的一种,并维持不变,即工作模式不随时间切换,即:
Figure BDA0002581600740000161
其中,P为概率密度;k代表采样时间;ξk为k时刻工作模式,ξk∈M={m1,m2,…,mz};M为工作模式集合,mi为第i个工作模式;z为工作模式的个数。
由于观测量同状态向量间的复杂的非线性关系,基于单脉冲比观测的跟踪滤波方法应选取非线性卡尔曼滤波方法,步骤S2的初始化方法将得到多个可能的目标初始状态,对于每一个单调段得到的一个目标初始状态,它都被认为是潜在的目标实际状态与高斯的观测噪声相加后的结果,因而可假定对应一个单调段的目标状态的后验概率密度是高斯分布的,所以对每一个单调段都将可应用不敏卡尔曼滤波器进行滤波;
步骤S3、对得到的目标状态向量中的每个目标初始状态值应用不敏卡尔曼滤波器进行滤波,具体方法为:
将每个可能的初始状态作为一个模式;
对各模式分别进行不敏卡尔曼滤波
Figure BDA0002581600740000162
模式j的状态一步预测
Figure BDA0002581600740000163
和一步预测协方差矩阵
Figure BDA0002581600740000164
分别为:
Figure BDA0002581600740000165
Figure BDA0002581600740000166
其中,P为目标状态向量的协方差矩阵;
生成模式j的Sigma点
Figure BDA0002581600740000167
为:
Figure BDA0002581600740000168
其中,nx为状态向量维数;γ为尺度参数;
模式j的Sigma点的量测预测
Figure BDA0002581600740000169
为:
Figure BDA00025816007400001610
模式j的量测预测
Figure BDA00025816007400001611
为:
Figure BDA0002581600740000171
其中,i表示第i个Sigma点;
Figure BDA0002581600740000172
为第i个Sigma点的均值加权所用的权值;
模式j的量测的预测协方差矩阵
Figure BDA0002581600740000173
和状态及量测之间的协方差矩阵
Figure BDA0002581600740000174
为:
Figure BDA0002581600740000175
Figure BDA0002581600740000176
其中,
Figure BDA0002581600740000177
为协方差加权所用权值;
模式j的滤波增益为:
Figure BDA0002581600740000178
模式j的状态更新和协方差矩阵更新分别为:
Figure BDA0002581600740000179
Figure BDA00025816007400001710
步骤S4、应用静态多模型估计器得到目标高度估计的具体方法为:
模式j的概率更新为:
Figure BDA00025816007400001711
其中,
Figure BDA00025816007400001712
N表示正态分布;
Figure BDA00025816007400001713
表示以
Figure BDA00025816007400001714
为均值,
Figure BDA00025816007400001715
为协方差矩阵的正态分布;
最终的目标状态估计结果为:
Figure BDA00025816007400001716
目标的协方差矩阵估计结果为:
Figure BDA00025816007400001717
本发明所述方法中步骤S3和S4涉及的静态多模型估计器的一种流程框架在图6中示出。
下面举例说明本发明实施例所述方法的效果。
采用单频率单脉冲雷达对某匀速巡航运动目标进行跟踪;雷达参数设置:雷达高度为100m,天线视轴固定0度,天线半功率波束宽度1.1度,两单脉冲天线偏离天线视轴0.55度;天线发射电磁波频率为6GHz,雷达采样数据率50Hz;目标从距雷达50km、距海面40m的高度处朝向雷达飞行;目标水平方向速度-300m/s,目标垂直方向速度为零,定义目标水平速度以远离雷达的情况为正值,速度在垂直方向上远离海面定义为正值;目标的整个运动过程都在天线主波束半功率波束宽度以内;目标运动过程中噪声标准差[σx,σy]=[1m/s2,0.1m/s2];径向距离量测噪声标准差为20m,单脉冲比量测噪声标准差为7e-6。由常规单脉冲法计算得到的目标高度值可看做是目标量测值,以此作为多径误差条件交互式多模型方法的量测量;而本发明方法的量测量是由常规单脉冲法得到的单脉冲比加上量测噪声得到。
常规单脉冲方法、多径误差条件交互式多模型方法(即多径误差条件IMM方法)、以及本发明方法的滤波效果如图7所示,小图是将前500步的跟踪结果放大呈现;其中本发明方法在目标状态四个维度的滤波误差如图8所示,图8中(a)、(b)、(c)、和(d)分别对应水平距离x、水平方向速度vx、垂直高度y、以及垂直方向速度vy四个维度。可以看出,本发明的方法估计出的目标轨迹逐渐向目标实际运动过程收敛;四个维度的滤波误差都随着跟踪步数增加而逐渐减小;常规单脉冲方法推导得到的目标高度会出现近似周期性尖峰,但本发明利用了单脉冲比构造目标状态函数,因此可以避免常规单脉冲方法在海面低空多径下的失效;而多径误差条件交互式多模型方法对于尖峰测量误差有很好的滤波效果,但始终有固定偏差没有消除,对固定偏差的滤波效果较差。
通过大量蒙特卡洛仿真实验计算各方法下目标高度估计的均方根误差,结果如图9所示。均方根误差定义为
Figure BDA0002581600740000191
其中,h为目标高度实际值;
Figure BDA0002581600740000192
为方法估计的目标高度值;N为蒙特卡洛仿真实验的次数。在图9中展示了多径误差条件交互式多模型方法和本发明所述方法滤波后的均方根误差,参数设置与图7基本相同,不同之处在于径向距离量测噪声标准差为0.1m,单脉冲比比值量测噪声标准差为7e-6;目标从距雷达43km、距海面40m的高度处朝向雷达水平巡航飞行70s;进行20次蒙特卡洛仿真实验。
从图9的小图中可以看到,本发明方法可以在跟踪一定时间后,RMSE趋于零,即根据多个时刻的量测信息以及目标运动模型估计等先验信息,逐步减小了高度模糊的影响,对目标跟踪的效果好于仅对尖峰测量误差有很好的滤波效果的多径误差条件交互式多模型方法。在34km至更近的距离中出现许多相距很近有重合的尖峰,不同于单次滤波的规律,这是因为每次蒙特卡洛仿真实验中加入了不同的目标运动过程噪声,对于超低空目标的位置有很小的不同就可能导致出现尖峰时刻的位置有大的改变。
对不同径向距离量测噪声标准差下的本发明方法的性能进行了比较,结果如图10所示,通过蒙特卡洛仿真得到目标高度估计的均方根误差。参数设置与图7基本相同,不同之处在于图10实验了不同标准差下的径向距离量测噪声,每一个径向距离量测标准差下都进行20次蒙特卡洛仿真实验。可以看到,整体趋势上目标高度估计的均方根误差随径向距离量测噪声标准差增大而增大;在径向距离量测噪声标准差较小时,本发明所述方法对高度模糊问题解决的效果较好;径向距离量测噪声标准差越小,用本发明所述方法去消除高度模糊问题所需要的跟踪步数相对越少。
综上所述,本发明提供了一种多径效应下的低空目标跟踪滤波方法,所述方法充分利用了多径效应中包含的目标位置信息,通过各模型单脉冲比估计误差实时比较逐步解决高度模糊问题,减小目标高度估计的误差,最终实现对目标的有效地跟踪。此外,本发明还结合多径场景中单脉冲比的特点在初始化方案中引入了分段二分法,使得上述方法的运行时间得以缩短,同时能够减少可能遗漏的目标状态。
最后应说明的是:以上实施例仅用以说明本发明的技术方案,而非对其限制;尽管参照前述实施例对本发明进行了详细的说明,本领域的普通技术人员应当理解:其依然可以对前述各实施例所记载的技术方案进行修改,或者对其中部分技术特征进行等同替换;而这些修改或者替换,并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (1)

1.一种多径效应下的低空目标跟踪滤波方法,其特征在于,包括:
建立基于单脉冲比观测的系统模型,所述系统模型中的观测方程为单脉冲比关于目标状态的函数;
应用基于分段二分法的初始化方法进行滤波器初始化,得到目标可能的初始状态;
对得到的每个可能的初始状态分别应用不敏卡尔曼滤波器进行滤波;
将所述不敏卡尔曼滤波器的滤波结果应用静态多模型估计器得到目标高度估计,实现低空目标跟踪;
所述的建立基于单脉冲比观测的系统模型包括:
建立低空目标跟踪的状态方程:
Xk=Fk-1Xk-1+Vk-1
其中,Xk为k时刻目标状态向量,Fk-1为k-1时刻状态转移矩阵;Vk-1为k-1时刻过程噪声,所述噪声为零均值高斯白噪声,协方差矩阵为Qk-1
对于恒速度模型目标,状态转移矩阵为:
Figure FDA0003647282740000011
其中,T为雷达采样时间间隔;
将所述目标状态向量作为所述观测方程的输入,将单脉冲和差波束比幅雷达获得的单脉冲比量作为所述观测方程的输出,建立k时刻观测方程为:
Zk=hk(Xk)+Wk
其中,Zk为k时刻观测向量:
Zk=[ρk rk]T
其中,ρk是k时刻径向距离量测值,rk是k时刻实单脉冲比量测值;
Wk为k时刻观测噪声,所述观测噪声为零均值高斯白噪声,协方差矩阵为Rk
hk()满足下式:
hk(Xk)=[hρ(Xk) hratio(Xk)]k T
其中,
Figure FDA0003647282740000021
其中,hr为雷达高度,xk为k时刻目标水平距离,yk为k时刻目标垂直高度;
hratio()是表征单脉冲比量测量与状态向量之间关系的函数;
所述的hratio()的获得方法包括:通过以下表达式由目标状态向量得到单脉冲比以表征此函数关系:
Figure FDA0003647282740000022
其中,FΣ和FΔ分别为雷达和波束与差波束的电压增益方向性函数,即天线方向图函数;θ为目标偏向于天线视轴的角度;
Figure FDA0003647282740000023
为镜像目标偏向于天线视轴的角度;θ0为天线视轴的偏置角度;θt为目标相对于海面水平线的真实俯仰角值;θm为镜像目标相对于海面水平线的俯仰角值;ν为海面总反射系数;Δρ为反射路径与直射路径之间的路径长度差;φl是由反射路径与直射路径之间的路径长度差造成的相位差;λ为电磁波波长;
所述的应用基于分段二分法的初始化方法进行滤波器初始化,得到目标可能的初始状态包括:
确定目标俯仰角的搜索范围;
按照每个单脉冲比波峰波谷点对应的高度值将所述目标俯仰角的搜索范围分段,得到若干段单脉冲比的单调曲线;
在单脉冲比的每一段单调曲线上用二分法搜索目标可能的初始状态,即:对于所述的每一段单调曲线,计算目标状态对应的单脉冲比量与单脉冲比量测量的差值,将第一个差值符合单脉冲比误差范围的目标状态作为该段单调曲线上目标可能的初始状态;
所述的确定目标俯仰角的搜索范围包括:
选取天线主波束半功率波束宽度内雷达能照射到目标的区域作为所述目标俯仰角的搜索范围;
若天线主波束半功率波束宽度下边界照射到海面以下,则将所述下边界调整为海面水平线处的俯仰角;
所述的单脉冲比波峰波谷点对应的高度值的确定方法包括:
当反射路径与直射路径之间的路径长度差是半波长的偶数倍时,对应为单脉冲比波峰点;当路径长度差是半波长的奇数倍时,对应为单脉冲比波谷点;由反射路径长度ρm即可得到目标此刻高度值h;
所述的在单脉冲比的每一段单调曲线上用二分法进行搜索目标可能的初始状态的具体方法为:
步骤S231、计算目标状态对应的单脉冲比量
Figure FDA0003647282740000031
与单脉冲比量测量r的差值,将其记作f(h):
Figure FDA0003647282740000032
对于任意一个单脉冲比单调段,其两个端点对应的目标高度值分别记作hmax和hmin,且hmax>hmin
若f(hmin)*f(hmax)<0且min(|f(hmin)|,|f(hmax)|)>e,e为允许误差,则执行步骤S232;否则,执行步骤S235;
步骤S232、计算hmax、hmin的中点值hhalf
hhalf=(hmin+hmax)/2
若|f(hhalf)|<e,则执行步骤S233;否则执行步骤S234;
步骤S233、输出hhalf对应的目标状态作为一个可能的目标初始状态,完成此单调段的搜索;
步骤S234、判断f(hhalf)*f(hmin)<0是否成立:
若成立,则令hmax=hhalf,然后返回至步骤S232,直至找到符合误差条件的目标状态,完成此单调段的搜索;
否则,令hmin=hhalf,然后返回至步骤S232,直至找到符合误差条件的目标状态,完成此单调段的搜索;返回至计算hhalf的步骤,然后重新执行后续步骤直至找到符合误差条件的目标状态,完成此单调段的搜索;
步骤S235、判断|f(hmin)|<|f(hmax)|是否成立;
若成立,则输出hmin对应的目标状态作为一个可能的目标初始状态,完成此单调段的搜索;
否则,输出hmax对应的目标状态作为一个可能的目标初始状态,完成此单调段的搜索;
所述的静态多模型估计器满足如下公式:
Figure FDA0003647282740000041
其中,P为概率密度;k代表采样时间;ξk为k时刻工作模式,ξk∈M={m1,m2,…,mz};M为工作模式集合,mi为第i个工作模式;z为工作模式的个数;
所述的对得到的每个可能的初始状态分别应用不敏卡尔曼滤波器进行滤波的具体方法为:
将每个可能的初始状态作为一个模式;
对各模式分别进行不敏卡尔曼滤波
Figure FDA0003647282740000045
模式j的状态一步预测
Figure FDA0003647282740000042
和一步预测协方差矩阵
Figure FDA0003647282740000043
分别为:
Figure FDA0003647282740000044
Figure FDA0003647282740000051
其中,P为目标状态向量的协方差矩阵;
生成模式j的Sigma点
Figure FDA0003647282740000052
为:
Figure FDA0003647282740000053
其中,nx为状态向量维数;γ为尺度参数;
模式j的Sigma点的量测预测
Figure FDA0003647282740000054
为:
Figure FDA0003647282740000055
模式j的量测预测
Figure FDA0003647282740000056
为:
Figure FDA0003647282740000057
其中,i表示第i个Sigma点;
Figure FDA0003647282740000058
为第i个Sigma点的均值加权所用的权值;
模式j的量测的预测协方差矩阵
Figure FDA0003647282740000059
和状态及量测之间的协方差矩阵
Figure FDA00036472827400000510
为:
Figure FDA00036472827400000511
Figure FDA00036472827400000512
其中,
Figure FDA00036472827400000513
为第i个Sigma点的协方差加权所用权值;
模式j的滤波增益为:
Figure FDA00036472827400000514
模式j的状态更新和协方差矩阵更新分别为:
Figure FDA00036472827400000515
Figure FDA00036472827400000516
所述的应用静态多模型估计器得到目标高度估计的具体方法为:
模式j的概率更新为:
Figure FDA00036472827400000517
其中,
Figure FDA00036472827400000518
N表示正态分布;
Figure FDA00036472827400000519
表示以
Figure FDA0003647282740000061
为均值、
Figure FDA0003647282740000062
为协方差矩阵的正态分布;最终的目标状态估计结果为:
Figure FDA0003647282740000063
目标的协方差矩阵估计结果为:
Figure FDA0003647282740000064
CN202010671798.8A 2020-07-13 2020-07-13 一种多径效应下的低空目标跟踪滤波方法 Active CN111708015B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010671798.8A CN111708015B (zh) 2020-07-13 2020-07-13 一种多径效应下的低空目标跟踪滤波方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010671798.8A CN111708015B (zh) 2020-07-13 2020-07-13 一种多径效应下的低空目标跟踪滤波方法

Publications (2)

Publication Number Publication Date
CN111708015A CN111708015A (zh) 2020-09-25
CN111708015B true CN111708015B (zh) 2022-06-21

Family

ID=72546313

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010671798.8A Active CN111708015B (zh) 2020-07-13 2020-07-13 一种多径效应下的低空目标跟踪滤波方法

Country Status (1)

Country Link
CN (1) CN111708015B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112529797A (zh) * 2020-12-04 2021-03-19 中国人民解放军63921部队 基于序列视轴指向矢量的目标轨迹确认方法
CN114417942B (zh) * 2022-03-28 2022-06-07 成都数之联科技股份有限公司 一种杂波识别方法及系统及装置及介质

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6078281A (en) * 1996-06-28 2000-06-20 Milkovich Systems Engineering Signal processing architecture which improves sonar and pulse Doppler radar performance and tracking capability
JP2011080794A (ja) * 2009-10-05 2011-04-21 Mitsubishi Electric Corp パルスレーダ装置
CN104777478A (zh) * 2015-04-16 2015-07-15 电子科技大学 一种相控阵雷达搜索捕获目标的方法
CN107526070A (zh) * 2017-10-18 2017-12-29 中国航空无线电电子研究所 天波超视距雷达的多路径融合多目标跟踪算法
CN108802721A (zh) * 2018-08-22 2018-11-13 哈尔滨工业大学 一种任意直线约束下目标跟踪方法
CN110716196A (zh) * 2019-11-04 2020-01-21 广东中科四创科技有限公司 一种多点低慢小空中目标跟踪临视系统
CN111273278A (zh) * 2020-02-06 2020-06-12 零八一电子集团有限公司 四通道毫米波数字和差单脉冲精密跟踪系统

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7046188B2 (en) * 2003-08-14 2006-05-16 Raytheon Company System and method for tracking beam-aspect targets with combined Kalman and particle filters
US11269047B2 (en) * 2017-12-06 2022-03-08 Invensense, Inc. Three dimensional object-localization and tracking using ultrasonic pulses with synchronized inertial position determination

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6078281A (en) * 1996-06-28 2000-06-20 Milkovich Systems Engineering Signal processing architecture which improves sonar and pulse Doppler radar performance and tracking capability
JP2011080794A (ja) * 2009-10-05 2011-04-21 Mitsubishi Electric Corp パルスレーダ装置
CN104777478A (zh) * 2015-04-16 2015-07-15 电子科技大学 一种相控阵雷达搜索捕获目标的方法
CN107526070A (zh) * 2017-10-18 2017-12-29 中国航空无线电电子研究所 天波超视距雷达的多路径融合多目标跟踪算法
CN108802721A (zh) * 2018-08-22 2018-11-13 哈尔滨工业大学 一种任意直线约束下目标跟踪方法
CN110716196A (zh) * 2019-11-04 2020-01-21 广东中科四创科技有限公司 一种多点低慢小空中目标跟踪临视系统
CN111273278A (zh) * 2020-02-06 2020-06-12 零八一电子集团有限公司 四通道毫米波数字和差单脉冲精密跟踪系统

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
"A fuzzy approach for tracking of low-altitude target in the presence of multipath propagation";Y.M. Chen 等;《Proceedings of the Fifth International Conference on Information Fusion》;20021231;第143-148页 *
"Route-Based Dynamics Modeling and Tracking With Application to Air Traffic Surveillance";Linfeng Xu 等;《IEEE Transactions on Intelligent Transportation Systems 》;20190116;第209-221页 *
"Tracking low elevation targets in the presence of multipath propagation";Bar-Shalom Y 等;《IEEE Transactions on Aerospace and IEEE Transactions on Aerospace and Electronic Systems》;19941231;第973-979页 *
"海面多径环境下雷达目标俯仰角测量提取的研究与应用";吕韶昱 等;《国防科技大学学报》;20071015;第48-53页 *

Also Published As

Publication number Publication date
CN111708015A (zh) 2020-09-25

Similar Documents

Publication Publication Date Title
CN107976660B (zh) 弹载多通道雷达超低空目标分析与多径回波建模方法
CN111708015B (zh) 一种多径效应下的低空目标跟踪滤波方法
CN107607943B (zh) 基于干涉相位辅助的延迟多普勒雷达高度表的测高方法
CN109655802B (zh) 一种基于clean算法的多目标粒子群长时间积累检测方法
CN108490443B (zh) 基于解析解及NUFFT的多子阵合成孔径声纳ωk成像算法
CN113391315B (zh) 基于电磁波抛物方程伴随模式的雷达回波资料反演大气波导的方法
CN105445718A (zh) 一种基于阵列重构的分布式多载舰超视距雷达的doa估计方法
Lapierre et al. New methods for handling the range dependence of the clutter spectrum in non-sidelooking monostatic STAP radars
EP3418770B1 (en) Synthetic-aperture-radar signal processing device
CN114545411A (zh) 一种基于工程实现的极坐标格式多模高分辨sar成像方法
CN112612006A (zh) 基于深度学习的机载雷达非均匀杂波抑制方法
CN111856466A (zh) 一种高效的复杂运动目标isar平动补偿方法
CN108107432B (zh) 基于时域扰动的高低轨双基sar保相成像方法
CN110879391B (zh) 基于电磁仿真和弹载回波仿真的雷达图像数据集制作方法
CN110261837B (zh) 一种基于航迹信息的复杂目标rcs计算方法
CN112415512B (zh) 基于进退法和黄金分割法的sar运动目标聚焦方法
CN109633585B (zh) 分布式机会阵雷达非合作目标动态回波的高精度计算方法
CN113376625B (zh) 目标物体的偏离角度获得方法、装置、电子设备及存储介质
CN113671485B (zh) 基于admm的米波面阵雷达二维doa估计方法
CN114442061A (zh) 一种基于距离选通和交替反演的折叠杂波的抑制方法
CN108983192B (zh) 基于gps辐射源的雷达运动目标参数估计方法
CN112835006A (zh) 一种基于帧间积累的跟踪雷达海上小目标检测方法及系统
RU2761955C1 (ru) Способ определения высоты полета низколетающей цели моноимпульсной РЛС сопровождения
CN117491988B (zh) 一种粒子滤波宽带多频低空测角方法
Machhour et al. Synthetic Sea-Clutter Modelling for STAP

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