[发明专利]一种可估计未知有效声速的自适应水下单信标定位方法

权利要求书

1.一种可估计未知有效声速的自适应水下单信标定位方法，其特征在于，包括以下步骤：

A.以定位区域内任意点为原点，东、北、天三个方向分别设为x，y，z轴，建立水下局部惯性坐标系；

B.通过水下航行器所搭载的GPS系统获取该所述水下航行器在水下局部惯性系当中的初始位置；

C.建立所述水下航行器的运动学模型以及观测模型并进行离散化；

所述运动学模型的建立方法为：

定义状态向量为：

x＝[x y v_cx v_cy]^T

其中：x，y为所述水下航行器在所述水下局部惯性坐标系中的水平位置；v_cx，v_cy为未知的海流速度；

对x求导并加入所述水下航行器运动模型噪声影响，得到所述水下航行器的运动学模型：

其中：v_wx为所述水下航行器x方向的对水速度；v_wy为所述水下航行器y方向的对水速度；v_wx及v_wy通过读取螺旋桨转速与电子罗盘测得的航行器艏向角计算得出；ω_x为所述水下航行器在x方向的位置不确定性；ω_y为所述水下航行器在y方向的位置不确定性；ω_cx为x方向的海流不确定性；ω_cy为y方向的海流不确定性；

v_wx及v_wy的计算公式为：式中v_w为根据所述螺旋桨转速得出的所述水下航行器对水速度，为所述电子罗盘测得的艏向角；

所述观测模型的建立方法为：

S1.建立水声信号传递时间的观测模型；

设所述水下航行器获得水声信标发射水声信号的时刻为T_e，所述水声信标在所述水下局部惯性坐标系中的空间位置坐标为X_Te，Y_Te，Z_Te，所述水下航行器接收到该水声信号的时刻为T_a，观测方程为：

其中：ν_t为对应的观测噪声；z为所述水下航行器的深度，由深度计精确测得，为已知量；v_e为有效声速；

将观测方程记作m＝h(x,v_e)，其中：

S2.建立海流流速观测模型；

根据多普勒测速仪测得的所述水下航行器的绝对速度v_g，结合所述电子罗盘测得的艏向角计算得到所述水下航行器绝对速度在局部惯性坐标系下的分量v_gx,v_gy；

根据v_gx,v_gy以及v_wx,v_wy，计算得到海流速度观测分量为：

海流观测量为线性，满足m_vc＝Hx+ν_vc；

其中：观测向量m_vc＝[m_cx m_cy]^T；m_cx，m_cy分别为x，y方向海流速度观测；ν_vc为海流观测噪声向量，ν_vc＝[ν_v,cx ν_v,cy]^T，其中ν_v,cx为x方向的海流不确定性；ν_v,cy为y方向的海流不确定性；H为海流观测矩阵，满足：

所述运动学模型以及所述观测模型离散化方法为：

S1.运动学模型离散化；

以下标k为时间索引，以Δt＝t_k+1-t_k为离散间隔，运动学模型离散为：

x_k+1＝A_kx_k+B_ku_k+w_k

其中：A_k为运动学方程，满足：

B_k为控制方程，满足：

u_k为控制向量，满足u_k＝[v_wx,k v_wy,k]^T，为已知量；

w_k过程噪声向量，满足w_k＝[ω_x,k ω_y,k ω_cx,k ω_cy,k]^T，对应各个状态变量的不确定性；将系统状态x_k,y_k,v_cx,k,v_cy,k的过程噪声建模为零均值Gauss分布，其协方差矩阵满足：

其中，σ_w为所述水下航行器对水速度观测不确定性的标准差；σ_c为海流不确定性的标准差；

S2.观测模型离散化；

所述水下航行器在k-1至k间接收到该水声信号，将其假设为在k时刻接收到该水声信号，即离散后的水声信号传递时间观测方程为：

其中，ν_t,k为观测噪声，假设其满足方差为R_t,k的Gauss分布；考虑到有效声速v_e,k的时变未知性，将v_e,k也看作随机变量；离散形式的观测方程写作：

m_k＝h_k(x_k,v_e,k)，

由于海流观测采样频率较高，假设在每一个离散时间点k处均可以得到海流观测，故离散后的海流观测方程为：

m_vc，k＝H_kx_k+ν_vc,k

其中，H_k为k时刻海流速度观测矩阵，满足：

ν_vc,k为k时刻海流观测噪声，为零均值Gauss分布，其观测噪声协方差矩阵满足：

其中，σ_vc,m为海流速度观测噪声的标准差；

D.建立有效声速及其不确定性统计参数的随机模型以及函数模型；

将有效声速初始先验分布建模为Gauss分布：

其中：N(x|μ,Σ)表示以μ为均值，Σ为协方差矩阵，满足Gauss分布的随机向量x；为初始声速估计均值；P_e,0为初始声速不确定性方差；

将k-1时刻有效声速的后验分布也建模为Gauss分布：

P_e,k-1|k-1分别为k-1时刻有效声速的后验均值与方差；

有效声速的动态扩散的函数模型记作：

v_e,k＝v_e,k-1+ω_e,k-1

其中：ω_e,k-1为有效声速的不确定性，将其建模为均值及方差均未知的Gauss分布：

其中：μ_k-1与σ_e,k-1分别为声速不确定性的均值与标准差；

据此，可以将有效声速的概率扩散方程记作：

结合上式，获得有效声速的预测模型为：

其中P_e,k|k-1分别为k时刻有效声速的先验均值与方差，其计算公式为：

由于真实的有效声速的不确定性均值和方差未知，以其名义值与进行声速预测，预测方程为：

其中：与表示有效声速先验统计参数，即均值与方差的名义值；

通过在线估计P_e,k|k-1来补偿σ_e,k-1设置误差的影响；将μ_k-1与P_e,k|k-1的先验分布建模为其共轭先验，即Gauss-inverse Gama，简写为GIG分布：

其中：表示μ_k-1与P_e,k|k-1的先验统计参数，GIG(a,A|τ,α,λ,ν)表示以τ,α,λ,ν为参数的GIG分布，可以被分解为：

GIG(a,A|τ,α,λ,ν)＝N(a|τ,αA)IG(A|λ,ν)

其中：IG(A|λ,ν)表示以λ,ν为参数的逆Gama分布；

μ_k-1与P_e,k|k-1的先验分布可以被分解为：

μ_k-1与P_e,k|k-1的先验估计需要与其名义值匹配，即μ_k-1与P_e,k|k-1的先验统计参数满足：

设计μ_k-1与P_e,k|k-1统计参数预测方程为：

其中：ρ_α与ρ_λ为调制参数；

E.水声信标周期性广播水声信号，水声信号发射时间及水声信标位置已知；所述水下航行器在未接收到水声信号时，通过自身配备的电子罗盘、深度计以及读取自身的螺旋桨转速信息进行航位推算，同时进行有效声速随机模型参数预测；所述水下航行器在接收到所搭载的多普勒测速仪测得的绝对速度观测后，通过读取所述螺旋桨转速信息及电子罗盘信息，构造海流速度观测量并通过Kalman滤波进行海流速度校正；

所述水下航行器利用自身配备的电子罗盘、深度计以及读取自身的螺旋桨转速信息进行航位推算的方法为：

根据Kalman滤波的预测环节，将系统状态的先验分布近似为Gauss分布：

其中：与P_k|k-1分别为k时刻的先验状态和先验方差，其计算方法为：

其中：与P_k-1|k-1分别为k-1时刻的后验状态和后验方差；

所述水下航行器接收到所述多普勒测速仪测得的绝对速度观测后，进行海流速度校正的方法为：

P_k|k＝P_k|k-1-K_kH_kP_k|k-1

其中：K_k为海流更新Kalman增益；

F.所述水下航行器接收到水声信号后，记录接收时刻，根据已知的水声信号发射时刻及水声信标位置坐标，并考虑水下声速及其不确定性统计参数的未知性，以此基于扩展Kalman滤波算法及变分贝叶斯近似，以水声信号传递时间为观测变量，进行所述水下航行器的位置更新；

S1.建立观测似然函数、有效声速后验分布及其不确定性统计参数后验概率密度分布模型；

根据离散形式的水声信号传递时间观测模型，得到水声信号传递时间的观测似然密度为：

p(m_k|v_e,k,x_k)＝N(m_k|h_k(x_k,v_e,k),R_k)

将有效声速后验分布建模为Gauss分布，即：

μ_k-1与P_e,k|k-1的后验分布同样建模为GIG分布，分解为：

S2.定义状态变量、有效声速及声速不确定性统计参数后验估计近似值；

通过变分贝叶斯近似，将状态变量、有效声速及声速不确定性统计参数的联合后验估计近似为：

p(x_k,v_e,k,μ_k-1,P_e,k|k-1|m_1:k)≈q(x_k)q(v_e,k)q(μ_k-1)q(P_e,k|k-1)

以最小化近似前与近似后两个概率密度函数之间的Kullback-Leibler散度(KLD)最小化为目标，得到近似解为：

其中：log(·)表示对数运算；θ表示x_k,v_e,k,μ_k-1,P_e,k|k-1中的任意元素；E_x[·]表示相对于x的期望；上标(-θ)表示整个集合当中除了θ以外的其它元素；c_θ表示与θ无关的常数；采用固定点迭代来求解q(θ)；

S3.求解状态变量、有效声速及声速不确定性统计参数后验估计近似值；

S3.1求解联合概率密度对数值；

其对数形式表示为：

S3.2求解声速不确定性统计参数P_e,k|k-1估计值；

令θ＝P_e,k|k-1，得到：

其中：

变量加上标(i)表示该变量在第i次迭代中的估计值；进而得到：

取P_e,k|k-1在第i+1次迭代中的估计值为：

S3.3求解声速不确定性统计参数μ_k-1估计值；

令θ＝μ_k-1，得到：

进而得到：

取μ_k-1在第i+1次迭代中的估计值为其期望值，即：

S3.4求解有效声速v_e,k估计值；

记选择为展开点，对非线性的水声信号传递时间观测方程m_k＝h_k(x_k,v_e,k)进行线性化，只保留一阶项，得到：

其中：

为观测方程的雅克比矩阵；

由此得到线性化后的水声信号传递时间观测方程为：

令θ＝v_e,k，得到：

其中：

根据Kalman滤波更新式，得到有效声速更新方程为：

其中：为第i+1次迭代中有效声速更新的Kalman增益；

S3.5求解系统状态x_k估计值；

以为展开点，对非线性的水声信号传递时间观测方程m_k＝h_k(x_k,v_e,k)进行线性化，只保留一阶项，得到：

其中：

为观测方程的雅克比矩阵；

由此得到线性化后的水声信号传递时间观测方程为：

令θ＝x_k，得到：

其中：

根据Kalman滤波更新式，得到系统状态更新方程为：

其中：为第i+1次迭代中系统状态更新的Kalman增益；

记总迭代次数为N，则最终有效声速均值及方差、系统状态均值及协方差矩阵的估计值分别为：

参数需要用于k+1时刻的预测过程，记其估计值为