Obstructive sleep apnea syndrome detection based on ballistocardiogram via machine learning approach

1.   Introduction

As one of the most important life activities of human being, sleep is the best way to restore energy. Sleep apnea syndrome is a common sleep-related breathing disorder. According to the survey, about 4%-7% of men in the world suffer from sleep apnea syndrome, and women account for about 2%-4% ^[1,2,3]. Medical studies have shown that patients with sleep apnea syndrome have a higher risk of cardiovascular diseases, such as hypertension, headache and heart attack.

Polysomnography (PSG) is the golden standard for sleep staging and sleep apnea syndrome (SAS) detection. PSG monitors a variety of physiological parameters during sleep, including brain activity (EEG), eye movement (EOG), muscle activity or skeletal muscle activation (EMG) and heart rhythm (ECG). PSG monitoring needs to connect 30+ electrodes on the human body, which may seriously affect the experience and cause serious disturbance during sleep. At the same time, the identification of sleep stages and sleep apnea requires a professional doctor and it is costly, so it is necessary to find low-cost and reliable alternative methods to monitor sleep ^[4,5,6]. To make things worse, in many countries e.g. in China, medical resources for PSG monitoring are quite inadequate. If patients want to have sleep monitoring in the hospital, they usually need to make an appointment 1-2 months in advance. In recent years, many studies have studied to use simpler or different methods for sleep staging and OSA detection, which can expand the popularity of sleep health monitoring. As alternative methods of PSG, A. R. Hassan et al. researched single EEG based sleep staging for sleep quality evaluation to ensure wear ability and portability ^{[7,8,9,10,11,12]}. M. M. Rahman et al. proposed using single-channel EOG for sleep stage classification ^[13]. With respect to sleep apnea syndrome detection, many studies have also attempted to use electrodes or less signals, such as ECG ^{[14,15,16,17]}, SaO2 ^[18], EMG ^[19], nasal airflow ^[20,21,22], EEG ^[23], pulse ^[24] and other combined signals to achieve near precision sleep apnea detection. These signals are analyzed by time-domain information, frequency-domain information or combined information. In order to alleviate binding, A. R. Hassan et al. have researched single-channel ECG based OSA detection ^[25,26,27]. Though above mentioned OSA detection methods are much simpler than traditional PSG scheme, they still need to use electrode(s) to keep in touch with the human body for signal collection. Therefore, the problem of electrode dropping remains during the detection process, which may lead to unstable signal and affects the detection results. Considering that the above methods can't solve the problem of high cost, detection process may interfere with people's sleep, we focus on the study of non-contact intelligent sleep monitoring system, which uses pressure sensors to capture the pressure signals of breathing and heartbeat during sleep.

In this paper, we extract cardiac impulse and respiratory rhythm from a mixture of Ballistocardiogram (BCG) signals and respiratory signals. Sleep apnea related features are then obtained by extracting the time-domain and frequency-domain components of BCG and respiratory signals over fixed time intervals. Then machine learning based classification algorithm is used to detect OSA and model fusion technology is adopted to achieve higher recognition accuracy. The proposed sleep apnea detection method is low cost and will not disturb human's normal sleep.

2.   Methods

2.1.   BCG signal extraction

The mixed signal collected by the sensor is ${{S}}{}_{mixed}(t)$, which consists of heart beat signal ${S_{BCG}}(t)$, respiratory signal ${S_{rep}}(t)$ and noise signal $n(t)$, which can be written as follows:

The primary task is to extract BCG signals from mixed signals and then to obtain heart rate variability. Because the heart rate of different people is usually different, the difference of individual heart rate may affect the detection accuracy of heartbeat interval. So when extracting BCG signal, the average heart rate is calculated every half hour, and the filter parameters of band-pass filter are adjusted according to the average heart rate, so that noise and interference signals can be filtered and more accurate heartbeat interval sequence is obtained. We use $H(t)$ to represent the response function of the band-pass filter and obtain the BCG signal as:

For the analysis of non-linear and non-stationary signals, Norden E. Huang proposed a Hilbert-Huang Transform (HHT) method based on empirical model decomposition and Hilbert spectral analysis ^[24]. Due to that ${{S}}{}_{BCG}(t)$is a non-stationary signal, in this paper we use HHT method to extract BCG envelope and decompose it into the sum of multiple signals. The decomposition process is as follows:

(1) Set a set of ${r_0}(t)$ to satisfy that ${r_0}(t) = {S_{BCG}}(t)$.

(2) Let $U(t)$ be the set of local maximum values of ${r_0}(t)$, $L(t)$ be the set of local minimum value of ${r_0}(t)$ and ${m_1}(t)$ be the average of $U(t)$ and $L(t)$, that is:

(3) Obtain ${h_1}(t)$ by subtracting ${m_1}(t)$ from ${r_0}(t)$:

Then, it is judged that ${h_1}(t)$ is an Intrinsic Mode Function (IMF), that is, it is the signal of all local maximum points or local minimum points. If so, proceed with step (4), otherwise repeat the above process until it becomes and IMF and we record it as ${c_1}(t)$.

(4) Obtain the residual signal, say ${r_1}(t)$, by subtracting ${c_1}(t)$ from the original signal ${r_0}(t)$.

(5) Repeat steps (1) to (3) until the differential signal ${r_n}(t)$ becomes a monotonic function. Then the original signal $x(t)$ is expressed as the sum of all IMF and ${r_n}(t)$, and the expression is as follows:

In the multi-level decomposed signals ${c_i}(t)$, we choose the most typical signal which can reflect the characteristics of heartbeat interval to calculate the envelope of ${{S}}{}_{BCG}(t)$. As shown in Figure 1, since the J peak of the BCG signal is the maximum point of the signal in a cycle, the peak value of envelope signal is the the same as that of BCG signal.

Figure 1.  Typical BCG signal waveform.

DownLoad: Full-Size Img PowerPoint

After obtaining the envelope of BCG signal, peak detection algorithm is used to detect the peak value of envelope signal. The heartbeat interval of BCG signal is the difference sequence of the peak value of the envelope signal.

As shown in the same figure of Figure 1, it is likely that multiple peaks are detected in the same heartbeat cycle (e.g. J and H peaks), which will lead to inaccurate BCG cycle detection. It may also occur that no J peak can be detected in a BCG cycle. Therefore, the time interval sequence should be calibrated according to the characteristics of heartbeat interval. Assume that the J-J interval sequence is ${{J}}{{{J}}_{{i}}}$, ${{i = 1, 2, }}...{{, K}}$, where ${{K}}$ is the total number of J-J intervals in a given period of time. The calibration algorithm is as follows:

(1) Calculate the mean value of J-J interval sequence in a determined time duration, and express it as$\overline {{{J}}{{{J}}_{{j}}}}$.

(2) For any J-J interval ${{J}}{{{J}}_{{i}}}$, if ${{J}}{{{J}}_{{i}}} > 1.5*\overline {{{J}}{{{J}}_{{j}}}}$, then it is split into two parts of ${{J}}{{{J}}_{{{i - }}1}}$ and ${{J}}{{{J}}_{{{i + }}1}}$, and we can get the calibrated J-J intervals, that are:

(3) For any J-J interval ${{J}}{{{J}}_{{i}}}$, if ${{J}}{{{J}}_{{i}}} < 0.5*\overline {{{J}}{{{J}}_{{j}}}}$, then it is split into two parts and each is merged with the previous J-J interval, and we can get the calibrated J-J intervals as follows:

(4) Repeat the process until all data is processed.

Figure 2 shows the BCG signal envelope extracted by the improved algorithm. It can be seen that the J peak of BCG signal can be accurately detected, so that the accurate sequence of heartbeat interval can be obtained.

Figure 2.  BCG signal and its envelope.

DownLoad: Full-Size Img PowerPoint

2.2.   Feature Extraction

In this part, we extract the features of BCG signals, including time-domain features, frequency-domain features and similarity features.

2.2.1.   Time-domain features

Heart rate variability (HRV) reflects changes in heart rate and contains information about neurohumoral factors regulating cardiovascular function. Medical studies of heart rate variability have found that statistical characteristics of heart rate cycles during sleep apnea can be used to evaluate sleep apnea. Through the analysis of mean square deviation, standard deviation and coefficient of variation the characteristics of heartbeat interval and time domain change can be obtained. The indices derived from the difference of heartbeat interval include SDNN (standard deviation of NN interval), SDSD (standard deviation of difference between adjacent NN intervals), RMSSD (square root of the mean of the sum of the squares of differences between adjacent NN intervals), NN50 (number of pairs of adjacent NN intervals difference by more than 50 ms in the entire recording), PNN50 (NN50 count divided by the total number of all NN intervals) and CV (heart rate variability coefficient), as shown in Table 1.

Table 1.  Time-domain features.

Items	Explanation
Mean	Mean =${\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 {{N}}}}\right.} \!\lower0.7ex\hbox{${{N}}$}}\sum\nolimits_{i = 1}^N {J{J_i}}$
Variance	Variance =${\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 {{N}}}}\right.} \!\lower0.7ex\hbox{${{N}}$}}\sum\nolimits_{i = 1}^N {{{\left({J{J_i}{{ - }}\overline {{{J}}{{{J}}_{{j}}}} } \right)}^2}}$
Max	Max = max{$J{J_i}$}, $1 \leqslant {{i}} \leqslant N$
Min	Min = min{$J{J_i}$}, $1 \leqslant {{i}} \leqslant N$
SDNN	SDNN =$\sqrt {{\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 {{N}}}}\right.} \!\lower0.7ex\hbox{${{N}}$}}\sum\nolimits_{i = 1}^N {{{\left({J{J_i}{{ - }}\overline {{{J}}{{{J}}_{{j}}}} } \right)}^2}} }$
SDSD	SDSD = $\sqrt {{\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 {{N}}}}\right.} \!\lower0.7ex\hbox{${{N}}$}}\sum\nolimits_{i = 2}^N {{{\left({J{J_i}{{ - }}J{J_{i{{ - }}1}}} \right)}^2}} }$
RMSSD	RMSSD = $\sqrt {{\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 {{N}}}}\right.} \!\lower0.7ex\hbox{${{N}}$}}{{\sum\nolimits_{i = 1}^N {\left[{{{\left({J{J_i}{{ - }}J{J_{i{{ - }}1}}} \right)}^2}{{ - }}\overline {J{J_i}{{ - }}J{J_{i{{ - }}1}}} } \right]} }^2}}$
PNN50	PNN50 =$\frac{{{{Count}}\left\{ {J{J_i} \geqslant 50{{ms}}} \right\}}}{{{{Count}}\left\{ {J{J_i}} \right\}}}$
CV	CV = SDNN/Mean

---

 | Show Table

DownLoad: CSV

2.2.2.   Frequency-domain features

Studies have shown that heart rate and heart rate variability vary with different sleep stages, and heart rate variability is affected by sleep apnea. In the appropriate frequency range, spectrum analysis can reflect the periodic change of heart rate. The pattern of bradycardia and tachycardia during apnea is attributed to effective parasympathetic control of heart rate during sleep ^[31]. Intermittent apnea is accompanied by sympathetic nerve termination or awakening. Spectral analysis of HRV can be used to establish and quantitatively evaluate the activation sympathetic and parasympathetic nerves to heart beat ^[29]. According to the physiological interpretation of very low frequency (VLF) components, low frequency (LF) components (0.04-0.15 Hz) reflect the sympathetic control of blood pressure baroreflex, and high frequency (HF) components (0.15-0.4 Hz) reflect the respiratory rhythm, which is considered to be related to the parasympathetic control of the heart ^[30]. Overall, total power assesses overall heart rate variability; LF represents sympathetic and parasympathetic activity, HF represents parasympathetic activity, LF/(LF+HF) quantitatively measures the sympathetic activity, HF/(LF+HF) represents parasympathetic activity, and LF/HF represents the activity balance of autonomic nerves. Table 2 lists the frequency domain characteristics used.

Table 2.  Frequency-domain features.

Items	Explanation
TP	The total power in 0-0.4Hz band
vLF	The power in 0-0.04Hz band
LF	The power in 0.04-0.15Hz band
HF	The power in 0.15-0.4Hz band
nLF	The normalized power in LF band, i.e.LF/(TP-vLF)
nHF	The normalized power in HF band, i.e.HF/(TP-vLF)
LF/HF	The ratio of power in LF and HF band, i.e. LF/HF

---

 | Show Table

DownLoad: CSV

2.2.3.   Similarity features

Echocardiographic signals may be affected by noise. If the signal amplitude within a fixed time (e.g. 1 second) is taken as the characteristic of heart rate variability, it can't accurately reflect the characteristics of echocardiographic signal. At the same time, if two signal segments are separated at fixed intervals, even if they come from the same sleep apnea stage, they may have low similarities. In order to determine the change characteristics for the signals in a fixed time and make them have finer granularity and reduce the complexity of the algorithm, a vector composed of normalized time-frequency features is used to identify the signal characteristics in a period of time. The time-domain and frequency-domain features are extracted to identify the signals within the time interval, and the similarity between the time interval and the similarity of 20 time intervals before and after the time interval is calculated as a measure of the similarity characteristics of sleep apnea. The similarity of any two time periods is measured by cosine similarity, Euclidean distance and linear combination of features.

The cosine similarity coefficient is expressed as:

where $partA$ and $partB$ are vectors representing the time-domain and frequency-domain features of two time intervals.

The Euclidean distance similarity coefficient is:

Euclidean distance similarity coefficient measures the amplitude similarity of two periods. From above formula, it can be seen that the smaller the similarity coefficient, the higher the similarity. Vice versa, the smaller the similarity coefficient, the lower the similarity.

In addition to cosine similarity coefficient and Euclidean similarity coefficient, the non-linear combination of time-domain and frequency-domain features is introduced to enhance the expressive ability of the model and better reflect the characteristics of each fixed time period. The linear combination relation can be expressed as:

2.3.   Decision tree based OSA detection model

It is a binary classification problem to divide BCG signals into OSA and non-OSA segments. This paper firstly uses decision tree classification model to evaluate and then uses model fusion technology to improve the evaluation results. As shown in Figure 3, the decision tree model is used as an OSA recognition model. According to the characteristics of sleep apnea syndrome, the sample points in the sample space of sleep apnea syndrome are divided effectively, and the classification of each subspace sample is finally obtained.

Figure 3.  Decision tree based classification model.

DownLoad: Full-Size Img PowerPoint

The training data set of sleep apnea is defined as $D{{ = \{ (}}{{{x}}_1}{{, }}{{{y}}_1}{{), (}}{{{x}}_2}{{, }}{{{y}}_2}{{), }}...{{, (}}{{{x}}_m}{{, }}{{{y}}_m}{{)\} }}$ and the feature set is $A = \{ {a_1}, {a_2}, ..., {a_d}\}$, the decision tree algorithm is designed as follows:

Step 1: Generate decision tree nodes;

Step 2: If the sleep apnea sample training set belongs to the same OSA category or non-OSA category, return to step 1, otherwise continue Step 3;

Step 3: If the training set completes traversing feature A (such as LF) or the sleep apnea sample training set belongs to the same OSA category or non-OSA category on feature A, is divided into the category with the largest number of samples in the OSA training set and return to step 1, otherwise it enters Step 4;

Step 4: Select the best segmentation feature, say ${a_i}$, from the sample feature set A of sleep apnea and traverse every value of ${a_i}$ to generate a branch of OSA classification. If the subset of $D$ in ${a_i}$ is empty, the branch nodes of OSA classification are marked as leaf nodes, whose categories are aligned with the parent nodes with the most samples of sleep apnea. Otherwise, traverse other attributes in A.

The purpose of selecting the optimal OSA partition feature is to make the sleep apnea samples contained in the branch nodes of decision tree belong to the same category as possible, that is, the nodes have high "purity". Here, the ID3 algorithm is used to select the best partitioning feature, that is:

In this way, the attributes with larger information gain will be closer to the root node, and the decision tree will select the attributes with larger information gain to partition.

2.4.   Model fusion

Model fusion can improve the classification performance compared with single model by by constructing a stronger classifier with several weak classifier. Besides decision tree classifier, logistic regression and support vector machine classifiers are integrated with stacking method. The implementation of model fusion is shown in Figure 4, that is:

Figure 4.  Schematic diagram of model fusion.

DownLoad: Full-Size Img PowerPoint

(1) The data set is divided into $K$ parts and $n$ models ($n = 3$ in this paper) are used in the training process;

(2) In the training phase, for each of the $n$ models, $K - 1$ parts data are set as the training set and the remaining part is set as the testing set. When the parameters of the model are adjusted by the training set, the testing set is input and the predicted output is taken as a new feature. In the testing stage, the predicted mean value of the model is used as new features.

(3) This is done k times, and the entire data set gets $n$ New Features. After $K$ times operation, we get $n$ new features on the whole data set.

(4) The new features and labels are trained as input of the new classifier, and then the feature of the testing set is input into the model to obtain the final prediction results.

3.   Results and discussion

3.1.   Evaluation method

The predicted results of the sleep apnea detection model fall into four cases: True Positive (TP), False Positive (FP), True Negative (TN) and False Negative (FN). Let ${{TP}}$, $FP$, $TN$ and $FN$ represent the number of samples of their corresponding categories respectively. Obviously, the total number of samples is equal to ${{TP}}$+$FP$+$TN$+$FN$.

The metrics used in the classification models include sensitivity, specificity and accuracy, which can be calculated from the confusion matrix. In the sleep apnea detection model, sensitivity indicates the proportion of actual sleep apnea samples detected to total sleep apnea samples; specificity indicates the proportion of actual non-sleep apnea samples detected to all non-sleep apnea samples; accuracy indicates the proportion of actual sleep apnea samples detected plus actual non-sleep apnea samples detected to the total number of samples.

Sensitivity and specificity are contradictory. When the sensitivity is high, the specificity of the model is often low. On the contrary, when the specificity is high, the sensitivity of the model is low. Therefore, the tradeoff between sensitivity and specificity should be considered when adjusting the parameters of the classification model.

3.2.   Result analysis

Firstly, the performance for different length of time period is analyzed. As shown in Figure 5, it has higher precision and sensitivity for 20-second segmentation period, but it has the lowest specificity. The 60-second segmentation period has higher specificity, but its precision and sensitivity are low. A longer interval is beneficial to the analysis of heart rate variability and the detection of sleep apnea syndrome. According to the analysis results, in the following analysis we choose 50 seconds as the segmentation interval to study the detection of sleep apnea.

Figure 5.  Performance under different threshold.

DownLoad: Full-Size Img PowerPoint

Figure 6 shows the relationship among sensitivity, specificity and accuracy under different thresholds. The results showed that with the increase of threshold, the sensitivity and accuracy increased, while the specificity decreased. When the threshold approached 0.4, the specificity decreased significantly, which indicated that the model was easier to distinguish between the segments of sleep apnea and those of non-sleep apnea. According to the results of the graph, 0.3 was chosen as the threshold value for detecting sleep apnea in the following analysis, taking sensitivity, specificity and accuracy into account.

Figure 6.  Performance under different threshold.

DownLoad: Full-Size Img PowerPoint

Based on the training results of the above classification model, stacking model fusion method was used to train the sleep apnea samples. Figure 7 shows the test results after model fusion. Under the same tag serial number, the left column is the result without model fusion, while the right column is the result of model fusion. Obviously, the sensitivity, specificity and accuracy all have been improved after model fusion.

Figure 7.  Performance with model fusion.

DownLoad: Full-Size Img PowerPoint

We compare the four metrics with other methods proposed in the literature and the result is shown in Table 3. As we can see from the table, our method of achieves higher accuracy, specificity and sensitivity while lower false alarm rate than other methods. The advantages of our algorithm are better reflected in the mixed data set, because the mixed data set has many different shapes, so it is more sensitive to the shape.

Table 3.  Metrics comparison with other methods.

Method	Accuracy (%)	Sensitivity (%)	Specificity (%)
Proposed	75.0	74.0	74.0
C-SVM with linear kernel C=0.5 [14]	89.83	88.35	90.72
PSD and SampEn [16]	72.9	72.2	73.3
LD +RR [17]	85.8	76.9	91.3
ECG qualityexclusion criteria [24]	73.33	62.5	85.71
TQWT and RUSBoost [25]	88.88	87.58	91.49
DT-CWT and LogitBoost [26]	84.4	90.38	74.84
TQWT [27]	87.33	81.99	90.72

---

 | Show Table

DownLoad: CSV

4.   Conclusion

In this paper, a sleep apnea detection algorithm based on electrocardiogram is studied. The detection of sleep apnea is transformed into a classification problem by extracting the physiological signal features collected in a fixed time, namely sleep apnea and non sleep apnea. The time-domain features, frequency-domain features and similarity feature extracted from the original signal are used as the input of machine learning model. The sensitivity, specificity and accuracy of the single model were 68%, 73% and 71%, respectively. After model fusion, sensitivity, specificity and accuracy were increased to 74%, 75% and 75% respectively. The sleep apnea detection model proposed in this paper has the characteristics of moderate complexity, high spatial complexity and high sensitivity. It can be used for sleep apnea screening and sleep health detection at home, so as to reduce the risk of cardiovascular diseases. In the future, we will consider studying new methods, e.g. deep learning, to further improve the detection performance.

Acknowledgments

This work is supported by National Key R & D Program of China under grant number SQ2018YFC200148-03.

Conflict of interest

The authors have declared that no competing interests exist.