Earth crust three-dimensional density modeling method and system based on earthquake constraint

Abstract

The invention provides an earth crust three-dimensional density modeling method and system based on seismic constraints, and relates to the field of earth crust three-dimensional density modeling and the field of sea area oil-gas exploration, and the method comprises the steps: obtaining a two-dimensional reflection seismic section and a two-dimensional refraction seismic velocity model in a research area, constructing constraint information of a main density horizon, and obtaining a two-dimensional reflection seismic section and a two-dimensional refraction seismic velocity model; performing anomaly extraction and depth inversion to obtain the depth of the main density interface; and according to the depth of the main density interface and the density constraint information, establishing a three-dimensional density reference model and an upper and lower limit constraint model conforming to the prior constraint information, performing three-dimensional constraint inversion by adopting a regularization inversion method, and constructing an earth crust three-dimensional density model. Two-dimensional constraint information is introduced into three-dimensional inversion, so that the reliability of an inversion result is improved.

Description

Technical Field

[0001] The invention relates to the technical field of crust three-dimensional density modeling and oil and gas exploration, and in particular to a method and system for crust three-dimensional density modeling based on seismic constraints. Background Art

[0002] As the ocean has become the main growth area of ​​global oil and gas, studying the crust structure of the sea to determine favorable oil and gas areas can play a huge role in promoting breakthroughs in oil and gas exploration. In addition, the formation of the ocean is an important product of the evolution of the earth. Studying the crust structure of the ocean will help reveal the internal mechanism of the earth's evolution.

[0003] Different marine geological units often have significant density differences, and the special geological effects of oil and gas accumulation will also cause changes in density, providing a physical basis for conducting marine crustal structure research and oil and gas exploration based on gravity data; satellite gravity observation technology provides regional data with high coverage and high accuracy of low-frequency components, which has certain advantages in crustal scale research; but due to the volume effect of potential field, poor vertical resolution and non-uniqueness of geophysical inversion problems, its application in quantitative inversion of crustal scale and direct determination of favorable oil and gas areas is limited.

[0004] In three-dimensional gravity inversion, two-dimensional seismic profiles provide rich constraint information for inversion; however, due to their sparsity, two-dimensional seismic profiles are difficult to apply to the construction of three-dimensional reference models; existing methods are only used indirectly to determine the optimal inversion parameters and to evaluate the reliability of the results after the inversion is completed. This method does not fully utilize the constraint effect of seismic data; how to introduce two-dimensional constraint information into three-dimensional inversion is the key to improving the reliability of inversion results. Summary of the invention

[0005] The present invention provides a method and system for three-dimensional crust density modeling based on earthquake constraints, which introduces two-dimensional constraint information into three-dimensional inversion to improve the reliability of inversion results.

[0006] In order to solve the above technical problems, the technical solution of the present invention is as follows:

[0007] In a first aspect, a method for three-dimensional crustal density modeling based on seismic constraints, the method comprising:

[0008] Obtaining a two-dimensional reflection seismic profile and a two-dimensional refraction seismic velocity model in the study area, and constructing constraint information of main density layers, wherein the constraint information of the main density layers includes geometric constraint information and density constraint information;

[0009] Based on the satellite gravity anomaly, the main density interface is extracted and depth inverted according to the geometric constraint information and the density constraint information to obtain the depth of the main density interface;

[0010] According to the depth of the main density interface and the density constraint information, a three-dimensional density reference model and upper and lower limit constraint models that conform to the prior constraint information are established;

[0011] Based on satellite gravity anomalies, according to the three-dimensional density reference model and the upper and lower limit constraint models, the regularized inversion method is used to perform three-dimensional constrained inversion and construct a three-dimensional density model of the crust.

[0012] Furthermore, obtain the 2D reflection seismic profile and 2D refraction seismic velocity model in the study area and construct the constraint information of the main density layers, including:

[0013] Obtain 2D reflection seismic profiles and 2D refraction seismic velocity models in the study area;

[0014] Interpret the stratigraphic structure of 2D reflection seismic profiles and 2D refraction seismic velocity models to identify the main interfaces, such as the basement interface between sedimentary layers and the crust and the Moho between the crust and the mantle;

[0015] The depths of the basement interface and the Moho surface are discretized and digitized, and the geometric constraint information of the basement interface and the Moho surface is established.

[0016] Furthermore, the two-dimensional reflection seismic profile and the two-dimensional refraction seismic velocity model in the study area are obtained to construct the constraint information of the main density layers, which also includes:

[0017] Extract seismic P-wave velocity information of sedimentary layers, crust and mantle in different tectonic regions based on the two-dimensional refraction seismic velocity model;

[0018] According to the velocity-density conversion relationship, seismic P-wave velocity information is converted into density parameters, and density constraint information of sedimentary layers, crust and mantle in different regions is established, wherein the density constraint information includes maximum value, minimum value and average value;

[0019] According to the density constraint information of sedimentary layer, crust and mantle, the lateral variable density constraint information of basement interface and the lateral variable density constraint information of Moho surface are established.

[0020] Furthermore, based on the satellite gravity anomaly, the main density interface is subjected to anomaly extraction and depth inversion according to the geometric constraint information and the density constraint information to obtain the depth of the main density interface, including:

[0021] Obtain satellite gravity anomaly data, seafloor bathymetry data, and sediment thickness data;

[0022] According to the seabed bathymetric data, the seawater layer is divided into a series of right hexahedrons, and the spatial domain method is used for forward modeling to obtain the seawater correction value;

[0023] According to the sediment thickness data, the sedimentary layer is divided into a series of right hexahedrons. Based on the compaction model, the spatial domain method is used for forward modeling to obtain the sedimentary layer correction value. The density at different depths is:

[0024]

[0025] Among them, g0 is the gravity anomaly of the Moho surface, ρz is the density of the sediment layer at depth z, ρf is the density of the fluid in the sediment, ρg is the density of the rock particles, is the initial porosity, c is the porosity depth attenuation coefficient;

[0026] The satellite gravity anomaly data is corrected according to the seawater correction value and the sediment layer correction value to obtain the crustal Bouguer gravity anomaly data;

[0027] Based on the geometric constraint information of Moho surface, the potential field of crustal Bouguer gravity anomaly data is separated to obtain Moho surface gravity anomaly data;

[0028] Based on the geometric constraint information of the Moho surface, the initial Moho surface depth estimation pmoho0 is obtained using the linear regression method;

[0029] According to the Moho surface gravity anomaly data, the Moho surface geometric constraint information is introduced in a soft constraint manner. Based on the Moho surface lateral variable density constraint information, the depth of the Moho surface is iteratively calculated. The formula is:

[0030] pmohok=pmohok-1+gmoho0-gmohok-12πGρmoho+Δpmoho,

[0031] Where gmoho0 is the Moho gravity anomaly, gmohok-1 is the forward gravity response of the k-1th Moho depth, pmo h o k and pmo h o k-1 are the results of the kth and k-1th iterations respectively, G is the gravity constant, ρmoho is the lateral variable density information of the Moho surface, and Δpmoho is the depth correction value obtained by linear regression based on the geometric constraint information of the Moho surface for the kth time.

[0032] Further, based on the satellite gravity anomaly, the main density interface is extracted and depth inverted according to the geometric constraint information and the density constraint information to obtain the depth of the main density interface, and further includes:

[0033] According to the depth of the Moho surface, forward modeling is performed to obtain the Moho surface correction value;

[0034] Based on the satellite gravity anomaly data and the seabed bathymetry data, the correlation analysis was carried out to obtain the optimal density parameters for seawater correction. The seawater layer was divided into a series of right hexahedrons using the seabed bathymetry data, and the seawater correction value was obtained by forward calculation using the spatial domain method.

[0035] The satellite gravity anomaly data is corrected according to the Moho correction value and the seawater correction value to obtain the sedimentary layer Bouguer gravity anomaly data;

[0036] Under the geometric constraint of the basement interface, the Bouguer gravity anomaly data of the sedimentary layer is separated into potential fields by wavelet transform to obtain the gravity anomaly data of the basement interface.

[0037] According to the geometric constraint information of the basement interface, the initial basement depth estimation pbass0 is obtained by using the two-dimensional interpolation method;

[0038] According to the gravity anomaly data of the basement interface, the basement interface depth constraint information is introduced in a hard constraint manner, and based on the lateral variable density constraint information of the basement interface, the depth of the basement interface is iteratively calculated. The formula is:

[0039] pbassk=pbassk-1+JTJ+μRTR-1ATgbass0-gbassk-1-μRTRpbassk-1,

[0040] in, p b ass k and p b ass k-1 are the results of the k-th and k-1-th substrate depth iterations, respectively. J is the Jacobian matrix, J=2πGρbassI (where I is the unit matrix), ρbass is the substrate lateral variable density information, G is the gravity constant, R is the first-order derivative operator, which controls the smoothness of the depth solution; μ is the regularization parameter, which controls the relationship between the model term and the data fitting term.

[0041] Furthermore, according to the depth of the main density interface and the density constraint information, a three-dimensional density reference model and upper and lower limit constraint models that conform to the prior constraint information are established, including:

[0042] According to the seabed bathymetric data, the depth of the basement interface and the depth of the Moho surface, the underground three-dimensional medium is divided into four layers, and a three-dimensional geometric constraint model is established. The four layers include seawater layer, sedimentary layer, crust and mantle;

[0043] In the three-dimensional geometric constraint model, the average density corresponding to the four layers is added to construct a three-dimensional density reference model;

[0044] In the three-dimensional geometric constraint model, the maximum densities corresponding to the four layers are added to construct a three-dimensional density upper limit model;

[0045] In the three-dimensional geometric constraint model, the minimum densities corresponding to the four layers are added to construct a three-dimensional density lower limit model.

[0046] Furthermore, based on satellite gravity anomalies, according to the three-dimensional density reference model and the upper and lower limit constraint models, a regularized inversion method is used to perform three-dimensional constrained inversion to construct a three-dimensional density model of the crust, including:

[0047] According to the three-dimensional density reference model and the three-dimensional density upper and lower limit constraint model, the satellite gravity anomaly data is fitted by iterative inversion, and a series of three-dimensional density inversion results are obtained by setting different fitting errors and smoothness. The formula is:

[0048]

[0049] Among them, A is the sensitivity matrix, d is the satellite gravity anomaly data, m is the three-dimensional density inversion result, m0 is the three-dimensional density reference model, Wd is the data weight matrix; Wm is the model weight matrix; μ is the regularization parameter; the logarithmic barrier method is used to constrain the range of m:

[0050]

[0051] Among them, mmax and mmin are the upper limit model and lower limit model of three-dimensional density respectively, and λ is the weight factor;

[0052] By comparing the similarity between the 3D density inversion result and the 2D reflection seismic profile, a similarity comparison result is obtained, and the result with the highest similarity is taken as the final 3D density inversion result;

[0053] The final three-dimensional density inversion results are visualized using three-dimensional visualization software to obtain a three-dimensional density model.

[0054] In a second aspect, a system for three-dimensional crustal density modeling based on seismic constraints comprises:

[0055] An acquisition module is used to acquire a two-dimensional reflection seismic profile and a two-dimensional refraction seismic velocity model in the study area, and construct constraint information of a main density layer, wherein the constraint information of the main density layer includes geometric constraint information and density constraint information;

[0056] The processing module is used to extract anomalies and perform depth inversion on the main density interface based on the satellite gravity anomaly according to the geometric constraint information and the density constraint information to obtain the depth of the main density interface; establish a three-dimensional density reference model and upper and lower limit constraint models that meet the prior constraint information according to the depth of the main density interface and the density constraint information; based on the satellite gravity anomaly, according to the three-dimensional density reference model and the upper and lower limit constraint models, use a regularized inversion method to perform three-dimensional constraint inversion to construct a three-dimensional density model of the crust.

[0057] According to a third aspect, a computing device includes:

[0058] one or more processors;

[0059] The storage system is used to store one or more programs. When the one or more programs are executed by the one or more processors, the one or more processors implement the above method.

[0060] In a fourth aspect, a computer-readable storage medium is provided, wherein a program is stored in the computer-readable storage medium, and when the program is executed by a processor, the method according to any one of claims 1 to 7 is implemented.

[0061] The above solution of the present invention includes at least the following beneficial effects:

[0062] A scheme for introducing two-dimensional constraint information into three-dimensional physical property inversion is proposed, namely: first, geometric constraints and physical property constraint information are constructed based on the two-dimensional seismic profile results, and the density interface inversion process is introduced to obtain density interface inversion results that conform to the prior information; secondly, a three-dimensional density reference model and density upper and lower limit models are constructed based on the density interface inversion results and the physical property constraint information, and a three-dimensional density inversion process is introduced to obtain three-dimensional density inversion results of the crust that conform to the prior information, thereby improving the reliability of the inversion results. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 It is a flow chart of a method for three-dimensional crust density modeling based on earthquake constraints provided by an embodiment of the present invention;

[0064] Figure 2 is a schematic diagram of a system for three-dimensional crustal density modeling based on earthquake constraints provided by an embodiment of the present invention;

[0065] Figure 3 It is a schematic diagram of constructing a geometric constraint model (d) by using seabed bathymetry (a), basement depth (b), and Moho depth (c) in a method for three-dimensional crust density modeling based on seismic constraints provided by an embodiment of the present invention;

[0066] Figure 4The method for three-dimensional density modeling of the earth's crust based on earthquake constraints provided by an embodiment of the present invention adds density constraint parameters to three-dimensional geometric constraints to obtain a three-dimensional density reference model (a), a three-dimensional density upper limit model (b) and a three-dimensional density lower limit model (c);

[0067] Figure 5 The method for three-dimensional density modeling of the crust based on earthquake constraints provided by an embodiment of the present invention utilizes satellite gravity anomalies to invert a three-dimensional density reference model, a three-dimensional density upper limit model and a three-dimensional density lower limit model to obtain a schematic diagram of a three-dimensional density model. DETAILED DESCRIPTION

[0068] The exemplary embodiments of the present disclosure will be described in more detail below with reference to the accompanying drawings. Although the exemplary embodiments of the present disclosure are shown in the accompanying drawings, it should be understood that the present disclosure can be implemented in various forms and should not be limited by the embodiments set forth herein. On the contrary, these embodiments are provided to enable a more thorough understanding of the present disclosure and to fully convey the scope of the present disclosure to those skilled in the art.

[0069] like Figures 1 to 5 As shown, an embodiment of the present invention provides a method for three-dimensional density modeling based on seismic constraints, the method comprising:

[0070] Step 11, obtaining a two-dimensional reflection seismic profile and a two-dimensional refraction seismic velocity model in the study area, and constructing constraint information of main density layers, wherein the constraint information of the main density layers includes geometric constraint information and density constraint information;

[0071] Step 12, based on the satellite gravity anomaly, performing anomaly extraction and depth inversion on the main density interface according to the geometric constraint information and the density constraint information to obtain the depth of the main density interface;

[0072] Step 13, establishing a three-dimensional density reference model and upper and lower limit constraint models that conform to the prior constraint information according to the depth of the main density interface and the density constraint information;

[0073] Step 14, based on satellite gravity anomalies, according to the three-dimensional density reference model and the upper and lower limit constraint models, a regularized inversion method is used to perform three-dimensional constrained inversion to construct a three-dimensional density model of the crust.

[0074] In the embodiment of the present invention, by identifying, digitizing and discretizing the main density interfaces in the obtained two-dimensional reflection seismic profile results and the two-dimensional refraction seismic velocity model to establish geometric constraint information, the scattered depth constraints of the main density interfaces such as the basement and Moho surface can be obtained; by extracting the seismic longitudinal wave velocities in different tectonic regions from the two-dimensional refraction seismic velocity model and converting them, the density constraint information of the main density layer and the lateral variable density constraint information of the main interface can be established; by introducing the Moho surface geometric constraint information and the Moho surface lateral variable density constraint information into the Moho surface inversion, the Moho surface depth that better meets the actual constraint conditions can be obtained, providing a basis for subsequent modeling; by introducing the basement geometric constraint information and the Moho surface lateral variable density constraint information into the Moho surface inversion, the Moho surface depth that better meets the actual constraint conditions can be obtained, providing a basis for subsequent modeling; By introducing the lateral variable density constraint information of the basement into the basement inversion, a basement depth that is more in line with the actual constraint conditions can be obtained, providing a basis for subsequent modeling; by combining the acquired seabed bathymetric data, the basement and Moho surface depths obtained by inversion, and the density constraint information, a three-dimensional reference model and a density upper and lower limit constraint model that conform to the prior constraint information can be established, so that multiple constraint information can be integrated into three-dimensional constraint information; by introducing the three-dimensional reference model and the density upper and lower limit constraint model into the three-dimensional density inversion, and using the regularized inversion method to fit the satellite gravity anomaly, a three-dimensional density model that is more in line with the prior constraint information can be constructed, so as to achieve effective prediction of underground density distribution and provide an important reference for geological exploration and resource assessment.

[0075] like Figures 1 to 5 As shown, step 11 is to obtain a two-dimensional reflection seismic profile and a two-dimensional refraction seismic velocity model in the study area, and construct constraint information of the main density layers, including:

[0076] Step 1111, obtaining a two-dimensional reflection seismic profile and a two-dimensional refraction seismic velocity model in the study area;

[0077] Step 1112, interpreting the stratigraphic structure of the two-dimensional reflection seismic profile and the two-dimensional refraction seismic velocity model, identifying the main interfaces, so as to obtain the basement interface between the sedimentary layer and the crust and the Moho surface between the crust and the mantle;

[0078] Step 1113, discretize and digitize the depth of the basement interface and the Moho surface, and establish basement interface geometric constraint information and Moho surface geometric constraint information.

[0079] In an embodiment of the present invention, by interpreting a two-dimensional reflection seismic profile and a two-dimensional refraction seismic velocity model, different stratigraphic interfaces are identified, and geometric information of underground stratigraphic interfaces, including the depth and shape of the interfaces, can be provided; according to the interpretation results, digitization and discretization are performed to obtain scattered depth information of the main interfaces, and geometric constraint information can be established; the geometric constraint information can be used in subsequent interface inversion to constrain the extraction of basement and Moho gravity anomalies and control the inversion process of the basement and Moho, which can improve the accuracy of the basement and Moho inversion results; establishing geometric constraint information is one of the key steps in three-dimensional density modeling based on seismic constraints; converting two-dimensional seismic structure prior information into geometric constraint information and introducing it into three-dimensional density modeling of the crust, so that the modeling results are more in line with the actual geological conditions; the establishment of geometric constraint information provides a basis for the subsequent construction of a three-dimensional density reference model, upper and lower density limit models, and three-dimensional constrained gravity inversion.

[0080] like Figures 1 to 5 As shown, step 11, obtaining a two-dimensional reflection seismic profile and a two-dimensional refraction seismic velocity model in the study area, and constructing constraint information of the main density layer, also includes:

[0081] Step 1121, extracting seismic P-wave velocity information of sedimentary layers, crust and mantle in different tectonic regions based on the two-dimensional refraction seismic velocity model;

[0082] Step 1122, according to the velocity-density conversion relationship, convert the seismic P-wave velocity information into density parameters, and establish density constraint information of sedimentary layers, crust and mantle in different regions, wherein the density constraint information includes maximum value, minimum value and average value;

[0083] Step 1123, based on the density constraint information of the sedimentary layer, the crust and the mantle, establish the lateral variable density constraint information of the basement interface and the lateral variable density constraint information of the Moho surface.

[0084] In the embodiment of the present invention, by digitizing and discretizing the two-dimensional refraction seismic velocity model, the seismic P-wave velocity information of the sedimentary layer, crust and mantle in different tectonic regions can be extracted; according to the velocity-density conversion relationship, the seismic P-wave velocity can be converted into a density parameter to obtain the maximum, minimum and average density of the sedimentary layer, crust and mantle in different regions, so as to establish the density constraint information of the main density layer; according to the density constraint information of the sedimentary layer, crust and mantle in different regions, the lateral variable density constraint information of the basement and Moho surface can be established; the lateral variable density constraint information of the basement and Moho surface can be used for the depth inversion of the basement and Moho surface to improve the reliability of the inversion result; the two-dimensional seismic velocity prior information is converted into density constraint information and introduced into the three-dimensional density modeling of the crust, so that the modeling result is more in line with the actual geological conditions; establishing density constraint information is one of the key steps of three-dimensional density modeling based on seismic constraints; the establishment of density constraint information provides a basis for subsequent interface inversion and three-dimensional constrained gravity inversion.

[0085] like Figures 1 to 5 As shown, step 12, based on the satellite gravity anomaly, according to the geometric constraint information and the density constraint information, anomaly extraction and depth inversion of the main density interface are performed to obtain the depth of the main density interface, including:

[0086] Step 1211, obtaining satellite gravity anomaly data, seafloor bathymetry data, and sediment thickness data;

[0087] Step 1212, based on the seabed bathymetric data, the seawater layer is divided into a first right hexahedron series, and a spatial domain method is used for forward modeling to obtain seawater correction values;

[0088] Step 1213, according to the sediment thickness data, the sediment layer is divided into a second right hexahedron series, and forward modeling is performed based on the compaction model using the space domain method to obtain the sediment layer correction value. The density at different depths is:

[0089]

[0090] Among them, g0 is the gravity anomaly of the Moho surface, ρz is the density of the sediment layer at depth z, ρf is the density of the fluid in the sediment, ρg is the density of the rock particles, is the initial porosity, c is the porosity depth attenuation coefficient;

[0091] Step 1214, correcting the satellite gravity anomaly data according to the seawater correction value and the sediment layer correction value to obtain the crustal Bouguer gravity anomaly data;

[0092] Step 1215, based on the geometric constraint information of the Moho surface, perform potential field separation on the crustal Bouguer gravity anomaly data to obtain the Moho surface gravity anomaly data;

[0093] Step 1216, based on the Moho geometric constraint information, use the linear regression method to obtain an initial Moho depth estimate pmo h o 0 ;

[0094] Step 1217, according to the Moho surface gravity anomaly data, the Moho surface geometric constraint information is introduced in a soft constraint manner, and based on the Moho surface lateral variable density constraint information, the depth of the Moho surface is iteratively calculated, and the formula is:

[0095] pmohok=pmohok-1+gmoho0-gmohok-12πGρmoho+Δpmoho,

[0096] Wherein, gmoho0 is the Moho gravity anomaly, gmohok-1 is the forward gravity response of the k-1th Moho depth, pmohok and pmohok-1 are the results of the k-th and k-1th iterations respectively, G is the gravity constant, ρmoho is the lateral variable density information of the Moho, and Δpmoho is the depth correction obtained by linear regression based on the geometric constraint information of the Moho for the kth time.

[0097] In an embodiment of the present invention, forward modeling is performed based on seabed bathymetric data to obtain seawater correction values. First, the seabed bathymetric data is used to determine the shape and height of the ocean bottom. Based on the forward modeling method, the gravitational effect of the ocean bottom is simulated. Forward modeling refers to simulating the propagation and response process of the geophysical field by numerical calculation. According to the forward modeling calculation results, the gravitational effect of the ocean bottom, that is, the seawater correction value, is obtained. The seawater correction value represents the contribution of the ocean to the gravity anomaly. According to the sediment thickness data, forward modeling is performed based on the compaction model to obtain the sediment layer correction value. First, the sediment thickness data is used to determine the thickness distribution of the sediment layer in the crust. Then, based on the compaction model and the forward modeling method, the sediment layer is simulated. Gravitational effect, compaction model refers to a model that calculates the gravitational effect of the sedimentary layer based on its physical parameters, such as density and thickness. According to the forward calculation results, the gravitational effect of the sedimentary layer, that is, the correction value of the sedimentary layer, is obtained. The sedimentary layer correction value represents the contribution of the sedimentary layer to the gravity anomaly; the satellite gravity anomaly data is corrected according to the seawater correction value and the sedimentary layer correction value to obtain the crustal Bouguer gravity anomaly data. First, the satellite gravity anomaly data is used to subtract the seawater correction value and the sedimentary layer correction value to obtain the gravity anomaly data caused by other factors in the crust (such as crustal density changes, tectonic deformation, etc.). The gravity anomaly data obtained is the crustal Bouguer gravity anomaly data, which reflects the effects of seawater and sedimentary layers on the crust. The influence of other factors in the crust outside the Moho surface; Based on the geometric constraint information of the Moho surface, the potential field separation of the crustal Bouguer gravity anomaly data is carried out to obtain the Moho surface gravity anomaly data. First, the Moho surface geometric constraint information is used to determine the depth and shape of the Moho surface. Then, according to the potential field separation method, the crustal Bouguer gravity anomaly data is decomposed into Moho surface gravity anomaly data and residual crustal gravity anomaly data; Moho surface gravity anomaly data reflects the gravity anomaly caused by the change of the Moho surface in the crust, which can be used to study the structure and tectonic structure of the crust; According to the geometric constraint information of the Moho surface, the linear regression method can be used to make an initial estimate of the depth of the Moho surface more reliably to obtain the average depth of the Moho surface. The inversion process of the Moho surface is iterative; the geometric constraint information of the Moho surface is introduced into the interface inversion process by means of soft constraints, which helps to eliminate the fuzzy imaging of the Moho surface by two-dimensional reflection seismic; based on the lateral variable density constraint, the interface distribution is obtained by iterative optimization method, making the result more consistent with the actual geological situation; the depth of the Moho surface is an important parameter of the underground structure and can be used to construct a three-dimensional density reference model and a density upper and lower limit model; the depth inversion of the Moho surface under geometric constraints and lateral variable density constraints is one of the key steps in the three-dimensional density modeling of the crust based on seismic constraints; by obtaining the depth of the Moho surface, it can provide a basis for the subsequent three-dimensional constrained gravity inversion and the construction of a three-dimensional density model.

[0098] like Figures 1 to 5As shown, step 12, based on the satellite gravity anomaly, according to the geometric constraint information and the density constraint information, anomaly extraction and depth inversion of the main density interface are performed to obtain the depth of the main density interface, and also includes:

[0099] Step 1221, performing forward modeling according to the depth of the Moho surface to obtain a Moho surface correction value;

[0100] Step 1222, performing correlation analysis based on satellite gravity anomaly data and seafloor bathymetry data to obtain optimal density parameters for seawater correction, and using the seafloor bathymetry data to divide the seawater layer into a series of right hexahedrons, and using the spatial domain method to forward calculate and obtain the seawater correction value;

[0101] Step 1223, correcting the satellite gravity anomaly data according to the Moho correction value and the seawater correction value to obtain the sedimentary layer Bouguer gravity anomaly data;

[0102] Step 1224, under the geometric constraint of the basement interface, the Bouguer gravity anomaly data of the sedimentary layer is subjected to potential field separation by using a wavelet transform method to obtain gravity anomaly data of the basement interface;

[0103] Step 1225, based on the geometric constraint information of the substrate interface, a two-dimensional interpolation method is used to obtain an initial substrate depth estimate. p b ass 0 ;

[0104] Step 1226, based on the gravity anomaly data of the basement interface, the basement interface depth constraint information is introduced in a hard constraint manner, and based on the lateral variable density constraint information of the basement interface, the depth of the basement interface is iteratively calculated, and the formula is:

[0105] pbassk=pbassk-1+JTJ+μRTR-1ATgbass0-gbassk-1-μRTRpbassk-1,

[0106] in, p b ass k and p b ass k-1 are the results of the k-th and k-1-th substrate depth iterations, respectively. J is the Jacobian matrix, J=2πGρbassI (where I is the unit matrix), ρbass is the substrate lateral variable density information, G is the gravity constant, R is the first-order derivative operator, which controls the smoothness of the depth solution; μ is the regularization parameter, which controls the relationship between the model term and the data fitting term.

[0107] In the embodiment of the present invention, forward modeling is performed according to the depth of the Moho surface to obtain the Moho surface correction value. First, the depth data of the Moho surface is obtained to determine the shape and position of the Moho interface. Then, based on the forward modeling method, the influence of the Moho interface on the gravity anomaly is simulated. Forward modeling refers to simulating the propagation and response process of the geophysical field by numerical calculation. According to the forward modeling result, the gravitational effect of the Moho surface, that is, the correction value of the Moho surface, is obtained. The Moho surface correction value represents the contribution of the Moho interface to the gravity anomaly; the seabed bathymetric data is used to divide the seawater layer into a series of right hexahedrons, and the position and size of each hexahedron are determined; the spatial domain method is used to perform forward modeling, and the seawater correction value is calculated according to the seawater density parameter and the distribution of the seawater layer; the calculated seawater correction value is calculated. The seawater correction value is analyzed and verified to ensure the reliability and accuracy of the results; the final seawater correction value can be used to correct the satellite gravity anomaly data to improve the accuracy and reliability of the data; firstly, the satellite gravity anomaly data and the seabed bathymetric data are correlated to determine the optimal density parameter for seawater correction. The correlation analysis is used to determine the degree of correlation between the two data sets. The forward modeling method is used to simulate the contribution of seawater to the gravity anomaly using the seabed bathymetric data and the optimal density parameter. The forward calculation is performed using an analytical method. According to the forward calculation results, the gravitational effect of seawater is obtained, that is, the seawater correction value. The seawater correction value represents the correction of seawater to the gravity anomaly; the satellite gravity anomaly data is corrected according to the Moho correction value and the seawater correction value. Correction is performed to obtain the Bouguer gravity anomaly data of the sedimentary layer. First, the satellite gravity anomaly data is corrected using the Moho correction value and the seawater correction value. The Moho correction value and the seawater correction value are subtracted from the original data to obtain the Bouguer gravity anomaly data of the sedimentary layer. The gravity anomaly data obtained reflects the contribution of other factors to the gravity anomaly except the Moho interface and seawater, which mainly comes from the sedimentary layer in the crust; the Bouguer gravity anomaly data of the sedimentary layer is processed by the wavelet transform method, and it is decomposed into wavelet coefficients of different scales and frequencies for subsequent analysis; according to the geometric constraints of the basement interface, the appropriate wavelet coefficients are selected for inverse transformation to obtain the gravity anomaly data of the basement interface; the obtained basement interface gravity anomaly data is analyzed and Verification ensures that the results meet expectations, and makes necessary corrections and adjustments; the final basement interface gravity anomaly data can be used for further geological interpretation and research, providing important information for geological structures and sedimentary processes; the basement interface gravity anomaly data reflects the gravity anomaly caused by changes in the crustal basement, which is of great significance for studying the structure and geological characteristics of the crust; the use of hard constraints to introduce basement geometry constraints into the basement interface inversion process can improve the reliability of the basement inversion results; based on the basement lateral variable density constraints, the basement interface distribution is obtained through iterative optimization methods, making the results more consistent with the actual geological conditions; the depth of the basement is an important parameter of the underground structure, which can be used to construct a three-dimensional density reference model and a density upper and lower limit model;The basement interface inversion under geometric constraints and lateral variable density constraints is one of the key steps in the three-dimensional density modeling of the crust based on seismic constraints; by obtaining the depth of the basement, it can provide a basis for the subsequent three-dimensional constrained gravity inversion and the construction of a three-dimensional density model. ;

[0108] like Figures 1 to 5 As shown, step 13, according to the depth of the main density interface and the density constraint information, establish a three-dimensional density reference model and upper and lower limit constraint models that conform to the prior constraint information, including:

[0109] Step 131, dividing the underground three-dimensional medium into four layers according to the seabed bathymetric data, the depth of the basement interface and the depth of the Moho surface, and establishing a three-dimensional geometric constraint model, wherein the four layers include a seawater layer, a sedimentary layer, a crust and a mantle;

[0110] Step 132, adding the average densities corresponding to the four layers to the three-dimensional geometric constraint model to construct a three-dimensional density reference model;

[0111] Step 133, adding the maximum densities corresponding to the four layers to the three-dimensional geometric constraint model to construct a three-dimensional density upper limit model;

[0112] Step 134, adding the minimum densities corresponding to the four layers respectively into the three-dimensional geometric constraint model to construct a three-dimensional density lower limit model.

[0113] In an embodiment of the present invention, according to the acquired seabed bathymetric data, the basement and Moho depths obtained by inversion, the underground three-dimensional medium can be divided into four layers, namely, seawater layer, sedimentary layer, crust and mantle, to establish a three-dimensional geometric constraint model; in the three-dimensional geometric constraint model, the average density, maximum density and minimum density of the four layers are added to construct a three-dimensional density reference model, a density upper limit model and a density lower limit model respectively; the three-dimensional density reference model is the basis of the three-dimensional constraint inversion, which provides the approximate distribution and characteristics of the underground structure and provides a reference for the subsequent density inversion; the physical property upper and lower limit constraint model is used to constrain the physical property parameter range of the density model, so that the density model is more consistent with the actual geological conditions; the three-dimensional reference model and the density upper and lower limit constraint model that conform to the prior constraint information can improve the reliability of modeling; establishing a three-dimensional reference model and a physical property upper and lower limit constraint model that conform to the prior constraint information is one of the key steps of three-dimensional density modeling based on seismic constraints; the reference model and the physical property upper and lower limit constraint model provide the basis for the subsequent regularized inversion.

[0114] like Figures 1 to 5 As shown, step 14, based on satellite gravity anomalies, according to the three-dimensional density reference model and the upper and lower limit constraint models, a regularized inversion method is used to perform three-dimensional constrained inversion to construct a three-dimensional density model of the crust, including:

[0115] Step 141, according to the three-dimensional density reference model and the three-dimensional density upper and lower limit constraint model, the satellite gravity anomaly data is fitted by iterative inversion, and a series of three-dimensional density inversion results are obtained by setting different fitting errors and smoothness. The formula is:

[0116]

[0117] Among them, A is the sensitivity matrix, d is the satellite gravity anomaly data, m is the three-dimensional density inversion result, m0 is the three-dimensional density reference model, Wd is the data weight matrix; Wm is the model weight matrix; μ is the regularization parameter; the logarithmic barrier method is used to constrain the range of m:

[0118]

[0119] Among them, mmax and mmin are the upper limit model and lower limit model of three-dimensional density respectively, and λ is the weight factor;

[0120] Step 142, by comparing the similarity of the three-dimensional density inversion result with the two-dimensional reflection seismic profile, to obtain a similarity comparison result, and taking the result with high similarity as the final three-dimensional density inversion result;

[0121] Step 143, using three-dimensional visualization software, visualize the final three-dimensional density inversion result to obtain a three-dimensional density model.

[0122] In an embodiment of the present invention, by introducing a three-dimensional reference model and a density upper and lower limit model into the inversion process, the reliability of the inversion result can be improved; the regularization method can balance the smoothness of the model and the data fitting degree, and can make the inversion result smoother; the two-dimensional reflection seismic profile is compared with the three-dimensional slice, the similarity in the underground structure is observed, and the characteristics such as the stratigraphic interface and lithological changes are compared. According to the comparison results, the one with the highest similarity is selected as the three-dimensional density inversion result, so that the inversion result can meet the prior constraint information; using three-dimensional visualization software, the three-dimensional density inversion result can be visualized, and a three-dimensional density model of the crust can be established; through visualization, the underground density distribution can be intuitively observed, which helps geologists and geophysicists better understand the distribution of underground density.

[0123] Optionally, in one embodiment of the present invention, density slices at different depths and thicknesses of different layers are extracted according to the three-dimensional density model of the earth's crust, and favorable oil and gas favorable areas are determined based on regional oil and gas accumulation rules, including:

[0124] According to the regional oil and gas accumulation rules, determine the main hydrocarbon-generating strata and traps;

[0125] According to the two-dimensional reflection seismic profile, the top and bottom interface depths of the hydrocarbon-generating layer and the depth of the trap are extracted;

[0126] According to the 3D density model, the average density of the top and bottom interfaces of the hydrocarbon-generating layer at the corresponding position of the 2D reflection seismic section is obtained;

[0127] According to the average density of the top and bottom interfaces of the hydrocarbon-generating layer, the three-dimensional depths of the top and bottom interfaces of the hydrocarbon-generating layer are extracted from the three-dimensional density model, and the thickness of the hydrocarbon-generating layer is calculated by subtracting the interface depths.

[0128] Based on the three-dimensional density model, density slices at different depths are extracted according to the depth of the trap;

[0129] Determine the depth and location of the hydrocarbon-generating stratum uplift based on the high-density features in the density slices;

[0130] The favorable oil and gas areas are determined based on the thickness of the hydrocarbon-generating strata and the depth and position of the uplift of the hydrocarbon-generating strata.

[0131] In an embodiment of the present invention, the regional oil and gas accumulation law is studied to determine the main hydrocarbon-generating layers and traps; according to the two-dimensional reflection seismic profile, the top interface and bottom interface depths of the hydrocarbon-generating layer and the depth of the trap are extracted, firstly, the reflection waveform on the reflection seismic profile is analyzed to identify the top interface and bottom interface of the hydrocarbon-generating layer, and then the depths of the top interface and bottom interface of the hydrocarbon-generating layer are marked on the profile using seismic interpretation software or other related tools; then, according to the depths of the top interface and bottom interface of the hydrocarbon-generating layer, the depth of the trap is determined, and the trap is composed of the area between the top interface and the bottom interface of the hydrocarbon-generating layer; according to the three-dimensional density model, the average density of the top interface and the bottom interface of the hydrocarbon-generating layer at the corresponding position of the two-dimensional reflection seismic profile is obtained, firstly, the depths of the top interface and the bottom interface of the hydrocarbon-generating layer in the two-dimensional reflection seismic profile are matched with the three-dimensional density model, and in the three-dimensional density model, the average density value of the position corresponding to the depth of the top interface and the bottom interface of the hydrocarbon-generating layer is extracted, and the average density value of the top interface and the bottom interface of the hydrocarbon-generating layer is recorded; according to the average density of the top interface and the bottom interface of the hydrocarbon-generating layer, the average density of the top interface and the bottom interface of the hydrocarbon-generating layer is obtained from The three-dimensional depth of the top and bottom interfaces of the hydrocarbon-generating layer is extracted from the three-dimensional density model, and the thickness of the hydrocarbon-generating layer is calculated by subtracting the interface depths. First, the average density values ​​of the top and bottom interfaces of the hydrocarbon-generating layer are used to determine the three-dimensional depths of the top and bottom interfaces of the hydrocarbon-generating layer in the three-dimensional density model, and the thickness of the hydrocarbon-generating layer is calculated by subtracting the three-dimensional depths of the top and bottom interfaces. According to the three-dimensional density model, density slices of different depths are extracted according to the depth of the closure. First, based on the depth range of the closure, density slices of corresponding depths are extracted in the three-dimensional density model. The density slice is a density distribution image within a specific depth range. According to the high-density features in the density slices, the depth and position of the uplift of the hydrocarbon-generating layer are determined. In the density slices, the high-density features are observed, and the depth and position of the uplift of the hydrocarbon-generating layer are determined by determining the depth and position of the high-density features. Combined with the thickness of the hydrocarbon-generating layer and the depth and position of the uplift, the range of the favorable oil and gas area is determined; according to the geological conditions and sedimentary environment, the influence of the thickness and uplift of the hydrocarbon-generating layer on the accumulation of oil and gas is considered to further determine the favorable oil and gas area.

[0132] like Figure 2 As shown, an embodiment of the present invention further provides a system 20 for three-dimensional crust density modeling based on seismic constraints, comprising:

[0133] An acquisition module 21 is used to acquire a two-dimensional reflection seismic profile and a two-dimensional refraction seismic velocity model in the study area, and construct constraint information of a main density layer, wherein the constraint information of the main density layer includes geometric constraint information and density constraint information;

[0134] The processing module 22 is used to extract anomalies and perform depth inversion on the main density interface based on the satellite gravity anomaly and according to the geometric constraint information and the density constraint information to obtain the depth of the main density interface; establish a three-dimensional density reference model and upper and lower limit constraint models that conform to the prior constraint information based on the depth of the main density interface and the density constraint information; based on the satellite gravity anomaly, according to the three-dimensional density reference model and the upper and lower limit constraint model, use the regularized inversion method to perform three-dimensional constraint inversion to construct a three-dimensional density model of the crust.

[0135] The embodiment of the present invention further provides a computing device, comprising: a processor, a memory storing a computer program, wherein when the computer program is executed by the processor, the method described above is executed. All implementations in the above method embodiment are applicable to this embodiment and can achieve the same technical effect.

[0136] The embodiment of the present invention also provides a computer-readable storage medium storing instructions, which, when executed on a computer, enable the computer to execute the method described above. All implementations in the above method embodiment are applicable to this embodiment and can achieve the same technical effect.

[0137] Those skilled in the art will appreciate that the units and algorithm steps of each example described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are performed in hardware or software depends on the specific application and design constraints of the technical solution. Professional and technical personnel can use different methods to implement the described functions for each specific application, but such implementation should not be considered to be beyond the scope of the present invention.

[0138] Those skilled in the art can clearly understand that, for the convenience and brevity of description, the specific working processes of the above-described systems, systems and units can refer to the corresponding processes in the aforementioned method embodiments and will not be repeated here.

[0139] In the embodiments provided by the present invention, it should be understood that the disclosed systems and methods can be implemented in other ways. For example, the system embodiments described above are only schematic. For example, the division of the units is only a logical function division. There may be other division methods in actual implementation, such as multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed. Another point is that the mutual coupling or direct coupling or communication connection shown or discussed can be through some interfaces, indirect coupling or communication connection of systems or units, which can be electrical, mechanical or other forms.

[0140] The units described as separate components may or may not be physically separated, and the components shown as units may or may not be physical units, that is, they may be located in one place or distributed on multiple network units. Some or all of the units may be selected according to actual needs to achieve the purpose of the solution of this embodiment.

[0141] In addition, each functional unit in each embodiment of the present invention may be integrated into one processing unit, or each unit may exist physically separately, or two or more units may be integrated into one unit.

[0142] If the functions are implemented in the form of software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, or the part that contributes to the prior art, or the part of the technical solution, can be embodied in the form of a software product. The computer software product is stored in a storage medium, including several instructions for a computer device (which can be a personal computer, server, or network device, etc.) to perform all or part of the steps of the method described in each embodiment of the present invention. The aforementioned storage medium includes: various media that can store program codes, such as USB flash drives, mobile hard disks, ROM, RAM, magnetic disks, or optical disks.

[0143] In addition, it should be noted that in the system and method of the present invention, it is obvious that each component or each step can be decomposed and / or recombined. These decompositions and / or recombinations should be regarded as equivalent schemes of the present invention. Moreover, the steps of performing the above-mentioned series of processing can naturally be performed in chronological order according to the order of description, but it is not necessary to perform them in chronological order, and some steps can be performed in parallel or independently of each other. For those of ordinary skill in the art, it is understood that all or any steps or components of the method and system of the present invention can be implemented in any computing system (including processors, storage media, etc.) or a network of computing systems in hardware, firmware, software or a combination thereof, which can be achieved by those of ordinary skill in the art using their basic programming skills after reading the description of the present invention.

[0144] Therefore, the purpose of the present invention can also be achieved by running a program or a group of programs on any computing system. The computing system can be a well-known general system. Therefore, the purpose of the present invention can also be achieved by simply providing a program product containing a program code that implements the method or system. That is to say, such a program product also constitutes the present invention, and the storage medium storing such a program product also constitutes the present invention. Obviously, the storage medium can be any well-known storage medium or any storage medium developed in the future. It should also be pointed out that in the system and method of the present invention, it is obvious that each component or each step can be decomposed and / or recombined. These decompositions and / or recombinations should be regarded as equivalent schemes of the present invention. In addition, the steps of performing the above-mentioned series of processing can naturally be performed in chronological order according to the order of description, but it is not necessary to perform them in chronological order. Some steps can be performed in parallel or independently of each other.

[0145] The above is a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A method for three-dimensional crustal density modeling based on earthquake constraints, characterized in that: The method comprises: Obtaining a two-dimensional reflection seismic profile and a two-dimensional refraction seismic velocity model in the study area, and constructing constraint information of main density layers, wherein the constraint information of the main density layers includes geometric constraint information and density constraint information; Based on the satellite gravity anomaly, the main density interface is extracted and depth inverted according to the geometric constraint information and the density constraint information to obtain the depth of the main density interface; According to the depth of the main density interface and the density constraint information, a three-dimensional density reference model and upper and lower limit constraint models that conform to the prior constraint information are established; Based on satellite gravity anomalies, according to the three-dimensional density reference model and the upper and lower limit constraint models, the regularized inversion method is used to perform three-dimensional constrained inversion and construct a three-dimensional density model of the crust.

2. The method for three-dimensional crustal density modeling based on earthquake constraints according to claim 1, characterized in that: Obtain 2D reflection seismic profiles and 2D refraction seismic velocity models in the study area and construct constraint information for the main density horizons, including: Obtain 2D reflection seismic profiles and 2D refraction seismic velocity models in the study area; Interpret the stratigraphic structure of 2D reflection seismic profiles and 2D refraction seismic velocity models to identify the main interfaces, such as the basement interface between sedimentary layers and the crust and the Moho between the crust and the mantle; The depths of the basement interface and the Moho surface are discretized and digitized, and the geometric constraint information of the basement interface and the Moho surface is established.

3. The method for crustal three-dimensional density modeling based on earthquake constraints according to claim 2, characterized in that: Obtain 2D reflection seismic profiles and 2D refraction seismic velocity models in the study area, construct constraint information for the main density horizons, and also include: Extract seismic P-wave velocity information of sedimentary layers, crust and mantle in different tectonic regions based on the two-dimensional refraction seismic velocity model; According to the velocity-density conversion relationship, seismic P-wave velocity information is converted into density parameters, and density constraint information of sedimentary layers, crust and mantle in different regions is established, wherein the density constraint information includes maximum value, minimum value and average value; According to the density constraint information of sedimentary layer, crust and mantle, the lateral variable density constraint information of basement interface and the lateral variable density constraint information of Moho surface are established.

4. The method for crustal three-dimensional density modeling based on earthquake constraints according to claim 3, characterized in that: Based on the satellite gravity anomaly, the main density interface is extracted and depth inverted according to the geometric constraint information and the density constraint information to obtain the depth of the main density interface, including: Obtain satellite gravity anomaly data, seafloor bathymetry data, and sediment thickness data; According to the seabed bathymetric data, the seawater layer is divided into a series of right hexahedrons, and the spatial domain method is used for forward modeling to obtain the seawater correction value; According to the sediment thickness data, the sedimentary layer is divided into a series of right hexahedrons. Based on the compaction model, the spatial domain method is used for forward modeling to obtain the sedimentary layer correction value. The density at different depths is: Among them, g0 is the gravity anomaly of the Moho surface, ρz is the density of the sediment layer at depth z, ρf is the density of the fluid in the sediment, ρg is the density of the rock particles, is the initial porosity, c is the porosity depth attenuation coefficient; The satellite gravity anomaly data is corrected according to the seawater correction value and the sediment layer correction value to obtain the crustal Bouguer gravity anomaly data; Based on the geometric constraint information of Moho surface, the potential field of crustal Bouguer gravity anomaly data is separated to obtain Moho surface gravity anomaly data; Based on the geometric constraint information of the Moho surface, the initial Moho surface depth estimation pmoho0 is obtained using the linear regression method; According to the Moho surface gravity anomaly data, the Moho surface geometric constraint information is introduced in a soft constraint manner. Based on the Moho surface lateral variable density constraint information, the depth of the Moho surface is iteratively calculated. The formula is: pmohok=pmohok-1+gmoho0-gmohok-12πGρmoho+Δpmoho, Where gmoho0 is the Moho gravity anomaly, gmohok-1 is the forward gravity response of the k-1th Moho depth, pmo h o k and pmo h o k-1 are the results of the kth and k-1th iterations respectively, G is the gravity constant, ρmoho is the lateral variable density information of the Moho surface, and Δpmoho is the depth correction value obtained by linear regression based on the geometric constraint information of the Moho surface for the kth time.

5. The method for three-dimensional crustal density modeling based on earthquake constraints according to claim 4, characterized in that: Based on the satellite gravity anomaly, the main density interface is extracted and depth inverted according to the geometric constraint information and the density constraint information to obtain the depth of the main density interface, and further includes: According to the depth of the Moho surface, forward modeling is performed to obtain the Moho surface correction value; Based on the satellite gravity anomaly data and the seabed bathymetry data, the correlation analysis was carried out to obtain the optimal density parameters for seawater correction. The seawater layer was divided into a series of right hexahedrons using the seabed bathymetry data, and the seawater correction value was obtained by forward calculation using the spatial domain method. The satellite gravity anomaly data is corrected according to the Moho correction value and the seawater correction value to obtain the sedimentary layer Bouguer gravity anomaly data; Under the geometric constraint of the basement interface, the Bouguer gravity anomaly data of the sedimentary layer is separated into potential fields by wavelet transform to obtain the gravity anomaly data of the basement interface. According to the geometric constraint information of the basement interface, the initial basement depth estimation pbass0 is obtained by using the two-dimensional interpolation method; According to the gravity anomaly data of the basement interface, the basement interface depth constraint information is introduced in a hard constraint manner, and based on the lateral variable density constraint information of the basement interface, the depth of the basement interface is iteratively calculated. The formula is: pbassk=pbassk-1+JTJ+μRTR-1ATgbass0-gbassk-1-μRTRpbassk-1, in, p b ass k and p b ass k-1 are the results of the k-th and k-1-th substrate depth iterations, respectively. J is the Jacobian matrix, J=2πGρbassI (where I is the unit matrix), ρbass is the substrate lateral variable density information, G is the gravity constant, R is the first-order derivative operator, which controls the smoothness of the depth solution; μ is the regularization parameter, which controls the relationship between the model term and the data fitting term.

6. The method for three-dimensional crust density modeling based on earthquake constraints according to claim 5, characterized in that: According to the depth of the main density interface and the density constraint information, a three-dimensional density reference model and upper and lower limit constraint models that conform to the prior constraint information are established, including: According to the seabed bathymetric data, the depth of the basement interface and the depth of the Moho surface, the underground three-dimensional medium is divided into four layers, and a three-dimensional geometric constraint model is established. The four layers include seawater layer, sedimentary layer, crust and mantle; In the three-dimensional geometric constraint model, the average density corresponding to the four layers is added to construct a three-dimensional density reference model; In the three-dimensional geometric constraint model, the maximum densities corresponding to the four layers are added to construct a three-dimensional density upper limit model; In the three-dimensional geometric constraint model, the minimum densities corresponding to the four layers are added to construct a three-dimensional density lower limit model.

7. The method for three-dimensional crustal density modeling based on earthquake constraints according to claim 6, characterized in that: Based on satellite gravity anomalies, according to the three-dimensional density reference model and the upper and lower limit constraint model, the regularized inversion method is used to perform three-dimensional constrained inversion and construct a three-dimensional crustal density model, including: According to the three-dimensional density reference model and the three-dimensional density upper and lower limit constraint model, the satellite gravity anomaly data is fitted by iterative inversion, and a series of three-dimensional density inversion results are obtained by setting different fitting errors and smoothness. The formula is: Among them, A is the sensitivity matrix, d is the satellite gravity anomaly data, m is the three-dimensional density inversion result, m0 is the three-dimensional density reference model, Wd is the data weight matrix; Wm is the model weight matrix; μ is the regularization parameter; the logarithmic barrier method is used to constrain the range of m: Among them, mmax and mmin are the upper limit model and lower limit model of three-dimensional density respectively, and λ is the weight factor; By comparing the similarity between the 3D density inversion result and the 2D reflection seismic profile, a similarity comparison result is obtained, and the result with the highest similarity is taken as the final 3D density inversion result; The final three-dimensional density inversion results are visualized using three-dimensional visualization software to obtain a three-dimensional density model.

8. A system for three-dimensional crustal density modeling based on earthquake constraints, characterized in that: include: An acquisition module is used to acquire a two-dimensional reflection seismic profile and a two-dimensional refraction seismic velocity model in the study area, and construct constraint information of a main density layer, wherein the constraint information of the main density layer includes geometric constraint information and density constraint information; The processing module is used to extract anomalies and perform depth inversion on the main density interface based on the satellite gravity anomaly according to the geometric constraint information and the density constraint information to obtain the depth of the main density interface; establish a three-dimensional density reference model and upper and lower limit constraint models that meet the prior constraint information according to the depth of the main density interface and the density constraint information; based on the satellite gravity anomaly, according to the three-dimensional density reference model and the upper and lower limit constraint model, use a regularized inversion method to perform three-dimensional constraint inversion to construct a three-dimensional density model of the crust.

9. A computing device, characterized in that include: one or more processors; A storage device, used for storing one or more programs, when the one or more programs are executed by the one or more processors, the one or more processors implement the method according to any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a program, which, when executed by a processor, implements the method according to any one of claims 1 to 7.

Images