CN109509207B - A method for seamless tracking of point and extended targets - Google Patents

A method for seamless tracking of point and extended targets Download PDF

Info

Publication number
CN109509207B
CN109509207B CN201811338016.8A CN201811338016A CN109509207B CN 109509207 B CN109509207 B CN 109509207B CN 201811338016 A CN201811338016 A CN 201811338016A CN 109509207 B CN109509207 B CN 109509207B
Authority
CN
China
Prior art keywords
target
measurement
state
covariance
poisson
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
CN201811338016.8A
Other languages
Chinese (zh)
Other versions
CN109509207A (en
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.)
University of Electronic Science and Technology of China
Original Assignee
University of Electronic Science and Technology of China
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 University of Electronic Science and Technology of China filed Critical University of Electronic Science and Technology of China
Priority to CN201811338016.8A priority Critical patent/CN109509207B/en
Publication of CN109509207A publication Critical patent/CN109509207A/en
Application granted granted Critical
Publication of CN109509207B publication Critical patent/CN109509207B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/20—Analysis of motion
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/20—Analysis of motion
    • G06T7/207—Analysis of motion for motion estimation over a hierarchy of resolutions
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/20—Analysis of motion
    • G06T7/277—Analysis of motion involving stochastic approaches, e.g. using Kalman filters

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Multimedia (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了一种对点目标和扩展目标进行无缝跟踪的方法，属于高分辨传感器目标跟踪技术领域。本发明所述方法用可变维度的高斯过程GP建模目标的可变外形轮廓，目标轮廓占据的各个传感器分辨单元记为量测源；若相对于传感器观测的目标越小，目标轮廓上的量测源个数越少，反之亦然。本发明利用量测源个数的估计值大小在线调整GP模型的半径数。本发明能适应ET的外形变化和ET与PT之间的相互转换，无缝地跟踪多个ET和PT，并保持较好的跟踪性能。当目标为ET时，ET‑GP‑PMHT跟踪和输出目标的位置和外形；当目标为PT时，ET‑GP‑PMHT仅跟踪和输出目标的位置。此外，该方法计算复杂度与量测数、半径数和目标数相关。当ET外形变小时，采用的GP模型半径数变少，计算复杂度减小。

Figure 201811338016

The invention discloses a method for seamless tracking of point targets and extended targets, belonging to the technical field of high-resolution sensor target tracking. The method of the present invention uses a variable-dimensional Gaussian process GP to model the variable contour of the target, and each sensor resolution unit occupied by the target contour is recorded as a measurement source; if the target observed relative to the sensor is smaller, the The less the number of measurement sources, and vice versa. The invention adjusts the number of radii of the GP model online by using the estimated value of the number of measurement sources. The invention can adapt to the shape change of ET and the mutual conversion between ET and PT, seamlessly track multiple ETs and PTs, and maintain better tracking performance. When the target is ET, ET‑GP‑PMHT tracks and outputs the position and shape of the target; when the target is PT, ET‑GP‑PMHT only tracks and outputs the position of the target. In addition, the computational complexity of this method is related to the number of measurements, the number of radii and the number of targets. When the ET shape becomes smaller, the number of GP model radii used becomes smaller, and the computational complexity is reduced.

Figure 201811338016

Description

Method for seamless tracking of point target and extended target
Technical Field
The invention belongs to the technical field of high-resolution sensor target tracking, and particularly relates to a method for seamlessly tracking a point target and an extended target.
Background
In existing tracking algorithms, the target is modeled as a point source. With the improvement of the radar resolution or the fact that the target is close to the sensor, the target occupies a plurality of resolution units, a plurality of measurements are generated, and a Point Target (PT) model is not applicable any more, so that the problem of extending the target (ET) is solved, and therefore more and more documents for researching an ET tracking (ETT) algorithm appear.
In existing ETT studies, the ET state is modeled as two parts, the motion state and the target profile. In recent years, various ET shape modeling methods have been proposed. The Random Matrix (RM) model is modeled into an ellipse by adopting a symmetrical positive definite matrix, the motion state of the ellipse meets Gaussian distribution, and the target shape meets inverse Wishart distribution. Modeling the outline of an object with only an ellipse is not suitable for all objects, so a non-elliptical model is more suitable for an ET of arbitrary shape.
An intuitive idea is to model the target as a combination of ellipses, another method is to use a star-convex shape method, which models the unknown extended target shape as a finite number of unknown radius functions, a random hyper-surface model (RHM) based star-convex shape method defines the radius functions in the frequency domain, and a fourier series expansion method is used to parameterize the radius. The maximum Expectation (EM) method can be used in the RHM framework, and the advantages of the EM based on the recursive Gaussian RHM method are also researched and verified. The star-convex shape method based on the Gaussian Process (GP) model can model a target radius function in a spatial domain, namely model the shape of a target, and maintains the uncertainty of the unobserved part of the target. It is flexible enough to be used to represent a variety of shapes. In the GP framework, Extended Kalman Filtering (EKF) and Particle Filtering (PF) can track a single extended target, and in order to track multiple targets, label-multi-bernoulli (LMB) filtering and gaussian mixture probability hypothesis density (GM-PHD) filtering are proposed in heterogeneous multi-sensor scenarios.
Multiple ETs in an actual ETT scenario may have different profiles, far or near the sensor, and the target may occupy one or more resolution cells appearing as PT or ET. The same target, near the sensor, may appear as ET and far as PT. In some monitoring intervals, ET and PT may be present simultaneously, and the size of the extended target relative to the sensor may change as the target splits, merges, deflects and changes in distance from the sensor. The existing GP model ETT algorithm adopts a GP model with a fixed external radius number, and is not suitable for a plurality of expansion targets with changing external shapes any more. Although existing methods can handle some weak variations in target size, no method has been proposed for the system to efficiently handle profile variations, particularly variations between ET and PT.
Disclosure of Invention
The invention aims to solve the problems, and provides a method for seamlessly tracking a point target and an extended target, so as to track the shape change of an ET (augmented reality) and simultaneously track a PT (potential Transformer) and the ET under the conditions of clutter and missing detection, and solve the problems of measurement and data association of shape points.
The key technology for realizing the invention is as follows: the variable dimension GP is used for modeling the variable outline of the target, the Poisson rate is introduced to estimate the number of measurement sources, and the number of the measurement sources reflects the number of outline points, so that the number of the radii of the GP model can be dynamically adjusted by the Poisson rate, and the model is adapted to the change of the outline of the target. And the PMHT algorithm is adopted to solve the data association problem, so that the tracking of a plurality of extended targets in the clutter is completed.
The technical problem proposed by the invention is solved as follows:
a method for seamless tracking of a point target and an extended target comprises the following steps:
step 1, initializing ET-GP-PMHT algorithm parameters:
step 1a, initializing state background parameters and GP model parameters: state transition matrix, state noise covariance, initial state covariance, hyper-parameters, etc.;
step 1b, initializing various parameters of the observation environment: observing noise variance R, clutter density, sampling interval Δ t, monitored space V, sensor position, detection probability Pd；
Step 1c, importing observation information: including T frame data, T in each sliding windowbFrame data, 1-T in sliding windowbFrame measurement data set Z, tth frame measurement data set ZtThe number of measurement sets M of the t-th framet，1≤t≤Tb；
Step 2. initialize T in sliding windowbSetting the current iteration times i as 1 for a frame data and measurement data set Z;
step 2a, removing the measurement related to the existing track in the measurement space, stacking the measurement with the Euclidean distance between the other measurements lower than the distance threshold in the measurement space, and initializing the new target track by using the two-point difference method of the measurement stack center of the previous 2 frames; if a new target track exists, initializing various parameters of the state environment: number of targets NtAnd initializing the number of GP model radiuses of the target track according to the result of comparing the number of elements in the measurement set obtained by each pile division in the first frame with the detection probability
Figure BDA0001860334790000021
Obtaining the target initial state by two-point difference
Figure BDA0001860334790000022
Initializing a state transition matrix, an initial state covariance, a state noise covariance, a poisson rate and a poisson parameter;
step 2b. introduction
Figure BDA0001860334790000023
Of a radius-corresponding outline point
Figure BDA0001860334790000024
Measuring the model;
step 3, constructing a posterior probability calculation formula of the T frame of the ET-GP-PMHT:
step 3a. poisson velocity vector:
Figure BDA0001860334790000025
wherein N is 1t，
Figure BDA0001860334790000026
The Poisson velocity vector λ is 1 to TbIn the poisson rate vector set of the frame, the invention assumes that the measured number of the target and the clutter number in the measurement interval meet poisson distribution, so the prior probability in the original PMHT can be replaced by the poisson rate, and the poisson rate can also reflect the average number of the measurement generated by the target; lambda [ alpha ]0，tThe number obeying mean value of the representative clutter is lambda0，tThe poisson distribution of (1), set as a constant in this invention; poisson rate λn，l，tObeys a distribution ofn，l，t|tAnd betan，l，t|tGamma distribution lambda as poisson parametern，l，t＝γ(λn，l，t；αn，l，t|t，βn，l，t|t)，αn，l，t|tAs a shape parameter, βn，l，t|tIs a scale parameter;
and 3b, calculating a likelihood according to the formula:
assuming that the clutter is spatially uniform, the likelihood value is:
Figure BDA0001860334790000031
wherein z isj，tIs the jth (j 1.., M) of the t framet) Measurement, xn，tIs the target state for the t frame n target,
Figure BDA0001860334790000032
is shown in
Figure BDA0001860334790000033
As a mean value, with Rn，l，tIs a gaussian probability density function of the covariance,
Figure BDA0001860334790000034
hl，t(. R) represents the metrology function of the metrology model corresponding to the ith topographical point of the target at time t, nn，l，tFor its covariance matrix corresponding to the metrology model,
Figure BDA0001860334790000035
the angle of the ith appearance point of the nth target of the tth frame on the global coordinate axis (as shown in FIG. 1) is the same, and the measurement models of different targets are the same;
step 3c, posterior probability formula:
Figure BDA0001860334790000036
wherein, ω isj，l，n，tMeasurement z at a time tj，tIs derived from the object xn，tThe posterior probability of the l contour point;
step 3d. poisson rate formula:
αn，l，t-1|t＝exp{-Δt/τ}αn，l，t|t βn，l，t-1|t＝exp{-Δt/τ}βn，l，t|t
Figure BDA0001860334790000037
βn，l，t|t＝βn，l，t|t-1+1
Figure BDA0001860334790000038
where exp is the exponential power, αn，l，t|t-1To predict the shape parameter, betan，l，t|t-1For predicting the scale parameter, tau is a time constant, which means the response speed of estimating the change of the evolution Poisson rate;
step 4, calculating the comprehensive measurement and the comprehensive covariance:
comprehensive measurement
Figure BDA0001860334790000039
And integrated covariance
Figure BDA00018603347900000310
Respectively as follows:
Figure BDA0001860334790000041
so far, the problem of fuzzy association between measurement and target appearance points in an extended target scene is solved, and only one comprehensive measurement and one comprehensive covariance exist for one appearance point of each target;
step 5, judging T as TbWhether the result is true or not, if so, executing the next step; otherwise, returning to execute the step 3, if t is t + 1;
step 6, extended Kalman smoothing:
because the measurement function of the extended target is nonlinear, the state is estimated by the extended Kalman smoothing algorithm. The measurement matrix, the comprehensive measurement and the comprehensive covariance can be stacked by adopting a stacking method, and then an extended Kalman smoothing algorithm is applied.
Because the measurement function is a nonlinear function, a Jacobian matrix is required to be obtained for the measurement function as a measurement matrix:
Figure BDA0001860334790000042
then, stacking the measurement matrix, the comprehensive measurement and the comprehensive covariance respectively to obtain:
Figure BDA0001860334790000043
Figure BDA0001860334790000044
Figure BDA0001860334790000045
wherein diag (·) represents a diagonalized matrix; finally, for the target xn，tThe algorithm steps of the executed extended Kalman smoothing algorithm are consistent with those of the traditional extended Kalman smoothing algorithm;
step 7, judging whether the iteration number i meets the loop iteration convergence condition, and returning to the step 3 if the iteration number i does not meet the loop iteration convergence condition; convergence is performed 8;
step 8, judging the track termination: defining an average estimated rate
Figure BDA0001860334790000046
If xi is smaller than threshold xiTHIf the flight path is finished, otherwise, the flight path continues;
step 9, self-adapting the dynamic target shape, and adjusting the number of target shape points:
step 9a, estimating the number of measurement sources by using the target Poisson rate:
Figure BDA0001860334790000047
estimating the difference between the number of the measurement sources and the number of the appearance points:
Figure BDA0001860334790000048
step 9b, if var is larger than 0, finding out the target radii with the front var poisson parameters being large, and adding new radii with the same poisson parameters beside the radii;
if var is less than 0, deleting the radius of the-var bar with the minimum Poisson parameter;
if var ≠ 0, update the state transition matrix, state noise covariance Qn，tState covariance Pn，t；
If var is 0, executing step 10;
step 10, judging whether the sliding window contains the last T of the T frame data setbFrame data, if not, the sliding window slides forward by TsAt one moment, a new in-window T is formedbReturning to execute the step 2 after the frame data and the measured data set Z are collected; otherwise the algorithm ends.
The method estimates the number of measurement sources by using the Poisson rate, adjusts the number of radii of the GP model, and records each sensor resolution unit occupied by the target contour as a measurement source. If the target observed relative to the sensor is smaller, the number of measurement sources on the target profile is smaller, the number of radii of the GP is estimated to be smaller, and vice versa. When the target is PT, ET-GP-PMHT tracks and outputs only the position of the target. The PMHT algorithm is adopted as the target tracking algorithm, and the algorithm is premised on the assumption that one target can generate a plurality of measurements and accords with the actual situation of an extended target.
The invention has the beneficial effects that:
(1) the ET-GP-PMHT can track the shape change of the ET and seamlessly track the mutual conversion of the ET and the PT, and when the shape is the ET with the unchanged form, the number of GP radiuses estimated by an algorithm is almost unchanged, so that the better tracking performance can be kept; when the target is PT, the ET-GP-PMHT tracks only the position of the target.
(2) The calculation complexity of the ET-GP-PMHT algorithm is related to the measurement number, the radius number and the target number. When the ET external deformation is small, the number of radii of the adopted GP model is reduced, and the calculation complexity is reduced.
Drawings
FIG. 1 is an ET of 16 outliers in the GP model, with two coordinates-local and global-the black points being 16 outliers;
FIG. 2 is a graph of estimated tracks and true targets for four targets in a single Monte Carlo simulation, where the true targets are represented by black lines, the estimated targets are represented by red lines, target trajectories and profiles are both shown, four plus signs represent the starting positions of the four target trajectories, and sensing is at the origin;
FIG. 3 is a RMSE of 100 Monte Carlo shots of four targets;
FIG. 4 is a legend to FIGS. 5, 6, 7, 8, where A1 represents ET-GP-PMHT, A2 represents ET-GP-PMHT-FBP26, A3 represents ET-GP-PMHT-FBP10, and A4 represents ET-RM-PMHT. ET-GP-PMHT-FBP26 is an ET-GP-PMHT-FBP algorithm with 26 appearance points;
FIG. 5 shows the real shape of the target and the estimated shapes of the four algorithms in the case of shape tracking at 121s, 181s, 241s and 691s, using ET-GP-PMHT-FBP26, ET-GP-PMHT-FBP10, ET-RM-PMHT and ET-GP-PMHT to track the target 1;
FIG. 6 is a graph of the profile tracking of target 2 at 95s, 125s, 335s and 635s using ET-GP-PMHT-FBP26, ET-GP-PMHT-FBP10, ET-RM-PMHT, ET-GP-PMHT;
FIG. 7 shows the case of tracking target 3 with ET-GP-PMHT-FBP26, ET-GP-PMHT-FBP10, ET-RM-PMHT, ET-GP-PMHT at 25s, 155s profile;
FIG. 8 shows the case of tracking the target 4 with ET-GP-PMHT-FBP26, ET-GP-PMHT-FBP10, ET-RM-PMHT, and ET-GP-PMHT at profiles of 31s and 91 s.
Detailed Description
The invention is further described below with reference to the figures and examples.
The embodiment provides a method for seamlessly tracking a point target and an extended target, wherein simulation is performed in a scene with clutter, as shown in fig. 1, four moving targets are tracked, and root mean square error RMSE (RMSE) verification algorithm performance is calculated. And comparing the ET-GP-PMHT with ET-GP-PMHT-FBP26, ETGP-PMHT-FBP10, and ET-RM-PMHT. The ET-GP-PMHT-FBP represents an ET-GP-PMHT algorithm with fixed number of appearance points, the ET-GP-PMHT-FBP26 represents that the number of the appearance points is 26, and the PMHT algorithm based on the RM model is marked as ET-RM-PMHT.
The state equation adopts a uniform linear model, the target 1 is close to the sensor from near to far and then close, the target is changed from an extended target to a point target and then to an extended target, and the motion time is 1-701 s; the target 2 is far from the sensor, the target is changed from a point target to an expanded target and then to a point target, and the movement time is 5-701 s. The objects 1 and 2 make uniform circular motion, i.e. turn, at 230-481s, and the angular velocity is 1/80 rad/s. The target 3 is far away from the sensor all the time and is represented as a point target, and the movement time is 21-230 s; the target 4 is always close to the sensor and is represented as an extended target with the same size, and the movement time is 31-230 s. Number of target real measurement sources
Figure BDA0001860334790000061
In relation to the position S between the target sensors:
Figure BDA0001860334790000062
the method of the embodiment comprises the following steps:
step 1, initializing ET-GP-PMHT algorithm parameters:
step 1a, initializing state background parameters and GP model parameters: state transition matrix, state noise covariance, initial state covariance, hyper-parameters, etc.;
the target state is
Figure BDA0001860334790000063
Figure BDA0001860334790000064
In the state of motion, the device is in motion,
Figure BDA0001860334790000065
in the form of the external shape state,
Figure BDA0001860334790000066
is formed by
Figure BDA0001860334790000067
The vector formed by the contour radii corresponding to the contour points,
Figure BDA0001860334790000068
Figure BDA0001860334790000069
representing the position of the target center in two-dimensional space,
Figure BDA0001860334790000071
indicating its corresponding speed, #n，tThe angle at which the target is rotated, i.e. the angle between the global and local coordinates (as shown in figure 1),
Figure BDA0001860334790000072
for the angular velocity, the angle and the angular velocity of the initialized target rotation are respectively 0rad and 0 rad/s; the state transition matrix is
Figure BDA0001860334790000073
Wherein the transition matrix of the motion state
Figure BDA0001860334790000074
A state transition matrix of a uniform velocity linear (CV) model and a shape state transition matrix
Figure BDA0001860334790000075
Figure BDA0001860334790000076
The dimension of expression is
Figure BDA0001860334790000077
Where the frame time interval Δ t is 1s, the forgetting parameter α of the state space＝0.0001；
State noise covariance
Figure BDA0001860334790000078
Noise of motion state
Figure BDA0001860334790000079
Is the state noise of CV model, in which the standard deviation of the state noise of position and angle is sigmaq＝0.05，σψ0.001, noise of shape state
Figure BDA00018603347900000710
Figure BDA00018603347900000711
Is a basic vector formed by the angles of the target outline points in the local coordinates, the angles of the outline points in the local coordinates are also called basic points (as shown in figure 1), the included angles between adjacent outline points are equal, namely, the outline points are distributed on the target outline in an angular average manner, and the outline points are distributed on the target outline in an equal manner
Figure BDA00018603347900000712
Covariance matrix for GP model:
Figure BDA00018603347900000713
covariance is a modified Squared Exponential (SE) function with a period of 2 pi, u and u' being arguments of the k function
Figure BDA00018603347900000714
The hyper-parameters of the GP model are set as: sigmar＝0.3，σfAnother hyper-parameter, the length scale, is adjusted according to the following rule:
Figure BDA00018603347900000715
the initial states of the four targets are respectively:
Figure BDA00018603347900000716
Figure BDA00018603347900000717
Figure BDA00018603347900000718
Figure BDA00018603347900000719
the state covariance of the four targets is
Figure BDA0001860334790000081
The motion state covariance is the diagonal matrix diag ([0.01,0.001,0.01,0.001,0.001, 0.0001)]) The covariance of the shape state is
Figure BDA0001860334790000082
Step 1b, initializing various parameters of the observation environment: observing noise variance R, clutter density, sampling interval Δ t, monitored space V, sensor position, detection probability Pd；
The sensor position is (0m,0m)TThe two-dimensional detection area x, y has a range of [0,450 ]]×[0，450]m2The clutter in the region is uniformly distributed, the number of the clutter in the region follows Poisson distribution, and the average clutter number at each moment is 20; the signal-to-noise ratio of the scene is 21dB, and the detection probability of the target is Pd＝0.92；
Measuring the noise as
Figure BDA0001860334790000083
exp { - Δ t/τ } -0.9; duration of each batch process T in PMHTbAt 3 frame times, a sliding length Ts2 frame times. The fixed cycle iteration times are 5 times in each batch of processing;
step 1c, importing observation information: including T frame data, T in each sliding windowbFrame data, 1-T in sliding windowbFrame measurement data set Z, tth frame measurement data set ZtThe number of measurement sets M of the t-th framet，1≤t≤Tb；
Step 2. initialize T in sliding windowbSetting the current iteration times i as 1 for a frame data and measurement data set Z;
step 2a, removing the measurement related to the existing track in the measurement space, dividing the measurement by a certain distance threshold value to judge whether the measurement is related to the existing track, stacking the measurement with the Euclidean distance between the rest measurements being lower than the distance threshold value by 15m in the measurement space, and initializing the new target track by using the center of the measurement stack of the previous 2 frames to carry out a two-point difference method; if a new target track exists, initializing various parameters of the state environment: number of targets NtAnd initializing the number of GP model radiuses of the target track according to the result of comparing the number of elements in the measurement set obtained by each pile division in the first frame with the detection probability
Figure BDA0001860334790000084
Obtaining the target initial state by two-point difference
Figure BDA0001860334790000085
Initial state covariance, state noise covariance, Poisson rate and Poisson parameter lambdan，l，t＝0.7，αn，l，t|t＝8，βn，l，t|t＝10；
Step 2b. introduction
Figure BDA0001860334790000086
Of a radius-corresponding outline point
Figure BDA0001860334790000087
Measuring the model;
after the ET-GP-PMHT algorithm environmental parameters are determined, an observation model is determined. Mapping from the state space to the observation space, and measuring the model of the outline point corresponding to the ith radius:
Figure BDA0001860334790000091
wherein, the direction vector is:
Figure BDA0001860334790000092
Figure BDA0001860334790000093
Figure BDA0001860334790000094
Figure BDA0001860334790000095
Figure BDA0001860334790000096
the angle of the ith contour point of the nth target on the local coordinate axis (as shown in FIG. 1) in the tth framet，jN (0, R) represents a Gaussian distribution with a mean of 0 and a covariance of R;
step 3, constructing a posterior probability calculation formula of the T frame of the ET-GP-PMHT:
step 3a. poisson velocity vector:
Figure BDA0001860334790000097
wherein N is 1t，
Figure BDA0001860334790000098
The Poisson velocity vector λ is 1 to TbIn the poisson rate vector set of the frame, the invention assumes that the measured number of the target and the clutter number in the measurement interval meet poisson distribution, so the prior probability in the original PMHT can be replaced by the poisson rate, and the poisson rate can also reflect the average number of the measurement generated by the target; lambda [ alpha ]0，tThe number obeying mean value of the representative clutter is lambda0，tThe poisson distribution of (1), set as a constant in this invention; poisson rate λn，l，tObeys a distribution ofn，l，t|tAnd betan，l，t|tGamma distribution lambda as poisson parametern，l，t＝γ(λn，l，t；αn，l，t|t，βn，l，t|t)，αn，l，t|tAs a shape parameter, βn，l，t|tIs a scale parameter;
and 3b, calculating a likelihood according to the formula:
assuming that the clutter is spatially uniform, the likelihood value is:
Figure BDA0001860334790000099
wherein z isj，tIs the jth (j 1.., M) of the t framet) Measurement, xn，tIs the target state for the t frame n target,
Figure BDA00018603347900000910
is shown in
Figure BDA00018603347900000911
As a mean value, with Rn，l，tIs a gaussian probability density function of the covariance,
Figure BDA00018603347900000912
hl，t(. R) represents the metrology function of the metrology model corresponding to the ith topographical point of the target at time t, nn，l，tFor its covariance matrix corresponding to the metrology model,
Figure BDA00018603347900000913
for the nth target of the t frame, the l outline point isAngles on the global coordinate axis (as in fig. 1), the measurement models of different targets are the same;
step 3c, posterior probability formula:
Figure BDA0001860334790000101
wherein, ω isj，l，n，tMeasurement z at a time tj，tIs derived from the object xn，tThe posterior probability of the l contour point;
step 3d. poisson rate formula:
αn，l，t-1|t＝exp{-Δt/τ}αn，l，t|t βn，l，t-1|t＝exp{-Δt/τ}βn，l，t|t
Figure BDA0001860334790000102
βn，l，t|t＝βn，l，t|t-1+1
Figure BDA0001860334790000103
where exp is the exponential power, αn，l，t|t-1To predict the shape parameter, betan，l，t|t-1For predicting the scale parameter, tau is a time constant, which means the response speed of estimating the change of the evolution Poisson rate;
step 4, calculating the comprehensive measurement and the comprehensive covariance:
comprehensive measurement
Figure BDA0001860334790000104
And integrated covariance
Figure BDA0001860334790000105
Respectively as follows:
Figure BDA0001860334790000106
so far, the problem of fuzzy association between measurement and target appearance points in an extended target scene is solved, and only one comprehensive measurement and one comprehensive covariance exist for one appearance point of each target;
step 5, judging T as TbWhether the result is true or not, if so, executing the next step; otherwise, returning to execute the step 3, if t is t + 1;
step 6, extended Kalman smoothing:
because the measurement function of the extended target is nonlinear, the state is estimated by the extended Kalman smoothing algorithm. Will slide the window inside TbThe 3s stacking method stacks the measurement matrix, the comprehensive measurement and the comprehensive covariance, and then uses the extended kalman smoothing algorithm.
Because the measurement function is a nonlinear function, a Jacobian matrix is required to be obtained for the measurement function as a measurement matrix:
Figure BDA0001860334790000107
then, stacking the measurement matrix, the comprehensive measurement and the comprehensive covariance respectively to obtain:
Figure BDA0001860334790000111
Figure BDA0001860334790000112
Figure BDA0001860334790000113
wherein diag (·) represents a diagonalized matrix; finally, for the target xn，tThe algorithm steps of the executed extended Kalman smoothing algorithm are consistent with those of the traditional extended Kalman smoothing algorithm;
step 7, judging whether the current iteration times i are equal to 5, and if not, returning to the step 3; if yes, executing 8;
step 8, judging the track termination: defining an average estimated rate
Figure BDA0001860334790000114
If xi is smaller than threshold xiTHIf the track is 0.2, the track is ended, otherwise, the track continues;
step 9, self-adapting the dynamic target shape, and adjusting the number of target shape points:
step 9a, estimating the number of measurement sources by using the target Poisson rate:
Figure BDA0001860334790000115
estimating the difference between the number of the measurement sources and the number of the appearance points:
Figure BDA0001860334790000116
step 9b, if var is larger than 0, finding out the target radii with the front var poisson parameters being large, and adding new radii with the same poisson parameters beside the radii;
if var is less than 0, deleting the radius of the-var bar with the minimum Poisson parameter;
if var ≠ 0, the transition matrix, state noise covariance Q is updatedn，tState covariance Pn，t；
If var is 0, executing step 10;
step 10, determining whether the sliding window contains a frame data set last T of 701sbFrame data, if not, the sliding window slides forward by TsForming new in-window T at 2s momentsbReturning to execute the step 2 after the frame data and the measured data set Z are collected; otherwise the algorithm ends.
In this example implementation, FIG. 2 shows a real target and an estimated target in a single Monte Carlo simulation, with the real and estimated target outlines drawn every 10 frames, and the point targets represented by five-pointed stars. FIG. 2 demonstrates that ET-GP-PMHT can initialize PT and ET, track PT and ET simultaneously, and seamlessly track the interconversion between PT and ET. When PT is converted into ET, the target measurement number is increased, and ET-GP-PMHT can continuously track the position of a target even the shape of the ET and can track the deflection of the target.
FIG. 3 shows the position RMSE of the four targets, with the peak in RMSE due to dynamic model mismatch when targets 1, 2 transition between ET and PT, and target 4 being a contour-invariant ET with RMSE less than PT target 3 because target 4 can detect more measurements.
For better performance of the checking algorithm, in addition to RMSE, we can also calculate the following performance index: the average track number of the target; average initialization time delay, namely the time difference between the start of tracking the target and the start of the real target; average track termination delay. As shown in the following table:
TABLE 1
Performance index Average number of tracks Average initialization time delay Average track end delay
Object
1 2.08 0s 0s
Object 2 1.49 1.62s 0.1s
Target 3 1.2 1.02s 10.88s
Target 4 1 0.02s 12.38s
ET-GP-PMHT was compared to ET-GP-PMHT-FBP26, ETGP-PMHT-FBP10, and ET-RM-PMHT. Target 1 has 10 measurement sources at 121s, as shown in fig. 5, ET-GP-PMHT and ET-GP-PMHT-FBP10 can better track the target profile, and the target profile estimated by ET-GP-PMHT-FBP26 is larger than the actual target profile; when the target 1 only has 2 measurement sources at 181s and 1 measurement source at 241s, the ET-GP-PMHT judges that the target is PT, and other algorithms still estimate the target shape of a closed curve; when the target measurement sources are more, such as 691s, ET-GP-PMHT and ET-GP-PMHT-FBP26 have better shape estimation.
Therefore, the ET-GP-PMHT-FBP is only suitable for tracking a target of a certain size, and the target shape estimated by the ET-RM-PMHT is elliptical regardless of the real ET shape even when the target is PT (see fig. 5, 6, 7, 8). The ET-GP-PMHT can simultaneously and stably track a large target, a small target and the PT, when the target shape in the graph 8 is not changed, the better tracking precision is kept, and the PT in the graph 7 can also be better tracked.
Finally, the above embodiments are only for illustrating the technical solutions of the present invention and are not limited, and all equivalent changes and modifications made in the claims of the present invention should be covered by the present invention.

Claims (5)

1. A method for seamless tracking of a point target and an extended target is characterized by comprising the following steps:
step 1, initializing ET-GP-PMHT algorithm parameters:
step 1a, initializing state background parameters and Gaussian process GP model parameters: state transition matrix, state noise covariance, initial state covariance, hyper-parameters, etc.;
step 1b, initializing various parameters of the observation environment: observing noise variance R, clutter density, sampling interval Δ t, monitored space V, sensor position, detection probability Pd；
Step 1c, importing observation information: including T frame data, T in each sliding windowbFrame data, 1-T in sliding windowbFrame measurement data set Z, tth frame measurement data set ZtThe number of measurement sets M of the t-th framet，1≤t≤Tb；
Step 2. initialize T in sliding windowbSetting the current iteration times i as 1 for a frame data and measurement data set Z;
step 2a, removing the measurement related to the existing track in the measurement space, stacking the measurement with the Euclidean distance between the other measurements lower than the distance threshold in the measurement space, and initializing the new target track by using the two-point difference method of the measurement stack center of the previous 2 frames; if a new target track exists, initializing various parameters of the state environment: number of targets NtAnd initializing the number of GP model radiuses of the target track according to the result of comparing the number of elements in the measurement set obtained by each pile division in the first frame with the detection probability
Figure FDA0003009393600000011
Obtaining the target initial state by two-point difference
Figure FDA0003009393600000012
Initializing a state transition matrix, an initial state covariance, a state noise covariance, a poisson rate and a poisson parameter;
step 2b. introduction
Figure FDA0003009393600000013
Of a radius-corresponding outline point
Figure FDA0003009393600000014
Measuring the model;
step 3, constructing a posterior probability calculation formula of the T frame of the ET-GP-PMHT:
step 3a. poisson velocity vector:
Figure FDA0003009393600000015
wherein N is 1t，
Figure FDA0003009393600000016
The Poisson velocity vector λ is 1 to TbThe method comprises the steps that a Poisson rate vector set of a frame replaces prior probability in original PMHT with Poisson rate, and the number of measured Poisson rate response targets is generated; lambda [ alpha ]0,tThe number obeying mean value of the representative clutter is lambda0,tPoisson distribution of (a); poisson rate λn,l,tObeys a distribution ofn,l,t|tAnd betan,l,t|tGamma distribution lambda as poisson parametern,l,t＝γ(λn,l,t；αn,l,t|t,βn,l,t|t)，αn,l,t|tAs a shape parameter, βn,l,t|tIs a scale parameter;
and 3b, calculating a likelihood according to the formula:
assuming that the clutter is spatially uniform, the likelihood value is:
Figure FDA0003009393600000021
wherein z isj,tIs the jth measurement of t frames, j 1t，xn,tIs the target state for the t frame n target,
Figure FDA0003009393600000022
is shown in
Figure FDA0003009393600000023
As a mean value, with Rn,l,tIs a gaussian probability density function of the covariance,
Figure FDA0003009393600000024
hl,t(. R) represents the metrology function of the metrology model corresponding to the ith topographical point of the target at time t, nn,l,tIs a covariance matrix corresponding to the metrology model,
Figure FDA0003009393600000025
the angle of the ith appearance point of the nth target of the t frame on the global coordinate axis is the same as the measurement models of different targets;
step 3c, posterior probability formula:
Figure FDA0003009393600000026
wherein, ω isj,l,n,tMeasurement z at a time tj,tIs derived from the object xn,tThe posterior probability of the l contour point;
step 3d. poisson rate formula:
αn,l,t-1|t＝exp{-Δt/τ}αn,l,t|t βn,l,t-1|t＝exp{-Δt/τ}βn,l,t|t
Figure FDA0003009393600000027
βn,l,t|t＝βn,l,t|t-1+1
Figure FDA0003009393600000028
where exp is the exponential power, αn,l,t|t-1To predict the shape parameter, betan,l,t|t-1τ is a time constant for predicting the scale parameter;
step 4, calculating the comprehensive measurement and the comprehensive covariance:
comprehensive measurement
Figure FDA0003009393600000029
And integrated covariance
Figure FDA00030093936000000210
Respectively as follows:
Figure FDA00030093936000000211
step 5, judging T as TbWhether the result is true or not, if so, executing the next step; otherwise, making t equal to t +1, and returning to execute the step 3;
step 6, extended Kalman smoothing:
and (3) solving a Jacobian matrix as a measurement matrix for the measurement function:
Figure FDA00030093936000000212
respectively stacking the measurement matrix, the comprehensive measurement and the comprehensive covariance to obtain:
Figure FDA0003009393600000031
Figure FDA0003009393600000032
Figure FDA0003009393600000033
wherein diag (·) represents a diagonalized matrix;
finally, for the target xn,tAn extended Kalman smoothing algorithm is executed;
step 7, judging whether the iteration number i meets the loop iteration convergence condition, and returning to the step 3 if the iteration number i does not meet the loop iteration convergence condition; convergence is performed 8;
step 8, judging the track termination: defining an average estimated rate
Figure FDA0003009393600000034
If x is less than threshold xTHIf the flight path is finished, otherwise, the flight path continues;
step 9, self-adapting the dynamic target shape, and adjusting the number of target shape points:
step 9a, estimating the number of measurement sources by using the target Poisson rate:
Figure FDA0003009393600000035
estimating the difference between the number of the measurement sources and the number of the appearance points:
Figure FDA00030093936000000314
step 9b, if var is greater than 0, finding out the radii of the previous var poisson parameters which are large, and adding new radii which are the same as the poisson parameters beside the radii;
if var is less than 0, deleting the radius of the-var bar with the minimum Poisson parameter;
if var10, the state transition matrix, state noise covariance Q, is updatedn,tState covariance Pn,t；
If var is 0, executing step 10;
step 10, judging whether the sliding window contains the last T of the T frame data setbFrame data, if not, the sliding window slides forward by TsAt one moment, a new in-window T is formedbReturning to execute the step 2 after the frame data and the measured data set Z are collected; otherwise the algorithm ends.
2. The method for tracking the point target and the extended target seamlessly according to claim 1, wherein the specific process of initializing the state background parameters and the GP model parameters in step 1a is as follows:
the target state is
Figure FDA0003009393600000036
Figure FDA0003009393600000037
In the state of motion, the device is in motion,
Figure FDA0003009393600000038
in the form of the external shape state,
Figure FDA0003009393600000039
is formed by
Figure FDA00030093936000000310
The radius values corresponding to the outline points form a vector,
Figure FDA00030093936000000311
is dynamically adjusted by the estimated number of measurement sources,
Figure FDA00030093936000000312
Figure FDA00030093936000000313
representing the position of the target center in two-dimensional space,
Figure FDA0003009393600000041
indicating its corresponding speed, #n,tWhich represents the angle of rotation of the target,
Figure FDA0003009393600000042
for the angular velocity, the angle and the angular velocity of the initialized target rotation are respectively 0rad and 0 rad/s; the state transition matrix is
Figure FDA0003009393600000043
Wherein the transition matrix of the motion state
Figure FDA0003009393600000044
A state transition matrix of a uniform linear model and a transition matrix of an appearance state
Figure FDA0003009393600000045
Figure FDA0003009393600000046
The dimension of expression is
Figure FDA0003009393600000047
The unit matrix of (1), wherein the frame time interval Δ t is 1s, and the forgetting parameter α of the state space is 0.0001;
state noise covariance
Figure FDA0003009393600000048
Noise of motion state
Figure FDA0003009393600000049
Is the state noise of CV model, in which the standard deviation of the state noise of position and angle is sigmaq＝0.05,σψ0.001, noise of shape state
Figure FDA00030093936000000410
Figure FDA00030093936000000411
Is a basic vector formed by the angles of the target outline points in local coordinates, and the included angles between the adjacent outline points are equal, i.e. the outline points are evenly distributed on the target outline in terms of angles
Figure FDA00030093936000000412
Covariance matrix for GP model:
Figure FDA00030093936000000413
the covariance is a modified square exponential function with a period of 2p, u and u' being arguments of the function k
Figure FDA00030093936000000414
The hyper-parameters of the GP model are set as: sigmar＝0.3,σfAnother hyper-parameter, the length scale q, is adjusted according to the following rule:
Figure FDA00030093936000000415
the initial states of the four targets are respectively:
Figure FDA00030093936000000416
Figure FDA00030093936000000417
Figure FDA00030093936000000418
Figure FDA00030093936000000419
the state covariance of the four targets is
Figure FDA00030093936000000420
Covariance of motion state
Figure FDA00030093936000000421
Is the diagonal matrix diag ([0.01,0.001,0.01,0.001,0.001, 0.0001)]) The covariance of the shape state is
Figure FDA0003009393600000051
3. The method for seamless tracking of the point target and the extended target according to claim 2, wherein the step 1b. initializing each parameter of the observation environment comprises:
the sensor position is (0m,0m)TThe two-dimensional detection area x, y has a range of [0,450 ]]×[0,450]m2The clutter in the region is uniformly distributed, the number of the clutter in the region follows Poisson distribution, and the average clutter number at each moment is 20; the signal-to-noise ratio of the scene is 21dB, and the detection probability of the target is Pd＝0.92；
Measuring the noise as
Figure FDA0003009393600000052
exp { - Δ t/τ } -0.9; duration of each batch process T in PMHTbAt 3 frame times, a sliding length Ts2 frame times; the number of iterations in each batch was 5 with a fixed loop.
4. The method for seamless tracking of a point target and an extended target according to claim 1, wherein the poisson rate and poisson parameters are: lambda [ alpha ]n,l,t＝0.7，αn,l,t|t＝8,βn,l,t|t＝10。
5. The method for seamless tracking of point target and extended target according to claim 3, wherein in step 2b
Figure FDA0003009393600000053
Of a radius-corresponding outline point
Figure FDA0003009393600000054
The measurement model is as follows:
mapping from the state space to the observation space, and measuring the model of the outline point corresponding to the ith radius:
Figure FDA0003009393600000055
wherein, the direction vector is:
Figure FDA0003009393600000056
Figure FDA0003009393600000057
Figure FDA0003009393600000058
Figure FDA0003009393600000059
Figure FDA00030093936000000510
is the angle of the ith target shape point of the t frame on the local coordinate axis, et,jN (0, R) represents a Gaussian distribution with a mean of 0 and a covariance of R.
CN201811338016.8A 2018-11-09 2018-11-09 A method for seamless tracking of point and extended targets Active CN109509207B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811338016.8A CN109509207B (en) 2018-11-09 2018-11-09 A method for seamless tracking of point and extended targets

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811338016.8A CN109509207B (en) 2018-11-09 2018-11-09 A method for seamless tracking of point and extended targets

Publications (2)

Publication Number Publication Date
CN109509207A CN109509207A (en) 2019-03-22
CN109509207B true CN109509207B (en) 2021-07-06

Family

ID=65747968

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811338016.8A Active CN109509207B (en) 2018-11-09 2018-11-09 A method for seamless tracking of point and extended targets

Country Status (1)

Country Link
CN (1) CN109509207B (en)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111077494A (en) * 2019-11-27 2020-04-28 江苏理工学院 A tracking method and device
DE102020121064A1 (en) 2020-08-11 2022-02-17 Valeo Schalter Und Sensoren Gmbh Method for operating an assistance system in a motor vehicle, computer program product, computer-readable storage medium and assistance system
CN113917504B (en) * 2021-09-30 2024-08-13 浙江工业大学 Position estimation method for controllable intelligent carrier
CN116736286B (en) * 2023-05-24 2024-02-06 兰州理工大学 Asymptotic Bayesian extended target tracking method and system based on stochastic hypersurface

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7212652B1 (en) * 2003-07-07 2007-05-01 The United States Of America As Represented By The Secretary Of The Navy Method for tracking targets with hyper-spectral data
CN101354786A (en) * 2007-07-23 2009-01-28 中国科学院计算技术研究所 A sports video event analysis method
CN103759732A (en) * 2014-01-14 2014-04-30 北京航空航天大学 Angle information assisted centralized multi-sensor multi-hypothesis tracking method
CN104237880A (en) * 2014-09-18 2014-12-24 中国人民解放军海军航空工程学院 Variable structure joint probability data interconnection formation target tracking method
US9851461B1 (en) * 2012-04-04 2017-12-26 The United States Of America As Represented By The Secretary Of The Navy Modular processing system for geoacoustic sensing
CN107526070A (en) * 2017-10-18 2017-12-29 中国航空无线电电子研究所 The multipath fusion multiple target tracking algorithm of sky-wave OTH radar
CN108363054A (en) * 2018-02-07 2018-08-03 电子科技大学 Passive radar multi-object tracking method for Single Frequency Network and multipath propagation

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7212652B1 (en) * 2003-07-07 2007-05-01 The United States Of America As Represented By The Secretary Of The Navy Method for tracking targets with hyper-spectral data
CN101354786A (en) * 2007-07-23 2009-01-28 中国科学院计算技术研究所 A sports video event analysis method
US9851461B1 (en) * 2012-04-04 2017-12-26 The United States Of America As Represented By The Secretary Of The Navy Modular processing system for geoacoustic sensing
CN103759732A (en) * 2014-01-14 2014-04-30 北京航空航天大学 Angle information assisted centralized multi-sensor multi-hypothesis tracking method
CN104237880A (en) * 2014-09-18 2014-12-24 中国人民解放军海军航空工程学院 Variable structure joint probability data interconnection formation target tracking method
CN107526070A (en) * 2017-10-18 2017-12-29 中国航空无线电电子研究所 The multipath fusion multiple target tracking algorithm of sky-wave OTH radar
CN108363054A (en) * 2018-02-07 2018-08-03 电子科技大学 Passive radar multi-object tracking method for Single Frequency Network and multipath propagation

Also Published As

Publication number Publication date
CN109509207A (en) 2019-03-22

Similar Documents

Publication Publication Date Title
CN109509207B (en) A method for seamless tracking of point and extended targets
CN103729859B (en) A kind of probability nearest neighbor domain multi-object tracking method based on fuzzy clustering
Bar-Shalom et al. Multisensor track-to-track association for tracks with dependent errors
CN104794735B (en) Extension method for tracking target based on variation Bayes's expectation maximization
CN105354860B (en) Extension target CBMeMBer trackings based on case particle filter
CN101719278B (en) Automatic tracking method for video microimage cells based on KHM algorithm
CN107633256A (en) Joint objective positioning and sensor registration method under a kind of multi-source ranging
CN114449452A (en) An Indoor Localization Algorithm for Heterogeneous Devices
CN107797106A (en) A kind of PHD multiple target tracking smooth filtering methods of the unknown clutter estimations of acceleration EM
CN111736144B (en) Maneuvering turning target state estimation method only by distance observation
CN116520281A (en) DDPG-based extended target tracking optimization method and device
CN111262556B (en) A multi-target tracking method for simultaneously estimating unknown Gaussian measurement noise statistics
CN113344970B (en) Irregular multi-extended target joint tracking and classifying method based on multiple Bernoulli
CN107391446A (en) Irregular shape based on random matrix extends target shape and method for estimating state more
CN104050686B (en) A kind of dense space method for tracking target
CN112836582B (en) Online detection method of structural change points of high-dimensional flow systems based on dynamic sparse subspace
CN111722213B (en) A pure distance extraction method for maneuvering target motion parameters
CN115877371A (en) Large-batch cluster target tracking method, system, equipment and medium based on serial GLMB
CN118759514A (en) A multi-target tracking method
CN108169722A (en) A kind of unknown disturbances influence the system deviation method for registering of lower sensor
CN115114985B (en) A Distributed Fusion Method for Sensor Systems Based on Set Theory
CN116338573B (en) Unmanned system positioning method taking noise error characteristics into consideration in closed space
CN111340853A (en) Multi-sensor GMPHD self-adaptive fusion method based on OSPA iteration
Huang et al. A bank of maximum a posteriori estimators for single-sensor range-only target tracking
CN115578425A (en) A Dynamic Tracking Method Based on Boundless Kalman Filter Applied in Fry Counter

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