A three-dimensional seismic exploration method for coalfield in thick loess plateau region

Abstract

The present application relates to geological exploration technical field, disclose a kind of thick loess plateau area coalfield 3D seismic exploration method, comprising the following steps: S101, artificial seismic source shooting record and superficial Rayleigh wave detector data are collected;S102, grid division is carried out to target area, and set to be inverted parameter;S103, vertical ray path time estimation, multilayer medium dispersion forward and poststack energy estimation are carried out;S104, according to first wave arrival time observation value, Rayleigh wave observation dispersion curve, theoretical ray path time, phase velocity and poststack energy generates multi-target error function;S105, target area is divided into three subdomains, differential evolution iteration is carried out, until local convergence;S106, establish a plurality of parallel annealing chain, set chain temperature and gradually decrease by chain, optimize to be inverted parameter, output optimal to be inverted parameter;S107, the static correction time difference and absorption compensation coefficient of each acquisition point are calculated.The present application improves coal seam structure imaging precision and exploration data reliability under complex surface conditions.

Description

Technical Field

[0001] This invention belongs to the field of geological exploration technology, specifically relating to a three-dimensional seismic exploration method for coalfields in thick loess plateau areas. Background Technology

[0002] In the thick loess plateau region, coalfield seismic exploration faces significant challenges due to its unique geological structure and complex sedimentary characteristics. The thick, loose, and unevenly moist loess overburden in this area, coupled with low wave velocities, easily leads to severe distortion in seismic wave propagation, manifesting as rapid attenuation of high-frequency energy, waveform broadening, difficulty in identifying the first wave, and low signal-to-noise ratio. Traditional two-dimensional seismic exploration methods have limited effectiveness in such areas, struggling to accurately identify coal seam occurrence and structural features. Particularly in static correction and high-frequency absorption compensation, the strong lateral heterogeneity of the strata results in significant errors in conventional static correction methods based on surface velocity models, affecting subsequent inversion accuracy. Furthermore, the thick loess layer exhibits a strong absorption effect on high-frequency seismic waves, resulting in weak underground reflection energy and complex dispersion characteristics, further increasing the difficulty of parameter inversion and coal seam interpretation. Therefore, a three-dimensional seismic exploration method adapted to the unique geological background of the thick loess plateau region is urgently needed to improve the ability to accurately identify coalfield geological structures and coal seams, providing more reliable data support for resource development and geological disaster early warning. Summary of the Invention

[0003] This invention provides a three-dimensional seismic exploration method for coalfields in thick loess plateau areas, solving the technical problems in related technologies such as severe distortion of seismic wave propagation signals, large static correction errors, inaccurate modeling of dispersion characteristics, and unstable quality factor inversion in thick loess-covered areas.

[0004] This invention provides a three-dimensional seismic exploration method for coalfields in thick loess plateau areas, comprising the following steps:

[0005] S101, collects artificial source shot-shoal records and shallow Rayleigh wave detector data, and extracts seismic gathers, first wave arrival time observations and Rayleigh wave dispersion curves at different frequencies corresponding to each source location.

[0006] S102, the target area is divided into regular grids according to the preset depth range, and the loess layer thickness, P-wave propagation velocity and quality factor are set as parameters to be inverted for each grid cell;

[0007] S103, based on the parameters to be inverted for each grid cell, vertical ray path time estimation, multilayer medium dispersion forward modeling, and post-stack energy estimation are performed respectively to obtain vertical ray path time, phase velocity, and post-stack energy;

[0008] S104, based on the first wave arrival time observation, Rayleigh wave observation dispersion curve, theoretical ray path time, phase velocity and post-stack energy, a multi-objective error function is generated by weighting according to preset weights;

[0009] S105, based on the thickness of the loess layer, divides the target area into three sub-domains: low, medium, and high. Differential evolution iteration is performed on the parameters to be inverted for the grid cells in each sub-domain until local convergence is achieved.

[0010] S106, establish multiple parallel annealing chains, set the chain temperature based on the loess layer moisture content collected in the field and decrease it chain by chain, optimize the parameters to be inverted through jump variation, energy dissipation constraint correction of quality factor, and inter-chain state exchange mechanism until convergence and output the optimal parameters to be inverted.

[0011] S107: Calculate the static correction time difference of each acquisition point based on the optimal parameters to be inverted, calculate the absorption compensation coefficient based on the parameters to be inverted, and perform static correction and absorption compensation on the seismic gathers.

[0012] Furthermore, the artificial seismic source shot record is the raw amplitude time series data recorded on the ground geophone;

[0013] The seismic gather is an amplitude sequence continuously recorded by the detector at a preset sampling frequency during a single source excitation process;

[0014] The first wave arrival time observation value was extracted by using a threshold discrimination method on the amplitude signal of the seismic trace collection;

[0015] The Rayleigh wave observation dispersion curve is obtained by performing time-spectrum analysis on multiple channels of the same excitation to identify the phase velocity corresponding to the peak value of the Rayleigh wave energy concentration. The window function used in the time-spectrum analysis is the Hanning window.

[0016] Furthermore, the step of estimating the vertical ray path time includes:

[0017] S201: Extract the grid cells through which the seismic waves pass along the vertical propagation path of the seismic waves and form a set of grid cells.

[0018] S202, calculate the ratio of loess layer thickness to P-wave propagation velocity for each grid cell to obtain the one-way travel time;

[0019] S203, sum the one-way travel time of all grid cells in the grid cell set to obtain the vertical ray path time.

[0020] Furthermore, the step of the multilayer dielectric dispersion forward modeling includes:

[0021] S301, for each frequency in the Rayleigh wave observation dispersion curve, construct the transfer matrix of each medium layer based on the Haskell-Thomson method;

[0022] S302, combine the transfer matrices of all layers into a total transfer matrix;

[0023] S303, construct the determinant equation with the determinant of the total transfer matrix equal to zero, and solve it to obtain the phase velocity corresponding to the frequency.

[0024] Furthermore, the seismic traces recorded by each detector under the same excitation are collected, and the amplitude values ​​corresponding to the zero offset time are summed after squaring to obtain the post-stack energy.

[0025] Furthermore, the steps for constructing the multi-objective error function include:

[0026] S401, for each observation pair consisting of the source location and the acquisition point, calculate the squared difference between the first arrival time of the acquired seismic wave and the path time of the vertical ray, and sum the squared differences of all observation pairs to obtain the sum of squared errors of the first wave arrival time;

[0027] S402, calculate the squared difference between the P-wave propagation velocity and the phase velocity at each frequency, and sum the squared differences for all frequencies to obtain the sum of squared Rayleigh wave dispersion errors;

[0028] S403 linearly combines the sum of squared errors of the first wave arrival time, the sum of squared errors of the Rayleigh wave dispersion, and the post-stack energy to obtain a multi-objective error function.

[0029] Furthermore, the parameters to be inverted for the grid cells within each subdomain are subjected to differential evolution iterations until local convergence. Specific steps include:

[0030] S501, in each subdomain, multiple candidate individuals are randomly generated; each candidate individual is composed of the loess layer thickness, P-wave propagation velocity and quality factor of all grid cells in that subdomain;

[0031] S502, For any candidate individual, generate a mutation vector based on the individual in the subdomain that is closest to the actual measured thickness and three random candidate individuals;

[0032] S503: Dynamically adjust the crossover probability based on the average quality factor of all candidate individuals in the subdomain, and cross the mutation vector with the original candidate individuals according to the probability in the corresponding grid cell parameter dimension to generate experimental candidate individuals;

[0033] S504. For each experimental candidate individual, calculate the sum of squares of the first wave arrival error and the sum of squares of the Rayleigh wave dispersion error, and perform a weighted summation to obtain the local fitness. If the local fitness of the experimental candidate individual is better than that of the original candidate individual, then replace it; otherwise, retain the original candidate.

[0034] S505, after the first preset number of iterations, select a certain proportion of the best candidate individuals from each subdomain and put them into the global elite library; for each elite individual in the elite library, if its Euclidean distance with any candidate individual in any subdomain is lower than the first preset threshold, then replace the candidate individual with the largest Euclidean distance in that subdomain with that elite individual.

[0035] S506: When the change in the local fitness of the best candidate individual in a subdomain is lower than the second preset threshold and the difference in loess layer thickness between adjacent grid cells in the subdomain is lower than the third preset threshold, the subdomain is determined to be locally converged, and the differential evolution in the subdomain ends.

[0036] Furthermore, the steps for generating the mutation vector include:

[0037] S601, three different candidate individuals are randomly selected, denoted as random individual one, random individual two and random individual three, and the candidate individual with the smallest deviation from the measured thickness is selected from all candidate individuals in the subdomain, denoted as the individual with the best thickness fit;

[0038] S602, using all the parameters to be inverted of random individual one as the initial reference, calculate the difference between the thickness fit optimal individual and the parameters to be inverted of random individual one in each grid cell, multiply the difference by the preset scaling factor and combine it with the initial reference;

[0039] S603, calculate the difference between the parameters to be inverted for random individual two and random individual three in each grid cell, multiply the difference by the preset scaling factor, and combine it with the result of S602;

[0040] S604 is the result of summing the results of S603 to obtain the mutation vector.

[0041] Furthermore, the specific steps in S106 include:

[0042] S701, establish multiple parallel annealing chains within the target area, with each annealing chain corresponding to a candidate individual;

[0043] S702, based on the moisture content of the loess layer collected on site, set the initial temperature of the first annealing chain, and the temperature of subsequent annealing chains decreases by a fixed ratio;

[0044] S703, for each candidate individual corresponding to an annealing chain, a perturbation vector is generated using Levy flight and superimposed on the candidate individual to obtain a mutated candidate individual;

[0045] S704, for the quality factor of each grid cell in the mutant candidate individuals, calculate the viscoelastic dissipation energy increment based on the loess layer thickness and P wave propagation velocity, and correct the quality factor accordingly.

[0046] S705, compare the error values ​​of the global multi-objective error function of the mutated candidate individual with those of the original candidate individual. If the error value of the mutated candidate individual is not higher than that of the original candidate individual, then the mutated candidate individual is directly accepted.

[0047] S706, every second preset iteration number, attempt to swap candidate individuals between two adjacent annealing chains. If the swap reduces the sum of errors of the two chains, then swap; otherwise, perform a forced swap with the first preset probability. When the error change of all annealing chains is lower than the fourth preset threshold in multiple consecutive rounds, output the optimal parameters to be inverted.

[0048] Furthermore, the static correction time difference is calculated based on the P-wave propagation velocity and loess layer thickness of each grid cell in the optimal parameters to be inverted; the absorption compensation coefficient is determined based on the attenuation relationship between the quality factor and the frequency; wherein, static correction is used to adjust the time axis of the seismic gather by accumulating the time difference on the path between each source and the receiver, and absorption compensation is achieved by applying an inverse attenuation function to the spectrum of the seismic gather.

[0049] The beneficial effects of this invention are as follows: This invention integrates multi-dimensional data such as first wave arrival time, Rayleigh wave dispersion curves, and post-stack energy, and improves the stability and accuracy of underground medium parameter inversion through joint optimization of multi-objective error functions; it divides subdomains based on loess layer thickness, combines differential evolution and parallel annealing algorithms to balance global search and local exploitation, adapts to the strong heterogeneity of thick loess plateau areas, and avoids getting trapped in local optima; it achieves static correction and time difference compensation, absorption compensation based on the parameters to be inverted, eliminates near-surface interference, broadens the frequency band, and improves seismic data resolution and imaging quality, providing a reliable data foundation for coalfield exploration. Attached Figure Description

[0050] Figure 1 This is a flowchart of a three-dimensional seismic exploration method for coalfields in thick loess plateau areas according to the present invention. Detailed Implementation

[0051] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, features described in some examples may be combined in other examples.

[0052] It should be noted that, unless otherwise defined, the technical or scientific terms used in one or more embodiments of the present invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in one or more embodiments of the present invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed after the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0053] like Figure 1 As shown, a three-dimensional seismic exploration method for coalfields in thick loess plateau areas includes the following steps:

[0054] S101, collects artificial source shot-shoal records and shallow Rayleigh wave detector data, and extracts seismic gathers, first wave arrival time observations and Rayleigh wave dispersion curves at different frequencies corresponding to each source location.

[0055] S102, the target area is divided into regular grids according to the preset depth range, and the loess layer thickness, P-wave propagation velocity and quality factor are set as parameters to be inverted for each grid cell;

[0056] S103, based on the parameters to be inverted for each grid cell, vertical ray path time estimation, multilayer medium dispersion forward modeling, and post-stack energy estimation are performed respectively to obtain vertical ray path time, phase velocity, and post-stack energy;

[0057] S104, based on the first wave arrival time observation, Rayleigh wave observation dispersion curve, theoretical ray path time, phase velocity and post-stack energy, a multi-objective error function is generated by weighting according to preset weights;

[0058] S105, based on the thickness of the loess layer, divides the target area into three sub-domains: low, medium, and high. Differential evolution iteration is performed on the parameters to be inverted for the grid cells in each sub-domain until local convergence is achieved.

[0059] S106, establish multiple parallel annealing chains, set the chain temperature based on the loess layer moisture content collected in the field and decrease it chain by chain, optimize the parameters to be inverted through jump variation, energy dissipation constraint correction of quality factor, and inter-chain state exchange mechanism until convergence and output the optimal parameters to be inverted.

[0060] S107: Calculate the static correction time difference of each acquisition point based on the optimal parameters to be inverted, calculate the absorption compensation coefficient based on the parameters to be inverted, and perform static correction and absorption compensation on the seismic gathers.

[0061] In one embodiment of the present invention, a controllable seismic source or an explosive seismic source is used to excite seismic waves with controllable energy. The acquisition points are arranged at equal intervals along the survey line to ensure coverage of the target area. High-sensitivity velocity detectors are densely arranged at equal intervals along the survey line to form a shallow observation array for capturing surface wave signals such as Rayleigh waves.

[0062] After a single source excitation, the geophone synchronously collects surface vibration signals, and the recording time covers 0.5 to 2 seconds after the source excitation, forming the original amplitude time series data, i.e., the artificial source shot-vibration record; the data is stored in real time in binary format, including metadata such as timestamp, geophone coordinates, and amplitude value.

[0063] Seismic gathers are the amplitude sequences continuously recorded by the geophone during a single source excitation process at a preset sampling frequency;

[0064] The first wave arrival time observation is extracted using a threshold discrimination method. First, a 10-80Hz bandpass filter is applied to the single-channel seismic gather to remove low-frequency noise and high-frequency interference. The root mean square amplitude of the background noise in the first 50ms of the preprocessed signal is calculated, and a threshold of 2.5×RMS is set. The time point when the amplitude first exceeds this threshold is defined as the first wave arrival time. False triggers are corrected by manual verification to ensure that the picking error does not exceed 0.005 seconds.

[0065] The Rayleigh wave observation dispersion curve is generated based on time-spectrum analysis of multichannel gathers. A Hanning window function is applied to the multichannel data of the same excitation to suppress spectral leakage. A time-frequency matrix is ​​generated by fast Fourier transform to identify the energy peak position at each frequency. The phase velocity is calculated using the spacing between adjacent detectors. Finally, the Rayleigh wave observation dispersion curve is obtained by fitting the frequency as the horizontal axis and the phase velocity as the vertical axis.

[0066] In one embodiment of the present invention, the target area has a preset depth range of 0 to 120 meters, covering the needs of surface to shallow and middle-layer stratigraphic structure detection. A regular grid is used to discretize the space, with horizontal spacing of 25 meters and vertical layer thickness of 25 meters. Each grid cell corresponds to independent parameters to be inverted. This grid system transforms the continuous subsurface medium into discrete parameter units, facilitating subsequent layer-by-layer velocity model construction and geological structure interpretation based on seismic gathers, first wave arrival times, and Rayleigh wave dispersion curves.

[0067] In one embodiment of the present invention, loess layer thickness, P-wave propagation velocity and quality factor are set as parameters to be inverted for each grid cell; before inversion, initial values ​​are set for each grid cell, wherein the loess layer thickness is set to 25 meters according to the vertical grid layer thickness, the P-wave propagation velocity is preset based on regional experience, and the initial value of the quality factor is set to 50.

[0068] In one embodiment of the present invention, the step of estimating the vertical ray path time includes:

[0069] S201: Extract the grid cells that the seismic waves pass through along the vertical propagation path of the seismic waves and form a set of grid cells. That is, take the vertical path starting point directly below the source point as the starting point, and extract all the grid cells that the path passes through layer by layer according to the preset rule grid division to form a set of grid cells.

[0070] S202, calculate the ratio of loess layer thickness to P-wave propagation velocity for each grid cell to obtain the one-way travel time;

[0071] S203, sum the one-way travel time of all grid cells in the grid cell set to obtain the vertical ray path time.

[0072] The formula for calculating the vertical ray path time is as follows:

[0073] ;

[0074] in, The vertical ray path time represents the overall propagation time of a seismic wave along a vertical path through the target area. N represents the number of grid cells in the grid cell set, and i represents the grid cell index. This represents the thickness of the i-th mesh cell. This represents the P-wave propagation speed of the i-th grid cell.

[0075] In one embodiment of the present invention, the step of the multilayer dielectric dispersion forward modeling includes:

[0076] S301, for each frequency in the Rayleigh wave observation dispersion curve, construct the transfer matrix of each medium layer based on the Haskell-Thomson method; specifically, for the j-th medium layer, its transfer matrix is: ,in, Let f represent the transfer matrix of the j-th medium layer, and let f represent the frequency. This represents the phase velocity to be solved. Indicates the thickness of the j-th layer. This represents the density of the j-th layer. Indicates angular velocity;

[0077] S302, combine the transfer matrices of all layers into a total transfer matrix, which is the product of the transfer matrices of each layer;

[0078] S303, construct the determinant equation with the determinant of the total transfer matrix equal to zero, and solve it to obtain the phase velocity corresponding to the frequency.

[0079] In one embodiment of the present invention, the seismic traces recorded by each detector under the same excitation are collected, and the amplitude values ​​corresponding to the zero offset time are summed after squaring to obtain the post-stack energy; wherein, the formula for calculating the post-stack energy is: Where E represents the stacked energy, P represents the total number of detectors, and k represents the detector index. The amplitude value of the gather at the zero offset time (i.e., the source excitation time t=0) is represented; the post-stack energy enhances the strength of the effective signal and suppresses random noise by superimposing the energy of multiple channels in the same excitation, and is used to evaluate the stability of the source excitation energy and the overall response characteristics of the underground medium to seismic waves.

[0080] In one embodiment of the present invention, the step of constructing the multi-objective error function includes:

[0081] S401, for each observation pair consisting of the source location and the acquisition point, calculate the squared difference between the first arrival time of the acquired seismic wave and the path time of the vertical ray, and sum the squared differences of all observation pairs to obtain the sum of squared errors of the first wave arrival time;

[0082] S402, calculate the squared difference between the P-wave propagation velocity and the phase velocity at each frequency, and sum the squared differences for all frequencies to obtain the sum of squared Rayleigh wave dispersion errors;

[0083] S403 linearly combines the sum of squared errors of the first wave arrival time, the sum of squared errors of the Rayleigh wave dispersion, and the post-stack energy to obtain a multi-objective error function.

[0084] The formula for calculating the multi-objective error function is as follows:

[0085] ;

[0086] Where F represents the error value of the multi-objective error function, , and These represent the first, second, and third weighting coefficients, respectively; Q represents the number of source excitations; and q represents the source excitation index. This represents the first arrival time of the seismic wave generated by the k-th detector and the q-th source. Let R represent the path time of the vertical ray excited by the k-th detector at the q-th source, where R represents the number of frequency points and r represents the frequency index. Represents the r-th frequency. This represents the propagation speed of the P-wave corresponding to the r-th frequency. This represents the phase velocity at the r-th frequency. This represents the sum of squared errors when the first wave arrives. This represents the sum of squares of Rayleigh wave dispersion error. The above multi-objective error function, through multi-dimensional error joint optimization, can simultaneously calibrate the first wave travel time, dispersion characteristics, and data energy, thereby improving the stability and accuracy of underground medium parameter inversion.

[0087] In one embodiment of the present invention, the parameters to be inverted for the grid cells within each subdomain are subjected to differential evolution iteration until local convergence is achieved. Specific steps include:

[0088] S501, in each subdomain, multiple candidate individuals are randomly generated; each candidate individual is composed of the loess layer thickness, P-wave propagation velocity and quality factor of all grid cells in that subdomain;

[0089] S502, For any candidate individual, generate a mutation vector based on the individual in the subdomain that is closest to the actual measured thickness and three random candidate individuals;

[0090] Specifically, the steps for generating mutation vectors include:

[0091] S601, three different candidate individuals are randomly selected, denoted as random individual one, random individual two and random individual three, and the candidate individual with the smallest deviation from the measured thickness is selected from all candidate individuals in the subdomain, denoted as the individual with the best thickness fit;

[0092] S602, using all the parameters to be inverted of random individual one as the initial reference, calculate the difference between the thickness fit optimal individual and the parameters to be inverted of random individual one in each grid cell, multiply the difference by the first preset scaling factor and combine it with the initial reference;

[0093] S603, calculate the difference between the parameters to be inverted for random individual two and random individual three in each grid cell, multiply the difference by the second preset scaling factor, and combine it with the result of S602;

[0094] The result of S604 and S603 is accumulated to obtain the mutation vector.

[0095] The formula for calculating the mutation vector is as follows: Where V represents the mutation vector. This represents the individual with the optimal thickness fit. , and These represent random individual one, random individual two, and random individual three, respectively. By introducing optimal individual information and random perturbation, the global search and local development capabilities are balanced.

[0096] S503, dynamically adjust the crossover probability according to the average quality factor of all candidate individuals in the subdomain, and cross the mutation vector with the original candidate individuals according to the probability in the corresponding grid cell parameter dimension to generate experimental candidate individuals. Specifically, the higher the average quality factor, the smaller the crossover probability, so as to retain high-quality genes.

[0097] S504, for each experimental candidate individual, calculate the sum of squares of the first wave arrival error and the sum of squares of the Rayleigh wave dispersion error, and perform a weighted summation to obtain the local fitness; if the local fitness of the experimental candidate individual is better than that of the original candidate individual, then replace it; otherwise, retain the original candidate; wherein, the formula for calculating the local fitness is:

[0098] ;

[0099] in, Indicates local fitness. and These represent the fourth and fifth weight coefficients, respectively, and their sum is 1. By employing a local fitness comparison mechanism, we ensure that each iteration evolves towards a better solution, thereby improving the accuracy of parameter inversion.

[0100] S505: After the first preset number of iterations, a certain proportion of the best candidate individuals from each subdomain are selected and placed into the global elite library. For each elite individual in the elite library, if its Euclidean distance with any candidate individual in any subdomain is lower than the first preset threshold, then the elite individual replaces the candidate individual with the largest Euclidean distance in that subdomain. By constructing a global elite library and implementing cross-subdomain individual replacement, the local search limitation is broken, information sharing between subdomains is promoted, and the replacement strategy is controlled by the Euclidean distance threshold, which effectively avoids the algorithm from getting stuck in local optima and accelerates global convergence.

[0101] S506: When the change in local fitness of the best candidate individual in a subdomain is lower than the second preset threshold and the difference in loess layer thickness between adjacent grid cells in the subdomain is lower than the third preset threshold, the subdomain is determined to be locally converged, and the differential evolution in the subdomain ends. By using the dual convergence criteria of local fitness change and difference in thickness between adjacent grid cells, the stability of the solution is guaranteed, the spatial continuity is constrained, the geological characteristics of the loess layer are consistent, overfitting is prevented, and the geological rationality of the inversion results is ensured.

[0102] In one embodiment of the present invention, step S106 specifically includes:

[0103] S701, establish multiple parallel annealing chains within the target area, each annealing chain corresponding to a candidate individual, i.e. a set of parameters to be inverted;

[0104] S702, based on the loess layer moisture content collected on site, sets the initial temperature of the first annealing chain, and the temperature of subsequent annealing chains decreases by a fixed proportion. The initial temperature is achieved through a linear mapping relationship, and the temperature of subsequent annealing chains decreases by a fixed proportion, forming a temperature gradient sequence. This setting makes the initial temperature of high moisture content areas higher, enhancing the algorithm's ability to explore complex geological conditions.

[0105] S703: For each candidate individual corresponding to an annealing chain, a perturbation vector is generated by Lévy flight and superimposed on the candidate individual to obtain a mutated candidate individual. By using Lévy flight to generate perturbation vectors, it is possible to effectively escape local optima.

[0106] S704, for the quality factor of each grid cell in the mutant candidate individuals, the viscoelastic dissipation energy increment is calculated based on the loess layer thickness and P-wave propagation velocity, and the quality factor is corrected accordingly; wherein, the formula for calculating the viscoelastic dissipation energy increment is: ,in, The value represents the increase in viscoelastic dissipated energy, and f represents the frequency. This represents the loess layer thickness in the i-th grid cell. Let A represent the quality factor of the i-th grid cell, and let A represent the seismic wave amplitude. Let represent the P-wave propagation velocity of the i-th grid cell. The formula for calculating the corrected quality factor is: ,in, This represents the corrected quality factor for the i-th grid cell. This represents the quality factor of the i-th grid cell before correction. Indicates reference energy; by suppressing non-physical high-frequency gain caused by excessively low quality factor, it conforms to the viscoelastic attenuation characteristics of loess layer;

[0107] S705. Compare the error values ​​of the global multi-objective error function between the mutated candidate individual and the original candidate individual. If the error value of the mutated candidate individual is not higher than the error value of the original candidate individual, the mutated candidate individual is directly accepted; otherwise, it is not accepted.

[0108] S706, every second preset iteration, attempt to exchange candidate individuals between two adjacent annealing chains. If the exchange can reduce the sum of errors of the two chains, then exchange; otherwise, perform a forced exchange with the first preset probability to maintain population diversity. When the error change of all annealing chains is lower than the fourth preset threshold in multiple consecutive rounds, output the optimal parameters to be inverted.

[0109] The above steps, through the temperature gradient setting of the parallel annealing chain, combined with the Lévy flight perturbation and physical constraints of the quality factor, effectively balance global search and local exploitation capabilities, making it particularly suitable for parameter inversion in the highly heterogeneous media of the Loess Plateau region. The cross-chain individual exchange mechanism further enhances information flow among the populations, preventing the algorithm from getting trapped in local optima. The viscoelastic dissipative energy correction strategy ensures the physical rationality of the quality factor and improves the fitting accuracy of the inversion results to high-frequency seismic signals.

[0110] In one embodiment of the present invention, the static correction time difference is calculated based on the P-wave propagation velocity and loess layer thickness of each grid cell in the optimal parameters to be inverted; the absorption compensation coefficient is determined based on the attenuation relationship between the quality factor and the frequency; wherein, static correction is used to adjust the time axis of the seismic gather by accumulating the time difference on the path between each source and the receiver, and absorption compensation is achieved by applying an inverse attenuation function to the spectrum of the seismic gather.

[0111] Specifically, the formula for calculating the static correction time difference is:

[0112] ;

[0113] in, Indicates static correction time difference. Indicates one-way travel time. Indicates the thickness of the reference surface. The reference velocity is represented by z, and the elevation difference between the sampling point and the reference surface is represented by z. It represents the average velocity of the strata above the reference plane; by accumulating the static correction time difference to the sampling time of the corresponding seismic trace, the gap is filled by linear interpolation to achieve the alignment correction of the time axis;

[0114] The formula for calculating the absorption compensation coefficient is: ,in, The coefficient of absorption compensation is represented by exp, which represents the exponential function. denoted by v, the horizontal distance from the source to the receiver is given; v represents the path average velocity; Q represents the quality factor; and f represents the frequency. The seismic gathers are converted to the frequency domain using a fast Fourier transform to obtain the amplitude spectrum. Then, each frequency component is multiplied by an absorption compensation coefficient. Finally, the inverse fast Fourier transform is performed back to the time domain to complete the inverse compensation for high-frequency attenuation.

[0115] Static correction, by accurately calculating travel time deviations caused by near-surface heterogeneity, effectively eliminates the interference of surface undulations and low-velocity zones on seismic wave arrival times, ensuring accurate alignment of phase axes within the gather and improving the quality of subsequent stacked imaging. Absorption compensation utilizes the quality factor obtained from inversion to construct a frequency domain gain function, specifically recovering high-frequency energy losses caused by viscoelastic attenuation in the loess layer, broadening the effective frequency band of the seismic signal, and improving profile resolution and the accuracy of geological interpretation. The combination of these two methods forms a complete closed loop from parameter inversion to data correction, significantly enhancing the processing accuracy of 3D seismic exploration data for coalfields in the thick loess plateau region.

[0116] It should be noted that the interval and threshold sizes are set for ease of comparison. The size of the threshold depends on the amount of sample data and the base number set by those skilled in the art for each set of sample data, as long as it does not affect the proportional relationship between the parameter and the quantized value. Furthermore, the above formulas are all dimensionless calculations, and the formulas are derived from software simulations using a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0117] The embodiments of the present invention have been described above, but the present invention is not limited to the specific embodiments described above. The specific embodiments described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms based on the guidance of the present embodiments, all of which are within the protection scope of the present embodiments.

Claims

1. A three-dimensional seismic exploration method for coalfields in thick loess plateau areas, characterized in that, Includes the following steps: S101, collects artificial source shot-shoal records and shallow Rayleigh wave detector data, and extracts seismic gathers, first wave arrival time observations and Rayleigh wave dispersion curves at different frequencies corresponding to each source location. S102, the target area is divided into regular grids according to the preset depth range, and the loess layer thickness, P-wave propagation velocity and quality factor are set as parameters to be inverted for each grid cell; S103, based on the parameters to be inverted for each grid cell, vertical ray path time estimation, multilayer medium dispersion forward modeling, and post-stack energy estimation are performed respectively to obtain vertical ray path time, phase velocity, and post-stack energy; S104, based on the first wave arrival time observation, Rayleigh wave observation dispersion curve, theoretical ray path time, phase velocity and post-stack energy, a multi-objective error function is generated by weighting according to preset weights; S105, based on the thickness of the loess layer, divides the target area into three sub-domains: low, medium, and high. Differential evolution iteration is performed on the parameters to be inverted for the grid cells in each sub-domain until local convergence is achieved. S106, establish multiple parallel annealing chains, set the chain temperature based on the loess layer moisture content collected in the field and decrease it chain by chain, optimize the parameters to be inverted through jump variation, energy dissipation constraint correction of quality factor, and inter-chain state exchange mechanism until convergence and output the optimal parameters to be inverted. S107: Calculate the static correction time difference of each acquisition point based on the optimal parameters to be inverted, calculate the absorption compensation coefficient based on the parameters to be inverted, and perform static correction and absorption compensation on the seismic gathers.

2. The method for three-dimensional seismic exploration of coalfields in thick loess plateau areas according to claim 1, characterized in that, The artificial seismic source shot record is the raw amplitude time series data recorded on the ground geophone; The seismic gather is an amplitude sequence continuously recorded by the detector at a preset sampling frequency during a single source excitation process; The first wave arrival time observation value was extracted by using a threshold discrimination method on the amplitude signal of the seismic trace collection; The Rayleigh wave observation dispersion curve is obtained by performing time-spectrum analysis on multiple channels of the same excitation to identify the phase velocity corresponding to the peak value of the Rayleigh wave energy concentration. The window function used in the time-spectrum analysis is the Hanning window.

3. The three-dimensional seismic exploration method for coalfields in thick loess plateau areas according to claim 1, characterized in that, The steps for estimating the vertical ray path time include: S201: Extract the grid cells through which the seismic waves pass along the vertical propagation path of the seismic waves and form a set of grid cells. S202, calculate the ratio of loess layer thickness to P-wave propagation velocity for each grid cell to obtain the one-way travel time; S203, sum the one-way travel time of all grid cells in the grid cell set to obtain the vertical ray path time.

4. The three-dimensional seismic exploration method for coalfields in thick loess plateau areas according to claim 1, characterized in that, The steps of the multilayer medium dispersion forward modeling include: S301, for each frequency in the Rayleigh wave observation dispersion curve, construct the transfer matrix of each medium layer based on the Haskell-Thomson method; S302, combine the transfer matrices of all layers into a total transfer matrix; S303, construct the determinant equation with the determinant of the total transfer matrix equal to zero, and solve it to obtain the phase velocity corresponding to the frequency.

5. The three-dimensional seismic exploration method for coalfields in thick loess plateau areas according to claim 1, characterized in that, The seismic traces recorded by each geophone under the same excitation are collected, and the amplitude values ​​corresponding to the zero offset time are summed after squaring to obtain the post-stack energy.

6. The method for three-dimensional seismic exploration of coalfields in thick loess plateau areas according to claim 3, characterized in that, The steps for constructing the multi-objective error function include: S401, for each observation pair consisting of the source location and the acquisition point, calculate the squared difference between the first arrival time of the acquired seismic wave and the path time of the vertical ray, and sum the squared differences of all observation pairs to obtain the sum of squared errors of the first wave arrival time; S402, calculate the squared difference between the P-wave propagation velocity and the phase velocity at each frequency, and sum the squared differences for all frequencies to obtain the sum of squared Rayleigh wave dispersion errors; S403 linearly combines the sum of squared errors of the first wave arrival time, the sum of squared errors of the Rayleigh wave dispersion, and the post-stack energy to obtain a multi-objective error function.

7. The method for three-dimensional seismic exploration of coalfields in thick loess plateau areas according to claim 1, characterized in that, Within each subdomain, the parameters to be inverted for the grid cells are subjected to differential evolution iterations until local convergence. The specific steps include: S501, in each subdomain, multiple candidate individuals are randomly generated; each candidate individual is composed of the loess layer thickness, P-wave propagation velocity and quality factor of all grid cells in that subdomain; S502, For any candidate individual, generate a mutation vector based on the individual in the subdomain that is closest to the actual measured thickness and three random candidate individuals; S503: Dynamically adjust the crossover probability based on the average quality factor of all candidate individuals in the subdomain, and cross the mutation vector with the original candidate individuals according to the probability in the corresponding grid cell parameter dimension to generate experimental candidate individuals; S504. For each experimental candidate individual, calculate the sum of squares of the first wave arrival error and the sum of squares of the Rayleigh wave dispersion error, and perform a weighted summation to obtain the local fitness. If the local fitness of the experimental candidate individual is better than that of the original candidate individual, then replace it; otherwise, retain the original candidate. S505, after the first preset number of iterations, select a certain proportion of the best candidate individuals from each subdomain and put them into the global elite library; for each elite individual in the elite library, if its Euclidean distance with any candidate individual in any subdomain is lower than the first preset threshold, then replace the candidate individual with the largest Euclidean distance in that subdomain with that elite individual. S506: When the change in the local fitness of the best candidate individual in a subdomain is lower than the second preset threshold and the difference in loess layer thickness between adjacent grid cells in the subdomain is lower than the third preset threshold, the subdomain is determined to be locally converged, and the differential evolution in the subdomain ends.

8. The three-dimensional seismic exploration method for coalfields in thick loess plateau areas according to claim 7, characterized in that, The steps for generating mutation vectors include: S601, three different candidate individuals are randomly selected, denoted as random individual one, random individual two and random individual three, and the candidate individual with the smallest deviation from the measured thickness is selected from all candidate individuals in the subdomain, denoted as the individual with the best thickness fit; S602, using all the parameters to be inverted of random individual one as the initial reference, calculate the difference between the thickness fit optimal individual and the parameters to be inverted of random individual one in each grid cell, multiply the difference by the preset scaling factor and combine it with the initial reference; S603, calculate the difference between the parameters to be inverted for random individual two and random individual three in each grid cell, multiply the difference by the preset scaling factor, and combine it with the result of S602; S604 is the result of summing the results of S603 to obtain the mutation vector.

9. A three-dimensional seismic exploration method for coalfields in thick loess plateau areas according to claim 1, characterized in that, The specific steps in S106 include: S701, establish multiple parallel annealing chains within the target area, with each annealing chain corresponding to a candidate individual; S702, based on the moisture content of the loess layer collected on site, set the initial temperature of the first annealing chain, and the temperature of subsequent annealing chains decreases by a fixed ratio; S703, for each candidate individual corresponding to an annealing chain, a perturbation vector is generated using Levy flight and superimposed on the candidate individual to obtain a mutated candidate individual; S704, for the quality factor of each grid cell in the mutant candidate individuals, calculate the viscoelastic dissipation energy increment based on the loess layer thickness and P wave propagation velocity, and correct the quality factor accordingly. S705, compare the error values ​​of the global multi-objective error function of the mutated candidate individual with those of the original candidate individual. If the error value of the mutated candidate individual is not higher than that of the original candidate individual, then the mutated candidate individual is directly accepted. S706, every second preset iteration number, attempt to swap candidate individuals between two adjacent annealing chains. If the swap reduces the sum of errors of the two chains, then swap; otherwise, perform a forced swap with the first preset probability. When the error change of all annealing chains is lower than the fourth preset threshold in multiple consecutive rounds, output the optimal parameters to be inverted.

10. A three-dimensional seismic exploration method for coalfields in thick loess plateau areas according to claim 1, characterized in that, The static correction time difference is calculated based on the P-wave propagation velocity and loess layer thickness of each grid cell in the optimal parameters to be inverted; the absorption compensation coefficient is determined based on the attenuation relationship between the quality factor and the frequency; wherein, static correction is used to adjust the time axis of the seismic gather by accumulating the time difference on the path between each source and the receiver, and absorption compensation is achieved by applying an inverse attenuation function to the spectrum of the seismic gather.

Images