Joint imaging method and system based on seismic wave full waveform and gravity data

Abstract

The invention discloses a joint imaging method and system based on seismic wave full waveform and gravity data, and belongs to the technical field of geological exploration. The method comprises the following steps: firstly, carrying out preliminary modeling on a research area to obtain prior information, then collecting data to construct a total energy loss function, constructing a density and wave velocity model based on a level set method, associating inversion parameters corresponding to different physical data through an interface structure, and fusing prior information to initialize model parameters; and calculating prediction data, determining a data proportion coefficient, judging energy loss, calculating an update gradient, and updating parameters by using an asynchronous iteration Adam algorithm to finish joint inversion. According to the method, multi-physical data are effectively fused, the data proportion is reasonably set, the imaging effect is improved by means of structural similarity, prior information is fully utilized, regularization terms are optimized, good expandability is achieved, the complex interface capturing capacity is high, the geologic structure imaging resolution can be improved, interference can be reduced, and the calculation cost can be reduced. And the precision, the reliability and the overall quality of an inversion result are improved.

Description

Technical Field

[0001] The present invention relates to the field of geological exploration technology, and more particularly to a combined imaging method and system based on seismic wave full waveform and gravity data. Background Art

[0002] In the in-depth exploration of geological exploration and resource development, accurate depiction of underground geological structure is the core element to improve exploration efficiency and accuracy. Gravity inversion technology, as a passive exploration method, uses surface gravity anomaly data to infer underground density distribution, and plays an important role in regional geological surveys. However, the accuracy of this technology is easily affected by many factors, such as the noise of observation data, the non-uniqueness of the model, and the uncertainty of density distribution, which has certain limitations in high-resolution applications.

[0003] Relatively speaking, full waveform inversion technology simulates the propagation characteristics of seismic waves and inverts the wave velocity model of underground geological structures, which theoretically has the advantages of high resolution and high precision. However, in practical applications, full waveform inversion faces many challenges, especially when dealing with complex geological structures, such as salt body models, where the scattering and reflection of seismic waves are significant, often causing the inversion results to deviate from the actual situation and making it difficult to meet the needs of high-precision exploration. In response to the limitations of full waveform inversion, researchers have explored improved technologies such as multi-scale inversion and adaptive regularization to enhance the resolution of complex geological structures. Although some progress has been made, there is still a lot of room for improvement in full waveform inversion technology, and more in-depth exploration and innovation are still needed to further improve the efficiency and accuracy of inversion.

[0004] In this context, joint inversion methods have gradually become an important way to solve the above problems. Among them, the joint inversion technology of gravity and seismic wave data uses the complementarity of the two data to significantly improve the imaging capabilities of deep areas and complex geological structures. The key to joint inversion is to establish the connection between different physical parameters so that the two data can complement each other. However, the existing joint inversion of gravity and seismic waves mainly links the wave velocity and density models through empirical formulas of physical properties, and cannot effectively utilize the structural similarity of the target area in the wave velocity and density models during the imaging process.

[0005] Therefore, how to propose a joint imaging method and system based on seismic wave full waveform and gravity data to improve the efficiency and accuracy of gravity-full waveform joint inversion is a problem that technical personnel in this field urgently need to solve. Summary of the invention

[0006] In view of this, the present invention provides a joint imaging method and system based on seismic wave full waveform and gravity data to improve the efficiency and accuracy of gravity-full waveform joint inversion.

[0007] In order to achieve the above object, the present invention adopts the following technical solution:

[0008] In one aspect, the present invention discloses a joint imaging method based on seismic wave full waveform and gravity data, comprising the following steps:

[0009] S1. Acquire prior information of the target area, gravity anomaly measurement data and seismic wave full waveform measurement data;

[0010] S2. Construct density model and wave velocity model by combining level set function, and connect the inversion parameters corresponding to different physical data through interface structure;

[0011] S3. Initializing the level set function, and initializing the parameters of the density model and the wave velocity model based on the prior information;

[0012] S4. Calculating predicted gravity anomaly data and predicted seismic wave full waveform data according to the parameters of the density model and the wave velocity model;

[0013] S5. construct an energy loss function using the gravity anomaly measurement data, the seismic wave full waveform measurement data, the predicted gravity anomaly data and the predicted seismic wave full waveform data; if the data fitting term in the energy loss function is less than the target value, execute S9, otherwise execute S6;

[0014] S6. Calculate the updated gradient;

[0015] S7. Using an asynchronous iterative Adam algorithm, updating the parameters of the level set function, the density model and the wave velocity model based on the updated gradient;

[0016] S8. Return to S4;

[0017] S9. Output the updated level set function, density model, and wave velocity model.

[0018] Preferably, S2 includes:

[0019] Define the level set function φ, and characterize the common interface structure through the zero level set {r|φ(r)=0};

[0020] The density model ρ and wave velocity model v are constructed by combining the level set function as follows:

[0021] ρ=H τ (φ)×ρ0;

[0022] v=v1×H τ (φ)+v2×(1-H τ (φ));

[0023] Where ρ0 is the relative density of the target medium (such as salt body) to the background area, v1 is the wave velocity of the target medium, v2 is the wave velocity outside the target medium area, and H τ (·) is the smooth approximation function of the Heaviside function.

[0024] Preferably, S3 includes:

[0025] The level set function φ is initialized as a signed distance function of an ellipse to the center of the target area;

[0026] The parameters of the density model and the wave velocity model are initialized based on prior information, including the relative density ρ0 of the target medium to the background area, the wave velocity v1 of the target medium, and the wave velocity v2 outside the target medium area.

[0027] Preferably, the predicted gravity anomaly data in S4 The calculation formula is as follows:

[0028]

[0029] in is the component of the gravity kernel function in the height direction, They are The height direction component of n represents the spatial dimension considered; is the position coordinate of the measurement point, r is the position coordinate of the point in the region, ρ(r) is the density at position r, Γ g Indicates the measurement boundary;

[0030] The calculation formula for predicting the full waveform data of seismic waves d(r,t) is as follows:

[0031] At the detector measurement position r, let d(r,t) = u(r,t);

[0032] where u(r,t) is obtained by solving the wave equation We can get v as the wave velocity function, δ(rr s ) represents the Dirac Delta function, r s is the known source location, and ψ(t) is the known time wavelet function.

[0033] Preferably, the energy loss function E in S5 is j as follows:

[0034] E j =w×E g +E d +L(φ,ρ0,v1,v2);

[0035] The data fitting term of the energy loss function is: E g +E d ;

[0036] In the formula, E g is the gravity data fitting term, g z,i is the predicted gravity anomaly data calculated based on the density model at the i measurement point, is the gravity anomaly measurement data at measurement point i;

[0037] E d is the full waveform data fitting term, N s 、N r 、N T are the number of sources, the number of geophones, and the total number of measured time steps, represents the full waveform measurement data generated by source i received by detector j at time t, d i,j,t is the corresponding prediction data calculated based on the wave velocity model;

[0038] w is the weight coefficient of the balanced gravity data and the seismic wave full waveform data: w = w1w2, where w1 is the iterative attenuation term and w2 is the level set function gradient ratio term. The specific form is as follows:

[0039] (i) w1(k) = w0e -λk , w0 and λ are positive constants, k represents the kth iteration;

[0040] (ii) or

[0041] L(φ,ρ0,v1,v2) represents the regularization term introduced to the parameters φ,ρ0,v1,v2, and its expression is:

[0042]

[0043] where λ i ,1≤i≤4 represents the regularization term weight, ||·||2 represents the 2-norm;

[0044] For parameters that have been determined based on prior information, there is no need to regularize the corresponding parameters: for example, when the relative density ρ0 of the target medium and the wave velocity v1 of the target medium are known based on prior information, let λ2=λ3=0.

[0045] Preferably, the updated gradient in S6 includes the gradients of parameters φ, ρ0, v1, and v2;

[0046]

[0047]

[0048] For parameters that have been determined based on prior information, there is no need to calculate their gradients.

[0049] Preferably, the iterative update formula obtained based on the asynchronous iterative Adam algorithm in S7 is as follows:

[0050]

[0051] Wherein, θ∈(φ,ρ0,v1,v2) represents the parameters of the density model and the velocity model, and the subscript k represents the k-th iteration. For the parameters determined according to the prior information, no iteration is required;

[0052] α θ is the asynchronous iteration coefficient. For different model parameters θ, different α can be selected. θ The value of

[0053] ∈ is the basic learning rate, δ is a small constant; Where m k =β1m k-1 +(1-β1)G k , v k =β2v k-1 +(1-β2)G k ⊙G k , The gradient of the parameter θ corresponding to the k-th iteration, ⊙ represents the multiplication of the corresponding components, β1 and β2 are constants in [0,1), representing the decay rates of the first-order moment and the second-order moment, respectively.

[0054] On the other hand, the present invention also discloses a joint imaging system based on the full waveform of seismic waves and gravity data, which is used to implement the above-mentioned joint imaging method based on the full waveform of seismic waves and gravity data, including an information acquisition module, a model building module, an initialization module, a data prediction module, a judgment module, a calculation module, and an iterative update module connected in sequence; at the same time, the data prediction module is connected to the iterative update module;

[0055] It also includes an output module, which is connected to the judgment module;

[0056] Information acquisition module, used to obtain prior information of the target area, gravity anomaly measurement data and seismic wave full waveform measurement data;

[0057] Model building module, used to build density model and wave velocity model in combination with level set function;

[0058] An initialization module, used to initialize the level set function and initialize the parameters of the density model and the wave velocity model based on the prior information;

[0059] A data prediction module, used for calculating and predicting gravity anomaly data and seismic wave full waveform data according to the parameters of the density model and the wave velocity model;

[0060] A judgment module, used to construct an energy loss function using the gravity anomaly measurement data, the amplitude measurement data, the predicted gravity anomaly data and the predicted amplitude data and calculate the energy loss function value; if the energy loss function value is less than the target value, enter the output module, otherwise enter the calculation module;

[0061] The calculation module is used to calculate the update gradient;

[0062] An iterative update module, configured to use an asynchronous iterative Adam algorithm to update the parameters of the level set function, the density model, and the wave velocity model based on the updated gradient, and apply the updated parameters to a data prediction module;

[0063] The output module is used to output the level set function, density model, and wave velocity model after parameter update.

[0064] It can be seen from the above technical solutions that the present invention discloses a joint imaging method and system based on seismic wave full waveform and gravity data, which has the following beneficial effects:

[0065] 1. Use the level set method to describe the interface structure of the target area, characterize the structural similarity by sharing the level set function, link the density model and the wave velocity model, and organically combine the multi-physical information of gravity data and seismic wave full waveform data. Linking the density model and the wave velocity model based on structural similarity is more universal than the empirical formula of physical properties.

[0066] 2. Through the implicit boundary representation of the level set, the problem of the difficulty in accurately depicting complex geological boundaries in traditional methods is effectively solved, thereby improving the resolution of the interface structure of the joint imaging of the seismic wave full waveform and gravity data.

[0067] 3. In joint imaging, the level set expression can make full use of the prior information of model parameters and fix the model parameters known according to the prior information, which effectively reduces the interference of non-uniqueness problems on the inversion results, making the accuracy and reliability of the inversion results significantly better than the existing technical solutions.

[0068] 4. In joint imaging, the discontinuous interface structure is characterized by the zero level set of the continuous level set function, avoiding the direct introduction of discontinuous model parameters in the velocity model and density model. In the inversion calculation of joint imaging, the model parameters can be regularized by the L2 norm, which increases the smoothness and numerical stability of the optimization process and is more convenient for calculation.

[0069] 5. In the construction of the energy loss function of joint imaging, a scheme of weight coefficients balancing gravity data and seismic wave full waveform data is creatively proposed, so that the two types of physical data can be used reasonably, greatly improving the results of joint imaging.

[0070] 6. The framework design based on the level set method has good scalability, and can be further combined with deep learning methods to automatically extract geological features and optimize the inversion process. Compared with traditional methods, it not only improves the adaptability of the algorithm, but also provides innovative ideas for dealing with larger-scale and more diverse geological problems. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the provided drawings without paying creative work.

[0072] Figure 1 A flow chart of the method provided by the present invention;

[0073] Figure 2 This is a schematic diagram of gravity anomaly measurement data;

[0074] Figure 3(a)-Figure 3(d) Schematic diagram of full waveform measurement data of seismic waves corresponding to the earthquake sources at positions (2.3, 0), (6.9, 0), (11.5, 0), and (13.34, 0);

[0075] Figure 4(a) is the initialized density model, and Figure 4(b) is the initialized wave velocity model;

[0076] Figure 5(a) and Figure 5(b) are the density model and wave velocity model finally obtained by joint imaging, respectively;

[0077] Figure 6 This is a system architecture diagram provided by the present invention. DETAILED DESCRIPTION

[0078] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0079] This embodiment provides a combined imaging display of a submarine salt dome model.

[0080] First, let’s analyze some nouns:

[0081] Heaviside function: H(x)=1,x≥0; H(x)=0,x<0;

[0082] In the numerical calculation, we choose

[0083]

[0084] Its derivative is

[0085]

[0086] Delta function: generally refers to the Dirac delta function. The Dirac delta function is a generalized function. In physics, it is often used to represent the density distribution of ideal models such as particles and point charges. The function is equal to zero at all points except zero, and its integral over the entire domain is equal to 1. The Dirac delta function is the derivative of the Heaviside function. In calculation, we need to approximate it with a continuous function. We use:

[0087]

[0088] The embodiment of the present invention discloses a joint imaging method based on the full waveform of seismic waves and gravity data. The specific operation steps of the embodiment are described as follows:

[0089] S1. Obtain prior information, gravity anomaly measurement data and seismic full waveform measurement data of the target area.

[0090] Target area Ω0 = [0, 13.4 km] × [0, 4 km], Δx = Δz = 0.04 km, imaging matrix size n y ×n x =101×336.

[0091] The acquisition configuration of the acquired gravity anomaly measurement data is as follows: The measurement point position of the gravity data is set at the measurement boundary Γ of x0∈[-13.00km,27.00km],z=-0.2km g On, Δx m = 0.5km, with a total of 81 measurement points. Gravity anomaly measurement data like Figure 2 shown.

[0092] The acquisition configuration of the seismic wave full waveform data is as follows: there are 30 seismic sources, evenly distributed on the z = 0 km horizontal line, the interval between seismic sources is about 0.46 km (11.5 grid points), ranging from x = 0 km to x = 13.34 km. The receivers are evenly distributed on the z = 0 km horizontal line, recording seismic wave signals with a frequency of 5 Hz; the sampling interval Δt = 0.004 s, a total of 1800 time steps are recorded, corresponding to a total time T = 1800 × 0.004 = 7.2 s. The seismic full waveform measurement data is measured with d * express. Figure 3(a)-Figure 3(d) The full waveform measurement data of earthquakes corresponding to the locations of 4 out of 30 earthquake sources are shown: s =(2.3, 0), (6.9, 0), (11.5, 0), (13.34, 0).

[0093] S2. Construct density model and wave velocity model by combining level set function;

[0094] a. Define the level set function φ(r), where the zero level set {r|φ(r)=0} represents the location of the salt interface;

[0095] b. Construct density model and wave speed model:

[0096] ρ=H τ (φ)×ρ0;

[0097] v=v1×H τ (φ)+v2×(1-H τ (φ));

[0098] In this embodiment, the relative density ρ0 of the salt body and the wave velocity v1 in the salt body are given by prior information:

[0099] ρ0=(1.80-z)·0.2g / cm 3 , z represents the depth direction coordinate; v1=4.482km / s.

[0100] The above expression links the density model ρ and the wave velocity model v of the target area through the level set function φ. The zero level set {r|φ(r)=0} characterizes their common interface structure, that is, the structural similarity. In the above expression, some parameters can be fixed based on the known prior information of the target area and the target medium.

[0101] S3. Initialize the level set function;

[0102] Initialize the level set function to a relative ellipse Signed distance function, where a=4, b=1, (x c ,z c)=(6.7,2.0), the unit is km; initialize v2=1.5km / s; the initial density model and wave speed model are shown in Figure 4(a) and Figure 4(b);

[0103] S4. Calculate and predict gravity anomaly data and full seismic wave waveform data based on the parameters of the density model and the wave velocity model.

[0104] a. Calculate gravity anomaly data according to the following formula:

[0105]

[0106] in is the component of the gravity kernel function in the height direction, is the position vector of the measurement point, and r is the position vector of the point in the region. The calculation formula is:

[0107] b. Solve the following equation based on the existing wave velocity model to calculate the predicted seismic wave full waveform data:

[0108]

[0109] At the geophone measurement position r, let d(r, t) = u(r, t) to obtain the predicted seismic wave full waveform data;

[0110] S5. Use gravity anomaly measurement data, seismic wave full waveform measurement data, predicted gravity anomaly data and predicted seismic wave full waveform data to construct an energy loss function and calculate the energy loss function value; if the data fitting term in the energy loss function is less than the target value, execute S9, otherwise execute S6.

[0111] In this embodiment, since v1,ρ0 are fixed by prior information, the loss function E j Such as:

[0112] E j =w×E g +E d +L(φ,v2);

[0113] In the formula, E g is the gravity data fitting term, E d is the full waveform data fitting term:

[0114]

[0115] The data fitting term of the energy loss function is: E g +E d ,In the subsequent iteration process, when the data fitting term is less than the target threshold, the iteration stops;

[0116] L(φ,v2) represents the regularization term introduced for the parameter φ,v2, and its expression is:

[0117]

[0118] In this embodiment, λ1=10 -4 ,λ2=10 -4 ;

[0119] w is the weight coefficient of the balanced gravity data and the seismic wave full waveform data, w = w1w2. The specific selection of weight is as follows:

[0120] 1.

[0121] in is the number of iterations;

[0122] The present invention finds that taking w0≥1 makes the contribution of gravity data dominant in the initial iteration, which helps to quickly obtain the imaging structure of the deep area; introducing e -λk Allowing the specific gravity to decay will help make the inversion results more detailed and accurate.

[0123] 2. or

[0124] The principle of designing this item of the present invention is to make the gravity data and the seismic wave data contribute equally to the iteration of the shared level set function φ, so as to balance the different dimensions of the two physical data.

[0125] S6. Calculate the updated gradient, including the updated gradient of the level set function φ and the wave velocity v2 outside the salt body.

[0126] For the level set function φ, its updated gradient consists of three parts:

[0127] The gravity data fitting part is:

[0128]

[0129] Full waveform data fitting item part:

[0130]

[0131] Regularization term part:

[0132]

[0133] For the wave velocity v2(r) outside the salt body, its update gradient consists of two parts:

[0134]

[0135] In the above formula, Through deepwave computing, This is accomplished by automatic differentiation.

[0136] S7. Use the asynchronous iterative Adam algorithm to update the parameters of the level set function, density model, and wave velocity model based on the updated gradient.

[0137] The asynchronous iterative Adam algorithm is used for the inversion parameter φ,v2, and the iterative update formula is:

[0138]

[0139] Where θ∈(φ,v2), α θ Represents the correction coefficient, which is used to balance the iteration rate between φ and υ2; ∈ is the basic learning rate, and δ is a small constant; Where m k =β1m k-1 +(1-β1)G k , v k =β2v k-1 +(1-β2)G k ⊙G k , Corresponding to the gradient of parameter θ at the k-th iteration; in this embodiment,

[0140] S8. Return to S4;

[0141] S9. Output the updated level set function, density model, and wave velocity model.

[0142] The density model and wave velocity model reconstruction images finally obtained in this embodiment are shown in Figure 5(a) and Figure 5(b) respectively.

[0143] On the other hand, the present invention also discloses a combined imaging system based on seismic wave full waveform and gravity data (refer to Figure 6 ), used to implement the above-mentioned joint imaging method based on seismic wave full waveform and gravity data, including an information acquisition module, a model building module, an initialization module, a data prediction module, a judgment module, a calculation module, and an iterative update module connected in sequence; at the same time, the data prediction module is connected to the iterative update module;

[0144] It also includes an output module, which is connected to the judgment module;

[0145] Information acquisition module, used to obtain prior information of the target area, gravity anomaly measurement data and seismic wave full waveform measurement data;

[0146] Model building module, used to build density model and wave velocity model in combination with level set function;

[0147] Initialization module, used to initialize the level set function and initialize the parameters of the density model and wave velocity model based on prior information;

[0148] A data prediction module is used to calculate and predict gravity anomaly data and full seismic wave waveform data based on the parameters of the density model and the wave velocity model;

[0149] A judgment module is used to construct an energy loss function and calculate the energy loss function value by using gravity anomaly measurement data, seismic wave full waveform measurement data, predicted gravity anomaly data and predicted seismic wave full waveform data; if the energy loss function value is less than the target value, it enters the output module, otherwise it enters the calculation module;

[0150] The calculation module is used to calculate the update gradient;

[0151] Iterative update module, used to use the asynchronous iterative Adam algorithm to update the parameters of the level set function, density model and wave velocity model based on the updated gradient, and apply the updated parameters to the data prediction module;

[0152] The output module is used to output the level set function, density model, and wave velocity model after parameter update.

[0153] In this specification, each embodiment is described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the embodiments can be referred to each other. For the device disclosed in the embodiment, since it corresponds to the method disclosed in the embodiment, the description is relatively simple, and the relevant parts can be referred to the method part.

[0154] The above description of the disclosed embodiments enables one skilled in the art to implement or use the present invention. Various modifications to these embodiments will be apparent to one skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to the embodiments shown herein, but rather to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A joint imaging method based on full seismic waveform and gravity data, characterized in that: The following steps are involved: S1. Acquire prior information of the target area, gravity anomaly measurement data and seismic wave full waveform measurement data; S2. Construct density model and wave velocity model by combining level set function, and connect the inversion parameters corresponding to different physical data through interface structure; S3. Initializing the level set function, and initializing the parameters of the density model and the wave velocity model based on the prior information; S4. Calculating predicted gravity anomaly data and predicted seismic wave full waveform data according to the parameters of the density model and the wave velocity model; S5. construct an energy loss function using the gravity anomaly measurement data, the seismic wave full waveform measurement data, the predicted gravity anomaly data and the predicted seismic wave full waveform data; if the data fitting term in the energy loss function is less than the target value, execute S9, otherwise execute S6; S6. Calculate the updated gradient; S7. Using an asynchronous iterative Adam algorithm, updating the parameters of the level set function, the density model and the wave velocity model based on the updated gradient; S8. Return to S4; S9. Output the updated level set function, density model, and wave velocity model.

2. A joint imaging method based on seismic wave full waveform and gravity data according to claim 1, characterized in that: The S2 includes: Define the level set function φ, and characterize the common interface structure through the zero level set {r|φ(r)=0}; The density model ρ and wave velocity model v are constructed by combining the level set function as follows: p=H τ (φ)×ρ0; v=v1×H τ (φ)+v2×(1-H τ (f)); Where ρ0 is the relative density of the target medium to the background area, v1 is the wave velocity of the target medium, v2 is the wave velocity outside the target medium area, and H τ (·) is the smooth approximation function of the Heaviside function.

3. A joint imaging method based on seismic wave full waveform and gravity data according to claim 2, characterized in that: S3 includes: The level set function φ is initialized as a signed distance function of an ellipse to the center of the target area; The parameters of the density model and the wave velocity model are initialized based on prior information, including the relative density ρ0 of the target medium to the background area, the wave velocity v1 of the target medium, and the wave velocity v2 outside the target medium area.

4. The joint imaging method based on seismic wave full waveform and gravity data according to claim 1, characterized in that: Predicted gravity anomaly data in S4 The calculation formula is as follows: in is the component of the gravity kernel function in the height direction, z are r is the height direction component, n represents the spatial dimension considered; is the position coordinate of the measurement point, r is the position coordinate of the point in the region, ρ(r) is the density at position r, Γ g Indicates the measurement boundary; The calculation formula for predicting the full waveform data of seismic waves d(r,t) is as follows: At the detector measurement position r, let d(r,t) = u(r,t); where u(r,t) is obtained by solving the wave equation We can get v as the wave velocity function, δ(rr s ) represents the Dirac Delta function, r s is the known source location, and ψ(t) is the known time wavelet function.

5. The joint imaging method based on seismic wave full waveform and gravity data according to claim 2, characterized in that: The energy loss function E described in S5 j as follows: E j =w×E g +E d +L(φ,ρ0,v1,v2); The data fitting term of the energy loss function is: E g +E d ; In the formula, E g is the gravity data fitting term, g z,i is the predicted gravity anomaly data calculated based on the density model at the i measurement point, is the gravity anomaly measurement data at measurement point i; E d is the full waveform data fitting term, N s 、N r 、N T are the number of sources, the number of geophones, and the total number of measured time steps, represents the full waveform measurement data generated by source i received by detector j at time t, d i,j,t is the corresponding prediction data calculated based on the wave velocity model; w is the weight coefficient of the balanced gravity data and the seismic wave full waveform data: w = w1w2, where w1 is the iterative attenuation term and w2 is the level set function gradient ratio term. The specific form is as follows: (i) w1(k) = w0e -λk , w0 and λ are positive constants, k represents the kth iteration; (ii) or L(φ,ρ0,v1,v2) represents the regularization term introduced for parameters φ,ρ0,v1,v2. For parameters that have been determined based on prior information, there is no need to regularize the corresponding parameters.

6. The joint imaging method based on seismic wave full waveform and gravity data according to claim 2, characterized in that: The updated gradient in S6 includes the gradients of parameters φ, ρ0, v1, and v2; For parameters that have been determined based on prior information, there is no need to calculate their gradients.

7. The joint imaging method based on seismic full waveform and gravity data according to claim 2, characterized in that: The iterative update formula based on the asynchronous iterative Adam algorithm in S7 is as follows: Wherein, θ∈(φ,ρ0,v1,v2) represents the parameters of the density model and the velocity model, and the subscript k represents the k-th iteration. For the parameters determined according to the prior information, no iteration is required; α θ is the asynchronous iteration coefficient. For different model parameters θ, different α can be selected. θ The value of ∈ is the basic learning rate, δ is a small constant; Where m k =β1m k-1 +(1-β1)G k , v k =β2v k-1 + ( 1-β2)G k ⊙G k , The gradient of the parameter θ corresponding to the k-th iteration, ⊙ represents the multiplication of the corresponding components, β1 and β2 are constants in [0,1), representing the decay rates of the first-order moment and the second-order moment, respectively.

8. A joint imaging system based on seismic wave full waveform and gravity data, characterized in that: It includes an information acquisition module, a model building module, an initialization module, a data prediction module, a judgment module, a calculation module, and an iterative update module which are connected in sequence; and the data prediction module is connected to the iterative update module; It also includes an output module, which is connected to the judgment module; Information acquisition module, used to obtain prior information of the target area, gravity anomaly measurement data and seismic wave full waveform measurement data; Model building module, used to build density model and wave velocity model in combination with level set function; An initialization module, used to initialize the level set function and initialize the parameters of the density model and the wave velocity model based on the prior information; A data prediction module, used for calculating and predicting gravity anomaly data and seismic wave full waveform data according to the parameters of the density model and the wave velocity model; A judgment module, used to construct an energy loss function and calculate the energy loss function value by using the gravity anomaly measurement data, the seismic wave full waveform measurement data, the predicted gravity anomaly data and the predicted seismic wave full waveform data; If the energy loss function value is less than the target value, it enters the output module, otherwise it enters the calculation module; The calculation module is used to calculate the update gradient; An iterative update module, configured to use an asynchronous iterative Adam algorithm to update the parameters of the level set function, the density model, and the wave velocity model based on the updated gradient, and apply the updated parameters to a data prediction module; The output module is used to output the level set function, density model, and wave velocity model after parameter update.

Images