Pre-stack wide azimuth anisotropy inversion method and device, storage medium and equipment

Abstract

The application provides a pre-stack wide azimuth anisotropy inversion method, device, storage medium and equipment, the method comprises the following steps: obtaining original pre-stack wide azimuth gather; performing azimuth AVO feature correction processing on the original pre-stack wide azimuth gather to obtain target pre-stack wide azimuth gather; sequentially reading a set number of seismic line data from the target pre-stack azimuth gather data to form a current inversion array, and performing anisotropy inversion based on the current inversion array to obtain a corresponding inversion result, until all the seismic line data in the target pre-stack azimuth gather data are obtained. Through azimuth AVO feature correction processing on the original data, the quality of the gather data is effectively improved, the problems of low pre-stack azimuth data quality and large data quantity are solved, and the subsequent inversion requirements are better met; through rolling extraction of seismic data to form a specific inversion array for step-by-step inversion, the calculation efficiency and stability of the inversion are improved.

Description

Technical Field

[0001] This invention relates to the field of oil and gas geophysical exploration, and in particular to a pre-stack wide-azimuth anisotropy inversion method, apparatus, storage medium and equipment. Background Technology

[0002] Due to tectonic stress, underground rock strata inevitably fracture, resulting in fractures. In oil and gas exploration, fractures not only serve as conduits for oil and gas migration but also as boundaries of fault-block oil and gas fields, effectively controlling their distribution. Fracture identification is fundamental to fracture research, and accurate fracture identification is crucial for oil and gas field exploration and development. The formation of fractures is controlled by multiple factors, exhibiting complex physical properties with significant lateral and longitudinal variations, demonstrating strong anisotropy. Therefore, obtaining stable anisotropic information is essential for detecting fracture characteristics.

[0003] Conventional fracture detection methods are broadly classified into post-stack and pre-stack methods. Post-stack fracture detection, based on post-stack seismic data, primarily identifies fractures by detecting discontinuities in seismic traces, such as coherence and similarity. This method has limited accuracy and struggles to identify small- to medium-scale fractures and cracks. Pre-stack fracture detection, based on pre-stack seismic data, mainly utilizes anisotropic information caused by variations in azimuth seismic amplitude and azimuth travel time for fracture prediction. Anisotropic inversion based on pre-stack azimuth data offers higher accuracy, and the fracture development intensity and direction obtained through inversion can better characterize fracture development features. However, current azimuth data is often of low quality due to various factors such as acquisition and processing, and the large data volume severely limits the application of pre-stack azimuth anisotropic inversion in fracture and crack detection. There is usually a lack of targeted correction methods for the original gathers. Directly inputting the original gathers for inversion will reduce the accuracy of the inversion. Conventional inversion algorithms also lack effective constraints. At the same time, the amount of data in the pre-stack gathers is large. Scanning the data head and recording the head information, and then locating the data volume by reading the head information file during the inversion process, requires frequent hard disk scanning, which results in slow speed. Therefore, most methods such as multi-trace stacking and reducing the order of the inversion formula are used to improve the efficiency and stability of the inversion results by sacrificing the accuracy of the inversion. Summary of the Invention

[0004] To address the aforementioned problems, embodiments of the present invention provide a pre-stack wide-azimuth anisotropy inversion method, apparatus, storage medium, and device.

[0005] In a first aspect, embodiments of the present invention provide a pre-stack wide-azimuth anisotropy inversion method, comprising:

[0006] Obtain the original pre-stack wide azimuth gather;

[0007] The original pre-stack wide azimuth gather is subjected to azimuth AVO feature correction processing to obtain the target pre-stack wide azimuth gather;

[0008] A set number of seismic survey lines are sequentially read from the target stack front-side gather data to form a current inversion array. Anisotropic inversion is then performed based on the current inversion array to obtain the corresponding inversion results, until all seismic survey lines in the target stack front-side gather data have obtained inversion results.

[0009] In some implementations, the step of performing azimuth AVO feature correction processing on the original pre-stack wide azimuth gather to obtain the target pre-stack wide azimuth gather includes:

[0010] Continuous time-varying static corrections are calculated for the original pre-stack wide azimuth gather, and the remaining travel time is estimated using the time-varying static corrections.

[0011] For each phase axis of the reflected wave, the remaining travel time is fitted using the residual NMO plane of the ellipse;

[0012] Calculate the superposition velocity and orientation of the ellipses;

[0013] The original pre-stack wide azimuth gather is corrected by applying the stacking velocity of the ellipse to obtain the target pre-stack wide azimuth gather.

[0014] In some implementations, the time-varying static correction is obtained through anisotropic residual velocity analysis.

[0015] In some implementations, the superposition speed of the ellipse includes a fast superposition speed and a slow superposition speed, and the orientation of the ellipse includes the orientation of the fast superposition speed.

[0016] In some implementations, the step of performing anisotropic inversion based on the current inversion array to obtain the corresponding inversion result includes:

[0017] Anisotropic inversion is performed using the preconditional conjugate gradient method on the current inversion array to obtain the corresponding inversion results.

[0018] In some implementations, during the anisotropic inversion process using the preconditional conjugate gradient method for the current inversion array, the multi-channel constraint value matrix of each surface data element in the current inversion array is calculated before each external iteration. At the same time, when inverting each surface data element, the corresponding multi-channel constraint value matrix is ​​added to the inversion process of that surface data element.

[0019] In some implementations, the calculation of the multi-channel constraint value matrix for each facet data element in the current inversion array includes:

[0020] Initialize the horizontal and vertical constraint terms to 0, and set the value of the dimensionless minimum quantity;

[0021] Based on the intermediate and inversion results of the preconditional conjugate gradient inversion obtained from the external iteration, the horizontal and vertical constraint terms are calculated.

[0022] Based on the horizontal and vertical constraint terms, calculate the sum of squares of the inversion results, the sum of the absolute values ​​of the horizontal constraint terms, and the sum of the absolute values ​​of the vertical constraint terms.

[0023] Based on the horizontal constraint terms, vertical constraint terms, the sum of squares of the inversion results, the sum of the absolute values ​​of the horizontal constraint terms, and the sum of the absolute values ​​of the vertical constraint terms, calculate the constraint values ​​corresponding to each facet data element in the current inversion array, and obtain the multi-channel constraint value matrix corresponding to each facet data element.

[0024] Secondly, embodiments of the present invention provide a pre-stack wide-azimuth anisotropy inversion apparatus, comprising:

[0025] The acquisition module is used to acquire the original pre-stack width-azimuth gather;

[0026] The correction module is used to perform azimuth AVO feature correction processing on the original pre-stack wide azimuth gather to obtain the target pre-stack wide azimuth gather;

[0027] The inversion module is used to sequentially read a set number of seismic survey lines from the target stack front-side gather data to form a current inversion array, and perform anisotropic inversion based on the current inversion array to obtain the corresponding inversion results, until all seismic survey lines in the target stack front-side gather data have obtained inversion results.

[0028] Thirdly, embodiments of the present invention provide a computer storage medium on which a computer program is stored, and when the computer program is executed by one or more processors, it implements the method described in the first aspect.

[0029] Fourthly, embodiments of the present invention provide an electronic device, including a memory and one or more processors, wherein the memory stores a computer program, and the computer program, when executed by the one or more processors, implements the method described in the first aspect.

[0030] One or more embodiments of the present invention can bring at least the following beneficial effects:

[0031] This invention performs azimuth AVO feature correction on the original pre-stack wide azimuth gathers, and then sequentially reads a predetermined number of seismic lines using a rolling extraction method to form the current inversion array. Anisotropic inversion is then performed based on this array to obtain the corresponding inversion results, until all seismic lines in the target pre-stack azimuth gather data have been inverted. Aazimuth AVO feature correction on the original data effectively improves the gather data quality, solving the problems of low quality and large data volume in pre-stack azimuth data, making it more suitable for subsequent inversion requirements. The rolling extraction of seismic data to form a specific inversion array for step-by-step inversion improves the computational efficiency and stability of the inversion. Furthermore, adding multiple constraint matrices during the inversion enhances the constraint capability of the inversion algorithm, further improving inversion accuracy and solving the problem of insufficient constraint force in the inversion algorithm. Attached Figure Description

[0032] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly described below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation on the scope.

[0033] Figure 1 This is a flowchart of a pre-stack wide-azimuth anisotropy inversion method provided by an embodiment of the present invention;

[0034] Figure 2a This refers to the number of times the front-side data of the study area is covered by the data provided in the embodiments of the present invention;

[0035] Figure 2b This is an azimuth-offset rose diagram of the front-stack position data of the study area provided in this embodiment of the invention;

[0036] Figure 3a This refers to the original azimuth gather data provided in the embodiments of the present invention;

[0037] Figure 3b This is the azimuth gather data after correction processing provided in the embodiments of the present invention;

[0038] Figure 4 This is a planar comparison diagram of crack development intensity obtained by post-stack coherence and pre-stack anisotropic inversion provided in an embodiment of the present invention;

[0039] Figure 5 This is a vector diagram of pre-stack anisotropic inversion crack development intensity and direction provided in an embodiment of the present invention;

[0040] Figure 6 This is a block diagram of a pre-stack wide-azimuth anisotropy inversion device provided in an embodiment of the present invention. Detailed Implementation

[0041] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations. Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0042] Due to tectonic stress, underground rock strata inevitably fracture, resulting in fractures. In oil and gas exploration, fractures not only serve as conduits for oil and gas migration but also as boundaries of fault-block oil and gas fields, effectively controlling their distribution. Fracture identification is fundamental to fracture research, and accurate fracture identification is crucial for oil and gas field exploration and development. The formation of fractures is controlled by multiple factors, exhibiting complex physical properties with significant lateral and longitudinal variations, demonstrating strong anisotropy. Therefore, obtaining stable anisotropic information is essential for detecting fracture characteristics.

[0043] Conventional fracture detection methods are broadly classified into post-stack and pre-stack methods. Post-stack fracture detection, based on post-stack seismic data, primarily identifies fractures by detecting discontinuities in seismic traces, such as coherence and similarity. This method has limited accuracy and struggles to identify small- to medium-scale fractures and cracks. Pre-stack fracture detection, based on pre-stack seismic data, mainly utilizes anisotropic information caused by variations in azimuth seismic amplitude and azimuth travel time for fracture prediction. Anisotropic inversion based on pre-stack azimuth data offers higher accuracy, and the fracture development intensity and direction obtained through inversion can better characterize fracture development features. However, current azimuth data is often of low quality due to various factors such as acquisition and processing, and the large data volume severely limits the application of pre-stack azimuth anisotropic inversion in fracture and crack detection. There is usually a lack of targeted correction methods for the original gathers. Directly inputting the original gathers for inversion will reduce the accuracy of the inversion. Conventional inversion algorithms also lack effective constraints. At the same time, the amount of data in the pre-stack gathers is large. Scanning the data head and recording the head information, and then locating the data volume by reading the head information file during the inversion process, requires frequent hard disk scanning, which results in slow speed. Therefore, most methods such as multi-trace stacking and reducing the order of the inversion formula are used to improve the efficiency and stability of the inversion results by sacrificing the accuracy of the inversion.

[0044] This invention addresses the problems of low quality, large data volume, and insufficient constraint of pre-stack azimuth anisotropy inversion data by proposing an optimized pre-stack wide azimuth anisotropy inversion method. It effectively improves the quality of gather data by performing azimuth AVO feature correction on the original data, making it better suited for subsequent inversion requirements. Furthermore, it enhances the constraint capability of the inversion algorithm by adding multiple constraint matrices during inversion, thereby improving inversion accuracy. Finally, it improves the computational efficiency and stability of the inversion by using a rolling extraction of seismic data to form a specific array for step-by-step inversion.

[0045] Example 1

[0046] Figure 1 A flowchart of a pre-stack wide-azimuth anisotropy inversion method is shown, as follows: Figure 1 As shown, the pre-stack wide-azimuth anisotropy inversion method provided in this embodiment includes:

[0047] Step S101: Obtain the original pre-stack wide azimuth gather.

[0048] The original pre-stack wide azimuth gathers used as input data have low data quality and large data volume. Direct inversion would reduce the accuracy of the inversion. In order to improve the efficiency and stability of the inversion results, this embodiment performs azimuth AVO (Amplitude Variation with Offset) feature correction processing on the original pre-stack wide azimuth gather data, thereby effectively improving the quality of the original pre-stack wide azimuth gather data and making it more suitable for subsequent inversion requirements.

[0049] Step S102: Perform azimuth AVO feature correction processing on the original pre-stack wide azimuth gather to obtain the target pre-stack wide azimuth gather.

[0050] In some implementations, the original pre-stack wide azimuth gather is subjected to azimuth AVO feature correction processing to obtain the target pre-stack wide azimuth gather, which may further include:

[0051] Step S102a: Calculate continuous time-varying static corrections for the original pre-stack wide azimuth gather, and use the time-varying static corrections to estimate the remaining travel time.

[0052] In some implementations, the time-varying static correction is obtained through anisotropic residual velocity analysis.

[0053] This embodiment performs anisotropic residual velocity analysis on the original pre-stack wide azimuth gather to obtain a continuous time-varying static correction on the original pre-stack wide azimuth gather. This correction is then used to estimate the residual travel time. After azimuth AVO feature correction, the phase axis of the reflected wave is flatter and the azimuth AVO feature is more obvious, which helps to improve the accuracy of anisotropic inversion of the pre-stack wide azimuth gather.

[0054] Step S102b: For each phase axis of the reflected wave, the residual travel time is fitted using the residual NMO (dynamic correction, also known as normal move out) plane of the ellipse.

[0055] This embodiment performs fitting based on the residual travel time estimated in step S102a on the elliptical residual NMO plane, and then calculates the superposition velocity and orientation of the ellipse.

[0056] Step S102c: Calculate the superposition velocity and orientation of the ellipse.

[0057] In some implementations, the superposition rate of the ellipses includes the fast superposition rate (V0). fast ) and slow superposition speed (V slow The orientation of the ellipse includes the orientation of the fast superposition velocity (β).

[0058] Step S102d: Apply the stacking velocity of the ellipse to correct the original pre-stack wide azimuth gather to obtain the target pre-stack wide azimuth gather.

[0059] By performing anisotropic residual velocity analysis on the original pre-stack wide azimuth gather, estimating the residual travel time and performing ellipse-based fitting, the stacking velocity and azimuth of the ellipse are calculated. After correcting the original pre-stack wide azimuth gather with the stacking velocity of the ellipse, a target pre-stack wide azimuth gather with a flatter phase axis of the reflected wave and more obvious AVO characteristics can be obtained, which is beneficial to improving the accuracy of subsequent anisotropic inversion.

[0060] Step S103: Sequentially read a set number of seismic survey lines from the target stack front gather data to form the current inversion array, and perform anisotropic inversion based on the current inversion array to obtain the corresponding inversion results, until all seismic survey lines in the target stack front gather data have obtained inversion results.

[0061] In some implementations, anisotropic inversion is performed based on the current inversion array to obtain the corresponding inversion results, including:

[0062] Step S103a: Perform anisotropic inversion using the preconditional conjugate gradient method on the current inversion array to obtain the corresponding inversion results.

[0063] In some implementations, during the anisotropic inversion process using the preconditional conjugate gradient method for the current inversion array, the multi-channel constraint value matrix of each surface data element in the current inversion array is calculated before each external iteration. At the same time, when inverting each surface data element, the corresponding multi-channel constraint value matrix is ​​added to the inversion process of that surface data element.

[0064] In some implementations, the multi-channel constraint value matrix of each face data element in the current inversion array is calculated, including the following process:

[0065] Step S103a-1: Initialize the horizontal and vertical constraint terms to 0, and set the value of the dimensionless minimum quantity;

[0066] Initialize the horizontal and vertical constraints:

[0067] Let dy j,l,k =0,dz j,l,k =0,

[0068] Among them, dy j,l,k This indicates the horizontal constraint term.

[0069] dz j,l,k This indicates a vertical constraint term.

[0070] j = 0, 1, 2, ..., temp_int-1

[0071] l = 0, 1, ..., xn-1

[0072] k = 0, 1, 2, ..., ns-1

[0073] j represents time.

[0074] temp_int-1 represents the maximum value of the time.

[0075] l represents the main survey line.

[0076] xn-1 represents the maximum value of the main survey line.

[0077] k represents the contact survey line,

[0078] ns-1 represents the maximum value of the connection survey line.

[0079] Initialize the value of the dimensionless minimum quantity:

[0080] Let p = 0.0001,

[0081] Where p represents the dimensionless minimum quantity;

[0082] Step S103a-2: Based on the intermediate and inversion results of the preconditional conjugate gradient inversion obtained from the external iteration, calculate the horizontal and vertical constraint terms.

[0083] Specifically, it can be calculated through the following process:

[0084] Calculate dy j,l,k =RI j,l,k -RI inv,l,k ,

[0085] Among them, RI j,l,k This is an intermediate result of preconditional conjugate gradient inversion.

[0086] RI j,l,k The initial model, initialized as input, updates its results through each external iteration of preconditional conjugate gradient inversion.

[0087] RI inv,l,k Indicates the inversion result;

[0088] calculate

[0089] Calculate dz j,l,k =RI j,l,k -RI inv,l,k-1 Where k = 1, 2, ..., ns-1;

[0090] calculate

[0091] Step S103a-3: Based on the horizontal and vertical constraint terms, calculate the sum of squares of the inversion results, the sum of the absolute values ​​of the horizontal constraint terms, and the sum of the absolute values ​​of the vertical constraint terms.

[0092] Calculate dm = sum(RI) j,l,k ·RI j,l,k ),

[0093] sum y =sum(fabs(dy) j,l,k )),

[0094] sum z =sum(fabs(dz) j,l,k )),

[0095] Here, fabs represents the absolute value operation.

[0096] `sum` represents the summation operation.

[0097] dm represents the sum of squares of the inversion results.

[0098] sum y This represents the sum of the absolute values ​​of the horizontal constraint terms.

[0099] sum z This represents the sum of the absolute values ​​of the vertical constraint terms.

[0100] Step S103a-4: Based on the horizontal constraint terms, vertical constraint terms, the sum of squares of the inversion results, the sum of the absolute values ​​of the horizontal constraint terms, and the sum of the absolute values ​​of the vertical constraint terms, calculate the constraint values ​​of each face data element in the current inversion array to obtain the multi-channel constraint value matrix corresponding to each face data element.

[0101] The multichannel constraint matrix W can be calculated using the following formula:

[0102]

[0103] Among them, w j,l,k Indicates constraint value,

[0104] α represents a set coefficient, 0 < α < 1.

[0105] Sequentially read R seismic survey line data from the target stack front gather data to form the current inversion array, and perform anisotropic inversion based on the current inversion array to obtain the corresponding inversion results. The number of lines R can be set according to the computing power of the electronic device executing this method, and this embodiment does not impose any limitations.

[0106] After performing anisotropic inversion based on the current inversion array, new R seismic line data are read from the target stack front-side gather data, and the new R seismic line data replace the R seismic line data in the current inversion array to form a new inversion array. Inversion is then performed based on this new inversion array to obtain the corresponding inversion results. This process is repeated, continuously replacing the seismic line data in the current inversion array with newly read seismic line data, and performing inversion based on the replaced inversion array, until the inversion of all seismic line data in the target stack front-side gather data is completed.

[0107] In some cases, R seismic lines are read from the target stack front gather data, starting with the seismic lines with smaller sequence numbers. For example, the first read is from line 1 to line R, and the next read starts from line R+1, reading R seismic lines. If less than R seismic lines remain during the last read, the remaining seismic lines are read to form the current inversion array for inversion.

[0108] This paper optimizes anisotropic inversion from three aspects, addressing the issues of low quality and large volume of pre-stack azimuth data, as well as insufficient constraint of the inversion algorithm. First, azimuth AVO feature correction is performed on the original input pre-stack wide azimuth data. Anisotropic residual velocity analysis makes the phase axis straighter and the AVO feature more prominent, thereby improving the quality of the gather data and thus increasing inversion accuracy. Second, multi-trace constraint matrices are added during the inversion process to optimize the inversion algorithm, enhancing its constraint capability, reducing the ambiguity of the inversion results, and further improving inversion accuracy. Third, a step-by-step inversion optimization strategy is implemented by rolling extraction of seismic data to form a specific inversion array, improving the efficiency and stability of the inversion. This results in a high-accuracy, high-efficiency, and stable inversion scheme, effectively enhancing the applicability of azimuth anisotropic inversion in fracture and crack prediction.

[0109] Example 2

[0110] This embodiment uses an oil and gas field as an example to detect fractures and cracks using the pre-stack broad azimuth anisotropy inversion method provided by this invention. Fractures and associated fractures are extremely important controlling factors for hydrocarbon accumulation and enrichment in this region. The fracture zones in the study area are long, have good vertical inheritance, and are well-developed, but exhibit large segmental differences in the plane. Conventional post-stack attribute analysis is insufficient to accurately identify fractures and associated fractures in this area.

[0111] Figures 2a-2b A rose diagram showing the number of times the overlay data is generated and the azimuth-offset distance is plotted for the study area, which is 1400 km². 2 The data storage capacity reaches 6.5T. By scanning the original front-stack azimuth data, the coverage count is greater than 250 times. At the same time, the azimuth-offset rose diagram obtained from the scan shows that the azimuth aspect ratio of the data in the study area is about 0.71, which meets the basic requirements for wide azimuth data and can be used for subsequent anisotropic inversion work.

[0112] Figures 3a-3b The original azimuth gather data and the corrected azimuth gather data are compared. The original azimuth gather data shows poor quality, low signal-to-noise ratio, and indistinct azimuth AVO characteristics. After correction processing, the gather quality is significantly improved. Not only is the signal-to-noise ratio increased, but the phase axis is also straighter, and the AVO characteristics of the target layer are more pronounced. Furthermore, the amplitude energy decreases with increasing angle, thus laying a good foundation for subsequent inversion.

[0113] Figure 4This is a planar comparison of fracture development intensity derived from conventional post-stack coherence properties and pre-stack anisotropic inversion. From the perspective of fracture morphology and distribution, the planar characteristics of the main fractures in both methods are roughly equivalent, which to some extent verifies the accuracy of the pre-stack inversion results. However, anisotropic fracture development intensity can reflect more detailed information, especially given the significant differences in planar fracture segmentation in this study area. Fracture development intensity effectively characterizes the segmentation within strike-slip fractures; the fracture development intensity is lower in the compression zone; the fracture development intensity is higher in the strike-slip segment than in the compression segment, but still not high overall; the fracture development intensity is highest in the extensional segment.

[0114] Figure 5 This is a vector diagram showing the intensity and direction of fracture development based on pre-stack anisotropic inversion. The lines in the diagram indicate the direction of fracture development, and the color intensity of the lines indicates the strength of fracture development. The base map overlays the fracture development intensity display. It can be seen from the diagram that the concentrated distribution area of ​​the lines coincides with the main fracture area indicated by the fracture development intensity. The fracture development intensity is high along the fracture zone, and the well location matches the fracture development zone.

[0115] Example 3

[0116] Figure 6 A block diagram of a pre-stack wide azimuth anisotropy inversion device is shown. The pre-stack wide azimuth anisotropy inversion device provided in this embodiment of the invention includes:

[0117] Acquisition module 301 is used to acquire the original pre-stack width-azimuth gather;

[0118] Correction module 302 is used to perform azimuth AVO feature correction processing on the original pre-stack wide azimuth gather to obtain the target pre-stack wide azimuth gather;

[0119] The inversion module 303 is used to sequentially read a set number of seismic survey lines from the target stack front gather data to form the current inversion array, and perform anisotropic inversion based on the current inversion array to obtain the corresponding inversion results, until all the seismic survey lines in the target stack front gather data have obtained inversion results.

[0120] The original pre-stack wide azimuth gathers used as input data have low data quality and large data volume. Direct inversion would reduce the accuracy of the inversion. In order to improve the efficiency and stability of the inversion results, this embodiment performs azimuth AVO (Amplitude Variation with Offset) feature correction processing on the original pre-stack wide azimuth gather data, thereby effectively improving the quality of the original pre-stack wide azimuth gather data and making it more suitable for subsequent inversion requirements.

[0121] In some implementations, the original pre-stack wide azimuth gather is subjected to azimuth AVO feature correction processing to obtain the target pre-stack wide azimuth gather, which may further include:

[0122] For the original pre-stack wide azimuth gather, a continuous time-varying static correction is calculated, and the remaining travel time is estimated using the time-varying static correction. For each reflected wave phase axis, the remaining travel time is fitted using the residual NMO (dynamic correction, also known as normal move out) plane of the ellipse. The stacking velocity and azimuth of the ellipse are calculated. The stacking velocity of the ellipse is applied to correct the original pre-stack wide azimuth gather to obtain the target pre-stack wide azimuth gather.

[0123] In some implementations, the time-varying static correction is obtained through anisotropic residual velocity analysis.

[0124] This embodiment performs anisotropic residual velocity analysis on the original pre-stack wide azimuth gather to obtain a continuous time-varying static correction on the original pre-stack wide azimuth gather. This correction is then used to estimate the residual travel time. After azimuth AVO feature correction, the phase axis of the reflected wave is flatter and the azimuth AVO feature is more obvious, which helps to improve the accuracy of anisotropic inversion of the pre-stack wide azimuth gather.

[0125] This embodiment performs fitting based on the estimated remaining travel time using the elliptical residual NMO plane, and then calculates the superposition velocity and orientation of the ellipse.

[0126] In some implementations, the superposition rate of the ellipses includes the fast superposition rate (V0). fast ) and slow superposition speed (V slow The orientation of the ellipse includes the orientation of the fast superposition velocity (β).

[0127] By performing anisotropic residual velocity analysis on the original pre-stack wide azimuth gather, estimating the residual travel time and performing ellipse-based fitting, the stacking velocity and azimuth of the ellipse are calculated. After correcting the original pre-stack wide azimuth gather with the stacking velocity of the ellipse, a target pre-stack wide azimuth gather with a flatter phase axis of the reflected wave and more obvious AVO characteristics can be obtained, which is beneficial to improving the accuracy of subsequent anisotropic inversion.

[0128] In some implementations, anisotropic inversion is performed based on the current inversion array to obtain the corresponding inversion results, including: performing anisotropic inversion using the preconditional conjugate gradient method on the current inversion array to obtain the corresponding inversion results.

[0129] In some implementations, during the anisotropic inversion process using the preconditional conjugate gradient method for the current inversion array, the multi-channel constraint value matrix of each surface data element in the current inversion array is calculated before each external iteration. At the same time, when inverting each surface data element, the corresponding multi-channel constraint value matrix is ​​added to the inversion process of that surface data element.

[0130] In some implementations, the multi-channel constraint value matrix of each face data element in the current inversion array is calculated, including the following process:

[0131] First, initialize the horizontal and vertical constraint terms to 0, and set the value of the dimensionless minimum quantity.

[0132] Initialize the horizontal and vertical constraints:

[0133] Let dy j,l,k =0,dz j,l,k =0,

[0134] Among them, dy j,l,k This indicates the horizontal constraint term.

[0135] dz j,l,k This indicates a vertical constraint term.

[0136] j = 0, 1, 2, ..., temp_int-1

[0137] l = 0, 1, ..., xn-1

[0138] k = 0, 1, 2, ..., ns-1

[0139] j represents time.

[0140] temp_int-1 represents the maximum value of the time.

[0141] l represents the main survey line.

[0142] xn-1 represents the maximum value of the main survey line.

[0143] k represents the contact survey line,

[0144] ns-1 represents the maximum value of the connection survey line.

[0145] Initialize the value of the dimensionless minimum quantity:

[0146] Let p = 0.0001,

[0147] Where p represents the dimensionless minimum quantity;

[0148] Secondly, based on the intermediate and inversion results of the preconditioned conjugate gradient inversion obtained from the external iteration, the horizontal and vertical constraint terms are calculated.

[0149] Specifically, it can be calculated through the following process:

[0150] Calculate dy j,l,k =RI j,l,k -RI inv,l,k ,

[0151] Among them, RI j,l,k This is an intermediate result of preconditional conjugate gradient inversion.

[0152] RI j,l,k The initial model, initialized as input, updates its results through each external iteration of preconditional conjugate gradient inversion.

[0153] RI inv,l,k Indicates the inversion result;

[0154] calculate

[0155] Calculate dz j,l,k =RI j,l,k -RI inv,l,k-1 Where k = 1, 2, ..., ns-1;

[0156] calculate

[0157] Next, based on the horizontal and vertical constraint terms, calculate the sum of squares of the inversion results, the sum of the absolute values ​​of the horizontal constraint terms, and the sum of the absolute values ​​of the vertical constraint terms.

[0158] Calculate dm = sum(RI) j,l,k ·RI j,l,k ),

[0159] sum y =sum(fabs(dy) j,l,k )),

[0160] sum z =sum(fabs(dz) j,l,k )),

[0161] Here, fabs represents the absolute value operation.

[0162] `sum` represents the summation operation.

[0163] dm represents the sum of squares of the inversion results.

[0164] sum y This represents the sum of the absolute values ​​of the horizontal constraint terms.

[0165] sum z This represents the sum of the absolute values ​​of the vertical constraint terms.

[0166] Finally, based on the horizontal constraint terms, vertical constraint terms, the sum of squares of the inversion results, the sum of the absolute values ​​of the horizontal constraint terms, and the sum of the absolute values ​​of the vertical constraint terms, the constraint values ​​corresponding to each facet data element in the current inversion array are calculated to obtain the multi-channel constraint value matrix corresponding to each facet data element.

[0167] The multichannel constraint matrix W can be calculated using the following formula:

[0168]

[0169] Among them, w j,l,k Indicates constraint value,

[0170] α represents a set coefficient, 0 < α < 1.

[0171] Sequentially read R seismic survey line data from the target stack front gather data to form the current inversion array, and perform anisotropic inversion based on the current inversion array to obtain the corresponding inversion results. The number of lines R can be set according to the computing power of the electronic device executing this method, and this embodiment does not impose any limitations.

[0172] After performing anisotropic inversion based on the current inversion array, new R seismic line data are read from the target stack front-side gather data, and the new R seismic line data replace the R seismic line data in the current inversion array to form a new inversion array. Inversion is then performed based on this new inversion array to obtain the corresponding inversion results. This process is repeated, continuously replacing the seismic line data in the current inversion array with newly read seismic line data, and performing inversion based on the replaced inversion array, until the inversion of all seismic line data in the target stack front-side gather data is completed.

[0173] In some cases, R seismic lines are read from the target stack front gather data, starting with the seismic lines with smaller sequence numbers. For example, the first read is from line 1 to line R, and the next read starts from line R+1, reading R seismic lines. If less than R seismic lines remain during the last read, the remaining seismic lines are read to form the current inversion array for inversion.

[0174] This embodiment optimizes the input data by performing azimuth AVO feature correction processing on the original pre-stack wide azimuth gather. By adding multi-channel constraint matrices during the anisotropic inversion process and extracting data in a rolling manner to form an inversion array for step-by-step inversion optimization, the inversion strategy is effectively improved, enhancing the accuracy, efficiency, and stability of azimuth anisotropy inversion and increasing the applicability of pre-stack wide azimuth anisotropy inversion in fracture and crack detection.

[0175] Those skilled in the art will understand that the above-described modules or steps can be implemented using general-purpose computing devices. They can be centralized on a single computing device or distributed across a network of multiple computing devices. Optionally, they can be implemented using computer-executable program code, thereby storing them in a storage device for execution by the computing device, or fabricating them separately as individual integrated circuit modules, or fabricating multiple modules or steps into a single integrated circuit module. This invention is not limited to any specific hardware and software combination.

[0176] Example 4

[0177] This embodiment provides a computer storage medium on which a computer program is stored. When the computer program is executed by one or more processors, it implements the method of Embodiment 1.

[0178] The computer-readable storage medium can be implemented by any type of volatile or non-volatile storage device or a combination thereof, such as Static Random Access Memory (SRAM), Electrically Erasable Programmable Read-Only Memory (EEPROM), Erasable Programmable Read-Only Memory (EPROM), Programmable Read-Only Memory (PROM), Read-Only Memory (ROM), magnetic storage, flash memory, magnetic disk, or optical disk.

[0179] The method implemented by the computer-readable storage medium in this embodiment includes:

[0180] Step S101: Obtain the original pre-stack wide azimuth gather.

[0181] The original pre-stack wide azimuth gathers used as input data have low data quality and large data volume. Direct inversion would reduce the accuracy of the inversion. In order to improve the efficiency and stability of the inversion results, this embodiment performs azimuth AVO (Amplitude Variation with Offset) feature correction processing on the original pre-stack wide azimuth gather data, thereby effectively improving the quality of the original pre-stack wide azimuth gather data and making it more suitable for subsequent inversion requirements.

[0182] Step S102: Perform azimuth AVO feature correction processing on the original pre-stack wide azimuth gather to obtain the target pre-stack wide azimuth gather.

[0183] In some implementations, the original pre-stack wide azimuth gather is subjected to azimuth AVO feature correction processing to obtain the target pre-stack wide azimuth gather, which may further include:

[0184] Step S102a: Calculate continuous time-varying static corrections for the original pre-stack wide azimuth gather, and use the time-varying static corrections to estimate the remaining travel time.

[0185] In some implementations, the time-varying static correction is obtained through anisotropic residual velocity analysis.

[0186] This embodiment performs anisotropic residual velocity analysis on the original pre-stack wide azimuth gather to obtain a continuous time-varying static correction on the original pre-stack wide azimuth gather. This correction is then used to estimate the residual travel time. After azimuth AVO feature correction, the phase axis of the reflected wave is flatter and the azimuth AVO feature is more obvious, which helps to improve the accuracy of anisotropic inversion of the pre-stack wide azimuth gather.

[0187] Step S102b: For each phase axis of the reflected wave, the residual travel time is fitted using the residual NMO (dynamic correction, also known as normal move out) plane of the ellipse.

[0188] Step S102c: Calculate the superposition velocity and orientation of the ellipse.

[0189] In some implementations, the superposition rate of the ellipses includes the fast superposition rate (V0). fast ) and slow superposition speed (V slow The orientation of the ellipse includes the orientation of the fast superposition velocity (β).

[0190] Step S102d: Apply the stacking velocity of the ellipse to correct the original pre-stack wide azimuth gather to obtain the target pre-stack wide azimuth gather.

[0191] By performing anisotropic residual velocity analysis on the original pre-stack wide azimuth gather, estimating the residual travel time and performing ellipse-based fitting, the stacking velocity and azimuth of the ellipse are calculated. After correcting the original pre-stack wide azimuth gather with the stacking velocity of the ellipse, a target pre-stack wide azimuth gather with a flatter phase axis of the reflected wave and more obvious AVO characteristics can be obtained, which is beneficial to improving the accuracy of subsequent anisotropic inversion.

[0192] Step S103: Sequentially read a set number of seismic survey lines from the target stack front gather data to form the current inversion array, and perform anisotropic inversion based on the current inversion array to obtain the corresponding inversion results, until all seismic survey lines in the target stack front gather data have obtained inversion results.

[0193] In some implementations, anisotropic inversion is performed based on the current inversion array to obtain the corresponding inversion results, including:

[0194] Step S103a: Perform anisotropic inversion using the preconditional conjugate gradient method on the current inversion array to obtain the corresponding inversion results.

[0195] In some implementations, during the anisotropic inversion process using the preconditional conjugate gradient method for the current inversion array, the multi-channel constraint value matrix of each surface data element in the current inversion array is calculated before each external iteration. At the same time, when inverting each surface data element, the corresponding multi-channel constraint value matrix is ​​added to the inversion process of that surface data element.

[0196] In some implementations, the multi-channel constraint value matrix of each face data element in the current inversion array is calculated, including the following process:

[0197] Step S103a-1: Initialize the horizontal and vertical constraint terms to 0, and set the value of the dimensionless minimum quantity;

[0198] Initialize the horizontal and vertical constraints:

[0199] Let dy j,l,k =0,dz j,l,k =0,

[0200] Among them, dy j,l,k This indicates the horizontal constraint term.

[0201] dz j,l,k This indicates a vertical constraint term.

[0202] j = 0, 1, 2, ..., temp_int-1

[0203] l = 0, 1, ..., xn-1

[0204] k = 0, 1, 2, ..., ns-1

[0205] j represents time.

[0206] temp_int-1 represents the maximum value of the time.

[0207] l represents the main survey line.

[0208] xn-1 represents the maximum value of the main survey line.

[0209] k represents the contact survey line,

[0210] ns-1 represents the maximum value of the connection survey line.

[0211] Initialize the value of the dimensionless minimum quantity:

[0212] Let p = 0.0001,

[0213] Where p represents the dimensionless minimum quantity;

[0214] Step S103a-2: Based on the intermediate and inversion results of the preconditional conjugate gradient inversion obtained from the external iteration, calculate the horizontal and vertical constraint terms.

[0215] Specifically, it can be calculated through the following process:

[0216] Calculate dy j,l,k =RI j,l,k -RI inv,l,k ,

[0217] Among them, RI j,l,k This is an intermediate result of preconditional conjugate gradient inversion.

[0218] RI j,l,k The initial model, initialized as input, updates its results through each external iteration of preconditional conjugate gradient inversion.

[0219] RI inv,l,k Indicates the inversion result;

[0220] calculate

[0221] Calculate dz j,l,k =RI j,l,k -RI inv,l,k-1 Where k = 1, 2, ..., ns-1;

[0222] calculate

[0223] Step S103a-3: Based on the horizontal and vertical constraint terms, calculate the sum of squares of the inversion results, the sum of the absolute values ​​of the horizontal constraint terms, and the sum of the absolute values ​​of the vertical constraint terms.

[0224] Calculate dm = sum(RI) j,l,k ·RI j,l,k ),

[0225] sum y =sum(fabs(dy) j,l,k )),

[0226] sum z =sum(fabs(dz) j,l,k )),

[0227] Here, fabs represents the absolute value operation.

[0228] `sum` represents the summation operation.

[0229] dm represents the sum of squares of the inversion results.

[0230] sum y This represents the sum of the absolute values ​​of the horizontal constraint terms.

[0231] sum z This represents the sum of the absolute values ​​of the vertical constraint terms.

[0232] Step S103a-4: Based on the horizontal constraint terms, vertical constraint terms, the sum of squares of the inversion results, the sum of the absolute values ​​of the horizontal constraint terms, and the sum of the absolute values ​​of the vertical constraint terms, calculate the constraint values ​​of each face data element in the current inversion array to obtain the multi-channel constraint value matrix corresponding to each face data element.

[0233] The multichannel constraint matrix W can be calculated using the following formula:

[0234]

[0235] Among them, w j,l,k Indicates constraint value,

[0236] α represents a set coefficient, 0 < α < 1.

[0237] Sequentially read R seismic survey line data from the target stack front gather data to form the current inversion array, and perform anisotropic inversion based on the current inversion array to obtain the corresponding inversion results. The number of lines R can be set according to the computing power of the electronic device executing this method, and this embodiment does not impose any limitations.

[0238] After performing anisotropic inversion based on the current inversion array, new R seismic line data are read from the target stack front-side gather data, and the new R seismic line data replace the R seismic line data in the current inversion array to form a new inversion array. Inversion is then performed based on this new inversion array to obtain the corresponding inversion results. This process is repeated, continuously replacing the seismic line data in the current inversion array with newly read seismic line data, and performing inversion based on the replaced inversion array, until the inversion of all seismic line data in the target stack front-side gather data is completed.

[0239] In some cases, R seismic lines are read from the target stack front gather data, starting with the seismic lines with smaller sequence numbers. For example, the first read is from line 1 to line R, and the next read starts from line R+1, reading R seismic lines. If less than R seismic lines remain during the last read, the remaining seismic lines are read to form the current inversion array for inversion.

[0240] Example 5

[0241] This embodiment provides an electronic device, including a memory and one or more processors. The memory stores a computer program, which, when executed by one or more processors, implements the method of Embodiment 1.

[0242] In practical applications, the processor can be implemented as an Application Specific Integrated Circuit (ASIC), Digital Signal Processor (DSP), Digital Signal Processing Device (DSPD), Programmable Logic Device (PLD), Field Programmable Gate Array (FPGA), controller, microcontroller unit (MCU), microprocessor, or other electronic components to execute the methods described in the above embodiments.

[0243] The method implemented by the computer program of this embodiment when executed by one or more processors includes:

[0244] Step S101: Obtain the original pre-stack wide azimuth gather.

[0245] The original pre-stack wide azimuth gathers used as input data have low data quality and large data volume. Direct inversion would reduce the accuracy of the inversion. In order to improve the efficiency and stability of the inversion results, this embodiment performs azimuth AVO (Amplitude Variation with Offset) feature correction processing on the original pre-stack wide azimuth gather data, thereby effectively improving the quality of the original pre-stack wide azimuth gather data and making it more suitable for subsequent inversion requirements.

[0246] Step S102: Perform azimuth AVO feature correction processing on the original pre-stack wide azimuth gather to obtain the target pre-stack wide azimuth gather.

[0247] In some implementations, the original pre-stack wide azimuth gather is subjected to azimuth AVO feature correction processing to obtain the target pre-stack wide azimuth gather, which may further include:

[0248] Step S102a: Calculate continuous time-varying static corrections for the original pre-stack wide azimuth gather, and use the time-varying static corrections to estimate the remaining travel time.

[0249] In some implementations, the time-varying static correction is obtained through anisotropic residual velocity analysis.

[0250] This embodiment performs anisotropic residual velocity analysis on the original pre-stack wide azimuth gather to obtain a continuous time-varying static correction on the original pre-stack wide azimuth gather. This correction is then used to estimate the residual travel time. After azimuth AVO feature correction, the phase axis of the reflected wave is flatter and the azimuth AVO feature is more obvious, which helps to improve the accuracy of anisotropic inversion of the pre-stack wide azimuth gather.

[0251] Step S102b: For each phase axis of the reflected wave, the residual travel time is fitted using the residual NMO (dynamic correction, also known as normal move out) plane of the ellipse.

[0252] Step S102c: Calculate the superposition velocity and orientation of the ellipse.

[0253] In some implementations, the superposition rate of the ellipses includes the fast superposition rate (V0). fast ) and slow superposition speed (V slow The orientation of the ellipse includes the orientation of the fast superposition velocity (β).

[0254] Step S102d: Apply the stacking velocity of the ellipse to correct the original pre-stack wide azimuth gather to obtain the target pre-stack wide azimuth gather.

[0255] By performing anisotropic residual velocity analysis on the original pre-stack wide azimuth gather, estimating the residual travel time and performing ellipse-based fitting, the stacking velocity and azimuth of the ellipse are calculated. After correcting the original pre-stack wide azimuth gather with the stacking velocity of the ellipse, a target pre-stack wide azimuth gather with a flatter phase axis of the reflected wave and more obvious AVO characteristics can be obtained, which is beneficial to improving the accuracy of subsequent anisotropic inversion.

[0256] Step S103: Sequentially read a set number of seismic survey lines from the target stack front gather data to form the current inversion array, and perform anisotropic inversion based on the current inversion array to obtain the corresponding inversion results, until all seismic survey lines in the target stack front gather data have obtained inversion results.

[0257] In some implementations, anisotropic inversion is performed based on the current inversion array to obtain the corresponding inversion results, including:

[0258] Step S103a: Perform anisotropic inversion using the preconditional conjugate gradient method on the current inversion array to obtain the corresponding inversion results.

[0259] In some implementations, during the anisotropic inversion process using the preconditional conjugate gradient method for the current inversion array, the multi-channel constraint value matrix of each surface data element in the current inversion array is calculated before each external iteration. At the same time, when inverting each surface data element, the corresponding multi-channel constraint value matrix is ​​added to the inversion process of that surface data element.

[0260] In some implementations, the multi-channel constraint value matrix of each face data element in the current inversion array is calculated, including the following process:

[0261] Step S103a-1: Initialize the horizontal and vertical constraint terms to 0, and set the value of the dimensionless minimum quantity;

[0262] Initialize the horizontal and vertical constraints:

[0263] Let dy j,l,k =0,dz j,l,k =0,

[0264] Among them, dy j,l,k This indicates the horizontal constraint term.

[0265] dz j,l,k This indicates a vertical constraint term.

[0266] j = 0, 1, 2, ..., temp_int-1

[0267] l = 0, 1, ..., xn-1

[0268] k = 0, 1, 2, ..., ns-1

[0269] j represents time.

[0270] temp_int-1 represents the maximum value of the time.

[0271] l represents the main survey line.

[0272] xn-1 represents the maximum value of the main survey line.

[0273] k represents the contact survey line,

[0274] ns-1 represents the maximum value of the connection survey line.

[0275] Initialize the value of the dimensionless minimum quantity:

[0276] Let p = 0.0001,

[0277] Where p represents the dimensionless minimum quantity;

[0278] Step S103a-2: Based on the intermediate and inversion results of the preconditional conjugate gradient inversion obtained from the external iteration, calculate the horizontal and vertical constraint terms.

[0279] Specifically, it can be calculated through the following process:

[0280] Calculate dy j,l,k =RI j,l,k -RI inv,l,k ,

[0281] Among them, RI j,l,k This is an intermediate result of preconditional conjugate gradient inversion.

[0282] RI j,l,k The initial model, initialized as input, updates its results through each external iteration of preconditional conjugate gradient inversion.

[0283] RI inv,l,k Indicates the inversion result;

[0284] calculate

[0285] Calculate dz j,l,k =RI j,l,k -RI inv,l,k-1 Where k = 1, 2, ..., ns-1;

[0286] calculate

[0287] Step S103a-3: Based on the horizontal and vertical constraint terms, calculate the sum of squares of the inversion results, the sum of the absolute values ​​of the horizontal constraint terms, and the sum of the absolute values ​​of the vertical constraint terms.

[0288] Calculate dm = sum(RI) j,l,k ·RI j,l,k ),

[0289] sum y =sum(fabs(dy) j,l,k )),

[0290] sum z =sum(fabs(dz) j,l,k )),

[0291] Here, fabs represents the absolute value operation.

[0292] `sum` represents the summation operation.

[0293] dm represents the sum of squares of the inversion results.

[0294] sum y This represents the sum of the absolute values ​​of the horizontal constraint terms.

[0295] sum z This represents the sum of the absolute values ​​of the vertical constraint terms.

[0296] Step S103a-4: Based on the horizontal constraint terms, vertical constraint terms, the sum of squares of the inversion results, the sum of the absolute values ​​of the horizontal constraint terms, and the sum of the absolute values ​​of the vertical constraint terms, calculate the constraint values ​​of each face data element in the current inversion array to obtain the multi-channel constraint value matrix corresponding to each face data element.

[0297] The multichannel constraint matrix W can be calculated using the following formula:

[0298]

[0299] Among them, w j,l,k Indicates constraint value,

[0300] α represents a set coefficient, 0 < α < 1.

[0301] Sequentially read R seismic survey line data from the target stack front gather data to form the current inversion array, and perform anisotropic inversion based on the current inversion array to obtain the corresponding inversion results. The number of lines R can be set according to the computing power of the electronic device executing this method, and this embodiment does not impose any limitations.

[0302] After performing anisotropic inversion based on the current inversion array, new R seismic line data are read from the target stack front-side gather data, and the new R seismic line data replace the R seismic line data in the current inversion array to form a new inversion array. Inversion is then performed based on this new inversion array to obtain the corresponding inversion results. This process is repeated, continuously replacing the seismic line data in the current inversion array with newly read seismic line data, and performing inversion based on the replaced inversion array, until the inversion of all seismic line data in the target stack front-side gather data is completed.

[0303] In some cases, R seismic lines are read from the target stack front gather data, starting with the seismic lines with smaller sequence numbers. For example, the first read is from line 1 to line R, and the next read starts from line R+1, reading R seismic lines. If less than R seismic lines remain during the last read, the remaining seismic lines are read to form the current inversion array for inversion.

[0304] This invention addresses the issue from three main aspects: first, by optimizing the input data and correcting the original azimuth data through azimuth AVO feature correction, effectively improving the quality of the gather data and making it more suitable for subsequent inversion; second, by optimizing the inversion algorithm by adding multiple constraint matrices to enhance the algorithm's constraint capabilities and further improve inversion accuracy; and third, by optimizing the inversion strategy by using a rolling extraction of seismic data to form a specific array for step-by-step inversion, thus solving the computational efficiency and stability issues caused by the large amount of data.

[0305] In the several embodiments provided in this invention, it should be understood that the disclosed systems and methods can also be implemented in other ways. The system and method embodiments described above are merely illustrative.

[0306] It should be noted that, in this document, the terms "first," "second," etc., used in the specification, claims, and accompanying drawings of this application are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. The terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element.

[0307] While the embodiments disclosed in this invention are as described above, the content is merely for the purpose of facilitating understanding of the invention and is not intended to limit the invention. Any person skilled in the art to which this invention pertains may make any modifications and variations in form and detail of the implementation without departing from the spirit and scope disclosed herein; however, the scope of patent protection for this invention shall still be determined by the scope defined in the appended claims.

Claims

1. A pre-stack wide-azimuth anisotropy inversion method, characterized in that, include: Obtain the original pre-stack wide azimuth gather; The original pre-stack wide azimuth gather is subjected to azimuth AVO feature correction processing to obtain the target pre-stack wide azimuth gather; the step of performing azimuth AVO feature correction processing on the original pre-stack wide azimuth gather to obtain the target pre-stack wide azimuth gather includes: calculating a continuous time-varying static correction for the original pre-stack wide azimuth gather, and using the time-varying static correction to estimate the remaining travel time. For each reflected wave phase axis, the residual travel time is fitted using the residual NMO plane of the ellipse; the stacking velocity and azimuth of the ellipse are calculated; the stacking velocity of the ellipse is applied to correct the original pre-stack wide azimuth gather to obtain the target pre-stack wide azimuth gather. A set number of seismic survey lines are sequentially read from the target stack front-side gather data to form a current inversion array. Anisotropic inversion is then performed based on the current inversion array to obtain the corresponding inversion results, until all seismic survey lines in the target stack front-side gather data have obtained inversion results.

2. The pre-stack wide-azimuth anisotropy inversion method according to claim 1, characterized in that, The time-varying static correction is obtained through anisotropic residual velocity analysis.

3. The pre-stack wide-azimuth anisotropy inversion method according to claim 1, characterized in that, The superposition speed of the ellipse includes a fast superposition speed and a slow superposition speed, and the orientation of the ellipse includes the orientation of the fast superposition speed.

4. The pre-stack wide-azimuth anisotropy inversion method according to claim 1, characterized in that, The process of obtaining the corresponding inversion result by performing anisotropic inversion based on the current inversion array includes: Anisotropic inversion is performed using the preconditional conjugate gradient method on the current inversion array to obtain the corresponding inversion results.

5. The pre-stack wide-azimuth anisotropy inversion method according to claim 4, characterized in that, In the process of anisotropic inversion using the preconditional conjugate gradient method for the current inversion array, the multi-channel constraint value matrix of each surface data element in the current inversion array is calculated before each external iteration. At the same time, when inverting each surface data element, the corresponding multi-channel constraint value matrix is ​​added to the inversion process of that surface data element.

6. The pre-stack wide-azimuth anisotropy inversion method according to claim 5, characterized in that, The calculation of the multichannel constraint value matrix for each facet data element in the current inversion array includes: Initialize the horizontal and vertical constraint terms to 0, and set the value of the dimensionless minimum quantity; Based on the intermediate and inversion results of the preconditional conjugate gradient inversion obtained from the external iteration, the horizontal and vertical constraint terms are calculated. Based on the horizontal and vertical constraint terms, calculate the sum of squares of the inversion results, the sum of the absolute values ​​of the horizontal constraint terms, and the sum of the absolute values ​​of the vertical constraint terms. Based on the horizontal constraint terms, vertical constraint terms, the sum of squares of the inversion results, the sum of the absolute values ​​of the horizontal constraint terms, and the sum of the absolute values ​​of the vertical constraint terms, calculate the constraint values ​​corresponding to each facet data element in the current inversion array, and obtain the multi-channel constraint value matrix corresponding to each facet data element.

7. A pre-stack wide-azimuth anisotropy inversion device, characterized in that, include: The acquisition module is used to acquire the original pre-stack width-azimuth gather; A correction module is used to perform azimuth AVO feature correction processing on the original pre-stack wide azimuth gather to obtain a target pre-stack wide azimuth gather. The azimuth AVO feature correction processing on the original pre-stack wide azimuth gather to obtain the target pre-stack wide azimuth gather includes: calculating a continuous time-varying static correction for the original pre-stack wide azimuth gather, and using the time-varying static correction to estimate the remaining travel time; for each reflected wave phase axis, fitting the remaining travel time using the residual NMO plane of the ellipse; calculating the stacking velocity and azimuth of the ellipse; and applying the stacking velocity of the ellipse to correct the original pre-stack wide azimuth gather to obtain the target pre-stack wide azimuth gather. The inversion module is used to sequentially read a set number of seismic survey lines from the target stack front-side gather data to form a current inversion array, and perform anisotropic inversion based on the current inversion array to obtain the corresponding inversion results, until all seismic survey lines in the target stack front-side gather data have obtained inversion results.

8. A computer storage medium, characterized in that, The computer storage medium stores a computer program, which, when executed by one or more processors, implements the method as described in any one of claims 1 to 6.

9. An electronic device, characterized in that, It includes a memory and one or more processors, wherein the memory stores a computer program that, when executed by the one or more processors, implements the method as described in any one of claims 1 to 6.

Images