CN104933291B - 基于网函数插值的卫星测高数据平均海面高产品制作方法 - Google Patents

Abstract

本发明公开了一种基于网函数插值的卫星测高数据平均海面高产品制作方法,具体步骤分为：数据预处理、网格大小及中心点的确定、网格化插值、精细插值。本发明便于计算机自动划分网格，具备明显的统计特征，在物理特性上能够反映海洋重力场的空间分布规律，而海洋重力场与平均海面高成正相关，充分考虑了海平面的物理特性，能够提高平均海面高产品的精度和空间分辨率。

Description

基于网函数插值的卫星测高数据平均海面高产品制作方法

技术领域

本发明涉及一种制作方法，尤其涉及一种基于网函数插值的卫星测高数据平均海面高产品制作方法。

背景技术

平均海面高是指某段时间内的水位的平均值，时间尺度要足够长，消除了所有的物理振动后稳定的平均海面。

平均海面高产品是指利用多源卫星测高数据处理得到的格网化的平均海面高数据集。

高精度高分辨率的平均海面高数据是研究海洋动力环境变化的重要基础，是研究全球海平面变化和气候变化的重要支撑数据，有利于精确研究预测平均海面的变化趋势，进而作出地球物理和环境变化的解释。

随着测高卫星的不断发射和数据的累积获取，卫星测高技术已成为全球平均海面高数据获取的主要手段。卫星测高数据本是离散点集，相邻轨道间距较大，而测高点之间距离小，另外，不同测高卫星轨道不同，测高点分布无序错乱。为了获得格网化的平均海面高数据，需要对离散测高数据进行网格化处理(即数据插值问题)，插值方法的选择和使用直接影响平均海面高数据的制作精度。因此，需考虑合适的网格化插值算法。

目前格网化的算法很多，如距离加权最小二乘法、曲面拟合法、Shepard算法及谐核函数法等，这些都是数学中的一般方法，并没有考虑海平面的物理特性，而只有与插值对象保持物理特性较好的函数(例如周期函数适用于波动现象)方能获得高精度的插值结果，因此这些算法在平均海面高产品的制作精度上还存在缺陷。

发明内容

为了解决上述技术所存在的不足之处，本发明提供了一种基于网函数插值的卫星测高数据平均海面高产品制作方法。

为了解决以上技术问题，本发明采用的技术方案是：一种基于网函数插值的卫星测高数据平均海面高产品制作方法，其具体步骤如下：

a、数据预处理

根据各测高卫星的数据手册，开展格式转换与数据编辑，进而开展共线处理、交叉点平差，形成预处理后的测高点；

b、网格大小及中心点的确定

网格大小决定数据采样的空间分辨率，取值越大对于提高计算速度是有帮助的，但会丢失大量的细节信息，对于多源卫星测高数据,由于各轨道分布的非规律性，为保证插值精度，取轨道间距最稀疏处的一半为插值网格宽度，用W表示，指定左上角点坐标(x_t，y_t)为第1个网格点，结合网格宽度其它网格中心点坐标计算方法如下：G_(ij)＝(x_t+(i-1)×W，y_t+(j-1)×W)，其中i、j表示从左上角起算第i列、第j行；

c、网格化插值

预处理后的卫星测高点是离散、不规则的，采用三角形区域网函数插值算法，首先要构建以网格点为三角形内部点的最小外接三角形；根据最邻近点算法，在2倍网格宽度即2W为搜索半径找出距离网格点最近的3个数据点，判断网格点是否落在这三个数据点组成的三角形内部；若判断失败，则继续搜寻第4个最邻近点，然后利用这4个点中任意组合三角形，判断网格点是否落在三角形内部，若网格点落在2个以上三角形内部，则计算网格点与三角形顶点的总距离，取距离最短的三角形；若再次判断失败，则继续寻找下1个最邻近点，重复上述步骤，若在搜索半径内，无法构成该三角形，则该网格点失败；

判断网格点是否在三角形内可以使用向量同向法：若点P在三角形ABC内，则三个向量AB×AP、AP×AC、PB×PC平行同向(它们也与向量AB×AC平行同向)，由于这三个向量均有可能为0，直接判断它们平行同向比较麻烦，但考虑到AB×AC不可能为0，直接判断向量：AB×AP、AP×AC、PB×PC均与AB×AC平行同向更为简单；

三角形确定后，利用三角形区域网函数插公式求取网格点的值；其中Q_i是直线M_iQ和M_jM_k之交点，i,j,k∈{1,2,3}，F(Q)代表插值后的数值，Δ表示为三角形的总面积，Δ_1、Δ_2、Δ₃分别代表三角形M₃QM₂、M₃QM₁、M₁QM₂的面积；

然后，依次搜寻下一个网格点的最小外接三角形，并进行插值运算，方法同上，直到所有网格点都遍历完毕，形成规则的海面高数据点阵；

d、精细插值

基于网格化的卫星测高数据,利用矩形区域网函数插值算法内插出网格内部某点的海面高度,进一步提高海面高数据的空间分辨率，矩形面积加权网函数插值算法如下：

P₁(x₀,y₀)，P₂(x₁,y₀)，P₃(x₁,y₁)，P₄(x₀,y₁)是矩形区域D的四个角点(即相邻的四个网格数据点)，Q_i(i＝1,2,3,4)是过插值点Q与矩形平行的两条直线在矩形边界上截得的四个点；A_i(i＝1,2,3,4)是这两条直线将矩形区域D分成四个小矩形的面积，其中A₁＝(x-x₀)(y-y₀),A₂＝(x₁-x)(y-y₀),A₃＝(x₁-x)(y₁-y)，A₄＝(x–x₀)(y₁-y)；A＝(x₁–x₀)(y₁-y₀)＝A₁+A₂+A₃+A₄，是矩形区域D的总面积；

则Q点的海面高度值

以插值n倍为例，插值点坐标的计算方法为：x＝x₀+W/n，y＝y₀+W/n，其中W为网格宽度；

确定n后，逐行遍历每个网格点，针对每个网格分别计算插值点Q(x,y)的海面高度值，直至所有网格点遍历结束，生成新的平均海面高数据集。

本发明基于网函数插值理论推导严密准确、形式多样的特点，便于计算机自动划分网格，具备明显的统计特征，在物理特性上能够反映海洋重力场的空间分布规律，而海洋重力场与平均海面高成正相关，充分考虑了海平面的物理特性，能够提高平均海面高产品的精度和空间分辨率。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1为本发明的流程图。

图2为向量同向法示意图。

图3为三角形区域网函数插值算法示意图。

图4为矩形区域网函数插值算法图形示意。

具体实施方式

本发明的制作流程图如图1所示，针对卫星测高数据的离散特点，结合卫星轨道间距分析确定合理的网格大小和中心点位置，开展基于三角形区域网函数插值算法进行网格化插值，形成排列规则的海面高数据点阵，在此基础上应用矩形区域网函数插值算法开展精细插值，制作更高空间分辨率的平均海平面数据集。

具体制作方法步骤如下：

a、数据预处理

根据各测高卫星的数据手册，开展格式转换与数据编辑，进而开展共线处理、交叉点平差，形成预处理后的测高点；

b、网格大小及中心点的确定

网格大小决定数据采样的空间分辨率，取值越大对于提高计算速度是有帮助的，但会丢失大量的细节信息，对于多源卫星测高数据,由于各轨道分布的非规律性，为保证插值精度，取轨道间距最稀疏处的一半为插值网格宽度，用W表示，指定左上角点坐标(x_t，y_t)为第1个网格点，结合网格宽度其它网格中心点坐标计算方法如下：G_(ij)＝(x_t+(i-1)×W，y_t+(j-1)×W)，其中i、j表示从左上角起算第i列、第j行；

c、网格化插值

预处理后的卫星测高点是离散、不规则的，采用三角形区域网函数插值算法，首先要构建以网格点为三角形内部点的最小外接三角形；根据最邻近点算法，在2倍网格宽度即2W为搜索半径找出距离网格点最近的3个数据点，判断网格点是否落在这三个数据点组成的三角形内部；若判断失败，则继续搜寻第4个最邻近点，然后利用这4个点中任意组合三角形，判断网格点是否落在三角形内部，若网格点落在2个以上三角形内部，则计算网格点与三角形顶点的总距离，取距离最短的三角形；若再次判断失败，则继续寻找下1个最邻近点，重复上述步骤，若在搜索半径内，无法构成该三角形，则该网格点失败；

判断网格点是否在三角形内可以使用向量同向法：若点P在三角形ABC内，则三个向量AB×AP、AP×AC、PB×PC平行同向(它们也与向量AB×AC平行同向)，由于这三个向量均有可能为0，直接判断它们平行同向比较麻烦，但考虑到AB×AC不可能为0，直接判断向量：AB×AP、AP×AC、PB×PC均与AB×AC平行同向更为简单，如图2所示。

三角形确定后，利用三角形区域网函数插公式求取网格点的值；其中Q_i是直线M_iQ和M_jM_k之交点，i,j,k∈{1,2,3}，F(Q)代表插值后的数值，Δ表示为三角形的总面积，Δ_1、Δ_2、Δ₃分别代表三角形M₃QM₂、M₃QM₁、M₁QM₂的面积，如图3所示。

然后，依次搜寻下一个网格点的最小外接三角形，并进行插值运算，方法同上，直到所有网格点都遍历完毕，形成规则的海面高数据点阵；

d、精细插值

基于网格化的卫星测高数据,利用矩形区域网函数插值算法内插出网格内部某点的海面高度,进一步提高海面高数据的空间分辨率，矩形面积加权网函数插值算法如下：

P₁(x₀,y₀)，P₂(x₁,y₀)，P₃(x₁,y₁)，P₄(x₀,y₁)是矩形区域D的四个角点(即相邻的四个网格数据点)，Q_i(i＝1,2,3,4)是过插值点Q与矩形平行的两条直线在矩形边界上截得的四个点；A_i(i＝1,2,3,4)是这两条直线将矩形区域D分成四个小矩形的面积，其中A₁＝(x-x₀)(y-y₀),A₂＝(x₁-x)(y-y₀),A₃＝(x₁-x)(y₁-y)，A₄＝(x–x₀)(y₁-y)；A＝(x₁–x₀)(y₁-y₀)＝A₁+A₂+A₃+A₄，是矩形区域D的总面积，如图4所示。

则Q点的海面高度值

以插值n倍为例，插值点坐标的计算方法为：x＝x₀+W/n，y＝y₀+W/n，其中W为网格宽度；

确定n后，逐行遍历每个网格点，针对每个网格分别计算插值点Q(x,y)的海面高度值，直至所有网格点遍历结束，生成新的平均海面高数据集。

Claims (1)

1.一种基于网函数插值的卫星测高数据平均海面高产品制作方法，其特征在于：其具体步骤如下：

a、数据预处理

根据各测高卫星的数据手册，开展格式转换与数据编辑，进而开展共线处理、交叉点平差，形成预处理后的测高点；

b、网格大小及中心点的确定

网格大小决定数据采样的空间分辨率，取值越大对于提高计算速度是有帮助的，但会丢失大量的细节信息，对于多源卫星测高数据,由于各轨道分布的非规律性，为保证插值精度，取轨道间距最稀疏处的一半为插值网格宽度，用W表示，指定左上角点坐标(x_t，y_t)为第1个网格点，结合网格宽度其它网格中心点坐标计算方法如下：G_(ij)＝(x_t+(i-1)×W，y_t+(j-1)×W)，其中i、j表示从左上角起算第i列、第j行；

c、网格化插值

预处理后的卫星测高点是离散、不规则的，采用三角形区域网函数插值算法，首先要构建以网格点为三角形内部点的最小外接三角形；根据最邻近点算法，在2倍网格宽度即2W为搜索半径找出距离网格点最近的3个数据点，判断网格点是否落在这三个数据点组成的三角形内部；若判断失败，则继续搜寻第4个最邻近点，然后利用这4个点中任意组合三角形，判断网格点是否落在三角形内部，若网格点落在2个以上三角形内部，则计算网格点与三角形顶点的总距离，取距离最短的三角形；若再次判断失败，则继续寻找下1个最邻近点，重复上述步骤，若在搜索半径内，无法构成该三角形，则该网格点失败；

判断网格点是否在三角形内可以使用向量同向法：若点P在三角形ABC内，则三个向量AB×AP、AP×AC、PB×PC平行同向，它们也与向量AB×AC平行同向，由于这三个向量均有可能为0，直接判断它们平行同向比较麻烦，但考虑到AB×AC不可能为0，直接判断向量：AB×AP、AP×AC、PB×PC均与AB×AC平行同向更为简单；

三角形确定后，利用三角形区域网函数插公式求取网格点的值；其中Q_i是直线M_iQ和M_jM_k之交点，i,j,k∈{1,2,3}，F(Q)代表插值后的数值，Δ表示为三角形的总面积，Δ_1、Δ_2、Δ₃分别代表三角形M₃QM₂、M₃QM₁、M₁QM₂的面积；

然后，依次搜寻下一个网格点的最小外接三角形，并进行插值运算，方法同上，直到所有网格点都遍历完毕，形成规则的海面高数据点阵；

d、精细插值

基于网格化的卫星测高数据,利用矩形区域网函数插值算法内插出网格内部某点的海面高度,进一步提高海面高数据的空间分辨率，矩形面积加权网函数插值算法如下：

P₁(x₀,y₀)，P₂(x₁,y₀)，P₃(x₁,y₁)，P₄(x₀,y₁)是矩形区域D的四个角点(即相邻的四个网格数据点)，Q_i(i＝1,2,3,4)是过插值点Q与矩形平行的两条直线在矩形边界上截得的四个点；A_i(i＝1,2,3,4)是这两条直线将矩形区域D分成四个小矩形的面积，其中A₁＝(x-x₀)(y-y₀),A₂＝(x₁-x)(y-y₀),A₃＝(x₁-x)(y₁-y)，A₄＝(x–x₀)(y₁-y)；A＝(x₁–x₀)(y₁-y₀)＝A₁+A₂+A₃+A₄，是矩形区域D的总面积；

则Q点的海面高度值 <mrow> <mi>F</mi> <mrow> <mo>(</mo> <mi>Q</mi> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>4</mn> </msubsup> <mfrac> <mrow> <msub> <mi>A</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>A</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> </mrow> <mi>A</mi> </mfrac> <msub> <mi>Q</mi> <mi>i</mi> </msub> <mo>-</mo> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>4</mn> </msubsup> <mfrac> <msub> <mi>A</mi> <mi>i</mi> </msub> <mi>A</mi> </mfrac> <msub> <mi>P</mi> <mi>i</mi> </msub> <mo>;</mo> </mrow> 1

以插值n倍为例，插值点坐标的计算方法为：x＝x₀+W/n，y＝y₀+W/n，其中W为网格宽度；

确定n后，逐行遍历每个网格点，针对每个网格分别计算插值点Q(x,y)的海面高度值，直至所有网格点遍历结束，生成新的平均海面高数据集。