A curvelet domain de-azimuth processing method based on reflection coefficient attribute volume

Abstract

This invention provides a curvewave domain-based method for removing strong axes based on reflection coefficient attribute volumes, relating to the field of seismic data processing technology. The method includes: performing spectral inversion processing on raw seismic data to generate a reflection coefficient volume; based on the reflection coefficient volume, performing multi-scale, multi-directional layer-by-layer analysis in the curvewave domain, i.e., extracting reflection coefficient curves for each analysis window, determining inflection point locations by calculating the zero point of the second derivative of the curves to identify abrupt changes in the reflection coefficient; simultaneously, constructing local discrete surfaces based on the reflection coefficient volume, identifying local geometric anomaly regions of the reflection coefficient volume by fitting surface patches and calculating Gaussian curvature; and extracting strong reflection coefficient components by combining the determined inflection point locations and the identified local geometric anomaly regions. This invention can achieve precise stripping of strong reflection energy and faithful protection of effective weak signals, improving the reliability of seismic data in identifying and analyzing subsurface lithology and geological structures.

Description

Technical Field

[0001] This invention relates to the field of seismic data processing technology, and in particular to a curve wave domain strong axis removal method based on reflection coefficient attribute volume. Background Technology

[0002] The core of seismic data interpretation is to obtain the lithological characteristics of subsurface strata and clarify the geological structure. The integrity and authenticity of seismic reflection signals are key to the reliability of interpretation results. In actual exploration, if there are thick layers of extremely high-velocity or extremely low-velocity lithology on the top plate of the target layer, the difference in wave impedance between them and the underlying strata will generate extremely strong seismic reflection waves; while the wavelength and sidelobe characteristics of the seismic wavelet will cause strong axial energy diffusion, strongly shielding the weak reflection signals from the underlying strata.

[0003] Existing methods for solving this problem have the following drawbacks: They struggle to accurately define the spatial distribution and intensity boundaries of strong reflection energy, and cannot selectively strip away this energy. Due to incomplete stripping, residual strong reflection components continue to mask the weak seismic responses of underlying strata, resulting in insufficient accuracy for critical tasks such as stratigraphic lithology identification and geological structure assessment. Furthermore, the boundary between strong reflections and weak signals is blurred; traditional methods are prone to over-stripping during strong reflection suppression, leading to the accidental deletion of effective weak signals from underlying strata, compromising the fidelity of seismic data, and consequently affecting the accuracy of subsequent geological structural interpretation and the delineation of favorable geological targets. For example, in seismic exploration projects containing thick layers of special lithology, the strong axes formed by these lithologies severely suppress weak reflection signals from underlying strata. Because the spatial range and intensity gradient of strong reflections cannot be accurately determined, either residual strong reflection energy continues to obscure the target strata response, leading to significant deviations in the delineation of geological target distribution; or effective weak signals are mistakenly deleted, resulting in misjudgments of stratum thickness or the omission of crucial geological information. Summary of the Invention

[0004] The technical problem to be solved by the present invention is to provide a curve domain strong axis removal processing method based on the reflection coefficient attribute volume, so as to achieve accurate stripping of strong reflection energy and effective protection of weak signals, thereby improving the reliability of seismic data in identifying and analyzing underground strata lithology and geological structure.

[0005] To solve the above-mentioned technical problems, the technical solution of the present invention is as follows: A method for de-strong-axis processing in the curve wave domain based on reflection coefficient property volumes, the method comprising: Step 1: Perform spectral inversion processing on the original seismic data to generate a volume of reflection coefficients; Step 2: Based on the reflection coefficient volume, perform multi-scale and multi-directional layer-by-layer analysis in the curve domain. Specifically, for each analysis window, extract the reflection coefficient curve and determine the inflection point by calculating the zero point of the second derivative of the curve to identify the abrupt change characteristics of the reflection coefficient. At the same time, construct a local discrete surface based on the reflection coefficient volume, and identify local geometric anomaly regions of the reflection coefficient volume by fitting surface patches and calculating Gaussian curvature. Step 3: Based on the determined inflection point locations and the identified local geometric anomaly regions, extract the strong reflection coefficient components; Step 4: Based on the strong reflection coefficient components, establish a discrete data lattice of the spatial distribution of the strong reflection coefficient components; based on the discrete data lattice, construct a wave field energy topology skeleton that reflects the spatial structure of the strong reflection energy; define several characteristic control points located inside, at the edge and outside of the strong reflection energy cluster in the wave field energy topology skeleton; based on the characteristic control points, fit to form an implicit energy isomorphic manifold that characterizes the spatial distribution of the strong reflection energy. Step 5: Extract the morphological properties and energy gradient distribution properties of the implicit energy isomorphic manifold, and calculate the strong reflection stripping optimization coefficient accordingly. Step 6: Adaptively adjust the strong reflection coefficient component using the strong reflection stripping optimization coefficient. Subtract the adjusted strong reflection coefficient component from the reflection coefficient volume to obtain the reflection coefficient volume after removing the strong axis. Step 7: Reconstruct the seismic data based on the reflection coefficients after removing the strong axis to complete the strong axis removal process.

[0006] The above-described solution of the present invention has at least the following beneficial effects: By employing techniques such as generating accurate reflection coefficient volumes through spectral inversion, combining curve domain multi-scale and multi-directional layer-by-layer analysis with second-derivative inflection point identification and Gaussian curvature geometric anomaly judgment, extracting strong reflection components from comprehensive features, constructing wavefield energy topological skeletons and implicit energy equivalent manifolds, calculating adaptive strong reflection stripping optimization coefficients, and then stripping strong reflection components and reconstructing seismic data after weighted adjustment, this approach overcomes the technical problem of existing methods failing to accurately define the spatial distribution range and intensity boundary of strong reflection energy, leading to incomplete or excessive stripping of strong reflection energy and damage to effective weak signals. This approach achieves accurate stripping of strong reflection energy and faithful protection of effective weak signals, effectively improving the reliability of seismic data in identifying and analyzing underground strata lithology and geological structures. Attached Figure Description

[0007] Figure 1 This is a schematic flowchart of a curve wave domain de-strong axis processing method based on the reflection coefficient attribute volume provided by an embodiment of the present invention. Detailed Implementation

[0008] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

[0009] like Figure 1 As shown, an embodiment of the present invention proposes a method for curvature domain de-strong axis processing based on the reflection coefficient property volume. The method includes the following steps: Step 1: Perform spectral inversion processing on the original seismic data to generate a volume of reflection coefficients; Step 2: Based on the reflection coefficient volume, perform multi-scale and multi-directional layer-by-layer analysis in the curve domain. Specifically, for each analysis window, extract the reflection coefficient curve and determine the inflection point by calculating the zero point of the second derivative of the curve to identify the abrupt change characteristics of the reflection coefficient. At the same time, construct a local discrete surface based on the reflection coefficient volume, and identify local geometric anomaly regions of the reflection coefficient volume by fitting surface patches and calculating Gaussian curvature. Step 3: Based on the determined inflection point locations and the identified local geometric anomaly regions, extract the strong reflection coefficient components; Step 4: Based on the strong reflection coefficient components, establish a discrete data lattice of the spatial distribution of the strong reflection coefficient components; based on the discrete data lattice, construct a wave field energy topology skeleton that reflects the spatial structure of the strong reflection energy; define several characteristic control points located inside, at the edge and outside of the strong reflection energy cluster in the wave field energy topology skeleton; based on the characteristic control points, fit to form an implicit energy isomorphic manifold that characterizes the spatial distribution of the strong reflection energy. Step 5: Extract the morphological properties and energy gradient distribution properties of the implicit energy isomorphic manifold, and calculate the strong reflection stripping optimization coefficient accordingly. Step 6: Adaptively adjust the strong reflection coefficient component using the strong reflection stripping optimization coefficient. Subtract the adjusted strong reflection coefficient component from the reflection coefficient volume to obtain the reflection coefficient volume after removing the strong axis. Step 7: Reconstruct the seismic data based on the reflection coefficients after removing the strong axis to complete the strong axis removal process.

[0010] In this embodiment of the invention, it should be noted that the aforementioned layer-by-layer analysis specifically refers to: after obtaining a curved wave domain coefficient volume by performing a curved wave transformation on the reflection coefficient volume, a three-dimensional analysis time window is set within the curved wave domain coefficient volume based on the spatial range of the target geological stratum, and feature extraction is performed; this method, through multi-scale and multi-directional time window division, can simultaneously capture the abrupt changes in the reflection coefficient along the strike and dip of the strata and local geometric anomalies; by generating an accurate reflection coefficient volume through spectral inversion, a solid foundation is laid for strong reflection identification; the multi-scale and multi-directional layer-by-layer analysis in the curved wave domain, combined with the identification of second-derivative inflection points and Gaussian curvature geometric anomalies... The method significantly improves the positioning accuracy of strong reflection areas. By constructing a wavefield energy topology framework based on discrete data points and fitting an implicit energy isomorphic manifold using feature control points, the spatial distribution of strong reflection energy is clearly presented. By extracting manifold attributes and calculating the strong reflection stripping optimization coefficient, adaptive adjustment of the strong reflection coefficient components is achieved. This not only accurately strips strong reflection energy but also effectively protects the underlying weak signals. Finally, the reconstructed seismic data truly reflects the seismic response characteristics of underground strata, effectively improving the reliability of seismic data in identifying and analyzing stratigraphic lithology and geological structures, and providing high-quality data support for subsequent related work.

[0011] In a preferred embodiment of the present invention, step 1 above may include: Step 1.1 involves preprocessing the raw seismic data to obtain a preprocessed seismic data volume. This includes: collecting raw seismic data within the target work area's survey line range; converting all data into a standardized format, ensuring the integrity of the trace header information for each trace; maintaining a sampling interval within a reasonable range of 1 to 4 ms to guarantee data format uniformity and compatibility with subsequent processing; conducting a comprehensive quality screening of the converted raw data, identifying and removing invalid data such as empty or bad traces by verifying the trace header information; marking abnormal traces with serious acquisition errors or missing local data, and supplementing them with reasonable data from adjacent valid traces through interpolation to ensure data integrity; selectively suppressing random and useless signals mixed in the data, strictly preserving the dominant frequency and phase characteristics of valid seismic reflection signals during processing, without destroying the original form of useful signals; and separating coherent interference signals that affect data quality, such as interference from industrial activities and surface waves, from valid signals through separation processing to reduce the adverse effects of such interference on subsequent analysis.

[0012] Subsequently, amplitude consistency adjustment was performed. Based on the differences in excitation conditions at different excitation points and the differences in reception at different receiving points, the data amplitude was calibrated channel by channel to eliminate amplitude imbalance caused by factors such as uneven excitation charge and differences in receiving equipment status. Then, time difference correction processing was carried out. First, based on relevant data collected in the field, the differences in surface elevation and the thickness of the near-surface low-velocity zone were calculated to correct the reference surface time difference of the data. Then, through multiple iterations, residual time difference deviations were corrected to ensure the continuity and integrity of the same seismic reflection phase axis in the lateral distribution. Finally, a preprocessed seismic data volume with significantly improved signal-to-noise ratio, balanced and stable amplitude distribution, and effective time difference correction was obtained.

[0013] Step 1.2: Based on the preprocessed seismic data volume, extract seismic wavelets. Specifically, this includes: accurately delineating a range of 50 to 100 ms above and below the target layer from the preprocessed seismic data volume. This range must completely cover the reflection signals of the target layer and adjacent strata, while avoiding areas with chaotic signals such as fault fracture zones and areas with strong interference. Select continuous data gathers with continuous reflection phase axes, stable amplitudes, and low proportion of interference signals as the basic data for wavelet extraction. Perform amplitude normalization processing on the selected data gathers to adjust the peak amplitude of each signal to the same standard, eliminating the influence of single-channel energy differences on the extraction results. Then, perform time-domain smoothing processing to weaken residual random noise and retain the core waveform characteristics of the reflected wave.

[0014] Effective waveform segments of the reflecting phase axis are extracted from each channel of the processed channel set. The time-domain distribution characteristics of all effective waveform segments are statistically analyzed. Then, these waveform segments are converted to the frequency domain, and the amplitude values ​​of each channel are statistically analyzed at each frequency point and the arithmetic mean is taken to obtain the average amplitude spectrum reflecting the overall reflection characteristics of the work area. At the same time, the average phase spectrum is calculated using a weighted average method. Specifically, for each frequency point, n data channels participating in the calculation are first determined, and weights w1, w2, ..., wn are assigned to each channel (where w1 + w2 + ... + wn = 1). The weight allocation is based on the signal strength and signal-to-noise ratio of the channel. The channel with the greater signal strength and the higher the signal-to-noise ratio has a larger weight value. Then, the average phase is calculated, as shown in the following formula: By calculating the average phase at this frequency point, and completing the calculation point by point, a final average phase spectrum is formed, thereby improving the reliability of the phase spectrum. Each of the n data channels involved in the average phase calculation at a specific frequency point has its own phase value at that frequency point.

[0015] The dominant frequency range of seismic reflections in the work area is determined based on the average amplitude spectrum. The center frequency of this range is used as the dominant frequency of the wavelet, and the wavelet length is calculated in combination with the dominant frequency. The wavelet is ensured to cover 2 to 3 complete main cycles. For example, when the dominant frequency is 20Hz, the wavelet length is set to 40ms, which preserves the core reflection information and avoids excessive extension of side lobes. Based on the average amplitude spectrum and phase spectrum, an initial wavelet is generated in the time domain. Then, the initial wavelet is subjected to dual stability verification: first, waveform stability verification, which calculates the cross-correlation coefficient between wavelets corresponding to different data channels, requiring the coefficient value to be no less than 0.9 to ensure the consistency of wavelet waveforms in each channel; second, phase stability verification, which compares the phase spectra of different wavelets, requiring the phase difference within the dominant frequency range to not exceed 10°, and checking whether the phase spectrum is smooth and continuous, eliminating abnormal wavelets with phase abrupt changes.

[0016] For the verified initial wavelet, an iterative optimization process is initiated: the wavelet is convolved with the original data gather to obtain the predicted seismic response, and the fitting error between the predicted response and the actual data is calculated for each trace; if the error exceeds a preset threshold, usually set at 5%, the amplitude and phase near the dominant frequency are corrected: for frequency points with large amplitude deviations, the corresponding amplitude values ​​are adjusted to fit the average amplitude spectrum; for areas with phase shifts, a smoothing iterative method is used to correct the phase curve until the fitting error is reduced below the threshold; through multiple iterative adjustments, a seismic wavelet with regular waveform, concentrated dominant frequency, stable phase, and high fit with the seismic data of the target work area is finally obtained, providing reliable support for subsequent spectral inversion to generate an accurate reflection coefficient volume.

[0017] Step 1.3: Combining the acquired geological data and well logging data of the work area, construct a low-frequency model for spectral inversion. This includes: comprehensively collecting geological and well logging data of the target work area that meet the requirements of strong reflection stripping. The geological data should include a 1:50,000 stratigraphic plan, lithological columnar profile, fault occurrence data, and stratigraphic contact relationship report, focusing on confirming the vertical thickness variation, lateral distribution boundary, and interface characteristics with the underlying strata of extremely high-velocity or extremely low-velocity lithologies; the well logging data should include complete collection of sonic and density logging curves from all wells, supplemented by natural gamma and resistivity logging curves. The sonic and density data provide support for the basic calculation of wave impedance at low frequencies, while the natural gamma curve can help determine the lithology type. Taking coal-bearing strata exploration as an example, for the thick coal and rock layers at the top of the target stratum, logging data from 12 wells were collected to determine the thickness range of 5 to 18 meters, the velocity range of 1800 to 2200 m / s, and the velocity difference threshold between the coal and rock sections and the underlying sandstone (3200 to 3800 m / s). The logging data underwent refined preprocessing. First, the logging data was verified meter by meter using the caliber curve. For sections with enlargement exceeding 10%, the trend extrapolation method of logging values ​​from adjacent effective sections was used for correction. Then, based on multiple representative wells covering different lithologies and structural areas in the work area, a unified calibration standard was established using high-quality logging data with stable signals, minimal interference, and conformity to geological patterns. This standard was used to correct the deviation of the sonic logging data, controlling the numerical deviation within ±50 m / s and eliminating systematic errors. Finally, segmented Z-score standardization was performed according to lithological categories such as coal, mudstone, sandstone, and carbonate rocks, eliminating systematic deviations between wells while preserving the differences in the physical properties of the lithology itself.

[0018] Next, based on the seismic line spacing and depth sampling rate of the work area seismic data, the three-dimensional space of the work area was divided into fine grid units of 20m×20m×0.5m. The horizontal grid matched the seismic trace spacing, and the vertical grid was less than 1 / 2 of the minimum lithological thickness to ensure accurate characterization of thin layers such as coal and rock. The coordinates of each stratigraphic interface were projected onto the grid according to the geological stratification data to clarify the stratigraphic affiliation of each grid unit. Then, the standardized well logging data was accurately assigned to the corresponding grid units according to the well coordinates. For grids without well logging data, initial assignment was performed based on the lithological distribution patterns of adjacent wells to form the initial framework of the model. An improved Kriging interpolation algorithm was used to construct the initial model. On the basis of traditional interpolation, a lithological weight factor and a structural constraint threshold were introduced: the weight factor for the same lithological region was set to 0.8 to 0.9 to enhance spatial continuity, and the weight factor for lithological interface regions was set to 0.3 to 0.5 to preserve attribute abrupt changes. Interpolation barriers were set for grid units on both sides of the fault based on fault attitude data to avoid stratigraphic attribute confusion caused by direct interpolation across faults. For the strong reflection interface area where coal and sandstone meet, an interpolation control point is set every 5m to increase the density of interpolation calculations and accurately characterize the velocity and density gradient at the interface.

[0019] Secondly, model training was conducted to improve accuracy. 70% of the wells in the work area were selected as the training set, covering different lithological distribution areas and structural locations to ensure that the training data fully reflects the geological characteristics of the work area. 30% of the wells were used as the validation set, including two wells near faults and three wells with abrupt changes in coal and rock thickness, to verify the model's performance under complex conditions. The measured acoustic and density data of the training set wells were used as model training labels, and the measured data of the validation set were used as accuracy evaluation standards. The label data corresponded point by point with the grid cells to ensure consistency between training and actual scenarios. A multi-dimensional error evaluation system was constructed, requiring not only an absolute velocity error ≤100m / s and an absolute density error ≤0.1g / cm3, but also a relative error ≤5% and a correlation coefficient ≥0.92. The correlation coefficient is used to evaluate the model's ability to represent abrupt changes in lithological interfaces, avoiding blurring of strong reflective interfaces due to excessive model smoothing.

[0020] Then, multiple rounds of iterative optimization and scenario-based verification were conducted. After the first iteration, adjustments were made to different error-exceeding areas: in lithologically homogeneous areas, the interpolation coverage was expanded, and the interpolation radius was adjusted from 500m to 800m, while the interpolation smoothness parameter was optimized to further improve the spatial continuity of formation attributes; in lithological interface areas, the interpolation coverage was reduced, and the interpolation radius was adjusted to 200m, while the numerical difference of lithological weight factors was increased, making the weights within the same lithology more concentrated and the weights of different lithological interfaces more distinct, thereby strengthening the local abrupt change characteristics at lithological interfaces; in areas near faults, the structural constraint weight was increased from 0.6 to 0.8, and the fault barrier boundary parameters were optimized to avoid cross-fault interpolation interference. After each iteration, the error index was recalculated. If the relative error change rate was ≤0.3% for three consecutive iterations and the correlation coefficient was ≥0.95, the model training was considered converged; if convergence was not achieved after 80 iterations, control points of logging data near faults were added and the iteration was restarted. For coal-bearing strata scenarios, additional strong reflection interface adaptability verification was conducted. Velocity gradient data of the coal-rock and sandstone interface in the model were extracted and matched with the measured strong reflection interface location of the seismic earthquake. The interface depth deviation was required to be ≤2m and the velocity gradient deviation to be ≤300m / s / m. For areas with excessive deviation, the lithology weight factor of the grid near the interface was locally adjusted to ensure that the lithology interface represented by the model corresponds accurately to the strong reflection interface of the seismic earthquake.

[0021] The low-frequency model of this invention has significant core advantages: it is an improvement on the traditional geostatistical interpolation model, incorporating a structural constraint enhancement module and a lithology-sensitive feature adaptation module, perfectly meeting the core requirements of strong reflection stripping and weak signal protection, and solving the shortcomings of traditional models that ignore structural abrupt changes and poor lithology adaptability; through refined grid division, segmented standardized preprocessing, improved interpolation algorithms, and multi-dimensional error evaluation, the characterization error of extremely high-velocity or extremely low-velocity lithologies is controlled within 5%, providing a stable and accurate low-frequency background for spectral inversion and effectively reducing the volumetric ambiguity of reflection coefficients; it is particularly effective for thick-layered extremely high-velocity / extremely low-velocity lithologies. The targeted processing of sexual development zones can accurately characterize the low-frequency attribute differences of strong reflection interfaces, helping spectral inversion to highlight the characteristics of strong reflection coefficients and avoid extraction bias. The regional parameter adjustment and multi-round iteration mechanism can adapt to work areas with different structural complexities and lithological combinations, effectively improving the applicability of the method of this invention. At the same time, the low-frequency background field output by the model is continuous and stable, and when combined with seismic wavelet constraints, it can improve the convergence speed of spectral inversion by more than 30%. The generated reflection coefficient volume can clearly distinguish between strong reflection components and weak signal components, laying a solid foundation of accurate attributes for subsequent multi-scale and multi-directional layer-by-layer analysis in the curve wave domain.

[0022] Step 1.4: Using the extracted seismic wavelet and the low-frequency model used for spectral inversion as constraints, perform constrained sparse pulse inversion iterative calculations on the preprocessed seismic data volume until convergence, obtaining the initial reflection coefficient volume. Specifically, this includes: using the extracted seismic wavelet and the constructed low-frequency model as core constraints, ensuring that their sampling rates are consistent with the preprocessed seismic data (1 to 4 ms), and that the trace head information accurately corresponds to the original seismic data, ensuring the adaptability of the constraints to the processing object. The preprocessed seismic data undergoes format adaptation processing, unifying the data storage structure and identification rules to meet the input requirements of the constrained sparse pulse inversion iterative calculation; initial parameters for the inversion iteration are set based on the stratigraphic reflection characteristics of the work area: the reflection coefficient sparsity threshold is adjusted to 0.05 to 0.1 according to the stratigraphic reflection density, with higher values ​​in dense reflection areas and lower values ​​in sparse reflection areas; the iteration step size is set to 0.01 to 0.03 to balance iteration efficiency and adjustment accuracy; the maximum number of iterations is set to 50 to 100 to avoid computational redundancy caused by excessive iteration.

[0023] After initiating iterative calculations, point-by-point convolution operations are first performed between the current reflection coefficient data for each trace and the seismic wavelet; the temporal correspondence between the reflection coefficient sequence and the seismic wavelet is then determined: let the reflection coefficient sequence of the current trace be... r represents the reflection coefficient sequence of the current trace, and n is the number of reflection coefficient sampling points, arranged in chronological order; the seismic wavelet sequence is... , w represents the seismic wavelet sequence, m is the wavelet length, containing 2 to 3 main periods. The weighted superposition formula for the convolution operation is: the predicted seismic response k-th sample value = r1×wk + r2×wk-1 + ... + rk×w1, where k ranges from 1 to n+m-1. When the index exceeds the sequence range, the corresponding value is 0. During the operation, according to the temporal distribution of the reflection coefficient sequence, each reflection coefficient value is weighted (reflection coefficient value multiplied by wavelet point value) and superimposed (all weighted results summed) with the corresponding point of the seismic wavelet in a sliding superposition manner, generating the corresponding predicted seismic response data for each trace point by point, ultimately ensuring that the number of traces and sampling points of the predicted data are completely consistent with the actual preprocessed seismic data.

[0024] Subsequently, corresponding values ​​of the predicted response data and the actual preprocessed seismic data are extracted point by point. First, the difference between the two values ​​is calculated as predicted value - actual value. Then, each difference is squared: squaring eliminates the mutual cancellation effect of positive and negative deviations, ensuring that all deviations are represented as positive values. Next, the squared differences of all samples in the current trace are summed, and then the sum of the squared differences of all traces in the entire work area is added together. Finally, this sum is divided by the total number of samples involved in the calculation = total number of traces in the entire work area × number of sampling points per trace, yielding the global mean square error. This provides a direct and accurate quantification of the overall deviation between the predicted results and the actual data. The reflection coefficient is dynamically adjusted based on the spatial distribution characteristics of the mean square error: for areas with errors greater than 0.005, the single correction amplitude of the reflection coefficient is increased to enhance adaptability to areas with deviations; for areas with errors less than 0.001, the correction amplitude is decreased to maintain the stability of the reflection coefficient. Simultaneously, sparsity constraints are incorporated into the adjustment process. By suppressing the intensity of minor reflection coefficients with smaller amplitudes, the strong reflection characteristics of the main stratigraphic interfaces are highlighted, making the reflection coefficient distribution more closely match the true reflection patterns of the underground strata.

[0025] The process of continuous convolution calculation, error calculation, and reflection coefficient adjustment is repeated until any of the following convergence conditions are met: First, the global mean square error drops below 0.001, and the error change rate over five consecutive iterations is less than 0.5%, indicating that the prediction results have stably matched the actual data; second, the number of iterations reaches the maximum set value, and the error no longer decreases significantly, indicating convergence. The resulting initial reflection coefficient volume clearly shows the acoustic impedance differences between different stratigraphic interfaces, especially highlighting the strong reflection interface characteristics at the contact between extremely high-velocity / extremely low-velocity lithology and the underlying strata, providing a clear identification basis for subsequent extraction of strong reflection components.

[0026] Step 1.5 involves standardizing the initial reflection coefficient volume to generate the final reflection coefficient volume. This includes: performing a comprehensive data quality check on the initial reflection coefficient volume obtained after iterative convergence; using the 3-standard-deviation method to identify outlier points, i.e., removing outliers greater than the mean reflection coefficient plus 3 standard deviations or less than the mean minus 3 standard deviations, to avoid interference from extreme points in subsequent strong reflection extraction; then standardizing the initial reflection coefficient volume by mapping all reflection coefficient values ​​to the [-1, 1] interval, strictly preserving the positive and negative polarities of the reflection coefficients, with positive values ​​corresponding to peaks and negative values ​​to corresponding troughs, ensuring consistent relative strength of reflection coefficients in different regions and eliminating numerical fluctuations caused by differences in inversion parameters; subsequently, the standardized... Spatial continuity of the reflection coefficient volume was verified by calculating the correlation coefficients of reflection coefficients of adjacent traces and adjacent depth points. Regions with correlation coefficients less than 0.8 were marked as continuity fault areas. Gaussian smoothing correction was performed on these areas using a 3×3×3 spatial filter window. During the smoothing process, the filter intensity was controlled to avoid over-smoothing and destroying the true abrupt changes in the reflection coefficients of the strata. Finally, the corrected reflection coefficient volume was standardized in format and output in SEGY format consistent with the trace header information of the original seismic data to ensure accurate correspondence of trace number, depth label, and other information, which can be directly called in subsequent curve domain analysis. The final reflection coefficient volume was generated that is numerically stable, spatially continuous, clearly polarized, and can truly reflect the distribution characteristics of the strata's reflection coefficients.

[0027] In a preferred embodiment of the present invention, step 2 above may include: Step 2.1 involves performing a curvelet transform on the final reflection coefficient volume to obtain the corresponding curvelet domain coefficient volume. This includes: firstly, performing refined preprocessing on the final reflection coefficient volume; secondly, using the 3x standard deviation method to iterate through all data points, removing outlier points whose values ​​exceed the mean ± 3x standard deviation; and thirdly, filling in the missing data after the removal to ensure the spatial continuity of the reflection coefficient volume. Simultaneously, verifying that the trace number and depth sampling rate of the data are completely consistent with the original seismic data, and that the survey line number, trace number, and depth markings in the trace head information are accurately matched, avoiding errors caused by data mismatch. To mitigate transformation bias, and considering the strong reflection characteristics of the target strata in the work area, such as the strong reflection dominant frequency of 20 to 40 Hz at the coal-rock and sandstone interface in coal-bearing strata, the core parameters of the curve wave transform were determined: the scale decomposition level was set to 4 to 6 layers to ensure coverage of the frequency range corresponding to the dominant frequency of 20 to 40 Hz, while also taking into account the preservation of low-frequency weak reflection signals; the number of directional decompositions was set to 8 to 16 directions to adapt to the stratum strike of the target strata (such as 30° northeast) and the possible extension angle of the strong reflection interface, ensuring that the details of the reflection coefficient changes in different directions can be fully captured.

[0028] The curvelet transform process is initiated to convert the spatial domain reflection coefficient volume to the curvelet domain. Through multi-scale decomposition, the reflection coefficient signal is decomposed into three components: high frequency, mid frequency, and low frequency. The high frequency component corresponds to the local strong reflection abrupt signal, the mid frequency component corresponds to the continuous reflection signal of the strata, and the low frequency component corresponds to the regional background signal. Through multi-directional decomposition, the reflection characteristics at different angles are separated, so that the coefficient in each direction can accurately correspond to the extension trend of the strata interface at a certain angle. Finally, a curvelet domain coefficient volume containing complete multi-scale and multi-directional information is obtained. This coefficient volume can effectively highlight the coefficient amplitude anomalies corresponding to strong reflections, while preserving the detailed information of weak reflection signals to the greatest extent, laying the foundation for subsequent accurate capture of strong reflection interfaces.

[0029] Step 2.2: Within the corresponding curvelet transform coefficient volume, set multi-scale, multi-directional analysis windows along the target stratum. Specifically, this includes: combining the 1:50000 geological stratification report of the work area and drilling stratigraphic correlation data to accurately locate the spatial range of the target stratum within the curvelet transform coefficient volume; calibrating the depth of the top and bottom interfaces of the target stratum using drilling data to confirm that the burial depth error of the top interface is ≤2m and the burial depth error of the bottom interface is ≤2m; simultaneously determining the lateral distribution boundary of the target stratum and eliminating interference from non-target areas such as fault fracture zones; and setting corresponding multi-scale analysis windows based on the scale decomposition levels of the curvelet transform: small-scale windows are adapted to high... For fine analysis of frequency intensity reflection signals, the size is set to 20m×20m×0.5m, i.e., horizontal×horizontal×vertical, with 3 seismic trace spacings matched horizontally and 50ms transition zone at the top and bottom of the target layer vertically. For analysis of mid-frequency continuous reflection signals, the size is set to 30m×30m×1m, with 5 seismic trace spacings covered horizontally and 80ms at the top and bottom vertically. For macroscopic analysis of low-frequency background signals, the size is set to 50m×50m×2m, with 8 seismic trace spacings covered horizontally and 100ms at the top and bottom vertically, ensuring that time windows of different scales can complement each other to capture reflection characteristics.

[0030] Meanwhile, based on the stratigraphic strike (e.g., 30° NE), dip angle (e.g., 15°), and the extension direction of the strong reflection interface of the target stratum, eight multi-directional analysis time windows are set, with directions of 0° NE, 45° NE, 90° NE, 135° NE, 180° NE, 225° NE, 270° NE, and 315° NE. The extension angle of the time window in each direction is precisely matched with the directional decomposition angle of the curve transform, and the edge of the time window is parallel to the stratigraphic strike, ensuring that the time window in each direction can completely capture the distribution characteristics of the strong reflection coefficient at the corresponding angle. All time windows are set with a 50% overlap rate to avoid missing strong reflection features due to time window intervals.

[0031] Step 2.3: For each analysis window, extract the reflection coefficient sequence within the window to construct a one-dimensional reflection coefficient curve; calculate the second derivative of the one-dimensional reflection coefficient curve, and determine the inflection point of the one-dimensional reflection coefficient curve by locating the zero point of the second derivative. The inflection point is used to identify the abrupt change characteristics of the reflection coefficient. Specifically, for each set analysis window, first take the center layer of the target layer within the window as the baseline, and along the strike or dip direction of the strata, select according to the dip angle: dip angle <10° along the strike, ≥10° along the dip, extract the reflection coefficient values ​​corresponding to the curve wave domain coefficient volume point by point according to the original sampling step size of the reflection coefficient volume. For each extracted value, mark its corresponding lateral survey line number, trace number and longitudinal depth coordinate, and then arrange them in the order of spatial coordinates to form a one-dimensional reflection coefficient curve. The sampling interval of the curve is completely consistent with the lateral or longitudinal grid of the reflection coefficient volume to ensure that the curve can accurately reflect the true change trend of the reflection coefficient along the extension direction of the strata within the window.

[0032] The extracted one-dimensional reflectance coefficient curve undergoes smoothing preprocessing. First, the application rules for the 5-point moving average window are determined: the window slides along the curve point by point from left to right. For the 3rd to 3rd-to-last data point on the curve, at each target data point, five reflectance coefficient values ​​are simultaneously selected: the target point, the two points before the target point, and the two points after the target point. The arithmetic mean of these five values ​​is calculated, and this mean is used as the smoothed value for that target data point. For the first two data points on the curve, a 3-point moving average is used. For each data point, only three reflection coefficient values ​​can be selected, including the data point itself and the two points following it. For the second data point, three reflection coefficient values ​​are selected, including the data point itself, the data point before it, and the data point after it. The arithmetic mean of these two sets of values ​​is calculated as the smoothing value for the corresponding data point. For the last two data points of the curve, the same three-point moving average method is used. For the second to last data point, three values ​​are selected, including the data point itself, the data point before it, and the data point after it. For the last data point, three values ​​are selected, including the data point itself and the two points before it. The arithmetic mean of these values ​​is calculated to complete the smoothing.

[0033] After smoothing, the smoothed one-dimensional reflectance coefficient curve is compared point-by-point with the original curve. The correlation coefficient between the two curves is calculated by statistically analyzing the correlation of all corresponding point values. The correlation coefficient should be ≥0.95. If it is not reached, the number of points in the moving average window is adjusted (e.g., temporarily changed to 4 points), and the process is repeated until the requirement is met. This weakens the small fluctuations caused by random noise and avoids over-smoothing masking the abrupt changes in large values ​​corresponding to strong reflectance. After smoothing preprocessing, the second derivative of the preprocessed one-dimensional reflectance coefficient curve is calculated using the central difference method, with the sampling step size of the reflectance coefficient volume as the difference step size. First, the first derivative is calculated... The first derivative is calculated by subtracting the reflection coefficient of the adjacent point to the right from the reflection coefficient of the adjacent point to the left, and then dividing the difference by twice the sampling step size. This is the core operation of the central difference method, which differs from single-direction difference. For the first endpoint of the curve, the forward difference method is used: the reflection coefficient of the first point to the right is subtracted from the value of the endpoint, and then divided by the sampling step size to calculate the first derivative. For the last endpoint, the backward difference method is used: the reflection coefficient of the endpoint is subtracted from the value of the last adjacent point to the left, and then divided by the sampling step size to calculate the first derivative. This process forms a complete sequence of first derivatives.

[0034] Next, repeat the central difference method operation described above for the first derivative sequence to calculate the second derivative value at each point: for non-endpoint points, subtract the left adjacent value from the right adjacent value in the first derivative sequence, and then divide by twice the sampling step size; the first and last endpoints of the first derivative sequence are calculated using forward and backward difference methods respectively to obtain the complete second derivative sequence, thereby reducing the calculation error caused by single-direction difference; after completing the calculation of the second derivative sequence, traverse the sequence point by point, and for each data point, synchronously retrieve the second derivative values ​​of the one preceding and one following points. Compare the signs (positive, negative, or zero) of the second derivatives at these three points to pinpoint the location where the sign changes from positive to negative or vice versa. If a sign reversal of the second derivative is found between two adjacent sampling points (one positive and one negative), it indicates that the zero point is located between these two sampling points. Based on the second derivative values ​​and corresponding spatial coordinates of these two sampling points, the precise spatial location of the zero point, the transverse survey line number, the track number, and the longitudinal depth are gradually calculated using linear interpolation. The coordinate information of all zero points is then fully recorded.

[0035] The validity of the located zero-value points was verified. The reflection coefficient values ​​of five data points in the neighborhood of each zero-value point were extracted, and the difference between the maximum and minimum reflection coefficients in the neighborhood was calculated. Based on the statistical characteristics of the strong reflection interfaces in the work area, the following judgment criteria were set: the difference between the coal-rock and sandstone interfaces ≥ 0.2, the difference between the mudstone and sandstone interfaces ≥ 0.15, and the difference between the carbonate rock and mudstone interfaces ≥ 0.25. If the neighborhood difference reaches the corresponding criterion and the rate of change of the reflection coefficient is ≥ 0.01 / ms, it is judged as a valid inflection point. If the neighborhood difference is < 0.05 or the rate of change is < 0.01 / ms, it is judged as a false inflection point caused by noise and is removed. The locations of these valid inflection points accurately identify the abrupt change characteristics of the reflection coefficient and can be directly associated with the lateral extension trajectory and vertical burial depth of the underground strong reflection interface.

[0036] Step 2.4: For each analysis window, a local discrete surface is constructed using the 3D spatial data points within the window. This includes: for each analysis window, using the 3D spatial boundary of the window as the extraction range, the reflection coefficient values ​​of the data points in the curve domain coefficient volume are read point by point, and the corresponding 3D spatial coordinates of each data point are recorded. The horizontal coordinate is based on the unified survey network coordinate system of the work area, with the survey line number corresponding to the X-axis and the track number corresponding to the Y-axis. The vertical coordinate is based on the depth value after the two-way travel time conversion of the seismic data, ensuring that the coordinates of each data point accurately correspond to the actual geographic space of the work area. Multi-dimensional quality screening is carried out on the extracted 3D spatial data points. First, outliers with reflection coefficient values ​​exceeding the physical range of [-1,1] are removed. Then, by comparing the coordinate deviations of adjacent data points, misaligned points with coordinate offsets exceeding 0.5 grid units are removed. Finally, the signal-to-noise ratio (SNR) of each data point is calculated as the ratio of signal amplitude to noise amplitude, and low-quality points with an SNR lower than 2 are removed to avoid invalid data interfering with subsequent surface construction.

[0037] The selected valid 3D spatial data points are mapped into a 3D grid according to the grid rules of the surveying coordinate system (20m×20m×0.5m). Each data point is accurately placed into the corresponding grid node according to its coordinates. For grid nodes without data points, the average of the three valid neighboring data points is used to fill in the missing data points, ensuring the integrity of the grid nodes. Based on the grid nodes, adjacent nodes are connected according to the natural extension trend of the strata. Spatial connections between nodes are constructed first along the strike and dip direction of the strata, and then auxiliary connections in the horizontal and vertical directions are added to form a local discrete surface composed of discrete nodes and connections. The node density of this surface is consistent with the grid density of the reflectance coefficient volume, which can truly restore the 3D distribution of the reflectance coefficient within the time window.

[0038] Step 2.5 involves fitting surface patches to the local discrete surface and calculating the Gaussian curvature of each patch to form Gaussian curvature distribution data. Specifically, this includes: based on the node distribution characteristics and reflection coefficient variation gradient of the local discrete surface, first identifying abrupt change lines in the reflection coefficient within the surface. The criterion for determining an abrupt change line is that the difference in reflection coefficient between adjacent nodes is ≥0.1. Then, the surface patches are divided according to the rule of 4×4 grid units as a basic unit, with adjacent units overlapping by one grid unit. The division avoids abrupt change lines in the reflection coefficient, ensuring that there are no significant abrupt changes in the reflection coefficient within a single surface patch. Grid units in overlapping areas simultaneously belong to two adjacent surface patches, providing data support for the subsequent fitting results. For each divided surface patch, the standard deviation of the reflection coefficient of all discrete nodes within the patch is calculated to determine the fitting method: if the standard deviation of the reflection coefficient within the patch is <0.05, indicating a smooth distribution, such as in stable thick mudstone areas, a quadratic polynomial fitting method is used, with the fitting polynomial set as z=ax². The formula is: +by2+cxy+dx+ey+f (where x and y are the horizontal coordinates of the nodes within the surface patch, z is the reflection coefficient value, and a to f are the fitting coefficients to be solved). Using the coordinates of all discrete nodes within the surface patch and the reflection coefficient values ​​as constraints, the polynomial coefficients are solved using the least squares method to minimize the sum of squared residuals between the fitted surface and the discrete nodes, with a sum of squared residuals ≤ 0.01. If the standard deviation of the reflection coefficient within the patch is ≥ 0.05, indicating a significant abrupt change, such as in a strong reflection interface region, the surface patch is first divided into 2 to 3 sub-regions along the reflection coefficient abrupt change line. The reflection coefficient distribution within each sub-region is relatively uniform. Then, a linear fit is performed on each sub-region individually, with the fitting formula being z = ax + by + c. After fitting, the fitting results of adjacent sub-regions are connected using a linear transition method. The node values ​​in the transition region are taken as the weighted average of the fitted values ​​of the two sub-regions. The weights are distributed according to the distance from the node to the abrupt change line; the closer the distance, the higher the weight, avoiding numerical jumps at the sub-region connection points.

[0039] After completing the surface fitting, the Gaussian curvature of each node on each surface patch is calculated: First, the first-order partial derivatives of the fitted surface at the node are calculated, along the x-direction and y-direction of the surface. Then, the second-order partial derivatives of the node are calculated, including the degree of change of the rate of change along the x-direction, the degree of change of the rate of change along the y-direction, and the degree of change of the rate of change along the mixed x and y directions. Based on these partial derivative results, the three core parameters of the first basic form representing the local length and angular characteristics of the surface are calculated, and the three core parameters of the second basic form representing the local curvature of the surface are also calculated. Combining the parameters of the first and second basic forms, the curvature of the node in all possible directions is calculated, and the maximum and minimum curvature are selected. These two values ​​are the two principal values ​​of the node. Curvature; Multiply the two principal curvatures to obtain the Gaussian curvature value of the node. A positive value indicates that the surface is convex outward at the node, a negative value indicates that the surface is concave inward at the node, and a zero value indicates that the surface is flat at the node. Integrate the Gaussian curvature values ​​of all surface patches into continuous Gaussian curvature distribution data according to the node coordinates. First, verify the rationality of each curvature value and remove abnormal curvature values ​​with an absolute value > 1, which are mostly caused by overfitting. Then, for the curvature values ​​corresponding to surface patches with fitting residuals > 5%, replace them with the mean of the three effective curvature values ​​in the neighborhood. At the same time, ensure that the curvature value deviation in the overlapping area of ​​adjacent surface patches is ≤ 0.05. Finally, Gaussian curvature distribution data that is spatially continuous, numerically reliable, and can accurately reflect the geometric shape of the reflection coefficient volume is formed.

[0040] Step 2.6: Based on the Gaussian curvature distribution data, identify local geometric anomaly regions in the reflection coefficient volume. Specifically, this includes selecting 3 to 5 recognized normal strata areas within the work area—stable, thick mudstone areas without faults or strong reflection interfaces—each area covering ≥20 grid cells. The mean and standard deviation of the Gaussian curvature within each normal area are calculated. The average of all normal area mean values ​​is taken as the baseline mean, and the maximum standard deviation of all normal area standard deviations is taken as the baseline standard deviation. The normal threshold interval is set to the baseline mean ± 2 times the baseline standard deviation to ensure... This interval can cover more than 95% of the normal stratum curvature values. The Gaussian curvature distribution data is compared with the normal threshold interval node by node, and nodes whose curvature values ​​exceed the upper or lower limit of the interval are marked as preliminary anomalies. Spatial correlation analysis is carried out on the preliminary anomalies. The 8-neighborhood connectivity rule is adopted, that is, a node has 8 neighboring nodes above, below, left, right, front, back and diagonal. Adjacent anomalies are classified into anomaly regions. At the same time, the minimum area threshold of the anomaly region is set to be ≥3 grid cells. Isolated anomalies with insufficient area are removed, which are mostly caused by noise or fitting error.

[0041] Based on the effective inflection points located in step 2.3, the filtered abnormal areas are verified. If the abnormal area contains ≥2 effective inflection points, and the spatial distribution direction of the inflection points is consistent with the extension direction of the abnormal area, and the abrupt change in the reflection coefficient at the inflection point is ≥80% of the average abrupt change in the strong reflection interface of the work area, then it is determined to be a real local geometric abnormal area. If there are no effective inflection points or the inflection points are irregularly distributed in the abnormal area, it is marked as a suspicious area, and removed after re-verification of the fitting process. The finally determined local geometric abnormal areas accurately correspond to the distribution range of the underground strong reflection interface, providing a clear basis for the subsequent extraction of strong reflection components.

[0042] In a preferred embodiment of the present invention, step 3 above may include: Step 3.1 involves spatially fusing the inflection point locations determined within each analysis window with the identified local geometric anomaly areas to generate a comprehensive feature distribution map. This includes: firstly, conducting a unified spatial coordinate calibration, using the measured coordinates of the formation interfaces of the three key wells in the work area as a benchmark, and then performing secondary corrections on the determined effective inflection point locations and the coordinates of the identified local geometric anomaly areas to ensure that both strictly match the unified survey network coordinate system of the work area. The X-axis corresponds to the survey line number, the Y-axis corresponds to the trace number, and the Z-axis corresponds to the depth. The final coordinate deviation is controlled within 0.2 grid cells, ensuring the accuracy of spatial fusion from the root. Secondly, using the complete spatial range of the curve wave domain coefficient volume as a template, a blank comprehensive feature distribution base map is constructed. The base map not only retains the same grid density as the reflection coefficient volume but also marks basic attributes such as the survey line number range, trace number range, depth scale, and key formation interface identifiers.

[0043] Subsequently, each local geometric anomaly region was precisely plotted on the base map according to its three-dimensional boundary coordinates, and the region range was marked with a specific color block, with the region number marked at the edge of the region. Then, all effective inflection points were placed one by one to the corresponding grid nodes on the base map according to their coordinates, and marked with special symbols. At the same time, the reflection coefficient change amplitude and the corresponding stratigraphic interface type, such as coal-sandstone interface, carbonate-mudstone interface, etc., were added next to the symbols. Inflection points falling within the local geometric anomaly region and within one grid unit around it were marked with key associations. The total number of inflection points, the lateral extension density of inflection points (unit: points / 10m), the vertical distribution depth range, and the average change amplitude of inflection points in each anomaly region were counted. These statistical information were associated with the corresponding anomaly region numbering as additional attributes. Finally, the consistency of the fusion results was checked to see if there were any spatial misalignments between inflection points and anomaly regions or mismatches in attribute information. Misaligned data were recalibrated and fused again, and finally, a comprehensive feature distribution map containing the anomaly region range, inflection point spatial distribution, single-point attributes, and regional statistical characteristics was generated, clearly presenting the spatial correlation pattern of strong reflection-related features.

[0044] Step 3.2: Based on the comprehensive feature distribution map, set a feature intensity threshold to screen out candidate areas that meet the strong reflectivity characteristics. Specifically, this includes: constructing a multi-dimensional feature intensity evaluation system based on the additional statistical information of the comprehensive feature distribution map, and determining three core evaluation indicators and their corresponding calculation methods: First, inflection point density = total number of effective inflection points N in the abnormal area ÷ total area S of the abnormal area, where S is converted from the number of grid cells to the actual area, unit: m2, used to quantify the density of inflection point distribution in the abnormal area; Second, the abrupt change amplitude of the average reflection coefficient of inflection points = (amplitude of the first inflection point change + amplitude of the second inflection point change + ... + amplitude of the Nth inflection point change) ÷ total number of effective inflection points N, used to characterize the intensity of the abrupt change in reflection coefficient corresponding to the inflection point in the area; Third, the deviation of Gaussian curvature = |average Gaussian curvature of the abnormal area Kavg - mean value of normal stratum benchmark K0| ÷ standard deviation of normal stratum benchmark σ, used to quantify the degree to which the curvature of the abnormal area deviates from the range of normal stratum.

[0045] Based on measured statistical data from three typical strong reflective blocks in the work area, such as the inflection point density and abrupt change amplitude of the coal-sandstone strong reflective interface and the characteristic intensity range of normal strata, the lower thresholds for each indicator were determined through comparative analysis: inflection point density ≥ 0.3 per grid cell (equivalent to actual density ≥ 0.00075 per m²), average reflectance coefficient abrupt change amplitude ≥ 0.15, and Gaussian curvature deviation ≥ 0.2. All three indicators must be met simultaneously to be considered as meeting the characteristics of strong reflectivity. To ensure the rationality of the thresholds, two known strong reflective areas and two normal strata areas were selected for thresholding. Value verification is performed, and the threshold parameters are fine-tuned based on the verification results, with the adjustment range not exceeding ±0.02. Subsequently, the three index values ​​of all abnormal regions in the comprehensive feature distribution map are compared with the set thresholds one by one, and abnormal regions that simultaneously meet the three threshold conditions are selected as strong reflection candidate regions. Finally, the spatial continuity of the candidate regions is checked, and a minimum effective area threshold of ≥5 grid units is set, corresponding to an actual area of ​​≥2000m2. Isolated small regions with insufficient area are eliminated, and adjacent candidate regions with similar feature attributes are merged to avoid incomplete coefficient extraction in the subsequent process due to region fragmentation.

[0046] Step 3.3, in the curve wave coefficient volume, locate the curve wave coefficients corresponding to the candidate regions. Specifically, this includes: establishing a precise spatial mapping relationship between the comprehensive feature distribution map and the curve wave coefficient volume; by compiling a coordinate mapping table, determining the three-dimensional coordinates of the survey line number, trace number, and depth of each grid node in the comprehensive feature distribution map, and the one-to-one correspondence rules between these coordinates and the scale index, direction index, and spatial index in the curve wave coefficient volume, where the scale index corresponds to the 4th to 6th levels of the curve wave transform, and the direction index corresponds to the 8th to 16th directions of the decomposed direction angles; for each selected strong reflection candidate region, extract the complete spatial coordinate information of each grid node according to its three-dimensional boundary coordinates to form a candidate region coordinate list, which includes the survey line number, trace number, depth, and corresponding coordinates of each node. The corresponding region number is determined; based on the coordinate mapping table, each node in the coordinate list is mapped and transformed to accurately locate the corresponding curvesweep coefficient position in the curvesweep coefficient volume; for the boundary nodes of the candidate region, since they may cross the scale or direction boundary of the curvesweep domain, linear interpolation is used to correct the mapping deviation to ensure that the coefficient position corresponding to the boundary node is accurate. At the same time, the scale level, direction angle and original amplitude information of each located curvesweep coefficient are recorded to form a coefficient location list, in which the mid-to-high frequency scale level 3 to 5 of strong reflection signal concentration and the direction angle consistent with the stratum strike are marked to provide clear guidance for subsequent targeted coefficient extraction and ensure that no curvesweep coefficient related to strong reflection in the candidate region is missed.

[0047] Step 3.4: Extract the located curvesweep domain coefficients to form strong reflection coefficient components. Specifically, this includes: based on the generated coefficient location list, extracting curvesweep domain coefficients corresponding to candidate regions in descending scale level and clockwise order of direction angle within the same scale; during extraction, setting differentiated validity verification standards for different scale levels: for mid-to-high frequency scales (layers 3 to 5), where strong reflection signal characteristics are obvious, eliminating small coefficients with amplitudes less than 0.08; for low-frequency scales (layers 1 to 2), where background signals are dominant, only retaining coefficients with amplitudes greater than 0.03 and closely related to the spatial structure of the strong reflection region, avoiding the extraction of noise coefficients or irrelevant background coefficients; for valid coefficients that pass verification, archiving them in an orderly manner according to their original scale level, direction angle, and spatial coordinates to establish a scale... A three-dimensional indexing system of direction and coordinates is used to ensure that the spatial distribution structure of the coefficients is completely consistent with the structure of the candidate region in the curve domain coefficient volume. Subsequently, the completeness of the archived effective coefficients is checked by comparing the coefficient location list with the actual number of effective coefficients extracted. If the missing rate exceeds 5%, the missing coefficients are relocated and extracted. Finally, all archived effective coefficients are integrated and summarized to generate an independent strong reflection coefficient component. At the same time, a component metadata report is compiled, which clarifies the key information such as the extraction area coordinate range of the strong reflection coefficient component, the scale level and direction angle included, the total number of effective coefficients and the average amplitude. This component accurately focuses on the curve domain coefficients related to strong reflection, laying a solid foundation for the subsequent accurate removal of strong reflection components and the maximum preservation of weak reflection effective signals.

[0048] In a preferred embodiment of the present invention, step 4 above may include: Step 4.1 involves three-dimensional spatial sampling of the strong reflection coefficient components to establish a discrete data matrix of the spatial distribution of these components. Specifically, this includes: first, using the measured coordinates of the strong reflection interfaces of three key wells in the work area as a reference, performing a secondary three-dimensional spatial coordinate calibration on the generated strong reflection coefficient components. This is achieved by calculating the deviation between the coordinates of each coefficient and the corresponding measured coordinates of the wells, and adjusting coordinate points with deviations exceeding 0.2 grid cells. Ultimately, this ensures that the spatial position of all coefficients is completely aligned with the original seismic data and the X-axis survey line number, Y-axis trace number, and Z-axis depth of the reflection coefficient volume's coordinate system, with a coordinate deviation ≤ 0.1 grid cells. The grid unit; combining the energy concentration characteristics of the strong reflection coefficient component (mainly mid-to-high frequency scale) and the spatial scale corresponding to the main frequency of the strong reflection signal in the work area, 20 to 40 Hz corresponds to a horizontal resolution of 20m and a vertical resolution of 0.5m, the three-dimensional spatial sampling parameters are determined: the sampling interval strictly matches the grid interval of the final reflection coefficient volume, with 20m in the horizontal X and Y directions and 0.5m in the vertical Z direction. The sampling range covers the entire spatial area of ​​the strong reflection coefficient component, and at the same time, it extends outward by 1 sampling interval in each of the X, Y, and Z directions to form the sampling range of the core area + transition area, so as to avoid missing weak energy information at the edge.

[0049] According to the set sampling interval, the strong reflection coefficient components are sampled in three-dimensional space in a row-by-row, column-by-column, and depth-by-depth manner. The three-dimensional coordinates (X, Y, Z), corresponding strong reflection energy amplitude, and scale level of each sampling point are recorded simultaneously. Referring to the three-dimensional indexing system in step 3.4, an initial set of sampling points is formed. Multi-dimensional quality screening is then performed on the initial set of sampling points: first, low-energy invalid points with energy amplitude < 0.03 are removed. This threshold is based on the energy amplitude statistics of the normal background area to ensure noise removal; then, through neighborhood coordinate comparison, points with coordinate offsets exceeding 0 are removed. Outliers were identified at 5 sampling intervals. For missing points encountered during sampling, such as those not covered in the edge transition zone, the data was supplemented using a weighted average of three neighboring valid sampling points. Weights were distributed inversely proportional to distance, with closer points having higher weights, and the total weights equal to 1. Finally, the integrity of the supplemented discrete data matrix was verified, requiring a data integrity rate ≥99% and a missing point ratio ≤1%. If these standards were not met, resampling and supplementation were performed. The final result was a spatially continuous, complete, and accurately energy-informed discrete data matrix that perfectly matched the spatial distribution of the strong reflection coefficient components. Step 4.2: Based on the established discrete data point matrix, construct the wavefield energy topology skeleton reflecting the spatial structure of strong reflection energy by calculating the spatial correlation and energy coupling relationship between points. Specifically, this includes: First, based on the established discrete data point matrix, determine the spatial correlation analysis range between points: Use a 3×3×3 three-dimensional cubic neighborhood window, that is, each sampling point is only correlated with sampling points adjacent to it in 8 directions (up and down, left and right, front and back, and diagonal). The size of this window is adapted to the local aggregation scale of strong reflection energy to avoid introducing irrelevant background points due to an excessively large correlation range; Calculate the core parameters of each sampling point and its effective neighboring sampling points: First, the spatial distance, calculated as the straight-line distance based on the coordinate difference of the measurement network coordinate system, with the unit uniformly in meters; Second, the energy coupling coefficient, calculated by multiplying the energy amplitude of two points by the square of the spatial distance. The smaller the distance and the larger the amplitude, the larger the coupling coefficient, accurately quantifying the tightness of the energy correlation between two points.

[0050] By statistically analyzing the coupling coefficient distribution of sampling points with no strong reflection in the normal background area of ​​the work area, a screening threshold of 0.6 was determined. Pairs of adjacent sampling points with a coupling coefficient > 0.6 were selected and considered as effective correlation links for strong reflection energy. Based on these effective correlation links, a topological framework was constructed using a process of clustering (such as density clustering DBSCAN, commonly used in this field), extraction, and connection: First, discrete sampling points were connected through links to form multiple local energy clusters. The total energy amplitude of each cluster was calculated, and clusters with a total energy > 3 times the background mean were retained. Then, the energy high-value centers of each effective cluster were extracted, with the high-value center determined by an amplitude > region mean + 1 standard deviation. This line served as the backbone of the wavefield energy topological framework, and the backbone line had to traverse the clusters. Core Range: For branch energy regions around the main trunk where the energy amplitude is greater than 1.5 times the background mean, secondary energy high-value lines are extracted as skeleton branches. The branches are connected to the main trunk using smooth curves to ensure the connectivity of the topology. Finally, the topology skeleton is optimized by smoothing: a 5-point moving average window is used to adjust the coordinate position of the skeleton nodes, and isolated branches with a length of less than 3 sampling points are removed, as these branches are mostly caused by noise. The connection angle between nodes is corrected to ensure that the angle deviation is less than 30° to avoid abrupt broken lines. After optimization, the consistency of the energy distribution of the skeleton and the discrete data point array is checked. It is required that the energy covered by the skeleton accounts for more than 80% of the total energy of the clustered region. Finally, a wave field energy topology skeleton that can accurately reflect the spatial extension trend of strong reflective energy clusters, with clear connectivity and a smooth structure is formed.

[0051] Step 4.3: Based on energy intensity and topological location, define several feature control points located inside, at the edge of, and outside of strong reflection energy clusters within the wavefield energy topology framework. Specifically, this includes: first, statistically analyzing the strong reflection energy amplitudes corresponding to the wavefield energy topology framework, calculating the mean and standard deviation of energy amplitudes for all sampling points, and determining the energy intensity classification criteria: regions with energy amplitudes greater than the mean + 1.5 times the standard deviation are considered inside strong reflection energy clusters; regions with energy amplitudes between the mean - 0.5 times the standard deviation and the mean + 1.5 times the standard deviation are considered at the edge of energy clusters; and regions with energy amplitudes less than the mean - 0.5 times the standard deviation are considered outside energy clusters. Based on this classification criterion and the topological location of the wavefield energy topology framework, define three types of feature control points: control points inside strong reflection energy clusters, selecting sampling points from the high-value center region of each energy cluster, requiring each inner... At least three evenly distributed points are selected in the primary region, with a distance of no less than two sampling intervals between them. For energy cluster edge control points, sampling points at locations of abrupt changes in energy amplitude gradient (gradient change rate ≥ 0.02 / sampling interval) are selected, evenly distributed along the edge, with a distance of no more than three sampling intervals between adjacent edge points. For energy cluster external control points, sampling points located far from the energy cluster and in the background region are selected, ensuring that each energy cluster has at least two external control points, with a distance of no less than five sampling intervals from the energy cluster edge. All selected feature control points are labeled with attributes to clarify their type (internal, edge, external), corresponding energy amplitude, and topological location information. The distribution of control points is then checked to ensure that each strongly reflective energy cluster is completely surrounded and covered by all three types of control points, ultimately forming a set of feature control points containing spatial coordinates, energy attributes, and topological information.

[0052] Step 4.4: Based on the aforementioned feature control points and their spatial constraints, an implicit energy isomorphic manifold representing the spatial distribution of strong reflection energy is generated through implicit surface fitting calculations. Specifically, this includes: constructing a dual-constraint mechanism centered on the defined feature control points: a hard spatial constraint requires the fitted surface to pass through or be close to all feature control points; the constraint conditions can be expressed as follows: (Internal or edge control points) or (External control points), among which Let i be the three-dimensional coordinates of the i-th feature control point. For implicit surface functions, This refers to the isosurface to be fitted; soft constraints on energy properties ensure that the energy properties of the surface at the control points match the actual values, and the constraint conditions are as follows: ,in This is a function for calculating the energy amplitude of the surface at the corresponding location. Let be the actual strong reflection energy amplitude at the i-th control point. A threshold for allowable deviation of energy amplitude is set; simultaneously, an extension trend constraint is introduced into the wavefield energy topology skeleton, through the direction vectors of the skeleton nodes. Restrict the surface extension direction, that is, the normal vector of the surface at the skeleton node and the included angle to ensure that the extension direction of the fitted surface is consistent with the skeleton.

[0053] Combined with the statistical characteristics of the strong reflection energy clusters in the work area, select the average energy amplitude of the edge control points of the energy cluster as the threshold of the fitted energy isosurface. The calculation formula is: ; where is the total number of edge control points, is the energy amplitude of the i-th edge control point. This threshold needs to be verified by the edges of 2 known strong reflection regions to ensure that the corresponding isosurface can accurately define the boundary of the strong reflection energy cluster, neither including the background low-energy region, and the background energy amplitude is usually <TE−0.03, nor missing the weak energy details at the edge.

[0054] Adopt the radial basis function (RBF) fitting to construct the initial implicit surface relationship. First, select the multiquadric radial basis function that adapts to the three-dimensional irregular distribution characteristics of the strong reflection energy. The expression is , where r is the straight-line distance from the spatial point to the control point, c is the shape parameter of the basis function, and the value is c = 0.5×Δh, where Δh is the three-dimensional space sampling interval; then construct an independent basis element centered on each characteristic control point, and the influence range of the basis element is dynamically adjusted according to the distribution density of the control points: the control points in the core area are dense, ; the control points in the edge or outer area are sparse, ; finally, the initial implicit surface function is obtained by linearly combining all basis elements: ; where N is the total number of characteristic control points, is the weight coefficient of the i-th basis element, is the linear polynomial term, which is used to ensure the overall continuity of the surface. By solving the system of equations (j = 1, 2,..., N) to determine and the values of a, b, c, d, and obtain the initial implicit surface that initially fits the spatial positions of all control points.

[0055] Define the residual function with the goal of minimizing the spatial distance residual between the surface and the control points: ; where is the implicit surface function in the k-th iteration. In each iteration, adjust the weight coefficient of the basis element and the polynomial term parameters a, b, c, d through the gradient descent method, and update the surface function ; the iteration stop condition is R < 0.5Δh and Δh is the three-dimensional spatial sampling interval, ensuring that the spatial distance deviation between the surface and the control point is within the allowable range; for surface wrinkles or discontinuous regions that appear during the iteration process, a 3×3×3 local smoothing window is used for correction. For any surface node x within the window, the current node value is updated by weighted averaging of the surface function values ​​of neighboring nodes. ; Where Ω represents the set of neighboring nodes within the window, and M represents the number of neighboring nodes. The weights are inversely proportional to the distance; the closer the distance, the higher the weight. This operation smooths the surface transition and eliminates local abrupt defects.

[0056] After fitting, multi-dimensional validity verification is performed: spatial distance verification calculates the spatial distance between the implicit energy isomorphic manifold, i.e., S(x)=0, and all feature control points. , ( For curved surfaces (normal vector at the location), requiring more than 95% of control points to satisfy... Boundary integrity verification involves checking whether the length of the energy cluster edge covered by the isomanifold (Lcover / Ltotal) is greater than or equal to 98% (Lcover is the coverage length, and Ltotal is the total edge length). This verifies whether the isomanifold completely encloses the strongly reflective energy cluster and whether it misses details such as edge transitions and branches. Energy property verification involves extracting M uniformly distributed sampling points on the isomanifold. Calculate its energy amplitude The deviation from the set threshold is required to meet the following conditions. For areas that fail verification, such as missing edge details or excessive energy deviation, the reasons are analyzed and targeted optimizations are made: if the problem is insufficient control points, edge control points are added to the missing areas, the Nedge increases after addition, TE is recalculated and fitted; if the problem is an unreasonable threshold, TE is fine-tuned (fine-tuning amplitude ±0.01) and then refitted; after multiple optimizations and verifications, a clear boundary and spatially continuous implicit energy isomorphic manifold S(x)=0 is finally generated. This manifold can accurately characterize the spatial distribution of strong reflection energy, including the core area range, edge contour, and extension direction, providing a reliable morphological basis for the accurate removal of strong reflection components in the future.

[0057] In a preferred embodiment of the present invention, step 5 above may include: Step 5.1 involves performing geometric analysis on the implicit energy equivalent manifold to extract morphological attributes, including the manifold's surface area, volume, and curvature distribution characteristics. Specifically, this includes: first, verifying the geometric integrity and validity of the implicit energy equivalent manifold by using a 3×3×3 three-dimensional neighborhood search window to traverse all grid nodes of the manifold and checking for morphological defects such as missing nodes causing fractures, overlapping nodes, or isolated small patches; for fractured areas, performing linear interpolation repair based on the node coordinates and morphological trends on both sides of the fracture; deduplicating and merging overlapping nodes; directly removing isolated patches with an area less than 3 grid cells to ensure that the repaired manifold has good spatial continuity and closure; and then, based on the manifold's three-dimensional spatial coordinate information, performing grid discretization to convert the continuous implicit surface into a discrete three-dimensional grid configuration. The grid cell size remains consistent with the sampling interval in Step 4.1, 20m horizontally and 0.5m vertically, ensuring that the original morphology of the manifold can be accurately restored after discretization.

[0058] The discretized 3D mesh model is further decomposed using a triangular patch method. Each triangular patch consists of three adjacent mesh nodes. During the decomposition process, the difference in side lengths of the patches is strictly controlled to ensure that the side length deviation between adjacent patches does not exceed 50%, and adjacent patches share a common edge to avoid gaps or overlaps. This ultimately forms a continuous set of triangular patches covering the entire manifold. The manifold surface area is calculated for the decomposed set of triangular patches. First, the actual area of ​​each triangular patch is calculated using the coordinates of its three vertices. The area error caused by planar projection is corrected by incorporating the coordinate projection type of the work area, such as UTM projection. Finally, the areas of all triangular patches are summed one by one. The total surface area of ​​the implicit energy equivalent manifold is obtained. Before accumulation, abnormally small patches with an area less than 0.1 m2 need to be removed to avoid affecting the calculation accuracy. When calculating the manifold volume, the minimum bounding box of the manifold is first constructed. By traversing the three-dimensional coordinates of all nodes of the manifold, the maximum and minimum X, Y, and Z coordinates of the bounding box are determined to form a cubic bounding box that can completely enclose the manifold. Then, the bounding box is subdivided into several small cubic meshes with the same sampling interval. The number of small cubes that are completely located inside the manifold is counted. The preliminary volume of the manifold-enclosed region is calculated by combining the volume of the small cubes. For small cubes located at the manifold boundary, the volume value is corrected by judging the proportion of manifold nodes inside them to improve the accuracy of volume calculation.

[0059] When extracting curvature distribution features, the three vertices of each triangular facet are first locally smoothed to avoid individual abnormal nodes affecting the curvature calculation results. Then, the Gaussian curvature and mean curvature of each vertex are calculated separately. Gaussian curvature is used to characterize the overall bending trend of the surface at the vertex, with positive values ​​for convex regions, negative values ​​for concave regions, and near zero values ​​for flat regions. Mean curvature is used to characterize the local bending degree of the surface at the vertex, with larger values ​​indicating more pronounced bending. The curvature values ​​of all vertices are collected, and a curvature distribution histogram is generated using statistical analysis methods. Core statistical indicators such as mean curvature, maximum curvature, minimum curvature, and standard deviation of curvature are extracted from the histogram. At the same time, a criterion for judging high curvature regions is set, namely, the absolute value of curvature is greater than the mean curvature plus 1.5 times the standard deviation. The area of ​​the triangular facets contained in this region is counted, and its proportion of the total surface area of ​​the manifold is calculated. Finally, the surface area, volume, curvature statistical indicators, and the proportion of high curvature regions of the manifold are integrated to form a complete set of morphological attributes.

[0060] Step 5.2 involves performing energy field analysis on the implicit energy isomorphic manifold to calculate the energy gradient distribution attributes inside and on the surface. Specifically, this includes: first, determining the analysis range of the energy field, using the implicit energy isomorphic manifold as the boundary; the internal analysis range is the core region of strong reflection energy enclosed by the manifold, extending externally to the background region three sampling intervals outside the manifold, ensuring complete coverage of the influence range of strong reflection energy; based on the established discrete data point matrix, extracting the strong reflection energy amplitude of all sampling points within the analysis range to construct three-dimensional energy field data matching the spatial range of the manifold; when calculating the internal energy gradient distribution attributes, analysis points are selected inside the manifold at densities of two sampling intervals, and the central difference method is used to calculate the energy gradient distribution attributes of each analysis point. The energy gradient components in the X, Y, and Z directions are integrated to obtain a three-dimensional energy gradient vector. The average gradient magnitude, maximum gradient value, and gradient direction distribution of all internal analysis points are statistically analyzed to determine the main direction of the energy gradient. When calculating the surface energy gradient distribution attributes, sampling points are uniformly selected along the surface of the implicit energy isomorphic manifold, with at least one sampling point selected for each triangular facet. The energy gradient of each surface sampling point along the tangent direction of the manifold is calculated. The amplitude distribution range of the surface gradient, the location and extension length of gradient abrupt change regions (gradient amplitude greater than the mean + 1 standard deviation) are statistically analyzed. Finally, an energy gradient distribution attribute set containing the statistical characteristics of the internal energy gradient and the distribution characteristics of the surface energy gradient is formed.

[0061] Step 5.3: Combine the extracted morphological attributes and the calculated energy gradient distribution attributes to construct a multi-dimensional feature vector. This includes: firstly, standardizing the extracted morphological attributes and the calculated energy gradient distribution attributes to eliminate dimensional differences between different attributes. The standardization method uses min-max normalization to map all attribute values ​​to the [0, 1] interval, ensuring balanced weighting of each attribute in subsequent feature fusion; secondly, identifying the core indicators for both types of attributes. For morphological attributes, four core indicators are selected: normalized surface area, normalized volume, normalized mean curvature, and normalized proportion of high curvature regions. For energy gradient distribution attributes, the normalized mean internal gradient and the normalized maximum internal gradient are selected. Four core indicators are identified: normalized value, normalized value of surface gradient magnitude range, and normalized value of gradient abrupt change region length. Following a fixed order of morphological attribute indicators first, followed by energy gradient attribute indicators, these eight standardized core indicators are arranged sequentially to construct an 8-dimensional feature vector. Each dimension corresponds to a specific attribute indicator, and each element of the vector is the standardized value of the corresponding indicator. The constructed multidimensional feature vectors are validated for validity, and invalid vectors caused by data anomalies are removed. If the value of a certain dimension exceeds the [0, 1] interval, the attributes of the manifold region corresponding to the invalid vector are re-extracted and a vector is constructed. Finally, a set of standardized multidimensional feature vectors corresponding one-to-one with each implicit energy equivalent manifold is formed.

[0062] Step 5.4: Based on the constructed multidimensional feature vector, the strong reflection stripping optimization coefficient is calculated through the set optimization mapping function. Specifically, this includes: first, combining the engineering requirements for strong reflection stripping in the work area and historical measured data, determining the core objective of the optimization mapping function, namely, accurately mapping the multidimensional feature vector to obtain the strong reflection stripping optimization coefficient. This coefficient needs to balance the thoroughness of strong reflection stripping and the fidelity of weak reflection signals, with a value range set to [0.6, 1.0]. A larger coefficient indicates higher stripping strength. The optimization mapping function is constructed using Support Vector Regression (SVR), and its mathematical... The expression is: f(V)=ω·ϕ(V)+b; where f(V) is the strong reflection stripping optimization coefficient, which is the function output value, and the value range is [0.6, 1.0]; V=[v1,v2,...,v8] is the constructed 8-dimensional normalized multidimensional feature vector, v1 to v4 are the normalized values ​​of surface area, volume, mean curvature, and proportion of high curvature region, respectively, and v5 to v8 are the normalized values ​​of mean internal gradient, maximum internal gradient, range of surface gradient magnitude, and length of gradient abrupt region, respectively. ϕ(·) is the kernel function, and the radial basis function (RBF) is selected to map low-dimensional feature vectors to high-dimensional feature space, solving the nonlinear mapping problem. ω is the weight vector in the high-dimensional feature space, used to measure the contribution of high-dimensional features to the stripping optimization coefficient. b is the bias term, used to adjust the overall output benchmark of the function to ensure that the mapping result fits the actual stripping requirements of the work area. A training sample set is constructed based on the measured data of three known strong reflection blocks in the work area. Each sample contains a set of 8-dimensional standardized multi-dimensional feature vectors Vi and the corresponding optimal stripping coefficient yi. The optimal stripping coefficient is determined by multiple field stripping experiments and must meet the dual requirements of strong reflection residue <5% and weak reflection signal amplitude retention rate ≥90% after stripping. Based on the feature vectors and optimal stripping coefficients of the training sample set, the weight vector ω and bias term b of the optimized mapping function are determined by statistical fitting. The fitting objective is to minimize the mean square error between the function's predicted stripping coefficient and the measured optimal coefficient, and finally ensure that the average absolute deviation between the function's predicted value and the measured value is less than 0.03.

[0063] Subsequently, the performance of the constructed mapping function was verified. Two strong reflection blocks that were not involved in the fitting were selected as verification samples. Their multidimensional feature vectors were input into the function to calculate the predicted stripping coefficients, and the rationality of the predicted coefficients was verified through actual stripping experiments. If the experimental results did not meet the standard of strong reflection residue <5% and weak reflection retention rate ≥90%, the parameters of the radial basis kernel function were readjusted and fitted again. If the experimental results met the standard, the mapping function was deemed to be verified as qualified. All the constructed multidimensional feature vectors were input one by one into the verified optimized mapping function to calculate the strong reflection stripping optimization coefficients corresponding to each implicit energy isomorphic manifold. At the same time, the contribution of each feature dimension to the final coefficients was decomposed based on the weight vector ω. The calculation method is the ratio of the weight value of that dimension to the sum of the absolute values ​​of all weight values. Finally, a complete set containing the stripping optimization coefficients, the contribution ratio of each dimension, and the accuracy verification results was formed, providing a quantitative and interpretable basis for the subsequent accurate stripping of strong reflection components.

[0064] In a preferred embodiment of the present invention, step 6 above may include: Step 6.1 involves using a strong reflection stripping optimization coefficient to weight and adjust the strong reflection coefficient components, obtaining the weighted strong reflection coefficient components. Specifically, this includes: first, performing precise spatial coordinate matching between the strong reflection stripping optimization coefficient and the strong reflection coefficient components. Using the unified measurement network coordinate system of the work area as a reference, verifying the stripping optimization coefficient corresponding to each implicit energy isomorphic manifold to ensure that its coverage is completely aligned with the three-dimensional coordinates of the corresponding strong reflection region in the strong reflection coefficient component, with coordinate deviation controlled within 0.2 grid cells; based on the matching results, adjusting the strong reflection coefficient components using a grid-by-grid node weighting method. The core weighting adjustment formula is: Weighted coefficient value of the core area of ​​the strong reflection region = Original coefficient value of the strong reflection coefficient component × Strong reflection stripping optimization coefficient; where the original coefficient value of the strong reflection coefficient component is the extracted coefficient of the corresponding grid node. The strong reflection stripping optimization coefficient is the quantization coefficient of the corresponding region calculated. For the edge transition zone of the strong reflection region, within 1 to 2 grid cells outside the implicit energy isomorphic manifold, a gradual weighting strategy is adopted. The corresponding weighting formula is: weighted coefficient value of the transition zone = original coefficient value of the strong reflection coefficient component × (strong reflection stripping optimization coefficient × gradual attenuation coefficient). The gradual attenuation coefficient gradually decreases linearly from 1 at the manifold boundary to 0 at the edge of the transition zone to avoid obvious numerical abrupt changes at the edge after weighting. After the weighting adjustment is completed, the validity of the results is verified. Abnormal nodes with an absolute value greater than 1 after weighting are removed, as they exceed the physical meaning range of the reflection coefficient. Abnormal nodes are replaced by the mean of the three effective weighted nodes in the neighborhood. Finally, the weighted strong reflection coefficient component is obtained that is spatially continuous, numerically reasonable, and accurately matches the strong reflection region.

[0065] Step 6.2 involves subtracting the weighted strong reflection coefficient components from their corresponding positions in the reflection coefficient volume to generate an intermediate reflection coefficient volume. This includes: establishing a spatial mapping relationship between the weighted strong reflection coefficient components and the original reflection coefficient volume (i.e., the final generated reflection coefficient volume); confirming the 3D coordinate position of the original reflection coefficient volume corresponding to each grid node in the weighted components using a coordinate index table to ensure complete spatial grid correspondence between the two; and traversing all grid nodes within the strong reflection region of the original reflection coefficient volume in a line-by-line, track-by-track, and depth-by-depth sequence, subtracting the corresponding weighted strong reflection coefficient from the original reflection coefficient value of each node. The strong reflection coefficient value is used to obtain the intermediate reflection coefficient value of the node. For nodes in the original reflection coefficient volume that do not belong to the strong reflection region, their original reflection coefficient values ​​are directly retained without subtraction. After the subtraction operation is completed, the integrity of the generated intermediate reflection coefficient volume is checked to ensure that all nodes in the strong reflection region have completed the subtraction process without omission or duplicate operation. At the same time, the rationality of the intermediate reflection coefficient value is checked, and values ​​exceeding the range of [-1, 1] are truncated to the boundary of the interval to avoid abnormal values ​​that do not conform to physical meaning. Finally, an intermediate reflection coefficient volume is formed that retains the weak reflection signal and initially removes the strong reflection.

[0066] Step 6.3 involves energy equalization processing of the intermediate reflectance coefficient volume to obtain the reflectance coefficient volume after removing the strong axis. Specifically, this includes: statistically analyzing the energy distribution characteristics of the intermediate reflectance coefficient volume; dividing the work area into units according to the survey lines; calculating the mean, maximum, and standard deviation of the energy amplitude of the reflectance coefficient within each survey line unit; identifying areas with uneven energy distribution and areas where the energy amplitude deviates from the overall mean ±0.1; and using a regional adaptive energy equalization method. For areas with excessively high energy (amplitude greater than the mean +0.1), the gain is adjusted through linear attenuation to compress the energy amplitude within the area to near the overall mean. For areas with excessively low energy (amplitude less than the mean -0.1), the gain is adjusted to compress the energy amplitude within the area to near the overall mean. 1. By linearly enhancing the gain adjustment, the energy amplitude of the weak reflection signal within the region is increased, ensuring that the energy levels of different regions tend to be balanced. After equalization, a 5-point moving average window is used to perform global smoothing on the intermediate reflection coefficient volume, weakening the local energy abrupt changes that may occur during equalization, while preserving the detailed features of the weak reflection signal. Finally, the processing results are verified by comparing the energy distribution histograms before and after equalization to ensure that the overall energy distribution is smoother, and the amplitude retention rate of the weak reflection signal is ≥90% and the residual amount of strong reflection is <5%. The final result is a reflection coefficient volume after removing the strong axis, with complete removal of the strong axis, high fidelity of the weak reflection signal, and balanced energy distribution.

[0067] In a preferred embodiment of the present invention, step 7 above may include: Step 7.1 involves extracting the seismic wavelet used in generating the reflection coefficient volume. Specifically, this includes: first, locating the original seismic data gather upon which the reflection coefficient volume is generated; selecting continuous data segments from the gather that are free from strong reflection interference, have a high signal-to-noise ratio, and exhibit stable strata as the wavelet extraction window; setting the window length to 2 to 3 times the dominant frequency period of the seismic wavelet to ensure that the extracted wavelet possesses complete phase and amplitude characteristics; and then processing the seismic data within the extraction window using multichannel statistical correlation analysis. By analyzing the degree of similarity and correlation between different data channels, the correlation between the multichannel data is determined. The combined results are averaged to suppress random noise and obtain the initial seismic wavelet. Subsequently, the initial seismic wavelet is subjected to phase correction and amplitude normalization to correct the phase distortion problem of the wavelet, so that the phase spectrum of the wavelet tends to zero phase or constant phase. At the same time, the maximum amplitude of the wavelet is normalized to 1 to ensure the consistency of the wavelet energy. Finally, the amplitude spectrum and phase spectrum characteristics of the wavelet are verified to ensure that they match the dominant frequency and phase characteristics of the seismic data in the work area, so as to avoid the influence of factor wave deviation on the subsequent reconstruction accuracy. Finally, a standard seismic wavelet that is completely consistent with the wavelet used to generate the reflection coefficient volume is obtained.

[0068] Step 7.2 involves performing a convolution operation between the extracted seismic wavelet and the reflection coefficient volume after removing the strong axis to generate the initial reconstructed seismic data volume. This includes: conducting comprehensive spatial and dimensional matching preparation between the standard seismic wavelet and the reflection coefficient volume after removing the strong axis; establishing a coordinate index comparison table between the two based on the unified coordinate system of the work area; clarifying the one-to-one correspondence between the three-dimensional coordinate survey line number, trace number, and depth of each grid node in the reflection coefficient volume after removing the strong axis and the time sampling points of the standard seismic wavelet; considering that the sampling dimension of the reflection coefficient volume is depth and the sampling of the seismic wavelet... The dimension is time, requiring precise conversion between depth and time. During the conversion process, the conversion error is corrected by combining the measured velocity parameters of each depth layer to ensure that the depth sampling interval of the reflection coefficient volume (0.5m) and the time sampling interval of the wavelet (e.g., 2ms) are perfectly matched, and the coordinate deviation after conversion is ≤0.1 grid cells. At the same time, the standard seismic wavelet is preprocessed to confirm the phase characteristics of the wavelet, maintain the zero phase or constant phase after the correction in step 7.1, and truncate the effective length of the wavelet, usually 3 to 5 main frequency cycles, and remove redundant sampling points with amplitudes close to zero at both ends to improve computational efficiency.

[0069] Subsequently, convolution operations were performed in the following order: first grouped by survey line, then sorted by track number, and finally progressively along the depth direction. During the operation, a single grid node of the reflection coefficient volume after removing the strong axis was used as the core. The center point of the standard seismic wavelet was aligned with this node. Then, using the effective length of the wavelet as the range, the amplitude of each sampling point of the wavelet was multiplied one by one with the coefficient value of the corresponding grid node of the reflection coefficient volume within that range. All product results were summed to obtain the convolution result value corresponding to that node. For boundary grid nodes of the reflection coefficient volume, such as the beginning and end of the survey line, the edge of the track number, and the upper and lower depths, a strategy combining wavelet truncation and amplitude attenuation was adopted: when the wavelet exceeded the range of the reflection coefficient volume, only the sampling points within the range were retained for the operation. At the same time, the amplitude of the wavelet sampling points near the boundary was linearly attenuated, with the attenuation coefficient gradually decreasing from 1 to 0.2 to avoid waveform distortion caused by boundary effects. Throughout the operation, the phase characteristics of the standard seismic wavelet were strictly maintained. By monitoring the wavelet phase spectrum and the phase value of the nodes after the operation in real time, it was ensured that there was no phase shift or phase reversal problem in the reflection phase axis.

[0070] After the calculation is completed, the initial reconstructed seismic data volume is first checked for dimensional integrity. The number of survey lines, traces, and depth time sampling points are compared with the original seismic data volume to ensure that there is no missing data, dimensional misalignment, or duplicate calculation. Then, the waveform continuity of the data volume is checked by visualization comparison. The focus is on checking whether the weak reflection phase axis in the strong reflection removal area is intact and without obvious breaks or abrupt changes. Finally, the phase accuracy is checked. Seismic traces corresponding to three key wells in the work area are selected, and the reflection phase of the reconstructed data and the original data at the target layer is compared to ensure that the phase deviation is ≤5°. Finally, an initial reconstructed seismic data volume with continuous waveform, accurate phase, precise spatial location, and complete preservation of weak reflection characteristics is generated.

[0071] Step 7.3 involves amplitude compensation and residual noise suppression on the initial reconstructed seismic data volume to obtain the final seismic data after strong axis removal, thus representing the completion of strong axis removal. Specifically, this includes: conducting full-range energy distribution statistics on the initial reconstructed seismic data volume; calculating the mean, maximum, and standard deviation of reflection amplitude at each depth level (each 50m is a layer unit); determining the amplitude attenuation law at different depth levels; and highlighting deep weak reflection areas (amplitude less than the overall mean -0.1) and shallow areas with excessively high amplitude (amplitude greater than the overall mean +0.1); and employing a regional adaptive amplitude compensation method to process the data. For deep, weakly reflective regions, the gain is adjusted by linear enhancement channel by channel: using the average normal amplitude of the adjacent upper layer as a benchmark, the gain intensity is gradually increased according to the depth increasing trend to ensure that the amplitude of deep, weakly reflective regions after enhancement is close to the overall average, while avoiding excessive enhancement and amplification of noise; for shallow regions with excessively high amplitude, a moderate linear compression strategy is adopted to compress the excessively high amplitude proportionally to near the overall average, while preserving the amplitude difference characteristics of effective reflection; amplitude changes are monitored in real time during the compensation process, and the amplitude range is checked every 10 channels of data to ensure that the amplitude of all regions after compensation is within a reasonable range and there are no abnormal amplitudes that exceed physical meaning.

[0072] Subsequently, residual noise suppression was performed using a multi-channel coherent filtering method. The specific procedure was as follows: a 3×3 multi-channel neighborhood window was selected, covering 3 channels horizontally and 3 sampling points vertically, traversing all channels of the initial reconstructed seismic data volume; for the seismic data within each window, the waveform coherence between channels was calculated. By comparing the continuity and amplitude changes of the reflection phase axis of adjacent channels, a coherence threshold of 0.7 was set, retaining effective reflection signals with coherence higher than the threshold, and suppressing incoherent residual noise; during filtering, the window size and coherence threshold were strictly controlled to avoid blurring weak reflection details due to an excessively large window, or missing effective signals due to an excessively high threshold; for strong reflections... In addition to the edge transition zone of the region, local weakening filtering is used to further suppress any residual noise in the transition zone, while protecting the continuity of the weak reflection phase axis. After noise suppression, global smoothing optimization is performed, using a 5-point moving average window (2 points horizontally + the current point + 2 points vertically) to smooth the data volume point by point: the amplitude values ​​of each point within the window are averaged as the smoothed amplitude value of the current point. This operation weakens local amplitude abrupt changes and improves the waveform smoothness of the data volume. During the smoothing process, the location of the weak reflection phase axis is given special protection. If the window contains weak reflection signals, the weight of the window edge points is appropriately reduced to ensure that the weak reflection details are not blurred after smoothing.

[0073] Finally, a comprehensive quality verification and evaluation was conducted: First, the signal-to-noise ratio (SNR) before and after processing was compared. By statistically analyzing the amplitude ratio of effective reflected signals to noise, the SNR of the processed data volume was ensured to be improved by more than 10% compared to the initial reconstructed data. Second, the residual amount of strong reflections was checked. The remaining amplitude in the original strong reflection area was statistically analyzed to ensure that the residual amount was <5%. Third, the continuity of weak reflection phase axes was checked through visual comparison. Seismic traces corresponding to 3 key survey lines and 5 wells were selected to confirm that the weak reflection phase axes were clearly identifiable and without obvious breaks. If the verification did not meet the standards, the gain parameters of amplitude compensation or the coherence threshold of noise suppression were adjusted retrospectively, and the data was reprocessed and verified again. If all indicators met the requirements, the final seismic data after strong axis removal was obtained, with high fidelity of weak reflection signals and continuous and smooth waveforms.

[0074] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for de-strong-axis processing in the curve wave domain based on the reflection coefficient property volume, characterized in that, The method includes: The original seismic data is subjected to spectral inversion processing to generate a volume of reflection coefficients; Based on the reflection coefficient volume, multi-scale and multi-directional layer-by-layer analysis is performed in the curve domain. Specifically, for each analysis window, the reflection coefficient curve is extracted, and the inflection point is determined by calculating the zero point of the second derivative of the curve to identify the abrupt change characteristics of the reflection coefficient. At the same time, a local discrete surface is constructed based on the reflection coefficient volume, and the local geometric anomaly region of the reflection coefficient volume is identified by fitting the surface patch and calculating the Gaussian curvature. By combining the determined inflection point locations and the identified local geometric anomaly regions, the strong reflection coefficient component is extracted; Based on the strong reflection coefficient components, a discrete data lattice is established to represent the spatial distribution of the strong reflection coefficient components. Based on the discrete data lattice, a wave field energy topology skeleton reflecting the spatial structure of the strong reflection energy is constructed. Several characteristic control points located inside, at the edge and outside of the strong reflection energy cluster are defined in the wave field energy topology skeleton. Based on the characteristic control points, an implicit energy isomorphic manifold representing the spatial distribution morphology of the strong reflection energy is fitted. The morphological properties and energy gradient distribution properties of the implicit energy isomorphic manifold are extracted, and the strong reflection stripping optimization coefficient is calculated accordingly. The strong reflection coefficient component is adaptively adjusted by using a strong reflection stripping optimization coefficient. The strong reflection coefficient component is then subtracted from the reflection coefficient volume to obtain the reflection coefficient volume after removing the strong axis. The seismic data is reconstructed based on the reflection coefficients after removing the strong axis, thus completing the strong axis removal process.

2. The method for removing strong axes in the curve wave domain based on the reflection coefficient property volume according to claim 1, characterized in that, The original seismic data undergoes spectral inversion processing to generate a volume of reflection coefficients, including: The raw seismic data is preprocessed to obtain the preprocessed seismic data volume; Seismic wavelets are extracted from the preprocessed seismic data volume. By combining the acquired geological data and well logging data of the work area, a low-frequency model for spectral inversion was constructed; Using the extracted seismic wavelet and the low-frequency model used for spectral inversion as constraints, the preprocessed seismic data volume is subjected to constrained sparse pulse inversion iterative calculations until convergence, thus obtaining the initial reflection coefficient volume. The initial reflection coefficient volume is normalized and output to generate the final reflection coefficient volume.

3. The method for removing strong axes in the curve wave domain based on the reflection coefficient property volume according to claim 2, characterized in that, Based on the final reflection coefficient volume, multi-scale, multi-directional layer-by-layer analysis is performed in the curve wave domain. Specifically, for each analysis window, the reflection coefficient curve is extracted, and the inflection point is determined by calculating the zero point of the second derivative of the curve, thus identifying abrupt changes in the reflection coefficient. Simultaneously, a local discrete surface is constructed based on the reflection coefficient volume. By fitting surface patches and calculating Gaussian curvature, local geometric anomaly regions of the reflection coefficient volume are identified, including: Perform a curvilinear transformation on the final reflection coefficient volume to obtain the corresponding curvilinear domain coefficient volume; In the corresponding curve domain coefficient volume, multi-scale, multi-directional analysis windows are set along the target layer; For each analysis window, the reflection coefficient sequence within the window is extracted to construct a one-dimensional reflection coefficient curve; the second derivative of the one-dimensional reflection coefficient curve is calculated, and the inflection point of the one-dimensional reflection coefficient curve is determined by locating the zero point of the second derivative. The inflection point is used to identify the abrupt change characteristics of the reflection coefficient. For each analysis window, a local discrete surface is constructed using the three-dimensional spatial data points within the window. The local discrete surface is fitted with surface patches, and the Gaussian curvature of each surface patch is calculated to form Gaussian curvature distribution data; Based on Gaussian curvature distribution data, identify local geometric anomaly regions in the reflectance coefficient volume.

4. The method for removing strong axes in the curve wave domain based on the reflection coefficient property volume according to claim 3, characterized in that, Based on the determined inflection point locations and identified local geometric anomaly regions, the strong reflection coefficient components are extracted, including: The inflection point locations determined within each analysis window are spatially fused with the identified local geometric anomaly regions to generate a comprehensive feature distribution map; Based on the comprehensive feature distribution map, a feature intensity threshold is set to filter out candidate regions that meet the strong reflectivity characteristics. In the curve wave domain coefficient volume, locate the curve wave domain coefficients corresponding to the candidate region; The curve domain coefficients of the location are extracted to form the strong reflection coefficient component.

5. The method for removing strong axes in the curve wave domain based on the reflection coefficient property volume according to claim 4, characterized in that, Based on the strong reflection coefficient components, a discrete data lattice is established to represent the spatial distribution of the strong reflection coefficient components; based on the discrete data lattice, a wave field energy topology framework reflecting the spatial structure of the strong reflection energy is constructed; and several characteristic control points located inside, at the edge and outside of the strong reflection energy cluster are defined in the wave field energy topology framework. Based on the aforementioned feature control points, an implicit energy isomorphic manifold characterizing the spatial distribution of strong reflective energy is formed, including: Three-dimensional spatial sampling is performed on the strong reflection coefficient components to establish a discrete data matrix of the spatial distribution of the strong reflection coefficient components; Based on the established discrete data lattice, a wave field energy topology skeleton reflecting the spatial structure of strong reflection energy is constructed by calculating the spatial correlation and energy coupling relationship between points. Based on energy intensity and topological position, several characteristic control points are defined in the wave field energy topology framework, located inside the strong reflection energy cluster, at the edge of the energy cluster, and outside the energy cluster. Based on the aforementioned feature control points and their spatial constraints, an implicit energy isomorphic manifold representing the spatial distribution of strong reflective energy is generated through implicit surface fitting calculations.

6. The method for removing strong axes in the curve wave domain based on the reflection coefficient property volume according to claim 5, characterized in that, The morphological properties and energy gradient distribution properties of the implicit energy isomorphic manifold are extracted, and the strong reflection stripping optimization coefficients are calculated accordingly, including: Geometric analysis is performed on implicit energy isovalued manifolds to extract morphological properties, including manifold surface area, volume, and curvature distribution characteristics. Energy field analysis is performed on implicit energy isomorphic manifolds to calculate the energy gradient distribution properties inside and on the surface; By combining the extracted morphological attributes with the calculated energy gradient distribution attributes, a multidimensional feature vector is constructed. Based on the constructed multidimensional feature vector, the strong reflection stripping optimization coefficient is calculated through the set optimization mapping function.

7. The method for removing strong axes in the curve wave domain based on the reflection coefficient property volume according to claim 6, characterized in that, An optimization coefficient for strong reflection is used to adaptively adjust the components of the strong reflection coefficient. By subtracting the adjusted strong reflection coefficient components from the reflection coefficient volume, the reflection coefficient volume after removing the strong axis is obtained, including: The strong reflection coefficient components are weighted and adjusted using a strong reflection stripping optimization coefficient to obtain the weighted strong reflection coefficient components. The weighted strong reflection coefficient component is subtracted from the corresponding position in the reflection coefficient volume to generate the intermediate reflection coefficient volume; The intermediate reflectance volume is subjected to energy equalization to obtain the reflectance volume after the strong axis is removed.

8. The method for removing strong axes in the curve wave domain based on the reflection coefficient property volume according to claim 7, characterized in that, Based on the reflection coefficient volume reconstructed after removing the strong axis, the strong axis removal process is completed, including: Extract the seismic wavelet used in the process of generating the reflection coefficient volume; The extracted seismic wavelet is convolved with the reflection coefficient volume after removing the strong axis to generate the initial reconstructed seismic data volume. Amplitude compensation and residual noise suppression are performed on the initial reconstructed seismic data volume to obtain the final seismic data after strong axis removal, which represents the completion of strong axis removal.

Images