A multi-scale three-dimensional re-shock joint frequency domain inversion method based on mixed constraints

Abstract

The present application relates to the technical field of geophysical exploration, in particular to a multi-scale three-dimensional re-shock combined frequency domain inversion method based on mixed constraints. The method comprises the following steps: using discrete Fourier transform to determine the observed gravity anomaly on the ground, and using the weighted average method of the values at multiple Gaussian points to discretize the observed gravity anomaly; calculating the spherical harmonics and spherical harmonics coefficients of the two-dimensional density distribution of each layer in the inversion range, and calculating the gravity anomaly on the observation surface based on the spherical harmonics and spherical harmonics coefficients; converting the data in the inversion range into a three-dimensional gravity field forward kernel matrix; constructing an inversion objective function of the three-dimensional gravity field mixed constraint; inputting the pre-acquired three-dimensional seismic velocity model data and converting it into a reference density model; calculating the gradient operator in the inversion objective function, and simplifying the inversion objective function based on the calculation of the gradient operator, and then solving the simplified inversion objective function by using the alternating direction multiplier method.

Description

Technical Field

[0001] This invention relates to the field of geological exploration technology, and more particularly to a multi-scale three-dimensional gravity-seismic joint frequency domain inversion method based on hybrid constraints in geophysical exploration technology. Background Technology

[0002] Gravity exploration, as a geophysical exploration method, has the advantages of simplicity, high cost, and large exploration depth, and is widely used in geological and geophysical problems at various scales. Examples include local-scale mineral and oil and gas resource exploration, geothermal resource development, and regional-scale studies of the Earth's lithosphere crust-mantle density structure, as well as the three-dimensional density structure and formation and evolution history of lunar mass-tumor basins. Gravity field inversion is a crucial step in gravity exploration, involving the inversion of subsurface three-dimensional density distribution based on observed two-dimensional gravity anomaly data.

[0003] However, it is evident that using two-dimensional data to infer a three-dimensional density distribution is an underdetermined problem in mathematics, meaning there are more unknowns than equations, and therefore the solution is not unique. Furthermore, due to the significant volume effect of gravitational fields—that is, according to the law of universal gravitation, observed gravitational anomalies are influenced by both the mass and distance of the anomaly—and since mass is related to density and volume, a shallowly buried but massive sphere and a deeply buried but small sphere may produce the same gravitational field at a single point on the ground. The volume effect further exacerbates the non-uniqueness of the inversion, resulting in two major problems: severe multiple solutions and low depth resolution in three-dimensional gravity inversion.

[0004] To address the issues of high ambiguity and low depth resolution in gravity inversion, various improvements have been made by those skilled in the art. These improvements can be broadly categorized into two approaches: The first approach focuses on leveraging the inherent potential of gravity observation data itself, adding constraints from multiple dimensions to reduce ambiguity. This includes various forms of depth weighting functions, upper and lower bound constraints on physical properties, norm systems, and self-constrained mathematical methods. The second approach involves joint inversion of multi-source data, most commonly joint inversion of gravity and earthquake data, joint inversion of gravity and magnetic data, and joint inversion of gravity gradient tensors. Since earthquake wave velocity is strongly correlated with density, joint inversion of gravity and earthquake data has a robust physical mechanism to support it. Joint inversion of gravity gradient tensors is simpler, based on the fundamental principle of gravity, but with different components complementing each other; it doesn't add unnecessary physical property information, thus its joint effect is less effective than that of joint inversion of gravity and earthquake data. However, these improvements only reduce ambiguity and improve depth resolution to a certain extent; further research is needed to solve the problem of ambiguity in gravity inversion. On the other hand, traditional gravity field inversion methods are limited by the enormous storage requirements of the Jacobian matrix for gravity field inversion and the significant time consumption of forward modeling, making it difficult to solve the problem of large-scale three-dimensional density inversion. Furthermore, most of these methods are based on Cartesian coordinates at the exploration scale, which are ill-suited to the needs of large-scale gravity field inversion in spherical coordinates. Summary of the Invention

[0005] To address the technical problems existing in the prior art, this invention provides a multi-scale three-dimensional gravity field inversion method based on hybrid constraints, which solves the problems of low depth resolution, multiple solutions, and low adaptability of traditional three-dimensional gravity field inversion in the current rectangular coordinate system, and solves the problems mentioned in the background art.

[0006] To achieve the above objectives, the embodiments of the present invention provide the following technical solutions:

[0007] In a first aspect, in one embodiment of the present invention, a multi-scale three-dimensional re-seismic joint frequency domain inversion method based on hybrid constraints is provided, the method comprising the following steps:

[0008] Discrete Fourier transform was used to determine the observed gravity anomalies on the ground, and a weighted average of values ​​at multiple Gaussian points was used to discretize the observed gravity anomalies. The inversion range data was determined based on the discretized observed gravity anomalies.

[0009] Calculate the spherical harmonics and spherical harmonic coefficients of the two-dimensional density distribution in each layer of the inversion range, and calculate the gravity anomaly on the observation surface based on the spherical harmonics and spherical harmonic coefficients;

[0010] The inversion range data is converted into a three-dimensional gravity field forward modeling kernel matrix, and the observed gravity anomaly vector is calculated based on the three-dimensional gravity field forward modeling kernel matrix.

[0011] Construct the inversion objective function for a three-dimensional gravity field with hybrid constraints;

[0012] Input the pre-acquired three-dimensional seismic velocity model data and convert it into a reference density model;

[0013] Calculate the gradient operator in the inversion objective function, and simplify the inversion objective function based on the calculated gradient operator. Then, solve the simplified inversion objective function using the alternating direction multiplier method.

[0014] Continue the above iterative calculations until the convergence criterion is met.

[0015] As a further aspect of the present invention, the expression for observing gravity anomalies is:

[0016]

[0017] In the formula, Δg(x) m ,y n (z0) represents the observed gravity anomaly; Represented as the spectrum; xm represents the m-th point in the x-direction, y n Let z0 represent the nth point in the y-direction, and z0 represent the height of the observation point; k x k y dk represents the wavenumber in the x and y directions. x dk y represents the differential of the wave number in the x and y directions; e is the natural constant with a value of 2.71828; π is the mathematical constant pi.

[0018] As a further aspect of the present invention, the expression for the inversion range data is:

[0019]

[0020] In the formula, e represents the natural exponent, and i represents the imaginary unit. and λ j Let Δk represent the k-th and j-th Gaussian point values ​​in the north-south and east-west directions, respectively. x Δk y k represents the discrete spectral interval. xp and k yq The p-th and q-th discrete wavenumbers are represented in the x and y directions, respectively; Mx and Ny represent the number of segments in the x and y directions, respectively; x m and y n These represent the positions of the m-th point in the x-direction and the n-th point in the y-direction, respectively.

[0021] As a further aspect of the present invention, the calculation of the spherical harmonic coefficients of the two-dimensional density distribution of each layer in the inversion range, and the calculation of gravity anomalies on the observation surface based on the spherical harmonic coefficients, includes:

[0022] Calculate the gravitational potential, and based on the gravitational potential, calculate the spherical harmonics and spherical harmonic coefficients of the observed gravitational anomaly;

[0023] Gravity anomalies on the observation surface are calculated based on spherical harmonics and spherical harmonic coefficients.

[0024] As a further aspect of the present invention, the gravity potential can be expressed as the following formula:

[0025]

[0026] In the formula, V(r,θ,φ) represents the gravitational potential, γ represents the gravitational constant, π represents pi, a and b represent the spherical harmonic order, and r′ represents the radius of the location of the subsurface anomaly. and Let θ represent the spherical harmonic coefficient of the density distribution in the r′th underground layer, where r represents the radius of the point, and θ and λ represent latitude and longitude, respectively.

[0027] As a further aspect of the present invention, the formula for calculating the observed gravity anomaly vector is as follows:

[0028] Kρ=Δg

[0029] In the formula, K represents the three-dimensional gravity field forward modeling kernel matrix, which is an Mx*Ny matrix with elements K ij Let ρ represent the gravitational field response of the j-th unit cell with a density of unit 1 to the i-th observation point, where ρ represents the density vector and Δg represents the observed gravity anomaly vector.

[0030] As a further aspect of the present invention, the calculation formula for the inversion objective function is as follows:

[0031]

[0032] In the formula, Ψ(ρ) is the inversion objective function, ||Δg-K·ρ|| 2 For the data fitting term, ||W ρ (ρ-ρ ref )|| 2 For smooth model constraints, For the gradient operator term, ρ ref This is the reference density model.

[0033] As a further aspect of the present invention, the calculation formula of the gradient operator is as follows:

[0034]

[0035] Among them, D x D y D z Let x, y, and z represent the differential operators in the x, y, and z directions, respectively.

[0036] As a further aspect of the present invention, the simplified formula for calculating the inversion objective function is as follows:

[0037] Ψ(ρ)=||Δg-K·ρ|| 2 +β|Φρ|.

[0038] Φ represents the combined term of the smooth model constraint and the gradient operator constraint.

[0039] As a further aspect of the present invention, the expression for Φ is as follows:

[0040] Φ=[I T ,(D x W) T ,(D y W) T ,(D z W) T ] T

[0041] In the formula, I represents the identity matrix, and the superscript T indicates the transpose operation of the matrix.

[0042] The technical solution provided by this invention has the following beneficial effects:

[0043] This invention is applicable to a novel multi-scale three-dimensional gravity field gravity-seismic joint inversion method in both rectangular and spherical coordinate systems. By constructing a novel inversion objective function, a multi-norm constraint system and a gravity-seismic joint inversion system based on smooth constraint norm and gradient operator focusing constraint are formed, which can effectively solve the problems of low depth resolution and strong multiple solutions in traditional three-dimensional inversion methods.

[0044] These or other aspects of the invention will become more apparent from the following description of embodiments. It should be understood that the foregoing general description and the following detailed description are exemplary and explanatory only, and are not intended to limit the invention. Attached Figure Description

[0045] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other embodiments can be obtained based on these drawings without creative effort.

[0046] Figure 1This is a flowchart of a multi-scale three-dimensional re-seismic joint frequency domain inversion method based on hybrid constraints, according to an embodiment of the present invention.

[0047] Figure 2 This is a comparison image of the actual model in the density inversion cross section of the 3D salt dome model.

[0048] Figure 3 The image shows the results of traditional L2 norm inversion based on smoothness constraints, compared to the cross-section of the density inversion of the 3D salt dome model.

[0049] Figure 4 The image shows the results obtained by the traditional hybrid constraint-based inversion method in comparison with the cross-section of the density inversion of the 3D salt dome model.

[0050] Figure 5 A three-dimensional display of the actual three-dimensional density distribution of salt domes in a three-dimensional display for the inversion and comparison of three-dimensional salt dome density.

[0051] Figure 6 This is a 3D visualization of traditional smooth constraint inversion in a 3D display of 3D salt dome density inversion comparison.

[0052] Figure 7 The three-dimensional display diagram is an inversion diagram of the method proposed in this invention for the three-dimensional salt dome density inversion comparison.

[0053] Figure 8 The depth slice map of the Qinghai-Tibet Plateau at a depth of 10km is obtained by gravity seismic combined three-dimensional density inversion using the method proposed in this invention.

[0054] Figure 9 The depth slice map of the Qinghai-Tibet Plateau at a depth of 40km is obtained by gravity seismic combined three-dimensional density inversion using the method proposed in this invention.

[0055] Figure 10 The depth slice map of the Qinghai-Tibet Plateau at a depth of 140km is obtained by gravity seismic combined three-dimensional density inversion using the method proposed in this invention.

[0056] Figure 11 This is a depth slice of a 3D seismic velocity model at a depth of 10km on the Qinghai-Tibet Plateau. This velocity model is converted into an initial density model and then used in the joint inversion of heavy earthquakes.

[0057] Figure 12 This is a depth slice of a 40km-deep 3D seismic velocity model for the Qinghai-Tibet Plateau. This velocity model is converted into an initial density model and then used in the joint inversion of heavy earthquakes.

[0058] Figure 13 This is a depth slice map of a 3D seismic velocity model at a depth of 140km on the Qinghai-Tibet Plateau. This velocity model is converted into an initial density model and then used in the joint inversion of heavy earthquakes.

[0059] Figure 14 This is a depth profile map of different sections obtained from the 3D gravity-seismic joint inversion of the Tibetan Plateau. The first row shows the horizontal position of the different sections, namely CC′, DD′, EE′, FF′, GG′, and HH′. The second to fourth rows are the density-depth profile maps of these six sections. The dashed lines in the figure represent the Moho depth on the section. Detailed Implementation

[0060] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0061] The flowchart shown in the attached diagram is for illustrative purposes only and does not necessarily include all content and operations / steps, nor does it necessarily have to be performed in the order described. For example, some operations / steps can be broken down, combined, or partially merged, so the actual execution order may change depending on the actual situation.

[0062] It should be understood that the terminology used in this specification is for the purpose of describing particular embodiments only and is not intended to limit the invention. As used in this specification and the appended claims, the singular forms “a,” “an,” and “the” are intended to include the plural forms unless the context clearly indicates otherwise.

[0063] Specifically, the embodiments of the present invention will be further described below with reference to the accompanying drawings.

[0064] Please see Figure 1 , Figure 1 This is a flowchart of a multi-scale three-dimensional re-seismic joint frequency domain inversion method based on hybrid constraints provided by an embodiment of the present invention, as shown below. Figure 1 As shown, the multi-scale three-dimensional re-seismic joint frequency domain inversion method based on hybrid constraints includes steps S10 to S70.

[0065] S10. Discrete Fourier transform was used to determine the observed gravity anomaly on the ground, and a weighted average of values ​​at multiple Gaussian points was used to discretize the observed gravity anomaly. Based on the discretized observed gravity anomaly, the inversion range data was determined.

[0066] In an embodiment of the present invention, the method of determining the observed gravity anomaly on the ground using discrete Fourier transform, and discretizing the observed gravity anomaly using a weighted average of values ​​at multiple Gaussian points, and determining the inversion range data based on the discretized observed gravity anomaly, includes:

[0067] Let the north-south range be 0–X, the east-west range be 0–Y, and the depth range be 0–Z. Divide the north-south, east-west, and depth directions into segments Nx, Ny, and Nz, respectively, with intervals Δx, Δy, and Δz. Assume the position indices for each segment are p = 1, 2, ..., Nx, q = 1, 2, ..., Ny, and l = 1, 2, ..., Nz. Assume the density of each subdivided block is constant, ρ(ξ, η, ζ), where (ξ, η, ζ) represents the coordinates of the subdivided block. The coordinates of the observation point are (x, y, z0), and the horizontal indexes are m = 1, 2, ..., Nx and n = 1, 2, ..., Ny. Then, for small-scale forward modeling of the gravity field in a Cartesian coordinate system, the influence of the Earth's curvature can be disregarded. Therefore, the observed gravity anomaly Δg on the ground can be expressed as:

[0068]

[0069] Wherein, Δg(x) m ,y n (z0) represents the observed gravity anomaly; Represented as the spectrum; x m Let m be the m-th point in the x-direction, and y be the m-th point in the x-direction. n Let z0 represent the nth point in the y-direction, and z0 represent the height of the observation point; k x k y dk represents the wavenumber in the x and y directions. x dk y Let represent the differentials of the wavenumbers in the x and y directions; e is the natural constant with a value of 2.71828; and π is the mathematical constant pi. The above equation is a continuous inverse Fourier transform, which cannot be implemented on a computer. After discretization, it can be expressed as:

[0070]

[0071] Where K and J represent the Gaussian points in the north-south and east-west directions, respectively, and k and j represent the indices of their positions. e represents the natural exponent, and i represents the imaginary unit. and λ j Let Δk represent the k-th and j-th Gaussian point values ​​in the north-south and east-west directions, respectively. x Δk y k represents the discrete spectral interval. xp and k yq The p-th and q-th discrete wavenumbers are represented in the x and y directions, respectively; Mx and Ny represent the number of segments in the x and y directions, respectively; x m and y n Let A represent the position of the m-th point in the x-direction and the position of the n-th point in the y-direction, respectively. kj for:

[0072]

[0073] Formulas (1)-(3) above realize the frequency domain forward modeling of the three-dimensional gravity field. Due to the use of discrete Fourier transform and the weighted average of values ​​at multiple Gaussian points, the calculation efficiency is improved and the calculation accuracy is increased compared with the traditional spatial domain forward modeling method.

[0074] S20. Calculate the spherical harmonic coefficients of the two-dimensional density distribution of each layer in the inversion range, and calculate the gravity anomaly on the observation surface based on the spherical harmonic coefficients.

[0075] In an embodiment of the present invention, the calculation of the spherical harmonic coefficients of the two-dimensional density distribution of each layer in the inversion range, and the calculation of gravity anomalies on the observation surface based on the spherical harmonic coefficients, includes:

[0076] Calculate the gravitational potential, and based on the gravitational potential, calculate the spherical harmonics and spherical harmonic coefficients of the observed gravitational anomaly;

[0077] Gravity anomalies on the observation surface are calculated based on spherical harmonics and spherical harmonic coefficients.

[0078] It should be noted that before calculating the gravity potential, it is also necessary to perform a subdivision step in the rectangular coordinate system, just like in step S10, and the subdivision step is the same as in step S10. However, the coordinates of the observation point are (r, θ, λ), where r represents the radius of the point, and θ and λ represent latitude and longitude, respectively.

[0079] It should be noted that the gravitational potential V(r,θ,φ) can be expressed by the following formula:

[0080]

[0081]

[0082] Where γ represents the gravitational constant, π represents pi, a and b represent the spherical harmonic order, and r′ represents the radius of the location of the subsurface anomaly. and Let θ represent the spherical harmonic coefficient of the density distribution in the r′th underground layer, where r represents the radius of the point, and θ and λ represent latitude and longitude, respectively. and The fully normalized spherical harmonic function can be expressed as:

[0083]

[0084]

[0085] in This represents the fully normalized associated Legendre function.

[0086] Taking the partial derivative of the gravitational potential with respect to the direction of r in the above formula (4), we can obtain the spherical harmonic expression for the observed gravitational anomaly Δg(r,θ,λ):

[0087]

[0088] in and The spherical harmonic coefficients of the gravitational anomaly are calculated using the following formula:

[0089]

[0090] Where i = 1, 2. Substituting formula (8) into formula (7), we obtain the spherical harmonic expression formula for gravity anomalies as follows:

[0091]

[0092] As can be seen from the above formula, for large-scale regional gravity field forward modeling or global gravity field forward modeling, the spherical harmonic domain method can be used. The specific calculation method is as follows: first, calculate the spherical harmonic coefficients of each layer of two-dimensional density distribution, and then use formula (9) to calculate the gravity anomaly on the observation surface.

[0093] S30. Convert the inversion range data into a three-dimensional gravity field forward modeling kernel matrix, and calculate the observed gravity anomaly vector based on the three-dimensional gravity field forward modeling kernel matrix.

[0094] The formula for calculating the observed gravity anomaly vector is as follows:

[0095] Kρ=Δg (10)

[0096] Where KA represents the three-dimensional gravity field forward modeling kernel matrix, which is an Mx*Ny matrix with elements K ij Let ρ represent the gravitational field response of the j-th unit cell with density of unit 1 at the i-th observation point; ρ represents the density vector, and Δg represents the observed gravitational anomaly vector. It should be noted that the three-dimensional gravitational field forward modeling kernel matrix A is a huge non-sparse matrix. Storing all elements would consume a large amount of storage space, making large-scale forward and inverse modeling calculations impossible on ordinary computers. Therefore, a method of dynamically calculating the gravitational anomaly response during inversion is adopted here to avoid the storage problem of the huge kernel matrix.

[0097] S40. Based on the inversion range data, gravity anomalies, and observed gravity anomaly vectors, construct the inversion objective function of the three-dimensional gravity field hybrid constraint.

[0098] The formula for calculating the inversion objective function Ψ(ρ) is as follows:

[0099]

[0100] Where ||Δg-K·ρ||2 This is a data fitting term designed to ensure that the retrieved three-dimensional density distribution produces a gravity anomaly similar to the observed data. ||W ρ (ρ-ρ ref )|| 2 The constraint for the smooth model aims to ensure that the inverted three-dimensional density model is smooth and continuous in all directions, conforming to most geological conditions. The gradient operator term aims to ensure that the inverted 3D density distribution can have a large rate of change in the x, y, and z directions, that is, to allow abrupt changes in the model density. ρ ref This is the reference density model.

[0101] S50. Input the pre-acquired three-dimensional seismic velocity model data and convert the three-dimensional seismic velocity model data into a reference density model.

[0102] In this embodiment of the invention, the pre-acquired seismic three-dimensional velocity model data is converted into a reference density model using the following formula:

[0103] ρ ref =1.6612V p -0.4721V p 2 +0.0671V p 3 -0.0043V p 4 +0.000106V p 5 (12)

[0104] Among them, V p The velocity is the longitudinal wave velocity, and the density unit is g / cm³. 3 The unit of velocity is km / s. If the velocity model for this region is given in terms of shear waves, then the shear wave velocity can be converted into the longitudinal wave velocity using the following formula:

[0105] V p =0.9409 + 2.0947V s -0.8206V s 2 +0.2683V s 3 -0.0251V s 4 (13)

[0106] S60. Calculate the gradient operator in the inversion objective function, and simplify the inversion objective function based on the calculated gradient operator. Then, solve the simplified inversion objective function using the alternating direction multiplier method.

[0107] The formula for calculating the gradient operator is as follows:

[0108]

[0109] Where D x D y D z Let D represent the differential operators in the x, y, and z directions, respectively. x It has the following form:

[0110]

[0111] It should be noted that the remaining D y and D z Having the same characteristics as D x Similar structures.

[0112]

[0113]

[0114] The simplified calculation process for the inversion objective function is as follows:

[0115] Substituting formulas (14) and (15) into the inversion objective function, we can simplify to obtain:

[0116] Ψ(ρ)=||Δg-K·ρ|| 2 +β|Φρ| (16)

[0117] The expression for Φ is as follows:

[0118] Φ=[I T ,(D x W) T ,(D y W) T ,(D z W) T ] T (17)

[0119] Where I represents the identity matrix, i.e., a matrix whose elements are all 1; the superscript T indicates the transpose operation of the matrix; and Φ represents the combined term of the smooth model constraint and the gradient operator constraint.

[0120] W represents the gravity field inversion depth weighting matrix, which is calculated using the following formula:

[0121]

[0122] Where, N d N represents the number of observation points. m K represents the number of subsurface three-dimensional subdivision units. i,1This represents the i-th element in the first column of the kernel matrix.

[0123] The method of solving the simplified inversion objective function using the alternating direction multiplier method includes:

[0124] The simplified inversion objective function is rewritten in the form of a Lagrange equation; that is, formula (16) is rewritten in the form of a Lagrange equation:

[0125] Ψ(ρ,χ,ν * ,β)=||Δg-K·ρ|| 2 +β|χ|+ν *T (χ-Φρ) (19)

[0126] Where χ=Φρ,ν * This represents a Lagrange vector.

[0127] To solve the above formula (19), we introduce the augmented Lagrange function:

[0128] Ψ(ρ,χ,ν * ,β,τ)=Ψ(ρ,χ,ν * ,β,τ)+0.5τ||χ-Φρ|| 2 (20)

[0129] Where τ represents any real number.

[0130] Furthermore, the problem in formula (20) is transformed into a minimization problem, that is, solving the following problem:

[0131] minimize:||Δg-K·ρ|| 2 +β|χ|+ν *T (χ-Φρ) (21)

[0132] The Lagrange vector is normalized, i.e., v = ν. * If / ρ, then the above minimization process can be further simplified to:

[0133] Ψ(ρ,χ,ν)=||Δg-K·ρ|| 2 +β|χ|+0.5τ||χ-Φρ+v|| 2 (twenty two)

[0134] Formula (22) is solved iteratively using the alternating direction multiplier method, which simultaneously satisfies the following equations:

[0135]

[0136] Suppose the solution to the above problem in the kth iteration is (ρ (k) ,χ (k) ,ν (k)), then the solution of the (k+1)-th alternating direction method of multipliers can be expressed as:

[0137]

[0138] wherein C τ =K T K+τΦ T Φ. S(x,y) is a sign function, when x>y, S=x-y; when -y<x<y, S=0; when x<y, S=x+y, wherein x and y represent two parts of the function.

[0139] S70, complete the above iterative calculation until the convergence criterion is satisfied: ||ρ (k+1) -ρ (k) || / |ρ (k) ||<ε, wherein ε is 10 -5 .

[0140] The present invention is applicable to a new multi-scale three-dimensional gravity and seismic joint method in Cartesian coordinate system and spherical coordinate system. By constructing a brand-new inversion objective function, a multi-norm constraint system based on smooth constraint norm and gradient operator focusing constraint and a gravity and seismic joint inversion system are formed, which can effectively solve the problems of low depth resolution and strong multiplicity of solutions of traditional three-dimensional inversion methods.

[0141] In the iterative solution, it is necessary to perform forward calculation on the gravity field response of the modified variable of the model. In the present invention, a frequency domain forward method based on FFT is adopted in the Cartesian coordinate system, and a spherical harmonic domain forward method is adopted in the spherical coordinate system, and the storage problem of a large Jacobian matrix is avoided. When solving for the minimization of the inversion objective function, the alternating direction method of multipliers is adopted to avoid the problem of traditional solving of large-scale linear equations.

[0142] Exemplarily, Embodiment 1:

[0143] In this embodiment, a salt dome model is constructed to verify the correctness and high depth resolution of the inversion method proposed by the present invention. The model has a size of 10 km in both the east-west direction and the north-south direction, and 4 km in the depth direction, and the height of the observation surface is 100 meters above the ground. The abnormal body is a salt dome with density decreasing with depth, and the density anomaly decreases from -1000kg / m 3 to 0 from the ground to the deep part. The specific position and shape are as shown in Figure 2-4 . The model is uniformly divided into 100, 100 and 40 segments in the depth direction, the latitude direction and the longitude direction respectively, the total number of models is 400000, and the number of observation data is 10000. The gravity anomaly generated by the model on the observation surface is calculated by forward modeling, which is used as the observation data for inversion.

[0144] The inversion results of the traditional inversion method based on L2 norm smooth constraint and the method proposed by the present invention are compared, as shown in Figure 2-7 As shown, the traditional L2 norm smooth constraint inversion method struggles to accurately determine the deep location of the salt dome, while the proposed method clearly inverts the true contour of the salt dome, and the density value of the anomaly is closer to the real model. Furthermore, since the proposed method operates in the frequency domain, its computational efficiency is nearly 100 times higher than traditional methods. Traditional inversion methods take 10 hours, while the proposed method takes only 6 minutes.

[0145] Example, Embodiment Two:

[0146] The above embodiments are synthetic inversion models. Next, the effectiveness of the proposed combined gravity and seismic inversion method will be tested. The study area is the Tibetan Plateau, with latitude ranging from 30° to 37° and longitude ranging from 82° to 95°. First, the three-dimensional seismic data of this area is mapped, as shown below. Figure 11-13 The figure shows velocity models at different depths. This velocity model is then converted into a three-dimensional density model, serving as a reference model for three-dimensional gravity inversion. Finally, the proposed hybrid-constraint gravity-seismic joint frequency domain inversion method is used to obtain the subsurface three-dimensional density distribution, as shown below. Figure 8-10 As shown. Comparison Figure 8-13 It can be observed that the results based on the combined inversion of heavy and light earthquakes, compared to the second-column seismic velocity model, contain more detailed information and improve the lateral resolution of the inversion. From Figure 14 The six density sections show that the inversion results are consistent with the geological structure information.

[0147] To demonstrate the correctness of the proposed method, a common salt dome model inversion example is provided, where the anomaly is a salt dome with density varying with depth. Comparison with traditional smooth-constrained inversion methods reveals that the proposed method significantly improves both inversion accuracy and depth resolution. Traditional methods cannot accurately invert the bottom interface of the salt dome, and the horizontal range is larger than the actual model, while the proposed method yields results close to the actual model. Next, the proposed inversion method is applied to the joint inversion of gravity and seismic data in the Tibetan Plateau. First, the seismic data is converted into an initial three-dimensional density distribution model. Then, the proposed inversion method is used to obtain the joint inversion of the three-dimensional density distribution. It can be seen that the density distribution obtained using the joint gravity and seismic inversion method has a certain correlation with the initial seismic data, proving the correctness of the method. Therefore, the proposed method is a multi-scale, high-precision three-dimensional gravity field joint inversion method suitable for both exploration and regional scales.

[0148] The method proposed in this invention has high versatility. In a rectangular coordinate system, it can be used in fields such as mineral resource exploration, oil and gas resource exploration, engineering environment exploration, and geological hazard investigation; while in a spherical coordinate system, it can be used for three-dimensional density inversion and dynamic research of the crust and mantle of the Tibetan Plateau, three-dimensional density imaging and tectonic evolution research of the shallow lithosphere of the Moon, and can also be used for resource exploration research on Mars and other celestial bodies.

[0149] It should be understood that although the above description follows a certain order, these steps are not necessarily executed in that order. Unless otherwise expressly stated herein, there is no strict order restriction on the execution of these steps, and they can be executed in other orders. Moreover, some steps in this embodiment may include multiple steps or multiple stages, which are not necessarily completed at the same time, but may be executed at different times. The execution order of these steps or stages is not necessarily sequential, but may be performed alternately or in turn with other steps or at least a portion of the steps or stages in other steps.

[0150] It should be understood that, as used herein, the singular form "a" is intended to include the plural form as well, unless the context clearly supports an exception. It should also be understood that, as used herein, "and / or" refers to any and all possible combinations of one or more of the associatedly listed items. The embodiment numbers disclosed above are for descriptive purposes only and do not represent the superiority or inferiority of the embodiments.

[0151] Those skilled in the art should understand that the discussion of any of the above embodiments is merely exemplary and is not intended to imply that the scope of the invention (including the claims) is limited to these examples. Within the framework of the invention, technical features of the above embodiments or different embodiments can be combined, and many other variations of different aspects of the invention exist, which are not provided in the details for the sake of brevity. Therefore, any omissions, modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the invention should be included within the protection scope of the invention.

Claims

1. A multi-scale three-dimensional re-seismic joint frequency domain inversion method based on hybrid constraints, characterized in that, The method includes: Discrete Fourier transform was used to determine the observed gravity anomalies on the ground, and a weighted average of values ​​at multiple Gaussian points was used to discretize the observed gravity anomalies. Based on the discretized observed gravity anomalies, the inversion range data was determined. Calculate the spherical harmonics and spherical harmonic coefficients of the two-dimensional density distribution in each layer of the inversion range, and calculate the gravity anomaly on the observation surface based on the spherical harmonics and spherical harmonic coefficients; The inversion range data is converted into a three-dimensional gravity field forward modeling kernel matrix, and the observed gravity anomaly vector is calculated based on the three-dimensional gravity field forward modeling kernel matrix. Construct the inversion objective function for a three-dimensional gravity field with hybrid constraints; The pre-acquired three-dimensional seismic velocity model data is converted into a reference density model, and the three-dimensional seismic velocity model data or reference density model is used for inversion; Calculate the gradient operator in the inversion objective function, simplify the inversion objective function based on the gradient operator, and then solve the simplified inversion objective function using the alternating direction multiplier method; Continue the above iterative calculations until the convergence criterion is met; The calculation formula for the inversion objective function is as follows: In the formula, For the inversion objective function, For data fitting terms, For smooth model constraints, For gradient operator terms, ρ ref As a reference density model, ρ This represents the density vector to be solved, and β is the weighting coefficient between the data fitting term and the model fitting term; The formula for calculating the gradient operator is as follows: in, Denotes the differential operator term in the x-direction. Denotes the differential operator term in the y-direction. Represents the differential operator term in the z-direction; The simplified formula for calculating the inversion objective function is as follows: in, This represents a combined term of smooth model constraints and gradient operator constraints.

2. The multi-scale three-dimensional re-seismic joint frequency domain inversion method based on hybrid constraints as described in claim 1, characterized in that, The expression for the observed gravity anomaly is: In the formula, This is represented as an observed gravity anomaly; Represented as a spectrum; x m express x Direction first m One point, y n express y Direction first n One point, z 0 represents the height of the observation point; k x , k y express x direction, y Wave number in direction; dk x , dk y express x direction, y The differential of the wave number in the direction; e π is a natural constant with a value of 2.71828; π is the ratio of a circle's diameter to its diameter.

3. The multi-scale three-dimensional re-seismic joint frequency domain inversion method based on hybrid constraints as described in claim 2, characterized in that, The expression for the inversion range data is: In the formula, e represents the natural exponent, and i represents the imaginary unit. and Let Δk represent the k-th and j-th Gaussian point values ​​in the north-south and east-west directions, respectively. x Δk y Indicates discrete spectral interval; k xp and k yq This represents the p-th and q-th discrete wavenumbers in the x and y directions, respectively. x m and y n Let N represent the position of the m-th point in the x-direction and the position of the n-th point in the y-direction, respectively; x N y M x M y These represent the upper bound of the discrete spectral wavenumber in the x-direction, the upper bound of the discrete spectral wavenumber in the y-direction, the lower bound of the discrete spectral wavenumber in the x-direction, and the lower bound of the discrete spectral wavenumber in the y-direction, respectively.

4. The multi-scale three-dimensional re-seismic joint frequency domain inversion method based on hybrid constraints as described in claim 3, characterized in that, The calculation of the spherical harmonic coefficients of the two-dimensional density distribution in each layer of the inversion range, and the calculation of gravity anomalies on the observation surface based on the spherical harmonic coefficients, includes: Calculate the gravitational potential, and based on the gravitational potential, calculate the spherical harmonics and spherical harmonic coefficients of the observed gravitational anomaly; Gravity anomalies on the observation surface are calculated based on spherical harmonics and spherical harmonic coefficients.

5. The multi-scale three-dimensional re-seismic joint frequency domain inversion method based on hybrid constraints as described in claim 4, characterized in that, The gravitational potential can be expressed by the following formula: In the formula, V(r,θ,φ) represents the gravitational potential. Represents the gravitational constant, π represents pi, a and b represent the spherical harmonic order, and r The radius representing the location of the underground anomaly. and Indicates the rth underground The density distribution spherical harmonic coefficient of the layer, r represents the radius of the calculation point, and θ and λ represent latitude and longitude, respectively. and Let represent the cosine and sine components of the fully normalized spherical harmonic function, respectively, and let n represent a positive integer.

6. The multi-scale three-dimensional re-seismic joint frequency domain inversion method based on hybrid constraints as described in claim 5, characterized in that, The formula for calculating the observed gravity anomaly vector is as follows: Kρ =Δg In the formula, K Let K denote the forward modeling kernel matrix of the three-dimensional gravity field, which is an Mx*Ny matrix with elements K. ij This represents the gravitational field response of the j-th unit cell to the i-th observation point. ρ Let denot be the density vector to be solved, and Δg denote the observed gravity anomaly vector.

7. The multi-scale three-dimensional re-seismic joint frequency domain inversion method based on hybrid constraints as described in claim 1, characterized in that, The The expression is as follows: In the formula, I Represents the identity matrix. W This represents the depth-weighted matrix for gravity field inversion, with superscript indicating the depth. T This represents the transpose operation on a matrix. D x , D y , D z Let x, y, and z represent the differential operators in the x, y, and z directions, respectively.

Images