An Automatic Extraction Method for Dispersion Curves of Micro-Moving Surface Waves Based on Connected Component Analysis

Abstract

This invention relates to the field of micro-motion surface wave exploration methods, and discloses an automatic extraction method for dispersion curves of micro-motion surface waves based on connected component analysis. The method includes cross-correlation extraction of virtual source surface wave shot gathers, frequency-wavenumber domain amplitude spectrum calculation, adaptive threshold binarization, connected component analysis, frequency-wavenumber domain surface wave mask generation, mask-amplitude matrix element multiplication, phase velocity picking, Hampel filtering, and automatic generation of dispersion curves. This invention utilizes image connected component analysis and other operations to process the amplitude spectrum of micro-motion surface waves, separating the surface wave signal region from the frequency-wavenumber domain amplitude spectrum image. It automatically identifies discrete phase velocity points from images containing only surface wave dispersion energy stripes and uses Hampel filtering to remove distortion points, ultimately achieving high-precision automated measurement of the dispersion curve.

Description

Technical Field

[0001] This invention relates to the field of micro-motion surface wave exploration methods, and more particularly to the processing and analysis of micro-motion surface wave data. Specifically, it relates to an automatic extraction method for micro-motion surface wave dispersion curves based on connected component analysis. Background Technology

[0002] Surface waves are seismic wave fields that propagate along the Earth's free surface. During propagation, different frequency components of the surface wave propagate at different speeds, exhibiting dispersion characteristics. Both theory and practice have shown that the dispersion characteristics of surface waves are controlled by the velocity structure of the medium beneath their propagation path. By utilizing the dispersion characteristics of surface waves, the internal structure of the medium can be obtained. Therefore, surface wave detection methods are widely used in various fields such as engineering exploration, resource exploration, and geological imaging.

[0003] Micromotion is an effective method for seismic surface wave exploration. Micromotions refer to continuous, minute vibrations on Earth that can be observed at any time and location. The sources of micromotion signals are very diverse, including ground vibrations caused by industrial production and transportation in human life, as well as surface vibrations caused by atmospheric circulation and ocean currents in nature. By collecting micromotion signals and extracting surface wave information from them, and then combining this with surface wave exploration technology, it is possible to detect underground structures.

[0004] The key to using micromotion surface waves to detect subsurface structures is obtaining reliable surface wave dispersion curves, which reflect the velocity and thickness information of various subsurface strata. Dispersion curve extraction typically employs the SPAC (spatial autocorrelation) method and the FK (frequency-wavenumber) method. The SPAC method requires the seismograph acquiring micromotion signals to be deployed in a multi-ring array, which is difficult to implement in the field and results in low resolution; its results reflect the combined values ​​of strata parameters within the ring area. The arrangement of seismograph arrays along direct seismic lines is relatively convenient. This patent mainly optimizes the micro-motion surface wave exploration method acquired through direct seismic lines. First, the single-shot surface wave record of the virtual source is recovered from the micro-motion signal acquired through direct seismic lines using cross-correlation technology. Then, the dispersion curve is extracted based on the FK method. However, due to the characteristics of micro-motion signals, the surface wave signal obtained from micro-motion often has a low signal-to-noise ratio. Therefore, the dispersion energy spectrum generated by the conventional FK method has poor resolution and is not easy to identify the dispersion curve. Moreover, when the data volume is large or the exploration area is large, there can be hundreds or even thousands of virtual source surface wave records. If the dispersion curve is still extracted from each set of traces by manual identification, not only is the identification efficiency low, but it is also prone to misjudgment. Therefore, how to identify the surface wave signal area from the dispersion energy spectrum image and suppress the energy of other noise and interfering wave fields in the image to obtain a simple and single surface wave dispersion energy strip is the key to achieving high-precision automatic extraction of micro-motion surface wave dispersion curves. Summary of the Invention

[0005] To address the aforementioned shortcomings of existing technologies, the present invention aims to provide an automatic extraction method for micro-motion surface wave dispersion curves based on connected component analysis.

[0006] To achieve the above objectives, the present invention adopts the following technical solution:

[0007] The present invention discloses an automatic extraction method for dispersion curves of micro-moving surface waves based on connected component analysis, comprising the following steps:

[0008] Step 1: Based on cross-correlation technology, extract the virtual source surface wave single-shot record gather from the multi-channel micro-motion records of the straight measurement line. The virtual source position is the position of the first detector at the end of the measurement line. If the original multi-channel micro-motion record is denoted as U(x,t), then the reconstructed virtual source surface wave gather record is U(x,t). s (x,t), where x represents the distance between the detector and the virtual source, and t represents the recording time of each detector;

[0009] Step 2: Record U for the virtual source surface wave gather s Perform an FK transform on (x,t) and plot the result in the frequency-wavenumber domain to generate an energy spectrum image E. s (f,k), where f represents the frequency and k represents the wave number;

[0010] Step 3: Analyze the energy spectrum image E in the frequency wavenumber domain. s Binarize (f,k) to obtain a binary image consisting only of black or white, denoted as B_E. s (f,k), the threshold for binarization is adaptively determined by the maximum inter-class variance method.

[0011] Step 4: Process the binary image B_E s The image (f,k) is processed by connected component analysis. Based on the connection relationship between each pixel and its neighboring pixels, the binary image B_E is transformed. s Find each connected region in (f,k) and label them E1, E2, E3, ... according to their area.

[0012] Step 5: Convert image B_E s The two largest connected regions E1 and E2 in (f,k) are retained, and the pixel values ​​in the regions E1 and E2 are set to 1, while the pixel values ​​in other regions of the image are set to 0. The resulting new image is denoted as M, which contains only two values, 0 and 1, and can be regarded as a mask.

[0013] Step Six: Compare the mask M with the frequency-wavenumber domain energy spectrum E obtained in Step Two. s Let (f,k) be two matrices with the same number of rows and columns. Multiply the elements at the same positions in the two matrices, and denote the result as a new matrix E. d(f,k), since M contains only two element values, 0 and 1, then in matrix E d The non-zero portion of the image formed by (f,k) is the surface wave energy stripe;

[0014] Step 7: Apply the formula v = f / k to E d A projection transformation is performed on (f,k), where v represents the phase velocity of the surface wave, f represents the frequency, and k represents the wavenumber, transforming the energy spectrum from the frequency-wavenumber domain to the frequency-velocity domain. The result is denoted as E. d (f,v), through E d Find the energy distribution of (f,v) at each frequency f. i The maximum energy point in the cross section is denoted as P(f). i ,v i ), P(f i ,v i ) is a two-dimensional array, with the first column being the frequency [f i The second column is [v] i ].

[0015] Step 8: Use the Hampel filter to adjust P(f) i ,v i The second column vector [v] in ) i Median filtering is performed to remove outlier values ​​from the data. The filtered column vector is denoted as [phv]. i Then use the moving average to analyze phv i Perform smoothing processing, and the result after processing is denoted as [phase_v] i ].

[0016] Step Nine: Place [f] i ] and [phase_v i Energy spectrum E plotted in the frequency-velocity domain d In the (f,v) image, where [f i The element f in ] corresponds to the frequency axis, [phase_v i The element phase_v in the diagram corresponds to the velocity axis. i That is, frequency f i The phase velocity value at a given point is obtained by connecting all the phase velocity points in the graph to obtain the dispersion curve.

[0017] Preferably, in step four, when performing connected component analysis, an 8-connectivity method is used to determine the connection relationship between a pixel and its surrounding neighboring pixels.

[0018] Furthermore, in step seven, E d (f,v) energy spectrum f i The frequency cross section can be considered as having a frequency component f. iThe amplitude curve of the signal energy as a function of velocity is given by the frequency f, where the velocity corresponding to the maximum value of the amplitude curve is the frequency f. i Phase velocity v at time i However, this speed contains errors and still needs to be processed in step eight.

[0019] Furthermore, in step eight, the window size of the Hampel filter processing parameters is set to 3, the standard deviation threshold parameter is set to 2, and the number of iterations is set to 2; the window width of the moving average processing is set to 5.

[0020] Compared with the prior art, the present invention has the following beneficial effects.

[0021] (1) This invention uses image connected component analysis to process the frequency wavenumber spectrum of micro-moving surface waves. It divides the regions where different signals are located in the image in the frequency-wavenumber domain, extracts the frequency wavenumber spectrum of the surface waves based on the energy characteristics of the micro-moving surface waves, and generates a dispersion energy spectrum containing only surface wave energy strips in the frequency-velocity domain through projection transformation, thereby realizing the automatic extraction of phase velocity and dispersion curve.

[0022] (2) The present invention utilizes Hampel filtering technology to further filter the automatically extracted phase velocity points, which can eliminate distorted points in the data and reduce the disturbance error caused by discrete FK transformation, thereby improving the accuracy of automated measurement of dispersion curves. Attached Figure Description

[0023] Other features, objects, and advantages of the present invention will become more apparent from the following detailed description of non-limiting embodiments with reference to the accompanying drawings:

[0024] Figure 1 This is a flowchart of the method of the present invention;

[0025] Figure 2 This is a schematic diagram of the micro-motion surface wave virtual source gather recording in the embodiment;

[0026] Figure 3 This is a schematic diagram of the amplitude spectrum in the frequency-wavenumber domain in the embodiment.

[0027] Figure 4 This is a schematic diagram of the mask in the embodiment;

[0028] Figure 5 This is a schematic diagram of the amplitude spectrum containing only the frequency-wavenumber domain of surface waves in the embodiment.

[0029] Figure 6 This is a schematic diagram of the surface wave dispersion energy spectrum and its automatically generated dispersion curve in the embodiment. Detailed Implementation

[0030] The present invention will now be described in detail with reference to specific embodiments. These embodiments will help those skilled in the art to further understand the present invention, but do not limit the invention in any way. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention. These all fall within the scope of protection of the present invention.

[0031] Example

[0032] like Figure 1 As shown, an automatic extraction method for surface wave dispersion curves in petroleum seismic data includes the following specific steps:

[0033] Step 1: Based on cross-correlation technology, extract virtual source surface wave shot gathers from the original direct-line multi-channel micro-motion records. Figure 2 For the extracted surface wave record U s (x,t), where x represents the distance between the detector and the virtual source, t represents the recording time of each detector, and the virtual source position is the 0th channel position at the end of the measurement line. The signals of channels 1-13 in the figure are the virtual source records collected by each detector. The channel spacing is 5m, the sampling rate of each channel is 500Hz, and the recording length is uniformly set to 1s. The obvious surface wave waveform information can be seen in the 0-0.4s range of the figure.

[0034] Step 2: Record U for the virtual source surface wave gather s Perform an FK transform on (x,t) and plot the result in the frequency-wavenumber domain to generate an energy spectrum image E. s (f,k), where f represents the frequency and k represents the wave number, such as Figure 3 The image shows a frequency-wavenumber spectrum, from which surface energy bands and energy clusters formed by other noise can be observed.

[0035] Step 3: Analyze the energy spectrum image E in the frequency wavenumber domain. s Binarize (f,k) to obtain a binary image consisting only of black or white, denoted as B_E. s (f,k), the threshold for binarization is adaptively determined by the maximum inter-class variance method.

[0036] Step 4: Process the binary image B_E s The image (f,k) is processed by connected component analysis. Based on the connection relationship between each pixel and its neighboring pixels, the binary image B_E is transformed. s Find each connected region in (f,k) and label them E1, E2, E3, etc. according to their area.

[0037] Step 5: Convert image B_E sThe two largest connected components E1 and E2 in (f,k) are preserved, and the pixel values ​​within E1 and E2 are set to 1, while the pixel values ​​in other regions of the image are set to 0. The resulting new image is denoted as M, which contains only 0 and 1 values ​​and can be considered as a mask. Figure 4 The white area in the mask represents the corresponding region of the surface wave signal in the frequency wavenumber spectrum.

[0038] Step Six: Compare the mask M with the frequency-wavenumber domain energy spectrum E obtained in Step Two. s Let (f,k) be two matrices with the same number of rows and columns. Multiply the elements at the same positions in the two matrices, and denote the result as a new matrix E. d (f,k), since M contains only two element values, 0 and 1, then in matrix E d The non-zero portion of the image formed by (f,k) represents the surface wave energy stripes, such as... Figure 5 This is a frequency wavenumber spectrum image containing only surface waves.

[0039] Step 7: Apply the formula v = f / k to E d A projection transformation is performed on (f,k), where v represents the phase velocity of the surface wave, f represents the frequency, and k represents the wavenumber, transforming the energy spectrum from the frequency-wavenumber domain to the frequency-velocity domain. The result is denoted as E. d (f,v), Figure 6 The frequency-velocity domain dispersed energy spectrum contains visible surface wave energy bands; via E d Find the energy distribution of (f,v) at each frequency f. i The maximum energy point in the cross section is denoted as P(f). i ,v i ), P(f i ,v i ) is a two-dimensional array, with the first column being the frequency [f i The second column is [v] i ], Figure 6 The black circles in the diagram represent the P(f) at the maximum amplitude corresponding to each frequency. i ,v i )Location.

[0040] Step 8: Use the Hampel filter to adjust P(f) i ,v i The second column vector [v] in ) i Median filtering is performed to remove outlier values ​​from the data. The filtered column vector is denoted as [phv]. i Then use the moving average to analyze phv i Perform smoothing processing, and the result after processing is denoted as [phase_v] i ].

[0041] Step Nine: Place [f]i ] and [phase_v i Energy spectrum E plotted in the frequency-velocity domain d In the (f,v) image, where [f i The element f in ] corresponds to the frequency axis, [phase_v i The element phase_v in the diagram corresponds to the velocity axis. i That is, frequency f i The phase velocity value at a given point is used to connect all the phase velocity points in the graph to obtain the dispersion curve. Figure 6 The black curve in the image is the dispersion curve automatically generated after the above steps.

[0042] The specific embodiments of the present invention have been described above. It should be understood that the present invention is not limited to the specific embodiments described above, and those skilled in the art can make various modifications or variations within the scope of the claims, which do not affect the essence of the present invention.

Claims

1. An automatic extraction method for dispersion curves of micro-moving surface waves based on connected component analysis, characterized in that, Includes the following steps: Step one: Based on cross-correlation technique, the virtual source plane wave gather is extracted from the direct line multi-channel microseismic record. The virtual source position is the first channel detector position at the end of the measuring line. If the original multi-channel microseismic record is denoted as U(x,t), the recovered virtual source plane wave channel record is denoted as U s (x,t), wherein x represents the distance of the detector from the virtual source, and t represents the time recorded by each detector. Step two: record U of virtual source surface wave gather s (x,t) is F-K transformed and the transformed result is plotted in the frequency-wave number domain to generate an energy spectrum image E s (f,k), where f represents frequency and k represents wave number; Step three: energy spectrum image E in frequency-wave number domain s (f, k) is binarized to obtain a binary image composed of only black or white, denoted as B_E s (f, k), and the threshold of the binarization is adaptively determined by the maximum inter-class variance method. Step four: carry out connected domain analysis processing on the binary image B_E s (f, k) to find each connected domain in the binary image B_E s (f, k) based on the connection relationship between each pixel point and its neighborhood pixels, and sequentially record each connected domain as E1, E2, E3, … according to the size of the area of each domain. Step five: the image B_E s The two largest connected domains E1 and E2 in (f, k) are retained, and the pixel values in the regions of E1 and E2 are set to 1, and the pixel values in other regions of the image are set to 0, to obtain a new image M, which contains only two values 0 and 1 and can be regarded as a mask. Step six: mask M and the frequency-wavenumber energy spectrum E obtained in step two s (f, k) are two matrices with equal number of rows and columns, and the elements in the same position of the two matrices are multiplied to obtain a new matrix E d (f, k), since M only contains two element values of 0 and 1, then in the matrix E d (f, k) constitute the image of the non-zero part of the surface wave energy band; Step 7: Apply the formula v = f / k to E d A projection transformation is performed on (f,k), where v represents the phase velocity of the surface wave, f represents the frequency, and k represents the wavenumber, transforming the energy spectrum from the frequency-wavenumber domain to the frequency-velocity domain. The result is denoted as E. d (f,v), through E d Find the energy distribution of (f,v) at each frequency f. i The maximum energy point in the cross section is denoted as P(f). i ,v i ), P(f i ,v i ) is a two-dimensional array, with the first column being the frequency [f i The second column is [v] i ]; Step 8: Use the Hampel filter to adjust P(f) i ,v i The second column vector [v] in ) i Median filtering is performed to remove outlier values ​​from the data. The filtered column vector is denoted as [phv]. i Then use the moving average to analyze phv i Perform smoothing processing, and the result after processing is denoted as [phase_v] i ]; Step Nine: Place [f] i ] and [phase_v i Energy spectrum E plotted in the frequency-velocity domain d In the (f,v) image, where [f i The element f in ] corresponds to the frequency axis, [phase_v i The element phase_v in the diagram corresponds to the velocity axis. i That is, frequency f i The phase velocity value at a given point is obtained by connecting all the phase velocity points in the graph to obtain the dispersion curve.

2. The automatic extraction method for dispersion curves of micro-moving surface waves based on connected component analysis according to claim 1, characterized in that, In step four, when performing connected component analysis, the 8-connectivity method is used to determine the connection relationship between a pixel and its surrounding neighboring pixels.

3. The automatic extraction method for dispersion curves of micro-moving surface waves based on connected component analysis according to claim 1, characterized in that, In step seven, E d (f,v) energy spectrum f i The frequency cross section can be considered as having a frequency component f. i The amplitude curve of the signal energy as a function of velocity is given by the frequency f, where the velocity corresponding to the maximum value of the amplitude curve is the frequency f. i Phase velocity v at time i However, this speed contains errors and still needs to be processed in step eight.

4. The automatic extraction method for dispersion curves of micro-moving surface waves based on connected component analysis according to claim 1, characterized in that, In step eight, the window size of the Hampel filter processing parameters is set to 3, the standard deviation threshold parameter is set to 2, and the number of iterations is set to 2; the window width of the moving average processing is set to 5.

Images