CN113743475A - 一种基于ukf的实时多源数据融合方法 - Google Patents
一种基于ukf的实时多源数据融合方法 Download PDFInfo
- Publication number
- CN113743475A CN113743475A CN202110912112.4A CN202110912112A CN113743475A CN 113743475 A CN113743475 A CN 113743475A CN 202110912112 A CN202110912112 A CN 202110912112A CN 113743475 A CN113743475 A CN 113743475A
- Authority
- CN
- China
- Prior art keywords
- data
- time
- measured
- matrix
- variable
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000007500 overflow downdraw method Methods 0.000 title claims abstract description 14
- 239000011159 matrix material Substances 0.000 claims abstract description 54
- 238000012360 testing method Methods 0.000 claims abstract description 46
- 238000000034 method Methods 0.000 claims abstract description 37
- 238000007781 pre-processing Methods 0.000 claims abstract description 31
- 238000001914 filtration Methods 0.000 claims abstract description 10
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 7
- 238000005259 measurement Methods 0.000 claims description 68
- 230000002159 abnormal effect Effects 0.000 claims description 51
- 238000012937 correction Methods 0.000 claims description 51
- 238000004364 calculation method Methods 0.000 claims description 25
- 238000005070 sampling Methods 0.000 claims description 22
- 230000014509 gene expression Effects 0.000 claims description 19
- 238000006243 chemical reaction Methods 0.000 claims description 18
- 238000012545 processing Methods 0.000 claims description 17
- 230000005540 biological transmission Effects 0.000 claims description 11
- 238000013213 extrapolation Methods 0.000 claims description 9
- 230000001133 acceleration Effects 0.000 claims description 5
- 239000000126 substance Substances 0.000 claims description 3
- 238000012546 transfer Methods 0.000 claims description 3
- 230000001960 triggered effect Effects 0.000 claims description 3
- 230000004927 fusion Effects 0.000 abstract description 17
- 230000000295 complement effect Effects 0.000 abstract description 2
- 238000004088 simulation Methods 0.000 description 4
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000007499 fusion processing Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000002203 pretreatment Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000002035 prolonged effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/25—Fusion techniques
- G06F18/251—Fusion techniques of input or preprocessed data
Landscapes
- Engineering & Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Artificial Intelligence (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明提出了一种基于UKF的实时多源数据融合方法，其步骤为：首先，数据的接收：通过数据接收线程接收测元；其次，数据预处理：根据测元类型的不同，分别采用不同的预处理方法对不同类型的测元进行处理，并将预处理后的测元存放在数据缓存区中；最后，数据解算：利用测元选择策略从数据缓存区中查找当前时刻可用测元，生成测元矩阵，并通过UKF滤波算法对测元矩阵进行求解，得到并显示目标运动参数。本发明将各测量设备的直接观测数据进行数据融合，利用各个测量设备之间的互补关系，提高测量设备利用率，实现对直接观测数据的高效利用，减少数据冗余同时通过测量设备的合理部署，延长测量弧段并提高测量精度。
Description
技术领域
本发明涉及航天测控和靶场测量技术领域，特别是指一种基于UKF的实时多源数据融合方法。
背景技术
实时多源数据融合在航天测控和靶场测量系统中是数据处理中的重要一步。各个测站测量设备在跟踪目标的同时生成观测数据，并通过网络传输实时传送给数据处理软件实现多源数据融合，生成目标的实时运动参数。
目前，数据处理中多源数据融合大多是在运动轨迹级别的融合，即对每个/类测量设备的观测数据单独进行处理解算得到运动参数，再对各设备得到的运动参数进行融合处理，而非从各测量设备的直接观测数据进行融合。运动轨迹级别的融合是在单一设备的测量信息下进行的，不能有效利用同一时间段其他测量设备的观测数据，从而造成同一时间段内数据源之间的信息损失。另一方面，由于各测量设备的测量原理不同和实际跟踪情况等客观因素造成的直接观测数据的测量精度等先验信息也不能直接利用运动轨迹级别的融合。
发明内容
针对上述背景技术中存在的不足，本发明提出了一种基于UKF的实时多源数据融合方法，解决了由于各测量设备的测量原理不同和实际跟踪情况等客观因素造成的直接观测数据的测量精度等先验信息也不能直接利用运动轨迹级别的融合的技术问题。
本发明的技术方案是这样实现的：
一种基于UKF的实时多源数据融合方法，其步骤如下：
步骤一：数据的接收：通过数据接收线程接收测元，其中，测元包括测速元素、测距元素和测角元素；
步骤二：数据预处理：根据测元类型的不同，分别采用不同的预处理方法对不同类型的测元进行处理，并将预处理后的测元存放在数据缓存区中；
步骤三：数据解算：利用测元选择策略从数据缓存区中查找当前时刻可用测元，生成测元矩阵，并通过UKF滤波算法对测元矩阵进行求解，得到并显示目标运动参数。
测速元素的预处理包括量纲转换、异常值修正、时间误差修正和电波折射修正；
量纲转换：将设备输出的二进制数据恢复成十进制物理量，并将直接物理量多普勒频率信息转换成距离和变化率；
异常值修正：利用数据的连续性判断观测量中的数据是否异常，采用四点外推拟合法对观测量中的异常值进行处理；
时间误差修正：采样时刻tk测量设备接收到来自目标tk'时刻转发的电波信号，电波信号在空间传输时延为
式中，R(tk')为tk'时刻目标与观测设备的径向距离，c为光速；将tk'时刻的观测量向后修正Δtk，得到tk时刻对应的观测量；
电波折射修正：采用H.S.Hopfield折射指数模型将电磁波在经过不均匀大气层时产生传播路径弯曲折射进行修正，缩小实际测量值与理论观测值的差异。
测距元素的预处理包括异常值修正、时间误差修正和电波折射修正；
异常值修正：利用数据的连续性判断观测量中的数据是否异常，采用四点外推拟合法对观测量中的异常值进行处理；
时间误差修正：采样时刻tk测量设备接收到来自目标tk'时刻转发的电波信号，电波信号在空间传输时延为
式中，R(tk')为tk'时刻目标与观测设备的径向距离，c为光速；将tk'时刻的观测量向后修正Δtk，得到tk时刻对应的观测量；
电波折射修正：采用H.S.Hopfield折射指数模型将电磁波在经过不均匀大气层时产生传播路径弯曲折射进行修正，缩小实际测量值与理论观测值的差异。
测角元素的预处理包括方位角和俯仰角的预处理；方位角预处理包括角坐标转换、方位角跨零跳处理、异常值修正和时间误差修正；俯仰角预处理包括角坐标转换、异常值修正、时间误差修正和电波折射修正；
角坐标转换：方位角和俯仰角的测量是在垂线测量坐标系下测得的，将测得的方位角和俯仰角转换成统一坐标系下；
方位角跨零跳处理：方位角跨零跳是指目标在I象限和IV象限之间移动时，直接观测量方位角会出现0°和360°之间的跳跃，导致方位角在数值上的不连续，此时对方位角进行差分计算，得到差分序列，如果差分序列中有结果大于180°，则对应位置的方位值加360°，如果差分结果小于-180°，则对应位置的方位值减360°，从而进行方位角跨零跳处理，使得方位角在数值上连续有效；
异常值修正：利用数据的连续性判断观测量中的数据是否异常，采用四点外推拟合法对观测量中的异常值进行处理；
时间误差修正：采样时刻tk测量设备接收到来自目标tk'时刻转发的电波信号，电波信号在空间传输时延为
式中，R(tk')为tk'时刻目标与观测设备的径向距离，c为光速；将tk'时刻的观测量向后修正Δtk，得到tk时刻对应的观测量；
电波折射修正：采用H.S.Hopfield折射指数模型将电磁波在经过不均匀大气层时产生传播路径弯曲折射进行修正，缩小实际测量值与理论观测值的差异。
所述利用测元选择策略从数据缓存区中查找当前时刻可用测元，生成测元矩阵的方法为：
当解算命令触发后，从数据缓存区中进行测元选择，测元选择具体步骤如下：
a)以当前解算时刻Tj为基准，在数据缓存区中查找该时刻的可用测元；
b)判断查找的可用测元是否为异常值，如果是，则可用测元不写入测元矩阵；如果不是，则可用测元写入测元矩阵；并重复步骤a)直至遍历数据缓存区；
c)判断测元矩阵中的非零测元个数，如果非零测元个数小于6个，则该测元矩阵直接参与解算；如果非零测元个数大于等于6个，则从最高优先级的测元开始由高到低对非零测元个数进行判断，如果当前优先级及之前的测元数量大于等于6个，则将低于该优先级的测元数据置零，否则继续遍历优先级较低的测元数据，直到有效测元个数大于等于6个或者遍历整个测元矩阵为止。
所述通过UKF滤波算法对测元矩阵进行求解的方法为：
I)U变换取样
第k个测量时刻的原始状态变量为：
原始状态方差为：
式中：
原始状态变量xk为n维列向量；其中，xk表示地心直角坐标系下X方向的位置，
表示地心直角坐标系下X方向的速度，
表示地心直角坐标系下X方向的加速度，zk表示地心直角坐标系下Z方向的位置，
表示地心直角坐标系下Z方向的速度，
表示地心直角坐标系下Z方向的加速度；
对原始状态变量xk进行U变换得到变量Xk，Xk为n×(2n+1)维矩阵，Xk具有如下表达式：
式中：Pk为xk的协方差矩阵；
系数λ＝α2(n+κ)-n，α和κ均为待选参数；
II)UKF的时间更新
变量Xk的状态随时间的更新表达式为式(1)，变量Xk对应的加权采样均值
的表达式为式(2)，变量Xk的协方差矩阵表达式为式(3)，变量Xk对应的采样观测量Zi,k的表达式为式(4)，采样观测量Zi,k对应的加权采样均值
的表达式为式(5)：
Xi,k＝F(Xi,k-1,wi,k-1) (1)
Zi,k＝H(Xi,k,vi,k) (4)
其中，Xi,k表示k时刻变量Xk的第i列，Xi,k-1表示k-1时刻变量Xk的第i列，wi,k-1表示k-1时刻变量状态量噪声矩阵wk-1的第i列，F(·)表示状态转移函数，Wi表示变量Xi,k的加权系数，Pk'表示变量Xk的协方差矩阵，Wi c表示变量Xi,k均值的加权系数，vi,k表示k时刻观测噪声矩阵的第i列，H(·)表示观测函数，i＝0,1,2,...,2n表示i的取值范围为0～2n的正整数，其中n表示原始状态变量xk的维数；
III)UKF的测量更新
采样观测量Zi,k的协方差矩阵表达式为式(6)，变量Xk与采样观测量Zi,k的互协方差表达式为式(7)，则卡尔曼增益Kk的表达式为式(8)，状态量及状态协方差更新表达式分别为式(9)和(10)：
其中，
为采样观测量Zi,k的协方差矩阵，
为变量Xk与采样观测量Zi,k的互协方差，
为
的逆矩阵，
为卡尔曼增益Kk的转置。
与现有技术相比，本发明产生的有益效果为：本发明将各测量设备的直接观测数据进行数据融合，利用各个测量设备之间的互补关系，提高测量设备利用率，实现对直接观测数据的高效利用，减少数据冗余同时通过测量设备的合理部署，延长测量弧段并提高测量精度。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明的实时多源数据融合流程图。
图2为本发明的不同测元预处理的流程图。
图3为本发明的测元选择策略的实现流程图。
图4为本发明模拟的目标轨迹地面投影与测量设备的位置分布图。
图5为本发明模拟解算的结果比对；其中，(a)为连续波测速设备单独解算结果比对，(b)为单脉冲设备单独解算结果比对，(c)为相控阵设备单独解算结果比对，(d)为三种设备联合解算结果比对。
图6为实时数据处理软件的参数设置界面。
图7为利用实时数据出炉实现对数据进行接收和融合结算的结果界面。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1所示，本发明实施例提供了一种基于UKF的实时多源数据融合方法，包括对各种测量数据(测速元素、测距元素和测角元素，以下简称测元)的预处理、测元选择策略和UKF滤波解算。具体步骤如下：
步骤一：数据的接收：通过数据接收线程接收测元，其中，测元包括测速元素、测距元素和测角元素；
步骤二：数据预处理：根据测元类型的不同，分别采用不同的预处理方法对不同类型的测元进行处理，并将预处理后的测元存放在数据缓存区中；
在测量设备实际工作中，难免会有测元数据异常，信号质量差，实测值和理论值不符的情况，以及由于传输时延等问题引起的数据时间不齐，所以，在解算目标运动参数前，需要进行一系列的预处理工作。预处理主要流程如图2所示。对于不同类型的测元，需要选择不同的处理方法，所以分别针对距离和变化率等测速元素以及距离测元和角度测元进行预处理。
S2.1、测速元素预处理方法
测速元素为连续波体制的测量元素距离和变化率
预处理包括量纲转换、异常值修正、时间误差修正和电波折射修正；
量纲转换：将设备输出的二进制数据恢复成十进制物理量，并将直接物理量多普勒频率信息转换成距离和变化率，以便后续处理；
异常值修正：由于设备故障或数据记录、判读过程异常，周围环境突变和干扰以及人为操作错误等原因，测量设备的观测数据量通常带有误差偏大的异常值，并影响后续融合解算。所以需要在预处理时将观测量中的异常值进行判别和处理，以保证测量数据的合理可信。因此，利用数据的连续性判断观测量中的数据是否异常，采用四点外推拟合法对观测量中的异常值进行处理；
时间误差修正：由于电磁波在空间传播的时间延迟，导致在离散时刻接收并记录的观测数据不能反映该采样时刻的目标运动状态，而是较早的某个时刻的状态。假设采样时刻tk测量设备接收到来自目标tk'时刻转发的电波信号，电波信号在空间传输时延为
式中，R(tk')为tk'时刻目标与观测设备的径向距离，c为光速；将tk'时刻的观测量向后修正Δtk，得到tk时刻对应的观测量；
电波折射修正：采用H.S.Hopfield折射指数模型将电磁波在经过不均匀大气层时产生传播路径弯曲折射进行修正，缩小实际测量值与理论观测值的差异。
S2.2、测距元素的预处理方法
测距元素为单脉冲体制的测量元素径向距离R，预处理包括异常值修正、时间误差修正和电波折射修正；
异常值修正：由于设备故障或数据记录、判读过程异常，周围环境突变和干扰以及人为操作错误等原因，测量设备的观测数据量通常带有误差偏大的异常值，并影响后续融合解算。所以需要在预处理时将观测量中的异常值进行判别和处理，以保证测量数据的合理可信。因此，利用数据的连续性判断观测量中的数据是否异常，采用四点外推拟合法对观测量中的异常值进行处理；
时间误差修正：由于电磁波在空间传播的时间延迟，导致在离散时刻接收并记录的观测数据不能反映该采样时刻的目标运动状态，而是较早的某个时刻的状态。假设采样时刻tk测量设备接收到来自目标tk'时刻转发的电波信号，电波信号在空间传输时延为
式中，R(tk')为tk'时刻目标与观测设备的径向距离，c为光速；将tk'时刻的观测量向后修正Δtk，得到tk时刻对应的观测量；
电波折射修正：采用H.S.Hopfield折射指数模型将电磁波在经过不均匀大气层时产生传播路径弯曲折射进行修正，缩小实际测量值与理论观测值的差异。
S2.3、测角元素的预处理方法
测角元素为单脉冲体制的方位角和俯仰角AE，需要对方位角和俯仰角进行分别处理；方位角预处理包括角坐标转换、方位角跨零跳处理、异常值修正和时间误差修正；俯仰角预处理包括角坐标转换、异常值修正、时间误差修正和电波折射修正；
角坐标转换：方位角和俯仰角的测量是在垂线测量坐标系下测得的，将测得的方位角和俯仰角转换成统一坐标系下，以参与融合解算；
方位角跨零跳处理：方位角跨零跳是指目标在I象限和IV象限之间移动时，直接观测量方位角会出现0°和360°之间的跳跃，导致方位角在数值上的不连续，并作为异常值识别出来，导致数据错判，影响后续处理，此时对方位角进行差分计算，得到差分序列，如果差分序列中有结果大于180°，则对应位置的方位值加360°，如果差分结果小于-180°，则对应位置的方位值减360°，从而进行方位角跨零跳处理，使得方位角在数值上连续有效；
异常值修正：由于设备故障或数据记录、判读过程异常，周围环境突变和干扰以及人为操作错误等原因，测量设备的观测数据量通常带有误差偏大的异常值，并影响后续融合解算。所以需要在预处理时将观测量中的异常值进行判别和处理，以保证测量数据的合理可信。因此，利用数据的连续性判断观测量中的数据是否异常，采用四点外推拟合法对观测量中的异常值进行处理；
时间误差修正：由于电磁波在空间传播的时间延迟，导致在离散时刻接收并记录的观测数据不能反映该采样时刻的目标运动状态，而是较早的某个时刻的状态。假设采样时刻tk测量设备接收到来自目标tk'时刻转发的电波信号，电波信号在空间传输时延为
式中，R(tk')为tk'时刻目标与观测设备的径向距离，c为光速；将tk'时刻的观测量向后修正Δtk，得到tk时刻对应的观测量；
电波折射修正：采用H.S.Hopfield折射指数模型将电磁波在经过不均匀大气层时产生传播路径弯曲折射进行修正，缩小实际测量值与理论观测值的差异。
步骤三：数据解算：利用测元选择策略从数据缓存区中查找当前时刻可用测元，生成测元矩阵，并通过UKF滤波算法对测元矩阵进行求解，得到并显示目标运动参数。
本发明以目前使用比较广泛的测量设备以及先验的测元精度设定，主要分析连续波测量雷达(距离和变化率
方位俯仰角(AE)L)，单脉冲雷达(径向距离、方位俯仰角(RAE)D)，相控阵雷达(距离和变化率、方位俯仰角(RAE)X)的测量数据相互融合的解算模式。
测元使用优先级取决于测元的精度高低，假定测元精度高低关系为
则对应的测元使用优先级设定也是用同样的关系表达。
当解算命令触发后，从数据缓存区中进行测元选择，如图3所示，测元选择具体步骤如下：
a)以当前解算时刻Tj为基准，在数据缓存区中查找该时刻的可用测元；
b)判断查找的可用测元是否为异常值，如果是，则可用测元不写入测元矩阵；如果不是，则可用测元写入测元矩阵；并重复步骤a)直至遍历数据缓存区；
c)判断测元矩阵中的非零测元个数，如果非零测元个数小于6个，则该测元矩阵直接参与解算；如果非零测元个数大于等于6个，则从最高优先级的测元开始由高到低对非零测元个数进行判断，如果当前优先级及之前的测元数量大于等于6个，则将低于该优先级的测元数据置零，否则继续遍历优先级较低的测元数据，直到有效测元个数大于等于6个或者遍历整个测元矩阵为止。
UKF滤波以UT变换为基础，通过采样策略逼近非线性函数分布，设计少量的σ点，计算这些σ点经由非线性函数的传播，从而获得滤波值基于非线性状态方程的更新。通过UKF滤波算法对测元矩阵进行求解的方法为：
I)U变换取样
第k个测量时刻的原始状态变量为：
原始状态方差为：
式中：
原始状态变量xk为n维列向量；其中，xk表示地心直角坐标系下X方向的位置，
表示地心直角坐标系下X方向的速度，
表示地心直角坐标系下X方向的加速度，zk表示地心直角坐标系下Z方向的位置，
表示地心直角坐标系下Z方向的速度，
表示地心直角坐标系下Z方向的加速度；
对原始状态变量xk进行U变换得到变量Xk，Xk为n×(2n+1)维矩阵，Xk具有如下表达式：
式中：Pk为xk的协方差矩阵；
系数λ＝α2(n+κ)-n，α和κ均为待选参数；
II)UKF的时间更新
变量Xk的状态随时间的更新表达式为式(1)，变量Xk对应的加权采样均值
的表达式为式(2)，变量Xk的协方差矩阵表达式为式(3)，变量Xk对应的采样观测量Zi,k的表达式为式(4)，采样观测量Zi,k对应的加权采样均值
的表达式为式(5)：
Xi,k＝F(Xi,k-1,wi,k-1) (1)
Zi,k＝H(Xi,k,vi,k) (4)
其中，Xi,k表示k时刻变量Xk的第i列，Xi,k-1表示k-1时刻变量Xk的第i列，wi,k-1表示k-1时刻变量状态量噪声矩阵wk-1的第i列，F(·)表示状态转移函数，Wi表示变量Xi,k的加权系数，Pk'表示变量Xk的协方差矩阵，Wi c表示变量Xi,k均值的加权系数，vi,k表示k时刻观测噪声矩阵的第i列，H(·)表示观测函数，i＝0,1,2,...,2n表示i的取值范围为0～2n的正整数，其中n表示原始状态变量xk的维数；
III)UKF的测量更新
采样观测量Zi,k的协方差矩阵表达式为式(6)，变量Xk与采样观测量Zi,k的互协方差表达式为式(7)，则卡尔曼增益Kk的表达式为式(8)，状态量及状态协方差更新表达式分别为式(9)和(10)：
其中，
为采样观测量Zi,k的协方差矩阵，
为变量Xk与采样观测量Zi,k的互协方差，
为
的逆矩阵，
为卡尔曼增益Kk的转置。
具体实例
模拟7个测量设备对180°射向的运动目标进行实时测量融合解算，7个测量设备以及目标地面轨迹投影的分布如图4所示。其中，1～4号测量设备为连续波测量设备，1号和2号为两个主站，3号、4号和7号为三个副站，5号为单脉冲测量设备，6号为相控阵测量设备。
利用MATLAB仿真工具，分别模拟七个测量设备的原始数据，其中，测速测元的随机差为0.02m/s，系统差0.0045m/s；测距测元的随机差为7m，系统差为5m；测角随机差和系统差均为0.0115°。模拟发送软件实现测量数据的实时发送，实时数据处理软件对数据进行接收和融合解算，并显示目标轨迹。
采用上述模拟条件分别模拟连续波测速测元、单脉冲和相控阵单独进行解算以及三种设备联合解算的结果，并分别将实时处理软件的解算结果作为比对目标与基准目标(理论轨迹)进行目标运动参数的比对，在25s～170s的时间段内，目标运动参数的位置和速度三个方向的比对残差如图5所示，统计比对残差的均方根误差如表1所示。
表1四种测元组成模式解算结果比对统计结果
连续波单独解算有效解算时间为60s～170s，单脉冲、相控阵以及三种融合解算的有效解算时间为30s～170s。并且从表1中可以看出多设备融合解算的位置和速度均方根误差均小于连续波、单脉冲和相控阵单独解算的结果。
连续波测速设备相较于其他设备从60s才输出有效目标轨迹参数，60s之前因为7号站位置较远，俯仰角低于保精度最低俯仰角限制导致数据无效，连续波设备没有足够的有效测元构成独立解算条件。当达到独立解算条件并且稳定后，在局部段落也能获得较高的的位置和速度精度。
单脉冲和相控阵设备则都是在25s开始输出有效运动参数，并且位置精度相对较高，但因为没有有效高精度的测速信息，所以解算的速度精度相比多测速解算结果较差；由于相控阵设备的测量精度低于单脉冲设备，所以相控阵单独解算的位置和速度结果也相比较差。
三种设备融合解算，能够从25s开始输出较好的解算轨迹。在解算初始阶段连续波设备不满足独立解算时，联合其他两种设备的测元进行融合解算，有效延长测量弧段，获取解算全程的高精度测量结果。
利用本发明提出的一种基于UKF的实时多源数据融合方法，通过Qt软件开发平台设计开发一款实时数据处理软件。利用实时数据处理实现对数据进行接收和融合解算，并显示目标轨迹。参数设置如图6所示，解算结果如图7所示。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (6)
1.一种基于UKF的实时多源数据融合方法，其特征在于，其步骤如下：
步骤一：数据的接收：通过数据接收线程接收测元，其中，测元包括测速元素、测距元素和测角元素；
步骤二：数据预处理：根据测元类型的不同，分别采用不同的预处理方法对不同类型的测元进行处理，并将预处理后的测元存放在数据缓存区中；
步骤三：数据解算：利用测元选择策略从数据缓存区中查找当前时刻可用测元，生成测元矩阵，并通过UKF滤波算法对测元矩阵进行求解，得到并显示目标运动参数。
2.根据权利要求1所述的基于UKF的实时多源数据融合方法，其特征在于，测速元素的预处理包括量纲转换、异常值修正、时间误差修正和电波折射修正；
量纲转换：将设备输出的二进制数据恢复成十进制物理量，并将直接物理量多普勒频率信息转换成距离和变化率；
异常值修正：利用数据的连续性判断观测量中的数据是否异常，采用四点外推拟合法对观测量中的异常值进行处理；
时间误差修正：采样时刻tk测量设备接收到来自目标tk'时刻转发的电波信号，电波信号在空间传输时延为
式中，R(tk')为tk'时刻目标与观测设备的径向距离，c为光速；将tk'时刻的观测量向后修正Δtk，得到tk时刻对应的观测量；
电波折射修正：采用H.S.Hopfield折射指数模型将电磁波在经过不均匀大气层时产生传播路径弯曲折射进行修正，缩小实际测量值与理论观测值的差异。
3.根据权利要求1或2所述的基于UKF的实时多源数据融合方法，其特征在于，测距元素的预处理包括异常值修正、时间误差修正和电波折射修正；
异常值修正：利用数据的连续性判断观测量中的数据是否异常，采用四点外推拟合法对观测量中的异常值进行处理；
时间误差修正：采样时刻tk测量设备接收到来自目标tk'时刻转发的电波信号，电波信号在空间传输时延为
式中，R(tk')为tk'时刻目标与观测设备的径向距离，c为光速；将tk'时刻的观测量向后修正Δtk，得到tk时刻对应的观测量；
电波折射修正：采用H.S.Hopfield折射指数模型将电磁波在经过不均匀大气层时产生传播路径弯曲折射进行修正，缩小实际测量值与理论观测值的差异。
4.根据权利要求3所述的基于UKF的实时多源数据融合方法，其特征在于，测角元素的预处理包括方位角和俯仰角的预处理；方位角预处理包括角坐标转换、方位角跨零跳处理、异常值修正和时间误差修正；俯仰角预处理包括角坐标转换、异常值修正、时间误差修正和电波折射修正；
角坐标转换：方位角和俯仰角的测量是在垂线测量坐标系下测得的，将测得的方位角和俯仰角转换成统一坐标系下；
方位角跨零跳处理：方位角跨零跳是指目标在I象限和IV象限之间移动时，直接观测量方位角会出现0°和360°之间的跳跃，导致方位角在数值上的不连续，此时对方位角进行差分计算，得到差分序列，如果差分序列中有结果大于180°，则对应位置的方位值加360°，如果差分结果小于-180°，则对应位置的方位值减360°，从而进行方位角跨零跳处理，使得方位角在数值上连续有效；
异常值修正：利用数据的连续性判断观测量中的数据是否异常，采用四点外推拟合法对观测量中的异常值进行处理；
时间误差修正：采样时刻tk测量设备接收到来自目标tk'时刻转发的电波信号，电波信号在空间传输时延为
式中，R(tk')为tk'时刻目标与观测设备的径向距离，c为光速；将tk'时刻的观测量向后修正Δtk，得到tk时刻对应的观测量；
电波折射修正：采用H.S.Hopfield折射指数模型将电磁波在经过不均匀大气层时产生传播路径弯曲折射进行修正，缩小实际测量值与理论观测值的差异。
5.根据权利要求1或4所述的基于UKF的实时多源数据融合方法，其特征在于，所述利用测元选择策略从数据缓存区中查找当前时刻可用测元，生成测元矩阵的方法为：
当解算命令触发后，从数据缓存区中进行测元选择，测元选择具体步骤如下：
a)以当前解算时刻Tj为基准，在数据缓存区中查找该时刻的可用测元；
b)判断查找的可用测元是否为异常值，如果是，则可用测元不写入测元矩阵；如果不是，则可用测元写入测元矩阵；并重复步骤a)直至遍历数据缓存区；
c)判断测元矩阵中的非零测元个数，如果非零测元个数小于6个，则该测元矩阵直接参与解算；如果非零测元个数大于等于6个，则从最高优先级的测元开始由高到低对非零测元个数进行判断，如果当前优先级及之前的测元数量大于等于6个，则将低于该优先级的测元数据置零，否则继续遍历优先级较低的测元数据，直到有效测元个数大于等于6个或者遍历整个测元矩阵为止。
6.根据权利要求5所述的基于UKF的实时多源数据融合方法，其特征在于，所述通过UKF滤波算法对测元矩阵进行求解的方法为：
I)U变换取样
第k个测量时刻的原始状态变量为：
原始状态方差为：
式中：
原始状态变量xk为n维列向量；其中，xk表示地心直角坐标系下X方向的位置，
表示地心直角坐标系下X方向的速度，
表示地心直角坐标系下X方向的加速度，zk表示地心直角坐标系下Z方向的位置，
表示地心直角坐标系下Z方向的速度，
表示地心直角坐标系下Z方向的加速度；
对原始状态变量xk进行U变换得到变量Xk，Xk为n×(2n+1)维矩阵，Xk具有如下表达式：
式中：Pk为xk的协方差矩阵；
系数λ＝α2(n+κ)-n，α和κ均为待选参数；
II)UKF的时间更新
变量Xk的状态随时间的更新表达式为式(1)，变量Xk对应的加权采样均值
的表达式为式(2)，变量Xk的协方差矩阵表达式为式(3)，变量Xk对应的采样观测量Zi,k的表达式为式(4)，采样观测量Zi,k对应的加权采样均值
的表达式为式(5)：
Xi,k＝F(Xi,k-1,wi,k-1) (1)
Zi,k＝H(Xi,k,vi,k) (4)
其中，Xi,k表示k时刻变量Xk的第i列，Xi,k-1表示k-1时刻变量Xk的第i列，wi,k-1表示k-1时刻变量状态量噪声矩阵wk-1的第i列，F(·)表示状态转移函数，Wi表示变量Xi,k的加权系数，Pk'表示变量Xk的协方差矩阵，Wi c表示变量Xi,k均值的加权系数，vi,k表示k时刻观测噪声矩阵的第i列，H(·)表示观测函数，i＝0,1,2,...,2n表示i的取值范围为0～2n的正整数，其中n表示原始状态变量xk的维数；
III)UKF的测量更新
采样观测量Zi,k的协方差矩阵表达式为式(6)，变量Xk与采样观测量Zi,k的互协方差表达式为式(7)，则卡尔曼增益Kk的表达式为式(8)，状态量及状态协方差更新表达式分别为式(9)和(10)：
其中，
为采样观测量Zi,k的协方差矩阵，
为变量Xk与采样观测量Zi,k的互协方差，
为
的逆矩阵，
为卡尔曼增益Kk的转置。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110912112.4A CN113743475B (zh) | 2021-08-10 | 2021-08-10 | 一种基于ukf的实时多源数据融合方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110912112.4A CN113743475B (zh) | 2021-08-10 | 2021-08-10 | 一种基于ukf的实时多源数据融合方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113743475A true CN113743475A (zh) | 2021-12-03 |
| CN113743475B CN113743475B (zh) | 2024-05-17 |
Family
ID=78730531
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110912112.4A Active CN113743475B (zh) | 2021-08-10 | 2021-08-10 | 一种基于ukf的实时多源数据融合方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113743475B (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116566450A (zh) * | 2023-07-10 | 2023-08-08 | 成都华兴大地科技有限公司 | 一种基于zynq的波束控制算法实现方法 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104252178A (zh) * | 2014-09-12 | 2014-12-31 | 西安电子科技大学 | 一种基于强机动的目标跟踪方法 |
| CN107300697A (zh) * | 2017-06-07 | 2017-10-27 | 南京航空航天大学 | 基于无人机的运动目标ukf滤波方法 |
| CN108759838A (zh) * | 2018-05-23 | 2018-11-06 | 安徽科技学院 | 基于秩卡尔曼滤波器的移动机器人多传感器信息融合方法 |
| CN110792430A (zh) * | 2019-11-20 | 2020-02-14 | 中国地质大学（北京） | 一种基于多传感器数据融合的随钻测斜方法及装置 |
| CN111708773A (zh) * | 2020-08-13 | 2020-09-25 | 江苏宝和数据股份有限公司 | 一种多源科创资源数据融合方法 |
| CN112068092A (zh) * | 2020-08-31 | 2020-12-11 | 西安工业大学 | 一种用于高精度弹道实时定轨的抗差加权观测融合平方根ukf滤波方法 |
-
2021
- 2021-08-10 CN CN202110912112.4A patent/CN113743475B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104252178A (zh) * | 2014-09-12 | 2014-12-31 | 西安电子科技大学 | 一种基于强机动的目标跟踪方法 |
| CN107300697A (zh) * | 2017-06-07 | 2017-10-27 | 南京航空航天大学 | 基于无人机的运动目标ukf滤波方法 |
| CN108759838A (zh) * | 2018-05-23 | 2018-11-06 | 安徽科技学院 | 基于秩卡尔曼滤波器的移动机器人多传感器信息融合方法 |
| CN110792430A (zh) * | 2019-11-20 | 2020-02-14 | 中国地质大学（北京） | 一种基于多传感器数据融合的随钻测斜方法及装置 |
| CN111708773A (zh) * | 2020-08-13 | 2020-09-25 | 江苏宝和数据股份有限公司 | 一种多源科创资源数据融合方法 |
| CN112068092A (zh) * | 2020-08-31 | 2020-12-11 | 西安工业大学 | 一种用于高精度弹道实时定轨的抗差加权观测融合平方根ukf滤波方法 |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116566450A (zh) * | 2023-07-10 | 2023-08-08 | 成都华兴大地科技有限公司 | 一种基于zynq的波束控制算法实现方法 |
| CN116566450B (zh) * | 2023-07-10 | 2023-09-12 | 成都华兴大地科技有限公司 | 一种基于zynq的波束控制算法实现方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113743475B (zh) | 2024-05-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN103383261B (zh) | 一种改进型无损卡尔曼滤波室内动目标定位方法 | |
| CN107315171B (zh) | 一种雷达组网目标状态与系统误差联合估计算法 | |
| CN110738275B (zh) | 基于ut-phd的多传感器序贯融合跟踪方法 | |
| CN111913484B (zh) | 一种变电站巡检机器人在未知环境下的路径规划方法 | |
| CN113743475B (zh) | 一种基于ukf的实时多源数据融合方法 | |
| CN109917373B (zh) | 运动补偿搜索的动平台雷达的动态规划检测前跟踪方法 | |
| CN110856104B (zh) | 一种联合最小二乘和三边定位的超宽带室内定位方法 | |
| WO2022241991A1 (zh) | 高超声速飞行器轨迹跟踪方法 | |
| CN111693051B (zh) | 一种基于光电传感器的多目标数据关联方法 | |
| CN113189541B (zh) | 一种定位方法、装置及设备 | |
| CN111679269B (zh) | 一种基于变分的多雷达融合航迹状态估计方法 | |
| CN116047495B (zh) | 一种用于三坐标雷达的状态变换融合滤波跟踪方法 | |
| CN114722455B (zh) | 一种联合全站仪和激光跟踪仪的三维工程控制网构建方法 | |
| CN114252871B (zh) | 一种基于机器学习的雷达测量精度补偿方法 | |
| CN116108314A (zh) | 一种基于多传感器组合测量空间目标轨道确定方法 | |
| CN116027320A (zh) | 一种基于多因素欧式距离关联的雷达与ais数据融合方法 | |
| CN115392117B (zh) | 一种水下高速机动平台高帧率无模糊声学导航方法 | |
| CN116956647B (zh) | 一种气动数据融合方法及系统 | |
| Han et al. | Deviation Registration of Radar Networking Based on Improved Cuckoo Search Algorithm | |
| CN116993832A (zh) | 一种基于cnn的飞行目标的位置预测方法及装置 | |
| CN112556722B (zh) | 一种基于自动选取优选源的系统误差补偿方法 | |
| CN118050722A (zh) | 一种考虑系统误差修正的双雷达多目标实时跟踪方法 | |
| He et al. | Dual-station direction finding location based on disturbance detection | |
| CN117630853A (zh) | 一种外辐射源雷达目标跟踪方法及计算机可读介质 | |
| Guo et al. | System Error Registration Method of Bearings-Only Passive Location System for Surface Ships |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |