CN102540180A - 天基相控阵雷达空间多目标定轨方法 - Google Patents

Abstract

本发明公开了一种天基相控阵雷达空间多目标定轨方法，主要解决现有技术对空间弱目标不能有效估计和目标定轨精度低的问题。其过程为：对目标回波数据采用调零保形算法进行处理，运用距离脉冲压缩原理得到距离先验信息，并据此对回波信号进行分段处理；根据目标数将回波信号分段，每段数据为目标信号及其临近单元信号；对各段数据运用和差多波束测角原理实现多目标二维角度估计；基于空间目标跟踪结果和探测星轨道信息，进行坐标转换；采用Laplace型迭代算法实现对空间远近不同目标定轨，并运用最小二乘算法提高定轨精度。本发明能够降低强信号对弱目标的影响，对弱目标参数进行准确估计，可用于空间态势感知、轨道资源管理等实际应用领域。

Description

天基相控阵雷达空间多目标定轨方法

技术领域

本发明涉及信号处理技术领域，涉及卫星定轨，特别是基于天基相控阵雷达对空间远近不同目标定轨方法，可用于空间多目标的检测、跟踪和定轨。

背景技术

对空间目标进行定轨，所获取的目标轨道信息能够为轨道数据库建立提供数据支持，也可以为空间目标识别提供先验信息，对于空间态势感知、轨道资源管理和保障在轨飞行器安全具有重要意义。

美国目前的空间监视网(SSN：Space Surveillance Network)主要由3类系统组成：(1)地基雷达监视系统，主要包括：FPS-85雷达、铺路爪(Pave Paws)雷达、海军空间监视系统的雷达等；(2)地基光电系统，主要由“地基光电深空监视系统”(GEODSS：ground-based electro-optical deep space surveillance)组成；(3)天基监视系统，包括2011年春季投入运行的“天基空间监视系统”(SBSS：Space-based SpaceSurveillance System)、计划2015年投入运行的“轨道深空成像仪”(ODSI：Orbit DeepSpace Imager)和其他天基计划。SBSS部署在低地球轨道(1100Km)，携带高灵敏度的天基可视传感器(SBV：Space Based Visible)和1部可探测地球静止轨道(36000Km)卫星的光学望远镜，其目的是研制和部署能用于深空和近地空间目标的发现、跟踪和定轨；ODSI系统由运行在地球静止轨道的成像卫星组成的卫星星座组成，其主要功能是提供地球静止轨道上三轴稳定卫星的图像，即获取空间目标的特征信息，进行空间目标识别；其他天基计划如为导弹防御系统研制的“天基红外系统”(SBIRS：SpaceBased Infrared System)和“空间跟踪与监视系统(STSS：Space Track and SurveillanceSystem)”也具有很强的天基空间监视能力。俄罗斯的空间目标监视系统(SSS：SpaceSurveillance System)由用于监视低空目标的地基雷达和探测中低轨目标的光电望远镜组成，如天窗0系统。除美国外，加拿大、英国、日本、德国等国家也开展了天基空间监视的相关研究：加拿大2003年发射了空间探测小卫星，在此基础上启动了空间目标监视(SOS：Surveillance of Space)计划，主要用于跟踪地球同步轨道上通信卫星和其他高轨卫星；英国SURRY大学于2000年发射了实验小卫星，获得了目标空间光学图像；日本于2003年曾发射了空间垃圾回收卫星，具备天基空间目标监视能力；德国波恩大学开展了天基空间目标监视有效载荷技术研究。我国目前的空间目标监视系统主要由地基系统组成，对空间目标的探测和定轨能力有限。可见，目前天基监视系统国内外主要以光电系统为主，天基相控阵雷达系统鲜见报导。

地基系统相比于天基系统虽具有技术成熟、成本低、好维护等优点，但受到诸多条件的限制：(1)观测时间受限：地基系统受地域限制严重，难以实现对天基目标的连续监视与跟踪；(2)探测距离有限：地基系统无论是光电或是雷达，探测目标经过大气衰减导致信噪比降低从而影响远距离高轨目标的探测。雷达相比光电系统，具有如下的优点：(1)全天时、全天候工作，雷达可24小时不间断工作，而光电系统只能在目标被照亮的情况下才能工作；(2)雷达能够快速检测、跟踪和定轨目标，由于主动探测可提供目标多维信息：方位角、俯仰角、距离和速度等，能够实时对目标定轨，而光学传感器是被动仅测角形式，只能测量目标的赤经和赤纬，无先验信息的条件下往往需要数小时或数天才能对空间目标定轨。

将相控阵雷达安置于天基平台，除具有波束控制灵活的优点外，还具备空间多目标同时观测的优点，因此，研究天基相控阵雷达空间多目标定轨方法具有十分重要的意义。基于喇叭天线的传统和差测角虽具有实现简单、测向精度高的优点，但只能针对单个目标；基于子空间的超分辨测向方法，如多重信号分类(MUSIC：Multiple SignalClasssification)法、信号参数旋转不变估计(ESPRIT：Estimation of signal parametersvia rotational invariance techniques)法等，虽能够实现同时多目标测向，但实际应用中存在运算量大的难题。

发明内容

针对上述问题，本发明提供了一种基于天基相控阵雷达的空间多目标定轨方法，以克服空间多目标的远近效应，解决远距离目标由于信噪比低而导致的测角精度差的问题，提高轨道估计精度。

实现本发明目的的技术思路是：采用调零保形对多目标进行空分，通过距离脉压估计目标距离信息，由距离分段克服“远近”目标效应；由和差多波束测角估计多空间目标的方位和俯仰信息；根据空间目标距离和角度信息，采用Laplace型迭代的最小二乘方法实现对空间目标的定轨。其具体步骤包括如下：

(1)将相控阵雷达布置在卫星平台上，按脉冲体制发射信号s_e(t)，并接收回波数据，运用展宽调零保形算法，计算第k个目标的最优波束权值w_k，k＝1，...，P，P为目标数；

(2)用第k个目标的最优波束权值w_k乘以回波数据得到波束输出数据y_k，并对该波束输出数据y_k进行脉冲压缩处理，得到第k个目标的距离时延τ_k，k＝1，...，P，P为目标数；

(3)按阵元位置将相控阵分为左上、左下、右上和右下四个象限，把第k个目标的普通波束权值w_{q_k}划分为左上、左下、右上和右下权值，分别为

和

用

和

分别乘以对应象限的回波数据后，得到左上、左下、右上和右下象限t时刻的波束输出分别为

和

由左上、左下、右上和右下的波束输出构造和波束输出

方位差波束输出和俯仰差波束输出

$y_{q\_k}^{\mathrm{\Σ}}\left(t\right)=y_{q\_k}^{\mathrm{LU}}\left(t\right)+y_{q\_k}^{\mathrm{LD}}\left(t\right)+y_{q\_k}^{\mathrm{RU}}\left(t\right)+y_{q\_k}^{\mathrm{RD}}\left(t\right),$

$y_{q\_k}^X\left(t\right)=y_{q\_k}^{\mathrm{RU}}\left(t\right)+y_{q\_k}^{\mathrm{RD}}\left(t\right)-y_{q\_k}^{\mathrm{LU}}\left(t\right)-y_{q\_k}^{\mathrm{LD}}\left(t\right),$

$y_{q\_k}^Y\left(t\right)=y_{q\_k}^{\mathrm{RU}}\left(t\right)+y_{q\_k}^{\mathrm{LU}}\left(t\right)-y_{q\_k}^{\mathrm{RD}}\left(t\right)-y_{q\_k}^{\mathrm{LD}}\left(t\right);$

(4)以步骤(2)中的第k个目标的距离时延τ_k为中心，从和波束输出

方位差波束输出

和俯仰差波束输出

中均截取2T_P范围数据，T_P为发射脉冲宽度，对所截取的数据分别作距离脉压，得到和距离脉压数据

方位差距离脉压数据

俯仰差距离脉压数据

$z_{q\_k}^{\mathrm{\Σ}}\left(\mathrm{\τ}\right)={\∫}_{{\mathrm{\τ}}_k-T_p}^{{\mathrm{\τ}}_k+T_p}{y}_{q\_k}^{\mathrm{\Σ}}\left(t\right)s_e^{*}(\mathrm{\τ}-t)\mathrm{dt},$

$z_{q\_k}^X\left(\mathrm{\τ}\right)={\∫}_{{\mathrm{\τ}}_k-T_p}^{{\mathrm{\τ}}_k+T_p}{y}_{q\_k}^X\left(t\right)s_e^{*}(\mathrm{\τ}-t)\mathrm{dt},$

$z_{q\_k}^Y\left(\mathrm{\τ}\right)={\∫}_{{\mathrm{\τ}}_k-T_p}^{{\mathrm{\τ}}_k+T_p}{y}_{q\_k}^Y\left(t\right)s_e^{*}(\mathrm{\τ}-t)\mathrm{dt},$

其中*表示共轭，k＝1，...，P，P为目标数，在

中寻找到和距离脉压峰值

max(·)表示求最大值，同样可获得方位差距离脉压峰值

和俯仰差距离脉压峰值

(5)利用上述三个波束输出的距离脉压峰值得到俯仰偏角估计

和方位偏角估计

${\mathrm{\ϵ}}_k^{\mathrm{\θ}}=\arctan \left(\mathrm{imag}\left(\frac{P_{q\_k}^Y}{P_{q\_k}^{\mathrm{\Σ}}}\right)\right)\×\frac{\mathrm{\λ}}{2\mathrm{\π}{L}_y\cos {\mathrm{\θ}}_k^B}$

其中arctan(·)表示反正切函数，imag(·)表示取复数虚部操作，

为俯仰差比和，λ为波长，L_x为X方向的等效相位中心坐标，L_y为Y方向的等效相位中心坐标，为波束指向的俯仰角，为波束指向的方位角，为方位差比和，k＝1，...，P，P为目标数；

(6)由俯仰和方位偏角得到目标俯仰角和方位角估计值：

k＝1，...，P，P为目标数；

(7)基于坐标转换原理，利用对空间目标角度和距离的估计，计算空间目标在J2000坐标系下的坐标；

(8)基于Laplace迭代算法，运用最小二乘原理估计空间目标的轨道参数。

本发明与现有技术相比，具有以下优点：

(1)本发明通过对最优波束权值的计算，克服了空间多目标的远近效应，可同时对空间多目标进行检测。

(2)本发明对接收数据进行分段处理，抑制了强目标信号对弱目标的影响，降低了弱目标参数估计的误差。

(3)本发明在传统定轨方法的基础上，运用最小二乘原理估计空间目标的轨道参数，提高了定轨精度。

附图说明

图1是本发明的流程图；

图2是本发明的相控阵雷达工作场景图；

图3是本发明的二维均匀圆阵结构图；

图4是用本发明对目标1和目标2的距离进行估计的结果图；

图5是用本发明对目标1和目标2的距离进行估计的误差图。

图6是用本发明对目标1方位角和俯仰角进行估计的结果图；

图7是用本发明对目标1方位角和俯仰角进行估计的误差图；

图8是用本发明对目标2方位角和俯仰角进行估计的结果图；

图9是用本发明对目标2方位角和俯仰角进行估计的误差图；

具体实施方式

本发明使用的工作场景图如图2所示，相控阵雷达以探测星为平台，将阵列布置于Y′O′Z′平面，其结构如图3所示。

根据图3建立信号模型，得到t时刻各阵元回波数据x(t)为：

x(t)＝As(t)+n(t)，    1)

其中P为目标数，θ_k和

分别为第k个目标的俯仰角和方位角，k＝1，...，P，

为第k个目标的导向矢量，N表示阵元数，T表示转置，s(t)＝[s₁(t)  s₂(t)…s_P(t)]^T为信号包络，n(t)＝[n₁(t)，n₂(t)，...，n_N(t)]^T为噪声矢量。

参照图1，本发明基于上述信号模型，即t时刻各阵元回波数据x(t)，所提出的一种天基相控阵雷达对空间远近不同目标定轨的方法，具体步骤如下：

步骤1：对各阵元回波数据x(t)运用调零保形算法，计算第k个目标的最优波束权值w_k，k＝1，...，P，P为目标数。

调零保形算法可描述为一个带约束的二次优化问题：

$\left\{\begin{array}{cc}\min & {\left|\right|w_k-w_{q\_k}\left|\right|}^2\\ s.t.& C_k^H{w}_k=f\end{array}\right.---2)$

其中

表示第k个目标的普通波束权值，h表示加窗函数，

为波束指向的导向矢量，

为第k个波束指向，f＝[0 0 0]^T表示约束向量，||·||²表示矢量的2范数，

表示第k个约束矩阵，H表示共轭转置，F＝[I_N D_x D_y]表示系数矩阵，I_N表示N维单位矩阵，D_x＝diag(x₁ x₂...x_N)表示由阵元位置x方向坐标(x₁ x₂ ... x_N)构成的对角矩阵，D_y＝diag(y₁ y₂ ... y_N)表示由阵元位置y方向坐标(y₁ y₂ ... y_N)构成的对角矩阵，D_x和D_y用于展宽零陷以提高算法的稳健性；

对上述二次优化问题的解即计算第k个目标的最优波束权值w_k，通过构造拉格朗日代价函数得到：

$w_k=w_{q\_k}-C_k{\left(C_k^H{C}_k\right)}^{-1}{C}_k^H{w}_{q\_k}-C_k{\left(C_k^H{C}_k\right)}^{-1}f.---3)$

步骤2：用第k个目标的最优波束权值w_k乘以回波数据得到波束输出数据y_k，并对该波束输出数据y_k进行脉冲压缩处理，得到第k个目标的距离时延τ_k，k＝1，...，P；

步骤3：按阵元位置将相控阵分为左上、左下、右上和右下四个象限，把第k个目标的普通波束权值w_{q_k}根据四个象限划分为左上、左下、右上和右下权值，左上权值为

左下权值为

右上权值为

右下权值为

k＝1，...，P。

步骤4：用

和

分别乘以对应象限的回波数据，得到左上、左下、右上和右下象限t时刻的波束输出，左上象限波束输出为

左下象限波束输出为

右上象限波束输出为右下象限波束输出为由这四个象限波束输出按对称取反方式构造和波束输出：

$y_{q\_k}^{\mathrm{\Σ}}\left(t\right)=y_{q\_k}^{\mathrm{LU}}\left(t\right)+y_{q\_k}^{\mathrm{LD}}\left(t\right)+y_{q\_k}^{\mathrm{RU}}\left(t\right)+y_{q\_k}^{\mathrm{RD}}\left(t\right),$ 方位差波束输出：

$y_{q\_k}^X\left(t\right)=y_{q\_k}^{\mathrm{RU}}\left(t\right)+y_{q\_k}^{\mathrm{RD}}\left(t\right)-y_{q\_k}^{\mathrm{LU}}\left(t\right)-y_{q\_k}^{\mathrm{LD}}\left(t\right)$ 和俯仰差波束输出：

$y_{q\_k}^Y\left(t\right)=y_{q\_k}^{\mathrm{RU}}\left(t\right)+y_{q\_k}^{\mathrm{LU}}\left(t\right)-y_{q\_k}^{\mathrm{RD}}\left(t\right)-y_{q\_k}^{\mathrm{LD}}\left(t\right),$ k＝1，...，P。

步骤5：以步骤2中的第k个目标的距离时延τ_k为中心，从和波束输出

方位差波束输出

和俯仰差波束输出

中均截取2T_P范围数据，T_P为发射脉冲宽度，对所截取数据分别作距离脉压，得到和距离脉压数据

方位差距离脉压数据

俯仰差距离脉压数据

$z_{q\_k}^{\mathrm{\Σ}}\left(\mathrm{\τ}\right)={\∫}_{{\mathrm{\τ}}_k-T_p}^{{\mathrm{\τ}}_k+T_p}{y}_{q\_k}^{\mathrm{\Σ}}\left(t\right)s_e^{*}(\mathrm{\τ}-t)\mathrm{dt},---4)$

$z_{q\_k}^X\left(\mathrm{\τ}\right)={\∫}_{{\mathrm{\τ}}_k-T_p}^{{\mathrm{\τ}}_k+T_p}{y}_{q\_k}^X\left(t\right)s_e^{*}(\mathrm{\τ}-t)\mathrm{dt},---5)$

$z_{q\_k}^Y\left(\mathrm{\τ}\right)={\∫}_{{\mathrm{\τ}}_k-T_p}^{{\mathrm{\τ}}_k+T_p}{y}_{q\_k}^Y\left(t\right)s_e^{*}(\mathrm{\τ}-t)\mathrm{dt},---6)$

其中*表示共轭，k＝1，...，P，在中寻找到脉压峰值

其中max(·)表示求最大值，同样可获得和

的脉压峰值分别记为

和

其中 $P_{q\_k}^X=\underset{\mathrm{\τ}}{\max }\left(z_{q\_k}^X\left(\mathrm{\τ}\right)\right),$ $P_{q\_k}^Y=\underset{\mathrm{\τ}}{\max }\left(z_{q\_k}^Y\left(\mathrm{\τ}\right)\right).$

步骤6：利用上述三个波束输出的距离脉压峰值得到俯仰偏角估计

和方位偏角估计

${\mathrm{\ϵ}}_k^{\mathrm{\θ}}=\arctan \left(\mathrm{imag}\left(\frac{P_{q\_k}^Y}{P_{q\_k}^{\mathrm{\Σ}}}\right)\right)\×\frac{\mathrm{\λ}}{2\mathrm{\π}{L}_y\cos {\mathrm{\θ}}_k^B}---7)$

其中arctan(·)表示反正切函数，imag(·)表示取复数虚部操作，为俯仰差比和，λ为波长，L_x为X方向的等效相位中心坐标，L_y为Y方向的等效相位中心坐标，

为波束指向的俯仰角，

为波束指向的方位角，为方位差比和，k＝1，...，P，P为目标数；

步骤7：由俯仰和方位偏角得到第k个目标的俯仰角和方位角估计值： ${\widehat{\mathrm{\θ}}}_k={\mathrm{\θ}}_k^B+{\mathrm{\ϵ}}_k^{\mathrm{\θ}},$ k＝1，...，P。

步骤8：基于坐标转换原理，利用对空间目标距离和角度的估计，计算空间目标在J2000坐标系下的坐标r_k。

8a)计算J2000坐标系下第k个目标相对于雷达的向量ρ_k，k＝1，...，P，即对目标距离和角度的估计是以雷达为参考的，在相控阵雷达工作场景图中，参照图2所示，XYZ坐标系是J2000坐标系，X′Y′Z′坐标系是以雷达为参考的坐标系，基于坐标转换原理，ρ_k可按如下公式计算得到：

ρ_k＝C[x′_T y′_T z′_T]^T    9)

其中C为坐标转换矩阵，可由X′Y′Z′坐标系和XYZ坐标系的基矢量之间的关系得到，

[x′_kT y′_kT z′_kT]^T为第k个目标在X′Y′Z′坐标系中的位置矢量，T表示转置，

为第k个目标的距离，其中c为光速，ρ_k也可写成如下形式：ρ_k＝|ρ_k|[α_k β_k γ_k]^T，|ρ_k|为矢量ρ_k长度，[α_k β_k γ_k]^T为矢量ρ_k的单位矢量。

8b)由图2相控阵雷达工作场景图中的探测星、空间目标以及地球之间的几何关系，按如下公式计算第k个空间目标在J2000坐标系下的坐标r_k：

r_k＝ρ_k+R    10)

其中R＝[x_R y_R z_R]^T表示已知的探测星在J2000坐标系下的坐标，x_R、y_R、z_R分别为坐标的三个分量。

步骤9：基于Laplace迭代算法，运用最小二乘原理估计空间目标的轨道参数。

9a)基于Laplace迭代算法进行迭代，计算第k个目标t₀时刻的速度矢量v_k0的粗略估值

k＝1，...，P，即将第k个目标在J2000坐标系下的坐标r_k在t₀时刻进行泰勒展开得：

$r_k=r_{k0}+v_{k0}\mathrm{\τ}+\frac{1}{2}{a}_{k0}{\mathrm{\τ}}^2+\frac{1}{3!}\frac{d^3{r}_k}{{\mathrm{dt}}^3}{|}_{t=t_0}{\mathrm{\τ}}^3+O\left({\mathrm{\τ}}^4\right),---11)$

其中τ＝t-t₀表示时间间隔，r_k0＝[x_k0 y_k0 z_k0]^T表示t₀时刻第k个目标在J2000坐标系下的坐标，x_k0、y_k0、z_k0为r_k0的三个分量，

为t₀时刻第k个目标在J2000坐标系下的速度矢量，

为t₀时刻第k个目标在J2000坐标系下的加速度矢量。

由11)式经推导，r_k可写成

r_k＝Fr_k0+Gv_k0，    12)

其中 $F=1-\frac{{\mathrm{\μ\τ}}^2}{2{\left|r_{k0}\right|}^3}+\frac{{\mathrm{\μ\τ}}^3(r_{k0}\·v_{k0})}{2{\left|r_{k0}\right|}^5}+O\left({\mathrm{\τ}}^4\right)$ 为r_k0的系数， $G=\mathrm{\τ}-\frac{{\mathrm{\μ\τ}}^3}{6{\left|r_{k0}\right|}^3}+O\left({\mathrm{\τ}}^4\right)$ 为v_k0的系数，μ表示中心天体的质心引力常数，|r_k0|为t₀时刻空间目标的距离。

由10)式和12)式知，

可按

进行迭代计算，迭代时令初值F⁽⁰⁾＝1，G⁽⁰⁾＝τ。

9b)利用最小二乘原理，计算精确的第k个目标在J2000坐标系下的坐标r′_k0和速度矢量v′_k0，k＝1，...，P。

由图2相控阵雷达工作场景图中的探测星、空间目标以及地球之间的几何关系式：

$\left[\begin{array}{cc}{r}_{k0}^{\′}& v_{k0}^{\′}\end{array}\right]\left[\begin{array}{cccc}{F}^{\left({\mathrm{\τ}}_1\right)}& F^{\left({\mathrm{\τ}}_2\right)}& \·\·\·& F^{\left({\mathrm{\τ}}_M\right)}\\ G^{\left({\mathrm{\τ}}_1\right)}& G^{\left({\mathrm{\τ}}_2\right)}& \·\·\·& G^{\left({\mathrm{\τ}}_M\right)}\end{array}\right]=\left[\begin{array}{cccc}{\mathrm{\ρ}}_k\left({\mathrm{\τ}}_1\right)+R\left({\mathrm{\τ}}_1\right)& {\mathrm{\ρ}}_k\left({\mathrm{\τ}}_2\right)+R\left({\mathrm{\τ}}_2\right)& \·\·\·& {\mathrm{\ρ}}_k\left({\mathrm{\τ}}_M\right)+R\left({\mathrm{\τ}}_M\right)\end{array}\right]$

得到矩阵[r′_k0 v′_k0]＝DB^H(BB^H)^-1，这里矩阵DB^H(BB^H)^-1的第一列为r′_k0，矩阵DB^H(BB^H)^-1的第二列为v′_k0，

其中D＝[ρ_k(τ₁)+R(τ₁)ρ_k(τ₂)+R(τ₂)...ρ_k(τ_M)+R(τ_M)]表示空间目标的位置矢量矩阵，

$B=\left[\begin{array}{cccc}{F}^{\left({\mathrm{\τ}}_1\right)}& F^{\left({\mathrm{\τ}}_2\right)}& \·\·\·& F^{\left({\mathrm{\τ}}_M\right)}\\ G^{\left({\mathrm{\τ}}_1\right)}& G^{\left({\mathrm{\τ}}_2\right)}& \·\·\·& G^{\left({\mathrm{\τ}}_M\right)}\end{array}\right]$ 表示系数矩阵，

表示时间间隔为τ_i时的F值，

表示时间间隔为τ_i时的G值，i＝1，2，...，M，其中M为雷达测量次数，

和

可由9a)中的

目标角度距离信息及时间间隔τ_i计算得到，ρ_k(τ_i)表示时间间隔为τ_i时的ρ_k值，R(τ_i)表示时间间隔为τ_i时的R值，i＝1，2，...，M。

9c)根据步骤9b)得到的精确的第k个目标在J2000坐标系下的坐标r′_k0和速度矢量v′_k0，k＝1，...，P，计算卫星轨道参数，即轨道半长轴a、轨道偏心率e、轨道倾角i、轨道升交点赤经Ω、轨道真近点角f、轨道近地点幅角w：

轨道半长轴：

式中，|r′_k0|为r′_k0的模值，|v′_k0|²＝v′_k0 ^Tv′_k0，T表示转置，G_g为万有引力常量，M为地球质量，G_gM＝3.986005×10¹⁴；

轨道偏心率：

其中，|h|为矢量h的模值，h＝r′_k0×v′_k0＝[h_x h_y h_z]^T，h_x、h_y和h_z分别为h的三个分量；

轨道倾角：

arccos(·)表示反余弦函数；

轨道升交点赤经： $\mathrm{\Ω}=\arctan (-\frac{h_x}{h_y});$

轨道真近点角： $f=\arccos \left\{\frac{1}{e}\right[\frac{a(1-e^2)}{\left|r_{k0}^{\′}\right|}-1\left]\right\};$

轨道近地点幅角：

arcsin(·)表示反正弦函数，z_k0为矢量r′_k0的Z坐标分量。

本发明的效果可通过以下仿真进一步说明：

仿真参数：选平面阵为36元同心圆阵，目标1和目标2的信噪比分别为10dB，30dB。

仿真内容及结果：

仿真1，采用本发明对目标1和目标2的距离进行估计，仿真结果如图4和图5所示，其中图4为用本发明对目标1和目标2的距离进行估计的结果图，图5为用本发明对目标1和目标2的距离进行估计的误差图，由图4和图5可见，两目标的测距误差小于一个距离门宽度，测距准确、可靠。

仿真2，采用本发明对目标1方位角和俯仰角进行估计，仿真结果如图6和图7，其中图6为用本发明对目标1方位角和俯仰角进行估计的结果图，图7为用本发明对目标1方位角和俯仰角进行估计的误差图，由图6和图7可见，目标1的角度估计值与真实值非常接近，角度估计的误差范围为-0.1°～0.08°，误差较小。

仿真3，采用本发明对目标2方位角和俯仰角进行估计，仿真结果如图8和图9，其中图8为用本发明对目标2方位角和俯仰角进行估计的结果图，图9为用本发明对目标2方位角和俯仰角进行估计的误差图，由图8和图9可见，目标2角度估计的误差范围为-0.012°～0.007°，测角效果良好。

从图6～9还可以看出，目标1的信噪比很低，受到目标2强信号的影响，测角精度比目标2低，但目标1的角度仍然可以有效估计出来。

仿真4，本发明利用角度和距离信息对目标1采用传统算法与最小二乘算法进行定轨，定轨结果如表1所示：

表1对目标1采用传统算法与最小二乘算法定轨结果

轨道参数	设定值	传统方法定轨	最小二乘算法定轨
				i	55°	54.8105°	54.9994°
Ω	81°	80.7071°	80.9990°
				ω	60°	60.5343°	60.0023°
a	7500km	7504.5km	7499.3km
				e	0.05	0.0507	0.0499
f	-9.4728°	-9.8306°	-9.4743°

从表1中可以看出，利用雷达对目标1角度和距离的估计值，采用最小二乘算法对目标1进行定轨，其精度比采用传统方法有所改善。

仿真5，本发明利用角度和距离信息对目标2采用传统算法与最小二乘算法进行定轨，定轨结果如表2所示：

表2对目标2采用传统算法与最小二乘算法定轨结果

轨道参数	设定值	传统方法定轨	最小二乘算法定轨
				i	65°	64.9917°	64.9998°
Ω	100°	99.9922°	99.9998°
				ω	45°	45.0870°	44.9995°
a	8000km	8000.5km	7999.9km
				e	0.1	0.1001	0.1000
f	-4.9375°	-5.0211°	-4.9369°

从表2中可以看出，采用最小二乘算法对目标2进行定轨，定轨精度比采用传统方法有一定提高。

从表1和表2还可以看出，目标1和目标2均能有效定轨，采用最小二乘原理后，其定轨精度均高于传统方法，但目标2的定轨精度高于目标1。这是由于目标1的信噪比低于目标2，其角度估计误差大于目标2，因此利用角度和距离信息对目标定轨时，目标2的定轨精度比目标1高。

Claims (3)

1.一种天基相控阵雷达空间多目标定轨方法，包括如下步骤：

(1)将相控阵雷达布置在卫星平台上，按脉冲体制发射信号s_e(t)，并接收回波数据，运用展宽调零保形算法，计算第k个目标的最优波束权值w_k，k＝1，...，P，P为目标数；

(2)用第k个目标的最优波束权值w_k乘以回波数据得到波束输出数据y_k，并对该波束输出数据y_k进行脉冲压缩处理，得到第k个目标的距离时延τ_k，k＝1，...，P，P为目标数；

(3)按阵元位置将相控阵分为左上、左下、右上和右下四个象限，把第k个目标的普通波束权值w_{q_k}划分为左上、左下、右上和右下权值，分别为

和用

和分别乘以对应象限的回波数据后，得到左上、左下、右上和右下象限t时刻的波束输出分别为

和

由左上、左下、右上和右下的波束输出构造和波束输出

方位差波束输出和俯仰差波束输出

$y_{q\_k}^{\mathrm{\Σ}}\left(t\right)=y_{q\_k}^{\mathrm{LU}}\left(t\right)+y_{q\_k}^{\mathrm{LD}}\left(t\right)+y_{q\_k}^{\mathrm{RU}}\left(t\right)+y_{q\_k}^{\mathrm{RD}}\left(t\right),$

$y_{q\_k}^X\left(t\right)=y_{q\_k}^{\mathrm{RU}}\left(t\right)+y_{q\_k}^{\mathrm{RD}}\left(t\right)-y_{q\_k}^{\mathrm{LU}}\left(t\right)-y_{q\_k}^{\mathrm{LD}}\left(t\right),$

$y_{q\_k}^Y\left(t\right)=y_{q\_k}^{\mathrm{RU}}\left(t\right)+y_{q\_k}^{\mathrm{LU}}\left(t\right)-y_{q\_k}^{\mathrm{RD}}\left(t\right)-y_{q\_k}^{\mathrm{LD}}\left(t\right);$

(4)以步骤(2)中的第k个目标的距离时延τ_k为中心，从和波束输出

方位差波束输出

和俯仰差波束输出

中均截取2T_P范围数据，T_P为发射脉冲宽度，对所截取的数据分别作距离脉压，得到和距离脉压数据

方位差距离脉压数据俯仰差距离脉压数据

$z_{q\_k}^{\mathrm{\Σ}}\left(\mathrm{\τ}\right)={\∫}_{{\mathrm{\τ}}_k-T_p}^{{\mathrm{\τ}}_k+T_p}{y}_{q\_k}^{\mathrm{\Σ}}\left(t\right)s_e^{*}(\mathrm{\τ}-t)\mathrm{dt},$

$z_{q\_k}^X\left(\mathrm{\τ}\right)={\∫}_{{\mathrm{\τ}}_k-T_p}^{{\mathrm{\τ}}_k+T_p}{y}_{q\_k}^X\left(t\right)s_e^{*}(\mathrm{\τ}-t)\mathrm{dt},$

$z_{q\_k}^Y\left(\mathrm{\τ}\right)={\∫}_{{\mathrm{\τ}}_k-T_p}^{{\mathrm{\τ}}_k+T_p}{y}_{q\_k}^Y\left(t\right)s_e^{*}(\mathrm{\τ}-t)\mathrm{dt},$

其中*表示共轭，k＝1，...，P，P为目标数，在

中寻找到和距离脉压峰值

max(·)表示求最大值，同样可获得方位差距离脉压峰值

和俯仰差距离脉压峰值

(5)利用上述三个波束输出的距离脉压峰值得到俯仰偏角估计

和方位偏角估计

${\mathrm{\ϵ}}_k^{\mathrm{\θ}}=\arctan \left(\mathrm{imag}\left(\frac{P_{q\_k}^Y}{P_{q\_k}^{\mathrm{\Σ}}}\right)\right)\×\frac{\mathrm{\λ}}{2\mathrm{\π}{L}_y\cos {\mathrm{\θ}}_k^B}$

其中arctan(·)表示反正切函数，imag(·)表示取复数虚部操作，

为俯仰差比和，λ为波长，L_x为X方向的等效相位中心坐标，L_y为Y方向的等效相位中心坐标，

为波束指向的俯仰角，

为波束指向的方位角，

为方位差比和，k＝1，...，P，P为目标数；

(6)由俯仰和方位偏角得到目标俯仰角和方位角估计值：

k＝1，...，P，P为目标数；

(7)基于坐标转换原理，利用对空间目标角度和距离的估计，计算空间目标在J2000坐标系下的坐标；

(8)基于Laplace迭代算法，运用最小二乘原理估计空间目标的轨道参数。

2.根据权利要求1所述的天基相控阵雷达空间多目标定轨方法，其特征在于步骤(1)所述的计算第k个目标的最优波束权值w_k，k＝1，...，P，P为目标数，按如下公式计算：

$w_k=w_{q\_k}-C_k{\left(C_k^H{C}_k\right)}^{-1}{C}_k^H{w}_{q\_k}-C_k{\left(C_k^H{C}_k\right)}^{-1}f$

式中，

表示第k个普通波束权矢量，h表示加窗函数，

表示波束指向的导向矢量，表示第k个波束指向，f＝[0 0 0]^T表示约束向量，

表示第k个约束矩阵，k＝1，...，P，P为目标数，H表示共轭转置，F＝[I_N D_x D_y]表示系数矩阵，I_N表示N维单位矩阵，N为阵元数，D_x＝diag(x₁ x₂...x_N)表示由阵元位置x方向坐标(x₁ x₂...x_N)构成的对角矩阵，D_y＝diag(y₁ y₂...y_N)表示由阵元位置y方向坐标(y₁ y₂...y_N)构成的对角矩阵，D_x和D_y用于展宽零陷。

3.根据权利要求1所述的天基相控阵雷达空间多目标定轨方法，其特征在于步骤(8)所述的运用最小二乘原理估计目标星的轨道参数，其过程如下：

8a)基于Laplace迭代算法进行迭代，计算第k个目标t₀时刻的速度矢量v_k0的粗略估值

k＝1，...，P。

将第k个目标在J2000坐标系下的坐标r_k在t₀时刻进行泰勒展开得

$r_k=r_{k0}+v_{k0}\mathrm{\τ}+\frac{1}{2}{a}_{k0}{\mathrm{\τ}}^2+\frac{1}{3!}\frac{d^3{r}_k}{{\mathrm{dt}}^3}{|}_{t=t_0}{\mathrm{\τ}}^3+O\left({\mathrm{\τ}}^4\right),$

其中τ＝t-t₀表示时间间隔，r_k0＝[x_k0 y_k0 z_k0]^T表示t₀时刻第k个目标在J2000坐标系下的坐标，x_k0、y_k0、z_k0为r_k0的三个分量，

为t₀时刻第k个目标在J2000坐标系下的速度矢量，

为t₀时刻第k个目标在J2000坐标系下的加速度矢量。

r_k的泰勒展开式整理后，可写成

r_k＝Frk₀₊Gv_k0，

其中 $F=1-\frac{{\mathrm{\μ\τ}}^2}{2{\left|r_{k0}\right|}^3}+\frac{{\mathrm{\μ\τ}}^3(r_{k0}\·v_{k0})}{2{\left|r_{k0}\right|}^5}+O\left({\mathrm{\τ}}^4\right)$ 为r_k0的系数， $G=\mathrm{\τ}-\frac{{\mathrm{\μ\τ}}^3}{6{\left|r_{k0}\right|}^3}+O\left({\mathrm{\τ}}^4\right)$ 为v_k0的系数，μ表示中心天体的质心引力常数，|r_k0|为t₀时刻空间目标的距离。

由r_k＝ρ_k+R和r_k＝Fr_k0+Gv_k0知，可按

进行迭代计算，迭代时令初值F⁽⁰⁾＝1，G⁽⁰⁾＝τ，其中ρ_k表示第k个目标相对于雷达的向量，R表示已知的探测星在J2000坐标系下的坐标。

8b)利用最小二乘原理，计算精确的第k个目标在J2000坐标系下的坐标t′_k0和速度矢量v′_k0，k＝1，...，P。

由图2相控阵雷达工作场景图中的探测星、空间目标以及地球之间的几何关系式：

$\left[\begin{array}{cc}{r}_{k0}^{\′}& v_{k0}^{\′}\end{array}\right]\left[\begin{array}{cccc}{F}^{\left({\mathrm{\τ}}_1\right)}& F^{\left({\mathrm{\τ}}_2\right)}& \·\·\·& F^{\left({\mathrm{\τ}}_M\right)}\\ G^{\left({\mathrm{\τ}}_1\right)}& G^{\left({\mathrm{\τ}}_2\right)}& \·\·\·& G^{\left({\mathrm{\τ}}_M\right)}\end{array}\right]=\left[\begin{array}{cccc}{\mathrm{\ρ}}_k\left({\mathrm{\τ}}_1\right)+R\left({\mathrm{\τ}}_1\right)& {\mathrm{\ρ}}_k\left({\mathrm{\τ}}_2\right)+R\left({\mathrm{\τ}}_2\right)& \·\·\·& {\mathrm{\ρ}}_k\left({\mathrm{\τ}}_M\right)+R\left({\mathrm{\τ}}_M\right)\end{array}\right]$

得到矩阵[r′_k0 v′_k0]＝DB^H(BB^H)^-1，这里矩阵DB^H(BB^H)^-1的第一列为r′_k0，矩阵DB^H(BB^H)^-1的第二列为v′_k0，

其中D＝[ρ_k(τ₁)+R(τ₁)ρ_k(τ₂)+R(τ₂)...ρ_k(τ_M)+R(τ_M)]表示目标星的位置矢量矩阵，

$B=\left[\begin{array}{cccc}{F}^{\left({\mathrm{\τ}}_1\right)}& F^{\left({\mathrm{\τ}}_2\right)}& \·\·\·& F^{\left({\mathrm{\τ}}_M\right)}\\ G^{\left({\mathrm{\τ}}_1\right)}& G^{\left({\mathrm{\τ}}_2\right)}& \·\·\·& G^{\left({\mathrm{\τ}}_M\right)}\end{array}\right]$ 表示系数矩阵，

表示时间间隔为τ_i时的F值，

表示时间间隔为τ_i时的G值，i＝1，2，...，M，其中M为雷达测量次数，

和

可由8a)中的

目标角度距离信息及时间间隔τ_i计算得到，ρ_k(τ_i)表示时间间隔为τ_i时的ρ_k值，R(τ_i)表示时间间隔为τ_i时的R值，i＝1，2，...，M；

8c)根据步骤8b)得到的精确的第k个目标在J2000坐标系下的坐标r′_k0和速度矢量v′_k0，k＝1，...，P，计算卫星轨道参数，即轨道半长轴a、轨道偏心率e、轨道倾角i、轨道升交点赤经Ω、轨道真近点角f、轨道近地点幅角w：

轨道半长轴：式中，|r′_k0|为r′_k0的模值，|v′_k0|²＝v′_k0 ^Tv′_k0，T表示转置，G_g为万有引力常量，M为地球质量，G_gM＝3.986005×10¹⁴；

轨道偏心率：其中，|h|为矢量h的模值，h＝r′_k0×v′_k0＝[h_x h_y h_z]^T，h_x、h_y和h_z分别为h的三个分量；

轨道倾角：

arccos(·)表示反余弦函数；

轨道升交点赤经： $\mathrm{\Ω}=\arctan (-\frac{h_x}{h_y});$

轨道真近点角： $f=\arccos \left\{\frac{1}{e}\right[\frac{a(1-e^2)}{\left|r_{k0}^{\′}\right|}-1\left]\right\};$

轨道近地点幅角：

arcsin(·)表示反正弦函数，z_k0为矢量r′_k0的Z坐标分量。

Images