Oceanic Mesoscale Eddy Fitting Using Legendre Polynomial Surface Fitting Model Based on Along-Track Sea Level Anomaly Data

Abstract

Exploring the spatial distribution of sea surface height involves two primary methodologies: utilizing gridded reanalysis data post-secondary processing or conducting direct fitting along-track data. While processing gridded reanalysis data may entail information loss, existing direct fitting methods have limitations. Therefore, there is a pressing need for novel direct fitting approaches to enhance efficiency and accuracy in sea surface height fitting. This study demonstrates the viability of Legendre polynomial surface fitting, benchmarked against bicubic quasi-uniform B-spline surface fitting, which has been proven to be a well-established direct fitting method. Despite slightly superior accuracy exhibited by bicubic quasi-uniform B-spline surface fitting under identical order combinations, Legendre polynomial surface fitting offers a simpler structure and enhanced controllability. However, it is pertinent to note that significant expansion of the spatial scope of fitting often results in decreased fitting efficacy. To address this, the current research achieves the precise fitting of sea surface height across expansive spatial ranges through a regional stitching methodology.

1. Introduction

Sea Level Anomaly (SLA) denotes the deviation of sea level from its long-term average and is a crucial metric for analyzing ocean dynamic processes and climate change. SLA data are essential for understanding ocean circulation patterns, which significantly influence the global climate system. By examining SLA variations, researchers can track the oscillations of ocean currents and their role in global heat distribution. As global temperatures rise, the resultant melting of glaciers and ice shelves contributes to sea level rise, with SLA data providing direct evidence of these phenomena. Long-term SLA records are invaluable for identifying trends in sea level rise, which is vital for assessing and mitigating the long-term effects of climate change [1].

The computation of SLA involves several steps. Initially, the actual sea surface height (SSH) is measured using satellite altimetry. Subsequently, the mean sea surface height (MSSH) is derived from long-term SSH data. SLA is then calculated as the difference between the SSH and MSSH. When integrating data from multiple satellite missions, it is imperative to perform data fusion and cross-calibration to ensure consistency.

SLA observations predominantly utilize satellite altimetry technology. Satellite altimeters ascertain sea level height by emitting microwave pulses toward the Earth’s surface and measuring the time delay of the reflected signals. Additionally, tide gauges installed along coastlines provide precise local sea level data, which complement satellite altimeter measurements. By synthesizing these observational techniques and employing advanced data processing methods, researchers can generate high-precision and high-resolution SLA datasets, thereby furnishing a reliable foundation for oceanographic and climate research [2].

Recent studies have demonstrated that the application of composite sea state bias (SSB) correction can significantly enhance data accuracy by mitigating SSB and correlation errors within the 1 Hz frequency range [3]. Methods for identifying mesoscale eddies using SLA data include closed SLA contour lines, the Okubo-Weiss number, the Winding Angle method, and flow vector techniques [4]. The utilization of SLA data, coupled with analytical methods such as linear regression, Winters exponential smoothing, and empirical mode decomposition, enables the elucidation of spatiotemporal characteristics pertaining to sea level changes [5].

Mesoscale eddies are a prominent and recurring phenomenon in oceanic systems, with frequent occurrence and abundance across mid-to-high latitudes and subtropical regions. Mesoscale eddies are thought to arise through various mechanisms, including direct wind forcing, baroclinic instability, and barotropic instability, with direct wind forcing being a predominant factor [6]. A fundamental characteristic of these eddies is their rotational speed or eddy intensity, which exhibits notable seasonal variability. In the Northern Hemisphere, eddy intensity peaks in spring and reaches a minimum in autumn, whereas the Southern Hemisphere displays an inverse pattern. This variability is especially pronounced within the tropical–subtropical transition zones (15°–30° latitude) in regions such as the western Pacific, northwestern Atlantic, and eastern Indian Ocean, primarily due to rapid changes in baroclinic instability [7]. Most mesoscale eddies in the upper ocean layer, whether cyclonic or anticyclonic, have lifespans shorter than 30 days. Mesoscale eddies in the ocean typically propagate westward at the phase velocity of non-dispersive baroclinic Rossby waves. Cyclonic and anticyclonic eddies tend to exhibit slight deviations towards the polar and equatorial regions, respectively. The majority of these eddies are characterized by nonlinear dynamics [8]. Additionally, the size of cyclonic eddies remains fairly constant with depth, whereas anticyclonic eddies exhibit a slight reduction in size [9]. The influence of these eddies on sea surface wind speeds and heat fluxes is particularly pronounced during the winter season, whereas their impact on precipitation rates is more substantial in the summer. The magnitude of sea surface temperature (SST) anomalies induced by mesoscale eddies is modulated by the strength of the background SST gradient field, with a stronger gradient in winter leading to larger anomalies and a weaker gradient in summer resulting in smaller anomalies [10]. Furthermore, mesoscale eddies play a significant role in the dynamics of monsoon rainfall, influencing the onset and cessation of the rainy season within the Indian monsoon region [11]. These eddies are characterized by horizontal scales ranging from tens to hundreds of kilometers and vertical penetration to depths exceeding a thousand meters below the thermocline [12]. As a result, they play a significant role in the transport of momentum and energy within the ocean [13,14]. Their lifespan typically spans from several days to weeks, with generation and dissipation processes intricately linked to oceanic circulation dynamics and the broader influences of climate change [15]. Governed by the quasi-geostrophic vorticity conservation equation, the motion of mesoscale eddies holds critical importance in altering the spatial distribution characteristics of oceanic elements [16,17,18].

Mesoscale eddies wield a profound influence on the spatial distribution of marine biota and chemical constituents, while also modulating the physical attributes of seawater as they traverse through the oceanic realm [19,20]. Their dynamic motion offers invaluable insights into diverse phenomena, encompassing the propagation of acoustic signals, as well as the conveyance of oceanic energy and material [19,21]. Consequently, a substantial body of research has been dedicated to achieving the precise delineation of these eddies. Since the 1990s, the availability of SLA data from satellite altimetry has revolutionized global eddy identification, enabling a comprehensive assessment of mesoscale features [22]. Leveraging technological advancements and machine learning methodologies, particularly deep learning, has emerged as a potent tool for automating and refining eddy detection processes [23]. Notably, the YOLOv3 algorithm, rooted in deep learning principles, has demonstrated efficacy in rapid and accurate eddy identification. Furthermore, recent innovations, exemplified by the Sectional Paths in Time series (SPT) method developed by the Guangzhou Institute of Oceanography, offer a synchronized approach to eddy detection, addressing challenges associated with their transient nature, small-scale spatial extent, and irregular morphology [24].

In current research, the identification of mesoscale eddies commonly relies on two primary methodologies leveraging satellite-derived data. The first approach involves the assimilation of satellite data collected along tracks into oceanic models. Subsequently, these datasets undergo processing, often utilizing advanced techniques such as the Data Unification and Altimeter Combination System (DUACS), employing optimal interpolation methods to generate gridded data. This facilitates the automated identification of eddies within the oceanic system. Conversely, the second methodology entails a direct fitting of the original satellite along-track data. This approach bypasses the need for data importation into oceanic models and, thus, mitigates potential information loss and error accumulation. Furthermore, it enables the real-time visualization of sea surface height data, enhancing the temporal resolution of eddy identification. Common techniques within this category include B-spline surface fitting and orthogonal polynomial surface fitting [13]. This paper introduces a novel approach utilizing Legendre polynomial surface fitting for eddy identification. Through a comparative analysis of computational efficiency and accuracy, contrasting it with the quasi-uniform bicubic B-spline surface fitting, we demonstrate the efficacy of the Legendre polynomial method for mesoscale eddy identification.

This paper is structured into several sections to present a comprehensive analysis. Section 2 delineates the origin of the SLA data employed in this study along with the research area’s geographical scope. In Section 3, we commence with a discussion on one-dimensional curve fitting, elaborating on the methodologies of B-spline surface fitting and Legendre polynomial surface fitting. Moving forward, Section 3.2 conducts ideal experiments utilizing simulated data to evaluate the performance of both fitting techniques, followed by a comparative analysis and ensuing discussion. Subsequently, practical experiments are undertaken employing real-world trajectory data to further assess the efficacy and potential fitting errors associated with Legendre polynomial surface fitting. These results are juxtaposed against those derived from gridded reanalysis data. Lastly, Section 4 encapsulates the primary findings and conclusions drawn from this research endeavor.

2. Data and Methodology

2.1. Study Region and Data

The satellite orbit point latitude and longitude data employed in this study originate from the Copernicus Marine Environment Monitoring Service (CMEMS), accessible at https://marine.copernicus.eu/, accessed on 30 March 2022. CMEMS aggregates data from six satellite altimeters: HY-2B, Cryosat-2 (in its new orbit), Sentinel-3A, Sentinel-3B, Jason-3, and Saral/AltiKa. Our research utilizes global near-real-time products derived from Level 3 initial along-track data, with detailed specifications outlined in

Table 1. Detailed introduction of satellite along-track data.

Sample	Regional Scope	Longitude and Latitude Range	Timer Series	Level
SLA	Antarctic Circumpolar Current	56°S–62°S 55°W–80°W	20220301000000—20220313235959	3

. For illustrative purposes, we have selected a dataset spanning 13 days, from 1 March 2022 to 13 March 2022. Our focus lies within the Antarctic Circumpolar Current region, delimited by latitudes ranging from 56°S to 62°S and longitudes from 55°W to 80°W, as delineated in Figure 1a. To facilitate ideal experimentation, a bounded area spanning 56°S to 62°S and 65°W to 75°W has been designated. In order to construct an ideal model representing mesoscale eddies, we employ a two-dimensional Gaussian function, computing corresponding function values for each grid point in latitude and longitude. These values are utilized to form a Gaussian distribution surface. During fitting experiments, when surface velocity equals zero, we extract function values from the model surface at satellite orbit coordinates, substituting them as ideal replacements for sea surface height data.

The flow velocity and sea surface height data used in this study were obtained from the HYCOM model (https://www.hycom.org/), accessed on 13 July 2024. Specifically, data spanning 31 days, from 18 February 2022 to 20 March 2022, were processed to compute monthly averages. The region is delimited by latitudes 56°S to 62°S and longitudes 65°W to 75°W, as illustrated in Figure 2. The arrow vector signifies the direction of the flow field, while its magnitude is proportional to the flow velocity.

2.2. Methods

2.2.1. Cubic B-Spline Curve (One-Dimensional)

Bezier curves were originally proposed by Bezier in 1962 and later optimized by Cox [25], De Boor, and others [26,27], leading to the development of B-spline curves, which addressed the drawback of limited local flexibility in Bezier curves [28]. The expression of B-spline curves is primarily determined by the degree, control points, and knot vectors. Given the existence of $n$ + 1 control points and $k$ + 1 knot vectors on a two-dimensional plane, where the control point sequence is denoted as $P_0, P_1, \cdots ,P_n$, and the knot vector is denoted as $U=\left\{{\mathrm{u}}_0,{\mathrm{u}}_1,\cdots ,{\mathrm{u}}_k\right\}$, the expression of a p-degree B-spline curve $J$ determined by these $n$ + 1 control points and $k$ + 1 knot vectors is as follows [29]:

$$J\left(u\right)={\sum }_{i=0}^n{N}_{i,p}\left(u\right)P_i,0\ll u\ll 1$$

(1)

where $N_{i,p}$ is the B-spline basis function, the subscript $i$ denotes the $i$-th knot interval, and it satisfies [30]:

$$N_{i,0}\left(u\right)=\left\{\begin{array}{l}1, u_i\ll u<u_{i+1}\\ 0, Others\end{array}\right.$$

(2)

The relationship among $n$, $k$, and $p$ in Equation (1) must satisfy $k$ = $n$+ $p$ + 1. As the degree $p$ of the B-spline curve increases, its smoothness improves, but the computational complexity also increases, leading to low computational efficiency. Therefore, in the experiments conducted in this paper, the degree of the B-spline curve is fixed at 3. The control of the change in the order of the B-spline surface is achieved by altering the number of knot vectors. In other words, conceptually, the order of the curve is equivalent to the number of control points, and the order of the curve $J$ is referred to as $n$ + 1.

B-spline curves are classified into three categories: uniform B-spline curves, quasi-uniform B-spline curves, and non-uniform B-spline curves [31]. Among them, uniform B-spline curves have poor fitting effects at the edges due to the lack of a specific structure in the distribution of knot vectors. Non-uniform B-spline curves have randomly distributed knot vectors. Therefore, quasi-uniform B-spline curves are chosen for this experiment.

The recursive formula for quasi-uniform B-spline basis functions $N_{i, p}$ is as follows [32]:

$$N_{i,p}\left(u\right)={\displaystyle \frac{u-u_i}{u_{i+p}-u_i}}{N}_{i,p-1}\left(u\right)+{\displaystyle \frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}}{N}_{i+1,p-1}\left(u\right)$$

(3)

Firstly, we define ${\displaystyle \frac{0}{0}}=0$. To calculate $N_{i,p}$, one needs to compute $N_{i,p-1}$ and $N_{i+1,p-1}$ first. In the programming implementation, one calculates the 0th B-spline basis functions for each knot interval, then solves for the 1st B-spline basis functions, and gradually recurses according to the formula to obtain all the basis functions.

2.2.2. Bicubic Quasi-Uniform B-Spline Surface (Two-Dimensional)

For a bicubic quasi-uniform B-spline curve fitting in each direction, we now extend to a bicubic quasi-uniform B-spline surface fitting with one independent variable in the east–west direction and another independent variable in the north–south direction.

For any two-dimensional field function $Z\left(x,y\right)$, suppose there are $N$ discrete point $\left(x_1,y_1\right),\left(x_2,y_2\right),\cdots ,\left(x_N,y_N\right)$ within the known interval, with corresponding function values $Z_1,Z_2,\cdots ,Z_N$. Assuming the order in the east–west direction is m + 1 and the order in the north–south direction is n + 1, meaning the number of control points grid $B$ is $(m+1)\times (n+1)$, then the corresponding fitting expression is as follows [13]:

$$S_{i,j}\left(u,v\right)=\sum _{i=0}^3\sum _{j=0}^3{B}_{i,j}{N}_{i,3}\left(u\right)N_{j,3}\left(v\right),i=\mathrm{0,1},\cdots ,m-3,j=\mathrm{0,1},\cdots ,n-3$$

(4)

$N_{i,3}\left(u\right)$ and $N_{j,3}\left(v\right)$ are the bicubic quasi-uniform B-spline basis functions in the east–west and north–south directions, respectively, which can be obtained from Equations (2) and (3).

2.2.3. Legendre Polynomials (One-Dimensional)

During his youth, Legendre utilized spherical harmonic equations based on Laplace’s equation, leading to the Legendre equation, and subsequently derived Legendre polynomials from the degeneration of the solutions to the equation [31].

Below is the general form of Legendre polynomial basis functions. For any $x$ value on the x-axis, the Legendre polynomial at this $x$ value is:

$$P_n\left(x\right)={\displaystyle \frac{1}{2^{n}n!}}{\displaystyle \frac{d^n}{d{x}^n}}{\left(x^2-1\right)}^n,n=\mathrm{1,2},3,\dots$$

(5)

The specific expression obtained by unfolding each order is:

$$\begin{gathered} P_0\left(x\right)=1 \\ {P}_1\left(x\right)=x \\ {P}_2\left(x\right)={\displaystyle \frac{1}{2}}\left({3x}^2-1\right) \\ \dots \dots \\ {P}_n\left(x\right)=\sum _{k=0}^{\left[{\displaystyle \frac{n}{2}}\right]}{\displaystyle \frac{{\left(-1\right)}^k\left(2n-2k\right)!}{2^{n}k!\left(n-k\right)!\left(n-2k\right)}}{x}^{n-2k},n=\mathrm{1,2},3,\dots \end{gathered}$$

(6)

Legendre polynomials possess the following properties:

Property 1: On the interval [−1,1], Legendre polynomials $\left\{P_n\left(x\right)\right\}$ form an orthogonal function series with respect to the weight function $\omega \left(x\right)\equiv 1$, exhibiting orthogonality, and satisfying the following relationship:

$$\left(P_n,P_m\right)={\int }_{-1}^1{P}_n\left(x\right)P_m\left(x\right)dx=\left\{\begin{array}{c}0 , n\ne m\\ {\displaystyle \frac{2}{2n+1}},n=m \end{array}\right.$$

(7)

Property 2: Legendre polynomials have a recursive relationship, satisfying the following formula:

$$\left(n+1\right)P_{n+1}\left(x\right)=\left(2n+1\right)x{P}_{n+1}\left(x\right)-{nP}_{n-1}\left(x\right),n=\mathrm{1,2},3,\dots$$

(8)

Orthogonality ensures the independence of each coefficient, thereby preventing errors from interfering with one another. Furthermore, orthogonality often leads to a sparse or diagonalized coefficient matrix, which significantly improves computational efficiency. Orthogonality mitigates the risk of ill-conditioned matrix equations and enhances the accuracy of the results.

The use of recursive relationships enhances the computational efficiency of higher-order Legendre polynomials. By utilizing the initial polynomial, subsequent higher-order polynomials can be computed iteratively through the recursive relationship, obviating the need for solving explicit expressions for each polynomial. This approach not only simplifies the computation but also reduces the accumulation of numerical errors that are common in direct calculations of high-order polynomials.

2.2.4. Legendre Polynomial Surfaces (Two-Dimensional)

Given the Legendre polynomial basis functions for fitting along one direction, we now extend to fitting Legendre polynomial surfaces with one independent variable in the east–west direction and another independent variable in the north–south direction.

For any two-dimensional field function $Z\left(x,y\right)$, suppose there are N discrete points $\left(x_1,y_1\right),\left(x_2,y_2\right),\cdots ,\left(x_N,y_N\right)$ within the known interval, with corresponding function values $Z_1,Z_2$, …, $Z_N$. Assuming the degrees in the east–west and north–south directions are respectively denoted as $M_0$ and $N_0$, and both $M_0$ and $N_0$ are much smaller than the number of discrete points $N$, then the corresponding fitting expression is [31]:

$$F\left(x_i,y_i\right)=\sum _{m=0}^{M_0}\sum _{n=0}^{N_0}{A}_{m,n}{\phi }_m\left(x_i\right){\psi }_n\left(y_i\right),i=\mathrm{1,2},\cdots ,N$$

(9)

${\phi }_m\left(x_i\right)$ and ${\psi }_n\left(y_i\right)$ are Legendre polynomials in the east–west and north–south directions, respectively, as shown in Equation (6).

Given the predetermined orders $m$ + 1 and $n$ + 1 in the east–west and north–south directions, respectively, i.e., the degrees of the Legendre polynomial surface are already known, the approach to surface fitting involves finding a Legendre polynomial surface with order $m$ + 1 in the east–west direction and order $n$ + 1 in the north–south direction. This surface minimizes the sum of squared errors between the set of known data points Z and the corresponding fitted values $\widehat{z}$ on the fitting surface, i.e., $min{\left(Z-\widehat{z}\right)}^2$.

In this ideal experiment, the data point set $Z$ represents the function values corresponding to the initial satellite orbit latitude and longitude coordinates on a Gaussian surface. Therefore, the selection of the optimal surface is transformed into a least squares problem, seeking to $min{\left(Z-\widehat{z}\right)}^2$ to determine the coefficients in the fitting surface expression, thereby establishing the best-fitting surface corresponding to the data point set.

The sum of squared errors is defined as:

$$\epsilon \left(A_{\mathrm{0,0}},A_{\mathrm{1,0}},\cdots ,A_{k,s},\cdots ,A_{K_0,S_0}\right)=\sum _{i=1}^N{\left[\sum _{m=0}^{M_0}\sum _{n=0}^{N_0}{A}_{m,n}{\phi }_m\left(x_i\right){\psi }_n\left(y_i\right)-Z_i\right]}^2$$

(10)

According to the principle of least squares, taking the partial derivative of each variable yields:

$${\displaystyle \frac{\partial \epsilon }{\partial A_{m,n}}}=2\sum _{i=1}^N\left[Z_i-\sum _{k=0}^{M_0}\sum _{s=0}^{N_0}{A}_{k,s}{\phi }_k\left(x_i\right){\psi }_s\left(y_i\right)\right]{\phi }_m\left(x_i\right){\psi }_n\left(y_i\right)=0$$

(11)

The definition of orthogonal polynomials includes:

$$\begin{array}{l}{\displaystyle \sum _{i=1}^N{\phi }_k(x_i){\phi }_m(y_i)=0,k\ne m},\hfill \\ {\displaystyle \sum _{i=1}^N{\psi }_s}(y_i){\psi }_n(y_i)=0,s\ne n,\hfill \end{array}$$

(12)

Therefore, the expansion coefficient is:

$$A_{m,n}={\displaystyle \frac{{\sum }_{i=1}^{N}Z\left(x_i,y_i\right){\phi }_m\left(x_i\right){\psi }_n\left(y_i\right)}{{\sum }_{i=1}^N{\phi }_m^2\left(x_i\right){\psi }_n^2\left(y_i\right)}}$$

(13)

Before conducting surface fitting, it is necessary to determine the orders $m$ + 1 and $n$ + 1 for polynomial fitting in the east–west and north–south directions, respectively. Different orders correspond to different fitting surfaces. To select the optimal orders, this experiment adopts the method of cross-validation to compare and choose among different combinations of orders.

2.2.5. Cross-Validation

The process of cross-validation involves dividing the data into two parts: one for model training and the other for evaluating the performance of the trained model. The former is commonly referred to as the training set, while the latter is known as the validation set. Since the data in the training and validation sets are different, and only the validation set data are used to calculate errors, ensuring a sufficient number of data samples ensures that the estimation of errors on the validation set is close to the true situation.

According to the differences in dataset partitioning, cross-validation methods can be classified into four types: leave-one-out [32], leave-P-out [33], hold-out [34], and K-fold cross-validation [35]. In this experiment, a relatively simple implementation of 10-fold cross-validation is adopted. The specific workflow of 10-fold cross-validation is as follows:

• Randomly divide the existing data samples into 10 subsets;
• Select one subset as the validation set and the remaining nine subsets as the training set;
• Vary the orders of Legendre polynomials in the east–west and north–south directions in the training set, and obtain models trained with different combinations of orders;
• Import the validation set data into each model and record the mean absolute errors (MAEs) of the validation set data under different combinations of orders of Legendre polynomials;
• Repeat the process by iteratively selecting each of the nine training subsets as the validation set and the original validation set as the training set, thus obtaining new validation and training sets;
• Repeat the cross-validation process 10 times in total. Calculate the average MAE for each combination of orders of Legendre polynomials obtained from the 10 cross-validation iterations.
• Compare the average MAE values for different combinations of orders, and select the combination with the smallest MAE as the optimal combination of orders.

3. Results

3.1. Comparison of Static Model Bicubic B-Spline Surface Fitting and Legendre Polynomial Surface Fitting

Figure 3a demonstrates the computation of function values corresponding to each grid point of latitude and longitude utilizing a two-dimensional Gaussian function, thereby constructing a Gaussian distribution surface as the foundational model for mesoscale eddies. The accompanying Figure 3b portrays the surface values aligned with the satellite nadir point orbit.

Through the manipulation of controlled variables and the application of cross-validation methodologies, the orders of the east–west and north–south directions were systematically varied from 5 to 13, while the duration of observational data remained consistent at 1 to 13 days. For each combination of orders, surface fitting was executed employing both Legendre polynomial and bicubic quasi-uniform B-spline methods. The evaluation of model performance was facilitated by the mean absolute error (MAE) as the designated metric for cross-validation. The resulting MAE values corresponding to distinct combinations of orders of both spatial directions were documented for fixed data durations. Subsequent modifications to the duration of observational data enabled comparison of the minimum MAE for each data duration set. Notably, due to the limited availability of observed data points during the initial 4 days, yielding elevated MAE values, these data points were excluded from the analysis. Conversely, employing data spanning more than 13 days could yield distorted fitting outcomes attributable to the substantial time gap between the earliest and latest observations, as well as errors arising from the motion of eddies.

Table 2. Data quantity for different days and optimal MAE values corresponding to two fitting methods (cm).

	Days
	5	6	7	8	9	10	11	12	13
Legendre	2.25	0.63	0.49	1.21	0.33	0.33	0.32	0.31	0.30
B-spline	1.90	0.15	0.11	0.90	0.04	0.05	0.04	0.04	0.04
data Volume	2608	3097	3439	3863	4427	5113	5695	5933	6454

provides a comprehensive overview, detailing the quantity of data points for days 5 through 13 alongside the corresponding optimal MAE values for both Legendre polynomial and bicubic quasi-uniform B-spline fitting procedures [13].

Analysis of Table 2 reveals a consistent trend wherein the MAE decreases gradually with the augmentation of data days, barring an anomalous surge observed when employing 8 days of data. However, beyond 8 days, the reduction in error demonstrates diminishing returns, particularly evident from 9 days onward. Consequently, for mesoscale eddies characterized by the static model, 9-day data all the way up to 13-day data are deemed optimal, with both fitting methodologies proving viable. Notably, under equivalent data day conditions, the B-spline fitting method exhibits a marginally lower error rate in the ideal scenario, indicative of superior fitting precision.

In real-world scenarios, data are often subject to various forms of interference and noise, rendering them inherently less pristine. Understanding the ramifications of noise on fitting procedures is crucial for comprehending the model’s performance under realistic conditions and accurately gauging its efficacy. Noise imbues data with heightened complexity and uncertainty, necessitating the evaluation of model resilience amidst such disturbances, rather than succumbing to overfitting.

Subsequently, varying degrees of random noise will be introduced to investigate the corresponding trend in MAE derived from 10-fold cross-validation for both fitting methodologies as the noise levels escalate. The ensuing findings are presented in Figure 4.

Analysis of

Table 3. Optimal MAE values corresponding to two fitting methods for 9 days’ data with different signal-to-noise ratios of noise added (cm).

	Signal-to-Noise Ratio
	0	1%	2%	3%	4%	5%	6%	7%	8%
Legendre	0.30	0.44	0.65	0.87	1.08	1.35	1.55	1.81	2.05
B-spline	0.04	0.27	0.51	0.75	1.03	1.29	1.51	1.72	2.04

underscores that in the absence of added noise, B-spline fitting exhibits an optimal MAE nearing 0, indicative of markedly superior performance compared to Legendre polynomial fitting. However, with escalating noise levels, the optimal MAEs for both methods proportionally increase. Notably, as the signal-to-noise ratio reaches 4%, the advantageous performance of B-spline fitting diminishes. Upon further amplification of the signal-to-noise ratio to 8%, the optimal MAEs for both methodologies converge to approximate equality.

Subsequently, employing the Legendre polynomial surface fitting technique, utilizing 9 days of data, the orders m and n are systematically varied in the east–west and north–south directions, respectively, for 10-fold cross-validation. The outcomes of the cross-validation, depicted in Figure 4, visually indicate that the error levels are around 1% when the orders in the east–west direction exceed or equal 11 and in the north–south direction exceed or equal 9.

Utilizing the Legendre polynomial fitting approach with a dataset spanning 9 days and employing orders (11, 9) for the east–west and north–south directions, respectively, the fitting process is illustrated in Figure 5a. The ensuing errors at each satellite nadir point are delineated in Figure 5b.

3.2. Comparison of Bicubic B-Spline Surface Fitting and Legendre Polynomial Surface Fitting for Actual Measurement Data

We followed the steps of the ideal experiment and used the along-trackdata. Through the manipulation of controlled variables and the application of cross-validation methodologies, the orders of the east–west and north–south directions were systematically varied from 5 to 13, while the duration of the observational data remained consistent at 1 to 13 days. For each combination of orders, surface fitting was executed employing both Legendre polynomial and bicubic quasi-uniform B-spline methods. The evaluation of model performance was facilitated by the mean absolute error (MAE) as the designated metric for cross-validation. The resulting MAE values corresponding to distinct combinations of orders of both spatial directions were documented for fixed data durations. Subsequent modifications to the duration of observational data enabled comparison of the minimum MAE for each data duration set. Notably, due to the limited availability of observed data points during the initial 4 days, yielding elevated MAE values, these data points were excluded from the analysis.

Table 4. Optimal MAE values corresponding to two fitting methods for data with different days (cm).

	Days
	5	6	7	8	9	10	11	12	13
Legendre	4.43	3.16	2.65	3.10	2.14	2.35	2.40	2.42	2.58
B-spline	4.03	2.23	2.15	2.83	2.16	2.24	2.30	2.35	2.51

furnishes a comprehensive summary, delineating the optimal MAE values associated with both Legendre polynomial and bicubic quasi-uniform B-spline fitting methodologies across a duration spanning 5 days to 13 days.

Analysis of Table 4 reveals a consistent trend across both the Legendre polynomial and B-spline fitting methods concerning the alteration in error with increasing data days. However, unlike the controlled ideal experiment scenario, the discrepancy in error between Legendre polynomial and B-spline fitting diminishes notably when utilizing real measured data, particularly when the data days exceed or equal to 9. Consequently, the performance of both methods converges, exhibiting nearly identical fitting outcomes. Beyond 9 days of data, the reduction in error becomes less pronounced. Hence, for real measured data, opting for data collected from days 9 to 13 proves advantageous for achieving optimal fitting results for eddies, regardless of the chosen method. The motion of mesoscale eddies introduces greater errors as the number of data days utilized increases. Subsequently, the Legendre polynomial surface fitting method was employed, utilizing data from 9 consecutive days, while systematically varying the orders m and n along the east–west and north–south directions, respectively, for a 10-fold cross-validation procedure. The outcomes of this cross-validation process are visually depicted in Figure 6, showcasing a point where the error improvement becomes negligible as the orders in both the east–west and north–south directions reach or exceed 10. Beyond this threshold, escalating the orders would entail a substantial increase in computational complexity without significant error reduction. Hence, a pragmatic decision was made to adopt Legendre polynomial fitting with an order combination of (10, 10), leveraging 9 days of data for the fitting process. The results of this fitting endeavor are presented in Figure 7a, accompanied by Figure 7b illustrating the corresponding errors at each satellite’s sub-satellite point.

The study employed a desktop computer featuring an Intel(R) Core(TM) i7-9700 processor and 32.0 GB (Lenovo Group, Beijing, China) of RAM. It is important to acknowledge that the runtime of the program may vary based on the configuration of the computer used. However, on the same computer, the time required to complete the fitting using the Legendre method is consistently less than that required by the bicubic quasi-uniform B-spline method. Following multiple repeated experiments, for the area depicted in Figure 7, the time required for Legendre polynomial surface fitting is approximately 74.5% of the time needed for quasi-uniform bicubic B-spline surface fitting.

3.3. Hypothesis Testing

In empirical oceanographic research, data often manifest notable uncertainties, presenting challenges to strict adherence to a normal distribution. Consequently, deeper analysis commonly involves the utilization of normal distribution histograms and Probability–Probability (P-P) plots. In a normal distribution histogram, the robustness of the data’s normality is assessed based on its resemblance to a bell curve, characterized by a prominent central peak and tapered tails. Conversely, in a P-P plot, the degree of alignment between the observed data points and the theoretical distribution line serves as an indicator of the extent to which the dataset adheres to a normal distribution, whether approximately or precisely [13].

Figure 8 and Figure 9 exhibit the frequency distribution histogram and P-P plot of the data errors, respectively.

The histogram of the normal distribution typically displays a bell-shaped curve, characterized by a peak at the center and diminishing frequencies towards the tails. Additionally, the P-P plot demonstrates a close alignment between the observed cumulative probability (P) and the expected normal cumulative probability (P), exhibiting an approximation to a linear relationship. Line represents the expected normal cumulative probability, curve represents the observed cumulative probability. Consequently, despite the potential deviation from strict adherence to a normal distribution, the error data can be reasonably regarded as approximately normal.

3.4. Comparative Analysis of Grid Data and Along-Track Data

The fitting resolution is 0.01° × 0.01°, contrasting with the coarser resolution of the existing gridded reanalysis data at 0.25° × 0.25°. Post-fitting, the results are expected to present enhanced spatial granularity, capturing intricate variations across the domain. This finer modeling facilitates clearer visualizations and enables a denser spatial representation of the data, which is instrumental for spatial analyses.

Figure 10 showcases the outcomes of Legendre polynomial surface fitting utilizing both gridded reanalysis data and along-track data. Specifically, Figure 10a,c illustrate the original datasets for the same region and day, with Figure 10b,d presenting the corresponding fitting results. Notably, while noticeable disparities exist between the spatial distributions of the gridded reanalysis and along-track data, the resultant Legendre polynomial fittings exhibit minimal deviations.

In Figure 11a, the disparity between the interpolated values from the gridded data and the measured values at the satellite observation points is depicted, while Figure 11b showcases the distinction between the results derived from fitting the gridded data and the results derived from fitting the actual along-track data. The mean bias for these two products is 0.34, while the standard deviation is 2.98. Notably, a comparison between Figure 10b and Figure 11b reveals a significant error at the center of the mesoscale vortex, whereas the error diminishes towards its periphery. This discrepancy can be attributed to the grid data’s tendency to inaccurately assess data with prominent gradients, consequently resulting in errors.

Vorticity is a crucial physical parameter that delineates the rotational dynamics of fluids, pivotal for comprehending ocean circulation and vortex structures. The calculation of vorticity facilitates the analysis of eddy characteristics across diverse oceanic scales, thereby elucidating dynamic processes and phenomena. Utilizing sea surface height data for vorticity computation provides foundational support for advanced research in ocean physics and vortex dynamics.

$$\left\{\begin{array}{c}u=-{\displaystyle \frac{1}{\rho f}}{\displaystyle \frac{\partial P}{\partial y}} \\ v= {\displaystyle \frac{1}{\rho f}}{\displaystyle \frac{\partial P}{\partial x}}\end{array}\right.$$

(14)

$${P=P}_0-\rho gZ$$

(15)

$$\xi ={\displaystyle \frac{\partial v}{\partial x}}-{\displaystyle \frac{\partial u}{\partial y}}$$

(16)

Initially, we utilized Formulas (14) and (15) to compute the geostrophic flow velocity, followed by Formula (16) to determine the vorticity. Here, $g$ denotes gravitational acceleration, $f$ represents the geostrophic parameter, $P$ stands for pressure, $\rho$ signifies density, and $Z$ denotes the SLA data.

Due to the absence of velocity data in along-track data, discrepancies arise between geostrophic flow calculations derived from SLA data and actual flow fields. Therefore, it is imperative to initially utilize reliable flow field data for accurate vorticity field computation, followed by quantitative analysis of vorticity fields obtained through different methods. In this part, data within the specified spatial range from the HYCOM dataset is selected, and a 31-day period centered on 9 days is averaged monthly. The reliable vorticity field is computed using Formula (16), as depicted in Figure 12. From Figure 2, it is apparent that areas characterized by higher sea surface heights relative to their surroundings exhibit positive vorticity, whereas regions with lower sea surface heights exhibit negative vorticity.

The sea surface elevation data from HYCOM is first interpolated to correspond with the latitude and longitude coordinates of satellite orbit observations, thereby establishing along-track sea surface elevation. Subsequently, two distinct methods are employed to extrapolate sea surface height values across the entire studied region. Using Formulations (14)–(16), the vorticity fields are then computed. Figure 13a presents the vorticity obtained the results of using Legendre polynomial surface fitting, demonstrating an average absolute error of $3.67\times {10}^{-6} s^{-1}$ and an average error of $-1.34\times {10}^{-7} s^{-1}$ relative to the vorticity field depicted in Figure 12. In addition, Figure 13b showcases the vorticity outcomes derived from the bicubic quasi-uniform B-spline surface fitting, exhibiting an average absolute error of $2.90\times {10}^{-6} s^{-1}$ and an average error of $-9.19\times {10}^{-8} s^{-1}$ compared to Figure 12. Velocity data is interpolated to match the latitude and longitude coordinates corresponding to the along-track data. Subsequently, Legendre polynomial surface fitting is applied to establish the velocity across the entire region. Formula (16) is then used to compute the vorticity. Figure 13c illustrates these results, showing an average absolute error of $2.99\times {10}^{-6} s^{-1}$ and an average error of $5.86\times {10}^{-8} s^{-1}$ compared to the vorticity field depicted in Figure 12. Similarly, velocity data is interpolated to align with the latitude and longitude coordinates corresponding to the satellite orbit data, followed by bicubic quasi-uniform B-spline surface fitting to determine the velocity across the entire region. Formula (16) is employed to calculate the vorticity. Figure 13d presents these results, with an average absolute error of $2.54\times {10}^{-6} s^{-1}$ and an average error of $9.10\times {10}^{-8} s^{-1}$ relative to the vorticity field shown in Figure 12.

The vorticity fields derived from the fitted sea surface height or sea surface velocity demonstrate strong consistency with the reference vorticity field, with errors falling within acceptable limits. This result confirms the feasibility of employing both Legendre polynomial surface fitting and bicubic quasi-uniform B-spline surface fitting methods to interpolate satellite orbit height data and compute vorticity. It is noteworthy that vorticity calculated from fitted data appears smoother than the reference vorticity field, yet maintains overall consistency. Small vorticity regions characterized by high local variability tend to merge into larger vorticity structures, thereby introducing localized distortions in the vorticity field.

Next, vorticity will be computed using along-track data and gridded reanalysis data. Figure 14 illustrates the vorticity computations using different methods as follows: (a) depicts vorticity calculated based on Legendre polynomial surface fitting of satellite along-track data, (b) shows vorticity calculated based on bicubic quasi-uniform B-spline surface fitting of satellite along-track data, (c) presents vorticity derived from gridded reanalysis data, and (d) displays vorticity computed using Legendre polynomial surface fitting of gridded reanalysis data. The vorticity calculated directly from gridded reanalysis data exhibits consistency with the vorticity derived from sea surface height data fitted using along-track data. Moreover, the vorticity computed from sea surface height data using different fitting methods also demonstrates strong agreement. Variations in the calculated vorticity arise from disparities between the HYCOM dataset and satellite orbit data. Nevertheless, spatial correspondence can be observed between the two datasets; for instance, a positive vorticity region is evident in the upper right corner of the study area. Moreover, several primary positive and negative vorticity regions in the study area exhibit corresponding features.

Expanding the fitting domain facilitates an increase in the sample size, thereby enriching the dataset’s comprehensiveness and completeness. This expansion is crucial for enhancing the accuracy and reliability of the fitting model. By broadening the fitting area, the resulting dataset encompasses a wider range of oceanic features, incorporating not only the original localized seas but also the surrounding environmental and topographical contexts. This holistic approach integrates diverse factors that potentially influence the distribution of oceanic elements, thereby facilitating a nuanced understanding of their overall characteristics. Moreover, it enables a more precise analysis of spatial distribution patterns and trends, revealing interrelations and underlying mechanisms governing oceanic elements.

Expanding the fitting range to encompass 56°S–62°S and 55°W–80°W, Figure 15a presents satellite along-track data for the specified region, juxtaposed with Figure 15b illustrating gridded reanalysis data for the same area. Direct fitting of the enlarged area proves insufficient in yielding satisfactory results. To enhance fitting precision and minimize errors, the region is partitioned based on the fitting approach applied to the initial 6° × 10° area. Each smaller segment of 6° × 10° undergoes individual fitting, followed by trimming and seamless integration of overlapping sections. The resulting fitting outcomes are depicted in Figure 15c, with Figure 15d delineating the associated fitting errors. Thus, this method enables the precise fitting of larger areas by leveraging the accurate fitting of smaller regions. Its applicability extends beyond the scope illustrated in Figure 15 to encompass larger spatial domains. The mean bias for these two products is 0.41, while the standard deviation is 13.23.

4. Summary and Conclusions

In the realm of physical oceanography, the spatial distribution of sea surface height stands as a critical parameter for understanding ocean dynamics. Researchers often confront the choice between utilizing gridded reanalysis data, subject to potential information loss during processing, or employing direct fitting methods with satellite along-track data. However, existing direct fitting methodologies are somewhat constrained, necessitating the exploration of novel approaches to enhance efficiency and accuracy in sea surface height estimation.

Previous studies, notably by Xu et al., have advocated for B-spline surface fitting as the direct fitting method [13]. This study diverges from prior work by investigating the feasibility of Legendre orthogonal polynomial surface fitting, juxtaposed with bicubic quasi-uniform B-spline surface fitting. Initially, an idealized static model of mesoscale eddies is established to evaluate the performance of both methods. Through rigorous 10-fold cross-validation, it is demonstrated that, under ideal conditions, bicubic quasi-uniform B-spline surface fitting exhibits superior accuracy compared to Legendre polynomial surface fitting for equivalent order combinations. However, when introducing noise, representative of real-world observational data, both methods demonstrate comparable performance.

Subsequent 10-fold cross-validation using actual observational data reaffirms this parity, while emphasizing the simplicity and controllability of Legendre polynomial surface fitting. Additionally, hypothetical error testing in Legendre polynomial surface fitting validates its viability for practical applications.

Moreover, a common challenge identified in the study is the decrease in accuracy when fitting over large spatial extents. To address this, a novel approach is proposed: the accurate fitting of sea surface height over expansive spatial domains through region stitching. This innovative method offers promising avenues for improving the accuracy and reliability of sea surface height estimation in oceanographic research.