CN103557872A - 一种rnp中的综合系统误差实时计算方法 - Google Patents

Abstract

本发明公开了一种RNP中的TSE实时计算方法，具体包括以下步骤：步骤一：通过导航方程计算出导航误差概率椭圆。步骤二：利用坐标系旋转算法将导航误差概率椭圆旋转为正椭圆。步骤三：求解正椭圆的外切曲线参数，本发明采用的曲线以外切直线（线切椭圆法）和外切圆（圆切椭圆法）。步骤四：利用外切曲线参数计算实时TSE。步骤五：将实时TSE与RNP规范阈值进行比较，输出告警结果。本发明提出的一种RNP中的TSE实时计算方法具有较小的计算代价和较准确地TSE估计，适用于RNP中的TSE实时计算。

Description

一种RNP中的综合系统误差实时计算方法

技术领域

本发明属于航空导航领域，具体地说，是指一种RNP中的综合系统误差实时计算方法。

背景技术

航空应用对安全水平以及其他需要满足的性能指标提出了极高的要求。在导航方面，RNP（Required Navigation Performance，所需导航性能）对这些性能指标具有明确规定，它主要包含四个性能指标，即精度，完好性，连续性和可用性。其中，精度的提高与实时报警对保障RNP飞行安全具有重要作用。如图1所示，RNP中的精度指的是在某空域或者航路上，机载导航系统应以95%以上的概率保证飞机偏离预计航迹不超过某一规范阈值x。在RNP中，精度记为综合系统误差（TSE，Total System Error），通过对飞机TSE进行实时估计，并将其与RNP规范中对应的阈值进行比较，以便在TSE超过规范阈值时及时报警，保障RNP飞行安全。

当前，TSE计算与报警技术已经被应用于民航RNP飞行中，其一般方法是对所有种类的误差进行组合。如图2所示，TSE主要包括航径定义误差（PDE，Path Definition Error，定义航径与预期航径误差），飞行技术误差（FTE，Fight Technical Error，定义航径与估计位置误差）和导航系统误差（NSE，Navigation System Error，估计位置与真实位置误差）三个部分，其中航径定义误差相对于另外两项较小，可以忽略不计。因此，TSE近似等效于飞行技术误差与导航系统误差的矢量求和。然而，由于在飞行过程中，真实位置是未知的，因此导航系统误差矢量无法获得，这将导致无法对TSE进行直接计算。因此，为了保障飞行安全并实现有效告警，RNP中的TSE实时计算方法的研究非常必要。

TSE计算方法需要对飞行器所有误差的组合进行准确的概率估计。如果估计值过于保守，则有可能造成虚警，反之，则将对飞行安全构成威胁。同时，在RNP中的TSE的计算方法必须高效以便满足航空飞行的实时性需求。目前，TSE的计算方法主要分为两类。第一类通常假设飞行过程中所有误差均服从独立的零均值高斯分布，此时，TSE的标准差为FTE的标准差与NSE的标准差的均方根。然而，这种方法采用先验的FTE与NSE的概率统计值进行计算，只能用于预测，而无法提供实时告警。另外一类方法则将FTE与NSE作为标量进行相加求和获得TSE。尽管这种标量求和法能够获得实时TSE，然而，该方法没有区分TSE的横向或纵向的分量，因此当FTE与NSE方向不同时，标量和将大于矢量和，从而导致计算结果较真实值趋于保守。当采用该计算值与RNP规范的阈值进行比较时，很可能造成虚警，从而影响正常飞行。

发明内容

针对现有RNP中TSE实时计算方法的不足，本发明提出了一种快速有效的RNP中的TSE实时计算方法。RNP的TSE主要包括FTE和NSE两个部分，其中重点考虑相对于垂直预计航迹方向的误差（简称横向误差）。RNP中FTE可以通过导航估计位置与预计航迹横向距离求得，而NSE为一个概率分布，因此需要通过概率论方法进行估计。因此，本发明在误差概率椭圆理论的基础上，提出了另外两种可选择的TSE实时计算方法用于RNP中TSE的实时计算，称作线切椭圆法和圆切椭圆法。结合现有的标量求和法，在卫导导航模式下，对比了这3种TSE实时计算方法的效果。实验结果表明，本发明提出的一种RNP中的TSE实时计算方法具有较小的计算代价和较准确地TSE估计，适用于RNP中的TSE实时计算。

本发明的一种RNP中的TSE实时计算方法，具体包括以下步骤：

步骤一：通过导航方程计算出导航误差概率椭圆。

步骤二：利用坐标系旋转算法将导航误差概率椭圆旋转为正椭圆。

步骤三：求解正椭圆的外切曲线参数，本发明采用的曲线以外切直线（线切椭圆法）和外切圆（圆切椭圆法）。

步骤四：利用外切曲线参数计算实时TSE。

步骤五：将实时TSE与RNP规范阈值进行比较，输出告警结果。

附图说明

图1为现有技术中RNP-x规范示意图；

图2为现有技术中TSE组成示意图；

图3为本发明基于概率椭圆的RNP中的TSE实时计算流程图；

图4为RNP中的TSE实时计算模型示意图；

图5为正椭圆下的TSE实时计算模型示意图；

图6为线切椭圆法示意图；

图7为圆切椭圆法示意图；

图8为三种方法TSE实时计算时间结果；

图9为三种方法TSE实时计算误差。

具体实施方式

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

本发明的一种RNP中的TSE实时计算方法，方法流程图如图3所示，具体包括以下步骤：

步骤一：通过导航方程计算出导航误差概率椭圆。

建立卫星导航非线性观测方程：

${\mathrm{\ρ}}_{\mathrm{\π}}=\left|\right|p_{\mathrm{\π}}^{\mathrm{ENU}}-p_r^{\mathrm{ENU}}\left|\right|+C_b---\left(1\right)$

其中：ρ_π是第π颗可见卫星伪距，

是该卫星在ENU坐标系（东北天坐标系）下的位置，是接收机在ENU坐标系下的位置，C_b是接收机钟差，||·||为欧式距离。通过对式（1）线性化，可获得N_s颗可见卫星观测模型方程：

Δρ＝HΔx    （2）

其中：Δρ是N_s×1伪距残差矢量，H是N_s×4观测矩阵，Δx是状态误差矢量，包括三维位置[Δx,Δy,Δz]^T，其中：Δx,Δy,Δz分别为东北天三个方向的位置误差，和接收机钟差ΔC_b。利用最小二乘法求解式（2），可解得：

Δx＝(H^TH)^-1H^TΔρ    （3）

由于伪距观测量存在误差，通常假定伪距残差矢量服从相同的独立高斯分布，即

其中：

为方差，

是大小为N_s×N_s的单位矩阵，Δρ～N(μ,Σ)表示Δρ服从均值为μ协方差矩阵为Σ的高斯分布。由此可求得Δx协方差为：

$\begin{array}{c}\mathrm{cov}\left(\mathrm{\&Delta;x}\right)=E<{\mathrm{\&Delta;x\&Delta;x}}^T>={\left(H^{T}H\right)}^{-1}{H}^T\mathrm{cov}\left(\mathrm{\&Delta;\&rho;}\right)H{\left(H^{T}H\right)}^{-1}\\ ={\left(H^{T}H\right)}^{-1}{H}^T{\mathrm{\&sigma;}}_p^2{I}_{N_s}H{\left(H^{T}H\right)}^{-1}={\mathrm{\&sigma;}}_p^2{\left(H^{T}H\right)}^{-1}\end{array}---\left(4\right)$

其中：cov(·)表示协方差矩阵。

将cov(Δx)写成展开形式，可得

$\mathrm{cov}\left(\mathrm{\&Delta;x}\right)=\left[\begin{array}{cccc}{\mathrm{\&sigma;}}_x^2& {\mathrm{\&sigma;}}_{\mathrm{yx}}^2& {\mathrm{\&sigma;}}_{\mathrm{zx}}^2& {\mathrm{\&sigma;}}_{C_{b}x}^2\\ {\mathrm{\&sigma;}}_{\mathrm{xy}}^2& {\mathrm{\&sigma;}}_y^2& {\mathrm{\&sigma;}}_{\mathrm{zy}}^2& {\mathrm{\&sigma;}}_{C_{b}y}^2\\ {\mathrm{\&sigma;}}_{\mathrm{xz}}^2& {\mathrm{\&sigma;}}_{\mathrm{yz}}^2& {\mathrm{\&sigma;}}_z^2& {\mathrm{\&sigma;}}_{C_{b}z}^2\\ {\mathrm{\&sigma;}}_{x{C}_b}^2& {\mathrm{\&sigma;}}_{y{C}_b}^2& {\mathrm{\&sigma;}}_{z{C}_b}^2& {\mathrm{\&sigma;}}_{C_b}^2\end{array}\right]---\left(5\right)$

其中：表示误差α的方差，

且α≠β)表示误差α与误差β的协方差。

由上可得状态误差适量Δx的东向和北向两个方向的分量x_EN的协方差 $\mathrm{\&Sigma;}=\left[\begin{array}{cc}{\mathrm{\&sigma;}}_x^2& {\mathrm{\&sigma;}}_{\mathrm{xy}}^2\\ {\mathrm{\&sigma;}}_{\mathrm{xy}}^2& {\mathrm{\&sigma;}}_y^2\end{array}\right],$ 即x_EN服从二维联合高斯分布：

$f\left(x_{\mathrm{EN}}\right)=\frac{1}{2\mathrm{\&pi;}\sqrt{\det \left(\mathrm{\&Sigma;}\right)}}\exp \{-\frac{1}{2}{(x_{\mathrm{EN}}-{\mathrm{\&mu;}}_{\mathrm{EN}})}^T{\mathrm{\&Sigma;}}^{-1}(x_{\mathrm{EN}}-{\mathrm{\&mu;}}_{\mathrm{EN}})\}---6$

其中，μ_EN为导航定位结果估计值的东向和北向两个方向的分量，det(·)为矩阵的行列式。

通过式（6）即可获得导航误差概率椭圆：

(x_EN-μ_EN)^TΣ^-1(x_EN-μ_EN)＝K²    （7）

其中：K为待定等概率椭圆常数。根据TSE的定义，本发明选取95%概率椭圆进行TSE计算，设定K=1.96。由此可获得RNP中的TSE实时计算模型，图如图4所示。飞行器以95%的概率落在以定位点o’为圆心的概率椭圆内，其中a和b为该椭圆的半长轴和半短轴。RNP中的精度要求指的是飞机以95%的概率落在以预计点O为中心的横向误差小于x的范围内。

步骤二：利用坐标系旋转算法将导航误差概率椭圆旋转为正椭圆。

为了方便对椭圆方程进行表示，避免协方差的引入，本发明将坐标系绕预计点进行旋转，将误差概率椭圆转为正椭圆。

将横向误差的协方差Σ进行对角化，可得：

$\mathrm{\&Sigma;}=\left[\begin{array}{cc}{\mathrm{\&sigma;}}_x^2& {\mathrm{\&sigma;}}_{\mathrm{xy}}^2\\ {\mathrm{\&sigma;}}_{\mathrm{xy}}^2& {\mathrm{\&sigma;}}_y^2\end{array}\right]={\mathrm{V\&Lambda;V}}^T---\left(8\right)$

其中：正交矩阵V满足VV^T＝I₂（I₂为2×2单位矩阵）记为 $V=\left[\begin{array}{cc}\cos \mathrm{\&theta;}& -\sin \mathrm{\&theta;}\\ \sin \mathrm{\&theta;}& \cos \mathrm{\&theta;}\end{array}\right],$ 其中θ为参数，可以看出V是一个二维坐标系旋转矩阵，其中θ为坐标系顺时针旋转角度，对角矩阵Λ＝diag(λ_a,λ_b)，其中λ_a和λ_b为矩阵Σ的两个特征值。

假设ENU坐标系经过顺时针旋转角度θ，得到新的坐标系，原坐标系下的向量x在新坐标系下的坐标x'为：

x＝Vx′    （9）

结合式(7),则原始误差椭圆在新坐标系下为正椭圆:

K²=(x_EN-μ_EN)^TΣ^-1(x_EN-μ_EN)=(x′_EN-μ′_EN)^TV^TΣ^-1V(x′_EN-μ′_EN)=(x′_EN-μ′_EN)^T(V^TΣV)^-1(x′_EN-μ′_EN)=(x′_EN-μ′_EN)^TΛ^-1(x′_EN-μ′_EN)    (10)

综上所述，将原ENU坐标系下的任意平面坐标点x经线性变换可得新坐标系下的坐标点x'＝V^Tx，且在新坐标系下，导航定位误差椭圆为正椭圆，如式（10）所示。如图5所示，经过坐标旋转后，预计航迹由y轴转变为直线y＝kx（k＝cotθ为直线斜率），而水平系统误差包络变为直线y＝kx+d和y＝kx-d（d为直线截距,满足

）。不失一般性，这里假设k＞0且d＞0。记误差椭圆的中心为经过坐标转化后的转变为μ'_EN＝(x₀,y₀)，该椭圆半长轴和半短轴分别为

和

由此可得正椭圆的标准方程：

$\frac{{(x-x_0)}^2}{a^2}+\frac{{(y-y_0)}^2}{b^2}=1---\left(11\right)$

步骤三：求解正椭圆的外切曲线参数；

本发明采用的曲线以外切直线（线切椭圆法）和外切圆（圆切椭圆法）为例。

A.线切椭圆法

如图6所示，线切椭圆法的特点采用平行于预计航迹的直线与椭圆最远点相切，从而获得椭圆上距离预计航迹最远的切线。根据椭圆参数，可得切点(x_t,y_t)满足落在椭圆线上的条件：

$\frac{{(x_t-x_0)}^2}{a^2}+\frac{{(y_t-y_0)}^2}{b^2}=1---\left(12\right)$

对式（12）的求全微分，得：

$\frac{x_t-x_0}{a^2}{\mathrm{dx}}_t+\frac{y_t-y_0}{b^2}{\mathrm{dy}}_t=0---\left(13\right)$

又切线与预计航迹平行，因此切线斜率满足：

$\frac{{\mathrm{dy}}_t}{{\mathrm{dx}}_t}=k---\left(14\right)$

i)当k＝0或∞时，椭圆的轴方向与预计航迹平行或垂直，由椭圆的性质易得：

TSE＝b+y₀或TSE＝a+x₀    （15）

此时，线切椭圆法与标量求和法结果相同。

ii)当k≠0和∞时，联立求解式（12）-式（14）可以得到切点：

$\left\{\begin{array}{c}{y}_t=y_0\&PlusMinus;\frac{b^2}{\sqrt{k^2{a}^2+b^2}}\\ x_t=x_0\stackrel{\&OverBar;}{+}\frac{{\mathrm{ka}}^2}{\sqrt{k^2{a}^2+b^2}}\end{array}---\left(16\right)\right.$

B.圆切椭圆法

如图7所示，圆切椭圆法采用预计航迹点（坐标原心）的最小半径圆O(r)对导航误差椭圆进行包络。因此，该包络圆应过椭圆上距离预计航迹点最远的点(x_f,y_f)，即：

$\begin{array}{c}\underset{(x_f,y_f)}{\max }{r}^2=x_f^2+y_f^2\\ s.t.\frac{{(x_f-x_0)}^2}{a^2}+\frac{{(y_f-y_0)}^2}{b^2}=1\end{array}---\left(17\right)$

其中：r为圆O的半径，s.t.为约束条件。式（17）表明在(x_f,y_f)满足约束条件

的条件下，求圆O的半径r的最大值。

为方便计算，将x_f与y_f表示成极坐标的形式，可得x_f＝x₀+acosφ，y_f＝y₀+bsinφ，其中φ为极坐标参数，将其代入式（17）可转化得：

$\underset{\mathrm{\&phi;}}{\max }{r}^2\left(\mathrm{\&phi;}\right)={(x_0+a\cos )}^2+{(y_0+b\sin )}^2---\left(18\right)$

求r²(φ)相对于φ的导数，且满足极大值必要条件

即

因此可以得到：

$\begin{array}{c}\frac{\&PartialD;r\left(\mathrm{\&phi;}\right)}{\&PartialD;\mathrm{\&phi;}}=-a(x_0+a\cos )\sin \mathrm{\&phi;}+b(y_0+b\sin \mathrm{\&phi;})\cos \mathrm{\&phi;}\\ =(b^2-a^2)\cos \mathrm{\&phi;}\sin \mathrm{\&phi;}-{\mathrm{ax}}_0\sin \mathrm{\&phi;}+{\mathrm{by}}_0\cos \mathrm{\&phi;}=0\end{array}---\left(19\right)$

i)当a＝b时，误差椭圆退化为正圆，此时式（19）可化简为：

-x₀sinφ+y₀cosφ＝0    （20）

将式（20）代入到式（18）中，可得：

$r\left(\mathrm{\&phi;}\right)=\sqrt{x_0^2+y_0^2}+a---\left(21\right)$

ii)当a≠b时，令 $t=\tan \left(\frac{\mathrm{\&phi;}}{2}\right),$ 则 $\sin \left(\mathrm{\&phi;}\right)=\frac{2t}{1+t^2},\cos \left(\mathrm{\&phi;}\right)=\frac{1-t^2}{1+t^2},$ 式（19）可转化为

$\frac{\&PartialD;r\left(\mathrm{\&phi;}\right)}{\&PartialD;\mathrm{\&phi;}}=(b^2-a^2)\frac{2t(1-t^2)}{{(1+t^2)}^2}-a{x}_0\frac{2t}{1+t^2}+b{y}_0\frac{1-t^2}{1+t^2}=0---\left(22\right)$

式（22）可转化为一元四次方程

by₀t⁴+2(ax₀-a²+b²)t³+2(ax₀+a²-b²)t-by₀＝0     （23）

ii.i)当y₀＝0时，式（23）退化为一元三次方程

(ax₀-a²+b²)t³+(ax₀+a²-b²)t＝0    （24）

易得式（23）的三个解分别为 $t_1=0,t_2=\sqrt{\frac{-a{x}_0-a^2+b^2}{a{x}_0-a^2+b^2}},t_3=-\sqrt{\frac{-a{x}_0-a^2+b^2}{a{x}_0-a^2+b_2}}.$

ii.ii)当y₀≠0时，将式（23）首项系数化为1，可得：

$t^4+\frac{2(a{x}_0-a^2+b^2)}{b{y}_0}{t}^3+\frac{2(a{x}_0+a^2-b^2)}{b{y}_0}t-1=0---\left(25\right)$

式（25）的解实际上为如下矩阵的实特征值

$A=\left[\begin{array}{cccc}-\frac{2(a{x}_0-a^2+b^2)}{b{y}_0}& 0& -\frac{2(a{x}_0+a^2-b^2)}{b{y}_0}& -1\\ 1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\end{array}\right]---\left(26\right)$

由于椭圆上极值点为一个最大值和一个最小值，因此，矩阵A有且仅有两个实特征值t₁＝λ₁和t₂＝λ₂。因此，式（19）的解为

φ_i＝2tan^-1(t_i)    （27）

其中i为式（19）的实数解的个数。

步骤四：利用外切曲线参数计算实时TSE。

对于线切椭圆法，通过求该切点(x_t,y_t)到预计航迹的距离可以得到：

$\begin{array}{c}\mathrm{TSE}=\frac{|y_t-k{x}_t|}{\sqrt{1+k^2}}\\ =\frac{|y_0\&PlusMinus;\frac{b^2}{\sqrt{k^2{a}^2+b^2}}-k{x}_0\&PlusMinus;\frac{k^2{a}^2}{\sqrt{k^2{a}^2+b^2}}|}{\sqrt{1+k^2}}\\ =\frac{|y_0-k{x}_0\&PlusMinus;\sqrt{k^2{a}^2{b}^2}|}{\sqrt{1+k^2}}\end{array}---\left(28\right)$

可见线切椭圆法将得到两个解，为获得有效包络，需要计算外切直线，因此取其中较大的解作为结果。

对于圆切椭圆法最后，将获得的参数θ带入到式（18）中，为获得外切圆，因此选r(θ)较大值作为结果，即得

$\mathrm{TSE}=\underset{i}{\max }r\left({\mathrm{\&phi;}}_i\right)---\left(29\right)$

如图8所示，通过对1000个采样点进行试验，可得线切椭圆法，圆切椭圆法与传统标量求和法的TSE实时计算消耗结果，标量求和法和线切椭圆法的计算时间几乎为0，圆切椭圆法计算时间约为16毫秒，因此这些方法均可满足RNP实时性的需求。图9为三种方法TSE实时计算误差，由图可以看出，相比于标量求和法而言，线切椭圆法与圆切椭圆法均能能够对TSE进行的更精确的估计，其中，线切椭圆法效果最好。

步骤五：将实时TSE与RNP规范阈值x进行比较，输出告警结果。

当TSE＜x时，表明综合系统误差小于RNP规范阈值，此时正常飞行。反之，给出RNP告警信息。

本实例表明，通过采用本发明提出的RNP中的TSE实时计算方法，能够准确有效地对TSE进行实时估计，从而保障RNP飞行的安全性和效率。

Claims (1)

1.一种RNP中的综合系统误差实时计算方法，RNP表示所需导航性能，具体包括以下步骤：

步骤一：通过导航方程计算导航误差概率椭圆；

建立卫星导航非线性观测方程：

${\mathrm{\&rho;}}_{\mathrm{\&pi;}}=\left|\right|p_{\mathrm{\&pi;}}^{\mathrm{ENU}}-p_r^{\mathrm{ENU}}\left|\right|+C_b---\left(1\right)$

其中：ρ_π是第π颗可见卫星伪距，

是该卫星在ENU坐标系下的位置，ENU坐标系表示东北天坐标系，

是接收机在ENU坐标系下的位置，C_b是接收机钟差，||·||为欧式距离；通过对式（1）线性化，获得N_s颗可见卫星观测模型方程：

Δρ＝HΔx    （2）

其中：Δρ是N_s×1伪距残差矢量，H是N_s×4观测矩阵，Δx是状态误差矢量，包括三维位置[Δx,Δy,Δz]^T和接收机钟差ΔC_b，其中：Δx,Δy,Δz分别为东北天三个方向的位置误差，利用最小二乘法求解式（2），可解得：

Δx＝(H^TH)^-1H^TΔρ    （3）

假定

其中：

为方差，是大小为N_s×N_s的单位矩阵，Δρ～N(μ,Σ)表示Δρ服从均值为μ协方差矩阵为Σ的高斯分布，则Δx协方差为：

$\begin{array}{c}\mathrm{cov}\left(\mathrm{\&Delta;x}\right)=E<{\mathrm{\&Delta;x\&Delta;x}}^T>={\left(H^{T}H\right)}^{-1}{H}^T\mathrm{cov}\left(\mathrm{\&Delta;\&rho;}\right)H{\left(H^{T}H\right)}^{-1}\\ ={\left(H^{T}H\right)}^{-1}{H}^T{\mathrm{\&sigma;}}_p^2{I}_{N_s}H{\left(H^{T}H\right)}^{-1}{\mathrm{\&sigma;}}_p^2{\left(H^{T}H\right)}^{-1}\end{array}---\left(4\right)$

其中：cov(·)表示协方差矩阵；

将cov(Δx)写成展开形式，可得

$\mathrm{cov}\left(\mathrm{\&Delta;x}\right)=\left[\begin{array}{cccc}{\mathrm{\&sigma;}}_x^2& {\mathrm{\&sigma;}}_{\mathrm{yx}}^2& {\mathrm{\&sigma;}}_{\mathrm{zx}}^2& {\mathrm{\&sigma;}}_{C_{b}x}^2\\ {\mathrm{\&sigma;}}_{\mathrm{xy}}^2& {\mathrm{\&sigma;}}_y^2& {\mathrm{\&sigma;}}_{\mathrm{zy}}^2& {\mathrm{\&sigma;}}_{C_{b}y}^2\\ {\mathrm{\&sigma;}}_{\mathrm{xz}}^2& {\mathrm{\&sigma;}}_{\mathrm{yz}}^2& {\mathrm{\&sigma;}}_z^2& {\mathrm{\&sigma;}}_{C_{b}z}^2\\ {\mathrm{\&sigma;}}_{x{C}_b}^2& {\mathrm{\&sigma;}}_{y{C}_b}^2& {\mathrm{\&sigma;}}_{z{C}_b}^2& {\mathrm{\&sigma;}}_{C_b}^2\end{array}\right]---\left(5\right)$

其中：

表示误差α的方差，α＝x,y,z,C_b，

表示误差α与误差β的协方差，α,β＝x,y,z,C_b且α≠β；

由上可得状态误差适量Δx的东向和北向两个方向的分量x_EN的协方差 $\mathrm{\&Sigma;}=\left[\begin{array}{cc}{\mathrm{\&sigma;}}_x^2& {\mathrm{\&sigma;}}_{\mathrm{xy}}^2\\ {\mathrm{\&sigma;}}_{\mathrm{xy}}^2& {\mathrm{\&sigma;}}_y^2\end{array}\right],$ 即x_EN服从二维联合高斯分布：

$f\left(x_{\mathrm{EN}}\right)=\frac{1}{2\mathrm{\&pi;}\sqrt{\det \left(\mathrm{\&Sigma;}\right)}}\exp \{-\frac{1}{2}{(x_{\mathrm{EN}}-{\mathrm{\&mu;}}_{\mathrm{EN}})}^T{\mathrm{\&Sigma;}}^{-1}(x_{\mathrm{EN}}-{\mathrm{\&mu;}}_{\mathrm{EN}})\}---6$

其中，μ_EN为导航定位结果估计值的东向和北向两个方向的分量，det(·)为矩阵的行列式；

通过式（6）获取导航误差概率椭圆：

(x_EN-μ_EN)^TΣ^-1(x_EN-μ_EN)＝K²    （7）

其中：K为待定等概率椭圆常数，设定K=1.96；

步骤二：利用坐标系旋转算法将导航误差概率椭圆旋转为正椭圆；

将横向误差的协方差Σ进行对角化，可得：

$\mathrm{\&Sigma;}=\left[\begin{array}{cc}{\mathrm{\&sigma;}}_x^2& {\mathrm{\&sigma;}}_{\mathrm{xy}}^2\\ {\mathrm{\&sigma;}}_{\mathrm{xy}}^2& {\mathrm{\&sigma;}}_y^2\end{array}\right]={\mathrm{V\&Lambda;V}}^T---\left(8\right)$

其中：正交矩阵V满足VV^T＝I₂，I₂为2×2单位矩阵，记为 $V=\left[\begin{array}{cc}\cos \mathrm{\&theta;}& -\sin \mathrm{\&theta;}\\ \sin \mathrm{\&theta;}& \cos \mathrm{\&theta;}\end{array}\right],$ 其中θ为坐标系顺时针旋转角度，对角矩阵Λ＝diag(λ_a,λ_b)，其中λ_a和λ_b为矩阵Σ的两个特征值；

假设ENU坐标系经过顺时针旋转角度θ，得到新的坐标系，原坐标系下的向量x在新坐标系下的坐标x'为：

x＝Vx′    （9）

结合式（7），则原始误差椭圆在新坐标系下为正椭圆：

K²=(x_EN-μ_EN)^TΣ^-1(x_EN-μ_EN)=(x′_EN-μ′_EN)^TV^TΣ^-1V(x′_EN-μ′_EN)=(x′_EN-μ′_EN)^T(V^TΣV)^-1(x′_EN-μ′_EN)=(x′_EN-μ′_EN)^TΛ^-1(x′_ENμ′_EN)    (10)

综上所述，将原ENU坐标系下的任意平面坐标点x经线性变换可得新坐标系下的坐标点x'＝V^Tx，且在新坐标系下，导航定位误差椭圆为正椭圆，正椭圆的标准方程：

$\frac{{(x-x_0)}^2}{a^2}+\frac{{(y-y_0)}^2}{b^2}=1---\left(11\right)$

其中：误差椭圆的中心为经过坐标转化后的转变为μ'_EN＝(x₀,y₀)，半长轴

半短轴为 $b=K\sqrt{{\mathrm{\&lambda;}}_b};$

步骤三：求解正椭圆的外切曲线参数；

采用线切椭圆法或者圆切椭圆法获取，具体为：

A.线切椭圆法

线切椭圆法为采用平行于预计航迹的直线与椭圆最远点相切，获得椭圆上距离预计航迹最远的切线；根据椭圆参数，切点(x_t,y_t)满足落在椭圆线上的条件：

$\frac{{(x_t-x_0)}^2}{a^2}+\frac{{(y_t-y_0)}^2}{b^2}=1---\left(12\right)$

对式（12）的求全微分，得：

$\frac{x_t-x_0}{a^2}{\mathrm{dx}}_t+\frac{y_t-y_0}{b^2}{\mathrm{dy}}_t=0---\left(13\right)$

又切线与预计航迹平行，因此切线斜率满足：

$\frac{{\mathrm{dy}}_t}{{\mathrm{dx}}_t}=k---\left(14\right)$

i)当k＝0或∞时，椭圆的轴方向与预计航迹平行或垂直，由椭圆的性质得：

TSE＝b+y₀或TSE＝a+x₀    （15）

此时，线切椭圆法与标量求和法结果相同；

ii)当k≠0和∞时，联立求解式（12）-式（14）得到切点：

$\left\{\begin{array}{c}{y}_t=y_0\&PlusMinus;\frac{b^2}{\sqrt{k^2{a}^2+b^2}}\\ x_t=x_0\stackrel{\&OverBar;}{+}\frac{{\mathrm{ka}}^2}{\sqrt{k^2{a}^2+b^2}}\end{array}---\left(16\right)\right.$

B.圆切椭圆法

圆切椭圆法采用预计航迹点的最小半径圆O(r)对导航误差椭圆进行包络，该包络圆应过椭圆上距离预计航迹点最远的点(x_f,y_f)，即：

$\begin{array}{c}\underset{(x_f,y_f)}{\max }{r}^2=x_f^2+y_f^2\\ s.t.\frac{{(x_f-x_0)}^2}{a^2}+\frac{{(y_f-y_0)}^2}{b^2}=1\end{array}---\left(17\right)$

其中：r为圆O的半径，s.t.为约束条件；式（17）表明在(x_f,y_f)满足约束条件的条件下，求圆O的半径r的最大值；

将x_f与y_f表示成极坐标形式，得到x_f＝x₀+acosφ，y_f＝y₀+bsinφ，其中φ为极坐标参数，将其代入式（17）可转化得：

$\underset{\mathrm{\&phi;}}{\max }{r}^2\left(\mathrm{\&phi;}\right)={(x_0+a\cos )}^2+{(y_0+b\sin )}^2---\left(18\right)$

求r²(φ)相对于φ的导数，且满足极大值必要条件

即因此得到：

$\begin{array}{c}\frac{\&PartialD;r\left(\mathrm{\&phi;}\right)}{\&PartialD;\mathrm{\&phi;}}=-a(x_0+a\cos )\sin \mathrm{\&phi;}+b(y_0+b\sin \mathrm{\&phi;})\cos \mathrm{\&phi;}\\ =(b^2-a^2)\cos \mathrm{\&phi;}\sin \mathrm{\&phi;}-{\mathrm{ax}}_0\sin \mathrm{\&phi;}+{\mathrm{by}}_0\cos \mathrm{\&phi;}=0\end{array}---\left(19\right)$

i)当a＝b时，误差椭圆退化为正圆，此时式（19）可化简为：

-x₀sinφ+y₀cosφ＝0    （20）

将式（20）代入到式（18）中，可得：

$r\left(\mathrm{\&phi;}\right)=\sqrt{x_0^2+y_0^2}+a---\left(21\right)$

ii)当a≠b时，令 $t=\tan \left(\frac{\mathrm{\&phi;}}{2}\right),$ 则 $\sin \left(\mathrm{\&phi;}\right)=\frac{2t}{1+t^2},\cos \left(\mathrm{\&phi;}\right)=\frac{1-t^2}{1+t^2},$ 式（19）可转化为

$\frac{\&PartialD;r\left(\mathrm{\&phi;}\right)}{\&PartialD;\mathrm{\&phi;}}=(b^2-a^2)\frac{2t(1-t^2)}{{(1+t^2)}^2}-a{x}_0\frac{2t}{1+t^2}+b{y}_0\frac{1-t^2}{1+t^2}=0---\left(22\right)$

式（22）转化为一元四次方程

by₀t⁴+2(ax₀-a²+b²)t³+2(ax₀+a²-b²)t-by₀＝0    （23）

ii.i)当y₀＝0时，式（23）退化为一元三次方程

(ax₀-a²+b²)t³+(ax₀+a²-b²)t＝0    （24）

则式（23）的三个解分别为 $t_1=0,t_2=\sqrt{\frac{-a{x}_0-a^2+b^2}{a{x}_0-a^2+b^2}},t_3=-\sqrt{\frac{-a{x}_0-a^2+b^2}{a{x}_0-a^2+b_2}};$

ii.ii)当y₀≠0时，将式（23）首项系数化为1，可得：

$t^4+\frac{2(a{x}_0-a^2+b^2)}{b{y}_0}{t}^3+\frac{2(a{x}_0+a^2-b^2)}{b{y}_0}t-1=0---\left(25\right)$

式（25）的解实际上为如下矩阵的实特征值

$A=\left[\begin{array}{cccc}-\frac{2(a{x}_0-a^2+b^2)}{b{y}_0}& 0& -\frac{2(a{x}_0+a^2-b^2)}{b{y}_0}& -1\\ 1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\end{array}\right]---\left(26\right)$

由于椭圆上极值点为一个最大值和一个最小值，因此，矩阵A有且仅有两个实特征值t₁＝λ₁和t₂＝λ₂；因此，式（19）的解为

φ_i＝2tan^-1(t_i)    （27）

其中i为式（19）的实数解的个数；

步骤四：利用外切曲线参数计算实时综合系统误差；

A：对于线切椭圆法，通过求该切点(x_t,y_t)到预计航迹的距离可以得到：

$\begin{array}{c}\mathrm{TSE}=\frac{|y_t-k{x}_t|}{\sqrt{1+k^2}}\\ =\frac{|y_0\&PlusMinus;\frac{b^2}{\sqrt{k^2{a}^2+b^2}}-k{x}_0\&PlusMinus;\frac{k^2{a}^2}{\sqrt{k^2{a}^2+b^2}}|}{\sqrt{1+k^2}}\\ =\frac{|y_0-k{x}_0\&PlusMinus;\sqrt{k^2{a}^2{b}^2}|}{\sqrt{1+k^2}}\end{array}---\left(28\right)$

其中，TSE表示综合系统误差；

线切椭圆法将得到两个解，取其中较大的解作为结果；

B：对于圆切椭圆法，将获得的参数θ带入到式（18）中，为获得外切圆，因此选r(θ)较大值作为结果，即得

$\mathrm{TSE}=\underset{i}{\max }r\left({\mathrm{\&phi;}}_i\right)---\left(29\right)$

步骤五：将实时综合系统误差与RNP规范阈值x进行比较，当TSE＜x时，表明综合系统误差小于RNP规范阈值，此时正常飞行；反之，给出RNP告警信息。

Images