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

Abstract

The application discloses a kind of based on the joint imaging method and system of seismic wave full waveform and gravity data, belong to geological exploration technical field.First, preliminary modeling is carried out to research area to obtain prior information, then data is collected to construct total energy loss function, density and wave velocity model are constructed based on level set method, interface structure is linked to different physical data corresponding inversion parameter, prior information is initialized model parameter, calculated prediction data, determine data proportion coefficient, judge energy loss, calculate update gradient, use asynchronous iteration Adam algorithm to update parameter, complete joint inversion.The application effectively fuses multiple physical data, reasonably sets data proportion, improves imaging effect by the aid of structure similarity, fully utilizes prior information and optimizes regularization term, has good scalability, has strong ability to capture complex interface, can improve geological structure imaging resolution, reduce interference, reduce calculation cost, improve the precision, reliability and overall quality of inversion result.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geological exploration, more particularly to a joint imaging method and system based on seismic wave full waveform and gravity data. BACKGROUND

[0002] In the deep exploration of geological exploration and resource development, accurately depicting the underground geological structure is the core element to improve the efficiency and accuracy of exploration. 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 survey. However, the accuracy of this technology is easily affected by various factors, such as noise of observation data, non-uniqueness of model and uncertainty of density distribution, etc., thus there are certain limitations in high-resolution applications.

[0003] Relatively, full waveform inversion technology simulates the propagation characteristics of seismic waves to invert the wave velocity model of underground geological structure, which theoretically has the advantages of high resolution and high accuracy. However, in practical applications, full waveform inversion faces many challenges, especially in dealing with complex geological structures such as salt body model, the scattering and reflection of seismic waves are significant, which often leads to the deviation of the inversion result from the true situation, and it is difficult to meet the demand of high-precision exploration. In view of the limitations of full waveform inversion, researchers have explored improved technologies such as multi-scale inversion and adaptive regularization to enhance the analytical ability of complex geological structures. Although some progress has been made, there is still a wide space for improvement in full waveform inversion technology, and further exploration and innovation are still needed in how to further improve the efficiency and accuracy of inversion.

[0004] Under this background, joint inversion method gradually becomes an important way to solve the above problems. Among them, the gravity and seismic wave data joint inversion technology uses the complementarity of the two kinds of data to significantly improve the imaging ability of deep region and complex geological structure. The key of joint inversion is to establish the connection between different physical parameters, so that the two kinds of data can complement each other. However, the existing gravity and seismic wave joint inversion mainly connects the wave velocity and density model through the physical property empirical formula, and cannot effectively utilize the structural similarity in the wave velocity and density model in 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 needs to be solved by the person skilled in the art. SUMMARY

[0006] Therefore, the present application 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 application adopts the following technical solutions:

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

[0009] S1. Obtain prior information, gravity anomaly measurement data and seismic wave full waveform measurement data of a target region;

[0010] S2. Construct a density model and a wave velocity model in combination with a level set function, and link the inversion parameters corresponding to different physical data through an interface structure;

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

[0012] S4. Calculate 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 a data fitting term in the energy loss function is less than a target value, execute S9, otherwise execute S6;

[0014] S6. Calculate an update gradient;

[0015] S7. Update the parameters of the level set function, the density model and the wave velocity model based on the update gradient using an asynchronous iteration Adam algorithm;

[0016] S8. Return to S4;

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

[0018] Preferably, S2 comprises:

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

[0020] Construct a density model ρ and a wave velocity model v in combination with the level set function as follows:

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

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

[0023] In the formula, ρ0 is the relative density of the target medium (e.g., a salt body) relative to the background region, v1 is the wave velocity of the target medium, v2 is the wave velocity outside the target medium region, and H... τ (·) is a smooth approximation of the Heaviside function.

[0024] Preferably, S3 includes:

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

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

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

[0028]

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

[0030] The formula for calculating the full waveform data d(r,t) of the predicted seismic wave is as follows:

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

[0032] Where u(r,t) is obtained by solving the wave equation We obtain v as the wave speed function, δ(rr) s ) represents the Dirac Delta function, r s Given the known location of the earthquake source, ψ(t) is a wavelet function with a known time.

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

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

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

[0036] E g is a gravity data fitting term, g z,i is predicted gravity anomaly data calculated according to the density model at the i measuring point, is gravity anomaly measured data at the i measuring point;

[0037] E d is a full waveform data fitting term, N s , N r , N T are the number of sources, the number of receivers and the total number of measured time steps, respectively, represents the full waveform measured data generated by the i source received by the j receiver at the t time, d i,j,t is the corresponding predicted data calculated according to the wave velocity model;

[0038] w is the proportionality coefficient of balancing the gravity data and the seismic wave full waveform data: w = w1w2, wherein w1 is an iterative attenuation term and w2 is a level set function gradient ratio term, and the specific form is as follows:

[0039] (i) w1(k) = w0e -λk , w0 and λ are normal numbers, and k represents the kth step of iteration;

[0040] (ii) or

[0041] L(φ,ρ0,v1,v2) represents a regularization term introduced for the parameters φ, ρ0, v1, v2, and the expression is as follows:

[0042]

[0043] wherein λ i , 1≤i≤4 represent regularization term weights, and ||·||2 represents a 2-norm;

[0044] For the parameters determined according to prior information, the regularization of the corresponding parameters is not required: for example, when the relative density ρ0 of the target medium and the wave velocity v1 of the target medium are known by prior information, then λ2 = λ3 = 0.

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

[0046]

[0047]

[0048] For the parameters determined according to the prior information, the gradient thereof does not need to be calculated.

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

[0050]

[0051] In the formula, θ∈(φ, ρ0, v1, v2) represents the parameters of the density model and the velocity model, and the subscript k represents the kth iteration step. For the parameters determined according to the prior information, iteration is not needed.

[0052] α θ is an asynchronous iteration coefficient, and different α θ values can be selected for different model parameters θ.

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

[0054] On the other hand, the application also discloses a joint imaging system based on seismic wave full waveform and gravity data, which is used for implementing the joint imaging method based on seismic wave full waveform and gravity data, and comprises an information acquisition module, a model construction module, an initialization module, a data prediction module, a judgment module, a calculation module, an iterative updating module connected in sequence; meanwhile, the data prediction module is connected with the iterative updating module.

[0055] The joint imaging system further comprises an output module connected with the judgment module.

[0056] The information acquisition module is used for acquiring prior information, gravity anomaly measurement data and seismic wave full waveform measurement data of a target region.

[0057] The model construction module is used for constructing a density model and a wave velocity model in combination with a level set function.

[0058] The initialization module is used for initializing the level set function, and initializing parameters of the density model and the wave velocity model based on the prior information.

[0059] a data prediction module configured to calculate predicted gravity anomaly data and predicted seismic wave full waveform data according to parameters of the density model and the wave velocity model;

[0060] a judging module configured to construct an energy loss function and calculate an energy loss function value by using gravity anomaly measurement data, the amplitude measurement data, the predicted gravity anomaly data and the predicted amplitude data; if the energy loss function value is less than a target value, the method enters an output module, otherwise, the method enters the calculating module;

[0061] the calculating module is configured to calculate an update gradient;

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

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

[0064] According to the technical solution, the application discloses a joint imaging method and system based on seismic wave full waveform and gravity data, which has the following beneficial effects:

[0065] 1. The level set method is used to describe the interface structure of the target region, the structure similarity is described by sharing the level set function, the density model and the wave velocity model are connected, and the multi-physical information of the gravity data and the seismic wave full waveform data is organically combined. The density model and the wave velocity model are connected based on the structure similarity, which is more universal than the physical property empirical formula.

[0066] 2. The implicit boundary representation of the level set effectively solves the problem that it is difficult to accurately describe the complex geological boundary in the traditional method, thereby improving the resolution of the interface structure in the joint imaging of the seismic wave full waveform and the gravity data.

[0067] 3. In the joint imaging, the level set expression can fully utilize the prior information of the model parameters, and the model parameters known according to the prior information are fixed, thereby effectively reducing the interference of the non-uniqueness problem on the inversion result, and the accuracy and reliability of the inversion result are obviously better than those of the prior art.

[0068] 4. In the joint imaging, the zero level set of the continuous level set function is used to describe the discontinuous interface structure, thereby avoiding directly introducing the discontinuous model parameters in the wave velocity model and the density model. In the inversion calculation of the joint imaging, the L2 norm regularization can be used for the model parameters, thereby increasing the smoothness and numerical stability of the optimization process, and the calculation is more convenient.

[0069] 5. In the joint imaging of energy loss function construction, a specific gravity coefficient scheme is creatively given to balance the gravity data and seismic wave full waveform data, so that the two kinds of physical data can be reasonably used, and the joint imaging result is greatly improved.

[0070] 6. The framework based on the level set method has good scalability, and subsequent deep learning methods can be further combined to automatically extract geological features and optimize the inversion process. Compared with traditional methods, the adaptability of the algorithm is improved, and an innovative idea is provided for processing larger-scale and more diversified geological problems. BRIEF DESCRIPTION OF DRAWINGS

[0071] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the drawings needed to be used in the embodiments or prior art description will be briefly introduced as follows. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor on the basis of the provided drawings.

[0072] Figure 1 The method flowchart provided by the present application is shown in the figure.

[0073] Figure 2 The gravity anomaly measurement data schematic diagram is shown in the figure.

[0074] Figures 3(a)-3(d) The seismic wave full waveform measurement data schematic diagram corresponding to the seismic source with position (2.3, 0), (6.9, 0), (11.5, 0), (13.34, 0) respectively is shown in the figure.

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

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

[0077] Figure 6 The system architecture diagram provided by the present application is shown in the figure. DETAILED DESCRIPTION

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

[0079] The present embodiment provides a joint imaging display of a submarine salt dome model.

[0080] First, several nouns are analyzed:

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

[0082] In numerical calculation, we select

[0083]

[0084] Its derivative is

[0085]

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

[0087]

[0088] The embodiment of the application discloses a joint imaging method based on seismic wave full waveform and gravity data, and 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 a target region.

[0090] The target region Omega0 = [0, 13.4km] x [0, 4km], Delta x = Delta z = 0.04km, and the imaging matrix size n y x n x = 101 x 336.

[0091] The obtained gravity anomaly measurement data is collected and configured as follows: the measurement point position of the gravity data is set in x0 ∈ [-13.00km, 27.00km], z = -0.2km, the measurement boundary Gamma g is above, Delta x m = 0.5km, and there are 81 measurement points. The gravity anomaly measurement data is shown as follows. Figure 2

[0092] ​The acquisition configuration of the acquired seismic full waveform data is as follows: 30 sources are uniformly distributed on the z=0 km horizontal line, the interval distance between the sources is about 0.46 km (11.5 grid points), and the range is from x=0 km to x=13.34 km. The receivers are uniformly distributed on the z=0 km horizontal line, and the recording frequency of the seismic wave signal is 5 Hz; the sampling interval is Δt=0.004 s, and 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 represented by d * . Figures 3(a)-3(d) The seismic full waveform measurement data corresponding to the positions of 4 sources in the 30 sources is shown: the source positions r s =(2.3, 0), (6.9, 0), (11.5, 0), (13.34, 0).

[0093] S2. Constructing a density model and a wave velocity model in combination with a level set function;

[0094] a. Defining a level set function φ(r), and the zero level set {r|φ(r)=0} represents the position of the salt body interface;

[0095] b. Constructing a density model and a wave velocity 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.2 g / cm 3 , z represents the depth direction coordinate; v1=4.482 km / s.

[0100] The above expressions link the density model ρ and the wave velocity model v of the target region through the level set function φ, and the zero level set {r|φ(r)=0} describes their common interface structure, that is, the structural similarity; in the above expressions, part of the parameters can be fixed based on the known prior information of the target region and the target medium.

[0101] S3. Initializing the level set function;

[0102] The level set function is initialized as a relative ellipse directed distance function, where a=4, b=1, (x c ,z c) = (6.7, 2.0) km; initialize v2 = 1.5 km / s; initial density model and wave velocity model are shown in Fig. 4(a), Fig. 4(b);

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

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

[0105]

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

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

[0108]

[0109] The predicted seismic wave full waveform data can be obtained by letting d(r, t) = u(r, t) at the measurement position r of the receiver;

[0110] 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, 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, p0 have been fixed by prior information, the loss function E j is of the form:

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

[0113] In the formula, E g is the gravity data fitting term, and 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 parameters φ, v2, and its expression is:

[0117]

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

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

[0120] 1.

[0121] wherein is the iteration number;

[0122] The present application finds that, by taking w0≥1, the contribution of gravity data is dominant in the initial iteration, which helps to quickly obtain the imaging structure of the deep region; by introducing e -λk , the proportionality decays, which helps to make the inversion result more refined and accurate.

[0123] 2. or

[0124] The principle of the present application is to make the iterative contributions of gravity data and seismic data to the level set function φ equivalent, so as to balance the different dimensions of the two kinds of physical data.

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

[0126] For the level set function φ, the update gradient thereof is composed of three parts:

[0127] wherein the gravity data fitting part is:

[0128]

[0129] The full waveform data fitting part is:

[0130]

[0131] The regularization part is:

[0132]

[0133] For the outer wave velocity v2(r) of the salt body, the update gradient thereof contains two parts:

[0134]

[0135] In the above formula, Realized by deepwave calculation, Realized by automatic derivation.

[0136] S7. Update the parameters of the level set function, density model and wave velocity model based on the updated gradient using the asynchronous iteration Adam algorithm.

[0137] Use the asynchronous iteration Adam algorithm to update the inversion parameters φ and v2, and the iterative update formula is:

[0138]

[0139] In the formula, θ∈(φ,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 , The gradient of the parameter θ corresponding to the kth 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 final density model and wave velocity model reconstruction maps obtained in this embodiment are shown in Figures 5(a) and 5(b) respectively.

[0143] On the other hand, the application also discloses a joint imaging system based on seismic wave full waveform and gravity data (refer to Figure 6 for realizing the joint imaging method based on seismic wave full waveform and gravity data, comprising an information acquisition module, a model construction module, an initialization module, a data prediction module, a judgment module, a calculation module and an iterative update module connected in sequence; the data prediction module is connected with the iterative update module;

[0144] Further comprising an output module, which is connected with the judgment module;

[0145] The information acquisition module is used to acquire the prior information, gravity anomaly measurement data and seismic wave full waveform measurement data of the target region.

[0146] The model construction module is configured to construct a density model and a wave velocity model in combination with a level set function.

[0147] The initialization module is configured to initialize the level set function, and initialize parameters of the density model and the wave velocity model based on prior information.

[0148] The data prediction module is configured to calculate predicted gravity anomaly data and predicted seismic wave full waveform data according to the parameters of the density model and the wave velocity model.

[0149] The judgment module is configured to construct an energy loss function and calculate an 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.

[0150] The calculation module is configured to calculate an update gradient.

[0151] The iterative update module is configured to update the level set function, the parameters of the density model and the wave velocity model based on the update gradient by using an asynchronous iteration Adam algorithm, and apply the updated parameters to the data prediction module.

[0152] The output module is configured to output the level set function, the density model and the wave velocity model after the parameters are updated.

[0153] The embodiments in the specification are described in a progressive manner, and each embodiment focuses on the difference from other embodiments, and the same or similar parts between the embodiments can be referred to each other.

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

Claims

1. A joint imaging method based on seismic wave full waveform and gravity data, characterized in that, The method comprises the following steps: S1. obtaining prior information of a target area, gravity anomaly measurement data and seismic wave full waveform measurement data; S2. constructing a density model and a wave velocity model in combination with a level set function, and connecting inversion parameters corresponding to different physical data through an interface structure; S3. initializing the level set function, and initializing 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. constructing an energy loss function 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 a data fitting term in the energy loss function is less than a target value, performing S9, otherwise performing S6; S6. calculating an update gradient; S7. updating the parameters of the level set function, the density model and the wave velocity model based on the update gradient using an asynchronous iteration Adam algorithm; S8. returning to S4; S9. outputting the level set function, the density model and the wave velocity model after parameter update.

2. The method of claim 1, wherein, The S2 comprises: defining a level set function by a zero level set depicting the common interface structure; Building density models in conjunction with level set functions and wave speed models As follows: ； ； wherein is the relative density of the target medium to the background region, is the wave speed of the target medium, is the wave speed outside the target medium region, is a smooth approximation function of the Heaviside function.

3. The method of claim 2, wherein, The S3 comprises: the level set function a distance function initialized to an ellipse centered on the target region; initializing parameters of the density model and the wave velocity model based on prior information, including relative density of the target medium to the background region , wave velocity of the target medium , and wave velocity outside the target medium region .

4. The method of claim 1, wherein, Predicting gravity anomaly data in S4 The formula for calculating is as follows: ； in It is the component of the gravity kernel function in the height direction. , They are respectively The height direction component, Indicates the spatial dimensions considered; These are the coordinates of the measurement point. These are the position coordinates of a point within the region. for r Density at location Indicates the measurement boundary; Predicting seismic wave full waveform data The formula for calculating is as follows: At the geophone measurement location Let ; where By solving the wave equation ( r - ) , is the wave velocity function, ( r - ) denotes the Dirac delta function, is the known source location, is the known time wavelet function.

5. The method of claim 2, wherein, The energy loss function described in S5 As follows: ； The data fitting term of the energy loss function is: ; wherein is a gravity data fitting term, , is a gravity data fitting term, predicted gravity anomaly data calculated from the density model at the measurement point, is a gravity data fitting term, gravity anomaly measurement data at the measurement point; For fitting the full waveform data, , These represent the number of seismic sources, the number of geophones, and the total number of measurement time steps, respectively. represent j Detector in t Received in real time i Full waveform measurement data generated by the earthquake source, These are the corresponding predicted data calculated based on the wave velocity model; is the specific gravity coefficient balancing the gravity data and the seismic wave full waveform data: wherein is an iterative damping term, is a level set function gradient ratio term, in particular as follows: (i) , and are normal numbers, denotes the th iteration; (ii) or ; denotes the regularization of the parameters The introduced regularization term is not needed for parameters for which the prior information has determined the parameter.

6. The method of claim 2, wherein, The gradients in S6 include gradients of parameters of the parameters; ； ； ； ； For parameters determined according to prior information, there is no need to calculate the gradient thereof.

7. The method of claim 2, wherein, The iteration update formula obtained in S7 based on the asynchronous iteration Adam algorithm is as follows: ； wherein denote parameters of the density model and the wave velocity model, the index denote the parameters of the density model and the wave velocity model, the index step iteration, for parameters determined from prior information no iteration is needed; For asynchronous iteration coefficients, different values of are chosen for different model parameters ; a base learning rate, a small constant; wherein , , corresponds to the gradient of the parameters at the and are constants within the interval [0, 1] representing the decay rates of the first and second moments, respectively.

8. A joint imaging system based on seismic wave full waveform and gravity data, characterized in that, The device comprises an information acquisition module, a model construction module, an initialization module, a data prediction module, a judgment module, a calculation module and an iteration update module connected in sequence; and the data prediction module is connected with the iteration update module. The device further comprises an output module connected with the judgment module. The information acquisition module is configured to acquire prior information of a target area, gravity anomaly measurement data and seismic wave full waveform measurement data. The model construction module is configured to construct a density model and a wave velocity model in combination with a level set function. The initialization module is configured to initialize the level set function, and initialize parameters of the density model and the wave velocity model based on the prior information. The data prediction module is configured to calculate predicted gravity anomaly data and predicted seismic wave full waveform data according to the parameters of the density model and the wave velocity model. The judgment module is configured to construct an energy loss function 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, and calculate an energy loss function value. If the energy loss function value is less than a target value, the output module is entered, otherwise the calculation module is entered. The calculation module is configured to calculate an update gradient. The iteration update module is configured to update the parameters of the level set function, the density model and the wave velocity model based on the update gradient using an asynchronous iteration Adam algorithm, and apply the updated parameters to the data prediction module. The output module is configured to output the level set function, the density model and the wave velocity model after parameter update.

Images