np.gradient处理非均匀间距平滑数据方差过高的解决方法咨询

处理相对平滑的数据时，使用np.gradient计算dx/dt（时间t为非均匀间距），结果方差较高，数据中的噪声被大幅放大。复现代码如下：

import numpy as np
import matplotlib.pyplot as plt
x = np.array([13.11149679, 13.2141427 , 13.37743691, 13.3934357 , 13.56163066,
       13.60207566, 13.69304133])
t = np.array([0.73065159, 0.74012055, 0.75911018, 0.7607452 , 0.77811468,
       0.78031837, 0.79046324])
x_grad = np.gradient(x, t)

plt.plot(t[1:], np.diff(x) / np.diff(t), 'xb', label="findiff")
plt.plot(t, x_grad, 'or', label = 'np.gradient')
plt.plot(t, x, '+g', label = "x")
plt.xlabel('t')
plt.ylabel('x')
plt.legend()
plt.show()

问题原因 #

np.gradient对非均匀间隔数据采用相邻点加权平均计算梯度，当时间间隔突变（比如你的t数组中0.7607到0.7781的间隔突然变大）时，微小的数据波动会被加权机制放大，导致梯度结果噪声显著。

解决方案 #

方法一：先平滑数据再求梯度

用滤波方法先抑制数据中的噪声，再计算梯度，常用的有Savitzky-Golay滤波、滑动窗口平均。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter

x = np.array([13.11149679, 13.2141427 , 13.37743691, 13.3934357 , 13.56163066,
       13.60207566, 13.69304133])
t = np.array([0.73065159, 0.74012055, 0.75911018, 0.7607452 , 0.77811468,
       0.78031837, 0.79046324])

# 用Savitzky-Golay滤波平滑x，窗口大小3，一阶多项式拟合
x_smoothed = savgol_filter(x, window_length=3, polyorder=1)
x_grad_smoothed = np.gradient(x_smoothed, t)

plt.plot(t[1:], np.diff(x)/np.diff(t), 'xb', label="findiff")
plt.plot(t, np.gradient(x, t), 'or', label='np.gradient (raw)')
plt.plot(t, x_grad_smoothed, '^g', label='np.gradient (smoothed)')
plt.plot(t, x, '+k', label="x")
plt.xlabel('t')
plt.ylabel('x / dx/dt')
plt.legend()
plt.show()

方法二：多项式拟合后求导

对数据做连续曲线拟合（比如三次样条），再对拟合曲线求导，这种方法对非均匀间隔数据更稳健。

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

x = np.array([13.11149679, 13.2141427 , 13.37743691, 13.3934357 , 13.56163066,
       13.60207566, 13.69304133])
t = np.array([0.73065159, 0.74012055, 0.75911018, 0.7607452 , 0.77811468,
       0.78031837, 0.79046324])

# 三次样条拟合x(t)，并求一阶导数
cs = CubicSpline(t, x)
t_fine = np.linspace(t.min(), t.max(), 100)
x_grad_cs = cs(t_fine, 1)

plt.plot(t[1:], np.diff(x)/np.diff(t), 'xb', label="findiff")
plt.plot(t, np.gradient(x, t), 'or', label='np.gradient (raw)')
plt.plot(t_fine, x_grad_cs, '-g', label='CubicSpline derivative')
plt.plot(t, x, '+k', label="x")
plt.xlabel('t')
plt.ylabel('x / dx/dt')
plt.legend()
plt.show()

方法三：自定义加权梯度

给相邻点分配与间隔成反比的权重，减少间隔突变对梯度的影响，替代np.gradient的默认加权逻辑。

import numpy as np
import matplotlib.pyplot as plt

x = np.array([13.11149679, 13.2141427 , 13.37743691, 13.3934357 , 13.56163066,
       13.60207566, 13.69304133])
t = np.array([0.73065159, 0.74012055, 0.75911018, 0.7607452 , 0.77811468,
       0.78031837, 0.79046324])

def weighted_gradient(x, t):
    grad = np.zeros_like(x)
    # 首尾点用前后向差分
    grad[0] = (x[1] - x[0])/(t[1] - t[0])
    grad[-1] = (x[-1] - x[-2])/(t[-1] - t[-2])
    # 中间点用间隔反比加权平均
    for i in range(1, len(x)-1):
        dt_left = t[i] - t[i-1]
        dt_right = t[i+1] - t[i]
        w_left = dt_right / (dt_left + dt_right)
        w_right = dt_left / (dt_left + dt_right)
        grad[i] = w_left*(x[i]-x[i-1])/dt_left + w_right*(x[i+1]-x[i])/dt_right
    return grad

x_grad_weighted = weighted_gradient(x, t)

plt.plot(t[1:], np.diff(x)/np.diff(t), 'xb', label="findiff")
plt.plot(t, np.gradient(x, t), 'or', label='np.gradient (raw)')
plt.plot(t, x_grad_weighted, 'sm', label='Weighted gradient')
plt.plot(t, x, '+k', label="x")
plt.xlabel('t')
plt.ylabel('x / dx/dt')
plt.legend()
plt.show()

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