CN105157703A - 一种重力辅助惯性导航适配区可导航性的评价方法 - Google Patents

Abstract

本发明提出一种重力辅助惯性导航适配区可导航性的评价方法,解决了目前无法评价所选适配区是否有效的技术问题。本发明以分形几何理论为基础，通过分析重力异常曲面的轮廓谱距和表面谱距的特性，提出一种重力场特征度量值——等方性系数，分别得出等方性好与差时，适配区的匹配概率与方向之间的定性关系，有助于提高重力场适配区可导航性分析的效率。

Description

一种重力辅助惯性导航适配区可导航性的评价方法

技术领域

本发明涉及一种重力辅助惯性导航适配区可导航性的评价方法，属于重力辅助惯性导航系统航迹规划领域。

背景技术

水下载体在航行过程中常用重力来辅助惯导进行位置校正，用重力来辅助惯性导航就是利用重力异常值或者重力梯度值等辅助信息来修正惯性导航随时间积累的定位误差，以实现载体的精确导航。重力辅助惯性导航的核心是匹配算法，其中序列匹配算法的匹配精度与局部重力场特征有很大的关系，局部区域重力异常数据离散程度越高，校正精度越高；反之校正精度较差，甚至会发散，因此保证载体可以长航时、高精度航行的关键是提取丰富的适配区使得惯导的定位误差可以得到及时的校正。

目前有很多文献已给出了提取适配区的方法，但其所选的适配区究竟能否有校辅助惯导进行校正，并没有给出具体的评价方法。在实际航行过程中，一片海域内的适配区是固定的，但载体进入海域的方向是不定的，当载体沿不同方向进入适配区时，参与匹配的重力场曲面剖面也不同，得到的匹配效果也不同，在适配区内不一定在所有方向上都能得到好的匹配效果，因此为了全面评价重力场适配区的可导航性，本发明引入重力异常三维曲面的等方性系数来对适配区进行评价，从而有效的利用重力辅助惯导进行位置校正。

发明内容

本发明的目的是为了克服已有技术存在的不足，提出一种重力辅助惯性导航适配区可导航性的评价方法,解决了目前无法评价所选适配区是否有效的技术问题。

本发明的基本原理是：以分形几何理论为基础，通过分析重力异常曲面的轮廓谱距和表面谱距的特性，提出一种重力场特征度量值——等方性系数，分别得出等方性好与差时，适配区的匹配概率与方向之间的定性关系，有助于提高重力场适配区可导航性分析的效率。

本发明的目的是通过下述技术方案实现的。

一种重力辅助惯性导航适配区可导航性的评价方法，包括如下步骤：

步骤一、根据重力场重力异常值确定适配区重力场三维曲面图，利用随机游走法计算其三维曲面图的分形维数；

进一步地，采用如下方法计算三维分形维数，步骤如下：

(1)首先将重力异常值看作随机游走的结果，设定一个间隔值R，计算每个点跟它上下左右相邻R的重力异常值差，即：

P＝g(x₂,y₂)-g(x₁,y₁)(1)

其中P为两坐标点的重力异常值差，g(x,y)为在重力异常网格图上坐标为(x,y)的重力异常值，||(x₂,y₂)-(x₁,y₁)||＝R；

(2)计算P的期望E(P)；

(3)计算E(P)和R的双对数log-log函数,取不同R得到一组对应的E(P)和R，进行线性拟合，得到斜率k,由D＝3-k求得重力异常三维曲面图的分形维数D。

步骤二、计算实际重力场重力异常方差，反推求得重力异常序列特征尺度系数；

用于重力场重力异常序列轮廓的表征，分形维数为D的W-M函数形式如下：

$g(x,y)=G^{D-1}\underset{n=n_1}{\overset{\∞}{\mathrm{\Σ}}}\frac{\cos 2\mathrm{\π}{\mathrm{\γ}}^{n}x}{{\mathrm{\γ}}^{(2-D)n}}---\left(2\right)$

g(x,y)为重力场重力异常数据；D为分形维数；G为特征尺度系数；γ为大于1的常数，γⁿ为随机过程的空间频率；n₁与重力异常面轮廓的最低阶段频率相对应。

进一步地，对于近似服从正态分布的重力异常三维曲面，γ取1.5。

相应于(2)式的连续功率谱为

$P\left(\mathrm{\&omega;}\right)=\frac{G^{2(D-1)}}{2\ln \mathrm{\&gamma;}}\&CenterDot;\frac{1}{{\mathrm{\&omega;}}^{(5-2D)}},{\mathrm{\&gamma;}}^{n_1}<\mathrm{\&omega;}<+\&infin;---\left(3\right)$

其中ω为频率；

根据(3)式得重力异常序列方差表示为

$m_0={\&Integral;}_{w_l}^{w_h}P\left(\mathrm{\&omega;}\right)\mathrm{d\&omega;}=\frac{G^{2(D-1)}}{2\ln \mathrm{\&gamma;}}\&CenterDot;\frac{1}{4-2D}\&CenterDot;[\frac{1}{{\mathrm{\&omega;}}_l^{(4-2D)}}-\frac{1}{{\mathrm{\&omega;}}_h^{(4-2D)}}]---\left(4\right)$

计算实际重力异常的方差得到m₀，根据(4)式求得特征尺度系数G，其中△为取样间隔，L为取样长度。

步骤三、计算重力异常三维曲面图的表面二阶谱距；

重力异常曲面是一个三维均匀的随机表面，γ阶轮廓谱距m_r和表面谱距m_r-q,r定义分别如下：

$m_{\mathrm{\&gamma;}}=\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}{\mathrm{\&omega;}}^{\mathrm{\&gamma;}}P\left(\mathrm{\&omega;}\right)\mathrm{d\&omega;}---\left(5\right)$

$m_{\mathrm{\&gamma;}-q,q}={\&Integral;}_{-\&infin;}^{+\&infin;}{\&Integral;}_{-\&infin;}^{+\&infin;}{\mathrm{\&omega;}}_1^{\mathrm{\&gamma;}-q}{\mathrm{\&omega;}}_2^{q}Q({\mathrm{\&omega;}}_1,{\mathrm{\&omega;}}_2)d{\mathrm{\&omega;}}_{1}d{\mathrm{\&omega;}}_2---\left(6\right)$

其中Q(ω₁,ω₂)为重力异常轮廓的二元频谱函数，根据Longuet-Higgins等式，设定一个参考方向x，与参考方向成θ角方向上的γ阶轮廓谱距m_r与表面谱距m_r-q,q有如下关系：

$m_{\mathrm{\&gamma;}}\left(\mathrm{\&theta;}\right)=m_{\mathrm{\&gamma;}-q,q}{C}_r^q{\cos }^{r-q}\mathrm{\&theta;}{\sin }^q\mathrm{\&theta;}---\left(7\right)$

由(7)式推出在θ方向上轮廓的二阶谱距m₂(θ)为

m₂(θ)＝m₂₀cos²θ+2m₁₁cosθsinθ+m₀₂sin²θ(8)

根据(8)式，分别在三个方向上测量重力异常数据，可以求得重力异常序列轮廓表面的二阶谱距，根据(9)式可以求得重力异常三维曲面的表面谱距：

$\left[\begin{array}{c}{m}_{20}\\ m_{11}\\ m_{02}\end{array}\right]=T_2^{-1}\left[\begin{array}{c}{m}_2\left({\mathrm{\&theta;}}_1\right)\\ m_2\left({\mathrm{\&theta;}}_2\right)\\ m_2\left({\mathrm{\&theta;}}_3\right)\end{array}\right]---\left(9\right)$

其中 $T_2=\left[\begin{array}{ccc}{\cos }^2{\mathrm{\&theta;}}_1& 2\cos {\mathrm{\&theta;}}_1\sin {\mathrm{\&theta;}}_1& {\sin }^2{\mathrm{\&theta;}}_1\\ {\cos }^2{\mathrm{\&theta;}}_2& 2\cos {\mathrm{\&theta;}}_2{\sin \mathrm{\&theta;}}_2& {\sin }^2{\mathrm{\&theta;}}_2\\ {\cos }^2{\mathrm{\&theta;}}_2& 2\cos {\mathrm{\&theta;}}_3\sin {\mathrm{\&theta;}}_3& {\sin }^2{\mathrm{\&theta;}}_3\end{array}\right],$ θ₁，θ₂，θ₃为任选的与参考方向成一定夹角的三个方向，表面谱距m₂₀和m₀₂为相互垂直的两个方向上的斜率方差，而m₁₁为两个方向上的协方差；

步骤四、根据重力异常三维曲面的表面谱距计算其等方性系数，进行适配区可导航性的评价，

测量坐标系旋转时，该随机表面各条轮廓高度方向的概率分布不变。轮廓的二阶谱矩，既可以反映轮廓幅度分布，又可以反映频率特性，而且一般情况下很强地依赖于被测轮廓的方向，因此可以用其来表征各向异性表面的等方性。如果三维表面为等方性的，则必定m₂(θ)与θ无关，此时m₂(θ)的轨迹为圆。

由m₂(θ)＝m₂₀cos²θ+2m₁₁cosθsinθ+m₀₂sin²θ可知轮廓的二阶谱距不仅依赖于θ。对于均一的随机表面，轮廓斜率方差最大值和最小值所在的方向是相互垂直的，这两个方向称为主方向，而斜率方向最大的那条轮廓曲线与x夹角称为主方向θ_p。

对(8)式求导求取最大值和最小值可以求得

$m_{2\max }=\frac{1}{2}[M_2+\sqrt{M_2^2-4{\mathrm{\&Delta;}}_2}],m_{2\min }=\frac{1}{2}[M_2-\sqrt{M_2^2-4{\mathrm{\&Delta;}}_2}]\backslash *\mathrm{MERGEFORMAT}---\left(10\right)$

其中M₂＝m₂₀+m₀₂,若m₂(θ)与θ无关，则m_2max＝m_2min，则

因此若三维表面是完全等方性的，则我们用下式评定各向异性表面的等方性

$\mathrm{\&Delta;}=2\frac{\sqrt{{\mathrm{\&Delta;}}_2}}{M_2}=\frac{2\sqrt{m_{20}\&CenterDot;m_{02}-m_{11}^2}}{m_{20}+m_{02}}\backslash *\mathrm{MERGEFORMAT}---\left(11\right)$

当m₁₁＝0,m₂₀＝m₀₂,△＝1时，表示重力异常三维曲面是完全等方性的，所以△越接近1，等方性越好。实际评价时，对于同一重力异常曲面，测量的网格数不仅仅代表三条轮廓线，这时△不为单一数值，等方性参数可取其平均值 ${\mathrm{\&Delta;}}_{\mathrm{average}}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_i.$

根据等方性系数较低和较高的情况下特定区域上的适配性度量值和方向之间的关系，在进行导航任务的预处理时，可以先计算某方向上适配区的等方性系数，如果该适配区等方性较好，则它的导航性能对方向不敏感，各个方向上的匹配概率相差很小，不需要再经过复杂的分析就可以全面得出该区域的可导航性；如果该区域等方性较差，则它的导航性能对方向敏感，适配区在各个方向上的可导航性能可能不一样，并且各个方向上的匹配概率起伏很大，只有选择较小的方向粒度作为步调再计算分析该区域在各个方向上的匹配性能指标，才能全面分析该适配区的可导航性，按照上述策略可以大大提高导航规划系统预处理的效率。

本发明的有益效果：

本发明方法给出了一种重力辅助惯性导航适配区可导航性的评价方法，利用分形几何的原理，通过分析重力异常曲面的轮廓谱距和表面谱距的特性，提出用等方性系数来衡量适配区的可导航性，分别得出等方性好与差时，适配区的匹配概率与方向之间的定性关系，有助于提高重力场适配区可导航性分析的效率。

具体实施方式

本实施例中，取10°×10°的海洋范围，分辨率为10'×10'，局部区域的网格取m＝n＝5，也即5×5的局部区域网格，重力场背景数据为重力异常值。其过程如下：

步骤一、根据重力场重力异常值确定适配区重力场三维曲面图，利用随机游走法计算其三维分形维数；

(1)首先将重力异常值看作随机游走的结果，设定一个间隔值R(如R＝3)，计算每个点跟它上下左右相邻R的重力异常值差，即：

P＝g(x₂,y₂)-g(x₁,y₁)(1)

其中P为两坐标点的重力异常值差，g(x,y)为在重力异常网格图上坐标为(x,y)的重力异常值，||(x₂,y₂)-(x₁,y₁)||＝R；

(2)计算P的期望E(P)；

(3)计算E(P)和R的双对数log-log函数，即log(E)＝(3-D)logR+C,取不同R得到一组对应的E(P)和R，进行线性拟合，得到斜率k,由D＝3-k求得重力异常三维曲面图的分形维数D；

步骤二、对比实际重力场重力异常的方差，求重力异常序列特征尺度系数

分形维数为D的W-M函数形式如下：

$g(x,y)=G^{D-1}\underset{n=n_1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\frac{\cos 2\mathrm{\&pi;}{\mathrm{\&gamma;}}^{n}x}{{\mathrm{\&gamma;}}^{(2-D)n}}---\left(2\right)$

g(x,y)为重力场重力异常数据；D为分形维数；G为特征尺度系数；γ为大于1的常数，对于近似服从正态分布的重力异常三维曲面，γ取1.5。γⁿ为随机过程的空间频率；n₁与重力异常面轮廓的最低阶段频率相对应。

相应于(2)式的连续功率谱为

$P\left(\mathrm{\&omega;}\right)=\frac{G^{2(D-1)}}{2\ln \mathrm{\&gamma;}}\&CenterDot;\frac{1}{{\mathrm{\&omega;}}^{(5-2D)}},{\mathrm{\&gamma;}}^{n_1}<\mathrm{\&omega;}<+\&infin;---\left(3\right)$

根据(3)式可得重力异常序列方差表示为

$m_0={\&Integral;}_{w_l}^{w_h}P\left(\mathrm{\&omega;}\right)\mathrm{d\&omega;}=\frac{G^{2(D-1)}}{2\ln \mathrm{\&gamma;}}\&CenterDot;\frac{1}{4-2D}\&CenterDot;[\frac{1}{{\mathrm{\&omega;}}_l^{(4-2D)}}-\frac{1}{{\mathrm{\&omega;}}_h^{(4-2D)}}]---\left(4\right)$

计算实际重力异常的方差可以得到m₀，根据(4)式可以求得特征尺度系数G,其中△为取样间隔，L为取样长度。

步骤三、计算重力异常三维曲面图的表面二阶谱距；

对于三维均匀的随机表面，γ阶轮廓谱距m_r和表面谱距m_r-q,r定义分别如下：

$m_r=\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}{\mathrm{\&omega;}}^{r}P\left(\mathrm{\&omega;}\right)\mathrm{d\&omega;}---\left(5\right)$

$m_{r-q,r}={\&Integral;}_{-\&infin;}^{+\&infin;}{\&Integral;}_{-\&infin;}^{+\&infin;}{\mathrm{\&omega;}}_1^{r-q}{\mathrm{\&omega;}}_2^{q}P({\mathrm{\&omega;}}_1,{\mathrm{\&omega;}}_2)d{\mathrm{\&omega;}}_{1}d{\mathrm{\&omega;}}_2---\left(6\right)$

根据Longuet-Higgins等式，θ方向上的γ阶轮廓谱距m_r与表面谱距m_r-q,r有如下关系：

$m_r\left(\mathrm{\&theta;}\right)=m_{r-q,r}{C}_r^q{\cos }^{r-q}\mathrm{\&theta;}{\sin }^q\mathrm{\&theta;}---\left(7\right)$

由(7)式可以推出在θ方向上轮廓的二阶谱距m₂(θ)为

m₂(θ)＝m₂₀cos²θ+2m₁₁cosθsinθ+m₀₂sin²θ(8)

根据(8)式，分别在三个方向上测量重力异常数据，可以求得重力异常序列轮廓表面的二阶谱距,根据(9)式可求得重力异常三维曲面二阶表面谱距：

$\left[\begin{array}{c}{m}_{20}\\ m_{11}\\ m_{02}\end{array}\right]=T_2^{-1}\left[\begin{array}{c}{m}_2\left({\mathrm{\&theta;}}_1\right)\\ m_2\left({\mathrm{\&theta;}}_2\right)\\ m_2\left({\mathrm{\&theta;}}_3\right)\end{array}\right]---\left(9\right)$

其中 $T_2=\left[\begin{array}{ccc}{\cos }^2{\mathrm{\&theta;}}_1& 2\cos {\mathrm{\&theta;}}_1\sin {\mathrm{\&theta;}}_1& {\sin }^2{\mathrm{\&theta;}}_1\\ {\cos }^2{\mathrm{\&theta;}}_2& 2\cos {\mathrm{\&theta;}}_2{\sin \mathrm{\&theta;}}_2& {\sin }^2{\mathrm{\&theta;}}_2\\ {\cos }^2{\mathrm{\&theta;}}_2& 2\cos {\mathrm{\&theta;}}_3\sin {\mathrm{\&theta;}}_3& {\sin }^2{\mathrm{\&theta;}}_3\end{array}\right],$ 表面谱距m₂₀和m₀₂为相互垂直的两个方向上的斜率方差，而m₁₁为两个方向上的协方差。

步骤四、根据重力异常三维曲面的表面谱距计算其等方性系数，进行适配区可导航性的评价；

随机表面的等方性为：测量坐标系旋转时，该随机表面各条轮廓高度方向的概率分布不变。轮廓的二阶谱矩，既可以反映轮廓幅度分布，又可以反映频率特性，而且一般情况下很强地依赖于被测轮廓的方向，因此可以用其来表征各向异性表面的等方性。如果三维表面为等方性的，则必定m₂(θ)与θ无关，此时m₂(θ)的轨迹为圆。

我们用下式评定各向异性表面的等方性

$\mathrm{\&Delta;}=2\frac{\sqrt{{\mathrm{\&Delta;}}_2}}{M_2}=\frac{2\sqrt{m_{20}\&CenterDot;m_{02}-m_{11}^2}}{m_{20}+m_{02}}---\left(11\right)$

当m₁₁＝0,m₂₀＝m₀₂,△＝1时，表示重力异常三维曲面是完全等方性的，所以△越接近1，等方性越好。实际评价时，对于同一重力异常曲面，测量的网格数不仅仅代表三条轮廓线，这时△不为单一数值，等方性参数可取其平均值 ${\mathrm{\&Delta;}}_{\mathrm{average}}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}_i.$

为了说明本发明的效果，利用本发明方法在事先选取的适配区上进行可导航性分析，得到适配区的等方性系数，在适配区内以经纬度最小值为坐标原点，在0°～90°上进行匹配试验得到不同的匹配率，在45°～60°范围内两种匹配方法可以得到较好的匹配效果，该适配区的可导航性能一般，在45°～60°度方向上适合进行匹配。

综上所述，本发明方法能有效的评价所选的适配区的可导航性。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进，或者对其中部分技术特征进行等同替换，这些改进和替换也应视为本发明的保护范围。

Claims (5)

1.一种重力辅助惯性导航适配区可导航性的评价方法，其特征在于，包括如下步骤：

步骤一、根据重力场重力异常值确定适配区重力场三维曲面图，利用随机游走法计算其三维曲面图的分形维数；

步骤二、计算实际重力场重力异常方差，反推求得重力异常序列特征尺度系数；

步骤三、计算重力异常三维曲面图的表面二阶谱距；

步骤四、根据重力异常三维曲面的表面谱距计算其等方性系数，进行适配区可导航性的评价。

2.如权利要求1所述的一种重力辅助惯性导航适配区可导航性的评价方法，其特征在于，进一步地，采用如下方法计算三维曲面图的分形维数，步骤如下：

(1)首先将重力异常值看作随机游走的结果，设定一个间隔值R，计算每个点跟它上下左右相邻R的重力异常值差，即：

P＝g(x₂,y₂)-g(x₁,y₁)(1)

其中P为两坐标点的重力异常值差，g(x,y)为在重力异常网格图上坐标为(x,y)的重力异常值，||(x₂,y₂)-(x₁,y₁)||＝R；

(2)计算P的期望E(P)；

(3)计算E(P)和R的双对数log-log函数,取不同R得到一组对应的E(P)和R，进行线性拟合，得到斜率k,由D＝3-k求得重力异常三维曲面图的分形维数D。

3.如权利要求2所述的一种重力辅助惯性导航适配区可导航性的评价方法，其特征在于，采用下述方法求得重力异常序列特征尺度系数；

用于重力场重力异常序列轮廓的表征，分形维数为D的W-M函数形式如下：

$g(x,y)=G^{D-1}\underset{n=n_1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\frac{\cos {2\mathrm{\&pi;\&gamma;}}^{n}x}{{\mathrm{\&gamma;}}^{(2-D)n}}---\left(2\right)$

g(x,y)为重力场重力异常数据；D为分形维数；G为特征尺度系数；γ为大于1的常数，γⁿ为随机过程的空间频率；n₁与重力异常面轮廓的最低阶段频率相对应；

相应于(2)式的连续功率谱为

$P\left(\mathrm{\&omega;}\right)=\frac{G^{2(D-1)}}{2\ln \mathrm{\&gamma;}}\&CenterDot;\frac{1}{{\mathrm{\&omega;}}^{(5-2D)}},{\mathrm{\&gamma;}}^{n_1}<\mathrm{\&omega;}<+\&infin;---\left(3\right)$

其中ω为频率；

根据(3)式得重力异常序列方差表示为

$m_0={\&Integral;}_{w_l}^{w_h}P\left(\mathrm{\&omega;}\right)\mathrm{d\&omega;}=\frac{G^{2(D-1)}}{2\ln \mathrm{\&gamma;}}\&CenterDot;\frac{1}{4-2D}\&CenterDot;[\frac{1}{{\mathrm{\&omega;}}_l^{(4-2D)}}-\frac{1}{{\mathrm{\&omega;}}_h^{(4-2D)}}]---\left(4\right)$

计算实际重力异常的方差得到m₀，根据(4)式求得特征尺度系数G，其中 △为取样间隔，L为取样长度。

4.如权利要求3所述的一种重力辅助惯性导航适配区可导航性的评价方法，其特征在于，进一步地，对于近似服从正态分布的重力异常三维曲面，γ取1.5。

5.如权利要求1或2或3或4所述的一种重力辅助惯性导航适配区可导航性的评价方法，其特征在于，根据等方性系数较低和较高的情况下特定区域上的适配性度量值和方向之间的关系，在进行导航任务的预处理时，先计算某方向上适配区的等方性系数，如果该适配区等方性较好，则它的导航性能对方向不敏感，各个方向上的匹配概率相差很小，不需要再经过复杂的分析就可以全面得出该区域的可导航性；如果该区域等方性较差，则它的导航性能对方向敏感，适配区在各个方向上的可导航性能可能不一样，并且各个方向上的匹配概率起伏很大，只选择较小的方向粒度作为步调再计算分析该区域在各个方向上的匹配性能指标，全面分析该适配区的可导航性，从而提高导航规划系统预处理的效率。