A method for evaluating water abundance of sandstone aquifer of coal seam roof based on three-dimensional seismic

Abstract

The application discloses a coal seam roof sandstone aquifer water enrichment evaluation method based on three-dimensional seismic, which comprises the following steps: S10, collecting drilling data, well logging data and high-precision three-dimensional seismic data of a research area; S20, inverting pre-stack seismic data to obtain a P-S wave velocity ratio of the roof sandstone aquifer; S30, performing time-frequency analysis on post-stack seismic data and performing spectral equalization on frequency domain seismic data; S40, inverting a longitudinal wave frequency dispersion attribute of a target layer according to a reflection coefficient approximate formula and the time-frequency spectrum of the post-stack seismic data after the spectral equalization; S50, using an entropy weight method to calculate weight values of the P-S wave velocity ratio and the longitudinal wave frequency dispersion attribute; and S60, constructing a water enrichment indication factor calculation formula and evaluating water enrichment of the sandstone aquifer according to the water enrichment indication factor. The method is based on the theoretical basis of seismic rock physics, has low cost, high operability, more statistical significance, can quickly and accurately evaluate water enrichment of the coal seam roof sandstone aquifer, and provides technical support for development and utilization of coal resources.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of seismic rock physics, and specifically discloses a coal seam roof sandstone aquifer water abundance evaluation method based on three-dimensional seismic. BACKGROUND

[0002] Mine roof water disaster seriously threatens the safety production of coal mines, and the water abundance evaluation of the coal seam roof sandstone aquifer as an important water source of mine water disaster is related to the prevention and control of mine water disaster. At present, the commonly used water abundance evaluation methods in the academic and industrial circles are mainly electromagnetic exploration, which has problems such as shallow detection depth, limited coverage range, insufficient precision, etc. With the increase of exploration depth, these methods cannot realize large-scale quantitative or semi-quantitative roof water disaster prediction. Seismic exploration is the most commonly used method to obtain underground medium information, and has the advantages of wide coverage and high precision, providing a new way for water abundance identification.

[0003] The ratio of P-wave and S-wave velocity is an important elastic parameter reflecting the property of stratum fluid. The existence of fluid will increase the P-wave velocity of sandstone, but has little effect on the S-wave velocity, so as to increase the ratio of P-wave and S-wave velocity. In addition, under the influence of fluid, the seismic wave shows significant dispersion and attenuation when passing through the sandstone aquifer, which lays a theoretical foundation for the fluid identification and quantitative evaluation of water abundance of the aquifer.

[0004] The entropy weight method is a typical objective weighting method. This method measures the information amount of each index by calculating the data dispersion degree of each index, and then determines the weight, so as to construct a more stable and reasonable weight system.

[0005] The water abundance prediction technology based on three-dimensional seismic is applied to the prediction of roof sandstone water abundance, which can consider the P-wave and S-wave velocity ratio and dispersion attribute information of the sandstone aquifer at the same time, avoid the limitation of single attribute interference by noise or other factors, and realize more reliable water abundance evaluation. SUMMARY

[0006] In view of the above problems existing in the prior art, the present application provides a coal seam roof sandstone aquifer water abundance evaluation method based on three-dimensional seismic. The above technical purposes can be achieved and other technical effects can be brought by adopting the following technical features.

[0007] To solve the above problems, the technical scheme adopted by the present application is as follows:

[0008] A coal seam roof sandstone aquifer water abundance evaluation method based on multi-source data, comprising the following steps:

[0009] S10: Collecting drilling data, logging data and three-dimensional seismic data of a specified study area, including but not limited to drilling lithology layering data, density logging data, pre-stack three-dimensional seismic data and post-stack three-dimensional seismic data.

[0010] S20: Inversion of pre-stack seismic data to obtain the P-S wave velocity ratio of the top sandstone aquifer;

[0011] S30: Spectral decomposition of post-stack seismic data by time-frequency analysis technique to obtain the time-frequency spectrum of post-stack seismic data, and spectral equalization of the time-frequency spectrum;

[0012] S40: Establishment of P-wave dispersion attribute inversion method according to the reflection coefficient approximation formula, inversion of frequency domain post-stack seismic data to obtain the P-wave dispersion attribute of the target layer;

[0013] S50: Calculation of entropy weight values of the P-S wave velocity ratio and P-wave dispersion attribute of the water-bearing sandstone layer in the study area according to the entropy weight method;

[0014] S60: Construction of water enrichment indicator calculation formula, and evaluation of the water enrichment of the sandstone aquifer according to the size of the water enrichment indicator.

[0015] Further, the pre-stack three-dimensional seismic data in the S10 step is subjected to static correction, dynamic correction, combined deconvolution, and pre-stack denoising processing; the post-stack seismic data is obtained by stacking and interpolation of the pre-stack seismic data, and is subjected to random noise attenuation, time-varying filtering, amplitude equalization, etc.

[0016] Further, the S20 step of inversion of pre-stack seismic data to obtain the P-S wave velocity ratio of the top sandstone aquifer comprises the following steps:

[0017] S201: Extraction of statistical seismic wavelet of pre-stack seismic data, well-to-seismic calibration, and establishment of corresponding relationship between logging depth domain and seismic time domain;

[0018] S202: Establishment of initial inversion model of P-wave velocity, S-wave velocity and density according to the calibrated density and velocity logging curves;

[0019] S203: Obtaining the P-S wave velocity ratio distribution of the top sandstone aquifer according to the initial inversion model and the pre-stack seismic data.

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

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

[0022] S301: Smooth pseudo Wigner-Ville distribution (SPWVD) of each seismic record of post-stack seismic data to obtain amplitude spectrum data;

[0023] S302: Selection of reference frequency f refThe amplitude spectrum data obtained in S301 is subjected to spectral equalization processing;

[0024] S303: repeating the steps of S301 and S302 until the time-frequency spectrum S(t,f) of all seismic trace sets after spectral equalization is obtained, and storing the time-frequency spectrum data into a matrix;

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

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

[0027]

[0028] wherein g(u) is a smoothing window function in the time direction, h(τ) is a smoothing window function in the delay direction, x(t) is the seismic record to be analyzed, x * (t) is the conjugate complex of the seismic record to be analyzed.

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

[0030]

[0031] wherein S(t,f m ) is the time-frequency spectrum of f m , S(t,f m ) is the time-frequency spectrum of f m , and W(f ) is the spectral equalization coefficient, which is calculated by the following formula:

[0032]

[0033] wherein Max[S(f ref )] is the maximum amplitude value at the reference frequency, and Max[S(f m )] is the maximum amplitude value at the frequency f m .

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

[0035]

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

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

[0038]

[0039] Since the reflection coefficient has frequency dependence, assuming V p With the change of frequency, the density is not affected by the change of frequency, then:

[0040]

[0041] Select the main frequency of the earthquake as the reference frequency, and perform first-order Taylor expansion on the above formula at the reference frequency f ref :

[0042]

[0043] Subtract the frequency-dependent reflection coefficient equation at frequency f from the frequency-dependent reflection coefficient equation at the reference frequency f ref :

[0044]

[0045] The P-wave dispersion attribute D p is defined as The above formula can be written in the form of a matrix:

[0046]

[0047] Since the seismic record obtained in actual seismic exploration is the frequency spectrum after the reflection coefficient is converted to the spectrum, the above formula can be rewritten as:

[0048]

[0049] The ridge regression algorithm is used to solve the above formula, and the dispersion attribute D p is obtained:

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

[0051] Where G is G T , λ is the damping parameter, and d is f m mth target frequency.

[0052] Further in step S50, the entropy weight method calculates the weight value of the sandstone P-wave and S-wave velocity ratio and P-wave dispersion attribute, including the following steps:

[0053] S501: standardize each index:

[0054]

[0055] Wherein, x ij represents the observation value of the i-th object under the j-th index;

[0056] S502: calculate the proportion value of each sample under each index:

[0057]

[0058] Wherein, m represents the number of evaluation objects;

[0059] S503: calculate the information entropy E j of each index:

[0060]

[0061] S504: calculate the weight D j of each index by entropy weight method:

[0062]

[0063] Wherein, n represents the number of indexes.

[0064] Further in step S60, the calculation formula of the water enrichment indication factor WI i of the i-th sample point is:

[0065]

[0066] Wherein, D j is the entropy weight value of the j-th index, and I ij is the j-th index value of the i-th sample point.

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

[0068] 1. The water enrichment prediction technology based on three-dimensional seismic is applied to the roof sandstone water enrichment prediction, which can consider the P-wave and S-wave velocity ratio and dispersion attribute information of the water-bearing sandstone layer at the same time, avoid the limitation of single attribute interference by noise or other factors, and realize more reliable water enrichment evaluation.

[0069] 2. The entropy weight method is adopted to calculate the data dispersion degree of each index to measure the information amount, and then the weight is determined, so that a more stable and reasonable weight system can be constructed.

[0070] 3、The method has low cost, strong operability, more statistical significance, and can quickly and accurately evaluate the water enrichment condition of the sandstone aquifer of the coal seam roof, thereby providing technical support for the development and utilization of coal resources. BRIEF DESCRIPTION OF DRAWINGS

[0071] Figure 1 It is a flowchart of the coal seam roof sandstone aquifer water enrichment evaluation method based on three-dimensional seismic of the application.

[0072] Figure 2 is three-dimensional seismic data, Figure 2(a) is pre-stack three-dimensional seismic data, and Figure 2(b) is post-stack three-dimensional seismic data.

[0073] Figure 3 It is a P-wave and S-wave velocity ratio distribution map obtained by pre-stack seismic inversion.

[0074] Figure 4 It is a spectrum decomposition profile of three-dimensional seismic data after spectrum equalization, wherein, 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 It is a frequency dispersion attribute plane distribution map obtained by inversion.

[0076] Figure 6 It is a water enrichment indication factor plane distribution map. DETAILED DESCRIPTION

[0077] In order to make the purpose, technical scheme and advantages of the technical scheme of the application more clear, the technical scheme of the embodiment of the application will be described clearly and completely in the following combined with the drawings of the specific embodiment of the application. It should be noted that the described embodiment is a part of the embodiment of the application, not all the embodiments. Based on the described embodiment of the application, all other embodiments obtained by those skilled in the art without creative labor belong to the protection scope of the application.

[0078] The specific embodiments of the application will be further described below combined with the drawings, as shown in Figure 1 The steps include:

[0079] S10: Collecting drilling data, logging data and three-dimensional seismic data of the specified study area, including but not limited to drilling lithology layering data, density logging data, pre-stack three-dimensional seismic data, post-stack three-dimensional seismic data.

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

[0081] S20: Inversion of the pre-stack seismic data to obtain the P-S wave velocity ratio of the top sandstone aquifer.

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

[0083] S201: Extraction of the statistical seismic wavelet of the pre-stack seismic data, well-seismic calibration, and establishment of the corresponding relationship between the logging depth domain and the seismic time domain;

[0084] S202: Establishment of the initial inversion model of P wave velocity, S wave velocity and density according to the calibrated density and velocity logging curves;

[0085] S203: Obtaining the P-S wave velocity ratio distribution of the top sandstone aquifer according to the initial inversion model and the pre-stack seismic data.

[0086] In particular, in S201, the statistical seismic wavelet of the pre-stack seismic data should be extracted by selecting seismic traces far from geological structures, and the wavelet estimation window should be set near the top sandstone aquifer.

[0087] Figure 2 is the 3D seismic data of the specified research area, Figure 2(a) is the pre-stack 3D seismic data, the horizontal coordinate is offset, the unit is meter (m), and the vertical coordinate is two-way travel time, the unit is millisecond (ms); Figure 2(b) is the post-stack 3D seismic data, the horizontal coordinate is trace number, and the vertical coordinate is two-way travel time, the unit is millisecond (ms).

[0088] Figure 3 The P-S wave velocity ratio distribution of the top sandstone aquifer of the specified research area obtained by inversion according to the above steps, wherein the horizontal coordinate is X coordinate, the unit is meter (m), the vertical coordinate is Y coordinate, the unit is meter (m), the color indicates the size of the P-S wave velocity ratio, and the P-S wave velocity ratio of the top sandstone is between 1.64 and 1.78. The P-S wave velocity ratio distribution presents obvious transition distribution characteristics, the regional numerical variability is not large, and the high value is mainly distributed in the north-central and southern parts, which may indicate the water-rich area.

[0089] S30: Spectral decomposition of the post-stack seismic data by time-frequency analysis technology to obtain the post-stack seismic time-frequency spectrum, and spectral equalization of the time-frequency spectrum;

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

[0091] S301: performing a smoothing pseudo Wigner-Ville distribution (SPWVD) on each seismic record of the stacked seismic data to obtain amplitude spectrum data, according to the following formula:

[0092]

[0093] wherein g(u) is a smoothing window function in the time direction, h(τ) is a smoothing window function in the delay direction, x(t) is the seismic record to be analyzed, x * (t) is the conjugate complex of the seismic record to be analyzed.

[0094] S302: performing spectral equalization on the amplitude spectrum data obtained in S301, according to the following formula:

[0095]

[0096] wherein S(t,f ) is the time-frequency spectrum of f m after spectral equalization, S(t,f m ) is the time-frequency spectrum of f m , and W(f m ) is a spectral equalization coefficient, calculated according to the following formula:

[0097]

[0098] wherein Max[S(f ref )] is the maximum amplitude value at the reference frequency, and Max[S(f m )] is the maximum amplitude value at the frequency f m .

[0099] In this embodiment, the main seismic frequency 55 Hz is selected as the reference frequency f ref according to the spectral analysis result.

[0100] S303: repeating the steps of S301 and S302 until the time-frequency spectrum S(t,f of all seismic gathers after spectral equalization is obtained.

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

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

[0103] Figure 4The spectral decomposition profile of the specified study area after spectral equalization is shown in the figure, where the horizontal coordinate is the trace number, the vertical coordinate is the two-way travel time, the unit is millisecond (ms), the color indicates the amplitude, and the arrow indicates the position of the top sandstone aquifer. Figure 4 (a), Figure 4 (b), Figure 4 (c), Figure 4 (d) respectively correspond to the spectral decomposition results of 30Hz, 40Hz, 50Hz and 60Hz. After spectral equalization, the amplitudes of each frequency band are unified to the same order of magnitude, and the amplitude difference at the target layer is not large.

[0104] S40: According to the approximate formula of reflection coefficient, the longitudinal wave dispersion attribute inversion method is established, the post-stack seismic data in frequency domain is inverted, and the longitudinal wave dispersion attribute of the target layer is obtained.

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

[0106]

[0107] Where R(θ) is the reflection coefficient at incidence angle θ, v p is the average longitudinal wave velocity of the media on both sides of the interface in the vertical direction, v s is the average shear wave velocity of the media on both sides of the interface in the vertical direction, ρ is the average 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 represent the reflection response of near-normal incidence, the incidence angle is small, and the above formula can be further simplified as:

[0109]

[0110] Since the reflection coefficient has frequency dependence, assuming V p changes with frequency, and the density is not affected by frequency, then:

[0111]

[0112] Selecting the seismic main frequency as the reference frequency, the first-order Taylor expansion of the above formula at the reference frequency f ref is:

[0113]

[0114] Subtracting the frequency-dependent reflection coefficient equation at frequency f from the frequency-dependent reflection coefficient equation at reference frequency f ref , we get:

[0115]

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

[0117]

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

[0119]

[0120] The ridge regression algorithm is used to solve the above formula, and the dispersion attribute D is obtained p :

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

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

[0123] Figure 5 The longitudinal wave dispersion attribute plane distribution of the water-bearing sandstone layer in the specified research area is shown, wherein the horizontal coordinate is the X coordinate, the unit is meter (m), the vertical coordinate is the Y coordinate, the unit is meter (m), and the color indicates the size of the dispersion attribute value; The longitudinal wave dispersion attribute value is between 0 and 1, the longitudinal wave dispersion attribute value is unevenly distributed in the whole research area, and the spatial difference is strong. Most of the attribute values are low (<0.3), the high value area is sporadic, and is mainly concentrated in the middle and east.

[0124] S50: According to the entropy weight method, the entropy weight values of the longitudinal and transverse wave velocity ratio and the longitudinal wave dispersion attribute of the water-bearing sandstone layer in the research area are calculated, and a calculation formula of the water enrichment indicating factor is constructed.

[0125] In particular, in step S50, the entropy weight method calculates the weight value of the sandstone longitudinal and transverse wave velocity ratio and the longitudinal wave dispersion attribute, including the following steps:

[0126] S501: Standardize each index:

[0127]

[0128] Where, x ij represents the observation value of the ith object under the jth index.

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

[0130]

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

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

[0133]

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

[0135]

[0136] wherein n represents the number of indexes.

[0137] In this embodiment, the entropy weight value of the P-wave / S-wave velocity ratio is 0.21, and the entropy weight value of the P-wave dispersion attribute is 0.79.

[0138] S60: Construct a water enrichment indicator calculation formula, and evaluate the water enrichment of the sandstone aquifer according to the size of the water enrichment indicator.

[0139] In particular, in step S60, the water enrichment indicator WI i of the i-th sample point is calculated according to the following formula:

[0140]

[0141] wherein D j is the entropy weight value of the j-th index, and I ij is the j-th index value of the i-th sample point.

[0142] Figure 6 The water enrichment indicator plane distribution of the target water-bearing sandstone layer is shown in the figure, wherein the horizontal coordinate is the X coordinate, the unit is meter (m), the vertical coordinate is the Y coordinate, the unit is meter (m), and the color indicates the size of the water enrichment indicator; the water enrichment indicator value is between 0 and 1, the water enrichment indicator plane distribution is similar to the P-wave dispersion attribute value distribution law, and the spatial difference is strong, the attribute value in most areas is low (<0.3), the high value area is sporadic, and the water enrichment indicator considers the P-wave / S-wave velocity ratio and the dispersion attribute at the same time by using the entropy weight method, and has stronger water enrichment indicating significance.

[0143] In the description of the specification, reference to terms "one embodiment", "an example", "a specific example" and so on is intended to indicate that a particular feature, structure, material or characteristic described in connection with the embodiment or example is included in at least one embodiment or example of the application. Descriptive expressions of the above terms in the specification do not necessarily refer to the same embodiment or example. Moreover, the described specific features, structures, materials or characteristics can be combined in any one or more embodiments or examples in a suitable manner.

[0144] The above is only an example and illustration of the application, and those skilled in the art can make various modifications or supplements to the described specific embodiments or replace them with similar ways, as long as they do not deviate from the application or exceed the scope defined by the claims.

Claims

1. A method for evaluating the water-bearing capacity of sandstone aquifers in coal seam roofs based on three-dimensional seismic data, characterized in that, Including the following steps: S10: Collect borehole data, well logging data and 3D seismic data for the designated study area, including but not limited to borehole lithological stratification 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 ratio of the top sandstone aquifer; S30: Perform spectral decomposition on post-stack seismic data using time-frequency analysis techniques to obtain the time spectrum of post-stack seismic data, and then perform spectral equalization on the time spectrum; S40: Based on the approximate formula for the reflection coefficient, a method for inverting the P-wave dispersion attribute is established to invert post-stack seismic data in the frequency domain and obtain the P-wave dispersion attribute 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 layers in the study area using the entropy weight method; S60: Construct a formula for calculating the water-bearing indicator factor, and evaluate the water-bearing capacity of sandstone aquifers based on the magnitude of the water-bearing indicator factor; In step S30, the time-frequency analysis of the post-stack seismic data includes the following steps: S301: Perform 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 The amplitude spectrum data obtained in S301 is subjected to spectral equalization processing. S303: Repeat steps S301 and S302 until the time spectrum after equalization of all seismic gather spectra is obtained. The time-spectrum data is stored in the matrix; S304: Extract the m spectral decomposition data of the target frequency band after spectral equalization from the time-spectrum results of S303; In step S60, the first i Water-rich indicator factors for each sample point The calculation formula is: (1) in, D j For the first j The entropy weight value of each indicator I ij For the first i The first sample point j Individual indicator values.

2. The method for evaluating the water-bearing capacity of sandstone aquifers in coal seam roofs based on three-dimensional seismic analysis 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 superimposing 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-bearing capacity of sandstone aquifers in coal seam roofs based on three-dimensional seismic analysis according to claim 1, characterized in that: In step S20, the inversion of pre-stack seismic data to obtain the P-wave and S-wave velocity ratio of the top sandstone aquifer includes the following steps: S201: Extract statistical seismic wavelets from pre-stack seismic data, perform well-seismic calibration, and establish the correspondence between the well logging depth domain and the seismic time domain; S202: Establish an initial inversion model for 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 top sandstone aquifer was obtained.

4. The method for evaluating the water-bearing capacity of sandstone aquifers in coal seam roofs based on three-dimensional seismic analysis 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 sandstone aquifer in the top plate.

5. The method for evaluating the water-bearing capacity of sandstone aquifers in coal seam roofs based on three-dimensional seismic analysis according to claim 1, characterized in that: In step S301, for a given seismic record x(t) The formula for calculating the smoothed pseudo-Wigner distribution (SPWVD) is as follows: (2) in, For the time-direction smoothing window function, For the smoothing window function in the delay direction, x(t) For the earthquake record to be analyzed ，x * (t) Let be the conjugate complex number of the earthquake record to be analyzed.

6. The method for evaluating the water-bearing capacity of sandstone aquifers in coal seam roofs based on three-dimensional seismic data according to claim 1, characterized in that: In step S302, the purpose of spectral equalization is to eliminate the influence of wavelet overprinting. The processing method is as follows: (3) in, After spectral equalization f m The time spectrum, for f m The time spectrum, The spectral equalization coefficient is calculated using the following formula: (4) in, The maximum amplitude value at the reference frequency. For frequency f m The maximum amplitude value.

7. The method for evaluating the water-bearing capacity of sandstone aquifers in coal seam roofs based on three-dimensional seismic analysis according to claim 1, characterized in that: In step S40, the approximate formula for the reflection coefficient at near-perpendicular incidence is: (5) in, R ( θ ( ) is the angle of incidence θ The reflectance coefficient below, v p The mean longitudinal wave velocity in the direction perpendicular to the medium on both sides of the interface is denoted as . v s The average transverse wave velocity in the direction perpendicular to the medium on both sides of the interface. ρ The average density of the media on both sides of the interface. This represents the parameter difference between the media on both sides of the interface; Since the post-stack data mainly represents the reflection response of approximately perpendicular incidence, with a relatively small incident angle, the above equation can be further simplified to: (6) Since the reflection coefficient has frequency dependence, assuming V p If the density is unaffected by frequency changes, then: (7) Choosing the dominant earthquake frequency as the reference frequency, the above equation is applied at the reference frequency. f ref Performing a first-order Taylor expansion at the point, we have: (8) frequency f The frequency-varying reflection coefficient equation and reference frequency f ref Subtracting the equations for the frequency-varying reflection coefficients, we get: (9) Longitudinal wave dispersion properties D p Defined as The above formula can be written in matrix form: (10) Since actual seismic exploration yields seismic records, converting the reflection coefficient into a spectrum after spectral equalization allows the above equation to be rewritten as: (11) The ridge regression algorithm is used to solve the above equation to obtain the dispersion property. D p : (12) in, G for , G T for G The transpose of the matrix, λ Let d be the damping parameter. , f m Indicates the first m Target frequency.

8. The method for evaluating the water-bearing capacity of sandstone aquifers in coal seam roofs based on three-dimensional seismic analysis according to claim 1, characterized in that: In step S50, calculating the weight values ​​of the P-wave / S-wave velocity ratio and the P-wave dispersion property of sandstone includes the following steps: S501: Standardize each indicator: (13) in, x ij Indicates the first i The object in the first j Observations under each indicator; S502: Calculate the proportion of each sample under each indicator: (14) in, m Indicates the number of evaluation objects; S503: Calculate the information entropy of each indicator E j : (15) S504: Calculate the weights of each indicator using the entropy weight method. D j : (16) in, n Indicates the number of indicators.

Images