使用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()

我的问题是:有没有更好的方法来解决这个问题？
我发现必须手动指定削波阈值以消除噪声和高斯滤波器中的西格玛,这大大削弱了这种方法,使其无法应用于其他类似的情况.

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()

插入和附加的大惊小怪是因为具有大导数的第一个指数应该被称为"峰值开始",并且最后一个这样的指数应该被称为"峰值终点",即使它们旁边没有合适的间隙(差距)是无限的).

分段线性拟合

这是一种替代(不一定更好)的方法,它不使用导数:拟合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中任何现成的分段线性模型选择工具.