Modified Multi-Mode Target Tracker for High-Frequency Surface Wave Radar

Abstract

In high-frequency surface wave radar (HFSWR), part of the radiation signal inevitably propagates upward and illuminates the target through the ionosphere due to the poor controllability of the antenna’s vertical pattern. As a result, a target may have several echoes from different propagation modes, which affect target detection and tracking and is usually unmanageable in existing HFSWR systems. Without the information of the elevation angle, it is difficult to distinguish the propagation mode of measurements during the target detection phase. This paper makes the first attempt to propose a multi-mode target tracker for HFSWR to solve this problem during the target tracking phase. The multipath probability data association (MPDA) tracker is capable of exploiting multipath target signatures of discrete propagation modes and has been widely used. Based on this, we construct a modified multi-mode probability data association tracker for HFSWR to suppress false tracks caused by multiple propagation modes. Numerical simulations demonstrate that this novel tracker can effectively and accurately track the target from the measurements under multiple propagation modes in HFSWR. The processing results of the actual data collected in Weihai, China indicate that this tracker is of great significance for practical applications.

1. Introduction

High-frequency surface wave radar (HFSWR) usually consists of a one-dimensional line array with weak control of the vertical pattern. Therefore, energy cannot be strictly radiated along the sea surface and partially radiates upwards and impinges on the ionosphere. In some cases, the incident skywave energy may be reflected or scattered by the target or ionosphere and return to the receiver [1]. In fact, HF radar signals propagate through the ionosphere in three ways: the first is near-vertical reflection; the second is to emit a wave at an angle less than vertical, reflected by the ocean or the target, and then back along the same path or sea surface; the third is back-scattering phenomenon caused by the irregularities and the fluctuations of the ionosphere [2]. Traditionally, in HFSWR, ionospheric echo is generally considered as clutter, and there have been many studies on ionospheric clutter suppression [2,3]. However, in the second scenario, ionospheric echo may carry the information of the target (such as a plane, island, or vessel).

It is worth noting that the presence of the ionospheric echo of the target can extend the detection range of the HFSWR, which has not been considered or analyzed in the previous HFSWR systems. This paper proposes a novel tracker to solve the problem of false tracks in long-range detection of HFSWR, making it possible to expand the surveillance area of HFSWR. This is very favorable for marine environmental monitoring, and it can be directly applied to existing HFSWR systems without adding additional equipment.

However, the existence of ionospheric echo results in several echoes corresponding to one target (the phenomenon of multiple propagation modes) in HFSWR. Irrespective of ionospheric stratification, there are four propagation modes, as shown in Figure 1: emit and receive the wave through the ground wave path—GG, emit/receive the wave through the ground wave path and receive/emit the wave from the ionosphere—GS/SG, emit the wave to the ionosphere and receive the wave from the ionosphere—SS [4].

The phenomenon of multiple propagation modes results in the following problems:

(a)

One target may have several measurements and one measurement may correspond to many targets, which may cause false detections;

(b)

The actual ground range calculation formula for each propagation modes are different, therefore, misjudgment of the measurement’s propagation mode will lead to a false estimation of the target’s true position; and

(c)

One target may form multiple tracks.

Note that, multipath phenomenon caused by ionospheric stratification also exists. In addition, multiple propagation modes’ phenomenon also brings some advantages. The skywave propagation mode propagates farther than the groundwave propagation mode because its attenuation is smaller and decays slower than the ground attenuation. Hence, if we analyze the echoes of multiple propagation modes together, the detection range of HFSWR can be expanded. However, we need to solve those problems caused by multiple propagation modes’ phenomenon at first. In addition, since the skywave propagation mode echo may exceed the unambiguity range of the existing HFSWR, this paper extends the unambiguity detection range of HFSWR through a range ambiguity resolution method to observe the echoes of all propagation modes more effectively.

The key of the above problems lies in how to judge the propagation mode of the measurements. The direct method is to obtain the echo’s elevation angle which, however, is difficult for HFSWR, due to one-dimensional line array is deployed there. Some researches on two-dimensional arrays has been proposed, for example, Ref. [5] uses L-shaped arrays to suppress ionospheric interference and Ref. [6] uses 4 × 4 square planar arrays to mitigate ionospheric clutter. However, it is still difficult to solve the multi-mode problems in the target detection phase, because mode SG (GS) cannot be distinguished with mode GG (SS) through the elevation angle. Therefore, in this paper, we consider to solve the multi-mode problems when tracking targets for a one-dimensional array HFSWR system.

Target tracking in HFSWR has been extensively studied, and it mainly consists of various filters to describe the target movement properly, such as the Kalman filter [7], deferred decision filter [8], extended Kalman filter [9], and unscented Kalman filter [10], etc. To enhance the tracking performance in clutter background, many data association algorithms are employed to HFSWR, such as near neighbor data association [11], probabilistic data association (PDA) [12] and joint probabilistic data association (JPDA) [13]. There are also some knowledge-based tracking algorithms for reducing the probability of the track breaking [14]. However, all those tracking method cannot solve the multi-mode problem and will result in false tracks.

Multipath phenomenon of over the horizon radar (OTHR) caused by ionosphere stratification also forms false tracks, and there are two common types of methods for tracking. The first type is to establish a target dynamical model based on the radar coordinate system, and applying a data association tracking algorithm to achieve tracking under the radar coordinate system. After tracking, one target may form multiple tracks. Then, the track fusion algorithm is applied to find which tracks correspond to the same target and calculate the real target track. The data association tracking algorithm includes Viterbi data association (VDA) [15], probabilistic data association (PDA) [16] and probabilistic multi-hypothesis tracking (PMHT) [17], etc. The track fusion algorithm includes a dynamic weighted fusion algorithm [18], a sequential track-to-track fusion algorithm [19], etc. The second type establishes the target dynamical model under the geographic coordinate system and completes the data association in the radar coordinate system. This type of method can obtain the real target track directly (only one track for one target). The corresponding data association algorithm includes multipath data association (MPDA) [20], multipath Viterbi data association (MVDA) [21], Markov chain Monte Carlo (MCMC) [22], and expectation maximization data association (EMDA) [23], etc. The tracking accuracy of MCMC and EMDA is high, but their computation is heavy. MPDA uses Markov chain to characterize the probability transfer of target state; MVDA uses the dynamic programming optimization framework to merge or delete data association hypotheses, which is sub-optimal. Those algorithms have both advantages and disadvantages.

The essential difference between the two types of methods is that the target dynamic model is based on different coordinate systems. The first type is based on the radar coordinate system, while the second type is under the geographic coordinate system. In addition, the first type of method does not need the prior knowledge of the ionosphere, and can work stably without relying on coordinate transformation. However, it will form many tracks correspond to one target, which needs track fusion to obtain the real target track. Therefore, this type of method has a higher track loss rate when the echoes of some modes are not detected, while the second type of method requires coordinate transformation and ionosphere status information, which will introduce errors and degrade tracking stability, but the track loss rate is lower. Certain modes’ echoes may not happen because of the path attenuation and the instability of the ionosphere, which will be analyzed in Section 2. Therefore, we choose the second type of method. Since this paper is the first attempt to propose the multi-mode target tracking for HFSWR, we chose to modify the MPDA tracker that has been widely used.

In Section 2, we analyze the multiple propagation modes’ phenomenon in HFSWR, and study the coverage of each propagation mode by analyzing path attenuation. Then, we show the processing results of the actual data collected in Weihai, China to illustrate the existence of multiple propagation modes’ phenomenon (false tracks). In Section 3, we construct a modified multi-mode probability data association tracker. MPDA establishes a target dynamical model and measurement model in the geographic coordinate and the radar coordinate, respectively. Thus, firstly, we show the target dynamical model and measurement model. The target dynamical model uses a linear discrete-time model described by range, velocity, azimuth, and azimuth rate. The measurement model is constructed through the geometrical relationship between target ground range and each propagation mode’s path length and it consists of path length, azimuth, and velocity. Obviously, the measurement model is nonlinear, so the Jacobian of the coordinate transformation is also needed. Moreover, a one-point initiation algorithm proposed in [20] is modified to apply to HFSWR. In Section 4, the simulation results show that the modified multi-mode probability data association tracker can suppress false tracks and track the real target correctly. Moreover, the multi-mode target tracker is successfully applied to the actual data and it found a trajectory, which should be a plane flying to Hohhot, China.

2. Multiple Propagation Modes’ Phenomenon Analysis

Considering that both the groundwave path and the skywave path exist in HFSWR, there should be four propagation modes irrespective of ionospheric stratification, labeled as GG, GS, SG, and SS, as shown in Figure 1. However, the plasma distribution of the ionosphere varies with height, time, solar activity, latitude and longitude, etc. According to the plasma distribution, the ionosphere is usually divided into three layers, the E-layer, the Es-layer, and the F-layer, which may be separated into the F1-layer and F2-layer [24]. The layer from which the wave is reflected is uncertain and is related to the frequency, angle of incidence, and site, etc. Therefore, considering the complexities of the ionosphere, there will be quite a large number of propagation modes. In order to reduce the complexity of the tracking algorithm when taking the ionospheric characteristics into account, we only analyze these five propagation modes in this paper—Mode 1-GG, Mode 2-GS1 (ionospheric height is $h_1$), Mode 3-SG2 (ionospheric height is $h_2$), Mode 4-SS1 (ionospheric height is $h_1$), and Mode 5-SS2 (ionospheric height is $h_2$).

Then, we calculate the attenuation of each mode through their corresponding groundwave attenuation and skywave attenuation. Here, groundwave attenuation is calculated by the groundwave-propagation program GRWAVE developed by Rotheram [25], and skywave attenuation is calculated according to the empirical formula proposed by CCIR (Consultative Committee of International Radio) [26]. Moreover, the path length of each mode is calculated through the geometric relationship, as shown in Figure 2 and its calculation formula is the same as the measurement model that will be given in Section 3.

Figure 3a shows that the path length (the length calculated through the echo’s delay time) of each of five propagation modes is a function of ground range (target’s ground range). Here, the ionospheric heights are $h_1=110 \mathrm{km}$ and $h_2=220 \mathrm{km}$, and the electromagnetic frequency is 5 MHz. Figure 3b shows the path attenuation (the electric field level) in five propagation modes changing with the path length. It can be seen that there is one propagation mode within 200 km, four modes in the range of 200 km to 400 km, and five modes over 400 km. In addition, as the distance increases, the path attenuation is so large that the target cannot be detected. In this case, the number of propagation modes may be less than the theoretical value. The multipath range is also related to electromagnetic frequency and the ionospheric state, which will not be analyzed in detail in this paper.

Next, the actual data collected in Weihai, China is processed to illustrate the existence of the phenomenon of multiple propagation modes’. The processing result is shown in Figure 4, where red dots represent the detection targets with positive speed (the target is approaching the radar), and blue dots represent the detection targets with negative speed (the target is moving away from the radar). As the illuminated area is inland, a large amount of clutter affects the detection performance. Within 900–1100 km, we determine two suspected tracks, represented by green (track A) and yellow (track B) squares, respectively. The measured values of the two suspected tracks are shown in

Table 1. Tracking parameters of Track A.

Time	Path Length/km	Ground Distance/km	Velocity m/s	Azimuth/Degrees
10:57:56	974	934	33.8	−27
10:58:03	973	933	34.1	−27
10:58:10	971	930	33.9	−28
10:58:17	971	931	33.6	−28
10:58:24	971	930	33.3	−27
10:58:31	969	928	33.3	−27

and

Table 2. Tracking parameters of Track B.

Time	Path Length/km	Ground Distance/km	Velocity m/s	Azimuth/Degrees
10:57:49	1013	932	33.6	−29
10:57:56	1013	932	33.2	−28
10:58:03	1010	929	34.4	−27
10:58:10	1013	932	33.8	−25

.

Firstly, we calculate the ground distance of every measurement of the two tracks with two different ionospheric height. We assume the ionospheric heights for tracks A and B are 140 km and 200 km, respectively, then we calculate the ground distance, as shown in the third column of Table 1 and Table 2. Then, we calculate the absolute value of the difference between the two tracks’ parameters at each moment, as shown in

Table 3. The absolute value of the difference between two tracks’ tracking parameters.

Time	Ground Distance/km	Velocity m/s	Azimuth/Degrees
10:57:49	-	-	-
10:57:56	2	0.6	1
10:58:03	4	0.3	0
10:58:10	1	0.1	3
10:58:17	-	-	-
10:58:24	-	-	-
10:58:31	-	-	-

. Using the accuracy of the measurements as a criterion. The accuracy of the distance measurements is 1 km, the accuracy of velocity is 0.6 m/s, and the accuracy of azimuth is five degrees. Since the ionospheric height is not accurate, the accuracy of the distance will increase, set to 5 km. We find that the differences between the two tracks’ ground distance, velocity, and azimuth are smaller than 5 km, 0.6 m/s, and five degrees, respectively, and, according to our criteria, the two tracks are corresponding to the same target. Therefore, there are multiple tracks (false tracks) for one target, i.e., the phenomenon of multiple propagation modes in HFSWR.

After the above analysis, we know that the two tracks should come from one target. The echoes of modes GG/SG1/GS2 of this target have not been detected because path attenuation of the target at about 900 km in these three modes is very high. However, this is also a case of the phenomenon of multiple propagation modes in HSFWR.

3. Modified Multi-Mode PDA Tracker for HFSWR

To suppress false tracks, this section constructs a modified multi-mode PDA tracker, which consists of target dynamical model, measurement model, track initiation, event probabilities and state estimator. For simplicity, we assume a non-maneuvering and constant-velocity aircraft tracking in an ideal ionospheric state.

3.1. Target Dynamical Model

The target state at time k is described in ground coordinates through ground range $L(k)$ (length of the red line 1 in Figure 2), ground range rate $\dot{L}(k)$, azimuth $\theta (k)$ (angle with the center of illuminating aperture angle), azimuth rate $\dot{\theta }(k)$. Therefore, the target status at time k is written as:

$$\mathit{x}(k)=[\begin{array}{cccc}L(k)& \dot{L}(k)& \theta (k)& \dot{\theta }(k)\end{array}]{}^{\prime },$$

(1)

We assume that the target moves along a straight line with a certain deviation which can be expressed as an additive noise term $\mathit{v}(k)$. Then, the equation of target status can be given by a recursion formula:

$$\mathit{x}(k+1)=\mathit{Fx}(k)+\mathit{v}(k),$$

(2)

where the matrix F is given by:

$$\mathit{F}=\left[\begin{array}{cccc}1& T& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& T\\ 0& 0& 0& 1\end{array}\right],$$

(3)

T is the revisit time, $\mathit{v}(k)$ is a white Gaussian sequence with zero-mean and covariance $\mathit{Q}(k)$ which is calculated according to [22]:

$$\mathit{Q}\left(k\right)=\left(\begin{array}{cccc}{T}^3{\delta }_p^2/3& T^2{\delta }_p^2/2& 0& 0\\ T^2{\delta }_p^2/2& T{\delta }_p^2& 0& 0\\ 0& 0& T^3{\delta }_b^2/3& T^2{\delta }_b^2/2\\ 0& 0& T^2{\delta }_b^2/2& T{\delta }_b^2\end{array}\right),$$

(4)

where ${\delta }_p^2$ and ${\delta }_b^2$ are the noise variance of range and bearing filter, respectively.

3.2. Measurement Model

The radar measurement consists of path length $Rs$ (delay time of echo), velocity $Rr$, and azimuth $Az$ (angle with the center of illuminating aperture angle). Hence, the measurement vector at time k is defined as follows:

$$\mathit{y}(k)=[\begin{array}{ccc}Rs(k)& Rr(k)& Az(k)\end{array}]{}^{\prime },$$

(5)

The target measurement is a non-linear function of target status $\mathit{x}(t)$ and associated with propagation mode. Assuming that there are $N_k$ propagation modes, their measurement equation is expressed as:

$${\mathit{y}}^m(k)={\mathit{H}}^m(\mathit{x}(k))+{\mathit{w}}^m(k),m=1,\cdots ,N_k,$$

(6)

where m is the index of propagation mode, ${\mathit{y}}^m(k)$ is the measurement vector of the $m\mathrm{th}$ mode, ${\mathit{w}}^m(t)$ is the measurement noise term of the $m\mathrm{th}$ mode with zero-mean and covariance ${\mathit{R}}^m(k)$. Assuming $\mathit{v}(k)$ and $\mathit{w}(k)$ is unrelated. ${\mathit{R}}^m(k)$ can be obtained as:

$${\mathit{R}}^m\left(k\right)=\mathrm{cov}\left({\mathit{w}}^m\left(k\right)\right)=\left(\begin{array}{ccc}{\delta }_R^2& 0& 0\\ 0& {\delta }_D^2& 0\\ 0& 0& {\delta }_A^2\end{array}\right),$$

(7)

where ${\delta }_R^2$, ${\delta }_D^2$, and ${\delta }_A^2$ are the error variances of range, velocity, and azimuth, respectively.

The measurements are obtained under radar coordinates, ${\mathit{H}}^m(\cdot )$ in measurement equation is the mapping from ground coordinates X to radar coordinates Y. From the geometric relations shown in Figure 2, the coordinate transformation formula from X to Y of five propagation modes 1-GG, 2-GS1, 3-SG2, 4-SS1, and 5-SS2 (

Table 4. Index of propagation modes.

Index	1	2	3	4	5
Mode	GG	GS	SG	SS	SS
$h$	-	$h_1$	$h_2$	$h_1$	$h_2$

) are expressed as follows

$$\{\begin{array}{l}R{s}^m=\mathit{a}(m)L+\mathit{b}(m)\mathit{K}\left(m\right)\\ R{r}^m=\mathit{a}(m)\dot{L}+\mathit{b}(m)\mathit{S}\left(m\right)\\ A{z}^m=\theta \end{array},m=1,\cdots ,5,$$

(8)

$$\mathit{a}=[\begin{array}{ccccc}1& 1/2& 1/2& 0& 0\end{array}],$$

(9)

$$\mathit{b}=[\begin{array}{ccccc}0& 1/2& 1/2& 1& 1\end{array}],$$

(10)

$$\mathit{h}=\left[\begin{array}{ccccc}0& h_1& h_2& h_1& h_2\end{array}\right],$$

(11)

$$\mathit{K}=2\sqrt{{\left(\mu \sin \left(\frac{L}{2\mu }\right)\right)}^2+{\left(\mathit{h}+\mu \left(1-\cos \left(\frac{L}{2\mu }\right)\right)\right)}^2},$$

(12)

$$\mathit{S}=\frac{2\left(\mathit{h}+\mu \right)}{Rg}\sin \left(\frac{L}{2\mu }\right)\dot{L},$$

(13)

where $\mu$ is the earth radius, $\mathit{h}$ contains two different ionospheric height—$h_1$ and $h_2$.

As a result, ${\mathit{H}}^m(\cdot )$ in the measurement equation is as follows:

$${\mathit{H}}^m(x)=\left[\begin{array}{c}\mathit{a}(m)L+\mathit{b}(m)\mathit{K}\left(m\right)\\ \mathit{a}(m)\dot{L}+\mathit{b}(m)\mathit{S}\left(m\right)\\ \theta \end{array}\right],m=1,\cdots ,5,$$

(14)

The Jacobian of transformation from ground coordinates to radar coordinates is defined as follows:

$${\mathit{J}}_m=\frac{\partial \mathit{y}}{\partial \mathit{x}}=\left[\begin{array}{cccc}\mathit{a}(m)+\mathit{b}(m)\mathit{A}\left(m\right)& 0& 0& 0\\ \mathit{b}(m)\mathit{B}\left(m\right)& \mathit{a}(m)+\mathit{b}(m)\mathit{A}\left(m\right)& 0& 0\\ 0& 0& 1& 0\end{array}\right],m=1,\cdots ,5,$$

(15)

$$\mathit{A}=\frac{2\left(\mathit{h}+\mu \right)}{Rs}\sin \left(\frac{L}{2\mu }\right),$$

(16)

$$\mathit{B}=\frac{2\left(\mathit{h}+\mu \right)\dot{L}}{R{s}^2}\left(\cos \left(\frac{L}{2\mu }\right)\frac{Rs}{2\mu }-\sin \left(\frac{L}{2\mu }\right)\frac{Rr}{\dot{L}}\right),$$

(17)

3.3. Initiation

In this paper, we apply a one-point initiation algorithm proposed by [20]. This algorithm forms the initial state estimates for all propagation modes, which can prevent certain modes from being missed. For each valuable measurement $\mathit{y}\left(k\right)=\left[\begin{array}{ccc}Rs& Rr& Az\end{array}\right]$ in scan k, this algorithm forms an initial state estimate ${\widehat{\mathit{x}}}_m\left(0|0\right)=\left[\begin{array}{cccc}L& \dot{L}& \theta & \dot{\theta }\end{array}\right]$ of every propagation mode m (GG, GS1, SG2, SS1, SS2) using the appropriate ionospheric height h, i.e., mapping measurements from Y to X. The initial covariance $\mathit{P}\left(0|0\right)$ is assigned based on the known measurement noise variances $\mathit{R}$. The formula for mapping measurements from Y to X is:

$$GG:\{\begin{array}{l}L=Rs\\ \dot{L}=Rr\\ \theta =Az\\ \dot{\theta }\approx 0\end{array},$$

(18)

$$GS1/SG2:\{\begin{array}{l}L=\frac{2}{h}\left(-Rs\mu +\sqrt{R{g}^2{\mu }^2-\mu h\left(h^2-R{s}^2\right)}\right)\\ \dot{L}=\frac{Rr}{\frac{1}{2}+\frac{h+\mu }{Rs}\sin \left(\frac{L}{2\mu }\right)}\\ \theta =Az\\ \dot{\theta }\approx 0\end{array},$$

(19)

$$SS1/2:\{\begin{array}{l}L=2\mu {\cos }^{-1}\left(\frac{h^2+2{\mu }^2+2\mu h-R{s}^2/4}{2{\mu }^2+2\mu h}\right)\\ \dot{L}=\frac{RrRs}{2\left(h+\mu \right)\sin \left(\frac{L}{2\mu }\right)}\\ \theta =Az\\ \dot{\theta }\approx 0\end{array},$$

(20)

where $h=h_1^{\prime }$ for mode GS1 and SS1, $h=h_2^{\prime }$ for mode SG2 and SS2. Moreover, $h_1^{\prime }$ and $h_2^{\prime }$ are assigned according to the predicted value of the ionospheric height.

The mapping from Y-$Rs$ to X-$L$ of mode GS1/SG2 needs to solve the transcendental equation. For its calculation, here we substitute the cosine with the first two terms of its Taylor’s expansion. Figure 5 shows the transformation errors of modes GS1/SG2 mapping from Y-$Rs$ to X-$L$. It can be seen that the ground range error increases rapidly with ground range, the error within the distance of interest (0 to about 1000 km) is acceptable for track initiation.

3.4. Event Probabilities

Tracking a single target determines the minimum mean-square-error estimator (MMSE) $\widehat{\mathit{x}}\left(k|k\right)$ of the state $\mathit{x}\left(k\right)$ based on the multi-mode measurement process.

The PDA tracker needs to exhaust all possible associations of valid measurements with propagation modes. In order to reduce the number of valid detections, a validation gate (or region) for each mode is needed at each scan k, as shown in Figure 6 (G¹ and G²). In heavy clutter, we define a validation gate for each target state estimate, “a validation region or ‘gate’ in the measurement space in which we expect to see a detection from the target at the next dwell. In this way, many detections that are statistically distant from the expected target location are eliminated from the processing. Detections falling within the validation region are referred to as validated or gated measurements for that target” [22]. For every validation gate $G^m(k)$ of mode m, we can obtain the validated measurement sets $\{{\mathit{y}}_1(k),\cdots ,{\mathit{y}}_{M_k}(k)\}$ by judging whether the measurement $\mathit{y}$ in the k-th scan satisfies:

$$\begin{array}{ll}{G}^m(k)\Rightarrow & \{\mathit{y}\in {\text{ℝ}}^3:[\mathit{y}-{\widehat{\mathit{y}}}^m(k+1|k)]{}^{\prime }{({\mathit{S}}^m(k+1))}^{-1}\\ & [\mathit{y}-{\widehat{\mathit{y}}}^m(k+1|k)]<{\gamma }^m\},m=1,\cdots ,N_k\end{array},$$

(21)

where $N_k$ is the number of propagation modes, the scalar constant ${\gamma }^m$ is determined by the validation gate probability ${P_G}^m(k)$ in the k-th scan based on hypothesis testing of ${\chi }^2(n)$, n is the measurement vector dimension (equal to three in this paper). ${\widehat{\mathit{y}}}^m(k+1|k)$ and ${\mathit{S}}^m(k+1)$ are the predicted measurement and corresponding covariance for mode m respectively, which can be calculated by the following formulas:

$$\widehat{\mathit{x}}(k+1|k)=\mathit{F}\widehat{\mathit{x}}(k|k),$$

(22)

$$\mathit{P}(k+1|k)=\mathit{FP}(k|k){\mathit{F}}^{\prime }+\mathit{Q}(k),$$

(23)

$${\widehat{\mathit{y}}}^m(k+1|k)={\mathit{H}}^m(\widehat{\mathit{x}}(k+1|k)),m=1,\cdots ,N_k,$$

(24)

$${\mathit{S}}^m(k+1)={\mathit{J}}^m\mathit{P}(k+1|k){\mathit{J}}^m(k+1){}^{\prime }+{\mathit{R}}^m(k),m=1,\cdots ,N_k,$$

(25)

where $\widehat{\mathit{x}}(k+1|k)$ is the predicted state that obtained by one-point initiation algorithm, $\mathit{P}(k+1|k)$ is the prediction covariance, and ${\mathit{J}}^m(k+1)$ is the Jacobian of the measurement equation ${\mathit{H}}^m(x)$ evaluated at the current predicted state $\widehat{\mathit{x}}(k+1|k)$ for mode m.

Then, we need to exhaust all possible associations. Let ${\theta }_{i,m}^j\left(t\right)$ denote the i-th possible association hypothesis for the $j-\mathrm{th}\left(j=1,\cdots ,M_k\right)$ validated measurements when the number of active modes at time k is $m\in \left\{0,\cdots ,\min \left(M_k,N_k\right)\right\}$. It satisfies:

$${\theta }_{i,m}^j=\{\begin{array}{ll}theassociatedmode& ,iftarget\\ 0& ,ifclutter\end{array},$$

(26)

Let $r_k\left(m\right)$ denote the number of association hypotheses at time k when the active modes’ number is m.

For example, if there are three valid measurements $\left\{{\mathit{y}}_1\left(k\right),{\mathit{y}}_2\left(k\right),{\mathit{y}}_3\left(k\right)\right\}$ in two validation gates $\left\{G^1\left(k\right),G^2\left(k\right)\right\}$ as shown in Figure 6, there are nine possible association hypotheses: when $m=0$, $r\left(0\right)=1$, $\left\{{\theta }_{1,0}^1=0,{\theta }_{1,0}^2=0,{\theta }_{1,0}^3=0\right\}$; when $m=1$, $r\left(0\right)=4$, $\left\{{\theta }_{1,1}^1=1,{\theta }_{1,1}^2=0,{\theta }_{1,1}^3=0\right\}$, $\left\{{\theta }_{2,1}^1=0,{\theta }_{2,1}^2=1,{\theta }_{2,1}^3=0\right\}$, $\left\{{\theta }_{3,1}^1=0,{\theta }_{3,1}^2=2,{\theta }_{3,1}^3=0\right\}$, $\left\{{\theta }_{4,1}^1=0,{\theta }_{4,1}^2=0,{\theta }_{4,1}^3=2\right\}$; when $m=2$, $r\left(0\right)=5$, $\left\{{\theta }_{1,2}^1=1,{\theta }_{1,2}^2=1,{\theta }_{1,2}^3=0\right\}$, $\left\{{\theta }_{2,2}^1=1,{\theta }_{2,2}^2=2，{\theta }_{2,2}^3=0\right\}$, $\left\{{\theta }_{3,2}^1=1,{\theta }_{3,2}^2=0,{\theta }_{3,2}^3=2\right\}$, $\left\{{\theta }_{4,2}^1=0,{\theta }_{4,2}^2=1,{\theta }_{4,2}^3=2\right\}$; $\left\{{\theta }_{5,2}^1=0,{\theta }_{5,2}^2=2,{\theta }_{5,2}^3=2\right\}$.

Next, for active modes $m=0,1,\cdots ,N_k$ and feasible association hypotheses $i=1,\cdots ,r_k(n)$ at time k, we compute the event probabilities $\overline{\beta }(k)$, ${\beta }_0(k)$ and $\{{\beta }_{i,n}(k)\}$.

$\overline{\beta }(k)$ represents the probability of the event that the target does not exist and all detections in the validation gate are clutter. ${\beta }_0(k)$ represents the probability of the event that the target exists and is observable by all modes, but all measurements in the validation gate are clutter. $\{{\beta }_{i,m}(k)\}$ represents the probability of the event that the target exists and is observable and the measurements fall in the validation gate of mode $m=1,\cdots ,N_k$ [20].

Denote $\mathrm{Pr}({\mathrm{E}}_k|Y^{k-1})$ and $\mathrm{Pr}({\overline{E}}_k|Y^{k-1})$ to represent the target existence possibility and inexistence possibility, respectively. The probabilities are governed by a homogeneous Markov chain with known transition matrix M:

$$\mathit{M}=\left(\begin{array}{cc}\mathrm{Pr}(E_k|E_{k-1})& \mathrm{Pr}({\overline{E}}_k|E_{k-1})\\ \mathrm{Pr}(E_k|{\overline{E}}_{k-1})& \mathrm{Pr}({\overline{E}}_k|{\overline{E}}_{k-1})\end{array}\right),$$

(27)

The propagation Markov chain is:

$$\left[\begin{array}{c}\mathrm{Pr}(E_k|Y^{k-1})\\ \mathrm{Pr}({\overline{E}}_k|Y^{k-1})\end{array}\right]={\mathit{M}}^{\prime }\left[\begin{array}{c}\mathrm{Pr}(E_{k-1}|Y^{k-1})\\ \mathrm{Pr}({\overline{E}}_{k-1}|Y^{k-1})\end{array}\right],$$

(28)

We choose the Poisson model to describe the clutter. Hence, the probability of n clutter measurements in the validation region $G\left(k\right)$ is:

$$g_c\left(n\right)=\frac{{\left(\lambda V_G\left(k\right)\right)}^n{e}^{-\lambda V_G\left(k\right)}}{n!},n=0,1,2,\cdots ,$$

(29)

where, $\lambda$ is the density of detections, $V_G$ is the volumes of validation region $G\left(k\right)$.

Then, the event probabilities are:

$$\overline{\beta }\left(k\right)={\delta }_k^{-1}\mathrm{Pr}\left({\overline{E}}_k|Y^{k-1}\right)\frac{{\lambda }^{M_k}}{M_k!}\exp \left\{-\lambda V_G\left(k\right)\right\},$$

(30)

$${\beta }_0\left(k\right)={\delta }_k^{-1}\mathrm{Pr}\left(E_k|Y^{k-1}\right){\left(1-P_D{P}_G\right)}^{N_k}\frac{{\lambda }^{M_k}}{M_k!}\exp \left\{-\lambda V_G\left(k\right)\right\},$$

(31)

$$\begin{gathered} \begin{array}{cc}\hfill {\beta }_{i,m}\left(k\right)& ={\delta }_k^{-1}\mathrm{Pr}\left(E_k|Y^{k-1}\right){\left(P_D\right)}^m{\left(1-P_D{P}_G\right)}^{N_k-m}\frac{1}{r_k\left(m\right)}\hfill \\ & \times {\displaystyle \prod _{j\in \{j:{\theta }_{i,m}^j\ne 0\}}{\left(2\pi \right)}^{-n_y/2}{\left(\det {\mathit{S}}_q\left(k\right)\right)}^{-1/2}\exp \left\{-\frac{1}{2}{d}_{jq}\right\}}\hfill \\ & \times \{\begin{array}{ll}1& c_m=0\\ \exp \left\{-\lambda V_G\left(k\right)\right\}\frac{{\lambda }^{M_k-m}}{\left(M_k-m\right)!}& c_m>0\end{array}\hfill \end{array} \\ m=1,\cdots ,\min \left\{M_k,N_k\right\},i=1,\cdots ,r_k\left(m\right), \end{gathered}$$

(32)

$$d_{jq}=\left({\mathit{y}}_j\left(k\right)-{\widehat{\mathit{y}}}_q\left(k|k-1\right)\right){}^{\prime }{\mathit{S}}_q{\left(k\right)}^{-1}\left({\mathit{y}}_j\left(k\right)-{\widehat{\mathit{y}}}_q\left(k|k-1\right)\right),$$

(33)

$$q={\theta }_{i,m}^j\left(k\right)1,$$

(34)

where ${\delta }_k$ is a normalization factor at time k, $n_y$ is the dimension of the measurement vector, $P_D$ is the probability of detection, and $P_G$ is the validation gate probability.

3.5. State Estimator

For every possible association hypothesis, the conditional state estimates ${\widehat{\mathit{x}}}_{i,m}(k+1|k+1)$ and covariance ${\mathit{P}}_{i,m}(k+1|k+1)$, $i=1,\cdots ,r_k\left(m\right)$, $m=1,\cdots ,\min \left\{M_k,N_k\right\}$ are calculated according to the following formulas:

$$\begin{array}{ll}{\widehat{\mathit{x}}}_{i,m}\left(k|k\right)& \triangleq \widehat{\mathit{x}}\left(k|k-1\right)+{\displaystyle \sum _{j\in \left\{j:{\theta }_{i,m}^j\ne 0\right\}}{\displaystyle \sum _{l\in \left\{l:{\theta }_{i,m}^l\ne 0\right\}}\mathit{P}\left(k|k-1\right)}}\\ & \times {\mathit{J}}_{qj}^{\prime }{\displaystyle {\sum }_{jl}\left(k\right)\left\{{\mathit{y}}_j\left(K\right)-{\widehat{\mathit{y}}}_{ql}\left(k|k-1\right)\right\}}\end{array},$$

(35)

$$\begin{array}{ll}{P}_{i,m}\left(k|k\right)& \triangleq \mathit{P}\left(k|k-1\right)-{\displaystyle \sum _{j\in \left\{j:{\theta }_{i,m}^j\ne 0\right\}}{\displaystyle \sum _{l\in \left\{l:{\theta }_{i,m}^l\ne 0\right\}}\mathit{P}\left(k|k-1\right)}}\\ & \times {\mathit{J}}_{qj}^{\prime }{\displaystyle {\sum }_{jl}\left(k\right){\mathit{J}}_{ql}\left(k\right)\mathit{P}\left(k|k-1\right)}\end{array},$$

(36)

Next, we form the iterative update equation for the state estimate as expressed in Equation (37) and the state estimation error covariance as expressed in Equation (38):

$$\widehat{\mathit{x}}(k|k)=\{\overline{\beta }(k)+{\beta }_0(k)\}\widehat{\mathit{x}}(k|k-1)+{\displaystyle \sum _{n=1}^{\min (M_k,N_k)}{\displaystyle \sum _{i=1}^{r_k(n)}{\beta }_{i,n}(k){\widehat{\mathit{x}}}_{i,n}(k|k)}},$$

(37)

$$\begin{array}{l}\mathit{P}(k|k)={\beta }_0(k)\mathit{P}(k|k-1)+\overline{\beta }(k)\overline{\mathit{P}}(k)+\text{{}\overline{\beta }(k)+{\beta }_0(k)\}\widehat{\mathit{x}}(k|k-1)\widehat{\mathit{x}}{}^{\prime }(k|k-1)\\ -\mathit{x}(k|k){\mathit{x}}^{\prime }(k|k)+{\displaystyle \sum _{n=1}^{\min (M_k,N_k)}{\displaystyle \sum _{i=1}^{r_k(n)}{\beta }_{i,n}(k)\{{\mathit{P}}_{i,n}(k|k)+{\widehat{\mathit{x}}}_{i,n}(k|k){\widehat{\mathit{x}}{}^{\prime }}_{i,n}(k|k)}}\end{array},$$

(38)

We repeat the above steps of Section 3.4 and Section 3.5 using the update state estimate $\widehat{\mathit{x}}(k|k)$ and the state estimation error covariance $\mathit{P}(k|k)$.

4. Numerical Simulation and Actual Data Processing

In this section, we give the simulation results of the modified multi-mode PDA tracker and then apply the tracker to process the actual data collected in Weihai, China.

Table 5. Simulation parameters.

Label	Parameter	Value
T	Revisit time	10 (s)
L	Total number of iterations	50
$\mathit{x}(0|0)$	The initial state of target	100 km, 0.15 km/s, 5 deg, 0.001 deg/s
${\delta }_p$	Range noise standard deviation	10−6 (km)
${\delta }_d$	Bearing filter noise standard deviation	5 × 10−5 (deg)
${\delta }_R^2$	Range error variance	0.0025 (km2)
${\delta }_D^2$	Velocity error variance	0.0001 ((m/s)2)
${\delta }_A^2$	Azimuth error variance	0.04 (deg2)

lists the simulation parameters of target state model and measurement model.

For the validation gate, the gate probability at the k-th scan $P_G(k)$ is set to 0.95, same for each propagation mode. According to the hypothesis testing of ${\chi }^2(3)$, the scalar constant $\gamma$ should be 7.8. The clutter density is assigned 0.002.

Then, we assume that the probability of target detection $P_D$ is the same for all propagation mode. Markov transition matrix M is:

$$\mathit{M}=\left[\begin{array}{cc}0.75& 0.25\\ 0.05& 0.95\end{array}\right],$$

(39)

When the probability of target detection $P_D$ is equal to 0.4, the tracking result is shown in Figure 7a and the tracking error is given in Figure 7b. The false tracks are suppressed and the target state is estimated correctly.

The actual data collected in Weihai, China has been analyzed in Section 2. Since the detection number is enormous and the suspected target tracks’ locations are known, we just apply the modified multi-mode PDA tracker to the measurements around the suspected two tracks as shown in Figure 8. According to the previous analysis, the two suspected tracks belong to the same target, that is, the multiple propagation modes’ phenomenon. If we use the traditional tracking method, two tracks appear. The tracking result of the proposed tracker is shown in Figure 9 and the arrows indicate the direction of target movement. We can determine that the tracking results are basically consistent with the analysis in Section 2. The tracking result has only one track without false tracks. Therefore, the processing result proves that the modified multi-mode PDA tracker is significant for HFSWR to suppress false tracks caused by the phenomenon of multiple propagation modes.

5. Conclusions

The ionospheric echo of the targets in the HFSWR can be used to detect the targets at further distances, which have not received sufficient attention in the past. However, this will lead to the problem of multi-mode transmission, which is the main issue of this paper. We propose the multi-mode target tracking to solve the multiple propagation modes’ problem causes by the existence of ionospheric echo in HFSWR. The multiple propagation modes’ problem caused several measurements for one target, which will affect target detection and tracking. This problem can be solved by the modified multi-mode PDA tracker. The simulation results show that the tracker can track targets without false tracks. In addition, the tracker is applied to actual data collected in Weihai, China, and the processing result indicates the tracker is significant for HFSWR to solve problem of the multiple propagation modes.