水声目标弱线谱软判决被动探测方法

1.本发明涉及一种基于检测前跟踪的水声目标弱线谱软判决被动探测方法，属于水声目标探测技术领域。


背景技术：

2.舰船等辐射噪声主要包括连续的宽带谱及离散的窄带线谱。低频的窄带线谱具有较高的稳定性，因此对于目标的被动探测具有重要意义。被动声纳通过接收目标的辐射噪声从而完成信号处理、目标检测跟踪识别等一系列工作，然而除了目标的辐射噪声外，被动声纳会接收到全部的水下环境噪声及非目标干扰的辐射噪声，因此当目标距离我们较远或目标本身的辐射噪声能量较低时，其线谱会淹没在其他噪声中，呈现出能量起伏或断断续续的状态。传统的线谱自动检测方法在对线谱进行检测时，通常采用一个固定或自适应门限，对整个谱图中的频点进行筛选。然而在低信噪比情况下，检测门限通常难以选取，门限过高容易导致较弱的频点漏检，门限过低又带来大量虚警。无论是漏检还是虚警，都为后续的目标识别，状态估计带来困难。


技术实现要素：

3.针对现有检测前跟踪门限检测出现的弱目标线谱检测概率低，虚警概率高的问题，本发明提供一种水声目标弱线谱软判决被动探测方法。
4.本发明的一种水声目标弱线谱软判决被动探测方法，所述方法包括：
5.s1、将被动声纳系统的接收信号进行短时傅里叶变换；
6.s2、结合线谱的状态建立状态方程及量测方程；
7.状态方程为：x
k
＝x
k
‑1+q
k
；
8.量测方程为：z
k
＝h(e
k
x
k
,w
k
)；
9.其中，x
k
表示k时刻z
k
的状态量x
k
＝[f
k
,a
k
]
t
，q
k
＝[q
f
,q
a
]
t
，f
k
和a
k
分别表示线谱k时刻的频率和幅值，q
f
和q
a
分别表示频率和幅值的状态噪声，z
k
表示k时刻的量测，函数h(.)表示s1的短时傅里叶变换，w
k
表示k时刻的噪声，e
k
表示标记状态，e
k
＝1表示跟踪的水声目标线谱存在，e
k
＝0表示跟踪的水声目标线谱不存在；
[0010]
s3、将s1得到的短时傅里叶变换结果结合s2中的状态方程及量测方程，基于粒子滤波，得到线谱的软决策检测结果，所述软决策检测结果为：
[0011][0012]
其中，p
e,k
表示在k时刻的线谱存在概率，m表示粒子滤波产生粒子的数量，ei
k
表示第i个粒子的e
k
。
[0013]
作为优选，所述s1包括：
[0014]
接收信号的离散时域信号y(n)为：
[0015]
y(n)＝s(n)+w(n)
[0016]
s(n)表示线谱窄带信号，w(n)表示宽带噪声；
[0017]
信号y(n)的n点stft变换，获得y
k
(f)：
[0018]
y
k
(f)＝s
k
(f)+w
k
(f)
[0019]
其中，f表示离散的归一化频率；f0表示信号所在的频率，n表示短时傅里叶变换的点数。
[0020]
作为优选，量测z
k
的预测值为：
[0021][0022]
作为优选，所述s3包括：
[0023]
当k＝0时，执行s31，当1≤k≤k时，执行s32至s37；
[0024]
s31、按照均匀分布产生m个粒子，得到第i个粒子的初始状态权值同时按照目标初始存在概率p
birth
为每一个粒子赋予标记状态
[0025]
p
birth
表示目标从不存在转移至存在的概率，p
death
表示目标从存在转移至不存在的概率；
[0026]
s32、粒子标记状态转移服从两步马尔科夫过程，即有：
[0027][0028][0029][0030][0031]
s33、对于所有标记状态为且的粒子，有状态量
[0032]
对于所有标记状态为且的粒子，按照均匀分布产生粒子，得到状态量
[0033]
s34、进行量测预测：
[0034][0035]
s35、对每个粒子的权值进行计算：
[0036][0037]
其中，
i0()是0阶修正贝塞尔函数，表示方差；
[0038]
s36、通过每个粒子的权值对粒子进行重采样：将权值大的粒子进行复制，权值小的粒子删除，粒子总数不变，每个粒子的权值变为
[0039]
s37、获取软决策检测结果：
[0040][0041]
中上标i表示粒子的序号，下标k表示时刻。
[0042]
本发明的有益效果：本发明将检测与跟踪看作一个整体，不预先对线谱进行检测，而是将所有的数据输入都看作量测值，沿着线谱的运动轨迹，在完成其状态估计的同时完成检测，同时不直接给出检测结果，而是给出轨迹上每一时刻对应线谱的软判决置信级，本发明与一般检测器的非0即1的判决结果相比，软判决探测的结果更加灵活直观。在目标检测过程中，尤其是对弱目标的检测过程中，可以给出一个目标的存在概率而非直接判决，可以有效地减少漏检，并通过系统的要求完成后续的处理。且利用本发明探测，检测概率明显提高且检测结果稳定。
附图说明
[0043]
图1为本发明的流程示意图；
[0044]
图2为当输出频域信噪比为3db时的lofar谱图；
[0045]
图3为本发明pf
‑
tbd的软判决存在概率历程图；
[0046]
图4为100次蒙特卡洛试验下不同信噪比(3db,6db,10db)pf
‑
tbd的软判决检测结果；
[0047]
图5为3db和6db下的cfar和pf
‑
tbdroc曲线对比图，图5(a)表示3db下，图5(b)表示6db下；
[0048]
图6为虚警概率为0.5％、1％、5％下的cfar和tbd的roc曲线；
[0049]
图7为目标运动轨迹；
[0050]
图8为目标和被动声纳之间的距离；
[0051]
图9为接收信噪比；
[0052]
图10为stft处理后的lofar谱图；
[0053]
图11为本发明pf
‑
tbd处理得到的软判决检测结果；
[0054]
图12为当虚警概率为1％时的cfar和pf
‑
tbd的检测结果，图12(a)表示cfar方法，图12(b)表示本发明的pf
‑
tbd。
具体实施方式
[0055]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0056]
需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。
[0057]
下面结合附图和具体实施例对本发明作进一步说明，但不作为本发明的限定。
[0058]
如图1所示，本实施方式的水声目标弱线谱软判决被动探测方法，包括：
[0059]
步骤一、将被动声纳系统的接收信号进行短时傅里叶变换stft；
[0060]
步骤二、结合线谱的状态建立状态方程及量测方程；
[0061]
状态方程为：x
k
＝x
k
‑1+q
k
；
[0062]
量测方程为：z
k
＝h(e
k
x
k
,w
k
)；
[0063]
其中，x
k
表示k时刻z
k
的状态量x
k
＝[f
k
,a
k
]
t
，q
k
＝[q
f
,q
a
]
t
，f
k
和a
k
分别表示线谱k时刻的频率和幅值，q
f
和q
a
分别表示频率和幅值的状态噪声，z
k
表示k时刻的量测，函数h(.)表示s1的短时傅里叶变换，w
k
表示k时刻的噪声，e
k
表示标记状态，e
k
＝1表示跟踪的水声目标线谱存在，e
k
＝0表示跟踪的水声目标线谱不存在；
[0064]
步骤三、将步骤一得到的短时傅里叶变换结果结合步骤二中的状态方程及量测方程，基于粒子滤波，得到线谱的软决策检测结果，所述软决策检测结果为：
[0065][0066]
其中，p
e,k
表示在k时刻的线谱存在概率，m表示粒子滤波产生粒子的数量，表示第i个粒子的e
k
。
[0067]
本实施方式在检测前将检测与跟踪看作一个整体，不预先对线谱进行检测，而是将所有的数据输入都看作量测值，沿着线谱的运动轨迹，在完成其状态估计的同时完成检测，利用检测前跟踪的方法对水下弱目标线谱进行跟踪，在每一步的检测完成后，不直接给出检测硬判决“0”“1”结果，而是给出一个0
‑
1之间的软判决检验统计量，为每一时刻的检测频点赋予一个置信级。
[0068]
本实施方式的目的是对频谱中目标线谱进行检测与跟踪，因此获得的量测值为信号经过stft得到的结果。
[0069]
本实施方式的步骤一、对接收信号进行stft变换，包括：
[0070]
针对于线谱检测问题，离散时域信号y(n)可以表示为线谱窄带信号与宽带噪声的叠加：
[0071]
y(n)＝s(n)+w(n)
ꢀꢀꢀ
公式1
[0072]
其中：w(n)为宽带噪声，服从均值为0，方差为的高斯分布；s(n)可表示为一个正弦信号
[0073]
s(n)＝a cos(2πf0n)
ꢀꢀꢀ
公式2
[0074]
f0表示信号所在的频率，a表示信号的幅值。信号y(n)的n点stft变换为：
[0075][0076]
其中f为离散的归一化频率，k表示时刻(k＝1,2,...,k)，r表示矩形窗。假设信号y(n)的长度为l且l＝kn，矩形窗的长度为n：
[0077][0078]
那么公式3可以写为：
[0079][0080]
由于stft的结果可以表示为信号stft与噪声stft变换的累加，因此y
k
(f)还可写为：
[0081]
y
k
(f)＝s
k
(f)+w
k
(f)
ꢀꢀꢀ
公式6
[0082]
公式6中，w
k
(f)可以写为：
[0083][0084]
公式7中w
r,k
(f)和w
i,k
(f)分别是w
k
(f)的实部与虚部。由于w
k
(f)可以看做是w(n)的线性组合，因此w
k
(f)的虚部与实部也服从高斯分布，若式e(
·
)表示均值d(
·
)表示方差，则有：
[0085]
e(w
r,k
(f))＝e(w
i,k
(f))＝0
ꢀꢀꢀ
公式8
[0086][0087]
同时w
r,k
(f)、w
i,k
(f)是相互独立的，因为有：
[0088][0089]
从而可得其中
[0090]
公式6中，s
k
(f)是s(n)的stft结果：
[0091][0092]
若只考虑f＞0的情况，则有：
[0093][0094]
那么通过w
k
(f)、s
k
(f)的变换结果可知，y
k
(f)服从复高斯分布：
[0095][0096]
在检测前跟踪过程中，假设量测值z
k
为y
k
(f)的包络值，则有：
[0097]
z
k
＝|y
k
(f)|,k＝1,2,...,k
ꢀꢀꢀ
公式14
[0098]
本实施方式的步骤二结合线谱的状态建立状态方程及量测方程，具体包括：
[0099]
在线谱检测跟踪过程中，k时刻线谱的状态为x
k
＝[f
k
,a
k
]
t
，由于短时间内线谱及其功率的变化较小，可假设其为0速模型：
[0100]
x
k
＝x
k
‑1+q
k
ꢀꢀꢀ
公式15
[0101]
其中q
k
＝[q
f
,q
a
]
t
，q
f
和q
a
分别为频率和幅值的状态噪声。
[0102]
同时定义线谱存在变量e
k
，表示标记状态，e
k
＝1表示目标线谱存在，e
k
＝0表示目标线谱不存在。存在变量的转移服从两步马尔科夫转移，转移参数为：
[0103]
p
birth
＝p(e
k
＝1|e
k
‑1＝0)
ꢀꢀꢀ
公式16
[0104]
p
death
＝p(e
k
＝0|e
k
‑1＝1)
ꢀꢀꢀ
公式17
[0105]
检测前跟踪过程中，获得到的量测值为未经门限处理的数据，而非先通过检测得到的点迹。因此，在k时刻我们获得到的量测数据为一段经过处理的数据。假设k时刻我们获得了长度为n的信号，那么k时刻的量测方程就可以用向量z
k
表示：
[0106]
z
k
＝h(e
k
x
k
,w
k
)
ꢀꢀꢀ
公式18
[0107]
针对于线谱检测与跟踪来说，函数h(.)形式取决于步骤一中信号stft的过程。w
k
表示k时刻的噪声。k时刻目标线谱不存在时，即e
k
＝0时，量测值为h(0,w
k
)；k时刻目标线谱存在时e
k
＝1，量测值为目标的幅值分布为h(x
k
,w
k
)。
[0108]
根据步骤一中公式12的stft结果,量测z
k
的预测值可写为:
[0109][0110]
本实施方式的步骤三是对弱目标检测前跟踪方法处理：
[0111]
将步骤一获得的量测值，结合步骤二中建立的状态方程、量测方程，通过具体的检测前跟踪方法进行处理，最终获得弱目标的检测、跟踪结果。
[0112]
本实施方式的步骤三是基于粒子滤波，线谱检测前跟踪软决策过程实现具体过程为：
[0113]
设共有k个时刻，当k＝0时进行步骤三一的初始化过程；当1≤k≤k时进行步骤三二至步骤三七的循环：
[0114]
步骤三一、初始化：在一定的位置、速度范围内，按照均匀分布产生m个粒子，得到第i个粒子的初始状态权值同时按照目标初始存在概率p
birth
为每一个粒子赋予标记状态
[0115]
步骤三二、粒子标记状态转移：
[0116]
粒子标记状态转移服从两步马尔科夫过程，利用p
birth
表示目标线谱从不存在转移至存在的概率，p
death
表示目标线谱从存在转移至不存在的概率，即有：
[0117][0118][0119]
[0120][0121]
步骤三三、状态预测：
[0122]
对于所有标记状态为且的粒子，有
[0123][0124]
对于所有标记状态为且的粒子，按照均匀分布产生粒子，得到状态量
[0125]
步骤三四、量测预测
[0126][0127]
步骤三五、权值计算：
[0128]
利用似然函数对每个粒子的权值进行计算，通过步骤二中获得的量测z
k
的性质，可知，当e
k
＝0时，z
k
服从瑞利分布，有
[0129][0130]
当e
k
＝1时，z
k
服从莱斯分布，有
[0131][0132]
其中i0()是0阶修正贝塞尔函数：
[0133][0134]
因此有：
[0135][0136][0137]
步骤三六、重采样：
[0138]
通过每个粒子的权值对粒子进行重采样，将权值大的粒子进行复制，权值小的粒子删除，粒子总数不变，即将粒子的权值由粒子数量的多少表示。每个粒子的权值变为粒子的后验分布由粒子的数量来表示。
[0139]
步骤三七、获取软决策检测结果：
[0140]
计算重采样后标记的粒子的数量。则k时刻的目标线谱软决策检测结果定义为
[0141][0142]
p
e,k
为k时刻的线谱存在概率，也为最终输出的软决策结果。
[0143]
仿真实验：
[0144]
a)基于本实施方式的探测方法进行仿真，将本实施方式的探测方法简称为pf
‑
tbd：
[0145]
现仿真一突然出现，突然消失的弱目标，通过观察其检验统计量的特性，来观察pf
‑
tbd实现软决策检测的可行性。
[0146]
仿真条件：假设观察总时长为400s，目标在100s
‑
300s期间发射频率为100hz的信号，时域信号噪声为高斯白噪声，均值为0，方差、stft时长1s，点数n＝15000，采样频率f
s
＝15khz，pf
‑
tbd参数：参与跟踪的粒子数m＝500，p
bi
＝p
death
＝0.05，跟踪间隔t＝1s。
[0147]
图2为当输出频域信噪比为3db时的lofar谱图。图3为pf
‑
tbd进行在频谱上进行检测跟踪后的存在概率历程图，图中，色棒的范围为0
‑
1，表示每个时刻频点上的目标存在概率，也就是我们得到的软决策结果。
[0148]
对3db、6db、10db的目标进行100次蒙特卡洛仿真之后将检验统计量进行平均，可以看出更为清晰的结果对比，如图4所示。图中可以明显的看出，通过pf
‑
tbd的检验统计量可以作为目标检测的一个软判决结果，在0
‑
100s，300
‑
400s目标不存在时，目标的存在概率在0.18左右，在100
‑
300s，目标出现时，存在概率明显上升，且目标的信噪比越高，存在概率越高。与一般检测器的非0即1的判决结果相比，这种软判决的结果更加灵活直观。在目标检测过程中，尤其是对弱目标的检测过程中，可以给出一个目标的存在概率而非直接判决，可以有效地减少漏检，并通过系统的要求完成后续的处理。
[0149]
b)检测性能对比
[0150]
通过画出roc曲线，来对比基于pf
‑
tbd方法与常用的基于奈曼皮尔逊准则的恒虚警(cfar)线谱检测器的检测性能。
[0151]
pf
‑
tbd的探测输出结果为公式31的存在概率软决策结果，因此在衡量roc曲线时可将存在概率当做检验统计量，设置不同的可调节门限，通过蒙特卡洛仿真统计出检测概率与虚警概率之间的关系。
[0152]
仿真条件：分别仿真频域输出信噪比为3db、6db时cfar检测器和pf
‑
tbd检测器的roc曲线。其中，时域信号噪声为高斯白噪声，均值为0，方差stft时长1s，点数n＝15000，采样频率f
s
＝15khz，目标的频率f0＝100hz。pf
‑
tbd参数：参与跟踪的粒子数m＝500，p
birth
＝p
death
＝0.05，跟踪间隔t＝1。共进行10000次蒙特卡洛仿真，得到的roc曲线如图5所示；
[0153]
从图5中可以看出，在相同虚警概率的条件下，基于pf
‑
tbd检测器的检测概率明显要高于cfar检测器。
[0154]
图6中给出了当虚警概率为0.5％、1％、5％的情况下，两种检测器检测概率随随信噪比的变化的roc曲线。从图6中曲线可知，与cfar检测器相比，pf
‑
tbd检测器可带来3
‑
6db的增益。
[0155]
c)运动目标的线谱被动检测仿真
[0156]
假设有一运动目标从远处先靠近再远离被动声纳，目标做匀速直线运动，速度为(
‑
3m/s,
‑
4m/s).目标的初始位置为(7.5km,5km)，被动声纳的位置为(0km,0km)，目标一共运动3500s.假设目标线谱为100hz,噪声为均值为0方差为1的高斯白噪声。通过仿真得到的目标运动轨迹如图7所示，目标和被动声纳之间的距离变化如图8所示。
[0157]
被动声纳方程可以写为：
[0158]
s
r(t)
＝sl
‑
tl
‑
(nl+10lgb)
ꢀꢀꢀ
公式32
[0159]
其中s
r
为接收信噪比，sl为声源级，tl为传播损失，nl为环境噪声。若用r表示目标和被动声纳之间的距离，则有：tl＝20lgr，公式32给出的是时域的接收信噪比s
r(t)
，通过stft处理后，可以得到频域的信噪比s
r(f)
(在每频点)：
[0160]
s
r(f)
＝s
r(t)
+10lgb+10lgt＝sl
‑
tl
‑
nl+10lgt
ꢀꢀꢀ
公式33
[0161]
b为处理带宽，t为stft处理时长。
[0162]
设stft处理时长t＝10s，采样频率f
s
＝15khz，处理带宽为全频段当sl＝135db，nl＝70db时，s
r(f)
可以通过公式32和公式33计算，接收信噪比s
r(f)
随时间的变化如图9所示；
[0163]
图10展示了stft处理后的lofar谱图，图11展示了pf
‑
tbd软判决检测后的结果，从图11中可以看出100hz处的存在概率在运动目标靠近被动声纳时逐渐升高，远离时随着信噪比减小而降低。接下来将传统的cfar检测器与pf
‑
tbd检测方法进行对比。当保证虚警概率为1％时，得到的检测结果如图12所示：
[0164]
从图12可以看出，当检测运动目标的线谱时，pf
‑
tbd的检测效果要更好。随着信噪比的升高(目标靠近被动声纳)，两种检测方法的检测效果都变得更好，但是cfar检测器的检测结果断断续续，pf
‑
tbd方法在600
‑
2100s的时候可以给出一个相对稳定的检测结果。
[0165]
虽然在本文中参照了特定的实施方式来描述本发明，但是应该理解的是，这些实施例仅仅是本发明的原理和应用的示例。因此应该理解的是，可以对示例性的实施例进行许多修改，并且可以设计出其他的布置，只要不偏离所附权利要求所限定的本发明的精神和范围。应该理解的是，可以通过不同于原始权利要求所描述的方式来结合不同的从属权利要求和本文中所述的特征。还可以理解的是，结合单独实施例所描述的特征可以使用在其他所述实施例中。