A denoising scheme for DSPI phase based on improved variational mode decomposition

Mechanical Systems and Signal Processing 110 (2018) 28–41

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

A denoising scheme for DSPI phase based on improved variational mode decomposition Qiyang Xiao, Jian Li ⇑, Zhoumo Zeng State Key Laboratory of Precision Measurement Technology and Instrument, Tianjin University, Tianjin 300072, China

a r t i c l e

i n f o

Article history: Received 1 November 2017 Received in revised form 6 March 2018 Accepted 9 March 2018

Keywords: Digital speckle pattern interferometry Phase map Denoised Improved variational mode decomposition

a b s t r a c t Aiming at reducing the noise interference of digital speckle pattern interferometry (DSPI) phase map, a DSPI phase denoised method based on improved variational mode decomposition (IVMD) is proposed in this paper. Firstly, the IVMD method is adopted to decompose DSPI phase map and obtain the optimal mode components; Then, according to the mode components characteristics, an adaptively mode threshold method is proposed to process the mode components. Finally, the denoised mode components are reconstructed to obtain the denoised DSPI phase map. The result indicates that the proposed method can effectively filter out noise interference, and the peak signal-to-noise ratio (PSNR) is higher than other denoised methods. Ó 2018 Elsevier Ltd. All rights reserved.

1. Introduction In recent years, different kinds of composite materials have been widely used in the field of aeronautics and astronautics. However, the components may undergo significantly declined performance and reduced service time due to defects of these composite materials, such as deformation and displacement generated during processing, manufacturing, or serving process. Therefore, these composite materials are required to be tested and evaluated for providing solutions for subsequent urgent repair and guaranteeing robust operation of the aerospace equipment [1,2]. Relevant researches have indicated that when deformation or displacement occurs on the surface of an object exposed to coherent light, the object surface deformation will cause phase changes of the speckle interferogram, which is called digital speckle pattern interferometry (DSPI). The digital speckle pattern interferometry is a nondestructive testing technique which is capable of measuring deformation. The technique can obtain the information on object deformation by measuring the phase changes of the speckle pattern before and after deformation [3,4]. In the DSPI phase maps, the noise interference leads to a low peak signal-to-noise ratio (PSNR) of the speckle interferogram and severely influences the subsequent phase unwrapping result and accuracy [5,6]. Hence, it is necessary to study noise reduction for the DSPI phase maps to filter out noise and improve the measuring accuracy. For the speckle noise of DSPI phase maps, numerous different de-noised approaches have been proposed in the past few years, such as average filtering method [7], the median filtering method [8], Fourier transformation method [9,10], partial differential equation (PDE) filtering method [11], multi-scale filtering method [12], etc. The average filtering method is a linear filtering method where the value of target pixel point on the map is substituted with an average value of several pixels around that target pixel, wherein target pixel is considered as a center. The average filtering method also damages the image details while filtering the noise. As a result, the processed image is blurred and information on abrupt phase change is lost ⇑ Corresponding author. E-mail address: [email protected] (J. Li). https://doi.org/10.1016/j.ymssp.2018.03.014 0888-3270/Ó 2018 Elsevier Ltd. All rights reserved.

Q. Xiao et al. / Mechanical Systems and Signal Processing 110 (2018) 28–41

29

[7]. The median filtering method is a non-linear noise reduction method where the median of various points in the target point domain is substituted with a value of that point. Compared with the average filtering method, the median filtering method can maintain the marginal information while reducing noise in a simple manner. However, this method also has some disadvantages including long computation time and poor smoothing effect [8]. The filtering method based on Fourier transform is to transform the image into the frequency domain, and then achieve the purpose of filtering by suppressing or enhancing the different frequency components, finally transform the image into the space domain to obtain the filtered image [9]. A window function is required to be selected when the Fourier transformation is used for noise reduction [10,13]. The basic idea of the PDE method is to treat image processing as a continuous rather than a discrete process, converting a given image into a form of partial differential equations to obtain a de-noised image by solving partial differential equations [14]. Compared with other methods, partial differential equation filtering method is simple, and the filtering effect is better. However, the results of the PDE method often have a great relationship with the filtering model [11]. Selecting different iterations and different steps in the filtering process will produce different filtering effects. These parameters often through a large number of experiments and repeated comparative analysis to determine [10,13]. The multi-scale filtering method includes wavelet transformation [14], empirical mode decomposition [15], variational mode decomposition [16], etc. The wavelet filtering method is used to sparsely decompose the image at different scales and reconstruct it using the threshold-free noise-free scale to obtain a de-noised image. Wavelet filter method is sparse, multi-resolution and other characteristics; therefore, this method has been widely used. Wavelet transform method need to select the wavelet basis function, in addition, once the wavelet basis function and decomposition scale, wavelet transform results will only be a function of the sampling frequency [17]. The empirical mode decomposition (EMD) is a non-stationary signal processing method which is capable of adaptively decomposing a complicated multi-component signal into several components [18,19]. The empirical mode decomposition method processes the signal according to the intrinsic characteristics of the signal itself without requiring artificial selection of a primary function. Hence, the empirical mode decomposition method has been widely used in actual situations although having some disadvantages such as mode mixing, end effect, and poor noise immunity [19]. The variational mode decomposition (VMD) is a new non-stationary signal processing method. By iteratively searching the optimal solution for the variational mode during decomposition, the method can determine the frequency center and bandwidth of the decomposed component and eventually realize the purpose of decomposing the nonstationary signal [20]. VMD decomposes and converts the signal into a non-recursive variational mode decomposition mode and rationally control the convergence conditions. Therefore, it can effectively eliminate mode mixing. The VMD method has been widely applied in a series of fields such as fault diagnosis, electroencephalogram signal processing, and laser speckle imaging [16]. However, this method requires the presetting of the decomposed mode number, which may lead to information loss or over decomposition [21]. To solve the problem of decomposed mode number, this paper proposes an improved variational mode decomposition (IVMD) method, which can determine appropriate mode number using the orthogonal index based on the noise statistics characteristics. Firstly, the IVMD can be used to decompose the DSPI phase map into several mode components; then, the mode components are analyzed to self-adaptively compute the mode threshold of the mode components. Secondly, the mode components are processed on the basis of the mode threshold. Thirdly, the noise is filtered out to obtain the denoised mode components. Finally, the denoised mode components are reconstructed to obtain the denoised DSPI phase map. The remainder of this paper is organized as follows. Section 2 presents the improved VMD method and the simulation. An adaptive mode threshold method is proposed to deal with mode components in Section 3. In Section 4, the DSPI phase map denoised scheme is explained. In Section 5, the proposed DSPI phase denoised scheme is experimentally validated, and compared with average filtering method and median filtering method. Finally, a brief conclusion is given in Section 6. 2. Improved variational mode decomposition 2.1. Brief introduction of VMD Variational mode decomposition is a new self-adaptive image processing method, which was proposed by Dragomiretskiy in 2015 [16]. This method achieves self-adaptive image decomposition through constructing and solving variational problems. Any image is decomposed into a series of band-limited intrinsic mode functions (BLIMFs) by VMD. Each mode has a limited spectral bandwidth: the first mode component contains the highest frequency while the last one contains the lowest frequency. The mathematical expression of complex decomposed image f ðxÞ is:

f ðxÞ ¼

X uk ðxÞ

ð1Þ

k

where uk ðxÞ denotes the band limited intrinsic mode functions (BLIMFs), f ðxÞ is a 2D signal. The 2D analytic signal of mode uk ðxÞ can be obtained using Hilbert transform method:

 uAS;k ðxÞ ¼ uk ðxÞ  dðhx; xk iÞ þ

 j dðhx; xk ; ?iÞ phx; xk i

ð2Þ

30

Q. Xiao et al. / Mechanical Systems and Signal Processing 110 (2018) 28–41

In Eq. (2), uAS;k ðxÞ is the 2D analytic signal,  denotes convolution, dðxÞ is the dirac function, xk represents the center frequencies set of components after decomposition. The VMD algorithm transfers the signal decomposition process into a variational framework. Hence, the VMD decomposition process is an optimal-solution processing for a constrained variational process. The constrained variational mode is defined as follows:

( ) X    jhxk; xi 2  min ak r uAS;k ðxÞe 2

fuk g;fxk g

ð3Þ

k

where ak is the penalty parameter. For constrained-variation problems, the augmented Lagrange function is introduced to transform the constrained variation problem into an unconstrained variation problem [20]. The mathematical expression of augmented Lagrange function is as follows:

2 *  +   X   X X    jhxk ;xi 2  Lðfuk g; fxk g; kÞ ¼ ak r uAS;k ðxÞe þ f ðxÞ  uðxÞ þ kðxÞ; f ðxÞ  uk ðxÞ 2   k k k

ð4Þ

2

where k is the Lagrange multiplier. In order to obtain an optimal solution, the alternate direction method of multipliers (ADMM) is used to calculate the saddle point of the augmented Lagrange function, i.e., obtaining the optimal solution of , the constrained-variation equation. The saddle-point problem is solved by the alternate renewal of components unþ1 k

xnþ1 , and Lagrange multiplier knþ1 , i.e., by solving the optimal-solution problem of the variational problem. The mathematk and xnþ1 is as follows: ical expression of unþ1 k k unþ1 k

8 2 9   = <   X 2  kðxÞ   jh x ;xi ¼ arg min ak r uAS;k ðxÞe k 2 þ f ðxÞ  ui ðxÞ þ   : 2  ; uk i

ð5Þ

2

n

  2 xnþ1 ¼ arg min ak r uAS;k ðxÞeihxk ;xi 2 k

o

ð6Þ

xk

2.2. Improved VMD In VMD, the ADMM method is used to seek the optimal solution for the variational constraint mode and decomposes the signal into K mode components according to the optimal solution. However, the decomposed mode number K is preset during the process of VMD decomposition. The mode number can influence the accuracy and effect of decomposition, and thus easily generate false components. A one-dimensional signal collected by the sensor comprises noise signal and useful signal. Wu et al. conducted a statistical analysis for the noise energy density of signal, and obtained the decomposed mode number based on the length of the one-dimensional signal [22,23]. Images are a two-dimensional signal. The processing method for one-dimensional signal is extended to the two-dimensional signal for analyzing its noise energy characteristics. The mode number K is estimated according to the image length and width. The orthogonal index can be used to measure the decomposition degree of the complicated signal. A smaller orthogonal value will lead to a lower level of correlation of the decomposed components and a better decomposition result [24,25]. Thus, the mode number is estimated according to the image size. The most suitable mode number is selected according to the orthogonal value. The specific steps are as follows: The dimension of the DSPI phase map is M  N. The number of components following VMD decomposition is K: (1) When the phase map meets M ¼ N, the number of decomposed mode components K is:

K ¼ log2 N  1

ð7Þ

(2) When the phase map meets M–N, the number of the decomposed mode components K:

if M 6 N; log2 M  1 6 K 6 log2 N  1

ð8Þ

if M P N; log2 N  1 6 K 6 log2 M  1

ð9Þ

The orthogonal index under different K values is computed. The K value corresponding to the minimum orthogonal value is the appropriate number of mode components. The VMD algorithm continuously updates the modes in frequency domain, and finally transforms them to time domain by the Fourier inversion transformation. The updating is executed iteratively until the convergent Eq. (10) is satisfied.

X  .  unþ1 ðxÞ  un ðxÞ2 un 2 6 e k k 2 k 2

ð10Þ

K

where unþ1 ðxÞ denotes the decomposed components, k position are shown in Fig. 1.

e is epsilon. The specific steps for improved variational mode decom-

Q. Xiao et al. / Mechanical Systems and Signal Processing 110 (2018) 28–41

31

Input phase map Estimate the mode number K accroding to map length Determin mode number K using the orthogonal value Constructing variational problems ADMM method was used to solve the optimal solution

N

Judge the termination conditions

Y Obtain a series of mode components Fig. 1. Flow chart for the improved VMD.

2.3. Simulation experiment To verify the effectiveness of improved VMD method mentioned above, simulation experiment was carried out. The original phase measured by the digital speckle pattern interferometry is the wrapping phase of the true phase with respect to 2p. Its mathematical expression is as follows:

u2 ¼ modðu1 ; 2pÞ

ð11Þ

u1 is the true phase. For the sake of simplicity, it was assumed that u1 varied along the width direction of the image and remained unchanged along the height direction of the map. The size of the phase map is 512  212. According to

Fig. 2. Simulation phase map.

32

Q. Xiao et al. / Mechanical Systems and Signal Processing 110 (2018) 28–41

Eq. (10), we can obtain a simulation wrapping phase map, as shown in Fig. 2.Based on the size of the phase map, the decomposed mode number was estimated in advance. As mentioned above, the computed number of mode components was 6 6 k 6 8. The orthogonal index could be used to measure the decomposition degree of the simulation phase map. A smaller orthogonal value would lead to a lower level of correlation of the decomposed components and a better decomposition result. Thus, the orthogonal index was employed to select the most suitable mode number decomposed by IVMD in this paper. The orthogonal values under different K values are shown in Fig. 3. As seen from Fig. 3, a lower level of correlation among different mode components leads to a lower degree of mode mixing, more thorough VMD algorithm decomposition, and a smaller orthogonal value. When K = 8, the decomposed orthogonal value is the minimum. Thus, K = 8 is selected as the appropriate decomposed mode number for the simulation phase map. The energy storage level is the contrast ratio of the energy before and after phase map decomposition, which can be used to measure the orthogonal index among different components after phase map decomposition [26]. The mathematical expression is as follows:

ESL ¼

M X N Kþ1 X X

!

u2i ðx; yÞ

ð12Þ

2

x¼1 y¼1

i¼1

½f ðx; yÞ  rðx; yÞ

f ðx; yÞ is the complex image, ui ðx; yÞ is the decomposed component, rðx; yÞ is a residual. A better orthogonal index would lead to a lower level of correlation among different components and less energy leakage during the process of simulation phase map decomposition. The energy storage levels of simulation phase map after IVMD decomposition were computed when K = 5, 6, 7, 8, 9, and 10. Table 1 shows the energy storage levels of the simulation phase map after IVMD decomposition under different K values. The energy storage level reaches the highest value when K = 8, indicating that the energy leakage is the least after decomposition, and the level of correlation among different mode components at this point is the lowest. The improved VMD algorithm can effectively compute the decomposed mode number K value, accurately decompose the simulation maps, and obtain the optimal mode components. 3. Denoising scheme thresholds 3.1. Adaptive calculation mode thresholds The original DSPI phase map was decomposed to obtain multiple mode components by IVMD. Some mode components were closely associated with the phase information, while some mode components, including a great deal of noise, were less closely correlated with the phase map. Thus, it was necessary to reduce noise of the mode components and improve the peak signal-to-noise ratio (PSNR) of the DSPI phase map. The traditional wavelet threshold denoising method is a time domain method with excellent denoising capacity. Using appropriate threshold, this method processes the wavelet decomposition coefficient, and then reconstructs the processed wavelet decomposition coefficient to obtain the noise reduction signal. The wavelet threshold denoising method can not only decompose the low-frequency part of the image, but also decompose the high-frequency part, so as to effectively extract useful information of various frequency bands [27,28]. The wavelet threshold denoising method has been widely used in a series of fields such as data noise reduction, vibration information extraction, and image processing. According to the principle of the wavelet threshold denoising method, a mode threshold-based de-noising method was used to filter out the noise inherent in mode components. IVMD was used to decompose the phase map, compute the appropriate mode-dependent threshold, analyze the mode components, and obtain the denoising component. The specific steps for noise reduction method based on IVMD and mode-dependent threshold are as follows:

 ui ðxÞ ¼

sgnðui ðxÞÞðjui ðxÞj  T i Þ; jui ðxÞj > T i 0;

jui ðxÞj 6 T i

;

i ¼ 1; 2; 3; . . . ; J

ð13Þ

where ui ðxÞ is the mode components containing noise, ui ðxÞ is the mode components after noise reduction, i is the number of noise mode components, T i is a mode-dependent threshold, which is associated with the mode components ui ðxÞ. The mathematical expression of T i is as follows:

T i ¼ ri

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 lnðM  NÞ

ð14Þ

ri ¼ MADi =0:6745

ð15Þ

MADi ¼ medianðv arðui ðxÞÞÞ

ð16Þ

MADi represents the variance medium value of the ui ðxÞ component. ri represents estimated value of noise of the ui ðxÞ component. The amplitude of the noise component of the actual signal is small. Hence, using mode-dependent threshold T i can effectively achieve noise reduction for various types of signals. This mode-dependent threshold de-noising method typically works well because the amplitude of the noise component is inherently lower in the signal. However, the distribution of

33

Q. Xiao et al. / Mechanical Systems and Signal Processing 110 (2018) 28–41 5 4.5

Orthogonal values(%)

4 3.5 3 2.5 2 1.5 1 0.5 0

5

6

7

8

9

10

Mode number k Fig. 3. Orthogonal values corresponding to different mode numbers.

Table 1 Energy storage levels corresponding to different mode numbers. K

5

6

7

8

9

10

Energy storage levels

0.9331

0.9023

0.9214

0.9765

0.9132

0.9471

Fig. 4. Simulated phase map with noise.

60

BLIMF energies

Threshold values

50

Energies

40

30

20

10

0

1

2

3

4

5

6

7

8

BLIMF components Fig. 5. The decomposed mode threshold and mode absolute value.

34

Q. Xiao et al. / Mechanical Systems and Signal Processing 110 (2018) 28–41

Fig. 6. The DSPI phase noise reduction methods: (a) IVMD, (b) VMD.

15

14

14

13

13

12

PSNR

11

11 10

10

9

9

8

8

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.0

2

3

4

5

6

7

8

9

The value of rho

The value of C

(a)

(b)

14 13 12

PSNR

PSNR

12

11 10 9 8

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

The value of beta

(c) Fig. 7. Peak signal-to-noise ratio under different parameters: (a) parameter C; (b) parameter q; and (c) parameter b.

DSPI phase map

Estimate the mode number

IVMD

Calculate Mode dependent threshold,

Obtain the denoised mode components

Fig. 8. Flow chart for the DSPI phase denoising method.

Reconstruct denoised mode components

10

Q. Xiao et al. / Mechanical Systems and Signal Processing 110 (2018) 28–41

VBS

laser Lens

35

Mirror

Concave lens

Fiber

Camera

IL

Aperture

Measur ing Target

Fig. 9. Schematic diagram for DSPI system.

mode components samples may not be gaussian with variance equal to the noise variance, due to the self-adaptive nature of the IVMD. In practice, the noise is colored because it has different energies in the mode components [29,30]. This results in multiples of the mode-dependent threshold, which are computed as follows.

Ti ¼ C Ek ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ek 2 lnðM  NÞ

E21 k q ; b

ð17Þ

k ¼ 1; 2; 3; . . . ; J

ð18Þ

In Eqs. (17) and (18), E1 is the energy value of the mode component u1 ðxÞ, C is a constant, b; q are the parameter during the process of IVMD decomposition, which mainly depend on the number of screened iterations. 3.2. Simulation noise reduction experiment For comparison, the proposed denoised method is used to process simulate Fig. 2. The random speckle noise with a mean value of 0 pixels and variance of 0.75 pixels is added to Fig. 2. The simulated phase map containing noise is shown in Fig. 4. Thus, the obtained phase map is similar to the phase map measured by DSPI. IVMD was used to decompose simulated phase map into eight decomposed components, according to the principle of noise reduction described above, the absolute value of the mode component and the noise threshold were calculated respectively, and the result is shown in Fig. 5. Fig. 5 shows that the absolute values of the mode components u3, u4 and u5 are greater than the noise threshold. These mode components contain the primary phase information; the other components are mainly noise components. Therefore, these three components are selected to reconstruct to obtain a denoised phase map, as shown in Fig. 6(a). For comparison,

Fig. 10. Tested object.

36

Q. Xiao et al. / Mechanical Systems and Signal Processing 110 (2018) 28–41

Fig. 11. Out-of-plane deformation DSPI phase map.

VMD is used to process Fig. 4 to filter out noise interference and obtain the phase map after noise reduction. The result is shown in Fig. 6(b). As seen from Fig. 6, the VMD-based noise reduction method is prone to false components, which cannot effectively filter out noise interference. The proposed method in this paper can obtain the optimal mode numbers, filter out the noise adaptively and extract the noise-free mode components. The peak signal-to-noise ratio (PSNR) is used to measure the effect of noise reduction. The peak signal-to-noise ratio reflects the retention of image detail. A higher peak signal-to-noise ratio indicated a higher level of detail preservation. For these constant values b, q, and C, different constant values were used to calculate the peak signal to noise ratio after noise reduction, taking the average of multiple experiments, the results are shown in Fig. 7. Fig. 7 shows the peak signal to noise ratio under different parameters, the peak signal to noise ratio changes with different parameter values. When the parameters b, q, and C were 0.719, 2.01 and 0.7 respectively, the peak signal to noise ratio was the largest, which showed that the noise reduction effect was the best. Therefore, the values of parameters b, q, and C are respectively 0.719, 2.01 and 0.7. 4. DSPI phase map denoising scheme The phase map was collected with a CCD camera. The decomposed mode number was estimated according to the length and width of the DSPI phase map. The mode number K corresponding to the minimum orthogonal value was regarded as the most suitable mode number. The DSPI phase map was processed to obtain K mode components using IVMD. The decomposed mode was analyzed, and the mode threshold was adaptively computed. The threshold was used to process the noise of decomposed components. The denoising mode components were reconstructed to obtain the denoised phase map. The specific process is shown in Fig. 8.

1.4

Orthogonal valuse(%)

1.2

1

0.8

0.6

0.4

0.2

0

6

7

8

9

10

Mode number K Fig. 12. Orthogonal values corresponding to different mode numbers.

37

Q. Xiao et al. / Mechanical Systems and Signal Processing 110 (2018) 28–41 Table 2 Energy storage levels corresponding to different mode numbers. K

6

7

8

9

10

Energy storage levels

0.9241

0.8996

0.9345

0.9831

0.9784

Fig. 13. IVMD decomposition result.

5. Experimental analyses 5.1. Experimental conditions A digital speckle pattern interferometry measuring system was established. In this experiment, Model STC-CL152A digital camera (SENTECH, Japan) was used. The CCD target surface size of the CCD camera was 1.27 cm (1/2 in.), h  v = 6.4 mm  4.8 mm. The lens (TAMRON, Japan) with focal distance f of 8 mm (lens dimension, 2/3 in.) were adopted as imaging lens in this experiment. A single longitudinal mode green laser with wavelength of 532 nm was used as the light source. A copper

38

Q. Xiao et al. / Mechanical Systems and Signal Processing 110 (2018) 28–41 100

Mode energies

Threshold values

90 80

Energies

70 60 50 40 30 20 10 0

1

2

3

4

5

6

7

8

9

Mode components Fig. 14. The decomposed mode threshold and mode absolute value.

sheet was used as the tested plate. The optical path for object deformation measurement is shown in Fig. 9. The tested object is shown in Fig. 10. The interferometry speckle patterns of the copper sheet before and after deformation were collected with the CCD camera. Then, an accurate deformation force was applied on the tested cooper sheet, allowing the occurrence of deformation. After that, by opening light source, the obtained interferometry speckle patterns were processed, resulting in the DSPI phase maps, i.e. wrapped phase maps for out-of plane deformation, as shown in Fig. 11. 5.2. IVMD decomposition The phase map of out-of-plane deformation was subject to IVMD decomposition. The size of the phase map for out-ofplane deformation is 900  1500. Then, the mode component after decomposition is 8 6 K 6 9. The orthogonal values corresponding to different K values were computed, as shown in Fig. 12.

Fig. 15. Five denoised mode components.

39

Q. Xiao et al. / Mechanical Systems and Signal Processing 110 (2018) 28–41

Fig. 16. The denoised deformation DSPI phase map for out-of-plane.

As an important index for measuring the image decomposition effect, orthogonal index plays an important role in image decomposition. A higher level of correlation among various image decomposition components would lead to a larger orthogonal value and poorer decomposition effect, and vice versa. As seen from Fig. 12, different K values correspond to different orthogonal values. The corresponding K value is 9 when the orthogonal value is the minimum. Therefore, K = 9 serves as the optimal mode number after decomposition of the DSPI phase map. The energy storage level was used to judge the orthogonal index of various decomposed mode components. The energy storage levels under different mode numbers were computed. The results are shown in Table 2. According to Table 2, when K = 9, the energy storage level after subjecting the image to IVMD decomposition is the highest, indicating the least energy leakage during decomposition, the lowest level of correlation among different mode components, and the best decomposition effect. The DSPI phase map was subject to IVMD decomposition based on the K value. The decomposition results are shown in Fig. 13.

5.3. DSPI phase map denoising As seen from Fig. 13, the improved VMD could effectively decompose the maps and obtain the appropriate mode number. Finally, 9 mode components were obtained. According to the computation principle for the mode threshold mentioned above, the absolute value and the mode threshold were computed, as shown in Fig. 14. As seen from Fig. 14 the absolute value and the mode threshold of the decomposed mode components were computed, respectively. It can be seen that the absolute values of the mode components of u2, u3, u4, u5, and u6 are larger than the mode thresholds. As the five components contain the primary phase information, they were subject to threshold processing. Fig. 15 shows the result of the five mode components subjecting to mode threshold processing. Through analysis of Fig. 15, the components subjecting to mode threshold processing were reconstructed to obtain the denoised phase map. The results of the denoised DSPI phase map for out-of-plane is shown in Fig. 16.

1

13

0.9

12

0.8

11

0.7

SI

PSNR

14

10

0.6

9 0.5 8 0.4 5

6

7

8

9

10

5

6

7

8

The K of numbers

The K of numbers

(a)

(b)

Fig. 17. PSNR and SI under different mode numbers: (a) PSNR and (b) SI.

9

10

40

Q. Xiao et al. / Mechanical Systems and Signal Processing 110 (2018) 28–41 15

0.9

14 0.8

13 0.7

11

SI

PSNR (db)

12

0.6

10 9

0.5

8 0.4

7 6 Proposed method

AFM

MFM

(a)

EMD

VMD

0.3 Proposed method

AFM

MFM

EMD

VMD

(b)

Fig. 18. Comparison of results obtained by three denoising methods in terms of: (a) PSNR and (b) SI.

The peak signal-to-noise ratio (PSNR) and the speckle index (SI) were used to measure the denoising effect. The peak signal-to-noise ratio reflected the map detail-preserving capacity. A higher peak signal-to-noise ratio indicated a higher level of detail preservation. The speckle index reflected the noise suppression capacity. A smaller speckle index indicated higher denoising capacity. For this reason, PSNR and SI under different mode numbers were computed. The results are shown in Fig. 17. Different mode numbers tend to generate false components and influence the accuracy of VMD decomposition, leading to different PSNR and SI. It can be seen that when K is 9, the optimal mode component of the phase map can be obtained, and the accurate decomposition of the phase map with best denoising effect can be achieved. Therefore, the PSNR value was the largest, and the SI value was the smallest. For further comparison, the average filtering method (AFM), median filtering method (MFM), VMD method and EMD methods were used to conduct multiple denoising experiments. The denoised PSNR and SI are shown in Fig. 18. The AFM damages the details of the map while filtering out the noise, then the processed maps become blurred and the information on abrupt phase changes is lost. The MFM has poor smoothing effect. Mode mixing occurs in EMD decomposition, making it impossible to accurately decompose the DSPI phase map. The proposed method in this paper can obtain the optimal mode component by self-adaptive decomposition of the maps, separate the noise components and the noise-free components using the mode threshold, process the noise components, reconstruct the components, and achieve signal noise reduction. Thus, the proposed denoising method processes the DSPI, realizing the largest PSNR and the smallest SI in the processed phase map. 6. Conclusions The imaging process of the digital speckle phase is inevitably subjected to substantial noise interference, leading to great measuring errors. Hence, reducing the noise interference plays a significant role in DSPI. Based on the simulation and actual testing results, the following conclusions are drawn. (1) The mode number is required to be preset before variational mode decomposition. For this reason, the paper proposes an improved VMD method, which can effectively decompose the phase map, obtain the most suitable mode number, and improve the decomposition accuracy and effect. (2) To deal with the shortcomings of the traditional denoised method, the paper proposes an IVMD-based method for reducing noise in the DSPI phase map. The method can self-adaptively compute the mode threshold of the mode components, analyze the mode components using the thresholds, filter out noise, and obtain the denoised DSPI phase map.

Acknowledgments The authors gratefully acknowledge the support provides for this research by the National Natural Science Foundation Project of China (Nos. 51675055 and 61374219), National Key Research and Development Program (2016YFF0101802). We sincerely thank Lianxiang Yang and Sijin Wu, they are from School of Instrumentation Science and Opto-electronics Engineering of Beijing Information Science and Technology University, and they provide the DSPI phase map of this experiment.

Q. Xiao et al. / Mechanical Systems and Signal Processing 110 (2018) 28–41

41

References [1] J. Gryzagoridis, D. Findeis, Impact damage detection on composites using optical NDT techniques, Insight 52 (5) (2010) 248–251. [2] B. Ma, Z. Zhou, Progress and development trends of composite structure evaluation using noncontact nondestructive testing techniques in aviation and aerospace industries, Acta Aeronaut. Astronaut. Sin. 35 (7) (2015) 1787–1803 (in Chinese). [3] L. Yang, X. Xin, L. Zhu, et al, Review of electronic speckle pattern interferometry (ESPI) for three dimensional displacement measurement, Chin. J. Mech. Eng. En. 27 (1) (2014) 1–13. [4] S. Wu, M. Dong, Y. Fang, et al, Universal optical setup for phase-shifting and spatial-carrier digital speckle pattern interferometry, J. Eur. Opt. Soc-Rapid. 12 (1) (2016) 14. [5] Q. Kemao, W. Gao, H. Wang, Windowed Fourier-filtered and quality-guided phase-unwrapping algorithm, Appl. Optics. 47 (29) (2008) 5420–5428. [6] S. Wu, L. Zhu, L. Yang, et al, Spatiotemporal three-dimensional phase unwrapping in digital speckle pattern interferometry, Opt. Lett. 41 (5) (2016) 1050–1053. [7] H. Hosseini, F. Hessar, F. Marvasti, Real-time impulse noise suppression from images using an efficient weighted-average filtering, IEEE. Signal. Proc. Let. 22 (8) (2015) 1050–1054. [8] V. Jayaraj, D. Ebenezer, A new switching-based median filtering scheme and algorithm for removal of high-density salt and pepper noise in images, Eurasip. J. Adv. Sig. Pr. 2010 (1) (2010) 690218. [9] R. Kulkarni, P. Rastogi, Patch-wise denoising of phase fringe patterns based on matrix enhancement, Opt. Laser. Eng. 87 (2016) 191–196. [10] Q. Xiao, J. Li, Z. Bai, et al, A small leak detection method based on VMD adaptive de-noising and ambiguity correlation classification intended for natural gas pipelines, Sensors 16 (12) (2016) 2116. [11] C. Tang, L. Wang, H. Yan, Overview of anisotropic filtering methods based on partial differential equations for electronic speckle pattern interferometry, Appl. Optics 51 (20) (2012) 4916–4926. [12] Y. Qin, J. Chen, H. Fan, The study and application of a new filtering method on electronic speckle pattern interferometry fringe, Opt. Laser. Eng. 39 (4) (2003) 449–456. [13] Q. Kemao, Two-dimensional windowed Fourier transform for fringe pattern analysis: principles, applications and implementations, Opt. Laser Eng. 45 (2) (2007) 304–317. [14] W. Xu, C. Tang, F. Gu, et al, Combination of oriented partial differential equation and shearlet transform for denoising in electronic speckle pattern interferometry fringe patterns, Appl. Optics 56 (10) (2017) 2843–2850. [15] Q. Xiao, J. Li, S. Wu, et al, Adaptive DSPI phase denoising using mutual information and two-dimensional variational mode decomposition, Meas. Sci. Technol. 29 (4) (2017) 1–14. [16] A. Humeau-Heurtier, G. Mahe, P. Abraham, Multi-dimensional complete ensemble empirical mode decomposition with adaptive noise applied to laser speckle contrast images, IEEE. Trans. Med. Imag. 34 (10) (2015) 2103–2117. [17] D. Zosso, K. Dragomiretskiy, A.L. Bertozzi, et al, Two-Dimensional Compact Variational Mode Decomposition, J. Math. Imag. Vis. 58 (2) (2017) 294–320. [18] G. Gu, K. Wang, Y. Wang, et al, Synchronous triple-optical-path digital speckle pattern interferometry with fast discrete curvelet transform for measuring three-dimensional displacements, Opt. Laser. Technol. 80 (2016) 104–111. [19] Y. Zhou, H. Li, A denoising scheme for DSPI fringes based on fast bi-dimensional ensemble empirical mode decomposition and BIMF energy estimation, Mech. Syst. Signal. Pr. 35 (1) (2013) 369–382. [20] C. Wang, Q. Kemao, F. Da, Regenerated phase-shifted sinusoids assisted EMD for adaptive analysis of fringe patterns, Opt. Laser Eng. 87 (2016) 176– 184. [21] A. Humeau-Heurtier, P. Abraham, S. Henni, Bi-dimensional variational mode decomposition of laser speckle contrast imaging data: a clinical approach to critical limb ischemia, Comput. Biol. Med. 86 (2017) 107–112. [22] S. Maheshwari, R.B. Pachori, V. Kanhangad, et al, Iterative variational mode decomposition based automated detection of glaucoma using fundus images, Comput. Biol. Med. 88 (2017) 142–149. [23] Y. Lei, N. Li, J. Lin, et al, Fault diagnosis of rotating machinery based on an adaptive ensemble empirical mode decomposition, Sensors. 13 (12) (2013) 16950–16964. [24] Z. Wu, N.E. Huang, Ensemble empirical mode decomposition: a noise-assisted data analysis method, Adv. Adap. Data Anal. 1 (01) (2009) 1–41. [25] Y. Lei, J. Lin, Z. He, et al, A review on empirical mode decomposition in fault diagnosis of rotating machinery, Mech. Syst. Signal. Pr. 35 (1) (2013) 108– 126. [26] Y. Li, M. Xu, Y. Wei, et al, An improvement EMD method based on the optimized rational Hermite interpolation approach and its application to gear fault diagnosis, Measurement 63 (2015) 330–345. [27] Q. Chen, N. Huang, S. Riemenschneider, et al, A B-spline approach for empirical mode decompositions, Adv. Comput. Math. 24 (1) (2006) 171–195. [28] P. Nguyen, M. Kang, J.M. Kim, et al, Robust condition monitoring of rolling element bearings using de-noising and envelope analysis with signal decomposition techniques, Expert. Syst. Appl. 42 (22) (2015) 9024–9032. [29] P. Nguyen, J.M. Kim, Adaptive ECG denoising using genetic algorithm-based thresholding and ensemble empirical mode decomposition, Inform. Sci. 373 (2016) 499–511. [30] Y. Kopsinis, S. McLaughlin, Development of EMD-based denoising methods inspired by wavelet thresholding, IEEE. Trans. Signal. Pr. 57 (4) (2009) 1351–1362.