Method, apparatus, device and storage medium for synthesizing seismic fault training data

Abstract

The present invention discloses a method, apparatus, device and storage medium for synthesizing seismic fault training data. The method includes: generating a horizontal layered velocity model; adding dip angles and fold structures to convert the horizontal layered velocity model into a velocity model containing dip angles and fold structures; adding fault structures to the velocity model containing dip angles and fold structures to obtain a fault velocity model with dip angles, folds and fault structures and corresponding fault labels; calculating a fault reflectivity model corresponding to the fault velocity model and adding random noise; establishing a point spread function library; convolving the fault reflectivity model with a randomly selected point spread function in the point spread function library to obtain seismic fault training data. By establishing a fault reflectivity model and a corresponding three-dimensional point spread function and performing convolution, the present invention can quickly synthesize real seismic fault data, providing a data basis for the quantitative interpretation of seismic faults and the intelligent interpretation of faults based on deep learning.

Description

Technical Field

[0001] The present invention relates to the field of intelligent interpretation of seismic faults, and particularly to a method, device, equipment and storage medium for synthesizing seismic fault training data. Background Art

[0002] Intelligent interpretation of seismic faults based on deep learning is one of the hot research topics in the field of exploration seismology. Among them, training data plays a crucial role in deep learning methods. Currently, there are mainly three common methods for obtaining fault training data: 1. Manually annotating actual seismic fault data to obtain labeled fault training data; 2. Establishing a fault model and obtaining corresponding seismic data through forward modeling and imaging; 3. Using code to establish a random fault model and artificially synthesizing fault training data through convolution of the model with a one-dimensional wavelet.

[0003] However, these methods all have certain limitations. For example, manually annotating three-dimensional seismic fault data is time-consuming, and the process is highly subjective; although forward modeling and imaging can ensure the authenticity and objectivity of the data, the process is very time-consuming and not suitable for obtaining a large amount of three-dimensional fault training data; the artificial synthesis method for obtaining fault training data is fast and has good results, and is currently the most commonly used method for obtaining seismic training data. However, the seismic data obtained by convolving with a one-dimensional wavelet has too high a resolution and has an obvious difference from actual seismic data, which will affect the effect of fault identification in practical applications.

[0004] Therefore, the technical problem that needs to be urgently solved is that the seismic data obtained by convolving with a one-dimensional wavelet has too high a resolution and has an obvious difference from actual seismic data. Summary of the Invention

[0005] In order to solve the technical problem that the seismic data obtained by convolving with a one-dimensional wavelet has too high a resolution and has an obvious difference from actual seismic data, the present invention proposes a method for synthesizing seismic fault training data based on three-dimensional convolution. This method can quickly synthesize real seismic fault data by convolving a fault seismic geological model with a corresponding three-dimensional point spread function, providing a data basis for quantitative interpretation of seismic faults and intelligent interpretation of faults based on deep learning.

[0006] According to a first aspect of the present invention, a method for synthesizing seismic fault training data includes the following steps:

[0007] Generate a horizontal layered velocity model;

[0008] Add dip and fold structures to convert the horizontal layered velocity model into a velocity model containing dip and fold structures;

[0009] On the velocity model containing dip angle and fold structure, add fault structure to obtain a fault velocity model with dip angle, fold and fault structures and the corresponding fault labels;

[0010] Calculate the fault reflectivity model corresponding to the fault velocity model and add random noise;

[0011] Establish a point spread function library;

[0012] Convolve the fault reflectivity model with a randomly selected point spread function in the point spread function library to obtain seismic fault training data.

[0013] Further, the generation of the horizontal layered velocity model includes:

[0014] Randomly generate a velocity sequence v(Z), where the magnitude and position of the velocity are both random, and the value distribution range is given according to the geological conditions of the actual work area. Then replicate and extend v(Z) horizontally to form a three-dimensional horizontal layered velocity model v(X, Y, Z).

[0015] Further, the size of the horizontal layered velocity model is 128×128×128.

[0016] Further, the addition of dip angle and fold structure to convert the horizontal layered velocity model into a velocity model containing dip angle and fold structure includes:

[0017] Add formation dip angle. The varying displacement field M1(X, Y, Z) is defined by a plane function, and the expression is:

[0018] M1(X, Y, Z) = aX + bY + c0

[0019] In the formula, a, b, and c0 are the coefficients of the plane function, which control the form of the horizontal formation. a and b control the dip angles in the X direction and Y direction respectively. To keep the middle trace (X C , Y C , Z) stationary, c0 = -aX C – bY C ;

[0020] Generate fold structure by adding a vertical displacement field M2(X, Y, Z). M2(X, Y, Z) is composed of N two-dimensional Gaussian functions and a linear attenuation function δ(Z), and the expression is:

[0021]

[0022] Among them, ± controls the type of fold. When it is +, the fold is an anticline; when it is -, the fold is a syncline. The kth Gaussian function is determined by the center position (c k , d k ), the half-length in the X direction σxk , the half length σ in the Y direction yk and the amplitude b k are defined; by summing N Gaussian functions, multiple different fold structures are generated; through the linear attenuation function δ(Z), the amplitude of the fold is attenuated in the vertical direction;

[0023] When the amplitude attenuates from the bottom to the top, the expression of δ(Z) is:

[0024]

[0025] When the amplitude attenuates from the top to the bottom, the expression of δ(Z) is:

[0026]

[0027] where Z and Z max are the depth and the maximum depth respectively. Through these two displacement fields M1(X, Y, Z) and M2(X, Y, Z) in the vertical direction, the horizontal layered velocity model v(X, Y, Z) is converted into a velocity model containing dip angles and fold structures satisfies:

[0028]

[0029] Furthermore, on the velocity model containing dip angles and fold structures, a fault structure is added to obtain a fault velocity model with dip angles, folds and fault structures and the corresponding fault labels, including:

[0030] a: Determine the fault reference plane:

[0031] Determine the fault reference plane defined by the reference point P(X0, Y0, Z0), the strike angle φ and the true dip angle θ, where the reference point P(X0, Y0, Z0) is the center point of the fault plane, φ ∈ [0°, 360°), the due north points to the negative direction of the Y axis, and θ ∈ [0°, 180°);

[0032] b: Coordinate transformation:

[0033] According to the reference point P, the strike angle φ and the true dip angle θ, the global coordinate system (X, Y, Z) is transformed into the local coordinate system (x, y, z), where the reference point P is the origin of the local coordinate system; the positive direction of the x-axis points to the strike; the positive direction of the y-axis points to the dip; the positive direction of the z-axis is the normal of the reference plane and points to the hanging wall part:

[0034]

[0035] where, R is the rotation matrix:

[0036]

[0037] When θ ∈ (0°, 90°), n = 2; when θ ∈ (90°, 180°), n = 1

[0038] c: Generating a fault surface by perturbation:

[0039] Randomly generate a series of perturbation points near the reference plane in the local coordinate system. According to these perturbation points, use bicubic spline interpolation to obtain the fault surface z f = f(x, y);

[0040] d: Calculating the fault displacement field:

[0041] Calculate the projected displacement field d(x, y; z = 0) of the fault displacement field D(x, y, z) on the fault reference plane, and use d(x, y; z = 0) to calculate the displacement field of each point;

[0042] Among them, d(x, y; z = 0) is defined by an ellipse formula, and its center point is located at the reference point P(X0, Y0, Z0), satisfying:

[0043]

[0044]

[0045] In the formula, d max is the maximum displacement; l x and l y are the semi - lengths of the fault surface in the x and y directions respectively;

[0046] Calculate the displacement of each point in the model using d(x, y; z = 0), that is, D(x, y, z). D(x, y, z) is a vector field, which is composed of the strike x - direction displacement field D x (x, y, z), the dip y - direction displacement field D y (x, y, z) and the z - direction displacement field D z (x, y, z), and is expressed as:

[0047] D(x, y, z) = (D x (x, y, z)u x , D y (x, y, z)u y , D z (x, y, z)u z )

[0048] Among them, u x , u y and u z are the unit vectors in the x, y and z directions respectively, and D x (x, y, z) = 0;

[0049] Define a boundary γ. When z reaches or exceeds the boundary γ, D y (x, y, z) = 0. γ is called the traction radius. For each point (x, y, z) in the hanging wall ((f(x, y) ≤ z ≤ γ + f(x, y))), D y The calculation formula of D

[0050] D y (x, y, z) = λ · d(x, y; z = 0) · α(x, y, z)

[0051] In the formula, λ represents the proportion of the displacement of the upper and lower plates in the projected displacement field d(x, y; z = 0), and its value range is between 0 and 1. When λ = 1, only the upper plate generates displacement and the lower plate does not move. When λ = 0, only the lower plate slides and the upper plate does not move. α(x, y, z) is a non - linear attenuation factor, which makes D y gradually decrease when (x, y, z) moves away from the fault surface along the normal direction. Its expression is as follows:

[0052]

[0053] For each point (x, y, z) in the footwall (f(x, y) - γ ≤ z ≤ f(x, y)), D y is:

[0054] D y (x, y, z) = (λ - 1) · d(x, y; z = 0) · α(x, y, z)

[0055] D z The calculation formula of D

[0056] D z (x, y, z) = f(x, y + D y (x, y, z)) - f(x, y)

[0057] e: Generating fault structures using the fault displacement field:

[0058] Add fault structures to the previously generated velocity model containing dip and fold structures using the calculated three - dimensional displacement field. It mainly consists of the following three steps:

[0059] Use the rotation matrix R to transform the coordinates in the global coordinate system to the local coordinate system (x, y, z):

[0060]

[0061] Apply the three - dimensional displacement field D(x, y, z) in the local coordinate system to simulate the generation of fault displacement:

[0062]

[0063] Convert the local coordinate system back to the global coordinate system:

[0064]

[0065] f: Repeat steps a - e for n times to realize adding multiple different layered fault structures to the velocity model with dip angle and fold structures, and obtain multiple fault velocity models with dip angle, fold and fault structures;

[0066] g: For each fault velocity model with dip angle, fold and fault structures, obtain the labels of the fault training data.

[0067] Furthermore, calculating the fault reflectivity model corresponding to the fault velocity model and adding random noise includes:

[0068] Set the density of each layer in the fault velocity model to be the same, calculate the fault reflectivity model corresponding to the velocity model according to the velocity of each layer in the velocity model, and add random noise to the fault reflectivity model. The calculation formula of the fault reflectivity model is:

[0069]

[0070] where v1 and v2 are the velocities of the upper layer and the lower layer respectively, and ρ1 and ρ2 are the densities of the upper and lower layers respectively.

[0071] Furthermore, establishing the point spread function library includes:

[0072] Use an observation system with a single source - receiver pair to detect a small target area U(x) near the subsurface x. The seismic signal emits from the source x s and then propagates to the target area and generates scattered waves; The seismic data recorded at the receiver x g is expressed as:

[0073]

[0074] In the formula, ω is the angular frequency; S(ω) is the source spectrum; k0 is the local wave number; x' is the local coordinate defined in the small target area U(x); V(x') = Δv(x') / v0(x) is the relative velocity perturbation, where V(x') is approximately the fault reflectivity model; G(x'; x s , ω) and G(x'; x g , ω) are the Green's functions from the source x s and the receiver x g to the target point x' respectively, and satisfy the cross - correlation imaging condition G(x'; xg , ω) = G(x g ; x', ω);

[0075] For an observation system with multiple source-receiver pairs, the seismic data I at x is expressed as:

[0076]

[0077] where * represents the complex conjugate operation, and G(x; x s , ω) and G(x; x g , ω) are the Green's functions from the source x s and the geophone x g to x, respectively; the broadband seismic data is expressed as:

[0078] I(x) = ∫ U(x) V(x')PSF(x, x')dx'

[0079] where the broadband point spread function PSF is defined as:

[0080]

[0081] According to the second aspect of the present invention, a seismic fault training data synthesis device includes the following modules:

[0082] A generation module for generating a horizontal layered velocity model;

[0083] A first addition module for adding dip and fold structures to convert the horizontal layered velocity model into a velocity model with dip and fold structures;

[0084] A second addition module for adding fault structures to the velocity model with dip and fold structures to obtain a fault velocity model with dip, fold, and fault structures and corresponding fault labels;

[0085] A calculation module for calculating the fault reflectivity model corresponding to the fault velocity model and adding random noise;

[0086] A building module for building a point spread function library;

[0087] A convolution module for convolving the fault reflectivity model with a random point spread function to obtain seismic fault training data.

[0088] According to the third aspect of the present invention, an electronic device includes a memory, a processor, and a computer program stored on the memory and executable on the processor, and when the processor executes the program, the steps of the seismic fault training data synthesis method are implemented.

[0089] According to a fourth aspect of the present invention, a storage medium stores a computer program thereon, and when the computer program is executed by a processor, the steps of the seismic fault training data synthesis method described above are implemented.

[0090] The technical solution provided by the present invention has the following beneficial effects:

[0091] 1. Since the present invention uses a three-dimensional point spread function PSF to perform three-dimensional spatial convolution with a randomly constructed fault reflectivity model, compared with the one-dimensional wavelet convolution method, this method is more in line with the seismic imaging theory. Therefore, the synthesized seismic fault training data is more realistic.

[0092] 2. The present invention uses an analytical method to calculate the Green's function and then calculate the point spread function. Compared with the numerical method, it can greatly improve the calculation efficiency and provide an efficiency guarantee for quickly obtaining training data.

[0093] 3. Since the application effect of the deep learning-based fault intelligent recognition network is largely affected by the training data, when the synthesized training data is highly similar to the actual seismic data, the application effect of the network will be better. The present invention can provide more realistic and reliable fault training data for deep learning. Therefore, the deep learning intelligent fault interpretation based on the present invention can obtain better application effects. BRIEF DESCRIPTION OF THE DRAWINGS

[0094] The present invention will be further described below with specific embodiments in conjunction with the drawings. In the drawings:

[0095] Figure 1 is a flowchart of the seismic fault training data synthesis method of the present invention;

[0096] Figure 2 is a flowchart of the fault training data synthesis of the present invention;

[0097] Figure 3 is a schematic diagram of the observation system of the present invention with a single shot gather;

[0098] Figure 4 is a partial example in the point spread function library PSF of the present invention;

[0099] Figure 5 is the seismic fault data obtained by convolving the fault reflectivity model of the present invention with the point spread function PSF;

[0100] Figure 6 is a schematic diagram of a partial synthesized three-dimensional fault sample and label of the present invention;

[0101] Figure 7 is a structural diagram of the seismic fault training data synthesis device of the present invention;

[0102] Figure 8It is a schematic structural diagram of an electronic device according to the present invention. Detailed implementation manners

[0103] For a clearer understanding of the technical features, objectives, and effects of the present invention, the detailed implementation manners of the present invention will now be described in detail with reference to the accompanying drawings.

[0104] Refer to Figure 1 , Figure 1 It is a flowchart of a method for synthesizing seismic fault training data according to the present invention. The method for synthesizing seismic fault training data based on three-dimensional convolution according to the present invention includes the following steps:

[0105] Step S1: Generate a horizontal layered velocity model.

[0106] First, a velocity sequence v(Z) is randomly generated. The magnitude and position of the velocity are both random, and the value distribution range is given according to the geological conditions of the actual work area. Then, v(Z) is replicated and extended horizontally to form a three-dimensional horizontal layered velocity model v(X, Y, Z). The size of the model can be set according to the actual situation. In the present invention, the model size is set to 128×128×128, as shown in Figure 2 (a).

[0107] Step S2: Add dip and fold structures.

[0108] Adding formation dip and fold structures to the horizontal layer is achieved by adding vertical displacement fields M1(X, Y, Z) and M2(X, Y, Z) to the horizontal layer.

[0109] First, add the formation dip. The varying displacement field M1(X, Y, Z) is defined by a plane function, as shown in Equation (1).

[0110] M1(X, Y, Z) = aX + bY + c0 (1)

[0111] In the formula, a, b, and c0 are the coefficients of the plane function, which control the form of the horizontal formation. a and b control the dip in the X direction and Y direction respectively; to keep the middle trace (X C , Y C , Z) stationary, c0 = -aX C – bY C .

[0112] Next, generate fold structures by adding the vertical displacement field M2(X, Y, Z). M2(X, Y, Z) is composed of N two-dimensional Gaussian functions and a linear attenuation function δ(Z), as shown in Equation (2).

[0113]

[0114] Among them, ± controls the type of fold. When it is +, the fold is an anticline; when it is -, the fold is a syncline. The k-th Gaussian function is defined by the center position (c k , d k ), the half-length σ xk in the X direction, the half-length σ yk in the Y direction, and the amplitude b k . To generate a relatively gentle fold structure, it is recommended to select a smaller amplitude value b k and larger half-lengths σ xk , σ yk . By summing N Gaussian functions, multiple different fold structures can be generated. Through the linear attenuation function δ(Z), the amplitude of the fold can be attenuated in the vertical direction (Z direction). When the amplitude attenuates from the bottom to the top, the expression of δ(Z) is:

[0115]

[0116] When the amplitude attenuates from the top to the bottom, δ(Z) is:

[0117]

[0118] where Z and Z max are the depth and the maximum depth respectively. Through these two displacement fields M1(X, Y, Z) and M2(X, Y, Z) in the vertical direction, the horizontal layered velocity model v(X, Y, Z) can be transformed into a velocity model containing dip angles and fold structures As Figure 2 (b) shows, where:

[0119]

[0120] Step S3: Add fault structures.

[0121] Fault structures are very important structures in geological structures. In the present invention, a three-dimensional fault displacement field is used to simulate and generate faults in the previous velocity model. The generation of fault structures mainly includes five steps: determining a fault reference plane, coordinate transformation, generating a fault surface by perturbation, calculating the fault displacement field, and generating a fault structure using the fault displacement field. The specific description of the generation of fault structures is as follows:

[0122] a: Determine the fault reference plane;

[0123] The first step in simulating the generation of fault structures is to determine a fault reference plane, which is defined by a reference point P(X0, Y0, Z0) (the center point of the fault plane), a strike angle φ (φ ∈ [0°, 360°), with the due north pointing in the negative Y-axis direction), and a true dip angle θ (θ ∈ [0°, 180°)).

[0124] The purposes of determining the fault reference plane are as follows:

[0125] (1) For the convenience of subsequent calculations, convert the global coordinate system (X, Y, Z) to the local coordinate system (x, y, z) with the reference point P(X0, Y0, Z0) as the origin;

[0126] (2) Since the actual fault plane is rarely an ideal plane but an uneven curved surface, perturbations can be applied to the fault reference plane to generate a fault surface;

[0127] (3) Determine the hanging wall and footwall of the fault to facilitate the calculation of the fault displacement fields of the upper and lower walls.

[0128] b: Coordinate transformation;

[0129] Using the reference point P, the strike angle φ, and the true dip angle θ, the global coordinate system (X, Y, Z) can be converted to the local coordinate system (x, y, z), where the reference point P is the origin of the local coordinate system; the positive direction of the x-axis points to the strike; the positive direction of the y-axis points to the dip; the positive direction of the z-axis is the normal of the reference plane and points to the hanging wall part.

[0130] According to the theoretical knowledge of three-dimensional space coordinate transformation, the global coordinate system can be converted to the local coordinate system using the following formula:

[0131]

[0132] where R is the rotation matrix:

[0133]

[0134] When θ ∈ (0°, 90°), n = 2; when θ ∈ (90°, 180°), n = 1

[0135] c: Generating the fault surface by perturbation;

[0136] To generate an uneven fault surface, a series of perturbation points are randomly generated near the reference plane in the local coordinate system. Taking these perturbation points as known points and using the biharmonic spline interpolation method to interpolate all (x, y) of the reference plane, the fault surface z f = f(x, y) can be obtained.

[0137] d: Calculating the fault displacement field;

[0138] Since the fault displacement field D(x, y, z) in the present invention is approximately an ellipsoid and the displacement amount of each point on it is different, which is not convenient for calculation. First, calculate the projected displacement field d(x, y; z = 0) of D(x, y, z) on the fault reference plane, and then use d(x, y; z = 0) to calculate the displacement field of each point.

[0139] In this embodiment, d(x, y; z = 0) is defined by an ellipse formula with its center point located at the reference point P(X0, Y0, Z0), and d(x, y; z = 0) is expressed as follows:

[0140]

[0141]

[0142] where d max is the maximum displacement; l x and l y are respectively the semi-lengths of the cross-section in the x and y directions.

[0143] Next, the displacement of each point in the model, i.e., D(x, y, z), is calculated using d(x, y; z = 0). D(x, y, z) is a vector field, which is composed of the displacement field D x (x, y, z) in the strike x direction, the displacement field D y (x, y, z) in the dip y direction, and the displacement field D z (x, y, z) in the z direction, and it can be expressed as:

[0144] D(x, y, z) = (D x (x, y, z)u x , D y (x, y, z)u y , D z (x, y, z)u z ) (8)

[0145] where u x , u y and u z are the unit vectors in the x, y, and z directions respectively.

[0146] In this embodiment, only dip-slip faults are considered and strike-slip faults are not considered. The displacement in the strike direction of the model is 0, i.e., D x (x, y, z) = 0.

[0147] Since D y (x, y, z) is parallel to the fault reference plane, the displacement field d(x, y; z = 0) calculated previously is used to calculate D y (x, y, z). Since D y (x, y, z) gradually decreases along the normal direction away from the fault surface, a boundary γ needs to be defined. When z reaches or exceeds the boundary γ, D y (x, y, z) = 0, and γ is called the traction radius.

[0148] For each point (x, y, z) on the hanging wall ((f(x, y) ≤ z ≤ γ + f(x, y))), D y (x, y, z) is calculated by the formula:

[0149] D y (x, y, z) = λ · d(x, y; z = 0) · α(x, y, z) (9)

[0150] In the formula, λ represents the ratio of the displacement of the upper and lower plates to the projection field d(x, y; z = 0), and its value range is between 0 and 1. When λ = 1, only the upper plate generates displacement and the lower plate does not move; when λ = 0, only the lower plate slides and the upper plate does not move. α(x, y, z) is a non-linear attenuation factor, which makes D y (x, y, z) gradually decrease when moving away from the fault surface along the normal direction. Its expression is as follows:

[0151]

[0152] For each point (x, y, z) on the footwall (f(x, y) - γ ≤ z ≤ f(x, y)), D y (x, y, z) is:

[0153] D y (x, y, z) = (λ - 1) · d(x, y; z = 0) · α(x, y, z) (11)

[0154] In this embodiment, only the case of seamless movement of the fault block along the fault surface is considered, that is, the points initially on the fault plane are still on the fault plane after the fault displacement, only the position on the fault plane has changed and will not break away from the fault plane. Therefore, the calculation formula of D z (x, y, z) is as follows:

[0155] D z (x, y, z) = f(x, y + D y (x, y, z)) - f(x, y) (12)

[0156] So far, the three components D x (x, y, z), D y (x, y, z) and D z (x, y, z) of the three-dimensional fault displacement field D(x, y, z) have all been calculated.

[0157] e: Generating fault structures using the fault displacement field

[0158] Adding fault structures to the previously generated fold model using the calculated three-dimensional displacement field, which mainly includes the following three steps:

[0159] (1) Use the rotation matrix R to transform the coordinates in the global coordinate system to the local coordinate system (x, y, z):

[0160]

[0161] (2) Apply the three-dimensional displacement field D(x, y, z) in the local coordinate system to simulate the generation of fault displacement:

[0162]

[0163] (3) Convert the local coordinate system back to the global coordinate system.

[0164]

[0165] Execute steps a - e, and then a fault structure can be added to the fold model. Repeatedly use this method on the same model multiple times, and multiple faults with different attitudes can be added to the model. As Figure 2 (c) shows, for each fault velocity model with dip, fold, and fault structures, the number and pattern of faults have been determined. Therefore, the labels of the fault training data can be obtained in advance at this stage.

[0166] Step S4, calculate the fault reflectivity model.

[0167] After obtaining the fault velocity model with dip, fold, and fault structures, calculate the fault reflectivity model according to formula (16). Assume that the density of each layer in the fault velocity model is the same. According to the velocity of each layer in the velocity model, the corresponding fault reflectivity model of the velocity model can be calculated. To make the synthesized data closer to the real data, random noise can be added to the reflectivity model.

[0168]

[0169] Among them, v1 and v2 are the velocities of the upper and lower layers respectively, and ρ1 and ρ2 are the densities of the upper and lower layers respectively.

[0170] Step S5, establish a point spread function library.

[0171] Currently, the most commonly used method for synthesizing training data is the convolution of the reflectivity model and the wavelet. However, the problem with this method is that the data resolution is too high, and it is difficult to be similar to the actual data even after adding noise. Taking fault identification as an example, this will result in the inability to accurately identify faults in some complex seismic data. Considering this problem, the present invention uses the convolution of the three-dimensional point spread function (PSF) to synthesize training samples. The relevant principle of PSF is as follows:

[0172] As Figure 3As shown in the figure, an observation system with a single source-receiver pair is used to detect a small target area U(x) near the subsurface x. The seismic signal is emitted from the source x s and then propagates to the target area, generating scattered waves. The seismic data recorded at the receiver x g can be expressed as:

[0173]

[0174] In Equation (17), ω is the angular frequency; S(ω) is the source spectrum; k0 is the local wavenumber; x' is the local coordinate defined in the target area U(x); V(x') = Δv(x') / v0(x) is the relative velocity perturbation, where V(x') can be approximated by the above fault reflectivity model; G(x'; x s , ω) and G(x'; x g , ω) are the Green's functions from the source x s and the receiver x g to the target point x', respectively, and satisfy the cross-correlation imaging condition G(x'; x g , ω) = G(x g ; x', ω).

[0175] For an observation system with multiple source-receiver pairs, the seismic data I at x can be expressed as:

[0176]

[0177] where * represents the complex conjugate operation. Substituting Equation (17) into Equation (18), the broadband seismic data can be expressed as:

[0178] I(x) = ∫ U(x) V(x')PSF(x, x')dx' (19)

[0179] where the broadband point spread function (PSF) is defined as:

[0180]

[0181] Therefore, according to Equation (19), seismic data can be considered as the convolution of the PSF and model parameters such as reflectivity. The PSF generally consists of four Green's functions, two of which consider the influence from the source, and the other two contain the acquisition aperture effect at the receiver end. To simulate accurate seismic data, the reflectivity model needs to be convolved with the PSF of the local area. Compared with the complete three-dimensional seismic data, a training sample of size 128×128×128 is a relatively small part. It can be assumed that the PSF varies little within this small range, that is, a reflectivity model is convolved with a single PSF instead of the spatially varying PSF to simulate a seismic data. By adopting this method, the computational cost of generating training samples can be significantly reduced. Although there is some loss in the imaging amplitude, the structural features in the seismic data are not greatly affected. As Figure 4 shown Figure 4 in,

[0182] Step S6: Three-dimensional convolution of the reflectivity model and the point spread function

[0183] Convolve the fault reflectivity model generated in Step 4 with a random PSF in the PSF library established in Step 5 ( Figure 5 ), and a large amount of seismic fault data can be obtained flexibly and stably. As Figure 6 shown,

[0184] The technical key points of the present invention are as follows:

[0185] 1. Establish an initial velocity model according to the seismic geological background of the actual work area, and form geological structures such as dips, folds, and faults. The fault training data is targeted;

[0186] 2. Use code to control the formation and fault parameters, and perform completely random modeling within a reasonable geological background to make the fault training data diverse;

[0187] 3. Use the three-dimensional point spread function (3D PSF) to convolve with the fault reflectivity model to synthesize more realistic seismic fault data, thereby improving the effect of intelligent interpretation of seismic faults using deep learning methods.

[0188] As Figure 7 shown, this embodiment also provides a device for synthesizing seismic fault training data for implementing the above-mentioned method for synthesizing seismic fault training data. The device specifically includes the following modules:

[0189] Generation module 1, used to generate a horizontal layered velocity model;

[0190] The first addition module 2 is used to add dip angle and fold structure, and convert the horizontal layered velocity model into a velocity model containing dip angle and fold structure;

[0191] The second addition module 3 is used to add fault structure to the velocity model containing dip angle and fold structure, and obtain a fault velocity model with dip angle, fold and fault structure and the corresponding fault label;

[0192] The calculation module 4 is used to calculate the fault reflectivity model corresponding to the fault velocity model and add random noise;

[0193] The establishment module 5 is used to establish a point spread function library;

[0194] The convolution module 6 is used to convolve the fault reflectivity model with a random point spread function to obtain seismic fault training data.

[0195] In addition, this embodiment also provides an electronic device, including a memory, a processor, and a computer program stored on the memory and executable on the processor. When the processor executes the program, it implements the steps of the above seismic fault training data synthesis method. For details, please refer to Figure 8 , which exemplifies a schematic diagram of the physical structure of an electronic device. The electronic device may include: a processor 610, a communication interface 620, a memory 630, and a communication bus 640. Among them, the processor 610, the communication interface 620, and the memory 630 complete mutual communication through the communication bus 640. The processor 610 can call the logical instructions in the memory 630 to execute the steps of the above seismic fault training data synthesis method, specifically including: generating a horizontal layered velocity model; adding dip angle and fold structure to convert the horizontal layered velocity model into a velocity model containing dip angle and fold structure; adding fault structure to the velocity model containing dip angle and fold structure to obtain a fault velocity model with dip angle, fold and fault structure and the corresponding fault label; calculating the fault reflectivity model corresponding to the fault velocity model and adding random noise; establishing a point spread function library; convolving the fault reflectivity model with a random point spread function in the point spread function library to obtain seismic fault training data.

[0196] In addition, when the logical instructions in the above-mentioned memory 630 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, in essence, or the part that contributes to the prior art, or a 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 can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods described in various embodiments of the present invention. The aforementioned storage medium includes: various media such as USB flash drives, mobile hard disks, read-only memories (ROM, Read-Only Memory), random access memories (RAM, Random Access Memory), magnetic disks, or optical discs that can store program codes.

[0197] In another aspect, an embodiment of the present invention also provides a storage medium, on which a computer program is stored. When the computer program is executed by a processor, it implements the steps of the above-mentioned seismic fault training data synthesis method, specifically including: generating a horizontal layered velocity model; adding dip angles and fold structures to convert the horizontal layered velocity model into a velocity model containing dip angles and fold structures; adding fault structures to the velocity model containing dip angles and fold structures to obtain a fault velocity model with dip angles, fold structures, and fault structures and corresponding fault labels; calculating the fault reflectivity model corresponding to the fault velocity model and adding random noise; establishing a point spread function library; convolving the fault reflectivity model with a random point spread function in the point spread function library to obtain seismic fault training data.

[0198] An embodiment of the present invention provides a seismic fault training data synthesis method, device, equipment, and storage medium. By generating a horizontal layered velocity model; adding dip angles and fold structures to convert the horizontal layered velocity model into a velocity model containing dip angles and fold structures; adding fault structures to the velocity model containing dip angles and fold structures to obtain a fault velocity model with dip angles, fold structures, and fault structures and corresponding fault labels; calculating the fault reflectivity model corresponding to the fault velocity model and adding random noise; establishing a point spread function library; convolving the fault reflectivity model with a random point spread function in the point spread function library to obtain seismic fault training data. The present invention can quickly synthesize real seismic fault data by establishing a fault reflectivity model and a corresponding three-dimensional point spread function and performing convolution, providing a data basis for the quantitative interpretation of seismic faults and the intelligent interpretation of faults based on deep learning.

[0199] It should be noted that in this text, the terms "include", "comprise" or any other variants thereof are intended to cover non-exclusive inclusion, such that a process, method, article or system comprising a series of elements not only includes those elements but also includes other elements not expressly listed, or further includes elements inherent to such process, method, article or system. Without further limitation, an element defined by the statement "comprising one..." does not exclude the presence of additional identical elements in the process, method, article or system comprising such element.

[0200] The serial numbers of the above embodiments of the present invention are for description only and do not represent the superiority or inferiority of the embodiments. In the unit claims listing several devices, several of these devices may be embodied by the same hardware item. The use of the words first, second, and third, etc. does not denote any order and these words may be construed as identifiers.

[0201] The above are only the preferred embodiments of the present invention, and do not limit the patent scope of the present invention accordingly. Any equivalent structure or equivalent process transformation made by using the content of the specification and drawings of the present invention, or directly or indirectly applied in other related technical fields, shall be equally included in the patent protection scope of the present invention.

Claims

1. A method for synthesizing seismic fault training data, characterized in that, It includes the following steps: Generate a horizontal layered velocity model; Add dip angle and fold structure, and convert the horizontal layered velocity model into a velocity model containing dip angle and fold structure; Add fault structure to the velocity model containing dip angle and fold structure to obtain a fault velocity model with dip angle, fold and fault structures and corresponding fault labels; Calculate the fault reflectivity model corresponding to the fault velocity model and add random noise; Establish a point spread function library; Convolve the fault reflectivity model with a random point spread function in the point spread function library to obtain seismic fault training data.

2. The method for synthesizing seismic fault training data according to claim 1, wherein, The generation of the horizontal layered velocity model includes: Randomly generate a velocity sequence v(Z), where the magnitude and position of the velocity are both random, and the value distribution range is given according to the geological conditions of the actual work area. Then, copy and extend v(Z) horizontally to form a three-dimensional horizontal layered velocity model v(X, Y, Z).

3. The method for synthesizing seismic fault training data according to claim 1, wherein The size of the horizontal layered velocity model is 128×128×128.

4. The method for synthesizing earthquake fault training data according to claim 1, wherein The addition of dip angle and fold structure to convert the horizontal layered velocity model into a velocity model containing dip angle and fold structure includes: Add a formation dip angle. The varying displacement field M1(X, Y, Z) is defined by a plane function, and the expression is: M1(X, Y, Z) = aX + bY + c0 Wherein, a, b, and c0 are the coefficients of the plane function, which control the shape of the horizontal formation. a and b respectively control the dips in the X and Y directions. To keep the middle trace (X C , Y C , Z) stationary, c0 = -aX C – bY C ; Generate a fold structure by adding a vertical displacement field M2(X, Y, Z). M2(X, Y, Z) is composed of N two-dimensional Gaussian functions and a linear attenuation function δ(Z), and the expression is: Among them, ± controls the type of fold. When it is +, the fold is an anticline; when it is -, the fold is a syncline. The k-th Gaussian function is determined by the center position (c k , d k ), the half-length σ xk in the X direction, the half-length σ yk in the Y direction, and the amplitude b k . By summing N Gaussian functions, multiple different fold structures are generated. Through the linear attenuation function δ(Z), the amplitude of the fold is attenuated in the vertical direction; When the amplitude decays from the bottom to the top, the expression of δ(Z) is: When the amplitude decays from the top to the bottom, the expression of δ(Z) is: where Z and Z max are the depth and the maximum depth respectively. Through the two displacement fields M1(X, Y, Z) and M2(X, Y, Z) in the vertical direction, the horizontal layered velocity model v(X, Y, Z) is converted into a velocity model containing dip angles and fold structures satisfy:

5. The method for synthesizing seismic fault training data according to claim 1, wherein The addition of fault structure to the velocity model containing dip angle and fold structure to obtain a fault velocity model with dip angle, fold and fault structures and corresponding fault labels includes: a: Determine the fault reference plane: Determine the fault reference plane defined by the reference point P(X0, Y0, Z0), strike angle φ and true dip angle θ, where the reference point P(X0, Y0, Z0) is the center point of the fault plane, φ ∈ [0°, 360°), the due north points to the negative direction of the Y axis, and θ ∈ [0°, 180°); b: Coordinate transformation: Convert the global coordinate system (X, Y, Z) to the local coordinate system (x, y, z) according to the reference point P, strike angle φ and true dip angle θ, where the reference point P is the origin of the local coordinate system; the positive direction of the x-axis points to the strike; the positive direction of the y-axis points to the dip; the positive direction of the z-axis is the normal of the reference plane and points to the hanging wall part: Among them, R is the rotation matrix: When θ ∈ (0°, 90°), n = 2; when θ ∈ (90°, 180°), n = 1 c: Perturb to generate the fault surface: A series of perturbation points are randomly generated near the reference plane in the local coordinate system. Based on these perturbation points, the fault surface z f = f(x, y); d: Calculate the fault displacement field: Calculate the projection displacement field d(x, y; z = 0) of the fault displacement field D(x, y, z) on the fault reference plane, and use d(x, y; z = 0) to calculate the displacement field of each point; Among them, d(x, y; z = 0) is defined by an ellipse formula, and its center point is located at the reference point P(X0, Y0, Z0), satisfying: where d max is the maximum displacement; l x and l y are the half lengths of the cross-section in the x- and y-directions, respectively; Calculate the displacement of each point in the model using d(x, y; z = 0), i.e., D(x, y, z). D(x, y, z) is a vector field, which is composed of the displacement field D x (x, y, z) in the strike x direction, the displacement field D y (x, y, z) in the dip y direction, and the displacement field D z (x, y, z) in the z direction, and is expressed as: D(x,y,z) = (D x (x,y,z)u x , D y (x,y,z)u y , D z (x,y,z)u z ) where u x , u y and u z are unit vectors in the x, y, and z directions respectively, and D x (x, y, z) = 0; Define a boundary γ. When z reaches or exceeds the boundary γ, D y (x, y, z) = 0, and γ is called the traction radius. For each point (x, y, z) on the upper plate ((f(x, y) ≤ z ≤ γ + f(x, y))), D y The calculation formula of (x, y, z) is as follows: D y (x,y,z) = λ·d(x,y; z = 0)·α(x,y,z) where λ represents the ratio of the displacement of the upper and lower plates to the projected displacement field d(x, y; z = 0), and its value ranges from 0 to 1. When λ = 1, only the upper plate moves and the lower plate remains stationary. When λ = 0, only the lower plate slides and the upper plate remains stationary. α(x, y, z) is a non-linear attenuation factor that causes D y (x, y, z) to gradually decrease as it moves away from the fault surface along the normal direction, and its expression is as follows: For each point (x, y, z) on the lower surface (f(x, y) - γ ≤ z ≤ f(x, y)), D y (x, y, z) is: D y (x, y, z) = (λ - 1)·d(x, y; z = 0)·α(x, y, z) D z (x, y, z) is calculated as follows: D z (x,y,z) = f(x,y + D y (x,y,z)) - f(x,y) e: Generate the fault structure using the fault displacement field: Adding fault structures to the previously generated velocity model containing dip angles and fold structures using the calculated three-dimensional displacement field mainly involves the following three steps: Use the rotation matrix R to transform the coordinates in the global coordinate system to the local coordinate system (x, y, z): Simulating fault displacements using the three-dimensional displacement field D(x, y, z) in the local coordinate system: Converting the local coordinate system back to the global coordinate system: f: Repeating steps a - e a total of n times to add multiple different layered fault structures to the velocity model containing dip angles and fold structures, obtaining multiple fault velocity models with dip angles, folds, and fault structures; g: For each fault velocity model with dip angles, folds, and fault structures, obtaining the labels of the fault training data.

6. The method for synthesizing seismic fault training data according to claim 1, wherein Calculating the fault reflectivity model corresponding to the fault velocity model and adding random noise, including: Setting the density of each layer in the fault velocity model to be the same, calculating the fault reflectivity model corresponding to the velocity model based on the velocity of each layer in the velocity model, and adding random noise to the fault reflectivity model. The calculation formula for the fault reflectivity model is: where v1 and v2 are the velocities of the upper and lower layers respectively, and ρ1 and ρ2 are the densities of the upper and lower layers respectively.

7. The method for synthesizing seismic fault training data according to claim 1, wherein Establishing the point spread function library, including: Use an observation system with a single source-receiver pair to detect a small target area U(x) near the subsurface x. The seismic signal is emitted from the source x s and then propagates to the target area and generates scattered waves; the seismic data recorded at the receiver x g is expressed as: where ω is the angular frequency; S(ω) is the source spectrum; k0 is the local wavenumber; x' is the local coordinate defined in the small target region U(x); V(x') = Δv(x') / v0(x) is the relative velocity perturbation, where V(x') is approximated as the fault reflectivity model; G(x'; x s , ω) and G(x'; x g , ω) are the Green's functions from the source x s and the geophone x g to the target point x', respectively, and satisfy the cross-correlation imaging condition G(x'; x g , ω) = G(x g ; x', ω); For an observation system containing multiple source-receiver pairs, the seismic data I located at x is expressed as: where * denotes the complex conjugate operation, G(x; x s , ω) and G(x; x g , ω) are the Green's functions from the source x s and the geophone x g to x, respectively; the broadband seismic data are obtained as follows: I(x) = ∫ U(x) V(x')PSF(x, x')dx' In the formula, the broadband point spread function PSF is defined as:

8. An earthquake fault training data synthesis device, characterized in that, Including the following modules: A generation module for generating a horizontally layered velocity model; A first addition module for adding dip angles and fold structures to convert the horizontally layered velocity model into a velocity model containing dip angles and fold structures; A second addition module for adding fault structures to the velocity model containing dip angles and fold structures to obtain a fault velocity model with dip angles, folds, and fault structures and the corresponding fault labels; A calculation module for calculating the fault reflectivity model corresponding to the fault velocity model and adding random noise; An establishment module for establishing the point spread function library; A convolution module for convolving the fault reflectivity model with a random point spread function to obtain seismic fault training data.

9. An electronic device, comprising a memory, a processor, and a computer program stored on the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the steps of the seismic fault training data synthesis method according to any one of claims 1 - 7.

10. A storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the steps of the seismic fault training data synthesis method according to any one of claims 1 - 7.

Images