A marine gravity comprehensive correction method and system based on multi-source auxiliary data

Abstract

This invention discloses a method and system for comprehensive ocean gravity correction based on multi-source auxiliary data. The method includes acquiring multi-source auxiliary data and gravity observation sequences and aligning them spatiotemporally; obtaining a gravity-corrected observation sequence after motion correction; dividing the data into an irregular grid and segmenting it spatiotemporally; calculating the grid value variation coefficient and drift coefficient; correcting the gravity value sequence; calculating the tidal correction based on third auxiliary data; obtaining the first and second standard gravity values; fusing them to obtain the standard gravity field; calculating the long-wave and short-wave noise contribution indices and constructing anomaly factors; identifying anomalous grid nodes; obtaining the long-wave corrected gravity component; calculating the short-wave signal-to-noise ratio; performing anomaly detection and correcting short waves; outputting the short-wave corrected component; and fusing the long-wave and short-wave corrected components to obtain the gravity correction value. This method provides more efficient and accurate gravity correction results, while also exhibiting good interpretability, and can be directly applied to comprehensive ocean gravity correction systems.

Description

Technical Field

[0001] This invention relates to the field of marine gravity measurement technology, and in particular to a method and system for comprehensive marine gravity correction based on multi-source auxiliary data. Background Technology

[0002] Ocean gravity fields have significant applications in geophysical exploration, resource surveys, and underwater navigation. Existing fusion methods lack collaborative constraints from multi-source auxiliary data and refined correction mechanisms for different frequency bands. Furthermore, the accuracy of motion correction is constrained by nonlinear interference from complex sea conditions, and anomaly detection and long-wave and short-wave noise reduction processes are relatively independent. Effectively fusing multi-source gravity data and achieving high-precision comprehensive correction has become a key technical problem that urgently needs to be solved in this field.

[0003] Current integrated ocean gravity correction techniques have significant limitations: First, motion correction relies on item-by-item compensation using physical formulas, making it difficult to adapt to the nonlinear error characteristics under complex sea conditions; second, grid calibration and drift correction are independent, lacking quantitative assessment of grid variation and spatial differences in drift, and anomaly identification is mostly based on a single threshold, failing to fully utilize long-wave and short-wave noise information and lacking spatial continuity verification; finally, long-wave and short-wave denoising lacks refined noise mode identification and classification strategies. Therefore, this invention proposes an integrated ocean gravity correction method and system based on multi-source auxiliary data. This method achieves high-precision reconstruction of the ocean gravity field driven by multi-source auxiliary data through motion correction of gravity observation sequences, calculation of grid variation coefficients and drift coefficients, calculation of noise contribution indices and construction of anomaly factors, and determination and correction of short-wave signal-to-noise ratio. This method combines multi-source auxiliary data fusion, anomaly grid identification, and physical constraint embedding, solving industry pain points such as low accuracy of ocean gravity data correction, error frequency band splitting, and multi-level error coupling, and is of great significance for improving the accuracy of ocean gravity field modeling and the reliability of underwater navigation. Summary of the Invention

[0004] The purpose of this invention is to provide a method and system for comprehensive marine gravity correction based on multi-source auxiliary data.

[0005] To achieve the above objectives, the present invention is implemented according to the following technical solution:

[0006] This invention includes the following steps:

[0007] Multi-source auxiliary data and gravity observation sequences were collected and spatiotemporally aligned. Motion corrections were applied to the gravity observation sequences to obtain gravity-corrected observation sequences.

[0008] The gravity correction observation sequence is divided into irregular grids based on seabed topography data. The gravity correction observation sequence is divided according to time and space. The grid value variation coefficient and drift coefficient are calculated and compared to obtain the gravity value sequence.

[0009] The tidal correction is calculated based on the third auxiliary data. The first standard gravity and the second standard gravity are obtained by subtracting them from the first and second auxiliary data, respectively. The standard gravity field is obtained by fusing them and long-wave noise and short-wave noise are extracted.

[0010] Calculate the long-wave noise contribution index and the short-wave noise contribution index and construct anomaly factors. Identify anomalous grid nodes in the gravity value sequence based on the anomaly factors.

[0011] The long-wave gravity component and the short-wave gravity component are extracted separately. The long-wave corrected gravity component is obtained by modal decomposition based on the reference gravity data. The short-wave signal-to-noise ratio is directly calculated and anomaly detection and short-wave correction are performed. The short-wave corrected gravity component is then output.

[0012] The gravity correction value is obtained by fusing the long-wavelength corrected gravity component and the short-wavelength corrected gravity component;

[0013] The multi-source auxiliary data includes first auxiliary data, second auxiliary data, third auxiliary data, navigation data, reference gravity data, and seabed topography data; the first auxiliary data is satellite altimetry gravity data; the second auxiliary data is dock baseline gravity data; and the third auxiliary data includes marine environmental observation data and tidal model data.

[0014] Furthermore, the method for performing spatiotemporal alignment includes:

[0015] Using the gravity observation sequence as a spatiotemporal reference, the multi-source auxiliary data is calibrated and preprocessed, and the spatial reference deviation and time delay of the multi-source auxiliary data relative to the gravity observation sequence are calculated.

[0016] Based on the spatial reference deviation, bilinear interpolation was used to spatially align the first auxiliary data, reference gravity data, seabed topography data, and marine environmental observation data; kriging interpolation was used to spatially align the second auxiliary data; and biharmonic spline interpolation was used to spatially align the tidal model data.

[0017] Based on the time delay, piecewise linear interpolation is used to time-align the second auxiliary data and the gravity observation sequence, and harmonic analysis matching is used to time-align the tidal model data and the gravity observation sequence.

[0018] Furthermore, the method for motion correction includes:

[0019] Perform six-degree-of-freedom motion compensation on navigation data, determine the six-degree-of-freedom motion parameter sequence and extract features, input the pre-trained support vector regression model, and output the motion disturbance correction amount;

[0020] The motion interference correction is removed from the gravity observation sequence to obtain the initial gravity observation sequence. The difference between the modified initial gravity observation sequence and the baseline gravity data is used to obtain the first residual.

[0021] A time-domain analysis is performed on the first residual to derive the autocorrelation function and the experimental semivariogram. The correlation length and the range of the variogram are calculated based on the autocorrelation function and the experimental semivariogram, respectively. The time correlation type of the first residual is determined based on the decay characteristics. The time correlation types include white noise-dominated, low-frequency trend-dominated, and periodic component-dominated types.

[0022] Based on the type of time correlation, select the Kriging interpolation method: if it is dominated by white noise, perform Kriging interpolation directly; if it is dominated by low-frequency trends, separate the trend term first and then perform Kriging interpolation; if it is dominated by periodic components, remove the periodic components first and then perform Kriging interpolation.

[0023] Kriging interpolation is used to spatially interpolate the first residual to obtain the second residual. The radius and lag distance parameters of the Kriging interpolation are adjusted according to the ratio of the correlation length to the range of the variogram. The second residual is added to the predicted gravity anomaly value for correction, resulting in the gravity-corrected observation sequence.

[0024] Furthermore, the method for obtaining the gravity value sequence includes:

[0025] Based on the topographic complexity of the seabed topography data, the gravity correction observation sequence is divided into an irregular grid.

[0026] The irregular grid is divided into early, middle and late segments according to the time dimension, and into flat terrain area and complex terrain area according to the spatial dimension, resulting in multi-temporal and spatial segmentation;

[0027] In a multi-temporal segment, grid value-drift coupling equations are established. The least squares adjustment method is used to calculate the grid value and drift rate of each temporal segment. The grid value variation coefficient and drift coefficient are calculated and compared to determine the error type. The gridded gravity correction observation sequence is then corrected according to the error type.

[0028] Furthermore, the method for calculating the tidal correction includes:

[0029] Instantaneous tide heights at different observation times are calculated based on tidal model data, and seawater density profile data are obtained based on marine observation environment data to construct a tide-gravity transfer function.

[0030] Substituting the instantaneous tide height into the tide level-gravity transfer function yields the tidal correction.

[0031] Furthermore, the method for identifying abnormal grid nodes includes:

[0032] The long-wave gravity component and short-wave gravity component in the standard gravity field are extracted. Allan variance decomposition is performed on the long-wave gravity component and short-wave gravity component respectively to obtain the variance estimates of long-wave noise and short-wave noise. The long-wave noise contribution index and short-wave noise contribution index of each grid node are calculated.

[0033] Based on the long-wave noise contribution index and the short-wave noise contribution index, anomaly factors for each grid node are constructed.

[0034] The abnormal factors are compared with the preset abnormal thresholds to determine candidate abnormal grid nodes;

[0035] Spatial continuity is assessed for candidate abnormal grid nodes. The absolute difference between the abnormality factor of a candidate abnormal grid node and the mean of the abnormality factor of its neighboring grid nodes is calculated. The quotient of the absolute difference and the standard deviation of the neighboring abnormal factors is used as the deviation degree, and a threshold is used to determine the abnormal grid nodes.

[0036] Furthermore, the method for obtaining the long-wavelength corrected gravitational component includes:

[0037] Long-wave gravity components are extracted by spectral decomposition, and two-dimensional empirical mode decomposition is performed on the long-wave gravity components to obtain several intrinsic mode function components and residual data, thereby determining the noise-dominant mode data.

[0038] Triangulation is performed on the noise-dominant mode data to extract local extrema points.

[0039] Local extreme points are triangulated and their geometric features are extracted. The area uniformity and side length variation coefficient are calculated, and the noise type is determined. The noise types include uniform noise, isolated noise, and strip noise.

[0040] For uniform noise, threshold filtering is used to correct the maximum and minimum values ​​of abnormal gravity values ​​of triangles; for isolated noise, extreme points corresponding to triangles with abnormal areas are removed; for strip noise, directional filtering is used to correct it according to the direction of triangle elongation.

[0041] The denoised modal data is obtained, and the denoised modal data, the remaining undenoised intrinsic mode function components, and the residual data are merged to obtain the long-wavelength corrected gravity component.

[0042] Furthermore, the method for anomaly detection and shortwave correction includes:

[0043] The shortwave signal-to-noise ratio (SNR) of each grid node is calculated based on the shortwave gravity component. If the SNR is greater than the SNR threshold, the shortwave gravity component is output and fused with the corrected longwave gravity component to obtain the final corrected gravity field; otherwise, the shortwave component is corrected. The step of correcting the shortwave component includes:

[0044] The gradient directions of the shortwave gravity component and the seabed topography data are aligned and compared. The gradient direction deviation is calculated, and the gradient direction of the shortwave gravity component is corrected according to the gradient direction deviation to obtain the gradient-corrected shortwave gravity component.

[0045] The second derivative of the seabed topography data is calculated. Based on the second derivative, the concave and convex regions of the seabed topography are determined. The local extreme value positions of the shortwave gravity component and the seabed topography data are extracted, aligned and compared, and the extreme value position deviation is calculated. Based on the convex and concave regions and combined with the extreme value position deviation, the extreme value anomaly correction is performed on the shortwave gravity component to obtain the extreme value corrected shortwave gravity component.

[0046] The gradient-corrected shortwave gravity component and the extremum-corrected shortwave gravity component are fused to obtain the shortwave corrected gravity component.

[0047] Secondly, a comprehensive ocean gravity correction system based on multi-source auxiliary data includes:

[0048] Data alignment module: used to acquire multi-source auxiliary data and gravity observation sequences and perform spatiotemporal alignment;

[0049] Data processing module: used to perform motion correction on gravity observation sequence to obtain gravity-corrected observation sequence, divide gravity-corrected observation sequence into irregular grids according to seabed topography data, divide gravity-corrected observation sequence according to time and space, calculate grid value variation coefficient and drift coefficient, compare and correct to obtain gravity value sequence;

[0050] Fusion module: used to calculate the tidal correction based on the third auxiliary data, subtract it from the first auxiliary data and the second auxiliary data respectively to obtain the first standard gravity and the second standard gravity, and fuse them to obtain the standard gravity field;

[0051] Anomaly identification module: used to calculate the long-wave noise contribution index and short-wave noise contribution index based on long-wave noise and short-wave noise and construct anomaly factors, and to identify abnormal grid nodes in the gravity value sequence based on the anomaly factors;

[0052] Correction module: It is used to extract the long-wave gravity component and the short-wave gravity component respectively, perform mode decomposition based on the reference gravity data to obtain the long-wave corrected gravity component, directly calculate the short-wave signal-to-noise ratio and perform anomaly detection and short-wave correction, output the short-wave corrected gravity component, and fuse the long-wave corrected gravity component and the short-wave gravity component to obtain the gravity correction value.

[0053] The beneficial effects of this invention are:

[0054] This invention is a method and system for comprehensive ocean gravity correction based on multi-source auxiliary data. Compared with the prior art, this invention has the following technical advantages:

[0055] This invention achieves motion correction that combines data-driven and spatial constraints by performing feature engineering on the six-degree-of-freedom motion parameters and inputting them into a support vector regression model, combined with Kriging residual interpolation correction.

[0056] This invention achieves accurate identification of anomalous grid nodes by calculating the noise contribution index of long-wave and short-wave gravity components and constructing anomaly factors, combined with spatial continuity judgment.

[0057] This invention achieves physical constraint correction of shortwave components by determining the shortwave signal-to-noise ratio and using seabed topographic data to correct gradient direction and extreme value anomalies, thereby improving the accuracy of the corrected shortwave components. Attached Figure Description

[0058] Figure 1 This is a flowchart illustrating the steps of a comprehensive marine gravity correction method based on multi-source auxiliary data according to the present invention. Detailed Implementation

[0059] The present invention will be further described below through specific embodiments. The illustrative embodiments and descriptions herein are used to explain the present invention, but are not intended to limit the present invention.

[0060] This invention relates to a method and system for comprehensive ocean gravity correction based on multi-source auxiliary data, comprising the following steps:

[0061] like Figure 1 As shown, this embodiment includes the following steps:

[0062] Multi-source auxiliary data and gravity observation sequences were collected and spatiotemporally aligned. Motion corrections were applied to the gravity observation sequences to obtain gravity-corrected observation sequences.

[0063] The gravity correction observation sequence is divided into irregular grids based on seabed topography data. The gravity correction observation sequence is divided according to time and space. The grid value variation coefficient and drift coefficient are calculated and compared to obtain the gravity value sequence.

[0064] The tidal correction is calculated based on the third auxiliary data. The first standard gravity and the second standard gravity are obtained by subtracting them from the first and second auxiliary data, respectively. The standard gravity field is obtained by fusing them and long-wave noise and short-wave noise are extracted.

[0065] Calculate the long-wave noise contribution index and the short-wave noise contribution index and construct anomaly factors. Identify anomalous grid nodes in the gravity value sequence based on the anomaly factors.

[0066] The long-wave gravity component and the short-wave gravity component are extracted separately. The long-wave corrected gravity component is obtained by modal decomposition based on the reference gravity data. The short-wave signal-to-noise ratio is directly calculated and anomaly detection and short-wave correction are performed. The short-wave corrected gravity component is then output.

[0067] The gravity correction value is obtained by fusing the long-wavelength corrected gravity component and the short-wavelength corrected gravity component;

[0068] The multi-source auxiliary data includes first auxiliary data, second auxiliary data, third auxiliary data, navigation data, reference gravity data, and seabed topography data; the first auxiliary data is satellite altimetry gravity data; the second auxiliary data is dock baseline gravity data; the third auxiliary data includes marine environmental observation data and tidal model data.

[0069] The gravity observation sequence consists of ship-based gravity data from the sea surface.

[0070] In this embodiment, the method for performing spatiotemporal alignment includes:

[0071] Using the gravity observation sequence as a spatiotemporal reference, the multi-source auxiliary data is calibrated and preprocessed, and the spatial reference deviation and time delay of the multi-source auxiliary data relative to the gravity observation sequence are calculated.

[0072] Based on the spatial reference deviation, bilinear interpolation was used to spatially align the first auxiliary data, reference gravity data, seabed topography data, and marine environmental observation data; kriging interpolation was used to spatially align the second auxiliary data; and biharmonic spline interpolation was used to spatially align the tidal model data.

[0073] Based on the time delay, piecewise linear interpolation is used to time-align the second auxiliary data and the gravity observation sequence, and harmonic analysis matching is used to time-align the tidal model data and the gravity observation sequence.

[0074] In this embodiment, the method for motion correction includes:

[0075] Perform six-degree-of-freedom motion compensation on navigation data, determine the six-degree-of-freedom motion parameter sequence and extract features, input the pre-trained support vector regression model, and output the motion disturbance correction amount;

[0076] The motion interference correction is removed from the gravity observation sequence to obtain the initial gravity observation sequence. The first residual is obtained by subtracting the initial gravity correction sequence from the reference gravity data.

[0077] A time-domain analysis is performed on the first residual to derive the autocorrelation function and the experimental semivariogram. The correlation length and the range of the variogram are calculated based on the autocorrelation function and the experimental semivariogram, respectively. The time correlation type of the first residual is determined based on the decay characteristics. The time correlation types include white noise-dominated, low-frequency trend-dominated, and periodic component-dominated types.

[0078] Based on the type of time correlation, select the Kriging interpolation method: if it is dominated by white noise, perform Kriging interpolation directly; if it is dominated by low-frequency trends, separate the trend term first and then perform Kriging interpolation; if it is dominated by periodic components, remove the periodic components first and then perform Kriging interpolation.

[0079] Kriging interpolation is used to spatially interpolate the first residual to obtain the second residual. The radius and lag distance parameters of the Kriging interpolation are adjusted according to the ratio of the correlation length to the range of the variogram. The second residual is added to the predicted gravity anomaly value for correction, resulting in the gravity-corrected observation sequence.

[0080] In this embodiment, the method for obtaining the gravity value sequence includes:

[0081] Based on the topographic complexity of the seabed topography data, the gravity correction observation sequence is divided into an irregular grid.

[0082] The irregular grid is divided into early, middle and late segments according to the time dimension, and into flat terrain area and complex terrain area according to the spatial dimension, resulting in multi-temporal and spatial segmentation;

[0083] In a multi-temporal segment, grid value-drift coupling equations are established. The least squares adjustment method is used to calculate the grid value and drift rate of each temporal segment. The grid value variation coefficient and drift coefficient are calculated and compared to determine the error type. The gridded gravity correction observation sequence is then corrected according to the error type.

[0084] In practical evaluation, the expression for the grid value-drift coupling equation is:

[0085]

[0086] in, for Gravity observations at time 10:00. For observation time, These are the gravimeter readings. These are the original readings from the gravimeter. This is the initial offset term for the gravimeter. The initial observation time, This represents the signed zero-point drift rate;

[0087] The expression for the drift coefficient is:

[0088]

[0089] in, To find the maximum value function, The drift coefficient, For the drift rate of complex terrain and , For flat terrain, the drift rate and ;

[0090] When the grid value of a specific spatiotemporal segment is greater than the preset grid value threshold, the current spatiotemporal segment grid value is removed, and the grid value and drift rate are calculated in the next spatiotemporal segment until the current grid value is less than or equal to the preset grid value threshold.

[0091] The error dominance type is determined by comparing the grid value variation coefficient and the drift coefficient. Correction is then made based on the error dominance type by comparing the grid value variation coefficient and the drift coefficient with the first preset threshold and the second preset threshold, respectively.

[0092] When the coefficient of variation of grid value First preset threshold, drift coefficient The second preset threshold is used to determine that grid value error is the dominant factor, and grid value iterative correction is performed.

[0093] When the coefficient of variation of grid value First preset threshold, drift coefficient The second preset threshold is used to determine that the drift error is dominant, and drift segment fitting correction is performed.

[0094] When the coefficient of variation of grid value First preset threshold, drift coefficient The second preset threshold is determined to be dominated by coupling error, and grid value-drift coupling iterative correction is performed.

[0095] When the coefficient of variation of grid value First preset threshold, drift coefficient If the second preset threshold is determined to be without significant error, the current grid gravity value will be directly output.

[0096] In this embodiment, the method for calculating the tidal correction includes:

[0097] Instantaneous tide heights at different observation times are calculated based on tidal model data, and seawater density profile data are obtained based on marine observation environment data to construct a tide-gravity transfer function.

[0098] Substituting the instantaneous tide height into the tide level-gravity transfer function, we obtain the tidal correction.

[0099] In actual assessments, the expression for the tide level-gravity transfer function is:

[0100]

[0101] in, For observation point exist The tide level at any given moment - gravity. Location of the observation point The total number of depth layers in the sea column. This is the depth layer index number. For the first The density of seawater at each depth layer, in units of , For the first The depth of each depth layer For the first The thickness of each depth layer For the first Seawater at the observation point Place The gravitational potential at a given moment. For the first The partial derivative of the change in seawater mass with respect to gravity, in units of , This refers to the instantaneous tide level.

[0102] The expression for the tidal correction is:

[0103]

[0104] in, For observation point exist Tidal correction at time, in units of , For observation point exist The instantaneous tide level at any given moment. For the first Seawater load at observation point Elastic deformation at the point of application due to gravity.

[0105] In this embodiment, the method for identifying abnormal grid nodes includes:

[0106] The long-wave gravity component and short-wave gravity component in the standard gravity field are extracted. Allan variance decomposition is performed on the long-wave gravity component and short-wave gravity component respectively to obtain the variance estimates of long-wave noise and short-wave noise. The long-wave noise contribution index and short-wave noise contribution index of each grid node are calculated.

[0107] Based on the long-wave noise contribution index and the short-wave noise contribution index, anomaly factors are constructed for each grid node.

[0108] The abnormal factors are compared with the preset abnormal thresholds to determine candidate abnormal grid nodes;

[0109] Spatial continuity is judged for candidate abnormal grid nodes. The absolute difference between the abnormality factor of the candidate abnormal grid node and the mean of the abnormality factor of the neighboring grid nodes is calculated. The quotient of the absolute difference and the standard deviation of the neighboring abnormal factors is used as the deviation degree. Threshold judgment is performed to determine the abnormal grid nodes.

[0110] In practical evaluation, the expression for the long-wave noise contribution index is:

[0111]

[0112] in, For grid nodes Longwave noise contribution index This is the latitudinal index of the grid node. For the meridional index of the grid node, This represents the latitudinal lag step size for the long-wavelength component. This represents the meridional lag step size of the long-wavelength component. For long-wavelength gravity components at spatial lag scales Allan variance contribution in two-dimensional space The contribution of long-wavelength signal components to the variance at spatial scale;

[0113] The expression for the shortwave noise contribution index is:

[0114]

[0115] in, For grid nodes Shortwave noise contribution index The latitudinal lag step size for the shortwave component is... This represents the meridional lag step size of the shortwave component. For shortwave gravity components at spatial lag scales Allan variance contribution in two-dimensional space The contribution of shortwave signal components to the variance at spatial scales;

[0116] The expression for the abnormal factor is:

[0117]

[0118] in, For grid nodes Abnormal factors, Contribute index weighting factors to long waves. Contribute index weighting factor to shortwave;

[0119] If abnormal factors If a preset anomaly threshold is set, the grid node will be marked as a candidate anomaly grid node.

[0120] If the deviation A preset deviation threshold is used to determine whether a candidate abnormal grid node is a normal grid node;

[0121] If the deviation A preset deviation threshold is set to determine the candidate abnormal grid node as an abnormal grid node, and the long and short wave correction process is initiated.

[0122] The calculation method of the two-dimensional spatial Allan variance is to select different spatial lag step lengths, perform difference calculation on a pair of grid points separated by a step length in the two-dimensional grid, sum the squares of all the differences that meet the conditions in the whole figure, and then take half of the average. The difference of grid values ​​separated by a given step length is compared point by point in the two-dimensional spatial domain.

[0123] In this embodiment, the method for obtaining the long-wavelength corrected gravity component includes:

[0124] Long-wave gravity components are extracted by spectral decomposition, and two-dimensional empirical mode decomposition is performed on the long-wave gravity components to obtain several intrinsic mode function components and residual data, thereby determining the noise-dominant mode data.

[0125] Triangulation is performed on the noise-dominant mode data to extract local extrema points.

[0126] Local extreme points are triangulated and their geometric features are extracted. The area uniformity and side length variation coefficient are calculated, and the noise type is determined. The noise types include uniform noise, isolated noise, and strip noise.

[0127] For uniform noise, threshold filtering is used to correct the maximum and minimum values ​​of abnormal gravity values ​​of triangles; for isolated noise, extreme points corresponding to triangles with abnormal areas are removed; for strip noise, directional filtering is used to correct it according to the direction of triangle elongation.

[0128] The denoised modal data is obtained, and the denoised modal data, the remaining undenoised intrinsic modal function components, and the residual data are merged to obtain the long-wavelength corrected gravity components.

[0129] In actual evaluation, the cumulative correlation coefficient between each intrinsic mode function component and the reference gravity data is calculated separately. Intrinsic mode function components whose cumulative correlation coefficient change rate is less than the preset change rate threshold are judged as low correlation mode data.

[0130] For low-correlation modal data, the intrinsic mode function components are accumulated sequentially in descending order of frequency, and the cumulative variance contribution rate after each accumulation is calculated.

[0131] When the cumulative variance contribution rate reaches the preset contribution rate threshold, record the current accumulation count K, and determine the first K intrinsic mode function components as noise-dominant mode data;

[0132] When the area distribution uniformity First threshold, coefficient of variation of side length The second threshold is used to determine uniform noise, and threshold filtering is used to correct the maximum and minimum values ​​of abnormal gravity values ​​of the triangle.

[0133] When the area distribution uniformity If an abnormal triangle with an area greater than the third threshold is found to be isolated noise, the extreme points corresponding to the abnormal triangle with an area greater than the third threshold are removed, the triangle mesh is reconstructed, and spatial interpolation is performed to fill it.

[0134] When the area distribution uniformity The first threshold and the triangles are distributed in a strip shape along the preset direction. The area uniformity does not meet the judgment conditions for uniform noise and isolated noise, so it is judged as strip noise. Directional filtering is performed according to the elongation direction of the triangle to suppress high-frequency changes perpendicular to the strip direction.

[0135] When the area uniformity does not meet the requirements for uniform noise, isolated noise, and strip noise, median noise reduction is used for processing.

[0136] In this embodiment, the method for anomaly detection and shortwave correction includes:

[0137] The shortwave signal-to-noise ratio (SNR) of each grid node is calculated based on the shortwave gravity component. If the SNR is greater than the SNR threshold, the shortwave gravity component is output and fused with the corrected longwave gravity component to obtain the final corrected gravity field; otherwise, the shortwave component is corrected. The step of correcting the shortwave component includes:

[0138] The gradient directions of the shortwave gravity component and the seabed topography data are aligned and compared. The gradient direction deviation is calculated, and the gradient direction of the shortwave gravity component is corrected according to the gradient direction deviation to obtain the gradient-corrected shortwave gravity component.

[0139] The second derivative of the seabed topography data is calculated. Based on the second derivative, the concave and convex regions of the seabed topography are determined. The local extreme value positions of the shortwave gravity component and the seabed topography data are extracted, aligned and compared, and the extreme value position deviation is calculated. Based on the convex and concave regions and combined with the extreme value position deviation, the extreme value anomaly correction is performed on the shortwave gravity component to obtain the extreme value corrected shortwave gravity component.

[0140] The gradient-corrected shortwave gravity component and the extremum-corrected shortwave gravity component are fused to obtain the shortwave corrected gravity component.

[0141] In practical evaluation, the expression for the shortwave gravity component after gradient correction is as follows:

[0142] ;

[0143] in, The shortwave gravity component after gradient correction. The original shortwave gravity component, To preset the correction adjustment coefficient, The gradient direction deviation angle. For grid nodes The gradient vector of the shortwave gravity component at that location. The mean value of the gradient vector of the short-wave gravity component of the full-grid network;

[0144] The expression for the shortwave gravity component after extreme value correction is:

[0145] ;

[0146] in, The shortwave gravity component after extreme value correction. To preset the correction strength coefficient, and These are the latitudinal and longitudinal deviations of the extreme value locations, respectively. The preset spatial attenuation coefficient, For seabed topographic data, For symbolic functions, This is the second derivative of the seabed topography data.

[0147] Furthermore, taking a certain sea area in the northern part of the South China Sea as an example, the longitude range is 112 degrees east. ~113 Latitude range 15°N ~16 The target area consists of 3600 grid points, with a uniform grid resolution set to 1'. 1´, approximately 1.85km;

[0148] Using the upper left corner of the target area as the origin, bilinear interpolation was used to spatially align the first auxiliary data, the baseline gravity data, and the seabed topography data, with the target grid point at longitude 112°E. 0'30", 15" North Latitude Taking 59'30" as an example, the satellite altimeter gravity values ​​at the four corner points are 37.31 mGal, 39.17 mGal, 35.83 mGal, and 37.51 mGal, respectively. After bilinear interpolation, the gravity value at that point is calculated to be 37.455 mGal. The seabed topography data within the same grid window ranges in depth from -1873 to -2633 m. The interpolation result smoothly connects with the original grid data. Seawater density profile data is provided in depth layers, with a surface density of 1023.7. As the depth increases to the bottom layer 1035.1 Bilinear interpolation is used to match each grid point;

[0149] The second auxiliary data source is the gravity data of the dock's base points. Ordinary kriging interpolation was used to normalize the data to a unified grid. Five discrete points on the ship were selected, with gravity values ​​ranging from 38.47 to 43.51 mGal. The spherical model parameters are: nugget value 0.73 mGal squared, abutment value 7.91 mGal. 0.17 Strokes Located at 112 degrees east longitude 15', 15' North Latitude Taking a 45' grid point as an example, the weights of the four known points involved in the interpolation are 0.273, 0.271, 0.229 and 0.227, respectively. The Kriging estimate is 40.69 mGal, the actual gravity measurement of the base point is 40.61 mGal, and the interpolation error is 0.08 mGal.

[0150] The tidal model data in the third auxiliary data uses biharmonic spline interpolation, with the target grid point at longitude 112°E. 30', 15' North Latitude Taking 30' as an example, the interpolated amplitude of the M2 tidal bore is 44.3 cm, and the lag angle is 226°. ;

[0151] Multi-source auxiliary data were unified to the WGS84 coordinate system and the EGM2008 geoid datum, and stored in a 1´ resolution grid. The range of satellite altimetry gravity anomaly was 34~41mGal, the range of wharf benchmark gravity anomaly was 38~44mGal, the range of water depth was -1873~-2633m, and the M2 tidal amplitude was 37 to 54cm.

[0152] Motion corrections are applied to the gravity observation sequence. Taking the observation time at 1800s as an example, feature engineering is performed on the six-degree-of-freedom motion parameter sequence at the observation time at 1800s to extract statistical features of each motion parameter with a sliding window length of 60s.

[0153] The motion state feature vector is input into a pre-trained support vector regression model. The training samples are derived from shipborne gravity measurement data from no less than five historical voyages in the target sea area and adjacent sea areas. Shipborne gravity outliers are removed by gross error elimination in overlapping areas with satellite altimetry gravity data (gross error elimination is based solely on the consistency of repeated survey lines within the shipborne survey data and the absolute accuracy verification of land reference points). Voyage data with a root mean square error of repeated survey lines better than 0.5 mGal are retained as training samples, and shipborne gravity outliers are used as label ground truth values. Accuracy is transferred through shore-based gravity benchmark connection measurements to ensure absolute accuracy better than 0.3 mGal.

[0154] Each voyage includes complete six-DOF motion parameters and corresponding ship-measured gravity anomalies. The sampling frequency is 1Hz. For each motion parameter, the mean, standard deviation, and root mean square value within a sliding window of 60s are extracted to form an 18-dimensional feature vector. Principal component analysis is used to reduce the dimensionality to 12 dimensions. The support vector regression model uses a radial basis function kernel with a penalty function parameter of 10 and a kernel function parameter of 0.05. Training uses 5-fold cross-validation and quality control. The optimal hyperparameters are selected with the goal of minimizing the root mean square error. At time 1800s, the motion disturbance correction output by the support vector regression model is 2.69mGal, the initial corrected gravity value is 37.51mGal, and the extracted reference gravity anomaly value at this time is 39.47mGal. The difference between the gravity correction observation sequence and the reference gravity anomaly value is used to obtain the first residual of 1.96mGal.

[0155] For the experimental semivariogram, all grid point pairs are divided into 20 lag intervals based on distance (starting from 0.01 degrees and increasing in increments of 0.01 degrees to 0.20 degrees). The experimental semivariogram value is calculated within each lag interval using the following formula:

[0156]

[0157] in, The lag distance is The experimental semivariogram at time , The lag distance, The total number of sample points. , Let each be any two points on the grid. For index number, , They are respectively in , The first residual value at the two points;

[0158] The scatter sequence of experimental semivariogram values ​​as a function of lag distance was obtained by calculating interval by interval. All 3600 grid points were paired up to obtain the experimental semivariogram scatter points. The theoretical model of the experimental semivariogram was fitted using an exponential model and the fitting method was weighted least squares (the weight of each lag distance interval was the reciprocal of the logarithm of the grid points in that interval, so that the interval with more pairs had a higher weight in the fitting).

[0159] The autocorrelation function and experimental semivariogram were derived. Autocorrelation function analysis showed that the autocorrelation coefficient decayed slowly and remained positive over a long lag time, indicating a significant low-frequency trend component in the first residual. Experimental semivariogram analysis showed that the nugget value was 0.73. The base value is 7.91. The range is 0.17. (Approximately 18.9 km) Based on the slowly decaying characteristic of the autocorrelation function, the first residual is determined to be dominated by a low-frequency trend, with the time-domain correlation length set at 27 s and the spatial-domain variogram range at 0.17. ;

[0160] Determined to be a low-frequency trend-dominant type, the trend term is first separated from the first residual, and then fitted using the least squares method to obtain the trend term expression: The trend term has a value of 2.15 mGal at 1800s, and the mean of the detrended residuals is -0.19 mGal. Based on a ship speed of 3.9... Convert the time-domain correlation length into spatial distance ,in For spatial distance, The length of the correlation in the time domain. Given the ship's speed, the spatial distance is approximately 0.00095. The ratio of the transformed spatial distance to the range of the spatial domain variogram was calculated and found to be 0.006, which is less than the preset threshold of 1.0. This indicates that the spatial correlation of the residual is significantly stronger than the temporal correlation. When performing Kriging interpolation, the range of the spatial domain variogram of 0.17 degrees was used as the main search radius reference, and the lag distance parameter was appropriately reduced to highlight the local spatial structure characteristics.

[0161] Dynamically adjust the search radius of the Kriging interpolation (from 0.30) Adjusted to 0.23 ) and hysteresis parameter (from 0.017 Adjusted to 0.013 The second residual is obtained by spatial interpolating the detrended residual using the adjusted Kriging interpolation parameters.

[0162] Located at 112 degrees east longitude 19', 15' North Latitude Taking a 47' grid point as an example, the distances between the four known points involved in the Kriging interpolation are 0.07. 0.09 0.11 and 0.13 The corresponding weights are 0.37, 0.31, 0.19 and 0.13, respectively. The second residual obtained by interpolation is -0.23mGal. The second residual is added to the predicted gravity anomaly value for correction, and the corrected gravity anomaly value for the grid point is 37.28mGal. After motion correction, the root mean square error of the gravity correction observation sequence is reduced from 2.37mGal before correction to 1.03mGal.

[0163] Based on the topographic complexity of the seabed topographic data, the gravity correction observation sequence was divided into an irregular grid. The sea area with a topographic standard deviation greater than 220m was identified as a complex topographic area, and the sea area with a standard deviation less than 150m was identified as a flat topographic area. The complex topographic area accounts for about 37% of the total area in this target area, mainly concentrated in the central-eastern part of the study area. The flat topographic area accounts for about 51% of the total area, distributed in the western and southern parts of the study area. The remaining 12% is a transitional area.

[0164] The irregular grid is divided into an early segment (0-2400s), a middle segment (2400-4800s), and a late segment (4800-7200s) according to the time dimension, and into a flat terrain region and a complex terrain region according to the spatial dimension, resulting in 6 spatiotemporal segments: flat terrain - early segment, flat terrain - middle segment, flat terrain - late segment, complex terrain - early segment, complex terrain - middle segment, and complex terrain - late segment;

[0165] In each spatiotemporal segment, a grid value-drift coupling equation is established. If the normal equation becomes ill-conditioned during the adjustment process, conventional methods such as ridge estimation are used to handle it. The least squares adjustment method is used to perform joint inversion of grid values ​​and drift rates for each spatiotemporal segment. The Helmert variance component estimation method is used to determine the weights of each segment. The inversion is performed on the six spatiotemporal segments in sequence to obtain the grid values ​​and drift rates of each segment.

[0166] Calculate the coefficient of variation for each grid value; the arithmetic mean of the six segmented grid values ​​is: The standard deviation of the grid values ​​is 0.0113. The calculated value is 0.00353;

[0167] The drift coefficient was calculated, and the average drift rate of the three segments in the complex terrain area was 0.093. The average drift rate of the three segments in the flat terrain area was 0.087. The calculated value is 0.0645;

[0168] The drift coefficient calculation results are further used to determine the drift direction: when and When the signs are the same, it is determined that they are drifting in the same direction, indicating that the drift trends of the two regions are consistent; when and When the signs are opposite, it is determined to be a reverse drift, which may indicate systematic bias or complex sea state interference. and The ratio of absolute values ​​(the one greater than 1) is used as the corrected drift coefficient. In subsequent judgments, if the corrected drift coefficient exceeds the preset threshold, the drift segment fitting correction is triggered first.

[0169] The first preset threshold is the upper limit of relative change of grid value determined based on the statistical results of grid value calibration of historical gravimeters of the same type; the second preset threshold is the upper limit of drift difference determined based on the statistical results of drift rate difference between flat terrain areas and complex terrain areas.

[0170] The grid value variation coefficient is compared with the first preset threshold (0.5%) and the drift coefficient is compared with the second preset threshold (10%). The grid value variation coefficient of 0.353% is less than 0.5% and the drift coefficient of 6.45% is less than 10%. Since neither of the two indicators exceeds the threshold, it is determined that there is no significant error and no correction needs to be triggered.

[0171] The tidal correction was calculated using the load Green's function to account for the elastic deformation and gravity effects of seawater loading. A one-dimensional Earth model was used, with seawater densities ranging from 1023.7 to 1035.1 at various depths. Target grid point 112 East longitude 30', 15' North Latitude The transfer function at position 30' is 37.7305. At 1800s, the instantaneous tidal level was -0.13m, the direct gravitational effect was -4.9050mGal, the elastic deformation effect was -0.0063mGal, the total tidal correction was -4.9113mGal, the tidal correction range for the entire region was -6.13~5.47mGal, and the average was -0.73mGal.

[0172] The long-wave noise contribution index uses the step size corresponding to the stability of the variance value in the Allan variance curve as the lag step size, while the short-wave noise contribution index uses the minimum step size corresponding to the significant increase of the variance value in the Allan variance curve as the lag step size.

[0173] Allan variance boundary handling employs an effective region constraint method: for a given latitudinal and meridional lag step, only if the sum of the latitudinal index and the latitudinal lag step of a grid node does not exceed the total number of rows in the grid, and the sum of the meridional index and the meridional lag step does not exceed the total number of columns in the grid, does the node participate in the difference calculation at that lag step; grid point pairs that do not meet the boundary conditions are directly discarded and do not participate in the variance calculation at that lag step.

[0174] A spatial Allan variance curve is plotted with the lag step size as the horizontal axis and the Allan variance value as the vertical axis. The curve shows different trends in the short-scale and long-scale intervals. The curve is fitted to the short-scale and long-scale intervals respectively. The fitted value in the short-scale interval represents the noise variance contribution, and the fitted value in the long-scale interval represents the signal variance contribution. The ratio of the two represents the proportion of frequency band noise.

[0175] Abnormal grid nodes were identified. Within a 5×5 window, the long-wave component ranged from 36.73 to 38.03 mGal, and the short-wave component ranged from -1.07 to 2.87 mGal. The long-wave noise contribution index of the center node in the window was 0.187, and the short-wave noise contribution index was 0.241. Based on experimental calibration, the weighting factor for the long-wave contribution index was determined to be 0.4, and the weighting factor for the short-wave contribution index was determined to be 0.6. The distribution range of abnormal factors for all grid nodes was 0.143 to 0.230. Nodes with an abnormal factor less than the preset abnormal threshold of 0.17 were identified as normal nodes. Four candidate anomalous nodes were marked with a threshold of 0.17, with anomalous factors of 0.176, 0.181, 0.187, and 0.230, respectively. The spatial continuity judgment deviations were 4.69, 1.14, 0.71, and 0.43, respectively. The preset deviation threshold was the upper limit of the deviation based on the assumption of normal distribution, which was set to 2.0. Only the node with an anomalous factor of 0.230 was determined to be an anomalous grid node. A total of 17 anomalous nodes were detected in the whole area, accounting for 0.47%. The average shortwave noise contribution index of 0.473 was 2.11 times that of the average of normal nodes of 0.224.

[0176] As a preferred approach, after calculating the anomaly factor, the grid nodes identified as candidate anomalies are subjected to a robust secondary verification: the anomaly factor of the node is calculated for three window sizes: 3×3, 5×5, and 7×7. If the anomaly factor of the node is greater than the preset anomaly threshold for all three window sizes, it is finally confirmed as an anomalous grid node. If it is greater than the threshold only for some window sizes, it is marked as a suspected anomaly node and enters the spatial continuity judgment step for further confirmation, thereby further enhancing the robustness of anomaly identification.

[0177] The spatial frequency division between long-wavelength and short-wavelength gravity components is based on a spatial wavelength of 50 km. Wavelengths greater than 50 km (spatial frequencies less than...) The component of gravity is defined as the long-wavelength gravity component, with a wavelength less than or equal to 50 km (spatial frequency greater than or equal to 50 km). The components of ) are defined as short-wave gravity components;

[0178] In long-wave correction, the screening termination condition for two-dimensional empirical mode decomposition is that the number of extreme points in the remaining residual data is less than the preset extreme point number threshold (it is impossible to continue extracting effective intrinsic mode function components) or the number of intrinsic mode function components obtained by decomposition reaches the preset upper limit of the decomposition layer. The decomposition of 7 intrinsic mode function components is terminated. The correlation coefficient ranges from 0.197 to 0.897. The cumulative variance contribution rate is calculated with the total variance of the original long-wave gravity components as the denominator. When the cumulative variance contribution rate threshold is 60%, the first 4 components are determined to be noise-dominant modes.

[0179] Statistical analysis was performed on historical gravity data of the target area to calculate the mean and standard deviation of the uniformity of triangle area distribution and the coefficient of variation of side length in the normal area. The mean minus 2 times the standard deviation was used as the lower limit of the first threshold, and the mean plus 2 times the standard deviation was used as the upper limit of the second threshold.

[0180] The first threshold is set to 0.60, the second threshold to 0.50, and the third threshold to twice the average triangle area. After triangulation, the triangle area uniformity is 0.631 and the side length variation coefficient is 0.412, which is determined to be uniform noise. The maximum value of the typical triangle (0.53 mGal) is corrected to 0.29 mGal and the minimum value (-0.37 mGal) is corrected to -0.23 mGal. The long-wave signal-to-noise ratio is improved from 28.3 dB to 34.7 dB, and the root mean square error is reduced from 0.37 mGal to 0.19 mGal.

[0181] In shortwave correction, the shortwave signal-to-noise ratio is calculated using the conventional signal-to-noise ratio calculation method. The ratio of the signal power to the noise power of the shortwave gravity component within the spatial neighborhood window is calculated, and the signal quality of the shortwave component is characterized in decibels.

[0182] Based on experimental calibration, the preset correction intensity coefficient was set to 0.2, the preset spatial attenuation coefficient was set to 1.5 grid spacings, the shortwave component range was -2.37 to 3.53 mGal, the average shortwave signal-to-noise ratio at anomalous nodes was 14.7 dB, the preset signal-to-noise ratio threshold was 18.0 dB, and the gradient direction deviation angle at typical anomalous nodes was 46°. Greater than the threshold 25 After correction, the shortwave component changed from -1.73mGal to -1.49mGal. The extreme value position deviation modulus of 1.12´ was greater than the threshold of 1.0´, triggering correction. The corrected value is -1.94mGal.

[0183] This embodiment uses a dynamic weighted average driven by the signal-to-noise ratio (SNR) to fuse the shortwave gravity component after gradient correction and the shortwave gravity component after extreme value correction: for grid nodes with a shortwave SNR higher than a preset SNR threshold, the gradient correction result accounts for 0.7 and the extreme value correction result accounts for 0.3; for grid nodes with a shortwave SNR lower than the preset SNR threshold, the gradient correction result accounts for 0.4 and the extreme value correction result accounts for 0.6.

[0184] The weighted fusion of the shortwave correction component of this node is -1.80mGal, the root mean square error of the outlier node is reduced from 1.07mGal to 0.53mGal, and the correlation between the shortwave component and the terrain is improved from 0.31 to 0.87.

[0185] After superposition and fusion, the final gravity correction value ranges from 33.97 to 45.23 mGal, with a root mean square error of 1.37 mGal. The root mean square error of 1.37 mGal was verified by 12 ship-based survey stations in the target area, which is about 42.2% higher than the initial value of 3.31 mGal.

[0186] Secondly, a comprehensive ocean gravity correction system based on multi-source auxiliary data includes:

[0187] Data alignment module: used to acquire multi-source auxiliary data and gravity observation sequences and perform spatiotemporal alignment;

[0188] Data processing module: used to perform motion correction on gravity observation sequence to obtain gravity-corrected observation sequence, divide gravity-corrected observation sequence into irregular grids according to seabed topography data, divide gravity-corrected observation sequence according to time and space, calculate grid value variation coefficient and drift coefficient, compare and correct to obtain gravity value sequence;

[0189] Fusion module: used to calculate the tidal correction based on the third auxiliary data, subtract it from the first auxiliary data and the second auxiliary data respectively to obtain the first standard gravity and the second standard gravity, and fuse them to obtain the standard gravity field;

[0190] Anomaly identification module: used to calculate the long-wave noise contribution index and short-wave noise contribution index based on long-wave noise and short-wave noise and construct anomaly factors, and to identify abnormal grid nodes in the gravity value sequence based on the anomaly factors;

[0191] Correction module: It is used to extract the long-wave gravity component and the short-wave gravity component respectively, perform mode decomposition based on the reference gravity data to obtain the long-wave corrected gravity component, directly calculate the short-wave signal-to-noise ratio and perform anomaly detection and short-wave correction, output the short-wave corrected gravity component, and fuse the long-wave corrected gravity component and the short-wave gravity component to obtain the gravity correction value.

[0192] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for comprehensive ocean gravity correction based on multi-source auxiliary data, characterized in that, Includes the following steps: S1. Collect multi-source auxiliary data and gravity observation sequences and perform spatiotemporal alignment. Perform motion correction on the gravity observation sequences to obtain gravity-corrected observation sequences. S2. Based on the seabed topography data, the gravity correction observation sequence is divided into irregular grids. The gravity correction observation sequence is divided according to time and space. The grid value variation coefficient and drift coefficient are calculated and compared to obtain the gravity value sequence. S3. Calculate the tidal correction based on the third auxiliary data, and subtract it from the first and second auxiliary data to obtain the first standard gravity and the second standard gravity, and then merge them to obtain the standard gravity field. S4. Calculate the long-wave noise contribution index and the short-wave noise contribution index and construct the anomaly factor. Identify the abnormal grid nodes in the gravity value sequence based on the anomaly factor. S5. Extract the long-wave gravity component and the short-wave gravity component respectively. Perform mode decomposition based on the reference gravity data to obtain the long-wave corrected gravity component. Directly calculate the short-wave signal-to-noise ratio and perform anomaly detection and short-wave correction. Output the short-wave corrected gravity component. S6. The gravity correction value is obtained by fusing the long-wavelength corrected gravity component and the short-wavelength corrected gravity component; The multi-source auxiliary data includes first auxiliary data, second auxiliary data, third auxiliary data, reference gravity data, navigation data, and seabed topography data; the first auxiliary data is satellite altimetry gravity data; the second auxiliary data is dock baseline gravity data; and the third auxiliary data includes marine environmental observation data and tidal model data.

2. The ocean gravity comprehensive correction method based on multi-source auxiliary data according to claim 1, characterized in that, The method for spatiotemporal alignment includes: Using the gravity observation sequence as a spatiotemporal reference, the multi-source auxiliary data is calibrated and preprocessed, and the spatial reference deviation and time delay of the multi-source auxiliary data relative to the gravity observation sequence are calculated. Based on the spatial reference deviation, bilinear interpolation was used to spatially align the first auxiliary data, reference gravity data, seabed topography data, marine environmental observation data, and gravity observation sequence; kriging interpolation was used to spatially align the second auxiliary data and gravity observation sequence; and biharmonic spline interpolation was used to spatially align the tidal model data and gravity observation sequence. Based on the time delay, the gravity observation sequence is time-aligned using piecewise interpolation, and the tidal model data and gravity observation sequence are time-aligned using harmonic analysis matching.

3. The ocean gravity comprehensive correction method based on multi-source auxiliary data according to claim 1, characterized in that, The method for motion correction includes: Perform six-degree-of-freedom motion compensation on navigation data, determine the six-degree-of-freedom motion parameter sequence and extract features, input the pre-trained support vector regression model, and output the motion disturbance correction amount; The motion interference correction is removed from the gravity observation sequence to obtain the initial gravity observation sequence. The first residual is obtained by subtracting the initial gravity observation sequence from the reference gravity data. A time-domain analysis is performed on the first residual to derive the autocorrelation function and the experimental semivariogram. The correlation length and the range of the variogram are calculated based on the autocorrelation function and the experimental semivariogram, respectively. The time correlation type of the first residual is determined based on the decay characteristics. The time correlation types include white noise-dominated, low-frequency trend-dominated, and periodic component-dominated types. Based on the type of time correlation, select the Kriging interpolation method: if it is dominated by white noise, perform Kriging interpolation directly; if it is dominated by low-frequency trends, separate the trend term first and then perform Kriging interpolation; if it is dominated by periodic components, remove the periodic components first and then perform Kriging interpolation. Kriging interpolation is used to spatially interpolate the first residual to obtain the second residual. The radius and lag distance parameters of the Kriging interpolation are adjusted according to the ratio of the correlation length to the range of the variogram. The second residual is added to the predicted gravity anomaly value for correction, resulting in the gravity-corrected observation sequence.

4. The ocean gravity comprehensive correction method based on multi-source auxiliary data according to claim 1, characterized in that, The method for obtaining the gravity value sequence includes: Based on the topographic complexity of the seabed topography data, the gravity correction observation sequence is divided into an irregular grid. The irregular grid is divided into early, middle and late segments according to the time dimension, and into flat terrain area and complex terrain area according to the spatial dimension, resulting in multi-temporal and spatial segmentation; In a multi-temporal segment, grid value-drift coupling equations are established. The least squares adjustment method is used to calculate the grid value and drift rate of each temporal segment. The grid value variation coefficient and drift coefficient are calculated and compared to determine the error type. The gridded gravity correction observation sequence is then corrected according to the error type.

5. The ocean gravity comprehensive correction method based on multi-source auxiliary data according to claim 1, characterized in that, The method for calculating the tidal correction includes: Instantaneous tide heights at different observation times are calculated based on tidal model data, and seawater density profile data are obtained based on marine observation environment data to construct a tide-gravity transfer function. Substituting the instantaneous tide height into the tide level-gravity transfer function yields the tidal correction.

6. The ocean gravity comprehensive correction method based on multi-source auxiliary data according to claim 1, characterized in that, The method for identifying abnormal grid nodes includes: The long-wave gravity component and short-wave gravity component in the standard gravity field are extracted. Allan variance decomposition is performed on the long-wave gravity component and short-wave gravity component respectively to obtain the variance estimates of long-wave noise and short-wave noise. The long-wave noise contribution index and short-wave noise contribution index of each grid node are calculated. Based on the long-wave noise contribution index and the short-wave noise contribution index, anomaly factors are constructed for each grid node. The abnormal factors are compared with the preset abnormal thresholds to determine candidate abnormal grid nodes; Spatial continuity is assessed for candidate abnormal grid nodes. The absolute difference between the abnormality factor of a candidate abnormal grid node and the mean of the abnormality factor of its neighboring grid nodes is calculated. The quotient of the absolute difference and the standard deviation of the abnormality factor of the neighboring grid nodes is used as the deviation degree, and a threshold is used to determine the abnormal grid nodes.

7. The ocean gravity comprehensive correction method based on multi-source auxiliary data according to claim 1, characterized in that, The method for obtaining the long-wavelength corrected gravity component includes: Long-wave gravity components are extracted by spectral decomposition, and two-dimensional empirical mode decomposition is performed on the long-wave gravity components to obtain several intrinsic mode function components and residual data, thereby determining the noise-dominant mode data. Triangulation is performed on the noise-dominant mode data to extract local extrema points. Local extreme points are triangulated and their geometric features are extracted. The area uniformity and side length variation coefficient are calculated, and the noise type is determined. The noise types include uniform noise, isolated noise, and strip noise. For uniform noise, threshold filtering is used to correct the maximum and minimum values ​​of abnormal gravity values ​​of triangles; for isolated noise, extreme points corresponding to triangles with abnormal areas are removed; for strip noise, directional filtering is used to correct it according to the direction of triangle elongation. The denoised modal data is obtained, and the denoised modal data, the remaining undenoised intrinsic mode function components, and the residual data are merged to obtain the long-wavelength corrected gravity component.

8. The ocean gravity comprehensive correction method based on multi-source auxiliary data according to claim 1, characterized in that, The method for anomaly detection and shortwave correction includes: The shortwave signal-to-noise ratio (SNR) of each grid node is calculated based on the shortwave gravity component. If the shortwave SNR is greater than the SNR threshold, the shortwave gravity component is output and fused with the corrected longwave gravity component to obtain the final corrected gravity field. Otherwise, the shortwave component is corrected. The step of correcting the shortwave component includes: The gradient directions of the shortwave gravity component and the seabed topography data are calculated and aligned. The gradient direction deviation is calculated, and the gradient direction of the shortwave gravity component is corrected according to the gradient direction deviation to obtain the gradient-corrected shortwave gravity component. The second derivative of the seabed topography data is calculated. Based on the second derivative, the concave and convex regions of the seabed topography are determined. The local extreme value positions of the shortwave gravity component and the seabed topography data are extracted, aligned and compared, and the extreme value position deviation is calculated. Based on the convex and concave regions and combined with the extreme value position deviation, the extreme value anomaly correction is performed on the shortwave gravity component to obtain the extreme value corrected shortwave gravity component. The gradient-corrected shortwave gravity component and the extremum-corrected shortwave gravity component are fused to obtain the shortwave corrected gravity component.

9. A marine gravity integrated correction system based on multi-source auxiliary data, used to perform the method according to any one of claims 1-8, characterized in that, include: Data alignment module: used to acquire multi-source auxiliary data and gravity observation sequences and perform spatiotemporal alignment; Data processing module: used to perform motion correction on gravity observation sequence to obtain gravity-corrected observation sequence, divide gravity-corrected observation sequence into irregular grids according to seabed topography data, divide gravity-corrected observation sequence according to time and space, calculate grid value variation coefficient and drift coefficient, compare and correct to obtain gravity value sequence; Fusion module: used to calculate the tidal correction based on the third auxiliary data, subtract it from the first auxiliary data and the second auxiliary data respectively to obtain the first standard gravity and the second standard gravity, and fuse them to obtain the standard gravity field; Anomaly identification module: used to calculate the long-wave noise contribution index and short-wave noise contribution index based on long-wave noise and short-wave noise and construct anomaly factors, and to identify abnormal grid nodes in the gravity value sequence based on the anomaly factors; Correction module: It is used to extract the long-wave gravity component and the short-wave gravity component respectively, perform mode decomposition based on the reference gravity data to obtain the long-wave corrected gravity component, directly calculate the short-wave signal-to-noise ratio and perform anomaly detection and short-wave correction, output the short-wave corrected gravity component, and fuse the long-wave corrected gravity component and the short-wave gravity component to obtain the gravity correction value.

Images