CN102508217A - 建立雷达测量误差标定模型的方法 - Google Patents

Abstract

本发明公开了一种建立雷达测量误差标定模型的方法，用于解决现有的雷达测量误差标定方法精度差的技术问题。技术方案是通过建立系统状态方程和观测方程，并引入系统偏差和测量偏差，建立完整的雷达测量误差标定模型，通过联合估计方法直接修正雷达的测量误差，减少积累误差的影响，提高了雷达测量误差标定方法的精度。

Description

建立雷达测量误差标定模型的方法

技术领域

本发明涉及一种建模方法，特别是涉及一种建立雷达测量误差标定模型的方法。 

背景技术

雷达测量误差是指目标参数测量值和真实值之间的偏差，它直接决定雷达的测量精度。对雷达测量误差的标定能够减少误差，改进雷达的测量精度，提高雷达检测和跟踪目标的能力，具有重要意义。 

雷达测量的主要参数为目标斜距、方位角、高低角、速度等。测量误差主要由热噪声、雷达自身不完善因素、电波传播因素、目标及其运动引起的误差等造成。主要包括：跟踪环路存在热噪声、伺服噪声和不平衡等而引起的雷达相关跟踪误差；雷达数据系统非线性、轴系不正交、结构变化等因素带来的雷达相关转换误差；目标回波的起伏和闪烁、目标的速度、加速度引起的目标相关的跟踪误差；以及由于对流层、电离层等导致雷达电波传播折射及不规则性带来的传播误差(王德纯等著，精密跟踪测量雷达技术.北京：电子工业出版社，2006)。雷达测量误差主要可分为两种类型：系统误差和随机误差。其中，随机误差可以通过各种滤波方法进行消除；系统误差是一种确定性误差，无法通过滤波方法去除，需要通过校准或标定技术进行修正。 

文献“胡波，梁星霞，练学辉.雷达系统误差的测量和修正方法.雷达与对抗，2005，2：12-15”公开了一种雷达测量误差标定方法，该标定方法 

1.静态标定：目的是为了确保雷达测量目标距离、方位的准确性。方法是：选择一个已知的孤立地标进行跟踪，录取数据。对数据进行处理，算出该地标的距离、方位，与给定的真值进行比较。若两者的均方差小于标定指标要求则完成静态标定；否则，根据两者一次差的分布，对雷达方位零位进行修正，然后再跟踪该地标，直到满足标定指标要求； 

2.动态标定：对三坐标雷达来说，还需要动态标定，其目的是保证雷达测量目标距离、方位角和高低角的准确性。方法为：跟踪飞机或挂有角反射器的气球，录取数据。计算这些数据与GPS或真值雷达录取数据的平均差和均方差，若满足标定指标要求则完成动态标定；否则，根据一次差的分布，对方位、距离和仰角零位进行分段修正，修正后再次跟踪，直到满足标定指标要求。 

3.系统误差的估值计算：标定试验后，对录取的数据进行处理、分析，可得到雷达的系统误差。通常采用分段方法来处理数据，在保证一定样本数的前提下，段内的数据应是平稳或近似平稳的随机函数。具体方法为：对每个航次，进入和远离均作为 一次进入计算每次进入数据与其对应真值之间的一次差，画出一次差相对于真值的曲线。汇总所有航次，根据一次差曲线的分布，将数据分段，根据下式计算每段数据的系统误差： 

${\mathrm{\Δx}}_p=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\Δx}}_i$

其中，Δx_p为每段数据的系统误差，Δx_i为测量值与真值之间的差值，N为数据段内的数据个数。 

上述雷达系统误差标定方法利用GPS或真值雷达测量数据作为参考值，与被标定雷达测量数据进行比较，计算参考值与测量值之间的差值，并将差值取平均值作为系统误差，对雷达系统误差进行修正。由于仅仅使用误差的平均值作为系统误差过于粗糙，误差修正精度较差。 

发明内容

为了克服现有的雷达测量误差标定方法精度差的不足，本发明提供一种建立雷达测量误差标定模型的方法。该方法通过建立系统状态方程和观测方程，并引入系统偏差和测量偏差，建立完整的雷达测量误差标定模型，通过联合估计方法直接修正雷达的测量误差，可以提高雷达测量误差标定方法的精度。 

本发明解决其技术问题所采用的技术方案是：一种建立雷达测量误差标定模型的方法，其特点是包括以下步骤： 

在原有的飞机运动方程中，引入飞机的斜距、方位角和高低角，得到 

式中，飞机飞行状态为平飞或盘旋，选择原点O和飞机质心重合而固连于飞机机体坐标系的动坐标系OXYZ，设u，v，w分别为质心沿动坐标系轴OX，OY，OZ的速度分量；雷达极坐标系原点位于雷达质心，R，A，B分别为在雷达极坐标系下的飞机斜距、高低角和方位角； 

为俯仰角， 

为滚转角，Ψ为偏航角，指向正北方向，p，q，r分别为飞机的滚转、俯仰、偏航角速率，n_x，n_y，n_z为过载分量，且 

n_x＝A_x/g，n_y＝A_y/g，n_z＝A_z/g 

式中，g为重力加速度，A_x，A_y，A_z为沿动坐标系轴OX，OY，OZ的加速度分量；设 

u_i＝[n_x n_y n_z p q r]^T

考虑到测量值中的系统偏差和尺度因子偏差，得 

$u_i=\left[\begin{array}{cccccc}I+{\mathrm{\λ}}_{\mathrm{nx}}& 0& 0& 0& 0& 0\\ 0& I+{\mathrm{\λ}}_{\mathrm{ny}}& 0& 0& 0& 0\\ 0& 0& {I+\mathrm{\λ}}_{\mathrm{nz}}& 0& 0& 0\\ 0& 0& 0& I+{\mathrm{\λ}}_p& 0& 0\\ 0& 0& 0& 0& I+{\mathrm{\λ}}_q& 0\\ 0& 0& 0& 0& 0& I+{\mathrm{\λ}}_r\end{array}\right]u_{\mathrm{im}}+\underset{\‾}{w}+b_{u_i}$

式中，u_im为u_i的测量值，w为过程噪声向量，即u_i的量测噪声， 

为状态向量的系统偏差向量，λ_nx、λ_ny、λ_nz、λ_p、λ_q和λ_r分别为n_x，n_y，n_z，p，q和r测量偏差的尺度因子；设 

$b_{u_i}={\left[\begin{array}{cccccc}{b}_{\mathrm{nx}}& b_{\mathrm{ny}}& b_{\mathrm{nz}}& b_p& b_q& b_r\end{array}\right]}^T$

式中，b_nx，b_ny，b_nz，b_p，b_q和b_r分别为n_x，n_y，n_z，p，q和r的系统偏差； 

飞机全面运动的观测方程为 

$V_0=\sqrt{u^2+v^2+w^2}$

$\mathrm{\α}={\tan }^{-1}\left(\frac{w-q{X}_{\mathrm{\α}}+p{Y}_{\mathrm{\α}}}{u}\right)$

$\mathrm{\β}={\tan }^{-1}\left(\frac{v+r{X}_{\mathrm{\β}}-p{Z}_{\mathrm{\β}}}{u}\right)$

$R=\sqrt{{111}^2{({\mathrm{\ζ}}_p-{\mathrm{\ζ}}_l)}^2+{111}^2{\cos }^2\mathrm{\δ}{({\mathrm{\σ}}_p-{\mathrm{\σ}}_l)}^2+h^2}$

${A=\tan }^{-1}\left(\frac{\cos \mathrm{\δ}({\mathrm{\σ}}_p-{\mathrm{\σ}}_l)}{{\mathrm{\ζ}}_p-{\mathrm{\ζ}}_l}\right)$

$B={\tan }^{-1}\left(\frac{h}{R}\right)$

式中，V₀为空速，X_α，Y_α，Z_α和X_β，Y_β，Z_β分别为气流迎角和侧滑角传感器在动坐标系中的坐标，α，β分别为气流迎角和侧滑角，h为高度，σ_l，ζ_l为雷达质心的经度和纬度坐标，σ_p，ζ_p为飞机的经度和纬度坐标， λ_α，λ_β，λ_h，λ_R，λ_A，λ_B分别为V₀，α，β，h，R，A，B测量偏差的尺度因子， 

b_α，b_β， 

b_Ψ，b_h，b_R，b_A，b_B和 

v _α，v _β， 

v _Ψ，v _h，v _R，v _A，v _B分别为V₀，α，β， 

Ψ，h，R，A，B的系统偏差和观测噪声，下标m代表测量值；V_0m，α_m，β_m，θ_m， 

Ψ_m，h_m，R_m，A_m，B_m分别为V_O，α，β，θ， 

Ψ，h，R，A，B的实际测量值； 

将观测方程写成矩阵形式，有 

y_m＝H_yy_c+b_y+v    (3) 

式中， 

为观测向量的计算值， 

为观测向量的实际测量值， 

$H_y=\mathrm{diag}[(1+{\mathrm{\λ}}_{V_0}),(1+{\mathrm{\λ}}_{\mathrm{\α}}),(1+{\mathrm{\λ}}_{\mathrm{\β}}),\mathrm{1,1,1},(I+{\mathrm{\λ}}_h),(I+{\mathrm{\λ}}_R),(I+{\mathrm{\λ}}_A),(I+{\mathrm{\λ}}_B)]$ 为观测向量的测量偏差矩阵， 

为观测向量的系统偏差向量， 

为观测噪声向量， 

取 $b_1={\left[\begin{array}{ccccccc}{\mathrm{\λ}}_{\mathrm{nx}}& {\mathrm{\λ}}_{\mathrm{ny}}& {\mathrm{\λ}}_{\mathrm{nz}}& {\mathrm{\λ}}_p& {\mathrm{\λ}}_q& {\mathrm{\λ}}_r& b_{u_i}^T\end{array}\right]}^T$

$b_2={\left[\begin{array}{cccccccc}{\mathrm{\λ}}_{V_0}& {\mathrm{\λ}}_{\mathrm{\α}}& {\mathrm{\λ}}_{\mathrm{\β}}& {\mathrm{\λ}}_h& {\mathrm{\λ}}_R& {\mathrm{\λ}}_A& {\mathrm{\λ}}_B& b_y^T\end{array}\right]}^T$

根据(1)式，系统的状态方程写为 

式中，x为9维状态向量，f为9维函数向量，b₁为12维状态偏差向量，Γ(x)为过程噪声的系数矩阵，w为零均值高斯白噪声序列，方差为Q； 

观测方程写为 

y＝к(x，u_i，b₂)+v    (5) 

式中， 

为10维观测向量，b₂为17维观测偏差向量，к为10维输出函数向量，v为零均值高斯白噪声序列，方差为R； 

将(4)式和(5)式按泰勒级数展开，并取一次项，得雷达测量误差的标定模型为 

$\left\{\begin{array}{c}x(k+1)=\underset{\‾}{A}\left(k\right)x\left(k\right)+\underset{\‾}{B}\left(k\right)b_1\left(k\right)+R_u\left(k\right)+\mathrm{\Γ}\left(k\right)\underset{\‾}{w}\left(k\right)\\ y(k+1)=H(k+1)x(k+1)+D(k+1)b_2(k+1)+R_y(k+1)+\underset{\‾}{v}(k+1)\end{array}\right.---\left(6\right)$

式中，x(k)为目标kT时刻的状态向量， 

b₁(k)为kT时刻的状态偏差向量， 

Γ(k)为kT时刻过程噪声的系数矩阵， 

w(k)为kT时刻的过程噪声向量， 

y(k)为kT时刻的观测向量， 

b₂(k)为kT时刻的观测偏差向量， 

v(k)为kT时刻的观测噪声向量， 

$\underset{\‾}{A}\left(k\right)=\frac{\∂f[x\left(k\right),b_1,k]}{\∂x}{|}_{b_1=b_1\left(k\right),x\left(k\right)=x\left(k\right|k)}$

$\underset{\‾}{B}\left(k\right)=\frac{\∂f[x\left(k\right),b_1,k]}{\∂b_1}{|}_{b_1=b_1\left(k\right),x\left(k\right)=x\left(k\right|k)}$

$H(k+1)=\frac{\∂\mathrm{\κ}[x(k+1),b_2,k+1]}{\∂x}{|}_{b_2=b_2\left(k\right),x\left(k\right)=x(k+1|k)}$

$D(k+1)=\frac{\∂\mathrm{\κ}[x(k+1),b_2,k+1]}{\∂b_2}{|}_{b_2=b_2\left(k\right),x\left(k\right)=x(k+1|k)}$

R_u(k)＝f[x(k|k)，b₁(k)，k]-A(k)x(k|k)-B(k)b₁(k) 

R_y(k)＝к[x(k+1/k)，b₂(k)，k+1]-H(k+1)x(k+1|k)-D(k+1)b₂(k) 

$\frac{\∂\mathrm{\κ}[x(k+1),b_2,k+1]}{\∂x}=\left[\begin{array}{ccccccccc}(1+{\mathrm{\λ}}_{V_0})\frac{u}{V_0}& (1+{\mathrm{\λ}}_{V_0})\frac{v}{V_0}& (1+{\mathrm{\λ}}_{V_0})\frac{w}{V_0}& 0& 0& 0& 0& 0& 0\\ h_{21}& 0& h_{23}& 0& 0& 0& 0& 0& 0\\ h_{31}& h_{32}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& (1+{\mathrm{\λ}}_R)& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& (1+{\mathrm{\λ}}_A)& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& (1+{\mathrm{\λ}}_B)\end{array}\right]$

$h_{21}=\frac{-(1+{\mathrm{\λ}}_{\mathrm{\α}})(w-q{X}_{\mathrm{\α}}+p{Y}_{\mathrm{\α}})}{u^2(1+{\tan }^2\mathrm{\α})}$

$h_{23}=\frac{(1+{\mathrm{\λ}}_{\mathrm{\α}})}{u(1+{\tan }^2\mathrm{\α})}$

$h_{31}=\frac{-(1+{\mathrm{\λ}}_{\mathrm{\β}})(v+r{X}_{\mathrm{\β}}-p{Z}_{\mathrm{\β}})}{u^2(1+{\tan }^2\mathrm{\β})}$

$h_{32}=\frac{(1+{\mathrm{\λ}}_{\mathrm{\β}})}{u(1+{\tan }^2\mathrm{\β})}$

x(k/k)为目标kT时刻状态的滤波值， 

x(k/k-1)为目标kT时刻状态的一步预测值， 

u_i(k)为kT时刻的输入， 

n_xm，n_ym，n_zm，p_m，q_m，r_m分别为n_x，n_y，n_z，p，q，r的实际测量值， 

$\frac{\∂\mathrm{\κ}[x(k+1),b_2,k+1]}{{\∂b}_2}=\left[\begin{array}{ccccccccccccccccc}{V}_0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& \mathrm{\α}& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& \mathrm{\β}& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& h& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& R& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& A& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& B& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]$

$f_{71}=\cos B\cos A\cos \mathrm{\θ}\cos \mathrm{\ψ}+\cos B\sin A\cos \mathrm{\θ}\sin \mathrm{\ψ}-\sin B\sin \mathrm{\θ}$

$f_{81}=-\sin B\cos A\cos \mathrm{\θ}\cos \mathrm{\ψ}-\sin A\sin B\cos \mathrm{\θ}\sin \mathrm{\ψ}-u\sin \mathrm{\θ}({\cos }^{2}A\cos B+\sin A\cos A\cos B)$

p_c＝(I+λ_p)p_m+b_p

q_c＝(I+λ_q)q_m+b_q

r_c＝(I+λ_r)r_m+b_r

p_c，q_c，r_c为p，q，r的计算值。 

本发明的有益效果是：由于通过建立系统状态方程和观测方程，并引入系统偏差和测量偏差，建立完整的雷达测量误差标定模型，通过联合估计方法直接修正雷达的测量误差，减少积累误差的影响，提高了雷达测量误差标定方法的精度。 

下面结合具体实施方式对本发明作详细说明。 

具体实施方式

在原有的飞机运动方程中，引入飞机的斜距、方位角和高低角，得到： 

其中，假设飞机飞行状态为平飞或盘旋，选择原点O和飞机质心重合而固连于飞机机体坐标系的动坐标系(广义的体轴系)OXYZ，设u，v，w分别为质心沿动坐标系轴OX，OY，OZ的速度分量；雷达极坐标系原点位于雷达质心，R，A，B分别为在雷达极坐 标系下的飞机斜距、高低角和方位角； 为俯仰角， 

为滚转角，Ψ为偏航角，指向正北方向，p，q，r分别为飞机的滚转、俯仰、偏航角速率，n_x，n_y，n_z为过载分量，且 

n_x＝A_x/g，n_y＝A_y/g，n_z＝A_z/g 

g为重力加速度，A_x，A_y，A_z为沿动坐标系轴OX，OY，OZ的加速度分量；设 

u_i＝[n_x n_y n_z p q r]^T

输入向量u_i的准确值无法知道，只能得到其测量值，考虑到测量值中的系统偏差和尺度因子偏差，可得 

$u_i=\left[\begin{array}{cccccc}I+{\mathrm{\λ}}_{\mathrm{nx}}& 0& 0& 0& 0& 0\\ 0& I+{\mathrm{\λ}}_{\mathrm{ny}}& 0& 0& 0& 0\\ 0& 0& {I+\mathrm{\λ}}_{\mathrm{nz}}& 0& 0& 0\\ 0& 0& 0& I+{\mathrm{\λ}}_p& 0& 0\\ 0& 0& 0& 0& I+{\mathrm{\λ}}_q& 0\\ 0& 0& 0& 0& 0& I+{\mathrm{\λ}}_r\end{array}\right]u_{\mathrm{im}}+\underset{\‾}{w}+b_{u_i}$

其中，u_im为u_i的测量值，w为过程噪声向量，即u_i的量测噪声， 

为状态向量的系统偏差向量，λ_nx、λ_ny、λ_nz、λ_p、λ_q和λ_r分别为n_x，n_y，n_z，p，q和r测量偏差的尺度因子；设 

$b_{u_i}={\left[\begin{array}{cccccc}{b}_{\mathrm{nx}}& b_{\mathrm{ny}}& b_{\mathrm{nz}}& b_p& b_q& b_r\end{array}\right]}^T$

其中，b_nx，b_ny，b_nz，b_p，b_q和b_r分别为n_x，n_y，n_z，p，q和r的系统偏差； 

飞机全面运动的观测方程为 

$V_0=\sqrt{u^2+v^2+w^2}$

$\mathrm{\α}={\tan }^{-1}\left(\frac{w-q{X}_{\mathrm{\α}}+p{Y}_{\mathrm{\α}}}{u}\right)$

$\mathrm{\β}={\tan }^{-1}\left(\frac{v+r{X}_{\mathrm{\β}}-p{Z}_{\mathrm{\β}}}{u}\right)$

$R=\sqrt{{111}^2{({\mathrm{\ζ}}_p-{\mathrm{\ζ}}_l)}^2+{111}^2{\cos }^2\mathrm{\δ}{({\mathrm{\σ}}_p-{\mathrm{\σ}}_l)}^2+h^2}$

${A=\tan }^{-1}\left(\frac{\cos \mathrm{\δ}({\mathrm{\σ}}_p-{\mathrm{\σ}}_l)}{{\mathrm{\ζ}}_p-{\mathrm{\ζ}}_l}\right)$

$B={\tan }^{-1}\left(\frac{h}{R}\right)$

其中，V₀为空速，X_α，Y_α，Z_α和X_β，Y_β，Z_β分别为气流迎角和侧滑角传感器在动坐标系中的坐标，α，β分别为气流迎角和侧滑角，h为高度，σ_l，ζ_l为雷达质心(NEU坐标系原点)的经度和纬度坐标，σ_p，ζ_p为飞机的经度和纬度坐标， 

λ_α，λ_β，λ_h，λ_R，λ_A，λ_B分别为V₀，α，β，h，R，A，B测量偏差的尺度因子， 

b_α，b_β， 

b_Ψ，b_h，b_R，b_A，b_B和  v _α，v _β， 

v _Ψ，v _h，v _R，v _A，v _B分别为V₀，α，β， 

Ψ，h，R，A，B的系统偏差和观测噪声，下标m代表测量值；V_0m，α_m，β_m，θ_m， 

Ψ_m，h_m，R_m，A_m，B_m分别为V_O，α，β，θ， 

Ψ，h，R，A，B的实际测量值； 

将观测方程写成矩阵形式，有 

y_m＝H_yy_c+b_y+v    (3) 

其中 

为观测向量的计算值， 

为观测向量的实际测量值， 

$H_y=\mathrm{diag}[(1+{\mathrm{\λ}}_{V_0}),(1+{\mathrm{\λ}}_{\mathrm{\α}}),(1+{\mathrm{\λ}}_{\mathrm{\β}}),\mathrm{1,1,1},(I+{\mathrm{\λ}}_h),(I+{\mathrm{\λ}}_R),(I+{\mathrm{\λ}}_A),(I+{\mathrm{\λ}}_B)]$ 为观测向量的测量偏差矩阵， 

为观测向量的系统偏差向量， 

为观测噪声向量， 

取 $b_1={\left[\begin{array}{ccccccc}{\mathrm{\λ}}_{\mathrm{nx}}& {\mathrm{\λ}}_{\mathrm{ny}}& {\mathrm{\λ}}_{\mathrm{nz}}& {\mathrm{\λ}}_p& {\mathrm{\λ}}_q& {\mathrm{\λ}}_r& b_{u_i}^T\end{array}\right]}^T$

$b_2={\left[\begin{array}{cccccccc}{\mathrm{\λ}}_{V_0}& {\mathrm{\λ}}_{\mathrm{\α}}& {\mathrm{\λ}}_{\mathrm{\β}}& {\mathrm{\λ}}_h& {\mathrm{\λ}}_R& {\mathrm{\λ}}_A& {\mathrm{\λ}}_B& b_y^T\end{array}\right]}^T$

根据(1)式，系统的状态方程可写成 

其中，x为9维状态向量，f为9维函数向量，b₁为12维状态偏差向量，Γ(x)为过程噪声的系数矩阵，w为零均值高斯白噪声序列，方差为Q； 

观测方程可写成 

y＝к(x，u_i，b₂)+v    (5) 

为10维观测向量，b₂为17维观测偏差向量，к为10维输出函数向量，v为零均值高斯白噪声序列，方差为R； 

将(4)式和(5)式按泰勒级数展开，并取一次项，可得雷达测量误差的标定模型为 

$\left\{\begin{array}{c}x(k+1)=\underset{\‾}{A}\left(k\right)x\left(k\right)+\underset{\‾}{B}\left(k\right)b_1\left(k\right)+R_u\left(k\right)+\mathrm{\Γ}\left(k\right)\underset{\‾}{w}\left(k\right)\\ y(k+1)=H(k+1)x(k+1)+D(k+1)b_2(k+1)+R_y(k+1)+\underset{\‾}{v}(k+1)\end{array}\right.---\left(6\right)$

式中，x(k)为目标kT 时刻的状态向量， 

b₁(k)为kT时刻的状态偏差向量， 

Γ(k)为kT时刻过程噪声的系数矩阵， 

w(k)为kT时刻的过程噪声向量， 

y(k)为kT时刻的观测向量， 

b₂(k)为kT时刻的观测偏差向量， 

v(k)为kT时刻的观测噪声向量， 

$\underset{\‾}{A}\left(k\right)=\frac{\∂f[x\left(k\right),b_1,k]}{\∂x}{|}_{b_1=b_1\left(k\right),x\left(k\right)=x\left(k\right|k)}$

$\underset{\‾}{B}\left(k\right)=\frac{\∂f[x\left(k\right),b_1,k]}{\∂b_1}{|}_{b_1=b_1\left(k\right),x\left(k\right)=x\left(k\right|k)}$

$H(k+1)=\frac{\∂\mathrm{\κ}[x(k+1),b_2,k+1]}{\∂x}{|}_{b_2=b_2\left(k\right),x\left(k\right)=x(k+1|k)}$

$D(k+1)=\frac{\∂\mathrm{\κ}[x(k+1),b_2,k+1]}{\∂b_2}{|}_{b_2=b_2\left(k\right),x\left(k\right)=x(k+1|k)}$

R_u(k)＝f[x(k|k)，b₁(k)，k]-A(k)x(k|k)-B(k)b₁(k) 

R_y(k)＝к[x(k+1/k)，b₂(k)，k+1]-H(k+1)x(k+1|k)-D(k+1)b₂(k) 

$\frac{\∂\mathrm{\κ}[x(k+1),b_2,k+1]}{\∂x}=\left[\begin{array}{ccccccccc}(1+{\mathrm{\λ}}_{V_0})\frac{u}{V_0}& (1+{\mathrm{\λ}}_{V_0})\frac{v}{V_0}& (1+{\mathrm{\λ}}_{V_0})\frac{w}{V_0}& 0& 0& 0& 0& 0& 0\\ h_{21}& 0& h_{23}& 0& 0& 0& 0& 0& 0\\ h_{31}& h_{32}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& (1+{\mathrm{\λ}}_R)& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& (1+{\mathrm{\λ}}_A)& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& (1+{\mathrm{\λ}}_B)\end{array}\right]$

$h_{21}=\frac{-(1+{\mathrm{\λ}}_{\mathrm{\α}})(w-q{X}_{\mathrm{\α}}+p{Y}_{\mathrm{\α}})}{u^2(1+{\tan }^2\mathrm{\α})}$

$h_{23}=\frac{(1+{\mathrm{\λ}}_{\mathrm{\α}})}{u(1+{\tan }^2\mathrm{\α})}$

$h_{31}=\frac{-(1+{\mathrm{\λ}}_{\mathrm{\β}})(v+r{X}_{\mathrm{\β}}-p{Z}_{\mathrm{\β}})}{u^2(1+{\tan }^2\mathrm{\β})}$

$h_{32}=\frac{(1+{\mathrm{\λ}}_{\mathrm{\β}})}{u(1+{\tan }^2\mathrm{\β})}$

x(k/k)为目标kT时刻状态的滤波值， 

x(k/k-1)为目标kT时刻状态的一步预测值， 

u_i(k)为kT时刻的输入， 

n_xm，n_ym，n_zm，p_m，q_m，r_m分别为n_x，n_y，n_z，p，q，r的实际测量值， 

$\frac{\∂\mathrm{\κ}[x(k+1),b_2,k+1]}{{\∂b}_2}=\left[\begin{array}{ccccccccccccccccc}{V}_0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& \mathrm{\α}& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& \mathrm{\β}& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& h& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& R& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& A& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& B& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]$

$f_{71}=\cos B\cos A\cos \mathrm{\θ}\cos \mathrm{\ψ}+\cos B\sin A\cos \mathrm{\θ}\sin \mathrm{\ψ}-\sin B\sin \mathrm{\θ}$

$f_{81}=-\sin B\cos A\cos \mathrm{\θ}\cos \mathrm{\ψ}-\sin A\sin B\cos \mathrm{\θ}\sin \mathrm{\ψ}-u\sin \mathrm{\θ}({\cos }^{2}A\cos B+\sin A\cos A\cos B)$

p_c＝(I+λ_p)p_m+b_p

q_c＝(I+λ_q)q_m+b_q

r_c＝(I+λ_r)r_m+b_r

p_c，q_c，r_c为p，q，r的计算值； 

建立上述雷达测量误差的标定模型后，作为实例，采用非线性系统状态估计和偏差辨识的分离算法(史忠科.最优估计的计算方法[M].北京：科学出版社，2001)对雷达测量误差进行修正： 

时间更新 

$\left.\begin{array}{c}x(k+1/k)=f[x\left(k\right|k),b_1\left(k\right),k]\\ P_x^{\%}(k+1|k)=A\left(k\right)P_x^{\%}\left(k\right|k)A^T\left(k\right)+\mathrm{\Γ}\left(k\right)Q\left(k\right){\mathrm{\Γ}}^T\left(k\right)\end{array}\right\}$

式中， 

表示0偏差估计误差的协方差阵， 

f为函数向量， 

x(k/k)为目标kT时刻状态的滤波值， 

x(k/k-1)为目标kT时刻状态的一步预测值， 

b₁(k)为kT时刻的状态偏差向量， 

Γ(x)为过程噪声的系数矩阵， 

Q(k)为过程噪声向量w的方差； 

状态测量更新为 

$\left.\begin{array}{c}x(k+1|k+1)=f[x\left(k\right|k),b_1\left(k\right),k]+\{G_x^{\%}(k+1)+V_b(k+1)G_b(k+1)\}\\ \×\{y(k+1)-\mathrm{\κ}[x(k+1|k),b_2\left(k\right),k+1\left]\right\}\\ G_x^{\%}(k+1)=P_x^{\%}(k+1|k+1)H^T(k+1)R^{-1}(k+1)\\ P_x^{\%}(k+1/k+1)=[1-G_x^{\%}(k+1)H(k+1)]P_x^{\%}(k+1|k)\end{array}\right\}---(7.a)$

式中，b₂(k)为kT时刻的观测偏差向量， 

R(k)为观测噪声向量v的方差， 

y(k+1)为观测向量， 

к为输出函数向量， 

偏差辨识的测量更新为 

$\left.\begin{array}{c}{b}_2(k+1)=b_2\left(k\right)+G_b(k+1)\{y(k+1)-\mathrm{\κ}[x(k+1|k),b_2\left(k\right),k+1\left]\right\}\\ G_b(k+1)=P_0(k+1|k+1){[H(k+1)V_b(k+1)+D(k+1)]}^T{R}^{-1}(k+1)\\ P_b(k+1|k+1)=P_0\left(k\right|k)+C^T(k+1)[H(k+1)P_x^{\%}(k+1|k)H^T(k+1)\\ +R(k+1){]}^{-1}C(k+1)\end{array}\right\}---(7.b)$

式中 

$\left.\begin{array}{c}C(k+1)=H(k+1)U\left(k\right)+D(k+1)\\ V_b(k+1)=U\left(k\right)-G_x^{\%}(k+1)C(k+1)\\ U\left(k\right)=A\left(k\right)V_b\left(k\right)+B\left(k\right)\end{array}\right\}---(7.c)$

初始条件为P₀(k|k)； 

式(7.a)～式(7.c)就构成了完整的非线性离散系统的分离算法。 

Claims (1)

1.一种建立雷达测量误差标定模型的方法，其特征在于包括以下步骤：

在原有的飞机运动方程中，引入飞机的斜距、方位角和高低角，得到

式中，飞机飞行状态为平飞或盘旋，选择原点O和飞机质心重合而固连于飞机机体坐标系的动坐标系OXYZ，设u，v，w分别为质心沿动坐标系轴OX，OY，OZ的速度分量；雷达极坐标系原点位于雷达质心，R，A，B分别为在雷达极坐标系下的飞机斜距、高低角和方位角； 

为俯仰角， 为滚转角，Ψ为偏航角，指向正北方向，p，q，r分别为飞机的滚转、俯仰、偏航角速率，n_x，n_y，n_z为过载分量，且

n_x＝A_x/g，n_y＝A_y/g，n_z＝A_z/g

式中，g为重力加速度，A_x，A_y，A_z为沿动坐标系轴OX，OY，OZ的加速度分量；设

u_i＝[n_x n_y n_z p q r]^T

考虑到测量值中的系统偏差和尺度因子偏差，得 

式中，u_im为u_i的测量值，w为过程噪声向量，即u_i的量测噪声， 

为状态向量的系统偏差向量，λ_nx、λ_ny、λ_nz、λ_p、λ_q和λ_r分别为n_x，n_y，n_z，p，q和r测量偏差的尺度因子；设

式中，b_nx，b_ny，b_nz，b_p，b_q和b_r分别为n_x，n_y，n_z，p，q和r的系统偏差；

飞机全面运动的观测方程为

式中，V₀为空速，X_α，Y_α，Z_α和X_β，Y_β，Z_β分别为气流迎角和侧滑角传感器在动坐标系中的坐标，α，β分别为气流迎角和侧滑角，h为高度，σ_l，ζ_l为雷达质心的经度和纬度坐标，σ_p，ζ_p为飞机的经度和纬度坐标， λ_α，λ_β，λ_h，λ_R，λ_A，λ_B分别为V₀，α，β，h，R，A，B测量偏差的尺度因子， b_α，b_β， 

b_Ψ，b_h，b_R，b_A，b_B和 

v _α，v _β， 

v _Ψ，v _h，v _R，v _A，v _B分别为V₀，α，β， 

Ψ，h，R，A，B的系统偏差和观测噪声，下标m代表测量值；V_0m，α_m，β_m，θ_m， 

Ψ_m，h_m，R_m，A_m，B_m分别为V_O，α，β，θ， 

Ψ，h，R，A，B的实际测量值；

将观测方程写成矩阵形式，有

y_m＝H_yy_c+b_y+v    (3)

式中，

为观测向量的计算值，

为观测向量的实际测量值，

为观测向量的测量偏差矩阵，

为观测向量的系统偏差向量，

为观测噪声向量，

取

根据(1)式，系统的状态方程写为

式中，x为9维状态向量，f为9维函数向量，b₁为12维状态偏差向量，Γ(x)为过程噪声的系数矩阵，w为零均值高斯白噪声序列，方差为Q；

观测方程写为

y＝к(x，u_i，b₂)+v    (5)

式中， 

为10维观测向量，b₂为17维观测偏差向量，к为10维输出函数向量，v为零均值高斯白噪声序列，方差为R； 

将(4)式和(5)式按泰勒级数展开，并取一次项，得雷达测量误差的标定模型为

式中，x(k)为目标kT时刻的状态向量，

b₁(k)为kT时刻的状态偏差向量，

Γ(k)为kT时刻过程噪声的系数矩阵，

w(k)为kT时刻的过程噪声向量，

y(k)为kT时刻的观测向量，

b₂(k)为kT时刻的观测偏差向量，

v(k)为kT时刻的观测噪声向量，

R_u(k)＝f[x(k|k)，b₁(k)，k]-A(k)x(k|k)-B(k)b₁(k)

R_y(k)＝к[x(k+1/k)，b₂(k)，k+1]-H(k+1)x(k+1|k)-D(k+1)b₂(k)

x(k/k)为目标kT时刻状态的滤波值，

x(k/k-1)为目标kT时刻状态的一步预测值，

u_i(k)为kT时刻的输入，

n_xm，n_ym，n_zm，p_m，q_m，r_m分别为n_x，n_y，n_z，p，q，r的实际测量值，

p_c＝(I+λ_p)p_m+b_p

q_c＝(I+λ_q)q_m+b_q

r_c＝(I+λ_r)r_m+b_r

p_c，q_c，r_c为p，q，r的计算值。