专利名称：用于海中垂直阵接收信号互谱矩阵估计中的噪声抑制方法
技术领域：
本发明涉及水声信号处理领域，更具体地说，涉及一种用于海中垂直阵接收信号互谱矩阵估计中的噪声抑制方法。
背景技术：
在海洋中应用的水声系统常使用多阵元垂直水听器阵列(以下简称垂直阵)来探测目标声源的信号。而水声系统在进行信号处理时，通常首先要将垂直阵中各阵元接收到的水声信号构造为一个信号互谱矩阵，并依据该互谱矩阵展开进一步的信号处理。为了抑制接收信号中的噪声，需要在信号互谱矩阵估计时采用一定的噪声抑制方法。
如图1示出了垂直阵探测目标声源信号的工作示意图，垂直阵1的阵元数为M，zs为目标声源2在海水中的深度(以下简称为声源深度)，也就是目标声源2距离海面3的距离。r为目标声源2与垂直阵1的间距(以下简称接收距离)，该距离r会随目标声源2的运动而变化。第m阵元4为垂直阵1中的一个示例性阵元，zm为第m阵元4的深度。用标号5表示海底。
将垂直阵1接收的水声信号按采样时间划分为长度为T的多个子段信号，将垂直阵1中的第m阵元在tn(＝t0+nT)时刻采样的子段信号记为sr(zm，t0+nT+t) m＝1......M、 n＝-N......N、0＜t＜T。(1)其中t0为某一采样时刻，用tn表示t0前后的第n个采样时刻。根据海洋声传播理论，公式(1)进一步具有如下的形式sr(zm，t0+nT+t)＝s0(t0+nT+t)*h(r(t0+nT)，zs，zm，r)+n(zm，t0+nT+t)(2)其中，s0(t0+nT+t)为目标声源信号(为了论述简单，假定只有一个目标)；*为卷积运算；h(r(t0+nT)，zs，zm，t)为信道冲击响应函数，与声源深度zs、接收距离r和接收深度zm有关，且假定接收距离r在时间T内不变；n(zm，t0+nT+t)为背景噪声(一般为加性高斯白噪声)。
将公式(2)变换到频域，有Sr(zm，t0+nT，ω)＝S0(t0+nT，ω)H(r(t0+nT)，zs，zm，ω)+N(ω，zm，t0+nT)，(3)其中，S0(t0+nT,ω)=∫0Ts0(t0+nT+t)exp(-jωt)dt,----(4)]]>H(r(t0+nT),zs,zm,ω)=∫0Th(r(t0+nt),zs,zm,t)exp(-jωt)dt,----(5)]]>N(zm,t0+nT,ω)=∫0Tn(zm,t0+nT+t)exp(-jωt)dt.----(6)]]>进一步计算垂直阵1的各阵元在tn时刻接收到的子段信号的互谱矩阵(M×M阶)的测量值R(t0+nT，ω)，由于在计算中只取各阵元在一个时刻接收到的一个子段信号，因此R(t0+nT，ω)又可称为互谱矩阵的单样本测量值。矩阵R(t0+nT，ω)的(m，k)元素为垂直阵中的第m和第k阵元接收信号频谱的共轭乘积，如公式(7)所示Rm,k(t0+nT,ω)=Sr(zm,t0+nT,ω)Sr*(zk,t0+nT,ω),m,k=1......M----(7)]]>将公式(3)代入公式(7)有，Rm,k(t0+nT,ω)=|S0(t0+nT,ω)|2H(r(t0+nT),zs,zm,ω)H*(r(t0+nT),zs,zk,ω)---(8)]]>+N~m,k(t0+nT,ω)]]>其中N~m,k(t0+nT,ω)=S0(t0+nT,ω)H(r(t0+nT),zs,zm,ω)N*(zk,t0+nT,ω)]]>+S0*(t0+nT,ω)H*(r(t0+nT),zs,zk,ω)N(zm,t0+nT,ω).----(9)]]>+N(zm,t0+nT,ω)N*(zk,t0+nT,ω)]]>公式(8)等号右边的第一项，为目标声源信号经过信道传播后的有效信号信道分量的垂直阵互谱； 为互谱噪声项，可近似视为白噪声，其期望值(时间平均值)为0，即E[N~m,k(t0+nT,ω)]=0.----(10)]]>在常规互谱矩阵统计中，有以下两点假设
(1)在统计时间内，公式(8)中有效信号信道分量的传递函数H的互谱位相不变，即H(r(t0+nT)，zs，zm，ω)H*(r(t0+nT)，zs，zk，ω)＝H(r(t0)，zs，zm，ω)H*(r(t0)，zs，zk，ω)。
(11)(2)目标声源功率谱在统计时间内平稳，存在期望值(时间平均值)E(|S0(t0+nT，ω)|2)＝|S0(ω)|2， (12)于是Rm，k(t0+nT，ω)的期望值为E[Rm，k(t0+nT，ω)]＝|S0(ω)|2H(r(t0)，zs，zm，ω)H*(r(t0)，zs，zk，ω)。
(13)基于上述的两点假设，现有技术中为了抑制如公式(8)所示的互谱矩阵R(t0+nT，ω)中的噪声成份，通常采用多样本统计的方法，即对t0时刻前后共2N+1个互谱矩阵单样本测量值R(t0+nT，ω)做统计平均，得到t0时刻互谱矩阵的估计值 其各个元素为R^m,k(t0,ω)=12N+1Σn=-NNRm,k(t0+nT,ω),m,k=1......M----(14)]]>根据现有的信号处理理论，较之单样本测量值，经过2N+1个样本统计后， 信噪比的增加量近似为10log10(2N+1)分贝。用于统计的样本数越多，噪声抑制效果越好，得到的 信噪比就越高。所以为了获得较高的信噪比，需要尽量增加统计样本的个数。
但是在实际应用中，用于统计的时间样本个数2N+1是有限的。这是因为当目标声源运动时，接收距离r随时间变化，统计时间过长时，公式(11)所示的假设将得不到满足。此时，长时间统计不但不能提高 的信噪比，反而可能使其信噪比下降。
对水声目标声源，尤其是水面舰船等，通常可假定其在较长时间内(1分钟以上)相对于垂直阵只做匀速水平运动，而深度基本不变。于是接收距离r＝r0+v(t-t0)，v为目标声源相对于垂直阵的径向速度。定义使下式成立的δr(r，zs，zm，zk，ω)为Rm，k(r，zs，ω)的水平相关半径(超过水平相关半径，Rm，k(r，zs，ω)相位变化超过 )，则Re{[H(r,zs,zm,ω)H*(r,zs,zk,ω)][(r+δ,zs,zm,ω)H*(r+δ,zs,zk,ω)]*}|H(r,zs,zm,ω)H*(r,zs,zk,ω)H(r+δ,zs,zm,ω)H*(r+δ,zs,zk,ω)|≥22,]]>|δ|＜δr(r，zs，zm，zk，ω)(15)对于特定目标接收信号的垂直阵空间互谱矩阵R(r，zs，ω)，其水平相关半径可定义为δ~r(r,zs,ω)=minm,k(δr(r,zs,zm,zk,ω))----(16)]]>若统计时间内目标声源的距离变化小于 公式(11)近似成立。故而可进行接收目标信号垂直阵空间互谱矩阵统计的时间长度不能超过 以30米深、泥沙底质、冬季水文浅海信道为例，假定声源位于10米深度，垂直阵从海面到海底间隔2米分布。声场模型数值计算表明，在1kHz以下的低频段，在接收距离较远时(大于1公里)， 约为100米。若目标航速为9米/秒(约18节)，垂直阵互谱矩阵常规统计方法的统计时间不能超过11秒。
综上所述，现有技术中用统计的方法对垂直阵接收信号互谱矩阵估计进行噪声抑制有如下缺陷首先，受到互谱矩阵水平相关半径的限制，统计时间一般较小，限制了信噪比的提高潜力；其次，实际应用中，由于缺少先验知识，该统计时间限制很难加以预报；一旦统计时间过长，会导致统计得到的互谱矩阵信噪比严重下降。
因此就需要有一种对互谱矩阵进行噪声抑制的方法，该方法不受水平相关半径的限制，增长统计时间，从而进一步抑制互谱矩阵估计中的噪声成份，提高信号的信噪比。
发明内容
本发明的目的在于克服现有技术在对垂直阵接收信号的互谱矩阵估计进行噪声抑制方法的不足，提供一种新的噪声抑制方法，该方法相比于现有技术可具有更长的统计时间，以便进一步抑制互谱矩阵估计中的噪声成份。
本发明人通过理论计算和实验发现，垂直阵空间互谱矩阵各个元素中有效信号信道分量(即公式(8)中右边的第一项)的传递函数互谱[H(r，zs，zm，ω)H*(r，zs，zk，ω)]在距离—频率空间具有稳定的结构，如下式所示[H(r,zs,zm,ω)H*(r,zs,zk,ω)]=]]>[H(r+δr,zs,zm,ω(1+βδrr))H*(r+δr,zs,zk,ω(1+βδrr))]----(17)]]>从公式(17)可看出，在接收距离r+δr处、频率 的垂直阵空间互谱矩阵中各个元素信道分量的传递函数互谱与接收距离r、频率ω的互谱矩阵相同，其中β为一个常数。这就是说，通过适当的频移可以补偿目标距离变化引起的垂直阵互谱矩阵中信道分量的变化。而且，即使δr大于公式(16)中定义的水平相关半径 公式(17)依然满足。
当统计时间较长时，公式(11)的条件难以满足，但公式(17)却可依然满足。定义一个目标运动影响因子α=βδrr,]]>称α为目标运动影响因子。利用公式(17)，公式(14)可以改写为R^m,k(t0,ω)=12N+1Σn=-NNRm,k(t0+nT,ω(1+α))----(18)]]>根据公式(18)，本发明提供一种用于海中垂直阵接收信号互谱矩阵估计中的噪声抑制方法，包括所述垂直阵在第一采样时刻以及与第一采样时刻相邻的多个采样时刻接收到目标声源发出的信号；用第一采样时刻垂直阵的接收信号构造出第一互谱矩阵，用与第一采样时刻相邻的多个采样时刻垂直阵的接收信号相应构造出多个第二互谱矩阵；将所述多个第二互谱矩阵相对于第一互谱矩阵作频移，使得每一个第二互谱矩阵中的有效信号信道分量都与第一互谱矩阵中的有效信号信道分量相等；将第一互谱矩阵与多个频移后第二互谱矩阵作统计平均得到第一采样时刻的第一互谱矩阵的估计。
其中，将所述多个第二互谱矩阵相对于第一互谱矩阵作频移是用一个目标运动影响因子α对第二互谱矩阵中的频率进行变换，与第一采样时刻相邻的多个采样时刻中的第n个采样时刻对应的第二互谱矩阵的频率ω′＝ω(1+nα)，其中ω为第一互谱矩阵的频率，ω′为第二互谱矩阵频移后的频率；当n为正数时，第n个采样在第一采样时刻之后，当n为负数时，第n个采样在第一采样时刻之前。
所述目标运动影响因子α是通过分析所述第一互谱矩阵和多个互谱矩阵在距离一频率空间的互谱条纹得到的。具体为，计算不同α时，公式(19)的值B(α)=∫ωLωHΣk,m=1k≠m|ΣnRm,k(tn,ω(1+nα))|2dω,----(19)]]>B(α)最大值对应的α值就是所需的目标运动影响因子，其中Rm，k(tn，ω(1+nα))为第n个采样时刻tn对应的第二互谱矩阵中的第m行第n列的元素，当n为零时对应于第一采样时刻；ωH、ωL为起止频率。
采用本发明的方法可以补偿目标声源的运动，公式(17)成立所需的水平相关半径远大于公式(11)，一般是其10倍以上。因此，在对目标声源位于r＝r0处的垂直阵互谱统计中，通过对互谱矩阵进行频移，公式(18)可以包括更多的样本，有效地抑制了互谱矩阵中的噪声成份，与常规统计方法相比，信噪比可提高10分贝以上。
图1是垂直阵探测目标声源信号的工作示意图；图2(a)是在一个实施例中垂直阵两阵元的互谱幅度在频率—距离空间的分布；图2(b)是在一个实施例中垂直阵两阵元的互谱相位在频率—距离空间的分布；图3是一个实施例中7米、17米深的两阵元接收信号互谱幅度随时间和频率的变化；图4(a)是B(α)随α和时间的分布；图4(b)目标运动影响因子αop(t)的搜索结果；图5(a)是采用本发明的方法进行目标运动补偿后长时间累积匹配定位结果；图5(b)是采用常规方法进行目标运动补偿后长时间累积匹配定位结果。
具体实施例方式
下面结合附图和具体实施方式
对本发明作进一步详细描述。
(1)垂直阵在第一采样时刻t0以及与第一采样时刻t0相邻的多个采样时刻接收到目标声源发出的信号。在一个实施例中，第一采样时刻t0相邻的多个采样时刻是指与t0相邻的前后各N个采样时刻。这些采样时刻用tn＝t0+nT表示，其中n的取值是从-N到N，当n为负数时，表示该采样时刻在t0之前，当n为正数时，表示该采样时刻在t0之后，当n为零时表示第一采样时刻t0。应该理解，在选定第一采样时刻t0后，其他采样时刻的选取可以有多种变化，并不影响本发明的实施。
将垂直阵M个阵元在时刻t0+nT接收到的长度为T的信号记录为sr(zm，t0+nT+t) m＝1......M、n＝-N......N、0＜t＜T。(20)(2)将垂直阵在各个采样时刻接收到的信号构造成互谱矩阵。
首先对各阵元的接收信号进行傅立叶变化，分析得到其频谱，Sr(zm,t0+nT,ω)=∫0Tsr(zm,t0+nT+t)exp(-jωt)dt.----(21)]]>然后计算各段信号的垂直阵互谱矩阵(M×M阶)单样本测量值R(t0+nT，ω)，矩阵R(t0+nT，ω)的各个元素为Rm,k(t0+nT,ω)=Sr(zm,t0+nT,ω)Sr*(zk,t0+nT,ω),m,k=1......M----(22)]]>其中，第一采样时刻t0的互谱矩阵为R(t0，ω)，称为第一互谱矩阵。其他采样时刻的互谱矩阵称为第二互谱矩阵，共2N个。
(3)相对于第一互谱矩阵对第二互谱矩阵作频移，使得第二互谱矩阵中的有效信号信道分量与第一互谱矩阵中的有效信号信道分量相等。由于第一互谱矩阵和2N个第二互谱矩阵中的有效信号信道分量满足公式(17)，所以可根据公式(19)分析第一互谱矩阵和2N个第二互谱矩阵，估算得到目标运动影响因子α=βδrr.]]>以30米深、泥沙底质、冬季水文浅海信道为例，当目标声源深度10米时，位于5米和15米深的两阵元信道分量的传递函数互谱为[H(r，10，5，ω)H*(r，10，15，ω)]，其在频率ω和距离r空间的幅度分布如图2(a)所示、相位分布如图2(b)所示。显然，图2(a)和图2(b)中的互谱满足公式(17)。
在实际应用中，最好能够通过计算的方式得到公式(17)中的频移参数β或者目标运动影响因子α=βδrr.]]>在一个实施例中，对t0时刻前后共2N+1个互谱矩阵单样本测量值R(t0+nT，ω)做分析，用类似波束形成的方法，搜索目标运动影响因子α的所有可能取值。计算不同α时，B(α)=∫ωLωHΣk,m=1k≠m|Σn=-NNRm,k(t0+nT,ω(1+nα))|2dω,----(23)]]>其中ωH、ωL为起止频率。B(α)最大值对应的α值就是目标运动影响因子的真实值，记为αop。
(4)根据αop，对第一采样时刻t0前后的2N个相邻的第二互谱矩阵R(t0+nT，ω)相对于第一互谱矩阵R(t0，ω)做频移，以补偿目标声源的运动影响，使得第二互谱矩阵的有效信号信道分量与第一互谱矩阵中的有效信号信道分量相等，得到R~m,k(t0+nT,ω)=Rm,k(t0+nT,ω(1+nαop)).----(24)]]>由于在公式(24)中，n为零时没有频移，因此第一互谱矩阵R(t0，ω)也可以归纳到公式(24)中。
(5)对公式(24)中的2N+1个频移的 做统计平均，统计得到t0时刻垂直阵空间互谱矩阵的估计值为R^m,k(t0,ω)=12N+1Σn=-NNR~m,k(t0+nT,ω).----(25)]]>通过该统计平均，抑制了t0时刻垂直阵空间互谱矩阵中的噪声成份。
为助于形象理解具体实施方式
，下面在一个被动匹配定位海试数据处理中应用本发明。在对水面舰船目标的垂直阵被动匹配定位中，对比分析本发明方法和常规(非目标运动补偿)统计方法的结果。
首先简单介绍一下水声匹配场定位，匹配场定位是现代水声定位系统常用的方法。将垂直阵接收目标声源信号互谱的估计值 作为测量值。同时，根据所在海域的海洋环境声参数，可以预报目标在不同假定位置时的垂直阵互谱矩阵理论值Rcopy(r，zs，ω)，其各个元素为Rm,kcopy(r,zs,ω)=H(r,zs,zm,ω)H*(r,zs,zk,ω).----(26)]]>对测量值 和理论值Rcopy(r，zs，ω)加以匹配，以常规BARTLETT匹配方法为例，得到不同目标假定位置(r，zs)的匹配值B(r,zs)=|Σm,k[R^m,k*(r,zs,ω)Rm,kcopy(r,zs,ω)]/Σm,k|R^m,k(r,zs,ω)|2Σm,k|Rm,kcopy(r,zs,ω)|2,----(27)]]>
匹配结果B(r，zs)最大值对应的空间位置(r，zs)即为目标位置估计值。在匹配定位处理中， 的信噪比越高，匹配定位效果越好(匹配结果中目标位置清晰、背景噪声低)。
试验中使用16阵元垂直阵接收目标噪声信号，垂直阵阵元从5米到35米以2米等间距布放。在约25分钟时间内水面舰船目标与垂直阵的间距从14公里变到2公里，相对垂直阵径向航速大小为8米/秒(约16节)。
图3为7米、17米深的两阵元接收信号互谱幅度随时间和频率的变化，图中的条纹通过视觉就可以观察出来。通过用公式(23)分析包括图3的多阵元互谱条纹，得到B(α)随α和时间的分布，如图4(a)所示，图中深颜色代表B(α)大，浅颜色代表B(α)小。搜索每个时刻B(α)最大值位置，得到目标运动影响因子αop(t)，如图4(b)所示。根据目标运动补偿量αop(t)，先对每时刻前后共60个样本的垂直阵互谱矩阵用目标运动影响因子αop(t)作如公式(24)所示的频移，以补偿目标运动。60个样本对应的信号时间长度为84秒，对应的目标径向运动距离672米，远远超过了实验海区的垂直阵互谱水平相关半径(约100米)。对目标运动补偿后的垂直阵互谱矩阵进行统计平均，得到中间时刻的互谱矩阵估计值，最后使用该估计值进行匹配定位。得到的匹配定位结果如图5(a)所示，图中的圆圈“O”为试验中GPS(卫星全球定位系统)记录给出的目标距离参考值，深颜色代表的匹配值高、浅颜色匹配值低。图中的目标匹配位置清晰、匹配结果背景噪声较低，给出的目标距离与GPS参考值基本一致。
作为对比，不补偿目标运动，对每时刻前后共60个样本的空间互谱矩阵直接统计平均，得到其估计值后，进行相同的匹配定位处理。最后的匹配定位结果如图5(b)所示，图中根本找不出目标的明显位置。
从上面的试验数据处理实例可以看出目标运动会导致垂直阵空间互谱矩阵的相关时间存在限制，不补偿目标运动的影响，盲目进行长时间多样本统计，难以在互谱矩阵估计中有效抑制噪声、提高信噪比，甚至会给出错误的互谱矩阵。而根据本发明给出的方法，可以实现目标运动补偿，达到长时间累积垂直阵空间互谱矩阵的目的，最终有效抑制噪声、提高信噪比，改善匹配定位的效果。
权利要求
1.一种用于海中垂直阵接收信号互谱矩阵估计中的噪声抑制方法，包括所述垂直阵在第一采样时刻以及与第一采样时刻相邻的多个采样时刻接收到目标声源发出的信号；用第一采样时刻垂直阵的接收信号构造出第一互谱矩阵，用与第一采样时刻相邻的多个采样时刻垂直阵的接收信号相应构造出多个第二互谱矩阵；将所述多个第二互谱矩阵相对于第一互谱矩阵作频移，使得每一个第二互谱矩阵中的有效信号信道分量都与第一互谱矩阵中的有效信号信道分量相等；将第一互谱矩阵与多个频移后第二互谱矩阵作统计平均得到第一采样时刻的第一互谱矩阵的估计。
2.根据权利要求1所述的噪声抑制方法，其特征在于，将所述多个第二互谱矩阵相对于第一互谱矩阵作频移是指用一个目标运动影响因子α对第二互谱矩阵中的频率进行变换，与第一采样时刻相邻的多个采样时刻中的第n个采样时刻对应的第二互谱矩阵的频率ω′＝ω(1+nα)，其中ω为第一互谱矩阵的频率，ω′为第二互谱矩阵频移后的频率；当n为正数时，第n个采样在第一采样时刻之后，当n为负数时，第n个采样在第一采样时刻之前。
3.根据权利要求2所述的噪声抑制方法，其特征在于，所述目标运动影响因子α是通过分析所述第一互谱矩阵和多个互谱矩阵在距离—频率空间的互谱条纹得到的。
4.根据权利要求2所述的噪声抑制方法，其特征在于，所述目标运动影响因子α通过如下方法得到计算不同α时B(α)=∫ωLωHΣk,m=1k≠m|ΣnRm,k(tn,ω(1+nα))|2dω]]>的值，B(α)最大值对应的α值就是所需的目标运动影响因子，其中Rm，k(tn，ω(1+nα))为第n个采样时刻tn对应的第二互谱矩阵中的第m行第n列的元素，当n为零时对应于第一采样时刻；ωH、ωL为起止频率。
全文摘要
本发明公开了一种用于海中垂直阵接收信号互谱矩阵估计中的噪声抑制方法。垂直阵在第一采样时刻和与其相邻的多个采样时刻接收到目标声源发出的信号。用第一采样时刻的接收信号构造出第一互谱矩阵，用与第一采样时刻相邻的多个采样时刻的接收信号相应构造出多个第二互谱矩阵。将多个第二互谱矩阵相对于第一互谱矩阵作频移，使得每一个第二互谱矩阵中的有效信号信道分量都与第一互谱矩阵中的有效信号信道分量相等。将第一互谱矩阵与多个频移后第二互谱矩阵作统计平均得到第一采样时刻的第一互谱矩阵的估计。本发明的噪声抑制方法采用频移后的互谱矩阵作多样本统计，可利用更多的样本进行统计，互谱矩阵估计的信噪比比常规统计方法提高十分贝以上。
文档编号H04R1/44GK1684142SQ20041003115
公开日2005年10月19日 申请日期2004年4月13日 优先权日2004年4月13日
发明者宫在晓, 张仁和, 黄晓砥 申请人:中国科学院声学研究所