Extended Target Tracking and Feature Estimation for Optical Sensors Based on the Gaussian Process

Abstract

A problem of tracking surface shape-shifting extended target by using gray scale pixels on optical image is considered. The measurement with amplitude information (AI) is available to the proposed method. The target is regarded as a convex hemispheric object, and the amplitude distribution of the extended target is represented by a solid radial function. The Gaussian process (GP) is applied and the covariance function of GP is modified to fit the convex hemispheric shape. The points to be estimated on the target surface are selected reasonably in the hemispheric space at the azimuth and pitch directions. Analytical representation of the estimated target extent is provided and the recursive process is implemented by the extended Kalman filter (EKF). In order to demonstrate the algorithm’s ability of tracking complex shaped targets, a trailing target characterized by two feature parameters is simulated and the two feature parameters are extracted with the estimated points. The simulations verify the validity of the proposed method with compared to classical algorithms.

1. Introduction

The extended target tracking (ETT) problem has been widely studied to estimate target motions and shape. With the development of the sensors, more information such as amplitude, color, etc. can be collected. This additional information should be fully utilized since it helps to improve the tracking precision and contains more abundant target characteristics. Therefore, new demand for an effective algorithm that can take advantage of additional information has emerged.

Many scholars have carried out research on ETT problems from different aspects such as shape modeling, measurement modeling and dynamic modeling. In this letter, we focus on the first aspect, i.e., the shape model method for an extended object. The current research achievements can be divided into the following three categories, corresponding to three levels of shape complexity. No extent: the simplest approach is to ignore the target extent. Under the assumptions that the measurements are distributed in the neighbor of target with a Poisson distributed number, Gaussian mixture probability hypothesis density (GM-PHD) filter is proposed [1,2]. Simple extent: a better modeling approach is to assume the target as a basic geometry, such as an ellipse, a circle or a rectangle. The most widely adopted method is based on a random matrix (RM), which represents the target extent as an ellipse [3,4,5,6,7]. For the target with non-elliptic shapes, a combination of multiple elliptical subobjects can be employed to approximate its real appearances [8,9,10]. Rectangle shape is usually adopted when tracking cars and other box-shaped objects [11,12,13]. Arbitrary extent: a more accurate and general approach is called random hypersurface models (RHMs) [14,15,16], which models the extent as a star-convex shape [17]. The RHMs can adapt to more flexible shapes in contrast with RM. The Gaussian process (GP) method can also be employed to estimate a target with highly nonlinear extent [18,19]. A more complete review for an ETT problem about shape model methods is presented in [20]. As mentioned in part Simple extent, the classical RM, RHMs and other variations only output the two-dimensional contour of the target, they do not have the ability to utilize extra information while the GP has the potential to realize that. The GP based method can theoretically estimate target of any dimension due to the fitting ability for highly nonlinear function.

GP has been adopted to learn the shape of a stationary three-dimensional object in [21,22], but none of them provide a recursive solution in a tracking application. A recursive solution for a cube-shaped object is provided in [23] by using point measurements collected from the target surface, but the target shape is relatively simple and there are no comparisons to the classical methods. In this letter, we extend the GP method to the problem of ETT on an image plane for optical sensors. Since the traditional RM, RHMs and GP are unable to utilize the amplitude information, we try to present a novel approach to overcome this problem. The surface of an extended target is described with solid radial function, rather than a two-dimensional counter. In order to adapt to the three-dimensional target tracking, the spherical difference was adopted rather than the angle difference of two two-dimensional angles. The RM and RHMs need to parameterize the target extent and estimate these parameters, while this letter estimates some selected points on the target surface. These selected points can represent the target distribution feature. Thus, the proposed method does not require any parametric model. The measurements with amplitude are converted into polar coordinates, and the measurements are used to estimate these points by the GP method with modified covariance function. In addition, we use the estimated points to extract the feature parameters of the target distribution, which are beneficial to the identification and classification of the target in the future work.

2. GP Theory

The GP is a kind of stochastic process in probability theory and mathematical statistics, which combines a series of random variables that follow normal distribution within an index set. Unknown complicated functions can be learned with training data. In this section, GP theory is directly applied to tracking application. More detail about the theory can be referred to [18].

2.1. Basic GP

GP is a typical stochastic process. it is always used to model spatially correlated measurements [24]. With t being input, the distribution of a GP $f(t)$ can be uniquely determined by its mean function $\mu (t)$ and covariance function $k(t,t^{\prime })$ as

$$\mu (t)=\mathrm{E}\left[f(t)\right],$$

(1)

$$k(t,t^{\prime })=\mathrm{E}\left[\left(f(t)-\mu (t)\right){\left(f(t^{\prime })-\mu (t^{\prime })\right)}^T\right].$$

(2)

Then, the GP can be marked as

$$f(t)\sim \mathcal{G}\mathcal{P}\left(\mu (t),k(t,t^{\prime })\right).$$

(3)

For the function value of the finite input $t_1\cdots t_N$, the Gaussian process is a generalization of the multivariate Gaussian probability distribution, i.e.,

$$\left[\begin{array}{c}f(t_1)\\ ⋮\\ f(t_N)\end{array}\right]\sim \mathcal{N}\left(\mu ,K\right),$$

(4)

where

$$\mu =\left[\begin{array}{c}\mu (t_1)\\ ⋮\\ \mu (t_N)\end{array}\right],K=\left[\begin{array}{ccc}k(t_1,t_1)& \cdots & k(t_1,t_N)\\ ⋮& & ⋮\\ k(t_N,t_1)& \cdots & k(t_N,t_N)\end{array}\right].$$

(5)

2.2. GP Regression

The distribution of a highly nonlinear function can be fitted with a Gaussian process regression (GPR) method. In this subsection, the GPR is adopted to model the target extent. In tracking application, the function value $\mathbf{f}={[f(t_1^f),\cdots ,f(t_N^f)]}^T$ for the corresponding input ${\mathbf{t}}^{\mathbf{f}}={[t_1^f,\cdots ,t_N^f]}^T$ is learned by using the measurement $\mathbf{z}={[z_1,\cdots z_M]}^T$ and the corresponding input $\mathbf{t}={[t_1,\cdots t_N]}^T$. For a target with amplitude information shown in Figure 1, the above measurement $\mathbf{z}$ is defined by $\{r_z,c_z,I_z\}$, where $r_z$, $c_z$ and $I_z$ represent the row, column, and amplitude of the measurement in grayscale image, respectively. In our application, the input $t=\psi$, where $\psi \stackrel{\Delta }{=}(\theta ,\phi )$ is the azimuth and pitch angle of a solid angle. The radius value $z_k$ and the function value $\mathrm{f}=f({\psi }_k)$ satisfy the joint Gaussian distribution

$$\left[\begin{array}{c}{z}_k\\ \mathbf{f}\end{array}\right]\sim \mathcal{N}\left(0,\left[\begin{array}{cc}k({\psi }_k,{\psi }_k)+R& K({\psi }_k,{\psi }^{\mathbf{f}})\\ K({\psi }^{\mathbf{f}},{\psi }_k)& K({\psi }^{\mathbf{f}},{\psi }^{\mathbf{f}})\end{array}\right]\right).$$

(6)

According to the GP theory, the likelihood and initial priori probabilities are given as follows:

$$p(\left.z_k\right|\mathbf{f})=\mathcal{N}(z_k;H_k^{\mathbf{f}}\mathbf{f},R_k^{\mathbf{f}}),$$

(7)

$$p(\mathbf{f})=\mathcal{N}(0,P_0^{\mathbf{f}}),$$

(8)

where

$$H^{\mathbf{f}}({\psi }_k)=K({\psi }_k,{\psi }^{\mathbf{f}}){\left[K({\psi }^{\mathbf{f}},{\psi }^{\mathbf{f}})\right]}^{-1},$$

(9)

$$R^{\mathbf{f}}({\psi }_k)=k({\psi }_k,{\psi }_k)+R-K({\psi }_k,{\psi }^{\mathbf{f}}){\left[K({\psi }^{\mathbf{f}},{\psi }^{\mathbf{f}})\right]}^{-1}K({\psi }^{\mathbf{f}},{\psi }_k),$$

(10)

$$\begin{array}{cc}{P}_0^{\mathbf{f}}=K({\psi }^{\mathbf{f}},{\psi }^{\mathbf{f}})& {\psi }^{\mathbf{f}}=({\theta }^{\mathbf{f}},{\phi }^{\mathbf{f}}).\end{array}$$

(11)

2.3. Recursive GPR

In order to realize the real-time recursive processing, the approximate method is usually used to solve the recursive problem. Bayes formula is applied to the posterior probability $p\left(\mathbf{f}\left|z_{1:N}\right.\right)$, which is

$$\begin{array}{cc}p\left(\mathbf{f}\left|z_{1:N}\right.\right)\hfill & \propto p\left(z_N\left|\mathbf{f},z_{1:N-1}\right.\right)·p\left(\mathbf{f}\left|z_{1:N-1}\right.\right),\hfill \\ & \propto \prod _{k=1}^{N}p\left(z_k\left|\mathbf{f},z_{1:k-1}\right.\right)·p\left(\mathbf{f}\right).\hfill \end{array}$$

(12)

Therefore, it can be approximated that $\mathbf{f}$ is independent from the previous measurement $z_{1:k-1}$, which means that $\mathbf{f}$ is completely statistically independent from the previous measurement, i.e.,

$$p\left(z_k\left|\mathbf{f},z_{1:k-1}\right.\right)=p\left(z_k\left|\mathbf{f}\right.\right).$$

(13)

2.4. Covariance Function Modification

In the covariance function used in Ref. [18], the squared Euler distance between the two plane angles is used to measure the correlation, i.e.,

$$k(\psi ,{\psi }^{\prime })={\sigma }_f^{2}exp\left(-\frac{d(\psi ,{\psi }^{\prime })}{2l^2}\right)+{\sigma }_r^2,$$

(14)

$$d(\psi ,{\psi }^{\prime })={∥\psi -{\psi }^{\prime }∥}^2,$$

(15)

where ${\sigma }_f^2$ denotes the variance of the measurement amplitude and l denotes the scale of the function; ${∥·∥}^2$ denotes the squared Euler distance. Since the target amplitude surface is considered as a convex hemisphere, it is obvious that the spherical distance is more reasonable to measure the correlation of the two solid angles, which is expressed as

$$d(\psi ,{\psi }^{\prime })=\arccos {[cos(\phi )cos({\phi }^{\prime })cos(\theta -{\theta }^{\prime })+sin(\phi )sin({\phi }^{\prime })]}^2.$$

(16)

3. EKF Model

Kalman filters such as EKF and Unscented Kalman filter (UKF) [25,26] are the most commonly used filters for nonlinear filtering. Since the EKF exerts lower computing complexity than the UKF, we employed EKF as the Bayesian approximation in this section.

3.1. Measurement Amplitude Modification

Since the measurement amplitude and target size are different in units as well as unequal in the values, in order to obtain a better extent estimation, we hope that the input $x^{\mathbf{f}}$ can be distributed as uniformly as possible on the target surface. Therefore, the measurement amplitude and target size should be close in order of magnitude, which means that the amplitude information needs to be reasonably modified. Let $S_k$ represent the number of pixels covered by the target, and let $I_{z_k}$ denote the amplitude value of all the target measurements collected from the optical image plane. Then, we define a modification factor $\beta$ as

$$\beta =\sqrt{S_k}/\phantom{\sqrt{S_k}max\left(I_{z_k}\right)}max\left(I_{z_k}\right).$$

(17)

With the amplitude value of the measurements multiplied by the modification factor $\beta$, the maximum value of the modified measurement $max(I_{z_k}^{\prime })=\beta ·max(I_{z_k})=\sqrt{S_k}$. It can be seen that, after modification, the amplitude value $I_{z_k}^{\prime }$ and the target scale $\sqrt{S_k}$ are at the same order of magnitude. Thus, the modified measurement is re-defined as

$${\tilde{z}}_k\stackrel{\Delta }{=}\{r_{z_k},c_{z_k},\beta I_{z_k}\},$$

(18)

where $r_{z_k}$, $c_{z_k}$ and $I_{z_k}$ represent the row, column, and amplitude of the measurement $z_k$. We now replace $z_k$ with the modified measurement ${\tilde{z}}_k$ for the rest of this letter.

3.2. Measurement Model with Amplitude

With the amplitude information considered, the target dimension is augmented to three, and the target extent is regarded as a hemispherical convex surface. Thus, the extent of the target can be represented by a function of the solid angle $\psi =(\theta ,\phi )$ and radius $r=f(\psi )$ in polar coordinates. Two coordinate systems are adopted. The sensor coordinate system is fixed to the sensor and local coordinate system is fixed to the body of the extended target, denoted by S and L, respectively. The origin of the local system ${\mathrm{o}}^L$ is combined with the mass center of the target and fixed on the image plane; the ${\mathrm{x}}^L$ axis lies in the image plane and points to the target moving direction, the ${\mathrm{z}}^L$ axis is perpendicular to the plane. ${\mathrm{x}}^L$, ${\mathrm{y}}^L$ and ${\mathrm{z}}^L$ follow the right-hand rule. The target extent is represented in the local coordinate system and the target center is presented in the sensor coordinate system.

As is shown in Figure 1, the target measurement ${\tilde{\mathbf{z}}}_{k,l}$ can be expressed in terms of the direction vector $\mathbf{p}({\mathbf{x}}_k^c)$ and radius $f({\psi }_{k,l}^L({\mathbf{x}}_{\mathrm{k}}^{\mathrm{c}},{\zeta }_{\mathrm{k}}))$ relative to the centroid ${\mathbf{x}}_k^c$, which is

$${\tilde{\mathbf{z}}}_{k,l}={\mathbf{x}}_k^c+\mathbf{p}({\mathbf{x}}_k^c)f({\psi }_{k,l}^L({\mathbf{x}}_k^c,{\zeta }_k))+{\overline{\mathbf{e}}}_{k,l},$$

(19)

where ${\overline{\mathbf{e}}}_{k,l}\sim \mathcal{N}\left(0,R\right)$, and $\mathbf{p}({\mathbf{x}}_k^c)=\mathbf{p}({\psi }_{k,l}^S({\mathbf{x}}_k^c))=\frac{{\tilde{\mathbf{z}}}_{k,l}-{\mathbf{x}}_k^c}{\left|\right|{\tilde{\mathbf{z}}}_{k,l}-{\mathbf{x}}_k^c\left|\right|}$ is the unit vector that points to the measurement ${\tilde{\mathbf{z}}}_{k,l}$ from centroid ${\mathbf{x}}_k^c$.

Reviewing the likelihood given in Equation (7), $H_k^{\mathbf{f}}\mathbf{f}$ is the mean of the random variable $z_k$ and $R_k^{\mathbf{f}}$ is the variance of it. Thus, the measurement $z_k$ can be expressed as

$$z_k=H_k^{\mathbf{f}}\mathbf{f}+e_k^{\mathbf{f}},e_k^{\mathbf{f}}\sim \mathcal{N}(0,R_k^{\mathbf{f}}).$$

(20)

In tracking application, the target extent ${\mathrm{x}}_k^{\mathbf{f}}$ is the instantiation of function value $\mathbf{f}$, i.e., $\mathbf{f}={\mathbf{x}}_k^{\mathbf{f}}$. Similar to Equation (20), the radius value on the ${\psi }_{k,l}^L({\mathbf{x}}_k^c,{\zeta }_k)$ direction can be expressed as

$$f({\psi }_{k,l}^L({\mathbf{x}}_k^c,{\zeta }_k))=H^{\mathbf{f}}({\psi }_{k,l}^L({\mathbf{x}}_k^c,{\zeta }_k)){\mathbf{x}}_k^{\mathbf{f}}+{\mathrm{e}}_{k,l}^{\mathbf{f}},$$

(21)

where $H^{\mathbf{f}}(·)$ denotes the observation model and $e_{k,l}^{\mathbf{f}}$ denotes the noise of it. Substitute Equation (21) into Equation (19), the measurement equation can be rewritten as

$$\begin{array}{cc}{\tilde{\mathbf{z}}}_{k,l}\hfill & ={\mathbf{x}}_k^c+{\mathbf{p}}_{k,l}({\mathbf{x}}_k^c)[H^{\mathbf{f}}({\psi }_{k,l}^L({\mathbf{x}}_k^c,{\zeta }_k)){\mathbf{x}}_k^{\mathbf{f}}+{\mathrm{e}}_{k,l}^{\mathbf{f}}]+{\overline{\mathbf{e}}}_{k,l}\hfill \\ & =\underset{\tilde{\mathbf{h}}({\mathbf{x}}_k^c,{\zeta }_k,{\mathbf{x}}_k^{\mathbf{f}})}{\underbrace{{\mathbf{x}}_k^c+{\tilde{H}}_l({\mathbf{x}}_k^c,{\zeta }_k){\mathbf{x}}_k^{\mathbf{f}}}}+\underset{{\tilde{\mathbf{e}}}_{k,l}}{\underbrace{{\mathbf{p}}_{k,l}({\mathbf{x}}_k^c){\mathbf{e}}_{k,l}^{\mathbf{f}}+{\overline{\mathbf{e}}}_{k,l}}}\hfill \\ & ={\tilde{\mathbf{h}}}_{k,l}({\mathbf{x}}_k)+{\tilde{\mathbf{e}}}_{k,l},\hfill \end{array}$$

(22)

where

$${\tilde{H}}_l({\mathbf{x}}_k^c,{\psi }_k)={\mathbf{p}}_{k,l}({\mathbf{x}}_k^c)H^{\mathbf{f}}({\psi }_{k,l}^L({\mathbf{x}}_k^c,{\zeta }_k)),$$

(23)

$$H^{\mathbf{f}}({\psi }_{k,l}^L({\mathbf{x}}_k^c,{\zeta }_k))=K({\psi }_{k,l}^L,{\psi }^{\mathbf{f}}){[K({\psi }^{\mathbf{f}},{\psi }^{\mathbf{f}})]}^{-1},$$

(24)

$${\mathbf{x}}_k\stackrel{\Delta }{=}({\mathbf{x}}_k^c,{\zeta }_k,{\mathbf{x}}_k^{\mathbf{f}}).$$

(25)

According to Equation (22), it is easy to obtain the new noise covariance as follows:

$${\tilde{R}}_{k,l}={\mathbf{p}}_{k,l}{R}_{k,l}^{\mathbf{f}}{\mathbf{p}}_{k,l}^T+{\overline{R}}_{k,l},$$

(26)

$$R_{k,l}^{\mathbf{f}}=R^{\mathbf{f}}({\psi }_{k,l}^L({\mathbf{x}}_k^c,{\zeta }_k)),$$

(27)

where $R_{k,l}^{\mathbf{f}}$ is the covariance of noise ${\mathbf{e}}_{k,l}^{\mathbf{f}}$, and ${\overline{R}}_{k,l}$ is the covariance of noise ${\overline{\mathbf{e}}}_{k,l}$.

3.3. Motion Model

We model the target state as ${\mathbf{x}}_k={({\overline{\mathbf{x}}}_k,{\mathbf{x}}_k^{\mathbf{f}})}^T$, where ${\overline{\mathbf{x}}}_k=({\mathbf{x}}_k^c,{\zeta }_k)$ is the centroid, ${\zeta }_k$ denotes the orientation of the extended target and ${\mathbf{x}}_k^{\mathbf{f}}$ is the extent. Because the target centroid is different from that of the extended, the centroid and the extent state are modeled separately. The dynamics of the target centroid ${\overline{\mathbf{x}}}_k$ is consistent with the constant velocity model:

$${\overline{\mathbf{x}}}_{k+1}=\overline{F}{\overline{\mathbf{x}}}_k+{\overline{\omega }}_k{\overline{\omega }}_k\sim \mathcal{N}(0,\overline{Q}),$$

(28)

where the transition matrix $\overline{F}$ and covariance matrix $\overline{Q}$ for target centroid ${\overline{\mathrm{x}}}_k$ are given as

$$\overline{F}=\left[\begin{array}{cc}1& T\\ 0& 1\end{array}\right]\otimes I_3,\overline{Q}=\left[\begin{array}{cc}{T}_s^3/\phantom{T_s^{3}3}3& T_s^2/\phantom{T_s^{2}2}2\\ T_s^2/\phantom{T_s^{2}2}2& T_s\end{array}\right]\otimes \left[\begin{array}{ccc}{\sigma }_q^2& 0& 0\\ 0& {\sigma }_q^2& 0\\ 0& 0& {\sigma }_{q^{\zeta }}^2\end{array}\right],$$

(29)

where T denotes the sampling period, ⊗ denotes the Kronecker product, $I_3$ denotes the third order identity matrix, ${\sigma }_q^2$ and ${\sigma }_{q^{\zeta }}^2$ represents the process variance of the centroid and orientation. The dynamics of target extent ${\mathbf{x}}_k^{\mathbf{f}}$ is shown as follows:

$${\mathbf{x}}_{k+1}^{\mathbf{f}}=F^{\mathbf{f}}{\mathbf{x}}_k^{\mathbf{f}}+{{\omega }^{\mathbf{f}}}_k{{\omega }^{\mathbf{f}}}_k\sim \mathcal{N}(0,Q^{\mathbf{f}}),$$

(30)

where the extent transition matrix and covariance matrix are presented as

$$F^{\mathbf{f}}=exp(-\alpha T),$$

(31)

$$Q^{\mathrm{f}}=(1-exp(-2\alpha T))·K({\psi }^{\mathrm{f}},{\psi }^{\mathrm{f}}),$$

(32)

where $\alpha$ is the forgetting factor.

Finally, transition model of the augmented state, i.e., the combination of dynamics of target state ${\overline{\mathbf{x}}}_k$ and extent ${\mathbf{x}}_k^{\mathbf{f}}$, is given as

$${\mathbf{x}}_{k+1}=F{\mathbf{x}}_k+{\omega }_k{\omega }_k\sim \mathcal{N}(0,Q),$$

(33)

where

$$\begin{array}{ccc}{\mathbf{x}}_k={[{\overline{\mathbf{x}}}_k,{\mathbf{x}}_k^{\mathbf{f}}]}^T& F=\left[\begin{array}{cc}\overline{F}& 0\\ 0& F^{\mathbf{f}}\end{array}\right]& Q=\left[\begin{array}{cc}\overline{Q}& 0\\ 0& Q^{\mathbf{f}}\end{array}\right].\end{array}$$

(34)

The initial state ${\mathbf{x}}_0$ is given as

$$\begin{array}{cc}{\mathbf{x}}_0\left({\mu }_0,P_0\right),{\mu }_0={[{\overline{\mu }}_0,0]}^T& P_0=\left[\begin{array}{cc}{\overline{P}}_0& 0\\ 0& P_0^{\mathbf{f}}\end{array}\right],\end{array}$$

(35)

where ${\overline{\mu }}_0$ and ${\overline{P}}_0$ denote the mean and covariance of ${\overline{\mathbf{x}}}_k$, $P_0^{\mathbf{f}}$ denotes the covariance of ${\mathbf{x}}_k^{\mathbf{f}}$, which can be calculated by Equation (22).

So far, the state space description is provided in Equations (22) and (33). Now, the estimation can be performed using a nonlinear filter and EKF is employed in this letter. The calculation of Jacobian matrix is deeply discussed in the appendix of Ref [18].

4. Simulation and Analysis

4.1. Simulation Setting

An extended target with grayscale information moves in the $[300\times 300]$ surveillance area. The total number of scans is set to 100 s and the scan cycle is set to $T_s=1$ s. The grayscale amplitude along the velocity direction of the target follows the type I extreme value distribution with the scale parameter ${\sigma }_{ex}$. The grayscale amplitude on a direction perpendicular to velocity obeys the Gaussian distribution with the standard deviation ${\sigma }_{gau}$. A view of the amplitude distribution model of the target is shown in Figure 2.

The standard deviations of the process noise on position and angle are set to be ${\sigma }_q=1$ and ${\sigma }_{q^{\zeta }}=0.1$, respectively, forgetting factor $\alpha =0.0001$. The standard deviations of the measurement noise on position and amplitude are set as ${\sigma }_r=1$ and ${\sigma }_I=0.1$. The hyper-parameter of the GP are set to ${\sigma }_f=2$, ${\sigma }_r=2$ and $l=\pi /\phantom{\pi 4}4$. Those pixels with amplitude value greater than $T{h}_{gray}=0.01$ are retained as target measurements z and the others are removed. Each simulation are tested within 50 Monte Carlos.

4.2. Feature Estimation Steps

In addition to the estimation of the target centroid ${\mathbf{x}}_k^c$, another important goal is to extract the target feature parameters based on the estimated values ${\mathbf{x}}_k^{\mathbf{f}}$ at all directions. In this study, the steps shown in Figure 3 are designed to extract the target distribution feature parameters ${\sigma }_{ex}$ and ${\sigma }_{gau}$. In this feature estimation method, the estimated points along velocity direction and vertical velocity direction are extracted. These points are fitted with function $f_1$ and $f_2$, which is given below. After the fitting step, the center estimation $({\widehat{x}}_c^L,{\widehat{y}}_c^L)$ and feature parameters ${\sigma }_{ex}$ and ${\sigma }_{gau}$ are easily attained.

For the estimation step, too many ${\psi }^{\mathbf{f}}$ will increase the calculation cost and cause the irreversibility of the covariance function $K({\psi }^{\mathbf{f}},{\psi }^{\mathbf{f}})$. However, if the number of input is small, the estimation accuracy of target extent cannot be guaranteed. As a compromise, we set the input number as follows. For the proposed GP with AI, eight points are uniformly sampled within the interval of $\left[0\sim 2\pi \right]$ as azimuth ${\theta }^{\mathbf{f}}$ and five points are uniformly sampled within the interval of $\left[0\sim \pi /2\right]$ as ${\phi }^{\mathbf{f}}$. For the traditional GP, eight points are uniformly sampled within the interval of $\left[0\sim 2\pi \right]$. Since the number of Fourier coefficients in RHMs has to be odd, and to ensure that the number of parameters to be evaluated is consistent with that of GP as much as possible, we set the coefficients number to 9.

As for the extraction and fitting step, the feature parameters ${\sigma }_{ex}$ and ${\sigma }_{gau}$ can be subsequently obtained by fitting the estimated points on velocity direction and vertical velocity direction. Since the distribution of the trailing target on velocity and vertical velocity direction follows the type I extreme value distribution and Gaussian distribution, the functions $f_1$ and $f_2$ are given as

$$f_1=A_1·\underset{type\mathrm{I}extremevaluedistribution}{\underbrace{\frac{1}{{\sigma }_{ex}}·exp\left(\frac{x-{\mu }_{ex}}{{\sigma }_{ex}}\right)·exp\left(-exp\left(\frac{x-{\mu }_{ex}}{{\sigma }_{ex}}\right)\right)}},$$

(36)

$$f_2=A_2·\underset{Gaussiandistribution}{\underbrace{\frac{1}{\sqrt{2\pi }{\sigma }_{gau}}·exp\left(\frac{{\left(x-{\mu }_{gau}\right)}^2}{{\sigma }_{gau}^2}\right)}}.$$

(37)

Subsequently, ${\widehat{x}}_c^L$ and ${\widehat{y}}_c^L$, corresponding to the maximum value on each fitting curve in Figure 1, denotes the target local coordinate position, respectively.

4.3. Shape Deformation

The target moves at a constant speed and the initial state of the target centroid is ${[31,23,1.2,1.1]}^{\mathrm{T}}$. The scale parameter ${\sigma }_{ex}=3.5$ and becomes 2.5 at the 50th second. The standard deviation ${\sigma }_{gau}=1.5$ and becomes 2 at the 50th second. The sum of the target grayscale is 200. The performance is measured with the tracking precision and target feature parameter estimation accuracy. Figure 4 shows a snapshot of the tracking process. For display convenience, we added a bias of 20 pixels on the y-axis of the estimation. It can be observed that the shape estimate changes with the measurement progressively.

Figure 5 compares the tracking precision of the RHMs, GP, and the proposed GP with AI. It can be seen from Figure 5 that the proposed method outperforms the RHMs and GP method. Because the RHMs and the traditional GP method can only use the position information of the measurements, the contribution of each pixel to the target centroid is the same, so traditional methods can only estimate the center of form. On the other hand, the proposed GP with AI can estimate the center of mass with the amplitude information introduced. The estimation of target feature ${\sigma }_{ex}$ and ${\sigma }_{gau}$ are given in Figure 6, the red dotted lines denote the ground truth of target features and the blue lines denote the estimation of target features. It can be seen that the proposed method estimates the ground truth quite well during the first 50 s. Even if the target extent mutates in the 50th second, the proposed method is able to approximate the ground truth gradually with the new measurements’ accumulation. The result demonstrates that the proposed method has the ability to adapt to the shape-shifting extended target. There exists deviation between the estimation and the ground truth value since the systematic error is brought in by the truncation of measurement amplitude by $T{h}_{gray}$ and the fitting step.

4.4. Motion Change

The scale parameter ${\sigma }_{ex}=3.5$ and standard deviation ${\sigma }_{gau}=1.5$. The sum of target grayscale is 200. The target motion model follows Equation (38) during the 20th to 70th second and constant speed model during the rest of the simulation time

$$F=\left[\begin{array}{cccc}1& 0& \sin (\omega T_s)/\omega & -(1-\cos (\omega T_s))/\omega \\ 0& 1& (1-\cos (\omega T_s))/\omega & \sin (\omega T_s)/\omega \\ 0& 0& \cos (\omega T_s)& -\sin (\omega T_s)\\ 0& 0& \sin (\omega T_s)& \cos (\omega T_s)\end{array}\right],$$

(38)

where $\omega$ is the turn rate. Figure 7 shows a snapshot of a target with turn rate $\omega =\pi /80$. Figure 8 shows the feature parameter estimation results of the targets moving with $\omega =\pi /80$, $\pi /100$ and $\pi /140$, respectively.

It can be seen that the feature parameter estimation performance goes worse with the increase of the target maneuverability.

4.5. Intensity Change

The scale parameter ${\sigma }_{ex}=3.5$ and standard deviation ${\sigma }_{gau}=1.5$. The sum of target grayscale changes from 200 to 500 at the 50th second. Figure 9 shows a snapshot of intensity change and Figure 10 shows the feature estimation results.

It can be seen that, even if the target parameters’ estimation has a large deviation as the intensity abrupt changes at the 50th second, it tends to converge to the ground truth as long as the intensity of the target holds steady.

5. Conclusions

In this letter, a GP based method for the tracking of an extended target by using amplitude information is proposed. Surface shape and kinematic state are simultaneously estimated. The superiority of utilizing the amplitude information is demonstrated compared with traditional RHMs and GP methods. As a result, the tracking error implemented by the proposed method reduces about 0.5 pixels compared with the RHMs and GP. Furthermore, by using the estimated points on the target surface, we also extract the target feature parameters, which can be used for target classification in further study. The estimation error of feature parameters ${\sigma }_{ex}$ and ${\sigma }_{gau}$ are able to converge to 0.17 and 0.1, respectively. Through challenging scenes in Section 4.4 and Section 4.5, the estimation performance can be limited, influenced by the mutations on target maneuvering and intensity changing, which proves the robustness of our proposed method. In fact, the tracking performance improvement is mainly due to the introduction of amplitude information. If the amplitude information is ignored, this method degrades to the traditional GP method, and its performance will be no different from that of the GP and RHM methods. Since the proposed method does not require any parametric model, it is flexible enough to estimate any convex hemisphere-shaped target.