Levee quality elastic wave tomography detection method

Abstract

The application discloses a kind of embankment quality elastic wave tomography detection methods, comprising: synchronous acquisition compactor vibration source at source end reference signal and the ground surface response signal of measured region surface, establish the unified time base of both;Correlation function of source end reference signal and ground surface response signal is calculated, and the coarse estimation time is determined based on peak value.Extract feature frequency component from signal and obtain analytic signal phase information, calculate the phase difference of ground surface response signal at coarse estimation time relative to source end reference signal, and use the phase difference to determine fine time correction amount.Use fine time correction amount to correct coarse estimation time, obtain final elastic wave time and generate elastic wave velocity distribution image by inversion.The application breaks through the limitation of sampling rate on measurement accuracy by collecting real source end waveform and combining phase refinement technology, and improves the resolution and reliability of compaction quality detection.

Description

Technical Field

[0001] This invention relates to the field of non-destructive testing and digital construction technology in geotechnical engineering, and in particular to a method for elastic wave tomography detection of embankment quality. Background Technology

[0002] The compaction quality of embankment subgrades determines the stability and service life of transportation infrastructure, making high-precision, full-coverage testing of compaction degree of significant engineering importance. Traditional testing methods, such as sand cone methods or nuclear density gauges, have limitations including high destructiveness, sparse measurement points, low efficiency, and radiation risks. Utilizing the continuous vibration generated by road rollers during construction as an elastic wave excitation source, combined with distributed intelligent sensor networks for tomographic imaging, allows for continuous monitoring of the underground structure in the compacted area without interrupting construction, providing crucial data support for digital construction quality control.

[0003] Current collaborative monitoring technologies primarily employ a passive source detection mode, which involves acquiring elastic wave signals through detector arrays deployed on the road surface and extracting propagation travel times using cross-correlation techniques. Existing solutions typically assume the roller's vibration signal is an ideal single-frequency sine wave, or simply generate a theoretical reference signal based on a set rated frequency. The processing system performs cross-correlation calculations between the acquired surface vibration signal and this theoretical reference signal, determines the arrival time of the elastic wave by searching for the peak time of the cross-correlation function, and then uses travel-time tomography algorithms to invert the underground wave velocity distribution.

[0004] However, when faced with complex construction environments, the accuracy of current measurement and the resolution of imaging are constrained by both signal source distortion and discrete sampling limitations. Specifically, the vibration of an actual road roller is not an ideal sine wave, but a complex signal containing mechanical transmission harmonics, frequency drift, and amplitude modulation. If only an ideal waveform is used as a reference, it will lead to broadening of cross-correlation peaks and severe phase mismatch. At the same time, the measurement resolution of traditional peak-picking-based methods is limited by the sampling frequency. For example, a 1-millisecond sampling interval corresponds to a distance error of about 150 mm, making it difficult to capture microsecond-level time differences. Furthermore, a single-frequency model cannot effectively eliminate multipath effects and dispersion interference, resulting in the final inverted compaction image failing to accurately reflect minute structural defects. Summary of the Invention

[0005] The purpose of this invention is to provide a method for detecting the quality of dikes using elastic wave tomography, in order to solve the aforementioned problems existing in the prior art.

[0006] Technical solution: A method for elastic wave tomography detection of levee quality, comprising:

[0007] Simultaneously acquire the source end reference signal at the vibration source of the road roller and the surface response signal of the measured area, and establish a unified time reference for the source end reference signal and the surface response signal;

[0008] Calculate the cross-correlation function between the source-end reference signal and the surface response signal, and determine the rough estimate to time based on the peak value of the cross-correlation function;

[0009] At least one characteristic frequency component is extracted from the source reference signal and the surface response signal, and the phase information of the analytical signal corresponding to each characteristic frequency component is obtained.

[0010] Based on the phase information of the analytical signal, the phase difference of the surface response signal relative to the source reference signal at the coarsely estimated time is calculated, and the fine time correction is determined by using the phase difference.

[0011] The coarse estimated arrival time is corrected by using a fine time correction to obtain the final elastic wave arrival time, and an elastic wave velocity distribution image of the measured area is generated by inversion based on multiple final elastic wave arrival times.

[0012] Beneficial effects: By acquiring real source waveforms and combining them with phase refinement technology, this invention overcomes the limitation of sampling rate on measurement accuracy, and improves the resolution and reliability of compaction quality detection. Attached Figure Description

[0013] Figure 1 This is a schematic diagram of the overall process of an elastic wave tomography detection method for dike quality provided in an embodiment of this application.

[0014] Figure 2 This is a schematic diagram of the process for establishing a unified time reference between the source reference signal and the surface response signal, provided in an embodiment of this application.

[0015] Figure 3 This is a schematic diagram of the process for determining the peak value and making a rough estimate based on the cross-correlation function, provided in an embodiment of this application.

[0016] Figure 4 This is a schematic diagram of the process for obtaining the analytical signal phase information corresponding to each characteristic frequency component provided in the embodiments of this application. Detailed Implementation

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

[0018] Example 1 details the complete execution flow of an elastic wave tomography detection method for levee quality, such as... Figure 1 As shown, by introducing a real reference waveform at the source end and combining it with a two-stage processing strategy of coarse estimation and phase refinement, the technical problem of traditional methods being limited by the sampling rate and unable to achieve high-precision imaging is solved. Taking the fundamental frequency vibration signal of a road roller as an example, the whole process from data acquisition to the final generation of velocity distribution image is shown.

[0019] Step 101: Simultaneously acquire the source end reference signal at the vibration source of the road roller and the surface response signal of the measured area, and establish a unified time reference for the source end reference signal and the surface response signal.

[0020] In this embodiment, the synchronous acquisition process is initiated first. An accelerometer at the vibration source of the road roller acquires the vibration signal of the drum in real time, i.e., the source-end reference signal r(t). Simultaneously, several inspection robots (detectors) distributed across the surface of the measured area synchronously acquire the vibration waveform transmitted from the ground, i.e., the ground response signal s(t). To ensure consistency of the time axes of both, the system employs a precise time synchronization mechanism based on the IEEE 1588 protocol, calibrating the local clock errors of all devices to the microsecond level. The received source-end reference signal is reconstructed based on data frames with absolute timestamps and mapped to the same local time coordinate system as the ground response signal, thus establishing a unified time reference.

[0021] Step 102: Calculate the cross-correlation function between the source reference signal and the surface response signal, and determine the coarse estimate time based on the peak value of the cross-correlation function.

[0022] After signal alignment, the processing unit first performs coarse-grained time estimation. DC removal and bandpass filtering preprocessing are applied to both the source-end reference signal r(t) and the surface response signal s(t) to suppress environmental noise. The normalized cross-correlation function Rt of both signals is then calculated using a Fast Fourier Transform algorithm. rs (τ). Within a preset physical search window, for example, a range determined by the minimum and maximum wave velocities, the system scans the cross-correlation function to find the maximum peak point. The time delay corresponding to this peak point is determined as a coarse estimate of the time τ. coarse Since the calculation process is based on discrete sampling points, the accuracy of this coarse estimate is limited by the sampling interval, such as 1 millisecond.

[0023] Step 103: Extract at least one characteristic frequency component from the source reference signal and the surface response signal respectively, and obtain the analytical signal phase information corresponding to each characteristic frequency component.

[0024] To overcome the accuracy limitations of the sampling rate, this step proceeds to the phase analysis stage. The system extracts the characteristic frequency component with the most concentrated energy from the broadband vibration signal. In the basic mode of this embodiment, this characteristic frequency component is selected as the fundamental frequency f0 of the road roller vibration, for example, 30 Hz. A narrowband filter with a center frequency of f0 is used to filter r(t) and s(t) to separate the fundamental frequency component. Hilbert transform is applied to the filtered signal to construct a complex analytic signal z(t). By calculating the arctangent values ​​of the imaginary and real parts of the analytic signal, the instantaneous phase φ of the source-end reference signal is obtained. r (t) and the instantaneous phase φ of the surface response signal s (t).

[0025] Step 104: Based on the phase information of the analytical signal, calculate the phase difference of the surface response signal relative to the source reference signal at the coarsely estimated time, and use the phase difference to determine the fine time correction amount.

[0026] This step utilizes phase information to calculate the sub-period-level time difference. The system reads the surface response signal and roughly estimates the time τ. coarse The observed phase φ at time s (τ coarse ), and the reference phase φ of the source-end reference signal at the start time. r (0). Calculate the difference between the two, Δφ = φ s (τ coarse )-φ r (0), and normalize this difference to the interval from -π to +π to obtain the phase difference. This phase difference represents a small deviation of less than one period between the actual propagation time and the coarsely estimated time. Using the formula Δτ... fine =Δφ / (2π*f0) converts the phase difference into a fine time correction Δτ in the time domain. fine , where Δτ fine The unit is seconds, and f0 is the fundamental frequency of vibration.

[0027] Step 105: Correct the coarsely estimated arrival time using fine time correction to obtain the final elastic wave arrival time, and generate an elastic wave velocity distribution image of the measured area based on multiple final elastic wave arrival times.

[0028] The system synthesizes the final measurement results. The coarse estimate obtained in step 102 is then applied to time τ. coarse The integer periodic portion and the fine time correction Δτ obtained in step 104 fine Adding them together yields a highly accurate final elastic wave time τ. finalThe system collects the final elastic wave arrival time data measured by all inspection robots within the measured area and constructs an observation vector. Combining the position coordinates of each robot, a ray path matrix is ​​constructed, and inversion calculations are performed using the Joint Iterative Reconstruction Technique (SIRT). By iteratively updating the slowness values ​​of the grid cells until the residuals between the theoretical and observed arrival times converge, an elastic wave velocity distribution image reflecting the compactness of the subsurface medium is finally generated to guide road compaction construction.

[0029] Example 2 describes how to obtain a high-fidelity source-end reference signal from the harsh working conditions of a road roller and transmit it to the receiving end in real time through a wireless link with limited bandwidth, that is, describe the hardware basis for realizing source-end phase enhancement.

[0030] Step 201: The source-end reference signal is a vibration acceleration signal. The source-end reference signal at the vibration source of the road roller is obtained by: acquiring the vibration acceleration signal in real time through an acceleration sensor installed at the bearing seat of the vibrating drum of the road roller, wherein the sampling frequency of the acceleration sensor is more than 10 times the fundamental frequency of the vibration source.

[0031] In this embodiment, to accurately capture the true waveform of the vibration source, the installation position and specifications of the sensor must be appropriately selected. The accelerometer is installed on the bearing housing of the vibrating drum of the road roller, and the distance between the installation position and the rotation axis of the vibrating drum should preferably be controlled within 200 mm. The reason for this position selection is that the bearing housing is the first rigid node for transmitting the vibration force to the frame, which can most directly reflect the actual force of the vibrating drum on the ground and avoid the attenuation of high-frequency signals by the frame damping system.

[0032] The preferred accelerometer is a triaxial MEMS accelerometer with a range of ±50g and a sensitivity better than 1mg. Regarding the selection of the sampling frequency, considering that the typical fundamental vibration frequency f0 of a road roller is between 25Hz and 50Hz, and that the mechanical transmission contains abundant second and third harmonic components, to prevent frequency domain aliasing and retain sufficient phase information, the sensor's sampling frequency f0 is selected. ref The frequency is set to 500 Hz or higher, which meets the requirement of being more than 10 times the fundamental frequency. The sensor continuously collects the vibration acceleration signal a of the vibrating drum. ref (t), forming the original data stream of the source-side reference signal.

[0033] The selection of the sampling frequency needs to consider the following factors: If the sampling frequency is too low, such as below 10 times the fundamental vibration frequency, higher harmonic components may be lost due to frequency domain aliasing, affecting the accuracy of subsequent multi-frequency phase analysis; if the sampling frequency is too high, such as exceeding 2000 Hz, it will significantly increase the burden on data storage and wireless transmission. In typical engineering applications, a sampling frequency between 500 Hz and 1000 Hz is recommended to achieve a balance between signal fidelity and system resource consumption.

[0034] Step 202: Establish a unified time reference for the source-end reference signal and the surface response signal, such as... Figure 2 As shown, the specific steps include: performing anti-aliasing filtering on the acquired source-end reference signal, and compressing the processed signal using differential pulse code modulation to obtain a compressed reference waveform; performing spectral analysis on the source-end reference signal to determine the current vibration fundamental frequency; broadcasting synchronization data frames by the road roller end at a preset period, the synchronization data frames containing the compressed reference waveform, the timestamp of the current reference signal, and the current vibration fundamental frequency parameters; receiving the synchronization data frames by the receiving end, parsing out the compressed reference waveform and decoding it to reconstruct the source-end reference signal, and mapping the time axis of the source-end reference signal to the local time system based on the timestamp.

[0035] Due to the limited bandwidth of wireless communication, directly transmitting high-sampling-rate raw waveform data would result in significant latency, impacting real-time performance. This step employs a signal preprocessing and compressed transmission strategy. The vibration controller performs anti-aliasing low-pass filtering on the raw reference signal. Specifically, a fourth-order Butterworth low-pass filter is used, with a cutoff frequency set to 200 Hz to filter out high-frequency mechanical noise. The filtered signal is then downsampled, for example, reducing the sampling rate to 250 Hz.

[0036] In the data compression stage, to minimize the data volume while preserving waveform characteristics, this embodiment employs Differential Pulse Code Modulation (DPCM) technology. First, the signal is normalized to obtain the normalized signal a. norm (t). Then calculate the difference between adjacent sampling points.

[0037] d(n)=a norm (n)-a norm (n-1);

[0038] In the formula, d(n) represents the difference value of the nth sampling point; a norm (n) represents the normalized amplitude of the nth sampling point; a norm(n-1) represents the normalized amplitude of the (n-1)th sampling point. The difference value d(n) is quantized and encoded using 4 bits. Compared to the original 16-bit or 24-bit sampled data, this method achieves a compression ratio of over 4:1, resulting in a compressed reference waveform.

[0039] To distribute this reference waveform to all inspection robots, the vibration controller broadcasts a synchronization data frame at fixed intervals, such as 1 second. This synchronization data frame uses a compact binary structure and specifically includes the following fields:

[0040] The frame synchronization header is fixed at 0xAA55AA55 and is used by the receiver to identify the start of the frame.

[0041] Frame sequence number, used to detect packet loss;

[0042] A timestamp records the absolute start time of the reference signal for that frame, accurate to microseconds;

[0043] The fundamental frequency f0 of the vibration has a resolution of 0.01 Hz; the vibration amplitude A0 has a resolution of 0.1 mm.

[0044] The data body, namely the aforementioned DPCM compressed reference waveform data, for example, 1 second of data is compressed to approximately 125 bytes;

[0045] Cyclic Redundancy Check (CRC-16) is used to ensure data integrity.

[0046] The broadcast data uses LoRa spread spectrum modulation technology, with a center frequency set at 470 MHz, a bandwidth of 125 kHz, and a spreading factor (SF) of 7. This combination of parameters can provide an effective data rate of approximately 5.5 kbps while ensuring coverage, which is sufficient to support the aforementioned compressed data stream.

[0047] After receiving the synchronization data frame, each inspection robot first performs a cyclic redundancy check (CRC) check. If the check passes, it proceeds with decompression. Decompression is the reverse of compression; by accumulating the decoded differential values ​​and recovering the amplitude, the source reference signal r(t) is reconstructed. The robot extracts the absolute timestamp within the frame and, combined with the clock offset between its local machine and the master node, precisely maps the time axis of the reference signal to its local time system. Thus, each robot obtains a source reference waveform copy that is strictly time-aligned with the surface response signal it acquired, providing an accurate template for subsequent cross-correlation analysis.

[0048] In some optional implementations, when a slave node detects the loss of several consecutive synchronization data frames, such as five, it enters a communication hold mode. In hold mode, the slave node can choose to pause the acquisition of new data or approximate the data using a reference waveform from the most recently successfully received synchronization data frame, while simultaneously sending a communication recovery request to the master node. Once the communication link is restored, the slave node re-executes the clock synchronization process and marks the data acquired during hold mode as having synchronization quality to be verified, giving it lower weight or excluding it in subsequent inversion calculations.

[0049] In some alternative implementations, if high-bandwidth communication conditions are available on-site, such as the deployment of Wi-Fi 6 or a 5G private network, the compression step described above can be bypassed, and the vibration controller can directly broadcast the uncompressed raw floating-point waveform data, avoiding the slight impact of quantization noise on the accuracy of subsequent phase analysis.

[0050] Example 3 describes how to overcome the sampling rate limitation and achieve sub-period-level high-precision time measurement using single-frequency phase information under normal operating conditions. The method in this example constitutes the basic algorithm layer of the entire high-precision imaging system. Through a two-stage strategy of coarse estimation and fine refinement, it solves the technical problem that the resolution of traditional cross-correlation methods is limited by the sampling interval.

[0051] Step 301: Calculate the cross-correlation function between the source-end reference signal and the surface response signal. Based on the peak value of the cross-correlation function, determine the coarse estimate of the time, such as... Figure 3 As shown, the process includes: removing the DC components of the source reference signal and the surface response signal respectively, and filtering out out-of-band noise using a bandpass filter; calculating the normalized cross-correlation function of the preprocessed source reference signal and the surface response signal; searching for the maximum value of the normalized cross-correlation function within a preset search time window, and determining the time delay corresponding to the maximum value as the coarse estimate time.

[0052] To accurately extract relevant features from the signal, preprocessing of the original signal is necessary. The source-side reference signal r(t) and the surface response signal s(t) often contain a DC bias introduced by hardware circuitry, which can affect the accuracy of related calculations. The system first calculates the mean of the signal within the time window and subtracts it from the original signal to eliminate the DC component. A bandpass filter is then used for filtering. The center frequency of the bandpass filter is set to the fundamental vibration frequency f0, and the bandwidth is preferably set to twice f0. For example, when the fundamental vibration frequency is 30 Hz, the passband range of the filter can be set from 0 Hz to 60 Hz. This step effectively filters out high-frequency mechanical noise and low-frequency drift interference outside the frequency band, significantly improving the signal-to-noise ratio.

[0053] After preprocessing, the system calculates the normalized cross-correlation function of the two signals. To improve computational efficiency, in practical engineering implementations, the Fast Fourier Transform (FFT) is typically used to perform the convolution operation in the frequency domain. The normalized cross-correlation function R... rs (τ) can characterize the similarity between two signals under different time delays τ, and its value ranges from -1 to +1.

[0054] R rs (τ)=[∫r'(t)*s'(t+τ)dt] / [sqrt(∫r' 2 (t)dt)*sqrt(∫s' 2 (t)dt)];

[0055] In the formula, R rs (τ) represents the normalized cross-correlation function; r'(t) and s'(t) represent the preprocessed source reference signal and the surface response signal, respectively; τ is the time delay variable; sqrt(...) represents the square root operation; ∫...dt represents the integral within the effective signal window.

[0056] The system searches for R within a preset search time window. rs The peak value of (τ). The range of the search time window [τ] min ,τ max Determined based on the priori range of source-detection distance and medium wave velocity. Specifically, τ min τ equals the minimum source-receiver distance divided by the maximum wave velocity estimate. max It equals the maximum source-receiver distance divided by the minimum wave velocity estimate. The location of the maximum value within this window is then determined as the coarse estimate up to time τ. coarse Simultaneously, the system will record the peak correlation coefficient ρ at that location. max If ρ max If the quality is below a preset threshold, such as 0.3, the data set will be marked as invalid and will not be included in subsequent fine-grained calculations.

[0057] In some alternative implementations, to further improve the robustness of the coarse estimate, a generalized cross-correlation (GCC) algorithm can be used, for example, by using phase transformation (PHAT) weighting to sharpen the correlation peaks and more accurately locate the peak positions in a strongly reverberant environment.

[0058] Step 302: Obtain the analytical signal phase information corresponding to each characteristic frequency component, such as... Figure 4 As shown, the process includes: using a narrowband filter whose center frequency corresponds to the characteristic frequency component to filter the source reference signal and the surface response signal respectively; performing a Hilbert transform on the filtered signal to construct an analytic signal in complex form; and calculating the instantaneous phase that changes with time based on the real and imaginary parts of the analytic signal, as the phase information of the analytic signal.

[0059] After obtaining a coarse time estimate, phase information of the signal needs to be extracted to further improve time accuracy. In this embodiment, the characteristic frequency component refers to the fundamental vibration frequency f0 of the road roller. The system uses a narrowband bandpass filter to perform secondary filtering on the signal, and the bandwidth B of this filter is... n The range is set relatively narrow, preferably 0.1 times f0, to filter out sidelobe interference to the maximum extent and extract the pure fundamental frequency component.

[0060] An analytic signal is constructed using the Hilbert transform. An analytic signal is a complex signal whose real part is the original real signal and whose imaginary part is the Hilbert transform of the original real signal. Through the analytic signal, the amplitude envelope and instantaneous phase of the vibration signal can be separated.

[0061] z(t) = x(t) + j*H[x(t)];

[0062] In the formula, z(t) represents the constructed analytic signal; x(t) represents the real signal after narrowband filtering; j represents the imaginary unit; H[...] represents the Hilbert transform operator.

[0063] Based on the real and imaginary parts of the analytic signal, the instantaneous phase that changes continuously with time can be calculated using the four-quadrant arctangent function atan2.

[0064] φ(t)=atan2(Im[z(t)],Re[z(t)]);

[0065] In the formula, φ(t) represents the instantaneous phase at time t, with a value ranging from negative π to positive π; Im[...] indicates taking the imaginary part; Re[...] indicates taking the real part. This step generates the phase function φ for the source reference signal and the surface response signal, respectively. r (t) and φ s (t).

[0066] Step 303: The characteristic frequency components only include the fundamental vibration frequency; based on the analytical signal phase information, calculate the phase difference of the surface response signal relative to the source reference signal at the coarsely estimated time, and use the phase difference to determine the fine time correction amount, specifically including:

[0067] The observed phase of the surface response signal at the coarsely estimated time is obtained, as well as the reference phase of the source reference signal at the time corresponding to the timestamp of the acquisition start time or the synchronization data frame; the difference between the observed phase and the reference phase is calculated, and the difference is normalized to one vibration period to obtain the phase difference; the fine time correction is calculated using the formula Δτ=Δφ / 2πf0, where Δτ is the fine time correction, Δφ is the phase difference, and f0 is the vibration fundamental frequency.

[0068] The system first reads the instantaneous phase of the source-end reference signal at the initial time t=0, and uses it as the phase reference φ. ref =φ r (0). At a rough estimate, τ... coarse At the corresponding moment, the instantaneous phase φ of the surface response signal is read. obs =φ s (τ coarse ).

[0069] Furthermore, the difference between the two is calculated and then unwound or normalized to fall within the range of negative π to positive π, thus obtaining the phase difference Δφ.

[0070] Δφ=φ obs -φ ref ;

[0071] If Δφ is greater than π, then subtract 2π; if Δφ is less than or equal to -π, then add 2π.

[0072] This phase difference represents the small delay less than one period in the propagation time of the elastic wave. Converting this phase difference into a time-domain correction:

[0073] Δτ fine =Δφ / (2π*f0);

[0074] In the formula, Δτ fine Δφ represents the fine time correction in seconds; Δφ represents the normalized phase difference in radians; and f0 represents the fundamental frequency of vibration in Hertz.

[0075] The final elastic wave arrives at τ final It consists of a full-cycle part and a fine-tuning part. The number of full cycles N is determined by the coarse estimate, i.e., N = round(τ). coarse *f0), where round(·) is the rounding function.

[0076] τ final =(N+Δφ / (2π)) / f0.

[0077] As an example, assuming the fundamental frequency f0 of the road roller vibration is 30 Hz, corresponding to a period of approximately 33.33 milliseconds, the time τ obtained by coarse estimation using cross-correlation is... coarse The time interval is 100 milliseconds. The phase difference Δφ measured by the system is π / 5, or 36 degrees. First, calculate the number of integer cycles N = round(100 * 30 / 1000) = 3, then calculate the fine time correction:

[0078] Δτ fine =(π / 5) / (2π*30)=1 / 300≈3.33 milliseconds.

[0079] The final high-precision time τfinal =3*33.33+3.33=103.33 milliseconds.

[0080] As can be seen, compared with the coarse estimation results, this method corrects the time deviation, and the accuracy is no longer limited to the sampling interval of 1 millisecond or 2 milliseconds, and can theoretically reach the microsecond level.

[0081] By employing the phase refinement method described above, this invention effectively solves the problem of measurement accuracy limitations caused by sampling frequency. Traditional cross-correlation peak picking methods are limited in temporal resolution by the sampling interval; for example, a 1-millisecond sampling interval corresponds to a distance error of approximately 150 millimeters. This invention, by introducing phase information, elevates the accuracy from the sampling interval level to the phase resolution level. This improvement enables the system to detect finer compaction differences and smaller-scale structural defects.

[0082] The analysis of the above embodiments shows that the method of the present invention has a significant improvement in measurement accuracy compared with the traditional pure cross-correlation peak picking method. The accuracy is improved from being limited by the sampling interval (millisecond level) to being determined by the phase resolution (microsecond level). Those skilled in the art will understand that the specific improvement may vary depending on the signal-to-noise ratio, propagation distance, and medium conditions at the scene.

[0083] Example 4: In view of the inherent physical characteristics of the vibration signal of the road roller containing rich harmonics, this example adopts the multi-frequency joint constraint method to further solve the periodic ambiguity problem that may exist in the single-frequency method, and uses the information redundancy of multiple frequency components to improve the noise resistance of the measurement.

[0084] The characteristic frequency components include the fundamental frequency of vibration and at least one higher harmonic component whose frequency is an integer multiple of the fundamental frequency of vibration; based on the phase information of the analytical signal, the phase difference of the surface response signal relative to the source-end reference signal at the coarsely estimated time is calculated, including:

[0085] Step 401: For each characteristic frequency component, calculate the difference between the observed phase of the surface response signal at the coarsely estimated time and the reference phase of the source reference signal at the initial time, to obtain a set of multi-frequency phase differences corresponding to different frequencies.

[0086] In this embodiment, the characteristic frequency components are no longer limited to the fundamental frequency, but include the fundamental frequency and its integer multiples thereof, as well as higher harmonics. Based on the spectral analysis results of the source-end reference signal, the system automatically determines a set of valid harmonic frequency sequences:

[0087] f k =k*f0, k=1,2,...,K;

[0088] Among them, f kf0 is the k-th harmonic frequency, in Hz; k is the harmonic order; K is the highest harmonic order being analyzed; f0 is the fundamental frequency. The value of K is determined based on the spectral analysis of the source-end reference signal, typically K=3, which means analyzing the fundamental frequency, second harmonic, and third harmonic components.

[0089] The system constructs a multi-channel filter bank containing multiple narrowband bandpass filters with center frequencies of f1, f2, and f3. The source reference signal and the ground response signal are separated into component signals of different frequency bands by this filter bank. For each frequency band k, the system extracts the analytic signal phase using the Hilbert transform, following the method described in the previous embodiment.

[0090] The system calculates the time τ for each frequency component k in a coarse estimate. coarse Phase difference Δφ at the point k Therefore, we obtain not just a single phase difference value, but a set of phase differences containing multiple elements:

[0091] {Δφ1,Δφ2,Δφ3}.

[0092] This set contains richer time constraint information than single-frequency measurements.

[0093] Step 402, using phase differences to determine the fine time correction amount, specifically includes: constructing a multi-frequency consistency cost function, which is used to characterize the degree of consistency between the estimated phase corresponding to the candidate time correction amount and each observation value in the multi-frequency phase difference set; generating a candidate time correction amount set within a preset range centered on the coarsely estimated time; searching the candidate time correction amount set; and determining the candidate time correction amount that makes the multi-frequency consistency cost function reach its minimum value as the fine time correction amount.

[0094] To comprehensively utilize the aforementioned multi-frequency phase information, this embodiment introduces the concept of a multi-frequency consistency cost function. The actual propagation time should ensure that the phases calculated from all frequency components match the observed phases. If the candidate time correction is incorrect, then even if the fundamental frequency phase happens to match, the phases of the high-frequency harmonics will usually show significant deviations.

[0095] The multi-frequency consistency cost function J(τ) is specifically defined as follows:

[0096] J(τ)=Σ[w k *sin 2 (π*f k *(τ-τ φ_k ))];

[0097] In the formula, J(τ) represents the inconsistency cost of the candidate time τ; Σ represents the summation over all effective harmonic orders k; w kIt is the weighting factor for the k-th frequency component; f k It is the kth harmonic frequency; τ φ_k The time delay value, τ, is calculated separately based on the phase difference of the k-th frequency. φ_k =Δφ k / (2π*f k sin 2 The (...) term is used to measure the candidate time τ and the measured value τ. φ_k The phase distance between them is 0 when they match, i.e., they differ by an integer number of periods.

[0098] Step 403: Based on the spectral analysis results of the source-end reference signal, determine the signal-to-noise ratio of each characteristic frequency component, and assign weighting factors to each phase deviation term in the multi-frequency consistency cost function based on the signal-to-noise ratio; construct a discrete set of candidate time corrections within a preset range with the vibration period corresponding to the fundamental vibration frequency as the step size; calculate the multi-frequency consistency cost function value corresponding to each element in the candidate time correction set in order to find the minimum value of the weighted consistency.

[0099] The allocation of weights is crucial when constructing the cost function. Since different harmonic components have varying energy intensities and are affected by noise to varying degrees, a one-size-fits-all approach cannot be applied. This step first calculates the signal-to-noise ratio (SNR) of each harmonic component in the source-end reference signal. k .

[0100] SNR k =10*log10(P signal_k / P noise_k );

[0101] In the formula, P signal_k P represents the power spectral density at the k-th harmonic frequency. noise_k This represents the noise floor power in the adjacent frequency band.

[0102] Weighting factor w k Normalized allocation based on signal-to-noise ratio:

[0103] w k =SNR k / (ΣSNR j' );

[0104] In this context, the summation index j' indicates that all valid harmonic indices have been traversed.

[0105] In other words, frequency components with higher signal-to-noise ratios have greater influence in the final decision-making process, thus ensuring the robustness of the algorithm.

[0106] Next, we will search for the optimal solution. The search range is set within the range of time τ (roughly estimated). coarse Within the neighborhood, for example [τ coarse-T0 / 2,τ coarse +T0 / 2], where T0 is the fundamental frequency period. Within this range, a series of discrete candidate time points are generated using the vibration period corresponding to the fundamental frequency or a finer grid as the step size. The system calculates the cost function value J(τ) corresponding to each of these candidate points and finds the point that minimizes J(τ) as the optimal candidate time point τ. opt This process utilizes the short-period characteristics of high-frequency harmonics to eliminate the ambiguity of low-frequency measurements, while simultaneously using the long-period characteristics of low frequencies to prevent cycle slip in high-frequency measurements, achieving a measurement effect similar to that of a vernier caliper.

[0107] Step 404: Calculate the normalized weights of each characteristic frequency component based on the signal-to-noise ratio; calculate the difference between the observed phase of each characteristic frequency component and the theoretical phase calculated based on the current candidate time correction when the multi-frequency consistency cost function reaches its minimum value, and use this difference as the residual phase deviation corresponding to each characteristic frequency component; use the normalized weights to perform a weighted summation of the residual phase deviations of each characteristic frequency component, and superimpose the weighted results to the candidate time correction when the minimum value is reached, to obtain the final fine time correction.

[0108] When the optimal candidate is found, τ opt Subsequently, this value is typically located on a discrete search grid and may not yet be a precise extreme point. To further explore the accuracy potential, this step performs local phase refinement.

[0109] The system first calculates at τ opt At each frequency component k, the residual phase deviation still exists, corresponding to the time delay δτ. k .

[0110] δτ k =[φ obs_k -(φ ref_k +2π*f k *τ opt )] / (2π*f k );

[0111] The part in parentheses represents the phase residual after deducting the integer period.

[0112] Specifically, for the k-th frequency component, its residual phase deviation is defined as: at the candidate time τ opt The measured phase φ of the surface response signal at that location obs_k Compared with τ opt The calculated theoretical phase, i.e., the reference phase φ of the source-side reference signal. ref_k Add 2π·f k ·τ opt The difference between them.

[0113] Using the aforementioned weight wk Take a weighted average of these residual delays:

[0114] Δτ residual =Σ(w k *δτ k );

[0115] The weighted average residual correction is then added to the optimal candidate time to obtain the final fine-grained time correction τ. final =τ opt +Δτ residual .

[0116] Using the above method, the final measurement result integrates the contributions of all harmonic components, which is equivalent to utilizing the effective energy of all frequency bands. Its statistical error variance will be significantly smaller than that of single-frequency measurement, achieving a precision enhancement effect of 1+1>2.

[0117] Example 5 describes a quality control mechanism based on multi-frequency consistency. In complex construction sites, strong noise interference or special geological structures, such as cavities or weak interlayers, can cause conventional time-of-arrival measurements to fail. By analyzing the residual information in the multi-frequency inversion process, not only can unreliable measurement data be automatically eliminated, but the dispersion characteristics of the medium can also be further identified, providing a deeper basis for identifying potential engineering hazards.

[0118] Step 501: The multi-frequency consistency cost function value corresponding to the determined fine time correction amount is used as the consistency residual; the consistency residual is compared with the preset quality threshold; if the consistency residual exceeds the quality threshold, the corresponding final elastic wave arrival time is marked as invalid data, or it is identified as abnormal data indicating that the medium has a dispersion effect; if the consistency residual does not exceed the quality threshold, the corresponding final elastic wave arrival time is marked as valid data and included in the subsequent inversion calculation.

[0119] After the multi-frequency joint inversion is completed, the system obtains an optimal fine-grained time correction that minimizes the multi-frequency consensus cost function J(τ). This step defines the minimum value as the consensus residual J. res The residual is not merely an optimization index; it also carries a clear physical meaning: it reflects the degree of inherent conflict among different frequency components when calculating the propagation time of the same physical path. In an ideal homogeneous medium, elastic waves of different frequencies should have the same group velocity, therefore J res It should approach zero.

[0120] The system presets a quality threshold J. th As a preferred embodiment, the threshold can be determined through statistical analysis of known high-quality samples, for example, by setting it to 0.1.

[0121] When the calculated Jres ≤J th When the measurement data is valid, the system determines that it is valid and includes it in the subsequent imaging database.

[0122] Conversely, if J res >J th This indicates a significant conflict in the phase information of each frequency component. In this case, the system initiates a tiered processing logic: checking the signal-to-noise ratio (SNR). If the conflict is due to random phase disturbances caused by low SNR, the data point is directly marked as invalid and discarded to prevent noise contamination of the imaging results. If the signal SNR is acceptable, but high consistency residuals still exist, it usually suggests that the propagation medium itself has complex physical characteristics, such as dispersion effects or nonlinear responses. In this case, the system marks the data point as an anomaly and retains it, transferring it to the anomaly analysis module for further diagnosis, rather than simply discarding it.

[0123] Step 502: For each characteristic frequency component, calculate the corresponding optimal candidate correction amount based on the phase difference of the frequency component; calculate the time difference between the optimal candidate correction amount corresponding to each higher harmonic component and the optimal candidate correction amount corresponding to the fundamental frequency; if the time difference shows a monotonically changing trend with the increase of the frequency order, and the absolute value of the time difference exceeds the preset dispersion threshold, then it is determined that there is a dispersion effect in the current propagation path.

[0124] This step details how to use the aforementioned anomalous data to diagnose the dispersion characteristics of a medium. Dispersion refers to the phenomenon where the propagation speed of a wave changes with frequency. In compaction testing, this often corresponds to specific geological structures, such as aquifers or loose interlayers.

[0125] The system extracts each characteristic frequency component, namely the fundamental frequency and each harmonic, and calculates the optimal candidate correction value separately, denoted as τ. opt_k Calculate the time difference Δτ between each higher harmonic and the fundamental frequency. disp_k .

[0126] Δτ disp_k =τ opt_k -τ opt_1 ;

[0127] In the formula, τ opt_k τ represents the optimal candidate correction for the k-th harmonic; opt_1 This represents the optimal candidate correction value corresponding to the fundamental frequency.

[0128] The system analyzes the trend of time difference as a function of frequency order k. If Δτ disp_k If the propagation path exhibits a monotonically increasing or monotonically decreasing trend as k increases, and its absolute value exceeds a preset dispersion threshold, such as 2 milliseconds, then the propagation path is determined to have a dispersion effect.

[0129] Specifically, in normally compacted soil, high-frequency components decay rapidly, but usually slightly faster than low-frequency components (normal dispersion). If an anomalous reverse dispersion phenomenon occurs, i.e., high-frequency components lag significantly, it may indicate the presence of a microcrack network or uncompacted loose clumps along the path. These structures exert a strong scattering delay effect on short-wavelength high-frequency waves. Through this logic, this invention transforms data that might otherwise be considered error into effective features for detecting hidden defects.

[0130] As a specific application example, suppose that in the detection of a certain propagation path, the system calculates:

[0131] The optimal candidate correction for the fundamental vibration frequency (30 Hz) is τ. opt_1 =52.3 milliseconds, the optimal candidate correction for the second harmonic (60 Hz) is τ. opt_2 =53.1 milliseconds, the optimal candidate correction for the third harmonic (90 Hz) is τ. opt_3 =53.8 milliseconds.

[0132] Based on this, the time difference Δτ between the second harmonic and the fundamental frequency is calculated. disp_2 =53.1-52.3=0.8 milliseconds, the arrival time difference Δτ between the third harmonic and the fundamental frequency. disp_3 =53.8 - 52.3 = 1.5 milliseconds. Observation shows that the time difference increases monotonically with increasing frequency order, and all values ​​exceed the preset dispersion threshold, for example, 0.5 milliseconds. Based on this, the system determines that there is a dispersion effect in this propagation path. Combining engineering experience, the inverse dispersion phenomenon where high-frequency components lag significantly usually indicates the presence of uncompacted loose clumps or microcrack networks along the path. The system automatically marks this area as a suspected defect area on the compaction distribution map and prompts on-site personnel to conduct focused verification or recompaction.

[0133] Example 6 describes how discrete time-delayed data can be transformed into an intuitive image of compaction distribution, ultimately serving construction quality control. It supplements the specific implementation formula of the SIRT algorithm of the joint iterative reconstruction technology, and constructs a complete mapping chain from physical parameters to engineering indicators.

[0134] Step 601: Based on the spatial coordinates corresponding to the source reference signal and the surface response signal, the measured area is divided into grids, and a ray path matrix is ​​constructed. The spatial coordinates are obtained through a positioning system or pre-calibrated. Using the Joint Iterative Reconstruction Technology (SIRT), the theoretical arrival time is calculated based on the initial velocity model. The residual between the final elastic wave arrival time and the theoretical arrival time is calculated, and the residual is back-projected to each grid cell using the ray path matrix to update the velocity model until the iteration converges, thus obtaining the elastic wave velocity distribution image.

[0135] To reconstruct the internal structure of the measured area, spatial discretization is first required. The system determines the boundary of the imaging area based on the coordinate range of the road roller and the inspection robot, and divides it into regular three-dimensional mesh units. For example, the mesh size can be set to 0.5 m × 0.5 m × 0.2 m. Let the total number of meshes be M, and the total number of rays be N.

[0136] Construct the ray path matrix L. Based on the high-frequency ray approximation theory, assuming the elastic wave propagates in a straight line (in the initial iteration), calculate the length l of the line segment through the j-th grid cell of the i-th ray. ij This forms an N-row, M-column sparse matrix.

[0137] The inversion process employs the Joint Iterative Reconstruction Technique (SIRT). First, an initial velocity model is established, typically a uniform model, for example, where the slowness (the reciprocal of the velocity) of all grid cells is S. j (0) All were initialized to 1 / 150 of a second per meter.

[0138] In the k-th iteration, the slow-speed model is updated according to the following core formula:

[0139] ΔS j (k) =λ*[Σ(Δt i (k) *l ij / L i )] / [Σ(l ij 2 / L i )];

[0140] S j (k+1) =S j (k) +ΔS j (k) ;

[0141] In the formula, ΔS j (k) Δt represents the slow update amount of the j-th grid cell in the k-th iteration; λ is the relaxation factor, which is preferably between 0.5 and 1.0 to ensure the stability of iteration convergence; the first Σ represents the summation of all rays i passing through grid cell j; Δt i (k) Let l be the residual between the observed time of the i-th ray and the theoretical arrival time calculated by the current model; ij L is the length of the i-th ray within the j-th grid. i Let be the total path length of the i-th ray.

[0142] The system repeatedly executes the above process of calculating the theoretical timeout, calculating the residuals, and backprojecting updates until the root mean square (RMS) of the residuals for all rays is less than a preset convergence threshold, such as 0.5 milliseconds, or until the maximum number of iterations is reached. The finally converged slow-speed model is then inversely taken and converted into an elastic wave velocity distribution image V=[v1,v2,...,v...]. M ].

[0143] As an alternative, if computational resources are limited, the algebraic reconstruction technique ART can be used instead of SIRT. ART updates the model using one ray at a time, resulting in faster convergence but sensitivity to noise; while the preferred SIRT in this embodiment updates the model using all rays simultaneously, offering better noise resistance and smoothness, and is more suitable for construction site data.

[0144] Step 602: Invoke the pre-built velocity-compaction mapping relationship to convert the elastic wave velocity distribution image into compaction distribution data; compare the compaction distribution data with the preset compaction qualification threshold to identify and mark the under-compaction area or defect location where the compaction degree is lower than the qualification threshold.

[0145] While physical velocity models reflect the underground structure, construction workers are more concerned with direct engineering indicators. This step requires establishing a mapping model between velocity v and compaction degree K. Before construction, technicians select typical areas on site, measure the compaction degree using the sand cone method or nuclear density meter, and simultaneously measure the wave velocity, fitting a linear regression equation:

[0146] K q =a*v q +b;

[0147] In the formula, K q v represents the compaction percentage of the q-th grid cell; q Let be the elastic wave velocity of the unit; a and b are calibration coefficients. For typical subgrade fill materials, a is usually a positive value, indicating that the higher the velocity, the better the compaction.

[0148] This equation is used to convert the velocity image point by point into compaction degree distribution data. A qualified threshold K is set according to engineering design specifications. pass For example, 96%, and the warning threshold K warn For example, 94%.

[0149] Traverse all grid cells to generate a visualized 3D chromatogram:

[0150] Compaction degree greater than or equal to K pass The area is rendered in green, indicating that it is qualified;

[0151] Between K warn With K passThe area between them is rendered in yellow, indicating that reinforcement is needed due to undervoltage;

[0152] Below K warn The areas are rendered in red, representing serious defects or holes.

[0153] In addition, the system will automatically extract the boundary coordinates of the red and yellow areas to generate a list of specific defect locations.

[0154] Furthermore, to achieve closed-loop control, the aforementioned defect location data is encapsulated into JSON (JavaScript Object Notation) format instruction packets and transmitted in real time to the roller's onboard navigation terminal via a wireless network. The roller automatically plans its compaction path based on the instructions, enabling an intelligent construction mode of simultaneous compaction and inspection with immediate feedback. This fully automated process eliminates the risk of human error in inspections and improves the quality control level of embankment construction.

[0155] According to one aspect of this application, the source-end reference signal acquisition and coordinated triggering can further be:

[0156] A triaxial MEMS accelerometer is installed in the bearing housing of the vibrating drum of the rolling mill. The sensor is positioned no more than 200 mm from the axis of rotation of the vibrating drum. The sensor has a range of ±50 g, a sensitivity better than 1 mg, and a frequency response range covering 0.5 Hz to 200 Hz.

[0157] The reference sensor uses a sampling rate f ref Continuously acquire vibration acceleration signal a of the vibrating drum ref (t), where f ref ≥500Hz.

[0158] The vibration controller preprocesses the original reference signal:

[0159] A 4th-order Butterworth low-pass filter is used, with a cutoff frequency f. c =200Hz, filtering out high-frequency noise and aliasing components.

[0160] The filtered signal is downsampled to 250Hz to reduce the data transmission bandwidth requirements.

[0161] Perform amplitude normalization on the signal:

[0162] a norm (t)=a ref (t) / max{|a ref (t)|};

[0163] In the formula, a norm (t) is the normalized reference signal; max{|a ref (t)|} represents the maximum absolute value of the reference signal within the acquisition window.

[0164] The normalized reference signal is compressed using differential pulse code modulation (DPCM).

[0165] Calculate the difference between adjacent sampling points:

[0166] d(n)=a norm (n)-a norm (n-1);

[0167] In the formula, d(n) is the difference value of the nth sampling point; a norm (n) represents the normalized amplitude of the nth sampling point.

[0168] The differential values ​​are quantized and encoded using 4 bits, resulting in a compression ratio of approximately 4:1.

[0169] The compressed data is encapsulated into fixed-length data blocks, with 250 sampling points per second, resulting in approximately 125 bytes after compression.

[0170] The vibration controller operates at a fixed period T. bc =1s broadcast synchronization data frame, the data frame structure is shown in Table 1 below.

[0171] Table 1:

[0172]

[0173] The broadcast uses LoRa spread spectrum modulation, with a center frequency of 470MHz, a bandwidth of 125kHz, a spreading factor of SF=7, and an effective data rate of approximately 5.5kbps.

[0174] After receiving the synchronization data frame, each inspection robot performs the following operations:

[0175] CRC check: Verifies the integrity of the data frame. If the check fails, the frame is discarded and the system waits for the next frame.

[0176] Frame sequence number check: Compare with the sequence number of the previous frame to detect if there are any dropped frames. If there are dropped frames, record them in the log.

[0177] Timestamp alignment: Converts the data frame timestamp from the master node time base to the local time base.

[0178] Reference signal decompression: Decode the DPCM compressed data to recover the normalized reference signal r(t).

[0179] Each inspection robot initiates local elastic wave acquisition within the time window corresponding to the reference signal, based on the data frame timestamp.

[0180] Calculate the start time of data acquisition: T start =T stamp +Δt delay Tstamp For the data frame timestamp, Δt delay This is a preset propagation delay margin, typically set to 50ms.

[0181] In T start The ADC is triggered at a specific time to start acquiring the surface elastic wave response signal s i (t).

[0182] The collection duration is T w After the data collection is completed, it is stored in the local cache.

[0183] According to one aspect of this application, the two-stage phase tracking time measurement is specifically as follows:

[0184] The source reference signal r(t) and the surface response signal s(t) are preprocessed separately:

[0185] (1) DC component removal:

[0186] r'(t) = r(t) - mean{r(t)};

[0187] s'(t) = s(t) - mean{s(t)};

[0188] In the formula, r'(t) and s'(t) are the surface response signals after removing the DC component; mean{·} represents the mean value within the time window.

[0189] (2) Bandpass filtering: A bandpass filter with a center frequency of f0 and a bandwidth of 2×f0 is used to suppress noise outside the frequency band.

[0190] Calculate the normalized cross-correlation function between the preprocessed reference signal and the acquired signal:

[0191] R rs (τ)=[∫r'(t)×s'(t+τ)dt] / [sqrt(∫r' 2 (t)dt)×sqrt(∫s' 2 (t)dt)];

[0192] In the formula, R rs (τ) is the normalized cross-correlation function, with a value range of [-1, 1]; τ is the time delay variable; the integration interval is the effective signal window; sqrt(·) is the square root operation, and the acquired signal is the surface response signal.

[0193] The actual calculation is implemented using Fast Fourier Transform (FFT):

[0194] R rs (τ)=IFFT{FFT{r'(t)}×conj[FFT{s'(t)}]};

[0195] In the formula, FFT{·} is the Fast Fourier Transform; IFFT{·} is the Inverse Fast Fourier Transform; and conj[·] is the complex conjugate operation.

[0196] Within the preset search range [τ] min ,τ max Peak value of cross-correlation function within the search:

[0197] τ coarse =argmax{R rs (τ)},τ∈[τ min ,τ max ];

[0198] In the formula, τ coarse For a rough estimate; argmax represents the independent variable that makes the function reach its maximum value; τ min =L min / v max For the minimum theoretical time; τ max =L max / v min For the maximum theoretical time; L min L max These are the minimum and maximum values ​​of the source-detector distance, respectively; v min v max These are the minimum and maximum estimates of the medium wave velocity, respectively.

[0199] Simultaneously record the peak correlation coefficient:

[0200] ρ max =R rs (τ coarse );

[0201] In the formula, ρ max This is the peak correlation coefficient, used for subsequent signal quality assessment.

[0202] If ρ max <ρ th If the threshold is typically set to 0.3, the quality of the currently acquired signal is determined to be unsatisfactory, and the robot data is marked as invalid and will not be included in subsequent imaging calculations.

[0203] Narrowband filtering is performed on the reference signal and the acquired signal to extract the fundamental frequency component:

[0204] r n (t)=BPF{r'(t),f0,B n};

[0205] s n (t)=BPF{s'(t),f0,B n};

[0206] In the formula, r n (t), s n (t) represent the narrowband filtered reference signal and the acquired signal, respectively; BPF{·,f0,B n} represents the center frequency f0 and bandwidth B n Bandpass filtering operation; B n =0.1×f0 is the bandwidth of the narrowband filter.

[0207] The filter uses a fourth-order Butterworth bandpass filter, which has a linear phase response and avoids introducing additional phase distortion.

[0208] Perform Hilbert transform on the narrowband filtered signal to construct an analytic signal:

[0209] r a (t)=r n (t)+j×H{r n (t)};

[0210] s a (t)=s n (t)+j×H{s n (t)};

[0211] In the formula, r a (t), s a (t) represents the analytic signals (complex form) of the reference signal and the acquired signal, respectively; H{·} is the Hilbert transform operator; j is the imaginary unit.

[0212] The Hilbert transform is implemented in the frequency domain:

[0213] H{x(t)}=IFFT{-j×sign(f)×FFT{x(t)}};

[0214] In the formula, sign(f) is the sign function of frequency, with positive frequencies taking +1, negative frequencies taking -1, and zero frequencies taking 0.

[0215] Extracting instantaneous phase from an analytical signal:

[0216] φ r (t)=atan2{Im[r a (t)],Re[r a (t)]};

[0217] φ s (t)=atan2{Im[s a (t)],Re[s a (t)]};

[0218] In the formula, φ r (t) represents the instantaneous phase of the reference signal; φs (t) represents the instantaneous phase of the acquired signal; atan2{y,x} is the arctangent function in the four quadrants, with a return value range of (-π,π]; Im[·] and Re[·] represent taking the imaginary and real parts of the complex number, respectively.

[0219] The instantaneous phase at the start of the reference signal is taken as the phase reference:

[0220] φ ref =φ r (0);

[0221] In the formula, φ ref Used as a reference phase standard.

[0222] When τ is roughly estimated coarse At the corresponding moment, read the instantaneous phase of the acquired signal:

[0223] φ obs =φ s (τ coarse );

[0224] In the formula, φ obs The observed phase of the acquired signal at a time when it is roughly estimated.

[0225] Calculate the difference between the observed phase and the reference phase:

[0226] Δφ=φ obs -φ ref ;

[0227] In the formula, Δφ is the phase difference, and the unit is radians.

[0228] Normalize the phase difference to the interval (-π, π):

[0229] If Δφ>π, then Δφ=Δφ-2π;

[0230] If Δφ≤-π, then Δφ=Δφ+2π.

[0231] Convert the normalized phase difference into a time correction:

[0232] Δτ fine =Δφ / (2π×f0);

[0233] In the formula, Δτ fine Δφ is the phase refinement time correction, in seconds; Δφ is the phase difference, in radians; f0 is the fundamental frequency of vibration, in Hz.

[0234] The number of complete vibration cycles included in the propagation delay is determined based on a rough estimate:

[0235] N=round(τ coarse ×f0-Δφ / (2π));

[0236] In the formula, N is the number of integer cycles; round(·) is the rounding function.

[0237] By combining the contribution of the whole period and the phase refinement correction, the final high-precision time is obtained:

[0238] τ final =(N+Δφ / (2π)) / f0;

[0239] In the formula, τ final This represents the final high-precision time estimate; N is the number of integer cycles; Δφ is the phase difference; and f0 is the fundamental frequency of vibration.

[0240] The measurement results are packaged into time data records, as shown in Table 2 below.

[0241] Table 2:

[0242]

[0243] In other optional application scenarios, the method proposed in this invention is also applicable to engineering occasions requiring large-area compaction quality testing, such as airport runway compaction, high-speed railway subgrade construction, and port yard foundation treatment. For different fill material types and design compaction requirements, such as cohesive soil, gravel, and weathered rock fill materials, on-site calibration tests can be conducted before construction to establish the corresponding velocity-compaction degree mapping relationship and determine the applicable acceptable threshold K for the project. pass and warning threshold K warn Furthermore, the fundamental vibration frequency and harmonic characteristics may differ for different models and specifications of road rollers. However, the method of this invention, which involves real waveform acquisition at the source end, coarse estimation of cross-correlation, and precise measurement of phase refinement, is still applicable. Only the filter settings and search window range need to be adjusted according to the actual vibration parameters.

[0244] This invention improves the accuracy of time measurement from the sampling interval level (millisecond level) of traditional cross-correlation methods to the phase resolution level (microsecond level) through phase refinement technology, theoretically increasing the accuracy by 1-2 orders of magnitude. It also overcomes the sampling rate limitation.

[0245] By collecting the actual waveform at the vibration source of the road roller as a reference signal, the phase mismatch problem caused by the ideal sine wave assumption is avoided, the sharpness and reliability of the cross-correlation peak are improved, and the influence of signal source distortion is eliminated.

[0246] By combining multi-frequency harmonic constraints and utilizing the information redundancy between different frequency components, the influence of random noise on phase measurements is effectively suppressed, and abnormal propagation paths can be automatically identified. This enhances measurement robustness.

[0247] By analyzing multi-frequency phase consistency residuals and dispersion characteristics, the system can not only eliminate invalid data but also proactively identify loose regions or structural defects in the medium, providing in-depth evidence for engineering quality control. This achieves intelligent defect diagnosis.

[0248] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A method for detecting the quality of a dike using elastic wave tomography, characterized in that, include: Simultaneously acquire the source-end reference signal at the vibration source of the road roller and the surface response signal of the measured area, and establish a unified time reference for the source-end reference signal and the surface response signal; Calculate the cross-correlation function between the source-end reference signal and the surface response signal, and determine the coarse estimate to time based on its peak value; At least one characteristic frequency component is extracted from the source reference signal and the surface response signal, and the phase information of the analytical signal corresponding to each characteristic frequency component is obtained. Based on the phase information of the analytical signal, the phase difference of the surface response signal relative to the source reference signal at the coarsely estimated time is calculated, and the fine time correction is determined by using the phase difference. The coarse estimated arrival time is corrected by using a fine time correction to obtain the final elastic wave arrival time, and an elastic wave velocity distribution image of the measured area is generated by inversion based on multiple final elastic wave arrival times. The source-end reference signal is a vibration acceleration signal; Acquire the source-end reference signal at the vibration source of the road roller, including: Vibration acceleration signals are collected in real time by an acceleration sensor installed at the bearing seat of the vibrating drum of the road roller. The sampling frequency of the acceleration sensor is more than 10 times the fundamental frequency of the vibration source. Obtain the analytical signal phase information corresponding to each characteristic frequency component, including: Narrowband filters whose center frequency corresponds to the characteristic frequency component are used to filter the source reference signal and the surface response signal respectively. Perform a Hilbert transform on the filtered signal to construct an analytic signal in complex form; Based on the real and imaginary parts of the analytic signal, calculate the instantaneous phase that changes with time, and use it as the phase information of the analytic signal; The characteristic frequency components include the fundamental frequency of vibration and at least one higher harmonic component whose frequency is an integer multiple of the fundamental frequency of vibration; Based on the phase information of the analytical signal, the phase difference of the surface response signal relative to the source reference signal at the coarsely estimated time is calculated, including: For each characteristic frequency component, the difference between the observed phase of the surface response signal at the coarsely estimated time and the reference phase of the source reference signal at the initial time is calculated to obtain a set of multi-frequency phase differences corresponding to different frequencies. Determining fine time corrections using phase differences includes: A multi-frequency consistency cost function is constructed to characterize the degree of consistency between the calculated phase corresponding to the candidate time correction and each observation in the multi-frequency phase difference set; Searching for candidate time corrections within a preset range centered on the coarse estimate time, the candidate time correction that minimizes the multi-frequency consistency cost function is determined as the fine time correction.

2. The method according to claim 1, characterized in that, Establish a unified time reference for the source-end reference signal and the surface response signal, including: The acquired source reference signal is subjected to anti-aliasing filtering, and the processed signal is compressed using differential pulse code modulation to obtain a compressed reference waveform. Perform spectrum analysis on the source reference signal to determine the current vibration fundamental frequency, and broadcast synchronization data frames according to a preset period. The synchronization data frames include the compressed reference waveform, the timestamp of the current reference signal, and the current vibration fundamental frequency parameters. It receives synchronous data frames, parses out compressed reference waveforms and decodes them to reconstruct the source reference signal, and maps its time axis to the local time system based on timestamps.

3. The method according to claim 1, characterized in that, Calculate the cross-correlation function between the source-end reference signal and the surface response signal, and determine the coarse estimate based on the peak value of the cross-correlation function, including: The DC components of the source reference signal and the ground response signal are removed respectively, and out-of-band noise is filtered out using a bandpass filter; Calculate the normalized cross-correlation function between the preprocessed source-end reference signal and the surface response signal; Within a preset search time window, the maximum value of the normalized cross-correlation function is searched, and the corresponding time delay is determined as the coarse estimate.

4. The method according to claim 1, characterized in that, The characteristic frequency components include the fundamental vibration frequency; Based on the phase information of the analytical signal, the phase difference of the surface response signal relative to the source reference signal at the coarsely estimated time is calculated, and the fine time correction is determined using the phase difference, including: Obtain the observed phase of the surface response signal at the coarsely estimated time, and the reference phase of the source-end reference signal at the initial time; Calculate the difference between the observed phase and the reference phase, normalize the difference to one vibration period, and obtain the phase difference; The fine time correction is calculated using the formula Δτ=Δφ / 2πf0, where Δτ is the fine time correction, Δφ is the phase difference, and f0 is the fundamental frequency of vibration.

5. The method according to claim 1, characterized in that, Construct a multi-frequency consistency cost function and search candidate time corrections, including: Based on the spectral analysis results of the source-end reference signal, the signal-to-noise ratio of each characteristic frequency component is determined, and weighting factors are assigned to each phase deviation term in the multi-frequency consistency cost function accordingly. Using the vibration period corresponding to the fundamental frequency as the step size, a discrete set of candidate time correction quantities is constructed within a preset range; Calculate the multi-frequency consistency cost function value corresponding to each element in the candidate time correction set in order to find the minimum value of the weighted consistency.

6. The method according to claim 1, characterized in that, The elastic wave velocity distribution image of the measured region is generated based on the time inversion of multiple final elastic wave arrivals, including: Based on the spatial coordinates corresponding to the source reference signal and the surface response signal, the measured area is divided into grids, and a ray path matrix is ​​constructed. Using joint iterative reconstruction technology, the theoretical time is calculated based on the initial velocity model; The residual between the final arrival time of the elastic wave and the theoretical arrival time is calculated, and the residual is back-projected to each grid cell using the ray path matrix to update the velocity model until the iteration converges, thus obtaining the elastic wave velocity distribution image.

Images