CN109085656B - 一种面向特征的高精度重力图构建与插值方法 - Google Patents

Abstract

本发明涉及一种面向特征的高精度重力图构建与插值方法，其技术特点在于：包括以下步骤：步骤1、根据重力场特征分布特点，构建包含精准特征信息的高精度重力图；步骤2、根据步骤1构建的高精度重力图，利用多次线性插值方法或径向重力异常非线性特征等效法求取区域内任意点P(r_P,θ_p)处的重力异常值。本发明采用直接描述重力场特征方式，包含区域重力异常峰值点与从峰值点起沿径向重力异常下降的非线性信息，对重力场特征的描述更高效更准确，同时减少了构建重力图的船载测量工作量，满足重力场特征定位技术对重力场特征信息的需求。

Description

一种面向特征的高精度重力图构建与插值方法

技术领域

本发明属于地球物理场导航定位技术领域，涉及地球重力场信息资源在自主导航定位技术中的应用，尤其是一种面向特征的高精度重力图构建与插值方法。

背景技术

地球重力场是反映地球物质分布特征的物理场，实际重力场与正常重力场之差称为扰动重力场，扰动重力场包括重力异常、垂线偏差等参量。重力场信息具有长期稳定性和独特的空间分布特性，与当地地理位置密切相关，给导航提供了天然的参考基准。重力场信息定位技术利用重力场时空分布特征与地理位置信息的相关性，匹配重力实时测量值与重力图值，确定载体精确位置。重力场定位信息可用于惯性导航系统误差校正，解决目前存在的水下作战或作业平台导航系统缺乏水下无源校正手段的技术难题，大幅提升深远海长航时高精度自主导航能力。重力场定位技术需应用于重力场特征明显区域，重力图是实现重力场定位的信息基础，需要精准反映重力场特征分布。

重力测量数据是构建重力图的基础，现有海洋重力场资料主要通过卫星重力、卫星测高反演海洋重力、船载海洋重力测量等技术手段获取。卫星测高数据覆盖全球海区，但精度和分辨率偏低。船载海洋重力测量具有高精度、高分辨率特点，但存下述在问题：一是测量成本高、周期长，数据覆盖范围有限；二是目前船载重力测量在测区内等间距沿经度方向布设测线，同时按一定比例在测线垂线方向配置检测线，如图1所示，其测线布置与重力场特征分布无关，在反映重力场特征方面效率和能力不足。

发明内容

本发明的目的在于克服现有技术中的不足，提供一种设计合理、高效精确且能够减少构建高精度重力图所需船载测量工作量的精准捕捉重力场特征信息的高精度重力图构建与插值方法。

本发明解决其现实问题是采取以下技术方案实现的：

一种面向特征的高精度重力图构建与插值方法，包括以下步骤：

步骤1、根据重力场特征分布特点，构建包含精准特征信息的高精度重力图；

步骤2、根据步骤1构建的高精度重力图，利用多次线性插值方法或径向重力异常非线性特征等效法求取区域内任意点P(r_P,θ_p)处的重力异常值。

而且，所述步骤1的具体步骤包括：

(1)基于卫星测高或其它途径得到的重力图信息，初步确定区域重力异常峰值点，记为O’；

(2)基于卫星测高或其它途径得到的重力图信息，从O’点出发向四周延伸，0°到360°范围内搜索并确定重力异常下降最快方向n1、次快方向n2、下降最慢方向m1与次慢方向m2共4个极值方向；

(3)设计一组测线，获得该组测线上海洋重力测量数据，确定精确重力异常峰值点O位置，并完成O点精确性检验；

(4)建立以O点为极点的极坐标系，以从O点出发指向东向的射线为极轴Ox；

(5)设计一组圆周多边形测线，测线间距由重力异常沿径向下降的线性度决定，获取圆周多边形测线上海洋重力测量数据；

(6)设计一组径向测线，获取径向测线上海洋重力测量数据；

(7)圆周多边形测线和径向测线的海洋重力异常测量数据包含了这个区域重力场分布信息与特征信息，构成了这个区域的高精度重力图。

而且，所述步骤2的具体步骤包括：

(1)当P点位于两条相邻圆周多线段测线之间时，使用线性插值法计算点P处的重力异常值G_p；

(2)当P点位于O点与圆周多边形测线之间时，基于径向测线上重力异常值计算点P处的重力异常值G_p。

而且，所述步骤1第(3)步的具体步骤包括：

①设计测线A，沿n1方向，中点为O’，实施海洋重力测量得到测线A上重力异常测量值，记录测线A上重力异常峰值点O1；

②设计测线B，沿n1的垂线方向n1_⊥，中点为O1，实施海洋重力测量得到测线B上重力异常测量值，记录测线B上重力异常峰值点即是区域重力异常峰值点O；

③重力异常峰值点O的精确性检验：设计测线E，方向任取，中点为O；实施海洋重力测量得到测线E上重力异常测量值，测线E上重力异常峰值处与点O距离，可以反映区域重力异常峰值点O的精确程度。

而且，所述步骤2第(1)步的具体步骤包括：

①求过P点的极径方向与P点相邻2条圆周测线交点F和E处的重力异常值G_F和G_E；

首先按下式计算则线段SF长度L_SF：

L_SF＝L_ST*sinα/(sinα+sinβ)

其中，S和T是F点所在圆周多边形边上的顶点，L_ST为线段ST长度，α为线段OP与OS夹角，β为线段OP与OT夹角；

根据L_SF，在圆周测线ST上重力异常值中插值得到F点重力异常值G_F。采用与G_F同样方法得到E点处重力异常值G_E。

②F点极径为r_F，E点极径为r_E，按下式计算P点的重力异常值G_P

G_P＝G_E*(r_F-r_P)/(r_F-r_E)+G_F*(r_P-r_E)/(r_F-r_E)

而且，所述步骤2第(2)步的具体步骤包括：

①过P点的极径方向与P点相邻圆周测线交于点E，按前述方法得到E点处重力异常值G_E。

②以m1与m2为界将区域分为两部分，若P点与n1位于同一部分，令方向n＝n1；若P点与n2位于同一部分，令方向n＝n2。

③获取过P点极径方向的位置序列r_j及其对应的重力异常序列G_j

首先按下式计算位置序列r_j：

r_j＝r_i*r_E/r_E′

其中，r_i为重力图沿n方向径向测线上重力异常序列G_i对应的位置序列，r_E为E点极径，r_E′为沿n方向径向测线上重力异常值为G_E的点E′的极径；由于径向测线上重力异常具有随极径增大而单调递减特性，极径r_E′可根据其对应重力异常值G_E在沿n方向径向测线上重力异常序列G_i及对应位置序列r_i中线性插值得到。

由于区域重力异常峰值O点附近的重力异常等值线形状具有高度相似性，因此将沿n方向径向测线上位置序列r_i对应的重力异常序列G_i可视为过P点极径方向上位置序列r_j对应的重力异常序列G_j，即G_j＝G_i。

④根据P点极径r_P，在过P点极径方向位置序列r_j和对应重力异常序列G_j中线性插值得到P点处重力异常值G_P。

本发明的优点和有益效果：

1、本发明提出的面向特征的高精度重力图构建与插值方法，直接搜索确定特征区域重力异常峰值信息，精度高，弥补了传统重力图构建方法中易于漏掉或难以精确捕捉重力异常峰值点这一重要特征信息的缺点，为重力场定位提供更有力支撑。

2.本发明建立了以重力异常峰值点为原点的极坐标系，采用圆周和径向测线方式，插值方法不需要对重力测量数据进行网格化处理，充分保留了重力场非线性信息，同时可根据重力异常随径向距离的线性度设计圆周测线间距，设计更为科学合理，在减少构建重力图所需船载测量工作量同时，可更为全面和精确的描述重力场特征，从而提高重力场特征定位的精度。

3.本发明提出的面向特征的高精度重力图构建与插值方法，根据目标点在特征区所处位置不同，采用多次线性插值或径向重力异常非线性特征等效法求取目标点重力异常值，与常用平面插值方法相比，减少了计算量，保留了非线性信息，提高了精度。

附图说明

图1是本发明的船测重力测线示意图；

图2是本发明的特征区重力场示意图；

图3是本发明的南海某区域重力场分布情况示意图；

图4是本发明的L1-L12测线分布示意图；

图5是本发明的L1-L6测线上重力异常图；

图6是本发明的L7-12测线上重力异常图；

图7是本发明的N1-N12测线分布示意图；

图8是本发明的N1-N12测线上重力异常图；

图9是本发明的重力特征及极坐标系示意图；

图10是本发明的半径为3海里处圆周上重力异常曲线图；

图11是本发明的确定重力场峰值O点的航迹示意图；

图12是本发明的基于极坐标的圆周多边形测线设计示意图；

图13是本发明的区域内任一点P位置示意图1；

图14是本发明的区域内任一点P位置示意图2；

图15是本发明的重力异常峰值点附近形状高度相似的重力等值线示意图；

图16是本发明的“基于现有重力图方法”和“基于特征重力图方法”构建的重力图及线段L1、L2、L3示意图；

图17是本发明的“基于现有重力图方法”和“基于特征重力图方法”得到的线段L1重力异常图；

图18是本发明的“基于现有重力图方法”和“基于特征重力图方法”得到的线段L1重力异常误差图；

图19是本发明的“基于现有重力图方法”和“基于特征重力图方法”得到的线段L2重力异常图；

图20是本发明的“基于现有重力图方法”和“基于特征重力图方法”得到的线段L2重力异常误差图；

图21是本发明的“基于现有重力图方法”和“基于特征重力图方法”得到的线段L3重力异常图；

图22是本发明的“基于现有重力图方法”和“基于特征重力图方法”得到的线段L3重力异常误差图；

具体实施方式

以下结合附图对本发明实施例作进一步详述：

一种面向特征的高精度重力图构建与插值方法，包括以下步骤：

步骤1、根据重力场特征分布特点，构建包含精准特征信息的高精度重力图；

所述步骤1的具体步骤包括：

(1)基于卫星测高或其它途径得到的重力图信息，初步确定区域重力异常峰值点，记为O’(Lono1，Lato1)。

(2)基于卫星测高或其它途径得到的重力图信息，从O’点出发向四周延伸，0°到360°范围内搜索并确定重力异常下降最快方向n1、次快方向n2、下降最慢方向m1与次慢方向m2共4个极值方向。

以图9所示区域重力场分布特征为例，距其重力异常峰值点O’处3海里的圆周上重力异常值如图10所示，从图10可得到重力异常下降最快方向n1、次快方向n2与下降最慢方向m1与次慢方向m2共4个极值方向，这4方向对应的角度分别为233°、56°、145.8°、325.2°。

(3)设计一组测线，获得该组测线上海洋重力测量数据，确定精确重力异常峰值点O位置，并完成O点精确性检验。具体步骤如下：

①设计测线A，沿n1方向，中点为O’，测线A长度建议取4-6海里。实施海洋重力测量得到测线A上重力异常测量值，记录测线A上重力异常峰值点O1。

②设计测线B，沿n1的垂线方向n1_⊥，中点为O1，如图11所示。实施海洋重力测量得到测线B上重力异常测量值，记录测线B上重力异常峰值点即是区域重力异常峰值点O。

③重力异常峰值点O的精确性检验：设计测线E，方向任取，中点为O，如图11所示；实施海洋重力测量得到测线E上重力异常测量值，测线E上重力异常峰值处与点O距离，可以反映区域重力异常峰值点O的精确程度。

(4)建立以O点为极点的极坐标系，以从O点出发指向东向的射线为极轴Ox。

(5)设计一组圆周多边形测线，测线间距由重力异常沿径向下降的线性度决定，获取圆周多边形测线上海洋重力测量数据。。

如图12所示，设计一组以O点为圆心不同半径r的圆周多边形测线，当r较小时，圆周多边形可取四边形，其顶点为半径为r的圆周与n1、n2、m1、m2方向交点；随着半径r增加，圆周多边形测线可取八边形，其顶点为半径为r的圆周与n1、n2、m1、m2方向以及n1、n2、m1、m2中各相邻两个方向的角平分线方向交点。

完成圆周多边形测线海洋重力测量，获得测线上重力异常测量数据。

(6)设计一组径向测线，获取径向测线上海洋重力测量数据。

设计2条径向测线，以O点为起点，分别沿n1、n2方向，长度由重力场特征的覆盖范围而定。

完成海洋重力测量，获得测线上重力异常测量数据。

(7)圆周多边形测线和径向测线的海洋重力异常测量数据包含了这个区域重力场分布信息与特征信息，构成这个区域的高精度重力图。

步骤2、根据步骤1构建的高精度重力图，利用多次线性插值方法或径向重力异常非线性特征等效法求取区域内任意点P(r_P,θ_p)处的重力异常值。

所述步骤2的具体步骤包括：

(1)当P点位于两条相邻圆周多线段测线之间时，如图13所示，使用线性插值法计算点P处的重力异常值G_p。

①求过P点的极径方向与P点相邻2条圆周测线交点F和E处的重力异常值G_F和G_E；

首先按下式计算则线段SF长度L_SF：

L_SF＝L_ST*sinα/(sinα+sinβ)

其中S和T是F点所在圆周多边形边上的顶点，L_ST为线段ST长度，α为线段OP与OS夹角，β为线段OP与OT夹角；

根据L_SF，在圆周测线ST上重力异常值中插值得到F点重力异常值G_F。采用与G_F同样方法得到E点处重力异常值G_E。

②F点极径为r_F，E点极径为r_E，按下式计算P点的重力异常值G_P

G_P＝G_E*(r_F-r_P)/(r_F-r_E)+G_F*(r_P-r_E)/(r_F-r_E)

(2)如果P点位于O点与圆周多边形测线之间，如图14所示，则基于径向测线上重力异常值按下述步骤计算点P处的重力异常值G_p。

①过P点的极径方向与P点相邻圆周测线交于点E，按前述方法得到E点处重力异常值G_E。

②以m1与m2为界将区域分为两部分，若P点与n1位于同一部分，令方向n＝n1；若P点与n2位于同一部分，令方向n＝n2。如图14所示情况n＝n2。

③获取过P点极径方向的位置序列r_j及其对应的重力异常序列G_j

首先按下式计算位置序列r_j：

r_j＝r_i*r_E/r_E′

其中，r_i为重力图沿n方向径向测线上重力异常序列G_i对应的位置序列，r_E为E点极径，r_E′为沿n方向径向测线上重力异常值为G_E的点E′的极径；由于径向测线上重力异常具有随极径增大而单调递减特性，极径r_E′可根据其对应重力异常值G_E在沿n方向径向测线上重力异常序列G_i及对应位置序列r_i中线性插值得到。

由于区域重力异常峰值O点附近的重力异常等值线形状具有高度相似性，如图15所示，因此将沿n方向径向测线上位置序列r_i对应的重力异常序列G_i可视为过P点极径方向上位置序列r_j对应的重力异常序列G_j，即G_j＝G_i。

④根据P点极径r_P，在过P点极径方向位置序列r_j和对应重力异常序列G_j中线性插值得到P点处重力异常值G_P。

在本实施例中，为验证本发明的优越性，进行了面向特征的重力图构建与插值方法效果的仿真验证：

对于图9所示区域，其一是利用现有高精度重力图构建方法构建重力图，再采用常用平面双线性插值方法获取区域内一组线段的重力异常值(简称为“基于现有重力图方法”)；其二是利用本发明提出的面向特征的高精度重力图构建方法构建重力图，再采用本发明提出的插值方法获取区域内一组线段的重力异常值(简称为“基于特征重力图方法”)。比较上述两种方法得到的同一组线段重力异常值的误差，验证本发明提出的面向特征的重力图构建与插值方法的效果。

基于现有重力图方法和基于特征重力图方法构建的高精度重力图如图16所示，现有重力图方法采用间距为1海里的网格数据，基于特征重力图方法包含极径为3海里、5海里与7海里的圆周多边形测线上重力异常测量序列和沿n1、n2方向径向测线上重力异常测量序列(可取m1、m2方向径向测线为检测线，则重力图中亦包含沿m1、m2方向径向测线上重力异常测量序列)。基于上述方法分别求取图16中所示线段L1、L2、L3上重力异常值，结果如图17和图18、图19和图20、图21和图22所示，从图中可知本发明提出的面向特征的高精度重力图构建和插值方法得到的L1、L2、L3这3个线段上重力异常误差明显小于现有高精度重力图构建与插值方法结果，验证了本发明提出的方法的先进性。

本发明的工作原理是：

图2所示区域的重力场分布呈现“山峰”特征，此类特征区域在各大洋中较为常见，图3为南海某区域重力图，包含了若干具有“山峰”特征的区域。相对于重力变化平缓的区域，此类区域重力变化剧烈，特征信息含量大，适于重力场导航定位。本发明提出了面向此类特征的高精度重力图构建与插值方法。对于此类重力场特征，如图1所示的传统测线布设方式构建的重力图难以精确捕捉重力场主要特征，因为多条测线上可能存在重力异常峰值，如图4所示的测线L1-L12，其上重力异常如图5、图6所示，重力异常曲线存在峰值，重力图网格数据极有可能丢失这些峰值信息，提高网格分辨率可以缓解这一问题，但将增加测量工作量且对重力变化线性度较好区域没有改善作用。

图2所示重力特征区，重力异常峰值点为O，从O点出发沿0°到360°方向均随着远离O点而单调下降，O点附近区域重力异常下降的非线性明显，随着远离O点，重力异常随距离下降的线性度变好。如图7所示，以O点为起点，向四周延伸得到测线N1至N12，测线N1至N12上重力异常曲线如图8所示，重力异常呈现单调下降特性，且随着远离O点，重力异常下降线性度变好，与如5、图6比较，图7方式更适合描述特征区重力场信息。

本发明提出的面向特征的高精度重力图构建方法，首先根据卫星测高数据得到的重力图对区域内重力特征进行分析，确定重力异常峰值点O′以及从O′出发沿各方向线段中重力异常随距离下降最快、最慢的方向；在O′基础上通过两个正交测线的重力异常测量搜索确定精确的重力异常峰值点O；建立以区域峰值点O为极点的极坐标系，根据特征分布设计圆周多边形测线和径向测线，各圆周多边形测线与径向测线的重力异常测量序列构成该特征区域高精度重力图，不需对重力异常测量数据进行网格化处理。

本发明提出了适于面向特征的高精度重力图的插值方法。假定求取任意点P处重力异常值，当P点距离O点较远时，径向重力异常随距离下将线性度较好，该方法首先通过圆周多边形测线重力异常测量序列插值获取过P点的极径方向与圆周多边形测线各交点的重力异常值，再根据各交点重力异常值插值得到P点的重力异常值。当P点距离O点较近时，径向重力异常随距离下降的非线性较为明显，此时利用O点附近区域等值线形状高度相似性，将重力图径向测线的重力异常测量序列等效到对P点极径方向，再对过P点极径方向重力异常序列插值得到P点处的重力异常值。

需要强调的是，本发明所述实施例是说明性的，而不是限定性的，因此本发明包括并不限于具体实施方式中所述实施例，凡是由本领域技术人员根据本发明的技术方案得出的其他实施方式，同样属于本发明保护的范围。

Claims (2)

1.一种面向特征的高精度重力图构建与插值方法，其特征在于：包括以下步骤：

步骤1、根据重力场特征分布特点，构建包含精准特征信息的高精度重力图；

步骤2、根据步骤1构建的高精度重力图，利用多次线性插值方法或径向重力异常非线性特征等效法求取区域内任意点P(r_P,θ_p)处的重力异常值；

所述步骤1的具体步骤包括：

(1)基于卫星测高或其它途径得到的重力图信息，初步确定区域重力异常峰值点，记为O’；

(2)基于卫星测高或其它途径得到的重力图信息，从O’点出发向四周延伸，0°到360°范围内搜索并确定重力异常下降最快方向n1、次快方向n2、下降最慢方向m1与次慢方向m2共4个极值方向；

(3)设计一组测线，获得该组测线上海洋重力测量数据，确定精确重力异常峰值点O位置，并完成O点精确性检验；

(4)建立以O点为极点的极坐标系，以从O点出发指向东向的射线为极轴Ox；

(5)设计一组圆周多边形测线，测线间距由重力异常沿径向下降的线性度决定，获取圆周多边形测线上海洋重力测量数据；

(6)设计一组径向测线，获取径向测线上海洋重力测量数据；

(7)圆周多边形测线和径向测线的海洋重力异常测量数据包含了这个区域重力场分布信息与特征信息，构成了这个区域的高精度重力图；

所述步骤2的具体步骤包括：

(1)当P点位于两条相邻圆周多边形测线之间时，使用线性插值法计算点P处的重力异常值G_p；

(2)当P点位于O点与圆周多边形测线之间时，基于径向测线上重力异常值计算点P处的重力异常值G_p；

所述步骤2第(1)步的具体步骤包括：

①求过P点的极径方向与P点相邻2条圆周测线交点F和E处的重力异常值G_F和G_E；

首先按下式计算线段SF长度L_SF：

L_SF＝L_ST*sinα/(sinα+sinβ)

其中，S和T是F点所在圆周多边形边上的顶点，L_ST为线段ST长度，α为线段OP与OS夹角，β为线段OP与OT夹角；

根据L_SF，在圆周测线ST上重力异常值中插值得到F点重力异常值G_F；采用与G_F同样方法得到E点处重力异常值G_E；

②F点极径为r_F，E点极径为r_E，按下式计算P点的重力异常值G_P

G_P＝G_E*(r_F-r_P)/(r_F-r_E)+G_F*(r_P-r_E)/(r_F-r_E)；

所述步骤2第(2)步的具体步骤包括：

①过P点的极径方向与P点相邻圆周测线交于点E，按前述方法得到E点处重力异常值G_E；

②以m1与m2为界将区域分为两部分，若P点与n1位于同一部分，令方向n＝n1；若P点与n2位于同一部分，令方向n＝n2；

③获取过P点极径方向的位置序列r_j及其对应的重力异常序列G_j

首先按下式计算位置序列r_j：

r_j＝r_i*r_E/r_E′

其中，r_i为重力图沿n方向径向测线上重力异常序列G_i对应的位置序列，r_E为E点极径，r_E′为沿n方向径向测线上重力异常值为G_E的点E′的极径；由于径向测线上重力异常具有随极径增大而单调递减特性，极径r_E′可根据其对应重力异常值G_E在沿n方向径向测线上重力异常序列G_i及对应位置序列r_i中线性插值得到；

由于区域重力异常峰值O点附近的重力异常等值线形状具有高度相似性，因此将沿n方向径向测线上位置序列r_i对应的重力异常序列G_i可视为过P点极径方向上位置序列r_j对应的重力异常序列G_j，即G_j＝G_i；

④根据P点极径r_P，在过P点极径方向位置序列r_j和对应重力异常序列G_j中线性插值得到P点处重力异常值G_P。

2.根据权利要求1所述的一种面向特征的高精度重力图构建与插值方法，其特征在于：所述步骤1第(3)步的具体步骤包括：

①设计测线A，沿n1方向，中点为O’，实施海洋重力测量得到测线A上重力异常测量值，记录测线A上重力异常峰值点O1；

②设计测线B，沿n1的垂线方向n1_⊥，中点为O1，实施海洋重力测量得到测线B上重力异常测量值，记录测线B上重力异常峰值点即是区域重力异常峰值点O；

③重力异常峰值点O的精确性检验：设计测线E，方向任取，中点为O；实施海洋重力测量得到测线E上重力异常测量值，测线E上重力异常峰值处与点O距离，可以反映区域重力异常峰值点O的精确程度。

Images