Method and program product for multiple constraint adaptive noise reduction of receiver functions at a fixed station

Abstract

The application discloses a fixed station receiver function multi-constraint adaptive noise reduction method and a program product, and the method comprises the following steps: grouping receiver functions corresponding to events of a target fixed station and adjacent stations; performing time-frequency transformation on receiver functions in each group; calculating local time-frequency signal-to-noise ratio by using a P-wave arrival time noise window based on original complex time-frequency coefficients obtained through the transformation; constructing joint reservation weights according to the local time-frequency signal-to-noise ratio and multiple constraint terms; performing time-frequency connected domain clustering on time-frequency units corresponding to the weights; and performing inverse time-frequency transformation on the original complex time-frequency coefficients according to soft weights of noise dominant clusters, effective clusters and uncertain protection clusters divided out, so as to obtain receiver functions after noise reduction. The application realizes faithful noise reduction of fixed station receiver functions by multi-constraint adaptive time-frequency clustering and soft weight reconstruction, while significantly suppressing noise, accurately protecting weak converted waves, multiple waves and anisotropic signals.

Description

Technical Field

[0001] This invention relates to the field of seismic signal processing technology, and in particular to a multi-constraint adaptive noise reduction method and program product for fixed station receiver functions. Background Technology

[0002] The receiver function method (RFM) studies the crustal and upper mantle structure beneath a station by using Ps-transformed waves and their multiples generated by teleseismic P-waves at velocity discontinuities below the station. Traditional RFM processing typically includes teleseismic event screening, three-component waveform preprocessing, coordinate rotation, deconvolution, filtering, RFM stacking, H-κ analysis, harmonic decomposition, or CCP imaging.

[0003] For fixed stations, receiver function (RSF) data are typically obtained from the same station under different teleseismic events, different inverse azimuth angles, and different ray parameters. This type of data has the following characteristics: 1. Significant differences in the signal-to-noise ratio (SNR) of single-event RSFs; 2. Real variations in the Ps converted phase, harmonic phase, and tangential components under different inverse azimuth angles; 3. Tangential RSFs may contain anisotropic, tilted interface, or three-dimensional structural information; 4. Neighboring stations near the target station can provide spatial consistency constraints; 5. Fixed bandpass filtering and artificial thresholds are difficult to adapt to different events, azimuths, frequencies, and noise levels.

[0004] Existing receiver function (RFC) denoising methods mainly include bandpass filtering, wavelet thresholding, empirical mode decomposition (EMD), variational mode decomposition (VMD), singular spectrum analysis (SMD), singular value decomposition (SVD), and deep learning denoising. While these methods can suppress noise to some extent, they still have the following shortcomings: First, fixed filtering parameters lack adaptability. Traditional RFC processing often uses fixed bandpasses, such as 0.05–1.5Hz or 0.03–2Hz. This type of fixed filtering method is difficult to adaptively adjust the denoising intensity for each event, frequency, azimuth, or noise level. Second, single-channel denoising stability is insufficient. If denoising is performed independently on a single RFC of the target station, random noise is easily misjudged as valid signals, or weak Ps transition phases, multiples, and deep interface responses are misjudged as noise. Third, tangential components are easily mistakenly deleted. Traditional RFC interpretation usually focuses on radial components, while tangential components are often regarded as noise or secondary information. However, when azimuth anisotropy, tilted interfaces, or three-dimensional scattering structures exist, the tangential component may contain key structural information; fourth, the inverse azimuth information is underutilized. If all inverse azimuth data are mixed and processed, the amplitude, polarity, and arrival time characteristics of the receiver function as a function of azimuth may be weakened, thus affecting anisotropy analysis and tilted interface identification; fifth, hard thresholding denoising is prone to causing loss of effective signal. Wavelet thresholding, spectral subtraction, or binary time-frequency masks often directly clear low-energy regions, which may cause the azimuth-related energy in weak transition phases, multiples, or tangential components to be mistakenly deleted, and may introduce phase distortion; sixth, there is a lack of interpretable quality control output. Existing denoising methods usually only output the denoised waveform, making it difficult to determine which time-frequency regions are preserved, which are suppressed, and whether the difference receiver function contains the mistakenly deleted effective phase.

[0005] Therefore, there is an urgent need to invent an adaptive noise reduction method for received function data from fixed stations, in order to solve the problems of limited signal-to-noise ratio improvement, insufficient noise reduction stability of single station and single channel in existing noise reduction methods, and the inability to retain the Ps converted phase, multiple wave, tangential component and azimuth anisotropy characteristics while improving the signal-to-noise ratio. Summary of the Invention

[0006] In view of this, embodiments of the present invention provide a multi-constraint adaptive noise reduction method and program product for fixed station receiving functions, which at least partially solves the problems existing in the prior art.

[0007] Other features and advantages of the invention will become apparent from the following detailed description, or may be learned in part by practice of the invention.

[0008] To achieve the above objectives, the embodiments of the present invention provide the following technical solutions:

[0009] According to a first aspect of the present invention, a multi-constraint adaptive noise reduction method for a fixed station receiving function is provided, the method comprising:

[0010] Acquire multi-event teleseismic three-component seismic records from the target fixed station and nearby stations;

[0011] For the event seismic records corresponding to each event at each station, data preprocessing, coordinate rotation, and deconvolution are performed on the event seismic records to obtain the event reception function, which includes a radial reception function and a tangential reception function.

[0012] The event receiving functions are grouped according to the reverse azimuth angle and ray parameters, and time-frequency transformation is performed on the event receiving functions in each event group to obtain the original complex time-frequency coefficients.

[0013] Based on the original complex time-frequency coefficients, the frequency-related noise level is estimated using the P-wave arrival time-precedence noise window, and the local time-frequency signal-to-noise ratio is calculated based on the frequency-related noise level.

[0014] The radial component score and tangential component score are calculated based on the local time-frequency signal-to-noise ratio. The radial component score and the tangential component score are used at least and fused with at least one of the following: radial / tangential local coherence, spatial consistency of neighboring stations, azimuth consistency and azimuth harmonic protection, to obtain a joint retention weight.

[0015] Time-frequency connected component clustering is performed on the time-frequency units corresponding to the joint retention weights to divide them into noise-dominated clusters, effective clusters, and uncertain protection clusters;

[0016] Based on the soft weights of the noise-dominant cluster, the effective cluster, and the uncertain protection cluster, the original complex time-frequency coefficients are subjected to inverse time-frequency transformation to obtain the noise-reduced receiver function.

[0017] Furthermore, multi-event teleseismic three-component seismic records from the target fixed station and neighboring stations were acquired, including:

[0018] According to the target fixed station The target fixed station is determined by the preset spatial neighborhood radius. The neighboring stations, obtain the set of neighboring stations. ,in, This represents the surface distance between stations. The preset spatial neighborhood radius;

[0019] Acquire multi-event teleseismic three-component seismic records from the target fixed station and neighboring stations.

[0020] Furthermore, the event receiving functions are grouped according to the inverse azimuth angle and ray parameters, and a time-frequency transformation is performed on the event receiving functions within each event group to obtain the original complex time-frequency coefficients, including:

[0021] For each station For each event e, the inverse azimuth angle corresponding to event e is used. and ray parameters Group the event receiving functions to obtain event groups. ;

[0022] The event group The expression is ,in, Center of the preset orientation window The preset orientation window half-width, Center the preset ray parameter window. This is the preset allowable range for ray parameters;

[0023] For each event group The event receiving function within the function performs time-frequency transformation to obtain the original complex time-frequency coefficients. The expression for the original complex time-frequency coefficients is as follows: Where c represents the radial component or the tangential component, Indicates radial or tangential receiver function. This indicates a time-frequency transform operator, which is either a continuous wavelet transform operator or an S-transform operator;

[0024] The original complex time-frequency coefficients include radial complex time-frequency coefficients and tangential complex time-frequency coefficients.

[0025] Further, based on the original complex time-frequency coefficients, the frequency-dependent noise level is estimated using a P-wave arrival time-premature noise window, and the local time-frequency signal-to-noise ratio is calculated according to the frequency-dependent noise level, including:

[0026] The time window preceding the arrival of the P-wave is used as the noise estimation window, and the expression for the noise estimation window is as follows: ,in, The time when the P-wave arrives at the station. and To preset the noise window parameters, ;

[0027] Within the noise estimation window, the frequency-dependent noise level is estimated using the original complex time-frequency coefficients. The expression for the frequency-related noise level is: ,in, These are the original complex time-frequency coefficients. ;

[0028] Calculate the local time-frequency signal-to-noise ratio using the frequency-dependent noise level. The calculation formula is: ,in, To prevent division by zero constant;

[0029] The local time-frequency signal-to-noise ratio includes the radial local time-frequency signal-to-noise ratio and the tangential local time-frequency signal-to-noise ratio.

[0030] Further, radial component scores and tangential component scores are calculated based on the local time-frequency signal-to-noise ratio. At least the radial component scores and tangential component scores are used, and fused with at least one of the following: radial / tangential local coherence, spatial consistency of neighboring stations, azimuth consistency, and azimuth harmonic protection, to obtain a joint retention weight, including:

[0031] Based on the radial local time-frequency signal-to-noise ratio, the radial component score is calculated using the following formula: ,in, For radial local time-frequency signal-to-noise ratio, Represents the sigmoid function;

[0032] Based on the tangential local time-frequency signal-to-noise ratio, the tangential component score is calculated using the following formula:

[0033] ,in, For tangential local time-frequency signal-to-noise ratio;

[0034] In the neighborhood of the preset time frequency Internal calculation of radial / tangential local coherence The calculation formula is: ,in, Represents the radial complex time-frequency coefficient. Represents the tangential complex time-frequency coefficients. Represents the complex conjugate of the tangential complex time-frequency coefficient;

[0035] Based on the radial component score and the tangential component score, the spatial consistency of neighboring stations is calculated. The calculation formula is: ,in, , For indicator functions, Set a preset candidate retention threshold;

[0036] azimuth consistency is calculated based on the radial component score and the tangential component score. The calculation formula is: ;

[0037] The azimuth harmonic protection term is obtained by harmonic decomposition of the event receiving function as a function of the reverse azimuth angle.

[0038] The joint retention weight is obtained by fusing at least the radial component score and the tangential component score, and combining them with at least one of the following: radial / tangential local coherence, spatial consistency of neighboring stations, azimuth consistency, and azimuth harmonic protection. .

[0039] Further, the radial component score and tangential component score are calculated based on the local time-frequency signal-to-noise ratio. At least the radial component score and the tangential component score are used, and fused with at least one of the following: radial / tangential local coherence, spatial consistency of neighboring stations, azimuth consistency, and azimuth harmonic protection, to obtain a joint retention weight. This also includes:

[0040] If there are no neighboring stations within the preset spatial neighborhood radius of the target fixed station, or if the data quality of neighboring stations does not meet the preset conditions, then the repeatability of the same station event is calculated using the set of similar reverse azimuth and similar ray parameter events within the target fixed station. Instead of spatial consistency at neighboring stations, the expression for the repetition of events at the same station is: ,in, This represents the set of events that have similar reverse azimuth angles and similar ray parameters to event e;

[0041] At least the radial component score and the tangential component score are used, and at least one of the following is fused: radial / tangential local coherence, repetitiveness of co-occurrence events, azimuth consistency, and azimuth harmonic protection terms, to obtain the joint retention weight. The fusion includes linear fusion, product-type fusion, or logistic fusion.

[0042] Furthermore, time-frequency connected component clustering is performed on the time-frequency units corresponding to the joint retention weights to divide them into noise-dominated clusters, effective clusters, and uncertain protection clusters, including:

[0043] Based on the joint retention weights, a candidate binary mask is constructed according to a preset clustering threshold, and the set of candidate time-frequency units to be discriminated is determined according to the candidate binary mask;

[0044] Based on the candidate binary mask, candidate time-frequency unit pairs with a preset adjacency relationship are connected to obtain candidate time-frequency clusters to be judged. The preset adjacency relationship is a four-neighborhood, an eight-neighborhood, a frequency-weighted neighborhood, a time-weighted neighborhood, or a graph-connected neighborhood.

[0045] For each candidate time-frequency cluster q, the confidence score corresponding to the candidate time-frequency cluster is calculated based on the cluster energy, tangential / radial energy ratio, neighbor station repeatability, cluster orientation consistency, and theoretical phase time window matching degree.

[0046] The formula for calculating the cluster energy is: ;

[0047] The formula for calculating the tangential / radial energy ratio is: ;

[0048] The formula for calculating the repeatability of neighboring stations is: ;

[0049] The formula for calculating the cluster orientation consistency is as follows: ;

[0050] The formula for calculating the theoretical phase time window matching degree is as follows: ,in, For the theoretical phase arrival time, The center time of the candidate time-frequency cluster, This represents the theoretical phase time window width.

[0051] Determine whether the confidence score of the candidate time-frequency cluster is greater than or equal to a preset effective cluster threshold;

[0052] If the confidence score of the candidate time-frequency cluster is greater than or equal to the preset effective cluster threshold, then the candidate time-frequency cluster is regarded as an effective cluster.

[0053] If the confidence score of the candidate time-frequency cluster is less than the preset effective cluster threshold, then it is determined whether the candidate time-frequency cluster meets at least one uncertain protection condition. The uncertain protection conditions include preset harmonic protection conditions, preset theoretical phase window matching conditions, preset neighbor station repeatability conditions, and co-station multi-event repeatability conditions.

[0054] If the candidate time-frequency cluster satisfies at least one uncertainty protection condition, then the candidate time-frequency cluster is designated as an uncertainty protection cluster.

[0055] If the candidate time-frequency cluster does not meet any of the uncertain protection conditions, then the candidate time-frequency cluster is regarded as the noise-dominant cluster.

[0056] Further, based on the soft weights of the noise-dominant cluster, the effective cluster, and the uncertain protection cluster, an inverse time-frequency transform is performed on the original complex time-frequency coefficients to obtain the denoised receiver function, including:

[0057] Constructing soft weight functions ,in, For effective cluster regions The corresponding soft weights, For uncertain protected cluster regions The corresponding soft weights, Noise-dominant cluster region The corresponding soft weights are as follows: the soft weight corresponding to the effective cluster region is greater than the soft weight corresponding to the uncertain protection cluster region, the soft weight corresponding to the uncertain protection cluster region is greater than the soft weight corresponding to the noise-dominant cluster region, and the soft weight corresponding to the noise-dominant cluster region is greater than or equal to zero.

[0058] The original complex time-frequency coefficients are subjected to inverse time-frequency transformation using the soft weighting function to obtain the noise-reduced receiving function. The inverse time-frequency transformation expression is as follows: Among them, the noise-reduced receiver function This includes the denoised radial receiver function and the denoised tangential receiver function. These are the weighted original complex time-frequency coefficients. , It is an inverse time-frequency transform operator.

[0059] Furthermore, the method also includes:

[0060] Using the noise-reduced receiving function and the original event receiving function Calculate the difference receiving function The calculation formula is: ;

[0061] Based on the noise-reduced receiver function and the original event receiver function, quality control indicators are calculated. These quality control indicators include Ps time offset, amplitude retention rate, tangential component energy retention rate, harmonic term retention rate, signal-to-noise ratio improvement rate, time-frequency energy retention rate, and effective phase residue index of the difference receiver function.

[0062] Determine whether the quality control indicators meet the preset quality threshold;

[0063] If the quality control indicators do not meet the preset quality threshold, the joint retention weight parameter, preset clustering threshold, preset effective cluster threshold, azimuth harmonic protection threshold, theoretical phase time window width, or soft weight level will be adjusted accordingly.

[0064] According to a second aspect of the present invention, a computer program product is provided, the computer program product comprising a computing program stored on a non-transitory computer-readable storage medium, the computer program comprising program instructions that, when executed by a computer, cause the computer to implement the steps of the fixed station receiving function multi-constraint adaptive noise reduction method as described in any of the preceding claims.

[0065] This invention provides a multi-constraint adaptive noise reduction method and program product for fixed station receiver functions. The method includes: grouping the receiver functions corresponding to event seismic records of the target fixed station and neighboring stations; performing time-frequency transformation on the receiver functions within each group; calculating the local time-frequency signal-to-noise ratio (SNR) using a P-wave arrival time-precedence noise window based on the original complex time-frequency coefficients obtained from the transformation; constructing joint retention weights based on the local time-frequency SNR and various constraint terms; performing time-frequency connected component clustering on the time-frequency units corresponding to the weights; and performing inverse time-frequency transformation on the original complex time-frequency coefficients based on the soft weights of the identified noise-dominant clusters, effective clusters, and uncertain protection clusters to obtain the denoised receiver function. This invention, through multi-constraint adaptive time-frequency clustering and soft-weight reconstruction, significantly suppresses noise while accurately protecting weak converted waves, multiple waves, and anisotropic signals, achieving fidelity-preserving noise reduction for fixed station receiver functions. Attached Figure Description

[0066] To more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are merely exemplary, and those skilled in the art can derive other embodiments based on the provided drawings without creative effort.

[0067] Figure 1 This is a flowchart illustrating a multi-constraint adaptive noise reduction method for a fixed station receiving function provided in an embodiment of the present invention.

[0068] Figure 2 A schematic diagram illustrating the process principle of a multi-constraint adaptive noise reduction method for a fixed station receiving function provided in an embodiment of the present invention;

[0069] Figure 3 This is a schematic diagram of the spatial distribution of stations provided in an embodiment of the present invention;

[0070] Figure 4 This is a schematic diagram of reverse azimuth angle grouping provided in an embodiment of the present invention;

[0071] Figure 5 This is an exemplary schematic diagram of the target fixed station receiving function waveform and time-spectrum noise reduction provided in an embodiment of the present invention;

[0072] Figure 6 This is a schematic diagram of the soft weight reconstruction, noise reduction result output, and quality control feedback process provided in an embodiment of the present invention. Detailed Implementation

[0073] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0074] It should be noted that the terms "first," "second," etc., in the specification, claims, and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of the invention described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover a non-exclusive inclusion; for example, a process, method, system, product, or apparatus that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0075] Figure 1 A flowchart of a multi-constraint adaptive noise reduction method for a fixed station receiving function according to an embodiment of the present invention is shown.

[0076] Figure 2 A schematic diagram illustrating the process principle of a multi-constraint adaptive noise reduction method for a fixed station receiving function according to an embodiment of the present invention is shown.

[0077] like Figure 1 As shown, the fixed station receiving function multi-constraint adaptive noise reduction method according to the embodiment of the present invention may include steps S100, S200, S300, S400, S500, S600 and S700.

[0078] In step S100, multi-event teleseismic three-component seismic records of the target fixed station and neighboring stations are acquired.

[0079] Specifically, the above steps include:

[0080] Based on the target fixed station The target fixed station is determined by the preset spatial neighborhood radius. The neighboring stations, obtain the set of neighboring stations. ,in, This represents the surface distance between stations. The preset spatial neighborhood radius.

[0081] Acquire multi-event teleseismic three-component seismic records from the target fixed station and neighboring stations.

[0082] Figure 3 A schematic diagram of the spatial distribution of stations provided in an embodiment of the present invention is shown.

[0083] Next, in step S200, for the event seismic records corresponding to each event at each station, data preprocessing, coordinate rotation and deconvolution are performed on the event seismic records to obtain the event receiving function, which includes a radial receiving function and a tangential receiving function.

[0084] Specifically, the above steps include:

[0085] For the event seismic records corresponding to each event at each station, the vertical components are obtained after preprocessing and coordinate rotation. radial component and tangential components .

[0086] In one implementation, the event reception function is calculated using frequency domain level factor deconvolution, and the radial reception function is expressed as: The tangential receiving function is expressed as: ,in, Angular frequency, Represents the complex conjugate of the vertical component spectrum. It is either a level factor or a stability factor.

[0087] In another implementation, the receiver function can also be obtained using an iterative deconvolution method, i.e., finding... and , making , , where * represents the convolution operation.

[0088] It should be noted that this invention is not limited to a specific deconvolution algorithm, and can be applied to frequency domain level factor deconvolution, time domain iterative deconvolution, multi-window spectrum deconvolution, or other receiver function calculation methods.

[0089] In step S300, the event receiving functions are grouped according to the reverse azimuth angle and ray parameters, and the event receiving functions in each event group are subjected to time-frequency transformation to obtain the original complex time-frequency coefficients.

[0090] Specifically, the above steps include:

[0091] For each station For each event e, the inverse azimuth angle corresponding to event e is used. and ray parameters Group the event receiving functions to obtain event groups. .

[0092] The above event group The expression is ,in, Center of the preset orientation window The preset orientation window half-width, Center the preset ray parameter window. This is the preset allowable range for ray parameters.

[0093] Figure 4 A schematic diagram of reverse azimuth angle grouping provided in an embodiment of the present invention is shown.

[0094] For sliding azimuth windows, preset window width and preset step size can be used; for data with a small number of events, a wider azimuth window can be used to ensure statistical stability. The specific window width and step size can be adaptively set according to the number of events, research objectives, and azimuth coverage.

[0095] For each event group The event receiving function within the function performs time-frequency transformation to obtain the original complex time-frequency coefficients, which include radial complex time-frequency coefficients and tangential complex time-frequency coefficients.

[0096] The expression for the above-mentioned primitive complex time-frequency coefficients is: Where c represents the radial component or the tangential component, Indicates radial or tangential receiver function. This represents a time-frequency transform operator, which can be a continuous wavelet transform operator or an S-transform operator.

[0097] The above continuous wavelet transform operator is defined as follows: Where a is the scale parameter and b is the time translation parameter. The mother wavelet. Its scale and frequency can be approximated as: ,in, The center frequency of the mother wavelet. The sampling interval is denoted as .

[0098] The reason for using time-frequency conversion in this invention is that the Ps conversion phase and multiple waves in the receiving function have non-stationary transient characteristics, and fixed bandpass filtering is difficult to adapt to both short-time high-frequency energy and long-time low-frequency energy at the same time.

[0099] In other implementations, the time-frequency transformation may also employ short-time Fourier transform, synchronous compressed wavelet transform, matched pursuit, or other transformation methods that yield a time-frequency representation.

[0100] In step S400, based on the original complex time-frequency coefficients, the frequency-dependent noise level is estimated using the P-wave arrival time-precedence noise window, and the local time-frequency signal-to-noise ratio is calculated based on the frequency-dependent noise level.

[0101] Specifically, the above steps include:

[0102] Using the time window preceding the arrival of the P-wave as the noise estimation window, the expression for the noise estimation window is as follows: ,in, The time when the P-wave arrives at the station. and To preset the noise window parameters, The specific value can be set according to the sampling rate, data length and event noise level.

[0103] Within the noise estimation window, the frequency-dependent noise level is estimated using the original complex time-frequency coefficients. The expression for the frequency-dependent noise level is: ,in, These are the original complex time-frequency coefficients. .

[0104] Calculate the local time-frequency signal-to-noise ratio using the frequency-dependent noise level. The aforementioned local time-frequency signal-to-noise ratio includes radial local time-frequency signal-to-noise ratio and tangential local time-frequency signal-to-noise ratio.

[0105] The formula for calculating the local time-frequency signal-to-noise ratio is as follows: ,in, To prevent the removal of the zero constant.

[0106] The noise estimation method used in this invention does not rely on manually specifying a global threshold and can adapt to noise levels under different events, locations, and frequencies.

[0107] In other implementations, noise estimation may also employ quantile estimation, local block variance estimation, robust root mean square estimation, or Bayesian noise estimation.

[0108] In step S500, the radial component score and the tangential component score are calculated based on the local time-frequency signal-to-noise ratio. At least the radial component score and the tangential component score are used, and at least one of the radial / tangential local coherence, spatial consistency of neighboring stations, azimuth consistency and azimuth harmonic protection terms are fused to obtain the joint retention weight.

[0109] Specifically, the above steps include:

[0110] The radial component score is calculated based on the radial local time-frequency signal-to-noise ratio, using the following formula: ,in, For radial local time-frequency signal-to-noise ratio, Representing the sigmoid function: , To indicate the steepness of the function, This is the center threshold for the signal-to-noise ratio.

[0111] Based on the tangential local time-frequency signal-to-noise ratio, the tangential component score is calculated using the following formula:

[0112] ,in, This represents the tangential local time-frequency signal-to-noise ratio.

[0113] In the neighborhood of the preset time frequency Internal calculation of radial / tangential local coherence The calculation formula is: ,in, Represents the radial complex time-frequency coefficient. Represents the tangential complex time-frequency coefficients. It represents the complex conjugate of the tangential complex time-frequency coefficient.

[0114] Spatial consistency of neighboring stations is calculated based on radial and tangential component scores. The calculation formula is: ,in, , For indicator functions, Set the preset candidate retention threshold.

[0115] If there are no neighboring stations within the preset spatial neighborhood radius of the target fixed station, then the repeatability of similar reverse azimuth and similar ray parameter events within the target fixed station is utilized. Instead of spatial consistency at neighboring stations, the expression for the repetition of events at the same station is: ,in, This represents the set of events that have similar reverse azimuth angles and similar ray parameters to event e.

[0116] azimuth consistency is calculated based on the radial component score and the tangential component score. The calculation formula is: .

[0117] If there are azimuth anisotropy, tilted interfaces, or three-dimensional non-uniform structures beneath a fixed station, the receiver function may exhibit periodic changes with the reverse azimuth angle. To avoid accidentally deleting weak azimuth-related signals, this invention introduces an azimuth harmonic protection mechanism.

[0118] Azimuth harmonic decomposition can be applied to the radial receiver function, the tangential receiver function, or the characteristic quantity obtained by combining the radial and tangential components.

[0119] For the target fixed station The change of the event receiving function with respect to the reverse azimuth angle can be expressed as: The first-order harmonic energy and the second-order harmonic energy are defined as follows: and The harmonic contribution ratio is defined as .

[0120] The radial component score, tangential component score, radial / tangential local coherence, spatial consistency of neighboring stations, azimuth consistency, and / or azimuth harmonic protection terms are normalized to the [0,1] interval before fusion. The joint retention weight is obtained by fusing at least one normalized constraint term. .

[0121] In one implementation, linear fusion is used for the constraint terms: and satisfy ,in, This is for azimuth harmonic protection.

[0122] If there are no neighboring stations within the preset spatial neighborhood radius of the target fixed station, then the repetition rate of the same station event is used. Replace the spatial consistency of adjacent stations The joint retention weight is then expressed as: And satisfy: .

[0123] In another implementation, the constraint terms can also be fused using a product-type approach, for example: .

[0124] In another implementation, Logistic fusion can also be used for the constraint terms, for example: ,in, , representing each normalization constraint term, , b represents parameters to be optimized or preset.

[0125] For fixed stations targeting anisotropy studies, the weights of tangential component scores, azimuth consistency, and harmonic protection terms can be increased. For fixed stations targeting Moho Ps arrival picking or H-κ superposition, the weights of radial component scores, spatial consistency with neighboring stations, and theoretical phase window correlation terms can be increased.

[0126] Preferably, the weight coefficients, exponential parameters, or logistic fusion parameters in the joint retention weights are automatically optimized based on the quality control objectives.

[0127] In one implementation, the following objective function can be used to optimize the weight coefficients, exponential parameters, or logistic fusion parameters: ,in, This represents the set of parameters to be optimized. The weights are the objective function weights.

[0128] Parameter optimization can employ methods such as grid search, Bayesian optimization, multi-objective optimization, cross-validation, or learning from a small number of labeled samples. In engineering implementation, empirical weights can be given first, and then iteratively fine-tuned based on quality control results.

[0129] In step S600, time-frequency connected component clustering is performed on the time-frequency units corresponding to the joint retention weights to divide them into noise-dominated clusters, effective clusters, and uncertain protection clusters.

[0130] Specifically, the above steps include:

[0131] Based on the joint retention weights, candidate binary masks for time-frequency units corresponding to the joint retention weights are constructed according to a preset clustering threshold: ,in, The threshold can be either a fixed threshold or an adaptive threshold.

[0132] Based on the candidate binary mask, candidate time-frequency unit pairs with a preset adjacency relationship are connected to obtain candidate time-frequency clusters. The aforementioned preset adjacency relationships are four-neighborhood, eight-neighborhood, frequency-weighted neighborhood, time-weighted neighborhood, or graph-connected neighborhood.

[0133] For each candidate time-frequency cluster q, calculate the cluster energy, tangential / radial energy ratio, neighbor station repeatability, cluster orientation consistency, and theoretical phase time window matching degree of the candidate time-frequency cluster.

[0134] The formula for calculating the cluster energy is as follows: .

[0135] The formula for calculating the tangential / radial energy ratio is as follows: .

[0136] The formula for calculating the repetition of adjacent stations is as follows: .

[0137] The formula for calculating the cluster orientation consistency is as follows: .

[0138] The formula for calculating the theoretical phase-time window matching degree is as follows: ,in, For the theoretical phase arrival times such as Ps, PpPs, PsPs+PpSs, etc. The center time of the candidate time-frequency cluster, This represents the width of the theoretical phase time window. In this invention, the theoretical phase time window serves as a positive protection constraint, not the sole basis for deletion. When the velocity model is uncertain, the time window can be widened or combined with data-driven time-of-arrival picking results.

[0139] When calculating the spatial consistency of neighboring stations, the recurrence of candidate clusters among neighboring stations allows for preset time tolerances, frequency tolerances, and azimuth tolerances. The time tolerance can be determined by theoretical phase transition time difference, station spacing, velocity model, or cross-correlation alignment results.

[0140] For a candidate time-frequency cluster q, its harmonic protection term is defined as: .like If the candidate time-frequency cluster is marked as a harmonic protection cluster, it will not be completely deleted even if the energy of the cluster is weak, but will be assigned a retention weight or a medium soft weight. This mechanism can protect the periodic polarity changes in the tangential component, the first or second order azimuth harmonic response, the tilted interface response, and the weak energy related to scattering of three-dimensional structures.

[0141] This invention does not simply delete weak clusters based on cluster area or energy, but instead employs a multi-condition retention rule constrained by seismology.

[0142] For a candidate time-frequency cluster q, its reconstruction weight is retained or increased if any of the following conditions are met: , , , , , ,in, Statistical values ​​representing radial / tangential local coherence within a candidate cluster, such as mean, median, or quantile.

[0143] Furthermore, the credibility score corresponding to the candidate time-frequency cluster is calculated using the following formula: In this case, all terms are normalized to the interval [0,1], and .

[0144] Determine whether the credibility score of the candidate time-frequency cluster is greater than or equal to the preset effective cluster threshold.

[0145] If the credibility score of a candidate time-frequency cluster is greater than or equal to the preset effective cluster threshold, then the candidate time-frequency cluster is considered an effective cluster.

[0146] If the credibility score of the candidate time-frequency cluster is less than the preset effective cluster threshold, it is determined whether the candidate time-frequency cluster meets at least one uncertain protection condition. The aforementioned uncertain protection conditions include preset harmonic protection conditions, preset theoretical phase window matching conditions, preset neighbor station repeatability conditions, and co-station multi-event repeatability conditions.

[0147] If a candidate time-frequency cluster satisfies at least one uncertainty protection condition, then the candidate time-frequency cluster is designated as an uncertainty protection cluster.

[0148] If a candidate time-frequency cluster does not meet any of the uncertain protection conditions, then the candidate time-frequency cluster is regarded as the noise-dominant cluster.

[0149] The core of the above criterion is that weak signals will not be directly deleted simply because they are weak, as long as they have seismic consistency in neighboring stations, azimuth, theoretical phase, radial / tangential coherence, or harmonic terms.

[0150] In step S700, the original complex time-frequency coefficients are subjected to inverse time-frequency transformation based on the soft weights of the noise-dominant cluster, the effective cluster, and the uncertain protection cluster to obtain the noise-reduced receiver function.

[0151] Specifically, the above steps include:

[0152] Traditional hard mask reconstruction is usually represented as: In this case, M can only be 0 or 1, and traditional methods may cause phase distortion and false deletion of weak phases.

[0153] This invention employs soft-weight reconstruction to construct a soft-weight function. ,in, For effective cluster regions The corresponding soft weights, For uncertain protected cluster regions The corresponding soft weights, Noise-dominant cluster region The corresponding soft weights, .

[0154] In one implementation, The above values ​​are for illustrative purposes only and can be adjusted adaptively based on research objectives, noise levels, and quality control results.

[0155] For anisotropy studies, the preferred method is... This is to avoid completely removing weak azimuth-related energy in the tangential component. Even if a candidate region is determined to be a noise-dominant region, the present invention still assigns it a soft weight greater than zero, so that the region is attenuated rather than completely cleared.

[0156] In another implementation, the soft weights can also be obtained by continuous mapping of cluster credibility scores, for example: ,in, , , To map the steepness, This is the center value of the score.

[0157] This invention applies soft weights to the original complex time-frequency coefficients, rather than simply reconstructing the time-frequency amplitudes. This preserves the phase information of the original receiver function and reduces arrival offset and waveform distortion caused by hard thresholding or phase reestimation.

[0158] Using a soft weighting function to perform an inverse time-frequency transform on the original complex time-frequency coefficients, we obtain the denoised receiver function. The expression for the inverse time-frequency transform is: Among them, the noise-reduced receiver function This includes the denoised radial receiver function and the denoised tangential receiver function. These are the weighted original complex time-frequency coefficients. , It is an inverse time-frequency transform operator.

[0159] In addition, embodiments of the present invention also include:

[0160] Using the noise-reduced receiver function and the original event receiving function Calculate the difference receiving function The calculation formula is: .

[0161] To ensure the seismological reliability of the noise reduction results, quality control indices are calculated based on the noise-reduced receiver function and the original event receiver function. These indices include Ps arrival time migration, amplitude retention rate, tangential component energy retention rate, harmonic term retention rate, signal-to-noise ratio improvement rate, time-frequency energy retention rate, and effective phase residual index of the difference receiver function.

[0162] The expression for the time offset of Ps above is: ,like If the set threshold is exceeded, a message will be displayed indicating that the noise reduction may have altered the effective phase.

[0163] The expression for the amplitude retention rate is as follows: This indicator is used to determine whether noise reduction is excessively suppressed or abnormally amplified by the Ps conversion phase.

[0164] The expression for the energy retention rate of the tangential component is as follows: ,in, The focus time window is set for the tangential component. For anisotropic studies, the energy retention rate of the tangential component should not be abnormally low. It should be noted that... This represents the quality control indicator, not the tangential receiving function. The tangential receiving function is uniformly denoted as... .

[0165] The expression for the harmonic term retention rate is as follows: This indicator is used to determine whether information related to azimuth anisotropy has been excessively weakened.

[0166] The expression for the above signal-to-noise ratio improvement rate is: This metric is used to evaluate the degree of improvement in the signal-to-noise ratio of the effective phase window relative to the noise window after noise reduction.

[0167] The expression for the above time-frequency energy retention rate is: This metric is used to evaluate the final proportion of time-frequency energy retained, avoiding excessive noise reduction.

[0168] The above-mentioned differential receiver function inspection index is defined as follows: Ideally, the differential receiver function should mainly exhibit low-amplitude random noise residuals. If clear Ps, multiple waves, or tangential energy with azimuth periodicity appear in the differential receiver function, it indicates that noise reduction may have mistakenly deleted effective signals.

[0169] The expression for the effective phase vestigial index of the above differential receiver function is: ,in, For the Ps-conversion phase or multiple wave theory phase window, if An excessively high value indicates the presence of effective phase residue in the differential receiver function, suggesting that the noise reduction is too strong.

[0170] Furthermore, it is determined whether the quality control indicators meet the preset quality thresholds.

[0171] If the quality control indicators do not meet the preset quality threshold, the joint retention weight parameters, preset clustering threshold, preset effective cluster threshold, azimuth harmonic protection threshold, theoretical phase time window width, or soft weight level will be adjusted accordingly.

[0172] Specifically, when the quality control indicators do not meet the preset quality thresholds, feedback adjusts at least one of the following parameters: 1. Weight coefficients in the joint retention weights; 2. Exponential parameters in product-type fusion; 3. Bias terms or weight terms in Logistic fusion; 4. Time-frequency clustering threshold; 5. Candidate cluster credibility score threshold; 6. Soft weight levels corresponding to effective clusters, uncertain protection areas, and noise-dominant areas; 7. Azimuth harmonic protection threshold; 8. Theoretical phase time window width; 9. Time tolerance, frequency tolerance, or azimuth tolerance in the repeatability judgment of neighboring stations. This feedback mechanism allows the quality control output not only to be displayed as a result but also as a basis for adaptive optimization of noise reduction parameters.

[0173] Figure 5 This diagram illustrates an exemplary schematic of the target fixed station receiving function waveform and time-spectrum noise reduction provided in an embodiment of the present invention, wherein... Figure 5 (a) shows the waveform of the original event receiving function. Figure 5 (b) is the original time spectrum. Figure 5 (c) is a soft-weighted time-frequency mask. Figure 5 (d) is the waveform of the receiver function after noise reduction. Figure 5 (e) is the difference receiving function. Figure 5 This demonstrates that soft-weight reconstruction can suppress noise while preserving the effective transition phase.

[0174] Figure 6 This diagram illustrates the soft weight reconstruction, noise reduction result output, and quality control feedback process provided in an embodiment of the present invention.

[0175] In addition, embodiments of the present invention also provide a computer program product, which includes computer program instructions that, when executed by a processor, implement the steps of a fixed station receiving function multi-constraint adaptive noise reduction method as described above.

[0176] The following examples are used to further illustrate the present invention, but do not constitute a limitation on the scope of protection of the present invention.

[0177] Example 1: Noise Reduction of Target Fixed Station Receiver Function

[0178] Select a fixed broadband seismic station It collects teleseismic events whose epicentral distance and magnitude meet preset conditions.

[0179] For each event, the three-component records are preprocessed by removing the mean, trend, instrument response, resampling, and bandpass, and then rotated to the ZRT coordinate system.

[0180] Using the vertical component as the source function, deconvolution is performed on the radial and tangential components to obtain: .

[0181] The event is divided into several groups based on the reverse azimuth angle and ray parameters. A continuous wavelet transform is performed on the receiver function within each group, and the noise level is estimated using a preset time window on the P-wavefront to calculate the local time-frequency signal-to-noise ratio.

[0182] Subsequently, a radial / tangential multi-component joint retention weight is constructed, and connected component clustering is performed in the time-frequency domain. High soft weights are assigned to effective clusters, medium soft weights to uncertain protection regions, and low soft weights (greater than zero) to noise-dominant regions. Finally, an inverse transform is performed to obtain the denoised radial and tangential receiver functions.

[0183] Example 2: Collaborative Noise Reduction between Target Fixed Station and Neighboring Stations

[0184] At the target station The surrounding area selected is less than the preset spatial neighborhood radius. Collection of neighboring stations For neighboring stations, the same preprocessing, rotation, and receiving function calculation parameters as the target station are used.

[0185] For a candidate time-frequency cluster in the target station, determine whether the cluster recurs in neighboring stations with similar reverse azimuth angles, similar ray parameters, and similar time-frequency positions. If the recurrence rate exceeds a threshold, increase the spatial consistency score of the cluster.

[0186] When it is determined that a neighboring station is repeated, preset time tolerance, frequency tolerance, and azimuth tolerance can be allowed, or corrections can be made based on theoretical phase transition time difference, velocity model, or cross-correlation alignment results.

[0187] This implementation method can improve the identification stability of weak Ps transition phase and multiple waves of the target station, and reduce the impact of single-event random noise on the noise reduction results.

[0188] Example 3: Protection Treatment for Anisotropic Fixed Stations

[0189] For fixed stations suspected of having azimuth anisotropy or tilted interfaces beneath them, harmonic decomposition is performed on the radial and tangential receiving functions at different reverse azimuth angles.

[0190] If a candidate time-frequency cluster makes a significant contribution to the first or second harmonic terms, it is marked as a harmonic protection cluster and given a retention weight or a medium soft weight, even if its time-frequency energy is weak.

[0191] This processing can prevent periodic polarity changes, weak azimuth responses, or tilted interface responses in the tangential component from being misjudged as noise.

[0192] Example 4: Adaptive Weight Learning for Different Research Objectives

[0193] When the research objective is Moho Ps arrival picking or H-κ superposition, weight learning can improve the contributions of radial component scores, theoretical phase window matching, and neighbor station consistency terms, for example, the following can be obtained: .

[0194] When the research objective is receiver function harmonic decomposition or anisotropy identification, weight learning can improve the contributions of tangential component scores, azimuth consistency, and harmonic protection terms. For example, it can yield: .

[0195] When the research objective is CCP imaging or fixed network regional structure analysis, weight learning can improve the contribution of spatial consistency of neighboring stations and repeatability terms within azimuth groups.

[0196] The weights mentioned above are for illustrative purposes only. Actual weights are automatically determined by the data noise level, number of events, distribution of neighboring stations, and research objectives.

[0197] Example 5: Quality Control Feedback Noise Reduction

[0198] After one noise reduction step, the system calculates the time offset of Ps, amplitude retention rate, tangential component energy retention rate, harmonic term retention rate, signal-to-noise ratio improvement rate, time-frequency energy retention rate, and effective phase residual index of the difference receiver function.

[0199] If the Ps arrival time offset exceeds the threshold, the amplitude retention rate is significantly low, the tangential component energy retention rate is abnormally low, the harmonic term retention rate is significantly reduced, or clear Ps, multiples, or azimuth periodic tangential energy appear in the difference receiver function, the system will reduce the clustering screening intensity, increase the soft weight of uncertain protection areas or noise-dominated areas, or increase the weight of tangential components and harmonic protection terms, and reconstruct the soft weights.

[0200] Example 6: Alternative treatment under the condition of no nearby stations

[0201] When there are no neighboring stations that meet the conditions around the target station, the spatial consistency term of neighboring stations can be replaced by the repeatability term of similar reverse azimuth and similar ray parameter events within the same station.

[0202] Specifically, the recurrence rate of candidate time-frequency clusters can be statistically analyzed in the same or adjacent azimuth groups to maintain the applicability of the scheme in sparse network or single-station research scenarios.

[0203] This invention is not limited to the above-described embodiments, and the following alternative methods may also be used:

[0204] 1. Time-frequency transformation can employ continuous wavelet transform, S-transform, short-time Fourier transform, synchronous compressed wavelet transform, or matched pursuit;

[0205] 2. Noise estimation can be performed using MAD, quantile estimation, local block variance estimation, robust root mean square estimation, or Bayesian noise estimation.

[0206] 3. Clustering methods can include four-neighbor connectivity, eight-neighbor connectivity, adaptive clustering based on connectivity constraints, density clustering, or graph connectivity clustering;

[0207] 4. The weighting coefficients or fusion index can be determined automatically through empirical setting, grid search, Bayesian optimization, multi-objective optimization, cross-validation, or a small number of labeled samples;

[0208] 5. If there are no available neighboring stations around the target station, the spatial consistency item can be replaced with the repetition of multiple events within the same station;

[0209] 6. The theoretical phase time window can be determined by a one-dimensional velocity model, a three-dimensional velocity model, a ray tracing method, or a data-driven time-of-arrival picking method;

[0210] 7. The receiver function can be calculated using frequency domain level factor deconvolution, time domain iterative deconvolution, multi-window spectrum deconvolution, or other deconvolution methods;

[0211] 8. Soft weights can be piecewise constant functions, continuous sigmoid functions, mapping functions based on cluster confidence scores, or weighting functions based on Bayesian posterior probabilities.

[0212] 9. Quality control indicators can be added, removed, or reweighted according to research objectives;

[0213] 10. The output results can be used for H-κ superposition, receiver function superposition, harmonic decomposition, CCP imaging, fixed network area imaging, or station quality assessment.

[0214] This invention does not perform independent filtering on a single receiver function. Instead, under a fixed station receiver function data structure, it incorporates the repeatability of multiple events at the same station, spatial consistency of neighboring stations, radial / tangential multi-component consistency, inverse azimuth grouping consistency, theoretical phase time window, azimuth harmonic protection, and quality control feedback into the time-frequency cluster screening and soft-weight reconstruction process.

[0215] Compared with traditional fixed-site receiver function processing methods, this invention has the following advantages: 1. Improves the signal-to-noise ratio of single-event receiver functions and grouped superimposed receiver functions at fixed stations; 2. Reduces the subjectivity caused by fixed bandpass filtering and manual thresholding; 3. Preserves Ps-converted phase, multiples, and deep weak interface responses; 4. Preserves anisotropy, tilted interface, and scattering information in tangential receiver functions; 5. Maintains true azimuth variation characteristics through inverse azimuth angle and ray parameter grouping; 6. Improves the stability of weak signal identification at target stations through spatial consistency of neighboring stations; 7. Maintains the applicability of the method by maintaining the repeatability of multiple events at the same station when there are no neighboring stations; 8. 9. Avoids the accidental deletion of weak azimuth periodic signals through azimuth harmonic protection mechanisms; 10. Reduces phase distortion and weak phase loss caused by hard thresholding through soft weight reconstruction; 11. Retains original complex phase information by weighting the original complex time-frequency coefficients; 12. Achieves noise attenuation rather than hard deletion through non-zero weighting of noise-dominant regions; 13. Outputs time-frequency cluster catalog, cluster credibility score, difference receiver function, and quality control indicators, making the noise reduction process auditable and interpretable; 14. Achieves adaptive adjustment of weight parameters, cluster screening threshold, and soft weight level through a quality control feedback mechanism; 15. Applicable to various seismological applications such as H-κ superposition, receiver function harmonic decomposition, CCP imaging, fixed network area imaging, and station quality assessment.

[0216] Although the present invention has been described in detail above with general descriptions and specific embodiments, modifications or improvements can be made to it, which will be obvious to those skilled in the art. Therefore, all such modifications or improvements made without departing from the spirit of the present invention fall within the scope of protection claimed by the present invention.

[0217] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Any simple modifications, equivalent changes, or alterations made by those skilled in the art using the disclosed technical content shall fall within the protection scope of the present invention.

Claims

1. A multi-constraint adaptive noise reduction method for a fixed station receiving function, characterized in that, The method includes: Acquire multi-event teleseismic three-component seismic records from the target fixed station and nearby stations; For the event seismic records corresponding to each event at each station, data preprocessing, coordinate rotation, and deconvolution are performed on the event seismic records to obtain the event reception function, which includes a radial reception function and a tangential reception function. The event receiving functions are grouped according to the reverse azimuth angle and ray parameters, and time-frequency transformation is performed on the event receiving functions in each event group to obtain the original complex time-frequency coefficients. Based on the original complex time-frequency coefficients, the frequency-related noise level is estimated using the P-wave arrival time-precedence noise window, and the local time-frequency signal-to-noise ratio is calculated based on the frequency-related noise level. The radial component score and tangential component score are calculated based on the local time-frequency signal-to-noise ratio. The radial component score and the tangential component score are used at least and fused with at least one of the following: radial / tangential local coherence, spatial consistency of neighboring stations, azimuth consistency and azimuth harmonic protection, to obtain a joint retention weight. Time-frequency connected component clustering is performed on the time-frequency units corresponding to the joint retention weights to divide them into noise-dominated clusters, effective clusters, and uncertain protection clusters; Based on the soft weights of the noise-dominant cluster, the effective cluster, and the uncertain protection cluster, the original complex time-frequency coefficients are subjected to inverse time-frequency transformation to obtain the noise-reduced receiver function.

2. The fixed station receiving function multi-constraint adaptive noise reduction method according to claim 1, characterized in that, Acquire multi-event teleseismic three-component seismic records from the target fixed station and neighboring stations, including: According to the target fixed station The target fixed station is determined by the preset spatial neighborhood radius. The neighboring stations, obtain the set of neighboring stations. ,in, This represents the surface distance between stations. The preset spatial neighborhood radius; Acquire multi-event teleseismic three-component seismic records from the target fixed station and neighboring stations.

3. The fixed station receiving function multi-constraint adaptive noise reduction method according to claim 2, characterized in that, The event receiving functions are grouped according to the inverse azimuth angle and ray parameters. A time-frequency transformation is then performed on the event receiving functions within each event group to obtain the original complex time-frequency coefficients, including: For each station For each event e, the inverse azimuth angle corresponding to event e is used. and ray parameters Group the event receiving functions to obtain event groups. ; The event group The expression is ,in, Center of the preset orientation window The preset orientation window half-width, Center the preset ray parameter window. This is the preset allowable range for ray parameters; For each event group The event receiving function within the function performs time-frequency transformation to obtain the original complex time-frequency coefficients. The expression for the original complex time-frequency coefficients is as follows: Where c represents the radial component or the tangential component, Indicates radial or tangential receiver function. This indicates a time-frequency transform operator, which is either a continuous wavelet transform operator or an S-transform operator; The original complex time-frequency coefficients include radial complex time-frequency coefficients and tangential complex time-frequency coefficients.

4. The fixed station receiving function multi-constraint adaptive noise reduction method according to claim 3, characterized in that, Based on the original complex time-frequency coefficients, the frequency-dependent noise level is estimated using the P-wave arrival time-precedence noise window, and the local time-frequency signal-to-noise ratio is calculated based on the frequency-dependent noise level, including: The time window preceding the arrival of the P-wave is used as the noise estimation window, and the expression for the noise estimation window is as follows: ,in, The time when the P-wave arrives at the station. and To preset the noise window parameters, ; Within the noise estimation window, the frequency-dependent noise level is estimated using the original complex time-frequency coefficients. The expression for the frequency-related noise level is: ,in, These are the original complex time-frequency coefficients. ; Calculate the local time-frequency signal-to-noise ratio using the frequency-dependent noise level. The calculation formula is: ,in, To prevent division by zero constant; The local time-frequency signal-to-noise ratio includes the radial local time-frequency signal-to-noise ratio and the tangential local time-frequency signal-to-noise ratio.

5. The fixed station receiving function multi-constraint adaptive noise reduction method according to claim 4, characterized in that, The radial component score and tangential component score are calculated based on the local time-frequency signal-to-noise ratio. At least the radial component score and the tangential component score are used, and at least one of the following is fused: radial / tangential local coherence, spatial consistency of neighboring stations, azimuth consistency, and azimuth harmonic protection, to obtain a joint retention weight, including: Based on the radial local time-frequency signal-to-noise ratio, the radial component score is calculated using the following formula: ,in, For radial local time-frequency signal-to-noise ratio, Represents the sigmoid function; Based on the tangential local time-frequency signal-to-noise ratio, the tangential component score is calculated using the following formula: ,in, For tangential local time-frequency signal-to-noise ratio; In the neighborhood of the preset time frequency Internal calculation of radial / tangential local coherence The calculation formula is: ,in, Represents the radial complex time-frequency coefficient. Represents the tangential complex time-frequency coefficients. Represents the complex conjugate of the tangential complex time-frequency coefficient; Based on the radial component score and the tangential component score, the spatial consistency of neighboring stations is calculated. The calculation formula is: ,in, , For indicator functions, Set a preset candidate retention threshold; azimuth consistency is calculated based on the radial component score and the tangential component score. The calculation formula is: ; The azimuth harmonic protection term is obtained by harmonic decomposition of the event receiving function as a function of the reverse azimuth angle. At least the radial component score and the tangential component score are used, and at least one of the following is fused: radial / tangential local coherence, spatial consistency of neighboring stations, azimuth consistency, and azimuth harmonic protection, to obtain the joint retention weight. .

6. The fixed station receiving function multi-constraint adaptive noise reduction method according to claim 5, characterized in that, The radial component score and tangential component score are calculated based on the local time-frequency signal-to-noise ratio. At least the radial component score and the tangential component score are used, and at least one of the following is fused: radial / tangential local coherence, spatial consistency of neighboring stations, azimuth consistency, and azimuth harmonic protection, to obtain a joint retention weight. The calculation also includes: If there are no neighboring stations within the preset spatial neighborhood radius of the target fixed station, or if the data quality of neighboring stations does not meet the preset conditions, then the repeatability of the same station event is calculated using the set of similar reverse azimuth and similar ray parameter events within the target fixed station. Instead of spatial consistency at neighboring stations, the expression for the repetition of events at the same station is: ,in, This represents the set of events that have similar reverse azimuth angles and similar ray parameters to event e; At least the radial component score and the tangential component score are used, and at least one of the following is fused: radial / tangential local coherence, repetitiveness of co-occurrence events, azimuth consistency, and azimuth harmonic protection terms, to obtain the joint retention weight. The fusion includes linear fusion, product-type fusion, or logistic fusion.

7. The fixed station receiving function multi-constraint adaptive noise reduction method according to claim 6, characterized in that, Time-frequency connected component clustering is performed on the time-frequency units corresponding to the joint retention weights to divide them into noise-dominated clusters, effective clusters, and uncertain protection clusters, including: Based on the joint retention weights, a candidate binary mask is constructed according to a preset clustering threshold, and the set of candidate time-frequency units to be discriminated is determined according to the candidate binary mask; Based on the candidate binary mask, candidate time-frequency unit pairs with a preset adjacency relationship are connected to obtain candidate time-frequency clusters to be judged. The preset adjacency relationship is a four-neighborhood, an eight-neighborhood, a frequency-weighted neighborhood, a time-weighted neighborhood, or a graph-connected neighborhood. For each candidate time-frequency cluster q, the confidence score corresponding to the candidate time-frequency cluster is calculated based on the cluster energy, tangential / radial energy ratio, neighbor station repeatability, cluster orientation consistency, and theoretical phase time window matching degree. The formula for calculating the cluster energy is: ; The formula for calculating the tangential / radial energy ratio is: ; The formula for calculating the repeatability of neighboring stations is: ; The formula for calculating the cluster orientation consistency is as follows: ; The formula for calculating the theoretical phase time window matching degree is as follows: ,in, For the theoretical phase arrival time, The center time of the candidate time-frequency cluster, This represents the theoretical phase time window width. Determine whether the confidence score of the candidate time-frequency cluster is greater than or equal to a preset effective cluster threshold; If the confidence score of the candidate time-frequency cluster is greater than or equal to the preset effective cluster threshold, then the candidate time-frequency cluster is regarded as an effective cluster. If the confidence score of the candidate time-frequency cluster is less than the preset effective cluster threshold, then it is determined whether the candidate time-frequency cluster meets at least one uncertain protection condition. The uncertain protection conditions include preset harmonic protection conditions, preset theoretical phase window matching conditions, preset neighbor station repeatability conditions, and co-station multi-event repeatability conditions. If the candidate time-frequency cluster satisfies at least one uncertainty protection condition, then the candidate time-frequency cluster is designated as an uncertainty protection cluster. If the candidate time-frequency cluster does not meet any of the uncertain protection conditions, then the candidate time-frequency cluster is regarded as the noise-dominant cluster.

8. A multi-constraint adaptive noise reduction method for a fixed station receiving function according to claim 7, characterized in that, Based on the soft weights of the noise-dominant cluster, the effective cluster, and the uncertain protection cluster, an inverse time-frequency transform is performed on the original complex time-frequency coefficients to obtain the denoised receiver function, including: Constructing soft weight functions ,in, For effective cluster regions The corresponding soft weights, For uncertain protected cluster regions The corresponding soft weights, Noise-dominant cluster region The corresponding soft weights are as follows: the soft weight corresponding to the effective cluster region is greater than the soft weight corresponding to the uncertain protection cluster region, the soft weight corresponding to the uncertain protection cluster region is greater than the soft weight corresponding to the noise-dominant cluster region, and the soft weight corresponding to the noise-dominant cluster region is greater than or equal to zero. The original complex time-frequency coefficients are subjected to inverse time-frequency transformation using the soft weighting function to obtain the noise-reduced receiving function. The inverse time-frequency transformation expression is as follows: Among them, the noise-reduced receiver function This includes the denoised radial receiver function and the denoised tangential receiver function. These are the weighted original complex time-frequency coefficients. , It is an inverse time-frequency transform operator.

9. A multi-constraint adaptive noise reduction method for a fixed station receiving function according to claim 8, characterized in that, The method further includes: Using the noise-reduced receiving function and the original event receiving function Calculate the difference receiving function The calculation formula is: ; Based on the noise-reduced receiver function and the original event receiver function, quality control indicators are calculated. These quality control indicators include Ps time offset, amplitude retention rate, tangential component energy retention rate, harmonic term retention rate, signal-to-noise ratio improvement rate, time-frequency energy retention rate, and effective phase residue index of the difference receiver function. Determine whether the quality control indicators meet the preset quality threshold; If the quality control indicators do not meet the preset quality threshold, the joint retention weight parameter, preset clustering threshold, preset effective cluster threshold, azimuth harmonic protection threshold, theoretical phase time window width, or soft weight level will be adjusted accordingly.

10. A computer program product, characterized in that, The computer program product includes computer program instructions, which, when executed by a processor, implement the steps of the multi-constraint adaptive noise reduction method for fixed station receiving functions as described in any one of claims 1 to 9.

Images