Содержание к диссертации
Введение
1. Использование переходных процессов для целей глобального моделирования 17
1.1.Введение 17
1.2.Методики оценки качества глобальных моделей при их реконструкции по различным участкам временного ряда 18
1.3. Примеры использования переходных процессов для улучшения качества глобальных моделей 20
1.3.1. Расширение области хорошей аппроксимации объекта моделью за счет использования переходного процесса (реконструкция уравнений неавтономного осциллятора Тода) 21
1.3.2. Потеря информации о структуре объекта при установлении движения (реконструкция дискретной многомодовой системы) 23
1.3.3. Улучшение качества аппроксимации при наличии «лишних» базисных функций (реконструкциая автономного осциллятора В ан-дер-Поля - Тода) 3 0
Использование переходного процесса для оптимизации набора базисных функций 31
2. Модификация процедуры выбора динамических переменных 38
2.1.Введение 38
2.2.Процедура тестирования набора переменных на возможность построения глобальной динамической модели 40
2.3. Модификация процедуры, возможность использования для тестирования на нелинейность 45
2.4.Примеры приложения процедуры тестирования 52
3. Оптимизация набора базисных функций 62
3.1.Введение 62
3.2. Оценка чувствительности значений коэффициентов перед базисными функциями к изменениям распределения тренировочных точек в фазовом пространстве 64
3.3 .Методика выбора базисных функций 71
3.4.Соотношение предложенного и ранее известного метода оптимизации 73
3.5.Тестовые примеры 75
4. Использование реконструкции модельных уравнений для анализа динамической стационарности временного ряда 81
4.1.Введение 81
4.2.Методы анализа динамической нестационарности 84
4.3. Анализ нестационарности при скачкообразном изменении параметров объекта: тестовые примеры 85
4.4.Использование простых моделей для описания сложных систем 94
4.4.1. Общие положения 94
4.4.2. Статистические оценки, к которым приводит построение динамических моделей 97
4.4.3. Быстрый метод обнаружения изменений многомерных распределений по скалярному временному ряду с помощью «плохих» моделей 99
4.5.Приложение к анализу внутричерепной ЭЭГ человека во время эпилептического припадка 103
Приложение. Моделирование автономных колебаний антенного жгутика москитов 120
Заключение 128
Список литературы
- Примеры использования переходных процессов для улучшения качества глобальных моделей
- Модификация процедуры, возможность использования для тестирования на нелинейность
- Оценка чувствительности значений коэффициентов перед базисными функциями к изменениям распределения тренировочных точек в фазовом пространстве
- Анализ нестационарности при скачкообразном изменении параметров объекта: тестовые примеры
Введение к работе
Возможности современной вычислительной техники позволили развить новые подходы к созданию эмпирических динамических моделей колебательных явлений. Если раньше речь шла в основном об аппроксимации наблюдаемых в эксперименте простых зависимостей, то сейчас - о получении модельных дифференциальных или разностных уравнений, описывающих сложные или даже хаотические движения. Повсеместное использование в современной измерительной технике аналого-цифровых преобразователей и становление концепции динамического хаоса привели к тому, что, начиная с 80-х годов особое внимание стало уделяться реконструкции1 уравнений по хаотическим временным рядам [1-7] -дискретной последовательностям чисел, полученным стробированием наблюдаемых величин. Временные ряды могут быть векторными или скалярными; последнее еще более осложняет задачу. С помощью реконструкции решаются задачи прогнозирования дальнейшего поведения [1] и бифуркаций [8], классификации систем [9,10], скрытой передачи информации [11]. В диссертации проблема построения моделей по временному ряду рассматривается на примерах эталонных для радиофизики динамических систем, а также для некоторых сложных биологических сигналов, однако она актуальна и во многих других областях исследований. В частности для задач астрофизики [12], лазерной физики [13], метеорологии [14,15], сейсмографии [16] финансов [17,18] и т.д.
Пик интереса в задаче глобальной реконструкции уравнений по временному ряду приходится на девяностые годы, после чего появились и
Здесь используется общепринятый термин "реконструкция" (восстановление), хотя он полностью адекватен лишь ситуации, когда временной ряд получен путем численного решения уравнений. В приложении к реальным объектам и явлениям, для которых нет единственной или "истинной" математической модели, уместнее было бы говорить о "конструировании", а не о "реконструкции".
2 Термин «глобальная» означает, что модельные уравнения, записанные в замкнутой форме, описывают поведение объекта во всем фазовом пространстве (глобально).
обобщающие материал работы [19-22]. Однако практическое применение разработанных методик выявило ряд принципиальных трудностей, возрастающих вместе с размерностью и степенью нелинейности моделируемых процессов. Это послужило причиной некоторого разочарования в возможностях культивируемых тогда универсальных подходов к глобальной реконструкции и продемонстрировало необходимость разработки направленных методик и приемов (технологий), учитывающих специфику достаточно узких классов моделируемых объектов. На разработку элементов таких технологий и направлена данная диссертационная работа. Суть проблем и необходимость постановки решаемых в работе задач поясним при описании процедуры реконструкции модели по временному ряду. При всем разнообразии подходов и практических ситуаций в ней можно выделить следующие основные этапы:
1-й этап: выбор типа модельных уравнений, например, разностные
x/+1=F(x,.) (1)
или дифференциальные
^W(x,). (2)
2-й этап: формирование ряда векторов состояния {х,}, где
х/ =(хн>"'>хш)> по РЯЛУ наблюдаемой величины {v,}/=1, которая скалярна, продеформирована измерительными приборами и искажена шумами. То есть выбрать количество D и способ получения динамических переменных. Наиболее часто употребляются методы временных задержек и последовательного дифференцирования, последовательного интегрирования [23], взвешенного суммирования [24], но в принципе, наблюдаемые и переменные могут быть связаны любым другим мыслимым способом, предугадать который без информации о принципах функционирования объекта весьма сложно. Неудачный выбор переменных затрудняет аппроксимацию неизвестных зависимостей F в модельных уравнениях выбранного вида (1-2) или вовсе делает их неоднозначными, а следовательно,
непригодными для динамического моделирования. Проблема оптимизации
выбора переменных уже привлекала к себе внимание [25-27], но известные
численные процедуры нуждаются в модернизации.
3-й этап: выбор вида аппроксимирующих функций F. В отсутствии
априорной информации или иных соображений о виде функции F для ее
представления используются различные универсальные формы. Обычно ее
выбирают в виде универсальной линейной комбинации некоторых базисных
функций:
м
F(x,) = |>4/4(x,), (3)
причем вид fk может оказать решающее влияние на результат
моделирования. Часто F берут в виде полинома, однако их использование (как и других универсальных аппроксиматоров) ведет к очень громоздким выражениям с большим числом лишних слагаемых, что является одной из причин неудач реконструкции в случае моделирования сложных систем. Удалив лишние базисные функции можно, сделав модель более компактной, существенно расширить область, в которой она работоспособна. Для этой цели известно несколько процедур, например, критерий Стьюдента или специальные процедуры [28,29]. Но их применимость ограничена, а расширение возможностей требует новых решений.
На этапах 1-3 фиксируется структура модели.
4-й этап: выбор тренировочного участка во временном ряде переменных, по которому затем оцениваются параметры модели (обычно методом наименьших квадратов). Задача выбора тренировочного ряда становится нетривиальной и требует решения в первую очередь при наличии в ряде нестационарности, а это типично, особенно для живых систем. Так неизвестно и требует изучения, следует ли включать в тренировочный ряд
Применимость этих методик обоснована, если неточности аппроксимации связаны только с присутствием шума. На практике типичны ситуации, когда источником ошибок кроме шума служит невозможность точного представления функций, стоящих в объекте, через выбранные базисные.
переходные процессы, как их учет повлияет на качество восстанавливаемых моделей. С другой стороны реконструкция параметров модели по нестационарному ряду может быть полезной для фиксации факта произошедших в объекте изменений (задача идеологически близкая прогнозу бифуркаций [8]). Информация о моментах изменения параметров объекта, выделение квазистационарных участков по хаотическим временным реализациям и их типологизация (классификация) очень важны, например, для решения задач нейрофизиологии.
5-й этап: оценка качества модели, после чего она используется или дорабатывается.
В работе рассмотрена возможность приложения разрабатываемых методик к анализу динамической нестационарности и выделению этапов протекания эпилептического приступа по записи ЭЭГ. Нестационарность ЭЭГ много анализировалась со статистических позиций - с точки зрения изменений статистических распределений и спектральных свойств (см., например, обзор [30]). Этапы протекания эпилептического приступа исследовались, в частности, с помощью спектрального анализа [31,32]. С динамических позиций этот вопрос никогда ранее не рассматривался, хотя, как оказалось, это позволяет получить дополнительную информацию по сравнению с результатами традиционных методик.
Предложенные и исследованные в работе приемы реконструкции колебательных моделей по наблюдаемым временным рядам отрабатывались на традиционных для радиофизики нелинейных динамических системах, в том числе при внесении в сигнал искажений и шума, апробировались на радиотехнических макетах и прилагались для задач физиологии и медицинской диагностики.
Цель работы
Разработка специальных подходов к выбору структуры уравнений и тренировочных участков, расширяющих возможности реконструкции
эмпирических динамических колебательных моделей. Приложение разработанного аппарата к исследованию реальных сигналов, а также к задачам медицинской диагностики.
Основные задачи, решаемые в работе
~ разработка метода оптимального для целей моделирования выбора тренировочного участка нестационарного временного ряда;
~ модернизация метода тестирования выбранных динамических переменных на соответствие предполагаемой структуре уравнений;
-- разработка методов оптимизации набора базисных функций в отсутствии априорной информации о структуре модели;
~ исследование возможностей приложения методов реконструкции модельных отображений для анализа структуры электроэнцефалограммы.
Научная новизна
продемонстрирована эффективность использования переходных процессов при построении моделей по временному ряду, предложен метод оптимизации набора базисных функций (удаления лишних) для аппроксимации функций в модели, основанный на использовании свойств переходного процесса;
модифицирована процедура тестирования набора реконструированных динамических переменных на возможность динамического описания (существование однозначной непрерывной зависимости между последовательными векторами переменных);
~ предложен новый метод выбора базисных функций для аппроксимации неизвестных зависимостей при построении модели по временному ряду;
~ показаны возможности приложения метода реконструкции модели по временному ряду к анализу динамической нестационарности. Метод приложен к анализу временной структуры эпилептического припадка на основе анализа стационарности внутричерепной ЭЭГ.
Достоверность полученных результатов
Достоверность полученных результатов основывается на воспроизводимости всех численных экспериментов, использовании отработанных численных методов, непротиворечивости с известными в литературе результатами, а также высокой точности совпадения данных получаемых на тестовых примерах, когда восстановление идет по временным рядам эталонных динамических систем.
Основные результаты и положения, выносимые на защиту
Включение в тренировочный ряд участков переходных процессов способствует получению глобальных моделей. При реконструкции модели конкретного установившегося движения, нацеленного на прогноз, нестационарная часть временного ряда в общем случае снижает прогностические возможности модели. Учет переходных процессов повышает эффективность предложенной процедуры опознания и удаления лишних коэффициентов многокомпонентных моделей по уровню их вариабельности при сдвиге окна реконструкции.
Модифицирован метод тестирования способа формирования динамических переменных модели по ряду наблюдаемой на соответствие заданной структуре глобальной динамической модели. Выбранный критерий качества переменных4 позволяет определять еще и наличие нелинейности в аппроксимируемых зависимостях и судить о недостаточности объема данных во временном ряде.
Предложен новый метод оптимизации реконструируемых многокомпонентных моделей, основанный на оценке зависимости значений коэффициента перед базисной функцией от изменения плотности распределения точек тренировочного ряда в фазовом пространстве. Метод
Характерный вид зависимости максимального разброса значений аппроксимируемой величины от размера окрестности точек в восстановленном фазовом пространстве.
дает лучшие результаты по сравнению с ранее известными методиками, если набор базисных функций неполон.
4. Расчет расстояний между векторами в пространстве коэффициентов моделей, построенных по коротким участкам временного ряда, при высокой точности аппроксимации позволяет обнаруживать и анализировать динамическую нестационарность (изменение модельного оператора эволюции). При недостаточной для динамического прогноза точности, реконструированные уравнения могут быть использованы для экспресс-анализа на статистическую нестационарность.
Научная и практическая значимость результатов
Результаты диссертационной работы развивают методы получения моделей по временным рядам, имеющие общедисциплинарное значение. Разработанные методики могут применяться для анализа практически важных сигналов. Динамическая модель, описывающая поведение объекта, может служить инструментом для контроля его состояния. Параметры реконструированной по ряду модели после дополнительного анализа могут быть связаны с не измеряемыми явно свойствами объекта. Отдельной практически важной задачей является построение моделей со структурой выбранной с учетом априорной информации о природе моделируемого объекта. В этом случае построение модели может служить основой для проверки различных гипотез о механизмах его функционирования.
На основе разработанных методик были написаны программы для анализа ЭЭГ, внедряемые в настоящее время в диагностическую и исследовательскую практику (9-я гор. больницы города Саратова и Институт высшей нервной деятельности и нейрофизиологии РАН). С их помощью исследуется возможность контроля эффективности действия противосудорожных лекарств, и устанавливаются механизмы, лежащие в основе некоторых видов эпилепсии.
Личный вклад автора
Автором произведено программирование и проведены численные исследования и эксперименты. Формулировка поставленных задач, выбор методов их решения, а также интерпретация полученных результатов произведена совместно с научным руководителем и соавторами.
Апробация работы и публикации
Основные результаты диссертации составили содержание докладов на следующих конференциях и школах
VII Международная школа «Хаотические автоколебания и образование структур». Саратов, Россия, 1-6 октября, 2004.
XII Научная школа «Нелинейные волны - 2004», конференция молодых ученых. Нижний Новгород, Россия, 29 февраля - 7 марта, 2004.
International symposium «Topical problems of nonlinear wave physics» (Международный симпозиум «Актуальные проблемы физики нелинейных волн»), Nizhny Novgorod, Russia, September 6-12, 2003.
The 5th European Congress on Epileptology, Madrid, 2002.
VI Научная конференция «Нелинейные колебания механических систем». 16-19 сентября, 2002.
XI Всероссийская научная школа «Нелинейные волны - 2002». Нижний Новгород, Россия, 2-9 марта, 2002.
VI Международная школа «Хаотические автоколебания и образование структур». Саратов, Россия, 2001.
Научная школа-конференция «Нелинейные дни в Саратове для молодых - 2000». Саратов, Россия, 16-20 октября, 2000.
II Международная конференция «Фундаментальные проблемы физики». Саратов, Россия, 2000.
The 6th International Symposium on Nonlinear Theory and its Applications NOLTA-2000 (VI Международный симпозиум по нелинейной теории и ее приложениям). Dresden, Germany, 2000.
- V Международная конференция «Нелинейные колебания механических систем». Нижний Новгород, Россия, 1999.
Работы были поддержаны грантами РФФИ (№99-02-17735, № 02-02-17578), части работы по этим проектам выполнялись как индивидуальные разделы, поддержанные молодежными грантами РФФИ (№02-02-06502, № 03-02-06860), а также Американским фондом гражданских исследований и развития для государств бывшего Советского Союза (грант № REC-006).
По теме диссертации было опубликовано 17 научных работ, из них 5 статей в реферируемых научных журналах, 1 статья в коллективной монографии, 11 тезисов и публикаций в сборниках трудов конференций.
Структура работы
Во введении обосновывается актуальность и практическая значимость рассматриваемых в работе проблем, формулируется цель работы, перечисляются основные задачи, формулируются положения, выносимые на защиту.
В первой главе рассматривается возможность, преимущества и недостатки использования переходных процессов при проведении глобальной реконструкции динамических моделей по временным рядам.
Для исследования роли переходных процессов использовались временные ряды, получаемые путем численного решения уравнений эталонных динамических колебательных систем (осциллятора Тода под гармоническим внешним воздействием, связанных квадратичных отображений, осциллятора Ван дер Поля - Тода). Уровень шума при этом незначителен, определяется лишь погрешностями численного метода и машинного округления. Все это позволяет уйти от решения сложных проблем связанных с выбором качественного вида модели и влиянием шума и сосредоточиться на решении поставленной задачи.
Процедура исследования заключалась в следующем. Фиксировалась некоторая ширина окна реконструкции (М точек): Хк = {Xj}"}*^, где т -номер начальной точки (при увеличении т окно сдвигается по временному
ряду в область установившихся движений). Модели восстанавливались при различных значениях т и сравнивались различным критериям качества (точность аппроксимации функций, стоящих в объекте, предсказательные возможности модели).
Выделены случаи, влияния учета переходного процесса на качество модели. В зависимости от цели моделирования (получение глобальной модели или прогноз установившихся движений) использование переходного процесса может как улучшать, так и ухудшать качество модели.
Для случая, когда аппроксимирующие функции содержат «лишние» слагаемые, предложена процедура их выделения и удаления из базиса. Эффективность использования этой процедуры повышается при моделировании по переходному процессу.
Во второй главе развивается и модернизируется метод выбора способа формирования динамических переменных по временному ряду, предложенные ранее в работе [17].
Чтобы реконструкция была возможна необходимо, чтобы между последовательными восстановленными векторами состояния существовала однозначная и непрерывная зависимость. Поэтому, решая, какой набор переменных выбрать, необходимо проверить наличие однозначной непрерывной зависимости между левыми и правыми частями восстанавливаемых уравнений. Для этого в [17] предлагалось разбить фазовое пространство из восстановленных переменных на ячейки размером S и рассчитав в каждой из них разброс є аппроксимируемых значений, стоящих в левых частях модельных уравнений, выбрать из них максимальный. Если зависимость этого максимального разброса єттх от 6
идет в 0 при уменьшении S, то можно утверждать, что однозначная и непрерывная зависимость существует и построение модели возможно.
В данной главе предложено заменить разброс в ячейках на разброс между точками, лежащими на расстоянии меньше 8. Это позволяет
устранить немонотонность зависимостей разброса єттх (8), что облегчает их интерпретацию. Кроме того, после введения такой модификации оказывается возможным использовать зависимость ^щах(^) Для тестирования исследуемой зависимости на нелинейность.
Проанализировано влияние конечности и дискретности имеющегося набора тренировочных точек на результаты тестирования. Выделен признак некорректной работы процедуры, вследствие нехватки тренировочных точек для разрешения структуры зависимости на малых масштабах.
В третьей главе предлагается, исследуется и апробируется новый метод оптимизации набора базисных функций. В его основе лежит предположение о том, что коэффициенты перед «лишними» базисными функциями будут чувствительны к малым изменениям распределения тренировочных точек. Эффекта аналогичного изменению этого распределения можно добиться вводя для каждой точки ее вес. При оценке значений коэффициентов по методу наименьших квадратов критерий их определения после введения весов будет выглядеть следующим образом:
( м Л
V k=\ J
Для случая ортогональных базисных функций получена аналитическая
оценка чувствительности коэффициентов к малым изменениям весовой
функции р —> р'' = р + р. В случае неортогонального набора достаточно
рассчитать эту чувствительность для проекции базисной функции, ортогональной всем остальным базисным функциям, так как только она в действительности улучшает качество аппроксимации.
Предлагается процедура совершенствования базиса, которая в данном случае состоит из выбора некоего начального универсального набора с последующим исключением тех функций, чувствительность которых к малым изменениям весовой функции тренировочных точек максимальна.
Показана связь новой процедуры оптимизации с ранее известной [18,19]. Процедура апробируется, и эффективность ее работы сравнивается с эффективностью ранее известной методики на ряде численных примеров.
В 4-й главе дан краткий обзор методов анализа нестационарности, включая динамическую (изменение оператора эволюции порождающего ряд объекта). Рассмотрена возможность применения для анализа динамической нестационарности методики построения моделей по временному ряду. На простых численных примерах рядов динамических систем со скачкообразным изменением параметра (логистическое отображение, отображение косинуса) показано, что точное определение момента произошедших изменений возможно, только если модель с высокой точностью описывает динамику системы. В противном случае ее коэффициенты оказываются зависящими от распределения тренировочных точек, и получаемая картина отражает нестационарность параметров этого распределения (то есть статистических, а не динамических свойств ряда).
Методика применена для анализа стационарности внутричерепных ЭЭГ во время эпилептического припадка. С ее помощью удается выделить несколько стадий в его протекании. Результаты анализа динамической нестационарности припадка с помощью глобальных моделей сравниваются с более традиционными методами - построением спектрограмм и вейвлет-спектров, а также результатами анализа стационарности одномерной функции распределения. Показано, что анализ динамической нестационарности дает информацию, дополняющую то, что получается с помощью перечисленных, ставших традиционными, методик.
В приложении представлены результаты моделирования автономных колебаний антенного жгутика москита, возникающих под воздействием определенного лекарства (DMSO), по временным рядам скорости его движения, предоставленным специалистами университета города Цюриха (Швейцария). Выбор в качестве переменных самой наблюдаемой величины
'/
xt =v(tt) и интеграла от нее с переменным верхним пределом у( - \v(t)dt
позволил получить восстановленную фазовую траекторию, не содержащую самопересечений. Для такого выбора переменных получена модель в виде осциллятора Ван дер Поля с нелинейным потенциалом, которая с высокой точностью воспроизводит динамику объекта.
Примеры использования переходных процессов для улучшения качества глобальных моделей
Исследование затрагивает лишь один из этапов построения модели по временному ряду - выбор места окна реконструкции в исходном ряде, обеспечивающего работоспособность модели в широкой области фазового пространства. Поэтому в качестве объектов используются эталонные дифференциальные уравнения и связанные отображения, и строится модель в соответствии с их структурой - она или полностью повторяет структуру объекта и при реконструкции требуется только оценить неизвестные коэффициенты, или при сохранении числа уравнений в модели используется другой вид аппроксимирующих функций. Исходные временные ряды (дискретные последовательности значений одной из динамических переменных объекта) получаются путем численного решения, при этом уровень шума незначителен - определяется лишь погрешностями численного метода и машинным округлением. Все это позволяет уйти от решения сложных проблем связанных с выбором качественного вида модели и влиянием шума и сосредоточиться на решении поставленной задаче5.
Ниже на ряде эталонных примеров показано, когда учет переходных процессов полезен для моделирования. Расширение области хорошей аппроксимации объекта моделью за счет использования переходного процесса (реконструкция уравнений неавтономного осциллятора Тода) Рассмотрим пример реконструкции дифференциальных уравнений. В качестве объекта используем неавтономный осциллятор Тода: = У (13) у = -0Л-у-\ + е х +cos(t). к }
Эта система представляет широкий класс нелинейных динамических систем, содержит экспоненциальную нелинейность и явную зависимость от времени. Скалярный временной ряд значений координаты х был получен путем численного интегрирования уравнений (1.3) методом Рунге-Кутта 4-го порядка с шагом 0.01. Начальные условия выбирались вне аттрактора, которым при данных значениях параметров является предельный цикл. Фазовая траектория, соответствующая исходному ряду показана на рис. 1.1, где цифрами приведены номера некоторых точек (начальная имеет номер 0). Аттрактор расположен в области, где точки траектории расположены плотнее. По рисунку видно, что движения вблизи аттрактора занимают незначительную часть площади, охватываемой траекторией.
Структура модельных уравнений была выбрана в виде (1.1) с функцией /: f{x,y,t) = Pn(x) + jy + Asin(t) + 5cos(0 (1.4) где Pn(x) — полином 5-й степени. Использовалась технология реконструкции, опирающаяся на стандартную изложенную во введении схему, модернизированную для неавтономных систем [34-36]. По временному ряду определялись значения параметров у, А, В и коэффициентов полинома (методом наименьших квадратов).
На рисунке 1.2 приведена зависимость нормированной ошибки аппроксимации (1.2) от длины тренировочного ряда. Видно, что существует оптимальная длина ряда (около 2000 точек). Дальнейшее увеличение этой длины, вовлекающее все больше точек с аттрактора, приводит только к ухудшению качества аппроксимации.
Чтобы сравнить качество моделей, реконструированных по переходному процессу и аттрактору, проводилась реконструкция в скользящем окне длиной 2000 тренировочных точек. На рисунке 1.3 приведена зависимость получившейся ошибки аппроксимации є от номера первой точки окна реконструкции т. Наилучшие модели получались при малых значениях т -при учете переходного процесса. Причем, график є(т) имеет характерный «ступенчатый» вид, иллюстрирующий ценность того или иного участка ряда для целей глобальной реконструкции. Так почти горизонтальные участки соответствуют выведению из тренировочного ряда точек из окрестности аттрактора, что свидетельствует об их малом влиянии на качество модели. Скачки между плоскими ступеньками на графике, на которых качество модели с ростом т быстро ухудшается, соответствуют исключению из тренировочного ряда точек, отмеченных на фазовой траектории (рис. 1.1) кружками. Эти точки принадлежат переходному процессу и наиболее удалены от аттрактора.
Для нелинейных диссипативных систем типична ситуация, когда аттрактор содержится в некотором подпространстве фазового пространства. В процессе установления колебаний трансверсальные этому подпространству возмущения затухают. Поэтому часть информации потенциально полезной для моделирования теряется. Внизу приведены увеличенные участки графиков в окрестности резких скачков ошибки аппроксимации. Область интегрирования для вычисления є такая же, как и для рис. 1.2. (однородного) режима в системе идентичных связанных маятников каждый элемент движется одинаково, как если бы других элементов не было - т.е. ряд не несет никакой информации о сложности системы. Временные ряды переходных процессов могут нести в себе информацию о затухающих неоднородных движениях и расширить возможности моделирования. Продемонстрируем это, используя в качестве объекта систему диссипативно связанных квадратичных отображений где х, у - динамические переменные, к — параметр связи, А - параметр неравновесности, п = 0,1,2,... - дискретное время. Эта эталонная система характеризуется бесконечным числом регулярных и хаотических видов колебаний (вариантов синхронизации колебаний подсистем) и мультистабильностью [37-44]. Выберем значения параметров Я=1,8, к = 0,4-0,5 при которых в системе устанавливается синфазный хаотический режим — когда переменные х, у меняются во времени хаотически, но одинаково (имеет место хаотическая синхронизация). На рис. 1.4 представлен фазовый портрет этой системы при параметре связи к = 0.4 и начальных условиях хо=ОЛ, у0=0А. Жирными точками показан переходный процесс, цифрами указан номера итераций. Видно, что движения подсистем быстро синхронизируются и трансверсальная аттрактору составляющая движения затухает. На рис. 1.5 показана временная реализация этой составляющей (хп —у„) для значений к=0Л (кружки) и к = 0.42 (крестики). Заметим, что меньшим значениям к из выделенного интервала соответствует большая длительность переходных процессов по сравнению со значениями, близкими к к=0.5. Это объясняется тем, что синфазный аттрактор при выбранном значении Я устойчив лишь при к 0.361 [40], а вблизи бифуркационных значений параметра возмущения затухают слабее.
Модификация процедуры, возможность использования для тестирования на нелинейность
Первая вторая и третья итерации для логистического отображения (2.1); (б) Графики максимального разброса єтйХ (S): быстрее осциллирующей зависимости соответствует более крутой график; (в) При наличии шума графики «подняты» вверх. Наиболее существенно усиливается разброс быстро осциллирующих зависимостей (соответствующие им графики «подняты» выше). x(t() и x(fy), и данный подход не применим, так как нет возможности исследовать локальные свойства зависимости у(х).
В экспериментальных данных неизбежно присутствуют погрешность измерений (связанная, например, с разрядностью АЦП) и шумы (влияние многочисленных факторов, которые не поддаются детерминированному описанию). Обозначим crnoise - их суммарный вклад в значения х и у. Когда становится меньше crnoise, больше не уменьшается даже при наличии определенной закономерности, связывающей у их. Так, если величина crnoise превышает размер участка большой крутизны, то график тах( 5) будет свидетельствовать о неоднозначности исследуемой зависимости (например, рис.2.2в, верхняя кривая). В этом случае согласно предлагаемому критерию рассматриваемые переменные не годятся для глобального моделирования.
Основное отличие описанной выше методики от 8-Е метода, предложенного в [45], в использовании локальных характеристик вместо средних (интегральных). В принципе, в качестве дополнительной характеристики можно использовать и величину среднего локального разброса?:
Если є — О при 8 — О и крутизна графика є(8) мала, то это может указывать на более «плавную» в среднем зависимость у(х), которую иногда (при прочих равных условиях) легче аппроксимировать гладкой функцией. В [27] показано, что сама по себе величина s не содержит достаточной информации для заключения о наличии однозначной непрерывной связи между у кх. Так, если Xх) имеет локализованный участок неоднозначности или разрыв, то они могут вносить малый вклад в среднее значение є и график є (8) будет выглядеть примерно так же, как и для гладкой однозначной зависимости.
Сходная задача о проверке наличия однозначной непрерывной связи между двумя временными рядами возникает в теории обобщенной хаотической синхронизации. Критерий ее проверки, сходный по своей идеологии с описанным здесь подходом предлагался в работе [47]. Отличие от меры однозначности, предложенной в [47], в первую очередь заключается в анализе зависимости єтах (8) во всем допустимом диапазоне значений 8, включая анализ асимптотического поведения при 8 — 0. В то время как в [47] предлагается использовать только одно значение єітх, соответствующее некоторому достаточно малому 8.
В предлагаемой процедуре на оценку пригодности набора переменных для моделирования влияет по крайне мере три величины (рис.2.3). Это, как в [47], значение максимального разброса при наименьшем возможном 8, наклон графика єітх (8) при малых 8 и пересечение линейного продолжения графика Emzx (8) с осью ординат.
Изложенная выше процедура тестирования - это методика с фиксированным набором ячеек. Она имеет тот недостаток, что вектор х, лежащий близко к границе ячейки, не сравнивается с близкими к нему векторами соседних ячеек, но может сравниваться с более удаленными векторами «своей ячейки» (см. рис. 2.4). При наличии шума это может привести к сильным колебаниям на графике smax (8) при малых 8 из-за того, что при одном значении 8 два вектора попадают в одну ячейку, а при другом, но близком, - в разные. Такая немонотонность зависимости затруднит оценку пригодности рассматриваемых переменных.
Участок восстановленного фазового пространства, разбитый на квадратные ячейки. Точками показаны полученные по временному ряду вектора переменных. При использовании немодифицированной методики (предложенной в [27]) рассчитывается величина разброса между значениями тестируемой зависимости в точках А и С, но не проверяется разброс между значениями в точках А и В. При уменьшении изменении размера ячеек точки А и В могут то попадать в одну ячейку, то не попадать. Это приводит к неоправданным колебаниям зависимости разброса от размера ячеек.
Пример такой ситуации имеет место на рис. 2.2в, где приведены результаты тестирования по зашумленным хаотическим рядам, генерируемым логистическим отображением (2.1). Наблюдаемая г]п =хп +%п, где „ -случайные величины, распределенные по гауссовскому закону со стандартным отклонение составляющим 1% от дисперсии ряда {хп}. (По ряду длиной 1000 точек тестировались зависимости, соответствующие первой, второй и третьей итерациям: rjn+l(7ju)t 7jn+2(?]n) и 7]п+і(т]п))
Этот недостаток методики можно устранить, если вместо максимального разброса при фиксированном наборе ячеек использовать максимальный разброс только между значениями у(х) в точках лежащих на расстоянии меньше заданного 5:
Оценка чувствительности значений коэффициентов перед базисными функциями к изменениям распределения тренировочных точек в фазовом пространстве
Традиционный универсальный подход к реконструкции уравнений по временным рядам не предполагает наличия какой-либо информации о структуре неизвестных функций, входящих в модельные уравнения. Для их аппроксимации используются различные универсальные формы. Например, в соответствии с теоремой Вейерштрасса любую непрерывную функцию можно аппроксимировать полиномом. На практике подобные модели зачастую оказываются неработоспособными - в подобном громоздком базисе почти всегда присутствуют лишние элементы. Строго говоря, лишней следует считать такую базисную функцию, коэффициент перед которой отличен от нуля только из-за присутствия шумов. Присутствие таких базисных функций немного уменьшает ошибку аппроксимации моделью тренировочного ряда, однако приводит к большим ошибкам во всех остальных областях восстановленного фазового пространства. Выбросив из базиса эти функции можно существенно расширить область фазового пространства, где построенная модель будет хорошо описывать поведение объекта.
Традиционный подход к устранению лишних базисных функций основан на том, что поскольку они описывают шумовую компоненту, их вклад в уменьшение ошибки аппроксимации тренировочного ряда относительно мал. Можно просто последовательно выбрасывать из базиса такие функции, избавление от которых приводит к наименьшему росту ошибки аппроксимации [28,29]. Оптимальным можно считать то количество базисных функций, которое будет соответствовать минимальной ошибке предсказания тестового ряда или какому-нибудь другому критерию качества модели.
Однако влияние шума не единственная причина, ограничивающая область фазового пространства, где модель дает хорошую аппроксимацию поведения объекта. Если набор базисных функций не полон, то есть не позволяет описать объект с абсолютной точностью, коэффициенты модели подгоняются таким образом, чтобы улучшить аппроксимацию в точках тренировочного ряда, а вдали от них аппроксимирующая функция может вести себя произвольным образом. Для получения в этом случае аппроксимации, работоспособной в максимально широкой области фазового пространства, в работе предлагается специализированная процедура оптимизации набора базисных функций.
Строго говоря, ни одну из базисных функций в этом случае нельзя назвать лишней, так как все они уменьшают ошибку аппроксимации детерминированной составляющей временного ряда. Однако нежелательными можно назвать те, которые уменьшают ее только в узкой области фазового пространства (в окрестности точек тренировочного ряда). Метод обнаружения таких функций основан на предположении о том, что коэффициенты перед такими функциями будут сильно зависеть от распределения точек в реконструированном фазовом пространстве. Когда коэффициенты перед базисными функциями будут слабо зависеть от распределения, то можно ожидать что такая аппроксимирующая функция будет описывать объект не только в точках тренировочного ряда, но и в некоторой их окрестности.
В этой главе получена оценка зависимости коэффициентов аппроксимирующей функции от распределения точек в фазовом пространстве. Затем предлагается на основе этой оценки исключать базисные функции с нестабильными коэффициентами. Эффективность такой оптимизации на тестовых примерах сравнивается с эффективностью традиционного алгоритма.
Пусть {х( ,-)}/=1 - ряд векторов состояния, сформированный по исходному временному ряду. Модельные уравнения могут быть реконструированы, например, в виде дифференциальных: i(0 = F(x(0,a), (3.1) или в виде разностных уравнений: x(f,+1) = F(x(f,),a), (3.2) где F — функции выбранного вида, зависящие от неизвестных коэффициентов а. Эти функции обычно ищутся в виде линейных комбинаций известных базисных функций gk (х): м F{x,d)= akgk(x). (3.3) k=i
Для расчета коэффициентов а обычно используют метод наименьших квадратов, то есть значения коэффициентов выбираются так чтобы минимизировать величину среднеквадратичной ошибки п N( М Л2 (3.4) 1=14 k=\ где у,- - полученные по временному ряду значения, стоящие в левой части модельных уравнений, например, производные векторов состояния x(t) в уравнениях (3.1) или следующий по времени вектор состояния х(//+1) в уравнениях (3.2).
Заметим, что при таком подходе результат аппроксимации будет зависеть от распределения р(х) векторов состояния х(ґ,) в реконструированном фазовом пространстве. Аппроксимация будет более точной в той области, где плотность точек наблюдаемого ряда наиболее высока. Эффекта, аналогичного изменению плотности распределения точек
Анализ нестационарности при скачкообразном изменении параметров объекта: тестовые примеры
Большинство сигналов в природе, и в особенности это относится к физиологическим данным, нестационарны, то есть меняют свои свойства с течением времени. К нестационарности, например, приводит наличие в системе процессов, характерный временной масштаб которых больше длительности интервала наблюдения, внешние воздействия, наличие переходных процессов и дрейфа параметров. Задачи фиксации самого факта нестационарности и разработки методов анализа таких ситуаций актуальны для многих областей исследований.
В теории случайных процессов (при статистическом подходе) дается определение стационарности в узком и в широком смыслах. Процесс называется стационарным в узком смысле, если все его многомерные плотности распределения вероятности не зависят от времени. То есть pn(xl,tl,...,xn,tn) = pn(xl,tl+T,...,xn,tn+T) справедливо для любых n,tltТ.
Поскольку по конечному временному ряду невозможно оценить все плотности вероятности, то для практических целей используют менее строгие определения. В частности процесс называют стационарным в широком смысле, если от времени не зависит его среднее, а корреляционная функция зависит только от разности времен 4Jx(tl,t2) = 4Jx(tl2). Таким образом, стационарность в широком смысле свидетельствует о постоянстве спектра мощности сигнала.
В последнее время в связи с развитием методов анализа временных рядов, основанных на теории нелинейных динамических систем, появилось и стало популярным еще одно понятие стационарности. Это динамическая стационарность, которая подразумевает постоянство оператора эволюции системы, генерирующей ряд [50-55]. Именно такой вид стационарности важен, когда для анализа временных рядов применяются методы, заимствованные из нелинейной динамики. Стационарность в широком смысле говорит только о неизменности линейных связей между точками ряда. Динамическую нестационарность можно рассматривать как обобщение этого понятия на произвольный вид связи.
Учитывая, что во временном ряде может быть как детерминированная, так и случайная компонента, динамическую стационарность иногда определяют как постоянство всех условных распределений вероятности Ро(хп+і \хп), где \хп} - последовательность векторов состояния,
восстановленных по временному ряду [55]. Заметим, что из стационарности в узком смысле в этом случае следует и динамическая стационарность. Обратное же неверно, и примером динамически стационарного ряда, который нестационарен в узком смысле может служить переходный процесс в динамической системе.
Важным примером нестационарного временного ряда является запись ЭЭГ. Изменения вида ЭЭГ происходят почти постоянно и могут быть связаны с переходами между различными физиологическими или эмоциональными состояниями [30,56-59]. Особенно заметные изменения имеют место во время эпилептического припадка [60]. Период времени, в течение которого статистические (в частности спектральные) свойства ЭЭГ можно считать неизменными, по разным оценкам составляет от 4 секунд до 1-й минуты [30,31,57,61].
Анализ нестационарности важен не только в приложениях к биологическим сигналам. Так в [8] рассматривается возможность предсказания бифуркаций по нестационарному временному ряду. Эта методика используется затем для прогноза появления озоновых дыр [63]. Моделирование по нестационарному ряду использовалось в [11] для скрытой передачи информации. В [50,64] рассматривалась проблема разделения временного ряда на квазистационарные участки с последующей их классификацией. В этом случае целью анализа может быть отслеживание изменения состояния генерирующей ряд системы. В данной главе рассмотрены возможности такого подхода при анализе нестационарности эпилептических припадков.
Статистическая нестационарность (в широком смысле) во время припадка известна и исследовалась, например, в [31] и [32]. В последней работе с помощь оконного преобразования Фурье специального вида [65] был проведен частотно-временной анализ ЭЭГ во время эпилептического припадка, были выделены несколько стадий его протекания, отличающихся спектральным составом.
В то же время на основе тестов с суррогатными данными было показано, что свойства ЭЭГ во время припадка, не могут быть полностью описаны только спектром мощности [66]. Присутствие нелинейности [56,66-69] а также существование почти периодических комплексов, позволяют предположить, что динамика ЭЭГ во время эпилептического припадка может быть описана с помощью низкоразмерной нелинейной детерминированной модели (см., например, [70,71]). Хотя это предположение до сих пор не получило полного эмпирического подтверждения, оно является основой для исследований многих научных групп. Проводимые ими работы направлены на предсказание эпилептических припадков [72-81], локализацию эпилептического фокуса [68,69] и анализ синхронизации в мозге [82-87]. В то же время вопросы о стационарности ЭЭГ с точки зрения нелинейной динамики практически не рассматривались. Можно указать только работы [53,88] где с помощью специальных статистических тестов было показано, что ЭЭГ невозможно рассматривать как динамически стационарный сигнал. Вопрос же о «структуре» динамической нестационарности (существовании и продолжительности квазистационарных участков) во время эпилептического припадка, насколько нам известно, никогда ранее не рассматривался.
В данной главе дан краткий обзор некоторых методов обнаружения и анализа динамической нестационарности. Особое внимание уделено методу, основанному на реконструкции глобальных динамических моделей по временным рядам. Работоспособность этого метода демонстрируется на тестовых примерах временных рядов, полученных от систем со скачкообразным изменением параметров. Затем методика применяется для анализа стационарности ЭЭГ во время эпилептического припадка.
Похожие диссертации на Анализ нестационарности и выбор структуры уравнений при моделировании по временным рядам
