一种基于GP迭代外推的单星定位方法与流程

本发明属于信号处理领域，具体涉及一种基于gp(gerchberg-papoulis)迭代外推的单星定位方法。

背景技术：

单星定位在军事目标定位和民用导航上有着广泛的应用，信号外推对于通信，雷达，语音与视频信号的恢复均有重要意义。在单星定位中，如果能通过卫星接收到的部分信号g(t)(g(t)为卫星经过辐射源上方附近接收到的完整信号f(t)中的部分信号)，求得整个定位过程的f(t)数据，将会有效缩减单星定位的时间，减少卫星的飞行成本。

信号处理领域中的大部分信号都能通过某种变换映射到某个域中，在变换后的域中处理信号可以检测和提取信号的某些性质。对于信号的外推问题，gerchberg-papoulis外推算法(简称gp算法)利用了长球函数在σ-带限空间内的完备正交性，通过反复应用傅里叶变换及反变换，并截断、替换来迭代出已知区间外的信号。sanz-huang理论是基于gp算法提出的离散估计理论，使得gp算法可以用dft实现。这种离散信号迭代外推的方法可以在有效次迭代内获得误差递减的外推信号，并且计算简便，能够有效提高外推的数据和信号质量。

对于单星定位问题，r.j.webster等人提出了利用单观测站测量目标辐射源来波信号的频率，实现对其定位的方法。实质就是单观测站通过测量信号频率中目标与观测站相对运动的多普勒信息进行定位。虽然采用多普勒定位方法能够完成单星定位，但是多普勒定位方法存在着所需要的定位时间长，卫星飞行时间长的问题，而卫星飞行时间长又会带来飞行成本高的问题。

技术实现要素：

本发明的目的是为解决采用传统多普勒定位方法在单星定位过程中所需要的定位时间长，卫星飞行时间长以及飞行成本高的问题，而提出了一种基于gp迭代外推的单星定位方法。

本发明基于gp外推算法的思想对卫星接收机接收到的部分观测信号进行合理外推，通过原始观测信号与外推得到的信号进行地表辐射源的定位，该信号外推定位算法流程如图1所示，为解决上述技术问题采取的技术方案是：

步骤一、卫星在地表辐射源上空附近(具体飞行位置需要根据轨道的设定和实际卫星的飞行情况确定)飞行时，接收到辐射源发出的部分信号g(t)，将g(t)作为卫星接收到的观测信号；

步骤二、对步骤一中卫星接收到的观测信号g(t)进行后补零处理，即，将接收到的观测信号g(t)后补零至需要的时长，补零段信号的长度可以选择为观测信号g(t)的1倍以上，将补零段的信号记为p(t)；

步骤三、初始化待外推信号为g0(t)，g0(t)由观测信号g(t)与补零段信号p(t)构成；初始化支撑集为φ0，φ0为包含且仅包含观测信号g(t)与补零段信号p(t)的支撑集；

预置迭代次数为n，对初始化的待外推信号g0(t)进行迭代外推，直至满足预置的迭代次数时停止，获得外推出的信号f(t)；

步骤四、对外推出的信号f(t)进行时频分析，得到信号f(t)的时间-频率函数；

再按照时间间隔δt对时间-频率函数进行m次多普勒频率的采样，分别得到在各采样时刻卫星接收信号的多普勒频率估计值；

步骤五、根据卫星星历获得卫星在各采样时刻的位置及速度，再利用卫星在各采样时刻的位置、速度，以及在各采样时刻卫星接收信号的多普勒频率估计值，求解出地表辐射源的位置。

进一步地，所述步骤三的具体过程为：

步骤三一、初始化待外推信号为g0(t)，初始化支撑集为φ0；

步骤三二、对步骤三一中的初始化待外推信号g0(t)进行p阶分数傅里叶变换，经过p阶分数傅里叶变换后获得的信号为b[g0(t)]；

步骤三三、对步骤三二获得的信号b[g0(t)]进行滤波，获得滤波后的信号b[f0(t)]；

步骤三四、对步骤三三得到的滤波后信号b[f0(t)]进行-p阶分数傅里叶变换，得到时域中的映射信号f1(t)；

步骤三五、在初始化支撑集φ0内，利用步骤三一中初始化待外推信号g0(t)对步骤三四获得的信号f1(t)的对应时域部分进行替代，得到时域中第一次迭代获得的信号g1(t)；

替代方法为：

式中，b代表在整个时域范围内；

同时，对初始化支撑集φ0进行更新获得第一次迭代获得的支撑集φ1，支撑集φ1为包含且仅包含信号g1(t)的支撑集；

步骤三六、将步骤三五中获得的信号g1(t)赋值给步骤三一中的初始化待外推信号，将支撑集φ1赋值给初始化支撑集，来重复步骤三二至步骤三五的过程，直至满足预置迭代次数n时停止迭代，得到第n次迭代获得的信号gn(t)，将第n次迭代获得的信号gn(t)作为外推出的信号f(t)。

进一步地，所述步骤五的具体过程为：

在包含辐射源的位置范围内划分网格，获得网格点集∑，(xk,yk,zk)∈∑，(xk,yk,zk)为网格点集∑内的网格点坐标；

其中，中间变量(xj,yj,zj)是卫星在第j个采样时刻的位置，是卫星在第j个采样时刻的速度，c代表光速，ωj代表在第j个采样时刻卫星接收信号的多普勒频率估计值，ω为辐射源发射的信号频率，εj为第j个采样时刻的测频误差，向量rj是向量的模，j＝1,2,…,m；

将公式(2)改写为如下形式：

将简记为gj(xk,yk,zk)，则

ωj＝ωgj(xk,yk,zk)+εj(4)

将全部m个采样点所对应的公式(4)整理成如下的矩阵形式：

ω＝gω+e(5)

式中，ω、g和e均为中间变量，且ω＝[ω1,ω2,…,ωm]^t，e＝[ε1,ε2,…,εm]^t，g＝[g1(xk,yk,zk),g2(xk,yk,zk),…,gm(xk,yk,zk)]^t，上角标t代表矩阵的转置；

为了保证||e||²最小，建立目标函数j(ω)：

j(ω)＝||ω-gω||²(5)

式中，||e||代表e的范数，令则

式中，上角标-1代表矩阵的逆，代表使得目标函数j(ω)取得极小值时的ω；

将网格点集∑中各个网格点的坐标分别代入公式(6)，依次求得各网格点所对应的将使得目标函数j(ω)取值达到最小的所对应的网格点坐标作为辐射源的位置坐标。

x轴，y轴，z轴是指地心地固坐标系的三个坐标轴，地心地固坐标系简称为地心坐标系，是以地心为原点o的坐标系，z轴与地轴平行指向北极点，x轴指向本初子午线与赤道的交点，y轴垂直于xoz平面，x轴，y轴，z轴构成右手坐标系。

进一步地，所述预置迭代次数为n，n的取值为1000次(本发明设定的是1000，可以自行设置)。

更进一步地，所述步骤四中，时间间隔δt的取值为1min。实验中根据实际情况设定。

本发明的有益效果是：本发明提出了一种基于gp迭代外推的单星定位方法，本发明对传统的gp迭代外推算法和多普勒单星定位算法进行优化。由于卫星的接收信号是非平稳信号，仅仅用傅立叶变换不足以分析其显著特征，而运用分数傅立叶变换可以选取信号最集中的角度去分析，通过对卫星已经接收到的信号采用分数傅里叶变换进行合理外推，相比传统傅里叶变换，使得信号在变换域上更加带限，有利于选取合适的滤波器和达到更加良好的外推效果。同时在保证一定的定位精度的要求下，应用信号外推获得了更多的多普勒信息，减缓了定位过程的复杂性问题，缩减了单星定位的时间以及卫星飞行时间，降低了卫星的飞行成本。

附图说明

图1是基于gp算法迭代外推的单星定位流程图；

图2是对sinc信号的第一次外推过程示意图；

图3是对sinc信号的第n次外推过程示意图；

图4是对卫星观测信号进行迭代外推的流程图；

以对sinc信号应用傅里叶变换外推为例，并类比至卫星的观测信号的外推。图2中将有限长观测信号g(t)进行补零，通过离散傅里叶变换dft映射到频域中，得到频域离散信号g(ω)。

之后对变换后的映射域中的信号通过一个低通滤波器进行截断，只保留有限带宽[-σ,σ]内信号f1(ω)；之后对f1(ω)进行离散傅里叶反变换idft，信号长度被展宽，获得信号f1(t)，即，观测区域(-t,t)外出现外推信号，实现过程为对k点频域离散信号f1(ω)进行idft得到时域离散信号f1(t)。

再以观测信号g(t)替代外推信号中相应位置，保留观测信号信息，同时获得观测区域(-t,t)之外的外推信号。从图4的流程知，卫星的接收信号的外推也是按照图2和图3的流程，需要将dft和idft替换成p阶分数傅里叶变换(frft)和-p阶分数傅里叶变换，p阶分数傅里叶变换的表达式为：

其中：

图5为仿真卫星接收信号的外推效果图；

按照图4所示的卫星接收信号的外推流程图，首先对卫星接收信号进行补零，并预置一定的迭代次数，每次迭代时，先对信号进行p阶的分数傅里叶变换frft，再通过低通滤波器h，然后对滤波后的信号进行-p阶的分数傅里叶变换frft，再将变换所得信号在支撑集φ内的部分替换为已知的观测信号，当前次迭代结束；通过不断迭代，重构出完整长度的信号。图5中，前一半为已知的观测信号，后一半为外推出的信号。

图6为对外推重构出的信号进行时频分析后的时间-频率效果图；

图6给出了外推重构出的信号的时间-频率函数关系，设置一定的时间间隔进行采样，可以获得由于卫星与辐射源之间相对运动所产生的多普勒频移ωj，作为多普勒单星定位算法的输入值之一。

图7是对辐射源预估位置第一次划分网格时每个网格点测频误差倒数；

图8是对辐射源预估位置第二次划分网格时每个网格点测频误差倒数；

在卫星接收辐射源发出的信号时，可以根据测向信息等得到辐射源的大致位置，将辐射源所在的大致位置进行网格划分，进行测频误差的解算，这里选择进行了两次划分，第二次的网格比第一次更加细密，由于在多普勒单星定位算法中，测频误差最小的位置就是定位出来的位置，所以在图7和图8中的峰值位置即是解算出的位置。

图9a)为辐射源的实际位置的示意图；

图9b)为解算出的辐射源的最终位置的示意图；

当网格密度满足一定的要求时，输出分别在卫星星下点轨迹两侧令测频误差的模值取得最小值的点，分别为(39.9299°n，116.3879°e)和(33.8964°n，120.6765°e)。将卫星仿真软件中设定的地面辐射源与经过解算得到的两个定位结果标记在地图上，结果如图9a)和图9b)所示，target1(39.9289°n,+116.3883°e)为辐射源实际位置，target3(39.9299°n，116.3879°e)是定位出的辐射源位置，target2(33.8964°n，120.6765°e)为镜像(模糊)点，需要去除。

经由以上分析可知，在多普勒单星定位方法中，引入gp算法迭代外推，在满足一定的定位精度的要求下，可以减少单星定位时间和卫星的飞行成本，根据外推的补零长度以及时间，可以确定整体算法节省的时间；并且可以减少定位过程的复杂性。

仿真过程的参数设置为：

原始信号g(t)的持续时间为在[0,t],t＝6min，观测位置为前一半，补零部分的持续时间为[t,2t],t＝6min；

辐射源实际位置为(39.9289°n,+116.3883°e)，第一次网格划分经纬度密度为1°，第二次网格划分经纬度密度为0.0001°，测频时间间隔为60秒，多普勒频率精度为1hz，卫星位置精度为10^-3m，卫星速度精度为10^-3m/s。

具体实施方式

具体实施方式一、结合图2、图3和图4说明本具体实施方式，基于gp算法迭代外推的单星定位方法，这里为了更加直观清晰说明，以sinc信号应用傅里叶变换外推为例，并类比至卫星的接收信号的外推，它包括以下步骤：

步骤一、对观测信号l进行补零至需要的时长，得到初始化外推信号g和支撑集φ；

步骤二、对g进行迭代外推。预置一定的迭代次数，每次迭代时，先对g进行dft，再通过低通滤波器h，然后对滤波后的信号进行idft，再将反变换所得信号在支撑集φ内的部分替换为已知的观测信号l，更新g，当前次迭代结束；通过不断迭代，重构出完整长度的信号f。

步骤三、由图4可知，卫星的接收信号的外推也是按照步骤二的流程，但是需要将dft和idft替换成p阶和-p阶的frft。然后将满足预置迭代条件的信号f进行时频分析，根据一定的时间间隔进行采样，得到不同时刻卫星接收信号的多普勒频率。

步骤四、根据卫星星历和测频时间得到卫星各时刻位置及速度，根据多普勒单星定位算法，使得测频误差最小的位置即为单星定位的地表辐射源最终位置。

通过卫星仿真软件得到的卫星星历信息如表1所示：

表1

具体实施方式二、结合图1说明本实施方式。本实施方式的一种基于gp算法迭代外推的单星定位方法，所述方法通过以下步骤实现：

a1.对观测信号l进行补零至需要的时长，得到初始化外推信号g和支撑集φ；

a2.预置总迭代次数，将初始化外推信号g通过外推模块，通过不断迭代，重构出完整长度的信号f。

对于观测信号的外推，第i次的实现过程为：

b1.对外推信号进行初始化得到gi；设定迭代次数，获取支撑集φi；

b2.第i次迭代时，对gi进行p阶分数傅里叶变换frft，得到分数域形式gi；

b3.将gi通过低通滤波器hd，经低通滤波得到信号

b4.对进行-p阶分数傅里叶变换frft；

b5.将b4所得信号在支撑集φi内的部分替换为gi，得到gi+1，当前次迭代结束；若迭代终止条件不满足，则返回到b2步，否则，输出重构信号f＝gi+1。

c1.将输出的重构信号f进行时频分析。设置时间间隔，进行m次多普勒频率的采样，作为下一步单星定位解算的输入值之一；

c2.在c1确定的m次采样时刻，根据卫星星历确定卫星在这些时刻的位置(xi,yi,zi)及速度得到

结合地球面参考模型方程可以得到

c3.根据多普勒单星定位算法，使得测频误差εi最小的位置即为单星定位的地表辐射源最终位置。

本发明的上述算例仅为详细地说明本发明的计算模型和计算流程，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。