CN112596106A - Method for gravity-seismic joint inversion of density interface distribution in spherical coordinate system - Google Patents

Method for gravity-seismic joint inversion of density interface distribution in spherical coordinate system Download PDF

Info

Publication number
CN112596106A
CN112596106A CN202011240182.1A CN202011240182A CN112596106A CN 112596106 A CN112596106 A CN 112596106A CN 202011240182 A CN202011240182 A CN 202011240182A CN 112596106 A CN112596106 A CN 112596106A
Authority
CN
China
Prior art keywords
interface
depth
gravity
inversion
density
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.)
Granted
Application number
CN202011240182.1A
Other languages
Chinese (zh)
Other versions
CN112596106B (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.)
Fifth Gold Detachment Of Chinese People's Armed Police Force
Original Assignee
Fifth Gold Detachment Of Chinese People's Armed Police Force
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 Fifth Gold Detachment Of Chinese People's Armed Police Force filed Critical Fifth Gold Detachment Of Chinese People's Armed Police Force
Priority to CN202011240182.1A priority Critical patent/CN112596106B/en
Publication of CN112596106A publication Critical patent/CN112596106A/en
Application granted granted Critical
Publication of CN112596106B publication Critical patent/CN112596106B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G01V2210/622—Velocity, density or impedance
    • G01V2210/6224—Density

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

The invention provides a method for gravity-seismic joint inversion of density interface distribution in a spherical coordinate system, which comprises the following steps: estimating the initial depth of the density interface, and calculating the corrected depth of the density interface
Figure DDA0002768102320000011
Post-iteration gravity anomaly matrix values
Figure DDA0002768102320000012
Calculating, namely calculating the difference value delta h between the known seismic constraint point and the inversion calculation pointiRoot mean square error analysis, multiple iterative calculation and iterative end point judgment. The invention has the beneficial effects that: aiming at regional large-scale density interface inversion research, the method approximates the earth to be a sphere, and uses spherical unit bodies to replace conventional prisms, so that the problem that the curvature of the earth is neglected in large-scale research is solved, the inversion result is closer to the real condition, and meanwhile, a depth weighting factor and the ground are introduced in the inversion processAnd the inversion accuracy is improved by seismic restraint, and the multi-solution and uncertainty of density interface inversion are reduced.

Description

Method for gravity-seismic joint inversion of density interface distribution in spherical coordinate system
Technical Field
The invention belongs to the technical field of geophysical exploration, and particularly relates to a method for gravity-seismic joint inversion of density interface distribution in a spherical coordinate system.
Background
In the field of geophysical exploration, density interface inversion has important significance in the aspects of deducing regional structure research, acquiring fluctuation of a crystal basal plane and a Mohol plane and the like. The density interface inversion refers to inverting the depth of the density interface according to the abnormal change of gravity. There are many inversion methods for the depth of the interface at the present stage, which can be broadly divided into a space-domain method and a frequency-domain method. The space domain mainly comprises a direct iteration method, a compressed texture and surface method and a nonlinear inversion method, and the space domain inversion interface generally has discontinuity and is slow in calculation speed; the frequency domain density interface algorithm introduces Fourier transform, and improves the calculation efficiency.
At present, a small-range region Dec density interface inversion method is mature, and the method divides an underground model into regular prisms in a Cartesian coordinate system. However, when regional large scale or global scale problems are involved (such as regional mojohn plane inversion, regional tectonic geological interpretation, etc.), the curvature of the earth needs to be taken into account (Anderson.1976; Grombein et al, 2013; Uieda at al, 2016, Aowei Hao et al.2019), and the result obtained by the method has a certain deviation from the real situation. And the existence of multiple solutions in the gravity interface inversion is also a technical problem which is difficult to solve by the prior art.
Disclosure of Invention
In order to make up for the defects of the prior art, the invention provides a method for gravity-seismic joint inversion of density interface distribution in a spherical coordinate system.
A method for gravity-seismic joint inversion of density interface distribution in a spherical coordinate system comprises the following steps:
(1) estimating the initial density interface depth: according to geology or related physical data, giving out initial density interface reference depth, namely depth value of gravity anomaly datum plane, and giving out g which is a gravity anomaly matrix under a spherical coordinate system to be invertedobs(q)While presenting known seismic-constrained point data
Figure BDA0002768102300000021
Dividing the subsurface density layer into M and N spherical unit bodies in longitude and latitude, wherein
Figure BDA0002768102300000022
The depth values of the ith ball unit body in the longitude direction and the jth ball unit body in the latitude direction relative to the gravity anomaly datum surface after the kth iteration are obtained; the observation matrix of the sphere unit body has M multiplied by N points in total, and the gravity anomaly of each measuring point is approximate to that generated by an infinite flat plate, so that the depth of an initial density interface can be estimated
Figure BDA0002768102300000023
Figure BDA0002768102300000024
Wherein
Figure BDA0002768102300000025
G is a universal gravitation constant, and delta rho is a density difference value of an upper interface and a lower interface of a density layer;
(2) calculating corrected density interface depth
Figure BDA0002768102300000026
According to the estimated initial density interface depth
Figure BDA0002768102300000027
Adopting a forward modeling method for gravity anomaly of a fluctuation density interface under a spherical coordinate system to obtain a gravity anomaly value of each measuring point under the spherical coordinate system, and further obtaining an initial gravity anomaly matrix of
Figure BDA0002768102300000028
Residual of gravity at this time
Figure BDA0002768102300000029
Depth at the initial density interface
Figure BDA00027681023000000210
Based on the depth weighting function w (delta r) of the seismic data, the depth weighting function w (delta r) is inverted by introducing a lower interface of a spherical coordinate, and known seismic point constraint information is introduced based on a seismic data constraint method
Figure BDA00027681023000000211
Respectively correcting the gravity iteration correction term and the seismic data correction term to obtain the density interface depth after the 1 st iteration
Figure BDA00027681023000000212
Figure BDA00027681023000000213
Similarly, the abnormal residual of gravity for the kth iteration
Figure BDA00027681023000000214
After correcting the gravity iteration correction term and the seismic data correction term, obtaining the density interface depth after the k iteration
Figure BDA00027681023000000215
Figure BDA00027681023000000216
Where α and β are convergence constraints for gravity and earthquake, respectively, and the general range in practical testing is α e (11.5), β e (0.751.25).
(3) Gravity anomaly matrix values after kth iteration
Figure BDA0002768102300000031
And (3) calculating: adopting a fluctuating density interface gravity anomaly forward modeling method under spherical coordinates to carry out the depth of the density interface after the k-th iteration
Figure BDA0002768102300000032
The gravity anomaly forward formula of the fluctuation density interface under the spherical coordinates is introduced, and the gravity anomaly matrix value after the kth iteration is obtained through updating
Figure BDA0002768102300000033
(4) Calculating the difference value delta h between the known seismic constraint point and the inversion calculation pointi: connecting the known seismic tie points
Figure BDA0002768102300000034
And inversely calculated depth points
Figure BDA0002768102300000035
Substituting the following formula for the difference Δ hiAnd (3) calculating:
Figure BDA0002768102300000036
wherein
Figure BDA0002768102300000037
And
Figure BDA0002768102300000038
respectively a known seismic restraint point and an inversion-calculated depth point, i belongs to n; wherein the inversion isCalculated depth point
Figure BDA0002768102300000039
Depth of density interface through each iteration
Figure BDA00027681023000000310
Obtaining;
(5) root mean square error analysis: computing gravity anomaly matrix values after each iteration
Figure BDA00027681023000000311
And g is a gravity anomaly matrix under a spherical coordinate system to be invertedobs(q)Root mean square error RMS of gravity fit therebetweengraAnd calculating the difference value delta h between the known seismic constraint point and the inversion calculation pointiFitting root mean square error RMS of seismic tie pointsseisThe formula is as follows
Figure BDA00027681023000000312
(6) Multiple iterative computation and iterative end point judgment: when the gravity is fitted to the root mean square error RMSgraAnd fitting root mean square error RMS to the seismic tie pointsseisWhen the negligible minimum value is reached or the iteration frequency reaches the upper limit, the iteration is terminated, and the density interface inversion result is output, wherein the result is obtained at the moment
Figure BDA00027681023000000313
Can be regarded as the density interface depth value of the ith ball unit body in the longitude direction and the jth ball unit body in the latitude direction relative to the gravity reference surface; otherwise, repeating the steps (2) to (5) to carry out iterative calculation until the gravity is fitted with the root mean square error RMSgraFitting root mean square error RMS to seismic tie pointsseisAnd when the negligible minimum value is reached or the iteration frequency reaches the upper limit, the iteration is terminated, and the density interface inversion result is output.
Further, the method for forward modeling gravity anomaly of the fluctuation density interface under the spherical coordinates in step (2) comprises the following steps: the subsurface density layer is divided into longitude anddividing the latitude into M and N sphere unit bodies, wherein a measuring point P is any measuring point on an observation surface, and the gravity anomaly g of the fluctuation density interface of the underground rock stratum under the sphere coordinate caused by the measuring point PpComprises the following steps:
Figure BDA0002768102300000041
wherein Forward (Δ r)ij) The gravity anomaly of the fluctuating density interface under the spherical coordinates caused by the ith spherical unit body in the longitudinal direction and the jth spherical unit body in the latitudinal direction is shown.
Further, the computing principle and method of the interface inversion depth weighting function w (Δ r) under the spherical coordinates in step (2) are as follows:
gravity anomaly is attenuated rapidly along with the deepening of depth, gravity signals are insensitive to the deep part, the depth weighting function is introduced to balance the weights of different depths for earth surface gravity observation, so that an inversion result is closer to a real underground model space, distortion caused by different depths of inversion is corrected, the inversion can adapt to the condition that interface fluctuation is severe, and the interface inversion depth weighting function is performed under the spherical coordinate of an interface deviation datum plane under a spherical coordinate system:
Figure BDA0002768102300000042
where η is a depth-weighted contraction factor, typically 1 or 2, H is the depth relative to the observation plane in spherical coordinates, r0The reference interface depth radius, r is the actual interface relief radius, Δ r is the difference between the actual interface relief radius and the reference interface depth radius, and is positive along the radius direction, where r is r0+ Δ r; the unit bodies with the same scale in the longitude and latitude directions only have difference in the radius direction, and the volume of the unit bodies can be approximate to the product of the edge of each tomb of the spherical unit bodies, namely
Figure BDA0002768102300000043
The formula can be simplified into
Figure BDA0002768102300000044
When Δ r is 0, the actual interface coincides with the reference interface, w (Δ r) ═ 1; when delta r is larger than 0, the actual interface is positioned above the reference interface, and w (delta r) is smaller than 1, and the interface bump of the shallow part is pressed; when Δ r < 0, the actual interface is located below the reference interface, and when w (Δ r) > 1, the deep recess is protruded.
Further, the seismic data constraint method in the step (2) is as follows:
assuming that the difference of the interface model and the difference of the gravity anomaly have a linear relationship:
Figure BDA0002768102300000045
wherein
Figure BDA0002768102300000051
For the correction between the actual interface and the actual calculation model, Δ gM×NA and b are coefficients of linear fitting for theoretical gravity anomaly and residual error of calculation gravity anomaly;
assuming that n seismic measurement points are used as depth constraint points in the region, the method can be known through least square normal regression:
Figure BDA0002768102300000052
Figure BDA0002768102300000053
Figure BDA0002768102300000054
Figure BDA0002768102300000055
wherein
Figure BDA0002768102300000056
And
Figure BDA0002768102300000057
respectively the observed gravity anomaly and the calculated gravity anomaly of the ith seismic restraint point,
Figure BDA0002768102300000058
and
Figure BDA0002768102300000059
known seismic depth points and inversion-computed depth points, respectively, where i ∈ n.
The invention has the beneficial effects that:
the method for gravity-seismic joint inversion of density interface distribution under the spherical coordinate system balances the weights of different depths for earth surface gravity observation by introducing a method for inverting a depth weighting function by an interface under the spherical coordinate system, so that an inversion result is closer to a real underground model space, distortion caused by inverting different depths is corrected, and inversion can adapt to the condition of severe interface fluctuation.
Secondly, the method for the gravity-seismic joint inversion of density interface distribution in the spherical coordinate system can effectively reduce the existence of multi-solution by introducing known geological and other related information to obtain constraint based on a seismic data constraint method. Compared with gravity, the seismic result can well invert underground interface information, inversion accuracy is improved, conventional known seismic depth point constraints adopt 'hard constraints', the range near a known point is assumed to be known depth, local deformation is easy to grow when seismic depth data are inaccurate, and 'soft constraints' are adopted, so that the accuracy of one or more constraint points is not pursued, and the optimization of an integral inversion model is pursued.
The invention provides a method for gravity-seismic joint inversion of density interface distribution in a spherical coordinate system, which is aimed at the inversion research of a large-scale density interface in an area, approximates the earth to be a sphere, uses a spherical unit body to replace a conventional prism body, overcomes the problem that the curvature of the earth is neglected in the large-scale research, and improves the inversion precision and reduces the multiple solution and uncertainty of density interface inversion, and the inversion result is closer to the real condition.
Drawings
Fig. 1 is a schematic two-dimensional cross-sectional view of a relief density layer divided into spherical unit bodies in a relief density interface gravity anomaly forward method under spherical coordinates in embodiment 1 of the present invention;
FIG. 2 is a flowchart illustrating steps of a method for joint-inversion of density interface distribution by gravity and seismic in a spherical coordinate system according to embodiment 1 of the present invention;
FIG. 3 is a schematic diagram of a test result of a method for joint inversion density interface distribution of heavy shocks in a spherical coordinate system according to embodiment 2 of the present invention; wherein (a) is a theoretical simulated mohuo interface profile; (b) obtaining a theoretical gravity anomaly distribution map by a fluctuation density interface gravity anomaly forward modeling method under spherical coordinates; (c) fitting root mean square error RMS for iteration number and gravitygraThe convergence graph of (a); (d) a RMSseis convergence plot for iteration number and seismic tie points; (e) after 7 iterations, comparing the difference between the interface result and the abnormal result of the seismic restraint point; (f) computing a residual δ g between an outlier and a theoretical gravity anomaly for a known gravitykA schematic diagram; (g) a residual error graph of a theoretical density interface distribution result and an inversion density interface result is obtained; (h) the density interface inversion result distribution result after 7 times of iteration is obtained.
Detailed Description
In order to make the objects, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention are clearly and completely described, and it is obvious that the described embodiments are a part of the embodiments of the present invention, but not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
Example 1
1. Principle and method for forward modeling gravity anomaly of fluctuation density interface under spherical coordinate
Dividing the underground density layer into M and N spherical unit bodies in longitude and latitude, wherein a measuring point P is any measuring point on an observation surface (the two-dimensional section schematic diagram of the fluctuation density layer divided into the spherical unit bodies is shown in figure 1), and the gravity anomaly g of the fluctuation density interface under the spherical coordinates caused by the underground rock layer at the measuring point PpComprises the following steps:
Figure BDA0002768102300000071
wherein Forward (Δ r)ij) The gravity anomaly of the fluctuating density interface under the spherical coordinates caused by the ith spherical unit body in the longitudinal direction and the jth spherical unit body in the latitudinal direction is shown.
2. Computing principle and method for interface inversion depth weighting function w (delta r) under spherical coordinates
Gravity anomaly is rapidly attenuated along with the deepening of depth, gravity signals are insensitive to the deep part, the skin effect can grow in conventional inversion, the depth weighting function is introduced to balance the weights of different depths for observing the gravity of the earth surface, so that the inversion result is closer to the space of a real underground model, the distortion caused by different depths of inversion is corrected, the inversion can adapt to the condition that the fluctuation of an interface is severe, and the depth weighting function is inverted by the interface under the spherical coordinate of the spherical coordinate system, wherein the interface deviates from the reference surface:
Figure BDA0002768102300000072
where η is a depth-weighted contraction factor, typically 1 or 2, H is the depth relative to the observation plane in spherical coordinates, r0The reference interface depth radius, r is the actual interface relief radius, Δ r is the difference between the actual interface relief radius and the reference interface depth radius, and is positive along the radius direction, where r is r0+ Δ r; the unit bodies with the same scale in the longitude and latitude directions only have the difference in the radius directionThe volume of a unit cell can be approximated as the product of the edges of a sphere unit cell, i.e.
Figure BDA0002768102300000073
The formula can be simplified into
Figure BDA0002768102300000074
When Δ r is 0, the actual interface coincides with the reference interface, w (Δ r) ═ 1; when delta r is larger than 0, the actual interface is positioned above the reference interface, and w (delta r) is smaller than 1, and the interface bump of the shallow part is pressed; when delta r is less than 0, the actual interface is positioned below the reference interface, w (delta r) > 1 at the time, the deep depression is protruded, shallow protrusion and deep depression during inversion can be suppressed through a depth weighting function under a spherical coordinate system, and the accuracy of interface inversion is improved.
3. Principles and methods for seismic data constraint
The gravity interface inversion has multi-solution, and the addition of known geological and other related information to obtain constraint can effectively reduce the existence of multi-solution. Compared with gravity, the seismic result can well invert the underground interface information, and inversion accuracy is improved. The conventional known seismic depth point constraint adopts 'hard constraint', the range near the known point is assumed to be the known depth, local deformation is easy to grow when the seismic depth data are inaccurate, and the 'soft constraint' is adopted at the time, so that the correctness of one or more constraint points is not pursued, and the optimization of the integral inversion model is pursued.
Assuming that the difference of the interface model and the difference of the gravity anomaly have a linear relationship:
Figure BDA0002768102300000081
wherein
Figure BDA0002768102300000082
For corrections between actual interfaces and actual calculation modelsValue,. DELTA.gM×NA and b are coefficients of linear fitting for theoretical gravity anomaly and residual error of calculation gravity anomaly;
assuming that n seismic measurement points are used as depth constraint points in the region, the method can be known through least square normal regression:
Figure BDA0002768102300000083
Figure BDA0002768102300000084
Figure BDA0002768102300000085
Figure BDA0002768102300000086
wherein
Figure BDA0002768102300000087
And
Figure BDA0002768102300000088
respectively the observed gravity anomaly and the calculated gravity anomaly of the ith seismic restraint point,
Figure BDA0002768102300000089
and
Figure BDA0002768102300000091
known seismic depth points and inversion-computed depth points, respectively, where i ∈ n.
4. A method for joint inversion of density interface distribution by using gravity and earthquake under a spherical coordinate system comprises the following steps (see a flow chart in the step of figure 2):
(1) estimating the initial density interface depth: according to geological or related physical property data, giving out initialThe density interface reference depth is the depth value of the gravity anomaly datum plane, and the gravity anomaly matrix is given as g under a spherical coordinate system to be invertedobs(q)While presenting known seismic-constrained point data
Figure BDA0002768102300000092
Dividing the subsurface density layer into M and N spherical unit bodies in longitude and latitude, wherein
Figure BDA0002768102300000093
The depth values of the ith ball unit body in the longitude direction and the jth ball unit body in the latitude direction relative to the gravity anomaly datum surface after the kth iteration are obtained; the observation matrix of the sphere unit body has M multiplied by N points in total, and the gravity anomaly of each measuring point is approximate to that generated by an infinite flat plate, so that the depth of an initial density interface can be estimated
Figure BDA0002768102300000094
Figure BDA0002768102300000095
Wherein
Figure BDA0002768102300000096
G is a universal gravitation constant, and delta rho is a density difference value of an upper interface and a lower interface of a density layer;
(2) calculating corrected density interface depth
Figure BDA0002768102300000097
According to the estimated initial density interface depth
Figure BDA0002768102300000098
Adopting a forward modeling method for gravity anomaly of a fluctuation density interface under a spherical coordinate system to obtain a gravity anomaly value of each measuring point under the spherical coordinate system, and further obtaining an initial gravity anomaly matrix of
Figure BDA0002768102300000099
Residual of gravity at this time
Figure BDA00027681023000000910
At the initial density interface depth
Figure BDA00027681023000000911
Based on the depth weighting function w (delta r) of the seismic data, the depth weighting function w (delta r) is inverted by introducing a lower interface of a spherical coordinate, and known seismic point constraint information is introduced based on a seismic data constraint method
Figure BDA00027681023000000912
Respectively correcting the gravity iteration correction term and the seismic data correction term to obtain the density interface depth after the 1 st iteration
Figure BDA00027681023000000913
Figure BDA00027681023000000914
Similarly, the abnormal residual of gravity for the kth iteration
Figure BDA00027681023000000915
After correcting the gravity iteration correction term and the seismic data correction term, obtaining the density interface depth after the k iteration
Figure BDA00027681023000000916
Figure BDA00027681023000000917
Where α and β are convergence constraints for gravity and earthquake, respectively, and the general range in practical testing is α e (11.5), β e (0.751.25).
(3) Gravity anomaly matrix values after kth iteration
Figure BDA0002768102300000101
And (3) calculating: adopting a fluctuating density interface gravity anomaly forward modeling method under spherical coordinates to carry out the depth of the density interface after the k-th iteration
Figure BDA0002768102300000102
The gravity anomaly forward formula of the fluctuation density interface under the spherical coordinates is introduced, and the gravity anomaly matrix value after the kth iteration is obtained through updating
Figure BDA0002768102300000103
(4) Calculating the difference value delta h between the known seismic constraint point and the inversion calculation pointi: known seismic restraint points
Figure BDA0002768102300000104
And inversely calculated depth points
Figure BDA0002768102300000105
Substituting the following formula for the difference Δ hiAnd (3) calculating:
Figure BDA0002768102300000106
wherein
Figure BDA0002768102300000107
And
Figure BDA0002768102300000108
respectively a known seismic restraint point and an inversion-calculated depth point, i belongs to n; wherein the inversely calculated depth points
Figure BDA0002768102300000109
Depth of density interface through each iteration
Figure BDA00027681023000001010
Obtaining;
(5) root mean square error differenceAnd (3) analysis: computing gravity anomaly matrix values after each iteration
Figure BDA00027681023000001011
And g is a gravity anomaly matrix under a spherical coordinate system to be invertedobs(q)Root mean square error RMS of gravity fit therebetweengraAnd calculating the difference value delta h between the known seismic constraint point and the inversion calculation pointiFitting root mean square error RMS of seismic tie pointsseisThe formula is as follows
Figure BDA00027681023000001012
(6) Multiple iterative computation and iterative end point judgment: fitting the root mean square error RMS as gravitygraFitting root mean square error RMS to seismic tie pointsseisWhen the negligible minimum value is reached or the iteration frequency reaches the upper limit, the iteration is terminated, and the density interface inversion result is output, wherein the result is obtained at the moment
Figure BDA00027681023000001013
Can be regarded as the density interface depth value of the ith ball unit body in the longitude direction and the jth ball unit body in the latitude direction relative to the gravity reference surface; otherwise, repeating the steps (2) to (5) to carry out iterative calculation until the gravity is fitted with the root mean square error RMSgraFitting root mean square error RMS to seismic tie pointsseisAnd when the negligible minimum value is reached or the iteration frequency reaches the upper limit, the iteration is terminated, and the density interface inversion result is output.
5. Density interface distribution test is carried out by adopting the gravity-seismic joint inversion density interface distribution method under the spherical coordinate system
Assuming that the research area ranges from 0 to 15 degrees E, from 0 to 10 degrees N, and the difference between the upper and lower density interfaces is 280kg/m3And constructing the depth of the theoretical Mohuo surface interface as shown in figure 3(a), wherein the change range is 30.2km to 50.6km, the fluctuation change of the model interface is integrally large, the Mohuo surface in the west of the model is deep, and the Mohuo surface in the east of the model is shallow. The grid distance of longitude and latitude of the model subdivision and gravity observation net is 0.1 degrees, and the calculation grid totally comprises 151 degreesThe average depth of the burial depth of a mojo surface is selected to be 40km according to a reference datum plane of gravity anomaly, a gravity observation surface is a curved surface with the height of 0, theoretical gravity anomaly obtained by utilizing a fluctuating density interface gravity anomaly forward method under a spherical coordinate in the gravity-seismic joint inversion density interface distribution method under the spherical coordinate system is shown in a figure 3(b), high anomaly is caused by lifting or invasion of a high-density rock mass at the bottom at a shallow position of the mojo surface, conversely, low anomaly is caused by descending of a low-density rock mass at the upper part of the mojo surface at a deep position of the mojo surface, and the range of the gravity anomaly is-93.43 mGal to 97.02 mGal.
The theoretical gravity anomaly is used as the observed gravity anomaly to carry out space domain iterative inversion, and the iteration times and the gravity are fitted with root mean square error RMSgraSee FIG. 3(c), the number of iterations and seismic tie points fit the root mean square error RMSseisSee fig. 3 (d). The graph shows the results as: converge rapidly at the beginning and stabilize after 6 iterations.
The experiment takes the result after 7 times of iteration (gravity fitting root mean square error RMS)graAt 0.032mGal, seismic tie points were fitted with root mean square error RMSseis87.5 m); FIG. 3(e) is a graph showing the comparison of the difference between the interface result and the abnormal result at the seismic tie points after 7 iterations, with Δ h on the horizontal axisiIs the residual error between the interface inversion result and the known theoretical value, the depth error range is within 250m, and the vertical axis is deltagiIt is the residual result of the gravity theory value and the calculated value of the known earthquake restraint point, and the range is within 0.4 mGal; FIG. 3(f) is a graph of residual δ g between a calculated anomaly value for known gravity and a theoretical gravity anomalykThe overall residual fitting result is better and can be controlled within 0.68 mGal; FIG. 3(g) is a residual plot of the theoretical density interface distribution results and the inverted density interface results, i.e.
Figure BDA0002768102300000111
The residual error can be controlled within 650 m; FIG. 3(h) is the density interface inversion result after 7 iterations of the method
Figure BDA0002768102300000121
The effectiveness of the gravity-seismic joint inversion density interface distribution method under the spherical coordinate system is verified through testing, and the density interface inversion result precision is high. In the inversion process, the inversion of a gravity interface is taken as a main point, and the soft constraint of the known seismic points is added, so that the multi-solution property of the inversion is reduced.
It should be understood that the above-described specific embodiments are merely illustrative of the present invention and are not intended to limit the present invention. Obvious variations or modifications which are within the spirit of the invention are possible within the scope of the invention.

Claims (4)

1. A method for joint inversion of density interface distribution by using gravity and earthquake under a spherical coordinate system is characterized by comprising the following steps:
(1) estimating the initial density interface depth: according to geology or related physical data, giving out initial density interface reference depth, namely depth value of gravity anomaly datum plane, and giving out g which is a gravity anomaly matrix under a spherical coordinate system to be invertedobs(q)While presenting known seismic-constrained point data
Figure FDA0002768102290000011
Dividing the subsurface density layer into M and N spherical unit bodies in longitude and latitude, wherein
Figure FDA0002768102290000012
The depth values of the ith ball unit body in the longitude direction and the jth ball unit body in the latitude direction relative to the gravity anomaly datum surface after the kth iteration are obtained; the observation matrix of the sphere unit body has M multiplied by N points in total, and the gravity anomaly of each measuring point is approximate to that generated by an infinite flat plate, so that the depth of an initial density interface can be estimated
Figure FDA0002768102290000013
Figure FDA0002768102290000014
Wherein
Figure FDA0002768102290000015
G is a universal gravitation constant, and delta rho is a density difference value of an upper interface and a lower interface of a density layer;
(2) calculating corrected density interface depth
Figure FDA0002768102290000016
According to the estimated initial density interface depth
Figure FDA0002768102290000017
Adopting a forward modeling method for gravity anomaly of a fluctuation density interface under a spherical coordinate system to obtain a gravity anomaly value of each measuring point under the spherical coordinate system, and further obtaining an initial gravity anomaly matrix of
Figure FDA0002768102290000018
Residual of gravity at this time
Figure FDA0002768102290000019
Depth at the initial density interface
Figure FDA00027681022900000110
Based on the depth weighting function w (delta r) of the seismic data, the depth weighting function w (delta r) is inverted by introducing a lower interface of a spherical coordinate, and known seismic point constraint information is introduced based on a seismic data constraint method
Figure FDA00027681022900000111
Respectively correcting the gravity iteration correction term and the seismic data correction term to obtain the density interface depth after the 1 st iteration
Figure FDA00027681022900000112
Figure FDA00027681022900000113
Similarly, the abnormal residual of gravity for the kth iteration
Figure FDA00027681022900000114
After correcting the gravity iteration correction term and the seismic data correction term, obtaining the density interface depth after the k iteration
Figure FDA00027681022900000115
Figure FDA00027681022900000116
Wherein α and β are convergence constraint terms of gravity and earthquake, respectively, and the general range in practical test is α e (11.5), β e (0.751.25);
(3) gravity anomaly matrix values after kth iteration
Figure FDA0002768102290000021
And (3) calculating: adopting a fluctuating density interface gravity anomaly forward modeling method under spherical coordinates to carry out the depth of the density interface after the k-th iteration
Figure FDA0002768102290000022
The gravity anomaly forward formula of the fluctuation density interface under the spherical coordinates is introduced, and the gravity anomaly matrix value after the kth iteration is obtained through updating
Figure FDA0002768102290000023
(4) Calculating the difference value delta h between the known seismic constraint point and the inversion calculation pointi: connecting the known seismic tie points
Figure FDA0002768102290000024
And inversely calculated depth points
Figure FDA0002768102290000025
Substituting the following formula for the difference Δ hiAnd (3) calculating:
Figure FDA0002768102290000026
wherein
Figure FDA0002768102290000027
And
Figure FDA0002768102290000028
respectively a known seismic restraint point and an inversion-calculated depth point, i belongs to n; wherein the inversely calculated depth points
Figure FDA0002768102290000029
Depth of density interface through each iteration
Figure FDA00027681022900000210
Obtaining;
(5) root mean square error analysis: computing gravity anomaly matrix values after each iteration
Figure FDA00027681022900000211
And g is a gravity anomaly matrix under a spherical coordinate system to be invertedobs(q)Root mean square error RMS of gravity fit therebetweengraAnd calculating the difference value delta h between the known seismic constraint point and the inversion calculation pointiFitting root mean square error RMS of seismic tie pointsseisThe formula is as follows
Figure FDA00027681022900000212
(6) Multiple iterative computation and iterative end point judgment: when the gravity fits the root mean square errorDifferential RMSgraAnd fitting root mean square error RMS to the seismic tie pointsseisWhen the negligible minimum value is reached or the iteration frequency reaches the upper limit, the iteration is terminated, and the density interface inversion result is output, wherein the result is obtained at the moment
Figure FDA00027681022900000213
Can be regarded as the density interface depth value of the ith ball unit body in the longitude direction and the jth ball unit body in the latitude direction relative to the gravity reference surface; otherwise, repeating the steps (2) to (5) to carry out iterative calculation until the gravity is fitted with the root mean square error RMSgraFitting root mean square error RMS to seismic tie pointsseisAnd when the negligible minimum value is reached or the iteration frequency reaches the upper limit, the iteration is terminated, and the density interface inversion result is output.
2. The method for joint inversion of density interface distribution by using gravity and earthquake in a spherical coordinate system as claimed in claim 1, wherein the method for performing gravity anomaly forward on the fluctuating density interface in the spherical coordinate system in step (2) comprises: dividing the underground density layer into M and N spherical unit bodies in longitude and latitude, wherein a measuring point P is any measuring point on an observation surface, and the fluctuation density interface gravity anomaly g caused by the underground rock layer under the spherical coordinates at the measuring point PpComprises the following steps:
Figure FDA0002768102290000031
wherein Forward (Δ r)ij) The gravity anomaly of the fluctuating density interface under the spherical coordinates caused by the ith spherical unit body in the longitudinal direction and the jth spherical unit body in the latitudinal direction is shown.
3. The method for joint heavy-seismic inversion of density interface distribution in a spherical coordinate system according to claim 2, wherein the computing principle and method of the spherical-coordinate lower interface inversion depth weighting function w (Δ r) in step (2) are as follows:
gravity anomaly is attenuated rapidly along with the deepening of depth, gravity signals are insensitive to the deep part, the depth weighting function is introduced to balance the weights of different depths for earth surface gravity observation, so that an inversion result is closer to a real underground model space, distortion caused by different depths of inversion is corrected, the inversion can adapt to the condition that interface fluctuation is severe, and the interface inversion depth weighting function is performed under the spherical coordinate of an interface deviation datum plane under a spherical coordinate system:
Figure FDA0002768102290000032
where η is a depth-weighted contraction factor, typically 1 or 2, H is the depth relative to the observation plane in spherical coordinates, r0The reference interface depth radius, r is the actual interface relief radius, Δ r is the difference between the actual interface relief radius and the reference interface depth radius, and is positive along the radius direction, where r is r0+ Δ r; the unit bodies with the same scale in the longitude and latitude directions only have difference in the radius direction, and the volume of the unit bodies can be approximate to the product of the edge of each tomb of the spherical unit bodies, namely
Figure FDA0002768102290000033
The formula can be simplified into
Figure FDA0002768102290000034
When Δ r is 0, the actual interface coincides with the reference interface, w (Δ r) ═ 1; when delta r is larger than 0, the actual interface is positioned above the reference interface, and w (delta r) is smaller than 1, and the interface bump of the shallow part is pressed; when Δ r < 0, the actual interface is located below the reference interface, and when w (Δ r) > 1, the deep recess is protruded.
4. The method for the gravity-seismic joint inversion of the density interface distribution in the spherical coordinate system according to claim 3, wherein the seismic data constraint method in the step (2) is as follows:
assuming that the difference of the interface model and the difference of the gravity anomaly have a linear relationship:
Figure FDA0002768102290000041
wherein
Figure FDA0002768102290000042
For the correction between the actual interface and the actual calculation model, Δ gM×NA and b are coefficients of linear fitting for theoretical gravity anomaly and residual error of calculation gravity anomaly;
assuming that n seismic measurement points are used as depth constraint points in the region, the method can be known through least square normal regression:
Figure FDA0002768102290000043
Figure FDA0002768102290000044
Figure FDA0002768102290000045
Figure FDA0002768102290000046
wherein
Figure FDA0002768102290000047
And
Figure FDA0002768102290000048
respectively the observed gravity anomaly and the calculated gravity anomaly of the ith seismic restraint point,
Figure FDA0002768102290000049
and
Figure FDA00027681022900000410
known seismic depth points and inversion-computed depth points, respectively, where i ∈ n.
CN202011240182.1A 2020-11-09 2020-11-09 Method for carrying out gravity and vibration joint inversion on density interface distribution under spherical coordinate system Active CN112596106B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011240182.1A CN112596106B (en) 2020-11-09 2020-11-09 Method for carrying out gravity and vibration joint inversion on density interface distribution under spherical coordinate system

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011240182.1A CN112596106B (en) 2020-11-09 2020-11-09 Method for carrying out gravity and vibration joint inversion on density interface distribution under spherical coordinate system

Publications (2)

Publication Number Publication Date
CN112596106A true CN112596106A (en) 2021-04-02
CN112596106B CN112596106B (en) 2024-04-26

Family

ID=75183428

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011240182.1A Active CN112596106B (en) 2020-11-09 2020-11-09 Method for carrying out gravity and vibration joint inversion on density interface distribution under spherical coordinate system

Country Status (1)

Country Link
CN (1) CN112596106B (en)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113740915A (en) * 2021-07-28 2021-12-03 中国人民武装警察部队黄金第五支队 Method for joint inversion of crust structure parameters by using gravity and receiving function of spherical coordinate system
CN113761457A (en) * 2021-09-09 2021-12-07 中国自然资源航空物探遥感中心 Method for extracting local gravity anomaly based on measured gravity anomaly data
CN117148457A (en) * 2023-08-29 2023-12-01 长安大学 Magnetic layer magnetization modulus calculation method and system
CN117607954A (en) * 2023-11-23 2024-02-27 长安大学 Free parameter and substrate simultaneous inversion method based on index variable density model
CN118467898A (en) * 2024-05-15 2024-08-09 中国自然资源航空物探遥感中心 Method for inverting glacier thickness by utilizing aviation gravity anomaly

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109444973A (en) * 2018-11-06 2019-03-08 中南大学 Gravity forward modeling accelerated method under a kind of spherical coordinate system
CN110045432A (en) * 2018-12-05 2019-07-23 中南大学 Gravitational field forward modeling method and 3-d inversion method under spherical coordinate system based on 3D-GLQ
CN110287620A (en) * 2019-06-28 2019-09-27 中国地震局地球物理研究所 Spherical coordinate system density interface forward modeling method and system suitable for surface observation face
CN111400654A (en) * 2020-03-13 2020-07-10 中南大学 Gravity field rapid forward modeling method and inversion method based on Toplitz nuclear matrix

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109444973A (en) * 2018-11-06 2019-03-08 中南大学 Gravity forward modeling accelerated method under a kind of spherical coordinate system
CN110045432A (en) * 2018-12-05 2019-07-23 中南大学 Gravitational field forward modeling method and 3-d inversion method under spherical coordinate system based on 3D-GLQ
CN110287620A (en) * 2019-06-28 2019-09-27 中国地震局地球物理研究所 Spherical coordinate system density interface forward modeling method and system suitable for surface observation face
CN111400654A (en) * 2020-03-13 2020-07-10 中南大学 Gravity field rapid forward modeling method and inversion method based on Toplitz nuclear matrix

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
LEONARDO UIEDA, VAL´ERIA C.F. BARBOSA: "Fast nonlinear gravity inversion in spherical coordinates with application to the South AmericanMoho", GEOPHYSICAL JOURNAL INTERNATIONAL, vol. 208 *
ZHIWEI LI , TIANYAO HAO, YA XU, YI XU: "An efficient and adaptive approach for modeling gravity effects in spherical coordinates", JOURNAL OF APPLIED GEOPHYSICS, vol. 73 *
张杰;杨光亮;谈洪波;吴桂桔;王嘉沛;: "基于接收函数约束的川滇地区莫霍面深度反演研究", 地球物理学报, vol. 63, no. 07 *
王祥;郭良辉;: "球坐标系密度界面反演方法及在华南大陆的应用", 物探与化探, vol. 44, no. 05 *
陈石;郑秋月;徐伟民;: "南北地震带南段地壳厚度重震联合最优化反演", 地球物理学报, vol. 58, no. 11 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113740915A (en) * 2021-07-28 2021-12-03 中国人民武装警察部队黄金第五支队 Method for joint inversion of crust structure parameters by using gravity and receiving function of spherical coordinate system
CN113740915B (en) * 2021-07-28 2024-04-19 中国人民武装警察部队黄金第五支队 Method for jointly inverting crust structure parameters by spherical coordinate system gravity and receiving function
CN113761457A (en) * 2021-09-09 2021-12-07 中国自然资源航空物探遥感中心 Method for extracting local gravity anomaly based on measured gravity anomaly data
CN113761457B (en) * 2021-09-09 2024-05-14 中国自然资源航空物探遥感中心 Method for extracting local gravity anomaly based on measured gravity anomaly data
CN117148457A (en) * 2023-08-29 2023-12-01 长安大学 Magnetic layer magnetization modulus calculation method and system
CN117607954A (en) * 2023-11-23 2024-02-27 长安大学 Free parameter and substrate simultaneous inversion method based on index variable density model
CN118467898A (en) * 2024-05-15 2024-08-09 中国自然资源航空物探遥感中心 Method for inverting glacier thickness by utilizing aviation gravity anomaly

Also Published As

Publication number Publication date
CN112596106B (en) 2024-04-26

Similar Documents

Publication Publication Date Title
CN112596106A (en) Method for gravity-seismic joint inversion of density interface distribution in spherical coordinate system
AU2020102025A4 (en) Method for Reconstructing Basin Paleogeomorphology
Van Dam et al. Displacements of the Earth's surface due to atmospheric loading: Effects on gravity and baseline measurements
WAN et al. An algorithm of fault parameter determination using distribution of small earthquakes and parameters of regional stress field and its application to Tangshan earthquake sequence
CN102944896B (en) The modelling static correcting method of surface survey data
CN112147709B (en) Gravity gradient data three-dimensional inversion method based on partial smoothness constraint
CN102901985B (en) A kind of Depth Domain interval velocity modification method being applicable to relief surface
CN109799540B (en) Volcanic rock type uranium deposit magnetic susceptibility inversion method based on geological information constraint
CN113031068B (en) Reflection coefficient accurate base tracking prestack seismic inversion method
CN109884710A (en) For the micro logging chromatography imaging method of excitation well depth design
CN110531410A (en) A kind of least square reverse-time migration gradient Preconditioning method based on through wave field
CN113671570A (en) Seismic surface wave travel time and gravity anomaly joint inversion method and system
CN115508908A (en) Seismic surface wave travel time and gravity anomaly joint inversion method and system
CN111025388B (en) Multi-wave combined prestack waveform inversion method
CN110737018B (en) Method for modeling anisotropy of VSP seismic data
CN109490978B (en) Frequency domain rapid high-precision forward modeling method for undulating stratum
CN113740915A (en) Method for joint inversion of crust structure parameters by using gravity and receiving function of spherical coordinate system
CN117572502A (en) Method and system for dynamically correcting seismic profile based on velocity iteration
Sato et al. Local effects on tidal strain measurements at Esashi, Japan
CN114740540B (en) Method and system for constructing magnetic anomaly map of middle ridge area in ocean based on direction constraint
CN109901221B (en) Seismic data anisotropy modeling method based on dynamic correction velocity parameter
CN111538080B (en) Method of seismic imaging
CN115236667A (en) Potential landslide volume estimation method fusing multi-source SAR data
CN108508479A (en) A kind of vacant lot well solid gravity and magnetic data collaboration target location inversion method
Asgari et al. Estimate the Crust Thickness using the Gravity Data for the KopehtDagh Region

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