A gravity inversion method and system applicable to the crustal structure of solid planets

Abstract

This application belongs to the field of planetary geophysical methods, specifically disclosing a gravity inversion method and system applicable to the crustal structure of solid planets. The method includes: combining a planetary gravity field spherical harmonic model and a topographic spherical harmonic model to obtain discrete coordinate points of surface topography and discrete coordinates of gravity observation points; calculating the gravity signal and topographic signal at the crust-mantle interface using the finite amplitude method or spherical grid forward modeling; rotating the observation coordinate system and converting radial gravity anomalies into Bouguer gravity anomalies using angular projection; reconstructing the surface topography using discrete coordinate points in the topographic observation coordinate system to obtain a three-dimensional gravity inversion model; inputting the position coordinates of the observation points in the gravity observation coordinate system and the vertical gravity anomaly data into the three-dimensional gravity inversion model, and obtaining the spatial distribution of planetary density by iteratively updating the gravity inversion model. This application can extract the internal geological structure of the crust by combining planetary gravity field and topographic data.

Description

Technical Field

[0001] This application belongs to the field of planetary geophysical methods and technology, and more specifically, relates to a gravity inversion method and system applicable to the crustal structure of solid planets. Background Technology

[0002] Combining gravity and topography is a crucial method for studying planetary internal structure. Local gravity inversion can reveal hidden mass nodules or magma chambers solidified by ancient volcanic activity, thus providing geophysical evidence for planetary geological evolution. Furthermore, gravity inversion can clarify regional density structures, providing technical guidance for probe landing site selection. Therefore, establishing a three-dimensional gravity inversion system suitable for planetary crustal structure is extremely important.

[0003] Currently, planetary gravity three-dimensional inversion studies mainly rely on the direct processing of Bouguer gravity anomaly data, and the removal of gravity signals generated at the crust-mantle interface is not considered in the inversion process. The potential problems include: First, the Bouguer gravity anomaly data used may not be based on the latest gravity spherical harmonic model. The latest models typically contain more detailed information, have higher resolution, and are more conducive to new scientific discoveries. Second, topographic factors were not considered in the inversion. Unlike Earth, exoplanets do not have relatively accurate global density models, making the reference density selection during planetary Bouguer gravity anomaly correction selective. These corrected densities may differ from the actual density of local surfaces, resulting in incomplete correction of local topographic signals. However, the grid cells that allow for topographic density variations were not retained in the actual inversion. Third, the inversion of Bouguer gravity anomaly data did not consider that the crust-mantle interface is also a source of Bouguer gravity anomalies, and the latter actually accounts for a significant proportion. Fourth, in specific studies, polynomial fitting was used to extract gravity anomalies from local major structures, making background removal subjective and not adhering to the objective facts of the observational data. Fifth, the influence of curvature was not considered when studying large areas. This is because planetary gravity models generally have low resolution, and the study area is large, therefore the influence of planetary curvature cannot be ignored. Since local inversion is usually performed in a Cartesian coordinate system, the Bouguer gravity anomaly in the original coordinate system is a radial gravity anomaly, while the Cartesian coordinate system shows a vertical gravity anomaly. Both are directional derivatives of the gravity potential, but they do not actually coincide and cannot be completely equivalent. Furthermore, the surface of the inverted model is no longer a straight line, and the curvature effect should be considered. Therefore, the observed surface should be processed using curved surfaces. Summary of the Invention

[0004] To address the shortcomings of existing technologies, the purpose of this application is to provide a gravity inversion method and system applicable to the crustal structure of solid planets. This aims to solve the problem that existing planetary gravity inversion methods do not consider the bias introduced by the crust-mantle interface and curvature topography to the numerical observation of the gravity field when inverting Bouguer gravity anomaly data.

[0005] The first aspect of this application relates to a gravity inversion method applicable to the crustal structure of solid planets, comprising the following steps:

[0006] Step 1: Perform terrain sampling on the sphere and obtain the discrete coordinates of the terrain sampling points using a terrain spherical harmonic model;

[0007] Step 2: Determine the gravity observation reference surface, sample gravity observation points on the observation surface, and obtain free air gravity anomaly data and discrete coordinates of gravity observation points by combining the gravity spherical harmonic model;

[0008] Step 3: Solve for the gravity anomalies caused by topography and the crust-mantle interface using the finite amplitude method or spherical mesh forward modeling. Subtract the gravity anomalies caused by topography and the crust-mantle interface from the free air gravity anomalies to obtain the radial gravity anomalies inside the crust.

[0009] Step 4: Determine a new observation coordinate system by rotating the discrete coordinates of the terrain sampling points and the discrete coordinates of the gravity observation points; and obtain the vertical gravity anomaly data in the new coordinate system by angular projection of the radial gravity anomaly within the shell.

[0010] Step 5: Reconstruct the surface topography using discrete coordinate points in the new coordinate system to obtain a three-dimensional gravity inversion model;

[0011] Step 6: Input the coordinates of the gravity observation points and the vertical gravity anomaly data in the new coordinate system into the three-dimensional gravity inversion model, and obtain the density distribution inside the planetary shell by iteratively updating the three-dimensional gravity inversion model.

[0012] In some implementations, the method for rotating the discrete coordinates of the terrain sampling points and the discrete coordinates of the gravity observation points in step four specifically includes the following steps:

[0013] A three-dimensional rotation matrix is ​​constructed using the center point of the study area, the geocenter, and other arbitrary coordinate points;

[0014] The coordinates after three-dimensional rotation are obtained by multiplying the discrete coordinates of the terrain sampling points / gravity observation points with the three-dimensional rotation matrix;

[0015] Using the point with the maximum z-value after 3D rotation as a reference, the coordinates after 3D rotation are translated along the z-axis. The discrete coordinates of the terrain sampling points / gravity observation points after 3D rotation are projected onto the horizontal plane to obtain a set of two-dimensional coordinate points.

[0016] Construct a two-dimensional rotation matrix, multiply it by the two-dimensional coordinate point set, and obtain the two-dimensional coordinates after the plane is rotated;

[0017] By combining the two-dimensional coordinates after plane rotation with the corresponding z-axis values ​​after translation, a new observation coordinate system after rotation and projection is obtained.

[0018] In some implementations, a method for solving gravity anomalies caused by terrain is used, specifically:

[0019] The reference sphere corresponding to the reference radius is meshed to generate grid nodes with equal latitude and longitude spacing; based on the latitude and longitude range of the study area, latitude and longitude sampling is performed at equal intervals on the gravity observation surface to form the location information of gravity observation points;

[0020] The radial distance from the grid node to the Earth's center is obtained by combining the latitude and longitude information of the grid node with the topographic spherical harmonic model to form the planetary surface. Discrete spherical grid cells are formed by the planetary surface, the reference sphere, and the grid nodes. The latitude and longitude information of the center of the discrete spherical grid cells is statistically analyzed.

[0021] Set the shell correction density, and set the density of discrete spherical grid cells on the planetary surface that is lower than the reference sphere as the correction density negative value, and vice versa as positive value; carry out spherical coordinate system gravity forward modeling based on gravity observation point location information, discrete spherical grid cell density and center latitude and longitude to obtain the topographic gravity anomaly.

[0022] In some implementations, a forward modeling method using spherical meshes is employed to solve for the gravity anomaly caused by the crust-mantle interface, specifically as follows:

[0023] Based on the spherical harmonic model of the crust-mantle interface, the radial length set of discrete grid nodes generated by sampling from the reference sphere is calculated, and its minimum and maximum values ​​are statistically analyzed. Two reference spheres are then formed, with the crust-mantle interface embedded within each reference sphere. Spherical grid cells are formed from the three-layer interface and the latitude and longitude information of the grid nodes, and the latitude and longitude information of the center of each spherical grid cell is statistically analyzed.

[0024] To obtain the corrected reference density of the shell-mantle interface, the known mantle density of the spherical harmonic model at the shell-mantle interface is subtracted from the shell corrected density. The density of spherical grid cells above the shell-mantle interface is set to a negative value, otherwise it is set to a positive value. The gravity anomaly at the shell-mantle interface is calculated by combining the location information of gravity observation points, the density of spherical grid cells, and the latitude and longitude information of the center.

[0025] In some implementations, the finite amplitude method is used to solve for gravity anomalies caused by terrain, specifically as follows:

[0026] Using the spherical surface corresponding to the planet's average radius as a reference surface, and combining planetary mass information and topographic spherical harmonic model, the data are input into the topographic gravity anomaly spherical harmonic coefficient calculation function of the open-source spherical harmonic analysis tool; the relative deviation of the power spectrum of the spherical harmonic coefficients of adjacent orders is calculated by increasing the finite amplitude order from the 2nd order; when the deviation is lower than the preset value, the corresponding topographic gravity spherical harmonic coefficient is obtained.

[0027] By combining the location information of gravity observation points and the spherical harmonic coefficient of topographic gravity, gravity anomaly data caused by topography are obtained.

[0028] In some implementations, the finite amplitude method is used to solve for the gravity anomaly caused by the crust-mantle interface, specifically as follows:

[0029] The latitude and longitude range of the joint gravity observation points and the latitude and longitude coordinates of the center point of the study area were used to select a spherical cap, and the average radius within the spherical cap area was calculated by combining the crust-mantle interface model;

[0030] Using the spherical surface corresponding to the average radius within the spherical cap region as the reference surface, the planetary mass, crust-mantle interface correction density, and calculation order of the spherical cap region are input into the gravity anomaly spherical harmonic coefficient calculation function of the open-source spherical harmonic tool; the relative deviation of the power spectrum of the spherical harmonic coefficients of adjacent orders is calculated by increasing the calculation order sequentially from order 2. When the deviation is lower than the preset value, the gravity spherical harmonic coefficient corresponding to the crust-mantle interface is obtained.

[0031] By combining the gravitational spherical harmonic coefficients corresponding to the crust-mantle interface and the latitude and longitude information of the gravity observation points, the gravity anomaly data caused by the crust-mantle interface is calculated. In some implementations, step six specifically involves: inputting the location information of the gravity observation points in the new observation coordinate system and the vertical gravity anomaly data into the three-dimensional gravity inversion model; constructing the target functional using the zero-order minimum entropy constraint combined with the regularized inversion method; solving the target functional using the Gauss-Newton optimization algorithm; obtaining the update direction of the three-dimensional gravity inversion model; determining the optimal update step size using a linear search method; updating the current three-dimensional gravity inversion model; repeating the iteration until the normalized fitting difference of the inversion result reaches a preset threshold or the number of iterations reaches the maximum value; and then outputting the inversion density result.

[0032] The second aspect of this application relates to a gravity inversion system applicable to the crustal structure of solid planets, comprising:

[0033] The first observation coordinate system module performs terrain sampling on the sphere and obtains the discrete coordinates of the terrain sampling points by combining the terrain spherical harmonic model.

[0034] The second observation coordinate system module determines the gravity observation reference surface, samples gravity observation points on the observation surface, and obtains free air gravity anomaly data and discrete coordinates of gravity observation points by combining the gravity spherical harmonic model.

[0035] The terrain correction module is used to solve the gravity anomaly caused by the terrain using the finite amplitude method or spherical mesh forward modeling, and to subtract the gravity anomaly caused by the terrain from the gravity anomaly caused by the terrain in the free air gravity anomaly.

[0036] The crust-mantle interface correction module is used to solve the gravity anomaly caused by the crust-mantle interface by using the finite amplitude method or spherical mesh forward modeling. It subtracts the gravity anomaly caused by the crust-mantle from the free air gravity anomaly after removing the gravity anomaly caused by the terrain to obtain the radial gravity anomaly inside the crust.

[0037] The observation projection correction module is used to determine a new observation coordinate system by rotating the discrete coordinates of the terrain sampling points and the discrete coordinates of the gravity observation points; and to obtain the vertical gravity anomaly data in the new observation coordinate system by angular projection of the radial gravity anomaly inside the shell.

[0038] The spatial subdivision module is used to reconstruct the surface topography using surface coordinate points in the new observation coordinate system and obtain a three-dimensional gravity inversion model.

[0039] The 3D inversion module is used to input the coordinates of gravity observation points in the new coordinate system and the vertical gravity anomaly data into the 3D gravity inversion model. By iteratively updating the 3D gravity inversion model, the density distribution inside the planetary shell is obtained.

[0040] In some implementations, the observation projection correction module includes:

[0041] The 3D rotation matrix construction unit is used to construct a 3D rotation matrix using the center point of the study area, the geocenter, and other arbitrary coordinate points.

[0042] The three-dimensional coordinate rotation unit is used to multiply the discrete coordinates of terrain sampling points / gravity observation points with a three-dimensional rotation matrix to obtain the coordinates after three-dimensional rotation;

[0043] The coordinate projection unit is used to translate the three-dimensional rotated coordinates along the z-axis, using the maximum z-value point as a reference, so that the coordinates are zeroed, and project the discrete coordinates of the three-dimensional rotated terrain sampling points / gravity observation points onto the horizontal plane to obtain a two-dimensional coordinate point set;

[0044] A two-dimensional coordinate rotation unit is used to construct a two-dimensional rotation matrix, which is then multiplied by a set of two-dimensional coordinate points to obtain the two-dimensional coordinates after the plane has been rotated.

[0045] The observation coordinate system reconstruction unit is used to combine the two-dimensional coordinates after plane rotation with the corresponding translated z-axis values ​​to obtain a new observation coordinate system after rotation and projection.

[0046] In some implementations, the shell-mantle interface correction module includes a reference spherical construction unit, a second mesh generation unit, a shell-mantle density correction unit, and a shell-mantle interface gravity anomaly calculation unit.

[0047] The reference spherical building element is used to calculate the radial length set of the shell-mantle interface spherical harmonic model at the discrete grid nodes, forming the shell-mantle interface; the minimum and maximum radial lengths of the shell-mantle interface are statistically analyzed to form two reference spheres, with the shell-mantle interface embedded in the two reference spheres.

[0048] The second grid generation unit is used to form spherical grid units by taking the reference sphere and the crust-mantle interface according to the terrain sampling points, and to collect the latitude and longitude information of the center of each spherical grid unit.

[0049] The shell-mantle density correction element is used to obtain the shell-mantle interface correction reference density by subtracting the known mantle density and the shell correction density from the spherical harmonic model of the shell-mantle interface. The density of spherical grid elements below the shell-mantle interface is set to a negative value, and vice versa.

[0050] The crust-mantle gravity anomaly calculation unit is used to calculate the gravity anomaly corresponding to the crust-mantle interface by using the spherical grid cell density, center latitude and longitude information, and gravity observation point location information of the second grid generation unit;

[0051] The shell-mantle interface correction module includes: an average radius calculation unit, a second spherical harmonic coefficient calculation unit, and a second spherical harmonic coefficient-to-gravity anomaly calculation unit;

[0052] The average radius calculation unit is used to select a spherical cap by combining the latitude and longitude range of the gravity observation point and the latitude and longitude coordinates of the center point of the study area, and to calculate the average radius within the spherical cap area by combining the shell-mantle interface model;

[0053] The second spherical harmonic coefficient calculation unit uses the spherical surface corresponding to the average radius within the spherical cap region as the reference surface, and inputs the planetary mass, crust-mantle interface correction density, and calculation order of the spherical cap region into the gravity anomaly spherical harmonic coefficient calculation function of the open-source spherical harmonic tool; according to the calculation order increasing sequentially from the 2nd order, the relative deviation of the power spectrum of the spherical harmonic coefficients of adjacent orders is calculated; when the deviation is lower than the preset value, the gravity spherical harmonic coefficient corresponding to its crust-mantle interface is obtained.

[0054] The second spherical harmonic coefficient conversion gravity anomaly calculation unit combines the gravity spherical harmonic coefficient corresponding to the crust-mantle interface with the latitude and longitude information of the gravity observation point to obtain gravity anomaly data caused by the crust-mantle interface.

[0055] Thirdly, this application provides an electronic device, comprising: at least one memory for storing a program; and at least one processor for executing the program stored in the memory, wherein when the program stored in the memory is executed, the processor is configured to execute the method described in the first aspect or any possible implementation thereof.

[0056] Fourthly, this application provides a computer-readable storage medium storing a computer program that, when run on a processor, causes the processor to perform the method described in the first aspect or any possible implementation thereof.

[0057] Fifthly, this application provides a computer program product that, when run on a processor, causes the processor to perform the method described in the first aspect or any possible implementation thereof.

[0058] It is understood that the beneficial effects of the second to fifth aspects mentioned above can be found in the relevant descriptions in the first aspect mentioned above, and will not be repeated here.

[0059] Overall, the technical solutions conceived in this application have the following beneficial effects compared with the prior art:

[0060] This application employs a gravity inversion method that removes crust-mantle interface signals. This is because Bouguer gravity anomalies not only include the residual mass of subsurface geological bodies but also the influence of undulating interfaces with density differences. The latter may constitute a significant proportion of the total Bouguer gravity anomaly data, and directly inverting the Bouguer gravity anomaly data could lead to interfering interpretations of the intracranial structure. Compared to existing techniques, related studies do not consider crust-mantle interface signals and attempt to use polynomial fitting to remove background and locate the main anomaly signals. However, since crust-mantle interface signals and intracranial gravity signals are coupled, this subjective processing method cannot completely eliminate crust-mantle interface interference, as it does not rigorously examine the short-wavelength characteristics of the crust-mantle interface signal from a physical perspective; the wavelengths of the two signals may be similar.

[0061] This application employs curvature correction in gravity inversion. In large-scale regional studies, the conversion between spherical and rectangular coordinate systems leads to inequivalence between radial and vertical gravity anomalies, and significant projection errors occur when the model deviates from the region's center. Simple calculations show that ignoring surface curvature and using a planar approximation for modeling can introduce errors of up to 50% when the surface angle is large. Compared to existing techniques, related studies ignore the effects of curvature, directly treating the surface as a plane, and do not consider the systematic biases caused by coordinate transformation and geometric approximation. These factors limit the application of related techniques in large-scale regional studies.

[0062] This application also provides a three-dimensional gravity inversion method based on entropy constraints improved by the Softmax function. Compared with the prior art, the introduction of the Softmax function enables it to quickly distinguish between anomalous and non-anomalous regions during inversion. In particular, it causes the gradient change of anomalous boundaries to decay rapidly, thereby enhancing the focusing performance of the minimum entropy regularization constraint and making it easier to obtain focused inversion results, thus alleviating the problem of gravity inversion results being prone to divergence.

[0063] This application proposes a comprehensive workflow for inverting planetary internal structures from a gravitational spherical harmonic model. Compared to existing technologies, it does not rely on the acquisition of publicly available Bouguer gravity anomaly data for planets, as such data is often based on corrections to early gravitational spherical harmonic models. This application, however, is based on rigorous terrain correction, fundamentally analyzing the correction process for planetary Bouguer anomalies and providing two feasible approaches: the finite amplitude method and spherical mesh element modeling. It can directly use the latest gravitational spherical harmonic models and establish a custom observation coordinate system through coordinate rotation. It rigorously analyzes terrain undulations and strictly considers terrain in the inversion process, offering greater flexibility and providing effective technical support for subsequent research on planetary internal structures. Attached Figure Description

[0064] Figure 1 This is a flowchart of a gravity inversion method applicable to the crustal structure of a solid planet provided in this application embodiment.

[0065] Figure 2(a) is a map showing the original coordinate system distribution of gravity observation points in the Mars example area provided in the embodiments of this application.

[0066] Figure 2(b) is a three-dimensional coordinate distribution diagram after three-dimensional rotation provided in the embodiment of this application.

[0067] Figure 2(c) is a two-dimensional coordinate distribution diagram of the projection onto the horizontal plane after three-dimensional rotation provided in the embodiment of this application.

[0068] Figure 2(d) is a coordinate distribution diagram of the measuring points and terrain points after two-dimensional coordinate rotation provided in the embodiment of this application.

[0069] Figure 2(e) is a schematic diagram of the coordinate correction of measuring points on the horizontal plane provided in the embodiment of this application.

[0070] Figure 2(f) is a three-dimensional coordinate distribution map of the Martian volcanic region topography and gravity observation points after processing by the observation coordinate system module, provided in the embodiment of this application.

[0071] Figure 3 This is a free-air gravity anomaly map of the mountain study area provided in the embodiments of this application.

[0072] Figure 4 These are the power spectra of the topographic gravity spherical harmonic coefficients calculated by different orders of the finite amplitude method provided in the embodiments of this application;

[0073] Figure 5 This is a schematic diagram of a reference model for estimating the mass of the shell between the shell and mantle interface and the undulating terrain, provided in an embodiment of this application.

[0074] Figure 6 This is an inversion model of an example volcanic terrain on Mars based on tetrahedral mesh partitioning, provided in this application embodiment.

[0075] Figure 7(a) is a simple anomaly model based on actual volcanic terrain provided in the embodiments of this application.

[0076] Figure 7(b) is a schematic diagram of the smooth inversion method results provided in the embodiments of this application.

[0077] Figure 7(c) shows the inversion result obtained by the traditional zero-order minimum entropy constraint inversion method provided in the embodiments of this application.

[0078] Figure 7(d) shows the inversion results of the entropy constraint method based on the improved softmax function provided in the embodiments of this application. Detailed Implementation

[0079] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0080] In this application, the term "and / or" describes the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent three cases: A existing alone, A and B existing simultaneously, and B existing alone. In this application, the symbol " / " indicates that the related objects are in an "or" relationship, for example, A / B means A or B.

[0081] In this application, the terms “first” and “second” are used to distinguish different objects, rather than to describe a specific order of objects.

[0082] In the embodiments of this application, the terms "exemplary" or "for example" are used to indicate that something is an example, illustration, or description. Any embodiment or design that is described as "exemplary" or "for example" in the embodiments of this application should not be construed as being more preferred or advantageous than other embodiments or design. Specifically, the use of the terms "exemplary" or "for example" is intended to present the relevant concepts in a specific manner.

[0083] In the description of the embodiments in this application, unless otherwise stated, "multiple" means two or more.

[0084] The embodiments of this application are described below with reference to the accompanying drawings.

[0085] This application provides a gravity inversion method applicable to the crustal structure of solid planets. In practical studies, the final gravity inversion observation data is directly obtained from a spherical harmonic gravity model, without relying on specific publicly available Bouguer gravity anomaly data. Two feasible methods are provided for removing gravity signals from topography and the crust-mantle interface. Furthermore, the influence of curvature in large-scale studies is analyzed in detail, and rigorous numerical correction of the gravity field is achieved through mathematical derivation, eliminating projection errors that may arise after transforming from spherical coordinates to a local rectangular coordinate system. Finally, an improved entropy-constrained inversion method is presented, constructing entropy constraints using a softmax function to accelerate and achieve three-dimensional focused gravity inversion.

[0086] This application provides a gravity inversion method applicable to the crustal structure of solid planets, specifically including the following steps:

[0087] Step S1: Load the terrain spherical harmonic model, determine the degree interval of terrain sampling, and make it as dense as possible. It is recommended that the interval be no less than [insert value here]. The scale can be determined based on the research scale, and then the discrete coordinates of the corresponding terrain sampling points are calculated;

[0088] Step S2: Select the gravity spherical harmonic model required for the study, determine the appropriate observation reference surface through terrain sampling point coordinate analysis, then determine the sampling accuracy of the gravity field observation points on the reference surface according to the research requirements, and then calculate the free air gravity anomaly and discrete coordinates of the gravity observation points based on the gravity spherical harmonic model.

[0089] Step S3: Solve for the gravity signal caused by the terrain using the finite amplitude method or spherical mesh forward modeling, and remove the terrain gravity contribution from the free air gravity anomaly data;

[0090] Step S4: Load the shell-mantle interface model, obtain the gravity anomaly caused by the shell-mantle interface through the finite amplitude method or spherical mesh forward modeling, and further subtract the gravity contribution caused by the shell-mantle interface from the gravity anomaly data obtained in S3.

[0091] Step S5: Discretize the coordinates of the terrain sampling points and gravity observation points, and construct a new observation coordinate system by coordinate rotation; at the same time, obtain the vertical gravity anomaly data in the new coordinate system by angular projection of the gravity anomaly data obtained in S4.

[0092] Step S6: Import the discrete point coordinate data of the ground surface under the new observation coordinate system into the multi-scale simulation modeling software to reconstruct the ground surface topography, and perform appropriate grid subdivision to obtain a three-dimensional gravity inversion model.

[0093] Step S7: Input the coordinates of the gravity observation points in the new observation coordinate system and the gravity anomaly data obtained in S5 into the three-dimensional gravity inversion model. Construct the target functional using the zero-order minimum entropy constraint method combined with the regularized inversion method, and solve the target functional using the Gauss-Newton optimization algorithm to obtain the update direction of the three-dimensional gravity inversion model. Further combine the linear search method to determine the optimal update step size, and then update the current inversion model. If the data normalization fitting difference is lower than the preset threshold, the inversion terminates; if the preset threshold is not reached and the current iteration number is less than the set number, continue to execute the above steps to iteratively update the inversion model parameters until the optimal inversion result is obtained, and output the final density space distribution.

[0094] More preferably, step S5, which involves constructing the terrain observation coordinate system and the gravity observation coordinate system using a coordinate rotation method, specifically includes the following steps:

[0095] Step S5.1: Determine the precision of the spherical mesh division, i.e. the latitude and longitude grid spacing, and solve the radial length of the terrain at each node using the terrain spherical harmonic model; determine the original coordinate system, with the planetary "North Pole" point and the Earth's center as the z-axis, and the line connecting the 0-degree meridian and the Earth's center as the x-axis. At the same time, determine the y-axis direction through the z-axis and x-axis. Solve the three-dimensional coordinates of the corresponding terrain sampling points in this coordinate system based on the latitude, longitude and longitude and radial length of each node, and obtain the discrete coordinates of the terrain sampling points;

[0096] Step S5.2: Statistically calculate the maximum radial length of the terrain within the study area, select a radial length slightly larger than this value as the reference radius, and construct a gravity observation surface; determine the latitude and longitude grid spacing of the gravity observation points on the observation surface, that is, use latitude and longitude to represent the spatial position of each gravity observation point, and further convert the spherical coordinates of the gravity observation points into the three-dimensional coordinates of each observation point in the rectangular coordinate system to obtain the discrete coordinates of the gravity observation points;

[0097] The same rotational projection transformation is applied to the discrete coordinates of the terrain sampling points and gravity observation points. The following section uses the rotational projection of the discrete coordinates of the terrain sampling points as an example:

[0098] Step S5.3: If the center point of the study area Earth's center and arbitrary coordinate points (Not coinciding with the other two points), construct the coordinate system for the first rotation. The rotation coordinate matrix for this process is:

[0099]

[0100] in for:

[0101]

[0102]

[0103] in for:

[0104]

[0105]

[0106]

[0107] in for:

[0108]

[0109] Discretize the coordinates of the terrain sampling points and the rotation matrix Multiplying them together yields the three-dimensional coordinates after rotation.

[0110] Step S5.4: Calculate the maximum z value of the three-dimensional coordinates obtained in step S5.3, record this point as the reference point, and reset it to 0. This process is equivalent to translation in the z direction. Perform the same process on other coordinate points after three-dimensional rotation, and obtain new three-dimensional coordinates based on the translation method of the reference point. The projection of these coordinates onto the horizontal plane can form a two-dimensional coordinate point set.

[0111] Step S5.5: Rotate the two-dimensional coordinate point set, select any two points and add the origin of the coordinate system to reconstruct the x and y axes. Let the coordinates of any two points be denoted as follows: and Then the y-axis direction vector can be obtained. :

[0112]

[0113] Then combine transpose vector This further leads to the x-axis direction vector. :

[0114]

[0115] Then the two-dimensional rotation matrix for:

[0116]

[0117] Step S5.6: By comparing the two-dimensional coordinates obtained in S5.4 with... Multiplying them yields the rotated two-dimensional coordinates;

[0118] Furthermore, the planar coordinate distribution of the gravity observation points can be finely adjusted. By selecting a suitable reference origin, its coordinates can be close to or coincide with the minimum coordinate value of the gravity observation points. The two-dimensional coordinates of other points are processed based on the newly determined reference origin. This step does not affect the z-axis values ​​of all these points. At this point, the transformation processing of all two-dimensional coordinates is completed.

[0119] Step S5.7: Combine the two-dimensional coordinates after the plane is rotated with the z-axis values ​​of these points to obtain the three-dimensional coordinates of all points after rotation and projection.

[0120] More preferably, the method for obtaining free air gravity anomaly data in step S2 is as follows:

[0121] By using the loaded spherical harmonic model of gravity, combined with the spherical harmonic expansion of the free-air gravity anomaly and the latitude and longitude of the gravity observation point, the free-air gravity anomaly corresponding to the gravity observation point can be solved. Its spherical harmonic expansion is specifically as follows:

[0122]

[0123] in The spherical harmonic coefficient is the coefficient of gravity. It is a normalized spherical harmonic function; and These are latitude and longitude, respectively. It is the order; For the number of times; It is the gravitational constant; Planetary mass; The average radius of the planet; The radius of the unfolded shape;

[0124] To ensure efficient data processing, and to rapidly acquire free-air gravity anomalies using open-source spherical harmonic analysis tools, it is necessary to guarantee the number of target sampling points in the latitudinal direction. The maximum spherical harmonic order of the gravitational field satisfies ;

[0125] In the tools, select the built-in Gaussler-Gande quadrature method. When calculating gravity anomalies, select to remove the normal field. If the observation surface is a sphere, enter the calculation radius. If it is an ellipsoid, you also need to enter parameters such as angular velocity and ellipsoid degrees.

[0126] More preferably, step S3 involves calculating the gravity anomaly corresponding to the terrain:

[0127] When setting the correction density and calculating using the finite amplitude method, if the spherical harmonic function of the undulating interface is... , indicating the first Step, The core method for determining the spherical harmonic coefficients involves polynomial expansion of the undulating interface, approximating the actual interface through the accumulation of secondary terrain features. The greater the difference in interface undulation, the higher the required expansion order. The definition of the unfolded terrain is as follows:

[0128]

[0129] in, For the undulating interface function in the spatial domain;

[0130] The purpose of the finite amplitude method is to minimize the topographic curvature and meet the gravity accuracy constraints of the observation to obtain the final result. This method can be implemented using the corresponding calculation function in an open-source spherical harmonics tool. Input parameters include the reference surface radius, planetary mass, and calculation order n. When calculating topographic spherical harmonics, the sphere corresponding to the planet's average radius is usually chosen as the reference surface. To evaluate the calculation accuracy, n is increased sequentially from order 2 onwards in the experiment, and the corresponding gravity field spherical harmonics are calculated. If the power spectrum deviation between the (n-1)th and nth data is less than a preset value, the final calculation is complete. The specific principle of the conversion from topographic spherical harmonics to gravity spherical harmonics in this process is as follows:

[0131]

[0132] in, and Positive and negative spherical harmonic coefficients representing the contribution of terrain gravity; and Positive and negative spherical harmonic coefficients representing terrain; This represents a complex mapping function that considers secondary terrain; further utilizing... and It can calculate the terrain gravity anomaly; by subtracting this part of the data from the free air gravity, the gravity anomaly data after terrain removal can be obtained;

[0133] When using spherical mesh elements for 3D modeling, if the terrain sampling interval is... ,generally To ensure sufficient discretization accuracy. Radial gravity anomaly is considered when taking surface topography into account. The integral can be expressed as:

[0134]

[0135] in and This represents the length, latitude, and longitude of the center of the discrete grid cell; That is the unit density; It is the gravitational constant; , , These represent the upper and lower limits of the integration in the radial length, latitude, and longitude directions, respectively. This study proposes to use a bilinear interpolation method to approximate the two-dimensional surface of the Earth's surface, in order to reduce numerical errors that may arise from inaccuracies in the processing of the surface topography. Using the bilinear interpolation method, the relevant parameters can be further expanded:

[0136]

[0137]

[0138]

[0139] In the above formula Representing the upper or lower surface of a spherical mesh element, the topological transformation of the integration domain corresponding to the variable integral range is applied to the interval. The Jacobian matrix of the transformation is written as:

[0140]

[0141] The final gravity integral discretization is:

[0142]

[0143] in This module calculates the Gaussler-Gande quadrature coefficients. The specific process includes setting a reference sphere, gridding the surface, and statistically analyzing the location information of discrete grid cells; setting the terrain correction density and assigning it to the corresponding grid cells; combining gravity observation points and spherical grid cell information with the spherical grid forward modeling method to quickly calculate terrain gravity anomalies; and removing this portion of the data from the free-air gravity data to obtain the terrain-removed gravity anomaly data.

[0144] More preferably, the method for removing gravity anomalies corresponding to the crust-mantle interface is as follows:

[0145] When using the finite amplitude method (FAM) for calculation, a spherical harmonic model of the crust-mantle interface is loaded. Using the center point of the study area as a reference, a suitable spherical cap is selected to cover the latitude and longitude distribution range of the study area. Simultaneously, the average radius corresponding to the crust-mantle interface covered by the spherical cap is calculated. The crust-mantle interface corrected density is obtained by subtracting the known mantle density from the corrected crust density from the loaded crust-mantle interface model. The crust-mantle interface corrected density, planetary mass, and other parameters are then input into the FAM calculation tool to obtain the corresponding order of the crust-mantle interface gravitational spherical harmonic coefficients. In this specific study, the order is gradually increased from 2, and the power spectrum deviation between adjacent orders is calculated. If the required accuracy is met, the corresponding order of the crust-mantle interface gravitational spherical harmonic coefficients is selected.

[0146] By further combining the information from gravity observation points and the gravity spherical harmonic coefficient of the crust-mantle interface, the gravity anomaly caused by the crust-mantle interface can be calculated. Subtracting the former from the gravity anomaly value after removing the topography will yield the gravity observation value before projection.

[0147] When using spherical mesh elements for 3D modeling, the establishment of the reference surface includes: loading a spherical harmonic model of the crust-mantle interface, calculating the radial length values ​​of the model at discrete mesh nodes to form the crust-mantle interface, and statistically analyzing its minimum and maximum values, while simultaneously forming two reference spheres, which are distinguished by the crust-mantle interface; further, based on the latitude and longitude sampling pattern of the terrain sampling points, and combining the three-layer interface for discretization, corresponding mesh elements can be obtained. Density negative values ​​are assigned to each mesh element; the density of regions lower than the crust-mantle interface is set as the crust-mantle interface correction density negative value, and vice versa; then, combined with gravity observation point information, the gravity anomaly values ​​of the crust-mantle interface are obtained using the spherical mesh forward modeling method. Subtracting the former from the gravity anomaly values ​​after terrain removal also yields the gravity observation values ​​before projection.

[0148] After the original coordinate system is converted to a new coordinate system through the observation coordinate system module in step S5, the observed gravity anomaly values ​​need to be corrected synchronously, specifically including:

[0149] Based on the transformation method of the observation coordinate system module, the three-dimensional coordinates of the gravity observation point can be obtained. If the coordinates of any point within this system... For example, calculate the three-dimensional coordinates of the original coordinate system origin (geocenter) after processing by the observation coordinate system module, denoted as... Then the coordinates of any observation point Rotation coordinates from origin This will form a vector:

[0150]

[0151] Calculate vectors The cosine of the angle between the new coordinate system and the z-axis is as follows:

[0152]

[0153] If the coordinates of the observation point The corresponding gravitational field value is Then the projection correction The value is:

[0154]

[0155] The spatial subdivision module, relying on the three-dimensional coordinates of the terrain discrete points processed by the observation coordinate system module, reconstructs the three-dimensional terrain through simulation modeling software; imports the terrain file composed of the three-dimensional coordinates of the terrain discrete points into the software, constructs the inversion mesh, and uses tetrahedral mesh to discretize the inversion space to form a three-dimensional inversion model; in the gravity inversion, the density values ​​of the preset initial model are all set to 0.

[0156] The 3D inversion module constructs the 3D inversion objective function:

[0157]

[0158] in , , The regularization parameters are as follows:

[0159]

[0160]

[0161]

[0162]

[0163] in, , , , These are the data weighting matrix, model weighting matrix, roughness matrix, and entropy constraint term matrix, corresponding to the respective objective functionals. , , , The first two items are as follows:

[0164]

[0165] in, For the sensitivity matrix, The constant is between 0 and 1; the sensitivity matrix will be calculated by the following formula, if the... The coordinates of the observation points are , No. The centroid coordinates of each tetrahedral element are Then we have:

[0166]

[0167] Wherein, the roughness matrix represents the target element and its adjacent elements. N The mathematical mapping between units is essentially a smoothing filter, specifically:

[0168]

[0169] in, For the target cell volume, For the target unit, the first of The volume of each adjacent unit For the target unit and the first Euclidean distance between the geometric centers of adjacent units;

[0170] The entropy constraint term is specifically represented as follows:

[0171]

[0172] in, It is expressed as the sum of the absolute values ​​of the current cell values, i.e. , It is a very small constant;

[0173] Using the Gauss-Newton optimization algorithm, the inversion objective function can be expanded as follows:

[0174]

[0175]

[0176] in The reference model is typically the zero vector; the regularization parameters are updated in each iteration according to the following criteria:

[0177]

[0178] Each time the model parameters are updated When updating, use the Softmax function. Numerical values, which involve calculations of The model parameters will be corrected, and the corrected parameters will be denoted as... It is only used for updating Calculations are performed without participating in the inversion model parameter update; if the number of grid cells is... Specifically, there are:

[0179]

[0180] After determining the model update direction, the optimal step size is selected through linear search. Then the model update can be expressed as:

[0181]

[0182] The inversion stops when the misfit is lower than the preset value or the current number of inversions is lower than the maximum number.

[0183] Example 1

[0184] like Figure 1 As shown, this application provides a precise gravity inversion method applicable to planetary crustal structures, specifically including the following steps:

[0185] Step S1: Load the topographic spherical harmonic model, determine the degree interval of the topographic sampling, which should be more dense than that of gravity sampling, and can be determined according to the research scale. Then calculate the discrete coordinates of the corresponding sampling points.

[0186] This application establishes an observation coordinate system, taking into full account that the observation points are usually located on a sphere. This is because the unfolded reference surface of the spherical harmonic model basically corresponds to a sphere or ellipsoid. Taking Mars as an example, its corresponding reference model is a reference ellipsoid. Secondly, when studying local areas, the inversion model is converted into a rectangular coordinate system.

[0187] Step S2: Select the gravity spherical harmonic model required for the study, determine the observation reference surface through terrain data analysis, and determine the calculation accuracy according to the research requirements, that is, the degree interval of the gravity field at the sampling points on the spherical surface. Then, calculate the free air gravity anomaly corresponding to the reference surface based on the gravity spherical harmonic model.

[0188] This application takes the study of a Martian volcano as an example. The study area's center latitude and longitude is (174°E, 8.8°S), extending 5° to the left and right as the target study area. The latitude and longitude spacing of gravity observation points is 0.5°. The free-atmosphere gravity anomaly is calculated as follows: Figure 3 As shown.

[0189] Step S3: Select a suitable reference density, solve the gravity signal caused by the terrain using the finite amplitude method or spherical mesh forward modeling, and remove the terrain gravity contribution from the free air gravity anomaly data;

[0190] Taking the finite amplitude method as an example, when performing terrain correction for Mars, the maximum order of the known gravity spherical harmonic model is only 120, and a shell correction density of 2900 is selected. The average radius of Mars is 3389.5 km; the air density is approximately 0, therefore the density difference between the two sides of the surface is 2900; using the finite amplitude method, the spherical harmonic coefficients of topographic gravity corresponding to orders 2 to 10 are calculated. Figure 4 The calculated power spectrum is shown, where the power spectrum Defined as a spherical harmonic function of each order Total contribution:

[0191]

[0192] Figure 4The dashed line shows the first-order approximation results. When the finite amplitude order is greater than 1, there are significant numerical differences between short-wavelength gravity signals. Further calculation of the relative error between power spectra of adjacent orders shows that when the order is greater than 5, the relative error between adjacent orders is only 0.01%, indicating that selecting a finite amplitude order of 5 is sufficient to meet the accuracy requirements. Further combining the gravity spherical harmonic coefficient at order 5 and the gravity observation point, the topographic gravity anomaly can be calculated. Finally, subtracting the topographic gravity anomaly from the free air gravity anomaly data obtained in S3 yields the de-topographic gravity anomaly data.

[0193] Step S4: Load the crust-mantle interface model, obtain the gravity anomaly caused by the crust-mantle interface through the finite amplitude method or spherical mesh forward modeling, and further subtract the gravity contribution caused by the crust-mantle interface from the gravity anomaly data obtained in S3; unlike the calculation in topographic correction, the average radius of the study area must first be defined according to the crust-mantle interface model. Similar to the reference radius of the Earth's surface This part can be referenced. Figure 5 The diagram shows a schematic of the interface model. Combining the known mantle density and the corrected crust density from the crust-mantle interface model, the density difference at the crust-mantle interface can be obtained, which is the corrected crust-mantle interface density. This is then combined with the known planetary mass and average radius. The finite amplitude method is used to calculate the gravitational spherical harmonic coefficients corresponding to the shell-mantle interface. The calculation process is similar to that described in step S3. Finally, the gravity anomaly data obtained in step S4 is subtracted from the gravity anomaly data at the shell-mantle interface.

[0194] Step S5: Discretize the coordinates of the terrain sampling points and gravity observation points, determine the new observation coordinate system by coordinate rotation, and convert the gravity anomaly obtained in step S4 into vertical gravity anomaly data in the new coordinate system by angular projection;

[0195] The original coordinate system perspective is not conducive to local 3D inversion modeling research, so coordinate rotation transformation is essential. Figures 2(a) to 2(f) show the corresponding schematic diagrams of each step after coordinate rotation. Figure 2(a) shows the original coordinate distribution, which is obtained by transforming the terrain discrete points and gravity observation points through rectangular coordinates. In this process, the z-axis is the line connecting the North Pole and the Earth's center, and the y-axis is the line connecting the intersection of the 0° longitude and 0° latitude line and the Earth's center. Figure 2(b) shows the result after the first 3D coordinate rotation. Figure 2(c) is the coordinate projection of Figure 2(b) on the horizontal plane. Figure 2(d) shows the result after the 2D coordinate rotation. Figure 2(e) shows all the observation coordinate points corrected by the previous steps, and the result after translation based on the smallest coordinate point projected on the horizontal plane by selecting the set of gravity observation points that are close to or coincident with the previous points. Figure 2(f) shows the final observation coordinate system, and the 3D coordinate distribution of the Martian example volcanic terrain and gravity observation points after processing by the observation coordinate module.

[0196] Step S6: Import the terrain sampling point coordinate data obtained in the new coordinate system in step S5 into the simulation modeling software to reconstruct the surface terrain, and perform appropriate tetrahedral meshing to obtain a three-dimensional gravity inversion model.

[0197] The three-dimensional coordinates of the discrete terrain points obtained through coordinate rotation are imported into simulation modeling software to create a terrain interpolation function, thereby reconstructing the undulating terrain. In this embodiment, a Martian example volcano model is obtained through modeling as follows: Figure 6 As shown; the curvature effect is also considered, which can be reflected in the inversion model; since the thickness of the Martian crust varies from a few kilometers to 150 kilometers, the maximum inversion depth is set to 200 kilometers. Based on the above analysis, the basis for the final division of the inversion calculation domain is obtained; in this embodiment, the inversion region is discretized by a tetrahedral grid, which is beneficial for approximating the surface topography.

[0198] Step S7: Input the coordinates of the gravity observation point and the gravity anomaly data obtained in step S5 into the three-dimensional gravity inversion model. Construct the target functional using the zero-order minimum entropy constraint combined with the regularized inversion method, and solve the target functional using the Gauss-Newton optimization algorithm to obtain the model update direction. Further combine linear search to determine the optimal update step size and update the current inversion model. If the data normalization fitting difference is lower than the preset threshold or the maximum number of iterations is reached, the inversion terminates; otherwise, continue to execute the above iterative process until the optimal inversion result and output model density distribution are obtained.

[0199] In the inversion module, input the gravity observation point information. First, calculate the sensitivity matrix of the observation points corresponding to the discrete grid. The calculation applies the point mass approximation method, and the selected geometric center can be obtained through the coordinates of the four nodes of the tetrahedral element. The calculated mean value of the coordinates of all coordinate nodes is... ,Right now:

[0200]

[0201] A certain element in sensitivity That is, the element corresponding to the j-th measuring point of the i-th unit can be calculated by the following equation:

[0202]

[0203] Before inversion, set upper and lower limits for model inversion parameters. and Then, transform it to the logarithmic space, as follows:

[0204]

[0205] The objective functional is solved using the Gauss-Newton optimization method, and the model update direction at each iteration is calculated based on the expansion. Then, the optimal step size is determined; a set of step sizes is set up with a distribution from 0 to 1 and an interval of 0.01, and the corresponding step size and model update amount are compared. Multiplying yields the current model parameter update, and forward modeling is used to evaluate the computational data at this step size. and the original data The optimal step size is selected based on whether the data fit difference decreases. The data fit difference is defined as:

[0206]

[0207] The improvement to the entropy constraint term introduces the softmax function, aiming to fully utilize its classification ability to distinguish subtle differences and thus reasonably weight the target units in the entropy constraint matrix. For comparison, Figures 7(a) to 7(d) show the inversion results of different methods, which differ only in the difference of the target functional. Among them, the smooth inversion of the entropy-constrained positive term differs from the traditional entropy constraint method in that the improved softmax function is used to weight the target units. Improvements in computation.

[0208] Example 2

[0209] The second aspect of this application relates to a gravity inversion system applicable to the crustal structure of solid planets, comprising:

[0210] The first observation coordinate system module is used to determine the spacing between terrain sampling points and to obtain the discrete coordinates of terrain sampling points.

[0211] The second observation coordinate system module is used to determine the gravity observation reference surface and sample gravity observation points to obtain free air gravity anomaly data and discrete coordinates of gravity observation points.

[0212] The terrain correction module is used to solve the gravity anomaly caused by the terrain using the finite amplitude method or spherical mesh forward modeling, and to subtract the gravity anomaly caused by the terrain from the gravity anomaly caused by the terrain in the free air gravity anomaly.

[0213] The crust-mantle interface correction module is used to solve the gravity anomaly caused by the crust-mantle interface by using the finite amplitude method or spherical mesh forward modeling, and to subtract the gravity anomaly caused by the crust-mantle from the gravity anomaly caused by the crust-mantle interface in the free air gravity anomaly.

[0214] The observation projection correction module is used to determine a new observation coordinate system by rotating the discrete coordinates of terrain sampling points and gravity observation points; and to obtain vertical gravity anomaly data in a rectangular coordinate system by in-angle projection of radial gravity anomalies.

[0215] The spatial subdivision module is used to reconstruct the surface topography using discrete coordinate points of the surface under the new observation coordinate system and obtain a three-dimensional gravity inversion model.

[0216] The 3D inversion module is used to input the location coordinates of the observation points in the new observation coordinate system and the vertical gravity anomaly data into the 3D gravity inversion model. By iteratively updating the 3D gravity inversion model, the density distribution inside the planetary shell is obtained.

[0217] In some implementations, the observation projection correction module includes:

[0218] The 3D rotation matrix construction unit is used to construct a 3D rotation matrix using the center point of the study area (terrain / gravity), the geocenter, and other arbitrary coordinate points.

[0219] The three-dimensional coordinate rotation unit is used to multiply the discrete coordinates of the terrain sampling points / gravity observation points with the three-dimensional rotation matrix to obtain the three-dimensional coordinates after rotation.

[0220] The coordinate projection unit is used to translate the three-dimensional coordinates along the z-axis based on the point with the maximum z-value of the three-dimensional coordinates, and to project the discrete coordinates of the three-dimensional rotated terrain sampling points / gravity observation points onto the horizontal plane to obtain a set of two-dimensional coordinate points;

[0221] A two-dimensional coordinate rotation unit is used to construct a two-dimensional rotation matrix. By multiplying the two-dimensional coordinate point set with the two-dimensional rotation matrix, the two-dimensional coordinates after plane rotation are obtained.

[0222] The observation coordinate system reconstruction unit is used to combine the two-dimensional coordinates after plane rotation with the corresponding translated z-axis values ​​to obtain the topographic observation coordinate system and gravity observation coordinate system after rotation and projection.

[0223] In some implementations, the terrain correction module includes a first grid generation unit, a first gravity observation point establishment unit, and a terrain gravity anomaly calculation unit;

[0224] The first mesh generation unit is used to combine the terrain spherical harmonic model and spherical meshing to generate spherical mesh unit information;

[0225] The first gravity observation point establishment unit is used to sample measurement points on the gravity observation surface, form gravity observation point information, and obtain free air gravity anomaly data of the observation points by combining the gravity spherical harmonic model.

[0226] The terrain gravity anomaly calculation unit is used to combine gravity observation point information with spherical grid cell location information, use the spherical grid forward modeling method to calculate the terrain gravity anomaly, and subtract the terrain gravity anomaly from the free air gravity anomaly.

[0227] The terrain correction module may include a first spherical harmonic coefficient calculation unit and a first spherical harmonic coefficient converted gravity anomaly calculation unit;

[0228] The first spherical harmonic coefficient calculation unit is used to take the spherical surface corresponding to the planet's average radius as the reference surface, and input the planet's mass and finite amplitude order into the gravity spherical harmonic function of the open-source spherical harmonic tool; according to the calculation order increasing sequentially from the second order, it calculates and selects the gravity field spherical harmonic coefficients that meet the accuracy requirements.

[0229] The first spherical harmonic coefficient conversion gravity anomaly calculation unit is used to calculate and obtain gravity anomaly data caused by the terrain by combining the gravity spherical harmonic coefficient corresponding to the terrain and the latitude and longitude information of the gravity observation point.

[0230] In some implementations, the shell-mantle interface correction module includes a reference spherical construction unit, a second mesh generation unit, and a shell-mantle gravity anomaly calculation unit.

[0231] Reference spherical building blocks are used to calculate the radial length values ​​of the spherical harmonic model of the shell-mantle interface at discrete grid nodes, forming the shell-mantle interface; the minimum and maximum radial length values ​​of the interface are statistically analyzed to form two reference spheres; the interlayer formed by the reference spheres is distinguished by the shell-mantle interface.

[0232] The second grid generation unit is used to form discrete grid units from the system consisting of the reference sphere and the crust-mantle interface according to the terrain sampling points, and to collect the position information of each unit. The crust-mantle interface correction reference density is obtained by subtracting the known mantle density and the crust correction density from the crust-mantle interface spherical harmonic model. The density of the grid unit is set to negative if it is lower than that of the crust-mantle interface grid unit, and positive if it is lower. The above finally forms the density and position information of the second grid generation unit.

[0233] The crust-mantle gravity anomaly calculation unit is used to calculate the gravity anomaly corresponding to the crust-mantle interface by combining the density and location information of the second grid generation unit and the coordinate information of the gravity observation point;

[0234] The shell-mantle interface correction module includes: an average radius calculation unit, a second spherical harmonic coefficient calculation unit, and a second spherical harmonic coefficient-to-gravity anomaly calculation unit;

[0235] The average radius calculation unit is used to select a suitable spherical cap by combining the latitude and longitude range of the gravity observation point and the latitude and longitude coordinates of the center point of the study area, and to calculate the average radius within the spherical cap area by combining the crust-mantle interface model;

[0236] The second spherical harmonic coefficient calculation unit is used to take the spherical surface corresponding to the planet's average radius as the reference surface, and input the remaining planetary mass, crust-mantle interface correction density and calculation order into the gravity anomaly spherical harmonic coefficient calculation function of the open-source spherical harmonic tool. The relative deviation of the power spectrum of the spherical harmonic coefficients of adjacent orders is calculated according to the calculation order increasing sequentially from the second order. When the deviation is lower than the preset value, the accuracy requirement is met, and the gravity spherical harmonic coefficient corresponding to the crust-mantle interface is obtained.

[0237] The second spherical harmonic coefficient conversion gravity anomaly calculation unit will combine the gravity spherical harmonic coefficient corresponding to the crust-mantle interface and the latitude and longitude information of the gravity observation point to calculate and obtain gravity anomaly data caused by the crust-mantle interface.

[0238] In some implementations, the 3D inversion module is used to construct a target functional using the zero-order minimum entropy constraint method combined with a regularized inversion method, and solve the target functional using the Gauss-Newton optimization algorithm to obtain the update direction of the 3D gravity inversion model. The optimal update step size is determined by combining a linear search method, and the current 3D gravity inversion model is updated. The iteration is repeated until the normalized fitting difference of the inversion result reaches a preset threshold or the number of iterations reaches the maximum value, and the inversion density result is output.

[0239] In summary, this application has the following advantages compared with the prior art:

[0240] This application provides an accurate gravity inversion method applicable to the crustal structure of solid planets. This method does not rely on publicly available Bouguer gravity anomaly data. This part has always been a challenge in planetary research due to the accuracy of data processing and the analysis of the correction principle of planetary Bouguer gravity anomalies. In particular, most Bouguer gravity anomaly data are based on early publicly available spherical harmonic gravity field models, which often do not match the latest gravity field models. This application fundamentally analyzes the correction process of planetary Bouguer anomalies, which is more flexible, and provides two feasible approaches: the finite amplitude method and spherical mesh element modeling and solving, providing effective technical support for subsequent research on the internal structure of planetary crusts.

[0241] This application considers removing the crust-mantle interface signal to accurately obtain the internal structure of the crust. Surface-observed Bouguer gravity anomalies not only include the remaining mass of subsurface geological bodies but also the influence of undulating interfaces with density differences. Currently, when studying the internal structure of the crust, the Bouguer gravity anomaly inversion does not consider the influence of this crust-mantle interface signal, which accounts for a significant proportion in numerical calculations. Therefore, ignoring the influence of the crust-mantle interface signal may lead to interfering structural interpretations. In particular, some studies rely on polynomial fitting to remove the background, which involves significant subjectivity and fails to remove the crust-mantle interface signal exhibiting short-wavelength characteristics. Furthermore, this method of removing the gravity anomaly background differs significantly from traditional research, which typically considers a universally existing background, such as a uniform field or linear superposition, and usually handles it simply without affecting the subjective morphology of the target anomaly. Polynomial fitting, however, relies on polynomial coefficients for control, which may result in incomplete separation of local background fields.

[0242] This application considers the curvature effect in local planetary inversion studies, especially for studies of local regions. When the spherical coordinate system is transformed into a rectangular coordinate system, the radial gravity anomaly and the vertical gravity anomaly are not completely equivalent, especially in areas far from the center of the study area, where the error of the projection transformation is greater. Secondly, it also considers the changes in the surface morphology in the modeling. If the curvature effect is not considered and a plane is used instead of a curved surface for study, the larger the angle between the curved surface and the corresponding surface, the more likely it is to cause an error of nearly 50%, which may lead to incorrect inversion results.

[0243] This application also provides a three-dimensional gravity inversion method based on entropy constraints improved by the Softmax function. In the inversion, it can quickly distinguish between anomalous and non-anomalous regions. In particular, it can rapidly decay the gradient changes of anomalous boundaries, thereby enhancing the focusing performance of the minimum entropy regularization constraint and making it easy to obtain focused inversion results, thus alleviating the problem of gravity inversion results being prone to divergence.

[0244] This application addresses the shortcomings of existing methods by providing a precise gravity inversion method applicable to the structure of planetary crust, which will provide reliable geophysical evidence for the study of planetary geological evolution.

[0245] Based on the methods in the above embodiments, this application provides an electronic device that may include a processor, a communications interface, a memory, and a communication bus, wherein the processor, communications interface, and memory communicate with each other via the communication bus. The processor may invoke logical instructions stored in the memory to execute the methods in the above embodiments.

[0246] Furthermore, the logical instructions in the aforementioned memory can be implemented as software functional units and, when sold or used as independent products, can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or a portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of this application.

[0247] Based on the methods in the above embodiments, this application provides a computer-readable storage medium storing a computer program that, when run on a processor, causes the processor to execute the methods in the above embodiments.

[0248] Based on the methods in the above embodiments, this application provides a computer program product that, when run on a processor, causes the processor to execute the methods in the above embodiments.

[0249] Those skilled in the art will readily understand that the above description is merely a preferred embodiment of this application and is not intended to limit this application. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this application should be included within the scope of protection of this application.

Claims

1. A gravity inversion method applicable to the crustal structure of solid planets, characterized in that, Includes the following steps: Step 1: Determine the spacing between terrain sampling points and obtain the discrete coordinates of the terrain sampling points using the terrain spherical harmonic model; Step 2: Determine the gravity observation reference surface, sample gravity observation points on the observation reference surface, and obtain free air gravity anomaly data and discrete coordinates of gravity observation points by combining the gravity spherical harmonic model; Step 3: Solve for the gravity anomalies caused by topography and the crust-mantle interface using the finite amplitude method or spherical mesh forward modeling. Subtract the gravity anomalies caused by topography and the crust-mantle interface from the free air gravity anomalies to obtain the radial gravity anomalies inside the crust. Step 4: Determine a new observation coordinate system by rotating the discrete coordinates of the terrain sampling points and the discrete coordinates of the gravity observation points; and obtain the vertical gravity anomaly data in the new coordinate system by angular projection of the radial gravity anomaly within the shell. Step 5: Reconstruct the surface topography using discrete coordinate points in the new coordinate system to obtain a three-dimensional gravity inversion model; Step 6: Input the coordinates of the gravity observation points and the vertical gravity anomaly data in the new coordinate system into the three-dimensional gravity inversion model, and obtain the density distribution inside the planetary shell by iteratively updating the three-dimensional gravity inversion model.

2. The gravity inversion method according to claim 1, characterized in that, Step four involves rotating the discrete coordinates of the terrain sampling points and the gravity observation points. This process includes the following steps: A three-dimensional rotation matrix is ​​constructed using the center point of the study area, the geocenter, and other arbitrary coordinate points; The coordinates after three-dimensional rotation are obtained by multiplying the discrete coordinates of the terrain sampling points / gravity observation points with the three-dimensional rotation matrix; Using the point with the maximum z-value after 3D rotation as a reference, the coordinates after 3D rotation are translated along the z-axis. The discrete coordinates of the terrain sampling points / gravity observation points after 3D rotation are projected onto the horizontal plane to obtain a set of two-dimensional coordinate points. Construct a two-dimensional rotation matrix, multiply it by the two-dimensional coordinate point set, and obtain the two-dimensional coordinates after the plane is rotated; By combining the two-dimensional coordinates after plane rotation with the corresponding z-axis values ​​after translation, a new observation coordinate system after rotation and projection is obtained.

3. The gravity inversion method according to claim 1, characterized in that, The method for solving the gravity anomaly caused by terrain using spherical mesh forward modeling is as follows: The reference sphere corresponding to the reference radius is meshed to generate grid nodes with equal latitude and longitude spacing; based on the latitude and longitude range of the study area, latitude and longitude sampling is performed at equal intervals on the gravity observation surface to form the location information of gravity observation points; The radial distance from the grid node to the Earth's center is obtained by combining the latitude and longitude information of the grid node with the topographic spherical harmonic model to form the planetary surface. Discrete spherical grid cells are formed by the planetary surface, the reference sphere, and the grid nodes. The latitude and longitude information of the center of the discrete spherical grid cells is statistically analyzed. Set the shell correction density, and set the density of discrete spherical grid cells on the planetary surface that is lower than the reference sphere as the correction density negative value, and vice versa as positive value; carry out spherical coordinate system gravity forward modeling based on gravity observation point location information, discrete spherical grid cell density and center latitude and longitude to obtain the topographic gravity anomaly.

4. The gravity inversion method according to claim 1, characterized in that, The method of solving the gravity anomaly caused by the crust-mantle interface using forward modeling with spherical meshes is as follows: Based on the spherical harmonic model of the crust-mantle interface, the radial length set of discrete grid nodes generated by sampling from the reference sphere is calculated, and its minimum and maximum values ​​are statistically analyzed. Two reference spheres are then formed, with the crust-mantle interface embedded within each reference sphere. Spherical grid cells are formed from the three-layer interface and the latitude and longitude information of the grid nodes, and the latitude and longitude information of the center of each spherical grid cell is statistically analyzed. To obtain the corrected reference density of the shell-mantle interface, the known mantle density of the spherical harmonic model at the shell-mantle interface is subtracted from the shell corrected density. The density of spherical grid cells above the shell-mantle interface is set to a negative value, otherwise it is set to a positive value. The gravity anomaly at the shell-mantle interface is calculated by combining the location information of gravity observation points, the density of spherical grid cells, and the latitude and longitude information of the center.

5. The gravity inversion method according to claim 1, characterized in that, The method for solving gravity anomalies caused by topography using the finite amplitude method is as follows: Using the spherical surface corresponding to the planet's average radius as a reference surface, and combining planetary mass information and topographic spherical harmonic model, the data are input into the topographic gravity anomaly spherical harmonic coefficient calculation function of the open-source spherical harmonic analysis tool; the relative deviation of the power spectrum of the spherical harmonic coefficients of adjacent orders is calculated by increasing the finite amplitude order from the 2nd order; when the deviation is lower than the preset value, the corresponding topographic gravity spherical harmonic coefficient is obtained. By combining the location information of gravity observation points and the spherical harmonic coefficient of topographic gravity, gravity anomaly data caused by topography are obtained.

6. The gravity inversion method according to claim 1, characterized in that, The method for solving gravity anomalies caused by the crust-mantle interface using the finite amplitude method is as follows: The latitude and longitude range of the joint gravity observation points and the latitude and longitude coordinates of the center point of the study area were used to select a spherical cap, and the average radius within the spherical cap area was calculated by combining the crust-mantle interface model; Using the spherical surface corresponding to the average radius within the spherical cap region as the reference surface, the planetary mass, crust-mantle interface correction density, and calculation order of the spherical cap region are input into the gravity anomaly spherical harmonic coefficient calculation function of the open-source spherical harmonic tool; the relative deviation of the power spectrum of the spherical harmonic coefficients of adjacent orders is calculated by increasing the calculation order sequentially from order 2. When the deviation is lower than the preset value, the gravity spherical harmonic coefficient corresponding to the crust-mantle interface is obtained. By combining the gravitational spherical harmonic coefficient corresponding to the crust-mantle interface and the latitude and longitude information of the gravity observation point, the gravity anomaly data caused by the crust-mantle interface are calculated.

7. The gravity inversion method according to any one of claims 1 to 6, characterized in that, Step six specifically involves: inputting the location information of gravity observation points in the new observation coordinate system and the vertical gravity anomaly data into the three-dimensional gravity inversion model; constructing the target functional using the zero-order minimum entropy constraint combined with the regularized inversion method; solving the target functional using the Gauss-Newton optimization algorithm; obtaining the update direction of the three-dimensional gravity inversion model; determining the optimal update step size using the linear search method; updating the current three-dimensional gravity inversion model; repeating the iteration until the normalized fitting difference of the inversion result reaches the preset threshold or the number of iterations reaches the maximum value; and then outputting the inversion density result.

8. A gravity inversion system suitable for solid planetary crustal structures, characterized in that, include: The first observation coordinate system module performs terrain sampling on the sphere and obtains the discrete coordinates of the terrain sampling points by combining the terrain spherical harmonic model. The second observation coordinate system module determines the gravity observation reference surface, samples gravity observation points on the gravity observation surface, and obtains free air gravity anomaly data and discrete coordinates of gravity observation points by combining the gravity spherical harmonic model. The terrain correction module is used to solve the gravity anomaly caused by the terrain using the finite amplitude method or spherical mesh forward modeling, and to subtract the gravity anomaly caused by the terrain from the gravity anomaly caused by the terrain in the free air gravity anomaly. The crust-mantle interface correction module is used to solve the gravity anomaly caused by the crust-mantle interface by using the finite amplitude method or spherical mesh forward modeling. It subtracts the gravity anomaly caused by the crust-mantle from the free air gravity anomaly after removing the gravity anomaly caused by the terrain to obtain the radial gravity anomaly inside the crust. The observation projection correction module is used to determine a new observation coordinate system by rotating the discrete coordinates of the terrain sampling points and the discrete coordinates of the gravity observation points; and to obtain the vertical gravity anomaly data in the new observation coordinate system by angular projection of the radial gravity anomaly inside the shell. The spatial subdivision module is used to reconstruct the surface topography using surface coordinate points in the new observation coordinate system and obtain a three-dimensional gravity inversion model. The 3D inversion module is used to input the coordinates of gravity observation points in the new coordinate system and the vertical gravity anomaly data into the 3D gravity inversion model. By iteratively updating the 3D gravity inversion model, the density distribution inside the planetary shell is obtained.

9. The gravity inversion system according to claim 8, characterized in that, The observation projection correction module includes: The 3D rotation matrix construction unit is used to construct a 3D rotation matrix using the center point of the study area, the geocenter, and other arbitrary coordinate points. The three-dimensional coordinate rotation unit is used to multiply the discrete coordinates of terrain sampling points / gravity observation points with a three-dimensional rotation matrix to obtain the coordinates after three-dimensional rotation; The coordinate projection unit is used to translate the three-dimensional rotated coordinates along the z-axis, using the maximum z-value point as a reference, so that the coordinates are zeroed, and project the discrete coordinates of the three-dimensional rotated terrain sampling points / gravity observation points onto the horizontal plane to obtain a two-dimensional coordinate point set; A two-dimensional coordinate rotation unit is used to construct a two-dimensional rotation matrix, which is then multiplied by a set of two-dimensional coordinate points to obtain the two-dimensional coordinates after the plane has been rotated. The observation coordinate system reconstruction unit is used to combine the two-dimensional coordinates after plane rotation with the corresponding translated z-axis values ​​to obtain a new observation coordinate system after rotation and projection.

10. The gravity inversion system according to claim 8, characterized in that, The crust-mantle interface correction module includes a reference spherical construction unit, a second mesh generation unit, a crust-mantle density correction unit, and a crust-mantle interface gravity anomaly calculation unit. The reference spherical building element is used to calculate the radial length set of the shell-mantle interface spherical harmonic model at the discrete grid nodes, forming the shell-mantle interface; the minimum and maximum radial lengths of the shell-mantle interface are statistically analyzed to form two reference spheres, with the shell-mantle interface embedded in the two reference spheres. The second grid generation unit is used to form spherical grid units by taking the reference sphere and the crust-mantle interface according to the terrain sampling points, and to collect the latitude and longitude information of the center of each spherical grid unit. The shell-mantle density correction element is used to obtain the shell-mantle interface correction reference density by subtracting the known mantle density and the shell correction density from the spherical harmonic model of the shell-mantle interface. The density of spherical grid elements below the shell-mantle interface is set to a negative value, and vice versa. The crust-mantle gravity anomaly calculation unit is used to calculate the gravity anomaly corresponding to the crust-mantle interface by using the spherical grid cell density, center latitude and longitude information, and gravity observation point location information of the second grid generation unit; The shell-mantle interface correction module includes: an average radius calculation unit, a second spherical harmonic coefficient calculation unit, and a second spherical harmonic coefficient-to-gravity anomaly calculation unit; The average radius calculation unit is used to select a spherical cap by combining the latitude and longitude range of the gravity observation point and the latitude and longitude coordinates of the center point of the study area, and to calculate the average radius within the spherical cap area by combining the shell-mantle interface model; The second spherical harmonic coefficient calculation unit uses the spherical surface corresponding to the average radius within the spherical cap region as the reference surface, and inputs the planetary mass, crust-mantle interface correction density, and calculation order of the spherical cap region into the gravity anomaly spherical harmonic coefficient calculation function of the open-source spherical harmonic tool; according to the calculation order increasing sequentially from the 2nd order, the relative deviation of the power spectrum of the spherical harmonic coefficients of adjacent orders is calculated; when the deviation is lower than the preset value, the gravity spherical harmonic coefficient corresponding to its crust-mantle interface is obtained. The second spherical harmonic coefficient conversion gravity anomaly calculation unit combines the gravity spherical harmonic coefficient corresponding to the crust-mantle interface with the latitude and longitude information of the gravity observation point to obtain gravity anomaly data caused by the crust-mantle interface.

Images