Cross-dimensional multi-source satellite altimeter sea level anomaly joint estimation method and system

Abstract

The application discloses a cross-dimension multi-source satellite altimeter sea level anomaly joint estimation method and system, and relates to the field of ocean remote sensing.The joint estimation method comprises the following steps: obtaining multi-source observation data of a target sea area about a target period, forming a corrected joint observation data set through cross-dimension data processing and reconstruction operation and bias correction operation, and mapping the corrected joint observation data set to a corresponding regular space-time grid of the target sea area; based on a preset space-time search radius, constructing an observation operator, a background covariance matrix and an observation error covariance matrix corresponding to each grid point in the regular space-time grid, and generating a regular gridded sea level anomaly product of the target sea area about the target period through an optimal interpolation method without performing explicit system error correction.The application can consider large-scale background consistency and small and medium-scale structure resolution capability on multiple spatial scales, and improves the spatial continuity, statistical consistency and overall stability of the sea level anomaly estimation result.

Description

Technical Field

[0001] This invention relates to the field of marine remote sensing technology, specifically to a method and system for joint estimation of sea level anomalies using multi-dimensional, multi-source satellite altimeters. Background Technology

[0002] Sea level anomaly (SLA) is an important physical quantity characterizing ocean circulation structure, mesoscale and sub-mesoscale dynamic processes, and has wide applications in marine scientific research, operational ocean monitoring, and numerical ocean model assimilation. Accurate characterization of the spatiotemporal distribution of SLA can provide key physical constraints for ocean circulation diagnosis, energy cascade analysis, and ocean forecasting.

[0003] Existing methods for inverting and reconstructing sea level anomalies mainly rely on satellite altimeter observation data. These methods typically employ spatial interpolation, optimal interpolation, or variational estimation to map discrete multi-source observations onto a regular grid, thereby obtaining a continuous sea level anomaly field. On longer temporal and spatial scales, these methods can, to some extent, meet the needs of operational applications.

[0004] Furthermore, traditional methods primarily rely on one-dimensional along-orbit satellite altimeter observations. Because these observations are distributed along satellite orbits, their spatial sampling exhibits significant directionality and discontinuity, with sparse sampling in the lateral direction. This results in limited statistical constraints on lateral structures during the inversion process, making it difficult to fully resolve small- to medium-scale ocean dynamic structures. Even with increased observation density through multi-satellite combinations, such observations still suffer from spatial dimensional limitations in a statistical sense.

[0005] With the development of interferometric altimeter technology, wide-swath two-dimensional interferometric satellite altimeters can acquire high spatial resolution two-dimensional sea surface height observations, providing a new observational tool for fine characterization of small- to medium-scale ocean processes. However, while introducing high spatial resolution information, two-dimensional interferometric observations also bring more complex error statistics. Their observation errors often exhibit significant correlation and anisotropy in different spatial directions, and their spatial sampling density differs significantly from that of one-dimensional along-orbit observations.

[0006] Because one-dimensional along-track observations and two-dimensional interferometric observations differ fundamentally in spatial dimension, sampling structure, and error statistics, existing techniques typically employ methods such as staged processing, empirical weighting fusion, or explicit systematic error correction to fuse different observation data. These methods often introduce master-slave relationships or empirical weighting between different observations, making it difficult to impose consistent and controllable statistical constraints on cross-dimensional observations under a unified state space and statistical objective function. This can easily lead to an imbalance in information contributions at different spatial scales, thus limiting the full realization of the observational potential of two-dimensional interferometric altimeters.

[0007] Therefore, there is an urgent need for an estimation method that can take into account the statistical characteristics of observation errors in different dimensions under a unified statistical framework, and effectively regulate the information contribution at different spatial scales, so as to achieve a coordinated improvement in the consistency between spatial structure and statistics in the process of sea level anomaly inversion. Summary of the Invention

[0008] To address the aforementioned technical problems, this invention provides a method and system for joint estimation of sea level anomalies using multi-dimensional, multi-source satellite altimeters.

[0009] To achieve the above objectives, the present invention adopts the following technical solution:

[0010] This invention provides a method for joint estimation of sea level anomalies from multi-dimensional, multi-source satellite altimeters. Following steps S1 to S5, under a unified state space and a unified statistical estimation framework, it jointly models and estimates one-dimensional along-orbit observation data and two-dimensional interferometric observation data to generate a regular gridded sea level anomaly product for the target sea area over the target time period.

[0011] Step S1: Using a preset number of target satellites equipped with different altimeters, observation data corresponding to each observation point in the target sea area within the target time period are obtained respectively, and a one-dimensional along-orbit observation dataset and a two-dimensional interferometric observation dataset are formed; by applying quality screening and geophysical correction to the one-dimensional along-orbit observation dataset, a standard one-dimensional along-orbit observation dataset is formed, and at the same time, the two-dimensional interferometric observation dataset is spatially aggregated to form a downsampled two-dimensional interferometric observation dataset. Step S2: Map the standard one-dimensional track-along observation dataset and the downsampled two-dimensional interferometric observation dataset to a preset spatial frame, and then use reconstruction and bias correction operations in sequence to form a corrected joint observation dataset. Step S3: Based on the preset temporal resolution and preset spatial resolution, construct a regular spatiotemporal grid corresponding to the target sea area, and map the calibrated joint observation dataset to the regular spatiotemporal grid. Further, based on the preset spatiotemporal search radius, determine the spatiotemporal neighborhood of each grid point in the regular spatiotemporal grid, as well as the subset of calibrated joint observation data corresponding to each grid point.

[0012] Step S4: For each grid point, based on its corresponding subset of calibrated joint observation data, construct the observation operator, background covariance matrix, and observation error covariance matrix corresponding to that grid point, and combine the background covariance matrix and the observation error covariance matrix into a joint covariance matrix.

[0013] Step S5: For each grid point, based on its corresponding observation operator and joint covariance matrix, construct the joint statistical estimation objective function corresponding to the grid point without explicit systematic error correction, and calculate the optimal estimate of the sea level anomaly of the grid point by the optimal interpolation method. Then, based on the optimal estimate of the sea level anomaly of each grid point, generate a regular gridded sea level anomaly product of the target sea area for the target time period.

[0014] Furthermore, the observation points mentioned in step S1 include one-dimensional along-track observation points and two-dimensional interferometric observation points, wherein the one-dimensional along-track observation points correspond to one-dimensional along-track observation data, and the two-dimensional interferometric observation points correspond to two-dimensional interferometric observation data. Both the one-dimensional along-track observation data and the two-dimensional interferometric observation data include observation time, observation location, sea level anomaly value, mass marker, and orbit marker.

[0015] Further, step S1 applies quality screening and geophysical correction to the one-dimensional orbital observation dataset according to the following steps to form a standard one-dimensional orbital observation dataset: Step S1.1: Based on the preset quality identification parameters and preset physical parameter thresholds, the one-dimensional track observation dataset is divided into a valid one-dimensional track observation data subset and an invalid one-dimensional track observation data subset, and the invalid one-dimensional track observation data subset is deleted. Step S1.2: Apply standard geophysical corrections to each one-dimensional orbital observation data in the effective one-dimensional orbital observation data subset, and normalize each one-dimensional orbital observation data after applying standard geophysical corrections using a preset mean sea level model, thereby forming a standard one-dimensional orbital observation dataset.

[0016] Further, step S1 spatially aggregates the two-dimensional interferometric observation dataset according to the following steps to form a downsampled two-dimensional interferometric observation dataset: Step S1.1: Divide the target sea area using a preset downsampling factor to form non-overlapping and non-interval spatial aggregation windows, and determine the two-dimensional interferometric observation data subset corresponding to each spatial aggregation window; Step S1.2: For each spatial aggregation window, based on the preset quality identification parameters and preset physical parameter thresholds, filter out each effective two-dimensional interferometric observation point contained in the spatial aggregation window, and determine whether the number of effective two-dimensional interferometric observation points in the spatial aggregation window is greater than the preset threshold. If yes, proceed to step S1.3; otherwise, discard the spatial aggregation window and return to step S1.2 until all spatial aggregation windows have been judged. Step S1.3: For each spatial aggregation window, based on its corresponding subset of two-dimensional interferometric observation data, the super observation point corresponding to the spatial aggregation window is generated by taking the arithmetic mean of the sea level anomaly values ​​corresponding to all valid two-dimensional interferometric observation points within the spatial aggregation window. Each super observation point corresponds to a super observation data, and each super observation data includes the center observation time, center observation location, center sea level anomaly value, lateral distance indicator, and track index indicator. Step S1.4: Based on the orbital index identifier of each super observation point, the mean value of the central sea level anomaly of all super observation points in the same orbit is processed to form a downsampled two-dimensional interferometric observation dataset.

[0017] Furthermore, each target satellite operates along its corresponding orbit, and step S2 performs deviation correction as follows:

[0018] Step S2.1: Based on the track identifier of each observation point, determine the observation points contained in each track, and based on the observation position of each observation point, construct each observation point on the same track into a spatial linear geometric object, and each spatial linear geometric object has a spatial index. Step S2.2: Based on the spatial index of each spatial linear geometric object, use geometric calculation to determine the orbital intersection points between different orbits or between different orbital segments of the same orbit, and calculate the sea level anomaly value of each orbit at each orbital intersection point through neighborhood linear interpolation. Step S2.3: With the target reference satellite orbiting along the target reference orbit as the center, construct a star-shaped intersection topology and select the orbit intersections containing the target reference orbit as relevant intersections, that is, only retain the orbit intersections between each orbit and the target reference orbit. Step S2.4: Based on the sea level anomaly values ​​of each orbit at each relevant intersection point and the sea level anomaly values ​​of the target reference orbit at each relevant intersection point, the deviation values ​​of each orbit and the target reference orbit at all relevant intersection points are obtained. Then, the deviation prediction value corresponding to each orbit is generated by using a multinomial fitting method with the distance along the orbit as the independent variable. Finally, by subtracting the deviation prediction value of the sea level anomaly value corresponding to the observation point from the deviation prediction value of the orbit it belongs to, a corrected joint observation dataset with consistent system deviation is obtained.

[0019] Further, step S4 generates observation operators corresponding to each grid point using a spatiotemporal covariance function, wherein the spatiotemporal covariance function includes a basic spatiotemporal covariance function and a dual-scale spatiotemporal covariance function, and the basic spatiotemporal covariance function includes a basic spatiotemporal covariance kernel and a two-dimensional interferometric observation direction-related covariance term.

[0020] For each grid point, when the number of observation points in the spatiotemporal neighborhood of that grid point exceeds a preset threshold, a dual-scale spatiotemporal covariance function is used to generate the correlation between that grid point and the first grid point. Observation operators between observation points:

[0021] ;

[0022] in, For grid points and the first Observation operators between observation points For large-scale weights, For small-scale weights, For large-scale components, Small-scale components;

[0023] For each grid point, when the number of observation points in the spatiotemporal neighborhood of the grid point is not greater than the preset threshold for the number of observation points, the observation operator corresponding to the grid point is generated using the basic spatiotemporal covariance function.

[0024] Specifically, when the first grid point in the spatiotemporal neighborhood... When the nth observation point is a one-dimensional orbital observation point, the grid point and the nth observation point are related. The observation operator between observation points is determined by the fundamental spatiotemporal covariance kernel:

[0025] ;

[0026] ;

[0027] ;

[0028] in, Based on the spatiotemporal covariance kernel, For spatial distance, The longitudinal characteristic scale; For latitudinal characteristic scale, For time scale, For the time difference, Here, L is the spatial attenuation parameter, and L is the comprehensive spatial characteristic scale. It is an exponential function;

[0029] When the first in the spatiotemporal neighborhood of that grid point When the observation point is a two-dimensional interferometric observation point, the grid point and the... The observation operator between observation points is jointly determined by the fundamental spatiotemporal covariance kernel and the two-dimensional interferometric observation direction correlation covariance term:

[0030] ;

[0031] ;

[0032] ;

[0033] in, For the horizontal covariance, For the covariance along the track, The horizontal distance. For horizontal correlation scale, For the relevant dimensions along the track, This represents the difference in the index along the track.

[0034] Furthermore, the observation error covariance matrix is ​​a block covariance matrix constructed based on the observation point type, wherein the diagonal elements are the effective observation error variance of each observation point in the spatiotemporal neighborhood of the grid point, and the off-diagonal elements are zero or the correlation error term between any two observation points in the spatiotemporal neighborhood of the grid point.

[0035] Specifically, when two observation points have the same observation point type and are located within the same orbital period, the off-diagonal element is the correlation error term between the two observation points; when two observation points are located within the same orbital period but have different observation point types, the off-diagonal element is zero; when two observation points are located within different orbital periods, the off-diagonal element is zero.

[0036] Furthermore, based on the observation point type, the effective observation error variance includes a one-dimensional effective observation error variance and a two-dimensional effective observation error variance:

[0037] When the observation point is a one-dimensional observation point along the track, the variance of the one-dimensional effective observation error is calculated using the following formula:

[0038] ;

[0039] in, The variance of the one-dimensional effective observation error. It is the instrument's random noise variance. The residual long-wavelength error variance;

[0040] When the observation point is a two-dimensional interferometric observation point, the variance of the two-dimensional effective observation error is calculated using the following formula:

[0041] ;

[0042] ;

[0043] in, The effective observation error variance in two dimensions. This is the dynamic variance inflation factor related to lateral distance. The width is half the width of the cut. Let α be the horizontal distance. scale The variance inflation factor is a scale-dependent factor. It is the basic variance inflation coefficient.

[0044] Furthermore, when both observation points are one-dimensional along-track observation points, the correlation error term between the two observation points is calculated using the residual long-wavelength error; when both observation points are two-dimensional interferometric observation points, the correlation error term between the two observation points is calculated by combining the long-wavelength correlation term with the Roll error along-track exponential decay term through joint modeling.

[0045] Specifically, the Roll error along the track exponential decay term is: ,in This is the Roll error variance magnitude parameter, used to characterize the baseline strength of the Roll-related error terms; The interval between two two-dimensional interferometric observation points on the same track is... The length scale of the correlation along the Roll error is used to characterize how quickly the correlation of the Roll error along the track direction decays.

[0046] Another aspect of the present invention provides a joint estimation system for sea level anomalies from multi-dimensional multi-source satellite altimeters, used to implement the aforementioned joint estimation method for sea level anomalies from multi-dimensional multi-source satellite altimeters, comprising:

[0047] The observation data acquisition module is used to acquire one-dimensional track-along observation data and two-dimensional interferometric observation data;

[0048] The state space construction module is used to reconstruct one-dimensional along-track observation data and two-dimensional interferometric observation data into a unified observation structure and map it to a regular spatiotemporal grid. The observation structure includes observation type identifier, lateral distance identifier, and along-track index identifier.

[0049] The observation operator construction module is used to construct observation operators corresponding to one-dimensional along-track observation data and two-dimensional interferometric observation data, respectively.

[0050] The error covariance modeling module is used to construct the observation error covariance matrix containing off-diagonal correlation terms for one-dimensional track-side observation data; at the same time, based on the lateral distance-related dynamic variance inflation factor and the scale-dependent variance inflation factor, it constructs the observation error covariance matrix containing off-diagonal correlation terms for two-dimensional interferometric observation data.

[0051] The joint estimation module is used to construct a joint statistical estimation objective function without explicit systematic error correction, thereby completing the optimal estimation of sea level anomalies and outputting a regular gridded product.

[0052] The beneficial effects of adopting the above technical solution are as follows:

[0053] (1) By constructing a unified joint statistical objective function, this invention incorporates satellite altimeter observation data of different spatial dimensions into the same statistical estimation framework for processing, enabling one-dimensional along-orbit observation data and two-dimensional interferometric observation data to participate in sea level anomaly estimation under unified constraints, thus avoiding the problem of fragmented observation information in the phased processing method.

[0054] (2) This invention characterizes the error correlation of two-dimensional interferometric observation data in different spatial directions through the covariance term of two-dimensional interferometric observation direction, so that the effective information at small and medium scales can be reasonably utilized, while suppressing the estimation bias caused by error anisotropy and improving the analytical capability of sea level anomaly structure at small and medium scales.

[0055] (3) The present invention adaptively adjusts the statistical weights of different observation data by using a dynamic variance inflation factor, thereby avoiding the problem of local information dominance caused by the sampling density of two-dimensional interferometric observation data being significantly higher than that of one-dimensional along-track observation data, thus effectively preventing the imbalance of joint estimation results in spatial scale.

[0056] (4) Without performing explicit error correction, this invention suppresses systematic errors in the joint estimation process through reasonable statistical structure design and error covariance modeling, reduces the dependence on empirical error correction models, and improves the robustness of sea level anomaly estimation. Attached Figure Description

[0057] Figure 1 This is a flowchart of the estimation method of the present invention;

[0058] Figure 2 This is a schematic diagram illustrating the downsampling principle of the present invention.

[0059] Figure 3 This is a schematic diagram of the state space-observation space mapping of the present invention;

[0060] Figure 4 This is a schematic diagram illustrating the calculation principle of the observation operator for two-dimensional interferometric observation points in this invention.

[0061] Figure 5 This is a partially enlarged cross-sectional view of the high-kinetic-energy sea area in the Northwest Pacific obtained by the present invention;

[0062] Figure 6This is a comparison chart of the power spectral density of the regular gridded sea level anomaly product of this invention and the original SWOT observation. Detailed Implementation

[0063] The technical solution of the present invention will be described in detail below with reference to the accompanying drawings.

[0064] Example 1:

[0065] refer to Figure 1 A cross-dimensional, multi-source satellite altimeter joint estimation method for sea level anomalies, following steps S1 to S5, performs joint modeling and estimation of one-dimensional along-orbit observation data and two-dimensional interferometric observation data within a unified state space and a unified statistical estimation framework, generating a regular gridded sea level anomaly product for the target sea area over the target time period:

[0066] Step S1 involves acquiring seven consecutive days of one-dimensional along-orbit observation data globally using Sentinel-6A, Sentinel-3A / B, SARAL, and HY-2B satellites, each equipped with an altimeter and orbiting along its own orbit, and forming a one-dimensional along-orbit observation dataset. Specifically, Sentinel-6A orbits at an altitude of approximately 1336 km with a revisit period of approximately 10 days; Sentinel-3A / B orbits at an altitude of approximately 814 km with a revisit period of approximately 27 days; SARAL orbits at an altitude of approximately 800 km; and HY-2B orbits at an altitude of approximately 971 km with a revisit period of approximately 14 days.

[0067] In this embodiment, the acquired one-dimensional along-track observation data are all high-precision point observations. That is, each satellite performs a measurement at its fixed measurement cycle directly below it (the nadir point), forming a one-dimensional along-track observation point to provide long-wave reference information for sea level anomalies. Each one-dimensional along-track observation point corresponds to one set of one-dimensional along-track observation data; in other words, the one-dimensional along-track observation data is a set of attributes of the one-dimensional along-track observation point. Specifically, the one-dimensional along-track observation data includes: observation time, observation location, sea level anomaly value, mass indicator, and orbital indicator. The observation location is determined by the geographic longitude and geographic latitude of the one-dimensional along-track observation point.

[0068] Simultaneously, the SWOT satellite, equipped with a KaRIn interferometric altimeter, acquires seven consecutive days of two-dimensional interferometric observation data globally, forming a two-dimensional interferometric observation dataset. In this embodiment, the acquired two-dimensional interferometric observation data covers a wide swath area on both sides of the satellite orbit, using regularly arranged pixels as the basic unit; that is, each pixel is a two-dimensional interferometric observation point. Specifically, the two-dimensional interferometric observation data includes observation time, observation location, sea level anomaly value, mass marker, and orbital marker, and the observation location is determined by the geographic longitude and geographic latitude of the two-dimensional interferometric observation point.

[0069] Furthermore, the one-dimensional orbital observation dataset is subjected to quality screening and geophysical corrections as follows to form a standard one-dimensional orbital observation dataset:

[0070] Step S1.1: Based on the quality label corresponding to each one-dimensional track observation point, and according to the quality label parameters of the altimeter product and the preset sea level anomaly threshold, the one-dimensional track observation dataset is divided into a valid one-dimensional track observation data subset and an invalid one-dimensional track observation data subset. The invalid one-dimensional track observation data subset is then deleted, that is, the one-dimensional track observation points affected by instrument malfunction, track malfunction or environmental interference are removed.

[0071] Specifically, the quality labeling parameters include:

[0072] (1) ice_flag: Ice cover flag, a value of 0 indicates no ice cover;

[0073] (2) surface_type: surface type, a value of 0 indicates open ocean;

[0074] (3) orb_state_flag: Orb state flag, a value of 0 indicates that the orbit is normal;

[0075] (4) Sea level anomaly threshold: |SLA| ≤ 2.0 m.

[0076] Step S1.2 involves applying standard geophysical corrections to each one-dimensional orbital observation data point in the effective one-dimensional orbital observation data subset, and then normalizing the corrected one-dimensional orbital observation data using the DTU21 mean sea level model to generate standard sea level anomalies for each one-dimensional orbital observation point, thus forming a standard one-dimensional orbital observation dataset. The standard geophysical corrections include: orbital correction, solid tide correction, load tide correction, polar tide correction, atmospheric backpressure correction, and sea state bias correction.

[0077] At the same time, refer to Figure 2 The two-dimensional interferometric observation dataset is spatially aggregated according to the following steps to form a downsampled two-dimensional interferometric observation dataset:

[0078] Step S1.1: Divide the target sea area using a preset downsampling factor to form non-overlapping and non-interval spatial aggregation windows, and determine the subset of two-dimensional interferometric observation data corresponding to each spatial aggregation window based on the observation position of each two-dimensional interferometric observation point.

[0079] Step S1.2: For each spatial aggregation window, based on the preset quality identification parameters and preset physical parameter thresholds, select each effective two-dimensional interferometric observation point contained in the spatial aggregation window, and determine whether the number of effective two-dimensional interferometric observation points in the spatial aggregation window is greater than the preset threshold. If yes, proceed to step S1.3; otherwise, discard the spatial aggregation window and return to step S1.2 until all spatial aggregation windows have been judged.

[0080] Step S1.3: For each spatial aggregation window, based on its corresponding subset of two-dimensional interferometric observation data, an arithmetic mean is calculated on the sea level anomalies corresponding to all valid two-dimensional interferometric observation points within the spatial aggregation window to generate a super observation point. Each super observation point corresponds to a super observation data set, which includes the center observation time, center observation location, center sea level anomaly value, lateral distance identifier, and orbital index identifier. In this embodiment, the super observation point retains the areal features and geometric attributes (lateral distance, orbital index) of the two-dimensional interferometric observation point, and its essence remains that of a two-dimensional interferometric observation point; the super observation data is two-dimensional interferometric observation data after downsampling processing, and its essence also remains that of two-dimensional interferometric observation data.

[0081] Step S1.4: Based on the track index identifier of each super observation point, the sea level anomaly values ​​of all super observation points within the same track are processed to remove the mean, forming a downsampled two-dimensional interferometric observation dataset.

[0082] In this embodiment, the downsampling parameters are as follows:

[0083] Table 1. Downsampling parameters

[0084] Step S2, refer to Figure 3 The standard one-dimensional orbital observation dataset and the downsampled two-dimensional interferometric observation dataset are mapped to a preset spatial frame, and then a corrected joint observation dataset is formed by sequentially using reconstruction and deviation correction operations.

[0085] In this embodiment, by mapping the standard one-dimensional along-orbit observation dataset and the downsampled two-dimensional interferometric observation dataset to a preset spatial framework, observation data from different satellites, orbits, and dimensions can be unified to the same spatiotemporal reference. Through reconstruction operations, the data structures of one-dimensional along-orbit observation data and two-dimensional interferometric observation data can be unified, and one-dimensional along-orbit observation points and two-dimensional interferometric observation points can be distinguished by observation type identifiers. At the same time, the lateral distance identifier and along-orbit index identifier representing the two-dimensional interferometric observation points are displayed.

[0086] Furthermore, the deviation correction operation specifically includes:

[0087] Step S2.1: Based on the track identifier of each observation point, determine the observation points contained in each track, that is, determine the track to which each observation point belongs; and based on the observation position of each observation point, connect the observation points on the same track to construct a spatial linear geometric object, and each spatial linear geometric object has a spatial index.

[0088] Step S2.2: Based on the spatial index of each spatial linear geometric object, use geometric calculation to determine the orbital intersection points between different orbits or between different orbital segments of the same orbit, and calculate the sea level anomaly value of each orbit at each orbital intersection point through neighborhood linear interpolation, that is, perform weighted calculation based on the sea level anomaly value of the nearest observation point on both sides of the orbital intersection point and its corresponding distance.

[0089] Step S2.3: Select Sentinel-6A as the reference satellite, construct a star-shaped intersection topology centered on it, and select the orbit intersections containing the Sentinel-6A orbit as relevant intersections, that is, only retain the orbit intersections between each orbit and the Sentinel-6A orbit.

[0090] Step S2.4: Based on the sea level anomaly values ​​of each orbit at each relevant intersection point and the sea level anomaly values ​​of the Sentinel-6A operating orbit at each relevant intersection point, the deviation values ​​of each orbit and the Sentinel-6A operating orbit at all relevant intersection points are obtained. Then, the deviation prediction value corresponding to each orbit is generated by using a multinomial fitting method with the distance along the orbit as the independent variable. Finally, by subtracting the deviation prediction value of the sea level anomaly value corresponding to the observation point from the deviation prediction value of the orbit, a corrected joint observation dataset with consistent system deviation is obtained.

[0091] Step S3, refer to Figure 3 Based on preset temporal and spatial resolutions, a regular spatiotemporal grid corresponding to the target sea area is constructed, and the calibrated joint observation dataset is mapped to the regular spatiotemporal grid. Furthermore, according to a preset spatiotemporal search radius, the spatiotemporal neighborhood of each grid point is determined, and each observation point is assigned to its corresponding spatiotemporal neighborhood based on its observation location. This determines the subset of calibrated joint observation data corresponding to each grid point. In other words, the spatiotemporal neighborhood of each grid point may contain one-dimensional along-orbit observation points from different satellites or orbits, and may also contain super-observation points.

[0092] In this embodiment, the spatial resolution of the regular spatiotemporal grid is 0.125°×0.125°, and the temporal resolution is 7 days.

[0093] Step S4: For each grid point, based on its corresponding subset of calibrated joint observation data, construct the observation operator C, background covariance matrix F, and observation error covariance matrix R corresponding to that grid point, and combine the background covariance matrix and the observation error covariance matrix into a joint covariance matrix A.

[0094] Step S5: For each grid point, based on its corresponding observation operator and joint covariance matrix, construct the joint statistical estimation objective function corresponding to the grid point without explicit systematic error correction, and calculate the optimal estimate of the sea level anomaly of the grid point by means of the optimal interpolation method. Then, based on the optimal estimate of the sea level anomaly of each grid point, generate a global sea level anomaly product with a spatial resolution of 0.125° and a time window of 7 days.

[0095] Preferably, in this embodiment, the joint covariance matrix A is obtained by adding the background covariance matrix F and the observation error covariance matrix R, and the weight vector W is obtained by solving the linear equation A·W=C. Thus, the construction and solution of the joint statistical estimation objective function are achieved without explicit systematic error correction, and the solution process implicitly implements the following constraints:

[0096] (1) One-dimensional along-track observation data has a high weight in large-scale sea level anomaly estimation due to its small error variance and stable long-wave accuracy;

[0097] (2) Two-dimensional interferometric observation data play a role in the analysis of small and medium-scale structures due to the direction-correlation covariance term and variance inflation factor, while suppressing the unreasonable dominance of the large-scale background.

[0098] Furthermore, this embodiment utilizes a spatiotemporal covariance function to generate observation operators corresponding to each grid point. The spatiotemporal covariance function includes a basic spatiotemporal covariance function and a dual-scale spatiotemporal covariance function, and the basic spatiotemporal covariance function includes a basic spatiotemporal covariance kernel and a two-dimensional interferometric observation direction-related covariance term.

[0099] For each grid point, when the number of observation points in the spatiotemporal neighborhood of that grid point exceeds a preset threshold, the following dual-scale spatiotemporal covariance function is used to generate the correlation between that grid point and the first... Observation operators between observation points:

[0100] ;

[0101] in, For grid points and the first Observation operators between observation points For large-scale weights, For small-scale weights, As a large-scale component, it is used to characterize the spatial-temporal correlation structure of the medium-to-large-scale background field, ensuring the smoothness and stability of the estimated field; As a small-scale component, it is used to characterize the spatial-temporal correlation structure of small-to-medium-scale structures (such as vortices and fronts), enhancing the ability to resolve details.

[0102] Preferably, the large-scale component is adjusted according to the longitudinal and latitudinal characteristic scales; the small-scale component is a fixed short scale, which is set to 80km in this embodiment.

[0103] For each grid point, when the number of observation points in the spatiotemporal neighborhood of the grid point is not greater than the preset threshold for the number of observation points, the observation operator corresponding to the grid point is generated using the basic spatiotemporal covariance function.

[0104] Specifically, when the first grid point in the spatiotemporal neighborhood... When the nth observation point is a one-dimensional orbital observation point, the grid point and the nth observation point are related. The observation operator between observation points is determined by the fundamental spatiotemporal covariance kernel:

[0105] ;

[0106] ;

[0107] ;

[0108] ;

[0109] ;

[0110] ;

[0111] in, Based on the spatiotemporal covariance kernel, For spatial distance, The longitudinal characteristic scale; For latitudinal characteristic scale, For time scale, For time difference; L is the spatial attenuation parameter, and L is the comprehensive spatial characteristic scale. The geographic latitude of the target grid point is used to determine the meridional characteristic scale, latitudinal characteristic scale, and time scale. It is an exponential function;

[0112] refer to Figure 4 When the first grid point in the spatiotemporal neighborhood When the observation point is a two-dimensional interferometric observation point, the grid point and the... The observation operator between observation points is jointly determined by the fundamental spatiotemporal covariance kernel and the two-dimensional interferometric observation direction correlation covariance term:

[0113] ;

[0114] ;

[0115] ;

[0116] in, For the horizontal covariance, For the covariance along the track, The horizontal distance. For horizontal correlation scale, For the relevant dimensions along the track, This is the difference in the index along the track. In this embodiment, =30000m, =400m.

[0117] Furthermore, the observation error covariance matrix is ​​a block covariance matrix constructed based on the observation point type, wherein the diagonal elements are the effective observation error variance of each observation point in the spatiotemporal neighborhood of the grid point, and the off-diagonal elements are zero or the correlation error term between any two observation points in the spatiotemporal neighborhood of the grid point.

[0118] Specifically, when two observation points have the same observation point type and are located within the same orbital period, the off-diagonal element is the correlation error term between the two observation points; when two observation points are located within the same orbital period but have different observation point types, the off-diagonal element is zero; when two observation points are located within different orbital periods, the off-diagonal element is zero.

[0119] Furthermore, when both observation points are one-dimensional along-track observation points, the correlation error term between the two observation points is calculated using the residual long-wavelength error; when both observation points are two-dimensional interferometric observation points, the correlation error term between the two observation points is calculated by combining the long-wavelength correlation term with the Roll error along-track exponential decay term through joint modeling.

[0120] Preferably, the Roll error exhibits exponential decay along the trajectory as follows: ,in This is the Roll error variance magnitude parameter, used to characterize the baseline strength of the Roll-related error terms; The interval between two two-dimensional interferometric observation points on the same track is... This is the length scale for the correlation of Roll error along the track, used to characterize the rate of decay of the correlation of Roll error along the track direction. In this embodiment... , =500m.

[0121] Based on the observation point type, the effective observation error variance includes a one-dimensional effective observation error variance and a two-dimensional effective observation error variance:

[0122] When the observation point is a one-dimensional observation point along the track, the variance of the one-dimensional effective observation error is calculated using the following formula:

[0123] ;

[0124] in, The variance of the one-dimensional effective observation error. It is the instrument's random noise variance. This represents the residual long-wavelength error variance. In this embodiment, , .

[0125] When the observation point is a two-dimensional interferometric observation point, the variance of the two-dimensional effective observation error is calculated using the following formula:

[0126] ;

[0127] ;

[0128] in, The effective observation error variance in two dimensions. This is the dynamic variance inflation factor related to lateral distance. The width is half the width of the cut. Let α be the horizontal distance. scale The variance inflation factor is a scale-dependent factor that is continuously adjusted based on the local observation density, and observations in sparse regions are further deweighted. The basic variance inflation coefficient is used to compensate for the advantages of two-dimensional interferometric observation data in terms of degrees of freedom, thereby achieving a balance between long-wavelength accuracy and small-scale resolution in the estimation results. In this embodiment, =2.0, W swath =60km.

[0129] Furthermore, the dynamic variance inflation factor and the scale-dependent variance inflation factor are not external empirical weights independent of the observation error covariance matrix. Instead, they serve as constituent parameters of the observation error covariance matrix, participating in its construction to determine the effective observation error variance and adjust the strength of related terms, thus statistically characterizing and suppressing systematic errors. Simultaneously, by utilizing the variance inflation factor, the statistical freedom advantage brought by the high sampling density of two-dimensional interferometric observation data can be suppressed without negating its high-resolution advantage, preventing it from exerting non-physical dominance on the joint estimation results in local regions.

[0130] This embodiment, through joint modeling of effective observation error variance and orbital correlation error, maintains the small-scale analytical capability of two-dimensional interferometric observation data while suppressing the statistical degree of freedom advantage caused by its high-density sampling and the influence of orbital systematic errors, effectively improving the balance between large-scale stability and small-scale resolution of the joint estimation results.

[0131] Figure 5 This is a magnified profile of a high-kinetic-energy sea area in the Northwest Pacific (142°E-160°E, 25°N-40°N) obtained using the joint estimation method of this embodiment. Figure 5 It can be observed that the regular gridded sea level anomaly products obtained by the joint estimation method in this embodiment do not have orbital seams, that is, the obtained regular gridded sea level anomaly products are continuous fields, not discrete blocks "pasted by orbit"; and there are no abrupt changes at the boundaries between different satellites or different orbital segments, and the splicing boundary is significantly weakened, which indicates that the weight allocation of different observation sources is more reasonable under a unified statistical framework; at the same time, the joint estimation method in this embodiment can clearly and accurately reconstruct the mesoscale vortex and front structure in terms of spatial morphology.

[0132] The technical effects of this embodiment will be further explained in detail below based on the aforementioned global sea level anomaly product and comparative experiments.

[0133] In this embodiment, the comparative experiment was implemented based on the following hardware and software environments:

[0134] Hardware environment: Multi-core CPU server (≥20 cores), memory ≥64GB, supporting parallel computing;

[0135] Operating system: Microsoft Windows 11;

[0136] Software environment: Python 3.x runtime environment, which depends on scientific computing libraries such as NumPy, SciPy, and Numba, as well as the NetCDF4 data interface library.

[0137] Using September 2025 as the experimental period, the baseline data was the SWOT KaRIn L3 level 2D Swath sea surface height product (variable name ssha_unfiltered) released by CNES / CLS. This baseline data covers 196 independently sampled orbit files worldwide and has not undergone any spatial filtering, thus preserving the most authentic sub-mesoscale physical characteristics and load noise floor.

[0138] During the verification process, a trilinear interpolation algorithm was used to interpolate the global sea level anomaly products obtained using the joint estimation method of this implementation to the original sampling geographic coordinates of the aforementioned 196 SWOT orbits. Furthermore, taking the key band at a small to medium scale of 50.0 km to 200.0 km as an example, the Welch power spectral density (PSD) estimation method was used to calculate the along-track power spectrum of the original SWOT observations and the fused estimation profiles, respectively, to achieve quantitative verification of the signal preservation capability. The experimental statistical results are as follows: Figure 6 As shown.

[0139] from Figure 6 It can be observed that, within the 50.0 km to 200.0 km band, the energy retention rate of the global sea level anomaly product obtained using the joint estimation method of this embodiment reaches 0.8113. Furthermore, without explicit systematic error correction, the global sea level anomaly product completely overlaps with the original SWOT observations in the long-wavelength band above 100 km, demonstrating that the direction-dependent observation error covariance model can effectively absorb and suppress systematic biases in two-dimensional interferometric observation data.

[0140] Example 2:

[0141] This embodiment provides a cross-dimensional multi-source satellite altimeter sea level anomaly joint estimation system to implement the cross-dimensional multi-source satellite altimeter sea level anomaly joint estimation method described in Embodiment 1, including:

[0142] The observation data acquisition module is used to acquire one-dimensional track-along observation data and two-dimensional interferometric observation data;

[0143] The state space construction module is used to reconstruct one-dimensional along-track observation data and two-dimensional interferometric observation data into a unified observation structure and map it to a regular spatiotemporal grid. The observation structure includes observation type identifier, lateral distance identifier, and along-track index identifier.

[0144] The observation operator construction module is used to construct observation operators corresponding to one-dimensional along-track observation data and two-dimensional interferometric observation data, respectively.

[0145] The error covariance modeling module is used to construct the observation error covariance matrix containing off-diagonal correlation terms for one-dimensional track-side observation data; at the same time, based on the lateral distance-related dynamic variance inflation factor and the scale-dependent variance inflation factor, it constructs the observation error covariance matrix containing off-diagonal correlation terms for two-dimensional interferometric observation data.

[0146] The joint estimation module is used to construct a joint statistical estimation objective function without explicit systematic error correction, thereby completing the optimal estimation of sea level anomalies and outputting a regular gridded product.

[0147] The above description is merely a preferred embodiment of the present invention and does not constitute any limitation on the present invention. Any simple modifications, alterations, or equivalent structural changes made to the above embodiments based on the technical essence of the present invention shall still fall within the protection scope of the present invention.

Claims

1. A cross-dimensional, multi-source satellite altimeter joint estimation method for sea level anomalies, characterized in that, Following steps S1 to S5, within a unified state space and a unified statistical estimation framework, jointly model and estimate one-dimensional along-track observation data and two-dimensional interferometric observation data to generate a regular gridded sea level anomaly product for the target sea area over the target time period: Step S1: Using a preset number of target satellites equipped with different altimeters, observation data corresponding to each observation point in the target sea area within the target time period are obtained respectively, and a one-dimensional along-orbit observation dataset and a two-dimensional interferometric observation dataset are formed; by applying quality screening and geophysical correction to the one-dimensional along-orbit observation dataset, a standard one-dimensional along-orbit observation dataset is formed, and at the same time, the two-dimensional interferometric observation dataset is spatially aggregated to form a downsampled two-dimensional interferometric observation dataset. Step S2: Map the standard one-dimensional track-along observation dataset and the downsampled two-dimensional interferometric observation dataset to a preset spatial frame, and then use reconstruction and bias correction operations in sequence to form a corrected joint observation dataset. Step S3: Based on the preset temporal resolution and preset spatial resolution, construct a regular spatiotemporal grid corresponding to the target sea area, and map the calibrated joint observation dataset to the regular spatiotemporal grid. Further, based on the preset spatiotemporal search radius, determine the spatiotemporal neighborhood of each grid point in the regular spatiotemporal grid, as well as the subset of calibrated joint observation data corresponding to each grid point. Step S4: For each grid point, based on its corresponding subset of calibrated joint observation data, construct the observation operator, background covariance matrix, and observation error covariance matrix corresponding to that grid point, and combine the background covariance matrix and the observation error covariance matrix into a joint covariance matrix. Step S5: For each grid point, based on its corresponding observation operator and joint covariance matrix, construct the joint statistical estimation objective function corresponding to the grid point without explicit systematic error correction, and calculate the optimal estimate of the sea level anomaly of the grid point by the optimal interpolation method. Then, based on the optimal estimate of the sea level anomaly of each grid point, generate a regular gridded sea level anomaly product of the target sea area for the target time period.

2. The joint estimation method for sea level anomalies from multi-dimensional multi-source satellite altimeters according to claim 1, characterized in that, The observation points mentioned in step S1 include one-dimensional along-track observation points and two-dimensional interferometric observation points. The one-dimensional along-track observation points correspond to one-dimensional along-track observation data, and the two-dimensional interferometric observation points correspond to two-dimensional interferometric observation data. Both the one-dimensional along-track observation data and the two-dimensional interferometric observation data include observation time, observation location, sea level anomaly value, mass marker, and orbit marker.

3. The joint estimation method for sea level anomalies from multi-dimensional multi-source satellite altimeters according to claim 2, characterized in that, Step S1 applies quality screening and geophysical correction to the one-dimensional orbital observation dataset as follows, to form a standard one-dimensional orbital observation dataset: Step S1.1: Based on the preset quality identification parameters and preset physical parameter thresholds, the one-dimensional track observation dataset is divided into a valid one-dimensional track observation data subset and an invalid one-dimensional track observation data subset, and the invalid one-dimensional track observation data subset is deleted. Step S1.2: Apply standard geophysical corrections to each one-dimensional orbital observation data in the effective one-dimensional orbital observation data subset, and normalize each one-dimensional orbital observation data after applying standard geophysical corrections using a preset mean sea level model, thereby forming a standard one-dimensional orbital observation dataset.

4. The joint estimation method for sea level anomalies from multi-dimensional multi-source satellite altimeters according to claim 2, characterized in that, Step S1 involves spatial aggregation of the two-dimensional interferometric observation dataset to form a downsampled two-dimensional interferometric observation dataset as follows: Step S1.1: Divide the target sea area using a preset downsampling factor to form non-overlapping and non-interval spatial aggregation windows, and determine the two-dimensional interferometric observation data subset corresponding to each spatial aggregation window; Step S1.2: For each spatial aggregation window, based on the preset quality identifier parameter and the preset physical parameter threshold, select each effective two-dimensional interferometric observation point contained in the spatial aggregation window, and determine whether the number of effective two-dimensional interferometric observation points in the spatial aggregation window is greater than the preset threshold. If yes, proceed to step S1.3; otherwise, discard the spatial aggregation window and return to step S1.2, until all spatial aggregation windows have been judged. Step S1.3: For each spatial aggregation window, based on its corresponding subset of two-dimensional interferometric observation data, the super observation point corresponding to the spatial aggregation window is generated by taking the arithmetic mean of the sea level anomaly values ​​corresponding to all valid two-dimensional interferometric observation points within the spatial aggregation window. Each super observation point corresponds to a super observation data, and each super observation data includes the center observation time, center observation location, center sea level anomaly value, lateral distance indicator, and track index indicator. Step S1.4: Based on the orbital index identifier of each super observation point, the mean value of the central sea level anomaly of all super observation points in the same orbit is processed to form a downsampled two-dimensional interferometric observation dataset.

5. The joint estimation method for sea level anomalies from multi-dimensional multi-source satellite altimeters according to claim 2, characterized in that, Each target satellite operates along its corresponding orbit, and step S2 performs deviation correction as follows: Step S2.1: Based on the track identifier of each observation point, determine the observation points contained in each track, and based on the observation position of each observation point, construct each observation point on the same track into a spatial linear geometric object, and each spatial linear geometric object has a spatial index. Step S2.2: Based on the spatial index of each spatial linear geometric object, use geometric calculation to determine the orbital intersection points between different orbits or between different orbital segments of the same orbit, and calculate the sea level anomaly value of each orbit at each orbital intersection point through neighborhood linear interpolation. Step S2.3: With the target reference satellite orbiting along the target reference orbit as the center, construct a star-shaped intersection topology and select the orbit intersections containing the target reference orbit as relevant intersections, that is, only retain the orbit intersections between each orbit and the target reference orbit. Step S2.4: Based on the sea level anomaly values ​​of each orbit at each relevant intersection point and the sea level anomaly values ​​of the target reference orbit at each relevant intersection point, the deviation values ​​of each orbit and the target reference orbit at all relevant intersection points are obtained. Then, the deviation prediction value corresponding to each orbit is generated by using a multinomial fitting method with the distance along the orbit as the independent variable. Finally, by subtracting the deviation prediction value of the sea level anomaly value corresponding to the observation point from the deviation prediction value of the orbit it belongs to, a corrected joint observation dataset with consistent system deviation is obtained.

6. The joint estimation method for sea level anomalies from multi-dimensional multi-source satellite altimeters according to claim 2, characterized in that, Step S4 generates observation operators corresponding to each grid point using a spatiotemporal covariance function, wherein the spatiotemporal covariance function includes a basic spatiotemporal covariance function and a dual-scale spatiotemporal covariance function, and the basic spatiotemporal covariance function includes a basic spatiotemporal covariance kernel and a two-dimensional interferometric observation direction-related covariance term. For each grid point, when the number of observation points in the spatiotemporal neighborhood of that grid point exceeds a preset threshold, a dual-scale spatiotemporal covariance function is used to generate the correlation between that grid point and the first grid point. Observation operators between observation points: ； in, For grid points and the first Observation operators between observation points For large-scale weights, For small-scale weights, For large-scale components, Small-scale components; For each grid point, when the number of observation points in the spatiotemporal neighborhood of the grid point is not greater than the preset threshold for the number of observation points, the observation operator corresponding to the grid point is generated using the basic spatiotemporal covariance function. Specifically, when the th grid point in the spatiotemporal neighborhood... When the nth observation point is a one-dimensional orbital observation point, the grid point is related to the nth observation point. The observation operator between observation points is determined by the fundamental spatiotemporal covariance kernel: ； ； ； in, Based on the spatiotemporal covariance kernel, For spatial distance, The longitudinal characteristic scale; For latitudinal characteristic scale, For time scale, For the time difference, Here, L is the spatial attenuation parameter, and L is the comprehensive spatial characteristic scale. It is an exponential function; When the first in the spatiotemporal neighborhood of that grid point When the observation point is a two-dimensional interferometric observation point, the grid point and the... The observation operator between observation points is jointly determined by the fundamental spatiotemporal covariance kernel and the two-dimensional interferometric observation direction correlation covariance term: ； ； ； in, For the horizontal covariance, For the covariance along the track, The horizontal distance. For horizontal correlation scale, For the relevant dimensions along the track, This represents the difference in the index along the track.

7. The joint estimation method for sea level anomalies from multi-dimensional multi-source satellite altimeters according to claim 2, characterized in that, The observation error covariance matrix is ​​a block covariance matrix constructed based on the observation point type, where the diagonal elements are the effective observation error variance of each observation point in the spatiotemporal neighborhood of the grid point, and the off-diagonal elements are zero or the correlation error term between any two observation points in the spatiotemporal neighborhood of the grid point. Specifically, when two observation points have the same observation point type and are located within the same orbital period, the off-diagonal element is the correlation error term between the two observation points; when two observation points are located within the same orbital period but have different observation point types, the off-diagonal element is zero; when two observation points are located within different orbital periods, the off-diagonal element is zero.

8. The joint estimation method for sea level anomalies from multi-dimensional multi-source satellite altimeters according to claim 7, characterized in that, Based on the observation point type, the effective observation error variance includes a one-dimensional effective observation error variance and a two-dimensional effective observation error variance: When the observation point is a one-dimensional observation point along the track, the variance of the one-dimensional effective observation error is calculated using the following formula: ； in, The variance of the one-dimensional effective observation error. It is the instrument's random noise variance. The residual long-wavelength error variance; When the observation point is a two-dimensional interferometric observation point, the variance of the two-dimensional effective observation error is calculated using the following formula: ； ； in, The effective observation error variance in two dimensions. This is the dynamic variance inflation factor related to lateral distance. The width is half the width of the cut. Let α be the horizontal distance. scale The variance inflation factor is a scale-dependent factor. It is the basic variance inflation coefficient.

9. The cross-dimensional multi-source satellite altimeter sea level anomaly joint estimation method according to claim 7, characterized in that, When both observation points are one-dimensional along-track observation points, the correlation error term between the two observation points is calculated using the residual long-wavelength error; when both observation points are two-dimensional interferometric observation points, the correlation error term between the two observation points is calculated by combining the long-wavelength correlation term with the Roll error along-track exponential decay term through joint modeling. Specifically, the Roll error along the track exponential decay term is: ,in This is the Roll error variance magnitude parameter, used to characterize the baseline strength of the Roll-related error terms; The interval between two two-dimensional interferometric observation points on the same track is... The length scale of the correlation along the Roll error is used to characterize how quickly the correlation of the Roll error along the track direction decays.

10. A cross-dimensional multi-source satellite altimeter sea level anomaly joint estimation system, used to implement the cross-dimensional multi-source satellite altimeter sea level anomaly joint estimation method according to any one of claims 1-9, characterized in that, include: The observation data acquisition module is used to acquire one-dimensional track-along observation data and two-dimensional interferometric observation data; The state space construction module is used to reconstruct one-dimensional along-track observation data and two-dimensional interferometric observation data into a unified observation structure and map it to a regular spatiotemporal grid. The observation structure includes observation type identifier, lateral distance identifier, and along-track index identifier. The observation operator construction module is used to construct observation operators corresponding to one-dimensional along-track observation data and two-dimensional interferometric observation data, respectively. The error covariance modeling module is used to construct the observation error covariance matrix containing off-diagonal correlation terms for one-dimensional track-side observation data; at the same time, based on the lateral distance-related dynamic variance inflation factor and the scale-dependent variance inflation factor, it constructs the observation error covariance matrix containing off-diagonal correlation terms for two-dimensional interferometric observation data. The joint estimation module is used to construct a joint statistical estimation objective function without explicit systematic error correction, thereby completing the optimal estimation of sea level anomalies and outputting a regular gridded product.

Images