Содержание к диссертации
Введение
1 Статистический подход к реконструкции динамических систем (ДС) 20
1.1 Введение 20
1.2 Частные случаи 24
1.3 Модификация Байесова подхода 29
1.4 Восстановление значений параметров известной ДС по зашумленному временному ряду (ВР): сравнение модифицированного и "классического" подходов 33
1.5 Классификация режимов поведения известной ДС по коротким зашум-ленным ВР 36
1.6 Заключение 41
2 Реконструкция слабонестационарных динамических систем по хаотическим временным рядам 43
2.1 Введение 43
2.2 Описание модели наблюдаемой динамики 44
2.3 Пример реконструкции: прогноз поведения ДС по хаотическому незашумленному ВР 48
2.4 Реконструкция ДС по зашумленному слабонестационарному ВР 65
2.5 Заключительные замечания 77
3 Восстановление профилей атмосферных характеристик по данным дистанционного зондирования 81
3.1 Введение 81
3.2 Статистический подход в приложении к задаче восстановления атмосферных профилей 85
3.3 Восстановление профиля атмосферного озона по данным миллиметровых измерений 90
3.4 Заключение 100
Заключение 102
- Частные случаи
- Классификация режимов поведения известной ДС по коротким зашум-ленным ВР
- Реконструкция ДС по зашумленному слабонестационарному ВР
- Восстановление профиля атмосферного озона по данным миллиметровых измерений
Введение к работе
Разработке методов решения обратных задач реконструкции динамических систем (ДС) на основе порожденных ими наблюдаемых процессов (временных рядов) посвящено в последние тридцать лет большое количество работ (см., например, [1-3] и цитируемую там литературу). Актуальность такого подхода к реконструкции связана с тем, что он не требует наличия полной и детальной априорной информации о процессах, протекающих в системе, т. к. не включает в себя процедуру построения моделей из первых принципов (уравнений движения среды или отдельных частиц, уравнений для силовых полей, переноса излучения, химической кинетики, тепло и массопереноса и пр.). Математическая модель исследуемой ДС при этом строится путем прямого анализа наблюдаемых данных, вообще говоря, без каких-либо допущений о природе изучаемого явления. Как правило, каждый из предложенных к настоящему времени методов такой реконструкции, включает в себя два основных шага: (1) реконструкция области фазового пространства, в которой система эволюционирует, и (2) построение модели, воспроизводящей поведение системы в этой области фазового пространства.
Фундаментальной основой методов реализации первого шага являются доказанные Такенсом теоремы [4], из которых следует, что непрерывная бесконечная и стационарная зависимость от времени всего лишь одной (наблюдаемой) динамической переменной системы является достаточной для восстановления топологической структуры области фазового пространства, соответствующей наблюдаемой эволюции системы. В реальности протяженность наблюдаемого временного ряда (ВР) всегда ограничена, что приводит к ошибке при реконструкции реализующегося в фазовом пространстве системы аттрактора, и, как следствие, неточно-
Введение
сти в определении его характеристик. Тем не менее, созданные к сего-. дняшнему дню алгоритмы позволяют достаточно аккуратно восстанавливать динамические свойства системы по конечному стационарному ВР. Кроме того, результаты многочисленных компьютерных экспериментов (начиная со статьи [5], инициировавшей упомянутую работу Такенса) с очень широким кругом моделей показали, что такая реконструкция может производится на основе дискретного ВР. Наконец, было показано [6], что теоремы Такенса могут быть обобщены на случай нестационарных ВР, т.е. нестационарность не налагает ограничений с точки зрения принципиальной возможности реконструировать наблюдаемый аттрактор.
Методы построения моделей, описывающих динамику исследуемой системы в восстановленном фазовом пространстве, можно разбить на две группы [1,2]. К первой относятся методы, направленные на реконструкцию локальной динамики системы в отдельно взятых элементарных ячейках фазового пространства. В рамках таких методов строятся модели, способные воспроизвести (предсказать) эволюцию малой окрестности выбранного вектора состояния системы в фазовом пространстве на времена порядка характерного масштаба изменения динамической переменной. Другими словами, для каждой интересующей точки фазового пространства строится локальный оператор эволюции (ОЭ), наилучшим образом предсказывающий дальнейшее поведение фазовых траекторий, попавших в малую окрестность данной точки, на определенный шаг по времени. В качестве функций, аппроксимирующих данный ОЭ, используются как полиномы различной степени [7] (в частности, полином нулевой степени - в этом случае предсказание заключается в простом усреднении по образам всех точек из выбранной окрестности), так и более сложные функции, например, системы радиальных базисных функций [8]. В случае хаотической динамики такие модели обеспечивают характерное время предсказания, обратно пропорциональное значению старшего ляпуновского показателя ВР, являющегося мерой разбегания изначально близких фазовых траекторий хаотического аттрактора [9]. Для успешной реконструкции поведения системы с помощью описанных
Введение 5
локальных моделей необходимо, прежде всего, чтобы окрестность каждой точки исследуемого аттрактора была хорошо посещаема восстановленной фазовой траекторией, т. е. протяженность наблюдаемого ВР (объем данных) должна быть достаточной для хорошего покрытия всего аттрактора. Другим требованием к наблюдаемым данным является стационарность исследуемого ВР (предполагающая постоянство управляющих параметров реконструируемой системы), поскольку описанные методы базируются на эргодической гипотезе. Главным недостатком локальных моделей является очень большое количество коэффициентов, требуемых для их описания, что очевидным образом снижает точность реконструируемых характеристик системы. Кроме того, высокая чувствительность к измерительному шуму сильно затрудняет их использование при работе с реальными данными.
Вторую группу образуют методы глобальной реконструкции динамических систем, направленные на построение модели оператора эволюции, действующего во всей области фазового пространства, соответствующей наблюдаемому ВР. Такие модели воспроизводят качественную структуру фазового пространства ДС, при этом, вообще говоря, можно восстанавливать различные характеристики исследуемого аттрактора, такие как вероятностное распределение, фрактальная размерность, ляпунов-ские показатели и др. Привлекательность такого подхода связана с тем, что весь набор имеющихся данных описывается непрерывной гладкой моделью с небольшим (по сравнению с локальными моделями) числом параметров. Более того, с помощью глобальных моделей можно отслеживать изменения управляющих параметров исходной системы, что находит применение в задачах реконструкции неавтономности системы (см. главу 2 и [10-14]), восстановления бифуркационных диаграмм [15-18], передачи информации [19-21] и т.д. Чаще всего глобальная модель строится в виде дискретного ОЭ, описывающего связь между состояниями системы, соответствующими соседним по времени пересечениям фазовой траектории с выбранной секущей аттрактора в восстановленном фазовом пространстве. При этом для аппроксимации ОЭ могут использовать-
Бведение О
ся различные функции, такие как, например, системы ортогональных (на инвариантной плотности аттрактора) полиномов [22], системы радиальных базисных функций [8,23,24], искусственные нейронные сети (ИНС) [12,13,25,26] и др. Кроме того, предложены методы построения потоковых глобальных моделей, представляющих собой системы дифференциальных уравнений заранее выбранного вида [27-32].
Особенностью глобальной реконструкции является то, что используемые модели являются параметризованными, т.е. представляют собой функции, зависящие от наборов параметров. При этом задача реконструкции ДС часто понимается, как поиск оптимального (обеспечивающего наилучшее соответствие модели и наблюдаемого процесса), с точки зрения выбранного критерия, вектора параметров модели. Математическим представлением такого критерия обычно является ценовая функция (ЦФ), определяющая меру близости наблюдаемых данных и модели, маленькие значения которой отвечают хорошему соответствию, а большие - плохому. Например, хорошо известной и наиболее широко используемой ЦФ является среднеквадратичная невязка, лежащая в основе метода наименьших квадратов (МНК) (см., например, [9,33,35]); в методе обобщенных наименьших квадратов (МОНК) используется ЦФ, наиболее эффективная в задачах аппроксимации данных, когда погрешность присутствует как в образах, так и в прообразах реконструируемого отображения [34,36,37]; при реконструкции параметров известной системы по ВР применяется метод множественной стрельбы [39,49], основанный на ЦФ, учитывающей "долгие" корреляции наблюдаемой динамической переменной. Кроме того, при реконструкции параметров простых систем, в основе которых лежат известные уравнения, может применяться метод статистических моментов [2,38], обоснованность которого очевидна для эргодических процессов. Однако, данный метод может применяться только в случае реконструкции систем простого априори известного вида, а эффективность его невысока [38].
Неотъемлемой чертой любого процесса реального наблюдения является наличие случайной составляющей в экспериментальных данных. Та-
Введение 7
кая случайность может быть следствием, например, шума измерительной аппаратуры, конечной точности измерений, дискретности ряда по времени и т.п. Другими словами, измеренные значения отличаются от тех, которые были сгенерированы изучаемой системой, и эти отличия являются случайными величинами. Таким образом, наблюдаемые в реальном эксперименте величины содержат стохастическую компоненту, и, следовательно, уместным является вероятностный подход к их описанию. Ясно, что подход к оценке ненаблюдаемых величин (параметров системы, скрытых шумом (латентных) динамических переменных и других искомых характеристик системы) по таким данным должен быть статистически обоснованным, т. е. максимально корректно учитывать вероятностное распределение наблюдаемых величин. Можно легко показать, что в противном случае, когда используемый подход основан на неверной информации о статистических свойствах переменных (в частности, использование перечисленных выше ценовых функций далеко не всегда оказывается статистически обоснованным, см. главы I и II), производимые оценки искомых характеристик системы оказываются неконтролируемым образом систематически смещенными (см., например, [40,41]). Это обстоятельство делает наиболее адекватным вероятностный (Байе-сов) подход к решению задачи реконструкции, в рамках которого реконструируемые ненаблюдаемые переменные предполагаются случайными величинами, для которых конструируется функция апостериорной плотности вероятности (АПВ). Кроме информации о шуме, такой подход в вероятностной форме задействует также всю имеющуюся априорную информацию о реконструируемой системе. В наиболее типичном случае, когда решаемая задача является некорректной, использование обоснованной априорной информации, обеспечивающее регуляризацию решения, является ключевым моментом при реконструкции.
Одной из основных целей данной работы является разработка и реализация эффективного метода глобальной реконструкции ДС, позволяющего извлекать из наблюдаемого ВР максимально полную информацию об исследуемой системе на основе Байесового подхода к решению обрат-
Введение 8
ных задач. Попытки построения методик глобальной реконструкции ДС, в которых в основе производимых оценок лежат статистические соображения, начали предприниматься сравнительно недавно. Первой работой, направленной на получение несмещенной оценки параметров известной ДС по зашумлешюму хаотическому ВР, является статья [40], в которой, во-первых, продемонстрировано растущее с увеличением уровня шума систематическое смещение оценок, полученных с помощью МНК и МОНК, и, во-вторых, предложена ценовая функция, учитывающая инвариантную плотность состояний в фазовом пространстве модели при различных значениях параметров1. На примере ВР, сгенерированного простейшим точечным отображением (logistic map), показано отсутствие систематической погрешности реконструкции параметра данной системы во всем представленном диапазоне уровня шума. Однако, как справедливо отмечено в последующих работах [38,41], предложенный метод построения ЦФ содержит ряд ошибок, связанных с неправильной статистической интерпретацией переменных, при этом несмещенность оценки достигается крайне неэффективным способом. В этих же работах (кроме того, см. [42-44]) также отмечается трудность практического использования апостериорной плотности вероятности ненаблюдаемых, построенной статистически корректным образом (в рамках "классического" Байесова подхода), в случае реконструкции динамической системы по зашумленному хаотическому ВР существенной протяженности. Основная проблема заключается здесь в практической невозможности учета "динамичности" системы в полной мере, что отражается в чрезвычайно сложной ("изрезанной") структуре соответствующей функции АПВ. Поскольку именно хаотические ВР, неся в себе информацию о значительной части фазового пространства системы, представляют особый интерес с точки зрения глобальной реконструкции ДС, преодоление данной трудности является очень существенным шагом, на что и были направлены дальнейшие усилия в разработке Байесовых методов. Так, работе [41]
1 Фактически, предложено исключить из ценовой функции МОНК латентные переменные, проинтегрировав ее на предельном множестве состояний модели.
Введение
предлагается подход, основанный на смягчении требований к динамичности исследуемой системы: при построении искомой АПВ предполагается наличие (кроме динамической) слабой стохастической связи между латентными переменными системы, т.е. по сути, в модель вводится динамический шум. Несмотря на то, что, как отмечено в работе [42], авторы фактически "неправомерно" подменили динамическую задачу стохастической, продемонстрировано, что такой подход позволяет получить несмещенные оценки на искомые характеристики системы. Можно, однако, показать [44], что при этом существенно снижается точность реконструкции из-за ослабления априорных требований к системе2. В другой работе [38] предложен способ построения "кусочной" функции правдоподобия, в рамках которого исследуемый ВР делится на сегменты одинаковой фиксированной протяженности, при этом налагается требование, чтобы модель "максимально хорошо" воспроизводила наблюдаемую траекторию на этих сегментах. Как параметры, так и латентные переменные считаются независимыми величинами в различных сегментах, их точечная оценка получается усреднением оптимальных (в смысле максимизации функции правдоподобия) значений по всем сегментам. Вопросы выбора оптимальной длины сегмента, а также способа корректной оценки точности реконструкции остаются при этом открытыми.
В диссертационной работе развивается способ модификации Байесо-ва подхода, направленный на построение функции АПВ, который корректно учитывает статистику шумов, и в то же самое время максимально возможным образом принимает во внимание динамичность исследуемой системы. Указанная выше трудность практического использования Байесова подхода для анализа хаотических рядов преодолевается при этом включением в процедуру реконструкции априорной информации
2Надо отметить, что использование стохастических моделей в рамках Байесовой реконструкции является довольно популярным, при этом основная проблема заключается в большом количестве латентных переменных, а, следовательно, и аргументов функции АПВ, что существенно затрудняет ее анализ. К настоящему времени предложены различные способы интегрирования АПВ по латентным переменным, например, с помощью метода Монте-Карло [45] (что является достаточно ресурсоемким) или приближенных рекуррентных оценок интегралов модифицированными фильтрами Кальмана [46,50]
Введение
о свойстве динамического хаоса, состоящем в потере информационной связи между отсчетами ряда с увеличением временного интервала между ними. Последнее позволяет получить факторизованное (модифицированное) апостериорное распределение ненаблюдаемых, пригодное для численного анализа. Предлагается также основанный на методе Монте-Карло метод статистического анализа построенной функции АПВ, необходимого для расчета доверительных интервалов искомых величин. Отметим, что предложенный метод позволяет, кроме всего прочего, производить оценку параметров распределения шума (например, дисперсию), вообще говоря, неизвестных априори и являющихся ненаблюдаемыми наряду с параметрами модели и латентными переменными.
Первая глава данной работы посвящена описанию и иллюстрированию на модельных примерах модифицированного Байесова подхода. В начале главы делается введение, содержащее общие положения Байесова подхода применительно к анализу ВР. Далее рассматриваются частные случаи, соответствующие ситуациям, когда измерительный шум мал и когда он существенен, вводится функция АПВ, соответствующая предлагаемому модифицированному Байесову подходу. Предлагается основанная на методе Metropolis-Hasting [51,52] реализация метода Монте-Карло, пригодная для численного анализа данной функции, необходимого для извлечения из него информации об искомых характеристиках системы. Описание соответствующего алгоритма вынесено в Приложение. Возможности предложенного подхода демонстрируются на двух примерах. В одном из них решается задача классификации неразличимых традиционными методами режимов поведения известной динамической системы по коротким (около 20 характерных периодов изменения динамической переменной системы) существенно зашумленным ВР. Другой пример посвящен реконструкции параметра логистического отображения по зашумленному ВР: проводится сравнительный анализ результатов реконструкции, полученных упомянутыми существующими и предлагаемыми в данной работе методами. В заключение главы обсуждаются возможные факторы, ограничивающие точность предложенного
Введение 11
метода реконструкции.
Вторая глава представленной работы посвящена разработке метода глобальной реконструкции неизвестной ДС в условиях, когда наблюдаемые данные представляют собой слабонестационарные хаотические процессы, т. е. соответствующие временные ряды имеют масштаб нестационарности много больший, чем характерные времена изменения динамических переменных. Такое поведение характерно, в частности, для различных систем, определяющих протекание важнейших процессов в атмосфере и гидросфере Земли (эволюцию озонного слоя [53], поведение концентраций химических составляющих атмосферы в приземном слое воздуха [54], крупномасштабных вариаций поверхностной температуры тропических вод Тихого океана (явление Эль-Ниньо) [55,56])3. Данное условие на масштаб нестационарности позволяет рассматривать систему, породившую данный ВР, как неавтономную с медленно изменяющимися во времени управляющими параметрами. Актуальность данной задачи обусловлена тем, что описанный тип неавтономности характерен для природных систем, т. к. они практически никогда не бывают замкнутыми, но зависят от изменяющихся с течением времени внешних условий, что может приводить к изменению во времени параметров, определяющих динамические свойства системы. Это означает, прежде всего, возможность смены типа поведения системы (бифуркации) в процессе ее эволюции, что влечет за собой существенные (иногда катастрофические) изменения количественных характеристик наблюдаемого процесса. Поэтому в описанной ситуации одной из главных задач, которую должна решать глобальная реконструкция, является задача прогноза качественного поведения ДС, в качестве первоочередной цели которого естественно определить предсказание возможных типов бифуркаций и моментов бифуркационных переходов. Разработанный в работе метод прогноза, основанный на Байесовой реконструкции, позволяет производить
3 Отметим, что круг реальных систем, демонстрирующих хаотическую динамику, уже сегодня простирается от перечисленных атмосферных и атмосферно-океанических процессов до различных систем в живых организмах [57,58] и тектонической активности [59] и имеет тенденцию к расширению по мере приложения современных методов анализа к новым базам данных различной природы.
Введение
вероятностное предсказание поведения исследуемой системы. В рамках этого метода задается явно зависящая от времени модель ОЭ системы, строится функция АПВ параметров этой модели, которая определяет ансамбль возможных сценариев поведения системы, по которому затем вычисляются вероятности искомых характеристик (такие как, например, тип динамического режима в интересующий момент времени или моменты времени бифуркаций) как на интервале времени наблюдения, так и в будущем (путем экстраполяции модели за пределы времени наблюдения).
В последние 10-15 лет в печати демонстрируется повышенный интерес к анализу нестационарных ВР, в частности, большое внимание уделяется проблемам установления факта нестационарности, а также определению характера и меры нестационарности процессов различной природы. Разработан ряд методов, основанных на анализе распределений временных интервалов между соседними векторами в фазовом пространстве [60-62], использовании меры взаимной предсказуемости динамики между различными участками ВР [63], исследовании зависимости от времени плотности вероятности процесса и его спектральной плотности [64] и др. В данной работе, по сути, предлагается метод восстановления нестационарности процесса на основе глобальной реконструкции неавтономной ДС, породившей данный процесс.
Во вводных разделах описываемой главы ставится задача реконструкции ДС по слабонестационарному ВР, а также излагаются некоторые особенности связанной с ней задачи прогноза качественного поведения ДС. Затем описывается предлагаемый способ аппроксимации неизвестного неавтономного ОЭ системы: в форме ИНС. Выбор такой параметризации обусловлен тем, что ИНС является универсальным аппроксима-тором [25,65], т. е. позволяет аппроксимировать любую функцию на выбранном интервале изменения аргумента с наперед заданной точностью. Для модели в виде ИНС предлагается метод использования априорной информации о системе путем внесения вероятностных ограничений на параметры модели, необходимость которых диктуется вырожденностью
Введение 13
пространства параметров модели. Далее отдельно рассматриваются два случая - когда реконструкция производится по ВР без шума измерений, и когда ВР содержит измерительный шум, существенно превосходящий ошибку, связанную с неточностью аппроксимации (дефект модели). Для обоих случаев строятся функции АПВ параметров модели, причем в последнем случае предлагается использовать АПВ, построенную в соответствии с модифицированным Байесовым подходом, описанным в первой главе.
Обстоятельством, усложняющим анализ вероятностного распределения в последнем случае, является наличие в нем степеней свободы, соответствующих латентным переменным системы, количество которых возрастает пропорционально протяженности исследуемого ВР. В рассматриваемой задаче прогноза поведения ДС используются реализации существенной протяженности (например, используемые в данной работе ВР включают в себя порядка 1000 отсчетов), поскольку требуется, чтобы процесс содержал в себе достаточную информацию о структуре фазового пространства восстанавливаемой системы. В результате размерность пространства ненаблюдаемых зачастую оказывается настолько большой, что численный анализ соответствующей функции АПВ становится практически невозможным. Для преодоления указанной трудности в работе предлагается метод приближенного интегрирования АПВ по латентным переменным, позволяющий получить АПВ только для параметров модели ОЭ.
Результаты прогноза с использованием построенных функций АПВ демонстрируются на примерах ВР, порожденных как точечным отображением (отображение Эно), так и системами обыкновенных дифференциальных уравнений (система Ресслера и система уравнений, описывающая химические процессы, протекающие в мезосфере Земли). В частности, в случае с зашумленными ВР приведенные результаты иллюстрируют высокую эффективность модифицированного Вайесова подхода с точки зрения как прогноза поведения системы, так и реконструкции наблюдаемой динамики. В заключении ко второй главе кратко излагаются
Введение 14
основные выводы из полученных результатов, обсуждаются ограничения и возможные усовершенствования разработанного метода.
Таким образом, в двух первых главах представлена основанная на статистическом подходе методика решения обратной задачи реконструкции системы по данным наблюдений, включающая в себя (1) статистически обоснованную постановку задачи, основанную на вероятностном представлении данных, (2) построение пригодной для анализа функции АПВ, включающей в себя информацию о шумовой составляющей сигнала, физически обоснованную априорную информацию о системе, а также параметризацию оператора эволюции системы, (3) методы извлечения из построенной функции информации об искомых характеристиках системы. Хотя предложенная методика изначально разрабатывалась для задачи реконструкции динамических систем по генерируемым ими наблюдаемым процессам, она оказывается эффективной при решении других (не "динамических") некорректных обратных задач. Это демонстрируется в третьей главе работы, где рассматривается задача восстановления высотных распределений (профилей) параметров атмосферы по данным наземного пассивного радиометрического зондирования, представляющим собой спектры излучения (поглощения) химических составляющих атмосферы.
Данная задача предполагает решение интегрального уравнения [66] (в общем случае нелинейного), при этом она является некорректной [67-70], поскольку неточность в исходных данных, обусловленная их дискретностью и существенной зашумленностью, приводит к бесконечному множеству возможных решений, в то время как существует точное и единственное решение невозмущенной задачи. Ключевым моментом при этом является включение в процедуру восстановления физически обоснованной априорной информации о профиле, ограничивающей класс возможных решений. Наиболее часто применяемые алгоритмы восстановления профилей атмосферных характеристик основаны на использовании метода наименьших квадратов, дополненного регуляризацией, содержащей необходимые априорные представления. Эти алгоритмы отличаются
Сведение
между собой преимущественно методами регуляризации, т. е. типами используемой априорной информации и способами ее включения. Например, в хорошо известном методе Тихонова [68] налагаются условия на гладкость решения, при этом жесткость ограничений задается условием равенства функции невязки и дисперсии шума измерений. Хотя сходимость решения к точному при стремлении уровня шума к нулю и доказана [68,69], такой определяемый шумом выбор априорных ограничений неизбежно приводит к систематической погрешности ("переглаженным" профилям), которая является неконтролируемой и может быть причиной принципиальных ошибок в случае не слишком малой зашумленности данных. Другой метод, предложенный и развитый в работах [71,72], использует для регуляризации априорную статистику, полученную из предыдущих измерений. Данный метод дает несмещенное решение только в том случае, когда в распоряжении исследователя имеется достаточно богатая статистика, учитывающая все физически возможные вариации профиля в данной географической точке, однако это условие выполняется далеко не всегда. Общим недостатком традиционных методов восстановления является использование кусочно-однородной (кусочно-линейной) аппроксимации профиля, что сильно ограничивает возможность включения в алгоритм априорной информации различного типа. Кроме того, большинство существующих методов направлены на решение линеаризованной задачи, в то время как связь между искомым профилем и наблюдаемым спектром часто описывается нелинейным интегральным уравнением. В ситуации, типичной для наземного зондирования, когда уровень шума достаточно высок, систематическая ошибка, вносимая линеаризацией, может быть существенна. Другим очень важным моментом при решении задачи восстановления является корректный расчет погрешности произведенных оценок. При использовании традиционных методов восстановления определение данной погрешности является отдельной весьма нетривиальной задачей даже в случае линейной связи измеряемой и восстанавливаемой характеристик (см., напри-мер, [73,74]).
Введение 16
В то время как существующие методы направлены на поиск единственного "оптимального" профиля, в данной работе предлагается метод [75,76], основанный на вероятностном представлении решения. В соответствии с описанным в первой главе подходом искомая величина интерпретируется как случайная, и для нее строится функция АПВ, включающая в себя как информацию о распределении шума измерительной аппаратуры, так и априорные ограничения, налагаемые на свойства восстанавливаемого высотного распределения. При этом исходные интегральные уравнения, лежащие в основе построения АПВ, могут быть нелинейными по отношению как к восстанавливаемой величине, так и к параметрам функции, аппроксимирующей профиль. Анализ функции АПВ методом Монте-Карло позволяет получить статистический ансамбль возможных профилей, по которому рассчитываются доверительные интервалы с заданным уровнем вероятности для искомой величины во всем зондируемом диапазоне высот. Таким образом, оценка погрешности восстановления является неотъемлемой частью предложенной процедуры и получается автоматически.
В первом параграфе третьей главы формулируется задача восстановления атмосферных профилей по данным дистанционного зондирования, а также подробно поясняются недостатки существующих методов. Далее описывается общий Байесов подход в приложении к рассматриваемой задаче, что в дальнейшем позволяет произвести прямое сравнение различных методов восстановления. Затем предлагается параметризация профиля функцией в виде ИНС. Показано, что использование аппроксимации в виде ИНС позволяет более эффективным образом (по сравнению с кусочно-линейной аппроксимацией) задействовать различную априорную информацию о профиле, такую как гладкость, диапазон вариаций концентрации озона, количество участков монотонности и т.д. В следующей части главы описывается конкретная задача - восстановление профиля озона по данным пассивного зондирования в миллиметровом диапазоне длин волн.Возможности метода демонстрируются сначала на модельном примере, имитирующем измерения радиационной температу-
Введение 17
ры атмосферы спектрометром ИПФ РАН в Арктике [77,78]. Показано, что предлагаемый метод позволяет восстановить профиль с достаточно резкими вариациями озона, моделирующими ситуацию с сезонным образованием арктической озонной дыры, в то время как метод, включающий в себя кусочно-линейную параметризацию и регуляризацию Тихонова, оказывается неэффективным. Далее представляются результаты восстановления профиля озона по реальным измерениям, проведенным в Апатитах зимой 2002-2003 гг. [78], которые затем сравниваются с данными модели GOME/ROSE [79]. В заключении главы излагаются выводы из полученных результатов.
В Заключении сформулированы основные результаты диссертации.
В Приложении изложен алгоритм реализации метода Metropolis-Hasting, разработанный для анализа используемых в работе апостериорных распределений.
Апробация работы. Изложенные в диссертации результаты докладывались на семинарах и конкурсах работ молодых ученых Института прикладной физики РАН, семинарах НИИ ПМК, ИФА РАН, кафедры математической статистики ВМК МГУ им. М.В. Ломоносова, на семинаре «Российская наука - XXI век» Минпромнауки РФ, в Лондонском Империал колледже (Великобритания), на конкурсах научных работ, на международных и общероссийских конференциях и совещаниях: 12-ой Генеральной ассамблее Международного союза геодезии и геофизики (1999 г., Бирмингем, Великобритания), Международной конференции, посвященной 100-летию со дня рождения А.А. Андронова "Прогресс в нелинейной науке" (2001 г., Нижний Новгород), Международной конференции Ха-ос'01 (2001 г., Саратов), 120-ой Фарадеевской дискуссии Королевского химического общества "Nonlinear Chemical Kinetics: Complex dynamics and Spatiotemporal Patterns" (Манчестер, Великобритания, 2001), XXXII Международном научно-методическом семинаре «Шумовые и деграда-ционные процессы в полупроводниковых приборах» (Москва, 2001), 4-ом и 5-ом Рабочих совещаниях программы "REACTOR" Европейского на-
Введение
учного фонда (2003 г., Будапешт, Венгрия; 2004г., Прага, Чехия), 35-ой Научной ассамблее COSPAR (2004г., Париж, Франция), Генеральной ассамблее Европейского союза наук о Земле (2005 г., Вена, Австрия), Международном симпозиуме "Topical Problems of Nonlinear Wave Physics" (2005 г., Санкт-Петербург - Нижний Новгород), 31-ом Международном симпозиуме по дистанционному зондированию окружающей среды (2005 г., Санкт-Петербург), Международной конференции "Динамические дни" (2006 г., Крит, Греция), 3-ой, 4-ой, 5-ой, 6-ой Нижегородской научной сессии молодых ученых (Н. Новгород, 1998, 1999, 2000, 2001 гг.), 2-й Всероссийской научной конференции студентов-радиофизиков (1998 г., С.-Петербург), 9-ой Всероссийской школе-семинаре "Волны -2004" (Москва), 6-ой, 7-ой, 8-ой, 9-ой и 10-ой Всероссийской школе-конференции молодых ученых "МАПАТЭ" (2000 г. и 2003г., Нижний Новгород; 2004г., Москва; 2005г., Борок; 2006г., Москва), Международном симпозиуме стран СНГ "Атмосферная радиация" (2002 и 2004 г., С.-Петербург), 11-ой, 12-ой и 13-ой Научной школе "Нелинейные Волны" (2002, 2004 и 2006 г.г., Нижний Новгород), 20-ой Всероссийской конференции по распространению радиоволн (2002 г., Н. Новгород), 15-ой научной сессии Совета по нелинейной динамике (2006г., Москва).
Основные результаты диссертации опубликованы в 5 статьях в реферируемых российских (Известия ВУЗов - Радиофизика) и международных (Faradey Discussions, Advances in Space Research, Physical Review E) научных журналах, 3 препринтах ИПФ РАН, 2 отчетах по программе фундаментальных исследований ОФН РАН, 6 сборниках трудов и 27 сборниках тезисов всероссийских и международных конференций.
Введение
Частные случаи
Рассмотрим два предельных случая, отражающих наиболее типичные ситуации, возникающие при реконструкции ДС по данным измерений. Первый случай относится к ситуации с "высококачественно" измеренными данными, т.е. когда шум, вносимый измерением, пренебрежимо мал. Во втором случае уровень шума измерений много больше погрешности модели (динамического шума) - такая ситуация имеет место при реконструкции неизвестной ДС по существенно зашумленным данным, когда используется достаточно качественная аппроксимация ОЭ системы. Нулевой шум измерений В работе [12] на нескольких примерах ВР, сгенерированных системами дифференциальных уравнений, продемонстрированы успешные результаты реконструкции этих систем в случае, когда шум измерений отсутствовал. Для поиска параметров /і ОЭ в этой работе использовался (нелинейный) метод наименьших квадратов (МНК): оптимальным счи- /. Статистический ПОДХОД К реконструкции динамических систем (ДС) тался вектор /І, соответствующий минимуму функции невязки Нетрудно показать, что в этом случае, в предположении, что случайная величина г] - дельта-коррелированная гауссова с нулевым средним (wrj(r],oz) ос -щг е 2аг ), использование функции (1.7) является ста-тистически обоснованным, что и предопределило успех такой реконструкции. Действительно, представив плотность вероятности шума измерений в виде нормального закона с дисперсией о\ и считая отсчеты некоррелированными, введем априорное ограничение на сг, отражающее информацию об отсутствии шума измерений (а = 0): зададим функцию априорной плотности вероятности для а\ в виде дельта-функции Ppr(crj) = 8(&Т)- Для априорной вероятности величин а , /і и и выберем неинформативное распределение Pw{a ,ii, u) = const. Функция (1.5) при этом принимает вид: Как видно, в показателе экспоненты получившейся функции стоит невязка МНК (1.7), при этом максимум плотности вероятности (1.9) соответствует минимуму (1.7). pps(fl, (7ЧХ, t) = ffppa(fJ., U, г, 0-2х, t)rf(jdu ОС Существенный шум измерений В случае, когда исследуемый сигнал содержит существенный измерительный шум , соображения, лежащие в основе построения функции (1.9), несправедливы. Ложные результаты реконструкции с использованием распределения (1.9) в этом случае будут продемонстрированы в разделе 2.4 на примерах модельных ВР. Получим статистически корректное апостериорное распределение вероятности ненаблюдаемых для случая с ненулевым шумом измерений.
Ограничим наше рассмотрение ситуацией, когда модель ОЭ является достаточно "хорошей", т. е. среднеквадратичный разброс ее случайных различий с ОЭ наблюдаемой системы можно считать много меньшим aj. Другими словами, будем далее использовать Аналогично выкладкам при выводе функции (1.9), зададим нормальный закон распределения для г), априорное распределение для дисперсии а в виде дельта-функции, а также неинформативное априорное распределение для и и aj. К вопросу априорных ограничений на параметры модели [і вернемся в разделе 2.2, пока же лишь обозначим их в виде функции априорной плотности вероятности Ррг{у). Затем проинтегрируем распределение (1.5) с описанными априорными сомножителями по ai и латентным переменным: К сожалению, применение описанного "классического" Байесова подхода к реконструкции ДС по зашумленным хаотическим данным встречает, как указано в [41,42,44] фундаментальные трудности. Увеличение протяженности ВР, желательное для уменьшения эффективного уровня шума измерений, делает невозможным использование АПВ в виде (1.11) для вычислений вероятностей требуемых характеристик. Очевидно, что выражение (1.11) в случае хаотического достаточно длинного ВР будет слишком сложным как для Монте-Карло анализа, так и для поиска наиболее вероятных значений параметров и начальных условий. Это обусловлено тем, что из-за фрактальной природы странного аттрактора увеличение Т приводит к чрезвычайному усложнению формы областей значений латентных переменных и параметров модели, обеспечивающих существование в фазовом пространстве модели траектории, лежащей в шумовой окрестности фазовой траектории, реконструированной по исходному ВР. Это означает, что зависимость функции (1.11) от своих аргументов становится мультимодальной (сильно изрезанной).
Описанная неприменимость "классического" Байесова подхода к реконструкции хаотических ВР продемонстрирована в работе [42] на примере логистического отображения, ОЭ которого является дискретным отображением первого порядка ип+\ = 1 — аи\. Данная система может демонстрировать в зависимости от значения управляющего параметра а как регулярное (периодическое), так и хаотическое поведение. Переход к динамическому хаосу при изменении а происходит через каскад удвоений периода. В [42] была рассмотрена реконструкция значения параметра а по хаотическому ВР, содержащему аддитивный шум измерений : ип = 1 - аи2п_х Рассматривался случай, когда 5-коррелированный шум имеет гауссово распределение с известной дисперсией сф /. Статистический подход к реконструкции динамических систем (ДС) Было показано, что в классической байесовой формулировке (1.11) характерный масштаб изрезанности условной плотности вероятности ррг(щ, а\х) становится меньше машинной точности даже в случае умеренного шума j = 0.13 и не очень протяженного ВР (Т=70). Для поиска правильного значения параметра среднеквадратичный разброс aUo, характеризующий ширину априорного распределения р(щ,а), должен быть меньше Ю-17. Другими словами, классический Байесов подход требует недостижимо точной информации о начальном состоянии системы. В работе [41] Meyer and Christensen предложили модификацию классической Байесовой формулировки, направленную на преодоление указанной проблемы (несовместимость статистического подхода с динамической природой исследуемой системы). Критикуя работу [40], в которой выбор ценовой функции не являлся статистически обоснованным, Meyer and Christensen предложили в рамках Байесова подхода использовать допущение о присутствии в системе малого динамического шума. При этом второе уравнение системы (1.10) принимает вид u +i = f(ut,/z) + щ, где щ- гауссова случайная величина. Формально, как отмечено в [42], такая постановка является некорректной: систему (1.10), про которую мы знаем, что она динамическая, мы подменяем системой стохастической. Однако с точки зрения байесова подхода такое допущение является вполне корректным и означает просто ослабление априорных требований к модели: детерминированная связь (1.3) латентных переменных априори заменяется более "мягкой" - вероятностной. В [41,46] демонстрируется, что полученная в результате плотность вероятности (которая теперь, кроме щ, включает в себя также и остальные латентные переменные системы) позволяет осуществлять статистический анализ постериорного распределения, например, МСМС методом, используя при этом неинформативные априорные распределения для ненаблюдаемых. Модификация Байесова подхода, предлагаемая в данной работе, основана на использовании другой, в некотором смысле противополож- 3Для ВР, рассмотренного в [42] (ВР, сгенерированный системой (1.12) при а = 1.85; UQ = 0.3), j = 0.1 соответствует отношению сигнал-шум примерно 6:1. I. Статистический ПОДХОД К реконструкции динамических систем (ДС) ной [41] идее: возможно полнее использовать априорную информацию о динамической природе системы.
Классификация режимов поведения известной ДС по коротким зашум-ленным ВР
Полученные в предыдущем разделе результаты позволяют поставить важную с практической точки зрения задачу классификации режимов поведения ДС. В качестве примера рассмотрим очень короткие (Т = 20) ВР, сгенерированные более сложной ДС - отображением Эно (Непоп). Оператор эволюции данной ДС представляет собой функцию последова-ния 2-го порядка: В зависимости от значений параметров а и Ьсистема (1.18) демонстрирует различные регулярные и хаотический режимы. Области параметров, отвечающие различным режимам поведения, выделены оттенками серого на бифуркационной диаграмме (рис. 1.4). В качестве "измеряемой" рассматривалась величина х, представлявшая собой аддитивно зашумленную (латентную) переменную u: xn = un + п. Шум измерений был -коррелирован и имел гауссово распределение с дисперсией Т. Анализировались ВР, сгенерированные при различных значениях параметров а и Ь. В отсутствие шума ВР, отвечающие различным периодическим и хаотическому режимам, легко различимы. На рисунках 1.5(а,б) показаны для примера ВР, порожденные отображением (1.18), соответствующие двухпериодическому (а = 0.85, Ь = —0.2) и хаотическому (а = 1.15, Ь — —0.2) режимам. При добавлении даже сравнительно небольшого шума (см. рисунки 1.5(в,г), а = 0.3) различить динамические режимы становится невозможным. Покажем, как модифицированный Байесов подход решает эту проблему. Решение задачи классификации режимов поведения ДС, породившей исходный зашумленный ВР, основано на исследовании статистического ансамбля иезашумленных ВР, отвечающего статистическому ансамблю параметров, распределенных в соответствие с апостериорной плотностью вероятности (1.15) (последний ансамбль получается применением МСМС алгоритма к апостериорному распределению (1.15)). Подсчиты-вается, сколько элементов ансамбля незашумленных ВР соответствуют тому или иному динамическому режиму, затем вычисляется вероятность каждого из зарегистрированных режимов. На рис. 1.6 показана зависимость вероятности "правильного" режима от уровня шума для хаотического и одного из периодических зашумленных ВР, сгенерированных отображением (1.18). Видно, что вплоть до уровня шума 25% в случае хаотического ВР и 50% в случае периодического, модифицированный алгоритм выделяет правильный режим с вероятностью не менее 0,5. Отметим, что близость к единице вероятности наиболее вероятного режима является удобной количественной характеристикой качества классификации. Обсудим зависимость данной характеристики от длины сегмента w.
Из обсуждавшегося выше ясно, что в случае хаотического режима с увеличением w эта вероятность будет, вообще говоря, приближаться к единице: из-за уменьшения с ростом w ширины апостериорной плотности распределения параметров вероятности ошибочных режимов будут почти всегда уменьшаться. Нетрудно понять, что этот вывод справедлив и в отношении зависимости от w вероятности "правильного" регулярного режима, хотя зависимость ширины распределения параметров от w становится несколько иной. Именно, для регулярных режимов данная зависимость не является экспоненциально убывающей, как в случае хаотического режима (см. (1.18)). Этот факт объясняется тем, что на аттракторе, соответствующем регулярному режиму динамической системы, отсутствуют экспоненциально быстрое разбегание изначально близких фазовых траекторий, и, как следствие, экспоненциальный рост чувствительности текущего состояния системы к начальным условиям с увеличением времени наблюдения5. Это приводит к более медленному, по сравнению с экспоненциальным, закону уменьшения ширины распределения (1.15) с увеличением w. Тем не менее, уменьшение ширины p(/i\x) с ростом w обеспечивает и для регулярных режимов поведения увеличение вероятности правильного режима при увеличении длины сегмента. Результаты восстановления системы (1.18), приведенные на рисунках 1.4 и 1.6, подтверждают сказанное. На рисунке 1.4 выделены доверительные области восстановления параметров системы (1.18) по соответствующим (зашумленным) ВР. Видно, что при w 4 качество классификации улучшается при увеличении длины сегмента как для хаотического, так и для регулярных режимов поведения. На рис. 1.7 показана экспоненциально уменьшающаяся зависимость от длины сегмента размера доверительной области, соответствующей хаотическому режиму, а также аналогичная, но немонотонная зависимость для регулярного (двухпериодического) режима. 5Другими словами, слагаемые в (1.15), содержащие Р достаточно высоких порядков, слабо зависят от конкретных значений латентных переменных. /. Статистический подход к реконструкции динамических систем (ДС) 40 В главе обсуждался статистический байесов подход к реконструкции динамических систем (ДС) по сгенерированным ими временным рядам (ВР). Именно такой подход наиболее, на наш взгляд, адекватен задаче реконструкции "реальных" ДС, которые, во-первых, подвержены нерегулярным внешним воздействиям в процессе генерации сигнала, и, во-вторых, сам сигнал испытывает случайные возмущения как при распространении, так и во время регистрации. Статистический подход имеет целью не только реконструкцию оператора эволюция ДС, но и построение плотности вероятности параметров этого оператора, рассматриваемых как случайные величины.
Одна из двух трудностей в реализации такого подхода возникает при реконструкции по зашумленным хаотическим рядам: экспоненциально быстрое разбегание фазовых траекторий на хаотическом аттракторе имеет следствием нарастающую с увеличением продолжительности ВР изрезанность плотности вероятности параметров. В результате уже для сравнительно коротких ВР построение плотности вероятности параметров становится невозможным. Нами предложена модификация классического байесова подхода, позволяющая преодолеть указанную трудность. Точнее, показано, что учет априорной информации о хаотическом характере генерируемого системой процесса позволяет получить статистически корректное выражение для апостериорного распределения, пригодное для практического использования при анализе ВР, содержащих измерительный шум (см. разделы 1.4, 1.5, 2.4 а также [43,44]). Длина сегмента w в выражениях (1.14) и (1.15), отражающая эту информацию, определяет максимальное число динамических связей между последовательными отсчетами исследуемого ВР, которое может быть учтено при реконструкции системы: чем больше это число, тем точнее реконструкция. Эффективность модифицированного подхода продемонстрирована на примере решения задач отыскания неизвестных значений параметров известной ДС и классификации режимов поведения ДС (как регулярных, Нередко (например, при исследовании "природных"динамических систем) наблюдаемый динамический процесс не является стационарным. Как отмечалось во Введении, были предложены несколько способов выявления нестационарности ВР (см., например, [60-64,83,84]). Еще один способ выявления нестационарности является частью излагаемого ниже алгоритма. Утверждение о нестационарности наблюдаемого ВР означает наличие у соответствующего процесса по меньшей мере двух сильно различающихся временных масштабов. В данной работе мы рассмотрим ситуацию, когда продолжительность наблюдаемого ВР существенно превышает меньший из этих двух масштабов, будучи при этом меньше большего из них. Такая ситуация характерна для многих долгоживущих природных систем, медленные тренды параметров которых обусловлены плавными изменениями внешних условий. В описанной ситуации естественными с физической точки зрения являются следующие две гипотезы.
Реконструкция ДС по зашумленному слабонестационарному ВР
Рассмотрим теперь ситуацию, в которой шум измерений , входящий в слабонестационарный хаотический ряд х, существенен настолько, что деффект модели г/ можно считать пренебрежимо малым, т.е. а » av. Малость деффекта модели всегда может быть обеспечена, как упоминалось выше, повышением качества аппроксимации ОЭ системы (в случае предлагаемой в работе параметризации - увеличением количества нейронов ИНС), поэтому такая ситуация типична при работе с зашумлен-ными данными. В разделе 1.3 в общем виде выведены соответствующие описанному случаю выражения для функции АПВ (1.14) и (1.15). Как отмечалось в 1.3, несимметричность функции (1.14) проявляется в ее вариациях при изменениях начального отсчета ВР, что особенно существенно при анализе коротких ВР. Там же приводится симметризованное выражение (1.15) для апостериорной плотности вероятности, устраняющее данный недостаток. Ясно, однако, что описанный эффект исчезающе мал, когда мы имеем дело с ВР большой протяженности. Использование выражения (1.15) в этом случае приводит только к увеличению числа латентных переменных в w-І раз, что нежелательно с точки зрения численного анализа функции АПВ. В данной главе мы рассматриваем задачу реконструкции неизвестного ОЭ, при этом наблюдаемый ВР должен нести достаточную информацию для детального восстановления фазового пространства системы, т.е. быть достаточно протяженным (в рассматриваемых ниже примерах используются ВР, содержащие порядка 1000 характерных периодов изменения динамических переменных). Руководствуясь изложенным выше соображением, будем в дальнейшем использовать для плотности вероятности более простое выражение (1.14). Итак, полученное распределение ненаблюдаемых имеет большое количество степеней свободы, определяемое суммой числа параметров модели, числа латентных переменных, и числа параметров распределения шума измерений. Кроме того, зависимость функции (1.14) от ее многомерного аргумента является мультимодалыюй вследствие наличия в ней рекуррентных последовательностей {fJ(ufc, //)} , причем количество локальных "всплесков" тем больше, чем большее значение длины сегмента w мы используем (напомним, что увеличение w в допустимых пределах повышает точность реконструкции). Таким образом, функция АПВ в виде (1.14) в рассматриваемом случае протяженного ВР оказывается чрезмерно сложна для численного анализа.
В частности, прямое применение метода МСМС к функции (1.14), успешно используемое при анализе коротких ВР (см. главу I), оказывается совершенно неэффективно. В следующем параграфе предлагается метод снижения размерности ненаблюдаемых, делающий возможным статистический анализ АПВ параметров модели. Далее, на примерах реконструкции ДС по слабонестационарным ВР будет продемонстрирована полезность предложенного метода для прогноза будущего поведения ДС, а также будет показана неэффективность более простой функции (2.5), пригодной в случае неза-шумленного ВР (см. 2.3), при работе с "испорченными" шумом ВР. II. Реконструкция слабонестационарных динамических систем по хаотическим временным рядам 00 Метод анализа функции апостериорной плотности вероятности (АПВ) В главе I (и в работе [44]) совместная апостериорная функция плотности вероятности латентных переменных и и параметров \i численно исследовалась методом Монте-Карло: с помощью метода Metropolis-Hasting (МН) [51] генерировалась марковская последовательность наборов переменных (u, fi)i, статистическое распределение которой сходилось к заданной функции вида (1.15). Легко видеть, что ансамбль {/ІІ}І=1, извлеченный из стационарной части этого процесса (установившейся в статистическом смысле последовательности), соответствовует интегральному апостериорному распределению параметров (1-6). В нашем случае (протяженных ВР) такой метод интегрирования функции (1.14) по латентным переменным практически не реализуем, т. к. размерность искомого распределения чрезвычайно велика из-за большого объема используемых для реконструкции данных и, как следствие, большого числа латентных переменных, в результате чего сходимость метода МН будет очень медленной (см. работу [93]). В данной главе предлагается другой метод исключения латентных переменных, основанный на приближенном интегрировании (1.14) по и. Пользуясь тем, что к-й сегмент функции (1.14) зависит только от одной векторной латентной переменной щ, представим интеграл от (1.14) по латентным переменным в виде произведения интегралов по переменным, отвечающим различным сегментам: рря{ц,at\x,t) ос Pprb JnX стоящих в показателях экспонент (2.10), от и имеет резко выраженные максимумы (как указано главе I, максимум тем резче, чем больше длина сегмента w), поэтому для интегрирования (2.10) применим метод пере- вала [95]. В каждом отдельно взятом сегменте разложим х2(/Л иг) в РЯД Тейлора в окрестности минимума и ограничимся квадратичными членами разложения: где uf и Cf - зависящие от \х значение латентной переменной, соответствующее максимуму показателя экспоненты, и матрица вторых производных х2(д?иг) п0 латентным переменным, соответственно. Таким образом, каждое подынтегральное выражение в функции (2.10) аппроксимируется гауссовым законом и затем аналитически интегрируется по и: Далее, полученная приближенная функция апостериорного распределения для jf и // (2.12) используется для генерации ансамбля параметров: на каждом шаге итерационной процедуры отыскиваются значения uk(w+i)+v Расчитываются матрицы С[, . . и методом МН выбрасывается следующий элемент искомого ансамбля {( 7f,/i)i}. Как и ранее, предложенная процедура позволяет не только извлекать информацию о распределении параметров ОЭ (1.6), но и оценивать дисперсию измерительного шума j в случае, когда данная величина неизвестна априори.
Примеры реконструкции: прогноз поведения ДС по хаотическим зашумленным ВР Проиллюстрируем предложенный метод на примерах ВР, сгенерированных динамическими системами с медленно меняющимися управляющими параметрами. Для этого мы будем использовать отображение Эно, задаваемое двухмерной функцией последования і = 0.2 + z[x - с) Используемый нами ВР, порожденный системой (2.13), представляет собой последовательность значений двухмерного вектора {х }А.=1, N = 3000. В процессе генерации этого ряда управляющий параметр а системы (2.13) линейно менялся со временем в диапазоне от 1,39 до 1,23. Нами исследовались два случая: в первом шум измерений отсутствовал, во втором к ВР добавлялся аддитивный белый гауссов шум с нулевым средним и среднеквадратичным отклонением а = 0.1. На рис. 2.7а горизонтальной скобкой выделен незашумленный ряд первой компоненты {a4}fc=1 наблюдаемого вектора, имеющий протяженность = 3000 (в единицах времени системы (2.13)). Вне выделенного интервала показано продолжение этого ряда "в будущее". Множества, образованные исследуемыми ВР на плоскости динамических переменных системы (2.13), показаны на рис. 2.9. Скалярный ВР y(tk), сгенерированный системой (2.14), имеет протяженность Т = 7500 (в единицах времени системы (2.14)), причем на этом временном отрезке управляющий параметр с линейно меняется от 5.13 до 4.41. Мы использовали для реконструкции два таких ВР: первый из них "регистрировался" с малым дискретом (шагом по времени At = 0.025), что соответствует случаю пренебрежимо малого шума измерений (далее этот ряд будем называть незашумленным); второй ряд содержал существенную шумовую составляющую, возникающую вследствие существенно дискрета At = 0.1. Для извлечения из такого ряда выборки данных, пригодной для построения дискретной модели вида (2.1), мы поступали следующим образом. Методом задержек (для данного ряда выбиралась задержка г = 5) по скалярному ВР восстанавливалась фазовая траектория {x(fc)}f 0, х Є Ud в пространстве размерности d = 3, затем фиксировались пересечения этой траектории с секущей, за-
Восстановление профиля атмосферного озона по данным миллиметровых измерений
Рассмотрим задачу восстановления вертикального профиля озона по данным пассивного наземного зондирования атмосферы в миллиметровом диапазоне длин волн. Экспериментальные данные для восстановления представляют собой спектр радиационной температуры атмосферы в полосе частот, соответствующей одной из линий собственного излучения озона. Данные проходят сначала предварительную обработку, в процессе которой спектр излучения озона отделяется от фона, обусловленного излучением других составляющих атмосферы (молекулярного кислорода и водяного пара). Описание прямой задачи Прямая задача (вычисление спектра излучения озона по известному профилю концентрации) решается следующим интегральным уравнением, моделирующим столкновительный механизм уширения линии [103,1041і: Здесь N(z) - зависимость концентрации озона от высоты, T(z) - температурный профиль, K(z, /) - эффективное сечение поглощения молекулы озона. Свойства функций K(z,fi), а также расстановка и количество частотных каналов /;, определяют чувствительность уравнения (а, значит, и процедуры восстановления) к различным вариациям профиля. На рис. 3.1 показаны высотные зависимости функции К для разных отстроек по частоте от центра линии, соответствующей вращательному переходу озона на частоте НО ГГц. По расположению максимумов можно увидеть, что решение уравнения (3.10) чувствительно к возмущениям выражение (ЗЛО) справедливо для высот ниже 70 км, на больших высотах столкновительный механизм сменяется доплеровским. III. Восстановление профилей атмосферных характеристик по данным дистанционного зондирования 90 Рис. 3.1. Эффективное сечение рассеяния молекулы озона в зависимости от высоты для различных частотных каналов спектрометра. профиля на высотах от 20 до 50 км. Таким образом, в этом диапазоне высот можно ожидать успешного восстановления. Модельная задача Для демонстрации преимущества использования нейронной сети над кусочно-линейной аппроксимацией решалась модельная задача восстановления. В рамках этой задачи в качестве исходных данных использовался спектр оптической толщины озона т: о Отметим, что данная задача является линейной относительно концентрации озона N2.
Для иллюстрации возможностей изложенного подхода был смоделирован спектр оптической толщины атмосферного озона в полосе частот, соответствующих линии излучения озона НО 2В некоторых работах (см., например, [103]) уравнение (3.11) используется для восстановления профиля по реально измеренным данным. При этом делается предположение об изотермичности атмосферы, что позволяет грубо расчитать оптическую толщину по яркостной температуре, используя алгебраические соотношения, полученные из уравнения (3.10). III. Восстановление профилей атмосферных характеристик по данным дистанционного зондирования 91 цедура моделирования предполагала весьма скромные, по современным представлениям, характеристики как чувствительности озонометра, так и ширины полосы анализа спектра, и заключалась в следующем. Сначала конструировалось модельное распределение концентрации озона по высоте. Затем по такому профилю в соответствии с соотношением (3.11) рассчитывался спектр оптической толщины озона. Количество и расстановка частотных каналов, а также уровень шума в данном спектре соответствовали параметрам измерений, проводимых в 1998 - 2003 годах в Апатитах (Murmansk, 67N,33E) спектрометром ИПФ РАН [77,78]. Далее к полученному спектру прибавлялся белый гауссов шум, после чего зашумленная линия использовалась для восстановления профиля концентрации озона. Среднеквадратичный разброс шумовой добавки соответствовал типичным условиям наблюдений при усреднении по 10-15 отсчетам и составлял примерно 50% от уровня сигнала в крайних (отстоящих на 110-130 МГц от центра линии) спектральных каналах. На рисунке 3.2а в качестве примера приведен один из модельных профилей озона, включающий как участки плавной зависимости концентрации озона от высоты, так и резкий "провал" в диапазоне высот от примерно 15 до 30 км. Такое распределение озона по высоте позволяет, во-первых, исследовать возможности предложенного подхода. Во-вторых, данное распределение моделирует реальную ситуацию, которая возникает в процессе формирования озоновой дыры в зимне-весенней полярной стратосфере [105]. Затем по такому профилю рассчитывался спектр оптической толщины озона. Рассчитанная в соответствии с (3.11) и зашумленная, как описано выше, линия показана на рис. 3.26.
Решение модельной задачи: сравнение кусочно-однородной и нейронно-сетевой параметризаций профиля Для восстановления модельного профиля озона нами использовались функционалы правдоподобия, соответствующие как методу Тихонова (использующему кусочно-однородную аппроксимацию профиля), так и Рис. 3.2. а) модельный профиль концентрации озона; б) спектр оптической толщины озона, рассчитанный по модельному профилю с помощью 3.11. Пунктирная линия - не зашумленный спектр, сплошная линия - спектр, содержащий шум со среднеквдратичным отклонением равным Ю-3 нп. предлагаемому в данной работе методу (выражения (З.б), (3.7) в предыдущем параграфе). При использовании кусочно-однородной параметризации искомый вектор параметров fi представляет собой вектор N значений концентрации озона на высотах (} гДе Н - количество однородных слоев, на которые условно разбивается атмосфера. С учетом уравнения (3.11), функционал правдоподобия (3.6) может быть записан в матрично-векторной форме: где К п = K(fi,zn), А - матрица, производящая из вектора N вектор производных профиля по высоте. Нетрудно видеть, что распределение вероятности величин N, соответствующее функции (3.12), является гауссовым, а его среднее и ковариационная матрица DN выражаются через соотношения: где а — т1—. Некорректность исходной нерегуляризовашюй задачи (случай а = 0) отражается в плохой обусловленности матрицы КТК, вследствие чего ошибка DN будет принимать бесконечно большие зна- чения при сколь угодно малом уровне шума сг. С ростом параметра регуляризации а, обусловленность задачи будет улучшаться (матрица КТК + аАтА перестает быть сингулярной), а случайная ошибка восстановления будет падать, при этом, однако, систематическое смещение решения будет расти. Таким образом, при заданном априорном распределении, все параметры нормального распределения (3.13) полностью определены, что позволяет вычислить доверительные интервалы величин N. Использование предлагаемой нами параметризации в виде ИНН (3.8) является более сложным алгоритмически, т.к. зависимость функционала правдоподобия (3.4) от параметров аппроксимации профиля JJ, существенно нелинейна, а значит и распределение Pps(ji\T) описывается законом, сильно отличающимся от нормального. Для анализа данного распределения нами применялась предложенная модификация метода Metropolis-Hasting, подробно описанная в приложении. С помощью этого метода генерировался статистический ансамбль параметров профиля, вероятностное распределение которого задавалось функцией (3.4).
