Tracking and Track Management of Extended Targets in Range-Doppler Using Range-Compressed Airborne Radar Data

Ship tracking facilitates a comprehensive insight into maritime traffic situations and ensures its safety and security. However, with the current operational surveillance systems, detecting several maritime threats is still a major challenge. In this article, we propose a supportive ship tracking concept using an airborne-based radar sensor. The proposed tracking algorithm is suitable for dense multitarget scenarios. Tracking is performed in the range-Doppler domain. The primary advantage of using the range-Doppler domain is that ships even with low radar cross section moving with certain line-of-sight velocity appear out of the clutter region, thus improving their detectability. In addition, a powerful track management system is also developed to handle false targets. The simulated and real experimental results from the DLR’s airborne radar sensors, F-SAR and DBFSAR, are presented to prove the concept.


I. INTRODUCTION
M ARITIME transport is considered as the backbone of international trade and the world economy. While the seaborne trade is in continuous expansion, the increasing number of vessels for fulfilling the demands of the global population has stimulated the development of a robust, reliable, and accurate maritime surveillance system [1].
State-of-the-art sensors that are used for maritime traffic situational awareness are the automatic identification system (AIS) [2] and marine radars [3]- [5]. These kinds of systems are appealing and serve as the dominating source for the localization and tracking of vessels [6]- [8]. Some of the major limitations of these systems are given as follows.
1) Not all ships, especially the smaller ones, are equipped with AIS transponders. 2) Marine radars suffer from limited visibility.
3) The targets of low radar cross section (RCS) and also those targets that are lying in the shadow and blind sector of the marine radar are difficult to detect. To deal with these limitations, spaceborne and airborne radars are the desirable choices. One distinguishing feature of airborne radars is their flexibility to collect data with very high resolution and with shorter revisit and longer observation times. The performance of these radar systems is mainly restricted by their endurance, which can be overcome in the future by installing them on high-altitude platforms (HAPs) or high-altitude pseudo satellites (HAPSs) that are flying in the stratosphere for several days, weeks, months, or even years [9], [10]. Unlike marine radars, with airborne radars, the targets are observed and tracked differently, as illustrated in Fig. 1.
As shown in Fig. 1(a), in the case of the airborne radar, the target is tracked as long as it is illuminated by the antenna beam that itself moves at the speed of the aircraft. The antenna in the shown example is nonsteerable, has wide beamwidth, and is typically designed for synthetic aperture radar (SAR) systems. On the other hand, in Fig. 1(b), the marine radar uses a long rotating antenna to transmit a very narrow horizontal beam [11]. The narrow antenna beam of marine radars allows them to measure directly the range and azimuth (bearing) of the target for a much longer time, limited mainly by the revolutions per minute (RPM) of the antenna.
Several tracking approaches are available on ship tracking using the marine radar [5], [12], [13]. However, current stateof-the-art tracking algorithms published in the open literature are not well explored for tracking ships with moving airborne radar platforms [14]. The available algorithms for such platforms were originally designed for tracking road vehicles [15]. They use low range-resolution data where the vehicles of interest, in most cases, occupy not more than a single resolution cell. For such cases, tracking can be executed easily at the cost of limited detection performance. For tracking extended targets, such as ships, additional efforts are needed (see Section III).
A potential overall processing concept of maritime surveillance using airborne radars is shown in Fig. 2.
As shown in Fig. 2, the major blocks of the processing flowchart are detection, tracking, inverse SAR (ISAR) imaging, and AIS data fusion. ISAR imaging and AIS data fusion are out of the scope of this article.
We use range-compressed (RC) radar data (single-channel or multichannel data) for airborne-based maritime This  surveillance [16]. With the airborne or future HAP or HAPS radars, the signal-to-noise ratio (SNR) in many cases is sufficiently large. Therefore, there is no need to use fully focused radar images. This saves the processing time and paves the way for future real-time capability.
Ships are detected in the range-Doppler domain after applying the azimuth fast Fourier transform (FFT) to the RC data. In the Doppler domain, ships even with low RCS, moving with sufficient line-of-sight (LOS) velocity, are shifted out of the clutter region, thus improving their detection capability [see Fig. 4(a)].
Detected ships are then tracked to reconstruct their trajectories. Like detection, tracking is also performed in the  range-Doppler domain. Tracking in range-Doppler has the benefit that overlapping target signals in the time domain can be separated in the Doppler domain if they are moving with different LOS velocities. This is because the targets moving with different radial velocities are shifted to different Doppler frequencies. This is very advantageous for tracking multiple targets, especially in dense multitarget scenarios. The target tracks in the range-Doppler domain are needed for the following.
1) Generating high-resolution ISAR image sequences after successively extracting all the required ship data (see Fig. 27), which later on aid to target recognition and identification [17], [18]. 2) Mapping the detections accurately on the ground after computing additionally the direction-of-arrival (DOA) angle of the target [19]. The ground tracks can later be compared with the AIS data for validation purposes (see Fig. 29). Therefore, tracking in range-Doppler is indispensable and a prerequisite for the aforementioned applications. This article proposes a novel range-Doppler-based tracking algorithm for tracking multiple extended targets using RC airborne radar data. Extended targets in the data may give rise to more than one detection per target per frame [5]. Therefore, the extracted centroids of the targets are tracked using a suitable state-space motion model. The overall multitarget tracking (MTT) system is designed using an SQLite [20] database structure. A powerful track management system is also developed, which runs simultaneously within the tracker to automatically recognize and terminate false or ghost targets. An additional challenge associated with tracking in range-Doppler is the Doppler aliasing (or back-folding) effect. In this article, the problem is briefly discussed with real experimental results. To validate the proposed range-Doppler-based tracking algorithm, the simulated and real linearly and circularly acquired experimental data from DLR's (German Aerospace Center) airborne F-SAR [21] and DBFSAR [22] system are presented.
The remainder of this article is organized in the following order. Details about the experimental airborne radar data are provided in Section II. Ship detection and clustering algorithms are discussed briefly in Section III. The concept of the extended target tracking in the range-Doppler domain is explained in Section IV. Track management is explained in Section V followed by the Doppler aliasing phenomena in Section VI. In Section VII, the results from the simulated and real experimental data are presented. Section VIII concludes this article with a short summary and discussion.
II. EXPERIMENTAL AIRBORNE RADAR DATA Since, in the following, results based on real data are used for the explanation of our proposed algorithm, we will discuss  I   RADAR SYSTEM AND ACQUISITION GEOMETRY PARAMETERS OF THE  SINGLE-CHANNEL X-BAND LINEAR AND CIRCULAR FLIGHT TRACKS in the next paragraphs when and how these experimental data were acquired. A two-day F-SAR flight campaign was conducted in June 2016 in the North Sea. All radar data, in a total of more than 1 TB, were acquired fully polarimetric and simultaneously in the X-and L-bands. On the first day of the campaign, the island Helgoland and the town Cuxhaven, including the coastal areas and ships of opportunity, were observed, mainly during linear flight tracks, but also during a circular track with the radar antenna pointing not to the circle center but to the opposite direction [see the red circle in Fig. 3 (top)]. On the second day, a dedicated experiment with a controlled ship operated by the German federal police was carried out. The ship moved with velocities of 0-20 knots (kn) between three different waypoints. The circular flight tracks, this time with the antenna pointing to the circle center, were flown with a radius of 5600 m, resulting in a total ship observation time of approximately 400 s (=6.7 min) per circle [see the red circle in Fig. 3 The acquisition geometry and the system parameters of the linearly and circularly acquired real single-channel airborne radar data are listed in Table I. Further details about the conducted flight campaigns can be found in [17].

III. DETECTION AND CLUSTERING
The methodology proposed in [16] is used to detect ships in the range-Doppler domain. In that paper, the authors have proposed a constant false alarm rate (CFAR)-based ship detection algorithm, which is adapted for the RC airborne data. After obtaining multiple pixel-based detections from a single ship, clustering is applied to group these detections as a single "physical object." A standard density-based spatial clustering of applications with noise (DBSCAN) algorithm [23] is used to form the ship clusters. An example of clustering a ship signal in the range-Doppler domain of real X-band HH polarized F-SAR data is shown in Fig. 4.
To perform tracking the center of the clustered ship, i.e., its Doppler frequency and range position [see cluster centroid in Fig. 4(b)] have to be estimated and tracked at successive times. Therefore, we have computed here the center of gravity (COG) of the cluster to be used for extended target tracking. Other metrics, such as the center corresponding to the maximum peak or the center of the bounding box, could also be used. For our investigated airborne radar data and for almost all the ships, we found COG to be more stable. It is worth mentioning that principally 2-D correlation-based methods can also be used for cluster center position estimation [24].

IV. TRACKING IN RANGE-DOPPLER
Target tracking is defined as the reconstruction of the motion of the targets with associated parameter estimates based on a sequence of noisy measurements at each time step. In our case, tracking is completely performed in the range-Doppler domain. The noisy measurements are the estimated cluster center positions (COG in Section III). Fig. 5 shows a schematic representation of the ship tracking concept.
As shown in Fig. 5, the RC data are first partitioned along the azimuth direction into small coherent processing intervals (CPIs). The length of the CPI is data and system-parameterdependent. For the F-SAR data, each CPI has 128 azimuth samples. Thus, with the pulse repetition frequency (PRF) of 2.4 kHz, each CPI has a duration of approximately 53 ms. After transforming individual CPIs to the Doppler domain, ship cluster center positions are estimated and tracked at each CPI.
The simplified processing steps of the proposed range-Doppler-based tracker are shown in Fig. 6. The tracking methodology can be summarized as follows.
1) Estimate the cluster center positions for each detected target (see Section III). 2) Store the target position and its related motion parameters in a database (see Section IV-A). 3) Associate the targets to already existing tracks for MTT (see Section IV-F). 4) Run track management after every t mng seconds for updating the tracks (see Section V). 5) Monitor the Doppler frequency of the target (see Section VI) for considering the Doppler back-folding. The details related to each processing block in Fig. 6 are sequentially presented in the following.

A. Database Structure
We have developed an SQLite database structure for storing the detection and clustering results and also for doing target tracking. The database has a table where each row of the table represents a detected target at each CPI. The column contains among others the target motion parameters. The database is designed for storing and tracking an arbitrary number of targets, limited only by the available memory and the SQLite limitations. An example of a database table is shown in Fig. 7.
The table shows that three exemplary targets were detected at CPI = 0. These targets are stored in the database and have their unique IDs (=unique row numbers), which increments automatically with subsequent detections in the following CPIs. Along the column of each row, parameters such as the positions and the extents of the targets and also the data patches belonging to the targets for ISAR imaging (dashed block in Fig. 6) are stored. The table also has additional columns, such as a relation and a predicted flag (PF). Their importance and need are explained in Sections IV-F and V.
The detected ships (i.e., the cluster centers) at each CPI are tracked by using a target motion model. Tracking can still be performed without using any motion model. This is valid only when the targets are detected at every time step, and the time steps are relatively small (53 ms in our case). However, without a motion model, gaps due to missing detections cannot be bridged. The ability of the motion model to fill the gaps in the detections is explained in Section IV-E.

B. State and Measurement Model
In this section, we define the state and measurement models used for tracking targets in the range-Doppler domain.
Assume that a target is moving with constant velocity (CV) on the ground. Such a target in RC data has an azimuth timedependent range and Doppler history [see Fig. 8(b) and (c)]. This is because the range between the moving radar platform and a specific point of the target on ground changes.
The target's CV motion on the ground with respect to time (in the Cartesian coordinates) is given as where x 0 and y 0 are the ground positions of the target at t = 0, and v x0 and v y0 are the along-track and across-track velocity components, respectively. If v p is the platform velocity, z t and z p are the target and platform altitude with respect to the ground, respectively, the target range r (t) as a function of time [see Fig. 8(b)] can be written as The corresponding Doppler history f a (t) as a function of time [see Fig. 8 where λ is the radar wavelength. Inserting (3) into (4), the first-order Taylor series approximation of the Doppler history [see Fig. 8(c)] can be written as [25] f where f DC and k a denote the Doppler shift and the Doppler slope, respectively. The terms f DC and k a are mathematically expressed as [25] f DC ≈ − 2 λr 10 where r 10 is the range between the antenna and the target at t = 0. Using (6) and (7) in (3), the range history of a target moving with CV on the ground [see Fig. 8(a)] can be approximated using Doppler parameters in the following way: For simplicity, we write (8) as r (t) ≈ r 10 + u r t + 1 2 a r t 2 (9) where u r and a r can be considered as the initial velocity and acceleration components of the range history and have values of −(λ/2) f DC and −(λ/2)k a , respectively. Inserting (5) into (8), the moving target range history as a function of the Doppler frequency can be written as [26] r ( f a ) ≈ r 10 Substitutingf a = f a − f DC , (10) can further be written as where v dop and a dop are the Doppler velocity and Doppler acceleration, respectively. Equation (11) is depicted in Fig. 8(d), where the quadratic behavior of the target's range history can be seen as a function of the Doppler frequency.
In this article, the target kinematics in the range-Doppler domain, which are assumed to evolve in time, are expressed by the following five-state vector Z (t k ) ∈ n , which is defined as: where t k is the azimuth time corresponding to the CPI k and the symbol is the definition sign. The estimates f a (t k ) andḟ a (t k ) are the Doppler frequency and its first derivative [see (5)], respectively, and r (t k ),ṙ (t k ), andr (t k ) are the range and its first-and second-order derivatives [see (8)], respectively. The componentsḟ a (t k ),ṙ (t k ), andr (t k ) of the target state vector in (12) are written aṡ To estimate the target kinematics shown in (12), the target's Doppler frequency and range can be approximated by CV and constant acceleration (CA) motion models [27], respectively. For the range history, a CV motion model can also be used assuming that the range is piecewise linear due to the very small-time intervals.
The received measurements (i.e., the cluster centers) are the Doppler frequencies ( f am ) and the range positions (r m ) (see Section III). The 2-D target measurement vector m(t k ) ∈ l is defined as The state-space motion model and the measurements discussed in this section are incorporated within the framework of the Kalman filter (KF) [28]. It is computationally fast since it uses the current measurement and the estimated state and its uncertainty from the previous CPI in order to estimate the true state at the current CPI. Nonlinear filters, such as the extended/unscented KFs [29], could also be applied [see (11)]; however, their investigations are out of the scope of this article.

C. Kalman Filter
The KF works in three stages. In the first stage, also known as the initialization stage, the initial state meanẐ (0|0) and the initial posterior error covariance of the state are set up. We assume here thatẐ (0|0) is the first detected position of the target because the true position of the target at k = 0 cannot be known in advance. The initialization of the error covariance matrix of the target state is explained in Section IV-D2.
In the prediction stage, the target dynamics [see (12)] are predicted by using a motion model described in Section IV-B. The KF equations for the prediction stage arê whereẐ (k + 1|k) is the predicted state at k + 1 given the state estimateẐ (k|k) at k. The matrices P(k + 1|k) and Q are the n × n covariance matrices of the predicted state and the process noise, respectively. The state transition model F applied to the previous statê Z (k|k) to predictẐ (k + 1|k) is written as where T CPI is the time interval corresponding to one CPI (see Fig. 5). In the last stage, i.e., the update stage, the predicted state and the covariance matrix from (17) and (18) are updated based on the received noisy measurements [see (16)]. The updated target stateẐ (k + 1 | k + 1) and the associated covariance matrix P(k + 1|k + 1) then becomê where I is the identity matrix and H is the l × n observation matrix that transforms the state space to the measurement space. Here, n and l are the dimensions of the target and the measurement states, respectively. The matrix K (k + 1) is the Kalman gain, and D(k + 1) is the innovation (or the residual vector), which are expressed as (22), commonly known as the measurement prediction  II   DETAILS OF THE REAL AIRBORNE RADAR DATA TAKES USED FOR  THE INITIALIZATION OF KF MATRICES. THE DATA WERE ACQUIRED  USING LINEAR FLIGHT TRACKS. NOTE THAT THE SHIP SHOWN IN  THE TABLE IS ALSO USED IN THE CIRCULAR  FLIGHT EXPERIMENT   TABLE III ESTIMATED R MATRIX FOR DIFFERENT DATA TAKES. FIRST AND SECOND VALUES OF THE DIAGONAL MATRIX CORRESPOND TO VARIANCES IN DOPPLER AND RANGE, RESPECTIVELY covariance matrix or innovation covariance matrix, is written as where R is the l × l measurement noise covariance matrix. The significance of the matrix S(k + 1) and the innovation vector D(k + 1) in terms of MTT are explained in Section IV-F.

D. Initialization of the Kalman Filter Matrices
Kalman filtering requires the initialization of the P(0|0), R, and Q matrices. They are set offline based on either simulated or already available real experimental data and are kept constant throughout the filtering process. This is because the online estimation of these matrices is difficult as we do not have the ability to observe the process that we are estimating [30]. However, a fairly good approximation of R can be made by using the cluster centers themselves. In this section, we give, in our case, the plausible initialization values of R, P(0|0), and Q matrices. To initialize these matrices, three linearly acquired real X-band HH polarized experimental RC data were investigated. Each data had a real single controlled moving ship. The specifications of the data (e.g., wavelength and PRF) are listed in Table I, and details about the experiments were already given in Section II. Additional details of the radar data and the ships in the data are given in Table II. We found that these three acquisitions were sufficient to retrieve suitable initialization parameters of the KF.

1) Initialization of R:
The R matrix is a diagonal matrix that stores the variance of the deviations between the received measurements and the true measurements (see Table III). The received measurements are the cluster centers that were computed at each CPI using the methodology presented in Section III. The true measurements in our case are assumed as the higher order polynomial fitting to the received noisy measurements. This is because the single true position of an extended target in high-resolution data cannot be known accurately because of the ship dimensions.
After having at each time instant the true and the received measurements from the ship, the variances are computed over azimuth. Fig. 9 shows the square root of the computed variances, i.e., the standard deviation plot of the measured Doppler and range positions of the ship in data take I (see Table II).
As shown in Fig. 9, the uncertainties for both Doppler and range vary significantly over time. This is mainly because of the unstable estimated cluster centers. The factors that contribute to the instability of the cluster centers are the ship extents (see Table II), their amplitude fluctuations in the data, and the Doppler and range bin sizes, which, for the investigated F-SAR data, are approximately 20 Hz and 0.3 m, respectively.
In order to set the initials for the R matrix, the mean standard deviations in Doppler and range are computed using the results shown in Fig. 9. The obtained values are then squared to get the mean variances that are listed in Table III for three different data takes.
Using Table III as a reference, the initial values of the R matrix are set as diag(350, 5). The initialization is, however, data-dependent and is usually recommended to be set higher as the obtained estimates.
2) Initialization of P(0|0): The initialization of P(0|0) is based on our knowledge about the initialization error of the target states. If it is assumed that the initial state is far from its true state, then P(0|0) should be set higher. We also assumed that we have no information about the typical values of P(0|0), and therefore, we can write P(0|0) = σ 2 p0 I n , where I n is the identity matrix and n is the number of target states (see (12),  Table II). The noisy cluster centers can be seen in the figure. The ship track is shown in the range-Doppler domain using different initializations of σ 2 p0 .
where n = 5). The only tuning parameter of P(0 | 0) is σ 2 p0 . Fig. 10 shows the impact of different initializations of σ 2 p0 on the trajectory of a real moving ship in the Doppler domain.
In Fig. 10, it is observed with lower initialization values (i.e., P(0|0) = 0.01I 5 ); the predictions are trusted more and more at the beginning. As a result, the cluster centers are ignored, and larger errors are observed. Therefore, we set P(0 | 0) = 1000I 5 in order to avoid the discrepancies at the beginning of the estimation. Although it is fixed initially, it gets updated with successive times minimizing the estimation error [see (21)].
3) Initialization of Q: The Q matrix represents the expected uncertainties in the state equations. The uncertainties are due to the modeling errors, measurement errors, and the approximations made in the derivations [see (5) and (8)]. The initialization of Q is not straightforward and is constructed intuitively. To do so, we have tabulated different initial values of Q after fixing R and P(0 | 0)(see Table III and Fig. 10) and checked how their initialization affects the Doppler frequency and range position accuracies. The true state to compute the RMSE (root mean square error) in the Doppler and range is again the higher order polynomial fit to the measurements (as done for Fig. 9). The results are shown in Table IV where we performed the investigations again for the same three different data takes.
From Fig. 10 and Table IV, it is clear that the best accuracy is achieved with a high value of P(0|0) and low value of Q.
In a concluding remark for the matrix initialization of the KF, it can be said that R should be estimated offline using either simulated or already available real experimental data. The matrix P(0|0) should be set higher, and Q should be set lower. From now on, for further investigations, we set R = diag(350, 5), P(0|0) = 1000I n , and Q = 0.01I n . 4) Filter Consistency: With the given initializations, the consistency of the KF for the real data with extended targets and the airborne acquisition geometry is checked by computing the normalized innovation squared (NIS) [31]. NIS is a statistical hypothesis test for detecting the model mismatch. If the initializations of the KF filter are properly set then for l = 2 degrees of freedom, NIS should be chi-square (χ 2 ) distributed, and within the two-sided 95% confidence interval, the NIS bound (or the acceptance region) is given by [0.05, 7.38]. If a lot of the computed NIS points lie outside the given bound, then there is a serious mismatch and the initializations have to be set properly. In this article, we show the NIS plot of data take I (see Table II) in Fig. 11. As shown in Fig. 11, the computed NIS values of a real ship in data take I are within the given bounds (horizontal red dashed lines in the figure). Here, only the upper bound is of interest because the lower bound is practically zero (0.05), and it is shown that the NIS is well below the upper bound [31]. Although not shown, the NIS for the other two data takes also satisfies the criteria. This means that the proposed model and the selected initializations can be used.

E. Maximum Manageable Gap as Motion Model Performance Measure
In real scenarios, it is expected that there will be missing measurements at certain azimuth times. This could be due to the lower backscatter received from the moving target or when the target is not illuminated by the antenna beam [see Fig. 25(a)] [17]. In such situations, the KF is able to give a predicted position [see (17)].  However, inaccurate predictions of the missed measurements may further lead to a track loss. Therefore, we have investigated and compared "CV only" and (CV +CA) motion models of the KF to assess their performance in the presence of missed measurements. In the "CV only"-based motion model, along with the Doppler, the range history is also modeled as a CV by assuming it as piecewise linear. The (CV+CA) motion model was already discussed in Section IV-B.
For the investigations, the gaps in the data were artificially introduced (see Fig. 12). The gap duration was successively increased to a point after which the target track is lost. The performance achieved using the motion models is compared with the one where no motion model is considered. The results are shown in Table V again for the same three data takes (see Table II). The position accuracies are also evaluated.
From Table V, it can be seen that, although, in terms of the position accuracy, all methods perform similarly, the (CV+CA)-based target motion model can still be considered as a better model based on the track duration (e.g., 9.5 s for data take III).
In the following, we explain the use of relations and unique IDs (see Fig. 7) to form the tracks of an arbitrary number of targets.

F. Relation Generation Using Database
For tracking targets, it is essential to establish a relation between the detected targets at a given CPI and the already existing tracks from the previous CPI. Fig. 13 shows an example of the concept of relation generation using a database structure for tracking three consecutively detected targets.
In Fig. 13, there are three detected targets at CPI = 0. These targets are assigned to unique IDs. The unique IDs are the row numbers of the database table, and for a newly added row, it is incremented automatically.
If the same three targets are detected in the next CPI, the unique IDs are again generated (i.e., new rows are added into the database). The relation columns of these newly added rows are now updated with the unique IDs of the same targets from the previous CPI (see the relation columns at CPI = 1 and the Unique ID columns at CPI = 0 in Fig. 13). After the tracking is over, a link between the unique rows belonging to the same target is established (arrows of the same color in Fig. 13). For instance, in the figure, it is shown that unique IDs 7-4-1 (green arrow) belong to the same target and likewise for the other two targets. A relation of −1 indicates that the target is detected for the first time and has no relation with any previously existing tracks.
1) Target to Track Association: Previously, it was shown that unique row numbers belonging to the same target at different CPIs are linked to create the target tracks. Generating such links is possible only when the target detections at the current CPI are associated with the tracks from the previous CPI. In tracking literature, this is termed data association. The concept of data association for a single target tracking is shown in Fig. 14.
In Fig. 14 at CPI = k − 1 (previous CPI), only one measurement is observed, which is denoted as m (0) (k − 1). The track is initialized for this target, and the motion parameters are stored in the database. At CPI = k (current CPI), two measurements are observed m (1) (k) and m (2) Fig. 14 . Data association is now performed as a two-step procedure.
First, the position of the previously detected target is predicted at the current CPI. The positionm (0) (k | k − 1) is the KF-based predicted measurement of m (0) (k − 1) [see (23)]. Second, a validation region centered around the predicted measurement (which is now an established track) is created. This is done in order to eliminate unlikely observation-to-track pairings. Detections falling within the gate of the track are considered for updating the track [see Fig. 14, where m (1) (k) lies in the gating region ofm (0) (k|k − 1)]. Detections lying outside the gate m (2) (k) in this case could either be a false target or a new target. The tracker will initiate a new track for such cases.
Rectangular gating shown in Fig. 14 with the fixed size is one of the possibilities. An ellipsoidal gating could also be used to select measurements for a track [32]. Its limitations in the context of the state-of-the-art tracking filters are explained in Section IV-F2. Furthermore, a discussion on measurement uncertainty i.e., when more than one measurement lies inside the gating region of a given track, is also presented in Section IV-F2.
The rectangular gating criteria shown in Fig. 14 are mathematically expressed as where |m f a (k) −m f a (k|k − 1)| is the offset between the detected and predicted Doppler frequency and |m r (k) − m r (k|k − 1)| is the offset between the detected and predicted range [see (16) and (17)]. The terms f thres and r thres are the width (along Doppler) and length (along range) of the rectangular region, respectively. From Fig. 14 and (25) and (26), it is clear that the extents of the rectangular search window have to be first determined to associate detections to their corresponding target tracks. Before determining the extents of the rectangular gate, it is first necessary to investigate how the offsets shown on the left-hand side of (25) and (26) vary so that reasonable values of gate extents can be set. These offsets are obtained directly from one of the KF equations shown in (23). In the equation, the disparity between the received measurement and the predicted measurement is clearly reflected. An example of these offsets (also known as innovations) estimated for a real moving ship in data take I (see Table II) is shown over the azimuth time in Fig. 15.
From Fig. 15, the maximum observed offsets in the Doppler and range are approximately 40 Hz (almost twice the Doppler bin size) and 4 m (almost 14 times the range bin size), respectively. The factors that cause such variations were already explained in Section IV-D1. For the relation generation, the size of the rectangular window must be larger than these estimated offsets. A wise choice for the extents of the search window can be set three times the maximum offsets. This will prevent the track loss even if the maximum offsets are a bit larger than what was observed in Fig. 15.
2) Measurement Uncertainty: There are situations when multiple measurements fall within the gate of a single-track causing measurement uncertainty. Two popular algorithms that are used for resolving such uncertainties and are applied in real-world problems are the global nearest neighbor (GNN) [33] and the joint probabilistic data association filter (JPDAF) [5].
Given a set of measurements and tracks, both GNN and JPDAF first compute the Mahalanobis distances between the track and the measurements and compare them to a gating threshold to select the candidate measurements for each track. GNN then uses the Hungarian algorithm [34] for assigning the most likely measurements to the existing tracks [33]. JPDAF, on the other hand, uses all the measurements within the validation region of the track, evaluates the measurementto-track association probabilities, and combines them to update the target state and covariance [32], [35].
For tracking an unknown number of targets, it is advisable to use a separate track management logic along with GNN and JPDAF [36].
The Mahalanobis distance used in GNN and JPDAF for selecting the potential measurements is sensitive to the innovation covariance of the track [see (24)]. If the covariance is set smaller, then the distance between the track and the measurement is more than the validation region of the track. As a result, the true measurement might fall outside its validation region [see Fig. 16(a)]. Especially, for our data, such a situation is very likely because the resolution in Doppler is very coarse due to coherent integration of a small number of azimuth samples, and the resolution in the range is very fine (<40 cm). As a result, the center of the ship is jumping around significantly over time. The larger the ship is, the higher could the jump be. Increasing the R matrix can increase the innovation covariance, but the states will deviate from the radar measurements because the reduced Kalman gain causes less dependence on the measurements [see (20)].
Furthermore, in Fig. 16(b), when a true target is not detected in several subsequent CPIs, its innovation covariance grows bigger. Due to high covariance, the false measurement is closer to the true target track, and therefore, the threat of false target assignment is very high. This situation can by far be avoided by setting a threshold in the covariance, but, again, this is something to be set manually [5]. In addition, with an  increasing number of targets and measurement uncertainties, the computation time of JPDAF also increases.
3) Proposed Method: We propose here a rather simple but fast and efficient MTT approach shown in Fig. 17. To explain the method, an example is shown in Fig. 18.
As shown in Fig. 17, to solve the measurement uncertainty, an assignment matrix is initially computed by providing the number of tracks and detections. The elements of the assignment matrix are the Euclidean distance offsets between the existing tracks and the available detections. The assignment matrix is computed separately for Doppler and range because both are measured in different units. Based on our knowledge of Doppler and range bin sizes, we create a rectangular gating region around each track for selecting potential target-originated measurements [see Fig. 14 and (25) and (26)]. After applying the gating, if there is more than one candidate measurement of a confirmed track (see Fig. 18 where there are three detections within the validation region of track T2), we then compute the Mahalanobis distances [37] only between the track and its validated measurements [see Fig. 19(b)]. The Mahalanobis distance is used because the dimensions are not equally weighted. To estimate the Mahalanobis distance, the covariance matrix is required, which, in this case, is the covariance of the innovation vector D [see (24)].
Assume that, in a CPI, there are N vm numbers of validated measurements lying within the rectangular validation region of a particular track. We then calculate the Mahalanobis distances {d Mi } N vm i=1 between the track and its validated measurements. The measurement corresponding to the minimum Mahalanobis distance is considered for the track update. For a given track, the minimum Mahalanobis distance d min can be mathematically expressed as (27) where d M is the Mahalanobis distance that is written as where D is the innovation vector and S is its covariance matrix [see (23) and (24)]. An example to resolve the measurement uncertainty is shown in Fig. 19, which is specifically solved for the situation shown in Fig. 18. As shown in Fig. 19(b), after using (27), the detections D1, D3, and D4 are assigned to tracks T1, T2, and T3, respectively. The unassociated detections D2 and D5 initiate their own new tracks. Unlike JPDAF, our approach is not sensitive to the innovation covariance for selecting the candidate measurements for the track because we use a Euclidean distance-based approach [see (25) and (26)]. We found this as an advantage not just to apprehend the true measurement but also to discard false target assignment to a confirmed track.
This simple method of data association is found sufficient for most of the investigated airborne radar data.

V. TRACK MANAGEMENT
Every unassigned detection initiates a tentative track. If the detection belongs to a real target, then the target is expected to be detected in several subsequent CPIs. Such a track then becomes a confirmed track. Once the target moves out of the antenna beam, it is no longer detected and, therefore, should be terminated. Also, if there exist some false targets (clutter), they should also be terminated as they are not originated from real targets. To do this, we employ a track management scheme for updating the confirmed tracks and also to terminate the finished tracks and false targets [38]. The concept of track management is illustrated in Fig. 20.
Track management runs periodically within the tracker. For the investigated F-SAR data, it is run after, e.g., every 2 s (can be defined by the user), which corresponds to approximately 40 CPIs as per the current F-SAR system configuration.
As an example, in Fig. 20, there are three targets. One of them is a real target, and the others are false targets. Each target initiates its own track in its very first detection (given by the green dots). A value of zero is then assigned in the PF column of each target, which implies that these are detected targets (see Fig. 7). In the next time step, only the real target is detected, and the rest are not detected. In the absence of detection, the target position is predicted by the KF (given by the blue dots). For the predicted position, we now assign a value of one in its PF column. As the tracking continues, at CPI = 40 (i.e., nearly after 2 s), the individual target tracks are extracted, and their track lengths are determined. As shown in the figure, each target should at least be tracked for 2 s. This is because the observed time of less than 2 s is too small to terminate the targets (see false targets in the figure). For the tracks that are equal to 2 s, their prediction percentages (PP) are estimated using where, for each extracted track, N pred is the total number of CPIs in which the target is only predicted and N L is the total extracted track length, i.e., 40 CPIs in the shown example. We set the PP threshold to 70%, which means that the targets that are "only predicted" for more than 70% of the time without having their corresponding detections are terminated. In the figure, the real target is also terminated after 10 s (it was last observed at 8 s). An example of tracking a single real ship target in data take III (see Table II) in the presence of false targets and with and without using track management is shown in Fig. 21.
As shown in Fig. 21,Target ID_0 belongs to a real ship target, whereas Target ID_1 and Target ID_2 are the false targets, which are also termed "ghost targets." With track management, it was possible to terminate such target tracks after a short time. VI. DOPPLER ALIASING Another major challenge while tracking targets in the range-Doppler domain is the effect of Doppler aliasing. It is a special condition that occurs mainly when the Doppler shift of the target is larger than the PRF limit of the radar system [39]- [41]. Larger Doppler shifts of the target are the consequences of either its high radial velocity [see (4)] or squinted radar acquisition geometry or too low PRF. As a result, the target appears back-folded (or aliased) in the range-Doppler domain.
To illustrate the aliasing effect, an exemplary range history of a single moving target in the range-Doppler plane is shown in Fig. 22. The target track has already been initialized (see CPI k 0 in Fig. 22), and the detections are tracked in subsequent CPIs. As the tracking continues, at CPI = k als , the KF predicts the Doppler frequency that is smaller than −PRF/2 [see (17)]. The target detection to be used to assign to this track is found approximately one PRF apart and nearly at the same range (aliased detection). This is because the target's Doppler frequency shifts that are smaller than −PRF/2 are back-folded and detected from PRF/2 onward (see the ambiguous trajectory in Fig. 22). As a result, the tracker initiates a new track for the aliased detection (see the orange curve in Fig. 22). However, in reality, there is only one target in the scene (complete blue curve in Fig. 22). Detecting and correcting the aliased Doppler frequency are essential in order to extract the true range-Doppler history of the target.
For resolving the Doppler aliasing, the Doppler frequency of the aliased detection must be updated during tracking/track management. If the target is moving toward −PRF/2 and the aliasing occurred in the right, as shown in Fig. 22, then the Doppler frequency of the aliased detection is subtracted by one PRF, and if the target is moving toward PRF/2 and the aliasing is detected on the left, then one PRF is added to the aliased Doppler frequency. By doing this, new track initiation from the aliased target can be avoided, and a single unambiguous target track in range-Doppler can be extracted. Fig. 23 shows the real single ship trajectory reconstruction before and after solving the Doppler aliasing. The PRF used is approximately 1.2 kHz. In Fig. 23(a), after tracking the ship for nearly 1 s, it is found aliased. That is why, in Fig. 23(a), the Target ID_0 dies after a certain time (see Section V), and a new target is initialized (Target ID_1). However, after the aliasing is corrected, only one target appears during the entire illumination time [see Target ID_0 in Fig. 23(b)].

VII. SIMULATION AND EXPERIMENTAL RESULTS WITH REAL DATA A. Simulation Results
In the simulation, we considered three moving targets. The targets have overlapping range histories in the time domain [see Fig. 24(b)]. The effects of missing target detections, false targets, and Doppler aliasing are artificially included in the simulation. Simulated radar data acquisition and moving target parameters are given in Table VI.
The tracking results of the three simulated moving targets are shown in Fig. 24. Tracking is implemented in the range-Doppler domain; only for visualization purposes, we show the target tracks in the time domain.
In Fig. 24(b), the gaps are artificially introduced for Target ID_0 (blue) and Target ID_2 (green) that are filled by the KF-based predictions [see (17)]. Moreover, they are also found aliased at 2.4 and 7.9 s, respectively, which, in the end, resulted in several ambiguous targets, as shown in Fig. 24(a). The reason is that the PRF is not sufficiently high (see Table VI where the PRF is 1.5 kHz) to capture the larger Doppler shifts caused due to their higher ground velocities (20 and 35 m/s for Target ID_0 and Target ID_2, respectively). In addition, without using track management, the artificially added false target tracks are also never terminated. However, from Fig. 24(b), it is clear  that, after using KF (see Section IV-C) and track management (running here approximately after every 2 s) and recognizing and correcting the Doppler back-folding (see Section VI), the number of targets is reduced from 14 to 3 (plus 2 artificially introduced false targets).

B. Real Experimental Results
Using Single-Channel Data 1) Tracking Results: Fig. 25(b) shows the tracking results of a real moving ship in circularly acquired X-band HH polarized F-SAR airborne radar data. Additional information related to the tracked targets is listed in Table VII. Details about the acquired circular radar data were already given in Table I. In the data, the detected and clustered real ship signal (see Section III) is tracked using the (CV+CA) motion model-based KF (see Section IV-C) and the data association method presented in Section IV-F. From detection to tracking, all the algorithms are implemented in the range-Doppler domain.
As shown in Fig. 25(b), there are total of nine tracked targets. It is also observed that, if the gaps are significantly larger (in the order of several minutes), multiple new tracks from a single target can be created (see Fig. 25(b) and Table VII where the single police ship has five different target tracks).
In addition, some false targets or clutter are also detected, which are very well handled by the track management system (see Section V). Since there were larger gaps in the data, the tracks were, therefore, checked after every 4 s (see Table VII where the false targets are tracked for more than 4 s).
Using the tracked range information of the ship, the data belonging to the ship signal are also successively extracted in  FIG. 25(B) the time domain and are shown in Fig. 25(c). The extracted ship signals are later on used for generating high-resolution ISAR image sequences. We point out here that the gaps in the radar data shown in Fig. 25(a) are due to the small 3-dB azimuth antenna beamwidth of the F-SAR X-band acquisition (≈8 • ). Since the antenna also cannot be steered electronically or mechanically, during a circular flight with a ship moving in the circle center, the crosswind may cause a significant yaw angle so that the ship is not always illuminated.
To evaluate the tracker performances, there exist several metrics, such as optimal subpattern assignment (OSPA) [42] or generalized OSPA (GOSPA) [43]. Metrics such as OSPA are already implemented on real X-band marine radar data [5], [8]. In order to evaluate such metrics, information such as the true ground position of the target and the time duration of its true trajectory in the region of interest must be known.
In this article, due to the use of only a single receiving channel (see Table I), we do not have the angle information of the target. As a result, the tracks cannot be projected on the ground, and therefore, comparison with the AIS-based ground-truth data is not possible. In addition, as shown in Fig. 1(a) and (c), due to the antenna pattern weighting, wider antenna beam, and the nonideal platform motion, it is also impossible to know exactly when the target will enter and exit the beam. Therefore, the true time on target cannot be determined.
In the simulations presented in Section VII-A, we have demonstrated that, in terms of the total number of tracked targets, disregarding the crucial components of a tracker can significantly worsen its performance [see Fig. 24(a)].
2) Target's Range Position Accuracy in Circular Data: In this section, we present the range position accuracy results of the controlled police ship using both X-and L-band real single-channel airborne radar data acquired simultaneously during the circular flight track. Details about the flight experiment and the police ship are given in Section II. The X-band circular radar data with the police ship were already shown in Fig. 25(a), and the details about the L-band circular radar data can be found in [17]. Fig. 26. Range position differences computed over the azimuth time for the controlled police ship in both X-and L-band real airborne radar data acquired simultaneously during the circular flight track. The corresponding X-band circular radar data were already shown in Fig. 25(a). Gaps marked in the figure, e.g., between 250 and 305 s, occur when the target is not illuminated by the radar antenna beam. In the figure, red and green vertical dashed lines with colored guard zones correspond to the time instants when the ship moving direction was perpendicular (either toward or away) and parallel or antiparallel with respect to the aircraft flight direction, respectively.
For the accuracy assessment, the measured range position of the ship is compared with its true range position. The measured range position of the ship is derived from the COG of the ship cluster (see Section III). Since the COG of the ship cluster is a single position and can be anywhere on the ship, the measured range position of the ship is biased by the ship size (see Table II where the ship's length and beam are 66 and 11 m, respectively).
The reference range position of the ship is calculated by using the GPS position of the ship and the position of the aircraft at the same azimuth time. Note here that the position of the GPS antenna on the ship is not known. Therefore, the reference range position of the ship is also biased by the ship size.
After having the measured and the reference range positions of the ship at each azimuth time, the range position differences are computed, and the results are shown in Fig. 26.
From Fig. 26, it can be seen that there are significant gaps in the computed range position differences in the X-band circular data compared to the L-band data. This is because the X-band 3-dB azimuth antenna beamwidth is in the order of 8 • (see Table I), which is smaller than the L-band antenna beamwidth (≈18 • ). Moreover, during a circular flight, the crosswind significantly changes the yaw angle of the aircraft, and the target is not always illuminated by the narrower antenna beam of the X-band radar.
Furthermore, it is also observed that the range position differences vary significantly during the total illumination time. Ship amplitude fluctuations and its constantly changing moving direction with respect to the aircraft (magenta curve in Fig. 26) are responsible for changing the measured and true range positions of the ship because of the bias caused by the ship size.
The minimum and maximum range position differences for the X-and L-band circular radar data are (0.002 m, 23.25 m) and (0.005 m and 32.22 m), respectively. For both the frequency bands, the observed maximum range position differences are within 66 m, which is the length of the ship. Fig. 27. Some ISAR images obtained from the extracted RC ship data shown in Fig. 25(c). The data were extracted after using the tracking information. The images correspond to different azimuth times, and the coherent integration time for each image is 1.7 s. Finally, the calculated RMSEs of the tracked ship range history in X-and L-band radars are 13.18 and 13.75 m, respectively. For the ship that has a length of 66 m and a beam of 11 m, the calculated RMSE can be considered very good.
3) Preliminary ISAR Imaging Results: One primary advantage of circular flights is that, due to the very long observation time, dozens or even hundreds of ISAR images of the ships can be obtained under different aspect angles and used later on for ship classification and recognition purposes. In Fig. 27, a small cut of such an ISAR image sequence is shown.
The images were generated from the extracted data patch depicted in Fig. 25(c) by applying the ISAR processor discussed in [17]. A coherent integration time of 1.7 s was used for each ISAR image. The achievable azimuth (crossrange) resolution is a function of the radar data wavelength, the coherent integration time, and the target motion itself [18].

C. Geocoding Results Using Multichannel Data
The airborne radar data discussed so far in this article were acquired using only a single receiving channel (see Section II). With only one receiving channel, it was not possible to do along-track interferometry that would allow the DOA angle estimation of the detected target and the projection of radar-based target detections on the ground.
To complete the processing chain of our airborne radar-based maritime surveillance concept, we have provided some preliminary geocoded results using real multichannel experimental radar data.
Multichannel data are exploited to compute additionally the DOA angle of the detected target. With the known DOA angle, the geographical aircraft position, and the Euler angles from the aircraft's inertial measurement unit (IMU), the target tracks in the range-Doppler domain can directly be mapped onto the ground so that the target's geographical coordinates are obtained.
In November 2019, a multichannel flight campaign using DLR's DBFSAR system was carried out in the North Sea near Cuxhaven, Germany. More details related to the radar sensor and the multichannel data can be found in [22]. The aircraft was flying at an altitude of approximately 2400 m above ground and was also equipped with a dual-channel AIS receiver. The multichannel DBFSAR system and acquisition parameters are listed in Table VIII. Parameters such as the radar wavelength and 3-dB azimuth antenna beamwidth were already given in Table I. Some preliminary results with the multichannel radar data from the North Sea flight campaign are shown in Fig. 28.
As shown in the figure, there are more than just two ships in the data as revealed by the binary detection map [see Fig. 28(d)]. The detected and clustered ships were tracked in the range-Doppler domain using the approach proposed in this article. In total, there were 21 tracked targets. Since there is track management running in parallel, the false tracks are terminated after a short time [see ID_4 and ID_6 in Fig. 28(e)]. Fig. 28(f) shows the ISAR image generated for ship ID_15. The ship name is HAM 316, which is a dredger of dimensions 129 m × 22 m. This ship contains several cranes and booms on the deck; one of them causes a strong multipath reflection clearly recognizable in the ISAR image.
The tracked range histories shown in Fig. 28(e) were then mapped on the ground after computing the DOA angle of the targets. The projected ground tracks were then compared with the available AIS ground truth data. For demonstration purposes, we have selected here only two ships namely HAM 316 and LANGELAND (ID_15 and ID_17 in Fig. 28(e), respectively). The geocoding results are shown in Fig. 29.
We found a good agreement between the geocoded target detections and the AIS tracks, giving an average position estimation accuracy of better than 30 m. Such an accuracy for an extended target can be considered very good for radar-based ship monitoring using an airborne radar sensor. A detailed investigation on ISAR imaging, multichannel data preprocessing, DOA estimation, and geocoding techniques is out of the scope of this article.
VIII. CONCLUSION This article proposed a novel range-Doppler-based ship tracking algorithm suitable for the RC airborne radar data. The algorithm is applicable for linear and circular flight tracks, as shown in this article. However, there exists also no restriction for arbitrary flight tracks if the aircraft position and the associated Euler angles are efficiently incorporated in the estimation, and the data are processed using small CPIs. Due to the algorithm structure and the use of mainly azimuth FFTs, the algorithm is expected to have the real-time capability. For sure, this requires an efficient implementation in a parallelized way on a multicore or multiprocessor computer, taking also into account graphical processing units. The algorithm includes a detector and clustering module followed by a tracking module. The former detects and calculates the COG of the ship cluster, which is tracked over azimuth. The latter has the following components: 1) a KF with motion models based on CV and CA to handle missed target detections; 2) a simple but efficient data association block to do MTT and also to resolve detection uncertainties; 3) a powerful track management system running simultaneously within the tracker to terminate the false and already finished tracks; 4) a Doppler aliasing block to extract the true target range history and, hence, the unambiguous Doppler history by identifying and considering the Doppler back-folding. The complete tracking module, including track management, is designed using an SQLite database as the core. This article also showed some preliminary ISAR imaging and geocoding results by using the tracking information of the ship. The proposed methodology was validated by using simulated and real experimental radar data acquired with DLR's airborne radar sensors, F-SAR and DBFSAR, from altitudes of up to 5600 m above ground. Although not covered and discussed in detail in this article, we want to point out again that the proposed algorithm can easily be adapted for a multichannel radar system as well. With multiple receiving antennas arranged in the along-track direction, the tracked target Doppler and range histories can be projected on the ground (i.e., the actual geographical position can be obtained) by using the estimated DOA angles, as shown in the example in Fig. 29. To the best of our knowledge, such a detailed and comprehensive investigation of ship tracking using moving airborne radar sensors is, up to now, not available in the open literature. From 2016 to 2017, he was involved in the Earth Observations Applications Mission (EOAM) Project at IIRS, where he worked on automatic water feature extraction by the synergistic use of optical and microwave data. Since 2017, he has been with the Microwaves and Radar Institute, German Aerospace Center (DLR-HR), Weßling, Germany, as a Ph.D. Student. He has authored or coauthored over 5 peer-reviewed journal articles and 15 conference papers. His research interests include the development of multitarget tracking algorithms, and maritime moving target indication techniques using single-channel and multichannel synthetic aperture radar (SAR) and inverse ISAR imaging.
Mr. Joshi has won one of five prizes during the Three Minute Thesis (3MT) Contest organized by the IEEE during the IEEE Radar Conference 2020 in Florence, Italy. In 2018, he received the Franz-Xaver-Erlacher-Förderpreis (Franz-Xaver Erlacher Promotion Award)-an award for promoting young scientists and Ph.D. students granted by the Society of Friends of DLR. He was also invited as a Guest Scientist in the Congreso Internacional ultidisciplinario de Ingenierías 2020 (CONIMI 2020) organized by the Universidad Politécnica Juventino Rosas, Mexico.