An improved algorithm for phase-resolved sea surface reconstruction from X-band marine radar images

We present a modified methodology for phase-resolved surface wave reconstruction from incoherent X-band marine radar images. The method is based on the linear wave theory and uses the linear dispersion relation to extract the valuable signals associated with gravity waves. A parameter optimization of the proposed modification is performed based on simulated synthetic radar images. The quantitative comparisons in the accuracy of the standard and modified reconstruction methods are made for both simulated and real radar images. The correlation coefficient between reconstructed and true wave elevations is improved up to 0.9–0.92 for the present modified method from 0.69 to 0.74 for the standard method for the simulated sea surfaces. The wave spectra reconstructed from the real X-band radar measurements are in good agreement with those obtained from the independent point measurement by Miros RangeFinder for both unimodal and bimodal seas.


Introduction
Ocean wave prediction is of critical importance in marine operations and maritime transport. It is essentially required in the development and application of advanced technologies that help significantly reduce the costs, increase the safety and thus enhance the operation envelop of a wide range of challenging operations under severe sea conditions. Statistical methods have played a vital role in modeling and predicting the sea state conditions. They provide the forecast of statistical parameters of ocean environments useful in relatively long-term planning of marine activities and maritime transport. Without phase information, however, the statistical modeling cannot provide a quantitative prediction of the specific times and/or spatial locations on the occurrence of individual extreme events, such as rogue waves and ship capsizing. In practice, the marine operations are indeed limited by these extreme events that may cause structural damages and/or marine safety issues. In addition to the statistical modeling, phase-resolved methods have been developed to predict the spatial-temporal evolution of detailed sea surface elevation maps. The application of phase-resolved modeling enables one to obtain a deterministic forecast of individual large-amplitude waves and corresponding vessel/platform responses. The forecasted phase-resolved wave and structure motion information can be used to improve maneuverability and control of vessel motion and operations, provide operational guidance for improved safety and performance, increase operational envelop and survivability in rough seas (e.g. lifting operations, helicopter landing, and landing of subsea modules into platform installation), and optimize ship deployment and multi-vessel cooperation (e.g. transporting cargoes between ships or mobile platforms). In practical applications, a major challenge in phase-resolved wave prediction is to obtain reliable and accurate phase-resolved wave field measurements that can be used to initialize the phase-resolved wave model.
There is a great deal of research and development on the use of X-band marine radars for ocean wave measurements. Incoherent X-band marine radars are widely used on vessels to monitor the statistical parameters of the ocean waves and for navigation aids. Incoherent X-band marine radars can scan a sea surface area within a radius of 2-5 km continuously in space and time with high resolution and produce backscatter images. However, the intensity images do not 1 3 represent the real sea surface wave shape since they are complexly influenced by a number of physical factors including hydrodynamic modulation, tilt modulation, shadowing, scanning range and direction, and wind and wave conditions. It has been shown on the basis of the linear gravity wave theory that the use of a three-dimensional (i.e. two dimensions in space and one dimension in time) Fourier analysis of the radar intensity images could yield valuable information on both ocean wave directionality and near surface currents (Young et al. 1985;Dankert and Rosenthal 2004;Nieto Borge et al. 1999;Nieto-Borge and Guedes-Soares 2000;Hessner et al. 2001). Nieto-Borge et al. (2004) proposed an extension of the existing inverse modeling technique to derive wave spectra from incoherent marine radar images. This algorithm accounts for the complex physical uncertainties (such as tilt modulation and shadowing effects) by multiplying the wave spectra (directly obtained from the radar intensity images) by an empirically determined modulation transfer function (MTF). Independent wave measurements are needed to calibrate the reconstruction and determine the optimal values of the empirical parameters. Despite the challenge, significant advances have been achieved in the use of incoherent X-band marine radars to recover surface wave information (Nieto Borge et al. 2004;Blondel and Naaijen 2012;Naaijen and Blondel 2012). Various schemes have been developed to improve the accuracy of the inversion schemes for the common marine radar, including the use of supplementary buoy measurements and coherent (Doppler) radars (e.g. Cornell et al. 2015). Since coherent radars provide the Doppler shift as well as the backscatter strength of the scatters and may operate at different polarizations, they offer more information to be included in any wave inversion scheme. Qi et al. (2016) improved the reconstruction method based on concurrent phase-resolved wave field simulations. By maximizing the correlation between the reconstructed and simulated wave fields over time, optimal values of the reconstruction parameters are obtained. Despite the meaningful improvements in the consistency and fidelity of the reconstruction, the optimal values of the parameters involved in this method vary with the wave conditions and sampling domain. Moreover, in practical applications, the radar return data may be degraded due to physical and nonphysical conditions. In these cases, the optimized reconstruction parameters may not be easy to obtain and/or physically reasonable. Despite the attractiveness of using the incoherent marine radar for wave measurements, it remains a challenging task to develop an effective inversion algorithm (that is also convenient to be implemented in practical applications) for the reconstruction of the phase-resolved wave field from the radar return intensity data.
In this paper, we propose effective improvements to the standard reconstruction method (Young et al. 1985;Nieto-Borge et al. 2004) for the wave inversion of incoherent X-band marine radar return data. The key modification to the existing method is to manipulate the radar return data by transforming the image intensity from the non-negative range (used in the standard method) to a "wave elevation"-like representation with a near zero mean. The optimized value of the parameter controlling the mean position is obtained by maximizing the Fig. 1 Illustration of the proposed modification: the radar image intensity ( , t) after range and azimuth correction (top) and its centered modification ∼ ( , t) (bottom) 1 3 consistency and fidelity of the reconstructed phase-resolved wave fields using the simulated synthetic radar images for a broad range of sea states. In addition, the use of zero-padding technique is employed to minimize the effect of spectral leakage associated with the non-periodicity in time and/or space of the radar return data in the Fourier analysis. Both synthetic and real radar data are used to verify the effectiveness and demonstrate the enhancements of the improved method verse the standard method in terms of the consistency and fidelity of the wave reconstruction.
The rest of the paper is organized as follows. In Sect. 2, we summarize the algorithms in the generation of simulated multi-directional irregular linear wave field and the associated synthetic radar images. We briefly introduce the standard reconstruction method and describe the proposed modifications in Sect. 3. The verifications and comparisons of the performance of the improved method in the reconstruction of phases-resolved wave fields and the associated spectra with the reference to the standard method are given in Sect. 4. Brief conclusions are drawn in Sect. 5.

Numerical simulation of linear wave field
Assuming the linear wave model, the simulated sea surface elevation ( , t) in the polar coordinate = (r, ) and temporal coordinate t , is specified by N × L propagating spectral wave components, containing N frequencies i , i = 1, … , N and L propagation directions j , j = 1, … , L: where Φ ij are the random phases drawn from a uniform distribution on the interval [0, 2 ] , = k x , k y is the wavenumber vector related to the corresponding wave frequency by the dispersion relation for linear waves: Here, H is the water depth, and g is the gravitational acceleration. The amplitudes of wave components (with frequency increment Δω and direction spreading increment Δθ) in (2) are obtained from a directional JONSWAP spectrum S ( , ) . The directional spectrum S ( , ) of the sea is modeled as the energy spectral density function multiplied by a normalized angular spreading function D( ) . There are many functional models in the literature to represent the directional spreading of the sea surface. Here, we model the directional distribution using a normal (1) frequency-independent spreading function of the following form (e.g. Sand 1984): where A 0 is a scaling parameter, 0 is the mean wave direction, is the standard deviation. In this paper, the largest value of spreading we consider is = 20° for wind-generated waves and = 5° for swell.

Shadowing simulation model
In this paper, synthetic images are constructed by considering only the effect of geometric shadowing which is the main effect in radar mechanism. A simple geometric criterion is applied: a sea surface elevation at position r is visible by the radar if the radar ray, which scans the surface elevation at r, if it does not intersect any surface elevation at all distances less than r (e.g. Wijaya 2016). In obtaining synthetic radar image intensity ( , t) (Qi 2012;Qi et al. 2016), we set the values of shadowing points equal to zero, change the range of not-shadowing points to a non-negative range and rescale the values into 256 Gy levels. Assuming that SAR image intensity fluctuations are proportional to the wave-height fluctuations (Vesecky et al. 1981) and the reconstruction algorithm consists of only linear operators and filters, the rescaling of the obtained non-negative values into 256 Gy levels will not affect the features of the result.

The standard reconstruction method and proposed modifications
For convenience in the application of FFT (Fast Fourier transform), we choose a rectangular subdomain D 0 of interest in the radar scanned area. We denote the image intensity in the subdomain as The objective is to reconstruct the time-varying sea surface wave elevation in the selected rectangular area D 0 from the sequence of X-band marine radar images ( , t).

Standard reconstruction method
The existing standard reconstruction method is based on the work of Young et al. (1985) and Nieto-Borge et al. (2004). The reconstruction process contains the following main steps of procedures.
Step 1: Pre-filtering of radar image for correction of range and azimuth dependences For an incoherent radar, the radar image intensity (r, , t) generally decreases with increasing the range r due to the decay law of the received electromagnetic energy with the distance and increasing deviation between the azimuth and the upwind direction w (Qi 2012;Qi et al. 2016). To account for these effects, the radar image intensity (r, , t) is multiplied by a calibration function C(r, ): In this study, the wave analysis is carried out in a set of rectangular windows near the up-wind or up-wave directions ( ≈ w ) and the range of sample image-window center is not large ( r∕R < 0.5 , where R is the maximum radius in the original radar image sequence). Under these assumptions, according to Qi (2012), the calibration function (5) takes the form: where is a positive parameter optimized by Qi (2012). We note that this step is not applied to the synthetic data.
Step 2: Extraction of wave-related signals from radar images The image sequence is transformed into a spectrum I k x , k y , by applying 3D FFT of C(r, ) ( , t) . The waverelated signal is extracted from the 3D image spectrum by the use of filtering. The high-pass filter is used to eliminate the long-range dependence modulation effects (Nieto-Borge et al. 2004): where c is the reconstruction parameter, and Δ is the angular frequency increment in the analysis. The reconstruction parameter c is determined empirically. The band-pass filter is applied to extract the wave-related components based on the linear dispersion relationship (3):  where b is also an empirical parameter and we use b = 2 in this study (Qi 2012).
Step 3: Multiplication of the modulation transfer function M( ) To minimize the complex non-physical effects inherent in the X-band radar measurements such as shadowing and tilt modulations, Nieto-Borge et al. (2004) proposed to multiply the resulting wave-related spectrum I k x , k y , by a simple modulation transfer function (MTF) The empirical parameter q = 0.5 is chosen in this study based on the optimization proposed by Wijaya (2016).
Step 4: Production of phase-resolved sea surface To obtain the unscaled reconstructed phase-resolved sea surface elevation ∼ R ( , t) , the inverse 3D FFT is applied on the filtered image spectrum. The calibrated (scaled) estimation of the sea surface elevation R ( , t) is then given by: where the scaling factor c 0 is defined by with ∼ R being the standard deviation of the reconstructed wave elevation ∼ R ( , t) . The significant wave height H s could be derived from the signal-to-noise ratio of the radar backscatter (Nieto-Borge et al. 2004), or based on the shadowing statistics proposed by Wenzel (1990), or using independent wave measurements (e.g. point buoy measurements as used in the present study).

Proposed modifications
To increase the reconstruction accuracy, we propose the following manipulation with the radar image intensity ( , t) . Right before starting Step 2 of the procedure, we replace the original image intensity ( , t) in the visible region by the modified image intensity ∼ ( , t) that takes the form: Here is the empirically estimated parameter; s(x, y, t) is the shadowing function that is equal to 0 (or 1) if the point (x, y) is shadowed (or visible). In the shadowed region, (12) is applied to center the image data around the zero mean  and to transform the radar image intensity from non-negative range [0, A] to a "wave elevation"-like representation (Fig. 1). Thus, the zero level of ∼ ( , t) closely corresponds to the mean sea level of the original wave elevation ( , t) . The optimal value of is determined based on the optimization of reconstruction accuracy for the synthetic radar data, as described in Sect. 4.1. It is important to emphasize that the modification of radar image intensity in (12) is applied to ( , t) in the visible region only, but not in the shallowed region.
Performing FFT in a discrete finite domain for a non-integer number of periods of a signal introduces artificial effects such as spectral leakage when the 'power' of the original frequencies will 'leak' into a number of different frequency bins (Bracewell 1986). In this case, the Discrete Fourier Transform (DFT) assumes that its input is one period of an infinitely periodic function and, if the values at the opposite edges are not the same, then these edge discontinuities introduce spurious frequencies into the spectrum. When we consider a finite sequence of radar images, the spectral leakage in the frequency domain occurs due to the truncation of the infinitely long signal to a finite one. Spectral leakage cannot be eliminated but it can be reduced using window functions or zero padding. In this study, we shall demonstrate the effectiveness of using zero padding for reconstruction accuracy improvement. In the zero-padding algorithm, the sequence of N t images is extended by adding 'zero'images (zero-padded) before it is transformed by the DFT. Zero-padding a signal in time domain does not increase the amount of information that is contained in the signal. But it shows an interpolation effect on the spectrum and allows to apply filtering more accurately since reducing the spacing between bins may put a bin closer to the true frequency and give a better estimation of the amplitudes and frequencies.
Generally, spectral leakage occurs in both time and space domains. In this study, zero-padding is applied only in the time domain, but not in the space domain since the observed edge effect (e.g. reconstruction error near the border of the domain) in the space domain is not as severe as in the time domain. This is due to the fact that with the use of the large spatial domain, the wavenumber resolution is sufficient for the reconstruction methods. In the case of relatively small spatial domains, apparent edge effects have been observed to propagate in the main wave direction. The use of windowing function (or zero-padding) technique can remove such edge effects in the spatial domain . We remark that the success of using the radar intensity images for phase-resolved wave field reconstruction has been demonstrated in a number of studies, the robust application of this approach under general conditions remains challenging. A major fundamental and technical issue with this approach is the complexity associated with the  parameterization of the empirical MTF used in Step 3, which depends on a large number of factors, ranging from wind/ wave conditions and radar parameters to selected analysis domain parameters, such as range, azimuth and domain size. There have been some efforts on improving the radar technology on wave measurements such as the application of Doppler radar, with which one obtains the motion of an area of sea surface and thus avoids the use of MTF (e.g. Connell et al. 2015;Støle-Hentschel et al. 2018). The purpose of the present study is to develop an improved data processing algorithm that can reduce the non-physical effect associated with the shallowing factor. This would reduce the sensitivity of MTF on the range parameter in wave reconstruction. Another important technical issue is related to the fact that the radar backscatter modulations depend on the look direction relative to the wave propagation (e.g. Lund et al. 2014). This makes it a challenge to accurately reconstruct a complex wave field with wave components propagating in orthogonal directions with the use of a single rectangular analysis domain, as adopted in this study. In this case, multiple analysis domains corresponding to different dominant wave/viewing directions may need to be selected and combined to obtain a wave reconstruction with reasonable accuracy. One other solution is to use the polar Fourier transform method in which each spectral component is estimated from an angular region centered on that component (Lyzenga 2017).

Wave field reconstruction from synthetic radar images
Synthetic linear wave fields have been generated in polar coordinates with a spatial ( Δr = 3.5 m) and temporal ( Δt = 2.0 s) resolution and domain size that is realistic for wave radar systems. The considered wave conditions are classical North-sea conditions according to NORA10 (Norwegian ReAnalysis 10 km) data for the period from September 1957 to September 2016. We consider one case of unimodal (wind waves) and three cases of bimodal sea state in which a combination of swell and wind-waves is present within the same area of observation. The sea-state parameters used in the generation of simulated wave fields are summarized in Table 1.
For the radar height, we take as example H r = 30 m above the still water level, which is a reasonable value for large ships as oil and gas tankers (Van Groesen et al. 2017). Figure 2 shows a sample frame of the simulated wave elevation ( , t) and the corresponding synthetic radar image. Since the radar image intensity is non-negative, the non-shadowing values of ( , t) have been scaled to a positive range while the values of shadowing points have been set equal to zero.
A rectangular subdomain of 1500 × 1500 m 2 in the dominant up-wave directions of swell and sea is selected from N t = 32 sequential simulated radar images. Since the examples considered are for the bi-modal seas with relatively small angles between dominant directions of swell and sea, satisfactory wave field reconstructions are obtained by the use of a single large analysis region. The intensity values at equally spaced N x × N y = 512 × 512 Cartesian grid points in the selected subdomain are obtained from the values at the equally spaced polar grids using the linear interpolation method.
In this study, we use the 2D correlation coefficient and normalized point-to-point absolute error as a performance The normalized point-to-point reconstruction error is defined as follows: Corr t i .
(18) Figure 3 shows the variation of the average correlation coefficient Corr mean with respect to the parameter . We can see that Corr mean gradually increases with increasing , and reaches its maximum at ≈ 0.87 and then decreases. The value = 0 corresponds to the standard reconstruction method. At = 0.87 , the zero level of the modified radar image intensity (12) approximately corresponds to the sea level of the original waves. The fact that the correlation coefficient reaches its maximum at < 1 could be explained by the fact that the average image intensity 0 obtained from (13)-(14) does not consider the effects of unknown shadowed points. It turns out that the inclusion of the sea level information to the radar return data is found to significantly improve the consistency and fidelity of wave field reconstruction.
To assist in understanding the fundamental effect of shifting the mean level of radar return intensity on the wave field reconstruction, we illustrate the influence of the shadowing function s(x, y, t) that is introduced in Sect. 3.2. The radar return signal corresponds to the surface elevation given by [ (x, y, t) + a]s (x, y, t) where the constant a = min{ (x, y, t) ⋅ s(x, y, t)} . Among the radar return signal, the term (x, y, t) ⋅ s(x, y, t) contains the wave information while the term a ⋅ s(x, y, t) contains the shadowing effect irrelevant to the waves. As an illustration, Fig. 4 compares the spectral densities of (x, y, t) ⋅ s(x, y, t) and a ⋅ s(x, y, t) as well as the benchmark wave field (x, y, t) for Case 3. It is seen that the spectrum of the shadowing Fig. 11 Comparisons of the instantaneous three-dimensional wave surfaces of the reconstructed and original wave fields. The results are for Case 4 function s(x, y, t) also contains the wave components that satisfy the dispersion relation, which are actually nonphysical. Depending on the value of a , these non-physical wave components may be comparable or even larger in magnitude than those physical wave components contained in (x, y, t) ⋅ s(x, y, t) . The wave field reconstructed by the standard method, based on the direct radar intensity data [ (x, y, t) + a]s(x, y, t) , would consist of the wave components from both the physical effect of (x, y, t) ⋅ s(x, y, t) and the non-physical effect of a ⋅ s(x, y, t) . By lowering the radar return intensity in visible regions by β 0 in (12), the modified method minimizes the non-physical effect of the shadowing term a ⋅ s(x, y, t) . As a result, the modified method can achieve a much-improved accuracy and fidelity of wave reconstruction. We note that since the shadowing depends on the range and wave conditions, the modified method may be further improved by considering the parameter β as a function of the range or the percentage of shadowed area in the wave field. This is under investigation.
In Figs. 5 and 6, we compare the reconstructed wave fields by the standard method (with = 0 ), the modified method (with the optimal value of = 0.85 ), and the benchmark original wave field for the sample wave profiles and surface at the different stages of wave reconstruction. The results shown are obtained for Case 3. The reconstructed wave fields by the modified method at all stages of reconstruction agree better with the original wave field Fig. 12 Comparisons of the wave profiles (along the y-axis at x = 900 m) of the reconstructed wave fields by the standard method (with = 0 ), the modified method (with optimal = 0.85 ), and the modified method (with optimal = 0.85 ) and zero-padding, and the benchmark original wave field ("Original"). The results are for Case 2 than those by the standard method. The comparisons in Figs. 5 and 6 demonstrate the improvement of accuracy in wave reconstruction by the modified method.
The new parameter is used to bring the mean level of radar return intensity down to the realistic mean sea level to reduce the non-physical shadowing effect in wave reconstruction. From the concern of this purpose, we think it should be more effective to consider as a function of the percentage of shadowed area in the wave field, instead of multiple parameters such as range, significant wave height, radar position etc.
In the implementation of the zero-padding technique for the reduction of spectral leakage associated with the finitetime truncation of the infinitely long time signal, the input data of 3D DFT takes the size of N x × N y × (N t + N 0 ) , where N 0 is the number of extended zero images. Figure 7a compares the correlation coefficients in the reconstruction of wave fields for a duration of 64 s (with 32 radar images) obtained by the standard method, and the modified method with and without zero-padding. Much lower values of the correlation coefficient are seen near the beginning and ending time of signal for all four cases considered. The use of zero-padding is shown to be quite effective in removing the edge effect in all the cases. Figure 7b shows the variation of the mean value of correlation coefficient Corr mean with the number of zero-padding images. It is seen that the minimum required length of zero-padding is about N 0 = 5 for all the synthetic wave fields tested. Table 2 summarizes the results of maximum and average correlation coefficients in the reconstruction of the four synthetic sea states by the standard method and presents modified method with and without zero-padding. Figure 8 displays the point-to-point reconstruction errors in the selected rectangular analysis domain at three representative instants (corresponding to the initial, middle, and end time of the image data) for Case 2, which are obtained by the standard method and present modified method with and without zero-padding. The similar results for Case 4 are shown in Fig. 9. The comparisons between the results by the different methods show that the modified method significantly improves the accuracy of the phase-resolved wave field reconstruction from X-band radar images from the standard method. In addition, the employment of the zero-padding technique effectively removes the spectral leakage associated with the endpoint discontinuity. To examine edge effects, the two-dimensional images of the input wave fields for selected frames are presented in Figs. 10 and 11 along with the estimated wave fields for the cases shown in Figs. 8 and  9. The examples of one-dimensional plots versus range corresponding to the initial and middle time of the image data are shown in Figs. 12 and 13. Fig. 13 Comparisons of the wave profile (along the y-axis at x = 900 m) of the reconstructed wave fields by the standard method (with = 0 ), the modified method (with optimal = 0.85 ), and the modified method (with optimal = 0.85 ) and zero-padding, and the benchmark original wave field ("Original"). The results are for Case 4 1 3

Wave field reconstruction from real radar images
All data used for the comparison in this case were collected at Solstad IMR subsea construction vessel Normand Ocean, operated by DeepOcean. At Normand Ocean, a Furuno X-band radar is mounted at a height of 20 m above the surface. The radar has a slightly more than 180° unobstructed view in the forward direction. Close to the radar, Miros RangeFinder is mounted looking downwards right ahead of the bulb of the ship. This gives an independent measurement of the incoming waves and mitigates the uncertainties in the significant wave height estimation. We consider a radar image sequence generated when the radar is working in short pulse. Figure 14a-d displays one representative frame of the raw radar images for each of the four different sea states measured. In these figures, the color scale indicates the strength of radar image intensity, and Hm0 denotes the significant wave height evaluated from the spectrum of the measured waves.
After estimating the directional spectra, the images were rotated to have the main energy propagation along the positive y-direction. The region with the distance to the radar (set as the origin) less than 600 m is the blind area where the radar backscatter is too high to be useful due to specular scattering from the sea surface. The upper border of the probe area is chosen to be 1000 m from the origin. Then we obtain a time series of consecutive radar images in the rectangular subdomain of 1500 × 1500 m 2 with N x = N y = 512 nodes in the x-and y-directions. For better accuracy in wave reconstruction, the rectangular analysis domain is chosen to face the dominant up-wave direction. The sampling time in the time series was assumed to be 2.5 s (the antenna rotational period was varied from 2.47 to 2.53 s). The optimal value of β = 0.85 is used for the modified method in the construction of wave fields based on these realistic radar intensity data. In the application of 3D FFT, the radar intensity image in the rectangular subdomain is assumed to correspond to the snapshot of sea surface. In practice, the intensity image is produced by scanning while the radar's antenna continuously rotates. With the antenna rotational period of about 2.5 s, the snapshot assumption would lead to a maximum of about 0.3 s time difference in the time-varying phase of waves in the subdomain, which is not expected to significantly affect the accuracy of wave reconstruction. Figure 15 shows the comparisons of the spectral density functions as a function of wave frequency obtained by the standard and modified reconstruction methods for the four real sea states. The spectral density functions by the standard and modified reconstruction methods are based on a 160 s of radar measurements on a large spatial area. The integration of the 2D wavenumber spectrum obtained from reconstruction gives the 1D frequency-dependent spectral density function. The result from the independent Miros sensor measurement at a single space point (vessel position) is also shown Limited by the time resolution of X-band radar measurement with ~ 2.5 s time interval per radar image, the reconstruction methods are unable to capture short waves with periods smaller than 4-5 s. Compared to the standard method, the modified method shows a significant improvement in reconstruction of wave spectra with relatively broader bandwidth, as the comparisons in Fig. 15a, d indicate. Also, the modified reconstruction method shows a much better agreement with the independent Miros point measurement for the bimodal sea state, as in Fig. 15b, which contains a combination of wind-sea and swell. In this case, the standard method largely overestimates the amplitudes of swell components. For the sea state from 20 April 2019 (Fig. 15c), when the frequency band is relatively narrow, the shape of the reconstructed wave spectrum by the standard method matches that by the modified method, both of which compare well with the independent Miros point measurement. Figures 16, 17, 18 and 19 display one sample frame of the raw data image and its corresponding reconstructed wave elevation in the selected subdomain for the four real measured sea states. Unfortunately, it is not possible to estimate the reconstruction accuracy for individual waves from

Conclusions
In this study, we develop and verify an improved phaseresolved wave inversion method for the reconstruction of wave elevation map and its spectral components from incoherent X-band marine radar images. The method extracts the gravity wave signals that satisfy the linearized dispersion relation. The key modification to the existing standard method is based on the radar image intensity transformation from the non-negative range to a "wave elevation"-like representation by centering the intensity data around the zero mean. A parameter optimization of the proposed modification is performed using the simulated synthetic radar images. The quantitative comparisons of phase-resolved wave reconstruction by the standard and modified methods are made for both simulated and real radar wave data. The correlation between the reconstructed and the true wave elevation by the present modified method is enhanced to the range of 0.9-0.92 for the simulated radar data compared to the value of about 0.7 by the existing standard method. The wave spectra obtained from the real radar images of sea surface by the modified method are in good agreement with that from the independent Miros RangeFinder point measurement in both unimodal and bimodal sea states. In particular, the modified method is capable of obtaining an accurate reconstruction of wave spectra of broad-banded and bimodal seas while the stand method seems to be effective only for relatively narrow-banded unimodal seas.