Developing a Quality Flag for the SAR Ocean Wave Spectrum Partitioning with Machine Learning

1. Introduction

Copernicus, the Earth observation component of the European Union’s Space program, plays a pivotal role in monitoring our planet and its environment for the benefit of European citizens [33]. Within this program, the S1 mission, operated by the European Space Agency (ESA), is equipped with SAR satellites. The S1 satellites, which use C-band SAR (central frequency 5.405 GHz), are capable of collecting data across multiple acquisition modes, with WV predominantly employing vertical-vertical polarization (VV) to observe open ocean conditions [35]. The launch of Sentinel-1A occurred on April 3, 2014, followed by Sentinel-1B on April 25, 2016. Both satellites operated in the same orbital plane with a phase shift of 180°. Sentinel-1B stopped delivering data in December 23, 2021, was decommissioned in 2023 and replaced by Sentinel-1C (launched in December 5, 2024). The data obtained by Sentinel-1 in WM are routinely processed into L2 OCN products and include directional ocean wave spectra generated using a quasi-linear inversion method [17]. This inversion process provides valuable information on the distribution of wave variance in frequency and direction. However, SAR-based inversion is mainly applicable to long ocean swells (with wavelengths greater than 200 m) [1,5,8] and to wave propagation in the marginal ice zone [3,25].

Although integral parameters from SAR-derived spectra can effectively characterize a single-swell system, they become less informative when multiple, overlapping wave systems are present (Figure 1). In such cases, it becomes essential to partition the ocean wave spectrum into independent wave systems to ensure accurate data analysis and quality control. Within the L2 OCN product, the ocean wave spectrum is divided into a maximum number of five distinct partitions [17] using the "watershed" algorithm, originally introduced by [4,13]. Each partition is described by its corresponding integral parameters, including the partitioned Significant Wave Height (Hs), the peak wavelength, and peak wave direction (Figure 1).

A direct comparison between the SAR-derived wave partitions and those obtained from the WW3 model [28] reveals discrepancies, particularly between the dominant and secondary swell systems. Differences in partitioned wavelengths exhibit positive clustering in the azimuth direction and negative clustering in the range direction [18,29,31]. Notably, a study by [15] in the Australian region shows Root-Mean-Square Errors ($RMSE$) of 0.70 m for Hs, 0.9 s for mean periods, and 30° for wave direction. These discrepancies are largely attributed to the motion of ocean waves with radial velocity components, which induce Doppler frequency shifts in radar echoes. As a result, SAR image points with radial velocity components are displaced along the azimuth direction, leading to non-linear effects when converting SAR data into ocean wave spectra. To address these challenges, ML techniques are commonly employed to handle such non-linearities. While existing methods have corrected some integrated SAR wave parameters (e.g., total Hs, as discussed by [21]), they have yet to correct individual ocean wave systems.

This study applies data-driven ML techniques to define a prior QF for each of the five wave partitions derived from Sentinel-1 WV OCN products. The QF is determined by learning partition error parameters, including relative errors in partitioned significant wave height ($\Delta \mathit{Hs}$),

$$\Delta Hs=2{\displaystyle \frac{\left(H{s}_{\mathrm{SAR}}-H{s}_{\mathrm{WW}3}\right)}{H{s}_{\mathrm{SAR}}+H{s}_{\mathrm{WW}3}}}$$

(1)

the relative error of the peak period ($\Delta \mathit{T}$),

$$T_{\mathrm{SAR}}=\sqrt{{\displaystyle \frac{w{l}_{\mathrm{SAR}}}{1.56}}},T_{\mathrm{WW}3}=\sqrt{{\displaystyle \frac{w{l}_{\mathrm{WW}3}}{1.56}}}$$

(2)

$$\Delta \mathit{T}=2{\displaystyle \frac{\left(T_{\mathrm{SAR}}-T_{\mathrm{WW}3}\right)}{T_{\mathrm{SAR}}+T_{\mathrm{WW}3}}}$$

(3)

and the absolute error of the peak wave direction ($\Delta \varphi$),

$$\Delta \varphi ={\varphi }_{\mathrm{SAR}}-{\varphi }_{\mathrm{WW}3}$$

(4)

based on collocated SAR and WW3 partition pairs. By transitioning from absolute to relative errors, this data-driven approach ensures a more balanced distribution of integral parameters across the QF classes. This adjustment mitigates the disproportionate influence of higher $\mathit{Hs}$ values on error estimation, thereby reducing bias and providing a more equitable representation of partition quality.

The primary contributions of this research are: (1) the application of ML algorithms to improve the qualification of S1 ocean wave spectrum partitioning, (2) an explanation of the model’s predictions and how they relate to the geophysical inversion process using SAR data, and (3) the establishment of a framework for future research aimed at advancing the retrieval of ocean wave spectra from SAR.

The manuscript is structured as follows: Section 2 presents the dataset utilized in this study. Section 3 outlines the full ML methodology used. In Section 4, we present the results and validate the model. Section 5 offers an in-depth discussion of the findings, and Section 6 concludes the paper with future research perspectives.

2. Data Set Collection

2.1. Sentinel-1 WV OCN OSW Products

The S1A/B WV mode operates by acquiring SAR images, or imagettes, spaced 100 km apart, with two alternating incidence angles, WV1 and WV2, along the satellite’s flight path using a "leap frog" acquisition strategy [24]. In each acquisition, the S1 WV mode can only operate in a single polarization, VV or HH, with VV polarization being the default for global ocean observations. Each acquired imagette covers an area of approximately 20 x 20 km, with a high spatial resolution of 5 meters. The ocean wave spectrum (OSW) component of the WV OCN product is derived from Level-1 (L1) Single Look Complex (SLC) SAR imagettes, processed through the ESA L2 spectral inversion unit [17].

The inversion process begins with the estimation and removal of non-linear contributions from the SAR data, which is essential for improving the accuracy of the final results. Once the non-linearity is treated, the cross-spectra technique is applied to resolve the 180° directional ambiguity in wave propagation, which is crucial for obtaining reliable directional information. This method allows the retrieval of the ocean wave spectrum with a directional resolution of 10° and 30 exponentially spaced wave numbers, spanning from 30 to 800 meters. The final step of the process involves partitioning the 2D ocean wave spectrum into five independent ocean wave systems using the watershed algorithm, after which the integral parameters (such as Significant Wave Height, peak wavelength, and peak wave direction) for each partition are computed.

For the purposes of this study, we use data acquired exclusively from the S1A satellite, covering both WV1 and WV2 modes with VV polarization, collected between December 2024 and February 2025. This time period was selected because it corresponds to the deployment of version 3.90 of the S1 L2 Instrument Processing Facility (IPF), which includes a new SAR Modulation Transfer Function (MTF) tuning. This upgrade significantly improves wave inversion performance from SAR imagery [17].

Additionally, to ensure the integrity of the ocean wave data, products acquired at latitudes above 55° North and South are filtered out to exclude any possible contamination from sea ice, which could distort wave measurements.

The publicly available S1 OCN WV data used in this study can be accessed through the Copernicus Open Access Hub [34].

2.2. WW3 Hindcasts

The WW3 Hindcast, part of the Integrated Ocean Waves for Geophysical and other Applications (IOWAGA) project led by IFREMER, is a third-generation wave model that resolves the spectral action density equilibrium equation for wavenumber-directional spectra. It is driven by data on 3-hour wind and ice concentration from the European Center for Medium-Range Weather Forecasts (ECMWF) product [2,22]. The model operates with spatial and temporal resolutions of 0.5° and 3 hours, respectively. The spectral data are organized into 32 frequency bins, exponentially spaced from 0.038 to 0.7159 Hz, and 24 directional bins, each separated by 15°. To partition the wave spectra, the model employs the “watershed” algorithm, which divides the spectra into five partitions corresponding to different swells, as well as one partition for wind seas. Wind seas are represented by Partition 0, while Partitions 1 through 5 correspond to various swell systems. The accuracy of the model’s wave parameter has been validated against buoy and altimeter measurements, showing strong agreement [2,9,23]. The WW3 Hindcast data is available for download through the IFREMER FTP server (ftp://ftp.ifremer.fr/ifremer/ww3/HINDCAST/).

2.3. Match-Up Partitions

The S1 OCN WV OSW products were collocated with the WW3 hindcasts based on the closest spatial and temporal grid points, within a window of 0.25° in latitude/longitude and 1.5 hours in time. To ensure the best matchups at the partitioned swell level, these spatio-temporal collocations were further refined by identifying the closest partition matches. For this partition cross-assignment, the spectral distance $D_{spec}$ between the partitions was defined, as proposed by [14], according to the following equation:

$$D_{\mathrm{spec}}={\displaystyle \frac{1}{q}}\left(\left|{\varphi }_1-{\varphi }_2\right|mod360+2·r·{\displaystyle \frac{\left|T_1-T_2\right|}{T_1+T_2}}\right)$$

(5)

In this equation, ${\varphi }_i$ and $T_i$ represent the peak propagation direction and peak period of each swell system, respectively. The constant r was set to 250, which scales the period error appropriately: a 20° directional error is treated as equivalent to an 8% error in wave period, reflecting the typical SAR swell measurement accuracy. Additionally, the value of $q=30$ was chosen to define the error thresholds to ensure 30° directional error and a 12% in period error to give a spectral distance equal to one [30].

No further filtering was applied to the dataset to preserve a comprehensive representation of the various errors inherent in the inversion and partitioning processes, as well as any non-geophysical contributions present at the ocean surface.

3. Method Details

3.1. Quality Flag Definition

The Normalized Radar Cross-Section (NRCS) is influenced by several error sources, including geometric distortions, sensor-specific acquisition errors, and the presence of non-geophysical factors such as atmospheric influences and surface roughness at the air-sea interface. Additionally, the quasi-linear inversion process employed to derive ocean wave spectra from SAR data introduces further inaccuracies. The partitioning algorithm, which divides the ocean wave spectrum into distinct wave systems, is also subject to its own set of errors. To address these complexities, the first critical step in our algorithmic approach is to independently quantify the errors in $\Delta \mathit{Hs}$, $\Delta \mathit{T}$, and $\Delta \varphi$ at the partition level with respect to the WW3 model data.

These errors are estimated using a supervised ML algorithm, which leverages a set of SAR-derived observables (i.e., features) to predict discrepancies in the ocean wave parameters [7]. The features used for the error estimation process are detailed in Table 1. The computation of some features is described in detail in [17].

Regardless of the acquisition mode (WV1 or WV2), the same set of features is used to estimate the three primary errors. This necessitates the training of six separate models tailored to the S1A mission, each aimed at improving the accuracy of partition-level error quantification and ultimately refining the ocean wave parameter retrieval process.

The second step of the algorithm focuses on computing a combined error metric for each individual acquisition mode (i.e., WV1 and WV2). This combined error, denoted as $C_{\mathrm{error}}$ and defined in Eq. 6, is calculated by multiplying the three primary error estimates—$\Delta \mathit{Hs}$, $\Delta \mathit{T}$, and $\Delta \varphi$—previously inferred using the trained models applied to the validation dataset, which consists of all S1A acquisitions from February 2025. The metric is expressed as:

$$C_{\mathrm{error}}=\Delta \mathit{Hs}·\Delta \mathit{T}·\Delta \varphi$$

(6)

Since the ultimate goal is to assess the quality of ocean swell systems in a manner that is independent of the SAR acquisition mode, a unified combined error metric is produced by merging the mode-specific error estimates from both WV1 and WV2. This approach is motivated by the need for consistent and homogeneous error characterization across different SAR configurations. Since both WV1 and WV2 data are used interchangeably in operational wave monitoring and analysis, having a mode-agnostic quality metric ensures comparability between partitions retrieved under different acquisition conditions. By assembling the combined error from both modes, we enable a seamless interpretation of partition-level quality, facilitating the identification and classification of ocean swell systems without bias introduced by the acquisition geometry.

This reference error is then partitioned into five equally probable intervals using the q-quintile method, commonly referred to as Quintile. These intervals represent distinct levels of error severity, ranging from low to high. Once this reference error distribution is established, the calculated $C_{error}$ for each partition is matched to the corresponding quintile range. Based on this comparison, each partition is assigned a quality label: “Very Good,” “Good,” “Medium,” “Low,” or “Poor.” This classification serves as an indication of the reliability and accuracy of the partition, facilitating informed decision-making in further analyses.

3.2. Machine Learning Data Set

In the context of algorithm development, we focused on the S1 WV OCN dataset, which was acquired over two months : December 2024 and January 2025. These data were assigned to the training phase, with an 80%-20% split for training and validation, respectively. This approach ensures that the model is rigorously tested on unseen data during training. Data from February 2025 were used exclusively for model evaluation, providing an independent test set to assess the model’s generalizing ability.

During the training process, no data normalization was performed, as the chosen machine learning algorithm is inherently robust to raw, unnormalized data, facilitating direct input handling without additional preprocessing steps and data transformation.

The targets for the machine learning models are : $\left|\Delta \mathit{Hs}\right|$, $\left|\Delta \mathit{T}\right|$, and $\left|\Delta \varphi \right|$. We selected the absolute error to emphasize the overall discrepancy between WW3 and the SAR observations, rather than focusing on the signed difference. Given that the WW3 does not account for all geophysical complexities, we place greater trust in the observational data as a more reliable reference, thereby enhancing the confidence in the results.

3.3. Machine Learning Modeling

In this study, we initially used a traditional Random Forest algorithm; however, it soon revealed limitations in terms of interpretability. Consequently, we shifted our focus to the eXtreme Gradient Boosting (XGBoost) supervised learning algorithm, which offered improved performance and explainability. XGBoost is a scalable, distributed machine learning library based on gradient-boosted decision trees (GBDT). A decision tree is a model that makes decisions by splitting data into branches based on simple conditions, forming a set of "true or false" statements. Since its debut in 2014, XGBoost has become one of the most widely adopted algorithms among data scientists and machine learning practitioners due to its high performance, particularly for structured data problems. The library is open-source [32] and is able to efficiently train and test models on large datasets.

One of the primary reasons for selecting XGBoost in our study is its regularization capabilities. Regularization helps to control overfitting, where the model fits too closely to the training data, thereby reducing its ability to generalize to unseen data. XGBoost applies L1 and L2 penalties to the weights and biases of each tree, helping to maintain model simplicity and prevent overfitting. Additionally, XGBoost is highly optimized, with features that make it more memory efficient, such as cache-awareness, which is crucial when working with large datasets.

XGBoost also stands out for its ability to handle missing data, eliminating the need for imputation, and its ability to work with data in their raw, unnormalized state, simplifying the preprocessing process. These advantages make XGBoost an attractive choice for dealing with complex data problems.

To fine-tune the XGBoost model and achieve optimal performance, hyperparameter tuning was performed. Hyperparameters are critical settings that influence how the model learns, and selecting the right ones can significantly improve model accuracy. Given the large number of hyperparameters, finding the best combination manually is impractical. Therefore, we used the GridSearch technique [19] to systematically explore different combinations of hyperparameters. In this process, multiple sets of hyperparameters are tested within a defined search space, and the performance of the model is evaluated on a validation dataset. Although this method is effective, it can be computationally expensive, particularly as the number of parameters and their possible values increases. A list of the hyperparameters most commonly tuned in XGBoost is provided in Table 2, with additional parameters available in [32].

Machine learning models must exhibit robustness, meaning they should minimize the impact of outliers and prioritize the influence of typical data points. In tasks such as parameter estimation, using a robust loss function (e.g., absolute error) is often preferred over a non-robust one (e.g., squared error) because it is less sensitive to large deviations. Common loss functions in regression tasks are squared loss, $L\left(x\right)=x^2$, and absolute loss, $L\left(x\right)=\left|x\right|$. While squared loss is highly sensitive to outliers, making it less reliable in such cases, absolute loss is more resilient as it focuses on the order of data rather than their absolute magnitude [20].

To leverage the advantages of both, the pseudo-Huber [6,10] loss function is often used. This hybrid function blends the properties of both quadratic and absolute loss while maintaining the smooth differentiability required for optimization. It is mathematically expressed as:

$$L_{Hp}\left(x\right)={\delta }^2\left(\sqrt{1+{\displaystyle \frac{x^2}{{\delta }^2}}}-1\right),$$

(7)

where $\delta >0$ is the threshold parameter that controls the transition between the quadratic and linear behavior of the loss.

For small values of x, the pseudo-Huber loss behaves like a quadratic function,

$$L_{Hp}\left(x\right)\approx {\displaystyle \frac{x^2}{2}},$$

while for large values of x, it approximates the absolute loss,

$$L_{Hp}\left(x\right)\approx \delta \left|x\right|-{\delta }^2.$$

In XGBoost, the model is optimized using gradient-based methods, such as gradient descent, to minimize the chosen loss function while incorporating regularization to prevent overfitting. The model’s objective function combines the loss function and a regularization term, as shown in Eq. 8:

$${\mathcal{L}}^{\left(t\right)}=\sum _{i=1}^{n}l\left(y_i,{\widehat{y}}_i^{(t-1)}+f_t\left(x_i\right)\right)+\Omega \left(f_t\right).$$

(8)

Here, l represents the loss function (e.g., pseudo-Huber), $f_t$ denotes the prediction from the t-th tree, and $\Omega \left(f_t\right)$ is the regularization term that helps avoid overfitting by penalizing overly complex models. The gradient descent process iteratively updates the model parameters to minimize this objective. Additional details on regularization can be found in [7].

4. Results

4.1. Metrics Definition

Model training was conducted using an 8GB NVIDIA RTX A2000 graphics card. Thanks to XGBoost’s ability to leverage parallel processing, the grid search and model training time for each model was approximately three hours. The performance of the XGBoost models in predicting the errors in S1 ocean wave spectrum partitioning parameters and controlling the final partitioning labeling was assessed using several evaluation metrics. These metrics include the $RMSE$, Normalized Root Mean Squared Error ($NRMSE$), standard deviation ($STD$), Median Absolute Error ($MAE$), Scatter Index ($SI$), coefficient of determination ($R^2$), and explained variance score ($EVS$).

The $RMSE$ is a widely used metric that measures the square root of the average of squared differences between predicted and observed values. It is sensitive to large errors, thus providing a strong indication of how far off the predictions are from the actual values. The RMSE is defined as:

$$RMSE=\sqrt{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{(y_i-\widehat{y_i})}^2}$$

(9)

where $\widehat{y_i}$ is the predicted value obtained with the model, and $y_i$ is the true value.

The NRMSE is a normalized version of the RMSE that scales the error by the range or the mean of the observed values. This allows for a comparison of error across different datasets with different ranges or units. It is calculated as:

$$NRMSE={\displaystyle \frac{RMSE}{\mathrm{Range}\mathrm{or}\mathrm{Mean}\mathrm{of}{y}_i}}$$

(10)

The STD measures the amount of variation or dispersion of a set of values. In the context of model evaluation, it is useful for understanding the variability in the errors and is given by:

$$STD=\sqrt{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{(y_i-\mu )}^2}$$

(11)

where $\mu$ is the mean of $y_i$.

The MAE measures the median of the absolute differences between the predicted and true values, making it a robust metric to outliers, as it gives equal weight to each error regardless of magnitude. The MAE is defined as:

$$MAE=\mathrm{median}\left(|y_i-\widehat{y_i}|\right)$$

(12)

The SI is a metric used to assess the consistency of the model’s predictions by calculating the ratio of the RMSE to the mean of the observed values. It provides insight into the relative error of the predictions, and is given by:

$$SI={\displaystyle \frac{RMSE}{\mathrm{mean}\mathrm{of}{y}_i}}$$

(13)

The R² is a statistical measure that indicates how well the predicted values match the actual values. It measures the proportion of the variance in the dependent variable that is predictable from the independent variables. An R² value of 1 means perfect predictions, while a value of 0 indicates no correlation between predicted and actual values. R² is computed as:

$$R^2=1-{\displaystyle \frac{{\sum }_{i=1}^n{(y_i-\widehat{y_i})}^2}{{\sum }_{i=1}^n{(y_i-\overline{y})}^2}}$$

(14)

The EVS is similar to R² but with the key difference that it does not penalize for systematic offsets in the predictions. It measures the proportion of variance in the target variable that is explained by the model, and is defined as:

$$EVS=1-{\displaystyle \frac{\mathrm{Var}(y-\widehat{y})}{\mathrm{Var}\left(y\right)}}$$

(15)

$$\mathrm{Var}\left(y\right)\mathrm{is}\mathrm{the}\mathrm{variance}\mathrm{of}\mathrm{the}\mathrm{true}\mathrm{values},\mathrm{and}\mathrm{Var}(y-\widehat{y})\mathrm{is}\mathrm{the}\mathrm{variance}\mathrm{of}\mathrm{the}\mathrm{residuals}.$$

4.2. Overall Models Performances

The results for the two acquisition modes, categorized by partition class (i.e., quality flag), are derived from the validation dataset using the partition-effective Hs, peak wavelength, and peak wave direction with respect to WW3, and are presented in Table 3, Table 4, and Table 5, respectively. This analysis provides a clear view of how the model distinguishes between good and bad partitions.

Performance metrics for each class, evaluated based on partition parameters and compared to the WW3 model, reveal insightful trends. For each parameter, the model demonstrates strong initial agreement with WW3. Specifically, for the ’very good’ class, the agreement is excellent, with an R² value of 0.77 (overall mean). However, as the class quality decreases towards ’poor,’ the agreement progressively weakens. In the ’poorest’ class, the distributions of SAR and WW3 diverge significantly, reflecting a clear mismatch between the two. This trend holds true for almost the entire range of partition wave parameters.

For the subsequent analysis, we will focus exclusively on the extreme classes—"very good" and "poor"—along with the "medium" class, as these categories represent the most relevant distinctions in partition quality for our purposes.

To showcase the robustness of the methodology in partition classification, from high-quality (Very Good) to low-quality (Poor) classes for both WV1 and WV2 across partition parameter ranges, scatter plots are presented in Figure 2 for the case of Hs. In addition to the clear class separation, the plots reveal a consistent trend of increasing error as partition quality declines, with a notable mismatch in the Poor class. This mismatch underscores the model’s ability to effectively capture variations in error distribution, reinforcing the reliability of the classification approach. These findings highlight the model’s strength in preserving accuracy across diverse data quality levels, even under challenging conditions. Further partition parameters, including peak wavelength and peak wave direction, are detailed in Appendix A and Appendix B.

4.3. Focused Analysis of Partition Classification Performance

The class distribution, depicted in Figure 3, demonstrates a consistent and largely homogeneous spread across all partitions, showing no significant bias towards any particular class. This indicates an overall well-balanced distribution.

While the class distribution remains balanced for both ascending and descending passes, a slight difference is noticeable in the medium and low classes when comparing the two modes. This discrepancy likely stems from combining the errors between these modes, which can shift the percentile limits compared to maintaining separate percentiles for each.

A closer inspection of Figure 4 reveals how the algorithm assigns quality classifications across the different wave partitions. The first partition—corresponding to the main swell—is consistently identified as the most significant and is predominantly labeled as "very good" or "good." This outcome aligns well with expectations, as the algorithm is designed to better identify the most relevant wave state information from the SAR data.

The classification of the second partition, typically associated with a secondary swell, exhibits a more uniform distribution across quality categories. This pattern is particularly evident in WV1 data, where it remains largely independent of satellite heading. In contrast, for WV2, the behavior of the second partition is slightly different, which could be attributed to the combined uncertainty introduced by merging the classification thresholds from the two acquisition modes.

The algorithm tends to classify these lower-ranked partitions as "low" or "poor" simply because they are less relevant from a geophysical perspective. This results in a clear distinction in partition quality, which can be especially useful for end users of SAR wave data when interpreting the most meaningful components of the sea state.

To further evaluate the performance of the quality classification, we examine the spatial distribution of quality classes ("Very Good," "Medium," and "Poor") across both WV1 and WV2 modes (Figure 5). These maps show spatial variations that appear broadly consistent between the two modes, suggesting that the classification patterns are more closely related to the underlying algorithmic process than to differences in acquisition mode.

Some of these spatial trends can be partially explained by the varying frequency of WV acquisitions across regions. For instance, the coverage near Europe is limited due to the use of alternative Sentinel-1 acquisition modes: Interferometric Wide Swath (IW) mode is primarily employed for land monitoring over the mainland, while Extra Wide Swath (EW) mode is more commonly used in areas such as the Azores. This results in a lower number of WV observations available for quality assessment in those zones.

A more detailed examination of the spatial distribution of quality classes highlights certain regional tendencies:

• Very Good: This class tends to be more frequent in mid-latitude. In contrast, it is less commonly observed in the northern Indian Ocean and in coastal regions, where environmental factors such as monsoon activity, coastal topography, and proximity to land can affect the quality of wave retrievals.
• Medium: This class is relatively evenly distributed across the globe, with increased presence in transitional zones near the equator and subpolar regions. These areas are characterized by more variable conditions that often result in intermediate-quality inversions.
• Poor: This class is more frequently observed at high latitudes in both hemispheres, where strong and variable wind speeds, along with complex atmospheric phenomena, can negatively impact SAR wave retrieval. Additionally, some chaotic offshore regions experience increased maritime traffic and environmental variability, which may introduce biases in the sea state inversion process. These factors contribute to lower quality classifications in these areas.

These spatial patterns provide end users of SAR wave data with clearer guidance on the expected quality of sea-state information, supporting a more informed use of the products in different geographic contexts.

5. Discussion

Accurately deriving ocean wave parameters from SAR data is inherently complex, influenced by both radar system limitations and environmental factors. A critical component in this process is the MTF, which defines how spatial frequencies in the SAR cross-spectrum are translated into the ocean wave spectrum. The MTF is particularly sensitive to wave motion along the radar’s range axis, whether moving toward or away from the radar. This motion introduces a Doppler effect that affects the radar backscatter signal and, consequently, the retrieval of the wave spectrum. These Doppler-induced shifts, combined with the inherent limitations of the MTF, add complexity in accurately interpreting SAR data for ocean wave analysis (Figure 6).

In addition to MTF-related challenges, environmental phenomena such as rainfall and other atmospheric or oceanic features further impact the NRCS, affecting wave retrieval. Rainfall, (Figure 7), can influence the NRCS by altering the surface roughness, potentially masking the geophysical signal. Oceanic fronts, characterized by sharp gradients in sea surface temperature and salinity, and atmospheric fronts (Figure 8), marked by strong near-surface wind variations, also introduce discontinuities in the NRCS, further complicating wave spectrum retrieval.

To go deeper into model explanation and interpretability, an illustrative example of Model interpretation using SHAP-based explanations (Figure 7 )is provided for estimating the errors in the three partition-integral parameters for the WV1 acquisition (see Figure 9, Figure 10, and Figure 11). The explanation for the WV2 acquisition is essentially the same and is therefore not shown. Each plot investiagtion holds two type of analysis one for the model expalanation : this explanation helps interpret the importance of each feature (Table 1) and how they affect the model’s predictions. A second plot based on features importance introduce the concept of feature clustering to help visualize redundancy among features. This clustering aims to group features that provide similar information to the model. This means that a model might be able to use either of two related features and achieve similar predictive performance. Traditional methods like correlation matrices can identify such relationships. In the context of SHAP’s feature clustering, "distance" between features is typically scaled between 0 and 1, where a distance of 0 means the features are perfectly redundant, and a distance of 1 means they are completely independent. In our case, we set the clustering cutoff at 0.5 : the bar plot will only show clusters (or groups) of features that have a clustering distance less than 0.5. This means only highly redundant features (those sharing more than 50% of their explanation power) will be visually grouped together. Less redundant features, even if they have some relationship, will be displayed individually.

The most impactful predictors for the $\left|\Delta \mathit{Hs}\right|$ prediction (as depicted in Figure 9 (a)) are the effective partition $\mathit{Hs}$ and the normalized variance of the SLC imagette. These features are physically meaningful, as they directly relate to wave energy and sea surface texture, and their influence is clearly reflected in the model’s behavior. Specifically, the model tends to associate higher effective $\mathit{Hs}$ values with lower predicted error, leading to a higher quality classification for those partitions. Conversely, lower $\mathit{Hs}$ values are typically linked to higher predicted errors and lower quality assessments. This indicates that the model prioritizes energetic sea states, while treating weaker partitions more conservatively.

A further insight comes from the feature clustering (Figure 9 (b)): statistical moments derived from the SLC imagette—namely skewness, kurtosis, and normalized variance—form a single cluster, indicating they collectively explain a significant portion (over 50%) of the model’s behavior. Another correlated group includes peak wave parameters such as peak direction and peak period. Most other features appear relatively independent in their contribution to predicting $\left|\Delta \mathit{Hs}\right|$, suggesting the model relies on a focused subset of physically meaningful indicators to make its predictions.

In the context of peak period error estimation $\left|\Delta \mathit{T}\right|$ (Figure 10 (a)), the predominant determinant of model fidelity is the energy ratio between the partition’s spectral peak and the maximum energy boundary. A diminished energy ratio corresponds to an escalation in relative error, underscoring its critical influence on the model’s uncertainty quantification. Conversely, elevated energy ratios are associated with reduced error magnitudes, indicative of enhanced precision in the peak period retrieval.

Consistent with the model’s objective, the peak period itself constitutes the secondary most salient feature impacting error prediction. Its influence manifests nonlinearly and is modulated by the swell propagation direction relative to the SAR azimuth, as well as hydrodynamic effects and tilt modulation intrinsic to the SAR wave retrieval mechanism.

In particular, while the partition index has a limited influence on the estimation of $\left|\Delta \mathit{Hs}\right|$, the model leverages the partition rank (i.e., the "p" feature) more substantially when predicting $\left|\Delta \mathit{T}\right|$. This suggests that the model attempts to refine its error estimation by accounting for the sequential ordering of swell systems, which aligns with the earlier observations regarding the distribution of quality classes across partition ranks. Such a strategy reflects the model’s recognition of the varying relevance and complexity of each partition in accurately retrieving peak period information.

Furthermore, we observe a similar behavior to $\left|\Delta \mathit{Hs}\right|$ in terms of feature redundancy (Figure 10 - (b)): the SLC imagette statistics are often redundant, leading the model to group them together. This redundancy extends to some other peak partition parameters computed in different geometries, which also appear to be redundant for the model.

The analysis of partition wave direction error $\left|\Delta \Phi \right|$ (Figure 11 (a)) reveals a pronounced sensitivity to the alignment between the dominant wave direction and the satellite heading. This confirms the fundamental influence of wave propagation dynamics relative to the SAR flight azimuth, which imposes intrinsic constraints on the inversion accuracy, as previously outlined. The elevated importance of the ambiguity factor—ranked third among predictors—further highlights its critical role in modulating directional uncertainty, distinguishing this error metric from others.

Consistent with patterns observed in the estimation of other wave parameters, the energy ratio once again proves to be a principal driver, ranking as the second most impactful feature. The peak wave direction $\Phi$ also contributes meaningfully, occupying the sixth position in the feature importance hierarchy, emphasizing its relevance in capturing directional variability.

Apart from the correlated statistical features derived from the SLC imagette (Figure 11 (b)), the model does not identify additional significant feature groupings influencing $\left|\Delta \Phi \right|$. This suggests that directional error arises from a more discrete set of physical and retrieval factors, underscoring the complex interplay between wave dynamics and SAR observation geometry. These insights pave the way for targeted improvements in SAR wave inversion methodologies, with the potential to enhance directional retrieval fidelity under challenging environmental conditions.

6. Conclusions

This study presents a novel approach to improving the reliability of SAR-based ocean wave retrieval through an advanced partition classification algorithm. By leveraging machine learning and explainability techniques, the method effectively addresses challenges that have long hindered the accuracy of wave spectrum interpretation—particularly those arising from system-induced distortions such as azimuth cutoff effects and the modulation transfer function.

What sets this work apart is not just the classification accuracy, but the intelligent selection of physically meaningful features that align with the behavior of ocean waves. This synergy between data-driven modeling and geophysical understanding enables a deeper insight into SAR performance under varying sea conditions. The model’s capacity to isolate the most trustworthy partitions creates opportunities for more targeted and confident use of SAR data in both scientific research and operational settings.

Beyond the technical contributions, the integration of interpretability adds a valuable layer of transparency. Rather than treating the algorithm as a black box, the use of SHAP-based analysis illuminates the factors influencing each decision, allowing users to validate outcomes and potentially refine the input data or retrieval strategies.

In a broader sense, this work contributes to a growing effort in Earth observation to make remote sensing tools not only more accurate but also more accountable and accessible. The proposed method offers a foundation upon which future SAR-based systems can be built—systems that are capable of adapting to complex marine environments while providing clear, interpretable results. As ocean monitoring becomes increasingly vital in the context of climate change and maritime activity, such tools will be essential for both research and real-world applications.