CN105785358B - 一种方向余弦坐标系下带多普勒量测的雷达目标跟踪方法 - Google Patents
一种方向余弦坐标系下带多普勒量测的雷达目标跟踪方法 Download PDFInfo
- Publication number
- CN105785358B CN105785358B CN201610339346.3A CN201610339346A CN105785358B CN 105785358 B CN105785358 B CN 105785358B CN 201610339346 A CN201610339346 A CN 201610339346A CN 105785358 B CN105785358 B CN 105785358B
- Authority
- CN
- China
- Prior art keywords
- eta
- gamma
- centerdot
- sigma
- pseudo
- 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
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
-
- 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/02—Systems using reflection of radio waves, e.g. primary radar systems; Analogous systems
- G01S13/50—Systems of measurement based on relative movement of target
- G01S13/58—Velocity or trajectory determination systems; Sense-of-movement determination systems
- G01S13/589—Velocity or trajectory determination systems; Sense-of-movement determination systems measuring the velocity vector
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
本发明涉及一种方向余弦坐标系下带多普勒量测的雷达目标跟踪方法,包括以下步骤:伪量测构造步骤,通过当前时刻k雷达获得的距离量测和多普勒量测的乘积构造伪量测;量测转换步骤,将方向余弦坐标系下的位置量测转换到直角坐标系下,获得转换位置量测;无偏一二阶矩计算步骤,计算转换位置量测误差的无偏一二阶矩和伪量测误差的无偏一二阶矩;笛卡尔状态信息提取步骤,利用转换位置量测和转换位置量测误差的无偏一二阶矩提取目标的笛卡尔状态信息;伪状态空间构造及伪状态信息提取步骤,利用伪量测的真值及其导数构造伪状态空间,并利用伪量测和笛卡尔状态信息提取伪状态信息;静态融合步骤,对伪状态信息和笛卡尔状态信息进行静态融合。
Description
技术领域
本发明涉及雷达目标跟踪,尤其涉及方向余弦坐标系中的雷达目标跟踪。
背景技术
雷达目标跟踪就是根据目标当前状态和雷达量测,对目标运动状态进行估计和预测。在目标跟踪领域,雷达提供的量测除了位置(距离和角度)外,还有多普勒信息,理论与实践已经证明,充分利用多普勒量测信息可以有效地提高目标跟踪精度。
由于传统机械扫描雷达受限于扫描速度,而相控阵雷达中的方向图和扫描特征又可以很方便的在方向余弦坐标系(COS)下表示,相控阵雷达的应用越来越广泛。但相控阵雷达COS坐标系下目标跟踪的研究相对较少,考虑多普勒量测的跟踪算法就更少了。现有文献中Y.Kosuge,H.Iwama and Y.Miyazaki,“A tracking filter for phased array radarwith range rate measurement,”Proceedings of 1991International Conference onIndustrial Electronics,Control and Instrumentation,pp.2555-2560,1991直接利用扩展卡尔曼滤波器(EKF)对COS坐标系下的位置量测和多普勒量测进行跟踪滤波,但由于多普勒量测的强非线性特点,加之COS坐标系下第三个方向余弦是前两个方向余弦的强非线性函数,在大的量测误差情况下,EKF的近似误差在滤波迭代中容易累计扩大,从而使滤波器性能恶化,会出现精度下降和滤波发散问题。文献B.Zhang,H.Qu and S.Li,“A newmethod for target tracking with debiased consisitent converted measurementsin direction cosines,”Chinese Journal of Electronics,Vol.19,no.3,pp.538-542,2010利用转换量测的方法对方向余弦坐标下的目标进行跟踪,但没有考虑包含目标速度信息的多普勒量测,属于次优的方案。
发明内容
本发明鉴于背景技术的以上问题提出,用于解决背景技术中存在的问题,至少是提供一种有益的选择。
为了实现以上目的,一种带多普勒量测的雷达目标跟踪方法,包括以下步骤:伪量测构造步骤,通过当前时刻k雷达获得的距离量测和多普勒量测的乘积构造伪量测;量测转换步骤,将方向余弦坐标系下的位置量测转换到直角坐标系下,获得转换位置量测;无偏一二阶矩计算步骤,利用伪量测和转换位置量测计算转换位置量测误差的无偏一二阶矩和转换多普勒量测误差的无偏一二阶矩;笛卡尔状态信息提取步骤,利用转换位置量测和转换位置量测误差的无偏一二阶矩提取目标的笛卡尔状态信息;伪状态空间构造及伪状态信息提取步骤,利用伪量测的真值(转换多普勒)及其导数构造伪状态空间,并利用伪量测和笛卡尔状态信息提取伪状态信息;静态融合步骤,对伪状态信息和笛卡尔状态信息进行静态融合。
根据本发明的一些实施方式,针将多普勒量测引入方向余弦坐标系下进行目标跟踪,在静态融合滤波器的框架下,推导了方向余弦坐标系下的量测转换和量测转换误差的无偏一二阶矩,并同时利用转换多普勒量测卡尔曼滤波器(CDMKF)和转换位置量测卡尔曼滤波器(CPMKF)线性的提取多普勒信息和目标笛卡尔状态信息,然后联合两者的输出在最小均方误差(MMSE)准则下估计出目标最终状态,从而能够提高目标跟踪的精度。
附图说明
结合附图,可以更好地理解本发明,但是附图仅仅是示例性的,不是对本发明的限制。
图1示出了方向余弦坐标系和直角坐标系的关系。
图2示出了一种本发明一种实施方式的带多普勒量测的雷达目标跟踪方法的示例流程。
图3示出了仿真情况下的RRMSE位置误差图。
图4示出了仿真情况下的RRMSE速度误差图。
具体实施方式
下面结合附图具体说明本发明的实施方式。所说明的实施方式仅是示例性的,不是对本发明的限制。在所做的说明中,各实施方式可以互相参照。
在陈述发明步骤之前,先介绍一下COS坐标系下带多普勒量测的目标跟踪的基本数学模型。
带多普勒量测的目标跟踪模型在直角坐标系中用离散时间状态方程表示为
X(k+1)=ΦX(k)+ΓV(k) (42)
其中,为目标运动状态,x(k),y(k)和z(k)分别为目标target在x,y和z方向上的三个位置分量,和为相应的速度分量,Φ,Γ分别为状态转移矩阵和过程噪声增益矩阵,V(k)是均值为0,方差为Q(k)的高斯过程噪声。
图1示出了方向余弦坐标系和直角坐标系的关系。如图1所示,量测方程可表示为
Zm(k)=f[X(k)]+W(k) (43)
其中
rm(k),αm(k),βm(k)和分别为径向距离、两个方向余弦和多普勒量测,r(k),α(k),β(k)和为相应的真值,和为相应的均值为0的高斯量测噪声,方差分别为和且和互不相关,和互不相关,和的相关系数为ρ。
方向余弦坐标系中带多普勒量测的雷达目标跟踪的目的,就是根据k时刻相控阵雷达对于目标的量测rm(k),αm(k),βm(k)和以及先验的量测偏差信息均值为0、方差分别为和的高斯白噪声和和的相关系数ρ,估计出目标当前时刻的运动状态
图2示出了一种本发明一种实施方式的带多普勒量测的雷达目标跟踪方法的示例流程。如图2所示,首先在步骤一,S101:通过当前时刻k雷达获得的距离量测rm(k)和多普勒量测的乘积构造伪量测
其中是直角坐标系中伪量测ηc(k)的转换误差。伪量测的真值为转换多普勒。
下标m是measurement的首字母,表明是测量值;下标c是convert的首字母,表明是转换值,在方向余弦和直角坐标系中,伪量测的数学形式一样,为了统一后面转换量测(位置、多普勒)卡尔曼滤波器的数学形式,统一用c表明转换量测值;η(k)是转换多普勒,是伪量测对应的真值。
然后在步骤二S102进行量测转换,将方向余弦坐标系下的量测转换到直角坐标系下。在一种实施方式中,可以如下进行
其中,xc(k),yc(k)和zc(k)分别为直角坐标系中x,y和z方向上的转换位置量测,和分别是直角坐标系中相应的转换位置量测误差,rm(k),αm(k)和βm(k)分别是当前时刻k雷达获得的距离量测和两个方向余弦量测,其中第三个方向余弦γm(k)为
接着在步骤三S103,计算转换位置量测误差和转换多普勒量测误差的无偏一二阶矩。转换位置量测误差和转换多普勒量测误差的均值和方差依次为(为简化起见,部分变量的索引时刻k给予省略)
其中
其中,rm(k),αm(k),βm(k)和分别是当前时刻k雷达获得的距离量测、两个方向余弦量测和多普勒量测,σr,σα,σβ和分别是距离量测、两个方向余弦和多普勒量测的测量偏差。ρ是距离和多普勒量测之间的相关系数。γm(k)同步骤二S102。
步骤四S104:提取目标的笛卡尔状态信息,在一种实施方式中,利用CPMKF来进行提取,其迭代过程如下
Pp(k+1,k+1)=[I-Kp(k+1)Hp]Pp(k+1,k) (71)
其中
该部分输出的笛卡尔状态信息即为和
步骤五S105:通过当前时刻转换多普勒η(k)及其导数构造伪状态空间,并利用CDMKF提取伪状态信息。
构造伪状态空间为
CDMKF的迭代过程如下
Pη(k+1,k+1)=[I-Kη(k+1)Hη]Pη(k+1,k) (78)
其中
其中T是雷达扫描周期,q是直角坐标系中各个坐标轴方向的过程高斯白噪声的方差,式(34)中Pp(k,k)由步骤五中的CPMKF提供。
步骤六S106:静态融合步骤四和步骤五的输出结果(为简化起见,部分变量的索引时刻k给予省略)。
1)计算伪状态估计和目标位置估计之间的互协方差
其中
2)计算目标状态和伪观测状态(将伪状态η(k)当做目标最终状态的一种观测状态,伪状态是目标最终状态的一个数学函数)之间的协方差
其中C是伪状态与目标状态之间的函数关系,定义如下
是函数C的Jacobin矩阵。
3)计算伪观测状态的方差
其中,ei是直角坐标系中第i个nη维偏置单位向量,是函数C的Jacobin矩阵,为函数C的第i个分量的Hessian矩阵。
4)计算目标的最终状态和状态估计方差
P=Pp-PXZ(PZZ)-1(PXZ)T (87)
其中
本发明的一些实施方式相对于一些其他方法的优势在于,将COS坐标系下的动态非线性估计问题转换成一个动态线性估计问题和一个静态非线性融合问题,避免了利用非线性滤波方法EKF直接同时处理位置量测和多普勒量测可能出现的精度下降和滤波发散问题,从而精确稳健的估计COS坐标系下的运动目标状态。
为了验证方向余弦坐标系下静态融合无偏转换量测卡尔曼滤波器的有效性,将本文算法(SF-CMKFcos)与仅考虑位置量测的CPMKF算法、同时处理位置和多普勒量测的SEKF算法和UKF算法进行仿真比较。
仿真情况设定相控阵雷达位于坐标原点,以1s的扫描间隔给出目标的斜距、两个方向余弦和多普勒量测信息,量测的标准差分别为σr=1000m,σα=σβ=0.01和过程噪声方差为q=0.01m/s2,目标做匀速运动,初始位置为(30km,30km,30km),初始速度为(20m/s,20m/s,20m/s)。评价指标为位置、速度的相对均方根误差(RRMSE),其定义为
对上述条件做500次蒙特卡洛实验的100次跟踪扫描Monte-Carlo仿真结果如图3和图4所示。
从仿真结果可以看出,带多普勒量测的三种滤波器(SEKF、UKF和SF-CMKFcos)的RMSE明显小于不带多普勒量测的CPMKF滤波器,这说明多普勒量测的引入,可以显著改善跟踪滤波器的性能;而SF-CMKFcos性能最优。这是因为大的量测误差导致了大的非线性近似误差,并且通过SEKF和UKF的动态非线性迭代引起了滤波器性能恶化。而SF-CMKFcos中采用两个线性最优滤波器(CPMKF和CDMKF)处理转换量测,大的非线性近似误差只是用来更新静态融合中的权系数,没有传播到下一步迭代过程中,因而保证了滤波器性能。
以上的实施方式都是示例性的,不是对本发明的限制,本领域技术获益于本发明作出的对本发明的各种变换和改进也在本发明的保护范围内。
Claims (4)
1.一种方向余弦坐标系下带多普勒量测的雷达目标跟踪方法,包括以下步骤:
伪量测构造步骤,通过当前时刻k雷达获得的距离量测rm(k)和多普勒量测的乘积构造伪量测;
量测转换步骤,将方向余弦坐标系下的位置量测转换到直角坐标系下,获得转换位置量测;
无偏一二阶矩计算步骤,利用所述伪量测和所述转换位置量测计算转换位置量测误差的无偏一二阶矩和转换多普勒量测误差的无偏一二阶矩;
笛卡尔状态信息提取步骤,利用所述转换位置量测和所述转换位置量测误差的无偏一二阶矩提取目标的笛卡尔状态信息;
伪状态空间构造及伪状态信息提取步骤,利用伪量测的真值及其导数构造伪状态空间,并利用伪量测和笛卡尔状态信息提取伪状态信息;
静态融合步骤,对所述伪状态信息和所述笛卡尔状态信息提取步骤所提取的笛卡尔状态信息进行静态融合;
在所述伪量测构造步骤中,利用转换多普勒量测卡尔曼滤波器提取伪状态信息,以及在笛卡尔状态信息提取步骤中利用转换位置量测卡尔曼滤波器线性地提取目标笛卡尔状态信息;
在所述伪量测构造步骤中,根据以下公式构造伪量测
其中是直角坐标系中伪量测ηc(k)的转换误差,rm(k)和分别为径向距离和多普勒量测,下标m表明是测量值;下标c表明是转换值,η(k)是转换多普勒,是伪量测对应的真值;
在所述量测转换步骤,依据以下公式将方向余弦坐标系下的量测转换到直角坐标系下
其中,rm(k),αm(k),βm(k)和分别为径向距离、两个方向余弦和多普勒量测,xc(k),yc(k)和zc(k)分别为直角坐标系中x,y和z方向上的转换后的位置量测,和分别是直角坐标系中相应的转换位置量测误差,rm(k),αm(k)和βm(k)分别是当前时刻k雷达获得的距离量测和两个方向余弦量测,其中第三个方向余弦γm(k)为
在所述无偏一二阶矩计算步骤中,利用以下公式计算转换位置量测误差和转换多普勒量测误差的无偏一二阶矩;
其中转换位置量测误差和转换多普勒量测误差的均值和方差依次为
其中
其中,σr,σα,σβ和分别是距离量测、两个方向余弦和多普勒量测的测量偏差。ρ是距离和多普勒量测之间的相关系数。
2.根据权利要求1所述的方法,其特征在于,在所述笛卡尔状态信息提取步骤中,由转换位置量测卡尔曼滤波器提取目标的笛卡尔状态信息,其迭代过程如下
Pp(k+1,k+1)=[I-Kp(k+1)Hp]Pp(k+1,k) (25)
其中
3.根据权利要求2所述的方法,其特征在于,在所述伪状态空间构造步骤,构造伪状态空间为
其中,利用转换多普勒量测卡尔曼滤波器的迭代提取伪状态信息,过程如下
Pη(k+1,k+1)=[I-Kη(k+1)Hη]Pη(k+1,k) (32)
其中
其中T是雷达扫描周期,q是直角坐标系中各个坐标轴方向的过程高斯白噪声的方差,式(34)中Pp(k,k)是由笛卡尔状态信息提取步骤提取的笛卡尔状态信息。
4.根据权利要求3所述的方法,其特征在于,在所述静态融合步骤中,采用以下的公式进行静态融合:
1)根据以下公式进行计算
其中
2)再根据以下公式进行计算
其中C是伪状态与目标状态之间的函数关系,定义如下
是函数C的Jacobin矩阵;
3)接着根据以下公式进行计算
其中,ei是直角坐标系中第i个nη维偏置单位向量,是函数C的Jacobin矩阵,为函数C的第i个分量的Hessian矩阵;
4)计算目标的最终状态和状态估计方差
P=Pp-PXZ(PZZ)-1(PXZ)T (41)
其中
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610339346.3A CN105785358B (zh) | 2016-05-19 | 2016-05-19 | 一种方向余弦坐标系下带多普勒量测的雷达目标跟踪方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610339346.3A CN105785358B (zh) | 2016-05-19 | 2016-05-19 | 一种方向余弦坐标系下带多普勒量测的雷达目标跟踪方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN105785358A CN105785358A (zh) | 2016-07-20 |
CN105785358B true CN105785358B (zh) | 2017-04-12 |
Family
ID=56380082
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201610339346.3A Active CN105785358B (zh) | 2016-05-19 | 2016-05-19 | 一种方向余弦坐标系下带多普勒量测的雷达目标跟踪方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN105785358B (zh) |
Families Citing this family (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106646453B (zh) * | 2016-11-17 | 2019-04-05 | 电子科技大学 | 一种基于预测值量测转换的多普勒雷达目标跟踪方法 |
CN106950562B (zh) * | 2017-03-30 | 2020-02-18 | 电子科技大学 | 一种基于预测值量测转换的状态融合目标跟踪方法 |
CN108279412A (zh) * | 2018-01-30 | 2018-07-13 | 哈尔滨工业大学 | 一种目的地约束下目标跟踪装置及方法 |
CN111077518B (zh) * | 2019-12-20 | 2020-11-06 | 哈尔滨工业大学 | 一种基于距离-多普勒量测的跟踪滤波方法及装置 |
CN114089288B (zh) * | 2022-01-12 | 2022-04-15 | 中国人民解放军空军预警学院 | 一种相控阵雷达抗干扰方法、装置及存储介质 |
CN117630993B (zh) * | 2024-01-15 | 2024-04-12 | 华中科技大学 | 一种基于sair多快拍的rfi源地理定位方法 |
Family Cites Families (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US5912640A (en) * | 1997-08-26 | 1999-06-15 | Lockheed Martin Corporation | Boost engine cutoff estimation in Doppler measurement system |
CN103048658A (zh) * | 2012-11-10 | 2013-04-17 | 中国人民解放军海军航空工程学院 | 基于径向加速度的RA-Singer-EKF机动目标跟踪算法 |
-
2016
- 2016-05-19 CN CN201610339346.3A patent/CN105785358B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN105785358A (zh) | 2016-07-20 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN105785358B (zh) | 一种方向余弦坐标系下带多普勒量测的雷达目标跟踪方法 | |
CN103278813B (zh) | 一种基于高阶无迹卡尔曼滤波的状态估计方法 | |
CN105954742B (zh) | 一种球坐标系下带多普勒观测的雷达目标跟踪方法 | |
CN106950562B (zh) | 一种基于预测值量测转换的状态融合目标跟踪方法 | |
CN102288978B (zh) | 一种cors基站周跳探测与修复方法 | |
CN106885576B (zh) | 一种基于多点地形匹配定位的auv航迹偏差估计方法 | |
CN103940433B (zh) | 一种基于改进的自适应平方根ukf算法的卫星姿态确定方法 | |
CN104330083A (zh) | 基于平方根无迹卡尔曼滤波的多机器人协同定位算法 | |
CN102795323B (zh) | 一种基于ukf的水下机器人状态和参数联合估计方法 | |
CN101853243A (zh) | 系统模型未知的自适应卡尔曼滤波方法 | |
CN104166989B (zh) | 一种用于二维激光雷达点云匹配的快速icp方法 | |
CN105022040A (zh) | 基于杂波数据联合拟合的阵元误差估计方法 | |
CN103900564A (zh) | 一种潜深辅助地磁异常反演测速/水下连续定位方法 | |
CN105572629A (zh) | 一种适用于任意阵列结构的低运算复杂度的二维测向方法 | |
CN116882134A (zh) | 一种多基准站网的gnss基线联合解算方法及计算机可读介质 | |
CN103616024A (zh) | 一种行星探测进入段自主导航系统可观测度确定方法 | |
CN104614751B (zh) | 基于约束信息的目标定位方法 | |
Vecherin et al. | Three-dimensional acoustic travel-time tomography of the atmosphere | |
CN104021311A (zh) | 一种基于Hermite函数约束的数据融合计算方法 | |
CN104022757B (zh) | 一种高阶矩匹配的多层无迹卡尔曼滤波器的线性扩展方法 | |
CN104330772B (zh) | 基于多向寻优的全跟踪式ukf滤波算法的双站定位方法 | |
Liu et al. | Algorithm for HF radar vector current measurements | |
Komaragiri et al. | A sag monitoring device based on a cluster of code based GPS receivers | |
CN107063195A (zh) | 一种基于递归位置估计的大规模水下网络定位方法 | |
CN102645646B (zh) | 一种多信源不确定融合定位方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |