Fracture prediction method for narrow-azimuth and low-coverage seismic data acquisition of land data

Abstract

The invention provides a crack prediction method for acquiring seismic data in a narrow-azimuth and low-coverage manner in allusion to land data. The method mainly solves the problems that the fracture prediction difficulty is large and the prediction result is uncertain due to the fact that ellipse fitting cannot be achieved for seismic data collected in a narrow direction in the prior art. The method comprises the following steps: S1, establishing a floating reference surface model of a work area time domain; s2, correcting the time domain CMP gather to a small smooth surface; s3, preparing and loading data; s4, filling in main offset parameters; s5, omnibearing local angle domain pre-stack depth migration imaging is carried out; s6, azimuth anisotropy correction processing is carried out on the reflection angle gather, and a processed reflection angle gather is generated; s7, dividing the azimuth angle group of the corrected reflection angle gather to obtain reflection angle gather data of the divided azimuth angle group; and S8, pre-stack crack prediction. According to the crack prediction method, on the basis of existing acquired data, the crack prediction precision can be essentially improved.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic data processing and imaging, and in particular to a crack prediction method for seismic data collected in narrow azimuth and low coverage on land data. Background Art

[0002] Fracture oil and gas reservoir exploration technology is of great significance to the development of the petroleum industry. The geological causes of fractures are complex, and their development is controlled by many factors. They have the characteristics of complex distribution, poor regularity, and difficulty in prediction and identification. At present, the pre-stack and post-stack fracture prediction theory and technical methods based on high-density and wide-azimuth seismic data have been perfected and widely used. However, due to the constraints of the actual acquisition environment and exploration accuracy, the data in some areas are collected at narrow azimuths and low coverage times. For seismic data collected at narrow azimuths, ellipse fitting cannot be achieved, which increases the difficulty of fracture prediction. Therefore, a complete technical system is needed to solve the problem of fracture prediction for data collected at narrow azimuths and low coverage times. Summary of the invention

[0003] The technical problem to be solved by the present invention is to overcome the problem that the existing methods in the background technology cannot achieve ellipse fitting for seismic data collected in narrow azimuths, which increases the difficulty of fracture prediction and the uncertainty of the prediction results, and provide a fracture prediction method for seismic data collected in narrow azimuths and low coverage of land data. The fracture prediction method for seismic data collected in narrow azimuths and low coverage of land data can substantially improve the accuracy of fracture prediction on the basis of existing collected data, and the seismic information contained is more comprehensive, providing more accurate basic data for fracture prediction.

[0004] The present invention solves the problem through the following technical solution: The fracture prediction method for narrow-azimuth and low-coverage seismic data collected from land data includes the following steps:

[0005] S1: Establish a floating datum model in the time domain of the work area;

[0006] S2: Correct the time domain CMP gathers to a small smooth surface;

[0007] S3: Data preparation and loading for all-round local angle domain prestack depth migration imaging;

[0008] S4: Fill in the main migration parameters into the migration imaging system;

[0009] S5: All-round local angle domain pre-stack depth migration imaging to generate dip gathers and reflection angle gathers;

[0010] S6: Perform azimuth anisotropy correction processing on the reflection angle gathers obtained in S5 to generate reflection angle gathers after azimuth anisotropy correction processing;

[0011] S7: Divide the azimuth angle groups of the reflection angle gather after azimuth anisotropy correction to obtain the reflection angle gather data with divided azimuth angle groups:

[0012] S8: Prestack fracture prediction.

[0013] Further, the method of establishing a floating reference surface model in the time domain of a work area in S1 includes:

[0014] Smooth the near-surface elevation value to create a small smooth surface;

[0015] Then, a small smooth surface is used to establish a floating reference surface model in the depth domain;

[0016] Then, the replacement speed of the corresponding block is applied to perform time-depth conversion, and the floating reference surface model in the depth domain is converted to the time domain to obtain the floating reference surface model in the time domain.

[0017] Furthermore, the method S2 for correcting the time domain CMP gather to a small smooth surface includes:

[0018] Load the pre-stack CMP gathers in the time domain into the migration imaging system;

[0019] Then, the time domain floating datum model established in step S1 is applied to update the floating datum trace header information in the time domain CMP gather, and the time domain CMP gather is corrected to the small smooth surface.

[0020] Furthermore, the data preparation and loading of the S3 omnidirectional local angle domain prestack depth migration imaging includes:

[0021] The small smooth surface established by S1 is used as the starting surface for prestack depth migration in the omnidirectional local angle domain;

[0022] Load the pre-stack depth migration velocity volume and TTI anisotropy parameters into the migration imaging system;

[0023] The TTI anisotropy parameters include tilt volume, azimuth volume, δ and ε.

[0024] Furthermore, the main offset parameters of S4 include offset aperture, anti-aliasing parameters, offset step size, shallow half-opening angle and deep half-opening angle.

[0025] Furthermore, the specific method of the S5 all-round local angle domain prestack depth migration imaging includes:

[0026] Perform one-dimensional ray tracing on each diffraction point of the CMP gather that has been corrected to a small smooth surface to form a fan-shaped ray area;

[0027] Ray pairs are formed by this fan-shaped ray region;

[0028] The ray pairs simulate the event axes of reflection and scattering, and each ray pair is used to form a ray beam within the Fresnel zone;

[0029] Then the information of all shot points and receiving points within the beam diameter is synthesized into a seismic trace to generate a de-migrated image and a dip gather and a reflection angle gather.

[0030] Furthermore, the method of performing azimuthal anisotropy correction processing on the obtained reflection angle gather in S6 comprises the following steps:

[0031] 61. Convert the depth domain reflection angle gathers obtained by pre-stack depth migration in the omni-directional local angle domain to the time domain;

[0032] 62. Perform excision and dynamic balancing on the time domain reflection angle gathers; perform stacking and denoising on the time domain reflection angle gathers and use them as external model channels;

[0033] Input the time domain reflection angle gathers and external model traces after excision and dynamic balancing to the seismic data processing system. Perform cross-correlation processing on each trace in the reflection angle gathers and the input external model traces within the defined time window. According to the maximum cross-correlation coefficient principle, obtain the relative time shift in each cross-correlation time window point, that is, the residual time difference of the midpoint of the time window. Then perform interpolation to obtain the residual time difference and correlation coefficient data of the entire trace data.

[0034] 63. Apply the correlation coefficient data generated in the previous step and the picked residual time difference data, input them into the prestack time migration velocity field as a reference velocity field, estimate the azimuthal anisotropy parameters, and generate the azimuthal anisotropy parameters; the azimuthal anisotropy parameters include fast wave velocity, slow wave velocity and slow wave velocity azimuth;

[0035] 64. Input the reflection angle gather, reference velocity field and the azimuthal anisotropy parameters generated previously into the seismic data processing system, perform azimuthal anisotropy residual time difference correction processing, and generate reflection angle gathers that have been corrected for azimuthal anisotropy.

[0036] Furthermore, in step S7, the scheme for dividing the azimuth angle groups of the reflection angle gathers after azimuth anisotropy correction is as follows:

[0037] 71. Input the reflection angle gathers processed by azimuth anisotropy correction into the crack prediction system and conduct azimuth group division test;

[0038] 72. According to the principle of ensuring that the number of coverages is as consistent as possible, when the number of coverages of an azimuth group is higher than that of other azimuth groups, its range is adjusted to reduce the range of the azimuth group; when the number of coverages of an azimuth group is lower than that of other azimuth groups, the range of this azimuth group is expanded. The division of the entire azimuth group presents a fan-shaped distribution rather than an equally divided symmetrical distribution, thereby obtaining the reflection angle gather data of the divided azimuth groups.

[0039] Furthermore, the step S8 of pre-stack crack prediction includes the following steps:

[0040] 81. Input the HTI medium model into the crack prediction system;

[0041] 82. Based on the HTI medium model, azimuth attribute calculation is performed: by applying the pre-stack azimuth gathers that have been divided into azimuth groups, five pre-stack azimuth attributes are calculated, and the five pre-stack azimuth attributes are: seismic energy attribute, seismic frequency attribute, percentage of specified frequency energy in total energy, frequency attenuation gradient and seismic wave impedance attribute;

[0042] 83. Azimuth attribute analysis and optimization: By applying block logging data and rock physics data, the geological model and rock physics model of the fracture reservoir at the key well location are established. The full wave simulation technology is applied to calculate the seismic response of seismic waves, and the variation characteristics of the seismic attributes caused by the fractures with the offset and azimuth are analyzed, so as to determine the seismic attribute body that can be used for fracture prediction;

[0043] 84. Orientation attribute ellipse fitting to obtain ellipse flattening:

[0044] Furthermore, the method for obtaining the ellipse flattening by azimuth attribute ellipse fitting includes:

[0045] 101. The cracks under HTI medium conditions are high-angle cracks with a crack inclination angle close to 90 degrees. When the longitudinal wave incident angle is small, the relationship between the longitudinal wave reflection amplitude and the crack direction in the anisotropic medium can be expressed as the Aki-Richards equation simplified by the Zoplitz equation:

[0046]

[0047] Where: R pp is the longitudinal wave reflection amplitude, G is the rate of change of amplitude with offset, P is the reflection amplitude when the longitudinal wave is vertically incident (constant); G iso is the rate of change of amplitude with offset, i.e., isotropic gradient; G ani is the rate of change of amplitude with azimuth, i.e., the anisotropy gradient; is the crack direction; θ is the incident angle of the seismic wave; is the azimuth from the shot point to the receiver point;

[0048] 102. Through the analysis of the above equation (1), when the incident angle θ is 0°, the reflection coefficient of the fractured formation is independent of the azimuth, indicating that the close offset distance cannot reflect the anisotropy of the fractured medium; as the incident angle increases, the relative change of the reflection coefficient at each azimuth increases, and the anisotropic characteristics of the underground medium can be observed more; at the azimuth When the angle is 0°, that is, vertical to the crack direction, the P-wave reflection coefficient changes most significantly with the incident angle. When the angle is 90°, that is, parallel to the crack direction, the P-wave reflection coefficient changes the least with the incident angle;

[0049] It can be seen that the sensitivity of seismic waves at different incident angles to cracks is different. Select a certain incident angle range with greater sensitivity to cracks. Usually, a smaller incident angle range is selected. Within this range, sin 2 θ does not change much and can be regarded as a constant C. The above formula (1) can be simplified to:

[0050]

[0051] Where: A=P+G iso C, B = G ani C;

[0052] 103. In the crack prediction system, by applying the wave impedance attribute body of different azimuth groups and performing plane fitting, an ellipse can be obtained. A and B can be solved by applying formula (2); A represents the minor axis of the ellipse, B represents the major axis of the ellipse, and the ellipse flattening B / A is used as a relative measure of the crack density to predict the crack density. When B = 0, the medium is an isotropic medium, that is, it is considered that there are no cracks, and the P-wave reflection coefficient is fitted as a circle with a radius of A; the larger the absolute value of B, the flatter the shape of the ellipse, indicating that the more developed the cracks are, the greater the crack density.

[0053] Compared with the above-mentioned background technology, the present invention can have the following beneficial effects: the core of the fracture prediction method for narrow azimuth and low coverage seismic data collected from land data is to perform all-round local angle domain pre-stack depth migration processing on the seismic data collected in narrow azimuth and low coverage before fracture prediction. This migration method starts from the real reflection point underground and performs ray scanning and integration. It is a truly all-round processing method. Compared with the traditional OVT domain pre-stack time migration method, the imaging accuracy is higher and the seismic information contained is more comprehensive, providing more accurate basic data for fracture prediction. Combined with the non-equal azimuth group division method, it can eliminate the interference of coverage times on energy caused by the traditional equal azimuth group division method, and finally obtain a more accurate fracture density prediction result. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 The present invention discloses a flow chart of a fracture prediction method for seismic data collected in narrow azimuth and low coverage on land;

[0055] Figure 2 It is a reflection angle gather obtained by OVT domain migration gather and omnidirectional local angle domain prestack depth migration in an embodiment of the present invention (the left figure is an OVT domain migration gather; the right figure is an omnidirectional local angle domain migration gather);

[0056] Figure 3 1. The rose diagram of the OVT domain data of the variable view position and the rose diagram of the omnidirectional local angle domain data of the embodiment of the present invention (the left figure is the rose diagram of the OVT domain data; the right figure is the rose diagram of the omnidirectional local angle domain data);

[0057] Figure 4 It is the reflection angle gather of the omnidirectional local angle domain prestack depth migration before and after the azimuth anisotropy correction of the embodiment of the present invention (the left picture is before correction; the right picture is after correction);

[0058] Figure 5 These are the crack prediction results of OVT domain data and the crack prediction results of omni-directional local angle domain prestack depth migration in an embodiment of the present invention (the left figure is the crack prediction results of OVT domain data; the right figure is the crack prediction results of omni-directional local angle domain prestack depth migration). Specific implementation method:

[0059] The present invention will be further described below in conjunction with the accompanying drawings and embodiments:

[0060] The invention discloses a fracture prediction method for seismic data collected in narrow azimuth and low coverage on land data.

[0061] The following steps are involved:

[0062] S1: Establish a floating datum model in the time domain of the work area: Smooth the near-surface elevation value to establish a small smooth surface, and then use the small smooth surface to establish a floating datum model in the depth domain. Then, apply the replacement speed of the corresponding block to perform time-to-depth conversion, and convert the floating datum model in the depth domain to the time domain to obtain a floating datum model in the time domain.

[0063] S2: Time domain CMP gathers corrected to a small smooth surface;

[0064] The pre-stack CMP gathers in the time domain are loaded into the migration imaging system, and then the floating reference surface model in the time domain established in step S1 is applied to update the floating reference surface header information in the time domain CMP gathers, and the time domain CMP gathers are corrected to the small smooth surface;

[0065] S3: Data preparation and loading for all-round local angle domain prestack depth migration imaging: Use the small smooth surface established in S1 as the starting surface for all-round local angle domain prestack depth migration; load the prestack depth migration velocity body and TTI anisotropy parameters into the migration imaging system; the TTI anisotropy parameters include dip body, azimuth body, δ and ε.

[0066] S4: Fill in the main migration parameters into the migration imaging system; the main migration parameters include the migration aperture, anti-aliasing parameters, migration step size, shallow half-opening angle and deep half-opening angle.

[0067] S5: All-round local angle domain pre-stack depth migration imaging: One-dimensional ray tracing is performed on each diffraction point of the CMP gather that has been corrected to a small smooth surface to form a fan-shaped ray area; ray pairs are formed from this fan-shaped ray area; ray pairs simulate the phase axes of reflection and scattering, and each ray pair is formed into a ray beam within the Fresnel zone. Then, the information of all shot points and receiving points within the ray beam diameter is synthesized into a seismic trace de-migration imaging to generate dip gathers and reflection angle gathers.

[0068] S6: Perform azimuth anisotropy correction on the reflection angle gathers obtained in S5.

[0069] (1) Convert the depth domain reflection angle gathers obtained by prestack depth migration in the omni-directional local angle domain to the time domain;

[0070] (2) Perform excision and dynamic balancing on the time domain reflection angle gathers; stack and denoise the time domain reflection angle gathers and use them as external model tracks;

[0071] Input the time domain reflection angle gathers and external model traces after excision and dynamic balancing to the seismic data processing system. Perform cross-correlation processing on each trace in the reflection angle gathers and the input external model traces within the defined time window. According to the maximum cross-correlation coefficient principle, obtain the relative time shift in each cross-correlation time window point, that is, the residual time difference of the midpoint of the time window. Then perform interpolation to obtain the residual time difference and correlation coefficient data of the entire trace data.

[0072] (3) Apply the correlation coefficient data generated in the previous step and the picked residual time difference data to the prestack time migration velocity field as the reference velocity field; estimate the azimuthal anisotropy parameters to generate azimuthal anisotropy parameters; the azimuthal anisotropy parameters include fast wave velocity, slow wave velocity and slow wave velocity azimuth.

[0073] (4) Input the reflection angle gather, the reference velocity field and the azimuthal anisotropy parameter generated previously into the seismic data processing system, perform azimuthal anisotropy residual moveout correction processing, and generate a reflection angle gather that has been corrected for azimuthal anisotropy.

[0074] S7: dividing the azimuth angle groups of the reflection angle gather after azimuth anisotropy correction to obtain reflection angle gather data with divided azimuth angle groups;

[0075] In order to analyze the anisotropic response in different azimuths and obtain more accurate fracture prediction results, the reflection angle gathers that have been corrected for azimuth anisotropy need to be divided into azimuths.

[0076] The reflection angle gathers processed by azimuth anisotropy correction are input into the crack prediction system to carry out the azimuth group division test. According to the principle of ensuring the coverage times are as consistent as possible, when the coverage times of one azimuth group are higher than those of other azimuth groups, its range is adjusted to narrow the range of the azimuth group. When the coverage times of one azimuth group are lower than those of other azimuth groups, the range of this azimuth group is expanded. The division of the entire azimuth group presents a fan-shaped distribution instead of an equally divided symmetrical distribution, thereby obtaining the reflection angle gather data of the divided azimuth groups.

[0077] S8: Pre-stack fracture prediction:

[0078] (1) Input the HTI medium model into the crack prediction system;

[0079] (2) Based on the HTI medium model, azimuth attribute calculation is performed: by applying the pre-stack azimuth gathers that have been divided into azimuth groups, five pre-stack azimuth attributes are calculated, namely, seismic energy attribute, seismic frequency attribute, percentage of specified frequency energy in total energy, frequency attenuation gradient, and seismic wave impedance attribute;

[0080] (3) Azimuth attribute analysis and optimization: By applying block logging data and rock physics data, the geological model and rock physics model of the fracture reservoir at the key well location are established. The full wave simulation technology is used to calculate the seismic response of seismic waves and analyze the variation characteristics of the seismic attributes caused by the fractures with the offset and azimuth, so as to determine the seismic attribute body that can be used for fracture prediction;

[0081] (4) Calculate the ellipse flattening by fitting the orientation attribute ellipse:

[0082] The cracks under HTI medium conditions are high-angle cracks with a crack inclination angle close to 90 degrees. When the longitudinal wave incident angle is small, the relationship between the longitudinal wave reflection amplitude and the crack parameters (crack direction) in the anisotropic medium can be expressed as the Aki-Richards equation simplified by the Zoplitz equation:

[0083]

[0084] Where: R pp is the longitudinal wave reflection amplitude, G is the rate of change of amplitude with offset, P is the reflection amplitude when the longitudinal wave is vertically incident (constant); Giso is the rate of change of amplitude with offset, i.e., isotropic gradient; G ani is the rate of change of amplitude with azimuth, i.e., the anisotropy gradient; is the crack direction; θ is the incident angle of the seismic wave; is the azimuth from the shot point to the receiver point.

[0085] Through the above equation analysis, when the incident angle θ is 0°, the reflection coefficient is independent of the azimuth, indicating that the near offset distance cannot reflect the anisotropy of the fracture medium. As the incident angle increases, the relative change of the reflection coefficient at each azimuth increases, and the anisotropic characteristics of the underground medium can be observed. When the angle is 0°, that is, vertical to the crack direction, the P-wave reflection coefficient changes most significantly with the incident angle. When the angle is 90°, that is, parallel to the crack direction, the P-wave reflection coefficient changes the least with the incident angle.

[0086] It can be seen that the sensitivity of seismic waves at different incident angles to cracks is different. Select a certain incident angle range with greater sensitivity to cracks (usually a smaller incident angle range is selected, within which sin 2 θ does not change much and can be regarded as a constant C), and the above formula can be simplified to:

[0087]

[0088] Where: A=P+G iso C, B = G ani C.

[0089] According to the wave impedance attribute values ​​of different azimuth groups at any point, a plane fitting is performed to obtain an ellipse, and A and B in the formula can be solved (A represents the minor axis of the ellipse, and B represents the major axis of the ellipse). When B = 0, the medium is an isotropic medium, that is, when it is assumed that there are no cracks, the P-wave reflection coefficient is fitted as a circle with a radius of A. The larger the absolute value of B, the more developed the cracks, and the flatter the shape of the ellipse, so we can use the circular flattening ratio B / A as a relative measure of the crack density to obtain the crack density.

[0090] Example 1

[0091] The present invention will be further described below in conjunction with the accompanying drawings and specific embodiments:

[0092] Basic information of the application block: The acquisition surface element of the work area is 20x20m, the detection line distance is 400m, the shot line distance is 480m, the coverage times are 80 times, the observation system has serious local changes, the shot line distance and the detection line distance are large, the acquisition coverage times are low, the maximum offset distance is between 5000-6000m, the minimum offset distance is partially missing, and the aspect ratio of the target layer Xujiahe Formation is less than 0.5.

[0093] according to Figure 1 The specific implementation steps are as follows:

[0094] Step S1: Establish a floating datum model for the work area: Smooth the near-surface elevation value according to the 500m smoothing parameter to establish a small smoothing surface, then use the small smoothing surface to establish a floating datum model in the depth domain, and then use the replacement speed of the corresponding block to perform time-to-depth conversion, convert the floating datum model in the depth domain to the time domain, and obtain a floating datum model in the time domain.

[0095] Step S2: Data preparation: Load the pre-stack CMP gathers in the time domain, apply the time domain floating datum model established in step S1 to update the floating datum trace header information in the time domain CMP gathers, and correct the time domain CMP gathers to the small smooth surface. Use the small smooth surface as the starting surface of pre-stack depth migration for migration imaging. Load the pre-stack depth migration velocity body and TTI anisotropy parameters; the TTI anisotropy parameters include dip body, azimuth body, δ and ε.

[0096] Step S3: Apply step S2 to complete the loaded data, perform all-round local angle domain pre-stack depth migration imaging, and define the main migration parameters including: migration aperture 10000mx10000m, anti-aliasing parameter 0.75, migration step length 5m, shallow half-opening angle 60°, deep half-opening angle 45°, output reflection angle gathers and dip gathers, and data coverage times 250 times.

[0097] Step S4: Apply the omnidirectional local angle domain prestack depth migration imaging in step S3 to generate omnidirectional reflection angle gathers. The comparison between the conventional OVT domain prestack time migration gathers and the reflection angle gathers obtained by omnidirectional local angle domain prestack depth migration imaging after conversion to time domain is shown in Figure 4. Figure 2 , the azimuthal anisotropy characteristics of the reflection angle gather are more obvious. Figure 3 For a certain variable observation position, a rose diagram is drawn using the offset and azimuth information of the OVT domain prestack time migration gather and the omnidirectional local angle domain prestack depth migration reflection angle gather. As shown in the figure, the offset and azimuth information of the omnidirectional reflection angle gather is more completely preserved, and the aspect ratio is close to 1:1, providing high-precision prestack seismic data for subsequent fracture prediction.

[0098] Step S5: To eliminate the influence of azimuthal anisotropy moveout, azimuthal anisotropy correction is performed on the reflection angle gathers generated by omnidirectional local angle domain migration.

[0099] (1) Input reflection angle gathers and prestack depth migration velocity volume, and convert the reflection angle gathers in the depth domain to the time domain;

[0100] (2) Perform excision and dynamic balancing on the time domain reflection angle gathers; stack and denoise the time domain reflection angle gathers and use them as external model tracks;

[0101] Input the time domain reflection angle gathers and external model traces after resection and dynamic balancing to the seismic data processing system, define the maximum static correction allowable value as 10ms, the start and end time of the processing time window as 500-2400ms, the rolling length of the cross-correlation time window as 100ms, the cross-correlation time window length as 120ms, the transition zone at the upper and lower boundaries of the processing time window as 20ms, and output the residual time difference data and correlation coefficient data;

[0102] (3) Input the prestack time migration velocity field, residual time difference data and correlation coefficient data. The correlation coefficient threshold is defined as 0.1, the power of the correlation coefficient is defined as 2, and the maximum ellipticity of the allowed azimuthal velocity ellipse is defined as 1.3, which is the ratio of the fast velocity domain to the slow velocity. The azimuthal anisotropy parameters are estimated, including fast, slow and slow azimuth.

[0103] (4) Input azimuth gathers, prestack time migration velocity field and azimuth anisotropy parameters, and perform azimuth anisotropy residual moveout correction. Define the maximum imaging reference velocity error percentage allowed by residual azimuth moveout correction as 0.5, the maximum allowable ratio of fast velocity to slow velocity as 1.06, and define the percentage of the maximum passband frequency of the low-pass interpolation filter to the Nyquist frequency as 60 during moveout correction. Finally, generate reflection angle gathers after azimuth anisotropy correction, see Figure 4 .

[0104] Step S6: In order to eliminate the influence of the number of coverages on energy, a non-equal azimuth angle scheme is adopted to make the number of coverages of the target layer uniform. The reflection angle gathers corrected for azimuth anisotropy in step S5 are loaded into the fracture prediction system, and the azimuth angles are divided according to the principle of uniform coverage. The specific division scheme is: -180° to -135°, -135° to -89°, -89° to -30°, -30° to 23°, 23° to 79°, 79° to 130°, 130° to 180°.

[0105] Step S7: Pre-stack fracture prediction:

[0106] (1) Input the HTI medium model into the crack prediction system;

[0107] (2) Azimuth attribute calculation: After determining that the azimuth group division is reasonable, the reflection angle data of the divided azimuth groups are used to calculate the azimuth attributes, which include five pre-stack azimuth seismic attributes, namely, seismic energy attribute, seismic frequency attribute, percentage of specified frequency energy in total energy, frequency attenuation gradient, and seismic wave impedance attribute;

[0108] (3) Azimuth attribute analysis and optimization: Based on the well compliance rate, the seismic wave impedance attribute body is optimized and regarded as equivalent to the seismic amplitude calibration to analyze the change of amplitude with azimuth.

[0109] (4) Ellipse fitting to obtain the ellipse flattening: Through ellipse fitting, the spatial variation of the azimuthal ellipse of the amplitude is obtained, thereby predicting the density of the cracks.

[0110] like Figure 5 ,By comparing the fracture density prediction results with the traditional OVT domain data, it can be seen that the fracture density prediction has a higher consistency rate for the wells.

[0111] The fracture prediction method based on all-round local angle domain offset imaging data disclosed in the present invention can effectively solve the fracture prediction problem of seismic data with narrow azimuth, low coverage and serious observation system variation of land data. In the imaging process of this technical method, all wave field data within the offset aperture are used to form true dip and azimuth information about the geological body, the azimuth and offset distance are kept intact, the aspect ratio of the target layer is more ideal, and the anisotropic characteristics are more obvious. Combined with the non-equal azimuth group division method, it can eliminate the interference of uneven coverage on energy caused by the traditional equal azimuth group division method, and finally obtain a more accurate fracture density prediction result, which has broad application prospects.

[0112] It should be noted that the embodiments described above are only used to explain the present invention and do not constitute any limitation to the present invention. The present invention has been described with reference to typical embodiments, but it should be understood that the words used therein are descriptive and explanatory words, rather than restrictive words. The present invention may be modified as specified within the scope of the claims of the present invention, and the present invention may be revised without departing from the scope and spirit of the present invention. Therefore, the content of the present invention is not limited to the examples listed, and any equivalent changes taken by ordinary technicians in the field of the present invention by reading the specification of the present invention to the technical solution of the present invention are all covered by the claims of the present invention.

Claims

1. A fracture prediction method for narrow-azimuth, low-coverage seismic data collected on land, characterized by: The following steps are involved: S1: Establish a floating datum model in the time domain of the work area; S2: Based on the time domain floating datum model established in S1, the time domain CMP gathers are corrected to a small smooth surface; S3: Data preparation and loading for all-round local angle domain prestack depth migration imaging; S4: Fill in the main migration parameters into the migration imaging system; S5: All-round local angle domain pre-stack depth migration imaging to generate dip gathers and reflection angle gathers; S6: Perform azimuth anisotropy correction processing on the reflection angle gathers obtained in S5 to generate reflection angle gathers after azimuth anisotropy correction processing; S7: dividing the azimuth angle groups of the reflection angle gather after azimuth anisotropy correction to obtain reflection angle gather data with divided azimuth angle groups; S8: Perform pre-stack fracture prediction.

2. The fracture prediction method for narrow-azimuth, low-coverage seismic data collected from land data according to claim 1 is characterized by: The method for establishing a floating datum surface model in the time domain of a work area in S1 includes: Smooth the near-surface elevation value to create a small smooth surface; Then, a small smooth surface is used to establish a floating reference surface model in the depth domain; Then, the replacement speed of the corresponding block is applied to perform time-depth conversion, and the floating reference surface model in the depth domain is converted to the time domain to obtain the floating reference surface model in the time domain.

3. The fracture prediction method for narrow-azimuth, low-coverage seismic data collected from land data according to claim 1, characterized in that: The method S2 for correcting the time domain CMP gather to a small smooth surface comprises: Load the pre-stack CMP gathers in the time domain into the migration imaging system; Then, the time domain floating datum model established in step S1 is applied to update the floating datum trace header information in the time domain CMP gather, and the time domain CMP gather is corrected to the small smooth surface.

4. The fracture prediction method for narrow-azimuth, low-coverage seismic data acquisition of land data according to claim 1, characterized in that: The data preparation and loading of the S3 all-round local angle domain prestack depth migration imaging includes: The small smooth surface established by S1 is used as the starting surface for prestack depth migration in the omnidirectional local angle domain; Load the pre-stack depth migration velocity volume and TTI anisotropy parameters into the migration imaging system; The TTI anisotropy parameters include tilt volume, azimuth volume, δ and ε.

5. The fracture prediction method for narrow-azimuth, low-coverage seismic data collected from land data according to claim 1, characterized in that: The main migration parameters of S4 include migration aperture, anti-aliasing parameters, migration step size, shallow half-opening angle and deep half-opening angle.

6. The fracture prediction method for narrow-azimuth, low-coverage seismic data collected from land data according to claim 1, characterized in that: The specific method of the S5 all-round local angle domain prestack depth migration imaging includes: Perform one-dimensional ray tracing on each diffraction point of the CMP gather that has been corrected to a small smooth surface to form a fan-shaped ray area; Ray pairs are formed by this fan-shaped ray region; The ray pairs simulate the event axes of reflection and scattering, and each ray pair is used to form a ray beam within the Fresnel zone; Then the information of all shot points and receiving points within the beam diameter is synthesized into a seismic trace to generate a de-migrated image and a dip gather and a reflection angle gather.

7. The fracture prediction method for narrow-azimuth, low-coverage seismic data collected from land data according to claim 1, characterized in that: The method of performing azimuthal anisotropy correction processing on the obtained reflection angle gathers in S6 comprises the following steps:

61. Convert the depth domain reflection angle gathers obtained by pre-stack depth migration in the omni-directional local angle domain to the time domain; 62. Perform excision and dynamic balancing on the time domain reflection angle gathers; perform stacking and denoising on the time domain reflection angle gathers and use them as external model channels; Input the time domain reflection angle gathers and external model traces after excision and dynamic balancing to the seismic data processing system. Perform cross-correlation processing on each trace in the reflection angle gathers and the input external model traces within the defined time window. According to the maximum cross-correlation coefficient principle, obtain the relative time shift in each cross-correlation time window point, that is, the residual time difference of the midpoint of the time window. Then perform interpolation to obtain the residual time difference and correlation coefficient data of the entire trace data.

63. Apply the correlation coefficient data generated in the previous step and the picked residual time difference data, input them into the prestack time migration velocity field as a reference velocity field, estimate the azimuthal anisotropy parameters, and generate the azimuthal anisotropy parameters; the azimuthal anisotropy parameters include fast wave velocity, slow wave velocity and slow wave velocity azimuth; 64. Input the reflection angle gather, reference velocity field and the azimuthal anisotropy parameters generated previously into the seismic data processing system, perform azimuthal anisotropy residual time difference correction processing, and generate reflection angle gathers that have been corrected for azimuthal anisotropy.

8. The fracture prediction method for narrow-azimuth, low-coverage seismic data collected from land data according to claim 1, characterized in that: The scheme for dividing the azimuth angle groups of the reflection angle gathers after azimuth anisotropy correction in step S7 is:

71. Input the reflection angle gathers processed by azimuth anisotropy correction into the crack prediction system and conduct azimuth group division test; 72. According to the principle of ensuring that the number of coverages is as consistent as possible, when the number of coverages of an azimuth group is higher than that of other azimuth groups, its range is adjusted to reduce the range of the azimuth group; when the number of coverages of an azimuth group is lower than that of other azimuth groups, the range of this azimuth group is expanded. The division of the entire azimuth group presents a fan-shaped distribution rather than an equally divided symmetrical distribution, thereby obtaining the reflection angle gather data of the divided azimuth groups.

9. The fracture prediction method for narrow-azimuth, low-coverage seismic data collected from land data according to claim 1, characterized in that: The step S8, pre-stack crack prediction, comprises the following steps:

81. Input the HTI medium model into the crack prediction system; 82. Based on the HTI medium model, azimuth attribute calculation is performed: by applying the pre-stack azimuth gathers that have been divided into azimuth groups, five pre-stack azimuth attributes are calculated, and the five pre-stack azimuth attributes are: seismic energy attribute, seismic frequency attribute, percentage of specified frequency energy in total energy, frequency attenuation gradient and seismic wave impedance attribute; 83. Azimuth attribute analysis and optimization: By applying block logging data and rock physics data, the geological model and rock physics model of the fracture reservoir at the key well location are established. The full wave simulation technology is applied to calculate the seismic response of seismic waves, and the variation characteristics of the seismic attributes caused by the fractures with the offset and azimuth are analyzed, so as to determine the seismic attribute body that can be used for fracture prediction; 84. Calculate the ellipse flattening by fitting the orientation attribute ellipse.

10. The fracture prediction method for narrow-azimuth, low-coverage seismic data collected from land data according to claim 9, characterized in that: The method for obtaining the ellipse flattening by orientation attribute ellipse fitting comprises:

101. The cracks under HTI medium conditions are high-angle cracks with a crack inclination angle close to 90 degrees. When the longitudinal wave incident angle is small, the relationship between the longitudinal wave reflection amplitude and the crack direction in the anisotropic medium can be expressed as: Where: R pp is the longitudinal wave reflection amplitude; G is the rate of change of amplitude with offset; P is the reflection amplitude when the longitudinal wave is vertically incident (constant); G iso is the rate of change of amplitude with offset; i.e., isotropic gradient; G ani is the rate of change of amplitude with azimuth, i.e., the anisotropy gradient; is the crack direction; θ is the incident angle of the seismic wave; is the azimuth from the shot point to the receiver point; 102. Through the analysis of the above equation (1), when the incident angle θ is 0°, the reflection coefficient of the fractured formation is independent of the azimuth, indicating that the close offset distance cannot reflect the anisotropy of the fractured medium; as the incident angle increases, the relative change of the reflection coefficient at each azimuth increases, and the anisotropic characteristics of the underground medium can be observed more; at the azimuth When the angle is 0°, that is, vertical to the crack direction, the P-wave reflection coefficient changes most significantly with the incident angle. When the angle is 90°, that is, parallel to the crack direction, the P-wave reflection coefficient changes the least with the incident angle; It can be seen that the sensitivity of seismic waves at different incident angles to cracks is different. Select a certain incident angle range with greater sensitivity to cracks. Usually, a smaller incident angle range is selected. Within this range, sin 2 θ does not change much and can be regarded as a constant C. The above formula (1) can be simplified to: Where: A=P+G iso C, B = G ani C; 103. In the crack prediction system, by applying the wave impedance attribute body of different azimuth groups and performing plane fitting, an ellipse can be obtained. A and B can be solved by applying formula (2); A represents the minor axis of the ellipse, B represents the major axis of the ellipse, and the ellipse flattening B / A is used as a relative measure of the crack density to predict the crack density. When B = 0, the medium is an isotropic medium, that is, it is considered that there are no cracks, and the P-wave reflection coefficient is fitted as a circle with a radius of A; the larger the absolute value of B, the flatter the shape of the ellipse, indicating that the more developed the cracks are, the greater the crack density.

Images