An adaptive frequency band extension method for non-stationary post-stack seismic data

Abstract

The application discloses a kind of adaptive band extension methods for non-stationary post-stack seismic data, it is related to the technical field of seismic data processing, after band-limited profile normalization, construct original amplitude, directional gradient, Laplace response and directional local details, input pre-trained multi-scale residual encoder-decoder;Characteristic is extracted by anisotropic convolution, attention, vertical Fourier hybrid and sample level and local level spectrum modulation, predict full-band correction residual, add with input and apply adaptive amplitude constraint, output same size wideband profile;Training sample is generated by the same reflectivity model and wide, narrow frequency band space-time variation wave field pairing, and the training target jointly uses mask Charbonnier reconstruction loss, vertical logarithmic amplitude spectrum loss and common band complex spectrum loss.The method is used to compensate low, high frequency end component and keep reliable public band amplitude correlation.

Description

Technical Field

[0001] This invention relates to the field of seismic data processing technology, and in particular to an adaptive frequency band extension method for non-stationary post-stack seismic data. Background Technology

[0002] The effective frequency band of back-stack seismic profiles is limited by factors such as seismic wavelets, formation absorption, noise, and processing procedures, often manifesting as insufficient low-frequency information, high-frequency energy attenuation, and local waveform distortion. Frequency band limitation reduces the reliability of thin-layer resolution, fault characterization, and reservoir interpretation. Real seismic data also exhibits significant non-stationarity; wavelets and degradation levels can vary with time, depth, lateral location, or both.

[0003] Existing methods such as deconvolution, inverse Q-filtering, and sparse reflectivity inversion typically rely on wavelet estimation, sparsity assumptions, or high signal-to-noise ratios. Under complex geological structures, lateral heterogeneity, and strong noise conditions, non-uniqueness, noise amplification, or spurious frequency components may occur. While existing deep learning methods can recover details through encoder-decoder approaches, if seismic profiles are processed as ordinary image enhancement methods, they easily confuse vertical waveform frequency variations with lateral geological structural continuity, and it is also difficult to simultaneously address sideband compensation, common band preservation, and local spectral non-stationarity.

[0004] Therefore, there is a need for a method that can adaptively adjust the compensation intensity for spectral differences between samples and within the profile without changing the spatial dimensions of the post-stack seismic profile, while simultaneously constraining the band extension of waveform, amplitude spectrum and common band complex spectrum. Summary of the Invention

[0005] To address the technical problems of simultaneous loss or attenuation of low and high frequencies in band-limited post-stack seismic data, difficulty in adapting to non-stationary spectral degradation, easy disruption of common frequency band amplitude-phase correlation, and easy generation of uncontrolled high-frequency textures, this invention discloses an adaptive frequency band extension method for non-stationary post-stack seismic data.

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

[0007] An adaptive bandwidth spreading method for non-stationary post-stack seismic data includes the following steps:

[0008] S1, obtain the band-limited post-stack seismic profile to be processed, normalize the amplitude of the band-limited post-stack seismic profile while keeping its time sampling number and trace number unchanged;

[0009] S2, construct a priori seismic features based on the zone-limited post-stack seismic profile, and project the a priori seismic features as initial features;

[0010] S3 uses a pre-trained multi-scale residual encoder-decoder to process initial features, and uses adaptive frequency band expansion residual blocks of different scales to extract vertical waveform, lateral in-phase axis continuity and frequency domain interaction information. It also fuses encoded and decoded features through skip connections.

[0011] S4. Inside the residual block in step S3, sample-level channel modulation weights are generated based on the vertical amplitude spectrum of the current feature, and position-related modulation weights are generated based on the local energy of the sample-level modulation feature. The frequency band compensation intensity is adjusted in two stages, and the modulated feature is returned to the residual block for further processing.

[0012] S5. Based on the final decoding features, predict the single-channel full-band correction residual, add it to the normalized band-limited post-stack seismic profile, and apply adaptive amplitude constraints based on the effective amplitude statistics or dynamic range of the input profile.

[0013] S6 restores the normalized scale of the output profile, outputting a broadband post-stack seismic enhancement profile with the same size as the band-limited post-stack seismic profile.

[0014] Furthermore, in step S2, the a priori seismic features include six categories: the original amplitude, vertical gradient, lateral gradient, Laplace response, vertical local detail, and lateral local detail of the band-limited stacked seismic profile. These features are then stitched together and projected according to the following formula:

[0015] ;

[0016] In the formula, This is the normalized single-channel band-limited input profile; and These are the vertical gradient and the lateral gradient, respectively. , These represent the difference operators in the corresponding directions, with subscripts indicating the subscripts. Indicates the horizontal direction; It is a two-dimensional discrete Laplace response. It is the Laplace operator; ; These are local details in the vertical and horizontal directions, respectively. , These are the local averaging operators for the vertical and horizontal directions, respectively; Concat indicates concatenation along the channel dimension; It is a six-channel prior feature;

[0017] Boundary padding is used to maintain consistent dimensions across all prior channels, followed by shallow convolutional projection. We obtain the initial features, where, This is a shallow convolution projection. Output it.

[0018] Further, in step S3, the multi-scale residual encoder-decoder includes an encoding path that progressively reduces the spatial size and increases the number of channels, a bottleneck layer at the coarsest scale, and a decoding path that progressively restores the spatial size. The decoding path fuses the upsampled features with the same-scale features of the encoding path at each corresponding scale. Its pre-training process is completed before executing step S1, including: convolving the same reflectivity model with wide and narrow frequency band sub-wavefields respectively to form paired training samples; obtaining the prediction profile by using network forward computation; updating the network parameters by weighting and summing the mask Charbonnier reconstruction loss, vertical logarithmic amplitude spectrum loss, and common frequency band complex spectrum loss. For the network that has been trained, the trained parameters are reused when processing different profiles to be processed.

[0019] In step S3, the adaptive band-expanded residual block undergoes anisotropic seismic convolution, axial attention and window attention, vertical full-band Fourier mixing, local spectral modulation, and gated convolution feedforward processing sequentially, following the methods of normalization, feature transformation, learnable residual scaling, and residual summation; wherein:

[0020] Anisotropic convolutions are fused in parallel according to the following formula: local 2D, vertical long kernel, horizontal long kernel, dilated, and depth-separable convolution branches.

[0021] ;

[0022] In the formula, This is the input feature tensor for the current processing stage; This is an anisotropic convolutional fusion output; , , These are two-dimensional convolution operators with kernel sizes of 3×3, 7×1, and 1×7, respectively, where the first dimension of the kernel size is vertical and the second dimension is horizontal. For dilated convolution operators; It is a depthwise separable convolution operator; This is a fusion projection performed after stitching together the branches along the channel dimension;

[0023] Axial attention establishes long-range dependencies along both the vertical and horizontal axes, while window attention aggregates structural information within a local two-dimensional window. Vertical Fourier mixing is implemented according to the following formula:

[0024] ;

[0025] In the formula, Indicates the vertical sampling direction. This is a one-dimensional Fourier transform along the vertical sampling axis of the feature tensor, while other dimensions remain unchanged; S is the vertical complex spectrum obtained by the transform;

[0026] ;

[0027] In the formula, The frequency domain mixing characteristics after inverse transformation; , These are the real and imaginary parts of the complex spectrum S, respectively; A mapping that allows for the learnable mixing of the real and imaginary parts and forms the output complex spectrum; This is the inverse Fourier transform along the vertical sampling axis;

[0028] Gated convolution feedforward processing is performed according to the following formula:

[0029] ;

[0030] In the formula, , , These are gated branch mapping, value branch mapping, and output projection, respectively. It represents element-wise multiplication between homogeneous tensors, broadcasting by channel or spatial location when necessary; For the output of the gated feedforward branch, the subscripts of each mapping name , , Identify the gating, value, and output respectively.

[0031] Further, step S4 includes generating channel weights based on the characteristic vertical amplitude spectrum and performing sample-level modulation, followed by local-level modulation based on local energy; the sample-level and local-level modulation execution processes are as follows:

[0032] ;

[0033] In the formula, Its vertical amplitude spectrum; For spectral pooling operators, statistical information in the frequency dimension and lateral position is pooled to describe the channel of each sample; MLP stands for channel multilayer perceptron. It is a sigmoid activation function. This is the independent variable of the function; The modulation weights for each sample and each channel are shared at the spatial location of that channel;

[0034] ;

[0035] In the formula, Features after sample-level modulation;

[0036] because Single-level modulation factor Located at (0.5, 1.5), centered at 1;

[0037] ;

[0038] ;

[0039] In the formula, From Extracted local energy description; A convolutional mapping that maps energy descriptions to location-dependent responses; For local modulation weights; The features are those after local modulation, with the subscript 1 indicating the local level;

[0040] Single-level factor Located at (0.5, 1.5), the overall scaling range of the two-level multiplication is (0.25, 2.25).

[0041] Further, in step S5, the residual prediction head output is used to compensate for the full-band correction residual of the missing or attenuated components at the low-frequency and high-frequency ends, and is added to the input according to the following formula:

[0042] ;

[0043] In the formula, The residual prediction head outputs a single-channel full-band correction residual based on the final decoded features; This is the broadband prediction profile after summing the residuals; all three have the same number of time samples and channels; amplitude constraints apply to... The constrained profile is then normalized and scaled to obtain the final output.

[0044] right When applying amplitude constraints, the lower limit is determined based on the input effective amplitude quantiles or dynamic range. and upper limit ,and Each predicted amplitude is truncated to this interval; , These represent the adaptive lower and upper limits, respectively. The amplitude statistics only use the effective region, and the constraint rules remain consistent during training and inference.

[0045] The beneficial effects of this invention are:

[0046] First, by using a six-channel seismic prior, directional convolution, and long-range attention in a coordinated process, the vertical waveform frequency information and the lateral structural continuity can be distinguished and modeled, which helps to maintain the reflection position, polarity, and lateral morphology when increasing the resolution.

[0047] Second, by using vertical full-band Fourier mixing and sample-level and local-level dual-scale spectral modulation, the compensation intensity is adjusted according to the overall bandwidth state and local energy distribution, thereby adapting to non-stationary wavelets and attenuation that change with time and lateral position.

[0048] Third, by using full-band residual learning and common-band complex spectrum loss, the amplitude and phase of the reliable common-band are jointly constrained while recovering the low-frequency and high-frequency components, and the risk of overcompensation is reduced by using adaptive amplitude constraints.

[0049] Fourth, by pairing wide and narrow band non-stationary convolutions under the same reflectivity model, as well as absorption and noise degradation, the coverage of training samples in terms of construction patterns, wavelets, and degradation modes is improved.

[0050] In summary, ablation and independent synthesis testing support the complementary role of the above-mentioned techniques; however, the field migration effect should still be used in conjunction with residual and amplitude consistency checks. Attached Figure Description

[0051] Figure 1 Prepare a relationship diagram for the overall process and training of this method;

[0052] Figure 2 This is a flowchart illustrating the process of constructing paired non-stationary training samples during the network training preparation in this embodiment of the invention.

[0053] Figure 3 This is a diagram showing the overall network structure and corresponding steps in an embodiment of the present invention;

[0054] Figure 4 This is a structural diagram of the adaptive bandwidth extension residual block in an embodiment of the present invention;

[0055] Figure 5 This is the non-stationary wavelet band inflection point parameter field for narrow-frequency input and wide-frequency target in this embodiment of the invention, where (a) and (b) are the low-frequency passband inflection points for narrow-frequency input and wide-frequency target, respectively. The parameter fields (c) and (d) represent the high-frequency passband inflection points of the two parameters, respectively. Parameter field;

[0056] Figure 6 These are the three loss convergence curves in this embodiment of the invention;

[0057] Figure 7 This is the band-normalized amplitude spectrum ratio curve of the processed independent synthetic test plot subset in this embodiment of the invention.

[0058] Figure 8 The images show a comparison of the Marmousi2 model profiles before and after processing in this embodiment of the invention. (a) is the band-limited input profile, (b) is the broadband enhanced profile after processing by the method of this invention, and (c) and (d) are magnified views of the two profiles in the fault-dense region (traces 900-2000 and 700-1500ms). Detailed Implementation

[0059] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0060] This invention discloses an adaptive bandwidth extension method for non-stationary post-stack seismic data. For example... Figure 1 As shown, this embodiment retains the main steps S1 to S6. First, the network training preparation process, which is used in step S3 and is actually executed before step S1, is introduced. Then, S1, S2, S3 and their internal steps S4, S5 and S6 are explained in sequence. Finally, experimental verification is given.

[0061] The pre-training process of the network used in this method is as follows:

[0062] a. Generate a reflectivity model;

[0063] like Figure 2 As shown, a two-dimensional reflectivity model is generated that includes layered interfaces, faults, folds, pinch-outs, and local anomalies. Each input and target pair shares the same model to ensure that the reflection position, polarity, and structural geometry correspond to each other. Training, validation, and testing are divided into units based on the reflectivity model.

[0064] b. Construct wide-band and narrow-band sub-wavefields;

[0065] A broadband target wavelet field and a narrow-band input wavelet field are constructed for the same model, ensuring that their phase and spatial variation trends correspond to each other, while their effective frequency band endpoints differ. The wavelets can adopt the Ricker, Ormsby, or Klauder family of wavelets, and their variation modes include stationary, time-varying, transverse position-varying, and spatiotemporally combined variations. The unified spectral representation of each local wavelet is as follows:

[0066] (1);

[0067] In the formula, For position The complex spectrum of a local wavelet; This is the non-negative amplitude spectrum of the wavelet; The phase spectrum is expressed in radians. Frequency, in Hz; To meet The imaginary unit; exp is the exponential function, and the semicolon is used to distinguish frequency variables. The position on which the wavelet parameters depend Imaginary unit here. With the discrete channel index below Distinguish them according to the context of their respective formulas.

[0068] When using Ricker wavelets, equation (2) is used to generate the amplitude spectrum for local dominant frequency control, and equation (3) is used to set the corresponding phase and local time shift. Substituting these values ​​into equation (1):

[0069] (2);

[0070] In the formula, The normalized amplitude spectrum of the Ricker wavelet; For position The local dominant frequency at a given location, measured in Hz; The subscript represents the square of the local main frequency. Indicates the main frequency control type.

[0071] The amplitude spectrum is in The value reaches its maximum value of 1 at that point.

[0072] (3);

[0073] In the formula, The phase spectrum of the main frequency control wavelet; This represents the local phase rotation, expressed in radians. This represents the local wavelet time shift, in seconds. Pi; This is the linear phase term caused by time shift.

[0074] This formula is used for the non-negative frequency end, while the negative frequency end is filled with conjugate symmetry to generate a real-valued wavelet; the band control wavelet also uses the same phase construction. and When all phases are zero, a zero-phase wavelet is obtained.

[0075] When using a frequency band-controlled wavelet, equation (4) specifies the four inflection points of the local frequency band. The low-frequency transition band rises smoothly, the main passband maintains unit amplitude, the high-frequency transition band falls smoothly, and the passband is cut off; the transition band can use a cosine gradient. The phase is still constructed according to equation (3), and this parameterization is used to directly control the effective frequency band of the paired wavelet:

[0076] (4);

[0077] In the formula, This is a local frequency band parameter vector; , , , In order, they are the low-frequency cutoff inflection point, the low-frequency passband inflection point, the high-frequency passband inflection point, and the high-frequency cutoff inflection point; the inequality in to Indicates the same The corresponding value at that location. to This is a low-frequency transition band. to Main passband, to It is a high-frequency transition band; For the Nyquist frequency, This refers to the time sampling interval. All frequencies mentioned above are in Hz and are called the inflection point frequencies of the frequency band, as opposed to angular frequencies, which are measured in rad / s.

[0078] The above-mentioned main frequency, frequency band inflection point, phase, and time shift are all generated by interpolation and smoothing of the control point values ​​according to equation (5) to avoid non-physical abrupt changes in adjacent channels and adjacent sampling positions:

[0079] (5);

[0080] In the formula, The continuous parameter field to be generated can be , any component, or ; Control points The parameter values ​​at the specified locations are i and j, which are the time direction and lateral control point indices, respectively. It is the set of all control point values; Indicates in Two-dimensional interpolation and smoothing operators for position evaluation.

[0081] After generating parameter fields for the four components of β respectively, the frequency order of equation (4) is guaranteed to hold.

[0082] Taking the Ormsby wavelet as an example, Figure 5 In the figure, (a) and (b) represent the low-frequency passband inflection points of the narrow-band input and the wide-band target, respectively. Parameter field, Figure 5 In the diagram, (c) and (d) represent the high-frequency passband inflection points of the two components, respectively. Parameter fields. Each parameter field changes smoothly with time and channel number, reflecting the spatiotemporal nonstationarity of the wavelet band; compared to narrow-band input, the wide-band target... reduce, The parameters are increased, thus extending the main passband to both low and high frequencies. Applying these parameter fields to a non-stationary convolution of the same reflectivity model generates input and target samples with corresponding locations but different frequency bands.

[0083] This embodiment can employ a generation configuration that extends the target frequency by approximately 3-8 Hz at the low end and approximately 12-18 Hz at the high end relative to the input. The specific endpoints are subject to the metadata of each sample. This configuration is not intended to replace... Figure 5 The actual sample shown , Parameter field; Figure 5 The difference in the main passband parameters at both ends of the graph can be directly read from this graph, and the upper frequency limit is always limited by the Nyquist frequency.

[0084] c, generates paired non-stationary convolutional responses;

[0085] Substituting the wideband and narrowband sub-wavefields generated in b into equation (6), and performing a profile-level nonstationary convolution with the same reflectivity model in a, we obtain the wideband target response and the narrowband response without additional degradation:

[0086] (6);

[0087] In the formula, This is the synthetic seismic record at the horizontal coordinate x and the two-way travel time t; Let τ be the reflectivity of the channel at the integration time. Centered at ,by For local wavelets relative to time; For superimposed noise; τ is the convolution integral variable. Let x be the integral time element. The units of t and τ are seconds, and x is the horizontal position.

[0088] Constructing a noise-free target and a reference timing before adding noise The noise is only superimposed once in the subsequent degradation stage. This stage makes... The target response serves as the basis for broadband targets, while the narrowband response is denoted as... This provides a band-limited reference and a noisy input for d generation.

[0089] d, apply amplitude, absorption, and noise degradation;

[0090] Narrowband response By applying amplitude modulation and frequency-dependent absorption according to equation (7), the quality factor is allowed to vary with time and lateral position:

[0091] (7);

[0092] In the formula, This is the narrowband response without additional absorption attenuation obtained by equation (6) at n=0; This represents a complete time series with a fixed horizontal position x, where the dot “·” represents the time independent variable to be transformed. The recording before noise is applied after amplitude modulation and absorption. It is the amplitude modulation factor; The quality factor for spatiotemporal changes; The two-way propagation time used at this sampling location, in seconds; and These are the Fourier transform and inverse transform along the time direction, respectively, ending with... This indicates that the value is evaluated at output time t; This is a frequency-dependent attenuation factor.

[0093] Fix its at each output position and Then, an inverse transform is performed to evaluate the expression, which is the operator expression for non-stationary absorption simulation, and does not mean that a fixed-frequency filter is used only once for the entire profile.

[0094] Equation (8) specifically generates the equation in equation (7). Used to simulate smooth local amplitude changes:

[0095] (8);

[0096] In the formula, The positive amplitude modulation factor in equation (7); For logarithmic magnitude field; Control points The logarithmic magnitude at position (). This is the logarithmic magnitude field after interpolation and smoothing.

[0097] Exponential mapping guarantees This factor is used to simulate smooth local energy changes.

[0098] The output of equation (7) is retained as the band-limited reference before noise addition, and then Gaussian white noise, migration profile noise, or non-uniform noise with lateral weights is superimposed to form the band-limited input. The reference already includes the amplitude modulation and absorption attenuation used in the input, but does not include the noise added later; the noise term of equation (6) is not reused.

[0099] e: Forming training, validation, and test samples;

[0100] For band-limited inputs, wideband targets, and pre-noise band-limited references, consistent clipping, translation, and lateral flipping are applied, and normalization is performed using the same paired amplitude scale, denoted as X, Y, and , respectively. Meanwhile, the effective area mask M is preserved. Geometric transformations, amplitude scales, and sampling intervals are kept consistent. Data subsets are divided according to reflectivity models to avoid clipping blocks of the same model from falling into different subsets.

[0101] According to the original experimental configuration, the training set, validation set, and test set contain 7000, 400, and 1000 pairs of profiles, respectively; the profile size is 512×512, the time sampling interval is 4ms, corresponding to a Nyquist frequency of 125Hz. Four types of samples—stationary, time-varying, laterally varying, and spatiotemporally combined varying—each account for 25%. For noisy samples, the target signal-to-noise ratio is 8dB, with migration profile noise accounting for 20% and Gaussian white noise accounting for 80%. Pruning, translation, and normalization are performed offline during construction, and random lateral flipping is possible during training. The above values ​​are for the example configuration.

[0102] f, Forward computation and multi-domain joint optimization;

[0103] In each iteration, the training input X undergoes the same forward computations S2, S3 (inclusive) and S5 to obtain the prediction. The three loss branches are calculated as follows: the waveform branch first calculates equation (9), the amplitude spectrum branch calculates equations (10) and (11) in sequence, and the common frequency band branch calculates equations (12) to (15) in sequence; finally, the network parameters are updated using equations (16) and (17).

[0104] First, within the effective region, compare the predictions according to equation (9). Point-by-point waveform amplitude of broadband target Y:

[0105] (9);

[0106] In the formula, The mask Charbonnier reconstruction loss, also known as the effective region reconstruction loss; b, i, and j are the in-batch sample index, time sampling index, and channel index, respectively; and These are the corresponding amplitudes of the predicted profile and the broadband target profile, respectively; For the effective area mask, a value of 1 indicates effective sampling, and a value of 0 indicates non-participation in reconstruction loss; the denominator is the total number of effective sampling points within the batch and should be greater than 0; >0 is the smoothing stability constant, which is taken as 10⁻³ in this embodiment.

[0107] Secondly, perform orthogonal normalized real Fourier transform on each channel of the prediction and target according to equation (10) to obtain the input spectrum required by equation (11). And S, then calculate the vertical logarithmic amplitude spectrum loss:

[0108] (10);

[0109] In the formula, S and S are the predicted profiles after masking. The vertical complex spectrum relative to the target profile Y; M is the effective area mask with the same shape as the profile; For each channel, an orthogonal normalized one-dimensional real Fourier transform is performed, preserving the non-negative frequency coefficients; This is an element-wise multiplication. This formula generates two sets of input spectra of formula (11), and formula (13) uses the same transformation and normalization rules.

[0110] (11);

[0111] In the formula, For vertical logarithmic amplitude spectrum loss; , , respectively, are the predicted and target complex spectrum coefficients of the b-th sample, the k-th non-negative frequency, and the j-th channel in equation (10); k is the non-negative frequency index; the absolute value of the inner layer is the complex modulus, and the absolute value of the outermost layer is the absolute value of the difference between real numbers; log is the natural logarithm; This represents the total number of spectral coefficients involved in the summation.

[0112] When the batch size is B and each section is hour, , For time samples, For the Tao, This is for rounding down.

[0113] Next, the reliable common frequency band of Equation (12) is determined based on the frequency band parameters of each sample; the complex spectrum of the prediction and the band-limited reference before noise addition is calculated using Equation (13); and the difference power is calculated according to Equation (14) on the frequency index selected by Equation (12). and reference power Finally, substitute both into equation (15). The four formulas sequentially complete the frequency selection, transformation, power statistics, and loss normalization:

[0114] (12);

[0115] In the formula, This represents the reliable common frequency band for the b-th sample. and These represent the lower and upper limits of the frequency band, respectively, in Hz; the superscript (b) is the sample number, not a power; the subscript c indicates a common frequency band. This frequency band is determined based on the frequency band parameters of the sample input and the target, using the overlapping portion of their reliable frequency bands, and avoiding the low-energy transition region.

[0116] (13);

[0117] In the formula, The vertical complex spectral coefficients of the predicted profile are defined in the same way as in equation (10); These are the corresponding complex spectral coefficients of the band-limited reference profile before adding noise; The noisy pre-band-limited reference for the b-th sample already includes the amplitude and absorption degradation used by the input, and uses the same geometric transformation and normalization scale as the input and the target; , These are the mask and prediction profile for the b-th sample, respectively; the subscripts k and j outside the parentheses indicate the coefficients extracted for the k-th frequency and j-th channel; the superscript c is only for reference and does not indicate exponentiation.

[0118] (14);

[0119] In the formula, Let b be the average complex spectral difference power of the b-th sample within the reliable common frequency band. The reference average power at the same location; For the selected set of frequency indices, The k-th discrete frequency; The number of spectral coefficients used in the calculation for this sample, | | represents the index number in the set. Let j be the number of channels; the summation over j covers all channels. It must be guaranteed that... Non-empty, avoid ; It represents the square of the modulus of a complex number.

[0120] (15);

[0121] In the formula, B represents the common frequency band complex spectrum loss; B is the number of samples in the current calculation batch. and Calculated by equation (14); Indicates taking Compared with the lower limit of reference power The larger value in; Used to stabilize the denominator, distinct from the integration time τ in equation (6); Used to smooth the square root, and by subtracting... The loss is zero when the spectrum is made completely identical. This embodiment takes... .

[0122] For spatiotemporally varying wavelets, the common frequency band uses a shared and reliable frequency range from the sample metadata, avoiding directly treating the passband at any location as the common frequency band of the entire sample. When using the four-inflection-point parameter, the common intersection of the input and target main passbands within the effective region can be obtained, and low-energy frequencies can be eliminated; the selected frequency band must be non-empty. In this embodiment, the complete profile M=1; when there are filled or invalid regions, the same mask and the same boundary processing are used on the input, target, and reference.

[0123] Finally, the total loss is calculated according to equation (16), and the frequency domain constraint weights are added according to equation (17):

[0124] (16);

[0125] In the formula, This represents the total training loss for the adaptive bandwidth expansion network. , , The losses are represented by equations (9), (11), and (15), respectively. , These are the amplitude spectrum loss weight and the common band complex spectrum loss weight corresponding to the e-th training round, respectively; e is the training round number starting from 0.

[0126] (17);

[0127] In the formula, This represents the smaller of e / 10 and 1; 0.05 and 0.03 are the final weights of the two frequency domain losses, respectively, and 10 is the warm-up round parameter. When e=0, the two weights are 0; when e=10, they reach 0.05 and 0.03 respectively, and remain unchanged thereafter. If the training rounds in the figure start from 1, the position where the final weights are reached corresponds to the 11th round in the figure.

[0128] This embodiment uses the AdamW optimizer, trained for 60 epochs, with an initial learning rate of 2×10⁻⁻⁻⁻⁴. 4 Weight decay 1×10⁻ 4 The gradient norm cutoff threshold is 1.0; the physical batch size is 1, and the effective batch size is 4 obtained through 4 gradient accumulations. B in the formula is still the actual batch size corresponding to each loss calculation. Learning first performs 10 rounds of linear warm-up, and then decays according to the cosine plan; during mixed precision training, the common frequency band complex spectrum loss is calculated using 32-bit floating-point numbers.

[0129] In each round, three losses are examined on the validation set and weighted according to a fixed weight. The model checkpoints are selected in accordance with the selection rules of the original comparative experiment; the test set is used only for independent evaluation. After network training is complete, the parameters are saved for reuse in S3.

[0130] An adaptive bandwidth extension method for non-stationary post-stack seismic data, comprising the following steps:

[0131] S1, Obtain and normalize the profile to be processed;

[0132] Obtain band-limited post-stack seismic profiles, retaining the effective area mask, original dimensions, sampling interval, and trace head information. Normalize the profiles according to the amplitude rules consistent with those used during training and record the scale. Normalization only changes the amplitude scale and does not add time sampling points or seismic traces.

[0133] When the sampling interval of the profile to be processed differs from the training data, it can be resampled to the training interval before S1, and the original sampling position can be restored in S6. When the profile is larger than the training window, an overlapping sliding window can be used, with normalized weights used for fusion in the overlapping area, and the original range can be cropped back in S6. The training window size and temporary padding do not change the number of time samples and channels in the final output.

[0134] S2, construct and project six-channel earthquake priors;

[0135] like Figure 3 As shown in the network entry, the gradients, Laplace responses, and local details in both directions are explicitly calculated from the normalized input X, and then concatenated with the original amplitude according to equation (18). Appropriate padding is applied to the local averaging and difference operations to maintain the spatial dimensions of each channel:

[0136] (18);

[0137] In the formula, X is the normalized single-channel band-limited input profile; and These are the vertical gradient and the lateral gradient, respectively. , These represent the difference operators in the corresponding directions, with the subscript h indicating the transverse direction; For a two-dimensional discrete Laplace response, here It is the Laplace operator; , These are local details in the vertical and horizontal directions, respectively. , For the local averaging operator in the corresponding direction, and the magnitude factor of equation (7) Different; Concat indicates splicing along the channel dimension; These are six-channel prior features.

[0138] Boundary padding is used to maintain consistent dimensions across all prior channels, followed by shallow convolutional projection. The initial features are obtained, where This is a shallow convolution projection. Its output. The six-channel prior... Initial features are obtained through shallow convolution projection. This allows vertical waveforms, lateral continuity, and local boundary information to be expressed at the network entry point.

[0139] S3, multi-scale coding, residual block processing and skip fusion decoding;

[0140] Using the trained parameters, enter Figure 3The diagram illustrates a multi-scale residual encoder-decoder. Taking a 512×512 input as an example, the number of channels for the four encoding scales can be 64, 128, 256, and 512, respectively, with spatial dimensions of 512×512, 256×256, 128×128, and 64×64, respectively. The bottleneck occurs at the coarsest scale where information is aggregated, and the original dimensions are recovered from the three decoding scales through upsampling and same-scale skip connections. This scale configuration is used as an example.

[0141] like Figure 4 As shown, each residual block consists of five sub-layers executed sequentially. Each sub-layer employs pre-normalization, feature transformation, learnable scaling, and residual summation: anisotropic convolution sub-layer, axial and windowed attention sub-layer, vertical Fourier mixing and sample-level modulation sub-layer, local modulation sub-layer, and gated feedforward sub-layer. The specific calculation of sample-level and local-level modulation is S4.

[0142] Anisotropic convolution uses five parallel branches of equation (19) to model local two-dimensional variations, vertical waveforms, lateral in-phase axes, and larger-scale structures, respectively:

[0143] (19);

[0144] In the formula, F is the input feature tensor of the current processing stage; This is an anisotropic convolutional fusion output; , , These are two-dimensional convolution operators with kernel sizes of 3×3, 7×1, and 1×7, respectively, where the first dimension of the kernel size is vertical and the second dimension is horizontal. For dilated convolution operators; It is a depthwise separable convolution operator; This is a fusion projection performed after stitching the branches along the channel dimension. Each branch maintains the same spatial size through padding.

[0145] Axial attention aggregates long-range information along two directions, while window attention supplements the local two-dimensional structure. Then, vertical Fourier mixing is performed according to equations (20) and (21), and S4 is called within the corresponding sublayer:

[0146] (20);

[0147] In the formula, This is a one-dimensional Fourier transform along the vertical sampling axis of the feature tensor, while other dimensions remain unchanged; Let be the input feature tensor of the current processing stage; S is the vertical complex spectrum obtained by transformation. Here, t represents the vertical sampling direction, and the frequency index of the multi-scale latent features is used for feature mixing; the common band loss defined in Hz is calculated on the output profile of the original sampling interval.

[0148] (twenty one);

[0149] In the formula, The frequency domain mixing characteristics after inverse transformation; , These are the real and imaginary parts of the complex spectrum S in equation (20), respectively; To perform a learnable mixing of the real and imaginary parts and form the output complex spectrum, the local wavelet spectrum of equation (1) is used. different; This is the inverse Fourier transform along the vertical sampling axis. It employs a real Fourier transform combined with the corresponding inverse transform, or maintains the conjugate symmetry of the complete spectrum, to achieve... It is a real-valued feature.

[0150] After sample-level and local-level modulation are completed, the gated feedforward sublayer filters and fuses features according to equation (22), and then outputs the learned residuals after scaling and summing:

[0151] (twenty two);

[0152] In the formula, This is the input feature tensor for the current processing stage; , , These are gated branch mapping, value branch mapping, and output projection, respectively. It represents element-wise multiplication between homogeneous tensors, broadcasting by channel or spatial location when necessary; This is the output of the gated feedforward branch. The subscripts g, v, and o of each mapping name identify the gate, value, and output, respectively.

[0153] The residual blocks mentioned above are reused at their respective scales in the encoding path, bottleneck, and decoding path. Skip fusion supplements the decoding features with shallow structural information at the same scale. After all scale processing is completed, the final decoded features enter S5.

[0154] S4 performs sample-level and local-level spectral modulation within the residual block of S3;

[0155] S4 corresponds to Figure 4 The sample-level and local-level modulation stages calculate and return the result on the current feature of the corresponding residual block each time. The sample-level branch first extracts the channel frequency description of the entire sample according to equation (23), and then modulates the channel response according to equation (24):

[0156] (twenty three);

[0157] In the formula, F is the input feature tensor of the current processing stage; | Its vertical amplitude spectrum; For spectral pooling operators, statistical information in the frequency dimension and lateral position is pooled to describe the channel of each sample; MLP stands for channel multilayer perceptron. Let u be the sigmoid activation function, and u be the independent variable of the function. The modulation weights for each sample and each channel are shared at the spatial location of that channel.

[0158] (twenty four);

[0159] In the formula, Features after sample-level modulation; The modulation weights for each sample and each channel are shared at the spatial location of that channel. Indicates will After being broadcast to the corresponding spatial location, they are multiplied element by element. Because Single-level modulation factor Located at (0.5, 1.5), centered at 1.

[0160] The local branch is based on the sample-level modulation features, and generates position-related weights according to equation (25) and performs local modulation:

[0161] (25);

[0162] In the formula, From The extracted local energy description can be achieved by using the local average of the characteristic squares; This is a convolutional mapping that maps the energy description to the position-dependent response; σ is the sigmoid activation function of equation (23); For local modulation weights; These are the features after local modulation; the subscript 1 indicates the local level. and Maintain the same spatial dimensions. Single-level factor. Located at (0.5, 1.5), the overall scaling range of the two-level multiplication is (0.25, 2.25).

[0163] Both levels of weights are obtained through training, and the response is adjusted based on the overall spectral state and local energy. Local energy only provides conditional information and is not directly equivalent to geological confidence. The modulated features are returned to their respective sublayers for further residual summation and subsequent processing.

[0164] S5, predicts full-band correction residuals and constrains amplitude;

[0165] The residual prediction head maps the final decoded features to single-channel full-band corrected residuals. Add the normalized input X according to equation (26):

[0166] (26);

[0167] In the formula, X is the normalized band-limited input profile; The residual prediction head outputs a single-channel full-band correction residual based on the final decoded features; This is the broadband prediction profile after summing the residuals; all three have the same number of time samples and channels. Amplitude constraints apply to... The constrained profile is then normalized and scaled to obtain the final output.

[0168] Determine the adaptive cutoff lower limit based on the amplitude quantiles or dynamic range of the input valid region. and upper limit point by point Constraint to The statistical rules for the upper and lower limits are consistent during the training and inference phases. This processing is used to suppress abnormal amplitudes and overcompensation; the single-channel residual simultaneously handles low- and high-frequency compensation without changing the profile size.

[0169] S6, restores the amplitude scale and outputs the enhanced profile;

[0170] The enhanced profile obtained in S5 is restored to the amplitude scale recorded in S1. When using an overlapping sliding window, each block is first restored to the same amplitude scale, then weighted and fused in the overlapping area, and temporary fill is removed; if preprocessing resampling exists, the result is restored to the original sampling position, and trace head information is also restored. The final output is a broadband post-stack seismic enhanced profile consistent with the original input time sample number and trace number.

[0171] Experimental verification process:

[0172] Select an Ormsby wavelet sample and generate it using equations (4) and (5). Figure 5 of , The parameter field changes smoothly with time and channel number. The spatial trends of the input and the target correspond to each other, and the frequency band parameters are different, indicating that the sample can represent a spatiotemporally non-stationary wavelet.

[0173] The training process for f is carried out for 60 rounds according to the multi-domain objective described above. Figure 6 In the figure, the mask Charbonnier reconstruction loss decreased from 0.032932 to 0.013408, the vertical logarithmic amplitude spectrum loss decreased from 0.035557 to 0.009102, and the common band complex spectrum loss decreased from 0.176998 to 0.028885; the three types of losses generally decreased and then tended to stabilize in the later stages. The three abbreviations in the figure correspond one-to-one with the three losses listed in the invention description.

[0174] Based on the processing results of S1 to S6, 300 samples are derived from 1000 independent synthetic test samples to form a plotting subset. According to the frequency band endpoints of each sample, the low-frequency extended band, common frequency band and high-frequency extended band are mapped to [−1,0], [0,1] and [1,2] respectively. Figure 7 The vertical axis represents the amplitude spectrum ratio relative to a broadband target, in dB. 0 dB indicates that the amplitude spectrum ratio is 1.

[0175] Figure 7 In the corresponding plotting subset, the average gaps at the low-frequency and high-frequency ends of the narrow-band input are approximately -12.90 dB and -18.92 dB, respectively. The output of this method is approximately -0.75 dB and -0.55 dB in the corresponding frequency bands, and approximately -0.04 dB in the common frequency band, indicating that the amplitude spectra at both ends are compensated and the average amplitude spectrum in the common frequency band is close to the target. The amplitude spectrum ratio itself does not directly prove phase consistency; phase preservation is jointly constrained by the complex spectrum loss in the common frequency band and the corresponding evaluation.

[0176] Independent tests also compared waveform correlation, low-frequency to high-frequency correlation, and logarithmic spectral error. Ablation comparisons showed that removing structural priors weakened structural modeling, removing dual-scale modulation affected non-stationary adaptation, and removing common band loss increased common band complex spectral error. Figures 5 to 7 The values ​​are used to illustrate the implementation results, without limiting the number of training rounds, sample size, or parameter values.

[0177] The on-site data lacks strict broadband true values, and it cannot be concluded that the missing frequency has been accurately recovered based solely on the increase in external frequency band energy; the processing results are combined with common frequency band back-filter correlation, phase deviation, residual profile and amplitude calibration for quality inspection.

[0178] The feasibility and effectiveness of the present invention will be further illustrated below with an example based on the Marmousi2 model.

[0179] Example 1

[0180] This embodiment uses the Marmousi2 elastic model to test the processing capability of the method of the present invention under complex structural conditions. The model includes multi-scale faults, steeply dipping interfaces, and deep arc-shaped reflections, and its structural complexity is significantly greater than that of the layered background used in the training samples. The P-wave reflectivity of the model is convolved with the Ricker wavelet with a dominant frequency of 100Hz to obtain a synthetic record. Then, the spectrum is reconstructed into a band-limited profile with a dominant frequency of 18Hz using the Gaussian spectral ratio method, which serves as the band-limited post-stack seismic profile to be processed in this embodiment. The time sampling interval of this profile is 4ms. After removing the top 248ms to avoid Fourier transform loop effects, the time range is 248-2240ms, with a total of 499 time sampling points and 2671 traces.

[0181] The above profile is processed according to steps S1 to S6: after normalization, a 512×512 overlapping sliding window is used for block-by-block processing, with an overlap width of 128 sampling points. The overlapping areas are weighted and fused, and each block reuses the same set of pre-trained network parameters. The final output is a broadband enhanced profile with the same size as the input. Since the model does not have a strict broadband ground truth, a no-reference index is used to evaluate the processing effect: the input profile energy is concentrated in 5-30Hz, which is defined as the common frequency band. 0-5Hz and 30-50Hz are used as the low-frequency end and high-frequency end observation frequency bands, respectively. The average amplitude spectrum gain within the frequency band is used to measure the frequency band expansion intensity, the correlation coefficient between the output after back-filtering in 5-30Hz and the input is used to measure the degree of common frequency band preservation, and the correlation coefficient between adjacent channels within the band is used to measure lateral continuity.

[0182] Figure 8 In the diagram, (a) is the band-limited input profile, (b) is the broadband enhanced profile processed by the method of this invention, and (c) and (d) are local magnified views of the two profiles in the fault-dense region (traces 900-2000 and 700-1500ms), respectively. It can be seen that the method of this invention maintains the main structural morphology, reflection position, and polarity of the Marmousi2 model, without significant structural misalignment; the processed profile exhibits significant wavelet compression, enhanced interlayer contrast, and significantly improved vertical resolution in weak reflection zones and fault-related areas; the location of steeply dipping interfaces and the lateral traceability of weak phase axes are also clearer.

[0183] Quantitative results show that the method of this invention achieves average amplitude spectrum gains of +11.86dB and +12.91dB at the low-frequency end (0-5Hz) and high-frequency end (30-50Hz), respectively, realizing balanced bi-lateral expansion at both low and high frequencies. In contrast, the high-frequency gains of the three comparative networks, U-Net, ResUNet, and SRResUNet, are +4.78dB, +4.98dB, and +0.58dB, respectively. The common band gain of the method of this invention is only +0.58dB, and the back-filter correlation coefficient for 5-30Hz is 0.9835, indicating that the amplitude-phase relationship of the common band is effectively maintained. The correlation coefficient between adjacent channels at the high-frequency end is 0.9511, higher than the 0.9431 of the input profile itself, indicating that the newly added high-frequency energy is organized laterally along the real structure without generating disordered texture. These results verify that the method of this invention still has good effectiveness and adaptability even on models with construction complexity exceeding the training background.

[0184] Of course, the above description is not intended to limit the present invention, and the present invention is not limited to the examples given above. Any changes, modifications, additions or substitutions made by those skilled in the art within the scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. An adaptive bandwidth spreading method for non-stationary post-stack seismic data, characterized in that, Includes the following steps: S1, obtain the band-limited post-stack seismic profile to be processed, normalize the amplitude of the band-limited post-stack seismic profile while keeping its time sampling number and trace number unchanged; S2, construct a priori seismic features based on the zone-limited post-stack seismic profile, and project the a priori seismic features as initial features; S3 uses a pre-trained multi-scale residual encoder-decoder to process initial features, and uses adaptive frequency band expansion residual blocks of different scales to extract vertical waveform, lateral in-phase axis continuity and frequency domain interaction information. It also fuses encoded and decoded features through skip connections. S4. Inside the residual block in step S3, sample-level channel modulation weights are generated based on the vertical amplitude spectrum of the current feature, and position-related modulation weights are generated based on the local energy of the sample-level modulation feature. The frequency band compensation intensity is adjusted in two stages, and the modulated feature is returned to the residual block for further processing. S5. Based on the final decoding features, predict the single-channel full-band correction residual, add it to the normalized band-limited post-stack seismic profile, and apply adaptive amplitude constraints based on the effective amplitude statistics or dynamic range of the input profile. S6 restores the normalized scale of the output profile, outputting a broadband post-stack seismic enhancement profile with the same size as the band-limited post-stack seismic profile.

2. The adaptive bandwidth extension method for non-stationary post-stack seismic data as described in claim 1, characterized in that, In step S2, the a priori seismic features include six categories: original amplitude, vertical gradient, lateral gradient, Laplace response, vertical local detail, and lateral local detail of the band-limited stacked seismic profile. These features are then stitched together and projected according to the following formula: ； In the formula, X is the normalized single-channel band-limited input profile; and These are the vertical gradient and the lateral gradient, respectively. , These represent the difference operators in the corresponding directions, with the subscript h indicating the transverse direction; It is a two-dimensional discrete Laplace response. It is the Laplace operator; ; These are local details in the vertical and horizontal directions, respectively. , These are the local averaging operators for the vertical and horizontal directions, respectively; Concat indicates concatenation along the channel dimension; It is a six-channel prior feature; Boundary padding is used to maintain consistent dimensions across all prior channels, followed by shallow convolutional projection. We obtain the initial features, where, This is a shallow convolution projection. Output it.

3. The adaptive bandwidth extension method for non-stationary post-stack seismic data as described in claim 2, characterized in that, In step S3, the multi-scale residual encoder-decoder includes an encoding path that progressively reduces the spatial size and increases the number of channels, a bottleneck layer located at the coarsest scale, and a decoding path that progressively restores the spatial size; the decoding path fuses the upsampled features with the same-scale features of the encoding path at each corresponding scale. Its pre-training process is completed before performing step S1, and includes: The same reflectivity model is convolved with wide-band and narrow-band sub-wavefields respectively to form paired training samples; The predicted profile is obtained by using network forward computation; The network parameters are updated by weighting the Charbonnier reconstruction loss, vertical logarithmic amplitude spectrum loss, and common band complex spectrum loss. For a network that has been trained, the trained parameters are reused when processing different profiles.

4. The adaptive bandwidth extension method for non-stationary post-stack seismic data as described in claim 3, characterized in that, In step S3, the adaptive band-expanded residual block undergoes anisotropic seismic convolution, axial attention and window attention, vertical full-band Fourier mixing, local spectral modulation, and gated convolution feedforward processing sequentially, following the methods of normalization, feature transformation, learnable residual scaling, and residual summation; wherein: Anisotropic convolutions are fused in parallel according to the following formula: local 2D, vertical long kernel, horizontal long kernel, dilated, and depth-separable convolution branches. ； In the formula, F is the input feature tensor of the current processing stage; This is an anisotropic convolutional fusion output; , , These are two-dimensional convolution operators with kernel sizes of 3×3, 7×1, and 1×7, respectively, where the first dimension of the kernel size is vertical and the second dimension is horizontal. For dilated convolution operators; It is a depthwise separable convolution operator; This is a fusion projection performed after stitching together the branches along the channel dimension; Axial attention establishes long-range dependencies along both the vertical and horizontal axes, while window attention aggregates structural information within a local two-dimensional window. Vertical Fourier mixing is implemented according to the following formula: ； In the formula, Indicates the vertical sampling direction. This is a one-dimensional Fourier transform along the vertical sampling axis of the feature tensor, while other dimensions remain unchanged; The vertical complex spectrum obtained by transformation; ； In the formula, The frequency domain mixing characteristics after inverse transformation; , These are the real and imaginary parts of the complex spectrum S, respectively; A mapping that allows for the learnable mixing of the real and imaginary parts and forms the output complex spectrum; This is the inverse Fourier transform along the vertical sampling axis; Gated convolution feedforward processing is performed according to the following formula: ； In the formula, , , These are gated branch mapping, value branch mapping, and output projection, respectively. It represents element-wise multiplication between homogeneous tensors, broadcasting by channel or spatial location when necessary; For the output of the gated feedforward branch, the subscripts of each mapping name , , Identify the gating, value, and output respectively.

5. The adaptive bandwidth extension method for non-stationary post-stack seismic data as described in claim 4, characterized in that, Step S4 includes generating channel weights based on the characteristic vertical amplitude spectrum and performing sample-level modulation, followed by local-level modulation based on local energy; the sample-level and local-level modulation execution processes are as follows: ； In the formula, Its vertical amplitude spectrum; For spectral pooling operators, statistical information in the frequency dimension and lateral position is pooled to describe the channel of each sample; MLP stands for channel multilayer perceptron. for Shape activation function, This is the independent variable of the function; The modulation weights for each sample and each channel are shared at the spatial location of that channel; ； In the formula, Features after sample-level modulation; because Single-level modulation factor Located at (0.5, 1.5), centered at 1; ； ； In the formula, From Extracted local energy description; A convolutional mapping that maps energy descriptions to location-dependent responses; For local modulation weights; The features are those after local modulation, with the subscript 1 indicating the local level; Single-level factor Located at (0.5, 1.5), the overall scaling range of the two-level multiplication is (0.25, 2.25).

6. The adaptive bandwidth extension method for non-stationary post-stack seismic data as described in claim 5, characterized in that, In step S5, the residual prediction head outputs a full-band correction residual used to compensate for missing or attenuated components at the low and high frequencies, which is added to the input according to the following formula: ； In the formula, The residual prediction head outputs a single-channel full-band correction residual based on the final decoded features; This is the broadband prediction profile after summing the residuals; all three have the same number of time samples and channels; amplitude constraints apply to... The constrained profile is then normalized and scaled to obtain the final output. right When applying amplitude constraints, the lower limit is determined based on the input effective amplitude quantiles or dynamic range. and upper limit ,and Each predicted amplitude is truncated to this interval; , These represent the adaptive lower and upper limits, respectively. The amplitude statistics only use the effective region, and the constraint rules remain consistent during training and inference.

Images