UwDAS phase unwrapping method based on 2ndDUI-VMD-EMD

Abstract

The invention discloses a uwDAS phase unwinding method based on 2ndDUI-VMD-EMD. The method comprises the following steps: setting a signal detected by a uwDAS sensing optical fiber as O0 (n); demodulating the O0 (n) based on a 3 * 3 coupler and an arc tangent algorithm to obtain a winding phase signal # imgabs0 # and # imgabs1 #, performing second-order difference on the winding phase signal # imgabs0 # and # imgabs1 #, performing common unwinding on the second-order difference signal # imgabs2 # and # imgabs3 #, performing second-order integration on the differential unwinding signal # imgabs4 # and # imgabs5 #, performing VMD decomposition with the decomposition layer number of k on the integral signal # imgabs6 # and # imgabs7 #, and obtaining k modal components ui (n); the kurtosis of each modal component is calculated, three modal components with the maximum kurtosis are selected for signal reconstruction, an output signal # imgabs8 # is used for carrying out EMD decomposition for extracting a single intrinsic mode function on an output signal # imgabs9 # obtained through VMD decomposition, the single intrinsic mode function IMF1 (n) is obtained, and the IMF1 (n) is the final unwrapping signal y (n). According to the method, the real-time performance of a demodulation system is ensured while the amplitude of an unwinding signal is increased; and the real form of the unwound signal can be more accurately recovered under the condition of low signal-to-noise ratio.

Description

Technical Field

[0001] The present invention relates to the technical field of uwDAS phase demodulation, and in particular to a method based on 2 nd uwDAS phase unwrapping method for DUI-VMD-EMD. Background Art

[0002] Distributed Acoustic Sensing (DAS) system (uwDAS) based on ultra weak fiber Bragg Grating (uwFBG) is widely used due to its advantages such as high sensitivity, long-distance passive measurement and anti-electromagnetic interference. uwDAS is based on phase-sensitive optical time domain reflectometry ( TimeDomain Reflectometry, ) to extract the phase information, The system can realize distributed detection of vibration, strain and other signals by demodulating the phase information of Fizeau interference of reflected light from adjacent uwFBGs. The phase demodulation methods of the system mainly include digital coherent demodulation method, demodulation method based on 3×3 coupler, phase generation carrier demodulation method and I / Q demodulation method.

[0003] During demodulation, the inverse tangent function is first used to limit the complex phase information to between (-π,π) to obtain the warped phase. Although the inverse tangent function can simplify the phase demodulation process, the warped phase output by it cannot directly reflect the real phase change. For example, when the real phase increases from π to 3π, the inverse tangent function will map it to -π, resulting in a phase jump. This jump will seriously affect the analysis and processing of the signal. Therefore, it is necessary to restore the warped phase to a continuous real phase through an unwrapping algorithm. The core idea of ​​the ordinary unwrap algorithm is to detect the jump point by taking the difference between the current sampling point and the previous sampling point. If the difference is within the range of (-π,π), there is no need to compensate the current sampling point; if the difference is greater than or equal to π, the current sampling point and the subsequent points are subtracted by 2π; if the difference is less than or equal to -π, the current sampling point and the subsequent points are added with 2π to obtain the unwrapped signal.

[0004] Based on the common unwrap algorithm, the DUI (Differential-Unwrapping-Integral) algorithm first differentiates the wrapped phase to reduce the signal amplitude, then performs the unwrap operation, and finally restores the original signal by removing the offset through integration and polynomial fitting to achieve unwrap. As the signal amplitude increases, the order of differentiation and integration also increases. For example, the second-order DUI (2 ndThe signal amplitude that can be unwrapped by DUI is larger than that of the first-order DUI. The DUI algorithm can greatly improve the signal amplitude that can be correctly unwrapped compared to the ordinary unwrap algorithm. In addition to polynomial fitting, the method of modal decomposition is also a good choice for removing the offset, such as Empirical Mode Decomposition (EMD) and Variational Modal Decomposition (VMD). While selecting the target signal from the decomposed modal components, it can reduce noise interference.

[0005] Defects and deficiencies of the prior art:

[0006] (1) Due to the limitations of the algorithm principle, the ordinary unwrap algorithm can only correctly recover signals with a phase change less than π between adjacent sampling points. In addition, in the case of low signal-to-noise ratio, due to the large interference of noise, there are many misjudged jump points, and wrong compensation is performed on the signal during unwrapping, resulting in the failure of final unwrapping. This limits the application of the uwDAS system in large dynamic range and low signal-to-noise ratio situations.

[0007] (2) Although the high-order DUI algorithm has improved in unwrapping amplitude, it comes at the cost of increasing the order of difference and integration. A large number of difference and integration operations greatly increase the computational complexity, which limits the real-time performance of system demodulation, and the integration operation will cause a rapid increase in the noise level. Therefore, it is necessary to achieve a larger unwrapping dynamic range and applicability to low signal-to-noise ratio based on the DUI algorithm with as low an order as possible.

[0008] (3) The integration operation in the DUI algorithm will introduce an offset. The existing methods for removing the offset include polynomial fitting. The order of the polynomial is the same as that of DUI. As the order of the polynomial increases, the computational complexity also increases. Moreover, polynomial fitting is sensitive to noise, and noise will cause the fitting curve to deviate from the true trend, resulting in the signal not being correctly unwrapped at low signal-to-noise ratio. Therefore, it is necessary to find a method with better trend fitting effect at low signal-to-noise ratio. Summary of the Invention

[0009] To solve the above technical problems, the present invention provides a uwDAS phase unwrapping method based on 2 nd DUI-VMD-EMD, which ensures the real-time performance of the demodulation system while increasing the unwrapping signal amplitude; the unwrapped signal can more accurately restore its true form in the case of low signal-to-noise ratio; it not only saves computational resources but also can ensure the effect of signal unwrapping.

[0010] The technical solution adopted by the present invention is:

[0011] Based on 2nd uwDAS Phase Unwrapping Method for DUI-VMD-EMD, including the following steps:

[0012] Step 1: Let the signal detected by the uwDAS sensing optical fiber be O 0 (n), where n is the sampling point index, 0 ≤ n ≤ N - 1, and N is the total number of sampling points;

[0013] Step 2: Demodulate O 0 (n) based on a 3×3 coupler and the arctangent algorithm to obtain the wrapped phase signal

[0014] Step 3: Take the second-order difference of the wrapped phase signal to obtain the second-order difference signal ), where n2 is the sampling point index of the second-order difference signal, 0 ≤ n2 ≤ N - 3;

[0015] Step 4: Perform ordinary unwrapping on the second-order difference signal to obtain the unwrapped difference signal

[0016] Step 5: Take the second-order integral of the unwrapped difference signal to obtain the integral signal

[0017] Step 6: Perform VMD decomposition with the decomposition layer number of k on the integral signal to obtain k modal components u i (n), where u i (n) is the i-th modal component, i is the serial number of the modal component, k is the number of VMD decomposition layers, 1 ≤ i ≤ k;

[0018] Step 7: Calculate the kurtosis of each modal component respectively, select the three modal components with the largest kurtosis for signal reconstruction, and output the signal

[0019] Step 8: Perform EMD decomposition for extracting a single intrinsic mode function on the output signal of the VMD decomposition to obtain the single intrinsic mode function IMF 1 (n), and IMF 1 (n) is the final unwrapped signal y(n).

[0020] In the said Step 2, demodulating O 0 (n) based on a 3×3 coupler and the arctangent algorithm includes the following steps:

[0021] S2.1: The 3×3 coupler outputs three signals O with a phase difference of 1(n), O 2 (n) and O 3 (n), the expressions are respectively:

[0022] O 1 (n) = cos(O 0 (n));

[0023]

[0024]

[0025] where: cos(·) is the cosine function, O 1 (n), O 2 (n), O 3 the mean value av(n) of (n) is:

[0026]

[0027] S2.2: Use the arctangent function to obtain the wrapped phase signal Its calculation formula is:

[0028]

[0029] where: arctan(·) is the four - quadrant arctangent function.

[0030] In the said step 3, for the wrapped phase signal perform a second - order difference, including the following steps:

[0031] S3.1: Perform a first - order difference on the wrapped phase signal to obtain a first - order difference signal Its calculation formula is:

[0032]

[0033] where: n1 is the sampling point index of the first - order difference signal, 0 ≤ n1 ≤ N - 2; represents the value of the (n1 + 1)-th sampling point of the wrapped signal, represents the value of the n1 - th sampling point of the wrapped signal.

[0034] S3.2: Perform a first - order difference on the first - order difference signal again to obtain a second - order difference signal Its calculation formula is:

[0035]

[0036] where: represents the value of the (n2 + 1)-th sampling point of the first - order difference signal, Represents the value of the n2-th sampling point of the first-order difference signal.

[0037] Step 4 includes the following steps:

[0038] S4.1: Calculate the phase compensation factor m, and its calculation formula is:

[0039]

[0040] Where: Represents the value of the (n2 - 1)-th sampling point of the second-order difference signal, and round(i) represents the rounding operation.

[0041] S4.2: Calculate the unwrapped signal of the difference from the phase compensation factor m Its calculation formula is:

[0042]

[0043] Step 5 includes the following steps:

[0044] S5.1: Perform a first-order integration on the unwrapped signal of the difference To obtain the first-order integration signal Its calculation formula is:

[0045]

[0046] Where: Represents the value of the initial sampling point of the first-order difference signal; Represents the value of the q1-th sampling point of the unwrapped signal of the difference; q1 is the index in the first-order integration summation operation, q1 = 0, 1,..., n1 - 1;

[0047] S5.2: Perform a first-order integration on the first-order integration signal To obtain the final second-order integration signal Its calculation formula is:

[0048]

[0049] Where: Represents the value of the initial sampling point of the wrapped phase signal, Represents the value of the q2-th sampling point of the first-order integration signal; q2 is the index in the second-order integration summation operation, q2 = 0, 1,..., n - 1;

[0050] Step 6 includes the following steps:

[0051] S6.1: Initialize the following parameters: correlation coefficient r = 0; VMD decomposition layer number k = 2;

[0052] S6.2: Initialize the spectrum of the modal components The central frequency of each modal component and the Lagrange multiplier where: is the initial value of the spectrum of the i-th modal component, ω is the frequency index, 0 ≤ ω ≤ N / 2; is the initial value of the central frequency of the i-th modal component, and initialize the number of iterations n of the VMD decomposition d = 1;

[0053] S6.3: Update and The calculation formula is:

[0054]

[0055]

[0056]

[0057] where: is the spectrum of the i-th modal component in the (n d + 1)-th iteration; is the frequency-domain representation of;

[0058] is the spectrum of the i'-th modal component in the (n d + 1)-th iteration; i' is the sequence number of the modal component when updating the mode, 1 ≤ i' ≤ k; is the Lagrange multiplier in the n d -th iteration; α is the balance parameter of the fidelity constraint; is the central frequency of the i-th modal component in the n d -th iteration; is the central frequency of the i-th modal component in the (n d + 1)-th iteration;

[0059] is the Lagrange multiplier in the (n d + 1)-th iteration; τ is the update parameter of the Lagrange multiplier; denotes integration with respect to ω from 0 to ∞; |·| denotes the absolute value operation;

[0060] S6.4: Judge whether the following convergence condition is reached:

[0061]

[0062] where: denotes the spectrum of the i-th modal component in the n d -th iteration; ε is the convergence tolerance; Denotes the square of the 2-norm. If the convergence condition is reached, the spectra of the current k modal components are obtained. And perform S6.5. Is the spectrum of the i-th modal component; otherwise, increment n d by 1 and repeat S6.3.

[0063] S6.5: Through Reconstruct the time-domain modal component u i (n), including the following steps:

[0064] S6.5.1: Rearrange the frequencies of to obtain Denotes the spectrum of the i-th modal component after frequency rearrangement, and its calculation formula is:

[0065]

[0066] Where: IFFTShift(·) means shifting the low-frequency part of the signal spectrum to both ends of the spectrum, so that the spectrum is centrosymmetric.

[0067] S6.5.2 Perform the inverse fast Fourier transform on and take the real part of the transform result to obtain the i-th state modal component u i (n), and its calculation formula is:

[0068]

[0069] Where: Re(·) means taking the real part operation, j is the imaginary unit, and e is the base of the natural logarithm.

[0070] S6.6: Calculate the trend trend(n) fitted by the VMD decomposition for the integral signal , and its calculation formula is:

[0071]

[0072] Where: u k (n) is the k-th modal component.

[0073] S6.7: Calculate the correlation coefficient r between the integral signal and the trend trend(n), and the calculation formula is:

[0074]

[0075] Where: E[·] means taking the mean operation.

[0076] S6.8: Judge whether the following loop condition is reached

[0077] r < th && k < kmax;

[0078] Where: th is the threshold of the correlation coefficient, kmax is the maximum number of layers of VMD decomposition, && represents the logical operator "AND". If the loop execution condition is met, the value of k is incremented by 1 and steps S6.2 to S6.7 are repeated; otherwise, step S7 is performed.

[0079] The said step 7 includes the following steps:

[0080] S7.1: Calculate the kurtosis Kurt(i) of the i-th mode component u i (n), and its calculation formula is:

[0081]

[0082] Where: u i (n) 4 represents the fourth power of the i-th mode component, and u i (n) 2 represents the square of the i-th mode component.

[0083] S7.2: Arrange the mode components in descending order of kurtosis, and select the sum of the three mode components corresponding to the three largest kurtosis values as the output signal of VMD decomposition Its calculation formula is:

[0084]

[0085] Where: I 1 、I 2 、I 3 are the indices of the serial numbers of the three mode components with the largest kurtosis respectively, and are the three mode components with the largest kurtosis respectively.

[0086] The said step 8 includes the following steps:

[0087] S8.1: Initialize the candidate intrinsic mode function and the number of iterations t of EMD decomposition to 1.

[0088] S8.2: Update the decomposed signal r t (n) in the t-th iteration, and its calculation formula is:

[0089] r t (n) = imf t-1 (n);

[0090] Where, imf t-1 (n) is the candidate intrinsic mode function in the (t - 1)-th iteration.

[0091] S8.3: Find the decomposed signal r tAll the maximum and minimum points of (n). The judgment condition for the maximum point is:

[0092] r t (q) > r t (q - 1) && r t (q) > r t (q + 1);

[0093] Where: q is the index of the sampling point in r t (n), 1 ≤ q ≤ N - 2; r t (q), r t (q - 1), r t (q + 1) are the values of the q-th, (q - 1)-th, and (q + 1)-th sampling points of the decomposed signal in the t-th iteration respectively; if r t (q) satisfies the above judgment condition, then r t (q) is a maximum point;

[0094] The judgment condition for the minimum point is:

[0095] r t (q) < r t (q - 1) && r t (q) < r t (q + 1);

[0096] If r t (q) satisfies the above judgment condition, then r t (q) is a minimum point.

[0097] S8.4: Use cubic spline interpolation to interpolate the maximum and minimum points to obtain the upper envelope e max (n) and the lower envelope e min (n). The calculation formula for the upper envelope e max (n) is:

[0098] e max (n) = a l (q - q l_begin ) 3 + b l (q - q l_begin ) 2 + c l (q - q l_begin ) + d l ;

[0099] Where: q l_begin is the starting point of the l-th upper envelope interpolation interval; l is the index of the upper envelope interpolation interval, l = 1, 2, …, L - 1, L, and L is the total number of upper envelope interpolation intervals; a l , b l , cl , d l is the interpolation coefficient of the l-th upper envelope interval.

[0100] Lower envelope e min (n) is calculated as follows:

[0101]

[0102] Where: is the starting point of the l d -th lower envelope interpolation interval; l d is the index of the lower envelope interpolation interval, l d = 1, 2, …, L d -1, L d , L d is the total number of lower envelope interpolation intervals; is the interpolation coefficient of the l d -th lower envelope interval.

[0103] S8.5: Calculate the local mean m t (n) of the upper and lower envelopes of the decomposed signal r t (n) at the t-th iteration, and its calculation formula is:

[0104]

[0105] S8.6: Calculate the candidate intrinsic mode function imf t (n) at the t-th iteration, and its calculation formula is:

[0106] imf t (n) = r t (n) - m t (n);

[0107] S8.7: Calculate the number of extreme points N t and the number of zero crossings N ext of imf zero (n), including the following steps:

[0108] S8.7.1: Calculate the number of maximum points N max , and its calculation formula is:

[0109]

[0110] Where: q' is the index of the sampling point in imf t (n), 1 ≤ q' ≤ N - 1; imf t (q'), imf t (q' - 1), imf t(q'+1) is the value of the q', q'-1, and q'+1 sampling points of the candidate intrinsic mode function in the t-th iteration; (·) is an indicator function that takes the value of 1 when the condition is satisfied and 0 otherwise.

[0111] S8.7.2: Calculate the number of minimum points N min The calculation formula is:

[0112]

[0113] S8.7.3: imf t (n) The number of extreme points N ext The calculation formula is:

[0114] N ext = N max + N min ;

[0115] S8.7.4: imf t (n) The number of zero-crossing points N zero The calculation formula is:

[0116]

[0117] S8.8: Calculate the upper envelope h t (n) and the lower envelope h max (n) of imf min The calculation formula for the upper envelope is:

[0118] h max (n) = A s (q'-q' s_begin ) 3 + B s (q'-q' s_begin ) 2 + C s (q'-q' s_begin ) + D s ;

[0119] Where: q' s_begin is the starting point of the s-th upper envelope interpolation interval; s is the index of the upper envelope interpolation interval, s = 1, 2,..., N max - 1; A s , B s , C s , D s are the interpolation coefficients of the s-th upper envelope interval.

[0120] The calculation formula for the lower envelope is:

[0121]

[0122] Wherein: is the starting point of the sth d lower envelope interpolation interval; s d is the index of the lower envelope interpolation interval, s d = 1, 2, …, N min - 1; is the interpolation coefficient of the sth d lower envelope interval.

[0123] S8.9: Calculate the local mean M t (n) of the upper and lower envelopes of the t-th iteration imf t (n), and its calculation formula is:

[0124]

[0125] S8.10: Determine whether the following condition holds:

[0126] |N ext - N zero | ≤ 1 && max(|M t (n)|) < δ;

[0127] Wherein: max(·) is the maximum value operation, δ is the calculation tolerance. If the condition holds, go to step S8.11; otherwise, increment the value of t by 1 and repeat steps S8.2 - S8.9.

[0128] S8.11: The expression of the single intrinsic mode function IMF 1 (n) obtained by performing EMD decomposition is:

[0129] IMF 1 (n) = imf t (n).

[0130] The uwDAS phase unwrapping method based on 2 nd DUI - VMD - EMD of the present invention has the following technical effects:

[0131] 1) The present invention uses 2 nd DUI, which reduces the amplitude of the unwrapping input signal through second - order difference, thereby expanding the dynamic range of the phase unwrapping amplitude, breaking the limitation that the ordinary unwrap algorithm can only correctly unwrap signals with a phase change less than π between adjacent sampling points. Moreover, compared with higher - order processing methods, the computational load it brings is within the capacity of a general demodulation system, ensuring the real - time performance of the demodulation system while increasing the amplitude of the unwrapped signal.

[0132] 2) The present invention combines VMD and EMD to remove the offset caused by integral operation. First, the optimal number of layers for VMD decomposition is determined through the correlation coefficient between the integral signal and the offset, and the integral signal is decomposed into multiple different modal components by VMD with the obtained optimal number of layers, which represent different frequency components in the signal. Then, the kurtosis index is used to select the required signal modes, which can remove most of the high-frequency noise and low-frequency offset. Using EMD decomposition can further remove noise and offset, enabling the unwound signal to more accurately restore its true form under low signal-to-noise ratio conditions.

[0133] 3) Since EMD is prone to modal aliasing, the present invention first uses VMD to process the integral signal before EMD. On the one hand, it can improve the quality of the input signal for EMD decomposition, thus effectively avoiding the modal aliasing phenomenon of EMD; on the other hand, due to the preprocessing of VMD, EMD only needs to perform decomposition to extract a single intrinsic mode function, that is, the decomposition result of EMD only has one intrinsic mode function and the offset to be removed, which not only saves computing resources but also can ensure the effect of signal unwinding.

[0134] 4) Both VMD and EMD are adaptive methods, and the present invention can be adjusted according to the characteristics of the signal. For complex signals, ordinary unwrap algorithms and the DUI algorithm using polynomial fitting to remove the offset of the integral signal often have poor effects, while the combination of VMD and EMD in the present invention can better adapt to these complex signals. BRIEF DESCRIPTION OF THE DRAWINGS

[0135] The present invention will be further described below in conjunction with the drawings and examples;

[0136] Figure 1 is a flowchart of the present invention.

[0137] Figure 2 is the modal decomposition result of VMD in Example 1.

[0138] Figure 3 is the signal waveform diagram reconstructed from the three modal components with the maximum kurtosis in Example 1.

[0139] Figure 4 is the modal decomposition result of EMD in Example 1.

[0140] Figure 5 is the comparison diagram of the original signal waveform, ordinary unwrap result, 2 nd DUI result and 2 nd DUI-VMD-EMD result in Example 1.

[0141] Figure 6 is the ordinary unwrap result, 2 nd DUI result and 2nd Curve graph showing the variation of the similarity between the DUI-VMD-EMD result and the original signal with the signal-to-noise ratio.

[0142] Figure 7 For the ordinary de - winding result in Example 2, 2 nd DUI result and 2 nd Curve graph showing the variation of the similarity between the DUI-VMD-EMD result and the original signal with the amplitude. Detailed implementation manners

[0143] Example 1:

[0144] The signal collected by the sensing optical fiber is simulated using a sine signal with an amplitude of 10 rad, a signal frequency of 1 kHz, a sampling frequency of 10 kHz, and a signal-to-noise ratio of 35 dB to observe the phase de - winding effect of the present invention. The parameters of the present invention are set as follows:

[0145] The balance parameter α of the fidelity constraint = 2000;

[0146] The update parameter τ of the Lagrange multiplier = 0;

[0147] The convergence tolerance ε = 1e - 6;

[0148] The threshold th of the correlation coefficient = 0.999999992;

[0149] The maximum number of layers kmax for VMD decomposition = 50;

[0150] The calculation tolerance δ = 1e - 6.

[0151] Use the 2 nd DUI-VMD-EMD method proposed by the present invention to de - wind the signal, and compare the de - winding results with those of the ordinary de - winding algorithm and the 2 nd DUI algorithm:

[0152] Observe the signal waveforms at some process nodes when the present invention de - winds the signal. It can be seen from Figure 2 that by the correlation coefficient of the integral signal and the offset, the optimal number of decomposition layers for VMD decomposition is determined to be 16 at this time, that is, decomposed into 16 different modal components. Observing these modal components, it can be seen that not all modes have a strong correlation with the collected original signal. Therefore, it is necessary to use kurtosis to screen these modes. According to the calculation result of kurtosis, select the 14th, 15th, and 16th modal components, namely u 14 (n), u 15 (n), u 16 (n) for reconstruction. It can be seen from Figure 3It can be seen that the reconstructed signal waveforms of the three modal components with the largest kurtosis selected are relatively close to the original signal, but there is a certain offset. Therefore, it is necessary to further use EMD to remove the offset. Since the present invention performs EMD decomposition to extract a single intrinsic mode function, as Figure 4 can be seen, the decomposition result of EMD is 1 intrinsic mode component and a residue. Among them, the intrinsic mode component removes the offset on the basis of the Figure 3 result, that is, removes the residue part in the EMD decomposition, making the signal waveform closer to the original signal. In addition, the residue result of the EMD decomposition does not show an obvious pattern, which is an effect that cannot be achieved by the polynomial fitting of the DUI algorithm. This also reflects the advantage of the combination of VMD and EMD in removing the offset.

[0153] As Figure 5 can be seen, under the same conditions, the results of ordinary deconvolution and 2 nd DUI are completely distorted, while the present invention can well achieve deconvolution. Although the result of the present invention is somewhat distorted compared with the original waveform, the waveforms of the two are relatively consistent, and compared with the complete distortion of the results of ordinary deconvolution and 2 nd DUI, it has been greatly improved. This shows that the deconvolution effect of the present invention is better than that of ordinary deconvolution and 2 nd DUI.

[0154] Example 2:

[0155] The original analog signal collected by the sensing optical fiber has the same parameter settings as those in Example 1 except for the amplitude and signal-to-noise ratio.

[0156] First, keep the amplitude of the analog signal always at 9.25 rad, and change the signal-to-noise ratio within the range of [22, 42] dB. The applicability of different deconvolution methods to the signal-to-noise ratio is studied through the similarity between the deconvolution result and the original analog signal. The closer the similarity is to 1, the better. As Figure 6 can be seen, the similarity between the result of the present invention and the original signal is generally higher than that of the other two deconvolution methods as the signal-to-noise ratio changes. For ordinary deconvolution, the amplitude of 9.25 rad at this time has exceeded the maximum amplitude that can be deconvolved, so no matter how the signal-to-noise ratio changes, it cannot be correctly deconvolved; for 2 nd DUI, the deconvolution fails when the signal-to-noise ratio is lower than 38 dB; when the signal-to-noise ratio of the present invention is 30 dB and above, its deconvolution result maintains a high similarity with the original signal. When the signal-to-noise ratio is 25 dB, the similarity is 0.7144, and when the signal-to-noise ratio is 28 dB, the similarity is 0.8103. Compared with the other two methods, the deconvolution effect is significantly improved, indicating that the present invention is more suitable for the case of low signal-to-noise ratio and broadens the applicable range of the signal-to-noise ratio.

[0157] Keep the signal-to-noise ratio of the original analog signal in simulation at 35 dB all the time, change the amplitude within the range of [3, 12] rad, and study the amplitude dynamic range of different unwrapping methods through the similarity between the unwrapping result and the original analog signal. From Figure 7 It can be seen that when the signal-to-noise ratio is 35 dB, the maximum amplitude that can be unwrapped by the present invention is 11.75 rad. At this time, the similarity between the unwrapping result and the original signal is 0.8896. 2 nd The maximum amplitude that DUI can unwrap is 8 rad, and the maximum amplitude that the ordinary unwrapping method can unwrap is about 4 rad. The upper limit of the amplitude that the present invention can unwrap is higher than the other two unwrapping methods, indicating that the amplitude dynamic range of the unwrapping of the present invention has been improved.

[0158] For the signal collected by the present invention, first use a 3×3 coupler and the arctangent algorithm to limit the phase between (-π, π) to achieve phase wrapping, then perform a second-order difference on the wrapped signal and perform ordinary unwrapping on the difference signal to expand the range of the phase unwrapping amplitude, perform a second-order integration on the unwrapped difference signal to restore the original signal, and use VMD and EMD to remove the offset introduced by the integration. Determine the optimal number of layers of VMD decomposition through the correlation coefficient between the integrated signal and the offset, and use kurtosis to select the required signal mode. Finally, perform EMD decomposition on the selected mode components to extract a single intrinsic mode function, and this single intrinsic mode function is the final unwrapped signal. The present invention has been greatly improved in terms of dynamic range and applicability to low signal-to-noise ratio.

Claims

1. Based on 2 nd The uwDAS phase unwrapping method of DUI-VMD-EMD is characterized by The following steps are involved: Step 1: Assume that the signal detected by the uwDAS sensing fiber is O0(n), where n is the sampling point index, 0≤n≤N-1, and N is the total number of sampling points; Step 2: Demodulate O0(n) based on the 3×3 coupler and the inverse tangent algorithm to obtain the warped phase signal Step 3: Phase the Wrapped Signal Do 2nd order difference to get 2nd order difference signal Wherein, n2 is the sampling point index of the second-order differential signal, 0≤n2≤N-3; Step 4: 2nd order differential signal Perform normal unwinding to obtain the differential unwinding signal Step 5: Unwrap the differential signal Perform 2nd-order integration to obtain the integral signal Step 6: Integrate the signal Perform VMD decomposition with a decomposition layer of k to obtain k modal components u i (n), where u i (n) is the i-th modal component, i is the ordinal number of the modal component, k is the number of VMD decomposition levels, 1≤i≤k; Step 7: Calculate the kurtosis of each modal component separately, select the three modal components with the largest kurtosis for signal reconstruction, and output the signal Step 8: Decompose the output signal of VMD Perform EMD decomposition to extract a single intrinsic mode function and obtain a single intrinsic mode function IMF1(n). IMF1(n) is the final deconvolved signal y(n).

2. According to claim 1 based on 2 nd The uwDAS phase unwrapping method of DUI-VMD-EMD is characterized by: In step 2, O0(n) is demodulated based on a 3×3 coupler and an inverse tangent algorithm, including the following steps: S2.1: The three phase differences of the 3×3 coupler output are The signals O1(n), O2(n) and O3(n) are expressed as: O1(n)=cos(O0(n)); Where: cos(·) is the cosine function, and the mean value av(n) of O1(n), O2(n), and O3(n) is: S2.2: Using the inverse tangent function to obtain the warped phase signal The calculation formula is: Where: arctan(·) is the four-quadrant inverse tangent function.

3. According to claim 1 based on 2 nd The uwDAS phase unwrapping method of DUI-VMD-EMD is characterized by: In step 3, the warped phase signal Doing second-order differences includes the following steps: S3.1: For the winding phase signal Do first-order difference to get first-order difference signal The calculation formula is: Where: n1 is the sampling point index of the first-order differential signal, 0≤n1≤N-2; Represents the value of the n1+1th sampling point of the wrap signal, Represents the value of the n1th sampling point of the winding signal; S3.2: For 1st order differential signal Then do the first-order difference to get the second-order difference signal The calculation formula is: in: Represents the value of the n2+1th sampling point of the 1st-order differential signal, Represents the value of the n2th sampling point of the first-order difference signal.

4. According to claim 1 based on 2 nd The uwDAS phase unwrapping method of DUI-VMD-EMD is characterized by: The step 4 comprises the following steps: S4.1: Calculate the phase compensation factor m, which is calculated as follows: in: It represents the value of the n2-1th sampling point of the second-order differential signal, and round(i) represents the rounding operation; S4.2: Calculate the differential unwrapped signal by the phase compensation factor m The calculation formula is:

5. According to claim 1 based on 2 nd The uwDAS phase unwrapping method of DUI-VMD-EMD is characterized by: Step 5 The following steps are involved: S5.1: Unwrapping the differential signal Perform 1st-order integration to obtain 1st-order integrated signal The calculation formula is: in: Represents the value of the initial sampling point of the first-order differential signal; Represents the value of the q1th sampling point of the differential dewrapped signal; q1 is the index in the 1st order integral summation operation, q1=0,1,…,n1-1; S5.2: For the 1st order integrated signal Do the first-order integration again to get the final second-order integrated signal The calculation formula is: in: Represents the value of the initial sampling point of the wrapped phase signal, Represents the value of the q2th sampling point of the 1st-order integrated signal; q2 is the index in the 2nd-order integral summation operation, q2 = 0, 1, …, n-1.

6. According to claim 1 based on 2 nd The uwDAS phase unwrapping method of DUI-VMD-EMD is characterized by: The step 6 comprises the following steps: S6.1: Initialize the following parameters: correlation coefficient r = 0; VMD decomposition level k = 2; S6.2: Initialize the spectra of the modal components The center frequency of each modal component and Lagrange multipliers in: is the initial value of the spectrum of the i-th modal component, ω is the frequency index, 0≤ω≤N / 2; is the initial value of the center frequency of the i-th modal component, and the number of iterations n of initializing VMD decomposition d =1; S6.3: Update and The calculation formula is: in: For nth d +1 frequency spectrum of the ith modal component of iteration; for Frequency domain representation of ; For nth d +1 iteration spectrum of the i'th modal component; i' is the sequence number of the modal component when updating the mode, 1≤i'≤k; For nth d The Lagrange multiplier of the iteration; α is the balance parameter of the fidelity constraint; For nth d The center frequency of the i-th modal component of the iteration; For nth d +1 iteration center frequency of the ith modal component; For nth d +1 iteration of the Lagrange multiplier; τ is the update parameter of the Lagrange multiplier; Indicates the integration of ω from 0 to ∞; |·| indicates the absolute value operation; S6.4: Determine whether the following convergence conditions are met: in: Indicates the nth d The spectrum of the ith modal component of the iteration; ε is the convergence tolerance; Represents the square of the 2-norm. If the convergence condition is met, the spectrum of the current k modal components is obtained. And proceed to S6.5, is the spectrum of the ith modal component; otherwise, n d The value of is increased by 1 and S6.3 is repeated; S6.5: Pass Reconstruct the time domain modal component u i (n); S6.6: Calculate the integrated signal by VMD decomposition The calculation formula of the fitted trend trend(n) is: Where: u k (n) is the kth modal component; S6.7: Calculate the integrated signal The correlation coefficient r with trend(n) is calculated as: Where: E[·] represents the averaging operation; S6.8: Determine whether the following loop conditions are met r <th&&k<kmax; Wherein: th is the threshold value of the correlation coefficient, kmax is the maximum number of layers of VMD decomposition, && represents the logical operator "and", if the loop condition is met, the value of k is increased by 1 and steps S6.2 to S6.7 are repeated, otherwise step S7 is performed.

7. According to claim 6 based on 2 nd The uwDAS phase unwrapping method of DUI-VMD-EMD is characterized by: In S6.5, by Reconstruct the time domain modal component u i (n), comprising the following steps: S6.5.1: Rearrange the frequencies to get represents the spectrum of the i-th modal component after frequency rearrangement, and its calculation formula is: Among them: IFFTShift(·) means moving the low-frequency part of the signal spectrum to both ends of the spectrum, so that the spectrum center is symmetrical; S6.5.2 Pair Perform an inverse fast Fourier transform and take the real part of the transform result to obtain the i-th state mode component u i (n), which is calculated as follows: Where: Re(·) represents the real part operation, j is the imaginary unit, and e is the base of the natural logarithm.

8. According to claim 1 based on 2 nd The uwDAS phase unwrapping method of DUI-VMD-EMD is characterized by: The step 7 comprises the following steps: S7.1: Calculate the i-th modal component u i The kurtosis Kurt(i) of (n) is calculated as follows: Where: u i (n) 4 represents the fourth power of the ith modal component, u i (n) 2 represents the square of the i-th modal component; S7.2: The modal components are arranged in descending order of kurtosis, and the sum of the modal components corresponding to the first three largest kurtosis is selected as the output signal of VMD decomposition. The calculation formula is: Among them: I1, I2, I3 are the indexes of the three modal components with the largest kurtosis, are the three modal components with the largest kurtosis.

9. According to claim 1 based on 2 nd The uwDAS phase unwrapping method of DUI-VMD-EMD is characterized by: The step 8 comprises the following steps: S8.1: Initialize candidate eigenmode functions and EMD decomposition iteration number t = 1; S8.2: Update the decomposed signal r in the tth iteration t (n), which is calculated as follows: r t (n)=imf t-1 (n); Among them, IMF t-1 (n) is the candidate eigenmode function in the t-1th iteration; S8.3: Find the decomposed signal r in the tth iteration t All the maximum and minimum points of (n), the judgment condition of the maximum point is: r t (q)>r t (q-1)&&r t (q)>r t (q+1); Where: q is r t (n) is the index of the sampling point, 1≤q≤N-2; r t (q), r t (q-1), r t (q+1) are the values ​​of the qth, q-1th, and q+1th sampling points of the decomposed signal in the tth iteration respectively; if r t (q) If the above judgment conditions are met, then r t (q) is the maximum point; The judgment conditions for the minimum point are: r t (q)<r t (q-1)&&r t (q)<r t (q+1); If r t (q) If the above judgment conditions are met, then r t (q) is the minimum point; S8.4: Use cubic spline interpolation to interpolate the maximum and minimum points to obtain the upper envelope e max (n) and the lower envelope e min (n), upper envelope e max The calculation formula for (n) is: e max (n)=a l (q-q l_begin ) 3 +b l (q-q l_begin ) 2 +c l (q-q l_begin )+d l ; Where: q l_begin is the starting point of the lth upper envelope interpolation interval; l is the index of the upper envelope interpolation interval, l = 1, 2, ..., L-1, L, L is the total number of upper envelope interpolation intervals; a l 、b l 、c l ,d l is the interpolation coefficient of the lth upper envelope interval; Lower envelope e min The calculation formula for (n) is: in: For the first d The starting point of the lower envelope interpolation interval; l d is the index of the lower envelope interpolation interval, l d =1,2,…,L d -1,L d , L d is the total number of interpolation intervals of the lower envelope; For the first d Interpolation coefficients for the lower envelope intervals; S8.5: Find the decomposition signal r at the tth iteration t (n) Local mean of upper and lower envelopes m t (n), which is calculated as follows: S8.6: Find the candidate intrinsic mode function imf for the tth iteration t (n), which is calculated as follows: imf t (n)=r t (n)-m t (n); S8.7: Find imf t (n) The number of extreme points N ext and the number of zero crossings N zero ; S8.8: Find imf t The upper envelope of (n) h max (n) and the lower envelope h min (n), the upper envelope is calculated as: h max (n)=A s (q'-q' s_begin ) 3 +B s (q'-q' s_begin ) 2 +C s (q'-q' s_begin )+D s ; Where: q' s_begin is the starting point of the sth upper envelope interpolation interval; s is the index of the upper envelope interpolation interval, s=1,2,…,N max -1; A s , B s , C s , D s is the interpolation coefficient of the sth upper envelope interval; The calculation formula for the lower envelope is: in: For the s d The starting point of the lower envelope interpolation interval; s d is the index of the lower envelope interpolation interval, s d =1,2,…,N min -1; For the s d Interpolation coefficients for the lower envelope intervals; S8.9: Find the tth iteration imf t (n) Local mean of upper and lower envelopes M t (n), which is calculated as follows: S8.10: Determine whether the following conditions are met: |N ext -N zero |≤1&&max(|M t (n)|)<δ; Where: max(·) is the maximum value operation, δ is the calculation tolerance, if the condition is met, proceed to step S8.11; otherwise, the value of t is increased by 1 and steps S8.2 to S8.9 are repeated; S8.11: The expression of the single intrinsic mode function IMF1(n) obtained by EMD decomposition is: IMF1(n)=imf t (n)。 10. According to claim 9 based on 2 nd The uwDAS phase unwrapping method of DUI-VMD-EMD is characterized by: In S8.7, find imf t (n) The number of extreme points N ext and the number of zero crossings N zero , including the following steps: S8.7.1: Find the number of maximum points N max , and its calculation formula is: Where: q' is imf t (n) is the index of the sampling point, 1≤q'≤N-1; imf t (q'), imf t (q'-1), imf t (q'+1) is the value of the q', q'-1, q'+1 sampling points of the candidate intrinsic mode function in the tth iteration; is an indicator function, which takes the value 1 when the condition is met, otherwise it takes the value 0; S8.7.2: Find the number of minimum points N min The calculation formula is: S8.7.3: imf t (n) Number of extreme points N ext The calculation formula is: N ext =N max +N min ; S8.7.4: imf t (n) Number of zero crossings N zero The calculation formula is:

Images