CN103235327B - 一种gnss/mins超深组合导航方法、***及装置 - Google Patents

Abstract

为了实现GNSS和MINS的超深层次组合，本发明提出一种GNSS/MINS超深组合导航方法，利用INS提供的位置、速度辅助信息，GNSS使用直接位置估计算法估计位置误差、速度误差、时钟误差，以此闭合环路，***直接以MIMU、基带相关器为传感器，在顶层组合算法中实现所有的导航功能，即该算法将GNSS基带相关器作为敏感空时定位场的传感器，将MIMU作为敏感惯性场的传感器，实现了GNSS与MINS的完整融合，在INS辅助下，GNSS联合多通道信息进行矢量鉴相和矢量定位。本发明还提出了应用于此方法的超深组合***和装置。本发明具有导航精度高、动态性能好、GNSS跟踪灵敏度高、GNSS抗干扰性能强、GNSS动态牵引范围广等优点，理论上本发明动态范围受MINS动态限制，MINS辅助下GNSS能实现-160dBm信号跟踪。

Description

一种GNSS/MINS超深组合导航方法、***及装置

技术领域

本发明属于组合导航技术领域，尤其涉及一种GNSS/MINS超深组合导航方法、***及装置。

背景技术

GNSS/MINS(Global Navigation Satellite System全球卫星导航***/Micro-Electro-Mechanical Systems Inertial Navigation System微机电惯性导航***)组合导航***是指以MIMU(MEMS Inertial Measurement Unit微机电惯性测量单元)作为传感器，由GNSS接收机输出位置/速度等参考信息校正MINS误差，由MINS输出位置/速度辅助GNSS接收机捕获跟踪的导航***。

GNSS/MINS组合导航***因其成本低、体积小、误差不随时间累积、高灵敏度、抗高动态应力等优点，在未来航空、车载应用中有着良好的前景。目前组合导航算法有松组合、紧组合、深组合三种，分别基于位置/速度、伪距/伪距率、I/Q累加值层面进行组合。

目前为止，无论松组合、紧组合还是深组合，GNSS和MINS两者之间相对独立。MINS以MIMU为传感器，附加捷联解算实现惯性导航功能；GNSS以基带相关器为传感器，附加导航解算算法实现卫星导航功能；GNSS和MINS再在顶层组合模块中实现融合。

用于组合导航的INS一般采用捷联解算算法。GNSS环路跟踪可分为标量环路和矢量环路(VDFLL)。标量环路中使用独立的DLL(延迟锁定环)和PLL(相位锁定环)跟踪每颗卫星，每颗卫星的观测值在输入导航解算模块，使用卡尔曼滤波器或递推最小二乘法进行PVT(位置、速度和时间)解算；矢量算法首先根据预测位置和速度计算卫星伪距、伪距率，再根据伪距、伪距率推算载波频率和码相位，鉴相器输出作为观测量，代表预测伪距、伪距率与真实伪距、伪距率之差，反馈到卡尔曼滤波器对用户位置再进行预测，从而闭合环路。可见，标量环路跟踪和解算分开进行，而矢量环路跟踪和解算在一个卡尔曼滤波器中完成。

标量算法各通道独立工作，无法利用各通道之间的冗余信息，因此跟踪灵敏度低且动态性能不好，但优点是各通道不会相互影响。

尽管矢量环路耦合多个通道信息进行跟踪和定位，其跟踪灵敏度优于标量环路，但矢量环路中使用的鉴相器为非线性鉴相器，且随载噪比下降，鉴相误差非线性增长，导致跟踪和定位性能下降，因此矢量环路不适于跟踪GNSS弱信号。

发明内容

本发明旨在至少在一定程度上解决上述技术问题之一或至少提供一种有用的商业选择。为此，本发明的第一目的在于提出一种具有融合程度好、导航精度高的GNSS/MINS超深组合导航方法。

根据本发明实施例的GNSS/MINS超深组合导航方法，以MIMU和GNSS基带相关器为传感器，进行组合导航计算，包括：步骤S01，组合导航***初始化参数设置；步骤S02，GNSS捕获卫星并获取导航电文，INS在GNSS辅助下完成大失准角初始对准；步骤S03，根据INS输出的位置和速度以及GNSS输出的初始钟偏和初始钟漂，进行第一网格划分，并完成从位置、速度、钟偏和钟漂到载相频率和码相位的映射，实现第二网格划分，对应每个第二网格中心点建立相关器；步骤S04，进行1比特电文长度的相干积分，相关器输出相关值；步骤S05，判断是否达到积分时间T，若未达到，进入步骤S06，反之进入步骤S07；步骤S06，***进行非相干累加以及捷联导航解算，并根据所述捷联导航解算结果实时调整所述第二网格划分中心位置，返回步骤S05；步骤S07，将每颗卫星的非相干累加相关值进行第二网格到第一网格的逆映射，然后进行最大值的搜索，获得T时间内的位置误差、速度误差、钟偏误差和钟漂误差；步骤S08，将步骤S07中得到的各误差作为组合卡尔曼滤波器的观测量，进行INS误差更新及时钟误差更新，然后将卡尔曼滤波器状态清零；步骤S09，搜索网格大小自适应调整；步骤S10，捷联导航解算更新用户位置、速度以后，进入下一算法循环。

在本发明的方法的一个实施例中，所述第一网格划分为在三维位置、三维速度、钟漂、钟偏上的八维网格划分，第二网格划分为在GNSS载波频率、伪码相位上的二维网格划分。

在本发明的方法的一个实施例中，所述超深组合导航方法利用气压高度表引入气压高度数据辅助导航。

在本发明的方法的一个实施例中，所述超深组合导航方法应用高精度晶振。

在本发明的方法的一个实施例中，所述超深组合导航方法使用最大似然估计方法，最大似然估计式为：

其中，是待估计的接收机位置、钟偏、速度和钟漂，r(t)是GNSS接收基带信号加噪声，s_k(t)S_k(t)是接收机产生的第k颗卫星基带信号，N是可视卫星数量，T₀是预检测积分起始时刻，T是预检测积分时间，t是接收机生成基带信号时间序列。假设S_k(t)跨越D比特导航电文，典型情况下D取值50～150；

在第n数据位，第k颗卫星本地基带信号s_k(t)为：

s_k(t)＝a_k,nC_k(t)exp[j(ω_kt+φ_k,n)]

其中a_k,n为信号幅度，与载噪比有关；C_k(t)为CA码，由用户位置、钟偏决定；ω_k为载波圆频率，由用户速度、钟漂决定；φ_k,n为载波相位；

将最大似然估计式展开，得到：

由于R只与接收信号有关，因此最小化[R-S]，等效于最大化[S]，即：

其中I_k,n是第k颗卫星第n数据位的时间跨度，内部求和号为非相干累加，外部求和号为各卫星通道融合累加；

CA码序列波形m_k(t)求解式为：

$C_k\left(t\right)\≅C_k[t-b\left(T_0\right)-{\mathrm{\τ}}_k\left(T_0\right)-(t-T_0){\mathrm{\τ}}_k^{\′}\left(T_0\right)]$

其中C_k(t_GNSS)为以GNSS时间为参数的卫星CA码产生函数，C_k(t_GNSS)已知；τ_k(T₀)为T₀到达的第k颗卫星信号传输延时；τ'_k(T₀)为由多普勒频率引起的延迟时间变化率；b(T₀)为钟偏，由导航电文给出，几个参数计算如下：

${\mathrm{\τ}}_k\left(T_0\right)\≅\frac{|p_k\left(T_0\right)-p\left(T_0\right)|}{c+v_k\left(T_0\right)\·u_k\left(T_0\right)}\sec$

${\mathrm{\τ}}_k^{\′}\left(T_0\right)\≅\frac{1}{c}\left\{\right[v_k\left(T_0\right)-v\left(T_0\right)]\·u_k\left(T_0\right)\}\sec /\sec$

其中p_k(T₀)为第k颗卫星位置；p(T₀)为用户位置；v_k(T₀)为第k颗卫星速度；u_k(T₀)为用户到第k颗卫星的投影矢量；v(T₀)为用户速度；c为光速；

基带多普勒频率ω_k计算式为：

${\mathrm{\ω}}_k=2\mathrm{\π}\left\{\frac{1}{\mathrm{\λ}}\right[v\left(T_0\right)-v_k\left(T_0\right)]\·u_k\left(T_0\right)-\stackrel{\·}{b}\left(T_0\right)f\}$

其中f为载波频率，λ为载波波长。可见，已知导航电文、用户速度、钟漂的情况下，可以求解基带多普勒频率ω_k。

在本发明的方法的一个实施例中，在步骤S07中，按照水平速度、垂直速度、钟漂、水平位置、高度、钟偏的顺序进行最大值的搜索，获得对应T时间内的水平速度误差、垂直速度误差、钟漂误差、水平位置误差、高度误差、钟偏误差。

在本发明的方法的一个实施例中，所述水平速度误差搜索方法为：假设第n数据位起始时刻INS输出位置、速度为钟偏和钟漂为以这组参数为中心进行8维网格划分，每个维度划分网格的数量依次为n₁～n₈，且n₁～n₈均为奇数，针对每颗卫星，在码相位、载波多普勒两个维度均分为9份：

$\left[\begin{array}{ccccc}{\tilde{p}}_{n,x}-\frac{n_1-1}{2}\mathrm{\Δ}{p}_x& \·\·\·& {\tilde{p}}_{n,x}& \·\·\·& {\tilde{p}}_{n,x}+\frac{n_1-1}{2}\mathrm{\Δ}{p}_x\\ {\tilde{p}}_{n,y}-\frac{n_2-1}{2}\mathrm{\Δ}{p}_y& \·\·\·& {\tilde{p}}_{n,y}& \·\·\·& {\tilde{p}}_{n,y}+\frac{n_2-1}{2}\mathrm{\Δ}{p}_y\\ {\tilde{p}}_{n,z}-\frac{n_3-1}{2}\mathrm{\Δ}{p}_z& \·\·\·& {\tilde{p}}_{n,z}& \·\·\·& {\tilde{p}}_{n,z}+\frac{n_3-1}{2}\mathrm{\Δ}{p}_z\\ {\tilde{b}}_n-\frac{n_7-1}{2}\mathrm{\Δb}& \·\·\·& {\tilde{b}}_n& \·\·\·& {\tilde{b}}_n+\frac{n_7-1}{2}\mathrm{\Δb}\end{array}\right]\⇒\left[\begin{array}{ccccc}{m}_1-4\mathrm{\Δ}{m}_1& \·\·\·& m_1& \·\·\·& m_1+4\mathrm{\Δ}{m}_1\\ \·& & \·& & \·\\ \·& \·\·\·& \·& \·\·\·& \·\\ \·& & \·& & \·\\ m_k-4\mathrm{\Δ}{m}_k& \·\·\·& m_k& \·\·\·& m_k+4\mathrm{\Δ}{m}_k\\ \·& & \·& & \·\\ \·& \·\·\·& \·& \·\·\·& \·\\ \·& & \·& & \·\\ m_N-4\mathrm{\Δ}{m}_N& \·\·\·& m_N& \·\·\·& m_N+4\mathrm{\Δ}{m}_N\end{array}\right]$

$\left[\begin{array}{ccccc}{\tilde{v}}_{n,x}-\frac{n_4-1}{2}\mathrm{\Δ}{v}_x& \·\·\·& {\tilde{v}}_{n,x}& \·\·\·& {\tilde{v}}_{n,x}+\frac{n_4-1}{2}\mathrm{\Δ}{v}_x\\ {\tilde{v}}_{n,y}-\frac{n_5-1}{2}\mathrm{\Δ}{v}_y& \·\·\·& {\tilde{v}}_{n,y}& \·\·\·& {\tilde{v}}_{n,y}+\frac{n_5-1}{2}\mathrm{\Δ}{v}_y\\ {\tilde{v}}_{n,z}-\frac{n_6-1}{2}\mathrm{\Δ}{v}_z& \·\·\·& {\tilde{v}}_{n,z}& \·\·\·& {\tilde{v}}_{n,z}+\frac{n_6-1}{2}\mathrm{\Δ}{v}_z\\ {\widetilde{\stackrel{\·}{b}}}_n-\frac{n_8-1}{2}\mathrm{\Δ}\stackrel{\·}{b}& \·\·\·& {\widetilde{\stackrel{\·}{b}}}_n& \·\·\·& {\widetilde{\stackrel{\·}{b}}}_n+\frac{n_8-1}{2}\mathrm{\Δ}\stackrel{\·}{b}\end{array}\right]\⇒\left[\begin{array}{ccccc}{\mathrm{\ω}}_1-4\mathrm{\Δ}{\mathrm{\ω}}_1& \·\·\·& {\mathrm{\ω}}_1& \·\·\·& {\mathrm{\ω}}_1+4\mathrm{\Δ}{\mathrm{\ω}}_1\\ \·& & \·& & \·\\ \·& \·\·\·& \·& \·\·\·& \·\\ \·& & \·& & \·\\ {\mathrm{\ω}}_k-4\mathrm{\Δ}{\mathrm{\ω}}_k& \·\·\·& {\mathrm{\ω}}_k& \·\·\·& {\mathrm{\ω}}_k+4\mathrm{\Δ}{\mathrm{\ω}}_k\\ \·& & \·& & \·\\ \·& \·\·\·& \·& \·\·\·& \·\\ \·& & \·& & \·\\ {\mathrm{\ω}}_N-4\mathrm{\Δ}{\mathrm{\ω}}_N& \·\·\·& {\mathrm{\ω}}_N& \·\·\·& {\mathrm{\ω}}_N+4\mathrm{\Δ}{\mathrm{\ω}}_N\end{array}\right]$

其中，m₁,m₂,...,m_N为第二网格中心网格对应的各颗卫星的码相位；ω₁,ω₂,...,ω_N为第二网格中心网格对应的各颗卫星的载波圆频率。

对于卫星k，以每个网格中心对应码相位、载波多普勒频率设置相关器，则每颗卫星需要设置81个相关器，预检测积分时间T后，得到81个相关值，进入搜索阶段；首先按照搜索中心固定位置、垂直速度、钟偏和钟漂，将卫星k的81个相关值逆映射到水平速度搜索二维平面内；依次将其他卫星的相关值逆映射到水平速度搜索二维平面内，之后按照各颗卫星载噪比权重相加，最大值点即对应水平速度误差。

在本发明的方法的一个实施例中，在水平速度搜索结束后，固定水平速度误差，按照相同的方法依次进行垂直速度、钟漂、水平位置、高度、钟偏误差搜索。

在本发明的方法的一个实施例中，卡尔曼滤波误差模型向量为：

其中，地理坐标系选为北东地坐标系，下标N、E、D分别表示北向，东向，地向，为惯导姿态误差，δV为惯导速度误差，为惯导陀螺常值漂移，为惯导加速度计常值漂移，δL、δλ和δh分别为惯导纬度误差、经度误差和高度误差，b为时钟偏差，为时钟漂移；

卡尔曼滤波状态方程为：

其中：

$F=\left[\begin{array}{cccccc}{F}_{11}& F_{12}& F_{13}& F_{14}& 0_{3\×3}& 0_{3\×2}\\ F_{21}& F_{22}& F_{23}& 0_{3\×3}& F_{25}& 0_{3\×2}\\ 0_{3\×3}& F_{32}& F_{33}& 0_{3\×3}& 0_{3\×3}& 0_{3\×2}\\ 0_{6\×3}& 0_{6\×3}& 0_{6\×3}& 0_{6\×3}& 0_{6\×3}& 0_{3\×2}\\ 0_{2\×3}& 0_{2\×3}& 0_{2\×3}& 0_{2\×3}& 0_{2\×3}& F_{46}\end{array}\right],W=\left[\begin{array}{c}-C_b^n{\mathrm{\ϵ}}_w^b\\ C_b^n{\▿}_w^b\\ 0_{11\×1}\end{array}\right]$

$F_{11}=(-{\mathrm{\ω}}_{\mathrm{in}}^n\×)=\left[\begin{array}{ccc}0& -({\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_E}{R_N+h}\tan L)& \frac{V_N}{R_M+h}\\ {\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_E}{R_N+h}\tan L& 0& {\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E}{R_N+h}\\ -\frac{V_N}{R_M+h}& -({\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E}{R_N+h})& 0\end{array}\right]$

$F_{12}=\left[\begin{array}{ccc}0& \frac{1}{R_N+h}& 0\\ -\frac{1}{R_M+h}& 0& 0\\ 0& -\frac{\tan L}{R_N+h}& 0\end{array}\right]$

$F_{13}=\left[\begin{array}{ccc}-{\mathrm{\ω}}_{\mathrm{ie}}\sin L& 0& 0\\ 0& 0& 0\\ -{\mathrm{\ω}}_{\mathrm{ie}}\cos L-\frac{V_E}{R_N+h}{\sec }^{2}L& 0& 0\end{array}\right]$

$F_{14}=-\left[\begin{array}{ccc}\cos \mathrm{\θ}\cos \mathrm{\ψ}& -\cos \mathrm{\γ}\sin \mathrm{\ψ}+\sin \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}& \sin \mathrm{\γ}\sin \mathrm{\ψ}+\cos \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}\\ \cos \mathrm{\θ}\sin \mathrm{\ψ}& \cos \mathrm{\γ}\cos \mathrm{\ψ}+\sin \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}& -\sin \mathrm{\γ}\cos \mathrm{\ψ}+\cos \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}\\ -\sin \mathrm{\θ}& \sin \mathrm{\γ}\cos \mathrm{\θ}& \cos \mathrm{\γ}\cos \mathrm{\θ}\end{array}\right]$

$F_{21}=(f^n\×)=\left[\begin{array}{ccc}0& -f_z^n& f_y^n\\ f_z^n& 0& -f_x^n\\ -f_y^n& f_x^n& 0\end{array}\right]$

$F_{22}=\left[\begin{array}{ccc}\frac{V_D}{R_M+h}& -2({\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_e}{R_N+h}\tan L)& \frac{V_N}{R_M+h}\\ {2\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_E}{R_N+h}\tan L& \frac{V_D}{R_N+h}+\frac{V_N}{R_N+h}\tan L& {2\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E}{R_N+h}\\ -\frac{{2V}_N}{R_M+h}& -2({\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E}{R_N+h})& 0\end{array}\right]$

$F_{23}=\left[\begin{array}{ccc}-({2V}_E{\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E^2}{R_N+h}{\sec }^{2}L)& 0& 0\\ {2V}_N{\mathrm{\ω}}_{\mathrm{ie}}\cos L-{2V}_E{\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_N{V}_E}{R_N+h}{\sec }^{2}L& 0& 0\\ {2V}_E{\mathrm{\ω}}_{\mathrm{ie}}\sin L& 0& -\frac{2g}{R_M}\end{array}\right]$

$F_{25}=\left[\begin{array}{ccc}\cos \mathrm{\θ}\cos \mathrm{\ψ}& -\cos \mathrm{\γ}\sin \mathrm{\ψ}+\sin \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}& \sin \mathrm{\γ}\sin \mathrm{\ψ}+\cos \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}\\ \cos \mathrm{\θ}\sin \mathrm{\ψ}& \cos \mathrm{\γ}\cos \mathrm{\ψ}+\sin \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}& -\sin \mathrm{\γ}\cos \mathrm{\ψ}+\cos \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}\\ -\sin \mathrm{\θ}& \sin \mathrm{\γ}\cos \mathrm{\θ}& \cos \mathrm{\γ}\cos \mathrm{\θ}\end{array}\right]$

$F_{32}=\left[\begin{array}{ccc}\frac{1}{R_M+h}& 0& 0\\ 0& \frac{1}{(R_N+h)\cos }& 0\\ 0& 0& -1\end{array}\right]F_{23}=\left[\begin{array}{ccc}0& 0& -\frac{V_E}{{(R_M+h)}^2}\\ \frac{V_E\tan L}{(R_N+h)\cos L}& 0& -\frac{V_E}{{(R_N+h)}^2\cos L}\\ 0& 0& 0\end{array}\right]F_{46}=\left[\begin{array}{cc}0& 1\\ 0& 0\end{array}\right]$

其中，为载体坐标系到导航坐标系的转换矩阵，θ为俯仰角，γ为横滚角，ψ为航向角，V_N、V_E和V_D分别为载体的北向、东向和地向速度，L为载***置的地理纬度，h为载体海拔高度，R_M为地球子午圈曲率半径，R_N为地球卯酉圈曲率半径，ω_ie为地球的自转角速度，和为惯导的比力，为陀螺的噪声，为加速度计的噪声；

卡尔曼滤波器观测向量为： $Z=[\mathrm{\Δx},\mathrm{\Δy},\mathrm{\Δz},{\mathrm{\Δv}}_N,{\mathrm{\Δv}}_E,{\mathrm{\Δv}}_D,\mathrm{\Δb},\mathrm{\Δ}\stackrel{\·}{b}];$

卡尔曼滤波器观测方程为：Z＝HX+V；

其中，

$H=\left[\begin{array}{cccc}{0}_{3\×3}& 0_{3\×3}& I_{3\×3}& 0_{3\×2}\\ 0_{3\×3}& I_{3\×3}& 0_{3\×3}& 0_{3\×2}\\ 0_{2\×3}& 0_{2\×3}& 0_{2\×3}& I_{2\×2}\end{array}\right]$ 为观测矩阵，V为观测噪声。

在本发明的方法的一个实施例中，在步骤S06中，捷联惯导解算速率与导航电文码率相等或呈倍数关系。

在本发明的方法的一个实施例中，其特征在于，在步骤07中，根据卡尔曼滤波器状态协方差矩阵自适应调整搜索网格大小。

在本发明的方法的一个实施例中，所述GNSS***为GPS、北斗卫星定位***、格洛纳斯***或伽利略导航***。

本发明的第二目的在于提出一种具有融合程度好、导航精度高的GNSS/MINS超深组合导航***。

根据本发明实施例的GNSS/MINS超深组合导航***，包括：GNSS***，通过天线(601)接收卫星信号，经射频前端(602)，将所述卫星信号下变频采样后输出中频信号，并接收卫星导航电文；MINS***，通过MIMU(603)获取测量数据，经捷联惯导***(604)解算后得到载***置、速度及姿态；载波频率及码相位计算单元(605)，用于将所述MINS***输出的速度和位置结合所述卫星导航电文计算GNSS信号载波频率和码相位，以产生本地信号；以及组合卡尔曼滤波器(609)，将状态误差反馈回MINS***进行误差修正；其特征在于，还包括：相关器组(606)，用于将所述本地信号与中频信号进行相关；投影映射单元(607)，将所述相关器组(606)的输出结果完成从第二网格到第一网格的投影映射，其中，所述第一网格为在三维位置、三维速度、钟漂、钟偏共八维变量上的网格划分，第二网格为在GNSS载波频率、伪码相位等二维变量上的网格划分；以及搜索单元(608)，对所述投影映射单元(607)的输出结果进行搜索，获得最大值，即位置误差、速度误差、钟偏误差和钟漂误差，作为所述组合卡尔曼滤波器的观测量。

在本发明的***的一个实施例中，所述搜索单元(608)按照水平速度、垂直速度、钟漂、水平位置、高度、钟偏的顺序进行最大值的搜索，获得对应T时间内的水平速度误差、垂直速度误差、钟漂误差、水平位置误差、高度误差、钟偏误差。

在本发明的***的一个实施例中，所述GNSS***为GPS、北斗卫星定位***、格洛纳斯***或伽利略导航***。

在本发明的***的一个实施例中，所述***还包括气压高度表，为***引入气压高度数据。

在本发明的***的一个实施例中，所述***还包括高精度晶振。

本发明的第三目的在于提出一种具有融合程度好、导航精度高的GNSS/MINS超深组合导航***。

根据本发明实施例的组合导航装置，包括：天线和射频模块，通过FPGA接口将GNSS中频数据引入***；MIMU，通过FPGA接口将MIMU的测量数据引入***；以及电源模块、DSP和上位机；其特征在于，所述组合导航装置还包括相关器组，所述DSP将计算码相位、载波频率和投影矢量，并通过EMIF总线写入FPGA内部，FPGA根据搜索网格大小自动划分网格并映射到所述相关器组，预检测积分时间T后，FPGA处理相关器组的输出值并找到最大值，即对应T时间内的位置误差、速度误差、钟偏和钟漂误差，所述DSP读入所述位置误差、速度误差、钟偏和钟漂误差观测值并对用户状态量进行更新。

在本发明的装置的一个实施例中，还包括气压高度表，通过FPGA接口将气压高度数据引入***。

在本发明的装置的一个实施例中，所述GNSS***为GPS、北斗卫星定位***、格洛纳斯***或伽利略导航***。

现有技术中多采用矢量算法，由于利用多颗卫星信号联合定位，其导航性能优于标量算法，但是当卫星信号载噪比较低时，矢量算法鉴相器误差急剧增加，导航性能下降。本发明主要采用直接位置估计算法(DPE)而不使用鉴相器，使用贝叶斯估计算法，理论上能达到最高的精度，实现更深层次GNSS/INS组合。本发明的GNSS/MINS超深组合导航方法、***及装置具有融合程度好，导航精度高的优点。

本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1为直接位置估计算法原理图。

图2为第一网格到第二网格映射示意图。

图3为实施例一的超深组合导航方法流程图。

图4为本实施例一种所使用的水平速度误差搜索方法的工作原理图。

图5为水平速度搜索时网格划分和映射的原理图。

图6为实施例二的超深组合导航***的结构图。

图7为实施例三的超深组合导航装置的结构图。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。

在本发明的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“长度”、“宽度”、“厚度”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”“内”、“外”、“顺时针”、“逆时针”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明的描述中，“多个”的含义是两个或两个以上，除非另有明确具体的限定。

在本发明中，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”、“固定”等术语应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。

在本发明中，除非另有明确的规定和限定，第一特征在第二特征之“上”或之“下”可以包括第一和第二特征直接接触，也可以包括第一和第二特征不是直接接触而是通过它们之间的另外的特征接触。而且，第一特征在第二特征“之上”、“上方”和“上面”包括第一特征在第二特征正上方和斜上方，或仅仅表示第一特征水平高度高于第二特征。第一特征在第二特征“之下”、“下方”和“下面”包括第一特征在第二特征正下方和斜下方，或仅仅表示第一特征水平高度小于第二特征。

实施例一

任何卫星定位***，其本质和目标都是确定用户的位置和速度。本发明正是基于这一共同点，提出基于位置、速度层面的定位场概念。基于这一概念，可以将INS与GNSS进行超深层次组合：利用INS提供的位置、速度辅助信息，GNSS使用直接位置估计(DPE)算法估计位置误差、速度误差、时钟误差，以此闭合环路。由于***直接以MIMU、基带相关器为传感器，在顶层组合算法中实现所有的导航功能，即该算法将GNSS基带相关器作为敏感空时定位场的传感器，将MIMU作为敏感惯性场的传感器，实现了GNSS与MINS的完整融合，在INS辅助下，GNSS联合多通道信息进行矢量鉴相和矢量定位，因此称为超深组合导航算法。

为了敏感卫星定位场，本设计采用相关器组作为传感元件，相关器输出值主要受用户位置、速度、钟漂、钟偏影响，因此，考虑在用户三维位置、三维速度、钟漂、钟偏8个维度上进行网格划分，针对每颗卫星，每一种网格组合对应一个相关器，预检测积分时间T后相关器组输出一组累加值，其中最大的相关器累加值对应正确的三维位置、三维速度、钟漂和钟偏。利用观测值修正用户位置、速度、钟偏和钟漂后，再开始下一轮搜索。

具体算法原理如图1所示。

T₀时刻用户位置、速度、钟偏和钟漂已知，以该初始状态为中心划分网格位置搜索格101、速度搜索格102、钟偏搜索网格、钟漂搜索网格，可以映射到各颗卫星载波多普勒和码相位搜索格103。假设有n颗可视卫星，将每个网格对应n颗卫星的相关器输出值进行累加，所有累加值中最大值对应正确的用户状态误差。误差更新后可进行下一轮搜索并更新。

假设在三维位置、三维速度、钟偏、钟漂8个维度划分搜索网格的数量分别为n₁～n₈，则划分的网格总数量为：

N＝n₁×n₂×n₃×n₄×n₅×n₆×n₇×n₈

对于典型的GNSS导航应用，考虑每个维度划分5个网格，即n₁＝n₂＝n₃＝n₄＝n₅＝n₆＝n₇＝n₈＝5，则此时N＝390625，即需要390625个相关器才能完成1颗卫星信号的跟踪，数量如此巨大的相关器组是无法实现的，限制了基于卫星定位场的直接位置估计算法在GNSS导航领域中的应用。

可以使用一些策略来减少相关器数量，优化算法复杂度。

使用高精度晶振，钟偏和钟漂2个维度搜索格子数量可以减少，典型情况每个维度选用3个搜索格。

高度通道对相关器累加值影响较小，因此可以引入气压高度表辅助导航，从而减少高度方向的搜索格数量，典型情况高度及其导数每个维度选用3个搜索格。

在INS辅助下，可以压缩用户水平方向位置、速度搜索格大小，典型情况下，位置搜索格大小为5m，范围30m，经过插值可以得到1m的精度；速度搜索格大小为2m/s，范围12m/s，经过插值可以得到小于0.5m/s的精度。即水平位置、速度4个维度各需要7个搜索格。

综上，网格总数为：N＝n₁×……×n₈＝3⁴×9⁴＝194481。可见数量还是过大，无法硬件实现。

对于GNSS信号而言，只要给出准确的载波频率，伪码相位，则可以保证稳定跟踪信号。因此，对于每一颗卫星而言，可以将用户三维位置、三维速度、钟偏、钟漂8维变量映射到载波频率、码相位2维变量上，从而减少相关器的数量。即三维位置、三维速度、钟偏、钟漂到载波频率、码相位的映射为多对一线性映射。

如图2所示，位于x₁y₁z₁区域内的网格，都将映射到投影矢量上的ρ₁区域，位于x₂y₂z₂区域内的网格，都将映射到投影矢量上的ρ₂区域，从而实现3维到1维的映射。

因此，在用户到卫星投影矢量已知的情况下，可以将三维速度、钟漂、三维位置、钟偏8维变量投影到载波频率、码相位两个维度，然后在载波频率、码相位两个维度划分网格，得到相关器输出后，再通过投影矢量逆运算映射到三维速度、钟漂、三维位置、钟偏网格中，从而将8维网格划分压缩为2维网格划分，大幅度减少相关器数量，保证硬件可实现性。

本实施例方法如图3所示，具体按照以下步骤：

步骤S01，组合导航***初始化参数设置；

步骤S02，GNSS捕获卫星并获取导航电文，INS在GNSS辅助下完成大失准角初始对准；

步骤S03，根据INS输出的位置和速度以及GNSS输出的初始钟偏和初始钟漂，进行第一网格划分，并完成从位置、速度、钟偏和钟漂到载相频率和码相位的映射，实现第二网格划分，对应每个第二网格中心点建立相关器；

步骤S04，进行1比特电文长度的相干积分，相关器输出相关值；

步骤S05，判断是否达到积分时间T，若未达到，进入步骤S06，反之进入步骤S07；

步骤S06，***进行非相干累加以及捷联导航解算，并根据所述捷联导航解算结果实时调整所述第二网格划分中心位置，捷联惯导解算速率与导航电文码率相等或呈倍数关系。返回步骤S05；

步骤S07，将每颗卫星的相关值进行第二网格到第一网格的逆映射，然后进行最大值的搜索，获得T时间内的位置误差、速度误差、钟偏误差和钟漂误差；；

步骤S08，将步骤S07中得到的各误差作为组合卡尔曼滤波器的观测量，进行INS误差更新及时钟误差更新，然后将卡尔曼滤波器状态清零；

步骤S09，搜索网格大小自适应调整；

步骤S10，捷联导航解算更新用户位置、速度以后，进入下一算法循环。

由于搜索网格数量有限，网格划分越大，卫星信号跟踪能力越强，但是跟踪精度越差。因此要根据卡尔曼滤波器状态协方差矩阵自适应调整搜索网格大小，保证稳健跟踪卫星信号的同时，尽可能提高跟踪精度。

在步骤S07进行搜索时，按照对相关器输出值影响大小排序，依次为水平速度、垂直速度、钟漂、水平位置、高度、钟偏。因此，可以按照初始状态固定三维位置、垂直速度、钟偏和钟漂，对水平速度进行搜索得到准确的水平速度误差，然后固定误差更新后的水平速度，对垂直速度进行搜索，依次类推。在计算资源允许的情况下，可以按照上面顺序进行迭代，直到算法收敛。该搜索策略将8维搜索划分为多个2维以下的搜索，有效降低算法复杂度。

矢量环路受限于鉴相器非线性以及鉴相器工作范围，预检测积分时间T一般为20～200ms，对弱信号跟踪能力有限，并且当载噪比降低时，鉴相器误差增大，容易使环路失锁。相较于矢量环路，直接位置估计(DPE)算法环路中没有使用鉴相器，因此不受鉴相器非线性及工作范围的影响，预检测积分时间T可选为1～3秒，并且在低载噪比情况下也有良好的性能。

在本实施例中，***模型的建立包括DPE算法和组合算法。

DPE算法使用最大似然估计方法，DPE方法等效于最优化：

其中，是待估计的接收机位置、钟偏、速度和钟漂，r(t)是GNSS接收基带信号加噪声，s_k(t)是接收机产生的第k颗卫星基带信号，N是可视卫星数量，T₀是预检测积分起始时刻，T是预检测积分时间，t是接收机生成基带信号时间序列。假设s_k(t)跨越D bit导航电文，典型情况下D取值50～150(1～3秒)。

在第n数据位，第k颗卫星本地基带信号s_k(t)为：

s_k(t)＝a_k,nC_k(t)exp[j(ω_kt+φ_k,n)]   (2)

其中a_k,n为信号幅度，与载噪比有关；C_k(t)为CA码，由用户位置、钟偏决定；ω_k为载波圆频率，有用户速度、钟漂决定；φ_k,n为载波相位。由于信号幅度可以通过载噪比估计值进行加权，载波相位估计要求网格划分过密而无法实现，因此下面的讨论中不考虑a_k,n和φ_k,n的影响。

将(1)式展开，得到：

由于R只与接收信号有关，因此最小化[R-S]，等效于最大化[S]，即：

其中I_k,n是第k颗卫星第n数据位的时间跨度，对GPS而言大小为20ms。(3)式中内部求和号为非相干累加，外部求和号为各卫星通道融合累加。

CA码序列波形C_k(t)求解式为：

$C_k\left(t\right)\≅C_k[t-b\left(T_0\right)-{\mathrm{\τ}}_k\left(T_0\right)-(t-T_0){\mathrm{\τ}}_k^{\′}\left(T_0\right)]---\left(4\right)$

其中C_k(t_GNSS)为以GNSS时间为参数的卫星CA码产生函数，C_k(t_GNSS)已知；τ_k(T₀)为T₀到达的第k颗卫星信号传输延时；τ'_k(T₀)为由多普勒频率引起的延迟时间变化率；b(T₀)为钟偏，由导航电文给出。几个参数计算如下：

${\mathrm{\τ}}_k\left(T_0\right)\≅\frac{|p_k\left(T_0\right)-p\left(T_0\right)|}{c+v_k\left(T_0\right)\·u_k\left(T_0\right)}\sec ---\left(5\right)$

${\mathrm{\τ}}_k^{\′}\left(T_0\right)\≅\frac{1}{c}\left\{\right[v_k\left(T_0\right)-v\left(T_0\right)]\·u_k\left(T_0\right)\}\sec /\sec ---\left(6\right)$

其中p_k(T₀)为第k颗卫星位置；p(T₀)为用户位置；v_k(T₀)为第k颗卫星速度；u_k(T₀)为用户到第k颗卫星的投影矢量；v(T₀)为用户速度；c为光速。可见，已知导航电文、用户位置、钟偏的情况下，根据(4)～(6)式可以解算CA码序列m_k(t)。

基带多普勒频率ω_k计算式为：

${\mathrm{\ω}}_k=2\mathrm{\π}\left\{\frac{1}{\mathrm{\λ}}\right[v\left(T_0\right)-v_k\left(T_0\right)]\·u_k\left(T_0\right)-\stackrel{\·}{b}\left(T_0\right)f\}---\left(7\right)$

其中f为载波频率，λ为载波波长。可见，已知导航电文、用户速度、钟漂的情况下，可以求解基带多普勒频率ω_k。

按照水平速度、垂直速度、钟漂、水平位置、高度、钟偏的顺序进行搜索，由于搜索方法类似，下面以水平速度搜索为例进行说明，工作原理如图4所示：

假设第n数据位起始时刻INS输出位置、速度为钟偏和钟漂为以这组参数为中心进行8维网格划分，每个维度划分网格的数量依次为n₁～n₈，且n₁～n₈均为奇数。

由于用户位置、钟偏到码相位，用户速度、钟偏到载波多普勒频率的映射为线性映射，因此可以根据n₁～n₈网格划分求出各颗卫星最大的码相位偏移和载波多普勒频率偏移，无论偏移大小，均在码相位、载波多普勒两个维度均分为9份，如图5所示：

$\left[\begin{array}{ccccc}{\tilde{p}}_{n,x}-\frac{n_1-1}{2}\mathrm{\Δ}{p}_x& \·\·\·& {\tilde{p}}_{n,x}& \·\·\·& {\tilde{p}}_{n,x}+\frac{n_1-1}{2}\mathrm{\Δ}{p}_x\\ {\tilde{p}}_{n,y}-\frac{n_2-1}{2}\mathrm{\Δ}{p}_y& \·\·\·& {\tilde{p}}_{n,y}& \·\·\·& {\tilde{p}}_{n,y}+\frac{n_2-1}{2}\mathrm{\Δ}{p}_y\\ {\tilde{p}}_{n,z}-\frac{n_3-1}{2}\mathrm{\Δ}{p}_z& \·\·\·& {\tilde{p}}_{n,z}& \·\·\·& {\tilde{p}}_{n,z}+\frac{n_3-1}{2}\mathrm{\Δ}{p}_z\\ {\tilde{b}}_n-\frac{n_7-1}{2}\mathrm{\Δb}& \·\·\·& {\tilde{b}}_n& \·\·\·& {\tilde{b}}_n+\frac{n_7-1}{2}\mathrm{\Δb}\end{array}\right]\⇒\left[\begin{array}{ccccc}{m}_1-4\mathrm{\Δ}{m}_1& \·\·\·& m_1& \·\·\·& m_1+4\mathrm{\Δ}{m}_1\\ \·& & \·& & \·\\ \·& \·\·\·& \·& \·\·\·& \·\\ \·& & \·& & \·\\ m_k-4\mathrm{\Δ}{m}_k& \·\·\·& m_k& \·\·\·& m_k+4\mathrm{\Δ}{m}_k\\ \·& & \·& & \·\\ \·& \·\·\·& \·& \·\·\·& \·\\ \·& & \·& & \·\\ m_N-4\mathrm{\Δ}{m}_N& \·\·\·& m_N& \·\·\·& m_N+4\mathrm{\Δ}{m}_N\end{array}\right]$

$\left[\begin{array}{ccccc}{\tilde{v}}_{n,x}-\frac{n_4-1}{2}\mathrm{\Δ}{v}_x& \·\·\·& {\tilde{v}}_{n,x}& \·\·\·& {\tilde{v}}_{n,x}+\frac{n_4-1}{2}\mathrm{\Δ}{v}_x\\ {\tilde{v}}_{n,y}-\frac{n_5-1}{2}\mathrm{\Δ}{v}_y& \·\·\·& {\tilde{v}}_{n,y}& \·\·\·& {\tilde{v}}_{n,y}+\frac{n_5-1}{2}\mathrm{\Δ}{v}_y\\ {\tilde{v}}_{n,z}-\frac{n_6-1}{2}\mathrm{\Δ}{v}_z& \·\·\·& {\tilde{v}}_{n,z}& \·\·\·& {\tilde{v}}_{n,z}+\frac{n_6-1}{2}\mathrm{\Δ}{v}_z\\ {\widetilde{\stackrel{\·}{b}}}_n-\frac{n_8-1}{2}\mathrm{\Δ}\stackrel{\·}{b}& \·\·\·& {\widetilde{\stackrel{\·}{b}}}_n& \·\·\·& {\widetilde{\stackrel{\·}{b}}}_n+\frac{n_8-1}{2}\mathrm{\Δ}\stackrel{\·}{b}\end{array}\right]\⇒\left[\begin{array}{ccccc}{\mathrm{\ω}}_1-4\mathrm{\Δ}{\mathrm{\ω}}_1& \·\·\·& {\mathrm{\ω}}_1& \·\·\·& {\mathrm{\ω}}_1+4\mathrm{\Δ}{\mathrm{\ω}}_1\\ \·& & \·& & \·\\ \·& \·\·\·& \·& \·\·\·& \·\\ \·& & \·& & \·\\ {\mathrm{\ω}}_k-4\mathrm{\Δ}{\mathrm{\ω}}_k& \·\·\·& {\mathrm{\ω}}_k& \·\·\·& {\mathrm{\ω}}_k+4\mathrm{\Δ}{\mathrm{\ω}}_k\\ \·& & \·& & \·\\ \·& \·\·\·& \·& \·\·\·& \·\\ \·& & \·& & \·\\ {\mathrm{\ω}}_N-4\mathrm{\Δ}{\mathrm{\ω}}_N& \·\·\·& {\mathrm{\ω}}_N& \·\·\·& {\mathrm{\ω}}_N+4\mathrm{\Δ}{\mathrm{\ω}}_N\end{array}\right]$

其中，m₁,m₂,...,m_N为第二网格中心网格对应的各颗卫星的码相位；ω₁,ω₂,...,ω_N为第二网格中心网格对应的各颗卫星的载波圆频率。

对于卫星k，以中心黑点对应码相位、载波多普勒频率设置相关器，则每颗卫星需要设置81个相关器，预检测积分时间T后，得到81个相关值，进入搜索阶段。首先按照搜索中心固定位置、垂直速度、钟偏和钟漂，将卫星k的81个相关值逆映射到水平速度搜索二维平面内。

依次将其他卫星的相关值逆映射到水平速度搜索二维平面内，之后按照各颗卫星载噪比权重相加，最大值点即对应水平速度误差，可作为组合卡尔曼滤波器的观测量。

水平速度误差搜索结束后，固定水平速度误差，搜索垂直速度，依次类推。所有搜索结束，得到速度误差、位置误差、钟偏误差、钟漂误差8个观测量，反馈回组合卡尔曼滤波器对***状态进行更新。

本实施例中，地理坐标系选为北东地坐标系(下标N、E、D分别表示北向，东向，地向)，误差模型状态向量为：

其中，为惯导姿态误差，δV为惯导速度误差，为惯导陀螺常值漂移，为惯导加速度计常值漂移，δL、δλ和δh分别为惯导纬度误差、经度误差和高度误差，b为时钟偏差，为时钟漂移。

卡尔曼滤波状态方程为：

其中：

$F=\left[\begin{array}{cccccc}{F}_{11}& F_{12}& F_{13}& F_{14}& 0_{3\×3}& 0_{3\×2}\\ F_{21}& F_{22}& F_{23}& 0_{3\×3}& F_{25}& 0_{3\×2}\\ 0_{3\×3}& F_{32}& F_{33}& 0_{3\×3}& 0_{3\×3}& 0_{3\×2}\\ 0_{6\×3}& 0_{6\×3}& 0_{6\×3}& 0_{6\×3}& 0_{6\×3}& 0_{3\×2}\\ 0_{2\×3}& 0_{2\×3}& 0_{2\×3}& 0_{2\×3}& 0_{2\×3}& F_{46}\end{array}\right],W=\left[\begin{array}{c}-C_b^n{\mathrm{\ϵ}}_w^b\\ C_b^n{\▿}_w^b\\ 0_{11\×1}\end{array}\right]$

$F_{11}=(-{\mathrm{\ω}}_{\mathrm{in}}^n\×)=\left[\begin{array}{ccc}0& -({\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_E}{R_N+h}\tan L)& \frac{V_N}{R_M+h}\\ {\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_E}{R_N+h}\tan L& 0& {\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E}{R_N+h}\\ -\frac{V_N}{R_M+h}& -({\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E}{R_N+h})& 0\end{array}\right]$

$F_{12}=\left[\begin{array}{ccc}0& \frac{1}{R_N+h}& 0\\ -\frac{1}{R_M+h}& 0& 0\\ 0& -\frac{\tan L}{R_N+h}& 0\end{array}\right]$

$F_{13}=\left[\begin{array}{ccc}-{\mathrm{\ω}}_{\mathrm{ie}}\sin L& 0& 0\\ 0& 0& 0\\ -{\mathrm{\ω}}_{\mathrm{ie}}\cos L-\frac{V_E}{R_N+h}{\sec }^{2}L& 0& 0\end{array}\right]$

$F_{14}=-\left[\begin{array}{ccc}\cos \mathrm{\θ}\cos \mathrm{\ψ}& -\cos \mathrm{\γ}\sin \mathrm{\ψ}+\sin \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}& \sin \mathrm{\γ}\sin \mathrm{\ψ}+\cos \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}\\ \cos \mathrm{\θ}\sin \mathrm{\ψ}& \cos \mathrm{\γ}\cos \mathrm{\ψ}+\sin \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}& -\sin \mathrm{\γ}\cos \mathrm{\ψ}+\cos \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}\\ -\sin \mathrm{\θ}& \sin \mathrm{\γ}\cos \mathrm{\θ}& \cos \mathrm{\γ}\cos \mathrm{\θ}\end{array}\right]$

$F_{21}=(f^n\×)=\left[\begin{array}{ccc}0& -f_z^n& f_y^n\\ f_z^n& 0& -f_x^n\\ -f_y^n& f_x^n& 0\end{array}\right]$

$F_{22}=\left[\begin{array}{ccc}\frac{V_D}{R_M+h}& -2({\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_e}{R_N+h}\tan L)& \frac{V_N}{R_M+h}\\ {2\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_E}{R_N+h}\tan L& \frac{V_D}{R_N+h}+\frac{V_N}{R_N+h}\tan L& {2\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E}{R_N+h}\\ -\frac{{2V}_N}{R_M+h}& -2({\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E}{R_N+h})& 0\end{array}\right]$

$F_{23}=\left[\begin{array}{ccc}-({2V}_E{\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E^2}{R_N+h}{\sec }^{2}L)& 0& 0\\ {2V}_N{\mathrm{\ω}}_{\mathrm{ie}}\cos L-{2V}_E{\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_N{V}_E}{R_N+h}{\sec }^{2}L& 0& 0\\ {2V}_E{\mathrm{\ω}}_{\mathrm{ie}}\sin L& 0& -\frac{2g}{R_M}\end{array}\right]$

$F_{25}=\left[\begin{array}{ccc}\cos \mathrm{\θ}\cos \mathrm{\ψ}& -\cos \mathrm{\γ}\sin \mathrm{\ψ}+\sin \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}& \sin \mathrm{\γ}\sin \mathrm{\ψ}+\cos \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}\\ \cos \mathrm{\θ}\sin \mathrm{\ψ}& \cos \mathrm{\γ}\cos \mathrm{\ψ}+\sin \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}& -\sin \mathrm{\γ}\cos \mathrm{\ψ}+\cos \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}\\ -\sin \mathrm{\θ}& \sin \mathrm{\γ}\cos \mathrm{\θ}& \cos \mathrm{\γ}\cos \mathrm{\θ}\end{array}\right]$

$F_{32}=\left[\begin{array}{ccc}\frac{1}{R_M+h}& 0& 0\\ 0& \frac{1}{(R_N+h)\cos }& 0\\ 0& 0& -1\end{array}\right]F_{23}=\left[\begin{array}{ccc}0& 0& -\frac{V_E}{{(R_M+h)}^2}\\ \frac{V_E\tan L}{(R_N+h)\cos L}& 0& -\frac{V_E}{{(R_N+h)}^2\cos L}\\ 0& 0& 0\end{array}\right]F_{46}=\left[\begin{array}{cc}0& 1\\ 0& 0\end{array}\right]$

其中，为载体坐标系到导航坐标系的转换矩阵，θ为俯仰角，γ为横滚角，ψ为航向角，V_N、V_E和V_D分别为载体的北向、东向和地向速度，L为载***置的地理纬度，h为载体海拔高度，R_M为地球子午圈曲率半径，R_N为地球卯酉圈曲率半径，ω_ie为地球的自转角速度，和为惯导的比力，为陀螺的噪声，为加速度计的噪声；

卡尔曼滤波器观测向量为：由DPE最大似然估计得到，对于速度误差观测量，需要由ECEF坐标系转换到导航坐标系。

卡尔曼滤波器观测方程为：Z＝HX+V

其中， $H=\left[\begin{array}{cccc}{0}_{3\×3}& 0_{3\×3}& I_{3\×3}& 0_{3\×2}\\ 0_{3\×3}& I_{3\×3}& 0_{3\×3}& 0_{3\×2}\\ 0_{2\×3}& 0_{2\×3}& 0_{2\×3}& I_{2\×2}\end{array}\right]$ 为观测矩阵，V为观测噪声。

利用INS提供的位置、速度辅助信息，GNSS使用直接位置估计(DPE)算法估计位置误差、速度误差、时钟误差，以此闭合环路。由于***直接以MIMU、基带相关器为传感器，在顶层组合算法中实现所有的导航功能，即该算法将GNSS基带相关器作为敏感空间定位场的传感器，实现了GNSS与MINS的完整融合。

本实施例可以兼容所有卫星定位***，包括GPS、北斗卫星定位***、格洛纳斯***和伽利略导航***。

实施例二

如图6所示为本实施例的***架构图，该GNSS/MINS超深组合导航***主要由GNSS***、MINS***和相关器组606构成。

其中，GNSS***通过天线601接收卫星信号，经射频前端602，将所述卫星信号下变频采样后输出中频信号，并接收卫星导航电文；MINS***通过MIMU603获取测量数据，经捷联惯导***604解算后得到载***置、速度及姿态。

该***还包括：

载波频率及码相位计算单元605，用于将所述MINS***输出的速度和位置结合所述卫星导航电文计算GNSS信号载波频率和码相位，以产生本地信号；

组合卡尔曼滤波器609，将状态误差反馈回MINS***进行误差修正；

相关器组606，用于将所述本地信号与中频信号进行相关；

投影映射单元607，将所述相关器组606的输出结果进行投影映射；

搜索单元608，对所述投影映射单元607的输出结果进行搜索，获得最大值，即位置误差、速度误差、钟偏误差和钟漂误差，作为所述组合卡尔曼滤波器的观测量。

所述搜索单元608按照水平速度、垂直速度、钟漂、水平位置、高度、钟偏的顺序进行最大值的搜索，获得对应T时间内的水平速度误差、垂直速度误差、钟漂误差、水平位置误差、高度误差、钟偏误差。

本实施例可以兼容所有卫星定位***，包括GPS、北斗卫星定位***、格洛纳斯***和伽利略导航***。

本实施例使用高精度晶振，使钟偏和钟漂2个维度搜索格子数量可以减少。高度通道对相关器累加值影响较小，因此可以引入气压高度表辅助导航，从而减少高度方向的搜索格数量。

实施例三

本实施例提出一种应用于实施例一的方法的装置，如图7所示。

该GNSS/MINS组合导航装置，包括：天线701和射频模块702，通过FPGA706的接口将GNSS中频数据引入***；气压高度表703，通过FPGA706的接口将气压高度数据引入***；MIMU704，通过FPGA706的接口将MIMU704)的测量数据引入***；以及电源模块705、DSP708和上位机707。

其中，所述组合导航装置还包括相关器组709，所述DSP708将计算码相位、载波频率和投影矢量，并通过EMIF总线写入FPGA706内部，FPGA706根据搜索网格大小自动划分网格并映射到所述相关器组709，预检测积分时间T后，FPGA706处理相关器组709的输出值并找到最大值，即对应T时间内的位置误差、速度误差、钟偏和钟漂误差，所述DSP708读入所述位置误差、速度误差、钟偏和钟漂误差观测值并对用户状态量进行更新。

本实施例可以兼容所有卫星定位***，包括GPS、北斗卫星定位***、格洛纳斯***和伽利略导航***。

本实施例使用高精度晶振，使钟偏和钟漂2个维度搜索格子数量可以减少。高度通道对相关器累加值影响较小，因此可以引入气压高度表辅助导航，从而减少高度方向的搜索格数量。

该硬件架构充分利用FPGA并行处理能力，使用FPGA完成大量的相关运算和搜索运算，搭配一定的搜索策略，可以保证超深组合导航***的嵌入式实现。

需要说明的是，流程图中或在此以其他方式描述的任何过程或方法描述可以被理解为，表示包括一个或更多个用于实现特定逻辑功能或过程的步骤的可执行指令的代码的模块、片段或部分，并且本发明的优选实施方式的范围包括另外的实现，其中可以不按所示出或讨论的顺序，包括根据所涉及的功能按基本同时的方式或按相反的顺序，来执行功能，这应被本发明的实施例所属技术领域的技术人员所理解。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。

Claims (20)

1.一种GNSS/MINS超深组合导航方法，以MIMU和GNSS基带相关器为传感器，进行组合导航计算，包括：

步骤S01，组合导航***初始化参数设置；

步骤S02，GNSS捕获卫星并获取导航电文，INS在GNSS辅助下完成大失准角初始对准；

步骤S03，根据INS输出的位置和速度以及GNSS输出的初始钟偏和初始钟漂，进行第一网格划分，并完成从位置、速度、钟偏和钟漂到载相频率和码相位的映射，实现第二网格划分，对应每个第二网格中心点建立相关器；

步骤S04，进行1比特电文长度的相干积分，相关器输出相关值；

步骤S05，判断是否达到积分时间T，若未达到，进入步骤S06，反之进入步骤S07；

步骤S06，***进行非相干累加以及捷联导航解算，并根据所述捷联导航解算结果实时调整所述第二网格划分中心位置，返回步骤S05；

步骤S07，将每颗卫星的非相干累加相关值进行第二网格到第一网格的逆映射，然后进行最大值的搜索，获得T时间内的位置误差、速度误差、钟偏误差和钟漂误差；

步骤S08，将步骤S07中得到的各误差作为组合卡尔曼滤波器的观测量，进行INS误差更新及时钟误差更新，然后将卡尔曼滤波器状态清零；

步骤S09，搜索网格大小自适应调整；

步骤S10，捷联导航解算更新用户位置、速度以后，进入下一算法循环。

2.根据权利要求1所述的组合导航方法，其特征在于，所述第一网格划分为在三维位置、三维速度、钟漂、钟偏上的八维网格划分，第二网格划分为在GNSS载波频率、伪码相位上的二维网格划分。

3.根据权利要求2所述的组合导航方法，其特征在于，所述超深组合导航方法利用气压高度表引入气压高度数据辅助导航。

4.根据权利要求2所述的组合导航方法，其特征在于，所述超深组合导航方法应用高精度晶振。

5.根据权利要求3或4所述的组合导航方法，其特征在于，所述超深组合导航方法使用最大似然估计方法，最大似然估计式为：

其中，是待估计的接收机位置、钟偏、速度和钟漂，r(t)是GNSS接收基带信号加噪声，s_k(t)是接收机产生的第k颗卫星本地基带信号，N是可视卫星数量，T₀是预检测积分起始时刻，T是预检测积分时间，t是接收机生成基带信号时间序列，假设s_k(t)跨越D比特导航电文，典型情况下D取值50～150；

在第n数据位，第k颗卫星本地基带信号s_k(t)为：

s_k(t)＝a_k，nC_k(t)exp[j(ω_kt+φ_k，n)]

其中a_k，n为信号幅度，与载噪比有关；C_k(t)为CA码，由用户位置、钟偏决定；j为虚数单位；ω_k为载波圆频率，由用户速度、钟漂决定；φ_k，n为载波相位；

将最大似然估计式展开，得到：

由于R只与接收信号有关，因此最小化[R-S]，等效于最大化[S]，即：

其中I_k，n是第k颗卫星第n数据位的时间跨度，内部求和号为非相干累加，外部求和号为各卫星通道融合累加；

CA码序列波形C_k(t)求解式为：

C_k(t)≌C_k[t-b(T₀)-τ_k(T₀)-(t-T₀)τ′_k(T₀)]

其中C_k(t_GNSS)为以GNSS时间为参数的卫星CA码产生函数，C_k(t_GNSS)已知；τ_k(T₀)为T₀到达的第k颗卫星信号传输延时；τ′_k(T₀)为由多普勒频率引起的延迟时间变化率；b(T₀)为钟偏，由导航电文给出，几个参数计算如下：

${\mathrm{\τ}}_k\left(T_0\right)\≅\frac{|p_k\left(T_0\right)-p\left(T_0\right)|}{c+v_k\left(T_0\right)\·u_k\left(T_0\right)}\sec$

${\mathrm{\τ}}_k^{\′}\left(T_0\right)\≅\frac{1}{c}\left\{\right[v_k\left(T_0\right)-v\left(T_0\right)]\·u_k\left(T_0\right)\}\sec /\sec$

其中p_k(T₀)为第k颗卫星位置；p(T₀)为用户位置；v_k(T₀)为第k颗卫星速度；u_k(T₀)为用户到第k颗卫星的投影矢量；v(T₀)为用户速度；c为光速；

基带多普勒频率ω_k计算式为：

${\mathrm{\ω}}_k=2\mathrm{\π}\left\{\frac{1}{\mathrm{\λ}}\right[v\left(T_0\right)-v_k\left(T_0\right)]\·u_k\left(T_0\right)-\stackrel{\·}{b}\left(T_0\right)f\}$

其中f为载波频率，λ为载波波长，可见，已知导航电文、用户速度、钟漂的情况下，可以求解基带多普勒频率ω_k。

6.根据权利要求5所述的组合导航方法，其特征在于，在步骤S07中，按照水平速度、垂直速度、钟漂、水平位置、高度、钟偏的顺序进行最大值的搜索，获得对应T时间内的水平速度误差、垂直速度误差、钟漂误差、水平位置误差、高度误差、钟偏误差。

7.根据权利要求6所述的组合导航方法，其特征在于，所述水平速度误差搜索方法为：

假设第n数据位起始时刻INS输出位置、速度为钟偏和钟漂为以这组参数为中心进行8维网格划分，每个维度划分网格的数量依次为n₁～n₈，且n₁～n₈均为奇数，针对每颗卫星，在码相位、载波多普勒两个维度均分为9份：

$\left[\begin{array}{ccccc}{\tilde{p}}_{n,x}-\frac{n_1-1}{2}{\mathrm{\Δp}}_x& \·\·\·& {\tilde{p}}_{n,x}& \·\·\·& {\tilde{p}}_{n,x}+\frac{n_1-1}{2}{\mathrm{\Δp}}_x\\ {\tilde{p}}_{n,y}-\frac{n_2-1}{2}{\mathrm{\Δp}}_y& \·\·\·& {\tilde{p}}_{n,y}& \·\·\·& {\tilde{p}}_{n,y}+\frac{n_2-1}{2}{\mathrm{\Δp}}_y\\ {\tilde{p}}_{n,z}-\frac{n_3-1}{2}{\mathrm{\Δp}}_z& \·\·\·& {\tilde{p}}_{n,z}& \·\·\·& {\tilde{p}}_{n,z}+\frac{n_3-1}{2}{\mathrm{\Δp}}_z\\ {\tilde{b}}_n-\frac{n_7-1}{2}\mathrm{\Δb}& \·\·\·& {\tilde{b}}_n& \·\·\·& {\tilde{b}}_n+\frac{n_7-1}{2}\mathrm{\Δb}\end{array}\right]\⇒\left[\begin{array}{ccccc}{m}_1-{4\mathrm{\Δm}}_1& \·\·\·& m_1& \·\·\·& m_1+{4\mathrm{\Δm}}_1\\ \·& & \·& & \·\\ \·& \·\·\·& \·& \·\·\·& \·\\ \·& & \·& & \·\\ m_k-{4\mathrm{\Δm}}_k& \·\·\·& m_k& \·\·\·& m_k+{4\mathrm{\Δm}}_k\\ \·& & \·& & \·\\ \·& \·\·\·& \·& \·\·\·& \·\\ \·& & \·& & \·\\ m_N-{4\mathrm{\Δm}}_N& \·\·\·& m_N& \·\·\·& m_N+{4\mathrm{\Δm}}_N\end{array}\right]$

$\left[\begin{array}{ccccc}{\tilde{v}}_{n,x}-\frac{n_4-1}{2}{\mathrm{\Δv}}_x& \·\·\·& {\tilde{v}}_{n,x}& \·\·\·& {\tilde{v}}_{n,x}+\frac{n_4-1}{2}{\mathrm{\Δv}}_x\\ {\tilde{v}}_{n,y}-\frac{n_5-1}{2}{\mathrm{\Δv}}_y& \·\·\·& {\tilde{v}}_{n,y}& \·\·\·& {\tilde{v}}_{n,y}+\frac{n_5-1}{2}{\mathrm{\Δv}}_y\\ {\tilde{v}}_{n,z}-\frac{n_6-1}{2}{\mathrm{\Δv}}_z& \·\·\·& {\tilde{v}}_{n,z}& \·\·\·& {\tilde{v}}_{n,z}+\frac{n_6-1}{2}{\mathrm{\Δv}}_z\\ {\widetilde{\stackrel{\·}{b}}}_n-\frac{n_8-1}{2}\mathrm{\Δ}\stackrel{\·}{b}& \·\·\·& {\widetilde{\stackrel{\·}{b}}}_n& \·\·\·& {\widetilde{\stackrel{\·}{b}}}_n+\frac{n_8-1}{2}\mathrm{\Δ}\stackrel{\·}{b}\end{array}\right]\⇒\left[\begin{array}{ccccc}{\mathrm{\ω}}_1-{4\mathrm{\Δ\ω}}_1& \·\·\·& {\mathrm{\ω}}_1& \·\·\·& {\mathrm{\ω}}_1+{4\mathrm{\Δ\ω}}_1\\ \·& & \·& & \·\\ \·& \·\·\·& \·& \·\·\·& \·\\ \·& & \·& & \·\\ {\mathrm{\ω}}_k-{4\mathrm{\Δ\ω}}_k& \·\·\·& {\mathrm{\ω}}_k& \·\·\·& {\mathrm{\ω}}_k+{4\mathrm{\Δ\ω}}_k\\ \·& & \·& & \·\\ \·& \·\·\·& \·& \·\·\·& \·\\ \·& & \·& & \·\\ {\mathrm{\ω}}_N-{4\mathrm{\Δ\ω}}_N& \·\·\·& {\mathrm{\ω}}_N& \·\·\·& {\mathrm{\ω}}_N+{4\mathrm{\Δ\ω}}_N\end{array}\right]$

其中，m₁，m₂，...，m_N为第二网格中心网格对应的各颗卫星的码相位，ω₁，ω₂，...，ω_N为第二网格中心网格对应的各颗卫星的载波圆频率，

对于卫星k，以每个网格中心对应码相位、载波多普勒频率设置相关器，则每颗卫星需要设置81个相关器，预检测积分时间T后，得到81个相关值，进入搜索阶段；首先按照搜索中心固定位置、垂直速度、钟偏和钟漂，将卫星k的81个相关值逆映射到水平速度搜索二维平面内；依次将其他卫星的相关值逆映射到水平速度搜索二维平面内，之后按照各颗卫星载噪比权重相加，最大值点即对应水平速度误差。

8.根据权利要求7所述的组合导航方法，其特征在于，在水平速度搜索结束后，固定水平速度误差，按照相同的方法依次进行垂直速度、钟漂、水平位置、高度、钟偏误差搜索。

9.根据权利要求3或4所述的组合导航方法，其特征在于，卡尔曼滤波误差模型向量为：

其中，地理坐标系选为北东地坐标系，下标N、E、D分别表示北向，东向，地向，为惯导姿态误差，δV为惯导速度误差，为惯导陀螺常值漂移，为惯导加速度计常值漂移，δL、δλ和δh分别为惯导纬度误差、经度误差和高度误差，b为时钟偏差，为时钟漂移；

卡尔曼滤波状态方程为：

其中：

$F=\left[\begin{array}{cccccc}{F}_{11}& F_{12}& F_{13}& F_{14}& 0_{3\×3}& 0_{3\×2}\\ F_{21}& F_{22}& F_{23}& 0_{3\×3}& F_{25}& 0_{3\×2}\\ 0_{3\×3}& F_{32}& F_{33}& 0_{3\×3}& 0_{3\×3}& 0_{3\×2}\\ 0_{6\×3}& 0_{6\×3}& 0_{6\×3}& 0_{6\×3}& 0_{6\×3}& 0_{3\×2}\\ 0_{2\×3}& 0_{2\×3}& 0_{2\×3}& 0_{2\×3}& 0_{2\×3}& F_{46}\end{array}\right],W=\left[\begin{array}{c}-C_b^n{\mathrm{\ϵ}}_w^b\\ C_b^n{\▿}_w^b\\ 0_{11\×1}\end{array}\right]$

$F_{11}=(-{\mathrm{\ω}}_{\mathrm{in}}^n\×)=\left[\begin{array}{ccc}0& -({\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_E}{R_N+h}\tan L)& \frac{V_N}{R_M+h}\\ {\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_E}{R_N+h}\tan L& 0& {\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E}{R_N+h}\\ -\frac{V_N}{R_M+h}& -({\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E}{R_N+h})& 0\end{array}\right]$

$F_{12}=\left[\begin{array}{ccc}0& \frac{1}{R_N+h}& 0\\ -\frac{1}{R_M+h}& 0& 0\\ 0& -\frac{\tan L}{R_N+h}& 0\end{array}\right]$

$F_{13}=\left[\begin{array}{ccc}-{\mathrm{\ω}}_{\mathrm{ie}}\sin L& 0& 0\\ 0& 0& 0\\ -{\mathrm{\ω}}_{\mathrm{ie}}\cos L-\frac{V_E}{R_N+h}{\sec }^{2}L& 0& 0\end{array}\right]$

$F_{14}=-\left[\begin{array}{ccc}\cos \mathrm{\θ}\cos \mathrm{\ψ}& -\cos \mathrm{\γ}\sin \mathrm{\ψ}+\sin \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}& \sin \mathrm{\γ}\sin \mathrm{\ψ}+\cos \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}\\ \cos \mathrm{\θ}\sin \mathrm{\ψ}& \cos \mathrm{\γ}\cos \mathrm{\ψ}+\sin \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}& -\sin \mathrm{\γ}\cos \mathrm{\ψ}+\cos \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}\\ -\sin \mathrm{\θ}& \sin \mathrm{\γ}\cos \mathrm{\θ}& \cos \mathrm{\γ}\cos \mathrm{\θ}\end{array}\right]$

$F_{21}=(f^n\×)=\left[\begin{array}{ccc}0& -f_z^n& f_y^n\\ f_z^n& 0& -f_x^n\\ -f_y^n& f_x^n& 0\end{array}\right]$

$F_{22}=\left[\begin{array}{ccc}\frac{V_D}{R_M+h}& -2({\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_E}{R_N+h}\tan L)& \frac{V_N}{R_M+h}\\ 2{\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_E}{R_N+h}\tan L& \frac{V_D}{R_N+h}+\frac{V_N}{R_N+h}\tan L& {2\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E}{R_N+h}\\ -\frac{{2V}_N}{R_M+h}& -2({\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E}{R_N+h})& 0\end{array}\right]$

$F_{23}=\left[\begin{array}{ccc}-({2V}_E{\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E^2}{R_N+h}{\sec }^{2}L)& 0& 0\\ 2V_N{\mathrm{\ω}}_{\mathrm{ie}}\cos L-2V_E{\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_N{V}_E}{R_N+h}{\sec }^{2}L& 0& 0\\ 2V_E{\mathrm{\ω}}_{\mathrm{ie}}\sin L& 0& -\frac{2g}{R_M}\end{array}\right]$

$F_{25}=\left[\begin{array}{ccc}\cos \mathrm{\θ}\cos \mathrm{\ψ}& -\cos \mathrm{\γ}\sin \mathrm{\ψ}+\sin \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}& \sin \mathrm{\γ}\sin \mathrm{\ψ}+\cos \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}\\ \cos \mathrm{\θ}\sin \mathrm{\ψ}& \cos \mathrm{\γ}\cos \mathrm{\ψ}+\sin \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}& -\sin \mathrm{\γ}\cos \mathrm{\ψ}+\cos \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}\\ -\sin \mathrm{\θ}& \sin \mathrm{\γ}\cos \mathrm{\θ}& \cos \mathrm{\γ}\cos \mathrm{\θ}\end{array}\right]$

$F_{32}=\left[\begin{array}{ccc}\frac{1}{R_M+h}& 0& 0\\ 0& \frac{1}{(R_N+h)\cos L}& 0\\ 0& 0& -1\end{array}\right]F_{33}=\left[\begin{array}{ccc}0& 0& -\frac{V_N}{{(R_M+h)}^2}\\ \frac{V_E\tan L}{(R_N+h)\cos L}& 0& -\frac{V_E}{{(R_N+h)}^2\cos L}\\ 0& 0& 0\end{array}\right]F_{46}=\left[\begin{array}{cc}0& 1\\ 0& 0\end{array}\right]$

其中，为载体坐标系到导航坐标系的转换矩阵，θ为俯仰角，γ为横滚角，ψ为航向角，V_N、V_E和V_D分别为载体的北向、东向和地向速度，L为载***置的地理纬度，h为载体海拔高度，R_M为地球子午圈曲率半径，R_N为地球卯酉圈曲率半径，ω_ie为地球的自转角速度，和为惯导的比力，为陀螺的噪声，为加速度计的噪声；

卡尔曼滤波器观测向量为： $Z=[\mathrm{\Δx},\mathrm{\Δy},\mathrm{\Δz},\mathrm{\Δ}{v}_N,\mathrm{\Δ}{v}_E,\mathrm{\Δ}{v}_D,\mathrm{\Δb},\mathrm{\Δ}\stackrel{\·}{b}];$

卡尔曼滤波器观测方程为：Z＝HX+V；

其中，

$H=\left[\begin{array}{cccc}{0}_{3\×3}& 0_{3\×3}& I_{3\×3}& 0_{3\×2}\\ 0_{3\×3}& I_{3\×3}& 0_{3\×3}& 0_{3\×2}\\ 0_{2\×3}& 0_{2\×3}& 0_{2\×3}& I_{2\×2}\end{array}\right]$ 为观测矩阵，V为观测噪声。

10.根据权利要求1所述的组合导航方法，其特征在于，在步骤S06中，捷联惯导解算速率与导航电文码率相等或呈倍数关系。

11.根据权利要求1所述的组合导航方法，其特征在于，在步骤S07中，根据卡尔曼滤波器状态协方差矩阵自适应调整搜索网格大小。

12.根据权利要求1所述的组合导航方法，其特征在于，所述GNSS***为GPS、北斗卫星定位***、格洛纳斯***或伽利略导航***。

13.一种GNSS/MINS超深组合导航***，包括：

GNSS***，通过天线(601)接收卫星信号，经射频前端(602)，将所述卫星信号下变频采样后输出中频信号，并接收卫星导航电文；

MINS***，通过MIMU(603)获取测量数据，经捷联惯导***(604)解算后得到载***置、速度及姿态；

载波频率及码相位计算单元(605)，用于将所述MINS***输出的速度和位置结合所述卫星导航电文计算GNSS信号载波频率和码相位，以产生本地信号；以及

组合卡尔曼滤波器(609)，将状态误差反馈回MINS***进行误差修正；

其特征在于，还包括：

相关器组(606)，用于将所述本地信号与中频信号进行相关；

投影映射单元(607)，将所述相关器组(606)的输出结果完成从第二网格到第一网格的投影映射，其中，所述第一网格为在三维位置、三维速度、钟漂、钟偏共八维变量上的网格划分，第二网格为在GNSS载波频率、伪码相位等二维变量上的网格划分；以及

搜索单元(608)，对所述投影映射单元(607)的输出结果进行搜索，获得最大值，即位置误差、速度误差、钟偏误差和钟漂误差，作为所述组合卡尔曼滤波器的观测量。

14.根据权利要求13所述的组合导航***，其特征在于，所述搜索单元(608)按照水平速度、垂直速度、钟漂、水平位置、高度、钟偏的顺序进行最大值的搜索，获得对应T时间内的水平速度误差、垂直速度误差、钟漂误差、水平位置误差、高度误差、钟偏误差。

15.根据权利要求13所述的组合导航***，其特征在于，所述GNSS***为GPS、北斗卫星定位***、格洛纳斯***或伽利略导航***。

16.根据权利要求13所述的组合导航***，其特征在于，所述***还包括气压高度表，为***引入气压高度数据。

17.根据权利要求13所述的组合导航***，其特征在于，所述***还包括高精度晶振。

18.一种应用于权利要求1-12任一项所述的组合导航方法，或者应用于权利要求13-17任一项所述的组合导航***中的组合导航装置，包括：

天线和射频模块，通过FPGA接口将GNSS中频数据引入***；MIMU，通过FPGA接口将MIMU的测量数据引入***；以及电源模块、DSP和上位机；其特征在于，所述组合导航装置还包括相关器组，所述DSP将计算码相位、载波频率和投影矢量，并通过EMIF总线写入FPGA内部，FPGA根据搜索网格大小自动划分网格并映射到所述相关器组，预检测积分时间T后，FPGA处理相关器组的输出值并找到最大值，即对应T时间内的位置误差、速度误差、钟偏和钟漂误差，所述DSP读入所述位置误差、速度误差、钟偏和钟漂误差观测值并对用户状态量进行更新。

19.根据权利要求18所述的组合导航装置，其特征在于，还包括气压高度表，通过FPGA接口将气压高度数据引入***。

20.根据权利要求18所述的组合导航装置，其特征在于，所述GNSS***为GPS、北斗卫星定位***、格洛纳斯***或伽利略导航***。