Gravity anomaly calculation method, device, equipment, storage medium and product

Abstract

The application discloses a gravity anomaly calculation method, device, equipment, storage medium and product, relates to the satellite altimetry data processing technical field, and the gravity anomaly calculation method comprises: obtaining sea surface height data observed by an altimetry satellite; according to the sea surface height data, a plurality of residual geoid gradients in different directions are calculated; the plurality of residual geoid gradients in different directions are fused by a Bayesian triangular cap method, and an east-west direction component and a south-north direction component of the residual geoid gradient are obtained; and based on the east-west direction component and the south-north direction component of the residual geoid gradient, the marine gravity anomaly of a target region is calculated by a gravity anomaly inversion algorithm. The application solves the technical problems of low utilization rate of multi-direction gradient information of a single altimetry satellite and insufficient component extraction precision, realizes mining and optimization of multi-direction gradient information of a single satellite, and improves the inversion precision of the marine gravity anomaly.

Description

Technical Field

[0001] This application relates to the field of satellite altimetry data processing technology, and in particular to gravity anomaly calculation methods, devices, equipment, storage media and products. Background Technology

[0002] Ocean gravity anomalies are crucial foundational data for marine geophysical exploration, seafloor topography inversion, and geoid refinement. Related technologies primarily utilize shipborne gravity measurements, airborne gravity measurements, or satellite altimetry to obtain ocean gravity anomalies. Among these, satellite altimetry, with its advantages of wide coverage, short observation periods, and relatively low cost, has become the main method for acquiring global ocean gravity field information.

[0003] Traditional satellite altimetry gravity inversion methods typically rely on along-orbit geoid gradient data and perform integral calculations using the inverse Vining-Meinez formula. However, the orbits of traditional altimetry satellites (such as the Jason series and CryoSat-2) are mostly oriented north-south, resulting in a one-dimensional along-orbit distribution of observation data. This leads to inconsistencies in the accuracy of the calculated geoid gradient in the north-south and east-west directions, with the east-west component often exhibiting larger errors. To compensate for this deficiency, related technologies usually require fusing altimetry data from multiple satellites with different orbital inclinations and jointly solving the gradient components on a regular grid using statistical methods such as least-squares collocation. While this multi-source data fusion approach improves accuracy to some extent, the data processing is complex and dependent on spatiotemporal registration between different satellite data, making it prone to introducing additional systematic errors.

[0004] In recent years, the emergence of new-generation altimeter satellites (such as SWOT) has enabled two-dimensional planar observations of sea surface height, directly acquiring geoid gradient information in multiple directions within the same area. This has made it possible to retrieve ocean gravity anomalies with high precision using a single satellite. However, how to fully utilize the multi-directional gradient information in the two-dimensional observation data, effectively reduce the uncertainty between gradients in different directions, and extract reliable east-west and north-south components remains a pressing technical problem to be solved. Summary of the Invention

[0005] The main purpose of this application is to provide a gravity anomaly calculation method, apparatus, equipment, storage medium and product, which aims to solve the technical problems of low utilization rate of multi-directional gradient information of a single altimeter satellite and insufficient component extraction accuracy.

[0006] To achieve the above objectives, this application proposes a gravity anomaly calculation method, the method comprising:

[0007] Acquire sea surface height data from altimetry satellites;

[0008] Based on the sea level data, residual geoid gradients in multiple different directions were calculated;

[0009] The residual geoid gradients in multiple different directions are fused using the Bayesian triangular hat method to obtain the east-west and north-south components of the residual geoid gradient.

[0010] Based on the east-west and north-south components of the residual geoid gradient, the marine gravity anomaly of the target area is calculated using a gravity anomaly inversion algorithm.

[0011] In one embodiment, the step of calculating the residual geoid gradient in multiple different directions based on the sea level data includes:

[0012] The sea surface height observation grid is extracted from the sea surface height data, and the pre-constructed reference geoid model is subtracted from the observation grid to obtain the residual geoid;

[0013] Using each grid point in the observation grid as a calculation center, determine the adjacent grid points around the calculation center in multiple preset directions;

[0014] For each preset direction, calculate the residual geoid height difference between the calculation center and the adjacent grid points in the corresponding direction;

[0015] Determine the spherical distance between the calculation center and the adjacent grid points;

[0016] The residual geoid gradient at the calculation center in this direction is obtained based on the ratio of the residual geoid height difference to the spherical distance.

[0017] In one embodiment, the step of determining adjacent grid points in multiple preset directions around the computing center includes:

[0018] Based on the location of the computing center in the observation grid, identify the surrounding grid points in the observation grid that are adjacent to the computing center;

[0019] Based on the geometric orientation of the surrounding grid points relative to the calculation center, the surrounding grid points are divided into different directional categories, and each directional category corresponds to a preset direction;

[0020] In each preset direction, the nearest surrounding grid point to the computing center is selected as the adjacent grid point in that direction.

[0021] In one embodiment, the step of fusing the residual geoid gradients in multiple different directions using the Bayesian triangular hat method to obtain the east-west and north-south components of the residual geoid gradient includes:

[0022] Based on the residual geoid gradients in multiple different directions, the north-south and east-west components of the residual geoid gradient corresponding to each direction are calculated to obtain a multi-component dataset. Each component dataset contains the north-south and east-west components of the same spatial location under multiple observation periods.

[0023] The Bayesian triangular hat method was applied to perform uncertainty analysis on the multi-component dataset, and the error statistics corresponding to each component dataset were calculated.

[0024] Based on the error statistics, the fusion weights of each component dataset are determined using Bayesian theory.

[0025] The north-south and east-west components of the multi-component dataset are weighted and combined using the fusion weights to obtain the fused residual geoid gradient north-south and east-west components.

[0026] In one embodiment, the step of applying the Bayesian triangular hat method to perform uncertainty analysis on the multi-component dataset and calculating the error statistic corresponding to each component dataset includes:

[0027] From the multi-component quantity dataset, a set is designated as a reference dataset, and the difference sequence between the remaining component quantity datasets and the reference dataset is calculated to construct a difference matrix;

[0028] The sample covariance matrix is ​​calculated based on the difference matrix, and the covariance matrix to be estimated, which characterizes the data noise, is introduced to establish the mathematical relationship between the covariance matrix to be estimated and the sample covariance matrix.

[0029] An objective function is constructed based on optimization theory, wherein the objective function takes the diagonal elements of the covariance matrix to be estimated as unknown variables to be solved, and sets constraints to determine a unique solution.

[0030] Set initial values ​​for iteration, and iteratively solve the objective function using a numerical optimization algorithm. When the objective function converges to its minimum value, output the diagonal elements of the covariance matrix to be estimated as the error statistic.

[0031] In one embodiment, the step of calculating the ocean gravity anomaly of the target area using a gravity anomaly inversion algorithm includes:

[0032] Based on the principle of the inverse Vening-Meinesz formula, an integral relationship is constructed between the north-south and east-west components of the residual geoid gradient and the residual gravity anomaly.

[0033] The residual gravity anomaly field is obtained by performing a global surface integral on the north-south and east-west components of the residual geoid gradient using the integral relationship.

[0034] In the surface integration process, a one-dimensional fast Fourier transform technique is used to accelerate the integration operation in order to meet the data processing requirements of regular grid forms.

[0035] To address the singularity problem generated by the integral kernel function in the neighborhood of the calculation point, the innermost integration region is defined and a separate numerical calculation is performed using a planar approximation method to compensate and correct the residual gravity anomaly field.

[0036] The corrected residual gravity anomaly field is superimposed with the reference gravity anomaly field corresponding to the reference geoid model that was subtracted when calculating the residual geoid gradient to obtain the final ocean gravity anomaly field.

[0037] Shipborne gravity measurement data or publicly available gravity field model data are obtained as a verification benchmark. The accuracy of the ocean gravity anomaly field is evaluated through statistical comparative analysis to verify the effectiveness of the method.

[0038] Furthermore, to achieve the above objectives, this application also proposes a gravity anomaly calculation device, which includes:

[0039] The acquisition module is used to acquire sea surface height data observed by altimetry satellites;

[0040] The calculation module is used to calculate the residual geoid gradient in multiple different directions based on the sea level data.

[0041] The fusion module is used to fuse the residual geoid gradients in multiple different directions using the Bayesian triangular hat method to obtain the east-west and north-south components of the residual geoid gradient.

[0042] The calculation module is used to calculate the marine gravity anomaly of the target area based on the east-west and north-south components of the residual geoid gradient using a gravity anomaly inversion algorithm.

[0043] In addition, to achieve the above objectives, this application also proposes a gravity anomaly calculation device, the device comprising: a memory, a processor, and a computer program stored in the memory and executable on the processor, the computer program being configured to implement the steps of the gravity anomaly calculation method as described above.

[0044] In addition, to achieve the above objectives, this application also proposes a storage medium, which is a computer-readable storage medium, on which a computer program is stored, and when the computer program is executed by a processor, it implements the steps of the gravity anomaly calculation method described above.

[0045] In addition, to achieve the above objectives, this application also provides a computer program product, which includes a computer program that, when executed by a processor, implements the steps of the gravity anomaly calculation method described above.

[0046] One or more technical solutions proposed in this application have at least the following technical effects:

[0047] Compared to related technologies that utilize only a single-directional gradient or simply average multi-directional gradients, this application acquires sea surface height data observed by altimetry satellites; calculates residual geoid gradients in multiple different directions based on the sea surface height data; fuses these residual geoid gradients using the Bayesian triangular hat method to obtain the east-west and north-south components of the residual geoid gradient; and calculates the marine gravity anomaly of the target area based on the east-west and north-south components of the residual geoid gradient using a gravity anomaly inversion algorithm. It is understood that this application employs the Bayesian triangular hat method, an uncertainty assessment and fusion technique. When faced with multi-directional gradient data, the Bayesian triangular hat method quantifies the error characteristics of each directional gradient, thereby determining the fusion weights of the data in each direction, thus achieving an optimized combination of gradient information from each direction. Therefore, the Bayesian triangular hat method can make full use of the multi-directional information in the two-dimensional observation data provided by a single altimeter satellite, avoid the information loss caused by uneven directional accuracy in traditional methods, and thus significantly improve the extraction accuracy of the east-west and north-south components. Finally, it can complete the inversion of high-precision ocean gravity anomalies and solve the technical problems of low utilization rate of multi-directional gradient information and insufficient component extraction accuracy of a single altimeter satellite. Attached Figure Description

[0048] The accompanying drawings, which are incorporated in and form part of this specification, illustrate embodiments consistent with this application and, together with the description, serve to explain the principles of this application.

[0049] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, for those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0050] Figure 1 This is a flowchart illustrating an embodiment of the gravity anomaly calculation method of this application.

[0051] Figure 2 This is a flowchart illustrating Embodiment 2 of the gravity anomaly calculation method of this application;

[0052] Figure 3 This is a flowchart illustrating Embodiment 3 of the gravity anomaly calculation method of this application;

[0053] Figure 4 This is a schematic diagram of the process for inverting marine gravity anomalies based on the Bayesian triangular hat method and fusing SWOT altimetry satellite multi-directional geoid gradients in this application embodiment;

[0054] Figure 5 This is a schematic diagram of the SWOT eight-direction geoid gradient in the embodiments of this application;

[0055] Figure 6 This is a lollipop diagram showing the root mean square error of the gravity anomaly grid and the shipborne survey line in the embodiments of this application.

[0056] Figure 7 This is a schematic diagram of the module structure of the gravity anomaly calculation device according to an embodiment of this application;

[0057] Figure 8 This is a schematic diagram of the hardware operating environment involved in the gravity anomaly calculation method in this application.

[0058] The purpose, features, and advantages of this application will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0059] It should be understood that the specific embodiments described herein are merely illustrative of the technical solutions of this application and are not intended to limit this application.

[0060] To better understand the technical solution of this application, a detailed description will be provided below in conjunction with the accompanying drawings and specific implementation methods.

[0061] The main solution in this application embodiment is:

[0062] Acquire sea surface height data from altimetry satellites;

[0063] Based on the sea level data, residual geoid gradients in multiple different directions were calculated;

[0064] The residual geoid gradients in multiple different directions are fused using the Bayesian triangular hat method to obtain the east-west and north-south components of the residual geoid gradient.

[0065] Based on the east-west and north-south components of the residual geoid gradient, the marine gravity anomaly of the target area is calculated using a gravity anomaly inversion algorithm.

[0066] In this embodiment, the application uses a gravity anomaly calculation device as the execution subject. For ease of description, it will be referred to as "device" in detail below.

[0067] Because existing gravity anomaly inversion methods based on traditional altimetry satellite data are limited by the single satellite orbit and the one-dimensional distribution of observation data along the orbit, the calculated geoid gradient has inconsistent accuracy in the north-south and east-west directions, with a large error in the east-west component. Although related technologies have attempted to improve this problem by fusing data from multiple satellites with different orbital inclinations, the multi-source data fusion process is complex and relies on spatiotemporal registration between different satellite data, which can easily introduce additional systematic errors. While the new generation of wide-swath interferometric altimetry satellites can provide two-dimensional planar observation data containing gradient information in multiple directions, how to effectively utilize this multi-directional gradient information, reduce uncertainties in each direction, and extract reliable gradient components remains a pressing technical problem to be solved.

[0068] This application provides a solution that fully utilizes multi-directional gradient information provided by a single two-dimensional observation satellite during satellite altimetry data processing. By introducing the Bayesian triangular hat method, the error characteristics of gradients in each direction are quantitatively analyzed and optimized for fusion, thereby significantly improving the extraction accuracy of the east-west and north-south components of the residual geoid gradient, and ultimately enhancing the inversion quality of ocean gravity anomalies. Specifically, this application acquires sea surface height data from altimetry satellites and calculates multi-directional residual geoid gradients. The Bayesian triangular hat method is then used to assess the uncertainty of gradients in each direction. Based on the assessment results, fusion weights are determined, ultimately yielding high-precision gradient components and retrieving ocean gravity anomalies.

[0069] Based on this, the embodiments of this application provide a gravity anomaly calculation method, referring to... Figure 1 , Figure 1 This is a flowchart illustrating the first embodiment of the gravity anomaly calculation method of this application.

[0070] In this embodiment, the gravity anomaly calculation method includes steps S10 to S40:

[0071] Step S10: Obtain sea surface height data from altimetry satellite observations;

[0072] It should be noted that altimeter satellites refer to satellite platforms equipped with radar altimeters or other remote sensing equipment capable of measuring the sea surface height relative to a reference ellipsoid. These include, but are not limited to, wide-swath interferometric altimeter satellites such as SWOT and traditional along-orbit altimeter satellites such as the Jason series and CryoSat-2. Sea surface height data are instantaneous sea surface height values ​​obtained from satellite observations, typically stored in grid form in various levels of data products. For example, the Level 3 L3_LR_SSH product provides a 2-kilometer resolution sea surface height observation grid, and the data can be sourced from public data centers such as AVISO.

[0073] Understandably, obtaining high-quality sea surface elevation data is a fundamental step in gravity anomaly inversion. By collecting raw observation data from altimetry satellites, reliable raw input is provided for subsequent calculations of the geoid gradient, ensuring the data coverage and accuracy of the inversion process.

[0074] Step S20: Based on the sea level data, calculate the residual geoid gradient in multiple different directions;

[0075] It should be noted that the residual geoid gradient refers to the rate of change of the residual portion obtained after subtracting the reference geoid model from the original sea level height data. The reference geoid model includes, but is not limited to, publicly available Earth gravity field models such as EGM2008 and EGM2020, used to remove long-wave gravity signals to highlight local gravity anomalies. Multiple different directions refer to the directions formed by connecting adjacent grid points in different orientations, centered on each grid point in the observation grid. These typically include eight directions: east, west, south, north, northeast, southeast, southwest, and northwest. The residual geoid gradient in each direction is obtained by calculating the difference in residual geoid height between the center point and adjacent points, divided by the spherical distance between the two points.

[0076] Understandably, by calculating multi-directional gradients, two-dimensional grid observation data can be transformed into multi-angle gradient information reflecting local changes in the gravitational field. This fully leverages the spatial structure advantages of wide-swath altimetry satellites, providing diverse observation samples for subsequent fusion and extraction of reliable gradient components.

[0077] Step S30: The residual geoid gradients in multiple different directions are fused using the Bayesian triangular hat method to obtain the east-west and north-south components of the residual geoid gradient.

[0078] It should be noted that the Bayesian triangular hat method is a multi-source data uncertainty assessment and fusion method based on Bayesian theory. It is used to estimate the error statistics of each dataset and determine the optimal fusion weights from multiple sets of noisy observation data. Specifically, in this technical solution, the residual geoid gradient in each direction is first decomposed into north-south and east-west components, resulting in multiple datasets corresponding to different directions. Then, the Bayesian triangular hat method is used to perform uncertainty analysis on each dataset, calculating the error statistics for each dataset. Next, based on the error statistics and Bayesian theory, the fusion weights for each dataset are determined. Finally, the north-south and east-west components of each dataset are weighted and combined according to the fusion weights to obtain the fused north-south and east-west components of the residual geoid gradient.

[0079] Understandably, by fusing multi-directional gradients using the Bayesian triangular hat method, the uncertainty of gradient data in each direction can be quantified, and different weights can be assigned according to the magnitude of uncertainty. This suppresses the adverse effects of directions with large errors on the final components, significantly improves the extraction accuracy of the east-west and north-south components, and provides a more accurate gradient field input for gravity anomaly inversion.

[0080] Step S40: Based on the east-west and north-south components of the residual geoid gradient, the marine gravity anomaly of the target area is calculated using a gravity anomaly inversion algorithm.

[0081] It should be noted that the gravity anomaly inversion algorithm refers to the computational process of recovering the gravity anomaly field from gradient data by utilizing the physical relationship between the geoid gradient and gravity anomalies, through methods such as integral transformation. Specifically, this includes, but is not limited to, establishing an integral relationship between gradient components and residual gravity anomalies based on the inverse Vening-Meinesz formula, obtaining the residual gravity anomaly field through global surface integration; using one-dimensional fast Fourier transform to accelerate computation during integration to adapt to regular grid data; delineating the innermost region and using a planar approximation method for separate compensation to address the singularity of the integral kernel function in the neighborhood of the calculation point; and finally, superimposing the residual gravity anomaly field with the reference gravity anomaly field corresponding to the reference geoid model subtracted during gradient calculation to obtain the final ocean gravity anomaly field. The target region refers to the specific sea area where ocean gravity anomalies need to be calculated, which can be global or a local sea area.

[0082] Understandably, by converting the fused high-precision gradient components into ocean gravity anomalies through gravity anomaly inversion algorithms, the gravity field information of the target area can be obtained. This information can then be used for applications such as marine geophysical research, seabed topography inversion, and resource exploration. This achieves a complete processing flow from satellite observation data to gravity anomaly products, verifying the effectiveness and reliability of this method.

[0083] This embodiment provides a gravity anomaly calculation method that employs the Bayesian triangular hat method, an uncertainty assessment and fusion technique. When faced with multi-directional gradient data, the Bayesian triangular hat method quantifies and analyzes the error characteristics of each directional gradient, thereby determining the fusion weights of the data in each direction and achieving an optimized combination of gradient information from each direction. Therefore, based on the Bayesian triangular hat method, the multi-directional information in the two-dimensional observation data provided by a single altimeter satellite can be fully utilized, avoiding information loss caused by uneven directional accuracy in traditional methods. This significantly improves the extraction accuracy of the east-west and north-south components, ultimately achieving high-precision inversion of ocean gravity anomalies and solving the technical problems of low utilization rate and insufficient component extraction accuracy of multi-directional gradient information from a single altimeter satellite.

[0084] In one feasible implementation, the step of calculating the residual geoid gradient in multiple different directions based on the sea level data includes:

[0085] The sea surface height observation grid is extracted from the sea surface height data, and the pre-constructed reference geoid model is subtracted from the observation grid to obtain the residual geoid;

[0086] Using each grid point in the observation grid as a calculation center, determine the adjacent grid points around the calculation center in multiple preset directions;

[0087] For each preset direction, calculate the residual geoid height difference between the calculation center and the adjacent grid points in the corresponding direction;

[0088] Determine the spherical distance between the calculation center and the adjacent grid points;

[0089] The residual geoid gradient at the calculation center in this direction is obtained based on the ratio of the residual geoid height difference to the spherical distance.

[0090] It should be noted that the sea surface height observation grid refers to the two-dimensional grid of sea surface height data arranged according to latitude and longitude rules in altimetry satellite data products. Each grid point corresponds to an observation value at a specific geographical location, and the grid spacing includes, but is not limited to, 2 kilometers or other resolutions. The calculation center refers to the reference grid point currently selected in the observation grid for calculating the residual geoid gradient. The gradient with surrounding grid points is calculated using this point as the center. Multiple preset directions refer to different pre-defined orientations around the calculation center, including but not limited to eight directions: east, west, south, north, northeast, southeast, southwest, and northwest, used to comprehensively capture anisotropic geoid changes. Adjacent grid points refer to grid points in the observation grid that are spatially adjacent to the calculation center, usually those directly adjacent to the calculation center in the latitude and longitude direction, and their selection is determined based on the preset directions. The residual geoid height difference refers to the difference in residual geoid height between the calculation center and adjacent grid points, reflecting the amount of change in the residual geoid between the two locations. Spherical distance refers to the shortest distance on the Earth's surface between the calculation center and adjacent grid points. It is calculated based on the latitude and longitude coordinates of the two points and is used to convert elevation differences into spatial rates of change. The residual geoid gradient refers to the rate of change of the residual geoid per unit spherical distance, specifically the residual geoid elevation difference in that direction divided by the spherical distance.

[0091] Understandably, by using each grid point as the computation center to calculate the residual geoid gradient in multiple preset directions, two-dimensional grid observation data can be transformed into multi-angle gradient information reflecting local gravity field changes. This process fully utilizes the spatial coverage advantage of wide-swath altimeter satellites, enabling the quantification of gradient values ​​in different directions. This provides a rich data foundation for subsequent Bayesian triangular cap fusion, thereby helping to improve the extraction accuracy of the east-west and north-south components and ensuring the reliability of gravity anomaly inversion results.

[0092] For example, refer to Figure 2 This is a schematic diagram illustrating the specific process of calculating residual geoid gradients in multiple directions based on the sea level data. First, a sea level observation grid is extracted from the acquired sea level data, and a pre-constructed reference geoid model is subtracted from the observation grid to obtain the residual geoid. Then, using each grid point in the observation grid as the calculation center, adjacent grid points in multiple preset directions are determined around this calculation center, such as east, west, south, north, and the four diagonal directions. For each preset direction, the residual geoid height difference between the calculation center and the adjacent grid points in the corresponding direction is calculated, and the spherical distance between the two points is also determined. Finally, the residual geoid gradient at the calculation center in that direction is obtained based on the ratio of the residual geoid height difference to the spherical distance. By traversing all grid points, multi-directional residual geoid gradient data covering the target area can be obtained.

[0093] 1. For example, the specific steps for calculating the residual geoid gradient in different directions using SWOT sea level data include:

[0094] Extract the 2 km sea surface height observation grid from the SWOT altimeter satellite Level 3 ocean data product L3_LR_SSH, and subtract the EGM2008 geoid model from it to extract the residual geoid;

[0095] Using SWOT 2D grid points as the basis for calculating the geoid gradient, and relating them to adjacent observation points, geoid gradients in multiple directions are calculated, such as... Figure 2 As shown, the geoid gradient The calculation formula is:

[0096]

[0097] in, N Q and N P For two consecutive grid points Q and P The geoid is high. d PQ for P Point and Q The spherical distance of a point.

[0098] In one feasible implementation, the step of determining adjacent grid points in multiple preset directions around the computing center includes:

[0099] Based on the location of the computing center in the observation grid, identify the surrounding grid points in the observation grid that are adjacent to the computing center;

[0100] Based on the geometric orientation of the surrounding grid points relative to the calculation center, the surrounding grid points are divided into different directional categories, and each directional category corresponds to a preset direction;

[0101] In each preset direction, the nearest surrounding grid point to the computing center is selected as the adjacent grid point in that direction.

[0102] It should be noted that the location of the computing center within the observation grid refers to the row and column index or geographic coordinates of that grid point in the latitude and longitude coordinate system, serving as a spatial reference for identifying surrounding grid points. Adjacent surrounding grid points are those grid points in the observation grid that are adjacent to the computing center in the row or column direction, including but not limited to grid points in the four directions (up, down, left, right) and the four diagonal directions, forming an eight-neighbor or four-neighbor adjacency relationship centered on the computing center. The geometric orientation of a surrounding grid point relative to the computing center refers to its azimuth angle relative to the computing center. For example, true north corresponds to an azimuth angle of 0 degrees, and true east corresponds to an azimuth angle of 90 degrees. The direction of a grid point can be determined through its geometric orientation. Direction category refers to the classification and grouping of surrounding grid points based on their geometric orientation. Grid points within the same azimuth interval are grouped into the same category, and each category corresponds to a preset direction. The nearest surrounding grid point in space refers to the surrounding grid point with the smallest actual distance to the computation center on the sphere within the same directional category. When there are multiple grid points in the same direction, the closest one is selected to ensure the locality and accuracy of gradient calculation.

[0103] Understandably, by identifying surrounding grid points based on the computation center and classifying them into directional categories according to their geometric orientation, and then selecting the nearest grid point in each direction as the neighbor, a standardized and unified neighborhood selection rule can be established. This rule ensures that each grid point has definite neighboring points available for computation in each preset direction, avoiding inconsistencies in gradient calculations caused by directional ambiguity or multi-point competition. This lays a reliable spatial sampling foundation for the accurate calculation of subsequent multi-directional residual geoid gradients, guaranteeing the integrity and directional representativeness of the gradient field data.

[0104] In one feasible implementation, the step of fusing the residual geoid gradients in multiple different directions using the Bayesian triangular hat method to obtain the east-west and north-south components of the residual geoid gradient includes:

[0105] Based on the residual geoid gradients in multiple different directions, the north-south and east-west components of the residual geoid gradient corresponding to each direction are calculated to obtain a multi-component dataset. Each component dataset contains the north-south and east-west components of the same spatial location under multiple observation periods.

[0106] The Bayesian triangular hat method was applied to perform uncertainty analysis on the multi-component dataset, and the error statistics corresponding to each component dataset were calculated.

[0107] Based on the error statistics, the fusion weights of each component dataset are determined using Bayesian theory.

[0108] The north-south and east-west components of the multi-component dataset are weighted and combined using the fusion weights to obtain the fused residual geoid gradient north-south and east-west components.

[0109] It should be noted that the north-south component of the residual geoid gradient for each direction refers to the projected component obtained by decomposing the gradient vector in a specific direction to the north-south direction (i.e., the meridian direction), reflecting the contribution of the gradient in that direction to the north-south variation. The east-west component of the residual geoid gradient for each direction refers to the projected component obtained by decomposing the same gradient vector to the east-west direction (i.e., the latitudinal direction), reflecting the contribution of the gradient in that direction to the east-west variation. A multi-component dataset refers to a series of data sets obtained after decomposing gradients in different directions. Each direction constitutes a group, and each group contains the north-south and east-west components of data for multiple observation periods under that direction. Multiple observation periods refer to multiple observation data obtained by altimeter satellites repeatedly passing through the same geographical location at different times. The time span and number of periods can be determined according to the actual situation of the data product, including but not limited to the number of days of the satellite's revisit cycle or the scientific mission phase. The same spatial location refers to grid points with the same geographical coordinates in the observation grid, ensuring that data from different periods and directions can be spatially aligned. The north-south component refers to the projected value of the gravity field change in the north-south direction, and the east-west component refers to the projected value of the gravity field change in the east-west direction. The Bayesian triangular hat method, applied here, refers to using this method to statistically infer the error characteristics of each component dataset. Specifically, it involves constructing the difference matrix and covariance matrix between datasets and introducing optimization to obtain the error estimates for each dataset. Uncertainty analysis refers to quantifying the combined impact of random and systematic errors in each dataset using mathematical methods to obtain indicators reflecting data reliability. Error statistics refer to the error values ​​of each dataset calculated through uncertainty analysis, including but not limited to standard deviation and variance, used to characterize the noise level of the data. Bayesian theory is a statistical inference framework based on prior information and the probability distribution of sample data updates, used in this step to calculate the fusion weights based on the error statistics. Fusion weights are weighting coefficients assigned to each component dataset; the weight is inversely proportional to the error statistics of that group of data, with smaller errors resulting in larger weights. Weighted combination refers to linearly superimposing the north-south and east-west components of each dataset according to the fusion weights, resulting in fused north-south and east-west components. The fused north-south component is the result of a weighted sum of the north-south components of each group, and the fused east-west component is the result of a weighted sum of the east-west components of each group.

[0110] Understandably, by decomposing the residual geoid gradient in multiple directions into north-south and east-west components, a multi-component dataset is constructed. The Bayesian triangular hat method is then introduced for uncertainty analysis and weight calculation, ultimately achieving optimized fusion of gradient information from each direction. This process fully utilizes redundant information in multi-directional observation data, quantifies the reliability differences between data in each direction, and assigns different weights based on error magnitude. This effectively suppresses the interference of noisy directions on the fusion results, significantly improving the extraction accuracy of the north-south and east-west components. This provides a high-quality basic gradient field input for subsequent gravity anomaly inversion based on the inverse Vening-Meinesz formula, ensuring the accuracy and stability of the final inversion results.

[0111] For example, the residual geoid gradients of multiple periods and directions at the same location are first calculated using SWOT multi-directional residual geoid gradients, forming multiple datasets.

[0112] Set any reference dataset, and the rest n -1 The difference between the dataset and the reference dataset constitutes m ×( n -1) order matrix Y .set up Y The covariance matrix is ​​( n -1)×( n -1) order matrix S And introduce unknown noise. n × n order covariance matrix R and( n -1)× n 1-th order matrix J satisfy covariance matrix R for:

[0113]

[0114] in, It is a dataset and Covariance between them, diagonal elements These are the unknowns that need to be calculated; to obtain a good estimate of the covariance matrix, the dataset length... m It must be much larger than the number of data sets. n However, at this point, the number of unknowns exceeds the number of equations, and they must be determined reasonably. n Only by having a few "free" parameters can a unique solution be obtained;

[0115] At this point, the Kuhn-Tak condition is introduced to solve the unique solution problem, and an objective function is introduced. F :

[0116]

[0117] in ; Set constraints and initial conditions, and calculate the diagonal elements by iteratively minimizing the objective function value. The error of each dataset is obtained. According to Bayesian theory, the weight of each dataset is defined as the normalized form of the cumulative product of the uncertainties of the remaining datasets after removing that dataset:

[0118]

[0119] In the formula, the numerator represents the removal of the first... i After the first dataset, the cumulative product of the uncertainties of all other datasets is calculated, with the normalization factor in the denominator, to ensure that the sum of all weights is 1; the multi-directional residual geoid gradient components are fused according to the weights of each dataset.

[0120] In one feasible implementation, the step of applying the Bayesian triangular hat method to perform uncertainty analysis on the multi-component dataset and calculating the error statistic corresponding to each component dataset includes:

[0121] From the multi-component quantity dataset, a set is designated as a reference dataset, and the difference sequence between the remaining component quantity datasets and the reference dataset is calculated to construct a difference matrix;

[0122] The sample covariance matrix is ​​calculated based on the difference matrix, and the covariance matrix to be estimated, which characterizes the data noise, is introduced to establish the mathematical relationship between the covariance matrix to be estimated and the sample covariance matrix.

[0123] An objective function is constructed based on optimization theory, wherein the objective function takes the diagonal elements of the covariance matrix to be estimated as unknown variables to be solved, and sets constraints to determine a unique solution.

[0124] Set initial values ​​for iteration, and iteratively solve the objective function using a numerical optimization algorithm. When the objective function converges to its minimum value, output the diagonal elements of the covariance matrix to be estimated as the error statistic.

[0125] It should be noted that the reference dataset refers to a set of data arbitrarily selected from the multi-component datasets, used as a benchmark for calculating the differences between other groups. The selection of the reference dataset does not affect the final error estimation result but does affect the convenience of the intermediate calculation process. The difference sequence refers to the data sequence obtained by subtracting the north-south and east-west components of each of the other component datasets from the reference dataset at the same spatial location and observation period, reflecting the relative deviation between data in different directions. The difference matrix is ​​a two-dimensional matrix composed of all difference sequences arranged according to dataset group and observation period. The rows of the matrix correspond to the observation period, and the columns correspond to the datasets other than the reference group. The sample covariance matrix is ​​a statistic calculated based on the difference matrix, used to describe the linear correlation and fluctuation characteristics between the datasets. The covariance matrix to be estimated, representing data noise, is an unknown covariance matrix. Its diagonal elements represent the true noise variance of each dataset, and the off-diagonal elements represent the noise covariance between different datasets; it is the target parameter that needs to be solved. The mathematical relationship refers to the linear constraint equation satisfied between the covariance matrix to be estimated and the sample covariance matrix, which links the observable sample statistics with the unknown parameters to be estimated. Optimization theory is the mathematical framework used to solve constrained extremum problems; in this step, it is used to construct an objective function with the diagonal elements of the unknown covariance matrix as variables. The objective function is a scalar function constructed based on optimization theory; its value reflects the degree of matching between the solution and the constraints, and minimizing the objective function corresponds to the optimal parameter estimation. The unknown variables to be solved are the diagonal elements of the covariance matrix to be estimated, i.e., the noise variance of each dataset, which are key parameters that need to be obtained through optimization calculations. Constraints are mathematical restrictions imposed to ensure the uniqueness and physical meaning of the solution, including but not limited to the positive definiteness of the covariance matrix and the non-negativity of the diagonal elements. Initial values ​​for iteration refer to the initial guesses assigned to the variables to be solved before the numerical optimization algorithm begins; reasonable initial values ​​help accelerate convergence and avoid getting trapped in local optima. Numerical optimization algorithms refer to computational methods used to iteratively find the minimum value of an objective function, including but not limited to gradient descent, Newton's method, and quasi-Newton methods. These algorithms gradually approach the optimal solution by continuously updating variable values. Convergence to the minimum value means that when the change in the objective function value is less than a preset threshold or the number of iterations reaches an upper limit, the algorithm is considered to have found the optimal solution that meets the accuracy requirements. Error statistics refer to the diagonal elements of the covariance matrix to be estimated in the final output, i.e., the estimated noise variance of each dataset. These statistics are used to quantify the error level of each data set and serve as the basis for subsequent weight calculations.

[0126] Understandably, by introducing a reference dataset to construct a difference matrix and calculating the sample covariance matrix, a mathematical relationship is established between the covariance matrix to be estimated and the sample covariance matrix. Then, optimization theory and numerical optimization algorithms are used to solve for the error statistics of each dataset. This process achieves rigorous quantification of the uncertainty in multi-component datasets. This solution method based on the Bayesian triangular hat method can fully utilize the statistical correlation between multiple sets of data, separating and estimating the errors of each set of data without the need for external truth references. This objectively evaluates the reliability of gradient data in different directions, providing a scientific basis for subsequently determining the fusion weights based on error statistics, and ensuring the rationality and accuracy of multi-directional gradient information fusion.

[0127] In one feasible implementation, the step of calculating the ocean gravity anomaly of the target area using a gravity anomaly inversion algorithm includes:

[0128] Based on the principle of the inverse Vening-Meinesz formula, an integral relationship is constructed between the north-south and east-west components of the residual geoid gradient and the residual gravity anomaly.

[0129] The residual gravity anomaly field is obtained by performing a global surface integral on the north-south and east-west components of the residual geoid gradient using the integral relationship.

[0130] In the surface integration process, a one-dimensional fast Fourier transform technique is used to accelerate the integration operation in order to meet the data processing requirements of regular grid forms.

[0131] To address the singularity problem generated by the integral kernel function in the neighborhood of the calculation point, the innermost integration region is defined and a separate numerical calculation is performed using a planar approximation method to compensate and correct the residual gravity anomaly field.

[0132] The corrected residual gravity anomaly field is superimposed with the reference gravity anomaly field corresponding to the reference geoid model that was subtracted when calculating the residual geoid gradient to obtain the final ocean gravity anomaly field.

[0133] Shipborne gravity measurement data or publicly available gravity field model data are obtained as a verification benchmark. The accuracy of the ocean gravity anomaly field is evaluated through statistical comparative analysis to verify the effectiveness of the method.

[0134] It should be noted that the inverse Vening-Meinesz formula is the theoretical foundation for describing the transformation relationship between the geoid gradient and gravity anomaly in physical geodesy. This principle establishes the integral transformation relationship between the gradient field and the anomaly field. The integral relationship refers to the mathematical expression constructed based on the inverse Vening-Meinesz formula, taking the north-south and east-west components of the residual geoid gradient as input, and outputting the residual gravity anomaly through integration operations including a kernel function. Global surface integration refers to integrating over all flow points within the entire computational region, accumulating the gradient component contribution at each flow point to the computational point, thus obtaining the complete residual gravity anomaly field. The residual gravity anomaly field refers to the spatial distribution of the remaining local gravity anomaly signal after removing the long-wavelength reference field. One-dimensional fast Fourier transform (FFT) technology refers to the method of converting spatial domain integration operations to frequency domain calculations, significantly reducing computational load by utilizing the periodicity and symmetry of regular grid data. The regular grid form refers to a two-dimensional matrix structure where data is arranged at equal latitude and longitude intervals, facilitating processing using fast algorithms such as Fourier transform. The integral kernel function is a weighting function in the integral formula that describes the geometric relationship between the calculation point and the flow point. Its function value is related to the spherical distance between the two points. The singularity problem refers to the situation where the integral kernel function becomes infinitely large or undefined when the calculation point and the flow point coincide, making direct integration impossible. The innermost integration region refers to a small circular neighborhood centered on the calculation point. Special processing methods are needed within this region to avoid singularities. The planar approximation method involves treating the sphere as a plane within a small area for approximate calculation, and then performing separate numerical integration of the contribution of the innermost region using analytical formulas. Compensation correction combines the results of the separate calculation of the innermost region with the global integration results to compensate for signal loss caused by avoiding singularity regions. The reference gravity anomaly field refers to the gravity anomaly value corresponding to the reference geoid model subtracted when calculating the residual geoid gradient. It is usually obtained by expanding this model to a certain order. Superposition involves adding the residual gravity anomaly field to the reference gravity anomaly field to recover the final ocean gravity anomaly field containing complete frequency band information. Shipborne gravity measurement data refers to gravity anomaly values ​​measured by a gravimeter mounted on a ship during navigation at sea. It possesses high local accuracy but has limited coverage. Publicly available gravity field model data refers to high-order global gravity field models published internationally, such as EGM2008 and the EIGEN series, which can serve as a benchmark for accuracy verification over a wide range. Statistical comparative analysis involves quantitatively evaluating the consistency between the inversion results and the verification benchmark by calculating statistical indicators such as root mean square error, correlation coefficient, and standard deviation.

[0135] Understandably, by establishing an integral relationship through the inverse Vening-Meinesz formula and performing a global surface integral, the conversion from the fused high-precision gradient components to the residual gravity anomaly field was achieved. The use of one-dimensional fast Fourier transform technology accelerated the calculation, significantly improving processing efficiency for regular grid data. To address the integral singularity problem, the innermost region was delineated for separate compensation, avoiding signal loss caused by avoiding singularities in traditional methods. The residual anomaly field was superimposed with the reference field to recover the complete gravity anomaly spectrum. Finally, accuracy was verified using shipboard measured data or publicly available models, forming a complete technical closed loop from data input and core processing to result evaluation. This series of steps ensures that the final ocean gravity anomaly field possesses both theoretical rigor and practical reliability and verifiability.

[0136] For example, refer to Figure 3 This is a schematic diagram illustrating the specific process of inverting marine gravity anomalies based on the fused residual geoid gradient components and verifying their accuracy. First, based on the inverse Vening-Meinesz formula, an integral relationship is constructed between the north-south and east-west components of the residual geoid gradient and the residual gravity anomaly. The gradient components are then subjected to a global surface integral using this relationship to obtain the residual gravity anomaly field. During the integration process, a one-dimensional fast Fourier transform technique is used to accelerate the integration operation to accommodate the data processing requirements of a regular grid. To address the singularity problem generated by the integration kernel function in the neighborhood of the calculation point, the innermost integration region is delineated and a separate numerical calculation is performed using a planar approximation method to compensate and correct the residual gravity anomaly field. Subsequently, the corrected residual gravity anomaly field is superimposed with the reference gravity anomaly field corresponding to the reference geoid model subtracted during the calculation of the residual geoid gradient to obtain the final marine gravity anomaly field. Finally, shipboard gravity measurement data or publicly available gravity field model data are obtained as a validation benchmark. Statistical comparative analysis is used to evaluate the accuracy of the ocean gravity anomaly field, thereby verifying the effectiveness of the proposed method. This workflow fully demonstrates the entire process from gradient components to the final product and its validation.

[0137] For example, the IVM formula in kernel function form is:

[0138]

[0139] in, Indicates the calculation point. Indicates the flow point. yes The residual gravity anomaly at the point Indicates normal gravity. and They represent The north-south and east-west components of the residual geoid gradient at a point. express Click The azimuth of the point, For kernel functions, it can be derived from... , The spherical distance between two points is represented; the resulting fused multi-directional residual geoid gradient north-south and east-west components are both regular grid data, therefore, the 1D FFT method can be used to accelerate the convolution calculation and rigorously implement surface integrals.

[0140]

[0141] in, , , , and They are latitude and longitude, respectively. and These are regular grid data and In the latitudinal and longitudinal grid spacing, It is a one-dimensional fast Fourier transform operator; the surface integral can be used to calculate the residual gravity anomaly from the residual geoid gradient component. Combined with the influence of the innermost circle effect on the residual gravity anomaly, and the gravity anomaly corresponding to the reference field subtracted when calculating the residual geoid gradient (i.e., the gravity anomaly model expanded to order 2160 in EGM2008), the final ocean gravity anomaly can be obtained.

[0142] It should be noted that the above examples are only for understanding this application and do not constitute a limitation on the gravity anomaly calculation method of this application. Any simple transformations based on this technical concept are within the protection scope of this application.

[0143] like Figure 4 The diagram shown illustrates the process of ocean gravity anomaly inversion based on the Bayesian triangular hat method fused with SWOT altimeter satellite multi-directional geoid gradients in an embodiment of the present invention, including the following processing steps:

[0144] Step 1: Collect sea surface height data from the SWOT altimeter satellite Level 3 ocean data product L3_LR_SSH for different latitude regions around the world based on the AVISO data center.

[0145] Step 2: Calculate the residual geoid gradient in different directions using SWOT sea level data.

[0146] Step 2 includes the following specific steps:

[0147] (1) Extract the 2 km sea surface height observation grid from the SWOT altimeter satellite Level 3 ocean data product L3_LR_SSH, and subtract the EGM2008 geoid model from it to extract the residual geoid;

[0148] (2) Using the SWOT 2D grid points as the basis for calculating the residual geoid gradient, when solving for the residual geoid gradient of the middle red point, the eight nearby observation points can be used to calculate the residual geoid gradient in eight directions, such as... Figure 2 As shown, the geoid gradient The calculation formula is:

[0149]

[0150] in, N Q and N P For two consecutive grid points Q and P The geoid is high. d PQ for P Point and Q The spherical distance of a point.

[0151] Step 3: Use the BTCH method to fuse the multi-directional residual geoid gradient and calculate its east-west and north-south components.

[0152] Step 3 includes the following specific steps:

[0153] (1) Using the SWOT multi-directional residual geoid gradient, calculate the north-south and east-west components of the residual geoid gradient at the same location in multiple periods and directions, and form multiple datasets.

[0154] Specifically, taking a single-direction, single-cycle calculation as an example, the formulas for calculating the north-south and east-west components of the residual geoid gradient are as follows:

[0155]

[0156] In the formula, and The north-south and east-west components of the residual geoid gradient are given by the above formula. The residual geoid gradients in other directions and periods are calculated according to the above formula.

[0157] (2) Uncertainty in calculating the east-west and north-south components of the residual geoid gradient in each direction based on the Bayesian triangular hat method:

[0158] Specifically, n Direct solution of geoid gradient nDirectional geoid gradient, including north-south and east-west components. Based on SWOT data. m The observation periods constitute a dataset of geoid gradient components. That is, n There are 1 dataset, each dataset has 10 datasets. m Each observation has a time step of one period. It is determined by truth value Sum of error terms composition:

[0159]

[0160] First, define n -1 dataset and one reference dataset (from The difference between (any choice) is:

[0161]

[0162] It is 1-th order matrix, The covariance matrix is ​​set as follows:

[0163]

[0164] Introducing unknown noise n × n order covariance matrix R and( n -1)× n 1-th order matrix J :

[0165]

[0166] satisfy covariance matrix R for:

[0167]

[0168] in, It is a dataset and Covariance between them, diagonal elements These are the unknowns that need to be calculated; to obtain a good estimate of the covariance matrix, the dataset length... m It must be much larger than the dataset size n. According to... The following relationship can be obtained:

[0169]

[0170] However, at this point, the number of unknowns in the above equation exceeds the number of equations, and they must be determined reasonably. n Only with a few "free" parameters can a unique solution be obtained. At this point, the Kuhn-Tak condition is introduced to solve the unique solution problem, and an objective function is introduced. F :

[0171]

[0172] Its constraints are:

[0173]

[0174] In the formula , , The initial conditions for the iteration are:

[0175]

[0176] The diagonal elements are calculated by setting constraints and initial conditions and then iteratively minimizing the objective function value. The error of each dataset is obtained. ;

[0177] (3) Calculate the corresponding weights based on Bayesian theory and the uncertainty of each dataset:

[0178] Specifically, according to Bayesian theory, for each dataset... weight It is defined as the normalized form of the cumulative product of uncertainties in the remaining datasets after removing this dataset:

[0179]

[0180] In the formula, the numerator represents the removal of the first... i After the first dataset, the cumulative product of the uncertainties of all other datasets is calculated, with the normalization factor in the denominator, to ensure that the sum of all weights is 1; the multi-directional residual geoid gradient components are fused according to the weights of each dataset.

[0181] (4) Based on the weighted fusion of each dataset, the final residual geoid gradient direction components are obtained:

[0182] Based on the weights of each dataset calculated in the above process, the final residual geoid gradient east-west and north-south components are obtained by weighted averaging.

[0183] Step 4: Use the east-west and north-south components of the merged residual geoid gradient to invert the ocean gravity anomaly using the inverse Vining-Meinez (IVM) formula.

[0184] Step 4 includes the following specific steps:

[0185] (1) Using the inverse Vening Meinesz formula to invert the residual gravity anomaly:

[0186] The specific kernel function form of the IVM formula is as follows:

[0187]

[0188] in, Indicates the calculation point. Indicates the flow point. yes The residual gravity anomaly at the point Indicates normal gravity. and They represent The north-south and east-west components of the residual geoid gradient at a point. express Click The azimuth of the point, For kernel functions, it can be derived from... , The spherical distance between two points is represented.

[0189] The north-south and east-west components of the fused multi-directional residual geoid gradient obtained by the method described in step 3 are both regular grid data. Therefore, the 1D FFT method can be used to accelerate the convolution calculation and strictly implement surface integration.

[0190]

[0191] in, , , , and They are latitude and longitude, respectively. and These are regular grid data and In the latitudinal and longitudinal grid spacing, It is a one-dimensional fast Fourier transform operator; the surface integral can be used to calculate the residual gravity anomaly from the residual geoid gradient component.

[0192] (2) Calculation of the innermost circle effect:

[0193] when , When the spherical distance between two points approaches zero, the kernel function becomes singular, and the azimuth angle cannot be determined. To avoid singularities, gravity anomalies are typically not calculated in the innermost region. To avoid this practice causing loss of gravity anomalies, the gravity anomaly in the innermost region needs to be calculated separately. The innermost region is usually assumed to be a circular area. Due to its small size, the spherical surface can be approximated as a plane, and the calculation can be performed using the following formula:

[0194]

[0195] in, yes exist derivative in direction, yes exist derivative in direction, Construct a rectangular coordinate system for the region. The axis points east. The axis points north. It is the radius of the innermost circle.

[0196] (3) Restore the reference field gravity anomaly:

[0197] By combining the influence of the innermost circle effect on the residual gravity anomaly, and adding the gravity anomaly corresponding to the reference field subtracted when calculating the residual geoid gradient (i.e., the gravity anomaly model expanded to order 2160 in EGM2008), the final ocean gravity anomaly can be obtained.

[0198] This invention uses the Bayesian triangular hat method to calculate the uncertainty of the east-west and north-south components of the residual geoid gradient extracted from SWOT sea surface height data. Based on Bayesian theory and the uncertainty of each dataset, the corresponding weights are calculated and fused to obtain the final residual geoid gradient directional components. Finally, the inverse Vining-Meinez formula is used to invert the ocean gravity anomaly.

[0199] This application also provides a gravity anomaly calculation system, including:

[0200] Data module: Used to collect and store sea surface height data of the L3_LR_SSH Level 3 ocean data product from the SWOT altimeter satellite for global ocean areas;

[0201] Residual Geoid Calculation Module: Used to calculate the residual geoid gradient in different directions of the SWOT 2D sea surface grid;

[0202] Multi-directional geoid fusion module: used to fuse the east-west and north-south components of the geoid gradient calculation.

[0203] Ocean gravity anomaly inversion module: used to invert ocean gravity anomalies based on the IVM formula;

[0204] Comparison and verification module: used to evaluate the effectiveness and accuracy of the method for inverting marine gravity anomalies using shipborne survey lines and gravity field models.

[0205] For example, refer to Figure 5 This diagram illustrates the geoid gradient in eight directions from the SWOT altimeter satellite. Using a center point in a two-dimensional grid as an example, the diagram shows eight adjacent grid points located in the east, west, south, north, northeast, southeast, southwest, and northwest directions. The residual geoid gradient at the center point is obtained by calculating the difference in residual geoid height between the center point and its adjacent points in each direction, divided by the spherical distance. This diagram clearly illustrates the spatial sampling method for multi-directional gradient calculation in this application, demonstrating the multi-directional information extraction process from wide-swath altimeter satellite two-dimensional observation data.

[0206] For example, refer to Figure 6 The figure presents a lollipop plot of the root mean square error (RMSE) between the gravity anomaly grid and shipborne survey line data. The plot displays the statistical results of the RMSE between the inversion results of this application and the measured gravity anomaly values ​​on different shipborne survey lines, presented in bar or lollipop format. The horizontal axis represents the survey line number, and the vertical axis represents the RMSE value. This plot is used to verify the inversion accuracy of the method proposed in this application. Comparative analysis shows that the RMSE of each survey line is within a small range, indicating that the marine gravity anomaly inverted by the method proposed in this application matches the measured data well, demonstrating high reliability and practical value.

[0207] This application also provides a gravity anomaly calculation device, please refer to... Figure 7 The gravity anomaly calculation device includes:

[0208] Module 10 is used to acquire sea surface height data observed by altimetry satellites;

[0209] The calculation module 20 is used to calculate the residual geoid gradient in multiple different directions based on the sea level data.

[0210] The fusion module 30 is used to fuse the residual geoid gradients in multiple different directions using the Bayesian triangular hat method to obtain the east-west and north-south components of the residual geoid gradient.

[0211] The calculation module 40 is used to calculate the marine gravity anomaly of the target area based on the east-west and north-south components of the residual geoid gradient using a gravity anomaly inversion algorithm.

[0212] And / or, the gravity anomaly calculation device includes:

[0213] The first extraction model is used to extract the sea surface height observation grid from the sea surface height data, and subtract the pre-constructed reference geoid model from the observation grid to obtain the residual geoid;

[0214] The first determining model is used to determine adjacent grid points in multiple preset directions around each grid point in the observation grid as the calculation center;

[0215] The first calculation model is used to calculate the residual geoid height difference between the calculation center and the adjacent grid points in the corresponding direction for each preset direction.

[0216] The second determining model is used to determine the spherical distance between the computing center and the adjacent grid points;

[0217] The second calculation model is used to obtain the residual geoid gradient at the calculation center in the direction based on the ratio of the residual geoid height difference to the spherical distance.

[0218] And / or, the gravity anomaly calculation device includes:

[0219] The first identification model is used to identify the surrounding grid points in the observation grid that are adjacent to the computing center, based on the position of the computing center in the observation grid.

[0220] The first partitioning model is used to divide the surrounding grid points into different directional categories according to the geometric orientation of the surrounding grid points relative to the calculation center, with each directional category corresponding to a preset direction;

[0221] The first selection model is used to select the nearest surrounding grid point to the computing center in each preset direction as the adjacent grid point in that direction.

[0222] And / or, the gravity anomaly calculation device includes:

[0223] The third calculation model is used to calculate the north-south component and east-west component of the residual geoid gradient corresponding to each direction based on the residual geoid gradient in the multiple different directions, so as to obtain a multi-component dataset. Each component dataset contains the north-south component and east-west component of the same spatial location under multiple observation periods.

[0224] The fourth calculation model is used to perform uncertainty analysis on the multi-component dataset using the Bayesian triangular hat method, and to calculate the error statistics corresponding to each component dataset.

[0225] The fifth computational model is used to determine the fusion weights of each component dataset based on Bayesian theory, according to the error statistics.

[0226] The first combination model is used to weight and combine the north-south and east-west components in the multi-component dataset using the fusion weights to obtain the fused residual geoid gradient north-south and east-west components.

[0227] And / or, the gravity anomaly calculation device includes:

[0228] The first model is used to specify a set as a reference dataset from the multi-component quantity dataset, calculate the difference sequence between the remaining component quantity datasets and the reference dataset, and construct a difference matrix;

[0229] The first step is to establish a model to calculate the sample covariance matrix based on the difference matrix, and to introduce the covariance matrix to be estimated, which represents the data noise, and to establish the mathematical relationship between the covariance matrix to be estimated and the sample covariance matrix.

[0230] The first model is set up to construct an objective function based on optimization theory, wherein the objective function takes the diagonal elements of the covariance matrix to be estimated as unknown variables to be solved and sets constraints to determine a unique solution.

[0231] The first output model is used to set the initial values ​​for iteration, and to iteratively solve the objective function through a numerical optimization algorithm. When the objective function converges to the minimum value, the diagonal element values ​​of the covariance matrix to be estimated are output as the error statistics.

[0232] And / or, the gravity anomaly calculation device includes:

[0233] The second model is used to construct an integral relationship between the north-south and east-west components of the residual geoid gradient and the residual gravity anomaly, based on the principle of the inverse Vening-Meinesz formula.

[0234] The first integral model is used to perform global surface integration of the north-south and east-west components of the residual geoid gradient using the integral relationship to obtain the residual gravity anomaly field.

[0235] The first adaptive model is used to accelerate the integration operation by employing one-dimensional fast Fourier transform technology during the surface integration process, so as to adapt to the data processing requirements of regular grid form.

[0236] The first correction model is used to address the singularity problem generated by the integral kernel function in the neighborhood of the calculation point. It delineates the innermost integration region and uses a plane approximation method to perform separate numerical calculations to compensate and correct the residual gravity anomaly field.

[0237] The first superposition model is used to superimpose the corrected residual gravity anomaly field with the reference gravity anomaly field corresponding to the reference geoid model that was subtracted when calculating the residual geoid gradient, so as to obtain the final ocean gravity anomaly field.

[0238] The first verification model is used to obtain shipborne gravity measurement data or publicly available gravity field model data as a verification benchmark. The accuracy of the ocean gravity anomaly field is evaluated through statistical comparative analysis to verify the effectiveness of the method.

[0239] The gravity anomaly calculation device provided in this application, employing the gravity anomaly calculation method in the above embodiments, can solve the technical problems of low utilization rate of multi-directional gradient information and insufficient component extraction accuracy of a single altimeter satellite. Compared with the prior art, the beneficial effects of the gravity anomaly calculation device provided in this application are the same as those of the gravity anomaly calculation method provided in the above embodiments, and other technical features in the gravity anomaly calculation device are the same as those disclosed in the methods of the above embodiments, and will not be repeated here.

[0240] This application provides a gravity anomaly calculation device, which includes: at least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, and the instructions are executed by the at least one processor to enable the at least one processor to perform the gravity anomaly calculation method in the first embodiment described above.

[0241] The following is for reference. Figure 8 The diagram illustrates a structural schematic of a gravity anomaly calculation device suitable for implementing the embodiments of this application. The gravity anomaly calculation device in the embodiments of this application may include, but is not limited to, mobile terminals such as mobile phones, tablets, laptops, digital broadcast receivers, PDAs (Personal Digital Assistants), PMPs (Portable Media Players), in-vehicle terminals (e.g., in-vehicle navigation terminals), and fixed terminals such as digital televisions and desktop computers. Figure 8 The gravity anomaly calculation device shown is merely an example and should not impose any limitations on the functionality and scope of use of the embodiments of this application.

[0242] like Figure 8As shown, the gravity anomaly calculation device may include a processing unit 1001 (e.g., a central processing unit, a graphics processing unit, etc.), which can perform various appropriate actions and processes according to a program stored in a read-only memory (ROM) 1002 or a program loaded from a storage device 1003 into a random access memory (RAM) 1004. The RAM 1004 also stores various programs and data required for the operation of the gravity anomaly calculation device. The processing unit 1001, ROM 1002, and RAM 1004 are interconnected via a bus 1005. An input / output (I / O) interface 1006 is also connected to the bus. Typically, the following systems can be connected to the I / O interface 1006: input devices 1007 including, for example, a touchscreen, touchpad, keyboard, mouse, image sensor, microphone, accelerometer, gyroscope, etc.; output devices 1008 including, for example, a liquid crystal display (LCD), speaker, vibrator, etc.; storage devices 1003 including, for example, magnetic tape, hard disk, etc.; and communication devices 1009. Communication device 1009 allows the gravity anomaly calculation device to communicate wirelessly or wiredly with other devices to exchange data. Although the figure shows gravity anomaly calculation devices with various systems, it should be understood that implementation or possession of all the systems shown is not required. More or fewer systems may be implemented alternatively.

[0243] Specifically, according to the embodiments disclosed in this application, the processes described above with reference to the flowcharts can be implemented as computer software programs. For example, embodiments disclosed in this application include a computer program product comprising a computer program carried on a computer-readable medium, the computer program containing program code for performing the methods shown in the flowcharts. In such embodiments, the computer program can be downloaded and installed from a network via a communication device, or installed from storage device 1003, or installed from ROM 1002. When the computer program is executed by processing device 1001, it performs the functions defined in the methods of the embodiments disclosed in this application.

[0244] The gravity anomaly calculation device provided in this application, employing the gravity anomaly calculation method described in the above embodiments, can solve the technical problems of low utilization rate of multi-directional gradient information from a single altimeter satellite and insufficient component extraction accuracy. Compared with the prior art, the beneficial effects of the gravity anomaly calculation device provided in this application are the same as those of the gravity anomaly calculation method provided in the above embodiments, and other technical features of this gravity anomaly calculation device are the same as those disclosed in the method of the previous embodiment, and will not be repeated here.

[0245] It should be understood that the various parts disclosed in this application can be implemented using hardware, software, firmware, or a combination thereof. In the description of the above embodiments, specific features, structures, materials, or characteristics can be combined in any suitable manner in one or more embodiments or examples.

[0246] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

[0247] This application provides a computer-readable storage medium having computer-readable program instructions (i.e., a computer program) stored thereon, which are used to execute the gravity anomaly calculation method in the above embodiments.

[0248] The computer-readable storage medium provided in this application may be, for example, a USB flash drive, but is not limited to, electrical, magnetic, optical, electromagnetic, infrared, or semiconductor systems, devices, or any combination thereof. More specific examples of computer-readable storage media may include, but are not limited to: electrical connections having one or more wires, portable computer disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fiber, portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination thereof. In this embodiment, the computer-readable storage medium may be any tangible medium containing or storing a program that can be used by or in conjunction with an instruction execution system, system, or device. The program code contained on the computer-readable storage medium may be transmitted using any suitable medium, including but not limited to: wires, optical cables, RF (Radio Frequency), etc., or any suitable combination thereof.

[0249] The aforementioned computer-readable storage medium may be included in the gravity anomaly calculation device; or it may exist independently and not assembled into the gravity anomaly calculation device.

[0250] The aforementioned computer-readable storage medium carries one or more programs, which, when executed by the gravity anomaly calculation device, cause the gravity anomaly calculation device to:

[0251] Acquire sea surface height data from altimetry satellites;

[0252] Based on the sea level data, residual geoid gradients in multiple different directions were calculated;

[0253] The residual geoid gradients in multiple different directions are fused using the Bayesian triangular hat method to obtain the east-west and north-south components of the residual geoid gradient.

[0254] Based on the east-west and north-south components of the residual geoid gradient, the marine gravity anomaly of the target area is calculated using a gravity anomaly inversion algorithm.

[0255] Computer program code for performing the operations of this application can be written in one or more programming languages ​​or a combination thereof, including object-oriented programming languages ​​such as Java, Smalltalk, and C++, and conventional procedural programming languages ​​such as the "C" language or similar programming languages. The program code can be executed entirely on the user's computer, partially on the user's computer, as a standalone software package, partially on the user's computer and partially on a remote computer, or entirely on a remote computer or server. In cases involving remote computers, the remote computer can be connected to the user's computer via any type of network—including a Local Area Network (LAN) or a Wide Area Network (WAN)—or can be connected to an external computer (e.g., via the Internet using an Internet service provider).

[0256] The flowcharts and block diagrams in the accompanying drawings illustrate the architecture, functionality, and operation of possible implementations of systems, methods, and computer program products according to various embodiments of this application. In this regard, each block in a flowchart or block diagram may represent a module, segment, or portion of code containing one or more executable instructions for implementing a specified logical function. It should also be noted that in some alternative implementations, the functions indicated in the blocks may occur in a different order than those indicated in the drawings. For example, two consecutively indicated blocks may actually be executed substantially in parallel, and they may sometimes be executed in reverse order, depending on the functions involved. It should also be noted that each block in the block diagrams and / or flowcharts, and combinations of blocks in the block diagrams and / or flowcharts, can be implemented using a dedicated hardware-based system that performs the specified function or operation, or using a combination of dedicated hardware and computer instructions.

[0257] The modules described in the embodiments of this application can be implemented in software or hardware. The names of the modules do not necessarily limit the functionality of the unit itself.

[0258] The readable storage medium provided in this application is a computer-readable storage medium that stores computer-readable program instructions (i.e., a computer program) for executing the above-described gravity anomaly calculation method. This solves the technical problems of low utilization rate of multi-directional gradient information from a single altimeter satellite and insufficient component extraction accuracy. Compared with the prior art, the beneficial effects of the computer-readable storage medium provided in this application are the same as those of the gravity anomaly calculation method provided in the above embodiments, and will not be elaborated upon here.

[0259] This application also provides a computer program product, including a computer program that, when executed by a processor, implements the steps of the gravity anomaly calculation method described above.

[0260] The computer program product provided in this application can solve the technical problems of low utilization rate of multi-directional gradient information from a single altimeter satellite and insufficient component extraction accuracy. Compared with the prior art, the beneficial effects of the computer program product provided in this application are the same as those of the gravity anomaly calculation method provided in the above embodiments, and will not be repeated here.

[0261] All acquisition of signals, information, or actions in this application are carried out in compliance with the relevant data protection laws and policies of the country where the application is located, and with the authorization of the relevant device owner.

[0262] The above description is only a part of the embodiments of this application and does not limit the scope of protection of this application. All equivalent structural transformations made under the technical concept of this application and using the content of this application specification and drawings, or direct / indirect applications in other related technical fields, are included in the scope of protection of this application.

Claims

1. A method for calculating gravity anomalies, characterized in that, The method includes: Acquire sea surface height data from altimetry satellites; Based on the sea level data, residual geoid gradients in multiple different directions were calculated; The residual geoid gradients in multiple different directions are fused using the Bayesian triangular hat method to obtain the east-west and north-south components of the residual geoid gradient. Based on the east-west and north-south components of the residual geoid gradient, the marine gravity anomaly of the target area is calculated using a gravity anomaly inversion algorithm. The step of fusing the residual geoid gradients in multiple different directions using the Bayesian triangular hat method to obtain the east-west and north-south components of the residual geoid gradient includes: Based on the residual geoid gradients in multiple different directions, the north-south and east-west components of the residual geoid gradient corresponding to each direction are calculated to obtain a multi-component dataset. Each component dataset contains the north-south and east-west components of the same spatial location under multiple observation periods. The Bayesian triangular hat method was applied to perform uncertainty analysis on the multi-component dataset, and the error statistics corresponding to each component dataset were calculated. Based on the error statistics, the fusion weights of each component dataset are determined using Bayesian theory. The north-south and east-west components of the multi-component dataset are weighted and combined using the fusion weights to obtain the fused residual geoid gradient north-south and east-west components.

2. The method as described in claim 1, characterized in that, The step of calculating the residual geoid gradient in multiple different directions based on the sea level data includes: The sea surface height observation grid is extracted from the sea surface height data, and the pre-constructed reference geoid model is subtracted from the observation grid to obtain the residual geoid; Using each grid point in the observation grid as a calculation center, determine the adjacent grid points around the calculation center in multiple preset directions; For each preset direction, calculate the residual geoid height difference between the calculation center and the adjacent grid points in the corresponding direction; Determine the spherical distance between the calculation center and the adjacent grid points; The residual geoid gradient at the calculation center in this direction is obtained based on the ratio of the residual geoid height difference to the spherical distance.

3. The method as described in claim 2, characterized in that, The step of determining adjacent grid points in multiple preset directions around the computing center includes: Using the location of the computing center in the observation grid as a reference, identify the surrounding grid points in the observation grid that are adjacent to the computing center; Based on the geometric orientation of the surrounding grid points relative to the calculation center, the surrounding grid points are divided into different directional categories, and each directional category corresponds to a preset direction; In each preset direction, the surrounding grid point that is closest to the computing center in space is selected as the adjacent grid point in that direction.

4. The method as described in claim 1, characterized in that, The steps of applying the Bayesian triangular hat method to perform uncertainty analysis on the multi-component dataset and calculating the error statistics corresponding to each component dataset include: From the multi-component quantity dataset, a set is designated as a reference dataset, and the difference sequence between the remaining component quantity datasets and the reference dataset is calculated to construct a difference matrix; The sample covariance matrix is ​​calculated based on the difference matrix, and the covariance matrix to be estimated, which characterizes the data noise, is introduced to establish the mathematical relationship between the covariance matrix to be estimated and the sample covariance matrix. An objective function is constructed based on optimization theory, wherein the objective function takes the diagonal elements of the covariance matrix to be estimated as unknown variables to be solved, and sets constraints to determine a unique solution. Set initial values ​​for iteration, and iteratively solve the objective function using a numerical optimization algorithm. When the objective function converges to its minimum value, output the diagonal elements of the covariance matrix to be estimated as the error statistic.

5. The method as described in claim 1, characterized in that, The steps for calculating the ocean gravity anomaly of the target area using the gravity anomaly inversion algorithm include: Based on the principle of the inverse Vening-Meinesz formula, an integral relationship is constructed between the north-south and east-west components of the residual geoid gradient and the residual gravity anomaly. The residual gravity anomaly field is obtained by performing a global surface integral on the north-south and east-west components of the residual geoid gradient using the integral relationship. In the surface integration process, a one-dimensional fast Fourier transform technique is used to accelerate the integration operation in order to meet the data processing requirements of regular grid forms. To address the singularity problem generated by the integral kernel function in the neighborhood of the calculation point, the innermost integration region is defined and a separate numerical calculation is performed using a planar approximation method to compensate and correct the residual gravity anomaly field. The corrected residual gravity anomaly field is superimposed with the reference gravity anomaly field corresponding to the reference geoid model that was subtracted when calculating the residual geoid gradient to obtain the final ocean gravity anomaly field. Shipborne gravity measurement data or publicly available gravity field model data are obtained as a verification benchmark. The accuracy of the ocean gravity anomaly field is evaluated through statistical comparative analysis to verify the effectiveness of the method.

6. A gravity anomaly calculation device, characterized in that, The device includes: The acquisition module is used to acquire sea surface height data observed by altimetry satellites; The calculation module is used to calculate the residual geoid gradient in multiple different directions based on the sea level data. The fusion module is used to fuse the residual geoid gradients in multiple different directions using the Bayesian triangular hat method to obtain the east-west and north-south components of the residual geoid gradient. The calculation module is used to calculate the marine gravity anomaly of the target area based on the east-west and north-south components of the residual geoid gradient using a gravity anomaly inversion algorithm. The gravity anomaly calculation device is also used to achieve: Based on the residual geoid gradients in multiple different directions, the north-south and east-west components of the residual geoid gradient corresponding to each direction are calculated to obtain a multi-component dataset. Each component dataset contains the north-south and east-west components of the same spatial location under multiple observation periods. The Bayesian triangular hat method was applied to perform uncertainty analysis on the multi-component dataset, and the error statistics corresponding to each component dataset were calculated. Based on the error statistics, the fusion weights of each component dataset are determined using Bayesian theory. The north-south and east-west components of the multi-component dataset are weighted and combined using the fusion weights to obtain the fused residual geoid gradient north-south and east-west components.

7. A gravity anomaly calculation device, characterized in that, The device includes: a memory, a processor, and a computer program stored in the memory and executable on the processor, the computer program being configured to implement the steps of the gravity anomaly calculation method as described in any one of claims 1 to 6.

8. A storage medium, characterized in that, The storage medium is a computer-readable storage medium, and a computer program is stored on the storage medium. When the computer program is executed by a processor, it implements the steps of the gravity anomaly calculation method as described in any one of claims 1 to 5.

9. A computer program product, characterized in that, The computer program product includes a computer program that, when executed by a processor, implements the steps of the gravity anomaly calculation method as described in any one of claims 1 to 5.

Images