Нелинейная динамика и анализ временных рядов – обзор алгоритма Recurrence plots

В этой статье я сделаю обзор относительно нового и довольно мощного метода нелинейной динамики – метода Recurrence plots или рекуррентного анализа в приложении к анализу временных рядов. А, кроме того, тут вы найдете код короткой программы на языке Matlab, которая реализует все нижеописанное. Итак, начнем. По долгу службы я занимаюсь нелинейной динамикой, обработкой видео и изображений, я бы даже сказал, довольно узкой частью нелинейной динамики – нелинейными колебаниями роторов. Как известно, вибросигнал представляет собой ничто иное, как временной ряд, где в качестве сигнала выступает значение амплитуды отклонения, ну например, ротора турбины самолета. Как известно, не только колебания ротора можно представить в таком виде. Колебания биржевых котировок, активность Солнца и множество других процессов описываются простым вектором чисел, выстроенным по времени. Скажу даже больше, все эти процессы объединяет один важный фактор – они нелинейны, а некоторые даже хаотичны, что означает на практике невозможность предсказать состояние в системе на сколь угодно большой отрезок времени даже зная точно закон ее движения в виде дифференциальных уравнений. А самое главное, в большинстве случаев мы не можем даже записать эти самые уравнения в каком-либо виде. И тут на помощь приходит эксперимент и нелинейная динамика.

Нелинейная динамика

Снимая показания с датчика (скачивая файл с котировкой валют), мы имеем на выходе одномерный сигнал, как правило, сложной формы. Хорошо еще если сигнал периодический. А если нет? В противном случае мы имеем дело со сложной системой, которая к тому же, может находиться в режиме нелинейных и хаотических колебаний. Под сложной системой в данном случае понимают систему с большим числом степеней свободы. Но мы-то получаем с датчика одномерный сигнал. Как узнать о других степенях свободы этой системы? Прежде чем мы переедем к главному, стоит упомянуть о еще одном мощном методе анализа нелинейной динамики, который плавно перекочевал в методы анализа временных рядов. Это так называемое восстановление траектории системы по сигналу от всего лишь одной степени свободы. К примеру, вращается по кругу шарик на нити, и если мы смотрим в одной плоскости на это процесс, то сможем снимать сигнал лишь с по одной оси, для нас шарик будет бегать туда-сюда. И вот, выяснилось, что можно полностью восстановить N – мерный сигнал всей системы. Но не в полном объеме, а лишь с сохранением топологических свойств системы (проще говоря ее геометрии). В данном случае, если шарик движется по периодической траектории, то и восстановленная траектория будет периодической. Сам метод основан на простой формуле (в случае двумерной системы): Y(i) = X(t), Y1(i) = X(t+n), где n – это временная задержка, а Y, Y1 это уже восстановленный сигнал. Доказывается это все в теореме Такенса. В примере пространство двумерно, а в реальности оно может представлять из себя пространство сколь угодно большего числа измерений. Метод оценки размерности пространства не входит в тему данного топика, упомяну лишь о том, что это может быть, к примеру, метод ложных соседей.

Метод Recurrence plots

Итак, мы получили траекторию системы, в которой сохранены топологические свойства реальной физической системы, что очень важно. Теперь мы можем натравить на нее целый арсенал методов распознавания образов для выявления закономерностей. Но все они так или иначе имеют свои недостатки вычислительные или просто сложны в реализации. И вот в 1987 году Экман с коллегами разработали новый метод, суть которого сводиться к следующему. Полученная выше траектория, представляющая собой набор векторов размерности N отображается на двумерную плоскость по формулам: R[i,j,m,e] = O(e[i], — ||S[i] – S[j] ||), где i,j – индексы точки на плоскости, m – размерность пространства вложения, O – функция Хевисайда ||…|| — норма или расстояние (например евклидово). Внешний вид полученной и визуализированной матрицы и даст нам представление о динамике системы, изначально представленной в виде временного ряда. Все вышесказанное проиллюстрируем в виде эксперимента с числами Вольфа, которые описывают активность Солнца (вместо этого ряда легко можно взять котировки валют).

Результаты

Это, собственно, сам сигнал:

Это восстановленная траектория в двумерном пространстве (мы берем двумерное пространство для простоты, хотя на самом оно содержит большее число измерений, но как показывает практика этого вполне достаточно для качественной оценки) Это так называемая матрица расстояний (по смыслу всего лишь расстояния от i-ой точки до j-ой). Часто эти диаграммы имеют замысловатый рисунок, что может быть использовано (и используется) в дизайне.

Ну вот, собственно, и сама рекуррентная диаграмма. Для описания всех методов анализа такой диаграммы не хватит и 400 страниц. Помимо качественного анализа, данный метод допускает и количественные показатели, что может с успехом использоваться в случае использования нейронных сетей. Ну а самое важное, это то, что уже бросив беглый взгляд на эту диаграмму, мы можем сказать о ней многое. Во-первых наличие полос, которые перпендикулярны главной диагонали говорит о наличии хаотического или стохастического процессов в системе (что бы сказать точно необходимы дополнительные исследования). Наличие неравномерно заполненных черными точками зон говорит о нестационарных процессах и позволяет точно определить границы этих процессов во времени.

Код в Matlab

clear; clc; %Построение рекуррентной диаграммы

[X1,X2,X3,X4]=textread(‘data.txt’,’%f %f %f %f’);

plot(X2,X4);grid;hold; ; ylabel(‘Wolf’);xlabel(‘Time’);

N = length(X1);

M = round(0.3*N);

M1 = N – M;

m = 2; %Размерность пространства вложения

t = 10; %Временная задержка

%Сюда будем писать восстановленный аттрактор

X(1,1) = 0; X(1,2) = 0;

j = 1;

%Здесь восстанавливаем траекторию

for i=M1:(N – t)

X(j,1) = X3(i); X(j,2) = X3(i + t); j = j + 1;

end

figure;

plot(X(:,1),X(:,2));grid;hold; ; ylabel(‘Wolf’);xlabel(‘Time’);

N1 = length(X(:,1));

D1(1,1) = 0; D2(1,1) = 0; e = 30;

%Строим RP – диаграмму

for i = 1:N1

for j = 1:N1

D1(i,j)=sqrt((X(i,1)-X(j,1))^2+(X(i,2)-X(j,2))^2);

if D1(i,j) < e D2(i,j) = 0;

else D2(i,j) = 1;

end;

end;

end;

figure; pcolor(D2) ; shading interp; colormap(pink); figure; pcolor(D1) ; shading interp; colormap(pink); hold on

PS. Это всего лишь элементарный пример, который, естественно, не охватывает всю мощь этого метода. В простейшем случае для автоматизации процесса вы можете использовать алгоритмы поиска фигур на плоскости. Такие системы могут быть актуальны в системах мониторинга турбин самолета, а так же системах мониторинга различных финансовых процессов. Я не затронул самую интересную часть данной тематики — вычисление количественных показателей этих диаграмм, что постараюсь компенсировать в ближайших топиках.

Литература

У нас есть похожие новости по этим темам:
Наверх