[发明专利]一种RNP中的综合系统误差实时计算方法

权利要求书

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

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

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

ρπ=||pπENU-prENU||+Cb---(1)]]>

其中：ρ_π是第π颗可见卫星伪距，是该卫星在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协方差为：

cov(Δx)=E<ΔxΔxT>=(HTH)-1HTcov(Δρ)H(HTH)-1=(HTH)-1HTσp2INsH(HTH)-1σp2(HTH)-1---(4)]]>

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

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

cov(Δx)=σx2σyx2σzx2σCbx2σxy2σy2σzy2σCby2σxz2σyz2σz2σCbz2σxCb2σyCb2σzCb2σCb2---(5)]]>

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

由上可得状态误差适量Δx的东向和北向两个方向的分量x_EN的协方差Σ=σx2σxy2σxy2σy2,]]>即x_EN服从二维联合高斯分布：

f(xEN)=12πdet(Σ)exp{-12(xEN-μEN)TΣ-1(xEN-μEN)}---6]]>

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

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

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

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

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

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

Σ=σx2σxy2σxy2σy2=VΛVT---(8)]]>

其中：正交矩阵V满足VV^T＝I₂，I₂为2×2单位矩阵，记为V=cosθ-sinθsinθcosθ,]]>其中θ为坐标系顺时针旋转角度，对角矩阵Λ＝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，且在新坐标系下，导航定位误差椭圆为正椭圆，正椭圆的标准方程：

(x-x0)2a2+(y-y0)2b2=1---(11)]]>

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

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

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

A.线切椭圆法

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

(xt-x0)2a2+(yt-y0)2b2=1---(12)]]>

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

xt-x0a2dxt+yt-y0b2dyt=0---(13)]]>

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

dytdxt=k---(14)]]>

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

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

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

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

yt=y0±b2k2a2+b2xt=x0+‾ka2k2a2+b2---(16)]]>

B.圆切椭圆法

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

max(xf,yf)r2=xf2+yf2s.t.(xf-x0)2a2+(yf-y0)2b2=1---(17)]]>

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

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

maxφr2(φ)=(x0+acos)2+(y0+bsin)2---(18)]]>

求r²(φ)相对于φ的导数，且满足极大值必要条件即因此得到：

∂r(φ)∂φ=-a(x0+acos)sinφ+b(y0+bsinφ)cosφ=(b2-a2)cosφsinφ-ax0sinφ+by0cosφ=0---(19)]]>

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

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

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

r(φ)=x02+y02+a---(21)]]>

ii)当a≠b时，令t=tan(φ2),]]>则sin(φ)=2t1+t2,cos(φ)=1-t21+t2,]]>式（19）可转化为

∂r(φ)∂φ=(b2-a2)2t(1-t2)(1+t2)2-ax02t1+t2+by01-t21+t2=0---(22)]]>

式（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）的三个解分别为t1=0,t2=-ax0-a2+b2ax0-a2+b2,t3=--ax0-a2+b2ax0-a2+b2;]]>

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

t4+2(ax0-a2+b2)by0t3+2(ax0+a2-b2)by0t-1=0---(25)]]>

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

A=-2(ax0-a2+b2)by00-2(ax0+a2-b2)by0-1100001000010---(26)]]>

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

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

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

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

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

TSE=|yt-kxt|1+k2=|y0±b2k2a2+b2-kx0±k2a2k2a2+b2|1+k2=|y0-kx0±k2a2b2|1+k2---(28)]]>

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

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

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

TSE=maxir(φi)---(29)]]>

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