CN103698793B - 基于软件接收机的gnss信号模拟器仿真角速度范围测量方法 - Google Patents

Abstract

本发明公开了一种基于软件接收机的GNSS信号模拟器仿真角速度范围的测量方法，属于卫星导航领域。本发明在对GNSS信号模拟器输出射频信号进行高保真度采样的基础上，通过软件接收机解算出观测量输出，并利用伪距变化率消除卫星相对运动引入的载波多普勒影响，最后得到角速度的最优估计值。本发明可以根据测试指标对软件接收机参数进行灵活配置，输出的丰富的高更新率高精度观测量，可同时满足测量时间分辨率和测量精度的要求；本方法数据处理原理清晰，不需要进行复杂的坐标转换，非线性最小二乘估计的迭代过程能够保证很快收敛，数据处理简单。

Description

基于软件接收机的GNSS信号模拟器仿真角速度范围测量方法

技术领域

本发明属于卫星导航领域，具体涉及一种基于软件接收机的GNSS信号模拟器仿真角速度范围的测量方法。

背景技术

随着我国北斗系统正式提供区域导航服务，目前已形成GPS、GLONASS、BeiDou和Galileo多系统并存和竞争的局面。为了支持新GNSS星座和信号的设计验证，以及相应的地面监测设备和用户接收机部件的测试评估，GNSS信号模拟器作为关键的模拟测试设备，近年来得到了空前高速的发展。但目前关于对GNSS信号模拟器自身技术性能指标的测试方法尚没有标准可循，各GNSS信号模拟器厂商所声明的技术指标缺乏第三方标准机构的测试和计量认证，模拟器实际所能达到的技术性能难以验证。

目前对于GNSS信号模拟器载体动态类指标的传统测试方法是利用通用的频率测量仪器（频率计数器、频谱分析仪等），在GNSS信号模拟器模拟载体动态场景下，连续测量GNSS信号模拟器输出单载波射频或中频信号的瞬时载波频率，根据瞬时载波多普勒求出导航卫星和所模拟载体接收机天线相位中心之间的瞬时相对径向速度v_r(t)，两者的关系如式(1)所示。

$v_r\left(t\right)=f_d\left(t\right)\·\frac{c}{f_{\mathrm{RF}}}---\left(1\right)$

其中f_RF为GNSS信号模拟器输出信号的标称射频载波频率，c为光速，f_d(t)为瞬时载波多普勒。

旋转载体角速度范围作为GNSS信号模拟器模拟旋转载体动态的重要指标，其测试方法未见公开文献发表，而GNSS信号模拟器的角速度范围指标非常高，现有最先进的通用仪器无法同时满足测量该指标所需的频率分辨率和时间间隔分辨率的要求。以英国Spirent公司GSS8000信号模拟器声明的角速度指标为例（60πrad/s，0.05m的接收机天线杆臂长），若采用频率计数器进行测量，即便每个旋转周期内只测量一次，也需要频率计数器的时间间隔分辨率小于33ms。表1显示了美国Agilent公司目前高端的频率计数器53230A在指定测量时间间隔下的频率分辨率位数，可以看到当测量时间间隔为30ms时，有效频率分辨率位数只有9.2位，无法满足GHz量级的频率测量精度要求。而频谱分析仪等其他通用仪器的频率测量精度更低。即便将射频信号下变频为60MHz的中频信号进行测量，有效测量精度最多也只能达到0.04Hz。

表1测量时间间隔与频率分辨率位数（Agilent53230A）

另外，通用仪器方法只能测量未调制的载波信号，无法体现导航信号不同调制体制的特点，并且需要设计特殊的仿真场景，不能消除卫星相对运动引入的载波多普勒影响。综上，通用仪器方法不能满足GNSS信号模拟器仿真角速度范围测量的需求。

发明内容

本发明的目的是提出一种可行的GNSS信号模拟器仿真角速度范围的测量方法，利用软件接收机可以灵活配置参数和输出丰富的高更新率高精度观测量的优势，通过先验已知的载波多普勒动态辅助跟踪来消除跟踪环路在大动态情况下的动态应力误差，接收机环路参数的设置兼顾测量时间间隔和频率测量精度要求，实现对GNSS信号模拟器仿真旋转载体角速度范围指标的精确测量。

本发明提出的基于软件接收机的GNSS信号模拟器仿真角速度范围测量方法，在对GNSS信号模拟器输出正常信号结构和卫星动态的射频信号进行高保真度采样的基础上，通过软件接收机捕获跟踪单颗导航卫星信号，解算输出导航卫星仰角、方位角、载波多普勒、卫星三维位置坐标等观测量，将载体质心的三维位置坐标作为先验条件，利用卫星速度在径向上的投影消除卫星相对运动引入的载波多普勒影响，利用角度转换和非线性最小二乘拟合方法，计算出角速度的最优估计值。具体是通过以下步骤实现的：

第一步，在GNSS信号模拟器中设定测试动态场景，并进行射频信号的数据采集。

第二步，根据角速度指标配置软件接收机参数，将采集数据送入软件接收机进行处理，并将GNSS信号模拟器输出的载波多普勒作为先验信息，采用载波多普勒外部辅助跟踪的方法消除大动态应力影响，软件接收机最终输出导航卫星仰角、方位角、载波多普勒、卫星三维位置坐标等观测量。

第三步，根据软件接收机输出的卫星三维位置坐标数据及已知的载体质心三维位置坐标，计算伪距变化率，消除卫星相对运动引入的载波多普勒影响；然后根据载波多普勒与旋转角速度的解析关系进行非线性最小二乘估计，最终得到角速度的最优估计值。

本发明提供的方法具有以下优点：

1、本方法可以根据测试指标对软件接收机参数进行灵活配置，利用软件接收机输出丰富的高精度观测量，可同时满足测量时间、分辨率和测量精度的要求。

2、本方法不需要在GNSS信号模拟器中设置卫星静止和单频载波输出的特殊测试场景，也不需要对射频信号进行下变频等额外处理，并且能够适应所有调制体制的导航信号。

3、本方法数据处理原理清晰，不需要进行复杂的坐标转换，非线性最小二乘估计的迭代过程能够保证很快收敛，数据处理简单。

附图说明

图1为本发明提供的基于软件接收机的GNSS信号模拟器仿真载体角速度范围测量方法的流程图；

图2为GNSS信号模拟器测试动态场景设置及数据采集流程示意图；

图3为载体坐标系示意图；

图4为旋转载体剖面示意图；

图5为预检积分时间与载波多普勒测量精度关系图；

图6为软件接收机输出的瞬时载波多普勒曲线；

图7为消除卫星相对运动引入的载波多普勒影响后的瞬时载波多普勒曲线；

图8为最小二乘估计的迭代次数与角速度估计误差的示意图。

具体实施方式

下面将结合附图和具体实施例对本发明作进一步的详细说明。

本发明提供一种基于软件接收机的GNSS信号模拟器仿真角速度范围测量方法，所述测量方法流程如图1所示，具体通过以下步骤实现：

步骤1：在GNSS信号模拟器中设定测试动态场景，并采用高保真度射频信号采集设备进行数据采集，具体过程如下：

首先选择待测的导航系统及频点，选择待测的GNSS信号，在待测GNSS信号模拟器中设置测试动态场景如下：

如图2所示，将GNSS信号模拟器（简称模拟器）的大功率输出端口与高保真度射频信号采集设备（如型号为TG-X3600-2TB的超宽带高速数据连续记录设备）相连接，采集设备与模拟器连接同一高稳定度的外部时钟源；首先设置模拟器仿真单颗可见卫星，正常信号结构和卫星动态（卫星按照正常的运行轨道进行移动），在仿真中关闭星历误差、星钟误差、电离层延迟、对流层延迟、多径等各距离误差项，信号输出功率为最大功率；其次设置载体质心为固定坐标，接收机天线相位中心距离载体质心的距离为接收机天线杆臂长，载体绕其旋转轴自旋，旋转轴指向大地北方向，载体旋转角速度固定且不低于GNSS信号模拟器的角速度范围指标，设置仿真时长；最后启动GNSS信号模拟器输出仿真的射频信号，同时按照预设仿真时长进行数据采集。

步骤2：根据角速度指标配置软件接收机参数，将采集数据送入软件接收机进行处理，采用外部载波多普勒动态辅助跟踪方法消除大动态应力影响，由软件接收机输出导航卫星仰角、方位角、载波多普勒、导航卫星三维位置坐标。具体步骤如下：

（2.1）明确载体旋转角速度与软件接收机输出观测量之间的定量解析关系。

根据载体质心位置建立如图3所示的东北天坐标系：原点O为载体质心；Y轴指向大地北方向（N）；Z轴与过原点O的地球参考椭球面法线重合，指向椭球面外，即天顶方向（U）；X轴与Y轴、Z轴构成右手坐标系，即指向正东方向（E）。图3中，载体旋转轴朝向大地北方向，接收机天线的相位中心（假设接收机天线几何中心与接收机天线相位中心重合）围绕着载体质心O旋转，XOZ平面与接收机天线的旋转剖面重合。

假设在t时刻，导航卫星发射天线的相位中心在S点，接收机天线的相位中心在C点。设导航卫星发射天线的相位中心S和载体质心O的连线OS在XOY平面的投影为OS₁，在XOZ平面的投影为OS₂；导航卫星发射天线的相位中心S和天线相位中心C的连线CS在XOZ平面的投影为CS₂。OS与OS₁的夹角θ₁(t)为导航卫星相对于载体的仰角，OS₁与Y轴的夹角(t)为导航卫星相对于载体的方位角。令OS与OS₂的夹角为α(t)，OS₂与Z轴的夹角为β(t)。

由图3可知，虽然接收机天线的相位中心与载体质心O不重合，但是由于接收机天线杆臂长OC（米级以下）远远小于导航卫星与载体质心的距离OS（数万公里量级），因此可以认为OS近似等于CS；同样可以认为θ₁(t)和(t)分别近似等于软件接收机解算输出的导航卫星相对于接收机天线相位中心的仰角θ(t)和方位角(t)。另外，很容易选择一颗可见星使得OS的投影OS₂长度在万米以上的量级，因此OS₂与CS₂可以近似看作是两条平行线。

图4为对应t时刻，载体的旋转剖面图。图中v₀(t)为接收机天线旋转的瞬时相对切向速度，ω为载体旋转的角速度，ωt则为t时刻接收机天线的相位中心与载体质心O的连线OC相对于Z轴旋转的角度。

结合图3和图4，可知在t时刻，接收机天线旋转的瞬时相对切向速度v₀(t)可以分解为两个速度分量，其中速度分量v₂(t)垂直于CS₂，因而垂直于CS，对接收机天线的瞬时径向速度没有影响；而速度分量v₁(t)平行于CS₂，投影到CS方向上即为此时接收机天线的瞬时相对径向速度v_r(t)，两者满足如下关系：

v_r(t)＝v₁(t)cos[α(t)]    （2）

由此得到对应t时刻，瞬时相对径向速度v_r(t)与接收机天线旋转的瞬时相对切向速度v₀(t)的关系式：

v_r(t)＝v₀(t)sin[ωt-β(t)]cos[α(t)]    （3）

结合公式(1)，进一步得到瞬时载波多普勒f_d(t)为：

f_d(t)＝v₀(t)sin[ωt-β(t)]cos[α(t)]·f_RF/c    （4）

式中α(t)、β(t)可以利用方向向量的几何关系由导航卫星的仰角θ(t)和方位角导出；f_RF为GNSS信号模拟器输出信号的标称射频载波频率，c为光速。

根据图3，在t时刻，在东北天坐标系中导航卫星与载体质心连线的单位入射方向向量用导航卫星仰角θ(t)和方位角表示为：

用α(t)和β(t)表示为：

{cos[α(t)]·sin[β(t)],sin[α(t)],cos[α(t)]·cos[β(t)]}

由此得到如下关系式成立：

对式(5)整理得到：

将式(6)代入式(4)，得到瞬时载波多普勒f_d(t)与导航卫星的仰角θ(t)、方位角和载体旋转角速度ω的关系式：

式中，r为接收机天线杆臂长。

（2.2）接收机的动态辅助跟踪；

由于相对于角速度为ω的正弦变化而言，卫星仰角和方位角均为缓变量，根据式(1)和式(7)得到瞬时相对径向速度为：

v_r(t)＝Aωrsin(ωt-B)    （8）

其中，对于接收机天线相位中心的径向速度而言均为缓变量。

由此得到瞬时相对径向加速度a_r(t)和瞬时相对径向加加速度j_r(t)分别为：

$a_r\left(t\right)=\frac{{\mathrm{dv}}_r\left(t\right)}{\mathrm{dt}}=A{\mathrm{\ω}}^{2}r\cos (\mathrm{\ωt}-B)---\left(9\right)$

$j_r\left(t\right)=\frac{d^2{v}_r\left(t\right)}{{\mathrm{dt}}^2}=-A{\mathrm{\ω}}^{3}r\sin (\mathrm{\ωt}-B)---\left(10\right)$

可见接收机天线相位中心径向速度的一阶和二阶变化率分别与角速度的平方和三次方成正比，并且以同样的角速度随时间呈正弦变化。当角速度值很大时，相应的高动态很容易导致接收机失锁，因此需要利用外部多普勒信息辅助跟踪环路的方法来消除大动态的影响。由于在GNSS信号模拟器的测试中，GNSS信号模拟器能够输出通过数学仿真计算得到的瞬时载波多普勒，因此能够作为已知的先验信息，采用与传统紧组合导航原理相同的载波环外部辅助方法，消除大动态对软件接收机跟踪环路的动态应力误差，使之与热噪声引起的颤动相比可以忽略。

（2.3）设置软件接收机载波跟踪锁相环（PLL，Phase Lock Loop）参数，以保证输出观测量满足测量的时间分辨率和频率精度的要求。

根据式(7)，旋转载体的载波多普勒呈正弦变化，周期为2π/ω。根据正弦信号的采样方法，在一个周期内至少采样4个点才能恢复出原信号，软件接收机输出的载波多普勒时间分辨率应小于π/2ω。而软件接收机锁相环的预检积分时间（也是载波环的更新周期，对应于载波多普勒测量值的时间分辨率）最小可设为1ms，完全能够满足测量要求。

为验证软件接收机锁相环跟踪结果能否同时满足测量精度的要求，做如下分析：在消除大动态条件造成的锁相环动态应力误差后，由于采用了与模拟器同源的高稳时钟基准，接收机载波多普勒的测量精度主要受到热噪声误差的支配。由热噪声造成的锁相环相位抖动误差为：

${\mathrm{\σ}}_{\mathrm{PLL}}=\frac{1}{2\mathrm{\π}}\sqrt{\frac{B_n}{C/N_0}(1+\frac{1}{2T_{\mathrm{coh}}\·C/N_0})}\left(\mathrm{cycle}\right)---\left(11\right)$

相应的载波多普勒测量误差为：

${\mathrm{\σ}}_f=\frac{1}{2\mathrm{\π}{T}_{\mathrm{coh}}}\sqrt{\frac{B_n}{C/N_0}(1+\frac{1}{2T_{\mathrm{coh}}\·C/N_0})}\left(\mathrm{HZ}\right)---\left(12\right)$

式中B_n为锁相环噪声带宽（单位：Hz），T_coh为预检积分时间（单位：s），C/N₀为载噪比（单位：Hz），并且有(C/N₀)_dB单位为dB-Hz。

实测GNSS信号模拟器仿真单星场景下最大信号功率输出时的信号载噪比(C/N₀)_dB可达80dB-Hz左右，一些低噪声模拟器甚至能够达到105dB-Hz以上。图5给出了不同载噪比条件下接收机载波多普勒测量误差（1σ）随预检积分时间的变化。其中PLL环路噪声带宽为B_n=10Hz。可见，在强载噪比条件下，即便取T_coh为1ms，也能保证载波多普勒的测量误差小于0.15Hz；在80dB-Hz的载噪比和10ms的预检积分时间下载波多普勒的测量误差为0.0071Hz（1σ）。

步骤3：根据模拟器所仿真载体质心的三维位置坐标，以及软件接收机所输出的导航卫星三维位置坐标数据，计算伪距变化率，消除卫星相对运动对载波多普勒的影响；然后根据载波多普勒与旋转角速度的解析关系进行非线性最小二乘估计，最终得到角速度的最优估计值。

（1）计算伪距变化率，消除导航卫星相对运动引入的载波多普勒；

在GNSS信号模拟器的正常仿真条件下，导航卫星沿轨道正常运行，软件接收机输出的瞬时载波多普勒中包含了导航卫星相对径向运动带来的影响，如式(13)所示：

f_doppler(t)＝f_d(t)+f_ds(t)    （13）

其中，f_doppler(t)表示软件接收机输出的t时刻对应的瞬时载波多普勒，f_d(t)表示在t时刻，由于载体自旋运动造成的瞬时载波多普勒，f_ds(t)表示在t时刻由于导航卫星相对运动造成的瞬时载波多普勒。

设模拟器所仿真载体质心的三维位置坐标为O:(X_r,Y_r,Z_r)；软件接收机每隔一个预检积分时间间隔，会输出接收时刻坐标系下对应发射时刻的卫星三维位置坐标，设为S:(X_s(t),Y_s(t),Z_s(t))。当在仿真中关闭所有误差项时，两个坐标之间的距离即为伪距，由此可以求出伪距变化率，并转换为载波多普勒。

伪距：

$\mathrm{\ρ}\left(t\right)=\sqrt{{(X_s\left(t\right)-X_r\left(t\right))}^2+{(Y_s\left(t\right)-Y_r\left(t\right))}^2+{(Z_s\left(t\right)-Z_r\left(t\right))}^2}---\left(14\right)$

由于导航卫星相对运动引入的载波多普勒为：

$f_{\mathrm{ds}}(t+T_{\mathrm{coh}})=\frac{-(\mathrm{\ρ}(t+T_{\mathrm{coh}})-\mathrm{\ρ}\left(t\right))}{T_{\mathrm{coh}}}\frac{f_{\mathrm{RF}}}{c}---\left(15\right)$

其中T_coh为预检积分时间。

将对应时刻软件接收机输出的瞬时载波多普勒减去由式(14)和式(15)计算得到的由于导航卫星相对运动引入的瞬时载波多普勒，就可以得到仅由载体自旋引起的瞬时载波多普勒。

（2）非线性最小二乘估计方法得到角速度的估计值；

消除卫星相对运动引入的载波多普勒影响后，式(12)改写为：

f_d(t)＝f_doppler(t)-f_ds(t)＝A(t)ωrsin[ωt-B(t)]f_RF/c    （16）

其中，

对应第i个观测时刻的离散形式为：

f_d(i)＝A(i)ωrsin[ω(t₀+iT_coh)-B(i)]f_RF/c    （17）

＝A(i)ωrsin[ωiT_coh+φ-B(i)]f_RF/c

其中φ表示首个观测时刻（t₀时刻）对应的初相位。

为表达简便，令常数C＝r·f_RF/c，并用下标i表示第i个观测时刻，则式(17)可以改写为：

$\frac{f_{\mathrm{di}}}{C\·A_i}=\mathrm{\ω}\sin (\mathrm{\ωi}{T}_{\mathrm{coh}}+\mathrm{\φ}-B_i)---\left(18\right)$

其中消除卫星相对运动引入载波多普勒影响后得到f_di序列，根据接收机输出的导航卫星仰角和方位角序列可以计算出A_i与B_i序列，未知变量为角速度ω和初始相位φ。

令那么软件接收机输出的N组测量结果可以组成一个二元非线性方程组。利用牛顿迭代法进行线性化处理，假设解(ω,φ)第k次迭代的初始值为(ω_k-1,φ_k-1)，则第k次迭代的线性化方程为：

v_i-ω_k-1sin(ω_k-1iT_coh+φ_k-1-B_i)＝c_i,k-1(ω-ω_k-1)+d_i,k-1(φ-φ_k-1)    （19）

其中：

$\begin{array}{c}{c}_{i,k-1}=\frac{\∂\left[{\mathrm{\ω}}_{k-1}\sin ({\mathrm{\ω}}_{k-1}i{T}_{\mathrm{coh}}+{\mathrm{\φ}}_{k-1}-B_i)\right]}{\∂\mathrm{\ω}}\\ =\sin ({\mathrm{\ω}}_{k-1}i{T}_{\mathrm{coh}}+{\mathrm{\φ}}_{k-1}-B_i)+{\mathrm{\ω}}_{k-1}i{T}_{\mathrm{coh}}\cos ({\mathrm{\ω}}_{k-1}i{T}_{\mathrm{coh}}+{\mathrm{\φ}}_{k-1}-B_i)\end{array}---\left(20\right)$

$d_{i,k-1}=\frac{\∂\left[{\mathrm{\ω}}_{k-1}\sin ({\mathrm{\ω}}_{k-1}i{T}_{\mathrm{coh}}+{\mathrm{\φ}}_{k-1}-B_i)\right]}{\∂\mathrm{\φ}}={\mathrm{\ω}}_{k-1}\cos ({\mathrm{\ω}}_{k-1}i{T}_{\mathrm{coh}}+{\mathrm{\φ}}_{k-1}-B_i)---\left(21\right)$

由所述的二元非线性方程组近似转化为线性方程组，用矩阵形式表示为：

G_N×2,k-1·Δx_2×1,k＝b_N×1,k-1    （22）

其中下脚标表示矩阵的维数，各矩阵分别表示为：

$G_{N\×2,k-1}=\left[\begin{array}{cc}{c}_{1,k-1}& d_{1,k-1}\\ c_{2,k-1}& d_{2,k-1}\\ ...& ...\\ c_{N,k-1}& d_{N,k-1}\end{array}\right]---\left(23\right)$

$\mathrm{\Δ}{x}_{2\×1,k}=\left[\begin{array}{c}\mathrm{\ω}-{\mathrm{\ω}}_{k-1}\\ \mathrm{\φ}-{\mathrm{\φ}}_{k-1}\end{array}\right]---\left(24\right)$

$b_{N\×1,k-1}=\left[\begin{array}{c}{v}_1-{\mathrm{\ω}}_{k-1}\sin ({\mathrm{\ω}}_{k-1}{T}_{\mathrm{coh}}+{\mathrm{\φ}}_{k-1}-B_1)\\ v_2-{\mathrm{\ω}}_{k-1}\sin ({\mathrm{\ω}}_{k-1}2T_{\mathrm{coh}}+{\mathrm{\φ}}_{k-1}-B_2)\\ ...\\ v_N-{\mathrm{\ω}}_{k-1}\sin ({\mathrm{\ω}}_{k-1}N{T}_{\mathrm{coh}}+{\mathrm{\φ}}_{k-1}-B_N)\end{array}\right]---\left(25\right)$

根据式(22)，矩阵Δx_2×1,k的最小二乘解为：

$\mathrm{\Δ}{x}_{2\×1,k}={\left(G_{N\×2,k-1}^T{G}_{N\×2,k-1}\right)}^{-1}{G}_{N\×2,k-1}^T{b}_{N\×1,k-1}---\left(26\right)$

根据接收机输出的观测量信息计算得到v_i、A_i与B_i序列，并假设一组(ω,φ)的初始值(ω₀,φ₀)，分别代入式(20)、式(21)可得到G_N×2,k-1矩阵，代入式(25)可得到b_N×1,k-1矩阵，最终可根据式(26)求出角速度ω和观测初始相位φ的最小二乘估计值(ω_k,φ_k)，一般经过6次迭代后，角速度估值的误差就能够达到稳定，并且远小于1e-6rad/s。将角速度ω的最小二乘估计值与GNSS信号模拟器预设角速度值进行比较，如果偏差小于模拟器的角速度准确度指标，即可判定GNSS信号模拟器所能仿真的角速度范围指标满足要求。

下面以具体测试实例进一步说明本发明。假设待测GNSS信号模拟器的标称角速度范围指标为38rad/s（接收机天线杆臂长0.2m），角速度准确度指标为0.005rad/s，被测信号为GPSL1C/A，利用本发明提供的方法进行角速度测量的具体步骤为：

步骤1：在GNSS信号模拟器中设定测试动态场景，并采用高保真度采集设备进行数据采集。

首先将GNSS信号模拟器大功率射频输出端口与高保真度采集设备相连接，采集设备与模拟器连接同一外部高稳定时钟源。其次选择待测的导航系统为GPS系统，频点为L1（1575.42MHz），选择待测信号为L1C/A信号，设置GNSS信号模拟器仿真单颗可见卫星，正常卫星动态，在仿真中关闭星历误差、星钟误差、电离层延迟、对流层延迟、多径等各距离误差项，信号输出功率为最大功率；然后设置GNSS信号模拟器载体动态场景为：载体质心为固定坐标，接收机天线相位中心距离载体质心的距离为0.2m，载体绕其旋转轴自旋，旋转轴朝向大地北方向，旋转角速度固定为38rad/s，仿真时长180s。最后启动GNSS信号模拟器输出仿真的射频信号，同时按照预设的180s仿真时长进行数据采集，采样率设为100MHz。

步骤2：根据角速度指标配置软件接收机参数，将采集数据送入软件接收机进行处理，采用外部载波多普勒动态辅助跟踪方法消除大动态应力影响，由软件接收机输出卫星仰角、方位角、载波多普勒、卫星三维位置坐标。

根据角速度指标（38rad/s），可知其旋转周期约为165ms。按照正弦信号一个周期内至少采集4个点的要求，软件接收机PLL环路的预检积分时间不大于41.3ms即可。设置环路噪声带宽B_n=10Hz，以采集信号载噪比(C/N₀)_dB等于80dB-Hz计算，不同预检积分时间下的载波多普勒测量精度如表2所示。

表2预检积分时间与载波多普勒测量精度

由表2可知，预检积分时间越大，载波多普勒的测量精度越高，但是一个旋转周期内的采集点数越少，因此为了兼顾测量精度和时间分辨率的需要，这里选取预检积分时间为20ms。

设置GNSS信号模拟器输出以20ms为间隔的载波多普勒数学仿真计算结果，作为软件接收机解算处理的先验辅助信息，软件接收机最终输出时间间隔为20ms的导航卫星仰角、方位角、载波多普勒、导航卫星三维位置坐标等观测量序列。

步骤3：根据模拟器所仿真载体质心的三维位置坐标数据，以及软件接收机输出的导航卫星三维位置坐标数据，计算伪距变化率，消除卫星相对运动对载波多普勒影响；然后根据载波多普勒与旋转角速度的解析关系进行非线性最小二乘估计，最终得到角速度的最优估计值。

根据步骤2得到的观测量序列，做如下数据处理：

（1）根据模拟器所仿真载体质心的三维位置坐标及软件接收机输出的导航卫星三维位置坐标数据，计算伪距变化率，消除卫星相对运动对载波多普勒的影响，得到修正后的载波多普勒序列；

图6所示为消除导航卫星相对运动引入的载波多普勒影响前，软件接收机输出的瞬时载波多普勒曲线，可以看出载波多普勒正弦振荡的包络按一定斜率变化，这就是由于导航卫星相对运动造成的；图7所示为消除导航卫星相对运动对载波多普勒影响后，由于载体自旋造成的瞬时载波多普勒曲线，可以看到，消除了导航卫星相对运动影响后，仅由载体自旋造成的瞬时载波多普勒呈正弦变化。

（2）利用软件接收机输出的导航卫星仰角和方位角序列计算得到A_i与B_i序列，结合修正后的载波多普勒序列计算得到v_i序列，并设初始值(ω₀,φ₀)为（38.01，4.1），利用牛顿迭代法求出角速度的非线性最小二乘估计值。

图8所示为角速度估计值随最小二乘估计迭代次数的变化示意图，表3则示出了各次迭代的角速度相对于仿真真值的估计误差。可见经过6次迭代后角速度估计值的误差已稳定，并小于1e-8rad/s。其中前三次迭代值出现围绕真值的波动是因为初相的迭代初值含较大误差。最终的角速度最小二乘估计值为37.999999989891265rad/s，与待测GNSS信号模拟器设置的角速度指标38rad/s相比较，偏差为-1.01087e-8rad/s，远小于0.005rad/s的误差限，因此可判定待测GNSS信号模拟器的角速度范围指标满足要求。

表3不同迭代次数下的角速度估计误差

本发明的基于软件接收机的GNSS信号模拟器角速度测量方法，通过软件接收机载波多普勒动态辅助跟踪，输出导航卫星的仰角、方位角、载波多普勒、导航卫星三维位置坐标等观测量；根据伪距变化率消除导航卫星相对运动对载波多普勒的影响，并利用非线性最小二乘估计方法解算得到角速度的最优估计值。这种方法可以同时满足测量时间分辨率和测量精度的要求，不需要设置特殊的测试场景，能够适应所有调制体制的导航信号，数据处理分析实现简单。本发明所提供的方法适合各类GNSS信号模拟器的角速度范围指标的测量。

Claims (2)

1.一种基于软件接收机的GNSS信号模拟器仿真角速度范围的测量方法，其特征在于：包括如下步骤：

第一步，在GNSS信号模拟器中设定测试动态场景，并进行射频信号的数据采集；

第二步，根据角速度指标配置软件接收机参数，将采集数据送入软件接收机进行处理，软件接收机最终输出导航卫星仰角、方位角、载波多普勒和导航卫星三维位置坐标；

根据GNSS模拟器中设定的角速度指标ω，可知载波多普勒周期为2π/ω，在一个周期内至少采样四个点恢复出原始信号，确定软件接收机输出的载波多普勒时间分辨率应小于π/2ω；

同时，在消除了大动态条件造成的锁相环应力误差后，由于采用了与模拟器同源的高稳时钟基准，接收机载波多普勒的测量精度受到热噪声误差的支配；由热噪声造成的锁相环相位抖动误差为：

${\mathrm{\σ}}_{\mathrm{PLL}}=\frac{1}{2\mathrm{\π}}\sqrt{\frac{B_n}{C/N_0}(1+\frac{1}{{2T}_{\mathrm{coh}}\·C/N_0})}\left(\mathrm{cycle}\right)---\left(11\right)$

相应的载波多普勒测量误差为：

${\mathrm{\σ}}_f=\frac{1}{2\mathrm{\π}{T}_{\mathrm{coh}}}\sqrt{\frac{B_n}{C/N_0}(1+\frac{1}{{2T}_{\mathrm{coh}}\·C/N_0})}\left(\mathrm{Hz}\right)---\left(12\right)$

式中B_n为锁相环噪声带宽，单位：Hz；T_coh为预检积分时间，单位：s，C/N₀为载噪比，单位：Hz，并且有(C/N₀)_dB单位为dB-Hz；

在强载噪比条件下，(C/N₀)_dB大于80dB-Hz，在1ms～20ms范围内选取软件接收机锁相环预检积分时间，能够兼顾测量精度和时间分辨率的需求；

第三步，根据软件接收机输出的导航卫星三维位置坐标数据及已知的载体质心三维位置坐标，计算伪距变化率，消除卫星相对运动对载波多普勒的影响；然后根据载波多普勒与旋转角速度的解析关系进行非线性最小二乘估计，最终得到角速度的最优估计值，具体为：

软件接收机输出的瞬时载波多普勒中包含了导航卫星相对径向运动带来的影响，如式(13)所示：

f_doppler(t)＝f_d(t)+f_ds(t)   (13)

其中，f_doppler(t)表示软件接收机输出的t时刻对应的瞬时载波多普勒，f_d(t)表示在t时刻，由于载体自旋运动造成的瞬时载波多普勒，f_ds(t)表示在t时刻，由于导航卫星相对运动造成的瞬时载波多普勒；

设模拟器所仿真载体质心的三维位置坐标为O:(X_r,Y_r,Z_r)；软件接收机每隔一个预检积分时间间隔，会输出接收时刻坐标系下对应发射时刻的卫星三维位置坐标，设为S:(X_s(t),Y_s(t),Z_s(t))；当在仿真中关闭所有误差项时，两个坐标之间的距离即为伪距，由此求出伪距变化率，并转换为载波多普勒：

伪距：

$\mathrm{\ρ}\left(t\right)=\sqrt{{(X_s\left(t\right)-X_r\left(t\right))}^2+{(Y_s\left(t\right)-Y_r\left(t\right))}^2+{(Z_s\left(t\right)-Z_r\left(t\right))}^2}---\left(14\right)$

由于导航卫星相对径向运动引入的载波多普勒为：

$f_{\mathrm{ds}}(t+T_{\mathrm{coh}})=\frac{-(\mathrm{\ρ}(t+T_{\mathrm{coh}})-\mathrm{\ρ}\left(t\right))}{T_{\mathrm{coh}}}\frac{f_{\mathrm{RF}}}{c}---\left(15\right)$

其中T_coh为预检积分时间；

将对应时刻软件接收机输出的瞬时载波多普勒减去由式(14)和式(15)计算得到的由于导航卫星相对运动引入的瞬时载波多普勒，就得到仅由载体自旋引起的瞬时载波多普勒：

f_d(t)＝f_doppler(t)-f_ds(t)＝A(t)ωrsin[ωt-B(t)]f_RF/c   (16)

其中，其中，θ(t)和分别为导航卫星仰角和方位角；

对应第i个观测时刻的离散形式为：

$\begin{array}{c}{f}_d\left(i\right)=A\left(i\right)\mathrm{\ω}r\sin [\mathrm{\ω}(t_0+i{T}_{\mathrm{coh}})-B\left(i\right)]f_{\mathrm{RF}}/c\\ =A\left(i\right)\mathrm{\ω}r\sin [\mathrm{\ωi}{T}_{\mathrm{coh}}+\mathrm{\φ}-B\left(i\right)]f_{\mathrm{RF}}/c\end{array}---\left(17\right)$

其中φ表示首个观测时刻t₀对应的初相位，f_RF为GNSS信号模拟器输出信号的标称射频载波频率，c为光速；

为表达简便，令常数C＝r·f_RF/c，并用下标i表示第i个观测时刻，则式(17)改写为：

$\frac{f_{\mathrm{di}}}{C\·A_i}=\mathrm{\ω}\sin (\mathrm{\ωi}{T}_{\mathrm{coh}}+\mathrm{\φ}-B_i)---\left(18\right)$

其中消除卫星相对运动引入载波多普勒影响后得到f_di序列，根据接收机输出的导航卫星仰角和方位角序列计算出A_i与B_i序列，未知变量为角速度ω和采样初始相位φ，r为接收机天线杆臂长；

令那么软件接收机输出的N组测量结果组成一个二元非线性方程组；利用牛顿迭代法进行线性化处理，假设解(ω,φ)第k次迭代的初始值为(ω_k-1,φ_k-1)，则第k次迭代的线性化方程为：

v_i-ω_k-1sin(ω_k-1iT_coh+φ_k-1-B_i)＝c_i,k-1(ω-ω_k-1)+d_i,k-1(φ-φ_k-1)   (19)

其中：

$\begin{array}{c}{c}_{i,k-1}=\frac{\∂\left[{\mathrm{\ω}}_{k-1}\sin ({\mathrm{\ω}}_{k-1}i{T}_{\mathrm{coh}}+{\mathrm{\φ}}_{k-1}-B_i)\right]}{\∂\mathrm{\ω}}\\ =\sin ({\mathrm{\ω}}_{k-1}i{T}_{\mathrm{coh}}+{\mathrm{\φ}}_{k-1}-B_i)+{\mathrm{\ω}}_{k-1}i{T}_{\mathrm{coh}}\cos ({\mathrm{\ω}}_{k-1}i{T}_{\mathrm{coh}}+{\mathrm{\φ}}_{k-1}-B_i)\end{array}---\left(20\right)$

$d_{i,k-1}=\frac{\∂\left[{\mathrm{\ω}}_{k-1}\sin ({\mathrm{\ω}}_{k-1}i{T}_{\mathrm{coh}}+{\mathrm{\φ}}_{k-1}-B_i)\right]}{\∂\mathrm{\φ}}={\mathrm{\ω}}_{k-1}\cos ({\mathrm{\ω}}_{k-1}i{T}_{\mathrm{coh}}+{\mathrm{\φ}}_{k-1}-B_i)---\left(21\right)$

由所述的二元非线性方程组近似转化为线性方程组，用矩阵形式表示为：

G_N×2,k-1·Δx_2×1,k＝b_N×1,k-1   (22)

其中下脚标表示矩阵的维数，各矩阵分别表示为：

$g_{N\×2,k-1}=\left[\begin{array}{cc}{c}_{1,k-1}& d_{1,k-1}\\ c_{2,k-1}& d_{2,k-1}\\ ...& ...\\ c_{N,k-1}& d_{N,k-1}\end{array}\right]---\left(23\right)$

${\mathrm{\Δx}}_{2\×1,k}=\left[\begin{array}{c}\mathrm{\ω}-{\mathrm{\ω}}_{k-1}\\ \mathrm{\φ}-{\mathrm{\φ}}_{k-1}\end{array}\right]---\left(24\right)$

$b_{N\×1,k-1}=\left[\begin{array}{c}{v}_1-{\mathrm{\ω}}_{k-1}\sin ({\mathrm{\ω}}_{k-1}{T}_{\mathrm{coh}}+{\mathrm{\φ}}_{k-1}-B_1)\\ v_2-{\mathrm{\ω}}_{k-1}\sin ({\mathrm{\ω}}_{k-1}2T_{\mathrm{coh}}+{\mathrm{\φ}}_{k-1}-B_2)\\ ...\\ v_N-{\mathrm{\ω}}_{k-1}\sin ({\mathrm{\ω}}_{k-1}N{T}_{\mathrm{coh}}+{\mathrm{\φ}}_{k-1}-B_N)\end{array}\right]---\left(25\right)$

根据式(22)，矩阵Δx_2×1,k的最小二乘解为：

${\mathrm{\Δx}}_{2\×1,k}={\left(G_{N\×2,k-1}^T{G}_{N\×2,k-1}\right)}^{-1}{G}_{N\×2,k-1}^T{b}_{N\×1,k-1}---\left(26\right)$

根据接收机输出的观测量信息计算得到v_i、A_i与B_i序列，并假设一组初始值(ω,φ)的初始值(ω₀,φ₀)，分别代入式(20)、式(21)得到G_N×2,k-1矩阵，代入式(25)式得到b_N×1,k-1矩阵，最终根据式(26)求出角速度ω和观测初始相位φ的最优二乘估计值(ω_k,φ_k)。

2.根据权利要求1所述的一种基于软件接收机的GNSS信号模拟器仿真角速度范围的测量方法，其特征在于：所述的测试动态场景如下：将GNSS信号模拟器的大功率输出端口与高保真度射频信号采集设备相连接，采集设备与模拟器连接同一个外部高稳定时钟源；首先设置模拟器仿真单颗可见卫星，正常信号结构和卫星动态，在仿真中关闭星历误差、星钟误差、电离层延迟、对流层延迟、多径等各距离误差项，信号输出功率为最大功率；其次设置载体质心为固定坐标，接收机天线相位中心距离载体质心的距离为接收机天线杆臂长，载体绕旋转轴自旋，旋转轴指向大地北方向，载体旋转角速度固定且不低于GNSS信号模拟器的角速度范围指标，设置仿真时长；最后启动GNSS信号模拟器输出仿真的射频信号，同时按照预设仿真时长进行数据采集。