CN110045342A

CN110045342A - 雷达相对系统误差估值有效性评价方法

Info

Publication number: CN110045342A
Application number: CN201910378283.6A
Authority: CN
Inventors: 高效; 戴乐; 李光伟; 郑虹; 陈超; 张丽; 刘东红; 邢姗姗; 吴志宏; 马铭岩; 刘文华; 张明; 刘冰
Original assignee: Institute Of Strategic Early Warning Academy Of Air Force Chinese People's Liberation Army
Current assignee: Institute Of Strategic Early Warning Academy Of Air Force Chinese People's Liberation Army
Priority date: 2019-05-08
Filing date: 2019-05-08
Publication date: 2019-07-23
Anticipated expiration: 2039-05-08
Also published as: CN110045342B

Abstract

本发明属于多雷达数据融合技术领域，具体涉及一种雷达相对系统误差估值有效性评价方法。本发明通过选取主、次站雷达对某个典型航路目标的一组观测数据，经中心统一直角坐标转换后，先后采用不加权直线航迹线参数估计模型和加权直线航迹线参数估计模型对两雷达观测航迹线、次站雷达系统误差修正航迹线进行直线参数迭代估计，通过计算次站观测航迹线系统误差修正前后与主站观测航迹线的相对熵，构建次站修正航迹线系统误差估值有效性指标，用以评价系统误差估值的有效性以及修正效果。本发明原理方法科学、实施步骤合理，与传统系统误差有效性评价方法相比，在提高准确性和可操作性的同时，简化了复杂度，方便了工程上的实现。

Description

雷达相对系统误差估值有效性评价方法

技术领域

本发明属于多雷达数据融合技术领域，具体涉及一种雷达相对系统误差估值有效性评价方法。

背景技术

雷达相对系统误差是指雷达在对目标进行测量时存在的相对于指定参照系的固定误差，包括雷达动态测量中的测向误差Δθ、测距误差Δρ、雷达站址标定误差(Δx,Δy)和雷达组网应用中存在的坐标转换误差等。在雷达联网观测中，动态目标的绝对位置多数情况下是不可知的，比较容易得到的是在同一时间段内不同雷达对同一飞行目标的离散观测值。而相对于网内其它雷达来说，总会有一部雷达的观测是精确的，此时可以以该雷达(命名为主站)的测量值作为对目标位置的真实描述，其它雷达(命名为次站)以此为参照，从而求得次站雷达相对于主站雷达的相对系统误差。

多雷达组网应用中，不同雷达之间存在的相对系统误差会造成同一目标观测结果的空间分裂，严重时就会妨碍来自同一目标的航迹关联，而又造成对应不同目标的航迹的错误关联。因此研究系统误差的规律性，尽可能准确地定位、估计系统误差，对于提高雷达网目标状态估计和多雷达数据融合的准确性极为重要。

估计、消除雷达测量中存在的系统误差有两种途径：设备校准和数据校准。通常，雷达完成阵地部署后进行的入网前阵地检飞校准多采用设备校准手段。设备校准一般依靠固定地物回波或飞机检飞方式进行，在工程实施上有可能引入新的误差且效费比低，对多雷达组网后的系统误差校准，情况就更为复杂。因此，必须寻求科学合理、易于实现的方法评价一部新入网雷达经过阵地检飞校准后系统误差估值的有效性。

为了克服设备校准的局限性，人们提出了具有普遍意义的数据校准方法，如实时精度控制法、最小平方法、最大似然法、广义最小平方法和最小二乘相对系统误差修正法等等。它从数据处理的角度出发，在选定的参照系中对雷达测量数据进行分析，对其中存在的系统误差进行估计，并以此估计结果为依据，反过来对测量数据进行修正，以此达到同一坐标系内观测结果的一致性。这是一种直接、有效的校准方法，尤其在雷达组网后的误差修正中，显示出不可替代的优越性。但是如何评价某种数据校准方法的有效性，也成为一个亟待解决的问题。

与雷达测量中存在的随机误差相比，系统误差不随时间改变，是相对固定的误差。对于一部已经完成阵地部署的雷达而言，在一段时期内，相对于指定参照系的系统误差不会发生较大改变。但是，雷达在不同气候环境(春、夏、秋、冬、云、雨等)、不同探测区域(按照方位、距离分格)的系统误差可能存在差异，在一个系统误差估值不能“包打天下”的情况下，就要具体问题具体分析，“因地制宜”地建立雷达相对系统误差估值表。在这一过程中，研究雷达相对系统误差估值有效性评价方法就显得十分必要了。

传统的系统误差估值有效性评价多采用“取点评价法”，就是将同一时刻，雷达测量位置点与目标在参照系中的位置点(这两个点又称为时间配准点)之间的直线距离作为基本评价依据。为了提高评价的可信度，往往需要取多个时间配准点，将多个时间配准点修正前后距离的平均值之比作为最终的系统误差估值有效性评价指标。这种方法原理上比较科学，但在工程实现上，由于雷达异步工作，“多个时间配准点”在原始测量中几乎不存在，虽然采用卡尔曼滤波等方法可以外推得到心怡时刻的目标位置点，但外推过程又不可避免的引入了新的计算误差，因此取点评价法的可操作性和准确性不高。

发明内容

(一)要解决的技术问题

本发明要解决的技术问题是：如何提出一种雷达相对系统误差估值有效性评价方法，以提高雷达阵地检飞、数据校准方法选择和雷达相对系统误差估值表计算过程中对系统误差估值评价的科学性、准确性和可操作性，有助于提升雷达网目标状态估计和多雷达数据融合的质量。

(二)技术方案

为解决上述技术问题，本发明提供一种雷达相对系统误差估值有效性评价方法，所述评价方法包括如下步骤：

步骤1：选取主站雷达、次站雷达对同一空中目标的一段直线航迹线观测数据；

步骤2：对主站雷达观测数据进行坐标变换，得到一组对应的统一直角坐标{(X_tzi,Y_tzi)，i＝1,2,…n}，n表示主站雷达观测数据的总点数；

步骤3：分别对次站雷达观测数据、次站雷达系统误差修正数据进行坐标变换，得到两组统一直角坐标{(X_tcj,Y_tcj)，j＝1,2,…m}和{(X_txj,Y_txj)，j＝1,2,…m}，m表示次站雷达观测数据的总点数；

步骤4：在统一直角坐标系中，使用加权直线航迹线参数估计模型，分别对主站雷达观测数据、次站雷达观测数据和次站雷达系统误差修正数据进行直线参数迭代估计，得到主站观测航迹线参数(k_z,d_z)、次站观测航迹线参数(k_c,d_c)和次站修正航迹线参数(k_x,d_x)；所述k_z、k_c、k_x分别表示直线的斜率，d_z、d_c、d_x分别表示直线的截距；

步骤5：分别计算次站观测航迹线、次站修正航迹线与主站观测航迹线的相对熵JS_c和JS_x；相对熵越小，表明两条航迹线越相似；

步骤6：次站修正航迹线系统误差估值有效性指标V表示为：

(三)有益效果

与现有技术相比较，本发明通过选取主、次站雷达对某个典型航路目标的一组观测数据，经中心统一直角坐标转换后，先后采用不加权直线航迹线参数估计模型和加权直线航迹线参数估计模型对两雷达观测航迹线、次站雷达系统误差修正航迹线进行直线参数迭代估计，通过计算次站观测航迹线系统误差修正前后与主站观测航迹线的相对熵，构建次站修正航迹线系统误差估值有效性指标，用以评价系统误差估值的有效性以及修正效果。本发明原理方法科学、实施步骤合理，与传统系统误差有效性评价方法相比，在提高准确性和可操作性的同时，简化了复杂度，方便了工程上的实现。本发明有助于提高雷达阵地检飞、数据校准方法选择和雷达相对系统误差估值表计算过程中对系统误差估值评价的科学性、准确性和可操作性，从而提升雷达网目标状态估计和多雷达数据融合的质量。本发明所提供的方法的时间复杂度和空间复杂度都很低，可操作性和实用性很强。

附图说明

图1为本发明技术方案中评价方法的主体流程示意图。

图2为本发明实施例中主、次站雷达观测数据和两组次站雷达系统误差修正数据在统一直角坐标系中的显示图。

具体实施方式

为使本发明的目的、内容、和优点更加清楚，下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。

本发明针对传统方法从局部(单点)入手带来的问题与不足，引入信息论中相对熵的概念，通过选取典型航路(目标保持一定高度沿直线飞行)上的一组主、从雷达测量数据，整体比较次站观测航迹线系统误差修正前后与主站观测航迹线的相似性来评价系统误差估值的有效性以及修正效果。这样在提高准确性和可操作性的同时，简化了评价方法的复杂程度，方便了工程上的实现。

具体而言，为解决现有技术问题，本发明提供一种雷达相对系统误差估值有效性评价方法，所述方法应用于雷达阵地检飞、数据校准方法选择和雷达相对系统误差估值表计算过程中，属于雷达网目标状态估计和多雷达数据融合系统的前期数据准备工作范畴；

所述评价方法包括如下步骤：

步骤3：分别对次站雷达(原始)观测数据、次站雷达系统误差修正数据进行坐标变换，得到两组统一直角坐标{(X_tcj,Y_tcj)，j＝1,2,…m}和{(X_txj,Y_txj)，j＝1,2,…m}，m表示次站雷达观测数据的总点数；

步骤6：次站修正航迹线系统误差估值有效性指标V表示为：

V∈[0,1]时，V越大，表示系统误差估值越有效，航迹修正效果越好；V<0时，表示系统误差估值不仅完全无效，反而加剧了主、次站雷达观测航迹的分裂程度。

其中，所述步骤1包括如下步骤：

步骤1.1：选取空中目标处于直线航迹线时，主站雷达、次站雷达同时段上报的一段目标观测数据；

步骤1.2：选取的主站雷达观测数据为：{(ρ_zi,θ_zi,t_zi)，i＝1,2,…n}表示t_zi时刻主站雷达测得的目标距离ρ_zi和方位θ_zi，其中θ_zi∈[0,360),单位为度，雷达正北方向为0度，正东为90度，顺时针为正；n表示主站雷达观测数据的总点数；

步骤1.3：选取的次站雷达观测数据为：{(ρ_cj,θ_cj,t_cj)，j＝1,2,…m}，表示t_cj时刻次站雷达测得的目标距离ρ_cj和方位θ_cj，且|t_c1-t_z1|≤T，|t_cm-t_zn|≤T，m表示次站雷达观测数据的总点数。

其中，所述步骤1.1中，每部雷达观测数据点数量一般不少于10点。

其中，所述步骤1.1中，所述同时段是指：主站雷达观测数据、次站雷达观测数据的首点和末点时间差均不大于1个雷达探测周期T。

其中，所述雷达探测周期T一般为10或20秒。

其中，所述步骤2包括如下步骤：

步骤2.1：将主站雷达观测数据中的(ρ_zi,θ_zi)，i＝1,2,…n，转换为以本站为中心的二维直角坐标(X_zi,Y_zi)：

X_zi＝ρ_zisinθ_zi

Y_zi＝ρ_zicosθ_zi

步骤2.2：将(X_zi,Y_zi)，i＝1,2,…n转换为中心统一直角坐标(X_tzi,Y_tzi)：

X_tzi＝X_zicosδ_xz-Y_zisinδ_xz+X_zx

Y_tzi＝X_zisinδ_xz+Y_zicosδ_xz+Y_zx

其中：(X_zx,Y_zx)为主站雷达在中心统一直角坐标系中的坐标；δ_xz为主站雷达站址与直角坐标系中心点的经度差，单位为弧度。

其中，所述步骤3包括如下步骤：

步骤3.1：将次站雷达观测数据中的(ρ_cj,θ_cj)，j＝1,2,…m转换为以本站为中心的二维直角坐标(X_cj,Y_cj)：

X_cj＝ρ_cjsinθ_cj

Y_cj＝ρ_cjcosθ_cj

步骤3.2：将(X_cj,Y_cj)，j＝1,2,…m转换为中心统一直角坐标(X_tcj,Y_tcj)：

X_tcj＝X_cjcosδ_xc-Y_cjsinδ_xc+X_cx

Y_tcj＝X_cjsinδ_xc+Y_cjcosδ_xc+Y_cx

其中：(X_cx,Y_cx)为次站雷达在中心统一直角坐标系中的坐标；δ_xc为次站雷达站址与直角坐标系中心点的经度差，单位为弧度；

步骤3.3：将次站雷达系统误差修正数据(ρ_xj,θ_xj)，j＝1,2,…m转换为以本站为中心的二维直角坐标(X_xj,Y_xj)：

X_xj＝ρ_xjsinθ_xj

Y_xj＝ρ_xjcosθ_xj

此处若假设次站系统误差估值为(Δρ,Δθ,Δx,Δy)，表示次站雷达在测距、测向和站址定位上存在的相对于主站的系统误差，则有如下关系存在：

ρ_xj＝ρ_cj-Δρ

θ_xj＝θ_cj-Δθ

步骤3.4：将(X_xj,Y_xj)，j＝1,2,…m转换为中心统一直角坐标(X_txj,Y_txj)：

X_txj＝X_xjcosδ_xc-Y_xjsinδ_xc+X_cx-Δx

Y_txj＝X_xjsinδ_xc+Y_xjcosδ_xc+Y_cx-Δy

其中：(X_cx,Y_cx)为次站雷达在中心统一直角坐标系中的坐标；δ_xc为次站雷达站址与直角坐标系中心点的经度差，单位为弧度；(Δx,Δy)表示次站相对于主站的站址定位系统误差。

其中，所述步骤4包括如下步骤：

步骤4.1：采用不加权直线航迹线参数估计模型粗略估计主站雷达观测到的直线航迹线y-k₁x-d₁＝0，其中k₁为直线的斜率，d₁为直线在Y轴上的截距，步骤4.1包括如下步骤：

步骤4.1.1：用主站雷达所有观测点{(X_tzi,Y_tzi),i＝1,2,...n}到该直线的距离l_i的平方和最小作为条件构造直线，计算在此条件下的这条直线的最佳参数(k′,d′)，其中k′为待求直线的斜率，d′为待求直线的截距；此处，主站雷达所有观测点{(X_tzi,Y_tzi),i＝1,2,...n}简记为：{(x_i,y_i)，i＝1,2,…n}；

对于(1)式，应有f′(k′,d′)分别对k′和d′求偏导数，并等于零，即有下式成立：

令：a′＝c′₀-a′₁b′₁，b′＝a′₂-b′₂-a₁′²+b₁′²，c′＝a′₁b′₁-c′₀，则(2)式的解为：

步骤4.1.2：按照距离最小原则确定方程的合理解；

(k′₁,d′₁)和(k′₂,d′₂)都是方程(2)的实根，且k′₁×k′₂＝-1，即解得的两条直线相互垂直；按照观测点{(x_i,y_i)，i＝1,2,…n}到所求直线的距离的平方和最小原则，确定合理的直线参数值；该问题可简化为：计算测量点(x₁,y₁)分别到直线y＝k′₁×x+d′₁和直线y＝k′₂×x+d′₂的距离l′₁,l′₂：

若|l′₁|<|l′₂|，则取(k′₁,d′₁)，否则取(k′₂,d′₂)作为所求直线的合理参数，记为(k₁,d₁)；

步骤4.2：采用加权直线航迹线参数估计模型迭代估计主站雷达观测到的直线航迹线y-k_zx-d_z＝0，其中k_z为直线的斜率，d_z为直线的截距，步骤4.2包括如下步骤：

步骤4.2.1：计算主站雷达各观测点(x_i,y_i)到直线y-k₁x-d₁＝0的距离之和；

n为观测点数；

步骤4.2.2：求主站雷达各观测点(x_i,y_i)到直线y-k^(s)x-d^(s)＝0的距离l_i；

式中s表示迭代次数，n表示观测点数；s初始值为1，即：

k⁽¹⁾＝k₁，d⁽¹⁾＝d₁；

步骤4.2.3：求|l_i|的倒数；

步骤4.2.4：求各点的权值v_i；

步骤4.2.5：求解加权直线航迹线参数估计模型；

用主站雷达所有观测点{(x_i,y_i)，i＝1,2,…n}到该直线的加权距离(v_i×l_i)的平方和最小作为条件构造直线，计算在此条件下的这条直线的最佳参数(k,d)；

对于(3)式，应有f(k,d)分别对k和d求偏导数，并等于零，即有下式成立：

令：

c″＝c₀+a₁b₁，则(4)式的解为：

步骤4.2.6：按照距离最小原则确定方程的合理解；

按照点{(x_i,y_i)，i＝1,2,…n}到所求直线的距离的平方和最小原则，确定合理的直线参数值；该问题可简化为：计算点(x₁,y₁)分别到直线y＝k₁×x+d₁和直线y＝k₂×x+d₂的距离l₁,l₂：

若|l₁|+|l₂|>L min，L min初值为10⁶，则输出(k^(s-1),d^(s-1))作为所求直线的合理参数，并简记为(k_z,d_z)，迭代过程结束；否则Lmin＝|l₁|+|l₂|；

s值加1，若|l₁|<|l₂|，则取(k₁,d₁)，否则取(k₂,d₂)作为所求直线的合理参数，记为(k^(s),d^(s))，s表示迭代次数；

步骤4.2.7：计算所有观测点到新直线y-k^(s)x-d^(s)＝0的加权距离之和f^(s)(k^(s),d^(s))；

式中s表示迭代次数，n表示观测点数；

步骤4.2.8：判别是否为最佳解；

若f^(s)(k^(s),d^(s))≥f^(s-1)(k^(s-1),d^(s-1))，则输出解(k^(s-1),d^(s-1))，并简记为(k_z,d_z)；否则重复步骤4.2.2至步骤4.2.7；式中s表示迭代次数；

步骤4.3：参照步骤4.1的过程，采用不加权直线航迹线参数估计模型粗略估计次站雷达观测到的直线航迹线；

步骤4.4：参照步骤4.2的过程，采用加权直线航迹线参数估计模型迭代估计次站雷达观测到的直线航迹线；

步骤4.5：参照步骤4.1的过程，采用不加权直线航迹线参数估计模型粗略估计次站雷达数据进行系统误差修正后的直线航迹线y-k_x1x-d_x1＝0；

步骤4.6：参照步骤4.2的过程，采用加权直线航迹线参数估计模型迭代估计次站雷达数据进行系统误差修正后的直线航迹线y-k_xx-d_x＝0。

其中，所述步骤5包括如下步骤：

步骤5.1：计算次站观测航迹线与主站观测航迹线的相对熵JS_c；求解步骤如下:

①计算f(x_i′)和g(x_i′)。

f(x_i′)＝k_zx_i′+d_z，g(x_i′)＝k_cx_i′+d_c，i′＝1,2,...n+m。

其中：x_i′∈{(X_tzi|i＝1,2,...n)∪(X_tcj|j＝1,2,...m)}，表示主站观测点与次站观测点横坐标的并集，n和m分别表示主站雷达、次站雷达的观测点数；

②计算KL(f(x)||g(x))和KL(g(x)||f(x))；

③计算JS_c；

步骤5.2：计算次站修正航迹线与主站观测航迹线的相对熵JS_x；求解步骤如下：

①计算f′(x′_i′)和g′(x′_i′)。

f′(x′_i′)＝k_zx′_i′+d_z，g′(x′_i′)＝k_xx′_i′+d_x，i′＝1,2,...n+m。

其中：x′_i′∈{(X_tzi|i＝1,2,...n)∪(X_txj|j＝1,2,...m)}，表示主站观测点与次站修正点横坐标的并集，n和m分别表示主、次站雷达的观测点数；

②计算KL′(f′(x)||g′(x))和KL′(g′(x)||f′(x))；

③计算JS_x；

下面结合具体实施例来详细描述本发明。

实施例1

本实施例具体描述一种雷达相对系统误差估值有效性评价方法，并进行了模拟对比试验过程。所述评价方法应用于雷达阵地检飞、数据校准方法选择和雷达相对系统误差估值表计算过程中，属于雷达网目标状态估计和多雷达数据融合系统的前期数据准备工作范畴。

所述评价方法包括如下步骤：

步骤1：按照表1中的模拟雷达与目标参数，仿真产生主、次站雷达对同一空中目标的一段直线航迹线观测结果，见表2，表3。

表1：基本参数

表2：主站雷达观测数据

序号	时间t	距离ρ(km)	方位θ(度)
				1	0	249.92	179.96
2	10	247.25	179.53
				3	20	244.62	178.92
4	30	241.87	178.28
				5	40	239.06	177.85
6	50	236.60	177.31
				7	60	234.07	176.59
8	70	231.59	175.90
				9	80	229.01	175.54
10	90	226.55	174.78
				11	100	223.96	174.10
12	110	221.38	173.62
				13	120	218.97	172.76
14	130	216.63	172.24
				15	140	214.20	171.43

表3：次站雷达观测数据

步骤2：对表2的主站雷达观测数据进行坐标变换，得到统一直角坐标{(X_tzi,Y_tzi)，i＝1,2,…15}。

步骤2.1：将主站雷达观测数据点(ρ_zi,θ_zi)，i＝1,2,…15转换为以本站为中心的二维直角坐标(X_zi,Y_zi)，计算公式为：

X_zi＝ρ_zisinθ_zi

Y_zi＝ρ_zicosθ_zi

计算结果见表4。

表4：以主站为中心的二维直角坐标

序号	X<sub>z</sub>(km)	Y<sub>z</sub>(km)
			1	0.1803	-249.9159
2	2.0088	-247.2406
			3	4.6122	-244.5810
4	7.2612	-241.7625
			5	8.9662	-238.8894
6	11.0985	-236.3348
			7	13.9182	-233.6584
8	16.5745	-230.9913
			9	17.8050	-228.3177
10	20.5925	-225.6165
			11	23.0172	-222.7762
12	24.5882	-220.0104
			13	27.5795	-217.2243
14	29.2639	-214.6478
			15	31.9351	-211.8052

步骤2.2：将(X_zi,Y_zi)，i＝1,2,…11转换为中心统一直角坐标(X_xzi,Y_xzi)，计算公式为：

X_xzi＝X_zicosδ_xz-Y_zisinδ_xz+X_zx

Y_xzi＝X_zisinδ_xz+Y_zicosδ_xz+Y_zx

计算结果见表5。

表5：主站测量点的中心统一直角坐标

序号	X<sub>tz</sub>(km)	Y<sub>tz</sub>(km)
			1	1500.1803	2300.0841
2	1502.0088	2302.7594
			3	1504.6122	2305.4190
4	1507.2612	2308.2375
			5	1508.9662	2311.1106
6	1511.0985	2313.6652
			7	1513.9182	2316.3416
8	1516.5745	2319.0087
			9	1517.8050	2321.6823
10	1520.5925	2324.3835
			11	1523.0172	2327.2238
12	1524.5882	2329.9896
			13	1527.5795	2332.7757
14	1529.2639	2335.3522
			15	1531.9351	2338.1948

步骤3：分别对次站雷达原始观测数据、次站系统误差修正数据进行坐标变换，得到两组统一直角坐标{(X_tcj,Y_tcj)，j＝1,2,…15}和{(X_txj,Y_txj)，j＝1,2,…15}。

步骤3.1：将次站雷达观测数据点(ρ_cj,θ_cj)，j＝1,2,…15转换为以本站为中心的二维直角坐标(X_cj,Y_cj)，计算公式为：

X_cj＝ρ_cjsinθ_cj

Y_cj＝ρ_cjcosθ_cj

计算结果见表6。

表6：以次站为中心的二维直角坐标

步骤3.2：将(X_cj,Y_cj)，j＝1,2,…15转换为中心统一直角坐标(X_tcj,Y_tcj)，计算公式为：

X_tcj＝X_cjcosδ_xc-Y_cjsinδ_xc+X_cx

Y_tcj＝X_cjsinδ_xc+Y_cjcosδ_xc+Y_cx

计算结果见表7。

表7：次站测量点的中心统一直角坐标

序号	X<sub>tc</sub>(km)	Y<sub>tc</sub>(km)
			1	1492.9361	2292.3919
2	1497.8942	2294.5754
			3	1499.4064	2297.9326
4	1500.7656	2298.8805
			5	1504.3988	2304.2131
6	1508.1147	2306.1978
			7	1508.0098	2308.6039
8	1512.5714	2311.2676
			9	1514.0366	2313.9060
10	1517.1070	2316.1215
			11	1520.0867	2317.5894
12	1521.9967	2319.2396
			13	1525.3349	2324.1517
14	1525.6704	2324.5845
			15	1530.8372	2328.8352

步骤3.3：为了便于比较本发明的应用效果，我们使用表1中通过不同途径获得的需要评价的两组次站雷达系统误差估值1和估值2对表3中的次站雷达原始测量数据进行系统误差修正，得到次站雷达经测距、测向系统误差(Δρ,Δθ)修正后的数据集1和数据集2。见表8。其中，次站雷达系统误差修正数据1(ρ_1xj,θ_1xj)和数据2(ρ_2xj,θ_2xj)与次站雷达原始观测数据(ρ_cj,θ_cj)之间有如下关系：

ρ_1xj＝ρ_cj-Δρ₁，θ_1xj＝θ_cj-Δθ₁

ρ_2xj＝ρ_cj-Δρ₂，θ_2xj＝θ_cj-Δθ₂

表8：次站雷达系统误差修正数据集1和数据集2

然后，将表8中的次站雷达系统误差修正数据1和数据2转换为以本站为中心的二维直角坐标(X_1xj,Y_1xj)和(X_2xj,Y_2xj)：

X_1xj＝ρ_1xjsinθ_1xj，Y_1xj＝ρ_1xjcosθ_1xj

X_2xj＝ρ_2xjsinθ_2xj，Y_2xj＝ρ_2xjcosθ_2xj

计算结果见表9。

表9：次站雷达系统误差修正数据的站直角坐标

步骤3.4：将(X_1xj,Y_1xj)，(X_2xj,Y_2xj)，j＝1,2,…15转换为中心统一直角坐标(X_1txj,Y_1txj)和(X_2txj,Y_2txj)：

X_1txj＝X_1xjcosδ_xc-Y_1xjsinδ_xc+X_cx-Δx₁

Y_1txj＝X_1xjsinδ_xc+Y_1xjcosδ_xc+Y_cx-Δy₁

X_21txj＝X_2xjcosδ_xc-Y_2xjsinδ_xc+X_cx-Δx₂

Y_2txj＝X_2xjsinδ_xc+Y_2xjcosδ_xc+Y_cx-Δy₂

计算结果见表10。

表10：次站雷达系统误差修正数据的统一直角坐标

主、次站雷达原始测量点和两组次站雷达系统误差修正点在统一直角坐标系中的显示如图2所示。

步骤4：在统一直角坐标系中使用加权直线航迹线参数估计模型分别对主、次站雷达观测数据和次站雷达系统误差修正数据进行直线参数迭代估计，得到主站观测航迹线参数(k_z,d_z)、次站观测航迹线参数(k_c,d_c)和两组次站修正航迹线参数(k_1x,d_1x)、(k_2x,d_2x)。

步骤4.1：采用不加权直线航迹线参数估计模型粗略估计主站雷达观测到的直线航迹线y-k₁x-d₁＝0(k₁为直线的斜率，d₁为直线在Y轴上的截距)，包括如下步骤：

步骤4.1.1：用主站雷达所有观测点{(X_tzi,Y_tzi),i＝1,2,...15}(简记为：{(x_i,y_i)，i＝1,2,…15})到该直线的距离的平方和最小作为条件构造直线，计算在此条件下的这条直线的最佳参数(k,d)。对(1)式的求解步骤如下。

①计算a′₁,a′₂,b′₁,b′₂,c′₀(n＝15)。

②计算a′,b′,c′。

a′＝c′₀-a′₁b′₁＝115.0896，

b′＝a′₂-b′₂-a₁′²+b₁′²＝41.8806，

c′＝a′₁b′₁-c′₀＝-115.0896。

③解方程，计算所有解。

d′₁＝b′₁-a′₁k′₁＝502.4071，

d′₂＝b′₁-a′₁k′₂＝3584.1047。

步骤4.1.2：按照距离最小原则确定方程的合理解。

计算测量点(x₁,y₁)分别到直线y＝k′₁×x+d′₁和直线y＝k′₂×x+d′₂的距离l′₁,l′₂：

因|l′₁|<|l′₂|，则取(k′₁,d′₁)＝(1.1984，502.4071)作为所求直线的合理参数，记为(k₁,d₁)。

步骤4.2：采用加权直线航迹线参数估计模型迭代估计主站雷达观测到的直线航迹线y-k_zx-d_z＝0，包括如下步骤：

步骤4.2.1：计算主站各观测点(x_i,y_i)到直线y-k₁x-d₁＝0的距离之和(n＝15)。

步骤4.2.2：求各观测点(x_i,y_i)到直线y-k^(s)x-d^(s)＝0的距离l_i。计算结果如表11所示。

式中s表示迭代次数，初始值为1，即：k⁽¹⁾＝k₁，d⁽¹⁾＝d₁。

步骤4.2.3：求|l_i|的倒数。计算结果如表11所示。

步骤4.2.4：求各点的权值v_i。计算结果如表11所示。

表11距离和权值计算结果

序号	l<sub>i</sub>	P<sub>i</sub>	v<sub>i</sub>
				1	-0.056308	17.759487	0.122775
2	0.253868	3.939055	0.027232
				3	-0.041009	24.384793	0.168577
4	-0.269043	3.716876	0.025696
				5	0.262666	3.807117	0.026319
6	0.262270	3.812870	0.026359
				7	-0.187980	5.319726	0.036776
8	-0.518577	1.928353	0.013331
				9	0.249562	4.007015	0.027701
10	-0.159922	6.253035	0.043229
				11	-0.201866	4.953779	0.034247
12	0.363962	2.747540	0.018994
				13	-0.147662	6.772220	0.046818
14	0.209848	4.765347	0.032944
				15	-0.019809	50.483350	0.349002

步骤4.2.5：求解加权直线航迹线参数估计模型。

用主站雷达所有观测点{(x_i,y_i)，i＝1,2,…15}到该直线的加权距离(v_i×l_i)的平方和最小作为条件构造直线，计算在此条件下的这条直线的最佳参数(k,d)。对(3)式的求解步骤如下。

①计算a₀,a₁,a₂,b₁,b₂,c₀(n＝15)。

②计算a,b,c。

a＝-c₀-a₁b₁＝189.5367，

c″＝c₀+a₁b₁＝-189.5367。

③解方程，计算所有解。

d₁＝b₁-a₁k₁＝500.5643，

d₂＝b₁-a₁k₂＝3599.1113。

步骤4.2.6：按照距离最小原则确定方程的合理解。

计算测量点(x₁,y₁)分别到直线y＝k₁×x+d₁和直线y＝k₂×x+d₂的距离l₁,l₂：

因|l₁|+|l₂|＝37.192271<L min(初值为10⁶)，则L min＝37.192271。s加1后值为2，且|l₁|<|l₂|，则取(k₁,d₁)＝(1.1995，500.5643)作为所求直线的合理参数，记为(k⁽²⁾,d⁽²⁾)。

步骤4.2.7：计算所有观测点到新直线y-k⁽²⁾x-d⁽²⁾＝0的加权距离之和f⁽²⁾(k⁽²⁾,d⁽²⁾)。

步骤4.2.8：判别是否为最佳解。

因f⁽¹⁾(k⁽¹⁾,d⁽¹⁾)＝3.204352，f⁽²⁾(k⁽²⁾,d⁽²⁾)＝0.083881,f⁽²⁾(k⁽²⁾,d⁽²⁾)<f⁽¹⁾(k⁽¹⁾,d⁽¹⁾)，则重复步骤4.2.2至步骤4.2.7。

经计算，s＝3、执行步骤4.2.6时，|l₁|+|l₂|＝739.048264，大于上一次迭代时的Lmin＝37.192271，因此输出解(k⁽²⁾,d⁽²⁾)，并简记为(k_z,d_z)＝(1.1995，500.5643)。

步骤4.3：采用不加权直线航迹线参数估计模型粗略估计次站雷达观测到的直线航迹线y-k_c1x-d_c1＝0。具体过程与步骤4.1类同。

步骤4.4：采用加权直线航迹线参数估计模型迭代估计次站雷达观测到的直线航迹线y-k_cx-d_c＝0。具体过程与步骤4.2类同。计算结果为：(k_c,d_c)＝(1.004519，791.838008)。

步骤4.5：采用不加权直线航迹线参数估计模型粗略估计次站雷达数据采用系统误差估值1修正后的直线航迹线y-k_1x1x-d_1x1＝0。具体过程与步骤4.1类同。

步骤4.6：采用加权直线航迹线参数估计模型迭代估计次站雷达数据采用系统误差估值1修正后的直线航迹线y-k_1xx-d_1x＝0。具体过程与步骤4.2类同。计算结果为：(k_1x,d_1x)＝(1.153464，570.324616)。

步骤4.7：采用不加权直线航迹线参数估计模型粗略估计次站雷达数据采用系统误差估值2修正后的直线航迹线y-k_2x1x-d_2x1＝0。具体过程与步骤4.1类同。

步骤4.8：采用加权直线航迹线参数估计模型迭代估计次站雷达数据采用系统误差估值2修正后的直线航迹线y-k_2xx-d_2x＝0。具体过程与步骤4.2类同。计算结果为：(k_2x,d_2x)＝(1.095545，655.384635)。

步骤5：分别计算次站观测航迹线、两组次站修正航迹线与主站观测航迹线的相对熵JS_c、JS_1x和JS_2x。相对熵越小，表明两条航迹线越相似。所述步骤5包括如下步骤：

步骤5.1：计算次站观测航迹线与主站观测航迹线的相对熵JS_c。

①计算f(x_i′)和g(x_i′)。

f(x_i′)＝k_zx_i′+d_z，g(x_i′)＝k_cx_i′+d_c，i′＝1,2,...30；

x_i′∈{(X_tzi|i＝1,2,...15)∪(X_tcj|j＝1,2,...15)}

计算结果如表12所示。

表12 f(x_i')与g(x_i')计算结果

②计算KL(f(x)||g(x))和KL(g(x)||f(x))；

③计算JS_c。

步骤5.2：计算第一条次站修正航迹线与主站观测航迹线的相对熵JS_1x。

①计算f′(x′_i′)和g′(x′_i′)。

f′(x′_i′)＝k_zx′_i′+d_z，g′(x′_i′)＝k_xx′_i′+d_x，i′＝1,2,...30。

其中：x′_i′∈{(X_tzi|i＝1,2,...15)∪(X_1txj|j＝1,2,...15)}，表示主站观测点与次站第一组修正点横坐标的并集。

计算结果如表13所示。

表13 f′(x′_i′)与g′(x′_i′)计算结果

②计算KL′(f′(x)||g′(x))和KL′(g′(x)||f′(x))。

③计算JS_1x。

步骤5.3：计算第二条次站修正航迹线与主站观测航迹线的相对熵JS_2x。具体过程与步骤5.2类同。计算结果为：JS_2x＝0.084058。

由相对熵JS_c>JS_2x>JS_1x的结果可以看出，JS_1x＝0.002087对应的次站修正航迹线与主站观测航迹线更加相似(接近)，JS_2x＝0.084058对应的次站修正航迹线次之；相比之下，JS_c＝0.188534对应的次站观测航迹线与主站观测航迹线相似度最差，与图2的观察结果一致。

步骤6：则次站两组修正航迹线系统误差估值有效性指标V₁和V₂计算结果分别为：

由V₁>V₂>0可以得出以下结论：

第一，V₁和V₂对应的系统误差估值均能有效减小主、次站雷达观测航迹的分裂程度。

第二，在减小主、次站雷达观测航迹的分裂程度上，V₁对应的系统误差估值(Δρ₁,Δθ₁,Δx₁,Δy₁)＝(1.4,4.0,0.4,0.1)要优于V₂对应的系统误差估值(Δρ₂,Δθ₂,Δx₂,Δy₂)＝(0.5,2.5,-0.3,0.1)，取得了较好的航迹修正效果,与图2的观察结果一致。

本发明通过选取主、次站雷达对某个典型航路目标的一组观测数据，经中心统一直角坐标转换后，先后采用不加权直线航迹线参数估计模型和加权直线航迹线参数估计模型对两雷达观测航迹线、次站雷达系统误差修正航迹线进行直线参数迭代估计，通过计算次站观测航迹线系统误差修正前后与主站观测航迹线的相对熵，构建次站修正航迹线系统误差估值有效性指标，用以评价系统误差估值的有效性以及修正效果。本发明原理方法科学、实施步骤合理，与传统系统误差有效性评价方法相比，在提高准确性和可操作性的同时，简化了复杂度，方便了工程上的实现。本发明有助于提高雷达阵地检飞、数据校准方法选择和雷达相对系统误差估值表计算过程中对系统误差估值评价的科学性、准确性和可操作性，从而提升雷达网目标状态估计和多雷达数据融合的质量。本发明所提供的方法的时间复杂度和空间复杂度都很低，可操作性和实用性很强。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。例如但不限于以下几点：

(1)步骤3.3中引用的次站系统误差估值(Δρ,Δθ,Δx,Δy)仅仅是系统误差的某一种表述方式。系统误差在不同表述模型中具有多种样式，如增加时间误差、高度误差、测距系数误差等，在采用GPS等精确测量仪器后可以忽略定位误差(Δx,Δy)等。无论采取何种系统误差模型，都是用以将次站测量数据转换为修正数据，最终转换成统一直角坐标，进而转入步骤4。所以，系统误差模型的多样性不影响本发明的有效性，而应看成是本发明的某种变形应用。

(2)一定要确保步骤2和步骤3中将极坐标转换为中心统一直角坐标时，横、纵坐标值均大于零，必要时需要对中心统一直角坐标系原点进行平移。

(3)步骤5中，当航迹线与X轴夹角大于75度时，对于f(x_i)＝k_zx_i+d_z、g(x_j)＝k_cx_j+d_c和g(x_j)＝k_xx_j+d_x，要分别替换为：

和

(4)本发明所述实施步骤2和步骤3是针对两坐标雷达的表述，而对于三坐标雷达测量数据本发明同样适用。只要将相应的极坐标到雷达站二维直角坐标转换公式在考虑目标高度测量值的情况下稍作调整即可。

(5)本发明所述步骤1选取主、次站雷达极坐标值(ρ_zi,θ_zi,t_zi)作为观测数据适用于普遍情况。如果在雷达检飞或雷达站具备ADS-B接收设备的情况下，可以获得来自GPS或ADS-B设备更加准确、数据率更高的目标三维位置信息(通常以地理坐标形式给出)。此时只要将步骤2和步骤3中的极坐标到统一直角坐标转换公式替换为地理坐标到统一直角坐标的转换公式即可。

Claims

1.一种雷达相对系统误差估值有效性评价方法，其特征在于，所述评价方法包括如下步骤：

步骤6：次站修正航迹线系统误差估值有效性指标V表示为：

2.如权利要求1所述的雷达相对系统误差估值有效性评价方法，其特征在于，所述步骤1包括如下步骤：

3.如权利要求2所述的雷达相对系统误差估值有效性评价方法，其特征在于，所述步骤1.1中，每部雷达观测数据点数量不少于10点。

4.如权利要求2所述的雷达相对系统误差估值有效性评价方法，其特征在于，所述步骤1.1中，所述同时段是指：主站雷达观测数据、次站雷达观测数据的首点和末点时间差均不大于1个雷达探测周期T。

5.如权利要求4所述的雷达相对系统误差估值有效性评价方法，其特征在于，所述雷达探测周期T为10或20秒。

6.如权利要求2所述的雷达相对系统误差估值有效性评价方法，其特征在于，所述步骤2包括如下步骤：

X_zi＝ρ_zisinθ_zi

Y_zi＝ρ_zicosθ_zi

X_tzi＝X_zicosδ_xz-Y_zisinδ_xz+X_zx

Y_tzi＝X_zisinδ_xz+Y_zicosδ_xz+Y_zx

7.如权利要求6所述的雷达相对系统误差估值有效性评价方法，其特征在于，所述步骤3包括如下步骤：

X_cj＝ρ_cjsinθ_cj

Y_cj＝ρ_cjcosθ_cj

X_tcj＝X_cjcosδ_xc-Y_cjsinδ_xc+X_cx

Y_tcj＝X_cjsinδ_xc+Y_cjcosδ_xc+Y_cx

X_xj＝ρ_xjsinθ_xj

Y_xj＝ρ_xjcosθ_xj

ρ_xj＝ρ_cj-Δρ

θ_xj＝θ_cj-Δθ

X_txj＝X_xjcosδ_xc-Y_xjsinδ_xc+X_cx-Δx

Y_txj＝X_xjsinδ_xc+Y_xjcosδ_xc+Y_cx-Δy

8.如权利要求7所述的雷达相对系统误差估值有效性评价方法，其特征在于，所述步骤4包括如下步骤：

令：a′＝c′₀-a′₁b′₁，c′＝a₁′b₁′-c′₀，则(2)式的解为：

步骤4.1.2：按照距离最小原则确定方程的合理解；

n为观测点数；

式中s表示迭代次数，n表示观测点数；s初始值为1，即：k⁽¹⁾＝k₁，d⁽¹⁾＝d₁；

步骤4.2.3：求|l_i|的倒数；

步骤4.2.4：求各点的权值v_i；

步骤4.2.5：求解加权直线航迹线参数估计模型；

令：

c″＝c₀+a₁b₁，则(4)式的解为：

步骤4.2.6：按照距离最小原则确定方程的合理解；

若|l₁|+|l₂|>Lmin，Lmin初值为10⁶，则输出(k^(s-1),d^(s-1))作为所求直线的合理参数，并简记为(k_z,d_z)，迭代过程结束；否则Lmin＝|l₁|+|l₂|；

式中s表示迭代次数，n表示观测点数；

步骤4.2.8：判别是否为最佳解；

9.如权利要求8所述的雷达相对系统误差估值有效性评价方法，其特征在于，所述步骤5包括如下步骤：

①计算f(x_i′)和g(x_i′)。

f(x_i′)＝k_zx_i′+d_z，g(x_i′)＝k_cx_i′+d_c，i′＝1,2,...n+m。

②计算KL(f(x)||g(x))和KL(g(x)||f(x))；

③计算JS_c；

①计算f′(x′_i′)和g′(x′_i′)。

②计算KL′(f′(x)||g′(x))和KL′(g′(x)||f′(x))；

③计算JS_x；

10.如权利要求1所述的雷达相对系统误差估值有效性评价方法，其特征在于，所述步骤6中，V∈[0,1]时，V越大，表示系统误差估值越有效，航迹修正效果越好；V<0时，表示系统误差估值不仅完全无效，反而加剧了主、次站雷达观测航迹的分裂程度。