Non-contact heart rate variability detection method based on improved variational mode decomposition and mode selection criterion

Abstract

The invention discloses an HRV estimation method of a millimeter wave radar. The method comprises the following steps: firstly, positioning a human body by using the range resolution characteristic of an FMCW radar, and preprocessing noise and interference existing in an echo signal; secondly, a grey wolf optimization algorithm is introduced, the maximum mutual information coefficient and the minimum reconstruction error serve as fitness functions, a traditional variational mode decomposition method is optimized, and self-adaptive separation of vital sign signals is achieved. Meanwhile, an optimal heartbeat component is selected from decomposed modals based on a modal selection criterion of a correlation coefficient. And finally, extracting a time parameter corresponding to the heartbeat wave crest, further obtaining IBI, and calculating the time-frequency domain characteristics of the HRV through the IBI. According to the method, the heartbeat signals are effectively extracted through the improved mode decomposition method and the mode selection criterion, and then high-precision estimation of the heart rate variability is achieved.

Description

Technical Field

[0001] The invention belongs to human vital sign detection technology, and specifically relates to a non-contact heart rate variability detection method based on improved variational mode decomposition and mode selection criteria. Background Art

[0002] Heart rate is the number of heartbeats per minute. A healthy heartbeat is not a regular beat like a metronome, but each beat has slight differences. Heart rate variability (HRV) reflects the difference in the time intervals between adjacent heartbeats, called the interbeat interval (IBI). HRV has been proven to be an important indicator for assessing overall heart health and the state of the autonomic nervous system responsible for regulating heart activity. In addition, HRV detection has been applied in sleep monitoring, stress detection, emotion recognition and other fields. How to accurately estimate HRV is a hot issue and research trend in the field of contactless vital signs detection using millimeter-wave radar.

[0003] Traditional HRV is usually estimated by obtaining IBI through contact devices such as electrocardiograms or photoplethysmography. Although contact detection has high stability and accuracy, long-term contact with the skin will reduce user comfort and bring inconvenience, especially for special groups such as burn patients and infants. Compared with traditional sensing devices, wireless signals do not require users to wear any sensors, and can complete monitoring without being noticed. This medical health monitoring system based on wireless signals will not interfere with the user's daily life, and effectively solves the problem of user compliance. Therefore, non-contact detection technology is particularly important. In the field of non-contact detection, compared with wireless devices such as Wi-Fi and video devices, frequency modulated continuous wave radar is more suitable for vital sign detection because of its suitable distance resolution, micro-Doppler detection capability, simple system structure and good anti-interference ability.

[0004] Accurately estimating HRV is more challenging than simply estimating HR because it requires accurate acquisition of the heartbeat waveform in order to calculate the time length of each heartbeat cycle. Therefore, the present invention discloses a HRV estimation method for millimeter wave radar. First, the distance resolution characteristics of the FMCW radar are used to locate the human body, and the noise and interference in the echo signal are preprocessed. Secondly, the combination of the maximum mutual information coefficient (MIC) and the minimum reconstruction error (MRE) is used as the fitness function, and the Grey Wolf Optimization (GWO) algorithm is used to search for the optimal influencing parameter combination of the variational mode decomposition algorithm (VMD) to determine the optimal value of the penalty factor α and the number of modes k for variational mode extraction. Then, the improved VMD algorithm is used to separate the vital sign signal, and a mode selection criterion based on the correlation coefficient is proposed to select the optimal heartbeat signal component from the decomposed mode. Finally, the time parameter corresponding to the peak value of the heartbeat wave is extracted, and the IBI is further obtained, and the time-frequency domain characteristics of HRV are calculated by the IBI. The patent of this invention effectively reduces the interference and noise in the radar signal by combining the improved modal decomposition method and the modal selection criterion, thereby achieving high-precision estimation of heart rate variability. Summary of the invention

[0005] The purpose of the present invention is to provide an accurate and robust method for detecting human heart rate variability, which utilizes millimeter wave radar to capture human vital sign signals, extracts heartbeat signals from radar signals through a series of signal processing algorithms, and further extracts time domain and frequency domain features of HRV from the heartbeat waveform, thereby realizing the detection of human heart rate variability.

[0006] The present invention discloses a non-contact heart rate variability detection method based on improved variational mode decomposition and mode selection criteria. The method comprises the following steps:

[0007] Step 1: Use frequency modulated continuous wave (FMCW) radar to collect human vital signs radar signals. Use the distance-FFT method to obtain human distance information and extract the corresponding echo signal. After performing arc tangent demodulation and a series of pre-processing noise reduction reconstruction on the echo signal, the phase signal is obtained.

[0008] Step 2: A Grey Wolf Optimization (GWO) algorithm that combines the maximum mutual information coefficient (MIC) and the minimum reconstruction error (MRE) as the fitness function is used to adaptively optimize the parameter penalty factor α and the number of modes k of the variational mode decomposition algorithm (VMD) parameters, determine the optimal α and k values, and then adaptively separate the vital sign signals.

[0009] Step 3: Use the improved variational mode decomposition (VMD) algorithm to decompose the original signal into k modal components u k (t), through step-by-step optimization, different modal components (IMF, Intrinsic Mode Function) are finally obtained, thereby effectively separating breathing, heartbeat signals and other noise components. Then, the modal selection criterion based on correlation coefficient is used to select the optimal heartbeat signal component from the decomposed IMF. Specifically, since the normal number of heartbeats in one minute is 48 to 120 times, the corresponding frequency range is 0.8-2.0Hz. Therefore, for the IMF decomposed by VMD, the energy x corresponding to the heartbeat frequency range (0.8-2.0Hz) in each IMF is first calculated, and then the correlation coefficient between each IMF energy and x is evaluated. Finally, the IMF with the highest correlation coefficient with x is selected as the optimal modal component to extract the required heartbeat signal.

[0010] Step 4: Extract the heartbeat waveform from the heartbeat signal, and obtain the time point corresponding to each peak by analyzing the peak of the heartbeat waveform. Then, the IBI sequence is obtained by calculating the difference between adjacent time points, thereby further obtaining the HRV time-frequency domain features. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] Figure 1 It is a flow chart of the present invention;

[0012] Figure 2 This is the phase signal diagram after preprocessing.

[0013] Figure 3 Time domain and frequency domain diagrams of the modal components of the signal decomposed by VMD.

[0014] Figure 4 The time domain and frequency domain diagrams of the heartbeat waveform extracted based on the modal selection criteria.

[0015] Figure 5Extract time domain plots for peak detection based IBI sequences.

[0016] Figure 6 This is the time domain result diagram of HRV. Specific implementation plan

[0017] The present invention will be further described below in conjunction with the accompanying drawings:

[0018] like Figure 1 A non-contact heart rate variability detection method based on improved variational mode decomposition and mode selection criteria is shown, which specifically includes the following steps:

[0019] Step 1: Use frequency modulated continuous wave (FMCW) radar to collect human vital signs radar signals. Use the distance-FFT method to obtain human distance information and extract the corresponding echo signal. After performing arc tangent demodulation and a series of pre-processing noise reduction reconstruction on the echo signal, the phase signal is obtained. The signal x(t) is as follows: Figure 2 shown.

[0020] Step 2: Use a GWO optimization algorithm that combines MIC and MRE as the fitness function to automatically optimize the VMD parameter penalty factor α and the number of modes k, determine the optimal α and k values, and then adaptively separate the vital sign signals. MIC originates from the concept of mutual information, which is defined as formula (1):

[0021]

[0022] Among them, X and Y are two random variables, p(X,Y) is the joint probability of X and Y, and p(X) and p(Y) are the marginal distributions of X and Y respectively. MIC divides the values ​​of the variables into grids and calculates the normalized maximum mutual information under a fixed grid:

[0023]

[0024] Where I(X,Y|dx,dy) is the mutual information of X and Y when divided into a grid of dx and dy, and B represents the upper limit of the grid size. log2min(dx,dy) normalizes the MIC value to ensure that the MIC value is between [0,1].

[0025] In order to alleviate modal aliasing, the MIC between modal component IMFs should be minimized. At the same time, in order to retain effective information, the MIC between modal component IMFs and the input signal x(t) should be maximized. On this basis, the minimum reconstruction error MRE indicator is added to measure the overall decomposition effect. Therefore, the fitness function is as follows:

[0026]

[0027] in, λ1 and λ2 are weights.

[0028] Next, the GWO optimization algorithm is used to adaptively optimize the VMD parameter penalty factor α and the number of modes k to determine the optimal α and k values. First, initialize the individual positions of a group of gray wolves, each position represents a solution, and the initial population matrix is:

[0029]

[0030] Among them, the size of the gray wolf population is N, and the dimension of the solution is 2, namely α and k. According to the fitness function fobj(k,α), the fitness value corresponding to the position of each individual is calculated, and the positions of α, β, and δ wolves are determined, f(X i )=fobj(k,α) is the fitness function, then: X α represents the position with the best fitness, X β Indicates the position with suboptimal fitness, X δ Indicates the third position in fitness. Update the positions of the remaining gray wolves according to the positions of α, β, and δ wolves. First, calculate the distance vector:

[0031]

[0032] Among them, C1, C2, C3 are random vectors in the range [0,2]: C i =2·r i , r i The modulus is a random number between [0,1]. Then update the position vector:

[0033]

[0034] Among them, A1, A2, A3 are coefficient vectors in the range [-a, a], and a is a parameter that decreases over time: A i =2·a·r i -a, a is the convergence factor that decreases linearly from 2 to 0 as the number of iterations increases, and the modulus of r is a random number between [0,1]. Finally, update the position of the gray wolf:

[0035]

[0036] Through the above position update process, the position of the wolf pack is continuously iterated and updated. The final position of the α wolf is the final solution to the optimization problem. According to its fitness function f(X α ) obtains the optimal solution for the penalty factor α and the number of modes k.

[0037] Step 3: Use the variational mode decomposition (VMD) algorithm to decompose the original signal into k modal components u k (t), through step-by-step optimization, different modal component IMFs are finally obtained, thereby effectively separating the breathing, heartbeat signals and other noise components. First, for each modal component u k (t) Perform Hilbert transform:

[0038]

[0039] Modulate the spectrum of each mode function to the corresponding baseband:

[0040]

[0041] The bandwidth of each IMF is estimated by the L2 norm of the squared signal gradient in formula (9), forming a constrained variational problem:

[0042]

[0043] Among them, u k (t) are all modal components, ω k is the center frequency corresponding to all modes.

[0044] In order to solve the above constrained optimization problem, the constrained variational problem is transformed into an unconstrained variational problem, and the augmented Lagrangian function is introduced by using the quadratic penalty term α and the Lagrangian multiplier λ(t).

[0045]

[0046] The alternating direction method of multipliers (ADMM) is used to solve the problem. The idea is to fix the other two variables and update one of them. and n+1 Find the Lagrangian saddle point and get The optimal solution effectively separates breathing, heartbeat signals and other noise components. The final decomposition mode diagram is as follows Figure 3 shown.

[0047] Secondly, since the normal number of heartbeats in one minute is 48 to 120 times, the corresponding frequency range is 0.8-2.0 Hz. Therefore, for the IMF decomposed by VMD, the energy x corresponding to the heartbeat frequency range (0.8-2.0 Hz) in each IMF is first calculated, and then the correlation coefficient between each IMF energy and x is evaluated. The mode selection criterion based on the correlation coefficient is used to select the optimal heartbeat signal component from the decomposed IMF, as shown in formula (12):

[0048]

[0049] Among them, IMF is the total energy of each modal component, x is the energy corresponding to the heartbeat component in each IMF, and E represents the mean value of the corresponding signal.

[0050] Finally, the IMF with the highest correlation coefficient with x is selected as the optimal modal component to extract the required heartbeat signal. The best heartbeat signal is selected as shown in formula (13). The result is as follows: Figure 4 shown.

[0051] IMF best =IMF k When k = argmax k (ρ k ) (13)

[0052] Step 4: Extract the heartbeat waveform from the heartbeat signal, and obtain the time point corresponding to each peak by analyzing the peak of the heartbeat waveform. Then, calculate the difference between adjacent time points to obtain the IBI sequence. The result is as follows: Figure 5 As shown, the HRV time-frequency domain characteristics are further obtained.

[0053] (1)HRV time domain characteristics, the results are as follows Figure 6 As shown:

[0054] ①Mean: The mean of the IBI sequence.

[0055]

[0056] Among them, N IBI is the estimated total number of IBIs.

[0057] ② Overall standard deviation (SDRR): the standard deviation of all IBIs.

[0058]

[0059] ③RMSSD: The mean of the sum of continuous differences and then its square root is used to measure the continuous IBI changes.

[0060]

[0061] ④SDSD: Square the difference between the sum of consecutive differences and the mean, calculate the mean, and then calculate the square root.

[0062]

[0063] ⑤pNN50: The percentage of consecutive IBIs with a difference of more than 50ms.

[0064]

[0065] (2) HRV frequency domain characteristics: The frequency domain analysis method of HRV is also called power spectrum analysis. The IBI sequence is decomposed into different components through fast Fourier transform or autoregressive model, and then the spectrum of these component data is estimated to obtain the HRV frequency domain characteristics. This paper adopts the Welch method, an improved method of fast Fourier transform, to perform spectrum analysis on HRV. The purpose of using this method is to reduce the noise interference in the signal by processing the signal by segmented windowing and averaging the power spectrum of each segment, thereby enhancing the accuracy of HRV frequency domain characteristics. The obtained frequency domain features include: low frequency (LF, 0.04-0.15Hz), high frequency (HF, 0.15-0.4Hz) and the ratio of low frequency to high frequency (LF / HF).

Claims

1. A non-contact heart rate variability detection method based on improved variational mode decomposition and mode selection criteria, characterized in that The following steps are involved: Step 1: Use frequency modulated continuous wave (FMCW) radar to collect human vital signs radar signals. Use the distance-FFT method to obtain human distance information and extract the corresponding echo signal. After performing arc tangent demodulation and a series of pre-processing noise reduction reconstruction on the echo signal, the phase signal x(t) is obtained. Step 2: A Grey Wolf Optimization (GWO) optimization algorithm that combines the Maximum Information Coefficient (MIC) and the Minimum Reconstruction Error (MRE) as the fitness function is used to automatically optimize the penalty factor α and the number of modes k of the variational mode decomposition (VMD) algorithm parameters, determine the optimal α and k values, and then adaptively separate the vital sign signals. MIC originates from the concept of mutual information, which is defined as formula (1): Among them, X and Y are two random variables, p(X,Y) is the joint probability of X and Y, and p(X) and p(Y) are the marginal distributions of X and Y respectively. MIC divides the values ​​of the variables into grids and calculates the normalized maximum mutual information under a fixed grid: Where I(X,Y|dx,dy) is the mutual information of X and Y when divided into a grid of dx and dy, and B represents the upper limit of the grid size. log2min(dx,dy) normalizes the MIC value to ensure that the MIC value is between [0,1]. In order to alleviate modal aliasing, the MIC between modal components (Intrinsic Mode Function, IMF) should be minimized. At the same time, in order to retain effective information, the MIC between the modal component IMF and the input signal x(t) should be maximized. On this basis, the minimum reconstruction error MRE indicator is added to measure the overall decomposition effect. Therefore, the fitness function is as follows: in, λ1 and λ2 are weights. Next, the GWO optimization algorithm is used to adaptively optimize the VMD parameter penalty factor α and the number of modes k to determine the optimal α and k values. First, initialize the individual positions of a group of gray wolves, each position represents a solution, and the initial population matrix is: Among them, the size of the gray wolf population is N, and the dimension of the solution is 2, namely α and k. According to the fitness function fobj(k,α), the fitness value corresponding to the position of each individual is calculated, and the positions of α, β, and δ wolves are determined, f(X i )=fobj(k,α) is the fitness function, then: X α represents the position with the best fitness, X β Indicates the position with suboptimal fitness, X δ Indicates the third position in fitness. Update the positions of the remaining gray wolves according to the positions of α, β, and δ wolves. First, calculate the distance vector: Among them, C1, C2, C3 are random vectors in the range [0,2]: C i =2·r i , r i The modulus is a random number between [0,1]. Then update the position vector: Among them, A1, A2, A3 are coefficient vectors in the range [-a, a], and a is a parameter that decreases over time: A i =2·a·r i -a, a is the convergence factor that decreases linearly from 2 to 0 as the number of iterations increases, and the modulus of r is a random number between [0,1]. Finally, update the position of the gray wolf: Through the above position update process, the position of the wolf pack is continuously iterated and updated. The final position of the α wolf is the final solution to the optimization problem. According to its fitness function f(X α ) obtains the optimal solution for the penalty factor α and the number of modes k. Step 3: Use the variational mode decomposition (VMD) algorithm to decompose the original signal into k modal components u k (t), through step-by-step optimization, different modal component IMFs are finally obtained, thereby effectively separating the breathing, heartbeat signals and other noise components. First, for each modal component u k (t) Perform Hilbert transform: Modulate the spectrum of each mode function to the corresponding baseband: The bandwidth of each IMF is estimated by the L2 norm of the squared signal gradient in formula (9), forming a constrained variational problem: Among them, u k (t) are all modal components, ω k is the center frequency corresponding to all modes. In order to solve the above constrained optimization problem, the constrained variational problem is transformed into an unconstrained variational problem, and the augmented Lagrangian function is introduced by using the quadratic penalty term α and the Lagrangian multiplier λ(t). The alternating direction method of multipliers (ADMM) is used to solve the problem. The idea is to fix the other two variables and update one of them. k n+1 ,ω k n+1 and n+1 Find the Lagrangian saddle point to obtain u k n+1 (ω),ω k n+1 (ω) The optimal solution effectively separates breathing, heartbeat signals and other noise components. Secondly, since the normal number of heartbeats in one minute is 48 to 120 times, the corresponding frequency range is 0.8-2.0 Hz. Therefore, for the IMF decomposed by VMD, the energy x corresponding to the heartbeat frequency range (0.8-2.0 Hz) in each IMF is first calculated, and then the correlation coefficient between each IMF energy and x is evaluated. The mode selection criterion based on the correlation coefficient is used to select the optimal heartbeat signal component from the decomposed IMF, as shown in formula (12): Among them, IMF is the total energy of each modal component, x is the energy corresponding to the heartbeat component in each IMF, and E represents the mean value of the corresponding signal. Finally, the IMF with the highest correlation coefficient with x is selected as the optimal modal component to extract the required heartbeat signal, as shown in formula (13). IMF best = IMF k When k = argmax k (ρ k ) (13) Step 4: Extract the heartbeat waveform from the heartbeat signal, and obtain the time point corresponding to each peak by analyzing the peak of the heartbeat waveform. Then, the IBI sequence is obtained by calculating the difference between adjacent time points, thereby further obtaining the HRV time-frequency domain features.

Images