潜标阵目标运动分析算法设计与仿真

冷兆龙，刘高峰，王子齐

(1.海军工程大学， 武汉 430033；2.中国人民解放军91213部队，山东 青岛 266000)

1 引言

潜标也被称为水下浮标。通过浮体和锚系设备，潜标可以停留在预设的作业位置和作业深度，持续工作数月甚至更长的时间[1-2]。以往，潜标主要应用于海洋环境监测等民用领域。为创新反潜作战样式，文献[3]提出了一种将潜标阵应用于反潜作战的构想，即在重要的海峡、航道或其他海域，前置部署由若干被动声探测潜标组成的潜标阵，用于反潜目标的侦察和指示。在潜标阵中，各潜标使用被动声探测方式，仅能获得目标方位，再通过方位对目标运动要素进行估计，这类问题可归属于纯方位目标运动分析(BOTMA)问题。文献[4]提出了一种战场威胁约束下的观测站运动轨迹优化算法；文献[5-7]分别从计算复杂度、计算精度、消除估计有偏性的角度对BOTMA算法进行了改进；文献[8]针对未知分布的噪音，构建了一种改进的自适应滤波器；文献[9]改进了伯努利滤波器，用于多目标情况下的纯方位目标跟踪；文献[10-12]针对滤波器初值选择问题进行了分析。已有的BOTMA研究，从研究内涵来看，大多聚焦于观测站运动轨迹、滤波器改良、滤波器初值选择等算法优化方面；从应用场景来看，主要有陆基监听站、舰艇、潜艇、反潜机等装备搭载平台。

由于潜标阵应用的特殊性，已有的BOTMA算法大多无法直接运用。本文建立了一种“单潜标数据独立处理(single submarine buoy data independent processing，SSBDIP)和多潜标数据集中融合(multi submarine buoy data centralized fusion，MSBDCF)”的潜标阵目标运动分析描述模型，采用最小二乘法，开展非机动/机动目标运动要素解算研究。

2 潜标阵目标运动分析问题描述

2.1 传统多站BOTMA问题算法特性

根据观测站数量，BOTMA问题可分为单站BOTMA和多站BOTMA。相关研究已表明，单站BOTMA问题在观测站不进行有效机动情况下，目标运动状态不具备完全可观测性。多站BOTMA问题中，观测数据维度增加了，因此对观测站是否机动没有要求。传统多站BOTMA问题的解决思路如图1所示。

图1 传统多站BOTMA问题的解决思路框图

先获取目标位置的粗略估计，然后利用粗略估计的位置优化目标运动模型参数，同时利用运动模型修正位置估计，通过循环反馈的方式，得到较高精度的估计结果。

传统多站BOTMA问题中，对位置的粗略估计通常采用方位线交汇定位的方式，如图2所示。图2(a)中，观测站数量为2时，取观测所获得方位线的交点作为位置的粗略估计；图2(b)中，观测站数量大于2时，由于观测误差影响，方位线往往不能交于一点，而是形成了一个多边形，这时可以按一定准则从多边形中选择某个点(如选择多边形的重心)作为位置的粗略估计。显然，方位线交汇定位需要不同观测站对目标同一时刻的观测数据，或者说需要各观测站的实时观测数据，这也是绝大多数多站BOTMA算法对观测站的要求。

图2 传统多站BOTMA问题中位置的粗略估计示意图

相对于传统以平台为中心的探潜方式，潜标阵在反潜目标探测应用中具有一定的特殊性，主要表现在以下方面：

1) 潜标阵远离本土部署在重要海域，与指挥所之间难以建立有线通信。为将潜标阵探测的目标信息回传到指挥所，潜标需要上浮，通过卫星链路与指挥所进行数据通信。

2) 各潜标观测的目标角度数据多，通过卫星链路回传时，传输数据量需进行一定的限制，否则会占用大量通信资源，进而影响通信的可靠性和稳定性。

3) 指挥所得到的目标信息数据是潜标回传之前的数据，存在一定时延。

4) 为保持潜标的隐蔽性，各潜标相互独立且彼此之间不通信。

2.2 潜标阵目标运动分析描述模型

潜标阵目标运动分析问题属于多站BOTMA问题,但不具备直接运用传统多站BOTMA问题相关算法的条件。或者说，基于传统多站BOTMA问题算法，目标对于潜标阵“部分可观测”而非 “完全可观测”。已有研究表明：对于静止单站被动观测，可以解算出目标运动的部分要素，该结论可以推广应用于潜标阵目标运动分析之中。基本思路是：首先，各潜标独立解算，得到目标运动的部分要素；其次，通过多潜标数据集中融合，解算出目标运动的全部要素。提出的潜标阵目标运动分析描述模型如图3所示。

图3 潜标阵目标运动分析描述模型框图

当目标经过潜标阵部署海域时，各潜标被动感知目标方位，并对方位数据自行处理，得到目标运动的部分要素。当潜标判断目标已离开该潜标探测范围后，自行上浮，将处理结果和潜标位置、潜标观测目标的起止时刻、运动分段起始时刻的角度观测值等数据，通过卫星链路回传至指挥所。指挥所的情报中心获得多个(2个以上)潜标数据，进行时空对准，然后解算出目标运动的全部要素。由于潜标在目标离开探测范围后才上浮回传数据，上浮过程可能需要耗时数分钟左右，因此解算结果是非实时的，存在时延。针对这一问题，可以通过预估目标实时位置(利用目标航向、航速)的方式来解决。尽管这样处理存在误差，考虑到潜艇的运动速度和机动能力都是不很强，预估的目标运动要素可认为基本足以支撑反潜作战指挥决策，也可以为反潜武器运用提供一定参考。

3 单潜标数据独立处理算法设计

单潜标数据独立处理(SSBDIP)算法目的在于利用最小二乘法，估计得到目标运动的部分要素。

3.1 非机动目标的数据处理算法

假设目标做匀速直线运动，如图4所示。

图4 单个潜标对匀速直线运动目标的观测示意图

任选一点作为坐标系原点，以正东为轴正方向，正北为x轴正方向，建立平面直角坐标系。目标进行匀速直线运动，航向为Km，速度在x轴和y轴方向上的分量为vx、vy，在任意时刻t目标的位置坐标为[rx(t),ry(t)]。潜标位置坐标为(x,y)，可以对周围半径为Dr的圆形区域范围进行被动探测。记潜标首次探测到目标的时刻为t0，此时目标相对潜标的方位为β(t0)，以y轴正方向为0°，顺时针方向为正。潜标持续探测目标，直到目标离开潜标探测区域，把这些目标方位数据分别记为：β(t1)，β(t2)，…，β(tn)，其中，t1，t2，…,tn为目标方位数据所对应的时刻。

根据三角函数关系，对于任意时刻ti，均有：

(1)

假设目标做匀速直线运动，有：

(2)

将式(2)代入式(1)，整理后可得：

cosβ(ti){[ry(t0)-y]tanβ(t0)-x+vx(ti-t0)}=

sin(ti)[ry(t0)-y+vy(ti-t0)],i∈[0,1,…,n]

(3)

已有研究表明：静止单站被动观测条件下，如果观测站(潜标)不在目标运动轨迹及其延长线上，目标就是部分可观测的。对于式(3)，只要获得了至少3个不同时刻的观测数据，便可利用最小二乘法对目标运动的2个要素[ry(t0)-y]/vy和vx/vy进行估计，估计结果为：

(4)

其中，

(5)

通过上述处理过程，得到了目标运动的2个要素：[ry(t0)-y]/vy和vx/vy。这样，潜标在利用卫星链路通信时，不需要将目标的全部观测方位信息回传，大大降低了卫星链路通信的数据量。同时，[ry(t0)-y]/vy和vx/vy均为非时敏参数，因此对实时性没有严格要求。

3.2 机动目标的识别处理

描述机动目标的运动模型有很多，在潜标阵机动目标识别过程中，选择了一种分段匀速直线运动模型，将目标运动描述为若干段相连的匀速直线运动。该模型适用范围广，可以描述转向、变速等不同类型的机动目标，而且模型参数较为简单，便于计算。

使用分段匀速直线运动模型对机动目标观测数据处理的过程中，关键问题在于识别和判断各段匀速直线运动的起止。潜标仅能获取目标方位角度信息，无法直接识别目标的运动状态是否发生变化。文献[13-14]中提出了纯方位观测条件下通过角度预测判断目标是否发生机动的方法。本文在这一方法基础上进行了改进，识别的基本过程为：

过程1) 假定目标处于匀速直线运动状态，利用式(4) 可获得[ry(t0)-y]/vy和vx/vy的估计结果。

(6)

过程3) 计算目标角度预测数据与目标角度观测数据之间的差Δβ(ta+Δtc),即:

(7)

过程4) 设置检测阈值εsingle，用于判断识别目标是否进入了机动段。

若Δβ(ta+Δtc)>εsingle，可认为在ta+Δtc时刻，目标已经进入了新的机动段。记录下[t0,ta]时间里单潜标数据处理得到的目标运动部分要素，然后以ta+Δtc作为新的运动段起点t0′，返回过程1)。

过程5) 重复上述循环处理过程，直至完成了所有目标观测数据的处理。

需要特别指出的是，过程2)中的时间窗Δtc和检测阈值εsingle不能随意选取，具体取值需根据潜标观测的精度等因素确定。如图5所示，目标在tZ时刻进入了新的运动状态。可以看出，在tZ之前，预测值与观测数据之间偏差值较小；在tZ之后，由于目标机动，目标真实位置与机动前的轨迹偏离越来越大，预测值与真实值之间的偏差也逐渐积累增大。如果Δtc设置过小，可能无法积累较大的偏差，进而导致无法检测出目标运动状态的变化；如果Δtc设置过大，会增大[ry(t0)-y]/vy和vx/vy的估计误差，并将其误判为目标运动状态的变化。同理，选择检测阈值εsingle时，过大可能导致无法识别出目标运动状态的变化，过小会将[ry(t0)-y]/vy和vx/vy的估计误差误判为目标运动状态的变化。

从图5中还可以看出，tZ与ta、(ta+Δtc)之间没有确定性关系。因此，单个潜标仅能判断目标在某一时间段里是否进入了新的运动段，无法计算目标机动时刻tZ，需要在多潜标数据集中融合过程中对tZ进行估计。

图5 潜标对机动目标的识别处理过程曲线

通过SSBDIP算法，每个潜标便可利用观测数据得到数组[ry(t0)-y]/vy和vx/vy，每组对应一段匀速直线运动。

4 多潜标数据集中融合算法设计

单潜标的目标观测数据通过SSBDIP算法处理后，得到了目标运动的2个非时敏要素：[ry(t0)-y]/vy和vx/vy，并通过卫星链路将目标运动的2个要素，以及运动段起始时刻t0、运动段结束时刻tend、运动分段起始时刻的角度观测值β(t0)、潜标位置坐标(x,y)等数据回传至指挥所。将单个潜标传回的一个运动分段的数据记为Xi(i用于区分不同潜标或者同一潜标观测到的多个运动段)，有：

(8)

4.1 非机动目标的多潜标数据集中融合算法

对于非机动目标，MSBDCF算法计算过程如下。

不同潜标估计的vx/vy之间存在偏差。一般来说，潜标对目标持续观测时间越长，对目标运动部分要素的估计精度就越高，可以通过时间加权平均的方式对vx/vy进行融合处理。设指挥所收到了N个潜标的数据，有：

(9)

式(9)中，Vx/Vy表示融合结果，可解算出目标航向Km：

Km=arctan(Vx/Vy)

(10)

2) 目标位置和速度的计算

由于各潜标相互独立且彼此之间不通信，各个潜标回传至指挥所的数据基于自身观测数据和其首次发现目标的时间，因而需要进行时空对准。

对准方式为：在接收到多个潜标回传数据Xi后，定义最先感知到目标信息的时刻作为统一的时间基点，并对各潜标首次观测到目标的时刻校准。校准后，用T1,T2,…,TN分别表示第1个、第2个，…，第N个潜标发现目标时刻，其中，T1=0。用x1,x2,…,xN和y1,y2,…,yN分别代表这些潜标位置的横坐标和纵坐标。

对准后，各潜标的[ry(t0)-y]/vy应改写为[ry(Ti)-yi]/vy。

利用式(2)对[ry(Ti)-yi]/vy转化，有：

(11)

用ai表示[ry(Ti)-yi]/vy，则式(11)可看作二元一次方程：

就知识主题最晚结束年级而言，中国在4个知识主题设置时间均明显前置．相较于六国平均水平，“未知数、变量的使用”“代数式概念”“代数式的运算”“代数式的证明”设置时间较之六国均前置两年．从一定程度上可以说明中国代数思维课程内容设置时间相对较短．

ry(T1)+vy(Ti-T1-ai)=yi

(12)

式(12)中仅有ry(T1)和vy2个未知数。当指挥所接收到2组以上的数据Xi，便可以利用最小二乘法进行求解。

另外，可能出现一种特殊情况，观测到目标的所有潜标正好处于一行，即y1=y2=…=yN。此时式(12)退化为：

ry(T1)+vy(Ti-T1-ai)=y1

(13)

这种情况下，即使获取了多组Xi，式(13)也无法求解。此时，需要利用式(1)对式(12)进行变换，整理后得：

(14)

式(14)中，仅有rx(T1)和vy等2个未知数。而且，对于不同位置的若干潜标，如果y1=y2=…=yN，则x1,x2,…=,xN必然互不相等。指挥所获得2个或以上Xi后，便可进行求解。

获得ry(T1)或rx(T1)和vy后，利用式(1)以及Vx/Vy，可以计算得到rx(T1)或ry(T1)和vx，再通过式(2)可以计算任意时刻目标位置。

4.2 机动目标的多潜标数据集中融合算法

对于机动目标，潜标对目标观测数据通过SSBDIP算法，可得到若干组数据Xi，每组Xi对应一段匀速直线运动。多潜标数据集中融合过程中，需要先剔除掉可能存在的严重偏差，再将不同潜标对同一运动段的Xi关联起来，最后利用4.1节公式解算各匀速直线运动段的全部运动要素。

严重偏差是一种特殊情况，如图6所示。单潜标对机动目标识别处理时，假定目标处于匀速直线运动，利用观测数据处理获得目标运动的2个要素：[ry(t0)-y]/vy和vx/vy，待运动要素估计收敛到一定精度范围后，进行目标角度预测。如果目标刚进入潜标探测范围后不久，在[ry(t0)-y]/vy和vx/vy收敛前就进行了机动(图6中实线)，那么潜标不仅无法识别此变化，反而会将其目标运动状态变化前后的轨迹拟合成一段匀速直线运动，导致了严重偏差估计的出现(图6中虚线)。

图6 严重偏差的示意图

显然，此时单潜标数据处理得到的Xi与目标运动是不匹配的，该组数据会影响整体的解算效果，应进行剔除。

考虑到潜标阵的主要探测对象为潜艇，通常不会进行频繁的机动。因此，可以采用以下办法剔除严重偏差：

指挥所将收到的各潜标传回的目标各个运动段的数据Xi，按运动段起始时间Ti进行排序。依次选取时间上相邻的2组数据，按照4.1节中的计算方法进行融合解算。将这些解算结果中的vy分别记为：vy(1,2),vy(2,3),…,vy(i-1,i),，括号中的数字下标表示参与解算数据Xi的序号。

过程1) 取vy(1,2),vy(2,3),…,vy(i-1,i)中相邻的2个数值，计算其相对平均偏差δk为：

(15)

过程2) 设置判断阈值εmulti用于评价的大小。如果δk≤εmulti，说明vy(k-1,k)和vy(k,k+1)大小接近，参与解算的3组数据Xk-1、Xk和Xk+1属于同一段匀速直线运动。返回到过程1)，计算δk+1。

过程3) 如果δk>εmulti，说明vy(k-1,k)和vy(k,k+1)大小相差较大，参与解算的3组数据Xk-1、Xk和Xk+1不属于同一段匀速直线运动。考虑到过程2)中循环处理的条件，不妨假设δk-1≤εmulti，那么Xk-1和Xk属于一段匀速直线运动分段。因此，Xk+1有2种可能情况：①Xk+1属于一段新的匀速直线运动；②Xk+1是严重偏差估计。

过程4) 计算δk+2并利用εmulti评价其大小。

如果δk+2≤εmulti，说明Xk+1、Xk+2和Xk+3属于一段新的匀速直线运动；

如果δk+2>εmulti，说明Xk+1、Xk+2和Xk+3不属于同一段匀速直线运动。对于Xk+1来说，没有与其属于同一段匀速直线的其他数据。考虑到潜艇目标不会进行频繁的机动，可认为Xk+1就是需要剔除的严重偏差估计。

过程5) 完成对Xk+1的判断后，回到过程1)，计算δk+3。

通过上述过程，完成了剔除严重偏差估计和数据关联工作。关联后的数据，按照4.1节的方法便可解算各段匀速直线运动的要素。

另外，单潜标目标观测数据通过SSBDIP算法，无法给出目标机动时刻tZ。解算出各个匀速直线运动分段的要素后，取相邻2段运动的轨迹延长线交点作为前一运动段的终点(同时也是后一运动段的起点)，计算目标到达该位置的时刻，作为tZ的估计结果。

5 仿真实验与分析

5.1 仿真实验想定

假设，某特定海域部署了5行×5列，共计25枚潜标组成的潜标阵，相邻两行、相邻两列间的距离均为10 km，如图7所示。潜标的探测能力，参考文献[15]，取探测范围半径Dr=8 km，观测周期T为1 s，观测误差假定为高斯分布，均方差取1.5°，均值为0。生成4条模拟潜艇航迹，航迹参数见表1。SSBDIP算法中，时间窗Δtc取200 s，检测阈值εsingle取2.5°；MSBDCF算法中，判断阈值εmulti取0.1。

表1 模拟潜艇航迹参数

5.2 仿真实验分析

利用Matlab软件，对每条航迹各进行50次蒙特卡洛仿真实验。图8～图11给出了算法对4条航迹的仿真航迹和位置估计误差变化曲线。

图8 航迹1的实验结果

图9 航迹2的实验结果

图10 航迹3的实验结果

图11 航迹4的实验结果

表2给出了算法对目标机动的识别结果以及误差，表2中数据四舍五入省略了小数部分，潜艇目标实际机动位置为(19 995 m,11 311 m)。

表2 目标机动识别结果以及误差

在大幅度转向运动目标(航迹2)情况下，算法的位置估计误差较大；在小幅度转向运动目标(航迹3)情况下，算法的位置估计误差较小；变速不变向运动目标(航迹4)情况下，算法的位置估计误差介于前面2种情况之间。具体分析如下：① 从图9可见，对大幅度转向运动目标，算法在目标机动时刻(2 000 s)附近的位置估计误差很大，是因为算法对目标机动时刻tZ的估计不够准确引起，之后呈逐步下降趋势。② 从图10可见，目标机动后，算法的位置估计误差稍逊于匀速直线运动状态下的目标，但位置估计误差的变化幅度没有超过400 m。③ 从图11可见，对于变速不变向运动目标，位置估计误差在机动时刻附近很大，但由于目标航向未发生变化，因此位置估计误差在一段时间后迅速降低。

6 结论

1) 对匀速直线运动目标的位置估计。在设定的仿真条件下，在潜标阵观测范围内，算法的位置估计误差在200 m以内。

2) 对机动目标的识别判断。3种不同运动方式的机动目标，其算法均能够识别出目标机动，但机动时刻的估计误差较大。实际上，减小时间窗Δtc可以降低tZ的估计误差，但对于小幅度转向机动目标等情况，可能无法识别出目标机动。

3) 通过仿真实验可见：建立的SSBDIP-MSBDCF的综合算法，实现了潜标阵纯方位目标运动分析从“部分可观测”到“完全可观测”的突破，综合算法的稳定性和精度良好。