0 引　言

基于纯方位的目标运动分析（Target Motion Analysis，TMA）因观测信息量少、原理简单，在军事领域应用广泛，但需要水下观测平台至少机动一次，在实际应用中该可观测性条件并不经常满足，使用条件相对苛刻。目标辐射噪声中常含有稳定的线谱，当观测平台不作机动时，纯方位解算失效，但是在测得目标方位的基础上联合波束域的线谱信息，仍然具备解算目标运动要素的条件，因此，在观测平台保持匀速直航时，目标方位联合频率的TMA方法具有较强的应用价值。

国内外学者作了诸多相关研究，Passerieux等^[1]较早进行了方位联合频率的目标运动分析研究；刘健等^{[2 − 3]}对基于方位频率TMA的可观测性进行分析，通过伪线性化模型，证明了此系统可观测，并对方位-频率的伪线性卡尔曼滤波算法进行分析；李关防等^[4]建立高精度瞬时频率跟踪提取与方位-频率的定位跟踪模型，直接利用最小二乘法求解，可以在解算收敛的同时得到目标辐射源频率；赵建昕等^[5]建立伪线性模型，利用最小二乘法构建了2个伪线性滤波器，然后将二者融合得到解算结果，其特点只需测得瞬时频率即可。近年来，一些文献中还出现了方位频率联合解算中采用粒子群算法、遗传算法，官善政等^[6]基于最小均方误差（Minimum Mean Squared Error，MMSE）准则建立参数估计方程，并依靠粒子群优化（Particle Swarm Optimization，PSO）算法确定一组可使均方误差函数最小的运动参数，实现对目标实时位置、航速、正横距离的精确估计；王潇^[7]基于多普勒频移效应的目标参数解算模型，提出了精英-量子遗传算法，有效地提高了算法的收敛速率。另外，利用多普勒频移信息，在机动目标跟踪^[8]、航路选择^[9]等方面研究也取得了进展^[10]。

以上文献的前提都是先获得准确的线谱频率，然后再联合方位和频率测量信息解算目标运动要素。值得注意的是，在海洋环境噪声干扰、平台自噪声干扰和多普勒运动的影响下，上述传统方法存在2个比较棘手的现实问题：一是目标波束域线谱常出现时有时无、断断续续的现象，无法获得稳定准确的线谱频率，尤其是在信噪比较低的情况下，对接收线谱频率的估计误差很大，导致目标运动要素解算时间很长，甚至常常发散无法收敛；二是目标线谱往往不是唯一存在的，测量平台自噪声中也含有线谱，在众多线谱中迅速确认目标线谱和干扰线谱也并不容易，这就无法保证解算结果的正确性和可靠性。针对上述问题，本文提出一种新的解决方法，在纯方位算法求解相对运动要素的基础上，通过对目标波束域数据进行多普勒搜索补偿和相干累积，随着累积时间增长，当目标要素与实际情况完全匹配时，只有目标线谱能得到相干叠加，而噪声和其他干扰线谱均无法得到增强，根据出现的相干峰即可迅速确定目标运动要素。

1 纯方位求解相对运动要素

经典纯方位TMA仅利用观测平台测量到的一组目标方位来估计目标的运动参数。假定目标和观测平台均做匀速直线运动，只考虑在其运动平面的二维情形，如图1 所示，以地理正北方向为 y 轴，地理正东方向为x轴，观测平台初始位于坐标系的原点o。

记t=0时刻目标距离为D₀(m)，目标航向为H_t(°)，航速为V_t(m /s)，观测平台航向为C_m(°)，航速为V_m(m /s)，(F₀, F ₁, ... , F _i ) 为观测到的目标方位序列。记V_tx为目标速度的x轴分量；V_ty为目标速度的y轴分量；x₀为目标初始位置的x轴分量；y₀为目标初始位置的y轴分量；V_mx为观测平台速度的x轴分量；V_my为观测平台速度的y轴分量；$ {V_{rx}} = {V_{tx}} - {V_{mx}} $为目标相对于观测平台速度的x轴分量；$ {V_{ry}} = {V_{ty}} - {V_{my}} $为目标相对于观测平台速度的y轴分量，则(F₀，D₀，V_t，H_t)可以唯一的确定目标的运动轨迹，根据相对运动几何关系，经过t_i时间后，目标方位与上述参数之间关系为：

$ \begin{split} {\text{tan}}({F_i}) =& \frac{{{x_0} + ({V_{tx}} - {V_{mx}}){t_i}}}{{{y_0} + ({V_{ty}} - {V_{my}}){t_i}}} = \frac{{{D_0}\sin {F_0} + {V_{rx}}{t_i}}}{{{D_0}\cos {F_0} + {V_{ry}}{t_i}}} , \end{split} $

(1)

式中对于i=0,1,…,k，t_i、F_i为已知量，可构成关于(D₀，V_rx，V_ry)的非线性方程组为：

$ \left[ \begin{gathered} \begin{array}{*{20}{c}} {\sin ({F_1} - {F_0})}&{ - \cos {F_1}{t_1}}&{\sin {F_1}{t_1}} \end{array} \\ \begin{array}{*{20}{c}} {\sin ({F_2} - {F_0})}&{ - \cos {F_2}{t_2}}&{\sin {F_2}} {t_2}\end{array} \\ \begin{array}{*{20}{c}} {}& \cdots &{\begin{array}{*{20}{c}} {}&{}&{\begin{array}{*{20}{c}} {}&{ \cdots \begin{array}{*{20}{c}} {}&{}& \cdots \end{array}} \end{array}} \end{array}} \end{array} \\ \begin{array}{*{20}{c}} {\sin ({F_k} - {F_0})}&{ - \cos {F_k}{t_k}}&{\sin {F_k}{t_k}} \end{array} \\ \end{gathered} \right]\left[ {\begin{array}{*{20}{c}} {{D_0}} \\ {{V_{rx}}} \\ {{V_{ry}}} \end{array}} \right] = 0 。$

(2)

当观测平台等速直航的情况下，理论研究表明，该系统不完全观测，即无法通过关于式(2)的一系列非线性方程组唯一的确定目标参数向量$ {X_0} = ( {{D_0},{V_{rx}}, {V_{ry}}} ) $，但此情况却部分可观测^[11]，若记$ {X_0} = \displaystyle\left( {\frac{{{V_{rx}}}}{{{D_0}}},\displaystyle\frac{{{V_{ry}}}}{{{D_0}}}} \right) $，则可通过式(2)得到关于$ {X_0} $的方程：

$ {\boldsymbol{A}}{X_0} = b 。$

(3)

其中：

$ \begin{split}&{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} {\displaystyle\sum\limits_{i = 1}^k {{t_i}^2{{\cos }^2}{F_i}} }&{ - \displaystyle\sum\limits_{i = 1}^k {{t_i}^2\sin {F_i}\cos {F_i}} } \\ { - \displaystyle\sum\limits_{i = 1}^k {{t_i}^2\sin {F_i}\cos {F_i}} }&{\displaystyle\sum\limits_{i = 1}^k {{t_i}^2{{\sin }^2}{F_i}} } \end{array}} \right] \text{，} \\ &b = \left[ {\begin{array}{*{20}{c}} {\displaystyle\sum\limits_{i = 1}^k {{t_i}\cos {F_i}\sin \left( {{F_i} - {F_0}} \right)} } \\ { - \displaystyle\sum\limits_{i = 1}^k {{t_i}\sin {F_i}\sin \left( {{F_i} - {F_0}} \right)} } \end{array}} \right] 。\end{split}$

当观测平台、目标和观察者在一条直线上运动或同向平行运动时，此时$ {F_1}{\text{ = }}{F_2} \cdots = {F_k} $，方位无变化，$ {X_0} $无法得到，其他情况，系统可观测，$ {X_0} = \left( {\displaystyle\frac{{{V_{rx}}}}{{{D_0}}},\displaystyle\frac{{{V_{ry}}}}{{{D_0}}}} \right) $可以唯一确定，目标相对速度与初距比值$ \displaystyle\frac{{{V_r}}}{{{D_0}}} $和相对航向为$ {H_r} $可表示为：

$ \left\{ \begin{gathered} k = \frac{{{V_r}}}{{{D_0}}} = \sqrt {{{\left( {\frac{{{V_{rx}}}}{{{D_0}}}} \right)}^2} + {{\left( {\frac{{{V_{rx}}}}{{{D_0}}}} \right)}^2}} ，\\ {H_r} = {\tan ^{ - 1}}(\frac{{{V_{rx}}}}{{{V_{ry}}}}) 。\\ \end{gathered} \right. $

(4)

可知，在观测平台等速直航条件下，通过纯方位算法只可以获得目标的相对运动要素，如果在观测平台保持匀速直航时求解目标的绝对运动要素，那就必须引入其他测量信息，常见的有噪声中的目标线谱、传播时延等。

2 目标运动要素快速解算原理 2.1 线谱相干累积原理

目标辐射噪声常含有线谱信息，由于观测平台和目标之间的相对运动，观测平台声呐接收到的线谱信息存在多普勒频移，如果从目标所在方位的波束域数据中提取到线谱，可以得到任意时刻$ {t_i} $的频率测量方程为：

$ {f_i} = {f_0}\left(1 - \frac{{{V_r}\cos {\theta _{ri}}}}{c}\right)。$

(5)

式中：$ {f_0} $为目标辐射噪声中的线谱频率；$ {V_r} $为目标的相对速度；$ {\theta _{ri}} $为目标方位角与目标相对航向的夹角；c为声速。

记$ {t_i} $时刻的多普勒因子：

$ {\hat \beta _i} = - {{{V_r}\cos {\theta _{ri}}} \mathord{\left/ {\vphantom {{{V_r}\cos {\theta _{ri}}} c}} \right. } c} 。$

(6)

受目标航行工况和水声传播环境影响，平台接收到的目标线谱噪声存在能量起伏，有时会被淹没在宽带噪声中。目标线谱与宽带噪声的瞬时信噪比过低导致线谱瞬时频率$ {f_i} $估计误差较大，对目标运动分析结果存在不利影响。为了充分利用线谱信息，可以对接收信号进行时域相干累积，提高信噪比。由于目标线谱瞬时频率随相对态势变化产生多普勒频移，为了实现相干累积，需要对接收信号进行实时多普勒补偿，将接收线谱瞬时频率$ {f_i} $补偿至$ {f_0} $。由于缺少目标线谱的先验信息，无法在频域直接补偿，故选择在时域进行数据重采样，达到多普勒补偿的目的。设接收数据的采样率为$ {f_s} $，$ {t_i} $时刻的重采样频率$ f_s^{'} $，则二者之间应满足如下关系^[12]：

$ f_s^{'} = {f_s}(1 + {\beta _i})。$

(7)

设观测平台声呐输出目标连续辐射噪声，按声呐方位输出时刻$ {t_i} $将对应的波束域数据分段为$ {x}_{i}\left(\tau \right), 0\leqslant \tau < T,i=1,2,\cdots N $，每段数据的时间长度T等于方位采样的时间间隔，根据式(7)对各段信号进行时域重采样处理，得到多普勒补偿后的信号$ {y_i}(\tau ) $，进行FFT处理得到其频谱$ {Y_i}(f) $。经过重采样后多普勒得到补偿，各段数据的线谱频率均为$ {f_0} $，此时对$ {Y_i}(f) $进行频域相干累积，计算公式为：

$ Y(f) = \sum\nolimits_{i = 1}^N {{Y_i}(f)\exp ( - 2{\text{π}} jf{t_i})}。$

(8)

式中$ Y(f) $将在$ f = {f_0} $处出现峰值。当进行多普勒补偿时代入的$ {\beta _i} $与实际情况有较大误差时，重采样后的线谱频率不能准确补偿至$ {f_0} $，线谱相干性降低，$ Y({f_0}) $的峰值能量随之降低。对各段数据同时计算不同多普勒补偿后$ Y({f_0}) $的值，也就是进行多普勒匹配搜索，然后经过频域相干累积，当$ Y({f_0}) $最大时对应的多普勒补偿量就是目标的多普勒因子。

然而，直接按照以上步骤处理会存在多普勒因子搜索范围无法确定的问题，若在很大的范围内搜索会导致效率很低，因此，快速估计多普勒因子是提高解算收敛速度的关键。

2.2 基于Radon变换的多普勒因子快速估计

根据多普勒效应，以观测平台为参考，多普勒因子实际上只与目标的径向速度有关。将接收信号可能的多普勒因子（对应径向速度）范围划分为L个区间，记为$ {\beta _l} $，$ l = 1,2, \cdots L $。多普勒因子快速估计（见图2）步骤为：

步骤1　利用插值进行多普勒补偿。记第i个声呐方位时刻为$ {t_i} $，声呐波束域数据为$ {x_i}\left( \tau \right) $，依据多普勒因子$ {\beta _l} $进行插值，得到L个信号$ {y_{il}}(\tau ) $及其频谱$ {Y_{il}}({f_k}) $，$ {f_k} $为线谱及其多普勒范围内的频率点。

步骤2　快速估计目标多普勒因子。将$ {Y_{il}}({f_k}) $按$ {f_k} $进行累加得到：

$ {Y_{il}} = \sum\nolimits_k {{Y_{il}}({f_k})}。$

(9)

$ {Y_{il}} $刻画了目标多普勒因子随时间的变化，在较短时间内，目标多普勒因子随时间近似线性变化，$ {Y_{il}} $的峰值在时间—多普勒二维平面上为一条直线，对$ {Y_{il}} $进行Radon变换，记为$ {R_y} $，则$ {R_y} $的峰值坐标$ (d,\theta ) $与$ {t_i} $时刻目标多普勒因子之间$ {\beta _i} $存在如下关系：

$ {\beta _i} \cdot \cos \theta + {t_i} \cdot \sin \theta = d 。$

(10)

根据式(10)可以估算目标的多普勒因子$ {\hat \beta _i} $的近似值，该值具有相对较高的可信度，下一步可在此值周围继续精确搜索。若目标存在多条线谱时，不同线谱的时间-多普勒变化规律相同，可以将所有线谱按式(9)进行累加。

2.3 目标运动要素快速匹配

目标多普勒因子取决于运动速度和航向2个未知量，目标在不同的航向航速下可能有相同的多普勒因子，在这2个维度上同时进行搜索补偿计算量较大。在此，提出一种基于径向相对速度和相对航向的分步搜索估计算法，即在每一轮搜索中只遍历径向相对速度或相对航向中的一个维度，缩小搜索范围并提高计算速率。时域相干累积原理图如图3所示。

根据上述估计的目标多普勒因子$ {\hat \beta _i} $计算目标的相对航向和相对航速，目标相对航速$ {V_r} $及相对舷角$ {\theta _{ri}} $应满足$ {\hat \beta _i} = - {{{V_r}\cos {\theta _{ri}}} \mathord{\left/ {\vphantom {{{V_r}\cos {\theta _{ri}}} c}} \right. } c} $。一般来说$ {\hat \beta _i} $的大小主要取决于目标径向速度，$ {\hat \beta _i} $的变化率主要取决于目标的航向，在解算初期目标舷角变化相对较小，多普勒变化不大，对目标径向速度的估计精度远大于目标航向的估计精度。因此，将估计出的目标径向相对速度$ {V_r}\cos {\theta _{ri}} $作为已知条件，对目标相对航向作为未知条件进行搜索，将目标相对航向的搜索范围划分为R个区间，每个区间的目标航速$ {V_r} $和航向$ {H_r} $均有：

$ {\hat{\beta }}_{0}=-{V}_{r}\mathrm{cos}({H}_{r}-{F}_{0})/c ,\;r=1,2,\;\cdots,\;R。$

(11)

式中：$ {\hat \beta _0} $为目标初始多普勒估计值；$ {F_0} $为目标初始方位。结合目标方位$ {F_i} $，计算目标多普勒因子随时间的变化$ {\beta }_{ir}=-{V}_{r}\mathrm{cos}({H}_{r}-{F}_{i})/c $。

根据$ {\beta _{ir}} $对接收信号$ {x_i}\left( \tau \right) $进行插值得到$ {y_{ir}}(\tau ) $及其频谱$ {Y_{ir}}({f_k}) $，将$ {Y_{ir}}({f_k}) $按时间积累：

$ {Y_r}({f_k}) = \sum\nolimits_i {{Y_{ir}}({f_k})}。$

(12)

式中：i为数据累积的时间序号。在$ {f_k} = {f_0} $处，$ {Y_r}({f_k}) $有极大值，且$ {Y_r}({f_0}) $最大值对应的$ {V_r} $和$ {H_r} $分别为目标相对航速和相对航向的匹配值，此时目标多普勒因子也得到了准确搜索。当目标存在多条线谱时，可对不同频率的结果进行加权平均。

根据搜索匹配的相对航速$ {V_r} $，进一步结合相对运动关系，可得到目标的运动速度$ {V_t} $：

$ {V_t} = \sqrt {V_r^2 + V_m^2 - 2{V_r}{V_m}\cos ({H_r} - {C_m})} 。$

(13)

目标航向$ {H_t} $：

$ \frac{{{\text{sin}}\left( {{H_t} - {C_m}} \right)}}{{{\text{sin}}({C_m} - {H_r})}} = \frac{{{V_r}}}{{{V_t}}}。$

(14)

再由式(4)获得解算开始时的目标初距：

$ {D_0} = \frac{{{V_r}}}{k} 。$

(15)

3 算法有效性验证

利用某次航行试验的实录船舶辐射噪声数据，结合仿真处理进行算法有效性验证。

基本情况：观测平台航速为4 m/s，航向为0°，目标初始距离为29 km，航速为5.2 m/s，航向为173.5°，目标方位数据1 s更新一次，实测目标方位标准差估计为1.5°左右。由于条件限制，目标噪声只在50 Hz处存在一根弱线谱，线谱频率随目标多普勒变化。为了增强对比性，仿真162.5 Hz的目标线谱和400 Hz的干扰线谱（模拟平台自噪声干扰线谱）添加到辐射噪声中，并按照相对运动生成接收信号。依照线谱频率由大到小的顺序，估算信噪比分别为0、6、3 dB，截取数据的时间长度为600 s，目标接收信号的Lofar（Low frequency analysis and recording）谱如图4所示。

在本实验条件下，若使用传统方法解算目标运动要素将面临2个问题：一是无法事先确定3根线谱是否都是目标线谱，尤其是对于实时信号处理，在较短时间内很难辨别出线谱真伪；二是低信噪比情况下提取的线谱频率误差大，尤其是对于本实验中50 Hz的低频线谱，频率误差可能都要比多普勒频移量大，直接使用此频率往往导致解算精度较低或者结果发散。

利用前100 s数据进行目标多普勒估计，多普勒搜索范围为相对速度2~10 m/s，搜索步长0.020 m/s，将目标信号40~700 Hz范围内的频谱叠加并进行Radon变化，结果如图5所示，顶点位置为（−16，−2.6°），图6显示了前100 s目标多普勒因子的仿真值和估计值，初始多普勒因子约0.381%，初始径向相对速度约为5.72 m/s。由此可见，仅利用了前100 s的数据，便可获得较好的目标多普勒因子估计值。

根据前述获得的目标径向相对速度5.72 m/s为参考，在舷角±80°范围以内由式(11)和式(12)进一步精确搜索匹配，根据不同线谱进行相干累积，图7显示了3条线谱的累积结果，并与无线谱时的相干累积结果进行了对比。可以看出，在相应舷角处，只有目标线谱随时间相干积累，峰值增大，而干扰线谱无随时间累积增强的现象。

根据相干累积匹配的相对航速和航向结果计算目标运动要素，图8～图10分别显示了根据2根目标线谱的累积结果估算的目标速度、航向和初距结果，由结果可知，线谱处的信噪比越高解算收敛时间越短，而且高频线谱由于多普勒相对明显，解算结果也更快更稳定。对不同线谱的计算结果进行加权平均可得到更优的综合估计结果，图11显示了2根线谱按信噪比加权平均后的估计结果，相比之下，加权平均后的收敛速度更快，在200 s时目标运动要素就已经稳定收敛。

综合以上分析处理过程可知，对于存在多线谱的情况，只有目标线谱才具有和实际情况一致的多普勒频移，只要目标运动要素匹配一致，目标线谱能量就会随处理时间延长而不断累积，体现为同频相干叠加，而干扰线谱则是随机变化的，无法得到持续累积增强，因此，随着处理时间的延长，信噪比条件会逐渐改善，同时目标线谱和干扰线谱也得以分辨。

值得注意的是，由于不同线谱的信噪比水平存在差异，通过对多根线谱进行加权融合综合处理，如表1所示，与单根线谱的解算结果相比，解算误差明显减小，可获得更好的解算精度。

4 结　语

当观测平台作匀速直线运动时，在低信噪比和线谱干扰条件下，针对传统解算目标运动要素无法获得稳定准确的线谱频率的问题，通过多普勒匹配搜索补偿和相干累积方法，可以快速匹配得到目标运动要素。本方法无需事先辨认目标和干扰线谱，且具有适应低信噪比条件和解算收敛速度快的特点，具有较强的应用价值。在实际情况下，由于海洋信道的频率选择性衰落，目标线谱往往较弱而干扰线谱则较强，在某些态势下，利用Radon变换进行多普勒因子快速估计时，可能会存在峰值模糊现象，影响多普勒因子初值的估计精度。