CN107315171A - 一种雷达组网目标状态与系统误差联合估计算法 - Google Patents
一种雷达组网目标状态与系统误差联合估计算法 Download PDFInfo
- Publication number
- CN107315171A CN107315171A CN201710529734.2A CN201710529734A CN107315171A CN 107315171 A CN107315171 A CN 107315171A CN 201710529734 A CN201710529734 A CN 201710529734A CN 107315171 A CN107315171 A CN 107315171A
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- mtd
- mover
- radar
- 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.)
- Granted
Links
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 43
- 230000009897 systematic effect Effects 0.000 title claims abstract description 20
- 239000002131 composite material Substances 0.000 title abstract 2
- 238000005259 measurement Methods 0.000 claims abstract description 62
- 239000011159 matrix material Substances 0.000 claims abstract description 35
- 238000000034 method Methods 0.000 claims abstract description 32
- 238000001914 filtration Methods 0.000 claims abstract description 19
- 230000008569 process Effects 0.000 claims abstract description 13
- 230000006855 networking Effects 0.000 claims description 29
- 238000013178 mathematical model Methods 0.000 claims description 14
- 230000001902 propagating effect Effects 0.000 claims description 6
- 239000013598 vector Substances 0.000 claims description 5
- 230000001360 synchronised effect Effects 0.000 claims description 4
- 230000003595 spectral effect Effects 0.000 claims description 3
- 230000009466 transformation Effects 0.000 abstract description 2
- 230000003190 augmentative effect Effects 0.000 abstract 2
- 230000000694 effects Effects 0.000 description 15
- 230000004927 fusion Effects 0.000 description 10
- 238000010586 diagram Methods 0.000 description 9
- 238000005516 engineering process Methods 0.000 description 8
- 230000008901 benefit Effects 0.000 description 4
- 238000003672 processing method Methods 0.000 description 4
- 238000004088 simulation Methods 0.000 description 4
- 238000006243 chemical reaction Methods 0.000 description 3
- 230000007547 defect Effects 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 238000007476 Maximum Likelihood Methods 0.000 description 2
- 239000000654 additive Substances 0.000 description 2
- 230000000996 additive effect Effects 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000000342 Monte Carlo simulation Methods 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
- 230000007704 transition 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
- 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
- G01S7/40—Means for monitoring or calibrating
-
- 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/87—Combinations of radar systems, e.g. primary radar and secondary radar
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)
- Medicines Containing Antibodies Or Antigens For Use As Internal Diagnostic Agents (AREA)
Abstract
本发明提供一种雷达组网目标状态与系统误差联合估计算法,步骤如下,根据机载雷达测量原理,结合机载平台的地理位置及姿态角,构建出含有系统误差、测量噪声以及目标状态的每部机载雷达量测数学模型;根据目标运动状态,联合每部机载雷达的系统误差,构建出每一时刻含有过程噪声的扩维后目标运动状态数学模型;根据上两步构建的数学模型,设置测量噪声、过程噪声、扩维后目标运动状态初值及其估计误差协方差矩阵初值,应用CKF滤波方法,实现对目标运动状态及每部机载雷达系统误差的实时同步估计。本发明解决了传统基于坐标投影的二维平面系统误差配准算法中存在的诸如数据变形、无法估计俯仰角系统误差及不适于远距离系统误差配准等问题。
Description
技术领域
本发明涉及雷达组网信息融合目标跟踪领域,特别是涉及量测高精度、系统高维数、复杂强非线性条件下的机载雷达多平台协同跟踪领域,具体涉及一种雷达组网目标状态与系统误差联合估计算法。
背景技术
雷达组网信息融合所能带来的巨大效益已经得到了世界各国的公认。尽管雷达组网信息融合技术依然在不断飞速发展,但是由于实际系统中各雷达探测目标误差的存在,系统实时融合效果的保障已成为雷达组网信息融合技术领域长期以来十分棘手的问题。实际应用表明,在多雷达组网跟踪系统中,雷达系统误差的存在会导致目标跟踪均方根误差比理论值要大。当雷达系统误差太大时,就会出现多部雷达融合跟踪甚至不如单部雷达目标跟踪效果的情况。最恶劣情况下,雷达系统误差会导致来自同一轨迹的多雷达量测互联失败,产生相对同一目标的多条航迹,这样本来是同一个目标的航迹,却由于相互偏差较大而可能被认为是不同的目标,从而给航迹关联及融合带来模糊和困难,使融合得到的系统航迹的性能下降。尤其是在目标密集、编队飞行等复杂的场景中更容易造成航迹关联混乱、融合精度降低,进而使整体系统融合失去意义,丧失了雷达组网系统本身应有的优点。因此,雷达组网系统误差配准技术是确保雷达网性能稳定所不可缺少的基础和关键技术,其目的是准确估计并消除组网雷达的系统误差。
从坐标转换的角度来看,常用的雷达组网系统误差配准算法主要可以分为两类:其中一类是基于球(极)投影的二维平面系统误差配准算法。该类算法是把各雷达的量测投影到与地球正切的局部载体坐标系上,然后再变换到二维平面来估计各雷达的系统误差。此类方法在工程上应用较多,具有算法简单、便于实现等优点,但是此类雷达系统误差配准技术也存在以下缺陷:
①球(极)投影法虽然利用高阶近似来提高精度,但由于地球是椭球而不是圆球,所以在投影时会给测量引入误差;
②球(极)投影法会使数据变形。球(极)保角投影,只保证方位角不变形,并不能保证斜距不变形,这样会导致系统误差不再是常数,而且与量测有关;
③在二维公共坐标系中只能估计方位角和径向距离的系统误差,无法估计俯仰角系统误差。
所以基于球(极)投影法的系统误差配准技术通常用于短距离雷达组网系统误差配准,而对于远距离雷达系统误差的配准,我们多采用基于地心地固(ECEF)坐标系的雷达系统误差配准技术。
从数据处理的角度来看,雷达组网系统误差配准算法又可分为离线处理法和在线处理法。离线处理法主要包括最小二乘算法(LS)、广义最小二乘算法(GLS)、最大似然算法(ML)和精确极大似然估计算法(EML)。该类算法需要对一段时间内的数据进行集中处理,因此计算量相对比较大。在线处理法主要是基于扩展卡尔曼滤波(EKF)和无迹卡尔曼滤波(UKF)的实时估计算法,计算量相对较小。但是EKF需要计算雅克比矩阵,在量测高精度、系统高维数、复杂强非线性的情况下往往不能很好的近似系统状态的后验概率密度函数,导致滤波精度大大降低,甚至造成滤波发散。虽然UKF避免了EKF的缺点,估计精度有一定的提高,但当UKF应用于高维非线性系统时,其确定性采样点的权值容易出现负值的情况,导致其矩积分中引入截断误差,无法保证数值稳定性和状态协方差矩阵的半正定性,造成滤波精度下降。
发明内容
本发明的目的在于提供一种雷达组网目标状态与系统误差联合估计算法,克服或减轻现有技术的至少一个上述缺陷。
本发明的目的通过如下技术方案实现:一种雷达组网目标状态与系统误差联合估计算法,包括如下步骤,
步骤一:根据机载雷达测量原理,结合机载平台的地理位置及姿态角,构建出含有系统误差、测量噪声以及目标状态的每部机载雷达量测数学模型;
步骤二:根据目标运动状态,联合每部机载雷达的系统误差,构建出每一时刻含有过程噪声的扩维后目标运动状态数学模型;
步骤三:根据步骤一、步骤二构建的数学模型,设置测量噪声、过程噪声、扩维后目标运动状态初值及其估计误差协方差矩阵初值,应用CKF滤波方法,实现对目标运动状态及每部机载雷达系统误差的实时同步估计。
优选地是,所述机载雷达为两部。
优选地是,所述机载雷达测量原理通过如下公式实施:
x、y、z表示为目标相对于每部机载雷达的位置在ECEF坐标系下每个坐标轴上的投影。
优选地是,所述机载平台的地理位置为机载平台所在大地坐标系下的经度、纬度及海拔高度;所述机载平台的姿态角为偏航角、纵摇角及横摇角。
优选地是,所述步骤一所构建的含有系统误差、测量噪声以及目标状态的每部机载雷达量测数学模型公式:
ri(k)为目标距离,θi(k)为目标方位角,ηi(k)为目标俯仰角;其中φi(k),αi(k)分别表示k时刻机载雷达i所在载体平台的偏航角、纵摇角及横摇角;A为由机载直角坐标系转换到END坐标系的坐标旋转矩阵表示为:
k时刻机载雷达i在大地坐标系的坐标为Xisp(k)=[Li(k) λi(k) Hi(k)]T,Li(k),λi(k),Hi(k)分别表示k时刻机载雷达i所在的纬度、经度及海拔高度;T为由END坐标系转换到ECEF坐标系的坐标旋转矩阵矩阵表示为:
Xt(k)=[xt(k) yt(k) zt(k)]T表示为k时刻目标在ECEF坐标系下的坐标;Xis(k)表示为k时刻机载雷达i在ECEF坐标系下的位置表示为:
其中,Eq表示赤道半径,e表示地球偏心率;bi(k)为系统误差,ni(k)为测量噪声。
优选地是,所述步骤二所构建的每一时刻含有过程噪声的扩维后目标运动状态数学模型公式:
XA(k+1)=FA(k)XA(k)+WA(k)
其中,FA(k)=diag(F(k),I6×6),F(k)=diag(Fx(k),Fy(k),Fz(k)),i=x,y,z;WA(k)为零均值高斯白噪声;XA(k)为联合目标状态X(k)和两部雷达的系统误差b1(k),b2(k)扩维后的目标运动状态,表示为XA(k)=[X(k)T b1(k)T b2(k)T]T。
优选地是,所述测量噪声表示为 分别表示目标距离、方位角和俯仰角测量误差的协方差;过程噪声表示为QA=diag(Q,06×6),Q=diag(Qx,Qy,Qz),i=x,y,z,q表示噪声的功率谱密度。
优选地是,所述CKF滤波方法包括时间更新与量测更新两个步骤,
A、时间更新:
Ⅰ、通过Cholesky分解系统状态估计误差协方差矩阵P(k|k);
P(k|k)=S(k|k)ST(k|k)
Ⅱ、根据Spherical-Radial Cubature准则计算Cubature点;
Ⅲ、通过状态方程传播Cubature点;
Ⅳ、计算系统状态的先验均值及先验协方差矩阵;
Ⅴ、计算系统状态的先验协方差矩阵;
B、量测更新:
Ⅰ、通过Cholesky分解系统状态先验协方差矩阵P(k+1|k);
P(k+1|k)=S(k+1|k)ST(k+1|k)
Ⅱ、根据Spherical-Radial Cubature准则计算Cubature点;
Ⅲ、通过量测方程传播Cubature点;
Zi(k+1|k)=hA(Xi(k+1|k))
Ⅳ、计算量测的先验均值;
Ⅴ、计算量测的先验协方差矩阵;
Ⅵ、计算量测和状态向量的互相关协方差矩阵;
Ⅶ、计算卡尔曼滤波增益;
Ⅷ、计算系统状态的后验均值;
Ⅸ、计算系统状态的后验协方差矩阵;
P(k+1|k+1)=P(k+1|k)-W(k+1)Pzz(k+1|k)WT(k+1)。
本发明所提供的一种雷达组网目标状态与系统误差联合估计算法的有益效果在于,首先采用基于ECEF坐标系的雷达系统误差配准技术,解决了传统的基于坐标投影的二维平面系统误差配准算法中存在的诸如数据变形、无法估计俯仰角系统误差及不适于远距离系统误差配准等固有问题。其次,该方法对系统状态进行扩维处理,把系统误差作为系统的未知、待估计状态,并引入容积卡尔曼滤波(CKF)算法,实现了对目标状态和系统误差的联合估计。本发明所提出的方法无需计算雅克比矩阵,而且在递推运算中具有很强的数值稳定性,解决了传统非线性高斯滤波器在量测高精度、系统高维数、复杂强非线性情况下对系统状态后验概率密度函数近似精度不高这一问题,提高了雷达组网条件下对系统误差的估计精度及目标的跟踪精度,有效地实现了雷达组网系统误差的实时配准,同时提升了目标跟踪的可靠性和稳定性。同时本发明具有良好的扩展性和适应性,可广泛应用于火控、监视、预警等战术功能,对稳定目标跟踪有较高要求的多平台有源/无源雷达信息融合跟踪系统,应用前景广阔,应用价值巨大。
附图说明
图1是本发明雷达组网目标状态与系统误差联合估计算法的流程图;
图2是根据本发明一个实施例中的仿真环境经纬图;
图3是根据本发明一个实施例中的第一雷达距离系统误差估计效果图;
图4是根据本发明一个实施例中的第二雷达距离系统误差估计效果图;
图5是根据本发明一个实施例中的第一雷达方位角系统误差估计效果图;
图6是根据本发明一个实施例中的第二雷达方位角系统误差估计效果图;
图7是根据本发明一个实施例中的第一雷达俯仰角系统误差估计效果图;
图8是根据本发明一个实施例中的第二雷达俯仰角系统误差估计效果图;
图9是根据本发明一个实施例中的第一雷达径向距离估计效果图;
图10是根据本发明一个实施例中的第二雷达径向距离估计效果图;
图11是根据本发明一个实施例中的第一雷达方位角估计效果图;
图12是根据本发明一个实施例中的第二雷达方位角估计效果图;
图13是根据本发明一个实施例中的第一雷达俯仰角估计效果图;
图14是根据本发明一个实施例中的第二雷达俯仰角估计效果图。
具体实施方式
为使本发明实施的目的、技术方案和优点更加清楚,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行更加详细的描述。在附图中,自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。所描述的实施例是本发明一部分实施例,而不是全部的实施例。下面通过参考附图描述的实施例是示例性的,旨在用于解释本发明,而不能理解为对本发明的限制。基于本发明中的实施例,本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
下面结合附图对本发明的雷达组网目标状态与系统误差联合估计算法做进一步详细说明。
一种雷达组网目标状态与系统误差联合估计算法,其核心在于机载雷达系统的数学建模及非线性滤波算法的设计。首先完成对机载雷达系统数学模型的构建;然后采用基于ECEF坐标系的系统误差配准技术,以解决传统基于坐标投影的二维平面系统误差配准算法中存在的诸如数据变形、无法估计俯仰角系统误差及不适于远距离系统误差配准等固有问题。最后对系统状态进行扩维处理,把系统误差作为系统未知、待估计的状态,并引入CKF滤波算法,实现对目标状态和系统误差的联合估计。见图1所示,具体包括以下几个步骤:
1)步骤一:根据机载雷达测量原理,结合机载平台的地理位置及姿态角,构建出含有系统误差、测量噪声以及目标状态的每部机载雷达量测数学模型。
本发明中的系统选择由两部三坐标机载雷达组成。机载雷达i(i=1,2)在极坐标系下对同一目标进行同步测量,存在距离系统误差方位角系统误差俯仰角系统误差且假设它们均为常量加性系统误差,表示为
机载雷达i在k时刻测量得到含有系统误差bi(k)和测量噪声ni(k)的目标距离ri(k)、目标方位角θi(k)和目标俯仰角ηi(k)。其中ni(k)为零均值高斯白噪声,对应的噪声协方差矩阵表示为 分别表示目标距离、目标方位角和目标俯仰角测量误差的协方差,且系统误差与测量噪声之间是相互独立的。
假设不含系统误差的真值量测表示为ri'(k),θi'(k),ηi'(k),则k时刻机载雷达i在极坐标系下的测量方程为:
将极坐标系下的测量转化到局部直角载体坐标系下,可得:
式中h-1(r,θ,η)=[r sinθcosη,r cosθcosη,r sinη]T。
假设k时刻机载雷达i所在的载体平台姿态角表示为其中φi(k),αi(k)分别表示k时刻机载雷达i所在载体平台的偏航角、纵摇角及横摇角。则根据载体坐标系和END坐标系间的转换关系,可把机载雷达i在载体坐标系的测量转换到END坐标系中:
式中,
假设k时刻机载雷达i在大地坐标系的坐标为Xisp(k)=[Li(k) λi(k) Hi(k)]T,Li(k),λi(k),Hi(k)分别表示k时刻机载雷达i所在的纬度、经度及海拔高度。则根据大地坐标系和ECEF坐标系间的转换关系,可把k时刻机载雷达i在大地坐标系的位置转换到ECEF坐标系中:
式中,Eq表示赤道半径,e表示地球偏心率。根据END坐标系和ECEF坐标系间的转换关系,将k时刻机载雷达i在END坐标系的量测转换到ECEF坐标系中:
式中
因此,k时刻机载雷达i在ECEF坐标系下的量测可以表示为:
假设k时刻目标在ECEF坐标系下的坐标为Xt(k)=[xt(k) yt(k) zt(k)]T,如果机载雷达i不存在系统误差,且不含测量噪声,那么此时机载雷达i在ECEF坐标系下的量测就应该等同于目标此时在ECEF坐标系下的坐标,即:
Xis(k)+T(Xisp(k))A(vi(k))h-1(ri'(k),θi'(k),ηi'(k))=Xt(k)
经过整理,可得机载雷达的量测方程为:
其中,x、y、z表示为目标相对于每部机载雷达的位置在ECEF坐标系下每个坐标轴上的投影。
2)步骤二:根据目标运动状态,联合每部机载雷达的系统误差,构建出每一时刻含有过程噪声的扩维后目标运动状态数学模型。
首先将系统的状态进行扩维,把系统误差作为系统未知、待估计的状态。联合目标状态X(k)和两部雷达的系统误差b1(k),b2(k),构建新的系统状态(即扩维后目标运动状态)为:
XA(k)=[X(k)T b1(k)T b2(k)T]T
由于假设机载雷达的系统误差为不变常量,因此可得:
bi(k+1)=I3×3bi(k)
扩维后,新的系统状态XA(k)的状态转移方程可表示为:
XA(k+1)=FA(k)XA(k)+WA(k)
其中,FA(k)=diag(F(k),I6×6),F(k)=diag(Fx(k),Fy(k),Fz(k)),i=x,y,z。WA(k)=diag(W(k),06×6),WA(k)为零均值高斯白噪声,其噪声协方差矩阵为QA=diag(Q,06×6),Q=diag(Qx,Qy,Qz),i=x,y,z,q表示噪声的功率谱密度。
3)根据步骤一、步骤二构建的数学模型,设置测量噪声、过程噪声、扩维后目标运动状态初值及其估计误差协方差矩阵初值,应用CKF滤波方法,实现对目标运动状态及每部机载雷达系统误差的实时同步估计。
首先将步骤二中构建的新的系统状态和机载雷达1、2的量测方程联合,可得:
ZA(k)=hA(XA(k))+nA(k)
式中, nA(k)仍为零均值高斯噪声,对应的噪声协方差矩阵表示为RA(k)=diag(R1(k),R2(k))。
接下来对扩维后的系统状态方程和量测方程应用CKF滤波方法,进而对系统状态及系统误差进行联合估计。其中,CKF滤波方法包括时间更新与量测更新两个步骤。
A、时间更新:
Ⅰ、通过Cholesky分解系统状态估计误差协方差矩阵P(k|k);
P(k|k)=S(k|k)ST(k|k)
Ⅱ、根据Spherical-Radial Cubature准则计算Cubature点;
Ⅲ、通过状态方程传播Cubature点;
Ⅳ、计算系统状态的先验均值及先验协方差矩阵;
Ⅴ、计算系统状态的先验协方差矩阵;
B、量测更新:
Ⅰ、通过Cholesky分解系统状态先验协方差矩阵P(k+1|k);
P(k+1|k)=S(k+1|k)ST(k+1|k)
Ⅱ、根据Spherical-Radial Cubature准则计算Cubature点;
Ⅲ、通过量测方程传播Cubature点;
Zi(k+1|k)=hA(Xi(k+1|k))
Ⅳ、计算量测的先验均值;
Ⅴ、计算量测的先验协方差矩阵;
Ⅵ、计算量测和状态向量的互相关协方差矩阵;
Ⅶ、计算卡尔曼滤波增益;
Ⅷ、计算系统状态的后验均值;
Ⅸ、计算系统状态的后验协方差矩阵;
P(k+1|k+1)=P(k+1|k)-W(k+1)Pzz(k+1|k)WT(k+1)。
需要说明的是,假设在加性噪声条件下,考虑如下离散形式的系统状态空间模型:
其中,表示k时刻系统状态向量(nx为状态维数),表示k时刻外部量测向量(nz为量测维数)。FA(·)表示状态转移函数,hA(·)表示量测函数。与分别表示系统噪声和量测噪声,二者互不相关且均为零均值高斯白噪声,噪声协方差矩阵分别为QA(k)和RA(k)。
下面结合一个实施例通过其数值仿真详细说明该发明所提出的方法。
假设系统由两部机载雷达1、2和一个空中目标构成。机载雷达平台设置为直升机平台,飞行速度较慢,飞行海拔高度为1km,目标设置为战斗机目标,飞行速度较快,飞行海拔高度为2km。机载雷达扫描周期为1s,目标沿经线飞行,具体仿真场景如图2所示。
各机载雷达的量测噪声ni(k)均为零均值高斯白噪声,目标距离、目标方位角和目标俯仰角测量误差的协方差分别为 雷达1机载平台姿态角变化规律可以描述为v1(k)=[0.002k,0.01+0.002k,0.01+0.002k]T,雷达2机载平台姿态角变化规律可以描述为v2(k)=[0.002k,0.001k,0.001k]T,各个雷达的系统误差都为bi=[1000m,0.0087rad,0.0087rad]T。下面对ECEF-CKF-ASR算法的估计性能进行仿真分析。对ECEF-CKF-ASR算法进行50次蒙特卡洛仿真结果如图3~图14所示。可以看出,随着ECEF-CKF-ASR算法估计的收敛,ECEF-CKF-ASR算法给出的目标状态估计基本消除了系统误差的影响。
表1 ECEF-CKF-ASR算法系统误差估计精度
从表1可以看出,ECEF-CKF-ASR算法对各个系统误差的估计精度基本上都达到95%以上,可实现算法收敛。因此ECEF-CKF-ASR算法对各个雷达的系统误差具有良好的估计效果,很好地解决了机动雷达的误差配准问题。
以上所述,仅为本发明的具体实施方式,但本发明的保护范围并不局限于此,任何熟悉本技术领域的技术人员在本发明揭露的技术范围内,可轻易想到的变化或替换,都应涵盖在本发明的保护范围之内。因此,本发明的保护范围应以所述权利要求的保护范围为准。
Claims (8)
1.一种雷达组网目标状态与系统误差联合估计算法,其特征在于,包括如下步骤,
步骤一:根据机载雷达测量原理,结合机载平台的地理位置及姿态角,构建出含有系统误差、测量噪声以及目标状态的每部机载雷达量测数学模型;
步骤二:根据目标运动状态,联合每部机载雷达的系统误差,构建出每一时刻含有过程噪声的扩维后目标运动状态数学模型;
步骤三:根据步骤一、步骤二构建的数学模型,设置测量噪声、过程噪声、扩维后目标运动状态初值及其估计误差协方差矩阵初值,应用CKF滤波方法,实现对目标运动状态及每部机载雷达系统误差的实时同步估计。
2.根据权利要求1所述的雷达组网目标状态与系统误差联合估计算法,其特征在于,所述机载雷达为两部。
3.根据权利要求2所述的雷达组网目标状态与系统误差联合估计算法,其特征在于,所述机载雷达测量原理通过如下公式实施:
<mrow>
<mi>h</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msqrt>
<mrow>
<msup>
<mi>x</mi>
<mn>2</mn>
</msup>
<mo>+</mo>
<msup>
<mi>y</mi>
<mn>2</mn>
</msup>
<mo>+</mo>
<msup>
<mi>z</mi>
<mn>2</mn>
</msup>
</mrow>
</msqrt>
<mo>,</mo>
<mi>arctan</mi>
<mrow>
<mo>(</mo>
<mfrac>
<mi>y</mi>
<mi>x</mi>
</mfrac>
<mo>)</mo>
</mrow>
<mo>,</mo>
<mi>arctan</mi>
<mrow>
<mo>(</mo>
<mfrac>
<mi>z</mi>
<msqrt>
<mrow>
<msup>
<mi>x</mi>
<mn>2</mn>
</msup>
<mo>+</mo>
<msup>
<mi>y</mi>
<mn>2</mn>
</msup>
</mrow>
</msqrt>
</mfrac>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
</mrow>
<mi>T</mi>
</msup>
</mrow>
x、y、z表示为目标相对于每部机载雷达的位置在ECEF坐标系下每个坐标轴上的投影。
4.根据权利要求2所述的雷达组网目标状态与系统误差联合估计算法,其特征在于,所述机载平台的地理位置为机载平台所在大地坐标系下的经度、纬度及海拔高度;所述机载平台的姿态角为偏航角、纵摇角及横摇角。
5.根据权利要求2所述的雷达组网目标状态与系统误差联合估计算法,其特征在于,所述步骤一所构建的含有系统误差、测量噪声以及目标状态的每部机载雷达量测数学模型公式:
<mrow>
<msub>
<mi>Z</mi>
<mrow>
<mi>i</mi>
<mi>p</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>r</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>&theta;</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>&eta;</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>=</mo>
<mi>h</mi>
<mo>(</mo>
<mo>&lsqb;</mo>
<msup>
<mi>A</mi>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mrow>
<mo>(</mo>
<msub>
<mi>v</mi>
<mi>i</mi>
</msub>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mo>&lsqb;</mo>
<msup>
<mi>T</mi>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mrow>
<mo>(</mo>
<msub>
<mi>X</mi>
<mrow>
<mi>i</mi>
<mi>s</mi>
<mi>p</mi>
</mrow>
</msub>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mo>&lsqb;</mo>
<msub>
<mi>X</mi>
<mi>t</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msub>
<mi>X</mi>
<mrow>
<mi>i</mi>
<mi>s</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mo>)</mo>
<mo>+</mo>
<msub>
<mi>b</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mi>n</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
ri(k)为目标距离,θi(k)为目标方位角,ηi(k)为目标俯仰角;其中φi(k),αi(k)分别表示k时刻机载雷达i所在载体平台的偏航角、纵摇角及横摇角;A为由机载直角坐标系转换到END坐标系的坐标旋转矩阵表示为:
k时刻机载雷达i在大地坐标系的坐标为Xisp(k)=[Li(k) λi(k) Hi(k)]T,Li(k),λi(k),Hi(k)分别表示k时刻机载雷达i所在的纬度、经度及海拔高度;T为由END坐标系转换到ECEF坐标系的坐标旋转矩阵矩阵表示为:
<mrow>
<mi>T</mi>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<mo>-</mo>
<mi>s</mi>
<mi>i</mi>
<mi>n</mi>
<mi>&lambda;</mi>
</mrow>
</mtd>
<mtd>
<mrow>
<mo>-</mo>
<mi>sin</mi>
<mi> </mi>
<mi>L</mi>
<mi> </mi>
<mi>c</mi>
<mi>o</mi>
<mi>s</mi>
<mi>&lambda;</mi>
</mrow>
</mtd>
<mtd>
<mrow>
<mi>cos</mi>
<mi> </mi>
<mi>L</mi>
<mi> </mi>
<mi>c</mi>
<mi>o</mi>
<mi>s</mi>
<mi>&lambda;</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>-</mo>
<mi>cos</mi>
<mi>&lambda;</mi>
</mrow>
</mtd>
<mtd>
<mrow>
<mo>-</mo>
<mi>sin</mi>
<mi> </mi>
<mi>L</mi>
<mi> </mi>
<mi>s</mi>
<mi>i</mi>
<mi>n</mi>
<mi>&lambda;</mi>
</mrow>
</mtd>
<mtd>
<mrow>
<mi>cos</mi>
<mi> </mi>
<mi>L</mi>
<mi> </mi>
<mi>s</mi>
<mi>i</mi>
<mi>n</mi>
<mi>&lambda;</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mrow>
<mi>cos</mi>
<mi> </mi>
<mi>L</mi>
</mrow>
</mtd>
<mtd>
<mrow>
<mi>sin</mi>
<mi> </mi>
<mi>L</mi>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
1
Xt(k)=[xt(k) yt(k) zt(k)]T表示为k时刻目标在ECEF坐标系下的坐标;
Xis(k)表示为k时刻机载雷达i在ECEF坐标系下的位置表示为:
<mrow>
<msub>
<mi>X</mi>
<mrow>
<mi>i</mi>
<mi>s</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mo>&lsqb;</mo>
<msub>
<mi>C</mi>
<mi>i</mi>
</msub>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
<mo>+</mo>
<msub>
<mi>H</mi>
<mi>i</mi>
</msub>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
<mo>&rsqb;</mo>
<mi>c</mi>
<mi>o</mi>
<mi>s</mi>
<mo>(</mo>
<msub>
<mi>L</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>)</mo>
<mi>c</mi>
<mi>o</mi>
<mi>s</mi>
<mo>(</mo>
<msub>
<mi>&lambda;</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>)</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>&lsqb;</mo>
<msub>
<mi>C</mi>
<mi>i</mi>
</msub>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
<mo>+</mo>
<msub>
<mi>H</mi>
<mi>i</mi>
</msub>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
<mo>&rsqb;</mo>
<mi>c</mi>
<mi>o</mi>
<mi>s</mi>
<mo>(</mo>
<msub>
<mi>L</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>)</mo>
<mi>s</mi>
<mi>i</mi>
<mi>n</mi>
<mo>(</mo>
<msub>
<mi>&lambda;</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>)</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>&lsqb;</mo>
<msub>
<mi>C</mi>
<mi>i</mi>
</msub>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
<mo>(</mo>
<mn>1</mn>
<mo>-</mo>
<msup>
<mi>e</mi>
<mn>2</mn>
</msup>
<mo>)</mo>
<mo>+</mo>
<msub>
<mi>H</mi>
<mi>i</mi>
</msub>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
<mo>&rsqb;</mo>
<mi>s</mi>
<mi>i</mi>
<mi>n</mi>
<mo>(</mo>
<msub>
<mi>&lambda;</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>)</mo>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
其中,Eq表示赤道半径,e表示地球偏心率;bi(k)为系统误差,ni(k)为测量噪声。
6.根据权利要求2所述的雷达组网目标状态与系统误差联合估计算法,其特征在于,所述步骤二所构建的每一时刻含有过程噪声的扩维后目标运动状态数学模型公式:
XA(k+1)=FA(k)XA(k)+WA(k)
其中,FA(k)=diag(F(k),I6×6),F(k)=diag(Fx(k),Fy(k),Fz(k)),i=x,y,z;WA(k)为零均值高斯白噪声;XA(k)为联合目标状态X(k)和两部雷达的系统误差b1(k),b2(k)扩维后的目标运动状态,表示为XA(k)=[X(k)T b1(k)T b2(k)T]T。
7.根据权利要求2所述的雷达组网目标状态与系统误差联合估计算法,其特征在于,所述测量噪声表示为 分别表示目标距离、方位角和俯仰角测量误差的协方差;过程噪声表示为QA=diag(Q,06×6),Q=diag(Qx,Qy,Qz),i=x,y,z,q表示噪声的功率谱密度。
8.根据权利要求2所述的雷达组网目标状态与系统误差联合估计算法,其特征在于,所述CKF滤波方法包括时间更新与量测更新两个步骤,
A、时间更新:
Ⅰ、通过Cholesky分解系统状态估计误差协方差矩阵P(k|k);
P(k|k)=S(k|k)ST(k|k)
Ⅱ、根据Spherical-Radial Cubature准则计算Cubature点;
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>&xi;</mi>
<mi>i</mi>
</msub>
<mo>=</mo>
<msqrt>
<mfrac>
<mi>m</mi>
<mn>2</mn>
</mfrac>
</msqrt>
<msub>
<mrow>
<mo>&lsqb;</mo>
<mn>1</mn>
<mo>&rsqb;</mo>
</mrow>
<mi>i</mi>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>w</mi>
<mi>i</mi>
</msub>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mi>m</mi>
</mfrac>
<mo>,</mo>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>...</mn>
<mi>m</mi>
<mo>,</mo>
<mi>m</mi>
<mo>=</mo>
<mn>2</mn>
<msub>
<mi>n</mi>
<mi>x</mi>
</msub>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mrow>
<msub>
<mi>X</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mi>S</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<msub>
<mi>&xi;</mi>
<mi>i</mi>
</msub>
<mo>+</mo>
<mover>
<mi>X</mi>
<mo>^</mo>
</mover>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
Ⅲ、通过状态方程传播Cubature点;
<mrow>
<msubsup>
<mi>X</mi>
<mi>i</mi>
<mo>*</mo>
</msubsup>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msub>
<mi>F</mi>
<mi>A</mi>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>X</mi>
<mi>i</mi>
</msub>
<mo>(</mo>
<mrow>
<mi>k</mi>
<mo>|</mo>
<mi>k</mi>
</mrow>
<mo>)</mo>
<mo>)</mo>
</mrow>
</mrow>
Ⅳ、计算系统状态的先验均值及先验协方差矩阵;
<mrow>
<mover>
<mi>X</mi>
<mo>^</mo>
</mover>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mi>m</mi>
</mfrac>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msubsup>
<mi>X</mi>
<mi>i</mi>
<mo>*</mo>
</msubsup>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
Ⅴ、计算系统状态的先验协方差矩阵;
<mrow>
<mi>P</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mi>m</mi>
</mfrac>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msubsup>
<mi>X</mi>
<mi>i</mi>
<mo>*</mo>
</msubsup>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<msubsup>
<mi>X</mi>
<mi>i</mi>
<mrow>
<mo>*</mo>
<mi>T</mi>
</mrow>
</msubsup>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mover>
<mi>X</mi>
<mo>^</mo>
</mover>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<msup>
<mover>
<mi>X</mi>
<mo>^</mo>
</mover>
<mi>T</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>Q</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
B、量测更新:
Ⅰ、通过Cholesky分解系统状态先验协方差矩阵P(k+1|k);
P(k+1|k)=S(k+1|k)ST(k+1|k)
Ⅱ、根据Spherical-Radial Cubature准则计算Cubature点;
<mrow>
<msub>
<mi>X</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mi>S</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<msub>
<mi>&xi;</mi>
<mi>i</mi>
</msub>
<mo>+</mo>
<mover>
<mi>X</mi>
<mo>^</mo>
</mover>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
Ⅲ、通过量测方程传播Cubature点;
Zi(k+1|k)=hA(Xi(k+1|k))
Ⅳ、计算量测的先验均值;
<mrow>
<mover>
<mi>Z</mi>
<mo>^</mo>
</mover>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mi>m</mi>
</mfrac>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msub>
<mi>Z</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
Ⅴ、计算量测的先验协方差矩阵;
<mrow>
<msub>
<mi>P</mi>
<mrow>
<mi>z</mi>
<mi>z</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mi>m</mi>
</mfrac>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msub>
<mi>Z</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<msubsup>
<mi>Z</mi>
<mi>i</mi>
<mi>T</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mover>
<mi>Z</mi>
<mo>^</mo>
</mover>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<msup>
<mover>
<mi>Z</mi>
<mo>^</mo>
</mover>
<mi>T</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>R</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
Ⅵ、计算量测和状态向量的互相关协方差矩阵;
<mrow>
<msub>
<mi>P</mi>
<mrow>
<mi>x</mi>
<mi>z</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mi>m</mi>
</mfrac>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msub>
<mi>X</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<msubsup>
<mi>Z</mi>
<mi>i</mi>
<mi>T</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mover>
<mi>X</mi>
<mo>^</mo>
</mover>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<msup>
<mover>
<mi>Z</mi>
<mo>^</mo>
</mover>
<mi>T</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
Ⅶ、计算卡尔曼滤波增益;
<mrow>
<mi>W</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msub>
<mi>P</mi>
<mrow>
<mi>x</mi>
<mi>z</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<msubsup>
<mi>P</mi>
<mrow>
<mi>z</mi>
<mi>z</mi>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msubsup>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
Ⅷ、计算系统状态的后验均值;
<mrow>
<mover>
<mi>X</mi>
<mo>^</mo>
</mover>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mover>
<mi>X</mi>
<mo>^</mo>
</mover>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>W</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<msub>
<mi>Z</mi>
<mi>A</mi>
</msub>
<mo>(</mo>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
<mo>)</mo>
<mo>-</mo>
<mover>
<mi>Z</mi>
<mo>^</mo>
</mover>
<mo>(</mo>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
</mrow>
<mo>)</mo>
<mo>)</mo>
</mrow>
</mrow>
Ⅸ、计算系统状态的后验协方差矩阵;
P(k+1|k+1)=P(k+1|k)-W(k+1)Pzz(k+1|k)WT(k+1)。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710529734.2A CN107315171B (zh) | 2017-07-02 | 2017-07-02 | 一种雷达组网目标状态与系统误差联合估计算法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710529734.2A CN107315171B (zh) | 2017-07-02 | 2017-07-02 | 一种雷达组网目标状态与系统误差联合估计算法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN107315171A true CN107315171A (zh) | 2017-11-03 |
CN107315171B CN107315171B (zh) | 2021-04-20 |
Family
ID=60179983
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201710529734.2A Active CN107315171B (zh) | 2017-07-02 | 2017-07-02 | 一种雷达组网目标状态与系统误差联合估计算法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN107315171B (zh) |
Cited By (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108872973A (zh) * | 2018-08-30 | 2018-11-23 | 北京电子工程总体研究所 | 一种弹道导弹目标定轨的ekf滤波方法 |
CN109001699A (zh) * | 2018-01-30 | 2018-12-14 | 哈尔滨工业大学 | 基于带噪声目的地信息约束的跟踪方法 |
CN110031802A (zh) * | 2019-04-04 | 2019-07-19 | 中国科学院数学与系统科学研究院 | 具有未知量测零偏的双红外传感器的融合定位方法 |
CN110426702A (zh) * | 2019-07-24 | 2019-11-08 | 中国人民解放军海军航空大学 | 基于多雷达的系统偏差与目标状态的序贯滤波方法及系统 |
CN110426689A (zh) * | 2019-07-02 | 2019-11-08 | 中国航空工业集团公司雷华电子技术研究所 | 一种基于em-cks的机载多平台多传感器系统误差配准算法 |
CN111027204A (zh) * | 2019-12-05 | 2020-04-17 | 中国人民解放军63620部队 | 航天发射光、雷、遥与导航卫星测量数据融合处理方法 |
CN111624594A (zh) * | 2020-05-12 | 2020-09-04 | 中国电子科技集团公司第三十八研究所 | 一种基于转换量测重构的组网雷达跟踪方法 |
CN112285701A (zh) * | 2020-10-22 | 2021-01-29 | 香港中文大学(深圳) | 三维组网雷达系统误差校正方法 |
CN113484832A (zh) * | 2021-07-29 | 2021-10-08 | 西安电子科技大学 | 一种地基雷达组网的系统误差配准方法 |
CN114461968A (zh) * | 2022-01-21 | 2022-05-10 | 中国船舶重工集团公司第七0九研究所 | 一种基于m估计的雷达系统误差鲁棒配准方法及系统 |
CN115267760A (zh) * | 2022-06-25 | 2022-11-01 | 中国人民解放军战略支援部队信息工程大学 | 一种地心地固坐标系下协同被动测向与主动雷达的运动目标定位方法 |
CN118243102A (zh) * | 2024-03-25 | 2024-06-25 | 哈尔滨工业大学 | 一种基于有限间断信息的目标位置估计方法 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103616036A (zh) * | 2013-11-29 | 2014-03-05 | 中国航空无线电电子研究所 | 一种基于合作目标的机载传感器系统误差估计与补偿方法 |
KR20160070383A (ko) * | 2014-12-10 | 2016-06-20 | (주)레코디아 | 위상배열 안테나 기반의 감시망 시스템 |
CN106597428A (zh) * | 2016-12-20 | 2017-04-26 | 中国航空工业集团公司雷华电子技术研究所 | 一种海面目标航向航速估算方法 |
-
2017
- 2017-07-02 CN CN201710529734.2A patent/CN107315171B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103616036A (zh) * | 2013-11-29 | 2014-03-05 | 中国航空无线电电子研究所 | 一种基于合作目标的机载传感器系统误差估计与补偿方法 |
KR20160070383A (ko) * | 2014-12-10 | 2016-06-20 | (주)레코디아 | 위상배열 안테나 기반의 감시망 시스템 |
CN106597428A (zh) * | 2016-12-20 | 2017-04-26 | 中国航空工业集团公司雷华电子技术研究所 | 一种海面目标航向航速估算方法 |
Non-Patent Citations (3)
Title |
---|
KUMAR PAKKI BHARANI CHANDRA ET AL.: "Square Root Cubature Information Filter", 《IEEE SENSORS JOURNAL》 * |
刘瑜等: "基于平方根容积卡尔曼滤波的目标状态与传感器偏差扩维联合估计算法", 《吉林大学学报(工学版)》 * |
胡振涛等: "基于CKF的系统误差与目标状态联合估计算法", 《光电子·激光》 * |
Cited By (19)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109001699A (zh) * | 2018-01-30 | 2018-12-14 | 哈尔滨工业大学 | 基于带噪声目的地信息约束的跟踪方法 |
CN108872973B (zh) * | 2018-08-30 | 2022-07-29 | 北京电子工程总体研究所 | 一种弹道导弹目标定轨的ekf滤波方法 |
CN108872973A (zh) * | 2018-08-30 | 2018-11-23 | 北京电子工程总体研究所 | 一种弹道导弹目标定轨的ekf滤波方法 |
CN110031802B (zh) * | 2019-04-04 | 2020-10-09 | 中国科学院数学与系统科学研究院 | 具有未知量测零偏的双红外传感器的融合定位方法 |
CN110031802A (zh) * | 2019-04-04 | 2019-07-19 | 中国科学院数学与系统科学研究院 | 具有未知量测零偏的双红外传感器的融合定位方法 |
CN110426689A (zh) * | 2019-07-02 | 2019-11-08 | 中国航空工业集团公司雷华电子技术研究所 | 一种基于em-cks的机载多平台多传感器系统误差配准算法 |
CN110426702A (zh) * | 2019-07-24 | 2019-11-08 | 中国人民解放军海军航空大学 | 基于多雷达的系统偏差与目标状态的序贯滤波方法及系统 |
CN110426702B (zh) * | 2019-07-24 | 2021-06-08 | 中国人民解放军海军航空大学 | 基于多雷达的系统偏差与目标状态的序贯滤波方法及系统 |
CN111027204A (zh) * | 2019-12-05 | 2020-04-17 | 中国人民解放军63620部队 | 航天发射光、雷、遥与导航卫星测量数据融合处理方法 |
CN111027204B (zh) * | 2019-12-05 | 2023-07-28 | 中国人民解放军63620部队 | 航天发射光、雷、遥与导航卫星测量数据融合处理方法 |
CN111624594A (zh) * | 2020-05-12 | 2020-09-04 | 中国电子科技集团公司第三十八研究所 | 一种基于转换量测重构的组网雷达跟踪方法 |
CN111624594B (zh) * | 2020-05-12 | 2022-09-23 | 中国电子科技集团公司第三十八研究所 | 一种基于转换量测重构的组网雷达跟踪方法及系统 |
CN112285701A (zh) * | 2020-10-22 | 2021-01-29 | 香港中文大学(深圳) | 三维组网雷达系统误差校正方法 |
CN112285701B (zh) * | 2020-10-22 | 2024-05-10 | 香港中文大学(深圳) | 三维组网雷达系统误差校正方法 |
CN113484832A (zh) * | 2021-07-29 | 2021-10-08 | 西安电子科技大学 | 一种地基雷达组网的系统误差配准方法 |
CN114461968A (zh) * | 2022-01-21 | 2022-05-10 | 中国船舶重工集团公司第七0九研究所 | 一种基于m估计的雷达系统误差鲁棒配准方法及系统 |
CN115267760A (zh) * | 2022-06-25 | 2022-11-01 | 中国人民解放军战略支援部队信息工程大学 | 一种地心地固坐标系下协同被动测向与主动雷达的运动目标定位方法 |
CN115267760B (zh) * | 2022-06-25 | 2023-08-15 | 中国人民解放军战略支援部队信息工程大学 | 一种地心地固坐标系下协同被动测向与主动雷达的运动目标定位方法 |
CN118243102A (zh) * | 2024-03-25 | 2024-06-25 | 哈尔滨工业大学 | 一种基于有限间断信息的目标位置估计方法 |
Also Published As
Publication number | Publication date |
---|---|
CN107315171B (zh) | 2021-04-20 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN107315171B (zh) | 一种雷达组网目标状态与系统误差联合估计算法 | |
CN110375730B (zh) | 基于imu和uwb融合的室内定位导航系统 | |
US20210373855A1 (en) | Rotation Matrix-Based Factor Graph Cooperative Localization Algorithm | |
CN107300697A (zh) | 基于无人机的运动目标ukf滤波方法 | |
Ullah et al. | Evaluation of Localization by Extended Kalman Filter, Unscented Kalman Filter, and Particle Filter‐Based Techniques | |
WO2022193106A1 (zh) | 一种通过惯性测量参数将gps与激光雷达融合定位的方法 | |
Xu et al. | An improved indoor localization method for mobile robot based on WiFi fingerprint and AMCL | |
CN111783307A (zh) | 一种高超声速飞行器状态估计方法 | |
Zhang et al. | High-precision, limited-beacon-aided AUV localization algorithm | |
CN103344946A (zh) | 一种地基雷达与空中移动平台雷达的实时误差配准方法 | |
Tang et al. | GNSS/inertial navigation/wireless station fusion UAV 3-D positioning algorithm with urban canyon environment | |
Xie et al. | An improved algorithm based on particle filter for 3D UAV target tracking | |
CN115900708A (zh) | 基于gps引导式粒子滤波的机器人多传感器融合定位方法 | |
Hu et al. | A reliable cooperative fusion positioning methodology for intelligent vehicle in non-line-of-sight environments | |
Tang et al. | Vehicle Heterogeneous Multi-Source Information Fusion Positioning Method | |
Tang et al. | Geomagnetic matching cooperative positioning method for unmanned boat cluster based on factor graph | |
CN107463871A (zh) | 一种基于角特征加权的点云匹配方法 | |
CN114236480A (zh) | 一种机载平台传感器系统误差配准算法 | |
CN111736144B (zh) | 一种仅用距离观测的机动转弯目标状态估计方法 | |
CN102707268A (zh) | 机动雷达组网批处理式误差配准器 | |
CN116047495B (zh) | 一种用于三坐标雷达的状态变换融合滤波跟踪方法 | |
CN113076634A (zh) | 一种多机协同无源定位方法、装置及系统 | |
Zhang et al. | UWB/INS-based robust anchor-free relative positioning scheme for UGVs | |
CN110426689A (zh) | 一种基于em-cks的机载多平台多传感器系统误差配准算法 | |
Fang et al. | Model-free approach for sensor network localization with noisy distance measurement |
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 |