球面网格之间的一种双线性插值方法

原文：An Alternative Bilinear Interpolation Method Between Spherical Grids

摘要

在地球科学研究中，传统的双线性插值已被广泛用于球体表面逻辑矩形网格之间的重新映射。最近，已经提出了包括测地线网格在内的各种球形网格系统来解决由传统经纬度网格引起的奇异性问题。我们建议使用一种替代预先存在的双线性插值方法来重新映射任何球面网格，甚至是球体上随机分布的点。该方法支持与目标点相邻的四个源点的任何几何配置进行插值，并提供与传统双线性方法等效的重映射精度。此外，为了有效搜索相邻源点，我们使用以立方球体作为参考系的链接单元映射方法。因此，计算成本与 NlogN 而不是 $N^2$（N 是网格点的数量）成正比，即使对于球体上随机分布的点的重新映射也是如此。

1.介绍

已经提出了各种球形网格系统，例如立方球体 [1] 和二十面体 [2]，以解决传统经纬度网格的奇异性问题。地球系统建模本质上需要在球形网格之间重新映射或重新网格化 [3]。有多种重新映射方法可以在球形网格之间进行插值，例如双线性 [4]、双三次 [4]、保守 [5,6] 和补丁恢复 [7] 方法。地球系统建模框架 (ESMF) 重新网格化包和 SCRIP（球面坐标重映射和插值包）是众所周知的重映射工具包，并支持这些不同的重映射方法。
双线性插值因其易于实现且计算成本低而被广泛用于逻辑矩形网格。尽管 ESMF 和 SCRIP 重映射工具包都支持双线性插值，但它们有一些限制。SCRIP 中的双线性重映射工具仅适用于逻辑矩形网格 [8]，ESMF包中的一个支持逻辑矩形网格和由多边形组成的非结构化网格，但不支持自相交单元； 例如，一个细胞扭曲成领结形状 [9]。
在这项研究中，我们提出了一种替代的双线性插值方法，该方法可用于任何球面网格，甚至包括随机分布的点，保持其实现的简单性，并且计算成本低。 特别是，我们提出了一种有效的搜索算法，以使用链接单元映射方法和立方球体作为参考系来查找相邻网格点。由于立方球体是其中元素均匀分布在球体上的球形网格之一，因此它可以为搜索提供合理大小的网格。 有关该方法的更多详细信息，请参见第 3 节。
我们通过测试最初使用纬度-经度系统、立方球体、斐波那契球体和随机分布点的球谐函数的重新映射来验证新建议的双线性插值方法。测试结果表明，该方法的精度与现有的双线性插值方法相当或略好。 此外，我们发现计算成本与 NlogN 成正比，而不是 $N^2$，其中 N 是网格点的数量。下一节将详细介绍这种双线性插值方法。第 3 节介绍了基于链接单元映射方法的搜索算法。我们描述了用于验证新建议的双线性插值方法的测试集，并在第 4 节中将其精度与经纬度系统、立方球体、斐波那契球体和随机分布点上的现有方法进行了比较。在最后一节，我们总结了新的插值方法和搜索算法的可行性。

2.一种可选的双线性插值方法

图 1 显示了使用围绕目标点的四个矩形点进行常规双线性插值的目标点的配置。该过程是简单的线性插值，首先沿 x 轴，然后沿 y 轴。以目标点为中心，插值的权重由四个矩形的面积与其面积总和的比值给出。这种加权方式与传统的双线性插值方法相同[4]。
$f(x,y)=w_1{f}_1+w_2{f}_2+w_3{f}_3+w_4{f}_4$,
其中：

图1：在由四个网格点 (x1,y1)、(x2,y2)、(x3,y3) 和 (x4,y4) 形成的矩形星座中，函数 f 在位置 (x, y) 上的双线性插值示例使用它们的 对应的函数 f1、f2、f3 和 f4。
同时，图 2 展示了使用四个不形成矩形的点对目标进行双线性插值的示例，并表明本研究中建议的双线性插值方法适用于目标点附近源点的任何几何配置。现在，我们解释这个假设和本研究中提出的双线性插值背后的想法。 插值函数由下式给出
$f(x,y)=a+bx+cy+dxy$,
图2：函数 f 在位置 (x, y) 上的替代双线性插值示例：(a) 四边形和 (b) 四个网格点 (x1,y1)、(x2,y2)、(x3,y3) 、(x4,y4) 使用它们对应的函数 f1、f2、f3 和 f4。
函数的选择是基于它尽可能简单并且如果四个源点形成一个矩形，则它与传统的双线性插值相同。这四个源点的坐标分别用(x1,y1)、(x2,y2)、(x3,y3)和(x4,y4)表示，它们对应的函数值为f1、f2、f3和f4， 分别。可以形成方程组并将其表示为矩阵和向量的乘法：

#一段python的程序，计算一下插值的效果，还可以，但是也没有那么理想
import numpy as np  #导入numpy的库函数
x1=2
y1=1
f1=0
x2=9
y2=1
f2=15
x3=1
y3=15
f3=20
x4=10
y4=15
f4=30
xt=5.5
yt=14
c=0
D=np.array([[1,x1,y1,x1*y1],[1,x2,y2,x2*y2],[1,x3,y3,x3*y3],[1,x4,y4,x4*y4]])
b=np.array([f1,f2,f3,f4])
x=np.linalg.solve(D,b)
print(x)
aaa=x[0]
bbb=x[1]
ccc=x[2]
ddd=x[3]
f=aaa+bbb*xt+ccc+yt+ddd*xt*yt
print(f)  //16.23

一旦方程组已经求解，就可以得到相应函数的系数 a、b、c 和 d，并且函数 f(x, y) 在这四个附近任意点的位置上的双线性插值 点可以继续。特别是，应检查以下两个关键点，以确保插值的精度与传统双线性插值的精度相当。
(1) 我们需要通过旋转四个源点的 x-y 轴来寻找一个行列式最大的 4 × 4 矩阵，以获得与插值相关的坐标。换句话说，这种双线性插值过程的精度取决于所选择的四个相邻点的坐标。根据用于方程组的 Cramer 规则 (Apostol, 1969)，系数可以由行列式 D1、D2、D3 和 D4 与行列式 D 的比率给出，它们的定义如下：
随后，方程组中的系数可以表示为：

例如，让我们考虑一个相对于原点和与原始矩阵的角度 θ 变换的旋转矩阵，给定目标点和前面提到的四个相邻点。 在此之后，行列式 D 和 D4 变为：

例如，让我们考虑一个相对于原点和与原始矩阵的角度 θ 变换的旋转矩阵，给定目标点和前面提到的四个相邻点。 在此之后，行列式 D 和 D4 变为
这里，D4 与四个源点的坐标变化无关，因此当我们以行列式变得最大的方式旋转源点的轴时，系数变得最小。插值的精度主要受系数 d 的影响，因为当 d 最小时，即当 D 最大时精度最高。 从几何上讲，当 d 最小时，该函数看起来是平坦的。