Coal seam roof sandstone aquifer water yield evaluation method based on three-dimensional earthquake

Abstract

The invention discloses a coal seam roof sandstone aquifer water yield evaluation method based on a three-dimensional earthquake. The method comprises the steps that S10, drilling data, logging data and high-precision three-dimensional earthquake data of a research area are collected; s20, inverting the pre-stack seismic data to obtain a longitudinal-transverse wave velocity ratio of the roof sandstone aquifer; s30, time-frequency analysis is carried out on the post-stack seismic data, and spectrum equalization is carried out on the frequency domain seismic data; s40, according to a reflection coefficient approximation formula and the time-frequency spectrum of the post-stack seismic data after spectrum equalization, inverting the longitudinal wave dispersion attribute of the target layer; s50, using an entropy weight method to calculate the weight values of the longitudinal and transverse wave velocity ratio and the longitudinal wave dispersion attribute; and S60, constructing a water yield indicator factor calculation formula, and evaluating the sandstone aquifer water yield according to the water yield indicator factor. The method is based on the theoretical basis of seismic rock physics, is low in cost and high in operability, has more statistical significance, can quickly and accurately evaluate the water-rich condition of the coal seam roof sandstone aquifer, and provides technical guarantee for development and utilization of coal resources.

Description

Technical Field

[0001] The invention relates to the technical field of seismic rock physics research, and specifically discloses a method for evaluating the water-richness of a coal seam roof sandstone aquifer based on three-dimensional seismic analysis. Background Art

[0002] Mine roof flooding poses a serious threat to coal mine safety and production. As a key source of water for mine flooding, the assessment of the water content of sandstone aquifers in coal seam roofs is crucial for mine flood prevention and control. Currently, electromagnetic exploration, a commonly used method for assessing water content in academia and industry, suffers from shallow detection depths, limited coverage, and insufficient accuracy. As exploration depths increase, these methods are unable to achieve large-scale quantitative or semi-quantitative prediction of roof flooding. Seismic exploration, the most common method for acquiring information about underground media, offers the advantages of wide coverage and high accuracy, providing a new approach for identifying water content.

[0003] The ratio of P-wave to S-wave velocities is an important elastic parameter reflecting the properties of formation fluids. The presence of fluids increases the P-wave velocity of sandstone while having little effect on the S-wave velocity, thereby increasing the P-wave to S-wave velocity ratio. Furthermore, due to the influence of fluids, seismic waves exhibit significant dispersion and attenuation when traveling through sandstone aquifers. This characteristic provides a theoretical foundation for fluid identification and quantitative assessment of water-rich aquifers.

[0004] The entropy weight method is a typical objective weighting method. This method measures the amount of information by calculating the degree of data dispersion of each indicator, and then determines the weight, which can build a more robust and reasonable weight system.

[0005] The water-richness prediction technology based on 3D seismic is applied to the prediction of water-richness of roof sandstone. It can simultaneously consider the P-wave and S-wave velocity ratios and dispersion attribute information of the water-bearing sandstone layer, avoid the limitation of a single attribute being interfered by noise or other factors, and achieve a more reliable water-richness evaluation. Summary of the Invention

[0006] In response to the above-mentioned problems existing in the existing technology, the present invention proposes a method for evaluating the water richness of coal seam roof sandstone aquifers based on three-dimensional seismic analysis. Due to the adoption of the following technical features, the above-mentioned technical objectives can be achieved and many other technical effects can be brought about.

[0007] In order to solve the above problems, the technical solutions adopted by the present invention are as follows:

[0008] A method for evaluating the water richness of a coal seam roof sandstone aquifer based on multi-source data includes the following steps:

[0009] S10: Collect borehole data, well logging data and 3D seismic data in the designated study area, including but not limited to borehole lithologic layering data, density logging data, pre-stack 3D seismic data, and post-stack 3D seismic data;

[0010] S20: Invert pre-stack seismic data to obtain the P-wave and S-wave velocity ratios of the roof sandstone aquifer;

[0011] S30: Perform spectral decomposition on the post-stack seismic data using time-frequency analysis technology to obtain the post-stack seismic time-frequency spectrum, and perform spectral equalization on the time-frequency spectrum;

[0012] S40: Based on the reflection coefficient approximation formula, a P-wave dispersion attribute inversion method is established to invert the frequency domain post-stack seismic data to obtain the P-wave dispersion attributes of the target layer;

[0013] S50: Calculate the entropy weights of the P-wave velocity ratio and P-wave dispersion properties of the water-bearing sandstone layer in the study area according to the entropy weight method;

[0014] S60: Construct a calculation formula for the water-richness indicator factor and evaluate the water-richness of the sandstone aquifer based on the size of the water-richness indicator factor.

[0015] Furthermore, the pre-stack 3D seismic data in step S10 undergoes static correction, dynamic correction, combined deconvolution, and pre-stack denoising; the post-stack seismic data is obtained by stacking and interpolating the pre-stack seismic data, and undergoes random noise attenuation, time-varying filtering, amplitude equalization, and other processing.

[0016] Furthermore, in step S20, inverting the pre-stack seismic data to obtain the P-wave and S-wave velocity ratios of the roof sandstone aquifer comprises the following steps:

[0017] S201: extracting statistical seismic wavelets from pre-stack seismic data, performing well-seismic calibration, and establishing a corresponding relationship between the well logging depth domain and the seismic time domain;

[0018] S202: establishing an initial inversion model of P-wave velocity, S-wave velocity and density based on the calibrated density and velocity logging curves;

[0019] S203: Based on the initial inversion model and pre-stack seismic data, the P-wave and S-wave velocity ratio distribution of the roof sandstone aquifer is obtained.

[0020] Furthermore, in S201, when extracting statistical seismic wavelets from pre-stack seismic data, seismic traces far away from geological structures should be selected, and the wavelet estimation window should be set near the roof sandstone aquifer.

[0021] Furthermore, in step S30, the time-frequency analysis of post-stack seismic data includes the following steps:

[0022] S301: performing smoothed pseudo-Wigner transform (SPWVD) on each seismic record of the post-stack seismic data to obtain amplitude spectrum data;

[0023] S302: Select reference frequency f ref, performing spectrum equalization processing on the amplitude spectrum data obtained in S301;

[0024] S303: Repeat steps S301 and S302 until the time-frequency spectrum S(t,f) of all seismic gathers after spectral equalization is obtained, and store the time-frequency spectrum data in a matrix;

[0025] S304: extracting m spectrum decomposition data after spectrum equalization in the target frequency band from the time-frequency spectrum result of S303.

[0026] Furthermore, in step S301, for a given earthquake record x(t), the calculation formula of the smoothed pseudo-Wigner distribution (SPWVD) is:

[0027]

[0028] Among them, g(u) is the smoothing window function in the time direction, h(τ) is the smoothing window function in the delay direction, x(t) is the earthquake record to be analyzed, and x * (t) is the conjugate complex number of the earthquake record to be analyzed.

[0029] Furthermore, in step S302, the purpose of spectral equalization is to eliminate the influence of wavelet superimposition, and the processing method is:

[0030]

[0031] in, f after spectral equalization m The time spectrum, S(t,f m ) is f m Time spectrum, W(f m ) is the spectrum equalization coefficient, which is calculated by the following formula:

[0032]

[0033] Among them, Max[S(f ref )] is the maximum amplitude value at the reference frequency, Max[S(f m )] is the frequency f m The maximum amplitude value under .

[0034] Furthermore, in step S40, the approximate formula of the reflection coefficient at near-vertical incidence is:

[0035]

[0036] Where R(θ) is the reflection coefficient at the incident angle θ, v p is the average longitudinal wave velocity in the vertical direction of the media on both sides of the interface, v sis the mean shear wave velocity in the vertical direction of the media on both sides of the interface, ρ is the mean density of the media on both sides of the interface, and Δ represents the parameter difference of the media on both sides of the interface;

[0037] Since the post-stack data mainly represents the reflection response of approximately vertical incidence and the incident angle is small, the above formula can be further simplified as:

[0038]

[0039] Since the reflection coefficient is frequency dependent, assuming V p The density is not affected by the frequency change, so:

[0040]

[0041] Select the earthquake main frequency as the reference frequency, and the above formula is ref The first-order Taylor expansion is:

[0042]

[0043] The frequency-dependent reflection coefficient equation at frequency f is compared with the reference frequency f ref Subtracting the frequency-dependent reflection coefficient equation under , we get:

[0044]

[0045] Longitudinal wave dispersion property D p is defined as The above formula can be written in matrix form:

[0046]

[0047] Since seismic records are obtained in actual seismic exploration, the reflection coefficient is converted into a spectrum after spectral equalization, and the above formula can be rewritten as:

[0048]

[0049] The ridge regression algorithm is used to solve the above equation and obtain the dispersion attribute D p :

[0050] D p =2(G T G+λI) -1 G T d (11)

[0051] Among them, G is G T is the transposed matrix of G, λ is the damping parameter, and d is f m represents the mth target frequency.

[0052] Furthermore, in step S50, the entropy weight method is used to calculate the weight values ​​of the P-wave velocity ratio and the P-wave dispersion attribute of the sandstone, including the following steps:

[0053] S501: Standardize each indicator:

[0054]

[0055] Among them, x ij Represents the observation value of the i-th object under the j-th indicator;

[0056] S502: Calculate the proportion of each sample under each indicator:

[0057]

[0058] Where m represents the number of evaluation objects;

[0059] S503: Calculate the information entropy E of each indicator j :

[0060]

[0061] S504: Calculate the weight D of each indicator using the entropy weight method j :

[0062]

[0063] Where n represents the number of indicators.

[0064] Furthermore, in step S60, the water-richness indicator factor WI of the i-th sample point i The calculation formula is:

[0065]

[0066] Among them, D j is the entropy weight value of the jth indicator, I ij is the jth index value of the i-th sample point.

[0067] Compared with the prior art, the present invention has the following beneficial effects:

[0068] 1. The 3D seismic-based water-richness prediction technology is applied to the prediction of water-richness of roof sandstone. It can simultaneously consider the P-wave velocity ratio and S-wave velocity dispersion attribute information of the water-bearing sandstone layer, avoid the limitation of a single attribute being interfered with by noise or other factors, and achieve a more reliable water-richness evaluation.

[0069] 2. The present invention uses the entropy weight method to calculate the data dispersion degree of each indicator to measure its information content, and then determine the weight, which can build a more robust and reasonable weight system.

[0070] 3. This method is low-cost, highly operational, and more statistically significant. It can quickly and accurately evaluate the water-richness of the sandstone aquifer in the coal seam roof, providing technical support for the development and utilization of coal resources. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] Figure 1 It is a schematic flow chart of the method for evaluating the water-richness of a coal seam roof sandstone aquifer based on three-dimensional seismic analysis of the present invention.

[0072] Figure 2 shows 3D seismic data, where Figure 2(a) shows pre-stack 3D seismic data and Figure 2(b) shows post-stack 3D seismic data.

[0073] Figure 3 This is the distribution diagram of the P-wave and S-wave velocity ratio obtained by pre-stack seismic inversion.

[0074] Figure 4 is the spectrum decomposition cross-section diagram of the three-dimensional seismic data after spectrum equalization, where: Figure 4 (a) is the 30HZ spectrum decomposition result; Figure 4 (b) is the 40HZ spectrum decomposition result; Figure 4 (c) is the 50HZ spectrum decomposition result; Figure 4 (d) is the 60HZ spectrum decomposition result.

[0075] Figure 5 This is the plane distribution diagram of the dispersion properties obtained by inversion.

[0076] Figure 6 It is the plane distribution diagram of water-richness indicator factor. DETAILED DESCRIPTION

[0077] In order to make the purpose, technical solution and advantages of the technical solution of the present invention clearer, the technical solution of the embodiment of the present invention will be clearly and completely described below in conjunction with the drawings of the specific embodiment of the present invention. It should be noted that the described embodiment is a part of the embodiment of the present invention, not all the embodiments. Based on the described embodiment of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0078] The specific embodiments of the present invention will be further described below with reference to the accompanying drawings. Figure 1 As shown, the following steps are included:

[0079] S10: Collect drilling data, logging data and 3D seismic data in the designated study area, including but not limited to drilling lithologic layering data, density logging data, pre-stack 3D seismic data, and post-stack 3D seismic data.

[0080] In particular, the collected pre-stack 3D seismic data undergoes static correction, dynamic correction, combined deconvolution, and pre-stack denoising; the post-stack seismic data is obtained by superposition and interpolation of the pre-stack seismic data, and undergoes random noise attenuation, time-varying filtering, amplitude equalization, and other processing.

[0081] S20: Invert pre-stack seismic data to obtain the P-wave and S-wave velocity ratios of the roof sandstone aquifer.

[0082] Specifically, in step S20, the time-frequency analysis of post-stack seismic data includes the following steps:

[0083] S201: extracting statistical seismic wavelets from pre-stack seismic data, performing well-seismic calibration, and establishing a corresponding relationship between the well logging depth domain and the seismic time domain;

[0084] S202: establishing an initial inversion model of P-wave velocity, S-wave velocity and density based on the calibrated density and velocity logging curves;

[0085] S203: Based on the initial inversion model and pre-stack seismic data, the P-wave and S-wave velocity ratio distribution of the roof sandstone aquifer is obtained.

[0086] In particular, in S201, when extracting statistical seismic wavelets from pre-stack seismic data, seismic traces far away from geological structures should be selected, and the wavelet estimation window should be set near the roof sandstone aquifer.

[0087] Figure 2 shows the 3D seismic data for the designated study area. Figure 2(a) shows the pre-stack 3D seismic data, where the horizontal axis represents the offset distance in meters (m) and the vertical axis represents the two-way travel time in milliseconds (ms). Figure 2(b) shows the post-stack 3D seismic data, where the horizontal axis represents the channel number and the vertical axis represents the two-way travel time in milliseconds (ms).

[0088] Figure 3 The P-wave velocity ratio distribution of the roof sandstone aquifer in the designated study area, obtained by inverting the above steps, is shown. The horizontal axis represents the X-axis, in meters (m), and the vertical axis represents the Y-axis, in meters (m). The color indicates the P-wave velocity ratio. The P-wave velocity ratio of the roof sandstone ranges from 1.64 to 1.78. The P-wave velocity ratio distribution shows a distinct transitional distribution with little regional variability. High values ​​are primarily located in the north-central and southern parts of the country, potentially indicating water-rich areas.

[0089] S30: Perform spectral decomposition on the post-stack seismic data using time-frequency analysis technology to obtain the post-stack seismic time-frequency spectrum, and perform spectral equalization on the time-frequency spectrum;

[0090] Specifically, in step S30, the time-frequency analysis of post-stack seismic data includes the following steps:

[0091] S301: Perform smoothed pseudo-Wigner transform (SPWVD) on each seismic record of the post-stack seismic data to obtain amplitude spectrum data. The formula is as follows:

[0092]

[0093] Among them, g(u) is the smoothing window function in the time direction, h(τ) is the smoothing window function in the delay direction, x(t) is the earthquake record to be analyzed, and x * (t) is the conjugate complex number of the earthquake record to be analyzed.

[0094] S302: Perform spectrum equalization processing on the amplitude spectrum data obtained in S301. The formula is as follows:

[0095]

[0096] in, f after spectral equalization m The time spectrum, S(t,f m ) is f m Time spectrum, W(f m ) is the spectrum equalization coefficient, which is calculated by the following formula:

[0097]

[0098] Among them, Max[S(f ref )] is the maximum amplitude value at the reference frequency, Max[S(f m )] is the frequency f m The maximum amplitude value under .

[0099] In this embodiment, according to the spectrum analysis results, the earthquake main frequency 55HZ is selected as the reference frequency f ref .

[0100] S303: Repeat steps S301 and S302 until the time-frequency spectrum of all seismic gathers after spectrum equalization is obtained. Store the time-frequency spectrum data into a matrix.

[0101] S304: extracting m spectrum decomposition data after spectrum equalization in the target frequency band from the time-frequency spectrum result of S303.

[0102] In this embodiment, according to the magnitude of the earthquake main frequency, 9 spectrum decomposition data after spectrum equalization in the range of 30 to 70 Hz (with a step size of 5 Hz) are extracted from the time-frequency spectrum result of S303.

[0103] Figure 4The figure shows the spectrum decomposition section of the 3D seismic data after spectral equalization in the designated study area. The horizontal axis is the channel number, the vertical axis is the two-way travel time, and the unit is milliseconds (ms). The color indicates the amplitude, and the arrow indicates the location of the roof sandstone aquifer. Figure 4 (a), Figure 4 (b) Figure 4 (c) Figure 4 (d) Spectral decomposition results corresponding to 30Hz, 40Hz, 50Hz and 60Hz respectively. After spectral equalization processing, the amplitudes of each frequency band are unified to the same order of magnitude, and the amplitude difference at the target layer is not much.

[0104] S40: Based on the reflection coefficient approximation formula, a P-wave dispersion attribute inversion method is established to invert the frequency domain post-stack seismic data to obtain the P-wave dispersion attributes of the target layer.

[0105] In particular, in step S40, the approximate formula for the reflection coefficient at near-vertical incidence is:

[0106]

[0107] Where R(θ) is the reflection coefficient at the incident angle θ, v p is the average longitudinal wave velocity in the vertical direction of the media on both sides of the interface, v s is the mean shear wave velocity in the vertical direction of the media on both sides of the interface, ρ is the mean density of the media on both sides of the interface, and Δ represents the parameter difference of the media on both sides of the interface.

[0108] Since the post-stack data mainly represents the reflection response of approximately vertical incidence and the incident angle is small, the above formula can be further simplified as:

[0109]

[0110] Since the reflection coefficient is frequency dependent, assuming V p The density is not affected by the frequency change, so:

[0111]

[0112] Select the earthquake main frequency as the reference frequency, and the above formula is ref The first-order Taylor expansion is:

[0113]

[0114] The frequency-dependent reflection coefficient equation at frequency f is compared with the reference frequency f ref Subtracting the frequency-dependent reflection coefficient equation under , we get:

[0115]

[0116] Longitudinal wave dispersion property Dp is defined as The above formula can be written in matrix form:

[0117]

[0118] Since seismic records are obtained in actual seismic exploration, the reflection coefficient is converted into a spectrum after spectral equalization, and the above formula can be rewritten as:

[0119]

[0120] The ridge regression algorithm is used to solve the above equation and obtain the dispersion attribute D p :

[0121] D p =2(G T G+λI) -1 G T d (11)

[0122] Among them, G is G T is the transposed matrix of G, λ is the damping parameter, and d is f m represents the mth target frequency.

[0123] Figure 5 The figure shows the planar distribution of the P-wave dispersion attribute for the water-bearing sandstone layer in the roof of the designated study area. The horizontal axis is the X-axis, in meters (m), and the vertical axis is the Y-axis, in meters (m). The color indicates the dispersion attribute value. P-wave dispersion attribute values ​​range from 0 to 1 and are unevenly distributed across the study area, with strong spatial variability. Most areas have low attribute values ​​(<0.3), while high-value areas are scattered, primarily concentrated in the central and eastern regions.

[0124] S50: The entropy weights of the P-wave and S-wave velocity ratios and P-wave dispersion properties of the water-bearing sandstone layers in the study area are calculated according to the entropy weight method, and a calculation formula for the water-richness indicator factor is constructed.

[0125] Specifically, in step S50, the entropy weight method for calculating the weight values ​​of the P-wave velocity ratio and the P-wave dispersion attribute of sandstone includes the following steps:

[0126] S501: Standardize each indicator:

[0127]

[0128] Among them, x ij Represents the observation value of the i-th object under the j-th indicator.

[0129] S502: Calculate the proportion of each sample under each indicator:

[0130]

[0131] Where m represents the number of evaluation objects.

[0132] S503: Calculate the information entropy E of each indicator j :

[0133]

[0134] S504: Calculate the weight D of each indicator using the entropy weight method j :

[0135]

[0136] Where n represents the number of indicators.

[0137] In this embodiment, the calculated entropy weight of the ratio of longitudinal and transverse wave velocities is 0.21, and the entropy weight of the longitudinal wave dispersion attribute is 0.79.

[0138] S60: Construct a calculation formula for the water-richness indicator factor and evaluate the water-richness of the sandstone aquifer based on the size of the water-richness indicator factor.

[0139] Specifically, in step S60, the water-richness indicator factor WI of the i-th sample point i The calculation formula is:

[0140]

[0141] Among them, D j is the entropy weight value of the jth indicator, I ij is the jth index value of the i-th sample point.

[0142] Figure 6 The figure shows the planar distribution of the water-rich indicator factor of the target water-bearing sandstone layer, where the horizontal axis is the X-coordinate, in meters (m), the vertical axis is the Y-coordinate, in meters (m), and the color indicates the size of the water-rich indicator factor. The water-rich indicator factor value ranges from 0 to 1. The planar distribution of the water-rich indicator factor is similar to the distribution pattern of the P-wave dispersion attribute value, with strong spatial differences. The attribute values ​​in most areas are low (<0.3), and the high-value areas are scattered. The water-rich indicator factor uses the entropy weight method to simultaneously consider the P-wave velocity ratio and the dispersion attribute, which has a stronger water-rich indicator significance.

[0143] Throughout this specification, references to terms such as "one embodiment," "example," or "specific example" indicate that the specific features, structures, materials, or characteristics described in conjunction with that embodiment or example are included in at least one embodiment or example of the present invention. In this specification, schematic representations of these terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in any one or more embodiments or examples.

[0144] The above contents are merely examples and explanations of the present invention. Those skilled in the art may make various modifications or additions to the described specific embodiments or replace them in similar ways. As long as they do not deviate from the invention or exceed the scope defined by the claims, they should all fall within the scope of protection of the present invention.

Claims

1. A method for evaluating the water richness of a coal seam roof sandstone aquifer based on 3D seismic data, characterized in that: Including steps: S10: Collect borehole data, well logging data and 3D seismic data in the designated study area, including but not limited to borehole lithologic layering data, density logging data, pre-stack 3D seismic data, and post-stack 3D seismic data; S20: Invert pre-stack seismic data to obtain the P-wave and S-wave velocity ratios of the roof sandstone aquifer; S30: Perform spectral decomposition on the post-stack seismic data using time-frequency analysis technology to obtain the post-stack seismic time-frequency spectrum, and perform spectral equalization on the time-frequency spectrum; S40: Based on the reflection coefficient approximation formula, a P-wave dispersion attribute inversion method is established to invert the frequency domain post-stack seismic data to obtain the P-wave dispersion attributes of the target layer; S50: Calculate the entropy weights of the P-wave velocity ratio and P-wave dispersion properties of the water-bearing sandstone layer in the study area according to the entropy weight method; S60: Construct a calculation formula for the water-richness indicator factor and evaluate the water-richness of the sandstone aquifer based on the size of the water-richness indicator factor.

2. The method for evaluating the water richness of a coal seam roof sandstone aquifer based on 3D seismic data according to claim 1, characterized in that: In step S10, the collected pre-stack 3D seismic data undergoes static correction, dynamic correction, combined deconvolution, and pre-stack denoising; the post-stack seismic data is obtained by stacking and interpolating the pre-stack seismic data, and undergoes random noise attenuation, time-varying filtering, amplitude equalization, and other processing.

3. The method for evaluating the water richness of a coal seam roof sandstone aquifer based on 3D seismic data according to claim 1, characterized in that: In step S20, inverting pre-stack seismic data to obtain the P-wave and S-wave velocity ratios of the roof sandstone aquifer includes the following steps: S201: extracting statistical seismic wavelets from pre-stack seismic data, performing well-seismic calibration, and establishing a corresponding relationship between the well logging depth domain and the seismic time domain; S202: establishing an initial inversion model of P-wave velocity, S-wave velocity and density based on the calibrated density and velocity logging curves; S203: Based on the initial inversion model and pre-stack seismic data, the P-wave and S-wave velocity ratio distribution of the roof sandstone aquifer is obtained.

4. The method for evaluating the water-richness of a coal seam roof sandstone aquifer based on 3D seismic data according to claim 3, characterized in that: In step S201, when extracting statistical seismic wavelets from pre-stack seismic data, seismic traces far from geological structures should be selected, and the wavelet estimation window should be set near the roof sandstone aquifer.

5. The method for evaluating the water-richness of a coal seam roof sandstone aquifer based on 3D seismic data according to claim 1, characterized in that: In step S30, the time-frequency analysis of post-stack seismic data includes the following steps: S301: performing a smooth pseudo-Wigner transform on each seismic record of the post-stack seismic data to obtain amplitude spectrum data; S302: Select reference frequency f ref , performing spectrum equalization processing on the amplitude spectrum data obtained in S301; S303: Repeat steps S301 and S302 until the time-frequency spectrum of all seismic gathers after spectrum equalization is obtained. Store the time-frequency spectrum data into a matrix; S304: extracting m spectrum decomposition data after spectrum equalization in the target frequency band from the time-frequency spectrum result of S303.

6. The method for evaluating the water-richness of a coal seam roof sandstone aquifer based on 3D seismic data according to claim 5, characterized in that: In step S301, for a given earthquake record x(t), the calculation formula of the smoothed pseudo-Wigner distribution (SPWVD) is: Among them, g(u) is the smoothing window function in the time direction, h(τ) is the smoothing window function in the delay direction, x(t) is the earthquake record to be analyzed, and x * (t) is the conjugate complex number of the earthquake record to be analyzed.

7. The method for evaluating the water-richness of a coal seam roof sandstone aquifer based on 3D seismic data according to claim 5, characterized in that: In step S302, the purpose of spectral equalization is to eliminate the influence of wavelet superposition, and the processing method is as follows: in, f after spectral equalization m The time spectrum, S(t,f m ) is f m The time spectrum of W(f m ) is the spectrum equalization coefficient, which is calculated by the following formula: Among them, Max[S(f ref )] is the maximum amplitude value at the reference frequency, Max[S(f m )] is the frequency f m The maximum amplitude value under .

8. The method for evaluating the water-richness of a coal seam roof sandstone aquifer based on 3D seismic data according to claim 1, characterized in that: In step S40, the approximate formula of the reflection coefficient at near vertical incidence is: Where R(θ) is the reflection coefficient at the incident angle θ, v p is the average longitudinal wave velocity in the vertical direction of the media on both sides of the interface, v s is the mean shear wave velocity in the vertical direction of the media on both sides of the interface, ρ is the mean density of the media on both sides of the interface, and Δ represents the parameter difference of the media on both sides of the interface; Since the post-stack data mainly represents the reflection response of approximately vertical incidence and the incident angle is small, the above formula can be further simplified as: Since the reflection coefficient is frequency dependent, assuming V p The density is not affected by the frequency change, so: Select the earthquake main frequency as the reference frequency, and the above formula is ref The first-order Taylor expansion is: The frequency-dependent reflection coefficient equation at frequency f is compared with the reference frequency f ref Subtracting the frequency-dependent reflection coefficient equation under , we get: Longitudinal wave dispersion property D p is defined as The above formula can be written in matrix form: Since seismic records are obtained in actual seismic exploration, the reflection coefficient is converted into a spectrum after spectral equalization, and the above formula can be rewritten as: The ridge regression algorithm is used to solve the above equation and obtain the dispersion attribute D p : D p =2(G T G+λI) -1 G T d (11) in, G T is the transposed matrix of G, λ is the damping parameter, and d is f m represents the mth target frequency.

9. The method for evaluating the water-richness of a coal seam roof sandstone aquifer based on 3D seismic data according to claim 1, characterized in that: In step S50, calculating the weight values ​​of the P-wave velocity ratio and the P-wave dispersion attribute of the sandstone includes the following steps: S501: Standardize each indicator: Among them, x ij Represents the observation value of the i-th object under the j-th indicator; S502: Calculate the proportion of each sample under each indicator: Where m represents the number of evaluation objects; S503: Calculate the information entropy E of each indicator j : S504: Calculate the weight D of each indicator using the entropy weight method j : Where n represents the number of indicators.

10. The method for evaluating the water-richness of a coal seam roof sandstone aquifer based on 3D seismic data according to claim 1, characterized in that: In step S60, the water-richness indicator factor WI of the i-th sample point i The calculation formula is: Among them, D j is the entropy weight value of the jth indicator, I ij is the jth index value of the i-th sample point.

Images