一种提高船舶自动识别系统测量精度的估计方法与流程

本发明属于海上船舶状态测量技术领域，具体涉及一种提高船舶自动识别系统测量精度的估计方法。

背景技术：

随着社会经济的快速发展，不断增长的船舶数量和交通密度严重影响着水域交通安全及海洋环境污染。频繁发生的海损事故，将会造成生命与财产的重大损失，而且对邻近水域的海洋生态环境会造成严重的破坏。为了确保海上航行安全，船与船、船与岸之间的实时数据交换以及目标信息识别等功能变得越来越重要。船舶自动识别系统(automaticidentificationsystem，ais)正是针对上述问题而发展起来的一种用于岸一船、船一岸以及船一船之间通信的新型助航系统。ais设备能将本船的动态信息如位置、速度、航向等和静态信息如海上移动业务标识码、船名、呼号等，经信息处理器处理后在vhf信道向外播发，本船周围船舶只要安装了ais设备就可以自动接收这些船舶信息并显示；同时，本船也可将其他船的识别信息和本船信息一起存储并向范围内其他目标船发送。然而，由于船舶模型不确定性、随机环境干扰与高频测量噪声等因素的存在，传统的ais系统得出的目标船舶运动状态信息存在较大的测量偏差，严重影响了ais系统对目标船舶的判别，这种测量偏差的存在，对船舶海上航行的安全性产生了巨大的威胁。

技术实现要素：

本发明的目的是为了解决水面船舶自动识别系统提供的船舶运动状态数据存在较大测量误差的问题，提供一种提高船舶自动识别系统测量精度的估计方法，该种方法基于水面船舶三自由度标称运动学模型与动力学模型，分析水面船舶三自由度模型中的模型不确定性与测量噪声，建立包含模型不确定性与测量噪声的船舶三自由度受扰模型，进而建立水面船舶三自由度测量方程，基于测量方程设计变结构扩展容积卡尔曼滤波器对船舶运动状态进行估计，利用ais系统对目标船舶进行跟踪，将估计出的船舶运动状态作为最终目标跟踪状态信息通过通信单元进行发送。

为实现上述目的，本发明采取的技术方案如下：

一种提高船舶自动识别系统测量精度的估计方法，所述的方法具体步骤如下：

步骤一：建立水面船舶三自由度标称运动学模型与动力学模型为：

其中，η＝[x,y,ψ]^t，ν＝[u,v,r]^t，x为北向位置，y为东向位置，ψ为船舶艏向角，u为纵荡速度，v为横荡速度，r为艏摇角速度，j(η)为坐标转换矩阵，τc为控制力和力矩的矩阵；

m0为船舶质量与惯性矩阵，其详细定义为：

其中，m为船体质量，iz为z轴转动惯量，xg为船体重心在船体坐标系下x轴的坐标，均是指相应自由度下的附加质量和附加静矩；

c0(ν)为科里奥利与向心力矩阵，其详细定义为：

是指相应自由度的附加质量和附加静矩；

d0(ν)为水动力阻尼矩阵，其详细定义为：

其中，xu、yv、yr、nv、nr分别为相应自由度的线性水动力阻尼；

步骤二：分析水面船舶三自由度标称动力学模型中的不确定性与测量噪声，建立包含模型不确定性与测量噪声的船舶三自由度受扰模型；

基于式(10)给出模型不确定性分析如下，设置以下变量：

m＝m0+mδ(11)

c(ν)＝c0(ν)+cδ(ν)(12)

d(ν)＝d0(ν)+dδ(ν)(13)

其中，带有下标0的表示矩阵中标称模型部分，带有下标δ的表示模型参数摄动部分，则含有模型不确定性的船舶三自由度受扰模型表示成以下形式：

其中，式(10)为带下标0的标称模型，式(14)为添加了模型不确定性的模型，含义不同；

由一阶波浪诱导产生的水面船舶三自由度高频运动模型由下式所示

其中，ξ为等效状态变量，aw为系统矩阵，ew为随机噪声矩阵，w1为零均值高斯白噪声向量，ηw＝cwξ代表一阶线性波浪诱导作用下的状态高频输出，cw为输出矩阵；

将公式(15)所描述的模型看作是一种附带非线性阻尼的二阶振荡环节，如下式所示：

其中，kωi为波浪强度度量值，s是为传递函数中拉氏变换的算子，阻尼比ζi的取值范围为0.05～0.2，ω0i为波浪p-m谱中的主导频率；

利用一阶高斯—马尔科夫过程模拟由海上环境因素引起的缓慢变化的环境扰动作用力，该慢变环境扰动模型的描述形式为：

其中，tb为二维对角阵，eb为随机噪声矩阵，w2为零均值高斯白噪声向量；

考虑到水面船舶导航与测量系统能够提供的是带有测量噪声的系统输出值，因此，测量方程由下式得出：

y＝η+ηw+v(18)

其中，v为零均值高斯白噪声向量；

综上，得到包含环境干扰与随机噪声特性的完整系统模型为：

为了便于步骤中变结构扩展容积卡尔曼滤波器的设计，提取式(19)中的第二与第四子方程离散化写成如下的一体化标准状态空间描述形式：

其中，xk＝[ηk,υk]^t，ηk和υk即表征离散后的第k时刻状态η和状态υ，下标k表征当前的k时刻，下标k-1表示k-1时刻，

uk-1＝τck-1，wk-1由于传感器高频噪声和随机海浪干扰产生的随机系统噪声，zk为k时刻系统的观测输出，h(·)＝η，vk为k时刻传感器测量噪声；

步骤三：基于测量方程设计变结构扩展容积卡尔曼滤波器对船舶运动状态进行估计；首先容积卡尔曼滤波ckf选取2n个cubature点及同等权值计算高斯权重积分，cubature点形式如下：

使用cubature点集得出ckf滤波算法，而后进行状态扩维，在每个滤波周期内进行时间更新和测量更新，n为自然数、ei为状态矩阵的第i个列元素；

对于基础型扩展ckf滤波器设计步骤包含以下几部分：

1.时间更新：

a.假设k-1时刻后验密度函数已知，通过cholesky分解误差协方差矩阵pk-1|k-1：

其中，n(·,·)高斯分布，为k-1时刻系统状态[η,υ]^t的估计值，sk-1|k-1为通过cholesky分解获取得到的正定对称阵，为sk-1|k-1的转置；

b.针对系统状态进行状态扩维(j＝1,2,...n,n＝2(n+n+m))：

n为状态的维数，m为噪声的维数，j是要进行取cubature点的数量

其中，为扩维后k-1时刻系统状态[η,υ]^t的估计值，q和r都为正定对称矩阵，即表示扩展之后的sk-1|k-1矩；

c.计算cubature点(j＝1,2,...n,n＝2(n+n+m))：

其中，xj,k-1|k-1为cubature点；

通过状态方程传播的cubature点

其中，r表示采样基准，x表示采样域，表示与系统状态维数相对应的cubature点，表示与状态噪声相对应的cubature点，uk-1为k-1时刻的控制输入；

d.估计k时刻的状态预测值如下：

其中，是在k时刻的预测值；

估计k时刻的状态误差协方差预测值pk|k-1如下：

2.测量更新：

a.通过观测方程传播的cubature点zi,k|k-1如下：

其中，v表示基准采样点，x表示采样域，为与系统状态维数相对应的cubature点，表示与系统观测噪声相对应的cubature点；

估计k时刻的观测预测值

zj,k|k-1为求平均数用的原始元素；

b.估计自相关协方差阵pzz,k|k-1：

估计互相关协方差阵pxz,k|k-1：

xj,k|k-1为cubature点的异步预报估计；ωj＝1/2n,j＝1,2,...,n；

估计卡尔曼增益wk：

c.k时刻的状态估计值

k时刻状态误差协方差估计值pk|k：

pk|k＝pk|k-1-wkpzz,k|k-1wk^t(35)

对于滑模变结构滤波器，

1.k时刻状态预测值：

其中，

2.估计k时刻的状态误差协方差预测值：

其中，即对函数f求取偏导数；qk-1代表k-1时刻的矩阵q，q为正定对称阵；uk就是k时刻的控制输入u；

3.则k时刻的观测预测值为：

其中，h(·)＝η；

设观测值z的估计误差变量ez,k|k-1为：

4.滤波器增益kk表示为：

其中，○为schur乘法算子，算子[·]⁺代表矩阵的广义逆，γ为标量，代表滤波器的收敛率增益参数，ez,k-1|k-1为k-1时刻系统观测误差变量，zk-1为k-1时刻的系统观测值，为k-1时刻系统输出的估计值；

滤波器增益的饱和函数定义为：

其中，为由滑动平面边界向量ψ导出的对角阵，带下标i的表征矩阵里面的某一个元素；不带i的表示整个矩阵；

ψm代表矩阵里面第m行第m列的那个元素中的分母，是滑动平面边界向量中的元素；

k时刻的状态估计值为：

其中，kk为滤波器增益参数，ez,k|k-1为估计误差；

k时刻的状态误差协方差估计值：

pk|k＝(i-kkhk-1)pk|k-1(i-kkhk-1)^t+kkrkkk^t(49)

其中，矩阵rk为待选择的正定对称阵；i为单位阵；pk|k-1代表正定对称阵p在k-1时刻的异步预报估计；

k时刻的观测估计值和观测误差估计值依次为：

定义滑动平面边界层矩阵ψ，该矩阵包含了不同滑动平面边界层之间的关系；

为了解决基于ψ的滑模面边界层随时间变化的问题，将式(49)和k时刻的状态误差协方差估计值方程联立可得：

其中，trace代表矩阵的迹；上式解定义为如下形式：

其中，sk和a分别定义为：

a＝(|ez,k|k-1|+γ|ez,k-1|k-1|)(55)

由可知，代表由a组成的对角阵，利用此对角阵计算得出的滑动平面边界方程，将扩展cubature卡尔曼滤波器与滑模变结构滤波器联合进行目标船舶三自由度状态估计；

步骤四：基于步骤一至三中对ais系统中目标水面船舶三自由度运动状态的最优估计，对目标船舶运动状态进行持续目标跟踪，将基于容积卡尔曼滤波器估计出的船舶运动状态作为最终目标跟踪状态信息通过通信单元进行发送，从而基于此判断本船是否通过危险区域。

本发明相对于现有技术的有益效果是：本发明采用一种基于扩展ckf的变结构滤波器作为船舶自动识别系统对目标船舶运动状态进行最优估计。该滤波器结合了ckf滤波器与滑模变结构滤波器的优点，既能提供含有高斯条件下测量噪声的二维系统状态信息的精确逼近，又对随机海浪干扰条件下水面船舶三自由度模型不确定性具有良好的鲁棒特性；由于利用状态估计策略对ais系统输出数据进行最优状态估计，因此提高了ais系统对目标船舶运动状态的检测精度。

附图说明

图1为基于变结构扩展ckf滤波的船舶自动识别系统工作流程图；

图2为变结构滤波器原理示意图。

具体实施方式

下面结合附图和实施例对本发明的技术方案作进一步的说明，但并不局限于此，凡是对本发明技术方案进行修正或等同替换，而不脱离本发明技术方案的精神范围，均应涵盖在本发明的保护范围之中。本发明中符号上带一个点的表示一阶导数，带两个点表示二阶导数。

具体实施方式一：本实施方式记载的是一种提高船舶自动识别系统测量精度的估计方法，所述的方法具体步骤如下：

步骤一：建立水面船舶三自由度标称运动学模型与动力学模型为：

其中，η＝[x,y,ψ]^t，ν＝[u,v,r]^t，x为北向位置，y为东向位置，ψ为船舶艏向角，u为纵荡速度，v为横荡速度，r为艏摇角速度，j(η)为坐标转换矩阵，τc为控制力和力矩的矩阵；

m0为船舶质量与惯性矩阵，其详细定义为：

其中，m为船体质量，iz为z轴转动惯量，xg为船体重心在船体坐标系下x轴的坐标，均是指相应自由度下的附加质量和附加静矩；

c0(ν)为科里奥利与向心力矩阵，其详细定义为：

是指相应自由度的附加质量和附加静矩；

d0(ν)为水动力阻尼矩阵，其详细定义为：

其中，xu、yv、yr、nv、nr分别为相应自由度的线性水动力阻尼；

步骤二：分析水面船舶三自由度标称动力学模型中的不确定性与测量噪声，建立包含模型不确定性与测量噪声的船舶三自由度受扰模型；

由于实际模型中存在诸多未建模动态和模型不确定性，因此需要基于式(10)给出模型不确定性分析如下，设置以下变量：

m＝m0+mδ(11)

c(ν)＝c0(ν)+cδ(ν)(12)

d(ν)＝d0(ν)+dδ(ν)(13)

其中，带有下标0的表示矩阵中标称模型部分，带有下标δ的表示模型参数摄动部分，则含有模型不确定性的船舶三自由度受扰模型表示成以下形式：

其中，式(10)为带下标0的标称模型，式(14)为添加了模型不确定性的模型，含义不同；

由一阶波浪诱导产生的水面船舶三自由度高频运动模型由下式所示

其中，ξ为等效状态变量，aw为系统矩阵，ew为随机噪声矩阵，w1为零均值高斯白噪声向量，ηw＝cwξ代表一阶线性波浪诱导作用下的状态高频输出，cw为输出矩阵，其数值的大小代表不同的海情状态；

将公式(15)所描述的模型看作是一种附带非线性阻尼的二阶振荡环节，如下式所示：

其中，kωi为波浪强度度量值，kωi越大，则波浪强度越大，s是为传递函数中拉氏变换的算子，阻尼比ζi的取值范围为0.05～0.2，ω0i为波浪p-m谱中的主导频率，与有义波高值有关；

利用一阶高斯—马尔科夫过程模拟由风浪流等海上环境因素引起的缓慢变化的环境扰动作用力，诸如二阶波浪漂移力、海流作用力以及海风等作用力。该慢变环境扰动模型的描述形式为：

其中，tb为二维对角阵，表征环境干扰力矩的时间常数，eb为随机噪声矩阵，w2为零均值高斯白噪声向量；

考虑到水面船舶导航与测量系统能够提供的是带有测量噪声的系统输出值，因此，测量方程由下式得出：

y＝η+ηw+v(18)

其中，v为零均值高斯白噪声向量，表征传感器系统测量噪声；

综上，得到包含环境干扰与随机噪声特性的完整系统模型为：

为了便于步骤中变结构扩展容积卡尔曼滤波器的设计，可以提取式(19)中的第二与第四子方程离散化写成如下的一体化标准状态空间描述形式：

其中，xk＝[ηk,υk]^t，ηk和υk即表征离散后的第k时刻状态η和状态υ，下标k表征当前的k时刻，下标k-1表示k-1时刻，uk-1＝τck-1，wk-1由于传感器高频噪声和随机海浪干扰产生的随机系统噪声，zk为k时刻系统的观测输出，h(·)＝η，vk为k时刻传感器测量噪声；

步骤三：基于测量方程设计变结构扩展容积卡尔曼滤波器对船舶运动状态进行估计；首先容积卡尔曼滤波ckf选取2n个cubature点及同等权值计算高斯权重积分，cubature点形式如下：

使用cubature点集得出ckf滤波算法，而后进行状态扩维，在每个滤波周期内进行时间更新和测量更新，n为自然数、ei为状态矩阵的第i个列元素；

对于基础型扩展ckf滤波器设计步骤包含以下几部分：

1.时间更新：

a.假设k-1时刻后验密度函数已知，通过cholesky分解误差协方差矩阵pk-1|k-1：

其中，n(·,·)高斯分布，为k-1时刻系统状态[η,υ]^t的估计值，sk-1|k-1为通过cholesky分解获取得到的正定对称阵，为sk-1|k-1的转置；

b.针对系统状态进行状态扩维(j＝1,2,...n,n＝2(n+n+m))：

n为状态的维数，m为噪声的维数，j是要进行取cubature点的数量

其中，为扩维后k-1时刻系统状态[η,υ]^t的估计值，q和r都为正定对称矩阵，即表示扩展之后的sk-1|k-1矩；

c.计算cubature点(j＝1,2,...n,n＝2(n+n+m))：

其中，xj,k-1|k-1为cubature点；

通过状态方程传播的cubature点

其中，r表示采样基准，x表示采样域，表示与系统状态维数相对应的cubature点，表示与状态噪声相对应的cubature点，uk-1为k-1时刻的控制输入；

d.估计k时刻的状态预测值如下：

其中，是在k时刻的预测值；

估计k时刻的状态误差协方差预测值pk|k-1如下：

2.测量更新：

a.通过观测方程传播的cubature点zi,k|k-1如下：

其中，v表示基准采样点，x表示采样域，为与系统状态维数相对应的cubature点，表示与系统观测噪声相对应的cubature点；

估计k时刻的观测预测值

zj,k|k-1为求平均数用的原始元素；

b.估计自相关协方差阵pzz,k|k-1：

估计互相关协方差阵pxz,k|k-1：

xj,k|k-1为cubature点的异步预报估计；ωj＝1/2n,j＝1,2,...,n；

估计卡尔曼增益wk：

c.k时刻的状态估计值

k时刻状态误差协方差估计值pk|k：

pk|k＝pk|k-1-wkpzz,k|k-1wk^t(35)

水面船舶在整个工作域范围内是一个典型的具有模型参数不确定性的强非线性动力系统，而滑模变结构滤波器对于该类系统表现出很好的稳定性和鲁棒性。滑模变结构滤波器的基本概念如图2所示。

对于滑模变结构滤波器，

1.k时刻状态预测值：

其中，

2.估计k时刻的状态误差协方差预测值：

其中，即对函数f求取偏导数；qk-1代表k-1时刻的矩阵q，q为正定对称阵；uk就是k时刻的控制输入u；

3.则k时刻的观测预测值为：

其中，h(·)＝η；

设观测值z的估计误差变量ez,k|k-1为：

4.滤波器增益kk表示为：

其中，○为schur乘法算子，算子[·]⁺代表矩阵的广义逆，γ为标量，代表滤波器的收敛率增益参数，ez,k-1|k-1为k-1时刻系统观测误差变量，zk-1为k-1时刻的系统观测值，为k-1时刻系统输出的估计值；

滤波器增益的饱和函数定义为：

其中，为由滑动平面边界向量ψ导出的对角阵，带下标i的表征矩阵里面的某一个元素；不带i的表示整个矩阵；

ψm代表矩阵里面第m行第m列的那个元素中的分母，是滑动平面边界向量中的元素；

k时刻的状态估计值为：

其中，kk为滤波器增益参数，ez,k|k-1为估计误差；

k时刻的状态误差协方差估计值：

pk|k＝(i-kkhk-1)pk|k-1(i-kkhk-1)^t+kkrkkk^t(49)

其中，矩阵rk为待选择的正定对称阵；i为单位阵；pk|k-1代表正定对称阵p在k-1时刻的异步预报估计；

k时刻的观测估计值和观测误差估计值依次为：

为了得到一个能够产生更精确估计值的滑动平面的边界方程，定义了如下的滑动平面边界层矩阵ψ，该矩阵包含了不同滑动平面边界层之间的关系；

为了解决基于ψ的滑模面边界层随时间变化的问题，将式(49)和k时刻的状态误差协方差估计值方程联立可得：

其中，trace代表矩阵的迹；上式解定义为如下形式：

其中，sk和a分别定义为：

a＝(|ez,k|k-1|+γ|ez,k-1|k-1|)(55)

式(54)、(55)中相关参数定义可参见前文定义；

由可知，代表由a组成的对角阵，利用此对角阵计算得出的滑动平面边界方程，可以将扩展cubature卡尔曼滤波器与滑模变结构滤波器联合进行目标船舶三自由度状态估计；

步骤四：基于步骤一至三中对ais系统中目标水面船舶三自由度运动状态的最优估计，对目标船舶运动状态进行持续目标跟踪，将基于容积卡尔曼滤波器估计出的船舶运动状态作为最终目标跟踪状态信息通过通信单元进行发送，从而基于此判断本船是否通过危险区域。