An Evolutionary Process-Embedded Spatiotemporal Interpolation Method for Marine Environmental Fields

1. Introduction

In the geographic world, there exists one kind of spatiotemporal phenomena that exhibit continuous and gradual changes, such as mesoscale ocean eddies, ocean fronts, wildfires, and the spread of epidemics. The highly spatiotemporally continuous and gradual characteristics of the marine dynamic environments make the study of ocean dynamics processes an essential part of digital twin ocean construction. It also plays a crucial role in promoting the development of smart ocean and addressing global climate change research [1,2,3,4]. Simulations, control, and optimization of ocean dynamic processes heavily rely on fine and accurate marine environmental fields. Marine environmental fields with a high temporal resolution play a key role in monitoring and forecasting ocean dynamics at regional weather scales [5,6,7,8]. However, due to the limitations of the temporal resolution, it is currently difficult to obtain the required grid data that supports a true representation of the ocean dynamics at the timescale of evolutionary, let alone meet the interpolation accuracy demands [6]. The discreteness of ocean spatiotemporal data severely restricts the development and application of modeling, continuous expression, and analysis in ocean dynamics processes [9]. Therefore, achieving high spatiotemporal continuity for marine environmental elements and ocean phenomenon data has become a critical issue [8,9,10]. The essence of transformation of marine environmental fields from discrete to continuous lies in the spatiotemporal reconstruction for missing data. Constructing appropriate geographic spatiotemporal interpolation methods is a vital tool for bridging the gap between the scale of observational data and the scale of geographical evolution, and it also provides the data foundation for storage, expression, simulation, and analysis in the digital twin of the oceans. Thus, there is an urgent need for effective spatiotemporal reconstruction methods to handle the creation of continuous spatiotemporal datasets for ocean environments.

Spatiotemporal interpolation refers to the process of constructing a mathematical model for the studied random variables to predict or estimate the unknown attribute values within a dataset [11]. Current spatiotemporal interpolation or reconstruction methods used in the field of marine information science can be broadly divided into three categories: physical model-based methods, geostatistics-based methods, and data-driven methods. Numerical physical model-based methods use a series of complex physical equations to describe the variation patterns of marine environmental factor fields. These equations are often intricate, resulting in high computational costs, and require substantial prior knowledge and assumptions, which makes them less adaptable to different ocean regions and time periods [12,13].

Statistical methods analyze the spatiotemporal characteristics of data and establish corresponding spatiotemporal dependencies for time, space, and spatiotemporal interpolation. Common methods include exponential smoothing, spline interpolation, inverse distance weighted interpolation, and kriging interpolation [14,15,16,17]. These methods are more flexible than physical models, but they require assumptions of spatiotemporal relationships and struggle to capture complex nonlinear relationships.

Data-driven methods use machine learning (ML) and deep learning (DL) techniques to capture nonlinear relationships in data through learning, establishing mappings that do not directly depend on physical models or explicit rules. This allows for the imputation, reconstruction, and prediction of missing data by discovering patterns and relationships from observational data and approximating target results. ML methods, such as random forests (RF) [18,19] and support vector machines (SVM) [20], directly discover patterns from historical data, though they are often less accurate because they cannot fully utilize large volumes of historical data to learn deep features of marine environmental factor data. DL methods, including convolutional neural networks (CNN) [21], long short-term memory networks (LSTM) [22,23,24,25], spatiotemporal graph neural networks (STGNN) [26,27,28,29], convolutional autoencoders (DINCAE) [30,31], and convolutional long short-term memory (ConvLSTM) networks [32,33,34,35,36], have been widely used for forecasting or interpolating marine environmental factor fields. These methods, which do not rely on theoretical or prior rules, have achieved state-of-the-art performance [37]. However, they still face challenges and limitations in effectively capturing both the temporal and spatial dependencies of marine environmental factor fields, as well as the variations in different patterns.

Despite notable advancements in marine data reconstruction, current models still exhibit two key limitations. First,few interpolation methods take into account the overall data trends, sample choices, and spatiotemporal dependencies [38]. Existing methods primarily focus on modeling from a single dimension (spatial or temporal), perspective, or model, failing to effectively integrate spatiotemporal dependencies. Whether using physical models, geostatistical methods, or AI techniques, each single model has its advantages and drawbacks in capturing spatiotemporal dependencies [39,40]. Particularly when marine environmental factors exhibit long-term trends and periodic variations, single models often mix multiple characteristic patterns, limiting their ability to precisely capture the spatiotemporal dependencies and dynamic evolutionary processes. As a result, it becomes difficult to fully extract spatial correlations and temporal dependencies, leading to inaccurate reflections of the complex evolution of marine environmental factors, which in turn reduces interpolation accuracy and limits the representation of true oceanic dynamical processes. Second, phenomena such as mesoscale eddies often display anomalies in marine environmental fluctuations, along with continuous gradual changes and moving wave characteristics in the spatiotemporal process. Existing marine spatiotemporal interpolation methods overlook the non-stationary, fluctuating, and moving spatiotemporal characteristics in the dynamic evolution of ocean phenomena. To overcome these two challenges, this paper introduces an Evolutionary Process-embedded Marine Spatiotemporal Interpolation Model (EPMSIM), based on the integration of multiple spatiotemporal modeling operators. It comprehensively considers the trend, periodic, and local spatiotemporal dynamic features of long time-series marine data, coupling geostatistical models with data-driven models.The main contributions of EPMSIM are threefold.

(1) This paper designs a spatiotemporal reconstruction model based on a convolutional bidirectional long short-term memory (ConvBiLSTM) network. By leveraging the bidirectional spatiotemporal feature extraction advantage of ConvBiLSTM, the model captures the high precision and stability of trend and seasonal features in long time series of marine environmental factors, addressing the interpolation issues related to periodic and trend characteristics in ocean spatiotemporal data.

(2) This paper proposes a process-based dynamic tracking interpolation method (PSDTIM). PSDTIM combines the maximum cross-correlation (MCC) from fluid motion estimation with spatiotemporal mixed inverse distance weighting (IDW) interpolation to capture the dynamic gradual-moving fluctuation characteristics in the evolutionary process at the local spatiotemporal scale, solving the interpolation problem for the evolutionary process component of ocean spatiotemporal data.

(3) A comparative evaluation of geostatistical and deep learning-based interpolation approaches is carried out to substantiate the enhanced performance of EPMSIM. Ablation experiments confirm the necessity and effectiveness of integrated coupling modeling and further validate the efficacy of PSDTIM when incorporating evolutionary process features for interpolation.

The rest of the paper is organized as follows: Section 2 reviews related works on marine spatiotemporal interpolation methods; Section 3 describes the model methodology in detail; Section 4 presents the experimental setup and the data; Section 5 presents the results analysis;Section 6 discusses the issue; Section 7 concludes the paper.

2. Related Work

In this section, we will describe three types of spatiotemporal interpolation methods in the field of marine information science: numerical physical models, statistical models, and data-driven models.

2.1. Numerical Physical Models

Methods based on numerical physical models rely on solving mathematical equations that describe physical processes. These models perform interpolation and reconstruction of spatiotemporal data using numerical simulation techniques or mathematical physical formulas. Relevant research methods include General Circulation Models (GCM) [12,13], Coupled Model Intercomparison Project (CMIP) models [41,42], and Regional Ocean Modeling Systems (ROMS) [43], among others. These physical models typically have high complexity and explainability of dynamic mechanisms. However, they face challenges such as high model complexity due to the interaction of multiple variables in physical processes, high computational resource consumption, regional and temporal applicability differences, uncertainty in parameter settings, and the need for substantial domain knowledge and prior assumptions [27,28,29].

2.2. Statistical Models

Statistical methods typically utilize spatiotemporal features of data to build appropriate statistical models that express spatiotemporal dependencies, enabling data imputation, reconstruction, and prediction. Typical methods include exponential smoothing, spline interpolation [14], inverse distance weighting interpolation [16], kriging interpolation [17], optimal interpolation [44], empirical orthogonal function data interpolation [45], and spatiotemporal non-integrated methods based on spatiotemporal mixing [11,38,39]. These methods are based on geostatistical principles or parameter models that require the establishment of complex spatiotemporal relationships or dependencies in advance. When interpolation is performed under conditions that do not meet these assumptions, the precision of the results is limited, and it is challenging to capture the nonlinear spatiotemporal relationships in complex marine environmental factor fields [38].

2.3. Data-Driven Models

Data-driven methods are models that make decisions and predictions by leveraging large datasets and recognizing underlying patterns. In contrast to traditional rule-based or theoretical models, data-driven models prioritize data itself over predefined assumptions or theoretical frameworks. The detailed related work of data-driven methods includes both basic machine learning (ML) and deep learning (DL) methods.

Traditional machine learning models mainly include Support Vector Machines (SVM) [20], Random Forests (RF) [18,19], and XGBoost models [46], which are regression models that primarily focus on extracting time features from marine environmental data. However, these methods are limited in capturing spatial features and fail to fully capture the deep nonlinear relationships of marine environmental factor sequence data.Deep learning-based methods mainly include the following three types of models.

Spatial deep learning models estimate values by utilizing the spatial dependencies of marine environmental factor fields in different regions. They primarily use CNN models to transform historical sea surface environmental factor data into image format for estimation [47,48], or use graph convolutional network (GCN) models to construct nodes and edges for estimating unknown locations [29].

Temporal deep learning models primarily rely on recurrent neural networks (RNN) and their variants to estimate values by capturing the temporal dependencies of marine environmental factor data. This category mainly includes RNN-based models [49], gated recurrent unit (GRU)-based models [50], and LSTM-based models [23,25]. Although these methods perform well in terms of estimation, they often focus only on temporal factors, neglecting the spatial relationships of marine environmental factor fields.

Spatiotemporal deep learning models introduce deep spatiotemporal data analysis methods, which can simultaneously extract both spatial and temporal features of marine environmental factor data. ConvLSTM was introduced specifically to address the limitation of fully connected LSTM (FC-LSTM) in preserving spatial correlations while processing spatiotemporal data [51], and it has been widely applied in ocean environment spatiotemporal data estimation [32,33,34,35,36]. Additionally, Hou et al. [52] developed the D2CL model, using dilated convolution networks and LSTM networks to learn spatial and temporal features. Li et al. [53]combined CNN with LSTM to propose a 3D CNN and LSTM-based SST prediction model. Zhao et al. [28] combined graph convolution with GRU to propose a spatiotemporal graph recurrent network to improve the performance of marine environmental factor data estimation. Zhou et al. [54] designed time and space transformer modules based on the Transformer model to extract both temporal and spatial correlations of input spatiotemporal sequences, improving the estimation performance of chlorophyll-a (Chla).

Although the deep learning models mentioned above have demonstrated strong performance in the estimation of marine environmental factor data, existing models still fail to fully address the problem of spatiotemporal data interpolation. First, most of these models are limited to a single model framework, making it difficult to simultaneously handle complex long-term spatiotemporal dependencies and fine-grained spatiotemporal pattern changes. Second, existing methods still have limitations in the application of data interpolation for ocean phenomena and other spatiotemporal ocean processes, failing to fully capture the dynamic changes in the complex evolutionary processes of ocean spatiotemporal data, especially in ocean phenomena. To address these shortcomings, this paper proposes a new spatiotemporal interpolation method, EPMSIM, which couples deep learning models with geostatistical models to fully integrate the advantages of multiple models. This approach captures complex long-term spatiotemporal dependencies, fine-grained pattern changes, and dynamic fluctuations in the evolutionary process, thereby improving the accuracy and stability of spatiotemporal process data interpolation in the ocean.

3. Methodology

The input data for EPMSIM is ocean time-series grid field data. The method includes four core steps: STL time-series decomposition, trend and seasonal feature interpolation based on ConvBiLSTM, process-based spatiotemporal dynamic tracking interpolation method (PSDTIM), and additive spatiotemporal coupling. The technical framework of EPMSIM is illustrated in Figure 1. The method mainly consists of the following components: First, time-series data is extracted from each spatial location of the ocean environment grid, and time-series decomposition is performed to obtain the seasonal component, trend component, and evolutionary process component (residual). The second part introduces the ConvBiLSTM model designed and built in this paper, which is used to integrate both spatial and temporal dimensions to capture long-term trends and seasonal patterns. The third part introduces the PSDTIM module, which considers the dynamic evolutionary process fluctuations at the local spatiotemporal scale. The fourth part involves the additive spatiotemporal coupling of the seasonal, trend, and evolutionary process components (residual) to obtain the final results. The detailed content is described in sections 3.1 to 3.5 below. Finally, using an experimental case of mesoscale ocean eddies based on sea surface temperature, the method's performance is validated by comparing it with multiple interpolation methods at daily and weekly scales.

3.1. Problem Statement

Spatiotemporal interpolation (reconstruction) of marine environmental factor fields typically divides the region of interest into a 2-D grid of rows (row) and columns (col) along latitude and longitude. Considering the influence of spatiotemporal sequence evolution on the variation of marine environmental factor fields, the modeling problem of marine environmental factor fields is described using the spatial expression (1):

$$\begin{array}{c}\mathrm{F}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}\left(\mathrm{E}\mathrm{l}\mathrm{e}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t}\right)=\left\{{\mathrm{f}}_{\mathrm{i},\mathrm{j}}\left(\mathrm{t}\right)\right|\mathrm{i}=\mathrm{0,1},\dots \mathrm{I}-1;\mathrm{j}=\mathrm{0,1},\dots \mathrm{J}-1,\mathrm{t}\in \mathrm{R}\}\end{array}$$

(1)

In this context, $\mathrm{F}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}\left(\mathrm{E}\mathrm{l}\mathrm{e}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t}\right)$denotes the marine environmental factor field, while ${\mathrm{f}}_{\mathrm{i},\mathrm{j}}\left(\mathrm{t}\right)$ describes the continuous variation function of the marine environmental factor at the (i, j) location within the grid. Additionally, I and J indicate the total number of rows and columns in the marine environmental factor grid, respectively.

At each specific time $\mathrm{t}$, the recorded marine environmental factor values of all grid regions form a matrix ${\mathrm{X}}_{\mathrm{t}}\in {\mathrm{R}}^{\mathrm{I}\times \mathrm{J}}$. The grid interpolation estimate for time t requires the data from previous time periods ${\mathrm{X}}_{\mathrm{t}-\mathrm{u}+1},{\mathrm{X}}_{\mathrm{t}-\mathrm{u}+2},...{\mathrm{X}}_{\mathrm{t}-1}$ and subsequent time periods ${\mathrm{X}}_{\mathrm{t}+1},...{\mathrm{X}}_{\mathrm{t}+\mathrm{u}-2},{\mathrm{X}}_{\mathrm{t}-\mathrm{u}-1}$.. Therefore, the interpolation estimate for ${\mathrm{X}}_{\mathrm{t}}$ is as shown in (2):

$$\begin{array}{c}{\mathrm{X}}_{\mathrm{t}}=F\left({\mathrm{X}}_{\mathrm{t}-\mathrm{u}+1},{\mathrm{X}}_{\mathrm{t}-\mathrm{u}+2},\dots {\mathrm{X}}_{\mathrm{t}-1},{\mathrm{X}}_{\mathrm{t}+1},\dots {\mathrm{X}}_{\mathrm{t}+\mathrm{u}-2},{\mathrm{X}}_{\mathrm{t}+\mathrm{u}-1}\right)\end{array}$$

(2)

where $\mathrm{F}$ represents a spatiotemporal interpolation model.

3.2. STL Time-Series Decomposition

The variation in marine environmental factor data is influenced by temporal factors, resulting in a strong correlation between the data. Time-series decomposition algorithms are effective at revealing these inherent correlations and identifying the regularities of marine environmental factor changes across different time scales.This can, in turn, assist in marine environmental factor modeling to achieve high accuracy [34,55].

The seasonal-trend decomposition using loess (STL) method, based on locally weighted regression (Loess) [56], is centered on using the loess technique for two iterative stages: an inner loop and an outer loop. First, it removes the trend component to extract the seasonal component. Then, it applies trend smoothing to the deseasonalized series. Finally, the time series is decomposed into three independent components: long-term trend, seasonal component, and evolutionary process component (residual). Ultimately, the time series is broken down into three distinct components: the long-term trend, the seasonal component, and evolutionary process component (residual). STL decomposition is an additive time-series decomposition method. Using STL decomposition, the extracted grid point time-series data can be broken down into three components, as described earlier, resulting in discrete subsequences. The decomposition expression is shown in (3):

$$\begin{array}{c}{\mathrm{T}}_{\mathrm{t}}={\mathrm{S}}_{\mathrm{t}}+{\mathrm{C}}_{\mathrm{t}}+{\mathrm{R}}_{\mathrm{t}}，\left(\mathrm{t}\right|0\le \mathrm{t}\le \left|\mathrm{T}\right|,\mathrm{t}\in \mathrm{Z})\end{array}$$

(3)

where ${\mathrm{T}}_{\mathrm{t}}$ represents the time-series data extracted from the grid point, ${\mathrm{S}}_{\mathrm{t}}$ is the seasonal component, ${\mathrm{C}}_{\mathrm{t}}$ is the trend component, and ${\mathrm{R}}_{\mathrm{t}}$ is the residual component; $\left|\mathrm{T}\right|$ is the length of the extracted time-series data.

Current STL decomposition techniques are currently focused primarily on time series forecasting of station-based data [57,58,59], with no applications in spatiotemporal interpolation of two-dimensional marine environmental fields. To effectively capture the changing characteristics of the marine environmental field, EPMSIM first extracts time-series data from each grid point within the ocean environmental field. It then applies the STL time series decomposition method to separate the data, extracting the long-term trend and periodic variation characteristics, as well as the evolving fluctuation (residual) changes, at each spatial location in the field.

3.3. ConvBiLSTM

3.3.1. Module Theory

The spatiotemporal interpolation problem of marine environmental factor time-series grids involves both temporal and spatial dimensions. Therefore, this paper integrates both dimensions and adopts ConvLSTM as the foundational building block for the model. Unlike fully connected LSTM, ConvLSTM determines the future state of a unit in the grid based on the input of that unit and the past state of its local neighbors. In this process, both spatial and temporal correlations are captured and utilized. The key difference between ConvLSTM and LSTM models is that, in the state-to-state and input-to-state transitions, convolution operations replace the matrix multiplication used in FC-LSTM. The working mechanism of ConvLSTM includes gates (input gate, forget gate, and output gate) and information flow. The following (4)-(8) represent the information state transmission equations in a single ConvLSTM unit.

$$\begin{array}{c}{\mathrm{i}}_{\mathrm{t}}=\sigma \left({\mathrm{W}}_{\mathrm{x}\mathrm{i}}\mathrm{*}{\mathrm{X}}_{\mathrm{t}}+{\mathrm{W}}_{\mathrm{h}\mathrm{i}}\mathrm{*}{\mathrm{h}}_{\mathrm{t}-1}+{\mathrm{W}}_{\mathrm{c}\mathrm{i}}\circ {\mathrm{c}}_{\mathrm{t}-1}+{\mathrm{b}}_{\mathrm{i}}\right)\end{array}$$

(4)

$$\begin{array}{c}{\mathrm{f}}_{\mathrm{t}}=\sigma \left({\mathrm{W}}_{\mathrm{x}\mathrm{f}}\mathrm{*}{\mathrm{X}}_{\mathrm{t}}+{\mathrm{W}}_{\mathrm{h}\mathrm{j}}\mathrm{*}{\mathrm{h}}_{\mathrm{t}-1}+{\mathrm{W}}_{\mathrm{c}\mathrm{f}}\circ {\mathrm{c}}_{\mathrm{t}-1}+{\mathrm{b}}_{\mathrm{f}}\right)\end{array}$$

(5)

$$\begin{array}{c}{\mathrm{c}}_{\mathrm{t}}={\mathrm{f}}_{\mathrm{t}}\circ {\mathrm{c}}_{\mathrm{t}-1}+{\mathrm{i}}_{\mathrm{t}}\circ \mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left({\mathrm{W}}_{\mathrm{x}\mathrm{c}}\mathrm{*}{\mathrm{X}}_{\mathrm{t}}+{\mathrm{W}}_{\mathrm{h}\mathrm{c}}\mathrm{*}{\mathrm{h}}_{\mathrm{t}-1}+{\mathrm{b}}_{\mathrm{c}}\right)\end{array}$$

(6)

$$\begin{array}{c}{\mathrm{o}}_{\mathrm{t}}=\sigma \left({\mathrm{W}}_{\mathrm{x}\mathrm{o}}\mathrm{*}{\mathrm{X}}_{\mathrm{t}}+{\mathrm{W}}_{\mathrm{h}\mathrm{o}}\mathrm{*}{\mathrm{h}}_{\mathrm{t}-1}+{\mathrm{W}}_{\mathrm{c}\mathrm{o}}\circ {\mathrm{c}}_{\mathrm{t}}+{\mathrm{b}}_{\mathrm{o}}\right)\end{array}$$

(7)

$$\begin{array}{c}{\mathrm{h}}_{\mathrm{t}}={\mathrm{o}}_{\mathrm{t}}\circ \mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left({\mathrm{c}}_{\mathrm{t}}\right)\end{array}$$

(8)

In this context, where ${\mathrm{i}}_{\mathrm{t}}$ represents the input gate, ${\mathrm{f}}_{\mathrm{t}}$ represents the forget gate, ${\mathrm{o}}_{\mathrm{t}}$ represents the output gate, ${\mathrm{c}}_{\mathrm{t}}$ represents the current state, ${\mathrm{c}}_{\mathrm{t}-1}$ represents the previous state, and ${\mathrm{h}}_{\mathrm{t}}$represents the final output,$\mathrm{W}$ represents the weight coefficients of the given gates, $\mathrm{b}$represents the corresponding bias coefficients, $\mathrm{*}$represents the convolution operator,and $\circ$represents the Hadamard product.$\mathsf{\sigma }$denotes the sigmoid activation function. Through convolution operations, spatial features of the data can be extracted, while LSTM captures the temporal variation features of the data.

3.3.2. Model Structure

The spatial and attribute changes between consecutive time points in the marine environmental field are intrinsically linked. The interpolation of unknown spatiotemporal points relies on the neighboring points in both the spatial domain and the temporal dimension. Therefore, in addition to the past spatiotemporal information, future spatiotemporal information is also significant for estimating the current interpolation point. However, a unidirectional ConvLSTM model can only consider the influence of past time points on future time points and does not account for the relationships with future time points. Therefore, this paper adopts a bidirectional ConvLSTM model structure. Unlike the unidirectional model, the ConvBiLSTM layer has two sets of ConvLSTM units, enabling bidirectional learning and optimization of the sequence grid data.

To better capture the spatiotemporal correlations between marine environmental factor fields, this paper constructs a spatiotemporal deep learning model by stacking ConvBiLSTM layers. The model consists of four layers: a feature extraction layer, two ConvBiLSTM layers, and a convolutional layer as the output layer. The input to the entire model is in the form of a five-dimensional tensor (batch size, time steps, width, height, features). The feature extraction layer takes the normalized time-series grid samples$\mathrm{X}=[{\mathrm{X}}_1,{\mathrm{X}}_2,...,{\mathrm{X}}_{\mathrm{t}}]$ as input, which are passed through two convolutional layers to extract the spatial features of the grid samples. The ConvBiLSTM layer uses bidirectional ConvLSTM to compute the time-series feature vector and combines the forward and backward output hidden states by concatenating the channel dimensions to propagate the output. The output layer uses convolutional mapping to output the model's normalized predictions, and the final results are obtained by reversing the normalization process, converting the output back to the original scale.

In this paper, the ReLU activation function is applied after the convolution operations in both the feature extraction and output layers to help the network learn more complex nonlinear features and further enhance the model's expressive power. Batch normalization is also applied to improve the model's generalization ability and to avoid problems such as gradient vanishing or explosion. After each ConvBiLSTM layer's output, a dropout regularization operation is applied to prevent overfitting. Finally, a convolutional layer is used for nonlinear mapping to output the results. Based on experimental experience, this paper sets the number of convolutional kernels in the two ConvBiLSTM layers to 8 and refers to the optimal performance parameters in [36], setting the size of all convolutional kernels to 3*3. The specific model composition is shown in Table 1, and the structural diagram of ConvBiLSTM is shown in Figure 2.

3.4. PSDTIM Module

Marine spatiotemporal processes refer to the continuous, gradual evolution of marine phenomena or entities throughout their entire life cycle. In the spatiotemporal evolution process, there are spatial and attribute differences, as well as intrinsic relationships, between consecutive time points. Therefore, interpolation methods must account for the temporal and spatial proximity correlations. Additionally, the dynamic process characteristics of marine phenomena—across their lifecycle (from emergence, expansion, stabilization, weakening, to extinction)—are continuous and gradual. This requires considering the motion and changes of targets within sequence images. Based on this, this paper proposes the PSDTIM module, which integrates spatiotemporal mixed IDW (Inverse Distance Weighting) to extract local spatiotemporal process features. The module continuously tracks sample points across consecutive time steps using fluid motion estimation methods. It then collects relevant sample points and computes sample values within surrounding spatial neighborhoods. Finally, a time series model with temporal inverse distance weighting is applied to assign weights to the tracked sample values, capturing the dynamic process characteristics of marine phenomena. The schematic diagram of the PSDTIM principle is shown in Figure 3.

The PSDTIM module is an algorithm developed based on fluid motion estimation and the spatiotemporal mixed IDW method. The fluid motion estimation method is an important approach for tracking the movement of targets in sequence images. Existing related research primarily focuses on cross-correlation methods and optical flow methods. Cross-correlation technology, also known as particle image velocimetry (PIV), is based on the principle of using maximum cross-correlation (MCC) to compare the correlation between two images, thereby determining the speed and relative displacement of fluid flow in the image. Additionally, the MCC method has been widely applied in the study of marine environmental factor field motion changes in the marine domain [60,61,62,63]. The core of the method is to extract the query window from the image pair and calculate the cross-correlation. By finding the direction of the matching window with the highest cross-correlation coefficient, the displacement vector within that window is obtained. The cross-correlation function operation is shown in (9). By solving for the highest correlation matching window in each neighboring image, an average displacement vector is obtained, forming a displacement vector field.

The displacement vectors generated can be used for motion estimation and tracking. In this study, the impact of spatiotemporal proximity is also considered, leading to the development of the PSDTIM algorithm. To improve computational speed, the study uses a Fast Fourier Transform (FFT)-based cross-correlation method to calculate the displacement field, which is then used in subsequent experimental calculations.

$$\begin{array}{c}{\mathrm{R}}_{\mathrm{A}\mathrm{B}}\left(\mathsf{\Delta }\mathrm{x},\mathsf{\Delta }\mathrm{y}\right)=\sum _{\mathrm{x},\mathrm{y}}{\mathrm{I}}_{\mathrm{A}}\left(\mathrm{x},\mathrm{y}\right)\times {\mathrm{I}}_{\mathrm{B}}\left(\mathrm{x}+\mathsf{\Delta }\mathrm{x},\mathrm{y}+\mathsf{\Delta }\mathrm{y}\right)\end{array}$$

(9)

where $(\mathrm{x},\mathrm{y})$ represents the pixel coordinates within the image window, and $(\mathsf{\Delta }\mathrm{x},\mathsf{\Delta }\mathrm{y})$ represents the displacement between the two image windows. ${\mathrm{I}}_{\mathrm{A}}$ and ${\mathrm{I}}_{\mathrm{B}}$ are the grayscale values of the corresponding pixel windows in the adjacent images A and B, respectively.

The specific process of the PSDTIM algorithm is as follows: ➀ Displacement Field Calculation: First, starting from the adjacent images before and after the target interpolation time, cross-correlation methods are applied between adjacent grids within the temporal neighborhood to estimate fluid motion and obtain the continuous displacement motion field. ➁ Sample Point Tracking: Starting from the target interpolation grid point, the displacement vector generated in step ➀ is used for continuous tracking within the temporal neighborhood time-series grid, thus obtaining sample points for both the forward and backward time points. ➂ Sample Value Calculation: Considering spatial correlation, for each tracked sample point at a given time, the spatial IDW formula is used to account for the influence of neighboring samples of the tracked sample point. The calculation of the sample value for a given time point is shown in (10)-(12). ➃ Sample Weighting: Generally speaking, entities that are closer to each other are more similar than those that are farther apart. Considering that the closer the temporal distance to the target interpolation time, the greater the influence, temporal IDW is used in the time series model to assign weights to the sample values obtained from the forward and backward time points in the temporal neighborhood, reflecting the temporal correlation. The formulas are shown in (13)-(14). ➄ Traversal of All Interpolation Grid Points: Repeat steps ➁, ➂, and ➃ for all grid points in the target interpolation grid until all interpolation points are completed.

$$\begin{array}{c}{\mathrm{f}}_{\mathrm{I}\mathrm{D}\mathrm{W}}\left(\mathrm{i}，\mathrm{j}\right)=0.5{\mathrm{f}}_{\mathrm{center}}+0.5{\displaystyle \frac{{\sum }_{\left({\mathrm{i}}_{\mathrm{k}}，{\mathrm{j}}_{\mathrm{k}}\right)\in \mathrm{N}}{\mathrm{w}}_{\mathrm{k}}\mathrm{f}\left({\mathrm{i}}_{\mathrm{k}},{\mathrm{j}}_{\mathrm{k}}\right)}{{\sum }_{\left({\mathrm{i}}_{\mathrm{k}}，{\mathrm{j}}_{\mathrm{k}}\right)\in \mathrm{N}}{\mathrm{w}}_{\mathrm{k}}}}\end{array}$$

(10)

$$\begin{array}{c}{\mathrm{w}}_{\mathrm{k}}={\displaystyle \frac{1}{{{\mathrm{d}}_{\mathrm{i},\mathrm{j}}^{{\mathrm{i}}_{\mathrm{k}},{\mathrm{j}}_{\mathrm{k}}}}^2}}\end{array}$$

(11)

$$\begin{array}{c}{\mathrm{d}}_{\mathrm{i},\mathrm{j}}^{{\mathrm{i}}_{\mathrm{k}},{\mathrm{j}}_{\mathrm{k}}}=\sqrt{(\mathrm{i}-{\mathrm{i}}_{\mathrm{k}}{)}^2+(\mathrm{j}-{\mathrm{j}}_{\mathrm{k}}{)}^2}\end{array}$$

(12)

where $\mathrm{f}(\mathrm{i},\mathrm{j})$ represents the value at position $(\mathrm{i},\mathrm{j})$, ${\mathrm{f}}_{\mathrm{center}}$ represents the value of the center point, $\mathrm{N}$ is the set of points within the neighborhood window, and the value of each neighboring point $({\mathrm{i}}_{\mathrm{k}},{\mathrm{j}}_{\mathrm{k}})$is $\mathrm{f}({\mathrm{i}}_{\mathrm{k}},{\mathrm{j}}_{\mathrm{k}})$, with the distance from the center point being ${\mathrm{d}}_{\mathrm{i},\mathrm{j}}^{{\mathrm{i}}_{\mathrm{k}},{\mathrm{j}}_{\mathrm{k}}}$.

$$\begin{array}{c}{\mathsf{\lambda }}_{\mathrm{i}}={\displaystyle \frac{\frac{1}{({\mathrm{d}}_{\mathrm{i}}{)}^2}}{{\sum }_{\mathrm{i}=\mathrm{s}}^{\mathrm{e}}\frac{1}{({\mathrm{d}}_{\mathrm{i}}{)}^2}}}\end{array}$$

(13)

$$\begin{array}{c}Z={\displaystyle \sum _{\mathrm{i}=\mathrm{s}}^{\mathrm{e}}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathrm{Z}}_{\mathrm{i}}\left({\mathrm{t}}_{\mathrm{i}}\right)\end{array}$$

(14)

where ${\mathrm{d}}_{\mathrm{i}}$represents the time interval between the $\mathrm{i}\_\mathrm{t}\mathrm{h}$ time point and the target time point; ${\mathsf{\lambda }}_{\mathrm{i}}$ is the sample weight at the $\mathrm{i}\_\mathrm{t}\mathrm{h}$ time point; ${\mathrm{Z}}_{\mathrm{i}}\left({\mathrm{t}}_{\mathrm{i}}\right)$ is the tracked sample value at the $\mathrm{i}\_\mathrm{t}\mathrm{h}$ time point; $\mathrm{s}$ is the sample start time; and e is the sample end time.See Error! Reference source not found. for details of PSDTIM.

Algorithm 1 PSDTIM (Inputs:Images,tTarget,SpatialRadius, TemporalRadius).

Result=InitializeGrid(Images.size(), Images[0].rows, Images[0].cols); tBefore = {t | t < tTarget, t ∈ Images.timePoints, tTarget - t ≤ TemporalRadius}; tAfter = {t | t > tTarget, t ∈ Images.timePoints, t - tTarget ≤ TemporalRadius}; // Displacement field calculation FOR i = 1 to tBefore.size() DO FOR j = 1 to tAfter.size() DO       DisplacementField[i][j]=CrossCorrelation(Images[tBefore[i]], Images[tAfter[j]]); END FOR END FOR FOR each GridPoint(x, y) in Result DO Samples = {}; // Sample point tracking FOR each time t in tBefore DO TrackPoint = (x, y); FOR timestep = tTarget-1 TO t STEP -1 DO Displacement = DisplacementField[timestep][timestep+1]; TrackPoint=TrackPoint - displacement[TrackPoint.x][TrackPoint.y]; END FOR // sample value calculation Neighbors=GetNeighbors(TrackPoint,Images[t], SpatialRadius); SampleValue =SpatialIDW(Neighbors, TrackPoint); Samples.Add(t, SampleValue); END FOR FOR each time t in tAfter DO TrackPoint = (x, y); FOR timestep = tTarget+1 TO t STEP 1 DO displacement=DisplacementField[timestep-1][timestep];       TrackPoint=TrackPoint+displacement[TrackPoint.x][TrackPoint.y]; END FOR // sample value calculation Neighbors = GetNeighbors(TrackPoint, Images[t], SpatialRadius); SampleValue = SpatialIDW(Neighbors,TrackPoint); Samples.Add(t, SampleValue); END FOR // Sample weighting FinalValue = TemporalIDW(Samples, tTarget); Result[x][y] = FinalValue; END FOR Return Result; End Algorithm

Where SpatialRadius is the spatial domain range, TemporalRadius is the temporal neighborhood range; CrossCorrelation is the cross-correlation calculation function, SpatialIDW is the spatial IDW function, and TemporalIDW is the temporal IDW function.

3.5. Additive Spatiotemporal Coupling

For each spatial location in the marine environmental field, the time series data can be decomposed into discrete subsequences. These subsequences can then be fitted into continuous functions. The continuous functions are coupled into one overall continuous function using additive operations. This is referred to as the time evolution function (15). Using this continuous function, interpolation can estimate missing or unknown time steps, transitioning from discrete to continuous. Based on this, the final EPMSIM interpolation result is obtained. This result is derived by comparing the interpolation outcomes from the additive spatiotemporal coupling trend ConvBiLSTM model, the periodic ConvBiLSTM model, and the PSDTIM interpolation model.

$$\begin{array}{c}{\mathrm{F}}_{\mathrm{row}\text{，}\mathrm{col}}\left(\mathrm{t}\right)={\mathrm{f}}_{\mathrm{s}}\left(\mathrm{t}\right)+{\mathrm{f}}_{\mathrm{c}}\left(\mathrm{t}\right)+{\mathrm{f}}_{\mathrm{r}}\left(\mathrm{t}\right),\left(\mathrm{t}\right|0\le \mathrm{t}\le \left|\mathrm{T}\right|)\end{array}$$

(15)

where ${\mathrm{F}}_{\mathrm{r}\mathrm{o}\mathrm{w},\mathrm{c}\mathrm{o}\mathrm{l}}\left(\mathrm{t}\right)$ represents the time evolution function at the grid point corresponding to the row i and column j,${\mathrm{f}}_{\mathrm{s}}\left(\mathrm{t}\right)$ is the seasonal function,${\mathrm{f}}_{\mathrm{c}}\left(\mathrm{t}\right)$ is the trend function, ${\mathrm{f}}_{\mathrm{r}}\left(\mathrm{t}\right)$ is the residual function, t is the continuous time, and $\left|\mathrm{T}\right|$ is the length of the extracted time-series data.

4. Experiments

This study conducts experiments using ocean mesoscale eddies based on sea surface temperature (SST) as a case study. The experiments focus on eddy and non-eddy regions at both daily and weekly time scales. Nine interpolation methods are selected for comparative analysis, including five traditional geostatistical spatiotemporal interpolation methods (linear interpolation, spline functions, exponential smoothing, spatiotemporal inverse distance weighting, and spatiotemporal kriging interpolation [64]), as well as four deep learning methods (BiLSTM [65], Conv3D [66], CNN-LSTM [67], and ConvLSTM [32]).

4.1. Dataset and Experimental Region Overview

The SST data used in this study comes from the Optimum Interpolation Sea Surface Temperature (OI-SST) MW_IR dataset product, which is a multisource sensor fusion product from the National Environmental Satellite, Data, and Information Service (NESDIS) of the U.S. National Oceanic and Atmospheric Administration (NOAA) [68]. This dataset covers the global range, with a spatial resolution of 0.09°, a temporal resolution of 1 day, and a time span from 2006 to 2020. The original time-series grid dataset is processed and merged to obtain daily-scale data, and weekly-scale data is derived by averaging the daily data. The daily dataset contains 5477 images, and the weekly dataset contains 784 images. The corresponding target interpolation time grids in both datasets are used as the ground truth for interpolation accuracy validation.

This study selects a typical mesoscale eddy case from 2006 for interpolation experimental validation. The eddy is a typical warm eddy that lasted from January 2 to February 21, 2006, with a lifecycle of about 50 days. During this period, the dynamic process of the eddy showed significant changes, with notable evolutionary process characteristics, making it suitable for an interpolation case study. The center of this eddy is located at 172-174.5°E and 28.5-31.5°N, and it reached its maximum area on January 29, making it a typical mesoscale eddy phenomenon. The interpolation case study region is selected to be between 172-175°E, 28-32°N. The study area is shown in Figure 4. The statistical table for the data at different scales for the entire year of 2006 in the study area is shown in Table 2. The time-series decomposition example for the daily-scale experimental data is shown in Figure 5.

4.2. Data Processing

In this study, the IDW interpolation method is used in the data preprocessing stage to fill in missing values in the time-series SST data, resulting in complete spatiotemporal image data. The input data is normalized using Z-score standardization, and the test data is normalized using the mean and standard deviation from the training data.

The dataset is divided using a sliding window method. Based on the backcasting window length proposed in [22], which is set to be at least four times the prediction length, this paper sets the window length to 4 and the step size to 1, with images from a certain time step used to predict the target time step's image. To ensure the independence of the training and testing sets, the test set is selected from the full year of 2006 and is not included in the training process. The remaining data is divided into an 80% training set and a 20% validation set. To prevent data leakage during time-series decomposition, the training and testing sets are processed independently during the time-series decomposition, ensuring they remain separate for subsequent interpolation experiments.

In terms of feature input, when performing interpolation for the trend component of the marine environmental factor field, the trend information obtained from the decomposition (trend) is used as the feature input for the trend model, with the feature count set to 1. For the seasonal component interpolation, considering the periodic position information, sine and cosine functions are used to handle the periodicity of time. The input features for the seasonal component model are (seasonal, sin_feature, cos_feature), so the feature count is set to 3. The formula for calculating the periodic position information is as follows.

$$\begin{array}{c}\mathrm{s}\mathrm{i}\mathrm{n}\_\mathrm{f}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{e}\left(\mathrm{\Phi }\right)=\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\Phi }\right)\end{array}$$

(16)

$$\begin{array}{c}\mathrm{c}\mathrm{o}\mathrm{s}\_\mathrm{f}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{e}\left(\mathrm{\Phi }\right)=\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\Phi }\right)\end{array}$$

(17)

$$\begin{array}{c}\mathrm{\Phi }={\displaystyle \frac{2\mathsf{\pi }\mathrm{*}\mathrm{i}\mathrm{d}\mathrm{x}}{\mathrm{T}}}\end{array}$$

(18)

where idx is the index of the current time step, and T is the length of the period, representing the number of time steps in a full cycle.

4.3. Experimental Parameter Settings

During the model training process, this paper uses Root Mean Squared Error (RMSE) as the loss function, with the number of epochs set to 200 and early stopping with a patience of 10. The AdamW optimizer is used for model training, with an initial learning rate of 0.001 and weight decay set to 0.01. The batch size is set to 32, and a dropout rate of 0.3 is used to prevent overfitting. Additionally, a cosine annealing learning rate scheduler is employed during training to dynamically adjust the learning rate. Specifically, the learning rate is adjusted every 10 epochs, with a minimum learning rate set to 1e-6. The entire model is trained in an end-to-end manner. Finally, the experiments are conducted on a PC with an Intel 12th i7 CPU, 32GB of memory, and an NVIDIA GeForce RTX 4060 GPU to train the deep learning networks involved in this study. To ensure fairness, all deep learning models used in the experiments are trained with the same training and validation datasets, as well as the same training-related hyperparameters (including loss function, training iterations, and batch size). In addition, for the PSDTIM model, the temporal neighborhood range is set to 4, and the spatial neighborhood range is set to 3*3.

4.4. Evaluation Metrics

To quantitatively evaluate the interpolation accuracy and morphological performance of the spatiotemporal interpolation models, this paper selects three evaluation metrics for accuracy assessment: root mean square error (RMSE), mean absolute error (MAE), and mean absolute percentage error (MAPE). The structural similarity index measure (SSIM) is chosen as the evaluation metric to assess the image morphological structural quality. The formulas are shown in (19)-(22).

$$\begin{array}{c}RMSE=\sqrt{{\displaystyle \frac{1}{\mathrm{n}}}{\mathsf{\Sigma }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}-{\widehat{\mathrm{y}}}_{\mathrm{i}}\right)}^2}\end{array}$$

(19)

$$\begin{array}{c}MAE={\displaystyle \frac{1}{\mathrm{n}}}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}}\left|{\mathrm{y}}_{\mathrm{i}}-{\widehat{\mathrm{y}}}_{\mathrm{i}}\right|\end{array}$$

(20)

$$\begin{array}{c}MAPE={\displaystyle \frac{100\mathrm{\%}}{\mathrm{n}}}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}}\left|{\displaystyle \frac{{\mathrm{y}}_{\mathrm{i}}-{\widehat{\mathrm{y}}}_{\mathrm{i}}}{{\mathrm{y}}_{\mathrm{i}}}}\right|\end{array}$$

(21)

$$\begin{array}{c}SSIM\left(\mathrm{x},\mathrm{y}\right)={\displaystyle \frac{\left[2{\mathsf{\mu }}_{\mathrm{x}}{\mathsf{\mu }}_{\mathrm{y}}+{\left({\mathrm{k}}_1\mathrm{L}\right)}^2\right]\left[2{\mathsf{\sigma }}_{\mathrm{x}\mathrm{y}}+{\left({\mathrm{k}}_2\mathrm{L}\right)}^2\right]}{\left[{\mathsf{\mu }}_{\mathrm{x}}^2+{\mathsf{\mu }}_{\mathrm{y}}^2+{\left({\mathrm{k}}_1\mathrm{L}\right)}^2\right]\left[{\mathsf{\sigma }}_{\mathrm{x}}^2+{\mathsf{\sigma }}_{\mathrm{y}}^2+{\left({\mathrm{k}}_2\mathrm{L}\right)}^2\right]}}\end{array}$$

(22)

where n is the total number of grid point samples, ${\mathrm{y}}_{\mathrm{i}}$ is the actual attribute value at the i-th time point, and ${\hat{\mathrm{y}}}_{\mathrm{i}}$ is the estimated value at the i-th target grid point. ${\mathsf{\mu }}_{\mathrm{x}}$ and ${\mathsf{\mu }}_{\mathrm{y}}$ are the means of the two images, ${\mathsf{\sigma }}_{\mathrm{x}}$ and ${\mathsf{\sigma }}_{\mathrm{y}}$ are the standard deviations of the two images, ${\mathsf{\sigma }}_{\mathrm{x}\mathrm{y}}$ is the covariance between the two images, L is the dynamic range of pixel values, and ${\mathrm{k}}_1,{\mathrm{k}}_2$ are two scalar constants, typically set as ${\mathrm{k}}_1$=0.01 and ${\mathrm{k}}_2$=0.03.

5. Discussion

5.1. Ablation Experiment on Integrated Coupling

To verify the advantages of time-series decomposition and integrated coupling, this study trained the ConvBiLSTM model using a dataset without time-series decomposition. The model parameters were set the same as those of the trend ConvBiLSTM model mentioned earlier. The model was applied to the daily-scale dataset, with the same interpolation times as in the previous daily-scale experiment, and compared with the EPMSIM method. The experimental results are shown in Figure 10 and Figure 11. The results indicate that, compared to the ConvBiLSTM model without time-series decomposition, EPMSIM achieved an average reduction of 9.57% in RMSE, 8.15% in MAE, and 8.24% in MAPE. This improvement is due to EPMSIM's integration of the advantages of different models, leveraging both data-driven and geostatistical models. By combining trend, seasonal, and evolutionary process features, EPMSIM achieves lower RMSE, MAE, and MAPE values, confirming the feasibility and effectiveness of coupling spatiotemporal deep learning models with geostatistical models through time-series decomposition.

5.2. Ablation Experiment on Evolutionary Process Component Reconstruction

To verify the advantages of PSDTIM in reconstructing the dynamic evolutionary process component, this study compares the following four scenarios at the same interpolation time points: using trend and seasonal ConvBiLSTM models to reconstruct the trend and seasonal components, while keeping these unchanged and using different combinations of Linear, SES, ConvLSTM, and PSDTIM to reconstruct the evolutionary process component. The purpose of these comparisons is to evaluate the performance of PSDTIM and its differences from other methods. As shown in Figure 12 and Figure 13, the experimental results indicate that, compared to the second-best combination, the PSDTIM combination achieves an average reduction of 23.96% in RMSE, 24.56% in MAE, and 24.42% in MAPE. Compared to the linear and SES methods, which reconstruct the evolutionary process component in the temporal dimension, PSDTIM integrates spatial information and dynamically moving evolutionary process information related to spatial correlations. In contrast to using ConvLSTM for reconstructing the evolutionary process component, PSDTIM takes into account the temporal relationships between past and future time points. Therefore, the interpolation results show lower RMSE, MAE, and MAPE values, confirming the feasibility and effectiveness of PSDTIM.

6. Conclusions

Facing the complex spatiotemporal dependencies and the spatiotemporal characteristics of dynamic fluctuations in the marine environmental fields, this paper addresses the limitations of single model-based spatiotemporal interpolation methods and proposes a spatiotemporal interpolation method by integrating the long-term trend and seasonal characteristics with the local spatiotemporal scale fluctuations, as well as the characteristics of the evolutionary process. As parallely coupling the deep learning model ConvBiLSTM with the geostatistical model PSDTIM, this method, named EPMSIM, outperforms the traditional geostatistical interpolation methods and deep learning-based baseline interpolation models. The main conclusions are as follows:

(1)

This paper proposes a spatiotemporal interpolation method that incorporates evolutionary process features. The core of the method is to use time-series decomposition to integrate and couple different interpolation models, constructing the trend and seasonal component reconstruction model (ConvBiLSTM) and the evolutionary process component reconstruction model (PSDTIM). These are then combined additively to reconstruct the spatiotemporal data. The method takes into account different modes in the spatiotemporal dimension, avoiding the limitations of using a single spatiotemporal interpolation model, thus improving interpolation accuracy.

(2)

A process-based spatiotemporal dynamic tracking interpolation method (PSDTIM) is designed. PSDTIM uses the spatiotemporal cross-correlation of residual fluctuations from the ocean spatiotemporal process decomposition to track sample points and perform weighted interpolation. This captures the dynamic, continuous, gradual changes and fluctuations of ocean phenomena, addressing the interpolation problem of the evolutionary process component in ocean phenomena. In the ablation experiment for evolutionary process component reconstruction, the PSDTIM combination outperforms the second-best combination, with average reductions of 23.96% in RMSE, 24.56% in MAE, and 24.42% in MAPE. This proves the effectiveness and feasibility of the algorithm.

(3)

In this study, four evaluation metrics were selected to evaluate and analyze the performance of five traditional geostatistical interpolation methods and four deep learning model methods at both the daily and weekly scales. The experimental results show that EPMSIM outperforms other interpolation methods in terms of RMSE, MAE, MAPE, and SSIM. Compared to the second-best model, RMSE, MSE, and MAPE were reduced by an average of 20.26%, 20.23%, and 19.68% in the daily-scale eddy experiment, 6.09%, 4.06%, and 5.11% in the daily-scale non-eddy experiment, 19.18%, 18.14%, and 18.90% in the weekly-scale eddy experiment, and 8.03%, 8.57%, and 7.98% in the weekly-scale non-eddy experiment. These results validate the effectiveness and feasibility of EPMSIM.

Our proposed method provides a reference for the reconstruction of geographic phenomena with continuous gradual process characteristics and has scalability for application to other geospatial datasets with continuous gradual features. This method could offer potential technical support and solutions for the reconstruction of other ocean-related satellite products, contributing to the study of marine ecological environments and ocean dynamic processes. Although the EPMSIM model demonstrates the best interpolation performance, the model still has some limitations. The EPMSIM model in this study has not yet implemented adaptive parameter adjustment and optimization. Future research will focus on improving the algorithm and optimizing the model. PSDTIM heavily relies on the accuracy of fluid motion estimation methods. Future work will consider using physical motion models for fluid motion estimation to improve the method. During the model construction, only single-element time-series environmental data was used, without considering the impact of external environmental factors on geographic phenomenon data interpolation. Future work will integrate environmental factors and covariates to enhance the spatiotemporal interpolation model.