A binary structure-constrained gravity inversion method based on seismic data

Abstract

A binary structure-constrained gravity inversion method based on seismic data, including using a high-resolution velocity model obtained by inverting seismic first-arrival travel-time data as a guiding model. Combining the fuzzy c-means clustering technique and the threshold technique is used to divide the background area and the target area of the guiding model, and the structural information of the velocity model can be accurately extracted. In the binary structure-constrained inversion algorithm, since the binary structure-constrained model only contains two elements, 0 and 1, the inversion result shows a certain sparsity sharpness and can better depict the boundary. At the same time, the inversion only calculates in the target area while the background area remains unchanged, so the inversion effectively reduces the dimension of the solution space and reduces the non-uniqueness of the inversion.

Description

Technical Field

[0001] The present invention relates to the technical field of two-dimensional gravity data inversion processing, and particularly provides a high-precision method for gravity inversion with binary structure constraint based on seismic data. Background Art

[0002] Gravity exploration has the advantages of being portable, fast, and low-cost, and has been widely used in the fields of studying the internal structure of the crust, exploring the distribution of solid minerals and oil and gas resources. Gravity inversion is one of the important means for quantitative interpretation of deep resource exploration, but there is a serious problem of non-uniqueness in inversion. Especially in terms of vertical resolution ability, it is difficult to obtain a high-resolution density reconstruction result. Fortunately, electromagnetic methods with rich frequency information, high-resolution seismic methods, and borehole prior information can provide high-resolution vertical structure information. For conventional image-guided methods, prior information on rock physical properties is required, which has certain limitations. However, for cross-gradient constraint joint inversion, two geophysical data need to be inverted simultaneously, and the calculation process is relatively complex. Moreover, when the convergence rates of different data differ too much or the resolution differences of the two data are too large, it is difficult for the cross-gradient constraint joint inversion to reconstruct a high-resolution density model. Therefore, how to incorporate this high-resolution structure information into gravity inversion has always been a challenging problem. Summary of the Invention

[0003] The purpose of the present invention is to provide a gravity inversion method with binary structure constraint based on seismic data to solve the problem of low resolution in existing gravity inversion.

[0004] As Figure 1 shown, the present invention includes the following steps:

[0005] 1) Input the observed data of seismic first arrival travel time, input the initial velocity model and reference model, and the initial regularization factor;

[0006] 2) Construct the two-dimensional seismic data inversion objective function Φ1(m1);

[0007] 3) Use the Gauss-Newton method to optimize and solve the two-dimensional seismic data inversion objective function, and use the obtained velocity model as the guiding model;

[0008] 4) Combine the fuzzy c-means clustering algorithm (FCM) and threshold technology to separate the target area and background area of the guiding model to obtain a binary structure constraint model;

[0009] 5) Input the gravity observed data, input the initial density model and reference model, and the initial regularization factor;

[0010] 6) Construct the two-dimensional gravity data inversion objective function Φ2(m2) based on binary structure constraint;

[0011] 7) Under the constraint of the binary structure constraint model, the Gauss-Newton method is used to optimize and solve the objective function of two-dimensional gravity data inversion to obtain the density model.

[0012] Advantages of the present invention:

[0013] 1. By combining the fuzzy c-means clustering technique and the threshold technique, the present invention can accurately divide the background area and the target area of the velocity model, and then effectively extract the structural information of the velocity model.

[0014] 2. The binary structure constraint model of the present invention is a matrix composed of two elements, 0 and 1. Therefore, the density model under the binary structure constraint will exhibit sparsity, which can better depict the interface of the anomaly body and avoid the ambiguity brought by the smooth constraint.

[0015] 3. The binary structure constraint gravity inversion algorithm of the present invention only performs calculations in the target area, and the background area remains unchanged, effectively reducing the dimension of the inversion solution space, and then reducing the non-uniqueness of the inversion. Compared with the separate gravity inversion and the joint inversion with cross-gradient constraint, this method obtains a more accurate spatial geometric shape and physical property parameter values of the anomaly body. Description of the drawings

[0016] Figure 1 is a flow chart of the binary structure constraint gravity inversion method based on seismic data;

[0017] Figure 2 is a schematic diagram of the density and velocity theoretical models;

[0018] Figure 3 is a schematic diagram of the binary structure constraint model;

[0019] Figure 4 is a schematic diagram of the density models obtained by different inversion methods. Detailed implementation manners

[0020] As Figure 1 shown, a binary structure constraint gravity inversion method based on seismic data includes the following steps:

[0021] 1) Input the observed data of the seismic first arrival travel time, the initial velocity model and the reference model, and the initial regularization factor;

[0022] 2) Construct the objective function Φ1(m1) of two-dimensional seismic data inversion;

[0023] 3) Use the Gauss-Newton method to optimize and solve the objective function of two-dimensional seismic data inversion, and use the obtained velocity model as the guiding model;

[0024] 4) Combine the fuzzy c-means clustering algorithm (FCM) and threshold technology to separate the target area and background area of the guidance model, and obtain a binary structure constraint model;

[0025] 5) Input gravity observation data, input the initial density model and reference model, and the initial regularization factor;

[0026] 6) Construct a two-dimensional gravity data inversion objective function Φ2(m2) based on binary structure constraints;

[0027] 7) Under the constraint of the binary structure constraint model, use the Gauss-Newton method to optimize and solve the two-dimensional gravity data inversion objective function to obtain a density model.

[0028] In step 1), input the seismic first arrival travel time observation data, input the initial velocity model and reference model, and the initial regularization factor. The seismic first arrival travel time observation data is obtained by calculating the forward response through establishing a velocity theoretical model.

[0029] In step 2), construct a two-dimensional seismic data inversion objective function Φ1(m1), and the expression is as follows:

[0030] Φ1(m1) = ||W d1 (d1 - f(m1))|| 2 + α1·||W m1 (m1 - m 1ref )|| 2 (1)

[0031] where f(m1) is the seismic forward response; d1 is the seismic observation data; W d1 is the diagonal inverse matrix of seismic data noise error; W m1 is the velocity model smoothing matrix; m 1ref is the velocity reference model; α1 is the regularization factor of the seismic method.

[0032] In step 3), use the Gauss-Newton method to optimize and solve the two-dimensional seismic data inversion objective function, and the obtained velocity model is used as the guidance model.

[0033] The expression of the velocity model obtained by using the Gauss-Newton method is as follows:

[0034]

[0035] where A is the Jacobian matrix of the seismic method, m 1k is the velocity model parameter, k is the number of inversion iterations, and T is the matrix transpose symbol.

[0036] Due to the inherent physical properties of gravity inversion, it exhibits low vertical resolution. To reduce the non-uniqueness of gravity inversion, it is necessary to impose prior structural information constraints from other geophysical inversions, drilling, or geological reference models. Compared with electromagnetic and magnetic data, seismic data has higher vertical resolution. Therefore, using the velocity model obtained from seismic data inversion as a guiding model has higher credibility.

[0037] In step 4), the fuzzy c-means clustering algorithm (FCM) and threshold technique are combined to separate the target area and background area of the guiding model, obtaining a binary structure constraint model.

[0038] The expression of the FCM objective function is as follows:

[0039]

[0040]

[0041] where m i represents the i-th velocity model, M represents the number of velocity models, C represents the number of clusters, v j is the clustering center of the j-th cluster, q is the fuzzy coefficient, usually set to 2, and u ij represents the membership degree of the i-th velocity model m i belonging to the j-th cluster, and the membership degree is between [0, 1].

[0042] By minimizing the objective function (Equation 1), the membership degree u ij and the clustering center v j are obtained, and the expressions are as follows:

[0043]

[0044]

[0045] Calculate the absolute value (|v j -m b |) of the difference between the clustering center v j and the background value m b of the velocity model, and take the minimum value as the threshold ξ. Then, when |m i -m b |<ξ (i = 1, 2,..., M), set Q i = 0, and for others set Q i = 1. Here, Q is called the binary structure constraint model, which is a matrix composed of only 0 and 1.

[0046] In step 5), input the gravity observation data, the initial density model and the reference model, and the initial regularization factor. The gravity observation data is obtained by calculating the forward response through establishing a velocity theoretical model.

[0047] In step 6), construct the two-dimensional gravity data inversion objective function Φ2(m2) based on binary structure constraints, and the expression is as follows:

[0048] Φ2(m2) = ||W d2 (d2 - f(m2)) 2 +α2·||W m2 (m2 - m 2ref )|| 2 (7)

[0049] Constraint condition: m2 = diag(R)·m2 (8)

[0050] Wherein, f(m2) is the forward response; d2 is the observed data; W d2 is the diagonal inverse matrix of data noise error; W m2 is the model smoothing matrix; m 2ref is the reference model; α2 is the regularization factor.

[0051] In step 7), under the constraint of the binary structure constraint model, use the Gauss-Newton method to optimize and solve the two-dimensional gravity data inversion objective function, and the expression of the (k + 1)-th density model obtained is as follows:

[0052]

[0053] m 2k+1 = diag(R)·m 2k+1 (10)

[0054] Wherein, A2 is the Jacobian matrix of the gravity method, and m 2k is the density model parameter of the k-th iteration. The gravity inversion is only calculated in the target area, and the background area remains unchanged, effectively reducing the dimension of the inversion solution space. Compared with the single gravity inversion, the binary structure constraint gravity inversion can obtain a higher-resolution density model.

[0055] Next, verify the binary structure constraint gravity inversion method based on seismic data provided by the present invention.

[0056] As Figure 2 shown, a rectangular combined model is designed. The model consists of two rectangular bodies of different sizes, and the density model is as Figure 2 shown in a, and the velocity model is as Figure 2 shown in b. The size of the background area is 8 km × 3 km, and the background velocity and density are 4000 m / s and 0 g / cm3 The target area contains two rectangular anomalies. The upper rectangle has a size of 2 km × 0.5 km, with velocities and densities of 3000 m / s and 0.7 g / cm respectively. 3 The lower rectangle has a size of 4 km × 0.7 km, with velocities and densities of 5000 m / s and 1 g / cm respectively. 3 There are 29 gravity observation points evenly distributed on the survey line from 0 - 8 km. The model is divided into 80×30 horizontal and vertical cells, each with a size of 100×100 m. This test directly uses the true velocity model as the guiding model, and then reconstructs the density model through separate inversion, cross-gradient joint inversion algorithm, and binary structure-constrained gravity inversion. In all inversion algorithm calculations, the initial density model is the background density model.

[0057] First, a separate inversion is performed on the gravity data set. The separate inversion stops after 5 iterations, and the reconstructed density model is as shown in Figure 3 (a). As expected, the reconstructed density model cannot clearly define the location and geometry of the target. The boundaries of the anomalies are not clear, and the upper and lower anomalies cannot be distinguished, which is usually due to the limited vertical resolution of the gravity data. In the cross-gradient joint inversion algorithm iteration, by fixing the velocity model at each iteration and using the cross-gradient function to structurally constrain the density model, this algorithm stops after 5 iterations, and the reconstructed density model is as shown in Figure 3 (b). It can be found that at the boundaries of the anomalies in the velocity model, the density model obtains a structural constraint effect, and the sharp boundaries of the two rectangular anomalies can be restored. However, the structural constraint effect is poor outside the target area. This is mainly because there is no gradient change in the velocity values of the velocity model outside the target area, and the cross-gradient value in this area is always zero, resulting in the cross-gradient function being unable to play a structural constraint role. Therefore, false anomalies appear outside the rectangular anomalies, showing a structural difference from the true density model. This phenomenon also indicates that the cross-gradient joint inversion algorithm needs to be used with caution when there are large differences in the convergence rates of different geophysical methods, and it is necessary to adjust the two methods to have the same convergence rate as much as possible to improve the joint effect. In the FCM space-constrained joint inversion algorithm, the target area of the true velocity model is extracted through FCM clustering technology and threshold technology. We set the number of clusters C = 2 and continuously update the cluster centers. Finally, the target area and the background area are divided into two categories. The red area is the target area, and the blue area is the background area. Figure 4 This is the obtained binary structure-constrained model. Under the binary structure constraint, the gravity inversion stops after 8 iterations, and the reconstructed density model is as shown in Figure 3(c). It can be found that due to the reduction of the inversion solution space, the multi-solution problem of gravity inversion is reduced, resulting in a reconstructed density model closer to the true model. Compared with the cross-gradient constraint joint inversion algorithm, a more accurate spatial geometric shape and physical property parameter values of the anomaly body are obtained, improving the resolution ability of gravity inversion.

Claims

1. A binary structure-constrained gravity inversion method based on seismic data, characterized in that: Including the following steps: 1) Input the observed data of the first arrival time of seismic waves, the initial velocity model, the reference model, and the initial regularization factor; 2) Construct the objective function for 2D seismic data inversion ; 3) Use the Gauss-Newton method to optimize and solve the inversion objective function of two-dimensional seismic data, and the obtained velocity model is used as the guiding model; 4) Combine fuzzy c The mean clustering algorithm and threshold technology are used to separate the target area and background area of the guidance model, and a binary structure constraint model is obtained; 5) Input the gravity observation data, the initial density model, the reference model, and the initial regularization factor; 6) Construct a two-dimensional gravity data inversion objective function based on binary structure constraints ; 7) Under the constraint of the binary structure constraint model, use the Gauss-Newton method to optimize and solve the inversion objective function of two-dimensional gravity data to obtain the density model.

Images