Micro-Motion Estimation of Maritime Targets Using Pixel Tracking in Cosmo-Skymed Synthetic Aperture Radar Data—An Operative Assessment

: In this paper, we propose a novel strategy to estimate the micro-motion (m-m) of ships from synthetic aperture radar (SAR) images. To this end, observe that the problem of motion and m-m detection of targets is usually solved using synthetic aperture radar (SAR) along-track interferometry through two radars spatially separated by a baseline along the azimuth direction. The approach proposed in this paper for m-m estimation of ships, occupying thousands of pixels, processes the information generated during the coregistration of several re-synthesized time-domain and not overlapped Doppler sub-apertures. Speciﬁcally, the SAR products are generated by splitting the raw data according to a temporally small baseline using one single wide-band staring spotlight (ST) SAR image. The predominant vibrational modes of different ships are then estimated. The performance analysis is conducted on one ST SAR image recorded by COSMO-SkyMed satellite system. Finally, the newly proposed approach paves the way for application to the surveillance of land-based industry activities.


Introduction
Synthetic aperture radar (SAR) provides high-resolution images of static ground scenes, whereas processing of data containing ground object motion results in a variation of focused areas. A special case of such motion is vibration where its pattern may be distinctly identifiable in focused SAR intensity images, as well as in a precise amplitude analysis. The X-band SAR can be well suited to image vibration because its wavelength is close to typical vibration amplitudes. This work extends the investigation technique observed in [1] where the pixel tracking has been used for the first time to detect small movements of marine targets. Before this work, several interesting research studies have been performed in the field of micro-motion (m-m) estimation. In [2] the effects of rotation and vibration in millimeter-waveforms have been assessed by processing simulated and real SAR data. Vibrating targets cause the modulation of the azimuth phase history of a SAR image. This phenomenon may be considered as a micro-Doppler (m-D) frequency modulation. In [3] the oscillation frequency and amplitude were computed from the time-frequency distributions and the findings were in accordance with the respective ground truths. In [4] a vibrating m-D signature for a bistatic SAR system with

Estimation Procedure and Processing Architecture
This section consists of three parts aimed at guiding the reader towards the design of the estimation procedure. More precisely, the first part provides a description of the underlying vibrational model from the physical point of view, whereas the second part explains the ideas behind the estimation procedure. Finally, the third part presents the entire processing chain and describes each operation performed in the course of the parameter estimation.

Computational Model
In this subsection, we present the model used to describe the vibrations generated by a ship. Since the whole ship can be considered as a complex body consisting of several interacting mechanical components, the analysis of its behavior represents a formidable task. For this reason, at the design stage, the vibration analysis is carried out by resorting to finite element techniques. In fact, a common practice consists in first dividing the ship into five parts: Then, the vibration contribution of each of these components are analyzed. The forces involved along the 3-dimensional reference system can be seen in Figure 1a. The forces (F1, F2, F3) are generated by the rotation movement of the propulsion propeller, amplified by the sail effect generated by the rudder. The forces (F4, F5) are instead due to the movement of the main axis of the transmission of motion generated by the machines. The coupling system and the engine revolutions-per-minute (rpm) reduction box generates the force F6 while the diesel engine will pull the forces (F7, F8). These are the main sources of vibration, generated by a ship, that propagate along all the structures and superstructures that make up the hull. Figure 1b represents the same ship observed in SAR coordinates. The point P1 is hit by the electromagnetic pulses and, during the image formation process, is projected on the point P2 located on the layover projection plane constituting the SAR image. If P1 vibrates along the three Euclidean dimensions, P2 exhibits a vibration in the range-azimuth SAR plane located into the new vector-space slant coordinates. Standard SAR-processing methods are based upon the assumption of a static scene. If targets are moving, their positions in the SAR-image are shifted in azimuth and range-azimuth plan and a defocusing may occur. Figure 2 shows the SAR acquisition geometry of the target P1 having a displacement of velocity (∆x, ∆y) with respect to the sensor. According to Figure 2 the azimuth SAR resolution is determined by the SAR aperture length L, the electromagnetic wavelength λ and the angle between the velocity vector and line of sight θ: The Doppler frequency shift of the stationary target is given by: where V t is the platform velocity and θ is the angle existing between the velocity vector and the target direction located at distance R 0 . On the other hand, the Doppler shift defines the azimuth position in the SAR image located on the following angle: If the target P1 has a velocity component v R along the range direction, it exhibits the following Doppler shift: In the course of the focusing process, this frequency error is transformed into an angular error denoted by δθ. Let us consider the differential version of (2) (5) and equal this variation to (4), then the differential of θ is given by:  This parameter produces an azimuth target miss-location due to a range velocity component with respect to the radar platform. The Doppler phase variation from a point reflector can be observed through the output of the matched filter whose expression is: where R(t) is the range variation due to the SAR antenna motion along its aperture. According to Figure 2, if the sensor is flying following the path r 1 (t) = (V t , 0, H) and the target is moving along r 2 = (x 0 + ∆x, y 0 + ∆y, 0) where ∆x = v x t + a x t 2 /2 and ∆y = v y t + a y t 2 /2 with (a x , a y ) ∈ R, it is possible to approximate the relative range variation R(t) = |r 1 (t) − r 2 (t)| by that of a stationary reflector modified by a correction term due to the target movement. If V t t << R 0 , then Replacing ∆x and ∆y with the respective expressions in (8) and neglecting the higher order terms, R(t) can be expressed as: where and a R = a x cos θ + a y cos γ sin φ + The target velocity components of (12) introduce an azimuth displacement as described by (6). The velocities and accelerations, a x , a y , v x , v y , in Formula (13), represent the contribution of the target to global acceleration and provide the variation of the Doppler frequency in the received signal. This result can be viewed as a smearing effect along azimuth and polarized in the same direction as the target's speed. Thus, the existence of nonzero target speeds and accelerations leads to a quality degradation of the image contrast since the signal to clutter ratio decreases. This degradation effect is quantified through the following equation: According to (7), where R(t) = a R t 2 2 , the corresponding phase error δφ(t) = 2πδa R t 2 λ returns a defocusing and the reduction of the target contrast in the SAR image. The resulting frequency is given by the phase derivative f (t) = 1 2π δφ(t) dt = 2δa R t λ and covers the following bandwidth over the time interval T = L/V t : Now, from the differentiation of (2), δ f d = −2δθV t sin θ λ and imposing that δ f d = B C , the azimuth resolution of the moving target becomes δθ b = δa R L v 2 sin θ . Using Equation (1), we obtain: By combining (16) and (14), the degraded azimuth resolution of moving targets can be estimated. The blurring significantly increases for fine-resolution mapping and highly squinted angles of observation. There is also a range resolution decrease due to the target range migration during the time of measurement. In the case of a non-accelerating target, the parameter V t is dominant and according to δa R = 2V t v x /R 0 and Formulas (16) and (1) yields: Recalling that the L V t = T, the target defocusing δR b = δθ b R 0 can be expressed as: which relates the defocusing length to the movement of the target during the time of measurement. Summarizing, moving targets produce significant position errors in the azimuth positioning due to the target range velocity component. Moreover, it leads to smearing effects on the focused signals in both range and cross-range directions. A pure displacement along the azimuth direction returns a pure defocusing along the same direction. All these effects are detected by two-dimensional normalized cross-correlation and tracked during the staring spotlight orbital time by the sub-pixel coregistration process. In the next subsection, this process is explained in detail.

Estimation Procedure
Sub-pixel offset tracking (SPOT) is a relevant technique to measure large-scale ground displacements in both range and azimuth directions. The technique is complementary to differential interferometric SAR and persistent scatterers interferometry when the radar phase information is unstable [16,17]. In this paper, we apply the above pixel tracking technique to a single spotlight image instead of multi-temporal interferometric images. Specifically, the single observation is divided into Doppler sub-apertures in order to investigate the fastest displacements of moving targets. Note that the acquisition duration is of the order of a few hundredths of a second.
With the above remarks in mind, we focus on maritime applications and estimate the displacement generated by the marine targets detected using a single SAR image focused in Doppler sub-apertures. In accordance with the temporal subdivision strategy of Figure 3, we observe the offset trend by computing the normalized cross-correlation once the image is partitioned into small patches. The estimation procedure consists in shifting the master for each temporal event and calculating the correlation between adjacent Doppler sub-apertures, according to a small-temporal baseline strategy. In order to provide a more operational approach to the technique we used the software Sarproz, implemented by Prof. D. Perissin. The software is versatile because it provides the use of parallel computing and allows the possibility of designing all the coregistration parameters.
For the proposed estimation procedure, we set the cross-correlation window size to 128 × 128 pixels and the oversampling factor to 64 in both the range and azimuth directions. All the coregistration parameters are reported in Table 2. The offset components of the sub-pixel normalized cross-correlation, according to [18,19] are described by the complex parameter D tot (range, azimuth) which is estimated by the following equation: Specifically, in (19) D displ is the offset component (18) generated by the earth displacement; D topo is the offset component generated by the earth displacement when located on highly sloped terrain; D orbit is the offset caused by residual errors of the satellite orbits; D control is the offset component generated by general attitude and control errors of the flying satellite trajectory; D atmosphere and D noise are the contributions generated by the electromagnetic aberrations due to atmosphere parameters, space and time variations, and general disturbances due to thermal and quantization noise, respectively. The atmospheric time-variation during the very short acquisition time interval has little influence on the temporal component of the last displacement parameters because of its low accuracy. All errors are compensated for by choosing only high energy and stable points and subtracting the initial offsets in order to retrieve the shift contributions only generated by the target displacement. As the concluding remark of this subsection, we note that the limits of the proposed procedure are mainly due to the maximum frequency of vibration observed and the minimum sensitivity in sensing vibration amplitudes. The first parameter is directly related to the sampling frequency. In fact, it is equal to the sampling frequency divided by two. For the specific case considered in the next section, 200 sub-apertures are produced over an acquisition time of 12 s. If we refer to this configuration, the time baseline is equal to 0.06 s. Thus, the maximum observation frequency is approximately 33 Hz, this means that it is possible to observe the vibrations generated by an endothermic engine running at 33-60 rpm. As for the minimum sensitivity, it is directly related to the spatial resolution of the sensor. In respect of the specific sensor used to assess the procedure, namely the SAR sensor of a CSK satellite, the spatial resolution is approximately one meter. In order to finalize a precise coregistration, we choose an over-sampling factor equal to 200. This value, in principle, allows dividing the single pixel range-azimuth into 200 × 200 sub-pixels and, thanks to this, the minimum estimated displacement was equal to five millimeters.

Computational Architecture Description
In this subsection, we describe in detail the processing chain used to perform the estimation of m-m. To this end, the block-scheme of the architecture is shown in Figure 4a. The input of the processing architecture is a single raw dataset (block 1). The computational block 2 consists of a temporal n-stage splitter that is responsible for generating n temporal sub-apertures. In the numerical examples, we generate 200 sub-products. The duration of the full-band Doppler Spotlight acquisition is about 12 s. This means that the observation time increment between two contiguous sub-apertures is approximately 0.06 s. The aim of block 3 is to focus all the n sub-apertures. The final results are a time series of SAR images in the SLC configuration having lower resolution in azimuth. This resolution loss is inversely proportional to the lost Doppler band fraction. Note that SAR products generated by block 3 are not yet coregistered, this operation is implemented in block 5, which performs the sub-pixel coregistration whose setting parameters are generated by block 4. Specifically, the SAR acquisition characteristics summary and the coregistration settings are shown in Tables 1 and 2 respectively. Moreover, the master of the multiple coregistration processes is selected according to the strategy referred to in the following as "small-temporal baseline" with the sliding master. The output of the co-recorder is represented by block 6. The result is a stack of maps of shifts in range and azimuth. Now, blocks 7 and 8 perform pixel by pixel analysis to estimate the oscillations of the ships. Figure 4b represents a more detailed explanation of computational blocks 7 and 8. The input of block 7 consists of a single column of shifts, that for the case at hand, has a dimensionality of 1 × 200. The vector of shifts is complex and after the coregistration process the displacements are estimated along both the range and azimuth direction. The computational block 7.1 interpolates the input data by a given factor. Next, block 7.

Experimental Results
In this section, we provide some experimental results by analyzing the vibrational modes of two ships at anchor. Moreover, we analyze in depth the frequency response and the trend over the time of the vibrations by selecting seven measurement points for the first study case and two measurement points for the second one. All points are distributed widely on the outer deck of the ships. The results also include the estimation of the distribution of the resonance frequencies distributed in the space where a comparative analysis of two vibration profiles extracted both along the keel and along the ship's deck is also carried out. Detailed information about the processed SAR image are summarized in Table 3.  Figure 3 shows the spotlight SAR acquisition geometry. The SAR observation takes about 12 s using approximately 25 kHz of the Doppler band. The whole image formation history was divided into 200 temporal azimuth sub-apertures of duration about 0.3871 s, according to a small temporal baselines strategy. Figure 5 is a schematic representation of the parameters estimated by the coregistration procedure. The square number one is a focused pixel of the master image and the square number two is the same pixel but located on the slave image. The parameters d and θ are the distance between the master and slave pixel centers and the angle respect to the horizontal axis respectively. Figure 6a depicts the magnitude of the ROI of the ship under test. Figure 6b represents the estimated vibration energy, consisting in parameter d of Figure 5 integrated through time. Figure 7a shows the average phase spanned in time of the pixel displacement. This parameter consists of the angle θ indicated in Figure 5, whereas Figure 7b contains the measurement points map of the ship.   It is important to observe that the estimation of slow-rate variation velocities is performed using the SAR investigation technique called along-track interferometry (ATI). Unfortunately, this configuration requires the use of two or more radars separated by an along-track-baseline of variable length from a few meters up to a few tens of meters. The longer the baseline, the greater the sensitivity of the Doppler variations of this bistatic SAR system. It is clear that in order to obtain a high sensitivity to vibrations, it is necessary to design high baselines ATI geometries. However, there are no satellite systems configured in this way and therefore it is necessary to find other solutions [1]. Decomposing the raw SAR data into multiple Doppler sub-apertures may not be a fully successful solution. This is because investigating the raw data generated by Earth observation satellites, which are usually never oversampled [1], will never be possible to estimate correlated ATI interferograms. In fact, in Figure 8a,b the module and phase of infra-chromatic coherence are visible. Inspection of the figures confirms that no information can be extrapolated. Another technique to estimate the m-m consists in evaluating the anomalies of the Doppler centroid during the focusing stage. This technique is also not immune to problems because a reliable estimation might occur if raw SAR images are available. However, such data are not always easy to obtain at this level of processing. Moreover, this solution is suitable for estimating vibrations only in the azimuth dimension. The proposed technique allows the estimation of the displacements, even if of lower intensity, in both range and azimuth directions, thanks to the highly accurate estimation of the shifts made by the co-registrator. As stated before, we consider seven measurement points and for each of them the vibration spectrum and the temporal trend of movements in the time domain are displayed in Figures 9-15. The location of each point is summarized in Table 4.    Figure 10a,b. Figure 10a shows the trend of displacements over time. In this case, the amplitude of oscillation is rather low compared to the previous case with the maximum equal to 0.2. Figure 10b shows the vibrational frequency spectrum. In this case, the spectral lines are well defined with a maximum at about 34 Hz. The results for measurement point number 3 are shown in Figure 11a,b. Figure 11a shows the trend of displacements over time. On this measurement point, the amplitude of oscillation is also lower with respect to vibrations observed on the measurement point 1 where an average value is located to 0.25 starting from 2 s to approximately 11 s. Figure 11b shows the vibrational frequency spectrum. In this case, the spectral lines are well defined with a maximum around 23 Hz and 36 Hz. The measurement point number 4 is the most particular because it is chosen to study the vibrations present at the end of the front master of the ship. As Figure 12a shows, the oscillations over time have numerous peaks at the maximum amplitude of 1 unit. The preponderance of these peaks is very stable at 25 Hz, as shown by Figure 12b, which represents the spectrum. The results concerning the time domain vibration pattern of measurement points 5 and 6 are reported in Figures 13a and 14a, respectively. The figures show a similarity with mean values around 0.4. The corresponding spectra are shown in Figures 13b and 14b, confirming resonance oscillations of about 12 and 42 Hz for point 5 and 32 Hz at measurement point 6. No significant oscillations are measured at this point. Figure 15a shows the time course of the displacement observed on the measurement point number 7 for which negligible amplitude values are measured (on average about 0.0030 compared to the maximum of 1 unit). The spectrum has a very low maximum frequency, located around 2 Hz (result showed in Figure 15b). A global view of the resonance frequencies distributed over all the pixels composing the ship is shown in Figure 16. The figure shows two reflectivity profiles marked with the numbers 1 and 2. The first reflectivity profile is shown in Figure 17a. This orientation was chosen in order to study the trend of keel vibrations. This representation in amplitude is given by the plus markers while the frequency value is represented by times markers. It is considered very important to underline that both physical values have been normalized to 1 so as to show the trend on the same graph. Finally, Figure 17b represents the trend of the resonance frequency (plus points) due to the vibration energy (times symbols). The observations are measured along with the vibrational profile 2. It is found that an increase in frequency corresponds to an increase in vibration energy, both on the keel (profile number 1) and on the transverse profile (profile number 2).

Study Case 2
In this second case, the object under investigation is a merchant ship specially chosen with many lateral lobes and disturbances due to defocusing effects that are markedly more visible in the Doppler domain. This worst case allows us to test the robustness of the estimation algorithm. Figure 18 is the SLC in magnitude of the ship. The optical image of the vessel concerning the second case study is contained in Figure 19a. The considered measurement points are two as summarized in Table 5.  Figure 19b represents the energy map of vibrations where a high energy level can be observed at the measurement point 1 and coinciding with the chimney stack. This structure clearly vibrates more than the rest of the ship and this anomaly is correctly detected. Figure 20a

Discussion and Future Assessments
The results obtained in this work have proved that the use of amplitude information alone from SAR Spotlight images is enough to estimate the m-m generated by ships. As a matter of fact, the exploitation of high-resolution SAR Spotlight data has allowed us to obtain displacement maps with millimetric precision and, hence, to appreciate the m-m. The vibrations have been estimated in the time and frequency domain considering two different ships. In the first domain, an analysis from the energy point of view can be performed, while in the second domain the main oscillation modes can be identified at frequencies ranging from a few Hz to a few tens of Hz. More importantly, it is very likely that the observed oscillation modes are due to the machines used to generate the electrical energy on board, namely generated by the diesel alternators. As a matter of fact, the propulsion system of large ships can be generally designed following two ways. In the first case, the propeller is connected directly to the endothermic engine. This is a case which is usually adopted on oil tankers, equipped with very large and slow engines. The resulting frequencies are very low, quantifiable at no more than 10 Hz. The other approach consists in installing a speed reducer between the diesel endothermic engine and the propeller shaft. This solution is widely used because it allows the engine to run at higher rates with respect to the rotation frequency of the propeller shaft. The resulting frequencies will necessarily be higher with respect to the previous situation and varying from approximately 10 to few tens of Hz. Thus, considering the estimated results, we are confident about the reliability of the estimation algorithms and we have observed ships belonging to the second category.
In addition, given the high quality of the vibration profiles, we believe that this work opens the way to the design of high-performance ship detection, classification, and recognition methods. For instance, the vibrational map of coded targets allows establishment of whether or not the ship engines are turned on, or to locate the position of the engine room within the naval silhouette. Moreover, in principle, it could be also possible to recognize the number and type of engines installed on board through the deep-learning approach. Finally, the herein developed approach can be further extended to estimate the vibrations generated by terrestrial industrial installations, coal-fired, and nuclear power plants.

Materials and Methods
In this section, we describe the software tools and the hardware used in the performance analysis in order to provide the detailed guidelines for the replication of the experiments. Specifically, the authors used a laptop equipped with am Intel core i7 and 16 GB RAM. The considered software to estimate the coregistration shifts is SARPROZ: https://www.sarproz.com/. Finally, MATLAB scripts have been created by the authors to generate the Doppler sub-apertures.

Conflicts of Interest:
The authors declare no conflict of interest. The authors whose names are listed immediately below the title certify that they have no affiliations with or involvement in any organization or entity with any financial interest (such as honoraria; educational grants; participation in speakers' bureaus; membership, employment, consultancies, stock ownership, or other equity interest; and expert testimony or patent-licensing arrangements), or non-financial interest (such as personal or professional relationships, affiliations, knowledge or beliefs) in the subject matter or materials discussed in this manuscript. No founders are involved during the designation and the life of this project.

Abbreviations
The following abbreviations are used in this manuscript: