一种GNSS模拟器中的通用三维载体运动轨迹生成方法

摘要

一种GNSS模拟器中的通用三维载体运动轨迹生成方法，它有六大步骤：一、预生成载体运动轨迹用户配置文件，设置载体运动分段与轨迹参数；二、在载体运动轨迹模拟开始时，以设置的初始位置为原点建立地理坐标系；三、载体运动分段为空间直线运动类型的实时轨迹模拟计算；四、载体运动分段为空间圆周运动类型的实时轨迹模拟计算；五、对于运动分段之间的衔接，以上一运动段的末状态为当前运动段的初状态；六、将载体即时位置转换到大地坐标系实时显示输出，并将载体即时位置、速度转换到ECEF坐标系支持模拟器进一步功能实现。本发明能建立通用的三维载体运动轨迹产生模型，它在卫星导航技术领域里具有较好的实用价值和广阔的应用前景。

说明书

(一)技术领域

本发明涉及全球卫星导航系统(GNSS)卫星信号模拟器中接收机载体的运动轨迹模拟产生方法，尤其涉及一种GNSS模拟器中的通用三维载体运动轨迹生成方法，属于卫星导航技术领域。

(二)背景技术

GNSS卫星信号模拟器通过产生卫星信号接收机载体的运动轨迹，实时计算出接收机的位置，并结合卫星的实时位置判断卫星可见性，计算卫星信号功率、伪距、多普勒、载波相位等接收机实时观测量，以此实现对导航数据的扩频码调制和载波调制，模拟产生实际应用场景中接收机接收到的卫星导航信号，成为卫星信号接收机研制开发与测试验证中的仿真信号源。

相比直接利用真实卫星信号，GNSS卫星信号模拟器能够提供精确可控、可复现的仿真环境以及非正常情况的测试条件，使接收机的研发效率得以保证。随着美国GPS系统现代化和俄罗斯GLONASS系统改进升级，以及欧洲Galileo系统和我国北斗卫星导航系统的建设(上述全球卫星导航系统可通称GNSS系统)，GNSS模拟器持续受到工业及军事部门的关注。

国外已研制出多种型号的GNSS卫星信号模拟器，如，美国CAST公司的GPS模拟器、英国Spirent公司的GSS系列模拟器、德国IFEN公司的GNSS模拟器，但大都功能固定，所能模拟的导航卫星系统、信号频点与结构、载体动态等严格受限。

模拟用户动态是GNSS模拟器的重要环节之一，精确模拟出接收机载体的运动模式，有利于对接收机的动态性能进行测量评定；而动态导航信号的模拟产生能够支持接收机的捕获跟踪性能测试，以及定位精度等性能的分析评估。接收机载体运动轨迹的模拟可以基于信号模拟器内置的载体运动模型，也可以实测轨迹数据文件的形式载入GNSS信号模拟器，但内置载体运动模型的实现方式更为灵活，在精确性与可重复性方面更为突出。

(三)发明内容

1、目的

在现有的GNSS模拟器产品中，有关接收机载体运动轨迹生成，通常首先划分若干载体类型，在载体类型之下再划分以载体运动特性区分的运动模式，而在运动模式之下支持用户进行与运动模式关联的运动轨迹参数自定义，依据上述流程最终完成载体轨迹的定义，模拟器据此生成载体运动轨迹。以Spirent GSS系列模拟器为例，其载体类型的用户选项分为：简单运动(仅包括平面圆周运动和平面矩形轨迹运动)、车辆、舰船、飞机。在上述载体类型的选项之下，用户需要进行与载体类型关联的运动模式选择，以飞机载体类型为例，依据不同飞行阶段的运动特性，其运动模式分为：爬升、加速直线、匀速直线、匀速转弯、加速转弯、盘旋等。以匀速转弯运动模式为例，用户可定义的运动轨迹参数包括：运动段内飞机航向角的变化量、向心加速度；以加速直线运动模式为例，用户可定义的运动轨迹参数包括：加速段的持续时间或持续距离，加速度或末速度。

扩展性是上述模拟器接收机载体运动轨迹生成技术的主要问题，即载体类型、运动模式、可定义的轨迹参数受到繁琐且特定的限制。由模拟器支持的接收机测试与试验，通常只是利用实际GNSS信号进行测试之前的关键步骤，从而面向特定的应用场景与典型的测试试验成为一般GNSS模拟器的特点之一，这是GNSS模拟器载体动态模式固定并采用概略轨迹模型的主要原因。对于仍在建设中的GNSS系统(如，Galileo系统，北斗卫星导航系统)，通常情况下利用实际的卫星信号并不可行，由于系统建设周期的影响，实际信号测试可能远滞后于模拟器测试，从而对相应的GNSS模拟器提出了更高的要求，在系统建成并具备运行能力之前更加全面的测试试验需要模拟器的载体运动轨迹模拟具备扩展能力、具备反映各种动态的通用性、具备模拟复杂轨迹的能力。

本发明的目的是为了克服现有GNSS模拟器载体运动轨迹生成技术的局限，提供一种GNSS模拟器中的通用三维载体运动轨迹生成方法，它能将任意复杂形式的载体运动轨迹分解为空间直线运动段与空间圆周运动段的组合进行近似，借助姿态描述及多个坐标系的相互转换，建立通用的三维载体运动轨迹产生模型。

2、技术方案

本发明的目的通过以下技术方案实现：一种GNSS模拟器中的通用三维载体运动轨迹生成方法，该方法具体步骤如下：

步骤一：预生成载体运动轨迹用户配置文件，配置文件定义了载体运动分段与轨迹参数设置；

载体初始位置设置为大地坐标系中的坐标(λ₀，h₀)，即经度、纬度和高程。载体初始速度设置为地理坐标系中的东向、北向和天向初速度(v_e，0，v_n，0，v_u，0)；

关于载体运动分段的轨迹参数设置，对于直线运动类型，需要设置的轨迹参数包括：直线运动类型标识，运动段持续时间，东向加速度，北向加速度，天向加速度；对于圆周运动类型，需要设置的轨迹参数包括：圆周运动类型标识，运动段持续时间，运动方向(顺时针或逆时针运动)，圆周半径，线加速度，圆周运动平面方向相对地理坐标系的方向角。

步骤二：在载体运动轨迹模拟开始时，以设置的初始位置(即经度、纬度和高程)为原点建立地理坐标系，初始位置(x₀，y₀，z₀)即为(0，0，0)；

地理坐标系的x轴(又称e轴)沿当地纬线切线指向东，y轴(又称n轴)沿当地经线切线指向北，z轴(又称u轴)沿当地地理垂线(地球参考椭球表面法线)指向天顶。对于动态范围在180km内的载体运动轨迹模拟，这个初始设置的地理坐标系可视为三个轴线方向维持不变，所引起的误差可忽略不计，则在此距离范围内初始的地理坐标系将被保持。在载体运动轨迹模拟结束之前，如果载体动态超出这一距离范围，则需要更新地理坐标系；

步骤三：获取载体运动第一分段的轨迹参数，如果为空间直线运动类型，在地理坐标系中，根据初速度(v_e，0，v_n，0，v_u，0)、加速度(a_e，a_n，a_u)的设置，可以确定载体的即时位置、速度；

运动段持续时间内，根据t_k-1仿真历元时刻载体的位置(x_k-1，y_k-1，z_k-1)、速度(v_e，k-1，v_n，k-1，v_u，k-1)来计算t_k时刻载体位置(x_k，y_k，z_k)、速度(v_e，k，v_n，k，v_u，k)的方法如下：

(a)计算时间间隔Δt

Δt＝t_k-t_k-1

(b)计算t_k-1时刻到t_k时刻的距离变化量(Δs_e，Δs_n，Δs_u)

$\left(\begin{array}{c}\Delta s_e=v_{e,k-1}·\mathrm{\Delta t}+\frac{1}{2}{a}_e·\Delta t^2\\ \Delta s_n=v_{n,k-1}·\mathrm{\Delta t}+\frac{1}{2}{a}_n·\Delta t^2\\ \Delta s_u=v_{u,k-1}·\mathrm{\Delta t}+\frac{1}{2}{a}_u·\Delta t^2\end{array}\right)$

(c)计算t_k时刻载体速度(v_e，k，v_n，k，v_u，k)

$\left(\begin{array}{c}{v}_{e,k}=v_{e,k-1}+a_e·\mathrm{\Delta t}\\ v_{n,k}=v_{n,k-1}+a_n·\mathrm{\Delta t}\\ v_{u,k}=v_{u,k-1}+a_u·\mathrm{\Delta t}\end{array}\right)$

(d)计算t_k时刻载体位置(x_k，y_k，z_k)

$\left(\begin{array}{c}{x}_k=x_{k-1}+\Delta s_e\\ y_k=y_{k-1}+\Delta s_n\\ z_k=z_{k-1}+\Delta s_u\end{array}\right)$

步骤四：如果载体运动第一分段为空间圆周运动类型，根据运动平面相对地理坐标系方向角γ(取运动平面方向与地理坐标系z轴的夹角，以z轴右侧为正角度)的设置，首先确定地理坐标系中圆周运动的二维平面，该平面的x轴即为地理坐标系x轴(e轴)，y轴与x轴垂直；

其次，根据初速度(v_e，0，v_n，0，v_u，0)、圆周半径R和线加速度a的设置，可以确定圆周运动平面内载体的即时位置、速度。运动段持续时间内，计算t_k仿真历元时刻圆周运动平面内载体位置速度(v_x，k，v_y，k)的方法如下：

(a)计算时间间隔Δt

Δt＝t_k-t_k-1

(b)计算t_k-1时刻到t_k时刻的转角变化量Δφ

$\mathrm{\Delta \phi }=-\mathrm{sign}·(\frac{\left|v_{k-1}\right|}{R}·\mathrm{\Delta t}+\frac{1}{2}·\frac{a}{R}·\Delta t^2)$

其中，|.|为绝对值算子，

而$v_{k-1}=\sqrt{v_{e,k-1}^2+v_{n,k-1}^2+v_{u,k-1}^2}$

(c)计算t_k时刻运动平面内载体方位角φ_k

φ_k＝φ_k-1+Δφ

而初始方位${\phi }_{k-1}=\left(\begin{array}{cc}\arctan \frac{v_{x,k-1}}{v_{y,k-1}},& \frac{v_{x,k-1}}{v_{y,k-1}}>0\\ \pi -\arctan \frac{v_{x,k-1}}{v_{y,k-1}},& \frac{v_{x,k-1}}{v_{y,k-1}}<0\end{array}\right)$

(d)计算t_k时刻运动平面内载体速度(v_x，k，v_y，k)

$\left(\begin{array}{c}{v}_{x,k}=\mathrm{sign}·\left|v_k\right|·\sin {\phi }_k\\ v_{y,k}=-\mathrm{sign}·\left|v_k\right|·\cos {\phi }_k\end{array}\right)$

而|v_k|＝|v_k-1|+a·Δt

(e)计算t_k时刻运动平面内载体位置

$\left(\begin{array}{c}{x}_k^p=\mathrm{sign}·R\cos {\phi }_k+\mathrm{sign}·\frac{v_{y0}·R}{\sqrt{v_{x0}^2+v_{y0}^2}}+x_0^p\\ y_k^p=\mathrm{sign}·R\sin {\phi }_k-\mathrm{sign}·\frac{v_{x0}·R}{\sqrt{v_{x0}^2+v_{y0}^2}}+y_0^p\end{array}\right)$

而$\left(\begin{array}{c}{x}_0^p\\ y_0^p\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \gamma & -\sin \lambda \end{array}\right)·\left(\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right)$

最后，将运动平面内载体的位置和速度转换到地理坐标系中的相应运动状态参数，转换方法如下：

(a)确定运动平面坐标系p到地理坐标系L的坐标旋转矩阵

$R_p^L=\left(\begin{array}{cc}1& 0\\ 0& \cos \gamma \\ 0& -\sin \gamma \end{array}\right)$

(b)计算t_k时刻地理坐标系载体速度(v_e，k，v_n，k，v_u，k)

$\left(\begin{array}{c}{v}_{e,k}\\ v_{n,k}\\ v_{u,k}\end{array}\right)=R_p^L·\left(\begin{array}{c}{v}_{x,k}\\ v_{y,k}\end{array}\right)$

(c)计算t_k时刻地理坐标系载体位置(x_k，y_k，z_k)

$\left(\begin{array}{c}{x}_k\\ y_k\\ z_k\end{array}\right)=R_p^L·\left(\begin{array}{c}{x}_k^p\\ y_k^p\end{array}\right)$

步骤五：根据直线或圆周运动类型的用户设置，按照步骤三、步骤四，依次处理后续运动分段。对于运动分段之间的衔接，以上一运动段的末状态为当前运动段的初状态，即上一运动段的载体末位置(x_m，y_m，z_m)、末速度(v_e，m，v_n，m，v_u，m)分别成为当前运动段的初始位置(x₀，y₀，z₀)、初始速度(v_e，0，v_n，0，v_u，0)，即

$\left(\begin{array}{c}{x}_m\\ y_m\\ z_m\end{array}\right)=\left(\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right),\left(\begin{array}{c}{v}_{e,m}\\ v_{n,m}\\ v_{u,m}\end{array}\right)=\left(\begin{array}{c}{v}_{e,0}\\ v_{n,0}\\ v_{u,0}\end{array}\right)$

步骤六：在模拟器的每个仿真历元，将地理坐标系中的载体即时位置转换到大地坐标系，用于用户界面实时显示输出；并将载体即时位置、速度转换到地球直角坐标系(又称ECEF坐标系)以进一步支持卫星可见性判断，伪距、多普勒等信号状态参数计算等模拟器功能实现；

(a)计算t_k时刻大地坐标系中相应的载体位置坐标(λ_k，h_k)

其中，a为地球参考椭球半长轴，e为地球参考椭球偏心率。

便于数值计算起见，考虑R_N和R_M的近似算法，即

而扁率$f=1-\sqrt{1-e^2}.$

(b)计算t_k时刻ECEF坐标系中相应的载体位置坐标(x_e，k，y_e，k，z_e，k)

其中，λ、为载体位置坐标(x_k，y_k，z_k)对应的经、纬度，而(x_oe，y_oe，z_oe)为地理坐标系原点(λ₀，h₀)在ECEF坐标系中的对应坐标，可通过下式计算得到：

而

(c)计算t_k时刻ECEF坐标系中相应的载体速度(v_x，k，v_y，k，v_z，k)

3、优点及功效

由上述技术方案可以看出，与国内外已有的GNSS模拟器载体运动轨迹生成方法相比，本发明提出的载体轨迹生成方法具有以下技术优点：

(1)载体运动轨迹模拟具有通用性，不区分载体类型，认为载体轨迹自身的动态特征即反映载体类型；以空间直线运动、空间圆周运动两种运动模式涵盖适用于GNSS模拟器的卫星信号接收机载体的运动模式。

(2)基于运动轨迹分段分解近似并引入姿态描述，具有适用于GNSS模拟器的复杂载体运动轨迹的模拟能力，对复杂轨迹的近似度仅取决于运动段划分的精细程度，支持载体动态细节的小尺度精确描述。

(3)采用载体运动轨迹用户配置文件方式，载体轨迹参数的配置灵活、便捷，支持用户可完全自定义的载体运动轨迹模拟产生，使载体轨迹模拟具有扩展性。

(四)附图说明

图1为本发明的通用三维载体运动轨迹生成方法流程框图

图2为本发明的载体运动轨迹用户配置文件格式示意图

图3为本发明的载体运动轨迹生成示例图

图中符号说明如下：

+e：东向        +n：北向            +u：天向

(五)具体实施方式

本发明一种适用于GNSS模拟器的通用三维载体运动轨迹生成方法，所述方法流程框图如图1所示。该方法通过GNSS模拟器的软件平台实现，在方法具体实施之前，基于GNSS模拟器的人机交互界面，选择所需仿真的导航系统及频点，设置仿真起始时间和仿真时间长度，选取星历源文件，设置接收机载体动态，选择是否引入将折算为测距误差的各种误差项，等等。当设置接收机载体动态时，需要选择载体运动轨迹用户配置文件，从而进入本发明所述方法的流程。

该方法具体通过以下步骤实现：

步骤一：预生成载体运动轨迹用户配置文件，配置文件定义了载体运动分段与轨迹参数设置。有关用户配置文件的格式说明，如图2所示。

载体初始位置设置为大地坐标系中的坐标(λ₀，h₀)，即经度、纬度和高程。载体初始速度设置为地理坐标系中的东向、北向和天向初速度(v_e，0，v_n，0，v_u，0)。

关于载体运动分段的轨迹参数设置，对于直线运动类型，需要设置的轨迹参数包括：直线运动类型标识，运动段持续时间，东向加速度，北向加速度，天向加速度；对于圆周运动类型，需要设置的轨迹参数包括：圆周运动类型标识，运动段持续时间，运动方向(顺时针或逆时针运动)，圆周半径，线加速度，圆周运动平面方向相对地理坐标系的方向角。

步骤二：在载体运动轨迹模拟开始时，以设置的初始位置(即经度、纬度和高程)为原点建立地理坐标系，初始位置(x₀，y₀，z₀)即为(0，0，0)。

地理坐标系的x轴(又称e轴)沿当地纬线切线指向东，y轴(又称n轴)沿当地经线切线指向北，z轴(又称u轴)沿当地地理垂线(地球参考椭球表面法线)指向天顶。对于动态范围在180km内的载体运动轨迹模拟，这个初始设置的地理坐标系可视为三个轴线方向维持不变，所引起的误差可忽略不计，则在此距离范围内初始的地理坐标系将被保持。在载体运动轨迹模拟结束之前，如果载体动态超出这一距离范围，则需要更新地理坐标系。

步骤三：获取载体运动第一分段的轨迹参数，如果为空间直线运动类型，在地理坐标系中，根据初速度(v_e，0，v_n，0，v_u，0)、加速度(a_e，a_n，a_u)的设置，可以确定载体的即时位置、速度。

运动段持续时间内，根据t_k-1仿真历元时刻载体的位置(x_k-1，y_k-1，z_k-1)、速度(v_e，k-1，v_n，k-1，v_u，k-1)来计算t_k时刻载体位置(x_k，y_k，z_k)、速度(v_e，k，v_n，k，v_u，k)的方法如下：

(a)计算时间间隔Δt

Δt＝t_k-t_k-1

(b)计算t_k-1时刻到t_k时刻的距离变化量(Δs_e，Δs_n，Δs_u)

$\left(\begin{array}{c}\Delta s_e=v_{e,k-1}·\mathrm{\Delta t}+\frac{1}{2}{a}_e·\Delta t^2\\ \Delta s_n=v_{n,k-1}·\mathrm{\Delta t}+\frac{1}{2}{a}_n·\Delta t^2\\ \Delta s_u=v_{u,k-1}·\mathrm{\Delta t}+\frac{1}{2}{a}_u·\Delta t^2\end{array}\right)$

(c)计算t_k时刻载体速度(v_e，k，v_n，k，v_u，k)

$\left(\begin{array}{c}{v}_{e,k}=v_{e,k-1}+a_e·\mathrm{\Delta t}\\ v_{n,k}=v_{n,k-1}+a_n·\mathrm{\Delta t}\\ v_{u,k}=v_{u,k-1}+a_u·\mathrm{\Delta t}\end{array}\right)$

(d)计算t_k时刻载体位置(x_k，y_k，z_k)

$\left(\begin{array}{c}{x}_k=x_{k-1}+\Delta s_e\\ y_k=y_{k-1}+\Delta s_n\\ z_k=z_{k-1}+\Delta s_u\end{array}\right)$

步骤四：如果载体运动第一分段为空间圆周运动类型，根据运动平面相对地理坐标系方向角γ(取运动平面方向与地理坐标系z轴的夹角，以z轴右侧为正角度)的设置，首先确定地理坐标系中圆周运动的二维平面，该平面的x轴即为地理坐标系x轴(e轴)，y轴与x轴垂直。

其次，根据初速度(v_e，0，v_n，0，v_u，0)、圆周半径R和线加速度a的设置，可以确定圆周运动平面内载体的即时位置、速度。运动段持续时间内，计算t_k仿真历元时刻圆周运动平面内载体位置速度(v_x，k，v_y，k)的方法如下：

(a)计算时间间隔Δt

Δt＝t_k-t_k-1

(b)计算t_k-1时刻到t_k时刻的转角变化量Δφ

$\mathrm{\Delta \phi }=-\mathrm{sign}·(\frac{\left|v_{k-1}\right|}{R}·\mathrm{\Delta t}+\frac{1}{2}·\frac{a}{R}·\Delta t^2)$

其中，|.|为绝对值算子，

而$v_{k-1}=\sqrt{v_{e,k-1}^2+v_{n,k-1}^2+v_{u,k-1}^2}$

(c)计算t_k时刻运动平面内载体方位角φ_k

φ_k＝φ_k-1+Δφ

而初始方位${\phi }_{k-1}=\left(\begin{array}{cc}\arctan \frac{v_{x,k-1}}{v_{y,k-1}},& \frac{v_{x,k-1}}{v_{y,k-1}}>0\\ \pi -\arctan \frac{v_{x,k-1}}{v_{y,k-1}},& \frac{v_{x,k-1}}{v_{y,k-1}}<0\end{array}\right)$

(d)计算t_k时刻运动平面内载体速度(v_x，k，v_y，k)

$\left(\begin{array}{c}{v}_{x,k}=\mathrm{sign}·\left|v_k\right|·\sin {\phi }_k\\ v_{y,k}=-\mathrm{sign}·\left|v_k\right|·\cos {\phi }_k\end{array}\right)$

而|v_k|＝|v_k-1|+a·Δt

(e)计算t_k时刻运动平面内载体位置

$\left(\begin{array}{c}{x}_k^p=\mathrm{sign}·R\cos {\phi }_k+\mathrm{sign}·\frac{v_{y0}·R}{\sqrt{v_{x0}^2+v_{y0}^2}}+x_0^p\\ y_k^p=\mathrm{sign}·R\sin {\phi }_k-\mathrm{sign}·\frac{v_{x0}·R}{\sqrt{v_{x0}^2+v_{y0}^2}}+y_0^p\end{array}\right)$

而$\left(\begin{array}{c}{x}_0^p\\ y_0^p\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \gamma & -\sin \lambda \end{array}\right)·\left(\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right)$

最后，将运动平面内载体的位置和速度转换到地理坐标系中的相应运动状态参数，转换方法如下：

(a)确定运动平面坐标系p到地理坐标系L的坐标旋转矩阵

$R_p^L=\left(\begin{array}{cc}1& 0\\ 0& \cos \gamma \\ 0& -\sin \gamma \end{array}\right)$

(b)计算t_k时刻地理坐标系载体速度(v_e，k，v_n，k，v_u，k)

$\left(\begin{array}{c}{v}_{e,k}\\ v_{n,k}\\ v_{u,k}\end{array}\right)=R_p^L·\left(\begin{array}{c}{v}_{x,k}\\ v_{y,k}\end{array}\right)$

(c)计算t_k时刻地理坐标系载体位置(x_k，y_k，z_k)

$\left(\begin{array}{c}{x}_k\\ y_k\\ z_k\end{array}\right)=R_p^L·\left(\begin{array}{c}{x}_k^p\\ y_k^p\end{array}\right)$

步骤五：根据直线或圆周运动类型的用户设置，按照步骤三、步骤四，依次处理后续运动分段。对于运动分段之间的衔接，以上一运动段的末状态为当前运动段的初状态，即上一运动段的载体末位置(x_m，y_m，z_m)、末速度(v_e，m，v_n，m，v_u，m)分别成为当前运动段的初始位置(x₀，y₀，z₀)、初始速度(v_e，0，v_n，0，v_u，0)，即

$\left(\begin{array}{c}{x}_m\\ y_m\\ z_m\end{array}\right)=\left(\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right),\left(\begin{array}{c}{v}_{e,m}\\ v_{n,m}\\ v_{u,m}\end{array}\right)=\left(\begin{array}{c}{v}_{e,0}\\ v_{n,0}\\ v_{u,0}\end{array}\right)$

步骤六：在模拟器的每个仿真历元，将地理坐标系中的载体即时位置转换到大地坐标系，用于用户界面实时显示输出；并将载体即时位置、速度转换到地球直角坐标系(又称ECEF坐标系)以进一步支持卫星可见性判断，伪距、多普勒等信号状态参数计算等模拟器功能实现。

(a)计算t_k时刻大地坐标系中相应的载体位置坐标(λ_k，h_k)

其中，a为地球参考椭球半长轴，e为地球参考椭球偏心率。

便于数值计算起见，考虑R_N和R_M的近似算法，即

而扁率$f=1-\sqrt{1-e^2}.$

(b)计算t_k时刻ECEF坐标系中相应的载体位置坐标(x_e，k，y_e，k，z_e，k)

其中，λ、为载体位置坐标(x_k，y_k，z_k)对应的经、纬度，而(x_oe，y_oe，z_oe)为地理坐标系原点(λ₀，h₀)在ECEF坐标系中的相应坐标，可通过下式计算得到：

而

(c)计算t_k时刻ECEF坐标系中相应的载体速度(v_x，k，v_y，k，v_z，k)

实施例

利用本发明提供的提供通用三维载体运动轨迹生成方法，以模拟产生如图3所示的载体运动轨迹为例，并以WGS-84地球坐标系示例，说明具体实现步骤。

步骤一：GNSS模拟器用户设置载体运动轨迹参数，生成载体运动轨迹用户配置文件。

例如，设置载体初始位置：经度116°，纬度39°，高程100m；设置地理坐标系中载体初始速度：东向100m/s，北向0m/s，天向0m/s；设置载体轨迹分段及参数，该轨迹由10个直线运动分段和圆周运动分段组成，以前5个分段为例说明，分别为：

(a)第一分段为直线运动，东向、北向及天向加速度均为0，历时5s；

(b)第二分段为圆周运动，运动平面相对地理坐标系的方向角为0，圆周半径为200m，顺时针运动方向，线加速度为0，历时20s；

(c)第三分段为直线运动，东向、北向及天向加速度均为0，历时5s；

(d)第四分段为圆周运动，运动平面相对地理坐标系的方向角为45°，圆周半径为200m，顺时针运动方向，线加速度为40m/s²，历时5s；

(e)第五分段为直线运动，东向加速度为65m/s²、北向加速度为34m/s²、天向加速度为-40m/s²，历时5s。

步骤二：读取载体运动轨迹用户配置文件，以设置的初始位置为原点建立地理坐标系。该地理坐标系的x轴(又称e轴)沿39°纬线切线指向东，y轴(又称n轴)沿116°经线切线指向北，z轴(又称u轴)与x和y轴所确定的平面垂直指向天顶。

步骤三：依据步骤一中的运动分段轨迹参数设置，实时模拟载体轨迹、计算载体即时速度及位置坐标。

对于第一分段直线运动，载体从前述地理坐标系原点(0，0，0)起始，以100m/s初速度东向匀速直线运动，历时5s(即Δt＝5s)，则在该运动段的结束时刻，载体的距离变化量、速度和位置坐标分别为：

$\left(\begin{array}{c}\Delta s_e=500\\ \Delta s_n=0\\ \Delta s_u=0\end{array}\right)\left(m\right),$$\left(\begin{array}{c}{v}_{e,m}=100\\ v_{n,m}=0\\ v_{u,m}=0\end{array}\right)(m/s),$$\left(\begin{array}{c}{x}_m=500\\ y_m=0\\ z_m=0\end{array}\right)\left(m\right)$

步骤四：对于第二分段圆周运动，载体从上一运动段的末位置(500，0，0)起始，以上一运动段的末速率为初速率顺时针匀速圆周运动，运动平面维持在前述地理坐标系的水平面内(即γ＝0)，历时20s(即Δt＝20s)，则在该运动段的起始时刻，载体的初始方位、初速率、初始位置分别为：

φ₀＝0.5πrad，$v_0=100m/s,\left(\begin{array}{c}{x}_0^p=500\\ y_0^p=0\end{array}\right)\left(m\right)$

而在该运动段的结束时刻，载体的转角变化量、方位角、运动平面速度和运动平面坐标分别为：

Δφ＝-10rad，φ_m＝-8.4292rad，$\left(\begin{array}{c}{v}_{x,m}=-83.9072\\ v_{y,m}=54.4021\end{array}\right)(m/s),$$\left(\begin{array}{c}{x}_m^p=391.1958\\ y_m^p=-367.8143\end{array}\right)\left(m\right)$

进一步实现至地理坐标系的转换，则坐标旋转矩阵、地理坐标系载体速度及位置坐标分别为：

$R_p^L=\left(\begin{array}{cc}1& 0\\ 0& 1\\ 0& 0\end{array}\right),\left(\begin{array}{c}{v}_{e,m}=-83.9072\\ v_{n,.m}=54.4021\\ v_{u,m}=0\end{array}\right)(m/s),$$\left(\begin{array}{c}{x}_m=391.1958\\ y_m=-367.8143\\ z_m=0\end{array}\right)\left(m\right)$

步骤五：按照步骤三、步骤四，依次处理后续运动分段。

(a)对于第三分段，载体从上一运动段的末位置(391.1958，-367.8143，0)起始，以(-83.9027，54.4021，0)初速度匀速直线运动，历时5s(即Δt＝5s)，则在该运动段的结束时刻，载体的距离变化量、速度和位置坐标分别为：

$\left(\begin{array}{c}\Delta s_e=-419.5358\\ \Delta s_n=272.0106\\ \Delta s_u=0\end{array}\right)\left(m\right),$$\left(\begin{array}{c}{v}_{e,m}=-83.9072\\ v_{n,m}=54.4021\\ v_{u,m}=0\end{array}\right)(m/s),$$\left(\begin{array}{c}{x}_m=-28.3400\\ y_m=54.4021\\ z_m=0\end{array}\right)\left(m\right)$

(b)对于第四分段，载体从上一运动段的末位置(-28.3400，-95.8038，0)起始，以上一运动段的末速率为初速率顺时针加速圆周运动，运动平面与前述地理坐标系的水平面呈45度方向角(即γ＝45)，历时5s(即Δt＝5s)，则在该运动段的起始时刻，载体的初始方位、初速率、初始位置分别为：

φ₀＝4.1372rad，v₀＝100m/s，$\left(\begin{array}{c}{x}_0^p=-28.3400\\ y_0^p=-95.8038\end{array}\right)\left(m\right)$

而在该运动段的结束时刻，载体的转角变化量、方位角、运动平面速度和运动平面坐标分别为：

Δφ＝-5rad，φ_m＝-0.8628rad，$\left(\begin{array}{c}{v}_{x,m}=-227.9064\\ v_{y,m}=-195.0864\end{array}\right)(m/s),$$\left(\begin{array}{c}{x}_m^p=210.5218\\ y_m^p=-79.9270\end{array}\right)\left(m\right)$

进一步实现至地理坐标系的转换，则坐标旋转矩阵、地理坐标系载体速度和位置坐标分别为：

$R_p^l=\left(\begin{array}{cc}1& 0\\ 0& 0.7071\\ 0& -0.7071\end{array}\right),$$\left(\begin{array}{c}{v}_{e,m}=-227.9064\\ v_{n,m}=-137.9469\\ v_{u,m}=137.9469\end{array}\right)(m/s),$$\left(\begin{array}{c}{x}_m=210.5218\\ y_m=-56.5169\\ z_m=56.5169\end{array}\right)\left(m\right)$

(c)对于第五分段，载体从上一运动段的末位置(210.5218，-56.5169，56.5169)起始，以(-227.9064，-137.9469，137.9469)初速度加速直线运动，历时5s(即Δt＝5s)，则在该运动段的结束时刻，载体的距离变化量、速度和位置坐标分别为：

$\left(\begin{array}{c}\Delta s_e=-327.0319\\ \Delta s_n=-264.7344\\ \Delta s_u=189.7344\end{array}\right)\left(m\right),$$\left(\begin{array}{c}{v}_{e,m}=97.0936\\ v_{n,m}=32.0531\\ v_{u,m}=-62.0531\end{array}\right)(m/s),$$\left(\begin{array}{c}{x}_m=-116.5101\\ y_m=-321.2514\\ z_m=246.2514\end{array}\right)\left(m\right)$

步骤六：在模拟器的每个仿真历元，将地理坐标系中的载体即时位置转换到大地坐标系，用于用户界面实时显示输出。由于以WGS-84地球坐标系示例，WGS-84中地球椭球的主要参数定义如下：半长轴a＝6378137m，扁率f＝1/298.257223563，则在上述运动分段的结束时刻，用户界面实时显示的大地坐标系载体位置坐标与地理坐标系载体速度，见表1；而用于支持模拟器进一步功能实现的ECEF坐标系中相应的载体位置坐标及速度，见表2。

表1各运动段结束时刻的大地坐标系载体位置与地理坐标系载体速度

表2各运动段结束时刻的ECEF坐标系载体位置及速度