CN104598971A - 基于径向基函数神经网络的单位脉冲响应函数提取方法 - Google Patents

Abstract

本发明公开了一种基于径向基函数神经网络的单位脉冲响应函数提取方法，其在求解隐层权值矩阵的过程中，先将隐层权值矩阵进行稀疏表示，然后设定稀疏度，并利用正交匹配追踪OMP算法求解得到隐层权值矩阵在二维离散余弦变换下的稀疏系数矩阵，再根据稀疏系数矩阵得到隐层权值矩阵，这个过程相比利用最小二乘法求逆来得到权值矩阵的方法，不会出现奇异矩阵，并且在较强噪声下也能准确地提取出单位脉冲响应函数，鲁棒性更好，可以有效地提高网络的泛化能力，减小误差；此外，这个过程中利用了正交匹配追踪OMP算法，可以有效地提高本发明方法的运算速度，而且可以有效地提高本发明方法的计算精度。

Description

基于径向基函数神经网络的单位脉冲响应函数提取方法

技术领域

本发明涉及一种脉冲响应函数提取技术，尤其是涉及一种基于径向基函数神经网络的单位脉冲响应函数提取方法。

背景技术

随着科学技术的快速发展，人工神经网络的应用领域越来越广，特别是径向基函数神经网络因其简单性而应用广泛。在振动信号中，脉冲响应函数提取是振动信号模态分析理论的重要组成部分，许多算法(如随机子空间算法SSI、特征系统实现算法ERA等)均以脉冲响应函数为基础进行振动信号模态参数提取，虽然目前有部分算法采用直接分析振动响应信号的思路，但是要求提取脉冲响应函数的输入输出系统具备更高精度更稳定的分析结果，以适用于结构健康监测的高精度要求场合。

提取脉冲响应函数的算法大致分为频域法、时域法和小波分解法三大类。有文献证明了小波分解法与时域法具有等价性，并提出采用伪逆矩阵求解的时域法来避免最小二乘法系数矩阵的奇异性，这种伪逆矩阵求解的时域法在低噪声强度下能够取得出色的脉冲响应函数提取效果，能够有效地应用于多种实际振动检测工程。然而，一旦观测噪声(包括环境噪声、采集噪声等)和非线性振动效应较强时，那么这种伪逆矩阵求解的时域法的提取精度将会急剧下降，导致提取出来的脉冲响应函数的误差较大。

发明内容

本发明所要解决的技术问题是提供一种基于径向基函数神经网络的单位脉冲响应函数提取方法，其在噪声较强的情况下也能准确地提取出单位脉冲响应函数，鲁棒性更好。

本发明解决上述技术问题所采用的技术方案为：一种基于径向基函数神经网络的单位脉冲响应函数提取方法，其特征在于包括以下步骤：

①令U表示一个已知的随机信号集，并以矩阵的形式将U表示为U＝[u₁ u₂...u_i...u_N]，然后将U作为径向基函数神经网络的输入信号集，其中，在此N表示U中包含的输入信号的总条数，1≤i≤N，u₁、u₂、u_i和u_N对应表示U中的第1条输入信号、第2条输入信号、第i条输入信号和第N条输入信号，U中的每条输入信号为M维的列向量；

②利用K-means聚类算法对U进行聚类，获得U的K个类中心，并确定径向基函数神经网络的隐层的节点数为K，其中，K∈[1,N]；

③令Y表示径向基函数神经网络的输出信号集，且Y已知，并以矩阵的形式将Y表示为Y＝[y₁ y₂...y_i...y_N]，然后对Y中的每条输出信号添加信噪比为R的噪声信号，将含噪声信号的输出信号集记为Y'，其中，Y为U与单位脉冲响应函数卷积得到，在此N表示Y中包含的输出信号的总条数，与U中包含的输入信号的总条数一致，1≤i≤N，y₁、y₂、y_i和y_N对应表示Y中的第1条输出信号、第2条输出信号、第i条输出信号和第N条输出信号，Y中的每条输出信号为M维的列向量，R表示信噪比；

④将U的每个类中心与U中的N条输入信号进行高斯函数求解，获得由径向基函数神经网络的隐层中的所有节点的隐层输出构成的维数为K×N的隐层输出矩阵，记为Ψ，Ψ中的第k行隐层输出为径向基函数神经网络的隐层中的第k个节点的隐层输出，将Ψ中的第k行隐层输出记为中的第i个隐层输出为径向基函数神经网络的隐层中的第k个节点与U中的第i条输入信号的隐层输出，将中的第i个隐层输出记为其中，1≤k≤K，1≤i≤N，exp()表示以自然基数e为底的指数函数，λ为一常数且λ≥0，符号“|| ||”为求欧氏距离符号，c_k表示U的第k个类中心；

⑤根据径向基函数神经网络的特点，获知Ψ、隐层权值矩阵W与Y'之间存在线性关系：W×Ψ＝Y'；然后将W的转置W^T稀疏表示为W^T＝D1×S×D2；接着设定S的稀疏度为m1，利用正交匹配追踪OMP算法求解得到S；再根据求解得到的S，求解W^T＝D1×S×D2得到W^T；之后对W^T进行转置得到W，此时径向基函数神经网络为一个已训练好的径向基函数神经网络；其中，W表示径向基函数神经网络的隐层权值矩阵，W的维数为M×K，D1和D2为两个不同维度方向的二维离散余弦变换基，S表示W^T在二维离散余弦变换下的稀疏系数矩阵，1≤m1≤G_n，G_n表示采集度；

⑥将任意一个单位脉冲信号δ(t)输入到已训练好的径向基函数神经网络中；在已训练好的径向基函数神经网络中，根据U的K个类中心，获取δ(t)的隐层输出矩阵，记为Ψ'，Ψ'＝exp(-λ||δ(t)-c_k||²)，其中，exp()表示以自然基数e为底的指数函数，λ为一常数且λ≥0，符号“|| ||”为求欧氏距离符号，c_k表示U的第k个类中心；接着根据Ψ'和求解得到的W，计算单位脉冲响应函数，记为h(t)，h(t)＝W×Ψ'。

所述的步骤②中U的K个类中心的获取过程为：

②-1、随机选取U中的K条输入信号，将随机选取的每条输入信号作为U的一个类中心，其中，K∈[1,N]；

②-2、将U中当前待处理的第i条输入信号定义为当前输入信号，其中，i的初始值为1，1≤i≤N；

②-3、计算当前输入信号与U的各个类中心之间的欧式距离，然后从当前输入信号对应的K个欧式距离中找出最小欧式距离，假设最小欧式距离对应的类中心为U的第k个类中心，则确定当前输入信号属于U的第k个类中心，其中，1≤k≤K；

②-4、令i＝i+1，将U中下一个待处理的输入信号作为当前输入信号，然后返回步骤②-3继续执行，直至U中所有的输入信号处理完毕，再执行步骤②-5，其中，i＝i+1中的“＝”为赋值符号；

②-5、对属于U的任一个类中心的所有输入信号求均值，然后用该类中心对应的均值替换该类中心，实现类中心的更新；

②-6、重复执行步骤②-2至步骤②-5多次，获得U的K个最终的类中心。

所述的步骤④中取λ＝1。

所述的步骤⑤中取G_n＝0.5K。

所述的步骤⑤的具体过程为：

⑤-1、根据径向基函数神经网络的特点，获知Ψ、隐层权值矩阵W与Y'之间存在线性关系：W×Ψ＝Y'，然后对W×Ψ＝Y'两边进行转置，得到Ψ^T×W^T＝Y'^T，其中，W表示径向基函数神经网络的隐层权值矩阵，W的维数为M×K，M表示Y中的每条输出信号的维数；

⑤-2、将W的转置W^T稀疏表示为W^T＝D1×S×D2，其中，D1和D2为两个不同维度方向的二维离散余弦变换基，S表示W^T在二维离散余弦变换下的稀疏系数矩阵；

⑤-3、根据Ψ^T×W^T＝Y'^T和W^T＝D1×S×D2，得到Ψ^T×D1×S×D2＝Y'^T，进一步得到Ψ^T×D1×S＝Y'^T×D2^T，其中，D2^T为D2的转置；然后令A＝Ψ^T×D1，并令B＝Y'^T×D2^T，将Ψ^T×D1×S＝Y'^T×D2^T简化为A×S＝B；

⑤-4、令G_n表示采集度，构造一个维数为G_n×N的随机采样矩阵G，并将G作为稀疏复原过程的观测矩阵；然后利用G对A进行随机采样，得到感知矩阵，记为Ar，Ar＝G×A，并利用G对B进行随机采样，得到由多个观测向量构成的测量矩阵，记为Br，Br＝G×B；

⑤-5、设定S的稀疏度为m1，根据A×S＝B及Ar＝G×A和Br＝G×B，得到Ar×S＝Br；然后利用正交匹配追踪OMP算法求解Ar×S＝Br得到S，其中，1≤m1≤G_n；

⑤-6、根据求解得到的S，由W^T＝D1×S×D2求解W^T，之后对W^T进行转置得到W。

所述的步骤⑤-5中利用正交匹配追踪OMP算法求解Ar×S＝Br得到S的具体过程为：

⑤-5a、以ks表示迭代次数，令pp表示维数为m1的向量，令B1表示维数为G_n×K的矩阵，其中，ks的初始值为1，1≤ks≤m1；

⑤-5b、计算Ar中的每一列与第ks-1次迭代的残差值r_ks-1之间的相关系数，共得到K个相关系数；然后将K个相关系数中的最大值对应的Ar中的一列的列号作为pp中的第ks个元素的值，并将K个相关系数中的最大值对应的Ar中的一列作为B1中的第ks列；再在Ar中将K个相关系数中的最大值对应的Ar中的一列剔除；

其中，Ar中的第j列与第ks-1次迭代的残差值r_ks-1之间的相关系数为Ar中的第j列与第ks-1次迭代的残差值r_ks-1之间的内积，1≤j≤K，当ks＝1时r_ks-1为Br中的第1列观测向量；

⑤-5c、利用最小二乘法，计算r_ks-1在B1的每一列上的稀疏投影系数，得到稀疏投影系数序列，记为a_ks；

⑤-5d、根据B1和a_ks，计算第ks次迭代的残差值，记为r_ks，r_ks＝r_ks-1-B1×a_ks；

⑤-5e、判断ks是否等于m1，如果是，则得到S，S中的非零值的位置为pp中的值，相应位置所对应的非零值为a_ks中的值，S中的其它位置的值均为零；否则，令ks＝ks+1，然后返回步骤⑤-5b继续执行，其中，ks＝ks+1中的“＝”为赋值符号；

所述的步骤⑥中取λ＝1。

所述的步骤⑤中D1的构建过程为：

a1、令D1表示一个K阶方阵，令g和h为两个变量，1≤g≤K,1≤h≤K；

a2、将D1中第g行第h列的元素记为D1(g,h)， $D1(g,h)=\sqrt{2/K}\cos [\frac{\mathrm{\π}(g-0.5)}{K}\×(h-0.5)];$

所述的步骤⑤中D2的构建过程为：

b1、令D2表示一个M阶方阵，令g'和h'为两个变量，1≤g'≤M,1≤h'≤M；

b2、将D2中第g'行第h'列的元素记为D2(g',h')， $D2(g^{\′},h^{\′})=\sqrt{2/M}\cos [\frac{\mathrm{\π}(g^{\′}-0.5)}{M}\×(h^{\′}-0.5)].$

与现有技术相比，本发明的优点在于：

1)本发明方法引入两个不同的二维离散余弦变换基来构建不同的稀疏基，对径向基函数神经网络的隐层权值矩阵进行稀疏表示，利用正交匹配追踪OMP算法来获得隐层权值矩阵，这种方式获得的隐层权值矩阵的稳定性能好。

2)本发明方法在求解隐层权值矩阵的过程中，先将隐层权值矩阵进行稀疏表示，然后设定稀疏度，并利用正交匹配追踪OMP算法求解得到隐层权值矩阵在二维离散余弦变换下的稀疏系数矩阵，再根据稀疏系数矩阵得到隐层权值矩阵，这个过程相比利用最小二乘法求逆来得到权值矩阵的方法，不会出现奇异矩阵，并且在较强噪声下也能准确地提取出单位脉冲响应函数，鲁棒性更好，可以有效地提高网络的泛化能力，减小误差；此外，这个过程中利用了正交匹配追踪OMP算法，可以有效地提高本发明方法的运算速度，而且可以有效地提高本发明方法的计算精度。

附图说明

图1为本发明方法的流程框图；

图2a为稀疏度m1＝42、信噪比为R＝10、隐层的节点数即聚类中心数K＝300、输入信号的条数N＝1000的情况下，精确的单位脉冲响应函数、采用本发明方法提取出的单位脉冲响应函数与采用经典最小二乘法获得的单位脉冲响应函数、采用规则化最小二乘法获得的单位脉冲响应函数的误差对比；

图2b为稀疏度m1＝100、信噪比为R＝10、隐层的节点数即聚类中心数K＝300、输入信号的条数N＝1000的情况下，精确的单位脉冲响应函数、采用本发明方法提取出的单位脉冲响应函数与采用经典最小二乘法获得的单位脉冲响应函数、采用规则化最小二乘法获得的单位脉冲响应函数的误差对比；

图2c为稀疏度m1＝42、信噪比为R＝5、隐层的节点数即聚类中心数K＝300、输入信号的条数N＝1000的情况下，精确的单位脉冲响应函数、采用本发明方法提取出的单位脉冲响应函数与采用经典最小二乘法获得的单位脉冲响应函数、采用规则化最小二乘法获得的单位脉冲响应函数的误差对比；

图2d为稀疏度m1＝42、信噪比为R＝10、隐层的节点数即聚类中心数K＝300、输入信号的条数N＝800的情况下，精确的单位脉冲响应函数、采用本发明方法提取出的单位脉冲响应函数与采用经典最小二乘法获得的单位脉冲响应函数、采用规则化最小二乘法获得的单位脉冲响应函数的误差对比；

图2e为稀疏度m1＝42、信噪比为R＝10、隐层的节点数即聚类中心数K＝500、输入信号的条数N＝1000的情况下，精确的单位脉冲响应函数、采用本发明方法提取出的单位脉冲响应函数与采用经典最小二乘法获得的单位脉冲响应函数、采用规则化最小二乘法获得的单位脉冲响应函数的误差对比。

具体实施方式

以下结合附图实施例对本发明作进一步详细描述。

本发明提出的一种基于径向基函数神经网络的单位脉冲响应函数提取方法，其流程框图如图1所示，其包括以下步骤：

①令U表示一个已知的随机信号集，并以矩阵的形式将U表示为U＝[u₁ u₂...u_i...u_N]，然后将U作为径向基函数神经网络的输入信号集，其中，在此N表示U中包含的输入信号的总条数，500≤N≤1500，如在本实施例中可取N＝1000，1≤i≤N，u₁、u₂、u_i和u_N对应表示U中的第1条输入信号、第2条输入信号、第i条输入信号和第N条输入信号，U中的每条输入信号为M维的列向量，M≥32，在本实施例中可取M＝32， $u_i=\left[\begin{array}{c}{u}_{i,1}\\ u_{i,2}\\ \·\\ \·\\ \·\\ u_{i,t}\\ \·\\ \·\\ \·\\ u_{i,M}\end{array}\right],$ 1≤t≤M，u_i,1、u_i,2、u_i,t和u_i,M对应表示u_i中的第1个信号值、第2个信号值、第t个信号值和第M个信号值。

在此，该随机信号集中的随机信号可以是由传感器采集的大桥斜拉索的振动信号。

②利用K-means聚类算法对U进行聚类，获得U的K个类中心，并确定径向基函数神经网络的隐层的节点数为K，其中，K∈[1,N]。

在本实施例中，步骤②中U的K个类中心的获取过程为：

②-1、随机选取U中的K条输入信号，将随机选取的每条输入信号作为U的一个类中心，其中，K∈[1,N]。

②-2、将U中当前待处理的第i条输入信号定义为当前输入信号，其中，i的初始值为1，1≤i≤N。

②-3、计算当前输入信号与U的各个类中心之间的欧式距离，然后从当前输入信号对应的K个欧式距离中找出最小欧式距离，假设最小欧式距离对应的类中心为U的第k个类中心，则确定当前输入信号属于U的第k个类中心，其中，1≤k≤K。

②-4、令i＝i+1，将U中下一个待处理的输入信号作为当前输入信号，然后返回步骤②-3继续执行，直至U中所有的输入信号处理完毕，再执行步骤②-5，其中，i＝i+1中的“＝”为赋值符号。

②-5、对属于U的任一个类中心的所有输入信号求均值，然后用该类中心对应的均值替换该类中心，实现类中心的更新。

②-6、重复执行步骤②-2至步骤②-5多次，获得U的K个最终的类中心，将U的第k个最终的类中心记为c_k，U的每个类中心为M维的列向量， $c_k=\left[\begin{array}{c}{c}_{k,1}\\ c_{k,2}\\ \·\\ \·\\ \·\\ c_{k,t}\\ \·\\ \·\\ \·\\ c_{k,M}\end{array}\right],$ 1≤t≤M，c_k,1、c_k,2、c_k,t和c_k,M对应表示c_k中的第1个值、第2个值、第t个值和第M个值。在此，重复执行步骤②-2至步骤②-5的次数可设定为大于或等于20次，如具体操作时取50次。

③令Y表示径向基函数神经网络的输出信号集，且Y已知，并以矩阵的形式将Y表示为Y＝[y₁ y₂...y_i...y_N]，然后对Y中的每条输出信号添加信噪比为R的噪声信号，将含噪声信号的输出信号集记为Y'，其中，Y为U与单位脉冲响应函数卷积得到，在此N表示Y中包含的输出信号的总条数，与U中包含的输入信号的总条数一致，1≤i≤N，y₁、y₂、y_i和y_N对应表示Y中的第1条输出信号、第2条输出信号、第i条输出信号和第N条输出信号，Y中的每条输出信号为M维的列向量， $y_i=\left[\begin{array}{c}{y}_{i,1}\\ y_{i,2}\\ \·\\ \·\\ \·\\ y_{i,t}\\ \·\\ \·\\ \·\\ y_{i,M}\end{array}\right],$ 1≤t≤M，y_i,1、y_i,2、y_i,t和y_i,M对应表示y_i中的第1个信号值、第2个信号值、第t个信号值和第M个信号值，R表示信噪比，在本实施例中取R＝10，Y'＝[y₁' y₂'...y_i'...y_N']，y₁'、y₂'、y_i'和y_N'对应表示Y中的第1条输出信号、第2条输出信号、第i条输出信号和第N条输出信号各自添加信噪比为R的噪声信号后的信号，Y'中的每条添加有信噪比为R的噪声信号的输出信号为M维的列向量， ${y_i}^{\′}=\left[\begin{array}{c}{y_{i,1}}^{\′}\\ {y_{i,2}}^{\′}\\ \·\\ \·\\ \·\\ {y_{i,t}}^{\′}\\ \·\\ \·\\ \·\\ {y_{i,M}}^{\′}\end{array}\right],$ 1≤t≤M，y_i,1'、y_i,2'、y_i,t'和y_i,M'对应表示y_i'中的第1个信号值、第2个信号值、第t个信号值和第M个信号值。

在此，对Y中的每条输出信号添加信噪比为R的噪声信号的目的是为了说明输出信号在有噪声的情况下，本发明方法也能较好的提取出单位脉冲响应函数，从而可说明本发明方法的抗噪性能好。

④将U的每个类中心与U中的N条输入信号进行高斯函数求解，获得由径向基函数神经网络的隐层中的所有节点的隐层输出构成的维数为K×N的隐层输出矩阵，记为Ψ，Ψ中的第i列隐层输出为U中的第i条输入信号与径向基函数神经网络的隐层中的所有节点的隐层输出，Ψ中的第k行隐层输出为径向基函数神经网络的隐层中的第k个节点与U中的所有输入信号的隐层输出，Ψ中的第k行隐层输出亦可直接理解为径向基函数神经网络的隐层中的第k个节点的隐层输出，将Ψ中的第k行隐层输出即将径向基函数神经网络的隐层中的第k个节点的隐层输出记为中的第i个隐层输出为径向基函数神经网络的隐层中的第k个节点与U中的第i条输入信号的隐层输出，将中的第i个隐层输出记为其中，1≤k≤K，1≤i≤N，exp()表示以自然基数e为底的指数函数，λ为一常数且λ≥0，在本实施例中取λ＝1，符号“|| ||”为求欧氏距离符号，c_k表示U的第k个类中心。在此，高斯函数作为隐层函数，可具体选择高斯径向基函数。

⑤根据径向基函数神经网络的特点，获知Ψ、隐层权值矩阵W与Y'之间存在线性关系：W×Ψ＝Y'；然后将W的转置W^T稀疏表示为W^T＝D1×S×D2；接着设定S的稀疏度为m1，利用正交匹配追踪OMP算法求解得到S；再根据求解得到的S，求解W^T＝D1×S×D2得到W^T；之后对W^T进行转置得到W，此时径向基函数神经网络为一个已训练好的径向基函数神经网络；其中，W表示径向基函数神经网络的隐层权值矩阵，W的维数为M×K，D1和D2为两个不同维度方向的二维离散余弦变换基，S表示W^T在二维离散余弦变换下的稀疏系数矩阵，1≤m1≤G_n，G_n表示采集度，G_n＝0.5K，稀疏度m1表示稀疏系数的非零值个数，以便求解稀疏系数矩阵。

在此具体实施例中，步骤⑤的具体过程为：

⑤-1、根据径向基函数神经网络的特点，获知Ψ、隐层权值矩阵W与Y'之间存在线性关系：W×Ψ＝Y'，然后对W×Ψ＝Y'两边进行转置，得到Ψ^T×W^T＝Y'^T，其中，W表示径向基函数神经网络的隐层权值矩阵，W的维数为M×K，M表示Y中的每条输出信号的维数，在本实施例中可取M＝32，W^T的维数为K×M，Ψ^T的维数为N×K，Y'^T的维数为N×M。

⑤-2、由于隐层权值矩阵W的转置W^T不一定稀疏，因此本发明先选用两个不同维度方向的二维离散余弦变换基(DCT)对W的转置W^T进行稀疏表示，即将W的转置W^T稀疏表示为W^T＝D1×S×D2，其中，D1和D2为两个不同维度方向的二维离散余弦变换基，S表示W^T在二维离散余弦变换下的稀疏系数矩阵；由于一般矩阵在二维离散余弦变换下的系数值集中在低频，因此可以假定S是稀疏的，由此把问题从求解一般非稀疏的W^T转变成求解稀疏的S。

在具体实施时，如果D1为X方向的二维离散余弦变换基，则D2为Y方向的二维离散余弦变换基；如果D2为X方向的二维离散余弦变换基，则D1为Y方向的二维离散余弦变换基。

在本实施例中，D1的构建过程为：

a1、令D1表示一个K阶方阵，令g和h为两个变量，1≤g≤K,1≤h≤K。

a2、将D1中第g行第h列的元素记为D1(g,h)， $D1(g,h)=\sqrt{2/K}\cos [\frac{\mathrm{\π}(g-0.5)}{K}\×(h-0.5)].$

在本实施例中，D2的构建过程为：

b1、令D2表示一个M阶方阵，令g'和h'为两个变量，1≤g'≤M,1≤h'≤M。

b2、将D2中第g'行第h'列的元素记为D2(g',h')， $D2(g^{\′},h^{\′})=\sqrt{2/M}\cos [\frac{\mathrm{\π}(g^{\′}-0.5)}{M}\×(h^{\′}-0.5)].$

⑤-3、根据Ψ^T×W^T＝Y'^T和W^T＝D1×S×D2，得到Ψ^T×D1×S×D2＝Y'^T，进一步得到Ψ^T×D1×S＝Y'^T×D2^T，其中，D2^T为D2的转置；然后令A＝Ψ^T×D1，并令B＝Y'^T×D2^T，将Ψ^T×D1×S＝Y'^T×D2^T简化为A×S＝B。

⑤-4、令G_n表示采集度，构造一个维数为G_n×N的随机采样矩阵G，并将G作为稀疏复原过程的观测矩阵；然后利用G对A进行随机采样，得到感知矩阵，记为Ar，Ar＝G×A，并利用G对B进行随机采样，得到由多个观测向量构成的测量矩阵，记为Br，Br＝G×B，其中，G_n＝0.5K。

⑤-5、设定S的稀疏度为m1，根据A×S＝B及Ar＝G×A和Br＝G×B，得到Ar×S＝Br；然后利用正交匹配追踪OMP算法求解Ar×S＝Br得到S，其中，1≤m1≤G_n。

在此，步骤⑤-5中利用正交匹配追踪OMP算法求解Ar×S＝Br得到S的具体过程为：

⑤-5a、以ks表示迭代次数，令pp表示维数为m1的向量，令B1表示维数为G_n×K的矩阵，其中，ks的初始值为1，1≤ks≤m1。

⑤-5b、计算Ar中的每一列与第ks-1次迭代的残差值r_ks-1之间的相关系数，共得到K个相关系数；然后将K个相关系数中的最大值对应的Ar中的一列的列号作为pp中的第ks个元素的值，并将K个相关系数中的最大值对应的Ar中的一列作为B1中的第ks列；再在Ar中将K个相关系数中的最大值对应的Ar中的一列剔除。

其中，Ar中的第j列与第ks-1次迭代的残差值r_ks-1之间的相关系数为Ar中的第j列与第ks-1次迭代的残差值r_ks-1之间的内积，1≤j≤K，当ks＝1时r_ks-1为Br中的第1列观测向量。

⑤-5c、利用最小二乘法，计算r_ks-1在B1的每一列上的稀疏投影系数，得到稀疏投影系数序列，记为a_ks。

⑤-5d、根据B1和a_ks，计算第ks次迭代的残差值，记为r_ks，r_ks＝r_ks-1-B1×a_ks。

⑤-5e、判断ks是否等于m1，如果是，则得到S，S中的非零值的位置为pp中的值，相应位置所对应的非零值为a_ks中的值，S中的其它位置的值均为零；否则，令ks＝ks+1，然后返回步骤⑤-5b继续执行，其中，ks＝ks+1中的“＝”为赋值符号。

⑤-6、根据求解得到的S，由W^T＝D1×S×D2求解W^T，之后对W^T进行转置得到W。

⑥将任意一个单位脉冲信号δ(t)输入到已训练好的径向基函数神经网络中；在已训练好的径向基函数神经网络中，根据U的K个类中心，获取δ(t)的隐层输出矩阵，记为Ψ'，Ψ'＝exp(-λ||δ(t)-c_k||²)，其中，exp()表示以自然基数e为底的指数函数，λ为一常数且λ≥0，在本实施例中取λ＝1，符号“|| ||”为求欧氏距离符号，c_k表示U的第k个类中心；接着根据Ψ'和求解得到的W，计算单位脉冲响应函数，记为h(t)，h(t)＝W×Ψ'。

为进一步说明本发明方法的可行性，对本发明方法进行仿真试验。

取随机的输入信号集U＝[u₁ u₂...u_i...u_N]，其中第i条输入信号为u_i＝0.05×randn(M,1)，randn(M,1)为一个维数为M×1的随机信号序列，在本实例中取M＝32，取N＝1000，表示随机的输入信号集中有1000条随机输入信号。在仿真时，输出信号集Y是输入信号集U和仿真的理想单位脉冲响应函数卷积获得。

已知一个理想的单位脉冲响应函数为h(t)，在本实例中构造的理想的单位脉冲响应函数是h(t)＝0.9×exp(-5t)×cos(2π×10t+0.1)+0.5×exp(-3t)×cos(2π×30t+0.05)。

图2a给出了稀疏度m1＝42、信噪比R＝10、隐层的节点数K＝300、输入信号的条数N＝1000的情况下，精确的单位脉冲响应函数、采用本发明方法提取出的单位脉冲响应函数与采用经典最小二乘法获得的单位脉冲响应函数、采用规则化最小二乘法获得的单位脉冲响应函数的误差对比；图2b给出了稀疏度m1＝100、信噪比为R＝10、隐层的节点数即聚类中心数K＝300、输入信号的条数N＝1000的情况下，精确的单位脉冲响应函数、采用本发明方法提取出的单位脉冲响应函数与采用经典最小二乘法获得的单位脉冲响应函数、采用规则化最小二乘法获得的单位脉冲响应函数的误差对比；图2c给出了稀疏度m1＝42、信噪比为R＝5、隐层的节点数即聚类中心数K＝300、输入信号的条数N＝1000的情况下，精确的单位脉冲响应函数、采用本发明方法提取出的单位脉冲响应函数与采用经典最小二乘法获得的单位脉冲响应函数、采用规则化最小二乘法获得的单位脉冲响应函数的误差对比；图2d给出了稀疏度m1＝42、信噪比为R＝10、隐层的节点数即聚类中心数K＝300、输入信号的条数N＝800的情况下，精确的单位脉冲响应函数、采用本发明方法提取出的单位脉冲响应函数与采用经典最小二乘法获得的单位脉冲响应函数、采用规则化最小二乘法获得的单位脉冲响应函数的误差对比；图2e给出了稀疏度m1＝42、信噪比为R＝10、隐层的节点数即聚类中心数K＝500、输入信号的条数N＝1000的情况下，精确的单位脉冲响应函数、采用本发明方法提取出的单位脉冲响应函数与采用经典最小二乘法获得的单位脉冲响应函数、采用规则化最小二乘法获得的单位脉冲响应函数的误差对比。

从图2a中可以看出，采用本发明方法提取出的单位脉冲响应函数的误差最小，仅为0.0322；图2b、图2c、图2d、图2e分别为稀疏度m1、信噪比R、隐层的节点数K及输入信号的条数N中几个变量中的一个变化时，采用本发明方法提取出的单位脉冲响应函数与采用经典最小二乘法和规则化最小二乘法获得的单位脉冲响应函数的误差对比，从图2a至图2e中可以看出采用经典最小二乘法获得的单位脉冲响应函数都会出现很大的误差，而规则化最小二乘法也会使得误差增大。在图2a中，当稀疏度m1＝42、信噪比R＝10、隐层的节点数K＝300、输入信号的条数N＝1000时，从整体上对比三种方法的准确性，从图2a中可获知采用本发明方法提取出的单位脉冲响应函数与仿真的精确单位脉冲响应函数间的误差仅为0.0322，而经典最小二乘法的误差为0.0748，正则化最小二乘法的误差为0.0507，都比本发明方法的误差要大，说明采用本发明方法提取出的单位脉冲响应函数更加准确，产生的误差小，具有优越性；在图2b中，所示的是仅当稀疏度m1变大即m1＝100其他条件不变时，采用本发明方法提取出的单位脉冲响应函数的误差为0.0493，相对于稀疏度m1＝42时的误差0.0322变大了，与采用规则化最小二乘法提取出的单位脉冲响应函数的误差为0.0486相差不大，比采用经典最小二乘法提取出的单位脉冲响应函数的误差0.1050小很多，说明本发明方法在稀疏度变大时，也是可行的；图2c所示的是只改变信噪比R的时候，发现在小信噪比R＝5时即噪声较大的时候，采用本发明方法提取出的单位脉冲响应函数的误差为0.0483，相比在R＝10的情况下的误差0.0322增大，却也比采用经典最小二乘法提取出的单位脉冲响应函数的误差为0.1774和采用规则化最小二乘法提取出的单位脉冲响应函数的误差为0.0820小，说明在噪声较大时，采用本发明方法提取出的单位脉冲响应函数更准确，抗噪性能好；图2d所示的只改变输入信号的条数N大小时，当N＝1000和N＝800时的误差变化并无很大区别，说明不需要1000条的大输入样本时也可获得较好的单位脉冲响应函数，且采用本发明方法提取出的单位脉冲响应函数的误差为0.0254，比采用经典最小二乘法提取出的单位脉冲响应函数的误差0.1522和采用规则化最小二乘法提取出的单位脉冲响应函数的误差0.1069小的多，说明采用本发明方法提取出的单位脉冲响应函数依旧更准确；图2e所示的只改变隐层的节点数K时，误差随之改变，当隐层的节点数K越大，误差变大，获得的单位脉冲响应函数越趋于不准确，说明选择合理的隐层的节点数很重要。综上所述，通过不同对比，采用本发明方法提取单位脉冲响应函数是可行的，并且相对于采用经典最小二乘法和规则化最小二乘法提取出的单位脉冲响应函数更准确，在噪声比较大时，采用本发明方法得到的单位脉冲响应函数更精确，误差更小，效果更好，泛化性能好，具有实用性，而采用经典最小二乘法提取出的单位脉冲响应函数的误差最大，极其不稳定，这是因为经典最小二乘法会产生奇异矩阵的问题。

Claims (8)

1.一种基于径向基函数神经网络的单位脉冲响应函数提取方法，其特征在于包括以下步骤：

①令U表示一个已知的随机信号集，并以矩阵的形式将U表示为U＝[u₁u₂...u_i...u_N]，然后将U作为径向基函数神经网络的输入信号集，其中，在此N表示U中包含的输入信号的总条数，1≤i≤N，u₁、u₂、u_i和u_N对应表示U中的第1条输入信号、第2条输入信号、第i条输入信号和第N条输入信号，U中的每条输入信号为M维的列向量；

②利用K-means聚类算法对U进行聚类，获得U的K个类中心，并确定径向基函数神经网络的隐层的节点数为K，其中，K∈[1,N]；

③令Y表示径向基函数神经网络的输出信号集，且Y已知，并以矩阵的形式将Y表示为Y＝[y₁y₂...y_i...y_N]，然后对Y中的每条输出信号添加信噪比为R的噪声信号，将含噪声信号的输出信号集记为Y'，其中，Y为U与单位脉冲响应函数卷积得到，在此N表示Y中包含的输出信号的总条数，与U中包含的输入信号的总条数一致，1≤i≤N，y₁、y₂、y_i和y_N对应表示Y中的第1条输出信号、第2条输出信号、第i条输出信号和第N条输出信号，Y中的每条输出信号为M维的列向量，R表示信噪比；

④将U的每个类中心与U中的N条输入信号进行高斯函数求解，获得由径向基函数神经网络的隐层中的所有节点的隐层输出构成的维数为K×N的隐层输出矩阵，记为Ψ，Ψ中的第k行隐层输出为径向基函数神经网络的隐层中的第k个节点的隐层输出，将Ψ中的第k行隐层输出记为 中的第i个隐层输出为径向基函数神经网络的隐层中的第k个节点与U中的第i条输入信号的隐层输出，将中的第i个隐层输出记为  其中，1≤k≤K，1≤i≤N，exp( )表示以自然基数e为底的指数函数，λ为一常数且λ≥0，符号“|| ||”为求欧氏距离符号，c_k表示U的第k个类中心；

⑤根据径向基函数神经网络的特点，获知Ψ、隐层权值矩阵W与Y'之间存在线性关系：W×Ψ＝Y'；然后将W的转置W^T稀疏表示为W^T＝D1×S×D2；接着设定S的稀疏度为m1，利用正交匹配追踪OMP算法求解得到S；再根据求解得到的S，求解 W^T＝D1×S×D2得到W^T；之后对W^T进行转置得到W，此时径向基函数神经网络为一个已训练好的径向基函数神经网络；其中，W表示径向基函数神经网络的隐层权值矩阵，W的维数为M×K，D1和D2为两个不同维度方向的二维离散余弦变换基，S表示W^T在二维离散余弦变换下的稀疏系数矩阵，1≤m1≤G_n，G_n表示采集度；

⑥将任意一个单位脉冲信号δ(t)输入到已训练好的径向基函数神经网络中；在已训练好的径向基函数神经网络中，根据U的K个类中心，获取δ(t)的隐层输出矩阵，记为Ψ'，Ψ'＝exp(-λ||δ(t)-c_k||²)，其中，exp( )表示以自然基数e为底的指数函数，λ为一常数且λ≥0，符号“|| ||”为求欧氏距离符号，c_k表示U的第k个类中心；接着根据Ψ'和求解得到的W，计算单位脉冲响应函数，记为h(t)，h(t)＝W×Ψ'。

2.根据权利要求1所述的基于径向基函数神经网络的单位脉冲响应函数提取方法，其特征在于所述的步骤②中U的K个类中心的获取过程为：

②-1、随机选取U中的K条输入信号，将随机选取的每条输入信号作为U的一个类中心，其中，K∈[1,N]；

②-2、将U中当前待处理的第i条输入信号定义为当前输入信号，其中，i的初始值为1，1≤i≤N；

②-3、计算当前输入信号与U的各个类中心之间的欧式距离，然后从当前输入信号对应的K个欧式距离中找出最小欧式距离，假设最小欧式距离对应的类中心为U的第k个类中心，则确定当前输入信号属于U的第k个类中心，其中，1≤k≤K；

②-4、令i＝i+1，将U中下一个待处理的输入信号作为当前输入信号，然后返回步骤②-3继续执行，直至U中所有的输入信号处理完毕，再执行步骤②-5，其中，i＝i+1中的“＝”为赋值符号；

②-5、对属于U的任一个类中心的所有输入信号求均值，然后用该类中心对应的均值替换该类中心，实现类中心的更新；

②-6、重复执行步骤②-2至步骤②-5多次，获得U的K个最终的类中心。

3.根据权利要求1或2所述的基于径向基函数神经网络的单位脉冲响应函数提取方法，其特征在于所述的步骤④中取λ＝1。

4.根据权利要求3所述的基于径向基函数神经网络的单位脉冲响应函数提取方法，其特征在于所述的步骤⑤中取G_n＝0.5K。

5.根据权利要求4所述的基于径向基函数神经网络的单位脉冲响应函数提取方法，其特征在于所述的步骤⑤的具体过程为：

⑤-1、根据径向基函数神经网络的特点，获知Ψ、隐层权值矩阵W与Y'之间存在线性关系：W×Ψ＝Y'，然后对W×Ψ＝Y'两边进行转置，得到Ψ^T×W^T＝Y'^T，其中，W表示径向基函数神经网络的隐层权值矩阵，W的维数为M×K，M表示Y中的每条输出信号的维数；

⑤-2、将W的转置W^T稀疏表示为W^T＝D1×S×D2，其中，D1和D2为两个不同维度方向的二维离散余弦变换基，S表示W^T在二维离散余弦变换下的稀疏系数矩阵；

⑤-3、根据Ψ^T×W^T＝Y'^T和W^T＝D1×S×D2，得到Ψ^T×D1×S×D2＝Y'^T，进一步得到Ψ^T×D1×S＝Y'^T×D2^T，其中，D2^T为D2的转置；然后令A＝Ψ^T×D1，并令B＝Y'^T×D2^T，将Ψ^T×D1×S＝Y'^T×D2^T简化为A×S＝B；

⑤-4、令G_n表示采集度，构造一个维数为G_n×N的随机采样矩阵G，并将G作为稀疏复原过程的观测矩阵；然后利用G对A进行随机采样，得到感知矩阵，记为Ar，Ar＝G×A，并利用G对B进行随机采样，得到由多个观测向量构成的测量矩阵，记为Br，Br＝G×B；

⑤-5、设定S的稀疏度为m1，根据A×S＝B及Ar＝G×A和Br＝G×B，得到Ar×S＝Br；然后利用正交匹配追踪OMP算法求解Ar×S＝Br得到S，其中，1≤m1≤G_n；

⑤-6、根据求解得到的S，由W^T＝D1×S×D2求解W^T，之后对W^T进行转置得到W。

6.根据权利要求5所述的基于径向基函数神经网络的单位脉冲响应函数提取方法，其特征在于所述的步骤⑤-5中利用正交匹配追踪OMP算法求解Ar×S＝Br得到S的具体过程为：

⑤-5a、以ks表示迭代次数，令pp表示维数为m1的向量，令B1表示维数为G_n×K的矩阵，其中，ks的初始值为1，1≤ks≤m1；

⑤-5b、计算Ar中的每一列与第ks-1次迭代的残差值r_ks-1之间的相关系数，共得到K个相关系数；然后将K个相关系数中的最大值对应的Ar中的一列的列号作为pp中的第ks个元素的值，并将K个相关系数中的最大值对应的Ar中的一列作为B1中的第ks列；再在Ar中将K个相关系数中的最大值对应的Ar中的一列剔除；

其中，Ar中的第j列与第ks-1次迭代的残差值r_ks-1之间的相关系数为Ar中的第j列与第ks-1次迭代的残差值r_ks-1之间的内积，1≤j≤K，当ks＝1时r_ks-1为Br中的第1列观测向量；

⑤-5c、利用最小二乘法，计算r_ks-1在B1的每一列上的稀疏投影系数，得到稀疏投影系数序列，记为a_ks；

⑤-5d、根据B1和a_ks，计算第ks次迭代的残差值，记为r_ks，r_ks＝r_ks-1-B1×a_ks；

⑤-5e、判断ks是否等于m1，如果是，则得到S，S中的非零值的位置为pp中的值，相应位置所对应的非零值为a_ks中的值，S中的其它位置的值均为零；否则，令ks＝ks+1，然后返回步骤⑤-5b继续执行，其中，ks＝ks+1中的“＝”为赋值符号。

7.根据权利要求6所述的基于径向基函数神经网络的单位脉冲响应函数提取方法，其特征在于所述的步骤⑥中取λ＝1。

8.根据权利要求7所述的基于径向基函数神经网络的单位脉冲响应函数提取方法，其特征在于所述的步骤⑤中D1的构建过程为：

a1、令D1表示一个K阶方阵，令g和h为两个变量，1≤g≤K,1≤h≤K；

a2、将D1中第g行第h列的元素记为D1(g,h)， 

所述的步骤⑤中D2的构建过程为：

b1、令D2表示一个M阶方阵，令g'和h'为两个变量，1≤g'≤M,1≤h'≤M；

b2、将D2中第g'行第h'列的元素记为D2(g',h')，