CN116256808B - Forward and backward modeling method and system for regional gravity field - Google Patents

Forward and backward modeling method and system for regional gravity field Download PDF

Info

Publication number
CN116256808B
CN116256808B CN202310517692.6A CN202310517692A CN116256808B CN 116256808 B CN116256808 B CN 116256808B CN 202310517692 A CN202310517692 A CN 202310517692A CN 116256808 B CN116256808 B CN 116256808B
Authority
CN
China
Prior art keywords
data
gravity
gravity field
terrain
wave
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202310517692.6A
Other languages
Chinese (zh)
Other versions
CN116256808A (en
Inventor
杨萌
冯伟
钟敏
蒋弥
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Sun Yat Sen University
Original Assignee
Sun Yat Sen University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Sun Yat Sen University filed Critical Sun Yat Sen University
Priority to CN202310517692.6A priority Critical patent/CN116256808B/en
Publication of CN116256808A publication Critical patent/CN116256808A/en
Application granted granted Critical
Publication of CN116256808B publication Critical patent/CN116256808B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V7/00Measuring gravitational fields or waves; Gravimetric prospecting or detecting
    • G01V7/02Details
    • G01V7/06Analysis or interpretation of gravimetric records
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

The invention discloses a forward and backward fused regional gravity field modeling method and system, comprising the following steps: obtaining ground gravity anomaly data through ground gravity measurement data calculation; solving according to a global gravitational field model to obtain discrete point long wave gravitational field data; establishing a global layered terrain model, and obtaining high-frequency terrain data through spherical harmonic analysis and spherical harmonic synthesis; converting the high-frequency topographic data into discrete point short wave gravity field data through forward model theory; removing the discrete point long wave gravity field data and the discrete point short wave gravity field data from the ground gravity anomaly data to obtain residual gravity anomaly data; inversion is carried out on the basis of residual gravity anomaly data to obtain residual density anomaly data, and the residual gravity field caused by residual density anomaly and residual topography is solved by combining high-frequency topography data; and fusing the long-wave gravity field grid obtained based on the global gravity field model and the short-wave gravity field grid with the residual gravity field to obtain regional gravity field data.

Description

Forward and backward modeling method and system for regional gravity field
Technical Field
The invention relates to the technical field of gravitational and physical geodetic measurement, in particular to a method and a system for modeling a forward and backward fused regional gravitational field, which are suitable for regional gravitational field refinement and modeling of a gravity data sparse region.
Background
The gravity field reflects the mass distribution, movement and change in the earth, and has important practical significance and scientific significance for establishing global and regional mapping references, researching urgent topics such as sea level rising, geological evolution, resource investigation and the like. Fine regional gravitational field modeling has been one of the important scientific goals of modern geodetics, marine surveying, polar surveying, and solid geophysics. Currently, the construction of cm-level ground level and mGlul-level gravitational field is proposed internationally, and is a major challenge facing the field of current domestic and foreign ground measurement.
In recent years, with rapid developments in gravity measurement technology, satellite gravity detection technology, satellite altimetry technology, and the like, human knowledge of the global gravitational field has made a trans-epoch progress. The global gravity field model in the present stage, such as EGM2008, EIGEN-6C4, etc., can provide the gravity field signal with the spatial resolution half wavelength of 9 km, and the global ground level accuracy is due to meter level. Remote sensing, photogrammetry, etc. provide rich terrain data. Based on forward model theory, the topography and corresponding density information are known, and the fine medium short wave gravity field signal can be recovered according to the Newton gravitational force formula.
The current gravity field modeling in the gravity sparse area mainly integrates a global gravity field model and a terrain short-wave gravity field signal. Because of the lack of high-quality and high-resolution topographic density information, it is generally assumed that the topographic density is the constant 2670 kg/m3, and therefore, the medium-short wave gravity field signal constructed through the forward model does not contain density anomaly information. In areas with larger density values deviating from constant density hypothesis constants, such as volcanoes, mining areas and the like, larger errors can be introduced by forward model theory under the assumption that the density values are purely dependent on the constant density. In addition, the traditional residual topography model technology also has the problems of inaccurate reconciliation and correction, spectrum filtering and the like, and can achieve image modeling accuracy.
Disclosure of Invention
The invention provides a forward and backward fused regional gravity field modeling method and a forward and backward fused regional gravity field modeling system, which can effectively improve the fineness of gravity field modeling, and solve the problem of poor regional gravity field modeling precision caused by defects of constant density assumption, harmonic correction, spectral filtering and the like in a ground data sparse region.
In order to achieve the above purpose of the present invention, the following technical scheme is adopted:
a method for modeling a forward and backward fused regional gravity field comprises the following steps:
solving gravity anomaly by acquiring ground gravity measurement data in the area to obtain ground gravity anomaly data;
according to the selected global gravity field model, obtaining ground discrete point long wave gravity field data based on spherical harmonic synthesis solution;
the heterogeneous DEM data and SRTM15+ data are fused to establish a global layered terrain model, and terrain spherical surface filtering is completed through spherical harmonic analysis and spherical harmonic synthesis to obtain high-frequency terrain data;
converting high-frequency topographic data into ground discrete point short wave gravity field data by adopting a airspace method through forward model theory;
removing the discrete point long wave gravity field data and the discrete point short wave gravity field data from the ground gravity anomaly data to obtain residual gravity anomaly data;
inversion is carried out by adopting an equivalent source method based on the residual gravity anomaly data to obtain residual density anomaly data; solving a residual gravity field caused by residual density abnormality and residual topography by adopting a forward model theory in combination with residual density abnormality data and high-frequency topography data;
based on the selected global gravity field model, obtaining a long-wave gravity field grid through spherical harmonic synthesis and calculation; and obtaining a short wave gravity field grid through forward model theory calculation;
and fusing the long-wave gravity field grid, the short-wave gravity field grid and the residual gravity field to obtain regional gravity field data, and completing regional gravity field modeling.
Preferably, when the ground gravity measurement data and the global gravity field model do not belong to the same tidal system, the global gravity field model is uniformly transferred to the tidal system to which the ground gravity measurement data belong.
Preferably, the formula expression of the global gravity field model is as follows:
wherein,,representing the gravitational field of the area>Representing global long wave gravitational field data, +.>Representing and reconciling topographic shortwave gravity field data, +.>Gravitational field data generated for density anomalies.
Preferably, the heterogeneous DEM data comprises bedrock elevation information and ice elevation information; the SRTM15+ data comprise ocean sounding data and ground elevation data;
the heterogeneous DEM data adopts Eigen-6c4 geodetic plane as a reference plane, and the spatial resolution is 500 m; SRTM15+ data adopts an EGM96 ground level surface as a reference surface, and the spatial resolution is 450 m;
the fusion unifies the spatial resolution and the elevation reference of the heterogeneous DEM data and the SRTM15+ data, and finally builds a global layered terrain model with the spatial resolution of 500 m.
Preferably, spherical surface filtering of the terrain is completed through spherical harmonic analysis and spherical harmonic synthesis to obtain high-frequency terrain data, specifically, the global layered terrain model is subjected to spherical harmonic analysis to obtain a spherical harmonic unfolding form of the terrain, the spherical harmonic coefficient of the terrain with the same order as that of the global gravity field model is deducted, only a high-order item is reserved, spherical harmonic synthesis is carried out, and high-frequency terrain information in the area is obtained;
in the process of ball harmonic analysis and ball harmonic synthesis, an extended range algorithm is introduced, and the number is kept within the range of the double-precision numerical value by extending the floating point number index.
Preferably, the high-frequency topographic data is converted into discrete point short wave gravity field data by adopting a airspace method through forward model theory, and the method specifically comprises the following steps:
in the airspace, based on a layered terrain forward model technology, a plurality of forward model algorithms of polyhedron, prism, spherical crown and point quality are fused, global terrain quality is integrated, and a full-wave gravity field signal generated by a global layered terrain model is calculated:
wherein,,is a constant of universal gravitation->Is of topography density->For the spherical distance of the terrain quality element relative to the calculation point, < >>Representing the quality of the terrain>Representing the radial coordinates of the calculation points>For the earth radius>And->Elevation of upper and lower boundary of terrain, respectively, +.>For calculating the distance between the point and the integral point, < +.>For calculating the angular distance between the point and the integration point, < >>Andfor integrating radius>Full wave gravitational field signals generated for terrain;
in order to improve the calculation efficiency, the integration area is divided into different areas; the central area is 40 to km from the calculated point range, 3' resolution terrain is adopted, and each integral unit is replaced by prism approximation; the middle area distance calculation point 660 and km adopts 30' resolution terrain, and each integral unit is approximately replaced by a spherical crown body; the remote area can cover global topography, 15' resolution topography is adopted, and each integral unit is replaced by a point quality element;
in the spectral domain, resolving the influence of the long-wave gravity field and spectral filtering on the resolving of the high-frequency gravity field through three-time spherical harmonic analysis and synthesis; selecting a reference gravitational field model, the highest orderThe method comprises the steps of carrying out a first treatment on the surface of the Calculating the spherical harmonic coefficient of the layered terrain model with the same resolution as the reference gravitational field model through terrain spherical harmonic analysis; ball harmonic synthesis is carried out to obtain a smooth surface elevation grid;
the secondary spherical harmonic analysis of the smooth surface elevation grid is performed, and the ultra-high order is developed to obtain Gao Chengqiu harmonic coefficientsAnd then the high Cheng Qiuxie coefficient is converted into the spherical harmonic coefficient of the terrain gravity field:
wherein,,for gravitational field spherical harmonic coefficient, n and M are the order and number of expansion of spherical harmonic coefficient, M is earth mass, < >>Gao Chengmi times, & gt>Is->Is (are) factorial (I)>Is the density;
the long wave gravitational field is obtained by spherical harmonic synthesis:
and->The gravity anomaly and the gravity position of the terrain long wave are respectively;
discrete point short wave gravity field data is full wave gravity field signalAnd a long-wave gravitational field signal->The difference is:
further, the FFT technology is adopted to realize the conversion from the space elevation-density function to the gravitational field spherical harmonic coefficient, and the spherical harmonic synthesis is adopted to obtain the spherical harmonic composite material which is larger than the gravitational field spherical harmonic coefficientHigh frequency correction of (i.e. spectral filter correction->And finally obtaining the blended discrete point short wave gravity field data:
still further, removing long discrete point long wave gravity field data and discrete point short wave gravity field data from the ground gravity anomaly data to obtain a residual gravity signal:
wherein,,representing the gravitational disturbance measurement.
Preferably, based on the residual gravity anomaly data, inversion is carried out by adopting an equivalent source method to obtain residual density anomaly data; based on residual density anomaly data and high-frequency topographic data, solving a residual gravity field caused by residual density anomaly and residual topographic by adopting a forward model theory, wherein the method comprises the following steps of:
the basic idea of the equivalent source method is as follows: dividing the underground space into two-dimensional or three-dimensional regular grids, regarding each grid as a uniform cuboid, and inverting the density of each cuboid by the gravity value of an observation surface; the gravity formula of the three-dimensional cuboid is as follows:
in the method, in the process of the invention,is a universal gravitation constant; />Is the density; />For integrating voxel coordinates, +.>The distance between the coordinate point and the origin is set; />The lower and upper integration limits in the x-direction, respectively; />And->In the y-direction and in the z-direction;
the gravitational response of a three-dimensional cuboid grid at a series of observation points is written as the following linear equation:
in the method, in the process of the invention,constructing unknown vector->The remainder constitutes the design matrix->;/>Establishing an inversion observation equation for the observation vector;
the inversion objective function is constructed by adopting a direct inversion method, namely, an observation equation:
wherein,,is an observation data weight matrix,/->Is a model weight matrix,/->Is a regularization factor; the single-layer equivalent source directly takes the minimum value of the objective function without considering the model weight, namely the corresponding +.>As a regularized density anomaly solution:
adding two parts of smooth weight and depth weight into the model weight matrix, wherein the expression can be written as:
wherein,,is a global smooth factor, +.>、/>、/>Is a smoothing factor in three directions, +.>、/>、/>Is a corresponding smooth weight matrix; />For depth weight matrix, the depth weight function is used>The function of the formed diagonal matrix is written as:
wherein,,for observing the height of the surface +.>Is an attenuation factor->Depth;
the regularization solution of the observation equation corresponds to an objective function minimum, that is, the first derivative of the objective function is equal to 0, and regardless of the initial condition, the regularization solution can be written as:
based on layered equivalent sources, recovering residual gravity disturbance of any other point on the ground or in the air, and constructing a nuclear matrix according to the coordinates of the point to be solved and the grid coordinates of the equivalent sourcesCarry in the obtained equivalent source density +.>The residual gravity disturbance of the point to be solved can be calculated:
thereby realizing the recovery of the regional residual gravity field.
A system for modeling a forward and backward fused regional gravity field comprises
The gravity anomaly resolving module is used for resolving gravity anomalies by acquiring ground gravity measurement data in the area to obtain ground gravity anomaly data;
the discrete point long wave gravity field resolving module is used for resolving based on spherical harmonic synthesis according to the selected global gravity field model to obtain discrete point long wave gravity field data;
the high-frequency terrain data processing module is used for establishing a global layered terrain model by fusing heterogeneous DEM data and SRTM15+ data, and completing terrain spherical surface filtering through spherical harmonic analysis and spherical harmonic synthesis to obtain high-frequency terrain data;
the discrete point short wave gravity field calculation module is used for converting high-frequency topographic data into discrete point short wave gravity field data through forward model theory by adopting a airspace method;
the residual gravity anomaly data processing module is used for removing the discrete point long-wave gravity field data and the discrete point short-wave gravity field data from the ground gravity anomaly data to obtain residual gravity anomaly data;
the residual gravity field resolving module is used for inverting the residual density abnormal data by adopting an equivalent source method based on the residual gravity abnormal data; solving a residual gravity field caused by residual density abnormality and residual topography by adopting a forward model theory in combination with residual density abnormality data and high-frequency topography data;
the regional gravitational field modeling module is used for obtaining a long-wave gravitational field grid through spherical harmonic synthesis and calculation based on the selected global gravitational field model; and obtaining a short wave gravity field grid through forward model theory calculation; and fusing the long-wave gravity field grid, the short-wave gravity field grid and the residual gravity field to obtain regional gravity field data, and completing regional gravity field modeling.
The beneficial effects of the invention are as follows:
the invention combines two geophysical methods of a forward model and an inversion model for gravity field inversion of a gravity data insufficient area, and also overcomes the problem of non-harmony and spectral filtering in a classical residual topography model method.
According to the invention, by fusing the forward theory and the inversion theory, the dependence of the fine gravity field research on the fine gravity data is reduced; the problem of inaccurate classical harmonic correction is overcome by the harmonic nature of the spherical harmonic analysis and spherical harmonic synthesis. And the spectral filtering problem of classical spectrometry gravitational field filtering is also solved by introducing super Gao Jieqiu harmonic analysis. Through the process, the precision of gravity field modeling of the gravity data sparse region is improved, so that the method can be widely applied to the research of gravity fields of gravity sparse regions in developing countries, mountain regions, deserts and the like.
Drawings
Fig. 1 is a flowchart of a forward and backward fused regional gravity field modeling method according to the present invention.
FIG. 2 is a solution flow for converting high frequency topographic data into discrete point short wave gravity field data in accordance with the present invention.
Detailed Description
The invention is described in detail below with reference to the drawings and the detailed description.
Example 1
As shown in fig. 1, a method for modeling a forward and backward fused regional gravity field includes the following steps:
solving gravity anomaly by acquiring ground gravity measurement data in the area to obtain ground gravity anomaly data;
according to the selected global gravity field model, discrete point long wave gravity field data are obtained based on spherical harmonic synthesis solution;
the heterogeneous DEM data and SRTM15+ data are fused to establish a global layered terrain model, and terrain spherical surface filtering is completed through spherical harmonic analysis and spherical harmonic synthesis to obtain high-frequency terrain data;
converting high-frequency topographic data into discrete point short wave gravity field data by adopting a airspace method through a forward model theory;
removing the discrete point long wave gravity field data and the discrete point short wave gravity field data from the ground gravity anomaly data to obtain residual gravity anomaly data;
inversion is carried out by adopting an equivalent source method based on the residual gravity anomaly data to obtain residual density anomaly data; solving a residual gravity field caused by residual density abnormality and residual topography by adopting a forward model theory in combination with residual density abnormality data and high-frequency topography data;
based on the selected global gravity field model, obtaining a long-wave gravity field grid through spherical harmonic synthesis and calculation; and obtaining a short wave gravity field grid through forward model theory calculation;
and fusing the long-wave gravity field grid, the short-wave gravity field grid and the residual gravity field to obtain regional gravity field data, and completing regional gravity field modeling.
In a specific embodiment, when the ground gravity measurement data and the global gravity field model do not belong to the same tidal system, the global gravity field model is uniformly transferred to the tidal system to which the ground gravity measurement data belong.
The embodiment evaluates and selects the optimal global gravity field model according to a small amount of gravity data in the research area. For example, if the ground gravity field model is a tide-free tidal system and the global gravity field model is a zero-tide tidal system, and the ground gravity field model and the global gravity field model do not belong to the same tidal system, the global gravity field model is unified to the tide-free tidal system, and the conversion formula is as follows:
wherein,,and->The 2 nd order 0 th order spherical harmonic coefficients under the tide-free system and the zero-tide system are respectively.
It should be noted that: this embodiment provides for the interconversion between the tide-free, zero-tide and mean-tide tidal systems and is therefore uniform for tidal systems for all gravity data.
In a specific embodiment, in areas where gravity data is sparse, the gravitational field model is formulated as follows:
wherein,,representing the gravitational field of the area>Representing global long wave gravitational field data, +.>Representing and reconciling topographic shortwave gravity field data, +.>Gravitational field data generated for density anomalies. Topography shortwave gravitational field->Usually adopt the normal density +.>And (5) calculating.
In a specific embodiment, different terrain types are considered to have different densities, and a global layered terrain model is built to comprehensively consider the gravity influence of near zone and far zone terrain quality on the calculation points.
The global layered terrain model is established, and the method comprises the following steps: and fusing the heterogeneous DEM data with the SRTM15+ data to establish a global layered terrain model.
The heterogeneous DEM data comprise bedrock elevation information and ice elevation information; the SRTM15+ data comprise ocean sounding data and ground elevation data;
the heterogeneous DEM data adopts Eigen-6c4 geodetic plane as a reference plane, and the spatial resolution is 500 m; SRTM15+ data adopts an EGM96 ground level surface as a reference surface, and the spatial resolution is 450 m;
the fusion unifies the spatial resolution and the elevation reference of the heterogeneous DEM data and the SRTM15+ data, and finally builds a global layered terrain model with the spatial resolution of 500 m.
In a specific embodiment, the spherical surface filtering of the terrain is completed through spherical harmonic analysis and spherical harmonic synthesis to obtain high-frequency terrain data, specifically, the global layered terrain model is subjected to spherical harmonic analysis to obtain a spherical harmonic unfolding form of the terrain, the spherical harmonic coefficient of the terrain with the same order as that of the global gravity field model is subtracted, only high-order items are reserved, and the spherical harmonic synthesis is carried out to obtain high-frequency terrain information in the area;
in the ball harmonic analysis and synthesis process, the ultra-high order association Legend function calculation is involved, an extended range algorithm (Extended Range Arithmetic) is introduced for ensuring the stability of calculation, and the number is kept within the common double-precision numerical range by extending the floating point number index, so that the overflow problem is eliminated.
In a specific embodiment, as shown in fig. 2, the high-frequency topographic data is converted into discrete point short wave gravity field data by adopting a space domain method through forward model theory, which is specifically as follows:
in the airspace, based on a layered terrain forward model technology, a plurality of forward model algorithms of polyhedron, prism, spherical crown and point quality are fused, global terrain quality is integrated, and a full-wave gravity field signal generated by a global layered terrain model is calculated:
wherein,,is a constant of universal gravitation->Is of topography density->For the spherical distance of the terrain quality element relative to the calculation point, < >>Representing the quality of the terrain>Representing the radial coordinates of the calculation points>For the earth radius>And->Elevation of upper and lower boundary of terrain, respectively, +.>For calculating the distance between the point and the integral point, < +.>For calculating the angular distance between the point and the integration point, < >>Andfor integrating radius>Full wave gravitational field signals generated for terrain;
in order to improve the calculation efficiency, the integration area is divided into different areas; the central area is 40 to km from the calculated point range, 3' resolution terrain is adopted, and each integral unit is replaced by prism approximation; the middle area distance calculation point 660 and km adopts 30' resolution terrain, and each integral unit is approximately replaced by a spherical crown body; the remote area can cover global topography, 15' resolution topography is adopted, and each integral unit is replaced by a point quality element;
in the spectral domain, resolving the influence of the long-wave gravity field and spectral filtering on the resolving of the high-frequency gravity field through three-time spherical harmonic analysis and synthesis; selecting a reference gravitational field model, the highest orderThe method comprises the steps of carrying out a first treatment on the surface of the Calculating the spherical harmonic coefficient of the layered terrain model with the same resolution as the reference gravitational field model through terrain spherical harmonic analysis; ball harmonic synthesis is carried out to obtain a smooth surface elevation grid; the reference gravitational field model can be ITSG, GO_CONS_GCF or XGM series. Because the terrain filtering is inconsistent with the gravitational field filtering, the smooth surface can generate ultrahigh frequency gravitational field information, and the part of high frequency information is ignored in the airspace method resolving process, so that correction is needed.
The secondary spherical harmonic analysis of the smooth surface elevation grid is performed, and the ultra-high order is developed to obtain Gao Chengqiu harmonic coefficientsAnd then the high Cheng Qiuxie coefficient is converted into the spherical harmonic coefficient of the terrain gravity field:
wherein,,for gravitational field spherical harmonic coefficient, n and M are the order and number of expansion of spherical harmonic coefficient, M is earth mass, < >>Gao Chengmi times, & gt>Is->Is (are) factorial (I)>Is the density;
the long wave gravitational field is obtained by spherical harmonic synthesis:
and->The gravity anomaly and the gravity position of the terrain long wave are respectively;
discrete point short wave gravity field data is full wave gravity field signalAnd a long-wave gravitational field signal->The difference is:
reconciliation correction problem: the problem of 'non-reconciliation' is a main problem of resolving a high-frequency gravitational field based on a residual terrain model technology, and a classical compression method is generally adopted for reconciling and correcting, but the method is influenced by infinite Bragg assumption, and larger errors are introduced in a region with larger terrain fluctuation. In the embodiment, the long-wave gravitational field is resolved by adopting a mode of spherical harmonic analysis and spherical harmonic synthesis, so that the gravitational field at any resolved point meets the reconciliation condition, and the reconciliation and correction problem of the high-frequency gravitational field resolved by the traditional airspace method is solved.
In a specific embodiment, the topographic gravity field filtering is usually developed based on spectroscopy, and the problem of inconsistent topographic filtering and gravity field filtering is ignored, so that a large error is introduced. The embodiment adopts FFT technology to realize the conversion from space elevation-density function to gravitational field spherical harmonic coefficient, and obtains the spherical harmonic compound to be larger than thatHigh frequency correction of (i.e. spectral filter correction->And finally obtaining the blended discrete point short wave gravity field data:
in a specific embodiment, long discrete point long wave gravity field data and short discrete point short wave gravity field data are removed from the ground gravity anomaly data, and a residual gravity anomaly signal is obtained:
wherein,,representing the gravitational disturbance measurement.
In a specific embodiment, inversion is performed by an equivalent source method based on residual gravity anomaly data to obtain residual density anomaly data; based on residual density anomaly data and high-frequency topographic data, solving a residual gravity field caused by residual density anomaly and residual topographic by adopting a forward model theory, wherein the method comprises the following steps of:
the basic idea of the equivalent source method is as follows: dividing an underground space into two-dimensional or three-dimensional regular grids, regarding each grid as a uniform cuboid, and inverting the density of each cuboid by the gravity value of an observation surface (ground); the gravity formula of the three-dimensional cuboid is as follows:
in the method, in the process of the invention,is a universal gravitation constant; />Is the density; />For integrating voxel coordinates, +.>The distance between the coordinate point and the origin is set; />The lower and upper integration limits in the x-direction, respectively; />And->In the y-direction and in the z-direction;
the gravitational response of a three-dimensional cuboid grid at a series of observation points is written as the following linear equation:
in the method, in the process of the invention,constructing unknown vector->The remainder constitutes the design matrix->;/>Establishing an inversion observation equation for the observation vector;
the inversion objective function is constructed by adopting a direct inversion method, namely, an observation equation:
wherein,,is an observation data weight matrix,/->Is a model weight matrix,/->Is a regularization factor; the single-layer equivalent source directly takes the minimum value of the objective function without considering the model weight, namely the corresponding +.>As a regularized density anomaly solution:
adding two parts of smooth weight and depth weight into the model weight matrix, wherein the expression can be written as:
wherein,,is a global smooth factor, +.>、/>、/>Is a smoothing factor in three directions, +.>、/>、/>Is a corresponding smooth weight matrix; />For depth weight matrix, the depth weight function is used>The function of the formed diagonal matrix is written as:
wherein,,for observing the height of the surface +.>Is an attenuation factor->Depth;
the regularization solution of the observation equation corresponds to an objective function minimum, that is, the first derivative of the objective function is equal to 0, and regardless of the initial condition, the regularization solution can be written as:
based on layered equivalent sources, recovering residual gravity disturbance of any other point on the ground or in the air, and constructing a nuclear matrix according to the coordinates of the point to be solved and the grid coordinates of the equivalent sourcesCarry in the obtained equivalent source density +.>The residual gravity disturbance of the point to be solved can be calculated:
thereby realizing the recovery of the regional residual gravity field.
In the embodiment, unconstrained inversion is performed by using a two-dimensional or three-dimensional gravity inversion technology, so that an underground shallow equivalent density layer or a three-dimensional equivalent density structure of a research area is obtained, and the underground shallow equivalent density layer or the three-dimensional equivalent density structure is used as an equivalent source of an area residual gravity field.
Based on the embodiment, the long-wave gravity field grid, the short-wave gravity field grid and the residual gravity field are fused to obtain regional gravity field data, and regional gravity field modeling is completed. The present embodiment also verifies the accuracy of the resolved regional gravity field based on the independent gravity data and GPS level data.
In the embodiment, based on the density anomaly inversion result recovered from a small amount of residual gravity data in the test area, the density anomaly result of each 3km depth is inverted, and the density anomaly difference in the volcanic area and the volcanic bag area is larger. Further, according to the residual gravity disturbance result calculated by the density anomaly, the residual gravity disturbance is changed within the range of-10 to 10 mGlul, and the experimental result of the checking point shows that the method improves the modeling accuracy of the volcanic area gravity field by about 30%.
Example 2
The embodiment also provides a system for forward and backward fused regional gravity field modeling method based on the forward and backward fused regional gravity field modeling method of the embodiment 1, which comprises
The gravity anomaly resolving module is used for resolving gravity anomalies by acquiring ground gravity measurement data in the area to obtain ground gravity anomaly data;
the first long-wave gravitational field resolving module is used for resolving based on spherical harmonic synthesis according to the selected global gravitational field model to obtain discrete point long-wave gravitational field data;
the high-frequency terrain data processing module is used for establishing a global layered terrain model by fusing heterogeneous DEM data and SRTM15+ data, and completing terrain spherical surface filtering through spherical harmonic analysis and spherical harmonic synthesis to obtain high-frequency terrain data;
the discrete point short wave gravity field calculation module is used for converting high-frequency topographic data into discrete point short wave gravity field data through forward model theory by adopting a airspace method;
the residual gravity anomaly data processing module is used for removing the discrete point long-wave gravity field data and the discrete point short-wave gravity field data from the ground gravity anomaly data to obtain residual gravity anomaly data;
the residual gravity field resolving module is used for inverting the residual density abnormal data by adopting an equivalent source method based on the residual gravity abnormal data; solving a residual gravity field caused by residual density abnormality and residual topography by adopting a forward model theory in combination with residual density abnormality data and high-frequency topography data;
the regional gravitational field modeling module is used for obtaining a long-wave gravitational field grid through spherical harmonic synthesis and calculation based on the selected global gravitational field model; and obtaining a short wave gravity field grid through forward model theory calculation; and fusing the long-wave gravity field grid, the short-wave gravity field grid and the residual gravity field to obtain regional gravity field data, and completing regional gravity field modeling.
Example 3
A computer device comprising a memory, a processor and a computer program stored on the memory and executable on the processor, said processor implementing the steps of the method according to embodiment 1 when said computer program is executed.
Where the memory and the processor are connected by a bus, the bus may comprise any number of interconnected buses and bridges, the buses connecting the various circuits of the one or more processors and the memory together. The bus may also connect various other circuits such as peripherals, voltage regulators, and power management circuits, which are well known in the art, and therefore, will not be described any further herein. The bus interface provides an interface between the bus and the transceiver. The transceiver may be one element or may be a plurality of elements, such as a plurality of receivers and transmitters, providing a means for communicating with various other apparatus over a transmission medium. The data processed by the processor is transmitted over the wireless medium via the antenna, which further receives the data and transmits the data to the processor.
Still further embodiments provide a computer readable storage medium having stored thereon a computer program which, when executed by a processor, implements the steps of the method as described in embodiment 1.
That is, it will be understood by those skilled in the art that all or part of the steps in implementing the methods of the embodiments described above may be implemented by a program stored in a storage medium, where the program includes several instructions for causing a device (which may be a single-chip microcomputer, a chip or the like) or a processor (processor) to perform all or part of the steps in the methods of the embodiments described herein. And the aforementioned storage medium includes: a U-disk, a removable hard disk, a Read-only memory (ROM), a random access memory (RAM, random Access Memory), a magnetic disk, or an optical disk, or other various media capable of storing program codes.
It is to be understood that the above examples of the present invention are provided by way of illustration only and not by way of limitation of the embodiments of the present invention. Any modification, equivalent replacement, improvement, etc. which come within the spirit and principles of the invention are desired to be protected by the following claims.

Claims (9)

1. A forward and backward modeling method of a fused regional gravity field is characterized by comprising the following steps: the method comprises the following steps:
solving gravity anomaly by acquiring ground gravity measurement data in the area to obtain ground gravity anomaly data;
according to the selected global gravity field model, obtaining ground discrete point long wave gravity field data based on spherical harmonic synthesis solution;
the heterogeneous DEM data and SRTM15+ data are fused to establish a global layered terrain model, and terrain spherical surface filtering is completed through spherical harmonic analysis and spherical harmonic synthesis to obtain high-frequency terrain data;
converting high-frequency topographic data into ground discrete point short wave gravity field data by adopting a airspace method through forward model theory;
removing the ground discrete point long-wave gravity field data and the ground discrete point short-wave gravity field data from the ground gravity anomaly data to obtain residual gravity anomaly data;
inversion is carried out by adopting an equivalent source method based on the residual gravity anomaly data to obtain residual density anomaly data; solving a residual error gravity field caused by residual error density abnormality and high-frequency topography by adopting a forward model theory in combination with residual error density abnormality data and high-frequency topography data;
based on the selected global gravity field model, obtaining a long-wave gravity field grid through spherical harmonic synthesis and calculation; and obtaining a short wave gravity field grid through forward model theory calculation;
fusing the long-wave gravity field grid, the short-wave gravity field grid and the residual gravity field to obtain regional gravity field data, and completing regional gravity field modeling;
the high-frequency topographic data is converted into ground discrete point short wave gravity field data by adopting a airspace method through forward model theory, and the method is concretely as follows:
in the airspace, based on a layered terrain forward model technology, a plurality of forward model algorithms of polyhedron, prism, spherical crown and point quality are fused, global terrain quality is integrated, and a full-wave gravity field signal generated by a global layered terrain model is calculated:
wherein,,is a constant of universal gravitation->Is of topography density->For the spherical distance of the terrain quality element relative to the calculation point, < >>Representing the quality of the terrain>Representing the radial coordinates of the calculation points>For the earth radius>And->Elevation of the lower and upper boundary of the terrain, respectively,/->For calculating the distance between the point and the integral point, < +.>For calculating the angular distance between the point and the integration point, < >>And->For integrating radius>Full wave gravitational field signals generated for terrain;
in order to improve the calculation efficiency, the integration area is divided into different areas; the central area is 40 to km from the calculated point range, 3' resolution terrain is adopted, and each integral unit is replaced by prism approximation; the middle area distance calculation point 660 and km adopts 30' resolution terrain, and each integral unit is approximately replaced by a spherical crown body; the remote area can cover global topography, 15' resolution topography is adopted, and each integral unit is replaced by a point quality element;
in the spectral domain, resolving the influence of the long-wave gravity field and spectral filtering on the resolving of the high-frequency gravity field through three-time spherical harmonic analysis and synthesis; selecting a reference gravitational field model, the highest orderThe method comprises the steps of carrying out a first treatment on the surface of the Calculating the spherical harmonic coefficient of the layered terrain model with the same resolution as the reference gravitational field model through terrain spherical harmonic analysis; ball harmonic synthesis is carried out to obtain a smooth surface elevation grid;
the secondary spherical harmonic analysis of the smooth surface elevation grid is performed, and the ultra-high order is developed to obtain Gao Chengqiu harmonic coefficientsAnd then the high Cheng Qiuxie coefficient is converted into the spherical harmonic coefficient of the terrain gravity field:
wherein,,is the spherical harmonic coefficient of the gravitational field, n andm is the order and the number of times the spherical harmonic coefficient is developed, M is the earth mass, ++>Gao Chengmi times, & gt>Is->Is (are) factorial (I)>Is the density;
the long wave gravitational field is obtained by spherical harmonic synthesis:
and->The gravity anomaly and the gravity position of the terrain long wave are respectively;
the ground discrete point short wave gravity field data is full wave gravity field signalAnd a long-wave gravitational field signal->The difference is:
2. the method for modeling the forward and backward fused regional gravity field according to claim 1, which is characterized by comprising the following steps of: when the ground gravity measurement data and the global gravity field model do not belong to the same tidal system, the global gravity field model is uniformly transferred to the tidal system to which the ground gravity measurement data belong.
3. The method for modeling the forward and backward fused regional gravity field according to claim 1, which is characterized by comprising the following steps of: the formula expression of the global gravity field model is as follows:
wherein,,representing the gravitational field of the area>Representing global long wave gravitational field data, +.>Representing harmonious topographic shortwave gravity field data +.>Gravitational field data generated for density anomalies.
4. The method for modeling the forward and backward fused regional gravity field according to claim 1, which is characterized by comprising the following steps of: the heterogeneous DEM data comprise bedrock elevation information and ice elevation information; the SRTM15+ data comprise ocean sounding data and ground elevation data;
the heterogeneous DEM data adopts Eigen-6c4 geodetic plane as a reference plane, and the spatial resolution is 500 m; SRTM15+ data adopts an EGM96 ground level surface as a reference surface, and the spatial resolution is 450 m;
the fusion unifies the spatial resolution and the elevation reference of the heterogeneous DEM data and the SRTM15+ data, and finally builds a global layered terrain model with the spatial resolution of 500 m.
5. The method for modeling the forward and backward fused regional gravity field according to claim 1, which is characterized by comprising the following steps of: the method comprises the steps of finishing terrain spherical filtering through spherical harmonic analysis and spherical harmonic synthesis to obtain high-frequency terrain data, specifically, obtaining a spherical harmonic unfolding form of the terrain through spherical harmonic analysis by the global layered terrain model, deducting the spherical harmonic coefficient of the terrain with the same order as the global gravity field model, only retaining a high-order item, and obtaining high-frequency terrain information in a region through spherical harmonic synthesis;
in the process of ball harmonic analysis and ball harmonic synthesis, an extended range algorithm is introduced, and the number is kept within the range of the double-precision numerical value by extending the floating point number index.
6. The method for modeling the forward and backward fused regional gravity field according to claim 1, which is characterized by comprising the following steps of: the FFT technology is adopted to realize the conversion from the space elevation-density function to the gravitational field spherical harmonic coefficient, and the spherical harmonic synthesis is adopted to obtain the spherical harmonicHigh frequency correction of (i.e. spectral filter correction->And finally obtaining the blended ground discrete point short wave gravity field data:
7. the method for modeling the forward and backward fused regional gravity field according to claim 1, which is characterized by comprising the following steps of: removing ground discrete point long-wave gravity field data and ground discrete point short-wave gravity field data from ground gravity anomaly data to obtain residual gravity anomaly data, wherein the specific formula is as follows:
wherein,,representing the ground gravity anomaly.
8. The method for modeling the forward and backward fused regional gravity field according to claim 1, which is characterized by comprising the following steps of: inversion is carried out by adopting an equivalent source method based on the residual gravity anomaly data to obtain residual density anomaly data; combining residual density anomaly data and high-frequency terrain data, solving a residual gravity field caused by residual density anomaly and high-frequency terrain by adopting a forward model theory, wherein the method comprises the following steps of:
the basic idea of the equivalent source method is as follows: dividing the underground space into three-dimensional regular grids, regarding each grid as a uniform cuboid, and inverting the density of each cuboid by the gravity value of an observation surface; the gravity formula of the three-dimensional cuboid is as follows:
in the method, in the process of the invention,is a universal gravitation constant; />Is the density; />For integrating voxel coordinates, +.>The distance between the coordinate point and the origin is set; />The lower and upper integration limits in the x-direction, respectively; />Is the lower and upper integral limit in the y-direction,/->Is the lower and upper integral limits in the z-direction;
the gravitational response of a three-dimensional cuboid grid at a series of observation points is written as the following linear equation:
in the method, in the process of the invention,constructing unknown vector->The remainder constitutes the design matrix->;/>Establishing an inversion observation equation for the observation vector;
the inversion objective function is constructed by adopting a direct inversion method, namely, an observation equation:
wherein,,is an observation data weight matrix,/->Is a model weight matrix,/->Is a regularization factor; the single-layer equivalent source directly takes the minimum value of the objective function without considering the model weight, namely the corresponding +.>As a regularized density anomaly solution:
adding two parts of smooth weight and depth weight into the model weight matrix, wherein the expression can be written as:
wherein,,is a global smooth factor, +.>、/>、/>Is a smoothing factor in three directions, +.>、/>、/>Is a corresponding smooth weight matrix; />For depth weight matrix, the depth weight function is used>The function of the formed diagonal matrix is written as:
wherein,,for observing the height of the surface +.>Is an attenuation factor->Depth;
the regularization solution of the observation equation corresponds to an objective function minimum, that is, the first derivative of the objective function is equal to 0, and regardless of the initial condition, the regularization solution can be written as:
based on layered equivalent sources, recovering residual gravity disturbance of any other point on the ground or in the air, and constructing a nuclear matrix according to the coordinates of the point to be solved and the grid coordinates of the equivalent sourcesSubstituting the obtained equivalent source density +.>The residual gravity disturbance of the point to be solved can be calculated:
thereby realizing the recovery of the regional residual gravity field.
9. A forward and backward fused regional gravity field modeling system is characterized in that: comprising the following steps:
the gravity anomaly resolving module is used for resolving gravity anomalies by acquiring ground gravity measurement data in the area to obtain ground gravity anomaly data;
the discrete point long wave gravity field resolving module is used for resolving based on spherical harmonic synthesis according to the selected global gravity field model to obtain ground discrete point long wave gravity field data;
the high-frequency terrain data processing module is used for establishing a global layered terrain model by fusing heterogeneous DEM data and SRTM15+ data, and completing terrain spherical surface filtering through spherical harmonic analysis and spherical harmonic synthesis to obtain high-frequency terrain data;
the discrete point short wave gravity field resolving module is used for converting high-frequency topographic data into ground discrete point short wave gravity field data through forward model theory by adopting a airspace method;
the residual gravity anomaly data processing module is used for removing the ground discrete point long-wave gravity field data and the ground discrete point short-wave gravity field data from the ground gravity anomaly data to obtain residual gravity anomaly data;
the residual gravity field resolving module is used for inverting the residual density abnormal data by adopting an equivalent source method based on the residual gravity abnormal data; solving a residual error gravity field caused by residual error density abnormality and high-frequency topography by adopting a forward model theory in combination with residual error density abnormality data and high-frequency topography data;
the regional gravitational field modeling module is used for obtaining a long-wave gravitational field grid through spherical harmonic synthesis and calculation based on the selected global gravitational field model; and obtaining a short wave gravity field grid through forward model theory calculation; fusing the long-wave gravity field grid, the short-wave gravity field grid and the residual gravity field to obtain regional gravity field data, and completing regional gravity field modeling;
the high-frequency topographic data is converted into ground discrete point short wave gravity field data by adopting a airspace method through forward model theory, and the method is concretely as follows:
in the airspace, based on a layered terrain forward model technology, a plurality of forward model algorithms of polyhedron, prism, spherical crown and point quality are fused, global terrain quality is integrated, and a full-wave gravity field signal generated by a global layered terrain model is calculated:
wherein,,is a constant of universal gravitation->Is of topography density->For the spherical distance of the terrain quality element relative to the calculation point, < >>Representing the quality of the terrain>Representing the radial coordinates of the calculation points>For the earth radius>And->Elevation of the lower and upper boundary of the terrain, respectively,/->For calculating the distance between the point and the integral point, < +.>For calculating the angular distance between the point and the integration point, < >>And->For integrating radius>Full wave gravitational field signals generated for terrain;
in order to improve the calculation efficiency, the integration area is divided into different areas; the central area is 40 to km from the calculated point range, 3' resolution terrain is adopted, and each integral unit is replaced by prism approximation; the middle area distance calculation point 660 and km adopts 30' resolution terrain, and each integral unit is approximately replaced by a spherical crown body; the remote area can cover global topography, 15' resolution topography is adopted, and each integral unit is replaced by a point quality element;
in the spectral domain, resolving the influence of the long-wave gravity field and spectral filtering on the resolving of the high-frequency gravity field through three-time spherical harmonic analysis and synthesis; selecting a reference gravitational field model, the highest orderThe method comprises the steps of carrying out a first treatment on the surface of the Calculating the spherical harmonic coefficient of the layered terrain model with the same resolution as the reference gravitational field model through terrain spherical harmonic analysis; ball harmonic synthesis is carried out to obtain a smooth surface elevation grid;
the secondary spherical harmonic analysis of the smooth surface elevation grid is performed, and the ultra-high order is developed to obtain Gao Chengqiu harmonic coefficientsAnd then the high Cheng Qiuxie coefficient is converted into the spherical harmonic coefficient of the terrain gravity field:
wherein,,for gravitational field spherical harmonic coefficient, n and M are the order and number of expansion of spherical harmonic coefficient, M is earth mass, < >>Gao Chengmi times, & gt>Is->Is (are) factorial (I)>Is the density;
the long wave gravitational field is obtained by spherical harmonic synthesis:
and->The gravity anomaly and the gravity position of the terrain long wave are respectively;
the ground discrete point short wave gravity field data is full wave gravity field signalAnd a long-wave gravitational field signal->The difference is:
CN202310517692.6A 2023-05-10 2023-05-10 Forward and backward modeling method and system for regional gravity field Active CN116256808B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310517692.6A CN116256808B (en) 2023-05-10 2023-05-10 Forward and backward modeling method and system for regional gravity field

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310517692.6A CN116256808B (en) 2023-05-10 2023-05-10 Forward and backward modeling method and system for regional gravity field

Publications (2)

Publication Number Publication Date
CN116256808A CN116256808A (en) 2023-06-13
CN116256808B true CN116256808B (en) 2023-08-01

Family

ID=86682804

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310517692.6A Active CN116256808B (en) 2023-05-10 2023-05-10 Forward and backward modeling method and system for regional gravity field

Country Status (1)

Country Link
CN (1) CN116256808B (en)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117147552A (en) * 2023-10-30 2023-12-01 北京交通大学 Rock slag grading analysis method for TBM tunnel

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102169178B (en) * 2010-12-21 2012-11-21 中国测绘科学研究院 Method for determining sea surface topographic structure based on barodynamics
CN108267792B (en) * 2018-04-13 2019-07-12 武汉大学 Building global gravitational field model inversion method
CN113505479B (en) * 2021-07-05 2023-05-12 中国科学院地质与地球物理研究所 Density inversion method and device and electronic equipment
CN113779797B (en) * 2021-09-13 2024-06-18 中国地质调查局自然资源综合调查指挥中心 Three-dimensional gravity forward modeling method and system based on residual density model

Also Published As

Publication number Publication date
CN116256808A (en) 2023-06-13

Similar Documents

Publication Publication Date Title
Rexer et al. Layer-based modelling of the Earth’s gravitational potential up to 10-km scale in spherical harmonics in spherical and ellipsoidal approximation
Denker Regional gravity field modeling: theory and practical results
Balmino et al. Spherical harmonic modelling to ultra-high degree of Bouguer and isostatic anomalies
CN116256808B (en) Forward and backward modeling method and system for regional gravity field
Kiamehr et al. Effect of the SRTM global DEM on the determination of a high-resolution geoid model: a case study in Iran
Zhu et al. SDUST2021GRA: Global marine gravity anomaly model recovered from Ka-band and Ku-band satellite altimeter data
Abdalla et al. A new gravimetric geoid model for Sudan using the KTH method
Shen et al. Improved geoid determination based on the shallow-layer method: a case study using EGM08 and CRUST2. 0 in the Xinjiang and Tibetan regions
Li et al. The contribution of the GRAV‐D airborne gravity to geoid determination in the Great Lakes region
CN111257956A (en) Matlab-based regional quasi-geoid surface refinement method
Bajracharya Terrain effects on geoid determination
Wu et al. High-resolution regional gravity field recovery from Poisson wavelets using heterogeneous observational techniques
Chen et al. A high speed method of SMTS
Sulzbach et al. Modeling gravimetric signatures of third-degree ocean tides and their detection in superconducting gravimeter records
Ji et al. On performance of CryoSat-2 altimeter data in deriving marine gravity over the Bay of Bengal
Feng et al. A hierarchical network densification approach for reconstruction of historical ice velocity fields in East Antarctica
Root et al. Benchmark forward gravity schemes: the gravity field of a realistic lithosphere model WINTERC-G
Lestari et al. Local geoid modeling in the central part of Java, Indonesia, using terrestrial-based gravity observations
CN110287620B (en) Spherical coordinate system density interface forward modeling method and system suitable for earth surface observation surface
Kim et al. Utility of Slepian basis functions for modeling near-surface and satellite magnetic anomalies of the Australian lithosphere
CN113885101B (en) Method for constructing gravity gradient reference map based on ellipsoidal model
Fairhead et al. The use of GPS in gravity surveys
Li et al. The SDUST2022GRA global marine gravity anomalies recovered from radar and laser altimeter data: Contribution of ICESat-2 laser altimetry
Catalão et al. Mapping the geoid for Iberia and the Macaronesian Islands using multi-sensor gravity data and the GRACE geopotential model
Koch Simple layer model of the geopotential in satellite geodesy

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant