Seismic event extraction method based on B-COSFIRE filter of differential evolution algorithm

Abstract

The present invention discloses a seismic event extraction method based on a B‑COSFIRE filter of a differential evolution algorithm. In the method of the present invention, only two values ​​are used to describe each pixel, one from a strip structure selection filter and the other from a strip structure end selection filter, and they are merged by summing. The B‑COSFIRE filter proposed by the present invention is universal because its selectivity is not predefined when implemented, but is determined by a prototype pattern given in an automatic configuration process. The weighted geometric mean used in the present invention is a nonlinear function that responds only when filters for all points of interest are present. This type of function is more robust to noise than other methods. The present invention can automatically select strip structures of different widths through a differential evolution algorithm, has a high degree of flexibility, and can pick up events of different widths. The present invention achieves significant extraction of event target contours, preparing for subsequent seismic interpretation, reservoir prediction, and the like.

Description

Technical Field

[0001] The present invention belongs to the technical field of seismic exploration, relates to a combined shift filter response model (B-COSFIRE) filter, and uses a differential evolution algorithm to determine a filter subset to maximize the performance of an event extraction task. The present invention is a seismic data event detection method based on ganglion cell receptive field characteristics (B-COSFIRE) and a differential evolution algorithm. Technical Background

[0002] In the process of seismic exploration image processing, the stratigraphic section is structurally interpreted, and the position description of the detection point of the seismic image is the content of the structural interpretation of the seismic section. Since the visual cognitive process in seismic data processing often follows the structural influence of the event picking of the seismic wave reflection section, the position of effective oil and gas in the stratigraphic layer in the seismic section is obtained by picking up the visual information of the event of the seismic section, and then the seismic data processing and interpretation is realized through the effective stratigraphic layer information to improve the efficiency of seismic exploration.

[0003] In seismic exploration, most of the information in the collected signal is contained in the event picking, so it is meaningful to discuss how to realize visual feature extraction, such as: the self-excitation time of the event can reflect the depth information of the interface to a certain extent, and the shape of the event reflects the propagation speed information of the seismic wave in the medium. Since the seismic wavelet changes of the event are derived from the conditions such as the seismic wave propagation process, the detection and picking of the event are very important for the processing and interpretation of seismic data. The usual way to obtain reflected seismic data is: seismic waves are excited by artificial seismic sources (commonly explosives). During the propagation of seismic waves, when seismic waves encounter rock interfaces with different medium properties, part of the energy will be reflected back to the ground. The computational model with visual features receives and obtains multiple seismic records (signals), and each seismic signal corresponds to a detector. When seismic exploration seismic waves reflect multiple reflection interfaces underground, each seismic signal will generate multiple reflected seismic wavelets. At this time, the line connecting the extreme values ​​(usually called peaks or valleys) of the vibration profile phase on the seismic record is called the event.

[0004] Discussing the event axis information in seismic data, various detection technologies have been developed for seismic exploration, mainly including:

[0005] (1) Edge detection method: It converts the amplitude value of each sampling point of the two-dimensional seismic signal into different grayscale values, and regards the seismic wave gather as a grayscale image, so that the edge detection method can be used in image processing to detect the phase axis in the image. The event axis picked up by this method is relatively fuzzy, the resolution is greatly affected by space and frequency interference, and is greatly affected by noise factors. The detection result is the envelope of the grayscale mutation area, with a very low resolution, and it is often not directly used as the event axis detection result.

[0006] (2) Neural network method: It uses known epochs as standard samples to simulate the neural network and uses the error return method to gradually modify the connection weights between neurons. The process of the neural network method requires sufficient sample support, but the selection of samples is more difficult, and the adoption of new samples will also affect the neural network. A large number of iterative operations will lead to longer processing time.

[0007] (3) Cross-correlation method: The basic principle is to use the waveform similarity characteristics of the event axes between seismic traces to extract the event axes that are easily affected by noise. The spatial resolution decreases with the increase of calculation traces, which increases the amount of calculation.

[0008] (4) Chain matching method: In this method, each waveform is represented by a link of peaks and valleys with multiple features, and then the picking of the in-phase axis is converted into the optimal matching problem between chains, that is, finding the smallest link path. This method is also affected by noise, space and frequency interference, so it is also difficult to solve complex in-phase axis processing problems.

[0009] (5) Tu et al. proposed an automatic event-picking system: This method includes three steps: two-dimensional matched filter, Kalman filter and flexible template.

[0010] In the process of visual cognition, according to the influence of event picking of seismic wave reflection profile on structural interpretation, the position of effective oil and gas in the stratigraphic layer is obtained by event picking of seismic profile, and then the processing and interpretation of seismic data are realized through effective stratigraphic layer information to improve the efficiency of seismic exploration. Simulating biological visual recognition mechanism to realize seismic data target contour extraction will greatly improve the precise connection between seismic exploration and artificial intelligence, which is the hot spot and difficulty of current research. Summary of the invention

[0011] In order to solve the problems existing in the prior art, the purpose of the present invention is to propose a seismic event extraction method based on the B-COSFIRE filter of the differential evolution algorithm, which is used for seismic data event picking technology. The picked events can be used for sequence stratigraphy research, reservoir prediction and storage feature extraction.

[0012] In order to achieve the above object, the technical solution adopted by the present invention is:

[0013] The seismic event extraction method of B-COSFIRE filter based on differential evolution algorithm includes the following steps:

[0014] Step 1): Acquire seismic data and preprocess the seismic data;

[0015] Step 2): Use the two-dimensional expression of the DoG model to simulate the receptive field model, i.e., the B-COSFIRE filter, wherein the DoG model can realize an inhibitory field with an excitatory central region and an inhibitory surround;

[0016] Step 3): setting two groups of B-COSFIRE filters, combining the responses of the two groups of B-COSFIRE filters together for event picking of seismic data, wherein one group of the two groups of B-COSFIRE filters is used for detecting the event, and the other group is used for detecting the end of the event;

[0017] Step 4): Use the preprocessed seismic data as the input of the B-COSFIRE filter, blur and shift the DoG response, and obtain the weighted geometric mean to obtain the output of the B-COSFIRE filter;

[0018] Step 5): Based on the output results of step 4), the differential evolution algorithm is used as a feature selection method to evaluate the impact on system performance, determine the minimum feature subset that maximizes the event depiction performance, and complete the seismic event extraction.

[0019] Furthermore, the seismic data is in sgy format.

[0020] Furthermore, the preprocessing is specifically: wavelet fusion of the envelope peak instantaneous frequency of the pre-stack seismic data and the inclined stack peak instantaneous amplitude.

[0021] Furthermore, the two-dimensional expression of the DoG model is as follows:

[0022]

[0023] The corresponding value of the central area with excitability is set to 0.5σ, where σ is the spatial distribution of the Gaussian function of the surrounding inhibitory domain. It is an On-centered receptive field, that is, the center is the excitation field and the periphery is the inhibition field. It is the Off center receptive field, defined as:

[0024] For a given location (x, y) and an intensity distribution I(x', y') of seismic data I, the response c of the DoG is σ (x,y) is calculated as follows:

[0025]

[0026] Among them, * is convolution, |·| + It is half-wave rectification.

[0027] Furthermore, the DoG response is blurred and shifted, and the weighted geometric mean is calculated to obtain the output of the B-COSFIRE filter, which is specifically:

[0028] With a length of ρ i The vector moves each ambiguous DoG response toward the center of the support region, and the associated change vector is (Δx i ,Δy i ),Δx i =-ρ i cosφ i , Δy i =-ρ i sinφ i , in φ i Calculate the tuple (σ in the reverse direction i ,ρ i ,φ i ) and the response of DoG and shift it as shown in the following formula:

[0029]

[0030] In the formula, σ i is the standard deviation of the DoG model, (ρ i ,φ i ) is the polar coordinate, G σ' (x', y') is the maximum weighted response of DoG, -3σ'≤x', y'≤3σ', where σ'=σ'0+αρ i , the values ​​of σ'0 and α are constants.

[0031] Furthermore, the weighted geometric mean in step 4) is specifically obtained as follows:

[0032]

[0033] In the formula, S is the output response of the B-COSFIRE filter. s The weighted geometric mean is an AND-type function, that is, a B-COSFIRE filter is only valid if all incoming blur and shift responses A response is achieved only when it is greater than zero.

[0034] Furthermore, the differential evolution algorithm performs iterative search using mutation operations, crossover operations and selection operations.

[0035] Furthermore, in the mutation operation, a local fitness function is used to guide the mutation direction of each sub-component.

[0036] Furthermore, in the crossover operation, a uniform crossover method is adopted.

[0037] Compared with the prior art, the present invention has at least the following beneficial effects:

[0038] The present invention proposes a seismic data event axis detection method based on the receptive field characteristics of ganglion cells (B-COSFIRE) and a differential evolution algorithm. The method is a trainable filter method that adopts an unsupervised approach. Compared with the supervised method based on machine learning technology, in the method of the present invention, we only use two values ​​to describe each pixel, one from the strip structure selection filter and the other from the strip structure end selection filter, and merge them by summing. This provides a very effective method for picking up event axes. However, as more operations need to be performed, the time complexity increases with the increase in the number of dimensions. The B-COSFIRE filter proposed by the present invention is universal because its selectivity is not predefined when implemented, but is determined by a prototype pattern given in the automatic configuration process. The weighted geometric mean used in the present invention is a nonlinear function that only responds when the filters of all points of interest exist. This type of function is more robust to noise than other methods. The present invention can automatically select strip structures of different widths through a differential evolution algorithm, has high flexibility, and can pick up event axes of different widths. The present invention realizes the significant extraction of the event axis target contour, making preparations for subsequent seismic interpretation, reservoir prediction, etc.

[0039] Furthermore, the seismic data is preprocessed before being processed by B-COSFIRE filter, wherein the preprocessing is wavelet fusion of the envelope peak instantaneous frequency and the tilt stack peak instantaneous amplitude of the prestack seismic data, which can provide the best quality data to the maximum extent.

[0040] Further, a new B-COSFIRE filter is configured at the end of the strip structure to achieve the end response. The present invention uses the point at the end of the strip structure as the point of interest. Compared with the symmetrical B-COSFIRE filter, the asymmetrical B-COSFIRE filter obtains a stronger response at the end of the strip structure.

[0041] Furthermore, the weighted geometric mean is an AND-type function, i.e., a B-COSFIRE filter can achieve a response only when the values ​​of all incoming blur and shift responses are greater than zero. The contribution of blur and shift responses decreases as the distance from the support center of the B-COSFIRE filter increases.

[0042] Furthermore, the versatility of B-COSFIRE filters lies in their trainability, which can automatically configure the optimal parameters for the filters according to the objectives and obtain the best picking results.

[0043] Furthermore, the B-COSFIRE filter proposed in the present invention is nonlinear because it achieves directional selectivity by multiplying the outputs of a set of Gaussian difference (DoG) filters whose supports are aligned in a colinear manner. It can tolerate changes in rotation and slight deformations. In addition, COSFIRE is a trainable filter method. This means that the selectivity of the filter is not predefined in the implementation, but is determined by the prototype mode (such as straight container, bifurcation or intersection) specified by the user during the automatic configuration process.

[0044] Furthermore, a differential evolution algorithm is used to construct a feature vector based on the response of the B-COSFIRE filter. The differential evolution algorithm can evaluate the combined contribution of many features and explore a larger solution space.

[0045] Furthermore, the seismic data of the present invention is in sgy format, which contains a large amount of information. The present invention is unique in that it applies visual cognition methods to the earthquake field. Instead of converting the seismic data into a grayscale image, the data itself is used. More attribute information can be extracted from the data. The present invention uses seismic data in sgy format and extracts attributes such as instantaneous frequency and instantaneous amplitude therefrom, and then performs a co-phase axis picking operation to maximize the resolution and retain information.

[0046] Furthermore, the present invention performs event picking based on pre-stack and post-stack seismic data. Pre-stack seismic data contains more abundant information than post-stack seismic data, including common shot point gathers, VSP data, etc. acquired from seismic data; and it has not been stretched and stacked by NMO, and the amplitude and frequency information are not damaged, and it can reflect some subtle stratigraphic features, but its disadvantages are low signal-to-noise ratio and high noise. Therefore, the present invention proposes a high-precision event picking method based on the low signal-to-noise ratio of pre-stack seismic data. Of course, the method of the present invention is also applicable to event picking of post-stack seismic data. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 Schematic diagram of Radon transform.

[0048] Figure 2 An example of a B-COSFIRE filter configuration using a horizontally synthesized strip structure, where (a) is the schematic diagram of the B-COSFIRE filter, (b) is the response of the DoG filter, (c) is the local maximum DoG response along the concentric circles of the target point, and (d) is a sketch of the final filter.

[0049] Figure 3 Schematic diagram of the B-COSFIRE filter.

[0050] Figure 4Examples of B-COSFIRE filter configurations using vertically synthesized stripe structure ends, where (a) is the prototype of the stripe structure end and (b) is the asymmetric B-COSFIRE filter configured using (a).

[0051] Figure 5 are actual seismic data, where (a) is the original pre-stack CSP data; (b) is the instantaneous amplitude profile of (a).

[0052] Figure 6 is the instantaneous frequency of the envelope peak.

[0053] Figure 7 It is the wavelet fusion profile of the envelope peak instantaneous frequency and the tilt-stacked peak amplitude.

[0054] Figure 8 is the dip profile obtained by Radon transform.

[0055] Fig. 9 This is the B-COSFIRE filter response result.

[0056] Fig.10 Results of picking the event for the B-COSFIRE filter. DETAILED DESCRIPTION

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

[0058] Aiming at extracting the contour and position features of the event axis of the seismic profile, the present invention proposes a seismic data event axis detection method based on the receptive field characteristics of optic ganglion cells (B-COSFIRE) and differential evolution algorithm. The method uses the two-dimensional expression of the DoG model to simulate the receptive field model, which can realize the function of an excitatory (i.e. positive) central area and an inhibitory surround (i.e. negative) inhibitory domain. The DoG response is blurred and shifted, and the weighted geometric average is obtained to obtain the B-COSFIRE filter response model. Finally, the differential evolution method is used to optimize the best feature selection method to determine the minimum feature subset that maximizes the event axis depiction performance. The model proposed by the present invention can effectively detect the edge of the event axis of the target contour of the seismic data with huge data volume and numerous noise interference. Compared with other traditional edge detection, the B-COSFIRE filter proposed by the present invention can pick up the events of seismic data of different widths, which is more in line with the biological vision mechanism. The events picked up by this method are more continuous and accurate than those picked up by the existing methods, providing a new research direction for the interpretation of seismic data in the future.

[0059] The line connecting the extreme values ​​(usually called peaks or valleys) of the phase of the vibration profile on the seismic record is called the event axis. Most of the information in the acquired signal in seismic exploration is contained in the event axis picking, so it is very meaningful to discuss how to realize visual feature extraction, such as: the self-excitation time of the event axis can reflect the depth information of the interface to a certain extent, and the shape of the event axis reflects the propagation speed information of the seismic wave in the medium. Since the seismic wavelet changes of the event axis are derived from the conditions such as the seismic wave propagation process, the detection and picking of the event axis are very important for the processing and interpretation of seismic data. Therefore, the present invention proposes a high-precision event axis picking method based on reflection seismic data.

[0060] The present invention proposes a set of B-COSFIRE filters to solve the problem of picking up thin and long structures, that is, picking up the phase axis of seismic data, and the filters can select phase axes of different widths. The B-COSFIRE filter method was originally used to delineate retinal blood vessels. Such filters have also been used to analyze computed tomography angiography (CTA) images. This shows that they are suitable for various applications of picking up and segmenting thin and long structures. The present invention combines the responses of two B-COSFIRE filters, one dedicated to detecting the phase axis and the other dedicated to detecting the end of the phase axis. The parameter selection of the phase axis end filter should ensure that the performance of both filters is maximized. This means that the configuration of the phase axis end filter depends on the phase axis filter. The parameter selection of the filter is very important, and the configuration parameters of each filter are selected to achieve the best performance at the common width of all strip structures. The present invention determines a subset of B-COSFIRE filters through information theory and machine learning, selects strip structures of different widths, and uses differential evolution algorithm as a feature selection method to evaluate the impact on system performance.

[0061] Firstly, the seismic data is preprocessed, which is the wavelet fusion of the envelope peak instantaneous frequency and the inclined stacking peak instantaneous amplitude of the prestack seismic data; then the B-COSFIRE filter principle and the feature selection process based on the differential evolution algorithm are introduced.

[0062] The preprocessing process is briefly described below. In the embodiment of the present invention, it is assumed that the source wavelet is approximated by the following constant phase wavelet with 4 parameters:

[0063]

[0064] In the formula, is the modulation angular frequency, δ is the energy attenuation factor, A' and φ are the amplitude and phase respectively, t is the time, Performing Fourier transform on both sides of equation (1) yields:

[0065]

[0066] Where ω is the angular frequency and φ is the phase.

[0067] According to Barnes' definition, the envelope peak instantaneous frequency (EPIF) of a constant phase wavelet is the average frequency of the wavelet amplitude spectrum weighted by the Fourier frequency, that is,

[0068]

[0069] Where τ is time, f is frequency, and S(τ,f) is the Fourier transform of the source wavelet at time τ.

[0070] Substituting the modulus of equation (2) into equation (3), we can obtain the EPIF of the source wavelet at time zero:

[0071]

[0072] The tilt stack peak amplitude is calculated based on the instantaneous amplitude. The calculation formula for the instantaneous amplitude is as follows:

[0073] A(t)=|s(t)+i·H[s(t)]| (5)

[0074] in, is the analytical part of the signal calculated by wavelet transform, A(t) is the modulus of the instantaneous amplitude, s(t) is a seismic data, S(b,a) is the wavelet transform of the seismic data, is the wavelet function, g(t) is the real part of Fourier transform, a is the scale factor, b is the translation factor, and the wavelet transform of s(t) with respect to g(t) is defined as:

[0075]

[0076] Where t, b∈R, a>0; g(t)∈L 1 (R,dt)∩L 2 (R,dt),L 1 Indicates L 1 norm, which is the sum of the absolute values ​​of the elements in the vector, is the complex conjugate of g(t).

[0077] The oblique superposition is a local Radon transform, and the present invention adopts a local linear Radon transform. The integral path of the linear Radon transform is linear, also known as the τ-p transform, and its discrete form in the frequency domain is:

[0078]

[0079] Where p is the slope of the line, M(f,p)=∫m(τ,p)e -j2πfτ dτ, is the instantaneous amplitude A(t,x m ), x m is the mth reference channel, nx is the number of channels to be superimposed, and m(τ,p) is the Radon transform in the time domain. The calculation process of the local linear Radon transform is as follows: Figure 1 As shown: Select a reference trace on the instantaneous amplitude profile, for a certain time intercept τ on the reference trace j , move several tracks near the reference track along n p with different slopes p j (j=1,2,…,n p ) are superimposed with straight lines (the slope is sampled at intervals of Δp), and the sum of the instantaneous amplitude at the time intercept along different directions is calculated, and the superimposed value is recorded at the corresponding position (τ) of the τ-p coordinate axis. j ,pj), when the selected stacking slope is close to or equal to the slope of the event axis, the stacking value recorded in the tx domain along the straight line is the largest. The average value of the maximum stacking value is placed at the corresponding position (τ) in the tx domain (time-distance domain) j ,x m ), a super gather can be constructed to increase the signal-to-noise ratio, and the profile is called the inclined stack peak amplitude profile. There are many methods for calculating the linear Radon transform, and the present invention adopts a high-precision frequency-space domain matrix multiplication method for calculation.

[0080] So far, the envelope peak instantaneous frequency and the tilt stack peak amplitude (both are two-dimensional images in cross-section form) have been obtained. Next, the two will be fused by wavelet. The two-dimensional image fusion based on wavelet transform is mainly divided into three steps:

[0081] (1) Perform wavelet decomposition on each original image to be fused and obtain decomposition coefficients in different directions.

[0082] (2) The coefficients in different directions obtained by decomposition are fused.

[0083] (3) Performing inverse wavelet transform on the fused coefficients can obtain the final result image.

[0084] Next, we will find the B-COSFIRE filter. The human visual nervous system has evolved to receive and recognize external scenes with extremely high resolution. Computer vision studies how to simulate and realize similar visual mechanisms. When the sensory organs receive light stimulation and excitement, the ganglion cells will transmit the sensory information to the central nervous system, and the area where the response is generated is called the classical receptive field. The classical receptive field is a concentric circle receptive field. Outside the concentric circle receptive field, a large periphery is found that has a modulating effect on the classical receptive field. The classical receptive field plus the large periphery is a non-classical receptive field. Neurophysiology divides the brain's visual system into three levels: the retina, the lateral geniculate body, and the visual cortex. The receptive fields include the receptive field of retinal ganglion cells, the receptive field of lateral geniculate body nerve cells, and the classical receptive field of visual cortex cells. Among them, the classical receptive field of visual cortex cells is divided into simple cell receptive field, complex cell receptive field, and super-complex cell classical receptive field. The structure of the classical receptive field is concentric, and the center and the periphery have mutual antagonism. This structure can be described by the Difference of Gaussian (DoG) model. Figure 2 (a) Schematic diagram of the B-COSFIRE filter, configured as a horizontal bar, which serves as the input center DoG filter response at a specific location in the center of its support area. Figure 2 (b) is the response of the DoG filter, Figure 2 (c) is the local maximum DoG response along the concentric circles of the target point (marked by the cross mark in the center), Figure 2 (d) is a sketch of the final filter, where the size of the spots corresponds to the standard deviation of the Gaussian blur function. The following is a one-dimensional expression of the DoG model:

[0085]

[0086] Among them, A1 and A2 are the sensitivities of the central excitation domain and its surrounding inhibitory domain, σ1 is the mean square error of the Gaussian function of the central domain, and σ2 is the mean square error of the Gaussian function of the inhibitory domain.

[0087] COSFIRE is a combination of shifted filter response models. The method that combines filter response combination with strip selection is called B-COSFIRE, which can realize the extraction of strip features. It is sensitive to directional changes and has strong irritability and adaptability to rotation and slight shape changes. The B-COSFIRE filter proposed in the present invention is nonlinear, and its nonlinear response is achieved by multiplying the output of a group of Gaussian difference (DoG) filters to achieve directional selectivity. These filter groups support alignment in a colinear manner, and it can tolerate rotation changes and slight deformations. The following formula gives a two-dimensional expression of the DoG model to simulate the receptive field model, which can realize the function of an inhibitory domain with an excitatory (ie positive) central area and an inhibitory surround (ie negative):

[0088]

[0089] Among them, the corresponding value of the central excitation domain is set to 0.5σ, σ is the spatial distribution of the Gaussian function of the surrounding inhibition domain, that is, the standard deviation. This type of function is a recognized computational model for some cells in the lateral geniculate nucleus (LGN) of the brain. It is an On-centered receptive field, that is, the center is the excitation field and the periphery is the inhibition field. It is the Off center receptive field, defined as:

[0090]

[0091] The above DoG model can simulate the concentric antagonistic structure of cells.

[0092] For a given location (x, y) and an intensity distribution I(x', y') of seismic data I, the response c of the DoG is σ (x,y) can be calculated as follows:

[0093]

[0094] Among them, * is convolution, |·| + Then, we consider the DoG responses around a given point of interest along concentric circles and select the DoG with a local maximum from them.

[0095] The direction of B-COSFIRE is the same as that of the stripe structure. We use three parameters to describe the i-th point: using the DoG model with a set σ, we analyze c within the range of N concentric circles with the specified point as the center. σ (x,y) response. According to the characteristics of the seismic data event axis hyperbola, the B-COSFIRE filter combination is set as follows Figure 3The center of the light gray circle is the center point of the B-COSFIRE filter combination, and the black point generated by the intersection of the middle vertical line and the concentric circle is the extreme point of the DoG response. The parameters of the filter combination are:

[0096] S={(σ i ,ρ i ,φ i )|i=1,…,n} (12)

[0097] S defines a B-COSFIRE filter with a selective preference for a given prototype. i is the standard deviation of the DoG filter, (ρ i ,φ i ) is the polar coordinate, indicating its position, and n is the number of DoG filters that need to be configured. The parameter values ​​set by the following formula are used to Figure 2 (c) is analyzed by the horizontal line form. Figure 2 (c) The input pattern is determined by automatic analysis, where ρ and φ are automatically generated by the points obtained previously:

[0098]

[0099] The input image is first convolved with a DoG function with a specified standard deviation in a tuple of the set S. Next, the DoG filter is blurred and shifted to allow some error in the preferred location of the relevant points. The blur operation is defined as computing the maximum value of a weighted threshold response of the DoG filter. The blur operation computes the DoG maximum weighted response G σ' The maximum DoG response in the local neighborhood of (x',y') is used to improve its fault tolerance. The weights are calculated using the DoG filter response and the Gaussian function c σ' The coefficients of (x', y') are multiplied together, and their standard deviation σ' is the distance from the filter support center ρ i Linear function: σ'=σ'0+αρ i , σ'0 and α are constants that we adjust according to the application, ρ i Indicates the distance between the preferred position of the relevant point and the filter support center.

[0100] Then, the present invention uses a length of ρ i The vector moves each ambiguous DoG response toward the center of the support region, which is φ i The associated change vector is (Δx i ,Δy i ),Δx i =-ρ i cosφ i , Δy i =-ρ i sinφi In φ i Calculate the tuple (σ in the reverse direction i ,ρ i ,φ i ) and the response of DoG and shift it as shown in the following formula:

[0101]

[0102] In the formula, -3σ'≤x',y'≤3σ'.

[0103] Next, calculate the B-COSFIRE filter response model. The B-COSFIRE filter output response r s (x, y) is the weighted geometric mean of all blurred and shifted DoG responses in the corresponding tuple S:

[0104]

[0105] In the formula, (ρ i ,φ i ) is the polar coordinate of the ith DoG response relative to the filter support center, |·| τ represents the threshold of the response at the maximum response at a fraction τ (0 ≤ τ ≤ 1). The weighted geometric mean is an AND-type function, i.e., a B-COSFIRE filter is only effective if all incoming blurred and shifted responses A response is achieved only when it is greater than zero.

[0106] The above program configures the B-COSFIRE filter, which can be used for horizontally oriented containers. In order to achieve multi-directional selectivity, multiple B-COSFIRE filters can be configured by using prototype patterns of different directions. The weighted geometric mean is an "AND function". When all responses satisfy ">0", the response value can be obtained. The farther the specified point is from the center point, the smaller the effect of blurring and shifting will be. It can be seen that the set strip direction is the premise of this filter. We manipulate the parameter ψ of each tuple and create a new set R ψ (S), here we integrate the set S and add R ψ (S) realizes rotational invariance, then:

[0107]

[0108] It represents a B-COSFIRE filter with a directional preference of ψ radians, which is offset from the directional preference of the original filter S. It can track bar targets in all directions, sum all the responses obtained, calculate the response value of each pixel at the (0°, 15°, 30°, ..., 165°) orientation, and select the maximum value as the feature of the point. The present invention selects a group of B-COSFIRE filters with different directional preferences, which can be expressed as:

[0109]

[0110] in, n r =12.

[0111] Because the single reflection wave and multiple reflection waves of seismic data are discontinuous and the symmetry characteristics are not obvious, it is difficult to fully detect the phase axis when picking it. B-COSFIRE is a "logical and" function, which will only respond when all inputs are activated. In principle, it cannot respond at discontinuities. Therefore, consider building a new filter at the single and multiple reflection waves. The filter that previously detects continuous phase axes and direct waves is called B-COSFIRE with a symmetrical structure, and the filter at discontinuities and reflection waves is called B-COSFIRE with an asymmetric structure. The filter will calculate the response within a range of 360°. In other words, the difference between this filter and the symmetrical filter is that it only targets the target point on one side. In actual application, the number of orientations should be doubled to obtain the endpoints of the strip structure in all directions. There is a possibility of target points in each direction, such as Figure 4 As shown, the center of the concentric circles is the response center, and the black dot is the response extreme value.

[0112] The width of the seismic event axis may vary from 1 pixel to several pixels, depending on the offset. The larger the offset, the greater the high-frequency loss, the more low-frequency, and the wider the waveform. Therefore, the present invention configures a set of B-COSFIRE filters, which consists of 21 strip structures {S1, ..., S 21} and 21 bar structure ends {S 22 ,…,S 42}, with optional strip structures of different thicknesses. Large-scale filters are selective for wide events and tolerant to background noise, but have small responses along narrow events. In contrast, small-scale filters are more selective for narrow events but less tolerant to background noise. Their responses are fused together, with the promise of better event extraction performance at various scales. We construct a pixel-level feature vector v(x,y) for each event position (x,y), which consists of the responses of the 42 BCOSFIRE filters in the filter bank, plus the intensity value g(x,y) of the instantaneous amplitude:

[0113]

[0114] in, BCOSFIRE filter S i The rotation tolerance response.

[0115] Prior to classification, the present invention applies an inverse hyperbolic sine transformation function to each element of the feature vector. It reduces the skewness in the data and is defined as follows:

[0116]

[0117] For i and for larger values ​​of θ>0, the application of the function is similar to the logarithmic transformation. Let θ→0, f(υ i ,θ)→υ i The Z-score is then calculated to normalize all 43 eigenvalues. The present invention applies the Z-score normalization process to each image separately to compensate for illumination variations between images.

[0118] The filter bank designed above may have many redundant filters. The present invention studies various feature selection methods to determine the minimum feature subset that maximizes the event delineation performance. The present invention uses training data consisting of a matrix of size N×43 as input, where N corresponds to the number of pixels randomly selected from the training data and the number of columns corresponds to the size of the filter bank plus the EPIF data.

[0119] Inversion optimization technology has been widely used in many fields. In seismic exploration, it can be used to invert the medium parameters of complex underground structures. There are many inversion methods, such as waveform inversion and tomography. Waveform inversion optimization technology includes local optimization inversion methods and global optimization inversion methods. The local optimization method has a fast calculation speed, but the disadvantage is that it is difficult to invert the matrix in high-dimensional parameter inversion, and it is easy to fall into local minima, and the inversion accuracy depends on the initial value. Compared with the local optimization method, the global optimization method has no dependence on the initial value and does not need to calculate the inverse of the matrix. The disadvantage is that the calculation is large and time-consuming. In recent years, due to the significant improvement in computer computing power, the global optimization method has been widely used.

[0120] Global optimization algorithms include simulated annealing, differential evolution, genetic algorithm, co-evolution, and improved differential evolution with collaborative selection (DE-CCS), differential evolution with collaborative mutation (DE-CCM), etc. Evolutionary algorithms are algorithms that use computers to simulate the optimal process of individuals (i.e., species) based on the laws of evolution of various species in nature, the principle of natural selection and survival of the fittest, and through the iteration of the group from generation to generation, genes that are not conducive to the continuation of the species are continuously eliminated, and good genes are passed on, ultimately making the individual optimal. The present invention uses a differential evolution algorithm (DE), uses a local fitness function to guide the selection of the optimal gene for each subcomponent, and then performs inversion to search for the best feature subset in a large number of possible combinations.

[0121] The DE method first gives an initial population X containing M individuals:

[0122] X=[x1,x2,…,x M ] (20)

[0123] The kth individual in the population is represented by:

[0124] x k =[x k,1 ,x k,2 ,…,x k,N ],k=1,…,M (21)

[0125] Where M is the initial number of individuals, which is set according to the actual situation, such as M = 400, and N is the dimension of 43 parameters to be estimated (i.e., EPIF and 42 B-COSFIRE filters). The fitness function calculates the average accuracy of the training data and the selected columns. The local fitness function of the kth individual defined by the L1 norm is:

[0126]

[0127] Where, d o,j (·) is the event axis data selected in the jth column, d k,j (·) is the jth column of training data calculated by the kth individual forward modeling, NR is the number of observed seismic data, and g(l) is the Gaussian window function, which is defined as:

[0128]

[0129] Where l is a variable and σ is the standard deviation, which is used to control the width of the window. The center of the Gaussian window is located at the peak of the wavelet envelope. For post-stack data, the center of the time window should be located at the center of the sub-component. For pre-stack data, as the offset increases, the center of the time window should slide along the hyperbola. kis the parameter of B-COSFIRE and the position of the seismic event axis is picked up. A local fitness value is calculated for each trace. Then the local fitness value vector corresponding to the kth individual is:

[0130] f k =[f k,1 ,f k,2 ,…,f k,N ] (twenty four)

[0131] In the formula, the element f k,1 ,f k,2 ,…,f k,N is the local fitness value of 1 to N channels.

[0132] Therefore, for a population X containing M individuals, the local fitness value corresponding to its forward data is an M×N matrix:

[0133]

[0134] The initial population X is also an M×N matrix.

[0135] The differential evolution algorithm performs iterative search using mutation, crossover, and selection operations. The three main operations are introduced as follows.

[0136] An effective mutation strategy can quickly improve the genes (parameters to be estimated) in individuals, and finally obtain the optimal individual after multiple generations of cycles. In order to find the correct direction of gene mutation, in the mutation operation, the present invention uses the local fitness function to guide the mutation direction of each subcomponent, thereby proposing a new mutation operator. The present invention arranges all variables in the jth column of population X in ascending order:

[0137]

[0138] In the formula, the superscript s represents the sorted sequence. The local fitness value corresponding to this sequence is:

[0139]

[0140] Let the local fitness function corresponding to the sequence before sorting be f j , Depend on The elements in correspond to f j The value in is composed of j Sort the elements in ascending order. Pick up the minimum value of the vector and the NP values ​​in the neighborhood on both sides of the minimum value. The analytical expression of these (2×NP+1) values ​​can be approximated by fitting a concave function polynomial with a minimum value:

[0141]

[0142] In the formula, A, B, r1, and r2 are the four parameters of the composite equation, which can be obtained by the inversion method. x is the variable after sorting the jth column, and y is the corresponding local fitness value. The variable corresponding to the minimum point of the concave function may be a better mutant gene. In order to determine the values ​​of the four parameters A, B, r1, and r2, the objective function is defined as:

[0143]

[0144] In the formula, is the calculated local fitness value, is the value of the polynomial to be fitted. Since there are 4 parameters to be estimated, the inverse of the Hessian matrix is ​​relatively easy to obtain, so the Gauss-Newton method can be used to solve the system of equations.

[0145] However, the global optimization algorithm requires a large number of iterations, and its calculation amount is relatively large. The Gauss-Newton method is an inversion method, and the calculation speed is relatively slow. In order to improve the calculation speed, the present invention adopts the following direct method for mutation. Let the variable corresponding to the minimum value of the local fitness function be x r,j , randomly select two variables x except the minimum value in the jth column of the population r1,j and x r2,j , then the mutation strategy of the jth element of the kth individual is:

[0146]

[0147] In the formula, F is the control parameter, and its value is a random variable between 0 and 1. All individuals cycle to update the genes in the jth column. The mutated genes are disturbed near the optimal value point. The values ​​of different individuals after mutation are slightly different. This is because the randomly selected variables and control parameters are different, which is in line with the laws of nature and reflects the randomness of mutation.

[0148] After all individuals are cycled, in order to make the deep parameter inversion accurate, the present invention constructs a local cycle process, that is, multiple individual cycles are performed. During this local iterative cycle, we observed a phenomenon that as the number of iterations increases, the local fitness value (i.e., error energy) is still very large, but the change between different local fitness values ​​is very small. This is because as the number of mutations increases, this column of parameters gradually approaches the optimal value. When the parameters have not reached the optimal value, the difference between them has become smaller and smaller, and the mutation lacks the motivation to break through, resulting in a small difference in the values ​​before and after the mutation. If no new mutation mechanism is introduced, all components of the jth column in the population will tend to be constants. In order to improve the inversion accuracy, the present invention proposes the following weighted mutation strategy. Given an accuracy threshold of the error energy, when the error energy is greater than the threshold and the difference between the error energies of two individual iterations is less than another threshold, the mutation disturbance amount is amplified by the weighted method of the local fitness value, and the formula is as follows:

[0149]

[0150] In the formula, λ is a control parameter, which can be selected according to the accuracy, the superscript M represents the value after mutation, and the weight coefficient p1 is the normalized local fitness value, which is calculated as follows:

[0151]

[0152] In the formula, f r,j ,f r1,j ,f r2,j is the gene x r,j ,x r1,j ,x r2,j The weighted mutation strategy utilizes the information of error energy, so that the mutation disturbance amount can be adaptively adjusted according to the error energy.

[0153] Crossover operation: Randomly select The corresponding variant Extract genes from the test individuals The two commonly used crossover methods are uniform crossover and exponential crossover. The present invention adopts the uniform crossover method, and the formula is as follows:

[0154]

[0155] where CR∈[0,1) is the crossover rate, n rand It is a random integer between 1 and D, which ensures that at least one gene in the test individual comes from the mutant individual.

[0156] The conventional selection operator is based on the target individual and experimental individuals The global fitness value of will select one of the individuals into the next generation as a whole. The mathematical expression is:

[0157]

[0158] In each iteration, we configure an SVM classifier with a linear kernel using 90% of the training set and apply it to the remaining 10%. After training once on all samples in the training set, we sort the individuals in descending order of fitness and keep only the top 10%. We use this group of elite individuals to generate the remaining 90% of new individuals by performing a crossover operation on randomly selected pairs of elite individuals. Each of the new individuals has a 10% probability of mutating (i.e., changing the code from 1 to 0 or from 0 to 1). We run these iterative steps until all individuals reach the optimum and no longer mutate. Finally, we select the filter corresponding to the position of the individual with the highest fitness value in the population.

[0159] We use the selected features to train a SVM classifier with a linear kernel. The SVM classifier is particularly suitable for binary classification problems because it finds an optimal separating hyperplane that maximizes the margin between classes.

[0160] Case Analysis

[0161] The present invention applies the proposed method to an actual common shot point gather (CSP), which has 595 channels, a minimum offset of 90 m, a distance between detectors of 10 m, 1250 sampling points, and 2 ms sampling. Figure 5 (a) is the original pre-stack seismic data. As can be seen from the figure, the original data has obvious events, but the background noise is also very large, especially the events at close offsets and shallow surface cannot be distinguished, and the peaks and troughs of the events affect each other, which is not conducive to picking up other attribute information such as dip. Figure 5 (b) is the instantaneous amplitude profile. Figure 6 It is the profile of the envelope peak instantaneous frequency (EPIF). As can be seen from the figure, the envelope peak instantaneous frequency reduces the mutual interference between peaks and troughs, and the information is further concentrated, but the deep noise is still very large. Figure 7 This is the wavelet fusion result of EPIF and tilt stack peak amplitude profile. It can be seen that the fused image highlights the event axis and weakens the influence of noise to a certain extent. Figure 8 Dip profile picked for Radon transform. Fig. 9 is the response of the B-COSFIRE filter. Fig.10 The results of event picking using the method of the present invention show that the event picking of the shallow layer is relatively complete, while the event picking of the deep layer is relatively large and the information is submerged in the noise, so fewer events are picked up. However, the event picking of the thick layer has been completely performed, which illustrates the effectiveness of the present invention.

[0162] The present invention uses a two-dimensional expression of the DoG model to simulate the receptive field model, and the model can realize the function of an excitatory (i.e. positive) central region and an inhibitory surround (i.e. negative) inhibitory domain. The DoG response is blurred and shifted, and its weighted geometric average is obtained to obtain the B-COSFIRE filter response model. Finally, the differential evolution method is used to optimize the best feature selection method to determine the minimum feature subset that maximizes the event depiction performance. The present invention constructs a seismic data event picking method based on the analysis of the information processing mechanism of the primary visual cortex of biological vision. The model proposed by the present invention can effectively detect the event edge of the target contour of seismic data with huge data volume and numerous noise interference. Compared with other traditional edge detection, the B-COSFIRE filter proposed by the present invention can pick up the events of seismic data of different widths, which is more in line with the biological vision mechanism. The events picked up by this method are more continuous and accurate than those picked up by the existing method, providing a new research direction for future seismic data interpretation.

Claims

1. A seismic event extraction method based on the B-COSFIRE filter of the differential evolution algorithm, characterized in that: The following steps are involved: Step 1): Acquire seismic data and preprocess the seismic data; Step 2): Use the two-dimensional expression of the DoG model to simulate the receptive field model, i.e., the B-COSFIRE filter, wherein the DoG model can realize an inhibitory field with an excitatory central region and an inhibitory surround; The two-dimensional expression of the DoG model is as follows: The corresponding value of the central area with excitability is set to 0.5σ, where σ is the spatial distribution of the Gaussian function of the surrounding inhibitory domain. It is an On-centered receptive field, that is, the center is the excitation field and the periphery is the inhibition field. It is the Off center receptive field, defined as: For a given location (x, y) and an intensity distribution I(x', y') of seismic data I, the response c of the DoG is σ (x,y) is calculated as follows: Among them, * is convolution, |·| + It is half-wave rectification; Step 3): setting two groups of B-COSFIRE filters, combining the responses of the two groups of B-COSFIRE filters together for event picking of seismic data, wherein one group of the two groups of B-COSFIRE filters is used for detecting the event, and the other group is used for detecting the end of the event; Step 4): Use the preprocessed seismic data as the input of the B-COSFIRE filter, blur and shift the DoG response, and obtain the weighted geometric mean to obtain the output of the B-COSFIRE filter; The DoG response is blurred and shifted, and its weighted geometric average is obtained to obtain the output of the B-COSFIRE filter, which is specifically: With a length of ρ i The vector moves each ambiguous DoG response toward the center of the support region, and the associated change vector is (Δx i ,Δy i ),Δx i =-ρ i cosφ i , Δy i =-ρ i sinφ i , in φ i Calculate the tuple (σ in the reverse direction i ,ρ i ,φ i ) and the response of DoG and shift it as shown in the following formula: In the formula, σ i is the standard deviation of the DoG model, (ρ i ,φ i ) is the polar coordinate, G σ' (x', y') is the maximum weighted response of DoG, -3σ'≤x', y'≤3σ', where σ'=σ'0+αρ i , the values ​​of σ'0 and α are constants; The weighted geometric mean is calculated as follows: In the formula, S is the output response of the B-COSFIRE filter. s The weighted geometric mean is an AND-type function, that is, a B-COSFIRE filter is only valid if all incoming blur and shift responses Response can only be achieved when it is greater than zero; Step 5): Based on the output results of step 4), the differential evolution algorithm is used as a feature selection method to evaluate the impact on system performance, determine the minimum feature subset that maximizes the event depiction performance, and complete the seismic event extraction.

2. The seismic event extraction method of B-COSFIRE filter based on differential evolution algorithm according to claim 1, characterized in that: The seismic data is in sgy format.

3. The seismic event extraction method of B-COSFIRE filter based on differential evolution algorithm according to claim 1, characterized in that: The preprocessing is specifically: wavelet fusion of the envelope peak instantaneous frequency of pre-stack seismic data and the inclined stack peak instantaneous amplitude.

4. The seismic event extraction method of B-COSFIRE filter based on differential evolution algorithm according to claim 1, characterized in that: The differential evolution algorithm performs iterative search using mutation operations, crossover operations and selection operations.

5. The seismic event extraction method of B-COSFIRE filter based on differential evolution algorithm according to claim 4, characterized in that: In the mutation operation, the local fitness function is used to guide the mutation direction of each sub-component.

6. The seismic event extraction method of B-COSFIRE filter based on differential evolution algorithm according to claim 4, characterized in that: In the crossover operation, a uniform crossover method is adopted.

Images