A method and system for three-dimensional crustal density modeling based on seismic constraints

Abstract

The present invention provides a method and system for three-dimensional crustal density modeling based on seismic constraints, which relates to the fields of three-dimensional crustal density modeling and offshore oil and gas exploration. The method includes: obtaining two-dimensional reflection seismic profiles and two-dimensional refraction seismic velocity models in the study area, constructing constraint information of the main density horizons, and then performing anomaly extraction and depth inversion to obtain the depths of the main density interfaces; establishing a three-dimensional density reference model and upper and lower limit constraint models that conform to the prior constraint information according to the depths of the main density interfaces and the density constraint information, and performing three-dimensional constrained inversion using a regularization inversion method to construct a three-dimensional crustal density model. The present invention introduces two-dimensional constraint information into three-dimensional inversion to improve the reliability of the inversion results.

Description

Technical Field

[0001] The present invention relates to the technical fields of three-dimensional crust density modeling and oil and gas exploration, and particularly to a method and system for three-dimensional crust 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 crustal structure of the sea area 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 earth's evolution. Studying the oceanic crustal structure helps to reveal the internal mechanism of the earth's evolution.

[0003] Different marine geological units often have significant density differences, and the special geological processes of oil and gas accumulation can also cause density changes, providing a physical property basis for studying the crustal structure of the sea area and oil and gas exploration based on gravity data; satellite gravity observation technology provides high-coverage regional data, and the low-frequency components have high accuracy, having certain advantages in crust-scale research; however, due to problems such as the volume effect of the potential field, poor vertical resolution, and the non-uniqueness of geophysical inversion problems, it restricts its application in crust-scale quantitative inversion and directly determining favorable oil and gas areas.

[0004] In gravity three-dimensional 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 indirectly used to determine the best inversion parameters and evaluate the reliability of the results after inversion, and 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 seismic constraints, introducing two-dimensional constraint information into three-dimensional inversion to improve the reliability of inversion results.

[0006] 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 crust density modeling based on seismic constraints, the method comprising:

[0008] Obtain two-dimensional reflection seismic profiles and two-dimensional refraction seismic velocity models in the study area, and construct constraint information for the main density horizons, wherein the constraint information for the main density horizons includes geometric constraint information and density constraint information;

[0009] Based on satellite gravity anomalies, perform anomaly extraction and depth inversion on the main density interfaces according to the geometric constraint information and density constraint information to obtain the depths of the main density interfaces;

[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, ρ z is the sediment density at depth z, ρ f is the density of the fluid in the sediment, ρ g is the density of 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 constraints of the Moho surface, the initial Moho depth estimate 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]

[0031] in, The Moho gravity anomaly. is the forward gravity response of the k-1th Moho depth, and are the results of the kth and k-1th iterations, respectively, G is the gravity constant, ρ moho is the lateral density variation information of the Moho surface, Δ pmoho is the depth correction obtained by linear regression based on the Moho geometric constraint information 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 satellite gravity anomaly data and bathymetric data, correlation analysis is carried out to obtain the optimal density parameter for seawater correction. Using bathymetric data, the seawater layer is divided into a series of upright hexahedrons, and the seawater correction value is 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 Bouguer gravity anomaly data of the sedimentary layer;

[0036] Under the geometric constraint of the basement interface, the potential field separation of the Bouguer gravity anomaly data of the sedimentary layer is carried out by the wavelet transform method to obtain the gravity anomaly data of the basement interface;

[0037] According to the geometric constraint information of the basement interface, the two-dimensional interpolation method is used to obtain the initial estimated basement depth p bass 0 ;

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

[0039]

[0040] where and are the results of the k-th and (k - 1)-th basement depth iterations respectively. J is the Jacobian matrix, J = 2πGρ bass I (where I is the identity matrix), ρ bass is the lateral variable density information of the basement, G is the gravitational constant, R is the first 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 bathymetric data, the depth of the basement interface and the depth of the Moho surface, the underground three-dimensional medium is divided into four horizons, and a three-dimensional geometric constraint model is established. The four horizons include the seawater layer, the sedimentary layer, the crust, and the mantle;

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

[0044] In the three-dimensional geometric constraint model, the maximum density corresponding to each of the four horizons is added to construct a three-dimensional density upper limit model;

[0045] In the three-dimensional geometric constraint model, the minimum density corresponding to four horizons is added to construct a three-dimensional lower density 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 three-dimensional constrained inversion is carried out using the regularization inversion method to construct a three-dimensional crustal density model, including:

[0047] According to the three-dimensional density reference model and the three-dimensional density upper and lower limit constraint models, the satellite gravity anomaly data is fitted by means of 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, and W d is the data weight matrix; W m is the model weight matrix; μ is the regularization parameter; the range of m is constrained by the logarithmic barrier method:

[0050]

[0051] Among them, m max and m min are the three-dimensional density upper limit model and the lower limit model respectively, and λ is the weight factor;

[0052] By comparing the similarity between the three-dimensional density inversion result and the two-dimensional reflection seismic profile, the similarity comparison result is obtained, and the result with high similarity is used as the final three-dimensional density inversion result;

[0053] Using three-dimensional visualization software, the final three-dimensional density inversion result is visualized to obtain a three-dimensional density model.

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

[0055] An acquisition module for acquiring two-dimensional reflection seismic profiles and two-dimensional refraction seismic velocity models in the study area and constructing constraint information of the main density horizons, where the constraint information of the main density horizons includes geometric constraint information and density constraint information;

[0056] A processing module, configured to perform anomaly extraction and depth inversion on a main density interface based on satellite gravity anomalies according to the geometric constraint information and density constraint information, so as 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 according to the depth of the main density interface and the density constraint information; and perform three-dimensional constrained inversion by using a regularization inversion method based on satellite gravity anomalies according to the three-dimensional density reference model and the upper and lower limit constraint models to construct a three-dimensional crust density model.

[0057] In a third aspect, a computing device includes:

[0058] One or more processors;

[0059] A storage system for storing one or more programs, which, when executed by the one or more processors, cause the one or more processors to implement the above method.

[0060] In a fourth aspect, a computer-readable storage medium stores a program, which, when executed by a processor, implements the above method.

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

[0062] A solution of introducing two-dimensional constraint information into three-dimensional physical property inversion is proposed, that is: first, construct geometric constraint and physical property constraint information according to the results of two-dimensional seismic profiles, introduce the density interface inversion process to obtain the density interface inversion result that conforms to the prior information; secondly, construct a three-dimensional density reference model and density upper and lower limit models according to the density interface inversion result and the physical property constraint information, introduce the three-dimensional density inversion process to obtain the three-dimensional crust density inversion result that conforms to the prior information, and improve the reliability of the inversion result. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 is a schematic flowchart of a method for three-dimensional crust density modeling based on seismic constraints provided by an embodiment of the present invention;

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

[0065] Figure 3 is a schematic diagram of constructing a geometric constraint model (d) by using 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 4Schematic diagram of the method for three-dimensional crustal density modeling based on seismic constraints provided by an embodiment of the present invention, where density constraint parameters are added 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 Schematic diagram of the three-dimensional density model obtained by inverting a satellite gravity anomaly based on a three-dimensional density reference model, a three-dimensional density upper limit model, and a three-dimensional density lower limit model by the method for three-dimensional crustal density modeling based on seismic constraints provided by an embodiment of the present invention. Detailed implementation manner

[0068] Hereinafter, exemplary embodiments of the present disclosure will be described in more detail with reference to the accompanying drawings. Although the exemplary embodiments of the present disclosure are shown in the 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 so that the present disclosure can be more thoroughly understood and the scope of the present disclosure can be completely conveyed to those skilled in the art.

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

[0070] Step 11: Obtain two-dimensional reflection seismic profiles and two-dimensional refraction seismic velocity models in the study area, and construct constraint information of the main density horizons, where the constraint information of the main density horizons includes geometric constraint information and density constraint information;

[0071] Step 12: Based on the satellite gravity anomaly, perform anomaly extraction and depth inversion on the main density interfaces according to the geometric constraint information and density constraint information to obtain the depths of the main density interfaces;

[0072] Step 13: According to the depths of the main density interfaces 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;

[0073] Step 14: Based on the satellite gravity anomaly, perform three-dimensional constrained inversion using a regularization inversion method according to the three-dimensional density reference model and upper and lower limit constraint models to construct a three-dimensional crustal density model.

[0074] In the embodiments of the present invention, by identifying, digitizing, and discretizing the main density interfaces in the obtained two-dimensional reflection seismic profile results and two-dimensional refraction seismic velocity model, geometric constraint information can be established, and the scatter 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 performing conversions, density constraint information of the main density layers and lateral variable density constraint information of the main interfaces 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, a Moho surface depth that more conforms to the actual constraint conditions can be obtained, providing a basis for subsequent modeling; by introducing the basement geometric constraint information and the basement lateral variable density constraint information into the basement inversion, a basement depth that more conforms to the actual constraint conditions can be obtained, providing a basis for subsequent modeling; by establishing a three-dimensional reference model and density upper and lower limit constraint models that conform to the prior constraint information according to the obtained bathymetric data, the inverted basement and Moho surface depths, and combining the density constraint information, multiple constraint information can be integrated into three-dimensional constraint information; by introducing the three-dimensional reference model and density upper and lower limit constraint models into the three-dimensional density inversion and using the regularization inversion method to fit the satellite gravity anomaly, a three-dimensional density model that more conforms to the prior constraint information can be constructed, realizing an effective prediction of the underground density distribution and providing an important reference for geological exploration and resource evaluation.

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

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

[0077] Step 1112, interpret the horizon structures of the two-dimensional reflection seismic profile and two-dimensional refraction seismic velocity model, identify 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 depths of the basement interface and the Moho surface, and establish the basement interface geometric constraint information and the Moho surface geometric constraint information.

[0079] In the embodiments of the present invention, by interpreting the two-dimensional reflection seismic profile and the two-dimensional refraction seismic velocity model, different stratigraphic interfaces can be identified, and geometric information of the underground stratigraphic interfaces can be provided, including the depth, morphology, etc. of the interfaces; according to the interpretation results, digitization and discretization are performed to obtain the scatter depth information of the main interfaces, and geometric constraint information can be established; the geometric constraint information can be used for subsequent interface inversion, for constraining the extraction of the basement and Moho gravity anomalies, and controlling the inversion processes of the basement and Moho, which can improve the accuracy of the inversion results of the basement and Moho; establishing the geometric constraint information is one of the key steps in seismic-constrained three-dimensional density modeling; transforming the two-dimensional seismic structural prior information into geometric constraint information and introducing it into the three-dimensional density modeling of the crust makes the modeling results more in line with the actual geological situation; the establishment of the geometric constraint information provides a basis for subsequent construction of the three-dimensional density reference model, density upper and lower limit models, and three-dimensional constrained gravity inversion.

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

[0081] Step 1121, according to the two-dimensional refraction seismic velocity model, extracting the seismic longitudinal wave velocity information of the sedimentary layer, crust and mantle in different tectonic regions;

[0082] Step 1122, according to the conversion relationship between velocity and density, converting the seismic longitudinal wave velocity information into density parameters, and establishing the density constraint information of the sedimentary layer, crust and mantle in different regions, where the density constraint information includes the maximum value, minimum value and average value;

[0083] Step 1123, according to the density constraint information of the sedimentary layer, crust and mantle, establishing the lateral variable density constraint information of the basement interface and the lateral variable density constraint information of the Moho interface.

[0084] In the embodiments of the present invention, through digitization and discretization of the two-dimensional refraction seismic velocity model, seismic longitudinal wave velocity information of sedimentary layers, the crust, and the mantle in different tectonic regions can be extracted; according to the velocity-density conversion relationship, the seismic longitudinal wave velocity can be converted into density parameters, and the maximum, minimum, and average values of the densities of sedimentary layers, the crust, and the mantle in different regions can be obtained to establish density constraint information for the main density layers; according to the density constraint information of sedimentary layers, the crust, and the mantle in different regions, lateral variable density constraint information for the basement and the Mohorovicic discontinuity can be established; the lateral variable density constraint information for the basement and the Mohorovicic discontinuity can be used for depth inversion of the basement and the Mohorovicic discontinuity to improve the reliability of the inversion results; converting two-dimensional seismic velocity prior information into density constraint information and introducing it into three-dimensional crustal density modeling makes the modeling results more consistent with the actual geological situation; establishing density constraint information is one of the key steps in 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] As Figures 1 to 5 shown, step 12, based on satellite gravity anomalies, perform anomaly extraction and depth inversion on the main density interfaces according to the geometric constraint information and density constraint information to obtain the depths of the main density interfaces, including:

[0086] Step 1211, obtain satellite gravity anomaly data, bathymetric data, and sediment thickness data;

[0087] Step 1212, according to the bathymetric data, divide the seawater layer into a first series of upright hexahedrons, and perform forward modeling using the spatial domain method to obtain the seawater correction value;

[0088] Step 1213, according to the sediment thickness data, divide the sedimentary layer into a second series of upright hexahedrons, and perform forward modeling based on the compaction model using the spatial domain method to obtain the sedimentary layer correction value. The densities at different depths are:

[0089]

[0090] where g 0 is the Mohorovicic discontinuity gravity anomaly, ρ z is the density of the sedimentary layer at a depth of z, ρ f is the density of the fluid in the sediment, ρ g is the density of the rock particles, is the initial porosity, and c is the porosity depth attenuation coefficient;

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

[0092] Step 1215: Based on the Mohorovicic discontinuity geometric constraint information, perform potential field separation on the crust Bouguer gravity anomaly data to obtain the Mohorovicic discontinuity gravity anomaly data;

[0093] Step 1216: Based on the Mohorovicic discontinuity geometric constraint information, use the linear regression method to obtain an initial estimate of the Mohorovicic discontinuity depth

[0094] Step 1217: According to the Mohorovicic discontinuity gravity anomaly data, introduce the Mohorovicic discontinuity geometric constraint information in a soft constraint manner. Based on the Mohorovicic discontinuity lateral variable density constraint information, iteratively calculate the depth of the Mohorovicic discontinuity. The formula is:

[0095]

[0096] where, is the Mohorovicic discontinuity gravity anomaly, is the forward gravity response of the (k - 1)-th Mohorovicic discontinuity depth, and are the results of the k-th and (k - 1)-th iterations respectively. G is the gravitational constant, ρ moho is the Mohorovicic discontinuity lateral variable density information, and Δp moho is the depth correction amount obtained by linear regression based on the Mohorovicic discontinuity geometric constraint information in the k-th iteration.

[0097] In the embodiments of the present invention, forward modeling is performed based on bathymetric data to obtain a sea water correction value. First, the 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 through numerical calculations. According to the forward modeling calculation results, the gravitational effect of the ocean bottom, that is, the sea water correction value, is obtained. The sea water correction value represents the contribution of the ocean to the gravity anomaly. Based on 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 gravitational effect of the sediment layer is simulated. The compaction model is a model that calculates the gravitational effect according to the physical properties of the sediment layer, such as density and thickness. According to the forward modeling calculation results, the gravitational effect of the sediment layer, that is, the sediment layer correction value, is obtained. The sediment layer correction value represents the contribution of the sediment layer to the gravity anomaly. The satellite gravity anomaly data is corrected according to the sea water correction value and the sediment layer correction value to obtain the crustal Bouguer gravity anomaly data. First, the satellite gravity anomaly data is used, and the sea water correction value and the sediment layer correction value are subtracted to obtain the gravity anomaly data caused by other factors in the crust (such as crustal density changes, tectonic deformations, etc.). The obtained gravity anomaly data is the crustal Bouguer gravity anomaly data, which reflects the influence of other factors in the crust except for the sea water and the sediment layer. Based on the Mohorovicic discontinuity geometric constraint information, potential field separation is performed on the crustal Bouguer gravity anomaly data to obtain the Mohorovicic discontinuity gravity anomaly data. First, the Mohorovicic discontinuity geometric constraint information is used to determine the depth and shape of the Mohorovicic discontinuity. Then, according to the potential field separation method, the crustal Bouguer gravity anomaly data is decomposed into the Mohorovicic discontinuity gravity anomaly data and the residual crustal gravity anomaly data. The Mohorovicic discontinuity gravity anomaly data reflects the gravity anomaly caused by the change of the Mohorovicic discontinuity in the crust and can be used to study the structure and tectonics of the crust. According to the Mohorovicic discontinuity geometric constraint information, the depth of the Mohorovicic discontinuity can be initially estimated more reliably by using the method of linear regression to obtain the average depth of the Mohorovicic discontinuity for the Mohorovicic discontinuity iterative inversion process. The Mohorovicic discontinuity geometric constraint information is introduced into the interface inversion process in a soft constraint manner, which helps to eliminate the blurred imaging of the Mohorovicic discontinuity by two-dimensional reflection seismic. Based on the lateral variable density constraint, the interface distribution is obtained through iterative optimization, making the result more in line with the actual geological situation. The depth of the Mohorovicic discontinuity 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 Mohorovicic discontinuity depth inversion under geometric constraint and lateral variable density constraint is one of the key steps in seismic-constrained three-dimensional density modeling of the crust. By obtaining the depth of the Mohorovicic discontinuity, it can provide a basis for subsequent three-dimensional constrained gravity inversion and the construction of a three-dimensional density model.

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

[0099] In step 1221, forward modeling is performed according to the depth of the Moho surface to obtain the Moho correction value;

[0100] In step 1222, correlation analysis is performed based on the satellite gravity anomaly data and bathymetric data to obtain the optimal density parameter for seawater correction. Using the bathymetric data, the seawater layer is divided into a series of upright hexahedrons, and the seawater correction value is obtained by forward calculation using the spatial domain method;

[0101] In step 1223, the satellite gravity anomaly data is corrected according to the Moho correction value and the seawater correction value to obtain the Bouguer gravity anomaly data of the sedimentary layer;

[0102] In step 1224, under the geometric constraint of the basement interface, the wavelet transform method is used to separate the potential field of the Bouguer gravity anomaly data of the sedimentary layer to obtain the gravity anomaly data of the basement interface;

[0103] In step 1225, according to the geometric constraint information of the basement interface, the two-dimensional interpolation method is used to obtain the initial estimated basement depth p bass 0 ;

[0104] In step 1226, according to the gravity anomaly data of the basement interface, the depth constraint information of the basement interface 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 obtained by iterative calculation. The formula is:

[0105]

[0106] Among them, and are the results of the k-th and (k - 1)-th basement depth iterations respectively. J is the Jacobian matrix, J = 2πGρ bass I (where I is the identity matrix), ρ bass is the lateral variable density information of the basement, G is the gravitational constant, R is the first 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 embodiments of the present invention, forward modeling is performed according to the depth of the Moho discontinuity to obtain the Moho correction value. First, the depth data of the Moho discontinuity is acquired to determine the shape and position of the Moho interface. Then, based on the forward modeling method, the influence of the Moho interface on gravity anomalies is simulated. Forward modeling refers to simulating the propagation and response process of the geophysical field through numerical calculations. According to the forward calculation results, the gravitational effect of the Moho discontinuity, that is, the Moho correction value, is obtained. The Moho correction value represents the contribution of the Moho interface to gravity anomalies; using bathymetric data, the seawater layer is divided into a series of upright hexahedrons, and the position and size of each hexahedron are determined; the forward calculation is performed using the spatial domain method, and according to the seawater density parameters and the distribution of the seawater layer, the seawater correction value is calculated; the calculated seawater correction value is analyzed and verified to ensure the reliability and accuracy of the results; the finally obtained seawater correction value can be used to correct satellite gravity anomaly data, improving the accuracy and reliability of the data; first, the satellite gravity anomaly data and the bathymetric data are subjected to correlation analysis to determine the optimal density parameter for seawater correction. Correlation analysis is used to judge the degree of correlation between two data sets. Through the forward modeling method, using the bathymetric data and the optimal density parameter, the contribution of seawater to gravity anomalies is simulated. The forward calculation is carried out by an analytical method. According to the forward calculation results, the gravitational effect of seawater, that is, the seawater correction value, is obtained. The seawater correction value represents the correction of seawater to gravity anomalies; the satellite gravity anomaly data is corrected according to the Moho correction value and the seawater correction value 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, and 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 obtained gravity anomaly data reflects the contribution of other factors to gravity anomalies except for the Moho interface and seawater, mainly from the sedimentary layer in the crust; the wavelet transform method is used to process the Bouguer gravity anomaly data of the sedimentary layer, and it is decomposed into wavelet coefficients of different scales and frequencies for subsequent analysis; according to the geometric constraints of the basement interface, appropriate wavelet coefficients are selected for inverse transformation to obtain the gravity anomaly data of the basement interface; the obtained gravity anomaly data of the basement interface is analyzed and verified to ensure that the results meet the expectations, and necessary corrections and adjustments are made; the finally obtained gravity anomaly data of the basement interface can be used for further geological interpretation and research, providing important information for geological structures and sedimentary processes; the gravity anomaly data of the basement interface reflects the gravity anomaly caused by the changes in the crustal basement and is of great significance for studying the structure and geological characteristics of the crust; introducing the basement geometric constraints into the basement interface inversion process in a hard constraint manner can improve the reliability of the basement inversion results; based on the lateral variable density constraint of the basement, the distribution of the basement interface is obtained through an iterative optimization method, making the results more in line with the actual geological situation; the depth of the basement 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 inversion of the basement interface under geometric constraints and lateral variable density constraints is one of the key steps in 3D crustal density modeling based on seismic constraints; by obtaining the depth of the basement, it can provide a basis for subsequent 3D constrained gravity inversion and the construction of a 3D density model.

[0108] As Figures 1 to 5 shown, in step 13, according to the depth of the main density interface and the density constraint information, a 3D density reference model and upper and lower limit constraint models that conform to the prior constraint information are established, including:

[0109] In step 131, according to the bathymetric data of the seabed, 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 3D geometric constraint model is established. The four layers include the seawater layer, the sedimentary layer, the crust, and the mantle.

[0110] In step 132, the average density corresponding to each of the four layers is added to the 3D geometric constraint model to construct a 3D density reference model.

[0111] In step 133, the maximum density corresponding to each of the four layers is added to the 3D geometric constraint model to construct a 3D density upper limit model.

[0112] In step 134, the minimum density corresponding to each of the four layers is added to the 3D geometric constraint model to construct a 3D density lower limit model.

[0113] In the embodiment of the present invention, according to the obtained bathymetric data of the seabed, the inverted basement and Moho depths, the underground three-dimensional medium can be divided into four layers: the seawater layer, the sedimentary layer, the crust, and the mantle, and a 3D geometric constraint model is established; in the 3D geometric constraint model, the average density, maximum density, and minimum density of the four layers are added to respectively construct a 3D density reference model, a density upper limit model, and a density lower limit model; the 3D density reference model is the basis for 3D constrained inversion, providing the general distribution and characteristics of the underground structure and providing a reference for subsequent density inversion; the physical property upper and lower limit constraint models are used to constrain the range of physical property parameters of the density model, making the density model more in line with the actual geological situation; the 3D reference model and density upper and lower limit constraint models that conform to the prior constraint information can improve the reliability of modeling; establishing a 3D reference model and physical property upper and lower limit constraint models that conform to the prior constraint information is one of the key steps in 3D density modeling based on seismic constraints; the reference model and physical property upper and lower limit constraint models provide a basis for subsequent regularization inversion.

[0114] As Figures 1 to 5 shown, in step 14, based on the satellite gravity anomaly, according to the 3D density reference model and upper and lower limit constraint models, a 3D constrained inversion is carried out using a regularization inversion method to construct a 3D crustal density model, including:

[0115] Step 141: According to the three-dimensional density reference model and the three-dimensional density upper and lower limit constraint model, use the iterative inversion method to fit the satellite gravity anomaly data, and obtain a series of three-dimensional density inversion results by setting different fitting errors and smoothness. The formula is as follows:

[0116]

[0117] where 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, and W d is the data weight matrix; W m is the model weight matrix; μ is the regularization parameter; the range of m is constrained by the logarithmic barrier method:

[0118]

[0119] where m max and m min are the three-dimensional density upper limit model and the lower limit model respectively, and λ is the weight factor;

[0120] Step 142: By comparing the similarity between the three-dimensional density inversion result and the two-dimensional reflection seismic profile, obtain the similarity comparison result, and take the result with high similarity as the final three-dimensional density inversion result;

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

[0122] In the embodiment of the present invention, by introducing the three-dimensional reference model and the 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; by comparing the two-dimensional reflection seismic profile with the three-dimensional slice, observing the similarity on the underground structure, comparing the characteristics such as the formation interface and lithology change, and selecting the one with the highest similarity as the three-dimensional density inversion result according to the comparison result, the inversion result can meet the prior constraint information; using three-dimensional visualization software, the three-dimensional density inversion result can be visualized to establish a three-dimensional density model of the crust; through visualization, the underground density distribution can be intuitively observed, which helps geologists and geophysicists better understand the underground density distribution.

[0123] Optionally, in an embodiment of the present invention, extract density slices at different depths and thicknesses of different horizons based on the three-dimensional density model of the crust, and determine favorable oil and gas areas based on the regional oil and gas accumulation law, including:

[0124] Determine the main hydrocarbon-generating horizons and traps according to the regional oil and gas accumulation law;

[0125] Extract the depths of the top and bottom interfaces of the hydrocarbon-generating horizons and the depth of the trap according to the 2D reflection seismic profile;

[0126] Obtain the average densities of the top and bottom interfaces of the hydrocarbon-generating horizons at the corresponding positions of the 2D reflection seismic profile according to the 3D density model;

[0127] Extract the 3D depths of the top and bottom interfaces of the hydrocarbon-generating horizons from the 3D density model according to the average densities of the top and bottom interfaces of the hydrocarbon-generating horizons, and calculate the thickness of the hydrocarbon-generating horizons by subtracting the interface depths;

[0128] Extract density slices at different depths according to the depth of the trap according to the 3D density model;

[0129] Determine the depth and position of the uplift of the hydrocarbon-generating horizons according to the high-density features in the density slices;

[0130] Determine the favorable areas for oil and gas according to the thickness of the hydrocarbon-generating horizons and the depth and position of the uplift of the hydrocarbon-generating horizons.

[0131] In the embodiments of the present invention, by studying the regional hydrocarbon accumulation law, the main hydrocarbon generation horizons and traps can be determined; according to the two-dimensional reflection seismic profile, the depths of the top and bottom interfaces of the hydrocarbon generation horizons and the depth of the trap are extracted. First, the reflection waveforms on the reflection seismic profile are analyzed to identify the top and bottom interfaces of the hydrocarbon generation horizons, and then a seismic interpretation software or other relevant tools are used to mark the depths of the top and bottom interfaces of the hydrocarbon generation horizons on the profile; then, according to the depths of the top and bottom interfaces of the hydrocarbon generation horizons, the depth of the trap is determined. The trap is composed of the area between the top and bottom interfaces of the hydrocarbon generation horizons; according to the three-dimensional density model, the average densities of the top and bottom interfaces of the hydrocarbon generation horizons at the corresponding positions of the two-dimensional reflection seismic profile are obtained. First, the depths of the top and bottom interfaces of the hydrocarbon generation horizons in the two-dimensional reflection seismic profile are corresponded to the three-dimensional density model. In the three-dimensional density model, the average density values at the positions corresponding to the depths of the top and bottom interfaces of the hydrocarbon generation horizons are extracted, and the average density values of the top and bottom interfaces of the hydrocarbon generation horizons are recorded; according to the average densities of the top and bottom interfaces of the hydrocarbon generation horizons, the three-dimensional depths of the top and bottom interfaces of the hydrocarbon generation horizons are extracted from the three-dimensional density model, and the thickness of the hydrocarbon generation horizon is calculated by subtracting the interface depths. First, the average density values of the top and bottom interfaces of the hydrocarbon generation horizons are used to determine the three-dimensional depths of the top and bottom interfaces of the hydrocarbon generation horizons in the three-dimensional density model, and the thickness of the hydrocarbon generation horizon is calculated by subtracting the three-dimensional depths of the top and bottom interfaces; according to the three-dimensional density model, density slices at different depths are extracted according to the depth of the trap. First, according to the depth range of the trap, the density slices at the corresponding depths are extracted from the three-dimensional density model. The density slice is an image of the density distribution within a specific depth range; according to the high-density features in the density slice, the depth and position of the uplift of the hydrocarbon generation horizon are determined. In the density slice, the high-density features are observed, and by determining the depth and position of the high-density features, the depth and position of the uplift of the hydrocarbon generation horizon are determined. Combining the thickness of the hydrocarbon generation horizon and the depth and position of the uplift, the range of the favorable hydrocarbon area is determined; according to the geological conditions and sedimentary environment, considering the influence of the thickness and uplift of the hydrocarbon generation horizon on hydrocarbon accumulation, the favorable hydrocarbon area is further determined.

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

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

[0134] A processing module 22, configured to perform anomaly extraction and depth inversion on a main density interface based on satellite gravity anomalies according to the geometric constraint information and density constraint information, so as 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 according to the depth of the main density interface and the density constraint information; and perform three-dimensional constrained inversion by using a regularization inversion method based on satellite gravity anomalies according to the three-dimensional density reference model and the upper and lower limit constraint models to construct a three-dimensional crustal density model.

[0135] An embodiment of the present invention further provides a computing device, including: a processor and a memory storing a computer program. When the computer program is run by the processor, the above-mentioned method is executed. All implementation manners in the above method embodiments are applicable to this embodiment and can achieve the same technical effects.

[0136] An embodiment of the present invention further provides a computer-readable storage medium storing instructions. When the instructions are run on a computer, the computer is caused to execute the above-mentioned method. All implementation manners in the above method embodiments are applicable to this embodiment and can achieve the same technical effects.

[0137] Those of ordinary skill in the art can realize that the units and algorithm steps of each example described in combination with the embodiments disclosed herein can be implemented by electronic hardware, or by a combination of computer software and electronic hardware. Whether these functions are executed in a hardware or software manner depends on the specific application and design constraints of the technical solution. Professional technicians can use different methods to implement the described functions for each specific application, but such implementation should not be considered to exceed 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 foregoing method embodiments and will not be elaborated herein.

[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 merely illustrative. For example, the division of the units is only a logical function division, and there may be other division methods in actual implementation. For example, 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 displayed or discussed couplings or direct couplings or communication connections to each other can be through some interfaces, and the indirect couplings or communication connections of the systems or units can be in electrical, mechanical, or other forms.

[0140] The unit described as a separation component may or may not be physically separated. The component presented as a unit may or may not be a physical unit, that is, it may be located in one place or may be distributed across multiple network units. Some or all of the units can be selected according to actual needs to achieve the purpose of the solution of this embodiment.

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

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

[0143] In addition, it should be noted that in the systems and methods of the present invention, obviously, each component or each step can be decomposed and / or recombined. These decompositions and / or recombinations should be regarded as equivalent solutions of the present invention. And, the steps of performing the above series of processes can naturally be executed in chronological order according to the described order, but it is not necessary to execute them in chronological order. Some steps can be executed in parallel or independently of each other. For those of ordinary skill in the art, it is understandable that all or any steps or components of the methods and systems of the present invention can be implemented in any computing system (including processors, storage media, etc.) or in a network of computing systems in the form of 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] Accordingly, the object of the present invention can also be achieved by running a program or a set of programs on any computing system. The computing system can be a well-known general-purpose system. Therefore, the object of the present invention can also be achieved only by providing a program product containing program code for implementing the method or system. That is to say, such a program product also constitutes the present invention, and a 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 noted that in the system and method of the present invention, obviously, each component or each step can be decomposed and / or recombined. These decompositions and / or recombinations should be regarded as equivalent solutions of the present invention. And, the steps of performing the above series of processes can naturally be executed in chronological order according to the described order, but it is not necessary to be executed in chronological order. Some steps can be executed in parallel or independently of each other.

[0145] The above is the preferred embodiment of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the principle of the present invention, several improvements and refinements can be made, and these improvements and refinements should also be regarded as the protection scope 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 upper and lower limit constraint models, the regularized inversion method is used to perform three-dimensional constrained inversion to construct a three-dimensional density model of the crust. 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, ρ z is the sediment density at depth z, ρ f is the density of the fluid in the sediment, ρ g is the density of 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 constraints of the Moho surface, the initial Moho depth estimate 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: in, The Moho gravity anomaly. is the forward gravity response of the k-1th Moho depth, and are the results of the kth and k-1th iterations, respectively, G is the gravity constant, ρ moho is the lateral density variation information of the Moho surface, Δp moho is the depth correction obtained by linear regression based on the Moho geometric constraint information for the kth time.

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 three-dimensional crustal density modeling based on earthquake constraints according to claim 1, 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 substrate interface, the initial substrate depth estimation p is obtained by using the two-dimensional interpolation method. bass 0 ; 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: in, and are the results of the kth and k-1th base depth iterations, respectively, J is the Jacobian matrix, J = 2πGρ bass I (where I is the identity matrix), ρ bass is the lateral variable density information of the basement, 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.

5. The method for three-dimensional crustal density modeling based on earthquake constraints according to claim 4, 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.

6. The method for three-dimensional crust density modeling based on earthquake constraints according to claim 5, 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, W d is the data weight matrix; W m is the model weight matrix; μ is the regularization parameter; the logarithmic barrier method is used to constrain the range of m: Among them, m max and m min 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.

7. A system for the method of three-dimensional crust density modeling based on seismic constraints as claimed in claim 1, 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 models, use a regularized inversion method to perform three-dimensional constraint inversion to construct a three-dimensional density model of the crust.

8. 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 6.

9. 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 6.

Images