Interference suppression method of VSP multi-component polarization and frequency-wave number domain in metal mine area

Abstract

The application belongs to the technical field of geophysical exploration, and relates to a VSP multi-component polarization and frequency-wavenumber domain interference suppression method, time-depth domain data is obtained through multi-component VSP data preprocessing, multi-component polarization parameters are extracted by using a dynamically adjusted moving time window, polarization differences of effective waves and interference waves are captured, a classification result is corrected based on polarization parameter clustering and combined with geological prior knowledge of a high impedance value interface, a time-depth domain interference mask is generated, data and the interference mask are synchronously mapped to a frequency-wavenumber domain through an overlapping window technology, a joint domain optimization target is constructed through polarization parameters and a frequency domain interference mask, filter parameters are dynamically adjusted, and while the interference suppression term suppresses frequency domain interference, the polarization fidelity term is used to constrain the polarization characteristics of effective waves from being damaged, synchronous achievement of interference suppression and polarization fidelity is realized, and finally, data signal-to-noise ratio and imaging quality are improved through inverse transformation reconstruction.

Description

Technical Field

[0001] This application belongs to the field of geophysical exploration technology, and more specifically, relates to a method for suppressing interference in the VSP multi-component polarization and frequency-wavenumber domain in metal mining areas. Background Technology

[0002] Vertical seismic profiling (VSP) technology is a key tool for exploring deep ore bodies in metal mining areas. By deploying three-component geophones in the well to collect seismic data, it can directly obtain the longitudinal and lateral geological information of the subsurface medium, providing an important basis for ore body location and reserve assessment. However, metal mining areas are generally characterized by complex geological structures and the development of high-impedance interfaces (such as the contact surface between the ore body and the surrounding rock, and fracture zones). These high-impedance interfaces can cause strong scattering interference. At the same time, seismic waves are easily affected by dispersion effects during propagation, resulting in the confusion of polarization characteristics and overlapping frequency domain distributions between effective waves and interference waves in the VSP multi-component data.

[0003] Existing VSP data interference suppression methods are mainly divided into two categories: one is frequency-wavenumber domain-based filtering methods, which suppress interference within a specific wavenumber range by mapping the data to the frequency domain. However, this type of method does not consider the polarization difference between the effective wave and the interfering wave, and is prone to damaging the polarization characteristics of the effective wave while suppressing the interference, resulting in distortion of the effective wave propagation information. The other type is interference suppression methods based on multi-component polarization analysis, which distinguish between the effective wave and the interfering wave by extracting polarization parameters. However, this type of method has poor adaptability to polarization distortion caused by strong scattering interference and is difficult to handle the frequency domain diffusion problem of interference caused by dispersion effects.

[0004] In practical applications in metal mining areas, both of the above methods have significant limitations: due to the easy overlap and frequency domain distribution of strong scattering interference from high-impedance interfaces and the polarization characteristics of effective waves, relying solely on frequency domain filtering or polarization analysis methods cannot achieve accurate suppression of interference and simultaneous preservation of the polarization characteristics of effective waves. This ultimately results in low signal-to-noise ratio and poor imaging quality of the processed data, making it difficult to meet the needs of fine exploration of deep ore bodies in metal mining areas. Therefore, how to effectively suppress strong scattering interference and accurately preserve the polarization characteristics of effective waves in VSP multi-component data processing in metal mining areas has become a core technical problem that urgently needs to be solved. Summary of the Invention

[0005] This invention provides a method for suppressing interference in the polarization and frequency-wavenumber domains of VSP multi-components in metal mining areas. It aims to solve the technical problem that existing technologies in the processing of VSP multi-component data in metal mining areas are unable to simultaneously achieve effective suppression of strong scattering interference and accurate preservation of effective wave polarization characteristics.

[0006] A method for suppressing VSP interference in metal mining areas using multiple polarization and frequency-wavenumber domains includes the following steps:

[0007] Step 1: Deploy a three-component acquisition system in the metal mining area to acquire raw multi-component VSP seismic data. After mean removal, time difference correction and bandpass filtering, preprocessed multi-component VSP data is obtained.

[0008] Step 2: Perform covariance matrix analysis on the preprocessed multi-component VSP data using a moving time window, and extract multi-component polarization parameters through eigenvalue decomposition. The length of the moving time window is dynamically adjusted according to the local signal-to-noise ratio of the preprocessed multi-component VSP data.

[0009] Step 3: Construct feature vectors based on the multi-component polarization parameters, perform initial classification of effective waves and interference waves in the preprocessed multi-component VSP data using a clustering algorithm, and correct the classification results by combining geological prior knowledge of interfaces in metal mining areas with higher than preset impedance values, and generate a time-depth domain interference mask.

[0010] Step 4: Using the overlapping window technique, perform two-dimensional Fourier transform on the preprocessed multi-component VSP data and the time-depth domain interference mask respectively, mapping from the time-depth domain to the frequency-wavenumber domain to obtain the frequency domain multi-component VSP data and the frequency domain interference mask.

[0011] Step 5: Based on the multi-component polarization parameters and the frequency domain interference mask, dynamically adjust the cutoff frequency and attenuation coefficient of the filter, construct a joint domain optimization objective that includes interference suppression and polarization fidelity terms, and filter the frequency domain multi-component VSP data to obtain the filtered frequency domain multi-component VSP data.

[0012] Step 6: Restore the filtered frequency domain multi-component VSP data to the time-depth domain using two-dimensional inverse Fourier transform, eliminate window boundary distortion using the overlap-addition method, calibrate the data amplitude and baseline, and obtain the reconstructed multi-component VSP data.

[0013] This invention obtains standardized time-depth domain data through multi-component VSP data preprocessing, laying the foundation for subsequent analysis. Then, a dynamically adjusted moving time window is used to extract multi-component polarization parameters, accurately capturing the polarization differences between effective and interfering waves and solving the problem of polarization feature confusion. Subsequently, based on polarization parameter clustering and classification, and combined with geological prior knowledge of interfaces above the preset impedance value, the classification results are corrected to generate a time-depth domain interference mask, achieving precise location of areas with strong scattering interference. Overlapping window technology is used to synchronously map the data and interference mask to the frequency-wavenumber domain, avoiding interference identification ambiguity caused by frequency domain diffusion. A joint domain optimization objective is constructed using polarization parameters and the frequency domain interference mask, and filter parameters are dynamically adjusted. While the interference suppression term specifically suppresses frequency domain interference, the polarization fidelity term constrains the effective wave polarization characteristics to prevent damage, achieving simultaneous interference suppression and polarization fidelity preservation. Finally, inverse transformation reconstruction improves the data signal-to-noise ratio and imaging quality.

[0014] Preferably, the dynamic adjustment of the moving time window includes the following steps:

[0015] Set minimum and maximum window length thresholds, perform moving average smoothing on the local signal-to-noise ratio of the preprocessed multi-component VSP data, calculate the real-time window length using an exponential decay function, and increase the window length according to the power law of the exponential decay function based on the decrease in signal-to-noise ratio, and the real-time window length is always between the set minimum and maximum window length thresholds.

[0016] Preferably, the covariance matrix analysis uses Gaussian weighting factor optimization. The standard deviation of the Gaussian weighting factor is proportional to the window length. The weight value is calculated by the distance between each data point in the window and the center of the window. The greater the distance, the smaller the weight value. That is, the weight is determined by combining the relationship between the distance from the data point to the center and the standard deviation through the natural exponential function.

[0017] Preferably, the clustering algorithm uses K-means unsupervised clustering, and the standardized feature vector is input into the clustering algorithm for classification; wherein the feature vector is composed of polarization angle, ellipticity and polarization intensity in multi-component polarization parameters;

[0018] The standardization is defined as subtracting the mean of each parameter and then dividing it by its own standard deviation.

[0019] The initial cluster centers of the clustering algorithm are determined by randomly selecting multiple sample points with different polarization characteristics. Then, the algorithm iterates through a set number of iterations. When the movement distance of the cluster centers reaches a preset threshold, the optimal number of clusters is determined by the contour coefficient. After clustering, the cluster labels are converted into an initial interference probability matrix, which serves as the basis for correcting geological prior knowledge.

[0020] Preferably, the geological prior knowledge-corrected classification result includes the following steps:

[0021] Obtain the time-depth range of interfaces with known impedance values ​​higher than the preset impedance value, construct a geological prior weight matrix consistent with the dimensions of the preprocessed multi-component VSP data, in the geological prior weight matrix, the weight coefficient of the interface with the preset impedance value within the influence range is greater than 1, and the closer to the interface, the larger the weight coefficient is, the weight coefficient of the interface that is not higher than the preset impedance value is set to 1, and when the distance exceeds the set range, it is completely restored to 1.

[0022] The obtained initial interference probability matrix is ​​multiplied point by point with the geological prior weight matrix to obtain the weighted interference probability;

[0023] The threshold is then dynamically adjusted based on the local signal-to-noise ratio of the preprocessed multi-component VSP data; the lower the signal-to-noise ratio, the higher the threshold.

[0024] The weighted interference probability is compared with the threshold. When the weighted interference probability is greater than or equal to the threshold, it is marked as an interference region; otherwise, it is marked as a non-interference region, thereby generating a time-depth domain interference mask.

[0025] Preferably, the overlap rate of the overlapping window technology is set according to a preset ratio, and the preprocessed multi-component VSP data is symmetrically extended before the two-dimensional Fourier transform, with the extension length having a fixed proportional relationship with the window length; after the two-dimensional Fourier transform, the window data is spliced ​​by a weighted average method, and the weights of the overlapping areas of adjacent windows are executed according to a linear allocation rule, with the weight of the previous window gradually decreasing and the weight of the next window gradually increasing, so as to achieve smooth splicing of the window data.

[0026] Preferably, the joint domain optimization objective is to balance the interference suppression term and the polarization fidelity term using weighting coefficients, which are determined by cross-validation.

[0027] The interference suppression term is calculated based on the amplitude weight of the frequency domain interference mask; the larger the amplitude, the higher the corresponding weight.

[0028] The polarization fidelity term is calculated based on the differences in multi-component polarization parameters in the effective wave region; the differences in multi-component polarization parameters are obtained by comparing the polarization angle, ellipticity, and polarization intensity in the multi-component polarization parameters with the polarization parameters obtained by inversion from the filtered data.

[0029] Preferably, after the filtering process in step 5, frequency domain phase calibration is also included:

[0030] The phase shift of the filtered frequency domain data is calculated based on the effective wave polarization angle, and then the phase shift is corrected by a linear phase compensation method. The compensation coefficient mask in the compensation method is proportional to the polarization angle difference. The filtered frequency domain data is calibrated point by point by a phase correction factor.

[0031] Preferably, the amplitude calibration includes the following steps:

[0032] The ratio of the maximum amplitude of the preprocessed multi-component VSP data to the maximum amplitude of the reconstructed multi-component VSP data is calculated. Each sampling point of the reconstructed multi-component VSP data is multiplied by the ratio to ensure that the amplitude range of the reconstructed multi-component VSP data is consistent with that of the preprocessed multi-component VSP data. Furthermore, the window parameters for the overlap-addition method elimination are completely consistent with the window parameters for the frequency-wavenumber transformation.

[0033] The beneficial effects of the invention include:

[0034] This invention obtains standardized time-depth domain data through multi-component VSP data preprocessing, laying the foundation for subsequent analysis. Then, a dynamically adjusted moving time window is used to extract multi-component polarization parameters, accurately capturing the polarization differences between effective and interfering waves and solving the problem of polarization feature confusion. Subsequently, based on polarization parameter clustering and classification, and combined with geological prior knowledge of interfaces above the preset impedance value, the classification results are corrected to generate a time-depth domain interference mask, achieving precise location of areas with strong scattering interference. Overlapping window technology is used to synchronously map the data and interference mask to the frequency-wavenumber domain, avoiding interference identification ambiguity caused by frequency domain diffusion. A joint domain optimization objective is constructed using polarization parameters and the frequency domain interference mask, and filter parameters are dynamically adjusted. While the interference suppression term specifically suppresses frequency domain interference, the polarization fidelity term constrains the effective wave polarization characteristics to prevent damage, achieving simultaneous interference suppression and polarization fidelity preservation. Finally, inverse transformation reconstruction improves the data signal-to-noise ratio and imaging quality. Attached Figure Description

[0035] To more clearly illustrate the technical solutions in the embodiments of this application, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0036] Figure 1 The overall flowchart provided for embodiments of the present invention.

[0037] Figure 2 The flowchart for step 3 provided in the embodiment of the present invention is shown.

[0038] Figure 3The flowchart for step 5 provided in this embodiment of the invention is shown.

[0039] Figure 4 This is a schematic diagram showing the variation of filter parameters with effective wave confidence level provided in an embodiment of the present invention.

[0040] Figure 5 This is a schematic diagram showing the comparison of frequency domain data before and after filtering, provided in an embodiment of the present invention.

[0041] Figure 6 This is a schematic diagram showing the comparison of the fidelity of effective wave polarization parameters before and after filtering, provided in an embodiment of the present invention. Detailed Implementation

[0042] To make the technical problems, technical solutions, and beneficial effects to be solved by this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and are not intended to limit the scope of this application.

[0043] See Figure 1 As shown, this embodiment provides a method for suppressing interference in the VSP multi-component polarization and frequency-wavenumber domains in metal mining areas, including the following steps:

[0044] Step 1: Acquire multi-component VSP data and perform preprocessing, specifically including the following steps:

[0045] Multi-component VSP data input and preliminary verification: A three-component VSP acquisition system was deployed in a metal mining area. The system includes east-west... , north-south ,vertical Three-component detectors, with a deployment depth range of 0-1000m and a time sampling interval. Total recording time Depth sampling interval Total depth number of channels Raw multi-component VSP seismic data were acquired using the three-component VSP acquisition system, and the data format was SEG-Y standard format.

[0046] Perform preliminary quality checks on the input data, calculate the signal-to-noise ratio (SNR) to assess data usability, and the SNR... The calculation formula uses a piecewise calculation method: ;

[0047] In the formula: The effective signal power is calculated by the square mean of the signal amplitude within the effective wave time window (e.g., 0.5-1.5s) determined by geological priors. To obtain the noise power, the mean square of the amplitude is calculated from the initial no-signal segment (0-0.1s) at the beginning of the recording.

[0048] The signal-to-noise ratio (SNR) threshold is set to 10dB. If the SNR of a certain component data is lower than 10dB, the data will be reacquired or marked as a low-quality data segment.

[0049] Data preprocessing and filtering: First, the calibrated multi-component data is subjected to mean removal to eliminate instrument offset and environmental baseline drift. Each component is processed independently. ;

[0050] In the formula: For the original time-domain data of the i-th component, The time mean of the i-th component. , This represents the total number of sampling points (N=1000 in this embodiment). The time for the kth sampling point; This represents the time-domain data of the i-th component after removing the mean.

[0051] Subsequently, bandpass filtering was performed to remove high-frequency instrument vibration noise (>100Hz) and low-frequency temperature drift (<10Hz). Specifically, a Butterworth linear phase filter was used (phase linearity avoids distortion of the effective wave signal and adapts to the phase propagation characteristics of seismic waves). The passband frequency range was set to 10-100Hz, and the filter order was 4th. The filtered data is denoted as... .

[0052] Step 2: Covariance matrix analysis is performed on the preprocessed multi-component VSP data using a moving-time window. Multi-component polarization parameters are extracted through eigenvalue decomposition. The length of the moving-time window is dynamically adjusted based on the local signal-to-noise ratio of the preprocessed multi-component VSP data. In this embodiment, due to the high scattering characteristics of metal mining areas, seismic wave polarization characteristics undergo drastic spatiotemporal changes. Calculating polarization parameters with a fixed window can lead to parameter ambiguity. Therefore, an adaptive window sliding calculation strategy is introduced, dynamically adjusting the window size based on the local signal-to-noise ratio. Combined with weighted covariance matrix optimization, this improves the accuracy of polarization parameter calculation, providing reliable feature indicators for subsequent interference classification. The specific steps are as follows:

[0053] Adaptive time window setting: First, the signal-to-noise ratio obtained in step 1 is smoothed by moving average to eliminate the influence of local fluctuations. The smoothing formula is: ;

[0054] In the formula: The number of sampling points within the smoothing window (100). The moving average window length (e.g., 0.2s). The time sampling interval is 0.002s. Indicates the current time; Let be the smoothed signal-to-noise ratio at time t; Indicates the current time;

[0055] The window length is dynamically calculated based on the smoothed signal-to-noise ratio, using an exponential decay function. ;

[0056] In the formula: Let be the real-time window length at time t; A minimum window length (e.g., 0.05 s) is used to ensure the capture of rapidly changing polarization characteristics; The maximum window length (e.g., 0.2s) is used to ensure that there are enough data points in the low signal-to-noise ratio region to calculate stable parameters; The attenuation coefficient (e.g., 0.5) is calibrated using actual data from the mining area to balance the window size and parameter accuracy; in this embodiment, the high signal-to-noise ratio region... The window length is approximately 0.06 seconds, in the low signal-to-noise ratio region. The window length is approximately 0.15 seconds, all within... and Between; where the window movement step size This ensures that adjacent windows have sufficient overlap to avoid abrupt changes in polarization parameters.

[0057] Weighted covariance matrix calculation: To suppress the interference of window edge data on the covariance matrix calculation (window edge data is easily affected by adjacent wave fields, leading to parameter deviation), a Gaussian weighting factor is introduced to optimize the covariance matrix. Gaussian weighting factor calculation: ;

[0058] In the formula: , which is the standard deviation and is in a fixed proportion to the window length, to ensure that the weights decay reasonably within the window; The current time is the center time of the window. The time of the kth data point within the window has a weight range of [0,1], with the center data of the window having the largest weight and the edge data having the smallest weight. This represents the Gaussian weight value for the k-th data point within the window.

[0059] Construct a weighted covariance matrix based on weighting factors. Formula for calculating matrix elements: ;

[0060] In the formula: , is the component index; Let be the element in the i-th row and j-th column of the weighted covariance matrix; , where is the number of sampling points within the window; For the i-th component after preprocessing Data at any given time; For the j-th component after preprocessing Data at any given time; A real symmetric matrix with dimensions equal to the square of the amplitude. In this embodiment, through weighted processing, the weight of the window edge data is reduced by more than 60% compared with the center data, which effectively reduces edge interference.

[0061] Instantaneous polarization parameter extraction: For the weighted covariance matrix Perform eigenvalue decomposition to obtain eigenvalues. and corresponding feature vectors The following polarization parameters are extracted based on eigenvalues ​​and eigenvectors:

[0062] polarization angle : Describes the angle between the wave propagation direction and the principal axis of the polarization ellipsoid, based on the largest eigenvalue. Corresponding feature vector calculate: ,in The largest eigenvalue Corresponding feature vector Components in the x-component (east-west); For feature vectors Components in the y-component;

[0063] Ellipticity E: Quantifies the flattening of the polarization ellipse, distinguishing between linear waves (such as P-waves) and elliptically polarized waves (such as S-waves). Calculation formula: The range is [0,1], and the closer E is to 0, the closer the polarization is to linear; This represents the smallest eigenvalue obtained after eigenvalue decomposition of the weighted covariance matrix.

[0064] Polarization intensity I: reflects the stability of polarization, calculated using the following formula: , dimensionless, range [0,1], the closer I is to 1, the more stable the polarization.

[0065] Based on the above feature extraction, multi-component polarization parameters, i.e., a two-dimensional polarization parameter matrix, are obtained. This serves as the feature vector for subsequent interference wave classification.

[0066] Step 3: See Figure 2As shown, due to the fundamental differences in polarization characteristics between the effective wave and the interfering wave (the effective wave has stable polarization and low ellipticity, while the interfering wave has chaotic polarization and abnormal ellipticity), clustering algorithms can initially separate the two. However, high-impedance interfaces in metal mining areas are prone to strong scattering interference, and purely data-driven clustering is prone to misjudgment. Therefore, the classification results are corrected by combining prior geological knowledge to generate an accurate time-depth domain interference mask, providing a clear identifier of the interference region for subsequent frequency domain filtering. The specific steps are as follows:

[0067] Polarization attribute normalization and eigenvector construction: First, the two-dimensional polarization parameter matrix output in step 2 is... Reshape into a one-dimensional vector (number of samples) , This represents the number of time sampling points; The depth (channel number) is used to eliminate dimensional differences; then, each parameter is Z-score normalized to eliminate differences in dimensions and numerical range, resulting in the normalized polarization angle. Ellipticity polarization intensity ;

[0068] The three standardized parameters of each sample are combined into a three-dimensional feature vector. This characterizes the polarization state of seismic waves at the corresponding time-depth point.

[0069] Initial classification of interference waves based on K-means: The K-means unsupervised clustering algorithm is used to classify the feature vectors. First, the optimal number of clusters K is determined by the silhouette coefficient, where the silhouette coefficient... The calculation formula is: ,in Let i be the average distance between sample i and other samples in the same cluster. Let i be the average distance between sample i and the nearest heterogeneous sample. The range is [-1, 1]. The closer the value is to 1, the better the clustering effect. In this embodiment, the average profile coefficient is 0.72 when K is 2 and 0.65 when K is 3. Therefore, the optimal number of clusters is determined to be K=2 (effective clusters and interfering clusters).

[0070] The initial cluster centers are determined by randomly selecting 10 sample points with different polarization characteristics. The maximum number of iterations is set to 50, and the iteration stops when the distance the cluster centers have moved is less than 10. The goal of the K-means algorithm is to minimize the within-cluster sum of squares (WCSS). ;

[0071] In the formula: For the k-th cluster; It is the cluster centroid; Let i be the i-th eigenvector; This represents the i-th eigenvector;

[0072] After clustering, the statistical characteristics of polarization parameters of each cluster were analyzed: the mean polarization intensity of the effective clusters was 0.82, the mean ellipticity was 0.15 (close to linear), and the standard deviation of the polarization angle was 0.12 (stable); the mean polarization intensity of the interfering clusters was 0.35, the mean ellipticity was 0.68 (close to circular polarization), and the standard deviation of the polarization angle was 0.85 (random). Based on this, an initial interference mask was generated. , Indicates the area of ​​interference. Indicates the effective wave region.

[0073] Based on prior geological knowledge, the classification results were revised: Known high-impedance interfaces (such as the interface between the ore body and the surrounding rock) in metal mining areas are prone to strong scattering interference, necessitating an increase in the interference identification weight for this region; firstly, the depth range of known high-impedance interfaces was determined. Corresponding time range Constructing a geological prior weight matrix (Consistent with data dimensions), if and ,but ,otherwise ;in , which is an empirical weighting coefficient, to enhance the probability of identifying interference near high impedance interfaces; The depth (the depth at which the geophone is deployed in the well, in meters, and the range of values ​​is the depth range of the geophone deployment, such as 0–1000m in the example). Indicates time (the recording time of the seismic data, with a value range of the total recording time interval, such as 0–2s in the example); The vertical depth range affected by the interface; This represents the time-related influence range of the interface; in the weight matrix, the closer to the interface, the closer the weight coefficient is to 1.5, and the farther away it is, the smaller the weight coefficient is to 1.5. and It is then restored to the baseline weight of 1.

[0074] The initial interference probability matrix (transformed from the initial cluster labels, with higher probabilities for sample points closer to the centroid of interference clusters) is multiplied point-by-point by the weight matrix to obtain the weighted interference probability. Subsequently, the threshold is dynamically adjusted based on the local signal-to-noise ratio. This allows for the use of a higher threshold in low signal-to-noise ratio regions, i.e.: ,in The base threshold is 0.5. The adjustment factor is 0.3.

[0075] The weighted interference probabilities are thresholded to generate the final interference mask. :when hour ,otherwise .

[0076] Step 4: In the time-depth domain, the dispersion characteristics of the effective wave and the interfering wave are not significantly different. By mapping the data to the frequency-wavenumber domain using a two-dimensional Fourier transform, separation can be achieved by utilizing the difference in their wavenumber-frequency distributions. The overlapping window technique and symmetrical boundary extension are employed to suppress spectral leakage and ensure the spatial correspondence between the frequency domain data and the interference mask. The specific steps are as follows:

[0077] Window parameter design: preprocessed multi-component data and interference mask All are 1000×500 two-dimensional matrices; considering the dominant frequency (10-100Hz) of seismic waves in the mining area and the characteristics of geological layer thickness, the window parameters are designed as follows: time window length. Depth window length Number of sampling points within the time window Number of sampling points within the depth window Overlap rate Adjacent windows overlap by 50 sampling points in the time dimension and 25 sampling points in the depth dimension to ensure smooth window stitching.

[0078] Overlapping window segmentation and data extension: Overlapping window segmentation is performed on multi-component data and interference masks according to design parameters, specifically, the sliding step size in the time dimension. A total of A time window; depth dimension sliding step size A total of A depth window.

[0079] Symmetrical boundary extension is performed on the sub-data within each window, with an extension length of 1 / 4 of the window length (25 sampling points for the time window and 12 sampling points for the depth window) to avoid boundary distortion after transformation. After extension, the window dimension becomes 150×75.

[0080] Two-dimensional Fourier transform: number of sub-windows for each extended version Perform a two-dimensional Fourier transform to map the time-depth domain data to the frequency-wavenumber domain. The transformation formula is as follows: In the formula: For the frequency domain data of the i-th component under the m-th time window and the n-th depth window ; The imaginary unit; For frequency; Wave number; For local time indexing within the window; The time sampling interval; For local depth index within the window; This represents the depth sampling interval.

[0081] The frequency domain data of all windows are stitched together using an overlap-addition method: for each frequency-wavenumber point... Extract all window transformation results covering the point, and apply a linear weighted average to the overlapping areas (the weight of the previous window linearly decreases from 1 to 0.5, and the weight of the next window linearly increases from 0.5 to 1), finally obtaining complete frequency domain multi-component data. , , .

[0082] Frequency domain interference mask mapping: Using parameters completely consistent with the data transformation (window size, overlap rate, continuation method), the time-depth domain interference mask is mapped. Perform a two-dimensional Fourier transform to obtain the frequency domain interference mask. Its amplitude This indicates the percentage of interference energy at the corresponding frequency-wavenumber point, providing an identifier for the interference region in subsequent filtering.

[0083] Step 5: Because traditional frequency domain filtering tends to damage the effective wave polarization characteristics while suppressing interference, see [link to relevant documentation]. Figure 3 As shown, this invention introduces the polarization parameters extracted in step 2 as dynamic constraints to construct a joint polarization-frequency-wavenumber domain optimization objective. By balancing the interference suppression term and the polarization fidelity term, the filter parameters are dynamically adjusted to achieve simultaneous precise interference suppression and effective wave polarization fidelity preservation; specifically as follows:

[0084] Filter constraint coefficient calculation: By fusing polarization parameters and interference mask information, a point-by-point effective wave confidence coefficient in the frequency-wavenumber domain is constructed. The formula for quantifying the reliability of effective waves is: ;

[0085] In the formula: Weighting coefficients (e.g.) ); This represents the polarization intensity mapped to the frequency domain. This is the ellipticity mapped to the frequency domain; Indicates a frequency domain interference mask; The range is [0,1]. The closer the value is to 1, the higher the confidence level of the effective wave.

[0086] right Global normalization is performed to obtain the normalized effective wave confidence coefficient. .

[0087] Construction of the adaptive filter response function: A 4th-order Butterworth filter is used, based on... Dynamically adjust the cutoff frequency and attenuation coefficient:

[0088] Dynamic cutoff frequency: ,in The base cutoff frequency is approximately 1.3 times the effective wave frequency of 30Hz. This indicates the adjustment factor (e.g., a value of 0.5); the region with high confidence in the effective wave. The cutoff frequency is extended to 60Hz, covering the entire effective wave band; interference area The cutoff frequency is reduced to 40Hz, compressing the interference band.

[0089] Dynamic attenuation coefficient: ,in , representing the maximum attenuation coefficient; see [link to relevant documentation] Figure 4 As shown, the interference area ( The attenuation coefficient is 40dB, strongly suppressing interference; the effective wave region ( The attenuation coefficient is close to 0dB, so it hardly attenuates the effective wave.

[0090] Constructing the filter response function : ;

[0091] In the formula: express The filter response function at position, Indicates frequency, Indicates wave number; Indicates the dynamic cutoff frequency; This represents the betweenness number of the filter, with a value of 4.

[0092] Optimization of the joint domain loss function of polarization-frequency-wavenumber: Constructing a joint domain loss function L to balance interference suppression and polarization fidelity: ;

[0093] In the formula: Indicates the interference suppression weight; Represents the polarization fidelity weight, satisfying This was determined through cross-validation; For the interference energy suppression term, This is the effective wave polarization fidelity term; in this embodiment, the focus is on interference suppression. The value is 0.6.

[0094] Interference energy suppression term: ;

[0095] In the formula: , These represent the number of sampling points for frequency and wavenumber, respectively. Indicates a frequency domain interference mask; Represents the original frequency domain data of the i-th component; the interference energy suppression term indicates the residual interference energy after quantization filtering, and the smaller the value, the better the interference suppression effect.

[0096] Effective wave polarization fidelity term: ;

[0097] In the formula: , , Polarization parameters retrieved from filtered data , , The values ​​represent the original polarization angle, polarization intensity, and ellipticity; the effective wave polarization fidelity term quantifies the degree of distortion of the effective wave polarization characteristics, with smaller values ​​indicating better fidelity.

[0098] Optimization using gradient descent method Iteration formula: ,in , where is the learning rate; Let be the filter response function for the current iteration; This is the updated filter response function; This represents the partial derivative of the total loss function L with respect to the current filter response function; the iteration termination condition is... The iteration count may reach 50 times; after optimization, the loss function L is reduced by 78% compared with the initial value, achieving a balance between interference suppression and polarization fidelity.

[0099] Filter Applications and Phase Calibration: Optimizing the filter response function Applied to frequency domain multi-component data: This represents the filtered frequency domain data;

[0100] Since filtering may introduce phase shift, phase calibration is performed based on the effective wave polarization angle. ;

[0101] in , is the phase correction amount; This is the phase correction factor; after calibration, the effective wave phase offset deviation is controlled within 0.02 rad to ensure that the effective wave propagation characteristics are not affected; the final output is filtered frequency domain data. See also Figure 5As shown, the left figure is the amplitude distribution of the frequency domain data before filtering, and the right figure is the amplitude distribution after joint domain optimization filtering. As can be seen from the figure, there is obvious strong interference energy in the 20-30Hz and 110-150Hz frequency bands before filtering (the maximum amplitude in the 20-30Hz frequency band before filtering reaches 1.8V·s). This type of interference mainly comes from the strong scattering effect and dispersion effect of the high impedance interface in the metal mining area. After the joint domain filtering process in step 5 of this application, the energy suppression rate of the above-mentioned interference frequency bands reaches 84%, and the amplitude of the effective wave band in the 10-100Hz frequency band does not have obvious attenuation, maintaining the energy integrity of the effective wave, realizing targeted suppression of the interference frequency band, and avoiding damage to the effective wave band.

[0102] See also Figure 6 As shown in the figure, the variation trends of three core parameters—polarization angle (rad), ellipticity, and polarization intensity—are illustrated. The figure reveals that before filtering, the polarization parameters in the interference region (e.g., 0.8-0.9s, corresponding to the time range of 400-420m at the high-impedance interface in the embodiment) are chaotic, with an average polarization intensity of only 0.35, an average ellipticity of 0.68 (close to circular polarization), and a standard deviation of 0.85 for the polarization angle. After joint domain filtering, the polarization parameters in the effective wave region remain stable, with an average polarization intensity of 0.82 (consistent with the statistical values ​​of effective wave cluster polarization intensity in the embodiment), an average ellipticity of 0.15 (close to linear polarization, consistent with the polarization characteristics of effective P-waves and S-waves), and a standard deviation of 0.12 for the polarization angle, highly consistent with the extracted original effective wave polarization parameter characteristics. Furthermore, the time-domain variation trend of the polarization parameters after filtering is continuous without abrupt changes, achieving accurate preservation of the effective wave polarization characteristics during strong interference suppression.

[0103] Step 6: The filtered frequency domain data needs to be converted back to the time-depth domain before it can be used for subsequent imaging. Using the same inverse transform parameters as the forward transform can avoid data distortion. Window boundary distortion is eliminated by the overlap-addition method, and amplitude calibration is used to ensure that the reconstructed data has the same dimensions as the original data, thus improving data reliability. Specifically:

[0104] Frequency domain data preprocessing and parameter verification: First, verify the dimensionality integrity of the filtered frequency domain data to ensure consistency with the dimensionality of the frequency domain data output in step 4; then perform amplitude preprocessing on the frequency domain data to correct the amplitude scaling effect of the Fourier transform. ,in This represents the total number of sampling points for the original data. This represents the total number of sampling points for the frequency domain data, ensuring that the amplitude of the frequency domain data matches that of the time domain data.

[0105] Single-window two-dimensional inverse Fourier transform: Following the window segmentation logic of step 4, the preprocessed frequency domain data is divided into 19×19 sub-frequency domain matrices. A two-dimensional inverse Fourier transform is performed on each sub-matrix, as shown in the following formula: ;

[0106] in The number of sampling points in the extended window is used to obtain the local time-depth data within the window after inverse transformation. ; This indicates the number of sampling points in the extended frequency window; This indicates the number of sampling points in the extended wavenumber window; This represents the frequency domain data of the i-th component within the m-th time window and the n-th depth window; For local time indexing within the window; This is a local depth index within the window.

[0107] Overlap-addition data stitching: Removes the boundary extensions of each local data set, restoring the original 100×50 window size sub-data. A linear weighted averaging strategy consistent with the forward transform is used to stitch adjacent windows together: the weights of overlapping time and depth regions are linearly transitioned to eliminate differences in window boundaries; the stitched data yields preliminary reconstructed data. .

[0108] Post-processing of reconstructed data: Mean removal is performed on the initially reconstructed data to eliminate baseline drift. ,in for The global mean;

[0109] Then perform amplitude normalization calibration to ensure that the amplitude range of the reconstructed data is consistent with that of the data after preprocessing in step 1: ;in This represents the maximum amplitude of the i-th component data after preprocessing in step 1. The maximum amplitude of the i-th component data after initial reconstruction; the final output is the reconstructed time-depth domain multi-component data. .

[0110] The above are merely preferred embodiments of this application and are not intended to limit this application. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this application should be included within the protection scope of this application.

Claims

1. A method for suppressing interference in the multi-component polarization and frequency-wavenumber domains of VSP in metal mining areas, characterized in that, Includes the following steps: Step 1: Deploy a three-component acquisition system in the metal mining area to acquire raw multi-component VSP seismic data. After mean removal, time difference correction and bandpass filtering, preprocessed multi-component VSP data is obtained. Step 2: Perform covariance matrix analysis on the preprocessed multi-component VSP data using a moving time window, and extract multi-component polarization parameters through eigenvalue decomposition. The length of the moving time window is dynamically adjusted according to the local signal-to-noise ratio of the preprocessed multi-component VSP data. Step 3: Construct feature vectors based on the multi-component polarization parameters, perform initial classification of effective waves and interference waves in the preprocessed multi-component VSP data using a clustering algorithm, and correct the classification results by combining geological prior knowledge of interfaces in metal mining areas with higher than preset impedance values, and generate a time-depth domain interference mask. The clustering algorithm uses K-means unsupervised clustering, and the standardized feature vectors are input into the clustering algorithm for classification; The eigenvectors are composed of polarization angle, ellipticity and polarization intensity from the multi-component polarization parameters. The standardization is defined as subtracting the mean of each parameter and then dividing it by its own standard deviation. The initial cluster centers of the clustering algorithm are determined by randomly selecting multiple sample points with different polarization characteristics, and then it is iterated through a set fixed number of iterations. When the moving distance of the cluster centers reaches a preset threshold, the number of clusters is determined by the silhouette coefficient. After clustering, the cluster labels are converted into an initial disturbance probability matrix, which serves as the basis for correcting geological prior knowledge. Step 4: Using the overlapping window technique, perform two-dimensional Fourier transform on the preprocessed multi-component VSP data and the time-depth domain interference mask respectively, mapping from the time-depth domain to the frequency-wavenumber domain to obtain the frequency domain multi-component VSP data and the frequency domain interference mask. Step 5: Based on the multi-component polarization parameters and the frequency domain interference mask, dynamically adjust the cutoff frequency and attenuation coefficient of the filter, construct a joint domain optimization objective that includes interference suppression and polarization fidelity terms, and filter the frequency domain multi-component VSP data to obtain the filtered frequency domain multi-component VSP data. Step 6: Restore the filtered frequency domain multi-component VSP data to the time-depth domain using two-dimensional inverse Fourier transform, eliminate window boundary distortion using the overlap-addition method, calibrate the data amplitude and baseline, and obtain the reconstructed multi-component VSP data.

2. The method for suppressing interference in the VSP multi-component polarization and frequency-wavenumber domain in metal mining areas according to claim 1, characterized in that, The dynamic adjustment of the moving time window includes the following steps: Set minimum and maximum window length thresholds, perform moving average smoothing on the local signal-to-noise ratio of the preprocessed multi-component VSP data, calculate the real-time window length using an exponential decay function, and increase the window length according to the power law of the exponential decay function based on the decrease in signal-to-noise ratio, and the real-time window length is always between the set minimum and maximum window length thresholds.

3. The method for suppressing interference in the VSP multi-component polarization and frequency-wavenumber domain in metal mining areas according to claim 1, characterized in that, The covariance matrix analysis employs Gaussian weighting factor optimization. The standard deviation of the Gaussian weighting factor is proportional to the window length. The weight value is calculated by the distance between each data point within the window and the center of the window. The greater the distance, the smaller the weight value. In other words, the weight is determined by combining the natural exponential function with the correlation between the distance from the data point to the center and the standard deviation.

4. The method for suppressing interference in the multi-component polarization and frequency-wavenumber domain of VSP in metal mining areas according to claim 1, characterized in that, The geological prior knowledge-corrected classification results include the following steps: Obtain the time-depth range of interfaces with known impedance values ​​higher than the preset impedance value, construct a geological prior weight matrix consistent with the dimensions of the preprocessed multi-component VSP data, in the geological prior weight matrix, the weight coefficient of the interface with the preset impedance value within the influence range is greater than 1, and the closer to the interface, the larger the weight coefficient is, the weight coefficient of the interface that is not higher than the preset impedance value is set to 1, and when the distance exceeds the set range, it is completely restored to 1. The obtained initial interference probability matrix is ​​multiplied point by point with the geological prior weight matrix to obtain the weighted interference probability; The threshold is then dynamically adjusted based on the local signal-to-noise ratio of the preprocessed multi-component VSP data; the lower the signal-to-noise ratio, the higher the threshold. The weighted interference probability is compared with the threshold. When the weighted interference probability is greater than or equal to the threshold, it is marked as an interference region; otherwise, it is marked as a non-interference region, thereby generating a time-depth domain interference mask.

5. The method for suppressing interference in the VSP multi-component polarization and frequency-wavenumber domain in metal mining areas according to claim 1, characterized in that, The overlap rate of the overlapping window technology is set according to a preset ratio, and the preprocessed multi-component VSP data is symmetrically extended before the two-dimensional Fourier transform. The extension length is proportional to the window length. After the two-dimensional Fourier transform, the window data is spliced ​​by a weighted average method. The weights of the overlapping areas of adjacent windows are assigned according to a linear distribution rule. The weights of the previous window gradually decrease and the weights of the next window gradually increase, so as to achieve smooth splicing of the window data.

6. The method for suppressing interference in the VSP multi-component polarization and frequency-wavenumber domain in metal mining areas according to claim 1, characterized in that, The joint domain optimization objective is to balance the interference suppression term and the polarization fidelity term by weighting coefficients, which are determined by cross-validation. The interference suppression term is calculated based on the amplitude weight of the frequency domain interference mask; the larger the amplitude, the higher the corresponding weight. The polarization fidelity term is calculated based on the differences in multi-component polarization parameters in the effective wave region; the differences in multi-component polarization parameters are obtained by comparing the polarization angle, ellipticity, and polarization intensity in the multi-component polarization parameters with the polarization parameters obtained by inversion from the filtered data.

7. The method for suppressing interference in the VSP multi-component polarization and frequency-wavenumber domain in metal mining areas according to claim 1, characterized in that, After the filtering process in step 5, frequency domain phase calibration is also included: The phase shift of the filtered frequency domain data is calculated based on the effective wave polarization angle, and then the phase shift is corrected by a linear phase compensation method. The compensation coefficient in the compensation method is proportional to the polarization angle difference. The filtered frequency domain data is calibrated point by point by a phase correction factor.

8. The method for suppressing interference in the VSP multi-component polarization and frequency-wavenumber domain in metal mining areas according to claim 1, characterized in that, The amplitude calibration includes the following steps: The ratio of the maximum amplitude of the preprocessed multi-component VSP data to the maximum amplitude of the reconstructed multi-component VSP data is calculated. Each sampling point of the reconstructed multi-component VSP data is multiplied by the ratio to ensure that the amplitude range of the reconstructed multi-component VSP data is consistent with that of the preprocessed multi-component VSP data. Furthermore, the window parameters for the overlap-addition method elimination are completely consistent with the window parameters for the frequency-wavenumber transformation.

Images