Method for extracting spatially regularized iterative deconvolution receiver functions based on dense array

Abstract

The present application relates to geophysical exploration and deep seismology data processing field, disclose a kind of spatial regularization iterative deconvolution receiver function extraction method based on dense array, comprising the following steps: step S1: the radial component and vertical component of each station of seismic array are overlapped and sliding sampling, construct multiple subarray data blocks with spatial topological relationship;Step S2: for each subarray data block, in the iterative process of deconvolution, based on the algebraic superposition of cross-correlation function of all stations in the subarray, coherent peak searching is carried out, and the pulse position of the current iteration is determined;Step S3: according to the pulse position, the pulse amplitude of each station is extracted, and the pulse amplitude is updated by using spatial smoothing operator.The present application introduces coherent peak searching mechanism, and the cross-correlation function of all stations in the subarray is algebraically superimposed, so that the random noise specific to single station cannot form effective peak value, thereby eliminating the generation of false seismic phase from the root.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration and deep seismic data processing, and in particular to a method for extracting spatially regularized iterative deconvolution receiver functions based on dense arrays. Background Technology

[0002] Receiver function (RFF) is a core method for studying the velocity structure of the Earth's crust and upper mantle. Traditional RRF extraction methods (such as leveling pool deconvolution and iterative deconvolution) are mostly performed on a single station. However, single-station processing is often severely affected by near-surface noise and local scattered waves, leading to difficulties in phase identification. Specifically, traditional single-station deconvolution has the following limitations:

[0003] Signal-to-noise ratio bottleneck: Single-station processing relies entirely on the seismic record of that station. Due to strong scattering from near-surface sedimentary layers, instrument noise, and environmental fluctuations, P-wave tails are often mixed with a large amount of random noise, causing key phases such as Ps and PpPs to be submerged.

[0004] Phase ambiguity: Single-station deconvolution can easily produce "pseudo-phases." In the absence of spatial constraints, researchers find it difficult to distinguish whether a wave crest is a genuine reflection from an underground interface or a local noise interference unique to a single station.

[0005] Spatial discontinuity: Independent processing of adjacent stations leads to jumps in the calculated Moho surface or breaks in the stratigraphic interface, which fails to reflect the continuous physical characteristics of the underground structure.

[0006] Poor numerical stability: Traditional time-domain iteration is prone to divergence when dealing with high-frequency components and is extremely sensitive to initial parameters (such as Gaussian factor), making it difficult to achieve automated processing of large-scale arrays.

[0007] Although single-station iterative algorithms demonstrate good stability in single-station processing, they have not been effectively extended to joint array processing under sliding window conditions. In existing technologies, researchers typically attempt to directly apply single-station algorithms to each station in the array, then simply average or stitch the results together. However, this fails to address the problem of spatial consistency failure—the receiver functions of adjacent stations often exhibit lateral jumps, and pseudo-phases are difficult to identify due to a lack of spatial constraints. Therefore, how to deeply integrate the stability of frequency-domain iteration with the spatial coherence of the array, and utilize neighboring station information in real-time to constrain pulse extraction during the inversion process, thereby obtaining a spatially continuous and phase-clean receiver function profile, is a pressing technical problem in this field, and existing technologies have not provided corresponding technical insights. Summary of the Invention

[0008] To address the technical problems mentioned in the background section, this invention provides a spatially regularized iterative deconvolution receiver function extraction method based on dense arrays.

[0009] This invention is achieved using the following technical solution: a spatially regularized iterative deconvolution receiver function extraction method based on dense arrays, comprising the following steps:

[0010] Step S1: Overlap sliding sampling is performed on the radial and vertical components of each station in the seismic array to construct multiple sub-array data blocks with spatial topological relationships;

[0011] Step S2: For each sub-array data block, during the deconvolution iteration process, coherent peak finding is performed based on the algebraic superposition of the cross-correlation functions of all stations within the sub-array to determine the pulse position of the current iteration.

[0012] Step S3: Extract the pulse amplitude of each station based on the pulse position, and use the spatial smoothing operator to update the pulse amplitude under constraints;

[0013] Step S4: Perform nonlinear weighted synthesis of multiple receiver function results obtained from different sub-array data blocks of the same station, and output the final receiver function waveform.

[0014] Furthermore, in step S1, the window width and sliding step size of the overlapping sliding sampling are automatically adjusted according to the station spacing.

[0015] Furthermore, in step S2, the cross-correlation function is obtained by calculating the frequency domain cross-correlation between the current residual signal and the vertical component, and then transforming it back to the time domain after Gaussian filtering.

[0016] Furthermore, in step S2, the coherent peak finding is achieved by solving for the time delay value that maximizes the absolute value of the algebraic superposition of the cross-correlation functions of all stations within the subarray.

[0017] Furthermore, in step S2, at least one of the following discrimination conditions is used when determining the pulse position:

[0018] An energy ratio threshold is set, and a valid phase is determined only when the coherence peak is significantly higher than the average energy level of the background noise within the search window.

[0019] Verify the polarity consistency of the cross-correlation function of each station in the sub-array at the pulse position;

[0020] Determine whether the pulse position has iterative evolution stability during the recursive iteration process.

[0021] Furthermore, in step S3, the spatial smoothing operator is a smoothing operator in the form of a Laplace matrix, and the lateral consistency of the received function waveforms of adjacent stations is controlled by adjusting the smoothing factor λ.

[0022] Furthermore, the Laplace matrix is ​​a second-order difference matrix, and for stations within the subarray, the corrected pulse amplitude satisfies:

[0023]

[0024] in The pulse amplitude before smoothing. This represents the smoothed pulse amplitude.

[0025] Furthermore, in step S4, the weights of the nonlinear weighted synthesis are calculated based on the fitting degree of the deconvolution iteration of each sub-array data block. The higher the fitting degree, the greater the corresponding weight.

[0026] Furthermore, the formula for calculating the weight is as follows: ,in σ represents the fitting percentage in the m-th calculation, and σ is the decay constant.

[0027] Furthermore, all time alignment operations in steps S1 to S4 are implemented using the frequency domain shift theorem, and the B field of the SAC header file and the IZTYPE reference time type are forcibly corrected during the storage phase using the underlying SAC protocol.

[0028] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0029] This invention introduces a coherent peak-finding mechanism to algebraically superimpose the cross-correlation functions of all stations within the subarray, preventing the unique random noise of a single station from forming an effective peak, thereby eliminating the generation of pseudo-phases at the source.

[0030] This invention utilizes the Laplace smoothing operator to spatially constrain the pulse amplitude, forcing the amplitude of the receiver function of adjacent stations to not undergo drastic changes, so that the velocity discontinuity interface presents a clear and continuous band on the receiver function profile, which greatly facilitates subsequent migration stacking imaging.

[0031] The proposed solution ensures that even if the data quality of a station in the array is poor, its result will be completed and corrected by adjacent high-quality stations through coherent weights during the sliding window processing. At the same time, the nonlinear weighted synthesis based on the fit degree can effectively suppress the contribution of the noise-affected window.

[0032] This invention is based on a frequency domain iterative algorithm and combines the Fourier translation theorem for time axis alignment, which avoids the cumulative error caused by time domain convolution and ensures stable extraction of high-frequency details.

[0033] The algorithm parameters of this invention (such as window size and step size) can be automatically adjusted according to the spacing between the arrays, making it suitable for rapid and automated processing of large-scale dense arrays. Attached Figure Description

[0034] Figure 1 This is a crustal velocity model used to synthesize seismic maps in an embodiment of the present invention;

[0035] Figure 2 The following is a comparison of the synthesized receiver function extraction results in the embodiments of the present invention (signal-to-noise ratio including 10% noise). (a) Theoretical receiver function benchmark; (b) Extraction results of traditional deconvolution method; (c)–(e) Extraction results of the method in this paper when the regularization parameter λ is 0.1, 0.2 and 0.3, respectively; (f) Root mean square error (RMS) curves under different methods and parameter settings, reflecting the degree of deviation of the extraction results from the theoretical benchmark (a);

[0036] Figure 3 Comparison of synthesized receiver function extraction results in this embodiment of the invention (signal-to-noise ratio including 20% ​​noise). (a) Theoretical receiver function benchmark; (b) Extraction results of traditional deconvolution method; (c)–(e) Extraction results of the method in this paper when the regularization parameter λ is 0.1, 0.2 and 0.3, respectively; (f) Root mean square error (RMS) curves under different methods and parameter settings, reflecting the degree of deviation of the extraction results from the theoretical benchmark (a);

[0037] Figure 4 The following is a comparison of ablation tests of Laplace smoothing factor and sliding step size in the embodiments of the present invention: (a) Laplace smoothing with sliding window superposition; (b) Laplace smoothing without sliding window superposition; (c) No Laplace smoothing with sliding window superposition; (d) No Laplace smoothing and no sliding window superposition; (e) Comparison of root mean square error curves of the four parameter combinations.

[0038] Figure 5 This is a map showing the actual location and distribution of the observation array and the distribution of teleseismic data in this embodiment of the invention;

[0039] Figure 6 This is a comparison chart of the results of the receiving function extracted by the conventional method in the embodiments of the present invention and the method proposed in the present invention. Detailed Implementation

[0040] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be noted that the specific parameters used in this embodiment (such as sub-array window width, sliding step size, number of iterations, smoothing factor, etc.) are merely illustrative and do not constitute a limitation on the scope of protection of this invention.

[0041] I. Data Preparation and Matrix Construction

[0042] In this embodiment, the three-component SAC format records observed by the seismic array are first read. The read data undergoes preprocessing, specifically including: instrument response removal, bandpass filtering, and coordinate system rotation. After the above processing, the radial components of each station are obtained. and vertical components .

[0043] To facilitate subsequent sliding slice processing, this embodiment arranges the records of all stations row by row, thereby constructing a full-array data matrix. (Radial component matrix) and (Vertical component matrix). Simultaneously, a path list is established. This is used to store the original SAC file paths corresponding to each station, so that the results can be indexed and stored in subsequent steps.

[0044] II. Signal Model and Objective Function Definition

[0045] In this embodiment, for the first in the array Each station, its observed radial component With vertical component The convolutional model is satisfied, as shown in Equation 1:

[0046] (1)

[0047] The objective of this invention is to find a set of receiving functions within the sub-array range. This makes the global residual functional Minimum, as shown in Formula 2

[0048] (2)

[0049] in It is an introduced spatial smoothness penalty term. It is a constraint factor. This objective function lays the theoretical foundation for subsequent inversion.

[0050] III. Slicing of the sliding sub-array

[0051] In this step, the present invention introduces a sliding window mechanism to spatially slice the entire array data. It is worth noting that in this embodiment, the width of the sub-array window is set to... (That is, each sub-window contains 3 stations), the sliding step size is (That is, slide back one station at a time). Use a loop iterator to iterate over the radial component matrix constructed in step one. and vertical component matrix Perform row slicing operations.

[0052] Through the aforementioned sliding slicing, this embodiment generates a series of sub-array data blocks with spatial topological relationships, denoted as... subscript Indicates the first Each sub-array window. This approach implicitly encodes the spatial continuity of the array into the subsequent inversion process.

[0053] IV. Coherent Constraint Joint Deconvolution of Sub-arrays

[0054] This is the core step of the invention. For each sub-array data block obtained through sliding slicing... This embodiment performs the following calculation process:

[0055] (1) Global coherence peak finding

[0056] First, it should be noted that this invention abandons the traditional method of finding peaks independently station by station, and instead defines a coherent cumulative energy function for the sub-array. In the iteration of the _____ Next, we need to find the optimal position for the next pulse. .

[0057] Specifically, this embodiment first calculates the residual signal. With vertical component The frequency domain cross-correlation. Let... The Fourier transform is calculated in the frequency domain according to the following formula 3;

[0058] (3)

[0059] in For a Gaussian filter, its expression is shown in Equation 4 below, which is used to control the bandwidth and ensure the numerical stability of the iteration process;

[0060] (4)

[0061] It is important to emphasize that the coherent peak finding here uses frequency domain cross-correlation superposition of the "residual-vertical component," rather than directly superimposing cross-correlation on the original radial components. This is because directly superimposing cross-correlation on the radial components is subject to severe noise interference. Within the framework of iterative deconvolution, each round deals with the residual signal, which represents uninterpreted energy. By cross-correlating the residual with the vertical component, we are essentially performing an "online search of deconvolution." Only those residual energies that match the source waveform (vertical component) are physically reliable converted wave phases, which is much more accurate than "blindly searching" for energy peaks directly on the radial component.

[0062] Subsequently, the frequency domain results were transformed back to the time domain to obtain the cross-correlation functions of each station. To extract the common phase within the subarray, this embodiment constructs a criterion function as shown in Formula 5 below:

[0063] (5)

[0064] This formula uses algebraic superposition (rather than absolute value or square superposition) to strengthen spatially coherent signals and cancel out randomly distributed noise.

[0065] Furthermore, the reason for using algebraic superposition instead of absolute value or square superposition is that the biggest advantage of algebraic superposition is that it utilizes polarity constraints. The actual seismic phases are not only close in position within the subarray, but their waveform polarities (positive and negative directions) should also be consistent.

[0066] Algebraic superposition enhances the energy of the same polarity while canceling out noise of opposite polarities, which is a form of "spatial polarity filtering". While absolute or square superposition can increase energy, it loses polarity information, causing noise of opposite phase to be added in, which can easily produce spurious peaks.

[0067] Based on this, in order to transform the above assumptions into an executable discrimination mechanism, this embodiment determines... The following three specific constraints were also introduced:

[0068] Significance determination: Set an energy ratio threshold. Only when the coherent peak value Amax after superposition is significantly higher than the average energy level of the background noise within the search window is it determined to be a valid phase.

[0069] Polarity consistency verification: determined by coherent peak finding At this point, the algorithm will further verify the polarity (positive or negative) of the original cross-correlation function of each station within the sub-array. True interface reflections should have a consistent polarity within the local array.

[0070] Iterative evolution stability: Because the algorithm uses recursive iteration, the actual seismic phases usually exhibit stable peak positions in the initial few iterations. If an energy peak appears only in specific iterations and cannot maintain spatial stability in subsequent residual updates, it will be considered an incoherent random error.

[0071] In summary, this is one of the advantages of this algorithm: traditional single-station-based extraction algorithms cannot distinguish between real signals and noise, while this invention can effectively distinguish between the two through the above-mentioned triple discrimination mechanism.

[0072] (2) Pulse amplitude extraction and spatial Laplace smoothing

[0073] Determining the pulse position Then, in this embodiment, the preliminary pulse intensity corresponding to each station is calculated according to Formula 6;

[0074] (6)

[0075] To ensure the spatial continuity of the receiver function profile and suppress singular amplitude values ​​caused by single-station site effects, this invention introduces a spatial Laplace smoothing operator. Let the intensity vector be as shown in Equation 7:

[0076] (7)

[0077] This embodiment calculates the smoothed intensity vector based on the following formula 8. ;

[0078] (8)

[0079] in It is a second-order difference matrix (Laplace matrix), and its form is shown in Equation 9:

[0080] (9)

[0081] After derivation, for stations within the subarray, the corrected amplitude satisfies the equation shown in Formula 10:

[0082] (10)

[0083] It is important to note that Laplace smoothing is applied to the pulse signals extracted in each iteration, not to the overall waveform of the final receiver function. Smoothing the final waveform after inversion (i.e., "post-processing smoothing") is essentially performing post-processing low-pass filtering, which inevitably smooths out the sharpness of the seismic phases along the time axis, reducing vertical resolution. The logic employed in this invention is to introduce regularization, or "intra-iterative constraints," during the model building phase. The physical significance of this approach is that it only constrains the continuity of the seismic phases in the transverse direction (dynamic characteristics), while fully preserving their instantaneous nature along the time axis (kinematic characteristics).

[0084] Furthermore, this "iterative internal constraint" is more valuable than "post-event smoothing" because it creates a dynamic recursive feedback: the pulses injected into the model in each round are spatially coherently corrected, resulting in cleaner residual signals that guide subsequent iterations to more accurately capture weak seismic phases. This method effectively suppresses local site effect interference while clearly preserving abrupt structural details such as the Moho step, which is the core advantage of this algorithm.

[0085] (3) Derivation of recursive update and time translation

[0086] After obtaining the pulse components of this round, construct the pulse function according to the following formula 11.

[0087] (11)

[0088] Subsequently, the receiver function estimate is updated according to Formula 12 below, and the residual signal is updated synchronously.

[0089] (12)

[0090] In this embodiment, the maximum number of iterations is set to... Next. In practical applications, this parameter can be selected based on the convergence situation. Repeat the above process (1) to (3) until the sum of squared residuals no longer decreases significantly, at which point the receiving function waveform matrix of the sub-array is output.

[0091] For ease of final imaging, the results need to be shifted over time. According to the Fourier translation theorem, as shown in Equation 13 below, transforming back to the time domain yields the final output waveform:

[0092] (13)

[0093] V. Result Registration and Overlapping Storage

[0094] This embodiment establishes a result registry. It iterates through the output results of all sub-arrays in step four. For each station in each sub-array, its original storage path is used as an index, and the waveform data, fit (Fit%), ray parameters, etc. obtained in this calculation are stored as a "record" in the list corresponding to that path.

[0095] A key advantage of this sliding window overlay process is that it allows the window to identify step changes in discontinuities. When the window crosses a tectonic anomaly zone, the in-phase behavior of stations within the window changes, and the response fluctuations of the same station in different contexts are fully recorded, providing a rich information basis for subsequent weighted synthesis.

[0096] VI. Multi-window weighted synthesis

[0097] After all window calculations are completed, the aggregation phase begins. This applies to stations with multiple calculation records in the repository. This embodiment abandons the simple arithmetic mean and instead uses a nonlinear weighted exponential algorithm based on the degree of fit for synthesis.

[0098] First, define the first The fit weights calculated in this step The calculation method is shown in Formula 14 below;

[0099] (14)

[0100] in This represents the percentage of fit obtained from the inversion of the secondary array.

[0101] Furthermore, in the actual programming implementation, this embodiment adopts the simplified form shown in Formula 15 below.

[0102] w=exp((Fit-100) / 10)(15)

[0103] The weighting function has the following explicit physical logic:

[0104] Monotonically increasing property: due to the percentage of fit The theoretical upper limit is 100, and the numerator term Always negative or zero. Exponential function. Monotonically increasing, as As the exponent increases and approaches 100, it will increase from negative infinity and approach 0, thus increasing the weight. exist The coefficient of fit increases monotonically within the interval. Therefore, the higher the coefficient of fit, the greater the corresponding weight.

[0105] Considerations for numerical stability: Using Instead of using directly It is a normalized anchor point that uses perfect fit (100%) as the weight (i.e. This effectively avoids the risk of numerical overflow that may occur when directly performing exponential operations on large positive numbers.

[0106] Nonlinear suppression mechanism: The exponential operator is introduced to construct a "soft threshold" screening mechanism. This is achieved by adjusting the decay constant. This allows for flexible control over the sensitivity of weights to fitting errors. In practical inversion, this nonlinear mapping can significantly suppress interference from sudden noise or non-convergent windows (low latency). The contribution of the final receiver function profile.

[0107] Finally, the high signal-to-noise ratio receiving function waveform of the station is synthesized according to the following formula 16.

[0108] (16)

[0109] By using the above weighted superposition, the boundary effect of the sub-array is effectively eliminated, and a unique set of high signal-to-noise ratio receiver function results for the entire array is obtained.

[0110] VII. Absolute Time Alignment and High-Fidelity Storage

[0111] Finally, this embodiment calls the underlying SACTrace library, based on the preset delay parameters. Force correction of the B field in the SAC header file Then, set IZTYPE to 'IB' (reference time type). After completing the above processing, write the synthesized data to disk, thus completing the entire process of extracting the receiving function of a dense array based on coherence constraints.

[0112] VIII. Performance Verification and Comparative Analysis

[0113] To evaluate the reliability and calculation accuracy of the present invention, this embodiment conducted synthetic seismic map testing and applied actual observation data.

[0114] (1) Simulation calculation test

[0115] This embodiment first obtains theoretical waveform records through forward modeling. Figure 1 A two-dimensional velocity model used for testing is presented. This model incorporates a laterally non-uniform Mohorovičić discontinuity (Moho), with a crustal thickness of 30 km on the left and increasing to 35 km on the right. The crustal medium is assumed to be isotropic, with a P-wave velocity (Vp) of 6.3 km / s and a S-wave velocity (Vs) of 3.6 km / s. A linear array of 20 seismic stations is uniformly distributed across the surface.

[0116] To evaluate the robustness of the algorithm in noisy environments, this embodiment adds 10% and 20% Gaussian random noise to the synthesized waveform, respectively. Figure 2 The comparison of extraction results at a 10% noise level is shown. Figure 2 (a) is the theoretical receiver function benchmark. Figure 2 (b) shows the results extracted using the traditional deconvolution method. Figure 2 (c)-(e) show the extraction results of the proposed method when the regularization parameter λ is 0.1, 0.2, and 0.3, respectively. Figure 2 (f) shows the root mean square error (RMS) curves under different methods and parameter settings. The results show that as the Laplace smoothing factor λ gradually increases from 0.1 to 0.3, the extracted receiver function waveform exhibits better stability and a significantly improved agreement with the theoretical benchmark.

[0117] To quantitatively evaluate the accuracy and robustness of the algorithm, this scheme calculates the root mean square error (RMS) between the extraction results of different methods and the theoretical receiver function, to characterize the degree of deviation between the two (e.g., ...). Figure 2(as shown in f). Analysis results show that the RMS error of the traditional method is significantly higher than that of the proposed method, reflecting its limitation in waveform reconstruction under noise interference. In contrast, the proposed method effectively reduces the error level by introducing spatial constraints. Furthermore, as the value of λ increases, the RMS error shows a monotonically decreasing trend, further demonstrating the crucial role of the Laplace smoothing operator in suppressing incoherent noise and improving the reliability of receiver function extraction.

[0118] To further verify the robustness of the algorithm under extreme noise conditions, we increased the noise level to 20% to simulate a low signal-to-noise ratio observation environment (such as...). Figure 3 As shown in the figure. Experimental results show that under high noise interference, the quality of the receiver function extracted by traditional methods deteriorates significantly: although the Moho conversion wave (Ps phase) is still barely discernible, multiple waves carrying deep structural information (such as PpPs and PsPs+PpSs phases) are completely submerged in background noise and difficult to identify. In contrast, the method presented in this paper exhibits superior noise reduction performance and phase recovery capability. As the smoothing factor λ increases to 0.3, not only is the signal-to-noise ratio of the Ps phase significantly improved, but the previously blurred multiple reflection phases also become clearly visible, providing a reliable data foundation for subsequent H-kappa stacking or crustal modeling. Furthermore, Figure 3 The trend shown by f is that the RMS error decreases as λ increases. Figure 2 Consistent and quantitative evidence confirms that the algorithm still exhibits good convergence and stability even in the context of strong noise.

[0119] To objectively evaluate the technical gains of each core functional module in the algorithm, this scheme designs an ablation test targeting the Laplacian smoothing term (λ) and the sliding window overlap mechanism. Using the controlled variable method, the inversion effects and the evolution characteristics of the root mean square error (RMS) under four different parameter combinations are compared. The results are as follows: Figure 4 As shown.

[0120] Analysis of ablation test results: From the perspective of the lateral continuity and convergence stability of the inverted waveform, the performance of the complete algorithm ( Figure 4 a): When Laplace smoothing and sliding window overlap are introduced simultaneously, the extracted receiver function profile has optimal lateral continuity, the Moho phase is clear and stable, and the pseudo-phases caused by random noise are suppressed to the greatest extent. Figure 4 The red curve in (e) shows that the root mean square error (RMS) converges the fastest and the final residual level is the lowest under this combination.

[0121] The contribution of the overlapping and superposition mechanism ( Figure 4b): When maintaining Laplace smoothing but eliminating sliding window overlap (step size equals window length), the stability of the receiver function decreases significantly. The waveform begins to exhibit subtle sawtooth-like jumps in space, and the corresponding RMS curve shows significant fluctuations. This indicates that multi-window nonlinear weighted superposition plays a crucial role in smoothing sub-array edge differences and eliminating spatial sampling distortion.

[0122] The core role of Laplace smoothing ( Figure 4 c): If the sliding window stacking is retained but Laplace smoothing is removed (lambda=0), the overall profile effect deteriorates further, with a large amount of incoherent anomalous energy being incorrectly extracted. A comparison shows that the performance loss is greater than that in Figure (b), which strongly demonstrates that the Laplace smoothing term contributes better to constraining the diffusion of lateral anomalous energy than the stacking mechanism, and is the core operator ensuring the stability of the inversion dynamics.

[0123] Basic performance benchmark ( Figure 4 d): When both are cancelled, the algorithm degenerates into a simple sliding inversion involving only coherent peak finding. In this case, the profile is cluttered, the RMS level is at its highest, and reliable geological interpretation cannot be obtained against a strong noise background.

[0124] In summary, both Laplace smoothing and sliding window weighted stacking significantly improve inversion accuracy. Laplace smoothing lays the foundation for spatial stability, while multi-window weighted stacking provides significant gains in the transition between local details and overall continuity. The synergistic effect of these two methods constitutes the core of the algorithm's high robustness.

[0125] Regarding the impact of the smoothing factor λ on the imaging results, this scheme establishes a protection mechanism of 'travel time independence and amplitude constraint' to prevent lateral abrupt changes in geological structures from being erroneously smoothed.

[0126] The smoothing term acts only on the impulse amplitude vector to suppress incoherent amplitude disturbances caused by local site responses. Interface undulations or fault jumps deep within the Earth's crust are primarily carried out by the kinematic travel time term τ. Because the algorithm allows for local travel time deviations during the peak-finding phase, the true tectonic jump information is absorbed into the travel time parameter, thus preserving physical fidelity during spatial smoothing.

[0127] Furthermore, the application of smoothing intensity is not indiscriminate coverage, but rather linked to the coherence of adjacent stations: in gentle interface regions with high waveform consistency, λ exerts a strong constraint effect to improve the signal-to-noise ratio; in tectonic abrupt change zones with severe waveform distortion, the spatial constraint automatically weakens as coherence decreases. This adaptive adjustment strategy ensures that the algorithm can accurately characterize fine geological features such as the Moho discontinuity steps while suppressing anomalous energy.

[0128] Experimental example:

[0129] Reference Figures 5-6 To further verify the effectiveness of our proposed method in processing real-world complex source data, we applied it to measured data from the Zhangbaling uplift area in the southern segment of the Tanlu Fault Zone. The experimental data relied on a linear array of 102 short-period seismographs deployed in this region, such as... Figure 5 The array spans the core region of the fault zone, with an average inter-station spacing of approximately 2 km, and continuous observation lasted for nearly a month. During the data preparation phase, we rigorously screened the teleseismic events recorded by the array, selecting high signal-to-noise ratio events with epicentral distances between 30° and 95° and magnitudes greater than 5.8. The inset in the upper left corner shows the distribution of teleseismic receiver functions; red indicates teleseismic events used in the calculations, while blue indicates actual events not used, mainly due to smaller magnitudes or lower signal-to-noise ratios.

[0130] Subsequently, we extracted the receiver function sequence using both the traditional deconvolution method and the algorithm presented in this paper. A comparison of their actual imaging results is shown below. Figure 6 Analysis results show that while the receiver functions obtained by traditional methods can identify obvious Ps converted waves at some stations, the lateral continuity of the phases is poor across the entire array profile, and the waveforms exhibit significant incoherent noise interference between stations. In contrast, the proposed method, by introducing spatial constraints, significantly enhances the consistency of receiver functions between adjacent stations, resulting in excellent spatial coherence of the Ps phases. Especially in the area of ​​stations 70-102, due to the strong multiple waves and near-surface scattering caused by deep sedimentary layers, the waveforms obtained by traditional methods appear disordered and the phases blurred; while the proposed method demonstrates superior low signal-to-noise ratio processing capabilities, effectively suppressing incoherent background noise and significantly improving the clarity of the receiver functions and imaging quality. This proves that the new algorithm has stronger robustness and more refined structural characterization capabilities when processing complex crustal structure data.

[0131] In summary, this embodiment fully demonstrates that the coherence-constrained dense array receiver function extraction method proposed in this invention has achieved significantly better results than existing technologies in terms of improving signal-to-noise ratio, ensuring spatial continuity, enhancing algorithm robustness, and realizing automated processing.

[0132] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Those skilled in the art can make various improvements and modifications without departing from the spirit and principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A spatially regularized iterative deconvolution receiver function extraction method based on dense arrays, characterized in that, Includes the following steps: Step S1: Overlap sliding sampling is performed on the radial and vertical components of each station in the seismic array to construct multiple sub-array data blocks with spatial topological relationships; Step S2: For each sub-array data block, during the deconvolution iteration process, coherent peak finding is performed based on the algebraic superposition of the cross-correlation functions of all stations within the sub-array to determine the pulse position of the current iteration. Step S3: Extract the pulse amplitude of each station based on the pulse position, and use the spatial smoothing operator to update the pulse amplitude under constraints; Step S4: Perform nonlinear weighted synthesis of multiple receiver function results obtained from different sub-array data blocks of the same station, and output the final receiver function waveform.

2. The method according to claim 1, characterized in that, In step S1, the window width and sliding step size of the overlapping sliding sampling are automatically adjusted according to the station spacing.

3. The method according to claim 1, characterized in that, In step S2, the cross-correlation function is obtained by calculating the frequency domain cross-correlation between the current residual signal and the vertical component, and then transforming it back to the time domain after Gaussian filtering.

4. The method according to claim 3, characterized in that, In step S2, the coherent peak finding is achieved by solving for the time delay value that maximizes the absolute value of the algebraic superposition of the cross-correlation functions of all stations in the subarray.

5. The method according to claim 1, characterized in that, In step S2, when determining the pulse position, at least one of the following discrimination conditions is also used: An energy ratio threshold is set, and a valid phase is determined only when the coherence peak is significantly higher than the average energy level of the background noise within the search window. Verify the polarity consistency of the cross-correlation function of each station in the sub-array at the pulse position; Determine whether the pulse position has iterative evolution stability during the recursive iteration process.

6. The method according to claim 1, characterized in that, In step S3, the spatial smoothing operator is a smoothing operator in the form of a Laplace matrix, and the lateral consistency of the received function waveforms of adjacent stations is controlled by adjusting the smoothing factor λ.

7. The method according to claim 6, characterized in that, The Laplace matrix is ​​a second-order difference matrix. For stations within the subarray, the corrected pulse amplitude satisfies: ； in The pulse amplitude before smoothing. This represents the smoothed pulse amplitude.

8. The method according to claim 1, characterized in that, In step S4, the weights of the nonlinear weighted synthesis are calculated based on the fitting degree of the deconvolution iteration of each sub-array data block. The higher the fitting degree, the greater the corresponding weight.

9. The method according to claim 8, characterized in that, The formula for calculating the weight is: ,in σ represents the fitting percentage in the m-th calculation, and σ is the decay constant.

10. The method according to claim 1, characterized in that, All time alignment operations in steps S1 to S4 are implemented using the frequency domain shift theorem, and the B field of the SAC header file and the IZTYPE reference time type are forcibly corrected during the storage phase using the underlying SAC protocol.

Images