Space Target Tracking with the HRRP Characteristic-Aided Filter via Space-Based Radar

Abstract

Approaching space target tracking is a typical and challenging mission in the space situational awareness (SSA) field. As the space-based radar is able to monitor the space targets of interest full-weather all-time, the space-based radar system is utilized in this paper. However, most multi-target tracking (MTT) filters in target tracking studies merely utilize the location or narrow measurements, and many potentially valuable electromagnetic scattering characteristics are missed, which leads to space target false tracking problems. The space-based radar transmits a wide-band signal, and the measured high-resolution range profile (HRRP) information is an effective characteristic for different target discrimination. Therefore, the HRRP characteristics of space targets are implemented into the update recursion of the MTT filter, which can be utilized to improve the tracking performance. Then, to predict the target HRRP sequence, the geometrical theory of diffraction (GTD) model is utilized. Additionally, a modified spatial spectrum method with a novel covariance matrix is designed to improve the scattering parameter estimation accuracy. Finally, an adapting threshold is devised for merging the Gaussian mixture (GM) components weights. The proposed threshold is on the basis of the proposed HRRP characteristic-aided probability hypothesis density (PHD) filter, and it can tackle the problem of space target discrimination. Simulation results validate the effectiveness and robustness of the proposed probability hypothesis density (HGI-PHD) filter aided by HRRP information and improved with GM weights.

1. Introduction

1.1. Background and Problem Statement

The radar system in this paper is the space-based radar, and the reasons for utilizing a space-based radar to detect or track space targets are analyzed as follows. Note that with the rapid growth of stealth technology and space activity, the existing ground-based equipment fails to monitor our space targets of interest all-weather full-time. Compared with ground-based equipment and space-based optical equipment, a space-based radar can provide better coverage, and it is not constrained by the weather or the curvature of the Earth. Moreover, it possesses a stronger early warning ability for high-speed space targets. Therefore, the space-based radar is utilized to track space targets in this paper, and some valuable problems are solved under space-based observing scenarios.

The problems that exist in the above-mentioned space-based observing scenarios are essential for perceiving the orbital state of the space targets of concern in real-time accurately, which is a challenging mission for space situational awareness (SSA) [1,2]. In general, the SSA surveillance tasks conclude the effective observation of our space targets of interest in the surveillance region, i.e., how the space target orbits change with varying time in order to estimate the individual characteristics or states (such as their positions and velocities in the reference coordinates). It should be mentioned that the number of space targets will change as time varies, and it is not known as a priori knowledge; therefore, the SSA tasks can be regarded as classical multiple target tracking (MTT) problems.

However, MTT will encounter some specific and typical problems when faced with the SSA surveillance tasks. For example, the space-based radar can only cover a limited tracking range due to its relatively smaller radar transmitting power compared to the ground-based radar system. Furthermore, it is challenging to track space targets accurately when they approach each other or the distance between them is close. What is more, all target state components sometimes are not provided thoroughly in a space-based radar system, and these characteristics make it more difficult and complicated for MTT problems.

An appealing and effective way to solve the aforementioned problems is using the random finite set (RFS) theory and finite set statistics (FISST). It is known that FISST is an elegant Bayesian formulation for the description of MTT based on the RFS theory [3,4]. Efficient solutions such as the probability hypothesis density (PHD) [5,6] and cardinalized PHD (CPHD) [7,8,9] have emerged in the MTT region in recent years. The core goal for the aforementioned filters is to update the first-order multi-target moments of the posterior probability density and update the global posterior probability density rather than recurrence, which is able to reduce the computational complexity to a large degree. Moreover, with the aim of alleviating poor cardinality estimation performance accuracy of the standard PHD filter, the CPHD filter is proposed to ease the above tension. It should be noted that the sequential Monte Carlo (SMC) model [10] and the Gaussian mixtures (GM) model [11] are the two most typical implementations of the above-mentioned Bayesian-based filters.

Note that among the various MTT filters, the PHD filter and the CPHD filter have concise formulation expressions and a slightly small computational burden, which is feasible for engineering applications. However, due to the “spooky action” in the process of PHD filter recursion, some limits exist as we propagate these two filters’ legacy PHD. We would like to mention that it is difficult to derive a theoretical mathematical expression for the posterior probability density. Therefore, to tackle this problem, the multi-target cardinality-balanced multi-Bernoulli (CBMeMer) filter [12] is developed by propagating the posterior PDF. But under high-clutter density scenarios, the tracking performance of the CBMeMer filter will deteriorate sharply. In recent years, the generalized labeled multi-Bernoulli (GLMB) filter [13] and its computational simplified form, the labeled multi-Bernoulli (LMB) filter [14], have been proposed sequentially. These two filters are able to distinguish various target trajectories more accurately. However, considering the moment approximate operation, it is hard for the LMB filter to effectively avoid missed detection.

The space-based radar system is able to transmit a wide-band signal, and, thus, it can provide targets additional electromagnetic information like the high-resolution range profile (HRRP). Motivated by the above, we designed a filter that can update the target HRRP pseudo-likelihood. In order to predict the HRRP sequence and employ it in update recursion, the geometrical theory of diffraction (GTD) model is utilized in this paper. Being a classical scattering center model, the GTD model is effective in describing the radar targets' electromagnetic characteristics at high frequencies. In recent years, the GTD scattering center model has been widely applied in lots of radar domains, such as automatic target recognition (ATR), radar cross-section (RCS) extrapolation and interpolation, radar target three-dimensional (3D) reconstruction, and so on. Once the GTD scattering center is reconstructed, we can obtain the targets’ electromagnetic scattering data, and the HRRP sequence can be predicted as well. Therefore, it is vitally important to construct a precise GTD scattering model and estimate the scattering parameters accurately. Thus, many effective methods, such as estimating signal parameters via the rational invariant technique (ESPRIT) algorithm [15,16,17,18], the multiple signal classification (MUSIC) algorithm [19,20,21,22], and the matrix enhancement and matrix pencil (MEMP) algorithm [23,24,25], have been proposed to estimate the scattering parameters from the back-scattered electromagnetic data.

Based on the aforementioned background as well as to deeply analyze the scientific problems in the process of tracking space targets based on a space-based radar, we separately analyzed the specific problems in the observation scene of a space-based radar and provided specific solutions so as to effectively improve the tracking accuracy of space targets.

1.2. Problem Analysis and Contributions

It is worth stressing that some electromagnetic scattering characteristics remain in radar signals, while most MTT filters fail to utilize them properly. Therefore, some drawbacks existed during the tracking recursion as follows:

• Generally, most MTT filters merely utilize target information such as distance and orientation but fail to use some targets’ additional characteristics that a wide-band radar can provide, for example, the HRRP information.
• It is difficult to track space targets accurately when they are closely spaced with each other. So how can we improve the tracking performance when space targets are in dense clutter environments (due to thermal noise in space-based platforms, space debris, or other space targets)?
• If the HRRP characteristics are utilized in the MTT tracking recursion process, then how can we predict the HRRP sequence, and how can we compute the HRRP pseudo-likelihood for different targets?
• A more accurate GTD model is important for providing precise electromagnetic information for tracking update procedures. However, due to the thermal noise in a space-based platform, the signal-to-noise ratio (SNR) will be lower than the ground-based radar, which will lead to inaccuracy in scattering center parameter estimation. Therefore, it is vitally important to design a novel scattering parameter estimation method that can perform well in low-SNR scenarios.

The standard PHD filters do not utilize the targets’ potential electromagnetic information abundantly, and they will encounter the problem of poor tracking performance under low-SNR scenarios or dense clutter environments. To improve the aforementioned problems, a modified PHD filter with HRRP characteristics is proposed in this paper, and the main contributions are shown as follows.

• To fully utilize the electromagnetic scattering characteristics of space targets, HRRP information is utilized in the update recursion of our proposed MTT filter. The HRRP for various targets and clutter differs significantly from each other as tracking time varies; thus, if the HRRP information is fully utilized, the difference between different targets or clutter will be enlarged, which is beneficial for improving the tracking accuracy.
• To solve the aforementioned second problem, an adaptive threshold for merging the Gaussian components is designed in this paper. The designed adaptive threshold is relevant to the noise covariance, and it is able to improve the tracking performance of the MTT filter. Furthermore, the discrimination ability between closely spaced targets can be enhanced at the same time.
• In order to predict the HRRP sequence, the GTD scattering center model is constructed in this paper. Then, the similarity rate between different HRRP sequences is denoted to describe the HRRP pseudo-likelihood.
• To solve the fourth problem mentioned above with the aim of improving the GTD scattering parameter estimation accuracy in low-SNR scenarios, a modified 3D-ESPRIT algorithm with a novel correlation matrix is proposed in this paper. Furthermore, by squaring the total covariance matrix, the parameter estimation accuracy can be further improved.

The rest of this paper is organized as follows. Section 2 investigates a review of the PHD filter and the GTD scattering center model. Furthermore, to obtain the time-varying observing angle, the observation geometry transformation for a space-based radar is also given in Section 2. Section 3 presents our novel scattering parameter estimation method and the analytic implementation of our proposed PHD (HGI-PHD) filter aided by HRRP characteristics and improved with Gaussian mixture components. The simulation and performance evaluation are illustrated in Section 4. Finally, conclusions are drawn in Section 5.

2. Technical Foundation

In this section, we first introduce the basic theory of PHD filtering, which provides a foundation for space target tracking. Then, to predict the HRRP sequence of space targets, the GTD scattering center model is demonstrated in this section. Finally, the observation geometry transformation for the space-based radar is derived to compute the time-varying observing angle.

2.1. PHD Filtering

In order to ease the computational intractability of the Bayesian multiple-target filter, the PHD filter is proposed by using the first-order moment of an RFS, which can be called the intensity function, and it is used to approximate the posterior multiple-target state density and propagate posterior intensity. The intensity of an RFS is defined as follows:

$$v_k^{}\left(x|Z_{1:k}^{}\right)={\displaystyle \int {\delta }_X(x)f\left(X|Z_{1:k}^{}\right)\delta X}$$

(1)

where ${\delta }_X(x)={\displaystyle \sum _{w\in X}{\delta }_w(x)}$ and ${\delta }_w(x)$ represent the Dirac delta function. $X$ and $Z$ represent the states and measurements, respectively. Here, $Z_{1:k}^{}$ is omitted from the intensity to notate it more simply. The integral of the PHD is not a constant, and the expected number of targets is expressed as

$${\displaystyle \int v_k^{}(x)\mathrm{d}x={\lambda }_k^{}}$$

(2)

The PHD filter approximates the density of multiple targets by using a Poisson process as follows:

$$f_{k|k}^{}(X)=e^{-{\lambda }_k^{}}\prod _{x\in X}{\lambda }_k^{}{v}_k^{}(x)$$

(3)

Then, the PHD prediction recursion is computed as

$$v_{k|k-1}^{}(x)={\displaystyle \int p_{S,k}(\zeta )f(x|\zeta )v_{k-1}^{}(\zeta )\mathrm{d}\zeta +{\gamma }_k(x)}$$

(4)

where $v_{k|k-1}^{}(x)$ represents the predicted PHD density of time $k$, $p_{S,k}(\zeta )$ is the target survival probability, ${\gamma }_k(x)$ is the target birth intensity function, and $f(x|\zeta )$ denotes the Markov transition density with state $x$.

Finally, the PHD update step can be summarized as

$$v_k(x)=\left[1-p_{D,k}(x)\right]v_{k|k-1}^{}(x)+{\displaystyle \sum _{z\in Z}\frac{p_{D,k}(x)g_z(z|x)v_{k|k-1}^{}(x)}{\kappa (z)+{\displaystyle \int p_{D,k}(\xi )g_z(z|\xi )v_{k|k-1}^{}(\xi )\mathrm{d}\xi }}}$$

(5)

where at time $k$, $p_{D,k}(x)$ is the detection probability of the target, $g_z(z|x)$ is the measurement likelihood of the target, and $\kappa (z)$ is the clutter RFS intensity function.

2.2. GTD Scattering Center Model

As one of the typical scattering center models, the geometric theory of diffraction (GTD) model enables us to describe the electromagnetic characteristics of radar targets effectively. Therefore, in a high-frequency region, the GTD scattering center model can be expressed as [26]

$$\begin{array}{l}E(f_m,{\theta }_n,{\phi }_k)={\displaystyle \sum _{i=1}^I{A}_i}{(\mathrm{j}\frac{f_m}{f})}^{{\alpha }_i}\exp [-4\pi \mathrm{j}{f}_m(x_i\cos \theta \cos \phi +y_i\sin \theta \cos \phi +z_i\sin \phi )/\mathrm{c}]+\omega (f_m,{\theta }_n,{\phi }_k)\\ ={\displaystyle \sum _{i=1}^I{A}_i}{(\mathrm{j}\frac{f_0+m\Delta f}{f})}^{{\alpha }_i}\exp [-4\pi \mathrm{j}{f}_m(x_i\cos \theta \cos \phi +y_i\sin \theta \cos \phi +z_i\sin \phi )/\mathrm{c}]+\omega (f_m,{\theta }_n,{\phi }_k)\end{array}$$

(6)

where $E(f_m,{\theta }_n,{\phi }_k)$ represents the back-scattering data of radar targets, I denotes the total scattering centers, and $\{A_i,{\alpha }_i,x_i,y_i,z_i\}$ denote the scattering intensity, scattering type, transversal position parameter, and longitudinal position parameter and vertical position parameter of the i-th scattering center, respectively. $f_m=f_0+m\Delta f$ represents the operating frequency, and $\{f_0,\Delta f,m\}$ represent the initial frequency, the frequency step, and the frequency index, respectively. In similarity, ${\theta }_n$ is equal to ${\theta }_0+n\Delta \theta$; ${\phi }_k$ is equal to ${\phi }_0+k\Delta \phi$; ${\theta }_0$ and ${\phi }_0$ are the initial azimuth angle and the initial elevation angle, respectively; and $n\Delta \theta$ and $k\Delta \phi$ are the radar rotation angles, which are relatively small. $\mathrm{c}=3\times {10}^8 \mathrm{m}/\mathrm{s}$ represents the electromagnetic wave propagation speed, and $\omega (f_m,{\theta }_n,{\phi }_k)$ denotes the Gaussian white noise. The type scattering parameter ${\alpha }_i$ of typical scattering structures is given in [26].

Then, according to Ref. [27], the scattering center parameters ${\alpha }_i$, $x_i$, $y_i$, and $z_i$ can be estimated as

$${\alpha }_i=\frac{(\left|P_{xi}\right|-1)f_0}{\Delta f}$$

(7)

$$x_i=\frac{-\mathrm{angle}(P_{xi})\times \mathrm{c}}{4\pi \Delta f_x}$$

(8)

$$y_i=\frac{-\mathrm{angle}(P_{yi})\times \mathrm{c}}{4\pi \Delta f_y}$$

(9)

$$z_i=\frac{-\mathrm{angle}(P_{zi})\times \mathrm{c}}{4\pi \Delta f_z}$$

(10)

where $P_{xi}=(1+{\alpha }_i\frac{\Delta f_x}{f_{x0}})\exp (-4\pi \mathrm{j}\Delta f_x{x}_i/c)$, $P_{yi}=\exp (-4\pi \mathrm{j}\Delta f_y{y}_i/c)$, and $P_{zi}=\exp (-4\pi \mathrm{j}\Delta f_z{z}_i/c)$.

Finally, the intensity parameters $\tilde{A}$ can be estimated by using the least square method as follows

$$\tilde{A}={(G^{\mathrm{H}}G)}^{-1}{G}^{\mathrm{H}}{E}_{\mathrm{k}}$$

(11)

where

$$G=[a_1,\dots ,a_I]$$

(12)

$$a_i={\left[\begin{array}{c}{a}_i(0,0,0),\dots ,a_i(M-1,0,0),a_i(0,1,0),\dots ,\\ a_i(M-1,1,0),\dots ,a_i(M-1,N-1,0),\\ a_i(0,0,1),\dots ,a_i(M-1,N-1,K-1)\end{array}\right]}^{\mathrm{T}}$$

(13)

$$a_i(m,n,k)=\mathrm{j}{(\frac{f_m}{f_0})}^{{\alpha }_i}\exp [\frac{-4\pi \mathrm{j}{f}_m}{\mathrm{c}}(x_i\cos {\theta }_n\cos {\phi }_k+y_i\sin {\theta }_n\cos {\phi }_k+z_i\sin {\phi }_k)]$$

(14)

$$E_k=\left[\begin{array}{c}E(f_0,{\theta }_0,{\phi }_0),\dots ,E(f_{M-1},{\theta }_0,{\phi }_0),a_i(f_0,{\theta }_1,{\phi }_0),\dots ,\\ a_i(f_{M-1},{\theta }_1,{\phi }_0),\dots ,a_i(f_{M-1},{\theta }_{N-1},{\phi }_0),\\ a_i(f_0,{\theta }_0,{\phi }_1),\dots ,a_i(f_{M-1},{\theta }_{N-1},{\phi }_{K-1})\end{array}\right]$$

(15)

${\left[\right]}^{\mathrm{T}}$ represents the transpose operation.

2.3. Observation Geometry Transformation for Space-Based Radar

In this section, the observation geometry transformation for a space-based radar is derived to obtain accurate observing angles, which are significant for GTD scattering center estimation and the HRRP prediction process.

As depicted in Figure 1, the Earth-centered, Earth-fixed (ECEF) coordinate system $O_e-X_e{Y}_e{Z}_e$ and the radial tangential normal (RTN) coordinate system $O_t-X_t{Y}_t{Z}_t$ are marked in black lines and in red lines, respectively. The origin of the ECEF coordinate system is at the center of the earth, and it is a non-rotating right-handed Cartesian coordinate system. The direction of its $X_e$ axis points to the equator, and the prime meridian intersection and the $Z_e$ axis pass through the Earth's north pole. As for the RTN coordinate system, the $X_t$ axis is parallel to the direction of the velocity, the $Z_t$ axis is along the main body of the space target, and the $Y_t$ axis is determined according to the right-hand theorem.

The orbit of the space target is determined by six orbital elements, which are derived from a priori two-line elements (TLEs). To be specific, these elements are expressed as $h$ (semi-major axis), $\alpha$ (inclination), $\Omega$ (right ascension of the ascending node), $e$ (eccentricity), $\omega$ (argument of perigee), and $\theta$ (true anomaly). Thereafter, the observation geometry transformation between the ECEF coordinate system and the RTN coordinate system can be derived according to the given six orbital elements.

Then, the positions of the space target and the space-based radar in the ECEF coordinate system are given as

$$\begin{array}{c}{O}_{t_1}{}_{ECEF}=R_z({\Omega }_1)R_x({\alpha }_1)R_z({\omega }_1)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}·{\left[\begin{array}{ccc}{h}_1\frac{e_1+\cos {\theta }_1}{1+e_1\cos {\theta }_1}& \left(h_1\frac{1-e_1{}^2}{1+e_1\cos {\theta }_1}\right)\sin {\theta }_1& 0\end{array}\hspace{1em}\right]}^{\mathrm{T}}\end{array}$$

(16)

$$\begin{array}{c}{O}_{t_2}{}_{ECEF}=R_z({\Omega }_2)R_x({\alpha }_2)R_z({\omega }_2)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}·{\left[\begin{array}{ccc}{h}_2\frac{e_2+\cos {\theta }_2}{1+e_2\cos {\theta }_2}& \left(h_2\frac{1-e_2{}^2}{1+e_2\cos {\theta }_2}\right)\sin {\theta }_2& 0\end{array}\hspace{1em}\right]}^{\mathrm{T}}\end{array}$$

(17)

where $R_x$, $R_y$, and $R_z$ represent the rotation matrix around the $x$ axis, the $y$ axis, and the $z$ axis, respectively.

Also, the rotation matrices can be expressed as follows:

$$R_z(\Omega )=\left[\begin{array}{ccc}\cos \Omega & -\sin \Omega & 0\\ \sin \Omega & \cos \Omega & 0\\ 0& 0& 1\end{array}\right]$$

(18)

$$R_x(\alpha )=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \alpha & -\sin \alpha \\ 0& \sin \alpha & \cos \alpha \end{array}\right]$$

(19)

$$R_z(\omega )=\left[\begin{array}{ccc}\cos \omega & -\sin \omega & 0\\ \sin \omega & \cos \omega & 0\\ 0& 0& 1\end{array}\right]$$

(20)

Thus, the coordinate transformation matrix between the ECEF coordinate system and the RTN coordinate system can be expressed as

$$M_{ECEF-OF}={\left[R_z(\Omega )R_x(\alpha )R_z(\omega +\theta )\right]}^{-1}$$

(21)

So, the positions of the space target and space-based radar in the RTN coordinate system are given as

$$O_{t_1}{}_{OF}=M_{ECEF-OF}\left({\Omega }_1,{\alpha }_1,{\omega }_1,{\theta }_1\right)O_{t_1}{}_{ECEF}$$

(22)

$$O_{t_2}{}_{OF}=M_{ECEF-OF}\left({\Omega }_1,{\alpha }_1,{\omega }_1,{\theta }_1\right)O_{t_2}{}_{ECEF}$$

(23)

Also, the line of sight (LOS) of the radar in the RTN reference coordinate system is computed as

$$LOS=O_{t_1}{}_{OF}-O_{t_2}{}_{OF}=\left(x_{\Delta },y_{\Delta },z_{\Delta }\right)$$

(24)

$${\theta }_{\Delta }=\arctan \left(\frac{x_{\Delta }}{y_{\Delta }}\right)$$

(25)

$${\phi }_{\Delta }=\arctan \left(\frac{z_{\Delta }}{\sqrt{x_{\Delta }^2+y_{\Delta }^2+z_{\Delta }^2}}\right)$$

(26)

where $x_{\Delta },y_{\Delta },z_{\Delta }$ represent the location vector of the LOS in the RTN coordinate system and ${\theta }_{\Delta }$ and ${\phi }_{\Delta }$ are the observing azimuth angle and elevation angle in the RTN coordinate system, respectively.

By deriving the observation geometry transformation for the space-based radar, the inherent motion of space targets can be eliminated from consideration in this paper.

3. Proposed Method

In Section 3.1, to predict the HRRP sequence more accurately, we first introduce a novel scattering parameter estimation method based on an improved 3D-ESPRIT algorithm. Then, in Section 3.2, HRRP characteristics obtained via the wide-band radar and HRRP pseudo-likelihood predicted by using the GTD scattering center model are employed in the PHD filter to fully utilize space targets' electromagnetic characteristics. Finally, the main steps of our proposed HGI-PHD filter are given in Section 3.3. To be specific, the flowchart of our proposed HGI-PHD filter is shown in Figure 2.

3.1. Parameter Estimation of Scattering Centers

In this section, in order to predict the HRRP sequence more accurately, an improved 3D-ESPRIT algorithm with reconstructed covariance matrices is proposed to estimate the 3D-GTD model parameters, and the main contribution of our proposed 3D-ESPRIT algorithm is shown as follows.

Firstly, by performing the spatial smoothing operation in the $x$ direction, a Hankel matrix pencil $X^x$ can be given by reference [27]

$$X^x=\left[\begin{array}{cccc}{X}_0^x& X_1^x& \cdots & X_{M-P}^x\\ X_1^x& X_2^x& \cdots & X_{M-P+1}^x\\ ⋮& ⋮& \cdots & ⋮\\ X_{P-1}^x& X_P^x& \cdots & X_{M-1}^x\end{array}\right]$$

(27)

where

$$X_m^x=\left[\begin{array}{cccc}x(m,0)& x(m,1)& \cdots & x(m,K-L)\\ x(m,1)& x(m,2)& \cdots & x(m,K-L+1)\\ ⋮& ⋮& \cdots & ⋮\\ x(m,L-1)& x(m,L)& \cdots & x(m,K-1)\end{array}\right]$$

(28)

$$x(m,k)=\left[\begin{array}{cccc}E(m,0,k)& E(m,1,k)& \cdots & E(m,N-Q,k)\\ E(m,1,k)& E(m,2,k)& \cdots & E(m,N-Q+1,k)\\ ⋮& ⋮& \cdots & ⋮\\ E(m,Q-1,k)& E(m,Q,k)& \cdots & E(m,N-1,k)\end{array}\right]$$

(29)

where $P\in [I+1,M-I+1]$, $Q\in [I+1,N-I+1]$, $L\in [I+1,K-I+1]$, and $M$, $N$, $K$, and $I$ represent the frequency steps, azimuth angle steps, elevation angle steps, and the scattering center numbers, respectively.

Then, a permutation matrix $J$ is defined as

$$J={\left[\begin{array}{cccc}0& \dots & 0& 1\\ ⋮& 0& 1& 0\\ 0& ⋰& \ddots & ⋮\\ 1& 0& \dots & 0\end{array}\right]}_{PQL\times PQL}$$

(30)

By combining the permutation matrix $J$ and the original back-scattering electromagnetic data $X^x$, a new matrix $E_{conj}$ containing the covariance information of the original scattered data $X^x$ is obtained by

$$E_{conj}=J·X^x$$

(31)

Then, we construct the following three covariance matrices to fully utilize the electromagnetic scattering characteristic of radar targets:

$$\{\begin{array}{c}{R}_{X_x{X}_x}=X_x{X}_x{}^{\mathrm{H}}\\ R_{E_{conj}{E}_{conj}}=E_{conj}{E}_{conj}{}^{\mathrm{H}}\\ R_{X_x{E}_{conj}}=X_x{E}_{conj}{}^{\mathrm{H}}\end{array}$$

(32)

Next, by averaging the above three covariance matrices, a novel covariance matrix $R$ is proposed in this paper, which is shown as follows:

$$R=\frac{R_{X_x{X}_x}+R_{E_{conj}{E}_{conj}}+R_{X_x{E}_{conj}}}{3}$$

(33)

Note that $R$ is a Hermitten matrix, and it satisfies $R=R^{\mathrm{H}}$. Thus, by squaring matrix $R$, we have the final covariance matrix $R_1$ as follows:

$$R_1=R{R}^{\mathrm{H}}=R^2$$

(34)

Also, the following relation satisfies

$$\{\begin{array}{c}{\lambda }_1{}_i={\lambda }_i{}^2\\ {\Lambda }_1{}_i={\Lambda }_i\end{array}\begin{array}{cc}& i=1,\dots ,I\end{array}$$

(35)

where ${\lambda }_1{}_i$ and ${\Lambda }_1{}_i$ are the eigenvalue and the eigenvector of $R_1$, respectively, and ${\lambda }_i$ and ${\Lambda }_i$ are the eigenvalue and the eigenvector of $R$, respectively.

Note that from (34), the eigenvalues of $R_1$ are twice those of $R$. Therefore, the differences between the noise eigenvalues and the signal eigenvalues can be broadened by constructing $R_1$.

To analyze our motivation for constructing the aforementioned covariance matrix mathematically, the estimated parameters variance is derived as follows:

$$E\left\{{(\hat{\xi }-\xi )}^2\right\}=\frac{\frac{{\sigma }^2}{2MNK}{\displaystyle \sum _{i=1}^I}-\frac{{\gamma }_i}{{\left({\sigma }^2-{\gamma }_i\right)}^2}{\left|{\mathit{G}}^{\mathrm{H}}·{\left(v_i\right)}^{\mathrm{H}}\right|}^2}{{\displaystyle \sum _{i=I+1}^{MNK}}-{\left|\left[d{\mathit{G}}^{\mathrm{H}}/d(z)\right]·{\left(v_i\right)}^{\mathrm{H}}\right|}^2}$$

(36)

where $E\left\{\right(\widehat{\xi }-\xi {)}^2\}$ represents variance of the estimated parameters; $\widehat{\xi }$ and $\xi$ represent the estimated parameter and the original parameter, respectively; ${\sigma }^2$ and ${\gamma }_i$ denote eigenvalues of noises and eigenvalues of signals, respectively; $v_i={\gamma }_{i}I-X^x$ represents the eigenmatrix of ${\gamma }_i$; and $I$ denotes a identify matrix.

It can be seen from (36) that $E\left\{\right(\widehat{\xi }-\xi {)}^2\}$ decreases when ${\sigma }^2$ and ${\gamma }_m$ differ greatly from each other, which will bring out more accurate scattering center parameters. So, constructing $R$ can broaden the differences between ${\sigma }^2$ and ${\gamma }_m$, which can estimate the GTD model parameters more accurately.

Then, the scattering parameters $\{{\alpha }_i,x_i,y_i,z_i\}$ are estimated according to reference [27], and the intensity parameters $\tilde{A}$ can be obtained by using the least square method in (27). It should be noticed that for the sake of readability, the detailed CRB derivation of the GTD scattering center parameter is given in Appendix A.

3.2. HRRP Characteristic-Aided Filter

Once we have estimated the scattering center model parameters, the HRRP sequence can be predicted at a relatively small observing angle. Thus, the predicted HRRP information can be utilized in the update recursion of the PHD filter. In this section, we will introduce the specific implementation of our proposed HRRP characteristic-aided filter.

Let us denote the surveillance region as three-dimensional; then, the target state can be expressed as

$$x_k={\left[l_{x,k},v_{x,k},l_{y,k},v_{y,k},l_{z,k},v_{z,k}\right]}^{\mathrm{T}}$$

(37)

where at time $k$, $l_{x,k}$, $l_{y,k}$, and $l_{z,k}$ are the three-dimensional locations of targets. $v_{x,k}$, $v_{y,k}$, and $v_{z,k}$ are the three-dimensional velocities of targets.

Thereafter, as the HRRP characteristics are employed, the measurement state of targets at time $k$ yields

$$z_k=\left[z_{l,k},z_{hrrp,k}\right]$$

(38)

where at time $k$, $z_{l,k}$ and $z_{hrrp,k}$ denote the location measurements and the HRRP measurements, respectively.

The specific prediction and update recursion of the proposed HGI-PHD filter are discussed as follows.

Prediction:

The prediction recursion of the proposed filter is written as

$$D_{k|k-1}(x_k)={\gamma }_k(x_k)+{\displaystyle \int \left(p_{s,k}(x_{k-1})f_{k|k-1}(x_k|x_{k-1})+b_{k|k-1}(x_k|x_{k-1})D_{k-1|k-1}(x_{k-1})\right)}d{x}_{k-1}$$

(39)

where at time $k$, ${\gamma }_k(x_k)$ denotes the intensity of the birth RFS, $p_{s,k}(x_{k-1})$ is the survival probability, $f_{k|k-1}(x_k|x_{k-1})$ is the Markov transmission density, $b_{k|k-1}(x_k|x_{k-1})$ represents the density of the spawned at time $k$ under the condition of the previous state at time $k-1$, and $D_{k-1|k-1}(x_{k-1})$ is the density of the PHD filter at time $k-1$.

Update:

After employing the HRRP characteristics in the standard PHD filter, the update recursion of the proposed HGI-PHD filter can be given by

$$D_{k|k}(x_k)=L_Z(x_k)·D_{k|k-1}(x_k)$$

(40)

Assuming the HRRP characteristic of targets is independent of location states, the likelihood for targets and for clutter can be denoted as

$$g(\mathit{z}|x)=g((z_{l,k},z_{hrrp,k})|x)=g_l(z_l|x)g_{hrrp}(hrrp)$$

(41)

$$c(\mathit{z}|x))=c_z(z|x)c_{hrrp}(hrrp)$$

(42)

where $g_l(z_l|x)$ and $g_{hrrp}(hrrp)$ are the state location likelihood and HRRP pseudo-likelihood functions for targets, respectively. $c_z(z|x)$ and $c_{hrrp}(hrrp)$ are the state location likelihood and HRRP pseudo-likelihood functions for clutter, respectively.

However, how do we predict the HRRP sequence as time varies, and how do we compute the HRRP pseudo-likelihood for various targets and clutter? These two problems are core points to be solved in this paper, and we give the detailed process as follows.

(1)

HRRP sequence prediction

Note that the radar is able to provide range and bearing measurements of the target’s center, and the observation model of kinematic measurement is denoted as follows:

$$M_k=\left[R_k,{\beta }_k\right]=h\left(l_k\right)+w_k$$

(43)

where at time $k$, $h(·)$ denotes the kinematic observation function; $R_k$ and ${\beta }_k$ represent the target range measurements and bearing measurements, respectively; and $w_k$ represents the zero-mean observation noise matrix.

For the space target, due to its previously known orbital information and low maneuverability, the velocity direction is almost aligned with the axial direction of the target's main body, which is shown in Figure 3.

Therefore, the observing angle can be computed as

$${\phi }_k={\theta }_k-{\beta }_k={\cos }^{-1}\left(\frac{x_k{v}_{x,k}+y_k{v}_{y,k}}{\sqrt{x_k^2+y_k^2}\sqrt{v_{x,k}^2+v_{y,k}^2}}\right)$$

(44)

where ${\theta }_k={\tan }^{-1}\left(\frac{v_{y,k}}{v_{x,k}}\right)$ denotes the heading angle for space targets.

Once the specific observing angles are given as time varies, different back-scattering echoes can be obtained via the 3D scattering center model in (6). Then, the frequency response of the m-th frequency point can be written as $E_m=E(f_m,{\theta }_n,{\phi }_k)$, and by taking an inverse discrete Fourier transform (IDFT) on frequency response sequence $\mathit{E}=E_0,E_1,\dots E_M$, the HRRP sequence of space targets can be predicted eventually, which is denoted as $hrr{p}_{k\_est}$ in this paper.

(2)

HRRP pseudo-likelihood evaluation

At time $k$, we can acquire the actual HRRP sequence of different space targets or clutter, which is denoted as $hrr{p}_{k\_act}$. Then, the HRRP pseudo-likelihood can be computed as

$$g_{hrrp,k}=\frac{\langle hrr{p}_{k\_act},hrr{p}_{k\_est}\rangle }{\Vert hrr{p}_{k\_act}\Vert ·\Vert hrr{p}_{k\_est}\Vert }$$

(45)

where $\langle a,b\rangle ={\displaystyle \int a(x)}b(x)dx$ represents the inner product operation and $\Vert h\Vert$ is the $l_2$ norm function.

Note that the posterior intensity at time $k-1$ is in the form of a Gaussian mixture; then, we have

$$v_{k|k-1}(x)={\displaystyle \sum _{i=1}^{J_{k|k-1}}{w}_{k|k-1}^i\mathcal{N}\left(x;m_{k-1}^i,P_{k-1}^i\right)}$$

(46)

where $w_{k|k-1}^i$ is the i-th Gaussian mixture weight at time $k-1$; $m_{k-1}^i$ and $P_{k-1}^i$ denote the mean vector and covariance vector of the i-th Gaussian mixture weight at time $k-1$, respectively; and $J_{k|k-1}$ represents the number of Gaussian mixture weights at time $k-1$.

Substituting (46) into (40), we have the posterior intensity at time $k$ as

$$D_{k|k}(x_k)=\left[1-p_d\right]D_{k|k-1}(x_k)+{\displaystyle \sum _{i=1}^{J_{k|k-1}}{\displaystyle \sum _{z_k\in {\rm Z}}{\tilde{w}}_{k|k-1}^i\mathcal{N}\left(x_k;m_{k|k-1}^i,P_{k|k-1}^i\right)}}$$

(47)

where the normalized GM weight is denoted as

$${\tilde{w}}_{k|k-1}^j(z_k)=\frac{p_{D,k}(x)g_z(z|x)g_{hrrp}(hrrp)w_{k|k-1}^j{q}_k^j(z_k)}{\kappa (z)c_{hrrp}(hrrp)+p_{D,k}(x)g_z(z|x)g_{hrrp}(hrrp){\displaystyle \sum _{i=1}^{J_{k|k-1}}{w}_{k|k-1}^i{q}_k^i(z_k)}}$$

(48)

where $\{\begin{array}{c}{q}_k^i(z_k)=N\left(z_k,{\eta }_{k|k-1}^i,S_{k|k-1}^i\right)\\ {\eta }_{k|k-1}^i=H_k{m}_{k|k-1}^i\\ S_{k|k-1}^i=H_k{P}_{k|k-1}^i{H}_k^{\mathrm{T}}+R_k\end{array}$.

3.3. Improved Gaussian Component Weight Merging Principal

By utilizing (43), the posterior intensity of targets can be obtained, and different intensities are represented by different constituent Gaussian component weights. With the aim of tracking space targets with high clutters effectively, label $\kappa$ is denoted in this paper. Then, we have the set $v_k=\left\{w_k,x_k^i,P_k^i,{\kappa }_k^i\right\}$, and ${\kappa }_k^i$ denotes the label of the $k\mathrm{th}$ target at time $i$.

Following the RFS scheme, the constituent Gaussian mixture component weights are modified as follows.

(1)

First, the index of the highest weight component in the target posterior density is selected as follows:

$$i_{\max }=\underset{i\in I}{\arg \max }(w_k^i)$$

(49)

where $I$ is the index set.

(2)

Then, the set of component orders that are most similar to the largest weighted component in the target posterior intensity at time k is represented as

$${\Phi }_k=\left\{i\mid \left({\left(x_k^i-x_k^{i_{\max }}\right)}^{\mathrm{T}}{\left({\mathit{P}}_k^i\right)}^{-1}\left(x_k^i-x_k^{i_{\max }}\right)\le U_{mm}\right)\right\}$$

(50)

$$U_{mm}=\left(1+w_k^i\right){\sigma }_s$$

(51)

where ${\sigma }_s$ is the variance of the measured Gaussian white noise.

(3)

If the Gaussian component satisfies Equation (50) and the sum of the component weights fails to exceed twice the state extraction threshold, then a novel set $v_k^{\tau }$, which contains components, can be obtained as follows:

$${\tilde{l}}_k^{\tau }=i_k^{i_{\max }},{\tilde{w}}_k^{\tau }={\displaystyle \sum _{i\in {\mathsf{\Phi }}_k}{w}_k^i},{\mathit{x}}_k^{\tau }=\frac{1}{{\tilde{w}}_k^{\tau }}{\displaystyle \sum _{i\in {\mathsf{\Phi }}_k}{w}_k^i{\mathit{x}}_k^i}$$

(52)

While ${\displaystyle \sum _{k=1}^{M_k}{\tilde{w}}_k^{\tau }}>1$, the sub-weights of the target are multiplied by a penalized factor except the highest weight $i_{\max }$. Therefore, the novel weights yield the following:

$$w_{{}_{novel}}^k=\{\begin{array}{c}{\tilde{w}}_{{}_k}^{\tau },if\hspace{0.17em}k=i_{\max }\\ \beta {\tilde{w}}_{{}_k}^{\tau },if\hspace{0.17em}k\ne i_{\max },k\in \left[1,\dots ,M_k\right]\end{array}$$

(53)

where $\beta =\alpha (1-w_{{}_k}^{i_{\max }})$ denotes the penalized factor and $0\le \alpha \le 1$ is a constant that can determine the penalized degree.

(4)

Substituting (53) into (48) yields

$${\overline{w}}_{{}_{novel}}^{i,n}=\{\begin{array}{c}\frac{w_{{}_k}^{i,n}}{{\kappa }_k({\lambda }_c)+w_k^{i_{\max }}\beta +{\displaystyle \sum _{i=1}^{N_{k|k-1}}{w}_k^{i,n}}},i\ne i_{\max },i\in [1,M_k],n\in [1,N_{k|k-1}]\\ \frac{w_{{}_k}^{i_{\max },n}\beta }{{\kappa }_k({\lambda }_c)+w_k^{i_{\max }}\beta +{\displaystyle \sum _{i=1}^{N_{k|k-1}}{w}_k^{i,n}}},i=i_{\max },n\in [1,N_{k|k-1}]\end{array}$$

(54)

where $w_{{}_k}^{n,m}=p_{D,k}(x)g_{rcs}(rcs)g_d(d)w_{k|k-1}^j{q}_k^j(z_k)$ and ${\kappa }_k({\lambda }_c)={\lambda }_c{c}_{rcs}(rcs)c_d(d)$.

(5)

Finally, the novel posterior density is computed as follows:

$$v(x)={\displaystyle \sum _{J_k}^{\tau =1}{w}_k^{\tau }N\left(x;{\mathit{x}}_k^{\tau },{\mathit{p}}_k^{\tau }\right)}$$

(55)

$${\tilde{\mathit{P}}}_k^{\tau }=\frac{1}{x_k^{\tau }}{\displaystyle \sum _{i\in {\mathsf{\Phi }}_k}{w}_k^i\left({\mathit{P}}_k^i+\left({\mathit{x}}_k^{\tau }-{\mathit{x}}_k^i\right){\left({\mathit{x}}_k^{\tau }-{\mathit{x}}_k^i\right)}^{\mathrm{T}}\right)}$$

(56)

We can observe from (51) that $U_{mm}$ is an adaptive threshold varying with the measured noise in tracking scenarios. Therefore, the threshold can be dynamically adjusted during the whole tracking process and has a better tracking performance for the closely spaced targets.

3.4. Key Steps of the Proposed HGI-PHD Filter

The main steps of our proposed methods are given as follows. For the sake of readability, the corresponding pseudo-code of our proposed HGI-PHD filter is provided in Appendix B.

Step 1 Calculating the time-varying line-of-sight (LOS) and observing angles.

Step 2 Estimating the GTD scattering parameter via our proposed improved 3D-ESPRIT algorithm and reconstructing the scattering center model.

Step 3 Combining the time-varying observing angles with the scattering center model to predict the real-time scattering echo of space targets and applying an IFFT transform to the real-time electromagnetic echo to obtain the predicted HRRP sequence.

Step 4 Computing the HRRP likelihood for different targets and clutter.

Step 5 Improving the constituent Gaussian component weights by using the method introduced in Section 3.3.

Step 6 Employing the HRRP likelihood and the improved Gaussian component weight merging principal (e.g., ${\overline{w}}_{{}_{novel}}^{i.n}$) in the update recursion of the PHD filter.

4. Simulation Results

In this paper, a merging PHD (HGI-PHD) filter aided by HRRP characteristics and improved with Gaussian mixture components is proposed, which can be applied to track space targets accurately. The time-varying HRRP pseudo-likelihood is utilized in the GM process of the standard PHD filter. The performances of our proposed HGI-PHD filter are validated by the following simulation experiments. Furthermore, the capability of the proposed HGI-PHD filter to discriminate approaching space targets is also demonstrated by utilizing the star-link satellites' real raw data.

4.1. Performance Evaluation of Scattering Center Estimation

In this section, the performance of the modified algorithm with respect to the SNR is compared with the other two methods. We set the scattering center parameters as shown in

Table 1. Scattering parameters.

Scattering Centers	$x_i(m)$	$y_i(m)$	$z_i(m)$	${\alpha }_i$	$A_i$
Scattering center 1	1.202	1.080	1.320	1.0	6.2
Scattering center 2	1.353	1.263	1.541	0.5	5.6
Scattering center 3	1.534	1.702	2.520	0	4.7
Scattering center 4	1.892	2.322	3.210	1.0	3.4

; the matrix beam parameters $P$, $Q$, and $L$ are all set as 6; the SNR varies from 0 dB to 30 dB with an interval of 1 dB; and 200 Monte Carlo trials are performed for each fixed SNR. To simplify, here, we just provide the mean RMSE of the scattering centers obtained by using different methods. Simulation results are depicted in Figure 4. It can be observed from Figure 4a–e that all of the scattering parameter estimation accuracies obtained via different methods increase as SNR increases, and our proposed improved 3D-ESPRIT method has a better parameter estimation performance than the traditional 3D-ESPRIT algorithm, the method in reference [15], and the method in reference [27]. Furthermore, note that the location parameter estimation accuracy is better than that of the type parameter and intensity parameter due to the latter two parameters being dependent on the aforementioned location parameters.

4.2. HRRP Sequence Reconstructed by the GTD Model

To evaluate the reconstruction and parameter estimation performance of the proposed improved 3D-ESPRIT algorithm with simulated data, the back-scattering electromagnetic echo of three typical targets is calculated using electromagnetic computing software. To be specific, the shape and geometrical size of these three typical targets are depicted in Figure 5. In this simulation, we set the frequency number as 101 points with a range of 15.7~17.7 GHz and a frequency step of 20 MHz. Here, the incident azimuth angle and the elevation angle are set as 0° and 90°, respectively.

It can be observed from Figure 6 that for these three typical targets, there are 3, 3, and 2 sharp peaks in the HRRP curve of the original data. Note that the reconstructed HRRP curve of these three targets obtained by using our improved 3D-ESPRIT algorithm matches the whole trend well, and it can also describe the fluctuation characteristics accurately. The above simulation results validate the effectiveness of our improved 3D-ESPRIT method and the feasibility of HRRP reconstruction via the GTD scattering center model.

4.3. HRRP Similarity Rate Comparison

Firstly, the following three-dimensional (3D) simulation scenario is set to obtain the dynamic HRRP of three typical space targets, as depicted in Figure 7. From Figure 7, it can be seen that there are three typical space targets in the surveillance tracking region of the space-based radar.

In most applications, LOS parameter sequences are acquired via different methods. In this paper, we solve satellite targets with TLE ephemeris to compute their positions in the SGP4 model, and then their positions relative to the space-based radar can be calculated eventually. Hence, the LOS parameters can be calculated by using the spacecraft orbit computing software, which is shown in Figure 8.

Table 2. Simulation parameters in the spacecraft orbit computing software.

Parameters	Value
Range of detection $R$	200 km
Orbit altitude of space-based radar $H_0$	500 km
Orbit altitude of target 1 $H_1$	400 km
Orbit altitude of target 2 $H_2$	500 km
Orbit altitude of target 3 $H_3$	480 km
Orbit inclination of space-based radar ${\vartheta }_0$	45°
Orbit inclination of target 1 ${\vartheta }_1$	48°
Orbit inclination of target 2 ${\vartheta }_2$	55°
Orbit inclination of target 2 ${\vartheta }_3$	60°

gives the specific simulation parameters of the space radar and the three typical satellite targets in the spacecraft orbit computing software. Finally, combining the predicted HRRP sequence and the theoretical HRRP sequence of different targets and clutter, we obtain the dynamic coefficients among different targets and clutter as shown in Figure 9. It can be obviously seen that the predicted HRRP and the theoretical HRRP of the same target have the highest similarity rate than those of different targets or clutter. That means the same target has the largest HRRP pseudo-likelihood with its theoretical HRRP sequence, which can be utilized to improve the update accuracy and space target tracking performance.

4.4. Comparisons of OSPA Errors

Suppose there are three space targets existing in the surveillance region with a detection range of $\left[0,100\right]\mathrm{km}\times \left[0,100\right]\mathrm{km}\times \left[0,100\right]\mathrm{km}$. As for these three space targets, all of them obey the variable acceleration motion model, and the kinematic state vector includes the positions, velocities, and accelerations as $z_k={\left[x,y,z,v_x,v_y,v_z,a_x,a_y,a_z\right]}^{\mathrm{T}}$. Moreover, it should be noticed that the measurement vector ${\tilde{x}}_k={\left[x_k,HRR{P}_k\right]}^{\mathrm{T}}$ includes the position measurements and HRRP characteristic measurements.

Let us set the duration for the whole surveillance simulation as 100 s and the survival probability for independent space targets as 0.98. As for the three space targets, they are born at 1s, and the second target vanishes at 70s, while the third one is present throughout the whole tracking process. Then, the pruning procedure is employed with varying time by utilizing the Gaussian weight threshold as $T_{th}={10}^{-5}$ and the maximum number of the mixture components as 100, and we fix the cardinality distribution over 100 terms. Furthermore, the detection probability is set as $P_d=0.98$ for all filters, $V=8\times {10}^6 {\mathrm{m}}^3$ represents the whole space surveillance region volume, and measurements for detected targets are immersed with clutter, which obeys a Poisson RFS distribution and has an intensity of $1.4\times {10}^{-5} {\mathrm{m}}^{-3}$ (i.e., 112 false alarms over space surveillance region per frame). Notice that the HRRP probability density functions of clutter are obtained by adding white Gaussian noise to the HRRP probability density functions of targets.

Here, the optimal sub-pattern assignment (OSPA) metric is employed for the detection performance evaluation metric, and its specific form is given as follows [28]:

$$D_c^p(X,Y)\triangleq {\left[\frac{1}{N}\left(\underset{\ell \in {\Pi }_N^M}{\min }{\displaystyle \sum _{i=1}^M{d}_c{(x_i,y_{\ell (i)})}^p+c^p(N-M)}\right)\right]}^{\frac{1}{p}}$$

(57)

where $d_c(x_i,y_{\ell (i)})\triangleq \min \left(c,\Vert x_i-y_{\ell (i)}\Vert \right)$ is the cut-off distance with parameter $c>0$, $\Vert ·\Vert$ denotes the Euclidean distance, and ${\Pi }_N^M$ represents the set of permutations of cardinality $M$ on $\left\{1,2,\dots ,N\right\}$. $p$ is the order parameter of the OSPA metric with a range of $1\le p\le \infty$, and the OSPA distance is given by $d_c^p(X,Y)=d_c^p(Y,X)$, which corresponds to the case of $M>N$. Then, the average OSPA results after 100 Monte Carlo simulations for detection probabilities with parameters set as $p=1$ and $c=200 \mathrm{m}$ are shown in Figure 10, which is able to demonstrate the performance metrics for the accuracy of target location estimation as well as target number estimation.

It can be seen from Figure 10a,b that the GM-PHD filter has the highest location error and average OSPA cardinality among the four methods in general. Note that for most of the observing time, the amplitude-aided GM-PHD (Am-GM-PHD) filter and the Doppler-aided GM-PHD (Do-GM-PHD) filter have better location estimation accuracy and cardinality estimation accuracy than the conventional GM-PHD filter. We would like to mention that the estimated OSPA location obtained by the Am-GM-PHD filter fluctuates sharply as the time sequence varies. Furthermore, it can be observed that our proposed HRRP-aided filter has lower estimation OSPA results than the other three methods as depicted in Figure 11. We can also notice that our proposed filter fits better with the ground truth data than the other three methods.

The explanation is that, as expected, the proposed filter utilizes the HRRP characteristics of space targets obtained via the wide-band space-based radar system, which is more delicate in describing the target scattering features than the aforementioned ones (e.g., amplitude, etc.). Simulation results show that the proposed HRRP-aided filter outperforms the GM-PHD filter, the Am-GM-PHD filter, and the Do-GM-PHD filter, especially in dense clutter detection scenarios.

4.5. Effective Tracking of Two Closely Spaced Targets

To further exploit the performance of our proposed HGI-PHD filter to track closely spaced multiple targets, two closely spaced targets are chosen for tracking in this subsection. Note that due to Space-X satellites having similar orbit characteristics, two star-link satellites with serial numbers 55713U and 55715U are selected as the closely spaced targets to validate the discrimination functionality of the proposed filter.

The TLE information and observing duration for the two selected star-link satellites are depicted in

Table 3. TLE information of the two selected star-link satellites.

Targets	TLE Information
Satellite 1	1 55713U 23026U 23066.91667824 .00087040 00000-0 95406-3 0 9990 2 55713 43.0022 235.7401 0001445 276.7526 56.7822 15.62262486 1301
Satellite 2	1 55715U 23026W 23066.91667824 .00053204 00000-0 58720-3 0 9995 2 55715 41.0058 235.6894 0002065 272.0442 84.7138 15.62266923 1303

, and the observing duration is from 3 Mar 2023 01:55:39.000 to 3 Mar 2023 01:57:19.000. Moreover, the 3D observation diagram in the spacecraft orbit computing software is shown in Figure 12.

After obtaining the trajectories of the two space targets, the tracking performance of various MTT filters can be compared. Figure 13a,b show the time-averaged OSPA distance and the time-averaged OSPA cardinality over 200 Monte Carlo trials, respectively. Figure 14 shows the mean deviations of the cardinality distributions. As shown in Figure 13 and Figure 14, the proposed HGI-PHD filter produces the best results in the simulation not only in OSPA distance but in cardinality estimation as well. In contrast, the worst performance is from the standard GM-PHD filter, i.e., both its OSPA distance and OSPA cardinality are high in the initial stage of our simulation. Furthermore, it can be also seen that the results obtained from the Am-GM-PHD filter and the Do-GM-PHD filter fluctuate more violently than those obtained from the HGI-PHD filter. Hence, we are convinced by comparing Figure 13 and Figure 14 that the proposed HGI-PHD filter outperforms the other three filters for resolving closely spaced targets.

4.6. Computational Time Analysis

To analyze the computational complexity of different filters, the computational time is compared with an Intel(R) Core(TM) i7-8700 CPU @ 3.20GHz PC with 16 GB RAM and MATLAB R2021b. The core number we have utilized is 6, and the operating system is Windows 10. And the average computational time of the four filters is given in

Table 4. Comparisons of Average Computational Time.

$\mathbf{Clutter} {\mathit{\lambda }}_{\mathit{c}}$	GM-PHD	Am-GM-PHD	Do-GM-PHD	The Proposed Filter
50	0.725 s	3.249 s	0.780 s	3.005 s
100	0.815 s	3.425 s	1.075 s	3.086 s

after 100 Monte Carlo trials. It can be observed the proposed filter has a little heavier computational burden than the Am-GM-PHD filter under two clutter scenarios, which are 3.005 s and 3.086 s corresponding to two different clutter rates. The Do-GM-PHD filter and the standard GM-PHD filter have similar computational costs. In fact, due to computing the prior knowledge of the time-varying HRRP pseudo-likelihood, the proposed filter has a similar computational cost to the Am-GM-PHD filter theoretically, and the results in Table 4 validate it quantitatively.

5. Conclusions

In the standard PHD filter, the tracking performance deteriorates sharply when clutter rates are high, detection probabilities are low, or targets approach each other. To deal with the aforementioned problems, an adapted PHD-based filter aided by HRRP characteristics and with an adapting threshold for the GM component merging procedure was proposed in this paper. By fully using the HRRP information, the electromagnetic characteristics of space targets were utilized in the update recursion of the standard PHD filter. We have also derived the mathematical proof for the GM component merging procedure.

Our proposed HGI-PHD method can be applied to the space situational awareness region. It was shown that our proposed HGI-PHD filter outperformed the GM-PHD filter, the Am-GM-PHD filter, and the Do-GM-PHD filter, especially in scenarios where the detection probabilities are low, clutter rates are high, or there are closely spaced targets. the robustness and effectiveness of the proposed HGI-PHD filter were validated through simulation experiments and using real trajectories. Future work includes the development of efficient real-time algorithms for SAA tasks and the improvement of tracking accuracy for space targets.