Seismic facies guided high-precision geological anomalous body identification method

Abstract

The invention discloses a seismic facies-guided high-precision geological anomalous body identification method, belongs to the technical field of seismic exploration data processing, and aims to consider the change of seismic facies in different regions, perform cross-correlation time shifting on layer-leveled seismic data channel by channel, improve the layer-leveled precision and obtain high-precision geological anomalous body identification factors of a target horizon. And a good foundation is laid for post-fault development law and reservoir identification and prediction. The method is simple in algorithm, high in flexibility, wide in applicability, good in result stability and high in precision.

Description

Technical Field

[0001] The present invention belongs to the technical field of seismic exploration data processing, and particularly relates to a method for identifying high-precision geological anomalies guided by seismic facies. Background Technique

[0002] Geological anomalies refer to geological bodies that cause geophysical anomalies, which are manifested as discontinuous and energy-abnormal reflections in seismic data, showing heterogeneous characteristics, and are therefore also called heterogeneous bodies. During multiple geological movements from ancient times to the present, a large number of geological anomalies such as salt domes, karst caves, faults, river channels, fractures, and volcanoes have been formed underground. These geological anomalies are usually both channels for oil and gas migration and main sites for oil and gas accumulation. Identifying geological anomalies in seismic data and determining the morphology and spatial distribution characteristics of geological anomalies are crucial for finding oil and gas traps and reservoir prediction. Therefore, it is necessary to strip continuous strong reflections to highlight the response characteristics of the target reservoir, laying a good foundation for improving the accuracy of reservoir prediction.

[0003] Conventional methods for identifying geological anomalies, such as coherence body and similarity techniques, are easily affected by complex seismic data and the accuracy of horizon tracking, are insensitive to large-dip strata, and it is difficult to guarantee the resolution. Summary of the Invention

[0004] In order to solve the technical problems existing in the background technique, the present invention aims to provide a method for identifying high-precision geological anomalies guided by seismic facies. Considering the changes in seismic facies in different regions, cross-correlation time shift is performed on the flattened seismic data channel by channel to improve the accuracy of flattening and obtain high-precision geological anomaly identification factors for the target horizon.

[0005] In order to solve the technical problems, the technical solution of the present invention is:

[0006] A method for identifying high-precision geological anomalies guided by seismic facies, the method comprising:

[0007] S1: Obtain three-dimensional seismic data, target horizon data, and sampling interval, perform flattening on the data to align the strata, and obtain the flattened seismic data;

[0008] S2: Take each channel data of the flattened seismic data as the target channel, extract the three-dimensional data subset around the target channel, and calculate the Euclidean distance between the center point of each channel data in the data subset and the center point of the target channel to obtain the three-dimensional data subset and Euclidean distance matrix corresponding to the target channel;

[0009] S3: Calculate the cross-correlation between each channel data in the data subset and the target channel, determine the time shift amount corresponding to each channel data, indicating the time difference between the horizon time of each channel data in the data subset and the time of the target channel, and obtain the time shift amount matrix of each channel data in the data subset relative to the target channel;

[0010] S4: Based on the time shift amount matrix, determine whether the absolute value of the time shift amount of all elements in the matrix is less than 1. If it is satisfied, stop the iteration. Otherwise, update the time shift result of the data subset, recalculate the cross-correlation, and continue the iteration until the time shift amount matrix meets the iteration termination condition or reaches the maximum number of iterations, and output the optimized final time shift amount matrix and the time shift result of the data subset;

[0011] S5: Calculate the correlation coefficient matrix of the target trace based on the time shift result of the final data subset, weight the correlation coefficient matrix according to the Euclidean distance matrix, sort all elements of the weighted correlation coefficient matrix from largest to smallest, and set a threshold for screening to remove elements with low correlation coefficients, obtaining the screened data elements and their quantities;

[0012] S6: Use the output result of step S5 to calculate the standard slope, fitting slope, reciprocal of the mean, and reciprocal of the median to obtain each recognition factor;

[0013] S7: Combine different recognition factors to calculate multiple combined recognition factors, obtaining multiple combined recognition factors;

[0014] S8: Traverse all data of the flattened seismic data, calculate the recognition factors of the target horizon in the entire area, and obtain the distribution map of the geological anomaly recognition factors of the target horizon in the whole area.

[0015] Further, in step S1, the flattened data Q is expressed as:

[0016] Q[n - T[m][l]][m][l] = P[n][m][l];

[0017] where n ∈ [0, N - 1], m ∈ [1, M], l ∈ [1, L].

[0018] Further, in step S2, the three-dimensional data subset A m,l and the Euclidean distance O m,l are:

[0019] A m,l [n 2 [x][y] = Q[n 1 + n 2 [m - I - 1 + x][l - J - 1 + y];

[0020]

[0021] where n 2 ∈[-N 2 / 2, N 2 / 2], x ∈ [1, 2I + 1], y ∈ [1, 2J + 1], where Δx and Δy are the spatial sampling intervals in the inline and crossline directions respectively.

[0022] Furthermore, in the step S3, the time shift amount data is calculated by the following formula:

[0023]

[0024] where n 3 ∈ [-N 2 + 1, N 2 - 1], and argmax() represents the function for calculating the position of the maximum value.

[0025] Furthermore, in the step S4, the sum W m,l of the first k time shift amounts is calculated by the following formula:

[0026] W m,l [x][y] = ∑D m,l,k [x][y];

[0027] According to W m,l the three-dimensional data subset A m,l after time shift is updated to:

[0028] A m,l [n 2 [x][y] = Q[n 1 + n 2 + W m,l [x][y]][m - I - 1 + x][l - J - 1 + y];

[0029] The final correlation coefficient C m,l is:

[0030]

[0031] Furthermore, in the step S5, the weighted correlation coefficient G m,l is calculated by the following formula:

[0032]

[0033] where α ≥ 0 is the weighting parameter.

[0034] Furthermore, in the step S6, the standard slope S of the geological anomaly identification factor is calculated by the following formula:

[0035]

[0036] The fitting slope F is calculated by the following formula:

[0037]

[0038] where k 1 ∈ [1, K m,l , and the mean reciprocal E is calculated by the following formula:

[0039]

[0040] The median reciprocal U is calculated by the following formula:

[0041]

[0042] The median reciprocal U, the mean reciprocal E, and the combined recognition factor are calculated by the following formula:

[0043]

[0044] This technology combines seismic facies analysis and mathematical optimization methods to achieve high-precision identification of geological anomalies, and has the following advantages:

[0045] Accurate horizon flattening processing to reduce the influence of formation undulations: Through the horizon flattening technology, seismic event axes are aligned, improving data comparability and enhancing the ability to identify weak anomalies.

[0046] Adaptive three-dimensional data subset extraction to enhance local feature analysis: An adaptive subset extraction method based on seismic facies is adopted to effectively retain the spatial information of geological structures and improve the identification accuracy.

[0047] Iteratively optimizing the time shift amount to improve the correlation matching degree: Combining cross-correlation calculation and iterative adjustment methods to ensure the minimum time shift error between data and improve the matching degree of seismic signals.

[0048] Euclidean distance weighting to improve spatial consistency: Introducing Euclidean distance weights when calculating the correlation coefficient, making neighboring data points contribute more, reducing the influence of noise, and improving the stability of identification.

[0049] Multi-level recognition factors to enhance the ability to characterize anomalies: Calculating multiple factors such as standard slope, fitting slope, mean reciprocal, and median reciprocal to comprehensively characterize the characteristics of geological anomalies and improve the accuracy and reliability of identification.

[0050] Combined recognition factor to improve the resolution of anomalies: By combining multiple recognition factors, the discrimination degree of anomaly identification is improved, which is suitable for anomaly detection in different geological environments.

[0051] Global traversal calculation to ensure the integrity of identification results: Traversing all trace points in turn to ensure the complete calculation of geological anomaly recognition factors in the entire study area and avoid missing important geological information.

[0052] Application value: This method is applicable to multiple fields such as oil and gas exploration, fault detection, volcanic rock identification, and groundwater resource investigation. It can effectively improve the accuracy of seismic data interpretation and provide strong technical support for the accurate identification of geological anomalies. Description of the drawings

[0053] Figure 1 It is a certain inline profile of the flattened seismic data of a certain work area in the embodiment;

[0054] Figure 2 In the embodiment Figure 1 The corresponding data subset (unfolded according to the line number) of the vertical black dotted line in it;

[0055] Figure 3a In the embodiment Figure 2 The initial time shift amount curve graph of each other trace and the target trace in the data subset;

[0056] Figure 3b In the embodiment Figure 2 The last iteration time shift amount curve graph of each other trace and the target trace in the data subset;

[0057] Figure 3c In the embodiment Figure 2 The final time shift data of the data subset;

[0058] Figure 4a In the embodiment Figure 3c The correlation coefficient curve graph of each other trace and the target trace in it;

[0059] Figure 4b In the embodiment Figure 3b The correlation coefficient curve graph after weighting the correlation coefficients and sorting them from large to small;

[0060] Figure 5a It is the horizon slice of the target horizon of the original 3D data in the work area of the embodiment;

[0061] Figure 5b It is the standard slope S of the geological anomaly identification factor of the target horizon of the original 3D data in the work area of the embodiment m,l ;

[0062] Figure 5c It is the fitting slope F of the geological anomaly identification factor of the target horizon of the original 3D data in the work area of the embodiment m,l ;

[0063] Figure 5d It is the reciprocal of the mean E of the geological anomaly identification factor of the target horizon of the original 3D data in the work area of the embodiment m,l ;

[0064] Figure 5eMedian reciprocal U of the geological anomaly identification factor for the target horizon in the original 3D data of the work area in the embodiment m,l . Detailed implementation manners

[0065] The following describes the detailed implementation manners of the present invention in conjunction with embodiments:

[0066] It should be noted that the structures, ratios, sizes, etc. shown in this specification are only used to cooperate with the content disclosed in the specification for those skilled in this technology to understand and read, and are not used to limit the limiting conditions under which the present invention can be implemented. Any modification of the structure, change of the proportional relationship or adjustment of the size, without affecting the effects that the present invention can produce and the purposes that can be achieved, should still fall within the scope covered by the technical content disclosed in the present invention.

[0067] At the same time, the terms such as "upper", "lower", "left", "right", "middle" and "one" cited in this specification are only for the convenience of clear narration, and are not used to limit the scope under which the present invention can be implemented. The change or adjustment of their relative relationship, without substantial change in the technical content, should also be regarded as the scope under which the present invention can be implemented.

[0068] Embodiment 1:

[0069] A high-precision geological anomaly identification method guided by seismic facies, the method comprising:

[0070] 1. Data input and preprocessing

[0071] Input: 3D seismic data (including multiple survey lines, each survey line includes multiple traces, and each trace includes multiple sampling points). Target horizon data (a two-dimensional matrix representing the time information of the target horizon). Sampling interval (the time or depth difference between adjacent sampling points).

[0072] Processing: According to the target horizon data, flatten the seismic data in layers to align the reflection signals of the same geological layer and reduce the influence of formation undulations.

[0073] Output: Seismic data flattened in layers.

[0074] 2. Initialization of iterative calculation

[0075] Input: Calculation parameters, including data subset size, cross-correlation calculation window size, maximum number of iterations, correlation coefficient threshold, Euclidean distance weight parameter, etc.

[0076] Processing: Set all elements of the initial time shift matrix to 0 according to the data subset size parameter.

[0077] Output: Calculation parameters and the initialized time shift matrix.

[0078] 3. Extract the local data subset of the target trace

[0079] Input: Flattened seismic data. Coordinate information (trace number and sample number) of the target trace. Calculation parameters (data subset size, cross - correlation calculation window size, etc.).

[0080] Processing: Extract the 3D data subset around the target trace, calculate the Euclidean distance between each trace in the data subset and the target trace, which serves as the basis for subsequent weight calculation.

[0081] Output: The 3D data subset corresponding to the target trace. The calculated Euclidean distance matrix.

[0082] 4. Calculate the relative time shift between the data subset and the target trace

[0083] Input: The 3D data subset of the target trace points. Calculation parameters (cross - correlation calculation window size).

[0084] Processing: Calculate the cross - correlation function between each trace in the data subset and the target trace, and find the time shift corresponding to the maximum value of the cross - correlation function, which represents the time difference of the target horizon of each trace in the data subset relative to the target trace.

[0085] Output: The time shift matrix of each trace in the data subset relative to the target trace.

[0086] 5. Iteratively optimize the time shift

[0087] Input: The time shift matrix, calculation parameters (maximum number of iterations, convergence criterion, etc.).

[0088] Processing: Judge whether the time shift of all points is less than 1. If satisfied, stop the iteration; otherwise, update the time shift data of the data subset and calculate the cross - correlation of the time shift data. Continue the iteration until the time shift converges or reaches the maximum number of iterations.

[0089] Output: The final optimized time shift matrix. The final time shift result corresponding to the data subset of the target trace.

[0090] 6. Calculate the weighted correlation coefficient

[0091] Input: The final time shift result corresponding to the data subset of the target trace, the Euclidean distance matrix of the target trace, calculation parameters (Euclidean distance weight).

[0092] Processing: Calculate the correlation coefficient matrix according to the final time shift result corresponding to the data subset of the target trace, and use the Euclidean distance matrix to weight the correlation coefficient matrix so that data points with closer distances have a greater impact. Sort the elements of the weighted correlation coefficient matrix from large to small and filter according to the threshold to remove data elements with lower correlation.

[0093] Output: The screened data elements and their quantities.

[0094] 7. Calculate the geological anomaly identification factor

[0095] Input: The screened data elements and their quantities.

[0096] Processing: Calculate the standard slope to measure the change trend of the correlation coefficient. Calculate the fitting slope to describe the change law of the correlation coefficient by the method of curve fitting. Calculate the reciprocal of the mean and the reciprocal of the median to measure the average distribution and the middle trend of the data respectively.

[0097] Output: The standard slope, the fitting slope, the reciprocal of the mean, the reciprocal of the median.

[0098] 8. Calculate the combined identification factor

[0099] Input: The calculated identification factors (standard slope, fitting slope, reciprocal of the mean, reciprocal of the median). Calculation parameters (weight factors, etc.).

[0100] Processing: Combine different identification factors to calculate multiple combined identification factors to improve the ability to identify geological anomalies.

[0101] Output: Multiple combined identification factors.

[0102] 9. Traverse all data points and calculate the identification factor of the complete area

[0103] Input: All parameters, data and identification factor calculation methods in the calculation process.

[0104] Processing: Traverse all survey lines and trace points in turn, and repeat steps 2 - 8 to ensure that the geological anomaly identification factors of the entire work area are calculated.

[0105] Output: The distribution map of the geological anomaly identification factor of the whole area.

[0106] Final output: The calculated geological anomaly identification factor, which can be used for geological analysis applications such as oil and gas exploration, fault detection, volcanic rock identification, etc.

[0107] Example 2:

[0108] The specific steps of the present invention are as follows:

[0109] S1. Obtain the 3D seismic data P of a certain work area, the 2D data T of the target horizon and the time sampling interval Δt. Regard the data P as a 3D matrix of N×M×L, where L is the total number of lines, M is the number of traces per line, and N is the number of sampling points per trace. Then T is a 2D matrix of M×L;

[0110] S2. Use the horizon data T to flatten the original data P to obtain data Q. Each sampling point of the flattened data has N 1 , and the corresponding in-phase axis after flattening is located at the sampling point n 1 ;

[0111] Figure 1 is an inline profile of the flattened seismic data in a certain work area in the embodiment. It can be seen that the in-phase axis is not completely flattened in some local parts, especially at the fault;

[0112] S3. Initialize the line number m = 1 and the trace number l = 1, set the window size N for cross-correlation coefficient calculation 2 , the initial cross-correlation iteration times k = 1, the correlation coefficient threshold λ, and the Euclidean distance weighting parameter α. Determine the number of lines J and the number of traces I on both sides included in the three-dimensional data subset of the target trace according to the seismic facies;

[0113] S4. Generation of three-dimensional data subset: Let the target trace be the data of the m-th trace on the l-th line of data Q, then the corresponding three-dimensional data subset in Q is A m,l , which is a three-dimensional matrix of N 2 ×(2I + 1)×(2J + 1). Calculate the Euclidean distance O m,l between each trace data in the three-dimensional data subset A m,l and the target trace data;

[0114] Figure 2 is the data subset corresponding to the vertical black dotted line in the embodiment Figure 1 , where I = J = 5, the data subset has a total of 121 trace data, and the target trace data is in the 61st trace;

[0115] S5. According to the set window size M 2 for cross-correlation coefficient calculation, calculate the time shift amount D m,l,k corresponding to the maximum cross-correlation value between each trace data in the k-th three-dimensional data subset and the target trace;

[0116] S6. Statistically record the sum of the time shift amounts of the first k times as W m,l , and then judge whether the time shift amount D m,l,k of the k-th time is all less than 1. If so, calculate the correlation coefficient C m,l between each trace data and the target trace when the time shift amount is W m,l , otherwise k = k + 1, and shift the three-dimensional data subset A m,l according to W m,l , and repeat step S5 to obtain the time shift amount D m,l,k at this time, until the time shift amount D m,l,k of the k-th time is all less than 1;

[0117] Figure 3a and Figure 3bThey are respectively the initial and final time shift amount curves of each other trace and the target trace in the data subset of Figure 2 . Figure 3c It is the result of the final time shift of the data subset in Figure 2 . It can be seen that the flattening accuracy of the time shift data is significantly improved, laying a good data foundation for calculating the accurate correlation coefficient.

[0118] S7. Use the Euclidean distance O m,l to weight the final correlation coefficient C m,l to obtain G m,l , and sort G m,l from large to small. Then, according to the correlation coefficient threshold, screen the correlation coefficients to obtain one-dimensional data H m,l . The number of the screened correlation coefficients is K m,l .

[0119] Figure 4a It is the correlation coefficient curve of each other trace and the target trace in Figure 3c . Figure 4b For Figure 3b , the correlation coefficient curve after weighting the correlation coefficients and sorting them from large to small (α = 0.0).

[0120] S8. According to K m,l and G m,l , calculate the standard slope S m,l , fitting slope F m,l , mean reciprocal E m,l and median reciprocal U m,l of the geological anomaly body recognition factor, as well as the combined recognition factors and

[0121] Figure 5a . 5b 5c, 5d, and 5e are respectively the target horizon of the original 3D data in this work area, and the standard slope S m,l , fitting slope F m,l , mean reciprocal E m,l and median reciprocal U m,l of the geological anomaly body recognition factor corresponding to the target horizon. It can be seen that all four recognition factors can identify the geological anomaly bodies on the target horizon, and the accuracy and resolution of the mean reciprocal E m,l and median reciprocal U m,l are significantly higher than those of the standard slope S m,l and fitting slope F m,l , which is beneficial to the identification and prediction of the reservoir on the target horizon in the next step.

[0122] S9. Sequentially increase the line number l and trace number m to complete the calculation of the geological anomaly body recognition factor for all traces of the layer-flattened data.

[0123] Further, the flattened data Q described in step S2 is expressed as:

[0124] Q[n - T[m][l]][m][l] = P[n][m][l];

[0125] where n ∈ [0, N - 1], m ∈ [1, M], l ∈ [1, L].

[0126] Further, the three - dimensional data subset A corresponding to the target trace in step S4 m,l and the Euclidean distance O m,l are:

[0127] A m,l [n 2 [x][y] = Q[n 1 + n 2 [m - I - 1 + x][l - J - 1 + y];

[0128]

[0129] where n 2 ∈ [-N 2 / 2, N 2 / 2], x ∈ [1, 2I + 1], y ∈ [1, 2J + 1].

[0130] Further, the time - shift amount data described in step S5 is calculated by the following formula:

[0131]

[0132] where n 3 ∈ [-N 2 + 1, N 2 - 1].

[0133] Further, the sum W of the time - shift amounts for the first k times described in step S6 m,l is calculated by the following formula:

[0134] W m,l [x][y] = ∑D m,l,k [x][y];

[0135] According to W m,l the three - dimensional data subset A after time - shift m,l is updated to:

[0136] A m,l [n 2 [x][y] = Q[n 1 + n 2 + W m,l[x][y]][m-I-1+x][l-J-1+y]

[0137] The final correlation coefficient C m,l is as follows:

[0138]

[0139] Furthermore, the data G described in step S7 m,l is calculated by the following formula:

[0140]

[0141] Furthermore, the standard slope S of the geological anomaly identification factor described in step S8 is calculated by the following formula:

[0142]

[0143] The fitting slope F is calculated by the following formula:

[0144]

[0145] where k 1 ∈[1, K m,l , and the reciprocal of the mean E is calculated by the following formula:

[0146]

[0147] The reciprocal of the median U is calculated by the following formula:

[0148]

[0149] The reciprocal of the median U, the reciprocal of the mean E, and the combined identification factor are calculated by the following formula:

[0150]

[0151] Those skilled in the art should understand that the embodiments of the present invention can be provided as a method, a system, or a computer program product. Therefore, the present invention can take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present invention can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0152] The present invention is described with reference to the flowcharts and / or block diagrams of methods, apparatuses (systems), and computer program products according to embodiments of the present invention. It should be understood that each flow and / or block in the flowcharts and / or block diagrams, and combinations of flows and / or blocks in the flowcharts and / or block diagrams can be implemented by computer program instructions. These computer program instructions can be provided to the processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing device to produce a machine, such that the instructions executed by the processor of the computer or other programmable data processing device generate means for implementing the functions specified in one flow Figure 1 one flow or multiple flows and / or blocks Figure 1 or multiple blocks.

[0153] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing device to work in a specific manner, such that the instructions stored in the computer-readable memory produce a manufactured article including instruction means that implement the functions specified in one flow Figure 1 one flow or multiple flows and / or blocks Figure 1 or multiple blocks.

[0154] These computer program instructions can also be loaded onto a computer or other programmable data processing device, such that a series of operational steps are executed on the computer or other programmable device to produce a computer-implemented process, so that the instructions executed on the computer or other programmable device provide steps for implementing the functions specified in one flow Figure 1 one flow or multiple flows and / or blocks Figure 1 or multiple blocks.

[0155] The preferred embodiments of the present invention are described in detail above, but the present invention is not limited to the above embodiments. Within the knowledge of those of ordinary skill in the art, various changes can be made without departing from the spirit of the present invention.

[0156] Many other changes and modifications can be made without departing from the concept and scope of the present invention. It should be understood that the present invention is not limited to specific embodiments, and the scope of the present invention is defined by the appended claims.

Claims

1. A high-precision geological anomaly identification method guided by seismic phases, characterized in that: The method comprises: S1: Acquire 3D seismic data, target layer data and sampling interval, flatten the data to align the layers, and obtain seismic data after flattening; S2: Taking each data track of the seismic data after layer flattening as the target track, extracting the 3D data subset around the target track, and calculating the Euclidean distance between the center point of each track in the data subset and the center point of the target track, to obtain the 3D data subset and Euclidean distance matrix corresponding to the target track; S3: Calculate the cross-correlation between each data track in the data subset and the target track, determine the time shift corresponding to each data track, represent the time difference of each data track layer time in the data subset relative to the target track, and obtain the time shift matrix of each data track in the data subset relative to the target track; S4: Based on the time shift matrix, determine whether the absolute values ​​of the time shifts of all elements in the matrix are less than 1. If so, stop the iteration; otherwise, update the time shift results of the data subset, recalculate the cross-correlation, and continue iterating until the time shift matrix meets the iteration termination condition or reaches the maximum number of iterations, and output the optimized final time shift matrix and the time shift results of the data subset; S5: Calculate the correlation coefficient matrix of the target track based on the time-shift result of the final data subset, weight the correlation coefficient matrix according to the Euclidean distance matrix, sort all elements of the weighted correlation coefficient matrix from large to small, set a threshold for screening, remove elements with low correlation coefficients, and obtain screened data elements and their number; S6: using the output result of step S5, calculate the standard slope, the fitting slope, the inverse of the mean and the inverse of the median to obtain each identification factor; S7: combining different recognition factors, calculating multiple combined recognition factors, and obtaining multiple combined recognition factors; S8: Traverse all the data of the layer-leveling seismic data, calculate the identification factor of the target layer in the whole area, and obtain the distribution map of the geological anomaly identification factor of the target layer in the whole area.

2. The high-precision geological anomaly identification method guided by seismic phase according to claim 1 is characterized in that: In step S1, the flattened data Q is expressed as: Q[nT[m][l]][m][l]=P[n][m][l]; Among them, n∈[0,N-1],m∈[1,M],l∈[1,L].

3. The high-precision geological anomaly identification method guided by seismic phase according to claim 1 is characterized in that: In step S2, the three-dimensional data subset A corresponding to the target point m,l and Euclidean distance O m,l for: A m,l [n2][x][y]=Q[n1+n2][m-I-1+x][l-J-1+y]; Among them, n2∈[-N2 / 2,N2 / 2],x∈[1,2I+1],y∈[1,2J+1], Δx and Δy are the spatial sampling intervals in the inline and crossline directions respectively.

4. The high-precision geological anomaly identification method guided by seismic phase according to claim 1 is characterized in that: In step S3, the time shift data is calculated by the following formula: Among them, n3∈[-N2+1,N2-1], argmax() represents the function of calculating the maximum value position.

5. The high-precision geological anomaly identification method guided by seismic phase according to claim 1 is characterized in that: In step S4, the sum of the time shifts of the previous k times W m,l Calculated by the following formula: W m,l [x][y]=∑D m,l,k [x][y]; According to W m,l Three-dimensional data subset A after time shift m,l Updated to: A m,l [n2][x][y]=Q[n1+n2+W m,l [x][y]][m-I-1+x][l-J-1+y]; The final correlation coefficient C m,l for:

6. The high-precision geological anomaly identification method guided by seismic phase according to claim 1 is characterized in that: In step S5, the weighted correlation coefficient G m,l Calculated by the following formula: Where α≥0 is the weighting parameter.

7. The high-precision geological anomaly identification method guided by seismic phase according to claim 1 is characterized in that: In step S6, the standard slope S of the geological anomaly identification factor is calculated by the following formula: The fitting slope F is calculated by the following formula: where k1∈[1,K m,l ], the inverse of the mean E is calculated by the following formula: The median reciprocal U is calculated by the following formula: The median inverse U, the mean inverse E and the combination identification factor are calculated by the following formula:

Images