Multi-scale three-dimensional gravitational-seismic joint frequency domain inversion method based on hybrid constraint

Abstract

The invention relates to the technical field of geophysical exploration, in particular to a multi-scale three-dimensional gravitational-seismic joint frequency domain inversion method based on hybrid constraints. The method comprises the following steps: determining observation gravity anomaly on the ground by adopting discrete Fourier transform, and discretizing the observation gravity anomaly by adopting a method of weighted average of values at a plurality of Gaussian points; calculating spherical harmonics and spherical harmonic coefficients of two-dimensional density distribution of each layer in the inversion range, and calculating gravity anomalies on the observation surface based on the spherical harmonics and the spherical harmonic coefficients; converting the inversion range data into a three-dimensional gravitational field forward kernel matrix; constructing a three-dimensional gravitational field hybrid constraint inversion objective function; inputting pre-collected earthquake three-dimensional velocity model data, and converting the earthquake three-dimensional velocity model data into a reference density model; and calculating a gradient operator in the inversion objective function, simplifying the inversion objective function based on the gradient operator obtained through calculation, and then solving the simplified inversion objective function by using an alternating direction multiplier method.

Description

Technical Field

[0001] The present invention relates to the technical field of geological exploration, and in particular to a multi-scale three-dimensional gravity-seismic joint frequency-domain inversion method based on hybrid constraints in geophysical exploration technology. Background Art

[0002] As one of the geophysical exploration methods, gravity exploration has the advantages of simple method, high economy and large detection depth, and is widely used in various scale geological and geophysical problems. For example, exploration of mineral and oil and gas resources at the local scale, development of geothermal resources, etc., and research on the density structure of the crust-mantle of the Earth's lithosphere at the regional scale, research on the three-dimensional density structure and formation and evolution history of lunar mass anomaly basins, etc. Gravity field inversion is one of the important links in gravity exploration, that is, based on the two-dimensional gravity anomaly data obtained by observation, the three-dimensional density distribution underground is inversed.

[0003] However, it is obvious that the method of inversing three-dimensional density distribution with two-dimensional data belongs to an underdetermined problem in mathematics, that is, the number of unknowns is more than the number of equations, so the solutions to such problems are not unique. In addition, due to the serious volume effect of the gravity field, that is, according to the law of universal gravitation, the observed gravity anomaly is affected by both the mass and distance of the anomaly body, and the mass is related to density and volume. Therefore, a sphere with a shallow burial depth but a large mass and a sphere with a deep burial depth but a small mass may produce the same gravity field at a point on the ground. The emergence of the volume effect further exacerbates the non-uniqueness of the inversion, resulting in two major problems in three-dimensional gravity inversion: serious multi-solution and low depth resolution.

[0004] To solve the problems of strong non-uniqueness and low depth resolution in gravity inversion, relevant technical personnel in this field have made various improvements. These improved methods can basically be divided into two ideas: The first improvement idea is to vigorously explore the potential of the observed gravity data itself, adding constraints from multiple dimensions to reduce non-uniqueness. Such improvements include: various forms of depth weighting functions, various forms of physical property upper and lower limit constraints, various forms of norm systems, various forms of self-constrained mathematical methods, etc.; The second improvement idea is the joint inversion of multi-source data, and the most common ones are the joint inversion of gravity and seismic data, the joint inversion of gravity and magnetic data, and the joint inversion of gravity gradient tensors. Since the seismic wave velocity has a good correlation with density to a certain extent, the joint inversion of gravity and seismic data has a high physical mechanism as a guarantee. The joint inversion of gravity gradient tensors is even simpler, all based on the basic principle of gravity, just the mutual complement of different components, and essentially no redundant physical property information is added, so its joint effect is not as good as that of gravity and seismic data. However, the above-mentioned improvement measures only reduce the existence of non-uniqueness to a certain extent and improve the depth resolution, and the problem of non-uniqueness in gravity inversion still needs to be further developed. On the other hand, restricted by problems such as the huge storage of the Jacobian matrix of gravity field inversion and the serious time-consuming of forward calculation, it is difficult for traditional gravity field inversion to solve the problem of large-scale three-dimensional density inversion. In addition, most of the above methods are proposed based on the rectangular coordinate system at the exploration scale and are difficult to meet the needs of large-scale gravity field inversion in the spherical coordinate system. Summary of the Invention

[0005] To solve the technical problems existing in the above-mentioned prior art, the present invention provides a multi-scale three-dimensional gravity-seismic joint frequency domain inversion method based on hybrid constraints, which has solved the problems of low depth resolution, strong non-uniqueness, and low adaptability in traditional three-dimensional gravity field inversion under the current rectangular coordinate system, and has solved the problems proposed in the background technology.

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

[0007] In the first aspect, in an embodiment provided by the present invention, a multi-scale three-dimensional gravity-seismic joint frequency domain inversion method based on hybrid constraints is provided, and the method includes the following steps:

[0008] The discrete Fourier transform is used to determine the observed gravity anomaly on the ground, and the method of weighted averaging of values at multiple Gaussian points is used to discretize the observed gravity anomaly, and the inversion range data is determined based on the discretized observed gravity anomaly;

[0009] Calculate the spherical harmonics and 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 harmonics and spherical harmonic coefficients;

[0010] Convert the inversion range data into a 3D gravity forward kernel matrix, and calculate the observed gravity anomaly vector based on the 3D gravity forward kernel matrix;

[0011] Construct an inversion objective function with mixed constraints for the 3D gravity field;

[0012] Input the pre-collected 3D seismic velocity model data and convert it into a reference density model;

[0013] Calculate the gradient operator in the inversion objective function, simplify the inversion objective function based on the calculated gradient operator, and then solve the simplified inversion objective function using the alternating direction method of multipliers;

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

[0015] As a further solution of the present invention, the expression of the observed gravity anomaly is:

[0016]

[0017] In the formula, Δg(x m ,y n ,z0) represents the observed gravity anomaly; represents the spectrum; xm represents the m-th point in the x direction, y n represents the n-th point in the y direction, z0 represents the height of the observation point; k x 、k y represent the wave numbers in the x direction and y direction; dk x 、dk y represent the differentials of the wave numbers in the x direction and y direction; e is the natural constant, and its value is 2.71828; π is the pi.

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

[0019]

[0020] In the formula, e represents the natural exponent, i represents the imaginary unit, and λ j represent the k-th and j-th Gaussian point values in the north-south and east-west directions, Δk x 、Δk y represent the discrete spectrum intervals; k xp and k yq represent the p-th and q-th discrete wave numbers in the x direction and y direction; Mx and Ny represent the number of segments divided in the x direction and y direction respectively; x m and y n 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 solution of the present invention, calculating the spherical harmonic coefficients of the two-dimensional density distribution of each layer within the calculation and inversion range, and calculating the gravity anomaly on the observation surface based on the spherical harmonic coefficients, includes:

[0022] Calculating the gravitational potential, and obtaining the spherical harmonic and spherical harmonic coefficients for observing the gravity anomaly based on the gravitational potential;

[0023] Calculating the gravity anomaly on the observation surface based on the spherical harmonic and spherical harmonic coefficients.

[0024] As a further solution of the present invention, the gravitational potential can be expressed by the following formula:

[0025]

[0026] In the formula, V(r,θ,φ) represents the gravitational potential, γ represents the gravitational constant, π represents the circumference ratio, a and b represent the spherical harmonic orders, r′ represents the radius of the location of the underground anomaly, and represents the spherical harmonic coefficients of the density distribution of the r′-th layer underground, r represents the radius of this point, and θ and λ represent the latitude and longitude respectively.

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

[0028] Kρ = Δg

[0029] In the formula, K represents the forward modeling kernel matrix of the three-dimensional gravity field, which is a matrix of Mx*Ny, and its element K ij represents the gravity field response generated by the j-th unit body with a density of 1 on the i-th observation point, ρ represents the density vector, and Δg represents the observed gravity anomaly vector.

[0030] As a further solution 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 is the data fitting term, ||W ρ (ρ - ρ ref )|| 2 is the smooth model constraint, is the gradient operator term, and ρ ref is the reference density model.

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

[0034]

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

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

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

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

[0039] As a further solution of the present invention, the expression of Φ 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 represents the transpose operation of the matrix.

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

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

[0044] These aspects or other aspects of the present invention will be more clearly understood in the following description of the embodiments. It should be understood that the above general description and the following detailed description are only exemplary and explanatory, and cannot limit the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the following drawings are only some embodiments of the present invention. For those of ordinary skill in the art, other embodiments can be obtained based on these drawings without creative efforts.

[0046] Figure 1 ​Flowchart of the multi-scale three-dimensional seismic and gravity joint frequency-domain inversion method based on hybrid constraints for an embodiment of the present invention.

[0047] Figure 2 True model diagram in the contrast diagram of the density inversion section of the three-dimensional salt dome model.

[0048] Figure 3 Result diagram obtained by the traditional L2 norm inversion based on smooth constraints in the contrast diagram of the density inversion section of the three-dimensional salt dome model.

[0049] Figure 4 Result diagram obtained by the traditional inversion method based on hybrid constraints in the contrast diagram of the density inversion section of the three-dimensional salt dome model.

[0050] Figure 5 True three-dimensional density distribution diagram of the salt dome in the three-dimensional stereoscopic display of the three-dimensional salt dome density inversion comparison.

[0051] Figure 6 Three-dimensional stereoscopic display diagram of the traditional smooth constraint inversion in the three-dimensional stereoscopic display of the three-dimensional salt dome density inversion comparison.

[0052] Figure 7 Three-dimensional stereoscopic display diagram of the inversion by the method proposed in the present invention in the three-dimensional stereoscopic display of the three-dimensional salt dome density inversion comparison.

[0053] Figure 8 Depth slice diagram of the seismic and gravity joint three-dimensional density inversion at a depth of 10 km in the Qinghai-Tibet Plateau using the method proposed in the present invention.

[0054] Figure 9 Depth slice diagram of the seismic and gravity joint three-dimensional density inversion at a depth of 40 km in the Qinghai-Tibet Plateau using the method proposed in the present invention.

[0055] Figure 10 Depth slice diagram of the seismic and gravity joint three-dimensional density inversion at a depth of 140 km in the Qinghai-Tibet Plateau using the method proposed in the present invention.

[0056] Figure 11 Depth slice diagram of the three-dimensional seismic velocity model at a depth of 10 km in the Qinghai-Tibet Plateau. This velocity model is converted into an initial density model and then participates in the seismic and gravity joint inversion.

[0057] Figure 12 Depth slice diagram of the three-dimensional seismic velocity model at a depth of 40 km in the Qinghai-Tibet Plateau. This velocity model is converted into an initial density model and then participates in the seismic and gravity joint inversion.

[0058] Figure 13 Depth slice diagram of the three-dimensional seismic velocity model at a depth of 140 km in the Qinghai-Tibet Plateau. This velocity model is converted into an initial density model and then participates in the seismic and gravity joint inversion.

[0059] Figure 14 It is the depth section map of different profiles in the three-dimensional combined gravity and seismic inversion in the Qinghai-Tibet Plateau. The first row is the horizontal position map of different profiles, namely the CC′, DD′, EE′, FF′, GG′, and HH′ sections; the second to fourth rows are the density-depth section maps of these 6 sections respectively. The dashed line in the figure represents the Moho depth of the section. Specific implementation manners

[0060] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0061] The flowchart shown in the accompanying drawings is only an example, and does not necessarily include all the contents and operations / steps, nor does it necessarily need to be executed in the described order. For example, some operations / steps can also be decomposed, combined, or partially merged, so the actual execution order may be changed according to the actual situation.

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

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

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

[0065] S10. The discrete Fourier transform is used to determine the observed gravity anomaly on the ground, and the method of weighted average of values at multiple Gauss points is used to discretize the observed gravity anomaly, and the inversion range data is determined based on the discretized observed gravity anomaly.

[0066] In the embodiment of the present invention, the discrete Fourier transform is used to determine the observed gravity anomaly on the ground, and the method of weighted average of values at multiple Gauss points is used to discretize the observed gravity anomaly, and the inversion range data is determined based on the discretized observed gravity anomaly, including:

[0067] Let the north-south range be 0–X, the east-west range be 0–Y, and the depth range be 0–Z. The north-south direction, east-west direction, and depth direction are each divided into Nx, Ny, and Nz segments, with the division intervals being Δx, Δy, and Δz respectively. Assume that the position index of each segment is p = 1, 2, …, Nx, q = 1, 2, …, Ny, and l = 1, 2, …, Nz. Assume that the density value of each block after subdivision underground is constant, which is ρ(ξ, η, ζ), where (ξ, η, ζ) is the coordinate position of the subdivided block. The coordinate position of the observation point is (x, y, z0), and the index in the horizontal direction is m = 1, 2, …, Nx, n = 1, 2, …, Ny. Then, for the forward calculation of the small-scale gravity field in the rectangular coordinate system, the influence of the earth's curvature can be ignored. Therefore, the observed gravity anomaly Δg on the ground can be expressed as:

[0068]

[0069] where, Δg(x m ,y n ,z0) represents the observed gravity anomaly; represents the spectrum; x m represents the m-th point in the x direction, y n represents the n-th point in the y direction, z0 represents the height of the observation point; k x , k y represent the wave numbers in the x direction and y direction; dk x , dk y represent the differentials of the wave numbers in the x direction and y direction; e is the natural constant, and its value is 2.71828; π is the pi. The above formula 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 number of Gauss points in the north-south and east-west directions, and k and j represent their position indices respectively. e represents the natural exponent, i represents the imaginary unit, and λ j represent the values of the k-th and j-th Gauss points in the north-south and east-west directions, Δk x , Δk y represent the discrete spectrum intervals; k xp and k yq represent the p-th and q-th discrete wave numbers in the x direction and y direction; Mx and Ny represent the number of segments divided in the x direction and y direction respectively; x m and y n represent the positions of the m-th point in the x direction and the n-th point in the y direction respectively. A kj is:

[0072]

[0073] The above formulas (1)-(3) achieve the frequency-domain forward calculation of the three-dimensional gravity field. Due to the adoption of the discrete Fourier transform and the method of weighted averaging of the values at multiple Gaussian points, compared with the traditional spatial-domain forward method, both the calculation efficiency and the calculation accuracy are improved.

[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 the 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 the gravity anomaly on the observation surface based on the spherical harmonic coefficients include:

[0076] Calculate the gravitational potential, and obtain the spherical harmonic and spherical harmonic coefficients for observing the gravity anomaly based on the gravitational potential;

[0077] Calculate the gravity anomaly on the observation surface based on the spherical harmonic and spherical harmonic coefficients.

[0078] It should be noted that before calculating the gravitational potential, it is also necessary to perform a meshing step in the rectangular coordinate system as in S10, and the meshing step is the same as that in step S10, but the observation point coordinates are (r, θ, λ), where r represents the radius of the point, and θ and λ represent the 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 the circumference ratio, a and b represent the spherical harmonic orders, r′ represents the radius of the location of the underground anomaly, and represent the spherical harmonic coefficients of the density distribution of the r′-th layer underground, r represents the radius of the point, and θ and λ represent the latitude and longitude respectively. and represent the fully normalized spherical harmonic function, which can be expressed as:

[0083]

[0084]

[0085] Where represents the fully normalized associated Legendre function.

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

[0087]

[0088] Where and represent the spherical harmonic coefficients of the gravity anomaly, and their calculation formulas are as follows:

[0089]

[0090] where i = 1, 2. Substituting formula (8) into formula (7), the spherical harmonic expression formula for the gravity anomaly can be obtained as:

[0091]

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

[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 calculation formula for 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 a matrix of Mx * Ny, and its element K ij represents the gravity field response generated by the j-th unit cell with a density of 1 on the i-th observation point; ρ represents the density vector, and Δg represents the observed gravity anomaly vector. It should be noted that the three-dimensional gravity field forward modeling kernel matrix A is a huge non-sparse matrix. If all elements are stored, it will consume a large amount of storage space and cannot perform large-scale forward and inverse calculations on ordinary computers. Therefore, here a method of dynamically calculating the gravity anomaly response during inversion is adopted to avoid the problem of storing the huge kernel matrix.

[0097] S40. Based on the inversion range data, gravity anomaly, and observed gravity anomaly vector, construct an inversion objective function with three-dimensional gravity field hybrid constraints.

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

[0099]

[0100] where ||Δg - K·ρ||2 is a data fitting term, aiming to ensure that the three-dimensional density distribution obtained by inversion can generate a gravity anomaly similar to the observed data. ||W ρ (ρ - ρ ref )|| 2 is a smooth model constraint, aiming to ensure that the three-dimensional density model obtained by inversion is smooth and continuous in all directions, conforming to most geological conditions. is a gradient operator term, aiming to ensure that the three-dimensional density distribution obtained by inversion can have a large rate of change in the x, y, and z directions, that is, allowing sudden changes in the model density. ρ ref is a reference density model.

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

[0102] In the embodiment of the present invention, for the input of the pre-collected seismic three-dimensional velocity model data and its conversion into a reference density model, the following formula is used:

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

[0104] where V p is the P-wave velocity, the density unit is g / cm 3 , and the velocity unit is km / s. If the velocity model in this area is given in terms of S-wave, then the S-wave velocity can be converted into the P-wave velocity through 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, simplify the inversion objective function based on the calculated gradient operator, and then use the alternating direction multiplier method to solve the simplified inversion objective function.

[0107] The calculation formula of the gradient operator is as follows:

[0108]

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

[0110]

[0111] It should be noted that the remaining D y and D z have a structure similar to that of D x .

[0112]

[0113]

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

[0115] Substituting formulas (14) and (15) into the inversion objective function and after simplification, we can get:

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

[0117] where the expression of Φ 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, that is, the matrix with all matrix elements being 1; the superscript T represents the transpose operation of the matrix, and Φ represents the comprehensive term of the smooth model constraint and the gradient operator term constraint.

[0120] W represents the gravity field inversion depth weighting matrix, and has the following calculation formula:

[0121]

[0122] where, N d represents the number of observation points, and N m represents the number of three-dimensional underground subdivision unit bodies. K i,1 ​Denotes the i-th element in the first column of the kernel matrix elements.

[0123] The use of the alternating direction method of multipliers to solve the simplified inversion objective function includes:

[0124] Rewrite the simplified inversion objective function into the form of a Lagrangian equation; that is: rewrite formula (16) into the form of a Lagrangian equation:

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

[0126] where χ = Φρ, ν * Denotes the Lagrangian vector.

[0127] To solve the above formula (19), an augmented Lagrangian function is introduced:

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

[0129] where τ represents an arbitrary real number.

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

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

[0132] Normalize the Lagrangian vector, that is v = ν * / ρ, then the above minimization process can be further simplified to:

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

[0134] Use the alternating direction method of multipliers to iteratively solve formula (22), that is, simultaneously satisfy the following equations:

[0135]

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

[0137]

[0138] where 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, where 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) || < ε, where ε is 10 -5 .

[0140] The present invention is applicable to a multi-scale three-dimensional gravity field gravity and seismic joint new method in a rectangular coordinate system and a spherical coordinate system. By constructing a new inversion objective function, a multiple norm constraint system based on a smooth constraint norm and a 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 multi-solution of traditional three-dimensional inversion methods.

[0141] When iteratively solving, it is necessary to forward calculate the gravity field response of the variable of the model. The present invention adopts a frequency-domain forward method based on FFT in a rectangular coordinate system and a spherical harmonic domain forward method in a spherical coordinate system, and avoids the problem of storing a large Jacobian matrix. When minimizing the inversion objective function, the alternating direction method of multipliers is adopted to avoid the problem of solving a large linear equation system traditionally.

[0142] Exemplarily, Example 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 is 10 km in the east-west direction and 10 km in the north-south direction, 4 km in the depth direction, and the height of the observation surface is 100 meters from the ground. The anomaly body is a salt dome with a density decreasing with depth, and the density anomaly decreases from -1000 kg / m 3 to 0 from the ground to the deep part. The specific position and shape are as Figures 2-4 shown. The model is uniformly divided into 100, 100, and 40 segments in the depth direction, latitude direction, and longitude direction respectively, and the total number of model elements is 400000, and the number of observation data is 10000. Forward calculate the gravity anomaly generated by the model on the observation surface, and use it as the observation data for inversion.

[0144] Compare the inversion results of the traditional L2-norm smooth constraint inversion method and the method proposed by the present invention, asFigures 2-7 As shown in the figure. It can be seen that it is very difficult to accurately obtain the deep position of the salt dome from the three-dimensional density distribution obtained by the traditional smooth constraint inversion method using the L2 norm. However, the method proposed in the present invention can clearly invert the true contour of the salt dome, and the density value of the anomaly body is closer to the true model. In addition, since the method proposed in the present invention is carried out in the frequency domain, the calculation efficiency is increased by nearly 100 times compared with the traditional method. The traditional method takes 10 hours for inversion, while the inversion method of the present invention takes 6 minutes.

[0145] Exemplarily, Example 2:

[0146] The above embodiment is a synthetic inversion model. Next, the effectiveness of the joint gravity-seismic inversion method proposed in the present invention is tested. The research area is the Qinghai-Tibet Plateau, with a latitude range of 30° to 37° and a longitude range of 82° to 95°. First, the three-dimensional seismic data of this area is mapped, as Figures 11-13 shown, which is the velocity model at different depths. The velocity model is converted into a three-dimensional density model as the reference model for three-dimensional gravity inversion. Then, the joint gravity-seismic frequency domain inversion method based on hybrid constraints proposed in the present invention is used to obtain the underground three-dimensional density distribution, as Figures 8-10 shown. By comparing Figures 8-13 , it can be found that the result of the joint gravity-seismic inversion has more detailed information and higher lateral resolution compared with the seismic velocity model in the second column. From Figure 14 the six density sections, it can be seen that the inversion result conforms to the geological structure information.

[0147] To prove the correctness of the method proposed in the present invention, a common salt dome model inversion example is given in the present invention, and the anomaly body is a salt dome with a density varying with depth. By comparing with the traditional smooth constraint inversion method, it can be found that the method proposed in the present invention has obvious improvements in both inversion accuracy and depth resolution. It can be seen that the traditional method cannot accurately invert the bottom interface of the salt dome, and the horizontal range is larger than the true model, while the inversion result of the method proposed in this paper is close to the true model. Next, the applied inversion method is used for the joint gravity-seismic inversion of the Qinghai-Tibet Plateau. First, the seismic data is converted into an initial model of three-dimensional density distribution, and then the inversion method proposed in the present invention is used to obtain the three-dimensional density distribution after joint inversion. It can be seen that the density distribution obtained by the joint gravity-seismic inversion method has a certain correlation with the initial seismic data, which proves the correctness of the method. Therefore, the method proposed in the present invention is a multi-scale high-precision three-dimensional gravity field joint inversion method applicable to exploration scale and regional scale.

[0148] The method proposed by the present invention has high universality. In the rectangular coordinate system, it can be used in the fields of mineral resource exploration, oil and gas resource exploration, engineering environment exploration, geological disaster detection, etc.; while in the spherical coordinate system, it can be used for three-dimensional density inversion and dynamic research of the crust and mantle in the Qinghai-Tibet Plateau, three-dimensional density imaging and tectonic evolution research of the shallow part of the lunar lithosphere, and can also be used for resource exploration research on other celestial bodies such as Mars.

[0149] It should be understood that although the above is described in a certain order, these steps are not necessarily executed in the above order successively. Unless there is a clear indication in this article, the execution of these steps has no strict order limit, and these steps can be executed in other orders. Moreover, a part of the steps of this embodiment may include multiple steps or multiple stages, and these steps or stages are not necessarily executed at the same moment, but can be executed at different moments, and the execution order of these steps or stages is not necessarily sequential, but can be executed alternately or in turn with at least a part of other steps or steps or stages in other steps.

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

[0151] Those of ordinary skill in the art should understand that: the discussion of any of the above embodiments is only exemplary and is not intended to imply that the scope (including the claims) of the disclosure of the embodiments of the present invention is limited to these examples; under the idea of the embodiments of the present invention, the technical features in the above embodiments or different embodiments can also be combined, and there are many other variations in different aspects of the embodiments of the present invention as above, and they are not provided in detail for the sake of brevity. Therefore, any omission, modification, equivalent replacement, improvement, etc. made within the spirit and principle of the embodiments of the present invention shall be included in the protection scope of the embodiments of the present invention.

Claims

1. A multi-scale three-dimensional seismic joint frequency-domain inversion method based on hybrid constraints, characterized in that The method includes: using the discrete Fourier transform to determine the observed gravity anomaly on the ground, and using the method of weighted averaging of values at multiple Gaussian points to discretize the observed gravity anomaly, and determining the inversion range data based on the discretized observed gravity anomaly; calculating the spherical harmonics and spherical harmonic 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 harmonic coefficients; converting the inversion range data into a three-dimensional gravity field forward kernel matrix, and calculating the observed gravity anomaly vector based on the three-dimensional gravity field forward kernel matrix; constructing an inversion objective function with a three-dimensional gravity field hybrid constraint; converting the pre-collected seismic three-dimensional velocity model data into a reference density model; calculating the gradient operator in the inversion objective function, and simplifying the inversion objective function based on the gradient operator, and then using the alternating direction multiplier method to solve the simplified inversion objective function; completing the above iterative calculation until the convergence criterion is satisfied.

2. The multi-scale three-dimensional seismic joint frequency-domain inversion method based on hybrid constraints according to claim 1, wherein The expression of the observed gravity anomaly is: where, Δg(x m , y n , z0) represents the observed gravity anomaly; represents the spectrum; xm represents the m-th point in the x direction, y n represents the n-th point in the y direction, z0 represents the height of the observation point; k x , k y represent the wave numbers in the x direction and y direction; dk x , dk y represent the differentials of the wave numbers in the x direction and y direction; e is the natural constant, with a value of 2.71828; π is the pi.

3. The multi-scale three-dimensional seismic joint frequency-domain inversion method based on hybrid constraints according to claim 2, wherein, The expression of the inversion range data is: where e represents the natural exponential, and i represents the imaginary unit, and λj represent the values of the k-th and j-th Gauss points in the north-south and east-west directions, respectively. Δk x , Δk y represent the discrete spectral intervals; k xp and k yq represent the p-th and q-th discrete wave numbers in the x- and y-directions, respectively; Mx and Ny represent the number of segments divided in the x- and y-directions, respectively; x m and y n represent the positions of the m-th point in the x-direction and the n-th point in the y-direction, respectively.

4. The multi-scale three-dimensional joint frequency-domain inversion method based on hybrid constraints according to claim 3, wherein, The calculating the spherical harmonic 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 harmonic coefficients includes: calculating the gravitational potential, and calculating the spherical harmonics and spherical harmonic coefficients of the observed gravity anomaly based on the gravitational potential; calculating the gravity anomaly on the observation surface based on the spherical harmonics and spherical harmonic coefficients.

5. The multi-scale three-dimensional seismic joint frequency-domain inversion method based on hybrid constraints according to claim 4, wherein The gravitational potential can be expressed as the following formula: In the formula, V(r,θ,φ) represents the gravity potential, γ represents the gravitational constant, π represents the pi, a and b represent the spherical harmonic degrees, r′ represents the radius of the location of the subsurface anomaly, and represents the spherical harmonic coefficient of the density distribution of the r′-th subsurface layer, r represents the radius of the point, and θ and λ represent the latitude and longitude respectively.

6. The multi-scale three-dimensional seismic joint frequency-domain inversion method based on hybrid constraints according to claim 5, characterized in that The calculation formula of the observed gravity anomaly vector is as follows: Kρ = Δg In the formula, K represents the forward modeling kernel matrix of the three-dimensional gravity field, which is an Mx*Ny matrix, and its element K ij represents the gravity field response generated by the j-th unit body with a density of 1 at the i-th observation point, ρ represents the density vector, and Δg represents the observed gravity anomaly vector.

7. The multi-scale three-dimensional seismic joint frequency-domain inversion method based on hybrid constraints according to claim 6, characterized in that, The calculation formula of the inversion objective function is as follows: In the formula, Ψ(ρ) is the inversion objective function, ||Δg - K·ρ|| 2 is the data fitting term, ||W ρ (ρ - ρ ref )|| 2 is the smooth model constraint, is the gradient operator term, ρ ref is the reference density model; β is the weight coefficient between the data fitting term and the model fitting term.

8. The multi-scale three-dimensional seismic joint frequency-domain inversion method based on hybrid constraints according to claim 7, wherein, The calculation formula of the gradient operator is as follows: Among them, D x , D y , D z respectively represent the differential operators in the x, y, and z directions.

9. The multi-scale three-dimensional seismic joint frequency-domain inversion method based on hybrid constraints according to claim 8, wherein, The calculation formula of the simplified inversion objective function is as follows: Ψ(ρ) = ||Δg - K·ρ|| 2 + β|Φρ| where, Φ represents the comprehensive term of the smooth model constraint and the gradient operator term constraint.

10. The multi-scale three-dimensional seismic joint frequency-domain inversion method based on hybrid constraints as claimed in claim 9, wherein The expression of Φ is as follows: Φ = [I T ,(D x W) T ,(D y W) T ,(D z W) T T ​ In the formula, I represents the identity matrix, W represents the gravity field inversion depth weighting matrix, and the superscript T represents the transpose operation of the matrix.

Images