Real-time prediction and multi-dimensional intelligent evaluation method for compaction quality of soft soil subgrade

Abstract

This invention belongs to the field of data processing and analysis technology, specifically relating to a real-time prediction and multi-dimensional intelligent evaluation method for compaction quality of soft soil subgrades. The method includes: Step 1: Deploying multiple vertical component engineering seismic detectors along the direction of travel on the soil surface behind the vibratory roller wheel to form a detector array, with all sensors connected to an edge computing gateway for synchronous data acquisition; Step 2: Performing filtering processing on the raw seismic waveform records acquired by the detector array; Step 3: Establishing an initial velocity model based on the Rayleigh wave dispersion curve and outputting a final shear wave velocity model; Step 4: Extracting the estimated compaction coefficient and compaction uniformity variation coefficient for each depth layer based on the final shear wave velocity model, and outputting the compaction quality grade according to preset judgment rules. This invention can transform quality inspection from traditional post-construction sampling to real-time monitoring during construction, helping to promptly identify and address compaction quality problems.

Description

Technical Field

[0001] This invention belongs to the field of intelligent sensing systems, and further to the field of data processing and analysis technology, specifically involving a method for real-time prediction and multi-dimensional intelligent evaluation of compaction quality for soft soil subgrades. Background Technology

[0002] Seismic wave method is a geophysical exploration method based on the principle of elastic wave propagation. It inverts the underground velocity structure by analyzing the propagation characteristics of seismic waves in soil. Surface wave exploration is an important branch of seismic wave method. It utilizes the dispersion characteristics of Rayleigh waves propagating along the Earth's surface to obtain the shear wave velocity profile of the soil by extracting the phase velocities of Rayleigh waves at different frequencies and performing inversion. Since shear wave velocity is directly related to the shear modulus of the soil, and the shear modulus increases with the degree of compaction, shear wave velocity can be used as an effective indicator to characterize the compaction state of the soil. Multichannel surface wave analysis is currently a widely used method for dispersion curve extraction. It involves deploying a detector array on the surface to receive surface wave signals excited by artificial sources, extracting dispersion curves using signal processing techniques such as frequency wavenumber analysis or phase shift transformation, and then inverting the shear wave velocity profile using linearization inversion or optimization methods such as genetic algorithms. However, traditional surface wave exploration methods require independent artificial seismic sources such as drop hammers or explosives, and construction work needs to be stopped during detection, making it impossible to achieve real-time monitoring of the compaction process; dispersion curve inversion usually adopts the assumption of a horizontal layered model, which makes it difficult to accurately characterize the lateral velocity changes caused by uneven compaction; the resolution of the inversion results is limited by the bandwidth of the dispersion curve and the length of the detector array, and the accuracy of identifying local compaction anomalies needs to be improved.

[0003] Full waveform inversion is a high-resolution seismic imaging method that has emerged in recent years. It reconstructs subsurface velocity models by minimizing the differences between observed and synthesized waveforms. Compared with traditional inversion methods that only utilize travel time or dispersion information, full waveform inversion fully utilizes the amplitude, phase, and waveform details contained in the seismic waveform, theoretically achieving spatial resolution close to half-wavelength scale. However, full waveform inversion faces many challenges in engineering applications: it requires enormous computation, necessitating multiple wavefield forward models and gradient calculations, and traditional implementation methods struggle to meet the timeliness requirements of real-time on-site detection; the inversion problem exhibits strong nonlinear characteristics and multiple solutions, is highly dependent on the initial model, and is prone to getting trapped in local minima; and the uniform grid discretization method uses the same spatial resolution across the entire computational domain, which can lead to wasted computational resources and insufficient resolution in cases where local velocity anomalies exist. Summary of the Invention

[0004] The main objective of this invention is to provide a real-time prediction and multi-dimensional intelligent evaluation method for compaction quality of soft soil subgrades. It utilizes the vibration of the roller itself to achieve synchronous detection of the compaction process without the need for an independent vibration source device. The adaptive mesh refinement strategy achieves high resolution in areas of uneven compaction while effectively controlling the computational burden. The multi-dimensional evaluation system comprehensively considers the degree of deep stratified compaction and overall uniformity, enabling a comprehensive and accurate evaluation of the compaction quality of soft soil subgrades.

[0005] To address the aforementioned technical problems, this invention provides a method for real-time prediction and multi-dimensional intelligent evaluation of compaction quality for soft soil subgrades. This method includes:

[0006] Step 1: Multiple vertical component engineering seismic detectors are deployed on the soil surface behind the vibratory roller wheel along the direction of travel to form a detector array. An accelerometer is installed on the bearing seat of the vibratory roller wheel to record the vibration characteristics of the vibratory roller wheel. All sensors are connected to the edge computing gateway for synchronous data acquisition.

[0007] Step 2: Filter the raw seismic waveform records acquired by the detector array, and use the phase shift transformation method to extract the Rayleigh wave dispersion curve from the filtered multi-channel waveform records;

[0008] Step 3: Establish an initial velocity model based on the Rayleigh wave dispersion curve. Through a cyclic iterative process of wavefield forward modeling, adjoint wavefield back propagation, sensitive core field calculation, identification of highly sensitive anomaly elements and adaptive mesh refinement, velocity model iterative update, and iterative convergence determination, the final shear wave velocity model is output.

[0009] Step 4: Extract the estimated compaction coefficient and compaction uniformity variation coefficient of each depth layer based on the final shear wave velocity model, and output the compaction quality grade according to the preset judgment rules.

[0010] Furthermore, in step one, the detector array includes 16 vertical component engineering seismic detectors. The first vertical component engineering seismic detector is 1 meter away from the rear edge of the vibrating wheel, and the spacing between adjacent vertical component engineering seismic detectors is 0.3 meters. One triaxial microelectromechanical accelerometer is installed on the left and right sides of the vibrating wheel bearing seat. The synchronous acquisition module of the edge computing gateway synchronously acquires data from all channels using a sampling frequency of 4000 Hz, and uses the zero-crossing rising edge of the vertical component of the left triaxial microelectromechanical accelerometer in each excitation cycle as the synchronous triggering reference for all channels.

[0011] Furthermore, in step two, a fourth-order Butterworth filter is used for filtering, with a passband range of 5 to 200 Hz. The specific process of the phase shift transformation method is as follows: a fast Fourier transform is performed on each waveform record to obtain a complex spectrum. Each analysis frequency is selected sequentially within the frequency range of 5 to 100 Hz, and each test phase velocity is selected sequentially within the phase velocity range of 50 to 350 m / s. For each combination of analysis frequency and test phase velocity, the theoretical propagation time is calculated based on the horizontal distance between the engineering seismic detector and the vibrating wheel for each vertical component and the current test phase velocity, and converted into a phase offset. The complex value at the current analysis frequency is extracted from the complex spectrum of each waveform record and multiplied by the complex exponential factor corresponding to the phase offset to complete the phase correction. All 16 phase-corrected complex values ​​are summed and the modulus is taken as the dispersion energy value. The test phase velocity corresponding to the maximum dispersion energy value at each analysis frequency is selected as the Rayleigh wave phase velocity at the current analysis frequency.

[0012] Furthermore, in step three, the initial velocity model is established as follows: the rectangular roadbed area directly below the vibrating wheel, with a horizontal width of 2.4 meters and a depth of 1.2 meters, is divided into an initial inversion grid. The initial inversion grid is uniformly divided into 8 columns of elements in the horizontal direction and 6 rows of elements in the depth direction, thus forming 48 rectangular grid elements. An initial shear wave velocity value is set for each rectangular grid element according to the Rayleigh wave dispersion curve. The setting method is as follows: calculate the depth value of the center point of each rectangular grid element, find the Rayleigh wave phase velocity value that matches the wavelength of the depth value of the current rectangular grid element in the Rayleigh wave dispersion curve, and multiply the current Rayleigh wave phase velocity value by a coefficient of 0.92 to obtain the initial shear wave velocity value of the current rectangular grid element. The correspondence between the depth value and the wavelength is converted according to the principle that the wavelength is equal to 2.5 times the depth value.

[0013] Furthermore, in step three, the specific process of wavefield forward modeling is as follows: using the current velocity model as the medium parameter field, the soil surface position directly below the center of the vibrating wheel as the source position, and the arithmetic mean of the vertical component records of the left and right triaxial MEMS accelerometers as the source time function, the elastic wave equation under two-dimensional plane strain conditions is solved using a rotating staggered grid finite difference scheme. The rotating staggered grid finite difference scheme uses a fourth-order precision central difference operator in the spatial direction and a second-order precision central difference operator in the time direction. During the wavefield propagation calculation, the vertical particle velocity values ​​at each time step of the engineering seismic detector position of each vertical component are extracted to form a synthetic seismic waveform record, and the wavefield snapshots of all spatial grid nodes at each time step are stored in the cache area of ​​the edge computing gateway.

[0014] Furthermore, in step three, the specific process of the accompanying wavefield backpropagation is as follows: For each vertical component engineering seismic detector location, the synthetic seismic waveform record of the current vertical component engineering seismic detector location is subtracted from the measured seismic waveform record at each sampling point in the time domain to obtain the residual waveform. The residual waveform is then flipped along the time axis to obtain the time-reversed residual waveform. Using all vertical component engineering seismic detector locations as the simultaneously acting virtual source locations, and using the time-reversed residual waveform corresponding to each vertical component engineering seismic detector as the source time function of each virtual source, the elastic wave equation is solved in reverse from the end of the calculation time to the beginning using the same rotating staggered grid finite difference scheme as the wavefield forward modeling calculation to obtain the accompanying wavefield. The accompanying wavefield snapshot of all spatial grid nodes at each time step is stored in the cache area of ​​the edge computing gateway.

[0015] Furthermore, in step three, the specific process of calculating the sensitive core field is as follows: The forward wave field snapshot and the accompanying wave field snapshot are read sequentially at each time step. For each time step, each rectangular grid cell in the current inversion grid is traversed. The forward wave field horizontal particle velocity component and the forward wave field vertical particle velocity component of all spatial grid nodes within the coverage area of ​​the current rectangular grid cell are extracted, and the forward wave field shear strain rate component is calculated. The accompanying wave field shear stress component within the coverage area of ​​the current rectangular grid cell is extracted. The forward wave field shear strain rate component and the accompanying wave field shear stress component are multiplied node by node within the current rectangular grid cell and then summed to obtain the sensitive core contribution value of the current rectangular grid cell at the current time step. For each rectangular grid cell, the sensitive core contribution values ​​of all time steps are accumulated and multiplied by the time step length to obtain the shear wave velocity sensitive core value of the current rectangular grid cell. The shear wave velocity sensitive core values ​​of all rectangular grid cells constitute the sensitive core field.

[0016] Furthermore, in step three, the specific process of identifying highly sensitive anomaly cells and adaptively refining the mesh is as follows: For each non-boundary rectangular mesh cell in the current inverted mesh, a neighborhood range centered on the current non-boundary rectangular mesh cell is determined. The neighborhood range includes four adjacent rectangular mesh cells: the left and right rectangular mesh cells adjacent to the current non-boundary rectangular mesh cell in the lateral direction, and the upper and lower rectangular mesh cells adjacent in the depth direction. The shear wave velocity sensitive kernel value of the current non-boundary rectangular mesh cell and the shear wave velocity sensitive kernel values ​​of the four adjacent rectangular mesh cells are extracted, and their arithmetic mean is calculated as the neighborhood sensitive kernel mean. The absolute value of the difference between the shear wave velocity sensitive kernel value of the current non-boundary rectangular grid cell and the mean value of the neighboring sensitive kernels is used as the local anomaly index of the current non-boundary rectangular grid cell. The maximum value of the local anomaly index of all non-boundary rectangular grid cells is selected as the anomaly benchmark value. The anomaly benchmark value is multiplied by a coefficient of 0.55 to obtain the anomaly judgment threshold. Non-boundary rectangular grid cells with a local anomaly index greater than the anomaly judgment threshold are marked as high-sensitivity anomaly cells. Quadtree mesh refinement is performed on each high-sensitivity anomaly cell to form 4 sub-cells. The shear wave velocity values ​​of the 4 sub-cells all inherit the shear wave velocity values ​​of the corresponding high-sensitivity anomaly cells.

[0017] Furthermore, in step three, the specific process of velocity model iterative update is as follows: traverse each rectangular grid cell in the current inversion grid, divide the shear wave velocity sensitive kernel value of the current rectangular grid cell by the maximum absolute value of the shear wave velocity sensitive kernel values ​​of all rectangular grid cells to obtain the normalized sensitive kernel value, multiply the normalized sensitive kernel value by the current shear wave velocity value and then multiply by the fixed iteration step size coefficient of 0.03 to obtain the velocity update increment, subtract the velocity update increment from the current shear wave velocity value to obtain the updated shear wave velocity value; apply physical feasibility constraints to the updated shear wave velocity value, if the updated shear wave velocity value is less than 70 m / s, it is forcibly corrected to 70 m / s, if the updated shear wave velocity value is greater than 380 m / s, it is forcibly corrected to 380 m / s.

[0018] The present invention provides a real-time prediction and multi-dimensional intelligent evaluation method for compaction quality of soft soil subgrades, which has the following advantages: The invention creatively uses the periodic vibration of the vibratory roller as an active seismic excitation source, obtaining high-quality seismic wave signals without the need for a separate artificial seismic source device. The vibratory roller continuously generates stable vibration output during normal compaction operations, with good repeatability in excitation frequency and excitation force amplitude, providing ideal source conditions for seismic wave detection. This source utilization method allows compaction quality detection to be carried out simultaneously with the compaction construction process, eliminating the operational interruption problem of traditional detection methods requiring machine shutdown for sampling or separate seismic source deployment, significantly improving detection efficiency and achieving real-time monitoring of compaction status. The full-waveform iterative inversion method based on hierarchical progressive grid refinement and local sensitive core focusing proposed in this invention breaks through the resolution limitations of traditional uniform grid inversion. By analyzing the spatial coherence characteristics of the sensitive core field, it automatically identifies regions with significant local anomalies in the velocity model, performing quadtree-style grid refinement only on these regions. This achieves high spatial resolution in compaction inhomogeneity regions requiring fine characterization, while maintaining coarser resolution in regions with gradual velocity changes. This adaptive grid strategy optimizes the allocation of computational resources, effectively controlling the computational burden while ensuring inversion accuracy, enabling the full-waveform inversion method to meet the timeliness requirements of real-time on-site detection. This invention establishes a multi-dimensional compaction quality evaluation framework that comprehensively considers depth stratification and overall uniformity. By dividing the inversion area into shallow, middle, and deep layers according to depth, and calculating compaction coefficient estimates for each layer, it can accurately identify potential deep-layer compaction deficiencies, overcoming the limitations of traditional detection methods in obtaining deep-layer compaction information. Simultaneously, the compaction uniformity variation coefficient is introduced as an overall uniformity evaluation index, which can identify potential localized weak areas even when the average compaction coefficient meets the standard. The combined judgment rules of the multi-dimensional indicators achieve a comprehensive and objective evaluation of compaction quality. Attached Figure Description

[0019] Figure 1 This is a schematic diagram of the overall layout scheme of the self-excited seismic wave synchronous acquisition system for road rollers provided in an embodiment of the present invention;

[0020] Figure 2 The calculation results of the Rayleigh wave dispersion energy spectrum and the schematic diagram of the dispersion curve extracted from it are provided for the embodiments of the present invention;

[0021] Figure 3 This is a schematic diagram of multichannel seismic waveform recordings acquired by a detector array according to an embodiment of the present invention;

[0022] Figure 4 This is a schematic diagram of the convergence curve showing the change of residual energy with the number of iterations during the full waveform iterative inversion process provided in an embodiment of the present invention. Detailed Implementation

[0023] The method of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0024] A method for real-time prediction and multi-dimensional intelligent evaluation of compaction quality in soft soil subgrades, comprising:

[0025] Step 1: Multiple vertical component engineering seismic detectors are deployed on the soil surface behind the vibratory roller wheel along the direction of travel to form a detector array. An accelerometer is installed on the bearing seat of the vibratory roller wheel to record the vibration characteristics of the vibratory roller wheel. All sensors are connected to the edge computing gateway for synchronous data acquisition.

[0026] This invention utilizes the periodic vibration generated by the vibratory drum of a road roller during operation as an active seismic excitation source. By deploying an array of seismic detectors on the soil surface to receive seismic wave signals propagating through the subgrade medium, it enables the detection of the internal structural characteristics of soft soil subgrades. Compared with traditional methods using manual hammering or drop-weight devices as the seismic source, using the vibration of the road roller itself as the excitation source has significant advantages: the excitation frequency and excitation force amplitude of the vibratory drum are relatively stable, providing seismic source conditions with good repeatability; the excitation process is synchronized with the compaction work, and detection can be completed without additional machine shutdown; the coupling state between the vibratory drum and the soil directly reflects the current compaction conditions, and the acquired seismic wave response information is intrinsically correlated with the compaction quality.

[0027] The geophone array is arranged linearly, with 16 vertical component geophones sequentially installed on the soil surface behind the vibratory roller along the direction of the roller's travel. These vertical component geophones are specifically designed to receive the vertical vibrations of seismic waves. For Rayleigh waves, the vertical component amplitude dominates in shallow soil, thus using vertical component geophones effectively captures Rayleigh wave energy. The first vertical component geophone in the array is installed 1 meter from the rear edge of the vibratory roller. This distance is chosen by considering both near-field effects and signal attenuation. If the geophone is too close to the seismic source, the received signal will contain a strong near-field wave component. The propagation characteristics of near-field waves differ from those of far-field waves, interfering with subsequent dispersion curve extraction. If the distance is too far, the signal will attenuate significantly after propagating through the soil, reducing the signal-to-noise ratio. A 1-meter offset ensures the geophone is in the far-field region while maintaining sufficient signal strength. The spacing between adjacent vertical component seismic geophones was set to 0.3 meters, with 16 geophones forming a receiving array with a total length of 4.5 meters. The selection of geophone spacing is related to the shortest wavelength to be detected. According to the spatial sampling theorem, the geophone spacing should be less than half of the shortest wavelength to avoid spatial aliasing. The 0.3-meter spacing corresponds to a minimum resolvable wavelength of 0.6 meters, which translates to a Rayleigh wave phase velocity of approximately 60 meters per second multiplied by the reciprocal of the frequency, meeting the acquisition requirements for high-frequency surface wave components in soft soil subgrades.

[0028] refer to Figure 1 In this system, the vibratory roller of the road roller serves as the active seismic excitation source. The seismic waves generated by its periodic vibration propagate through the roadbed medium and are received by a geophone array deployed on the soil surface. The vibratory roller, represented by a circle, is located on the left side of the system, with its center approximately 0.6 meters above the ground surface. A triaxial microelectromechanical accelerometer is installed on the left and right sides of the vibratory roller's bearing housing to record the vibration characteristics of the roller. The symmetrical arrangement of the two sensors aims to eliminate any potential eccentric vibration components from the vibratory roller. The geophone array consists of 16 vertical component engineering seismic geophones, arranged linearly along the direction of the road roller's travel on the soil surface behind the vibratory roller. The first geophone is 1.0 meter horizontally from the rear edge of the vibratory roller, and the spacing between adjacent geophones is 0.3 meters. The 16 geophones form a receiving array with a total length of 4.5 meters. Each geophone is represented by an inverted triangle symbol, with its pointed conical base inserted below the soil surface to form a rigid coupling. The diagram uses concentric arcs to represent the wavefront of Rayleigh waves propagating outward from the vibrating wheel. The wavefront's energy gradually decreases as the propagation distance increases. The soil profile is divided into three layers based on depth: the shallow layer (0-0.35 meters), the middle layer (0.35-0.75 meters), and the deep layer (0.75-1.2 meters). Different layers are distinguished by different shades of tan. All sensors are connected via signal lines to an edge computing gateway installed in the roller's cab. The edge computing gateway integrates a multi-channel synchronous acquisition module, synchronously acquiring data from all 18 channels at a uniform sampling frequency of 4000 Hz. The zero-crossing rise time of the vertical component of the left-side triaxial MEMS accelerometer during each excitation cycle is used as the synchronous trigger reference for all channels.

[0029] In practical engineering applications, the method of fixing the geophone has a significant impact on signal quality. It is recommended to use an insertion-type geophone with a metal conical base, inserting the conical base completely 3 to 5 cm below the soil surface to ensure good rigid coupling between the geophone and the soil. If the soil surface is too soft, resulting in insufficient insertion depth, fine sand can be covered around the geophone and lightly compacted to improve coupling. The natural frequency of the geophone should be lower than the lowest frequency component of interest in the acquired signal. This invention recommends using an engineering seismic geophone with a natural frequency of 4.5 Hz, whose sensitivity maintains a flat response within the 10 to 200 Hz range, meeting the frequency band requirements for roadbed compaction monitoring.

[0030] The vibration characteristics of the vibratory roller are recorded using accelerometers mounted on the bearing housing. One triaxial MEMS accelerometer is installed on each of the left and right sides of the bearing housing; this symmetrical arrangement aims to eliminate any potential eccentric vibration components from the vibratory roller. The triaxial MEMS accelerometers can simultaneously measure acceleration components in three directions: horizontal, longitudinal, and vertical. The vertical component directly reflects the excitation characteristics of the vibratory roller on the soil. Compared to traditional piezoelectric accelerometers, MEMS accelerometers offer advantages such as smaller size, lower power consumption, and direct digital signal output, making them suitable for long-term stable operation in harsh engineering machinery environments. The sensor's range should be selected based on the vibration parameters of the vibratory roller. For conventional vibratory rollers, the vibration roller acceleration amplitude is typically in the range of 50 to 150 m / s², and a sensor with a range of ±200 m / s² is recommended to allow for a margin of error.

[0031] All sensors are connected to an edge computing gateway installed in the roller's cab via shielded cables. The shielded cables effectively suppress electromagnetic interference from the roller's engine and hydraulic system, ensuring signal transmission quality. The edge computing gateway integrates a multi-channel synchronous acquisition module, a data preprocessing module, and a wireless communication module. The synchronous acquisition module uses a uniform 4000 Hz sampling frequency for analog-to-digital conversion across all 18 channels. Based on the Nyquist sampling theorem, the 4000 Hz sampling frequency allows for distortion-free acquisition of signal components below 2000 Hz, fully covering the effective frequency band of roadbed seismic wave signals.

[0032] The implementation of multi-channel synchronous acquisition relies on a precise clock synchronization mechanism. This invention uses the zero-crossing rising edge of the vertical component signal from the left-side triaxial MEMS accelerometer as the synchronization trigger reference for all channels in each excitation cycle. Under normal operating conditions, the vibrating wheel reciprocates at a fixed frequency, and its vertical acceleration exhibits an approximately sinusoidal waveform, with each vibration cycle corresponding to a complete sinusoidal waveform. The zero-crossing rising edge refers to the moment when the acceleration value changes from negative to positive after crossing zero. This moment is unique and repeatable in each vibration cycle, making it suitable as a synchronization trigger reference. When the edge computing gateway detects the zero-crossing rising edge of the vertical component signal from the left-side triaxial MEMS accelerometer, it immediately triggers all 18 channels to simultaneously begin acquisition, continuously acquiring data for a complete excitation cycle. This synchronous triggering method ensures that each acquired data corresponds to the same vibration phase state of the vibrating wheel, facilitating the superposition and processing of data from multiple subsequent excitation cycles to improve the signal-to-noise ratio.

[0033] As an optional implementation, the number of detector arrays can be adjusted according to the actual detection depth requirements. If only shallow compaction within 0.5 meters needs to be detected, a shorter array of 8 detectors can be used; if detection of compaction at depths greater than 1.5 meters is required, the number of detectors can be increased to 24 or 32. The detector spacing can also be adjusted within the range of 0.2 to 0.5 meters according to the required detection resolution. Furthermore, the sampling frequency of the edge computing gateway can be selected within the range of 2000 to 8000 Hz based on signal characteristics. Higher sampling frequencies can acquire richer high-frequency information but increase the data storage and processing burden.

[0034] Step 2: Filter the original seismic waveform records acquired by the detector array, and use the phase shift transformation method to extract the Rayleigh wave dispersion curve from the filtered multi-channel waveform records.

[0035] The raw seismic waveform records acquired by the detector array need to undergo preprocessing to improve signal quality before entering the dispersion curve extraction process. The preprocessing process includes two steps: mean removal and bandpass filtering.

[0036] The purpose of mean removal processing is to eliminate any DC bias components that may exist in the original signal. DC bias may originate from factors such as sensor zero-point drift, analog-to-digital converter offset, or cable interference. If not eliminated, it will affect the accuracy of subsequent spectrum analysis. The specific operation of mean removal processing is as follows: for each waveform record, calculate the arithmetic mean of the amplitudes of all sampling points, and then subtract the calculated arithmetic mean from the amplitude of each sampling point to obtain waveform data with zero mean.

[0037] Bandpass filtering is used to preserve the frequency components of interest in the signal while suppressing noise interference. A fourth-order Butterworth filter is employed, whose amplitude-frequency response is maximally flat within the passband, avoiding amplitude distortion of the effective signal. The fourth-order design provides a roll-off rate of approximately 24 dB per octave, effectively suppressing out-of-band noise without introducing an excessively long filter impulse response length. The passband range is set from 5 to 200 Hz. The lower cutoff frequency of 5 Hz is used to suppress low-frequency mechanical vibration interference such as from roller movement and hydraulic pump operation, while the upper cutoff frequency of 200 Hz is used to suppress high-frequency random noise. This frequency range covers the main energy distribution range of Rayleigh waves in soft soil subgrades.

[0038] Rayleigh waves are surface waves that propagate along a free surface. Their propagation speed is slightly lower than the shear wave velocity of the medium, and a relatively stable conversion relationship exists between the two. A key characteristic of Rayleigh waves is their dispersion, meaning that Rayleigh wave components of different frequencies propagate with different phase velocities. In layered stratigraphic structures, the energy of high-frequency Rayleigh waves is mainly concentrated in the shallow layers, and their phase velocity is primarily controlled by the properties of the shallow medium. Low-frequency Rayleigh waves penetrate deeper, and their phase velocity reflects the comprehensive characteristics of the deeper layers. This correspondence between frequency and depth allows for the inversion of subsurface velocity structures by analyzing the phase velocities of Rayleigh waves at different frequencies. The dispersion curve represents the relationship between Rayleigh wave phase velocity and frequency.

[0039] This invention employs the phase shift transform method to extract the Rayleigh wave dispersion curve from 16 filtered waveform records. The basic principle of the phase shift transform method is as follows: if a Rayleigh wave of a certain frequency propagates at a specific phase velocity, there is a definite phase difference relationship between the frequency components received at different detector positions. This phase difference is related to the detector spacing and the phase velocity. By performing phase correction on each channel in the frequency domain, when the assumed experimental phase velocity matches the actual phase velocity, the corrected components of the same frequency in each channel will have the same phase, and their superposition will yield the maximum amplitude. When the assumed experimental phase velocity does not match the actual phase velocity, the corrected phases of each channel will be inconsistent, and their superposition will cancel each other out, resulting in a reduced amplitude. Therefore, by scanning different experimental phase velocities, the phase velocity that maximizes the superposition amplitude is the actual phase velocity of the Rayleigh wave at that frequency.

[0040] refer to Figure 2 The horizontal axis represents the analysis frequency, ranging from 5 to 100 Hz; the vertical axis represents the experimental phase velocity, ranging from 50 to 350 m / s. The dispersion energy spectrum is presented as a heatmap, with colors grading from white through yellow, orange to deep red. The darker the color, the greater the dispersion energy value at the combination of that frequency and phase velocity. A distinct high-energy band exists in the energy spectrum, corresponding to the dispersion characteristics of the Rayleigh fundamental mode. At each analysis frequency... At that point, the experimental phase velocity corresponding to the maximum energy value in the dispersion energy spectrum is selected as the Rayleigh wave phase velocity at that frequency. The Rayleigh wave phase velocities corresponding to all analyzed frequencies form a dispersion curve, plotted as a solid blue line on the energy spectrum. The dispersion curve exhibits typical normal dispersion characteristics, i.e., the phase velocity decreases with increasing frequency. This is because the deep soil in the soft soil subgrade has a higher degree of compaction due to prior compaction and natural consolidation, while the shallow newly filled soil has not been fully compacted. At approximately 10 Hz in the low-frequency range, the Rayleigh wave phase velocity is approximately 230 m / s, reflecting the velocity characteristics of the higher compaction of the deeper soil layers; at approximately 90 Hz in the high-frequency range, the Rayleigh wave phase velocity decreases to approximately 140 m / s, reflecting the velocity characteristics of the lower compaction of the shallow soil layers. Several sampling points are marked with black dots on the dispersion curve, representing the phase velocity values ​​extracted at discrete frequencies. The dispersion energy is calculated using the phase shift transformation method for each analyzed frequency. Phase velocity of the test The combination is based on the horizontal distance of each detector from the vibrating wheel. Calculate theoretical propagation time Converting theoretical propagation time into phase offset The dispersion energy value is obtained by performing phase correction on the complex spectrum of each waveform and then superimposing and taking the modulus. .

[0041] The specific implementation process of the phase shift transform method is as follows: First, a 2048-point Fast Fourier Transform is performed on each filtered waveform record to obtain a complex spectrum. The 2048-point transform length corresponds to a time window length of 0.512 seconds at a sampling frequency of 4000 Hz, with a frequency resolution of approximately 1.95 Hz, which meets the accuracy requirements for dispersion curve extraction. The Fast Fourier Transform converts the time-domain signal into a frequency-domain representation, and the transform result is in complex form, containing the amplitude and phase information of each frequency component.

[0042] In the dispersion energy calculation stage, it is necessary to traverse the two-dimensional parameter space composed of the analysis frequency and the test phase velocity. The scanning range of the analysis frequency is set to 5 to 100 Hz, with each analysis frequency selected sequentially at 0.5 Hz intervals, for a total of 191 frequency points. The frequency range of 5 to 100 Hz corresponds to the main energy distribution range of Rayleigh waves in soft soil subgrades. Components with frequencies that are too low are difficult to be effectively acquired by a detector array of finite length, while components with frequencies that are too high have weak energy due to soil absorption and attenuation. The scanning range of the test phase velocity is set to 50 to 350 m / s, with each test phase velocity selected sequentially at 2 m / s intervals, for a total of 151 phase velocity points. This phase velocity range covers the typical shear wave velocity range of soft soil from an uncompacted loose state to a fully compacted dense state.

[0043] refer to Figure 3The horizontal axis represents the recording time, ranging from 0 to 150 milliseconds; the vertical axis represents the detector number and its corresponding offset. The figure displays the waveforms of eight detectors (numbers 1, 3, 5, 7, 9, 11, 13, and 15) at intervals, with offsets increasing from 1.0 meter to 5.2 meters. Each waveform is plotted as a black curve, with the positive half-cycle filled in gray for enhanced visibility. The propagation characteristics of Rayleigh waves are clearly observed in the waveform records: the main vibration phase of each waveform gradually delays with increasing offset, reflecting the travel time difference required for the seismic wave to propagate from the epicenter to each detector location. The theoretical arrival time line of Rayleigh waves is marked with a red dashed line in the figure. This arrival time line is calculated based on an average phase velocity of approximately 180 meters per second, and the main energy envelope of each waveform matches the theoretical arrival time line well. The waveform recordings also reveal dispersion, where different frequency components propagate at varying speeds, causing stretching and distortion during propagation. High-frequency components propagate slower and arrive later, while low-frequency components propagate faster and arrive earlier. The waveform at longer offsets exhibits more pronounced dispersion broadening compared to that at shorter offsets. This is because the longer the propagation distance, the more significant the cumulative difference in travel time between different frequency components. The time window of the waveform recordings covers the arrival periods of the main energy components of the Rayleigh wave, providing a complete data foundation for subsequent dispersion curve extraction and full waveform inversion.

[0044] For each analysis frequency Phase velocity of the test The process of calculating the dispersion energy value based on the combination of factors is as follows: First, calculate the theoretical propagation time corresponding to each geophone based on the horizontal distance between the geophone and the vibrating wheel for each vertical component. Let the first geophone be... The horizontal distance between the vertical component seismic detector and the center of the vibrating wheel is: Then the theoretical propagation time corresponding to the detector is Calculate according to the following formula: ;in This represents the horizontal distance, in meters, from the center of the vibrating wheel to the i-th vertical component seismic detector. This indicates the current experimental phase velocity, expressed in meters per second. This represents the calculated theoretical propagation time, expressed in seconds.

[0045] The theoretical propagation time is then converted into a phase offset. For the analysis frequency... Theoretical dissemination time Corresponding phase offset Calculate according to the following formula: ;in This indicates the current analysis frequency, in Hertz (Hz). Indicates the first The theoretical propagation time corresponding to each vertical component engineering seismic detector, in seconds; This represents the calculated phase offset, expressed in radians.

[0046] Extract the current analysis frequency from the complex spectrum of each waveform record. Let the complex value at the position be . Waveform recordings at the analysis frequency The complex spectrum value at that point is The process of performing phase correction on it is as follows: convert the complex spectrum values Multiply by phase offset The corresponding complex exponential factor is obtained as a phase-corrected complex value. The calculation formula is: ;in Indicates the first Waveform recordings at the analysis frequency Complex spectral values ​​at; Indicates the first Phase offset corresponding to each vertical component of the engineering seismic detector; Represents the imaginary unit; This represents the complex exponential factor corresponding to the phase offset; This represents the complex value after phase correction.

[0047] After completing phase correction for all 16 waveform records, the complex values ​​after phase correction are summed to obtain a superimposed complex value, and its magnitude is taken as the current analysis frequency. Phase velocity of the test The corresponding dispersion energy value of the combination The calculation formula is: ;in Indicates the first The waveform records the complex values ​​after phase correction; the summation symbol indicates that the 16 data points are accumulated; the vertical bar symbol indicates that the modulus of the complex number is taken; This represents the dispersion energy value corresponding to the combination of the current analysis frequency and the experimental phase velocity.

[0048] By traversing all 191 analytical frequencies and 151 experimental phase velocities, a two-dimensional dispersion energy spectrum containing 28,841 dispersion energy values ​​was calculated. At each analytical frequency, all experimental phase velocities were traversed, and the experimental phase velocity corresponding to the maximum dispersion energy value was selected as the Rayleigh wave phase velocity for that current analytical frequency. The location of the maximum dispersion energy value corresponds to the situation where the experimental phase velocity and the true phase velocity are consistent; at this point, the complex values ​​after phase correction of each channel have the same phase, and the vector superposition reaches its maximum amplitude. The Rayleigh wave phase velocities corresponding to all 191 analytical frequencies constitute a Rayleigh wave dispersion curve, which describes the variation of the Rayleigh wave phase velocity with frequency.

[0049] In typical cases of soft soil subgrades, Rayleigh wave dispersion curves usually exhibit normal dispersion characteristics, with phase velocity decreasing as frequency increases. This is because the deep soil has undergone prior compaction and natural consolidation, resulting in a high degree of density, while the shallow newly filled soil has not yet been fully compacted. If the dispersion curve exhibits abnormal dispersion characteristics, with phase velocity increasing as frequency increases, it may indicate problems such as shallow soft interlayers or uneven compaction, requiring attention from construction quality management personnel.

[0050] As an optional implementation, the scanning range and interval of the analysis frequency and test phase velocity can be adjusted according to specific working conditions. For subgrades with high compaction, the upper limit of the test phase velocity can be increased to 400 m / s; for soft soils with high water content, the lower limit of the test phase velocity can be reduced to 30 m / s. Refining the scanning interval can improve the resolution accuracy of the dispersion curve, but it will correspondingly increase the computation time. In addition, when there are multiple maxima in the dispersion energy spectrum, it may correspond to the presence of higher-order mode Rayleigh waves. This invention prioritizes extracting the dispersion curve of the most energetic fundamental mode for subsequent inversion, and the dispersion information of higher-order modes can be used as auxiliary constraints in more advanced inversion algorithms.

[0051] Step 3: Establish an initial velocity model based on the Rayleigh wave dispersion curve. Through a cyclic iterative process of wavefield forward modeling, adjoint wavefield back propagation, sensitive core field calculation, identification of highly sensitive anomaly elements and adaptive mesh refinement, velocity model iterative update, and iterative convergence determination, the final shear wave velocity model is output.

[0052] Full waveform inversion is a method for reconstructing subsurface velocity structures by minimizing the differences between observed and synthesized seismic waveforms. Compared to traditional inversion methods that only utilize travel time information or dispersion curves, full waveform inversion fully leverages the amplitude, phase, and waveform details contained in the seismic waveform, enabling the acquisition of higher-resolution velocity models. This invention addresses the specific needs of compaction quality testing in soft soil subgrades by proposing a full waveform iterative inversion method based on hierarchical progressive grid refinement and local sensitive kernel focusing. By adaptively adjusting the inversion grid resolution, it achieves a fine characterization of unevenly compacted areas while ensuring computational efficiency.

[0053] Establishing an initial velocity model is a crucial prerequisite for successful full-waveform inversion. Due to the strong nonlinear characteristics of the full-waveform inversion problem, if the initial model differs significantly from the actual model, the iterative process is prone to getting trapped in local minima and failing to converge to the global optimum. This invention utilizes Rayleigh wave dispersion curves to construct an initial velocity model. These dispersion curves reflect the average velocity characteristics of soil layers at different depths, providing a reasonable initial estimate for the inversion.

[0054] The initial inversion mesh was set to a rectangular subgrade area directly below the vibratory roller, with a lateral width of 2.4 meters and a depth of 1.2 meters. The 2.4-meter lateral width is slightly larger than the width of a conventional vibratory roller's drum, covering the direct impact area of ​​the roller and the affected areas on both sides. The 1.2-meter depth range covers the depth range typically considered in subgrade compaction testing; for assessing the compaction state of deeper layers, the inversion depth can be extended as needed. The initial inversion mesh was uniformly divided into 8 columns of cells laterally and 6 rows of cells in the depth direction, forming 48 rectangular mesh cells. Each rectangular mesh cell has a lateral dimension of 0.3 meters and a depth dimension of 0.2 meters. This initial mesh size strikes a balance between computational efficiency and model resolution; subsequent adaptive refinement mechanisms gradually increase the resolution in areas requiring fine characterization.

[0055] The method for setting the initial shear wave velocity value for each rectangular grid cell based on the Rayleigh wave dispersion curve is as follows: Calculate the depth value at the center point of each rectangular grid cell, and find the Rayleigh wave phase velocity value from the Rayleigh wave dispersion curve that matches the wavelength corresponding to the depth value at the center point of the current rectangular grid cell. The penetration depth of the Rayleigh wave is closely related to its wavelength. Empirical studies have shown that the Rayleigh wave energy is mainly concentrated in the depth range of approximately one-third to one-half of the wavelength. This invention adopts the correspondence between depth value and wavelength by converting the wavelength to 2.5 times the depth value. Let the depth value at the center point of the rectangular grid cell be... The corresponding wavelength Calculate according to the following formula: ;in This represents the depth value of the center point of a rectangular grid cell, in meters. This represents the Rayleigh wave wavelength corresponding to the depth value, in meters.

[0056] According to wavelength The relationship between Rayleigh wave propagation speed and frequency, and the corresponding frequency. Calculate according to the following formula: ;in This represents the phase velocity of Rayleigh waves, measured in meters per second. This represents the Rayleigh wave wavelength, measured in meters. This represents the corresponding frequency, measured in Hertz (Hz). Since the dispersion curve describes the relationship between phase velocity and frequency, the depth value is found from the dispersion curve using iterative interpolation. Matched Rayleigh wave phase velocity values .

[0057] Multiplying the current Rayleigh wave phase velocity value by a coefficient of 0.92 yields the initial shear wave velocity value for the current rectangular mesh element. The calculation formula is: ;in This represents the Rayleigh wave phase velocity value obtained from the dispersion curve, in meters per second; 0.92 is the conversion factor between Rayleigh wave phase velocity and shear wave velocity. This represents the initial shear wave velocity value, in meters per second. The conversion factor between Rayleigh wave phase velocity and shear wave velocity depends on the Poisson's ratio of the medium. For typical soil materials with a Poisson's ratio in the range of 0.25 to 0.35, the conversion factor varies between 0.90 and 0.94. This invention takes 0.92 as the representative conversion factor for soft soil subgrades.

[0058] After the initial velocity model is established, a hierarchical and progressive full-waveform iterative inversion process is initiated. The iterative inversion at each grid level includes six sub-steps: forward wavefield modeling, backward propagation of the adjoint wavefield, sensitive core field calculation, identification of highly sensitive anomalous elements and adaptive mesh refinement, iterative update of the velocity model, and iterative convergence determination. Each sub-step is executed cyclically until the convergence condition is met or the maximum number of iterations is reached.

[0059] The purpose of wavefield forward modeling is to simulate the propagation of seismic waves in a medium under given velocity model conditions and to calculate the synthetic seismic waveform records at each detector location. The forward modeling uses the current velocity model as the medium parameter field and the location of the soil surface directly below the center of the vibrating wheel as the source location. The source time function is the arithmetic mean of the vertical component records from the left and right triaxial MEMS accelerometers. Averaging the records from both sensors eliminates the influence of possible eccentric vibration components from the vibrating wheel, resulting in a purer vertical excitation signal.

[0060] Forward modeling employs a rotated staggered mesh finite difference scheme to solve the elastic wave equations under two-dimensional plane strain conditions. Compared to traditional staggered meshes, the rotated staggered mesh exhibits better numerical stability when handling interfaces with strong velocity contrasts, reducing numerical dispersion and spurious reflections caused by mesh discretization. In the rotated staggered mesh finite difference scheme, stress and velocity components are defined at different mesh node locations, achieving second-order accuracy time progression through spatial staggering. The rotated staggered mesh finite difference scheme uses a fourth-order accurate central difference operator in the spatial direction and a second-order accurate central difference operator in the temporal direction. The fourth-order spatial difference operator has a smaller numerical dispersion error than the second-order difference operator, enabling more accurate wave field simulation results under the same mesh size.

[0061] The spatial step size of the computational grid is set to one-quarter of the smallest element size of the current inversion grid. When solving the wave equation using the finite difference method, the spatial step size needs to be sufficiently small to accurately represent the spatial variation characteristics of the wave field; typically, at least 8 to 10 grid nodes are required within each shortest wavelength range. Setting the spatial step size to one-quarter of the smallest element size of the inversion grid ensures both the accuracy of the wave field calculation and an integer multiple relationship with the inversion grid, facilitating spatial integration of the sensitive core field. The time step size needs to satisfy the Courant-Friedrich-Levy stability conditions. In this invention, the time step size is set as the spatial step size divided by the maximum shear wave velocity value in the current velocity model, multiplied by a coefficient of 0.3. This coefficient of 0.3 provides a safety margin while satisfying the stability conditions. Let the spatial step size be... The maximum shear wave velocity is Then the time step Calculate according to the following formula: ;in This indicates the spatial step size of the computational grid, in meters. This represents the maximum shear wave velocity value of all rectangular grid elements in the current velocity model, in meters per second. This indicates the time step of the calculation, in seconds.

[0062] The computation time was set to 0.5 seconds, which is sufficient to allow seismic waves to propagate from the hypocenter to all geophone locations and include the main reflected and surface wave components. During the wavefield propagation calculation, the vertical particle velocity values ​​at each time step of the engineering seismic geophone locations were extracted to form a synthetic seismic waveform record. Simultaneously, a wavefield snapshot of all spatial grid nodes at each time step was stored in the cache area of ​​the edge computing gateway. The wavefield snapshot includes the particle velocity and stress components of each node and is used for subsequent sensitive core field calculations.

[0063] The adjoint wave field backpropagation is the core step in calculating gradient information during full waveform inversion. The mathematical basis of the adjoint method comes from the variational principle and the Lagrange multiplier method. By introducing the adjoint wave field, the gradient calculation of the objective function with respect to the model parameters is transformed into the cross-correlation operation of two wave fields, thereby reducing the computational complexity of gradient calculation from being proportional to the number of model parameters to requiring only two wave field propagation calculations.

[0064] The specific process of the accompanying wavefield backpropagation is as follows: For each vertical component engineering seismic detector location, the synthetic seismic waveform record at the current vertical component engineering seismic detector location is subtracted from the measured seismic waveform record at each sampling point in the time domain to obtain the residual waveform. The residual waveform reflects the difference between the synthetic waveform and the measured waveform under the current velocity model conditions, and these differences are precisely the targets that need to be eliminated by adjusting the velocity model. The residual waveform is flipped along the time axis to obtain the time-reversed residual waveform. The time-reversal operation allows the accompanying wavefield to propagate in the opposite direction to the forward wavefield, thereby achieving the meeting and cross-correlation of the two wavefields in the spatiotemporal domain.

[0065] Using the locations of all 16 vertical component engineering seismic detectors as the simultaneously acting virtual source locations, and the time-reversed residual waveforms corresponding to each vertical component engineering seismic detector as the source time function of each virtual source, the elastic wave equation is solved backwards from the end of the computation time to the beginning using the same rotating staggered grid finite difference scheme as the wavefield forward modeling to obtain the adjoint wavefield. Backpropagation starts from the end of the time frame because the time-reversed source function needs to be excited from the end of the original record, allowing the adjoint wavefield to be cross-correlated with the stored forward wavefield snapshot at the same physical moment. During the adjoint wavefield propagation calculation, the adjoint wavefield snapshots of all spatial grid nodes at each time step are stored in the cache area of ​​the edge computing gateway.

[0066] Sensitive kernel field calculations reveal the degree of sensitivity of the objective function to velocity parameters at various spatial locations. Regions with large sensitive kernel values ​​have a greater impact on waveform fitting residuals, and prioritizing the adjustment of velocity values ​​in these regions can more effectively reduce residuals. The specific process of sensitive kernel field calculation is as follows: sequentially read the forward wavefield snapshot and the accompanying wavefield snapshot at each time step, and perform sensitive kernel contribution value calculation for each rectangular grid cell in the current inversion grid for each time step.

[0067] For the current rectangular grid cell, extract the positive wave field horizontal particle velocity components of all spatial grid nodes within the coverage area of ​​the current rectangular grid cell. and the vertical particle velocity component of the positive wave field The shear wave velocity sensitive nucleus is related to the shear strain rate, and the shear strain rate component of the positive wave field... Calculate according to the following formula: ;in This represents the particle velocity component in the horizontal direction of the positive wave field, with units of meters per second. This represents the vertical particle velocity component of the positive wave field, measured in meters per second. This represents the partial derivative of the horizontal velocity component of a particle with respect to the depth direction. This represents the partial derivative of the vertical particle velocity component with respect to the horizontal direction; This represents the shear strain rate component of the positive wave field, measured in seconds.

[0068] Extract the shear stress components of the associated wave field within the coverage area of ​​the current rectangular grid cell. The sensitive core contribution value of the current rectangular grid element at the current time step is obtained by multiplying the shear strain rate component of the positive wave field and the shear stress component of the accompanying wave field node by node within the current rectangular grid element and then summing them. : ;in This represents the shear strain rate component of the positive wave field; This represents the shear stress component of the accompanying wave field, in Pascals. This represents the summation of all spatial grid nodes within the coverage area of ​​the current rectangular grid cell; This represents the sensitive kernel contribution value of the current rectangular mesh cell at the current time step.

[0069] For each rectangular mesh cell, the sensitivity kernel contribution values ​​of all time steps are summed and multiplied by the time step size to obtain the shear wave velocity sensitivity kernel value of the current rectangular mesh cell. : ;in Indicates the time step, in seconds; This represents the summation over all time steps; This represents the sensitive kernel contribution value at each time step; This represents the shear wave velocity sensitive core value of the current rectangular grid cell. The shear wave velocity sensitive core values ​​of all rectangular grid cells constitute the sensitive core field.

[0070] The key innovations of this invention, distinguishing it from traditional full-waveform inversion methods, lie in the identification of highly sensitive anomaly cells and adaptive mesh refinement. Traditional full-waveform inversion employs uniform mesh partitioning, using the same mesh resolution across the entire inversion region. This approach incurs unnecessary computational overhead in regions with gentle velocity changes, while the resolution may be insufficient in regions with drastic velocity variations. This invention identifies regions in the velocity model requiring fine characterization by analyzing the spatial coherence characteristics of sensitive core fields, and refines the mesh only in these regions, thus optimizing the allocation of computational resources.

[0071] The basis for performing local spatial coherence analysis on the sensitive core field is that in regions where the velocity model closely approximates the true model, the sensitive core field exhibits smooth and continuous spatial variation characteristics; while in regions where the velocity model has significant deviations or where there are uncaptured velocity anomalies, the sensitive core field shows local anomalies, manifested as significant differences in sensitive core values ​​between adjacent grid cells. By identifying these local anomaly regions and refining the mesh, the true velocity structure details can be gradually approximated in subsequent iterations.

[0072] For each non-boundary rectangular grid cell in the current inversion grid, determine the neighborhood range centered on the current non-boundary rectangular grid cell. Boundary cells are excluded from the analysis range due to the lack of complete neighborhood information. The neighborhood range includes four adjacent rectangular grid cells: the left and right rectangular grid cells adjacent to the current non-boundary rectangular grid cell in the lateral direction, and the upper and lower rectangular grid cells adjacent in the depth direction. Extract the shear wave velocity sensitive kernel value of the current non-boundary rectangular grid cell. and the shear wave velocity sensitive core value of 4 adjacent rectangular grid cells , , , The arithmetic mean of the five sensitive kernel values ​​is calculated as the mean of the neighborhood sensitive kernel. : ;in This represents the shear wave velocity sensitive kernel value of the current non-boundary rectangular mesh cell; , , , These represent the shear wave velocity sensitive core values ​​of four adjacent rectangular grid cells; This represents the mean of the neighborhood sensitive kernels.

[0073] The absolute value of the difference between the shear wave velocity sensitive kernel value of the current non-boundary rectangular mesh cell and the mean value of the neighboring sensitive kernel is calculated as the local anomaly index of the current non-boundary rectangular mesh cell. : ;in This represents the shear wave velocity sensitive kernel value of the current non-boundary rectangular mesh cell; This represents the mean of the neighborhood sensitive kernels; This represents the local anomaly index. The local anomaly index reflects the degree of deviation between the current cell's sensitivity kernel value and the average level of its neighborhood. The larger the value, the more significant the local features of the velocity model in that region, which require more detailed characterization.

[0074] The maximum value among all local anomaly indices of non-boundary rectangular grid cells is selected as the anomaly baseline value. The anomaly threshold is obtained by multiplying the anomaly baseline value by a coefficient of 0.55. : ;in This represents the maximum value of the local anomaly index for all non-boundary rectangular grid cells; 0.55 is the anomaly determination coefficient. This represents the threshold for anomaly determination. The coefficient 0.55 was selected through numerical experiments, which verified that it can identify the main anomaly areas while avoiding a significant increase in computational burden caused by excessive refinement.

[0075] Non-boundary rectangular mesh elements with local anomaly indices exceeding the anomaly threshold are marked as highly sensitive anomaly elements. Quadtree mesh refinement is performed on each highly sensitive anomaly element, dividing it horizontally into left and right halves and vertically into upper and lower halves, forming four sub-elements. The shear wave velocity values ​​of these four sub-elements inherit those of the corresponding highly sensitive anomaly element. Rectangular mesh elements not marked as highly sensitive anomaly elements retain their original size and velocity values. Quadtree mesh refinement results in a non-uniform structure in the inverted mesh, providing higher resolution in areas requiring fine detail and maintaining coarser resolution in areas with gradual velocity changes, thus achieving a balance between computational efficiency and model accuracy.

[0076] The velocity model iteratively updates the shear wave velocity values ​​of each rectangular grid cell based on the gradient direction indicated by the sensitive core field. It iterates through each rectangular grid cell in the current inversion grid, dividing the shear wave velocity sensitive core value of the current rectangular grid cell by the maximum absolute value of the shear wave velocity sensitive core values ​​of all rectangular grid cells to obtain the normalized sensitive core value. : ;in This represents the shear wave velocity sensitive kernel value of the current rectangular mesh cell; This represents the maximum absolute value of the shear wave velocity sensitive kernel value for all rectangular grid cells; This represents the normalized sensitivity kernel value. The normalization operation maps the sensitivity kernel value to the range of -1 to +1, eliminating the influence of the magnitude difference in sensitivity kernel values ​​between different iterations, so that the fixed iteration step size coefficient can maintain a stable update amplitude throughout the inversion process.

[0077] Multiply the normalized sensitive kernel value by the current shear wave velocity value. Multiply by the fixed iteration step size coefficient of 0.03 to obtain the speed update increment. : ;in Indicates the normalized sensitive kernel value; This represents the current shear wave velocity value, in meters per second; 0.03 is the fixed iteration step size coefficient. This represents the velocity update increment, in meters per second. The updated shear wave velocity value is obtained by subtracting the velocity update increment from the current shear wave velocity value. : ;in This indicates the current shear wave velocity value, in meters per second. This indicates the speed update increment, in meters per second; This represents the updated shear wave velocity value, in meters per second.

[0078] Physical feasibility constraints are applied to the updated shear wave velocity values ​​to ensure that the inversion results fall within the reasonable velocity range for soft soil subgrade materials. If the updated shear wave velocity value is less than 70 m / s, it is forcibly corrected to 70 m / s, which corresponds to the lower limit of the typical shear wave velocity for extremely soft soil with high water content. If the updated shear wave velocity value is greater than 380 m / s, it is forcibly corrected to 380 m / s, which corresponds to the upper limit of the typical shear wave velocity for fully compacted soil. Physical feasibility constraints can prevent extreme velocity values ​​that do not conform to reality from appearing during the inversion process, thereby improving the stability and reliability of the inversion results.

[0079] Iterative convergence determination is used to determine whether the iteration at the current grid level has converged and whether it is necessary to proceed to the next grid level for refinement inversion. The method for calculating the waveform fitting residual energy of the current iteration is as follows: calculate the sum of squares of each sampling point of the residual waveform for each of the 16 synthetic and 16 measured seismic waveform records, and then sum the sums of squares of the 16 residual waveforms to obtain the total residual energy. Let the first... The first residual waveform The value of each sampling point is The total residual energy is calculated according to the following formula: ;in Indicates the first The first residual waveform Values ​​at each sampling point; This represents the summation of the 16 residual waveforms; This represents the summation of all sampling points of each residual waveform; This represents the total residual energy.

[0080] Compare the total residual energy of the current iteration with the total residual energy of the previous iteration, and calculate the percentage decrease in residual energy. : ;in This represents the total residual energy from the previous iteration; This represents the total residual energy in the current iteration; This indicates the percentage decrease in residual energy.

[0081] If the percentage decrease in residual energy is less than 1.5%, the current grid level is considered converged. A decrease in the magnitude of the residual energy decrease indicates that the velocity model is close to optimal at the current grid resolution, and further iterations are unlikely to yield significant improvement. After convergence, the lateral dimensions of all rectangular grid cells in the current inversion grid are checked, and the minimum value is selected as the current minimum cell size. If the current minimum cell size is greater than 0.04 meters, the iteration continues at the next grid level. In the next grid level, high-sensitivity anomaly cell identification and grid refinement are re-executed to further improve the resolution of anomaly regions. If the current minimum cell size is less than or equal to 0.04 meters, the inversion at all levels is considered complete. The 0.04-meter minimum size limit ensures that the inversion resolution does not exceed the accuracy required for actual engineering, avoiding excessive refinement that leads to wasted computational resources and numerical instability. The current velocity model is then output as the final shear wave velocity model.

[0082] If the residual energy decreases by a percentage greater than or equal to 1.5%, the current grid level is considered not yet converged, and the velocity model still has room for improvement. The iteration count counter is incremented by 1. If the iteration count counter does not exceed 30, the process returns to the wavefield forward modeling step to continue the next iteration of the current grid level. If the iteration count counter has exceeded 30, the current grid level is forcibly considered converged, and the converged processing flow is executed. The maximum number of iterations of 30 prevents excessive computation time due to slow convergence or getting trapped in local minima.

[0083] As an optional implementation, the iteration step size coefficient can be adaptively adjusted according to the decreasing trend of residual energy. When the percentage decrease in residual energy remains at a high level after multiple iterations, the step size coefficient can be appropriately increased to accelerate the convergence speed; when the residual energy oscillates or the decrease slows down, the step size coefficient can be appropriately decreased to improve convergence stability. Furthermore, the anomaly determination coefficient can be adjusted within the range of 0.4 to 0.7 according to the characteristics of the roadbed material and the required detection accuracy. Smaller coefficient values ​​will identify more highly sensitive anomaly units, resulting in a finer mesh but increasing the computational burden, while larger coefficient values ​​have the opposite effect.

[0084] refer to Figure 4 The horizontal axis represents the number of iterations, ranging from 1 to 30; the vertical axis uses a logarithmic scale to represent the percentage of normalized residual energy. Residual energy is defined as the sum of the squares of the differences between all 16 synthetic seismic waveforms and the measured seismic waveforms, calculated using the following formula: ,in Indicates the first The first residual waveform The sampling point values ​​are shown. The convergence curve is plotted with a solid blue line connecting the circular marker points, showing an exponential decay trend, indicating that the velocity model gradually approaches the true model during iteration. In the initial iteration phase, the residual energy decreases rapidly, from 100% to approximately 35% in the first 5 iterations; in the middle iteration phase, the residual energy continues to decrease steadily, reaching approximately 18% by the 15th iteration; in the later iteration phase, the decrease in residual energy flattens and gradually converges, eventually stabilizing at approximately 3.5%. The convergence threshold is indicated by a red horizontal dashed line in the figure, representing the percentage decrease in residual energy between two consecutive iterations. When the residual energy is less than 1.5%, the current mesh level is considered to have converged. The figure also uses green vertical dashed lines to mark the three mesh refinement points, occurring after the 8th, 16th, and 24th iterations, respectively. Each mesh refinement corresponds to the turning point after the current level has converged, leading to the next level of refinement and inversion. The light blue filled area below the curve visually shows the cumulative decrease in residual energy, reflecting the effectiveness and convergence stability of the full waveform inversion algorithm.

[0085] Step 4: Extract the estimated compaction coefficient and compaction uniformity variation coefficient of each depth layer based on the final shear wave velocity model, and output the compaction quality grade according to the preset judgment rules.

[0086] The final shear wave velocity model reflects the shear wave velocity distribution at various spatial locations within the subgrade, and the shear wave velocity is directly related to the compaction state of the soil. During compaction, the porosity of the soil decreases, the particle contact area increases, and the skeleton stiffness improves, macroscopically manifested as an increase in the shear modulus. Since the shear wave velocity is proportional to the square root of the shear modulus, the higher the degree of compaction, the greater the shear wave velocity. Based on this physical relationship, this invention extracts multi-dimensional evaluation indicators from the final shear wave velocity model to comprehensively determine the compaction quality grade of soft soil subgrades.

[0087] Based on the depth values ​​of the center points of each rectangular grid cell in the final shear wave velocity model, all rectangular grid cells are divided into three depth layers: shallow, intermediate, and deep. The necessity of layered evaluation lies in the fact that roadbed compaction employs a layered filling and compaction process, and the number of compaction passes and compaction effects may differ at different depths. Furthermore, the vibration load from the road roller decreases with depth in the soil, with shallow soil experiencing stronger compaction than deeper soil. Therefore, different depth layers need to be evaluated separately to accurately identify potential weak points.

[0088] The specific method for dividing the depth layers is as follows: rectangular grid cells with a depth value in the range of 0 to 0.35 meters are classified as shallow layers; rectangular grid cells with a depth value in the range of 0.35 to 0.75 meters are classified as intermediate layers; and rectangular grid cells with a depth value in the range of 0.75 to 1.2 meters are classified as deep layers. The layer boundaries of 0.35 meters and 0.75 meters match the loose paving thickness commonly used in roadbed filling construction, and each depth layer roughly corresponds to the compaction state of one to two filling layers.

[0089] For each depth layer, calculate the arithmetic mean of the shear wave velocity values ​​of all rectangular mesh elements within that depth layer as the average shear wave velocity for that depth layer. Suppose a certain depth layer contains... There are rectangular grid cells, and the shear wave velocity values ​​of each cell are as follows: , ... The average shear wave velocity at the current depth level Calculate according to the following formula: ;in This indicates the number of rectangular grid cells within the current depth level; Indicates the first Shear wave velocity values ​​for each rectangular grid cell, in meters per second; This represents the average shear wave velocity at the current depth, expressed in meters per second.

[0090] The estimated compaction coefficient for each depth layer is obtained by referring to a pre-calibrated table of the relationship between shear wave velocity and compaction coefficient. This table is established through field calibration tests. Before practical engineering application, comparative tests of compaction coefficient and shear wave velocity measurements on representative soil samples are necessary to establish a conversion relationship applicable to the current engineering conditions. For typical soft soil subgrade materials, a shear wave velocity ranging from 120 to 180 meters per second corresponds to a compaction coefficient ranging from 0.90 to 0.96. The specific conversion relationship varies depending on factors such as soil type and moisture content.

[0091] Compaction uniformity is another important dimension for evaluating the compaction quality of roadbeds. Even if the average compaction coefficient at each depth layer meets the requirements, the existence of areas with insufficient compaction can still lead to uneven settlement or localized defects in the roadbed. The coefficient of variation of compaction uniformity reflects the dispersion of shear wave velocity distribution throughout the entire inversion region. It is calculated by determining the standard deviation of the shear wave velocity values ​​of all rectangular grid elements in the final shear wave velocity model. Divide by the arithmetic mean of the shear wave velocity values ​​of all rectangular grid cells. Obtain the coefficient of variation of compaction uniformity .

[0092] Assume the final shear wave velocity model includes There are rectangular grid cells, and the shear wave velocity values ​​of each cell are as follows: , ... The arithmetic mean of all element shear wave velocity values Calculate according to the following formula: ;in This represents the total number of rectangular mesh elements in the final shear wave velocity model; Indicates the first Shear wave velocity values ​​for each rectangular grid cell, in meters per second; This represents the arithmetic mean of the shear wave velocity values ​​for all rectangular grid cells, expressed in meters per second. Calculate according to the following formula:

[0093] ;

[0094] in This represents the total number of rectangular mesh elements in the final shear wave velocity model; Indicates the first Shear wave velocity values ​​for each rectangular grid cell, in meters per second; This represents the arithmetic mean of the shear wave velocities of all rectangular grid cells, in meters per second. The standard deviation of the shear wave velocity values ​​for all rectangular grid cells, in meters per second. Coefficient of variation for compaction uniformity. Calculate according to the following formula: ;in The standard deviation of the shear wave velocity values ​​of all rectangular grid cells is expressed in meters per second. This represents the arithmetic mean of the shear wave velocities of all rectangular grid cells, in meters per second. The coefficient of variation for compaction uniformity is a dimensionless quantity. A smaller coefficient of variation indicates a more uniform distribution of shear wave velocity within the subgrade and better consistency in the compacted material quality.

[0095] The determination of compaction quality grade comprehensively considers two indicators: the estimated compaction coefficient at each depth and the coefficient of variation of compaction uniformity. The determination rule employs a multi-condition combination approach for layer-level assessment.

[0096] If the estimated compaction coefficients of the shallow, middle, and deep layers are all greater than or equal to 0.96, and the coefficient of variation for compaction uniformity is less than or equal to 0.07, then the compaction quality grade is determined to be excellent. An excellent grade indicates that all layers at each depth of the subgrade have achieved a high degree of compaction, and the overall compaction state is uniform and consistent, meeting the compaction quality requirements of high-grade highways or important projects.

[0097] Otherwise, if the estimated compaction coefficients of the shallow, middle, and deep layers are all greater than or equal to 0.93, and the coefficient of variation for compaction uniformity is less than or equal to 0.11, then the compaction quality grade is determined to be good. A good grade indicates that the overall compaction quality of the subgrade is relatively good, meeting the compaction quality requirements of general highway engineering, with slight room for improvement that does not affect normal use.

[0098] Otherwise, if the estimated compaction coefficients of the shallow, middle, and deep layers are all greater than or equal to 0.90, and the coefficient of variation for compaction uniformity is less than or equal to 0.15, then the compaction quality grade is determined to be qualified. A qualified grade indicates that the subgrade compaction quality meets the basic requirements and can satisfy the compaction quality requirements of low-grade highways or temporary projects. It is recommended to appropriately increase the number of compaction passes to further improve the compaction effect.

[0099] Otherwise, the compaction quality grade is determined to be unqualified. An unqualified grade indicates that the subgrade compaction quality has not met the basic requirements, and there are problems such as excessively low compaction coefficients in one or more depth layers or poor overall compaction uniformity. It is necessary to carry out additional compaction treatment and retest until the qualified standard is met.

[0100] The compaction quality grade, estimated compaction coefficients for shallow, intermediate, and deep layers, and the coefficient of variation for compaction uniformity are displayed in real-time on the human-machine interface screen inside the roller's cab. The HMI screen uses an intuitive graphical interface, with different colors indicating the compaction quality grade for quick driver identification. All numerical indicators are displayed precisely in digital form for detailed analysis by technicians. When the compaction quality grade is unacceptable, the HMI screen simultaneously displays a prompt to guide the driver in performing additional compaction.

[0101] As an optional implementation method, the threshold parameters for determining the compaction quality grade can be adjusted according to the project category and technical specifications. For high-grade projects such as highways, the compaction coefficient threshold can be increased and the coefficient of variation threshold can be decreased to implement stricter quality control; for low-grade projects such as rural roads, the threshold requirements can be appropriately relaxed to improve construction efficiency. In addition to the compaction quality grade, spatial location information of weak compaction areas can also be output. By marking the locations of rectangular grid cells with shear wave velocity values ​​below the set threshold in the final shear wave velocity model on the roadbed plan or profile, construction personnel can be guided to perform targeted compaction in specific areas.

[0102] While specific embodiments of the present invention have been described above, those skilled in the art should understand that these specific embodiments are merely illustrative. Those skilled in the art can omit, substitute, and modify the details of the above methods and systems in various ways without departing from the principles and essence of the present invention. For example, combining the above method steps to perform substantially the same function and achieve substantially the same result according to substantially the same method falls within the scope of the present invention. Therefore, the scope of the present invention is defined only by the appended claims.

Claims

1. A method for real-time prediction and multi-dimensional intelligent evaluation of compaction quality for soft soil subgrades, characterized in that, The method includes: Step 1: Multiple vertical component engineering seismic detectors are deployed on the soil surface behind the vibratory roller wheel along the direction of travel to form a detector array. An accelerometer is installed on the bearing seat of the vibratory roller wheel to record the vibration characteristics of the vibratory roller wheel. All sensors are connected to the edge computing gateway for synchronous data acquisition. Step 2: Filter the raw seismic waveform records acquired by the detector array, and use the phase shift transformation method to extract the Rayleigh wave dispersion curve from the filtered multi-channel waveform records; Step 3: Establish an initial velocity model based on the Rayleigh wave dispersion curve. Through a cyclic iterative process of wavefield forward modeling, adjoint wavefield back propagation, sensitive core field calculation, identification of highly sensitive anomaly elements and adaptive mesh refinement, velocity model iterative update, and iterative convergence determination, the final shear wave velocity model is output. Step 4: Extract the compaction coefficient estimate and compaction uniformity variation coefficient of each depth layer based on the final shear wave velocity model, and output the compaction quality grade according to the preset judgment rules; In step three, the specific process of wavefield forward modeling is as follows: using the current velocity model as the medium parameter field, the position of the soil surface directly below the center of the vibrating wheel as the source position, and the arithmetic mean of the vertical component records of the left and right triaxial MEMS accelerometers as the source time function, the elastic wave equation under two-dimensional plane strain conditions is solved using the rotating staggered grid finite difference scheme. The rotating staggered grid finite difference scheme uses a fourth-order precision central difference operator in the spatial direction and a second-order precision central difference operator in the time direction. During the wavefield propagation calculation, the vertical particle velocity values ​​at each time step of the engineering seismic detector position of each vertical component are extracted to form a synthetic seismic waveform record, and the wavefield snapshots of all spatial grid nodes at each time step are stored in the cache area of ​​the edge computing gateway. In step three, the specific process of the accompanying wavefield back propagation is as follows: For each vertical component engineering seismic detector location, the synthetic seismic waveform record of the current vertical component engineering seismic detector location is subtracted from the measured seismic waveform record at each sampling point in the time domain to obtain the residual waveform. The residual waveform is then flipped along the time axis to obtain the time-reversed residual waveform. Using all vertical component engineering seismic detector locations as the simultaneously acting virtual source locations, and using the time-reversed residual waveform corresponding to each vertical component engineering seismic detector as the source time function of each virtual source, the elastic wave equation is solved in reverse from the end of the calculation time to the beginning using the same rotating staggered grid finite difference scheme as the wavefield forward modeling calculation to obtain the accompanying wavefield. The accompanying wavefield snapshot of all spatial grid nodes at each time step is stored in the cache area of ​​the edge computing gateway. In step three, the specific process of calculating the sensitive core field is as follows: The forward wave field snapshot and the accompanying wave field snapshot are read sequentially at each time step. For each time step, each rectangular grid cell in the current inversion grid is traversed. The forward wave field horizontal particle velocity component and the forward wave field vertical particle velocity component of all spatial grid nodes within the coverage area of ​​the current rectangular grid cell are extracted, and the forward wave field shear strain rate component is calculated. The accompanying wave field shear stress component within the coverage area of ​​the current rectangular grid cell is extracted. The forward wave field shear strain rate component and the accompanying wave field shear stress component are multiplied node by node within the current rectangular grid cell and then summed to obtain the sensitive core contribution value of the current rectangular grid cell at the current time step. For each rectangular grid cell, the sensitive core contribution values ​​of all time steps are accumulated and multiplied by the time step length to obtain the shear wave velocity sensitive core value of the current rectangular grid cell. The shear wave velocity sensitive core values ​​of all rectangular grid cells constitute the sensitive core field.

2. The method according to claim 1, characterized in that, In step one, the detector array contains 16 vertical component engineering seismic detectors. The first vertical component engineering seismic detector is 1 meter away from the rear edge of the vibrating wheel, and the spacing between adjacent vertical component engineering seismic detectors is 0.3 meters. One triaxial microelectromechanical accelerometer is installed on the left and right sides of the vibrating wheel bearing seat. The synchronous acquisition module of the edge computing gateway synchronously acquires data from all channels using a sampling frequency of 4000 Hz. The zero-crossing rising edge of the vertical component of the left triaxial microelectromechanical accelerometer in each excitation cycle is used as the synchronous triggering reference for all channels.

3. The method according to claim 1, characterized in that, In step two, a fourth-order Butterworth filter is used for filtering, with a passband range of 5 to 200 Hz. The specific process of the phase shift transform method is as follows: a fast Fourier transform is performed on each waveform record to obtain a complex spectrum. Each analysis frequency is selected sequentially within the frequency range of 5 to 100 Hz, and each test phase velocity is selected sequentially within the phase velocity range of 50 to 350 m / s. For each combination of analysis frequency and test phase velocity, the theoretical propagation time is calculated based on the horizontal distance between the engineering seismic detector and the vibrating wheel for each vertical component and the current test phase velocity, and converted into a phase offset. The complex value at the current analysis frequency is extracted from the complex spectrum of each waveform record and multiplied by the complex exponential factor corresponding to the phase offset to complete the phase correction. All 16 phase-corrected complex values ​​are summed and the modulus is taken as the dispersion energy value. The test phase velocity corresponding to the maximum dispersion energy value at each analysis frequency is selected as the Rayleigh wave phase velocity at the current analysis frequency.

4. The method according to claim 1, characterized in that, In step three, the initial velocity model is established as follows: The rectangular roadbed area directly below the vibrating wheel, with a horizontal width of 2.4 meters and a depth of 1.2 meters, is divided into an initial inversion grid. The initial inversion grid is uniformly divided into 8 columns of elements in the horizontal direction and 6 rows of elements in the depth direction, thus forming 48 rectangular grid elements. An initial shear wave velocity value is set for each rectangular grid element according to the Rayleigh wave dispersion curve. The setting method is as follows: Calculate the depth value of the center point of each rectangular grid element, find the Rayleigh wave phase velocity value that matches the wavelength of the depth value of the current rectangular grid element in the Rayleigh wave dispersion curve, and multiply the current Rayleigh wave phase velocity value by a coefficient of 0.92 to obtain the initial shear wave velocity value of the current rectangular grid element. The correspondence between the depth value and the wavelength is converted according to the principle that the wavelength is equal to 2.5 times the depth value.

5. The method according to claim 4, characterized in that, In step three, the specific process of identifying highly sensitive anomaly cells and adaptively refining the mesh is as follows: For each non-boundary rectangular mesh cell in the current inversion mesh, a neighborhood range centered on the current non-boundary rectangular mesh cell is determined. The neighborhood range includes four adjacent rectangular mesh cells: the left and right rectangular mesh cells adjacent to the current non-boundary rectangular mesh cell in the lateral direction, and the upper and lower rectangular mesh cells adjacent in the depth direction. The shear wave velocity sensitive kernel value of the current non-boundary rectangular mesh cell and the shear wave velocity sensitive kernel values ​​of the four adjacent rectangular mesh cells are extracted, and their arithmetic mean is calculated as the neighborhood sensitive kernel mean. The absolute value of the difference between the shear wave velocity sensitive kernel value of the current non-boundary rectangular grid cell and the mean value of the neighboring sensitive kernels is calculated as the local anomaly index of the current non-boundary rectangular grid cell. The maximum value of the local anomaly index of all non-boundary rectangular grid cells is selected as the anomaly benchmark value. The anomaly benchmark value is multiplied by a coefficient of 0.55 to obtain the anomaly judgment threshold. Non-boundary rectangular grid cells with a local anomaly index greater than the anomaly judgment threshold are marked as highly sensitive anomaly cells. Quadtree mesh refinement is performed on each highly sensitive anomaly cell to form 4 sub-cells. The shear wave velocity values ​​of the 4 sub-cells all inherit the shear wave velocity values ​​of the corresponding highly sensitive anomaly cells.

6. The method according to claim 5, characterized in that, In step three, the specific process of velocity model iterative update is as follows: traverse each rectangular grid cell in the current inversion grid, divide the shear wave velocity sensitive kernel value of the current rectangular grid cell by the maximum absolute value of the shear wave velocity sensitive kernel values ​​of all rectangular grid cells to obtain the normalized sensitive kernel value, multiply the normalized sensitive kernel value by the current shear wave velocity value and then multiply by the fixed iteration step size coefficient 0.03 to obtain the velocity update increment, and subtract the velocity update increment from the current shear wave velocity value to obtain the updated shear wave velocity value; Physical feasibility constraints are imposed on the updated shear wave velocity values. If the updated shear wave velocity value is less than 70 m / s, it is forcibly corrected to 70 m / s. If the updated shear wave velocity value is greater than 380 m / s, it is forcibly corrected to 380 m / s.

Images