CN114895298A - 一种雷达慢速弱机动目标量测修正伯努利滤波方法及装置 - Google Patents

一种雷达慢速弱机动目标量测修正伯努利滤波方法及装置 Download PDF

Info

Publication number
CN114895298A
CN114895298A CN202210357684.5A CN202210357684A CN114895298A CN 114895298 A CN114895298 A CN 114895298A CN 202210357684 A CN202210357684 A CN 202210357684A CN 114895298 A CN114895298 A CN 114895298A
Authority
CN
China
Prior art keywords
target
value
moment
time
current
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
Application number
CN202210357684.5A
Other languages
English (en)
Other versions
CN114895298B (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.)
Xidian University
Original Assignee
Xidian University
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 Xidian University filed Critical Xidian University
Priority to CN202210357684.5A priority Critical patent/CN114895298B/zh
Publication of CN114895298A publication Critical patent/CN114895298A/zh
Application granted granted Critical
Publication of CN114895298B publication Critical patent/CN114895298B/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
    • G01S13/00Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/66Radar-tracking systems; Analogous systems
    • G01S13/72Radar-tracking systems; Analogous systems for two-dimensional tracking, e.g. combination of angle and range tracking, track-while-scan radar
    • G01S13/723Radar-tracking systems; Analogous systems for two-dimensional tracking, e.g. combination of angle and range tracking, track-while-scan radar by using numerical data
    • G01S13/726Multiple target tracking
    • 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/40Means for monitoring or calibrating
    • G01S7/4004Means for monitoring or calibrating of parts of a radar system
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/16Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明涉及一种雷达慢速弱机动目标量测修正伯努利滤波方法及装置,滤波方法包括:获取目标航迹当前时刻的前M个时刻的状态估计值,得到邻近时刻目标的航迹历史值;利用航迹历史值对目标运动模型进行拟合,得到前M个时刻中任一时刻的运动模型参数;利用任一时刻的运动模型参数对当前时刻的目标量测值进行修正,得到当前时刻目标的量测修正值;根据航迹历史值和运动模型对每个时刻目标的量测修正值计算目标量测的协方差矩阵;采用高斯和伯努利滤波器,结合当前时刻目标的量测修正值和协方差矩阵计算目标当前时刻的存在概率,并根据当前时刻的存在概率对目标状态进行更新。该滤波方法能够针对低分辨雷达对弱机动目标的跟踪提供高精度的滤波值。

Description

一种雷达慢速弱机动目标量测修正伯努利滤波方法及装置
技术领域
本发明属于雷达数据处理技术领域,具体涉及一种雷达慢速弱机动目标量测修正伯努利滤波方法及装置。
背景技术
雷达目标跟踪是对回波信号处理之后形成的量测数据进行处理,获得目标的位置和运动状态等参数,并形成目标的运动轨迹,从而是将雷达下发的点迹集合处理输出为目标航迹集合的过程,是雷达数据处理的一个重要环节。跟踪滤波算法的本质是对目标的运动状态进行估计,即从包含噪声的量测数据中提取出目标的状态信息,雷达系统根据每一时刻的状态形成目标运动轨迹,对目标进行评估与预测。因此跟踪滤波算法的优劣直接影响着雷达的安全效能。
现有的雷达目标跟踪系统中,对目标进行跟踪滤波时,主要有以下两种方式。第一种是传统的卡尔曼滤波方法。将状态估计问题视为概率的推理问题,不直接求解目标的状态估计,而是求解目标后验概率密度,将后验概率密度p(xk|yk)取最大值时目标所处的位置xk视作当前时刻目标状态的最优估计。该方法将滤波分为两个过程,第一个过程为预测过程,分别对状态向量和误差协方差矩阵进行一步预测;第二个过程为更新过程,使用当前时刻获得的量测值对状态向量和误差协方差矩阵的预测进行更新。该方法采用最小方差准则进行估计,通过递推的方式实现对目标的更新。但是该方法没有考虑虚警、杂波和目标的漏检,且不适用于非线性非高斯的系统。第二种方法是基于随机有限集的伯努利滤波方法。将目标的出现和消失以及存在杂波等物理现象统一以集合的形式描述,目标消失时为空集,目标出现时为概率密度为pk的单元素集,在杂波环境下也能获得较好的跟踪效果。但是该方法相较于第一种方法只能提高杂波环境下目标跟踪的准确性,无法提高目标的跟踪精度。
低分辨雷达测量精度低,在对慢速运动目标进行跟踪时,由于慢速目标运动速度慢,在探测周期内的位移距离难以超出低分辨雷达的误差范围,使得难以捕捉慢速目标的运动趋势,进一步增大了低测量精度对目标跟踪精度和航迹稳定性的影响,使用常规滤波方法进行跟踪,会导致较大的跟踪误差,甚至可能出现跟丢目标的情况,即常规方法的稳定性不高。
综上,在量测误差较大的情况下对慢速目标进行跟踪时,目标的运动趋势难以捕捉,直接对目标的量测数据进行滤波,或只是简单对量测数据的可靠性进行判断后决定丢弃还是直接进行处理,会导致目标跟踪的精度较低甚至跟丢目标。
发明内容
为了解决现有技术中存在的上述问题,本发明提供了一种雷达慢速弱机动目标量测修正伯努利滤波方法及装置。本发明要解决的技术问题通过以下技术方案实现:
本发明实施例提供了一种雷达慢速弱机动目标量测修正伯努利滤波方法,包括步骤:
S1、获取目标航迹当前时刻的前M个时刻的状态估计值,得到邻近时刻目标的航迹历史值;
S2、利用所述航迹历史值对目标的运动模型进行拟合,得到所述前M个时刻中任一时刻的运动模型参数,所述运动模型为匀速运动模型;
S3、利用所述任一时刻的运动模型参数对所述当前时刻获得的目标量测值进行修正,得到当前时刻目标的量测修正值;
S4、根据所述航迹历史值和所述运动模型对每个时刻目标的量测修正值计算目标量测的协方差矩阵;
S5、采用高斯和伯努利滤波器,结合所述当前时刻目标的量测修正值和所述协方差矩阵计算目标当前时刻的存在概率,并根据所述当前时刻的存在概率对目标状态进行更新。
在本发明的一个实施例中,步骤S2包括:
S21、利用总体最小二乘方法,根据所述航迹历史值中目标邻近时刻航迹的纬度坐标和运动速度计算目标的北向速度估计值和任一时刻的纬度估计值;
S22、利用总体最小二乘方法,根据所述航迹历史值中目标邻近时刻航迹的经度坐标和运动速度计算目标的东向速度估计值和任一时刻的经度估计值。
在本发明的一个实施例中,所述北向速度估计值和所述任一时刻的纬度估计值为:
Figure BDA0003587990560000031
其中,
Figure BDA0003587990560000032
为任一时刻的纬度估计值,
Figure BDA0003587990560000033
为北向速度估计值,AB为系数矩阵,σB为增广矩阵B'=[-B,AB]的最小奇异值,B为目标邻近时刻航迹的纬度坐标,B=[Bk-M,Bk-M+1,…,Bk]T,k为当前时刻,M为当前时刻向后回溯的时刻数,I为单位矩阵;
系数矩阵AB为:
Figure BDA0003587990560000034
其中,T为相邻时刻的时间间隔,k为当前时刻,RB=6356752m为地球短轴半径。
在本发明的一个实施例中,所述东向速度估计值和所述任一时刻的经度估计值为:
Figure BDA0003587990560000035
其中,
Figure BDA0003587990560000036
为任一时刻的经度估计值,
Figure BDA0003587990560000037
为东向速度估计值,AL为系数矩阵,σL为增广矩阵L'=[-L,AL]的最小奇异值,L为目标邻近时刻航迹的经度坐标,L=[Lk-M,Lk-M+1,…,Lk]T,k为当前时刻,M为当前时刻的向后回溯的时刻数,I为单位矩阵;
系数矩阵AL为:
Figure BDA0003587990560000041
其中,T为相邻时刻的时间间隔,k为当前时刻,
Figure BDA0003587990560000042
为目标所处的纬度圆半径,
Figure BDA0003587990560000043
R为地球平均半径,Bk-m为目标k-m时刻所处的纬度。
在本发明的一个实施例中,步骤S3包括:
S31、利用所述北向速度估计值、所述任一时刻的纬度估计值、所述东向速度估计值和所述任一时刻的经度估计值计算当前时刻目标量测坐标的估计值:
Figure BDA0003587990560000044
其中,
Figure BDA0003587990560000045
为当前时刻的纬度坐标,
Figure BDA0003587990560000046
为当前时刻的经度坐标,
Figure BDA0003587990560000047
为任一时刻的纬度估计值,
Figure BDA0003587990560000048
为任一时刻的经度估计值,m为任一时刻距离当前时刻的时刻数,ΔT为每一帧数据的时间间隔,
Figure BDA0003587990560000049
为北向速度估计值,
Figure BDA00035879905600000410
为东向速度估计值,RB为地球平均半径,RL为目标所处的纬度圆半径;
S32、利用所述北向速度估计值和所述东向速度估计值计算当前时刻目标量测速度的估计值:
Figure BDA00035879905600000411
其中,
Figure BDA00035879905600000412
为北向速度估计值,
Figure BDA00035879905600000413
为东向速度估计值。
S33、由所述当前时刻目标量测坐标的估计值和所述当前时刻目标量测速度的估计值得到所述当前时刻目标的量测修正值:
Figure BDA00035879905600000414
其中,
Figure BDA0003587990560000051
为当前时刻目标的量测修正值。
在本发明的一个实施例中,所述目标量测的协方差矩阵为:
Figure BDA0003587990560000052
其中,
Figure BDA0003587990560000053
为目标量测的协方差矩阵,
Figure BDA0003587990560000054
为量测误差均值,
Figure BDA0003587990560000055
E为期望函数,zm为k-M时刻至k时刻的航迹历史值,
Figure BDA0003587990560000056
为k-M时刻至k时刻的量测修正值。
在本发明的一个实施例中,步骤S5包括:
S51、利用上一时刻目标的存在概率预测当前时刻目标的存在概率,利用上一时刻所有高斯成分的均值预测当前时刻每个高斯成分的均值,利用上一时刻所有高斯成分的协方差矩阵预测当前时刻每个高斯成分的协方差矩阵,利用上一时刻存活高斯成分的权值预测当前时刻存活高斯成分的权值和新生高斯成分的权值;
S52、利用所述当前时刻目标的量测修正值和预测的所述当前时刻每个高斯成分的均值更新每个高斯成分的均值,利用预测的所述当前时刻每个高斯成分的协方差矩阵更新每个高斯成分的协方差矩阵,利用所述当前时刻目标的量测修正值和预测的所述当前时刻每个高斯成分的权值更新每个高斯成分的权值,根据预测的所述当前时刻目标的存在概率和更新后的所述每个高斯成分的权值更新目标的存在概率:
Figure BDA0003587990560000057
其中,rk|k-1为预测的当前时刻目标的存在概率,
Figure BDA0003587990560000058
为更新后的高斯成分的权值,PF为雷达的虚警概率;
S53、判断更新后的所述目标的存在概率是否满足预设条件,若满足,则根据更新后的所述每个高斯成分的均值和更新后的所述每个高斯成分的权值设置当前时刻所述目标状态,如不满足,则将当前时刻的所述目标状态更新为上一时刻目标状态的预测。
在本发明的一个实施例中,所述预设条件为:更新后的所述目标的存在概率大于0.5。
在本发明的一个实施例中,根据更新后的所述每个高斯成分的均值和更新后的所述每个高斯成分的权值设置当前时刻所述目标状态,包括,
将所有所述高斯成分的权值中权值最大的所述高斯成分的均值设置为所述当前时刻所述目标状态。
本发明的另一实施例提供了一种雷达慢速弱机动目标量测修正伯努利滤波装置,包括:
航迹历史值获取模块,用于获取目标航迹当前时刻的前M个时刻的状态估计值,得到航迹历史值;
运动模型拟合模块,用于利用所述航迹历史值对目标的运动模型进行拟合,得到所述前M个时刻中任一时刻的运动模型参数,所述运动模型为匀速运动模型;
目标量测值修正模块,用于利用所述任一时刻的运动模型参数对所述当前时刻获得的目标量测值进行修正,得到当前时刻目标的量测修正值;
协方差矩阵量测模块,用于根据所述航迹历史值和所述运动模型对每个时刻目标的量测修正值计算目标量测的协方差矩阵;
目标状态更新模块,用于采用高斯和伯努利滤波器,结合所述当前时刻目标的量测修正值和所述协方差矩阵计算目标当前时刻的存在概率,并根据所述当前时刻的存在概率对目标状态进行更新。
与现有技术相比,本发明的有益效果:
本发明的滤波方法利用弱机动目标邻近时刻运动状态的一致性,使用邻近时刻目标的航迹历史值对当前时刻目标的量测值进行修正,使用修正后的量测值进行滤波,使得目标跟踪的准确性和稳定性得以大幅提升,能够针对低分辨雷达对弱机动目标的跟踪提供高精度的滤波值,能有效提高目标的跟踪精度。
附图说明
图1为本发明实施例提供的一种雷达慢速弱机动目标量测修正伯努利滤波方法的流程示意图;
图2为本发明实施例提供的另一种雷达慢速弱机动目标量测修正伯努利滤波方法的流程示意图;
图3为本发明实施例提供的一种伯努利滤波方法仿真对目标航迹进行滤波处理的结果示意图;
图4为本发明实施例提供的一种雷达慢速弱机动目标量测修正伯努利滤波方法对目标航迹进行滤波处理的结果示意图;
图5为本发明实施例提供的一种伯努利滤波方法和雷达慢速弱机动目标量测修正伯努利滤波方法的跟踪误差对比图;
图6为本发明实施例提供的一种伯努利滤波方法仿真对实测数据进行滤波处理的结果示意图;
图7为本发明实施例提供的一种雷达慢速弱机动目标量测修正伯努利滤波方法对实测数据进行处理的结果示意图;
图8为本发明实施例提供的一种雷达慢速弱机动目标量测修正伯努利滤波装置的示意图。
具体实施方式
下面结合具体实施例对本发明做进一步详细的描述,但本发明的实施方式不限于此。
实施例一
请参见图1和图2,图1为本发明实施例提供的一种雷达慢速弱机动目标量测修正伯努利滤波方法的流程示意图,图2为本发明实施例提供的另一种雷达慢速弱机动目标量测修正伯努利滤波方法的流程示意图。该滤波方法适用于低分辨雷达对慢速弱机动目标的跟踪,该方法包括步骤:
S1、获取目标航迹当前时刻的前M个时刻的状态估计值,得到邻近时刻目标的航迹历史值。
具体的,根据高斯和伯努利滤波方法从目标航迹当前时刻k向后回溯M个时刻,获取前M个时刻的状态估计值,得到邻近时刻即目标前M步的航迹历史值{Yk-m},其中,m=1,2,…,M,k-m表示时刻,M的取值范围可以为4~10。航迹历史值{Yk-m}包括目标历史时刻的经度Lk-m、目标历史时刻的纬度Bk-m和运动速度vk-m
S2、利用所述航迹历史值对目标的运动模型进行拟合,得到所述前M个时刻中任一时刻的运动模型参数,所述运动模型为匀速运动模型。
本实施例中,前M个时刻中任一时刻的运动模型参数均可用于对当前时刻的量测值进行修正,即运动模型参数可以为任一时刻的参数。进一步的,运动模型参数包括任一时刻目标的经度坐标、任一时刻目标的纬度坐标和任一时刻目标的运动速度;由于运动模型为匀速运动模型,因此,k-M到k时刻的运动速度是相同的。
本实施例以k-M时刻为例进行说明,具体包括步骤:
S21、利用总体最小二乘方法,根据所述航迹历史值中目标邻近时刻航迹的纬度坐标和运动速度计算目标的北向速度估计值和任一时刻的纬度估计值。
具体的,航迹历史值{Yk-m}中,其目标邻近时刻航迹的纬度坐标为B=[Bk-M,Bk-M+1,…,Bk]T,则目标的北向速度估计值
Figure BDA0003587990560000081
和k-M时刻的纬度估计值
Figure BDA0003587990560000082
由下式计算得出:
Figure BDA0003587990560000083
其中,
Figure BDA0003587990560000084
为k-M时刻的纬度估计值,
Figure BDA0003587990560000085
为北向速度估计值,AB为系数矩阵,σB为增广矩阵B'=[-B,AB]的最小奇异值,B为目标邻近时刻航迹的纬度坐标,B=[Bk-M,Bk-M+1,…,Bk]T,k为当前时刻,M为当前时刻向后回溯的时刻数,I为单位矩阵。
系数矩阵AB为:
Figure BDA0003587990560000091
其中,T为相邻时刻的时间间隔,k为当前时刻,RB为地球短轴半径。
S22、利用总体最小二乘方法,根据所述航迹历史值中目标邻近时刻航迹的经度坐标和运动速度计算目标的东向速度估计值和任一时刻的经度估计值。
具体的,航迹历史值{Yk-m}中,其目标邻近时刻航迹的经度值为L=[Lk-M,Lk-M+1,…,Lk]T,则目标的东向速度估计值
Figure BDA0003587990560000092
和k-M时刻的经度估计值为:
Figure BDA0003587990560000093
其中,
Figure BDA0003587990560000094
为k-M时刻的经度估计值,
Figure BDA0003587990560000095
为东向速度估计值,AL为系数矩阵,σL为增广矩阵L'=[-L,AL]的最小奇异值,L为目标邻近时刻航迹的经度坐标,L=[Lk-M,Lk-M+1,…,Lk]T,k为当前时刻,M为当前时刻向后回溯的时刻数,I为单位矩阵。
系数矩阵AL为:
Figure BDA0003587990560000096
其中,T为相邻时刻的时间间隔,k为当前时刻,
Figure BDA0003587990560000097
为目标所处的纬度圆半径,
Figure BDA0003587990560000098
R为地球平均半径,Bk-m为k-m时刻目标所处的纬度值。
S3、利用所述任一时刻的运动模型参数对所述当前时刻获得的目标量测值进行修正,得到当前时刻目标的量测修正值。
本实施例中,当前时刻目标的量测修正值包括:目标的经度坐标估计值、目标的纬度坐标估计值和目标运动速度估计值。
具体包括步骤:
S31、利用所述北向速度估计值、所述任一时刻的纬度估计值、所述东向速度估计值和所述任一时刻的经度估计值计算当前时刻目标量测坐标的估计值。
具体的,以任一时刻为k-M时刻为例,利用步骤S2估计得到的运动模型参数:k-M时刻的纬度估计值
Figure BDA0003587990560000101
k-M时刻的经度估计值
Figure BDA0003587990560000102
北向速度估计值
Figure BDA0003587990560000103
东向速度估计值
Figure BDA0003587990560000104
计算当前时刻目标量测坐标的估计值:
Figure BDA0003587990560000105
其中,
Figure BDA0003587990560000106
为当前时刻的纬度坐标,
Figure BDA0003587990560000107
为当前时刻的经度坐标,
Figure BDA0003587990560000108
为k-M时刻的纬度估计值,
Figure BDA0003587990560000109
为k-M时刻的经度估计值,M为向后回溯的时刻数,ΔT为每一帧数据的时间间隔,
Figure BDA00035879905600001010
为北向速度估计值,
Figure BDA00035879905600001011
为东向速度估计值,RB为地球平均半径,RL为目标所处的纬度圆半径。
S32、利用所述北向速度估计值和所述东向速度估计值计算当前时刻目标量测速度的估计值。
具体的,利用北向速度估计值
Figure BDA00035879905600001012
东向速度估计值
Figure BDA00035879905600001013
计算当前时刻目标量测速度的估计值:
Figure BDA00035879905600001014
其中,
Figure BDA00035879905600001015
为北向速度估计值,
Figure BDA00035879905600001016
为东向速度估计值。
S33、由所述当前时刻目标量测坐标的估计值和所述当前时刻目标量测速度的估计值得到所述当前时刻目标的量测修正值。
具体的,当前时刻目标的量测修正值即为当前时刻目标量测坐标的估计值
Figure BDA00035879905600001017
和当前时刻目标量测速度的估计值
Figure BDA00035879905600001018
形成的矩阵:
Figure BDA00035879905600001019
其中,
Figure BDA00035879905600001020
为当前时刻目标的量测修正值。
进一步的,通过上述方法可以获得每个时刻目标的量测修正值
Figure BDA0003587990560000111
S4、根据所述航迹历史值和所述运动模型对每个时刻目标的量测修正值计算目标量测的协方差矩阵。
具体的,首先,由邻近时刻目标的航迹历史值{Yk-m}和运动模型对每个时刻目标的量测修正值
Figure BDA0003587990560000112
计算得到所有时刻的量测误差均值
Figure BDA0003587990560000113
Figure BDA0003587990560000114
其中,E为期望函数,zm为k-M时刻至k时刻的航迹历史值,
Figure BDA0003587990560000115
为k-M时刻至k时刻的量测修正值。
然后,由量测误差均值
Figure BDA0003587990560000116
计算所有时刻目标量测的协方差矩阵
Figure BDA0003587990560000117
Figure BDA0003587990560000118
其中,
Figure BDA0003587990560000119
为目标量测的协方差矩阵,
Figure BDA00035879905600001110
为量测误差均值。
S5、采用高斯和伯努利滤波器,结合所述当前时刻目标的量测修正值和所述协方差矩阵计算目标当前时刻的存在概率,并根据所述当前时刻的存在概率对目标状态进行更新。
具体的,高斯和伯努利滤波器将目标状态视为多个服从高斯分布状态的叠加,目标的概率密度由L个高斯成分表示。高斯和伯努利滤波器分两个步骤实现,分别为预测步和更新步,其具体实现步骤为:
S51、预测步。
具体的,利用上一时刻目标的存在概率rk-1|k-1预测当前时刻目标的存在概率:
rk|k-1=pb(1-rk-1|k-1)+psrk-1|k-1 (10)
其中,rk|k-1为当前时刻目标的存在概率,pb为目标的新生概率,rk-1|k-1为上一时刻目标的存在概率,ps为目标的存活概率。
利用上一时刻所有高斯成分的均值
Figure BDA00035879905600001111
预测当前时刻每个高斯成分的均值,其中,所有高斯成分是指L个高斯成分。预测得到的当前时刻每个高斯成分的均值为:
Figure BDA0003587990560000121
其中,l=1,2,…,L,Fk-1为状态转移矩阵。
利用上一时刻所有高斯成分的协方差矩阵
Figure BDA0003587990560000122
预测当前时刻每个高斯成分的协方差矩阵。其中,所有高斯成分是指L个高斯成分。预测得到的当前时刻每个高斯成分的协方差矩阵为:
Figure BDA0003587990560000123
其中,l=1,2,…,L,Fk-1为状态转移矩阵,Qk-1为状态噪声协方差矩阵,F′k-1为状态转移矩阵的转置。
利用上一时刻存活高斯成分的权值预测当前时刻存活高斯成分的权值和新生高斯成分的权值。预测的存活高斯成分的权值为:
Figure BDA0003587990560000124
其中,ps为高斯成分的存活概率,rk-1为目标上一时刻的存在概率,
Figure BDA0003587990560000125
为上一时刻存活高斯成分的权值。
预测的新生高斯成分的权值为:
Figure BDA0003587990560000126
其中,pR为高斯成分的新生概率,rk-1为目标上一时刻的存在概率,
Figure BDA0003587990560000127
为上一时刻存活高斯成分的权值。
S52、更新步。
利用所述当前时刻目标的量测修正值
Figure BDA0003587990560000128
和预测的所述当前时刻每个高斯成分的均值更新每个高斯成分的均值。具体的,第l个高斯成分的均值更新为:
Figure BDA0003587990560000129
其中,
Figure BDA00035879905600001210
为预测的当前时刻第l个高斯成分的均值,
Figure BDA00035879905600001211
为卡尔曼增益,
Figure BDA00035879905600001212
为预测得到的当前时刻第l个高斯成分的协方差矩阵,
Figure BDA0003587990560000131
为量测协方差矩阵,
Figure BDA0003587990560000132
Rk为量测噪声协方差矩阵,
Figure BDA0003587990560000133
为当前时刻目标的量测修正值,
Figure BDA0003587990560000134
Hk为量测矩阵。
利用预测的所述当前时刻每个高斯成分的协方差矩阵更新每个高斯成分的协方差矩阵。具体的,高斯成分的协方差矩阵更新为:
Figure BDA0003587990560000135
其中,
Figure BDA0003587990560000136
为预测的当前时刻第i个高斯成分的协方差矩阵,
Figure BDA0003587990560000137
为卡尔曼增益,Hk为量测矩阵。
利用所述当前时刻目标的量测修正值和预测的所述当前时刻每个高斯成分的权值更新每个高斯成分的权值。具体的,高斯成分权值更新为:
Figure BDA0003587990560000138
其中,pD为雷达的检测概率,
Figure BDA0003587990560000139
为量测的高斯似然函数,
Figure BDA00035879905600001310
为预测的当前时刻第l个高斯成分的权值。
高斯似然函数的表达式为:
Figure BDA00035879905600001311
其中,
Figure BDA00035879905600001312
为当前时刻目标的量测修正值,Sk|k-1为量测修正值对应的误差协方差矩阵,
Figure BDA00035879905600001313
根据预测的所述当前时刻目标的存在概率和更新后的所述每个高斯成分的权值更新目标的存在概率:
Figure BDA00035879905600001314
其中,rk|k-1为预测的当前时刻目标的存在概率,
Figure BDA00035879905600001315
为更新后的高斯成分的权值,PF为雷达的虚警概率。
在滤波过程中,由于高斯成分会不断的增加,因此需要对高斯成分进行以下处理:首先,对权值低于一定阈值δ的高斯成分进行修剪以将权值较小的高斯成分去除掉,修剪阈值δ一般取1e-5;然后,对修剪后的高斯成分进行合并,具体的,对距离小于一定阈值ε的高斯成分进行合并以保留两个距离较近中的一个高斯成分,合并阈值ε一般取4,第l个高斯和第i个高斯成分之间的距离定义为
Figure BDA0003587990560000141
最后,经过合并后只保留一定数目且具有较高权值的高斯成分。
S53、判断更新后的所述目标的存在概率是否满足预设条件,若满足,则根据更新后的所述每个高斯成分的均值和更新后的所述每个高斯成分的权值设置当前时刻所述目标状态,如不满足,则将当前时刻的所述目标状态更新为上一时刻目标状态的预测。
具体的,所述预设条件为:更新后的所述目标的存在概率大于0.5。进一步的,当更新后的目标的存在概率>0.5,则将所有高斯成分的权值中权值最大的高斯成分的均值设置为当前时刻所述目标状态,即目标当前时刻的更新滤波状态xk|k为权值最大的高斯高斯成分均值
Figure BDA0003587990560000142
当更新后的目标的存在概率≤0.5,则目标当前时刻的更新滤波状态为上一时刻目标状态xk|k-1的预测xk|k=Fkxk|k-1,其中,Fk为状态转移矩阵。
本实施例滤波方法的效果通过以下对仿真和实测数据的处理进一步说明。
1、实验环境与条件
数据率:40s/帧;目标数:1个;目标运动速度:14m/s;测距误差:10km;数据帧数:80帧。
2、实现内容与结果
实验1,用伯努利滤波方法仿真对目标航迹进行滤波处理,结果如图3所示,图3为本发明实施例提供的一种伯努利滤波方法仿真对目标航迹进行滤波处理的结果示意图。
实验2,用本实施例方法仿真对目标航迹进行滤波处理,结果如图4所示,图4为本发明实施例提供的一种雷达慢速弱机动目标量测修正伯努利滤波方法对目标航迹进行滤波处理的结果示意图。
比较图3和图4可见,伯努利滤波方法的滤波值随量测值波动较大,目标运动趋势不明显,跟踪效果较差。而本实施例的滤波值随量测值得波动很小,目标运动趋势明显,跟踪效果得到了明显得改善。
实验3,进行500次蒙特卡洛实验,对比伯努利滤波方法和本实施例方法的跟踪误差,结果如图5所示,图5为本发明实施例提供的一种伯努利滤波方法和雷达慢速弱机动目标量测修正伯努利滤波方法的跟踪误差对比图。从图中可以看出针对大量测误差下对慢速目标的跟踪,相比于伯努利滤波,本实施例的滤波误差有明显的减小,能够明显改善目标跟踪的精度。
实验4,用伯努利滤波方法对实测数据进行处理,结果如图6所示,图6为本发明实施例提供的一种伯努利滤波方法仿真对实测数据进行滤波处理的结果示意图。
实验5,用本发明方法对实测数据进行处理,结果如图7所示,图7为本发明实施例提供的一种雷达慢速弱机动目标量测修正伯努利滤波方法对实测数据进行处理的结果示意图。
图6和图7中目标的真实值来源于船舶自动识别系统,可作为真实值的参考,对比两图可以看出,在实际应用中本实施例也能够获得较好的跟踪效果,本实施例相较于伯努利滤波方法,拥有更好的跟踪效果,跟踪稳定性得到了明显的改善。
综上所述,本实施例的滤波方法利用弱机动目标邻近时刻运动状态的一致性,使用邻近时刻目标的航迹历史值对当前时刻目标的量测值进行修正,使用修正后的量测值进行滤波,使得目标跟踪的准确性和稳定性得以大幅提升,能够针对低分辨雷达对弱机动目标的跟踪提供高精度的滤波值,能有效提高目标的跟踪精度。
实施例二
在实施例一的基础上,请参见图8,图8为本发明实施例提供的一种雷达慢速弱机动目标量测修正伯努利滤波装置的示意图。该滤波装置包括:航迹历史值获取模块、运动模型拟合模块、目标量测值修正模块、协方差矩阵量测模块和目标状态更新模块。
具体的,航迹历史值获取模块用于获取目标航迹当前时刻的前M个时刻的状态估计值,得到航迹历史值。运动模型拟合模块用于利用所述航迹历史值对目标的运动模型进行拟合,得到所述前M个时刻中任一时刻的运动模型参数,所述运动模型为匀速运动模型。目标量测值修正模块用于利用所述任一时刻的运动模型参数对所述当前时刻获得的目标量测值进行修正,得到当前时刻目标的量测修正值。协方差矩阵量测模块用于根据所述航迹历史值和所述运动模型对每个时刻目标的量测修正值计算目标量测的协方差矩阵。目标状态更新模块用于采用高斯和伯努利滤波器,结合所述当前时刻目标的量测修正值和所述协方差矩阵计算目标当前时刻的存在概率,并根据所述当前时刻的存在概率对目标状态进行更新。
上述各个模块的具体实施方式请参见实施例一,本实施例不再赘述。
以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明,不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说,在不脱离本发明构思的前提下,还可以做出若干简单推演或替换,都应当视为属于本发明的保护范围。

Claims (10)

1.一种雷达慢速弱机动目标量测修正伯努利滤波方法,其特征在于,包括步骤:
S1、获取目标航迹当前时刻的前M个时刻的状态估计值,得到邻近时刻目标的航迹历史值;
S2、利用所述航迹历史值对目标的运动模型进行拟合,得到所述前M个时刻中任一时刻的运动模型参数,所述运动模型为匀速运动模型;
S3、利用所述任一时刻的运动模型参数对所述当前时刻获得的目标量测值进行修正,得到当前时刻目标的量测修正值;
S4、根据所述航迹历史值和所述运动模型对每个时刻目标的量测修正值计算目标量测的协方差矩阵;
S5、采用高斯和伯努利滤波器,结合所述当前时刻目标的量测修正值和所述协方差矩阵计算目标当前时刻的存在概率,并根据所述当前时刻的存在概率对目标状态进行更新。
2.根据权利要求1所述的雷达慢速弱机动目标量测修正伯努利滤波方法,其特征在于,步骤S2包括:
S21、利用总体最小二乘方法,根据所述航迹历史值中目标邻近时刻航迹的纬度坐标和运动速度计算目标的北向速度估计值和任一时刻的纬度估计值;
S22、利用总体最小二乘方法,根据所述航迹历史值中目标邻近时刻航迹的经度坐标和运动速度计算目标的东向速度估计值和任一时刻的经度估计值。
3.根据权利要求2所述的雷达慢速弱机动目标量测修正伯努利滤波方法,其特征在于,所述北向速度估计值和所述任一时刻的纬度估计值为:
Figure FDA0003587990550000011
其中,
Figure FDA0003587990550000012
为任一时刻的纬度估计值,
Figure FDA0003587990550000013
为北向速度估计值,AB为系数矩阵,σB为增广矩阵B'=[-B,AB]的最小奇异值,B为目标邻近时刻航迹的纬度坐标,B=[Bk-M,Bk-M+1,…,Bk]T,k为当前时刻,M为当前时刻向后回溯的时刻数,I为单位矩阵;
系数矩阵AB为:
Figure FDA0003587990550000021
其中,T为相邻时刻的时间间隔,k为当前时刻,RB为地球短轴半径。
4.根据权利要求2所述的雷达慢速弱机动目标量测修正伯努利滤波方法,其特征在于,所述东向速度估计值和所述任一时刻的经度估计值为:
Figure FDA0003587990550000022
其中,
Figure FDA0003587990550000023
为任一时刻的经度估计值,
Figure FDA0003587990550000024
为东向速度估计值,AL为系数矩阵,σL为增广矩阵L'=[-L,AL]的最小奇异值,L为目标邻近时刻航迹的经度坐标,L=[Lk-M,Lk-M+1,…,Lk]T,k为当前时刻,M为当前时刻向后回溯的时刻数,I为单位矩阵;
系数矩阵AL为:
Figure FDA0003587990550000025
其中,T为相邻时刻的时间间隔,k为当前时刻,
Figure FDA0003587990550000026
为目标所处的纬度圆半径,
Figure FDA0003587990550000027
R为地球平均半径,Bk-m为k-m时刻目标所在的纬度。
5.根据权利要求2所述的雷达慢速弱机动目标量测修正伯努利滤波方法,其特征在于,步骤S3包括:
S31、利用所述北向速度估计值、所述任一时刻的纬度估计值、所述东向速度估计值和所述任一时刻的经度估计值计算当前时刻目标量测坐标的估计值:
Figure FDA0003587990550000031
其中,
Figure FDA0003587990550000032
为当前时刻的纬度坐标,
Figure FDA0003587990550000033
为当前时刻的经度坐标,
Figure FDA0003587990550000034
为任一时刻的纬度估计值,
Figure FDA0003587990550000035
为任一时刻的经度估计值,m为任一时刻距离当前时刻的时刻数,ΔT为每一帧数据的时间间隔,
Figure FDA0003587990550000036
为北向速度估计值,
Figure FDA0003587990550000037
为东向速度估计值,RB为地球短轴半径,RL为目标所处的纬度圆半径;
S32、利用所述北向速度估计值和所述东向速度估计值计算当前时刻目标量测速度的估计值:
Figure FDA0003587990550000038
其中,
Figure FDA0003587990550000039
为北向速度估计值,
Figure FDA00035879905500000310
为东向速度估计值;
S33、由所述当前时刻目标量测坐标的估计值和所述当前时刻目标量测速度的估计值得到所述当前时刻目标的量测修正值:
Figure FDA00035879905500000311
其中,
Figure FDA00035879905500000312
为当前时刻目标的量测修正值。
6.根据权利要求1所述的雷达慢速弱机动目标量测修正伯努利滤波方法,其特征在于,所述目标量测的协方差矩阵为:
Figure FDA00035879905500000313
其中,
Figure FDA00035879905500000314
为目标量测的协方差矩阵,
Figure FDA00035879905500000315
为量测误差均值,
Figure FDA00035879905500000316
E为期望函数,zm为k-M时刻至k时刻的航迹历史值,
Figure FDA00035879905500000317
为k-M时刻至k时刻的量测修正值。
7.根据权利要求1所述的雷达慢速弱机动目标量测修正伯努利滤波方法,其特征在于,步骤S5包括:
S51、利用上一时刻目标的存在概率预测当前时刻目标的存在概率,利用上一时刻所有高斯成分的均值预测当前时刻每个高斯成分的均值,利用上一时刻所有高斯成分的协方差矩阵预测当前时刻每个高斯成分的协方差矩阵,利用上一时刻存活高斯成分的权值预测当前时刻存活高斯成分的权值和新生高斯成分的权值;
S52、利用所述当前时刻目标的量测修正值和预测的所述当前时刻每个高斯成分的均值更新每个高斯成分的均值,利用预测的所述当前时刻每个高斯成分的协方差矩阵更新每个高斯成分的协方差矩阵,利用所述当前时刻目标的量测修正值和预测的所述当前时刻每个高斯成分的权值更新每个高斯成分的权值,根据预测的所述当前时刻目标的存在概率和更新后的所述每个高斯成分的权值更新目标的存在概率:
Figure FDA0003587990550000041
其中,rk|k-1为预测的当前时刻目标的存在概率,
Figure FDA0003587990550000042
为更新后的高斯成分的权值,PF为雷达的虚警概率;
S53、判断更新后的所述目标的存在概率是否满足预设条件,若满足,则根据更新后的所述每个高斯成分的均值和更新后的所述每个高斯成分的权值设置当前时刻所述目标状态,如不满足,则将当前时刻的所述目标状态更新为上一时刻目标状态的预测。
8.根据权利要求7所述的雷达慢速弱机动目标量测修正伯努利滤波方法,其特征在于,所述预设条件为:更新后的所述目标的存在概率大于0.5。
9.根据权利要求7所述的雷达慢速弱机动目标量测修正伯努利滤波方法,其特征在于,根据更新后的所述每个高斯成分的均值和更新后的所述每个高斯成分的权值设置当前时刻所述目标状态,包括,
将所有所述高斯成分的权值中权值最大的所述高斯成分的均值设置为所述当前时刻所述目标状态。
10.一种雷达慢速弱机动目标量测修正伯努利滤波装置,其特征在于,包括:
航迹历史值获取模块,用于获取目标航迹当前时刻的前M个时刻的状态估计值,得到航迹历史值;
运动模型拟合模块,用于利用所述航迹历史值对目标的运动模型进行拟合,得到所述前M个时刻中任一时刻的运动模型参数,所述运动模型为匀速运动模型;
目标量测值修正模块,用于利用所述任一时刻的运动模型参数对所述当前时刻获得的目标量测值进行修正,得到当前时刻目标的量测修正值;
协方差矩阵量测模块,用于根据所述航迹历史值和所述运动模型对每个时刻目标的量测修正值计算目标量测的协方差矩阵;
目标状态更新模块,用于采用高斯和伯努利滤波器,结合所述当前时刻目标的量测修正值和所述协方差矩阵计算目标当前时刻的存在概率,并根据所述当前时刻的存在概率对目标状态进行更新。
CN202210357684.5A 2022-04-08 2022-04-08 一种雷达慢速弱机动目标量测修正伯努利滤波方法及装置 Active CN114895298B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210357684.5A CN114895298B (zh) 2022-04-08 2022-04-08 一种雷达慢速弱机动目标量测修正伯努利滤波方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210357684.5A CN114895298B (zh) 2022-04-08 2022-04-08 一种雷达慢速弱机动目标量测修正伯努利滤波方法及装置

Publications (2)

Publication Number Publication Date
CN114895298A true CN114895298A (zh) 2022-08-12
CN114895298B CN114895298B (zh) 2024-07-12

Family

ID=82715046

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210357684.5A Active CN114895298B (zh) 2022-04-08 2022-04-08 一种雷达慢速弱机动目标量测修正伯努利滤波方法及装置

Country Status (1)

Country Link
CN (1) CN114895298B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117724087A (zh) * 2024-02-07 2024-03-19 中国人民解放军海军航空大学 雷达多目标跟踪双标签多伯努利滤波算法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104766320A (zh) * 2015-04-02 2015-07-08 西安电子科技大学 阈值化量测下的多伯努利滤波弱目标检测与跟踪方法
CN109633599A (zh) * 2019-01-29 2019-04-16 中国人民解放军空军预警学院 一种机载预警雷达多目标跟踪方法
CN109655826A (zh) * 2018-12-16 2019-04-19 成都汇蓉国科微系统技术有限公司 一种低慢小目标轨迹滤波方法及装置
WO2020007487A1 (en) * 2018-07-06 2020-01-09 Bayerische Motoren Werke Aktiengesellschaft Object tracking based on multiple measurement hypotheses
CN113300690A (zh) * 2021-02-05 2021-08-24 苏州经贸职业技术学院 一种有色量测噪声下鲁棒泊松多伯努利滤波方法
CN113866755A (zh) * 2021-07-20 2021-12-31 桂林电子科技大学 基于多伯努利滤波的雷达微弱起伏目标检测前跟踪算法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104766320A (zh) * 2015-04-02 2015-07-08 西安电子科技大学 阈值化量测下的多伯努利滤波弱目标检测与跟踪方法
WO2020007487A1 (en) * 2018-07-06 2020-01-09 Bayerische Motoren Werke Aktiengesellschaft Object tracking based on multiple measurement hypotheses
CN109655826A (zh) * 2018-12-16 2019-04-19 成都汇蓉国科微系统技术有限公司 一种低慢小目标轨迹滤波方法及装置
CN109633599A (zh) * 2019-01-29 2019-04-16 中国人民解放军空军预警学院 一种机载预警雷达多目标跟踪方法
CN113300690A (zh) * 2021-02-05 2021-08-24 苏州经贸职业技术学院 一种有色量测噪声下鲁棒泊松多伯努利滤波方法
CN113866755A (zh) * 2021-07-20 2021-12-31 桂林电子科技大学 基于多伯努利滤波的雷达微弱起伏目标检测前跟踪算法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
WEI LI-XING: "A Labeled GM-CBMeMBer Filter with Adaptive One-Step Track Initiation", 《ELECTRONICS OPTICS & CONTROL》, 8 August 2019 (2019-08-08) *
ZHOU, XQ: "An Experimental Multi-Target Tracking of AM Radio-Based Passive Bistatic Radar System via Multi-Static Doppler Shifts", 《SENSORS》, 8 October 2021 (2021-10-08) *
冯燕: "基于标签多伯努利滤波器的非线性多机动目标跟踪算法研究", 《中国优秀硕士学位论文全文数据库 信息科技辑》, 15 January 2020 (2020-01-15) *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117724087A (zh) * 2024-02-07 2024-03-19 中国人民解放军海军航空大学 雷达多目标跟踪双标签多伯努利滤波算法
CN117724087B (zh) * 2024-02-07 2024-05-28 中国人民解放军海军航空大学 雷达多目标跟踪双标签多伯努利滤波算法

Also Published As

Publication number Publication date
CN114895298B (zh) 2024-07-12

Similar Documents

Publication Publication Date Title
CN110823217B (zh) 一种基于自适应联邦强跟踪滤波的组合导航容错方法
CN110061716B (zh) 一种基于最小二乘和多重渐消因子的改进kalman滤波方法
CN110161543B (zh) 一种基于卡方检验的部分粗差抗差自适应滤波方法
CN109425876B (zh) 一种提高定位精度的改进卡尔曼滤波方法
CN108801131B (zh) 北斗高频变形监测数据的滤波方法及系统
CN111257865B (zh) 一种基于线性伪量测模型的机动目标多帧检测跟踪方法
CN109031229B (zh) 一种杂波环境下目标跟踪的概率假设密度方法
CN110376582B (zh) 自适应gm-phd的机动目标跟踪方法
CN114895298B (zh) 一种雷达慢速弱机动目标量测修正伯努利滤波方法及装置
CN110889862A (zh) 一种网络传输攻击环境中多目标跟踪的组合测量方法
CN110501686A (zh) 基于一种新型自适应高阶无迹卡尔曼滤波的状态估计方法
CN111259332B (zh) 一种杂波环境下的模糊数据关联方法及多目标跟踪方法
CN113486960A (zh) 基于长短时记忆神经网络的无人机跟踪方法、装置、存储介质及计算机设备
CN110233608B (zh) 一种基于权值自适应的粒子滤波方法和雷达系统
CN112328959A (zh) 一种基于自适应扩展卡尔曼概率假设密度滤波器的多目标跟踪方法
CN113191427B (zh) 一种多目标车辆跟踪方法及相关装置
CN109188422B (zh) 一种基于lu分解的卡尔曼滤波目标跟踪方法
CN115291253B (zh) 一种基于残差检测的车辆定位完好性监测方法及系统
CN112305513A (zh) 一种传感器测量参数修正方法及系统
CN115420290B (zh) 非合作目标机动检测方法、装置、设备及计算机存储介质
CN116520311A (zh) 一种基于glmb的自适应航迹起始方法
CN110007298B (zh) 一种目标超前预测跟踪方法
CN110716219A (zh) 一种提高定位解算精度的方法
CN110515069B (zh) 一种用于分布式目标跟踪的自适应一致性信息滤波方法
CN109117965A (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