CN103973383B - 基于Cholesky分解与特征值的协作频谱检测方法 - Google Patents

Abstract

本发明公开了一种基于Cholesky分解与特征值的协作频谱检测方法，主要解决现有的基于特征值协作频谱检测方法的检测门限难以确定的问题。其实现步骤是：(1)各个检测用户根据所要检测的频段，采集该频段的数据，并上传到处理中心；(2)处理中心根据上传的数据，构造归一化协方差矩阵，并对其进行Cholesky分解；(3)利用该分解结果计算检测统计量，分析其概率分布；(4)根据检测统计量的概率分布计算目标虚警概率下的判决门限；(5)处理中心对检测统计量与检测门限进行比较，判决主用户信号是否存在。本发明具有检测门限精准，检测性能高的优点，可用于无线通信。

Description

基于Cholesky分解与特征值的协作频谱检测方法

技术领域

本发明属于无线通信技术领域，涉及一种频谱检测方法，可用于感知网络和认知无线电系统中的频谱检测。

背景技术

随着感知网络的发展及其在日常生活中的普及，面向感知网络的关键技术的研究得到了广泛关注。其中，频谱检测技术是感知网络中的重要关键技术，通过检测频段的空闲与否来决定是否利用当前频段。同时，随着无线通信和移动通信的发展，人们对通信技术的要求越来越高，所需的服务越来越多，这使得有限的频谱资源变得日益稀缺。在目前的静态频谱分配框架下，很多频谱资源被分配给一些特定的服务，使得频谱的利用率低下。为了改善频谱利用率低下的现状，J.Mitola等人提出了认知无线电的概念，其主要思想是在已授权的频段内寻找空闲频段，在不影响授权用户正常通信的前提下，允许检测用户能够检测并接入到当前空闲的频段，从而大幅提高频谱利用率。为了能接入到空闲的频段，检测用户必须准确地检测其周围的频谱占用情况，因此频谱检测技术在认知无线电中具有关键的作用。该项技术包括频谱检测方法和协作频谱检测方法。

现有的频谱检测方法，主要有三种：

1)能量检测。检测用户通过计算接收到信号的能量，根据信号的能量大小来确定主用户是否存在。该方法实现简单，容易确定检测门限。然而，在低信噪比情况下，由于受到深度衰落和多径衰落等因素的干扰，该方法不能有效地正常工作。而且，该方法受噪声不确定度的干扰，在实际应用中受限。

2)基于循环平稳的检测。检测用户利用主用户信号在循环频率处所表现的峰值特性与噪声在循环频率处没有峰值特性来检测主用户是否存在。该方法抗噪性能好但需要主用户的先验信息且复杂度高。在认知无线电中其系统效率低，在实际应用中受限。

3)基于特征值分解的检测。检测用户利用主用户信号的相关性，通过对协方差矩阵进行特征值分解，来构造检测统计量。该方法能够抵抗噪声不确定度问题，且性能优于能量检测。然而，该方法只能采用无限采样点以确定近似的检测门限，其检测性能随之降低。

协作频谱检测是各检测用户通过协作来确定主用户存在与否。现有的协作频谱检测方法，包括：

1)基于能量检测的协作检测。各检测用户利用能量检测方法，通过协作确定主用户存在与否。该方法虽然实现简单，但是易受噪声不确定的干扰，而且，在低信噪比情况下，容易受深度衰落和多径衰落等因素的干扰，不能有效地正常工作。

2)基于循环平稳的协作检测。各检测用户利用基于循环平稳的检测方法，通过协作确定主用户存在与否。该方法抗造性能好，但复杂度高，且需要主用户信号的先验知识，无法做到盲检测，在实际中受限。

3)基于特征值分解的协作检测。各检测用户利用基于特征值分解的检测方法，通过协作确定主用户存在与否。该方法能够抵抗噪声不确定度的影响，检测性能好，但难以精确地确定检测门限，给实际应用造成了局限。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于Cholesky分解的分布式协作频谱检测方法，以提高检测性能，降低检测复杂度，精确确定判决门限，提高对主用户信号的检测性能。

为实现上述目的，本发明的技术方法包括如下步骤：

(1)将占用当前频段的用户信号定义为主用户，将通过检测当前频段上主用户是否存在，以试图占用该频段的用户信号定义为检测用户，将通过融合和分析各个检测用户采集的数据，以确定当前频段主用户信号是否存在的设备定义为处理中心；

(2)各个检测用户根据所要观察的频段，采集该频段的数据x_i(n)，其中n＝1,…,N；i＝1,…,M，N为各个检测用户的采样点数，M为检测用户数，各个检测用户将采集到的数据x_i(n)上传到处理中心；

(3)处理中心根据各检测用户上传的数据x_i(n)，构建检测统计量T_ξ：

(3.1)处理中心根据各个检测用户上传的数据x_i(n)，构建数据矩阵X和协方差矩阵R_x，其中数据矩阵X为：

$X=\left[\begin{array}{cccc}{x}_1\left(1\right)& x_1\left(2\right)& ...& x_1\left(N\right)\\ x_2\left(1\right)& x_2\left(2\right)& ...& x_2\left(N\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ x_M\left(1\right)& x_M\left(2\right)& ...& x_M\left(N\right)\end{array}\right],$

协方差矩阵为：

$R_x=\frac{1}{N}X{X}^H,$

其中(·)^H为Heimitian转置；

(3.2)处理中心根据协方差矩阵R_x，计算归一化协方差矩阵R′_x：

$R_x^{\text{'}}=\frac{N}{{\mathrm{\σ}}_w^2}{R}_x,$

其中，为处理中心设置的噪声方差，N为各个检测用户的采样点数；

(3.3)处理中心对归一化协方差矩阵R′_x进行Cholesky分解，得到分解后的上三角矩阵，即：

R′_x=L^TL，

其中，L为上三角矩阵，其表示为：

$L=\left[\begin{array}{cccc}{l}_{11}& l_{12}& ...& l_{1M}\\ 0& l_{22}& ...& l_{2M}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 0& 0& ...& l_{\mathrm{MM}}\end{array}\right],$

其中，l_ij为上三角矩阵L的第i行第j列元素，i＝1,…,M,j＝1,…,M；

(3.4)处理中心根据分解后得到的上三角矩阵L，构建检测统计量T_ξ：

T_ξ＝λ_ξ，

其中λ_ξ为上三角矩阵L的最大特征值；

(4)处理中心根据检测统计量T_ξ，计算检测门限γ_ξ：

${\mathrm{\γ}}_{\mathrm{\ξ}}=u_{{\mathrm{\λ}}_{\mathrm{\ξ}}}+{\mathrm{\σ}}_{{\mathrm{\λ}}_{\mathrm{\ξ}}}{\mathrm{\Φ}}^{-1}\left(P_{\mathrm{fa}}\right),$

其中，P_fa为虚警概率，取值范围为(0,1)，Φ^-1(·)为标准正态分布的累积量分布函数Φ(·)的逆函数，表示如下：

$\mathrm{\Φ}\left(x\right)=\underset{-\∞}{\overset{x}{\∫}}\frac{1}{2\mathrm{\π}}{e}^{-\frac{u^2}{2}}\mathrm{du},$

其中，自变量x的取值为(-∞,+∞)；

和分别为最大特征值λ_ξ的均值和标准差，分别表示如下：

$\left\{\begin{array}{c}{u}_{{\mathrm{\λ}}_{\mathrm{\ξ}}}=\underset{-\∞}{\overset{+\∞}{\∫}}{\mathrm{yF}}_{{\mathrm{\λ}}_{\mathrm{\ξ}}}^{\′}\left(y\right)\mathrm{dy}\\ {\mathrm{\σ}}_{{\mathrm{\λ}}_{\mathrm{\ξ}}}=\sqrt{\underset{-\∞}{\overset{+\∞}{\∫}}{(y-u_{{\mathrm{\λ}}_{\mathrm{\ξ}}})}^2{F}_{{\mathrm{\λ}}_{\mathrm{\ξ}}}^{\′}\left(y\right)\mathrm{dy}}\end{array}\right.,$

其中，表示最大特征值λ_ξ的概率分布的累积量分布函数，表示为：

$F_{{\mathrm{\λ}}_{\mathrm{\ξ}}}\left(y\right)=\underset{i=1}{\overset{M}{\mathrm{\Π}}}{\∫}_0^{y^2}\frac{x^{\frac{N-i-1}{2}}{e}^{-\frac{x}{2}}}{2^{\frac{N-i+1}{2}}\mathrm{\Γ}\left(\frac{N-i+1}{2}\right)}\mathrm{dx},$

其中，N为每个检测用户的采样点数，M为检测用户数，Γ(·)为伽玛函数；

(5)将步骤(3.4)中得到的检测统计量T_ξ与步骤(4)中得到的检测门限γ_ξ进行比较，当T_ξ≥γ_ξ时，判决为主用户存在，即当前频段频谱已被某用户占用，否则，判决为主用户不存在，即当前频段频谱为空闲状态，允许检测用户利用。

本发明具有以下优点：

1、本发明利用主用户信号的相关性进行检测，检测性能优于基于最大特征值的协作检测方法和基于最大和最小特征值的协作检测方法。

2、本发明是一种全盲检测方法，不需要任何有关主用户，信道和噪声的先验信息。

3、本发明基于有限随机矩阵的Cholesky分解和最大特征值提取，根据随机矩阵理论，得到检测门限的闭式表达式，能在任意的采样点数下获得与目的虚警概率对应的精确的检测门限。

4、本发明能快速确定检测门限，降低了频谱检测复杂度，可以在实际中广泛应用。

附图说明

图1是本发明的实现流程图；

图2是本发明在主用户信号不存在情况下，理论检测门限的累积量分布曲线与仿真检测门限的累积量分布曲线对比图；

图3是本发明在主用户信号不存在情况下，理论虚警概率检测门限曲线与仿真虚警概率检测门限曲线对比图；

图4是本发明与现有频谱检测方法的性噪比检测概率曲线对比图。

具体实施方式

参照图1，本发明的实现步骤如下：

步骤1，各个检测用户采集数据，并上报给处理中心。

1.1)各个检测用户根据所要检测信号的频段，用相应的滤波器滤出该频段的信号；

1.2)在满足采样定理的前提条件下，对该频段的数据进行采集，得到采集到的数据x_i(n)，其中n＝1,…,N；i＝1,…,M，其中N为各个检测用户的采样点数，M为检测用户数；

1.3)将采集的数据x_i(n)上报给处理中心。

步骤2，处理中心根据上报的数据x_i(n)，得到归一化协方差矩阵R'_x。

2.1)处理中心根据各检测用户上报的数据，构建M行N列数据矩阵X：

$X=\left[\begin{array}{cccc}{x}_1\left(1\right)& x_1\left(2\right)& ...& x_1\left(N\right)\\ x_2\left(1\right)& x_2\left(2\right)& ...& x_2\left(N\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ x_M\left(1\right)& x_M\left(2\right)& ...& x_M\left(N\right)\end{array}\right],$

其中，N为各个检测用户的采样点数，M为检测用户数；

2.2)处理中心根据构建的数据矩阵X，计算协方差矩阵R_x：

$R_x=\frac{1}{N}X{X}^H,$

其中(·)^H为Heimitian转置；

2.3)处理中心根据协方差矩阵R_x，计算归一化协方差矩阵R′_x：

$R_x^{\text{'}}=\frac{N}{{\mathrm{\σ}}_w^2}{R}_x,$

其中，为处理中心设置的噪声方差，取值为(0,+∞)。

步骤3，处理中心根据归一化协方差矩阵R'_x，按如下公式进行Cholesky分解：

R′_x=L^TL，

其中，L为上三角矩阵，其表示为：

$L=\left[\begin{array}{cccc}{l}_{11}& l_{12}& ...& l_{1M}\\ 0& l_{22}& ...& l_{2M}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 0& 0& ...& l_{\mathrm{MM}}\end{array}\right],$

其中，l_ij为上三角矩阵L的第i行j列元素，i＝1,…,M,j＝1,…,M。

步骤4，处理中心根据上三角矩阵L计算检测统计量T_ξ。

处理中心根据上三角矩阵L计算检测统计量T_ξ，其表示为：

T_ξ＝λ_ξ，

其中λ_ξ为上三角矩阵L的最大特征值。

步骤5，处理中心根据检测统计量T_ξ的表达式，分析主用户信号不存在的情况下的检测统计量T_ξ的概率分布。

5.1)上三角矩阵L的对角元素的概率分布：

在主用户信号不存在的情况下，归一化协方差矩阵R'_x为Wishart矩阵，上三角矩阵L的对角元素l_ii相互独立，且非负数，其中i＝1,…,M，且服从自由度为N-i+1的卡方分布，用公式表示为：

$\left\{\begin{array}{c}{l}_{\mathrm{ii}}\≥0\\ l_{\mathrm{ii}}^2~{\mathrm{\χ}}_{N-i+1}^2\end{array}\right.,$

其中，表示自由度为N-i+1的卡方分布；

5.2)上三角矩阵L的最大特征值的累积量分布函数：

由上三角矩阵的特性可得，上三角矩阵L的特征值即为L的对角线元素，由步骤5.1)的分析可得，上三角矩阵L的对角元素为非负数且相互独立，因此，上三角矩阵L的特征值也非负数且相互独立，

根据以上分析，可得最大特征值λ_ξ的概率分布的累积量分布函数用公式表示为：

$F_{{\mathrm{\λ}}_{\mathrm{\ξ}}}\left(y\right)=\underset{i=1}{\overset{M}{\∏}}P(l_{\mathrm{ii}}^2\≤y^2)=\underset{i=1}{\overset{M}{\mathrm{\Π}}}{\∫}_0^{y^2}\frac{x^{\frac{N-i-1}{2}}{e}^{-\frac{x}{2}}}{2^{\frac{N-i+1}{2}}\mathrm{\Γ}\left(\frac{N-i+1}{2}\right)}\mathrm{dx},$

其中，为最大特征值λ_ξ的概率分布的累积量分布函数，其中y取值为(0,+∞)，N为每个检测用户的采样点数，M为检测用户数，Γ(·)为伽玛函数。

5.3)检测统计量T_ξ的高斯近似分布：

根据中心极限定理，可将最大特征值λ_ξ的概率分布近似为高斯分布，用表示最大特征值λ_ξ的均值，表示最大特征值λ_ξ的方差，其公式分别表示为：

$\left\{\begin{array}{c}{u}_{{\mathrm{\λ}}_{\mathrm{\ξ}}}=\underset{-\∞}{\overset{+\∞}{\∫}}y{F}_{{\mathrm{\λ}}_{\mathrm{\ξ}}}^{\′}\left(y\right)\mathrm{dy}\\ {\mathrm{\σ}}_{{\mathrm{\λ}}_{\mathrm{\ξ}}}^2=\underset{-\∞}{\overset{+\∞}{\∫}}{(y-u_{{\mathrm{\λ}}_{\mathrm{\ξ}}})}^2{F}_{{\mathrm{\λ}}_{\mathrm{\ξ}}}^{\′}\left(y\right)\mathrm{dy}\end{array}\right.,$

其中，为步骤5.2)得到的最大特征值λ_ξ的概率分布的累积量分布函数。

根据上述分析，可得检测统计量T_ξ服从均值为方差为的高斯分布。

步骤6，处理中心根据检测统计量T_ξ的高斯近似分布计算判决门限γ_ξ。

根据步骤5得到的检测统计量T_ξ的概率分布，计算检测门限γ_ξ如下：

${\mathrm{\γ}}_{\mathrm{\ξ}}=u_{{\mathrm{\λ}}_{\mathrm{\ξ}}}+{\mathrm{\σ}}_{{\mathrm{\λ}}_{\mathrm{\ξ}}}{\mathrm{\Φ}}^{-1}\left(P_{\mathrm{fa}}\right),$

其中P_fa表示检测的虚警概率，取值为(0,1)，Φ^-1(·)为标准正态分布的累积量分布函数Φ(·)的逆函数。Φ(·)表示如下：

$\mathrm{\Φ}\left(x\right)=\underset{-\∞}{\overset{x}{\∫}}\frac{1}{2\mathrm{\π}}{e}^{-\frac{u^2}{2}}\mathrm{du},$

其中，自变量x的取值为(-∞,+∞)。

步骤7，处理中心根据检测统计量T_ξ和检测门限γ_ξ，判决主用户信号是否存在。

处理中心通过对步骤4中计算得到的检测统计量T_ξ与步骤6中计算得到的检测门限γ_ξ进行比较，判决主用户信号是否存在：当T_ξ≥γ_ξ时，判决为主用户存在，即当前频段频谱已被某用户占用，否则，判决为主用户不存在，即当前频段频谱为空闲状态，允许检测用户利用。

本发明的频谱检测效果可以通过以下仿真进一步说明：

A、仿真条件

主用户信号为BPSK信号，采用的噪声为均值是0，方差是1的高斯白噪声，仿真方法为10000000次的蒙特卡洛仿真。对于仿真1，检测用户数和采样点数分别设置为10和20，20和40以及40和100三种情况，虚警概率设置为0.1。对于仿真2，检测用户数和采样点数分别设置为10和20，20和40以及40和100三种情况，信噪比设置为0dB。对于仿真3，信噪比设置为从-10dB到2dB，检测用户数和采样点数分别设置为40和100，虚警概率设置为0.1。

B、仿真内容

仿真1：在主用户信号不存在情况下，对理论检测门限的累积量分布曲线与仿真检测门限的累积量分布曲线进行对比，结果如图2所示，其中，“仿真CDF”表示本发明的最大特征值的实验累积量分布曲线，“近似CDF”表示本发明的最大特征值的理论累积量分布曲线。“10和20”，“20和40”，“40和100”分别表示三种不同的检测用户数和采样点数的设置组合。

仿真2：在主用户信号不存在情况下，对理论虚警概率检测门限曲线和仿真虚警概率检测门限曲线进行对比，结果如图3所示，其中，“仿真值”表示仿真结果,“理论值”表示理论推导结果。

仿真3：对本发明和现有频谱检测方法的性噪比检测概率曲线进行对比，结果如图4所示。其中，“最大特征值算法”表示基于最大特征值的协作检测方法，“最大最小特征值算法”表示基于最大和最小特征值的协作检测方法，“提出的方法”表示本发明方法。

C、仿真结果分析

由图2可见，本发明在检测用户数和采样点数较少的情况下，得到的最大特征值的实验累积量分布曲线和理论累积量分布曲线基本吻合，因此，本发明的检测统计量的概率分布精准性高，且需要检测用户数和采样点数少，可以在实际中广泛应用。

由图3可见，本发明在检测用户数和采样点数较少的情况下，得到的理论虚警概率判决门限曲线和仿真虚警概率判决门限曲线基本吻合，因此，本发明的检测门限精准性高，且需要检测用户数和采样点数少，可以在实际中广泛应用。

由图4可见，当性噪比在-10dB到2dB间时，本发明的检测性能好于基于最大特征值的协作检测方法以及基于最大和最小特征值的协作检测方法，这说明本发明在实际中可以得到更广泛的应用。

综合上述仿真结果和分析可得，本发明所需检测用户数和采样点数少，复杂度低，检测门限精确性高，检测性能比现有的基于最大特征值的协作检测方法和基于最大和最小特征值的协作检测方法的性能好，可在实际中得到更好地应用。

Claims (2)

1.一种基于Cholesky分解的分布式协作频谱检测方法，包括如下步骤：

(1)将占用当前频段的用户信号定义为主用户，将通过检测当前频段上主用户是否存在，以试图占用该频段的用户信号定义为检测用户，将通过融合和分析各个检测用户采集的数据，以确定当前频段主用户信号是否存在的设备定义为处理中心；

(2)各个检测用户根据所要观察的频段，采集该频段的数据x_i(n)，其中n＝1,...，N；i＝1,...,M，N为各个检测用户的采样点数，M为检测用户数，各个检测用户将采集到的数据x_i(n)上传到处理中心；

(3)处理中心根据各检测用户上传的数据x_i(n)，构建检测统计量T_ξ：

(3.1)处理中心根据各个检测用户上传的数据x_i(n)，构建数据矩阵X和协方差矩阵R_x，其中数据矩阵X为：

$X=\left[\begin{array}{cccc}{x}_1\left(1\right)& x_1\left(2\right)& ...& x_1\left(N\right)\\ x_2\left(1\right)& x_2\left(2\right)& ...& x_2\left(N\right)\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ x_M\left(1\right)& x_M\left(2\right)& ...& x_M\left(N\right)\end{array}\right],$

协方差矩阵为：

$R_x=\frac{1}{N}{\mathrm{XX}}^H,$

其中(·)^H为Heimitian转置；

(3.2)处理中心根据协方差矩阵R_x，计算归一化协方差矩阵R'_x：

$R_x^{\′}=\frac{N}{{\mathrm{\σ}}_w^2}{R}_x,$

其中，为处理中心设置的噪声方差，N为各个检测用户的采样点数；

(3.3)处理中心对归一化协方差矩阵R'_x进行Cholesky分解，得到分解后的上三角矩阵，即：

R'_x＝L^TL，

其中，L为上三角矩阵，其表示为：

$L=\left[\begin{array}{cccc}{l}_{11}& l_{12}& ...& l_{1M}\\ 0& l_{22}& ...& l_{2M}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ 0& 0& ...& l_{MM}\end{array}\right],$

其中，l_ij为上三角矩阵L的第i行第j列元素，i＝1,...,M,j＝1,...,M；

(3.4)处理中心根据分解后得到的上三角矩阵L，构建检测统计量T_ξ：

T_ξ＝λ_ξ，

其中λ_ξ为上三角矩阵L的最大特征值；

(4)处理中心根据检测统计量T_ξ，计算检测门限γ_ξ：

${\mathrm{\γ}}_{\mathrm{\ξ}}=u_{{\mathrm{\λ}}_{\mathrm{\ξ}}}+{\mathrm{\σ}}_{{\mathrm{\λ}}_{\mathrm{\ξ}}}{\mathrm{\Φ}}^{-1}\left(P_{fa}\right),$

其中，P_fa为虚警概率，取值范围为(0,1)，Φ^-1(·)为标准正态分布的累积量分布函数Φ(·)的逆函数，表示如下：

$\mathrm{\Φ}\left(x\right)=\underset{-\mathrm{\∞}}{\overset{x}{\∫}}\frac{1}{2\mathrm{\π}}{e}^{-\frac{u^2}{2}}du,$

其中，自变量x的取值为(-∞，+∞)；

和分别为最大特征值λ_ξ的均值和标准差，分别表示如下：

其中，表示最大特征值λ_ξ的概率分布的累积量分布函数，表示为：

$F_{{\mathrm{\λ}}_{\mathrm{\ξ}}}\left(y\right)=\underset{i=1}{\overset{M}{\Π}}{\∫}_0^{y^2}\frac{x^{\frac{N-i-1}{2}}{e}^{-\frac{x}{2}}}{2^{\frac{N-i+1}{2}}\mathrm{\Γ}\left(\frac{N-i+1}{2}\right)}dx,$

其中，N为每个检测用户的采样点数，M为检测用户数，Γ(·)为伽玛函数；

(5)将步骤(3.4)中得到的检测统计量T_ξ与步骤(4)中得到的检测门限γ_ξ进行比较，当T_ξ≥γ_ξ时，判决为主用户存在，即当前频段频谱已被某用户占用，否则，判决为主用户不存在，即当前频段频谱为空闲状态，允许检测用户利用。

2.根据权利要求1所述的基于Cholesky分解的分布式协作频谱检测方法，其特征在于步骤(3.3)所述的处理中心对归一化协方差矩阵R'_x进行Cholesky分解，按如下公式进行：

$\left\{\begin{array}{cc}{l}_{ii}={({R_{ii}}^{\′}-\underset{k=1}{\overset{i-1}{\Σ}}{l}_{ik}^2)}^{1/2}& i=j\\ l_{ij}=\frac{{R_{ij}}^{\′}-\underset{k=1}{\overset{j-1}{\Σ}}{l}_{ik}{l}_{jk}}{l_{jj}}& j>i\end{array}\right.,$

其中R_ij'为归一化协方差矩阵R'_x的第i行j列元素，i＝1,...,M,j＝1,...,M，l_ij为上三角矩阵L的第i行j列元素，i＝1,...,M,j＝1,...,M。