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

摘要

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

主权项

1.一种RNP中的综合系统误差实时计算方法，RNP表示所需导航性能，具体包括以下步骤：步骤一：通过导航方程计算导航误差概率椭圆；建立卫星导航非线性观测方程：${\rho }_{\pi }=\left|\right|p_{\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{\sigma }_p^2{I}_{N_s}H{\left(H^{T}H\right)}^{-1}{\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}{\sigma }_x^2& {\sigma }_{\mathrm{yx}}^2& {\sigma }_{\mathrm{zx}}^2& {\sigma }_{C_{b}x}^2\\ {\sigma }_{\mathrm{xy}}^2& {\sigma }_y^2& {\sigma }_{\mathrm{zy}}^2& {\sigma }_{C_{b}y}^2\\ {\sigma }_{\mathrm{xz}}^2& {\sigma }_{\mathrm{yz}}^2& {\sigma }_z^2& {\sigma }_{C_{b}z}^2\\ {\sigma }_{x{C}_b}^2& {\sigma }_{y{C}_b}^2& {\sigma }_{z{C}_b}^2& {\sigma }_{C_b}^2\end{array}\right]---\left(5\right)$其中：表示误差α的方差，α＝x,y,z,C_b，表示误差α与误差β的协方差，α,β＝x,y,z,C_b且α≠β；由上可得状态误差适量Δx的东向和北向两个方向的分量x_EN的协方差$\Sigma =\left[\begin{array}{cc}{\sigma }_x^2& {\sigma }_{\mathrm{xy}}^2\\ {\sigma }_{\mathrm{xy}}^2& {\sigma }_y^2\end{array}\right],$即x_EN服从二维联合高斯分布：$f\left(x_{\mathrm{EN}}\right)=\frac{1}{2\pi \sqrt{\det \left(\Sigma \right)}}\exp \{-\frac{1}{2}{(x_{\mathrm{EN}}-{\mu }_{\mathrm{EN}})}^T{\Sigma }^{-1}(x_{\mathrm{EN}}-{\mu }_{\mathrm{EN}})\}---6$其中，μ_EN为导航定位结果估计值的东向和北向两个方向的分量，det(·)为矩阵的行列式；通过式（6）获取导航误差概率椭圆：(x_EN-μ_EN)^TΣ^-1(x_EN-μ_EN)＝K²    （7）其中：K为待定等概率椭圆常数，设定K=1.96；步骤二：利用坐标系旋转算法将导航误差概率椭圆旋转为正椭圆；将横向误差的协方差Σ进行对角化，可得：$\Sigma =\left[\begin{array}{cc}{\sigma }_x^2& {\sigma }_{\mathrm{xy}}^2\\ {\sigma }_{\mathrm{xy}}^2& {\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 \theta & -\sin \theta \\ \sin \theta & \cos \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{{\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\pm \frac{b^2}{\sqrt{k^2{a}^2+b^2}}\\ x_t=x_0\overline{+}\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{\phi }{\max }{r}^2\left(\phi \right)={(x_0+a\cos )}^2+{(y_0+b\sin )}^2---\left(18\right)$求r²(φ)相对于φ的导数，且满足极大值必要条件即因此得到：$\begin{array}{c}\frac{\partial r\left(\phi \right)}{\partial \phi }=-a(x_0+a\cos )\sin \phi +b(y_0+b\sin \phi )\cos \phi \\ =(b^2-a^2)\cos \phi \sin \phi -{\mathrm{ax}}_0\sin \phi +{\mathrm{by}}_0\cos \phi =0\end{array}---\left(19\right)$i)当a＝b时，误差椭圆退化为正圆，此时式（19）可化简为：-x₀sinφ+y₀cosφ＝0    （20）将式（20）代入到式（18）中，可得：$r\left(\phi \right)=\sqrt{x_0^2+y_0^2}+a---\left(21\right)$ii)当a≠b时，令$t=\tan \left(\frac{\phi }{2}\right),$则$\sin \left(\phi \right)=\frac{2t}{1+t^2},\cos \left(\phi \right)=\frac{1-t^2}{1+t^2},$式（19）可转化为$\frac{\partial r\left(\phi \right)}{\partial \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\pm \frac{b^2}{\sqrt{k^2{a}^2+b^2}}-k{x}_0\pm \frac{k^2{a}^2}{\sqrt{k^2{a}^2+b^2}}|}{\sqrt{1+k^2}}\\ =\frac{|y_0-k{x}_0\pm \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({\phi }_i\right)---\left(29\right)$步骤五：将实时综合系统误差与RNP规范阈值x进行比较，当TSE＜x时，表明综合系统误差小于RNP规范阈值，此时正常飞行；反之，给出RNP告警信息。