Реконструкция уравнений динамики и диагностика взаимодействия нелинейных колебательных систем по временным рядам

ОГЛАВЛЕНИЕ
ВВЕДЕНИЕ...........................................................7
ЧАСТЬ I. Реконструкция уравнений динамики.........................25
Глава 1. Оценка параметров хаотических систем.....................29
1.1. Оценка параметров одномерных отображений...................30
1.1.1. Известные методы.......................................31
1.1.2. Оценка с использованием обратного отображения..........35
1.2. Оценка параметров в случае скрытых переменных..............47
1.2.1. Известные методы.......................................48
1.2.2. Модифицированный метод множественной стрельбы..........51
Основные результаты главы 1.....................................60
Глава 2. Восстановление нелинейных характеристик..................61
2.1. Оптимизация модельной функции..............................61
2.2. Моделирование автоколебаний аристы дрозофилы...............65
2.2.1. Объект исследования и данные...........................66
2.2.2. Предварительный анализ и выбор структуры модели........67
2.2.3. Результаты моделирования...............................69
2.2.4. Обсуждение.............................................73
2.3. Реконструкция уравнений неавтономных осцилляторов..........75
2.3.1. Особенности процедуры реконструкции....................76
2.3.2. Примеры эталонных осцилляторов.........................77
2.4. Восстановление характеристик нелинейных элементов электрических цепей.............................................79
2.4.1. Контур с переключаемыми конденсаторами.................79
2.4.2. Контур с полупроводниковым диодом......................82
Основные результаты главы 2.....................................88
Глава 3. Реконструкция в случае черного ящика.....................90
3.1. Известные методы реконструкции.............................90
3.1.1. Реконструкция фазовой траектории.......................90
2
3.1.2. Аппроксимация функций многих переменных..................97
3.2. Подбор динамических переменных...............................98
3.3. Описание внешних воздействий................................102
3.4. Исследование динамической нестационарности..................105
3.3.1. Стационарность и нестационарность, эталонный пример 105
3.4.2. Приложение для анализа электроэнцефалограмм.............109
Основные результаты главы 3......................................112
ЧАСТЬ II. Методы диагностики взаимодействия систем.................114
Глава 4. Обобщения причинности но Грейнджеру.......................116
4.1. Причинность по Грейнджеру и динамический эффект.............116
4.1.1. Причинность по Г рейнджеру..............................117
4.1.2. Трудности интерпретации.................................121
4.1.3. Динамический эффект воздействия.........................127
4.2. Оценка связей между радиофизическими автогенераторами 131
4.2.1. Описание экспериментальной установки....................131
4.2.2. Результаты оценки связей по временным рядам.............136
4.3. Оценка связей между системами с переключениями..............140
4.3.1. Описание методов........................................141
4.3.2. Результаты численных экспериментов......................145
Основные результаты главы 4......................................148
Глава 5. Оценка связей между осцилляторами
по коротким временным рядам...............................150
5.1. Метод моделирования фазовой динамики........................151
5.1.1. Определение фазы........................................152
5.1.2. Теоретические характеристики связи......................153
5.1.3. Оценки характеристик связи..............................155
5.2. Трудности в случае коротких рядов...........................155
5.3. Новые оценки связи для коротких рядов.......................161
5.3.1. Вывод формул............................................161
3
5.3.2. Иллюстрации в численном эксперименте....................164
5.4. Пределы применимости в зависимости от свойств осцилляторов.. 168
5.4.1. Методика исследования...................................169
5.4.2. Влияние свойств шума....................................170
5.4.3. Влияние индивидуальной нелинейности осцилляторов 171
5.4.4. Влияние силы связи между осцилляторами..................173
5.4.5. Влияние погрешностей при расчете фаз....................174
5.5. Пределы применимости в зависимости от длины ряда.............177
5.5.1. Объект и методика исследования..........................177
5.5.2. Результаты численных экспериментов......................178
Основные результаты главы 5.......................................181
Глава 6. Сравнительный анализ методов оценки связей.................183
6.1. Моделирование фазовой динамики
и анализ в пространствах состояний............................183
6.1.1. Метод анализа в пространствах состояний.................183
6.1.2. Тестовые объекты........................................186
6.1.3. Методика исследования...................................188
6.1.4. Результаты численных экспериментов......................189
6.2. Моделирование фазовой динамики
и частная направленная когерентность..........................197
6.2.1. Метод оценки частной направленной когерентности.........197
6.2.2. Тестовые объекты........................................199
6.2.3. Методика исследования...................................205
6.2.4. Результаты численных экспериментов......................207
Основные результаты главы 6.......................................219
Глава 7. Оценка связей в ансамблях осцилляторов.....................220
7.1. Выявление структуры связей в ансамблях осцилляторов..........220
7.1.1. Описание метода ........................................222
7.1.2. Результаты численных экспериментов......................229
4
7.2. Оценка нелинейных связей между осцилляторами................243
7.2.1. Описание метода.........................................243
7.2.2. Результаты численных экспериментов......................248
7.3. Интервальная оценка времени запаздывания связи;.............259
7.3.1. Точечная оценка времени запаздывания....................260
7.3.2. Интервальная оценка времени запаздывания................262
7.3.3. Результаты численных экспериментов......................265
Основные результаты главы 7......................................268
ЧАСТЬ III. Приложения методов диагностики взаимодействии 270
Глава 8. Приложения в нейрофизиологии..............................270
8.1. Взаимодействие между глубокими структурами мозга
и конечностями при паркинсоновском треморе...................271
8.1.1. Экспериментальные данные................................272
8.1.2. Оценки связей на основе моделирования фазовой динамики ..275
8.1.3. Оценки связей на основе причинности по Грейнджеру.......283
8.2. Взаимодействие между структурами мозга крыс - генетических моделей абсанс-эпилепсии.........................................291
8.2.1. Экспериментальные данные................................292
8.2.2. Оценки связей на основе причинности по Грейнджеру 293
Основные результаты главы 8......................................297
Глава 9. Приложения в климатологии.................................299
9.1. Оценки воздействия различных факторов на глобальную
приповерхностную температуру.................................300
9.1.1. Используемые данные.....................................300
9.1.2. Индивидуальная модель вариаций ГПТ......................302
9.1.3. Модель вариаций ГПТ с учетом солнечной активности.......305
9.1.4. Модель вариаций ГПТ с учетом вулканической активности...307
9.1.5. Модель вариаций ГПТ с учетом содержания ССЬ.............308
9.1.6. Четырехкомпонентная модель вариаций ГПТ.................309
5
9.2. Оценки взаимодействия процессов Эль-Ниньо/Южное колебание
и Северо-Атлантическое колебание.........................312
9.2.1. Используемые данные.................................313
9.2.2. Результаты анализа..................................314
9.3. Оценки взаимодействия между Эль-Ниньо в Тихом океане
и его аналогом в экваториальной Атлантике................323
9.3.1. Используемые данные.................................324
9.3.2. Результаты анализа..................................325
9.4. Оценки взаимодействия процесса Эль-Ниньо/Южное колебание
и индийского муссона.....................................334
9.4.1. Используемые данные.................................335
9.4.2. Результаты анализа................................ 337
Основные результаты главы 9..................................347
ЗАКЛЮЧЕНИЕ.....................................................349
БЛАГОДАРНОСТИ..................................................352
ЛИТЕРАТУРА.....................................................353
6
Введение
За почти вековую историю своего развития радиофизика и се идейная база - теория колебаний — постоянно расширяли круг объектов исследования и приложений своих методов. Если сначала речь шла о способах генерации и преобразования регулярных сигналов для радиосвязи и локации, то в последние десятилетия в центре внимания находятся сложные и хаотические сигналы, нелинейные колебательные процессы в системах различной природы [1-13]. Модели и методы радиофизики активно используются в различных областях - от физики до биологии, медицины и наук о Земле.
Отличительной чертой современного периода исследования колебаний является также цифровая форма представления и обработки информации: сигналы на выходе большинства современных измерительных приборов имеют вид временных рядов - дискретных последовательностей значений наблюдаемых величин. Подходы и идеи теории колебаний, нелинейной динамики, статистической радиофизики оказываются особенно плодотворными при анализе сигналов, демонстрирующих колебательный характер, нерегулярность, признаки нелинейности. В этом круге задач выделяется проблема реконструкции уравнений динамики (построения математической модели) по временным рядам, поскольку при ее успешном решении полученные уравнения могут использоваться для целого ряда приложений - от прогноза [14-17] до диагностики взаимодействия исследуемых систем [18-23]. Последнее очень востребовано в физике и химии [18,19], кардиологии [20,21], нейрофизиологии [22], климатологии [23]. Этими обстоятельствами обусловлена важность и актуальность темы диссертации, что подробнее обосновывается ниже.
Актуальность работы. Построение математических моделей по временным рядам получило название «идентификации систем» в математической статистике [24] и «реконструкции динамических систем» - в нелиней-
7
ной динамике [25-27]. Предшественницами современных задач реконструкции были задачи аппроксимации, и статистического исследования зависимостей между наблюдаемыми величинами. Долгое время наблюдаемые процессы моделировались с помощью явных функций времени 77 = /(О, аппроксимирующих множество экспериментальных точек на плоскости (/, 77), где / -время, 77 — наблюдаемая. Целью моделирования был прогноз будущего развития процесса или сглаживание зашумленных данных. В начале XX века серьезный шаг в развитии методов эмпирического моделирования нерегулярных процессов был сделан в математической статистике, когда было предложено использовать линейные стохастические модели [28]. Этот подход практически не имел альтернатив в течение полувека (1920-е - 1970-е) и нашел многочисленные приложения, особенно для прогноза и автоматического управления [29]. Успехи нелинейной динамики, обосновавшей возможность хаотического поведения нелинейных динамических систем малой размерности, и развитие вычислительной техники открыли перспективы успешного моделирования сложных процессов на основе нелинейных разностных и дифференциальных уравнений [30,31].
Обсуждая современное состояние проблемы, воспользуемся сложившейся типовой схемой реконструкции уравнений движения, которая включает в себя несколько этапов. На первом получают временной ряд; на втором -выбирают структуру модельных уравнений, т.е. все, кроме конкретных значений параметров; на третьем - оценивают параметры (подгонка модели); в итоге проверяют, удовлетворительно ли полученная модель описывает наблюдаемый процесс. Предложены теоретические идеи, обосновывающие алгоритмы действия на каждом этапе. При выборе структуры уравнений — это теоремы Такенса для реконструкции вектора состояния [32] и обобщенная аппроксимационная теорема для задания оператора эволюции [33]; при оценке параметров - минимизация различных целевых функций; при проверке
адекватности модели - расчет метрических и топологических характеристик
8
аттрактора [34] и т.д. Тем не менее, на практике получить удовлетворительную модель с помощью'существующих универсальных подходов зачастую не удается. Одним из широко известных препятствий является так называемое «проклятие размерности» - резко возрастающие с ростом размерности модели трудности аппроксимации и требования к объему и качеству данных. Но и кроме этого на каждом этапе процедуры моделирования остаются' нерешенные проблемы, связанные с оценкой параметров хаотических систем, оптимизацией структуры модели, выбором динамических переменных. Сложившаяся ситуация требует развития практически эффективных методов для реконструкции уравнений нелинейных колебательных систем в реалистичных условиях коротких временных рядов, хаотичности наблюдаемых процессов, нестационарности и зашумленности. Одним из перспективных направлений представляется разработка методов, ориентированных на избранные классы систем и задач, учитывающих особенности этих классов при выборе структуры модели, расширяющих возможности физической интерпретации результатов моделирования.
Среди практических приложений реконструкции в настоящее время выделяется обнаружение и оценка связей (диагностика взаимодействия) между нелинейными системами по данным наблюдений. Эта задача является обратной по отношению к исследованию динамики в ансамблях нелинейных колебательных сйстем. Во многих работах [35-39] показано, что динамика существенно зависит от структуры связей между элементами ансамбля. В частности, эти связи определяют виды синхронизации, простоту или сложность динамики, образование различных пространственно-временных структур. Поэтому оценка связей дает важную информацию при исследовании различных объектов, а методы решения этой задачи представляют значительный практический интерес. При этом особенное внимание уделяется возможности получения оценок по коротким и нестационарным сигналам, широко распространенным в физике, биологии, науках о Земле.
9
Существуют «непосредственные» (не опирающиеся на построение моделей) методы выявления зависимостей между процессами, включая корреляционный и спектральный анализ, теоретико-информационные и нелинейно-динамические характеристики. Однако для оценки направленных связей, т.е. для ответа на вопрос о том, влияет ли один процесс на другой и с какой «силой», наиболее подходящим инструментом оказывается реконструкция уравнений. С 1960-х гг. широко применяется линейная оценка «причинности по Грейнджеру», основанная на построении линейных авторегрессионных моделей и расчете улучшения прогноза одного процесса при учете в модели данных о другом процессе [40]. Она позволяет с надежностью установить факт наличия связей между линейными системами, но сталкивается с трудностями при решении более сложных задач. Во-первых, физическая интерпретация ее количественных характеристик не очевидна: не ясно, в какой степени они отражают влияние связей на те или иные свойства наблюдаемой динамики. Во-вторых, при анализе нелинейных систем линейная оценка причинности по Грейнджеру часто оказывается не пригодной даже для выявления связей. Попытки ее нелинейных обобщений, использующие реконструкцию уравнений в универсальном виде, сталкиваются с упомянутыми выше общими трудностями нелинейного моделирования по временным рядам и, как следствие, с недостоверностью выводов о наличии связей. Более перспективный нелинейный метод оценки связей, предложенный для исследования колебательных систем, допускающих введение фазы, основан на моделировании фазовой динамики [41]. Однако его известный вариант ориентирован на случай длинных стационарных сигналов. Для коротких или нестационарных временных рядов получаемые с его помощью результаты не надежны, т.к. не исследованы статистические свойства итоговых оценок связи и не предусмотрена оценка статистической значимости выводов. Таким образом, для исследования взаимодействий между сложными колебательными процессами различной природы на практике требуются новые подходы, позво-
10
ляющие получать физически интерпретируемые характеристики связей и применимые для исследования нелинейных систем в реалистичных постановках достаточно коротких и зашумленных временных рядов, нелинейных и запаздывающих связей, многих взаимодействующих систем.
Наконец, следует отметить, что теоретические идеи нелинейной динамики, на которых основаны реконструкция динамических систем, диагностика связей, решение ряда других задач, и соответствующие методы обычно тестируются на эталонных динамических системах в численном эксперименте. При этом количество работ, посвященных их приложениям для анализа реальных процессов относительно мало. Поэтому в последние годы все настойчивее интерес научного сообщества к практическим приложениям идей и методов теории колебаний и нелинейной динамики в биофизике, физиологии, науках о Земле. Таким образом, актуальность темы диссертации связана с необходимостью развития практически эффективных методов для реконструкции уравнений и диагностики взаимодействия нелинейных систем в условиях дефицита данных и сложности наблюдаемой динамики, а также их востребованностью при исследовании нелинейных колебательных процессов различной природы.
Нелыо диссертационной работы является развитие методов реконструкции уравнений и диагностики взаимодействия нелинейных колебательных систем в условиях дефицита данных и сложности наблюдаемой динамики и их приложения для анализа колебательных процессов различной природы.
Для достижения цели решались следующие основные задачи.
1. Разработка комплекса методов, повышающих эффективность процедуры реконструкции уравнений динамики по временным рядам, включая оценку параметров хаотических систем при дефиците данных, оптимизацию модельных функций, подбор динамических переменных. Апробация развитых методов при моделировании радиотехнических и биологических систем.
2. Разработка и апробация комплекса методов для диагностики взаимодействия в ансамблях нелинейных колебательных систем, ориентированного на получение достоверных выводов о наличии связей и физически интерпретируемых характеристик связи по относительно коротким временным рядам.
3. Приложения развитых методов диагностики взаимодействия для исследования колебательных процессов в реальных системах различной природы, включая радиотехнические, нейрофизиологические и климатические.
На защиту выносятся следующие положении н результаты.
1) Разработанные методы оценки параметров хаотических динамических систем по зашумленным временным рядам дают более точные результаты, чем известные подходы. А именно, модифицированный метод множественной стрельбы, допускающий разрывы фазовой траектории на интервале наблюдения, позволяет эффективно использовать сколь угодно длинные ряды и смягчает требования к стартовым догадкам для искомых параметров. Для одномерных отображений использование обратного отображения при расчете целевой функции дает оценки параметров, среднеквадратическая- погрешность которых в типичном случае уменьшается с ростом длины временного ряда N как 1//У в отличие от закона 1/4Й для подходов, основанных на сегментировании ряда.
2) Разработанный комплекс методов для реконструкции уравнений динамики расширяет возможности моделирования нелинейных колебательных систем по временным рядам. Он включает в себя приемы подбора динамических переменных на основе тестирования аппроксимируемых зависимостей на однозначность и непрерывность, оптимизации модельных уравнений за счет исключения лишних слагаемых, описания внешних воздействий за счет использования многочленов с переменными коэффициентами. По сравнению с известными подходами это позволяет получать модели с меньшим числом динамических переменных, воспроизводящие наблюдаемую динамику в более широкой области фазового пространства.
12
3) Предложенный метод оценки динамического эффекта воздействий по временным рядам позволяет количественно охарактеризовать, в какой степени различные свойства одного процесса зависят от других наблюдаемых процессов (факторов). Он основан на построении эмпирической модели и анализе ее динамики при искусственных изменениях рассматриваемых факторов. Это дополняет широко используемые характеристики причинности по Грейнджеру, которые позволяют выявить наличие связей между исследуемыми системами, но не дают возможности оценить степень влияния этих связей на динамику.
4) Предложенный модифицированный метод моделирования фазовой динамики, основанный на учете корреляционных свойств фазовых шумов, позволяет с заданной доверительной вероятностью делать выводы о наличии связей между двумя нелинейными колебательными системами по временным рядам длиной от двадцати характерных периодов колебаний. Метод становится более чувствительным к слабой связи, чем известные подходы (оценка частной направленной когерентности и статистика ближайших соседей в пространствах состояний), при уменьшении коэффициентов диффузии фазы исследуемых систем и длины временного ряда.
5) Предложенные обобщения метода моделирования фазовой динамики позволяют выявлять структуру связей в ансамблях колебательных систем, получать физически интерпретируемые характеристики взаимодействий, оценивать связи, характеризующиеся нелинейностью произвольно высокого порядка, получать интервальные оценки времени запаздывания связей.
6) По эмпирическим данным за период 1856 - 2005 гг. выявлено влияние вариаций солнечной и вулканической активности и содержания углекислого газа в атмосфере на вариации глобальной приповерхностной температуры (ГИТ) с помощью оценки причинности по Грейнджеру. Согласно оценке динамического эффекта воздействий рост ГПТ в последние десятилетия объяс-
13
нястся эмпирической моделью только при учете в ней вариаций содержания
ео2.
7) По эмпирическим данным за период с 1870 г. до начала XXI в. с помощью метода моделирования фазовой динамики и оценки причинности по Грейнджеру выявлены связи процесса Эль-Ниньо — Южное колебание (ЭНЮК) в Тихом океане с процессами в других регионах. А именно, обнаружены и количественно охарактеризованы воздействие ЭНЮК на Северо-Атлантическое колебание, воздействие экваториальной атлантической моды на ЭНЮК и двунаправленная связь между ЭНЮК и индийским муссоном.
8) При анализе записей локальных потенциалов с глубинных электродов и сигналов.акселерометров с помощью метода моделирования фазовой динамики и линейной и нелинейной оценок причинности но Грейнджеру выявлена двунаправленная связь между активностью субталамического ядра (структуры мозга в базальных ганглиях) и колебаниями противоположной руки при паркинсоновском треморе. Влияние колебаний руки на активность субталамического ядра обнаруживается и линейным, и нелинейными методами. Обратное воздействие выявляется только нелинейными методами и характеризуется запаздыванием от 200 до 400 миллисекунд.
Достоверность научных выводов обусловлена теоретическим обоснованием разработанных методов реконструкции уравнений и оценки связей с позиций нелинейной динамики и математической статистики, тестированием методов на эталонных системах в численных экспериментах и установлением эмпирических критериев их применимости, согласованием результатов численных расчетов и физических экспериментов, совпадением ряда результатов с результатами других авторов. При анализе физических, физиологических и климатических процессов достоверность подтверждается также проверкой адекватности эмпирических моделей, выполнением критериев применимости методов, высокой статистической значимостью выводов согласно аналитическим оценкам и воспроизводимостью результатов.
• Научная новизна результатов работы состоит в следующем.
Предложен оригинальный метод оценки параметров хаотических одномерных отображений по зашумленным временным рядам. Впервые показано, что модифицированный метод множественной стрельбы позволяет эффективно использовать для оценки параметров сколь угодно длинные хаотические ряды. Предложен и апробирован комплекс методов, повышающих эффективность реконструкции уравнений динамики нелинейных колебательных систем по временным рядам по сравнению с известными подходами.
Предложен метод оценки динамического эффекта воздействий различных факторов на исследуемый процесс, дополняющий широко используемую идею причинности по Грейнджеру. Разработан оригинальный комплекс методов для оценки связей между нелинейными колебательными системами в условиях относительно коротких временных рядов и существенных шумов.
На основе анализа направленных связей по эмпирическим данным выявлено и охарактеризовано воздействие различных факторов на вариации глобальной приповерхностной температуры. По эмпирическим данным выявлены связи процесса Эль-Ниньо/Южное колебание с другими крупномасштабными климатическими процессами: Северо-Атлантическим колебанием, экваториальной атлантической модой и индийским муссоном.
Впервые по эмпирическим данным выявлено воздействие субталамиче-ского ядра на колебания конечностей у пациентов с паркинсоновским тремором и получена оценка времени запаздывания этого воздействия. Показано, что в начале пик-волнового разряда у крыс линии WAG/Rij связь между ядрами таламуса и лобной корой (двунаправленная и асимметричная) резко усиливается. Показано, что автономные колебания аристы Drosophila melanogaster, вызванные введением диметилсульфоксида, адекватно описываются с помощью уравнений автогенератора с одной степенью свободы, и получены характеристики диссипации и возвращающей силы для этой системы.
15
Научное н практическое значение результатов работы.
В распространенной на практике ситуации, когда структура модели полностью известна из физических соображений, а все параметры и переменные имеют физический смысл, но не могут быть измерены, предложенные методы позволяют получить оценки параметров хаотических систем и восстановить временные ряды скрытых переменных.
В ситуации, когда структура модельных уравнений- отчасти известна, а неизвестны входящие в них функции (характеристики объекта), которые не могут быть непосредственно измерены, предложенный метод оптимизации структуры уравнений позволяет восстановить эти характеристики по временным рядам. Он апробирован при восстановлении эквивалентных характеристик нелинейных элементов электрических цепей. Получен патент на изобретение.
В ситуации, когда структура модели неизвестна, предложенный комплекс методов для реконструкции уравнений позволяет получать модели, точнее воспроизводящие наблюдаемую динамику в широкой области фазового пространства по сравнению с известными подходами, что имеет универсальное значение для моделирования колебательных процессов различной природы.
Метод оценки динамического эффекта воздействий позволяет установить, в какой степени те или иные свойства исследуемого процесса обусловлены различными факторами, т.е. охарактеризовать важность воздействий и предсказать изменения в наблюдаемой динамике при вариации упомянутых факторов.
Предложенный комплекс методов оценки связей между нелинейными колебательными системами позволяет решать практически важные задачи диагностики структуры связей в ансамбле и оценки нелинейности и запаздывания связей по относительно коротким временным рядам с контролируемой
16
. надежностью.-Последнее обстоятельство-важно-при нестационарности сигналов и дефиците данных, что типично в нейрофизиологии и геофизике.
Модели, автономных колебании аристы Drosophila melanogaster и полученные нелинейные характеристики значимы для биофизических исследований. Поскольку многие свойства спонтанных отоакустических излучений насекомых и позвоночных аналогичны, эти результаты должны быть, востребованы при дальнейших исследованиях физики усилителя в улитке уха позвоночных.
Выявленное с помощью реконструкции уравнений деление эпилептического припадка на квазистационарные сегменты может стать основой для дальнейшего нейрофизиологического анализа, направленного на изучение и классификацию таких сегментов при различных формах эпилепсии.
Результаты анализа связей между колебаниями конечностей и активностью подкорковых структур мозга во время спонтанного паркинсоновского тремора важны для дальнейшего развития нейрофизиологических представлений о его механизме. Результаты соответствуют гипотезе о том, что исследуемые структуры мозга играют активную роль в генерации тремора. Однако представления о петле обратной связи, содержащей только линии проведения нервных сигналов в обоих направлениях, по-видимому, следует обновить.
Результаты анализа климатических данных расширяют сведения о характере связей между крупномасштабными климатическими процессами, что важно в теории климата для усовершенствования и исследования моделей земной климатической системы и построения прогнозов.
Разработанные в диссертации методы и программы используются при проведении научных исследований в Институте физики атмосферы им. А.М. Обухова РАН (г. Москва). Результаты работы внедрены в учебный процесс на факультете нано- и биомедицинских технологий Саратовского госунивер-| ситета. Они составили основу монографий [340, 341], содержащих образова-
тельную компоненту, предназначенную для студентов и аспирантов.
!
t,
«
£ ‘
I
Личный вклад соискателя. Соискатель лично сформулировал, рассмотренные в диссертации задачи, разработал новые методы реконструкции уравнений и оценки связей и компьютерные программы для- их реализации, провел тестирование и. сравнительный анализ этих методов, выполнил основную часть работ по анализу экспериментальных данных. Выбор направления исследований и интерпретация ряда результатов осуществлялись совместно с проф. Б.П. Безручко. Некоторые результаты по тестированию методов и анализу экспериментальных данных получены совместно с Т.В. Ди-каневым, М.Б. Бодровым, И.В. Сысоевым, A.C. Караваевым, B.C. Власкиным, С.С. Козленко, ГТ.И. Наконечным. Лабораторные радиотехнические макеты для тестирования развитых методов в физическом эксперименте разработаны Е.П. Селезневым и В.И. Пономаренко. Постановки задач анализа нейрофизиологических и климатических данных и интерпретация соответствующих результатов- осуществлены совместно со специалистами в этих областях (И.И. Мохов, P. Tass, Е.Ю. Ситникова, Е. van Luijtelaar, R. Stoop, M. Goepfert, J.-L. Perez Velazquez, R. Wennberg).
Апробация работы н публикации. Основные результаты диссертации составили содержание докладов на всероссийских школах «Нелинейные волны» (Нижний Новгород, 2002; 2004; 2006; 2008; 2010); международных школах «Хаотические автоколебания и образование структур» - ХАОС (Саратов, 1998; 2001; 2004; 2007); международных симпозиумах «Topical problems of nonlinear wave physics» (Нижний Новгород, 2003; 2005; 2008), «Topical problems of Biophotonics» (Нижний Новгород, 2007; 2009), «Nonlinear Theory and its Applications» - NOLTA (Дрезден, Германия, 2000); международных семинарах «Nonlinear Dynamics of Electronic Systems» - NDES (Борнхольм, Дания, 1999; Дслфт, Нидерланды, 2001; Рапперсвиль, Швейцария, 2009), «WE Heraues Seminar» (Бад Хонеф, Германия, 2001), «PASCAL» (Лавин, Швейцария, 2005); международных конференциях «Нелинейные колебания механических систем» (Нижний Новгород, 1999; 2002; 2005; 2008), «Фундамснталь-
18
ные проблемы физики» (Саратов, 2000), «Современные проблемы электроники СВЧ и радиофизики» (Саратов, 2001), «Synchronization of chaotic and stochastic oscillations» (Саратов, 2002), European Congress on Epilcptology (Мадрид, Испания, 2002), IEEE on Circuits and Systems for Communications (Москва, 2004), «Dynamic Days» (Пальма-де-Майорка, Испания, 2004), «Идентификация систем и проблемы управления» - SICPRO (Москва, 2005), «Biomagnetism» (Ванкувер, Канада, 2006), «Nonlinear Dynamics, Chaos, and Applications» (Меллас, Крым, Украина, 2006); научно-практических конференциях «Системный анализ в проектировании и управлении» (Санкт-Петербург, 2004), «Новые технологии в экспериментальной биологии и медицине» (Ростов-на-Дону, 2007); научно-технических конференциях «Радиотехника и связь» (Саратов, 2005), «Физика и радиоэлектроника в медицине и экологии» (Суздаль, 2006); всероссийских симпозиумах «Медленные колебательные процессы в организме человека. Теоретические и прикладные аспекты нелинейной динамики в физиологии и медицине» (Новокузнецк, 2005; 2007); русско-японском семинаре по нейробиологии и нейродинамике «Bridging non-linear dynamics with cellular and molecular neuroscience» (Токио, Япония, 2008); всероссийских школах-семинарах «Волновые явления в неоднородных средах» (Звенигород, 2008; 2009); научных конференциях «Нано-электроника, нанофотоника и нелинейная физика» (Саратов, 2006; 2007; 2008; 2009) и «Состав атмосферы. Атмосферное электричество. Климатические процессы» (Борок, 2008); школах-семинарах «Методы компьютерной диагностики в биологии и медицине» (Саратов, 2007; 2008; 2009); конкурсах работ молодых ученых им. И.В. Анисимкина (ИРЭ им. В.А. Котельникова РАН, Москва, 2005; 2006; 2007).
Результаты работы многократно обсуждались на научных семинарах Саратовского филиала Института радиотехники и электроники им. В.А. Котельникова РАН, Института физики атмосферы им. А.М. Обухова (г. Москва), Института медицины и Института вычислений им. Дж. фон Неймана Ис-
следовательского центра Юлих (Германия), Центра по анализу данных и моделированию Фрайбургского университета (Германия), научной группы нелинейной динамики Потсдамского университета (Германия), факультета эпилептологии Боннского университета (Германия), Института познания и информации Наймегенского университета (Нидерланды), университета Торонто (Канада), кафедры динамического моделирования и биомедицинской инженерии Саратовского государственного университета.
Работы были поддержаны грантом Президента РФ для молодых кандидатов наук, Российским фондом фундаментальных исследований, Фондом содействия отечественной науке, Американским фондом гражданских исследований и развития для государств бывшего Советского Союза (CRDF), Немецким научным фондом (DFG), ФЦПТП «Исследования и разработки по приоритетным направлениям развития науки и техники на 2002-2006 годы», программами РАН и Министерства образования и науки РФ.
По теме диссертации опубликовано 67 работ (без учета тезисов докладов), в том числе 2 монографии, 40 статей в журналах, рекомендованных ВАК, 5 статей в прочих журналах, монографиях и научных сборниках, 1 патент, 15 статей в сборниках трудов научных конференций, 4 учебнометодических пособия [340-406].
Структура н объем диссертации. Диссертация состоит из введения, девяти содержательных глав, сгруппированных в три части, заключения и списка литературы. Объем диссертации - 390 страниц, включая 111 рисунков, 5 таблиц и список литературы из 406 наименований.
Во введении обосновывается актуальность рассматриваемых в работе проблем, определяются цели исследования, ставятся основные задачи, формулируются положения и результаты, выносимые на защиту, раскрывается научная новизна и научно-практическое значение полученных результатов, их достоверность и личный вклад соискателя, кратко описывается содержание работы.
20
В первой части, состоящей из трех глав, разрабатываются и апробируются методы реконструкции уравнений динамики. Представлена общая схема процедуры моделирования, дана систематизация постановок задач по объему априорной информации о структуре уравнений, обсуждается проблема возможной некорректности задач. Распределение материала первой части по главам подчиняется принципу последовательного усложнения ситуации, т.е. сокращения объема априорной информации об объекте.
В первой главе рассматривается наиболее простая постановка задачи, когда структура модельных уравнений полностью известна и процедура реконструкции начинается с этапа оценки параметров. В разделе 1.1 предложен метод, позволяющий повысить точность оценок параметров одномерных отображений х„+1 = /(;с„,с). В разделе 1.2 рассматривается более общая и практически важная ситуация, когда оценка параметров осложняется тем, что некоторые из компонент вектора состояния могут быть не наблюдаемыми (скрытые переменные). Предложена модификация метода множественной стрельбы, расширяющая возможности оценок параметров по хаотическим временным рядам.
Во второй главе рассматривается более сложная постановка задачи, когда С'груктура уравнений частично известна, но в них входят неизвестные функции (компоненты функции f в модельных уравнениях dx/dt = t‘(x)), а не только неизвестные параметры. В разделе 2.1 предложен метод оптимизации функции f, основанный на исключении лишних слагаемых. В разделе 2.2 проведено моделирование динамики биологического объекта (колебания аристы в органе слуха Drosophila melanogaster). В разделе 2.3 процедура восстановления модельных функций адаптирована для случая гармонически возбуждаемых систем. В разделе 2.4 она апробируется в физическом эксперименте на неавтономных электрических контурах для получения эквивалентных характеристик нелинейных элементов в режимах больших амплитуд и хаоса.
21
В третьей главе рассматривается наиболее сложная постановка задачи моделирования, когда информация о структуре модельных уравнений отсутствует. В разделе 3.1 представлен обзор известных методов, включающий в себя изложение теорем Такенса и методов реконструкции фазовой траектории, способов аппроксимации многомерных функций и других элементов процедуры моделирования. Обсуждаются трудности моделирования при использовании известных универсальных подходов. В разделе 3.2 предложен метод подбора динамических переменных модели, основанный на тестировании аппроксимируемых модельных зависимостей на однозначность и непрерывность. В разделе 3.3 показана возможность описания внешних воздействий за счет использования многочлена с переменными коэффициентами в модельных уравнениях. В разделе 3.4 представлено приложение методов реконструкции для исследования динамической нестационарности: проведена сегментация электроэнцефалограмм пациентов с височной эпилепсией во время эпилептического припадка.
Во второй части диссертации развиваются методы диагностики взаимодействия между исследуемыми системами на основе реконструкции уравнений. Дана систематизация известных методов оценки зависимостей между процессами по данным наблюдений. Большую группу составляют «непосредственные» методы, т.е. не опирающиеся на построение моделей. Обсуждаются трудности их использования для выявления направленных связей. Обосновывается продуктивность «опосредованных» методов, основанных на реконструкции уравнений.
В четвертой главе предлагается и апробируется метод, дополняющий широко используемую идею «причинности по Грейнджеру», которая, как показано в работе, сталкивается с трудностями при анализе нелинейных систем и физической интерпретации результатов (раздел 4.1). В разделе 4.2 идея причинности по Грейнджеру реализована при оценке связей между нелинейными системами (радиотехническими генераторами) в физическом экспери-
22
менте. В разделе 4.3 предложен метод оценки связей между системами с переключениями, основанный на учете характерных свойства динамики этих нелинейных систем.
В пятой главе развивается метод диагностики, взаимодействия между нелинейными колебательными системами по коротким временным рядам. В разделе 5.1 описан известный подход, основанный на моделировании фазовой динамики (МФД). В разделе 5.2 показано аналитически и в численных экспериментах, что при небольшой длине ряда достоверно определить характер связи с помощью известного метода невозможно, т.к. оценки связей не снабжены доверительными интервалами и оказываются смещенными. В разделе 5.3 предложен модифицированный метод, основанный на учете корреляционных свойств фазовых шумов, который позволяет устранить смещение оценок связи и получить доверительные интервалы. В разделе 5.4 исследуются пределы применимости предложенных интервальных оценок связи при изменении свойств шума и фазовой нелинейности на примерах эталонных осцилляторов. В разделе 5.5 исследуется применимость предложенных оценок при уменьшении длины ряда, что важно при практическом применении к нестационарным нейрофизиологическим и коротким климатическим сигналам.
В шестой главе проводится сравнительный анализ предложенного выше метода, основанного на моделировании фазовой динамики с другими методами оценки связей, которые широко используются на практике и часто оказываются эффективными. Методы сравниваются по их чувствительности (способности обнаружения существующих слабых связей) в численном эксперименте. Сравнение с нелинейным анализом ближайших соседей в пространствах состояний представлено в разделе 6.1, с линейной оценкой частной направленной когерентности - в разделе 6.2.
В седьмой главе метод оценки связей между колебательными системами, развитый в главе 5 для случая двух исследуемых систем, связи между ко-
23
торыми не запаздывающие и описываются функциями с нелинейностью невысокого порядка, обобщается на более широкий круг ситуаций. В разделе 7.1 предложен метод диагностики, связей в ансамбле из М осцилляторов, где М> 2. В разделе 7.2 предложен метод, позволяющий использовать функции связи с произвольно высоким порядком нелинейности в модельных уравнениях. В разделе 7.3 рассмотрен случай запаздывающих связей и на основе известной точечной оценки времени запаздывания предложена соответствующая интервальная оценка, что важно при анализе коротких рядов для получения надежного заключения о наличии ненулевого запаздывания.
В третьей части диссертации развитые методы диагностики взаимодействия нелинейных систем применяются для решения практически важных задач нейрофизиологии и климатологии.
В восьмой главе представлены нейрофизиологические приложения. В разделе 8.1 исследуется взаимодействие между активностью глубоких структур мозга и колебаниями конечностей при паркинсоновском треморе. В разделе 8.2 исследуется взаимодействие между различными структурами мозга (лобная кора, таламус) крыс - генетических моделей абсанс-эпилепсии до, во время, и после пик-волнового (эпилептического) разряда.
В девятой главе представлены климатологические приложения развитых методов. В разделе 9.1 исследуется влияние солнечной и вулканической активности и содержание СОг в атмосфере на вариации глобальной приповерхностной температуры (ГПТ). В разделах 9.2 - 9.4 исследуется взаимодействие между крупномасштабными процессами в азиатско-тихоокеанском регионе и северной и экваториальной Атлантике. А именно, проводится анализ связей процесса Эль-Ниньо/Южное колебание с Северо-Атлантическим колебанием в разделе 9.2, с экваториальной Атлантической модой в разделе 9.3, с индийским муссоном в разделе 9.4.
Основные выводы и результаты работы представлены в заключении.
24
Часть I: Реконструкция.уравненийдинамики
В первой части,, состоящей из трех глав, разрабатываются и апробируются методы реконструкции уравнений динамики нелинейных колебательных систем по временным рядам. Во вводном разделе (ниже) представлена схема процедуры моделирования и систематизация постановок задач по объему априорной информации о структуре уравнений, обсуждается проблема возможной некорректности задач. Распределение материала первой части по главам подчиняется принципу последовательного усложнения ситуации, т.е. сокращения объема априорной информации об объекте.
Типовая схема процедуры моделирования по временным рядам. Несмотря на безграничное число ситуаций, объектов и целей, вносящих в процесс свое специфическое, можно выделить основные этапы процедуры моделирования, показанные в виде схемы на рис.1. Работа начинается с рассмотрения имеющейся информации об объекте (экспериментальных данных о по-добных объектах; теорий, разработанных для описания исследуемого класса объектов и т.д.), получения и предварительного анализа рядов наблюдаемых величин, а заканчивается использованием полученной модели для решения конкретной задачи. Этот процесс сопровождается неоднократными повторениями, последовательными приближениями к «хорошей» модели.
На ключевом в схеме этапе 2 формируется структура модели. В теории идентификации этот этап называют структурной идентификацией [42]. Согласно [24] «структура модели - это параметризованное множество моделей». Во-первых, выбирается тип уравнений. В диссертационной работе речь идет о моделях в виде отображений х„+1 =Г(х„,с) или обыкновенных дифференциальных уравнений с/х/б// = Г(х,с). Во-вторых, задается вид входящих в них функций (компонент функции I). В-третьих, устанавливается связь динамических переменных (компонент вектора х) с наблюдаемыми величинами ц. В качестве переменных могут выступать сами наблюдаемые, но в общем
/ \
модель удовлетворительна не удовлетворительна
У
Использование модели
Рис.1. Типовая схема процедуры реконструкции уравнений (построения математической модели) по временному ряду.
26
случае эту связь задают в виде т] = Ь(х), где Ь называют измерительной функцией. Чтобы учесть измерительный шум, вводят случайную добавку £: Я = Ь(х) + ?. Чтобы сделать модель еще более реалистичной, случайную добавку вводят и в сами уравнения - так называемый динамический шум.
Формирование структуры модели - это наиболее сложный этап. После него остается только определить значения параметров с (этап 3). Здесь говорят об оценке параметров, «подгонке» модели, или параметрической / ненараметрической идентификации [42]. Для этого, как правило, проводится поиск экстремального значения некоторой целевой функции, например, минимизируется сумма квадратов отклонений наблюдаемых данных от решения модельных уравнений. При необходимости проводятся предварительные преобразования наблюдаемого ряда: фильтрация, численное дифференцирование или интегрирование и т.п.
В финале (этап 4) проводится проверка адекватности (верификация) модели в отношении интересующих исследователя свойств объекта или пригодности модели для решения поставленной практической задачи. Если модель признана удовлетворительной, то процедура заканчивается, иначе модель возвращается на доработку на любой из этапов схемы рис.1.
Особенности задач эмпирического моделирования. Задачи построения модельных уравнений по наблюдаемой временной реализации принадлежат к классу обратных задач и часто являются некорректно поставленными, т.е. решение не существует, не единственно или не непрерывно зависит от входных данных [43]. В практике моделирования по временному ряду некорректность постановки задачи проявляется наиболее явно на этапе выбора структуры модели: множество (достаточно сложных) моделей с различной структурой могут одинаково хорошо воспроизвести наблюдаемый конечный набор данных. Эта некорректность чаще всего снимается за счет дополнительных предположений или априорной информации о структуре модели,
27
что позволяет сузить класс моделей, в котором ищется-решение задачи. Этот подход широко используется и в диссертационной работе. Он представляет собой вариант построения квази-решения - известного подхода из теории некорректных задач [43].
В зависимости от происхождения временных рядов складываются две качественно различные ситуации при постановке задачи моделирования. Первая - когда наблюдения представляют собой решение каких-либо уравнений, полученное численными методами. Именно для нее полностью уместен термин «реконструкция» уравнений. Здесь можно сформулировать и теоретические условия эффективности методов реконструкции для различных классов систем. Во втором, качественно ином случае, когда временной ряд получен в результате измерений реального процесса, «единственно правильной» модели не существует и успех моделирования гарантировать нельзя. Использование слова «реконструкция» здесь не вполне обосновано и больше подходит термин «построение модели». По в качестве жаргона термин «реконструкция» употребляется и в отношении практических задач.
Систематизация задач по объему априорной информации. Фон на рис.1 меняется от черного до белого, что отражает степень априорной неопределенности. Наименее благоприятна ситуация, получившая название «черный ящик», когда информация о подходящей структуре модели отсутствует, и начинать приходится с самого верха схемы. Чем больше известно о том, как должна выглядеть модель, т.е. чем ниже «стартовая точка» на схеме, тем вероятнее успех. С практической точки зрения удобно упорядочить различные постановки задачи моделирования по степени осведомленности исследователя о подходящей структуре модели.
1) Структура модельных уравнений полностью записана из физических соображений, неизвестны только значения параметров. Здесь могут возникать существенные трудности, связанные с большим числом неизвестных параметров и скрытых переменных. Эта постановка рассматривается в главе 1.
2) Структура уравнений в значительной степени известна из физических соображений. Неизвестны только некоторые компоненты функции Г.которые-могут представлять собой важные физические характеристики объекта. Этой постановке посвящена глава 2.
3) Реконструкция, в случае «черного ящика» - наиболее тяжелая ситуация. Разработка и апробация практически эффективных методов реконструкции в.такой постановке проводится в главе 3.
Глава 1. Оценка параметров хаотических систем
Здесь рассматривается наиболее простая постановка, когда структура модельных уравнений полностью известна и процедура реконструкции начинается с этапа оценки параметров. Временной ряд переменной-^ - набор значений Г1^х\г1^2 )...,77(7л.) - моменты наблюдений, N - длина ряда)
- генерируется дискретным отображением
х«.1=г(х».с) (1)
или системой обыкновенных дифференциальных уравнений
с1х/А = {(х, с), (2)
где х - £>-мерный вектор состояния, f - вектор-функция, с - Р-мерный вектор параметров, п - дискретное время, / - непрерывное время. Наблюдаемая переменная 7 =/г(х) + £*, где к - измерительная функция, С - измерительный шум (шум наблюдений). Временной ряд получен при значении с = с0, соответствующем хаотическому режиму. Нужно получить статистическую оценку д величины с0.
Рассматриваемая постановка имеет самостоятельную ценность. Как правило, в этой постановке все переменные и параметры модели имеют физический смысл, а прямого пути их измерения часто не существует, так что про-
цедура реконструкции выступает как замена измерительного прибора. Такие задачи встречаются в радиофизике (моделирование электрических цепей с
29
сегнетоэлектрическими или полупроводниковыми нелинейностями [44,45]), лазерной физике (оценка скоростей перехода между уровнями рабочего вещества газового лазера [3]), аэродинамике (оценка коэффициентов аэродинамических сил при летном эксперименте с самолетом [46]); клеточной биологии (оценка параметров-модели нейрона [47], модели сигнального пути клеток [48])..
Несмотря на относительную простоту, даже в этой постановке возникают трудности, связанные с использованием« хаотических временных рядов: такие сигналы богаты информацией о системе, так что теоретически можно ожидать получения очень точных, оценок параметров [49], но для практического достижения этой цели необходимы соответствующие методы, которые и развиваются в диссертационной работе. Раздел 1.1 посвящен повышению точности оценок [49-54] в специальном случае одномерных отображений. В разделе 1.2 рассматривается ситуация дефицита данных, когда по имеющемуся ряду наблюдаемой р не удается сформировать ряды всех динамических переменных хк,к — 1,т.е. некоторые переменные являются «скрытыми» [55-59].
1.1. Оценка параметров одномерных отображений
Одномерные отображения описывают сложную динамику многих реальных систем [60,61]. Так, отображение Пуанкаре сильно диссипативных нелинейных потоковых систем часто с высокой точностью одномерно (см., например, системы Лоренца и Ресслера [8]). Поэтому задача оценки параметров одномерного отображения специально рассматривалась в ряде работ [50-54]. В простом варианте задача сводится к проведению кривой через экспериментальные точки на плоскости (лп>т]}М). Есть и другие методики (см п.1.1.1), но общей особенностью практически используемых вариантов является то, что среднеквадратическая погрешность оценки параметров уменьшается с ростом длины ряда не быстрее, чем 1/л//У.
30
Особенностью хаотической динамики является высокая чувствительность фазовой траектории модели к значениям параметров. Это позволяет надеяться на получение оценок, погрешность которых гораздо быстрее уменьшается с ростом длины ряда. Такое повышение точности действительно проиллюстрировано на нескольких эталонных примерах, наблюдалось даже экспоненциальное уменьшение дисперсии оценок [49,54]. Однако это имеет место только, в очень узком диапазоне длин ряда. Практически же получить такие оценки невозможно из-за принципиальных трудностей с поиском глобального минимума целевой функции.
В п. 1.1.2 предложен метод повышения точности оценок, при котором их среднеквадратическая погрешность уменьшается с ростом от длины ряда в типичном случае как 1/Лг. Он основан оценке параметров обратного отображения, что ранее применялось для нелинейной фильтрации шума [62] и разделения суммы хаотических сигналов на составляющие [63].
1.1.1. Известные методы
Методы иллюстрируются на примере оценки параметра квадратичного отображения
Хп+1 = /(хп>°),
Лп ~ Хп £п’
Л
где функция /'(*„ ус) = 1 — с • хп, значение параметра с = с0 соответствует хаотическому режиму, С, - нормальный белый шум с дисперсией <7*.
В отсутствие измерительного шума (о^ =0) имеет место Лп~хп и экс~ периментальные точки на плоскости (г/п,т]п+,) лежат точно на квадратичной параболе (рис. 1.1,а). Поиск величины с0 сводится к решению алгебраического уравнения. Итоговая оценка выражается как с = (1 “*„+,)/*« (здесь и далее оценки, полученные по временному ряду, обозначаются «крышечками»).
1.2 - ^ а) р > 1.2 - б) 1.2 •
0.4 / \ / * 0.4 } Ч / \ 0.4
-0.4 $ X % 9 О -0.4 * 1 ? -0.4
-1.2 - 1“ »■ 1 Т ■ -т Г-Х ' Ч—<—1 -1.2 - Л —■ Г—1—г —і—■—1—■—I 1 1 -1.2
п„+, о„ч, в)
*4''' •*>
ьо Ло
/ V
•» А
■ *© \°
О*. О
л„
—1---:--1--;---1--1--1---1--1---і--1--.
-1.2 -0.4 0.4 1.2 -1.2 -0.4 0.4 1.2
-1.2 -0.4 0.4 1.2
Рис. 1.1. Оценка параметров на примере квадратичного отображения (1.3) при с*о = 1.85, кружки - наблюдаемые значения: а) нет шума, пунктир — исходная парабола; б) только динамический шум, пунктир - модельная парабола, полученная минимизацией среднего квадрата вертикальных расстояний (некоторые показаны жирными линиями); в) только измерительный шум, пунктир - модельная парабола, полученная минимизацией среднего квадрата ортогональных расстояний; г) только измерительный шум, ромбики — реализация модели, наиболее близкая к наблюдаемому ряду в смысле (1.6).
32
Достаточно использовать любые два измерения хп,хп+{ при хп *0. Оценка 6 совпадает с истинным значением с0 с точностью до погрешности вычислений. При наличии шумов вместо точного значения-параметра ищут его статистическую оценку.
Динамический шум, измерительного шума нет. В методических целях рассмотрим для иллюстрации известных методов случай отсутствия измерительного шума (сг^=0) и наличия динамического- шума £
дгя+1 = 1 -С'Х2п + %п, где £п - последовательность случайных величин, статистически независимых и одинаково распределенных с плотностью /?(£). Для оценки параметра с часто используют метод максимального правдоподобия (ММП), который наиболее эффективен при достаточно общих условиях [64,49]. Функция правдоподобия имеет вид:
Лг-1
1пДс) = 1пр(Л1,Г1г Х>/^(т7„+1 ~/(Я„’<■'))■ (1-4)
И=|
Чаще всего рассматривают нормально распределенный шум. Тогда максимизация (1.4) эквивалентна методу наименьших квадратов (МНК), т.е. минимизации суммы квадратов невязок
£(с) = £(77„+1 -/(я„,с))2 ->гшп. (1.5)
л=!
Геометрически это означает, что график модельной функции на плоскости (т7;1,77,1+|) должен пройти так, чтобы сумма квадратов вертикальных расстояний от экспериментальных точек до него была минимальна (рис. 1.1,6). В примере (1.3) функция /линейна по с, поэтому целевая функция Л' квадратична по с и имеет единственный глобальный минимум, который легко находится решением линейного алгебраического уравнения.
Ошибка оценки д уменьшается с ростом длины ряда N. В рассматриваемой постановке задачи МНК дает асимптотически несмещенные и состоятельные оценки. Дисперсия оценок убывает как #-1, т.к. слагаемые в (1.5)
33
стационарны по и, т.е. частные информационные количества Фишера ограничены [49] (см: «асимптотические свойства оценок» в п. 1.1.2).
Измерительный шум* динамического шума нет. Наличие измерительного шума £ усложняет задачу. Эго связано с тем, что в.искомой зависимости ;с/|+| от хп «зашумлены» теперь не только значения «зависимой» переменной. *я+|, но и «независимой» - хп. При значительном уровне шума
смещение оценок, полученных с помощью-МНК (1.5), оказывается велико для сколь угодно длинного ряда [51, 52]. Уменьшить систематическую ошибку оценки в случае значительного измерительного шума отчасти удается при использовании метода полных наименьших квадратов [50], когда минимизируется сумма квадратов ортогональных расстояний от точек (rjni7jn+]) до графика f (хп ,с), рис. 1.1 ,в. Но это улучшение не принципиально.
Если записать функцию правдоподобия с учетом измерительного шума, то при нормальном шуме задача сводится к минимизации суммы квадратов отклонений реализации модели от наблюдаемого ряда (рис. 1.1,г):
S(c,x,)= ]T(?7„+i -P'\xuC)f -> min, (1.6)
/т=0
где /(,,) - п-я итерация отображения хя+1 = f(xn,c), /(0)(л',с) = дт, и в число оцениваемых величин включено начальное состояние х,.
Траектория хаотической системы очень чувствительна к начальным условиям и параметрам. Теоретически, дисперсия оценок (1.6) убывает очень быстро с ростом Ny иногда даже экспоненциально [49, 54]. Однако это достигается на практике только при условии, что удается находить глобальный минимум (1.6). При минимизации используют итерационные методы [65], согласно которым выбираются стартовые догадки для искомых величин и с
и по шагам происходит «спуск» в минимум (зачастую, в локальный). Обычно пробуют большое число М стартовых догадок и выбирают наиболее глубо-
34
кий из найденных минимумов. Догадки берутся из некоторой области, в которой предположительно должны находиться истинные значения оцениваемых величин. Часто используется параллелепипед с осями длины Ал по начальному условию и А, по параметрам ci (назовем Дх и А,- интервалами' стартовой неопределенности). Даже при умеренном N график функции S в случае хаотической системы становится сильно «изрезанным» (рис.1.2,а), так что найти глобальный минимум численными- методами [65] можно только при астрономически большом М-или очень малых Лх и А., что практически не реально. Для этого потребовалось бы очень «удачное» задание стартовых догадок для с и лг,. Об асимптотических свойствах оценок говорить трудно, т.к. в пределе N —> оо целевая функция становится негладкой [49].
В связи с описанными трудностями разрабатываются модификации ММП в приложении к задаче оценки параметров по хаотическому ряду. Так, часто применяемый кусочный подход [49] состоит в том, что фиксируются некоторые «практически разумные» значения Му АхУ А,-, а исходный ряд делится на сегменты длины L, для каждого из которых удается найти глобальный минимум (1.6) за «разумное» время. Итоговую оценку получают как эмпирическое среднее оценок, полученных по отдельным сегментам, обозначим ее dj- (от английского “forward” - прямо). Оценки по каждому из сегментов могут быть сильно смещенными, тогда итоговая оценка тоже смещена. Причем ее смещение не устраняется при увеличении N. Кроме того, дисперсия оценки убывает с ростом длины ряда всего лишь как N~l. Увеличивать длину сегмента L бесполезно, т.к. не удастся найти глобальный минимум (1.6). Это принципиальное ограничение описанного подхода.
1.1.2. Оценка с использованием обратного отображения
В диссертационной работе предложено использовать обратное отображение при расчете целевой функции, т.е. минимизировать величину:
Рис Л .2. Целевые функции 5(с,х,) и 3(с,хы) для квадратичного отображения (1.3) при N = 20, с0 = 1.85, х] = 0.3: а) при расчете целевой функции используются итерации отображения в прямом времени (1.6); б) в обратном времени (1.7). Графики соответствуют фиксированным (истинным) значениям х1 и Лдг.
36
S(c,xN) = “ /( "4xN,c))2 -> min, (1.7)
л=0
где /(w,) - л-ая итерация: отображения x,l+1 = f(xtt9c) в обратном времени. При'обращении времени ляпуновский показатель становится отрицательным -и чувствительность траектории к.«начальному» условию. xN пропадает. Поэтому следует ожидать, что целевая функция (4.7) будет иметь мало локальных минимумов (рис. 1.2,6). Минимизация (Г.7)> также должна проводиться численно при некоторых Му Ах, А) у но; как показывает численный эксперимент, предложенные оценки оказываются не чувствительны к этим величинам. Поэтому ниже используется М- 1, а Ах, А,- - такие же, как для дf.
Поскольку хаотический, режим может демонстрироваться только немо-нотонным одномерным - отображением, то возникает проблема неоднозначности обратного отображения-. При итерациях в обратном-времени нужно решать вопрос,.какой из корней уравнения х„+1 -/(х,с) = 0 брать в качестве хпу если хп+1 - значение, полученное на предыдущей итерации в обратном времени. В работе предложено брать тот корень, который ближе всего к соответствующему значению наблюдаемой rjn. Итоговая оценка обозначается сь от английского “backward”. Применимость метода ограничивается одномерными отображениями, т.к. для многомерной диссипативной системы траектория в обратном времени еще более чувствительна к параметрам.
Численный эксперимент. Погрешность оценки одного параметра д обычно определяется в математической статистике как средний квадрат ее
отклонения от истинного значения. Обозначим ее с2 =^(^-с0)2^, где угловые скобки означают математическое ожидание. В случае оценивания нескольких параметров (размерность вектора параметров Р > 1), обозначим е2
р
- погрешность оценки с{, as2— суммарную погрешность s2 = £s2.
/-!
37