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

Question

我正在尝试用 Python 找出在一系列计划的航天器机动中发生的俯仰旋转的通用方法。你可以将其视为 移位检测 问题的一个特例。

让我们考虑我的测量数据中的 solar_elevation_angle 变量，它标识了从航天器仪器测量到的太阳高度角。对于那些想要使用数据的人，我将 solar_elevation_angle.txt 文件保存在 这里。

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，会严重削弱这种方法，使其无法应用于其他类似情况。

- stm4tt

Answer 1

首先的改进是使用Savitzky-Golay滤波器以更少的噪声方式找到导数。例如，它可以对某个大小的每个数据片段拟合一个抛物线（在最小二乘意义下），然后取该抛物线的二阶导数。结果比仅使用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()

插入和追加的问题是因为具有大导数的第一个索引应该被视为“峰值开始”，而最后一个这样的索引应该被视为“峰值结束”，即使它们旁边没有合适的间隙（间隙是无限的）。

分段线性拟合

这是一种替代方法（不一定更好），它不使用导数：拟合一阶平滑样条（即分段线性曲线），并注意其节点的位置。

首先，将数据标准化（我称之为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中是否有任何现成的分段线性模型选择工具。

- user6655984

---

谢谢，这是个很好的解决方案！我之前不太熟悉Savitzky-Golay滤波器。结果证明它非常适用于平滑处理数据。我喜欢你能够识别每次音高转换的开始和结束，这将非常有用。我需要更多时间来完全理解你的另一种方法，看起来也非常有趣。 - stm4tt

我没有尝试过但至少应该提及的方法是将 find_peaks_cwt 应用于二阶导数；它与简单的 argrelmax 不同，具有滤波功能。 - user6655984