2.1. Track tree formation

TOMHT uses the structure of track trees to manage targets. A track tree corresponds to a target, each node on the track tree corresponds to a measurement (including dummy measurements) and each path from the root node to the leaf node corresponds to a track hypothesis, of which only one track corresponds to the real track of the target. At each frame, a new track tree is constructed for each measurement to account for a new target. The previously existing track is extended by spawning a separate branch for each measurement within its gating area. If no measurements are associated with the track, an empty node is spawned.

2.2. Track score

The track score is used to evaluate each track in MHT, and it can be defined as the cumulative Log-Likelihood Ratio (LLR). According to Bayesian theory, the score of the j-th track at frame j can be expressed in the following recursive form,

(1)$$S^j(k)=S^j(k-1)+\Delta S^j(k)$$

where ΔS ^j(k) is the increment of the track score, which is calculated as

(2)$$\Delta S^j(k) = \left\{ {\matrix{ {\ln \left( {\displaystyle{{P_{\rm D}} \over {{(2\pi )}^{q/2}(\rho _{\rm n} + \rho _{\rm f})\sqrt { \vert \sum\limits_k^j \vert } }}} \right) - \displaystyle{{d_{\rm M}^2 } \over 2},} \hfill & {i > 0} \hfill \cr {\ln {\mkern 1mu} (1 - P_{\rm D}),} \hfill & {i = 0} \hfill \cr } } \right.$$

where P _D is the probability of detection, ρ_n and ρ_f are the spatial density of new targets and clutter, respectively, $\sum^{j}_{k}$ is the covariance matrix which is estimated by a Kalman filter for the j-th track at frame k, q is the state dimension and d _M is the Mahalanobis distance between an observation location and a predicted location. i > 0 means the j-th track is associated with the i-th measurement at frame k while i = 0 means the j-th track is associated with the dummy measurement. Each measurement can initialise a new track, and the initial score of each track is S ^j(1) = ρ _n/ρ _f.

2.3. Global hypothesis formation and pruning

The global hypothesis formation is equivalent to the Multidimensional Assignment (MDA) or the Maximum Weight Independent Set (MWIS) problem (Papageorgiou and Salpukas, Reference Papageorgiou and Salpukas2009) and the MWIS problem is relatively easy to solve. An undirected graph G(V, B) is constructed by assigning each track to a graph vertex g _j, g _j ∈ V. A weight w _j of the vertex g _j corresponds to its track score S ^j(k). An edge (i, j) ∈ B connects two vertices g _i and g _j if the two tracks share common observations. An independent set is a set of vertices with no edges in common. Hence, finding MWIS is equivalent to finding the set of compatible tracks that maximises the total track score. This can be expressed as an integer linear programming problem,

(3)$$\eqalign{\max_{g_i} &\quad \sum_i g_iw_i \cr \hbox{s.t.}&\quad g_i+g_j\le 1,\ \forall(i,j)\in B\quad g_i\in \{0,1\},\ \forall i }$$

where s.t. is the abbreviation of “subject to”. Because of the exponential increase in the number of track hypotheses over time (Kim et al., Reference Kim, Li, Ciptadi and Rehg2015), it is necessary to apply some strategies to reduce computational complexity. There are many pruning mechanisms in a classical MHT implementation. In this paper, we use N-scan pruning and global pruning based on the posterior probability of each track. The N-scan pruning makes the data association decisions at frame k − N + 1 using old measurements from frame k − N + 1 to k. After N-scan pruning is completed, global track deletion is performed based on the track probability. All tracks with a probability lower than a threshold will be deleted. For more details about pruning mechanisms, see Blackman and Popoli (Reference Blackman and Popoli1999).

3.1. Motion model

Ships' tracks mainly include great-circle and rhumb-line navigation tracks (Sölver and Marcus, Reference Sölver and Marcus1958). The former is the shortest track, but in reality, piecewise rhumb-line tracks are normally used to approximate the great-circle track due to its frequent changes of course. Therefore, the latter is studied in this paper. Latitude and longitude are usually used to describe position and can be transformed into plane coordinates by projections, such as the Transverse Mercator (TM) projection and the Universal TM (UTM) projection. The revisit interval of an electronic intelligence satellite network is so long that the target state has a deformation in the projection plane, which affects the precision of prediction of the target state. Therefore, we describe the target state in geographical coordinates, that are defined in frame k as

(4)$${\bi x}_k=(\lambda_k, \dot{\lambda}_k, \varphi_k,\dot{\varphi}_k)^{\rm T}$$

where λ_k, φ_k and $\dot{\lambda}_{k}$, $\dot{\varphi}_{k}$ represent the position and velocity components along longitude and latitude directions, respectively. The basic unit of latitude, longitude, and course is the degree (°). The units of distance, time and velocity are nautical miles (nm), hours (h), knots (kn, kn = nm/h), respectively. The locations in the northern and southern hemisphere are denoted as positive and negative latitudes, respectively, while locations in the eastern and western hemisphere are denoted as positive and negative longitudes, correspondingly. Given the speed v _k and the course θ_k over ground (relative to true north), it can be obtained that

(5)$$\left\{ {\matrix{ {{\dot \lambda }_k = v_k\sin \theta _k} \cr {{\dot \varphi }_k = v_k\cos \theta _k} \cr } } \right.$$

Middle latitude sailing is simpler and more convenient to use in the dead reckoning of a rhumb-line track (Sölver and Marcus, Reference Sölver and Marcus1958), so we use it for track prediction. The state transition model from one point (λ₁, φ₁) to the other point (λ₂,φ₂) is

(6)$$\left\{ \matrix{\varphi _2 = \varphi _1 + \dot \varphi T/a \hfill \cr \lambda _2 = \lambda _1 + \dot \lambda T{\mkern 1mu} {\rm sec}(\varphi ^{\rm m})/a \hfill} \right.$$

where T is the time interval, a is the conversion factor which equals 60 nm/° since 1° of longitude along the equator is approximately 60 nm. φ^m is the middle latitude at which the arc length of the parallel separating the meridians passing through two specific points is exactly equal to the departure in proceeding from one point to the other. φ^m is approximately the mean latitude in sailing (Bowditch, Reference Bowditch2011),

(7)$$\varphi^{\rm m}=(\varphi_1+\varphi_2)/2=\varphi_1+\dot{\varphi}T/(2a)$$

As can be seen from Equation (6), λ₂ is prone to larger errors at high latitudes than at low latitudes when the errors of φ^m are the same due to the secant function. λ₂ tends to infinity when $\varphi^{\rm m}=\pm90^{\circ}$, so Equation (6) cannot be used accurately in the high north or the far south. If we want to track ships at high latitudes (latitude higher than 80°), we should use other coordinate systems, such as universal polar stereographic coordinates which cover the Earth's polar regions (Hager et al., Reference Hager, Behensky and Drew1989). In this paper, the region in which we are interested is at low and middle latitudes. The motion model can be described as

(8)$${\bi x}_k=f_{k-1}({\bi x}_{k-1})+{\bi G}_{k-1}{\bi v}_{k-1}$$

where f _k−1( · ) is the state transition function, and G_k−1 and v_k−1 are the process noise matrix and process noise, respectively. We define ${\bi v}_{k-1}=\lpar \ddot{\lambda}_{k-1}\comma \; \ddot{\varphi}_{k-1}\rpar ^{\rm T}$ and v_k−1 represents a multivariate Gaussian with zero mean. The standard deviations of $\ddot{\lambda}_{k-1}$ and $\ddot{\varphi}_{k-1}$ are both set to σ_v. Derived from Equation (6), the state transition equation is

(9)$$\matrix{ {x_k} \hfill & { = \left[ {\matrix{ {\lambda _k} \cr {{\dot \lambda }_k} \cr {\varphi _k} \cr {{\dot \varphi }_k} \cr } } \right] = \left[ {\matrix{ {\lambda _{k - 1} + {\dot \lambda }_{k - 1}T_{k - 1}{\rm sec}(\varphi _{k - 1}^{\rm m} )/a} \cr {{\dot \lambda }_{k - 1}} \cr {\varphi _{k - 1} + {\dot \varphi }_{k - 1}T_{k - 1}/a} \cr {{\dot \varphi }_{k - 1}} \cr } } \right]} \hfill \cr {} \hfill & {\quad + \left[ {\matrix{ {{(T_{k - 1})}^2{\mkern 1mu} {\rm sec}(\varphi _{k - 1}^{\rm m} )/(2a)} & 0 \cr {T_{k - 1}} & 0 \cr 0 & {{(T_{k - 1})}^2/(2a)} \cr 0 & {T_{k - 1}} \cr } } \right]\left[ {\matrix{ {{\ddot \lambda }_{k - 1}} \cr {{\ddot \varphi }_{k - 1}} \cr } } \right]} \hfill \cr } $$

with T _k−1 = t _k − t _k−1, $\varphi^{\rm m}_{k-1}=\varphi_{k-1}+\dot{\varphi}_{k-1}T_{k-1}/\lpar 2a\rpar $, where t _k−1 and $\varphi^{\rm m}_{k-1}$ are the time stamp and the middle latitude at frame k − 1, respectively.

3.2. Measurement model

By means of geolocation processing (multi-satellite time difference of arrival technique, etc.) in the short receiving time available, we can get the longitude and latitude of an emitter, which is used as the position measurement z_k. The measurement equation for position information at frame k is

(10)$${\bi z}_k={\bi H}_k{\bi x}_k+{\bi w}_k,{\bi w}_k\sim N(\bf{0},{\bi R}_k)$$

(11)$${\bi H}_k = \left[ {\matrix{ 1 & 0 & 0 & 0 \cr 0 & 0 & 1 & 0 \cr } } \right]$$

(12)$${\bi R}_k = \left[ {\matrix{ {\sigma _\lambda ^2 } & 0 \cr 0 & {\sigma _\varphi ^2 } \cr } } \right]$$

where H_k is the observation matrix, w_k is the position noise and R_k is the error covariance matrix which dominates the spread of measurements (spatial extension or shape). σ _λ and σ _φ are the standard deviations of positioning in longitude and latitude, respectively. Different circles of latitude have different lengths, and the Equator is the longest. To correspond to Equation (6), we set

(13)$$\sigma_\varphi=\sigma_{\rm p},\sigma_\lambda=\alpha\sigma_{\rm p},\alpha\ge 1$$

where σ_p is the standard deviation of positioning at the Equator, α is the weight in longitude and it should increase with an increase of latitude. To simplify processing, α is set to a constant in the surveillance region. Since an emitter is regarded as an extended target, we assume that the number of measurements from each emitter obtained by the satellite receiver is in a Poisson distribution.

An electronic intelligence satellite can not only acquire position information of emitters, but also obtain rich feature information such as Centre Frequency (CF), Pulse Repetition Interval (PRI) and Pulse Width (PW) by the processing of emitter signals. The CF is the carrier frequency of the electromagnetic wave and it is more able to reflect the characteristics of an emitter compared with other features.

Therefore, in order to simplify the feature model, only the CF of an emitter is used as the feature information in this paper. In reality, the CF often changes for anti-jamming and anti-reconnaissance, and the changes of CF of different emitters vary. For example, the CF range of a navigation radar is much smaller than that of an early-warning radar, and the frequency of some phased array radars can range over a few hundred Mega-Hertz (MHz). We assume that each ship is equipped with one or more emitters, and the CF of each emitter is a one-dimensional Gaussian distribution with different means and variances. As the measure accuracy of frequency is high (generally less than 1 MHz), the measurement error of the CF can be negligible compared with the frequency range. The measurement equation for frequency information is defined as follows,

(14)$$f_k=e^i_k,e^i_k\sim N(\mu_i,\sigma^2_i)$$

where f _k denotes the frequency measurement, $e^{i}_{k}$ is the true CF of the i-th emitter of a ship and μ_i and σ_i are the mean and the standard deviation of CF, respectively.

In this paper, both position information and feature information are used in a measurement and the augmented measurement is denoted as $\tilde{{\bi z}}_{k}=\lcub {\bi z}_{k}\comma \; f_{k}\rcub $ when satellites revisit the target at frame k. As shown in Figure 2, by simulation, a ship is equipped with three emitters with different means and variances of CF, and each emitter can produce multiple measurements during the short time when a satellite passes over the ship. The positions of measurements (left part of Figure 2) are distributed around the ground truth position of the ship. Points in different colours represent the measurements of different emitters. The centre frequencies of measurements (right part of Figure 2) are distributed around different means.

Figure 2.

A simulated example of measurement distribution of a ship.

4.1. Pre-processing of measurements

In a complex electromagnetic environment, an electronic intelligence satellite can not only receive emitter signals, but also receive clutter resulting in a large number of false positioning measurements. Measurements of an extended target cannot be directly used for tracking. Therefore, a pre-processing of measurements, which mainly includes two clustering operations, is performed to feed the tracking algorithm. The first clustering operation is to cluster the measurements of the same emitter to get a new measurement and reduce clutter, and the second operation is to cluster different emitters from the same ship, so as to track ships instead of tracking emitters and improve the accuracy of tracking.

To solve the clustering problem with the number of targets unknown, we use agglomerative hierarchical clustering (Jain et al., Reference Jain, Murty and Flynn1999) with the location and frequency information. The new measurements obtained after the first and the second clustering are represented as a dot and two dots on the top of measurements, respectively. The process is as follows:

(1) First clustering operation.

Input: Measurement set $\tilde{{\bi Z}}_{k}=\lcub \tilde{{\bi z}}^{i}_{k}\comma \; 1\le i\le N_{k}\rcub $, $\tilde{{\bi z}}^{i}_{k}=\lcub {\bi z}^{i}_{k}\comma \; f^{i}_{k}\rcub $, where N _k indicates the total number of measurements at the k-th revisit time.

Initialisation: Select each measurement as the initial cluster which is defined as $W_{i}=\lcub \tilde{{\bi z}}^{i}_{k}\rcub $. The maximum distance between the clusters is used to measure similarity, that is

(15)$$d(W_i,W_j)=\max_{\tilde{{\bi z}}^m_k\in W_i,\tilde{{\bi z}}^n_k\in W_j} d(\tilde{{\bi z}}_k^m,\tilde{{\bi z}}_k^n)$$

The distance between two measurements is defined as

(16)$$d(\tilde{{\bi z}}_k^i ,\tilde{{\bi z}}_k^j ) = \left\{ {\matrix{ {\displaystyle{{ \vert f_k^i - f_k^j \vert } \over {\Delta f_k}},} \hfill & {{\rm if }d_{\rm M}({\bi z}_k^i ,{\bi z}_k^j ) \le \sqrt {\delta _{{\rm P}_{\rm G}}} } \hfill \cr {1,} \hfill & {{\rm otherwise}} \hfill \cr } } \right.$$

(17)$$d_{\rm M}({\bi z}_k^i ,{\bi z}_k^j ) = \sqrt {{({\bi z}_k^i - {\bi z}_k^j )}^T{\left[ {\matrix{ {2{(\sigma _{\rm p})}^2} & 0 \cr 0 & {2{(\alpha \sigma _{\rm p})}^2} \cr } } \right]}^{ - 1}({\bi z}_k^i - {\bi z}_k^j )} $$

where Δ f _k is the range of frequency measurement at frame k, $\Delta f_{k}=\max\limits_{i} f_{k}^{i}-\min\limits_{i} f_{k}^{i}$ and $d_{\rm M}\lpar {\bi z}_{k}^{i}\comma \; {\bi z}_{k}^{j}\rpar $ denotes the Mahalanobis distance between two measurements. Ignoring slow movement of ships during the time that satellites are overhead, d _M² is a χ² distribution with two degrees of freedom. $\delta_{{\rm P}_{\rm G}}$ is a threshold which can be computed using inverse cumulative χ² distribution for a given probability P _G, where $\delta_{{\rm P}_{\rm G}}=\hbox{invchi}2\lpar P_{\rm G}\rpar $, and invchi2( · ) is the function of inverse cumulative χ² distribution. If $\delta_{{\rm P}_{\rm G}}$ is small, measurements from the same emitter cannot be clustered, which would cause a lot of false targets. If $\delta_{{\rm P}_{\rm G}}$ is too large, measurements from different emitters (relatively close to each other with equal frequencies) may be clustered. However, the impact of wrong clustering caused by a large $\delta_{{\rm P}_{\rm G}}$ on tracking can be reduced by multi-frame data association. We set $\delta_{{\rm P}_{\rm G}}=8$ for $P_{G}\approx98\percnt $ to ensure measurements of the same emitter can belong to the same cluster and reduce false targets.

Loop: Find the nearest two clusters W _m and W _n, and merge them to obtain a new cluster W _t, $W_{t}=W_{m}\cup W_{n}$. When the distance between clusters is greater than a threshold f _T, we stop clustering. Taking into account the frequency range of the same emitter, f _T is set to 200/Δ f, which means frequency measurements within 200 MHz will be regarded as one cluster.

Output: Cluster set {W _i, 1 ≤ i ≤ N _m}, where N _m is the number of clusters. Considering that more than one measurement of an emitter can be obtained, W _i is saved when the number of measurements is greater than one, otherwise it is regarded as clutter to be removed. The new cluster set is defined as $\lcub \dot{W}_{t}\comma \; 1\le i\le \dot{N}_{\rm m}\rcub $, where $\dot{N}_{\rm m}$ is the number of saved clusters. We define a new measurement set as $\dot{\tilde{{\bi Z}}}=\lcub \dot{\tilde{{\bi z}}}^{i}_{k}\comma \; 1\le i\le \dot{N}_{m}\rcub $, where $\dot{\tilde{{\bi Z}}}^{i}_{k}=\lcub \dot{{\bi z}}^{i}_{k}\comma \; F^{i}_{k}\rcub $, $F^{i}_{k}=\lsqb \mu^{i}_{k}\comma \; \sigma^{i}_{k}\comma \; L^{i}_{k}\rsqb $. $\dot{\tilde{{\bi Z}}}^{i}_{k}$ is the mean of position measurements in $\dot{W}_{i}$, $\mu^{i}_{k}$, $\sigma^{i}_{k}$ and $L^{i}_{k}\lpar L^{i}_{k}\ge 2\rpar $ are the mean, the standard deviation and the number of frequency measurements in $\dot{W}_{i}$, respectively. The standard deviations of positioning for a new position measurement $\dot{{\bi z}}^{i}_{k}$ in latitude and longitude are approximately $\sigma_{\rm p}/\sqrt{L^{i}_{k}}$ and $\alpha\sigma_{\rm p}/\sqrt{L^{i}_{k}}$, respectively.

(2) Second clustering operation.

According to the above clustering method, we merge measurements in $\dot{\tilde{{\bi Z}}}_{k}$ if their Mahalanobis distance $d_{\rm M}\lpar \dot{{\bi z}}^{i}_{k}\comma \; \dot{{\bi z}}^{j}_{k}\rpar $ is below the threshold $\sqrt{\delta_{{\rm P}_{\rm G}}}$,

(18)$$d_{\rm M}({\dot{{\bi z}}^{i}_{k}},{\dot{{\bi z}}^{j}_{k}}) = \sqrt {{({\dot{{\bi z}}^{i}_{k}} - {\dot{{\bi z}}^{j}_{k}})}^{\rm T}{\left[ {\matrix{ {{(\sigma _{\rm p})}^2(1/L_k^i + 1/L_k^j )} & 0 \cr 0 & {{(\alpha \sigma _{\rm p})}^2(1/L_k^i + 1/L_k^j )} \cr } } \right]}^{ - 1}({\dot{{\bi z}}^{i}_{k}} - {\dot{{\bi z}}^{i}_{k}})} $$

We also define a new measurement set as $\ddot{\tilde{{\bi Z}}}_{k}=\lcub \ddot{\tilde{{\bi z}}}^{i}_{k}\comma \; 1\le i\le \ddot{N}_{\rm m}\rcub $ after the second clustering operation, where $\ddot{\tilde{{\bi z}}}^{i}_{k}=\lcub \ddot{\tilde{{\bi z}}}^{i}_{k}\comma \; E^{i}_{k}\rcub $, $E^{i}_{k}=\lcub F^{i}_{n\lpar k\rpar }\vert 1\le n\le N^{i}\rcub $. N̈ _m is the number of clusters, $\ddot{{\bi z}}^{i}_{k}$ is the mean of position measurements in the i-th cluster, and $E^{i}_{k}$ is the frequency parameter set which contains multiple frequency measurements $F^{i}_{k}$ and $N^{i}_{k}$ is the number of $F^{i}_{k}$. In this paper, $\ddot{\tilde{{\bi Z}}}_{k}$ is used as the measurement set feeding the tracking. To simplify processing for tracking, we take $\sigma_{\rm p}/\sqrt{2}$ (the maximum after the first and the second clustering) in latitude and $\alpha\sigma_{\rm p}/\sqrt{2}$ in longitude as the unified standard deviations of positioning for all new measurements in the form of Equation (13).

4.2. Track score combined with features

The frequency measurement and the frequency state can be both regarded as Gaussian distributions. To better describe the feature similarity between emitters, Jensen-Shannon Divergence (JSD) (Lin, Reference Lin1991) is used to measure the distance of two frequency distributions. Assuming two Gaussian distributions p ₁(x) = N(x; m₁, P₁) and p ₂(x) = N(x; m₂, P₂), the JSD distance between them is defined as

(19)$$JSD(p_1\ \vert p_2)=\displaystyle{{1}\over{2}}(KLD(p_1\ \vert p_3)+KLD(p_2\ \vert p_3)),p_3=\displaystyle{{p_1+p_2}\over{2}}$$

(20)$$KLD(p_1\ \vert p_3)=\displaystyle{{1}\over{2}}\{\ln(\hbox{det}({\bi P}_3)/\hbox{det}({\bi P}_1))+ \hbox{tr}({\bi P}_1{\bi P}_3^{-1})-\hbox{dim}({\bi x})+({\bi m}_1-{\bi m}_3)^{\rm T} {\bi P}_3^{-1}({\bi m}_1-{\bi m}_3)\}$$

(21)$$KLD(p_2\ \vert p_3)=\displaystyle{{1}\over{2}}\{\ln(\hbox{det}({\bi P}_3)/\hbox{det}({\bi P}_2))+ \hbox{tr}({\bi P}_2{\bi P}_3^{-1})-\hbox{dim}({\bi x})+({\bi m}_2-{\bi m}_3)^{\rm T} {\bi P}_3^{-1}({\bi m}_2-{\bi m}_3)\}$$

(22)$${\bi m}_3=\displaystyle{{{\bi m}_1+{\bi m}_2}\over{2}}$$

(23)$${\bi P}_3=\displaystyle{{{\bi P}_1+{\bi m}_1^T{\bi m}_1+{\bi P}_2+{\bi m}_2^T{\bi m}_2}\over{2}}-{\bi m}_3^T{\bi m}_3$$

where KLD is Kullback-Leibler Divergence (KLD) (Duchi, Reference Duchi2007) and dim(·) represents the dimensionality of features. We define $\hat{E}^{j}_{k-1}$ as the frequency state of the j-th track at frame k − 1, and $E^{i}_{k}$ as the i-th frequency measurement at frame k,

(24)$$\hat{E}^j_{k-1}=\{\hat{F}^j_{m(k-1)} \vert 1\le m\le M^j\},\hat{F}^j_{m(k-1)}=[\hat{\mu}^j_{m(k-1)},\hat{\sigma}^j_{m(k-1)},I^j_{m(k-1)}]$$

(25)$$E^j_k=\{F^i_{n(k)} \vert 1\le n\le N^i\},F^i_{n(k)}=[\mu^i_{n(k)},\sigma^i_{n(k)},L^i_{n(k)}]$$

where M ^j and N ⁱ may not be equal. We define the distance between two sets as

(26)$$d(E^i_k,\hat{E}^j_{k-1})=\min_{1\le n\le N^i}\min_{1\le m\le M^j} JSD(N(x;\mu^i_{n(k)},(\sigma^i_{n(k)})^2)\ \vert\vert N(x;\hat{\mu}^j_{n(k-1)},(\hat{\sigma}^j_{n(k-1)})^2))$$

To convert the distance into likelihood, we use Fuzzy C-Means (FCM) clustering (Jain et al., Reference Jain, Murty and Flynn1999) to calculate the membership degree of feature information. This clustering algorithm, which assumes each sample can belong to multiple categories, is a soft partition method based on an optimal function. Providing that there are N _s samples from N _c categories given a weighted power b, the FCM clustering problem can be expressed as

(27)$$\eqalign{ \min &\quad \sum^{N_{\rm s}}_{i=1}\sum^{N_{\rm c}}_{j=1}(u_{ij})^b(d_{ij})^2 \cr \hbox{s.t.}&\quad u_{ij}\in [0,1], i=1,2,\ldots,N_{\rm s}\quad j=1,2,\ldots,N_{\rm c} \cr &\quad \sum^{N_{\rm c}}_{j=1}u_{ij}=1,\forall i\quad 0< \sum^{N_{\rm s}}_{i=1}u_{ij}< N_{\rm s},\forall j }$$

where u _ij and d _ij represent the membership degree and the distance between the i-th sample and the j-th category respectively. Using the Lagrange multiplier method, the fuzzy membership degree is given by

(28)$$u_{ij}=\displaystyle{{1}\over{\displaystyle\sum^{N_{\rm c}}_{k-1}\left(\displaystyle{{d_{ij}}\over{d_{ik}}}\right)^{\displaystyle{{2}\over{b-1}}}}}=\displaystyle{{(d^2_{ij})^{\displaystyle{{1}\over{1-b}}}}\over{\sum^{N_{\rm c}}_{k-1}\lpar d_{ij}^2\rpar ^{\displaystyle{{1}\over{1-b}}}}}$$

Through the above formula, we can obtain u _ij between $E^{i}_{k}$ and $\hat{E}^{j}_{k-1}$ within the gating area. Taking into account the measurement error, the error of the estimated position and the maximum distance, the gating distance between the current measurement position and estimated position at a previous time is set to satisfy the constraint,

(29)$$(\ddot{{\bi z}}^i_k-{\bi H}_{k-1}{\bi x}^j_{k-1})^T (\ddot{{\bi z}}_k^i-{\bi H}_{k-1}{\bi x}^j_{k-1})\le {\bi v}_{\rm max}{\bi T}_{k-1}/a+\sqrt{2}\alpha\sigma_{\rm p}$$

where v_max is the maximum speed. The increment of the i-th track's score with feature information can be defined as

(30)$$\Delta S^j(k) = \left\{ {\matrix{ {\omega \Delta S_{{\rm kin}}^j (k) + (1 - \omega )\Delta S_{{\rm fea}}^j (k),} \hfill & {i > 0} \hfill \cr {\ln (1 - P_{\rm D}),} \hfill & {i = 0} \hfill \cr } } \right.$$

where $\Delta S^{j}_{\rm kin}\lpar k\rpar $ and $\Delta S^{j}_{\rm fea}\lpar k\rpar $ are the increments of the kinematic and feature scores at frame k respectively, ω is the weight of kinematic information and $\Delta S^{j}_{\rm kin}\lpar k\rpar $ is calculated by Equation (2).

(31)$$\Delta S^j_{\rm fea}(k)=\ln (u_{ij}/p_e)$$

where p _e denotes the possibility of a measurement belonging to an extraneous (false or new) target.

4.3. Prediction and update

The prediction and update of the track information include kinematic and feature information. As the state transition equation is nonlinear, we use a first-order Extended Kalman Filter (EKF) to predict and update the motion state. For an EKF, it is important to calculate the state transition matrix F_k−1, which is a Jacobi matrix obtained by using a partial derivative of the state vector in Equation (9), and it is obtained as

(32)$${\bi F}_{k-1} = \left[ {\matrix{ 1 & {\displaystyle{{T_{k - 1}{\rm sec}(\varphi _{k - 1}^{\rm m} )} \over a}} & {\displaystyle{{{\dot \lambda }_{k - 1}T_{k - 1}\tan (\varphi _{k - 1}^{\rm m} )} \over {a\cos (\varphi _{k - 1}^{\rm m} )}}} & {\displaystyle{{{\dot \lambda }_{k - 1}{(T_{k - 1})}^2\tan (\varphi _{k - 1}^{\rm m} )} \over {2a^2\cos (\varphi _{k - 1}^{\rm m} )}}} \cr 0 & 1 & 0 & 0 \cr 0 & 0 & 1 & {T_{k - 1}/a} \cr 0 & 0 & 0 & 1 \cr } } \right]$$

The detailed process of prediction and update can refer to a standard EKF (Jazwinski, Reference Jazwinski2007), which is not described here. The frequency information of the track is updated but not predicted, and the update can be expressed as

(33)$$\hat{E}^j_k=J(\hat{E}^j_{k-1}, F^i_{n(k)})\cup F^i_{l(k)}$$

where $F^{i}_{n\lpar k\rpar }$ and $F^{i}_{l\lpar k\rpar }$ represent the frequency measurement used for estimation of an existing frequency and a new frequency, respectively, that satisfy

(34)$$ \vert \hat{\mu}^j_{m(k-1)}-\mu^i_{n(k)} \vert < \Delta f_1\quad \hbox{and}\quad \vert \hat{\sigma}^j_{m(k-1)}-\sigma^i_{n(k)} \vert < \Delta\sigma_1$$

(35)$$ \vert \hat{\mu}^j_{m(k-1)}-\mu^i_{l(k)} \vert >\Delta f_2,\forall m$$

That means the frequency measurement within tolerances Δ f ₁ and Δσ ₁ is used to update, and it is added as a new frequency if there is a difference of frequency excess Δ f ₂ compared with all frequency states. J(·) denotes that $F^{i}_{n\lpar k\rpar }$ is used to update $\hat{E}^{j}_{k-1}$ and the update of frequency information of the m-th emitter of the j-th ship is

(36)$$\hat{\mu}^j_{m(k)}=\displaystyle{{I^j_{m(k-1)}\hat{\mu}^j_{m(k-1)}+L^i_{n(k)}\mu^i_{n(k)}}\over{I^j_{m(k-1)}+L^i_{n(k)}}}$$

(37)$$(\hat{\sigma}^j_{m(k)})^2=\displaystyle{{I^j_{m(k-1)}[(\hat{\mu}^j_{m(k-1)})^2+(\hat{\sigma}^j_{m(k-1)})^2]+ L^i_{n(k)}[(\mu^i_{n(k)})^2+(\sigma^i_{n(k)})^2]}\over{I^j_{m(k-1)}+L^i_{n(k)}}}-(\hat{\mu}^j_{m(k)})^2$$

(38)$$I^j_{m(k)}=I^j_{m(k-1)}+L^j_{n(k)}$$

We use a logic-based method to perform track extraction. Tracks are confirmed by M-of-N test (we use 3-of-4 test) and are terminated after K (we set K = 3) missed detections. The new frequency of target is confirmed and outputted after the times of update is not less than a certain number (we use two times) in the whole process of tracking.