本发明公开了一种基于水下地形高程数据库的水下航行器辅助导航定位方法，其主要目的在于解决水下航行器长时间航行时由于惯性导航定位系统的误差会随着时间的积累逐渐增大的固有缺陷以至于不能满足航行器长时间精确导航定位需要的问题。本发明的主要步骤包括：主惯导系统指示航迹序列和实测高程序列的获取、匹配定位预处理、一级匹配定位、二级匹配定位以及修正过程。本发明可以解决水下航行器长时间航行所带来的初始位置误差过大致使匹配失败问题以及匹配过程中由于惯性器件的误差所引起的匹配误差问题，从而修正主惯导系统，提高导航定位精度。
1.一种基于水下地形高程数据库的水下航行器辅助导航定位方法，其特征为：步骤1：获取主惯导系统指示航迹序列和实测高程序列
当水下航行器驶入匹配区时，在主惯性导航系统INS的导航下航行一段距离，与此同时每隔一个时间段Times，Times的典型值为10分钟，通过航行器的高程测量装置测量高程，由此得到惯性导航系统给出的一个惯导系统指示航迹序列 i＝1,2…N，N代表序列点的个数，取值为20， 代表惯导系统指示航迹序列Pi的纬度， 代表惯导系统指示航迹序列Pi的经度，和一个由高程测量装置给出的对应于惯导系统指示航迹序列Pi的实测高程序列Di,i＝1,2…N，考虑高程测量装置的误差可得：其中 表示理想高程序列，γ表示高程测量装置的测量噪声序列，服从均值为0标准差为1的正态分布，
在加载地形高程数据库中，选择以惯导系统指示航迹序列Pi为中心、以惯性导航系统6倍误差为边长的区域，将此区域的高程数据存入一个预处理二维矩阵，预处理二维矩阵行标为row,1≤row≤rowmax，rowmax为所述区域行标的最大值，预处理二维矩阵列标为col,1≤col≤colmax，colmax为所述区域列标的最大值，矩阵中的值存储为高程数据；所述区域的起始纬度为lati0，起始经度为longti0，矩阵行列交点row，col处的经度、纬度分别为：lati0+(row-1)×grid，longti0+(col-1)×grid，grid为地形数据库中相邻高程数据的位置间距，考虑地形高程数据库建库时测量和量化所产生的误差可得：其中row为预处理二维矩阵的行标，col为预处理二维矩阵的列标，DB(row,col)表示预处理二维矩阵中行标为row、列标为col的高程数据， 表示对应于DB(row,col)的理想高程数据，η表示地形高程数据库建库时的测量和量化噪声，服从均值为0标准差为1的正态分布；
将实测高程序列Di和在DB(row,col)中所提取的高程序列作相关运算，所述的相关运算采用均方差MSD算法，使用去均值后的数据进行相关运算：其中，Ral(row,col)表示MSD相关运算的结果，DB(row,col)表示惯导系统指示航迹Pi平移至预处理二维矩阵row，col对应位置lati0+(row-1)×grid，longti0+(col-1)×grid处高程数据序列，DBA(row,col)为DB(row,col)中所提取的高程序列的均值，Di表示实测高程序列，DA表示实测序列的均值；
运算结束后，取得使Ral(row,col)最小的row，col，则将惯导系统指示航迹Pi平移至预处理二维矩阵row，col对应位置lati0+(row-1)×grid，longti0+(col-1)×grid，得到预处理位置序列为 表示得到预处理位置序列Pi′的纬度， 表示得到预处理位置序列Pi′的经度；
将地形高程数据库中以预处理位置序列Pi′为中心，以惯性导航系统3倍误差为边长的一级匹配区域内的高程数据存入一个一级匹配二维矩阵，并利用双线性插值算法在该一级匹配区域提取高程为实测高程Di的一级等高线Ci；
惯导系统给出的航向误差为ψINS，纬度和经度误差分别为tINSlati和tINSlongti，分别在-3ψINS～3ψINS，-3tINSlati～3tINSlati，-3tINSlongti～3tINSlongti范围内分别随机取一个航向误差、一个纬度和一个经度且值分别为：ψINSrand，tINSlatirand，tINSlongtirand；对预处理位置序列Pi′作一级随机旋转和一级随机平移变换，得到一级随机旋转和平移位置序列P′irand：步骤3.2：
通过一级随机旋转和平移位置序列P′irand向一级等高线Ci作垂线，得到一个一级匹配垂足序列Yi，求取一个刚性变换Trans且所述刚性变换Trans包括旋转矩阵和平移矢量 φICCP为一级匹配旋
转角，tICCPLati代表纬度偏移量，tICCPLongti代表经度偏移量，使得一级随机旋转和平移位置序列P′irand逐渐逼近一级匹配垂足序列Yi，也就是使以下目标函数最小：其中Object为目标函数，iteration代表第iteration次迭代，iteration＝0表示迭代未开始，weighti代表序列点的权值，Trans为刚性变换，Dist(P′ irand,Yi)为P′irand和Yi之间的距离，Distmax为P′ irand和Yi序列点之间距离的最大值；
qua1为四元数quaternion的第一个元素，qua2为四元数quaternion的第二个元素，qua3为四元数quaternion的第三个元素，qua4为四元数quaternion的第四个元素，φICCP为一级匹配旋转角；
Mat11为矩阵MatICCP第一行第一列的元素，Mat12为矩阵MatICCP第一行第二列的元素，Mat21为矩阵MatICCP第二行第一列的元素，Mat22为矩阵MatICCP第二行第二列的元素， 和分别为一级匹配垂足序列的质心、随机旋转和平移位置序列的质心，weighti代表序列点的权值，T代表转置；
为求解旋转矩阵Rotation而构造一个对称的哈密尔顿矩阵：求解并得到哈密尔顿矩阵Hamilton的四个特征值Υ1,Υ2,Υ3,Υ4：Υ1＝[(Mat11+Mat22)2+(Mat21-Mat12)2]1/2Υ2＝-[(Mat11+Mat22)2+(Mat21-Mat12)2]1/2Υ3＝[(Mat11-Mat22)2+(Mat21+Mat12)2]1/2Υ4＝-[(Mat11-Mat22)2+(Mat21+Mat12)2]1/2记最大的特征值为Υm，最大的特征值为Υm可使下式成立：
(Mat11+Mat22-Υm)qua1+(Mat21-Mat12)qua4＝0由此可得一级匹配旋转角为：
φICCP＝2×arctan((Mat11+Mat22-Υm)/(Mat12-Mat21))所以旋转矩阵Rotation为：
利用步骤3.2得到的包括旋转矩阵Rotation和平移矢量translation的刚性变换Trans，对一级随机旋转和平移后的位置序列P′irand作P′ irand＝Trans·P′ irand处理，此时，若迭代次数iteration大于最大迭代次数iterationmax，表明收敛速度过慢，舍弃此次迭代的结果，返回步骤3.1；若iteration小于等于iterationmax并且|Objectiteration-Objectiteration-1|＞τ，则返回步骤3.2，若iteration小于等于iterationmax并且|Objectiteration-Objectiteration-1|≤τ，则迭代终止并记录Objectitetation和Trans·P′irand的值，iterationmax-6值为100，τ值为10 ；
重复3.1-3.3步骤，得到50个Objectitetation(Yi,Trans·P′irand)值后，找出其中最小的Objectitetation(Yi,Trans·P′irand)对应的Trans·P′irand，得到一级匹配位置序列LatiP″代表一级匹配位置序列的纬度，LontiP″代表一级匹配序列的经度；
选择以加载地形高程数据库中以一级匹配位置序列Pi″为中心、3×grid为边长的区域，将此区域的高程数据存入一个二级匹配二维矩阵，并利用双线性插值算法在该区域提取高程为实测高程Di的二级等高线Ci'，在二级等高线Ci'上寻找一级匹配位置序列Pi″的对应垂足得到一个二级垂足序列 代表二级垂足序列的纬度，代表二级垂足序列的经度；对二级垂足序列Yi′到一级匹配位置序列Pi″作仿射变换G且所述仿射变换包括在纬度方向上的平移因子taffinex，在经度方向上的平移因子taffiney，缩放因子αaffine，旋转因子θaffine，ei为二级垂足序列Yi′与二级匹配定位的输入序列Pi″的经仿射变换后的坐标的差：令仿射变换求解向量为taffine，仿射变换求解矩阵为
将使ESDaffine最小的taffine记为 对ESDaffine求导可得：令 可得
当ESDaffine对taffine的二阶导数 时， 就是使ESDaffine最小的解：
二级垂足序列的纬度之和， 代表一级匹配位置序列的经度的之和，代表二级垂足序列的经度之和， 代表
一级匹配位置序列和二级垂足序列的点积之和， 代表一
级匹配位置序列和二级垂足序列的叉积之和， 代表一级匹配位置序列的模的平方之和， 代表矩阵 的模；
其中 为求解所得的平移因子， 为求解所得的缩放因子， 为求解所得的旋转因子；
对一级匹配位置序列作仿射变换得二级匹配位置序列Pi″′：代表二级匹配位置序列的纬度， 代表二级匹配位置序列的经度；
步骤5：将二级匹配位置序列Pi″′中的最新位置信息P20″′作为观测信息，利用卡尔曼滤波器中的间接滤波法估计并修正主惯性导航系统位置误差，卡尔曼滤波器接收惯性导航系统的当前位置信息和最新位置信息P20″′的差值，经过滤波计算，估计出惯性导航系统位置误差量，用估计出的位置误差去修正惯性导航系统输出的位置信息，作为导航位置信息输出。
2.根据权利要求1所述的一种基于水下地形高程数据库的水下航行器辅助导航定位方法，惯性导航系统位置误差量的估计采用如下方法：
步骤5.1：建立惯性导航系统状态方程，选取东北天坐标系OENU作为导航坐标系，其中OE为东向轴，ON为北向轴，OU表示天向轴，载体坐标系Oxyz,其中Ox轴沿航行器横轴指向右舷，Oy轴沿航行器纵轴指向前，Oz轴垂直于Ox与Oy轴所确定的平面构成右手坐标系，选取状态向量XINS为：其中，δVE、δVN、δVU分别是东向、北向和天向速度误差；φE、φN、φU分别是东向、北向和天向失准角；δLati、δLongti、δh分别是纬度、经度和天向误差；▽bx、▽by、▽bz分别是Ox、Oy、Oz向的加速度计偏置；εbx、εby、εbz分别是Ox、Oy、Oz向的陀螺漂移；
其中F11为矩阵FINS的第1行第1列元素，F12为矩阵FINS的第1行第2列元素，F13为矩阵FINS的第1行第3列元素，F17为矩阵FINS的第1行第7列元素，F19为矩阵FINS的第1行第9列元素，F21为矩阵FINS的第2行第1列元素，F27为矩阵FINS的第2行第7列元素，F29为矩阵FINS的第2行第9列元素，F31为矩阵FINS的第3行第1列元素，F39为矩阵FINS的第
3行第9列元素，F45为矩阵FINS的第4行第5列元素，F46为矩阵FINS的第4行第6列元素，F54为矩阵FINS的第5行第4列元素，F64为矩阵FINS的第6行第4列元素，F67为矩阵FINS的第6行第7列元素，cij为姿态转移矩阵 的元素，FINS为系统状态转移矩阵；WINS为系统随机干扰向量，ωie为地球自转角速度，Radiusn、Radiuse分别为沿子午圈和卯酉圈的主曲率半径，VE、VN及VU为载体真实东向、北向和天向速度，Longti、Lati及h表示载体所处的真实经度、纬度及高程；卡尔曼滤波用离散化模型来描述系统，采用离散化后的差分方程描述连续系统，离散化后的系统状态方程为：其中Xk为k时刻的系统状态向量，Φk,k-1为k-1到k时刻的一步转移矩阵， 为系统状态向量的估计，Wk-1为k-1时刻的系统噪声，且
I为单位矩阵，ΔPer为将滤波周期Per分成σ个时间间隔ΔPer，Fj-1为FINS在时刻k-1+(j-1)ΔPer时的值；
步骤5.2：建立观测方程，将二级匹配位置序列Pi″′中的最新位置信息P20″′与惯性导航系统输出的经度LongtiINS、纬度LatiINS位置信息的差值，高程测量装置实测高程Di的最新实测高程D20与主惯导系统输出高程hINS的差值构成位置量测向量：其中： 为观测矩阵，υk为观测噪声；
T代表转置， 为对系统状态向量的估计，Pk为状态估计误差的协方差阵，Kk为滤波增益， 表示对量测值的一步预测， 为新息序列，Qk-1、Rk分别为Wk-1、υk的方差阵；给定初始条件X0和P0，由Zk即可得到每一步的状态最优估计步骤5.4利用当前卡尔曼滤波器得出的误差量最优估计向量 中的位置误差量，修正惯性导航系统位置信息，得到修正后的位置信息
一种基于水下地形高程数据库的水下航行器辅助导航定位
方法
技术领域
[0001] 本发明涉及水下导航技术领域，具体地说是设计一种能够满足水下潜器长时间航行要求的导航方法。
背景技术
[0002] 水下航行器在未来的军用和民用领域都日显其重要性，因而成为了世界军事和科技强国重点发展的领域。水下航行器要求具备高度的隐蔽性，能保持长时间的精确航行，而实现这一要求的关键技术之一是其导航技术。水下航行器的导航方式有别于陆路，GPS、无线电、天文导航等导航方式不能满足水下航行器隐蔽性的要求。惯性导航因其具有自主性、隐蔽性等优点，成为了水下潜器导航手段的首选。但惯性导航有固有的缺陷，惯性导航定位系统的误差会随着时间的积累逐渐增大以至于发散，因而在惯性导航定位系统运行一段时间后必须采用外来信息对惯性导航定位系统进行修正。随着导航技术特别是无源导航技术的发展，应用地形、重力、地磁等信息来辅助惯性导航为水下器导航研究的提供了一种新途径。由于地形高程信息测量较为方便，且陆路上的地形匹配技术已应用于实际，可为水下地形匹配导航提供借鉴，因而应用地形高程信息来辅助水下航行器的导航具备可行性。
发明内容
[0003] 发明目的：本发明提供一种基于水下地形高程数据库的能够实现长时间高精度导航定位的水下航行器辅助导航定位方法。
[0004] 本发明的技术方案具体如下：
[0005] 一种基于水下地形高程数据库的水下航行器辅助导航定位方法，其特征为：
[0006] 步骤1：获取主惯导系统指示航迹序列和实测高程序列
[0007] 当水下航行器驶入匹配区时，在主惯性导航系统INS的导航下航行一段距离，与此同时每隔一个时间段Times，Times的典型值为10分钟，通过航行器的高程测量装置测量高程，由此得到惯性导航系统给出的一个惯导系统指示航迹序列 i=1，2…N，N代表序列点的个数，取值为20， 代表惯导系统指示航迹序列Pi的纬度， 代表惯
导系统指示航迹序列Pi的经度，和一个由高程测量装置给出的对应于惯导系统指示航迹序列Pi的实测高程序列Di，i=1，2…N，考虑高程测量装置的误差可得：
[0008]
[0009] 其中 表示理想高程序列，γ表示高程测量装置的测量噪声序列，服从均值为0标准差为1的正态分布。
[0010] 步骤2：匹配定位预处理
[0011] 在加载地形高程数据库中，选择以惯导系统指示航迹序列Pi为中心、以惯性导航系统6倍误差为边长的区域，将此区域的高程数据存入一个预处理二维矩阵，预处理二维矩阵行标为row，1≤row≤rowmax，rowmax为所述区域行标的最大值，预处理二维矩阵列标为col，1≤col≤colmax，colmax为所述区域列标的最大值，矩阵中的值存为储高程数据；所述区域的起始纬度为lati0，起始经度为longti0，矩阵行列交点row，col处的经度、纬度分别为：lati0+(row-1)×grid，longti0+(col-1)×grid，grid为地形数据库中相邻高程数据的位置间距，考虑地形高程数据库建库时测量和量化所产生的误差可得：
[0012]
[0013] 其中row为预处理二维矩阵的行标，col为预处理二维矩阵的列标，DB(row，col)表示预处理二维矩阵中行标为row、列标为col的高程数据， 表示对应于DB(row，col)的理想高程数据，η表示地形高程数据库建库时的测量和量化噪声，服从均值为0标准差为1的正态分布；
[0014] 将实测高程序列Di和在DB(row，col)中所提取的高程序列作相关运算，所述的相关运算采用均方差MSD算法，使用去均值后的数据进行相关运算：
[0015]
[0016] 其中，Ral(row，col)表示MSD相关运算的结果，DB(row，col)表示惯导系统指示航迹Pi平移至预处理二维矩阵row，col对应位置lati0+(row-1)×grid，
longti0+(col-1)×grid处高程数据序列，DBA(row，col)为DB(row，col)中所提取的高程序列的均值，Di表示实测高程序列，DA表示实测序列的均值；
[0017] 运算结束后，取得使Ral(row，col)最小的row，col，则将惯导系统指示航迹Pi平移至预处理二维矩阵row，col对应位置lati0+(row-1)×grid，longti0+(col-1)×grid，得到预处理位置序列为 表示得到预处理位置序列P′i的纬度， 表示得到预处理位置序列P′i的经度；
[0018] 步骤3：一级匹配定位
[0019] 将地形高程数据库中以预处理位置序列P′i为中心，以惯性导航系统3倍误差为边长的一级匹配区域内的高程数据存入一个一级匹配二维矩阵，并利用双线性插值算法在该一级匹配区域提取高程为实测高程Di的一级等高线Ci；
[0020] 步骤3.1：
[0021] 惯导系统给出的航向误差为ψINS，纬度和经度误差分别为tINSlati和tINSlongti，分别在-3ψINS～3ψINS，-3tINSlati～3tINSlati，-3tINSlongti～3tINSlongti范围内分别随机取一个航向误差、一个纬度和一个经度且值分别为：ψINSrand，tINSlatirand，tINSlongtirand；对预处理位置序列P′i作一级随机旋转和一级随机平移变换，得到一级随机旋转和平移位置序列P′ irand：
[0022]
[0023] 步骤3.2：
[0024] 通过一级随机旋转和平移位置序列P′irand向一级等高线Ci作垂线得到一个一级匹配垂足序列Y，求取一个刚性变换Trans且所述刚性变换Trans包括旋转矩阵
[0025] 和平移矢量 φICCP为一级匹配旋转角，tICCPLati代表纬度偏移量，tICCPLongti代表经度偏移量，使得一级随机旋转和平移位置序列P′irand逐渐逼近一级匹配垂足序列Yi，也就是使以下目标函数最小：
[0026]
[0027]
[0028] 其中Object为目标函数，iteration代表第iteration次迭代，iteration=0表示迭代未开始，weighti代表序列点的权值，Trans为刚性变换，Dist(P′irand，Yi)为P′irand和Yi之间的距离，Distmax为P′ irand和Yi序列点之间距离的最大值；
[0029] 求取刚性变换Trans采用四元数法；
[0030] 设四元数为quaternion：
[0031]
[0032] qua1为四元数quaternion的第一个元素，qua2为四元数quaternion的第二个元素，qua3为四元数quaternion的第三个元素，qua4为四元数quaternion的第四个元素，φICCP为一级匹配旋转角；
[0033] 构造矩阵MatICCp：
[0034]
[0035] Mat11为矩阵MatICCP第一行第一列的元素，Mat12为矩阵MatICCP第一行第二列的元素，Mat21为矩阵MatICCP第二行第一列的元素，Mat22为矩阵MatICCP第二行第二列的元素，和 分别为一级匹配垂足序列的质心、随机旋转和平移位置序列的质心，weighti代表序列点的权值，T代表转置；
[0036] 为求解旋转矩阵Rotation而构造一个对称的哈密尔顿矩阵：
[0037]
[0038] 求解并得到哈密尔顿矩阵Hamilton的四个特征值
[0039]
[0040]
[0041]
[0042]
[0043] 记最大的特征值为 最大的特征值为 可使下式成立：
[0044]
[0045] 由此可得一级匹配旋转角为：
[0046]
[0047] 所以旋转矩阵Rotation为：
[0048]
[0049] 平移量translation为：
[0050]
[0051] 步骤3.3：
[0052] 利用步骤3.2得到的包括旋转矩阵Rotation和平移矢量translation的刚性变换Trans，对-级随机旋转和平移后的位置序列P′irand作P′ irand=Trans·P′irand处理，此时，若迭代次数iteration大于最大迭代次数iterationmax，表明收敛速度过慢，舍弃此次迭代的结果，返回步骤3.1；若iteration小于等于iterationmax并且
[0053] |Ogjectiteration-Ogjectiteration-1|＞τ，则返回步骤3.2，若iteration小于等于iterationmax并且|Objectiteration-Objectiteration-1|≤τ，则迭代终止并记录Objectiteration和Trans·P′irand的值，iterationmax值为100，τ值为10-6；
[0054] 步骤3.4：
[0055] 重复3.1-3.3步骤，得到50个Objectiteration 值后，找出其中最小的Objectitetation 对应的 得到一级匹配位置序列
代表一级匹配位置序列的纬度， 代表一级匹配序
列的经度；
[0056] 步骤4：二级匹配定位
[0057] 选择以加载地形高程数据库中以一级匹配位置序列P″i为中心、3×grid为边长的区域，将此区域的高程数据存入一个二级匹配二维矩阵，并利用双线性插值算法在该区域提取高程为实测高程Di的二级等高线C′ i，在二级等高线C′i上寻找一级匹配位置序列P″i的对应垂足得到一个二级垂足序列 代表二级垂足序列的纬度，代表二级垂足序列的经度；对二级垂足序列Y′i到一级匹配位置序列P″ i作仿
射变换G且所述仿射变换包括在纬度方向上的平移因子taffine，在经度方向上的平移因子taffine，缩放因子αafine’旋转因子θaffine，ei为二级垂足序列Y′ i与二级匹配定位的输入序列P″i的经仿射变换后的坐标的差：
[0058]
[0059] 令仿射变换求解向量为taffine，仿射变换求解矩阵为
[0060]
[0061]
[0062] 则
[0063] 欧式平方距离ESDaffine为：
[0064]
[0065] T代表转置；
[0066] 将使ESDaffine最小的taffine为 对ESDaffine求导可得：
[0067]
[0068] 令 可得
[0069]
[0070] 当ESDaffine对taffine的二阶导数 时， 就是使ESDaffine最小的解：
[0071]
[0072] 其中： 代表一级匹配位置序列的纬度的之和， 代表二级垂足序列的纬度之和， 代表一级匹配位置序列的经度的之和，
代表二级垂足序列的经度之和， 代表一
级匹配位置序列和二级垂足序列的点积之和， 代表一级
匹配位置序列和二级垂足序列的叉积之和， 代表一级匹配位置序
列的模的平方之和， 代表矩阵 的模；
[0073] 由此就可以得到：
[0074]
[0075] 其中 为求解所得的平移因子， 为求解所得的缩放因子， 为求解所得的旋转因子；
[0076] 对一级匹配位置序列作仿射变换得二级匹配位置序列P″′i：
[0077]
[0078] 代表二级匹配位置序列的纬度， 代表二级匹配位置序列的经度；
[0079] 步骤5：将二级匹配位置序列P″′i中的最新位置信息P″′20作为观测信息，利用卡尔曼滤波器中的间接滤波法估计并修正主惯性导航系统位置误差，卡尔曼滤波器接收惯性导航系统的当前位置信息和最新位置信息P″′20的差值，经过滤波计算，估计出惯性导航系统位置误差量，用估计出的位置误差去修正惯性导航系统输出的位置信息，作为导航位置信息输出。
[0080] 与现有技术相比，本发明具有如下优点：
[0081] 1)本发明对主惯导系统进行了匹配定位预处理，匹配定位预处理所采用的算法在大初始位置误差情况下定位速度快，但定位精度较低，一级匹配所采用的算法则是在匹配定位预处理的定位结果上进行匹配，具有较小的搜索范围，因而能加快定位速度并提高定位精度；
[0082] 2)本发明二级匹配定位对主惯导系统的定位误差建模采用了仿射模型，仿射模型中引入了旋转因子和缩放因子，可更加真实地反映航迹匹配过程中惯性仪器陀螺仪和加速度计的误差特性，因而比现有的匹配算法有更高的精度；
[0083] 3)本发明采用两级匹配定位方式，在求解仿射参数过程中，由于一级匹配定位的航迹已在理想航迹的附近，因此并不需要采用复杂的最优搜索算法，可直接采用简单的最小二乘法来求解仿射参数，这种求解方式大幅降低了算法的复杂度，同时从实施例的结果得出本发明的匹配方法精度比现有匹配方法精度大约提高20％。
附图说明
[0084] 图1为本发明所描述的地形辅助导航定位的原理图；
[0085] 图2为双线性插值算法示意图；
[0086] 图3为等高线提取方法示意图；
[0087] 图4为本发明所描述的仿射变换示意图；
[0088] 图5为实施例中所采用的地形图；
[0089] 图6为实施例中TERCOM/ICCP算法匹配和本发明的匹配结果图；
[0090] 图7为实施例中TERCOM／ICCP算法匹配和本发明匹配误差的对比图。
具体实施方式
[0091] 下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域普通技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。
[0092] 一种基于水下地形高程数据库的水下航行器辅助导航定位方法：
[0093] 步骤1：获取主惯导系统指示航迹序列和实测高程序列
[0094] 当水下航行器驶入匹配区时，在主惯性导航系统INS的导航下航行一段距离，与此同时每隔一个时间段Times，Times的典型值为10分钟，通过航行器的高程测量装置测量高程，由此得到惯性导航系统给出的一个惯导系统指示航迹序列
[0095] i＝1，2…N，N代表序列点的个数，取值为20， 代表惯导系统指示航迹序列Pi的纬度， 代表惯导系统指示航迹序列Pi的经度，和一个由高程测量装置给出的对应于惯导系统指示航迹序列Pi的实测高程序列D，i＝1，2…N，考虑高程测量装置的误差可得：
[0096]
[0097] 其中 表示理想高程序列，γ表示高程测量装置的测量噪声序列，服从均值为0标准差为1的正态分布。
[0098] 步骤2：匹配定位预处理
[0099] 在加载地形高程数据库中，选择以惯导系统指示航迹序列Pi为中心、以惯性导航系统6倍误差为边长的区域，将此区域的高程数据存入一个预处理二维矩阵，预处理二维矩阵行标为row，1≤row≤rowmax，rowmax为所述区域行标的最大值，预处理二维矩阵列标为col，1≤col≤colmax，colmax为所述区域列标的最大值，矩阵中的值存为储高程数据；所述区域的起始纬度为lati0，起始经度为longti0，矩阵行列交点row，col处的经度、纬度分别为：lati0+(row-1)×grid，longti0+(col-1)×grid，grid为地形数据库中相邻高程数据的位置间距，考虑地形高程数据库建库时测量和量化所产生的误差可得：
[0100]
[0101] 其中row为预处理二维矩阵的行标，col为预处理二维矩阵的列标，DB(row，col)表示预处理二维矩阵中行标为row、列标为col的高程数据， 表示对应于DB(row，col)的理想高程数据，η表示地形高程数据库建库时的测量和量化噪声，服从均值为0标准差为1的正态分布；
[0102] 将实测高程序列Di和在DB(row，col)中所提取的高程序列作相关运算，所述的相关运算采用均方差MSD算法，使用去均值后的数据进行相关运算：
[0103]
[0104] 其中，Ral(row，col)表示MSD相关运算的结果，DB(row，col)表示惯导系统指示航迹Pi平移至预处理二维矩阵row，col对应位置lati0+(row-1)×grid，
longti0+(col-1)×grid处高程数据序列，DBA(row，col)为DB(row，col)中所提取的高程序列的均值，Di表示实测高程序列，DA表示实测序列的均值；
[0105] 运算结束后，取得使Ral(row，col)最小的row，col，则将惯导系统指示航迹Pi平移至预处理二维矩阵row，col对应位置lati0+(row-1)×grid，longti0+(col-1)×grid，得到预处理位置序列为 表示得到预处理位置序列P′i的纬度， 表示得到预处理位置序列P′i的经度；
[0106] 步骤3：一级匹配定位
[0107] 将地形高程数据库中以预处理位置序列P′i为中心，以惯性导航系统3倍误差为边长的一级匹配区域内的高程数据存入一个一级匹配二维矩阵，并利用双线性插值算法在该一级匹配区域提取高程为实测高程Di的一级等高线Ci；
[0108] 双线性插值算法示意图如图2所示，其计算方法为：
[0109] 坐标轴Oposx和坐标轴Oposy构成平面直角坐标系，坐标轴Oele垂直于坐标轴Oposx和坐标轴Oposy构成平面直角坐标系构成右手坐标系在坐标posy方向：
[0110]
[0111] scale2=posm--pos2
[0112] 在坐标posx方向：
[0113]
[0114] scale1=posn-pos1
[0115] posn，posm的高程值为：
[0116] 其中pos1，pos2为高程数据库中的网格点，scale2代表坐标posm方向上的比例，scale1代表坐标posn方向上的比例，posn，posm为pos1，pos2、pos1+1，pos2、pos1+1，pos2+1、pos1，pos2+1四个网格点所围成的矩形区域内的任意一点，ele(pos1，posm)表示pos1，posm处的高程值，ele(pos1+1，posm)表示pos1+1，posm处的高程值，ele(posn，pos2)表示posn，pos2处的高程值，ele(posn，pos2+1)表示posn，pos2+1处的高程值，ele(posn，posm)表示posn，posm处的高程值，ele(pos1，pos2)表示pos1，pos2处的高程值，ele(pos1+1，pos2)表示pos1+1，pos2处的高程值，ele(pos1，pos2+1)表示pos1，pos2+1处的高程值，ele(pos1+1，pos2+1)表示pos1+1，pos2+1处的高程值；
[0117] 提取等高线的方法为：
[0118] 设实测高程为Dtemp，四个网格点中部分高程小于Dtemp，而其它大于Dtemp时等高线穿过该网格单元，当等高线穿过该网格单元时，四个网格点所组成的矩形的四条边上高程为Dtemp的点有2个点或者是4个点；当是2个点的情况时，将两点所连直线近似为该网格单元内高程为Dtemp的等高线；当是如图3所示4个点的情况时，网格的四条边上都有出入点，所以首先必须得确定成对的出入点。
[0119] 令标记高程eleflag为：
[0120]如图3
所示，pA，pB，pC，pD为网格的四条边上出入点，当Dtemp＞eleflag时pA，pB成对，pC，pD成对；当Dtemp＜eleflag时pA，pD成对，pB，pC成对；当Dtemp=eleflag时pA，pC成对，pB，pD成对；两条成对点所连直线近似为该网格单元内高程为Dtemp的等高线。
[0121] 步骤3.1：
[0122] 惯导系统给出的航向误差为ψINS，纬度和经度误差分别为tINSlati和tINSlongti，分别在-3ψINS～3ψINS，-3tINSlati～3tINSlati，-3tINSlongti～3tINSlongti范围内分别随机取一个航向误差、一个纬度和一个经度且值分别为：ψINSrand，tINSlatirand，tINSlongtirand；对预处理位置序列P′i作一级随机旋转和一级随机平移变换，得到一级随机旋转和平移位置序列P′ irand：
[0123]
[0124] 步骤3.2：
[0125] 通过一级随机旋转和平移位置序列P′irand向一级等高线Ci作垂线得到一个一级匹配垂足序列Yi。
[0126] 获取垂足序列的方法为：由等高线的提取方法可知等高线是由一系列线段组成。求一点到线段的垂足有三种情况：第一种是该点落在线段上，则该点就是垂足。第二种是由该点向线段作垂线，且垂线和线段相交在线段之内，则相交点就是垂足。第三种是由该点向线段作垂线，但垂线和线段的延长线相交，或者该点落在线段的延长线上，在这种情况下，本发明选取线段的两个端点中距离该点较近的点作为垂足。由于等高线由一系列线段组成，则一级随机旋转和平移位置序列P′irand中的点在其对应等高线上可以求得不止一个的垂足，本发明中采用的策略是在垂足中选取与一级随机旋转和平移位置序列P′irand中的点的距离最近点组成一级匹配垂足序列Yi。
[0127] 求取一个刚性变换Trans且所述刚性变换Trans包括旋转矩阵
[0128] 和平移矢量 φICCP为一级匹配旋转角，tICCPLati代表纬度偏移量，tICCPLongti代表经度偏移量，使得一级随机旋转和平移位置序列P′irand逐渐逼近一级匹配垂足序列Yi，也就是使以下目标函数最小：
[0129]
[0130]
[0131] 其中Object为目标函数，iteration代表第iteration次迭代，iteration=0表示迭代未开始，weighti代表序列点的权值，Trans为刚性变换，Dist(P′irand，Yi)为P′irand和Yi之间的距离，Distmax为P′ irand和Yi序列点之间距离的最大值；
[0132] 求取刚性变换Trans采用四元数法；
[0133] 设四元数为quaternion：
[0134]
[0135] qua1为四元数quaternion的第一个元素，qua2为四元数quaternion的第二个元素，qua3为四元数quaternion的第三个元素，qua4为四元数quaternion的第四个元素，φICCP为一级匹配旋转角；
[0136] 构造矩阵MatICCP：
[0137]
[0138] Mat11为矩阵MatICCP第一行第一列的元素，Mat12为矩阵MatICCP第一行第二列的元素，Mat21为矩阵MatICCP第二行第一列的元素，Mat22为矩阵MatICCP第二行第二列的元素，和 分别为一级匹配垂足序列的质心、随机旋转和平移位置序列的质心，weighti代表序列点的权值，T代表转置；
[0139] 为求解旋转矩阵Rotation而构造一个对称的哈密尔顿矩阵：
[0140]
[0141] 求解并得到哈密尔顿矩阵Hamilton的四个特征值 ：
[0142]
[0143]
[0144]
[0145]
[0146] 记最大的特征值为 最大的特征值为 可使下式成立：
[0147]
[0148] 由此可得一级匹配旋转角为：
[0149]
[0150] 所以旋转矩阵Rotation为：
[0151]
[0152] 平移量translation为：
[0153]
[0154] 步骤3.3：
[0155] 利用步骤3.2得到的包括旋转矩阵Rotation和平移矢量translation的刚性变换Trans，对一级随机旋转和平移后的位置序列P′irand作P′ irand=Trans·P′irand处理，此时，若迭代次数iteration大于最大迭代次数iterationmax，表明收敛速度过慢，舍弃此次迭代的结果，返回步骤3.1；若iteration小于等于iterationmax并且|O勿ectiteration-Objectiteration-1|＞τ，则返回步骤3.2，若iteration小于等于iterationmax并且|Objectiteration-Objectiteration-1|≤τ，则迭代终止并记录Objectitetation和Trans·P′irand的值，iterationmax-6值为100，τ值为10 ；
[0156] 步骤3.4：
[0157] 重复3.1-3.3步骤，得到50个Objectitetation(Yi，Trans·P′irand)值后，找出其中最小的Objectitetation(Yi，Trans·P′irand)对应的Trans·P′irand，得到一级匹配位置序列代表一级匹配位置序列的纬度， 代表一级匹配序列的经度；
[0158] 步骤4：二级匹配定位
[0159] 加载地形高程数据库中以一级匹配位置序列P″i为中心，3×grid为边长区域，将此区域的高程数据存入一个二级匹配二维矩阵，并利用双线性插值算法在该区域提取高程为实测高程Di的二级等高线Ci′，在二级等高线Ci′上寻找一级匹配位置序列P″i的对应垂足得到一个二级垂足序列 代表二级垂足序列的纬度， 代表二级垂足序列的经度；对二级垂足序列Y′i到一级匹配位置序列P”i作仿射变换G且所述仿射变换包括在纬度方向上的平移因子taffinx，在经度方向上的平移因子taffinx，缩放因子αaffine，旋转因子θaffine，ei为二级垂足序列Y′ i与二级匹配定位的输入序列P″ i的经仿射变换后的坐标的差：
[0160]
[0161] 令仿射变换求解向量为taffine，仿射变换求解矩阵为
[0162]
[0163]
[0164] 则
[0165] 则欧式平方距离ESDaffine为：
[0166]
[0167] T代表转置；
[0168] 将使ESDDaffine最小的taffine记为 对ESDaffine求导可得：
[0169]
[0170] 令 可得
[0171]
[0172] 当ESDaffine对taffine的二阶导数 时， 就是使ESDaffine最小的解：
[0173]
[0174] 其中： 代表一级匹配位置序列的纬度的之和， 代表二级垂足序列的纬度之和， 代表一级匹配位置序列的经度的之和，
代表二级垂足序列的经度之和， 代表一
级匹配位置序列和二级垂足序列的点积之和， 代表一级
匹配位置序列和二级垂足序列的叉积之和， 代表一级匹配位置序
列的模的平方之和， 代表矩阵 的模；
[0175] 由此就可以得到
[0176]
[0177] 其中 为求解所得的平移因子， 为求解所得的缩放因子， 为求解所得的旋转因子；
[0178] 对一级匹配位置序列作仿射变换得二级匹配位置序列P″′i：
[0179]
[0180] 代表二级匹配位置序列的纬度， 代表二级匹配位置序列的经度；
[0181] 步骤5：将二级匹配位置序列P″′i中的最新位置信息P″′20作为观测信息，利用卡尔曼滤波器中的间接滤波法估计并修正主惯性导航系统位置误差，卡尔曼滤波器接收惯性导航系统的当前位置信息和最新位置信息P″′20的差值，经过滤波计算，估计出惯性导航系统位置误差量，用估计出的位置误差去修正惯性导航系统输出的位置信息，作为导航位置信息输出。
[0182] 惯性导航系统位置误差估计修正步骤如下：
[0183] 步骤5.1：建立惯性导航系统状态方程，选取东北天坐标系OENU作为导航坐标系，其中OE为东向轴，ON为北向轴，OU表示天向轴，载体坐标系Oxyz，其中Ox轴沿航行器横轴指向右舷，Oy轴沿航行器纵轴指向前，Oz轴垂直于Ox与Oy轴所确定的平面构成右手坐标系，选取状态向量XINS为：
[0184]
[0185] 其中，δVE、δVN、δVU分别是东向、北向和天向速度误差；φE、φN、φU分别是东向、北向和天向失准角；δLati、δLongti、δh分别是纬度、经度和天向误差；Δbx、Δby、Δbz分别是Ox、Oy、Oz向的加速度计偏置；εbx、εby、εbz分别是Ox、Oy、Oz向的陀螺漂移；
[0186] 建立惯性导航系统线性状态方程为：
[0187]
[0188] 其中：
[0189]
[0190]
[0191]
[0192]
[0193]
[0194]
[0195]
[0196]
[0197]
[0198]
[0199] 其中F11为矩阵FINS的第1行第1列元素，F12为矩阵FINS的第1行第2列元素，F13为矩阵FINS的第1行第3列元素，F17为矩阵FIVS的第1行第7列元素，F19为矩阵FINS的第1行第9列元素，F21为矩阵FIVS的第2行第1列元素，F27为矩阵FINS的第2行第7列元素，F29为矩阵FINS的第2行第9列元素，F31为矩阵FINS的第3行第1列元素，F39为矩阵FINS的第3行第9列元素，F45为矩阵FINS的第4行第5列元素，F46为矩阵FINS的第4行第6列元素，F54为矩阵FINS的第5行第4列元素，F64为矩阵FINS的第6行第4列元素，F67为矩阵FINS的第6行第7列元素，cij为姿态转移矩阵 的元素，FINS为系统状态转移矩阵；WINS为系统随机干扰向量，ωie为地球自转角速度，Radiusn、Radiuse分别为沿子午圈和卯酉圈的主曲率半径，VE、VN、及VU为载体真实东向、北向和天向速度，Longti、Lati及h表示载体所处的真实经度、纬度及高程；卡尔曼滤波用离散化模型来描述系统，采用离散化后的差分方程描述连续系统，离散化后的系统状态方程为：
[0200]
[0201] 其中Xk为k时刻的系统状态向量，Φk，k-1为k-1到k时刻的一步转移矩阵， 为系统状态向量的估计，Wk-1为k-1时刻的系统噪声，且
[0202]
[0203] I为单位矩阵，ΔPer为将滤波周期Per分成σ个时间间隔ΔPer，Fj-1为FINS在时刻k-1+(j-1)ΔPer时的值；
[0204] 步骤5.2：建立观测方程，将二级匹配位置序列P″′i中的最新位置信息P″′ 20与惯性导航系统输出的经度LongtiINS、纬度LatiINS位置信息的差值，高程测量装置实测高程Di的最新实测高程D20与主惯导系统输出高程hINS的差值构成位置量测向量：
[0205]
[0206] 其中： 为观测矩阵，υk为观测噪声；
[0207] 步骤5.3：卡尔曼滤波器误差估计：
[0208] 步骤5.3.1时间更新：
[0209]
[0210] 步骤5.3.1量测更新：
[0211]
[0212] T代表转置， 为对系统状态向量的估计，Pk为状态估计误差的协方差阵，Kk为滤波增益， 表示对量测值的一步预测， 为新息序列，Qk-1、Rk分别为Wk-1、uk的方差阵；给定初始条件X0和P0，由Zk即可得到每一步的状态最优估计
[0213] 步骤5.4利用当前卡尔曼滤波器得出的误差量最优估计向量 中的位置误差量，修正惯性导航系统位置信息，得到修正后的位置信息
[0214]
[0215] 其中 分别是卡尔曼滤波器输出的纬度、经度误差。
[0216] 本发明所描述的地形辅助导航定位的原理图如图1所示；
[0217] 实施例：
[0218] 仿真实验中地形高程数据库的范围为(北纬38.0°，东经120.0°)至(北纬38.5°，东经120.5°)的矩形区域，数据库地形的精度为0.075′×0.075′，在该区域内设置一条理想航迹位置序列，位置序列点数为20点，起始位置点为(北纬38.2°，东经
120.2°)，初始航向为38.7°，航向变化范围为±7°，初始航速为11.5节，航速变换范围为±1.3节。水下航行器惯导系统的陀螺仪漂移为0.02°/h，加速度计零偏为0.0005m／
2
s，初始位置误差在经度方向取为0.75′，纬度方向取为0.15′，采样时间为30s，测深装置的测量误差取为幅值为0.2的白噪声。表1为现有的ICCP迭代最近等值线匹配定位算法，现有的TERCOM／ICCP地形轮廓／迭代最近点匹配定位算法和本发明的匹配定位仿真的统计结果：
[0219] 表1：
[0220]参数 参数值 单位
ICCP匹配纬度误差均值误差 436.2 m
现有的TERCOM／ICCP匹配纬度均值误差 69.2 m
本发明的匹配纬度均值误差 47.6 m
ICCP匹配经度误差均值误差 448.1 m
现有的TERCOM／ICCP匹配经度均值误差 50.1 m
本发明匹配经度均值误差 28.9 m
ICCP匹配时间 203.5 s
现有的TERCOM／ICCP匹配时间 7.568 s
本发明匹配时间 7.594 s
[0221] 从仿真结果来可以看出，在较大的初始位置误差下，ICCP匹配的误差很大，属于误匹配，这是由于初始位置过大而不满足ICCP匹配算法正确执行所需满足的初始航迹在理想航迹附近的假设前提。由于TERCOM算法对初始位置不敏感，因而采用TERCOM／ICCP的方法进行匹配时，可以得到较为准确的匹配结果。而采用本发明的方法，在匹配时间相当的情况下，相比于TERCOM／ICCP匹配方法拥有更高的匹配精度。因此，本发明可以解决水下航行器长时间潜航所产生大初始位置情况下易出现误匹配的问题，同时以极小时间代价提高了导航定位的精度。
