Improved CEEMD-based correction method for low-frequency shock response spectrum in large dual-wave shock tester devices

2.1. SRS calculation

The SRS, which represents the maximum response values of a series of SDOF systems under a specific shock load as a function of the system’s natural frequency, is widely used in shock-related research [29].

Fig. 1SDOF systems

If a specific shock load is applied to the SDOF, as shown in Fig. 1, then the system’s mathematical model can be described by the following second-order differential equation:

where $m$, $c$, and $k$ represent the mass, damping coefficient, and stiffness of the SDOF system, respectively; $y\text{'}\text{'}$ represents the input shock excitation acceleration signal, and $x\text{'}\text{'}$ represents the response of the system to the shock signal.

The following can be derived from Eq. (1) by defining the relative displacement of the mass to the base as $z=x-y$:

The undamped intrinsic frequency is:

where ${\omega }_n$ is the intrinsic frequency in rad/s. The system damping ratio is:

Importantly, $\xi$ is commonly represented by the amplification factor $Q$, where $Q=1/\left(2\xi \right)$ is referred to as the “quality factor”.

Substituting Eq. (3) and (4) into Eq. (2) yields the equation of motion for the relative response:

Solving the differential Eq. (5) yields the relationship between the shock response’s peak and the SDOF system’s natural frequency, enabling conversion of the time-domain shock acceleration signal into a frequency-domain SRS.

2.2. Principles of CEEMD algorithm

The CEEMD algorithm enhances EEMD by “completely” ensembling white noise to more accurately mitigate the impact of noise on decomposition results. By adding a pair of auxiliary white noise signals with opposite signs to the signal and averaging the noise to nearly zero, CEEMD effectively avoids the mode-mixing problem of EEMD and achieves higher decomposition accuracy [30].

Assuming $a\left(t\right)$ is the signal to be decomposed, with a length of $N$, the CEEMD decomposition process $a\left(t\right)$ unfolds as follows.

(1) Determine the number of iterations, then add a pair of white noise sequences with equal amplitude and opposite signs, $n_i^{+}\left(t\right)$ and $n_i^{-}\left(t\right)$, to the original signal $a\left(t\right)$. This generates two sets of noise-contaminated signals, $a_i^{+}\left(t\right)$ and $a_i^{-}\left(t\right)$:

(2) Decompose $a\left(t\right)$ with CEEMD to obtain two sets of IMF components, $\mathrm{I}\mathrm{M}{\mathrm{F}}_{si}^{+}$ and $\mathrm{I}\mathrm{M}{\mathrm{F}}_{si}^{-}$. Here, $\mathrm{I}\mathrm{M}{\mathrm{F}}_{si}^{+}$ represents the ensemble mean of the IMFs after adding positive noise; $\mathrm{I}\mathrm{M}{\mathrm{F}}_{si}^{-}$ represents the ensemble mean of the IMFs after adding negative noise.

(3) Calculate the mean values of $\mathrm{I}\mathrm{M}{\mathrm{F}}_{si}^{+}$ and $\mathrm{I}\mathrm{M}{\mathrm{F}}_{si}^{-}$ to determine $\mathrm{I}\mathrm{M}{\mathrm{F}}_{si}$:

(4) Obtain the mean values for all $\mathrm{I}\mathrm{M}{\mathrm{F}}_{si}$ components after $m$ decompositions to determine the final decomposition result.

(5) Decompose the signal into all IMF components and the residual, which represents the non-decomposable trend components of the signal.

The final decomposition results consist of all IMF components and the residual. The IMF components represent the signal’s different frequency or scale components, while the residual corresponds to the low-frequency trend or non-periodic elements of the signal.

2.3. Principles of cross-correlation analysis

This study employs cross-correlation analysis as a critical criterion for selecting effective IMF components. Cross-correlation analysis is a method used to measure the similarity between two sets of signals [31-33]. It is commonly applied in signal processing, image analysis, and time series analysis to detect similarities, time lags, and signal alignments, among other tasks.

The cross-correlation function $R_{xy}\left(\tau \right)$ of two signals $x\left(t\right)$ and $y\left(t\right)$ is defined as follows:

where $\tau$ represents the delay time, indicating the time shift of one signal relative to the other.

The cross-correlation coefficient, obtained by normalizing the cross-correlation function, provides a convenient way to compare the similarity between two signals. Assuming two discrete signals $x$ and $y$, each of length $N$, the calculation process for the cross-correlation coefficient unfolds over five steps.

(1) Calculate the mean value. First, compute the mean values of signals $x$ and $y$:

(2) Remove the mean value. Subtract the mean value of each signal from the signal itself to obtain a signal with a 0 mean:

(3) Calculate the cross-correlation value. For any time delay $\tau$, the unnormalized value of cross-correlation is:

This represents the correlation between signal $x$ and signal $y$ at time delay $\tau$.

(4) Calculate auto-correlation values. The auto-correlation values, defined as:

are employed to normalize the cross-correlation values to a range between [–1, 1].

(5) Calculate the normalized cross-correlation coefficient:

where the cross-correlation coefficient ${\rho }_{xy}\left(\tau \right)$ falls within the range [−1, 1]. A value of $\left|{\rho }_{xy}\left(\tau \right)\right|$ closer to 1 indicates a higher degree of correlation between the two signals, while a value closer to 0 indicates a lower degree of correlation. After a series of comparative tests, a cross-correlation coefficient threshold of 0.1 is appropriate for selecting effective IMF components.

2.4. Principles of K-means clustering

K-means clustering is employed to further validate the rationality of the selected IMF components. This method groups data based on distance, initially by randomly selecting $K$ points as cluster centers, then assigning each data point to the nearest such center to form clusters. After assignment, the cluster centers are recalculated based on the current data points. This iterative process continues until one of the following termination criteria is met: no data points are reassigned, the cluster centers no longer change, or the sum of squared errors reaches a local minimum [34].

The specific steps of K-means clustering are as follows:

(1) Normalize data points and filter out outliers.

(2) Select $K$ initial center points, denoted as $z_1,z_2,\cdot \cdot \cdot ,z_k$.

(3) Define the loss function as:

where $y_i$ represents the $i$-th data point, $c_i$ denotes the cluster assigned to $y_i$, $z_{c_i}$ is the center point of the corresponding cluster, and $M$ is the total number of data points.

(4) Let $t=\mathrm{0,1},2\cdots$ represent the iteration step. During the $t$-th iteration, $c_i^t$ is the cluster assigned to the $i$-th point, and $z_i^t$ represents the $K$ center points. Assign each data point $x_i$ to the nearest center point:

Recalculate the center points to minimize the loss function, and repeat these steps until a termination criterion is met.

3.1. Shock environment measurement process

As shown in Fig. 2, this study employs a 12-t vertical dual-wave shock tester device developed by our institute. The improved CEEMD method eliminates trend term errors in acceleration signals, thereby restoring the actual SRS.

Fig. 2Vertical dual-wave shock tester device

During shock testing, the load platform is lifted to a preset height along the guide rod using the control system. Once the lifting clamping device is released, the platform slides down along the guide rods under the influence of gravity, impacting the waveform generator and producing a positive half-sine acceleration shock wave. Following the shock, a hydraulic damping cylinder buffers and restricts the platform’s motion as it rebounds upward. This creates a negative half-sine shock wave. The full sequence ultimately forms a dual-wave acceleration waveform.

A piezoelectric acceleration sensor is used to measure the shock load during testing. The sensor has a frequency response range of 0.1 Hz-12,600 Hz, a sensitivity of 1 pC/ms², and a maximum measurement range of 60,000 ms^-2. It is configured with a sampling frequency of 1280 Hz and mounted on the surface of the loading platform using a mechanical filter.

3.2. Analysis of original SRS

Shock tests were conducted under four different conditions using the dual-wave shock tester device. For all conditions (I-IV), the acceleration peaked at approximately 0.08 s with maximum values of 2149.1 m/s², 2491.7 m/s², 2438.5 m/s², and 2374.4 m/s², respectively. Taking Condition I as an example, Fig. 3 shows the acceleration signal and its time-domain integral curve. The acceleration signal does not immediately return to zero after the shock ends but continues to oscillate for some time. Consequently, the integrated velocity derived from the acceleration is non-zero; the integrated displacement increases over time, with the curve exhibiting significant drift. This indicates that the acceleration signal is affected by the trend term error.

The SRS of the acceleration signal containing trend term error under Condition I was calculated using the method described in Section 2.1, as shown in Fig. 4. According to our references [35], the standard slope of the low-frequency band in the four-coordinate pseudo-velocity shock spectrum is 6 dB/oct. However, the trend term error in the original acceleration signal causes a drift in the low-frequency band of the spectrum, resulting in a slope of –4.878 dB/oct, much smaller than the standard 6 dB/oct.

Using this SRS as the baseline for shock resistance testing would severely compromise the evaluation of low-frequency shipborne equipment. Therefore, removing the trend term error from the acceleration signal is essential before calculating the SRS to accurately reflect the actual shock environment.

Fig. 3Original acceleration signal and its integrated velocity and integrated displacement

a) Acceleration signal

b) Integral velocity

c) Integral displacement

Fig. 4SRS calculated from original acceleration signal

3.3. Correcting SRS based on improved CEEMD

3.3.1. Acceleration signal decomposition

The CEEMD algorithm was employed to decompose the acceleration signal obtained from the dual-wave shock tester device with 10 as the maximum number of iterations, a white noise coefficient of 0.2, and the number of added noises set to 100 [36].

Fig. 5 shows all IMF components and the residual $\delta$ obtained from CEEMD decomposition of the acceleration signal under Condition I. The time scales within each IMF component are consistent, which effectively resolves the mode-mixing problem encountered during EMD decomposition. Thus, the CEEMD algorithm allows for a more accurate decomposition of the shock acceleration signal into its effective IMF components, representing the actual acceleration and ensuring that error components are associated with the trend term.

Fig. 5CEEMD decomposition results of acceleration signal for Con. I

3.3.2. Effective IMF component selection

3.3.2.1. Selection criterion for cross-correlation analysis

Table 1 lists the cross-correlation coefficients of each IMF component under all four conditions. IMF components with cross-correlation coefficients exceeding the threshold were retained as effective components, while those below the threshold were identified as trend term error components and removed. Fig. 6 illustrates the variations in cross-correlation coefficients for each IMF component under Condition I. The first five IMF components have cross-correlation coefficients exceeding the threshold of 0.1, qualifying them as effective IMF components.

Based on the information shown Table 1 and Fig. 6, under Condition I, the first five IMF components are effective and primarily represent the actual acceleration signal. Starting from the 6^th IMF component, the cross-correlation coefficients all fall below 0.1 and gradually approach zero, indicating that they are dominated by trend term error.

However, cross-correlation analysis is sensitive to noise. When the signal contains significant noise, some effective IMFs may exhibit high correlations with the noise components, resulting in their misclassification as trend components. In such cases, relying solely on cross-correlation analysis to select effective IMFs may lead to unreliable results. Therefore, it is essential to evaluate each IMF more comprehensively by considering both its spectral and temporal characteristics.

Fig. 6Cross-correlation coefficients of each IMF for Con. I

Table 1Cross-correlation coefficients of each IMF

	Coefficient
Condition	IMF1	IMF2	IMF3	IMF4	IMF5	IMF6	IMF7	IMF8	IMF9
Con. I	0.4079	0.4817	0.5915	0.4087	0.1973	0.0864	0.0562	0.0091	–0.0002
Con. II	0.3930	0.4458	0.7158	0.3343	0.2046	0.0658	0.0579	0.0175	0.0002
Con. III	0.3650	0.4883	0.5448	0.3035	0.1407	0.0624	0.0529	0.0013	–0.0012
Con. IV	0.3817	0.4734	0.7440	0.3330	0.2176	0.0891	0.0582	0.0095	0.0012

3.3.2.2. Selection criterion for K-means clustering

The acceleration signals of the dual-wave shock tester are complex oscillatory signals. The K-means clustering method was employed to prevent over-removal of the effective components selected through cross-correlation coefficients, ensuring reliable results.

Table 2 lists the peak moments of each IMF component under the four conditions. The peak moments were divided into effective components and trend components based on K-means clustering, with the number of centroids $K=$2. Fig. 7 illustrates the results after applying K-means clustering to divide the peak moments of IMF components. The first five IMF components fall into the first category; the 6th to 9th IMF components comprise the second category. The peak moments of the first five IMF components are nearly identical, concentrated around 0.08 s.

Table 2Peak moments of each IMF

	Peak moments
Conditions	IMF1	IMF2	IMF3	IMF4	IMF5	IMF6	IMF7	IMF8	IMF9
Con. I	0.080	0.080	0.080	0.078	0.075	0.2175	0.331	1.393	0.544
Con. II	0.081	0.080	0.080	0.079	0.077	0.2653	0.320	1.352	2.065
Con. III	0.081	0.080	0.080	0.076	0.075	0.3114	0.383	0.599	1.106
Con. IV	0.081	0.080	0.080	0.079	0.069	0.2366	0.360	0.7083	1.593

As shown in Table 2 and Fig. 7, the first five IMF components were identified as the primary constituents of the actual acceleration signal. The peak moments show significant deviations from the 6th IMF component onward, indicating that these components are dominated by trend term errors. By combining the results of cross-correlation analysis and K-means clustering, the 6th to 9th IMF components and the residual can be identified as trend term errors caused by the measurement system and external environmental factors. These components should be excluded during the reconstruction of the acceleration signals to ensure accuracy.

Fig. 7K-means clustering division results at the peak moments of each IMF for Con. I

4.1. Correction effects of different decomposition algorithms

EMD is well-suited for analyzing and processing non-stationary and nonlinear shock signals. However, the mode-mixing effect persists during EMD decomposition, leading to inaccuracies in each IMF component’s frequency distribution. EEMD reduces mode mixing by adding white noise to the original signal before decomposition. However, residual white noise remains after averaging in EEMD and cannot be entirely offset. CEEMD addresses this limitation through further refinement.

As shown in Fig. 10, under Condition I, EMD, EEMD, and CEEMD were applied to decompose and reconstruct an original acceleration signal to obtain the SRS. The maximum iteration number was set to 10 for both EMD and EEMD, which is consistent with the parameters of the CEEMD algorithm. Cross-correlation analysis and K-means clustering were also employed to select effective IMF components.

Fig. 8Corrected acceleration signal and its integrated velocity and integrated displacement

a) Corrected acceleration signal

b) Integral velocity

c) Integral displacement

Fig. 9SRS before and after correction

Using EMD for decomposition and reconstruction led to an SRS with mode-mixing issues, causing overlap between IMF frequency bands and distortion in the low-frequency region. As the figure also shows, the resulting slope is 2.037 dB/oct and the low-frequency SRS values deviate significantly from the actual values. When EEMD is applied, the decomposition introduces residual white noise with random characteristics that create noticeable fluctuations in the low-frequency SRS.

CEEMD, however, not only mitigates mode mixing but also effectively eliminates the influence of residual noise. Combined with cross-correlation analysis and K-means clustering, it ensures precise decomposition and selection of effective IMF components. The corrected SRS closely aligns with the low-frequency equal displacement characteristics, achieving a slope of 5.731 dB/oct. In practice, this would significantly enhance the accuracy of shock resistance evaluation tests.

The correction results indicate that a method combining CEEMD decomposition, cross-correlation analysis, and K-means clustering effectively removes trend term errors in acceleration signals. This approach accurately reflects actual shock environments and could play a crucial role in improving the precision of onboard equipment shock evaluations.

Fig. 10Comparison of correction effect of different methods for Con. I

4.2. Low-frequency spring oscillators verification tests

Based on the fundamental principles of the SRS outlined in Section 2.1, an SDOF linear sprin-mass system (referred to as “spring oscillators”) was employed during a shock test to measure the maximum displacement responses at three fixed frequencies. These results were compared with the low-frequency SRS obtained after correction to further verify the accuracy of the proposed method [37].

4.2.1. Structure of spring oscillators

The spring oscillator set used in this study consists of double helical compression springs with natural frequencies of 2 Hz, 10 Hz, and 20 Hz. Each spring oscillator exists independently without mutual interference. The three spring oscillators are mounted on the same base to ensure exposure to identical shock loads during testing. The assembly diagram is shown in Fig. 11.

Fig. 11Assembly of low-frequency spring oscillator

The spring oscillator primarily consists of springs, mass blocks, guiding devices, a shell, and connectors. The springs serve as the core component, providing elastic restoring force. Two pre-compressed springs were used due to the non-linear characteristics of compression springs under tension. The guiding devices ensure that the mass block vibrates in a specified direction, preventing unnecessary displacement or rotation while offering anti-friction and dust-proof functions. Spring seats integrated into the mass blocks and shell position the springs securely.

A magnetostrictive non-contact displacement sensor was employed to record the displacement of mass blocks in the low-frequency oscillator. The sensor has a sensitivity of 100 mV/mm and a range of 200 mm. Compared to contact-type sensors, it reduces damping in the low-frequency oscillator and improves measurement accuracy. Calibration revealed that the actual natural frequencies of the three low-frequency oscillators were 1.937 Hz, 9.687 Hz, and 19.37 Hz, with deviations from the design frequencies of 3.15 %, 3.13 %, and 3.15 %, respectively. The damping ratios were only 0.0371, 0.0145, and 0.0111, while the displacement measurement accuracies were 98.89 %, 96.72 %, and 96.48 %. These results indicate that the low-frequency oscillator set is highly precise and has low damping ratios, making it suitable for accurate measurement in low-frequency shock environments.

4.2.2. Experimental validation results

During the experiment, the low-frequency spring oscillator set was rigidly connected to the shock table surface with bolts. The set was mounted adjacent to the acceleration sensor to ensure consistency with the sensor’s measurements, as shown in Fig. 12. Table 3 presents the maximum displacement responses of the oscillator set under Conditions I-IV.

Fig. 12Shock test site layout

Table 3Maximum displacement response of low-frequency spring oscillator

Conditions	2 Hz	10 Hz	20 Hz
Con. I	51.42 mm	49.60 mm	50.38 mm
Con. II	58.56 mm	61.42 mm	57.78 mm
Con. III	30.94 mm	28.74 mm	30.55 mm
Con. IV	57.66 mm	60.13 mm	58.62 mm

Again, taking Condition I as an example, the low-frequency spectral lines of the corrected SRS were compared to the maximum displacement responses of the spring oscillator at three frequencies. This visual comparison highlights the effectiveness of the correction applied to the low-frequency SRS of the dual-wave shock tester.

Fig. 13 shows the measured displacement responses of the low-frequency spring oscillator set compared to the corrected low-frequency SRS spectral line. The errors at the three frequencies are 3.09%, 5.06%, and 3.94%. The slope of the corrected low-frequency spectral line is 5.6316 dB/oct, which falls within a reasonable range. These results demonstrate that the proposed method effectively eliminates trend term errors in acceleration signals and accurately reflects the actual shock environment.

Fig. 13Measured displacement of spring oscillator versus position of SRS

Fig. 14Correction results of SRS for Con. II to Con. IV

a) Condition II

b) Condition III

c) Condition IV

Table 4Correction errors for four conditions

Conditions	2 Hz	10 Hz	20 Hz
Con. I	3.09 %	5.06 %	3.94 %
Con. II	5.03 %	4.72 %	2.31 %
Con. III	8.45 %	6.34 %	6.17 %
Con. IV	5.82 %	4.36 %	6.24 %

Corrected SRS curves under Conditions II-IV were obtained using the proposed method as shown in Fig. 14. The slopes of the low-frequency spectral lines improved from –6.041 dB/oct, –12.912 dB/oct, and –12.873 dB/oct to 5.535 dB/oct, 5.672 dB/oct, and 5.604 dB/oct after the correction. Table 4 presents the correction results for the low-frequency SRS under all four conditions, comparing them with the measured displacement responses of the low-frequency spring oscillator at the same three frequencies. The maximum errors are 5.06 %, 5.03 %, 8.45 %, and 6.24 %.

These findings indicate that the improved CEEMD method proposed in this paper effectively eliminates trend term errors in acceleration signals, reflects the actual shock environment accurately, and significantly enhances the effectiveness of shock resistance evaluations for onboard equipment.