CN112363236B - Gravity field data equivalent source continuation and data type conversion method based on PDE - Google Patents

Gravity field data equivalent source continuation and data type conversion method based on PDE Download PDF

Info

Publication number
CN112363236B
CN112363236B CN202011103444.XA CN202011103444A CN112363236B CN 112363236 B CN112363236 B CN 112363236B CN 202011103444 A CN202011103444 A CN 202011103444A CN 112363236 B CN112363236 B CN 112363236B
Authority
CN
China
Prior art keywords
gravity
data
equivalent source
pde
representing
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
CN202011103444.XA
Other languages
Chinese (zh)
Other versions
CN112363236A (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.)
China University of Geosciences
Original Assignee
China University of Geosciences
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 China University of Geosciences filed Critical China University of Geosciences
Priority to CN202011103444.XA priority Critical patent/CN112363236B/en
Publication of CN112363236A publication Critical patent/CN112363236A/en
Application granted granted Critical
Publication of CN112363236B publication Critical patent/CN112363236B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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

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 provides a gravity field data equivalent source continuation and data type conversion method based on PDE, which comprises the following steps: inputting gravity data and establishing a topographic relief curved surface; determining a space range of mesh generation according to the elevation information of the fluctuating observation curved surface and the preset inversion maximum depth, performing continuous structural non-uniform mesh generation, and determining an equivalent source inversion mesh space; according to the gravity background field, carrying out PDE three-dimensional inversion calculation with depth normalization factors, positive value constraint terms and normalization terms on gravity field data based on equivalent source inversion grid space to obtain a multilayer equivalent source model; different observation parameters are given, and the model is utilized to carry out three-dimensional forward modeling calculation on the gravity field of the nonlinear PDE so as to obtain one or more of gravity abnormal plumb data, gravity abnormal field data and gravity gradient tensor data. The invention has the beneficial effects that: the method can perform self-adaptive and high-precision continuation and data type conversion calculation on the underground gravity anomaly data in the complex environment.

Description

Gravity field data equivalent source continuation and data type conversion method based on PDE
Technical Field
The invention relates to the technical field of geophysical surveying, in particular to a gravity field data equivalent source continuation and data type conversion method based on PDE.
Background
In gravity detection, people usually detect the vertical first derivative of a gravity field, but in actual data interpretation, the data often needs to be converted into required parameters and types, such as conversion of different observation heights (continuation) and gradient tensor data, and the like.
The prior art mainly comprises the following steps: (1) the equivalent source method utilizes a virtual field source to simulate the actual measurement abnormity and can be used for space continuation (including curved surface continuation), gradient calculation, component conversion and the like of the bit field data; selecting a single layer of equivalent sources and placing them near the surface is a key feature of the equivalent source method.
(2) A single-layer equivalent source is used for analyzing gravity field data, and the single-layer equivalent source is used for denoising gravity gradient tensor data; in order to reconstruct the potential field with high precision while ensuring the calculation efficiency, the multilayer equivalent source method is a relatively reasonable choice,
(3) the document "plum ends; aging; fipronil; carrying out roof beam green treatment; a velvet sea dragon; a gravity data processing method [ C ]//2017 annual meeting of China Union of geoscience, 0. "based on multilayer equivalent sources provides a method for dividing underground equivalent sources into a shallow layer, a middle layer and a deep layer to realize the conversion of gravity field data,
(4) the document "Chentao, Zhang Guibin, equivalent source technology [ J ] for calculating gravity gradient by using gravity anomaly.
The application of these technologies enables people to perform polarization and data type conversion operation of ground gravity data through equivalent source technology, but there are the following problems: 1) the construction of equivalent source grids is only three layers at most, and the grid division is discontinuous, in the prior art, the number of the set layers of equivalent sources is limited, and the grid division of an underground space is unreasonable, so that the precision and the speed of gravity field continuation and data type conversion are influenced; 2) the depth position of each equivalent layer needs to be estimated independently and then placed independently; 3) the equivalent source is calculated by adopting a linear forward modeling method and a linear model inversion method, and the anti-interference capability and the calculation accuracy are poor.
Disclosure of Invention
In view of the above, the technical problems to be solved by the present invention are: how to perform equivalent source extension and data type conversion on the gravity data of the irregular measurement position d and improve the calculation accuracy;
therefore, the invention provides a gravity data equivalent source continuation and data type conversion method based on PDE, which adopts continuous structural non-uniform mesh subdivision, introduces a depth normalization factor based on a forward and backward theoretical framework of nonlinear Differential PDE (Partial Differential Equations, and can also comprise finite elements, finite volumes and finite difference methods), directly determines the depth and distribution of an equivalent source in the inversion process, and can calculate the gravity field data of the whole three-dimensional free space at one time.
The invention provides a gravity data equivalent source continuation and data type conversion method based on PDE, which comprises the following steps:
s1, acquiring gravity field data d on the undulation observation curved surface0Establishing a topographic relief curved surface according to the topographic height information of the area where the gravity field data is located; the undulation observation curved surface is a recorded known value during observation;
s2, determining a space range of mesh subdivision according to the elevation information of the undulation observation surface and the preset inversion maximum depth, and performing continuous structural non-uniform mesh subdivision on the space range according to the undulation surface of the terrain to determine an equivalent source inversion mesh space;
s3, according to the gravity background field U0Inverting the grid space versus gravity field data d based on the equivalent source0Carrying out PDE three-dimensional inversion calculation with depth normalization factors, positive value constraint terms and normalization terms to obtain a multilayer equivalent source model of the gravity abnormal body;
s4, giving different observation parameters, and performing three-dimensional forward calculation of the gravity field of the nonlinear PDE by using the multilayer equivalent source model to obtain gravity abnormal plumb data and gravity abnormal field data U generated after the gravity abnormal body is extended at different observation positionss' and gravity gradient tensor data.
The technical scheme provided by the invention has the beneficial effects that: the technical scheme provided by the invention can be used for carrying out high-precision polar and data type conversion calculation on the gravity anomaly data generated by the underground gravity anomaly body obtained at the irregular measurement position on the undulating observation curved surface; the invention adopts continuous mesh subdivision, the number of equivalent source layers is usually more than 3, the calculation precision is higher, meanwhile, a depth normalization factor is added, the depth and the range of the equivalent layer do not need to be independently estimated, the gravity anomaly can be subjected to self-adaption, the required data after type conversion and continuation data can be rapidly, efficiently and accurately generated, and the stability and the precision are higher.
Drawings
FIG. 1 is a flowchart of a method for equivalent source extension and data type conversion based on gravity data of a PDE according to an embodiment of the present invention;
fig. 2 is a schematic diagram illustrating an effect of non-uniform mesh generation according to an embodiment of the present invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, embodiments of the present invention will be further described with reference to the accompanying drawings.
Referring to fig. 1, an embodiment of the present invention provides a PDE-based gravity data equivalent source extension and data type conversion method, including the following steps:
s1, acquiring gravity field data d on the undulation observation curved surface0Establishing a topographic relief curved surface according to topographic height information of the area where the gravity field data is located; the gravity field data d0The data may be gravity abnormally vertical data or gravity gradient tensor data, and the gravity abnormally vertical data is taken as an example in this embodiment.
S2, determining a space range of mesh generation according to altitude information of the aerial survey undulation observation surface and the set inversion maximum depth, performing continuous structural non-uniform mesh generation on the space range according to the aerial survey undulation observation surface, and further determining an equivalent source inversion mesh space.
Specifically, the space range of the mesh generation includes an upper top surface and a lower bottom surface, wherein the upper top surface is determined according to the maximum height of the undulation observation curved surface, then the maximum depth of the gravity anomaly which can exist is estimated based on the existing detection technology or actual experience, and the inversion maximum depth is set according to the maximum depth, so that the lower bottom surface is determined;
after the spatial range of mesh generation is determined, dividing the spatial range by the lowest point of the topographic relief surface, carrying out uniform mesh generation on the spatial range above the lowest point to obtain a fine mesh, and carrying out non-uniform expanded mesh generation on the spatial range below the lowest point to obtain an expanded mesh; preferably, if the vertical edge of the fine mesh is 1 length unit (the specific value of the length unit may be set according to the size of the observation area space, for example, 100m is 1 length unit), the vertical edge of the extended mesh is increased at a speed 1.2 times that of the vertical edge of the fine mesh, and the maximum speed increase is set to 1.5 times, so that the calculation amount is reduced on the basis of ensuring a certain inversion accuracy. Please refer to fig. 2, which is a schematic diagram of a result of non-uniform mesh generation in the present embodiment, wherein a spatial range from a relief surface to a bottom surface forms an equivalent source inversion mesh space.
S3, according to the gravity background field U0Inversion of grid space gravity based on equivalent sourcesField data d0And carrying out PDE three-dimensional inversion calculation with depth normalization factors, positive value constraint terms and normalization terms to obtain a multilayer equivalent source model of the gravity abnormal body, wherein the multilayer equivalent source model specifically means that an inversion model solving space comprises a plurality of model depth surfaces, and based on the scheme, the number of model depth surfaces obtained by solving is usually more than 3 layers.
Preferably, the objective function of the PDE-based three-dimensional inversion calculation is:
Figure GDA0003514466110000041
wherein the content of the first and second substances,
Figure GDA0003514466110000042
Us=F(U0,m)
m≥0
where phi denotes an optimization objective (i.e., error),
Figure GDA0003514466110000051
a numerical constraint representing the objective function is represented,
Figure GDA0003514466110000052
representing model constraint of a target function, wherein m represents a density matrix of a multilayer equivalent source model to be solved, and considering that the physical property of an object leads the density to be a positive value, carrying out positive value constraint, namely m is more than or equal to 0; f (-) represents the three-dimensional PDE forward calculation of the multilayer equivalent source model, and can adopt a finite volume or finite element method, UsRepresenting gravity abnormal field data obtained by forward operation, and T (-) representing a conversion function from the gravity abnormal field data to the gravity abnormal plumb data; qx、Qy、QzRespectively representing interpolation functions in the north direction, the east direction and the vertical direction, wherein the interpolation functions comprise observation position information, and a Krigin interpolation function and the like can be selected; u shape0Representing the gravity background field intensity vector, U0x、U0y、U0zRespectively showing its north, east and vertical directionsA directional component; beta represents a regularization factor added according to actual requirements, and if not required, beta is 1; m isrefDensity matrix, W, representing a reference modelrRepresents the depth normalization factor:
Figure GDA0003514466110000053
wherein z represents the distance from the equivalent source to the relief surface, z0Representing the relief surface height, and r represents the depth factor, typically 3.
In this embodiment, the forward computation adopts a finite volume PDE method, and it should be noted that, in order to satisfy the finite volume solution condition, the equivalent source inversion grid space needs to be expanded, please refer to fig. 2, and the grid expansion is performed on the horizontal space below the set inversion maximum depth, above the undulation observation curved surface, and around the equivalent source inversion grid space according to the finite volume method, so as to meet the solution requirement.
In the forward calculation, the gravity field PDE equation is:
Figure GDA0003514466110000054
wherein U is U ═ U0+UsGamma denotes a gravitational constant, and ρ denotes the density of the multi-layered equivalent source model; based on the equation, gravity abnormal field data U is obtainedsForward calculation formula of (c):
Us=F(U0,m);
and (3) iteratively solving an objective function, namely minimizing the error phi, wherein a new density matrix m (namely the distribution of the density rho) is obtained after each iteration is completed and is used for fitting the measured gravity field data d0And finally obtaining a density matrix m of the multilayer equivalent source model. The forward calculation obtains three-component data, and the three-component data can be directly converted into plumb data or gradient tensor data in the inversion process so as to meet the requirements of different observation data types.
S4, given different observation parameters, using multiple layersThe equivalent source model performs three-dimensional forward calculation of the gravity field of the nonlinear PDE to obtain gravity abnormal plumb data and gravity abnormal field data U generated after the gravity abnormal body is extended at different observation positionssAnd one or more of gravity gradient tensor data.
Performing forward calculation according to the density matrix m determined in the step S3 to obtain extended gravity field data Us′:
Us′=F(U0′,m);
The formula for calculating the gravitational field tensor conversion is as follows:
Figure GDA0003514466110000061
wherein the content of the first and second substances,
Figure GDA0003514466110000062
represents a gradient operator, U ═ Ux,Uy,Uz]TAnd representing the gravity field data subjected to the gravity field tensor conversion, wherein factors in a rightmost matrix of the equation are different tensors obtained by converting the gravity field data.
While the present invention has been described with reference to the embodiments shown in the drawings, the present invention is not limited to the embodiments, which are illustrative and not restrictive, and it will be apparent to those skilled in the art that various changes and modifications can be made therein without departing from the spirit and scope of the invention as defined in the appended claims.

Claims (7)

1. A gravity data equivalent source continuation and data type conversion method based on PDE is characterized by comprising the following steps:
s1, acquiring gravity field data d on the undulation observation curved surface0Establishing a topographic relief curved surface according to the topographic height information of the area where the gravity field data is located; the undulation observation curved surface is a recorded known value during observation;
s2, determining a space range of mesh subdivision according to the elevation information of the undulation observation surface and the preset inversion maximum depth, and performing continuous structural non-uniform mesh subdivision on the space range according to the undulation surface of the terrain to determine an equivalent source inversion mesh space;
s3, according to the gravity background field U0Inverting the grid space versus gravity field data d based on the equivalent source0Carrying out PDE three-dimensional inversion calculation with depth normalization factors, positive value constraint terms and normalization terms to obtain a multilayer equivalent source model of the gravity abnormal body;
s4, giving different observation parameters, and performing three-dimensional forward calculation of the gravity field of the nonlinear PDE by using the multilayer equivalent source model to obtain gravity abnormal plumb data and gravity abnormal field data U generated after the gravity abnormal body is extended at different observation positionss' and gravity gradient tensor data;
the number of layers of the model depth surface of the multilayer equivalent source model is more than 3;
in step S3, the objective function of the PDE three-dimensional inversion calculation is shown in formula (1):
Figure FDA0003514466100000011
in the formula (1), the reaction mixture is,
Figure FDA0003514466100000012
Us=F(U0,m)
m≥0
phi denotes an optimization objective, i.e. an error;
Figure FDA0003514466100000013
a numerical constraint representing an objective function;
Figure FDA0003514466100000014
model constraints representing the objective function, m representing the value ofA density matrix of a multilayer equivalent source model solved by iterative optimization; m isrefA density matrix representing a reference model; f (-) represents the three-dimensional PDE forward calculation of the multilayer equivalent source model, UsRepresenting gravity abnormal field data obtained by forward operation, and T (-) representing a conversion function from the gravity abnormal field data to the gravity abnormal plumb data; qx、Qy、QzRespectively representing interpolation functions in the north direction, the east direction and the vertical direction, wherein the interpolation functions comprise observation position information; u shape0Representing the gravity background field intensity vector, U0x、U0y、U0zRespectively representing the north, east and vertical components; beta represents a preset weight parameter;
the depth normalization factor is as shown in formula (2):
Figure FDA0003514466100000021
in the formula (2), z represents the distance from the equivalent source to the relief surface, and z0Representing the relief surface height and r representing the depth factor.
2. The PDE-based gravity data equivalent source extension and data type conversion method according to claim 1, wherein in step S2, the spatial range of the mesh subdivision includes an upper top surface and a lower bottom surface, wherein the upper top surface is a plane determined by a maximum height of the undulation observation surface, and the lower bottom surface is a plane determined by a preset inversion maximum depth.
3. The PDE-based gravity data equivalent source continuation and data type conversion method according to claim 2, wherein in step S2, according to the topographic relief surface, continuous structured non-uniform mesh subdivision is performed on the spatial range to determine an equivalent source inversion mesh space, specifically:
dividing the space range according to the lowest point of the relief curved surface;
carrying out uniform mesh subdivision on the space range above the lowest point to obtain a fine mesh;
carrying out non-uniform mesh subdivision on the space range below the lowest point to obtain an expanded mesh;
and the spatial range from the topographic relief surface to the lower bottom surface, namely the fine grid and the expanded grid form an equivalent source inversion grid space.
4. The PDE-based gravity data equivalent source continuation and data type conversion method according to claim 3, wherein said vertical edge of said extended grid is alpha to said vertical edge of said fine grid1The speed is increased by multiple times and the maximum speed increase is set to alpha2Wherein α is21>1。
5. The method for PDE-based gravity data equivalent source extension and data type conversion according to claim 1, wherein in step S3, the density matrix m and the gravity background field U are generated iteratively during optimization0Performing PDE forward calculation to obtain gravity anomaly field data generated by the gravity anomaly body:
Us=F(U0,m)。
6. the method for PDE-based gravity data equivalent source extension and data type conversion according to claim 1, wherein in step S4, forward calculation is performed according to the density matrix m determined in step S3 to obtain extended gravity anomaly field data Us′:
Us′=F(U0′,m)。
7. The PDE-based gravity data equivalent source extension and data type conversion method according to claim 1, wherein the calculation of the gravity field tensor data is as follows (3):
Figure FDA0003514466100000031
wherein the content of the first and second substances,
Figure FDA0003514466100000032
represents a gradient operator, U ═ Ux,Uy,Uz]TAnd (3) representing the gravity field data subjected to the gravity field tensor conversion, wherein the factors in the rightmost matrix of the formula (3) are different tensor data obtained by converting the gravity field data.
CN202011103444.XA 2020-10-15 2020-10-15 Gravity field data equivalent source continuation and data type conversion method based on PDE Active CN112363236B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011103444.XA CN112363236B (en) 2020-10-15 2020-10-15 Gravity field data equivalent source continuation and data type conversion method based on PDE

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011103444.XA CN112363236B (en) 2020-10-15 2020-10-15 Gravity field data equivalent source continuation and data type conversion method based on PDE

Publications (2)

Publication Number Publication Date
CN112363236A CN112363236A (en) 2021-02-12
CN112363236B true CN112363236B (en) 2022-04-01

Family

ID=74506805

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011103444.XA Active CN112363236B (en) 2020-10-15 2020-10-15 Gravity field data equivalent source continuation and data type conversion method based on PDE

Country Status (1)

Country Link
CN (1) CN112363236B (en)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112949133B (en) * 2021-03-08 2024-03-22 中国地质大学(武汉) Heavy magnetic joint inversion method based on PDE
CN113587921B (en) * 2021-05-31 2024-03-22 中国人民解放军61540部队 Gravity gradient field and gravity anomaly field submersible vehicle fusion positioning method and system
CN113627051B (en) * 2021-07-23 2024-01-30 中国地质科学院地球物理地球化学勘查研究所 Gravity anomaly field separation method, system, storage medium and electronic equipment
CN114167511B (en) * 2021-11-26 2023-08-11 兰州大学 Bit field data rapid inversion method based on continuous expansion downward continuation
CN114814967B (en) * 2022-04-26 2024-04-02 中国人民解放军61540部队 High-resolution submarine topography nonlinear method for inversion of local sea area disturbance gravity data
CN116047617B (en) * 2023-03-10 2023-06-27 中国地质科学院地球物理地球化学勘查研究所 Method and device for identifying geological features between wells

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP3033689B1 (en) * 2013-04-24 2021-09-15 The University of British Columbia A penalty method for pde-constrained optimization
WO2019013906A1 (en) * 2017-07-13 2019-01-17 Exxonmobil Upstream Research Company Visco-pseudo-elastic tti fwi/rtm formulation and implementation
CN107402409A (en) * 2017-09-26 2017-11-28 西南石油大学 A kind of three-dimensional irregular stratum fluctuating interface gravity forward modeling method
CN109375280B (en) * 2018-12-10 2020-03-31 中南大学 Gravity field rapid high-precision forward modeling method under spherical coordinate system
CN110221344B (en) * 2019-06-17 2020-08-28 中国地质大学(北京) Seismic full-waveform and gravity joint inversion method for three-dimensional density structure of earth crust
CN110389391B (en) * 2019-08-01 2020-12-15 自然资源部第二海洋研究所 Heavy magnetic bit field analytic extension method based on spatial domain
CN111142169A (en) * 2020-02-25 2020-05-12 中国地质大学(北京) Submarine topography inversion method based on gravity gradient data

Also Published As

Publication number Publication date
CN112363236A (en) 2021-02-12

Similar Documents

Publication Publication Date Title
CN112363236B (en) Gravity field data equivalent source continuation and data type conversion method based on PDE
CN105549106B (en) A kind of gravity multiple solutions inversion method
CN112147709B (en) Gravity gradient data three-dimensional inversion method based on partial smoothness constraint
CN111856599B (en) Magnetic measurement data equivalent source pole and type conversion method based on PDE
US9915742B2 (en) Method and system for geophysical modeling of subsurface volumes based on label propagation
CN111856598B (en) Magnetic measurement data multilayer equivalent source upper extension and lower extension method
WO2009092992A1 (en) Geophysical data processing systems
WO2014099204A1 (en) Method and system for geophysical modeling of subsurface volumes based on computed vectors
WO2014099201A1 (en) Geophysical modeling of subsurface volumes based on horizon extraction
CN109902315B (en) Method for delineating deep boundary of hidden granite rock mass
CN112346139B (en) Gravity data multilayer equivalent source continuation and data conversion method
CN114943178A (en) Three-dimensional geological model modeling method and device and computer equipment
CN112526625B (en) Computing device for abnormal value of Bragg gravity of aviation gravity measurement point
CN115437027A (en) Method and device for calculating bump gravity anomaly by using geological information variable density forward modeling
CN111859251B (en) Magnetic measurement data equivalent source extension and extension method based on PDE
CN109444973B (en) Gravity forward modeling acceleration method under spherical coordinate system
US20160334537A1 (en) System and method for two dimensional gravity modeling with variable densities
CN107942374A (en) Diffracted wave field extracting method and device
CN112666614B (en) Debris flow source static reserve calculation method based on electrical prospecting and digital elevation model
EP3281044B1 (en) Method for estimating elastic parameters of subsoil
CN111208558B (en) Method and device for establishing ultra-deep low-amplitude three-dimensional geological structure
CN112748471B (en) Gravity-magnetic data continuation and conversion method of unstructured equivalent source
WO2012021938A1 (en) A method of analysing data obtained using a gravity gradiometer
CN107730586B (en) Method and system for modeling stratum
CN114966878B (en) Three-dimensional gravity field inversion method based on mixed norm and cross-correlation constraint

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