An ice cover snow grain density estimation method and system based on a vehicle-mounted noise source

Abstract

The application provides an ice cover granular snow layer density estimation method and system based on a vehicle-mounted noise source, relates to the technical field of physical exploration and ice cover detection, and the method comprises the following steps: according to the exact start and end time period of each vehicle-mounted noise event, a corresponding multi-channel vehicle-mounted noise event data is obtained by intercepting from an original continuous noise record; the mutual coherence processing and linear superposition of the intercepted multi-channel vehicle-mounted noise event data are carried out to obtain a virtual shot record with improved signal-to-noise ratio; P wave and surface wave in the record are jointly inverted to obtain a P wave and S wave velocity structure, and on the basis, a double-parameter granular snow layer density formula fitting is carried out. The application is suitable for realizing the quantitative estimation of the regional scale ice cover granular snow layer density structure by using passive source seismic signals, and provides parameter support for ice cover engineering activity site evaluation and environmental evolution research.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration and ice sheet detection technology, and in particular to a method and system for estimating the density of ice sheet granular snow layer based on a vehicle-mounted noise source. Background Technology

[0002] The density of the granular snow layer on ice sheets is a core parameter reflecting their structural evolution and snow-ice transformation process. It directly determines the contribution of the granular snow column to the change in the surface elevation of the ice sheet, and at the same time provides key basic data support for the assessment of engineering activity sites and environmental change research in some polar regions.

[0003] Currently, most technologies for obtaining the density of ice sheet granular snow layers, especially in the crucial step of converting seismic wave velocity information into density, face the following technical shortcomings: existing empirical relationships between velocity and density mostly rely on a single wave velocity constraint, failing to simultaneously incorporate the propagation characteristics of both P-wave velocity (Vp) and S-wave velocity (Vs) for consistent constraint.

[0004] The above-mentioned shortcomings mean that a single wave velocity constraint cannot fully integrate the physical information of the subsurface medium carried by P-waves and S-waves respectively. P-wave velocity mainly reflects the volume compression characteristics of the medium, while S-wave velocity mainly reflects the shear characteristics of the medium. Empirical relationships constructed based on only one wave velocity are insufficient to comprehensively characterize the densification physical process of granular snow layers, resulting in insufficient accuracy and robustness of density estimation results. In scenarios where only a small amount of core data provides constraints, the lack of empirical relationships that constrain the consistent propagation characteristics of P-waves and S-waves makes it difficult to reliably extend the local density information of sparse cores to the regional scale. This leads to a technical bottleneck in the quantitative estimation of the density structure of granular snow layers at the regional scale, failing to meet the needs of engineering assessments and environmental studies for large-scale density distribution data. Summary of the Invention

[0005] The technical problem to be solved by the present invention is to provide a method and system for estimating the density of ice sheet granular snow layer based on vehicle-mounted noise sources. This method is applicable to the quantitative estimation of the density structure of ice sheet granular snow layer at the regional scale using passive source seismic signals, and provides parameter support for site assessment and environmental evolution research of ice sheet engineering activities.

[0006] To solve the above-mentioned technical problems, the technical solution of the present invention is as follows:

[0007] In a first aspect, a method for estimating the density of ice sheet granular snow layer based on vehicle-mounted noise sources, the method comprising:

[0008] Vehicle noise events are identified and extracted from the original continuous noise records. The energy ratio of the short-term average to the long-term average is calculated and compared with a preset threshold. Finally, the exact start and end time periods of each vehicle noise event are output.

[0009] Based on the exact start and end time periods of each vehicle noise event, corresponding multi-channel vehicle noise event data are extracted from the original continuous noise record; the extracted multi-channel vehicle noise event data are subjected to mutual interference processing and linear superposition to obtain a virtual gun set record with improved signal-to-noise ratio.

[0010] The dispersion energy of Rayleigh surface wave and P-guided wave is extracted from the virtual shot gather record after the signal-to-noise ratio is improved, and the corresponding Rayleigh surface wave dispersion curve and P-guided wave dispersion curve are generated by dispersion analysis.

[0011] The Rayleigh surface wave dispersion curve and the P-guided wave dispersion curve are jointly inverted. By constructing and iteratively optimizing the objective function to update the velocity model synchronously, a convergent one-dimensional velocity structure model of the P-wave velocity and the S-wave velocity is finally obtained.

[0012] By combining the measured density data of known borehole ice cores at multiple depth points with the P-wave velocity and S-wave velocity obtained by inversion at the corresponding depth, a power-law type empirical formula for density with P-wave velocity and S-wave velocity as independent variables is established.

[0013] Substituting the P-wave velocity and S-wave velocity values ​​at each depth in the one-dimensional velocity structure model into the established power-law density empirical formula, the density is calculated layer by layer, and finally the density distribution of the ice sheet granular snow layer at the regional scale is output.

[0014] Furthermore, vehicle noise events are identified and extracted from the original continuous noise records. The energy ratio of the short-term average to the long-term average is calculated and compared with a preset threshold. Finally, the exact start and end time periods of each vehicle noise event are output, including:

[0015] Obtain the original continuous noise record containing background noise and vehicle noise;

[0016] Calculate the average energy of the short-time window and the average energy of the long-time window of the original continuous noise record;

[0017] An energy ratio sequence is generated based on the calculated ratio of the average energy of the short time window to the average energy of the long time window.

[0018] The energy ratio sequence is compared with a preset trigger threshold to determine the start time of the vehicle noise event;

[0019] After determining the start time of the vehicle noise event, the system continuously monitors the time when the energy ratio sequence drops to a preset detriggered threshold to determine the end time of the vehicle noise event.

[0020] Based on the determined start and end times of the vehicle noise events, the exact start and end time periods of each vehicle noise event are output.

[0021] Furthermore, based on the exact start and end time periods of each vehicle-mounted noise event, corresponding multi-channel vehicle-mounted noise event data are extracted from the original continuous noise record; the extracted multi-channel vehicle-mounted noise event data undergoes mutual interference processing and linear superposition to obtain a virtual shot set record with improved signal-to-noise ratio, including:

[0022] Based on the exact start and end time periods of each output vehicle noise event, the corresponding data segments are extracted from the original continuous noise records to obtain multichannel vehicle noise event data.

[0023] Multichannel vehicle noise event data are subjected to mutual interference processing to obtain the corresponding mutual interference results;

[0024] The mutual interference results of multiple vehicle noise events are linearly superimposed to suppress random noise components, resulting in the superimposed mutual interference results.

[0025] Based on the superimposed mutual coherence results, a virtual shot gather record with clear surface wave and guided wave signals is reconstructed.

[0026] Furthermore, the dispersion energy of Rayleigh surface waves and P-guided waves is extracted from the virtual shot gather record after the signal-to-noise ratio improvement, and the corresponding Rayleigh surface wave dispersion curves and P-guided wave dispersion curves are generated through dispersion analysis, including:

[0027] The virtual shot gather record is preprocessed to enhance the Rayleigh surface wave and P-guided wave signals contained therein, and the preprocessed virtual shot gather record is obtained.

[0028] Rayleigh surface wave signal and P-guided wave signal were identified and extracted from the preprocessed virtual shot gather record;

[0029] Dispersion analysis was performed on the extracted Rayleigh surface wave signal to obtain the dispersion energy distribution of the Rayleigh surface wave;

[0030] The extracted P-guided wave signal was subjected to dispersion analysis to obtain the dispersion energy distribution of the P-guided wave.

[0031] Based on the dispersion energy distribution of Rayleigh surface waves, the corresponding Rayleigh surface wave dispersion curves are generated.

[0032] Based on the dispersion energy distribution of the P-guided wave, the corresponding P-guided wave dispersion curve is generated.

[0033] Furthermore, the Rayleigh surface wave dispersion curve and the P-guided wave dispersion curve are jointly inverted. By constructing and iteratively optimizing the objective function to synchronously update the velocity model, a convergent one-dimensional velocity structure model of the P-wave velocity and S-wave velocity is finally obtained, including:

[0034] Based on the generated Rayleigh surface wave dispersion curve and the generated P-guided wave dispersion curve, an initial one-dimensional velocity model including P-wave velocity and S-wave velocity parameters is set.

[0035] Construct a joint inversion objective function that includes a data fitting term and a model regularization term; wherein, the data fitting term is based on the difference between the dispersion curve obtained by forward modeling of the initial one-dimensional velocity model and the observed dispersion curve;

[0036] The joint inversion objective function is optimized and iterated using a weighted least squares iterative algorithm, and the P-wave velocity and S-wave velocity parameters in the one-dimensional velocity model are updated synchronously.

[0037] After each optimization iteration, calculate the dispersion curve fitting error corresponding to the current model and determine whether the preset convergence condition has been met.

[0038] When the preset convergence condition is met, the iteration stops and the finally converged one-dimensional velocity structure model of the P-wave velocity and S-wave velocity is output.

[0039] Furthermore, by combining measured density data from known borehole ice cores at multiple depths with the P-wave and S-wave velocities obtained from inversion at the corresponding depths, a power-law empirical formula for density with P-wave and S-wave velocities as independent variables is established, including:

[0040] Obtain the P-wave velocity and S-wave velocity at each depth in the finally converged one-dimensional velocity structure model, as well as the measured density data of the known borehole ice core at the corresponding depth points;

[0041] Establish a correspondence between the P-wave velocity and S-wave velocity at the same depth point and the measured density data to form a data sample set for fitting.

[0042] A power-law model is constructed with P-wave velocity and S-wave velocity as independent variables and density as the dependent variable.

[0043] Based on the data sample set, a regularized fitting method is used to determine the undetermined parameters in the power law relationship model, and a power law relationship model with determined parameters is generated.

[0044] The output parameter-defined power-law relationship model serves as a power-law empirical formula for density estimation.

[0045] Furthermore, the P-wave velocity and S-wave velocity values ​​at each depth in the one-dimensional velocity structure model are substituted into the established power-law density empirical formula to calculate the density layer by layer. The final output is the density distribution of the ice sheet granular snow layer at the regional scale, including:

[0046] For each depth point in the one-dimensional velocity structure model, read the corresponding longitudinal wave velocity value and transverse wave velocity value at that depth point.

[0047] Substitute the P-wave velocity and S-wave velocity values ​​of each depth point into the power-law density empirical formula to calculate the density estimate of the depth point.

[0048] The depth information of all depth points and their corresponding density estimates are collected to form a distribution data set that characterizes the density variation with depth;

[0049] Output a set of distribution data representing the density variation with depth, as the density distribution of ice sheet granular snow layer with depth at the regional scale.

[0050] Secondly, an ice sheet snow layer density estimation system based on vehicle-mounted noise sources includes:

[0051] The noise event extraction and processing module is used to identify and extract vehicle noise events from the original continuous noise records. By calculating the energy ratio of the short-term window average to the long-term window average and comparing it with a preset threshold, the module finally outputs the exact start and end time periods of each vehicle noise event.

[0052] The virtual gun set reconstruction module is used to extract the corresponding multi-channel vehicle noise event data from the original continuous noise record based on the exact start and end time period of each vehicle noise event; the extracted multi-channel vehicle noise event data is subjected to mutual interference processing and linear superposition to obtain the virtual gun set record with improved signal-to-noise ratio.

[0053] The multi-wave dispersion curve extraction module is used to extract the dispersion energy of Rayleigh surface wave and P-guided wave from the virtual shot gather record after the signal-to-noise ratio is improved, and generate the corresponding Rayleigh surface wave dispersion curve and P-guided wave dispersion curve through dispersion analysis.

[0054] The P-wave and S-wave velocity joint inversion module is used to jointly invert the Rayleigh surface wave dispersion curve and the P-guided wave dispersion curve. By constructing and iteratively optimizing the objective function, the velocity model is updated synchronously, and finally a converged one-dimensional velocity structure model of P-wave velocity and S-wave velocity is obtained.

[0055] The empirical formula construction module is used to establish a power-law type empirical formula for density with P-wave velocity and S-wave velocity as independent variables by combining the measured density data of known borehole ice cores at multiple depth points with the P-wave velocity and S-wave velocity obtained by inversion at the corresponding depth.

[0056] The density calculation and output module is used to substitute the P-wave velocity and S-wave velocity values ​​at each depth in the one-dimensional velocity structure model into the established power-law density empirical formula to calculate the density layer by layer, and finally output the density distribution of the ice sheet granular snow layer at the regional scale.

[0057] The above-described solution of the present invention has at least the following beneficial effects:

[0058] Because it employs vehicle-mounted noise event extraction and virtual shot set reconstruction, Rayleigh surface wave and P-guided wave joint inversion, and regularized power-law fitting techniques for P-wave and S-wave velocities and densities, it overcomes the technical problems of insufficient estimation accuracy and robustness caused by the reliance on a single wave velocity constraint in existing velocity-density empirical relationships, and the difficulty in extending them to regional scales. This achieves the technical effects of improving the estimation accuracy of ice sheet granular snow layer density, realizing accurate acquisition of regional-scale density distribution with depth, and providing reliable basic parameters for ice sheet engineering assessment and environmental evolution research. Attached Figure Description

[0059] Figure 1 This is an embodiment of the present invention, which provides a method for estimating the density of ice sheet particles and snow based on vehicle-mounted noise sources, and the received waveforms and time-spectrum characteristics of different channels of vehicle-mounted noise signals.

[0060] Figure 2 The following are examples of an embodiment of the present invention: (a) a reconstructed virtual shot gather record; (b) a dispersion image of the seismic record; and (c) and (d) dispersion energy maps of Rayleigh surface waves and P-guided waves corresponding to the seismic record enclosed by the orange and purple boxes in Figure (a).

[0061] Figure 3 This is an embodiment of the present invention, which provides the surface wave and P-guided wave dispersion curve inversion results and joint inversion fitting process diagram of an ice sheet snow layer density estimation method based on vehicle-mounted noise source;

[0062] Figure 4 This is a 1D P-wave and S-wave velocity and mean plot, as well as a standard deviation and coefficient of variation plot, of an ice sheet snow layer density estimation method based on vehicle-mounted noise source provided by an embodiment of the present invention.

[0063] Figure 5 This is an experimental result diagram of the regularization parameter of an ice sheet snow layer density estimation method based on vehicle-mounted noise source provided by an embodiment of the present invention;

[0064] Figure 6 This is a density-depth variation curve of an ice sheet snow layer density estimation method based on vehicle-mounted noise source provided by an embodiment of the present invention;

[0065] Figure 7 This is a flowchart illustrating a method for estimating the density of ice sheet particles based on a vehicle-mounted noise source, provided by an embodiment of the present invention.

[0066] Figure 8 This is a schematic diagram of an ice sheet snow layer density estimation system based on a vehicle-mounted noise source, provided by an embodiment of the present invention. Detailed Implementation

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

[0068] like Figure 7 As shown, an embodiment of the present invention proposes a method for estimating the density of ice sheet granular snow layer based on vehicle-mounted noise sources. The method includes the following steps:

[0069] Step 1: Identify and extract vehicle noise events from the original continuous noise records. Calculate the energy ratio of the short-term average to the long-term average and compare it with a preset threshold. Finally, output the exact start and end time periods of each vehicle noise event.

[0070] Step 2: Based on the exact start and end time periods of each vehicle noise event, extract the corresponding multi-channel vehicle noise event data from the original continuous noise record; perform mutual interference processing and linear superposition on the extracted multi-channel vehicle noise event data to obtain the virtual gun set record with improved signal-to-noise ratio.

[0071] Step 3: Extract the dispersion energy of Rayleigh surface wave and P-guided wave from the virtual shot gather record after the signal-to-noise ratio is improved. The P-guided wave refers to the longitudinal wave guided wave, which is a longitudinal wave guided wave mode that propagates in a layered medium and exhibits dispersion characteristics. The corresponding Rayleigh surface wave dispersion curve and P-guided wave dispersion curve are generated through dispersion analysis.

[0072] Step 4: Jointly invert the Rayleigh surface wave dispersion curve and the P-guided wave dispersion curve, construct and iteratively optimize the objective function to update the velocity model synchronously, and finally obtain a converged one-dimensional velocity structure model of the P-wave velocity and the S-wave velocity.

[0073] Step 5: By combining the measured density data of known borehole ice cores at multiple depth points with the P-wave velocity and S-wave velocity obtained by inversion at the corresponding depth, a power-law type empirical formula for density with P-wave velocity and S-wave velocity as independent variables is established.

[0074] Step 6: Substitute the P-wave velocity and S-wave velocity values ​​at each depth in the one-dimensional velocity structure model into the established power-law density empirical formula to calculate the density layer by layer, and finally output the density distribution of the ice sheet granular snow layer at the regional scale.

[0075] In this embodiment of the invention, by utilizing vehicle-mounted noise as a passive source, active excitation is not required, effectively reducing detection costs and operational difficulty. Vehicle-mounted noise events are accurately identified using the STA / LTA ratio method, and the signal-to-noise ratio of the virtual shot ensemble is reasonably improved through mutual interference processing and linear superposition, ensuring clear differentiation between Rayleigh surface wave and P-guided wave signals. Based on the joint inversion of Rayleigh surface waves and P-guided waves, P-wave and S-wave velocities are simultaneously obtained, overcoming the limitations of single wave velocity constraints and making the velocity model more closely match the physical characteristics of the snow layer. Regularized fitting is used to construct power-law empirical formulas for P-wave and S-wave velocities and density, effectively solving the ill-conditioned problems caused by variable correlation, resulting in stable and reliable fitting results. Combined with sparse borehole ice core data, regional-scale density estimation is achieved, covering a wide range and with high estimation accuracy, providing high-quality density parameter support for ice sheet engineering assessment and environmental evolution research.

[0076] In a preferred embodiment of the present invention, step 1 above may include:

[0077] Step 1.1: Obtain the original continuous noise record containing background noise and vehicle noise. Specifically, this includes: selecting a broadband seismic detector array, with the number of detectors set to 110 based on the detection area and resolution requirements, and the spacing between adjacent detectors uniformly set to 10m to ensure comprehensive capture of the propagation characteristics of vehicle noise; firmly coupling the detectors to the ice sheet surface to avoid signal distortion caused by environmental vibrations; configuring a data acquisition instrument with a sampling rate of 1000Hz and a sampling duration of 13 days to obtain the original continuous noise record containing low-amplitude, stable, randomly distributed background noise and vehicle noise generated by snowmobile driving. The time series of the original vibration signal received by each detector is represented by s(n), where n is the sampling point number, and each sampling point corresponds to a time interval of 0.001s. The original continuous noise record completely preserves the time and frequency domain characteristics of various vibration signals in the environment.

[0078] Step 1.2 calculates the short-term average energy and long-term average energy of the original continuous noise records. Specifically, for each detector's received original record s(n), within a preset candidate time period (set according to the expected passage time of the snowmobile in the detection plan, usually covering 4 to 12 hours), the records of all detectors corresponding to the channel are superimposed point by point at the same sampling time to obtain the array total amplitude energy sequence x(n) that changes continuously with time. This sequence cancels the random fluctuations of the background noise through the averaging effect, making the energy peak characteristics of the vehicle noise event more prominent. Since the superimposed energy curve still has occasional spike interference and short-period fluctuations of 1 to 2 seconds, a one-dimensional sliding median filter is applied to x(n). The filter window length is strictly set to 2 seconds. Based on a sampling rate of 1000Hz, it corresponds to 2000 continuous sampling points. The sliding method is to move the window point by point. Every time a sampling point is moved, the median of the 2000 energy values ​​in the current window is calculated and used to replace the energy value at the center of the window, finally obtaining a smooth and interference-free energy envelope curve.

[0079] Subsequently, based on the energy envelope curve, the short-term average energy STA is calculated using formula (1), and the long-term average energy LTA is calculated using formula (2), as follows:

[0080] (1)

[0081] (2)

[0082] In the formula, t is the current calculation time, x(tk) is the energy value of the energy envelope curve at time tk, and N s The short time window length ranges from 1 to 5 seconds, corresponding to 1000 to 5000 sampling points. It is adjusted based on the initial energy rise rate of the vehicle noise; a fast rise takes 1 to 2 seconds, and a slow rise takes 3 to 5 seconds. N l The time window length is set to 1-2 hours, based on the actual observation duration of vehicle noise of 200 to 300 seconds. This is to avoid the background energy estimation being contaminated by vehicle noise events, while ensuring that the dynamic changes of energy can be tracked quickly.

[0083] Step 1.3: Based on the calculated ratio of the short-term average energy to the long-term average energy, an energy ratio sequence is generated. Specifically, this includes: using a sliding step of one sampling point, based on a 1000Hz sampling rate, with a step size corresponding to 0.001s to ensure real-time capture of energy changes, calculating the ratio of the short-term average energy (STA) to the long-term average energy (LTA) at each calculated moment; that is, generating a real-time updated energy ratio sequence using the STA / LTA method. The core function of this sequence is to amplify the energy mutation characteristics when a vehicle noise event occurs. When there is no vehicle noise interference, the STA and LTA values ​​are close, and the energy ratio sequence R(t) is stable around 1.0 with a fluctuation range of no more than ±0.2. When a snowmobile approaches or passes through the receiving array, the energy enhancement caused by vehicle noise will cause the STA to rise rapidly in a short period of time, while the LTA remains relatively stable due to the smoothness of the long time window. This causes the energy ratio sequence R(t) to rise rapidly to 1.5 to 2.0 or above, forming a significant peak, thereby achieving the initial identification and location of vehicle noise events.

[0084] In a preferred embodiment of the present invention, step 2 above may include:

[0085] Step 2.1: Based on the exact start and end time periods of each output vehicle noise event, extract the corresponding data segments from the original continuous noise record to obtain multi-channel vehicle noise event data. Specifically, this includes: retrieving the output vehicle noise event structured information table; accurately associating each event's unique identifier, such as Event_001, with the complete record of the corresponding 110 detectors in the original continuous noise record; and based on the sampling point range n of each event in the information table... start up to n end The signal segment within the interval is extracted from each raw record using a data extraction tool, forming a single-event multi-channel data set with the number of channels exactly matching the number of detectors. After extraction, each data channel is zero-mean normalized by subtracting its arithmetic mean to eliminate DC offset interference. Simultaneously, the data amplitude is checked for any abnormal peaks exceeding three standard deviations; if any are found, linear interpolation is used to replace the outliers. Finally, the number of sampling points for each data channel is verified to be equal to n. end -n start +1 ensures that the length of each data channel perfectly matches the duration of the event, with no data loss, over-sampling, or data distortion, ultimately obtaining high-quality multi-channel vehicle noise event data corresponding to each vehicle noise event.

[0086] Step 2.2 involves preprocessing the multichannel vehicle noise event data to suppress non-stationary noise and instrument drift, improve the stability of subsequent interference processing, and obtain the preprocessed multichannel vehicle noise events. Specifically, this includes: performing mean-removal and de-trending processing on each channel record to eliminate DC components and slow drift; applying bandpass filtering to the records according to the target frequency band to retain the effective frequency range relevant to surface wave and guided wave dispersion analysis; optionally, to weaken the dominant role of amplitude transient events on the interference results, time-domain normalization (e.g., One-bit) is performed on the records.

[0087] Step 2.3 involves performing mutual interference processing on the preprocessed noise data to obtain the corresponding mutual interference results. Specifically, this includes: first, preprocessing the multi-channel vehicle noise event data for a single event to filter out low-frequency background interference; then, using Fast Fourier Transform (FFT) to convert the time-domain signal into a frequency-domain signal, with 1024 FFT points and a 1000Hz sampling rate to obtain the frequency-domain expression for each detector observation point. , Let A and B be any two detectors, and ω be the angular frequency in rad / s. Based on the mutual coherence method of seismic interferometry, substitute the following formula:

[0088] ;

[0089] Calculate the mutual interference of detector pair (A, B), where... for The conjugate of complex numbers, The cross-average power spectrum of the signals at points A and B was calculated using the Welch method (segment length 2500 points, overlap rate 0%). , The power spectra of the signals at points A and B are shown below. ε is a regularization parameter, fixed at 0.001 to avoid the denominator approaching zero, and |·| represents the complex modulus. Following the above procedure, cross-correlation calculations are sequentially performed on all pairwise detector combinations in the multichannel data to obtain the full-gather cross-correlation result corresponding to the vehicle noise event. The cross-correlation results of the multichannels are arranged according to the offset to obtain a pseudo-shot gather record containing all shots. This result fully preserves the signal propagation phase relationship and amplitude characteristics, laying the foundation for subsequent virtual shot gather reconstruction.

[0090] Step 2.3 involves linearly superimposing the cross-coherence results of multiple vehicle noise events to suppress random noise components, resulting in a superimposed cross-coherence result. Specifically, this includes collecting verified valid vehicle noise events and processing the superimposed full-gather cross-coherence results. An inverse fast Fourier transform (IFFT) is performed to convert the frequency domain signal back to the time domain signal, resulting in a multi-channel interference record in the time-offset domain. During the superposition process, the coherent signal energy of the vehicle noise is enhanced by superposition in the same direction, while the random noise (incoherent component) in each event cancels each other out due to the irregular phase. After superposition, the signal-to-noise ratio is verified by calculating the ratio of signal energy to noise energy. It is confirmed that the signal-to-noise ratio is significantly improved from 1.22 dB to 5.18 dB in a single event, and finally stabilizes at around 7.55 dB, achieving effective suppression of random noise.

[0091] After reconstruction, the virtual shot gather record was quality verified. Dispersion analysis confirmed that the dispersion energy of the surface wave and the guided wave was continuous and stable, meeting the requirements for subsequent dispersion curve extraction and joint inversion.

[0092] In a preferred embodiment of the present invention, step 3 above may include:

[0093] Step 3.1: Identify and extract Rayleigh surface wave signals and P-guided wave signals from the post-processed virtual shot gather records. Specifically, this includes: accurately identifying and extracting signals based on the differences in propagation characteristics between Rayleigh surface waves and P-guided waves (P-guided waves are volume waveguides with a propagation speed of 3000 to 4000 m / s; Rayleigh surface waves are surface waves with a speed of 1400 to 2000 m / s, and P-guided waves arrive first in the time domain); first, sort the virtual shot gather records by offset Δx from smallest to largest, and for each offset channel, calculate the dynamic time window based on the typical velocity range of the snow layer: the P-guided wave window is t p ±0.2+0.0001Δx,t p =Δx / V p V p Using 2500 m / s as the reference velocity, the window is widened slightly for larger offsets to avoid signal overflow. The Rayleigh surface wave window is t. r ±0.5+0.0002Δx,t r =Δx / V r V r 1500 m / s was taken as the reference velocity. After the two signals were separated by time window, the correlation between the extracted signal and the original channel signal was calculated to ensure that the correlation coefficient was ≥0.85 and to verify the integrity of the signal. At the same time, the two windows were ensured to have no overlap in the time domain with an interval of ≥0.3s to avoid signal crosstalk. Finally, the pure Rayleigh surface wave signal and P-guided wave signal corresponding to each offset channel were extracted.

[0094] Step 3.2 involves performing dispersion analysis on the extracted Rayleigh surface wave signal to obtain its dispersion energy distribution. Specifically, this includes generating a frequency-phase velocity dispersion energy map using the phase-shift method. For each Rayleigh surface wave time-domain signal with an offset of x... A one-dimensional Fourier transform along the time axis yields the corresponding complex spectrum in the frequency domain. , of which spectrum It is a complex spectrum, containing amplitude and phase information. ,in, Represents the phase spectrum. This represents the amplitude spectrum. The phase velocity range is scanned within a selected frequency range, for example (...). Hz, corresponding =2 f) and scanning the phase velocity c within the candidate phase velocity range, for each group The complex spectrum of each channel in the frequency domain is multiplied by a phase compensation term to achieve phase alignment at the candidate phase velocity, and then superimposed along the spatial gather to obtain the superimposed amplitude corresponding to the frequency and the phase velocity. To reduce the impact of inter-channel amplitude differences on the spectrum, it is possible to optionally adjust the... Normalization is performed to obtain the dispersion energy value of Rayleigh surface waves, and a frequency-phase velocity dispersion energy distribution map is formed accordingly.

[0095] Step 3.3, similarly, performs dispersion analysis on the extracted P-guided wave signal to obtain the dispersion energy distribution of the P-guided wave. Specifically, this includes: generating a frequency-phase velocity dispersion map of the P-guided wave signal using the phase-shifting method; performing a one-dimensional Fourier transform along the time axis on the time-domain signal of the P-guided wave for each offset x to obtain the corresponding frequency domain complex spectrum; and scanning the phase velocity range within a selected angular frequency range, for example (...). Hz, corresponding =2 f) Multiply the complex spectrum of each channel in the frequency domain by a phase compensation term to achieve phase alignment under the candidate phase velocity, and superimpose them along the spatial gather to obtain the superposition amplitude corresponding to the frequency and the phase velocity, thereby obtaining the dispersion energy value of the Rayleigh surface wave, and forming a frequency-phase velocity dispersion energy distribution map accordingly.

[0096] Step 3.4: Based on the dispersion energy distribution of the Rayleigh surface wave, generate the corresponding Rayleigh surface wave dispersion curve. Specifically, this includes: performing hierarchical screening and optimization of the Rayleigh surface wave dispersion energy distribution, first targeting each center frequency f. c The energy values ​​corresponding to all phase velocities at this frequency are statistically analyzed, and 70% of the maximum energy value is taken as the effective energy threshold for this frequency. Experimental verification shows that this can eliminate more than 80% of noise interference, retaining only the phase velocity range with energy above the threshold. For each f... c Within the effective range, select the phase velocity v corresponding to the maximum energy value. ϕ (f c This is used as the dispersion characteristic value for that frequency; anomaly judgment rules are set: if the dispersion value of a certain frequency differs from the average dispersion value of the two adjacent frequencies by more than 20%, and such fluctuations occur at three consecutive frequency points, then it is judged as an abnormal segment, and linear interpolation method is used. Fill the gap Center frequency f c The Rayleigh surface wave phase velocity obtained by down-interpolation is used to correct outliers or missing values; f c−1 f c The previous adjacent center frequency, f c+1 f c The next adjacent center frequency, and both of them with f c The intervals are all 1Hz; f c−1 The corresponding Rayleigh surface wave phase velocity, f c+1The corresponding Rayleigh surface wave phase velocities are all selected effective dispersion characteristic values. All effective frequencies from 5 to 50 Hz and their corresponding dispersion values ​​are sorted in ascending order of frequency, and the points are connected by a smooth curve. The generated Rayleigh surface wave dispersion curve must satisfy the physical law that low frequencies correspond to deeper layers and the velocity gradually decreases as the frequency increases, and the average deviation from the theoretical dispersion curve should not exceed 10%.

[0097] Step 3.5: Based on the dispersion energy distribution of the P-guided wave, generate the corresponding P-guided wave dispersion curve. Specifically, this includes: performing targeted processing on the P-guided wave dispersion energy distribution, first for each center frequency f... c The energy distribution is statistically analyzed, and 25% of the maximum value is taken as the effective threshold. Based on the characteristics of P-guided wave signals, which have a higher signal-to-noise ratio and more concentrated energy, this threshold can accurately preserve the effective signal; for each f c Extract the phase velocity v corresponding to the energy peak. ϕp (f c The dispersion characteristic value is used as the outlier. The outlier handling rule is as follows: if the dispersion value of a certain frequency differs from the dispersion value of the adjacent frequency by more than 30%, and the difference exceeds twice the standard deviation of the dispersion value of the frequency band, it is judged as an outlier. The cubic spline interpolation method is used for correction. The interpolation nodes are selected from two effective frequency points before and after the outlier to ensure that the overall trend of the dispersion curve is not changed. All effective frequencies from 40 to 120 Hz and their corresponding dispersion values ​​are arranged in ascending order of frequency and connected to form a continuous and stable P-guided wave dispersion curve. Finally, the stability of the curve is verified: different offset ranges, such as 0 to 500 m and 500 to 1000 m, are selected to extract dispersion curves. The dispersion value deviation of the corresponding frequency points does not exceed 15%, ensuring that the curve can truly reflect the longitudinal velocity structure of the snow layer and can be directly used for subsequent joint inversion.

[0098] In a preferred embodiment of the present invention, step 4 above may include:

[0099] Step 4.1: Based on the generated Rayleigh surface wave dispersion curve and the generated P-guided wave dispersion curve, an initial one-dimensional velocity model including P-wave velocity and S-wave velocity parameters is set. Specifically, this includes: based on the generated Rayleigh surface wave dispersion curve and the generated P-guided wave dispersion curve, combined with the physical evolution law of shallow loose and deep dense granular snow layer, and referring to prior information of this region, the initial one-dimensional velocity model is accurately set; the vertical stratification of the model adopts non-equidistant division, and after setting and adjusting based on the prior information of this region, the specific stratum thickness information is: 6m, 13m. The depths are 22m, 32m, 42m, 52m, 66m, and 90m, followed by a half-space. This stratification method fully adapts to the actual distribution characteristics and evolution of the local granular snow layer, taking into account the resolution requirements and computational efficiency of each depth layer. The initial shear wave velocity Vs is based on the velocity measurement data of a small number of rock cores in the reference area and existing granular snow layer research literature. For example, the shallow layer is set to 800 to 1200 m / s, the middle layer to 1200 to 1800 m / s, and the deep layer to 1800 to 2000 m / s. The specific model velocity parameters are based on the geological information of the test area. The initial P-wave velocity Vp is uniformly set to twice Vs. This ratio has been verified by previous regional tests and is consistent with the elastic characteristics of the local granular snow layer. The reference model parameter vector mref preferentially adopts the velocity profiles corresponding to 1 to 2 existing boreholes in the detection area. If none are available, a general velocity model of the same type of ice sheet granular snow layer is selected. The model parameter vector m is organized in the order of layer 1 Vp, layer 1 Vs, layer 2 Vp, layer 2 Vs, ..., layer n Vs, with a vector dimension of 2 × layer number (including parameters corresponding to half space). Finally, a clear initial one-dimensional velocity model that fits the local reality is formed, providing a reliable starting point for joint inversion.

[0100] Step 4.2: Construct a joint inversion objective function that includes a data fitting term and a model regularization term. The data fitting term is based on the difference between the dispersion curve obtained from the forward modeling of the initial one-dimensional velocity model and the observed dispersion curve. Specifically, this includes constructing a joint inversion objective function that includes both a data fitting term and a model regularization term, as shown in the formula: ; It is the joint inversion objective function. The square of the L2 norm is used to quantify the overall size of vectors such as residual vectors and model bias vectors. f(m) is the forward dispersion curve vector, and d is the observed dispersion data vector, which is composed of all the eigenvalues ​​of the extracted Rayleigh surface wave dispersion curve and P-guided wave dispersion curve, and serves as the reference benchmark for inversion. The data variance in W specifically refers to the sampling variance of the observed dispersion data of the Rayleigh surface wave and P-guided wave in step 3, reflecting the degree of dispersion of the observed data. The smaller the variance, the greater the weight and the stronger the constraint. The weight matrix W of the data fitting term is constructed as follows: first, the observed data variance of the obtained Rayleigh surface wave and P-guided wave dispersion curves is calculated, based on... The residual statistics of the scatter points are used, and the reciprocal of the variance is taken as the basic weight value. When updating the shear wave velocity Vs, the basic weight value of the Rayleigh surface wave is multiplied by 1.5, while the weight value of the P-guided wave is kept at 1.0. When updating the longitudinal wave velocity Vp, the basic weight value of the P-guided wave is multiplied by 1.5, while the weight value of the Rayleigh surface wave is kept at 1.0, to ensure a strong constraint match between the target parameters and the corresponding wave type. The regularization coefficient λ (λ is used to balance the weight coefficients of the data fitting term and the model regularization term. It is initially set to 0.01. The larger the value, the smoother the model. The smaller the value, the closer it is to the observed data) is 0.01. This value has been verified by trial calculation to avoid the model from being overly smooth or fluctuating drastically.

[0101] The smoothness matrix L (L is a matrix obtained by discretizing the Laplacian operator, with the same dimension as the model parameter vector m; longitudinal continuity of the model is controlled by constraining the difference in parameters between adjacent layers) adopts the discrete form of the Laplacian operator. Longitudinal continuity of the model is ensured by constraining the velocity difference between adjacent layers to not exceed 10% of the velocity of the previous layer; the reference model parameter vector m... ref If it is a general model, it needs to be scaled and adjusted according to the thickness and average density of the snow layer in the detection area to ensure that the prior information fits the local reality. The objective function minimizes the difference between the observed value and the forward value through the data fitting term, and constrains the physical rationality of the model through the regularization term, so as to achieve the goal of accurate fitting and stable and reliable inversion.

[0102] Step 4.3: The joint inversion objective function is optimized and iterated using a weighted least squares iterative algorithm, simultaneously updating the P-wave and S-wave velocity parameters in the one-dimensional velocity model. Specifically, this includes: First, based on the current velocity model m, the generalized ray method, suitable for layered media, is selected for dispersion forward modeling. This method can efficiently calculate the dispersion curves f(m) of Rayleigh surface waves and P-guided waves in layered granular snow media, and the calculation accuracy meets the inversion requirements. Then, the integral Jacobian matrix J is constructed, where J is... The sensitivity matrix, with dimensions N×M, represents the total number of observed dispersion data points, and M is the dimension of the model parameter vector m. Each element reflects the degree of influence of the corresponding model parameter on the corresponding dispersion point. It is the core matrix for solving the parameter update amount. During construction, a small perturbation of 1% is applied to each layer of model parameters Vp and Vs. This perturbation ratio has been verified to capture the sensitivity of the parameters to the dispersion curve without causing distortion of the forward modeling results. The change in the dispersion curve after each perturbation is recorded to form the sensitivity matrix. Based on the gradient information of the objective function, it is calculated using the following formula:

[0103] ;

[0104] Calculate the model parameter update amount Δm, with the same dimension as m, and store the adjustment values ​​of Vp and Vs for each layer. Positive values ​​indicate that the parameters need to be increased, and negative values ​​indicate that they need to be decreased. J T W is the transpose of the Jacobian matrix J, used for dimension matching in matrix operations. T It is the transpose of the weight matrix W, and J T In conjunction with completing the mapping of residual information to the model parameter space, L T Let be the transpose of the smoothness matrix L, which, together with λ, forms the matrix form of the regularization term, constraining the smoothness of the model during parameter update. (d−f(m)) is the residual vector, i.e., the difference between the observed dispersion data and the forward dispersion data, reflecting the fitting deviation between the current model and the observed data, and is the core driving force for parameter update. -1 It is a matrix inversion operation used to solve a system of linear equations and obtain an analytical solution for parameter updates.

[0105] A dynamic step size factor α is used when updating the model. α is an iteration step size control coefficient, ranging from 0.3 to 0.5, used to adjust the parameter update amplitude and avoid divergence due to excessively rapid updates during iteration: 0.3 is used for the first 10 iterations to prevent divergence in the initial stage, and 0.5 is used after 10 iterations to accelerate convergence. Update parameters, m k This is the model parameter vector after the k-th iteration, storing the Vp and Vs values ​​for each layer in the current iteration; m k+1 Let m be the model parameter vector after the (k+1)th iteration. kThe updated result after superimposing the correction amount; Vp and Vs are strictly updated synchronously in each iteration, and the Vp / Vs ratio is forced to be kept between 1.6 and 2.0, which is a reasonable range for the elastic parameters of the snow layer. If the ratio exceeds the range after the update, the Vs value is adjusted proportionally to ensure that the model parameters conform to the physical characteristics of the snow layer.

[0106] Step 4.4: After each optimization iteration, calculate the dispersion curve fitting error corresponding to the current model and determine whether the preset convergence condition has been met. Specifically, this includes: after each optimization iteration, using the error function... Calculate the fitting error of the dispersion curve of the current model; where N is the total number of observed dispersion data points, which are summed after being statistically analyzed for Rayleigh surface wave dispersion points and P-guided wave dispersion points to ensure statistical completeness; j is the index of the dispersion data, ranging from 1 to N, corresponding to the j-th specific dispersion point; It is the dispersion data of the j-th observation point; It is the dispersion data of the j-th calculation point; It is a summation operator, which means summing the normalized squared residuals of all dispersion points from j=1 to j=N; The average of the summed total squared residuals is obtained by averaging the squared residuals.

[0107] The three preset convergence conditions are further refined: First, the fitting error MF < 0.05, i.e., 5%. This threshold is based on the engineering accuracy requirements for estimating the density of ice sheet granular snow layer; a value below 5% is sufficient for subsequent density calculations. Second, the change in MF over three consecutive iterations. ΔMF represents the relative change in fitting error between two consecutive iterations, and is dimensionless; M Fk Let M be the fitting error of the k-th iteration. Fk+1 The fitting error of the (k+1)th iteration is used to determine whether the error tends to stabilize. If the error no longer decreases significantly, the inversion tends to stabilize. 3. The number of iterations reaches 50. In most scenarios, convergence can be achieved in 20 to 30 iterations. 50 iterations is the upper limit for extreme cases to avoid infinite iteration and consuming computing power. After each iteration, the error threshold is checked first, then the amount of change is checked, and finally the number of iterations is checked. As long as one of them is satisfied, it is immediately marked as convergence and the iteration is ready to be terminated.

[0108] Step 4.5: When the preset convergence condition is met, stop the iteration and output the finally converged one-dimensional velocity structure model of P-wave velocity and S-wave velocity. Specifically, this includes: when the preset convergence condition is met, stop the iteration and output the finally converged one-dimensional velocity structure model; the model output is further refined: including the precise depth range of each layer (e.g., 0 to 6.0m, 6.0 to 13.0m, etc.), the mean of Vp for each layer (rounded to one decimal place) and the sample standard deviation, the mean of Vs for each layer (rounded to one decimal place) and the sample standard deviation; simultaneously outputting the statistical results across the entire depth range: single-point values ​​of Vp and Vs for all depth points, as well as the overall mean ± standard deviation and coefficient of variation C. V = Standard Deviation / Mean × 100%; Additional physical validity verification is performed after output: First, the longitudinal velocity trend is checked to ensure that the velocity increases overall with depth, and the velocity gradient of each layer does not exceed 50 m / s, consistent with the physical process of gradual densification of the granular snow layer; Second, the upper limit of velocity is verified, with Vp in the deep layer not exceeding the longitudinal wave velocity of ice (approximately 4000 m / s) and Vs not exceeding the transverse wave velocity of ice (approximately 2000 m / s); Third, transverse non-uniformity is verified, with the CV value controlled between 1% and 3%, corresponding to the typical characteristics of layered granular snow layers; After successful verification, a one-dimensional velocity structure model with both accuracy and physical validity is finally output, providing accurate and reliable velocity parameter support for the subsequent construction of density empirical formulas.

[0109] In a preferred embodiment of the present invention, step 5 above may include:

[0110] Step 5.1: Obtain density information for key depth points that can be used to establish empirical density relationships based on borehole logging ice core data. Specifically, this includes: acquiring logging data and ice core sampling records from at least one borehole within the target detection area. The logging data is used to determine the ice core sampling depth, depth benchmark, and layer interface, among other depth information. Within the ice core sampling depth range, select several representative key depth points and record the measured density values ​​of the ice cores at these key depth points. Preferably, key depth points include ice core sample points near layer interfaces and in areas with significant density gradient changes. The depth identifier of each key depth point and its corresponding measured density value form a set of known density points. The set of known density points is then output as constraint data for subsequent fitting of empirical density relationships.

[0111] Step 5.2 involves statistically fusing the multiple velocity structures obtained from the joint inversion and establishing a one-to-one correspondence between the P-wave velocity and S-wave velocity at key depth points to form a ternary sample set for fitting. Specifically, this includes: obtaining the joint inversion results output in Step 4, which contain multiple one-dimensional velocity structure profiles; reading the P-wave velocity and S-wave velocity values ​​at the layer interfaces for each one-dimensional velocity structure profile; performing depth alignment on the multiple velocity structures using depth as a unified benchmark; and statistically averaging the depth-aligned one-dimensional velocity structures to obtain a representative total one-dimensional S-wave velocity structure. Preferably, the S-wave velocity at each depth point... Calculate the average shear wave velocity It can also simultaneously calculate dispersion indices (such as standard deviation) for quality control; similarly, the average P-wave velocity can be obtained. To improve the stability of the fit, it is optional to first... , Perform a logarithmic transformation. Bind the set of known ground density points obtained in step 5.1 to the total one-dimensional velocity structure: based on depth. As the unique correlation identifier, the logarithmic values ​​of the P-wave velocity and S-wave velocity at the corresponding depth are obtained from the total one-dimensional velocity structure; when When the velocity values ​​do not completely coincide with the formation depth grid points, interpolation between adjacent depth points is preferred to determine the corresponding velocity values. A ternary sample is constructed for each point with known density, in the form of a... ,in , This represents the total shear wave velocity structure corresponding to the key depth points. The measured density of the ice core at this depth is given. All ternary samples are compiled to form a data sample set for the next step of parameter fitting calculation.

[0112] Step 5.3: Construct a power-law model with P-wave velocity and S-wave velocity as independent variables and density as the dependent variable. Specifically, this includes: combining the physical evolution characteristics of the ice sheet granular snow layer, constructing a power-law model with P-wave velocity Vp and S-wave velocity Vs as independent variables and density ρ as the dependent variable. The formula remains simple and intuitive, i.e. Without adding complexity, the core reason for choosing the power-law model can be further refined as follows: the densification process of the granular snow layer is a nonlinear evolution process. As the depth increases, the granular snow particles are gradually compacted and the porosity continuously decreases. The longitudinal wave velocity and transverse wave velocity will increase nonlinearly. This nonlinear relationship between density and wave velocity cannot be accurately fitted by a linear model. The power-law model can match this characteristic of increasing wave velocity and power-law growth in density very well. Moreover, as verified by previous experiments, its goodness of fit is significantly higher than that of the linear model and the exponential model, and it is more in line with the actual physical state of the granular snow layer.

[0113] The physical meaning of each parameter in the model is further clarified: ρ is the measured or estimated density of the snow pellet layer at a certain depth point on the ice sheet, with the unit uniformly set to kg / m³. 3 This facilitates comparison with borehole measurement data; Vp is the P-wave velocity at this depth point, and Vs is the S-wave velocity, both in m / s, consistent with the units of the velocity model parameters output in step 4; a is the proportionality coefficient, reflecting the basic compactness of the granular layer, its value is related to the basic density of the shallow granular layer, and usually ranges from 0.02 to 0.06; b is the sensitivity coefficient of P-wave velocity to density, dimensionless, the larger the value, the more significant the effect of small changes in P-wave velocity on density, which fits the dominant influence of Vp on density during the densification process of the granular layer; c is the sensitivity coefficient of S-wave velocity to density, dimensionless, with a value less than b, reflecting the auxiliary influence of Vs on density, both are positive, and the value of b+c ranges from 0.8 to 1.2, consistent with the existing research conclusions on the relationship between granular layer wave velocity and density, ensuring the physical rationality of the model.

[0114] Step 5.4: Based on the data sample set, a regularized fitting method is used to determine the undetermined parameters in the power-law relationship model, generating a power-law relationship model with determined parameters. Specifically, this includes fitting the model based on the obtained effective data sample set. If there are three known ice core prior information points, theoretically, three sample points can achieve 100% accurate fitting, but ill-conditioned solutions and abnormal condition numbers are prone to occur. Therefore, the ridge regression L2 regularized fitting method is used to determine the undetermined parameters a, b, and c in the power-law model. Complex formulas are not used throughout the process; the focus is on clarifying the operational procedures and core requirements. First, the regularization parameters are determined by selecting parameters from 10... -8 Up to 10 0 Multiple trials were conducted to ensure optimal fitting results. Secondly, the natural logarithm of both sides of the power-law model was taken to transform the nonlinear model into a linear form, facilitating parameter solving via linear regression. During the transformation, Vp, Vs, and ρ of the three sample points were simultaneously transformed to obtain linear input data of lnVp, lnVs, and lnρ. Outliers (such as data that are meaningless after logarithmic transformation) were recorded and directly removed, then added to the outlier record in step 5.2. Subsequently, ridge regression with L2 regularization was used to solve for the undetermined parameters, effectively avoiding ill-conditioned solutions and abnormal condition numbers. The fitting results were stable, and the final goodness of fit R0 was [value missing]. 2 The value is 0.993. After fitting, verify the rationality of the parameters: ensure that b and c are both positive numbers, and that b+c is between 0.8 and 1.2, and a is between 0.02 and 0.06. If the parameters exceed this range, the fitting process needs to be optimized again (the value of the regularization parameter can be adjusted), eliminate data or method errors, and finally determine the specific values ​​of a, b, and c, retaining 3 decimal places, to generate a power-law relationship model with determined parameters.

[0115] Step 5.5: Output the power-law relationship model with determined parameters as the power-law empirical formula for density estimation. Specifically, this includes: after fitting, output the power-law relationship model with determined parameters, and determine it as the power-law empirical formula for subsequent density estimation. Focus on refining the output content, applicable scope, and usage instructions to ensure practicality and operability; confirm the complete form of the output empirical formula, and label the specific values ​​of parameters a, b, and c (keeping 3 decimal places).

[0116] Finally, in step 6, the specific values ​​of Vp and Vs at a certain depth point in the one-dimensional velocity structure model are extracted to ensure that the values ​​are within the applicable range. They are then directly substituted into the empirical formula to calculate the estimated value of the snow layer density at that depth point. At the same time, supplementary verification is provided: the three drilling core points used in step 5 are selected as prior information for error calculation and physical property comparison.

[0117] In a preferred embodiment of the present invention, step 6 above may include:

[0118] Step 6.1: For each depth point in the one-dimensional velocity structure model, read the corresponding P-wave velocity and S-wave velocity values. Specifically, this includes: retrieving the output final converged one-dimensional velocity structure model, where each depth point is bound to a unique depth identifier and the corresponding effective values ​​of P-wave velocity Vp and S-wave velocity Vs; reading data from each depth point in ascending order of depth value, strictly associating the depth identifier during reading, and synchronously recording the depth point number, depth value, and corresponding Vp and Vs values; after reading, validating each set of velocity data to confirm whether it falls within the applicable velocity range of the power-law density empirical formula determined in step 5.5, i.e., Vp is between 2000 and 4000 m / s and Vs is between 1500 and 2000 m / s. Velocity data outside this range are marked as points to be checked and will not participate in subsequent density calculations, ensuring that the input velocity data meets the requirements of the formula application.

[0119] Step 6.2: Substitute the P-wave velocity and S-wave velocity values ​​read from each depth point into the power-law density empirical formula obtained above to calculate the density estimate for that depth point. Specifically, this includes substituting the read and verified P-wave velocity and S-wave velocity values ​​from each depth point into the power-law density empirical formula output in step 5.5 to calculate the density estimate. During the calculation, the power terms are calculated first, followed by the product terms. The density estimate for each depth point is rounded to three decimal places. After the calculation is complete, the results are checked for reasonableness to confirm whether the density estimate falls within the reasonable density range of 400 to 900 kg / m³ for the snow layer. 3 This ensures that reliable density estimation results are obtained for each qualified depth point.

[0120] Step 6.3: Collect the density estimation results from each depth point to form a density-depth distribution dataset, and complete consistency adjustments and supplementary verification to enhance the reliability of the results. Specifically, this includes: using depth identifiers as the core association basis, binding all depth point information with the corresponding calculated density estimates one by one. Each depth point generates a dataset containing a depth number, depth value (m), P-wave velocity value (m / s), S-wave velocity value (m / s), and density estimate (kg / m³). 3 The data records are multi-dimensional and arranged in ascending order of depth value. Furthermore, to verify the reliability and physical consistency of the output results, leave-one-out (LOO) cross-validation is used to evaluate the generalization error of the empirical formula, and the cross-validation error is used as a reference for the uncertainty of density estimation. Elastic parameters such as the elastic modulus and Poisson's ratio are derived based on the multi-dimensional data record set, and consistency is checked in conjunction with the geological background information of the experimental area to verify the rationality of the velocity-density relationship at the elastic physics level. The two-parameter power-law empirical formula of this invention is compared with the single-parameter density wellbore formula proposed by predecessors. If the constraint results of the two in the density profile are consistent or more stable in trend and magnitude, the effectiveness of the established relationship is further proven.

[0121] Step 6.4 outputs a distribution dataset representing the density variation with depth, serving as the density distribution result of the ice sheet frit layer at the regional scale. This includes: standardizing the resulting density distribution dataset, determining a unified format for the output data, listing the core fields of depth value and density estimate, and including supplementary data information such as the total number of data points, the number of valid data points, the depth coverage range, and the mean and standard deviation of the density estimates; formally outputting the standardized distribution dataset as the density distribution result of the ice sheet frit layer at the regional scale, simultaneously indicating the applicable range of the results, i.e., the 0 to 100m ice sheet frit layer depth range consistent with velocity models and empirical formulas, and specifying the application boundary of the results, i.e., beyond this depth range, the frit layer approaches dense ice, the density variation pattern changes, and the results are no longer applicable; the output distribution results can be directly used for subsequent work such as thickness analysis, porosity evolution research, and mass estimation of the ice sheet frit layer, providing accurate density distribution data support for the study of the physical properties of the ice sheet frit layer at the regional scale.

[0122] In the specific implementation process of this embodiment, the following processes and methods are also included:

[0123] The specific implementation method and process mainly consists of three parts: First, the extraction and processing of vehicle-mounted noise to reconstruct virtual artillery records; second, the joint inversion of Rayleigh surface waves and P-guided waves to calculate the longitudinal and transverse wave velocities of the snow layer; and third, the fitting of the power-law function relationship between velocity and density to obtain a new method for estimating the density of the underground medium.

[0124] First, virtual shot set reconstruction of vehicle-mounted noise is performed: This involves analyzing the vehicle-mounted noise generated by some polar snowmobiles, such as... Figure 1 When no vehicles are passing, the signal exhibits low amplitude and a stable, randomly distributed background noise. When a snowmobile passes the receiver array, a typical spindle-shaped vibration is detected, with its energy peak time systematically delayed as the offset distance increases, reflecting the relative motion between the vehicle and the detector array. To identify and extract vehicle noise events from long-term sequences, we use the following process: First, the original records s(n) received by each detector are superimposed on the records of each channel within the candidate time period to obtain the array total amplitude energy sequence x(n) that changes over time. This curve averages out the random fluctuations of the background noise, making the vehicle noise events more prominent. Since the superimposed energy curve still contains spikes and short-period fluctuations, we apply a one-dimensional median filter to it to obtain a relatively smooth energy envelope. The time window length of the median filter can generally be selected as 2s. The ratio of the average of the short time window and the long time window is calculated using formulas (1) and (2), respectively. In order to avoid the background estimation being contaminated by the event energy, and at the same time to track the rapid changes in energy, it is recommended that the long time window length be 2 to 3 times greater than the length of the complete vehicle noise, and the short time window length be 1 to 5s. When the STA / LTA ratio first exceeds the threshold, it is recorded as the start time of traffic noise. When the ratio drops to the de-triggering threshold, it is defined as the end of the event.

[0125]

[0126] Where N s N represents the short time window length. l This indicates the length of the long time window; the obtained trigger time information is used to truncate the original data, retaining the time period of vehicle noise events. Since traffic noise often has strong coherence, seismic interferometry using the mutual coherence method is employed to reconstruct a virtual shot gather record similar to the active source seismic record.

[0127] (3)

[0128] in , Let A and B represent the frequency domain Green's functions at observation points A and B, respectively, where ω is the angular frequency. It is a regularization parameter. Represents the modulus of a complex number. This indicates the calculation of the average power spectrum. The results show the mutual interference between detectors A and B. The signal-to-noise ratio (SNR) of the pseudo-shot gather record for a single event in the seismic interferometric reconstruction is approximately 1.22 to 5.18 dB. By linearly superimposing the mutual interference results of multiple vehicle-mounted noise events, random noise can be effectively suppressed, and the SNR of the virtual shot gather is increased to 7.55 dB. Figure 2 In a, Rayleigh surface waves and P-guided waves can be clearly identified in the recording, and the surface waves and guided waves have continuous and clear dispersive energy.

[0129] Second, the next step is the joint inversion of P-guided waves and surface waves: a comprehensive inversion method using the dispersion curves of P-guided waves and surface waves is employed for the reconstructed pseudo-shot gather record to simultaneously estimate the Vp and Vs models. This method requires constructing an integral Jacobian matrix to perform a global inversion of the dispersion curves of guided waves and surface waves. The objective function Ф(m) includes data fitting terms and model regularization terms.

[0130] (4)

[0131] f(m) represents the dispersion curve obtained from forward modeling, d represents the observed dispersion data vector, and m represents the model parameter vector. ref denoted as the reference model parameter vector containing prior information, W represents the weight matrix composed of the inverses of the data variances, when updating the Vs results, the surface wave dispersion has a higher weight, and when updating the Vp results, the P guided wave dispersion has a higher weight. Let L represent the square of the L2 norm; L represents the smoothness matrix formed by the discrete form of the Laplacian operator. After multiple trials, a reasonable initial model was set, and the values ​​of the model parameters were updated through the least squares inversion model, so that the inversion process converged rapidly within a finite number of iterations. At the same time, the convergence of the inversion was judged using the error function (5):

[0132] (5)

[0133] in and These are the dispersion data of the j-th observation point and the calculation point, respectively. Iterative inversion stops when one of the following three stopping conditions is met: 1. MF is less than a predefined threshold; 2. MF has converged and no longer decreases significantly; 3. The number of iterations reaches the predefined maximum number of iterations. After joint inversion, a dispersion curve with good fitting can be obtained, such as... Figure 3 Obtain the 1D P-wave and S-wave velocity structure of the survey line area, such as... Figure 4In the figure, the light-colored solid circles represent the single-point values ​​at various depths of the 15 profiles, while the hollow symbols with error bars represent the mean and standard deviation. The standard deviation at each depth is approximately on the order of tens of m / s, which is very small relative to the average velocity of P-waves and S-waves. Longitudinal stratification is evident, with velocity increasing overall with depth, and the velocity gradient corresponds to the physical process of gradual densification of granular snow with depth. Lateral non-uniformity is relatively weak overall (CV is mainly between 1% and 3%), corresponding to a layered granular snow layer structure.

[0134] Third, we then began constructing an empirical formula for density calculation based on Vp-Vs: In the next step of deriving the relationship between velocity and density, we were inspired by the power-law relationship between P-wave velocity and density proposed by Gardner et al. (1974). We assumed that density could be expressed as a power function of Vp and Vs, while introducing constraints on density from P-wave and S-wave velocities, and using limited core data for local calibration:

[0135] (6)

[0136] Where a, b, and c are undetermined coefficients; this form empirically allows for different power-law sensitivities of density to P-wave and S-wave velocities. To transform the nonlinear relationship into a linear problem, we take the logarithm of Equation 6 to obtain:

[0137] (7)

[0138] make , , , Then for the i-th calibration point, we have:

[0139] (8)

[0140] Written in matrix form:

[0141] (9)

[0142] Theoretically, Equation 7 can provide an exact solution when there are three sets of sample points. However, Vp and Vs may exhibit a strong correlation in the granular snow medium, causing the eigenvalues ​​of the design matrix X to approach 0 and the condition number to increase abnormally. This makes ordinary least squares estimation highly sensitive to noise, resulting in ill-formed solutions. To improve the stability of the solution and avoid non-physical solutions, we introduce ridge regression as a regularization constraint. The parameters to be estimated are obtained by minimizing the following equation:

[0143] (10)

[0144] Where λ is the regularization parameter. Its analytical solution is:

[0145] (11)

[0146] Where I is the identity matrix. By adjusting λ, the ill-conditioned problems caused by the limited sample size and strong correlation of independent variables can be effectively alleviated, resulting in smooth and physically reasonable fitting parameters. After testing, as shown... Figure 5 When λ is 10 -3 At that time, it was located at the inflection point of the L-shape, and the parameter norm, condition number, and error term were all within a stable range. The inverted velocity information was then compared with the density information at three known drilling depths (6m: 550 kg / m³). 3 51m: 830kg / m 3 89m: 901kg / m 3 By fitting the data and using formula 11, the exponents Vp and Vs are calculated to be 0.545 and 0.407, respectively. The final relationship is as follows:

[0147] (12)

[0148] We take the density of ice as ρ ice =915kg / m 3 Based on the polycrystalline ice modulus K=8.90GPa and μ=3.46GPa, the calculated result is taken as... , Combining the velocity-density empirical formulas proposed by Kohnen and Diez in 1972 and 2014 respectively, density-depth variation curves derived from the three empirical formulas were plotted, as follows: Figure 6 The trend of Equation 12 with depth is a combination of the curves of the other two. In the shallow layer above 20m, it is closer to the derivation result of Diez's Vs-density formula, in the middle layer it is closer to the derivation result of Kohnen's Vp-density formula, and in the deep layer it has good similarity with both, which verifies the effectiveness of the method.

[0149] like Figure 8 As shown, embodiments of the present invention also provide an ice sheet snow layer density estimation system based on vehicle-mounted noise sources, comprising:

[0150] The noise event extraction and processing module is used to identify and extract vehicle noise events from the original continuous noise records. By calculating the energy ratio of the short-term window average to the long-term window average and comparing it with a preset threshold, the module finally outputs the exact start and end time periods of each vehicle noise event.

[0151] The virtual gun set reconstruction module is used to extract the corresponding multi-channel vehicle noise event data from the original continuous noise record based on the exact start and end time period of each vehicle noise event; the extracted multi-channel vehicle noise event data is subjected to mutual interference processing and linear superposition to obtain the virtual gun set record with improved signal-to-noise ratio.

[0152] The multi-wave dispersion curve extraction module is used to extract the dispersion energy of Rayleigh surface wave and P-guided wave from the virtual shot gather record after the signal-to-noise ratio is improved, and generate the corresponding Rayleigh surface wave dispersion curve and P-guided wave dispersion curve through dispersion analysis.

[0153] The P-wave and S-wave velocity joint inversion module is used to jointly invert the Rayleigh surface wave dispersion curve and the P-guided wave dispersion curve. By constructing and iteratively optimizing the objective function, the velocity model is updated synchronously, and finally a converged one-dimensional velocity structure model of P-wave velocity and S-wave velocity is obtained.

[0154] The empirical formula construction module is used to establish a power-law type empirical formula for density with P-wave velocity and S-wave velocity as independent variables by combining the measured density data of known borehole ice cores at multiple depth points with the P-wave velocity and S-wave velocity obtained by inversion at the corresponding depth.

[0155] The density calculation and output module is used to substitute the P-wave velocity and S-wave velocity values ​​at each depth in the one-dimensional velocity structure model into the established power-law density empirical formula to calculate the density layer by layer, and finally output the density distribution of the ice sheet granular snow layer at the regional scale.

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

Claims

1. A method for estimating the density of ice sheet granular snow layer based on vehicle-mounted noise sources, characterized in that, The method includes: Vehicle noise events are identified and extracted from the original continuous noise records. The energy ratio of the short-term average to the long-term average is calculated and compared with a preset threshold. Finally, the exact start and end time periods of each vehicle noise event are output. Based on the exact start and end time periods of each vehicle noise event, corresponding multi-channel vehicle noise event data are extracted from the original continuous noise record; the extracted multi-channel vehicle noise event data are subjected to mutual interference processing and linear superposition to obtain a virtual gun set record with improved signal-to-noise ratio. The dispersion energy of Rayleigh surface wave and P-guided wave is extracted from the virtual shot gather record after the signal-to-noise ratio is improved, and the corresponding Rayleigh surface wave dispersion curve and P-guided wave dispersion curve are generated by dispersion analysis. The Rayleigh surface wave dispersion curve and the P-guided wave dispersion curve are jointly inverted. By constructing and iteratively optimizing the objective function to update the velocity model synchronously, a convergent one-dimensional velocity structure model of the P-wave velocity and the S-wave velocity is finally obtained. By combining the measured density data of known borehole ice cores at multiple depth points with the P-wave velocity and S-wave velocity obtained by inversion at the corresponding depth, a power-law type empirical formula for density with P-wave velocity and S-wave velocity as independent variables is established. Substituting the P-wave velocity and S-wave velocity values ​​at each depth in the one-dimensional velocity structure model into the established power-law density empirical formula, the density is calculated layer by layer, and finally the density distribution of the ice sheet granular snow layer at the regional scale is output.

2. The method for estimating ice sheet snow density based on vehicle-mounted noise sources according to claim 1, characterized in that, Vehicle noise events are identified and extracted from raw continuous noise records. The energy ratio of the short-term average to the long-term average is calculated and compared with a preset threshold. Finally, the exact start and end time periods of each vehicle noise event are output, including: Obtain the original continuous noise record containing background noise and vehicle noise; Calculate the average energy of the short-time window and the average energy of the long-time window of the original continuous noise record; An energy ratio sequence is generated based on the calculated ratio of the average energy of the short time window to the average energy of the long time window. The energy ratio sequence is compared with a preset trigger threshold to determine the start time of the vehicle noise event; After determining the start time of the vehicle noise event, the system continuously monitors the time when the energy ratio sequence drops to a preset detriggered threshold to determine the end time of the vehicle noise event. Based on the determined start and end times of the vehicle noise events, the exact start and end time periods of each vehicle noise event are output.

3. The method for estimating the density of ice sheet granular snow layer based on vehicle-mounted noise sources according to claim 2, characterized in that, Based on the exact start and end time periods of each vehicle noise event, the corresponding multi-channel vehicle noise event data is extracted from the original continuous noise records. The captured multichannel vehicle noise event data are subjected to mutual interference processing and linear superposition to obtain a virtual shot set record with improved signal-to-noise ratio, including: Based on the exact start and end time periods of each output vehicle noise event, the corresponding data segments are extracted from the original continuous noise records to obtain multichannel vehicle noise event data. Multichannel vehicle noise event data are subjected to mutual interference processing to obtain the corresponding mutual interference results; The mutual interference results of multiple vehicle noise events are linearly superimposed to suppress random noise components, resulting in the superimposed mutual interference results. Based on the superimposed mutual coherence results, a virtual shot gather record with clear surface wave and guided wave signals is reconstructed.

4. The method for estimating ice sheet snow density based on vehicle-mounted noise sources according to claim 3, characterized in that, The dispersion energy of Rayleigh surface waves and P-guided waves was extracted from the virtual shot gather record after signal-to-noise ratio enhancement. Dispersion analysis was then used to generate the corresponding Rayleigh surface wave dispersion curves and P-guided wave dispersion curves, including: The virtual shot gather record is preprocessed to enhance the Rayleigh surface wave and P-guided wave signals contained therein, and the preprocessed virtual shot gather record is obtained. Rayleigh surface wave signal and P-guided wave signal were identified and extracted from the preprocessed virtual shot gather record; Dispersion analysis was performed on the extracted Rayleigh surface wave signal to obtain the dispersion energy distribution of the Rayleigh surface wave; The extracted P-guided wave signal was subjected to dispersion analysis to obtain the dispersion energy distribution of the P-guided wave. Based on the dispersion energy distribution of Rayleigh surface waves, the corresponding Rayleigh surface wave dispersion curves are generated. Based on the dispersion energy distribution of the P-guided wave, the corresponding P-guided wave dispersion curve is generated.

5. The method for estimating the density of ice sheet granular snow layer based on vehicle-mounted noise sources according to claim 4, characterized in that, By jointly inverting the Rayleigh surface wave dispersion curve and the P-guided wave dispersion curve, and constructing and iteratively optimizing the objective function to synchronously update the velocity model, a convergent one-dimensional velocity structure model for the P-wave velocity and S-wave velocity is finally obtained, including: Based on the generated Rayleigh surface wave dispersion curve and the generated P-guided wave dispersion curve, an initial one-dimensional velocity model including P-wave velocity and S-wave velocity parameters is set. Construct a joint inversion objective function that includes a data fitting term and a model regularization term; wherein, the data fitting term is based on the difference between the dispersion curve obtained by forward modeling of the initial one-dimensional velocity model and the observed dispersion curve; The joint inversion objective function is optimized and iterated using a weighted least squares iterative algorithm, and the P-wave velocity and S-wave velocity parameters in the one-dimensional velocity model are updated synchronously. After each optimization iteration, calculate the dispersion curve fitting error corresponding to the current model and determine whether the preset convergence condition has been met. When the preset convergence condition is met, the iteration stops and the finally converged one-dimensional velocity structure model of the P-wave velocity and S-wave velocity is output.

6. The method for estimating ice sheet snow density based on vehicle-mounted noise sources according to claim 5, characterized in that, By combining measured density data from known borehole ice cores at multiple depths and fitting them with the P-wave and S-wave velocities obtained from inversion at the corresponding depths, a power-law empirical formula for density with P-wave and S-wave velocities as independent variables is established, including: Obtain the P-wave velocity and S-wave velocity at each depth in the finally converged one-dimensional velocity structure model, as well as the measured density data of the known borehole ice core at the corresponding depth points; Establish a correspondence between the P-wave velocity and S-wave velocity at the same depth point and the measured density data to form a data sample set for fitting. A power-law model is constructed with P-wave velocity and S-wave velocity as independent variables and density as the dependent variable. Based on the data sample set, a regularized fitting method is used to determine the undetermined parameters in the power law relationship model, and a power law relationship model with determined parameters is generated. The output parameter-defined power-law relationship model serves as a power-law empirical formula for density estimation.

7. The method for estimating ice sheet snow density based on vehicle-mounted noise sources according to claim 6, characterized in that, Substituting the P-wave velocity and S-wave velocity values ​​at each depth in the one-dimensional velocity structure model into the established power-law density empirical formula, the density is calculated layer by layer. The final output is the density distribution of the ice sheet granular layer at the regional scale, including: For each depth point in the one-dimensional velocity structure model, read the corresponding longitudinal wave velocity value and transverse wave velocity value at that depth point. Substitute the P-wave velocity and S-wave velocity values ​​of each depth point into the power-law density empirical formula to calculate the density estimate of the depth point. The depth information of all depth points and their corresponding density estimates are collected to form a distribution data set that characterizes the density variation with depth; Output a set of distribution data representing the density variation with depth, as the density distribution of ice sheet granular snow layer with depth at the regional scale.

8. A system for estimating the density of ice sheet granular snow layer based on vehicle-mounted noise sources, the system implementing the method as described in any one of claims 1 to 7, characterized in that, include: The noise event extraction and processing module is used to identify and extract vehicle noise events from the original continuous noise records. By calculating the energy ratio of the short-term window average to the long-term window average and comparing it with a preset threshold, the module finally outputs the exact start and end time periods of each vehicle noise event. The virtual gun set reconstruction module is used to extract the corresponding multi-channel vehicle noise event data from the original continuous noise record based on the exact start and end time period of each vehicle noise event; the extracted multi-channel vehicle noise event data is subjected to mutual interference processing and linear superposition to obtain the virtual gun set record with improved signal-to-noise ratio. The multi-wave dispersion curve extraction module is used to extract the dispersion energy of Rayleigh surface wave and P-guided wave from the virtual shot gather record after the signal-to-noise ratio is improved, and generate the corresponding Rayleigh surface wave dispersion curve and P-guided wave dispersion curve through dispersion analysis. The P-wave and S-wave velocity joint inversion module is used to jointly invert the Rayleigh surface wave dispersion curve and the P-guided wave dispersion curve. By constructing and iteratively optimizing the objective function, the velocity model is updated synchronously, and finally a converged one-dimensional velocity structure model of P-wave velocity and S-wave velocity is obtained. The empirical formula construction module is used to establish a power-law type empirical formula for density with P-wave velocity and S-wave velocity as independent variables by combining the measured density data of known borehole ice cores at multiple depth points with the P-wave velocity and S-wave velocity obtained by inversion at the corresponding depth. The density calculation and output module is used to substitute the P-wave velocity and S-wave velocity values ​​at each depth in the one-dimensional velocity structure model into the established power-law density empirical formula to calculate the density layer by layer, and finally output the density distribution of the ice sheet granular snow layer at the regional scale.

Images