CC BY
163
ОПРЕДЕЛЕНИЕ СТАБИЛЬНОСТИ ТРАЕКТОРИИ ПРОЦЕССА В ФАЗОВОМ ПРОСТРАНСТВЕ ПРИ ПОМОЩИ РЕКУРРЕНТНОГО
АНАЛИЗА
В.Б. Киселев
Научный руководитель - к. т.н., доцент Б.А. Крылов
Предложен способ количественной оценки стабильности фазовой траектории изучаемого процесса при помощи диаграмм расстояний (подвид рекуррентных диаграмм) и их количественного анализа. Приведены примеры применения способа к модельным и реальным системам, дана оценка возможностей использования.
Введение
Исследования сложных систем, как природных, так и искусственных, показали, что в их основе лежат нелинейные процессы, тщательное изучение которых необходимо для понимания и моделирования сложных систем. В последние десятилетия набор традиционных (линейных) методик исследования был существенно расширен нелинейными методами, полученными из теории нелинейной динамики и хаоса. Однако большинство методов нелинейного анализа требуют либо достаточно длинных, либо стационарных рядов данных, которые довольно трудно получить из природы. Более того, было показано, что данные методы дают удовлетворительные результаты в основном для идеализированных моделей реальных систем [1].
Рекуррентный анализ [2-4] - достаточно молодой и динамично развивающийся подход к анализу сложных систем, не требующий длинных или стационарных временных рядов. Рекуррентные диаграммы позволяют судить о характере протекающих в системах процессов, наличии и влиянии шума, дрейфа, наличии состояний повторения и замирания (ламинарность), совершении экстремальных событий, наличии скрытой периодичности и цикличности. Количественный анализ рекуррентных диаграмм позволяет сопоставить диаграмме некоторые численные меры, основанные на плотности рекуррентных точек, диагональных и вертикальных (горизонтальных) линий. Следует отметить, что пока не создано удовлетворительной теории применения рекуррентных диаграмм и их количественных мер; этот метод сам по себе представляет собой поле для исследований.
В настоящей работе предложен способ количественной оценки стабильности фазовой траектории изучаемого процесса при помощи диаграмм расстояний (подвид рекуррентных диаграмм) и их количественного анализа. Приведены примеры применения способа к модельным и реальным системам, дана оценка возможностей использования.
Рекуррентные диаграммы
Рекуррентные диаграммы были предложены для отображения траектории х е Ят (г = 1...N) в т -мерном фазовом пространстве на двумерную двоичную матрицу размером N х N. Единица в ячейке матрицы соответствует повторению состояния (проход траектории через одну и ту же точку фазового пространства) при некотором времени г в некоторое другое время ], а обе координатные оси диаграммы являются осями времени. Математически это выражается следующим образом:
Кт/' = -
э( - хг -х] ) х еШт, г,] = 1.N, (1)
где N - количество рассматриваемых состояний х^, вг- - размер окрестности точки х в момент г, II ■ II - норма (расстояние), ©( • ) - функция Хэвисайда. Графически рекур-
рентная диаграмма может быть представлена монохромным изображением, где единице соответствует черная точка.
При наличии только одномерного ряда ut эквивалентная траектория в m -мерном фазовом пространстве может быть восстановлена по методу временных задержек Такен-са [5] Xc(() = (uj, Uj+т,..., Uj+(m-i)), где m - размерность вложения, т - временная задержка (реальная временная задержка определяется как т • At ). Топологические структуры восстановленной траектории сохраняются, если m < 2d +1, где d - размерность аттрактора. На практике оказывается, что в большинстве случаев аттрактор может быть восстановлен и при m < 2d [6]. Задержка т, как правило, выбираются априорно.
Размер окрестности Sj определяет радиус окрестности в фазовом пространстве с центром в точке Xj. Если точка Xj попадает внутрь данной окрестности, то такое состояние считается подобным состоянию Xj, и на диаграмме устанавливается Rj j = 1. Радиус Sj может быть постоянным для всех Xj либо определяться для каждой точки индивидуально, чтобы в получаемую окрестность всегда попадало определенное количество подобных состояний. В этой работе используется постоянное значение Sj, что приводит к получению симметричной рекуррентной диаграммы относительно линии Rj j = 1 ( = j). Для вычисления расстояний между точками фазового пространства
используется максимальная норма Lœ = sup|vk\.
Рисунок рекуррентной диаграммы отображает поведение процесса во времени и позволяет сделать выводы о его характере по ее топологии и текстуре.
Збилут (Zbilut) и Вебер (Webber) разработали количественный анализ рекуррентных диаграмм (recurrence quantification analysis, RQA) [7] для определения численных показателей рекуррентной диаграммы. Они предложили меры, использующие плотность рекуррентных точек и диагональные структуры диаграммы: показатель подобия (RR), детерминизм (DET), максимальная длина диагональных линий (L), энтропия (ENTR), тренд (TREND).
Рис. 1. Фрагменты диаграмм расстояний: а - рассматриваемая в работе модельная система, б - две компоненты уравнения Ван-дер-Поля для стандартных параметров
Если исключить проверку попадания точек в s -окрестность других точек, получим диаграмму расстояний, каждая точка которой несет информацию о расстоянии между точками Xj и Xj траектории [8, 9]:
D?j =||Xj -Xj||. (2)
Такая диаграмма удобна тем, что позволяет проводить исследования для различных значений s. Графически она может быть изображена при помощи отображения на некоторую цветовую палитру, например уровней серого (рис. 1). Также диаграммы расстояния могут быть интересны с эстетической точки зрения, так как сложнопериодиче-ские процессы позволяют получить довольно интересные изображения.
Количественный анализ диаграмм. Показатель подобия
Существует несколько общеизвестных мер количественного анализа рекуррентных диаграмм, основанных на подсчете плотности точек, диагональных и вертикальных (горизонтальных) линий, позволяющие численно выразить структуры на рекуррентной диаграмме. Удовлетворительная теория использования этих мер пока не разработана. Одна из этих мер - показатель подобия (recurrence rate) [7] - представляет собой меру плотности рекуррентных точек на диаграмме: 1 N
RR = Л IRy. (3)
N ',j=1
Очевидно, что для одной и той же траектории показатель подобия является функцией от радиуса окрестности s :
rr(s) = RRfa), sk = so +k'As' sk <max(Dj,j* (4)
k = 0,1,..., i, j = Ok N -1,
где so - начальное значение радиуса окрестности, As - шаг увеличения, Dj,j - диаграмма расстояний (2) размером N х N. Значение sk изменяется с шагом As в пределах от so до максимального расстояния между точками траектории в фазовом пространстве. Скорость изменения RR(s) равна
RR'(s) = RR(sk)-RR(sk-1), k = 1.K -1, (5)
где K - количество значений sk, использованное для вычисления RR. Очевидно, что при ЯЯ0. к - 2 ~ RR ~ const происходит строго линейное нарастание RR(s), что говорит о равномерном распределении точек временного ряда по траектории (RR(s) при каждом приращении sk увеличивается на RR' = const ).
Строго говоря, такая ситуация возможна, во-первых, если процесс обладает гладкой, без флуктуаций траекторией, сходящейся к точке, и, во-вторых, если задержка временного ряда At ^ 0 .
Периодические процессы, встречающиеся на практике, порождают траектории различной степени сложности, зачастую отягощенные шумами, случайными флуктуа-циями и разного рода искажениями. Если принять, что траектория гладкая, задержка временного ряда At ^ 0, то кривая роста RR(s) будет представлять собой несколько сопряженных прямых отрезков, количество которых будет зависеть от формы траектории. Шумы и случайные флуктуации системы приводят к увеличению скорости возрастания RR(s) за счет появления на рекуррентной диаграмме отдельно стоящих рекуррентных точек, вертикальных и горизонтальных линий так, что контуры рисунка диаграммы размываются.
Таким образом, по изменению размаха скорости ARR' = max(RR')- min(RR') = STAB можно судить о стабильности исследуемой траектории.
Моделирование
Проследим динамику изменения ЯЯ(г) и для траектории системы
Г х = соб(( )-г , ч
\ ){ , г = |-В1П((-п + п-|-RND, (7)
[ У = 81П(()• г
где п - уровень равномерно распределенного шума ЯКО, вносящего нестабильность в траекторию, | = 0.2 - коэффициент толщины окружности.
Фазовые портреты системы (7) для значений п = 0.0, п = 0.2 и п = 0.7 представлены на рис. 2. Временной ряд сгенерирован с временной задержкой At = 0.061 длиной порядка 3000 точек.
а б в
Рис. 2. Фазовые портреты: а - п = 0.0, б - п = 0.2 , в - п = 0.7
Система (7) порождает довольно сложную периодическую траекторию; видны искажения (рис. 2б, в), вносимые равномерно распределенным шумом. На рис. 1а показан фрагмент диаграммы расстояний для системы (7) при п = 0.0 .
Рассмотрим фрагменты рекуррентных диаграмм для системы (7) при п = 0.0 (рис. 3) и п = 0.7 (рис. 4) для разных s . Из диаграмм видно, что плотность рекуррентных точек в случае зашумленной системы растет сильнее за счет появления отдельно стоящих рекуррентных точек и вертикальных и горизонтальных линий; края относительно гладкой структуры (рис. 3) становятся «мохнатыми» (рис. 4), размываются. Визуально заметно появление отдельных точек и линий, отображающих шумовую составляющую.
Рассмотрим графики изменения RR(s) для системы (7) при значениях п = 0.0, 0.2, 0.7 . Вычисление RR(s) проводилось, начиная со значения S0 = 0.01, с шагом As = 0.05 . На рис. 5 представлены графики RR(s) и RR'(s) для системы (7) при значениях п = 0.0, 0.2, 0.7 . Видно, что график RR(s) разбит на две части - практически линейный рост в диапазоне s^ 0.2k 1, затем - переход к более крутому росту после s = 1.1. Такая форма графика объяснима формой фазовой траектории системы (рис. 2а) - в районе значения s « 1 в окрестность i -й точки траектории начинают попадать точки из другой половины окружности, что и обуславливает более интенсивный рост значений RR.
а б в
Рис. 3. Рекуррентные диаграммы модельной системы при п = 0.0:
а - s = 0.3 , б -s = 0.7 , в - s = 1.1
а б в
Рис. 4. Рекуррентные диаграммы модельной системы при п = 0.7:
а - s = 0.3 , б - s = 0.7 , в - s = 1.1
Ч»М <|>М
а б в
Рис. 5. Графики изменения RR (сверху) и RR' (снизу) для модельной системы
при (s0 = 0.01, As = 0.05): а - п = 0.0, б - п = 0.2, в - п = 0.7
Резкие выпады на графиках RR'(s) (рис. 5б, в) обусловлены несовершенством модели и дискретностью временного ряда и ряда 8^ . Для нивелирования выпадов целесообразно проводить аппроксимацию значений RR'. На приведенном графике аппроксимированная кривая изображена толстой пунктирной линией; для ее получения использовалась функция аппроксимации сплайнами, входящая в стандартную библиотеку системы МайаЬ.
Следует отметить влияние шума на вид графиков RR(8), выражающееся в более крутом росте значений RR(8), особенно во второй половине графиков. Для наглядно-
сти на всех трех графиках проведена линия пересечения кривой с прямой, соответствующей sk = 16 . Видно, что при п = 0.7 эта линия расположена выше.
Значения размаха (6) ARR' = STAB составили:
STABn=0.0 = 0.0400619;
STABn=0.2 = 0.0413922;
STAB^0.7 = 0.0432716.
Очевиден рост размаха с ростом нестабильности траектории. Таким образом, меньшие значения STAB будут соответствовать большей стабильности траектории.
Проверим зависимость STAB от длины временного ряда. Для этого используем временные ряды системы Лоренца (рис. 6).
а б в
Рис. 6. Фазовые портреты системы Лоренца (стандартные параметры r = 28 , а = 10, b = 8/3) для временных рядов длиной: а - 2000, б - 4000, в - 8000 точек
Известно, что увеличение длины временного ряда улучшает достоверность результатов. Однако может оказаться, что мера STAB обладает зависимостью от длины временного ряда. На рис. 7 представлены графики RR(s) и RR'(s) для рассматриваемых рядов системы Лоренца.
а б в
Рис. 7. Графики изменения RR (сверху) и RR' (снизу) временных рядов системы Лоренца длиной: а - 2000, б - 4000, в - 8000 точек
Значения размаха (6) ARR' = STAB составили:
STABn=2000 = 0.0686575 ;
STABn=4000 = 0.0682356;
STABn=8000 = 0.0682186.
Видно, что значения STAB практически равны; наблюдается лишь незначительное уменьшение по мере увеличения длины ряда. Это объясняется, в частности, хаотическим характером системы Лоренца. Также можно говорить о некотором улучшении достоверности получаемого значения. В целом можно сделать вывод о независимости
STAB от длины временного ряда, что подтверждается также анализом других временных рядов.
Различия между STAB и DET
Мера детерминизма (DET) [7] учитывает диагональные линии рекуррентной диаграммы с использованием частотного распределения длин диагональных линий
Ps/()={{; i = 1k Ni }, где Ni - абсолютное количество диагональных линий (каждая
линия считается только один раз). Процессы со стохастическим поведением могут порождать очень короткие диагонали либо вообще не порождать их, в то время как детерминистские процессы дают длинные диагонали и малое количество отдельных рекуррентных точек. Таким образом, отношение рекуррентных точек, составляющих диа-
IP s(i)
гональные структуры, к их общему количеству DET =-№-, называется мерой
i, jR i,j
детерминизма (determinism, DET) или предсказуемости системы. Следует отметить, что эта мера не имеет значения реального детерминизма процесса. Пороговое значение минимальной длины lmin исключает диагональные линии, образованные тангенциальным движением траектории в фазовом пространстве. Очевидно, что lmjn = 1 ^ DET = 1.
Отличие детерминизма от предлагаемой меры STAB заключается в том, что первая оценивает периодичность траектории в смысле близкого прохождения участков траектории в фазовом пространстве и зависит, прежде всего, от характера изучаемой системы. Так как пороговое значение lmin позволяет отсечь шумовую составляющую диаграммы - т. е. отдельно стоящие точки и короткие линии - то при достаточно малом уровне шума значение DET будет постоянным при одном и том же s; мало того, в случае существенного увеличения зашумленности значение DET будет повышаться при сохранении s за счет высокой плотности (вплоть до объединения в достаточно длинные диагональные линии) порождаемых шумовой составляющей рекуррентных точек. В то же время STAB является мерой стабильности траектории и позволяет оценить влияние шумов и случайных флуктуаций вне зависимости от формы траектории (т.е. характера изучаемой системы).
Применение STAB
Особенность предложенной меры заключается в том, что она может быть использована только для сравнения рядов, полученных при разных измерениях одних и тех же переменных состояния одной и той же системы или одинаковых систем. Обусловливается это тем, что другой набор переменных состояния приведет к совершенно иному фазовому портрету, как следствие - совершенно иной рекуррентной диаграмме и поведению
RR(s) . Данная мера может быть использована:
• для оценки стабильности системы с целью, например, подстройки ее параметров (полученное значение сравнивается с эталонным);
• для оценки влияния шумов и искажений в результате некоторых преобразований;
• для оценки стабильности связанных систем, одна из которых некоторым образом определяет поведение другой.
Рассмотрим пример. Солнечная активность, количественно выражаемая числами Вольфа, оказывает существенное влияние на геомагнитное поле Земли. Рекуррентная диаграмма временного ряда чисел Вольфа за период с 01.1868 г. по 12.2000 г. представлена на рис. 8а.
а б
Рис. 8. Рекуррентные диаграммы чисел Вольфа (а) и аа-индекса (б). Временные ряды взяты за период с 01.1868 по 12.2000 с шагом в 1 месяц
Для характеристики общепланетарной возмущенности геомагнитного поля Земли сконструированы различные индексы, характеризующие возмущенность поля на различных широтах. Одним из таких индексов является аа-индекс, получаемый усреднением локальных индексов, полученных в геомагнитных обсерваториях на средних широтах. Рекуррентная диаграмма ряда аа-индекса за период с 01.1868 г. по 12.2000 г. представлена на рис. 8б.
пигпЬвг* ¿л-|п4«1
а б
Рис. 9. Графики изменения RR (сверху) и RR' (снизу) для чисел Вольфа (а) и аа-индекса (б) для диаграмм временных рядов с 01.1868 по 12.2000
Известно, что солнечная активность оказывает влияние на геомагнитное поле Земли. Этот факт отчетливо виден из построенных для одного временного промежутка диаграмм на рис. 8. 11-летний цикл солнечной активности явно просматривается в динамике аа-индекса. Также на обеих диаграммах видна область относительно низкой солнечной активности, завершающаяся в районе 1937 года (граница фазового перехода). Так как оба временных ряда, очевидно, отображают связанные процессы, один из которых (солнечная активность) модулирует другой (геомагнитная активность), вычисление STAB для этих рядов позволит количественно оценить, насколько сильно влияние солнечной активности на геомагнитное поле на данном промежутке времени
(рис. 9). Из изложенного выше понятно, что уровень влияния будет определяться близостью STABaa к STABW.
Рассмотрим полученные значения STAB для рядов чисел Вольфа (STABw ) и аа-индекса (STABaa) на полном диапазоне дат и на промежутках до 12.1936 года и после 01.1937 года.
Даты 01.1868 - 12.2000 01.1868 - 12.1936 01.1937 - 12.2000
STABW 0.0644426 0.0534079 0.0308646
STABaa 0.0947625 0.0573547 0.0527186
HW 0.93 0.873 0.909
H aa 0.92 0.908 0.748
Таблица. Значения STAB и показателя Харста, полученные для чисел Вольфа
и аа-индекса
Близкие значения STAB на промежутке 01.1868-12.1936 свидетельствуют о большей связи между солнечной и геомагнитной активностями (что также видно и по рекуррентным диаграммам - см. рис. 8). В то же время с 1937 года, несмотря на возросшую солнечную активность и возросшую стабильность ряда чисел Вольфа, влияние последних на магнитосферу Земли меньше.
Полученные результаты хорошо согласуются с приведенными в [10] расчетами показателя Харста (см. в таблице значения Hw и Haa) для тех же интервалов. По-видимому, это связано со структурой потока солнечной плазмы. При умеренной и слабой солнечной активности числа Вольфа определяют структуру потока солнечной плазмы, а взаимодействие с магнитосферой Земли почти линейно. В периоды высокой солнечной активности числа Вольфа уже не определяют структуру и параметры потоков солнечной плазмы. Так, по данным авторов работы [11] и цитируемой в ней литературы, скорость солнечного ветра не меняется в фазе с числами Вольфа, а в отдельных циклах их изменения существенно различны. Это подтверждается опубликованными на сайте «Энциклопедия Кругосвет» [12] данными измерения при помощи космических аппаратов IMP-8 и Voyager-2 среднего (за 300 дней) динамического давления солнечного ветра в районе орбиты Земли (на 1 АЕ) в течение одного 11-летнего солнечного цикла с 1978 по 1991 гг.
Заключение
В статье предложена новая мера количественного анализа рекуррентных диаграмм - мера стабильности траектории STAB . Достоинством данной меры является отсутствие необходимости подбора радиуса окрестности s . Рассмотрены выражения для вычисления этой меры. Приведены примеры получения STAB при исследовании модельных систем, показано влияние нестабильности фазовой траектории на значение меры. Дана оценка влияния длины временного ряда на значение STAB . Рассмотрены различия между стандартной мерой детерминизма DET и STAB .
Рассмотрены некоторые из возможных областей применения данной меры, приведены требования к временным рядам.
Литература
1. Manuca R., Savit R. Stationarity and nonstationarity in time series analysis. // Physica D 99 (2-3). 1996. Р. 134-161.
2. Eckmann J.-P., Kamphorst S.O., Ruelle D. Recurrence Plots of Dynamical Systems. // Europhysics Letters 5. 1987. Р. 973-977.
3. Киселев В.Б. Некоторые методы нелинейного анализа // Научно-технический вестник СПбГУ ИТМО. 2005. В. 20. С. 172-180.
4. Киселев В. Б. Рекуррентный анализ - теория и практика // Научно-технический вестник СПбГУ ИТМО. 2006. В. 29. С. 118-127.
5. Takens F. Detecting Strange Attractors in Turbulence. Berlin% Springer, 1981. Р. 366381.
6. Киселев Б.В., Волобуев Д.М. Реконструкция аттракторов в трехмерном фазовом пространстве и моделирование геомагнитных пульсаций. // Вопросы геофизики. СПб: Изд. С.-Петербургского унив., 1998. В. 35. С. 338-348.
7. 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.
8. Iwanski J.S., Bradley E. Recurrence plots of experimental data: To embed or not to embed? // Chaos. 8 (4). 1998. Р. 861-871.
9. Webber Jr.C.L. Recurrence Quantification Analysis, 2003. URL: http://homepages.luc.edu/~cwebber
10. Киселев Б.В., Киселев В.Б. Различия в динамике солнечной и геомагнитной активности. // Вопросы геофизики. Вып. 39 СПб: Изд-во С.-Петербургского ун-та, 2006.
11. Вальчук Т.Е., Фельдштейн Я.И. Корреляционные и регрессионные соотношения между индексом aa геомагнитной активности и характеристиками околоземного космического пространства. / Солнечный ветер, магнитосфера и геомагнитное поле. М., Наука, 1983. 127 с.
12. http://www.krugosvet.ru/articles/110/1011048/1011048a2.htm
