A Multi-Stage Vessel Tracklet Association Method for Compact High-Frequency Surface Wave Radar

Abstract

A compact high-frequency surface wave radar, used for target detection, suffers from a low signal-to-noise ratio, low detection probability, a high false alarm rate, and low positioning accuracy; this is due to its low transmit power and the reduced aperture size of the receiving antenna array. When target tracking algorithms are applied to compact high-frequency surface wave radar data, track fragmentation often occurs and a long track may be broken into several track segments (a.k.a. tracklets), which degrade the tracking continuity for a maritime surveillance system. We present a multi-stage vessel tracklet association method, based on bidirectional prediction and optimal assignment, to associate the broken tracklets belonging to the same target, and connect them to form one continuous track in a multi-target tracking scenario. Firstly, two global motion parameters, i.e., the average heading and average speed, were, respectively, extracted from the newly initiated and terminated tracklets as features for a rough tracklet association, then k-means clustering was used to produce the preliminary tracklet pairs. Subsequently, the temporal and spatial constraints on the initiated and terminated tracklets were considered to refine the preliminary tracklet pairs, to obtain the candidate tracklet pairs. Finally, the tracklet association costs were calculated using Doppler velocity, range, and azimuth to determine the similarity between tracklets in the candidate tracklet pairs, and an association cost matrix was obtained. Then an optimal assignment method based on Jonker–Volgenant–Castanon algorithm was applied to the association matrix to achieve optimal tracklet matching by minimizing the total association costs. Tracklet association experiments with both simulated and field data were conducted; experimental results show that, compared with existing track segment association methods, the association accuracy of the proposed method is significantly improved with better tracking continuity.

1. Introduction

A compact high-frequency surface wave radar (HFSWR) uses vertically-polarized electromagnetic waves of 3–30 MHz to continuously monitor sea surface moving vessels. Due to its advantages, e.g., flexibility in deployment and maintenance, as well as over-the-horizon and all-weather detection capabilities, it has become an indispensable ocean remote sensing sensor [1,2,3]. There are two representative compact HFSWR systems, as introduced in [4,5,6,7,8,9,10,11]; this article focuses on the newly developed compact over-the-horizon radar for maritime surveillance (CORMS) system that uses an eight-element linear receiving antenna array.

However, system miniaturization leads to many problems for target detection and tracking, e.g., a low transmitting power limits its detection range and results in a low signal-to-noise ratio, which reduces its detection probability. The reduced aperture size of the receiving antenna array causes wide beamwidth and coarse azimuth resolutions. Thus, its positioning accuracy drops. A long coherent integration time (CIT) for vessel detection leads to a low data rate. Therefore, the compact HFSWR system suffers from a low signal-to-noise ratio, low detection probability, low data rate, low positioning accuracy, as well as a high false alarm rate. These disadvantages bring about great challenges in target detection and tracking. When target tracking algorithms are applied to the data obtained by compact HFSWR, track fragmentation, which means that the track of a target is broken into several track segments (a.k.a. tracklets), often occurs and leads to poor tracking continuity. The objective of this paper was to develop a tracklet association method to address the track fragmentation problem, and improve the tracking continuity for compact HFSWR systems.

Track fragmentation (or breakage) is mainly caused by an incorrect plot-to-track association during the tracking procedure. An incorrect plot-to-track association may be due to missed detections because of low detection probability, target motion model mismatch, etc. Thus far, many researches have attempted to address the track fragmentation problem by proposing various track segment association (TSA) methods. Existing tracklet association methods can be categorized into three types. The first type of methods uses statistical methods to calculate the target state similarity between tracklets. For example, a track segment association method based on statistical weighting was proposed in [12]. In that method, a backward prediction was exerted on an initiated tracklet until the last state of a terminated tracklet, then an association condition relying on the Mahalanobis distance between the predicted state of the initiated tracklet, and the last state of the terminated tracklet was tested to determine the relevance. Similarly, target identification features as described in [13] were introduced into the tracklet association method to improve the tracklet association performance. Moreover, a turning model [14], the interactive multiple models [15], and a multi-model global nearest neighbor method [16] were employed to repair the broken tracks caused by target maneuvering. Moreover, incorrectly associated measurements were released first, then a multi-frame assignment-based TSA method was proposed to estimate the track during the breakage period using both unassociated and released measurements [17]. The second type includes fuzzy mathematics-based methods. The membership matrix [18] or similarity matrix [19] were calculated first using different target features, and then clustering methods [18,19], fuzzy K-nearest neighbor, and fuzzy C-means methods [20] were applied to achieve tracklet association. The third category involves artificial intelligence-based methods, e.g., an extreme learning machine was used for the HFSWR tracklet association and it achieved improved accuracy [21].

Due to the disadvantages of target detection with compact HFSWR, as well as target maneuvering and the complexity of tracking environments, existing TSA methods usually suffer from large errors in feature extraction, backward prediction, etc., when they are directly applied to compact HFSWR data; thus, they may not obtain satisfactory tracklet association results. To improve the tracklet association accuracy for a compact HFSWR system, a multi-stage tracklet association method was proposed in this article. In this method, the global target motion characteristics, temporal and spatial constraints, and kinematic parameters were employed to enhance the tracklet discrimination. K-means clustering and bidirectional prediction were used to improve the efficiency and accuracy. Moreover, the Jonker–Volgenant–Castanon (JVC) algorithm was employed to achieve optimality in multi-target tracking scenarios. Experiments with both simulated and field data were conducted to verify the accuracy and efficiency of the proposed method. The remainder of this paper is organized as follows. Tracklet representations as well as average heading and speed estimation methods are described in Section 2. The proposed multi-stage tracklet association method is introduced in detail in Section 3. In Section 4, experimental results are presented and analyzed; conclusions are drawn in Section 5.

2. Preliminaries

2.1. Tracklet Representation

The compact HFSWR transmits a linear frequency-modulated interrupted continuous wave (FMICW) to illuminate the sea surface within the coverage area. The backscattered echoes are received by a linear array of antennas. The signal received by each antenna is digitally processed to attain the range and Doppler velocity information; thus, a range–Doppler map is obtained. Then, a constant false alarm rate (CFAR) algorithm is applied to the range–Doppler map to achieve target detection, and a direction of arrival (DOA) estimation method, such as multiple signal classification (MUSIC) or digital beamforming (DBF) is used to obtain the azimuth of the detected target. Therefore, a compact HFSWR locates a target in terms of range r and azimuth $\theta$ under the polar coordinate system with the radar site as its origin. Moreover, as a Doppler radar, it can measure the target velocity component $v_r$ along the radar radial direction (a.k.a. Doppler velocity). Thus, it represents a target with a state vector ${\left[v_{r}r\theta \right]}^T$. HFSWR continuously observes the sea surface, at each sampling instant k, it acquires a frame of data containing plenty of “plots” from echoes of different targets and interferences. After several consecutive sampling periods, multiple frames of plot data can be collected. Then, a multi-target tracking algorithm can be applied to the obtained plot data sequence to produce target tracks.

In general, a multi-target tracking algorithm consists of three steps, i.e., track initiation, track maintenance, and track termination. Track maintenance includes state prediction, measurement-to-track association, and state estimation. The converted measurement Kalman filter (CMKF) and the minimal cost data association method were combined for target tracking in this paper and are described as follows.

(1) The dynamic and observation models.

The converted measurement Kalman filter operates with a dynamic model and an observation model. The dynamic model of moving vessels can be defined in a Cartesian coordinate system as

$${\mathbf{s}}_t={\mathbf{Hs}}_{t-1}+{\mathsf{\omega }}_t,$$

(1)

where ${\mathbf{s}}_t={[x_t,v_{x_t},y_t,v_{y_t}]}^T$ is the target’s true state vector at time t, $x_t$ and $y_t$ are the target’s position components, $v_{x_t}$ and $v_{y_t}$ denote the target’s true velocity components along the x and y directions, respectively. $\mathbf{H}$ is the state transition matrix, defined as

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

(2)

and T is the sampling time. ${\mathsf{\omega }}_{\mathbf{t}}$ denotes the Gaussian process noise with a mean of zero and covariance matrix ${\mathbf{Q}}_t$.

The observation model is also defined in the Cartesian coordinate system as

$${\tilde{\mathbf{s}}}_t=\mathbf{M}{\mathbf{s}}_t+{\mathbf{v}}_t,$$

(3)

where ${\tilde{\mathbf{s}}}_t={[{\tilde{x}}_t,{\tilde{v}}_{x_t},{\tilde{y}}_t,{\tilde{v}}_{y_t}]}^T$ is the target’s measured state vector at time t, ${\tilde{x}}_t$ and ${\tilde{y}}_t$ represent the measured target’s position components, ${\tilde{v}}_{x_t}$ and ${\tilde{v}}_{y_t}$ denote the corresponding measured velocity components along x and y directions. $\mathbf{M}$ is the measurement matrix and ${\mathbf{v}}_t$ represents measurement noise following Gaussian distribution with a mean of zero and covariance matrix ${\mathbf{R}}_t$.

(2) Track initiation.

Potential tracks are initiated using the logic method with the M-of-N rule [22]. If there are more than M plots connected in the most recent N frames, the track is successfully initiated, and it will be added for track maintenance; otherwise, it will be discarded.

(3) Track maintenance.

A. State prediction. For each initiated or maintained track, denote ${\widehat{\mathbf{s}}}_{t-1}=[{\widehat{x}}_{t-1},{\widehat{v}}_{x_{t-1}},$${\widehat{y}}_{t-1},{\widehat{v}}_{y_{t-1}}{]}^T$ as the estimated target state at time $t-1$, the predicted state ${\widehat{\mathbf{s}}}_{t\mid t-1}=[{\widehat{x}}_{t\mid t-1},$${\widehat{v}}_{x_{t\mid t-1}},{\widehat{y}}_{t\mid t-1},{\widehat{v}}_{y_{t\mid t-1}}{]}^T$ at time t can be obtained by ${\widehat{\mathbf{s}}}_{t\mid t-1}=\mathbf{H}{\widehat{\mathbf{s}}}_{t-1.}$ In addition, the corresponding state prediction covariance matrix ${\mathbf{P}}_{t\mid t-1}$ is calculated by ${\mathbf{P}}_{t\mid t-1}=\mathbf{H}{\mathbf{P}}_{t-1}\mathbf{H}$$+{\mathbf{Q}}_t.$

B. Coordinate conversion. For the subsequent measurement-to-track association procedure, the predicted state ${\widehat{\mathbf{s}}}_{t\mid t-1}$ is converted from the Cartesian coordinates to polar coordinates as ${\left[\begin{array}{ccc}{R}_t^p\hfill & {\theta }_t^p\hfill & v_t^p\hfill \end{array}\right]}^T$, in which

$$R_t^p=\sqrt{{\widehat{x}}_{t\mid t-1}^2+{\widehat{y}}_{t\mid t-1}^2},$$

(4)

$${\theta }_t^p=arctan\left(\frac{{\widehat{y}}_{t\mid t-1}}{{\widehat{x}}_{t\mid t-1}}\right),$$

(5)

$$v_t^p=\frac{{\widehat{x}}_{t\mid t-1}{\widehat{v}}_{x_{t\mid t-1}}+{\widehat{y}}_{t\mid t-1}{\widehat{v}}_{y_{t\mid t-1}}}{\sqrt{{\widehat{x}}_{t\mid t-1}^2+{\widehat{y}}_{t\mid t-1}^2}},$$

(6)

where $R_t^p$, ${\theta }_t^p$, and $v_t^p$ denote the predicted range, azimuth, and radial velocity, respectively.

C. Measurement-to-track association. The minimal cost criterion is utilized to find the most likely measurements ${\left[\begin{array}{ccc}{R}_t^m\hfill & {\theta }_t^m\hfill & v_t^m\hfill \end{array}\right]}^T$ for the current track at time t within a predefined validation gate [22]. If the measurement is associated with a track, go to step D; otherwise, go to step F.

D. Measurement conversion. The associated measurement ${\left[\begin{array}{ccc}{R}_t^m\hfill & {\theta }_t^m\hfill & v_t^m\hfill \end{array}\right]}^T$ is converted from polar coordinates to Cartesian coordinates to obtain the measured target state ${\tilde{\mathbf{s}}}_t={\left[{\tilde{x}}_t,{\tilde{v}}_{x_t},{\tilde{y}}_t,{\tilde{v}}_{y_t}\right]}^T$ by

$${\tilde{x}}_t=R_t^{m}cos{\theta }_t^m,$$

(7)

$${\tilde{y}}_t=R_t^{m}sin{\theta }_t^m,$$

(8)

$${\tilde{v}}_{x_t}=\left({\tilde{x}}_t-{\tilde{x}}_{t-1}\right)/T,$$

(9)

$${\tilde{v}}_{y_t}=\left({\tilde{y}}_t-{\tilde{y}}_{t-1}\right)/T.$$

(10)

E. State estimation. The estimated target state ${\widehat{\mathbf{s}}}_t$ at time t and the state estimation covariance matrix ${\mathbf{P}}_t$ are updated by

$${\mathbf{K}}_t={\mathbf{P}}_{t\mid t-1}{\mathbf{M}}^T\left(\mathbf{M}{\mathbf{P}}_{t\mid t-1}{\mathbf{M}}^T+{\mathbf{R}}_t\right),$$

(11)

$${\widehat{\mathbf{s}}}_t={\widehat{\mathbf{s}}}_{t\mid t-1}+{\mathbf{K}}_t\left({\tilde{\mathbf{s}}}_t-\mathbf{M}{\widehat{\mathbf{s}}}_{t\mid t-1}\right),$$

(12)

$${\mathbf{P}}_t={\mathbf{P}}_{t\mid t-1}-{\mathbf{K}}_k\mathbf{M}{\mathbf{P}}_{t\mid t-1},$$

(13)

where ${\mathbf{K}}_t$ is the Kalman gain at time t. Then the estimated target state ${\widehat{\mathbf{s}}}_t$ is used to update the current track.

F. Determine if the track termination conditions are satisfied. If the conditions are met, the track will be terminated; otherwise, t is increased by 1 and go to A.

(4) Track termination.

A maintained track will be terminated if one of the following conditions occurs:

A. There are no associated measurements in the past K frames out of the most recent L frames.

B. The estimated velocity reaches an unrealistic value $v_{max}$.

(5) Track smoothing.

The obtained tracks by the above tracking procedure usually fluctuate significantly and deviate from their true positions due to a low positioning accuracy of the compact HFSWR. Denote the position data sequence in longitudes and latitudes of an estimated track with a length of n as $\left\{\left({lon}_i,{lat}_i\right)\mid i=1,2,\dots ,n\right\}$, the track can be smoothed by a moving average filter with a window length of m as

$$\left\{\begin{array}{c}lo{n}_i^{\prime }=\frac{1}{m}\sum _{k=i-m/2}^{i+m/2}lo{n}_k\hfill \\ la{t}_i^{\prime }=\frac{1}{m}\sum _{k=i-m/2}^{i+m/2}la{t}_k\hfill \end{array}\right..$$

(14)

Then a smooth track can be generated. It should be noted that the main objective of this article is not target detection and tracking, but a target tracklet association method that is directly applied to the tracklets provided. Target detection and tracking are dependent on signal-to-noise ratio (SNR), the signal-to-clutter ratio (SCR), etc. The effects of these factors can be mitigated in target tracking, which involves state filtering and smoothing and is reflected in target parameter measurement errors. Therefore, investigating the influence of target parameter measurement errors on tracklet associations are more meaningful.

Due to the aforementioned shortcomings of target detection using compact HFSWR, the obtained tracks usually fluctuate and deviate from their true positions and are even fragmented into several short track segments. To improve the continuity of target tracking, several tracklets belonging to each (same) target should be associated and connected. Therefore, once a track is initiated, it is necessary to determine whether it comes from a new target or is a continuation of an existing target track. Two types of tracklets are defined as follows.

(1) Terminated tracklet. It represents a track that meets the termination condition [23] and stops updating its state. The terminated tracklet set $\mathbf{T}\_\mathbf{old}$ is defined as

$$\mathbf{T}\_\mathbf{old}=\left\{\mathbf{T}\_\mathbf{old}\right(j\left)\right|j=1,2,\dots ,N\},$$

(15)

where N denotes the number of terminated tracklets, $\mathbf{T}\_\mathbf{old}\left(j\right)=\{P_{old}\left(1\right),P_{old}\left(2\right),\dots ,P_{old}\left(n\right)\}$ is the jth terminated tracklet that contains n plots. It should be noted that a terminated tracklet could be a fully completed track or a portion of a track that is interrupted.

(2) Initiated tracklet. It represents a new track that satisfies the track initiation condition [22] and is defined as

$$\mathbf{T}\_\mathbf{new}=\left\{\mathbf{T}\_\mathbf{new}\right(i\left)\right|i=1,2,\dots ,M\},$$

(16)

where M is the number of newly initiated tracklets, $\mathbf{T}\_\mathbf{new}\left(i\right)=\{P_{new}\left(1\right),P_{new}\left(2\right),\dots ,P_{new}\left(l\right)\}$ denote the ith initiated tracklet with a length of l. It should be noted that an initiated tracklet could be an independent new track or a track portion that can be associated with an existing terminated tracklet.

2.2. Average Heading and Average Speed Calculation

Heading and speed are two important motion characteristics of a moving target. Compact HFSWR can only provide a coarse azimuth resolution; thus, the positions of the measured plots may deviate from their true values. Therefore, the instantaneous heading and speed cannot be accurately obtained using adjacent target positions. Fortunately, average heading and speed can be robustly estimated and reflect the overall motion characteristics of a moving target. An illustrative comparison between instantaneous and average headings is shown in Figure 1.

In Figure 1, the instantaneous and average headings are depicted in solid and dot dash lines, respectively, for a terminated tracklet $\mathbf{T}\_\mathbf{old}\left(j\right)$ and its corresponding initiated tracklet $\mathbf{T}\_\mathbf{new}\left(i\right)$. It can be seen that the instantaneous headings at different sampling times change abruptly, while the average headings for $\mathbf{T}\_\mathbf{old}\left(j\right)$ and $\mathbf{T}\_\mathbf{new}\left(i\right)$ are almost the same. Therefore, the average heading is a more stable characteristic; the same for the average speed.

Taking the ith initiated tracklet $\mathbf{T}\_\mathbf{new}\left(i\right)$ with a length of l as an example, the average heading and speed can be calculated as follows.

(1) Average heading.

The instantaneous heading ${\phi }_{T\_new}^i\left(k\right)$ of the tracklet $\mathbf{T}\_\mathbf{new}\left(i\right)$ at time k is defined as

$${\phi }_{T\_new}^i\left(k\right)=arctan(y/x),$$

(17)

where

$$\begin{array}{cc}\hfill & x=sin({lon}_{T\_new}^i\left(k\right)-{lon}_{T\_new}^i(k-1))cos\left({lat}_{T\_new}^i\left(k\right)\right),\hfill \\ \hfill & y=cos\left({lat}_{T\_new}^i(k-1)\right)sin\left({lat}_{T\_new}^i\left(k\right)\right)\hfill \\ \hfill & -sin\left({lat}_{T\_new}^i(k-1)\right)cos\left({lat}_{T\_new}^i\left(k\right)\right)cos({lon}_{T\_new}^i\left(k\right)-{lon}_{T\_new}^i(k-1)),\hfill \end{array}$$

(18)

$L_{T\_new}\left(k\right)=\left({lon}_{T\_new}^i\left(k\right),{lat}_{T\_new}^i\left(k\right)\right)$ represents the target position of the tracklet $\mathbf{T}\_\mathbf{new}\left(i\right)$ at time k in longitude and latitude, and it is determined by the measured range $r_k$, azimuth ${\theta }_k$ at time k as well as the radar site. Based on the above definition, the average heading of the tracklet $\mathbf{T}\_\mathbf{new}\left(i\right)$ is calculated as

$${\overline{\phi }}_{T\_new}^i=\frac{1}{l-1}\sum _{k=2}^l{\phi }_{T\_new}^i\left(k\right).$$

(19)

In order to further verify the feasibility of using the average heading as the track feature, two track segments were selected for validation, as shown in Figure 2a. These two tracklets can be associated with the same automatic identification system (AIS) track using the track-to-track association method [24], i.e., it is confirmed that they are derived from the same target. The instantaneous headings of the two tracklets were calculated separately using Equation (17) and are shown in Figure 2b, which illustrates that the instantaneous headings of the two tracklets fluctuate severely. In contrast, the average heading of tracklet 1 calculated by Equation (19) is 121.69${}^{\circ }$ and that of tracklet 2 is 120.98${}^{\circ }$, showing that the average headings have better consistencies.

(2) Average speed.

The instantaneous speed $v_{T\_new}^i\left(k\right)$ of the tracklet $\mathbf{T}\_\mathbf{new}\left(i\right)$ at time k is defined as

$$v_{T\_new}^i\left(k\right)=d(L_{T\_new}\left(k\right),L_{T\_new}(k-1))/T,$$

(20)

where $d(L_{T\_new}\left(k\right),L_{T\_new}(k-1))$ is the geodesic distance between adjacent target positions $L_{T\_new}\left(k\right)$ and $L_{T\_new}(k-1)$, and T is the radar sampling interval. Based on the above definition, the average speed of the tracklet $\mathbf{T}\_\mathbf{new}\left(i\right)$ is calculated as

$${\overline{v}}_{T\_new}^i=\frac{1}{l-1}\sum _{k=2}^l{v}_{T\_new}^i\left(k\right).$$

(21)

A motion vector containing the average heading and average speed can be denoted as

$${\mathbf{X}}_{T\_new}^i={\left[{\overline{\phi }}_{T\_new}^i{\overline{v}}_{T\_new}^i\right]}^T.$$

(22)

Similarly, the motion vector of the terminated tracklet $\mathbf{T}\_\mathbf{old}\left(j\right)$ can be represented as

$${\mathbf{X}}_{T\_old}^j={\left[{\overline{\phi }}_{T\_old}^j{\overline{v}}_{T\_old}^j\right]}^T.$$

(23)

3. A Multi-Stage Tracklet Association Method

3.1. Rough Tracklet Association Based on K-Means Clustering

The motion characteristics of different vessels are usually different. To determine the possible tracklet pair set, the clustering method was used first to roughly associate the initiated and terminated tracklets. Clustering is an unsupervised machine learning method that can divide the unlabeled data into several classes. The data are similar within a class but different between classes [25]. The k-means clustering method, due to its simple, effective, and real-time characteristics, is (here) used for a preliminary tracklet association. The motion vectors of the initiated tracklets were selected as the initial clustering centers to reduce the number of tracklet pairs to be associated, and to improve the effectiveness of the association procedure. The procedure for the rough tracklet association using k-means clustering is summarized in Algorithm 1.

Algorithm 1 Rough tracklet association using k-means clustering

Input:  An initiated tracklet set $\mathbf{T}\_\mathbf{new}$ that contains M tracklets, a terminated tracklet set $\mathbf{T}\_\mathbf{old}$ that contains N tracklets. Output:  A preliminary tracklet pair set $\mathbf{T}=\left\{\left[\mathbf{T}\_\mathbf{new}\right(i),\mathbf{T}\_\mathbf{old}(j\left)\right]|i\in [1,M],\right.$$\left.j\in [1,N]\right\}$. 1: Equations (20)–(23) are applied to the initiated tracklets $\mathbf{T}\_\mathbf{new}\left(i\right),i=1,2,\dots ,M$ and terminated tracklets $\mathbf{T}\_\mathbf{old}\left(j\right),j=1,2,\dots ,N$, respectively, to obtain the motion vectors ${\mathbf{X}}_{T\_new}^i$ and ${\mathbf{X}}_{T\_old}^j$. Initialize ${\mathbf{X}}_{T\_new}^i,i=1,2,\dots ,M$ as the clustering centers. 2: The similarities ${\mathbf{D}}_{ij}s$ between the motion vectors of the terminated tracklet and the clustering centers are calculated using Equation (24), and a terminated tracklet is assigned to the cluster with the highest similarity. $${\mathbf{D}}_{ij}=\sqrt{{\Delta }_{ij}^T{\mathbf{Q}}^{-1}{\Delta }_{ij}},$$ (24) where ${\Delta }_{ij}={\mathbf{X}}_{T\_new}^i-{\mathbf{X}}_{T\_old}^j$, $\mathbf{Q}$ represents the covariance matrix of average heading and speed estimation error. 3: Equation (25) is used to recalculate the mean values of the motion vectors of the terminated tracklets in each cluster to obtain the updated clustering centers ${\mu }_h$ as $${\mu }_h=\frac{1}{K_h}\sum _{T\_old\left(j\right)\in {\mathbf{C}}^h}{\mathbf{X}}_{T\_old}^j,h=1,2,\dots ,M,$$ (25) where ${\mathbf{C}}^h$ denotes the $h^{th}$ cluster set and $K_h$ is its cardinality. 4: The loss value is calculated as $$\mathbf{J}({\mu }_1,\dots ,{\mu }_M)=\sum _{h=1}^M\sum _{{\mathbf{X}}_{T\_old}^j\in {\mathbf{C}}^i}{({\mathbf{X}}_{T\_old}^j-{\mu }_h)}^2.$$ (26) 5: Repeat steps 2–4 until the loss values or the positions of the cluster centers do not change any more. Then a preliminarily associated tracklet pair set $\mathbf{T}=\left\{\left[\mathbf{T}\_\mathbf{new}\right(i),\mathbf{T}\_\mathbf{old}(j\left)\right]|i\in [1,M],j\in [1,N]\right\}$ is obtained. It should be noted that the values of i and j may not take all the values in $[1,M]$ and $[1,N]$, respectively.

3.2. Tracklet Pair Set Refinement by Spatiotemporal Constraints

There is a chronological relationship between the initiated tracklet and terminated tracklet of the same target, i.e., the time of the first plot of the initiated tracklet always lags behind that of the last plot of the terminated tracklet. In addition, the distance between the last plot of the terminated tracklet and the first plot of the initiated tracklet is constrained by the target velocity and time gap between two tracklets. Therefore, the preliminary tracklet pair set produced via the rough tracklet association procedure may be refined using the temporal and spatial constraints. Assume the terminated tracklet $\mathbf{T}\_\mathbf{old}\left(j\right)$ and initiated tracklet $\mathbf{T}\_\mathbf{new}\left(i\right)$ belong to one target track, the temporal constraint can be formulated as

$$t_{T\_new\left(i\right)}^{start}-t_{T\_old\left(j\right)}^{end}\ge 0,$$

(27)

where $t_{T\_old\left(j\right)}^{end}$ denotes the time of the last plot of the terminated tracklet $\mathbf{T}\_\mathbf{old}\left(j\right)$, while $t_{T\_new\left(i\right)}^{start}$ represents the time of the first plot of the initiated tracklet $\mathbf{T}\_\mathbf{new}\left(i\right)$.

The spatial constraint is formulated as

$$|d_1-d_2|\le d_{max},$$

(28)

where $d_1$ denotes the target traveling distance at an average speed $\frac{({\overline{v}}_{T\_new}^i+{\overline{v}}_{T\_old}^j)}{2}$ within the time interval $\left[t_{T\_old\left(j\right)}^{end},t_{T\_new\left(i\right)}^{start}\right]$, $d_2$ is the measured distance between the first plot of the initiated tracklet and the last plot of the terminated tracklet, $d_{max}$ represents the distance threshold.

The tracklet pairs in the preliminary tracklet pair set can be further screened using the constraints in Equations (27) and (28).

3.3. Optimal Tracklet Assignment Based on a Bidirectional Prediction

In conventional tracklet association methods, the target states are backward predicted from the last plot in an initiated tracklet to the instant of the last plot of the terminated tracklets for association. However, due to the low positioning accuracy of the compact HFSWR, the predicted states may deviate from the true trajectory after a few predictions, which may degrade the association performance. To address this problem, a bidirectional prediction method involving initiated tracklet backward prediction and terminated tracklet forward prediction is proposed.

Firstly, an intermediate common time $t^{ij}=\frac{t_{T\_new\left(i\right)}^{start}-t_{T\_old\left(j\right)}^{end}}{2}$ is calculated. Then, the initiated tracklet $\mathbf{T}\_\mathbf{new}\left(i\right)$ is reversely filtered from $t_{T\_new\left(i\right)}^{end}$ to $t^{ij}$ by the Kalman filter to obtain the predicted state ${\mathbf{S}}_{T\_new\left(i\right)}^p\left(t^{ij}\right)$, which can be expressed as

$${\mathbf{S}}_{T\_new\left(i\right)}^p\left(t^{ij}\right)={\left[v_{T\_new\left(i\right)}^p\left(t^{ij}\right)r_{T\_new\left(i\right)}^p\left(t^{ij}\right){\theta }_{T\_new\left(i\right)}^p\left(t^{ij}\right)\right]}^T.$$

(29)

Meanwhile, the terminated tracklet $\mathbf{T}\_\mathbf{old}\left(j\right)$ is forwardly filtered from $t_{T\_old\left(j\right)}^{start}$ to $t^{ij}$ by the Kalman filter to obtain the predicted state ${\mathbf{S}}_{T\_old\left(j\right)}^p\left(t^{ij}\right)$, which can be expressed as

$${\mathbf{S}}_{T\_old\left(j\right)}^p\left(t^{ij}\right)={\left[v_{T\_old\left(j\right)}^p\left(t^{ij}\right)r_{T\_old\left(j\right)}^p\left(t^{ij}\right){\theta }_{T\_old\left(j\right)}^p\left(t^{ij}\right)\right]}^T.$$

(30)

3.3.1. Tracklet Association Cost Calculation

To achieve optimal assignments of tracklets, a tracklet association cost function is defined as

$$cost=1-({sim}_v+{sim}_r+{sim}_{\theta }),$$

(31)

where ${sim}_v$, ${sim}_r$, and ${sim}_{\theta }$ represent the similarities between predicted states ${\mathbf{S}}_{T\_new\left(i\right)}^p\left(t^{ij}\right)$ and ${\mathbf{S}}_{T\_old\left(j\right)}^p\left(t^{ij}\right)$ in the Doppler velocity, range, and azimuth, respectively, and are formulated as

$${sim}_v=W_v\times exp(-{|v_{T\_old\left(j\right)}^p\left(t^{ij}\right)-v_{T\_new\left(i\right)}^p\left(t^{ij}\right)|}^2/{\sigma }_v^2),$$

(32)

$${sim}_r=W_r\times exp(-{|r_{T\_old\left(j\right)}^p\left(t^{ij}\right)-r_{T\_new\left(i\right)}^p\left(t^{ij}\right)|}^2/{\sigma }_r^2),$$

(33)

$${sim}_{\theta }=W_{\theta }\times exp(-{|{\theta }_{T\_old\left(j\right)}^p\left(t^{ij}\right)-{\theta }_{T\_new\left(i\right)}^p\left(t^{ij}\right)|}^2/{\sigma }_{\theta }^2),$$

(34)

where ${\sigma }_v^2$, ${\sigma }_r^2$, and ${\sigma }_{\theta }^2$ denote the variances of Doppler velocity, range, and azimuth, respectively. $W_v$, $W_r$, and $W_{\theta }$ represent the corresponding weights of three kinematic parameters and are set according to the radar measurement accuracy to satisfy

$$W_v+W_r+W_{\theta }=1.$$

(35)

The larger the values of ${sim}_v$, ${sim}_r$, and ${sim}_{\theta }$ are, the lower the association cost is, and the higher the probability that both tracklets come from the same target.

3.3.2. Optimal Tracklet Assignment

After the bidirectional state prediction and tracklet association cost calculations for the tracklet pairs in the refined tracklet pair set, one or more terminated tracklets may satisfy the association conditions with an initiated tracklet. Then an optimal assignment method is employed to achieve optimal allocation among possible tracklet pairs. The commonly used optimal assignment methods include the Munkres algorithm, auction algorithm, and JVC algorithm [26]. As the JVC algorithm has a better balance between performance and efficiency [27], it is used in this article for tracklet assignment. It should be pointed out that the word “optimal” here means the method with global optimal assignment.

Suppose there are m initiated tracklets and n terminated tracklets in the refined tracklet pair set. After bidirectional prediction, the association costs between the predicted states are calculated using Equation (31); an association cost matrix $\mathbf{D}$ can be obtained as

$$\mathbf{D}={\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m2}& \cdots & a_{mn}\end{array}\right]}_{m\times n},$$

(36)

where $a_{ij}$ is the association cost between the predicted states ${\mathbf{S}}_{T\_new\left(i\right)}^p\left(t^{ij}\right)$ and ${\mathbf{S}}_{T\_old\left(j\right)}^p\left(t^{ij}\right)$. The total association cost is calculated as

$$\lambda =\sum _{i=1}^m\sum _{j=1}^n{a}_{ij}{\delta }_{ij},$$

(37)

subject to

$$\sum _{i=1}^m{\delta }_{ij}=1,j=1,2,\dots ,n,$$

(38)

$$\sum _{j=1}^n{\delta }_{ij}=1,i=1,2,\dots ,m,$$

(39)

where ${\delta }_{ij}$ is a binary function and defined as

$${\delta }_{ij}=\left\{\begin{array}{c}1,if\mathbf{T}\_\mathbf{new}\left(i\right)isassociatedwith\mathbf{T}\_\mathbf{old}\left(j\right)\hfill \\ 0,if\mathbf{T}\_\mathbf{new}\left(i\right)isnotassociatedwith\mathbf{T}\_\mathbf{old}\left(j\right).\hfill \end{array}\right.$$

(40)

Equation (38) indicates that any initiated tracklet $\mathbf{T}\_\mathbf{new}\left(i\right)$ can only be associated with, at most, one terminated tracklet $\mathbf{T}\_\mathbf{old}\left(j\right)$, while Equation (39) indicates that any terminated tracklet $\mathbf{T}\_\mathbf{old}\left(j\right)$ can only be associated with, at most, one initiated tracklet $\mathbf{T}\_\mathbf{new}\left(i\right)$.

The best association cost combination $a_{ij}s$ is determined by minimizing the total association cost using the JVC algorithm and is formularized as

$$\underset{a_{ij}}{argmin}\left(\lambda \right).$$

(41)

The final tracklet pair set ${\mathbf{T}}_{final}=\left\{\left[\mathbf{T}\_\mathbf{new}\right(i),\mathbf{T}\_\mathbf{old}(j\left)\right]|i\in [1,m],j\in [1,n]\right\}$ is obtained according to the indices i and j of $a_{ij}s$.

The flowchart of the proposed multi-stage tracklet association algorithm is shown in Figure 3.

4. Experiment Results

To evaluate the association performance, tracklet association tests with both simulated and field data were conducted using the proposed multi-stage tracklet association method and the results are compared with that of the TSA method presented in [12].

4.1. Experiments with Simulated Data

Plot data of five different targets were simulated with initial kinematic parameters listed in

Table 1. Initial parameters of simulated targets.

	Initial Range (km)	Initial Azimuth ()	Initial Doppler Velocity (km/h)
Target 1	107	−10	22
Target 2	106	−11	21
Target 3	74	−7	−19
Target 4	112	−8	−23
Target 5	93	−7	−19

. Negative azimuth values indicate they are on the left side of the radar boresight, while negative Doppler velocities mean the targets move away from the radar. Some frames of data were removed intentionally to simulate missed detections.

According to the statistics reported in [22], the measurement errors in range, azimuth, and radial velocity of the developed radar are usually less than 4 km, 3, and 1 km/h, respectively. Thus, Gaussian white noise with a mean of zero and standard deviations of 4 km, 3, and 1 km/h were added to the target range, azimuth, and radial velocity data sequences, respectively, to generate simulated target data. The CMKF [28] was applied to the simulated data to produce target tracks; the obtained tracking results are shown in Figure 4. The red dots denote the starting positions of the tracklets, the tracklets of the same target are depicted using the same color and marker. It should be noted that target 1, target 2, and target 3 have two tracklets each, while target 4 and target 5 only have one tracklet. The kinematic parameters of target 1 and target 2, as well as target 5 and target 3 are set similar to challenge the tracklet association methods.

Both the proposed multi-stage tracklet association method and the track segment association method in [12] were applied to these simulated data and the results are shown in Figure 5. The associated tracklet pairs are depicted in the same color. Since target 4 and target 5 have only one tracklet and both tracklets are not associated with any other tracklet, the tracklets of target 4 and target 5 are not shown and analyzed.

It can be seen from the association results in Figure 5a that the TSA method in [12] correctly associated the tracklets of target 3, but failed for the tracklets of target 1 and target 2. The proposed multi-stage tracklet association method achieves successful tracklet association for all three targets. The association costs between different tracklet combinations of target 1 and target 2 were calculated by Equation (17) and the results are listed in

Table 2. Association costs between different tracklet combinations of target 1 and target 2.

	Terminated Tracklet	Terminated Tracklet 1	Terminated Tracklet 2
Initiated Tracklet
Initiated Tracklet 1		0.1117	0.0228
Initiated Tracklet 2		0.7020	0.5933

.

According to the association costs listed in Table 2, the TSA method in [12] will associate the initiated tracklet of target 1 with the terminated tracklet of target 2 as they have smaller association costs (0.0228); then the initiated tracklet of target 2 can only be connected to the terminated tracklet of target 1, i.e., false tracklet association occurs. On the contrary, the proposed multi-stage tracklet association method conducts tracklet association by minimizing the total association cost. In this case, the total association cost reaches the minimum value of 0.7050 = 0.1117 + 0.5933, and it corresponds to the correct tracklet pairs of target 1 and target 2.

4.2. Experiments with Field Data

To demonstrate the effectiveness of the proposed method in a real scenario, target tracking followed by tracklet association experiments were conducted using the data collected by a newly developed compact HFSWR–CORMS, which was deployed at the shore of the North China Sea on 18 January, 2019. The radar system uses an eight-element receiving antenna array with an antenna aperture of 105 m, an operating frequency of 4.7 MHz, and a data rate of 1 frame/min. A total of 266 frames of data were collected from 11:04 a.m. to 3:29 p.m.

Firstly, the CMKF method was applied to produce the target tracks, then the broken tracklets from four different targets were selected for tracklet association tests, marked as terminated tracklets 1–4 and initiated tracklets 1–4 in Figure 6, where the red dots indicate the starting plots.

4.2.1. Analysis of Track Fragmentation Cause

In order to better apply the tracklet association method, the causes of track fragmentation were analyzed using the tracklets from target 1 and target 2, as shown in Figure 7. The track-to-track association method in [24] was used to associate these tracklets with corresponding AIS tracks to confirm that the two tracklets in Figure 7a were from one target and those in Figure 7b were from another target. The measured plots within the breakage period are depicted in blue.

(1) Missed detection.

As can be seen from Figure 7a, terminated tracklet 1 ended at 12:24 and started the initiated tracklet 1 at 12:34. During 12:25–12:33, there should have been nine frames of data, but only six frames of data were acquired. The missed detections may result in track fragmentation.

(2) Clutter interference.

In Figure 7b, terminated tracklet 2 ended at 1:59 p.m. and started the initiated tracklet 2 at 2:04 p.m. During 2:00–2:03 p.m., the measured plots are scattered without obvious trend. By analyzing the track data of target 2, it was found that the target Doppler velocities during the breakage period were similar to that of the Bragg wave, i.e., the target echoes may be masked by the first-order sea clutter and cannot be detected. Thus, the measured plots do not form a smooth interpolation. The false alarms of the high-frequency surface wave radar were mainly caused by various clutters, such as sea clutter, ionospheric clutter, radio frequency interference, etc. The track, affected by the false plots, cannot associate the correct target plots; thus, track fragmentation occurs.

The above analysis is consistent with the discussion in Section 1. In addition to the above causes, a low data rate and low detection accuracy may also lead to track fragmentation.

4.2.2. Analysis of Tracklet Association Results

The selected tracklets in Figure 6 were used for tracklet association tests. Firstly, the preliminary tracklet association results using k-means clustering are shown in Figure 8a, and the preliminary tracklets that are classified into one cluster are depicted using the same color and marker. The results show that the initiated tracklets 1 and 2 and the terminated tracklets 1 and 2 are classified into the same cluster, while the initiated tracklets 3 and 4 and terminated tracklets 3 and 4 are placed into the same category.

Then the spatial and temporal constraints were applied to the preliminary associated tracklets to produce the candidate tracklet pairs, as shown in Figure 8b. After this step, the initiated tracklet 1 is associated with the terminated tracklet 1, while the initiated tracklet 2 is associated with the terminated tracklet 2; they are successfully distinguished and depicted in green and blue, respectively. However, the association results of the initiated tracklets 3 and 4 and terminated tracklets 3 and 4 remain the same.

Subsequently, bidirectional prediction was implemented on the candidate tracklets to produce the predicted states, the obtained results are shown as red lines in Figure 9, and the prediction results for the fragmentation parts are zoomed in for a better view. It can be seen that the initiated tracklet 1 and terminated tracklet 1, as well as the initiated tracklet 2 and terminated tracklet 2, can be easily associated, correspondingly. However, association ambiguities still exist in the initiated tracklets 3 and 4 and terminated tracklets 3 and 4 due to their similar headings and speeds.

In order to resolve the association ambiguity of the initiated tracklets 3 and 4 and terminated tracklets 3 and 4, the association costs between possible tracklet pairs were calculated using Equation (31) and listed in

Table 3. Association costs between tracklet pairs of target 3 and target 4.

	Terminated Tracklet	Terminated Tracklet 3	Terminated Tracklet 4
Initiated Tracklet
Initiated Tracklet 3		0.6052	0.1577
Initiated Tracklet 4		0.7173	0.1124

.

Finally, the optimal assignment method based on the JVC algorithm was applied to the association cost matrix of the candidate tracklet pairs to produce the final tracklet association results, as shown in Figure 10. The associated tracklet pairs for the same target are shown in the same colors. It is shown that the tracklets are correctly matched for both target 3 and target 4 with a minimum total association cost of 0.7176 = 0.6052 + 0.1124.

4.2.3. Analysis of Association Accuracy and Computational Complexity

To test the computational efficiency and association accuracy of the proposed method, the correct association rate $R_t$, false association rate $R_f$, and missing association rate $R_n$ are defined in Equations (42)–(44), respectively.

$$R_t=\frac{n_t}{num}$$

(42)

$$R_f=\frac{n_f}{num}$$

(43)

$$R_n=\frac{n_n}{num}$$

(44)

where $num$ is the total number of tracklets, $n_t$, $n_f$, and $n_n$ denote the numbers of correctly associated tracklets, incorrectly associated tracklets, and missed associated tracklets, respectively.

The converted measurement Kalman filter was applied to the field HFSWR data introduced in Section 4.2, first to produce target tracklets, and 51 pairs of tracklets were selected for tracklet association tests. It was verified that each tracklet pair was from the same target by associating them with the same AIS track using the track-to-track association method in [24]. Then, both the proposed multi-stage tracklet association method and the track segment association method in [12] were applied to these 51 tracklet pairs to produce the tracklet association results. The correct association rate $R_t$, false association rate $R_f$, and missing association rate $R_n$ were calculated and the obtained results are listed in

Table 4. Performance comparisons between two tracklet association methods.

	${\mathit{R}}_{\mathit{t}}$ (%)	${\mathit{R}}_{\mathit{f}}$ (%)	${\mathit{R}}_{\mathit{n}}$ (%)	Average Running Time (s)
TSA method in [12]	63.4	19.5	17.1	10.4
Proposed Method	93.5	4.3	2.2	3.3

. Moreover, in order to evaluate the computational efficiency, the proposed method and the track segment association method in [12] were run 200 times, and the average running time was calculated, respectively; the corresponding results are also listed in Table 4.

It can be seen from the results in Table 4 that, compared with the TSA method in [12], the correct association rate of the proposed method is improved by 30.1%, while its false association rate and missing association rate are reduced by 15.2% and 14.9%, respectively. The average running time of the proposed method is 7.1 s less than that of the TSA method in [12].

5. Discussion

From the above analysis and experimental results, it can be summarized that:

(i) Track fragmentation often occurs during target tracking with compact HFSWR due to its physical limits. Tracklet association methods can connect the broken track segments belonging to a same target, enhancing the tracking continuity.

(ii) The proposed multi-stage tracklet association method provides a “coarse-to-fine” way to match the same target’s tracklets, both accurately and efficiently, which has been verified by experimental results using both simulated and field data with respect to the correct association rate, false association rate, missing association rate, and average running time.

(iii) The rough tracklet association using k-means clustering and tracklet pair refinement by spatiotemporal constraint help to find the possible tracklet pairs. As the tracklets that cannot be associated with any other tracklet have been excluded for further analysis, on the one hand, the computational burden is reduced; on the other hand, the association accuracy can be improved.

(iv) Taking the low positioning accuracy of compact HFSWR into consideration, bidirectional prediction and a global optimal assignment based on JVC resolve the association ambiguities for close tracklets and improve the association accuracy.

6. Conclusions

Track fragmentation is a common problem for target tracking with a compact high-frequency surface wave radar since it degrades the tracking continuity and maritime surveillance performance. In this article, the tracklet association problem for compact HFSWR was investigated and a multi-stage tracklet association method, which consists of rough association using k-means clustering, refinement by spatiotemporal constraint, bidirectional prediction followed by optimal assignment based on JVC algorithm, was developed to repair the broken tracklets. The advantages of the proposed two-stage method are two-fold. In the first stage, k-means clustering can find the possible tracklet pairs efficiently, while the spatiotemporal constraint excludes the tracklet pairs that cannot be associated; thus, the computational efficiency for tracklet association is significantly improved. In the second stage, bidirectional state prediction and the optimal assignment based on the JVC algorithm are able to enhance the tracklet association accuracy under low positioning accuracy and multi-target tracking scenarios. Experimental results with both simulated and field data demonstrate that the proposed method is effective and able to resolve association ambiguities for similar tracklets and improve the track integrity; the association performance is superior to that of the existing track segment association method with a correct association rate improved by 30%.

Low detection probability and low spatial resolution are two main factors leading to track fragmentation for the compact high-frequency surface wave radar. On the one hand, more advanced tracklet association methods should be developed with detection characteristics of a compact high-frequency surface wave radar fully considered. On the other hand, a networking observation with multiple radars is a potential way to improve the target detection probability and positioning accuracy. Moreover, the proposed method should be extended to be applicable to the case involving maneuvering targets.