生成5000条过指定点的平滑地球层界面曲线的高效方法咨询

需要生成约5000条不同的平滑曲线，所有曲线必须经过一组预设点（用于表示地球内部层界面）。核心要求：

• 相邻预设点之间不能有明显振荡或突变（符合地球内部结构渐变的特点）
• 曲线的y值限制在25-40之间
• 现有方法是在预设点间随机生成点再拟合，效率较低，希望找到更高效的实现方式

附现有代码示例：

import numpy as np
from scipy import interpolate

x = [7.81, 21.65, 186.29, 246.62, 402.74, 446.24, 572.99, 585.11, 613.57, 762.97]
y = [26.27, 26.41, 26.46, 27.36, 34.32, 34.50, 37.94, 39.05, 39.23, 39.88]

# y can only range between 25 and 40

fig = plt.figure(figsize=(20,5))
xnew = np.linspace(min(x), max(x), num=50)
oned = interpolate.CubicSpline(x, y)
yfit = oned(xnew)
plt.scatter(x, y, color='red')
plt.plot(xnew, yfit)
plt.xlim(0, 800)
plt.ylim(20,45)
plt.gca().invert_yaxis()
plt.show()

示例曲线：

一、无振荡插值方法推荐 #

1. 单调三次插值（Pchip） #

scipy的interpolate.PchipInterpolator（或pchip）专为单调数据设计，能保证曲线在相邻点间保持单调性，完全避免振荡，适配地球内部结构渐变的需求。它会自动调整插值斜率，确保曲线平滑且无过冲。

2. 分段线性插值+平滑滤波 #

先做分段线性插值，再用scipy.ndimage.gaussian_filter1d等低通滤波工具平滑曲线。这种方法简单直观，可通过调整滤波参数控制平滑度，同时严格保证曲线经过预设点。

3. 约束样条插值 #

使用scipy.interpolate.LSQUnivariateSpline，通过设置平滑因子s控制曲线平滑度，还可添加y值范围、导数限制等约束条件，避免振荡。需注意调整s的值，确保曲线严格经过预设点。

二、高效生成5000条曲线的思路 #

无需在预设点间随机生成点再拟合，直接在插值参数空间引入可控随机扰动，同时保证曲线经过预设点且平滑：

• Pchip斜率扰动：计算原始Pchip的各点斜率，对内部点斜率添加小范围随机扰动（需保证扰动后曲线单调、y值在限制范围内），用新斜率生成Pchip曲线。
• 约束样条权重调整：给每个预设点添加微小随机权重（权重不为0），微调平滑因子s，生成不同平滑曲线。
• 分段二次曲线偏移：将相邻预设点间的曲线设为二次函数，在保证经过端点的前提下，随机调整二次项系数（控制弯曲程度，不超出y值范围）。

三、代码实现示例（基于Pchip的高效生成） #

import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate

# 预设点
x = np.array([7.81, 21.65, 186.29, 246.62, 402.74, 446.24, 572.99, 585.11, 613.57, 762.97])
y = np.array([26.27, 26.41, 26.46, 27.36, 34.32, 34.50, 37.94, 39.05, 39.23, 39.88])
y_min, y_max = 25, 40

# 生成5000条平滑曲线的函数
def generate_smooth_curves(x, y, num_curves=5000, num_points=50, slope_perturb=0.05):
    curves = []
    # 计算原始Pchip的斜率
    pchip_orig = interpolate.PchipInterpolator(x, y)
    slopes_orig = pchip_orig.derivative()(x)
    
    for _ in range(num_curves):
        slopes_perturbed = slopes_orig.copy()
        # 仅扰动中间点的斜率
        for i in range(1, len(x)-1):
            # 计算当前点允许的斜率范围，保证曲线不越界且单调
            max_slope = (y_max - y[i]) / min(x[i+1]-x[i], x[i]-x[i-1])
            min_slope = (y_min - y[i]) / max(x[i+1]-x[i], x[i]-x[i-1])
            # 添加随机扰动并限制在允许范围内
            delta = np.random.uniform(-slope_perturb, slope_perturb)
            slopes_perturbed[i] = np.clip(slopes_orig[i] + delta, min_slope, max_slope)
        
        # 生成新的Pchip曲线
        pchip_new = interpolate.PchipInterpolator(x, y, der=slopes_perturbed)
        x_new = np.linspace(x.min(), x.max(), num_points)
        y_new = pchip_new(x_new)
        # 最终确保y值在限制范围内
        y_new = np.clip(y_new, y_min, y_max)
        curves.append((x_new, y_new))
    
    return curves

# 生成曲线并展示5条示例
curves = generate_smooth_curves(x, y, num_curves=5000)

fig = plt.figure(figsize=(20,5))
plt.scatter(x, y, color='red', zorder=10, label='预设点')
for i in range(5):
    plt.plot(curves[i][0], curves[i][1], alpha=0.6, label=f'曲线{i+1}')
plt.xlim(0, 800)
plt.ylim(20,45)
plt.gca().invert_yaxis()
plt.legend()
plt.show()

四、注意事项 #

• 随机扰动需设置合理范围，避免曲线超出y值限制或出现突变
• 可通过multiprocessing并行生成曲线，进一步提升5000条曲线的生成效率
• 需严格单调性约束时优先选Pchip；允许轻微非单调但无振荡时，约束样条更灵活

内容的提问来源于stack exchange，提问作者Thoram