捷联惯性/天文组合导航(SINS/CNS)作为一种完全自主导航系统，以捷联惯导为主导航方式，辅以天文导航的量测信息实现惯导误差的估计校正，从而获取高精度的导航信息，其组合方式分为简单组合模式、基于陀螺仪漂移校正的组合模式、深组合模式以及基于全面最优校正的组合模式4种。

传统的惯性/天文组合导航采用捷联惯组和星敏感器两者相组合的方法，利用星敏感器高精度定姿特性，通过卡尔曼滤波原理，抑制因陀螺漂移而造成的姿态角发散误差，同时估计出陀螺的常值漂移并进行补偿。文献[1]提出一种基于弹道导弹主动段状态转移阵的误差估计方法，有效地解决了关机点时刻速度和位置误差较大的问题，但由于星敏感器无法从根本上消除加速度计常值偏置等因素对系统的影响，因此无法保证后续飞行中速度、位置误差收敛。针对该问题近些年涌现出一些新方法，旨在改善定位误差，比如文献[2-3]在全面最优组合校正模式下对基于星光折射间接敏感地平的组合导航方法开展研究，实现了高精度定位，但就如何保证星敏感器始终扫描平流层获取折射星的问题工程上尚未实现。文献[4]通过直接敏感地平观测信息并利用高度差法，设计了一种气压高度表辅助的天文/惯性位置组合导航方案，虽然可以提高机载组合导航系统的整体定位精度，但定姿精度远不如星敏感器且不适用于弹道导弹。

雷达高度表(RA)作为一种主动式距离传感器，可分为激光雷达与无线电雷达。随着制导导航与控制技术在深空探测、航天航空的迅猛发展，雷达测距技术己经成功应用于轨道测量、弹道修正、卫星定位、交会对接、巡航制导、航天着陆、大地测量、微波遥感等方面^[5-6]。与合成孔径雷达(SAR)适用于末制导微波成像不同^[7]，RA结构简单，可为稳定飞行平台提供精准的高度测量信息，适用于弹道导弹中段及末段飞行环境。本文针对弹道导弹中段飞行过程中导航精度偏低的问题，在SINS/CNS组合的基础上，提出一种基于雷达高度表的量测模型，引入新的位置量测量，在姿态机动的条件下，利用非线性扩展卡尔曼滤波(EKF)构建SINS/CNS/RA组合导航系统，以获取更高的速度、位置精度。

1 SINS与RA高程差量测模型

导弹在飞行过程中，通过向地球反射表面(海平面或地表平面)发射一系列雷达脉冲，由弹上雷达信号接收机测得回波信号后，可根据时间间隔Δt得到雷达高度表测量高度H′为^[8]

式中：c午磁波传播速度。

为了获得雷达高度表测得的海拔观测高度H_o，需在实测值H′的基础上进行高度补偿，即

式中：补偿值Δh由当地地表海拔等因素决定^[9]。弹道导弹的惯性导引通常以发射惯性系(以下简称l系)作渭航坐标系，如图 1所示，R_e伪地地球半径，R_e0为发射点至地心O_e的距离。

在l系O_l-X_lY_lZ_l下，位置矢量与发射点水平面夹角θ₀满足：

式中：r=; x、y、z分别为SINS位置解算结果在l系下的坐标。则该时刻SINS测得的海拔计算高度H_c为^[10]

H_c与H_o之差ΔH满足如下关系：

式中：

w_h为误差项，受导航误差及量测噪声影响；δx、δy、δz分别为捷联解算与标称弹道的位置误差。

2 弹道环境分析

弹道导弹飞行过程一般分为主动段、中段与再入段。主动段由于发动机点火工作，导弹时刻处于大振动环境且质量变化较大；而再入段中，高速飞行的弹头将会受到强大气动力作用^[11]。因此对于依赖相对静态环境的星敏感器而言，通常将中段作为工作段，此时空间气动力可以忽略，飞行稳定后三轴角速率较小^[12]。由于中段约占全弹道的80%~90%，因此如何提高中段导弹的导航精度对于弹道导弹来说显得尤为重要^[13]。在没有施加姿态机动的情况下，导弹飞行全过程如图 2所示。

星敏感器通过星光成像并获得在成像阵列上星像质心位置，由星图识别及定姿算法，通过安装矩阵，可以获得导弹本体系(以下简称b系)相对地心赤道惯性坐标系(以下简称i系)下的实时姿态^[14]。同时l系到i系的坐标转换矩阵C_lⁱ与发射时刻状态有关，其大小满足：

式中：A为发射时导弹相对于地理北向的发射方位角；A₀=λ₀+A_Υ0，λ₀、φ₀为初始发射位置的地理经纬度，A_Υ0为发射时刻春分点格林时角。进而得到b系相对l系的坐标转换矩阵

式中：ϕ、ψ、θ为绝对姿态角。当弹体在真空环境中无角速率处于绝对静止状态时，C_l^b保持不变。

以导弹质心为坐标原点建立当地铅垂坐标系O_n-X_nY_nZ_n(即北-东-地坐标系，以下简称n系)，进而得到b系相对n系的相对姿态角ϕ_r、ψ_r、θ_r，其坐标转换矩阵为

在n系下根据雷达高度表的捷联安装方式及相应测高算法对导弹进行实时姿态机动，以确保雷达有效测高^[8]。

弹道导弹再入段飞行导航精度取决于地球坐标系(以下简称e系)。通过解算中段飞行位置的地理经纬度λ、φ，有利于导弹的实时导航与控制，从而实现末段精确制导。i系到e系的坐标转换矩阵为

式中：A_Υ伪前时刻春分点格林时角^[15]。则由式(11)即可解算出当前λ、φ的计算主值：

式中：R_l为在l系下的投影向量; R= 。

3 SINS/CNS/RA组合方案 3.1 总体方案

由捷联惯组获得解算绝对姿态角、当前位置地理经纬度及海拔计算高度，由星敏感器获得高精度绝对姿态信息用于补偿惯导陀螺漂移，由雷达高度表结合数字地图获得海拔观测高度用于校正定位误差，设计非线性扩展卡尔曼滤波器，建立SINS/CNS/RA组合导航新算法。整体方案流程图如图 3所示。

图中X为状态变量，ϕ_c、ψ_c、θ_c为捷联解算获得的姿态角，ϕ_o、ψ_o、θ_o为星敏感器观测获得的姿态角，ξ为绝对姿态误差角，η为数学平台失准角，且η=M×ξ^[16]。

在姿态机动条件下，陀螺仪输出角速率ω由机动角速率ω_m、陀螺漂移ε及随机噪声w_g三部分构成，满足

3.2 状态方程

以捷联惯导三轴数学平台失准角η_x、η_y、η_z，l系三轴速度误差δv_x、δv_y、δv_z和位置误差δx、δy、δz，3个陀螺仪漂移ε_x、ε_y、ε_z和3个加速度计偏置∇_x、∇_y、∇_z建立15维状态变量：

式中：陀螺漂移ε包含陀螺随机常值漂移ε_b和一阶马尔可夫过程ε_r，即ε=ε_b+ε_r；加速度计偏置∇包含加速度计常值漂移∇_b和一阶马尔可夫过程∇_r，即∇=∇_b+∇_r^[17]。

考虑J₂项引力摄动项的影响^[6]，摄动状态方程为

式中：a_l为加速度计所敏感的视加速度矢量在l系下的投影；w_g、w_a分别为陀螺仪及加速度计模型的高斯白噪声；v为速度；w_rg、w_ra为对应一阶马尔可夫过程的驱动高斯白噪声；τ_g、τ_a为相关时间常数，当时间常数较大时一阶马尔可夫过程可近似为高斯白噪声；导航坐标系下重力加速度g满足

式中：μ为万有引力常数与地球质量的乘积；R_E呜球赤道半径；R_ix、R_iy、R_iz分别为在i系3个轴的投影坐标。从而建立系统状态方程：

式中：W为系统噪声。

3.3 观测方程

构建4维系统观测量，包括由星敏感器与捷联惯导测得的ξ以及由雷达高度表测得的ΔH两部分，即

于是，观测方程为

式中：

V=[V_ϕ  V_ψ  V_θ  w_h]^T为观测噪声。

3.4 EKF滤波设计

虽然捷联惯导严格意义上属于非线性系统，但为了简化模型便于研究，在姿态误差角为小量的条件下，采用卡尔曼滤波进行线性化处理，通常能够满足精度要求。不过，随着量测时间间隔的增加以及持续噪声的影响，SINS误差积累会导致系统非线性增强。同时由地球椭球模型所引起的重力加速度偏差也对系统摄动状态方程产生重要影响。鉴于此，本文采用EKF算法进行处理。

对于采样输出率高的惯导系统，其连续非线性状态方程近似泰勒展开并略去二阶以上项得

式中：X_k=X(t_k)为k时刻状态变量；A(X)=；为k时刻截断误差；T为系统采样周期。对于离散线性观测方程则采用卡尔曼滤波方法进行处理，最终得到滤波递推方程如下：

式中：P_k、Q_k、R_k、K_k分别为k时刻卡尔曼滤波状态估计阵、系统噪声阵、观测噪声阵、增益矩阵；Φ_{k+1, k}=I+A(X_k)T为k时刻状态传递函数。

由于星敏感器与雷达高度表的采样周期不同，造成观测量Z₁、Z₂数据的采集不同步，因而采用集中式滤波结构，在不断对系统进行时间更新的基础上，适时间断进行局部量测更新。

4 实验仿真 4.1 仿真条件设置

为验证滤波算法，设计标称弹道轨迹发生器作为信息源，其仿真参数如下：

1) 发射点时刻t₀(t₀=0)

(φ₀, λ₀)=(40°N, 100°E)，海拔高度10 m，方位角42°，A_Υ0=0°。

2) 主动段(推力加速度a_t=40 m/s²)

① 垂直飞行段t₀~t₁(t₁=10 s) θ_r=90°

② 程序飞行段t₁~t₂(t₂=30 s)

③ 瞄准飞行段t₂~t₃(t₃=205 s) θ_r=20°

3) 中段(a_t=0)

① 无控调整段t₃~t₄(t₄=250 s) θ_r ≈ 20°

② 姿态控制段

a.t₄~t₅(t₅=400 s)

b.t₅~t₆(t₆=1 810 s)

θ_r= 0°，ϕ_r=0°，弹道最高点海拔为1 017 km。

4) 再入点时刻t₆

(φ_r, λ_r)=(40°6′36″ N, 129°1′45″ W)，海拔高度为70 km，射程为9 828 km。

标称弹道轨迹曲线如图 4所示。

t₀时刻，位置误差5 m，速度误差0.1 m/s，方向失准角6′，水平失准角10″。导航器件精度如表 1所示，σ为误差标准差。

4.2 仿真结果分析

在标称弹道上分别对SINS、SINS/CNS与SINS/CNS/RA进行仿真，结果如图 5和图 6所示。

由图 5可以看出，受主动段影响导弹在关机点时刻已具有较大误差，经过1 810 s时间的飞行，当SINS单独工作时各项误差指标均明显发散；SINS/CNS组合后位置、速度误差发散得到抑制；SINS/CNS/RA组合后，通过一维位置信息的观测校正，经滤波处理，位置、速度误差趋于振荡收敛，导航位置、速度精度显著提高。如表 2所示。

由图 6可知，1 000 s后经滤波处理，SINS/CNS/RA组合姿态角误差分别减小为4.70″、3.83″、4.97″(RMS)。

为了更好验证雷达高度表在地心方向上的位置校正效果，对组合导航系统n系下各分量位置误差仿真图 7所示。

在姿态机动的条件下，雷达高度表可实时获得Z_n轴方向的观测信息，通过n系相对惯性系l系的坐标转换，经由滤波器耦合解算实现对导弹位置信息的校正。从图中可以看出，再入点位置处3轴定位误差分别为1 090、526、6 m，Z_n轴方向的定位精度明显高于其他方向。

独立重复试验20次，统计用本方法得到的再入点估计位置相对标称位置的误差，计算圆概率误差(CEP)。结果如图 8、表 3所示。

由统计特性可以看出，惯导级精度捷联惯组经校正后再入点时刻位置误差保持在2 km以内，CEP为1.2 km。

5 结论

1) 本文方法能够确保导航系统高精度定姿，姿态误差分别为4.70″、3.83″、4.97″(RMS)，趋近于星敏感器安装误差。

2) 在原有SINS/CNS两者组合的基础上引入高度测量信息，逐步校正了主动段造成的累积误差，提高了导航系统定位精度，使得速度、位置误差由原先的发散变为振荡收敛，再入点速度、位置误差分别减小了76.1%和65.0%。

考虑到弹体动态特性以及匹配地形图高度数据精度对系统测量和滤波精度的影响，本文采用精度为50 m(1σ)的雷达高度表测定高程，随着测量仪器精度的提高及测高算法的改进，最终定位精度仍有提升空间。