使用径向基函数在球面上插值函数
首先，给大家一点背景知识：
我正在使用球面调和函数作为一个例子，这个函数是在球体表面上的，就像这张图片中的前面那些球体：
我制作了其中一个球体，并根据球面上各个点的调和函数值给它上了色。我一开始是在很多点上进行计算，这样我的函数就非常准确。我把这个球体称为我的 fine 球体。
现在我有了我的 fine 球体，我在球体上选择了相对较少的点。这些点是我想要进行插值的点，也就是我的训练数据，我称它们为 interp 点。下面是我的 interp 点，按照它们的值上色，并绘制在我的 fine 球体上。
这个项目的目标是使用这些 interp 点来训练一个 SciPy 径向基函数，以便在球体上插值我的函数。我用以下代码实现了这个目标：
# Train the interpolation using interp coordinates
rbf = Rbf(interp.phi, interp.theta, harmonic13_coarse)
# The result of the interpolation on fine coordinates
interp_values = rbf(fine.phi, fine.theta)
这段代码生成了这个插值，绘制在球体上：
希望通过最后这张图片，你能看到我的问题。注意插值中间的那条线吗？这是因为插值数据有边界。边界的产生是因为我使用球坐标训练了径向基函数（边界在 [0,pi] 和 [0,2pi]）。
rbf = Rbf(interp.phi, interp.theta, harmonic13_coarse)
我的目标，以及我发帖的原因，是希望使用球体上数据的 x,y,z 笛卡尔坐标来插值我的函数。这样，由于球体没有边界，我就不会像使用球坐标时那样出现边界错误。然而，我就是搞不清楚该怎么做。
我尝试直接给 Rbf 函数提供 x,y,z 坐标和函数值。
rbf=Rbf(interp.x, interp.y, interp.z, harmonic13_coarse)
interp_values=rbf(fine.x,fine.y,fine.z)
但是 NumPy 给我抛出了一个奇异矩阵错误。
numpy.linalg.linalg.LinAlgError: singular matrix
有没有办法让我用笛卡尔坐标提供 Rbf 我的数据点，以及每个点的函数值，并让它像处理球坐标那样工作，但没有那些边界呢？从 Rbf 的文档来看，有一个 norm 属性可以定义不同的距离规范，我是否需要使用球面距离来让这个工作？
我对此感到很困惑。如果你有任何关于如何在没有球坐标边界的情况下在球体上插值我的函数的想法，请告诉我。
这是我完整的代码：
import matplotlib.pyplot as plt
from matplotlib import cm, colors
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from scipy import special
from scipy.interpolate import Rbf
from collections import namedtuple
from mayavi import mlab
# Nice aliases
pi = np.pi
cos = np.cos
sin = np.sin
# Creating a sphere in Cartesian and Sphereical
# Saves coordinates as named tuples
def coordinates(r, n):
phi, theta = np.mgrid[0:pi:n, 0:2 * pi:n]
Coor = namedtuple('Coor', 'r phi theta x y z')
r = r
x = r * sin(phi) * cos(theta)
y = r * sin(phi) * sin(theta)
z = r * cos(phi)
return Coor(r, phi, theta, x, y, z)
# Creating a sphere
# fine is coordinates on a fine grid
# interp is coordinates on coarse grid for training interpolation
fine = coordinates(1, 100j)
interp = coordinates(1, 5j)
# Defining finection to colour sphere
# Here we are using a spherical harmonic
def harmonic(m, n, theta, phi):
return special.sph_harm(m, n, theta, phi).real
norm = colors.Normalize()
# One example of the harmonic function, for testing
harmonic13_fine = harmonic(1, 3, fine.theta, fine.phi)
harmonic13_coarse = harmonic(1, 3, interp.theta, interp.phi)
# Train the interpolation using interp coordinates
rbf = Rbf(interp.phi, interp.theta, harmonic13_coarse)
# The result of the interpolation on fine coordinates
interp_values = rbf(fine.phi, fine.theta)
rbf=Rbf(interp.x, interp.y, interp.z, harmonic13_coarse)
interp_values=rbf(fine.x,fine.y,fine.z)
#Figure of harmoinc function on sphere in fine cordinates
#Points3d showing interpolation training points coloured to their value
mlab.figure()
vmax, vmin = np.max(harmonic13_fine), np.min(harmonic13_fine)
mlab.mesh(fine.x, fine.y, fine.z, scalars=harmonic13_fine, vmax=vmax, vmin=vmin)
mlab.points3d(interp.x, interp.y, interp.z, harmonic13_coarse,
scale_factor=0.1, scale_mode='none', vmax=vmax, vmin=vmin)
#Figure showing results of rbf interpolation
mlab.figure()
vmax, vmin = np.max(harmonic13_fine), np.min(harmonic13_fine)
mlab.mesh(fine.x, fine.y, fine.z, scalars=interp_values)
# mlab.points3d(interp.x, interp.y, interp.z, scalars, scale_factor=0.1, scale_mode='none',vmax=vmax, vmin=vmin)
mlab.show()
2 个回答
你可以改变标准的距离计算方式：
def euclidean_norm(x1, x2):
return np.sqrt( ((x1 - x2)**2).sum(axis=0) )
可以使用球面距离（比如，可以参考这个问题 Python中的哈弗辛公式（计算两个GPS点之间的方位和距离））。
你看到的边界是因为你把一个封闭的表面（S2）映射到了一个开放的表面（R2）。无论怎样，你都会遇到边界。虽然这些表面的局部特性是兼容的，所以在大部分球面上是可以工作的，但在整体上就不行了，你会得到一条线。
解决这个问题的方法是使用一个图集，而不是单一的图表。图集是多个重叠图表的集合。在重叠的区域，你需要定义权重，也就是一个平滑的函数，它在每个图表上从0到1变化。（抱歉，可能你并不想听到微分几何的内容）。
如果你不想深入到这个问题，可以注意到你原来的球面上有一个赤道，那里方差是最小的。然后你可以旋转你的球面，使它和那条线重合。这样虽然不能完全解决你的问题，但肯定能减轻一些。