python - 使用 numpy/scipy 识别数字信号的斜率变化？

Question

我正尝试在 Python 中提出一种通用方法来识别在一组计划的航天器机动过程中发生的俯仰旋转。您可以将其视为 shift detection 的特例问题。

让我们考虑一下我的测量集中的 solar_elevation_angle 变量，确定从航天器仪器测量的太阳仰角。对于那些可能想要玩转数据的人，我保存了 solar_elevation_angle.txt 文件 here .

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import gridspec
from scipy.signal import argrelmax
from scipy.ndimage.filters import gaussian_filter1d

solar_elevation_angle = np.loadtxt("solar_elevation_angle.txt", dtype=np.float32)

fig, ax = plt.subplots()    
ax.set_title('Solar elevation angle')
ax.set_xlabel('Scanline')
ax.set_ylabel('Solar elevation angle [deg]')
ax.plot(solar_elevation_angle)
plt.show()

扫描线是我的时间维度。斜率变化的四个点确定了航天器的俯仰旋转。

如您所见，航天器机动区域外的太阳高度角随时间的变化几乎是线性的，对于这个特定的航天器来说应该始终如此(重大故障除外)。

请注意，在每次航天器操纵期间，斜率变化显然是连续的，尽管在我的一组角度值中是离散的。这意味着:对于每次机动，尝试定位发生机动的单个扫描线实际上没有意义。我的目标是为每个操作识别扫描线范围内的“代表性”扫描线，定义操作发生的时间间隔(例如中间值或左边界)。

一旦我获得一组“代表性”扫描线索引，其中所有机动都已发生，我就可以使用这些索引粗略估计机动持续时间，或自动在绘图上放置标签。

到目前为止，我的解决方案是:

1. 使用以下方法计算太阳高度角的二阶导数 np.gradient。
2. 计算结果的绝对值和剪裁 曲线。剪辑是必要的，因为我认为是 线性段中的离散化噪声，这将严重影响点 4 中“真实”局部最大值的识别。
3. 对生成的曲线应用平滑处理，以去除多个峰值。我正在使用 scipy 的 1d 高斯滤波器和试错 sigma 值。
4. 确定局部最大值。

这是我的代码:

fig = plt.figure(figsize=(8,12))
gs = gridspec.GridSpec(5, 1) 

ax0 = plt.subplot(gs[0])
ax0.set_title('Solar elevation angle')
ax0.plot(solar_elevation_angle)

solar_elevation_angle_1stdev = np.gradient(solar_elevation_angle)
ax1 = plt.subplot(gs[1])
ax1.set_title('1st derivative')
ax1.plot(solar_elevation_angle_1stdev)

solar_elevation_angle_2nddev = np.gradient(solar_elevation_angle_1stdev)
ax2 = plt.subplot(gs[2])
ax2.set_title('2nd derivative')
ax2.plot(solar_elevation_angle_2nddev)

solar_elevation_angle_2nddev_clipped = np.clip(np.abs(np.gradient(solar_elevation_angle_2nddev)), 0.0001, 2)
ax3 = plt.subplot(gs[3])
ax3.set_title('absolute value + clipping')
ax3.plot(solar_elevation_angle_2nddev_clipped)

smoothed_signal = gaussian_filter1d(solar_elevation_angle_2nddev_clipped, 20)
ax4 = plt.subplot(gs[4])
ax4.set_title('Smoothing applied')
ax4.plot(smoothed_signal)

plt.tight_layout()
plt.show()

然后我可以使用 scipy 的 argrelmax 函数轻松识别局部最大值:

max_idx = argrelmax(smoothed_signal)[0]
print(max_idx)
# [ 689 1019 2356 2685]

它正确识别了我正在寻找的扫描线索引:

fig, ax = plt.subplots()
ax.set_title('Solar elevation angle')
ax.set_xlabel('Scanline')
ax.set_ylabel('Solar elevation angle [deg]')
ax.plot(solar_elevation_angle)
ax.scatter(max_idx, solar_elevation_angle[max_idx], marker='x', color='red')
plt.show()

我的问题是:有没有更好的方法来解决这个问题？
我发现必须手动指定限幅阈值以消除高斯滤波器中的噪声和 sigma 会大大削弱这种方法，从而无法将其应用于其他类似情况。

Answer 1

第一个改进是使用 a Savitzky-Golay filter 以噪音较小的方式找到导数。例如，它可以将抛物线(在最小二乘意义上)拟合到每个特定大小的数据切片，然后对该抛物线求二阶导数。结果比仅采用 gradient 的二阶差分要好得多。这是窗口大小 101:

savgol_filter(solar_elevation_angle, window_length=window, polyorder=2, deriv=2) 

其次，与其使用 argrelmax 寻找最大值点，不如寻找二阶导数较大的地方；例如，至少是其最大尺寸的一半。这当然会返回许多索引，但我们可以查看这些索引之间的差距，以确定每个峰值的开始和结束位置。然后很容易找到峰的中点。

这是完整的代码。唯一的参数是窗口大小，设置为 101。大小 21 或 201 给出基本相同的结果(它必须是奇数)。

from scipy.signal import savgol_filter
window = 101
der2 = savgol_filter(solar_elevation_angle, window_length=window, polyorder=2, deriv=2)
max_der2 = np.max(np.abs(der2))
large = np.where(np.abs(der2) > max_der2/2)[0]
gaps = np.diff(large) > window
begins = np.insert(large[1:][gaps], 0, large[0])
ends = np.append(large[:-1][gaps], large[-1])
changes = ((begins+ends)/2).astype(np.int)
plt.plot(solar_elevation_angle)
plt.plot(changes, solar_elevation_angle[changes], 'ro')
plt.show()

insert 和 append 的大惊小怪是因为具有大导数的第一个索引应该符合“peak begins”的条件，而最后一个这样的索引应该符合“peak ends”的条件，即使它们旁边没有合适的间隙(差距是无限的)。

分段线性拟合

这是一种替代(不一定更好)的方法，它不使用导数:拟合 smoothing spline of degree 1 (即分段线性曲线)，并注意它的节点在哪里。

首先，将数据(我称之为 y 而不是 solar_elevation_angle)归一化为标准差 1。

y /= np.std(y)

第一步是构建一条分段线性曲线，该曲线与数据的偏差最多为给定阈值，任意设置为 0.1(此处没有单位，因为 y 已标准化)。这是通过重复调用 UnivariateSpline 来完成的，从一个大的平滑参数开始并逐渐减小它直到曲线适合。 (不幸的是，不能简单地传入所需的统一误差范围)。

from scipy.interpolate import UnivariateSpline
threshold = 0.1

m = y.size
x = np.arange(m)
s = m
max_error = 1
while max_error > threshold: 
  spl = UnivariateSpline(x, y, k=1, s=s)
  interp_y = spl(x)
  max_error = np.max(np.abs(interp_y - y))
  s /= 2
knots = spl.get_knots()
values = spl(knots)

到目前为止，我们找到了节点，并记下了这些节点处的样条值。但并非所有这些结都非常重要。为了测试每个结的重要性，我将其移除并在没有它的情况下进行插值。如果新的插值与旧的有很大不同(加倍误差)，结点被认为是重要的并被添加到发现的斜率变化列表中。

ts = knots.size
idx = np.arange(ts)
changes = []
for j in range(1, ts-1):
  spl = UnivariateSpline(knots[idx != j], values[idx != j], k=1, s=0)
  if np.max(np.abs(spl(x) - interp_y)) > 2*threshold:
    changes.append(knots[j])
plt.plot(y)
plt.plot(changes, y[np.array(changes, dtype=int)], 'ro')
plt.show()

理想情况下，将分段线性函数拟合给定数据，增加结数直到再增加一个不会带来“实质性”改进。以上是 SciPy 工具的粗略近似值，但远非最佳。我不知道 Python 中有任何现成的分段线性模型选择工具。