Research on fractal interpolation of gravity anomaly random field for underwater aided navigation
-
摘要: 针对重力异常数据匮乏的地区，本文提出采用分形布朗运动(fractional Brownian motion，FBM)分形插值方法对重力异常随机场进行重构，实现可用于水下重力辅助导航的重力异常基准图生成. 实验结果表明：分形插值方法作为一种粗插值是可行的，可为水下重力辅助导航系统的实现提供理论参考.
-
关键词:
- 水下重力辅助导航 /
- 重力异常 /
- 随机场重构 /
- 分形布朗运动(FBM) /
- 分形插值
Abstract: For areas with scarce gravity anomaly data, this paper proposes to use fractal Brownian motion (FBM) fractal interpolation method to reconstruct the gravity anomaly random field, achieving the generation of gravity anomaly reference maps that can be used for underwater gravity assisted navigation. The experimental results indicate that the fractal interpolation method is feasible as a rough interpolation method and can provide theoretical reference for the implementation of underwater gravity assisted navigation systems. -
0. 引 言
地球物理场和惯性导航系统(inertial navigation system，INS)联合的辅助导航技术可有效遏制INS误差累积问题，始终是国内外水下导航领域的研究难点和热点. 水下辅助导航技术主要包括地磁匹配导航、地形匹配导航和重力匹配导航等. 而水下重力匹配导航通过测量重力信息与地球重力场进行匹配实现导航定位，在重力辅助导航系统运行时，采集重力特征显著区域的重力信息与重力基准图对比得到位置信息，实现天空海一体化水下潜器惯性/重力组合导航系统重调[1-4].
地球形状的不规则性和密度的不均匀性导致地球各点的重力场并不相同，表现为空间位置(经度、纬度、高度)的函数[5]. 因此潜航器可在航行过程中，通过重力测量仪器采集航线途经的重力数据，并与预先存储的重力数据相匹配，获得潜航器当前定位信息，进而对INS累积的位置误差进行校正[6-7]. 重力辅助导航系统在测量重力场数据过程中，潜航器无需露出或接近水面，测量仪器也无需向外部发射信号或接收外部信号，系统可进行无源、隐蔽的导航定位，潜航器在卫星、无线电定位系统失效或遭受破坏的特殊情况下仍可达到自主隐蔽导航的目的[8-10].
随着重力测量仪器及空间测量技术的发展[11]，在全球范围内快速、准确地获取重力数据成为现实，使得重力辅助导航系统具备校正INS积累位置误差的能力[12-13]. 重力辅助导航系统主要有两种代表性系统，均由Bell Aerospace研发，分别为于1990年研制出的重力梯度导航系统[14]，以及于1991年研制出的重力辅助INS [15]. 重力辅助导航系统由INS、重力测量仪器、数字重力基准图[16]和匹配定位算法等四个主要部分组成，各组成部分的性能差异对重力辅助定位系统的性能有重要影响. 数字重力基准图(重力模型)是重力辅助定位的基础[17-18]，重力基准图描述是否准确，包含的重力特征是否丰富，分辨率是否满足需求，都会影响重力辅助导航系统的性能[19].
为减小线性化误差，对重力异常参考地图进行预处理，选择合适的计算模型对地图进行插值逼近，得到高分辨率、小格网地图，以减小线性化误差对位置观测精度的影响. 在重力异常数据匮乏的地区，随着相邻点间的插值数的增加，拟合函数和Kriging方法无法满足逼近精度要求. 其中，拟合函数法跟插值的规模、控制点数目、控制点统计特征和地图统计特征差异、地图统计特征相关，不适宜重力异常数据匮乏的区域，而Kriging方法不适合大尺度条件下的重力插值，因此采用分形重构模拟真实的重力环境，为系统提供良好的仿真环境. 分形通常被用在模拟真实地表和其他自然形态的描述中，它用分维数表达数据场的统计特征. 与Kriging方法的空间变异函数相似，低分辨率地图不利于长距离的统计特征描述[20-22]. 本文考虑在这种条件下，提出利用分形插值的方法对重力异常随机数据场进行重构，重构模拟将为水下重力辅助INS提供更接近真实的重力环境[23].
分形就是对不规则和支离破碎的自然形状的描述，分形理论的发展促进了利用它对各种自然形态描述的研究. Goodchild在1980年考虑用分形理论分析地文表面，后来Clark和Schweitzer定义了具有更好鲁棒性的对地形表面分形估计. 从1991年开始，Kenichi Arakawa和Eric Krotkov对“利用分形重构自然地形”进行了研究，利用分形布朗运动(fractal Brownian motion, FBM)在不规则采样点集上估计自然地形的分形结构，并以此为基础构造出自然地形，Naokazu等给出基于FBM模型的插值方法并不断进行完善[24].
本文提出采用FBM分形插值方法对重力异常随机场进行重构，构造可用于水下重力辅助导航系统的重力异常基准图，通过实验验证分形插值方法的可行性，将为水下重力辅助导航系统的实现提供理论参考.
1. FBM函数及分维数估计
布朗运动(Brownian Motion，BM)是1982年英国植物学家R. Brown发现的，1923年德国数学家N. Wiener建立了BM的数学模型. BM是一种随机过程
$ B\left(t\right) $ ，满足以下两个条件：$$ B\left(0\right)=C $$ (1) $$ B(b)-B\left(a\right)\sim\mathrm{N}\left(0,\left(b-a\right){\sigma }^{2}\right),a\leqslant b $$ (2) 式中：
$ C $ 为常数；随机过程的任意两点之间($ B\left(b\right)- B\left(a\right) $ 和$ B\left(c\right)-B\left(b\right),a\leqslant b\leqslant c $ ) 的增量是相互独立的.$ B\left(t\right) $ 被写成$$ B\left(rt\right)={r}^{\frac{1}{2}}B\left(t\right) $$ (3) 显然，它具有自紧密性.
Mandelbrot提出FBM函数，由它产生BM过程.
$$ {B}_{H}\left(0\right)=C $$ (4) $$ \begin{split} {B}_{H}(t)-{B}_{H}(0)=&C\left[\int_{0}^{\infty}({(t-s)}^{H-0.5}-{(-s)}^{H-0.5})\mathrm{d}B(S)\right.\\&+\left.{\int_{0}^{\infty}}{(t-s)}^{H-0.5}\mathrm{d}B(S)\right] \end{split} $$ (5) Penland定义了FBM函数统计自相似性特征
$$ {P}_{r}\left\{\frac{f\left(t+\Delta t\right)-f\left(t\right)}{{\Delta t}^{H}} < x\right\}=g\left(x\right) $$ (6) 式中：
$ f\left(t\right) $ 为FBM函数；$ g(x) $ 为累计分布函数($ g(x)\sim {\mathrm{N}}(0,{\sigma }^{2}) $ )；$ H $ 为自仿射参数，它与$ f(t) $ 的分维数$ D $ 之间的关系为$$ D=n+1-H $$ (7) $ n $ 是几何拓朴维数，在重力异常地图中，$ n=2 $ .由式(6)可得
$$ E\left[\frac{f\left(t+\Delta t\right)-f\left(t\right)}{{\|\Delta t\|}^{H}}\right]=E\left[g\left(x\right)\right] $$ (8) 又
$$ E\left[\left|g\left(x\right)\right|\right]=2{\int }_{0}^{\infty }x\frac{1}{{\left(2{\text{π}}\right)}^{1/2}\sigma }{{\mathrm{e}}}^{-{x}^{2}/2{\sigma }^{2}}\mathrm{d}x=\frac{2\sigma }{\sqrt{2{\text{π}}}}=C $$ (9) 可以得到
$$ E\left[\left|f\left(t+\Delta t\right)-f\left(t\right)\right|\right]=C\cdot {\|\Delta t\|}^{H} $$ (10) 两边取对数
$$ \mathrm{log}\;E\left[\left|f\left(t+\Delta t\right)-f\left(t\right)\right|\right]=H\;\mathrm{log}\|\Delta t\|+\mathrm{log}\;C $$ (11) 以
$ \mathrm{log}\left[\left|f\left(t+\Delta t\right)-f\left(t\right)\right|\right] $ 为应变量，$ \mathrm{log}\|\Delta t \| $ 为自变量，则式(11)为直线方程($ H $ 和$ C $ 均为常数).在对二维随机数据场进行分维数估计时，
$ \Delta t $ 用地图上两点距离$ \Delta X=\sqrt{\Delta {x}^{2}+\Delta {y}^{2}} $ 代替，则$$ \mathrm{log}\;E\left[\left|f\left(X+\Delta X\right)-f\left(X\right)\right|\right]=H\;\mathrm{log} \|\Delta X \|+\mathrm{log}\;C $$ (12) 不同距离下的增量的
$ E\left[\left|f\left(X+\Delta X\right)-f\left(X\right)\right|\right] $ 可以通过统计计算得到，最后，按最小二乘法拟合出直线方程即可通过简单变换得到地图的$ H $ 和$ \sigma $ 值.2. 重力异常地图分维数估计
重力异常地图的分维数的估计也是以式(11)为基础进行的.
在计算时需要得到在不同
$ \Delta X $ 下的$ E[|f(X+\Delta X)- f(X)|] $ ，显然在重力异常地图阵列的基础上，可以得到的$ \Delta X $ 是离散的，它的取值与地图网格单位距离$ \left[\Delta {x}_{u},\Delta {y}_{u}\right] $ 和地图阵列的规模$ \left[{m}_{x},{m}_{y}\right] $ 有关，如图1所示.能取的距离为
$ \Delta X $ $$ \Delta X=\sqrt{{(i\cdot {u}_{x})}^{2}+{(j\cdot {u}_{y})}^{2}} $$ (13) 式中，
$ i=\mathrm{1,2},\cdots ,{m}_{x} $ ；$ j=\mathrm{1,2},\cdots ,{m}_{y} $ ；这里仅讨论$ {u}_{x}= {u}_{y}=u $ 的情形.显然，
$ \Delta X $ 增加，在计算$ E\left[\left|f\left(X+\Delta X\right)-f\left(X\right)\right|\right] $ 时，可参与统计点数目会急剧减少. 基于该思路，设定距离容差，可以更好地反映地图的统计特征.即，设
$ {\varepsilon }_{x} $ 为距离容差，当地图上两点距离满足$$ \left|\Delta {X}_{P}-\Delta X\right| < {\varepsilon }_{x} $$ (14) 那么，两点重力异常增量可以计入与
$ \Delta X $ 对应的$ E\left[\left|f\left(X+\Delta X\right)-f\left(X\right)\right|\right] $ 计算中. 其中，$ \Delta {X}_{P} $ 是两点之间的实际网格距离.按此方法分别对2个地块区域的重力异常地图进行了分维数计算，结果如图2所示.
用最小二乘法拟合直线，得到的结果如表1所示.
表 1 分形特征计算结果地块区域 拟合距离/
kmH Df $ \mathrm{log}\;C $ $ \sigma $ 残差/km 161×79 15 0.976 2.024 −0.803 0.197 0.019 60×60 5 0.870 2.130 0.975 11.821 0.017 3. FBM分形插值方法
用FBM分形插值，生成新的小尺度重力异常地图的基本步骤如图3所示.
分形特征一个有限的尺度范围内才能充分体现，这个范围称为无标度区间
$ \left[\Delta {X}_{\mathrm{m}\mathrm{i}\mathrm{n}},\Delta {X}_{\mathrm{m}\mathrm{a}\mathrm{x}}\right] $ .由地图得到的分形特征无标度区间的下限由地图间隔决定
$$ \Delta {X}_{\mathrm{m}\mathrm{i}\mathrm{n}}=u $$ (15) 其上限为分形特征曲线线性拟合区域的上限.
分形插值后得到地图的最小尺度
$ {u}_{c} $ 小于$ u $ ，因此分形插值所能得到的最小尺度不能由原始地图确定. 本文所讨论的重力异常地图分形插值只能用于重力异常地图的粗插值.分形插值采用随机中点位移法，即插入点附近区域的点理解为一个符合
$ (0,\sigma^{2}) $ 的随机数据与区域的平均值的和，它的实现如图4所示.在图4中，填充点为原始地图的控制点，虚线是插值点；插值点分为两类，标号为1的插值点
$ (x,y) $ ，由其周围的4个原地图控制点为参照进行插值$$ \begin{split}g(x,y)= & \frac{1}{4}[g(x-u_c,y-u_c)+g(x+u_c,y-u_c) \\ & +g(x-u_c,y+u_c)+g(x+u_c,y+u_c)] \\ & +\sqrt{1-2^{2\mathrm{\mathit{H}}-2}}\cdot u_c^H\cdot\sigma\cdot\mathrm{Gauss}\end{split} $$ (16) 标号为2的插值点，所用到的参照点由两个原地图控制点和标号为1的两个插值点构成
$$ \begin{split} g(x,y)=&\frac{1}{4}[g(x-{u}_{c},y)+g(x,y-{u}_{c})+g(x+{u}_{c},y)\\&+g(x,y+{u}_{c})]+\sqrt{1-{2}^{2H-2}}\\&\cdot {\left(\frac{{u}_{c}}{\sqrt{2}}\right)}^{H}\cdot \sigma \cdot {\mathrm{Gauss}} \end{split} $$ (17) 其中，
$ {\mathrm{Gauss}} $ 是符合均值为0，标准差为1的随机信号.4. 实验结果与分析
4.1 与原始地图比较
针对某块
$ 60\times 60 $ 的重力异常地图进行插值，在原始地图上进行等间隔抽样得到一块$ 30\times 30 $ 的地图，然后用上述分形插值方法进行插值，并对插值地图和原始地图进行比较. 抽样后的地图分形特征如表2所示.表 2 分形特征参数表地块区域 拟合距离/km $ H $ $ D_f$ $ \mathrm{log}\;C $ $ \sigma $ 残差/km 30×30 6 0.751 2.249 1.026 13.309 0.0150 119×119 5 0.840 2.160 0.992 12.311 0.0280 插值后的重力异常地图与原始地图的比较如图5所示.
插值后的精度分析如表3所示，采用误差平均值(
$ E\left({g}_{e}\right) $ )、绝对误差平均值($ E\left(\left|{g}_{e}\right|\right) $ )、绝对误差最大值($ \mathrm{m}\mathrm{a}\mathrm{x}\left(\left|{g}_{e}\right|\right) $ )、误差均方值($ E\left({{g}_{e}}^{2}\right) $ )和绝对误差标准差($ D\left({g}_{e}\right) $ )等参数来评定.表 3 分形插值精度分析$ E\left({g}_{e}\right) $/km $ E\left(\left|{g}_{e}\right|\right) $/km $ \mathrm{m}\mathrm{a}\mathrm{x}\left(\left|{g}_{e}\right|\right) $/km $ E\left({{g}_{e}}^{2}\right) $/km2 $ D\left({g}_{e}\right) $/km 0.014 6.470 53.419 98.316 7.514 4.2 与原重力异常地图比较
在原重力异常地图上进行插值，分析插值后地图的分形特征的变化. 图6是插值后的等值线图比较. 图7是插值后地图分形特征曲线. 其特征参数见表2第2行数据.
显然，插值后的地图与插值前地图的分形特征变化很小. 分形插值的精度与标准差相关，标准差越大，精度越低；分形插值引入随机信号，因此插值点的重力异常值具有不确定性；采用中点位移法的插值，插值点位置固定，即
$ {u}_{c}=u/2 $ .分形插值作为地形自然表面的重构方法，具有很强的真实感. 重力异常与地形之间的关联性，产生了利用分形对重力异常进行插值的联想，为得到可靠的重力异常地图的分维数，在分维数定义中加入距离容差，改善了分维数在长空间距离下的统计特征描述. 分形插值的精度虽然低，但插值后的地图仍保持原重力异常参考地图的统计特征，该方法为仿真提供了很好的环境模拟. 由于条件的限制，虽然没有对插值结果进行实际测量的验证，但从计算机仿真实验结果看，分形插值作为一种粗插值是可行的.
本文是针对水下重力数据匮乏区域的重力场仿真研究，由于数据异常和条件差，分形插值的实验结果仅与原始数据的计算结果进行了对比研究. 由于重力数据的匮乏，采用其他插值方法的效果不太理想，限于篇幅和数据条件，未进行详细分析，但后续仍需要通过获取大量实测数据来进行反复验证与探讨.
5. 结束语
本文针对重力异常数据匮乏的地区，提出了采用FBM分形插值方法对重力异常随机场进行重构，实现可用于水下重力辅助导航的重力异常基准图生成. 实验结果表明，分形插值后的地图与插值前地图的分形特征变化很小. 分形插值的精度与标准差相关，标准差越大，精度越低；分形插值引入随机信号，因此插值点的重力异常值具有不确定性；采用中点位移法的插值，插值点位置可固定. 分形插值作为地形自然表面的重构方法，具有很强的真实感. 从数值模拟结果来看，分形插值方法作为一种粗插值是可行的和有效的，可为水下重力辅助INS提供更接近真实的重力环境，为水下重力辅助导航系统的实现提供理参考.
-
表 1 分形特征计算结果
地块区域 拟合距离/
kmH Df $ \mathrm{log}\;C $ $ \sigma $ 残差/km 161×79 15 0.976 2.024 −0.803 0.197 0.019 60×60 5 0.870 2.130 0.975 11.821 0.017 
下载: 导出CSV
表 2 分形特征参数表
地块区域 拟合距离/km $ H $ $ D_f$ $ \mathrm{log}\;C $ $ \sigma $ 残差/km 30×30 6 0.751 2.249 1.026 13.309 0.0150 119×119 5 0.840 2.160 0.992 12.311 0.0280
下载: 导出CSV
表 3 分形插值精度分析
$ E\left({g}_{e}\right) $/km $ E\left(\left|{g}_{e}\right|\right) $/km $ \mathrm{m}\mathrm{a}\mathrm{x}\left(\left|{g}_{e}\right|\right) $/km $ E\left({{g}_{e}}^{2}\right) $/km2 $ D\left({g}_{e}\right) $/km 0.014 6.470 53.419 98.316 7.514
下载: 导出CSV
-
[1] 郑伟, 付滕达, 李钊伟, 等. 潜器水下重力匹配导航研究进展[J]. 科学技术与工程, 2022, 22(18): 7721-7734. DOI: 10.3969/j.issn.1671-1815.2022.18.001 [2] LI D, XU J N, HE H Y, et al. An underwater integrated navigation algorithm to deal with DVL malfunctions based on deep learning[J]. IEEE access, 2021, 9(1): 82010-82020. DOI: 10.1109/ACCESS.2021.3083493
[3] MU X K, HE B, WU S Y, et al. A practical INS/GPS/DVL/PS integrated navigation algorithm and its application on autonomous underwater vehicle[J]. Applied ocean research, 2021, 106(1): 102441. DOI: 10.1016/J.APOR.2020.102441
[4] WANG C L, WANG B, DENG Z H, et al. A delaunay triangulation based matching area selection algorithm for underwater gravity-aided inertial navigation[J]. IEEE/ASME transactions on mechatronics, 2021, 26(2): 908-917. DOI: 10.1109/tmech.2020.3012499
[5] 黄玉龙, 张勇刚, 赵玉新. 自主水下航行器导航方法综述[J]. 水下无人系统学报, 2019, 27(3): 232-253. DOI: 10.11993/j.issn.2096-3920.2019.03.002 [6] 张胜军, 李建成, 孔祥雪. 基于Laplace方程的垂线偏差法反演全球海域重力异常[J]. 测绘学报, 2020, 49(4): 452-460. DOI: 10.11947/j.AGCS.2020.20190108 [7] 冯浩. 未来水下重力辅助惯性导航微系统中的理论、模型和仿真研究[D]. 北京: 北京邮电大学, 2004. [8] 罗孝宽, 郭绍雍. 应用地球物理教程—重力磁法[M]. 北京: 地质出版社, 1991. [9] HAN Y, WANG B, DENG Z, et al. A combined matching algorithm for underwater gravity aided navigation[J]. IEEE/ASME transactions on mechatronics, 2018, 23(1): 233-241. DOI: 10.1109/tmech.2017.2774296
[10] 张育林, 王兆魁, 刘红卫. 地球重力场天基测量理论及其内编队实现方法[M]. 北京: 科学出版社, 2018: 156-162. [11] 陈俊勇, 文双江, 程鹏飞. 中国大地测量科学发展的若干问题[J]. 地球科学进展, 2001(5): 681-688,9. DOI: 10.3321/j.issn:1001-8166.2001.05.012 [12] 陈秋杰, 沈云中, 张兴福, 等. 基于GRACE卫星数据的高精度全球静态重力场模型[J]. 测绘学报, 2016, 45(4): 396-403. DOI: 10.11947/j.AGCS.2016.20150422 [13] 丁剑, 许厚泽, 章传银. 基于大地水准面经典定义的地球重力场模型评价[J]. 中国科学院大学学报, 2016, 33(4): 528-536. DOI: 10.7523/j.issn.2095-6134.2016.04.014 [14] AFFLECK C A, JIRCITANO A. Passive gravity gradiometer navigation system[C]//IEEE Symposium on Position, Location and Navigation. Las Vegas, USA, 1990: 60-66. DOI: 10.1109/PLANS.1990.66158
[15] DOSCH D, JIRCITANO A. Gravity aided inertial navigation system (GAINS): US19910805544[P]. 1994-08-23.
[16] 李正心. 重力场变化与地震[M]. 北京: 科学出版社, 2020: 34-45. [17] 王凯, 张传定, 吴星, 等. 地形/均衡重力场模型构制及其拟稳分析[J]. 地球物理学报, 2014, 57(3): 760-769. DOI: 10.6038/cjg20140307 [18] 徐翰, 周强波. 卫星重力梯度数据重力异常的精度分析[J]. 测绘科学, 2016, 41(11): 17-24. DOI: 10.16251/j.cnki.1009-2307.2016.11.004 [19] 闫利, 吴华玲. 重力场基准图的多尺度分析研究[J]. 测绘通报, 2009(2): 10-13. [20] 徐遵义, 姜玉祥, 赵亮, 等. 改进的Shepard算法及其在重力异常插值中的应用[J]. 武汉大学学报(信息科学版), 2010, 35(4): 477-480. [21] 杨金玉, 张训华, 张菲菲, 等. 应用多种来源重力异常编制中国海陆及邻区空间重力异常图及重力场解读[J]. 地球物理学报, 2014, 57(12): 3920-3931. DOI: 10.6038/cjg20141206 [22] 聂琳娟, 超能芳, 金涛勇, 等. 卫星测高反演海洋重力异常的精度分析[J]. 测绘科学, 2016, 41(9): 1-6,37. DOI: 10.16251/j.cnki.1009-2307.2016.09.001 [23] 潘宗鹏, 王伟祖, 安然. 多源重力场信息融合处理方法[J]. 测绘科学技术学报, 2013, 30(3): 236-240. DOI: 10.3969/j.issn.1673-6338.2013.03.005 [24] 齐东旭. 分形及其计算机生成[M]. 北京: 科学出版社, 1994.
计量
- 文章访问数: 220
- HTML全文浏览量: 119
- PDF下载量: 39
