CC BY
138
НЕКОТОРЫЕ МЕТОДЫ НЕЛИНЕЙНОГО АНАЛИЗА
В.Б. Киселев
Кратко рассмотрен перспективный метод анализа процессов — метод рекуррентных диаграмм. Уделено внимание реконструкции фазового пространства из одиночного ряда измерений, построению рекуррентных диаграмм и количественному анализу, основанному на их геометрических структурах.
Введение
Развитие науки и техники ставит перед человечеством все увеличивающееся число вопросов, подлежащих всестороннему изучению. Большинство из уже стоящих перед современной наукой проблем требуют более глубокого понимания зависимостей, имеющихся в Природе. Исследования сложных систем, как природных, так и искусственных, показали, что в их основе лежат нелинейные процессы, тщательное изучение которых необходимо для обеспечения понимания и моделирования сложных систем.
В последние десятилетия набор традиционных (линейных) методик исследования был существенно расширен нелинейными методами, полученными из теории нелинейной динамики и хаоса; многие исследования были посвящены оценке нелинейных характеристик и свойств (скейлинг, фрактальная размерность) процессов, протекающих в Природе. Однако большинство методов нелинейного анализа требуют либо достаточно длинных, либо стационарных рядов данных, которые довольно трудно получить при наблюдении Природы. Более того, было показано, что данные методы дают удовлетворительные результаты для идеализированных моделей реальных систем [1]. Эти факторы потребовали разработки новых методик нелинейного анализа данных.
В последние 15 лет развился оригинальный метод, основанный на фундаментальном свойстве диссипативных динамических систем — рекуррентности (повторяемости) состояний, выражающейся в том, что, несмотря на то, что даже самое малое возмущение в сложной динамической системе может привести систему к экспоненциальному отклонению от ее состояния, через некоторое время система стремится вернуться к состоянию, некоторым образом близкому к предыдущему, и проходит при этом подобные этапы эволюции. Такое рекуррентное поведение можно изобразить рекуррентными диаграммами (recurrence plots, RP), предложенными Экманом (Eckmann) и др. [2]. Практическое использование этого метода показало его применимость к коротким и нестационарным данным, что повысило его популярность у исследователей. В 1992 году Збилут (Zbilut) и Веббер (Webber) опубликовали работу [3], в которой описали элементы количественного анализа, основанного на рекуррентных диаграммах (recurrence quantification analysis, RQA). В настоящее время количество англоязычных публикаций на тему рекуррентного анализа превысило пять сотен. Следует отметить, что рекуррентные диаграммы имеют не только сугубо исследовательское, но и эстетическое значение, привлекая внимание так же, как и фрактальные рисунки.
В 1998 и 1999 гг. независимо [4,5] были предложены кросс-рекуррентные диаграммы (cross recurrence plots, CRP), позволяющие сравнивать две рассматриваемые системы, выявляя состояния подобия.
Траектории в фазовом пространстве
Состояние природных или искусственных систем обычно изменяется во времени. Изучение этих, как правило, сложных, процессов - важная задача во многих дисциплинах, позволяющая понять и описать сущность этих процессов, например, с распространенной целью прогнозирования их состояния во времени. Целью таких исследований
является нахождение математических моделей, которые бы достаточно соответствовали реальным процессам и могли бы быть использованы для решения поставленных задач.
Состояние системы может быть описано ее переменными состояния
X1 (t), X 2 (t), ... , x d (t), (1)
где верхний индекс - номер переменной. Набор из d переменных состояния во время t составляет вектор состояния X(t) в d -мерном фазовом пространстве. Данный вектор перемещается во времени в направлении, определяемом его вектором скорости:
X (t ) = d tX(t ) = F (x). (2)
Последовательность векторов X (t) образует траекторию в фазовом пространстве, причем поле скорости F (x ) касательно этой траектории. Эволюция данной траектории описывает динамику системы и ее аттрактор. Зная F (x) , можно получить информацию о состоянии системы в момент t путем интегрирования выражения (2). Так как форма траектории позволяет судить о характере процесса - периодические или хаотические процессы имеют характерные фазовые портреты - то для определения состояния системы не обязательно проводить интегрирование, достаточно построить графическое изображение траектории (подробнее см. [6, 7]).
При исследовании сложных систем зачастую либо нет информации обо всех переменных состояния, либо не все из них возможно измерить. Как правило, имеется единственное наблюдение u(t), проведенное через дискретный временной интервал At. Таким образом, измерения записываются в виде ряда ui (t), где t = i At. Интервал
At может быть постоянным, однако это не всегда возможно и создает проблемы для применения стандартных методов анализа данных, требующих равномерной шкалы наблюдений.
Взаимодействия и их количество в сложных системах таковы, что даже по одной переменной состояния можно судить о динамике всей системы в целом (данный факт был установлен группой американских ученых при изучении турбулентности). Таким образом, эквивалентная фазовая траектория, сохраняющая структуры оригинальной фазовой траектории, может быть восстановлена из одного наблюдения или временного ряда по теореме Такенса (Takens) методом временных задержек:
X(t)=(u,, Uг+т, Uг+{m-1)т), (3)
где m - размерность вложения, т - временная задержка (реальная временная задержка определяется как т-At). Топологические структуры восстановленной траектории сохраняются, если m > 2 • d +1, где d - размерность аттрактора [8]. Оба параметра вложения - размерность m и задержка т - выбираются априорно.
Существует несколько подходов к выбору минимально достаточной размерности m ; высокую эффективность показали методы, основанные на концепции ложных ближайших точек (false nearest neighbours, FNN). Суть ее заключается в том, что при уменьшении размерности вложения происходит увеличение количества ложных точек, попадающих в окрестность любой точки фазового пространства. Вытекающий отсюда простейший метод - определение количества FNN как функции от размерности, в пределе равной 0. Существуют и другие методы, основанные на этой концепции - например, определение отношений расстояний между одними и теми же соседними точками при разных m .
Так как низкая точность измерений и случайные ошибки могут привести к появлению линейной зависимости между последовательными векторами X(t), то задержку т следует выбирать так, чтобы минимизировать эту зависимость.
В последние годы были разработаны эвристические методы определения лучших размерности и задержки, с использованием RP. Вопросы определения размерности вложения также обсуждаются в работе [9]
Рекуррентные диаграммы (recurrence plots, RP)
Процессы в Природе обладают ярко выраженным рекуррентным поведением, таким, как периодичность или иррегулярная цикличность. Более того, рекуррентность (повторяемость) состояний в смысле прохождения последующей траектории достаточно близко к предыдущей является фундаментальным свойством детерминированных динамических систем. Это свойство, отмеченное еще в 80-х годах XIX века французским математиком Пуанкаре (Poincare), впоследствии было сформулировано в виде «теоремы рекуррентности», опубликованной в 1890 году:
Если система сводит свою динамику к ограниченному подмножеству фазового пространства, то система почти наверняка, т.е. с вероятностью, практически равной 1, сколь угодно близко возвращается к какому-либо изначально заданному режиму.
Как правило, размерность фазового пространства сложной динамической системы более 3, что делает практически неудобным его рассмотрение напрямую (при помощи проекций в 2- и 3-мерные пространства). В 1987 году Экман и др. предложили [2] способ отображения m -мерной фазовой траектории состояний x(t) на двумерную квадратную двоичную матрицу размером N х N, в которой 1 (черная точка) соответствует повторению состояния при некотором времени i в некоторое другое время j, а обе координатные оси являются осями времени. Такое представление было названо рекуррентной диаграммой (recurrence plot, RP). Математически вышесказанное описывается как
Rj = 0S-||Xi - Xj|) x е Rm, i, j = 1k N, (4)
где N - количество рассматриваемых состояний xf, si - размер окрестности точки X в момент i, 11 • || - расстояние и 0( • ) - функция Хэвисайда.
Поскольку Rи = 1 (( = 1 к N) по определению, то RP всегда содержит черную диагональную линию - линию идентичности (line of identity, LOI), под углом п/4 к осям координат. Произвольно взятая рекуррентная точка (i, j) не несет какой-либо информации о состояниях во времена i и j. Только вся совокупность рекуррентных точек позволяет восстановить свойства системы.
Непрактично и, как правило, невозможно обнаружить полную рекуррентность в смысле X = Xj (состояние хаотической системы не повторяется полностью эквивалентно начальному состоянию, а подходит к нему сколь угодно близко). Таким образом, рекуррентность определяется как достаточная близость состояния x j состоянию
X . Иными словами, рекуррентными являются состояния Xj, попадающие в m -мерную
окрестность с радиусом s и центром в X . Эти точки Xj называются рекуррентными
точками (recurrence points) - см. рис. 1.а, на котором показан увеличенный сегмент аттрактора Лоренца. В выражении (4) сказанное описывается функцией Хэвисайда 0 и ее параметром si .
Существует два подхода к определению параметра si, выбор которого зависит от целей и условий применения:
1. s =s, Vi; Rt j = R ji ;
2. s * s ; Ru * R}1.
Первый является наиболее используемым и приводит к симметричной RP. Второй означает, что для каждого состояния радиус определяется индивидуально, что приводит к несимметричной RP, но позволяет получить одинаковую плотность рекуррентности каждого столбца RP. В этом случае параметр используется для задания плотности рекуррентности — например, при s = 0,2, выбранном для момента i, радиус
si покроет 20% всех векторов фазового пространства. Такой подход называется подходом с постоянным числом соседних точек (fixed amount of nearest neighbours, FAN).
20 -40
а
sco
400
8L.....58.
к, %
К ЙЩ
100 200 ЗШ 400 ьоо
Рис. 1. а - сегмент фазового портрета системы Лоренца (ряд получен при г = 28, & = 10, ь = 8/3); б - рекуррентная диаграмма (т = 1, е = 0,15 , норма Ьт); в - диаграмма расстояний, на которой строится расстояние до ближайшей рекуррентной точки
Рис. 2. Основные нормы с одинаковым радиусом вокруг точки фазового пространства (черная точка), изображенные для двумерного фазового пространства: а - ь1, б - Ь2,
в - L
б
а
в
Перед расчетом КР также следует уделить внимание выбору нормы. Наиболее известные в данном случае нормы - Ц, Ь2 (евклидова норма), (максимальная норма). Границы этих норм имеют разные фигуры (см. рис. 2), что применительно к подходу с постоянным е означает, что норма обеспечивает нахождение наибольшего, Ь2 -среднего, а Ь1 - минимального количества соседних точек. Как правило, при расчете КР используется норма , так как, во-первых, она независима от размерности фазового пространства, во-вторых, наиболее просто вычисляема и, в-третьих, позволяет изу-
чать ЯР теоретически, поскольку аналитические выражения для решаются проще, нежели для любой другой нормы.
Независимость от размерности находит применение при анализе временных рядов с различными размерностями вложения, т.е. ЯР разных вложений могут сравниваться напрямую, в то время как для остальных норм требуется масштабирование.
Существует некоторое количество разновидностей ЯР. • ЯР с областью соседства в виде коридора \апп, £оШ ] (такие диаграммы неприменимы для количественного анализа) [9]:
Т> Lsin ,sout. R i, j
= ©(|X - X\-Sn Ms,
- V* - V
out pi Л j
Вместо изображения рекуррентности в виде черных точек можно изображать расстояния между состояниями на диаграмме расстояний (distance plot) [9, 10], которые отображаются на некоторую цветовую палитру (см. рис. 1.в):
Dm■=1X - j.
Оконная рекуррентная диаграмма (windowed RP) была предложена как средство изучения воздействия внешней силы или нестационарности в системе [11, 12]. Данная диаграмма получается путем покрытия исходной RP квадратами размером w х w (окнами) и усреднения значений рекуррентных точек внутри этих окон [11]. Таким образом, оконная RP представляет собой Nw х Nw -матрицу, размер которой Nw вычисляется как округленное в большую сторону отношение N/w; данная матрица состоит из значений, не ограниченных 0 и 1 и может быть изображена цветными точками. Эти значения соответствуют кросскорреляционной сумме
1
cms =-т Z R
w
m,s
i+( I -1) w, j+(J-1) w'
i, j=1
I, J = 1k N
w
между секциями длиной w и начинающимися в (I - 1)w +1 и (J - 1)w +1; • Разновидности RP не ограничиваются приведенными. Более обширный их список представлен, например, в работе [13].
Изначальное предназначение RP - зрительный анализ траекторий в фазовых пространствах высоких размерностей; внешний вид диаграммы может дать представление об эволюции этих траекторий во времени. Можно выделить два основных класса структур RP: топология (typology), соответствующая крупномасштабным структурам, и текстура (texture), соответствующая мелкомасштабным структурам.
ттм
'f'l.
г;:;.;
а
llffffi : М -IfЩ •.
эьшШШШ' 'Г'""
б
в
г
Рис. 3. Характерные топологии РР: а - однородная (нормально распределенный шум); б - периодическая (генератор Ван дер Поля); в - дрейф (отображение Икеды с наложенной линейно растущей последовательностью); г - контрастные области или полосы (броуновское движение)
Топология дает общее представление о характере процесса; выделяют четыре различных класса:
• однородные RP типичны для стационарных и автономных систем, в которых время релаксации мало по сравнению с длиной ряда (рис. 3.а);
• периодические повторяющиеся структуры (диагональные линии, узоры в шахматном порядке) соответствуют различным осциллирующим системам (рис. 3.б);
• дрейф соответствует системам с медленно изменяющимися параметрами, что делает белыми левый верхний и правый нижний углы RP (рис. 3.в);
• резкие изменения в динамике системы, равно как и экстремальные ситуации, обуславливают появление белых областей или полос (рис. 3.г). RP упрощает выявление экстремальных и редких событий.
Пристальное изучение RP позволяет выявить мелкомасштабные структуры - текстуру, которая составляется из простых точек, диагональных, горизонтальных и вертикальных линий. Комбинации вертикальных и горизонтальных линий формируют прямоугольные кластеры рекуррентных точек.
• Одинокие, отдельно стоящие рекуррентные точки появляются в том случае, когда соответствующие состояния редки, либо неустойчивы во времени, либо вызваны сильной флуктуацией. При этом они не являются признаками случайности или шума.
• Диагональные линии Ri+k j+k = 1 (при k = 1.../, где l - длина диагональной линии)
появляются в случае, когда сегмент траектории пролегает параллельно другому сегменту, т. е. траектория повторяет саму себя, возвращаясь в одну и ту же область фазового пространства в разное время. Длина таких линий определяется временем, в течение которого сегменты траектории остаются параллельными; направление (угол наклона) линий характеризует внутренне время подпроцессов, соответствующих данным сегментам траектории. Прохождение линий параллельно LOI (под углом ж/ 4 к осям координат) свидетельствует об одинаковом направлении сегментов траектории, перпендикулярность к LOI - о противоположном («отраженные» сегменты), что может также являться признаком несоответствующего вложения. Нерегулярное появление диагональных линий является признаком хаотического процесса.
• Вертикальные (горизонтальные) линии Ri j+k = 1 (при k = 1. v, где v - длина вертикальной линии) выделяют промежутки времени, в которые состояние не изменяется или изменяется незначительно (система как бы «заморожена» на это время), что является типичным поведением при ламинарных состояниях.
Рис. 3 показывает, насколько различными могут быть мелкомасштабные структуры. Большое количество отдельных точек и небольшое количество линий - признаки сильных флуктуаций, которые характерны для некоррелированного шума (рис. 3.а). Очевидна зависимость между периодическими повторяющимися структурами и осцилляторами (рис. 3.б). Строгая рекуррентная динамика вызывает появление длинных диагоналей на одинаковом расстоянии друг от друга. Нерегулярное появление коротких и длинных диагональных линий характерно для хаотических процессов (рис. 3.в); при этом нерегулярное появление обширных черных кластеров и обширных белых областей соответствует нерегулярному поведению системы (которая может быть, например, коррелированным шумом, рис. 3.г).
Перечисленные мелкомасштабные структуры лежат в основе количественного анализа рекуррентных диаграмм. Несмотря на достаточно точную классификацию, на практике сигналы отягощены шумами, нередко - весьма значительными, в результате чего анализ RP требует некоторого опыта.
Количественный анализ рекуррентных диаграмм (recurrence quantification analysis, RQA)
Очевидно, что образуемые структуры рекуррентных диаграмм можно некоторым образом анализировать численно. Впервые элементы количественного анализа рекуррентных диаграмм были предложены в [2] - авторами был разработан инструментарий вычисления мер сложности структур RP, используя плотность рекуррентных точек и диагональных структур. Ряд исследований показал, что эти меры могут использоваться для выявления точек бифуркации и переходов хаос-порядок.
Позднее были предложены меры, использующие плотность вертикальных структур, позволяющие выявлять переходы хаос-хаос и нечувствительные к стационарности и длине исследуемых рядов.
Для вычисления мер используются рекуррентные диаграммы, построенные с фиксированным значением порога s .
Мера рекуррентность (recurrence rate, RR) 1 N
RR = -V У RT (5)
N2 j г, 1
iy z,j =1
показывает плотность рекуррентных точек, просто подсчитывая их, включая LOI. В пределе
1 N
P%= lim-V У Rm,s
• N2 j ^1
ly z, 1=1
данная мера показывает вероятность нахождения рекуррентной точки в RP (вероятность повторения состояния).
Следующая мера рассматривает диагональные линии. Частотное распределение длин l диагональных линий в RP Ps(l)={l;i = 1...Nl}, где Nl - абсолютное количество диагональных линий (каждая линия считается только один раз). Процессы со стохастическим поведением могут порождать очень короткие диагонали либо вообще не порождать их, в то время как детерминистские процессы дают длинные диагонали и малое количество отдельных рекуррентных точек. Таким образом, отношение рекуррентных точек, составляющих диагональные структуры, к общему количеству рекуррентных точек
yw lPs(l)
DET = ^Mss-li (6)
z, 1 R z,1
называется мерой детерминизма (determinism, DET) или предсказуемости системы. Следует отметить, что эта мера не имеет значения реального детерминизма процесса. Пороговое значение lmin исключает диагональные линии, образованные тангенциальным движением траектории в фазовом пространстве. При lmin = 1 DET = RR .
Диагональные структуры показывают время, в течение которого участок траектории подходит достаточно близко к другому участку траектории. Таким образом, эти линии позволяют судить о расхождении элементов траектории. Средняя длина диагональных линий
У ipsIi)
L = ^Nh---(7)
yl=i ps()
^^^ min
- это среднее время, в течение которого два участка траектории проходят близко один к другому, и может рассматриваться как среднее время предсказуемости. В RQA, однако, используется максимальная длина диагональных структур, либо ее инверсия - дивергенция (divergence, DIV):
Lmax = max(; i = 1...N,}), DIV = (8)
max
Установлено, что длины диагональных линий соотносятся с наибольшим положительным показателем Ляпунова, если он существует для рассматриваемой системы [2]. Разными авторами были предложены методы оценки максимального положительного показателя Ляпунова с использованием длин диагональных линий.
Мера энтропии (entropy, ENTR) соотносится с энтропией Шеннона (Shannon) частотного распределения длин диагональных линий
N PS(] )
ENTR = -£p(()lnp(), где p(()= If ^ , (9)
'^ P ()
и отражает сложность детерминистской составляющей в системе.
Тренд (trend, TREND) - мера, представляющая линейный коэффициент убывания частной плотности рекуррентных точек RR„ диагоналей, параллельных к LOI, как функцию времени (расстояния) между этими диагоналями и LOI
YN(i - N/2I -( жг) )
TREND = —^A—Ц^—'Л. (10)
Eli - n/2)2 ( )
Тренд характеризует нестационарность процесса, особенно - дрейф. Выражение исключает границы RP (N < N ). Выбор значения N зависит от изучаемого процесса.
Рассмотрим точку траектории xi и множество ассоциированных с ней рекуррентных точек R' = |x;. : Rt j = 1; j g [l...N]}. Определим подмножество этих рекуррентных точек r = & g R' : (ra • Rhk + ) + (r,k • Ra-1 ) > 0; l g [1,..., NJ R,0 = R^+1 = 0}, которое содержит рекуррентные точки, формирующие вертикальные линии на RP в колонке i. Определим теперь длины v всех соединенных в si подмножеств {x £ si; xj+1,...,xj+v g si; xj+v+1 £ si}. Обозначим p(v) = {vk; k = 1,2,...,K} множество длин
присутствующих в si подмножеств, и из выражения U=1 P' (v) определяем распределение длин вертикальных линий Ps(v ) для всего RP в целом. Мера замирания (laminarity, LAM) TN vPs(v )
LAM = -— (11)
',j R'J
определяется отношением количества рекуррентных точек, образующих вертикальные линии, к общему количеству рекуррентных точек. LAM характеризует наличие состояний замирания системы (т. е. когда движение системы по фазовой траектории останавливается или продвигается очень медленно). Средняя длина вертикальных структур
TN vPs{v)
TT = -— (12)
ZN=v. PS(v)
y vmin
называется мерой времени остановки (trapping time, TT) и характеризует среднее время, которое система может провести в определенном состоянии. Мера максимальной длины вертикальных структур
Vmax = max(k; k = 1.K}) (13)
аналогична мере максимальной длины диагональных структур.
Подробнее меры, основанные на анализе вертикальных структур, рассмотрены в
Заключение
В данной статье была вкратце рассмотрена методика рекуррентных диаграмм (recurrence plots, RP) и количественного анализа на основе их геометрических структур. Не было уделено внимание расширению RP, носящему название кросс-рекуррентных диаграмм (cross-recurrence plots, CRP) [4,5]. Если RP строятся путем оценки расстояний между точками траектории в фазовом пространстве, то при построении CRP оцениваются расстояния между точками X и y двух разных траекторий в одном фазовом пространстве. Аналогично, количественный анализ CRP призван определить меру подобия двух систем.
В последние годы появилось множество работ, посвященных, помимо теоретического изучения RP, еще и вопросам их практического применения. RP активно применяются в геофизике, палеонтологии, медицине (в частности, для оценки сердечной деятельности с целью выявления болезней и предсказания их возможного развития), экономике. Преимущества, предоставляемые данным методом исследователям - одинаковая работа с данными любой размерности, некритичность к источнику и качеству данных, маргинальность (метод не вносит дополнительную информацию «от себя»), наглядность результатов - позволяют применять его достаточно широко при анализе и изучении сложных динамических систем совершенно различной природы.
Следует отметить, что RP сами по себе предоставляют поле для теоретического изучения как их самих, так и получаемых на их основе количественных показателей, так как удовлетворительная теория использования этих показателей при анализе еще не разработана.
Литература
1. Manuca R., Savit R. Stationarity and nonstationarity in time series analysis. // Physica D, 99(2-3), 134-161 (1996).
2. Eckmann J.-P., Kamphorst S. O., Ruelle D. Recurrence Plots of Dynamical Systems. // Europhysics Letters 5, 1987, 973-977.
3. Zbilut J. P., Webber Jr. C. L. Embeddings and delays as derived from quantification of recurrence plots. // Physics Letters A 171 (3-4), 1992, 199-203.
4. Zbilut J. P., Giuliani A., Webber Jr. C. L. Detecting deterministic signals in exceptionally noisy environments using cross-recurrence quantification. // Physics Letters A 246 (1-2), 1998, 122-128.
5. Marwan N. Untersuchung der Klimavariabilitat in NW Argentinien mit Hilfe der quantitativen Analyse von Recurrence Plots. Masters Thesis, Dresden University of Technology, 1999.
6. Берже П., Помо И., Видаль К. Порядок в хаосе. О детерминированном подходе к турбулентности. М.: Мир, 1991.
7. Мун Ф. Хаотические колебания: Вводный курс для научных работников и инженеров. М.: Мир, 1990.
8. Takens F. Detecting Strange Attractors in Turbulence. Р. 366-381, Springer, Berlin, 1981.
9. Iwanski J. S., Bradley E. Recurrence plots of experimental data: To embed or not to embed? // Chaos, 8(4), 861-871 (1998).
10. Webber Jr. C. L. Recurrence Quantification Analysis, 2003. http://homepages.luc.edu/~cwebber
11. Casdagli M. C. Recurrence plots revisited. // Physica D, 108(1-2), 12-44 (1997).
12. Manuca R., Savit R. Stationarity and nonstationarity in time series analysis. // Physica D, 99(2-3), 134-161 (1996).
13. Marwan N. Encounters With Neighbours - Current Developments Of Concepts Based On Recurrence Plots And Their Applications, PhD Thesis, University of Potsdam, ISBN: 300-012347-4 (2003). http://pub.ub.uni-potsdam.de/2003/0026/marwan.pdf
14. Marwan N., Wessel N., Meyerfeldt U., Schirdewan A., Kurths J. Recurrence-plot-based measures of complexity and their application to heart-rate-variability data. // Physical Review, E66, 026702 (2002).
