如何用插值法平滑卫星通道权重函数曲线？遇重复值报错求解

问题背景 #

绘制卫星部分通道的权重函数时，尝试用scipy.interpolate.make_interp_spline实现平滑拟合效果，运行代码触发了重复x值的错误。未做平滑处理时已能生成多通道阶梯状曲线（纵轴为倒置的气压0-1000hPa，横轴为归一化权重）。

原始代码 #

from scipy.interpolate import make_interp_spline, BSpline
import numpy as np
import matplotlib.pyplot as plt
import example_data as ex

# 读取透过率数据
tau_level = np.squeeze(np.load('/home/swadhin/rttov/wrapper/tau_level.npy'))
print(np.max(tau_level))

# 计算不同层的透过率差值
ans = []
for i in range(len(tau_level)):
    l = []
    for j in range(len(tau_level[i])-1):
        a = tau_level[i][j]-tau_level[i][j+1]
        l.append(a)
    ans.append(l)
print(np.max(ans))

p = []
p_plt = []
s = ex.p_ex[0]  # 从数据文件导入气压廓线
for x in range(len(s)-1):
    # 计算连续气压层的对数差绝对值
    y = np.abs(np.log(s[x]) - np.log(s[x+1]))
    # 取两个气压层的中点作为该层的绘图代表值
    p1 = (s[x]+s[x+1])/2
    p.append(y)
    p_plt.append(p1)
print(np.shape(p))

# 计算不同通道的权重函数
plot = np.array([i / j for i, j in zip(ans, p)])

# 绘制未平滑的权重函数
x_plt_store = []
for i in plot:
    x_plt = np.array(i/max(i))
    plt.plot(x_plt, p_plt)
    x_plt_store.append(x_plt)
x_plt_store = np.nan_to_num(x_plt_store)  # 处理nan/inf值

# 插值实现平滑绘图（原报错部分）
x_new_store = []
y_spl = []
for i in x_plt_store:
    x_new = np.linspace(i.min(),i.max(),500)
    x_new_store.append(x_new)
    spl = make_interp_spline(i, p_plt, k=3)
    y_spl.append(spl)

p_smooth = y_spl(x_new_store)
plt.plot(x_new_store, p_smooth)
plt.yticks(np.arange(0, 1000, 100))
plt.axis([0, 1, 0, 1000])
plt.gca().invert_yaxis()
plt.show()
plt.savefig('/home/swadhin/rttov/wrapper/wfn_insat.png', format='png', dpi=1000)

报错信息 #

ValueError                                Traceback (most recent call last)
/home/swadhin/rttov/wrapper/wfn.py in <cell line: 50>()

     51     x_new = np.linspace(i.min(),i.max(),500) # 500 represents number of points to make between i.min and i.max
     52     x_new_store.append(x_new)
---&gt; 53     spl = make_interp_spline(i, p_plt, k=3)
     54     y_spl.append(spl)
     56 p_smooth = y_spl(x_new_store)

File ~/anaconda3/envs/rttov/lib/python3.9/site-packages/scipy/interpolate/_bsplines.py:1067, in make_interp_spline(x, y, k, t, bc_type, axis, check_finite)
   1065     raise ValueError("Expect x to be a 1-D sorted array_like.")
   1066 if np.any(x[1:] == x[:-1]):
-&gt; 1067     raise ValueError("Expect x to not have duplicates")
   1068 if k &lt; 0:
   1069     raise ValueError("Expect non-negative k.")

ValueError: Expect x to not have duplicates

解决方法 #

报错原因 #

make_interp_spline要求输入的自变量x必须是无重复的单调序列，但原代码把归一化权重值作为x输入，而权重函数中存在多个相同的归一化值（比如多个层的权重为0，归一化后仍为0），导致重复x值触发错误。

修正思路 #

权重函数是随气压变化的，气压p_plt是单调递减的序列（无重复值），完全符合插值要求。因此应该以p_plt为自变量x，归一化权重为因变量y，对y进行平滑插值。

修正后的完整代码 #

from scipy.interpolate import make_interp_spline, BSpline
import numpy as np
import matplotlib.pyplot as plt
import example_data as ex

# 读取透过率数据
tau_level = np.squeeze(np.load('/home/swadhin/rttov/wrapper/tau_level.npy'))

# 计算不同层的透过率差值
ans = []
for i in range(len(tau_level)):
    l = []
    for j in range(len(tau_level[i])-1):
        a = tau_level[i][j]-tau_level[i][j+1]
        l.append(a)
    ans.append(l)

p = []
p_plt = []
s = ex.p_ex[0]  # 从数据文件导入气压廓线
for x in range(len(s)-1):
    # 计算连续气压层的对数差绝对值
    y = np.abs(np.log(s[x]) - np.log(s[x+1]))
    # 取两个气压层的中点作为该层的绘图代表值
    p1 = (s[x]+s[x+1])/2
    p.append(y)
    p_plt.append(p1)
p_plt = np.array(p_plt)  # 转换为numpy数组，方便后续操作

# 计算不同通道的权重函数
plot = np.array([i / j for i, j in zip(ans, p)])

plt.figure(figsize=(8,6))
# 绘制未平滑的权重函数（可选，作为对比）
for i in plot:
    x_plt = np.array(i/max(i))
    plt.plot(x_plt, p_plt, alpha=0.3)

# 插值实现平滑绘图（修正部分）
x_new = np.linspace(p_plt.min(), p_plt.max(), 500)  # 以气压为自变量生成平滑点
for i in plot:
    x_plt = np.array(i/max(i))
    # 以气压p_plt为x，归一化权重为y，生成样条插值
    spl = make_interp_spline(p_plt, x_plt, k=3)
    x_smooth = spl(x_new)
    plt.plot(x_smooth, x_new)

plt.yticks(np.arange(0, 1000, 100))
plt.axis([0, 1, 0, 1000])
plt.gca().invert_yaxis()
plt.savefig('/home/swadhin/rttov/wrapper/wfn_insat.png', format='png', dpi=1000)
plt.show()

关键修改点 #

1. 将插值的自变量从归一化权重改为气压p_plt，彻底避免重复值问题
2. 统一使用numpy数组操作，提升代码稳定性
3. 平滑点生成基于气压范围，确保曲线覆盖整个气压区间
4. 保留未平滑曲线作为对比（可通过注释移除）

内容的提问来源于stack exchange，提问作者The Emerging Star