Weak signal enhancement for small drill condition monitoring in PCB drilling process by using adaptive multistable stochastic resonance

The online tool condition monitoring is demanded to detect the tool wear and to ensure the hole drilling process of printed circuit boards (PCB) goes on smoothly. However, due to the impact of ambient noise caused by the limited size of the small drill and the laminated material of PCB, the tool wear signal features are too weak to extract. The stochastic resonance (SR) method has been proven to be effective in enhancing weak signals among various weak signal extractions. In this paper, an adaptive multistable stochastic resonance is presented to improve the performance of the SR method and process the tool wear signals for PCB drilling. The differential evolution (DE) algorithm is applied to adaptively optimize potential parameters and compensation factor, which makes the SR method suitable for high-frequency signals. Moreover, tool wear experiments with different drill wear are carried out to verify the effectiveness of the proposed method. The results indicate that the proposed method improves the signal-to-noise ratio and has great potential in enhancing weak signals for small drill condition monitoring in the PCB drilling process.


Introduction
The hole drilling process is one of the key processes for PCB production [1]. And the small drills (Φ0.1 ~ 0.3 mm) with high rotating speed (≥ 10,000 rpm) are used for the hole drilling on the PCB substrates of fiber-reinforced laminated composite materials. The combination of composite material, small drill, and high rotating and feeding speed easily leads to tool wear, which gradually deteriorates the drilling quality [2]. The online tool wear monitoring has particular significance for the hole drilling quality of PCB before severe breakage of tool occurs [3]. However, the extraction of weak wear signals in the PCB hole drilling process is very difficult, with the high ambient noise.
In the field of tool condition monitoring, vibration signal analysis has been proven as an effective and reliable method [4][5][6]. The features of the PCB hole drilling process determine the weakness of the collected tool wear signals, which are submerged in loud ambient noise, thereby resulting in a low signal-to-noise ratio (SNR). Hence, the weak signal processing and enhancing method plays a vital role in extracting tool wear features, which influences the accuracy of tool condition monitoring. Many weak signal processing and enhancing methods have been investigated and applied to machinery fault diagnosis, such as singular value decomposition [7], empirical mode decomposition [8], adaptive filtering [9], spectral kurtosis [10], and wavelet transform [11], most of which extract the weak signal features by noise suppression [12]. Nevertheless, noise suppression will decrease the target signal features involved in the noise bandwidth. Different from the above methods, stochastic resonance (SR) is a weak signal enhancing method that utilizes noise energy to enhance the useful signal, which has great potential in detecting bearing and gear fault signals [13]. For instance, Wang et al. proposed a noise-assisted signal feature enhancement and stochastic resonance for early fault diagnosis of rolling bearing [14]. Lu et al. developed an adaptive stochastic resonance method to enhance sound and vibration signals for bearing fault detection [15]. Lei et al. proposed an adaptive SR method for planetary gearbox fault diagnosis [16]. These studies promoted the application of the SR method in rotating machine fault detection. However, the SR method has not been used to process the tool failure weak signal.
In the SR model, the potential function is the primary factor that affects the detection results. In addition to classical bistable models, there are lots of improved and novel potential functions of the SR model reported, such as monostable SR, improved bistable SR, and multistable SR. Zhang et al. proposed an SR method based on a joint Woods-Saxon potential and Gaussian potential and found that the novel model had a smooth bottom and a steep barrier to achieve a better output [17]. Qiao et al. established an unsaturated potential comprised of piecewise bistable potentials to overcome the disadvantage of inherent output saturation of the classical bistable model, which produced higher output SNR [18]. Li et al. investigated the tristable potential with three stable states so as to detect weak signals with lower input SNR and improve the output SNR [19]. Besides the potential function, SR model parameters significantly affect the output SNR. Although manually adjusting the parameters is a conventional and frequently used way, it will be easily trapped in local optimization and thereby limit the detection effect. Many researchers employed intelligent optimization algorithms, e.g., genetic algorithms (GA) [20], grey wolf optimizer algorithm (GWOA) [21], and particle swarm optimization algorithms (PSOA) [22], to optimize the SR model parameters. Therefore, it is of great significance to establish the potential function and optimize the SR model parameters in order to enhance weak signals for small drill condition monitoring.
This research proposes an adaptive multistable SR method to enhance weak signals for small drill condition monitoring based on differential evolution algorithm and multistable SR. Firstly, the structural characteristics and performance advantages of multistable SR are analyzed and compared with that of the classical bistable SR. Then, parameter compensation is performed to improve the multistable SR model that can be applied to high-frequency signal detection. Taking the output mean signal-to-noise ratio (MSNR) as the objective function, the model parameters of multistable SR including potential and compensation parameters are adaptively optimized through the differential evolution algorithm. The experimental results demonstrate that the proposed method is effective to enhance weak signals for small drill condition monitoring.

Theory of SR
Stochastic resonance (SR) describes a physical phenomenon that the motion of particles in the potential well is driven by the interaction of periodic force, noise, and the nonlinear system potential model. Different from other weak signal enhancement methods, SR utilizes the inherent noise in the system to enhance weak fault features and improve the output SNR of the system. A nonlinear system commonly used in the SR model is an overdamped system, as shown in Fig. 1, which can be described by the following Langevin equation: where U(x) is the potential function. s(t) = Acos(2πft + φ) is the input signal with the amplitude A, the frequency f, and the phase φ. η(t) is the system noise and is usually idealized as where ξ(t) is a zero-mean Gaussian white noise and D stands for noise intensity.
The principle of SR is to realize the detection of the weak signal through the noise transfer mechanism under the action of the nonlinear system. The potential U(x) of the nonlinear system plays an important role in the SR method. In the study of stochastic resonance, the most commonly used model is the classic bistable SR model and the classical bistable potential is as follows: where a and b are parameters of classical bistable potential.
According to Eq. (3), the classic bistable SR has two potential wells, as shown in Fig. 2(a). The two potential wells are located at ±x m = ± √ a∕b , respectively, and the barrier height is ∆U = a 2 /4b. Under the combined action of periodic signals and noise, the motion of the particle becomes a periodic transition between bistable potential wells. However, the energy is not sufficient for the particle to jump into the other potential well in the case of small input periodic signal amplitude. With the assistance of the noise, the output of the system will switch between potential wells according to the modulation frequency of the input signal. Therefore, when the input signal, noise, and SR system reach a good synergy, the weak periodic signal component can be enhanced so that it can be identified and detected. Nevertheless, the weak signal processing effect is limited when detecting lower SNR signals because of the classical bistable

Multistable SR
The noise and the useful periodic signal of a specific signal acquired by sensors are restricted in industrial applications. Hence, the potential model of the nonlinear system will greatly influence the output SNR of the SR system. Many studies on different SR potential models such as monostable potential and improved bistable potential have been conducted in order to obtain better weak signal detection results, but few studies on multistable SR models applied to industrial and mechanical fault diagnosis have been carried out. Multistable SR potential has more potential wells and barriers to further improve the noise utilization and the output SNR of the SR system, which means multistable SR can detect weak signals with lower SNR compared to monostable and bistable SR. The multistable SR potential function is given as follows: where a, b, and c are parameters of multistable potential. The shapes of the multistable potential with different parameters are shown in Fig. 3(a). The potential function can be converted to a bistable or monostable SR model by adjusting the parameters. According to Fig. 3(b), the typical multistable potential has three wells located at x min and two barriers located at x max , and the barrier heights on both sides and the middle position are ∆U S and ∆U M ; the functions of which are given as follows: According to Eqs. (5) and (6), the parameters a, b, c of multistable SR potential affect the shape of the potential function and energy barrier for particle transition and consequently have great influences on the performance of the SR method. Better weak signal enhancement results can be obtained through the optimization of potential parameters.

Numerical simulation
In many studies, commonly used SR measurement indicators include signal-to-noise ratio (SNR), power spectral density, cross-correlation coefficient, kurtosis, approximate entropy, and residence time distribution. In this paper, we use SNR, which is the most intuitive indicator, to measure the output performance of the system. A larger SNR indicates the better weak signal detection performance. The specific formula of single frequency signal SNR is as follows: where S(f), N(f), and P are the power spectrum intensity, average noise power intensity, and total power intensity of the output signal with target signal frequency f. The calculation equations of S(f) and N(f) are as follows: where k 0 is the sequence number of the spectral peak at the signal frequency f 0 , so there is f 0 = [k 0 /(N-1)]f s . In the field of mechanical fault diagnosis, f s and N are sampling frequency and length of the input signal, respectively. x(n), n = 1, 2, …, N is the discrete form of the output signal x(t), and X(k) is frequency-domain value after Fourier transform of x(n). In addition, the actual signals collected by sensors are discrete signals, so the SR model in Eq. (1) is discretized to calculate the numerical solutions by using the fourth-order Runge-Kutta equation as follows: where y(n), n = 1, 2, …, N is the practical discrete form of the input signal y(t) sampled by sensors including useful signal s(t) and noise η(t), and h is the calculation step. According to Eq. (9), the output of the multistable SR model can be expressed as In one case of the numerical simulations, the parameters of multistable SR in Eq. (4) are a = 3, b = 3, c = 30, and the useful signal is s(t) = 0.5cos (0.02πt) with the amplitude A = 0.5 and frequency f = 0.01 Hz, which means the periodic signal amplitude is smaller than the side barrier height of the potential. The sampling frequency is 5 Hz, and the length of the data is 4096. The noise density of added Gaussian white noise is 6. The time-domain waveform of noisy signal and power spectrum are shown in Fig. 4(a), with the input SNR calculated as − 7.34 dB. The spectrum peak at 0.01 Hz is not obvious in the spectrum. The multistable SR output and power spectrum are shown in Fig. 4(b). It can be seen that the weak periodic signal is enhanced with the output SNR of 7.52 dB compared to that in Fig. 4(a). Figure 5 shows the trends of the SNR of a multistable SR method with different potential parameters of a, b, and c. The trend of the SNR with variable parameters of a = 1 ~ 5 and the fixed parameters of b = 3 and c = 30 is shown in Fig. 5(a). The trend of the SNR with variable parameters of b = 2 ~ 5 and the fixed parameters of a = 3 and c = 30 is shown in Fig. 5(b). The trend of the SNR with variable parameters of c = 10 ~ 50 and the fixed parameters of a = 3 and b = 3 is shown in Fig. 5(c). The results show that SR performance can be induced by adjusting any parameter in the multistable SR system. At the same time, it can be found that under the coupled action of different parameters and noise, there may be more than one parameter interval for SR to achieve effective detection of weak signals. Hence, it is of great significance to optimize the potential parameters simultaneously in order to obtain the best SR output performance. Figure 6 shows the trend of the SNR of a multistable SR method with the simply optimal parameters of a = 3, b = 5, c = 30 and a classical bistable SR method with a = 5 and b = 1.5 with different noise intensities. The output SNR of the multistable SR method is greater than that of the classical bistable SR method with lower noise intensity. In addition, the output of the classical bistable SR method diverges when the noise intensity is larger than 10, which means the classical bistable SR method cannot process the signals with lower SNR. The numerical simulation results demonstrate that the multistable model can transfer more noise energy into useful signals, and the enhancement effect of the noise on the periodic signal in the multistable SR method is larger than that in the classical bistable SR method for the same input signal. Moreover, the multistable SR method is able to process and enhance noisy signals with lower SNR. Therefore, the multistable SR model has more advantages in improving the output SNR and achieves better weak signal enhancement performance.  Figure 7 shows the time-domain waveform, the spectrum of the noisy signal, and the output of multistable SR with the frequency of the input signal increased to f = 0.3 Hz and other parameters unchanged. It can be seen that the spectrum peak at 0.3 Hz of the output signal processed by SR becomes lower, and the weak signal enhancement is not realized. The simulation results show that the response frequency band of the SR model is mainly in the low-frequency region, which means the SR model shows low-frequency characteristics and is suitable for useful signals with small parameters. To solve the above problems, the parameter compensation method is adopted to detect the useful signals with high frequency. The system model is as follows:

Parameter compensation for large-parameter system
Both sides of Eq. (11) are integrated to get where K is the compensation parameter. According to Eq. (12), when the value of K is equal to 2πf i , the amplitude A i of the SR output signal remains the same. The value of K significantly affects the output performance and should be larger than 2πf i because the amplitude A i of the SR output signal is enlarged. Hence, the value of the compensation parameter K can be optimized to improve the output SNR of the multistable SR system.

Specific process of the method
The actual signals are characteristically multi-frequency, and the mean signal-to-noise ratio (MSNR) is used to judge the SR output performance. Hence, Eq. (7) and (8) become as follows: where n is the number of characteristic frequencies and f i is the characteristic frequency. It is known from the above analysis and discussion that factors (potential parameters and the compensation parameter) of multistable SR have different influencing mechanisms on SR output performance. The conventional way is to adjust these parameters manually in order to achieve the optimal SR output results. However, the change of one factor needs the sequential adaptation of other factors to maintain the output of the whole system in the desired range or even the optimal solution, which is a complicated process. To improve the efficiency of parameters adjustment, an adaptive optimization method based on differential evolution (DE) was developed to realize parameters' adaptive adjustment, which is one kind of genetic algorithm and has high computational efficiency. It uses a unique mutation mechanism to avoid falling into local optimization. The specific algorithm flow is as follows.
Step 1. Initialization: Each individual corresponds to a solution of the problem , where x i,1 represents the i-th individual in the population, and the initial population (population size is NP) is randomly generated within the given constraint boundary, namely where u k and l k are searching upper and lower bounds of the k-th dimension variable, d is the dimensionality of the problem, and rand() represents random number uniformly distributed between 0 and 1.
Step 2. Mutation: In each generation search, a target individual t i (g) is generated as each individual x i (g) in the current population through mutation operation in the population, where g represents the evolutionary algebra, and the following mutation mechanism is adopted: where r 1 , r 2 , r 3 , r 4 ∈ {1, 2, ⋯ , N} are different integers that are not equal to i, x best (g) is the best individual in the g-th generation population, and F is a scaling factor between 0 and 1.
Step 3. Cross: The parts of the variables of the current individual x i are replaced by the corresponding variables of the target individual t i to generate the test individual v i . The two-phase crossover method is adopted: where rnd is an integer uniformly distributed between 1 and d.
Step 4. Selection: Using the greedy selection method, for the test individual and the current individual, the better individual is selected to enter the next-generation search, namely The greedy selection mechanism ensures that the individuals in the next-generation population are at least no worse than those in the current population so that the average performance of the population is improved and the optimal solution is gradually reached.
The procedure of the adaptive multistable SR method is shown in Fig. 8 and described as follows.
(1) Signal preprocessing. Since mechanical sensor signals such as vibration signals are generally modulated, they are firstly demodulated by Hilbert Transform [22] to obtain their envelope signals. The characteristic frequencies are calculated based on the practical processing state. The large-parameter signals are preprocessed using the parameter compensation method. (2) Parameter initialization. The preprocessed signal is put into the adaptive multistable SR model, and the model parameters including SR potential parameters and the compensation parameter are optimized by the DE method according to Eq. (14)- (17). (3) Output calculation. The output signal is achieved by solving the Eq. (1) using the above-mentioned discrete Fourthorder Runge-Kutta method. Then, the power spectrum of the output signal is obtained through Fourier transform, and the MSNR is calculated in terms of Eq. (13). (4) Output assessment. The optimal matching parameters are searched corresponding to the maximal output MSNR. The optimal output signal of the multistable SR model is achieved through the substitution of the optimal matching parameters.
x i (g), others

Experimental setup
To verify the effect of the proposed method in small drill condition monitoring, a small drill condition monitoring experimental setup is established, as shown in Fig. 9, which mainly consists of the motion platform, the high-speed electric spindle, the vibration signal acquisition system, and the open CNC system. The maximum speed of the Z-axis is 30 mm/s, and the maximum rotating speed of the spindle is 80krpm, which meets the requirement of PCB small holes drilling. The CNC system is composed of an industrial personal computer and a programmable multi-axis controller (PMAC). The vibration signal acquisition system consists of two three-axis acceleration sensors (356B11, PCB, USA), vibration signal acquisition device (PXIe-1082, PXIe-6124, PXI-4498, NI, USA), and signals acquisition software written based on Labview. The signal frequency range of the sensor is 2 ~ 10,000 Hz, and the maximum sampling frequency of the data acquisition card is 204.8 kS/s. Three-axis acceleration sensors are installed on the spindle and workpiece. The workpiece material is a multi-layer printed circuit board (S1000-T1.60 mm, SYTECH, China), and the drilling tool is a small drill with two cutting edges and a diameter of 0.8 mm (A129UCQM0.80, Jinzhou, China).

Weak signal enhancement of tool wear features
In the PCB drilling process, the rotating speed of the spindle is n = 65 krpm and the number of the cutting edge is z = 2, so the spindle rotation frequency and the tooth passing frequency are calculated as follows: The vibration signals of the drilling process are continuously collected with the sampling frequency of 50 kHz by the signal acquisition system and divided into signals of drilling Fig. 8 The proposed adaptive multistable SR procedure to enhance weak signals for small drill condition monitoring one hole. With the number of drilling holes increasing, the drilling tool flank wear area with the triangle shape observed by digital microscope (VHX-600, Keyence Corp., Japan) is increased. The pictures of drilling tool flank wear with different degrees are shown in Fig. 10(a). It can be seen that the new tool has two sharp blades, and the flank wear area of the drilling tool is the triangle distribution from the outer edge point to the drill core and is increased with the drilled holes number increasing, which is caused by the severe friction between the flank and printed circuit board. The increasing flank wear influences the dynamic characteristics of the drilling tool. The midterm wear tool with drilling holes of 2500 has a half-worn flank face, and the severe wear tool with drilling holes of 4000 has a completely worn flank face. The spectrums of the vibration data from the sensor mounted on the spindle with the new tool, midterm wear tool, and severe wear tool are shown in Fig. 10(b), respectively. It is clear that the spectrum has two peak values at frequencies of 1090 Hz and 2180 Hz, which are approximately consistent with the characteristic frequencies including spindle rotation frequency f = 1083.3 Hz and tooth passing frequency f z = 2166.6 Hz, for there is a small error between the actual speed and the set speed of the electric spindle. However, the calculated characteristic frequencies are in the same frequency band as the ambient noise, which leads to a low signal-to-noise ratio. Because of the strong ambient noise, the spectrums show that other interference frequencies are in the range of 5 ~ 10 kHz and part of their spectrum values are even larger than characteristic frequencies, which interferes with the characteristic frequencies and may have negative effects on the extraction of drilling tool wear features.
In order to enhance the characteristic weak signals of the drilling process and further extract the tool wear features, the proposed adaptive multistable SR method is adopted to process the vibration signals acquired in the experiments. Due to the coupled effects of multistable SR factors (potential parameters and the compensation parameter) on the output performance, the ranges of parameters are achieved by manual adjustment. Hence, the proposed adaptive multistable SR model can be considered as a single objective optimization, which is expressed as follows: (19)  where K is the compensation parameter, a, b, c are SR potential parameters, and R represents the value of calculation step h. To solve the Eq. (19), the iteration process of MSNR is illustrated in Fig. 11. The MSNR of DE is about − 76.7 dB at the initial stage. As the number of iterations increases, the value of MSNR is increased and reaches the maximum of 46.45 dB at the 87th generation. The result indicates that the proposed method is effective to improve output performance.
The time-domain waveforms, the spectrum of the raw signal, and the adaptive multistable SR model with midterm tool wear are shown in Fig. 12. According to the results of Eq. (19), the optimal corresponding parameters are K = 9987, a = 0.1676, b = 0.069, c = 25.4953, R = 1. As shown in Fig. 12(a), the raw signal waveform of drilling one hole has two interference peaks at the beginning and end of the drilling because the printed circuit board is covered by an Aluminum sheet to avoid inlet burrs in the practical drilling process. In addition, interference frequencies cannot be ignored in the spectrum. As shown in Fig. 12(b), the time-domain waveform of the adaptive multistable SR model has no extra interference peaks, and the spectrum shows that most of the interference frequencies are submerged and characteristic frequencies become obvious. The values of the spectrum at characteristic frequencies of 1090 Hz and 2180 Hz increase from 0.039 and 0.015 to 0.054 and 0.033, respectively. The calculated MSNR of the output signal is increased to 46.4536 dB from 25.89 dB of the raw signal.
The spectrum values of three drill wear states at spindle rotation frequency are shown in Fig. 13. Because of the increased friction force due to flank wear of the drilling tool, the spectral energy at the spindle rotation frequency during the drilling process shows an increasing trend. As shown in Fig. 13, the spectrum value of the midterm wear tool raw signal at 1090 Hz is smaller than that of the new tool and severe wear tool raw signal. This means the spectrum value of the raw signal is not linearly related to the tool wear value. However, the spectrum values of output signals are increased with the tool wear area increasing. It is consistent with vibration signal laws with different drill wear states because axial force and vibration acceleration are increased with the tool wear area increasing. In addition, the spectrum values of raw signals are suppressed by the ambient noise. The weak features are enhanced while improving the SNR by adopting the proposed method. The results demonstrate that the spectrum value at the spindle rotation frequency processed by the proposed method can be selected as the tool wear feature for small drill condition monitoring. According to the results, the spectrum values of output signals at tooth passing frequency f z = 2180 Hz are not linearly related to the tool wear value compared with those at the spindle rotation frequency, which demonstrates that tool wear of the two blades and flank is uniform and the spectral energy difference of different drill wear states mainly reflects in the spindle rotation frequency. In conclusion, the proposed

Conclusions
In this paper, in order to enhance weak signals for small drill condition monitoring, an adaptive multistable SR model is proposed, which is constructed by using multistable SR, parameter compensation method, and differential evolution algorithm. The main conclusions can be drawn as follows.
(1) The SR characteristics of classical bistable SR and multistable SR indicate that the SR model has an important influence on the energy distribution of the system's output by the periodic force acting on the particle and the parameters of the SR potential significantly affect the output performance. (2) By means of numerical simulations, the structural characteristics of the multistable SR model and the detection capability of the multistable SR method for weak signals are analyzed. It is found that the multistable SR method increases the output SNR and improves the detection effect of SR, which shows a better effect than the classical bistable SR method in weak signal enhancement. (3) The parameter compensation method is adopted to overcome the problem that the SR model has low-frequency characteristics. The proposed adaptive multistable SR model uses the DE algorithm to adaptively adjust the parameters (compensation parameter and SR potential parameter) of the SR model.  The spectrum values of three wear states (5) The output signal spectrum values of three drill wear states at spindle rotation frequency are increased with the tool wear area increasing, which indicates that the weak features are enhanced by adopting the proposed method. The proposed method is verified to be feasible in weak signal enhancement for small drill condition monitoring. Availability of data and material The datasets used or analyzed during the current study are available from the corresponding author or the first author on reasonable request.

Code availability
The codes used during the current study are available from the corresponding author or the first author on reasonable request.

Declarations
Ethics approval All the authors declare that no animals or human participants are involved in this research. There are no ethical issues to declare.

Consent to participate
All the authors declare that no human participants are involved in this research. Only authors participate in the research work of this paper.

Consent for publication
All authors agree to publish the research results in this paper.

Competing interests
The authors declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work. This manuscript has not been published or presented elsewhere in part or in entirety and is not under consideration by another journal. I have read and understood your journal's policies, and I believe that neither the manuscript nor the study violates any of these.