本发明涉及一种超声回波信号预处理方法，特别是一种基于奇异谱分析的超声回波信号在线预处理方法。
背景技术：
：超声检测技术利用超声换能器发射超声脉冲与被测对象作用，通过接收并研究反向散射的超声回波信号，可确定被测目标的形状尺寸、力学结构等物理特性，属于非侵入式检测技术，被广泛用于设备故障诊断、无损探伤等领域。超声回波信号属于典型的时变、非线性、非平稳信号，良好的信号预处理对于其后续检测意义重大，特别是对于在线检测设备而言。具体来说，一是因为在线检测设备往往需要兼顾实用性和便捷性，其系统结构较之实验级的精密设备已被简化，测量中往往容易引入噪声而影响准确性和精度；二是由于检测现场噪声源的不确定性，超声回波在传播过程中可能会受到多种因素的影响。综上可知，超声在线检测系统对于回波信号的预处理方法提出了较高要求，其既要能适用于超声信号非平稳、非线性的特点，并且对于不同特性的噪声应具备较好的鲁棒性，去噪时不能改变一些重要的波形特征，如峰位、峰宽、波形等，另外还应避免选用那些过于复杂、难以现场实现的算法，尽量减少算法中需要人工干预的部分。目前对于超声回波信号的数据预处理方法很多，切入点也各自有所不同。近年来，诸如小波变换、小波包变换、经验模态分解等时-频域分析信号方法在超声回波信号去噪、特征提取等方面应用广泛，但对于小波变换、小波包变换的最优小波基选择、最佳分解层数和阈值的确定，以及经验模态分解的固有模态函数筛选、端点效应和冗余模态处理等问题尚有待进一步研究，目前大都仍依赖实验加以确定，同时这些方法也存在复杂度高、计算速度慢的问题。奇异谱分析技术(singularspectrumanalysis，SSA)是一种非线性信号处理方法，较之其他时-频域信号处理技术，如小波变换、小波包变换，它具有以下显著优点：首先，SSA是一种非参数的信号处理方法，它无需对信号和噪声特性做出先验假设，因此即使是对非线性、非平稳的信号也能获得较好的效果，具备滤波频域范围的普适性。其次，SSA是一种数据驱动的信号分析方法，可将一维超声信号序列依据其自身的局部特征信息自适应的分解为若干相互独立的子序列(SSA主成分重构序列)的加和形式，分解所得的子序列则表征了原始信号的局部特征信息，可能具有一定的物理意义，如趋势分量、周期或类周期分量、噪声分量。但与小波变换和小波包变换不同，SSA中用于信号分解-重构运算中的基函数是从原始信号中直接生成，不再需要人为的预先定或选择，可适应数据在线预处理的需求。第三，SSA是能够有效地提取隐含在超声信号中的频域特性(周期性分量、类周期性分量)，可以视为一个广义频域滤波器，与傅里叶分析相比，基于SSA的周期性特征识别/提取方法具有更强的稳定性和识别能力。第四，SSA自1986年提出以来已经过国内外学者的多次优化，在算法计算量和运行实时性方面都大为改善，目前已成功用于在线数据处理方面并显示出其优越性，例如电气设备检测信号自适应去噪、胎儿心电图实时分析、机械振动信号在线处理等。技术实现要素：发明目的：本发明所解决的技术问题在于提供一种超声回波数据在线预处理方法，能够自动实现回波信号的噪声去除和不同频段信号分量的分离提取，改善超声回波测量信号的信噪比、提高在线设备的检测能力。实现本发明目的的技术解决方案为：一种基于奇异谱分析的超声回波信号在线预处理方法，包括以下步骤：步骤1：构造超声回波信号轨道矩阵；步骤2：对超声回波信号轨道矩阵X进行奇异值分解，获得超声回波信号轨道矩阵的奇异值谱；步骤3：求得超声回波信号轨道矩阵X的特征值谱和奇异值差分谱；步骤4：根据超声回波信号轨道矩阵的奇异值谱、特征值谱和奇异值差分谱，采用自适应重构算法确定奇异谱分析的重构分组数r，根据前r个奇异谱分析SSA主成分重构序列得到去噪后的超声回波信号序列；步骤5：前r个奇异谱分析SSA主成分重构序列，分别对应超声回波信号中不同频段下的信号分量，对于每个奇异谱分析SSA主成分重构序列进行傅里叶变换FFT(FastFourierTransformation)得到其对应的频谱，从而完成超声回波信号在线预处理。步骤1包括：对于{t1,t2,...,tN}时刻对应的实测超声回波信号序列xsignal＝{x(t1),x(t2),...,x(tN)}＝{x1,x2,...,xN}，构造如下L×K阶Hankel型的超声回波信号轨道矩阵X：X=x(t1)x(t2)...x(tK)x(t2)x(t3)...x(tK+1)............x(tL)x(tL+1)...x(tN)=x1x2...xKx2x3...xK+1............xLxL+1...xN,]]>其中，L、K分别为矩阵X的行数和列数，x(tN)＝xN表示tN时刻对应的实测超声回波信号，L也被称为轨道矩阵嵌入窗口长度，L为正整数且满足L≤K≤N/2，采用嵌入窗口长度选择算法确定轨道矩阵嵌入窗口长度L，K＝N-L+1，N为超声回波信号序列长度，{t1,t2,...,tN}为超声回波采样时刻序列。步骤1中所述采用嵌入窗口长度选择算法确定轨道矩阵嵌入窗口长度L，包括如下步骤：步骤1-1，预设一个超声回波信号序列的长度门限值N0；步骤1-2，当实测超声回波信号序列xsignal＝{x1,x2,...,xN}的长度N≤N0时，直接指定轨道矩阵X的嵌入窗口长度当实测超声回波信号序列xsignal的长度N>N0时，执行步骤1-3，其中符号表示向下取整；步骤1-3，根据经验公式确定轨道矩阵嵌入窗口长度的下限L0＝(lnN)c，其中参数c∈[1.5,2.5]；步骤1-4，通过如下公式，计算当嵌入窗口长度值L取不同值时，超声回波信号序列xsignal的自相关测试系数Corr(L)：Corr(L)=Σi=1N(xi+L*-xave)(xi-xave)Σi=1N(xi-xave)2,]]>依次选取不同的L，即得到的Corr(L)实际是关于L的函数，其中，xave为实测的超声回波信号序列xsignal的算术平均值，变量满足：当1≤i≤N-L时当N-L<i≤N时xi、xi+L和xi+L-N分别表示超声回波信号序列xsignal＝{x1,x2,...,xN}中第i个元素、第i+L个元素和第i+L-N个元素，i＝1,2,...,N；步骤1-5，基于已确定的轨道矩阵嵌入窗口长度的下限L0，选取满足L>L0且自相关测试系数Corr(L)曲线第一次穿越零点时对应的嵌入窗口长度集合{L1,L2}，即有Corr(L1)Corr(L2)≤0且L2＝L1+1>L0；步骤1-6，通过如下公式确定轨道矩阵X的嵌入窗口长度L，：步骤2包括：采用如下公式对超声回波信号轨道矩阵X进行奇异值分解：X=UΛVT=Σj=1LσjujvjT,]]>其中矩阵U和矩阵V为K×L阶正交矩阵，U＝[u1,u2,...,uL]，V＝[v1,v2,...,vL]，uj和vj分别表示矩阵U和V的第j个列向量，对角阵Λ＝diag(σ1,σ2,...,σL)为L×L阶方阵，其主对角元σ1,σ2,...,σL为超声回波信号轨道矩阵X的L个奇异值，序列{σ1,σ2,...,σL}即超声回波信号矩阵X的奇异值谱，且满足σ1≥σ2≥...≥σL，σj表示奇异值谱的第j个元素，j＝1,2,...,L。在计算方式上，向量u1,u2,...,uL是矩阵XXT的L个特征向量，L个奇异值的平方是轨道矩阵XXT的L个非零特征值，对于向量v1,v2,...,vL有vj＝XTuj/σj，j＝1,2,...,L；步骤3中，根据步骤2获得的超声回波信号轨道矩阵X的L个奇异值，求得超声回波信号轨道矩阵X的特征值谱{δ1,δ2,...,δj,...,δL}和奇异值差分谱{d1,d2,...,dp,...,dL-1}，δj表示超声回波信号轨道矩阵X的第j个特征值，j＝1,2,...,L，dp表示轨道矩阵X奇异值差分谱的第p个元素，p＝1,2,...,L-1，其中dp＝d(σp)＝σp-σp+1，σj、σp和σp+1分别表示矩阵X奇异值谱的第j个元素、第p个元素和第p+1个元素。步骤4包括如下步骤：步骤4-1，通过如下方式构建表示奇异值差分谱{d1,d2,...,dp,...,dL-1}中所有峰值点的集合Ψ：对于奇异值差分谱中的元素dp，若dp>dp-1且dp>dp+1则令dp∈Ψ，反之则令再通过如下方式构建阈值候选集Φ：对于奇异值差分谱中的元素dp，若dp∈Ψ且dp＝(σp-σp+1)>dave则令其对应的奇异值σp∈Φ，反之则令阈值候选集Φ表示集合Ψ中大于dave的点所对应的奇异值集合，其中p＝1,2,...,L-1，dave为判定奇异值差分谱中一个峰是否显著的一个经验参数，dave定义为奇异值差分谱中所有峰值的算术平均，即dave＝mean(Ψ)；步骤4-2，通过如下方式构建表示特征值谱{δ1,δ2,...,δj,...,δL}中大于δave的点所对应的奇异值集合Ω：对于特征值谱中的元素δj，若则令其对应的奇异值σj∈Ω，反之则令其中j＝1,2,...,L，δave为判定一个奇异谱分析SSA主成分是否显著的一个经验参数，δave定义为所有特征值的算术平均，即步骤4-3，选取集合Ω中的最小元素σs＝argmin(Ω)作为信号、噪声的近似边界；步骤4-4，若σs∈Φ则直接将σs确定为最终的去噪阈值参量σr，否则找出集合Φ中最接近σs的元素并将其作为最终的去噪阈值参量σr，r表示重构分组数；步骤4-5，通过对角平均化分别将前r个奇异谱分析SSA主成分轨道矩阵{X(1),X(2),...,X(q),...,X(r)}转化为对应的重构序列{g(1),g(2),...,g(q),...,g(r)}，其中序列过程如下公式所示：gi(q)=1i+1Σm=1i+1xm,i-m+2(q),i∈[1,L]1LΣm=1Lxm,i-m+2(q),i∈[L,K]1N-iΣm=i-K+2N-K+1xm,i-m+2(q),i∈[K,N],]]>其中，X(q)表示第q个奇异谱分析SSA主成分轨道矩阵，为矩阵X(q)中第m行第n列元素，表示重构序列g(q)的第i个元素，q＝1,2,...,r，i＝1,2,...,N；步骤4-6，确定重构分组数r后，将前r个奇异谱分析SSA主成分重构序列{g(1),g(2),...,g(r)}进行加合，得到去噪后的重构超声回波信号序列其中ysignal＝{y(t1),y(t2),...,y(tN)}＝{y1,y2,...,yN}，g(q)表示第q个奇异谱分析SSA主成分重构序列，有益效果：本发明与现有技术相比，其显著优点为：该方法集成了奇异谱分析(SSA)技术非参数、数据驱动的优点，对于信号平稳性、噪声特性没有特别要求，同时具备自适应的参数选取策略，使得数据处理过程中无需人工干预，满足超声回波信号在线预处理的使用要求。本发明利用奇异谱分析技术，可自动实现超声检测中回波信号的噪声去除和不同频段信号分量的分离提取，具有对信号平稳性和噪声特性无特殊要求、数据处理过程中无需人工干预等优点，能够满足超声回波信号在线预处理的使用要求。附图说明下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述和/或其他方面的优点将会变得更加清楚。图1为本发明的一种基于奇异谱分析的超声回波信号在线预处理方法的流程图。图2为本发明中超声回波信号轨道矩阵的嵌入窗口长度选择算法的流程图。图3为本发明中超声回波信号自适应重构算法的流程图。图4a为实施实例中仿真生成的不含噪的超声回波信号时域信息。图4b为实施实例中仿真生成的不含噪的超声回波信号频域信息。图5a为信噪比SNR＝8时实施实例中仿真生成的加噪超声回波信号。图5b为信噪比SNR＝15时实施实例中仿真生成的加噪超声回波信号。图6是本发明实施实例中嵌入窗口长度选择算法的自相关测试系数Corr(L)。图7a是实施实例中SNR＝8情况下的奇异值谱和奇异值差分谱。图7b是实施实例中SNR＝8情况下的特征值谱。图8a为信噪比SNR＝8时实施实例中超声回波信号去噪前后的对比效果。图8b为信噪比SNR＝15时实施实例中超声回波信号去噪前后的对比效果。图9a为实施实例中SNR＝8情况下第1-2个、第3-4个SSA主成分重构序列时域信息。图9b为实施实例中SNR＝8情况下第1-2个、第3-4个SSA主成分重构序列频域信息。具体实施方式结合图1本发明的一种基于奇异谱分析的超声回波信号在线预处理方法，包括以下步骤：步骤1：超声回波信号轨道矩阵构造，假设{t1,t2,...,tN}时刻对于超声回波信号的采样序列为xsignal＝{x(t1),x(t2),...,x(tN)}＝{x1,x2,...,xN}，将其构造为L×K阶Hankel型超声回波信号轨道矩阵X，X=x(t1)x(t2)...x(tK)x(t2)x(t3)...x(tK+1)............x(tL)x(tL+1)...x(tN)=x1x2...xKx2x3...xK+1............xLxL+1...xN]]>其中，L、K分别为矩阵X的行数和列数，L也被称为轨道矩阵嵌入窗口长度，L为正整数且满足L≤K≤N/2，采用嵌入窗口长度选择算法确定，K＝N-L+1，N为超声回波序列长度；步骤2：对超声回波信号轨道矩阵X进行奇异值分解其中矩阵U和矩阵V为K×L阶正交矩阵，U＝[u1,u2,...,uL]，V＝[v1,v2,...,vL]，uj和vj分别表示矩阵U和V的第j个列向量，对角阵Λ＝diag(σ1,σ2,...,σL)为L×L阶方阵，其主对角元σ1,σ2,...,σL为轨道矩阵X的L个奇异值，序列{σ1,σ2,...,σL}也即矩阵X的奇异值谱，且满足σ1≥σ2≥...≥σL。在计算方式上，向量u1,u2,...,uL是矩阵XXT的L个特征向量，L个奇异值的平方是轨道矩阵XXT的L个非零特征值，对于向量v1,v2,...,vL有vj＝XTuj/σj，j＝1,2,...,L；步骤3：求得超声回波信号轨道矩阵X的特征值谱{δ1,δ2,...,δj,...,δL}和奇异值差分谱{d1,d2,...,dp,...,dL-1}，其中dp＝d(σp)＝σp-σp+1，j＝1,2,...,L，p＝1,2,...,L-1，，σj、σp和σp+1分别表示矩阵X奇异值谱的第j个元素、第p个元素和第p+1个元素；步骤4：根据矩阵X的奇异值谱、特征值谱和奇异值差分谱，采用自适应重构算法确定重构分组数r，则去噪后的超声回波序列ysignal为前r个SSA主成分重构序列的加和其中ysignal＝{y(t1),y(t2),...,y(tN)}＝{y1,y2,...,yN}，为第q个SSA主成分重构序列，q＝1,2,...,r；步骤5：前r个SSA主成分重构序列{g(1),g(2),...,g(r)}分别对应超声回波信号ysignal中不同频段下的信号分量，且g(q)对于ysignal的方差贡献率依次递减，q＝1,2,...,r，对g(q)进行傅里叶变换(FFT，FastFourierTransformation)可得到其对应的频谱。结合图2，本发明的嵌入窗口长度选择算法的具体步骤为：(1)根据系统运算能力，人工预设一个超声回波序列的长度门限值N0；(2)当实测超声回波序列xsignal＝{x1,x2,...,xN}的长度N≤N0时，直接指定轨道矩阵X的嵌入窗口长度其中符号表示向下取整；(3)当实测超声回波序列xsignal的长度N>N0时，为避免因构造的轨道矩阵X维数过大而造成的运算复杂度增加，需选取一个适当小的嵌入窗口长度L，为此首先根据经验公式确定嵌入窗口长度的下限L0＝(lnN)c，其中参数c∈[1.5,2.5]；(4)计算xsignal的自相关测试系数该参数表征了以不同的嵌入窗口长度值L构造轨道矩阵X时各SSA主成分间的相关性(可分离性)，其中xave为序列xsignal的算术平均值变量满足：当1≤i≤N-L时当N-L<i≤N时xi、xi+L和xi+L-N分别表示超声回波信号序列xsignal＝{x1,x2,...,xN}中第i个、第i+L个和第i+L-N个元素，i＝1,2,...,N。特别的，当L的取值满足使Corr(L)趋于零时，说明奇异值分解后各SSA主成分间的相关性很小(可分离性很大)，奇异值分解的效果较好，也更利于进行后续的信号分析和预处理；(5)基于已确定的轨道矩阵嵌入窗口长度的下限L0，选取满足L>L0且Corr(L)曲线第一次穿越零点时对应的嵌入窗口长度集合{L1,L2}，即Corr(L1)Corr(L2)≤0且L2＝L1+1>L0；(6)通过如下公式确定轨道矩阵X的嵌入窗口长度L：：结合图3，本发明的自适应重构算法的具体步骤为：(1)通过如下方式构建集合Ψ，表示奇异值差分谱{d1,d2,...,dp,...,dL-1}中所有峰值点的集合：对于奇异值差分谱中的元素dp，若dp>dp-1且dp>dp+1则令dp∈Ψ，反之则令随后通过如下方式构建阈值候选集Φ，表示Ψ中大于dave的点所对应的奇异值集合：对于奇异值差分谱中的元素dp，若dp∈Ψ且dp＝(σp-σp+1)>dave则令其对应的奇异值σp∈Φ，反之则令其中p＝1,2,...,L-1，dave为判定奇异值差分谱中某个峰是否显著的一个经验参数，定义为奇异值差分谱中所有峰值的算术平均，即dave＝mean(Ψ)；(2)通过如下方式构建集合Ω，表示特征值谱{δ1,δ2,...,δL}中大于δave的点所对应的奇异值集合：对于特征值谱中的元素δj，若则令其对应的奇异值σj∈Ω，反之则令其中j＝1,2,...,L，δave为判定某个SSA主成分是否显著的一个经验参数，定义为所有特征值的算术平均，即(3)选取集合Ω中的最小元素σs＝argmin(Ω)作为信号、噪声的近似边界；(4)若σs∈Φ则直接将σs确定为最终的去噪阈值参量σr，否则找出集合Φ中最接近σs的元素并将其作为σr，r表示重构分组数；(5)确定重构分组数r后，通过对角平均化(diagonalaveraging)，分别将前r个SSA主成分轨道矩阵{X(1),X(2),...,X(q),...,X(r)}转化为对应的重构序列{g(1),g(2),...,g(q),...,g(r)}，其中序列该过程可表述为gi(q)=1i+1Σm=1i+1xm,i-m+2(q),i∈[1,L]1LΣm=1Lxm,i-m+2(q),i∈[L,K]1N-iΣm=i-K+2N-K+1xm,i-m+2(q),i∈[K,N]]]>其中为第q个SSA主成分轨道矩阵，为矩阵X(q)中第m行第n列元素，表示重构序列g(q)的第i个元素，q＝1,2,...,r，i＝1,2,...,N；(6)将前r个SSA主成分重构序列{g(1),g(2),...,g(r)}进行加合，得到去噪后的重构超声回波信号序列下面结合实施例对本发明做进一步详细的描述：结合一个对于超声回波仿真信号的预处理过程，对本发明作进一步详细描述。对于一个超声检测系统，根据超声回波的物理特性，超声换能器接收到的含噪回波信号模型可表征为x(t)＝s(t)+n(t)其中，t表示采样时间，x(t)表示实测的含噪超声回波信号序列，n(t)表示加性高斯噪声,s(t)表示回波中不含噪的信号分量，且满足以下形式s(t)=βe-α(t-τ)2cos[2πfc(t-τ)+φ]]]>其中，α为回波带宽系数，τ为回波到达时间，fc为回波中心频率，φ为回波相位，β为回波幅度系数。仿真实验中，取x(t)＝s1(t)+s2(t)+n(t)，其中α1＝α2＝1.5(MHz)2，τ1＝2μs，τ2＝7μs，fc,1＝fc,2＝4MHz，φ1＝φ2＝1.8rad，β1＝0.5，β2＝0.2。检测系统对于超声回波信号采样时长(0-10)μs，回波采样序列长度N＝500，采样频率50MHz。加性噪声n(t)特性为高斯白噪声，分别测试本发明对于信噪比SNR＝8和SNR＝15两种情况下的预处理效果。不含噪的超声回波信号s(t)＝s1(t)+s2(t)的时域和频域信息分别如图4a和图4b所示，两种加噪条件下的超声回波信号x(t)分别如图5a和图5b所示。关于预设的模型参数，本实施例中嵌入窗口长度选择算法取序列长度门限值N0＝200，嵌入窗口长度下限的经验公式中参数c＝2。需要说明的是，上述参数选取对于本发明的预处理方法来说是鲁棒的，即参数N0和c的选择不会影响最终回波信号预处理的效果。采用图2所示的嵌入窗口长度选择算法，考虑到回波序列长度N>N0，计算出嵌入窗口长度下限L0＝38，自相关测试系数Corr(L)如图6所示，经算法确定L＝41。据此，分别构造两种不同SNR时的超声回波信号轨道矩阵。随后对超声回波信号轨道矩阵进行奇异值分解，以SNR＝8情况为例，轨道矩阵对应的奇异值谱、奇异值差分谱和特征值谱如图7a和图7b所示。采用图3所示自适应重构算法，两种SNR条件下的最佳重构分组数r均为4，进而可确定去噪后的重构超声回波信号序列ysignal。图8a和图8b分别显示了两种SNR条件下超声回波信号去噪前后的对比效果。将图8a和图8b与图4a、图4b对比，可以看出算法的去噪效果较为理想，去噪前后回波信号的波形特性(如波峰、波谷的位置、强度等特种)得到了较好的保留。仍以SNR＝8情况为例，进一步分析预处理算法对于回波信号中各频域分量的提取效果。结合图7a和图7b可知，在第2、第4个奇异值处奇异值谱的幅值存在较大突变，对应奇异值差分谱则有两个较大的峰。鉴于信号、噪声分量相关性的差异会导致其对应奇异值幅值的不同，在奇异值谱中表征信号、噪声的边界附近应存在明显的趋势变化现象(急速降低)，反映在奇异值差分谱中则表现为较明显的峰。因此，分析后认为：分解后的SSA主成分重构序列大致可分为三类：第1-2个SSA主成分、第3-4个SSA主成分，其余的SSA主成分，其对于ysignal的方差贡献率η如表1所示。可见，前4个SSA主成分的方差贡献率之和已经占到88.96％，可以基本表征信号分量所含的特征信息。表1各SSA主成分重构序列对于ysignal的方差贡献率SSA主成分数12345678其余方差贡献率(％)42.9039.393.762.920.410.410.410.409.39分别比较图4a和图9a、图4b和图9b，第1-2个SSA主成分重构序列已能基本表征信号s(t)的时域、频域特征，第3-4个SSA主成分重构序列可视为对其部分细节特征的补充修正，使重构信号ysignal与s(t)趋于一致。本发明提供了一种基于奇异谱分析的超声回波信号在线预处理方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本
技术领域：
的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。当前第1页1 2 3