A video heart rate detection method for removing irregular motion artifacts

Abstract

The application discloses a video heart rate detection method for removing irregular motion artifacts, which is characterized by the following steps: firstly, extracting the time sequence of green channel pixel mean value in the face region of interest in the video, mainly including the positioning, division and screening of the region of interest, tracking the region of interest and extracting the time sequence of green channel pixel mean value; secondly, extracting the heart rate signal, including signal preprocessing, using the method combining non-negative matrix factorization (NMF) and independent vector analysis (IVA) to perform blind source separation on the regular heart rate source signal, selecting the blood volume pulse signal; finally, calculating the heart rate, after performing Fourier transform on the blood volume pulse signal, the heart rate value is obtained through screening and calculation. The application extracts the heart rate signal from the video image, is used for heart rate detection under the condition of irregular motion artifacts, can improve the accuracy of heart rate detection under the condition of irregular motion, and thus provides a new method for non-contact heart rate detection research.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of video-based non-contact physiological parameter detection, and particularly relates to a video heart rate detection method for removing irregular motion artifacts by combining non-negative matrix factorization and independent vector analysis. BACKGROUND

[0002] Heart rate is an important physiological indicator in the diagnosis of cardiovascular diseases. A commonly used non-invasive detection method is a pulse oximeter. The pulse oximeter uses the PhotoPlethysmoGraph (PPG) technology to extract heart rate from the blood flow of the capillary tissue bed under the skin using an optical sensor. However, the optical sensor must be worn on the human skin, which often causes inconvenience in actual measurement. Blood flow changes can cause corresponding changes in facial skin color. Although such changes are not visible to the naked eye, they can be captured by a consumer-level camera, providing a low-cost and widely applicable non-contact detection method for heart rate detection. This video-based heart rate detection technology can be referred to as Remote PPG (RPPG).

[0003] Although the RPPG method can perform non-contact heart rate detection, the signal amplitude captured by the consumer-level camera is very weak, and is also disturbed by noise such as light and motion artifacts, causing the captured blood volume pulse (BVP) signal to be distorted. Therefore, improving noise robustness, especially motion noise robustness, is the key to ensuring the accuracy of the RPPG technology. In order to ensure the accuracy of the detected physiological parameters, traditional noise removal methods can be roughly divided into two categories: optical model-based methods. Because of the changes in the distance or angle between the face and the camera, motion denoising can be modeled as an optical model, and a color difference model and a plane orthogonal skin reflection model are proposed to solve the motion denoising problem. Blind signal separation-based methods, such as independent component analysis, independent vector analysis, and principal component analysis, separate mixed source signals by assuming the independence between source signals.

[0004] But in practical application, the subject often accompanies irregular head movement in the detection process, and the accuracy of heart rate signal detection will decrease with the movement, and the movement and light change will also affect the signal quality of the region of interest of the face. The method based on the optical model will cause the relative displacement of the position and angle between the region of interest of the face and the camera due to the movement, thereby affecting the accuracy of the result. The method based on principal component analysis and independent component analysis cannot meet the precondition of source signal independence and linear combination in some high amplitude and severe movement scenes, and also makes the heart rate detection result deviate. Independent vector analysis (IVA) solves the arrangement problem of source decomposition of independent component analysis, and is a better and more stable blind source separation method than independent component analysis and principal component analysis, but for regular signals such as heart rate signals, the target signal cannot be accurately extracted. SUMMARY

[0005] The present application overcomes the deficiencies of the prior art, and provides a video heart rate detection method for removing irregular motion artifacts, so as to improve the accuracy of heart rate detection in irregular motion scenes, and thus provides a new method for non-contact heart rate detection.

[0006] In order to achieve the above-mentioned application purposes, the present application adopts the following technical solutions:

[0007] The video heart rate detection method for removing irregular motion artifacts has the characteristics that the following steps are performed:

[0008] Step 1, a high-quality region of interest is selected from the face region of the face video image of the subject, and the pixel mean time series information of the green channel is extracted;

[0009] Step 1.1, the jth frame face region is identified from the jth frame face video image of the subject by using a face feature point detection method, and the jth frame face region is divided into N×N regions of interest, and a feature point tracking algorithm is used to identify and locate the face region in each frame of face video image, so as to obtain N×N regions of interest of each frame of face video image; j = 1, 2,..., J; J represents the total number of frames of face video image;

[0010] Step 1.2, the pixel mean of each region of interest in the jth frame face video image is calculated, so as to obtain the pixel mean of each region of interest in J frames of face video image, and to form a pixel mean time series of each region of interest; let the pixel mean time series of any single region of interest include: the red channel time series R = [R1, R2,..., RJ] of J frames of video image, the green channel time series G = [G1, G2,..., GJ] of J frames of video image, and the blue channel time series B = [B1, B2,..., BJ] of J frames of video image; j = 1, 2,..., J; J represents the total number of frames of face video image. j..., R J ] and the time series of the green channel G = [G1, G2,..., G j ..., G J ] in the J-frame video image; and the time series of the blue channel B = [B1, B2,..., B j ..., B J ] in the J-frame video image; wherein R j represents the pixel mean value of the red channel in the jth frame in a single region of interest, G j represents the pixel mean value of the green channel in the jth frame in a single region of interest, and B j represents the pixel mean value of the blue channel in the jth frame in a single region of interest.

[0011] Step 1.3, calculate the illumination intensity index, illumination change index and green signal signal-to-noise ratio index of a single region of interest; and sort the illumination intensity indexes of the N×N regions of interest in descending order, sort the illumination change indexes of the N×N regions of interest in ascending order, and sort the green signal signal-to-noise ratios of the N×N regions of interest in descending order; select the top Q regions of interest from the three sorted index sequences, and select the M regions of interest that exist in all the top Q regions of interest from the three sorted index sequences as the final selected high-quality regions of interest; M < Q < N×N.

[0012] Step two, extraction of heart rate signal under irregular motion scene.

[0013] Step 2.1, taking the green channel signals of the M high-quality regions of interest as the multi-channel signal C = [G [1] ,G [2] ,..., G [m] ,..., G [M] ] T ; wherein G [m] represents the green channel signal of the mth region of interest, i.e. the mth channel signal, and wherein, represents the pixel mean value of the jth frame of the green channel of the mth high-quality region of interest, m = 1, 2,..., M.

[0014] Step 2.2, set the detrending parameter as λ, and perform detrending processing on the multi-channel signal C to obtain the detrended multi-channel signal denoted as wherein, represents the mth channel signal after detrending.

[0015] Step 2.3, perform wavelet transform on the detrended multi-channel signal each channel signal in the filtered multi-channel signal is up-sampled to obtain an up-sampled transposed multi-channel signal matrix wherein, denotes the mth channel signal after filtering;

[0016] Step 2.4, each channel signal in the filtered multi-channel signal is up-sampled to obtain an up-sampled transposed multi-channel signal matrix wherein, denotes the mth channel signal after up-sampling transposition, and is a 1xJ matrix;

[0017] Step 2.5, the signal of the mth channel in the multi-channel signal matrix X is decomposed by using a non-negative matrix factorization method to obtain an mth feature matrix T [m] ∈R 1×L and an mth gain matrix V [m] ∈R L×J , and the mth feature matrix T [m] is decomposed by using a non-negative matrix factorization method to obtain an mth feature matrix T l [m] , the mth gain matrix V [m] is decomposed by using a non-negative matrix factorization method to obtain an mth feature matrix T l,j [m] , L is the number of feature vectors, and l = 0, 1, 2,..., L;

[0018] Step 2.6, a loss function Q of the improved independent vector analysis method is established by using formula (1):

[0019]

[0020] In formula (1), y j [m] denotes the jth source component estimation vector of the mth channel, W denotes a demixing matrix, and W = [w [1] ,w [2] ,...,w [m] ,...,w [M] ] H , wherein w [m] denotes a component vector of the demixing matrix of the mth channel; H denotes a conjugate transpose;

[0021] Step 2.7, the loss function Q is optimized and solved to obtain the source component estimation signal Y [m] of the mth channel;

[0022] Step three, calculation of heart rate;

[0023] ​The power spectral density distribution of the mth channel source component estimation signal Y [m] is calculated by using fast Fourier transform, and the source component estimation signal with the most frequency in the set range is determined as the blood volume pulse signal, so as to calculate the main frequency f max of the blood volume pulse signal for estimating the heart rate value.

[0024] The video heart rate detection method for removing irregular motion artifacts according to the present application is characterized in that the step 2.7 comprises:

[0025] Step 2.7.1, defining the current iteration number as inter, and initializing inter = 1;

[0026] Let W inter be a full one vector in the interth iteration; initialize the non-negative vector t [m] in the mth feature matrix T l in the interth iteration as t [m],inter l [m] ; initialize the non-negative vector v [m] in the mth gain matrix V l,j [m],inter in the interth iteration as v l,j [m] ;

[0027] Step 2.7.2, calculating the jth estimation variance r j [m],inter in the interth iteration by using formula (2):

[0028]

[0029] Step 2.7.3, calculating the intermediate variable Z [m],inter in the interth iteration by using formula (3):

[0030]

[0031] Step 2.7.4, calculating the component vector w' [m],inter of the unmixed matrix of the mth channel in the interth iteration by using formula (4):

[0032] w' [m],inter = (WZ [m],inter ) -1 e [m] (4)

[0033] In formula (4), e [m] represents the mth unit vector;

[0034] Step 2.7.5, compute the component vector w of the demixing matrix of the mth channel after the mth update at the interth iteration using formula (5) [m],inter :

[0035]

[0036] Step 2.7.6, compute the estimate vector y of the jth source component of the mth channel at the interth iteration using formula (6) j [m],inter :

[0037] y j [m],inter = (w [m],inter ) H X (6)

[0038] Step 2.7.7, compute the non-negative vector t of the lth eigenmatrix at the inter+1th iteration using formula (7) l [m ],inter+1 :

[0039]

[0040] Step 2.7.8, compute the non-negative vector v of the lth gain matrix at the inter+1th iteration using formula (8) lj [m],inter+1 :

[0041]

[0042] Step 2.7.9, after assigning inter+1 to inter, return to Step 2.8 to sequentially execute until inter>inter_max, thus obtaining the estimate vector y of the jth source component at the inter_maxth iteration j [m,inter_max and as the jth estimate vector y of the final source component of the mth channel j [m] ; thus obtaining the source component estimate signal Y of the mth channel [m] = [y1 [m] ,y2 [m] ,…,y j [m] ,…,y J [m] ].

[0043] The electronic device comprises a memory and a processor, and is characterized in that the memory is used for storing a program supporting the processor to execute the video heart rate detection method, and the processor is configured to execute the program stored in the memory.

[0044] The computer readable storage medium stores a computer program, and the computer program is characterized in that when the computer program is run by a processor, the steps of the video heart rate detection method are executed.

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

[0046] 1. The traditional blind source separation methods such as independent component analysis and principal component analysis cannot meet the prerequisite conditions of source signal independence and linear combination in the presence of large amplitude motion artifacts, which will cause large errors in heart rate measurement results. The method proposed in the present application has good blind source separation ability for regular heart rate signals, and even in the environment with large amplitude irregular motion artifacts, it can still realize the separation of quasi-periodic blood volume pulse signals, and has better motion robustness.

[0047] 2. Independent vector analysis has a wide application in the field of joint blind source separation, but when multiple data sets contain the same rigid motion artifact, the motion artifact will also be extracted as a common source component vector (SCV), which will affect the heart rate detection result. After the joint blind source separation method, it is also more complex to determine the target blood volume pulse signal, and usually needs spectral clustering or setting a certain criterion to achieve it. Although the non-negative matrix factorization (NMF) algorithm can decompose the characteristics of the heart rate signal, it is difficult to achieve the effect of blind source separation due to the lack of means to cluster the characteristics. The present application creatively combines non-negative matrix factorization with independent vector analysis method (hereinafter referred to as NMF-IVA algorithm), considers the non-negative matrix factorization cost function and the independent vector analysis cost function when determining the loss function of the demixing matrix, fully utilizes the source signal characteristic matrix obtained by non-negative matrix factorization, makes it sensitive to regular heart rate signals, and at the same time utilizes the source separation ability of independent vector analysis, finally realizes the separation of quasi-regular blood volume pulse signals under irregular motion artifact conditions, improves the motion robustness of video heart rate detection, and thus provides ideas for the practicalization process of this technology. BRIEF DESCRIPTION OF DRAWINGS

[0048] Figure 1 The flowchart of the present application;

[0049] Figure 2 The schematic diagram of the region of interest division and screening of the present application;

[0050] Figure 3 Comparison of the heart rate signal waveform obtained by the present application with the waveform obtained by IVA and the reference signal waveform;

[0051] Figure 4 Comparison chart of the heart rate signal spectrogram obtained by the present application and the spectrogram obtained by IVA. DETAILED DESCRIPTION

[0052] In this embodiment, the video heart rate detection method for removing irregular motion artifacts mainly uses non-negative matrix decomposition and independent vector analysis to detect the heart rate of the subject in the video under the motion scene. The J-frame video image is positioned, tracked and screened in the region of interest. The green channel time sequence in the color channel time sequence in the screened region of interest is selected for pretreatment. The noise is preliminarily removed through detrending, normalization and band-pass filtering. Then the non-negative matrix decomposition algorithm is used to decompose the characteristics of the pretreated multi-channel signal. The independent vector analysis is used to extract the heart rate signal. The heart rate pulse signal is screened to complete the calculation of the heart rate. As shown in the figure, the video heart rate detection method comprises the following specific steps: Figure 1

[0053] Step one, screen out high-quality sub-regions of interest from the face region of the face video image of the subject for extracting the pixel mean time sequence information of the green channel;

[0054] Step 1.1, identify the jth frame face region from the jth frame face video image of the subject by using the face feature point detection method, and divide the jth frame face region into N×N sub-regions of interest. The face region in each frame of face video image is identified and positioned by using the feature point tracking algorithm, so as to obtain N×N sub-regions of interest of each frame of face video image; j = 1, 2,..., J; J represents the total frame number of the face video image; in order to ensure the stability of the region of interest, in this embodiment, 60% of the detected face region without hair interference is selected as the divided face region. The detected region is divided into N×N regions as sub-regions of interest, N = 4, a total of 16 sub-regions of interest. The KLT tracking algorithm is used to track the detected sub-regions of interest to position the sub-regions of interest of each frame of picture.

[0055] Step 1.2, calculate the pixel mean of each sub-region of interest in the jth frame of face video image, so as to obtain the pixel mean of each sub-region of interest in J frames of face video image, and compose the pixel mean time sequence of each sub-region of interest. Let the pixel mean time sequence of any single sub-region of interest include: the red channel time sequence R = [R1, R2,..., RJ] of J frames of video image, the green channel time sequence G = [G1, G2,..., GJ] of J frames of video image and the blue channel time sequence B = [B1, B2,..., BJ] of J frames of video image. j ​..., R J ..., G j ..., G J ] and the time series of the blue channel B = [B1, B2,..., B j ..., B J ] in the J-frame video image; wherein R j represents the pixel mean value of the red channel in the jth frame in a single sub-region of interest, G j represents the pixel mean value of the green channel in the jth frame in a single sub-region of interest, and B j represents the pixel mean value of the blue channel in the jth frame in a single sub-region of interest.

[0056] Step 1.3, calculate the illumination intensity index, illumination change index and green signal signal-to-noise ratio index of a single sub-region of interest; and sort the illumination intensity indexes of the N×N sub-regions of interest in descending order, sort the illumination change indexes of the N×N sub-regions of interest in ascending order. Research shows that the green channel signal contains the best quality Blood Volume Pulse (BVP) signal, so the green signal signal-to-noise ratios of the N×N sub-regions of interest are also sorted in descending order. From the three sorted index sequences, select the top Q sub-regions of interest from each of the three sorted index sequences, and select M sub-regions of interest that exist in all the top Q sub-regions of interest from the three sorted index sequences as the final selected high-quality sub-regions of interest; M < Q < N×N.

[0057] In this example, the division of the region of interest is shown in Figure 2 , and the number of regions of interest M = 3.

[0058] Step two, extraction of heart rate signal under irregular motion scene;

[0059] Step 2.1, take the green channel signals of the M high-quality sub-regions of interest as multi-channel signals C = [G [1] ,G [2] ,..., G [m] ,..., G [M] ] T ; wherein G [m] represents the green channel signal of the mth sub-region of interest, i.e. the mth channel signal, and wherein, represents the pixel mean value of the jth frame of the green channel of the mth high-quality sub-region of interest, m = 1, 2,..., M; in this example, M = 3.

[0060] Step 2.2, Since the collected signals have different degrees of fluctuations, the uneven waveform of the signal will affect the final calculation result, so the signal C is de-trended. Set the de-trend parameter as λ, and de-trend the multi-channel signal C to obtain the de-trended multi-channel signal , where represents the de-trended mthchannel signal.

[0061] Step 2.3, Since the normal heart rate range is 0.75-3Hz (corresponding to heart rate of 45-180bpm), each channel signal in the de-trended multi-channel signal is normalized and band-pass filtered to obtain a filtered multi-channel signal , where represents the filtered mthchannel signal.

[0062] Step 2.4, Each channel signal in the filtered multi-channel signal is up-sampled to obtain an up-sampled transposed multi-channel signal matrix , where represents the up-sampled transposed mthchannel signal, and is a 1xJ matrix.

[0063] Step 2.5, The mthchannel signal in the multi-channel signal matrix X is decomposed by using a non-negative matrix factorization method to obtain an mthfeature matrix T [m] ∈R 1×L and an mthgain matrix V [m] ∈R L×J , and the lthnon-negative vector in the mthfeature matrix T [m] is denoted as t l [m] , and the non-negative vector in the lthrow and jthcolumn of the mthgain matrix V [m] is denoted as v l,j [m] , L is the number of feature vectors, and l=0, 1, 2, …, L; in a specific example, the number of feature vectors is L=10.

[0064] Step 2.6, The loss function Q of the improved independent vector analysis method is established by using formula (1):

[0065]

[0066] In formula (1), y j [m] represents the jthsource component estimation vector of the mthchannel, W represents the demixing matrix, and W=[w [1] , w[2] ..., w [m] ..., w [M] ] H where w [m] represents the component vector of the unmixing matrix of the mth channel; H represents the conjugate transpose;

[0067] Step 2.7, define the current iteration number as inter, and initialize inter = 1;

[0068] Let W inter be an all-one vector; initialize the non-negative vector in the mth eigenmatrix T [m] of the interth iteration as t l [m,inter = t l [m] ; initialize the non-negative vector in the lth row and jth column of the mth gain matrix V [m] of the interth iteration as v l,j [m,inter = v l,j [m] ;

[0069] Step 2.8, calculate the jth estimated variance r j [m],inter of the interth iteration by using formula (2):

[0070]

[0071] Step 2.9, calculate the intermediate variable Z [m],inter of the interth iteration by using formula (3):

[0072]

[0073] Step 2.10, calculate the component vector w' [m],inter of the unmixing matrix of the mth channel after updating of the interth iteration by using formula (4):

[0074] w' [m],inter = (WZ [m],inter ) -1 e [m] (4)

[0075] In formula (4), e [m] represents the mth unit vector;

[0076] Step 2.11, calculate the component vector w" [m],inter of the unmixing matrix of the mth channel after updating again of the interth iteration by using formula (5):

[0077]

[0078] Step 2.11: Calculate the estimated vector y of the j-th source component of the m-th channel in the inter-th iteration using equation (6). j [m],inter :

[0079] y j [m],inter =(w” [m],inter ) H X (6)

[0080] Step 2.12: Calculate the non-negative vector t of the l-th characteristic matrix in the (inter+1)-th iteration using equation (7). l [m ],inter+1 :

[0081]

[0082] Step 2.13: Calculate the non-negative vector v of the l-th gain matrix in the (inter+1)-th iteration using equation (8). lj [m ],inter+1 :

[0083]

[0084] Step 2.14: After assigning inter+1 to inter, return to step 2.8 and execute sequentially until inter > inter_max, thus obtaining the estimated vector y of the j-th source component under the inter_max-th iteration. j [m,inter_max And as the j-th estimated vector y of the final source component of the m-th channel. j [m] Thus, the source component estimation signal Y of the m-th channel is obtained. [m] =[y1 [m] ,y2 [m] ,…,y j [m] ,…,y J [m] ].

[0085] Step 3: Calculate heart rate;

[0086] The source component estimation signal Y of the m-th channel is calculated using the Fast Fourier Transform. [m] The power spectral density distribution is determined, and the source component with the most frequencies within a set range is identified as the blood volume pulse signal, thereby calculating the dominant frequency f of the blood volume pulse signal. maxtime x f max wherein time is the unit time.

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

[0088] In this embodiment, a computer readable storage medium has a computer program stored thereon, and the computer program, when executed by a processor, performs the steps of the above method.

[0089] To verify the effectiveness of the method, the performance of the method is verified on the public databases UBFC-RPPG and UBFC-PHYS. The UBFC-RPPG database contains a small amount of slight motion, and the UBFC-PHYS database is divided into three sub-datasets. The UBFC-PHYS1 is a static case, the UBFC-PHYS2 is an interview case, and the UBFC-PHYS3 is a digital game case. The heads of the subjects in the PHYS2 and the PHYS3 contain a large amount of motion. The method uses conventional evaluation indicators to evaluate the heart rate detection performance, including the mean absolute error HR mae (Mean Absolute Error, MAE), the root mean square error HR rmse (root mean square error, RMSE), the Pearson correlation coefficient (Pearson’s Correlation Coefficient, R), and the standard deviation HR sd (Standard Deviation).

[0090] Table 1 UBFC-RPPG database comparison test results

[0091] Method HR mae (bpm) HR rmse (bpm) HR sd (bpm) R CHROM 8.63 10.79 6.39 0.86 POS 6.60 8.20 4.87 0.88 ICA 6.47 11.74 9.79 0.83 SCF 4.06 5.52 3.37 0.95 IVA 4.06 5.58 3.82 0.95 NMF-IVA 2.81 4.19 3.10 0.97

[0092] Table 2 UBFC-PHYS1 database comparison test results

[0093] Method HR mae (bpm) HR rmse (bpm) HR sd (bpm) R CHROM 7.62 10.35 7.00 0.78 POS 7.98 10.50 6.83 0.74 ICA 2.83 5.41 4.60 0.91 SCF 5.27 10.45 9.02 0.69 IVA 5.23 10.30 8.88 0.70 NMF-IVA 1.16 2.08 1.76 0.98

[0094] Table 3 UBFC-PHYS2 database comparison test results

[0095]

[0096]

[0097] Table 4 UBFC-PHYS3 database comparison test results

[0098] Method HR mae (bpm) HR rmse (bpm) HR sd (bpm) R CHROM 12.15 17.27 12.28 0.53 POS 18.59 25.15 16.94 -0.02 ICA 11.30 15.29 10.30 0.59 SCF 12.70 21.46 17.30 0.27 IVA 11.66 16.87 12.22 0.57 NMF-IVA 5.32 6.75 4.16 0.91

[0099] Table 1 presents the results of the proposed method and several other methods on the UBFC-RPPG database. As can be seen from Table 1, the proposed method achieved the best results in all four metrics, specifically HR... mae 2.81 bpm, HR rmse Reduced to 4.19 bpm, HR sd The heart rate was reduced to 3.10 bpm, and the R-squared value increased to 0.97. Under static or minimally dynamic conditions, the algorithm combining NMF and IVA demonstrates strong separation of regular, periodic heart rate signals, enabling it to extract heart rate signals more accurately than model-based methods and traditional blind source separation methods. Overall, the proposed method achieved near-best results on the UBFC-RPPG database, demonstrating its superior performance on the relatively static UBFC-RPPG database.

[0100] Tables 2 to 4 present the results of the proposed method and several other methods on the UBFC-RPPG database. Table 2 shows the results from the PHYS1 database, representing the static case; Table 3 shows the results from the PHYS2 database, representing the case with significant head movement; and Table 4 shows the results from the PHYS3 database, representing the case with minimal head movement. The results from the three tables show that the Single Channel Filtering (SCF) method achieves good results in static cases, but its performance deteriorates significantly compared to blind source separation methods and model-based methods when motion interference is present. The model-based method (CHROM, POS) outperforms the traditional blind source separation ICA algorithm when there is significant motion, indicating that the combination of motion noise and pulse signals does not conform to a linear model when the motion amplitude is large. Independent Vector Analysis (ICA), as a joint blind source separation method, performs well in various situations because it can extract pulse information from different regions of interest. Compared to the IVA algorithm, the method proposed in this invention can not only extract pulse information from different regions of interest and has a good suppression effect on irregular motion noise, but also extract regular heart rate signals, greatly increasing the accuracy of heart rate calculation. Experimental results prove that the method proposed in this invention is superior in removing irregular motion artifacts with large motion amplitude.

[0101] Figure 3 The presentation compares the waveform of the BVP signal obtained by the proposed method as an improvement to IVA and the IVA algorithm when there is significant head movement. Figure 4The power spectrum of the signals obtained by the two algorithms is compared.From the figure, it can be seen that the fitting degree of the result obtained by the framework proposed in the application to the reference value is higher than that of the IVA method.In the case of large motion amplitude, the waveform recovered by the IVA algorithm has a larger error compared with the reference signal, but the method proposed in the application can recover a waveform with a high degree of approximation to the reference waveform.Comparing the power spectrum, it can be found that compared with IVA, the signal obtained by the method proposed in the application has a clearer main frequency, and the heart rate calculated through the main frequency is also closer to the true value.

[0102] In summary, the experiments on the two databases prove that the method proposed in the application is more suitable for heart rate detection under the presence of irregular motion artifacts.The motion-robust non-contact heart rate detection method proposed in the application improves the accuracy of video heart rate detection in irregular motion scenes and has good motion robustness.

Claims

1. A video heart rate detection method for removing irregular motion artifacts, characterized in that, is performed according to the following steps: Step one, screening high-quality sub-regions of interest from the face region of the face video image of the subject for extracting the pixel mean time series information of the green channel; Step 1.1, using a face feature point detection method to identify the jth frame face region from the jth frame face video image of the subject, and dividing the jth frame face region into N×N sub-regions of interest, using a feature point tracking algorithm to identify and locate the face region in each frame of face video image, thereby obtaining N×N sub-regions of interest in each frame of face video image; j=1, 2,..., J; J represents the total number of frames of face video image; Step 1.2, calculating the pixel mean of each sub-region of interest in the jth frame of face video image, thereby obtaining the pixel mean of each sub-region of interest in J frames of face video image, and forming the pixel mean time series of each sub-region of interest; Let the time series of pixel mean value of any single sub-region of interest include: the time series of red channel in J frame video image R = [R1, R2,..., R j ,...,R J ], the time series of green channel in J frame video image G = [G1, G2,..., G j ,...,G J ] and the time series of blue channel in J frame video image B = [B1, B2,..., B j ,...,B J ]; wherein R j represents the pixel mean value of the jth frame of red channel in a single sub-region of interest, G j represents the pixel mean value of the jth frame of green channel in a single sub-region of interest, and B j represents the pixel mean value of the jth frame of blue channel in a single sub-region of interest; Step 1.3, calculating the illumination intensity index, illumination change index and green signal signal-to-noise ratio index of a single sub-region of interest; and sorting the illumination intensity index of N×N sub-regions of interest in descending order, sorting the illumination change index of N×N sub-regions of interest in ascending order, and sorting the green signal signal-to-noise ratio of N×N sub-regions of interest in descending order, selecting the top Q sub-regions of interest from the sorted three index sequences, and selecting M sub-regions of interest that exist in all the top Q sub-regions of interest from the sorted three index sequences as the final high-quality sub-regions of interest; M<Q<N×N; Step two, extraction of heart rate signal under irregular motion scene; Step 2.1, taking the green channel signals of M high-quality region-of-interest subdomains as a multi-channel signal C = [G [1] ,G [2] ,...,G [m] ,...,G [M] ] T ; wherein G [m] represents the green channel signal of the mth region-of-interest subdomain, i.e., the mth channel signal, and wherein, represents the pixel mean value of the jth frame of the green channel of the mth high-quality region-of-interest subdomain, m = 1, 2,..., M; Step 2.2, set the detrending parameter as λ, and detrend the multi-channel signal C to obtain a detrended multi-channel signal denoted as C wherein, denotes the mth channel signal after detrending. Step 2.

3. Normalized band-pass filtering each channel signal in the de- trended multi-channel signal to obtain a filtered multi-channel signal denoted as wherein and denotes the m-th channel signal of the filtered multi-channel signal. Step 2.

4. upsampling each of the filtered multi-channel signals to obtain an upsampled transposed multi-channel signal matrix wherein represents the mth channel signal after upsampled transposition, and is a 1 x J matrix;​ Step 2.5, decompose the signal of the mth channel in the multi-channel signal matrix X by using the non-negative matrix factorization method to obtain the mth feature matrix T [m] ∈R 1×L and the mth gain matrix V [m] ∈R L×J , and the mth feature matrix T [m] The lth non-negative vector in the mth feature matrix T l [m] The non-negative vector in the lth row and jth column of the mth gain matrix V [m] l,j [m] L is the number of feature vectors, and l = 0, 1, 2, …, L.​​ Step 2.6, using formula (1) to establish the loss function Q of the improved independent vector analysis method: In formula (1), y j [m] represents the jth source component estimation vector of the mth channel, W represents the demixing matrix, and W = [w [1] ,w [2] ,...,w [m] ,...,w [M] ] H where w [m] represents the component vector of the demixing matrix of the mth channel; H represents the conjugate transpose; Step 2.

7. Solve the optimization problem for the loss function Q to obtain the source component estimate signal Y for the mth channel [m] ; Step three, calculation of heart rate; The power spectral density distribution of the mth channel source component estimation signal Y [m] is calculated using fast Fourier transform, and the source component estimation signal with the most frequency within the set range is determined as the blood volume pulse signal, so as to calculate the main frequency f max of the blood volume pulse signal for estimating the heart rate value.

2. The video heart rate detection method of removing irregular motion artifacts according to claim 1, characterized in that, The step 2.7 includes: Step 2.7.1, defining the current iteration number as inter, and initializing inter=1; Let W inter be the all-ones vector; initialize the m-th feature matrix T [m] at the inter-th iteration as the non-negative vector t l [m],inter l [m] ; initialize the m-th gain matrix V [m] at the inter-th iteration as the non-negative vector v l,j [m],inter l,j [m] ;​​ Step 2.7.2, computing the jth estimated variance r at the inter iteration using formula (2) j [m],inter : Step 2.7.

3. Compute the intermediate variable Z at the inter-th iteration using formula (3) [m],inter : Step 2.7.4, calculating the component vector w' of the updated demixing matrix for the mth channel at the inter iteration using formula (4) [m],inter : w' [m]，infer = (Wz [m]，infer)-1 e [m] (4) In formula (4), e [m] represents the mth unit vector; Step 2.7.5, calculating the component vector w of the demixing matrix of the mth channel after the update again at the inter iteration [m],inter : Step 2.7.

6. Compute the estimate vector y of the jth source component of the mth channel at the inter iteration using formula (6) j [m],inter : y j [m]，inter = (w" [m]，inter ) H x (6) Step 2.7.

7. Compute the non-negative vector t for the lth feature matrix at the (inter+1)th iteration using formula (7) l [m],inter+1 : Step 2.7.

8. Compute the non-negative vector v of the lth gain matrix at the (inter+1)th iteration using formula (8) lj [m ],inter+1 : Step 2.7.9, after inter+1 is assigned to inter, return to step 2.8 to be executed sequentially until inter > inter_max, so as to obtain the jth source component estimation vector y of the inter_maxth iteration j [m],inter_max and the jth estimation vector y of the final source component of the mth channel j [m] ; so as to obtain the source component estimation signal Y of the mth channel [m] = [y1 [m] , y2 [m] ,..., y j [m] ,..., y J [m] ].

3. 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 video heart rate detection method of claim 1 or 2, and the processor is configured to execute the program stored in the memory.

4. 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 video heart rate detection method of claim 1 or 2.

Images