Micro-Motion Extraction for Marine Targets by Multi-Pulse Delay Conjugate Multiplication and Layered Tracking

Abstract

The detection and recognition of marine targets can be improved by utilizing the micro-motion induced by ocean waves. However, distinguishing and tracking overlapping targets is challenging when multiple extended targets overlap in the range dimension of the radar echo. In this paper, we propose a multi-pulse delay conjugate multiplication and layered tracking (MDCM-LT) algorithm for micro-motion trajectory tracking. The MDCM method is first applied to obtain the conjugate phase from the radar echo, which enables high-precision micro-motion extraction and overlapping state identification of extended targets. Then, the LT algorithm is proposed to track the sparse scattering points belonging to different extended targets. In our simulation, the root mean square errors of the distance and velocity trajectories were better than 0.277 $\mathrm{m}$ and 0.016 $\mathrm{m}/\mathrm{s}$, respectively. Our results demonstrate that the proposed method has the potential to improve the precision and reliability of marine target detection through radar.

1. Introduction

Micro-motion plays an important role in radar target detection and classification [1]. The vibration of a marine target caused by wave slapping is a unique feature that can be exploited for target identification and voyage monitoring [2]. Extracting vibration parameters can be challenging when performing long-distance radar detection, since the projection of micro-motion on the radar line-of-sight may be weak. The phase-derived measurement technique has been employed to achieve high-precision range estimation and motion feature extraction [3,4,5,6], which utilizes the principle that a half-wavelength translation of a target induces a $2\pi$ phase delay in radar echo.

In the range domain of wideband radar echo, the marine vessel can be regarded as an extended target composed of strong scattering points. Extended target detection is studied in [7,8,9]. In order to capture the target motion variation over time, the target tracking method was introduced. The Bayesian filtering framework is commonly used, including the Kalman filter (KF) [10] and the particle filter [11]. Moreover, the tracking problem can be divided into sub-problems of association between adjacent points and be solved through the dynamic programming method [12,13]. In long-distance scenes, the tangential coverage range of the radar beam may reach several hundreds of meters, and there may exist multiple targets in the radar sight. The main challenge in multiple extended target tracking is correctly associating different targets. For the multi-target tracking issue, the pairing between tracks and measurements is usually achieved through data association algorithms, such as nearest neighbor [14] and joint probabilistic data association [15]. In [16], the multiple points were separated into several groups, and the centers of the groups were tracked through the probability hypothesis density filter. In [17], the multi-Bernoulli filter algorithm was adopted to track the centers of groups.

However, when tracking the overlapping extended targets, current methods rely on motion model predictions or prior knowledge of target movements. It is difficult to directly acquire the overlapping state between targets. Additionally, existing extended target tracking methods primarily focus on tracking isolated points or the center of entire targets, without fully utilizing the motion correlation between scattering points, which reflects the structure of the extended target and offers new insights for target characteristic extraction. There is a pressing need to explore how to use the micro-motion of scattering points to enhance extended target tracking.

To address the problems, we proposed an extended target tracking algorithm based on the multi-pulse delay conjugate multiplication (MDCM) and layered tracking (LT). First, we analyzed the micro-motion features of marine targets and established the radar echo model. Then, the MDCM method was applied to calculate the conjugate phase, which was further used to differentiate the coincident points (superimposed by multiple scattering points) and independent points. With the assistance of independent points, the conjugate phases of the coincident points were corrected, and LT was finally achieved.

The reminder of this paper is organized as follows. Section 2 describes the movement model and the radar echo simulation of the marine target. Section 3 presents the MDCM-LT method. In Section 4, the proposed algorithm is verified through the simulation and the experiment. Furthermore, Section 5 concludes this paper.

2. Marine Target Micro-Doppler Modeling

2.1. Micro-Motion Model

According to sea-keeping theory, the stability and seaworthiness of a marine target are related to its size, shape, and weight distribution. Furthermore, various targets may exhibit distinct motion characteristics under identical sea conditions.

Linear strip theory [18] assumes the amplitude of the marine target motion is proportional to the amplitude of the ocean wave. The ocean wave forms from a superposition of regular waves as

$$\xi \left(t\right)=\sum _{i=1}^{\infty }{\xi }_{i}cos({\omega }_{i}t+{\phi }_i)$$

(1)

where ${\xi }_i$, ${\omega }_i$, and ${\varphi }_i$ denote the amplitude, the angular frequency, and the phase of the i-th wave, respectively.

The micro-motion of a marine target encompasses three translational and three rotational degrees of freedom, and it can be expressed as

$$a_u\left(t\right)=\sum _{i=1}^{\infty }{\mathrm{Rao}}_u\left({\omega }_i\right){\xi }_{i}cos({\omega }_{i}t+{\phi }_i)$$

(2)

where u ranges from one to six, thereby representing the target movements in six degrees of freedom. Furthermore, the response amplitude operator, ${\mathrm{Rao}}_u\left({\omega }_i\right)$, denotes the ratio between the movement amplitude to the regular wave amplitude at frequency ${\omega }_i$. In sea states of medium to low intensity, the movements along the longitudinal direction (heaving and pitching) are particularly influential on the micro-motion characteristics of the marine target. To simplify the vessel motion modeling, we focused on heaving and pitching movements.

We adopted the microtome section theory [19] to create the micro-motion model of the marine vessel. Two coordinate systems are established: the global coordinate system O-$XYZ$ and the vessel-fixed coordinate system $O_v$-$xyz$ (shown in Figure 1). The radar is positioned at point O, while the longitudinal axis of the vessel is oriented along the $+X$ axis. We assume that the coordinate of the vessel center, $O_v$, is $(R,0,0)$. $h\left(t\right)$ and $\theta \left(t\right)$ denote the height variation caused by heaving and the pitch angle of pitching, which are expressed as

$$\begin{array}{c}h\left(t\right)={\displaystyle \sum _{i=1}^{\infty }}{\mathrm{Rao}}_h\left({\omega }_i\right){\xi }_{i}cos({\omega }_{i}t+{\phi }_i)\hfill \\ \theta \left(t\right)={\displaystyle \sum _{i=1}^{\infty }}{\mathrm{Rao}}_{\theta }\left({\omega }_i\right){\xi }_{i}cos({\omega }_{i}t+{\phi }_i)\hfill \end{array}$$

(3)

where ${\mathrm{Rao}}_h\left({\omega }_i\right)$ and ${\mathrm{Rao}}_{\theta }\left({\omega }_i\right)$ are the response amplitude operators of heaving and pitching, respectively.

The vessel is separated into multiple sections along its longitudinal direction, and the coordinate of each section is calculated as

$$\begin{array}{c}{x}^{\prime }=R+\mathsf{\Delta }xcos\theta \left(t\right)\hfill \\ y^{\prime }=0\hfill \\ z^{\prime }=h\left(t\right)+\mathsf{\Delta }xsin\theta \left(t\right)\hfill \end{array}$$

(4)

where $\mathsf{\Delta }x$ is the distance between the section and $O_v$.

Then, the distance between the section and the radar is calculated by using the following equation:

$$\begin{array}{c}{R}^{\prime }\left(t\right)={\left\{{[R+\mathsf{\Delta }xcos\theta \left(t\right)]}^2+{\left[h\right(t)+\mathsf{\Delta }xsin\theta (t\left)\right]}^2\right\}}^{\frac{1}{2}}\hfill \end{array}$$

(5)

By calculating the derivative of Equation (5) with respect to time, the velocity of the section is

$$v\left(t\right)=\frac{1}{b}\left\{\left[h\left(t\right)-\mathsf{\Delta }xsin\theta \left(t\right)\right]\frac{dh\left(t\right)}{dt}+\mathsf{\Delta }x\left[h\left(t\right)cos\theta \left(t\right)-Rsin\theta \left(t\right)\right]\frac{d\theta \left(t\right)}{dt}\right\}$$

(6)

where $b=\sqrt{{(R+\mathsf{\Delta }xcos\theta )}^2+{(h\left(t\right)+\mathsf{\Delta }xsin\theta )}^2}$.

Due to $\mathsf{\Delta }x$ and $h\left(t\right)$ being much smaller than R, b could be approximated as R. Based on Equation (6), the velocity difference between two sections with interval $\mathsf{\Delta }{x}_1$ and $\mathsf{\Delta }{x}_2$ is calculated as follows

$$\begin{array}{c}\begin{array}{cc}\hfill \mathsf{\Delta }v& =v_{\mathsf{\Delta }{x}_1}\left(t\right)-v_{\mathsf{\Delta }{x}_2}\left(t\right)\hfill \\ \hfill & \approx \frac{1}{R}\left[h\left(t\right)cos\theta \left(t\right)\frac{d\theta \left(t\right)}{dt}-Rsin\theta \left(t\right)\frac{d\theta \left(t\right)}{dt}-sin\theta \left(t\right)\frac{dh\left(t\right)}{dt}\right]\left(\mathsf{\Delta }{x}_1-\mathsf{\Delta }{x}_2\right)\hfill \end{array}\end{array}$$

(7)

In Equation (7), $\theta \left(t\right)$ and $h\left(t\right)$ are variables associated with the overall vessel movement, and do not change with $\mathsf{\Delta }x$. Therefore, the difference in velocities between sections depends linearly on their relative position.

2.2. Radar Echo Model

The radar waveform is the chirp pulse signal with carrier frequency $f_c$. We regard the vessel sections as scattering points and construct the radar echo by superimposing the echoes from strong scattering points as follows

$$s_r\left(t\right)=\sum _{i=1}^I{\sigma }_i\mathrm{rect}\left[\frac{t-2R_i^{\prime }\left(t\right)/c}{T_w}\right]exp\left\{j2\pi \left[f_c\left(t-\frac{2R_i^{\prime }\left(t\right)}{c}\right)+\frac{1}{2}K{\left(t-\frac{2R_i^{\prime }\left(t\right)}{c}\right)}^2\right]\right\}$$

(8)

where I denotes the total number of the scattering points, $\mathrm{rect}(·)$ is the rectangular function, $T_w$ is the pulse width, K is the modulation ratio, c is the microwave speed, and ${\sigma }_i$ and $R_i^{\prime }\left(t\right)$ represent the scattering amplitude and the distance of the i-th scattering point, respectively.

The matched filtering is applied to the radar echo to obtain the range profile of the vessel. The reference signal in the frequency domain is denoted as $H\left(k\right)$, which is calculated by applying Fourier transform on the transmitting signal. The output of the matched filter is calculated as follows

$$\begin{array}{c}\begin{array}{cc}\hfill x\left(t\right)& =\mathrm{IFFT}\left[S_r\left(k\right)·H^{*}\left(k\right)\right]\hfill \\ \hfill & =\sum _{i=1}^I{A}_i\mathrm{sinc}\left[B\left(t-\frac{2R_i^{\prime }\left(t\right)}{c}\right)\right]exp\left[-j\frac{4\pi R_i^{\prime }\left(t\right)}{\lambda }\right]\hfill \end{array}\end{array}$$

(9)

where $\mathrm{IFFT}(·)$ denotes the inverse fast Fourier Transform process, $S_r\left(k\right)$ is the Fourier transforms of $s_r\left(t\right)$, $H^{*}\left(k\right)$ is the conjugate of $H\left(k\right)$, $A_i$ is the amplitude of the output result, B is the bandwidth of the chirp signal, and $\lambda$ is the wavelength.

The Doppler frequency of the i-th scattering point could be aquired by calculating the phase derivation of the range profile in Equation (9) as

$$f_{d,i}\left(t\right)=-\frac{2}{\lambda }\frac{d{R}_i^{\prime }\left(t\right)}{dt}$$

(10)

3. Micro-Motion Tracking for Extended Target

3.1. Micro-Motion Feature Extraction through MDCM

To simplify the analysis, we focused on the micro-motion of a single scattering point. The amplitude and envelope terms in Equation (9) are denoted as $A_r$, and the range profile of the m-th pulse can be expressed as

$$x_i(m,t_n)=A_{r}exp\left[-j\frac{4\pi R_i^{\prime }\left(t_m\right)}{\lambda }\right]+n\left(t_m\right)=A_{r}exp\left[-j\frac{4\pi R_i^{\prime }(m{T}_r+t_n)}{\lambda }\right]+n(m{T}_r+t_n)$$

(11)

where $t_n$ is the fast time in pulse, $t_m$ represents the sampling time in the m-th pulse, $T_r$ is the pulse repetition time, $R^{\prime }\left(t_m\right)$ is the distance between the point and the radar at time $t_m$, and $n\left(t_m\right)$ denotes the echo of the environment clutter.

The approximation of $R^{\prime }\left(t_m\right)$ in Equation (5) can be calculated using the first-order Taylor expansion as

$$\begin{array}{cc}\hfill R^{\prime }\left(t_m\right)& ={\left\{R^2+h{\left(t_m\right)}^2+\mathsf{\Delta }{x}^2+2Rcos\theta \left(t_m\right)\mathsf{\Delta }x+2h\left(t_m\right)sin\theta \left(t_m\right)\mathsf{\Delta }x\right\}}^{\frac{1}{2}}\hfill \\ \hfill & \approx R+\frac{1}{2R}\left[2Rcos\theta \left(t_m\right)\mathsf{\Delta }x+h{\left(t_m\right)}^2+\mathsf{\Delta }{x}^2+2h\left(t_m\right)sin\theta \left(t_m\right)\mathsf{\Delta }x\right]\hfill \end{array}$$

(12)

Since $\mathsf{\Delta }x$ and $h\left(t_m\right)$ are much smaller than R, the last three terms in the square brackets in Equation (12) can be neglected for analysis. Then, $R^{\prime }\left(t_m\right)$ can be approximated as

$$R^{\prime }\left(t_m\right)\approx R+cos\theta \left(t_m\right)\mathsf{\Delta }x$$

(13)

where $\theta \left(t_m\right)$ represents the pitch angle at time $t_m$.

$\theta \left(t_m\right)$ varies slowly over pulses. Therefore, we assume $\theta \left(t_m\right)$ to be linear during processing time with a velocity $v_{\theta }$. Substitute Equation (13) into Equation (11), and calculate the l-pulse MDCM as

$$\begin{array}{cc}\hfill y(m+l,T_n)& =x_r(m+l,t_n)·x_r^{*}(m,t_n)\hfill \\ \hfill & =A_{c}exp\left\{j\frac{8\pi \mathsf{\Delta }x}{\lambda }sin\frac{l{v}_{\theta }{T}_r}{2}sin\left[\theta \left(t_m\right)+\frac{v_{\theta }{T}_r}{2}\right]\right\}\hfill \end{array}$$

(14)

where $A_c$ is the amplitude, and l is the number of delayed intervals.

In Equation (14), the phase term is named as the conjugate phase. The range of the phase after the MDCM process is constrained by $sin\left(l{v}_{\theta }{T}_r/2\right)$, which is named as the reduction ratio $\eta$. By changing the delayed interval l, we can manipulate the value of $\eta$ and decrease the phase variation caused by micro-motion, thereby unwrapping the ambiguity phase. $l{v}_{\theta }{T}_r$ approaches zero when l does not exceed a few hundred, and $T_r$ is at the microsecond level. In such circumstances, $\eta$ is approximately equivalent to $l{v}_{\theta }{T}_r/2$, and is linearly proportional to l.

3.2. Layered Tracking of Extended Target

In this subsection, the tracking of the extended target was divided into two layers: the center tracking in the upper layer, and the scattering points tracking in the lower layer. In each time step, the measurement points $\left\{z_i(k+1)\right\}$ were initially partitioned into sets belonging to different extended targets. The wrong conjugate phases in $\left\{z_i(k+1)\right\}$ will be corrected on the basis of upper layer state prediction, ${\widehat{X}}_c(k+1|k)$. Then, the average of $z_i(k+1)$ will be adopted as the measurement values of upper layer, $Z_c(k+1)$. The tracking process of the lower layer and the upper layer were both implemented using the KF method. The framework of the proposed LT method is depicted in Figure 2. The details of the LT algorithm are explained in the following parts.

3.2.1. Measurement Partition

The appropriate partitioning of measurements is the prerequisite for effective tracking. We employed the distance partitioning principle [20] to divide the measurements at the initial frame. The Mahalanobis distance between two measurements, $z_i$ and $z_j$, is calculated as follows

$$d_{i,j}={\left[{\left(z_i-z_j\right)}^T\mathbf{R}\left(z_i-z_j\right)\right]}^{\frac{1}{2}}$$

(15)

where $\mathbf{R}$ is the distance weight matrix. When $d_{i,j}$ is smaller than the empirical threshold, the two measurements will be considered to belong to the same extended target.

When the scattering points of two extended targets coincide in the range dimension, the measurement partition becomes difficult, as more than one scattering points merge into one measurement, which we name the coincident point. To simplify the analysis, we assumed that a coincident point is superimposed by two scattering points with constant amplitudes. The phase of the coincident point, ${\varphi }_{\mathsf{\Sigma }}$, could be calculated as

$${\varphi }_{\mathsf{\Sigma }}=d{\varphi }_A+e{\varphi }_B$$

(16)

where ${\varphi }_A$ and ${\varphi }_B$ are the phases of the two points, respectively. d and e represent the proportion of ${\varphi }_A$ and ${\varphi }_B$ in ${\varphi }_{\mathsf{\Sigma }}$, respectively.

The MDCM under two delayed intervals, $l_1$ and $l_2$, is applied to the coincident point, and the conjugate phases are denoted as ${\varphi }_{\mathsf{\Sigma },1}$ and ${\varphi }_{\mathsf{\Sigma },2}$, respectively. Here, $l_1$ is small enough that the conjugate phases of two points are unambiguous, while $l_2$ is larger, and the phases are aliased. Let ${\varphi }_{A,2}$ and ${\varphi }_{B,2}$ denote the conjugate phases of the two scattering points under interval $l_2$; ${\varphi }_{\mathsf{\Sigma },2}$ can be calculated as

$${\varphi }_{\mathsf{\Sigma },2}=d{\varphi }_{A,2}+e{\varphi }_{B,2}$$

(17)

Based on the linear variation of the conjugate phase with l, ${\varphi }_{\mathsf{\Sigma },2}$ can be rewritten as

$$\begin{array}{cc}\hfill {\varphi }_{\mathsf{\Sigma },2}& =d\left(\frac{l_2}{l_1}{\varphi }_{A,1}-2k_A\pi \right)+e\left(\frac{l_2}{l_1}{\varphi }_{B,1}-2k_B\pi \right)\hfill \\ \hfill & =\frac{l_2}{l_1}\left(d{\varphi }_{A,1}+e{\varphi }_{B,1}\right)-\left(k_{\mathsf{\Sigma }}\right)\pi \hfill \\ \hfill & =\frac{l_2}{l_1}{\varphi }_{\mathsf{\Sigma },1}-\left(k_{\mathsf{\Sigma }}\right)\pi \hfill \end{array}$$

(18)

where ${\varphi }_{A,1}$ and ${\varphi }_{B,1}$ denote the conjugate phases of A and B, respectively, under interval $l_1$, $k_A$ and $k_B$ are the phase ambiguity integers, and $k_{\mathsf{\Sigma }}=2d{k}_A+2e{k}_B$. By choosing moderate values of $l_1$ and $l_2$, $k_{\mathsf{\Sigma }}$ is no longer a multiple of two. In this case, ${\varphi }_{\mathsf{\Sigma },2}$ can not be unwrapped correctly based on ${\varphi }_{\mathsf{\Sigma },1}$, while the single scatters still satisfies the linearity among different delayed intervals.

We therefore calculate the difference between the conjugate phases under $l_1$ and $l_2$ as

$$\mathsf{\Delta }\varphi =\left|{\varphi }_2-\left(\frac{l_2}{l_1}{\varphi }_1-2k\pi \right)\right|$$

(19)

where ${\varphi }_1$ and ${\varphi }_2$ denote the conjugate phases of the scattering point under $l_1$ and $l_2$, respectively. When $\mathsf{\Delta }\varphi$ is higher than the threshold ${\eta }_{\varphi }$, the point will be considered to be a coincident point and be assigned to multiple extended targets.

3.2.2. State Prediction and Update

The scattering points are detected from the range profile through the constant false alarm rate detection, which are input into the lower layer as measurements. The state vector of each measurement in frame k contains its distance $r\left(k\right)$ and conjugate phase $\omega \left(k\right)$ and is denoted as

$$X\left(k\right)={\left[r\left(k\right),\omega \left(k\right),\dot{r}\left(k\right),\dot{\omega }\left(k\right)\right]}^T$$

(20)

where $\dot{r}\left(k\right)$ and $\dot{\omega }\left(k\right)$ denote the first derivative of $r\left(k\right)$ and $\omega \left(k\right)$, respectively.

The state transition model is written as

$$\begin{array}{cc}\hfill \widehat{X}(k+1|k)& =FX\left(k\right|k)+V(k)\hfill \\ \hfill \widehat{P}(k+1|k)& =FP\left(k\right|k)F^T+Q\left(k\right)\hfill \end{array}$$

(21)

where $\widehat{X}(k+1|k)$ and $\widehat{P}(k+1|k)$ represent the predicting state and covariance matrix in frame $k+1$, respectively, $X\left(k\right|k)$ and $P\left(k\right|k)$ are the state matrix and covariance matrix after frame k, respectively, F is the transition matrix under a constant velocity model, and $V\left(k\right)$ and $Q\left(k\right)$ are the transition noise matrix and the process noise matrix, respectively.

Then, the tracker is updated as

$$\begin{array}{c}K(k+1)=\widehat{P}(k+1|k)H^T{\left[H\widehat{P}(k+1|k)H^T+R(k+1)\right]}^{-1}\\ X(k+1|k+1)=\widehat{X}(k+1|k)+K(k+1)\left[Z(k+1)-H\widehat{X}(k+1|k)\right]\\ P(k+1|k+1)=[I-K(k+1)H]\widehat{P}(k+1|k)\end{array}$$

(22)

where $K(k+1)$ is the Kalman gain in frame $k+1$, H is the measurement vector, $R(k+1)$ is the observation noise, $Z(k+1)$ is the measurement in frame $k+1$, I is the identity matrix, $X(k+1|k+1)$ and $P(k+1|k+1)$ are the state matrix and the covariance matrix after updating, respectively.

3.2.3. State Correction

Based on Equations (6) and (13), the averages of the radar distance and velocity of the scattering points in an extended target are calculated as

$$\begin{array}{c}\hfill {\overline{R}}^{\prime }\left(t\right)=R+cos\theta \left(t\right)\sum _{i=1}^n\frac{\mathsf{\Delta }{x}_i}{n}\\ \hfill \overline{v}\left(t\right)=\frac{h\left(t\right)}{a}\frac{dh\left(t\right)}{dt}+d\sum _{i=1}^n\frac{\mathsf{\Delta }{x}_i}{n}\end{array}$$

(23)

where n is the total number of scatters, $\mathsf{\Delta }{x}_i$ is the distance between each scatter and the target center, and

$$d=\frac{1}{a}\left[h\left(t\right)cos\theta \left(t\right)\frac{d\theta \left(t\right)}{dt}-sin\theta \frac{dh\left(t\right)}{dt}-Rsin\theta \left(t\right)\frac{d\theta \left(t\right)}{dt}\right]$$

(24)

When the distribution of scatters remains unchanged, ${\displaystyle \sum _{i=1}^n}\mathsf{\Delta }{x}_i/n$ is constant. As a result, ${\overline{R}}^{\prime }\left(t\right)$ and $\overline{v}\left(t\right)$ are only related to the overall motion of the extended target. Therefore, the state of the tracking center $(r_c,{\omega }_c)$ is regarded as the average state of scattering points and is calculated as

$$r_c=\sum _{i=1}^n\frac{r_i}{n},{\omega }_c=\sum _{i=1}^n\frac{{\omega }_i}{n}$$

(25)

where $r_i$ and ${\omega }_i$ denote the radar distance and conjugate phase of the i-th scattering point, respectively.

The center of the extended target is tracked using the same KF framework as in Section 3.2.2. During the tracking process, the prediction values of the target center, ${\widehat{r}}_c$ and ${\widehat{\omega }}_c$, were used to correct the conjugate phases of the coincident points. The correction process was based on the linear variation of micro-motion with $\mathsf{\Delta }x$, and each coincident point will be corrected to fit each extended target, respectively.

We assumed a coincident point to be distributed in the tracking gates of two extended targets at the same time and superimposed by two scattering points. For each extended target, if there exist more than two scattering points with correct conjugate phases, the states of the two nearest points, $\left(r_1,{\omega }_1\right)$ and $\left(r_2,{\omega }_2\right)$, will be carried out to calculate the revised conjugate phase as follows

$${\omega }_r={\omega }_1+\left(r-r_1\right)\frac{{\omega }_2-{\omega }_1}{r_2-r_1}$$

(26)

where ${\omega }_r$ is the revised phase.

If only one scattering point has the correct phase, we calculate the difference as ${\delta }_{\omega }={\widehat{\omega }}_0-{\omega }_0$, where ${\widehat{\omega }}_0$ and ${\omega }_0$ denote the prediction and the measurement in the last frame, respectively. Then, the revised phase is calculated as ${\omega }_r=\widehat{\omega }+{\delta }_{\omega }$, where $\widehat{\omega }$ and $\omega$ denote the prediction and the measurement in the current frame, respectively. If no single scattering points exist, the predicted phases of the corresponding lower layer trackers will be used as the phases of the coincident points directly.

In this way, the conjugate phases of the coincident points were corrected and then used for state update process.

4. Results and Discussion

4.1. Simulation Results of Micro-Motion Extraction

A vessel with no translation was considered in this simulation. The initial distance between the radar and the vessel center was 500 $\mathrm{m}$. The vessel was composed of three scattering points, and $\mathsf{\Delta }x$ of the three points were −10 $\mathrm{m}$, 0 $\mathrm{m}$ and 10 $\mathrm{m}$, respectively. The motion parameters were set as follows: the heaving amplitude was 0.5 $\mathrm{m}$, with a period of 10 $\mathrm{s}$; the pitching amplitude was 1.5°, with a period of 10 $\mathrm{s}$. The carrier frequency of the chirp pulse was 16 $\mathrm{GHz}$, the pulse repetition time was 100 $\mathsf{\mu }\mathrm{s}$, the number of pulses in one frame was 512, and the bandwidth was 200 $\mathrm{MHz}$. In this simulation, a sea clutter generation method considering the sea texture distribution, speckles, and sea spikes was adopted to generate the sea clutter echoes $n\left(t_m\right)$ [21]. The theoretical velocity and the Doppler-time diagram of the vessel under a signal-to-clutter ratio (SCR) of 20 $\mathrm{dB}$ are shown in Figure 3a and Figure 3b, respectively. In Figure 3b, we noticed that the velocity variation caused by heaving and pitching was too small to be observed under the velocity resolution of 0.18 $\mathrm{m}/\mathrm{s}$.

The initial phases of scattering points in each time step were extracted from the radar diagram and illustrated in Figure 4a. To unwrap the ambiguous phases, a 10-pulse MDCM was applied to the signal, and the conjugate phases were obtained, which were unwrapped and transformed into velocity values. The velocities derived from the conjugate phases of the simulated data with an SCR of 20 $\mathrm{dB}$ are plotted in Figure 4b. According to Figure 4b, the phase unwrapping was achieved effectively through the MDCM method, and the estimated velocity demonstrated good agreement with the theoretical velocity, with a root mean square error (RMSE) of 3.98 × 10${}^{-4}$$\mathrm{m}/\mathrm{s}$.

Figure 5 presents the RMSEs of the velocity estimation with an SCR ranging from 8 $\mathrm{dB}$ to 23 $\mathrm{dB}$. The phase measurement accuracy was directly influenced by the signal quality, which affected the velocity estimation results. Therefore, as the SCR increased, the RMSEs of all three points decreased. The RMSE of Point 3 was noticeably higher than that of Point 1 and Point 2, which was attributed to its lower scattering amplitude. When the SCR was above 15 $\mathrm{dB}$, the RMSEs of all three points were less than 0.01 $\mathrm{m}/\mathrm{s}$, which is comparable to the wavelength and adequate for high-precision micro-motion extraction. We excluded the results with SCRs below 8 dB, because the RMSEs were extremely large, and the results could not reflect the real target velocity.

4.2. Simulation Results of Extended Target Tracking

We considered a radar detection scenario involving two extended vessels. The motion parameters of the two targets are listed in

Table 1. The parameters of simulated extended targets.

Parameter Name	Target 1	Target 2
Translational velocity (m/s)	0	5
Center position (m)	5000	4955
$\mathsf{\Delta }x$ of scattering points (m)	(−10, 0, 20, 30)	(−15, 0, 15)
Amplitude of heaving (m)	0.9	0.9
Period of heaving (s)	15	5
Amplitude of pitching (°)	3	3
Period of pitching (s)	10	5

, while the theoretical velocities of the two targets are displayed in Figure 6. As is evident from the figure, the velocities of the scattering points of each target varied in a small range, which were difficult to discern from the Doppler-time profile with a decimeter-scale velocity resolution.

The MDCM was applied to calculate the conjugate phases of the scattering points. We extracted the conjugate phases under 10-pulse delay (${\varphi }_{10}$) and 160-pulse delay (${\varphi }_{160}$), then calculated the difference between the linear converted phases ${\widehat{\varphi }}_{160}$ and ${\varphi }_{160}$ as follows

$$\begin{array}{cc}\hfill \mathsf{\Delta }\varphi & =\left|{\widehat{\varphi }}_{160}-{\varphi }_{160}\right|\hfill \\ \hfill & =\left|\left(16{\varphi }_{10}-2k\pi \right)-{\varphi }_{160}\right|\hfill \end{array}$$

(27)

where k is the phase ambiguity integer. $\mathsf{\Delta }\varphi$ was employed for coincident point classification, and the selected points are demonstrated in red in Figure 7a. As is seen in the figure, the linear variation of the conjugate phases among different delayed intervals l provides a novel way for identifying coincident points.

Subsequently, the tracking of the scattering points was accomplished by utilizing the proposed LT method. The tracking result of two targets moving near to each other is shown in Figure 7b. The black lines correspond to the four scattering points of extended target 1 (1-1, 1-2, 1-3, 1-4), while the red lines represent the three scattering points of extended target 2 (2-1, 2-2, 2-3). As can be seen from the figure, the shapes of the tracks conform to the theoretical velocity curves (Figure 6).

The tracking results of the proposed method were compared under various scenarios. The translational velocities of target 1 and target 2 were denoted by $v_1$ and $v_2$, respectively. In this simulation, $v_1$ was set to zero, while $v_2$ took values of −2, 0, 2, and 5 $\mathrm{m}/\mathrm{s}$. The RMSEs of the simulated scenes are shown in

Table 2. The RMSEs of range and velocity tracks.

${\mathit{v}}_1$$(\mathbf{m}/\mathbf{s})$	${\mathit{v}}_2$$(\mathbf{m}/\mathbf{s})$	RMSEs of Range Tracks (m)	RMSEs of Velocity Tracks (m)
0	−2	0.470	0.023
0	0	0.590	0.019
0	2	0.765	0.020
0	5	0.277	0.016

. The RMSE values indicate that the estimate accuracies of the ranges and velocities did not significantly change when the targets were moving near or away from each other, even when the scattering points of the two targets were coincident in the range dimension (when $v_2=$ 2 $\mathrm{m}/\mathrm{s}$ or 5 $\mathrm{m}/\mathrm{s}$). These results demonstrate the effectiveness of the state correction process in the LT algorithm, which utilizes the motion correlation between scattering points and achieves the velocity estimation for the points without valid conjugate phase information.

The performance of the proposed MDCM-LT method was compared with the KF-ridge path regrouping (KF-RPRG) method [22] and the centroid group tracking (CGT) method [23]. The KF-RPRG method tracks individual scattering points and connects the interrupted tracks based on the gradient direction. The CGT method focuses on target centroids and classifies the target merging condition based on the number of measurements falling into the tracking gate of each target. The tracking results of the KF-RPRG and CGT are presented in Figure 8. In Figure 8a, a-b represents the b-th scattering points in the a-th extended target. It should be noted that the KF-RPRG does not differentiate which scattering point belongs to which extended target, and the black and red lines were only used for visual clarity in the plot. According to the figure, the KF-RPRG tracked the scattering points, but the fluctuations in the conjugate phases caused fluctuations in the movement trajectories. According to Figure 8b, the CGT offered stable tracking of the target centroids in both the range and velocity dimensions, but it was unable to track the micro-motion of each point.

The comparison between the three methods is listed in

Table 3. The comparison between the KF-RPRG method, the CGT method, and the proposed MDCM-LT method.

Item	KF-RPRG	CGT	The Proposed Method
Group tracking	×	✓	✓
Scattering point tracking	✓	×	✓
Real-time output	×	✓	✓
RMSE of velocity ($\mathrm{m}/\mathrm{s}$)	0.062	0.028 (center)	0.016
RMSE of distance ($\mathrm{m}$)	0.276	1.529 (center)	0.277
Processing time ($\mathrm{s}$)	3.26	1.07	3.06

. The regrouping process of the KF-RPRG could only be performed once the tracking was completed, which limited the real-time trajectory output. Due to the CGT only tracking the centroids of the extended targets, its processing time was the shortest among the three methods. The proposed MDCM-LT method utilized the micro-motion disparities among scattering points to rectify the conjugate phases, which led to reduced fluctuations in the trajectory tracking results. Furthermore, the RMSEs of the extracted velocity and distance trajectories were held within 0.016 $\mathrm{m}/\mathrm{s}$ and 0.277 $\mathrm{m}$, respectively, thereby indicating the effectiveness of the proposed method.

4.3. Experimental Results of Extended Target Tracking

A Ku-band radar was used for data acquisition in this experiment. The height of the radar placement was 60 $\mathrm{m}$, and the distance between the radar and the targets was approximately 5600 $\mathrm{m}$. In our experiment, extended target 1 was an anchored vessel. Extended target 2 consisted of two small ships connected by a rope, and the relative position between the two ships remained almost unchanged throughout the experiment. The experimental scene, the schematic diagram of target locations, and the range–time profile of the experimental radar data are illustrated in Figure 9.

Using delayed intervals of 6 and 10, the conjugate phases of the scattering points of two extended targets were extracted and used to classify the coincident points (as shown in Figure 10a). The LT method was then implemented for distance and conjugate phase tracking, with the tracking results presented in Figure 10b. In this figure, the black lines and red lines indicate the tracks of target 1 and target 2, while target 1 consists of three strong scattering points (1-1, 1-2, 1-3) and target 2 consists of two strong scattering points (2-1, 2-2). The conjugate phases were transformed into velocities, with the average velocities of target 1 and target 2 being −0.03 $\mathrm{m}/\mathrm{s}$ and 1.56 $\mathrm{m}/\mathrm{s}$, respectively. The average velocity of the vessel scattering points was nearly zero, which agrees with its anchored state. The average velocity of small ships is in accordance with the values calculated from the translation in the range domain. In addition, the small ships showed shorter motion periods and larger motion amplitudes than the vessel, which is attributed to their lower mass. Consequently, the proposed method enabled us to successfully track the micro-motion of the extended marine targets when the targets overlapped in the range dimension of the radar echo.

5. Conclusions

In this paper, we proposed a MDCM-LT algorithm for extracting the trajectories of rigid marine targets from radar echoes. Using the proposed method, we extracted the high-precision micro-motion trajectories of overlapping targets on the wavelength scale. The LT took advantage of the motion association within rigid extended targets to improve the tracking performance. Additionally, the MDCM method offers a novel approach to classify overlapping extended targets, even when they have similar translational velocities. This approach could potentially be used for target classification in complex marine conditions, indicating potential avenues for future research in this field.