一种星载实时掩星预报方法与流程

本发明涉及大气掩星探测技术领域，具体涉及一种星载实时掩星预报方法，利用星载大气掩星探测设备和初始化硬件设备，提前预报掩星事件发生，为快速捕获掩星事件提供初始条件的方法。

背景技术：

gps/met实验的成功，验证了gnss无线电掩星技术对于地球大气探测的具有重大价值。随着gnss无线电掩星技术的发展，掩星技术被越来越多国家重视。掩星事件是指gnss卫星的导航信号经过地球表面的电离层和大气层而产生弯折的物理变化现象。掩星技术是指掩星事件发生时，掩星信号被掩星接收机捕获并跟踪，从而通过分析采集到的观测量来反演信号传播路径中的一些参量的技术。

星载掩星接收机和掩星相对速度较大，因此掩星现象通常持续时间短，因此对于星载掩星探测仪捕获掩星事件，就需快速捕获，利用掩星预报算法，可以提高掩星捕获速度，另外对于不同的轨道掩星事件又具有随机性，因此结合天线指向，利用惯导设备获取当前姿态，通过精准预报，也能够合理设置通道，节省硬件资源。

现有掩星预报算法主要是运用仿真软件计算，然后提前置入设备中，未考虑实际应用中的可见性和实时性，因此准确率不足，从而导致浪费硬件资源。

技术实现要素：

本发明要解决的技术问题是提供一种星载实时掩星预报方法，克服现有技术的不足，通过轨道动力学方程以及坐标系转换，在保证gnss导航星可视的情况下，计算连续弧段内gnss导航星是否满足掩星事件要求，解决了星载掩星探测仪实时精确预报的需求，有效的提高计算效率和预报准确度。

为解决上述技术问题，本发明提供了一种星载实时掩星预报方法，步骤如下：

(1)获取航天器在wgs-84坐标系下的轨道参数

通过掩星探测仪的定位设备获取航天器的位置速度p＝[x，y，z]，v＝[vx，vy，vz]和当前时间t，并根据该位置速度计算航天器的轨道六根数s；

所述轨道六根数s由以下公式：

构建中间变量，

o＝[coss[1]sins[1]0]

od＝[-coss[0]*o[1]coss[0]*o[0]sins[0]]

构建中间变量，

s[4]＝phi-s[2]*sinphi-s[3]*cosphi

给出，其中r，v为位置速度的模；

rr为法向矢量，并将rr单位化，计算公式如下：

rr＝p×v

其中p＝[x，y，z]，v＝[vx，vy，vz]，为航天器的位置速度。

(2)获取航天器和所有gnss导航星在10s后的位置速度

航天器10s后的位置速度获取方法：利用步骤(1)中获取的轨道六根数s，根据轨道动力学建立航天器运动方程，以当前时间t+10s带入，计算10s后的轨道六根数s10，再利用轨道根数转直角坐标系公式，可以获得10s后的航天器位置速度p10，v10；

各gnss导航星10s后的位置速度获取方法：通过掩星探测仪的定位设备获取的各gnss导航星的广播电文，根据各星座公布icd接口文件获得所有导航星在当前时间t之后10s的位置速度pgnss10，vgnss10。

(3)获取所有gnss导航星的俯仰角elv、方位角及切点高度h

方位角和俯仰角ely，计算公式如下：

式中xrelnew，yrelnew，zrelnew是10s后gnss导航星相对于航天器的相对矢量rrel10转化为东北天坐标下的相对矢量rrelnew的各轴分量，由以下转化公式计算而来：

rrelnew＝rrel10*c

式中c为转换矩阵，用来将wgs-84坐标系下航天器的位置速度计算转换到东北天坐标系；rrel10为gnss导航星相对于航天器的相对矢量，由如下公式计算得到：

rrel10＝pgnss10-p10

式中p10为步骤(2)中得到的航天器10s后的位置速度；pgnss10为步骤(2)中得到的各gnss导航星在当前时间t之后10s的位置速度；

切点高度h的计算，通过如下的面积公式来计算：

h*||rrel||＝||p10||*||pgnss10||*sin(θ)

式中θ为leo-地心-gnss之间的夹角，通过三角形余弦定理计算获得；

再计算切点-地心-leo之间的夹角运用的公式如下：

以角度旋转gnss位置坐标到切点位置，计算切点坐标rtanget，并转化为经纬高[lat，lon，hight]。

(4)判断是否符合门限

利用步骤(3)中得到的俯仰角elv和切点高度hight判断是否满足掩星事件的门限；如果满足掩星事件门限，则再利用惯导姿态信息计算航天器天线相对于东北天的方位角，判断导航星是否满足天线方位角约束。

(5)如果步骤(4)中的天线方位角满足约束，则外推多个固定时点，分别继续计算步骤(1)～(3)在选定的固定时点时是否满足步骤(4)；若全部满足，则判定该导航星为掩星。

所述步骤(4)中的掩星事件的门限为：elvmin≤elv≤elvmax，hightmin≤hight≤hightmax；其中：elvmax为掩星的俯仰角门限上限，elvmin为掩星的俯仰角门限下限：所述elvmax及elvmin由公式elvmax＝w/2+ae和elvmin＝elvtanget-2计算，式中w为天线的幅角，ae为天线安装角俯仰角，az为方位角，elvtanget为掩星事件切地球的最低俯仰角。所述最低俯仰角elvtanget由公式计算，式中切点最低高度a为地球半径，h为航天器的高度。

hightmin和hightmax则根据所需要的观测的大气高度选取。如果关注底层大气，可设置为[1，10]km，如果关注中层大气，可设置为[10，80]km。

所述步骤(4)中的天线方位角约束具体为：天线方位角az应该在天线幅角w范围内，即掩星实际的方位角和天线的安装角方位角的差dangle的绝对值小于w/2，天线的安装方位角由天线的安装角和惯导提供的姿态角相加获得。

所述步骤(5)中的固定时点为30s和1min。

本发明与现有技术相比的有益效果：

(1)本发明采用连续弧段判断掩星事件，避免了单次判断造成的误判或者漏判；

(2)本发明通过惯导设备来判断导航星的可见性，增强了预报结果的准确性，方法简单易行，对于实时预报具有重要应用价值；

(3)本发明综合考虑了卫星实际的姿态和掩星天线的安装方向，增强了掩星事件预报的准确性。

附图说明

下面结合附图对本发明的具体实施方式作进一步详细地描述，其中：

图1为本发明在wgs-84坐标系下的几何示意图。

图2为本发明星载实时掩星预报方法在水平面的投影图。

图3为本发明星载实时掩星预报方法的流程图。

图4为本发明的实施场景设计。

具体实施方式

下面结合实施例附图，对本发明的技术方案进行清楚、完整的描述，显然，描述的实施例仅仅是本发明的一个具体的实施例，不是全部的实施例。下述实施例是说明性的，不是限制性的，不能以下述实施例来限定本发明的保护范围。

以下结合图1到图3阐述本发明星载实时掩星预报方法，详细步骤如下：

(1)获取航天器在wgs-84坐标系下的轨道参数

通过掩星探测仪的定位设备获取航天器的位置速度p＝[x，y，z]，v＝[vx，vy，vz]和当前时间t，并根据该位置速度计算航天器的轨道六根数s；

所述轨道六根数s由以下公式：

构建中间变量，

o＝[coss[1]sins[1]0]

od＝[-coss[0]*o[1]coss[0]*o[0]sins[0]]

构建中间变量，

s[4]＝phi-s[2]*sinphi-s[3]*cosphi

给出，其中r，v为位置速度的模；

rr为法向矢量，并将rr单位化，计算公式如下：

rr＝p×v

其中p＝[x，y，z]，v＝[vx，vy，vz]，为航天器的位置速度；

(2)获取航天器和所有gnss导航星在10s后的位置速度

航天器10s后的位置速度获取方法：利用步骤(1)中获取的轨道六根数s，根据轨道动力学建立航天器运动方程，以当前时间t+10s带入，计算10s后的轨道六根数s10，再利用轨道根数转直角坐标系公式，可以获得10s后的航天器位置速度p10，v10；

各gnss导航星10s后的位置速度获取方法：通过掩星探测仪的定位设备获取的各gnss导航星的广播电文，根据各星座公布icd接口文件获得所有导航星在当前时间t之后10s的位置速度pgnss10，vgnss10；

(3)获取所有gnss导航星的俯仰角elv、方位角azm及切点高度h

方位角和俯仰角elv，计算公式如下：

式中xrelnew，yrelnew，zrelnew是10s后gnss导航星相对于航天器的相对矢量rrel10转化为东北天坐标下的相对矢量rrelnew的各轴分量，由以下转化公式计算而来：

rrelnew＝rrel10*c

式中c为转换矩阵，用来将wgs-84坐标系下航天器的位置速度计算转换到东北天坐标系；rrel10为gnss导航星相对于航天器的相对矢量，由如下公式计算得到：

rrel10＝pgnss10-p10

式中p10为步骤(2)中得到的航天器10s后的位置速度；pgnss10为步骤(2)中得到的各gnss导航星在当前时间t之后10s的位置速度；

切点高度h的计算，通过如下的面积公式来计算：

h*||rrel||＝||p10||*||pgnss10||*sin(θ)

式中θ为leo-地心-gnss之间的夹角，通过三角形余弦定理计算获得；

再计算切点-地心-leo之间的夹角

以角度旋转gnss位置坐标到切点位置，计算切点坐标rtanget，并转化为经纬高[lat，lon，hight]；

(4)判断是否符合门限

利用步骤(3)中得到的俯仰角elv和切点高度hight判断是否满足掩星事件的门限；如果满足掩星事件门限，则再利用惯导姿态信息计算航天器天线相对于东北天的方位角，判断导航星是否满足天线方位角azm约束；

(5)如果步骤(4)中的天线方位角满足约束，则外推30s和1min，分别继续计算步骤(1)～(3)在30s和1min时是否满足步骤(4)；若全部满足，则判定该导航星为掩星。

所述步骤(4)中的掩星事件的判断门限为：elvmin≤elv≤elvmax，hightmin≤hight≤hightmax；

其中：elvmax为掩星的俯仰角门限上限，elvmin为掩星的俯仰角门限下限：所述elvmax及elvmin由如下公式计算：

elvmax＝w/2+ae

e1vmin＝elvtanget-2

其中：w为天线的幅角，ae为天线安装角俯仰角，az为方位角，elvtanget为掩星事件切地球的最低俯仰角；所述最低俯仰角elvtanget由如下公式计算：

其中：切点最低高度a为地球半径，h为航天器的高度；

其中hightmin和hightmax根据所需要的观测的大气高度选取，对于底层大气关注的话，可设置为[1，10]km，对于中层大气关注的话，可设置为[10，80]km。

该方法中，步骤(5)中的天线方位角约束具体为：方位角az应该在天线幅角w范围内，即掩星实际的方位角和天线的安装角方位角的差dangle的绝对值小于w/2，天线的安装方位角由天线的安装角和惯导提供的姿态角相加获得。

实施例

仿真场景如图4所示设计，卫星设计轨道高度400km，轨道倾角设计为55度，姿态设计为对地定向，三轴稳定，utc时间设为2017年8月5号4:00:00至2017年8月6号4:00:00，导航星座以gps星座为例，星数选区prn1，航天器接收天线幅角w为40度，天线安装俯仰角为-20，掩星切角为-30，则天线俯仰角门限为[-22，0]，由于场景设计为三轴稳定，所以天线和运动方向保持重合，因此天线方位角门限为[-20，20]，如果天线有姿态机动，需引入惯导设备，加上姿态机动导致的天线方位变化。以时间60s为步长，通过本发明计算全天prn1号gps卫星掩星事件如下表所示：

通过半实物仿真，得出以上数据，可以证实掩星探测仪可以根据预报结果快速、准确跟踪掩星事件。

本发明提供的星载实时掩星预报方法，增强了预报结果的准确性，方法简单易行，对于掩星实时预报具有重要应用价值。

上面结合附图对本发明的实施方式作了详细说明，本发明说明书中未作详细阐述的内容属于本领域专业技术人员的公知技术。