具体实施方式

下面结合附图对本发明的具体实施方式做进一步说明。

本申请公开了一种基于主被动联合声纳阵列的目标路径跟踪方法，该方法利用多个声纳构成的主被动联合声纳阵列来进行目标路径跟踪，主被动联合声纳阵列中包括若干个排列形成预定阵列形式的主动声纳和被动声纳，如图1所示。其中，每个被动声纳对应一个全向波束、能够测量360角度范围内的目标。每个主动声纳对应至少两个探测角度范围不同的有向波束，每个有向波束只能探测相应的探测角度范围内的目标，且在同一个时刻，一个主动声纳对应的所有的有向波束中最多只有一个在工作，有向波束之间的切换时间较短可以忽略不计、可以认为是在瞬时完成的。

在主被动联合声纳阵列中，主动声纳在被跟踪目标进入其探测范围后开始发射声信号输出测量信息，而被跟踪目标是否进入主动声纳的探测范围，则由被动声纳来判断，主动声纳和被动声纳共同协调工作实现目标路径跟踪。但是主动声纳在工作时需要发射声信号实现对目标的探测，这种信号是被动声纳工作时的干扰源，同时不同主动声纳之间也会彼此构成干扰源，因此在主被动联合声纳阵列工作过程中，需要对各个声纳的开关机状态进行规划使得彼此间可以协调工作，本申请的协调过程如下：

在初始时刻、主被动联合声纳阵列刚开始工作时，所有被动声纳都处于开机状态，所有主动声纳都处于关机状态，当被动声纳发现被跟踪目标时可以获取到其距离信息和方位信息，然后根据主被动联合声纳阵列的排布结构，将该被跟踪目标分配给目标主动声纳的目标波束并打开目标主动声纳切换至目标波束，目标主动声纳是与被跟踪目标距离最近的一个主动声纳，目标波束是主动声纳的探测角度范围与被跟踪目标的方位信息匹配的一个有向波束。然后根据波束干扰关系关闭与目标主动声纳的目标波束存在波束干扰的波束，与目标主动声纳的目标波束存在波束干扰的波束可能是被动声纳的全向波束也可能是主动声纳的其中一个有向波束。一旦被跟踪目标分配给相应的主动声纳的有向波束，则被动声纳不再处理与该被跟踪目标的相关信息。

主被动联合声纳阵列中各个波束之间的波束干扰关系可以根据排布结构预先确定，若两个波束同时工作时彼此会造成干扰、则这两个波束存在波束干扰，否则两个波束不存在波束干扰、可以同时工作，且根据上述介绍，可以认为同一个主动声纳对应的各个有向波束之间都是存在波束干扰的。

在主被动联合声纳阵列的工作过程中，被跟踪目标往往不只有一个，当有多个被跟踪目标时，就需要考虑在探测资源冲突的情况下优先探测哪些被跟踪目标，因此需要确定各个被跟踪目标的处理优先级。被跟踪目标的处理优先级与被跟踪目标的目标属性对应，则在通过被动声纳第一次检测到被跟踪目标时，根据被跟踪目标的频率域信息确定被跟踪目标的目标属性，目标属性常见的包括舰船和蛙人等等，不同目标属性对应的处理优先级通常预先配置。可选的，被跟踪目标的处理优先级还与被跟踪目标的运动速度对应，被跟踪目标的运动速度越快、则对应的处理优先级越高，这是因为运动速度较慢的被跟踪目标可以用外推数据替代量测数据，而运动速度较快的被跟踪目标则需要更多的量测数据才能实现准确跟踪。实际处理时，被跟踪目标的处理优先级可以通过目标属性和运动速度的任一项来确定，较为优选的，本申请的一个实施例中，将目标属性作为第一层次的优先级的划分依据，将运动速度作为第二层次的优先级的划分依据，同时通过这两个参数来对各个被跟踪目标的处理优先级进行排序。

然后基于各个被跟踪目标的处理优先级来调整各个声纳的工作状态，具体的：将所有被跟踪目标按照处理优先级依次进行重新分配，所有被跟踪目标包括通过被动声纳检测新出现的被跟踪目标以及已经出现过的被跟踪目标，新出现的被跟踪目标未被分配给任何主动声纳，已经出现过的被跟踪目标则已经被分配给对应的目标主动声纳。在对每个被跟踪目标进行重新分配时，具体做法与上述单个被跟踪目标的方法类似，本申请不再详细赘述。

每隔预定时长根据探测结果重新调整主被动联合声纳阵列中各个声纳的工作状态，而在每两次调整之间的时间段内，主被动联合声纳阵列中各个声纳的工作状态保持不变，本申请将两次调整之间的时间段定义为一个检测时段，则在每个检测时段内、各个声纳的工作状态保持不变，在当前检测时段结束时，重新调整各个声纳的工作状态进入下个检测时段。

如上所述，每个被跟踪目标在进入主被动联合声纳阵列的范围内后，会被分配给至少一个声纳的波束，也即每个被跟踪目标对应至少一个声纳，主被动联合声纳阵列中的各个声纳可以获取到对应的被跟踪目标的相关信息记为量测变量，包括如上提到的频率域信息、距离信息和方位信息等等。但是各个声纳的量测变量对应的时间是不具备统一的时间戳的，而且浅海声传播信道的复杂特性导致主动、被动声纳测量噪声很难用确定性统计参数描述，而是呈现出高非线性、强时变性特征，因此这些不同步的量测变量需要进一步融合计算才能实现路径跟踪，本申请利用基于先验概率的卡尔曼融合估计滤波器来实现融合计算。

被跟踪目标的运动方程定义为如下形式X_i＝F_iX_i-1+W_i，X_i＝[x_p,i,x_v,i,y_p,i,y_v,i]^T是描述被跟踪目标在第i个检测时刻的运动状态的向量也即状态变量，x_p,i表示被跟踪目标以坐标系x轴为参考的位置，x_v,i表示被跟踪目标以坐标系x轴为参考的速度，y_p,i表示被跟踪目标以坐标系y轴为参考的位置，y_v,i表示被跟踪目标以坐标系y轴为参考的速度，坐标系x轴通常选取所在地的地理东向，坐标系y轴选取地理北向。F_i是由运动学方程确定的状态转移矩阵，W_i为零均值高斯白噪声，其方差矩阵定义为Q_i。在第i个检测时刻，将通过主被动联合声纳阵列中对应的声纳共获取被跟踪目标的N_i个量测变量记为其中第k个量测变量表示被跟踪目标与第k个声纳之间的距离，k≥1。在一个检测时段内，可以有多个检测时刻。

量测模型写成其中H(X_i,X^k)＝||X_i-X^k||₂，X^k表示第k个声纳的坐标值。V_i ^k为零均值高斯白噪声，其方差矩阵定义为由于测量精度与距离有相关性，所以的特征值与所述第k个量测变量的大小相关，也即是距离的函数，可以记为

本申请利用基于先验概率的滤波器进行量测变量和状态变量的更新，具体方法实现如下：

定义第i个检测时刻获取到的N_i个量测变量记为从第一个检测时刻到第i个检测时刻的全部量测变量记为Z₁:_i，第i-1个检测时刻的先验概率密度函数记为p(X_i-1|Z_1:i-1)，第i-1个检测时刻到第i个检测时刻的状态转移概率密度函数记为p(X_i|X_i-1)，状态变量X_i与第k个量测变量之间的似然函数为在上述定义下，滤波过程分为三部分：

(1)初始化：设定初始时刻的概率分布函数p(X₀)。

(2)预测：p(X_i|Z_1:i-1)＝∫p(X_i|X_i-1)p(X_i-1|Z_1:i-1)dX_i-1。

(3)量测更新：

(4)量测空间更新：

在上述算法中，为了提高概率密度函数计算效率，采用高斯噪声特性进行线性变换，具体的：

(1)预测环节的状态变量预测计算公式为：

P_i|i-1＝F_iP_i-1|i-1F_i ^T+Q_i；

(2)量测更新环节的状态变量与量测变量的联合概率密度计算公式为：

其中，

(3)量测空间更新环节的状态变量与量测变量的联合概率密度计算公式为：

其中，

联合高斯概率密度服从以下分布：其中：

第k个预测值及其方差矩阵分别由如下两个公式计算：

状态变量与量测变量之间的第k个耦合方差矩阵为：

其中，H^k(X_i)＝||X_i||₂且X_i表示状态变量， 为零均值高斯白噪声的方差矩阵，为预设初始值。

基于上述算法原理得到采用高斯噪声特性对基于先验概率的滤波器进行线性变换后的卡尔曼融合估计滤波器，在第i个检测时刻的量测更新过程为：

(1)确定第i个检测时刻的状态变量的初始值及其均方根矩阵的初始值除此之外还会定义方差矩阵P_i|i。

其中，F_i是由运动学方程确定的状态转移矩阵，为预设的值、通常可以定义为大于0的实数。方差矩阵P_i|i通常可以定义为一个正定矩阵。S_i-1|i-1是P_i-1|i-1的均方根矩阵，P_i-1|i-1是第i-1个检测时刻的方差矩阵且通常为一个正定矩阵、S_Q,i是零均值高斯白噪声W_i的方差矩阵Q_i的均方根矩阵。

(2)根据第k-1个状态变量及其对应的第k-1个均方根矩阵计算得到第k个预测值计算公式为其中，N_x为系数且2N_x为[0,k]区间内最大的偶数，r为系数，

(3)确定状态变量与量测变量之间的第k个耦合方差矩阵

本申请根据第k个预测值第k-1个状态变量和第k-1个均方根矩阵来计算计算公式为：

(4)更新得到第k个状态变量及其对应的第k个均方根矩阵

在更新时，主要利用第k个量测变量第k个预测值和第k个耦合方差矩阵计算公式为：

其中，函数QR()用于对矩阵进行QR分解，是零均值高斯白噪声的方差矩阵的均方根矩阵，表示对应的增益矩阵，由测量值的高斯分布以及估计值的高斯分布的乘积确定，实际计算过程中，一般可近似地将取值为单位矩阵。

(5)若k<N_i，则令k＝k+1并再次执行根据第k-1个状态变量及其对应的第k-1个均方根矩阵计算得到第k个预测值的步骤，也即继续确定第k+1个状态变量。

(6)若k＝N_i，则利用得到的N_i个状态变量对被跟踪目标进行路径跟踪。

到第i+1个检测时刻再利用上述循环进行处理以及路径跟踪。

在上述方法中，如上提到的特征值与所述第k个量测变量的大小相关，而计算过程中真实距离很难获得，如果未能保持连续的估计，估计器容易出现发散的情况，为了计算的连续性，提高估计器的自适应特性和鲁棒性，该实施例还会利用与相关的第k个自适应伸缩因子补偿具体的：

计算 是第i个检测时段之前的所有量测变量的均值，的公式见上已公开内容。

检测是否满足其中，为预设参数。

若则基于和按照上述方法进行计算。

若则利用与相关的第k个自适应伸缩因子补偿得到并修正得到然后基于和进行计算，后续具体的方法都是类似的。

在上述过程中，自适应伸缩因子的计算公式为其中 是的均方根矩阵。

以上所述的仅是本申请的优选实施方式，本发明不限于以上实施例。可以理解，本领域技术人员在不脱离本发明的精神和构思的前提下直接导出或联想到的其他改进和变化，均应认为包含在本发明的保护范围之内。