A high-resolution frequency-wavenumber method for modal separation

Abstract

This paper provides a high-resolution frequency-wavenumber method capable of modal separation. Based on the frequency, wavenumber, and phase velocity characteristics of multimodal Rayleigh surface waves, a method for enhancing high-modal energy by weighting the F-K spectrum with a two-dimensional Gaussian window is proposed. This algorithmic approach enhances high-modal Rayleigh wave information, resulting in simple computation and clear physical meaning. Using probability distribution statistics, a frequency-phase velocity probability distribution is plotted. All extreme points are statistically analyzed, and the median of their phase velocity intervals is taken as the phase velocity of the frequency point in that mode. This more accurately separates and extracts multimodal dispersion curves.

Description

Technical Field

[0001] The present invention relates to the field of seismic exploration, and in particular to a method for extracting phase velocity from micro-motion or seismic background noise imaging. The technology proposes a new method for enhancing and separating high-mode energy of Rayleigh surface waves to extract multi-modal Rayleigh surface wave dispersion curves. Background Art

[0002] The high-resolution frequency-wavenumber method is an algorithm for extracting surface wave phase velocity from microtremor or seismic background noise. It obtains information about the dominant Rayleigh wave dispersion from an energy perspective by calculating the FK spectrum. The core of the high-resolution frequency-wavenumber method is to solve the FK spectrum. Estimating the FK spectrum essentially means estimating the time delay t between stations. Assuming that the wavefront passes through the array as a straight line without any curvature, the time delay t can be calculated from the relative positions of the detectors. When the wavenumber vector k is unknown, a grid search can be performed on all possible wavenumber vectors. For each wavenumber in the grid, the current wavenumber is considered the dominant wavenumber and its beam power is calculated. The maximum energy in the beam power generally corresponds to the wavenumber vector k of the incident dominant Rayleigh wave. Combining the data collected by each detector, the time delay t can be calculated by cross-correlation. Since the relative positions of the detectors determine the array response function, the FK spectrum can be derived from the cross-spectral matrix and the array response function.

[0003] The existing HRFK method extracts the dominant mode Rayleigh wave information with the maximum energy, so what is extracted is often the fundamental wave dispersion curve. However, the extreme points of its FK spectrum contain the dispersion information of higher modes that are suppressed and difficult to extract.

[0004] In the research on extracting high-modal Rayleigh surface wave dispersion curves, for active source surface wave exploration with equally spaced measuring points and a fixed wave propagation direction, the resolution of each mode of the traditional FK algorithm can be improved by decomposing the wavefield response into a phase spectrum and an amplitude spectrum. However, this method is not suitable for exploration based on micro-motion background noise. Another processing method is to enhance the high modes on the final frequency-phase velocity image. For example, a weighted sliding average window is used to process the amplitude spectrum in the frequency-phase velocity domain, or the lambda power of the amplitude value of each point on the dispersion energy spectrum is calculated. However, this method does not enhance the high modes in principle. Summary of the Invention

[0005] The object of the present invention is to provide a high-resolution frequency wavenumber method capable of mode separation.

[0006] The object of the present invention is achieved in that a high-resolution frequency-wavenumber method with modal separation comprises the following steps:

[0007] Step 1: Calculation of FK spectrum

[0008] The FK spectrum is obtained by convolving the cross-spectral matrix and the array response function.

[0009] The cross-spectral matrix reflects the correlation of data collected by different stations. The phase difference is analyzed based on the cross-correlation in the frequency domain, which is expressed as a time delay t in the time domain. The formula for calculating the cross-spectral matrix Φ is:

[0010]

[0011] In the formula: *---conjugate operation;

[0012] F i ---Fourier transform of the current sampling point of station i.

[0013] Beamforming technology is the core of array signal processing. Its essence is to make the acquisition array gain constant in the direction of the desired signal through certain weighting, thus achieving the purpose of spatial filtering. This spatial filter is the array response function, which is expressed as

[0014]

[0015] Where: X n ---The coordinate position of station n.

[0016] The high-resolution frequency-wavenumber method estimates the FK spectrum as

[0017]

[0018] Where: w---angular frequency;

[0019] φ nm ---The values ​​of station n and station m on the cross-spectral matrix Φ;

[0020] k x ---The component of the wave number vector k in the x direction.

[0021] According to the above process, the original FK spectrum P solved by the high-resolution frequency-wavenumber method is obtained.

[0022] Step 2: High-modality enhancement

[0023] A two-dimensional Gaussian window is introduced to weight the FK spectrum, thereby achieving a relative enhancement of the energy extreme points of the high-mode Rayleigh surface waves in the central area of ​​the FK spectrum. The expression of the two-dimensional Gaussian window is:

[0024]

[0025] Where: k max ---The maximum wave number searched by the high-resolution frequency-wave number method, half the width of the Gaussian window;

[0026] a---High modal enhancement parameter, which determines the range and effect of central wave number enhancement.

[0027] The enhancement of high-mode by Gaussian window is a direct weighting of the FK spectrum calculated by formula (8). The weighted power spectrum density P′ G for:

[0028] P′ G (k x , k y , w) = P′(k x , k y ,w)G(k x , k y )

[0029] Step 3: Search for extreme points and calculate phase velocity

[0030] Search and summarize all extreme points on the FK spectrum and calculate their corresponding phase velocities. Based on the calculated FK spectrum, the maximum extreme point found is the wave number vector k of the dominant surface wave at the current frequency. For the wave number vector k, its wave number magnitude |k| and phase angle θ are:

[0031]

[0032] When the angular frequency w is constant, the phase velocity V of the dominant surface wave at the current frequency point can be obtained according to the following conversion formula between wave number and phase velocity: R .

[0033]

[0034] Step 4: Probability distribution statistics

[0035] All extreme points extracted from the FK spectra are weighted and analyzed to create a frequency-phase velocity probability distribution diagram. The probability distribution statistics process is divided into two steps: single-frequency phase velocity weight distribution statistics and full-band probability distribution statistics.

[0036] Single frequency phase velocity weight distribution statistics: Summarize the multi-window extreme points under a single frequency point, according to V R (f) = 2πf / |k|, which is converted into the corresponding phase velocity point. The number of phase velocity points at a specific phase velocity value is regarded as the phase velocity weight. The maximum weight at this frequency point is used for normalization to obtain the phase velocity weight distribution at a single frequency point.

[0037] Full-band probability distribution statistics: Summarize the weight distribution of all frequency points, linearly map the weights from 1 to 0 to the color system, and draw a two-dimensional probability distribution graph, where each column represents the phase velocity weight distribution of a frequency point.

[0038] Step 5: Picking the dispersion curve

[0039] The dispersion curve is picked from the frequency-phase velocity probability distribution diagram obtained in the probability distribution step. For the multi-modal velocity band in the probability distribution diagram, a phase velocity interval is divided for each mode, and the median value is taken as the Rayleigh surface wave phase velocity value of that mode.

[0040] The advantages of the present invention are: 1. Based on the frequency, wavenumber, and phase velocity characteristics of multimodal Rayleigh surface waves, a processing method is proposed that uses a two-dimensional Gaussian window to weight the FK spectrum, thereby achieving high-modal energy enhancement. This algorithm enhances high-modal Rayleigh wave information based on algorithmic principles, is computationally simple, and has clear physical meaning. 2. A method using probability distribution statistics is proposed to plot a frequency-phase velocity probability distribution diagram, perform statistical analysis on all extreme points, and take the median of their phase velocity intervals as the phase velocity of the frequency point in that mode, thereby more accurately achieving the separation and extraction of multimodal dispersion curves. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] In order to more clearly illustrate the implementation methods of this specification or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the implementation methods or the description of the prior art. Obviously, the drawings described below are only some implementation methods recorded in this specification. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative labor.

[0042] Figure 1 This is the implementation flow chart of this patent for high modal enhancement and multi-modal separation.

[0043] Figure 2 This is a schematic diagram of the high-modal enhancement method used in this patent. Figure (a) shows the original FK spectrum, obtained by convolving the cross-spectral matrix with the array response function; Figure (b) shows a two-dimensional Gaussian window weighted according to the size and range of the FK spectrum; and Figure (c) shows the FK spectrum after high-modal enhancement.

[0044] Figure 3 This is a comparison of the probability distribution before and after the high-modal enhancement of this patent. Figure (a) is the original frequency-phase velocity probability distribution diagram, where the black line is the picked fundamental wave dispersion curve; Figure (b) is the frequency-phase velocity probability distribution diagram after enhancement, where the black line is the picked fundamental wave dispersion curve, and the gray line is the picked first-order wave dispersion curve; Figure (c) is the frequency-phase velocity probability distribution diagram of the original first-order wave region; Figure (d) is the frequency-phase velocity probability distribution diagram of the first-order wave region after enhancement.

[0045] Figure 4This is a comparison chart of the inversion results of this patent and the drilling results. The left side shows the shear wave velocity layer structure obtained by inverting only the fundamental wave dispersion curve, and the right side shows the shear wave velocity layer structure obtained by inverting the fundamental wave and first-order wave extracted by this method. The composition of the media in each formation is as follows:

[0046] Stratum 1: Mainly dark brown and brown silty mudstone, with small amounts of muddy siltstone, siltstone, fine sandstone and gravelly sandstone.

[0047] Stratum 2: dark gray-brown and gray-brown silty mudstone, mainly mudstone, with a small amount of gray-brown and light gray silty mudstone.

[0048] Stratum 3: Dark gray-brown mudstone, primarily silty mudstone, interbedded with thin layers of brown-gray calcareous sandstone. A few thin layers or bands of brown-gray or gray-brown silty marlstone are present in the upper portion. Occasional secondary gypsum is observed.

[0049] Stratum 4: dark gray-brown mudstone, mainly silty mudstone, interspersed with brown-gray, light brown calcareous siltstone and carbonated siltstone bands or thin layers. DETAILED DESCRIPTION

[0050] The following will be combined with the drawings in the embodiments of this specification to clearly and completely describe the technical solutions in the embodiments of this specification. Obviously, the embodiments described are only part of the embodiments of this specification, not all of the embodiments. Based on the embodiments of this specification, all other embodiments obtained by ordinary technicians in this field without making creative efforts should fall within the scope of protection of this application.

[0051] A high-resolution frequency-wavenumber method for mode separation comprising the following steps:

[0052] Step 1: Calculation of FK spectrum

[0053] The FK spectrum is obtained by convolving the cross-spectral matrix and the array response function.

[0054] The calculation formula of the cross-spectral matrix Φ is:

[0055]

[0056] In the formula: *---conjugate operation;

[0057] F i ---Fourier transform of the current sampling point of station i.

[0058] The expression of the array response function is:

[0059]

[0060] Where: X n---The coordinate position of station n.

[0061] The high-resolution frequency-wavenumber method estimates the FK spectrum as

[0062]

[0063] Where: w---angular frequency;

[0064] φ nm ---The values ​​of station n and station m on the cross-spectral matrix Φ;

[0065] k x ---The component of the wave number vector k in the x direction.

[0066] According to the above process, the original FK spectrum P solved by the high-resolution frequency-wavenumber method is calculated.

[0067] Step 2: High-modality enhancement

[0068] A two-dimensional Gaussian window is introduced to weight the FK spectrum, and the energy extreme points of the high-mode Rayleigh surface waves in the central area of ​​the FK spectrum are relatively enhanced. The expression of the two-dimensional Gaussian window is:

[0069]

[0070] Where: k max ---The maximum wave number searched by the HRFK method, half the width of the Gaussian window;

[0071] a---High modal enhancement parameter, which determines the range and effect of central wave number enhancement.

[0072] According to the maximum wave number range observable by the array, a reasonable high-mode enhancement parameter a is selected, and the original FK spectrum is directly weighted. The weighted power spectrum density P′ G for:

[0073] P′ G (k x , k y , w) = P′(k x , k y ,w)G(k x , k y ) Step 3: Search for extreme points and calculate phase velocity

[0074] Search and summarize all extreme points on the FK spectrum and calculate their corresponding phase velocities. Based on the calculated FK spectrum, the maximum extreme point found is the wave number vector k of the dominant surface wave at the current frequency. For the wave number vector k, its wave number magnitude |k| and phase angle θ are:

[0075]

[0076] When the angular frequency w is constant, the phase velocity V of the dominant surface wave at the current frequency point can be obtained according to the following conversion formula between wave number and phase velocity: R .

[0077]

[0078] Step 4: Probability distribution statistics

[0079] All extreme points extracted from the FK spectra are weighted and analyzed to create a frequency-phase velocity probability distribution diagram. The probability distribution statistics process is divided into two steps: single-frequency phase velocity weight distribution statistics and full-band probability distribution statistics.

[0080] (1) Single frequency phase velocity weight distribution statistics: Summarize the multi-window extreme points under a single frequency point, according to V R (f) = 2πf / |k|, which is converted into the corresponding phase velocity point. The number of phase velocity points at a specific phase velocity value is regarded as the phase velocity weight. The maximum weight at this frequency point is used for normalization to obtain the phase velocity weight distribution at a single frequency point.

[0081] (2) Full-band probability distribution statistics: Summarize the weight distribution of all frequency points, linearly map the weights from 1 to 0 to the color system, and draw a two-dimensional probability distribution diagram, where each column represents the phase velocity weight distribution of a frequency point.

[0082] Step 5: Pick the dispersion curve.

[0083] The dispersion curve is picked from the frequency-phase velocity probability distribution diagram obtained from the probability distribution statistics link. For the multi-modal velocity band in the probability distribution diagram, the phase velocity interval is divided for each mode, and the dispersion curve is picked using the median fitting method: at a specific frequency point, the start and end intervals of the given phase velocity are given, and the phase velocity weights within the range are summed, recorded as W sum , the weights are accumulated from the low phase velocity point to the high phase velocity point, and the accumulated result is recorded as W i , when W i ≥W sum When the corresponding phase velocity point V Ri That is the phase velocity fitting value of the frequency point in this velocity range.

[0084] Step 6: Joint inversion of multi-modal Rayleigh wave dispersion curves

[0085] Based on the extracted multi-modal Rayleigh surface wave dispersion curves, a joint inversion is performed to obtain the shear wave velocity layered structure, deduce the layered information of the underground medium structure, and realize the exploration of the underground medium.

[0086] The present invention is based on the formation and propagation principles of multimodal Rayleigh surface waves, combined with the distribution characteristics of the phase velocity and wave number of multimodal Rayleigh surface waves on the FK spectrum at specific frequencies, and provides a processing method for achieving high-modal Rayleigh surface wave energy enhancement through two-dimensional Gaussian window weighting. The probability distribution statistics of the processed and extracted multi-extreme points are performed, and a frequency-phase velocity probability distribution diagram is drawn, thereby realizing the separation and extraction of the multimodal Rayleigh surface wave dispersion curve.

[0087] Next, we combine actual data collection with the method described in this case to separate and extract multimodal dispersion curves and invert the layered structure of the underground medium.

[0088] Based on the background noise method, the micro-seismic detector is used to collect micro-seismic data at a single measuring point and multiple stations. The high-resolution frequency wave number method with modal separation proposed in this case is used to process and calculate the collected data at a single measuring point, and invert the accurate and reliable underground medium layered structure. The process of this patented technology is shown in Figure 1 , the high modal enhancement diagram is shown in Figure 2 , the probability distribution comparison before and after high modal enhancement is shown in Figure 3 , the comparison between the inversion effect and the drilling result is shown in Figure 4 .

[0089] Combining the parameters in the attached figure with the calculation method of this case, the dispersion curves of the multi-mode Rayleigh surface waves are separated. The specific implementation steps of the project are as follows:

[0090] 1. Deploy a geophone array to collect micromotion data.

[0091] 2. Input the collected micro-motion data into the calculation program and perform time domain processing (counting bad points, frequency points and frequency windowing).

[0092] 3. According to the array layout, set the maximum wave number vector k to search max In this process, the wave number range k max It determines the range of the maximum search wave number and the width of the Gaussian window. The k value of the micro-motion data is calculated based on the observation array response function. max is 0.033;

[0093] The original FK spectrum is calculated according to the following formula:

[0094]

[0095] Where: w---angular frequency;

[0096] φ nm ---The values ​​of station n and station m on the cross-spectral matrix Φ;

[0097] k x---The component of the wave number vector k in the x direction.

[0098] 4. Based on the estimation of the initial formation conditions and the maximum wave number vector k max , set the high-modal enhancement parameter a. Parameter a determines the range and effect of the center wavenumber enhancement. For the micro-motion data actually implemented this time, in order to control the enhanced wavenumber within a reasonable range, a is set to 0.08. The two-dimensional Gaussian window used in this processing is calculated according to the following formula:

[0099]

[0100] The original FK spectrum is weighted using a two-dimensional Gaussian window to obtain the FK spectrum after high modality enhancement as follows:

[0101] P′ G (k x , k y , w) = P′(k x , k y ,w)G(k x , k y )

[0102] 5. Search for extreme points on the FK spectrum after high modal enhancement and convert the extreme points into phase velocity points according to the following formula:

[0103] 6. Search and summarize all extreme points on the FK spectrum and calculate their corresponding phase velocities. Based on the calculated FK spectrum, the maximum extreme point found is the wave number vector k of the dominant surface wave at the current frequency. For the wave number vector k, its wave number magnitude |k| and phase angle θ are:

[0104]

[0105] When the angular frequency w is constant, the phase velocity V of the dominant surface wave at the current frequency point can be obtained according to the following conversion formula between wave number and phase velocity: R .

[0106]

[0107] 7. Perform weighted analysis on all extreme points extracted from all FK spectra and plot the frequency-phase velocity probability distribution. The probability distribution statistics process is divided into two steps: single-frequency phase velocity weighted distribution statistics and full-band probability distribution statistics.

[0108] (1) Single frequency phase velocity weight distribution statistics: Summarize the multi-window extreme points under a single frequency point, according to V R(f) = 2πf / |k|, which is converted into the corresponding phase velocity point. The number of phase velocity points at a specific phase velocity value is regarded as the phase velocity weight. The maximum weight at this frequency point is used for normalization to obtain the phase velocity weight distribution at a single frequency point.

[0109] (2) Full-band probability distribution statistics: Summarize the weight distribution of all frequency points, linearly map the weights from 1 to 0 to the color system, and draw a two-dimensional probability distribution diagram, where each column represents the phase velocity weight distribution of a frequency point.

[0110] 8. Pick up the dispersion curve from the frequency-phase velocity probability distribution diagram obtained from the probability distribution statistics link. For the multi-modal velocity band in the probability distribution diagram, divide the phase velocity interval for each mode and pick up the dispersion curve using the median fitting method: at a specific frequency point, given the start and end intervals of the phase velocity, sum the phase velocity weights within the range, denoted as W sum , the weights are accumulated from the low phase velocity point to the high phase velocity point, and the accumulated result is recorded as W i , when W i ≥W sum When the corresponding phase velocity point V Ri That is the phase velocity fitting value of the frequency point in this velocity range.

[0111] 9. Based on the extracted multi-modal Rayleigh surface wave dispersion curves, a joint inversion is performed to obtain the shear wave velocity layered structure, deduce the layered information of the underground medium structure, and realize the exploration of the underground medium.

[0112] Although the present application is described above through embodiments, those skilled in the art will appreciate that the present application may be subject to many modifications and variations without departing from the spirit of the present application, and the appended claims include these modifications and variations without departing from the inventive content of the present application.

Claims

1. A high-resolution frequency-wavenumber method with modal separation, characterized in that: It includes the following steps: Step 1: Calculation of FK spectrum The FK spectrum is obtained by convolving the cross-spectral matrix and the array response function; The cross-spectral matrix reflects the correlation of data collected by different stations. The phase difference is analyzed based on the cross-correlation in the frequency domain, which is expressed as a time delay t in the time domain. The cross-spectral matrix Φ is calculated as follows: In the formula: *---conjugate operation; F i ---Fourier transform of the current sampling point of station i; Beamforming technology is the core of array signal processing. Its essence is to achieve spatial filtering by using certain weights to make the acquisition array gain constant in the direction of the desired signal. This spatial filter is the array response function, which is expressed as follows: Where: X n ---The coordinate position of station n; The high-resolution frequency-wavenumber method estimates the FK spectrum as follows: Where: ω---angular frequency; Φ nm ---The values ​​of station n and station m on the cross-spectral matrix Φ; k x ---The component of the wave number vector k in the x direction; According to the above process, the original FK spectrum P solved by the high-resolution frequency-wavenumber method is obtained; Step 2: High-modality enhancement A two-dimensional Gaussian window is introduced to weight the FK spectrum, thereby achieving relative enhancement of the energy extreme points of the high-mode Rayleigh surface waves in the central area of ​​the FK spectrum. The expression of the two-dimensional Gaussian window is: Where: k max ---The maximum wave number searched by the high-resolution frequency-wave number method, half the width of the Gaussian window; a---High mode enhancement parameter, which determines the range and effect of central wave number enhancement; The enhancement of the high mode by the Gaussian window is a direct weighting of the FK spectrum estimation. The weighted power spectrum density P′ G for: P′ G (k x ,k y ,ω)=P′(k x ,k y ,ω)G(k x ,k y ) Step 3: Search for extreme points and calculate phase velocity All extreme points on the FK spectrum are searched and summarized, and their corresponding phase velocities are calculated. Based on the calculated FK spectrum, the maximum extreme point found is the wave number vector k of the dominant surface wave at the current frequency point. For the wave number vector k, its wave number magnitude |k| and phase angle θ are: When the angular frequency ω is constant, the phase velocity V of the dominant surface wave at the current frequency point can be obtained according to the following conversion formula between wave number and phase velocity: R ; Step 4: Probability distribution statistics Perform weight analysis on all extreme points extracted from all FK spectra and draw a frequency-phase velocity probability distribution diagram; Step 5: Picking the dispersion curve The dispersion curve is picked from the frequency-phase velocity probability distribution diagram obtained from the probability distribution link; for the multi-modal velocity band in the probability distribution diagram, the phase velocity interval is divided for each mode, and the median value is taken as the Rayleigh surface wave phase velocity value of the mode.

2. The high-resolution frequency-wavenumber method with modal separation according to claim 1, characterized in that: The specific process of the probability distribution statistics in step 4 includes dividing the probability distribution statistics process into two steps: single frequency point phase velocity weight distribution statistics and full frequency band probability distribution statistics; the single frequency point phase velocity weight distribution statistics method includes: summarizing the multi-window extreme value points under the single frequency point, according to V R (f)=2πf / |k|, convert it into the corresponding phase velocity point, regard the number of phase velocity points at a specific phase velocity value as the phase velocity weight, and use the maximum weight at the frequency point for normalization to obtain the phase velocity weight distribution at a single frequency point; the method for full-band probability distribution statistics includes: summarizing the weight distribution of all frequency points, linearly mapping the weight from 1 to 0 to a color system, and drawing a two-dimensional probability distribution diagram, in which each column represents the phase velocity weight distribution of a frequency point.

3. The high-resolution frequency-wavenumber method with modal separation according to claim 1, characterized in that: The step 5 of picking the dispersion curve includes the following steps: For the multi-modal velocity band in the probability distribution diagram, the phase velocity interval is divided for each mode, and the dispersion curve is picked up by the median fitting method: at a specific frequency point, the start and end intervals of the phase velocity are given, and the phase velocity weights within the range are summed, which is recorded as W sum , the weights are accumulated from the low phase velocity point to the high phase velocity point, and the accumulated result is recorded as W i , when W i ≥W sum When the corresponding phase velocity point V Ri That is the phase velocity fitting value of the frequency point in this velocity range.

Images