CN103584840B - 基于脑电、心率变异及两者相干性的自动睡眠分期方法 - Google Patents

Abstract

本发明公开了一种基于脑电、心率变异及两者相干性的自动睡眠分期方法，所述方法包括以下步骤：提取被测试者的脑电信号和心率变异信号；从脑电信号、心率变异信号中共提取9个特征向量；对9个特征向量进行主成分分析；特征提取后，使用支持向量机分类器对特征进行识别，进行自动睡眠分期。该方法可有效地提高自动睡眠分期系统准确性和简便性，并获得可观的社会效益和经济效益。最佳实施方案拟采用专利转让、技术合作或产品开发。由于该技术操作简单，敏感性强，基于该技术开发的产品可应用于睡眠监测、睡眠科研等多种场景等。

Description

基于脑电、心率变异及两者相干性的自动睡眠分期方法

技术领域

本发明涉及自动睡眠分期领域，特别涉及一种基于脑电、心率变异及两者相干性的自动睡眠分期方法。

背景技术

睡眠是人体重要的生理活动，人类生命约有1/3的时间在睡眠中度过。睡眠是机体自我修复和自我恢复的过程，而良好的睡眠是维持身心健康的基本需要。进行睡眠分期是探讨睡眠期间生理调节变化规律、进行睡眠质量评价的基础和前提，更重要的是，作为睡眠监护仪最核心的技术模块，睡眠分期为针对睡眠的多模式生理学研究提供了一种可靠的量化指标。

睡眠分期就是根据睡眠期间人体生理信号的不同将睡眠分为不同的阶段。主要分为非快速眼动期（NREM）和快速眼动期(REM)。非快速眼动期从夜间入睡开始，又细分为4期，NREM-Ⅰ为入睡期，NREM-Ⅱ为浅睡期，NREM-Ⅲ为中度睡眠期，NREM-Ⅳ为深度睡眠期。在这个阶段人的呼吸变浅、变慢，心率变慢、血压下降，全身肌肉松弛，无明显的眼球运动。快速眼动期的特征是眼球快速转动。人体的感觉功能进一步减退，肌肉也更加松弛，肌腱反射消失。这个阶段，体内各种代谢功能都显著增加，以保证大脑组织蛋白的合成和消耗物质的补充，使神经系统正常发育，并为第二天的活动积蓄能量。图1是正常人的睡眠周期图。

目前，基于各种生理信号的睡眠监测方法主要有多导睡眠图（PSG）、EEG（脑电信号）监测法、ECG（心电信号）监测法、血氧饱和度监测法和身体活动变化记录法等。（1）PSG是进行睡眠医学研究的经典方法，也称金方法。它同时记录EEG、EOG（眼电信号）和EMG（肌电信号），在此基础上根据R&K规则^[1]，进行睡眠状态和睡眠阶段的判定，但此方法需要佩戴至少十枚电极，被试者经常感觉不适，对正常睡眠干扰较大。而且一般采用对脑电波形进行人工判读的方法，但这需要很大的工作量，且有赖于判读者的个人经验，往往夹带有主观因素，缺乏一致性客观标准。（2）EEG监测法、ECG监测法、血氧饱和度监测法和身体活动变化记录法等多采用单一参数或将各生理参数简单融合作为自动睡眠分期的指标，正确率较低。

发明内容

本发明提供了一种基于脑电、心率变异及两者相干性的自动睡眠分期方法，本发明提高了睡眠的质量，且提高了自动睡眠分期的质量，详见下文描述：

一种基于脑电、心率变异及两者相干性的自动睡眠分期方法，所述方法包括以下步骤：

（1）提取被测试者的脑电信号和心率变异信号；

（2）从脑电信号、心率变异信号中共提取9个特征向量；

（3）对9个特征向量进行主成分分析；

（4）特征降维后，使用支持向量机分类器对特征进行识别，进行自动睡眠分期。

所述从脑电信号、心率变异信号中共提取9个特征向量的操作具体为：

1）通过对脑电信号进行小波变换、希尔伯特黄变换和奇异值分解获取脑电信号的小波熵、Hilbert-Huang熵，以及奇异值第一主成分，分别记为特征向量P₁,P₂,P₃；

2）计算极低频VLF、低频LF和高频HF的频谱熵，并通过小波变换计算心率变异信号的分形维数，分别记为特征向量P₄,P₅,P₆,P₇；

3）对脑电信号以及心率变异信号进行功率谱密度分析，并进行归一化处理，然后分别计算脑电信号delta频段与心率变异信号LF、HF的相干系数，分别记为特征向量P₈,P₉。

所述对脑电信号以及心率变异信号进行功率谱密度分析，并进行归一化处理，然后分别计算脑电信号delta频段与心率变异信号LF、HF的相干系数，分别记为特征向量P₈,P₉的操作具体为：

1）截取同一时间段的EEG、HRV信号，进行预处理，包括变参考、降采样、去噪声与干扰，预处理后的EEG信号记为x，HRV信号记为y；

2）采用Welch算法计算EEG、HRV的功率谱密度以及两者的互功率谱密度，分别记为Pxx(f₁)，Pyy(f₂)，Pxy(f₁,f₂)；

3）用改进的相干函数计算EEG信号delta频段与HRV信号LF、HF的相干系数。

所述对9个特征向量进行主成分分析的操作具体为：

（1）对原始数据P_n*9进行标准化处理，矩阵中的元素减去所在列的均值，然后除以所在列的标准差，使得每个变量的均值为0，方差为1，得到矩阵P_n*9 ^*；

（2）然后求其协方差矩阵C_9*9,P_n*9 ^*中任两列之间可以计算两变量间的协方差；

（3）对协方差矩阵C_9*9进行特征根分解，得到特征根矩阵及特征向量；

（4）求原始数据P_n*9在新的矢量空间中的投影，即主成分向量组F_n*9：

（5）求累积贡献率；

（6）选定预设的累积贡献率，使前d个主成分F_n*d作为新的数据进行模式识别。

本发明提供的技术方案的有益效果是：本发明通过计算EEG信号delta频段与HRV信号LF、HF的相干系数，与提取的EEG、HRV多种频域特征和非线性特征融合作为睡眠分期的特征矩阵，PCA去冗余信息后作为后续模式识别的输入，从而准确、客观的进行自动睡眠分期。该项发明可有效地提高睡眠分期的准确性和简便性，并获得可观的社会效益和经济效益；为睡眠监护仪等提供了最核心的技术支持，给实际应用带来便捷，并且可应用于多种工作场景。

附图说明

图1为睡眠周期图；

图2为自动睡眠分期方法流程图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。

睡眠作为一种生理状态，必然要在身体的多个部位及多个生理过程中反映出来。EEG是中枢神经信号，是睡眠分期最敏感最准确的生理信号，心电信号是自主神经信号，提取的心率变异信号（HRV）反映交感神经系统和副交感神经系统的活动水平，而自主神经系统的调控能力随着睡眠程度的不同而改变。有研究表明，自主神经活动与睡眠在皮层下有共同的调节中枢，心电信号作为主要的自主神经信号，与脑电信号的有效融合，势必能够提高自动睡眠分期的正确率。

101：提取被测试者的脑电信号和心率变异信号；

EEG包含大量的生理信息，经过分析处理后，可对睡眠进行分期。相对于其他生理参数，EEG反映睡眠各期特征最为明显，是目前睡眠分期的“金标准”，所以近年来通过脑电信号自动睡眠分期的准确性和客观性研究一直受到人们的关注，研究也最为成熟。EEG通常包含5种基本节律：delta波（0.5-4Hz），theta波（4-8Hz），aleph波（8-13Hz），beta波（13-30Hz），gamma波（30-50Hz）。

心率变异即逐次窦性心跳间期之间的微小涨落，直接受自主神经系统中的交感神经和副交感神经的双重支配。心率变异性是自主神经系统功能最直接、有效的反映，而自主神经系统的调控能力随着睡眠程度的不同而改变，研究表明HRV是进行自动睡眠分期的很好指标。HRV典型的频谱可有三个峰，大致位于0.04Hz以下，0.05-0.15Hz和>0.15Hz，分别称为极低频(VLF)、低频(LF)和高频(HF)峰。

在本方法中，由于该两个参数测量方法成熟，与其他生理参数相比能更准确、敏感的反应各睡眠阶段，并且心脑电信号在不同的睡眠阶段相关程度不同，所以采用这两种生理参数有效融合来进行自动睡眠分期在实际应用中具有优越性。

102：从脑电信号、心率变异信号中共提取9个特征向量；

1）EEG特征提取：通过对脑电信号进行小波变换、希尔伯特黄(Hilbert-Huang)变换和奇异值分解获取脑电信号的小波熵^[2]、Hilbert-Huang熵^[3]，以及奇异值第一主成分（奇异谱中的最大值）^[4]，分别记为特征向量P₁,P₂,P₃。

2）HRV特征提取：计算极低频（VLF）、低频（LF）和高频（HF）的频谱熵，并通过小波变换计算HRV的分形维数^[5][6]，分别记为特征向量P₄,P₅,P₆,P₇。

3）EEG信号delta频段与HRV参数的相干系数：对EEG信号以及HRV信号进行功率谱密度分析，并进行归一化处理，然后利用改进的相干函数分别计算EEG信号delta频段与HRV信号LF、HF的相干系数，分别记为特征向量P₈,P₉。

相干分析是在频域上描述两个信号相关程度的实值函数。频域相干性是对两个信号在不同频率下的相位一致性的量度，当两者的某一频率成分fi相位一致时,

趋于1，可以用它来判断两个信号是否在某一频率发生了固定相位的振荡。但是该相干系数评价参数存在一个很大的局限性就在于它只能计算某一相同频段的相干系数。为了克服这种局限性，本发明采用改进的相干性评价参数Coh_xy(f₁f₂)，其定义为：

${\mathrm{Coh}}_{\mathrm{xy}}\left(f_1{f}_2\right)=\frac{{\left|\mathrm{Pxy}\left(f_1{f}_2\right)\right|}^2}{\mathrm{Pxx}\left(f_1\right)\mathrm{Pyy}\left(f_2\right)}---(1-1)$

它表示脑电信号x在频率分量为f1，心率变异信号y在频率分量为f2处的分量振幅乘积的标准化均值，其取值区间为[0,1]，反映了2个信号的相关程度。相干谱越接近1，表明2个信号越相关。相干系数为1，表明2个信号之间高度相关，且一个信号是另一个信号的倍数；相干系数为0，表明2个信号完全无关。

本发明中，求EEG信号delta频段与HRV信号低频段（LF）、高频段（HF）的相干系数步骤如下：

1）在睡眠分期标准中，时间分辨率为30s，与其一致，截取同一时间段30s的EEG、HRV信号，进行预处理，包括变参考、降采样、去噪声与干扰，预处理后的EEG信号记为x，HRV信号记为y。

2）采用Welch算法计算EEG、HRV的功率谱密度以及两者的互功率谱密度，分别记为Pxx(f₁)，Pyy(f₂)，Pxy(f₁,f₂)。

$\mathrm{Pxx}\left(f_1\right)=\frac{1}{\mathrm{MU}}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\left|\underset{n=0}{\overset{M-1}{\mathrm{\Σ}}}{x}_i\left(n\right)d_2\left(n\right)e^{-j2\mathrm{\π}{f}_{1}n}\right|}^2---(1-2)$

$\mathrm{Pyy}\left(f_2\right)=\frac{1}{\mathrm{MU}}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\left|\underset{n=0}{\overset{M-1}{\mathrm{\Σ}}}{y}_i\left(n\right)d_2\left(n\right)e^{-j2\mathrm{\π}{f}_{2}n}\right|}^2---(1-3)$

$\mathrm{Pxy}\left(f_1{f}_2\right)=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{X}_i\left(f_1\right)Y_i\left(f_2\right)---(1-4)$

其中U是归一化因子，d₂(n)是高斯窗函数，L为将该段数据所分段数，M为每段的长度。x_t(n)是x(EEG信号)的第i段数据，y_t(n)是y(HRV信号)的第i段数据。 $X_i\left(f_1\right)=\underset{-\∞}{\overset{\∞}{\mathrm{\Σ}}}{x}_i\left(n\right)e^{-j2\mathrm{\π}{f}_{1}n},Y_i\left(f_2\right)=\underset{-\∞}{\overset{\∞}{\mathrm{\Σ}}}{y}_i\left(n\right)e^{-j2\mathrm{\π}{\mathrm{fn}}_2},$ j为虚单位， $j=\sqrt{\text{-1}}.$

3）用改进的相干函数计算EEG信号delta频段与HRV信号LF、HF的相干系数

${\mathrm{Coh}}_{\mathrm{xy}}\left(f_1{f}_2\right)=\frac{{\left|\mathrm{Pxy}\left(f_1{f}_2\right)\right|}^2}{\mathrm{Pxx}\left(f_1\right)\mathrm{Pyy}\left(f_2\right)}---(1-5)$

然后计算EEG信号的f₁在delta频段[0.5,4]，HRV信号的f₂在LF[0.05,0.15]范围内的平均相干系数Coh_xyLF，作为EEG信号delta频段与HRV信号LF的相干系数，再计算f₁在delta频段[0.5,4]，f₂在HF[0.15,0.4]范围的平均相干系数Coh_xyHF，作为delta频段与HF的相干系数。

${\mathrm{Coh}}_{\mathrm{xyLF}}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{Coh}}_{\mathrm{xy}}\left(f_1{f}_2\right),0.5\≤f_1\≤\mathrm{4,0.05}\≤f_2\≤0.15---(1-6)$

${\mathrm{Coh}}_{\mathrm{xyLF}}=\frac{1}{m}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{Coh}}_{\mathrm{xy}}\left(f_1{f}_2\right),0.5\≤f_1\≤\mathrm{4,0.15}\≤f_2\≤0.4---(1-7)$

其中，n为0.5≤f₁≤4,0.05≤f₂≤0.15范围内Coh_xy(f₁f₂)的点数，m为0.5≤f₁≤4,0.15≤f₂≤0.4范围内Coh_xy(f₁f₂)的点数。

103：对9个特征向量进行主成分分析（PCA）；

各个参数所含的信息之间具有一定的重叠性和相关性，若直接将它们用于模式识别，会造成模型参数的过度拟合而降低分类的准确性和可靠性，且会因为数据量过大而降低分类的速度。因此在模式分类之前，本发明首先利用PCA对得到的特征向量进行降维处理。

PCA根据方差最大化原理，用一组线性无关且相互正交的新向量表征原来的数据矩阵的行（或列），达到压缩变量个数，剔除冗余信息，最大化保存有效信息的目的。原始向量组为（P₁，P₂，…，P₉），主成分向量组记为（F₁，F₂，…,F_m），通常m小于9。则主成分与原始向量组的关系为：

$\left\{\begin{array}{c}{F}_1=a_{11}{P}_1+a_{12}{P}_2+...+a_{19}{P}_9\\ F_2=a_{21}{P}_1+a_{22}{P}_2+...+a_{29}{P}_9\\ ...\\ F_k=a_{k1}{P}_1+a_{k2}{P}_2+...+a_{k9}{P}_9\\ ...\\ \mathrm{Fm}=a_{m1}{P}_1+a_{m2}{P}_2+...+a_{m9}{P}_9\end{array}\right.---(1-8)$

其中F₁蕴含信息量最多，具有最大方差，称为第一主成分，F₂，…，F_m依次递减，称为第二主成分、…、第m主成分。因此主成分分析的过程可以看作是确定权重系数a_ik(i=1,…,m;k=1,…9)的过程。

在本发明中，首先对该9个变量进行n次观测，得到的观测数据可用下面的矩阵表示

$P_{n*9}=\left[\begin{array}{cccccc}{P}_{11}& P_{12}& ...& P_{1h}& ...& P_{19}\\ P_{21}& P_{22}& ...& P_{2h}& ...& P_{29}\\ .& .& & .& & .\\ .& .& ...& .& ...& .\\ .& .& & .& & .\\ P_{b1}& P_{b2}& ...& P_{\mathrm{bh}}& ...& P_{b9}\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ P_{n1}& P_{n2}& ...& P_{\mathrm{nh}}& ...& P_{n9}\end{array}\right]---(1-9)$

其中，P_bh为第b次观测（即第b个样本）的第h个特征。

用PCA进行特征降维的求解过程如下：

（1）对原始数据P_n*9进行标准化处理，矩阵中的元素减去所在列的均值，然后除以所在列的标准差，使得每个变量的均值为0，方差为1，得到矩阵P_n*9 ^*。

P_n*9 ^*=[y_bh]_n*9，b=1,2,…,n;h=1,2,…,9     （1-10）

$y_{\mathrm{bh}}=\[P_{\mathrm{bh}}-\stackrel{\‾}{P_h}\]/S_h---(1-11)$

其中 $\stackrel{\‾}{P_h}=\frac{1}{n}\underset{b=1}{\overset{n}{\mathrm{\Σ}}}{P}_{\mathrm{bh}},S_h=\sqrt{\frac{1}{n-1}\underset{b=1}{\overset{n}{\mathrm{\Σ}}}{(P_{\mathrm{bh}}-\stackrel{\‾}{P_h})}^2}.$

（2）然后求其协方差矩阵C_9*9,P_n*9 ^*中任两列之间可以计算两变量间的协方差，于是得到协方差矩阵：

$C_{9*9}=\left[\begin{array}{cccc}{S_1}^2& \mathrm{cov}\left(\mathrm{1,2}\right)& ...& \mathrm{cov}\left(\mathrm{1,9}\right)\\ \mathrm{cov}\left(\mathrm{2,1}\right)& {S_2}^2& ...& \mathrm{cov}\left(\mathrm{2,9}\right)\\ .& .& & .\\ .& .& ...& .\\ .& .& & .\\ \mathrm{cov}\left(\mathrm{9,1}\right)& \mathrm{cov}\left(\mathrm{9,2}\right)& ...& {S_9}^2\end{array}\right]---(1-12)$

（3）对协方差矩阵C₉*₉进行特征根分解，得到特征根矩阵Λ_9*9及特征向量U_9*9。

C_9*9=U_9*9Λ_9*9U_9*9′     （1-13）

其中特征向量U_9*9作为主成分的坐标轴，构成新的矢量空间，

其中特征根λ_r(r=1,2,…9)的大小代表第r个主成分蕴含的信息量。U_9*9′是U_9*9的转置矩阵。

（4）求原始数据P_n*9在新的矢量空间中的投影，即主成分向量组F_n*9：

F_n*9=P_n*9U_9*9     （1-14）

（5）求累积贡献率。每个主成分的特征根大小代表其蕴含信息量的多少。求前k（k=1,…,9）个主成分的累积贡献率。

${\mathrm{pre}}_k=\frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\λ}}_i}{\underset{i=1}{\overset{9}{\mathrm{\Σ}}}{\mathrm{\λ}}_i}---(1-15)$

其中，λ_i是求出的第i个特征根。

（6）选定预设的累积贡献率，使前d个主成分F_n*d作为新的数据进行模式识别。

例如：一共得出8个主成分。第一个主成分F1的贡献率为48%，F2的贡献率为32%，F3的贡献率为15%，F4,F5,F6,F7,F8总共的贡献率为5%（8个主成分的贡献率一共是100%）。那么前三个主成分（F1,F2,F3）的累积贡献率达到了95%，也就是说，前三个主成分蕴含了8个主成分95%的信息，那么，选择这三个主成分作为新的数据进行模式识别，在保证信息量的同时降低了特征矩阵的维数。

105：特征提取后，使用支持向量机（Support Vector Machine，SVM）分类器^[7]对特征进行识别，进行自动睡眠分期。

用支持向量机做模式识别时，将由PCA去除冗余信息后的特征参数作为训练支持向量机的输入参数，并以睡眠分期为输出，经过训练得到基于EEG、HRV及其相干性的自动睡眠分期预测模型，然后进行自动睡眠分期。

综上所述，本方法通过计算EEG与HRV参数的相干系数，并与所提取的EEG，HRV信号的频域、非线性特征融合作为特征矩阵，继而利用PCA主成分分析去除冗余信息，从而准确、客观、简便的进行自动睡眠分期。该方法可有效地提高自动睡眠分期系统准确性和简便性，并获得可观的社会效益和经济效益。最佳实施方案拟采用专利转让、技术合作或产品开发。由于该技术操作简单，敏感性强，基于该技术开发的产品可应用于睡眠监测、睡眠科研等多种场景等。

参考文献

[1]Rechtschaffen A,Kales A.A Manual of standardized terminology,techniques and scoring systemfor sleep stages of human subjects[M].Washington D C:Government Printing Office,PublicHealth Service,1968:3-7

[2]封洲燕.应用小波熵分析大鼠脑电信号的动态变化特性[J].生物物理学报,2002,18(3):325-330.

[3]李小俚,崔素媛,Sleigh J W2.基于希尔伯特黄熵的麻醉深度估计[J].中国生物医学工程学报,2008,27(5):689-694.

[4]Broomhead DS,King GP.Extracting qualitative dynamics from experimental data[J].Physica D,1986,20:217-236.

[5]Scherpers HE,Von Beek JHGM,Bassingthwaighte JB.Four methods to estimate the fractaldimension from selfaffine signals[J].EngMedBiol,1992,11(6):57.

[6]Togo F,Yamamoto Y.Decreased fractal component of human heart rate variability duringnon-REM sleep[J].Am J Physiol Heart Circ Physiol.2000,280:H17–H21.

[7]Bennett K P,Campbell C.Support vector machines:hype or hallelujah?[J].ACM SIGKDDExplorations Newsletter,2000,2(2):1-13.

本领域技术人员可以理解附图只是一个优选实施例的示意图，上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种基于脑电、心率变异及两者相干性的自动睡眠分期方法，其特征在于，所述方法包括以下步骤：

(1)提取被测试者的脑电信号和心率变异信号；

(2)从脑电信号、心率变异信号中共提取9个特征向量；

(3)对9个特征向量进行主成分分析；

(4)特征降维后，使用支持向量机分类器对特征进行识别，进行自动睡眠分期；

其中，所述从脑电信号、心率变异信号中共提取9个特征向量的操作具体为：

1)通过对脑电信号进行小波变换、希尔伯特黄变换和奇异值分解获取脑电信号的小波熵、Hilbert-Huang熵，以及奇异值第一主成分，分别记为特征向量P₁,P₂,P₃；

2)计算极低频VLF、低频LF和高频HF的频谱熵，并通过小波变换计算心率变异信号的分形维数，分别记为特征向量P₄,P₅,P₆,P₇；

3)对脑电信号以及心率变异信号进行功率谱密度分析，并进行归一化处理，然后分别计算心率变异信号delta频段与心率变异信号LF、HF的相干系数，分别记为特征向量P₈,P₉。

2.根据权利要求1所述的一种基于脑电、心率变异及两者相干性的自动睡眠分期方法，其特征在于，所述对脑电信号以及心率变异信号进行功率谱密度分析，并进行归一化处理，然后分别计算脑电信号delta频段与心率变异信号LF、HF的相干系数，分别记为特征向量P₈,P₉的操作具体为：

1)截取同一时间段的EEG、HRV信号，进行预处理，包括变参考、降采样、去噪声与干扰，预处理后的EEG信号记为x，HRV信号记为y；

2)采用Welch算法计算EEG、HRV的功率谱密度以及两者的互功率谱密度，分别记为Pxx(f₁)，Pyy(f₂)，Pxy(f₁,f₂)；

3)用改进的相干函数计算EEG信号delta频段与HRV信号LF、HF的相干系数。

3.根据权利要求1所述的一种基于脑电、心率变异及两者相干性的自动睡眠分期方法，其特征在于，所述对9个特征向量进行主成分分析的操作具体为：

(1)对原始数据P_n*9进行标准化处理，矩阵中的元素减去所在列的均值，然后除以所在列的标准差，使得每个变量的均值为0，方差为1，得到矩阵P_n*9 ^*；

(2)然后求其协方差矩阵C_9*9,P_n*9 ^*中任两列之间可以计算两变量间的协方差；

(3)对协方差矩阵C_9*9进行特征根分解，得到特征根矩阵及特征向量；

(4)求原始数据P_n*9在新的矢量空间中的投影，即主成分向量组F_n*9：

(5)求累积贡献率；

(6)选定预设的累积贡献率，使前d个主成分F_n*d作为新的数据进行模式识别。