中国指挥与控制学会会刊 
目前,基于卫星平台的对地定位方法一般分为有源定位和无源定位[1]。其中,有源定位通常运用GNSS系统。该系统由卫星发射定位信号,接收机接收信号并根据信号解算自身位置,现行系统有美国GPS、俄罗斯GLONASS、欧洲Galileo和中国北斗导航定位系统[2]。这种定位方法通常只能对合作目标进行定位。为解决非合作目标的定位问题,学者们提出无源定位方法,通过接收目标发射的电磁波对目标进行定位[3]。这种方法一定程度上解决了非合作目标的定位问题。
无源定位一般采用基于参数测量的两步定位法,即在正确分选信号的前提下,首先估计信号参数,然后根据信号参数估计辐射源位置[4⇓-6]。这种方法受到诸多制约,主要包括:特征参数估计误差,比如天线阵误差等造成的角度测量偏差[7],时间基准不一致等导致的时差测量偏差[8],频率基准不一致等导致的频差测量偏差[9];卫星导航、定轨误差,如卫星轨道位置误差、速度误差、时延误差[10-11]。此外,卫星与目标的几何构型变化、电磁波在大气中的折射、系统本身存在的热噪声、空间环境造成的随机误差、观测方程的线性化近似处理等也是产生误差的原因。而且,当存在多组定位观测量时,各观测量通常是分别独立进行估计的,忽略了所有参数必须对应于相同位置的约束,因此,基于参数测量的两步定位法是次优的[12]。
1 信号模型
为简化分析,假设辐射源发射第一个脉冲的起始时刻为0时刻,令各脉冲的初始相位均为零,且脉冲之间是相参的。对于单个辐射源l个接收机的情况,辐射源位于p0,接收机保持时频同步,并于pl,n处接收到第n个脉冲信号。对于短脉冲辐射源而言,多普勒频移可以忽略,因此仅考虑接收信号的时延和相位。接收脉冲信号模型可表示为
rl,n(t)=bl,nsn(t-τl,n) +wl,n(t),0<t<T
其中
τl,n=‖pl,n-p0‖/c
rl,n(t)表示第l个接收机接收到的第n个脉冲信号,bl,n表示接收信号的信道衰减,sn为第n个发射的脉冲信号,f0为载波频率,wl,n(t)为传输噪声。这里假设wl,n(t)为零均值高斯白噪声,分布服从wl(t)~N(0, )。τl,n为信号时延。由于一般线性调频信号的带宽远小于载频,信号载频的延时近似等于信号的相移。
假设信号的中心频率f0已通过载频估计得到,则经过接收机下变频采样处理后,上式可化为
=bl,nDl,n +
其中,Dl,n表示在时域上将 延迟[τl,n·fs]个采样时间,fs为采样率, 、 、 代表采样后的序列,展开表示为:
以上完成了信号时域模型的建立。对式(3)进行离散傅里叶变换,可以得到信号频域模型
(k)=bl,n (k) + (k),0<k<N-1
写成矩阵形式为
rl,n=bl,nal,nsn+wl,n
其中
sn=
wl,n=
rl,n=
al,n=diag
2 相参累积定位算法
2.1 算法流程
定位算法基于极大似然估计原理,首先根据信号的时差和相位差构建代价函数,再在目标区域划分网格点,计算每两个网格点之间的时差、相位差的估计值,用估计值补偿实际接收到的信号。将补偿后的信号代入代价函数中计算,补偿后的信号相似程度越高,代价函数的输出值越大,因此会在辐射源所在网格点处出现峰值。最后通过对网格点的遍历寻找到峰值点,实现对辐射源的定位。步骤如图1 所示。
2.2 代价函数
下面对定位算法中的代价函数进行推导。式(6)可以写为
wl,n=rl,n-bl,nal,nsn
假设各通道接收信号的噪声在时间、空间上是独立同分布的零均值高斯白噪声,即wl,n(t)~N(0, ),也就相当于rl,n-bl,nal,nsn~N(0, )。由此可以得到一个协方差矩阵rl,n-bl,nal,nsn的极大似然估计函数的表达式为
f(,b,s,p)= exp
上式两边取自然对数得
ln(f(,b,s,p))=-ln(π )- ‖rl,n-bl,nal,nsn‖2
对式(10)中的噪声方差 求偏导,然后使偏导数等于0,则噪声方差估计值为
${\hat{ }\!\!\sigma\!\!\text{ }}_{l}^{2}=\overset{i}{\mathop{\underset{n=1}{\mathop \sum }\,}}\,\overset{j}{\mathop{\underset{l=1}{\mathop \sum }\,}}\,\|{{r}_{l}}_{,n}-{{b}_{l}}_{,n}{{a}_{l}}_{,n}{{s}_{n}}{{\|}^{2}}$
代入式(10)中可得代价函数C(p)
C(p)= ‖rl,n-bl,nal,nsn‖2
其中,b为信道衰减,对sn进行归一化,再求式(12)的最小值,则b的极大似然估计值为
${{{\hat{b}}}_{l,n}}={{[{{({{a}_{l}}_{,n}{{s}_{n}})}^{H}}{{a}_{l}}_{,n}{{s}_{n}}]}^{-1}}{{({{a}_{l}}_{,n}{{s}_{n}})}^{H}}{{r}_{l}}_{,n}$
再代入式(12)得
C1(p)= ‖rl,n‖2-
其中,rl,n是独立参数,因此可以将最小化式(14)转化为最大化式(15)
C2(p)= = Qnsn
其中,定义N×N的Hermitian矩阵Qn为:
Qn≜Vn
Vn≜[r1,n,…, rj,n]
对于未知信号的情况,最大化代价函数式(15)可以通过最大化每一个关于sn的二次型来实现。因此要选择一组向量sn使得代价函数最大,也就是需要sn作为Qn的最大特征值λmax{Qn}对应的特征向量。因此代价函数可以转化为
C3(p)= λmax{Qn}
但是Qn的维度是N×N的,会随采样点的增加而增加,导致求特征值的计算量过大。根据矩阵理论,对于一个给定的矩阵X,XXH与XHX的非零特征值相等。因此,对于式(18),可以将N×N维的Qn替换为l×l维的 ,即
≜ Vn
则式(18)等价于
C3(p)= λmax{ }
而辐射源位置就位于代价函数的最大值处,即
${\hat{p}}=\underset{p}{\mathop{\text{argmax}}}\,~({{C}_{3}}(p))$
3 数值仿真
本节主要通过数值仿真,验证定位方法在双星平台的可行性,并且探究其定位性能。
3.1 场景模型
1)卫星动力学计算。人造卫星的运动主要受地球引力影响,此外,在低轨情况下,还受到大气阻力、光压等的影响。这里为简化分析,忽略其他各种摄动力的影响,将卫星绕地球的运动抽象为两个质点在万有引力作用下的运动学问题,即二体问题。
三维空间中,唯一确定物体轨迹需要6个参数,如位置矢量(x,y,z)和速度矢量(vx,vy,vz)可共同确定物体轨迹。此外,用6个轨道根数也可描述。通常的轨道六根数指的是半长轴a、偏心率e、轨道倾角i、近心点幅角ω、升交点赤经Ω和真近点角ϕ。双星的轨道六根数见表1 所示。
表1 双星轨道六根数Tab.1 Orbital elements of double-satellite |
| 轨道六根数 | SAT1 | SAT2 |
|---|---|---|
| a/km | 6 900.687 7 | 6 900.696 2 |
| e | 0.001 769 6 | 0.001 656 1 |
| i/rad | 1.702 292 9 | 1.702 073 0 |
| ω/rad | 2.715 988 9 | 2.762 334 4 |
| Ω/rad | 3.123 603 5 | 3.117 768 9 |
| φ/rad | 3.569 344 6 | 3.522 999 1 |
根据轨道六根数,通过轨道动力学方程运算,可以求出卫星任意时刻的坐标和速度[18]。
2)网格坐标转换。为计算与表示方便,辐射源和划分网格点位置采用基于WGS-84地球模型的经纬高坐标系(LLA)。卫星位置和定位计算过程采用地心球面固连坐标系(ECEF)。因此,需通过式(22)将经纬高坐标系转换到地心球面固连坐标系下进行定位计算。
其中,f为基准椭球体的极扁率,a为基准椭球体的长半轴,具体数值根据WGS-84地球模型数据获得。lon、lat、alt分别代表辐射源的经度、纬度和高程。假设地面辐射源的经度为68°E,纬度为20°N,高程为0,坐标表示为(lon=-68,lat=20,alt=0)。通过对卫星的轨道动力学计算和辐射源坐标转换,得到两颗卫星及辐射源在ECEF和LLA坐标系下的坐标分别如表2 和表3 所示。
表2 卫星和辐射源的初始位置与速度(ECEF)Tab.2 Initial position and velocity of satellites and radiation source (ECEF) |
| 坐标 | SAT1 | SAT2 | 辐射源 |
|---|---|---|---|
| x/km | 2 050.970 | 2 423.374 | 2 246.080 |
| y/km | -6 014.457 | -6 016.345 | -5 559.243 |
| z/km | 2 718.960 | 2 386.824 | 2 167.697 |
| vx/(km/s) | -2.412 047 | -2.386 594 | 0 |
| vy/(km/s) | 2.304 762 | 1.841 402 | 0 |
| vz/(km/s) | 6.903 510 | 7.051 100 | 0 |
表3 卫星和辐射源的初始位置与速度(LLA)Tab.3 Initial position and velocity of satellites and radiation source (LLA) |
| 坐标 | SAT1 | SAT2 | 辐射源 |
|---|---|---|---|
| lon/(°) | -71.170 28 | -68.060 50 | -68 |
| lat/(°) | 23.293 58 | 20.318 53 | 20 |
| alt/km | 536.980 1 | 535.724 2 | 0 |
3)信号数据生成。辐射源发射信号为典型的线性调频脉冲信号,设置载波频率为4 GHz,脉冲重频为5 kHz,空占比为50%,调频带宽为20 MHz。接收机采样率为50 MHz,接收信号为一个完整脉冲。图3 显示了发射信号的时域波形。
3.2 单脉冲直接定位仿真分析
理想条件(无噪声影响,信号未知)下,单脉冲定位仿真得到代价函数分布如图4 所示,代价函数在接近真实发射位置处有一个明显的峰值。
增加噪声使得信噪比(SNR)=-5/(-10)dB,接收脉冲数依旧为1,代价函数计算结果如图5 、图6 所示。
可见,在SNR=-10 dB时,代价函数中的峰值已经被噪声淹没,导致无法成功定位。
3.3 脉冲相参累积直接定位仿真分析
上文中,低信噪比条件下定位遇到困难,因此考虑利用脉冲相参累积提升算法在低信噪比条件下的定位能力。仿真实验中,将脉冲累积数增加到10。考虑卫星飞行速度较快,不能假设卫星始终静止,因此采取“动-停”仿真,即假设卫星在接收脉冲时静止,在其余时间段内运动。通过脉冲相参累积,得到仿真结果如图7 所示。
由仿真结果可以发现,通过累积10个脉冲,在同样SNR=-10 dB的条件下,相比与上一节中峰值被噪声淹没的情况,代价函数中的峰值较为明显。说明通过脉冲相参累积可以有效提升直接定位的效果。
3.4 脉冲相参累积量对定位精度影响仿真分析
在验证脉冲相参累积可以提升直接定位效果后,本小节考虑研究脉冲累积量与信噪比、定位精度之间的关系。在之前同样的仿真场景下,将接收信号的信噪比设置为从-25 dB到15 dB,脉冲累积量分别为1、10和70,蒙特卡洛仿真次数为30次,由此可以得到在不同脉冲累积量下定位误差与信噪比的关系,如图8 所示。
通过统计定位误差可以发现,随着接收信号信噪比的降低,定位精度也随之降低。而在相同信噪比条件下,通过脉冲累积量的增加,可以明显提升定位精度。仿真结果有效说明了脉冲相参累积对定位性能的贡献度。
4 结束语
本文针对卫星对未知地面辐射源定位,提出脉冲相参累积直接定位算法,实现了对地面未知辐射源的直接定位,通过仿真验证了脉冲相参累积直接定位算法的有效性,并进一步分析了脉冲累积量的变化在不同信噪比下对辐射源定位精度的影响。仿真结果表明,在低信噪比条件下,通过脉冲相参累积可以减少定位误差。下一步考虑将算法由对单一辐射源的定位拓展至对多辐射源的联合定位,并提升算法对不同信号的适应性和定位精度。