logo资料库

AIS航迹校正.docx

第1页 / 共10页
第2页 / 共10页
第3页 / 共10页
第4页 / 共10页
第5页 / 共10页
第6页 / 共10页
第7页 / 共10页
第8页 / 共10页
资料共10页,剩余部分请下载后查看
1.任务背景
2.相关技术介绍
2.1AIS数据格式分析
2.1.1AIS设备配备规定
2.1.2AIS数据质量分析
2.2卡尔曼滤波算法基本方程
2.2.1卡尔曼滤波算法的系统状态方程
2.2.2卡尔曼滤波算法的动态系统测量方程
2.2.3无迹变换
2.2.4 UKF算法
3 基于卡尔曼滤波算法的船舶AIS轨迹建模
3.1 场景设置
3.2 参数设置
3.3 仿真结果评估
3.4 仿真结果和分析
4结论
1. 任务背景 近年来,随着我国经济的迅速发展,国内外贸易迅猛发展,海上交通运输量不断增长, 特 别 是 某 些 重 要 的 水 道 , 船 舶 交 通 非 常 繁 忙 ; 随 着 装 载 船 舶 自 动 识 别 系 统 (AutomaticIdentificationSystem,AIS)设备的船舶数量越来越多,充分利用船舶 AIS 数据进行船舶交通研究已是焦点问题。实时获取准确的船舶 AIS 信息显得至关重要,只有获 取实时的船舶 AIS 数据,才能有效地监控某些水域通航状况,及时发现通航水域中存在碰撞 以及搁浅等海事事故的区域,这样有助于岸基人员更好地协调和管理通航水域的安全。尽管 在船舶和岸基管理已经安装了大量的 AIS 设备,但是由于 AIS 船站和岸站时隙拥堵和网络传 输,造成 AIS 数据位置报告更新不及时等情况。为了弥补 AIS 数据堵塞等原因导致更新数据 不及时,造成船舶轨迹的不准确或者误差较大的问题,本文提出利用对卡尔曼滤波算法进行 适当的修改,引入系统噪声和测量噪声,利用 AIS 船舶观测节点数据对系统状态做最小二乘 法估计,对船舶轨迹进行平滑和预测处理,能够比较正确地估计出船舶轨迹。 2. 相关技术介绍 2.1AIS 数据格式分析 2.1.1AIS 设备配备规定 目前,按照国际海事组织的规定,国际航行的 500 总吨及以上的船舶以及所有客船均需 强制配备安装 AIS。我国海事主管机关对中国籍 100 总吨以上的沿海航行船舶以及内河航行 船舶都配备安装 AIS 的要求。 2.1.2AIS 数据质量分析 船舶位置报告的更新延时长短船舶轨迹的分析研究产生较大的影响,图 1 是选取了上海 港外高桥水域内自 2011 年 4 月 27 日 0805—4 月 28 日 0725 时间段内,每隔 35min 统计一次 该段时间内船舶位置报告的更新延时情况。 图 1 上海港历史 AIS 数据更新率对比图 从图 1 分析可得,2min 未更新位置报告的船舶比例为 30%左右,3min 未更新位置报 告的船舶比例为 15%左右,5min 未更新位置报告的船舶比例为 10%左右,15min 未更新位置 报告的船舶比例约为 2%。按照 AIS 设备规范的规定,左边是 A 类船载 AIS 设备,右边是 B 类船载 AIS 设备,船载 AIS 设备的报告间隔。通过与船载 AIS 设备的标准时间间隔对比,实 测 AIS 数据的位置报告间隔确有一定的延时,这种定位信息的延时会对后期的船舶位置的同 步产生一定的误差,同时 AIS 位置报告数据中船舶位置信息也客观存在着误差,其主要源于
船载定位设备的定位误差,基站数据接收和记录延迟误差。 通过上面数据分析可知,5min 以内有 90%船舶 AIS 数据都进行了更新,仅有 2%的船舶 定位信息更新延时会超过 15min,考虑到船舶航行速度相对较慢,在研究中期对船舶数据更 新延时在 5min、10min 的数据进行处理和分析,可以满足船舶轨迹预测和估计的精度要求。 2.2 卡尔曼滤波算法基本方程 2.2.1 卡尔曼滤波算法的系统状态方程 系统状态方程为:X(k)=Φ(k,k-1)X(k-1)+Γ(k,k-1)U(k-1)+W(k-1)(1) 式中:X(k)为 n 维状态向量,反映的是系统 k 时刻的状态;Φ(k,k-1)是系统的 n×n 维状态转移矩阵,描述系统从 k-1 时刻状态转移到 k 时刻状态的参数矩阵;Γ(k,k-1)是 n×m 维输入控制项矩阵;U(k-1)是 k-1 时刻的控制输入向量;W(k)为 n 维均值为零的系统过程白 噪声序列,用于描述从一个状态转移到另一个状态的噪声或者误差。 2.2.2 卡尔曼滤波算法的动态系统测量方程 动态系统测量方程: Z(k)=H(k)X(k)+V(k)(2) 式中:Z(k)为 m 维观测向量,表示动态系统在 k 时刻被测量到的状态;X(k)为 n 维状态 向量,反映的是系统的 k 时刻的状态;V(k)为 m 维观测噪声序列,用于描述动态系统从一个 状态转移到另一个状态的噪声或者误差;H(k)为 m×n 维观测矩阵,观测噪声 V(k)为零均值 白噪声序列。一般情况下,假定系统噪声向量 W(k)和动态系统测量噪声向量 V(k)均为均值 为零的高斯白噪声,他们的自相关矩阵分别为: E(W(n)W(k)={Q(n),n=k 0,≠0 E(V(n)V(k)={P(n),n=k 0,≠0 (3) 假设状态的初始值 X(0)与 W(k),V(k)≥0 的时刻互不相关,并且噪声向量 W(k)与 V(k) 互不相关,即有:E(W(k)(k)=0。利用动态系统测量数据向量 Z(1),Z(2),Z(3),⋯,Z(k)求 系统状态向量 X(i)分量进行最小二乘估计,根据 i 和 n 的不同关系,卡尔曼滤波可以分为 滤波、平滑和预测 3 种情况。 2.2.3 无迹变换 Simon Julier 和 Jeffrey K. Uhlmann 在 1996 年提出了一种非常新颖的求解非线性变换 后的分布的期望和方差的思想——无迹变换 UT(Unscented Transformation)。 无迹变换是基于这样的一个思想:用固定数量的参数去近似一个高斯分布比近似任意的 非线性函数更容易。其实现原理为:在原先状态分布中按某一规则取一些点,使这些点的均 值和协方差等于原状态分布的均值和协方差;将这些点代入非线性函数中,相应得到非线性 函数值点集,通过这些点求取变换后的均值和协方差。如图 1 所示。由于无迹变换没有对非 线性函数进行线性化、没有忽略其高阶项,因而得到的后验分布的均值和协方差的估计比用 线性截断法得到的要精确。 设正态 n 维随机向量 X N X P  ( , 图 2 无迹变换原理 ) ,m 维随机向量Y 为 X 的某一非线性函数,即
Y  ( g X ) , X 的统计特性 X 、 XP 通过非线性函数 ( )g  进行传播得到Y 的统计特性 Y , YP 。无迹变换求Y 、 YP 的基本步骤为:根据 X 、 XP 设计一系列的点 i,i=0,1,⋯ , L,称其为点;基于设定的点(i=0,1,⋯ ,L)计算其经过 ( )g  的传播 i (i=0,1, ⋯ ,L);然后基于 i(i=0,1,⋯ ,L)计算Y 、 YP 。通常点的数量取为 2n+1,即 L=2n。 无迹变换可具体描述如下: 计算 2n+1 个点及其权系数 0 X   i  X  ( ( n )   ) , P X i i ,   1, n    i n X  ( ( n )   ) , P X i i ,   1, n ) ( mw 0  / ( n  )  ( ) cw 0  / ( n  )  (1   2 )    m ) ( w i  ( ) c w i  0.5 / ( n  ),  i   1, ,2 n )    2(   n 其中 n  ,决定点的散布程度,通常取一小的正值(如 0.01), ) 通常取为 0;用来描述 X 的分布信息(正态情况下的最优值为 2); ( 示矩阵平方根第 i 列(提示采用 matlab 函数中[V, D]=eig(A)函数实现,其中 V,D 分别为 矩阵 A 的特征向量和特征值矩阵,A 矩阵的平方根阵列向量可用 ch=V*(D.^0.5)求出,ch 矩阵中的第 i 列即为所求矩阵平方根第 i 列- ( iw (i=0,1,⋯ ,2n) P P )X i )X i ); 表 n n ( ( )m ( ) 为求一阶统计特性时的权系数; ( )c iw (i=0,1,⋯ ,2n)为求二阶统计特性时的权系数。 注意:点采样策略则采用比例采样修正算法,该算法中需要确定,和共 3 个 参数。查阅有关文献,给出了上述参数的一般范围:确定 x 周围 Sigma 点的分布,通常设 4 [1 e ,1) 一个 区间内较小的正数,这里取 0.01,对于高斯分布 2 为最优,而是一个比 例参数,通常设置为 0 或3 n ,这里仿真实验中取 0。  点 通 过 非 线 性 函 数 ( )g  的 传 播 为  i 从而可得,  i g ( ), i   0,1 ,2 n
Y 2 n   i  0 ) m ( w  i i P Y  2 n  i  0 ( ) c w i (  i  Y )(  i  Y T ) P XY  2 n  i  0 ( ) c w i (  i  X )(  i  Y T ) 无迹变换初看似乎是一种 Monte-Carlo 方法,其实不然。它只是取少数确定的点,而 不是随机地从给定分布中进行抽样。另外一点,它不是通常意义上的加权方法,其权系数应 该说是“广义”上的权系数,不一定都为正,其和也不一定为 1,因而不能将其理解为抽样 统计。 2.2.4 UKF 算法 传统的非线性滤波的方法主要是扩展卡尔曼滤波算法( EKF) ,但是该算法存在着精度不 高、稳定性差、对目标机动反应迟缓等缺点. 近年来,文献提出了一种非线性滤波算法- Unscented 卡尔曼滤波(UnscentedKalman Filter,即 UKF). 它是根据 Unscented 变化(无味 变换)和卡尔曼滤波相结合得到的一种算法. 这种算法主要运用卡尔曼滤波的思想,但是在 求解目标后续时刻的预测值和量测值时,则需要应用采样点来计算. UKF 通过设计加权点δ, 来近似表示 n 维目标采样点,计算这些δ点经由非线性函数的传播,通过非线性状态方程获 得更新后的滤波值 ,从而实现了对目标的跟踪. UKF 有效地克服了扩展卡尔曼滤波的估计精 度低、稳定性差的缺陷. 卡尔曼最初提出的滤波理论只适用于线性系统,Bucy,Sunahara 等人提出并研究了扩展 卡尔曼滤波(Extended Kalman Filter,简称 EKF),将卡尔曼滤波理论进一步应用到非线性 领域。EKF 的基本思想是将非线性系统线性化,然后进行卡尔曼滤波,因此 EKF 是一种次优 滤波。其后,多种二阶广义卡尔曼滤波方法的提出及应用进一步提高了卡尔曼滤波对非线性 系统的估计性能。二阶滤波方法考虑了 Taylor 级数展开的二次项,因此减少了由于线性化 所引起的估计误差,但大大增加了运算量,因此在实际中反而没有一阶 EKF 应用广泛。 在状态方程或测量方程为非线性时,通常采用扩展卡尔曼滤波(EKF)。EKF 对非线性函 数的 Taylor 展开式进行一阶线性化截断,忽略其余高阶项,从而将非线性问题转化为线性, 可以将卡尔曼线性滤波算法应用于非线性系统中。这样一来,解决了非线性问题。EKF 虽然 应用于非线性状态估计系统中已经得到了学术界认可并为人广泛使用,然而该种方法也带来 了两个缺点,其一是当强非线性时 EKF 违背局部线性假设,Taylor 展开式中被忽略的高阶 项带来大的误差时,EKF 算法可能会使滤波发散;另外,由于 EKF 在线性化处理时需要用雅 克比(Jacobian)矩阵,其繁琐的计算过程导致该方法实现相对困难。所以,在满足线性系统、 高斯白噪声、所有随机变量服从高斯(Gaussian)分布这 3 个假设条件时,EKF 是最小方差准 则下的次优滤波器,其性能依赖于局部非线性度。 无迹卡尔曼滤波是一种新型的滤波估计算法。UKF 以 UT 变换为基础,摒弃了对非线性 函数进行线性化的传统做法,采用卡尔曼线性滤波框架,对于一步预测方程,使用无迹(UT) 变换来处理均值和协方差的非线性传递,就成为 UKF 算法。UKF 是对非线性函数的概率密度 分布进行近似,用一系列确定样本来逼近状态的后验概率密度,而不是对非线性函数进行近 似,不需要求导计算 Jacobian 矩阵。UKF 没有线性化忽略高阶项,因此非线性分布统计量
的计算精度较高。基于上述优点,UKF 被广泛应用于导航、目标跟踪、信号处理和神经网络 学习等多个领域。 在 UKF 算法中,由于具有噪声项,需要对状态进行扩维处理,针对目标运动方程和定义 观测方程的系统,令 x a T [ x T w T T ] v ,UKF 算法的具体步骤如下: 初始条件为 ˆ x ) 0 P 0 x ( 0 x (( E E   0  ˆ x 0 )( x 0  ˆ x 0 T ) ) 状态的初始条件扩维,即 ˆ x ;0;0] E ˆ x x   ( ) [ a 0 a 0 0 P a 0  E (( x a 0  ˆ x a 0 )( x a 0  ˆ x a T ) ) 0 1)Sigma 点采样       P 0 0 0 0 Q 0 0 0 R      i k kχ a 采用某种采样策略,得到 k 时刻状态估计的点集 , | , 1,  i , L ,其中 L 为所采用的 iχ 的 iχ 为粒子 a x 采样策略的采样 Sigma 点个数。需要注意的是,此时状态维数为   n q m , 前 n 维组成的列向量, w iχ 为粒子 a iχ 的 1n 维到 n q 维组成的列向量, v iχ 为粒子 a iχ 的 n q 维到     n q m 维组成的列向量。 1 ˆ k kX 及预测误差协方差 1|k 2)运用无迹变换求基于状态方程的一步预测 1| kP 。 χ 计算 Sigma 点 , | ,( i i k k 1,   , L ) 通过状态方程的传播, χ x , i k   1| k F χ ( , x i k k k , |  1) ˆ X k 1|  k 1   L i  0 W i ( m ) χ x , i k 1|  k P k 1|  k  1  L i  0 ( ) c W i ( χ x , i k 1|  k ˆ X  )( χ x , i k 1|  k k 1|  k  ˆ X T ) k 1|  k ˆ k kX 、 1|k 3)求 1| kP 通过观测方程的传播 z , i k 1|  k  H χ ( x , i k , k k 1|   1) ˆ Z k 1|  k 1   L i  0 W i ( m ) z , i k 1|  k
P , ZZ k 1   P , X Z k 1   W i ( ) c ( z , i k 1|  k  ˆ Z )( z k 1|  k , i k 1|  k  ˆ Z T ) k 1|  k 1  L i  0 W i ( ) c ( χ , i k 1|  k  ˆ X )( z k 1|  k , i k 1|  k  ˆ Z ) T k 1|  k 1  L i  0 4)滤波更新 K k 1  P , XZ k P 1 , ZZ k  1  1  ˆ X k 1|  k 1   ˆ X  k 1|  k K Z ( 1  k  ˆ Z ) k 1|  k k 1  P k 1|  k 1   P k 1|  k  K P 1  k , ZZ k K 1  T k 1  采用无迹卡尔曼滤波(UKF)时,由于无需计算雅可比矩阵求导,因此,可以将这种滤波 方法直接应用于大地坐标系下的航迹滤波,同时可避免传统的“平台为中心”技术体制下复 杂繁琐的坐标转换过程。 3 基于卡尔曼滤波算法的船舶 AIS 轨迹建模 3.1 场景设置 根据对舰艇运动动力学的分析研究,可以认为舰艇机动是由切向加速度和法向加速度两 部分共同作用形成的,其中切向加速度改变舰艇速度的大小,法向加速度改变速度方向,在 推算下一时刻舰艇位置时,须把两个方向的加速度都考虑在内。舰艇在航行过程中的机动不 外乎是匀速、加速、减速、转弯、S 型机动。水面舰艇机动实际可以看作是上述几种机动模 型的任意组合,在仿真场景中加入各种机动模式的控制机制,在不同机动模式之间任意切换, 从而产生不同舰艇的各种机动情况,但要符合舰艇实际机动的一些限制,例如一般舰艇机动 的速度不超过 60 节,舰艇的加速度不超过 2g 等等。 给出的典型仿真场景中目标运动曲线和速度变化曲线分别如图 2、图 3 所示,其中采样 周期为 10s,初始角度为 45 度,运动过程如下: 第 01 周期到第 35 周期,目标匀速直线运动,速度为 10 米/秒; 第 35 周期到第 45 周期,目标做匀加速直线运动,加速度为 0.1 米/秒; 第 45 周期到第 60 周期,目标做匀减速直线运动,加速度为-0.1 米/秒; 第 60 周期到第 90 周期,目标做匀速右转运动,转动角度为 90 度,此时速度为 5 米/秒 第 90 周期到第 120 周期,目标接着做匀速左转运动,转动角度为 90 度; 第 120 周期到第 130 周期,目标做匀加速直线运动,加速度为 0.1 米/秒; 第 130 周期到第 180 周期,目标做匀速直线运动,此时速度为 15 米/秒; 目标初始位置经纬度坐标为( 100.0 ,40.0 ),该仿真场景包括匀速直线运动、匀加速、 匀减速直线运动,还有由匀速右转和左转所构成的蛇形机动(或称之为 S 型机动),覆盖了 CV、CA 和 CT 几种常用的运动模型,可以较好的检验滤波算法的性能。 0 0
40.1 40.09 40.08 40.07 40.06 40.05 40.04 40.03 40.02 40.01 ) 度 ( 度 纬 40 100 100.02 100.04 典型仿真场景 20 18 16 14 12 10 8 6 / ) 秒 米 ( 度 速 100.06 100.08 100.1 经度(度) 100.12 100.14 100.16 4 0 200 400 600 1000 800 时间(秒) 1200 1400 1600 1800 图 2 典型仿真场景运动曲线图 图 3 典型仿真场景速度曲线图 3.2 参数设置 仿真中,采用大地坐标系下二维 CT 模型来描述水面舰艇目标机动方式,状态向量为 X ]T [ v  , , , , ,初始状态滤波初值取为 X 0 [100 ,40 ,10 / m s 0 0 0 ,0 ,0 ]s 0/ ,动力学噪 声方差矩阵为 Q k diag 3.3 仿真结果评估 6  [10 ,10 ,100,1,1]  6 ,观测误差方差矩阵为 R k diag 2 2 [0.1 ,0.1 ] 。 滤波算法对目标的跟踪性能主要采用均方根误差(Root-mean-square error,RMSE)参 数来衡量,探测距离误差( S )、探测方位()、速度( v )RMSE 定义如下: ( ) RMSE S  1 N N  1 i  ˆ S i ( 2 S ) RMSE ( )   ; 1 N N  ˆ ) (   i i 1  2 ; ( ) RMSE v  1 N N  1 i  ˆ v i ( 2 v ) 其中 S 、、 v 为真实值, 结果通过图形的方式表现出来。 3.4 仿真结果和分析 ˆ iS 、 ˆi、 ˆiv 为第i 次滤波值。 N 为仿真次数。最后将计算
图 4 目标真实轨迹 图 5 目标观测航迹
分享到:
收藏