Detection and Tracking of Weak Exoatmospheric Target with Elliptical Hough Transform

Abstract

An elliptical Hough transform (EHT) algorithm is proposed in the framework of track-before-detect (TBD) for joint detection and tracking of weak exoatmospheric targets. The new approach exploits the fact that when restricted to a two-body problem, the exoatmospheric target often follows an elliptical orbit, and thus the Hough transform integrated with orbital geometry information would have better detection performance. The relationship between the original radar measurements in data space and the elliptical parameters in parameter space is explicitly derived with multiple steps of coordinate transformation. It is found that the data points mapping into the parameter space essentially represent a quartic curve. An EHT-based algorithm is then designed, and orbit planarity is also taken into account to reduce the effect of noise accumulation. The influences of primary and secondary thresholds and the signal-to-noise ratio (SNR) on the detection performance are compared by simulations. Additionally, a real radar tracking dataset from a scientific satellite on 28 May 2017 is used to investigate the efficiency of the method. By adding some imaginary clutter to the raw orbit, the results indicate that it is very effective in detecting the real satellite trajectory in a low signal-to-noise ratio (SNR) environment. The advantage of the new method lies in it can not only simultaneously detect and track weak exoatmospheric targets but also can predict the trajectory by using these available detected parameters.

1. Introduction

An important issue in the radar tracking community that arouses wide attention is the detection and tracking of weak (dim) targets in very low signal-to-noise ratio (SNR) environments [1,2,3,4,5]. A key difficulty for solving this lies in the inherent relationship between target detection probability and false alarm probability. To achieve a certain value of detection probability sufficient to maintain the track, one must lower the detection threshold, which in effect inevitably increases the false alarm probability. On the other hand, as the false alarm rate increases (due to clutter, jamming, or thermal noise), traditional recursive trackers, e.g., the Kalman filter (KF) combined with some association scheme [6,7], might have degraded tracking accuracy or even lead to track losses due to the appearance of massively and undesirably validated points.

Ballistic target tracking (BTT) is a good example for illustrating this. During the midcourse (exoatmospheric flight), many countermeasures are adopted to improve the survival capability of the reentry vehicles [8]. For instance, applying radar absorbing material (RAM) [9], shaping the surface, releasing chaff [10], and deployment of repeater jammers [11,12,13] are some notable representatives. Among these, stealth techniques, e.g., shaping, and RAM might reduce the radar cross-section (RCS) of the reentry vehicles to a considerably lower level. Typically, the RCS of an X-band cone-shaped warhead can be reduced to a lower bound of 0.0001 m² or even smaller [8]. When scanning a predetermined large surveillance volume, this will make the radar face a very low SNR environment. In some cases, even with the help of apriori early-warning indications (e.g., handover of the early-warning satellite), the tracking radar still cannot have the capability to find targets of interest in a certain surveillance volume. As a result, the attacker would decrease the range at which the reentry vehicles would be expected to be tracked well by surveillance radar or even lead to frequent loss of its tracks during the exoatmospheric flight.

Instead of awkward working, an alternative way for the victim radar is to employ some sophisticated algorithms, and to change into a track-before-detect (TBD) [14] processing mode, such that the target’s detection and tracking can be made simultaneously at a much greater surveillance range, or at least within the predetermined acquisition range. TBD approaches are quite effective algorithms for detecting weak targets in dense clutter environments in that they process multiple frames of raw data directly instead of single-frame data alone. For TBD processing, the judgment is carried out according to some accumulation metric. When the accumulation parameters are detected, the target tracks are picked up simultaneously. Since the detection is executed in the parameter space, even if the target measurements are not continuous and some points are lost due to signal decay, TBD may still work with a good performance. Therefore, TBD indeed has the advantages of providing a higher probability of detection at the same level of probability of false alarms, while circumventing the troublesome data association problem. TBD has a large family, and many algorithms can fall into this category, such as the Hough transform (HT) [15,16,17,18], maximum likelihood approach (MLA) [19,20], sequential hypothesis testing [21,22], maximum likelihood-probabilistic data association (ML-PDA) [3,4], [23,24], particle filter track-before-detect (PF-TBD) [25,26,27,28,29], 3D matched filtering [30,31], 4D-TBD [32], dynamic programming (DP) [33,34,35], etc. During the past several decades, many extensive studies of TBD have been published in the available literature. Generally, these techniques have some common characteristics. They typically use unthresholded data or thresholded data with significantly lower thresholds and operate on measurement data over several scans in a batch-processing fashion. As a consequence, this usually requires more computational load and has the drawback of not being real-time compared to that of recursive tracking families. Since “a sophisticated software can save sensor hardware” [6], to detect and track low SNR targets at ranges greater than that of a traditional recursive processing mode, it is very important to employ some sophisticated algorithms to extend radar performance under stressful circumstances.

There have been many TBD approaches highlighting detecting weak infrared targets and aircraft targets; however, not many of them have been applied for exoatmospheric applications, especially for those Hough transform ones. The Hough transform is a classical method that has been widely used in radar tracking fields, such as in track initiation [36,37] and joint detection and tracking [15,16,17,18], but very few results have been reported on exoatmospheric applications. The reason partly lies in that the exoatmospheric target often follows a parabolic curve (or more precisely elliptical curve); therefore, the standard Hough transform (SHT) [38] which is intended to detect straight lines, cannot be expected to have better performance. If the accumulation time is sufficiently short, a straight line is a good approximation of an exoatmospheric target trajectory, and the SHT works well. However, it is not the case for long time accumulation (such as the entire midcourse time); in the latter case, the SHT may work badly due to the inherent mismatch of the model. Therefore, one must resort to some innovations wherein the HT can have better detection performance. Note that when restricted to a two-body problem, the motion of an exoatmospheric target is essentially a part of an elliptical orbit, governed by the Keplerian motion equation [39]. What is more natural and rational, as well as practice, is to develop TBD methods specific for exoatmospheric targets that more fully exploit the inherent geometrical or motional characteristics of this kind of object.

The approach we presented is the elliptical Hough transform (EHT) [40]; we extend this algorithm and test it with field data. Since this approach coincides with the theoretical motion model, it is expected that it will have better performance than that of the simple SHT. It must be noted that the EHT we presented differs vastly from those classical ones used in the detection of ellipses in the image processing field, although their names are somewhat similar. For the detection of ellipses in an image, the methods commonly used are the five-dimensional parameter method [41,42], randomized Hough transform (RHT) [43,44], and geometric symmetry method [45,46]. These approaches cannot be directly applied to exoatmospheric applications. For instance, due to the radar line-of-sight limitation, the trajectory of an exoatmospheric target is usually only a part of an elliptical orbit, in contrast to an entire ellipse; therefore, the geometric symmetry method cannot account for this case.

This study is focused on the design of a new Hough detection algorithm. The approach we presented can not only detect a weak reentry vehicle in an earlier midcourse phase but also can predict the orbit with these available detected parameters. More specifically, it has the advantage of the integration of detection, tracking, and prediction capability. Nonetheless, to limit the scope of this investigation, the algorithm should be confined to some hypotheses. One of these hypotheses is the apriori information of early warning, such that the radar can search only limited surveillance space, improve the scan data rate, and reduce the high search dimension of the parameters by a track-while-scan (TWS) mode.

The main contributions of this work are as follows:

(i)

The relationships between the elliptical parameters and the radar raw measurements are established with multiple steps of coordinate transformation. Interestingly, a single data point mapping into the parameter space represents a quartic curve.

(ii)

The EHT detection algorithm is systematically designed, and the planarity of the track is also taken into account to reduce the cumulative effect of noise or clutter.

(iii)

The estimation performances related to the primary and secondary thresholds, measurement error, and SNR are analyzed and verified both by simulations and field data. Thus, a guideline is obtained as to what parameters should be used in similar scenarios.

Specifically, we use the real radar tracking data set of a scientific satellite on the day 28 May 2017 to investigate the efficiency of the method. By adding some imaginary (artificial) clutter to imitate a very low SNR environment, the results indicate that it is very effective in detecting the satellite trajectory in heavy clutter.

The remainder of this paper is organized as follows. Section 2 opens with the elliptical characteristics of exoatmospheric targets. Section 3 establishes the analytical relationship between the raw radar measurements and the defined elliptical parameters. Section 4 presents the EHT algorithm steps in detail. This is followed in Section 5 by simulations and field data analysis. Concluding remarks and recommendations for follow-on studies are provided in the last section.

2. Elliptical Characteristic of Exoatmospheric Targets

The entire trajectory of a ballistic target can be divided into three basic portions: boost, ballistic (coast, midcourse, or exoatmospheric), and reentry [47]. Among these, the boost and reentry phases are relatively short compared with the Earth’s radius, and take place only in close vicinity of Earth. On the contrary, the ballistic phase is an exoatmospheric, free-flight motion, continuing until the atmosphere is reached again. Usually, a typical intercontinental ballistic missile (ICBM) can last about 30 min in the midcourse, which constitutes the main flight time and range of a ballistic target’s motion. Therefore, the defense radar must detect and track this type of target earlier. In particular, it is expected that the target will be detected and tracked during the entire exoatmospheric phase. Under these circumstances, the defense system will have more time to intercept a target and also evaluate its performance of the engagement.

During exoatmospheric flight, the predominant force acting on the target is Earth’s gravity, while the effect of atmospheric drag, Earth non-sphericity, and the gravitational forces due to other celestial bodies can be ignored. Under these assumptions, the dynamics (motion) models of the target are governed by the laws of Keppler [39]. Since the target has a negligible mass relative to the Earth, the gravitational acceleration of the target is the solution of a so-called restricted two-body problem, obtained by Newton’s inverse-square gravity law as

$$a_G(r)=-\frac{\mu }{r^2}{u}_r=-\frac{\mu }{r^3}r$$

(1)

where $r$ is the vector from the Earth’s center to the target, $r\triangleq \left|r\right|$ is its length, $u_r\triangleq r/r$ is the unit vector in the direction of $r$, and $\mu =3.986005\times {10}^{14}\hspace{0.17em}{\mathrm{m}}^3{/\mathrm{s}}^2$ is the Earth gravitational constant.

Using Equation (1) combined with the conservation laws of mechanical energy and moment of momentum, it can be demonstrated that the trajectory follows an elliptical orbit in essence. For a comprehensive treatment of the theories and solutions to the two-body problem, the reader is referred to [39].

To illustrate this, a new coordinate system (CS) is set up (we denote it as THIS-CS), and the target-Earth-radar geometry is shown in Figure 1. In this scenario, the curve $\Gamma$ denotes the target trajectory in the exoatmospheric phase, with $\mathrm{A}$ the thrust cutoff or burnout point, $\mathrm{C}$ the apogee, and $\mathrm{B}$ the reentry point. The inside sphere ${\Omega }_{\mathrm{R}}$ denotes the Earth’s surface, while the outside sphere denotes the atmospheric boundary (in this paper, we simply assume that this boundary coincides with the turnoff and reentry points. The THIS-CS ${\mathrm{O}}_{\mathrm{e}}-xyz$ is a right-handed inertial system, with the origin ${\mathrm{O}}_{\mathrm{e}}$ at Earth’s center, axis ${\mathrm{O}}_{\mathrm{e}}x$ pointing in the AB direction, axis ${\mathrm{O}}_{\mathrm{e}}y$ pointing in the ${\mathrm{O}}_{\mathrm{e}}\mathrm{C}$ direction, and axis ${\mathrm{O}}_{\mathrm{e}}z$ perpendicular to the orbit plane. In this user-defined CS, the elliptical orbit of the target in the exoatmospheric phase can be simply formulated as

$$\Gamma :\hspace{1em}\frac{x^2}{b^2}+\frac{{[y-{(a^2-b^2)}^{1/2}]}^2}{a^2}=1,\hspace{1em}|x|\le (r_{\mathrm{e}}+H_{\mathrm{C}})\sin \left(\theta /2\right),\hspace{0.17em}y\ge 0,\hspace{0.17em}z=0$$

(2)

where a and b are the semi-major and semi-minor axes of the ellipse, respectively, i.e., $a=\mathrm{OC}$ and $b=\mathrm{OD}$. Note that wherein $\mathrm{O}$ is the real ellipse center, while the Earth center is only a focus, with half focal length ${\mathrm{OO}}_{\mathrm{e}}=\sqrt{a^2-b^2}$. Furthermore, $r_{\mathrm{e}}=6378.11\times {10}^3\hspace{0.17em}\mathrm{m}$ is the average Earth radius, $H_{\mathrm{C}}={\mathrm{O}}_{\mathrm{e}}\mathrm{A}-r_{\mathrm{e}}$ is the height of turnoff, $\theta$ is the range angle of the entire midcourse, i.e., $\theta =\angle {\mathrm{AO}}_{\mathrm{e}}\mathrm{B}$. According to classical astrodynamics, these elliptical parameters can be determined by the turnoff parameters [39]

$$\{\begin{array}{l}a=-\frac{\mu r_k}{r_k{V_k}^2-2\mu }\\ b=r_k\cos {\theta }_k{[{\upsilon }_k/(2-{\upsilon }_k)]}^{1/2}\end{array}$$

(3)

$$\theta =2\arcsin \left[\frac{{\upsilon }_k{\cos }^2{\theta }_k\tan {\theta }_k}{\sqrt{1+{\upsilon }_k({\upsilon }_k-2){\cos }^2{\theta }_k}}\right]$$

(4)

where $r_k$, $V_k$, and ${\theta }_k$ are the initial reference (or thrust cutoff) range, velocity, and tilt angle, respectively; ${\upsilon }_k$ is the energy parameter, defined as ${\upsilon }_k\triangleq V_k^2/(\mu /r_k)$. Let $H$ be the radar altitude, note that $H\ll H_{\mathrm{C}}$, then the radar is located somewhere on Earth’s surface, i.e.,

$$x_{\mathrm{R}}^2+y_{\mathrm{R}}^2+z_{\mathrm{R}}^2={(r_{\mathrm{e}}+H)}^2,x_{\mathrm{R}},y_{\mathrm{R}},z_{\mathrm{R}}\in {\Omega }_{\mathrm{R}}$$

(5)

where ${\Omega }_R$ is a spherical region, related to many limitations, such as geopolitics and the natural environment.

From Figure 1, we can draw out some constraint relationships between the elliptical parameters $a$, $b$, and $\theta$.

The reentry point B has the following coordinates

$$\{\begin{array}{l}{x}_{\mathrm{B}}=(r_{\mathrm{e}}+H_{\mathrm{C}})\sin (\theta /2)\\ y_{\mathrm{B}}=(r_{\mathrm{e}}+H_{\mathrm{C}})\cos (\theta /2)\end{array}$$

(6)

Substituting Equation (6) into Equation (2), we obtain

$$\frac{{(r_{\mathrm{e}}+H_{\mathrm{C}})}^2{\sin }^2\left(\theta /2\right)}{b^2}+\frac{{[(r_{\mathrm{e}}+H_{\mathrm{C}})\cos \left(\theta /2\right)-\sqrt{a^2-b^2}]}^2}{a^2}=1$$

(7)

Expansion of Equation (7) and collection of the terms yield the following constraint

$$\begin{array}{ll}\Lambda (a,b,\theta )& \triangleq (a^2-b^2){(r_{\mathrm{e}}+H_{\mathrm{C}})}^2{\cos }^2\left(\theta /2\right)\\ & +2b^2\sqrt{a^2-b^2}(r_{\mathrm{e}}+H_{\mathrm{C}})\cos \left(\theta /2\right)\\ & +b^4-a^2{(r_{\mathrm{e}}+H_{\mathrm{C}})}^2\\ & =0\end{array}$$

(8)

Equation (8) is a quadratic equation related to $\cos \left(\theta /2\right)$. Denote the following symbols as

$$\begin{array}{l}{A}_1\triangleq (a^2-b^2){(r_{\mathrm{e}}+H_{\mathrm{C}})}^2\\ B_1\triangleq 2b^2\sqrt{a^2-b^2}(r_{\mathrm{e}}+H_{\mathrm{C}})\\ C_1\triangleq b^4-a^2{(r_{\mathrm{e}}+H_{\mathrm{C}})}^2\end{array}$$

(9)

Note that we use the subscript to differentiate between symbols in Figure 1. Thus, we have the classic solution as

$$\begin{array}{ll}\cos \left(\theta /2\right)& =\frac{-B_1\pm \sqrt{B_1^2-4A_1{C}_1}}{2A_1}\\ & =\frac{-b^2\sqrt{a^2-b^2}(r_{\mathrm{e}}+H_{\mathrm{C}})\pm \sqrt{a^2(a^2-b^2){(r_{\mathrm{e}}+H_{\mathrm{C}})}^4}}{(a^2-b^2){(r_{\mathrm{e}}+H_{\mathrm{C}})}^2}\\ & =\frac{-b^2}{(r_{\mathrm{e}}+H_{\mathrm{C}})\sqrt{a^2-b^2}}\pm \frac{a}{\sqrt{a^2-b^2}}\end{array}$$

(10)

Next, it can be demonstrated that $\cos \left({\theta }_1/2\right)=\frac{-b^2-a(r_{\mathrm{e}}+H_{\mathrm{C}})}{(r_{\mathrm{e}}+H_{\mathrm{C}})\sqrt{a^2-b^2}}$ is not a valid solution. In fact, the expansion of ${\cos }^2\left({\theta }_1/2\right)$ yields

$$\begin{array}{ll}{\cos }^2\left({\theta }_1/2\right)& =\frac{{\left[-b^2-a(r_{\mathrm{e}}+H_{\mathrm{C}})\right]}^2}{{(r_{\mathrm{e}}+H_{\mathrm{C}})}^2(a^2-b^2)}\\ & =1+\frac{b^2\left[{(r_{\mathrm{e}}+H_{\mathrm{C}})}^2+2a(r_{\mathrm{e}}+H_{\mathrm{C}})+b^2\right]}{{(r_{\mathrm{e}}+H_{\mathrm{C}})}^2(a^2-b^2)}>1\end{array}$$

(11)

Obviously, Equation (11) is a contradiction, since we have a common sense that $0\le {\cos }^2\left({\theta }_1/2\right)\le 1$. It is not difficult to prove that another solution is a feasible one. Thus, the range angle $\theta$ has the following expression:

$$\theta =2\mathrm{acos}\left[\frac{-b^2+a(r_{\mathrm{e}}+H_{\mathrm{C}})}{(r_{\mathrm{e}}+H_{\mathrm{C}})\sqrt{a^2-b^2}}\right]$$

(12)

Similarly, treating a as a variable, and solving Equation (8), and ruling out unreasonable solutions, we can obtain the expression of a as

$$\begin{array}{ll}{a}^2& =-b^2{\mathrm{ctan}}^2(\theta /2)+\frac{b^4\left[1+{\cos }^2(\theta /2)\right]}{{(H_{\mathrm{C}}+r_{\mathrm{e}})}^2{\sin }^4(\theta /2)}\\ & -\frac{2b^3\cos (\theta /2)\sqrt{b^2-{(H_{\mathrm{C}}+r_{\mathrm{e}})}^2{\sin }^2(\theta /2)}}{{(H_{\mathrm{C}}+r_{\mathrm{e}})}^2{\sin }^4(\theta /2)}\end{array}$$

(13)

Using a similar approach, we have the expression $b$ as (Actually, we have four solutions, but only one is valid. For the sake of brevity, we do not list them here)

$$\begin{array}{ll}{b}^2& =a(H_{\mathrm{C}}+r_{\mathrm{e}})-\frac{1}{2}{\cos }^2(\theta /2){(H_{\mathrm{C}}+r_{\mathrm{e}})}^2\\ & +\frac{1}{2}(H_{\mathrm{C}}+r_{\mathrm{e}})\cos (\theta /2)\sqrt{4a^2-4a(H_{\mathrm{C}}+r_{\mathrm{e}})+{\cos }^2(\theta /2){(H_{\mathrm{C}}+r_{\mathrm{e}})}^2}\end{array}$$

(14)

Thus, Equations (13) and (14) have built up the relationship between $a$, $b$, and $\theta$, which is the constraint condition for the presented EHT method.

3. Relationship between Measurements and Elliptical Parameters

In practice, the measurements are usually available in the radar spherical CS at discrete time instants. Next, we establish the relationship between the radar raw measurements and the above-mentioned elliptical parameters through multiple steps of coordinate transformation.

Let $Z\triangleq {(R,A,E)}^{\mathrm{T}}$ be the raw radar measurements in the radar-centered spherical CS (range, azimuth, and elevation), and denote $r_{\mathrm{ENU}}\triangleq {(x_{\mathrm{ENU}},y_{\mathrm{ENU}},z_{\mathrm{ENU}})}^{\mathsf{{\rm T}}}$ the vector in the radar-centered East-North-Up coordinate system (ENU-CS) [47], then the first-step coordinate transformation is

$$\{\begin{array}{l}{x}_{\mathrm{ENU}}=R\cos E\cos A\hfill \\ y_{\mathrm{ENU}}=R\cos E\sin A\hfill \\ z_{\mathrm{ENU}}=R\sin E\hfill \end{array}$$

(15)

Let $L$, $B$, and $H$ be the known radar longitude, latitude, and altitude, respectively, and denote $r_{\mathrm{ECF}}={(x_{\mathrm{ECF}},y_{\mathrm{ECF}},z_{\mathrm{ECF}})}^{\mathrm{T}}$ the vector in the Earth-centered fixed (ECF) CS; then, the coordinate transformation between the ENU-CS and the ECF-CS can be obtained as

$$r_{\mathrm{ECF}}=T_{\mathrm{ENU}}^{\mathrm{ECF}}(r_{\mathrm{ENU}}+{\xi }_1)$$

(16)

where ${\xi }_1={(0,0,r_{\mathrm{e}}+H)}^{\mathrm{T}}$ is a known vector. $T_{\mathrm{ENU}}^{\mathrm{ECF}}$ is a $3\times 3$ orthogonal transformation matrix from the ENU-CS to the ECF-CS, i.e.,

$$T_{\mathrm{ENU}}^{\mathrm{ECF}}={\left[T_{\mathrm{ECF}}^{\mathrm{ENU}}\right]}^{-1}={\left[T_{\mathrm{ECF}}^{\mathrm{ENU}}\right]}^{\mathrm{T}}$$

(17)

Note that $T_{\mathrm{ECF}}^{\mathrm{ENU}}$ is related to the radar position (height $H$, longitude $L$, and latitude $B$) and can be expanded by the product of three rotation matrices

$$\begin{array}{ll}{T}_{\mathrm{ENU}}^{\mathrm{ECF}}& =R_{\mathrm{Z}}(\mathsf{\pi }/2-L)R_{\mathrm{Y}}(-\mathsf{\pi })R_{\mathrm{X}}(\mathsf{\pi }/2+B)\\ & =\left[\begin{array}{l}-\sin (L)\hspace{1em}-\cos (L)\sin (B)\hspace{1em}\cos (L)\cos (B)\\ \cos (L)\hspace{1em}-\sin (L)\sin (B)\hspace{1em}\hspace{0.17em}\sin (L)\cos (B)\\ 0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{1em}\hspace{1em}\cos (B)\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{1em}\sin (B)\end{array}\right]\end{array}$$

(18)

Note that in Figure 1 the ellipse is with reference to an inertial CS, whereas the measurement is obtained in the non-inertial CS. We should build up the relationship between the non-inertial CS and the inertial CS. The coordinate transformation between the Earth-centered inertial CS (ECI-CS) and the ECF-CS is

$$r_{\mathrm{ECI}}=T_{\mathrm{ECF}}^{\mathrm{ECI}}{r}_{\mathrm{ECF}}=R_{\mathrm{Z}}(-\omega \mathrm{t})r_{\mathrm{ECF}}$$

(19)

where $\omega =7.292115 \mathrm{rad}/\mathrm{s}$ is the Earth’s rotation rate, $t$ is the time with $t=0$, indicating the time at thrust-cutoff point $\mathrm{A}$. The expansion of $T_{\mathrm{ECF}}^{\mathrm{ECI}}$ is as

$$T_{\mathrm{ECF}}^{\mathrm{ECI}}=\left[\begin{array}{l}\cos (\omega t)\hspace{1em}-\sin (\omega t)\hspace{1em}0\\ \sin (\omega t)\hspace{1em} \cos (\omega t)\hspace{1em}\hspace{0.17em} 0\\ \hspace{0.17em}\hspace{0.17em}0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}0\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{1em}1\end{array}\right]$$

(20)

The coordinate transformation between the ECI-CS and the launching inertial CS (LCI-CS) is

$$r_{\mathrm{LCI}}=T_{\mathrm{ECI}}^{\mathrm{LCI}}(r_{\mathrm{ECI}}-{\xi }_2)$$

(21)

where ${\xi }_2$ is the thrust cutoff point described in the ECF-CS, i.e., point A in Figure 1, denoting $A_0$ the launching geodetic azimuth (defined as the angle anticlockwise from the launching direction to the meridian line at the thrust cutoff point), $L_{\mathrm{C}}$ and $B_{\mathrm{C}}$ the longitude and latitude of the thrust cutoff, respectively, then ${\xi }_2$ has the following expression when restricted to a spherical Earth model

$${\xi }_2=\left[\begin{array}{l}{x}_{\mathrm{C}}\\ y_{\mathrm{C}}\\ z_{\mathrm{C}}\end{array}\right]=\left[\begin{array}{l}(r_{\mathrm{e}}+H_{\mathrm{C}})\cos B_{\mathrm{C}}\cos L_{\mathrm{C}}\\ (r_{\mathrm{e}}+H_{\mathrm{C}})\cos B_{\mathrm{C}}\sin L_{\mathrm{C}}\\ (r_{\mathrm{e}}+H_{\mathrm{C}})\sin B_{\mathrm{C}}\end{array}\right]$$

(22)

The transformation matrix $T_{\mathrm{ECI}}^{\mathrm{LCI}}$ involved in Equation (21) can be expanded as

$$\begin{array}{ll}{T}_{\mathrm{ECI}}^{\mathrm{LCI}}& =R_{\mathrm{Y}}(-\mathsf{\pi }/2-A_0)R_{\mathrm{X}}(B_{\mathrm{C}})R_{\mathrm{Z}}(L_{\mathrm{C}}-\mathsf{\pi }/2)\\ & =\left[\begin{array}{l}-\sin (A_0)\sin (L_{\mathrm{C}})\hspace{1em}\hspace{1em}0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}0\\ \hspace{1em}\hspace{1em}\hspace{1em}0\hspace{1em}\hspace{1em}\cos (B_C)\sin (L_{\mathrm{C}})\hspace{1em}\hspace{0.17em}\hspace{1em}0\\ \hspace{1em}\hspace{1em}\hspace{1em}0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}0\hspace{1em}\hspace{1em}-\cos (B_C)\sin (A_0)\end{array}\right]\end{array}$$

(23)

The last step of coordinate transformation from the LCI-CS to the THIS-CS defined in Figure 1 is

$$r_{\mathrm{THIS}}=T_{\mathrm{LCI}}^{\mathrm{THIS}}(r_{\mathrm{LCI}}+{\xi }_3)=R_{\mathrm{Z}}(-\theta /2)(r_{\mathrm{LCI}}+{\xi }_3)$$

(24)

where ${\xi }_3={(0,r_{\mathrm{e}}+H_{\mathrm{C}},0)}^{\mathrm{T}}$. Note that the difference between the LCI-CS and THIS-CS is the translation of the origin by range $r_{\mathrm{e}}+H_{\mathrm{C}}$ and the rotation of the CS by angle $\theta /2$.

Figure 2 gives the flow chart of the above-mentioned coordinate transformations, which is relatively complicated. The rudimentary operations are translations and rotations.

Substitution of Equations (16)–(23) into Equation (24) and collections of terms yield the following equation:

$$\begin{array}{ll}{r}_{\mathrm{THIS}}& =T_{\mathrm{LCI}}^{\mathrm{THIS}}\{T_{\mathrm{ECI}}^{\mathrm{LCI}}[T_{\mathrm{ECF}}^{\mathrm{ECI}}{T}_{\mathrm{ENU}}^{\mathrm{ECF}}(r_{\mathrm{ENU}}+{\xi }_1)-{\xi }_2]+{\xi }_3\}\\ & =T_{\mathrm{LCI}}^{\mathrm{THIS}}{T}_{\mathrm{ECI}}^{\mathrm{LCI}}{T}_{\mathrm{ECF}}^{\mathrm{ECI}}{T}_{\mathrm{ENU}}^{\mathrm{ECF}}{r}_{\mathrm{ENU}}+T_{\mathrm{LCI}}^{\mathrm{THIS}}{T}_{\mathrm{ECI}}^{\mathrm{LCI}}{T}_{\mathrm{ECF}}^{\mathrm{ECI}}{T}_{\mathrm{ENU}}^{\mathrm{ECF}}{\xi }_1\\ & \hspace{1em}-T_{\mathrm{LCI}}^{\mathrm{THIS}}{T}_{\mathrm{ECI}}^{\mathrm{LCI}}{\xi }_2+T_{\mathrm{LCI}}^{\mathrm{THIS}}{\xi }_3\\ & \triangleq T_{\mathrm{ENU}}^{\mathrm{THIS}}(r_{\mathrm{ENU}}+{\xi }_1)-T_{\mathrm{ECI}}^{\mathrm{THIS}}{\xi }_2+T_{\mathrm{LCI}}^{\mathrm{THIS}}{\xi }_3\end{array}$$

(25)

where the transformation matrices are expanded as follows

$$\begin{array}{ll}{T}_{\mathrm{ENU}}^{\mathrm{THIS}}& =R_{\mathrm{Z}}(-\theta /2)R_{\mathrm{Y}}(-\mathsf{\pi }/2-A_0)R_{\mathrm{X}}(B_{\mathrm{C}})\\ & \times R_{\mathrm{Z}}(L_{\mathrm{C}}-L-\omega t)R_{\mathrm{Y}}(-\mathsf{\pi })R_{\mathrm{X}}(\mathsf{\pi }/2+B)\end{array}$$

(26)

$$T_{\mathrm{ECI}}^{\mathrm{THIS}}=R_{\mathrm{Z}}(-\theta /2)R_{\mathrm{Y}}(-\mathsf{\pi }/2-A_0)R_{\mathrm{X}}(B_{\mathrm{C}})R_{\mathrm{Z}}(L_{\mathrm{C}}-\mathsf{\pi }/2)$$

(27)

$$T_{\mathrm{LCI}}^{\mathrm{THIS}}=R_{\mathrm{Z}}(-\theta /2)$$

(28)

The full expansions of the above matrices are

$$T_{\mathrm{ENU}}^{\mathrm{THIS}}=\left[\begin{array}{l}\cos (\theta /2)\cos (L-L_{\mathrm{C}}+\omega t)\sin A_0\hspace{1em}\hspace{1em}\hspace{0.17em}0\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}0\\ 0\hspace{1em}-\cos B_{\mathrm{C}}\cos (\theta /2)\cos (L-L_{\mathrm{C}}+\omega t)\sin B\hspace{1em}0\\ 0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-\cos B_{\mathrm{C}}\sin A_0\sin B\end{array}\right]$$

(29)

$$T_{\mathrm{ECI}}^{\mathrm{THIS}}=\left[\begin{array}{l}-\cos (\theta /2)\sin A_0\sin L_{\mathrm{C}}\hspace{1em}0\hspace{1em}0\\ 0\hspace{1em}\cos B_{\mathrm{C}}\cos (\theta /2)\sin L_{\mathrm{C}}\hspace{1em}\hspace{0.17em}0\\ 0\hspace{1em}\hspace{1em}\hspace{1em}0\hspace{1em}\hspace{1em}\hspace{0.17em}-\cos B_{\mathrm{C}}\sin A_0\end{array}\right]$$

(30)

$$T_{\mathrm{LCI}}^{\mathrm{THIS}}=\left[\begin{array}{l}\cos (\theta /2)\hspace{1em}-\sin (\theta /2) \hspace{1em}\hspace{0.17em}0\\ \sin (\theta /2)\hspace{1em}\cos (\theta /2)\hspace{1em}\hspace{1em}\hspace{0.17em}0\\ 0\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{1em}\hspace{1em}0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1\end{array}\right]$$

(31)

Up to this point, Equation (25) has built up the interrelationship between raw radar measurements $Z$ and the elliptical vector $r_{\mathrm{THIS}}={(x,y,z)}^{\mathrm{T}}$. If all the relevant parameters are correct, the projection of the raw observations onto the THIS-CS defined in this paper will be a part of an elliptical curve and in a plane. The full expansion of $r_{\mathrm{THIS}}$ is considerably complicated, but it can still be reformulated as the following form

$$\begin{array}{ll}{r}_{\mathrm{THIS}}& \triangleq \Phi (Z;L_{\mathrm{C}},B_{\mathrm{C}},H_{\mathrm{C}},A_0,\theta ,a,b)\\ & =\Phi ({[R,A,E]}^{\mathrm{T}};L_{\mathrm{C}},B_{\mathrm{C}},H_{\mathrm{C}},A_0,\theta ,a,b)\end{array}$$

(32)

where $\Phi (.)$ is a nonlinear vector function. Or equivalently, we have the implicit equation

$$\Gamma :\phi (Z;a,b,\theta ,A_0,L_{\mathrm{C}},B_{\mathrm{C}},H_{\mathrm{C}})=0,\hspace{0.17em}\hspace{1em}[R,A,E]\in D$$

(33)

Note that in Equation (33) the measurement is $Z={(R,A,E)}^{\mathrm{T}}$, which is usually limited to a preconditioned surveillance region D, and can be obtained directly by the radar with some measurement errors, whereas the parameters $\{a,b,\theta ,A_0,L_C,B_C,H_C\}$ are usually unknown and need to be determined. This is also the focus of our investigation. A simple way to implement the EHT is to map the 3-D raw data into a 7-D parameter space and to make an exhaustive search. Since higher search dimensionality (7-D) requires much more processing load and data storage, and it is not easy to implement in real time, one must resort to some heuristic approaches to reduce the search dimensionality in the first place.

These parameters can be divided roughly into two groups: elliptical parameters $a,b,\theta$ and reference parameters $A_0,L_C,B_C,H_C$. Elliptical parameters are those, in essence, related to the thrust-cutoff parameters $r_k,V_k,{\theta }_k$ and are the parameters we must search via a voting scheme in the parameter space, whereas the reference parameters $A_0,L_C,B_C,H_C$ are not essential and can be roughly obtained by an early-warning radar (or by an early-warning satellite) as apriori.

If the early-warning information is not available, some heuristic methods can also be developed to estimate these parameters. For instance, $A_0$, in essence, determines that $z=0$ in the THIS-CS. Considering the factor of noise, we can search $A_0$ by minimizing the sum of the absolute value $z$. Furthermore, note that the thrust-cutoff point $\mathrm{A}(L_{\mathrm{C}},B_{\mathrm{C}},H_{\mathrm{C}})$ is only a reference point, which is not unique and dependent on the altitude $H_{\mathrm{C}}$ at the thrust cutoff point. For a common practice, $H_{\mathrm{C}}=80 \mathrm{km}$ is a reasonable reference value. Thus, the final reference parameters might be reduced to two-dimensions, i.e., $L_{\mathrm{C}}$ and $B_{\mathrm{C}}$. If the early-warning information is not available, three independent standard polynomial Hough transforms can be used to detect in the separate $R-t$, $A-t$, and $E-t$ planes. The detection results are then combined and transformed to the geodetic CS; the extrapolated value at altitude $H_{\mathrm{C}}=80 \mathrm{km}$ is the requisite reference thrust cutoff point.

In the following sections, for the sake of conciseness and due to space limitations, we assume that the early-warning information $\{A_0,L_{\mathrm{C}},B_{\mathrm{C}},H_{\mathrm{C}}\}=\{A_0,{\xi }_2\}$ is available, but with some estimation errors. Thus, the main work is focused on the search scheme for the elliptical parameters $a,b$, and $\theta$. Actually, according to Equations (13) and (14), it can be seen that $a,b$, and $\theta$ are not independent. Known any of the two, then another parameter is obtained. However, it is difficult to use Equations (13) and (14), since in most cases $\theta$ is usually applied at first. Equations (13) and (14) can be used as constraints to greatly reduce the huge search space.

4. Elliptical Hough Transform

4.1. Data Space and Parameter Space

Note that before transformation Equation (24), all the data expressed in the LCI-CS are explicitly known, except for some measurement errors and early-warning parameter-related errors. However, the transformation Equation (24) is still ambiguous since θ is an undetermined parameter, which must be searched over all possible discrete values, together with the search of parameters $a$ and $b$.

In fact, if we denote

$$r_{\mathrm{LCI}0}\triangleq r_{\mathrm{LCI}}+{\xi }_3\triangleq {(x_{\mathrm{LCI}0},y_{\mathrm{LCI}0},z_{\mathrm{LCI}0})}^{\mathrm{T}}$$

(34)

as a new vector, then the LCI0-CS can be viewed as the LCI-CS, translating from the cutoff point A to the Earth center ${\mathrm{O}}_{\mathrm{e}}$. The vector $r_{\mathrm{LCI}0}$ is also a known quantity related to the measurements.

Using $r_{\mathrm{LCI}0}$ to represent Equation (2) and expanding the rotation matrix $R_{\mathrm{Z}}(-\theta /2)$, we obtain the three-parameter elliptical equation as

$$\begin{array}{l}\Gamma :\hspace{1em}\frac{{\left[\cos \left(\frac{\theta }{2}\right)x_{\mathrm{ECI}0}-\sin \left(\frac{\theta }{2}\right)y_{\mathrm{ECI}0}\right]}^2}{b^2}+\frac{{\left[\sin \left(\frac{\theta }{2}\right)x_{\mathrm{ECI}0}+\cos \left(\frac{\theta }{2}\right)y_{\mathrm{ECI}0}\right]}^2}{a^2}=1\\ \mathrm{s}.\mathrm{t}.\hspace{0.17em}\{\begin{array}{l}\hspace{0.17em}|\cos \left(\frac{\theta }{2}\right)x_{\mathrm{ECI}0}-\sin \left(\frac{\theta }{2}\right)y_{\mathrm{ECI}0}|\le (r_{\mathrm{e}}+H_{\mathrm{C}})\sin \left(\frac{\theta }{2}\right)\\ \sin \left(\frac{\theta }{2}\right)x_{\mathrm{ECI}0}+\cos \left(\frac{\theta }{2}\right)y_{\mathrm{ECI}0}\ge 0\\ z_{\mathrm{ECI}0}=0\end{array}\end{array}$$

(35)

Thus, the EHT maps the points $r_{\mathrm{LCI}0}$ in the converted data space into curves in a $a-b-\theta$ three-parameter space. Note that, when discussed here, $(a,b,\theta )$ are the state variables, whereas $(x_{\mathrm{LCI}0},y_{\mathrm{LCI}0},z_{\mathrm{LCI}0})$ are the parameters. According to Equation (35), it is difficult to judge the curve type in the Hough parameter space. Instead, we use the vector in the THIS-CS to analyze its properties. Using the coordinate relationship

$$\{\begin{array}{l}x=\cos \left(\frac{\theta }{2}\right)x_{\mathrm{ECI}0}-\sin \left(\frac{\theta }{2}\right)y_{\mathrm{ECI}0}\\ y=\sin \left(\frac{\theta }{2}\right)x_{\mathrm{ECI}0}+\cos \left(\frac{\theta }{2}\right)y_{\mathrm{ECI}0}\\ z=z_{\mathrm{ECI}0}\end{array}$$

(36)

to replace Equation (35), it follows from that Equation (35) can be rewritten as

$$\frac{a^2}{y^2}+\frac{b^2}{x^2}-\frac{a^2{b}^2}{y^2{x}^2}=0,\hspace{0.17em}\hspace{0.17em}\theta \ge 2\arcsin \left(\frac{\left|x\right|}{r_{\mathrm{e}}+H_{\mathrm{C}}}\right)$$

(37)

Thus, for a given value $\theta$, the vector $r_{\mathrm{THIS}}$ in THIS-CS can be explicitly computed by using Equation (36); then, the data mapping Equation (37) results in a planar quartic curve, instead of an elliptic curve.

It can be seen that any one of the quartic curves in the Hough parameter space corresponds to the set of all possible elliptical curves in the data space in the THIS-CS through the corresponding data point. If an elliptical curve of points does exist in the THIS-CS, the curve is represented in Hough parameter space as the point of intersection of all of the mapped quartic curves, i.e., in the parameter space, multiple quartic curves will gather at a point, and the parameters at this point are the desirable ones we are searching for.

4.2. Algorithm Design

Next, we give full steps of the algorithm design for the presented EHT. The basic idea is illustrated in Figure 3, and it has four main parts.

(i)

Primary Threshold Processing: The radar operates in a TBD mode, searching a limited number of sectors (surveillance region) $D$. The search is implemented via a sequential mode, i.e., every resolution cell in the region is scanned only once during a search frame time $\mathsf{\Delta }t$. The radar raw pulse compressed signals (range-azimuth-elevation channels) are through a much lower primary threshold $\eta$. In this way, we obtain the original measurement points $\{t,R,A,E\}$, although this dataset contains quite a lot of clutter.

(ii)

Elliptical Hough Transform: Firstly, the values of the reference parameters, as well as the rough bound of the elliptic parameters, are obtained from the apriori information. The parameter space is then quantized in the $a-b-\theta$ dimensions. Any $R-A-E-t$ cell with a value exceeding $\eta$ is mapped into parameter space via multiple steps of coordinate transformation, and its power is reset to count 1. When a primary threshold crossing is mapped into the parameter space, its power (note, we use a binary integrator; hence, its power is 1) is added into the $a-b-\theta$ cells that intersect the corresponding quartic curve in the parameter space. In this way, the accumulator cell at the intersection of multiple quartic curves will have a high value; this is called voting accumulation.

(iii)

Secondary Threshold Processing: A secondary threshold can be applied in the parameter space to declare detections of targets of interest. Note that if the secondary threshold is set too low, there might be multiple detection scatters in neighboring parameter cells. In fact, there indeed may be only one target; it is the quantization of parameter space that leads to the false alarm. To reduce false alarms and meanwhile obtain fine estimation accuracy of the desired parameters, a clustering method can serve this purpose.

(iv)

Inverse Coordinate Transformations: The detected elliptical parameters $a-b-\theta$ are back-transformed into the measurement space $R-A-E-t$ by a series of inverse coordinate transformations, and both detection and tracking results are announced to show the time history and current position of the target. Meanwhile, it is very easy to use the detected parameters to predict the long-term evolutionary behaviors of the trajectory. For instance, the intersection between the detected trajectory and the Earth’s surface is the landing point.

The detailed algorithm to implement the EHT is listed as follows:

Step 1: Search the surveillance space $D$ sequentially and store the $N$ frame raw measurement data. If the power of one $R-A-E-t$ cell exceeds the primary threshold $\eta$, then its position $Z$ is recorded. Thus, during the time interval $N\mathsf{\Delta }t$, we obtain the total measurement data $Z^N\triangleq \left\{Z_{1:N}\right\}$, wherein the $i\mathrm{th}$ frame has $m_i$ measurements, and the total number of measurements in the processing period is $m={\displaystyle \sum _{i=1}^N{m}_i}$.

Step 2: Discretize the three-dimensional parameter space $a\times b\times \theta \in [a_{\mathrm{L}},a_{\mathrm{U}}]\times [b_{\mathrm{L}},b_{\mathrm{U}}]\times [{\theta }_{\mathrm{L}},{\theta }_{\mathrm{B}}]$ (the parameter discretization region can be obtained according to the early warning or determined by experience), and obtain $n_a\times n_b\times n_{\theta }$ cells. The parameter resolutions are

$$\{\begin{array}{l}\mathsf{\Delta }a=(a_{\mathrm{U}}-a_{\mathrm{L}})/n_a\\ \mathsf{\Delta }b=(b_{\mathrm{U}}-b_{\mathrm{L}})/n_b\\ \mathsf{\Delta }\theta =({\theta }_{\mathrm{U}}-{\theta }_{\mathrm{L}})/n_{\theta }\end{array}$$

(38)

and the parameter values can be obtained as

$$\{\begin{array}{l}a(i)=a_{\mathrm{L}}+(i-1/2)\mathsf{\Delta }a, \hspace{1em}i=1,\cdots ,n_a\\ b(j)=b_{\mathrm{L}}+(j-1/2)\mathsf{\Delta }b,\hspace{1em}j=1,\cdots ,n_b\\ \theta (k)={\theta }_{\mathrm{L}}+(k-1/2)\mathsf{\Delta }\theta ,\hspace{1em}k=1,\cdots ,n_{\theta }\end{array}$$

(39)

Set the initial accumulation variable as

$$c(i,j,k)=0,\hspace{0.17em}1\le i\le n_a,\hspace{0.17em}1\le j\le n_b,\hspace{0.17em}1\le k\le n_{\theta },\hspace{0.17em}i,j,k\in \mathbb{N}$$

Step 3: For every $\theta (k),1\le k\le n_{\theta }$ and every data $Z$ total number m), use the apriori information of early warning ${\widehat{A}}_0,{\widehat{L}}_{\mathrm{C}},{\widehat{B}}_{\mathrm{C}},{\widehat{H}}_{\mathrm{C}}$ to calculate the converted data $r_{\mathrm{THIS}}$ via Equation (25). Thus, we have

$${\widehat{r}}_{\mathrm{THIS}}(l,k)={[\widehat{x}(l,k),\widehat{y}(l,k),\widehat{z}(l,k)]}^{\mathrm{T}},\hspace{1em}1\le l\le m;\hspace{1em}1\le k\le n_{\theta }$$

(40)

Step 4: For every ${\widehat{r}}_{\mathrm{THIS}}(l,k)$ and every $a(i)$, compute the semi-minor axis $b(i,k,l)$

$$b(i,k,l)=a(i)\sqrt{\frac{{\widehat{x}}^2(l,k)}{a^2(i)-{\widehat{y}}^2(l,k)}},1\le i\le n_a,1\le k\le n_{\theta },\hspace{0.17em}1\le l\le m$$

(41)

Step 5.1: Scheme I (Binary integration): For every $b(i,k,l)$ and every $b(j)$, if $b(i,k,l)$ subjects to

$$\begin{array}{l}|b(j)-b(i,k,l)|\le \mathsf{\Delta }b/2\\ \{\begin{array}{l}1\le i\le n_a\\ 1\le j\le n_b\\ 1\le k\le n_{\theta }\\ 1\le l\le m\end{array}\end{array}$$

(42)

Then, the count of the accumulation cell $c(i,j,k)$ is added to 1, i.e.,

$$c(i,j,k)=c(i,j,k)+1$$

(43)

Step 5.2: Scheme II (Plane-constrained integration): For every $b(i,k,l)$ and every $b(j)$, if $b(i,k,l)$ subjects to Equation (42), then the accumulation steps change to

$$c(i,j,k)=c(i,j,k)+e^{-q|\widehat{z}(l,k)|}$$

(44)

where q is a design coefficient.

Step 6: Set the secondary threshold ${\eta }^{\mathrm{Hough}}$ in the Hough parameter space, if $c(i,j,k)\ge {\eta }^{\mathrm{Hough}}$, we obtain a detected cell, and the estimation parameters are obtained as $\widehat{a}=a(i)$, $\widehat{b}=b(j)$, and $\widehat{\theta }=\theta (k)$. In certain cases, the detections might be multiple and scatter in several adjacent regions. A k-means clustering method or hierarchical agglomerative clustering method [48] can be employed to combine the neighboring detected cells into one cell.

Step 7: Suppose ultimately that we have detected $n$ targets with parameters $[\widehat{a}(k),\widehat{b}(k),\widehat{\theta }(k)],1\le k\le n$. Using a series of inverse coordinate transformations indexed by Equation (37) back to Equation (15), we can retrieve the positions of targets of interest in any desirable CS. In this way, the presented EHT method declares the detection and tracking results simultaneously. Furthermore, it is also possible to use the estimated parameters to predict the long-term trajectory of the targets to be detected.

Note that, in Equation (44), we have considered the fact that the orbit of an exoatmospheric target in the midcourse usually follows a planar trajectory, hence $\widehat{z}(l,k)$ might be sufficiently small if the target turns out to be the true one. However, for random thermal noises or clutter, they would uniformly intersperse in the surveillance space and might have large values for $\widehat{z}(l,k)$. Therefore, for Scheme II, it indeed has the advantage of improving the accumulation performance for true targets, as well as reducing the false-alarm rate for noises. When in the case $q\to 0$, $e^{-q|\widehat{z}(l,k)|}\to 1$, Scheme II will reduce to Scheme I, whereas when $q\to \infty$ we will have $e^{-q|\widehat{z}(k)|}\to 0$. For the latter case, both targets and noises have no accumulation effect. Because remains a design parameter, in this investigation, we find that when $q=0.01\sim 0.0001$ it might have better accumulation performance.

Now we analyze the computational complexity of the EHT. The points in the measurement space correspond to a quadratic curve in the parameter space. Using apriori information to partition the parameter search space, we obtain a total of $n_a\times n_b\times n_{\theta }$ quantized parameters. Assuming that the total number of measurements is $m$ during N-frame scan, then in theory, the total of $n_a\times n_b\times n_{\theta }\times m$ cumulative calculations will be performed. In fact, since $a,b,\theta$ have the constraint relationship Equation (37), b does not need to be quantified, and the total computation is $n_a\times n_{\theta }\times m$. It follows that the more the number of targets, the denser the clutter, the longer the accumulation time, and the higher the accuracy of parameter quantization, the more time consuming the algorithm is.

Furthermore, in the above steps, we have used a three-dimensional search parameter space $a-b-\theta$ instead of two-dimensional space $a-b$ indicated by Equations (13) and (14). Since Equations (13) and (14) are relatively computationally sophisticated, we only use them as a constraint to reduce the number of search points. That is to say, when implementing the step of discretization Equation (39), we select only the points that satisfy Equations (13) or (14) to implement the mapping. Using Equation (14) as an example, the selection is judged as

$$|b(i,k,l)-\widehat{b}(i,k)|\le \mathsf{\Delta }b$$

(45)

where $b(i,k,l)$ is the measurement-related computed value of $b$, $\widehat{b}(i,k)$ is the theoretical computed value of $b$ according to Equation (14), and $\mathsf{\Delta }b$ is introduced to account for all types of measurement errors. With the steps of preselection, unlike pairs of $a(i)-b(j)-\theta (k)$ can be greatly reduced.

4.3. Detection and False Alarm Probabilities

Next, we analyze the statistical performance of the EHT.

The thermal noise-induced false alarm is considered as Rayleigh distributed in amplitude or exponential in power. Assuming a square law cell detector, then the normalized power probability of distribution (pdf) of noise is

$$f(x)={\mathrm{e}}^{-x}$$

(46)

For the sake of simplicity, we assume the exoatmospheric target follows a Swerling II distribution. Note that in reality, it is rarely the case that a reentry vehicle follows a classical Swerling model; instead, it depends on a lot of elements. Herein, since the Swerling II model has the simplest analytical expression of the probability of detection $P_{\mathrm{d}}$, we use it only to give a clear picture of the statistics of the method. Note that any other target fluctuation models are also applicable to this analysis.

The power pdf of the target is

$$f(x|S)=\exp \left(-x/(1+S)\right)/(1+S)$$

(47)

where $S$ is the SNR. It is easy to obtain the $P_{\mathrm{d}}$ for single $R-A-E-t$ detection cell as

$$P_{\mathrm{d}}={\displaystyle {\int }_{\eta }^{\infty }f(x|S)}\mathrm{d}x=\exp \left(-\eta /(1+S)\right)$$

(48)

where $\eta \triangleq -\ln (P_{\mathrm{fa}})$ is the primary threshold, with $P_{\mathrm{fa}}$ the probability of false alarms for noises. Obviously, the larger the value of $P_{\mathrm{fa}}$ is, the lower the threshold $\eta$ will be. To guarantee the subsequent Hough detection of a target, the primary threshold $\eta$ should be set relatively lower.

To obtain the final probability of detection $P_{\mathrm{d}}^{\mathrm{Hough}}$ per Hough parameter space cell, we follow the same way originated from Carlson [17]. Let n be the total number of $R-A-E-t$ cells that can contribute to a given parameter space cell. There are total $m$ data points exceeding the primary threshold $\eta$. Thus, the probability of this event is $P_{\mathrm{d}}^k{(1-P_{\mathrm{d}})}^{m{n}_a{n}_{\theta }-n}$, where $m{n}_a{n}_{\theta }$ is in essence the maximum accessible number for a given Hough parameter cell. According to [17], we directly give the probability $P_{\mathrm{d}}^{\mathrm{Hough}}$ that a given cell in parameter space exceeds the secondary threshold ${\eta }^{\mathrm{Hough}}$ as follows:

$$\begin{array}{ll}{P}_{\mathrm{d}}^{\mathrm{Hough}}& =P\{c>{\eta }^{\mathrm{Hough}}\}\\ & ={\displaystyle \sum _{n=1}^{m{n}_a{n}_{\theta }}{C}_{m{n}_a{n}_{\theta }}^k}{P}_{\mathrm{d}}^n{(1-P_{\mathrm{d}})}^{m{n}_a{n}_{\theta }-n}P\{y={\displaystyle \sum _{i=1}^n{c}_i}>{\eta }^{\mathrm{Hough}}\}\end{array}$$

(49)

where $c$ denotes the accumulation power for this given parameter cell and $c_i$ is the power of the $i\mathrm{th}$ data point. Actually, for binary integration (Scheme I), $c_i\equiv 1,\hspace{0.17em}1\le n\le m{n}_a{n}_{\theta }$; whereas for plane-constrained integration (Scheme II), $c_i=e^{-q|\widehat{z}(i)|},\hspace{0.17em}1\le n\le m{n}_a{n}_{\theta }$. Obviously, for Scheme I, the final probability of detection is reduced to

$$P_{\mathrm{d}}^{\mathrm{Hough}}=P\{c>{\eta }^{\mathrm{Hough}}\}={\displaystyle \sum _{n=\lfloor {\eta }^{\mathrm{Hough}}\rfloor }^{m{n}_a{n}_{\theta }}{C}_{m{n}_a{n}_{\theta }}^n}{P}_{\mathrm{d}}^n{(1-P_{\mathrm{d}})}^{m{n}_a{n}_{\theta }-n}$$

(50)

where $\lfloor .\rfloor$ denotes the rounding operation. The $P_{\mathrm{d}}^{\mathrm{Hough}}$ for Scheme II is obtained as

$$\begin{array}{ll}{P}_{\mathrm{d}}^{\mathrm{Hough}}& =P\{c>{\eta }^{\mathrm{Hough}}\}\\ & ={\displaystyle \sum _{n=1}^{m{n}_a{n}_{\theta }}{C}_{m{n}_a{n}_{\theta }}^n}{P}_{\mathrm{d}}^n{(1-P_{\mathrm{d}})}^{m{n}_a{n}_{\theta }-n}P\{y={\displaystyle \sum _{i=1}^n{e}^{-q|{\widehat{z}}_{\mathrm{t}}(i)|}}>{\eta }^{\mathrm{Hough}}\}\end{array}$$

(51)

Note that it is difficult to obtain the analytical expression $P\{y={\displaystyle \sum _{i=1}^k{e}^{-q|{\widehat{z}}_{\mathrm{t}}(i)|}}>{\eta }^{\mathrm{Hough}}\}$ since the accumulation statistic $y$ has a sophisticated distribution. For the sake of conciseness, we do not go into the details of the detection and false alarm probabilities in this investigation. Once $P_{\mathrm{d}}^{\mathrm{Hough}}$ is available, the false alarm probability can be obtained by limiting $S\to 0$.

Actually, the detection performance relates to many elements, such as the primary threshold $\eta$, secondary threshold ${\eta }^{\mathrm{Hough}}$, SNR, resolution cell of discretization $\mathsf{\Delta }a\times \mathsf{\Delta }b\times \mathsf{\Delta }\theta$, total processing Plane-constrained integration NΔt, radar measurement errors, and estimation accuracy of early-warning information, etc. Most of these statistics do not have analytical expressions and need massive Monte Carlo simulations to give a quantitative evaluation.

5. Simulations

In this section, we consider an example of a single search radar scenario and investigate the performances of the presented EHT algorithm.

5.1. Scenario Description

The simulation scenario is depicted as follows.

Target: We assume a single reentry vehicle flying in the exoatmospheric phase. The initial reference velocity, altitude, and tilt angle at the thrust cutoff point are 2500 m/s, 80 km, and 40.3070°, respectively, which subsequently determine the total exoatmospheric flight time of approximately 370 s. At the reference time $t=0$, the target is located at 0° N latitude 0°E longitude, pointing to the direction of East. The apogee is 0°N latitude 3.03485°E longitude, and the reentry point is 0° N latitude 6.0697° E longitude.

Radar: The radar is positioned at 2°N latitude 4.5°E longitude with a scan period $\mathsf{\Delta }t=0.5 \mathrm{s}$. The measurement accuracies of range, azimuth, and elevation are ${\sigma }_R=50 \mathrm{m}$ and ${\sigma }_A={\sigma }_E=1 \mathrm{mrad}$, respectively. The radar has received some apriori early warning information from the satellite; thus, the surveillance space is restricted to a limited region of $D=230\times 30\times 30$ km³. The processing time interval is [150, 250] s. Thus, the accumulation frame number is $N=200$, which accounts for a total data time of about 100 s.

False alarm: The number of false alarm detections is Poisson distributed with a spatial density $\lambda =8.69\times {10}^{-5}$/km³, which determines an average number of 18 false alarms of a single scan in the surveillance space D. Note that $\lambda$ is essentially a function of the cell false alarm probability $P_{\mathrm{fa}}$ of the radar and is independent across frames.

Early warning information: The true reference parameters provided by the early warning are $A_0=\mathsf{\pi }/2$, $L_{\mathrm{C}}=0°$, $B_{\mathrm{C}}=0°$, and $H_{\mathrm{C}}=80 \mathrm{km}$. The indicated errors are assumed to be Gauss distributed with standard deviations ${\sigma }_{A_0}=1 \mathrm{mrad}$ and ${\sigma }_{L_{\mathrm{C}}}={\sigma }_{L_{\mathrm{B}}}=0.05°$.

EHT specification: The Hough search space is restricted to $[a_{\mathrm{L}},a_{\mathrm{U}}]=[6600 \mathrm{km}, 6800 \mathrm{km}]$, $[b_{\mathrm{L}},b_{\mathrm{U}}]=[1000 \mathrm{km}, 2000 \mathrm{km}]$, and $[{\theta }_{\mathrm{L}},{\theta }_{\mathrm{U}}]=[5°,7°]$, with a resolution $\mathsf{\Delta }a=2000 \mathrm{m}$, $\mathsf{\Delta }b=2500 \mathrm{m}$, and $\mathsf{\Delta }\theta =0.01°$. The total quantization points are approximate 8 × 10⁶. The judgment of the successful detection of a target is set as ${[{(a-\widehat{a})}^2+{(b-\widehat{b})}^2]}^{1/2}\le 5\times {10}^3 \mathrm{m}$ and $|\theta -\widehat{\theta }|\le 0.01°$. Note that the threshold of judgment remains a design parameter in practice, adjusted/tuned mostly based on engineering experience and intuition. Actually, the optimal threshold is difficult to derive as it is related to specific radar-target geometry and related to many issues.

5.2. Raw Measurements

The radar measurements are available in radar spherical CS at discrete time instants. Figure 4 shows the true target and radar raw measurements in the preconditioned surveillance space D, with Figure 4a the measurements in the radar raw PPI display and Figure 4b the converted measurements in THIS-CS defined in this paper. It can be seen that, while the trajectory of the target in the PPI display is a continuous and complicated curve, in the THIS-CS it indeed is a part of the elliptical curve and symmetrical with the Earth center.

Note that in Figure 4 we do not consider the SNR. The gray-scale plot of all 200-frame primary threshold detections is shown in Figure 5, with darker points denoting higher power. Actually, the threshold $\eta =2.3026$ means an average SNR of 10 dB. It can be seen that the target is more visible as an elliptical curve, but it scatters in a cluttered environment. Conventional recursive trackers usually need an average SNR of 20 dB sufficient to maintain the track, and in this low SNR case, they may lose tracks due to frequent misdetections of the target erroneously associated with unwanted clutter points.

5.3. Integration Effects

Figure 6 and Figure 7 show the accumulation effects in the Hough parameter space, with Figure 6 for Scheme I, Binary integration, and Figure 7 for Scheme II, Plane-constrained integration, $q=0.001$. The target SNR is still 10 dB and the primary threshold is $\eta =2.3026$. For the sake of compactness, we do not show the accumulation effect in the $\theta$ dimension, since it has similar behaviors. It can be seen from Figure 6a and Figure 7a, in the parameter space, we did observe that multiple quartic curves gather together.

The two schemes can concentrate both on the true parameters $(a,b,\theta )=(3397.9 \mathrm{km},1135.4 \mathrm{km},6.0697°)$, but their accumulation effects have slight differences. It can be observed that Scheme II has a better accumulation effect than Scheme I. The reason partly lies in that Scheme II has sufficiently considered orbit planarity as a constraint, wherein the noises located outside the orbit plane have been largely suppressed. For this reason, in the following analysis, unless stated explicitly, we will use Scheme II, i.e., Plane-constrained integration to implement the EHT.

5.4. Parameter Estimation

With the K-means clustering, we can obtain the ultimate estimations of Hough parameters for the EHT. Figure 8 shows the estimated parameters of the EHT with different secondary thresholds ${\eta }^{\mathrm{Hough}}=10$ and ${\eta }^{\mathrm{Hough}}=3$, wherein the symbol “★” denotes the clustering center and “·” denotes the secondary threshold crossings. It can be seen that, when the Hough threshold ${\eta }^{\mathrm{Hough}}$ is relatively high, only a few points can pass the threshold, and the estimation accuracy of parameters will be satisfactory. However, this is not the case for a lower threshold. In Figure 8b, the secondary threshold is set to ${\eta }^{\mathrm{Hough}}=3$, and we can observe that the clustering center deviates vastly from the true value.

Once the parameters in the Hough parameter space are obtained, they can be mapped back into the measurement space R − A − E − t or any other CSs, to show the time history and current position of the target. Figure 9 shows the corresponding detected curve, true trajectory, and raw measurements in a single plot, with Figure 9a the detection threshold ${\eta }^{\mathrm{Hough}}=10$ and Figure 9b ${\eta }^{\mathrm{Hough}}=3$. It can be seen from Figure 9a that, the presented EHT works well as the detection result approximates highly to the nominal true trajectory. In contrast to Figure 9a, Figure 9b shows that it has worse detection performance. Because its secondary threshold is relatively low, multiple crossings will be recorded and this will reduce the clustering performance. However, in reality, the secondary threshold ${\eta }^{\mathrm{Hough}}$ is still a key design parameter and needs to be determined by experience or massive Monte Carlo simulations.

5.5. Performance Comparison

Next, we investigate the detection performance due to the influence of the primary threshold, secondary threshold, radar measurement error, and SNR. This is achieved by massive Monte Carlo runs. In the following simulations, the accumulation time interval is still restricted to [150, 250] s. To assess the performance of the proposed method, 10,000 Monte Carlo numerical simulations were performed in Matlab^® on a computing platform characterized by a 3.6 GHz Intel Core I7–4790 processor with 32 GB RAM.

Figure 10 shows the detection performance concerning the primary threshold $\eta$ and SNR. Figure 11 shows the detection performance concerning the secondary threshold ${\eta }^{\mathrm{Hough}}$ and radar measurement errors.

From Figure 10 and Figure 11, some conclusions are obvious. We briefly summarize them as follows:

(i)

The larger the SNR is, the higher the detection probability obtained;

(ii)

The lower the primary threshold is, the higher the detection probability obtained;

(iii)

The higher the secondary threshold is, the higher the detection probability obtained;

(iv)

The higher the radar measurement accuracy is, the higher the detection probability obtained.

Note that in Figure 11, it does not mean that the higher the secondary threshold, the better. It can be observed that, when ${\eta }^{\mathrm{Hough}}$ is set too high to cause no secondary threshold crossings, it may result in missing alarms, and consequently decrease the probability of detection. Alternatively, one can set no secondary threshold and assume that the peak value corresponds to the desired target. However, this may also induce false alarms, since even if there is no target, the detector will still declare a target.

In fact, the detection performance relates to many elements, such as the primary threshold $\eta$, secondary threshold ${\eta }^{\mathrm{Hough}}$, SNR, resolution cell of discretization $\mathsf{\Delta }a\times \mathsf{\Delta }b\times \mathsf{\Delta }\theta$, total processing time interval $N\mathsf{\Delta }t$, radar measurement errors, and the estimation accuracy of early-warning information. Most of the statistics do not have analytical expressions and need massive Monte Carlo simulations to give a quantitative evaluation.

Table 1. Computational Load with Different Parameters.

False Alarm Rate (Point Number)	nanbnθ = 2003	nanbnθ = 4003	nanbnθ = 6003
Pf = 0.01 (924)	0.3814 s	0.9579 s	1.7943 s
Pf = 0.1 (3557)	1.3723 s	3.9687 s	7.5323 s
Pf = 1 (7042)	2.7169 s	7.8832 s	14.9521 s

shows the average running time of the EHT method with different parameters. The primary false alarm probability is varied in the range $P_{\mathrm{f}}=\{0.01,0.1,1\}$ and the parameter space quantization in the range $n_a{n}_b{n}_{\theta }=\{{200}^3,{400}^3,{600}^3\}$. The target accumulation time interval is still restricted to [150, 250] s. As seen in Table 1, the computational load of the algorithm does increase as the primary detection threshold decreases and the number of quantization units increases. Nonetheless, this processing time is still within the total accumulated time. The total running time will also be greatly reduced if a parallel processing strategy and hardware computational platform are used.

Next, we investigate the performance in a very low SNR environment. Figure 12 gives the detected curves in the THIS-CS with an average SNR of 0 dB, which means that the primary threshold is very low, and most of the clutter points will contribute to the Hough space accumulation, while that of the target points will be greatly decreased. The other parameters are Scheme II, Plane-constrained integration, with $q=0.001$, time interval (150–250) s, primary threshold η = 0, and secondary threshold ${\eta }^{\mathrm{Hough}}=10$.

Transforming Figure 12 to the radar raw R-A-E measurement space with inverse coordinate transforms, we obtain the detection and tracking results. As shown in Figure 13a, the detected trajectory of ETH highly coincides with the true value. The ETH is a TBD method since detection and tracking are announced simultaneously. In fact, although the primary threshold is lowered and the clutter becomes massive, this also means that the probability of detection of the target becomes higher. In order to compare the performance, Figure 13b shows the tracking results of traditional data processing via a detect-before-track (DBT) mode. In Figure 13b, the track initialization is three-point differentiation, the association algorithm is nearest neighbor (NN), and the filtering algorithm is interacting multiple model (IMM). The filter parameters are adjusted to the best. Nevertheless, as can be seen from Figure 13b, the traditional DBT method tracks a large number of temporary tracks. Specifically, the tracks of the target are often intermittent. In fact, there are a total of 18 tracks that belong to the same target, and the total number of temporary tracks reaches 253. This not only wastes radar tracking resources, but also poses a challenge for subsequent target discrimination.

Figure 14 gives the accumulation effect of multiple targets, with Figure 14a the accumulation effect in the Hough parameter and Figure 14b the detection result in THIS-CS. The parameters are still Scheme II, Plane-constrained integration, with $q=0.001$, time interval [150, 250] s, primary threshold $\eta =0$, and secondary threshold ${\eta }^{\mathrm{Hough}}=10$. The difference between these two targets is the difference in release altitude of 10 km. As can be seen, our method can also detect multiple targets at the same time. In fact, the Hough transform is based on a voting mechanism where the accumulation between targets is independent; thus, the Hough transform has implied parallelism, as well as multi-target detection capability. We also find through simulation that our method is able to detect the same target even if the trajectory passes through the surveillance area and then returns to the surveillance area.

5.6. Field Data Analysis

To evaluate the effectiveness of the presented EHT method, we also use filed data to investigate the overall performance. The radar is a phased-array radar, which is manufactured by the China Electronics Technology Group Corporation (CETC) and is operated at the X-band. The data we analyzed involves a constellation comprising three identical scientific satellites of the European Space Agency (ESA), which are designed to provide the best survey of the geomagnetic field and its temporal evolution of the Earth system.

To meet the need for the detection accuracy of the magnetic field of the Earth system, the satellite constellation consists of three circular orbits, with satellites A and C approximately 450 km above the ground, and satellite B 530 km above the ground. Because these data span different periods, we use only the real tracking data of satellite A to investigate the performance. The shape of satellite A is shown on the left of Figure 15.

The orbital elements of satellite A are listed in

Table 2. Parameters of Satellite A.

Symbol	Quantity	True Value
a	semi-major axis	6811.44871 km
e	eccentricity	0.0003368
i	inclination	87.3681°
Ω	right ascension of ascending node	173.0810°
ω	argument of periapsis	96.9899°
f	true anomaly	263.1742°

. It can be seen from Table 2 that, the trajectory is approximately a part of a circular orbit, as the eccentricity $e=0.003368$ is close to zero. According to the basic knowledge of astronomy, we can calculate the true elliptical parameters (a,b,θ) = (6811.4487 km, 6811.4483 km, 0°). Since the target is a satellite and it has a circular orbit, the range angle $\theta$ indeed can be of any value.

Figure 16 shows the raw measurements of the satellite when tracked in the radar centered spherical CS. The data is collected on 8 May 2017. The total observation time is approximately 272 s of its exoatmospheric flight, with the range in the interval [1296 km, 1950 km], azimuth in the interval [225.7384°, 300.6767°], and elevation in the interval [4.4123°,14.6599°]. Since the radar is a phased-array radar and has the ability of adaptive scheduling, the maximum revisiting period is 3.94 s, the minimum period is 0.067 s, and the average is 0.5 s. The measurement accuracies of the range, azimuth, and elevation of the radar are approximate ${\sigma }_R=200 \mathrm{m}$, ${\sigma }_A=0.00228 \mathrm{rad}$, and ${\sigma }_E=0.0011769 \mathrm{rad}$, respectively.

With a series of coordination transformations, the satellite trajectory in the presented THIS-CS can be obtained. As can be seen in Figure 17, the trajectory of the orbit is indeed a part of a circular orbit, with an average height of approximately 449 m. Due to the radar line-of-sight limitation, the radar can only observe a small part of the entire circular orbit.

Since the raw data set is the independent tracking data containing the raw measurements of range, azimuth, and elevation, and for some reason lacking the SNR information, obviously, it cannot be used directly to evaluate the performance of the elliptical Hough transform, as the Hough transform usually needs target data, as well as clutter data.

To investigate the algorithm performance in stressful environments, we add some imaginary (non-real and man-made) clutter data, combining the raw satellite data to constitute a simulated cluttered environment.

The number of clutter detections is still assumed to be Poisson distributed with spatial density $\lambda =8.69\times {10}^{-5}$/km³, which determines an average number of 18 false alarms of a single scan in the surveillance space.

Thus, the field data are preprocessed with the following steps: adding clutter data (imaginary), implementation of the elliptical Hough transform, and inversion to the original CS.

Figure 18 shows the accumulation effects in the Hough parameter space of the synthetic data. The parameters are Scheme II, Plane-constrained integration with $q=0.001$, time interval [0, 272] s, and average SNR = 10 dB. The other parameters are assumed to be the same as earlier. It can be seen from Figure 16 that the corresponding quartic curves intersect at some point and have a high value in the parameter space. In fact, by setting the secondary threshold ${\eta }^{\mathrm{Hough}}=10$ and using the K-means clustering method, we obtain the detected elliptical parameters $(\widehat{a},\widehat{b},\widehat{\theta })=(6808.3778 \mathrm{km}, 6833.7457 \mathrm{km},0.003523°)$, which are now very close to the true values $(a,b,\theta )=(6811.4487 \mathrm{km},6811.4483 \mathrm{km},0°)$.

Once the parameters in the Hough parameter space are obtained, they can be converted back into the measurement CS or any other CSs to show the time history and current position of the target. Figure 19a shows the corresponding detected curve, assumed (imaginary) clutter, and raw satellite measurements in the THIS-CS, while Figure 19b shows the same result in the radar-centered ENU-CS. It can be seen from these two figures that the EHT works well, as the detected trajectory approximates highly to the true measurements of the satellite. An interesting fact in Figure 19 is that there seems to be a fixed deviation in the detection process. Even without adding any clutter (using the original tracking data), this fixed deviation will still exist and cannot be ruled out. The reason partly lies in the fact that a spherical Earth model is employed, which might result in systematic errors and cannot reflect the truth. Perhaps adopting a more precisely ellipsoidal Earth model would mitigate this effect, and this needs further investigation.

6. Conclusions

When tracking low observable targets in the presence of clutter, more challenges will occur. The TBD is a quite effective scheme to deal with low observable targets, as it exchanges the time for obtaining better SNR according to some accumulation metric. In this paper, still within the TBD framework, we present a new Hough detector, i.e., an elliptical Hough transform (EHT) detector, to detect and track weak exoatmospheric targets in the presence of heavy clutter. Both simulations and field data verify the validity of the algorithm. The advantage of the new algorithm lies in that it is specifically designed for exoatmospheric targets; hence, it indeed has the effect of integrating detection, tracking, and even prediction functions. Nonetheless, for multiple target detection, there are still some issues. One problem is that the number of k-means clustering centers still needs to be adjusted by experience. Another issue is that the secondary threshold also has an impact on the multi-target detection performance, especially in the presence of both large and small targets. All these problems need further investigation.

This method can also be extended to other situations. For instance, multisensor integration, data fusion, polarizations, and Doppler measurements can also be added to the approach to improve performance. Furthermore, in addition to using geometrical information of elliptical orbit, other motional information can also be used to design a TBD algorithm to detect an exoatmospheric target. For instance, many exoatmospheric targets have conservation characteristics (conservation laws of mechanical energy and moment of momentum); thus, a simple linear Hough transform detector may be used to detect a physically existent target and discriminate between false targets (e.g., active decoys, which are physically non-existent). The previous work has gained some interesting results concerning joint tracking and discrimination [49,50,51]; future work will be focused on TBD processing in joint detection and tracking and discrimination of weak exoatmospheric targets in clutter and jamming environments. This will be an intellectually motivating field and needs further investigation.