检测信号的爆发和重大变化 -凯发k8网页登录

通过累积和检测爆发

在许多实际应用中，您需要监控数据，并且希望在底层过程发生变化时尽快收到警报。通常可以使用累积和 (cusum) 控制图来实现这一点。

为了说明 cusum 的工作原理，我们首先查看 2014 年西非埃博拉疫情爆发的总报告病例， 有详细记录。

load westafricanebolaoutbreak2014
plot(whoreportdate, [totalcasesguinea totalcasesliberia totalcasessierraleone],'.-')
legend('guinea','liberia','sierra leone')
title('total suspected, probable, and confirmed cases of ebola virus disease')

通过观察几内亚疫情第一次爆发的前沿，您会发现前一百个病例是在 2014 年 3 月 25 日左右报告的，之后显著增加。值得注意的是，虽然利比里亚在 3 月份也有一些疑似病例，但病例数量一直相对可控，直到大约 30 天后才急剧增加。

为了了解新发现的发病率，我们来绘制从 2015 年 3 月 25 日发病开始总病例数的相对日变化。

dayssinceoutbreak = datetime(2014, 3, 24 (0:400));
cases = interp1(whoreportdate, totalcasesliberia, dayssinceoutbreak);
dayoverdaycases = diff(cases);
plot(dayoverdaycases)
title('rate of new cases (per day) in liberia since march 25, 2014');
ylabel('change in number of reported cases per day');
xlabel('number of days from start of outbreak');

如果放大前一百天的数据，您可以看到，虽然最初涌现大量病例，但其中许多病例都发生在第 30 天之前，在第 30 天，变化率暂时降至零以下。您还会看到，在第 95 天和第 100 天之间呈现显著上升趋势，达到每天 7 个新病例的速度。

xlim([1 101])

对输入数据执行 cusum 测试可以快速确定疫情爆发的时间。cusum 跟踪两个累积和：检测局部均值何时向上移动的上和，以及检测均值何时向下移动的下和。积分方法使 cusum 能够忽略发病率中的大（瞬态）尖峰，但仍对发病率中更稳定的小变化保持敏感。

带默认参数调用 cusum 将检查前 25 个采样的数据，并在初始数据中的均值偏离超过 5 倍标准差时发出警报。

cusum(dayoverdaycases(1:101))
legend('upper sum','lower sum')

请注意，cusum 在第 30 天（第 33 天）发现误报病例，并在第 80 天（第 90 天）发现疫情的首次爆发。如果您将这些结果与之前的图进行仔细比较，可以看到 cusum 能够忽略第 29 天的假上升，但仍在第 95 天开始的大幅上升趋势之前五天触发警报。

如果您调整 cusum，使其目标均值为零例/天，目标为加减三例/天，则您可以忽略第 30 天的错误警报，并在第 92 天发现疫情爆发：

climit = 5;
mshift = 1;
tmean = 0;
tdev = 3;
cusum(dayoverdaycases(1:100),climit,mshift,tmean,tdev)

发现方差的显著变化

检测统计量突变的另一种方法是通过变化点检测，它将一个信号分成若干相邻的段，在每个段中有一个统计量（例如，均值、方差、斜率等）为常量。

下一个示例是分析公元 622 年至 1281 年间尼罗河的年最低水位，该水位在开罗附近的洛达水位表上测得。

load nilometer
years = 622:1284;
plot(years,nileriverminima)
title('yearly minimum level of nile river')
xlabel('year')
ylabel('level (m)')

公元 715 年左右，开始建造一种更新、更精确的测量设备。在此之前的信息我们知道的不多，但在进一步的检查中，您可以看到在大约 722 年后变异性要小得多。为了找到新设备开始运行的时间段，您可以在执行按元素微分以去除任何缓慢变化趋势后搜索均方根水位的最佳变化。

i = findchangepts(diff(nileriverminima),'statistic','rms');
ax = gca;
xp = [years(i) ax.xlim([2 2]) years(i)];
yp = ax.ylim([1 1 2 2]);
patch(xp,yp,[.5 .5 .5],'facealpha',0.1)

按采样微分是一种简单的去趋势方法，还有其他更复杂的方法来检验更大规模的方差。有关如何使用此数据集通过小波执行变化点检测的示例，请参阅 (wavelet toolbox)。

检测输入信号的多个变化

下一个示例是关于 cr-cr 4 速传输模块的 45 秒仿真，采样间隔为 1 毫秒。汽车发动机 rpm 和扭矩的仿真数据如下所示。

load simcarsig
subplot(2,1,2)
plot(cartorquenm)
xlabel('samples')
ylabel('torque (n m)')
title('torque')
subplot(2,1,1);
plot(carenginerpm)
xlabel('samples')
ylabel('speed (rpm)')
title('engine speed')

在此示例中，汽车加速，换挡三次，换到空档，然后刹车。

由于发动机转速天然适合建模为一系列线性段，因此您可以使用 findchangepts 来找到汽车换挡的采样。

figure
findchangepts(carenginerpm,'statistic','linear','maxnumchanges',4)
xlabel('samples')
ylabel('engine speed (rpm)')

在此处，您可以看到四个变化（在五个线性段之间），它们发生在 10000、20000、30000 和 40000 采样标记附近。放大波形的空闲部分：

xlim([18000 22000])

直线拟合密切跟踪输入波形。然而，拟合可进一步改进。

观察信号间共享的多阶段事件的变化

为了观察改进情况，请将变化点的数量增加到 20 个，并观察采样编号 19000 处换挡附近的变化

findchangepts(carenginerpm,'statistic','linear','maxnumchanges',20)
xlim([18000 22000])

请注意，发动机转速在采样 19035 处开始下降，在采样 19550 处稳定下来之前进行了 510 次采样。由于采样间隔为 1 毫秒，这是约 0.51 秒的延迟，是换挡后的典型时间量。

现在看看同一区域内发动机扭矩的变化点：

findchangepts(cartorquenm,'statistic','linear','maxnumchanges',20)
xlim([19000 20000])

请注意，在发动机转速稳定后 55 毫秒，发动机扭矩在采样 19605 处完全传递给车轴。此时间与发动机进气冲程和扭矩产生之间的延迟有关。

为了找到离合器在什么时间变为啮合状态，您可以进一步放大信号。

xlim([19000 19050])

离合器在采样 19011 处踩下，并在大约 30 个采样（毫秒）后才完全松开。