Non-contact method for evaluating cardiopulmonary coupling

Abstract

The present invention discloses a non-contact cardiopulmonary coupling evaluation method, which correctly detects the pulse signal and breathing signal of the subject by shooting a video image of the upper body of a human body in a quiet state, and performs time domain and frequency domain analysis of cardiopulmonary coupling to evaluate and measure the cardiopulmonary coupling strength. The present invention shoots a video of the upper body of a subject through a camera, performs Gaussian pyramid decomposition and Laplace pyramid decomposition on the video, amplifies the pulse-related color changes and upper body movement information caused by breathing contained in the video, automatically extracts the head and face regions of interest containing pulse information and the chest and abdomen regions of interest containing breathing movement, calculates the brightness information in the regions of interest respectively, thereby extracting the pulse signal and the breathing signal, and performs time domain and frequency domain analysis of cardiopulmonary coupling through synchronization processing of the two, and performs cardiopulmonary coupling evaluation according to the coupling strength.

Description

Technical Field

[0001] The present invention relates to a non-contact cardiopulmonary coupling evaluation device and method, which monitors the pulse and respiratory signals of a subject in a non-contact manner, and measures and analyzes the cardiopulmonary coupling strength, thereby realizing early warning of respiratory and cardiovascular diseases and determining the respiratory status. Background Art

[0002] Cardiopulmonary coupling analysis was proposed by CK Peng of Harvard Medical School in 2005. It is mainly used in sleep structure and quality analysis and sleep apnea identification. It is usually based on ECG signals, and jointly analyzes heart rate and respiratory signals to obtain various time-frequency domain parameters, and then analyzes the relationship between each parameter and physiological mechanism. It is used in the study of sleep apnea, cardiovascular disease, etc.

[0003] Cardiopulmonary coupling analysis is used to evaluate the coupling strength between heart rate and respiratory signals. There are many analysis methods, including time domain analysis, frequency domain analysis, nonlinear analysis and other methods. In time domain analysis, the number of heartbeats in a period of time in each respiratory cycle is often counted to obtain the heartbeat rate, and respiratory diseases are diagnosed based on the variability of the heartbeat rate. In frequency domain analysis, the strength of cardiopulmonary coupling is usually expressed by the product of the square of the cross-power spectrum of the RR interval sequence and the respiratory signal in the ECG signal and the coherence. The cross-power spectrum of the RR interval sequence and the respiratory duration sequence signal refers to the product of the power of the two at the same frequency. The cross-power spectrum can be used to measure whether the two signals have a relatively large amplitude synchronously at this power. The coherence of the RR interval sequence and the respiratory duration sequence signal reflects the synchronization of the oscillations of the two signals, which is calculated by the cross-correlation function of the two signals.

[0004] Traditional cardiopulmonary coupling analysis is based on ECG signals, which can be obtained by contact measurement methods such as clinical electrocardiography and polysomnography. Although these methods can obtain high signal-to-noise ratio and high stability of heart rate and respiratory signals, thus providing a basis for cardiopulmonary coupling analysis, the measurement can easily cause discomfort to patients because it must be measured by contact. Especially in the diagnosis of respiratory diseases such as sleep apnea, the measurement process and measurement device may affect the patient's sleep, which increases the difficulty of diagnosis and reduces the diagnostic efficiency.

[0005] Cardiopulmonary coupling analysis often requires long-term measurements of more than 30 minutes. In this case, contact measurement will seriously affect the experience and cause serious interference to the subjects during sleep. At the same time, the electrode may fall off during the detection process, resulting in unstable signals and reducing the accuracy and stability of the test results.

[0006] On the other hand, the above-mentioned measuring equipment all needs to be carried out in a hospital environment, requires professional technicians to operate, and requires professional physicians to manually score the obtained measurement values, so the examination fee is relatively expensive.

[0007] Therefore, non-contact, low-load measurement devices can better meet the conditions of long-term, low-interference cardiopulmonary coupling measurement, and realize the measurement of relevant physiological signals and evaluation of cardiopulmonary coupling strength in a relaxed and comfortable environment such as home, thereby achieving higher diagnostic accuracy and effectiveness.

[0008] In response to the above needs, people have proposed many methods, such as using Doppler radar to measure chest and abdominal movements, which can achieve non-contact breathing measurement, but the accuracy of pulse or heart rate measurement is low and it is difficult to meet the requirements of cardiopulmonary coupling measurement. In addition, Doppler radar has a certain amount of electromagnetic radiation and is not suitable for long-term home use.

[0009] In addition, a method for detecting the respiration or pulse of a subject using a video analysis method is also known. However, due to the different acquisition sites and acquisition parameters required for video analysis of the subject's respiration and pulse signals, there is still a lack of research on simultaneously acquiring the subject's respiration and pulse through video. On the other hand, traditional cardiopulmonary coupling analysis is based on ECG signals, which can naturally synchronize respiration and heart rate signals; however, there is still a lack of relevant research on the synchronization method of respiration and pulse signals obtained by video analysis. Therefore, although the video-based analysis method provides a very promising technical method for non-contact cardiopulmonary coupling analysis, there are still many technical difficulties and bottlenecks, which make it difficult to meet the needs of cardiopulmonary coupling evaluation. Summary of the invention

[0010] To solve the above problems, the present invention proposes a non-contact cardiopulmonary coupling evaluation method, which uses video monitoring of the head, face and upper body images of the subject to correctly detect the pulse signal and respiratory signal, and can perform time domain and frequency domain analysis of cardiopulmonary coupling to achieve the evaluation of the cardiopulmonary coupling strength.

[0011] The non-contact cardiopulmonary coupling evaluation method of the present invention comprises the following specific steps:

[0012] Step 1: Obtain a video of the subject in a quiet state, including the subject's head, face and upper body.

[0013] Step 2: Amplify small changes in the video, including small color changes and small movement changes; amplify the color changes of the facial image in the video to extract the pulse signal, and amplify the small movements of the upper body image in the video to extract the breathing signal.

[0014] Step 3: Determine the region of interest of the enlarged video image, including the region of interest of the pulse signal and the region of interest of the respiratory signal.

[0015] Step 4: Extract pulse and respiratory signals from the image sequences in the regions of interest of the pulse and respiratory signals respectively, and synchronize the pulse and respiratory signals through cross-correlation operations.

[0016] Step 5: Perform cardiopulmonary coupling analysis on the synchronized breathing and pulse signals, including time domain analysis and frequency domain analysis, to achieve cardiopulmonary coupling measurement.

[0017] Based on the above method, the present invention proposes a non-contact cardiopulmonary coupling evaluation device, including a video image acquisition device, a video magnification device, a region of interest selection device, a video signal processing device, and a cardiopulmonary signal synchronization and analysis device.

[0018] The video image acquisition device is a camera, which is used to capture the video of the head, face and upper body of the subject.

[0019] The video magnification device comprises a pyramid decomposition unit, a time domain filtering unit and an amplification unit; wherein the pyramid decomposition unit amplifies the color change and contour change in the video image by Gaussian video amplification method and Laplace video amplification method in parallel. The time domain filtering unit performs time domain bandpass filtering on the brightness information of each layer of the pyramid decomposition image; outputs the frequency band information of interest of pulse and the frequency band information of interest of breathing; the amplification unit amplifies the frequency band information of interest of each layer of the pyramid, performs pyramid reconstruction on the amplified signal and adds it back to the original image, and outputs the amplified pulse video frame image and the breathing video frame image.

[0020] The region of interest selection device selects regions of interest for the Gaussian video magnified image sequence and the Laplace video magnified image sequence, respectively, and outputs regions of interest for respiratory and pulse signals.

[0021] The video signal processing device extracts the respiratory signal and the pulse signal from the respiratory and pulse signal regions of interest respectively.

[0022] The cardiopulmonary signal synchronization and analysis device is used to extract and synchronously process the pulse signal and the respiratory signal, and then perform time domain analysis and frequency domain analysis of the cardiopulmonary coupling.

[0023] The advantages of the present invention are:

[0024] 1. The non-contact cardiopulmonary coupling evaluation method of the present invention uses video to simultaneously obtain the subject's breathing and pulse for analysis, which can be monitored for a long time, avoiding the discomfort of long-term sticking of electrodes in traditional measurement methods, and also avoiding interference with the subject during sleep monitoring.

[0025] 2. The non-contact cardiopulmonary coupling evaluation method of the present invention can measure the cardiopulmonary coupling function without contact, and is suitable for the evaluation of respiratory diseases and respiratory function, especially for some highly contagious respiratory diseases, the non-contact measurement can greatly reduce the risk of infection.

[0026] 3. The non-contact cardiopulmonary coupling evaluation method of the present invention can evaluate the autonomic nervous function of cardiopulmonary coupling without load, and is suitable for scenarios such as mental fatigue evaluation, tension level evaluation, and lie detection.

[0027] 4. Compared with other non-contact measurement methods, the non-contact cardiopulmonary coupling evaluation method of the present invention has the advantage of low cost and is conducive to market promotion. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] Figure 1 It is a module diagram of the non-contact cardiopulmonary coupling evaluation method of the present invention;

[0029] Figure 2 It is a schematic diagram of the principle and results of the Gaussian pyramid amplification algorithm of the non-contact cardiopulmonary coupling evaluation method of the present invention;

[0030] Figure 3 It is a schematic diagram of the principle and results of the Laplace pyramid amplification algorithm in the non-contact cardiopulmonary coupling evaluation method of the present invention;

[0031] Figure 4 It is a flow chart of the pulse signal processing algorithm in the non-contact cardiopulmonary coupling evaluation method of the present invention;

[0032] Figure 5 It is a flow chart of the respiratory signal processing algorithm in the non-contact cardiopulmonary coupling evaluation method of the present invention;

[0033] Figure 6 It is a schematic diagram of a method for calculating time domain parameters of cardiopulmonary coupling in the non-contact cardiopulmonary coupling evaluation method of the present invention;

[0034] Figure 7 It is a schematic diagram of a method for calculating frequency domain parameters of cardiopulmonary coupling in the non-contact cardiopulmonary coupling evaluation method of the present invention;

[0035] Figure 8 Schematic diagram of the experimental scene of the non-contact cardiopulmonary coupling evaluation method of the present invention. DETAILED DESCRIPTION

[0036] The present invention is further described in detail below with reference to the accompanying drawings.

[0037] The non-contact cardiopulmonary coupling evaluation method of the present invention comprises the following specific steps:

[0038] Step 1: Take a video of the subject in a quiet state

[0039] A video of the subject in a quiet state under natural lighting conditions is captured using a camera (webcam or industrial camera), and the video image includes the subject's head, face and upper body.

[0040] Step 2: Zoom in on small changes in the video shot in step 1, including small color changes and small movement changes.

[0041] The video image (video) shot in step 1 is input into the video magnification device, and two techniques are used in parallel to magnify the video: one is the Gaussian video magnification method, which is used to magnify the color changes in the video image; the other is the Laplace video magnification method, which is used to magnify the contour changes in the video image. Finally, the video magnification device outputs the magnified video frame image (Mag frames).

[0042] First, the video magnification device in the device adopts two improved Euler video magnification methods, such as Figure 1 As shown, the video is magnified by Gaussian video and Laplace video respectively.

[0043] Among them, the Gaussian video magnification method is also called the Gaussian pyramid decomposition and magnification method. Its principle is to extract and amplify the image area or signal components of human interest through spatial filtering and sampling of different resolutions for a video image. Since blood absorbs different wavelength spectra in visible light, the pulse can cause the absorption and reflection of visible light on the head, face and other skin surfaces with rich blood vessels to show the same change pattern as the pulse frequency. However, the change of this reflected light is very weak in the video, which is not conducive to the extraction of the pulse signal. Therefore, the Gaussian video magnification method is used to preprocess the video, and the amplitude of the change of the pulse-related signal in the amplified pulse signal is more obvious, which is conducive to better extraction of the pulse signal from the video. The specific method of Gaussian pyramid decomposition and amplification is:

[0044] 1) Gaussian pyramid decomposition

[0045] Spatial filtering is performed on the image sequence that constitutes the video to obtain video images at different spatial resolutions; this helps to reduce noise. Video images present different SNRs (signal-to-noise ratios) at different spatial resolutions. According to research on human vision, the best signal-to-noise ratio can be achieved under specific spatial resolution conditions; therefore, in order to prevent distortion, these video images with different spatial resolutions should use different magnification factors. For objects or areas of interest, it is often necessary to perform image analysis such as feature extraction at a specific spatial resolution to achieve the best effect. Therefore, when amplifying images of a specific spatial resolution during the magnification process, the magnification factor should be increased, while for images of other spatial resolutions, the magnification factor should be reduced.

[0046] The principle of Gaussian pyramid decomposition is as follows Figure 2 As shown, let level j be the input image. First, the input image is subjected to Gaussian low-pass filtering, and the filtering function is:

[0047]

[0048] In the formula, D(u,v) represents the distance between the pixel point (u,v) and the origin of Fourier transform, u represents the column coordinate of the image pixel, v represents the row coordinate of the image pixel, and δ represents the bandwidth of the Gaussian filter, which determines the degree of smoothing. The larger δ is, the wider the frequency band corresponding to the filter is, and the better the degree of smoothing is.

[0049] Then downsampling is performed. Downsampling can be achieved by discarding the even rows and even columns in the image, so that the length and width of the image are reduced to half, and the area is reduced to one-fourth of the original image. The downsampling result at this time is the j-1 level approximate image. Repeating the above process with the j-1 level image as input can obtain the Gaussian pyramid.

[0050] 2) Frequency domain filtering

[0051] In the images of different spatial resolutions obtained by Gaussian pyramid decomposition, the color, brightness and other information of each pixel comes from the reflected light entering the camera, which is affected by factors such as the incident light and the reflector. Assuming that the incident light remains unchanged, the changes in the spectrum and brightness of the reflected light reflect the changes in the properties of the reflector. When shooting a video of the human head and face, the pulsation of blood causes changes in the filling degree of subcutaneous capillaries, and the blood flow in the capillaries has a strong absorption of specific spectra in the incident light, causing the spectrum of the reflected light to change periodically with the pulsation of the blood. Therefore, the pulse-related signals in the reflected light can be extracted through frequency domain filtering. The specific method is:

[0052] Frequency domain bandpass filtering is performed on each pixel of the image at each spatial resolution in order to extract the frequency band signal of interest (PLR). According to the pulse frequency range of adults, which is 30 to 240 times / min, the passband frequency of the bandpass filter is set to 0.5 to 4 Hz. If only the pulse frequency of a specific object needs to be extracted, a narrower passband can be selected to remove more motion interference and obtain a better extraction effect. The result of the above-mentioned bandpass filtering is the approximation of the change signal of interest. In the amplification of the pulse signal, the corresponding change is the change in image color.

[0053] 3) Video zoom

[0054] Amplify the signal of the frequency band of interest. Here we take amplifying a one-dimensional signal as an example. The two-dimensional image signal is similar. Assuming I(x,t) represents the brightness value at any position x and any time t, and δ(t) is a displacement function, then we have:

[0055] I(x,t)=f(x+δ(t)), I(x,0)=f(x) represents the brightness distribution corresponding to time 0.

[0056] The purpose of the final video amplification is to synthesize the signal:

[0057]

[0058] Among them, α represents the magnification,

[0059] By performing a first-order Taylor series expansion on the changed signal f(x+δ(t)) at time t, we can approximately obtain:

[0060]

[0061] Let B(x, t) represent the frequency domain bandpass filtering result of I(x, t) at each x position. Assuming that the motion signal of interest δ(t) is contained in the passband range, we have:

[0062]

[0063] Amplify the frequency domain filtering result B(x,t) by α times and add it back to the original signal, and we get:

[0064]

[0065] Combining equations (3), (4), and (5) yields:

[0066]

[0067] Then, according to the first-order Taylor series expansion, we can approximate:

[0068]

[0069] The above formula gives the result (filtered) after the spatial brightness signal δ(t) of the original signal f(x) is amplified by α times.

[0070] 4) Output the enlarged pyramid layer image

[0071] The traditional Gaussian video amplification technology includes a video reconstruction step, in which the amplified frequency band signal (filtered) is subjected to pyramid reconstruction and added back to the original image to obtain the amplified image sequence (Magframes), which is composed of the amplified image sequence (Mag frames) and output as a result. However, the reconstruction step is relatively time-consuming to calculate, and the reconstructed video causes image information of different spatial resolutions to be aliased again, reducing the signal-to-noise ratio, which is not conducive to the application of real-time pulse signal extraction. Therefore, the present invention adopts direct output of the amplified image frame for subsequent analysis, abandoning the video reconstruction step of the traditional Gaussian video amplification method to improve the signal-to-noise ratio and calculation rate, and save computer storage space. So far, the Gaussian pyramid video amplification process is completed.

[0072] The Laplacian pyramid video magnification method is similar in principle to the Gaussian pyramid process, but there are the following differences:

[0073] 1) Laplacian pyramid decomposition is different from Gaussian pyramid decomposition in pulse signal amplification. Laplacian pyramid decomposition is based on Gaussian pyramid decomposition. The j-1 level approximate image is subjected to the inverse process of pyramid decomposition, that is, it is upsampled with the same step size as the pyramid decomposition and then passed through the same Gaussian low-pass filter to make a difference with the original j level image to obtain the j level image prediction residual. At this time, a layer of the Laplacian pyramid image is obtained. Then the j-1 level image is used as the input image and the above steps are repeated to obtain the complete Laplacian pyramid, as shown in Figure 3 As shown in Figure 2, each frame of the original video is decomposed into a spatial Laplace pyramid as the baseband output (PL).

[0074] 2) Frequency domain filtering

[0075] Laplace pyramid amplification is mainly used to extract respiratory motion, so the frequency band of human respiratory motion is used in frequency domain filtering. The passband frequency of the bandpass filter is set to 0.05-1Hz based on the fact that the respiratory rate of an adult at rest is about 14-20 times / min. If only the respiratory rate of a specific object needs to be extracted, a narrower passband can be selected to remove more motion interference and obtain better extraction effect.

[0076] The other steps in the Laplacian pyramid upscaling method are the same as those in the Gaussian pyramid upscaling method.

[0077] Step 3: Determine the region of interest, including the region of interest of the pulse signal and the region of interest of the respiratory signal.

[0078] The video image amplified by the above method is used as input to quickly determine the region of interest where the pulse signal appears in the video. The specific method is:

[0079] Generally speaking, the region of interest of the pulse signal appears on the human head and face, because the human head and face have the richest capillary network. Therefore, when selecting the region of interest of the pulse signal, a face recognition method based on a neural network is used to determine the human head and face area. For the image sequence in the head and face area, the inter-frame difference method is used to calculate the difference Dg in the green channel signal intensity between frames. The advantages of using the inter-frame difference method are: the algorithm is simple to implement and the program design complexity is low; it is not very sensitive to scene changes such as light, can adapt to various dynamic environments, and has good stability. Let the green channel signal intensity of the nth frame and the 1st frame in the image sequence be g n and g1, the grayscale value of the corresponding pixel in the two frames is recorded as g n (u, v) and g1(u, v), subtract the grayscale values ​​of the pixels corresponding to the two frames of images according to the following formula to obtain the difference image D g (u,v), where u and v represent the column and row coordinates of the image pixel respectively:

[0080] D g (u,v)=|g n (u,v)-g1(u,v)| (8)

[0081] According to the average color (R m , G m , B m ), calculate the difference between the color of each pixel in the face area and the average color In the formula, R, G, and B represent the red channel signal, green channel signal, and blue channel signal of the image pixel respectively; R m , G m , B m They are the average values ​​of the red, green and blue channel signals of all pixels in the face area.

[0082] Calculate D for each pixel g / D c , which is a function of time; the maximum value of this function within the time of video acquisition is obtained as the inter-frame change intensity of each pixel. According to the inter-frame change intensity of all pixels in the face area, they are sorted from large to small, and the ratio of a certain pixel ranking to the number of all pixels is recorded as λ, and the inter-frame change intensity of the pixel at the sorting position is set as the inter-frame change intensity threshold T λThe value of λ is generally 0.3-0.5, which can be adjusted according to parameters such as lighting conditions and spatial resolution of the image sequence. When the value of λ is small, the region of interest of the pulse signal is small, which is suitable for conditions with high image resolution and sufficient lighting; when the value of λ is large, the region of interest of the pulse signal is large, which is suitable for conditions with low image resolution and insufficient lighting. For the frame-to-frame change intensity greater than the threshold T λ The pixels with strong pulse signals are considered as the pixels with strong pulse signals. After performing appropriate closing operations on these pixels, regions with good continuity are generated, which can reduce the influence of local noise, especially the salt and pepper noise. These regions are regarded as the regions of interest of pulse signals.

[0083] The area of ​​interest for respiratory motion generally appears in the chest area of ​​the human body. In order to specifically determine these areas, the image sequence after Laplace video amplification is input, and the color image sequence is first converted into a grayscale image sequence. According to the brightness value of the grayscale image, the local Otsu method is used to select a suitable brightness threshold to convert the image sequence into a binary image sequence. The method is:

[0084] Set the threshold T, and binarize the pixels one by one according to the following formula to obtain the binary image R:

[0085]

[0086] Among them, the point with a grayscale value of 255 is the foreground (moving target point), and the point with a grayscale value of 0 is the background point. Subtract the first frame image from each frame of the binary image R to generate a differential image sequence, and extract the region with the most dramatic change (ROI) in the differential image sequence based on the change in the area of ​​the motion region, that is, the connected domain with the largest area. According to the artificial neural network, the trunk features in the original video are located, including the chest and abdominal areas, and the extracted area with the most intense respiratory movement (ROI) is intersected with the located chest and abdominal areas to obtain the ROI of the area where the chest movement of the subject is significantly changed. K And the ROI of the area with obvious abdominal movement changes F , which are the two regions of interest for respiratory signals.

[0087] Step 4: Extract pulse and respiratory signals from the image sequences in the regions of interest of the pulse and respiratory signals respectively, and synchronize the pulse and respiratory signals through cross-correlation operations.

[0088] like Figure 4As shown in the figure, for pulse signal extraction, the red channel signal R, the green channel signal G, and the blue channel signal B are obtained by calculating the mean values ​​of the three color channel intensities in the region of interest of the pulse signal. Independent component analysis is performed on the three color channel signals. The principle is to assume that the collected signal is a linear superposition of several noises, and each signal and noise is independent of each other, which is called an independent component, as shown in the following formula:

[0089]

[0090] In the formula, I R ,I G ,I B They are the intensities of the red, green and blue channel signals respectively; C1, C2 and C3 are three independent components respectively. The red, green and blue channel signals in the color image sequence can be regarded as the results of multi-channel observation of the pulse signal. Each observation result is a linear superposition of the signal and noise with different weights. Through the independence assumption, the independent signals and noise can be separated from the multi-channel observation results, thereby obtaining the desired signal from the multi-channel observation results. The present invention uses the independent component analysis method to calculate the three independent components of C1, C2 and C3 from the three color channel signals. According to the above assumptions, one of the independent components must be a signal related to the pulse, and the other two independent components are noise. However, among the three independent components obtained by the independent component analysis, it is impossible to determine which one is the signal and which is the noise. Studies have shown that the blood in the capillaries absorbs the green spectrum most obviously, so the pulsation of blood in the capillaries has the most drastic effect on the green channel signal, which can cause the green channel signal to produce periodic fluctuations. Therefore, the three independent component results are correlated with the green channel signal G, and the result with the highest correlation is considered to be the pulse-related signal (IPPG). The IPPG signal is then analyzed in the frequency domain, and the frequency corresponding to the spectrum energy peak in 0.5-4Hz is selected as the pulse frequency fp. A narrowband bandpass filter is further constructed with a passband frequency of fp±0.15Hz, and the IPPG signal is narrowband filtered to obtain the pulse signal Sp.

[0091] like Figure 5 As shown, for the respiratory signal, two respiratory signal regions of interest (ROI K and ROI F )Internal difference image brightness and (IROI K and IROI F ), where the value of the differential image pixel brightness value may be 0 (indicating that the pixel brightness values ​​of the two frames are the same), +1 (indicating that the pixel brightness value of the subsequent image is greater than that of the first frame), and -1 (indicating that the pixel value of the subsequent image is less than that of the first frame); and the differential image brightness and (IROIK and IROI F ) were used to extract chest and abdominal respiratory intensity (R K and R F The respiratory intensity of the chest and abdomen is weighted averaged, and the weight of the chest respiratory signal is set to 1, while the weight of the abdominal respiratory signal is the sum of the differential image brightness in the chest and abdominal regions of interest (IROI K and IROI F ), the calculated signal strength is the respiratory signal S r , the expression is:

[0092] S r =R K +R F *f ass (IROI K ,IROI F ) (11)

[0093] In the formula, f ass (IROI K ,IROI F ) represents the correlation between the differential image brightness and (IROIK and IROIF) within the thoracic and abdominal regions of interest.

[0094] Step 5: Perform cardiopulmonary coupling analysis on the synchronized pulse signal Sp and respiratory signal Sr, including time domain analysis and frequency domain analysis, to achieve cardiopulmonary coupling measurement.

[0095] First, the pulse signal Sp and the respiratory signal Sr are synchronized. The pulse signal Sp and the respiratory signal Sr are synchronized by calculating the cross-correlation function of the pulse signal Sp and the respiratory signal Sr in a respiratory cycle Tr:

[0096]

[0097] Wherein, τ and t are any moments in a respiratory cycle Tr.

[0098] If the cross-correlation function R pr (t) reaches the maximum value at time t = tm, then the pulse signal Sp(t) is time-shifted by tm, that is, the new pulse signal S p ′(t)=S p (tt m), the synchronization of cardiopulmonary signals can be achieved. The principle of synchronization is that the traditional cardiopulmonary coupling analysis is based on ECG signals. Since the ECG signal contains R waves that can describe ventricular depolarization, the amplitude of the ECG signal also changes with the respiratory movement. When the inhalation movement occurs, the chest increases, and the voltage difference between the ECG electrodes also increases, which is manifested as an increase in the amplitude of the ECG signal. Therefore, the synchronization of heart rate and respiratory signals can be achieved according to the changes in the position and amplitude of the R wave in the ECG signal. However, in non-contact measurement, ECG activity cannot be obtained, so innovative synchronization methods are needed. In this method, the extreme points in the pulse signal Sp indicate the strongest or weakest moment of the pulse's spectral absorption, corresponding to the maximum contraction moment and the end of diastole in the pulse, respectively; and the extreme points in the respiratory signal Sr indicate the corresponding moments of the end of inspiration and the end of exhalation in the respiratory movement. Therefore, the synchronization of pulse and respiration can be achieved by cross-correlation operation of the pulse signal Sp and the respiratory signal Sr.

[0099] After completing the synchronization of pulse and respiratory signals, cardiopulmonary coupling analysis is performed, including cardiopulmonary coupling time domain analysis and frequency domain analysis. The specific methods are as follows:

[0100] Cardiopulmonary coupling time domain analysis:

[0101] like Figure 6 As shown, for each cycle of the respiratory signal Sr, calculate the Sr signal gradient, and record the time when the gradient maximum value appears as the inhalation start time ti; delay a period of time ts based on ti, and calculate the number of times hit when the pulse signal Sp maximum value appears within a period of time Δt before and after ts (i.e., within the time range of ti+ts±Δt); continuously record the number of hits within 30 minutes, and divide the number of hits by the number of breaths within 30 minutes to obtain the heart rate. Change ts and draw an incidence variation curve with ts as the horizontal axis and heart rate as the vertical axis. Extract the maximum value Imax and minimum value Imin of the heart rate from the heart rate variation curve, and calculate the heart rate time variability index I v for:

[0102]

[0103] Stroke Time Variability Index I v It can be used to characterize the strength of cardiopulmonary coupling, I v The closer it is to 1, the higher the synchronization of cardiopulmonary activity. v The closer it is to 0, the lower the synchronization of cardiopulmonary activity.

[0104] Cardiopulmonary coupling analysis Frequency domain analysis:

[0105] like Figure 7As shown, the peak positions of the pulse signal Sp and the respiratory signal Sr are extracted respectively, and the time difference between each peak is obtained. With the peak occurrence time as the horizontal axis and the time difference between the peaks as the vertical axis, the pulse cycle sequence signal a(t) and the respiratory duration sequence signal b(t) are constructed; the cross-correlation function of the two signals is calculated:

[0106]

[0107] Then the cross-power spectrum S of the two signals is calculated by the cross-correlation function ab (ω), that is, the Fourier transform of the cross-correlation function is as follows:

[0108]

[0109] The autocorrelation function R of the pulse cycle sequence a(t) and the breathing duration sequence signal b(t) aa (t), R bb (t) and the auto-power spectral density S aa (ω), S bb (ω) can be calculated using the following formulas:

[0110]

[0111]

[0112]

[0113]

[0114] Then the coherence coefficient C of the two sequences can be calculated by the following formula: ab (ω):

[0115]

[0116] Finally, the coupling strength CPC(ω) of the two signals is calculated by the cross power spectrum square and coherence of the two sequence signals:

[0117]

[0118] By calculating the power of different frequency bands of the cardiopulmonary coupling strength CPC (ω), we can obtain the total power (TPC) in the frequency range of 0-0.4 Hz, the power in the frequency range of 0.04-0.1 Hz (low frequency power, LFC), the low frequency power divided by the total power (normalized low frequency power, LFC-norm), the power in the frequency range of 0.1-0.4 Hz (high frequency power, HFC), the high frequency power divided by the total power (normalized high frequency power, HFC-norm), and the ratio of low frequency power to high frequency power (MFC). The parameters of different frequency bands can be used to evaluate the cardiopulmonary coupling strength.

[0119] The measurement and evaluation of cardiopulmonary coupling can be achieved based on the heart rate variability index Iv in the time domain analysis results of cardiopulmonary coupling and the cardiopulmonary coupling strength CPC(ω) in the frequency domain analysis results.

[0120] Example

[0121] use Figure 1 The structure diagram of the cardiopulmonary coupling measurement device shown in the figure, first, use the video acquisition device to capture the video of the subject, fix the camera within a suitable distance from the subject, and ensure that the image contains the head, face and upper body of the subject. The scene diagram of video acquisition in this embodiment is as follows Figure 8 The video capture process can use sufficient natural light or other stable light sources.

[0122] Secondly, the captured video image (video) is input into the video magnification device. The number of spatial pyramid decomposition layers is determined according to the size of the input video frame image. Downsampling and upsampling are performed with a step size of 2. For example, if the input image size is 2048×2048, the number of spatial pyramid decomposition layers is 11 layers, forming Gaussian pyramid and Laplace pyramid video image sequences. Then, the two pyramid image sequences constructed for each image are respectively subjected to time domain filtering. The filter is designed as an infinite impulse response bandpass filter. The passband range can be selected by the user. The pulse signal passband range used in this embodiment is 0.5-4Hz, and the respiratory signal passband range is 0.05-0.5Hz. The filtered frequency band of interest is then amplified. The pulse video amplification factor of this embodiment is α. G =10, breathing video magnification α L The amplified frequency bands of interest of the pulse and respiratory signals are obtained respectively, and then the amplified frequency bands of interest are reconstructed by high-level pyramid and Laplace pyramid respectively, and the reconstructed images are added back to the original video images, and the amplified pulse video frame images and respiratory video frame images are output respectively.

[0123] Then, the amplified pulse video frame image and breathing video frame image are input, and the inter-frame difference method is used to determine the area in the video where the color change caused by the pulse is the most significant and the area where the respiratory movement change is the most obvious, and distinguish the area of ​​interest of the head and face pulse signal and the area of ​​interest of the chest and abdomen respiratory signal.

[0124] Next, for each frame of the pulse signal region of interest, the average signal intensities IR, IG, and IB of the red, green, and blue colors of the image are extracted respectively. After preprocessing such as detrending, whitening, smoothing, and denoising, the independent component analysis method is used to extract the pulse-related IPPG signal from IR, IG, and IB, and the IPPG signal is analyzed in the frequency domain to determine the pulse frequency. A bandpass filter is constructed with the pulse frequency as the center frequency, and the passband frequency is fp±0.15Hz. The IPPG signal is bandpass filtered to obtain the pulse signal Sp.

[0125] Two regions of interest (ROI) for respiratory signals in the chest and abdomen K and ROI F ) The sum of the brightness information of the pixel values ​​within (I ROI K and I ROI F ), and select the frequency band corresponding to the respiratory signal of interest to sum the brightness information of the region of interest (I ROI K and I ROI F ) were band-pass filtered, and the passband frequency of the band-pass filter was set to 0.05-1 Hz according to the breathing rate of an adult at rest, which is about 14-20 times / min. The breathing signal curve (R K and R F ). The respiratory signal curve R of the chest and abdomen K and R F Perform weighted averaging to obtain the respiratory signal Sr.

[0126] The cardiopulmonary signal synchronization and analysis device synchronizes the pulse signal Sp and the respiratory signal Sr, and realizes the time domain analysis and frequency domain analysis of the cardiopulmonary coupling, and obtains the heart beat time variability index Iv and the cardiopulmonary coupling strength CPC (ω) respectively. By calculating the power of different frequency bands of the cardiopulmonary coupling strength CPC (ω), the total power (TPC), low frequency power (LFC), normalized low frequency power (LFC-norm), high frequency power (HFC), normalized high frequency power (HFC-norm), and the ratio of low frequency power to high frequency power (MFC) can be obtained. The cardiopulmonary coupling strength is evaluated according to the power of different frequency bands, and early warning of respiratory diseases and cardiovascular diseases is realized, or the respiratory state, such as respiratory obstruction and apnea, is determined.

Claims

1. A non-contact cardiopulmonary coupling evaluation method, characterized in that: This is accomplished by: Step 1: Obtain a video of the subject in a quiet state, including the subject's head, face and upper body; Step 2: Amplify small changes in the video, including small color changes and small movement changes; amplify the color changes of the facial image in the video to extract the pulse signal, and amplify the small movements of the upper body image in the video to extract the breathing signal; The method to amplify small changes in the video is: The video image captured in step 1 is input into a video magnification device, and two techniques are used in parallel to perform video magnification: one is a Gaussian video magnification method, which is used to magnify color changes in the video image; The other is a Laplace video magnification method, which is used to magnify the contour changes in the video image, and finally the video magnification device outputs the magnified video frame image; First, the specific process of Gaussian video magnification method is: 1) Gaussian pyramid decomposition The image sequence constituting the video is spatially filtered to obtain video images at different spatial resolutions; video images present different SNRs at different spatial resolutions, and the best signal-to-noise ratio is achieved under specific spatial resolution conditions; therefore, video images with different spatial resolutions use different magnifications; For the target object or area of ​​interest, feature extraction needs to be performed at a specific spatial resolution. Therefore, when a specific spatial resolution image is amplified during the zooming process, the magnification factor is increased, while for other spatial resolution images, the magnification factor is reduced; In the Gaussian pyramid decomposition method, let level j be the input image. First, the input image is subjected to Gaussian low-pass filtering, and the filtering function is: In the formula, D(u,v) represents the distance between the pixel point (u,v) and the origin of Fourier transform, u represents the column coordinate of the image pixel point, v represents the row coordinate of the image pixel point, and δ represents the bandwidth of the Gaussian filter; Then downsampling is performed. Downsampling is achieved by discarding the even rows and even columns in the image. The length and width of the image are reduced to half, and the area is reduced to one quarter of the original image. The downsampling result at this time is a j-1 level approximate image. The j-1 level image is used as input and the above process is repeated to obtain a Gaussian pyramid. 2) Frequency domain filtering From the images of different spatial resolutions obtained by Gaussian pyramid decomposition, the pulse-related signals in the reflected light are extracted by frequency domain filtering. The specific method is as follows: Perform frequency domain bandpass filtering on each pixel of the image at each spatial resolution to extract the frequency band signal of interest; set the passband frequency of the bandpass filter according to the pulse frequency range of adults; if only the pulse frequency of a specific object needs to be extracted, a narrow passband is selected; the result of the above frequency domain bandpass filtering is the approximation of the change signal of interest; in the amplification of the pulse signal, it corresponds to the change of image color; 3) Video zoom Amplify the signal of the frequency band of interest; let I(x,t) represent the brightness value at any position x and any time t, and δ(t) is a displacement function, then, I(x,t)=f(x+δ(t)), I(x,0)=f(x) represents the brightness distribution corresponding to time 0; The purpose of the final video amplification is to synthesize the signal: Among them, α represents the magnification, By performing a first-order Taylor series expansion on the changed signal f(x+δ(t)) at time t, we can approximately obtain: Let B(x, t) represent the frequency domain bandpass filtering result of I(x, t) at each x position. Assuming that the motion signal of interest δ(t) is contained in the passband range, we have: Amplify the frequency domain filtering result B(x,t) by α times and add it back to the original signal, and we get: Combining equations (3), (4), and (5) yields: Then, according to the first-order Taylor series expansion, we can approximate: Thus, the spatial brightness signal δ(t) of the original signal f(x) is amplified by α times; 4) Output the enlarged pyramid layer image The amplified image frames are directly output for subsequent analysis, and the video reconstruction step of the traditional Gaussian video amplification method is abandoned, thereby improving the signal-to-noise ratio and calculation rate. The Laplacian pyramid video upscaling method is similar to the Gaussian pyramid process, with the following differences: 1) Laplacian pyramid decomposition is different from Gaussian pyramid decomposition in pulse signal amplification. Laplacian pyramid decomposition is based on Gaussian pyramid decomposition. The j-1 level approximate image is subjected to the inverse process of pyramid decomposition, that is, it is upsampled with the same step size as the pyramid decomposition and then passed through the same Gaussian low-pass filter to obtain the j-level image prediction residual. At this time, a layer of image of the Laplacian pyramid is obtained. The j-1 level image is used as the input image and the above steps are repeated to obtain a complete Laplacian pyramid. Each frame of the original video is subjected to spatial Laplacian pyramid decomposition as the baseband output; 2) Frequency domain filtering Laplace pyramid amplification is mainly used to extract respiratory motion, so the frequency band where human respiratory motion is located is used in frequency domain filtering; the passband frequency of the bandpass filter is set according to the respiratory frequency of an adult at rest; if only the respiratory frequency of a specific object needs to be extracted, a narrow passband is selected; Step 3: Determine the region of interest of the enlarged video image, including the region of interest of the pulse signal and the region of interest of the respiratory signal; Among them, the method for determining the region of interest of the pulse signal is: The video image amplified by the above method is used as input to quickly determine the region of interest where the pulse signal appears in the video. When selecting the region of interest of the pulse signal, a face recognition method based on a neural network is used to determine the head and face area of ​​the human body. For the image sequence in the head and face area, the frame difference method is used to calculate the difference Dg of the green channel signal intensity between frames. The green channel signal intensity of the nth frame and the first frame in the image sequence is recorded as g n and g1, the grayscale value of the corresponding pixel in the two frames is recorded as g n (u, v) and g1(u, v), subtract the grayscale values ​​of the pixels corresponding to the two frames of images according to the following formula to obtain the difference image D g (u,v), where u and v represent the column and row coordinates of the image pixel respectively: D g (u,v)=|g n (u,v)-g1(u,v)| (8) According to the average color (R m , G m , B m ), calculate the difference between the color of each pixel in the face area and the average color In the formula, R, G, and B represent the red channel signal, green channel signal, and blue channel signal of the image pixel respectively; R m , G m , B m are the average values ​​of the red, green and blue channel signals of all pixels in the face area; Calculate D for each pixel g / D c , which is a function of time; the maximum value of the function within the time of the video acquisition is obtained as the inter-frame change intensity of each pixel; according to the inter-frame change intensity of all pixels in the face area, they are sorted from large to small, and the ratio of a certain pixel ranking to the number of all pixels is recorded as λ, and the inter-frame change intensity of the pixel at the sorting position is set as the threshold value T of the inter-frame change intensity. λ When the value of λ is small, the region of interest of the pulse signal is small, which is suitable for conditions with high image resolution and sufficient illumination; when the value of λ is large, the region of interest of the pulse signal is large, which is suitable for conditions with low image resolution and insufficient illumination; for the case where the inter-frame change intensity is greater than the threshold value T λ The pixels of are the pixels with strong pulse signals; after closing operation on these pixels, a region with good continuity is generated; Step 4: Extract pulse and respiratory signals from the image sequences in the regions of interest of the pulse signal and the respiratory signal, respectively, and synchronize the pulse signal and the respiratory signal through cross-correlation operation; For pulse signal extraction, the red channel signal R, green channel signal G, and blue channel signal B are obtained by calculating the mean values ​​of the three color channel intensities in the region of interest of the pulse signal respectively; independent component analysis is performed on the three color channel signals. The principle is to assume that the collected signal is a linear superposition of several noises, and each signal and noise are independent of each other, which are called independent components, as shown in the following formula: In the formula, I R ,I G ,I B are the intensities of the red, green and blue channel signals respectively; C1, C2 and C3 are three independent components respectively; the three independent components C1, C2 and C3 are calculated from the three color channel signals. According to the above assumptions, one of the independent components must be a pulse-related signal, while the other two independent components are noise; the three independent component results are subjected to correlation analysis with the green channel signal G, and the result with the highest correlation is considered to be the pulse-related signal IPPG; then the IPPG signal is subjected to frequency domain analysis, and the frequency corresponding to the peak of the spectrum energy in 0.5-4Hz is selected as the pulse frequency fp; a narrowband bandpass filter is further constructed, and the passband frequency is fp±0.15Hz, and the IPPG signal is subjected to narrowband filtering to obtain the pulse signal Sp; Step 5: Perform cardiopulmonary coupling analysis on the synchronized breathing and pulse signals, including time domain analysis and frequency domain analysis, to achieve cardiopulmonary coupling measurement; The pulse and respiration are synchronized by performing a cross-correlation operation on the pulse signal Sp and the respiration signal Sr. The pulse signal Sp and the respiration signal Sr are synchronized by calculating the cross-correlation function of the pulse signal Sp and the respiration signal Sr within a respiration cycle Tr: In the formula, τ and t are any moments in a respiratory cycle Tr; If the cross-correlation function R pr (t) reaches the maximum value at time t = tm, then the pulse signal Sp(t) is time-shifted by tm, that is, the new pulse signal S′ p (t) = S p (tt m ) can achieve synchronization of cardiopulmonary signals.

2. The non-contact cardiopulmonary coupling evaluation method according to claim 1, characterized in that: In step 3, the respiratory signal region of interest is selected: the enlarged respiratory video image frame is binarized, each frame is subtracted from the first frame to generate a differential image, and the area with obvious changes in the subject's body movement is extracted from the differential image, and it is divided into the chest and abdomen and other areas. Finally, the chest and abdomen areas with obvious respiratory movement are selected as two regions of interest ROI respectively K and ROI F .

3. The non-contact cardiopulmonary coupling evaluation method according to claim 1, characterized in that: In step 4, the breathing video signal is extracted: Calculate two regions of interest (ROIs) of the respiratory signal in the ternary difference image respectively K and ROI F Differential image brightness and IROI within K and IROI F ; And select the frequency band corresponding to the respiratory signal of interest to calculate the difference image brightness and IROI of the region of interest K and IROI F After bandpass filtering, the respiratory signal curve R is extracted. K and R F , and the respiratory signal curve R K and R F Perform weighted averaging, setting the weight of the chest respiratory signal to 1, and the weight of the abdominal respiratory signal to the difference image brightness and IROI in the chest and abdominal regions of interest. K and IROI F The calculated signal intensity is the respiratory signal Sr.

4. The non-contact cardiopulmonary coupling evaluation method according to claim 1, characterized in that: In step 5, the cardiopulmonary coupling analysis method is: Cardiopulmonary coupling time domain analysis: In each cycle of the respiratory signal Sr, the Sr signal gradient is calculated, and the time when the gradient maximum value appears is recorded as the inhalation start time ti; a delay of a period of time ts is made based on ti, and the number of times hit when the pulse signal Sp maximum value appears in a period of time Δt before and after ts is calculated; the number of hits within 30 minutes is continuously recorded, and the number of hits is divided by the number of respirations within 30 minutes, and recorded as the heartbeat rate; the delay time ts is changed to obtain the heartbeat rate curve; the heartbeat time variability index is calculated according to the maximum and minimum values ​​in the curve, and the cardiopulmonary coupling strength is evaluated according to the heartbeat time variability index; Frequency domain analysis of cardiopulmonary coupling: Extract the peak positions of the pulse signal Sp and the respiratory signal Sr respectively, obtain the time difference between each peak, and form the pulse cycle sequence signal a(t) and the respiratory duration sequence signal b(t); calculate the cross-correlation function R of the two signals ab (t), and the cross-power spectral density function S of the two signals is calculated by the cross-correlation function ab (ω); Similarly, the autocorrelation function and autopower spectrum density R of the pulse cycle sequence signal a(t) and the breathing duration sequence signal b(t) are calculated respectively aa (t), R bb (t), S aa (ω), S bb (ω); then according to the cross power spectrum S ab (ω) and the auto-power spectrum S aa (ω), S bb (ω), calculate the coherence coefficient C ab (ω); finally, the coupling strength CPC(ω) of the two signals was calculated by the square of the cross-power spectrum and coherence of the two sequence signals; the cardiopulmonary coupling strength CPC(ω) was divided into different frequency bands, and the total power in the frequency range of 0-0.4 Hz, the low-frequency power in the frequency range of 0.04-0.1 Hz, the low-frequency power divided by the total power, the high-frequency power in the frequency range of 0.1-0.4 Hz, the high-frequency power divided by the total power, and the ratio of the low-frequency power to the high-frequency power were calculated respectively; the cardiopulmonary coupling strength was evaluated by the parameters of the above different frequency bands.

5. A cardiopulmonary coupling evaluation device according to the non-contact cardiopulmonary coupling evaluation method of claims 1 to 4, characterized in that: It includes a video image acquisition device, a video magnification device, a region of interest selection device, a video signal processing device, and a cardiopulmonary signal synchronization and analysis device; The video image acquisition device is a camera, which is used to capture the video of the head, face and upper body of the subject; The video magnification device comprises a pyramid decomposition unit, a time domain filtering unit and an amplification unit; wherein the pyramid decomposition unit amplifies the color change and the contour change in the video image by Gaussian video amplification method and Laplace video amplification method in parallel; the time domain filtering unit performs time domain band pass filtering on the brightness information of each layer of the image decomposed by the pyramid; outputs the frequency band information of interest of the pulse and the frequency band information of interest of the breathing; the amplification unit amplifies the frequency band information of interest of each layer of the pyramid, performs pyramid reconstruction on the amplified signal and adds it back to the original image, and outputs the amplified pulse video frame image and the breathing video frame image; The region of interest selection device selects regions of interest for the Gaussian video magnified image sequence and the Laplace video magnified image sequence, respectively, and outputs regions of interest for respiratory and pulse signals; The video signal processing device extracts the respiratory signal and the pulse signal from the respiratory and pulse signal regions of interest, respectively; The cardiopulmonary signal synchronization and analysis device is used to extract and synchronously process the pulse signal and the respiratory signal, and then perform time domain analysis and frequency domain analysis of the cardiopulmonary coupling.

Images