Decentralized Approach for Translational Motion Estimation with Multistatic Inverse Synthetic Aperture Radar Systems

Abstract

This paper addresses the estimation of the target translational motion by using a multistatic Inverse Synthetic Aperture Radar (ISAR) system composed of an active radar sensor and multiple receiving-only devices. Particularly, a two-step decentralized technique is derived: the first step estimates specific signal parameters (i.e., Doppler frequency and Doppler rate) at the single-sensor level, while the second step exploits these estimated parameters to derive the target velocity and acceleration components. Specifically, the second step is organized in two stages: the former is for velocity estimation, while the latter is devoted to velocity estimation refinement if a constant velocity model motion can be regarded as acceptable, or to acceleration estimation if a constant velocity assumption does not apply. A proper decision criterion to select between the two motion models is also provided. A closed-form theoretical performance analysis is provided for the overall technique, which is then used to assess the achievable performance under different distributions of the radar sensors. Additionally, a comparison with a state-of-the-art centralized approach has been carried out considering computational burden and robustness. Finally, results obtained against experimental multisensory data are shown confirming the effectiveness of the proposed technique and supporting its practical application.

1. Introduction

As well known, radar images of man-made targets can be obtained by means of Inverse Synthetic Aperture Radar (ISAR) techniques, [1,2]. Compared to other methods, ISAR imaging has the advantage of generating all-day, all-weather maps of target reflectivity. This strength makes it an important system for target classification and recognition applications.

Nevertheless, conventional (i.e., single channel and single sensor) ISAR systems show the following inherent drawbacks: $\left(\mathrm{i}\right)$ conventional ISAR images provide a two-dimensional projection of targets having a 3D structure; $\left(\mathrm{ii}\right)$ achievable cross-range resolution depends on the intrinsic characteristics of the target motion and therefore, depending on the specific conditions, can be very poor; and $\left(\mathrm{iii}\right)$ the formation of an ISAR image requires the knowledge of the target motion parameters (both translation and rotation). Since generally, the targets are non-cooperative, these parameters must be estimated directly from the received signal.

Regarding the 3D target reconstruction, point ($\mathrm{i}$), interferometric ISAR (InISAR) systems offer a solution by applying interferometric technology to the ISAR images of multi-antenna systems. The traditional InISAR imaging usually employs three antennas to construct two perpendicular baselines [3]. However, in practical applications, it is not a simple matter to guarantee the complete orthogonality of the baselines due to the existence of factors such as system errors or environmental constraints. In [4], a practical maneuvering target 3D imaging algorithm based on the InISAR of an arbitrary three-antenna configuration is investigated.

The use of multistatic, or even Multiple Input Multiple Output (MIMO), ISAR systems with multiple sensors observing the same target with different acquisitions geometries and the joint exploitation of the acquired data allow much richer information to be garnered, overcoming limitations (ii) and (iii) imposed by conventional systems due to their single perspective. Over the last decade, researchers investigated the potentialities of multistatic/MIMO ISAR configurations for different tasks, see [5] and references therein.

In particular, distributed ISAR techniques have been also proposed to address the issues specified in $\left(\mathrm{ii}\right)$. For example, to enhance the cross-range resolution of rotating targets, multiple transmitters-receivers can synthesize a wider observation angle than the single sensor, bringing to a higher resolution [6,7,8]. Such approaches rely on a proper coherent combination of the received signals and, therefore, are effective in the case of limited angular separations among the individual sensors so that a stable behavior of the scattering mechanism can be assumed. In the case of a wider separation among the sensors, such an assumption does not hold. As the target e.m. response varies when observed by significantly different observation angles, incoherent combination approaches must be considered (e.g., [9] shows that enhanced target classification can be obtained by incoherent summation of single-sensor images acquired by spatially distributed sensors).

Finally, point $\left(\mathrm{iii}\right)$, even if motion compensation and image scaling are possible also with a conventional single sensor (limitations may be present for scaling), multistatic systems can be exploited for an enhanced estimation of the motion of the target. Specifically, in this respect, the main advantage of the multistatic systems compared to the monostatic (conventional single sensor) ones is the possibility to recover the complete target motion, namely:

• Velocity and acceleration ($x, y$) components for translation motion while in contrast, at the most, a single sensor could allow us to estimate the radial velocity and the modulus of the cross-range velocity (i.e., indeterminate sign).
• Roll, pitch, and yaw rates for rotation motion while only the overall effective rotation rate, or at the most, the vertical and horizontal rotation components could be estimated by single-sensor techniques.

It could be noticed that complete information regarding translation is of paramount importance for maritime surveillance purposes and, in general, estimated kinematic parameters could be exploited also for classification/recognition. A few contributions can be found in the literature for rotational motion estimation (relevant for image scaling, if needed) [10,11] and for the estimation of the target trajectory (relevant also for motion compensation) [12,13,14,15]. This work fits into the latter category, addressing the estimation of the target translational motion capitalizing on the spatial diversity offered by multistatic ISAR configurations.

In this framework, in general, two different data fusion strategies can be followed, namely centralized or decentralized. In the former, the kinematic parameters are obtained directly from the fusion of the multi-sensor signals. This can be obtained in two different ways as follows:

(a)

By fusing at image level, i.e., single-sensor images in a single multi-sensor image. Along this line, a preliminary proof of concept for the case of two platforms with constrained geometry, considering only the slow-time domain, was shown in [14,15].

(b)

By fusing at the cost function level, i.e., a combination of multiple cost functions, each evaluated at the single-sensor level. Fusion at the cost function level was considered in [12,13] for multistatic autofocus.

However, decentralized data fusion architectures are often preferred in multistatic radar systems because of their larger robustness and scalability under different operative conditions. Noticeably, unlike centralized procedures, they do not require wideband communication links between the sensors [16]. Therefore, as an alternative to centralized approaches, we propose here a decentralized technique to accomplish the translational motion estimation task in multiplatform imaging systems; some preliminary results along this line were previously described in [17,18]. The proposed decentralized multistatic technique estimates the kinematic parameters via a two-step procedure: (1) first estimate the single-sensor signal parameters (Doppler centroid and Doppler rate); (2) then jointly exploit the estimated signal parameters to estimate the kinematic of the target. The second step is organized into two stages: the first stage is for velocity estimation, while the second one is for velocity estimation refinement (if a constant velocity model motion can be regarded as acceptable or, for acceleration estimation, if a constant velocity assumption does not apply). The selection between the two motion models is ruled by a decision criterion properly defined. It could be noticed that a decentralized approach could be suitable as a stand-alone technique (providing the estimates of the target motion parameters) or to initialize centralized techniques (providing the initial guess of the target motion parameters that can then be refined by centralized estimation approaches). The theoretical performance is analytically derived in terms of the covariance matrices of the estimates parameters (i.e., the single-sensor target signal parameters, Doppler frequency and Doppler rate, and the target kinematic parameters, velocity, and acceleration); theoretical results are completed and confirmed by using synthetic data. Performance is then assessed against a number of varying conditions and compared to the performance of alternative approaches from the literature. Particularly, the impact of the spatial diversity among the sensors on the achievable performance is investigated. Regarding the comparison, the performance of the proposed technique is also compared to that achievable by a centralized approach at the cost function level. This alternative approach has been selected as it can be suitable for both co-located and widely separated sensors, as the proposed decentralized technique, whereas the centralized technique at the image level is suitable for the co-located case only. Finally, to further prove the effectiveness and demonstrate the practical applicability of the proposed techniques, results obtained by applying them to live multistatic ISAR data are also shown.

The remainder of the paper is organized as follows: in Section 2, the multi-sensor system geometry and the echo model are introduced; in Section 3, the proposed multistatic estimation techniques are presented and their performance is theoretically analyzed in Section 4 and assessed in Section 5; Section 6 shows the results achieved against live multistatic data sets; and Section 7 concludes the paper. Analytical details are reported in the appendices.

2. Geometry and Signal Model

In this work, we mainly refer to maritime scenarios. Figure 1 shows a pictorial view of the selected reference scenario comprising a coastal multistatic ISAR system and a vessel sailing in its field of view. However, it is worth underlining that the work could be easily generalized to cope with other scenarios and/or different targets.

Specifically, a formation of $N$ sensors is considered; sensor 1 is a monostatic active radar system transmitting and receiving while the remaining $N-1$ are receiving only devices. In this case, $N$ acquisitions are provided by the sensor network; one monostatic acquisition, from sensor 1, and $N-1$ bistatic acquisitions, each arising from the transmission from sensor 1 and reception from sensor $i \left(i=2,\dots ,N\right)$. The ship is modeled as a rigid body; therefore, as usual in ISAR literature, its motion is decomposed as the translation of an arbitrary reference point and the rotation of the body around that point (hereafter named target fulcrum). Figure 2 shows the considered acquisition geometry in a fulcrum-centered reference system.

In such a system, $r_{\mathit{i}}\left(\mathrm{t}\right)$ represents the position vector, changing with time, of sensor $i$, being ${\psi }_{\mathrm{i}}\left(t\right)$ and ${\zeta }_{\mathrm{i}}\left(t\right)$, respectively, which are the grazing angle and the angle between the projection on the $X–Y$ plane of the line-of-sight (LOS) and the y-axis (measured counter-clockwise from the y-axis). At image time (i.e., the synthetic aperture center here assumed as $t=0$ without loss of generality), the position of sensor $i$ is specified by:

$${\mathit{r}}_{\mathrm{i}}\left(0\right)=\left[\begin{array}{c}{x}_{\mathrm{i}}^0\\ y_{\mathrm{i}}^0\\ z_{\mathrm{i}}^0\end{array}\right]=r_{\mathrm{i}}^0\left[\begin{array}{c}\sin {\zeta }_{\mathrm{i}}^0\cos {\psi }_{\mathrm{i}}^0\\ -\cos {\zeta }_{\mathrm{i}}^0\cos {\psi }_{\mathrm{i}}^0\\ \sin {\psi }_{\mathrm{i}}^0\end{array}\right]$$

(1)

The target translation motion is taken into account by introducing the velocity vector $\mathit{v}$ and the acceleration vector $\mathit{a}$. As we are dealing with ship targets, the vertical component of both velocity and acceleration is assumed negligible so that we can write $\mathit{v}={\left[v_x v_y 0\right]}^T$ and $\mathit{a}={\left[a_x a_y 0\right]}^T$, however, the approach could be easily generalized to cope with situations where a vertical motion component also has to be included. As usually performed, the distance of the target fulcrum from sensor $i$ as a function of the slow-time $t$ can be approximated in second order as:

$$r_{\mathrm{i}}\left(t\right)\approx r_{\mathrm{i}}^0+{\dot{r}}_{\mathrm{i}}^{0}t+{\ddot{r}}_{\mathrm{i}}^0\frac{t^2}{2}$$

(2)

where ${\dot{r}}_{\mathrm{i}}^0$ and ${\ddot{r}}_{\mathrm{i}}^0$ are the distance first and second derivatives evaluated at $t=0$, and $\lambda$ is the wavelength. It is easy to show that, for the assumed geometry and target motion, these derivatives can be written as [19]:

$${\dot{r}}_i^0=\cos {\psi }_{\mathrm{i}}^0\left(-v_x\sin {\zeta }_{\mathrm{i}}^0+v_y\cos {\zeta }_{\mathrm{i}}^0\right)$$

(3)

$${\ddot{r}}_i^0=-\cos {\psi }_{\mathrm{i}}^0\left(a_x\sin {\zeta }_{\mathrm{i}}^0-a_y\cos {\zeta }_{\mathrm{i}}^0\right)+\frac{1}{ r_{\mathrm{i}}^0}{\sin }^2{\psi }_{\mathrm{i}}^0 \left(v_x^2+v_y^2\right)+\frac{1}{ r_{\mathrm{i}}^0}{\cos }^2{\psi }_{\mathrm{i}}^0 {\left(v_x\cos {\zeta }_{\mathrm{i}}^0+v_y\sin {\zeta }_{\mathrm{i}}^0\right)}^2$$

(4)

On this basis, the two-way Doppler centroid and Doppler rate at the i-th sensor are written as $f_i=-\left(\frac{{\dot{r}}_1^0}{\lambda }+\frac{{\dot{r}}_i^0}{\lambda }\right)$ and ${\dot{f}}_i=-\left(\frac{{\ddot{r}}_1^0}{\lambda }+\frac{{\ddot{r}}_i^0}{\lambda }\right)$, respectively.

3. Multistatic Translational Motion Estimation Technique

Considering the availability of a multistatic ISAR system, in this work, we address the estimation of the target translation motion parameters by means of a two-step procedure:

• Step 1: The first step is aimed at estimating the signal parameters (basically Doppler centroid and Doppler rate) at the single-sensor level;
• Step 2: The second step is aimed at estimating the target motion parameters by inverting their analytical relationship with the target signal parameters. This step comprises two possibilities as the specific analytical relation depends on the assumed model for the target motion (the choice between the two is driven by a proper target motion model selection criterion).

The complete decentralized processing scheme is shown in Figure 3. Details concerning the two steps and the adopted selection criterion are provided in the following sub-sections.

3.1. Single-Sensor Signal Parameters Estimation Technique

From the geometry and signal model introduced in Section 2, it is apparent that, depending on the values of the integration time ($T$) of the vessel-specific motion, the transmitted signal bandwidth ($B$), and the relative position of the sensor with respect to the target, both the range and Doppler cell migration might be observed in the acquired single-sensor raw data; particularly, range migration is highly dependent from the Doppler centroid (range walk component), while Doppler migration is basically due to the Doppler rate. Therefore, these two single-sensor signal parameters, Doppler centroid and Doppler rate, are thus estimated as those values providing the single-sensor target image with the highest quality. In this work, as common in the ISAR literature, we resort to contrast maximization [20] but also other cost functions could be adopted (for example, entropy minimization [21]).

Since, for the application under consideration, the target size is limited to say about one hundred meters and we are not dealing with very high resolution, we consider here a processing technique that exactly compensates the migration through the range and Doppler cells for the ship fulcrum and then it images the target via Fourier transform (i.e., rectangular format). Specifically, starting from single-sensor data in the range-compressed, slow-time domain, the range migration correction is performed by compensating at slow time instant ($t$) for a fast-time delay given by $\frac{v_{r}t}{c}$, where $v_r$ is the generic radial velocity, and $c$ is the speed of light. By expressing the generic radial velocity in terms of the generic value for the Doppler centroid ($\tilde{f}),$ such a delay can be rewritten as $\frac{\tilde{f}\lambda t}{c}=-\frac{\tilde{f}t}{f_c}$, with $f_c$ being the carrier frequency. Therefore, the correction is performed in the fast-frequency and slow-time domain by multiplying the transformed data by:

$${\Phi }_1\left(f_{\tau },t,\tilde{f} \right)=exp\left(-j2\pi f_{\tau }\frac{\tilde{f}}{f_c}t\right)$$

(5)

with $f_{\tau }$ being the fast-frequency. The following step is a range inverse Fourier transform to map back the data in the range-compressed and slow-time domain, where they are multiplied with a reference chirp signal (i.e., data dechirping):

$${\Phi }_2\left(t,\tilde{\dot{f}} \right)=exp\left(-j\pi \tilde{\dot{f}}{t}^2\right)$$

(6)

with $\tilde{\dot{f}}$ being the generic value for the Doppler rate. Finally, an azimuth Fourier transform is performed to obtain the data in the image domain. Let $I_i\left(\tilde{f},\tilde{\dot{f}}\right)$ be the intensity value of this $i$-th single-sensor complex image by defining the image contrast as:

$$I{C}_i\left(\tilde{f},\tilde{\dot{f}}\right)=\frac{\sqrt{E\left\{{\left[I_i\left(\tilde{f},\tilde{\dot{f}}\right)-E\left(I_i\left(\tilde{f},\tilde{\dot{f}}\right)\right)\right]}^2\right\}}}{E\left(I_i\left(\tilde{f},\tilde{\dot{f}}\right)\right)}$$

(7)

where the operator $E\{·\}$ represents the image spatial mean. The final Doppler and Doppler rate estimates at sensor $i$ are obtained by maximizing the contrast.

$$\left({\widehat{f}}_i,{\widehat{\dot{f}}}_i\right)=argma{x}_{\tilde{f},\tilde{\dot{f}}}\left\{I{C}_i\left(\tilde{f},\tilde{\dot{f}}\right)\right\}$$

(8)

3.2. Kinematic Parameters Estimation Technique

The second step is organized into two stages: the first stage is for velocity estimation, while the second one is for velocity estimation refinement. If a constant velocity model motion can be regarded as acceptable, or for acceleration estimation, a constant velocity assumption does not apply. The decision between the two motion models is performed according to the criterion described in Section 3.3.

3.2.1. Kinematic Parameters Estimation Technique—Stage 1

Starting from the first stage, considering the formation of $N$ sensors and using a matrix notation and results in Section 2, we can write:

$$\mathit{f}=\mathit{G}\mathit{v}$$

(9)

where $\widehat{\mathit{f}}$ is the $N\times 1$ vector collecting the $N$ estimated Doppler centroids, $\mathit{v}={\left[v_x v_y\right]}^T$, and $\mathit{G}$ is a $N\times 2$ matrix taking into account the specific acquisition geometry:

$$\mathit{G}=\frac{1}{\lambda }\left[\begin{array}{cc}2\cos {\psi }_1^0\sin {\zeta }_1^0& -2\cos {\psi }_1^0\cos {\zeta }_1^0\\ \cos {\psi }_1^0\sin {\zeta }_1^0+\cos {\psi }_2^0\sin {\zeta }_2^0& -\cos {\psi }_1^0\cos {\zeta }_1^0-\cos {\psi }_2^0\cos {\zeta }_2^0\\ \dots & \dots \\ \cos {\psi }_1^0\sin {\zeta }_1^0+\cos {\psi }_N^0\sin {\zeta }_N^0& -\cos {\psi }_1^0\cos {\zeta }_1^0-\cos {\psi }_N^0\cos {\zeta }_N^0\end{array}\right]$$

(10)

Exploiting the Doppler frequencies estimated in step 1, and assuming the knowledge of the acquisition geometry, the target velocity is obtained as the least square (LS) solution of the system in Equation (9), thus obtaining:

$${\widehat{\mathit{v}}}_{\mathit{f}}={\left({\mathit{G}}^T\mathit{G}\right)}^{-1}{\mathit{G}}^T\widehat{\mathit{f}}={\mathit{G}}^{\#}\widehat{\mathit{f}}$$

(11)

where ${(\ast )}^{\#}$ denotes the pseudo-inverse operator.

3.2.2. Kinematic Parameters Estimation Technique—Stage 2

As stated above, in the second stage, if a constant velocity motion model can be regarded as acceptable, the velocity estimated at the previous step is refined by incorporating the Doppler rate information, otherwise, the velocity estimated at stage 1 is retained as the final estimate, and the Doppler rate information is exploited to estimate the acceleration.

Starting from the case of constant velocity, the exploitation of the Doppler rate measurements requires the solution of a system of $N$ non-linear equations in two unknowns:

$$\{\begin{array}{c}{\widehat{\dot{f}}}_i=-\frac{2}{r_1^0\lambda }\left[{\mathit{sin}}^2{\psi }_1^0\left(v_x^2+v_y^2\right)+{\mathit{cos}}^2{\psi }_1^0{\left(v_x\mathit{cos}{\zeta }_1^0+v_y\mathit{sin}{\zeta }_1^0\right)}^2\right]\\ ⋮\\ {\widehat{\dot{f}}}_i=-\frac{1}{r_1^0\lambda }\left[{\mathit{sin}}^2{\psi }_1^0\left(v_x^2+v_y^2\right)+{\mathit{cos}}^2{\psi }_1^0{\left(v_x\mathit{cos}{\zeta }_1^0+v_y\mathit{sin}{\zeta }_1^0\right)}^2\right]\\ \hspace{1em}\hspace{1em}-\frac{1}{r_i^0\lambda }\left[{\mathit{sin}}^2{\psi }_i^0\left(v_x^2+v_y^2\right)+{\mathit{cos}}^2{\psi }_i^0{\left(v_x\mathit{cos}{\zeta }_i^0+v_y\mathit{sin}{\zeta }_i^0\right)}^2\right]\\ ⋮\\ {\widehat{\dot{f}}}_N=-\frac{1}{r_1^0\lambda }\left[{\mathit{sin}}^2{\psi }_1^0\left(v_x^2+v_y^2\right)+{\mathit{cos}}^2{\psi }_1^0{\left(v_x\mathit{cos}{\zeta }_1^0+v_y\mathit{sin}{\zeta }_1^0\right)}^2\right]\\ \hspace{1em}\hspace{1em}-\frac{1}{r_N^0\lambda }\left[{\mathit{sin}}^2{\psi }_N^0\left(v_x^2+v_y^2\right)+{\mathit{cos}}^2{\psi }_N^0{\left(v_x\mathit{cos}{\zeta }_N^0+v_y\mathit{sin}{\zeta }_N^0\right)}^2\right]\end{array}$$

(12)

The system in Equation (12) can be linearized by using a first-order Taylor series approximation around a generic target tentative velocity, ${\mathit{\vartheta }}^0={\left[v_x^0 v_y^0\right]}^T$:

$$\widehat{\dot{\mathit{f}}}-\dot{\mathit{f}}\left({\mathit{\vartheta }}^0\right)=\mathit{H}·\left(\mathit{\vartheta }-{\mathit{\vartheta }}^0\right)$$

(13)

where $\widehat{\dot{\mathit{f}}}$ is the $N\times 1$ vector collecting the $N$ estimated Doppler rates, $\dot{\mathit{f}}\left({\mathit{\vartheta }}^0\right)$ is the $N\times 1$ vector collecting the $N$ Doppler rates evaluated for the tentative velocity, $\mathit{\vartheta }={\left[v_x v_y\right]}^T$, and a matrix notation has been adopted based on the following definition:

$$\mathit{H}={\left[\begin{array}{cc}\frac{\partial \dot{f_1}}{\partial v_x}& \frac{\partial \dot{f_1}}{\partial v_y}\\ ⋮& ⋮\\ \frac{\partial \dot{f_N}}{\partial v_x}& \frac{\partial \dot{f_N}}{\partial v_y}\end{array} \right]}_{|\vartheta ={\vartheta }^0}$$

(14)

The initial tentative value in this second stage could be chosen equal to the estimate provided by the first stage (i.e., initial ${\mathit{\vartheta }}^0={\widehat{\mathit{v}}}_{\mathit{f}}$), and the overall system to be solved is obtained by augmenting Equation (9) with the one in Equation (13), namely:

$$\widehat{\dot{\mathit{F}}}-\dot{\mathit{F}}\left({\mathit{\vartheta }}^0\right)=\left[\begin{array}{c}{\mathbf{0}}_{N\times 1}\\ \widehat{\dot{\mathit{f}}}-\dot{\mathit{f}}\left({\mathit{\vartheta }}^0\right)\end{array}\right]=\mathit{U}·\left(\mathit{\vartheta }-{\mathit{\vartheta }}^0\right)$$

(15)

where ${\mathbf{0}}_{N\times 1}$ is a vector of $N\times 1$ with all elements set to 0 and $\mathit{U}$ is a $2N\times 2$ block matrix equal to $\mathit{U}={\left[{\mathit{G}}^T{\mathit{H}}^T\right]}^T$. It is worth underlining that this joint use allows us to refine the estimate of the target velocity using both measured Doppler centroids and rates. The LS solution is thus given by:

$$\widehat{\mathit{\vartheta }}={\mathit{\vartheta }}^0+{\mathit{U}}^{\#}\left(\dot{\mathit{F}}\left(\mathit{\vartheta }\right)-\dot{\mathit{F}}\left({\mathit{\vartheta }}^0\right)\right)$$

(16)

Then, through Equation (16), the target kinematic parameters are updated with respect to the tentative value (${\mathit{\vartheta }}^0)$, and the procedure is reiterated until, at the generic iteration, the displacement $\sqrt{\delta v_x^2+\delta v_y^2}$ is within the requirements on the velocity accuracy or the maximum admitted number of iterations is reached. In the hypothesis of error-free measurements, this algorithm converges to the true target velocity components.

In case a constant velocity assumption does not apply, the second stage is instead devoted to the estimation of the acceleration, again by exploiting the Doppler rate information estimated in the first step. Using previous results and definitions and adopting, again, a matrix notation, we can write:

$$\widehat{\dot{\mathit{f}}}-\dot{\mathit{f}}\left({\widehat{\mathit{v}}}_{\mathit{f}}\right)=\mathit{G}·\mathit{a}$$

(17)

where $\mathit{a}={\left[a_x,a_y\right]}^T$. The target acceleration is then estimated as:

$$\widehat{\mathit{a}}={\mathit{G}}^{\#}·\left(\dot{\mathit{f}}-{\dot{\mathit{f}}}_v\left({\widehat{\mathit{v}}}_{\mathit{f}}\right)\right)$$

(18)

3.3. Automatic Motion Model Selection Criterion

The automatic selection of the motion model is based on a comparison between the Doppler rate measurements and the Doppler rate values evaluated, assuming the velocity is equal to that estimated at stage 1 and null acceleration. Similar values imply the acceptable hypothesis of null acceleration so that the constant velocity motion model is selected. The other schema is selected in the opposite case. The similarity is assessed by means of the Mahalanobis distance by evaluating the distance between the error vector, collecting the differences of the Doppler rate values measured and theoretically evaluated, and a Normal distribution characterized by zero mean value and proper covariance matrix. Specifically, let $\mathsf{\Delta }\dot{\mathit{f} }$ be the error Doppler rate vector defined as:

$$\mathsf{\Delta }\dot{\mathit{f} }=\widehat{\dot{\mathit{f}}}-\dot{\mathit{f}}\left({\widehat{\mathit{v}}}_{\mathit{f}},\mathit{a}=\mathbf{0}\right)$$

(19)

This error vector has zero mean. Additionally, since the Doppler rate measurements and ${\widehat{\mathit{v}}}_{\mathit{f}}$ are independent variables, it is possible to write the covariance matrix of the Doppler rate error vector as follows:

$${\Sigma }_{\mathsf{\Delta }\dot{\mathit{f}}}={\Sigma }_{\widehat{\dot{\mathit{f}}}}+{\Sigma }_{\dot{\mathit{f}}\left({\widehat{\mathit{v}}}_{\mathit{f}},\mathit{a}=\mathbf{0}\right) }$$

(20)

where the two matrix components are specified in the following Equations (23) and (25) for ${\Sigma }_{\widehat{\dot{\mathit{f}}}}$, and Equation (B8) for ${\Sigma }_{\dot{\mathit{f}}\left({\widehat{\mathit{v}}}_{\mathit{f}},\mathit{a}=\mathbf{0}\right) }$. For decision-making purposes, the Mahalanobis distance is compared to a threshold $K$, i.e.,

$$\mathsf{\Delta }{\dot{\mathit{f}}}^{\mathit{T}}{\left({\Sigma }_{\mathsf{\Delta }\dot{\mathit{f}}}\right)}^{-1}\mathsf{\Delta }\dot{\mathit{f}} \frac{<}{>} K$$

(21)

and the constant velocity model is selected in case the distance is below the decision threshold. The K value is set in order to achieve a fixed probability value that, in case of actual constant velocity, the error Doppler rate vector lies within the corresponding ellipsoid. In Section 4.3, it is explained in detail how to set the value of K.

4. Theoretical Performance Analysis

In this section, the theoretical performance of the proposed technique is analytically derived. Furthermore, theoretical derivations are validated by comparison with simulated analysis. Performance achievable at the single-sensor level in the estimate of the signal parameters is considered first (Section 3.1). Following this, the performance in the estimate of the target kinematics parameters (Section 3.2) and the impact of the motion model selection (Section 3.3) are addressed.

4.1. Theoretical Performance Analysis—Step 1: Single-Sensor Signal Parameters

The performance in the estimate of the signal parameters is derived under the following assumptions: (i) statistical independence among the sensors; (ii) target described by a single dominant scatterer globally, taking into account the overall target Radar Cross Section (RCS); and (iii) negligible rotation motion. Particularly, assumption (ii) allows a closed-form performance derivation, and assumption (iii) allows to derive a benchmark for the achievable performance. These two limitations are then removed in section V where performance is assessed for an extended target and in the presence of 3D rotation motions.

Starting from the Doppler frequency measurements collected in the $N\times 1$ vector $\widehat{\mathit{f}}$, it can be demonstrated (see Appendix A) that the error of the estimates $e_{{\widehat{f}}_i}$ can be modeled as a zero-mean Gaussian random variable, i.e., $e_{{\widehat{f}}_i}~\mathcal{N}\left(0,{\sigma }_{{\widehat{f}}_i}^2\right)$, with a standard deviation equal to:

$${\mathsf{\sigma }}_{{\widehat{f}}_i}=\frac{72}{{\left(2\pi \right)}^2{T}^{2}SN{R}_{in{t}_i}{\left(\frac{B}{f_c}\right)}^2}=\frac{6}{{\left(2\pi \right)}^2{T}^{2}SN{R}_{in{t}_i}}·\frac{12}{{\left(\frac{B}{f_c}\right)}^2}$$

(22)

where $T$ is the aperture time, $B$ is the bandwidth of the transmitted signal, $f_c$ is the carrier frequency, and S$N{R}_{in{t}_i}$ is the integrated signal-to-noise ratio (SNR) at sensor $i$ (i.e., SNR evaluated in a single-sensor image domain). Interestingly, it can be noticed that Equation (22) can be decoupled as the product of two terms, the first one $\left(\frac{6}{{\left(2\pi \right)}^2{T}^{2}SN{R}_{in{t}_i}}\right)$ representing the Cramer Rao bound for the performance achievable in the estimate of the frequency of complex signals having constant amplitude and polynomial phase, [22], and the second one $\left(12/\left(\raisebox{1ex}{$B$}\!\left/ \!\raisebox{-1ex}{$f_c$}\right.\right){ }^2\right)$ representing the impact on the performance of the estimation performed via optimization of the range migration correction.

Moving to the Doppler rate measurements collected in the $N\times 1$ vector $\widehat{\dot{\mathit{f}}}$, the performance can be easily evaluated by recalling that the contrast optimization technique can achieve the Cramer Rao bound in selected conditions (i.e., when a single scatterer is exploited or when the multiple exploited scatterers share a similar SNR value) and in any case shows an efficiency very close to one, [23]. On this basis, the error in the estimate of the Doppler rate $e_{{\widehat{\dot{f}}}_{\mathrm{i}}}$ can be modeled as a zero-mean Gaussian random variable, i.e., $e_{{\widehat{\dot{f}}}_{\mathrm{i}}}~\mathcal{N}\left(0,{\mathsf{\sigma }}_{{\widehat{\dot{f}}}_{\mathrm{i}}}^2\right)$, with a standard deviation equal to (Appendix A)

$${\mathsf{\sigma }}_{{\widehat{\dot{f}}}_i}=\frac{1}{\pi }\sqrt{\frac{90}{SN{R}_{in{t}_i}{T}^4}}$$

(23)

Furthermore, since the measurements at the different sensors are independent, the Doppler and Doppler rate covariance matrix are, respectively:

$${\Sigma }_{\widehat{\mathit{f}}}=diag\left({{\mathsf{\sigma }}_{{\widehat{f}}_1}}^2,\dots ,{{\mathsf{\sigma }}_{{\widehat{f}}_N}}^2\right)$$

(24)

$${\Sigma }_{\widehat{\dot{\mathit{f}}}}=diag\left({{\mathsf{\sigma }}_{{\widehat{\dot{f}}}_1}}^2,\dots ,{{\mathsf{\sigma }}_{{\widehat{\dot{f}}}_N}}^2\right)$$

(25)

and due to the decoupling between Doppler frequency and Doppler rate [22], the covariance matrix of the measurements as a whole has a diagonal block structure and can be easily derived from the two above matrices.

To verify the theoretical derivations, a comparison against synthetic data is shown here. For this purpose, we define the following reference scenario used from now on (any modification will be duly mentioned): a formation of two sensors is considered with distances $r_1^0=r_2^0=10$ km and aspect and grazing angles equal to ${\zeta }_1^0=2°$, ${\zeta }_2^0=-6°$, ${\psi }_1^0={\psi }_2^0=0°$. In agreement with the study cases presented in [14,15], the active system is assumed to transmit a bandwidth of $B=100 \mathrm{MHz}$ with a center frequency of f_c = 10 GHz, a Pulse Repetition Frequency of $PRF=600 \mathrm{Hz}$, and 2048 pulses in the Coherent Processing Interval (CPI). The target moves with translational motion according to velocity $\mathit{v}={\left[8 4 0\right]}^T \mathrm{m}/\mathrm{s}$, negligible acceleration, and a constant yaw rotation motion ${\mathsf{\omega }}_y=1 deg/s$. The SNR is set to 33 dB (image domain).

Figure 4 shows the obtained results for 1000 independent trials. Particularly, the red and green ‘×’ markers represent the Doppler and Doppler rate measurements at sensor 1 and sensor 2, respectively. For comparison, the figure shows also the theoretical ellipse as achievable from Equations (22)–(25), and the ellipse derived from the measurements, both associated with a probability value equal to 0.99. A very good agreement is observed between theoretical and simulated results for both measured signal parameters, thus validating the analytical results.

4.2. Theoretical Performance Analysis for Step 2

4.2.1. Step 2—Stage 1: Velocity Estimation

The covariance matrix of the velocity estimated at stage 1 can be derived by combining Equations (11) and (24):

$${\mathsf{\Sigma }}_{{\widehat{\mathit{v}}}_{\mathit{f}}}={\mathit{G}}^{\#}{\mathsf{\Sigma }}_{\widehat{\mathit{f}}}{\mathit{G}}^{{\#}^T}$$

(26)

Again, Figure 5 shows the estimated velocity components for 1000 independent trials and the same study case in Figure 4. Particularly, the blue markers represent the velocities as estimated at stage 1, the continuous black ellipse is the theoretical iso-level curve, and the blue one is the ellipse derived from the simulations, both referred to 0.99 probability. The statistics (the mean value and the standard deviation in $x$ and $y$-directions, respectively) are shown in

Table 1. Mean and standard deviation—Stage 1 velocity.

Statistics	Target Kinematic Parameters	Theoretical Values	Decentralized Approach Estimates
mean	$v_x \left[m/s\right]$	8.000	8.001
	$v_y \left[m/s\right]$	4.000	3.999
std	$v_x \left[m/s\right]$	0.269	0.289
	$v_y \left[m/s\right]$	0.009	0.010

. Firstly, we observe, once again, a very good agreement between the theoretical predictions and the simulated results. Moreover, we notice how the accuracy in the estimation is higher for the y-component than for the x-component. This is a consequence of the exploitation of the Doppler frequency measurements that, in the considered case study (limited angular diversity between the two sensors), are highly dependent on $v_y$, but scarcely sensible to $v_x$.

4.2.2. Step 2—Stage 2: Velocity Refinement

The covariance matrix of the velocity estimated at stage 2 can be proven equal to (see Appendix B):

$${\Sigma }_{\widehat{\mathit{\vartheta }}\left(k\right)}=\mathit{A}\left(k\right)\left[{\mathbf{0}}_{2N\times 2}{\left[{\mathbf{0}}_{NxN} {\Sigma }_{\widehat{\dot{\mathit{f}}}}\right]}^T\right]\mathit{A}{\left(k\right)}^T+B{\mathsf{\Sigma }}_{{\widehat{\mathit{v}}}_{\mathit{f}}}{\mathit{B}}^T$$

(27)

where matrices $\mathit{A}$ and $\mathit{B}$ are defined as:

$$\mathit{A}\left(k\right)=\sum _{i=0}^{k-1}{\left(\mathit{I}-{{\mathit{U}}^{\#}}_{|\mathit{\vartheta }={\mathit{\vartheta }}^0}\left[\begin{array}{c}{\mathbf{0}}_{N\times 2}\\ {\mathit{H}}_{|\mathit{\vartheta }={\mathit{\vartheta }}^0}\end{array}\right]\right)}^i{{\mathit{U}}^{\#}}_{|\mathit{\vartheta }={\mathit{\vartheta }}^0}$$

(28)

$$\mathit{B}=\left(\mathit{I}-{{\mathit{U}}^{\#}}_{|\mathit{\vartheta }={\mathit{\vartheta }}^0}\left[\begin{array}{c}{\mathbf{0}}_{N\times 2}\\ {\mathit{H}}_{|\mathit{\vartheta }={\mathit{\vartheta }}^0}\end{array}\right]\right)$$

(29)

$\mathit{I}$ is the identity matrix, $k$ is the iteration at which the refinement is stopped, and the initial velocity (${\mathit{\vartheta }}_{\mathbf{0}}$) is the value provided by Equation (11). In order to analyze the impact of the refinement process on the estimation accuracy, Figure 6 shows the area of the theoretical 1σ ellipse associated with the covariance matrix in Equation (27) as a function of the iteration number for three geometry configurations: (i) narrow angular separation in blue (${\zeta }_1^0=2°$ and ${\zeta }_2^0=-6°$), corresponding to the case previously discussed, (ii) medium angular separation in red (${\zeta }_1^0=6°$ and ${\zeta }_2^0=-18°$), and (iii) wide angular separation in green (${\zeta }_1^0=12°$ and ${\zeta }_2^0=-36°$).

From the figure, for all configurations, it can be observed a rapid improvement in the first iterations until a minimum value (identified in the figure with a triangle and named $k_{min}$ in the following) is reached. After this value, it is observed a slight increase. This characteristic is used to define a further stop criterion: based on Equation (27), the refinement process finishes after $k_{min}$ iterations.

The red crosses in Figure 5 represent the refined estimates corresponding to points coming from stage 1 (blue crosses). These results have been obtained by stopping the iterative procedure at $k_{min}$.

Table 2. Mean and standard deviation of the estimated stage 2 velocities–Point-like target.

Statistics	Target Kinematic Parameters	Theoretical Values	Decentralized Approach Estimates
mean	$v_x \left[m/s\right]$	8.000	8.001
	$v_y \left[m/s\right]$	4.000	3.999
std	$v_x \left[m/s\right]$	0.038	0.040
	$v_y \left[m/s\right]$	0.009	0.010

compares the theoretical mean value and standard deviation of the stage 2 estimated velocity components to the values achieved by simulations. Again, the simulated results are well in line with the theoretical predictions. As expected, the exploitation of the Doppler rate measurements (highly sensitive to the velocity x-component) greatly increases the performance with respect to stage 1 for the x-component, whereas almost unvaried performance is obtained for the y-component. Overall, the joint exploitation of Doppler centroid and Doppler rate measurements allows the achievement of considerably high accuracy for both components.

4.2.3. Step 2—Stage 2: Acceleration Estimation

The covariance matrix of the acceleration estimated at stage 2 can be proven equal to (see Appendix B):

$${\Sigma }_{\mathit{a}}={\mathit{G}}^{\#}\left({\Sigma }_{\widehat{\dot{\mathit{f}}}}+{\Sigma }_{\dot{\mathit{f}}\left({\widehat{\mathit{v}}}_{\mathit{f}},\mathit{a}=\mathbf{0}\right) } \right){\mathit{G}}^{{\#}^T}$$

(30)

with ${\Sigma }_{\dot{\mathit{f}}\left({\widehat{\mathit{v}}}_{\mathit{f}},\mathit{a}=\mathbf{0}\right) }=\mathit{C}{ \mathsf{\Sigma }}_{{\widehat{\mathit{v}}}_{\mathit{f}}}{\mathit{C}}^T$ and $\mathit{C}$ being the $N\times 2$ matrix defined in the Appendix. In this case, to compare the theoretical predictions to the simulations, the previous reference scenario is maintained, but an acceleration $\mathit{a}={\left[0.45 0.225 0\right]}^T \mathrm{m}/{\mathrm{s}}^2$ is also included. For this case, Figure 7 shows the estimated velocity (Figure 7a) and acceleration (Figure 7b) components with the superimposed corresponding ellipses, while

Table 3. Mean and standard deviation of the estimated velocity and acceleration–Point-like target.

Statistics	Target Kinematic Parameters	Theoretical Values	Decentralized Approach Estimates
mean	$v_x \left[m/s\right]$	8.000	7.986
	$v_y \left[m/s\right]$	4.000	4.001
	$a_x \left[m/s^2\right]$	0.450	0.450
	$a_y \left[m/s^2\right]$	0.225	0.225
std	$v_x \left[m/s\right]$	0.269	0.289
	$v_y \left[m/s\right]$	0.009	0.010
	$a_x \left[m/s^2\right]$	0.002	0.002
	$a_y \left[m/s^2\right]$	0.000	0.000

compares the theoretical and simulated mean values and standard deviations. The shown results once again confirm the validity and correctness of the theoretical derivations. Moreover, we can appreciate the high accuracy achievable by the proposed approach in the estimate of the acceleration components. This very good result is enabled by the exploitation of the Doppler rate measurements that are directly proportional to the acceleration showing high sensitivity with respect to its variations.

4.3. Automatic Motion Model Selection

As stated in Section 3.3, the motion model is selected by means of the decision rule in Equation (21). The constant velocity model is selected in case the distance is below the decision threshold, otherwise, an acceleration is estimated. The $K$ value is set in order to achieve a fixed probability value that, in case of actual constant velocity, such a condition is correctly detected so that the velocity undergoes the refinement foreseen at stage 2. In particular, when the target is actually uniformly moving, the distance in Equation (21) is a random variable following a chi-square probability density function with $N$ degrees of freedom. Therefore, the probability of correctly detecting this condition is provided by the chi-square cumulative distribution function:

$$P_d\left(K, N\right)=\frac{\gamma \left(N/2, K/2\right)}{\mathsf{\Gamma }\left(N/2\right)}$$

(31)

where $\gamma \left(\ast , \ast \right)$ is the lower incomplete gamma function and $\mathsf{\Gamma }\left(N/2\right)$ is the gamma function, [24]. For the special case of two sensors $\left(N=2\right)$, the $P_d$ simplifies as follows:

$$P_d\left(K, N=2\right)=1-e^{-\frac{K}{2}}$$

(32)

The threshold value ($K$) has to be set in order to guarantee a reasonably high $P_d$ value, while at the same time keeping a good sensitivity to the presence of small acceleration values. In order to prove the effectiveness of the proposed approach, for the study case comprising two sensors, Figure 8 shows the probability of selecting an accelerated model for a target moving with a velocity equal to [8 4 0] m/s and variable acceleration $\left(a_x,a_y\right)$ having set $K$ so that $P_d=0.99$. Figure 9 shows the probability of deciding on a uniform motion as a function of the acceleration in (a) x-direction and (b) y-direction for different $P_d\left(K\right)$ values. As evident from the results, the proposed approach is quite effective at discriminating between the two conditions; particularly in Figure 8 where the presence of an acceleration is correctly detected even in the presence of low values. Concerning the transition between the two decision regions, as expected, the higher the $P_d$ value, the larger the transition region is as it can be appreciated by inspecting Figure 9 which shows the cuts of Figure 8 along the $a_x$ and $a_y$ directions around $\mathit{a}=0$ m/s². Nevertheless, it is easy to observe that even in the presence of $P_d$ values close to 1 (i.e., 0.99), the decision test maintains a very good performance, assuring a high sensitivity.

5. Performance Assessment

5.1. Performance Assessment with Respect to SNR Conditions and Spatial Diversity

The performance for different SNR conditions is first investigated. Specifically, Figure 10a shows the standard deviation of stage 1 velocity in the x-direction (red curve) and the y-direction (blue curve) as a function of integrated (i.e., integrated over both fast and slow time and over the target scatterers) SNR values. When a constant velocity model is acceptable, the refined velocity of stage 2 exhibits the accuracy illustrated by the dashed curves in Figure 10a. Conversely, if a constant acceleration model is selected, the estimated acceleration demonstrates the performance indicated by the curves in Figure 10b.

As expected, higher integrated SNR values lead to improved estimation accuracy. Nevertheless, even for low/medium values of the SNR (we have to recall that the considered values refer to an integrated SNR and the ship has been already detected, so 20 dB is quite a low value for imaging purposes), good performance is achieved also by the stage 1 velocity estimation (that represents the initial value when the refinement in stage 2 is applied).

The performance of different sensor positions is now investigated. Particularly, for the study cases previously introduced, Figure 11a,b refers to the uniform motion case and shows the theoretical standard deviation of the estimated velocity as a function of the angular separation in aspect between the two sensors for different separations in grazing. Both stage 1 and stage 2 are shown. Figure 11c regards instead the performance in estimating the accelerations for the non-uniform case.

Starting from the case of uniform motion, from the curves shown in Figure 11a,b, it is possible to make the following observations: (i) in general the performances improve as the angular diversity in aspect increases and remain unchanged as the grazing diversity varies. This is a consequence of having assumed the target motion on the (x, y) plane (actually, the angular diversity in grazing would play a fundamental role in the estimation of the velocity along z); (ii) the improvement is marked in particular for stage 1, which, by exploiting only the Doppler information, fails to adequately estimate the cross velocity (x-direction) in the case of small diversity; and (iii) the estimate of the component in y reaches minimum standard deviation when both sensors observe the same radial component. In the case of accelerated motion, the performance improves as the diversity increases for both components. Consistently with the results shown above, the performance in estimating the component in y is in any case better than the estimation of the component in x as a consequence of the direct impact of component in y on the radial acceleration. In any case, even in the presence of limited angular diversity, good performance is obtained.

5.2. Performance Assessment with Respect to Extended Targets

The assumption of a target described by a single dominant scatterer globally, taking into account the overall target RCS, is now removed and the analysis is generalized considering the case of an extended ship target characterized by a number of dominant scatterers. The considered target model is shown in Figure 12. As in the previous analysis, the SNR is set to 33 dB (image domain) and the same motion conditions have been considered. Particularly, concerning SNR, to allow a comparison with the results for a point-like target, the contribution from each scatterer is scaled so that the integration among the different scatterers would result, again, in the same SNR value. This is obtained by considering the single scatterer reflectivity equal to a fraction (1/20) of the point-like reflectivity.

For the first tested case (uniform motion, velocity $v={\left[8 4 0\right]}^T \mathrm{m}/\mathrm{s}$, and negligible acceleration), Figure 13 shows the obtained results for 1000 independent trials. As in the previous analyses, the blue crosses represent the estimated velocities at stage 1, while the red ones represent the estimates after the refinement at stage 2. For comparison, the same figure shows also, for both stages, the theoretical iso-level curves corresponding to 0.99 probability level and the corresponding ellipses containing 99% of the estimates. In addition,

Table 4. Mean and standard deviation of the estimated velocity–Extended target and uniform motion.

Statistics	Target Kinematic Parameters	Stage 1 Velocity			Stage 2 Velocity
		Theoretical Values (Point like Target)	Simulated (Extended Target with Constant Yaw Rotation)	Simulated (Extended Target with 3D Rotation)	Theoretical Values (Point like Target)	Simulated (Extended Target with Constant Yaw Rotation)	Simulated (Extended Target with 3D Rotation)
mean	$v_x \left[m/s\right]$	8.000	8.009	7.689	8.000	8.001	8.047
	$v_y \left[m/s\right]$	4.000	3.999	4.033	4.000	3.999	4.033
std	$v_x \left[m/s\right]$	0.269	0.416	0.403	0.038	0.056	0.071
	$v_y \left[m/s\right]$	0.009	0.014	0.014	0.009	0.014	0.014

compares the theoretical (point-like target) and simulated (extended target) mean value and standard deviation for both velocity components and both stages. From the shown results, it is easy to observe that, when moving from the ideal point-like to the extended target case (keeping fixed the overall SNR), the estimates are still unbiased but a slight degradation in terms of error standard deviation is observed. This is a consequence of more dispersed estimates of the Doppler frequency and Doppler rate at the single-sensor level, due to the interferences and also to the approximated procedure used in distributing the overall RCS among the scatterers. A similar degradation is observed on both the components.

For the second tested case (velocity $v={\left[8 4 0\right]}^T \mathrm{m}/\mathrm{s}$ and acceleration $a={\left[0.45 0.225 0\right]}^T \mathrm{m}/{\mathrm{s}}^2$), Figure 14 shows the obtained results for 1000 independent trials. Particularly, Figure 14a represents the velocity components as estimated by the first stage, while Figure 14b shows the acceleration components as estimated using the first stage velocity and exploiting the Doppler rate information. Additionally, the theoretical ellipses are superimposed on the figures to compare the point-like and extended target performance.

Table 5. Mean and standard deviation of the estimated velocity and acceleration–Extended target and accelerated motion.

Statistics	Target Kinematic Parameters	Theoretical Values	Decentralized Approach Estimates (Extended Target with Constant Yaw Rotation)	Decentralized Approach Estimates (Extended Target with 3D Rotation)
mean	$v_x \left[m/s\right]$	8.000	8.007	7.611
	$v_y \left[m/s\right]$	4.000	3.999	4.025
	$a_x \left[m/s^2\right]$	0.450	0.451	0.448
	$a_y \left[m/s^2\right]$	0.225	0.225	0.226
std	$v_x \left[m/s\right]$	0.269	0.373	0.397
	$v_y \left[m/s\right]$	0.009	0.013	0.014
	$a_x \left[m/s^2\right]$	0.002	0.002	0.003
	$a_y \left[m/s^2\right]$	0.000	0.001	0.001

shows the estimated statistics. From the shown results, also in this case, a small degradation is observed and similar considerations as those already performed for the uniform case can be applied.

Finally, the impact of the presence of a 3D rotation motion is analyzed. For this purpose, again, a simulated analysis has been carried out by also including sinusoidal pitch and roll motion with amplitudes of 0.25° and 1.25° and frequencies of 0.178 Hz and 0.091 Hz, respectively, [25]. The obtained results are reported in Table 4 for the case of uniform motion and in Table 5 for the case of accelerated motion. From the shown results, we can observe that the inclusion of a three-dimensional rotation introduces a slight polarization in the velocity and acceleration estimates, while the standard deviations remain almost unchanged. Anyway, despite the small degradation, in both cases we observe that very good performance is maintained.

5.3. Performance Assessment with Respect to Centralized Approach

The performance of the proposed decentralized technique is now compared to that of a centralized approach at the cost function level that can be derived along the same line in [12,13]. In this case, the fusion of the multi-sensor signals is implemented at the cost-function level. Specifically, the single-sensor image contrast values are combined in a new multistatic cost function defined as the product of the individual cost functions, as illustrated in Figure 15. This centralized approach has been selected for performance comparison being, as the proposed one, suitable without restrictions with respect to the angular separation among the sensors.

The optimization problem to be solved for the multistatic case is therefore:

$$\widehat{\mathit{v}},\widehat{\mathit{a}}=argma{x}_{\mathit{v},\mathit{a}}\left\{\prod _{i=1}^{N}I{C}_i\left(\mathit{v},\mathit{a}\right)\right\}$$

(33)

where $I{C}_i$ denotes the contrast of the image at the $i$-th sensor when motion is compensated according to the generic ($\mathit{v},\mathit{a}$).

A drawback of this approach resides in the increased dimension of the space over which optimization has to be carried out. In our case, four parameters have to be estimated (two velocity and two acceleration components), thus resulting in a 4D optimization (2D when a constant velocity model is acceptable and 6D in a more general case, if also the vertical direction has to be included). This can represent a severe issue since the multistatic cost function can be expected to be multimodal, which generally represents a problem for the optimization algorithms to converge on the true solution.

To analyze the performance of the centralized and decentralized approaches, results against synthetic data are shown here. The same scenario explained in Section 4.1 is considered by testing two angular separations: (a) narrow angular separation (NAS), with aspect and grazing angles equal to ${\zeta }_1^0=2°$, ${\zeta }_2^0=-6°$, ${\psi }_1^0={\psi }_2^0=0°$; (b) wide angular separation (WAS), with aspect and grazing angles equal to ${\zeta }_1^0=10°$, ${\zeta }_2^0=-30°$, ${\psi }_1^0={\psi }_2^0=0°$.

The extended target is assumed to be moving again with velocity being $v={\left[8 4 0\right]}^T \mathrm{m}/\mathrm{s}$ and acceleration being $a={\left[0.45 0.225 0\right]}^T \mathrm{m}/{\mathrm{s}}^2$. The multistatic cost function maximization in the centralized approach is implemented using the deterministic algorithm Nelder–Mead [26]. This algorithm requires an initial point as input and it is not guaranteed to converge to a global minimum/maximum. To deal with this problem, first, a range of possible velocity and acceleration values are defined. In our case, $\left|v_x\right|=\left|v_y\right|<15 \left[\mathrm{m}/\mathrm{s}\right]$ and $\left|a_x\right|=\left|a_y\right|<0.5 \left[\mathrm{m}/{\mathrm{s}}^2\right]$. Then, random velocity and acceleration values are generated within the allowed range and used as input parameters by the Nelder–Mead algorithm to maximize the multistatic cost function in Equation (33), obtaining the estimated target kinematic parameters as a result. The process is repeated for 200 randomly generated points. The final estimated values are selected as those that generate the highest value of the contrast product. Additionally, to also study the bound of performance achievable by means of the centralized approach, the case of an initial point equal to the true value of velocity and acceleration is also considered.

Table 6. Comparison between centralized and decentralized approaches–Estimates statistics.

Statistics	Target Kinematic Parameters	Narrow Angular Separation-NAS			Wide Angular Separation-WAS
		Centralized Approach (Random Initial Points)	Centralized Approach (Real Value Initial Point)	Decentralized Approach	Centralized Approach (Random Initial Points)	Centralized Approach (Real Value Initial Point)	Decentralized Approach
mean	$v_x \left[m/s\right]$	8.007	8.007	8.007	7.953	7.950	7.950
	$v_y \left[m/s\right]$	3.999	3.999	3.999	4.007	4.007	4.007
	$a_x \left[m/s^2\right]$	0.451	0.451	0.451	0.450	0.450	0.450
	$a_y \left[m/s^2\right]$	0.225	0.225	0.225	0.225	0.225	0.225
std	$v_x \left[m/s\right]$	0.376	0.376	0.374	0.815	0.078	0.078
	$v_y \left[m/s\right]$	0.013	0.013	0.013	0.143	0.014	0.014
	$a_x \left[m/s^2\right]$	0.002	0.002	0.002	0.021	0.001	0.001
	$a_y \left[m/s^2\right]$	0.001	0.001	0.001	0.005	0.000	0.000

compares the mean value and the standard deviation of the estimated values for the decentralized technique, centralized technique with randomly selected initial point, and centralized with initial point coincident with the true value.

From the results for both velocity and acceleration, the following observations apply: (i) in all considered conditions, the decentralized technique achieves the same performance as the centralized one initialized with the true motion parameters; (ii) in WAS conditions, the decentralized approach outperforms the centralized technique with randomly selected initial values (this likely happens when the optimization over the 4D space, despite the high number of random initial points, does not converge to the good solution). These results confirm the effectiveness and robustness of the proposed approach which also has lower requirements from an implementation point of view (N 2D optimizations instead of a 4D optimization over an N factor function).

6. Experimental Results

In order to demonstrate the practical effectiveness of the proposed technique and validate the previous theoretical performance analysis, we apply both the decentralized and centralized techniques to multistatic ISAR data acquired in an experimental campaign carried out in an anechoic chamber at the SELEX Galileo facility in Caselle (Turin, Italy) [11].

The setup is based on the use of a single reflector Compact Range System which generates a planar wave front in the test zone. This includes the parabolic reflector, the feed system, and the positioner of the target under tests (TUT). The reflector (Figure 16a) is an offset parabolic reflector P/N 5755 made by Scientific Atlanta that, when illuminated by a spherical wave front transmitted by an antenna located into its focus (feed), generates a planar wave front on the TUT. Two HP 83622As are used for the transmitter and receiver signal sources, while the measurement instrumentation is based on an HP 8510C Network Analyzer. Chamber walls, ceilings, and floors are covered with absorber material (pyramid and wedge) to minimize unwanted reflections and diffractions. In this manner, a stray radiation level within the test zone below –35 dB is ensured.

The system transmits a stepped frequency waveform in the Ku-band ($16.5 \mathrm{GHz}$) with a fixed pulse-to-pulse frequency increment of $\mathsf{\Delta }f=3.75 \mathrm{MHz}$, resulting in an overall bandwidth of $3 \mathrm{GHz}$. The second receiver is located $60 \mathrm{cm}$ from the transmitter, resulting in a bistatic channel having an aspect angle separation of $\mathsf{\Delta }\zeta =4.3°$ with respect to the monostatic link. The turntable rotation yields an angular separation burst-to-burst of $\delta \theta =0.07°$. The considered TUT is an ATR42 aircraft model (shown in Figure 16b) as representative of an extended target.

An overall illumination angle and bandwidth of $\mathsf{\Delta }\mathsf{\theta }\simeq 5.18°$ and $B=1.5 \mathrm{GHz}$ have been selected, respectively, resulting in a cross-range resolution (${\rho }_{cr}=\frac{\mathsf{\lambda }}{2\mathsf{\Delta }\mathsf{\theta }}=10 \mathrm{cm}$) and a range resolution (${\rho }_r=c/2B=10 \mathrm{cm}$). Referring to Figure 2, ${\zeta }_1^0=\frac{1}{4}\mathsf{\Delta }\zeta =1.075°$, ${\zeta }_2^0=-\frac{3}{4}\mathsf{\Delta }\zeta =-3.225°$ and ${\psi }_1^0={\psi }_2^0=0°$. The experimental facility produces only a yaw rotational motion. Therefore, an emulated translational motion is superimposed on the acquired data according to velocity ($\mathit{v}=\left[8 1 0\right] \mathrm{m}/\mathrm{s}$) and acceleration ($\mathit{a}=\left[0.45 0 0\right] \mathrm{m}/{\mathrm{s}}^2$). Furthermore, appropriately scaled noise is applied to achieve an integrated SNR of $33 dB$. Figure 17 shows the resulting images. The defocusing effect, due to uncompensated translation motion, is quite visible.

These data are provided in input to the processing chain in Figure 3.

Table 7. Mean and standard deviation of Doppler and Doppler rate estimates–Experimental data.

Statistics	Target Parameters	Theoretical Values	Decentralized Approach Estimates
mean	${\widehat{f}}_1 \left[Hz\right]$	−93.535	−93.565
	${\widehat{f}}_2 \left[Hz\right]$	−126.489	−126.544
	$\widehat{{\dot{f}}_1} \left[Hz/s\right]$	0.221	0.211
	$\widehat{{\dot{f}}_2} \left[Hz/s\right]$	−1.629	−1.641
std	${\widehat{f}}_1 \left[Hz\right]$	0.064	0.083
	${\widehat{f}}_2 \left[Hz\right]$	0.064	0.078
	$\widehat{{\dot{f}}_1} \left[Hz/s\right]$	0.003	0.003
	$\widehat{{\dot{f}}_2} \left[Hz/s\right]$	0.003	0.003

shows the mean and standard deviation of the estimates of the Doppler centroid and of the Doppler rate obtained at step 1 of the decentralized approach for 500 independent trails for both sensors and compares them to the corresponding theoretical performance. It is easy to observe a very good agreement between the theoretical predictions and the results obtained against the experimental data.

Figure 18 shows the estimated kinematic parameters. The blue (Figure 18a) and green (Figure 18b) crosses represent the velocity and acceleration components estimated by applying the decentralized approach, respectively. Additionally, the same parameters are estimated for the same 500 independent trails, but by applying the selected centralized approach with the optimization started from the actual velocity and acceleration values. The corresponding results are represented by red (Figure 18a) and violet (Figure 18b) triangles.

Table 8. Mean and standard deviation of estimates, aircraft target.

Statistics	Target Kinematic Parameters	Aircraft Target
		Theoretical Values	Centralized Approach (Real Value Initial Point)	Decentralized Approach
mean	$v_x \left[m/s\right]$	8.000	8.006	8.006
	$v_y \left[m/s\right]$	1.000	1.000	1.000
	$a_x \left[m/s^2\right]$	0.450	0.450	0.450
	$a_y \left[m/s^2\right]$	0.000	9 × 10−5	9 × 10−5
std	$v_x \left[m/s\right]$	0.022	0.028	0.028
	$v_y \left[m/s\right]$	4 × 10−4	5 × 10−4	5 × 10−4
	$a_x \left[m/s^2\right]$	9 × 10−4	1 × 10−3	1 × 10−3
	$a_y \left[m/s^2\right]$	4 × 10−5	5 × 10−5	5 × 10−5

shows the mean and standard deviation as predicted by the theory and as obtained from the experimental data via the decentralized and centralized approach. We can observe that the achieved performance against the experimental data is well in line with the theoretical expectations. Additionally, in agreement with the results presented in Section 5.3, the proposed decentralized approach provides the same performance as that provided by the centralized fed with the actual kinematic parameters. Qualitatively similar results could be shown also for different (and in particular worst) SNR conditions (not reported here for the sake of compactness). These results further validate the effectiveness of the proposed technique and support its applicability in practical environments. Finally, Figure 19 shows an example of target images after motion compensation by using the estimated kinematic parameters, where the aircraft shape can be easily recognized in both the monostatic and bistatic ISAR images.

7. Conclusions

In this work, multi-sensor translational motion estimation techniques have been devised. A network formed by an active sensor and some passive devices has been considered and a decentralized two-step scheme proposed where the single-sensor, signal parameters are first estimated, and then this information is used in the second step to estimate the kinematic of the target. Particularly, Doppler measurements are used to make an initial velocity estimate, and when a constant velocity model is acceptable, Doppler rate information is exploited to provide a more accurate velocity estimate. If instead, via a proper decision rule, it is assessed that an acceleration has to be taken into account, such a decentralized approach estimates the acceleration by exploiting the initial velocity jointly with the Doppler rate information.

The performance of the proposed technique has been first analytically derived and then analyzed and assessed with respect to different conditions and parameters. Additionally, the proposed technique was compared with a centralized approach. The obtained results indicate that the decentralized technique provides the same performance as the centralized one operating in ideal conditions (i.e., the initial point for the optimization coinciding with the actual values of the parameters), and at the same time, lowering the complexity of the optimization problem. Finally, the techniques have been tested against experimental multi-sensor datasets and the results prove the suitability of the proposed approach for practical applications.