CN105242262A - 一种基于天线周期扫描时间差无源定位方法 - Google Patents

Abstract

一种基于天线周期扫描时间差无源定位方法，本发明涉及基于天线周期扫描时间差无源定位方法。本发明是为了解决目前采用的无源定位方法中存在跟踪收敛速度慢、精度低以及需要一定机动的问题。本发明建立了平台和目标之间坐标系，并根据所建立坐标系及转换关系，获得目标位置的一种计算方法，利用探测范围的限制对所得目标位置的解进行选取，最后对带有测量误差的定位结果采用递推最小二乘法进行优化，得到了优化后位置估计。本发明仅利用天线扫描时间观测值及平台状态即可完成对目标的定位，且定位精度随时间测量精度提高而提高，是一种新的定位方式，对目前的无源定位体制进行了有效补充。本发明应用于无源定位领域。

Description

一种基于天线周期扫描时间差无源定位方法

技术领域

本发明涉及基于天线周期扫描时间差无源定位方法。

背景技术

无源雷达在应用中具有隐蔽性好，探测隐身和低空目标能力较强的优点，但由于其无源特性导致无法对静默目标有效定位，同时也带来了测量精度的问题。一般研究的无源定位方法有基于到达时间法和测向线交叉定位法。基于到达时间法中虽然对到达时间的精确测量能够反映目标径向距离的变化，但是该方法存在着速度慢、精度低的弊端达时间；测向线交叉定位法则要求观测平台有一定的机动飞行，这为该技术的应用带来诸多的不便，此外,跟踪收敛速度慢、精度低也是该方法的严重缺点。

雷达天线扫描周期是一项重要的雷达技术参数，天线扫描的测量方法较为简单，过去一般通过测量相邻两次雷达照射的时间差来确定天线扫描周期。当目标相对观测平台进行有一定规则的运动时，时间间隔将产生相应的规则变化，因此通过天线周期扫描的时间差可求得雷达状态信息，但目前利用此项方法进行定位的研究比较少见。

发明内容

本发明是为了解决目前采用的无源定位方法中存在跟踪收敛速度慢、精度低以及需要一定机动的问题，而提出的一种基于天线周期扫描时间差无源定位方法。

一种基于天线周期扫描时间差无源定位方法按以下步骤实现：

步骤一：建立坐标系并得到各坐标系间的转换关系；

建立观测平台坐标系，以观测平台运动方向为y轴方向，垂直y轴为x轴，y轴正方向右侧为x正方向；

建立目标坐标系，将观测平台的速度等效到目标上，采用等效后目标运动方向为y'轴方向，垂直y'轴为x'轴，y'轴正方向右侧为x'正方向，此时的坐标系相当于目标雷达对观测平台的定位坐标系，可以得到观测平台相对于目标的位置(x',y')，根据观测平台坐标系与目标坐标系的转换关系得到坐标变换x＝x'和y＝y'，得到目标相对于观测平台在平台坐标系中的位置(x,y)；

建立大地坐标系，以正北方位Y轴方向，正东方为X轴方向；当目标坐标需要转换到大地坐标下时，根据平台坐标系与大地坐标系的转换关系，即公式(1)-(3)，可得到大地坐标系下的目标坐标；

γ是观测平台坐标系y方向和大地坐标系Y方向的夹角；是目标在平台坐标系中的方位角，X_T为目标在大地坐标系中X方向位置，Y_T为目标在大地坐标系中Y方向位置，x₀为观测平台在大地坐标系中X方向位置，y₀为观测平台在大地坐标系中Y方向位置；

步骤二：根据步骤一所述的平台坐标系与目标坐标系的转换关系以及大地坐标系与平台坐标系的转换关系，可以将目标坐标系中所求观测平台坐标转换为大地坐标系中目标坐标，直接在目标坐标系中对观测平台位置进行求解，等效后目标进行匀速直线运动时，可以得到求解方程中x,y,T三个未知量需要的方程为：

$\{\begin{array}{c}2\mathrm{\π}({\mathrm{\τ}}_i-T)\[x^2+y^2-(y_{i+1}+y_i)y+y_{i+1}{y}_i\]-xT(y_{i+1}-y_i)=0\\ 2\mathrm{\π}({\mathrm{\τ}}_{i+1}-T)\[x^2+y^2-(y_{i+2}+y_{i+1})y+y_{i+2}{y}_{i+1}\]-xT(y_{i+2}-y_{i+1})=0\\ 2\mathrm{\π}({\mathrm{\τ}}_{i+2}-T)\[x^2+y^2-(y_{i+3}+y_{i+2})y+y_{i+3}{y}_{i+2}\]-xT(y_{i+3}-y_{i+2})=0\end{array}---\left(9\right)$

所述T为天线扫描周期，τ_i＝t_i+1-t_i，t_i为第i时刻天线扫描时间测量值，t_i+1为第i+1时刻天线扫描时间测量值，y_i为目标坐标系中第i时刻平台的y'方向坐标，y_i+1为目标坐标系中第i+1时刻平台的y'方向坐标，y_i+2为目标坐标系中第i+2时刻平台的y'方向坐标，y_i+3为目标坐标系中第i+3时刻平台的y'方向坐标；

步骤三：求解公式(9)的方程组得到6组数学解，解的选取利用探测范围的限制，所述探测范围为100～350km，在探测范围内的解记为可行解，进一步筛选可行解集合时，若有唯一解则将其作为解，若无可行解则将解计为零，有多组解可以取各解的平均值；

步骤四：求解递推最小二乘法的参数R_i、和P₀，将遗忘因子与R_i相乘得到R_i'，将R_i'、和P₀带入递推最小二乘的基本公式，利用通过时间测量值计算得到的z_i估计目标的真实位置得到估计值所述z_i为第i时刻的观测向量(x_i,y_i,T_i)^T，为第i时刻的估计值矩阵T_i为第i时刻天线扫描周期计算值，为第i时刻天线扫描周期估计值，为目标在x'方向坐标的估计值，为目标在y'方向坐标的估计值。

本发明提出了一种远距离静止目标在仅有其雷达天线周期扫描时间差并且观测平台处于匀速直线运动的条件下位置估计的方法，先对已知条件进行分析，建立了平台和目标之间坐标系并获得目标位置的一种计算方法，同时对带有测量误差的定位结果采用递推最小二乘法进行优化，得到了位置估计，本发明利用天线周期扫描时间差进行定位，与目前常规无源定位方法相比不需要角度测量值，本发明在平台匀速直线运动时即可实现定位，不需要平台进行特殊形式机动，仅利用天线扫描时间观测值及平台状态即可完成对目标的定位，且定位精度随时间测量精度提高而提高，是一种新的定位方式，对目前的无源定位体制进行了有效补充。

附图说明

图1为目标Tar和观测器O的几何关系示意图；

图2为观测平台坐标系图；

图3为目标坐标系图；

图4为观测平台坐标系和目标坐标系转换关系图；

图5为大地坐标系与观测平台坐标系几何关系图；

图6为解的选取流程图；

图7为目标x方向计算结果图(T＝20S)；

图8为目标y方向计算结果图(T＝20S)；

图9为天线扫描周期计算结果图(T＝20S)；

图10为目标x方向计算误差图(T＝20S)；

图11为目标y方向计算误差图(T＝20S)；

图12为天线扫描周期计算误差图(T＝20S)；

图13为目标x方向计算结果图(增加噪声后T＝20S)；

图14为目标y方向计算结果图(增加噪声后T＝20S)；

图15为天线扫描周期计算结果图(增加噪声后T＝20S)；

图16为目标x方向计算误差图(增加噪声后T＝20S)；

图17为目标y方向计算误差图(增加噪声后T＝20S)；

图18为天线扫描周期计算误差图(增加噪声后T＝20S)；

图19为目标x方向计算结果图(T＝10S)；

图20为目标y方向计算结果图(T＝10S)；

图21为天线扫描周期计算结果图(T＝10S)；

图22为目标x方向计算误差图(T＝10S)；

图23为目标y方向计算误差图(T＝10S)；

图24为天线扫描周期计算误差图(T＝10S)；

图25为间隔一点取点法所得目标x方向计算结果图(T＝10S)；

图26为间隔一点取点法所得目标y方向计算结果图(T＝10S)；

图27为间隔一点取点法所得天线扫描周期计算结果图(T＝10S)；

图28为间隔一点取点法所得目标x方向计算误差图(T＝10S)；

图29为间隔一点取点法所得目标y方向计算误差图(T＝10S)；

图30为间隔一点取点法所得天线扫描周期计算误差图(T＝10S)

图31为间隔两点取点法所得目标x方向计算结果图(T＝10S)；

图32为间隔两点取点法所得目标y方向计算结果图(T＝10S)；

图33为间隔两点取点法所得天线扫描周期计算结果图(T＝10S)；

图34为间隔两点取点法所得目标x方向计算误差图(T＝10S)；

图35为间隔两点取点法所得目标y方向计算误差图(T＝10S)；

图36为间隔两点取点法所得天线扫描周期计算误差图(T＝10S)；

图37为间隔一点取点法所得目标x方向计算结果图(增加噪声T＝10S)；

图38为间隔一点取点法所得目标y方向计算结果图(增加噪声T＝10S)；

图39为间隔一点取点法所得天线扫描周期计算结果图(增加噪声T＝10S)

图40为间隔一点取点法所得目标x方向计算误差图(增加噪声T＝10S)；

图41为间隔一点取点法所得目标y方向计算误差图(增加噪声T＝10S)；

图42为间隔一点取点法所得天线扫描周期计算误差图(增加噪声T＝10S)

图43为递推最小二乘法得到目标x方向计算结果图(增加噪声T＝20S)；

图44为递推最小二乘法得到目标y方向计算结果图(增加噪声T＝20S)；

图45为递推最小二乘法得到天线扫描周期计算结果图(增加噪声T＝20S)；

图46为递推最小二乘法得到目标x方向计算误差图(增加噪声T＝20S)；

图47为递推最小二乘法得到目标y方向计算误差图(增加噪声T＝20S)；

图48为递推最小二乘法得到天线扫描周期计算误差图(增加噪声T＝20S)；

图49为间隔一点取点递推最小二乘法得到目标x方向计算结果图(增加噪声T＝10S)；

图50为间隔一点取点递推最小二乘法得到目标y方向计算结果图(增加噪声T＝10S)；

图51为间隔一点取点递推最小二乘法得到天线扫描周期计算结果图(增加噪声T＝10S)；

图52为间隔一点取点递推最小二乘法得到目标x方向计算误差图(增加噪声T＝10S)；

图53为间隔一点取点递推最小二乘法得到目标y方向计算误差图(增加噪声T＝10S)；

图54为间隔一点取点递推最小二乘法得到天线扫描周期计算误差图(增加噪声T＝10S)；

T为天线扫描周期。

具体实施方式

具体实施方式一：一种基于天线周期扫描时间差无源定位方法包括以下步骤：

考虑二维单站被动跟踪问题，图1给出了目标Tar和观测器O的几何关系示意图，图中T_ar(x,y)表示目标Tar在观测平台坐标系的坐标，(X_T,Y_T)表示目标在大地坐标系中的坐标，X(E)表示X轴为正东方，Y(N)表示Y轴为正北方。

目标Tar处于静止状态，观测平台O沿y方向进行匀速直线运动，观测器被动接收目标Tar辐射的信号，利用对该信号的被动测量对目标进行定位跟踪。目标雷达天线以固定周期进行圆周扫描，观测器则对目标雷达的扫描信号进行探测并记录时间，由于观测平台并非朝向目标运动，观测器所测得的目标雷达天线两次扫描的时间间隔将不等于其扫描周期，其中包含了目标的位置信息。本文研究的即是基于天线周期扫描时间差对目标进行定位的问题。

步骤一：建立坐标系并得到各坐标系间的转换关系；

建立观测平台坐标系如图2所示，以观测平台运动方向为y轴方向，垂直y轴为x轴，y轴正方向右侧为x正方向；

为简化求解，建立目标坐标系如图3所示，将观测平台的速度等效到目标上，采用等效后目标运动方向为y′轴方向，垂直y′轴为x′轴，y′轴正方向右侧为x′正方向，此时的坐标系相当于目标雷达对观测平台的定位坐标系，可以得到观测平台相对于目标的位置(x',y')，根据观测平台坐标系和目标坐标系转换关系如图4所示，得到坐标变换x＝x'和y＝y'，得到目标相对于观测平台在观测平台坐标系中的位置(x,y)；

大地坐标系与观测平台坐标系几何关系如图5所示，以正北方位(Y)y轴方向，正东方(E)为x轴方向；目标坐标系中当需要转换到大地坐标下时，根据观测平台坐标与大地坐标的相对关系即可得到大地坐标下的计算方式：

γ是观测平台坐标系y方向和大地坐标系Y方向的夹角；是目标在平台坐标系中的方位角，X_T为目标在大地坐标系中X方向位置，Y_T为目标在大地坐标系中Y方向位置，x₀为观测平台在大地坐标系中X方向位置，y₀为观测平台在大地坐标系中Y方向位置；

步骤二：根据步骤一所述的平台坐标系与目标坐标系的转换关系以及大地坐标系与平台坐标系的转换关系，可以将目标坐标系中所求观测平台坐标转换为大地坐标系中目标坐标，直接在目标坐标系中对观测平台位置进行求解，等效后目标进行匀速直线运动(即将观测平台的匀速直线运动等效到目标上)，且观测平台不朝向目标方向运动时，可以得到求解方程中x,y,T三个未知量需要的方程为：

$\{\begin{array}{c}2\mathrm{\π}({\mathrm{\τ}}_i-T)\[x^2+y^2-(y_{i+1}+y_i)y+y_{i+1}{y}_i\]-xT(y_{i+1}-y_i)=0\\ 2\mathrm{\π}({\mathrm{\τ}}_{i+1}-T)\[x^2+y^2-(y_{i+2}+y_{i+1})y+y_{i+2}{y}_{i+1}\]-xT(y_{i+2}-y_{i+1})=0\\ 2\mathrm{\π}({\mathrm{\τ}}_{i+2}-T)\[x^2+y^2-(y_{i+3}+y_{i+2})y+y_{i+3}{y}_{i+2}\]-xT(y_{i+3}-y_{i+2})=0\end{array}---\left(9\right)$

所述T为天线扫描周期，τ_i＝t_i+1-t_i，t_i为第i时刻天线扫描时间测量值，t_i+1为第i+1时刻天线扫描时间测量值，y_i为目标坐标系中第i时刻平台的y'方向坐标，y_i+1为目标坐标系中第i+1时刻平台的y'方向坐标，y_i+2为目标坐标系中第i+2时刻平台的y'方向坐标，y_i+3为目标坐标系中第i+3时刻平台的y'方向坐标，τ_i+1＝t_i+2-t_i+1，t_i+2为第i+2时刻天线扫描时间测量值，τ_i+2＝t_i+3-t_i+2，t_i+3为第i+3时刻天线扫描时间测量值；

步骤三：求解公式(9)的方程组得到若干组(一般为6组)数学解，应根据给定条件对解进行选取。对于选取后无解、唯一解和多解的情况应进行区分；

解的选取利用探测范围的限制，当解超过最大探测范围或者小于最小探测范围时记为无效解，在探测范围内的解记为可行解，所述探测范围为100～350km，具体操作如图6所示：

对于选取可行解集合应当予以处理，有唯一解则将其作为解，若无可行解则将解计为零，有多组解可以取各解的平均值。

步骤四：求解递推最小二乘法的参数R_i、和P₀，将遗忘因子与R_i相乘得到R_i'，将R_i'、和P₀带入递推最小二乘的基本公式，对定位数据进行处理，利用通过时间观测值计算得到的z_i估计目标的真实位置得到估计值所述z_i为第i时刻的观测向量(x_i,y_i,T_i)^T，为第i时刻的估计值矩阵T_i为第i时刻天线扫描周期计算值，为第i时刻天线扫描周期估计值，为目标在x'方向坐标的估计值，为目标在y'方向坐标的估计值。

具体实施方式二：本实施方式与具体实施方式一不同的是：所述步骤二中得到求解方程中x,y,T三个未知量需要的三个方程的具体过程为：

在图3θ_i所在三角形中，有如下表达式成立：

${\mathrm{\θ}}_i=arctan\left(\frac{x}{y-y_{i+1}}\right)-arctan\left(\frac{x}{y-y_i}\right)---\left(4\right)$

所述θ_i为第i时刻与第i+1时刻目标与观测平台连线间的夹角；

根据两角和与差的正切公式可以得到：

$tan\left({\mathrm{\θ}}_i\right)=\frac{x(y_{i+1}-y_i)}{x^2+y^2-(y_{i+1}+y_i)y+y_i{y}_{i+1}}---\left(5\right)$

由于目标距离足够远(100～350km)，此时θ_i的值远小于1，θ_i满足的条件，因此可以取：

tan(θ_i)≈θ_i(6)

目标雷达天线以固定周期扫描，可以得到公式(7)：

${\mathrm{\τ}}_i=t_{i+1}-t_i=T+T\frac{{\mathrm{\θ}}_i}{2\mathrm{\π}}---\left(7\right)$

联立公式(6)和公式(7)，消去θ_i后可以得到方程(8)：

2π(τ_i-T)[x²+y²-(y_i+1+y_i)y+y_i+1y_i]-xT(y_i+1-y_i)＝0(8)

求解方程中x,y,T三个未知量需要三个方程，利用连续四个点，每相邻两点为一组，分别代入公式(8)，得到三个方程，方程组为公式(9)。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：所述步骤四中求解递推最小二乘法的参数R_i、和P₀的具体过程为：

R_i的求解过程为：

由于误差在时间测量值t上，对公式(8)两侧求时间t的偏导，可得：

$\begin{array}{c}\frac{dx}{dt}\[4\mathrm{\π}x({\mathrm{\τ}}_i-T)-T(y_{i+1}-y_i)\]+\frac{dy}{dt}\[4\mathrm{\π}y({\mathrm{\τ}}_i-T)-2\mathrm{\π}({\mathrm{\τ}}_i-T)(y_{i+1}+y_i)\]+\\ \frac{dT}{dt}\{-2\mathrm{\π}\[x^2+y^2-(y_{i+1}+y_i)y+y_{i+1}{y}_i\]-x(y_{i+1}-y_i)\}=2\mathrm{\π}({\mathrm{\τ}}_i-T)(2v_{y^{\′}i}y-v_{y^{\′}i}{y}_{i+1}-y_{y^{\′}i}{y}_i)\end{array}---\left(12\right)$

取连续四点，每相邻两点为一组，分别代入公式(12)，得到三个方程，改写为矩阵形式如下：

$\begin{array}{c}\left[\begin{array}{ccc}4\mathrm{\π}x({\mathrm{\τ}}_i-T)-T(y_{i+1}-y_i)& 4\mathrm{\π}y({\mathrm{\τ}}_i-T)-2\mathrm{\π}({\mathrm{\τ}}_i-T)(y_{i+1}+y_i)& -2\mathrm{\π}\[x^2+y^2-(y_{i+1}+y_i)y+y_{i+1}{y}_i\]-x\left(y_{i+1}\right)\\ 4\mathrm{\π}x({\mathrm{\τ}}_{i+1}-T)-T(y_{i+2}-y_{i+1})& 4\mathrm{\π}y({\mathrm{\τ}}_{i+1}-T)-2\mathrm{\π}({\mathrm{\τ}}_{i+1}-T)(y_{i+2}+y_{i+1})& -2\mathrm{\π}\[x^2+y^2-(y_{i+2}+y_{i+1})y+y_{i+2}{y}_{i+1}\]-x(y_{i+2}-y_{i+1})\\ 4\mathrm{\π}x({\mathrm{\τ}}_{i+2}-T)-T(y_{i+3}-y_{i+2})& 4\mathrm{\π}y({\mathrm{\τ}}_{i+2}-T)-2\mathrm{\π}({\mathrm{\τ}}_{i+2}-T)(y_{i+3}+y_{i+2})& -2\mathrm{\π}\[x^2+y^2-(y_{i+3}+y_{i+2})y+y_{i+3}{y}_{i+2}\]-x(y_{i+3}-y_{i+2})\end{array}\right]\\ \left[\begin{array}{c}\frac{dx}{dt}\\ \frac{dy}{dt}\\ \frac{dT}{dt}\end{array}\right]=\left[\begin{array}{c}2\mathrm{\π}({\mathrm{\τ}}_i-T)(2v_{y^{\′}i}y-v_{y^{\′}i}{y}_{i+1}-v_{y^{\′}i}{y}_i)\\ 2\mathrm{\π}({\mathrm{\τ}}_{i+1}-T)(2v_{y^{\′}i}y-v_{y^{\′}i}{y}_{i+2}-v_{y^{\′}i}{y}_{i+1})\\ 2\mathrm{\π}({\mathrm{\τ}}_{i+2}-T)(2v_{y^{\′}i}y-v_{y^{\′}i}{y}_{i+3}-v_{y^{\′}i}{y}_{i+2})\end{array}\right]\end{array}---\left(13\right)$

公式(13)可表示为A_3×3d_3×1＝b_3×1，则有：

$d_{3\×1}={\left(A_{3\×3}^T{A}_{3\×3}\right)}^{-1}{A}_{3\×3}^T{b}_{3\×1}---\left(14\right)$

根据观测误差协方差定义可知，R_i计算式为：

$R_i=\left[\begin{array}{ccc}{\left(\frac{dx}{dt}\right)}^2{\mathrm{\σ}}_t^2& 0& 0\\ 0& {\left(\frac{dy}{dt}\right)}^2{\mathrm{\σ}}_t^2& 0\\ 0& 0& {\left(\frac{dT}{dt}\right)}^2{\mathrm{\σ}}_t^2\end{array}\right]---\left(15\right)$

由公式(14)和公式(15)可得到R_i，为时间测量值的误差的方差；

和P₀的求解过程为：

所述为初始的状态估计值，P₀为初始的状态协方差，若没有给出和P₀，则可以先解公式(9)所示方程，以解得结果作为初始值和P₀，或可以将设置为一充分小的向量，P₀设置为充分大的矩阵，等待最小二乘法自行收敛；

为加快收敛速度，采用前一方法。

$P_0={\left({\stackrel{\‾}{H}}^T{\stackrel{\‾}{R}}^{-1}\stackrel{\‾}{H}\right)}^{-1}---\left(17\right)$

${\hat{X}}_0={\left({\stackrel{\‾}{H}}^T{\stackrel{\‾}{R}}^{-1}\stackrel{\‾}{H}\right)}^{-1}{\stackrel{\‾}{H}}^T{\stackrel{\‾}{R}}^{-1}\stackrel{\‾}{Z}---\left(18\right)$

其中和由起始天线扫描时间测量数据和平台状态数据，所述状态包括位置和速度，所得的相应矩阵作为分块矩阵组合而成，即：

${\stackrel{\‾}{H}}_0={\[H_{3\×3},H_{3\×3},...H_{3\×3}\]}^T---\left(19\right)$

$\stackrel{\‾}{Z}={\left[\begin{array}{cccc}{z}_0^T& z_1^T& ...& z_n^T\end{array}\right]}^T---\left(21\right)$

所述H_3×3表示三阶单位阵；

利用公式(17)-(21)可求得最小二乘递推初始值，即和P₀。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：所述步骤四中得到R_i'的具体过程为：

考虑递推计算过程中,协方差矩阵P(k)随着递推的进程将衰减很快,此时算法中矩阵也急剧衰减,使得新数据失去对估计值的修正能力，因此需要考虑修正方案,以保持新数据对参数估计值一定的修正能力,使得能得到更准确的估计值或适应目标位置的慢速变化。一般修正方法有遗忘因子法和协方差重调法，对于定位在一定误差时数据即失去修正能力的情况，本发明采用遗忘因子进行处理。

采用遗忘因子进行修正处理，将遗忘因子作用在R_i上时如式(16)所示：

R′_i＝e^(l-i)αR_i(16)

其中l是选择用于位置估计的R阵总个数,α是根据经验选取的常数本文中选取为0.05，R_i是用公式(15)计算所得结果，进行递推时将R′_i带入递推最小二乘的基本公式；

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：所述步骤四中递推最小二乘的基本公式为公式(10)和(11)：

由于直接求解的结果仍有一定抖动，而且处理过程中测量数据具有时间顺序，考虑实际应用中对存储容量的限制，本发明采用递推最小二乘法对定位数据进行处理，利用通过时间观测值计算得到的(x_i,y_i,T_i)估计目标的真实位置(x,y,T)，x_i为目标坐标系中第i时刻平台的x'方向坐标，x'为目标坐标系中的横坐标。

${\hat{X}}_{i+1}={\hat{X}}_i+P_i{H}_{i+1}^T{\[H_{i+1}{P}_i{H}_{i+1}^T+R_{i+1}\]}^{-1}\[z_{i+1}-H_{i+1}{X}_i\]---\left(10\right)$

$P_{i+1}=P_i-P_i{H}_{i+1}^T{\[H_{i+1}{P}_i{H}_{i+1}^T+R_{i+1}\]}^{-1}{H}_{i+1}{P}_i---\left(11\right)$

其中表示第i时刻的估计值矩阵z_i表示第i时刻的观测向量(x_i,y_i,T_i)^T，P_i表示第i时刻的协方差矩阵，H_i表示第i时刻的测量矩阵在此处即是3×3的单位阵，R_i表示第i时刻的观测误差协方差；

实施例一：

1、目标位置求解准确性验证

本部分实验用于验证步骤一到步骤三建立的目标位置求解方法的正确性和在带有误差情况下的稳定性，以及讨论在不同条件下的可观测性。

1.1仿真环境设置

仿真时设置目标在距离观测平台300km处，与观测平台坐标系x轴成60°夹角，目标扫描周期设置为20s和10s，观测平台速度为550km/h并处于匀速直线运动状态，时间测量值精度为1μs，定位时采用公式(9)。

仿真时所用时间测量值可通过下式求出：

$arctan\left(\frac{y-{\mathrm{vt}}_i}{x}\right)=\frac{\mathrm{\π}}{3}-\frac{2\mathrm{\π}}{T}{t}_i+2\mathrm{\π}(i-1)---\left(28\right)$

其中i即表示第i次定位，i＝1,2,3...。

根据仿真设置可知目标距离很远，每次扫描到观测平台的时间与上次扫描时间的差值均近似为扫描周期T，因此可以合理简化求解过程，对t_i仅在t_i-1+T附近求解。

由于公式(28)中涉及到超越函数形式求解较为复杂，为简化时间测量值的获取可以采用如下方法：

$t=iT+arctan\left(\frac{x}{y-ivT}\right)\frac{T}{2\mathrm{\π}}-arctan\left(\frac{x}{y}\right)\frac{T}{2\mathrm{\π}}---\left(29\right)$

其中i表示探测次数，公式(29)利用目标每次扫描到观测平台的时间均在扫描周期倍数附近的规律，在每个周期平台所在位置点处求一次扫描时间，将其累加到总时间后减去初始时间即可得到时间测量值。由于求时间测量值时公式(28)右侧变量均为已知，故直接带入求解即可，运算量大幅减小。

然而该方法忽略了观测平台应当不断运动的情况，因此将带来额外的误差，对定位结果将带来一定影响。

综上所述，虽然公式(28)的精度受限但其计算量较小，适合仿真实验时采用，因此本部分中实验所需时间测量值均使用公式(28)所示方法得到。

1.2仿真结果分析

为便于比较分析，仿真结果图均在统一坐标尺度下绘制。

图7至图18是取扫描周期T＝20S时的仿真结果。

图7至图12是在获取时间测量值后直接带入求解所得到的定位结果。由图可知，虽然定位仍有抖动，初始定位的存在误差，但随着时间推移定位误差逐渐减小，解算结果趋近于真值，且结果的均值也接近真值，说明了定位公式的正确性。定位误差主要来源于精度受限的时间测量值，另一方面也和观测性有关，随着时间推移，两次扫描之间时间差增大使得观测性变好，从而有误差随时间减小的现象。

图13至图18是在时间测量值基础上叠加了方差为1us的高斯白噪声后再带入求解所得到的定位结果。从图中可以发现初始定位误差较大，而误差随时间增加而减小，这和未加入高斯白噪声时是相似的，最终定位结果趋于真值，说明定位公式的准确性。

对比两组结果可以发现，加入白噪声后定位误差增加明显，且收敛时间增长，这和利用扫描时间差定位方法本身有关。在该问题中有效信息蕴含在两次扫描时间差与扫描周期的差值中，目标距离远，时间差接近扫描周期，有效信息少，因此在此基础上的微量时间测量值误差也会导致较大的定位误差。

图19至图36是取扫描周期T＝10S时的仿真结果。

图19至图24是在获得时间测量值后直接带入求解所得到的定位结果。由图可知在初始之后相当长一段时间内定位误差都较大且有偏向，定位结果收敛较弱。在80～100点时定位点值有趋于平缓的趋势，但误差仍有20～40km，说明此时定位不准确。

图25至图30是在获得时间测量值后隔一点取一点带入公式求解所得到的结果。从图中可以发现，虽然初始定位有一定的误差但能以较快的速度收敛，最终收敛到真值附近，说明此时定位准确。相比于未间隔取点的结果，间隔取点结果的误差分布更为均匀，误差较小而且收敛快速，在60点左右已经降至10km以内，这是由于间隔取点后利用了连续两点所积累的信息，相当于扫描周期T＝20S的结果。因此，未采用间隔点方法时的定位不准，原因在于时间间隔中所含信息不足导致可观测性不好，而非定位公式问题。

将图25至图30与图7至图12相比较可以看到周期10S间隔取点定位的误差和收敛速度均逊于周期为20S时的结果，这是由于扫描周期短所造成的，因为周期短则会在更短的时间扫过相同的角度，使得扫描时间间隔稍小，降低了目标的可观测性。

图31至图36是进行隔两点取一点操作后进行定位得到的结果，此时利用的总信息为连续三点的积累信息。由图可以看出此时定位精度以及误差均优于隔一点取点的结果，和T＝20S相比定位也更为准确，说明间隔取点的确可以改善可观测性，也证明定位公式的正确。

图31至图42是在时间测量值基础上叠加了方差为1us的高斯白噪声后再隔一点取一点带入求解所得到的定位结果。可以看到在相当长时间内误差都较大但收敛明显，在80点左右误差已降至20km，趋势与未加噪声时相似。最终定位结果收敛在真值附近，说明定位准确，定位公式正确。

2、定位方法可行性实验

2.1仿真环境设置：

基本环境与1相同，设置目标在距离观测平台300km处，与观测平台坐标系x轴成60°夹角，目标扫描周期设置为20S和10S，观测平台速度为550km/h，时间测量值精度为1us。

时间观测值t_i采用公式(29)计算，P₀应用公式(17)与公式(18)计算，R_i用公式(15)和(16)计算，最终利用公式(10)和公式(11)进行递推运算。

2.2结果分析：

比较直接解算和采用递推最小二乘法处理后的结果可以发现，直接计算的结果收敛时间相对较长且收敛结果在真值附近仍有一定程度的抖动，而用递推最小二乘法处理后收敛加快且收敛结果相对更加稳定，只是结果相对真值有一较小恒定偏差，但偏差在可接受范围内。因此递推最小二乘法是优化根据扫描时间差无源定位结果的一种较为实用可行的方法。

对于收敛结果相对真值有一定偏差可以通过在递推最小二乘法中加入遗忘因子实现修正，但同时将会导致定位结果产生一定幅度抖动。因此实际中应权衡定位稳定性与准确性决定是否加入遗忘因子以及遗忘因子的取值。

Claims (5)

1.一种基于天线周期扫描时间差无源定位方法，其特征在于，所述定位方法包括以下步骤：

步骤一：建立坐标系并得到各坐标系间的转换关系；

建立观测平台坐标系，以观测平台运动方向为y轴方向，垂直y轴为x轴，y轴正方向右侧为x正方向；

建立目标坐标系，将观测平台的速度等效到目标上，采用等效后目标运动方向为y'轴方向，垂直y'轴为x'轴，y'轴正方向右侧为x'正方向，此时的坐标系相当于目标雷达对观测平台的定位坐标系，可以得到观测平台相对于目标的位置(x',y')，根据观测平台坐标系与目标坐标系的转换关系得到坐标变换x＝x'和y＝y'，得到目标相对于观测平台在平台坐标系中的位置(x,y)；

建立大地坐标系，以正北方位Y轴方向，正东方为X轴方向；当目标坐标需要转换到大地坐标下时，根据平台坐标系与大地坐标系的转换关系，即公式(1)-(3)，可得到大地坐标系下的目标坐标；

γ是观测平台坐标系y方向和大地坐标系Y方向的夹角；是目标在平台坐标系中的方位角，X_T为目标在大地坐标系中X方向位置，Y_T为目标在大地坐标系中Y方向位置，x₀为观测平台在大地坐标系中X方向位置，y₀为观测平台在大地坐标系中Y方向位置；

步骤二：根据步骤一所述的平台坐标系与目标坐标系的转换关系以及大地坐标系与平台坐标系的转换关系，可以将目标坐标系中所求观测平台坐标转换为大地坐标系中目标坐标，直接在目标坐标系中对观测平台位置进行求解，等效后目标进行匀速直线运动时，可以得到求解方程中x,y,T三个未知量需要的方程为：

$\left\{\begin{array}{c}2\mathrm{\π}({\mathrm{\τ}}_i-T)\[x^2+y^2-(y_{i+1}+y_i)y+y_{i+1}{y}_i\]-xT(y_{i+1}-y_i)=0\\ 2\mathrm{\π}({\mathrm{\τ}}_{i+1}-T)\[x^2+y^2-(y_{i+2}+y_{i+1})y+y_{i+2}{y}_{i+1}\]-xT(y_{i+2}-y_{i+1})=0\\ 2\mathrm{\π}({\mathrm{\τ}}_{i+2}-T)\[x^2+y^2-(y_{i+3}+y_{i+2})y+y_{i+3}{y}_{i+2}\]-xT(y_{i+3}-y_{i+2})=0\end{array}\right.---\left(9\right)$

所述T为天线扫描周期，τ_i＝t_i+1-t_i，t_i为第i时刻天线扫描时间测量值，t_i+1为第i+1时刻天线扫描时间测量值，y_i为目标坐标系中第i时刻平台的y'方向坐标，y_i+1为目标坐标系中第i+1时刻平台的y'方向坐标，y_i+2为目标坐标系中第i+2时刻平台的y'方向坐标，y_i+3为目标坐标系中第i+3时刻平台的y'方向坐标；

步骤三：求解公式(9)的方程组得到6组数学解，解的选取利用探测范围的限制，所述探测范围为100～350km，在探测范围内的解记为可行解，进一步筛选可行解集合时，若有唯一解则将其作为解，若无可行解则将解计为零，有多组解可以取各解的平均值；

步骤四：求解递推最小二乘法的参数R_i、和P₀，将遗忘因子与R_i相乘得到R_i'，将R_i'、和P₀带入递推最小二乘的基本公式，利用通过时间测量值计算得到的z_i估计目标的真实位置得到估计值所述z_i为第i时刻的观测向量(x_i,y_i,T_i)^T，为第i时刻的估计值矩阵T_i为第i时刻天线扫描周期计算值，为第i时刻天线扫描周期估计值，为目标在x'方向坐标的估计值，为目标在y'方向坐标的估计值。

2.根据权利要求1所述的一种基于天线周期扫描时间差无源定位方法，其特征在于所述步骤二中得到求解方程中x,y,T三个未知量需要的三个方程的具体过程为：

θ_i的表达式为：

${\mathrm{\θ}}_i=arctan\left(\frac{x}{y-y_{i+1}}\right)-arctan\left(\frac{x}{y-y_i}\right)---\left(4\right)$

所述θ_i为第i时刻与第i+1时刻目标与观测平台连线间的夹角；

根据两角和与差的正切公式可以得到：

$tan\left({\mathrm{\θ}}_i\right)=\frac{x(y_{i+1}-y_i)}{x^2+y^2-(y_{i+1}+y_i)y+y_i{y}_{i+1}}---\left(5\right)$

目标距离为100～350km，θ_i满足的条件，因此取：

tan(θ_i)≈θ_i(6)

目标雷达天线以固定周期扫描，可以得到公式(7)：

${\mathrm{\τ}}_i=t_{i+1}-t_i=T+T\frac{{\mathrm{\θ}}_i}{2\mathrm{\π}}---\left(7\right)$

联立公式(6)和公式(7)，消去θ_i后可以得到方程(8)：

2π(τ_i-T)[x²+y²-(y_i+1+y_i)y+y_i+1y_i]-xT(y_i+1-y_i)＝0(8)

求解方程中x,y,T三个未知量需要三个方程，利用连续四个点，每相邻两点为一组，分别代入公式(8)，得到三个方程，组成方程组为公式(9)。

3.根据权利要求2所述的一种基于天线周期扫描时间差无源定位方法，其特征在于所述步骤四中求解递推最小二乘法的参数R_i、和P₀的具体过程为：

R_i的求解过程为：

对公式(8)两侧求时间t的偏导，可得：

$\begin{array}{c}\frac{dx}{dt}\[4\mathrm{\π}x({\mathrm{\τ}}_i-T)-T(y_{i+1}+y_i)\]+\frac{dy}{dt}\[4\mathrm{\π}y({\mathrm{\τ}}_i-T)-2\mathrm{\π}({\mathrm{\τ}}_i-T)(y_{i+1}+y_i)\]+\\ \frac{dT}{dt}\{-2\mathrm{\π}\[x^2+y^2-(y_{i+1}+y_i)y+y_{i+1}{y}_i\]-x(y_{i+1}-y_i)\}=2\mathrm{\π}({\mathrm{\τ}}_i-T)(2v_{y^{\′}i}y-v_{y^{\′}i}{y}_{i+1}-v_{y^{\′}i}{y}_i)\end{array}---\left(12\right)$

取连续四点，每相邻两点为一组，分别代入公式(12)，得到三个方程，改写为矩阵形式如下：

$\begin{array}{c}\left[\begin{array}{ccc}4\mathrm{\π}x({\mathrm{\τ}}_i-T)-T(y_{i+1}-y_i)& 4\mathrm{\π}y({\mathrm{\τ}}_i-T)-2\mathrm{\π}({\mathrm{\τ}}_i-T)(y_{i+1}+y_i)& -2\mathrm{\π}\[x^2+y^2-(y_{i+1}+y_i)y+y_{i+1}{y}_i\]-x(y_{i+1}-y_i)\\ 4\mathrm{\π}x({\mathrm{\τ}}_{i+1}-T)-T(y_{i+2}-y_{i+1})& 4\mathrm{\π}y({\mathrm{\τ}}_{i+1}-T)-2\mathrm{\π}({\mathrm{\τ}}_{i+1}-T)(y_{i+2}+y_{i+1})& -2\mathrm{\π}\[x^2+y^2-(y_{i+2}+y_{i+1})y+y_{i+2}{y}_{i+1}\]-x(y_{i+2}-y_{i+1})\\ 4\mathrm{\π}x({\mathrm{\τ}}_{i+2}-T)-T(y_{i+3}-y_{i+2})& 4\mathrm{\π}y({\mathrm{\τ}}_{i+2}-T)-2\mathrm{\π}({\mathrm{\τ}}_{i+2}-T)(y_{i+3}+y_{i+2})& -2\mathrm{\π}\[x^2+y^2-(y_{i+3}+y_{i+2})y+y_{i+3}{y}_{i+2}\]-x(y_{i+3}-y_{i+2})\end{array}\right]\\ \left[\begin{array}{c}\frac{dx}{dt}\\ \frac{dy}{dt}\\ \frac{dT}{dt}\end{array}\right]=\left[\begin{array}{c}2\mathrm{\π}({\mathrm{\τ}}_i-T)(2v_{y^{\′}i}y-v_{y^{\′}i}{y}_{i+1}-v_{y^{\′}i}{y}_i)\\ 2\mathrm{\π}({\mathrm{\τ}}_{i+1}-T)(2v_{y^{\′}i}y-v_{y^{\′}i}{y}_{i+2}-v_{y^{\′}i}{y}_{i+1})\\ 2\mathrm{\π}({\mathrm{\τ}}_{i+2}-T)(2v_{y^{\′}i}y-v_{y^{\′}i}{y}_{i+3}-v_{y^{\′}i}{y}_{i+2})\end{array}\right]\end{array}---\left(13\right)$

公式(13)可表示为A_3×3f_3×1＝b_3×1，则有：

$f_{3\×1}={\left(A_{3\×3}^T{A}_{3\×3}\right)}^{-1}{A}_{3\×3}^T{b}_{3\×1}---\left(14\right)$

根据观测误差协方差定义可知，R_i计算式为：

$R_i=\left[\begin{array}{ccc}{\left(\frac{dx}{dt}\right)}^2{\mathrm{\σ}}_t^2& 0& 0\\ 0& {\left(\frac{dy}{dt}\right)}^2{\mathrm{\σ}}_t^2& 0\\ 0& 0& {\left(\frac{dT}{dt}\right)}^2{\mathrm{\σ}}_t^2\end{array}\right]---\left(15\right)$

由公式(14)和公式(15)可得到R_i，所述R_i表示第i时刻的观测误差协方差，为时间测量值的误差的方差；

和P₀的求解过程为：

$P_0={\left({\stackrel{\‾}{H}}^T{\stackrel{\‾}{R}}^{-1}\stackrel{\‾}{H}\right)}^{-1}---\left(17\right)$

${\hat{X}}_0={\left({\stackrel{\‾}{H}}^T{\stackrel{\‾}{R}}^{-1}\stackrel{\‾}{H}\right)}^{-1}{\stackrel{\‾}{H}}^T{\stackrel{\‾}{R}}^{-1}\stackrel{\‾}{Z}---\left(18\right)$

其中和由起始天线扫描时间测量数据和平台状态数据，所得的相应矩阵作为分块矩阵组合而成，即：

$\stackrel{\‾}{H}={\[H_{3\×3},H_{3\×3},...H_{3\×3}\]}^T---\left(19\right)$

$\stackrel{\‾}{Z}={\left[\begin{array}{cccc}{z}_0^T& z_1^T& ...& z_n^T\end{array}\right]}^T---\left(21\right)$

所述H_3×3表示三阶单位阵；

利用公式(17)-(21)可求得最小二乘递推初始值，即和P₀。

4.根据权利要求3所述的一种基于天线周期扫描时间差定位方法，其特征在于所述步骤四中得到R_i'的具体过程为：

R_i'＝e^(l-i)αR_i(16)

其中所述l是选择用于位置估计的R_i阵总个数，α是为常数，R_i为公式(16)计算所得结果，进行递推时将R_i'带入公式(10)和公式(11)中的R_i即可。

5.根据权利要求4所述的一种基于天线周期扫描时间差无源定位方法，其特征在于所述步骤四中递推最小二乘的基本公式为：

${\hat{X}}_{i+1}={\hat{X}}_i+P_i{H}_{i+1}^T{\[H_{i+1}{P}_i{H}_{i+1}^T+R_{i+1}\]}^{-1}\[z_{i+1}-H_{i+1}{X}_i\]---\left(10\right)$

$P_{i+1}=P_i-P_i{H}_{i+1}^T{\[H_{i+1}{P}_i{H}_{i+1}^T+R_{i+1}\]}^{-1}{H}_{i+1}{P}_i---\left(11\right)$

其中表示第i时刻的估计值矩阵z_i表示第i时刻的观测向量(x_i,y_i,T_i)^T，P_i表示第i时刻的协方差矩阵，H_i表示第i时刻的测量矩阵在此处即是3×3的单位阵。