分形插值和网函数插值结合的卫星测高重力数据与船测重力数据融合方法

摘要

本发明公开了一种分形插值和网函数插值结合的卫星测高重力数据与船测重力数据融合方法，基本步骤为：利用网函数插值对卫星测高重力数据插值；计算插值结果与船测重力数据的差值；对差值进行聚类分析，将整体区域划分成多个子区域；在每个子区域各纵向和横向网格边上对差值进行分形插值，在子区域内部进行网函数插值，得到重力数据改正值；以重力数据改正值为中心，根据卫星测高重力值融合计算得到每个子区域数据融合产品，最终得到整个区域数据融合产品。本发明提供的方法科学合理，实现过程简单易懂，在对船测重力数据和卫星测高重力数据进行融合后能够得到高精度高分辨率的重力数据，同时充分保持了多源重力数据局部特征与整体特征相似性。

说明书

技术领域

本发明涉及数据融合技术领域，尤其涉及多源海洋重力数据融合的应用领域，具体说是一种分形插值和网函数插值结合的卫星测高重力数据与船测重力数据融合方法。

背景技术

高精度高分辨率的海洋重力场有助于研究地球内部地质构造，同时也是研究海洋动力环境变化的重要基础，此外，高精度高分辨率的重力异常对于勘探开发海洋矿产资源、保障航天器和远程武器的发射也具有重要意义。传统意义上的海洋重力场探测主要依靠海底布设或海面舰船搭载的重力测量仪器，其精度高但成本高昂且作业周期长，重力异常数据的分辨率较低，且由于各种因素的制约，无法获得敏感区域的重力异常数据。随着航天技术的发展，利用卫星测高反演海洋重力异常技术已经成熟，同时其也成为获取大范围高分辨率海洋重力异常的唯一手段，但其不足之处在于反演重力异常精度较低。在此情形下，如何对船测高精度数据和卫星测高反演的高分辨率数据进行融合成为一大焦点性问题。

目前对于两种数据的融合主要有统计法和解析法两种，统计法主要是基于统计学原理，应用最小二乘配置法进行数据融合，但这种方式需要大量观测数据构建协方差函数，且不同函数模型的选择没有明确标准，总体而言，统计法操作和计算都较为复杂，不易实现。解析法重在数据插值，但目前其插值方式只是依赖于数学中传统插值的一般形式，都是使插值曲线或曲面朝着越来越光滑的方向发展，其对于离散点、极不规则曲线或曲面而言有较大不足。由于许多自然界中的景物形态极不规则，传统意义上的插值难以对其进行刻画，而近年来发展起来的分形几何理论是对其进行描述的有力工具。有鉴于此，本文提出了一种分形插值与网函数插值结合的卫星测高重力数据与船测重力数据融合方法。

发明内容

(一)要解决的技术问题

提出一种分形插值和网函数插值结合的卫星测高重力数据与船测重力数据融合方法，充分考虑多源重力数据局部特征与整体的相似性，通过利用分形插值与网函数插值相结合和使用聚类分析划分数据融合子区域的方式，融合得到高精度高分辨率的重力数据。

(二)技术方案

本发明包含以下步骤：

(1)确定n个船测重力数据g_k(1≤k≤n)所在的卫星测高重力网格，利用网函数将网格角点的卫星测高重力数据G_i(1≤i≤4)插值到船测重力数据g_k(1≤k≤n)处，得到船测重力点的卫星测高重力插值结果G_k'(1≤k≤n)；

(2)计算卫星测高重力插值结果G_k'(1≤k≤n)与船测重力数据g_k(1≤k≤n)的差值ΔG_k(1≤k≤n)；

(3)对ΔG_k(1≤k≤n)进行聚类分析，将船测重力点区域划分成t个子区域，每个子区域有s条纵向网格边、r条横向网格边和s·r个船测重力点；

(4)在第p个子区域分别沿s条纵向网格边和r条横向网格边对ΔG_k(1≤k≤s·r)进行分形插值，得到第p个子区域中每条纵向和横向网格边上的重力数据改正值ΔGA_p'(1≤p≤t)；

(5)根据步骤(4)中得到的第p个子区域每条纵向和横向网格边上的重力数据改正值ΔGA_p'(1≤p≤t)，利用1-网函数插值计算第p个子区域内的重力数据改正值ΔGB_p'(1≤p≤t)，得到第p个子区域所有重力数据改正值ΔG_p'＝{ΔGA_p'，ΔGB_p'}；

(6)以第p个子区域重力数据改正值ΔG_p'(1≤p≤t)为中心，卫星测高重力网格分辨率的2～3倍为搜索半径，确定第p个子区域在搜索范围内的u个卫星测高重力数据>pi(1≤p≤t,1≤i≤u)；

(7)根据第p个子区域重力数据改正值ΔG_p'和步骤(6)中确定的卫星测高重力数据>pi(1≤p≤t,1≤i≤u)得到第p个子区域数据融合结果

(8)根据每个子区域中的数据融合结果得到t个子区域中所有的数据融合产品

进一步，所述步骤(1)中网函数插值类型为1-网函数插值，公式为

其中G₁(L₀,B₀)，G₂(L₁,B₀)，G₃(L₁,B₁)，G₄(L₀,B₁)是船测重力点g_k(1≤k≤n)所在卫星测高重力网格的四个数据点，船测重力点g_k(1≤k≤n)坐标为(L,B)，Q_i(i＝1,2,3,4)是过船测重力点g_k(1≤k≤n)与卫星测高重力网格平行的两条直线在网格边界上截得的四个点，>i(i＝1,2,3,4)是两条直线将网格分成四个小矩形的面积，A是网格总面积，>1+A₂+A₃+A₄，其中

A₁＝(L-L₀)(B-B₀)，A₂＝(L₁-L)(B-B₀)

A₃＝(L₁-L)(B₁-B)，A₄＝(L-L₀)(B₁-B)

进一步，所述步骤(2)中的计算公式为：

ΔG_k＝G_k'-g_k(1≤k≤n)

其中，g_k为船测重力数据，ΔG_k为卫星测高重力数据插值结果与船测重力数据的差值。

进一步，所述步骤(3)中划分子区域过程为

a.分别以所有船测重力点为中心，船测重力网格分辨率的2倍为搜索半径，找到搜索范围内count₁个船测重力点处的ΔG_k(1≤k≤count₁)，计算ΔG_k＞0的船测重力点数量count₂，当时认为该点为正值，当时，认为该点为负值，否则认为该点为正负交界的边界点；

b.针对所有非边界点，采用DBSCAN密度聚类算法对其进行聚类分析，得到t类船测重力点；

c.在t类船测重力点分布的区域中分别选择面积最大的矩形作为子区域，共得到t个子区域，每个子区域中有s条纵向网格边、r条横向网格边和s·r个船测重力点。

进一步，所述步骤(4)中分形插值在第p个子区域的每条纵向和横向网格边上进行，分形插值计算过程为

a.构造迭代函数系统{R²,ω_n,P_n,n＝1,2,…,N}如下

其中ω_n为压缩仿射变换，x为自变量，y为因变量，N为自变量个数，迭代函数系统中P_n满足：P_n＝1/N；

(1)式中常数a_n,c_n,d_n,e_n,f_n满足如下条件条件

即

其中{d_n}为纵向压缩因子，且满足0≤d_n＜1；

b.在第p个子区域中分别以经度L_i(1≤i≤s)和纬度B_j(1≤j≤r)为自变量，ΔG_k(1≤k≤s·r)为因变量，利用(1)式进行分形插值得到第p个子区域每条纵向和横向网格边的重力数据改正值ΔGA_p'(1≤p≤t)。

进一步，所述步骤(7)中计算公式为：

其中G_pi为第p个子区域中的卫星测高重力数据，ΔG_p'为第p个子区域中的重力数据改正值。

进一步，所述步骤(8)中计算公式为：

(三)有益效果

本方法以网函数插值和分形插值为基础，利用聚类算法从统计学角度划分不同融合子区域进行数据融合，在对船测重力数据和卫星测高重力数据进行融合时充分保持其局部特征与整体特征相似性，实现过程简单易懂，两种数据的融合能够得到高精度高分辨率的重力异常。

附图说明

图1为本发明实施的步骤流程图

图2为网格内1-网函数插值示意图

图3为船测重力点分为两类时选择面积最大矩形示意图

图4为子区域有3*3个船测重力点的分形插值示意图

具体实施方式

为使本发明的目的、内容和优点更加清楚，下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述：

参照图1，本发明的实施步骤为：

(1)确定n个船测重力数据g_k(1≤k≤n)所在的卫星测高重力网格，利用网函数将网格角点的卫星测高重力数据G_i(1≤i≤4)插值到船测重力数据g_k(1≤k≤n)处，得到船测重力点的卫星测高重力插值结果G_k'(1≤k≤n)，网函数插值类型为1-网函数插值，公式为：

1-网函数插值如图2所示，其中G₁(L₀,B₀)，G₂(L₁,B₀)，G₃(L₁,B₁)，G₄(L₀,B₁)是船测重力点g_k(1≤k≤n)所在卫星测高重力网格的四个数据点，船测重力点g_k(1≤k≤n)坐标为(L,B)，Q_i(i＝1,2,3,4)是过船测重力点g_k(1≤k≤n)与卫星测高重力网格平行的两条直线在网格边界上截得的四个点，A_i(i＝1,2,3,4)是两条直线将网格分成四个小矩形的面积，A是网格总面积，A＝A₁+A₂+A₃+A₄，其中

A₁＝(L-L₀)(B-B₀)，A₂＝(L₁-L)(B-B₀)

A₃＝(L₁-L)(B₁-B)，A₄＝(L-L₀)(B₁-B)

(2)计算卫星测高重力插值结果G_k'(1≤k≤n)与船测重力数据g_k(1≤k≤n)的差值ΔG_k(1≤k≤n)，计算公式为：

ΔG_k＝G_k'-g_k(1≤k≤n)

其中，g_k为船测重力数据，ΔG_k为卫星测高重力插值结果与船测重力数据的差值。

(3)对ΔG_k(1≤k≤n)进行聚类分析，将船测重力点区域划分成t子区域，每个子区域有s>

a.分别以所有船测重力点为中心，船测重力网格分辨率的2倍为搜索半径，找到搜索范围内count₁个船测重力点处的ΔG_k(1≤k≤count₁)，计算ΔG_k＞0的船测重力点数量count₂，当时认为该点为正值，当时，认为该点为负值，否则认为该点为正负交界的边界点；

b.针对所有非边界点，采用DBSCAN密度聚类算法对其进行聚类分析，得到t类船测重力点；

c.在t类船测重力点分布的区域中分别选择面积最大的矩形作为子区域，共得到t个子区域，每个子区域中有s条纵向网格边、r条横向网格边和s·r个船测重力点。当船测重力点分为两类时，选择面积最大矩形如图3所示。

(4)在第p个子区域分别沿s条纵向网格边和r条横向网格边对ΔG_k(1≤k≤s·r)进行分形插值，得到第p个子区域中每条纵向和横向网格边上的重力数据改正值ΔGA_p'(1≤p≤t)，分形插值在子区域的每条纵向和横向网格边上进行，分形插值计算过程为

a.构造迭代函数系统{R2,ω_n,P_n,n＝1,2,…,N}如下

其中ω_n为压缩仿射变换，x为自变量，y为因变量，N为自变量个数，迭代函数系统中P_n满足：P_n＝1/N；

(1)式中常数a_n,c_n,d_n,e_n,f_n满足如下条件条件

即

其中{d_n}为纵向压缩因子，且满足0≤d_n＜1；

b.在第p个子区域分别以经度L_i(1≤i≤s)和纬度B_j(1≤j≤r)为自变量，ΔG_k(1≤k≤s·r)>p'(1≤p≤t)。当子区域有3*3个船测重力点时，每条纵向和横向网格边上分形插值如图4所示。

(5)根据步骤(4)中得到的第p个子区域每条纵向和横向网格边上的重力数据改正值ΔGA_p'(1≤p≤t)，利用1-网函数插值计算第p个子区域内的重力数据改正值ΔGB_p'(1≤p≤t)，得到第p个子区域所有重力数据改正值ΔG_p'＝{ΔGA_p'，ΔGB_p'}；

(6)以第p个子区域重力数据改正值ΔG_p'(1≤p≤t)为中心，卫星测高重力网格分辨率的2～3倍为搜索半径，确定第p个子区域在搜索范围内的u个卫星测高重力数据>pi(1≤p≤t,1≤i≤u)；

(7)根据第p个子区域重力数据改正值ΔG_p'和步骤(6)中确定的卫星测高重力数据>pi(1≤p≤t,1≤i≤u)得到第p个子区域数据融合结果计算公式为：

其中G_pi为第p个子区域中的卫星测高重力数据，ΔG_p'为第p个子区域中的重力数据改正值。

(8)根据每个子区域中的数据融合结果得到t个子区域中所有的数据融合产品，公式为：

。