A bayesian inversion method for extracting wide-swath altimetry data balanced signals

Abstract

This invention discloses a Bayesian inversion extraction method for the equilibrium signal of wide-span altimeter data. The method involves acquiring and preprocessing sea surface height anomaly sequences from radar interferometers and nadir altimeters. A normalized sine square window function is applied for windowing, and multidimensional spatial averaging is used to estimate the one-dimensional wavenumber power spectrum. A piecewise power law-based equilibrium signal spectrum model and a noise spectrum model constrained by dynamic sea state are constructed, and the set of spectral parameters is extracted through logarithmic domain weighted least squares fitting. A set of spatial covariance matrices is constructed using cosine integral transform and Abelian forward and inverse transforms. A graphics processor is scheduled to perform batch matrix decomposition and singular fault-tolerant regularized inversion to solve for the posterior mean vector and posterior covariance matrix of the target equilibrium signal. Window fusion and index mapping are applied to fill the gaps in nadir observations. Geostrophic dynamics parameters are calculated, uncertainty quantification is performed based on the linear error propagation law, and the knowledge base is updated based on the exponential moving average algorithm. This invention achieves suppression of observation noise and physical filling of observation gaps, improving the adaptability of the inversion system to environmental changes while preserving non-Gaussian dynamic characteristics.

Description

Technical Field

[0001] This invention relates to the field of marine remote sensing and satellite altimetry data processing technology, specifically a Bayesian inversion extraction method for balanced signals from wide-swath altimetry data. Background Technology

[0002] Satellite altimetry is an important means of acquiring global-scale sea surface height observation data. With the evolution of wide-swath radar interferometry, satellites can provide high-resolution sea surface height topographic data with two-dimensional spatial coverage characteristics, providing a foundation for detailed research on ocean dynamic processes.

[0003] In the field of marine science, wide-swath altimetry data is widely used to study mesoscale and sub-mesoscale eddy dynamics. Current data processing workflows typically employ spatial low-pass filtering or conventional optimal interpolation algorithms to suppress noise in the original observation sequence, and then deduce dynamic parameters such as geostrophic velocity and relative vorticity from the reconstructed sea surface height field.

[0004] In practical applications, sea surface height observation data acquired by wide-swath radar interferometers generally contain complex non-stationary red noise. Existing spatial smoothing or conventional filtering methods typically assume that the observed noise is statistically stable white noise, without fully considering the influence of changes in the observed noise level with significant wave height and environmental sea conditions such as wind speed. This approach, which ignores the spatiotemporal variability of noise, makes it difficult for the system to accurately identify and remove background noise components in strong wind and wave environments. Often, inaccurate estimation of parameters for complex noise characteristics leads to energy loss in the extracted equilibrium signal, or incomplete noise suppression causes numerical divergence in subsequent dynamic parameter estimation processes. Summary of the Invention

[0005] This invention aims to solve the technical problems of high-frequency measurement noise interference on the balanced ocean dynamics signal in wide-swath satellite altimetry data, and the breakage of the two-dimensional physical space ocean dynamics structure caused by the gap in nadir point observation.

[0006] The first aspect of this invention provides a Bayesian inversion extraction method for balanced signals of wide-amplitude height measurement data, comprising the following processing steps:

[0007] The sea surface height anomaly sequences from radar interferometers and nadir altimeters are acquired, and preprocessing including tidal and benchmark corrections, as well as spatiotemporal mean removal, is performed. The preprocessed sea surface height anomaly sequences are then divided into data segments along the track, and a sinusoidal squared window function satisfying variance normalization constraints is applied for spatial windowing. Complex spectra are obtained through discrete Fourier transform, and the modulus squared values ​​at each discrete wavenumber point are calculated to obtain the power spectral density sequence. A multidimensional spatial arithmetic mean is performed on the power spectral density sequence within the effective observation pixels and observation period to estimate the one-dimensional wavenumber power spectrum of the radar interferometer and the nadir altimeter.

[0008] A parameterized spectral model balancing ocean dynamics signals and space observation noise is constructed. The one-dimensional wavenumber spectrum of the balanced ocean dynamics signal is defined as following a piecewise power-law parameterization form including spectral amplitude, transition wavelength, and spectral slope parameters. A parameterized model of space observation noise, including red and white noise, is constructed, superimposed with a corrected spectral model characterizing spaceborne smoothing and subsampling aliasing effects. The significant wave height and near-surface wind speed parameters attached to the observation data are extracted, and a priori values ​​of dynamic noise amplitude constrained by sea state are constructed as the starting point for fitting. A cutoff interval for the least squares parameter search is defined. Historical spectral parameters are retrieved from a global spectral parameter knowledge base as initial fitting values. The one-dimensional wavenumber power spectrum and the parameterized spectral model are mapped to logarithmic coordinate space, and a weighted least squares fitting objective function is constructed using the reciprocal of the ordinary wavenumber as a weighting factor. When the calculated residual of the objective function exceeds the limit or the noise spectral slope parameter exceeds the physically reasonable range, a spatial distance-inverse weighted backoff mechanism is activated. Parameters from known successfully fitted historical spatial windows are extracted and weighted average substitution is performed to obtain the balanced signal spectral parameter set and the noise spectral parameter set.

[0009] By applying Wiener-Khinchin's theorem, the parameterized spectral model of the equilibrium signal is transformed into a fundamental equilibrium covariance function through cosine integral transform. Abelian forward and inverse Abelian transform operators are used to establish a mapping relationship between the one-dimensional orbital wavenumber spectrum and the two-dimensional isotropic wavenumber spectrum. Combined with spatial physical geometric distance, a set of spatial covariance matrices containing both autocovariance and cross-covariance matrices is constructed.

[0010] A discrete checklist of covariance values ​​corresponding to physical geometric distances is pre-constructed offline and linear interpolation addressing is performed during online computation. A joint observation vector is constructed and cut into multiple overlapping local computation windows along the track direction. Lower triangular matrix decomposition is performed on the covariance matrix of the observation dataset. When a non-positive definite error signal is triggered during decomposition, a basic regularization coefficient and a fault-tolerant iteration round counter are introduced to perform regularization correction and retry. If decomposition still fails after reaching the upper limit of retries, the corresponding computation window is marked as a failure state, and the computation results of the neighboring successful decomposition windows are extracted to perform linear extrapolation as a fallback. The posterior mean vector and posterior covariance matrix of the target equilibrium signal are solved by matrix inversion and product operations.

[0011] Within the spatially overlapping area of ​​adjacent calculation windows, a weighted arithmetic mean is calculated on multiple posterior estimates according to the associated Hanning window function weights. A network index mapping relationship is established between the posterior mean vector and the full-coverage two-dimensional spatial coordinate system, reshaping the one-dimensional posterior mean vector into a two-dimensional balanced sea surface height reconstruction value, thus completing the filling of the nadir observation gaps.

[0012] Coriolis parameters are calculated based on geographic latitude coordinates. A second-order spatial central difference scheme is applied to internal grid nodes, and a one-sided second-order spatial difference scheme is applied to boundary grid nodes to calculate the geostrophic velocity component and the geostrophic relative vorticity scalar field. A pixel-by-pixel division operation is performed on the geostrophic relative vorticity scalar field using the Coriolis parameter matrix to obtain the normalized geostrophic vorticity field. A physical constraint mask is set to mark the geostrophic velocity calculation results as invalid values ​​when the absolute value of the geographic latitude coordinate is less than 5°. Based on linear Taylor expansion and the error propagation law, the uncertainty error distribution matrix of the geostrophic velocity field and the normalized geostrophic vorticity field is calculated based on the absolute uncertainty standard deviation of the posterior covariance matrix.

[0013] Based on the extracted set of legitimate noise spectrum parameters, the exponential moving average algorithm is applied to iteratively update the recorded values ​​in the global spectrum parameter knowledge base. The Pearson correlation coefficient and root mean square error between the two-dimensional balanced sea surface height reconstruction map and the external multi-source fusion satellite altimetry grid product are calculated. When the Pearson correlation coefficient or root mean square error triggers a preset anomaly judgment rule threshold, the state update operation of the global spectrum parameter knowledge base is stopped, and the data corresponding to the abnormal orbital segment is added to the pending review queue.

[0014] A second aspect of the present invention provides a Bayesian inversion extraction system for balanced signals of wide-amplitude height measurement data, comprising:

[0015] The data preprocessing unit is used to acquire the sea surface height anomaly sequence and perform tidal, benchmark correction and spatiotemporal mean removal processing;

[0016] The spectral feature estimation unit is used to estimate the one-dimensional wavenumber power spectrum through spatial windowing, discrete Fourier transform, and multidimensional spatial arithmetic mean operation.

[0017] The parameterized modeling unit is used to construct a spectral model that balances the signal and the observation noise. It constructs a priori values ​​of dynamic noise amplitude based on the significant wave height and near-sea wind speed, extracts the set of spectral parameters through weighted least squares fitting, and applies a spatial distance-inverse weighted backoff mechanism to replace parameters when the residuals exceed the limits.

[0018] The covariance mapping unit is used to convert a parameterized spectral model into a set of spatial covariance matrices through cosine integral transform and Abelian integral transform.

[0019] The parallel inversion solution unit is used to generate matrix elements based on discrete check tables, decompose the spatial covariance matrix, and solve the posterior mean vector and posterior covariance matrix through regularization correction and linear extrapolation fallback mechanism.

[0020] Spatial reconstruction units are used to perform weighted fusion of the calculation results of overlapping windows, and to reconstruct two-dimensional balanced sea surface height values ​​through coordinate mapping to fill the observation gaps.

[0021] The dynamic evaluation unit is used to calculate the geostrophic dynamic parameters and uncertainty error distribution with mask constraints, and to perform global spectral parameter knowledge base updates based on the multi-source cross-validation calculation results.

[0022] This invention provides a Bayesian inversion extraction method for balanced signals from wide-amplitude height measurement data. It has the following beneficial effects:

[0023] 1. This invention constructs a joint noise model that includes satellite-borne smoothing correction, subsampling aliasing correction, and correction based on auxiliary variables of effective wave height and sea surface wind speed. It incorporates the non-stationary red noise characteristics generated by wide-span radar interferometers and data fluctuations caused by sea state into a unified parameterized representation framework. By introducing dynamic noise amplitude prior constraints in the least squares fitting process, it effectively suppresses the parameter estimation divergence problem of traditional smoothing algorithms in strong wind and wave environments. It achieves accurate removal of complex background noise and complete preservation of balanced signal energy in wide-span altimetry data.

[0024] 2. This invention establishes a mathematical mapping relationship between a one-dimensional wavenumber spectrum and a two-dimensional physical field by applying Abelian forward and inverse transforms, and performs joint Gaussian process inference by combining the spatial covariance matrix generated based on the Wiener-Khinchin theorem transformation. It uses effective observation data from radar interferometers on both sides to perform numerical reconstruction of the nadir point gap region under physical constraints. This solves the technical bottleneck of physical distortion and signal discontinuity in the observation fault zone of traditional interpolation methods, and realizes consistent and seamless filling of the sea surface height field across the entire raft width, as well as pixel-by-pixel uncertainty quantification of the reconstruction results.

[0025] 3. This invention retains the piecewise power-law decay characteristics that conform to the laws of ocean dynamics in wavenumber spectrum modeling, and introduces an online update mechanism for the global spectral parameter knowledge base based on the exponential moving average algorithm and a multi-source cross-validation feedback alarm mechanism. While dynamically tracking changes in ocean physical state, it effectively avoids the excessive smoothing of asymmetric dynamic structures by conventional spatial low-pass filtering through Bayesian posterior mean solution. This achieves high-fidelity extraction of the non-Gaussian statistical distribution characteristics of mesoscale and sub-mesoscale ocean processes, ensuring that the output geostrophic velocity field and normalized geostrophic vorticity field have extremely high dynamic reliability. Attached Figure Description

[0026] Figure 1 This is a flowchart illustrating the overall process of the method of the present invention.

[0027] Figure 2 This is a diagram showing the signal-to-noise separation results of the SWOT test data of the present invention;

[0028] Figure 3 This is a graph showing the wavenumber spectrum fitting results of the present invention;

[0029] Figure 4 This is a diagram illustrating the vorticity statistics and uncertainty analysis of the present invention. Detailed Implementation

[0030] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0031] Reference Figure 1 , Figure 1 The present invention provides a Bayesian inversion extraction method for balanced signals from wide-swath altimetry data. This method comprises five logical stages: data preprocessing, spectral analysis and modeling, Bayesian inversion and reconstruction, dynamic quantity calculation and output, and experimental verification. The overall workflow is based on two-dimensional sea surface height anomaly data observed by a spaceborne Ka-band radar interferometer. Balanced ocean dynamic signals are separated through mathematical inversion, filling the gaps in nadir point observations.

[0032] The data preprocessing stage comprises steps S1 to S3. Step S1 involves data acquisition and quality screening, reading Level-2 data products containing Ka-band radar interferometer low-rate sea surface height data and nadir altimeter data. Based on the quality marker fields provided in the data products, contaminated or abnormally interfered observation points are removed, as are ripples with missing data rates and spatial variance anomalies. The filtered observation dataset is then output. Step S2 involves tidal and benchmark correction, restoring high-resolution empirical intratidal corrections to the observation dataset and deducting geoid error components that have undergone time averaging and spatial filtering. Step S3 involves spatiotemporal mean removal, subtracting the spatial distribution mean of ripple data within the same observation period from the observation dataset, and detrending the observation dataset along the satellite flight trajectory. The preprocessed spatial sequence data of sea surface height anomalies is then output and input into the spectral analysis and modeling stage.

[0033] Step S1 specifically includes the following sub-steps:

[0034] Step S101: Read the raw altimeter data. Receive the wide-swath satellite altimeter data product from the satellite data center as the basic input. The wide-swath satellite altimeter data product is a Level-2 data product issued by the Surface Water and Ocean Topography mission. Analyze the wide-swath satellite altimeter data product to extract the radar interferometer low-rate sea surface height data and nadir altimeter observation data contained within it. The radar interferometer low-rate sea surface height data provides two-dimensional swath measurement information with a spatial resolution of 2km. The nadir altimeter observation data provides one-dimensional nadir point measurement information with a spatial resolution of 7km. For each spatial observation point, extract the sea surface height anomalies associated with the spatial observation point and the quality marker field used to indicate the data status.

[0035] Step S102: Perform pixel-level spatial quality screening. Traverse all resolved spatial observation points and perform data status judgment based on the quality marker field. Observation points with non-zero quality marker values ​​are removed. A non-zero quality marker indicates that the spatial observation point is affected by rain signal attenuation, radar instrument hardware malfunction, or external electromagnetic measurement interference. Observation points with non-zero quality marker values ​​are determined as invalid observation points and cleared from the data memory, retaining only valid spatial observation points with a quality marker of zero.

[0036] Step S103: Perform overall consistency screening at the raft level. After completing the pixel-level screening, structural integrity and numerical stability are assessed using a single spatial raft as an independent statistical unit. Calculate the proportion of observation points removed due to quality defects within each independent raft to the total number of observation points that should exist in the independent raft, and obtain the data missing rate. If the data missing rate is greater than 20%, the independent raft is identified as a region with severe spatial structural fractures and is removed entirely. Simultaneously calculate the total spatial variance of sea surface height anomalies within each retained raft. Summarize the sea surface height anomalies of all retained rafts within the same batch of tracks and calculate the overall sea surface height anomaly variance as a baseline level. Compare the total spatial variance of a single raft with the baseline level of the overall sea surface height anomaly variance. When the total spatial variance of a single raft is greater than 10 times the baseline level, the single raft is identified as having a systematic anomaly, and the single raft is removed entirely.

[0037] Step S104: Construct the effective spatial observation dataset. All spatial observation points remaining after pixel-level spatial quality cleaning and swath-level overall removal are compiled and combined to construct the effective observation dataset. The effective observation dataset consists of discrete observation column vectors, specifically including the effective observation vectors from the radar interferometer. and the effective observation vector of the nadir altimeter The effective observation dataset contains records of the radar interferometer's effective observation vectors. and the effective observation vector of the nadir altimeter The latitude and longitude spatial coordinates corresponding to each sea surface height anomaly scalar value within the data support subsequent spatial covariance matrix derivation and inversion extraction of equilibrium ocean dynamics signals.

[0038] Step S2 specifically includes the following sub-steps:

[0039] Step S201: Recover the high-resolution empirical tidal internal tide correction. In the system's default processing flow, the sea surface height anomaly values ​​have already been pre-subtracted by the high-resolution empirical tidal internal tide correction. Read the high-resolution empirical tidal internal tide correction values ​​bound to the valid observation dataset. Add the read high-resolution empirical tidal internal tide correction values ​​back to the sea surface height anomaly values ​​according to the spatial observation points. Perform the recovery operation to avoid accidentally deleting the low wavenumber balanced ocean dynamics signal variance in the subsequent signal-to-noise separation stage.

[0040] Step S202: Calculate the time-averaged sea surface height field. Collect sea surface height observation data from approximately 90 revisits at the same spatial location during the rapid repeating orbit phase. Perform an arithmetic mean operation on the sea surface height observation data from the 90 revisits along the time dimension, based on the coordinates of the spatial observation points, to generate a time-averaged sea surface height field covering the entire observation area.

[0041] Step S203 involves spatial filtering and subtraction of geoid error. A Gaussian high-pass filter with a cutoff scale of 100 km is applied to process the time-averaged sea level height field to extract the geoid error signal components. Specifically, the Gaussian high-pass filter is implemented by subtracting the convolution result of the time-averaged sea level height field and a two-dimensional Gaussian smoothing kernel from the time-averaged sea level height field.

[0042] After acquiring the geoid error signal, subtract each geoid error signal from the sea surface height anomaly values ​​recovered after internal tide correction, according to their corresponding spatial coordinates. Output the spatial sequence data of sea surface height anomalies after tidal and benchmark correction, and input it into the subsequent spatiotemporal mean removal step.

[0043] Step S3 specifically includes the following sub-steps:

[0044] Step S301: Perform large-scale spatial mean removal. Obtain the spatial sequence data of sea surface height anomalies after tidal and baseline correction. Calculate the spatial arithmetic mean of the radar interferometer swath data within a single observation period. The spatial arithmetic mean is calculated by summing the sea surface height anomaly values ​​of all valid radar interferometer observation pixels within a single observation period, and then dividing the sum by the total number of valid observation pixels involved in the summation.

[0045] Step S302: Perform detrending processing along the track direction. For the sea surface height anomaly sequence after deducting the spatial arithmetic mean, apply detrending processing along the track direction to remove any potentially residual large-scale gradients. Output the preprocessed sea surface height anomaly sequence. as well as Preprocessed sea surface height anomaly sequence as well as Input the data into the spectral analysis and modeling module as the basic input data for wavenumber spectrum estimation.

[0046] The spectral analysis and modeling stage comprises steps S4 to S6. Step S4 involves wavenumber spectrum estimation, applying windowing and discrete Fourier transform to the preprocessed sea surface height anomaly spatial sequence data, and calculating the wavenumber power spectrum along the track. Step S5 involves fitting the equilibrium signal spectrum model, using a weighted least squares method to fit a piecewise power-law parameterized equilibrium signal spectrum model based on the along-track wavenumber power spectrum, determining the spectral amplitude, transition wavelength, and spectral slope parameters of the equilibrium signal spectrum model. Step S6 involves fitting the noise spectrum model, constructing red noise and white noise models respectively, introducing significant wave height and near-sea surface wind speed parameters to perform adaptive noise amplitude correction to constrain the fitting search interval, and triggering a spatial distance inverse weighting mechanism to backtrack to historical parameters when the fitting residual exceeds the limit. The output is a set of basic spectral parameters containing the statistical characteristics of the equilibrium signal and noise, which is then used in the Bayesian inversion and reconstruction stage.

[0047] Step S4 specifically includes the following sub-steps:

[0048] Step S401: Perform along-orbit data segmentation and windowing. Receive the pre-processed radar interferometer sea surface height anomaly sequence and nadir altimeter sea surface height anomaly sequence. Divide the radar interferometer and nadir altimeter sea surface height anomaly sequences into continuous along-orbit data segments of fixed length along the satellite's flight orbit. Set the physical length of each along-orbit data segment to 790 km. Apply a sinusoidal squared window function to perform spatial windowing on each along-orbit data segment.

[0049] Step S402: Perform Discrete Fourier Transform and Complex Spectrum Calculation. Perform Discrete Fourier Transform on the windowed sea surface height anomaly data segment to obtain the complex spectrum. The Discrete Fourier Transform maps the wave signal sequence in the spatial distance domain to a complex sequence in the spatial frequency domain, outputting the complex real and complex imaginary parts at each discrete spatial wavenumber position according to the ordinary wavenumber interval.

[0050] Step S403: Perform power spectral density estimation and multidimensional spatial averaging. Calculate the squared modulus of the complex spectrum at each discrete wavenumber point. Apply a discrete Fourier transform length scaling factor to the squared modulus values ​​for numerical normalization to obtain the power spectral density sequence. For the data generated by radar interferometer observations, traverse all cross-track pixels and all observation periods corresponding to the radar interferometer data, and perform an arithmetic mean operation on the calculated power spectral density sequence to obtain the smoothed one-dimensional wavenumber power spectrum of the radar interferometer. For data generated by a nadir altimeter, an arithmetic mean was performed on the power spectral density sequence over all effective observation periods to obtain a smoothed one-dimensional wavenumber power spectrum from the nadir altimeter. .

[0051] One-dimensional wavenumber power spectrum of radar interferometer One-dimensional wavenumber power spectrum of nadir altimeter All using ordinary wavenumber The ordinary wavenumber is an independent variable. The unit is period per kilometer, ordinary wavenumber The value of is the reciprocal of the value of the wavelength in physical space. Estimate the output one-dimensional wavenumber power spectrum of the radar interferometer. One-dimensional wavenumber power spectrum of nadir altimeter All satisfy the global variance normalization constraint, and the mathematical formula for spectral normalization is expressed as:

[0052]

[0053] In the formula, This represents the sea surface height anomaly sequence values ​​input before the windowing process. This represents the statistical variance corresponding to the numerical values ​​of the sea surface height anomaly sequence. For statistical expectation operators, Represents ordinary wavenumber The corresponding power spectral density value at that location. This indicates a continuous integration operation performed within the positive half-wavenumber interval. This represents the wavenumber integral infinitesimal element. One-dimensional wavenumber power spectrum of a radar interferometer. One-dimensional wavenumber power spectrum of nadir altimeter Transmitted to the equilibrium signal spectrum model fitting module.

[0054] Step S5 specifically includes the following sub-steps:

[0055] Step S501: Construct a parameterized spectral model for the equilibrium ocean dynamics signal. Receive the one-dimensional wavenumber power spectrum of the radar interferometer output from the spectral analysis and modeling module. Assume the one-dimensional wavenumber spectrum of the equilibrium ocean dynamics signal follows a piecewise power-law parameterized form, and construct a parameterized spectral model for the equilibrium ocean dynamics signal.

[0056] The mathematical expression for the parameterized spectral model is:

[0057]

[0058] In the formula, This represents the theoretical estimate of the power spectral density of the balanced signal at ordinary wavenumber k. Spectral amplitude parameter (unit) ), indicating the height of the low wavenumber platform; The transition wavelength parameter (in km) represents the characteristic scale at which the spectrum transitions from a plateau to a power law, physically corresponding to the size of the dominant mesoscale vortex. is the spectral slope parameter (dimensionless), representing the power-law exponent in the high wavenumber region. In the low wavenumber region... )hour, It presents a plateau; at high wavenumbers ( )hour, It exhibits power-law decay. A weighted least squares method based on log-log space is used for fitting, with weighting factors of [value missing]. This prevents densely sampled high-wavenumber data from dominating the fitting results. A typical fitting result is as follows: , , The output is the balanced signal spectrum parameters. .

[0059] Step S502: Define the objective function for logarithmic domain weighted least squares fitting. Then, use the one-dimensional wavenumber power spectrum of the radar interferometer and the parameterized spectral model. Mapping to logarithmic coordinate space for nonlinear fitting. Introducing ordinary wavenumber. reciprocal The weighted least squares fitting objective function is constructed using the weighting factors, based on the ordinary wavenumber. reciprocal The weighting factor constrains the dominant role of densely sampled high-wavenumber observation frequency band data in the global fitting results.

[0060] Step S503: Solve for the parameters of the balanced signal spectrum model. Execute a nonlinear optimization algorithm to fit the objective function using weighted least squares. Minimize. Extract the weighted least squares fitting objective function. The set of equilibrium signal spectrum parameters when the minimum value is reached.

[0061] The balanced signal spectrum parameter set includes the balanced signal spectrum amplitude parameter. Balanced signal transition wavelength parameters and the equilibrium signal spectrum slope parameter The extracted set of balanced signal spectrum parameters is transmitted to the noise spectrum model fitting module as input conditions for calculating the signal-to-noise separation transition wavelength.

[0062] Step S6 specifically includes the following sub-steps:

[0063] Step S601: Construct a parameterized model for space observation noise. Construct a red noise model for the radar interferometer and a white noise model for the nadir altimeter, respectively.

[0064] The noise spectrum of the radar interferometer is in Matérn form, and its mathematical expression is as follows:

[0065]

[0066] In the formula, Noise spectral amplitude (unit) ); The noise transition wavelength is fixed at 100km because it has little impact on the fitting. This represents the slope of the noise spectrum.

[0067] For data from a nadir altimeter, the noise is white noise, and the spectral model is:

[0068]

[0069] In the formula, Represents ordinary wavenumber Theoretical value of the power spectral density model of the altimeter at the nadir point; This represents the numerical value of the parameterized spectral model of the balanced signal; The spatial sampling interval of the nadir altimeter; This represents the statistical standard deviation of the sea surface height measurement at the sub-satellite point.

[0070] Step S602 involves joint correction of onboard smoothing and aliasing effects. The sea surface height anomaly data transmitted by the radar interferometer is subject to onboard smoothing effects at the satellite hardware level. A two-dimensional Gaussian filter in the frequency domain is constructed to characterize the onboard smoothing effect.

[0071] The response function of a two-dimensional Gaussian filter is expressed as:

[0072]

[0073] In the formula, Indicates the ordinary wavenumber along the track Compared with transorbital ordinary wavenumber The two-dimensional wave vector mode length, This represents the filter scale constant set based on the pixel size of the radar interferometer.

[0074] The red noise model of the radar interferometer is superimposed with the parameterized spectrum model of the balanced signal. The inverse Abelian transform is applied to map the superposition result to a two-dimensional wavenumber space. The result is then multiplied with the response function of a two-dimensional Gaussian filter to achieve smooth attenuation. Finally, the forward Abelian transform is applied to map the attenuation result back to a one-dimensional wavenumber space.

[0075] The joint mathematical expression for the radar interferometer observation spectrum model is:

[0076]

[0077] In the formula, This represents the theoretical value of the radar interferometer observation spectrum model after spaceborne smoothing correction. This represents the Abelian positive transformation operator. This represents the inverse Abelian transform operator.

[0078] A subsampling aliasing compensation term is introduced simultaneously to construct an aliasing-corrected observation spectrum model. The mathematical expression of the aliasing-corrected observation spectrum model is:

[0079]

[0080] In the formula, This is the index for the folding order, with values ​​ranging from 0, 1, to 2. This is the KaRIn observation spectral model after spaceborne smoothing correction. By superimposing the aforementioned third-order folded energies, the spectral energy aliasing caused by subsampling can be effectively corrected, ensuring the accuracy of the spectral model in the high wavenumber region.

[0081] Step S603: Perform adaptive noise amplitude prior correction based on sea state auxiliary variables. Extract the effective wave height auxiliary variable Hs and near-sea surface wind speed auxiliary variable from the same orbit of the wide-swath satellite altimetry data product. An adaptive noise amplitude prior model is constructed to dynamically adapt to time-varying sea conditions.

[0082] The calculation expression for the adaptive noise amplitude prior model is as follows:

[0083]

[0084] In the formula, The prior value of the dynamic noise amplitude is constrained by sea state; Baseline amplitude; The effective wave height scalar value recorded synchronously at the altimeter observation point, in units of The wind speed scalar value at a height of 10m above the sea surface, recorded synchronously at the altimeter observation point, is expressed in units of... The empirical sensitivity coefficient is determined by offline engineering calibration, and typically takes the following values: .

[0085] Prior value of dynamic noise amplitude This is set as the initial iteration starting point for the nonlinear fitting optimization of the radar interferometer's red noise model. Based on the prior value of the dynamic noise amplitude... The absolute value is used to set the upper and lower limits of the least squares parameter search range with a fluctuation ratio of ±30%, which limits the parameter search space of the least squares optimization algorithm and prevents parameter estimation divergence caused by strong winds and waves.

[0086] Step S604: Solve for the final spectral parameters and execute the fitting anomaly tolerance backoff mechanism. Within the set upper and lower limit cutoff intervals, the system first queries the historical spectral parameters of the corresponding geographic grid from the global spectral parameter knowledge base as initial fitting values. Then, it substitutes the observed data into the aliasing correction observation spectral model and the nadir point white noise model to perform fitting calculations.

[0087] Noise spectrum amplitude parameters were extracted by fitting data from radar interferometer spectral observations. Noise spectrum slope parameter And the set of equilibrium signal spectrum parameters. With the equilibrium signal spectrum parameter set fixed, the statistical standard deviation of nadir measurements is extracted from the altimeter's spectral observation data. .

[0088] Determine whether the fitting status of a single track-side window of data triggers the anomaly tolerance condition. The set anomaly tolerance conditions include a calculated residual value greater than 0.2 for the weighted least squares fitting objective function and the slope parameter of the noise spectrum in the fitting output. It is not within the physically reasonable range of 1.0 to 2.5.

[0089] When any of the fault tolerance conditions is triggered, the fitting results of the window along which the anomaly occurred are discarded, and the spatial distance inverse weighted backoff mechanism is activated. The mathematical calculation logic of the spatial distance inverse weighted backoff mechanism is as follows: extract the known successfully fitted results around the window along which the anomaly occurred. The set of valid parameters for each historical spatial window is used to construct a weighted average based on the reciprocal of the physical distance between the center of the historical spatial window and the center of the abnormal window.

[0090] After the parameter substitution is completed, the output contains the final noise spectrum parameter set containing the valid parameters and the corresponding parameter quality validity flag.

[0091] The Bayesian inversion and reconstruction stage comprises steps S7 to S9. Step S7 involves constructing the covariance matrix. Based on the fundamental spectral parameter set, an integral transform is used to generate a spatial covariance matrix set encompassing the autocorrelation of measured data, the cross-correlation between nadir points and ripple observation points, and the spatial correlation between target reconstruction grid points and observation points. Step S8 performs the Bayesian inversion solution. Based on the spatial covariance matrix set and the actual observation dataset, under the Gaussian process prior assumption, the graphics processor is scheduled to perform block inversion calculations along a sliding overlapping window to solve for the posterior mean vector providing the optimal estimate of the equilibrium signal and the posterior covariance matrix quantifying the estimation uncertainty. Step S9 performs nadir point gap filling and reconstruction. The one-dimensional posterior mean vector is mapped to a two-dimensional grid array using spatial latitude and longitude coordinates. Spatial covariance constraints are used to achieve linear unbiased filling of the physical observation gaps, outputting a two-dimensional equilibrium sea surface height reconstruction map for the dynamics calculation and output stage.

[0092] Step S7 specifically includes the following sub-steps:

[0093] Step S701: Determine the fundamental mapping relationship of the spatial covariance function. Assume that both the equilibrium ocean dynamics signal and the high-frequency observation noise follow a Gaussian distribution, and that they are statistically independent. Based on the Wiener-Khinchin theorem, the spectral model in the frequency domain is converted into a covariance function in the physical space domain.

[0094] The mathematical expression for the cosine integral transform of the spatial covariance function is defined as follows:

[0095]

[0096] In the formula, Represents the corresponding power spectral density model The numerical value of the spatial covariance function; This represents a spatial physical distance variable, with the unit of measurement being km. Represents the ordinary wavenumber as the independent variable; This represents the cosine basis integral term in the inverse Fourier transform.

[0097] Based on the defined cosine integral transform expression, the parameterized spectral model of the equilibrium signal is transformed into the fundamental equilibrium covariance function. .

[0098] Step S702: Construct the autocovariance matrix and cross-covariance matrix of the observation dataset. Calculate the autocovariance matrix between the observation points of the nadir altimeter. The mathematical expression for the elements of the autocovariance matrix of the nadir altimeter observation points is:

[0099]

[0100] in For point and points The distance between them For the Kronecker function (when The value is 1 if the condition is met, and 0 otherwise. KaRIn-KaRIn covariance (considering onboard smoothing):

[0101]

[0102] KaRIn-nadir cross-covariance:

[0103]

[0104] For the target location (reconstructed grid point), the target-target covariance is: The target-nadir covariance is The target-KaRIn covariance is The calculation formulas are similar. In numerical implementation, the Abel transform uses piecewise linear interpolation integration, and the cosine transform uses Discrete Cosine Transform Type I (DCT-I). The wavenumber grid ranges from 0 to... ,spacing Typical , The output is a complete set of covariance matrices. .

[0105] In the formula, The first element in the autocovariance matrix of the radar interferometer is... The observation location of the first cut area and the first The covariance values ​​between the observation locations of each ripple; Indicates the first The observation location of the first cut area and the first The physical geometric distance between the observation locations of each ripple.

[0106] In the underlying computer discretization solution, the Abelian integral transform employs a piecewise linear interpolation integration method to perform discretization calculations to suppress numerical instability caused by the differential term. The cosine integral transform is solved using a discrete cosine transform computation engine.

[0107] The lower limit of the discrete grid interval for the one-dimensional ordinary wavenumber is set to zero, and the upper limit of the discrete grid interval for the ordinary wavenumber is set to... Set the frequency spacing between adjacent ordinary wavenumber grid nodes to . Set constant parameters The value is set to 5000km, and constant parameters are defined. With a value of 100000, based on the assumption that the equilibrium signal and observation noise are independent and both follow a Gaussian distribution, a joint prior distribution of radar interferometer observations, nadir altimeter observations, and target equilibrium signals is constructed:

[0108]

[0109] in, This represents the effective observation vector of the radar interferometer. This is the effective observation vector of the nadir altimeter. Let be the target equilibrium signal vector to be solved. This represents a multivariate normal distribution with a mean of 0, and its covariance matrix is ​​the block matrix described above.

[0110] After completing the numerical integral transformation and matrix element assembly, the set of spatial covariance matrices is output. The set of spatial covariance matrices contains... , , , , as well as The entire input is fed into the Bayesian inversion solution module to support the joint Gaussian inference process.

[0111] Step S8 specifically includes the following sub-steps:

[0112] Step S801: Construct the Gaussian process prior model and joint observation vector. Combine the effective observation vectors from the radar interferometer. Effective observation vector of the nadir altimeter Perform vertical concatenation to generate a joint observation vector. Assume the target equilibrium signal vector to be extracted. With joint observation vector They jointly follow a multivariate Gaussian distribution with a mean of zero.

[0113] Step S802: Offline construction of covariance checklist and online interpolation addressing. Extract discrete numerical pairs from the spatial covariance matrix set, using physical geometric distance as a single independent variable, set the discrete step length parameter of the one-dimensional physical geometric distance, generate a one-dimensional distance grid sequence based on the set discrete step length parameter, calculate the covariance value of the corresponding one-dimensional distance grid sequence, and pre-construct a discrete checklist of the covariance values ​​corresponding to the one-dimensional distance.

[0114] Step S803: Perform track-side sliding overlapping window slice calculation. Set the physical length and window overlap rate of the track-side calculation window, where the physical length of the calculation window is typically set to 200km, and the overlap length of adjacent calculation windows in the track direction is typically set to 40km. Combine the observation vectors... The target reconstruction grid points are cut into multiple overlapping local data slices along the satellite's flight trajectory.

[0115] Within each local data slice, the autocovariance matrix of the observed dataset for that slice is generated using a covariance check table interpolation. The cross-covariance matrix between the target reconstruction point and the observation point and the autocovariance matrix of the target reconstruction point itself. .

[0116] Step S804: Invoke the graphics processor to perform batch Cholesky decomposition and singularity-tolerant adaptive regularization inversion. This involves calculating the autocovariance matrix of the observation dataset from all local data slices. The data is transferred to the graphics processor's video memory. The parallel computing cores of the graphics processor then process the autocovariance matrix of the observation dataset from multiple local slices. Simultaneously perform the lower triangular matrix decomposition operation.

[0117] The equation for the decomposition of a lower triangular matrix is: , in the formula This represents the lower triangular matrix resulting from the decomposition. This represents the transpose of the lower triangular matrix. It determines whether the lower triangular matrix decomposition operation triggers a non-positive definite error signal.

[0118] Upon detecting a non-positive definite error signal, a singular fault-tolerant adaptive regularization mechanism is triggered. The mathematical expression for the singular fault-tolerant adaptive regularization mechanism is:

[0119]

[0120] In the formula, This represents the autocovariance matrix after regularization correction; This represents the base regularization coefficient, set to 10⁻¹. 0 ; This represents a fault-tolerant iteration round counter, which starts at zero and increments by one after each non-positive definite error signal is triggered. Represents the identity matrix.

[0121] The regularized autocovariance matrix is ​​re-inputted into the lower triangular matrix decomposition operation, and the process is retried up to 5 times. If the decomposition still fails after 5 iterations, the calculation window along the track is marked as a failure, and linear extrapolation is performed to provide a numerical fallback by extracting the calculation results of the neighboring successful decomposition windows. After successful decomposition, the posterior mean vector and posterior covariance matrix of the target equilibrium signal are solved.

[0122] The mathematical equation for solving the posterior mean vector is:

[0123]

[0124] The mathematical equation for solving the posterior covariance matrix is:

[0125]

[0126] In the formula, Represents the posterior mean vector of a local slice; This represents the posterior covariance matrix of a local slice.

[0127] Step S805: Apply Hanning weights to perform multi-window overlapping area result fusion. Extract the posterior mean vector of all local data slices. With the posterior covariance matrix For the spatially overlapping regions of adjacent local data slices along the satellite's flight orbit, a Hanning window function is introduced to calculate the fusion weights.

[0128] Within spatially overlapping regions, for multiple posterior mean vector estimates and multiple posterior covariance matrix estimates from different local data slices at the same spatial physical coordinate location, a weighted arithmetic mean is calculated based on the associated Hanning window function weights. The globally optimal posterior mean vector and global posterior covariance matrix are then concatenated and transmitted to the subsequent reconstruction module for gap filling.

[0129] Step S9 specifically includes the following sub-steps:

[0130] Step S901: Construct a full-coverage two-dimensional spatial reconstruction grid. Receive the globally optimal posterior mean vector and global posterior covariance matrix output from the multi-window fusion stage. Read the spatial latitude and longitude range of the original wide-swath altimeter data product. Set the spatial resolution parameters of the two-dimensional spatial reconstruction grid to generate a full-coverage two-dimensional spatial coordinate system covering the left radar interferometer swath observation area, the data-free area in the nadir point observation gap, and the right radar interferometer swath observation area. Extract the total number of cross-track grid nodes and the total number of along-track grid nodes in the full-coverage two-dimensional spatial coordinate system.

[0131] Step S902: Perform coordinate mapping from the one-dimensional posterior vector to the two-dimensional physical space. Establish a position index mapping relationship between the discrete data points in the globally optimal posterior mean vector and the fully covered two-dimensional spatial coordinate system. Extract the scalar values ​​from the globally optimal posterior mean vector and write them one by one into the corresponding grid nodes of the fully covered two-dimensional spatial coordinate system according to the position index mapping relationship.

[0132] Step S903 achieves seamless physical filling of the nadir observation gaps. Within the framework of the joint Gaussian process, the grid node values ​​falling into the nadir observation gap region in the fully covered two-dimensional spatial coordinate system constitute the conditional mathematical expectation inferred based on the observation dataset and spatial covariance matrix constraints.

[0133] The reconstructed two-dimensional equilibrium sea surface height values ​​in the gap region of the nadir observation are generated by performing a linear weighted integral operation on the cross-covariance matrix and the effective observation data of the radar interferometers on both sides, thus realizing the numerical filling of the original observation space fault zone.

[0134] The output is a continuous two-dimensional balanced sea level reconstruction map including the left ripple, the filling gap, and the right ripple.

[0135] Step S904: Extract and construct a pixel-by-pixel uncertainty quantization distribution map. Extract the set of main diagonal elements of the global posterior covariance matrix. The main diagonal elements of the global posterior covariance matrix represent the statistical variance of the optimal estimate of each grid node in the fully covered two-dimensional spatial coordinate system.

[0136] Perform square root operations on the main diagonal elements of the global posterior covariance matrix to calculate the standard deviation of the absolute uncertainty at spatial nodes.

[0137] After calculating the standard deviation of the absolute uncertainty, the standard deviation is reconstructed into a two-dimensional spatial uncertainty distribution map based on a full-coverage two-dimensional spatial coordinate system. The two-dimensional spatial uncertainty distribution map quantifies and marks the reconstructed confidence boundaries of the observation gap region at the nadir and the lateral reaping regions. The continuous two-dimensional equilibrium sea surface height reconstruction map and the two-dimensional spatial uncertainty distribution map are then transferred to the dynamics calculation and output module.

[0138] The dynamics calculation and output stage includes steps S10 to S12. Step S10 calculates the geostrophic velocity by obtaining the spatial partial derivative of the reconstructed two-dimensional equilibrium sea surface height map in the swath coordinate system, and calculating the geostrophic flow velocity components along the track and across the track. Step S11 calculates the geostrophic vorticity by performing spatial differencing on the geostrophic flow velocity components and calculating the curl of the flow field velocity gradient to obtain the geostrophic vorticity scalar field. Step S12 performs uncertainty assessment and output by comparing the mean covariance with the uncertainty posterior covariance to determine the effective resolution wavelength of the spatial filter, and statistically deriving the skewness and kurtosis of the geostrophic vorticity field. Verified spectral parameters are extracted and updated in the locally stored spectral parameter knowledge base using exponential moving averages. The correlation coefficient and root mean square difference between the reconstructed results and independent multi-source satellite altimetry observation data are calculated to trigger an anomaly alarm mechanism. The final output is a digital product set containing the filled sea surface height field, geostrophic flow field parameters, and point-state uncertainty quantification indicators.

[0139] The experimental verification phase was conducted to confirm the effectiveness of the system based on satellite-borne measured sea surface height products. The above-mentioned system processing was sequentially applied to the observation data of the mid-latitude boundary current region. The standard deviation of the equilibrium signal, the standard deviation of the separation noise, and the signal-to-noise ratio parameters were statistically extracted. The uncertainty increment of the gap region was calculated. By comparing the vorticity skewness values ​​before and after processing, it was confirmed that the system completely preserved the non-Gaussian statistical characteristics of ocean dynamics while removing high-frequency noise and reconstructing the fracture gap.

[0140] Step S10 specifically includes the following sub-steps:

[0141] S1001, extracts geostrophic dynamics calculation parameters. Receives the continuous two-dimensional equilibrium sea surface height reconstruction map and two-dimensional spatial uncertainty distribution map output by the gap filling and reconstruction module. Extracts the geographic latitude coordinates of each grid node in the continuous two-dimensional equilibrium sea surface height reconstruction map. In the wiring coordinate system, along the track direction is... The direction of the cross-track is The Earth's rotational velocity component is calculated using the following formula:

[0142]

[0143] in, It is the acceleration due to gravity. Coriolis parameters, The angular velocity of Earth's rotation. For latitude. Spatial derivatives are calculated using a second-order central difference scheme, with interior points using... The boundary points use the corresponding one-sided second-order scheme.

[0144] Calculate the two-dimensional geostrophic velocity field. Based on the geostrophic equilibrium relationship, calculate the geostrophic velocity components in the trans-rail direction and along the rail direction using a continuous two-dimensional equilibrium sea surface height reconstruction map. In the numerical solution, a second-order spatial central difference scheme is applied to replace the continuous spatial partial derivative calculations.

[0145]

[0146] In the formula, This represents the sum of the geostrophic velocity amplitudes. Amplitude summation is performed pixel-by-pixel, outputting a two-dimensional geostrophic velocity amplitude field matrix.

[0147] S1002, since the geostrophic velocity is a linear function of SSH, its uncertainty can be propagated from the posterior covariance through linear error propagation. Calculations show that the velocity uncertainty at the center of the cut is approximately 6.8 cm / s, and approximately 7.9 cm / s at the nadir gap. It should be noted that when the Coriolis parameter... Approaching 0 (i.e., latitude) In low-latitude regions, the Geostrophic relationship is no longer applicable; therefore, the calculated Geostrophic velocity for that region is marked as invalid (NaN) and noted in the output quality label. This method is mainly applicable to mid-to-high latitude sea areas (…). In these regions, the Coriolis parameter is sufficiently large, and the geostrophic equilibrium assumption holds. The output is the geostrophic velocity field. And its uncertainties.

[0148] Step S11 specifically includes the following sub-steps:

[0149] Step S1101: Receive geostrophic dynamics input data. Receive the two-dimensional geostrophic velocity field matrix output by the geostrophic velocity calculation module. This matrix includes the geostrophic velocity component matrix in the trans-orbital direction and the geostrophic velocity component matrix in the along-orbital direction. Simultaneously receive the spatial discrete grid step size in the trans-orbital direction of the reconstructed grid. Spatial discrete grid step size along the track direction , and the Coriolis parameter matrix containing the values ​​of each grid node.

[0150] Step S1102: Calculate the geostrophic relative vorticity scalar field. In fluid mechanics, relative vorticity is defined as the curl of the fluid velocity field gradient, used to characterize the rotational motion of fluid particles. Geostrophic vorticity can be calculated from the spatial derivative of the geostrophic velocity field, or it can be directly solved using the Laplace operator balancing sea surface height.

[0151] The formula for calculating geostrophic vorticity is:

[0152]

[0153] In the formula, Geostrophic vorticity; The component of the geodynamic velocity along the track direction; The trans-rail direction velocity component; It is the acceleration due to gravity; Coriolis parameters It is a two-dimensional Laplace operator; This is a two-dimensional equilibrium sea surface height reconstruction map.

[0154] In the numerical solution program, a second-order spatial central difference scheme is used to calculate the spatial derivative of the two-dimensional geostrophic velocity field matrix. For the boundary grid nodes of the two-dimensional geostrophic velocity field matrix, the corresponding one-sided second-order spatial difference scheme is used to perform differentiation operations to avoid memory access errors. The processing logic for boundary nodes along the track direction is consistent with that in the cross-track direction. Spatial differentiation is performed pixel by pixel to output the geostrophic relative vorticity scalar field matrix.

[0155] Step S1103 Perform geostrophic relative vorticity normalization processing

[0156] Mesoscale and sub-mesoscale ocean dynamics analysis relies on the numerical quantification of the ratio of fluid motion inertial force to Coriolis force. The calculated geostrophic relative vorticity scalar field matrix is ​​divided pixel-by-pixel by the Coriolis parameter matrix to obtain the normalized geostrophic vorticity field. The Coriolis parameter is commonly used for geostrophic vorticity. Normalization is represented as This dimensionless quantity directly reflects the intensity relative to planetary vorticity.

[0157] Output the normalized geostrophic vorticity field matrix.

[0158] Step S1104: Calculate the geostrophic vorticity uncertainty error distribution matrix.

[0159] Calculating geostrophic vorticity from the sea surface height field requires second-order spatial differentiation, which significantly amplifies high-frequency measurement errors. The uncertainty of geostrophic vorticity is derived from the posterior covariance matrix reconstructed from sea surface height using the linear error propagation law. Typical results show that the uncertainty of vorticity at the center of the swath is approximately... , nadir gap approximately Despite some uncertainty, the vorticity can reach [a certain level]. The above-mentioned strong vortex characteristics still significantly exceed the uncertainty level.

[0160] The geostrophic relative vorticity scalar field matrix, the normalized geostrophic vorticity field matrix, and the two-dimensional geostrophic vorticity uncertainty distribution matrix are transmitted to the uncertainty assessment and feedback learning module as the basic data for statistical feature verification.

[0161] Step S12 specifically includes the following sub-steps:

[0162] Step S1201: Perform spatial effective resolution wavelength assessment. Obtain a continuous two-dimensional equilibrium sea surface height reconstruction map and a two-dimensional spatial uncertainty distribution map. Calculate the one-dimensional spatial signal power spectrum of the continuous two-dimensional equilibrium sea surface height reconstruction map along the satellite flight trajectory. Based on the point-by-point posterior variance values ​​contained in the two-dimensional spatial uncertainty distribution map, calculate the corresponding one-dimensional error power spectrum along the satellite flight trajectory. Both the one-dimensional spatial signal power spectrum and the one-dimensional error power spectrum use the ordinary wavenumber as the independent variable. Compare the numerical values ​​of the one-dimensional spatial signal power spectrum and the one-dimensional error power spectrum at the same wavenumber. Locate the first intersection point of the one-dimensional spatial signal power spectrum curve and the one-dimensional error power spectrum curve using a logarithmic spatial linear interpolation algorithm. Define the wavenumber coordinate value corresponding to the first intersection point as the cross wavenumber node. Calculate the reciprocal of the value corresponding to the cross wavenumber node to obtain the effective resolution wavelength parameter of the spatial filtering process.

[0163] Step S1202: Calculate the non-Gaussian statistical characteristic skewness and kurtosis. Extract the scalar values ​​of all valid grid points contained within the normalized geostrophic vorticity field matrix.

[0164] Step S1203: Perform a spectral parameter knowledge base update based on exponential moving average. Extract the set of valid spectral parameters that have not triggered the fitting anomaly tolerance condition. The set of valid spectral parameters includes the equilibrium signal spectral amplitude parameters. Balanced signal transition wavelength parameters Balanced signal spectrum slope parameter Noise spectrum amplitude parameters and noise spectrum slope parameter The system reads numerical records from the global spectral parameter knowledge base, organized by latitude and longitude grid (typically 2°×2°). After each successful inversion, the system iteratively integrates the set of valid spectral parameters into the global spectral parameter knowledge base using an exponential moving average algorithm to perform an online update.

[0165] The mathematical equation for the exponential moving average update algorithm is:

[0166]

[0167] In the formula, This represents the record value written to the global spectral parameter knowledge base after a single iteration update operation; This represents the numerical values ​​corresponding to the set of valid spectral parameters extracted in this operation. This indicates the historical data values ​​extracted and saved before performing update operations within the global spectral parameter knowledge base; The learning rate is typically set to a value of [value to be filled in]. By performing the above iterative update operations, the global spectral parameter knowledge base can be used to track the time-varying ocean physical state online.

[0168] Step S1204: Perform multi-source cross-validation feedback and anomaly alarm freezing. Receive external multi-source fused satellite altimeter grid products independent of the wide-swath satellite altimeter mission system. Extract the reference sea surface height matrix of the external multi-source fused satellite altimeter grid products within the corresponding physical space region. Perform spatial grid node alignment mapping between the continuous two-dimensional balanced sea surface height reconstruction map and the reference sea surface height matrix.

[0169] After completing the cross-validation operator operation, a preset threshold for the feedback anomaly detection rule is set. This threshold includes the Pearson correlation coefficient. Values ​​less than 0.7, and root mean square error The absolute value is greater than 3cm. The system executes Boolean logic judgment. When any of the set feedback anomaly judgment rule thresholds is triggered, an anomaly alarm signal code is generated for the computing system, and the execution of the global spectral parameter knowledge base state update operation based on exponential moving average is suspended to prevent anomaly parameters from contaminating prior knowledge. At the same time, the data corresponding to the anomaly trajectory segment is added to the pending review queue for subsequent investigation. In the state where no anomaly alarm signal is triggered, the continuous two-dimensional equilibrium sea surface height reconstruction map, two-dimensional spatial uncertainty distribution map, geostrophic velocity field matrix, and normalized geostrophic vorticity field matrix are serialized, encapsulated, and output.

[0170] To verify the technical effectiveness of this invention, a verification experiment was conducted using measured data from the SWOT satellite. This experiment was based on the filtering results of the SWOT satellite L3 level product, covering the spectral separation and noise modeling steps in the method of this invention. The experimental data came from the SWOT sea surface height product, selecting observation data from a typical orbit in the mid-latitude western boundary current region.

[0171] Experiment 1: Data Sources and Regional Characteristics

[0172] The experiment used SWOT sea surface height anomaly data. The selected region is a typical mid-latitude region with strong currents, and the SSHA standard deviation reaches 16.45 cm, representing a typical equilibrium ocean dynamic signal. The data patch structure includes a left patch of 30 pixels, a nadir gap of 8 pixels (approximately 14 km), and a right patch of 31 pixels, with a spatial resolution of approximately 2 km.

[0173] Experiment 2: Verification of signal-to-noise separation effect

[0174] Reference Figure 2 , Figure 2 This is a graph showing the signal-to-noise separation results of the SWOT measured data according to the present invention. The signal-to-noise separation effect of the present invention was verified by applying the signal-to-noise separation method of the present invention to the SWOT observation data. Figure 2 (a) Original observations and Figure 2 (b) As shown in the comparison of the extracted balanced signals, the system effectively strips out, as Figure 2 The observation noise is shown in (c). The quantitative statistical results are as follows:

[0175] (1) The standard deviation of the extracted balanced signal is 16.39 cm.

[0176] (2) The standard deviation of the separated noise is 0.63 cm.

[0177] (3) The signal-to-noise ratio (SNR) is 25.92.

[0178] Experiment 3: Wavenumber Spectrum Fitting Verification

[0179] Reference Figure 3 , Figure 3The figure shows the wavenumber spectrum fitting results of this invention. The parametric spectral model of this invention was applied to the observed data for wavenumber spectrum fitting verification. As shown in the logarithmic domain distribution curve, the gray line represents the observed spectrum. The system successfully separated the noise component in the high wavenumber range from the equilibrium component in the low wavenumber range. The specific parameter fitting results are as follows: cm / cpkm, km, Noise spectral model The parameter fitting results are as follows: cm / cpkm, Spectral model goodness of fit .

[0180] Experiment 4: Nadir gap filling and uncertainty verification

[0181] Reference Figure 4 , Figure 4 This is the vortex statistics and uncertainty analysis diagram of the present invention. Part (b) of the diagram shows the cross-track uncertainty distribution profile, where the spatial correlation between the ripple data and the gap region is used to fill the nadir gap. The results show that, according to step S9, the spatial correlation between the ripple data and the gap region is used to fill the nadir gap.

[0182] (1) The uncertainty of the cutting area is about 0.63 cm.

[0183] (2) The uncertainty of the nadir gap region is about 0.75cm.

[0184] (3) The uncertainty at the gap increases by about 19.2%, which is within a reasonable range.

[0185] (4) After filling, the vortex structure remains continuous at the gaps without obvious boundary artifacts.

[0186] Experiment 5: Verification of Geostrophic Dynamics

[0187] Reference Figure 4 The histogram of vorticity probability density distribution in part (a) shows that the skewness of the original observed vorticity field (gray area) is 0.222. After extraction by this system, the skewness of the reconstructed equilibrium vorticity field (blue area) is significantly improved to 0.754, indicating that the non-Gaussian statistical characteristics were successfully preserved. Geostrophic velocity and vorticity were calculated from the reconstructed equilibrium SSH field. The results are as follows:

[0188] (1) The typical flow velocity in the geostrophic velocity field reaches tens of cm / s, and both cyclonic and anticyclonic structures can be identified.

[0189] (2) The original observed vorticity field skewness was 0.222, and the reconstructed equilibrium vorticity field skewness was 0.754.

[0190] Experimental conclusion:

[0191] (1) The standard deviation of the balanced signal is 16.39 cm, the standard deviation of the noise is 0.63 cm, and the SNR is 25.92.

[0192] (2) Goodness of fit of the spectral model It is 0.86.

[0193] (3) The uncertainty of the reaming area is 0.63 cm, the uncertainty of the nadir gap area is 0.75 cm, and the gap area is relatively increased by 19.2%.

[0194] (4) The original observed vorticity skewness was 0.222, and the reconstructed equilibrium vorticity skewness was 0.754.

Claims

1. A Bayesian inversion extraction method for balanced signals of wide-amplitude height measurement data, characterized in that, Includes the following steps: The radar interferometer sea surface height anomaly sequence and the nadir altimeter sea surface height anomaly sequence are acquired, and tidal and reference correction and spatiotemporal mean removal preprocessing are performed. The preprocessed sea surface height anomaly sequence is divided into data segments along the track and windowed. The one-dimensional wavenumber power spectrum of the radar interferometer and the one-dimensional wavenumber power spectrum of the nadir altimeter are estimated by discrete Fourier transform and multidimensional spatial averaging. A parameterized spectral model balancing ocean dynamics signals and space observation noise is constructed. Based on logarithmic domain weighted least squares fitting and adaptive correction of sea state auxiliary variables, the set of balanced signal spectral parameters and the set of noise spectral parameters are extracted. By applying the Wiener-Khinchin theorem and the Abelian integral transform, the parameterized spectral model is converted into a set of spatial covariance matrices; Construct a joint observation vector, perform matrix decomposition and singular fault-tolerant adaptive regularization inversion on the set of spatial covariance matrices, and solve for the posterior mean vector and posterior covariance matrix of the target equilibrium signal; The Hanning window function weight is applied to perform multi-window overlap area result fusion, and the coordinate mapping from one-dimensional posterior vector to two-dimensional physical space is performed to complete the filling of the gaps in the nadir point observation and the full-coverage two-dimensional spatial reconstruction. Geostrophic parameters were calculated and uncertainties were quantified based on the continuous two-dimensional equilibrium sea surface height reconstruction map.

2. The Bayesian inversion extraction method for balanced signals of wide-amplitude height measurement data according to claim 1, characterized in that, The specific method for estimating the one-dimensional wavenumber power spectrum is as follows: Set the physical length of the data segment along the track, and apply a sinusoidal squared window function that satisfies the global variance normalization constraint to perform spatial windowing operation; The complex spectrum is obtained by using discrete Fourier transform, and the power spectral density sequence is obtained by calculating the squared modulus at each discrete wavenumber point. An arithmetic mean operation is performed on the power spectral density sequence across all valid observation pixels and observation periods to obtain the one-dimensional wavenumber power spectrum of the radar interferometer and the one-dimensional wavenumber power spectrum of the nadir altimeter.

3. The Bayesian inversion extraction method for balanced signals of wide-amplitude height measurement data according to claim 1, characterized in that, The specific method for constructing a parameterized spectral model of the equilibrium ocean dynamics signal and extracting parameters is as follows: The one-dimensional wavenumber spectrum of the equilibrium ocean dynamics signal is set to follow a piecewise power-law parameterization form, which includes the equilibrium signal spectrum amplitude parameter, the equilibrium signal turn wavelength parameter, and the equilibrium signal spectrum slope parameter. The one-dimensional wavenumber power spectrum and the parameterized spectral model are mapped to logarithmic coordinate space, and the reciprocal of the ordinary wavenumber is introduced as a weighting factor to construct a weighted least squares fitting objective function to perform parameter fitting.

4. The Bayesian inversion extraction method for balanced signals of wide-amplitude height measurement data according to claim 3, characterized in that, The specific method for constructing the parameterized spectral model of space observation noise is as follows: Construct a parameterized model of space observation noise that includes a red noise model following the Matrn distribution and a nadir white noise model; By introducing a two-dimensional Gaussian filter response function and a folding order summation operator, a corrected observation spectrum model is constructed to superimpose and characterize the onboard smoothing effect and the subsampling aliasing effect.

5. The Bayesian inversion extraction method for balanced signals of wide-amplitude height measurement data according to claim 1, characterized in that, The specific method for extracting the noise spectrum parameter set based on adaptive correction of sea state auxiliary variables is as follows: Extract the effective wave height auxiliary variable and the near-sea surface wind speed auxiliary variable within the same orbit, construct the dynamic noise amplitude prior value constrained by sea state, and use the prior value as the fitting starting point to set the cutoff interval for least squares parameter search. The historical spectral parameters of the corresponding geographic grid are retrieved from the global spectral parameter knowledge base as initial values ​​for fitting, and then the fitting operation is performed. The fitting status of the data along the track window is determined. If the calculated residual of the weighted least square fitting objective function is greater than the set threshold, or the noise spectrum slope parameter exceeds the physically reasonable range, the spatial distance inverse weighted backoff mechanism is activated, and the parameters of the historically successfully fitted spatial windows in the surrounding area are extracted and weighted average replacement is performed.

6. The Bayesian inversion extraction method for balanced signals of wide-amplitude height measurement data according to claim 1, characterized in that, The specific method for constructing the set of spatial covariance matrices is as follows: Based on the Wiener-Khinchin theorem, the parameterized spectral model of the equilibrium signal is transformed into the basic equilibrium covariance function through cosine integral transform. A mapping between one-dimensional orbital wavenumber spectrum and two-dimensional isotropic wavenumber spectrum is established by applying the Abelian forward and inverse Abelian transform operators. Combined with spatial physical geometric distance, a set is constructed that includes the autocovariance matrix of the radar interferometer, the autocovariance matrix of the nadir altimeter, the cross-covariance matrix, the autocovariance matrix of the target reconstruction grid points, and the cross-covariance matrix between the target reconstruction grid points and the observation points.

7. The Bayesian inversion extraction method for balanced signals of wide-amplitude height measurement data according to claim 1, characterized in that, The specific method for performing matrix decomposition and singular fault-tolerant adaptive regularization inversion is as follows: A discrete lookup table of one-dimensional physical geometric distance corresponding to the covariance value is pre-constructed offline, and linear interpolation addressing is performed in the online calculation process to obtain the covariance value; Perform lower triangular matrix decomposition on the autocovariance matrix of the observation dataset; if the decomposition triggers a non-positive definite error signal, then introduce the basic regularization coefficient and the fault-tolerant iteration round counter to perform regularization correction and retry up to 5 times; If the decomposition still fails after 5 iterations, the calculation window along the track is marked as a failure, and the numerical fallback is performed by extracting the calculation results of the neighboring successful decomposition windows and performing linear extrapolation.

8. The Bayesian inversion extraction method for balanced signals of wide-span height measurement data according to claim 1, characterized in that, The specific method for applying the Hanning window function weights to perform multi-window overlap region result fusion and coordinate mapping is as follows: The joint observation vector is cut into multiple overlapping local data slices for calculation windows along the track direction, and a weighted arithmetic mean is calculated for multiple posterior estimates at the same spatial physical coordinate position within the spatially overlapping region, according to the associated Hanning window function weight values. Establish a grid index mapping relationship between the globally optimal posterior mean vector and the full-coverage two-dimensional spatial coordinate system, and reshape the one-dimensional posterior mean vector into a two-dimensional balanced sea surface height reconstruction value.

9. The Bayesian inversion extraction method for balanced signals of wide-amplitude height measurement data according to claim 1, characterized in that, The specific method for performing geostrophic dynamics parameter calculations and uncertainty quantification assessments is as follows: Coriolis parameters are calculated based on geographic latitude coordinates. A second-order spatial central difference scheme is applied to the internal grid nodes, and a one-sided second-order spatial difference scheme is applied to the boundary grid nodes. The geostrophic velocity components and the geostrophic relative vorticity scalar field are calculated. The normalized geostrophic vorticity field is obtained by performing a pixel-by-pixel division operation using the Coriolis parameter matrix. Introducing a physical constraint mask for judgment: When the absolute value of the geographic latitude coordinate is less than 5°, the corresponding geostrophic velocity calculation result will be forcibly marked as invalid. Based on the linear Taylor expansion and the error propagation law, the uncertainty error distribution matrix of the geostrophic velocity field and the normalized geostrophic vorticity field is calculated based on the absolute uncertainty standard deviation of the posterior covariance matrix.

10. The Bayesian inversion extraction method for balanced signals of wide-amplitude height measurement data according to claim 1, characterized in that, Following the steps of calculating geostrophic dynamic parameters and quantifying uncertainties, the process also includes updating the global spectral parameter knowledge base based on a feedback mechanism. Specifically: Based on the extracted set of legal noise spectrum parameters, the exponential moving average algorithm is applied to iteratively update the stored global spectrum parameter knowledge base records. A multi-source cross-validation feedback mechanism is introduced to calculate the Pearson correlation coefficient and root mean square error between the continuous two-dimensional balanced sea surface height reconstruction map and the external multi-source fused satellite altimetry grid product. When the Pearson correlation coefficient or root mean square error triggers the preset anomaly judgment rule threshold, the state update operation for the global spectral parameter knowledge base is stopped, and the data corresponding to the abnormal orbital segment is added to the pending review queue.

Images