A non-contact heart rate detection method based on weak head vibration

Abstract

The application discloses a non-contact heart rate detection method based on weak head vibration, which comprises the following steps: 1, extracting the motion track of the feature points in the region of interest, including face detection, selection of the region of interest, tracking the feature points in the region of interest to obtain the motion track; 2, reconstructing the motion track of the feature points into the motion track in the front-back direction, including calculating the distance between adjacent feature points, taking the time sequence of the distance as the motion track in the front-back direction, and removing unstable feature points; 3, motion detection and motion noise removal, including signal preprocessing, motion detection, motion noise removal and signal dimension reduction; 4, heart rate estimation, including fusing signals and performing spectrum analysis to obtain the heart rate value. The application tracks the weak head vibration caused by heartbeat from a video, reconstructs the tracking track, removes motion noise, thereby improving the accuracy of detecting the heart rate through the weak head vibration, and provides a new solution for non-contact heart rate detection.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of video-based non-contact physiological parameter detection, and particularly relates to a non-contact video heart rate detection method based on weak head vibration. BACKGROUND

[0002] With the continuous development of biomedical engineering field and the progress of technology, remote health monitoring has become an important part of modern medical services. Its emergence not only brings more convenient and efficient medical services for the majority of patients, but also provides more perfect and accurate medical security for medical institutions. In the field of remote health monitoring, heart rate is one of the most important and commonly used indicators. Because the heart rate of a person can not only reflect the health status of the body, but also reflect the emotional state, stress and other psychological states of a person. Traditionally, heart rate measurement needs to be achieved through body surface electrodes or blood pressure meters and other sensors. Although this method has high accuracy, it needs to be operated invasively, which brings inconvenience and discomfort to patients.

[0003] In recent years, with the rapid development of computer vision, image processing and other technologies, video-based non-contact heart rate measurement has become a popular research direction. At present, the most common two methods are remote photoplethysmography (rPPG) and imaging ballistocardiography (iBCG). Among them, the rPPG technology uses the micro changes of facial color to detect the blood volume changes caused by heartbeat, so as to obtain the heart rate estimation; while the iBCG technology measures the heart rate from the weak head vibration caused by the periodic movement of blood. Compared with rPPG, iBCG has the advantage of not using skin area, and iBCG is less sensitive to different types of skin color or environmental light changes.

[0004] However, iBCG heart rate measurement also faces some major difficulties. First, the source of iBCG signal is the weak vibration of head caused by heartbeat, which is very easy to be affected by motion artifacts caused by voluntary shaking of the subject's head, facial expressions and speaking, which will cause the increase of iBCG signal noise, thereby affecting the accurate measurement of heart rate. The current research mainly focuses on reducing the interference of motion noise, by removing low-quality video frames in advance, the motion noise contained in the iBCG signal can be reduced; through filtering, discrete cosine transform, independent component analysis and other technologies, the motion noise in the iBCG signal can be separated. Second, the weak mechanical movement of the head caused by the heartbeat will show different characteristics when the subject is in different postures. For the common application scenarios of non-contact heart rate monitoring, such as office, home or interrogation environment, the subject is mostly in a sitting posture, and the vibration amplitude of the human body in the front-back direction caused by heartbeat is larger. At present, the depth camera that can directly collect high-precision front-back direction motion information is mostly expensive and difficult to apply in general scenarios. How to extract reliable iBCG signal and remove motion noise is the main challenge of iBCG heart rate monitoring. SUMMARY

[0005] The present application is to solve the above-mentioned problems of the prior art, and proposes a non-contact heart rate detection method based on weak head vibration, in order to improve the accuracy and robustness of the heart rate detection affected by motion artifacts when detecting heart rate through weak head vibration, thereby providing a feasible scheme for the practical application of non-contact heart rate detection based on weak head vibration.

[0006] The present application is to solve the above-mentioned problems of the prior art, and proposes a non-contact heart rate detection method based on weak head vibration, in order to improve the accuracy and robustness of the heart rate detection affected by motion artifacts when detecting heart rate through weak head vibration, thereby providing a feasible scheme for the practical application of non-contact heart rate detection based on weak head vibration.

[0007] The present application is to solve the above-mentioned problems of the prior art, and proposes a non-contact heart rate detection method based on weak head vibration, in order to improve the accuracy and robustness of the heart rate detection affected by motion artifacts when detecting heart rate through weak head vibration, thereby providing a feasible scheme for the practical application of non-contact heart rate detection based on weak head vibration.

[0008] Step one: extracting feature point motion trajectory of the face region of interest from the T-frame video image of the subject;

[0009] Step 1.1, positioning of the region of interest;

[0010] The face rectangle region containing the head in the first frame of video image is obtained by using the face detection algorithm, and K regions are selected as the region of interest from the face rectangle region according to the proportional relationship of the face;

[0011] Step 1.2, feature point detection and tracking;

[0012] The N [k] feature points of the kth region of interest in the first frame of video image are determined by using the feature point detection method; the N [k]tracking the feature points on the T-th frame of the video images, thereby obtaining the horizontal motion trajectory X [k] and the vertical motion trajectory Y [k] of the N [k] feature points of the k-th region of interest on the T-th frame of the video images [k] [k] [k] [k] [k] [k] [k] [k] , wherein x [k] (t) and y [k] (t) represent the horizontal coordinate and the vertical coordinate of the n-th feature point of the k-th region of interest on the t-th frame of the video images, respectively, and represents the horizontal coordinate of the n-th feature point of the k-th region of interest on the t-th frame of the video images, represents the vertical coordinate of the n-th feature point of the k-th region of interest on the t-th frame of the video images.

[0013] Step two, reconstruct the motion trajectories in the vertical and horizontal directions into the motion trajectories in the front-back direction, and remove the unstable feature points.

[0014] Step 2.1, calculate the distance d n，j [k] (t) between the n-th feature point of the k-th region of interest on the t-th frame of the video images and other j-th feature points by using the Euclidean distance, thereby obtaining the time sequence of the distance between the n-th feature point and all other feature points, and taking it as the front-back direction trajectory of the n-th feature point on the t-th frame of the video images.

[0015] Select the P feature points closest to the n-th feature point of the k-th region of interest on the first frame of the video images, and calculate the average value of the distance between the n-th feature point and the closest P feature points on the t-th frame of the video images , thereby obtaining the average time sequence of the distance between the n-th feature point and the closest P feature points on the t-th frame of the video images, and taking it as the front-back direction motion trajectory z [k] (t) of the N [k] feature points on the t-th frame of the video images; and further obtaining the front-back direction motion trajectory Z [k] of the N [k] feature points of the k-th region of interest on the T-th frame of the video images.​​​​​​​[k] =[z [k] (1), z [k] (2), ..., z [k] (t), ..., z [k] (T)];

[0016] Step 2.2: Remove unstable feature points within each region of interest;

[0017] Step 2.2.1: Calculate the ratio of the standard deviation to the mean of the forward and backward trajectories of all feature points in the k-th region of interest on the T-frame video image. in, This represents the ratio of the k-th region of interest in a T-frame video image. std(z n [k] ) represents the standard deviation of the nth feature point in the kth region of interest across T frames of video image, mean(z) n [k] ) represents the average value of the nth feature point in the kth region of interest across T frames of video image; z n [k] This represents the forward and backward motion trajectory of the nth feature point in the kth region of interest on the T-frame video image.

[0018] Step 2.2.2: Use a statistical histogram to plot r [k] The elements in the dataset are divided into H intervals; after counting the number of elements in each interval, the elements are sorted in descending order and then summed until the sum of the element counts is greater than N. [k] When ρ% is reached, the accumulation stops, and the average of the midpoints of all intervals involved in the accumulation is used as the reference value for the k-th region of interest. Where ρ represents the threshold;

[0019] Step 2.2.3, Definition The two deviation limits include: the lower limit and upper limit Remove r [k] Median value exceeds The feature points are obtained by identifying stable feature points within the k-th region of interest; where σ min σ represents the range of downward deviation from the reference value. max This indicates the range of upward deviation from the reference value;

[0020] Step 3: Motion noise detection and removal;

[0021] Step 3.1, Signal Preprocessing:

[0022] Set the detrending parameter as λ, and detrend the front-back direction trajectory Z [k] , horizontal direction trajectory X [k] , and then filter the detrended trajectory signal with a Butterworth filter to obtain the preprocessed front-back direction signal Z' [k] and horizontal direction signal X' [k] ;

[0023] Step 3.2, motion detection:

[0024] Calculate the average value δ [k] of the standard deviation of the vertical motion trajectory Y [k] of all feature points in the kth region of interest, and compare it with the pre-set threshold μ. If δ [k] > μ, it is determined that the kth region of interest contains motion noise, and step 3.3 is executed. Otherwise, it is indicated that the kth region of interest does not contain motion noise, and the preprocessed front-back direction signal Z' [k] is taken as the pure signal set A [k] in the kth region of interest, and step 3.4 is executed.

[0025] Step 3.3, motion noise removal:

[0026] The preprocessed front-back direction signal Z' [k] in the kth region of interest is filtered one by one with a normalized least mean square adaptive filter to obtain the pure signal set A [k] in the kth region of interest.

[0027] Step 3.4, applying principal component analysis to the pure signal set A [k] in the kth region of interest to obtain the eigenvalues of the principal components, and arranging them in descending order, and retaining the first J principal components as the signal set B [k] of the kth region of interest containing heart rate information; thereby obtaining the signal set B = {B [1] ,..., B [k] ,..., B [K]} of the K regions of interest.

[0028] Step four, fuse the signals of the K regions of interest to obtain the target typical variable and calculate the heart rate value.

[0029] Step 4.1, apply multivariate canonical correlation analysis to B to obtain K sets of canonical correlation variables arranged in descending order of correlation coefficients wherein, represents the kth set of canonical correlation variables, represents the jth canonical correlation variable in the kth set of canonical correlation variables.

[0030] Step 4.2, spectrum analysis is performed on the KxJ typical correlation variables, and a typical correlation variable with the highest normalized power peak value is selected as a target signal, a peak frequency f of the target signal is taken as a final heart rate frequency, and thus an estimated value of the heart rate is fxtimes, wherein times represents a unit time.

[0031] The electronic device comprises a memory and a processor, and the memory is used to store a program supporting the processor to execute the non-contact heart rate detection method, and the processor is configured to execute the program stored in the memory.

[0032] The computer readable storage medium stores a computer program, and when the computer program is executed by a processor, the steps of the non-contact heart rate detection method are executed.

[0033] Compared with the prior art, the beneficial effects of the present application are embodied in that:

[0034] 1. The current iBCG method directly uses the vertical trajectory of the head to extract the iBCG signal, that is, the displacement of the head-to-foot direction of the human body. When the subject is in different postures, the amplitudes of the iBCG motion in the three directions in space are not the same. For most remote heart rate monitoring application scenarios, such as office, home or interrogation environment, the subject is mostly in a sitting posture. In the sitting state, studies have shown that the amplitude of the BCG motion in the front-back direction of the human body is greater than that in the head-to-foot direction. The present application proposes a new method according to the perspective law, extracts the iBCG signal in the front-back direction in the RGB video through the change of the distance between the feature points, improves the quality of the iBCG signal, and thus obtains a more accurate heart rate estimate.

[0035] 2. Due to the relatively weak motion caused by periodic heart activity, it is easy to be dominated by random head movement or facial expression. When the head movement amplitude of the target is slightly large, it is difficult to extract the true heart rate from the iBCG signal affected by the motion noise. Since the motion noise has the characteristics of strong randomness and non-stationarity, the present application adopts an adaptive filter that adjusts the filter coefficients according to the normalized least mean square algorithm, uses the horizontal direction motion trajectory as a noise reference to adjust the filter coefficients, tracks the motion noise signal, and removes it from the BCG signal, thereby greatly improving the quality of the iBCG signal.

[0036] 3、The iBCG signal is fully considered to be derived from the head vibration caused by heartbeat, and is globally correlated with the head. In order to fully apply the information of the head correlation of the iBCG signal, typical correlation analysis is adopted to extract common components from the region of interest, improve the heart rate component of the iBCG signal, further improve the accuracy and robustness of the video-based non-contact heart rate detection, and provide reliable guarantee for the application of the technology. BRIEF DESCRIPTION OF DRAWINGS

[0037] Figure 1 The flowchart of the present application.

[0038] Figure 2 The region of interest division schematic diagram of the present application.

[0039] Figure 3 The signal quality comparison schematic diagram of the reconstructed signal proposed by the present application. DETAILED DESCRIPTION

[0040] In the embodiment, a non-contact heart rate detection method (Z-MAR) based on weak head vibration is a method for positioning the region of interest of a video image, determining and tracking the feature points in the region of interest to obtain a motion trajectory containing heart rate information, and reconstructing the motion trajectory into a front-back direction motion trajectory. Then, adaptive filtering and typical correlation analysis are applied to remove motion noise in the signal and extract the heart rate component to obtain the final heart rate estimation. Specifically, as shown in the figure, the method comprises the following steps: Figure 1

[0041] Step 1: Extract the feature point motion trajectory of the face region of interest from the T-frame video image of the subject;

[0042] Step 1.1, positioning of the region of interest;

[0043] The Viola-Jones face detection algorithm is used to obtain the face rectangular region containing the head in the first frame of video image, and K regions are selected as the region of interest from the face rectangular region according to the proportional relationship of the face;

[0044] In the specific embodiment, the left top vertex of the face rectangular region and the width and height are denoted as [x, y, w, h], K=2, and the forehead region [x+0.2Sw, y, 0.5w, 0.2h] and the nose region [x+0.25w, y+0.5h, 0.5w, 0.2h] are selected as the region of interest according to the proportional relationship, as shown in the figure. Figure 2

[0045] Step 1.2, feature point detection and tracking;

[0046] ​​N feature points of the kth region of interest in the first frame of video image are determined by using a good-feature-to-track (GFT) feature point detection method. [k] The N feature points of the kth region of interest are tracked on T frames of video image by using a Kande Lucas-Tomasi (KLT) feature point tracking method, so as to obtain horizontal direction motion trajectories X [k] and vertical direction motion trajectories Y [k] of the N feature points of the kth region of interest in the T frames of video image, wherein X [k] = [x [k] (1), x [k] (2),..., x [k] (t),..., x [k] (T)] and Y [k] = [y [k] (1), y [k] (2),..., y [k] (t),..., y [k] (T)], x [k] (t) and y [k] (t) represent horizontal coordinates and vertical coordinates of the N feature points of the kth region of interest in the tth frame of video image, respectively, and x (n, t) represents the horizontal coordinate of the nth feature point of the kth region of interest in the tth frame of video image, y

[0047] (n, t) represents the vertical coordinate of the nth feature point of the kth region of interest in the tth frame of video image.

[0048] Step two, the motion trajectories in vertical and horizontal directions are reconstructed into motion trajectories in front-back directions, and unstable feature points are removed;

[0049] Step 2.1, the distance d n，j [k] (n, t) between the nth feature point of the kth region of interest in the tth frame of video image and other jth feature point is calculated by using Euclidean distance, so as to obtain a time sequence of the distance between the nth feature point and all other feature points, and the time sequence is taken as the front-back direction trajectory of the nth feature point of the kth region of interest in the tth frame of video image; P feature points closest to the nth feature point of the kth region of interest in the first frame of video image are selected, and the average value of the distance between the nth feature point of the kth region of interest in the tth frame of video image and the closest P feature points is calculated.This yields the average time series of the distances between the nth feature point of the kth region of interest on the tth frame of the video image and its nearest P feature points, and serves as the N-th feature point on the tth frame of the video image. [k] The trajectory of the feature point in the forward and backward directions z [k] (t); and then obtain the N of the k-th region of interest on the T-frame video image. [k] The forward and backward movement trajectory Z of each feature point [k] =[z [k] (1), z [k] (2), ..., z [k] (t), ..., z [k] (T)]; In this embodiment, P is 2.

[0050] Step 2.2: Remove unstable feature points within each region of interest;

[0051] Step 2.2.1: Calculate the ratio of the standard deviation to the mean of the forward and backward trajectories of all feature points in the k-th region of interest on the T-frame video image. in, This represents the ratio of the k-th region of interest in a T-frame video image. std(z n [k] ) represents the standard deviation of the nth feature point in the kth region of interest across T frames of video image, mean(z) n [k] ) represents the average value of the nth feature point in the kth region of interest across T frames of video image; z n [k] This represents the forward and backward motion trajectory of the nth feature point in the kth region of interest on the T-frame video image.

[0052] Step 2.2.2: Use a statistical histogram to plot r [k] The elements in the dataset are divided into H intervals; after counting the number of elements in each interval, the elements are sorted in descending order and then summed until the sum of the element counts is greater than N. [k] When ρ% is reached, the accumulation stops, and the average of the midpoints of all intervals involved in the accumulation is used as the reference value for the k-th region of interest. Where ρ represents the threshold;

[0053] Step 2.2.3, Definition The two deviation limits include: the lower limit and upper limit Remove r [k] Median value exceeds The feature points are obtained by identifying stable feature points within the k-th region of interest; where σ minσ represents the range deviating from the reference value downward max σ represents the range deviating from the reference value upward.

[0054] In this embodiment, the value of H is 10, the value of p is 50, and the value of σ min and σ max are respectively 0.25 and 5.

[0055] Step three, motion noise detection and motion noise removal;

[0056] Step 3.1, signal preprocessing:

[0057] The detrending parameter is set as λ, and the front and rear direction trajectories Z [k] and the horizontal direction trajectory X [k] of the k-th region of interest are subjected to detrending processing, and then the Butterworth filter is used to filter the trajectory signals after detrending, to obtain the preprocessed front and rear direction signal Z' [k] and the horizontal direction signal X' [k] .

[0058] In a specific embodiment, the value of λ is 15, and the cutoff frequency of the Butterworth filter is [0.75hz, 3hz] to filter out components other than the heart rate frequency band.

[0059] Step 3.2, motion detection:

[0060] The average value δ [k] of the standard deviation of the vertical motion trajectory Y [k] of all feature points in the k-th region of interest is calculated, and compared with the pre-set threshold μ, if δ [k] > μ, it is determined that the k-th region of interest contains motion noise, and step 3.3 is executed; otherwise, it is indicated that the k-th region of interest does not contain motion noise, and the preprocessed front and rear direction signal Z' [k] is taken as the pure signal set A [k] in the k-th region of interest, and step 3.4 is executed; in a specific embodiment, the value of μ is 2.5.

[0061] Step 3.3, motion noise removal:

[0062] The normalized least mean square adaptive filter is used to filter the preprocessed front and rear direction signal Z' [k] in the k-th region of interest one by one, to obtain the pure signal set A [k] in the k-th region of interest.

[0063] For a p-order adaptive filter, at the t-th moment, the filter coefficient is m t . The sliding window is used to calculate the average value of the filter coefficient in the last n time points, and the average value is taken as the filter coefficient at the t-th moment. The segment with the intercept length p is taken as the input signal of the adaptive filter at the tth time The filter output y(t) at the tth time is:

[0064]

[0065] The The error e(t) between the adaptive filter output y(t) and the desired signal at the tth time is:

[0066]

[0067] The filter coefficient m t is updated by formula (5):

[0068]

[0069] In a specific embodiment, p is set to 50, μ is set to 0.003, and ∈ is a number close to 0, set to 10 -6 .

[0070] Solving formula (3), formula (4), and formula (5) performs adaptive filtering. The error e(t) between the input signal and the desired signal in adaptive filtering is taken as the pure signal after removing motion noise. After adaptive filtering is performed on all signals in the kth region of interest, all pure signals constitute a new signal set A [k] .

[0071] Step 3.4, applying principal component analysis to the pure signal set A [k] in the kth region of interest, the characteristic values of the principal component components are obtained after processing, and the first J principal components are retained as the signal set B of the kth region of interest containing heart rate information [k] ; thereby obtaining the signal set B = {B [1] ,..., B [k] ,..., B [K]} of the K regions of interest; in a specific embodiment, J is set to 5, and the first five principal components contain most of the information related to heart rate.

[0072] Step four, fusing the signals of the K regions of interest to obtain a target canonical variable and calculating a heart rate value;

[0073] Step 4.1, applying multivariate canonical correlation analysis to B to obtain K sets of canonical correlation variable sets arranged in descending order of correlation coefficients wherein, Bk represents the kth set of canonical correlation variable sets, represents the jth canonical correlation variable in the kth group of canonical correlation variables.

[0074] Step 4.2, spectral analysis is performed on the KxJ canonical correlation variables, and a canonical correlation variable with the highest normalized power peak value is selected as a target signal, a peak frequency f of the target signal is taken as a final heart rate frequency, so that an estimated value of the heart rate is fxtimes, wherein time represents a unit time.

[0075] In this embodiment, an electronic device includes a memory for storing a program supporting a processor to execute the above method, and the processor is configured to execute the program stored in the memory.

[0076] In this embodiment, a computer readable storage medium has a computer program stored thereon, and the computer program is executed by a processor to perform the steps of the above method.

[0077] To verify the effectiveness of the method, the performance of the method (Z-axismotion artifact removal, Z-MAR) is verified on the public database UBFC-RPPG. The heart rate detection performance is evaluated by using conventional evaluation indexes, specifically including mean absolute error HR mae (mean absolute error, MAE), standard deviation HR sd (standard deviation, sd), Pearson's correlation coefficient (r), and root mean square error HR rmse (root mean square error, RMSE).

[0078] Table 1 UBFC database experimental results

[0079] HR mae (bpm) HR sd (bpm) r HR rmse (bpm) iBCG-Ori 24.36 24.28 0.15 34.39 iBCG-TLF 23.25 23.33 0.11 32.94 RiBCG 5.43 8.6 0.85 10.22 Y-MAR 4.19 9.35 0.84 10.25 Z-MAR 2.94 5.95 0.93 6.64

[0080] Table 2 UBFC database motion window experimental results

[0081] HR mae (bpm) HR sd (bpm) r HR rmse (bpm) iBCG-Ori 36.45 23.68 0.09 43.46 iBCG-TLF 30.66 24.25 0.27 39.09 RiBCG 8.15 10.94 0.77 13.64 Y-MAR 6.74 11.64 0.79 13.45 Z-MAR 5.18 7.16 0.90 8.83

[0082] Table 1 presents the results of the proposed method Z-MAR and other methods on the UBFC-RPPG database, including the original iBCG method (iBCG original, iBCG-Ori), the iBCG method using a three-layer filter (iBCG three layer filter, iBCG-TLF), and the motion-robust iBCG method (Robust iBCG, RiBCG). To verify the effectiveness of the proposed method, the ablation experiment group Y-MAR is added, which uses the commonly used Y-axis trajectory to replace the Z-axis trajectory in the proposed Z-MAR as the iBCG source signal, thereby illustrating the contribution of the adaptive filter and the Z-axis trajectory applied in the present application to the accuracy of heart rate measurement. The results show that the proposed method Z-MAR performs best compared with other methods, achieving HR mae of 2.94 bpm, HR rmse of 6.64 bpm, and a correlation coefficient r of 0.93. By comparing Y-MAR with other methods, it can be seen that the adaptive filter used in the present application can achieve better performance in removing motion artifacts, even when the same Y-axis trajectory as other methods is used as the source signal, HR mae of 2.94 bpm is still reduced by 1.24 bpm compared with RiBCG. The comparison between the Z-MAR and Y-MAR methods shows that the reconstructed Z-axis trajectory can further reduce HR mae by 1.25 bpm compared with the commonly used Y-axis trajectory. This finding highlights the superior quality of the reconstructed anterior-posterior iBCG signal compared with the commonly used vertical iBCG signal, and proves that the proposed method can calculate reliable anterior-posterior trajectories for heart rate estimation. To further prove the performance of the proposed method in removing motion artifacts, a performance comparison is made when motion is detected.

[0083] Table 2 presents the results of the proposed method and other methods on all subjects in the UBFC-RPPG database who were detected to have significant motion. All methods show a performance decline, which highlights the sensitivity of the iBCG signal to motion artifacts and emphasizes the importance of motion artifact removal techniques for enhancing the motion robustness of the iBCG method. The proposed method still achieves the best results under the influence of motion, with HR mae of 5.18 bpm, HR sd of 7.16, a correlation coefficient r of 0.9, and HR rmse of 8.83 bpm. Compared with other methods, the proposed method has the smallest decline in accuracy under the influence of motion artifacts, with a Pearson correlation coefficient r of only 0.03, which further demonstrates the robustness of the proposed method.

[0084] Figure 3The X, Y, Z axis trajectories are compared with the reference PPG signal. The Z axis trajectory reconstructed by the proposed method is closer to the reference signal than the X, Y axis trajectories, which again proves that the Z axis trajectory is beneficial to iBCG heart rate measurement.

[0085] In summary, the proposed method can achieve the best performance compared with other methods. This study provides a new feasible scheme for the practical application of iBCG heart rate measurement.

Claims

1. A non-contact heart rate detection method based on weak head vibration, characterized in that, The method comprises the following steps: Step one: extracting feature point motion trajectory of the face region of interest from the T-frame video image of the subject; Step 1.1, positioning of the region of interest; A face rectangle region containing the head in the first frame of video image is obtained by using a face detection algorithm, and K regions are selected from the face rectangle region as the regions of interest according to the proportional relationship of the face; Step 1.2, feature point detection and tracking; The k-th region of interest in the first frame of the video image is determined using a feature point detection method. [k] N feature points; the feature point tracking method is used to track the N features of the k-th region of interest. [k] N feature points are tracked on T frames of video image to obtain the N features of the k-th region of interest in the T frames of video image. [k] The horizontal movement trajectory X of the feature points [k] =[x [k] (1), x [k] (2), ..., x [k] (t), ..., x [k] (T)] and vertical trajectory Y [k] =[y [k] (1), y [k] (2), ..., y [k] (t), ..., y [k] (T)], where x [k] (t) and y [k] (t) represent N values ​​of the k-th region of interest in the t-th frame of the video image. [k] The horizontal and vertical coordinates of each feature point, and This represents the horizontal coordinate of the nth feature point in the t-th frame of the video image of the k-th region of interest. This represents the vertical coordinate of the nth feature point of the kth region of interest in the tth frame of the video image; Step two, reconstructing the motion trajectory in the vertical direction and the horizontal direction into the motion trajectory in the front-back direction, and removing unstable feature points; Step 2.1, calculate the distance d between the n-th feature point of the k-th region of interest on the t-th frame of video image and other j-th feature point by using the Euclidean distance n，j [k] (t), thereby obtaining the time series of the distance between the n-th feature point and all other feature points, and taking it as the front-back direction trajectory of the n-th feature point of the k-th region of interest on the t-th frame of video image; In the k-th region of interest, select the P feature points closest to the n-th feature point on the first frame of the video image, and calculate the average distance between the n-th feature point of the k-th region of interest and the P nearest feature points on the t-th frame of the video image. This yields the average time series of the distances between the nth feature point of the kth region of interest on the tth frame of the video image and its nearest P feature points, and serves as the N-th feature point on the tth frame of the video image. [k] The trajectory of the feature point in the forward and backward directions z [k] (t); and then obtain the N of the k-th region of interest on the T-frame video image. [k] The forward and backward movement trajectory Z of each feature point [k] =[z [k] (1), z [k] (2), ..., z [k] (t), ..., z [k] (T)]; Step 2.2, removing unstable feature points in each region of interest; Step 2.2.1, calculating the ratio of the standard deviation and the mean of the front-back direction trajectory of all feature points of the kth region of interest on the T-frame video image wherein, represents the ratio of the kth region of interest on the T-frame video image, std(z n [k] ) represents the standard deviation of the nth feature point of the kth region of interest on the T-frame video image, mean(z n [k] ) represents the mean of the nth feature point of the kth region of interest on the T-frame video image; z n [k] represents the front-back direction motion trajectory of the nth feature point of the kth region of interest on the T-frame video image; Step 2.2.2, divide the elements in r [k] into H intervals using statistical histogram; and count the number of elements in each interval, then sort the number of elements in descending order and accumulate until the accumulated value of the number of elements is greater than ρ% of N [k] , stop the accumulation, and then take the average of the midpoints of all intervals involved in the accumulation as the reference value of the kth region of interest where ρ represents the threshold value; Step 2.2.3, definition of two deviation limits, including: lower limit and upper limit r [k] remove feature points whose values exceed , thereby obtaining stable feature points in the kth region of interest; wherein σ min represents the range of deviation from the reference value, and σ max represents the range of deviation from the reference value; Step three, motion noise detection and motion noise removal; Step 3.1, signal preprocessing: Set the detrending parameter as λ, and perform detrending on the front-back direction trajectory Z [k] , horizontal direction trajectory X [k] , and then filter the trajectory signal after detrending using a Butterworth filter to obtain the preprocessed front-back direction signal Z′ [k] and horizontal direction signal x′ [k] ; Step 3.2, motion detection: Calculate the vertical motion trajectory Y of all feature points in the kth region of interest [k] the average value of the standard deviation of the vertical motion trajectory Y of all feature points in the kth region of interest [k] and compare it with a pre-set threshold μ, if δ [k] > μ, it is determined that the kth region of interest contains motion noise, and step 3.3 is executed; otherwise, it is indicated that the kth region of interest does not contain motion noise, and the pre-processed front-back direction signal Z′ [k] is taken as the pure signal set A [k] in the kth region of interest, and step 3.4 is executed; Step 3.3, motion noise removal: The normalized least mean square error adaptive filter is used to filter the preprocessed front and back direction signals Z' in the kth region of interest [k] The filtering is performed one by one to obtain the pure signal set A in the kth region of interest [k] ; Step 3.

4. Apply principal component analysis to the pure signal set A in the kth region of interest [k] After processing, the eigenvalues of the principal component components are arranged in descending order, and the first J principal components are retained as the signal set B of the kth region of interest containing heart rate information [k] ; thereby obtaining the signal set B of the K regions of interest = {B [1] ,..., B [k] ,..., B [K]} ; Step four, fusing the signals of the K regions of interest to obtain target typical variables and calculating the heart rate value; Step 4.1, apply multivariate canonical correlation analysis to B to obtain K sets of canonical correlation variables in descending order of correlation coefficients wherein, denotes the set of canonical correlation variables of the kth group, denotes the jth canonical correlation variable in the set of canonical correlation variables of the kth group; Step 4.2, performing spectrum analysis on the K×J typical correlation variables, selecting the typical correlation variable with the highest normalized power peak value as the target signal, taking the peak frequency f of the target signal as the final heart rate frequency, and thus obtaining the estimated value of the heart rate as f×time, wherein time represents a unit of time.

2. An electronic device comprising a memory and a processor, characterized in that The memory is used to store a program supporting the processor to execute the non-contact heart rate detection method of claim 1, and the processor is configured to execute the program stored in the memory.

3. A computer-readable storage medium having stored thereon a computer program, characterized in that, The computer program is executed by the processor to perform the steps of the non-contact heart rate detection method of claim 1.

Images