一种长时波束相位统计特性的线谱检测方法
阅读说明：本技术 一种长时波束相位统计特性的线谱检测方法 (Line spectrum detection method for long-time beam phase statistical characteristics ) 是由 吴艳群 张兵兵 郭微 朱家华 王俊 张文 徐国军 胡正良 于 2021-07-13 设计创作，主要内容包括：本发明属于水声工程、海洋工程和声呐技术等领域,具体涉及一种长时波束相位统计特性的线谱检测方法,先通过水平阵列的波束形成获得水平各个方位各个频点输出的复波束序列,在对准目标方位的输出序列信噪比得到了空间增强；进而以第一快拍的复波束序列为参考,计算其它时刻的复波束序列与参考复波束序列的互谱,获得互谱序列。本发明的特点在于首先对水平阵列进行波束形成,在增强线谱信号的信噪比的同时一定程度上抑制了其它方向上的强干扰；同时,基于长时互谱序列的相位统计特性来区分线谱与各类宽带噪声,其检测门限易于确定且相对稳定。(The invention belongs to the fields of underwater acoustic engineering, ocean engineering, sonar technology and the like, and particularly relates to a line spectrum detection method for long-time beam phase statistical characteristics, wherein a complex beam sequence output by each frequency point in each horizontal direction is obtained through beam forming of a horizontal array, and the signal-to-noise ratio of the output sequence aligned to a target direction is spatially enhanced; and further taking the complex beam sequence of the first snapshot as a reference, calculating the cross spectrum of the complex beam sequence at other moments and the reference complex beam sequence, and obtaining a cross spectrum sequence. The method is characterized in that firstly, the horizontal array is subjected to beam forming, so that the signal-to-noise ratio of a line spectrum signal is enhanced, and strong interference in other directions is inhibited to a certain extent; meanwhile, line spectrum and various broadband noises are distinguished based on the phase statistical characteristics of the long-term cross-spectrum sequence, and the detection threshold is easy to determine and relatively stable.)
技术领域
本发明属于水声工程、海洋工程和声呐技术等领域，具体涉及一种长时波束相位统计特性的线谱检测方法。
背景技术
舰船辐射噪声通常由连续谱和线谱叠加而成，其中线谱噪声由辅机机械的往复运动引起，其频率低，可远距离传播，且谱级比连续谱高10-25dB，因此线谱是水下目标远程探测和识别的重要特征。线谱检测主要采用功率谱估计、现代谱估计等方法，二者都以离散傅里叶变换(Discrete Fourier Transform，DFT)为基础采用假设检验类线谱检测方法，依据信号的统计特性构建检测器。低频线谱分析历程图(Low Frequency Analysis AndRecording，LOFAR)利用线谱信号具有时间连续性的特性，通过对一段接收数据连续进行短时傅里叶变换(Short-Time Fourier Transform，STFT)计算得到。LOFAR图中线谱检测的主要原理是利用人眼的视觉累积效应，对一定时间内连续的线谱信号进行频率的自动跟踪。
为了提升线谱检测能力，美国海军实验室的Wagstaff等人利用线谱信号的相位特征具有长时稳定性，而噪声信号的相位长时间观测具有随机性，发明了Wagstaff’sintegration silence processor(WISPR)及一系列派生处理技术，可对相位缓慢变化的线谱信号采用相位差分的方法进行相位补偿从而实现线谱信号的长时间积分，以获得足够的处理增益(参见:Exploiting phase fluctuations to improve temporal coherence,IEEE J.Ocean.Eng.,2004,29(2):498–510.以及The AWSUMfilter:A 20-dB gainfluctuation-based processor，IEEE J.Ocean.Eng.,1997,22(1):110–118)。解恺等人结合空域阵列的输出增益，通过对波束域输出信号的起伏相位差分对齐，有效提高目标线谱分量的检测信噪比增益的分辨能力和波达方位的估计能力(参见：基于目标辐射噪声的信号起伏检测算法研究，《电子与信息学报》，2013，35(4):844-851)。
但是现有检测方法存在如下问题：一方面，海洋环境背景噪声起伏或者是其它舰船干扰将使得线谱的信噪比发生起伏，当线谱的某些部分因信噪比起伏，导致信噪比过低而淹没在干扰之中时，该部分线谱将无法进行检测。虽然使用阵列通过波束形成技术可以抑制水平面上来自其他方向的干扰信号，在一定程度上提升了期望目标方向上波束域信号的信噪比，但是强干扰存在时仍然可能使期望目标线谱淹没在噪声中；另一方面，由于海洋波导的干涉效应，目标运动时接收器接收到的线谱功率也会发生起伏，起伏范围能达到十几甚至二十分贝，这导致线谱在能量高时能够被检测到，在能量低时无法被检测到。考虑到噪声信号相位随机抖动，长时间观测相位方差大，而线谱信号相位相对稳定，长时间观测相位方差非常小，即使存在由于目标或者接收器相对运动带来的多普勒频移，其线谱信号相位方差在一定的时间内仍然能保持一个较小值。
发明内容
本发明所要解决的技术问题是，针对现有线谱检测抗干扰能力不足且检测门限难以确定的缺点，提供了一种长时波束相位统计特性的线谱检测方法。该方法将目标方向上波束域输出频域信号的长时间相位二次统计特性(即相位方差)作为检测器的检验量，在获得了阵列增益提高输出信噪比的同时，能够根据线谱和噪声的相位统计特性来区分线谱与各类宽带噪声(包括泄露到线谱方向上的宽带干扰能量)，具有检测门限易于确定且抗干扰能力强的优点。
本发明所采用的技术方案是：
一种长时波束相位统计特性的线谱检测方法，利用一条由M个水听器阵元构成的水听器阵列将接收到的声信号进行频域波束形成，并求解目标波束方位各个频率点上长时间的互谱相位方差作为门限值，提升线谱检测的性能。具体分为以下步骤：
S1，对水听器阵列M个阵元接收到的声信号矩阵沿着每一列分窗口进行傅里叶变换，获得频谱矩阵P(fi)；具体步骤如下：
S1.1将水听器阵列M个阵元接收到的声信号表示为一个NT×M维的矩阵p：
其中矩阵p的每一行p1(t) p2(t) … pM(t)表示水听器阵列的M个阵元在t时刻接收到的声信号，t＝1,2,…,NT；设时间采样率为fs，总观测时间为T，则每个阵元总采样点数为NT＝floor[Tfs]，这里floor表示向下取整操作。
S1.2任取一个不大于NT的正整数Nwin作为窗口长度，并任取一个小于Nwin的正整数Noverlap作为窗口重叠长度，将矩阵p的每一列数据划分为K＝floor[(NT-Nwin)/(Nwin-Noverlap)+1]个窗口；
S1.3选取矩阵p每一列的第k个窗口，即第(k-1)(Nwin-Noverlap)+1行到第(k-1)(Nwin-Noverlap)+Nwin行的数据进行Nwin点的离散傅里叶变换，获得频谱矩阵P(fi)，k＝1,2,...,K：
设为傅里叶变换后的频率轴向量，阵列最低工作频率fmin和最高工作频率fmax在fFFT向量中对应的编号分别为imin和imax，imin＝floor[fminNwin/fs]+1，imax＝floor[fmaxNwin/fs]+1；记录第m列第k个窗口傅里叶变换后的第i个频点fi处的值为Pm(fi,k)，其中fi＝(i-1)fs/Nwin，m＝1,2,...,M，i＝imin,imin+1,...,imax；依次遍历第m列的K个窗口，然后遍历所有列，获得一个K×M维的频谱矩阵P(fi)：
其中矩阵P(fi)第k行的数据P(fi,k)称为第k个快拍的数据：P(fi,k)＝[P1(fi,k),...,PM(fi,k)]，k称为快拍数。
S2，对频谱矩阵P(fi)的第k行数据P(fi,k)进行波束形成，获得波束形成后的复波束序列b(fi,θ,k)，然后计算水听器阵列工作频段内宽带波束功率输出结果B(θ,k)及其极大值对应的可能目标方位角；具体步骤如下：
S2.1对P(fi,k)进行波束形成，获得复波束序列
b(fi,θ,k)＝wH(fi,θ)PT(fi,k)
式中上标H表示矩阵的复共轭转置，上标T表示矩阵的转置；w(fi,θ)为远场平面波模型下水听器阵列的导向向量；θ为水平方位搜索角，角度搜索间隔为Δθ，其最小值和最大值分别为θmin和θmax。
S2.2根据S2.1获得的复波束序列b(fi,θ,k)，求解水听器阵列工作频段内宽带波束功率输出结果：
S2.3记录所有使B(θ,k)达到极大值时对应的方位角θsk，s＝1,2,...,S，S表示B(θ,k)极大值的个数，即第k个快拍下可能的目标个数。
S3，计算在第s个极大值方位角θsk处、频点为fi时，复波束序列b(fi,θ,k)在第k个快拍的序列值与第1个快拍的序列值的复互谱值C(fi,s,k)及其相位值φ(fi,s,k)：
C(fi,s,k)＝b(fi,θsk,k)b*(fi,θs1,1)
φ(fi,s,k)＝arg{C(fi,s,k)}
其中上标*表示向量的复共轭，arg{·}表示对复数取相位的操作。
S4，重复S3直至遍历所有K个快拍数，获得长度为K的长时统计复互谱相位序列φ(fi,s)，然后计算该序列的方差δφ(fi,s)。
长时统计复互谱相位序列可以表示为
φ(fi,s)＝arg{[C(fi,s,1),C(fi,s,2),...,C(fi,s,K)]T}
其方差为
其中为第i个频点的平均相位。
S5，将第imin至imax共计imax-imin+1个频点处所有相位方差δφ(fi,s)的倒数的平均值作为检测门限值，然后将不同频点处δφ(fi,s)的倒数逐个与比较：超出门限值的视为线谱，低于门限的值视为随机噪声。其中门限值可表示为
S6，在第s个极大值方位角θsk处，重复S3至S5直至遍历所有imax-imin+1个频点，得到第s个可能目标的线谱检测结果。
S7，重复S6直至遍历所有s个可能目标，即可完成线谱检测。
本发明基于以下原理：线谱信号的相位具有长时稳定性，而噪声信号相位的长时间序列是随机的，本发明先通过水平阵列的波束形成获得水平各个方位各个频点输出的复波束序列，在对准目标方位的输出序列信噪比得到了空间增强。进而以第一快拍的复波束序列为参考，计算其它时刻的复波束序列与参考复波束序列的互谱，获得互谱序列。对于线谱信号，该互谱序列的相位反映了目标波束方向上线谱信号的相位长时稳定性。对于噪声或者其它宽带信号，该互谱序列的相位长时间统计结果是随机的。本发明的特点在于首先对水平阵列进行波束形成，在增强线谱信号的信噪比的同时一定程度上抑制了其它方向上的强干扰；同时，基于长时互谱序列的相位统计特性来区分线谱与各类宽带噪声(包括泄露到线谱方向上的宽带干扰能量)，其检测门限易于确定且相对稳定。
与现有技术相比，本发明所具有的有益效果为：
本发明通过对水平阵列进行波束形成并计算互谱相位方差获得了线谱检测的门限值。与传统的方法相比，本发明的门限获取方法具有自适应性，且计算量小，稳健性好，便于工程实现；相比传统的功率谱法能够获得更高的线谱信噪比。
附图说明
图1是本发明实施例的总体流程图；
图2是阵列波束形成方位历程结果图；
图3是对准期望目标方向波束输出的LOFAR结果图；
图4是使用目标波束方向上的LOFAR输出检测结果(a)和本发明实施例的波束互谱输出的相位统计特性检测结果(b)比对图，其中五角星表示声源发射的线谱信号。
具体实施方式
下面结合附图和某次海试数据处理的具体实施方式，对本发明进行详细阐述。本发明提出的一种长时波束相位统计特性的线谱检测方法，通过对波束输出互谱相位进行统计分析，当检测到长时相位异常判别为线谱。附图1给出本发明实施例的总体流程图，具体将结合某次海上试验得到的数据进行详细说明。
在此次海上试验中，采用水听器个数M＝28的不等间距阵列，其阵列形状接近线列阵。以阵列第1个阵元为原点，建立xoy直角坐标系，其第m个阵元的坐标为xm＝[xm,ym]T。水平方位搜索角θ定义为与x正轴的夹角，角度搜索范围从0度到359度，角度搜索间隔Δθ为1度。
S1，对水听器阵列M个阵元接收到的声信号矩阵沿着每一列分窗口进行傅里叶变换，获得频谱矩阵P(fi)；具体步骤如下：
S1.1将水听器阵列M个阵元接收到的声信号表示为一个NT×M维的矩阵p：
其中矩阵p的每一行p1(t) p2(t) … pM(t)表示水听器阵列的M个阵元在t时刻接收到的声信号，t＝1,2,…,NT；设时间采样率为fs，总观测时间为T，则每个阵元总采样点数为NT＝floor[Tfs]，这里floor表示向下取整操作。
S1.2任取一个不大于NT的正整数Nwin作为窗口长度，并任取一个小于Nwin的正整数Noverlap作为窗口重叠长度，将矩阵p的每一列数据划分为K＝floor[(NT-Nwin)/(Nwin-Noverlap)+1]个窗口；
S1.3选取矩阵p每一列的第k(k＝1,2,...,K)个窗口，即第(k-1)(Nwin-Noverlap)+1行到
第(k-1)(Nwin-Noverlap)+Nwin行的数据进行Nwin点的离散傅里叶变换，获得频谱矩阵P(fi)：
设为傅里叶变换后的频率轴向量，阵列最低工作频率fmin和最高工作频率fmax在fFFT向量中对应的编号分别为imin和imax，那么imin＝floor[fminNwin/fs]+1，imax＝floor[fmaxNwin/fs]+1；记录第m列(m＝1,2,...,M)第k个窗口傅里叶变换后的第i个频点fi处的值为Pm(fi,k)，其中fi＝(i-1)fs/Nwin(i＝imin,imin+1,...,imax)；依次遍历第m列的K个窗口，然后遍历所有列，可获得一个K×M维的频谱矩阵P(fi)，即
其中矩阵P(fi)第k行的数据称为第k个快拍的数据，k称为快拍数。
本实施例中，时间采样率fs＝3276.8Hz，选取总观测时间为T＝5min，每个水听器总采样点数为NT＝983040。进行傅里叶变换时采用的时间窗口点数为Nwin＝65536个点，即每个时间窗口的观测时长为20s。每次选取数据时与上一次数据的重叠点数Noverlap＝45875，因此总采样时间时间T内，一共可以划分为K＝47个窗口。因此，矩阵P(fi)的维数为47×28。由于时间窗口点数为65536个点，故傅里叶变换后频点数为65536个。阵列工作频段为40Hz到410Hz，在频率轴向量fFFT上对应的频点编号分别为imin＝802和imax＝8202，因此总的频率点数为imax-imin+1＝7401个。
S2，对频谱矩阵P(fi)的第k行数据进行波束形成，获得波束形成后的复波束序列b(fi,θ,k)和其工作频段内的宽带波束功率输出结果B(θ,k)及其极大值对应的可能目标方位角；具体步骤如下：
S2.1将P(fi,k)＝[P1(fi,k),...,PM(fi,k)]记为频谱矩阵P(fi)的第k行数据，即水听器阵列M个阵元第k个快拍的声信号经过傅里叶变换后得到第i个频点的频谱信号。对P(fi,k)进行波束形成，获得复波束序列
b(fi,θ,k)＝wH(fi,θ)PT(fi,k)
式中上标H表示矩阵的复共轭转置，T表示矩阵的转置；θ为水平方位搜索角，角度搜索间隔为Δθ，其最小值和最大值分别为θmin和θmax，因此θ＝[θmin,θmin+Δθ,...,θmax]T；不失一般性地，设θmin＝-180°，θmax＝180°，在实际操作中，可根据阵列对目标的有效探测方位角范围这一先验信息缩小水平方位搜索角的范围，例如阵列由于存在左右模糊只能进行[-90°,90°]入射角范围的波束形成，可以将θmin和θmax分别设置为-90°和90°；w(fi,θ)为远场模型下平面波对应的水听器阵列导向向量。通常，在远场模型下平面波波束形成对应的水听器阵列的导向向量w(fi,θ)可表示为
为方便分析，后续内容将以上述导向向量为例。类似地，w(fi,θ)也可以为非等间距阵列或者其它几何构型的阵列在远场模型下平面波波束形成对应的水听器阵列的导向向量，不同类型的水听器阵列导向向量的表达式参见文献：R.O.Nielsen.Sonar SignalProcessing.Artech House,London,1991。
S2.2根据S2.1获得的复波束序列，进一步求解其阵列工作频段内宽带功率输出结果
S2.3记录所有使B(θ,k)达到极大值对应的方位角θsk(s＝1,2,...,S)，这里S表示B(θ,k)极大值的个数，即第k个快拍下可能的目标个数。优选地，本领域技术人员可根据期望目标方位的先验信息(例如通过船舶自动识别系统(AIS)获得期望目标方位的先验信息)减少可能的目标个数(如假设该步骤中记录到的极大值有3个，分别对应方位角θ1k、θ2k、θ3k分别为79度、80度和100度，而根据实验过程的先验信息表明可能的目标只会存在于80度附近，当方位角排除门限设为5度时，可排除100度方位角对应的极大值，从而将可能的目标个数减少为2个)。
本实施例中，水平方位角空间扫描角度为1到180度，角度间隔为1度。因此，获得了不同频率点fi和不同快拍k情况下，180个角度方向上的复波束序列输出b(fi,θ,k)。按照波束序列的功率输出B(θ,k)，得到了附图2所示的目标方位历程图。附图2给出了5分钟内B(θ,k)的输出结果。根据试验条件中的先验信息，期望目标位于80度附近。从附图2可以看到，在5分钟内，强干扰方位在8度左右，期望目标角度位于78度基本不变，因此θ11＝...＝θ1K＝78度。根据附图1，干扰能量高于目标。如附图3所示，期望目标方向上输出波束信号特征，即LOFAR分析结果|b(fi,θ1k,k)|2。比对附图2和附图3可知，虽然干扰方向和目标方向差异很大，但是强干扰的能量仍然通过旁瓣泄露到了目标方向，导致此目标方位上LOFAR图的低频部分存在很强的舰船干涉能量，而且很多线谱能量不稳定，无法从LOFAR图中观察到。
S3，计算在第s个极大值方位角θsk处、频点为fi时，复波束序列在第k个快拍的序列值与第1个快拍的序列值的复互谱值C(fi,s,k)及其相位值φ(fi,s,k)：
C(fi,s,k)＝b(fi,θsk,k)b*(fi,θs1,1)
φ(fi,s,k)＝arg{C(fi,s,k)}
其中上标*表示向量的复共轭，arg{}表示对复数取相位的操作。
本实施例中，由于数据点数被划分为K＝47个窗口，总的频率点数为imax-imin+1＝7401个，因此共获得47个复互谱相位值φ(fi,k)，即每47个快拍进行后续的线谱检测。因此对阵列工作频带范围内的7401个频点，获取的互谱相位φ(fi,k)均为47个。
S4，重复S3直至遍历所有K个快拍数，获得长度为K的长时统计复互谱相位序列，然后计算该序列的方差δφ(fi,s,k)。
长时统计复互谱相位序列可以表示为
φ(fi,s,k)＝arg{[C(fi,s,1),C(fi,s,2),...,C(fi,s,K)]T}
其方差为
其中为第i个频点的平均相位。
S5，将第imin至imax共计imax-imin+1个频点处所有相位方差δφ(fi,s)的倒数的平均值作为检测门限值，然后将不同频点处δφ(fi,s)的倒数逐个与比较：超出门限值的视为线谱，低于门限的值视为随机噪声。其中门限值可表示为
S6，在第s个极大值方位角θsk处，重复S3至S5直至遍历所有imax-imin+1个频点，得到第s个可能目标的线谱检测结果。
S7，重复S6直至遍历所有s个可能目标，即可完成线谱检测。
本实施例中，对阵列工作频带范围内的每个频点，获得期望目标的互谱相位φ(fi,1,k)的均值和方差δφ(fi,1)。附图4(b)给出了所有频点的相位方差δφ(fi,1)。为了验证本发明的有效性，图4给出了传统的功率谱法(附图4(a))和本发明的结果(附图4(b))比对。传统的功率谱法通过将附图3中所示的期望目标方向上输出波束信号特征进行多个快拍的平均，即在本实例中存在宽带强干扰目标且信干比低，此时传统的功率谱法难以通过信号能量检测到线谱。即使对附图4(a)通过去趋势平滑功率谱的趋势项，也无法提高线谱的信噪比，依然难以检测到信号。本发明利用线谱信号的相位稳定性的统计规律与宽带噪声随机相位的统计规律，将线谱信号与噪声分离，而不是采用信号能量特征，具有更好的检测效果。同时，不同噪声频点的相位方差非常接近，即噪声相位方差谱非常平坦，这为门限确定提供了便利，便于系统实现。
最后应说明的是：以上所述仅为本发明的较佳实施例而已，并非用于用于限定本发明的保护范围，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
- 上一篇：一种医用注射器针头装配设备
- 下一篇：一种用于物理教学实验的声速测量装置