Исследование по методу Прони динамики систем на основе временных рядов Текст научной статьи по специальности «Математика»

УДК 541.182

Б.Г. Кухаренко

Институт машиноведения им. А.А. Благонравова РАН Московский физико-технический институт (государственный университет)

Исследование по методу Прони динамики систем на основе временных рядов

Рассматриваются локально-переходные временные ряды для переменных динамических систем. Описан метод определения спектральных характеристик локально-переходных временных рядов, обеспечивающий поиск подобных сегментов. Показано, что нелинейные временные ряды, описывающие хаотическую динамику в моделях стохастического резонанса, Росслера и Лоренца, являются локально-переходными. Поиск подобных сегментов нелинейных временных рядов позволяет исследовать локальную во времени динамику хаотических систем.

Ключевые слова: динамика систем, временные ряды, спектральный анализ, разложение сингулярных чисел, поиск подобных сегментов, стохастический резонанс, модель Росслера, модель Лоренца.

I. Временные ряды для переменных динамических систем

1.1. Модель локально-переходных временных рядов

Динамика большинства систем представляется нестационарными временными рядами, и это создаёт существенную проблему при их анализе. При моделировании — это значения функций х = х(Ь) в дискретные моменты времени гДЬ (ДЬ — шаг дискретизации времени) и численные решения моделей динамических систем, а при мониторинге — записи колебаний различных переменных, характеризующих состояние этих систем. Поскольку эти временные ряды (которые могут быть многомерными) нестационарные, их спектральные характеристики меняются со временем. Стандартные процедуры локального (во времени) анализа временного ряда х[1 : N0] используют его сегментирование посредством сдвига временного окна фиксированной длины N ^ N0 NДЬ определяет временной масштаб для локального анализа этого нестационарного временного ряда). При этом во многих исследованиях неявным является предположение о локальной стационарности временных рядов [1-5]. Предполагается, что зна-

чение N может быть выбрано так, чтобы обеспечить квазистационарность последовательных сегментов х[п : п + N — 1],] п : п + N — 1] С [1 : N0]. В настоящей работе используется модель локально-переходного изменения, которая не столь широко распространена для определения локальной во времени динамики систем на основе нестационарных временных рядов [6-9]. Модель локально-переходного изменения характеризует временные интервалы, в которых нарушается локальная стационарность временного ряда. Это может быть результатом присутствия быстро затухающих мод, внезапных изменений структуры динамической системы и быстрых локальных изменений модальных параметров. Цель настоящего исследования — показать, что среди нестационарных временных рядов, представляющих временные зависимости переменных динамических систем, локально-переходные временные ряды являются более представительным классом временных рядов, чем локально-стационарные временные ряды, и установить, какие режимы колебаний (паттерны) характерны для локально-переходных временных рядов. Таким образом, для нестационарных временных рядов, исследуемых в настоящей работе, квазистационарность сегментов этих (нестационарных) временных рядов является скорее исключением. Локально-переходные временные ряды для переменных динамиче-

ских систем характеризуются колебательным изменением. Это значит, что после удаления среднего значения х = ^^ х [г], то есть

[г],г = 1 : —^ — х)4 = 1 : Щ,

либо удаления линейного или нелинейного тренда f (¿) в результате регрессии

\

N0

(х[г] — f (гАЬ))2 — шш

г=1

и преобразования х И , г = 1: N0 - (х [г] — f (г • АЬ)), г = 1 : N0,

для преобразованного временного ряда х[1 : N0] может быть определена последовательность локальных максимумов

х[гк], к =1: М, М < N,

где

[гк,к = 1 : М] С [1 : N0],

таких, что для каждого локального максимума существует е-окрестность:

х[г] < х[гк],если \г — гк\ < е.

Ряды, элементы которых после удаления тренда f (Ь) меняют знак при изменении индекса на единицу, считаются приводимыми к локально-стационарным временным рядам и не рассматриваются. Определенная выше последовательность локальных максимумов х[гк], к = 1 : М — это (дискретная с неравномерным шагом) огибающая временного ряда х[г], г = 1 : N0. Она позволяет определить класс локально переходных временных рядов с колебательным изменением, рассматриваемых в настоящей работе. Для этих временных рядов х[1 : N0] имеется временной масштаб NАЬ, такой, что в последовательных сегментах

х [из,из + N — 1], из = N • .1,

3 = [0 : round(No/N)]

их огибающая х[гк ], к = 1 : М, меняется монотонно. Такие сегменты (а значит,

и сам ряд локально во времени) характеризуются переходным колебательным изменением. Имеются примеры, когда это переходное колебательное изменение характеризует временной ряд в целом (то есть глобально во времени). Рассматриваемый класс временных рядов с локальным во времени переходным колебательным изменением может и дальше конкретизироваться (с привлечением априорного знания), чтобы обеспечить решение задачи идентификации различных режимов колебаний (паттернов) конкретной динамической системы на основе её временного ряда. В исследованиях нестационарных временных рядов рассматриваются задачи, решение которых может способствовать разработке методов идентификации режимов колебаний (паттернов) локально-переходных временных рядов. Первая задача — это поиск подобия сегментов временных рядов. Используя временные окна различной длины, делается попытка найти во временном ряде сегменты, повторяющиеся с точностью до некоторой меры близости в пространстве размерности, равной длине используемого временного окна. Следует отметить работу [10] по разработке методов для сокращения размерности пространства поиска. В этих методах при проверке сходства сегментов используются преобразования этих сегментов, например, дискретное преобразование Фурье, а для сокращения размерности — теорема Парсеваля. Однако используемые преобразования сегментов нестационарных временных рядов носят формальный характер, и не предполагается информативность (физическая интерпретация) их результатов. Вторая задача — это обнаружение (локализация во времени) внезапных изменений зависимости элементов ряда от индекса (временной зависимости, паттерна) временного ряда, в том числе внезапных изменений его спектра частот [1-3].

1.2. Метод Прони

Требования к методу спектрального анализа основываются на используемой модели локально-переходного изменения временных рядов. То есть сегменты ряда — это не квазистационарные части нестационарного временного ряда (которые для наблюдаемого ряда могли бы

оказаться очень малыми по длине), а переходные, то есть с неизменным знаком демпфирования главных спектральных составляющих сегмента, что отражает специфику рассматриваемых в настоящей работе (гладких) локально-переходных временных рядов. Использование спектрального анализа, соответствующего модели временного ряда, делает его результаты информативными (то есть имеющими физическую интерпретацию). При спектральном анализе размеры сегментов ряда могут существенно варьироваться, в результате чего меняется временной масштаб локального анализа временного ряда. Поскольку предстоит анализировать локально-переходные временные ряды, амплитуда колебаний которых убывает или возрастает во временном окне ограниченной продолжительности, интерес представляет метод Прони (Ргопу И,., 1796) [11]. Исторически — это первый метод спектрального анализа временных рядов, но он не имеет непрерывного аналога. Его широкое использование откладывалось из-за отсутствия устойчивого с вычислительной точки зрения алгоритма этого метода. Декомпозиция Прони сегмента

временного ряда х[1 : N] имеет вид р

х[г]=^ гЩгЩ- + п[г], г =1: N, 1=1

где р — число определяемых полюсов сегмента временного ряда; г[1] = ехр (8Щ + з2п/[/]) ,1 = [1 : р] — полюса, которые определяются для сегмента временного ряда, где $[/] и / [/] — соответственно демпфирование и частота; г[1] = а[1] • ехр(з'^[/]) ,1 = [1 : р] — вычеты в этих полюсах, где а[1] и (р[1] — соответственно амплитуда и фаза; п[/] — аддитивный шум. Создание устойчивого с вычислительной точки зрения алгоритма метода Прони является результатом работ ряда авторов [12-17]. В этих работах решается задача спектрального сверхразрешения частот. Матричный алгоритм метода Прони использует ганкелеву матрицу данных У размерности N — V) х (Ь + 1):

У [: ,1] = х [I : (I + N — Ь)]Т , I = [1 : (Ь + 1)] ,

где L — размерность вложения. Показано, что z = z[l],l = [1 : p],p < L < N — p понижает ранг матричной формы (matrix pencil)

(у(i) — zY,

где

Y(0) = Y [: ,1 : L] ,

Y(i) = Y [: ,2: (L + 1)] .

Критерий, позволяющий установить число p главных полюсов сегмента временного ряда, использует сингулярное разложение

Y(v) = U(v) • S(v) • (V(v))Tv =[0:1] , где

S(v) = diag(s(v)[1],s(v)[2], ...), s(v) [1] > s(v) [2] > ... > 0

и

(U(V) )T U(V) = I(N-L)x(N-L),

(V(V) )TV(V) = I(L+i)x(L+i), а верхний индекс «T» обозначает транспонирование. Определяются усеченные матрицы

S(V) = diag (s(v) [1], ...s(v)[p], ..., 0) , v =[0:1] .

Использование матриц S^, v = [0:1]

позволяет определить матрицы Yr , v = [0:1] с фиксированным рангом p. Определение полюсов zi, l = [1 : p] как обобщённых собственных чисел матриц ,v = [0:1] гарантирует вычислительную стабильность определения этих полюсов в пределе временного ряда без шума, поскольку матрицы Y^v =[0:1] имеют фиксированный ранг p. Аппроксимация с минимальной размерностью p =2 полезна при спектральном оценивании среднего периода гладкого временного ряда. Работы по совершенствованию алгоритмов на основе метода Прони для временных рядов, представляющих переменные динамических систем с целью их идентификации, начались в 1990-х годах и продолжаются до сих пор [18-20]. В современных алгоритмах метода Прони отсутствуют какие-либо ограничения на выбор параметров N и p. Однако посредством выбора параметров N и p можно обеспечить информативность результатов спектрального оценивания конкретного локально-переходного

временного ряда. Основной вопрос, который требуется решить при спектральном анализе временного ряда по методу Про-ни, — это выбор размера N временного окна и выбор величины сдвига этого временного окна (с перекрытием последовательных сегментов или без него). Правильный выбор размера временного окна позволяет существенно сократить объём информации по сравнению с исходными сегментами временного ряда и в результате оптимально решить задачу поиска подобия в локально-переходном временном ряду (то есть построить абстракцию режимов колебаний (паттернов) этого временного ряда). Специальный интерес представляет оценка момента возникновения неустойчивости с целью прогнозирования (на некоторый временной интервал) дальнейшей эволюции динамической системы, представляемой анализируемым локально-переходным временным рядом. С оформлением разведывания данных (Data Mining) как самостоятельной дисциплины стало общепринятым демонстрировать, что предлагаемые алгоритмы способны анализировать произвольные временные ряды. В настоящей работе метод Прони применяется с учётом особенностей (конкретной модели) локально-переходных временных рядов (то есть используется априорное знание). Именно поэтому получаемые результаты являются информативными (то есть имеют физическую интерпретацию). Например, при обезличенном подходе отсутствует возможность выбора естественного (для анализируемого временного ряда) масштаба времени.

I.3. Нелинейные временные ряды

Ниже показано, что нелинейные временные ряды, представляющие хаотические колебания, являются локально-переходными временными рядами. Проведено детальное исследование режимов колебаний (паттернов) в численных решениях классических диссипативных моделей (модель стохастического резонанса, модель Росслера, модель Лоренца) в области значений параметров, соответствующей хаотической динамике, и установлены закономерности локальной во времени эволюции этих хаотических численных решений, в том числе бифуркации изменения периода

колебаний. Что касается особенностей истинных траекторий нелинейных динамических систем, то установлено, что нелинейная система с хаотическим поведением характеризуется двумя основными свойствами [21-23]. Во-первых, в пространстве состояний нелинейной системы имеется бесконечное число нестабильных периодических орбит, вложенных в странный аттрактор. Во-вторых, динамика в этом странном аттракторе обладает свойством эргодичности. Это означает, что в процессе временной эволюции фазовая траектория системы посещает малую окрестность любой точки у каждой неустойчивой периодической орбиты, вложенной в аттрактор. Главная причина, почему периодические орбиты трудно обнаружить, — это их нестабильность: фазовые траектории в окрестности нестабильной периодической орбиты отталкиваются от неё [23-25]. При произвольно выбранных начальных условиях фазовой траектории в среднем требуется большое число итераций, чтобы попасть в е-окрестность целевой периодической орбиты. Более того, распределение вероятностей за большой отрезок времени, генерируемое типичной фазовой траекторией на странном аттракторе, в общем случае имеет сингулярный характер. Фазовая траектория с произвольно выбранными начальными условиями в области притяжения посещает различные части аттрактора с заметно отличающимися вероятностями. В этом смысле странный аттрактор обладает мультифрактальной структурой [26]. Действительно, если нелинейная система описывается моделью с конечным числом степеней свободы, то в её пространстве состояний всегда имеется несколько положений равновесия. Тогда в этом пространстве переменных нелинейной системы любое хаотическое движение представляется последовательностью локальных движений возле нескольких её положений равновесия. Причём каждое из локальных движений возле одного из положений равновесия может завершиться в результате перескока к другому положению равновесия. В результате в пространстве состояний нелинейной модели любое хаотическое движение представляется в виде последовательности неустойчивых локальных периодических орбит.

Все временные ряды, представляющие хаотические движения, считаются гладкими. Если временной ряд не гладкий, его можно сделать гладким с использованием разложения сингулярных чисел ганке-левой матрицы данных. Как и в матричном алгоритме метода Прони, для у[1 : Щ, N ^ N0, представляющего собой некоторый сегмент полного временного ряда у[1 : N0], определяется ганкелева матрица данных У размерности N — Ь) х (Ь + 1):

У [: ,1] = у [I :(1 + N — Ь)]Т ,

I = [1 : (Ь + 1)] ,

где Ь — размерность вложения. Этот анализ основан на сингулярном разложении матрицы У, а именно:

У = иБУ

т

1 ^ к ^ N.

В последнее время описанная процедура сглаживания временных рядов получила название сингулярного спектрального анализа [27-29]. Из этого описания следует, что в описываемом ниже анализе хаотических временных рядов не требуется их предварительное сглаживание посредством сингулярного спектрального анализа. Матричный алгоритм метода Прони сам основан на использовании матриц данных У(0) = У [: ,1 : Ь] и У(1) = У [: ,2 : (Ь + 1)], ранг которых фиксируется посредством разложения сингулярных чисел [15-20].

II. Поиск подобия в нелинейных временных рядах

где

и

в = ^(в [1] ,5 [2] , ...) 5 [1] >5 [2] > ... > 0

иТ и =

N-ь)хN-ь)1

у у = 1(Ь+1)х(Ь+1),

а верхний индекс «Т» обозначает транспонирование. Определяется усеченная матрица

Б^ = dгag (5 [1], ..., в [р] ,0, ..., 0).

Использование в сингулярном разложении матрицы Б^ позволяет определить матрицу У& ранга р. Усреднение антидиагоналей матрицы У^ даёт Ганкелеву матрицу, по которой может быть восстановлен сглаженный временной ряд (считается, что Ь < N — Ь — 1):

к

У^ Угг [т,к — т + 1],

т=1

1 ^ к ^ Ь + 1,

Ь+1

1

У^ Угг [т,к — т + 1],

Ь +1 ,

т=1

Ь + 1 ^ к ^ N — Ь,

Ь+1

У^ Угг [т,к — т + 1]

m=k-N+Ь+1

N — к + 1 N — Ь < к < N

П.1. Использование теоремы Парсеваля для поиска подобных сегментов

Считается, что у всех частотных компонент сегмента {х[к],к = 1 : N} их временные зависимости симметричны относительно начала отсчёта времени (то есть одинаково затухают при Ь — +с и при Ь — —сю). Тогда спектр мощности

м

м

1=1

аЩбЩещАМ])

+ с.с.

Б(со) = -п

,=| и и , и и

(1)

где а[1] и (р[1] — амплитуды и фазы вычетов, $[/] и и[/] = 2пf [/],/ = 1 : М — соответственно фактор демпфирования и круговые частоты полюсов, а с.с. — комплексно-сопряженный компонент. С учётом (1) для разности двух сегментов со спектрами, подобными по норме Ьж, оценка для полного спектра мощности при условии, что фазы различных частотных составляющих считаются некоррелированными:

N

Е(х[к] — У[к])2

к=1

2 п

dиx

х

м

м

1=1

И ехрОУа

5Х[1\2 + Ш] - ^У2

2

1

к

[l] exP(jVy

Sy [l]2 + (Uy [l] —

+ c.c.

M

M

ax[l] ay [l]

1

x

8[l]3duj [6[1]2 + (u[l]2 - u)2]2'

x

(2)

В (2)

п

S[l]3du

2 + (u[l] — и)2)2 2'

поэтому аналог теоремы Парсеваля для разложения Прони разности двух подобных сегментов имеет вид

N

£№] — уШ

k=i

M

ax[l] ay [l]

-V,

2tt\VW] л/Ш,

(3)

Поскольку в методе Прони изначально (a priori) отсутствует некоторая предопределённая решётка частот, поиск сегмента, подобного исходному, происходит в два этапа:

1) поиск сегментов, спектры частот которых близки спектру частот исходного сегмента (при этом используется мера близости частот L— она не зависит от нормировки временного ряда);

2) поиск сходства по норме (3) в пространстве с координатами ja[l]/ (S[l])^2 , l = 1 : M}, представляющими собой набор фурье-амплитуд всех частотных составляющих для последовательных сегментов временного ряда (используется мера близости L2 — она зависит от нормировки временного ряда).

Таким образом, если мы оцениваем сходство по M комплексным составляющим сегментов, размерность соответствующего прони-пространства поиска оказывается 2M. Но тот факт, что метод Прони позволяет определять фактические, а не формальные спектральные составляющие сегмента, позволяет использовать редуцированное (сокращенное) прони-простран-ство поиска размерности M' < M. Это

приводит к тому, что при поиске проверяется /х[1] & /у[1],1 = 1 : М' только для главных частотных составляющих, а из (3)

N

5>и — УИГ

>

>

k=i

.. M'

-У

2

i=i

ax [l] ay [l]

VЩГ] VW]

что гарантирует отсутствие ложных (ошибочных) пропусков подобных сегментов. Как обычно, ложные идентификации удаляются при последующем сравнении во временной области. Ниже описан алгоритм идентификации почти периодических орбит на основе определения минимального временного масштаба ряда по результату состоявшейся одночастотной аппроксимации этого временного ряда по методу Прони (с использованием ганкелевой матрицы данных с размерностью вложения Ь = N/2 со снижением ранга до 2).

П.2. Алгоритм сегментирования нелинейных временных рядов по методу Прони

У нелинейного временного ряда, представляющего хаотическое движение, имеется характерный средний период колебаний То (средний период неустойчивых периодических орбит в пространстве состояний (фазовом пространстве) его динамической системы). Спектр частот гладкого квазипериодического нелинейного временного ряда, определённый за его текущий период, должен быть квазигармоническим. Однако фактически периодичность либо квазипериодичность временного ряда проявляется за два его последовательных периода [30]. Поэтому под спектром за период ниже будет подразумеваться дискретный спектр частот, определённый для двухпериодного сегмента, состоящего из текущего однопериодного сегмента и предшествующего ему. Именно такой средний спектр, определённый для сегмента из двух последовательных однопериодных сегментов, является квазигармоническим с фундаментальной частотой, усреднённой за два периода локальной орбиты в пространстве состояний (фазовом пространстве). Поэтому минимальный (атомный) сегмент, используемый для отождествления (установления е-подобия) по методу

2

2

2

1

1

2

Прони, имеет длину 2T0/At. Перекрывающиеся (наполовину) двухпериодные сегменты временного ряда хк[1 : N0],K = 1 : P определяются следующим образом:

Xjk[1 : N] = хк[nj : nj+2], nj = NJ/2,

J = 0 : R, (4)

где N < N0,R = round (2N0/N). Считается, что в интервале времени

Atj = At[nj : nj+2],

где nj определено в (4), изменение частоты нелинейного временного ряда происходит непрерывно. Спектр частот сегментов Xjk [1 : N] (нестабильных орбит периода 1) с непрерывной зависимостью спектра частот от номера J сегмента, оцениваемый в среднем за два последовательных периода AtJ нелинейного временного ряда хк [1 : N0], имеет вид

fjK[l] « l • fjK[1],

l = 1 : Mj, J = 0 : R, K = 1 : P,

где fjK [1] ~ foj — фундаментальная частота нестабильной орбиты периода 1 (усреднённая за два последовательных периода At j ). В случае временного ряда, представляющего одночастотное колебание, в этом квазигармоническом спектре только фундаментальная частота имеет физический смысл (в пространстве состояний или фазовом пространстве). Для произвольного осциллирующего переходного временного ряда хк [1 : N0] с гладким изменением

fjK[l] <>l • fjK[1],

l = 1 : MJ, J = 0 : R, K = 1 : P.

Поскольку в методе Прони точность определения спектра частот fjK [l] ~ l • fjK [1] ,l = 1 : Mj не зависит от частоты дискретизации (At)-1, мерой квазипериодичности двухпериодного сегмента Xjk [1 : N] может служить сумма квадратов отклонений нормированных частот от натурального ряда чисел

Mj

E((fjK[l]/fjK[1]) - l)2.

1=1

Таким образом, все частоты спектра двух-периодного сегмента больше или равны фундаментальной частоте ^ & 1/Т0. Для рассматриваемых нелинейных временных рядов это свойство спектра частот является общим. Представляют интерес специфические спектральные характеристики осциллирующих нелинейных временных рядов. Это достигается использованием в алгоритме сегментирования по методу Про-ни стратегии, состоящей в слиянии минимальных сегментов длиной 2Т0/АЬ в более крупные сегменты. Так последовательно увеличивается временной масштаб спектрального анализа.

Для каждого значения 3 = 0 : (Я — 3) определяется последовательность расширенных интервалов времени являющихся результатом слияния последовательных двухпериодных интервалов:

Ь

{ТЬ = иАЬ/,Ь =(3 + 2): (Я — 1)}. (5) I=з

Каждому временному интервалу Т Ь (5) соответствует расширенный сегмент, для которого по методу Прони определяются его спектры амплитуд, фаз, факторов демпфирования и частот:

{аЬЩ,0Ь[1],%ШЬ[1],1 = 1: МЬ}. (6)

То есть спектры (6) могут быть определены для частично перекрывающихся расширенных интервалов времени (5). Спектры (6) расширенных сегментов временного ряда х[1 : N0] в принципе могут содержать более низкие частоты, чем средняя частота ^ = Т-1 локальной орбиты периода 1 в пространстве состояний (фазовом пространстве). Это связано с тем, что минимальная частота спектра, определяемая по методу Прони для сегмента, соответствующего временному интервалу продолжительностью Т, оказывается порядка Т-1. Пусть для 3 = Я1 : Я2, где 0 < Я1 < Я2 < Я, ни сами временные интервалы Тз = АЬ[из : из+2], где из = N3/2, ни их границы Ьз = АЬ • из не содержат моменты времени = АЬ • т1,1 = 1 : ф, в которые происходит резкое изменение частоты нелинейного временного ряда (численные оценки

для резкого изменения частоты нелинейного временного ряда приведены ниже). Тогда спектр частот последовательных двух-периодных сегментов (4) меняется непрерывно, оставаясь при этом почти гармоническим. Поэтому при усреднении спектров частот последовательных двухпери-одных сегментов по временному интервалу : N^2] по формуле

и(х) под воздействием однородного высокочастотного поля имеет вид

1

R2

R2 — Ri

Efj [1

mjq)]

J =Ri

будет получен почти гармонический спектр частот. Оказывается, что (средний) спектр частот сегмента x[NR1 : NR2 ], соответствующего интервалу времени At[NR1 : NR2], также является почти гармоническим. Но поскольку происходит кратное увеличение длины анализируемого сегмента (по сравнению с двух-периодным сегментом), в спектре могут появиться дополнительные энергетически значимые частоты < T—1. Таким образом, энергетически значимая часть спектра частот сегментов, кратных двухпериодно-му сегменту, формируется исключительно в относительно низкочастотной области < T-1. Это и гарантирует конечную прони-размерность сегментов нелинейных временных рядов, представляющих хаотическую динамику.

III. Модели динамических систем

Ш.1. Модель стохастического резонанса

В работе [31] показано, что роль случайного поля может играть однородное поле с быстрой гармонической зависимостью от времени, то есть сильно нелинейная динамика типа стохастического резонанса имеет место для бистабильного потенциала и(х) = —^х2 + |х4 с возмущением вида

и (х,г) = и (х) + х(7 сов(2п/0£) +

+е сов(2п/1*)), е < 1, /1 < /о,

где 7 и /о — соответственно амплитуда и круговая частота высокочастотного гармонического воздействия. Уравнение движения частицы в бистабильном потенциале

d x _ dx 3

w+6dt-x+x

Y cos (2nfot), (7)

где 8 — линейное демпфирование. Однородное высокочастотное поле порождает хаотические движения частицы в биста-бильном потенциале и (х). В фазовом пространстве системы с хаотической динамикой «стохастическое море», принадлежащее одной хаотической траектории, содержит узкие области — островки, куда эта траектория не может попасть из области «стохастического моря» и наоборот [32]. Если такая узкая область характеризуется сильным притяжением по супер- либо субгармонике, то в ней могут существовать регулярные движения [21]. В действительности, каждый островок «стохастического моря» имеет такую же сложную структуру, как и само «море». Таким образом, фазовое пространство системы с хаотической динамикой имеет фрактальную структуру (сосуществование областей устойчивости и областей хаоса в фазовом пространстве создаёт значительные трудности в изучении динамики хаоса). В пространстве параметров динамической системы (7) также обнаруживаются (локализуются) области регулярной и хаотической динамики. Для этого при различных значениях параметров системы определяются её численные решения, которые идентифицируются как регулярные либо хаотические. То есть, как и в работе [33], дискретные траектории х[г],г = 1 : N при различных значениях амплитуды 7 и частоты /0 внешнего возбуждения и линейного демпфирования 8 представляются численными решениями уравнения (7) с постоянным шагом интегрирования. Существуют диапазоны значений параметров 7, /0 и 8, в которых аддитивное воздействие 7 сов(2п/0) генерирует хаотический аттрактор, отражающий апериодические перескоки между двумя устойчивыми неподвижными точками х± = ±1 невозмущённой системы. Поскольку рассматриваемый бистабильный осциллятор сильно демпфирован, траектория системы проводит большую часть времени вблизи состояний х+ и х-. При отсутствии возбуждения третья неподвижная точка х0 = 0 ведёт себя как седло. Ее стабильное и нестабильное подмножество

играет важную роль в структуре этого хаотического аттрактора. Но в отличие от результатов работы [33] оказывается, что хаотические численные траектории могут описывать такой переходный процесс, который при определённых значениях параметров 7, f0 и 8 устанавливается и оказывается чисто периодическим [34] (рис. 1).

X

0.8 0 6 0.4 0.2 0 -0.2 -0.4 -0.6

.0.8 и-1-1-1-1-1-1-1-1-1- t

О 20 40 60 ВО 100 120 140 160 1В0 200

Рис. 1. Временная зависимость х = х(Ь) для бистабильного осциллятора

dx/dt

осциллятора

На рис. 2 представлена эта численная траектория при 7 = 0,316,^ = 0,1750 и 8 = 0,316 в фазовом пространстве системы (7). Ниже показано, как идентифицируется субгармонический режим этого численного решения. Определяется атомный двухпериодный сегмент этого численного решения (период Т0 = 1/^) и для последовательности расширяющихся сегментов вычисляется их прони-размерность, пока она не перестает меняться для временных интервалов больших 16Т0 (рис. 3) [35]. Так, установлено, что это численное решение модели бистабильного осциллятора

(7) представляет субгармоническое движение порядка 1/7. Инструментальная частота метода идентификации равна частоте возбуждения ^ = 0,1750. Системная частота субгармонических колебаний равна ^/7 = 0,0250. На рис. 4 показано отличие спектра частот f [/] ,1 = 1 : 10 от нечетных гармоник частоты f [1] = ^/7.

X

0.S -

-0.8 1-1-1-1-1-1- f

40 60 80 100 120 140 160

Рис. 3. Сегмент временного ряда для временного интервала больше 16Т0

Г

от нечетных гармоник частоты f [1] = fo/7

Если в правую часть уравнения (7) добавить медленное гармоническое воздействие е cos (2nf1t) ,е ^ 1 с частотой f1 = 0,0250, то будет иметь место резонанс. Эта резонансная частота f1 является 7-й субгармоникой частоты f0 быстрого гармонического воздействия 7 cos (2nf0t) в уравнении (7). Таким образом, частота f1 медленного гармонического воздействия не определяет резонансную частоту, как это имело бы место в случае стохастического резонанса. Она сама определяется частотой f0 быстрого гармонического воздействия, которое при стохастическом

резонансе выполняла бы роль случайного фактора — высокочастотного шума, обеспечивающего энергетическую накачку перескоков между положениями равновесия х± = ±1. Поэтому для бистабильного осциллятора (7) при выбранных значениях параметров 7, /0 и 8 имеет место обычный резонанс, а не стохастический. Ранее предполагалось, что стохастический резонанс в системе (7) всегда имеет хаотическую природу. Однако в работе [36] теоретически доказан квазипериодический характер стохастического резонанса. Описанный выше подход позволяет подтвердить эту теорию.

ПХ.2. Модель Росслера с хаотической динамикой

Модель Росслера имеет вид дх/дЬ = —у — г,

ду/дЬ = х + ау, (8)

йг/йЬ = Ь + г(х — т),

где а, Ь, с — параметры. Для модели (8) координаты второго седлового фокуса индекса 1 и первого седлового фокуса индекса 2 имеют вид

х0 = т/2 ± [т2/4 — аЬ]1/2,

уо = —хо/а, го = —уо-

Движение в модели с хаотической динамикой принято изображать в пространстве её координат (пространстве состояний). При этом получается усреднённая по времени картина колебаний, и проявляются их статистические закономерности. Однако на примере численных решений модели Росслера показано, что существуют нелинейные временные ряды, в основном представляющие стационарные субгармонические колебания. Отличительной особенностью модели Росслера является то, что при каждом значении параметра этой динамической модели возможен единственный почти периодический переходный режим колебаний и единственный субгармонический переходный режим колебаний, амплитуда которого со временем стабилизируется [37]. Ниже рассматривается модель Росслера с параметрами: а =1/5, Ь = 1/5ит — варьируемым параметром модели. При начальных условиях хо = 0,1; уо = 0; го = 0 численные решения модели (8) с начальными условиями

хо = 0,1; уо = 0; го = 0 получены с постоянным шагом интегрирования ДЬ = 0,04. При т = 4,23 завершается последовательность бифуркаций удвоения периода. Для т = 4,23 численное решение модели (8) представлено на рис. 5. Орбита порядка 12 в фазовом пространстве системы Росслера соответствует стационарному субгармоническому колебанию порядка 1/12. Спектр частот является гармоническим с фундаментальной частотой, равной субгармонике порядка 1/12 от частоты локальной орбиты.

г

15

10 -10

Рис. 5. Орбита периода 12 при т = 4,23

л

Х{1],1 = 1 : 12 от ряда целых чисел I = 1 : 12

На рис. 6 показано отличие нормированных частот А[/] = /[/]// [1] ,1 = 1 : 12 от ряда целых чисел I = 1 : 12. Спектр амплитуд имеет максимум для частоты локальной орбиты. При увеличении параметра т аттрактор Росслера возникает как

предел последовательности стационарных субгармонических движений понижающегося порядка для частоты нестабильной орбиты периода 1. Для т > 5 система Рос-слера действительно демонстрирует хаотическое поведение. Тем не менее, как показано в [38], почти гармонический спектр частот продолжает характеризовать и временную эволюцию траектории для аттрактора Росслера. Для широкого интервала значений параметра т нестабильные орбиты периода 1 идентифицируются с высокой точностью. Если аттрактор Россле-ра существует, то стабилизация системы Росслера — это выделение орбиты периода 1 из аттрактора Росслера в процессе его временной эволюции. Сходство частотных спектров субгармонических движений системы Росслера без и в присутствии аттрактора Росслера во временной эволюции приводит к тому, что орбита периода 1 выделяется из стационарных субгармонических движений системы Росслера посредством тех же методов управления, которые обеспечивают выделение орбит периода 1 из аттрактора Росслера.

Ш.Э. Структура траекторий модели Лоренца в пространстве состояний

Модель Лоренца представляет собой систему уравнений

dx/dt = —вх + ух, dy/dt = —ау + ах, dz/dt = —ху + ру — г.

(9)

В оригинальной работе Лоренца система (9) исследуется при значениях параметров: в = 8/3,а =10. При таком выборе параметров стационарное состояние теряет устойчивость при р = 470/19 & 24,74 [39]. Установлено, что во временной области все гладкие переходные колебательные процессы, соответствующие непрерывным последовательностям локальных орбит в пространстве состояний модели Лоренца, завершаются в результате резкого изменения частоты локальных орбит [40].

Обнаружение внезапных изменений в нелинейных временных рядах. Стандартные процедуры локального (во времени) анализа временного ряда х[1 : N0] используют его сегментирование посредством сдвига временного окна фиксированной длины N < N0 NАЬ опре-

деляет временной масштаб для локального анализа этого нестационарного временного ряда). При обнаружении резких изменений предполагается, что значение N может быть выбрано так, что последовательные сегменты х[и : (и + N — 1) ],] и : (и + N — 1)] С [1 : N0] временного ряда являются квазистационарными. Стандартный подход к определению моментов резких изменений основан на построении для базового сегмента х[и : (и + N — 1)] тестовой ганкелевой матрицы данных У(п) размерности N—Ь)х(Ь+1) (Ь — размерность вложения):

Т

У(и) [: ,3] = х[(и + 3) : (и + 3 + N — Ь)] 3 = [1 : (Ь + 1)]

со столбцами

Х(п) = х[(и + 3) : (и + 3 + N — Ь)]Т, 3 = [1 : (Ь +1)],

принадлежащими пространству Я1 Ь. Определение главных компонентов ковариационной матрицы (У(п))Т • У(п) связано с разложением сингулярных чисел матрицы У(п):

У(п) = и{п) Б(п)(У(п))т.

Главные компоненты — это главные левые сингулярные векторы

и(п) [1 : М], М<Ь ^ N — Ь

(эмпирические ортогональные функции). Изменения обнаруживаются по изменению эвклидова расстояния между М главными эмпирическими ортогональными функциями и столбцами:

Х(п) = х[(и + 3) : (и + 3 + N — Ь)]т,

3 = [(Р + 1) : (р + Я + 1)], Р> 1,

ганкелевых матриц данных для последующих сегментов (которые могут частично перекрываться с базовым сегментом) [41]. Определяется подпространство ЯМ С Я1-Ь с базисными векторами и(п) [1 : М]. Поскольку эти базисные векторы ортонормированные, квадрат эвклидова расстояния между произвольным вектором Z € Я1-Ь и подпространством Ям

(п)

имеет вид

^Г — ||и(п)[1 : М]т • Z

2

= ZT Z — ZT • Uw[1 : M] • U{n) [1 : M]T • Z.

Квадрат расстояния

p+1+Q

Dn;p,Q

Z x

j =p+i

<n)

T

j —

— (Xn))T • U(n) [1 : M] • U(n) [1 : M]T • X^

— это сумма квадратов расстояний всех столбцов

Xj(

(n)

x[(n + j) : (n + j + N — L)]

T

j = [(p +1):(p + Q + 1)], P> 1,

тестовой ганкелевой матрицы данных. Квадрат расстояния Dn;p,Q начинает расти, если сегмент временного ряда, для которого построена тестовая ганкелева матрица данных, содержит резкое изменение по сравнению с базовым сегментом (то есть для обнаружения изменений определяется порог). Ускорение процедуры сравнения основано на использовании других мер расстояния и подпространства Крылова [42]. С оформлением разведывания данных (Data Mining) как самостоятельной дисциплины стало общепринятым демонстрировать, что предлагаемые алгоритмы способны обрабатывать произвольные временные ряды. Однако обезличенный подход обеспечивает только констатацию изменений. Ниже описан метод, позволяющий предсказывать внезапные изменения. Это связано с тем, что метод использует модель временного ряда, позволяющую выбрать естественный (для временного ряда) масштаб времени.

Обнаружение внезапных изменений фундаментальной частоты численного решения модели Лоренца. В качестве примера рассматривается численное решение модели (9) при значениях параметров: в = 8/3, а = 10, р = 28 и начальных условиях x (0) = 0,y (0) = 0,z (0) = 10"3 (рис. 7). При этих значениях параметров существует аттрактор Лоренца [35]. Прямые линии на рис. 7 показывают координаты (27; ±8,5; ±8,5) неустойчивых положений равновесия. Для колебаний в окрестности положения равновесия частота с ростом их амплитуды меняется непрерывно

вплоть до момента, предшествующего завершению этих гладких колебаний. Как показывает спектральный анализ по методу Прони, двухпериодные сегменты этого численного решения характеризуются почти гармоническим спектром частот. Для временной зависимости г = г(Ь) численного решения модели Лоренца на рис. 8 показано отличие частот /[1],1 = 1 : 10 двух-периодного сегмента xJ[1 : Щ = :

^+2],^ = N • 3/2, с N = 56,3 = 14 от гармоник I = 1 : 10 фундаментальной частоты / = /14 [1] (сплошная линия). Интервал времени ДЬ8 = [20 : 21,5]. Моменты перескоков между неустойчивыми положениями равновесия могут быть определены непосредственно на рис. 7 по скачкообразному изменению переменной г = г(Ь). Например, через два периода после Ь = 21,5 при Ь & 23,5 на рис. 7 происходит перескок к другому неустойчивому положению равновесия модели (9).

10 20 30 40 50 60 i

Рис. 7. Координатно-временные зависимости модели (eq9)

"I 2 3 4 5 Б 7 В 910

Рис. 8. Отличие спектра частот от гармоник фундаментальной частоты /14 [1] при I = 1:10

Однако до момента перескока при Ь & 22,5 состояние колебаний в окрестности неустойчивого положения равновесия (27; —8,5; —8,5) становится нестабильным, и это сопровождается резким изменением частоты. Для значений N = 56,3 = 15 у двухпериодного сегмента хз [1 : Щ обнаруживается две фундаментальные частоты. Оказывается, что в момент времени Ь & 22,5 происходит изменение фундаментальной частоты колебания от & 1,510 до & 1,355. На рис. 9 для г = г(Ь) показано:

1 — отличие частот ¡15[1],1 = 1:2:9 и

2 — частот ¡15 [/],/ = 2 : 2 : 10 от гармоник I = 1:5 фундаментальных частот ¡15 [1] & 1,510 и ¡15 [2] & 1,355 соответственно (сплошные линии). При этом относительное изменение частоты составляет & 0,1. Это позволяет оценить порог, выше которого по методу Прони идентифицируется изменение частоты временного ряда. Для сравнения, во всем временном интервале непрерывного изменения фундаментальной частоты нестационарного колебания при сдвиге на N = 56 точек шкалы дискретизации времени относительное из-

менение этой частоты составляет

0,02.

1 1.5 2 2.5 3 3.5 4 4 5 5 1

Рис. 9. Отличие спектра частот: ^1) — ¡15[1],1 = 1:2:9 и («£) - ¡15[1],1 = 2 : 2 : 10 от гармоник I = 1:5 фундаментальных частот ¡15 [1] и ¡15 [2] соответственно

Скелетные кривые переходных

процессов. Как уже отмечалось, последовательность наполовину перекрывающихся двухпериодных сегментов нелинейных временных рядов хк[1 : N0],К = 1 : Р определяется следующим образом:

хзк [1 : N] = хк [пз : Пз+2], пз = N3/2, 3 = 0: К.

Пусть в интервале времени

АЬз = АЬ[пз : Пз+2]

изменение частоты нелинейного временного ряда происходит непрерывно. Для гладкого квазипериодического нелинейного временного ряда спектр частот, определённый для однопериодного сегмента, должен быть квазигармоническим. Поскольку фактически периодичность либо квазипериодичность временного ряда проявляется за два его последовательных периода, ниже рассматриваются спектры двухпери-одных сегментов. Именно такой спектр является квазигармоническим с фундаментальной частотой, усреднённой за два периода локальной орбиты (например, в пространстве состояний модели (9)). Оценка медленно меняющегося спектра частот ¡зк [1 : Мзк] локальной орбиты периода 1, усреднённого за два последовательных периода, имеет вид

¡зк[1] & I • ¡з, I = 1: Мзк, 3 =0: К,

где ¡з — фундаментальная частота локальной орбиты. Определим общий спектр частот синхронных сегментов хзк [1 : N],К = 1 : Р нелинейных временных рядов хк[1 : Щ],К = 1 : Р за период Тз = АЬ[п.з : пз+г] как

р

{¡з[/],/ = 1: Мз} = р| ¡зк[1 : Мзк]. к=1

(10)

Спектр частот (10) может быть определён на основе только одного из сегментов хзк[1 : Щ,К = 1 : Р с помощью специального критерия. Установлено, что для каждого сегмента хзк [1 : N] с фиксированным значением индекса К множество его полюсов гзк [1 : Мзк] разделяется на два подмножества {г,зк[1],1 £ !ж} и

[z3K[l],l Е IJK}, J U J = [1 : Mjk], в соответствии с критерием

max \6jk[1] \ < min \6jk[1]\-

r(N)

lei

(Q)

lei

(N)

Полюса {г,зк[1],1 £ Iзк} являются специфическими для сегмента хзк [1 : N] с фиксированным значением индекса К. Полюса {гзк[1],1 £ 1$»} характеризуют все сегменты хзк[1 : N],К = 1 : Р. Таким образом,

{¡з[/],/ = 1: Мз} = {¡зк[1],1 £ з}.

Частоты полюсов типа Q являются гармониками фундаментальной частоты /J [1] с той точностью, с которой получен сегмент XJK [1 : Щ. Пронумеруем полюса типа Q в порядке возрастания их частот:

{1,1 Е 1™} = {¡г,г =[1: МК]}. Тогда

Лк] & ФМ* = 1: М^К.

При фиксированном 3 точность, с которой полюса типа Q являются общими для всех синхронных сегментов хл^[1 : N],К = 1 : Р, зависит от этого индекса 3 = 0 : Л Преимущество метода Прони — это возможность выбора временных интервалов ДЬJ для спектрального анализа временного ряда хк[1 : N0],K = 1 : Р в соответствии с характерным масштабом NДЬ его изменения. Поэтому этот метод обеспечивает простую интерпретацию спектральных характеристик.

Пусть в моменты времени Ь/ = ДЬ • т/, I = 1 : Q фундаментальная частота нелинейного временного ряда меняется скачком. Определим Кх = топид(т/ /N). Зафиксируем индекс I. Тогда для 3 = [К/ : (К/+1 — 1)] (табличная) функциональная зависимость фундаментальной частоты / от индекса 3 является гладкой. Поэтому могут быть определены табличные амплитудно-частотные зависимости:

ак/,3 =[К/ : (Л/+1 — 1)]},

К = 1: Р. (11)

Зафиксируем индекс К. Если соответствующая табличная амплитудно-частотная зависимость (11) является достаточно гладкой, она частично представляет скелетную кривую а к = а к (/) для К-й переменной динамической модели. Тогда полное табличное представление этой скелетной кривой а = а (/) для всего временного ряда хк [1 : имеет вид

{аJK ,/J,3 = 1 : С} =

Я

= и ^К/,3 = [К/ : (Л/+1 — 1)]}, (12) /=1

где С = £?=1 ¡видЬк([К/ : (К/+1 — 1)]). Предполагается, что объединение множеств (12) отсортировано по возрастанию значения фундаментальной частоты

fj,3 = 1 : G. Отметим, что скелетная кривая ак = ак(f) действительно существует, если степень гладкости её полного табличного представления (12) выше или такая же, как степень гладкости всех её частичных табличных представлений (11). Численные эксперименты показывают, что для численных решений модели Лоренца скелетные кривые действительно существуют. Каждое гладкое табличное представление (12) скелетной кривой может быть аппроксимировано некоторым аналитическим выражением. Это аналитическое выражение используется для экстраполяции. Зафиксируем индекс K. В настоящей работе для аппроксимации обратной табличной зависимости (12) используется набор экспонент с вещественными показателями. Ошибка такой аппроксимации минимизируется с помощью алгоритма Нелдера-Мида [43]. Полученная аппроксимация f = f (ак), которая является обратной функцией для скелетной кривой а к = а к (f), используется вплоть до значения а к = 0.

Поскольку для численного решения модели Лоренца (9) при фиксированных параметрах данной модели скелетные кривые являются универсальными характеристиками всех непрерывных последовательностей локальных орбит в окрестности как первого, так и второго неустойчивого равновесия в пространстве состояний (x,y,z) модели (9), эти скелетные кривые являются универсальными характеристиками численного решения. Экстраполяция зависимостей f = f (а), обратных скелетным кривым а = а(f) на рис. 10, даёт три независимые оценки для фундаментальной частоты малой локальной неустойчивой орбиты периода 1 в окрестности неустойчивого положения равновесия. Соответствующие оценки для фундаментальной частоты: fox ~ 0,16151, f0y и 0,16152 и foz ~ 0,16153. Для оценки фундаментальной частоты fo малой локальной орбиты в окрестности одного из двух нестабильных положений равновесия (фиксированных точек) линеаризованные уравнения модели (9) могут быть решены для u = x — x0, v = y — y0 и w = z — z0, где Xo, yo и Zo — координаты неустойчивых фиксированных точек С = (27,6\/2,6\/2) и С' = (27,-6\/2,-6\/2). Для рас-

сматриваемых значений параметров в = 8/3,а = 10,р = 28 матрица А линеаризованных уравнений модели (9) имеет вид

A

-8/3 0 1

-10 -10 -6у/2 -1

Матрица собственных чисел матрицы A:

e-

-13,8546

0

+ j 10,1945

0 0

- j 10,1945

откуда фундаментальная частота бесконечно малой расходящейся орбиты /о = 1,6225. Оценки частот /хо ,/уо,/го, полученные в результате экстраполяции скелетных кривых, отличаются на величину & 0,0074, то есть относительная ошибка оценки фундаментальной частоты /о по скелетным кривым составляет & 0,3%. Метод определения скелетных кривых переходных процессов для трёх координатно-временных зависимостей модели Лоренца (рис. 7) составляет основу алгоритма разделения области переходных последовательностей неустойчивых локальных орбит порядка 1 и области неустойчивых глобальных орбит в пространстве состояний модели Лоренца (9).

I

125 1 3 1 35 1 4 1 45 1 5 1 55 1 6 1 65 '

Рис. 10. Скелетные кривые переходных процессов для координатно-временных зависимостей модели Лоренца: 1 — x = x(t); 2 — y = y(t); 3 — z = z(t)

Литература

1. Basseville M. Edge detection using sequential methods for change in level. Part ii: Sequential detection of change in mean

// IEEE Transactions on Acoustics, Speech, and Signal Processing. — 1981. — V. 29. — P. 32-50.

2. Basseville M., Benveniste A. Sequential detection of abrupt changes in spectral changes of digital signals // IEEE Transactions on Information Theory. 1983. — V. 29. — P. 709-724.

3. Basseville M., Nikiforov I.V. Detection of Abrupt Changes: Theory and Application. — Prentice-Hall, Inc, 1993.

4. Appel U., Brandt A.V. Adaptive sequential segmentation of piecewise stationary time-series // Information Science. — 1985. — V. 29, N.1. — P. 27-56.

5. Djuric P. Segmentation of nonstationary signals // Proceedings of ICASSP 92. — 1992. — V. 5. — P. 161-164.

6. Porat B., Friedlander B. Adaptive detection of transient signals // IEEE Transactions on Acoustics, Speech and Signal Processing. — 1986. — V. 34, N.6. — P. 1410-1418.

7. Porat B., Friedlander B. Performance analysis of a class of transient detection algorithms: a unified framework // IEEE Transactions on Signal Processing. — 1992. — V. 40, N. 10. — P. 2536-2545.

8. Thornburg H., Gouyon F. A flexible analysis-synthesis method for transients. — Stanford University: CCRMA, 2000.

9. Manders E.-J., Biswas G. Transient detection and analysis for diagnosis of abrupt faults in continuous dynamic systems // Proceedings of the IEEE International Workshop on Intelligent Signal Processing. — Budapest, Hungary: IEEE Computer Society, 2001. — P. 15-20.

10. Agrawal R., Faloutsos C., Swami A. Efficient similarity search in sequence database // Lecture Notes in Computer Science 730. — Springer-Verlag, 1993. — P. 69-84.

11. Prony R. Essai experimental et analytique // J. l'Ecole Polytechnique. — 1796. — V. 1, I. 2. — P. 24-79.

12. Pisarenko V.F. On the estimation of spectra by means of nonlinear function of covariance matrix // Geophysics Journal of Royal Astronomical Society. — 1970. — V. 28. — P. 511-531.

13. Pisarenko V.F. The retrieval of harmonics from a covariance function // Geophysics Journal of Royal

Astronomical.Society. — 1973. — V. 33. — P. 347-366.

14. Kumaresan R., Tufts D.W. Estimating the parameters of exponentially damped sinusoids and pole-zero modeling in noise // IEEE Transactions on Acoustics, Speech, and Signal Processing. — 1982. — V. ASSP-30, N.4. — P. 833-840.

15. Tufts D.W., Kumaresan R. Singular value decomposition and improved frequency estimation using linear prediction // IEEE Transactions on Acoustics, Speech, and Signal Processing. — 1982. — V. ASSP-30, N.4. — P. 671-675.

16. Hua Y., Sarkar T.K. Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise // IEEE Transactions on Acoustics, Speech, and Signal Processing. — 1990. — V. ASSP-38, N.5. — P. 814-824.

17. Hua Y, Sarkar T.K. On SVD for estimating generalized eigenvalues of singular matrix pencil in noise // IEEE Transactions on Acoustics, Speech, and Signal Processing. — 1991. — V. ASSP-39, N.4. — P. 892-900.

18. Sarkar T.K., Pereira O. Using the matrix pencil method to estimate the parameters of a sum of complex exponentials // IEEE Antenna and Propagation Magazine. — 1995. — V. 37, N.1.

P. 48-55.

19. Kukharenko B.G. Use of the Prony method for modal identification of slow-evolutionary linear structures // Journal of Structural Control. — 2000. — V. 7, N.2. -P. 203-218.

20. Кухаренко Б.Г. Технология спектрального анализа на основе быстрого преобразования Прони // Информационные технологии. — 2008. — № 4. — С. 38-42.

21. Lichtenberg A.J., Liberman A.P. Regular and stochastic motion / Applied Mathematical Science. — V. 38. — New York, Heidelberg, Berlin: Springer-Verlag, 1983.

22. Universality of Chaos: 2nd edition / ed. by Cvitanovic P. — Copenhagen, Denmark: Niels Bohr Institute, 1989.

23. Cvitanovic P., Garpard P., Schreiber T. Investigating of the Lorenz gas in term of periodic orbits // Chaos: An Interdisciplinary Journal of Nonlinear Science. — 1992. — V. 2, I.1. — P. 85-90.

24. Lim E, Mareels I.M.Y. Locating periodic orbits from experimental data //

Proceedings of 37th Conference on Decision and Control (CDC'98) / ed. by Castanon D. — Tampa, USA: IEEE Computer Society. 1998. P.2449-2454.

25. Dhamala M., Lai Y.-C., Kostelich E.J. Detecting unstable periodic orbits from transient chaotic time-series // Physical Review E. — 2000. — V. 61. — P. 6485-6489.

26. Skowronski J.M., Guttalu R.S. Realtime attractors // Dynamical Systems. —

1990. — V. 5, I.4. — P. 255-279. "

27. Cadzow J.A.Wilkes D.M. Enhanced rational signal modeling // Signal Processing. — 1991. — V. 25. — P. 171-182.

28. Vautard R, Yiou P., Ghil M. Singular-spectrum analysis: A toolkit for short, noisy, chaotic signals // Physica D. — 1992. — V. 58. — P. 95-103.

29. Golyandina N., Nekrutkin V., Zhigliavsky A. Analysis of Time Series Structure: SSA and Related Techniques. — Boca Raton, FL: Chapman & Hall/CRC, 2002.

30. Kukharenko B.G. On real-time estimation of system dynamics // Proceedings of the 1st International Conference «Control of Oscillation and Chaos»: Piscataway, NJ: IEEE Computer Society. — 1997. — V. 2. — P. 315-316.

31. Nicolis G., Nicolis C., McKernan D. Stochastic resonance in chaotic dynamics // Journal of Statistical Physics. — 1993. — V. 70, N.1/2. — P. 125-139.

32. Заславский Г.М., Сагдеев Р.З., Усиков Д.А., Черников А.А. Слабый хаос и квазирегулярные структуры. — М.: Наука,

1991.

33. Szemplincka-Stepnicka W., Rudovski J. Steady states in the twin-well potential oscillation. Computer simulations and approximate analytical studies // Chaos. — 1993. — V. 3. — P. 375-385.

34. Kukharenko B.G. Spectral identification of periodic trajectories for the model of interwell jumps in a bistable oscillator // Nonlinear Control Systems 2001 / ed. by Fradkov A.L., Kurzhanski A.B. — Oxford: Elsevier Science, 2002. — V. 2. — P. 951-956.

35. Kukharenko B.G. Use of the spanning Prony dimension under the similarity search of non-linear time-series // Proceedings of V international conference «System Identification and Control

Problems (SICPRO'06)». — Moscow: The V.A. Trapeznikov Institute of Control Sciences RAS, 2006. — P. 582-596.

36. Freidlin M.I. Quasi-deterministic approximation, metastability and stochastic resonance // Physica D. — 2000. — V. 37. -P. 333-352.

37. Kukharenko B.G. Detecting orbits of high periods from numerical solutions of the Rossler system // Proceedings of International Conference «Physics and Control» / ed. by A.L. Fradkov, A.N. Churilov. — St. Petersburg: IEEE Computer Society, 2005. — P. 866-871.

38. Farmer D., Crutchfield J., Froehing H., Packard N., Show R. Power spectra and mixing properties of strange attractors // Annals of N.Y. Academy of Science. — 1980. — V. 357. — P. 453-472.

39. Sparrow C. The Lorenz Equations. — New York: Springer-Verlag, 1982.

40. Kukharenko B.G. Numerical estimates of local and global motions

of the Lorenz attractor // Proceedings of International Conference «Physics and Control» / ed. by A.L. Fradkov, A.N. Churilov. — Washington, DC: IEEE Computer Society, 2003. — V. 2. P. 632-637.

41. Moskvina V., Zhigljavsky A.A. An algorithm based on Singular-Spectrum Analysis for change-point detection / / Communication in Statistics. Statistics and Simulations. — 2003. — V. 32. — P. 319-352.

42. Ide T., Tsuda K. Change-point detection using Krylov subspace learning // Proceedings of 2007 SIAM International Conference on Data Mining (SDM 07). -Omaha, NE, 2007. — P. 515-520.

43. Nelder J.A., Mead R. A simplex method for function minimization // Computer Journal. — 1965. — V. 7. — P. 308-313.

Поступила в редакцию 02.07.2009.