A medium change monitoring algorithm based on multiple scattered wave interference

Abstract

The application aims to provide a medium change monitoring algorithm based on multiple scattering wave interference, comprising the following steps: signal acquisition, signal segmentation, similar information acquisition, time delay calculation; medium change monitoring. Compared with the Stretching method, the application can realize the evaluation of nonlinear time delay and local medium change, and does not need strict time synchronization between signals; compared with the WCC method, the time resolution is greatly improved, and the period jump is effectively suppressed; compared with the traditional DTW method, the signal is segmented and processed, and a new bending window is used, so that the width of the bending window after truncation processing can be reduced by several times, thereby greatly improving the calculation efficiency. The horizontal offset is calculated based on the waveform characteristics, the signal is truncated, the occurrence of ill-conditioned matching is effectively suppressed, the constraint condition value is reasonably selected, the calculation efficiency is not reduced, the calculation accuracy is effectively improved, the influence of waveform distortion is weakened, and the calculation accuracy is greatly improved.

Description

TECHNICAL FIELD

[0001] The present application relates to an elastic wave signal processing method, in particular to a medium change monitoring method. BACKGROUND

[0002] Monitoring medium changes with high sensitivity and high precision can provide technical support for volcano early warning, earthquake research, material detection, and engineering safety evaluation. In a propagation medium with complex geometry or heterogeneous material, elastic wave energy quickly converts from a direct wave field to a multiple scattering wave field. Compared to the direct wave field, multiple scattering waves sample the medium multiple times during their long-range propagation, so they are sensitive to changes in medium properties. By interfering with the acoustic signals corresponding to the multiple scattering waves, high sensitivity and high precision monitoring of medium changes can be achieved, so it is widely used in earthquake research and ultrasonic nondestructive testing.

[0003] The evaluation focus of current multiple scattering wave interference analysis methods is on the micro changes of the elastic wave velocity of the propagation medium. Traditional analysis methods include the Windowed Cross-Correlation (WCC) method and the Stretching method. The WCC method is simple to implement, but due to differences between its theoretical assumptions and physical scenarios and the lack of cycle skipping suppression in the algorithm itself, it results in insufficient precision and stability in medium change evaluation. The Stretching method has improved stability and evaluation precision, but it is limited by the linear time delay assumption and the high-precision time synchronization requirement for signal collection equipment, resulting in an algorithm that is not suitable for local velocity change evaluation and increased application costs. In summary, the existing multiple scattering wave interference algorithm has unsatisfactory calculation precision, and its reliability cannot be guaranteed when faced with nonlinear time delay caused by local velocity changes, low signal-to-noise ratio, or cycle skipping caused by excessively large time delay.

[0004] The Dynamic Time Warping (DTW) method is a time-domain acoustic nonlinear optimization method developed in the field of speech recognition. It can compensate for differences in audio signals by using a point-by-point matching approach. The DTW method has nonlinear analysis potential and some ability to suppress cycle skipping, but it cannot be directly applied to multiple scattering wave interference analysis due to the huge differences between application scenarios and propagation processes.

[0005] In summary, it is necessary to develop a new medium change monitoring algorithm based on multiple scattering wave interference by absorbing the point-by-point matching processing idea of the DTW method to meet the needs of medium change monitoring in the presence of waveform distortion, low signal-to-noise ratio, and nonlinear time delay scenarios. SUMMARY

[0006] The purpose of this invention is to provide a high-precision and high-efficiency algorithm for monitoring medium changes based on multiple scattering wave interferometry, which has the ability to achieve ill-condition matching and periodic jump suppression.

[0007] The objective of this invention is achieved as follows:

[0008] This invention discloses a medium change monitoring algorithm based on multiple scattering wave interferometry, characterized by:

[0009] (1) Signal acquisition;

[0010] (2) Signal segmentation;

[0011] (3) Obtain similarity information and, based on cross-correlation processing, obtain the horizontal offset α between segmented signals. i Based on α i Perform truncation to obtain similar information;

[0012] (4) Delay calculation: The error matrix D within the curved window is calculated based on the Euclidean distance. The cumulative error matrix C is calculated from the error matrix D and the constraint condition b. The optimal regularization path P is determined from the cumulative error matrix C and the constraint condition b. The average offset β is then obtained from the optimal regularization path P. i In order to obtain the time delay τ i ;

[0013] (5) Monitoring of medium changes: Calculate the decorrelation coefficient of the two signals under the optimal regularized path, based on the time delay τ in step (4). i The velocity change is calculated with the elapsed time, and the signals before and after the change of the medium are continuously received and acquired. The monitoring of the change of the medium is realized by steps (1)-(5).

[0014] The present invention may also include:

[0015] 1. The specific steps of signal acquisition in step (1) are as follows:

[0016] The time-domain signals of multiple scattered waves before and after the change in medium are obtained and denoted as u0(t) and u1(t).

[0017] 2. The specific steps of signal segmentation in step (2) are as follows:

[0018] Based on the total time window t all The window length t of the moving time window mw and the step distance t step Obtain k segmented signals u 0，i (t) and u 1，i (t), 1≤i≤k.

[0019] 3. The specific method for obtaining similar information in step (3) is as follows:

[0020] Horizontal offset acquisition: for segmented signal u0，i (t) and u 1，i (t) Perform cross-correlation, calculate the maximum value of the cross-correlation function and the position corresponding to the maximum value. This position is the piecewise signal u. 0，i (t) and u 1，i Horizontal offset α between (t) i ;

[0021] Truncation: The offset α obtained from the calculation i Align the segmented signals, and align the segmented signals u 0，i (t) and u 1，i (t) is truncated to the portion that almost overlaps within the shift time window; the segmented signal u after truncation processing. 0，i (t) and u 1，i (t) has a length of m i =t mw -|α i |

[0022] 4. The specific delay calculation in step (4) is as follows:

[0023] Calculate the error matrix: Construct a matrix of size m i ×m i The error matrix D is obtained, with all elements inside the matrix being positive infinity. Then, the error matrix D within the curved window is recalculated based on the Euclidean distance formula, where D(j, l) = (u 0，i (j)-u 1，i (l)) 2 ;

[0024] Calculate the cumulative error matrix: Constraint b will control the calculation of the cumulative error matrix C within the curved window. When b = 1, C(j, l) = D(j, l) + min{C(j-1, l), C(j-1, l-1), C(j, l-1)}; when b ≥ 2, the cumulative error matrix C is calculated using the following formula.

[0025]

[0026] Finding the optimal regularization path: In the cumulative error matrix C, determine the movement path based on the constraint values, and find the optimal regularization path P = [P1, P2, ... P]. K ], m i ≤K≤2m i -1; from P K =(m i m i Starting from C(m) i -b,m i -b+1), C(m) i -1,m i -1) and C(m i-b+1,m i -b) The row and column value where the minimum of the three values ​​is located is the next point P. K-a = The value of (row, column), a = 1 or a = b is determined by the path, and then according to P K-a Find the next point P K-2a Repeat this operation until P1(1,1) is found;

[0027] Time delay estimation: For each point P in the regular path P c =(i c j c ) to perform (i c -j c The operation is performed where 1 ≤ c ≤ K. The average of the results is the average offset β. i In addition, the offset α calculated in step (3) i That is, the time delay τ within the moving time window. i .

[0028] 5. Step (5) Monitoring of medium changes specifically involves:

[0029] Decorrelation coefficient calculation: based on u 0，i (i c ) and u 0，i (j c ), 1≤c≤K, calculate the cross-correlation coefficient C(i) between signals, and remove the correlation coefficient D. i =1-C i ;

[0030] Velocity change calculation: Based on steps (3) and (4), all segmented signals u are calculated and obtained. 0，i (t) and u 1，i The time delay τ of (t) i 1≤i≤k and the step distance are used to calculate the velocity change Δv i That is, Δv i =(τ i+1 -τ i ) / t step , 1≤i≤k-1.

[0031] The advantages of this invention are:

[0032] 1. Compared with the Stretching method, it can evaluate nonlinear time delay and local medium changes, and does not require strict time synchronization between signals;

[0033] 2. Compared with the WCC method, the time resolution is greatly improved due to the use of point-by-point matching, and the period jump is effectively suppressed;

[0034] 3. Compared to the traditional DTW method, this approach, tailored to the characteristics of multiple scattered wave signals, segments the signal and utilizes a novel curved window. The width of the curved window is significantly reduced after truncation, greatly improving computational efficiency. The horizontal offset is calculated based on waveform characteristics, and signal truncation effectively suppresses ill-conditioned matching. Furthermore, by appropriately selecting constraint values, computational accuracy is improved without sacrificing efficiency, mitigating the impact of waveform distortion and significantly enhancing computational precision. Attached Figure Description

[0035] Figure 1 This is a flowchart of the present invention;

[0036] Figure 2 This is a schematic diagram of the curved window of the present invention;

[0037] Figure 3 This is a calculation diagram of the optimal regularized path of the present invention;

[0038] Figure 4 This is a simulation signal diagram of the present invention;

[0039] Figure 5 This is a diagram of the segmented signal before truncation processing;

[0040] Figure 6 This is a diagram of the segmented signal after truncation.

[0041] Figure 7 The nonlinear time delay of this invention is 10. -2 The estimation results of four methods are shown in the figure.

[0042] Figure 8 The nonlinear time delay of this invention is 10. -2 CDF plots for the four methods. Detailed Implementation

[0043] The invention will now be described in more detail with reference to the accompanying drawings:

[0044] Combination Figures 1-8 The present invention specifically includes the following steps:

[0045] Step 1 Signal Acquisition: This invention uses a simulated signal u0(t) and performs linear or nonlinear stretching on the simulated signal to construct the signal u1(t) after the medium change;

[0046] Step 2: Signal Segmentation. Signals u0(t) and u1(t) are segmented. This is because multiple scattered waves are long-duration signals, and interpolation is often used when the medium change is very weak, leading to excessively long signal lengths and significantly reduced computational efficiency. Furthermore, it reduces the probability of ill-conditioned matching; the maximum allowable offset between segments is limited to the length of the segmented signal. First, the total time window t is determined. all This ensures that the multiple scattered wave signal dominates within the total time window, and that the signal-to-noise ratio is high. The window length t of the moving time window is set. mw The moving time window refers to a window with a specified length t. mw The time window, included within the total time window, is calculated according to the set movement step t. step The signal is segmented by sliding the moving time window along the time axis. It's important to note that the moving time window needs to contain sufficient information, meaning the waveforms must be sufficiently complex. However, if the window length is too large, the calculation speed will decrease; therefore, a reasonable window length needs to be set. Simultaneously, to ensure the continuity of the analysis results, the moving step size needs to be small, allowing for signal overlap within the moving time window. Assume a total time window t... all =[t1-t w , t k +t w The window length of the moving time window is t. mw =2t w The center position of the sliding time window is t. i 1≤i≤k, where k is the number of segmented signals. This invention makes the waveform change within the moving time window include five cycles, and the moving step size is set to 1 / 15 of the window length, i.e., t step =t i -t i-1 =1 / 15×t mw Obtain the segmented signal u 0，i (t) and u 1，i (t), 1≤i≤k;

[0047] Step 3: Obtain similar information:

[0048] 3.1: Horizontal Offset Acquisition. Based on cross-correlation processing and considering its waveform characteristics, the horizontal offset between segmented signals is calculated. Assume the time delay variable t... s In [-2t] w 2t w Within the range, obtain the i-th segment signal u. 0，i (t) and u 1，i (t), t∈[t i -t w , t i +t w The cross-correlation coefficient R(t)s The formula is shown below, where the time delay variable t corresponds to the maximum value of the cross-correlation coefficient. s That is, the segmented signal u 0，i (t) and u 1，i Horizontal offset α of (t) i , where k refers to the number of segmented signals;

[0049]

[0050] When α i ≥t mw When i≥1, consider the entire signal, set the number of segments k to 1, and calculate the segmented signal u. 0，i (t) and u 1，i The time delay α of (t) i , 1≤i≤k.

[0051] 3.2: Truncation Processing. Based on the segmented signal u obtained in step 3.1 0，i (t) and u 1，i Horizontal offset α between (t) i The segmented signals are truncated to ensure that the microstructure information of the propagation medium contained in the segmented signals within the time window is basically consistent, thereby reducing the error caused by the determinism of DTW boundary conditions and effectively suppressing ill-conditioned matching. If the time delay α... i ≥0, then u 0，i The i-th time window of (t) is adjusted to [t]. i -t w , t i +t w -α i ], i≥1, u 1，i The time window of the i-th segment of (t) is adjusted to [t]. i -t w +α i , t i +t w Conversely, if the time delay α i <0 then u 0，i The time window of the i-th segment of (t) is adjusted to [t]. i -t w +|α i |,t i +t w ], u 1，i The time window of the i-th segment of (t) is adjusted to [t]. i -t wi , t i +t w -|α i |]. After truncation, the segmented signal u 0，i (t) and u 1，iThe length of (t) is m. i =t mw -|α i |;

[0052] Step 4: Delay Calculation

[0053] 4.1: Error Matrix Calculation. To reduce computational complexity and improve efficiency, only the error matrix within the curved window is calculated for the i-th segment of the truncated signal. This invention proposes a novel curved window, as shown in the attached figure. Figure 2 As shown, by combining the high efficiency of Sakoe-Chiba band with the high precision of Itakura-Paralleogram, Itakura-Paralleogram is used at the beginning and end of the sequence pair, while Sakoe-Chiba band is used in the middle. This not only further improves computational efficiency but also ensures the reliability of the results. Figure 2 middle This represents rounding down to the nearest integer, i.e., the largest integer not exceeding x. First, define several parameters... k2 = 2k1 + 1, The entire window is divided into 7 regions along the vertical axis by the range of the gray box, which are (1) 1≤x≤k1+1; (2) k1+2≤x≤k2; (3) k2+1≤x≤k2+k3; (4) k2+k3+1≤x≤m i -(k2+k3); (5)m i -(k2+k3)+1≤x≤m i -k2; (6)m i -k²+1≤x≤m i -k1-1; (7)m i -k1≤x≤m i The range of the gray box on the vertical axis is calculated according to the formula in the figure based on the divided regions. For example, for a 2000×2000 sequence, performing 100 Monte Carlo calculations and then taking the average time, with the width of the curved window set to 100, the traditional DTW algorithm takes 132.3ms, the Sakoe-Chuba band takes 55.8ms, the Itakura-Paralleogram takes 69.7ms, while this method takes the shortest time of 54.4ms, representing an improvement of 58.9% compared to the traditional DTW algorithm.

[0054] First, construct m i ×m iThe identity matrix is ​​1 ≤ i ≤ k. This identity matrix is ​​then multiplied by positive infinity, making all elements within the matrix positive infinity. This is denoted as the error matrix D. Typically, the width of the curved window is set to 10% of the segmented signal length. However, after truncation, the width of the curved window can be reduced by a factor of two. A smaller width results in higher computational efficiency, but too small a width can lead to the inability to correctly find the optimal regularization path. Therefore, it can be adjusted according to actual needs. In this invention, the width n of the curved window is set to 5% of the segmented signal length, i.e., n = m. i ×0.05. Recalculate the error matrix D within the curved window, using Euclidean distance to calculate D(j,l)=(u 0，i (j)-u 1，i (l)) 2 , refers to u 0，i (t) The j-th element and u 1，i The square of the difference of the l-th element of (t) is the value of the j-th row and l-th column of the error matrix D;

[0055] 4.2: Calculate the cumulative error matrix. The slope constraint is used to limit the variation of the regularized path. The larger the constraint value, the smoother the matching path, which can reduce the impact of waveform distortion and isolated point interference, but the computational complexity will also increase. Construct a matrix of the same size as the error matrix D, denoted as the cumulative error matrix C, with all elements being positive infinity. Recalculate the cumulative error matrix C within the curved window based on the set constraint value b and the error matrix D calculated in step 4.1. In this invention, the constraint value is set to 2, and the calculation formula is shown below;

[0056] C 1，1 =D 1，1

[0057] C i，j =D i，j +min{C i-1，j-1 C i-1，j C i，j-1}i,j≤n

[0058]

[0059] 4.3: Finding the optimal normalization path. In the cumulative error matrix C, find the optimal normalization path P = [P1, P2, ... P...]. K ], m i ≤K≤2m i -1, path P represents u 0，i (t) and u 1，i The matching relationships between points in (t). The following three conditions must be met when calculating the regularized path:

[0060] ① Monotonicity: If Pk = (i, j), P k+1 = (i′, j′), then i′≥i, j′≥j, the purpose of which is to ensure that the points on the path must be monotonically increased over time;

[0061] ② Continuity: If P k = (i, j), P k+1 = (i′, j′), then i′≥i+1, j′≥j+1. The purpose is to restrict the points on the path to be adjacent points and not allow any jumps.

[0062] ③ Boundary condition determinism: P1 = (1, 1), P K =(m i m i );

[0063] Based on the set value of constraint b, from P K =(m i m i Find P1 = (1, 1) by reverse direction. Given P K =(m i m i ), determine C(m) i -b,m i -b+1), C(m) i -1,m i -1) and C(m i -b+1,m i -b) The row and column position of the minimum value of the three, if position X = (m i-1 m i-1 If we move one unit to the bottom left corner, then let a = 1; otherwise, let a = b. Thus, P... K-a = (x, y), where x is the row containing the minimum value and y is the column containing the minimum value. (For example, see attached...) Figure 3 Substituting b = 2, then C(m) i -2,m i -1), C(m) i -1,m i -1) and C(m i -1,m i -2) The minimum value among the three is C(m) i -2,m i -1), not moving one unit along the bottom left corner, therefore a = b = 2, i.e., P K-2 =(m i -2,m i -1). Based on the constraint b value and P... K-2 The coordinates of P can determine the position of P. K-1 =(m i -1,m i). And so on, until P1 = (1, 1);

[0064] 4.4: Time Delay Estimation. Optimal Regularized Path P in P c =(i c j c ) represents u 0，i (i c ) and u 1，i (j c ) match, therefore, (i c -j c ) represents the i-th segment signal u 0，i (t) and u 1，i (t) at time point i c The offset at time. For each point P in the optimal regular path. j Repeat the above operation for 1≤j≤K to obtain the segmented signal u. 1，i (t) relative to the segmented signal u 0，i (t) The offset at all time points, and the average offset β is calculated by averaging the results. i , i.e. β i =[(i1-j1)+,...,+(i K -j K )] / K, plus the offset α calculated in step 3 i That is, the time delay τ within the moving time window. i , τ i =β i +α i ;

[0065] Step 5: Monitoring changes in the medium

[0066] 5.1: Calculate the correlation coefficient; set time point i c j c Substitute the segmented signals u respectively 0，i (t) and u 1，i (t), i.e., u 0，i (i c ) and u 0，i (j c ), where 1≤c≤K. The signal u is obtained using the following formula. 0，i (i c ) and u 0，i (j c The cross-correlation coefficient C) i , where μ 0，i and σ 0，i Indicates signal u 0，i (i c The mean and variance of μ; while μ 1，i and σ 1，i Indicates signal u1，i (j c The mean and variance of ); according to D i =1-C i Calculate the decorrelation coefficient D i ;

[0067]

[0068] 5.2: Velocity Change Calculation: Based on steps 3 and 4, calculate and obtain all segmented signals u. 0，i (t) and u 1，i The time delay τ of (t) i And 1≤i≤k and the movement step distance are used to calculate the speed change Δv i That is, Δv i =(τ i+1 -τ i ) / t step , 1≤i≤k-1;

[0069] The following simulation examples will provide further explanation:

[0070] Figure 4 The black dashed line and gray solid line represent the signals u0(t) and u1(t) before and after the change in medium, respectively. For ease of calculation, time is expressed in time points, and the conversion relationship is time in seconds multiplied by the time sampling rate. The total time window t is set. all = [501, 2500], the window length t of the moving time window mw =150, step size t step =10, then the number of segmented signals k is 184, that is, the number of segmented signals u is obtained. 0，1 (t) and u 1，1 (t), u 0，2 (t) and u 1，2 (t), ..., u 0，184 (t) and u 1，184 (t).

[0071] The following section uses the 80th segment signal as an example to explain steps 3 and 4 in detail. Figure 5 As shown, u at this time 0，80 (t)=u0(1291:1440) and u 1，80 (t) = u1(1291:1440), and its horizontal offset α is calculated from step 3.1. 80 =-2, because α 80 <0, u 0，80 (t) is truncated to u 0，80 (t)=u0(1291:1438,u 1，80 (t) is truncated to u 1，80(t) = u1(1292:1440), the segmented signal after truncation is as follows Figure 6 As shown, compared to Figure 5 It can be observed that the segmented signals are highly similar after truncation, and the horizontal offset between the segmented signals is basically compensated. Therefore, the setting of the curved window in the subsequent calculation of the error matrix and cumulative error matrix can be reduced significantly, greatly improving computational efficiency. The signal length m after truncation is... 80 =t mw -|α 80 |=148. Since the time delay between segmented signals is very small, interpolation is performed on the segmented signals with an interpolation interval of 0.1. After interpolation, the segmented signal u... 0，80 (t) and u 1，80 The length of (t) is m 80 =1471. From m 80 Determine the length n of the curved window, and round down to the nearest integer. Construct m 80 ×m 80 The error matrix D of positive infinity is used to recalculate the curved window using the Euclidean distance formula, as shown in the attached figure. Figure 2 The error matrix within the gray box, for example Construct a matrix of the same size as the error matrix D, with internal element values ​​of positive infinity, denoted as the cumulative error matrix C. Set the slope constraint b to 2. Based on the error matrix D and the slope constraint b, recalculate the cumulative error matrix C within the bending window according to the formula in step 3.2, for example, from k1 = C. 99，148 +2×D 100，149 +D 100，150 =0.03447, Hall 2 = C 98，149 +2×D 99，150 +D 10，0150 =0.03896 and k3=C 99，149 +2×D 10，0150 The minimum value of the three factors (0.03672, 0.03672, 0.03447) determines the value of C(100, 150), i.e., C(100, 150) = 0.03447. Similarly, calculate the elements of the remaining cumulative error matrix C within the curved window. The optimal regular path P is determined by the cumulative error matrix C and the slope constraint, given that P is known. K =(m 80 m 80Let P1 = (1, 1) and P1 = (1, 1). The value of K is not known beforehand and can be considered an unknown with a fixed value. The value of K can only be determined after the optimal regularization path is determined. The path position is determined based on the row and column positions of the minimum values ​​of C(1469, 1470) = 0.008347, C(1470, 1470) = 0.008332, and C(1470, 1469) = 0.008261. Since C(1470, 1469) is the minimum and not in a symmetrical position (i.e., the minimum value is not C(1470, 1470), a = b = 2. Therefore, P1 = (1, 1). K-a =P K-2 = (1470, 1469), based on the slope constraint value b and P K-2 The coordinates determine its path direction by first moving downwards once, then moving along the symmetrical direction, thus determining P. K-1 = (1471, 1470). Continue in this manner until P1 = (1, 1). Then, we obtain K = 1713, and thus P = [P1, ..., P...]. c , ..., P 1713 ], where P c =(i c j c ), 1≤c≤K and 1≤i c ≤m 80 For example, for P K-2 = (1470, 1469) is used as the row coordinate, and the ordinate is subtracted from the ordinate to obtain the offset y. K-2 = (1470-1469) = 1. Similarly, for each point P in the regular path P... j For 1≤j≤K, perform the same process to obtain offsets y1, y2, ..., y K Finally, regarding the y1, y2, ..., y K Perform mean averaging and calculate the average offset β. 80 = (y1+y2+,...,+y K ) / K = -0.4776, compared with the horizontal offset α calculated in step 3. 80 Add the two segments together to obtain the time delay τ between the 80th segment and the 2-3 segments. 80 , τ 80 =β 80 +α 80 = -2.4776. Similarly, process the piecewise signal u in the same way. 0，1 (t) and u 1，1 (t), u 0，2 (t) and u 1，2 (t), ..., u 0，184 (t) and u 1，184 (t), obtain the time delays τ1, τ2, ..., τ between the segmented signals.184 .

[0072] Calculate the cross-correlation coefficient C of signals Sig1 and Sig2 under the optimal normalized path using the formula in step 5.1. 80 =0.9767, where Sig1(c) = u 0，80 (i c Sig2(c) = u 1，80 (j c ), that is, the decorrelation coefficient D 80 =0.0233. Similarly, obtain the decorrelation coefficients D1, ..., D for all segmented signals. 184 From the formula Δv i =(τ i+1 -τ i ) / t step 1≤i≤k-1, calculate the velocity change.

[0073] Finally, to analyze the performance of the medium change monitoring algorithm based on multiple scattering waves of the present invention, the performance differences between the method mentioned in the present invention and the WCC method, Stretching method and traditional DTW algorithm are considered under two scenarios: linear and nonlinear time delay changes, and under different time delay orders. Figure 6 The nonlinear time delay is on the order of 10. -2 The graphs showing the estimation results of the four methods reveal that only this method can accurately and reliably estimate the nonlinear time delay at this order of magnitude. While WCC, when fitted, can roughly reflect the trend of nonlinear time delay changes, it cannot accurately estimate the true value. The result of the traditional DTW method approaches zero, indicating that traditional DTW has reached its limit and exhibits insufficient resolution. The Stretching method cannot estimate the nonlinear time delay changes; its result shows a linear trend. Figure 7 This further demonstrates intuitively that this method has the smallest error and the highest calculation accuracy. Figure 7The CDF plot shows the probability distribution function, with the horizontal axis representing the estimation error and the vertical axis representing the cumulative distribution function, indicating the probability of occurrence under a certain estimation error. A smaller horizontal axis indicates a smaller estimation error and higher computational accuracy. Furthermore, combining the root mean square error (RMSE) of the four methods under different orders of magnitude of linear and nonlinear time delays in Table 1, a smaller RMSE indicates better method performance. It can be observed that the method mentioned in this invention performs slightly worse than the Stretching method for linear time delay estimation, while exhibiting the best performance in nonlinear time delay estimation, achieving high-precision estimation of time delay changes. Moreover, regardless of whether the time delay is linear or nonlinear, the RMSE of this method is an order of magnitude lower than that of the traditional DTW algorithm, meaning its performance is superior to the traditional DTW method. In summary, the medium change monitoring algorithm based on multiple scattering waves mentioned in this paper is applicable not only to linear but also to nonlinear medium changes, providing a more comprehensive and accurate description of medium changes.

[0074] Table 1. RMSE results of four linear and nonlinear methods

[0075]

[0076]

Claims

1. A medium change monitoring algorithm based on multiple scattered wave interference, characterized by: (1) Signal acquisition; acquire the time domain signals of multiple scattered waves before and after the change of medium, denoted as and ; (2) signal segmentation; obtain k segmented signals and , , the window length of the moving time window when the signal is segmented is ; (3) Based on cross-correlation processing, the horizontal offset between the segmented signals is obtained , based on truncation processing is performed to obtain similar information; correlating the segmented signals and calculating the maximum value of the cross-correlation function and the position corresponding to the maximum value, which is the horizontal offset between the segmented signals and ;​ the calculated offset aligning the segmented signal, segmenting the signal and truncating to portions that almost overlap within the moving time window, the segmented signal after truncation processing and both have a length of ; (4) Definition , , , the curved window is constructed, the Itakura-Paralleogram is used at the beginning and end of the sequence pair, the Sakoe-Chiba band is used in the middle part, n is the width of the curved window, the curved window divides the horizontal coordinate into 7 regions, which are (1) ; (2) ; (3) ; (4) ; (5) ; (6) ; (7) ; calculating the error matrix D in the curved window based on the Euclidean distance: wherein the error matrix D has a size of ; calculating the cumulative error matrix C from the error matrix D and the constraint condition b: If , ; If , ; determining the optimal regularization path P from the cumulative error matrix C and the constraint condition b: ； From start, the row and column value where the minimum of , and is located is taken as the value of the next point, or is determined by the path, and the next point is found according to , and this operation is repeated until is found; For each point in the best regularized path P perform an operation , take the average of the results, which is the average offset , plus the offset calculated in step (3) , to obtain the latency under this moving time window ; (5) Medium change monitoring, two signals under the best regular path With Calculate the cross-correlation coefficient between signals , and then get the decorrelation coefficient ; Based on the latency of step (4) With the elapsed time calculating its speed change : continuously receiving signals before and after the medium change, and realizing medium change monitoring by steps (1)-(5).

2. The algorithm for monitoring changes in medium based on multiple scattered wave interference according to claim 1, characterized in that: The signal segmentation in step (2) is specifically: based on the total time window , the window length of the moving time window and the moving step obtaining k segmented signals and .

3. The algorithm for monitoring changes in a medium based on multiple scattered wave interference according to claim 1, characterized in that: In step (4), an error matrix D of size is constructed, with all elements in the matrix being positive infinity. Then, the error matrix D in the bending window is recalculated based on the Euclidean distance formula, .

Images