这取决于您认为合适的插值方法和物理特性。
从你所写的内容来看，我倾向于假设与波长误差相比，空间维度沿着误差可以忽略不计。如果是这样的话，N维插值可能是错误的，因为像素信息应该是独立的。相反，你需要做的是对所有像素进行一维插值。
最简单（也是最快）的插值形式是最近邻插值。现在，如果新的波长可以用np.round(decimals=-1)计算，则数据已经插值，您只需要更新波长值。
如果新波长不是旧波长的四舍五入值，或者如果不需要或不想要最近邻插值，则需要使用不同的方法，在某些点上，这将涉及循环通过像素。
SciPy提供了scipy.interpolate.interp1d()，它以一种矢量化的方式（即，通过像素的循环被推到Python帧之外）完成了这一任务，并提供了各种插值方法。
例如，如果samples是测量的波长，new_samples包含新的波长，arr包含最后一个轴穿过波长的堆叠图像：

import scipy.interpolate

def interp_1d_last_sp(arr, samples, new_samples, kind="linear"):
    interpolator = scipy.interpolate.interp1d(samples, arr, axis=-1, kind=kind, fill_value="extrapolate")
    return interpolator(new_samples)

至少对于线性插值，可以使用更快的np.interp()函数，通过显式循环手动计算类似的结果：

import numpy as np

def interp_1d_last_np(arr, samples, new_samples):
    shape = arr.shape
    k = shape[-1]
    arr = arr.reshape((-1, k))
    result = np.empty_like(arr)
    n = arr.shape[0]
    for i in range(n):
        result[i, :] = np.interp(new_samples, samples, arr[i, :])
    return result.reshape(shape)

该函数不支持多维输入，但可以使用Numba：

import numba as nb

@nb.njit(parallel=True)
def interp_1d_last_nb(arr, samples, new_samples):
    shape = arr.shape
    k = shape[-1]
    arr = arr.reshape((-1, k))
    result = np.empty_like(arr)
    n = arr.shape[0]
    for i in nb.prange(n):
        result[i, :] = np.interp(new_samples, samples, arr[i, :])
    return result.reshape(shape)

虽然它们都得到了类似的结果（它们处理外推的方式略有不同），但时间可能不同：

np.random.seed(0)
a = np.random.random([240, 320, 30])
w_min = 400
w_max = 700
w_step = 10
w = np.arange(w_min, w_max, w_step)
nw = w_step * (np.random.random(w.size) - 0.5)

funcs = interp_1d_last_sp, interp_1d_last_np, interp_1d_last_nb
base = funcs[0](a, w + nw, w)
for func in funcs:
    res = func(a, w + nw, w)
    # the sum of the absolute difference with the non-interpolated array is reasonably the same
    is_good = np.isclose(np.sum(np.abs(base - a)), np.sum(np.abs(res - a)))
    print(f"{func.__name__:>12s}  {is_good!s:>5}  {np.sum(np.abs(res - a)):16.8f}  ", end="")
    %timeit -n 2 -r 2 func(a, w + nw, w)

# interp_1d_last_sp   True   140136.05282911  90.8 ms ± 21.7 ms per loop (mean ± std. dev. of 2 runs, 2 loops each)

# interp_1d_last_np   True   140136.05282911  349 ms ± 14.2 ms per loop (mean ± std. dev. of 2 runs, 2 loops each)

# interp_1d_last_nb   True   140136.05282911  61.5 ms ± 1.83 ms per loop (mean ± std. dev. of 2 runs, 2 loops each)

编辑了解决方案，使其完全矢量化，没有循环。可能更多的内存密集型虽然。
注意，矢量化的版本比仅仅循环通过波长要慢。

import numpy as np
import scipy.interpolate as interpolate

a = np.random.random([240, 320, 30])

wavelengths1 = np.array([395.13, 408.62, 421.63, 434.71, 435.64, 453.39, 456.88, 471.48,
       484.23, 488.89, 497.88, 513.35, 521.38, 528.19, 539.76, 548.39,
       557.78, 568.06, 577.64, 590.22, 598.63, 613.13, 618.87, 632.75,
       637.5 , 647.47, 655.6 , 672.66, 681.88, 690.1 ])
wavelengths2 = np.array([400, 410, 420, 430, 440, 450, 460, 470, 480, 490, 500, 510, 520,
       530, 540, 550, 560, 570, 580, 590, 600, 610, 620, 630, 640, 650,
       660, 670, 680, 690])

interp = interpolate.RegularGridInterpolator((np.arange(a.shape[1]), np.arange(a.shape[0]), wavelengths1), a.transpose((1, 0, 2)))

# Create new grid based on the integer wavelengths desired above

X, Y = np.meshgrid(np.arange(320), np.arange(240))

X2 = np.repeat(X[:, :, None], len(wavelengths2), axis=-1)

Y2 = np.repeat(Y[:, :, None], len(wavelengths2), axis=-1)

a2 = interp((Y2, X2, wavelengths2[None, None, :]))