CN106371070B - 一种改进的基于小波分析的欠定盲源分离源数估计方法 - Google Patents

Abstract

本发明属于阵列信号处理技术领域，公开了一种改进的基于小波分析的欠定盲源分离源数估计方法，对观测信号应用小波变换进行分解；对分解得到的小波系数进行单层重构；重构观测信号矩阵，构造该矩阵的协方差矩阵；求解矩阵的特征值；计算源数目；本发明能够将白化后的观测信号和其部分小波分解系数重构的信号组合成为新的多维观测信号，能够在混合信号信噪比变化范围较大的情况下准确地估计得到源信号个数，原理简单，可操作性较强。

Description

一种改进的基于小波分析的欠定盲源分离源数估计方法

技术领域

本发明属于阵列信号处理技术领域。

背景技术

盲源分离(blind source separation,BSS)是指仅从若干观测到的混合信号中提取、分离出无法直接观测到的原始信号的技术^[1-2]，源数估计是盲源分离算法应用的前提，实现源数估计的方法有：基于信息论的方法、奇异值分解法以及盖尔圆法，但这些方法都要求观测信号数目大于或等于源信号数目，对于观测信号数目少于源信号数目的欠定源数估计问题，还没有十分成熟有效的方法，是盲源分离算法和信号参数估计领域的一个难点^[3-4]。

与本发明相关的现有技术

目前欠定盲源分离源数估计主要有两种方案：第一种方案是基于源信号的稀疏性实现源数估计；第二种方案是对观测信号升维实现源数估计。

尽管很多信号在时域不满足稀疏性要求，但可以通过稀疏变换工具，如FFT变换、小波变换、短时傅里叶变换等，获得信号的更稀疏表示。若不考虑噪声的影响，采样数据呈现线性聚类特征，直线的条数即为稀疏源信号的个数；如果考虑少量噪声的影响，观测信号的采样数据会分布在各条直线的附近，根据这一特性实现源数估计，相关实现算法如K均值聚类的方法、搜索重构观测信号采样点的方法、基于FFT的源数估计方法等^[3-5]，这些算法都严重依赖于源信号在时域或变换域的稀疏性，当源信号稀疏性差时，算法性能急剧下降。

对观测信号升维的方法主要有：时空法、基于小波或小波包变换的方法、基于经验模态分解(EMD)或总体经验模态分解(EEMD)的方法[6-7]。时空法通过对观测信号延时获得更多的观测信号，该方法中对于延迟参数的设置主要依靠工程经验和先验条件，且假设源信号为平稳、非重频信号，在应用中有很大的局限性。基于小波或小波包变换的方法，将每一路观测信号分解为一系列低频和高频小波系数，然后对这些小波系数进行重构得到小波系数重构信号，选择全部或部分重构信号，与原始观测信号组合成为新的观测信号，从而实现升维的目的，该方法在应用小波或小波包变换时，小波基函数和小波分解层数的选择缺乏充分理论分析和指导，后续实现原始观测信号扩维时如何选择小波系数重构信号，没有充分理论分析，在欠定源数估计中效果不佳。基于EMD或EEMD的方法，将每一路原始观测信号分解为一系列本征模态函数和一个趋势项，选择其中全部或部分分解得到的信号作为原始观测信号的补充，EMD本身存在模态混叠的缺陷，EEMD抑制了模态混叠，但需要多次重复EMD分解求均值，计算量大，如何选择补充信号缺乏理论指导，需要源信号的先验信息，导致该方法局限性较大，在欠定盲源分离源数估计方面实用性不强。

对观测信号升维实现源数估计的方法主要有：时空法、基于小波或小波包变换的方法、基于经验模态分解或总体经验模态分解的方法。时空法延迟参数的设置主要依靠工程经验和先验条件，且要求源信号为平稳、非重频信号；基于小波或小波包变换的方法，对于如何选择重构信号、小波基函数和小波分解层数缺乏充分理论分析；基于EMD的方法本身存在模态混叠的缺陷，基于EEMD抑制了模态混叠，但需要多次重复EMD分解求均值，计算量大，如何选择补充信号缺乏理论指导，需要源信号的先验信息，导致该方法局限性较大。由于以上原因，三种方法在欠定盲源分离源数估计方面效果不佳。

其中的稀疏信号：在绝大多数采样点上为零或者较小，在少数采样点幅值较大的信号，非零值采样点越少，信号越稀疏^[1]。

参考文献

[1]何昭水,谢胜利,傅予力.信号的稀疏性分析[J].自然科学进展,2006,16(9):1167-1173.

[2]张贤达.现代信号处理(第三版)[M].北京:清华大学出版社,2015.

[3]S.G.Kim and C.D.Yoo.Underdetermined Blind Source Separation Basedon Subspace Representation[J].IEEETrans.SignalProcess.,2009,57(7):2604-2614.

[4]S.L.Xie,L.Yang,J.M.Yang,et al.Time-Frequency Approach toUnderdetermined Blind Source Separation[J].IEEE Trans.NeuralNetworksandLearningSystems,2012,23(2):306-316.

[5]李爱丽.欠定盲源分离混合矩阵估计算法的研究[D].西安:西安电子科技大学学位论文,2014.

[6]毋文峰,陈小虎,苏勋家,等.机械振动源数估计的小波方法[J].机械科学与技术,2011,30(10):1679-1682.

[7]高鹏.基于单通道盲源分离理论的故障特征提取技术[D].西安:长安大学学位论文,2015.

发明内容

本发明要解决的技术问题是：针对观测信号个数少于源信号个数的欠定盲源分离问题，基于小波分析和矩阵奇异值分解方法，提出了一种改进的基于小波分析的欠定盲源分离源数估计方法，能够在混合信号信噪比变化范围较大的情况下较准确地估计出源信号个数，原理简单，可操作性较强。

为实现上述发明目的，本发明采用的技术方案是：

一种改进的基于小波分析的欠定盲源分离源数估计方法，对观测信号应用小波变换进行分解；对分解得到的小波系数进行单层重构；重构观测信号矩阵，构造该矩阵的协方差矩阵；求解矩阵的特征值；计算源数目，其具体实现步骤如下：

1)n部辐射源发射信号表示为s(t)＝[s₁(t),s₂(t),…,s_n(t)]^T，T₀为信号采样点数，m个接收阵列接收的混合信号表示为x(t)＝[x₁(t),x₂(t),…,x_m(t)]^T。欠定盲源分离模型表示为：x(t)＝As(t)+n(t),t＝1,2,…,T₀，n(t)＝[n₁(t),n₂(t),…,n_m(t)]^T表示接收信号中叠加的m个加性噪声；

2)对原始观测数据x(t)进行去均值和白化处理；

3)对处理后的观测数据[x₁(t),x₂(t),…,x_m(t)]^T进行小波分解，第k路观测信号x_k(t)分解得到的多路小波系数为a_kj，d_kj，j＝1,2…i，k＝1,2…m，a_kj为低频小波系数，d_kj为高频小波系数；

4)对每一路原始观测信号的小波分解系数进行重构。对于第k路原始观测信号的低频小波系数a_kj重构为x_k0(t)，j＝1,2…i，最高层高频小波系数d_ki重构为x_ki(t)，k＝1,2…m；

5)将观测信号和部分小波重构信号组合成为新的观测信号，x_w(t)＝[x_k(t),x_k0(t),x_ki(t)]^T，k＝1,2…m；

6)计算多维观测信号x_w(t)的协方差矩阵R_w；

7)对协方差矩阵R_w进行奇异值分解，得到其特征值Λ_s＝{λ₁,λ₂,…λ_M}；

8)采用信号与噪声判断阈值公式计算判断阈值；

9)通过特征值与判断阈值比较，完成源信号个数估计。

由于采用如上所述的技术方案，本发明具有如下优越性：

(1)一种改进的基于小波分析的欠定盲源分离源数估计方法，能够在混合信号信噪比变化范围较大的情况下准确地估计得到源信号个数，原理简单，可操作性较强。

(2)给出了一种欠定条件下信号与噪声的判断阈值计算公式。

(3)将白化后的观测信号和其部分小波分解系数重构的信号组合成为新的多维观测信号。

附图说明

图1为改进的基于小波分析的源数估计方法流程图；

具体实施方式

以下将结合图1及实施例详细说明本发明的实施方式，本发明的实施步骤如下：

1)n部辐射源发射信号表示为s(t)＝[s₁(t),s₂(t),…,s_n(t)]^T，T₀为信号采样点数，m个接收阵列接收的混合信号表示为x(t)＝[x₁(t),x₂(t),…,x_m(t)]^T。欠定盲源分离模型表示为：x(t)＝As(t)+n(t),t＝1,2,…,T₀，n(t)＝[n₁(t),n₂(t),…,n_m(t)]^T表示接收信号中叠加的m个加性噪声；

2)对原始观测数据进行白化处理，使其每一个观测信号零均值，互不相关，且具有单位方差；

3)对白化后的观测信号x(t)的小波分解，设小波分解层数为i，x₁(t)分解得到的多路小波系数分别为a_1j，d_1j，j＝1,2…i，i为最大小波分解层数，上述小波系数的下标中，第一个“1”表示第一个原始观测信号x₁(t)，其中a_1j表示第1个观测信号的低频小波系数，d_1j表示第1个观测信号的高频小波系数；其他原始观测信号分解层数及小波分解方法按照上述方法进行操作；

4)小波系数重构，将m个观测信号分解得到的小波系数序列信号分别进行小波重构，以第一个观测信号x₁(t)分解得到的多路小波系数a_1j，d_1j，j＝1,2…i为例，将低频小波系数重构为x₁₀(t)，每一层高频小波系数重构为x_1j(t)，j＝1,2…i，其他观测信号分解得到的小波系数重构方法及表示方法与此相同；

5)将白化后的观测信号和其小波重构信号中的部分信号组合成为新的多维观测信号，以第一路观测信号x₁(t)为例说明其组合方法。选取x₁(t)，x₁₀(t)，x_1i(t)作为第一路观测信号扩维后的信号，x₁₀(t)为由x₁(t)进行i层小波系数分解后得到的低频系数重构的信号，x_1i(t)为第i层高频系数重构后的信号，其他层高频系数重构的信号并未应用。其他观测信号扩维的方法与上述相同。经过扩维后的观测信号表示为：

x_w(t)＝[x₁(t),x₁₀(t),x_1i(t),x₂(t),x₂₀(t),x_2i(t),…,x_m(t),x_m0(t),x_mi(t),]^T

6)计算多维观测信号x_w(t)的协方差矩阵R_w，T₀为观测信号采样点数。

7)协方差矩阵R_w的奇异值分解。计算矩阵R_w的奇异值分解得到其特征值Λ_s＝{λ₁,λ₂,…λ_M}，其中M＝3×m为升维后的观测信号的维数。

8)计算信号与噪声判断阈值。欠定盲源混合模型中，给出一种信号与噪声判断阈 值公式

9)源数估计，令λ_s为一个M维的数据集合，统计集合中大于1的特征值的个数，即大于1的物理意义是该特征值为有用信号与噪声共同作用产生的，即通过特征值与判断阈值比较，完成源信号个数估计，得到源的个数。

仿真实验

本发明可以通过以下5个仿真实验进一步说明。

仿真实验1：四个机械源信号的混合及个数估计

1.1)仿真条件设置

仿真实验的机械信号表达式及参数为：s₁(t)＝cos(2πf₁t+10)；s₂(t)＝0.5sin(2πf₂t)；s₃(t)＝cos(2πf_b1t)[1+β₁cos(2πf_r1t)]；s₄(t)＝sin(2πf_b2t)[1+β₂cos(2πf_r2t)]。式中，f₁＝100Hz，f₂＝220Hz，f_r1＝10Hz，f_b1＝50Hz，f_r2＝25Hz，f_b2＝180Hz，f_r1＝10Hz，β₁＝0.95，β₂＝1.5，采样频率为2000Hz，采样点数T₀＝1024。混合矩阵由随机函数生成，A＝1-2*rand(m,4)。混合信号叠加噪声由Matlab自带的Agwn函数进行，令小波分解层数为2，小波基函数为db4。

1.2)仿真结果

观测信号个数m取2，在叠加噪声后混合信号信噪比为[-25dB,25dB]范围内，应用本发明估计得到的源信号个数估计结果如表1所示：

表1机械源信号的混合及源数估计实验结果

仿真实验2：四个振动信号的混合及源个数估计

2.1)仿真条件设置

四个典型振动信号表达式及参数为：s₁(t)＝sin(2π·20t)；s₂(t)＝sin(2π·100t²)；s₃(t)＝sin(2π·100t)·sin(2π·5t)；s₄(t)＝2e^-5tsin(2π·50t)。其中，s₁为简单正弦信号，s₂为线性调频扫描信号，s₃为拍波信号，s₄为指数衰减正弦信号。采样频率为1024Hz，采样点数为5000。混合矩阵A＝[0.97900.20770.61320.8432；0.06290.15690.02650.0904]，小波基函数为db4，小波分解层数为2。

2.2)仿真结果

观测信号个数m取2，在叠加噪声后混合信号信噪比为[-25dB,25dB]范围内，应用本发明估计得到的源信号个数结果如表2所示：

表2振动信号的混合及源数估计实验结果

仿真实验3：四个雷达信号的混合及源数估计

3.1)仿真条件设置

四个雷达信号表达式如下：s₁为AM信号，信号为0.5+0.5cos(2πf₁t)sin(2πf₀t)，f₀＝10MHz，f₁＝0.2MHz；s₂为LFM信号，调频带宽为20MHz，载波频率为6MHz，脉冲重复频率为0.1MHz；s₃为NLFM信号，脉宽为7μs，脉内信号为cos[2π(f₀t+a₁t²/2+a₂t³/2)]，载频f₀为10MHz，a₁＝5·10¹¹Hz/s，a₂＝10¹⁶Hz/s²；s₄为FSK信号，子脉冲宽度为(15/13)μs，载波频率变化范围为[25MHz,50MHz]，编码序列(十六进制)为1F35，混合矩阵A＝[0.97900.20770.61320.8432；0.06290.15690.02650.0904]，混合信号叠加噪声由Matlab里面的Agwn函数进行，采样频率为200MHz，采样点数为2990，令小波分解层数为2，小波基函数为db4。

3.2)仿真结果

观测信号个数m取2，在叠加噪声后混合信号信噪比为[-25dB,25dB]范围内，应用本发明得到的源数估计结果如表3所示：

表3雷达信号的混合及源数估计实验结果

仿真实验4：复杂体制雷达信号的混合及源数估计

4.1)仿真条件设置

五部复杂体制雷达，信号样式分别为：

(1)简单体制脉冲信号，载频为23MHz，脉宽为80μs，PRI为400μs，脉冲幅度为1V；

(2)重频参差信号，载频为23MHz，脉宽为80μs，相邻脉冲的PRI分别为430μs、670μs、310μs，脉冲幅度为1V；

(3)频率分集信号，载频1为23MHz，载频2为330MHz，载频3为250MHz，脉宽为80μs，PRI为400μs，脉冲幅度为1V；

(4)PRI跳变信号，载频为33MHz，脉宽为80μs，PRI集合为{430μs，442μs，422μs，398μs，410μs，470μs，453μs，465μs，438μs，458μs}，脉冲幅度为1V；

(5)PRI滑变信号，载频为23MHz，脉宽为80μs，PRI滑变规律为PRI(m)＝PRI0+(m-1)×0.01，PRI0＝430μs，m为从1到6的自然数，脉冲幅度为1V；

上述体制雷达信号脉内调制方式为单载频正弦信号，公式为0.5×sin(2π×载频×t)。采样频率为100MHz，采样点数为100001。观测天线阵列数为3，混合矩阵为A＝[0.71630.9256 0.8251 0.2760 0.8507；0.1556 0.7811 0.8784 0.6557 0.5178；0.2976 0.65350.10990.8975 0.4768]。小波基函数为db4，小波分解层数为3。

4.2)仿真结果

在叠加噪声后混合信号信噪比为[-25dB,25dB]范围内，应用本发明得到的源信数估计结果如表4所示：

表4复杂体制雷达信号的混合及源数估计实验结果

仿真实验5：两个正弦调制信号和一个符号信号的混合及源数估计

5.1)仿真条件设置

两个正弦调制信号和一个符号信号设置为：s₁(t)＝sin[2π·300t/f₁-6cos(2π·60t/f₂)]；s₂(t)＝sin(2π·9t/f₁)·sin(2π·300t/f₁)；s₃(t)＝sign[cos(2π·155t/f₁)]。式中，f₁＝10kHz，f₂＝5kHz，采样点数为3000，观测信号个数m取2，混合矩阵A＝[-0.76180.9839-0.1129；0.64780.17860.9929]，小波基函数为db4，小波分解层数为3。

5.2)仿真结果

在叠加噪声后混合信号信噪比为[-25dB,25dB]范围内，应用本发明得到的源信数估计结果如表5所示：

表5正弦调制信号和符号信号的混合及源数估计实验结果

仿真结果分析

(1)源信号样式的影响

仿真实验1与仿真实验2源信号样式为机械振动信号，仿真实验3为多种调制方式雷达信号的混合实验，具体信号波形为脉冲与连续波，且4个信号中有2个信号出现同频情况，仿真实验4为复杂体制雷达信号的混合实验，仿真实验5为正弦调制信号与符号信号的混合实验。通过仿真发现，应用本发明对于不同的信号样式的混合，绝大部分情况下都正确地估计出了源信号个数，表明本发明对于多种样式的辐射源信号的混合及源数估计具有较强的适用性。

(2)混合矩阵的影响

仿真实验1所用混合矩阵由rand函数生成，将4个源信号混合成2个观测信号；仿真实验2和仿真实验3为一固定矩阵，将4个源信号混合成2个观测信号；仿真实验4所用混合矩阵为一固定矩阵，将5个源信号混合成3个观测信号，仿真实验5所用混合矩阵为一固定矩阵，将3个源信号混合成2个观测信号。分析可知，上述实验中混合矩阵的选取具有一定的普遍性，实验结果证明本发明对于欠定盲源分离的混合矩阵具有较普遍的适用性。

(3)混合信号信噪比的影响

5个仿真实验，混合信号的信噪比变化范围都是[-25dB,25dB]，实际中传感器能够检测到的信号信噪比通常不会有如此大的变化范围。应用本发明，绝大部分情况下都正确地估出了源个数，可以认为本发明对于混合信号的信噪比具有较好的鲁棒性，能够满足工程精度的要求。

(4)小波基函数及分解层数的影响

在小波分析理论中，小波基函数及分解层数的选择是小波分析应用的关键步骤。symN系列小波基函数和dbN系列小波基函数具有优良的特性，在信号的小波分析中效果比较突出，因此通常被选为小波基函数。分解层数对于信号的降噪，是一个重要的问题，一般分解2、3层，就会取得较好的应用效果。在上述仿真实验中，小波基函数选择了db4，分解层数选择2层或3层，实际上针对多个算例，本发明在仿真中所设置分解层数最高到了5层，应用本发明估计源数目，结果普遍比较准确，基本能够满足工程应用精度要求。

Claims (1)

1.一种改进的基于小波分析的欠定盲源分离源数估计方法，其特征是：对观测信号应用小波变换进行分解；对分解得到的小波系数进行单层重构；重构观测信号矩阵，构造该矩阵的协方差矩阵；求解矩阵的特征值；计算源数目，其具体实现步骤如下：

1)n部辐射源发射信号表示为s(t)＝[s₁(t),s₂(t),…,s_n(t)]^T，T₀为信号采样点数，m个接收阵列接收的混合信号表示为x(t)＝[x₁(t),x₂(t),…,x_m(t)]^T；

欠定盲源分离模型表示为：x(t)＝As(t)+n(t),t＝1,2,…,T₀，n(t)＝[n₁(t),n₂(t),…,n_m(t)]^T表示接收信号中叠加的m个加性噪声；

2)对原始观测数据x(t)进行去均值和白化处理；

3)对处理后的观测数据[x₁(t),x₂(t),…,x_m(t)]^T进行小波分解，第k路观测信号x_k(t)分解得到的多路小波系数为a_kj，d_kj，j＝1,2…i，k＝1,2…m，a_kj为低频小波系数，d_kj为高频小波系数；

4)对每一路原始观测信号的小波分解系数进行重构；对于第k路原始观测信号的低频小波系数a_kj重构为x_k0(t)，j＝1,2…i，最高层高频小波系数d_ki重构为x_ki(t)，k＝1,2…m；

5)将观测信号和部分小波重构信号组合成为新的观测信号，x_w(t)＝[x_k(t),x_k0(t),x_ki(t)]^T，k＝1,2…m；

6)计算多维观测信号x_w(t)的协方差矩阵R_w；

7)对协方差矩阵R_w进行奇异值分解，得到其特征值Λ_s＝{λ₁,λ₂,…λ_M}；

8)采用信号与噪声判断阈值公式M＝3×m，计算判断阈值；

9)通过特征值与判断阈值比较，完成源信号个数估计。