In this work, high-precision in situ observation profiles acquired from the UK Met Office Hadley Centre’s EN4 temperature and salinity profiles dataset version 4.2.1 (subsequently referred to as EN4-profiles) were used as model training labels. The World Ocean Database (WOD), the Arctic Synoptic Basin-wide Observation (ASBO), the Global Temperature and Salinity Profile Program (GTSPP), and the Argo Global Data Assembly Centers (GDACs) provide the fundamental observations in EN4-profiles (Good et al., 2013). EN4-profiles are widely used to evaluate model simulations as “ ground truth” (Kumar et al., 2017). We selected high-quality temperature profiles through quality flags.

EN4-profiles have a discontinuous and irregular spatial and temporal distribution and need to be interpolated into the previously mentioned target grid. First, we utilized linear interpolation to interpolate EN4-profiles to 23 standard levels. The processed profiles were then interpolated level by level onto the previously described 0.25° horizontal grid. Because the observation profiles are extremely sparse, a spatial-temporal weighted interpolation method (Zeng and Levy, 1995) was used to increase the number of samples and improve the interpolated data accuracy. For each horizontal objective grid to be interpolated, a spatial-temporal domain with a spatial radius R_s and a temporal radius R_t are specified as the interpolation neighborhood centered on the target grid. The objective grid will be null if there are no observation profiles in the interpolation neighborhood. The monthly temperature T_obj.i of the level i of the objective grid is computed as

where N represents the total number of observations in the neighborhood, T_k represents the k-th temperature observation in the neighborhood, and w_k represents the interpolation weight of T_k. The w_k is calculated as

Where x_k, y_k, t_k represent the longitude, latitude, and time corresponding to T_k, x₀, y₀ and t₀ represent the longitude, latitude, and time corresponding to T_obj,i, respectively. In this work, R_s is 0.48°, R_t is 15 days, and t₀ is the 16th day of each month.

To yield sea surface information for ocean data fusion, multi-source ocean satellite observation and analysis data were used as model training input. The surface variables include sea surface temperature (SST), sea level anomaly (SLA), and sea surface wind (SSW). Three satellite SST analysis datasets were collected, including NOAA’s Optimum Interpolation Sea Surface Temperature (OISST, version 2) (Reynolds et al., 2007), the Extended Reconstructed Sea Surface Temperature (ERSST, version 5) (Huang et al., 2017), and the Hadley Centre Global Sea Ice and Sea Surface Temperature (HadISST) (Rayner et al., 2003). The OISST data is daily with a horizontal resolution of 0.25°, the ERSST data is monthly with a horizontal resolution of 2°, and the HadiSST data is monthly with a horizontal resolution of 1°. The satellite SLA is a daily Aviso-SLA (version 4.0) dataset from Copernicus Marine Environment Monitoring Service (CMEMS) with a horizontal resolution of 0.25°. NASA’s monthly Cross-Calibrated Multi-Platform (CCMP, Version 2) wind data with a horizontal resolution of 0.25° is provided by the satellite SSW (Atlas et al., 2011). All sea surface data were collected between 1994 and 2017.

To generate subsurface information for ocean data fusion, numerical model data, objective analysis data, and reanalysis data were used. The addition of numerical model data and reanalysis data increased the physical rationality of the 3D ODF-ST dataset. Monthly historical simulation data (r1i1p1f1) from NCAR’s CESM2 Earth system model (Danabasoglu et al., 2020) and monthly historical simulation data (r1i1p1f1) from CMA’s BCC-CSM2-HR climate system model (Wu et al., 2021) were included in numerical model data. The Hadley Center’s EN4 monthly objective analysis data (version 4.2.1, subsequently referred to as EN4-analysis) is used for the objective analysis, which has a horizontal resolution of 1°. (Good et al., 2013). SODA (version 3) monthly ocean reanalysis data from the University of Maryland (Carton et al., 2018), ECCO (version 4) monthly reanalysis dataset with a horizontal resolution of 0.5° from NASA (Forget et al., 2015), and HYCOM daily reanalysis data with a horizontal resolution of 1/12° from the US Naval Research Laboratory are the reanalysis datasets used (Chassignet et al., 2007). All of these 3D gridded datasets were collected between 1994 and 2017.

To produce the monthly average dataset for sea surface data, daily OISST and AVSIO SLA analysis data were averaged. Using bilinear interpolation, all sea surface data were uniformly interpolated to the previously described 0.25° horizontal grid. The daily HYCOM reanalysis data for 3D gridded data were averaged to generate a monthly average dataset. All 3D gridded data were linearly interpolated to vertical standard levels before being uniformly interpolated to the previously defined horizontal grid with 0.25° resolution using a bilinear interpolation method.

The ODF-Net, shown in Figure 2A , is a variant of the 1D-ResNet model. There are three distinct sections to this design. To accomplish the task of extracting shallow features from multi-source ocean data, block A contains a 1D convolutional layer that uses a 3×1 kernel in addition to two GAMs. Blocks B and C, depicted in Figures 2B, C , respectively, are composed mostly of 1D convolutional layers with a 3×1 kernel, dropout layers (dropout probabilities of 0.5), and skip-connection structures to harvest the deep features. Block D is a decoding block that fully integrates sea temperature data from many sources through the use of a combination of shallow and deep features.

In the ODF-Net, we choose to use the more advanced Adaptively Parametric Rectifier Linear Unit (APReLU) activation function rather than the more common ReLU activation function used in the 1D-ResNet. When the original features are less than zero, APReLU runs each sample through a small fully connected network to produce matching weights, which are then used as coefficients of the original features to provide a more flexible method of nonlinear transformation (Zhao et al., 2020).

Improving the model’s interpretability assists with both understanding the deep learning model’s complex decision-making foundation and guaranteeing the model’s reliability (Xu et al., 2021). Giving neural networks an attention mechanism improves their ability to learn by focusing on the relevant important feature while discarding the rest. In order to emphasize the interaction of multiple sources of information at different depth levels and to enable the model to capture important features in both dimensions, we modify the GAM (Liu et al., 2021) and redesign the sub-modules. Figure 3A depicts the GAM’s overall attention mechanism process, which sequentially combines channel attention and spatial attention. Given input features F₁ the intermediate state F₂ and output F₃ are defined as

where M_c and M_s are the channel attention and spatial attention weights, respectively, indicating the model’s degree of attention to distinct channels and depth levels of input features. ⊗ denotes the multiplication operation by element. For the first GAM, in particular, channel attention reflects the importance of different sources, while spatial attention reflects the model’s concentration on diverse physical depth levels.

A major part of the channel attention submodule is shown in Figure 3B . When extracting information in two dimensions, 2D permutation is employed, and then a two-layer neuron network is used to magnify the dependence between channels and depth levels across dimensions. Each channel’s weights are then calculated using the sigmoid function. Figure 3C depicts the spatial attention submodule, which uses two convolutional layers to aggregate information from different depths and a sigmoid function to determine the weights of each depth.

In order to integrate prior knowledge such as the vertical structure of sea temperature into the model and reduce the spurious high-frequency signal of fusion data in the vertical direction, the vertical gradient and integral of the sea temperature profile were added to the loss function, which is calculated as

where Loss_{ST_Grad} is the sea temperature gradient constraint, Loss_Cumsum is the sea temperature profile integral constraint. α and β are hyperparameters, which were determined as 0.004 and 0.02 respectively through comparison experiments. The Loss_{ST_Grad} and Loss_Cumsum are defined as

where ST_Grad and Cumsum represent the vertical gradient and the vertical integral of the sea temperature profile, respectively. T_i denotes the sea temperature value of the i-th level and Z_i denotes the vertical depth of the i-th level.

Meanwhile, to determine whether the loss function with gradient and integral constraints may improve the model’s fusion performance, a comparative experiment was run with the loss function Loss₂=RMSE, and the experimental results are provided in Table 1 . R² is the coefficient of determination, a value ranging from 0 to 1 that indicates how effectively a statistical model predicts an outcome. The closer a model’s R² is to 1, the better it is at making predictions. Equation 8 gives the calculation of R².

where f_i and y_i represent the i-th prediction and label, respectively. y ¯ denotes the averaged value of all labels. The model with Loss₁ has nearly the same root mean square error (RMSE) as the model with Loss₂ but the sea temperature profile gradient error (STPGE) is reduced by 4%, indicating that the addition of physics-based prior knowledge constraints in the loss function has an enhancement effect on the vertical structure of the fusion sea temperature.

Four sets of experiments (Exp0–Exp3) were designed to validate the favorable impacts of GAM and APReLU on the fusion model. Exp0 is the baseline, which does not include GAM or APReLU. Exp1 includes GAM, and Exp2 includes APReLU. Exp3 includes both GAM and APReLU, i.e., the ODF-Net. In addition, we designed Exp4 to do the same fusion task using Transformer (Vaswani et al., 2017), a state-of-the-art model for sequence-to-sequence learning, in order to validate the ODF-Net’s performance in comparison to other models. The transformer’s hyperparameters were tuned, and the encode dimension, number of attention heads, number of identical layers, query vector length, and key vector length were all set to 128, 8, 6, 16, and 16, respectively.

Table 2 shows a comparison of test results from the five sets of experiments. Comparing Exp1 and Exp0, GAM reduces RMSE by 2.87% and STPGE by 4.28%; comparing Exp2 and Exp0, APReLU reduces RMSE by 4.33% and STPGE by 5.32%; comparing Exp3 and Exp0, GAM and APReLU reduce RMSE by 6.98% and STPGE by 6.86%. As a result, GAM and APReLU considerably increase model performance. The comparison of Exp4 and Exp3 shows that the ODF-Net outperforms the Transformer in all metrics, including RMSE, STPGE, and R². The quantities of trainable parameters are listed in the third column of Table 2 with M representing 10⁶, the higher the value, the more complicated the model. The amount of trainable parameters in the ODF-Net is only about 60% of the Transformer, demonstrating that the ODF-Net we developed is lightweight and high-performance.

To verify the accuracy of the fusion model, we compared the RMSE, STPGE, and R² of all eight data sources, including ODF-Net predictions, six sets of fused data, and their ensemble average data (EAD), with in situ observations. Table 3 summarizes the accuracy of different datasets over the test set using the metrics mentioned above. As Table 3 suggests, ODF-Net predictions are optimal on all evaluation metrics. ODF-Net predictions have an average RMSE of 0.74°C, while other data sources have RMSEs greater than 0.9°C. ODF-Net predictions have an average STPGE of 0.042°C/m, while other data sources have average STPGEs greater than 0.046°C/m. The ODF-Net improves both the accuracy and the vertical structure of the fusion sea temperature.

The spatial distribution of RMSE is depicted in Figures 4A, B . ODF-Net predictions are more accurate than EAD’s in almost all regions, especially in areas with large gradients such as the Gulf Stream, the Kuroshio Extension, and the West Wind Drift, where the improvement is noticeable. Figure 4C shows the spatial distribution of the percentage improvement in R² of ODF-Net predictions with observation profiles versus R² of EAD predictions with observation profiles. Statistical examination of the test set reveals that R² of ODF-Net predictions is greater than EAD for 78.51% of the profiles. The improvement is more significant in areas with large gradients, such as the boundary current regions, which are similar to the spatial distribution of RMSE, suggesting that the deep learning model learns more correct information from multiple datasets.

The distribution and variation of sea temperatures in different ocean areas and depth levels exhibit distinct characteristics due to the effects of several factors, such as solar radiation, land-sea distribution, ocean currents, and monsoons (Chen et al., 2002; Li et al., 2020). To examine the fusion effect, we evaluated the accuracy of ODF-Net predictions at different depth levels by region. The global ocean was divided into five oceans, which are the Pacific Ocean, the Atlantic Ocean, the Indian Ocean, the Arctic Ocean, and the Southern Ocean. ODF-Net performance in the five oceans and global regions is then discussed.

Here, we examined the time average RMSEs of eight data sources at different depth levels in the five oceans and global regions, including CESM2 (green line), BCC-CSM2-HR (red line), ECCO (purple line), SODA (pink line), EN4-analysis (orange line), HYCOM (brown line), EAD (gray line), and ODF-Net predictions (blue line) ( Figure 5 ). The time average RMSEs in the vertical direction of eight data sources varied in each of the five oceans, but the RMSE of ODF-Net predictions at different depth levels is significantly lower than that of other data sets in each ocean as well as the global region. This indicates that the proposed ODF-Net performs better on a global scale. The accuracy improvement in ODF-Net predictions above the thermocline is greater than in other deeper layers, particularly at the thermocline with the largest vertical gradient. This might be attributed to the large number of thermocline observations, which enables the ODF-Net to learn the bias between other data sources and observations.

The Taylor diagram incorporates numerous assessment measures that are commonly used to evaluate model performance, including the correlation coefficient (COEF), root mean square error (RMSE), and standard deviation (STD) (Taylor, 2001). Taylor diagrams were utilized to more thoroughly and objectively analyze the statistical connections between various data points and observations in this work. Taylor diagrams of sea temperatures in the five oceans and global regions based on nine data sources, including projections (black dots) and observations (red dots) from ODF-Net are shown in Figure 6 . Figure 6 shows that the results from several datasets vary widely, while the predictions generated by ODF-Net consistently and comprehensively outperform those generated by any other dataset. ODF-Net predictions are unbiased with a high level of correlation. The ODF-Net also has different outcomes depending on location. The performance of the ODF-Net in the Pacific, Atlantic, and Indian Oceans was comparable to that of the global region ( Figures 6A–C ). Figures 6D, E illustrate that the performance of the ODF-Net model in the Northern and Southern Oceans still needs to be improved. The main reason is that the performance of ODF-Net is highly dependent on high-quality observations, which are much less available in the Arctic and Southern Oceans than in other regions.

Attention weights, as an intermediate output of the network model, can be used as a convenient tool to explain model decisions. Many studies have discussed the model interpretation ability of the attention weight distribution for neural network models based on attention mechanisms (Pruthi et al., 2019; Serrano and Smith, 2019; Wiegreffe and Pinter, 2019). In this study, we utilized attention weight distribution to analyze the contribution of multi-source ocean data as well as spatiotemporal information to the ODF-Net fusion process. Since multi-source data were sent directly to the first GAM in the ODF-Net, the attention weights could represent the contributions of the original ocean data as well as of the spatiotemporal information. We obtained the channel attention weight M_c and the spatial attention weight M_s of the first GAM, and since the GAM used a combination of channel and spatial attention serially, the global attention weight M_global is defined as:

Figure 7 shows a heat map of the global attention weight M_global. The top six rows illustrate the contribution of each of the 3D gridded datasets to the final fusion task, while each column represents the contribution of ocean data at various depth levels. As shown in Figure 7 , the top three 3D gridded data contributors, in order, are EN4-analysis, SODA, and HYCOM, whereas ECCO, CESM2, and BCC-CMS2-HR contribute relatively little. The distribution of attention weights is rather reasonable, as the error between EN4-analysis and observed EN4-profiles is the smallest on the test set, followed by SODA and HYCOM, whereas the average errors are larger for ECCO, CESM2, and BCC-CMS2-HR. This implies that the ODF-Net has given more attention to high-precision data. In the spatial dimension, the attention weights of the shallow levels above 100 m are greater than those of the deeper levels below 100 m for EN4-analysis and SODA, the two datasets that received the most attention, and for other datasets, the weights of the deep levels are greater than those of the shallow levels. This is due to the fact that sea temperatures are more stable at deeper levels, and 3D gridded datasets at deeper levels are more accurate than those at shallow levels. Due to the extremely high correlation of sea temperature variations with latitude, latitude has the greatest weight in the last channel. In contrast, the attention weights of temporal information and sea surface data are not significantly different. The analysis of attention weights shows that the ODF-Net pays more attention to the data sources that are more accurate. At the same time, information from different depth levels of the same data source makes different contributions to the fusion process.

Since the purpose of this study is to accomplish accuracy, high resolution, and spatiotemporal continuous intelligent fusion of multi-dimensional, multi-source heterogeneous ocean data, we used the ODF-Net to generate a 3D sea temperature fusion dataset (ODF-ST dataset). ODF-ST dataset spans the years 1994 to 2017, with a spatial extent of the global ocean (180°W–180°E, 60°S–65°N), a monthly temporal resolution, and a spatial resolution of 0.25°. To assess the spatial rationality of the ODF-ST dataset, we compared the fusion SST with that of OISST, the spatial distribution of fusion sea temperature profiles, and sea temperatures at different levels with WOA18 (Boyer et al., 2018). To assess the temporal rationality of the ODF-ST dataset, we compared the fusion sea temperature profiles with Tropical Atmosphere Ocean Array (TAO) monthly observation profiles. Finally, we evaluated the ENSO index time series calculated with fusion sea temperature.

The global climatic SSTs (averaged from 1994 to 2017) of the ODF-ST dataset and the OISST are shown in Figure 8 , where both SSTs have a “low-high-low” distribution from north to south, with a notable high-value area in the mid-western Pacific Ocean and a center SST of nearly 29°C. The transition of 25°C isotherms in the east-central Pacific Ocean (red frame area) is similar. The differences in SSTs are in the range of 0.5°C in most regions, indicating that the distribution of ODF-ST SST is acceptable and trustworthy. ODF-ST SST is significantly higher than OISST in Hudson Bay, the Mediterranean Sea, the southwestern coast of Africa, and the western Okhotsk Sea, but markedly lower than OISST in the northwestern and southwestern Atlantic Ocean and the southeastern Okhotsk Sea, possibly due to more complex changes in nearshore currents and sparse observations.

The distribution of global climatic sea temperatures (averaged from 1995 to 2017) between the ODF-ST dataset and WOA18 at 500 m ( Figures 9A–C ) and 900 m ( Figures 9D–F ) shows that the ODF-ST sea temperature is very similar to WOA18. From north to south, both 500 m sea temperatures exhibit a “low-high-low-high-low” distribution, with noticeable high-value areas in the northwest Atlantic Ocean, the Mediterranean Sea, the southwest Indian Ocean, and the northwest Pacific Ocean, with the center sea temperature of the high-value area in the northwest Atlantic Ocean being around 17°C. Both 900 m sea temperatures have obvious high-value areas in the mid-eastern Atlantic Ocean, the Mediterranean Sea, and the Gulf of Aden, with the center sea temperature of the high-value area in the Mediterranean Sea being around 15°C. The isotherm trends are also strikingly similar, with homologous transitions of 500 m sea temperatures in the northern Atlantic Ocean and the southern Indian Ocean(red frame area). The differences between 500 m and 900 m sea temperatures of the ODF-ST dataset and WOA18 in most regions are lower than 0.25°C. Differences in 500 m sea temperatures have noticeable positive and negative oscillations in the Gulf Stream, the Kuroshio Extension, the North Pacific Current, and the West Wind Drift. Differences in 900 m sea temperatures have noticeable positive and negative oscillations in the Nansha Islands and the West Wind Drift.

To evaluate the rationality of ODF-ST sea temperature in the vertical direction, we averaged sea temperatures from the ODF-ST dataset and WOA18 in the longitudinal and latitudinal directions and then compared their sea temperature profiles along the latitudinal ( Figures 10A–C ) and longitudinal (Figures 10D–F) directions. The profile along the latitudinal direction demonstrates that the ODF-ST dataset and WOA18 both have a high sea temperature value of 28°C in-depth levels over 100m near 5°N. The isotherms exhibit evident grooves in the northern and southern hemispheres’ mid-latitudes. The difference in sea temperature is mostly less than 0.25°C. ODF-ST sea temperature is about 0.5°C higher than WOA18 in depth over 80 m in the Northern Hemisphere.

The profile along the longitudinal direction shows that the ODF-ST dataset’s isotherm change is essentially consistent with that of WOA18. Both have three regions with strong gradient variations located at 60°W–80°W, 0°–50°E, and 110°E–150°E, respectively. The differences in sea temperature are also higher in these regions, while differences in other regions are less than 0.25°C. The highly consistent ODF-ST sea temperature with WOA18 in the vertical direction implies that the spatial distribution of ODF-ST sea temperature is reasonable.

We used monthly TAO observation profiles from 1994 to 2017 to conduct a comparative analysis of temporal correlation in order to assess the temporal rationality of the ODF-ST sea temperature. Figure 11 depicts, for each observation site, the spatial distribution of the temporal correlation between the ODF-ST sea temperature and TAO observation profiles.The Pearson correlation coefficient has a range of 0.9755 to 0.9995 and a mean value of 0.996.Therefore, in the sea area where the TAO array is deployed, the temporal distribution of ODF-ST sea temperature is reasonable.

To verify the ODF-ST dataset’s ability to capture the ENSO signal, Figure 12 illustrates the ENSO index time series in the Nino3.4 region (5°S–5°N, 170°W–120°W) of multiple sources from January 1994 to December 2017. Variations of the ENSO index in the ODF-ST dataset are generally consistent with those of ERSST and HadiSST and can reflect the significant El Niño years (1995, 1998, 2003, 2007, 2010, 2015) and La Niña years (1999, 2000, 2008, 2011, 2012, 2017).

We compared the Pearson correlation coefficient and RMSE of the anomalies in the Nino3.4 region of average sea temperatures above 100 m from different sources with the ENSO index provided by the U.S. Climate Prediction Center (NOAA/CPC) and discovered that the Pearson correlation coefficients of the ODF-ST dataset, HYCOM, ECCO, and SODA were all above 0.9 or higher, and the RMSEs of them were all below 0.5°C, indicating that the ODF-ST dataset’s sea temperature field could fairly reflect the ENSO signal.