A Feasibility Study of Nearshore Bathymetry Estimation via Short-Range K-Band MIMO Radar

: This paper provides an assessment of a 24 GHz multiple-input multiple-output radar as a remote sensing tool to retrieve bathymetric maps in coastal areas. The reconstruction procedure considered here exploits the dispersion relation and has been previously employed to elaborate the data acquired via X-band marine radar. The estimation capabilities of the sensor are investigated firstly on synthetic radar data. With this aim, case studies referring to sea waves interacting with a constant and a spatially varying bathymetry are both considered. Finally, the reconstruction procedure is tested by processing real data recorded at Bagnoli


Introduction
Bathymetry measurements play a significant role in effectively managing the marine environment because they allow us to understand the underwater landscape, study marine ecosystems, locate potential hazards, identify areas at risk, and plan marine infrastructures.Additionally, marine sediment dynamics significantly affects the coastline geomorphology and, therefore, continuous monitoring of the seabed allows us to update flood and coastal sea protection schemes [1].
Acoustic sensors, such as multi-beam sonar and single-beam echo sounder, are commonly employed to reconstruct the bathymetry.These sensors are typically mounted under the ship's hull and can effectively measure depths, at fixed locations or along transects, from as shallow as a few meters up to several thousand meters with a high spatial resolution [2].Despite these systems providing a high-resolution reconstruction of the seabed, their use in nearshore waters is limited by the possible hazards for navigation, especially at low tides, and by the high operational costs [3].As a result, remote sensing tools such as satellite altimeter, synthetic aperture radar (SAR), X-band marine radar (MR), light detection and ranging (LiDAR), and optical cameras, are now gaining interest as efficient and cost-effective alternatives for reconstructing the bathymetry in the open sea as well as in shallow waters [4,5].According to [4], remote sensing technologies for bathymetry reconstruction can be classified into non-imaging and imaging methods.
LiDAR and satellite altimeter belong to the first group of methods.Specifically, LiDAR retrieves the water depth from the propagation delay of electromagnetic pulses emitted in the infrared and green wavelength ranges.LiDAR is usually mounted on aerial platforms and measures water depths in the range 1.5-60 m, based on the water transparency [6][7][8].A satellite altimeter measures the distance between the sensor and the sea (or ocean) surface waves and can indirectly estimate the bathymetry.Specifically, variations in the sea surface slope are related to changes in water depth, and this allows us to estimate the bathymetry [9].
On the other hand, imaging methods exploit the shoaling effect, which is observable in optical/radar images.In particular, one of the most common way to retrieve bathymetry from remotely sensed images is the use of dispersion relation, which allows the description of the angular frequency-wavenumber (ω − k) features of sea gravity waves propagating at different water depths [10].
SAR and X-band MR are microwave sensors capable of detecting and characterizing the sea surface depending on the wind speed in the investigated area.Specifically, a minimum speed of 3 m/s is necessary to form the ripples that are responsible for the backscattering of microwave signals [11,12].The bathymetry map can be achieved from SAR or MR data by partitioning the investigated scene into several subareas [13][14][15] to estimate a spatially variable seabed topography.Regarding the SAR-based technique, prior knowledge of sea wavelength and wave period is required.The wavelength is retrieved from the wavenumber spectrum of the considered subarea whilst the wave period is assumed to be known (e.g., via an external sensor).Depending on the sea-state conditions, the SAR-based technique can estimate depths ranging from about 70 m up to 10 m.However, as shown in [16], the presence of wave breaking and wave reflection in coastal areas (i.e., depth < 10 m) affects the accuracy of bathymetry estimation.
With regard to X-band MR, the first contribution to the field of local bathymetry estimation using MR can be attributed to Bell [17].In this work, the local wave celerity was determined by analyzing the coherence between consecutive radar images.On the other hand, one of the most common methods for retrieving the bathymetry map is achieved by performing a 3D Fourier transform of a raw image time sequence for each subarea.Then, the radar spectrum is fitted via the theoretical dispersion relation, with the bathymetry and surface current being the unknown parameters [18,19].A novel approach, introduced in a recent study by Chernyshov [20], adopts a 2D continuous wavelet transform.This method involves the tracking of the spectral peak in the wavenumber domain as waves undergo shoaling and refraction while approaching the shore.
Furthermore, bathymetric field estimates can be achieved by processing optical satellite [27] or airborne [28] data based on the exploitation of sunlight reflection and the observation of hydrodynamic processes, respectively.Each of the above-mentioned approaches has its pros and cons depending on factors such as the water depth, resolution requirements, and survey goal [29].Therefore, the combination of the different sensing techniques is crucial for achieving a more accurate and reliable bathymetry reconstruction of the investigated area [28,30].
Recently, short-range K-band (SRK-band) radars have been proposed to estimate the sea state in very nearshore regions [31][32][33][34].These systems have a short-range coverage, which typically extends up to a few hundred meters.Within this context, this work investigates the capabilities of an SRK-band radar sensor based on frequency modulated continuous wave (FMCW) multiple-input multiple-output (MIMO) technology to reconstruct the bathymetry field in very shallow waters, i.e., depths lower than 4 m.The sensing device considered in this work has recently been proposed to determine the sea surface currents and directional wave spectra [34] and perform marine target detection and tracking [35].The radar signal processing involves two main steps: radar imaging and bathymetry field reconstruction.In the first step, beamforming is performed to achieve a focused radar image of the investigated area from the baseband raw radar signals [36].In the second step, the temporal sequence of radar images is processed to estimate the seafloor topography.With this aim, the dispersion relation describing the sea wave spectrum is exploited [14,22,25,34].
This manuscript deals with a first assessment of the capabilities of the proposed SRKband radar for bathymetry estimation, by processing simulated radar images and real data.It should be stressed that the proposed data processing approach has been already assessed [14,22,25,34] for both X-and SRK-band radar.Herein, we demonstrate for the first time the possibility of reconstructing the seafloor near the shoreline and in partially enclosed areas (e.g., bays or harbors) using an SRK-band FMCW MIMO radar.
With regard to the numerical simulations, sea waves propagating first on a constant depth seabed and later on a spatially varying bathymetry are considered.For the constant sea depth scenario, a synthetic sea wave simulator based on the fast Fourier transform (FFT) [37,38] and the Joint North Sea Wave Observation Project (JONSWAP) model [39] for the sea wave spectrum is implemented.As for the variable depth scenario, a hydrodynamic solver for coastal dynamics is used to describe wave propagation over the reference bathymetry [40][41][42].Once the numerical wave field has been generated, synthetic radar images are produced by means of the 3D FMCW MIMO radar data simulator presented in [34].Finally, an experimental testing of the system performance is provided by processing data measured at the seafront in Bagnoli, Naples, South Italy.
The manuscript is organized as follows.Section 2 presents the numerical sea-wave models.The radar system and the synthetic data generation and processing are detailed in Section 3. The numerical and experimental validation is presented in Sections 4 and 5, respectively.Discussions and conclusions are reported in Section 6.

Numerical Sea-Wave Model
The models adopted for generating the synthetic sea wave field (SWF) are described in this section.Specifically, a SWF over a constant bathymetry is first produced by adopting the Fourier domain approach [38].Then, random waves are generated and propagated over a planar-sloping beach according to the model for coastal dynamics described in [40].For both the wave models, the JONSWAP sea wave spectrum is assumed [39].

Wave Spectrum Model
The JONSWAP spectrum is often used in engineering and coastal studies to estimate the wave conditions and design criteria for offshore structures, coastal protection, and coastal engineering projects.Such a spectrum provides a convenient way to estimate the expected wave characteristics based on wave frequency and direction data.The JONSWAP spectrum is based on the assumption that sea waves can be described as the summation of many wave components with different frequencies and amplitudes.
The 1D JONSWAP wave model is defined as [39]: where α = 0.0076 is the energy scale factor, ω p = 22 is the spectral peak angular frequency, g = 9.81 m/s 2 is the gravitational acceleration, and σ is the peak shape parameter.Additionally, W 10 and L f are the wind speed at a height of 10 m above the mean sea level (MSL) and the fetch length, respectively.A key parameter in Equation ( 1) is the peak-enhancement factor γ, which ranges from 1 to 7, with γ = 3.3 representing the standard sea wave spectrum.

Fourier Domain Approach
Synthetic sea images for a constant bathymetry scenario are generated in a very straightforward and computationally efficient manner by using the FFT algorithm.According to [38], the complex Fourier amplitudes of a wave elevation field at time zero are given by: where ξ Re and ξ Im are independent Gaussian random variables with zero mean and unity standard deviation [38], and k = k x , k y is the wave vector with amplitude The function P JS (k) accounts for the spectral properties of the SWF under the 2D JONSWAP model [39].
According to the dispersion relation [10], the 2D wavenumber directional spectrum is obtained from Equation (1) and expressed as: where k p is the peak wave number, and S JS (θ) is the directional spreading function defined as [40]: where θ is the propagation direction, and θ p is the direction of the dominant wave (i.e., wind direction).
According to the linear wave theory, sea surface waves undergo a dispersion phenomenon that is governed by the dispersion relation [10].
where d and U = U x , U y are the bathymetry and the sea surface current vector, respectively.Given Equation ( 8), the time-dependent amplitude of the sea-wave height at time t is expressed as: Equation ( 9) conserves the complex conjugate property and, based on it, the wave elevation is evaluated as a real quantity thanks to an inverse FFT [38], i.e., (10) in which r = (x, y) is a point in the plane z = 0 that is assumed to be coincident with the MSL.

Wave Resolving Model
Wave propagation over non-planar bathymetries is obtained by using the depth semiaveraged scheme described in [41,42].This model consists of a subset of depth averaged quantities (i.e., the total water depth h = η + d, the depth-averaged velocity V, and the generalized mass flux M) and of a Poisson equation for the variable Y that is the integral of the vertical component of the local velocity field from a certain quote over the depth.Succinctly, the subset of depth-averaged quantities is written as: where ∇ is the nabla operator in two dimensions, F is the flux tensor, and S is the source term (more details can be found in [41]).Finally, the Poisson equation for Y is: where ∆ is the Laplace operator in three dimensions.The subset of depth-averaged equations is represented via the finite volume method, and the fluxes across volumes are modeled through the MUSCL-Hancock scheme described in [43].Conversely, the Poisson equation is discretized by adopting second-order finite differences.The overall system is integrated in time via a fourth-order Adams-Bashforth-Moulton predictor/corrector scheme.We direct the reader to [41] for more details about the model and its numerical implementation.

Sensor and Data Processing
This section describes the SRK-band radar platform used for bathymetry estimation, as well as the procedure to generate synthetic raw radar images from sea images, by taking into account the tilt modulation and shadowing phenomenon.Moreover, the data-processing pipeline for bathymetry retrieval is also described.

Radar Platform and Signal Model
SRK-band radar is a compact, lightweight, and low power 24 GHz FMCW MIMO system developed by Inras Gmbh for automotive applications [35].The main features of the radar system, which operates within the frequency range of 24-24.25 GHz, are widely detailed in [34,35], and herein are briefly summarized.Specifically, Figure 1 displays the arrangement of antennas, which consists of two transmitting (Tx) and eight receiving (Rx).These antennas comprise linearly polarized patch arrays, each containing eight radiating elements, exhibiting a narrow vertical beam width (±6.4 • ) and a wide horizontal beam width (±75 • ).The 2 × 8 MIMO array is equivalent to a monostatic virtual array composed of 16 antennas (channels), where an overlap is present between the eighth and ninth elements {n = 0, . . .15, n ̸ = 8} [35].
Remote Sens. 2024, 16, x FOR PEER REVIEW 6 of 17 The readers can find more detailed information in [34,35] to delve deeper into the radar platform and signal model.

Radar Data Simulation (Forward Model)
Synthetic radar data are obtained by exploiting a 3D radar signal simulator applied to the SWF.Originally developed in [34], this simulator accounts for the main physical phenomena involved in the radar imaging process, such as shadowing and tilt modulation, and it is briefly described here for the reader's convenience.
The sea surface scattering is described using a facet model and the radar signal is evaluated as a summation of the contributions related to each facet  inside the radar FoV.In addition, the presence on the sea surface of long waves, which modulate the radar signal, and capillary waves, which are responsible for the main backscattering contribution (resonant Bragg scattering), is taken into account by exploiting a two-scale model [44][45][46].It is worth noticing that in comparison to the X-band radar, the K-band Bragg waves undergo significant modulation at low wind conditions [47] For the sake of simplicity, focusing on a single Tx-Rx antenna pair related to a generic radar channel n, Figure 2 shows the tilt modulation and shadowing phenomenon that occur in the vertical plane.The tilt modulation accounts for the slope of the observed face, which depends on the angle  between the incidence direction ̂ and the normal  to the facet p defined by The shadowing phenomenon occurs when the radar antenna does not receive any signal from the shadowed facet of the sea surface, i.e., when the following condition holds: As for the resolution, the SRK-band radar has a fine range resolution (∆r = 0.6 m) and an angular resolution that varies with the target's direction from the boresight.In other words, the best angular resolution (c.a.7.6 • ) is achieved at the boresight, whereas, as the target moves away from the boresight, the angular resolution progressively worsens, becoming approximately 30 • at the edge of the azimuth field of view (FoV).
Equation ( 14) reports the transmitted signal that is a chirp with duration T c and bandwidth B w = f max − f min , i.e., where A t is the amplitude, t is the fast time, and ε = B w /T c is the chirp rate.The received signal over channel n corresponds to a time-delayed version of (14) and is demodulated at baseband.By performing a filtering of the double frequency terms, the intermediate frequency (IF) signal can be written as: where τ n is the time delay over channel n, and ∼ A r is an amplitude factor.The raw data are represented by the baseband signals y n t with n = 0, . . .15, n ̸ = 8 for every frame m = 1, . . ., M.
The readers can find more detailed information in [34,35] to delve deeper into the radar platform and signal model.

Radar Data Simulation (Forward Model)
Synthetic radar data are obtained by exploiting a 3D radar signal simulator applied to the SWF.Originally developed in [34], this simulator accounts for the main physical phenomena involved in the radar imaging process, such as shadowing and tilt modulation, and it is briefly described here for the reader's convenience.
The sea surface scattering is described using a facet model and the radar signal is evaluated as a summation of the contributions related to each facet p inside the radar FoV.In addition, the presence on the sea surface of long waves, which modulate the radar signal, and capillary waves, which are responsible for the main backscattering contribution (resonant Bragg scattering), is taken into account by exploiting a two-scale model [44][45][46].It is worth noticing that in comparison to the X-band radar, the K-band Bragg waves undergo significant modulation at low wind conditions [47].
For the sake of simplicity, focusing on a single Tx-Rx antenna pair related to a generic radar channel n, Figure 2 shows the tilt modulation and shadowing phenomenon that occur in the vertical plane.The tilt modulation accounts for the slope of the observed face, which depends on the angle β p between the incidence direction ŝp and the normal np to the facet p defined by Remote Sens. 2024, 16, x FOR PEER REVIEW 7 of 17

Data-Processing Approach (Inverse Model)
The strategy adopted to process the data frame sequence has been already presented and discussed in detail in [34] and involves two phases: radar imaging and bathymetry retrieval.
Specifically, beamforming based on a 2D FFT (see [35] for more details) is performed in the first step to obtain a time sequence of focused intensity images at polar coordinates, i.e., (, Θ, ).The intensity images are subsequently linearly interpolated over a Cartesian grid, thus obtaining the image sequence (, , ), which is given as the input to the second block of the data processing chain in Figure 3.To this point, the bathymetry field is reconstructed by partitioning each radar image into  partially overlapping subareas.It is important to note that in each subarea, the spatial homogeneity and temporary stationarity of the sea-state parameters and bathymetric conditions are assumed [21][22][23][24][25].The shadowing phenomenon occurs when the radar antenna does not receive any signal from the shadowed facet of the sea surface, i.e., when the following condition holds: where ρ p is the projection of the slant range r p on the plane at the MSL, and ψ p is the grazing angle defined by the propagation direction ŝp and the horizontal direction ρ.The radar signal scattered by the sea surface is given by the superposition of the echoes from every illuminated facet p, i.e., In Equation ( 18), Γ is the set of illuminated facets, µ p = 1 r 2 p is the wave attenuation in free-space, A p is the area of the p-th facet, σ 0p and ζ p are the normalized radar cross-section (NRCS) and phase shift due to the backscattering of the facet.To this point, a raw data frame y n t is achieved at any slow time t = t 1 , t 2 , . . ., t M .

Data-Processing Approach (Inverse Model)
The strategy adopted to process the data frame sequence has been already presented and discussed in detail in [34] and involves two phases: radar imaging and bathymetry retrieval.
Specifically, beamforming based on a 2D FFT (see [35] for more details) is performed in the first step to obtain a time sequence of focused intensity images at polar coordinates, i.e., I(ρ, Θ, t).The intensity images are subsequently linearly interpolated over a Cartesian grid, thus obtaining the image sequence I(x, y, t), which is given as the input to the second block of the data processing chain in Figure 3.To this point, the bathymetry field is reconstructed by partitioning each radar image into N s partially overlapping subareas.It is important to note that in each subarea, the spatial homogeneity and temporary stationarity of the sea-state parameters and bathymetric conditions are assumed [21][22][23][24][25].

Data-Processing Approach (Inverse Model)
The strategy adopted to process the data frame sequence has been already presented and discussed in detail in [34] and involves two phases: radar imaging and bathymetry retrieval.
Specifically, beamforming based on a 2D FFT (see [35] for more details) is performed in the first step to obtain a time sequence of focused intensity images at polar coordinates, i.e., (, Θ, ).The intensity images are subsequently linearly interpolated over a Cartesian grid, thus obtaining the image sequence (, , ), which is given as the input to the second block of the data processing chain in Figure 3.To this point, the bathymetry field is reconstructed by partitioning each radar image into  partially overlapping subareas.It is important to note that in each subarea, the spatial homogeneity and temporary stationarity of the sea-state parameters and bathymetric conditions are assumed [21][22][23][24][25].The bathymetry is retrieved by exploiting the normalized scalar product (NSP) technique, which maximizes the scalar product between the 3D image spectrum in each sub-area  8)).It is important to highlight that different values of sea surface current and bathymetry can lead to changes in the dispersion relation as well as the radar image spectrum [10].
Based on this effect, the NSP procedure is able to jointly retrieve both parameters [22].However, since the purpose of this study is bathymetry estimation, the sea surface current components are assumed to be negligible, as in a number of previous works [14,18,20,25,26], and we assume this a priori information in the reconstruction.Accordingly, the local bathymetry value ďj is found as where G(•) is the characteristic function based on Equation ( 8), ⟨ •, • ⟩ is the scalar product of the functions F (3) j and G, and P F j is the spectrum power.

Numerical Validation
In this section, the data processing approach described in Section 3 is assessed by considering SWFs that propagate on constant and variable bathymetry.

SWFs with Constant Bathymetry
Three SWFs (i.e., SWF 1 , SWF 2 , and SWF 3 ) are generated by using custom codes developed in the MATLAB 2019 programming environment.These wave models are characterized by bathymetry values equal to 1.5 m, 1 m, and 0.8 m, respectively.
The 2D JONSWAP spectrum associated with each SWF is specified via various spectral parameters.Specifically, ω p , k p , and L f are assigned different values, as seen in Table 1, while θ p = 0 • , W 10 = 15 m/s, and γ = 3.3.According to the above parameters, the significant wave height H s of the generated SWFs varies within the range of 0.71-1.21m.According to the Douglas scale [48], this range corresponds to a rough sea condition.A temporal sequence of M = 256 raw radar images is generated with a uniform spacing ∆x = ∆y = 0.5 m and time interval ∆t = 0.5 s between two consecutive images, leading to a total observation time equal to T w = 128 s.With regard to the radar data simulation, the electrical and geometrical radar parameters are listed in Table 2.A range coverage from ρ min = 34 m to ρ max = 149 m is achieved by considering the antenna located at a quota H = 10 m above the MSL and tilted about 10 • in the vertical plane.It is assumed that the radar boresight points in the arrival direction θ p of the sea wave.As demonstrated in [34] regarding the retrieval of surface currents, this condition yields a large region of FoV in which reliable reconstructions can be achieved.Of course, the reconstruction procedure works also when the radar boresight does not point in the direction of the incoming SWF, but a smaller angular region where reliable reconstructions are achieved (see [34]).The radar images are produced via spatial domain beamforming on a polar grid, taking into account the specified range interval and azimuth FoV (±75 • ).Subsequently, a linear interpolation on a Cartesian grid is performed, with spatial discretization steps equal to ∆x = ∆y = 0.5 m.A snapshot of the SWF 1 elevation and its corresponding radar image at t = 64 s are shown in Figures 4a and 4b, respectively.To this point, the estimation procedure described in Section 3.

MRPE = ∑
× 100 An additional figure of merit considered here is the signal-to-noise ratio (SNR), which measures the strength of the desired signal compared to the level of noise in the 3D radar spectrum, i.e.,  = 10 log In the above expression,  is the power of the radar spectrum in the dispersion shell, and  is the noise power evaluated over the complementary spectral region.More specifically, it has been found out that a low SNR is correlated with the spatial regions in the radar FoV where the estimates are less reliable.Therefore, a threshold value for the minimum acceptable SNR is introduced to define a reliability region for the retrieved bathymetry map.
The latter concept is clearly understood by analyzing the results shown in Figure 5. Here, Figure 5a,b display the SNR map and the retrieved bathymetry map referred to SWF1, respectively.As can be observed in Figure 5a, the SNR values decrease significantly far from the boresight.The decay of the radar spectrum energy away from the arrival direction of the sea wave and the reduced angular resolution makes the data less reliable, as previously pointed out in [34].Accordingly, we assume that the spatial region where the estimates can be considered reliable is the one where SNR ≥ 3 dB.From this point on, To this point, the estimation procedure described in Section 3.3 and displayed in Figure 3 is applied to the whole temporal image sequence.Specifically, the local estimation of bathymetry involves partitioning radar images into N s subareas with a size of 40 m × 40 m, with a 10 m overlap in both the x and y directions.The accuracy of the bathymetric reconstruction is assessed by evaluating the mean relative percentage error (MRPE).

MRPE =
1 An additional figure of merit considered here is the signal-to-noise ratio (SNR), which measures the strength of the desired signal compared to the level of noise in the 3D radar spectrum, i.e., SNR = 10log 10 In the above expression, P s is the power of the radar spectrum in the dispersion shell, and P n is the noise power evaluated over the complementary spectral region.More specifically, it has been found out that a low SNR is correlated with the spatial regions in the radar FoV where the estimates are less reliable.Therefore, a threshold value for the minimum acceptable SNR is introduced to define a reliability region for the retrieved bathymetry map.
The latter concept is clearly understood by analyzing the results shown in Figure 5. Here, Figure 5a,b display the SNR map and the retrieved bathymetry map referred to SWF 1 , respectively.As can be observed in Figure 5a, the SNR values decrease significantly far from the boresight.The decay of the radar spectrum energy away from the arrival direction of the sea wave and the reduced angular resolution makes the data less reliable, as previously pointed out in [34].Accordingly, we assume that the spatial region where the estimates can be considered reliable is the one where SNR ≥ 3 dB.From this point on, it is implicitly assumed that the bathymetry results are obtained by adopting the SNR ≥ 3 dB criterion (see Figure 5b).The estimated bathymetry maps relevant to the SWF 2 and SWF 3 are illustrated in Figures 5c and 5d, respectively.
Table 3 provides a summary of the MRPE values obtained for the considered SWFs, revealing that successful reconstructions are achieved in all cases.

Wave Propagation over Planar Sloping Beach
As shown in Figure 6a, we consider a planar sloping beach with slope tan(ϑ) = −0.0175 in the y direction, where ϑ is the angle between the seabed and the y-axis.The numerical domain ranges between y = 0 and y = 200 m, and the still water depth increases from 0.5 m up to 4 m, respectively.Random waves are generated along the boundary at y = 200 m (inflow boundary) using a JONSWAP spectrum.Then, waves propagate over the bathymetry and leave the domain across an open boundary at y = 0 m.Three different SWFs (i.e., SWF PSB1 , SWF PSB2 , and SWF PSB3 in Table 4), characterized by their peak frequency and significant wave height, propagate along a sloping beach.The numerical domain is discretized through a Cartesian grid with spatial steps ∆x = ∆y = 0.1 m in the horizontal plane.Conversely, a stretched grid is adopted over the water depth going from ∆z = 0.1 m at z = 0 m up to ∆z = 0.4 m at z = −3 m.The simulation starts at t = 0 s, with a time discretization step equal to 0.005 s, and the wave evolution is saved any ∆t = 0.5 s for t > 60 s.The latter time interval is needed to allow waves to propagate all over the numerical domain.After that, 256 snapshots of the wave evolution are recorded.
Remote Sens. 2024, 16, x FOR PEER REVIEW 10 of 17 it is implicitly assumed that the bathymetry results are obtained by adopting the SNR ≥ 3 dB criterion (see Figure 5b).The estimated bathymetry maps relevant to the SWF2 and SWF3 are illustrated in Figure 5c and Figure 5d, respectively.Table 3 provides a summary of the MRPE values obtained for the considered SWFs, revealing that successful reconstructions are achieved in all cases.

Wave Propagation over Planar Sloping Beach
As shown in Figure 6a, we consider a planar sloping beach with slope tan() = −0.0175 in the y direction, where  is the angle between the seabed and the y-axis.The  4), characterized by their peak frequency and significant wave height, propagate along a sloping beach.The numerical domain is discretized through a Cartesian grid with spatial steps ∆ = ∆ = 0.1 m in the horizontal plane.Conversely, a stretched grid is adopted over the water depth going from ∆ = 0.1 m at z = 0 m up to ∆ = 0.4 m at z = −3 m.The simulation starts at t = 0 s, with a time discretization step equal to 0.005 s, and the wave evolution is saved any ∆ = 0.5 s for  > 60 s.The latter time interval is needed to allow waves to propagate all over the numerical domain.After that, 256 snapshots of the wave evolution are recorded.Figure 6a shows a snapshot of the generated sea-wave elevation at  = 64 s.As can be observed, long-crested waves propagate over a seabed with linear bathymetry ranging from 4 m at the inflow up to 0.5 m at the outflow.The radar image corresponding to the sea-wave elevation snapshot is depicted in Figure 6b.The radar settings and data-processing parameters to obtain the radar image from the synthetic wave field are those previously reported in Table 2.As for the radar data processing, the bathymetry is retrieved by processing a time sequence of 256 radar images with a sub-area size equal to 40 × 40 m and a 10 m overlap along both the  and  direction.
Figure 7a-c show the estimated bathymetry for each considered case based on the SNR ≥ 3 dB criterion.As can be observed, the spatial variability of the bathymetry is quite well reconstructed, and the good agreement is also confirmed by the scatterplots of the depth values reported in Figure 8a-c.Furthermore, a statistical analysis is conducted, and the results in terms of the MPRE and correlation coefficient R are presented in Table 5.   Figure 6a shows a snapshot of the generated sea-wave elevation at t = 64 s.As can be observed, long-crested waves propagate over a seabed with linear bathymetry ranging from 4 m at the inflow up to 0.5 m at the outflow.The radar image corresponding to the sea-wave elevation snapshot is depicted in Figure 6b.The radar settings and dataprocessing parameters to obtain the radar image from the synthetic wave field are those previously reported in Table 2.As for the radar data processing, the bathymetry is retrieved by processing a time sequence of 256 radar images with a sub-area size equal to 40 × 40 m and a 10 m overlap along both the x and y direction.
Figure 7a-c show the estimated bathymetry for each considered case based on the SNR ≥ 3 dB criterion.As can be observed, the spatial variability of the bathymetry is quite well reconstructed, and the good agreement is also confirmed by the scatterplots of the depth values reported in Figure 8a-c.Furthermore, a statistical analysis is conducted, and the results in terms of the MPRE and correlation coefficient R are presented in Table 5.
from 4 m at the inflow up to 0.5 m at the outflow.The radar image corresponding to the sea-wave elevation snapshot is depicted in Figure 6b.The radar settings and data-processing parameters to obtain the radar image from the synthetic wave field are those previously reported in Table 2.As for the radar data processing, the bathymetry is retrieved by processing a time sequence of 256 radar images with a sub-area size equal to 40 × 40 m and a 10 m overlap along both the  and  direction.
Figure 7a-c show the estimated bathymetry for each considered case based on the SNR ≥ 3 dB criterion.As can be observed, the spatial variability of the bathymetry is quite well reconstructed, and the good agreement is also confirmed by the scatterplots of the depth values reported in Figure 8a-c.Furthermore, a statistical analysis is conducted, and the results in terms of the MPRE and correlation coefficient R are presented in Table 5.

Experimental Testing
This section describes a first test of the radar system for bathymetry reconstruction in the marine environment.An experiment was carried out at Bagnoli Bay, South Italy, on 26 July 2023 at 11:30 a.m.local time.Figure 9 depicts a satellite image of the surveyed area taken from Google Earth Pro, with the radar location highlighted in an inset overlaid on a radar image.The tags A1 and A2 in Figure 9 show the location of two anemometric stations within the bay and close to the measurement location, from which information about the wind speed and direction was retrieved [49].To this point, Figure 11a shows the bathymetry field of the entire investigated area and, since there are no external reference data available for the bathymetry as well as the wave amplitude in this area, it is not possible to verify the accuracy of the estimation.However, as suggested via numerical simulations, the accuracy of the estimations is expected to grow in the direction of the incoming waves, which differs by a few degrees from the radar boresight, as seen from the data.The SNR ≥ 3 dB criterion is adopted to delimit the spatial region and confirms that the estimations are more reliable in the neighborhood of the arrival direction of the sea waves (see the black curve in Figure 11a).Finally, the dependable bathymetric estimations have been georeferenced and are shown in Figure 11b, confirming the smoother behavior of the bathymetry field with depth values falling within the approximate range of 1-2.5 m.Moreover, for a qualitative analysis of Bagnoli Bay is the subject of many studies by researchers, since it belongs to the Campi Flegrei volcanic area and it is an important underwater archeological park [50,51].Considering this, a digital terrain model (DTM) of the seabed was retrieved to monitor the "Bradyseism" phenomenon [50] and retrieve the main archaeological features of the Bay [51].However, the DTM shows a minimum depth of about 5 m reached at a distance of about 300 m from the coastline.This limitation arises because the acoustic sensors installed under the ship's hull can be endangered as the ship approaches shallow waters.
The geographical coordinates of the radar were determined as 40 • 48 ′ 58.99 ′′ N, 14 • 9 ′ 37.29 ′′ E by using the GNSS receiver of a smartphone.Additionally, the direction of the radar antenna boresight was established as 215 • SW thanks to the magnetic compass of a smartphone.To ensure protection against water damage, the radar unit was housed within a waterproof case and securely positioned on a tripod.As shown in Figure 10a, the system was controlled by a PC and located on a rock at a height of about 10 m above the MSL.The radar data were gathered on a cloudy day and, based on the anemometric information, with a fresh breeze (Beaufort scale [52]) from the W direction and a speed equal to 10 m/s that induced a moderate sea (Douglas scale [48]).Consequently, considering this information along with the image depicted in Figure 10a, it can be inferred that the sea surface during radar acquisition exhibited rough conditions.The radar setup configuration is outlined in Table 2, save for the values of ρ min and ρ max , both set at 20 m and 210 m, respectively.Additionally, the Cartesian grid interpolation of the radar images is performed with a pixel spacing ∆x = ∆y = 0.5 m.To this point, Figure 11a shows the bathymetry field of the entire investigated area and, since there are no external reference data available for the bathymetry as well as the wave amplitude in this area, it is not possible to verify the accuracy of the estimation.However, as suggested via numerical simulations, the accuracy of the estimations is expected to grow in the direction of the incoming waves, which differs by a few degrees from the radar boresight, as seen from the data.The SNR ≥ 3 dB criterion is adopted to delimit the spatial region and confirms that the estimations are more reliable in the neighborhood of the arrival direction of the sea waves (see the black curve in Figure 11a).Finally, the dependable bathymetric estimations have been georeferenced and are shown in Figure 11b, confirming the smoother behavior of the bathymetry field with depth values falling within the approximate range of 1-2.5 m.Moreover, for a qualitative analysis of the reconstruction performance, the reference bathymetry line at depth 2 m is overlaid on The radar image sequence is partitioned into subareas of the size 40 m × 40 m and with an overlap of 10 m along both the x and y directions.Figure 10b illustrates a cross-section of the 3D radar spectrum within a specific subarea located at coordinates x = 0 m and y = 98 m, where the NSP strategy is employed to estimate the bathymetry value.
To this point, Figure 11a shows the bathymetry field of the entire investigated area and, since there are no external reference data available for the bathymetry as well as the wave amplitude in this area, it is not possible to verify the accuracy of the estimation.However, as suggested via numerical simulations, the accuracy of the estimations is expected to grow in the direction of the incoming waves, which differs by a few degrees from the radar boresight, as seen from the data.The SNR ≥ 3 dB criterion is adopted to delimit the spatial region and confirms that the estimations are more reliable in the neighborhood of the arrival direction of the sea waves (see the black curve in Figure 11a).Finally, the dependable bathymetric estimations have been georeferenced and are shown in Figure 11b, confirming the smoother behavior of the bathymetry field with depth values falling within the approximate range of 1-2.5 m.Moreover, for a qualitative analysis of the reconstruction performance, the reference bathymetry line at depth 2 m is overlaid on Figure 11b.This isobath is extracted from the nautical chart available on the website [53], which is derived from the official publications of the Hydrographic Institute-Italian Navy.It is observed that the estimated bathymetry values are in good agreement with the contour line.
It is important to emphasize that, even in this scenario, the sea surface current components are neglected in the dispersion relation, and this may lead to bias in the bathymetry estimate [54].
which is derived from the official publications of the Hydrographic Institute-Italian Navy.It is observed that the estimated bathymetry values are in good agreement with the contour line.
It is important to emphasize that, even in this scenario, the sea surface current components are neglected in the dispersion relation, and this may lead to bias in the bathymetry estimate [54].

Discussion and Conclusions
An SRK-band radar based on an FMCW MIMO system has been proposed for the reconstruction of bathymetry in nearshore areas.The system capabilities have been assessed first on simulated data under different sea state conditions in terms of significant wave height ranging from 0.25 m to 1.21 m.To this end, constant and variable seabed scenarios have been simulated via the Fourier domain approach and an ad hoc developed

Discussion and Conclusions
An SRK-band radar based on an FMCW MIMO system has been proposed for the reconstruction of bathymetry in nearshore areas.The system capabilities have been assessed first on simulated data under different sea state conditions in terms of significant wave height ranging from 0.25 m to 1.21 m.To this end, constant and variable seabed scenarios have been simulated via the Fourier domain approach and an ad hoc developed waveresolving model, respectively.After, radar images have been generated with the aid of a 3D radar simulator, taking into account the tilt modulation and shadowing phenomenon.The numerical tests have confirmed that the radar provides reliable estimates with errors lower than 8.5% in an angular interval where SNR ≥ 3 dB, and a correlation coefficient greater than 0.90 for the wave fields propagating over a planar slope bathymetry.A preliminary field trial carried out at the seafront of Bagnoli Bay has confirmed that the radar prototype and a related signal processing pipeline can reconstruct the bathymetry field in the nearshore area, providing results in agreement with the limited available data.
However, it should be stressed that the herein presented results were achieved with the use of a linear dispersion relation in shallow waters, which may introduce some estimation errors.In fact, the linear dispersion relation in shallower area may become not very accurate due to the presence of wave breaking leading to an error exceeding 30% for the estimated bathymetry [55,56].
Furthermore, the sea state homogeneity assumption within each subarea contradicts the inherent characteristics of the shoaling wave field, particularly in the presence of abrupt bathymetric variations [14,20].
Therefore, future activities will be focused on overcoming the above limitations.
As another future activity, we will consider a multi-sensor approach based on the joint exploitation of sensors such as wave buoys, marine drone-mounted echo sounders, an acoustic Doppler current profile, and an anemometric station.This comprehensive setup will be useful to calibrate and validate bathymetry measurements conducted via SRK-band radar and determine the SRK-band modulation transfer function necessary for the subsequent estimation of the significant wave height, period, and wavelength.

Figure 1 .
Figure 1.RF front end of the radar system.

Figure 1 .
Figure 1.RF front end of the radar system.

Figure 2 .
Figure 2. Representation of the shadowing and tilt modulation effect in the vertical plane.The Tx/Rx antenna pair is located at quota H above the MSL. and  delimit the radar FoV.

Figure 2 .
Figure 2. Representation of the shadowing and tilt modulation effect in the vertical plane.The Tx/Rx antenna pair is located at quota H above the MSL.ρ min and ρ max delimit the radar FoV.

Figure 2 .
Figure 2. Representation of the shadowing and tilt modulation effect in the vertical plane.The Tx/Rx antenna pair is located at quota H above the MSL. and  delimit the radar FoV.

Figure 3 .Figure 3 .
Figure 3. Data processing chain for bathymetry map estimation.The bathymetry is retrieved by exploiting the normalized scalar product (NSP) technique, which maximizes the scalar product between the 3D image spectrum in each subarea  ( ) (, ) j=1,...,N s and a characteristic function G(•) based on the theoretical dispersion relation (see Equation (
3 and displayed in Figure 3 is applied to the whole temporal image sequence.Specifically, the local estimation of bathymetry involves partitioning radar images into  subareas with a size of 40 m × 40 m, with a 10 m overlap in both the x and y directions.The accuracy of the bathymetric reconstruction is assessed by evaluating the mean relative percentage error (MRPE).

Figure 5 .
Figure 5. (a) SNR map for SWF1 with, superimposed upon it, the black solid line enclosing the area where SNR ≥ 3dB.(b-d) Reconstructed bathymetry field for SWF1, SWF2, and SWF3, respectively.
numerical domain ranges between y = 0 and y = 200 m, and the still water depth increases from 0.5 m up to 4 m, respectively.Random waves are generated along the boundary at y = 200 m (inflow boundary) using a JONSWAP spectrum.Then, waves propagate over the bathymetry and leave the domain across an open boundary at y = 0 m.Three different SWFs (i.e., SWFPSB1, SWFPSB2, and SWFPSB3 in Table

Figure 5 .
Figure 5. (a) SNR map for SWF 1 with, superimposed upon it, the black solid line enclosing the area where SNR ≥ 3dB.(b-d) Reconstructed bathymetry field for SWF 1 , SWF 2 , and SWF 3 , respectively.

Figure 6 .
Figure 6.Synthetic SWFPSB1 at time  = 64 s.(a) Sea waves propagating over the linear bathymetry.(b) Radar image.Figure 6. Synthetic SWF PSB1 at time t = 64 s.(a) Sea waves propagating over the linear bathymetry.(b) Radar image.

Figure 6 .
Figure 6.Synthetic SWFPSB1 at time  = 64 s.(a) Sea waves propagating over the linear bathymetry.(b) Radar image.Figure 6. Synthetic SWF PSB1 at time t = 64 s.(a) Sea waves propagating over the linear bathymetry.(b) Radar image.

Figure 8 .
Figure 8.Comparison between estimated bathymetry (dots) and ground truth (red dashed line).The equation of the regression line (black dashed line) is provided in the legend.(a-c) Scatterplot for SWF PSB1 , SWF PSB2 , and SWF PSB3 , respectively.
Remote Sens. 2024,16,  x FOR PEER REVIEW 13 of 17The radar image sequence is partitioned into subareas of the size 40 m × 40 m and with an overlap of 10 m along both the x and y directions.Figure10billustrates a crosssection of the 3D radar spectrum within a specific subarea located at coordinates  = 0 m and  = 98 m, where the NSP strategy is employed to estimate the bathymetry value.

Figure 9 .
Figure 9. Satellite image of the test site.The inset displays the radar's positioning overlaid on a radar image.A1 and A2 show the location of two anemometric stations.

Figure 9 .
Figure 9. Satellite image of the test site.The inset displays the radar's positioning overlaid on a radar image.A1 and A2 show the location of two anemometric stations.

Figure 9 .Figure 10 .
Figure 9. Satellite image of the test site.The inset displays the radar's positioning overlaid on a radar image.A1 and A2 show the location of two anemometric stations.

Figure 11 .
Figure 11.(a) Bathymetry field estimated at Bagnoli Bay.The black solid line encloses the area where the SNR ≥ 3 dB.(b) Georeferenced SNR bathymetry field based on SNR ≥ 3dB criterion.The white line is the isobath at depth 2 m.

Figure 11 .
Figure 11.(a) Bathymetry field estimated at Bagnoli Bay.The black solid line encloses the area where the SNR ≥ 3 dB.(b) Georeferenced SNR bathymetry field based on SNR ≥ 3dB criterion.The white line is the isobath at depth 2 m.

Table 1 .
Main characteristic of SWFs.

Table 2 .
Parameters defining the radar configuration.

Table 2 .
Parameters defining the radar configuration.
x Max.range (above MSL) 149 m

Table 3 .
MRPE values for estimated bathymetry with the SNR ≥ 3 dB criterion.

Table 3 .
MRPE values for estimated bathymetry with the SNR ≥ 3 dB criterion.

Table 4 .
Main characteristic of SWF PSB .

Table 5 .
MRPE and R values for estimated bathymetry using the SNR ≥ 3 dB criterion.

Table 5 .
MRPE and R values for estimated bathymetry using the SNR ≥ 3 dB criterion.