автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Весовые алгоритмы статистического моделирования переноса поляризованного излучения и решение задачи восстановления индикатрисы рассеяния

О0340и^о

На правах рукописи

Чимаева Анна Сергеевна

ВЕСОВЫЕ АЛГОРИТМЫ СТАТИСТИЧЕСКОГО МОДЕЛИРОВАНИЯ ПЕРЕНОСА ПОЛЯРИЗОВАННОГО ИЗЛУЧЕНИЯ И РЕШЕНИЕ ЗАДАЧИ ВОССТАНОВЛЕНИЯ ИНДИКАТРИСЫ РАССЕЯНИЯ

05.13.18 - математическое моделирование, численные методы и комплексы программ

Автореферат

диссертации на соискание ученой степени кандидата физико-математических наук

Новосибирск - 2009

2 2 О К 7 ?

003480281

Работа выполнена в Институте вычислительной математики и математической геофизики СО РАН.

Научный руководитель: кандидат физико-математических наук,

старший научный сотрудник, Ухинов Сергей Анатольевич

Официальные оппоненты: доктор физико-математических наук,

профессор,

Иванов Михаил Самуилович

кандидат физико-математических наук, доцент,

Соболева Ольга Николаевна Ведущая организация: Институт теплофизики СО РАН

Защита состоится 10 ноября 2009 года в 16 часов на заседании диссертационного совета Д 003.061.02 при Институте вычислительной математики и математической геофизики СО РАН (630090, Новосибирск, проспект Лаврентьева, 6).

С диссертацией можно ознакомиться в библиотеке Института вычислительной математики и математической геофизики СО РАН.

Автореферат разослан 8 октября 2009 года.

Ученый секретарь диссертационного совета д.ф.-м.н.

С.Б.Сорокин

Общая характеристика работы

Актуальность темы.

Имеется целый ряд физических проблем, требующих достаточно точного решения задач теории переноса излучения в атмосфере с учетом поляризации. Это, прежде всего, задачи интерпретации оптических наблюдений.

Вариации яркости неба главным образом определяются изменчивостью аэрозольной компоненты атмосферы, которая обладает сильной рассеивающей способностью. Теоретическое исследование проблем переноса излучения в такой среде невозможно без знания ее оптических параметров. Одним из них является аэрозольная индикатриса рассеяния.

Диссертационная работа посвящена исследованию дисперсии стандартной векторной оценки метода Монте-Карло, применяющейся при решении задач переноса поляризованного излучения, а также разработке и обоснованию методов решения задачи восстановления индикатрисы аэрозольного рассеяния атмосферы по наземным наблюдениям "яркости неба" в альмукантарате Солнца с учетом поляризации излучения.

Математическое описание распространения поляризованного излучения предоставляет нам удобный инструмент для исследования этого процесса. Базовыми при этом являются вектор-функции Сгок-са, характеризующие свойства излучения в каждой конкретной точке фазового пространства, и интегро-дифференциальное уравнение переноса описывающее процесс переноса.

Интенсивность и состояние поляризации излучения определяются четырехкомпонентной вектор-функцией Стокса, компоненты которой определяют в совокупности интенсивность, степень поляризации, плоскость поляризации и степень эллиптичности излучения. Процесс переноса излучения в этом случае можно описать интегральным уравнением второго рода, оператор которого, в силу физических особенностей задачи, оставляет инвариантным множество вектор-функций Стокса.

Рассматриваемая математическая модель позволяет ставить достаточно большое множество практически интересных задач, для решения которых может быть эффективно применен метод Монте-Карло. Традиционный способ его использования заключается в следующем. Рассматривается некоторый линейный функционал J от

решения уравнения переноса, для него строится стандартная весовая оценка статистического моделирования математическое ожидание которой и дает нам искомое значение функционала.

Конкретный вид функционала ,/, разумеется, зависит от поставленной задачи. Так, например, для определения характеристик поляризованного излучения "в точке" нами используются "локальные оценки". Отметим, что алгоритмы численного статистического моделирования естественным образом распараллеливаются путем распределения численных статистических испытаний по отдельным процессорам, поэтому, в связи с ростом мощностей вычислительных систем, их исследование приобретает особое значение.

В силу особенностей рассматриваемого подхода для применения метода Монте-Карло к вычислению линейного функционала J необходимо на интегральный оператор, описывающий перенос излучения, наложить некоторые ограничения, обеспечивающие существование математического ожидания оценки ее несмещенность и конечность дисперсии. В диссертационной работе изучается вопрос о конечности дисперсии стандартной векторной оценки метода Монте-Карло. Однако, даже если дисперсия оценки £ конечна, она может оказаться довольно большой и полученный алгоритм окажется практически неприменимым. В этом случае нами применяются специальные весовые модификации моделирования переноса излучения с поляризацией, приводящие к уменьшению дисперсии оценок метода Монте-Карло.

Основные цели работы.

• Исследование спектрального радиуса оператора, определяющего ковариационную матрицу, т.е. и дисперсию, векторной оценки метода Монте-Карло, построение и реализация алгоритма для его численной оценки.

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

• Разработка алгоритма вычисления матриц Якоби для рассматриваемых итерационных методов и его численная реализация.

• Проведение модельных тестовых расчетов.

Методы исследования базируются на теории переноса излучения, теории интегральных уравнений второго рода и теории весовых методов статистического моделирования.

Научная новизна.

1. Проведено исследование спектрального радиуса р матричио-интегрального оператора Кр, определяющего дисперсию стандартной векторной оценки метода Монте-Карло в задаче о переносе поляризованного излучения.

2. На основе теории положительных операторов аналитически вычислена величина р(5р) = ро, где 5Р оператор, определяющий ковариационную матрицу для задачи переноса в бесконечной однородной среде.

3. Получены теоретические выводы о конечности дисперсий оценок функционалов при использовании различных весовых модификаций, которые являются основанием для применимости метода Монте-Карло.

4. Проведены численные эксперименты по исследованию поведения дисперсий векторных оценок метода Монте-Карло при значениях коэффициента поглощения в среде близких к критическим, т.е. для которых теоретически дисперсия оценки бесконечна.

5. Построен алгоритм оценки спектрального радиуса оператора Кр методом Монте-Карло на основе итераций соответствующей резольвенты.

6. С помощью расчетов разработанным методом, а также приближенно аналитически, на основе возмущения специального функционала, показано, что величина р{Кр) для ограниченной среды приближенно равна произведению значения ро на спектральный радиус оператора, соответствующего переносу излучения без поляризации.

7. Для решения задачи восстановления индикатрисы рассеяния атмосферы по наземным наблюдениям "яркости неба" в "альмукантарате" Солнца предложен новый итерационный метод, представляющий собой комбинацию аддитивного и мультипликативного методов и учитывающий отражение от подстилающей поверхности. Рассматриваемые методы распространены на

случай поляризованного излучения. Дается обоснование предложенного метода.

8. С помощью предложенного комбинированного метода, а также двух ранее известных - аддитивного и мультипликативного методов, проведены тестовые расчеты для однородного плоского слоя. Сравнение результатов, полученных разными методами показало, что

• при увеличении оптической толщины слоя г аддитивный метод расходится, а комбинированный и мультипликативный методы сходятся при сравнительно больших т;

• погрешность восстановления индикатрисы комбинированным методом меньше, чем мультипликативным, если присутствует отражение от подстилающей поверхности;

• учет поляризации излучения молекулами воздуха повышает точность методов.

9. Разработан алгоритм вычисления матриц Якоби для рассматриваемых итерационных методов. С помощью построенного алгоритма проведены расчеты матриц Якоби, позволяющие обосновать сходимость для различных параметров среды.

Достоверность. Достоверность полученных выводов подтверждается анализом разработанных алгоритмов и проведением численных экспериментов на модельных задачах с известным точным решением, а также контрольными расчетами альтернативными методами.

Практическая значимость работы. Исследование спектрального радиуса матрично-интегрального оператора, определяющего ковариационную матрицу стандартной векторной оценки метода Монте-Карло в задаче о переносе поляризованного излучения, позволило сделать выводы о конечности дисперсии этой оценки, которые дают основания к ее применению для решения различных практических задач. Предложен новый эффективный метод для решения задачи восстановления индикатрисы аэрозольного рассеяния по наземным наблюдениям "яркости неба" в альмукантарате Солнца, а также показано, что учет поляризации при решении этой задачи повышает точность восстановления.

Личный вклад соискателя заключается в построении и обосновании вычислительных алгоритмов, разработке на их основе ком-

плекса программ, проведении численных экспериментов и анализе их результатов. Все результаты, представленные в диссертации, получены автором самостоятельно или при ее непосредственном участии.

Апробация работы. Результаты, изложенные в диссертации, докладывались и обсуждались на семинаре Отдела статистического моделирования в физике ИВМиМГ СО РАН (2005 - 2009 гг.), на семинаре Отдела математических задач геофизики ИВМиМГ СО РАН (2009 г.), а также на следующих конференциях:

- Конференции молодых ученых ИВМиМГ СО РАН (2004, 2006,

2008 гг.)

- Всероссийская конференция по вычислительной математике КВМ-2007 (г. Новосибирск, 2007 г.)

- Всероссийская конференция по вычислительной математике КВМ-2009 (г. Новосибирск, 2009 г.)

- Международная конференция по математическим методам в геофизике ММГ-2008 (г. Новосибирск, 2008 г.)

- Шестой Санкт-Петербургский международный семинар по стохастическому моделированию Simulation (г. Санкт-Петербург,

2009 г.).

Публикации. По тематике диссертации автором опубликовано 7 работ (в том числе одна работа без соавторов), среди которых 4 работы в изданиях из списка ВАК [1-4]. Список опубликованных работ помещен в конце автореферата.

Структура и объём работы. Диссертация состоит из введения, четырех глав, разбитых на разделы, заключения и списка литературы. Диссертация изложена на 96 страницах, включает библиографический список из 37 наименований работ, 8 рисунков, 11 таблиц.

КРАТКОЕ СОДЕРЖАНИЕ РАБОТЫ

Во Введении обоснована актуальность темы диссертации, цель и задачи исследований, дается обзор литературы по изучаемым в диссертации вопросам, краткое содержание диссертации по главам и параграфам, приведен перечень результатов, выносимых на защиту.

Первая глава посвящена вопросам теории и алгоритмов метода Монте-Карло и исследованию стандартной векторной оценки метода Монте-Карло.

В разделе 1.1 приводится вводная информация. Вектор-параметр Стокса интенсивности света имеет вид I = (7, <5, ГУ, V) в четырехмерном функциональном пространстве. При этом для параметров Стокса справедливы следующие соотношения: / > 0, /2 > д2 4- и2 ->- V2. Для стационарных задач компоненты стоксовских вектор-функций зависят от пяти переменных: I = 1(х), х € X, где X = Я х 0,7? С И3 - пространство координат, П = {ш 6 И.3 : = 1} С И2 - пространство единичных векторов направлений.

Традиционной математической моделью процесса переноса излучения является интегро-дифференциальное уравнение переноса

{и>\7)Ф + аФ= [ а8Р(ш',ш)Ф(г,ш')сЬ' + Jn

или в операторном виде Ь Ф + сг Ф = Я Ф + ^, где Ф{х) — (4^1 (ж), Ф2{х), Фз{х), Фц{х))т = 1(ж) - вектор-функция плотности потока частиц ("векторных фотонов"), иначе - интенсивности излучения; Р(и>', ш) - матричная функция рассеяния, а = а(г) - полное сечение, а = а8 + ас, ас - сечение поглощения, аа - сечение рассеяния; ^ = {/о1\/о2\1о\1о1))Т ~ вектор-функция плотности распределения источника частиц, /л = (ш,о/) - косинус угла рассеяния. Рассматриваемый вектор-параметр Стокса I после рассеяния преобразуется согласно формуле

1(г,ш)=Р(иЛа;,г)-1(г,сУ),

где Р{ш',ш, г) - фазовая матрица рассеяния, которая определяется соотношением Р(ш',и) — Ь(ж — г2)Д(/и)1<(—г1)/27г. Здесь Ь(-) - специальная матрица поворота, Д(-) - матрица рассеяния, которая, в общем случае, представляется как средневзвешенное от матриц молекулярного и аэрозольного рассеяния:

е( ч = Я„(/х, г)о-д(г) + Дт(м)(тт(г)

где аа и ат - соответственно коэффициенты аэрозольного и молекулярного рассеяния; в - ось локальной сферической системы координат; ¿1 - угол между плоскостью а/, в и плоскостью рассеяния а;,и/; ¿2 - угол между плоскостью рассеяния н плоскостью и

Вектор-функция плотности столкновений Ф(г, ш) с вектор-функцией плотности потока связана соотношением: Ф = (аФ\, а<&2, иФз, &$4)Т = (ФъФ2,Фз, $4)T- Известно, что она удовлетворяет векторно-интегральному уравнение переноса с учетом поляризации:

Ф(х) = У К(х', x)<$>{x')dx' + F(x), или Ф = КФ + F, х

со следующим матричным ядром:

, q(r')e-^'^a(r)P(n,r) , г - г' ^

*(*.*) =-^тгттр-

Алгоритмы метода Монте-Карло основаны на представлении решения приведенного выше векторно-интегралыюго уравнения рядом Неймана. Приведено (несколько более точное, чем в [Марчук Г. И., Михайлов Г. А., Назаралиев М. А. и др., 1976]) обоснование сходимости ряда Неймана на основе мажорантного свойства первой компоненты вектора Стокса.

В разделе 1.2 описано построение векторных оценок метода Монте-Карло для функционалов от решения рассматриваемого векторно-интегрального уравнения второго рода вида

1я = (Ф,Я) = (^,Ф*),

где Ф* - решение сопряженного уравнения.

В разделе 1.3 приводится уравнение для ковариационной матрицы векторной оценки решения сопряженного уравнения. Для оценки этого решения методом Монте-Карло строится векторная случайная величина такая, что Е£х = Ф*(х). Если почленное осреднение ряда для допустимо, то ковариационная матрица E(£x£j) = ip(x) удовлетворяет матричио-интегральному уравнению [Михайлов Г. А., 1987]

J р{х,у)

где А = НФ*Т + Ф*НТ — ННТ, а р(х,у) - переходная плотность моделируемой цепи Маркова.

Обозначим матрично-иптегральпый оператор из этого уравнения через Кр. Если спектральный радиус р(Кр) < 1, то указанное выше осреднение матрицы допустимо и, следовательно, дисперсия

оценки конечна. Через ,5'р обозначен оператор, получаемый из Кр при переходе к бесконечной однородной среде.

В разделе 1.4 дан детализированный вывод системы уравнений для определения величины р(5р), которая приведена в [Михайлов Г. А., 2003] фактически без обоснования. Собственная матрица г/>(°) (а;) оператора 5Р представляется в виде диагональной матрицы с неотрицательными элементами (1,01,01,02) на диагонали. Получена система уравнений на Ао, 01,02:

сц + c21oj = А0 Ci2 + (fi22 + Сзз)01 + С43О2 = 2A0Oi C34G1 + С44О2 = AQ02

где Ci j= / —~—rdu, и доказано, что если существует решение систе-

J Ым) -1

мы с неотрицательными компонентами Ао, ai, 02 , то p(Sp) = Ао- Приводится общий алгоритм решения полученной системы. Здесь ]>2(ц) - моделируемая плотность распределения угла рассеяния, rfj(fi) -элементы матрицы рассеяния R(fi).

В разделе 1.5 сформулирована и доказана теорема об оценке спектрального радиуса оператора Кр через спектральный радиус оператора Sp. Если матрица R не зависит от г, то

р(Кр) < p{Sp)qQ.

где go — sup /--^-гДе QiQi ~ физическая и моделируемая вероятности выживания частицы, px(l] г, uj),p^\l\ г, uj) - физическая и моделируемая плотности распределения длины свободного пробега. Для прямого моделирования qo соответствует верхней границе выживания в среде. Выведены условия на коэффициент поглощения, обеспечивающие конечность дисперсии оценок функционалов, для различных способов моделирования процесса переноса.

Проведенные расчеты показывают, что для реальных типов рассеяния p(Kp)/p(Sp) ~ р{Ьр), где Lp - скалярный оператор, соответствующий переносу излучения без поляризации. Полученное приближенное равенство обосновано приближенно аналитически на ос-

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

В разделе 1.6. приводится алгоритм вычисления спектрального радиуса оператора Кр методом Монте-Карло на основе итераций резольвенты.

Вторая глава посвящена рассмотрению задачи восстановления индикатрисы рассеяния по наземным наблюдениям яркости неба в альмукантарате Солнца и методой ее решения.

Раздел 2.1 содержит вводную информацию о задаче восстановления индикатрисы аэрозольного рассеяния атмосферы по наземным наблюдениям "яркости неба" в альмукантарате Солнца, то есть в различных направлениях, составляющих с зенитом тот же угол что п направление на Солнце. Поле излучения в атмосфере формирз'ется одно- и многократно рассеянным солнечным светом, а также отраженным от подстилающей поверхности излучением. В "приближении однократного рассеяния" наблюдаемые значения яркости пропорциональны соответствующим значениям индикатрисы.

В разделе 2.2 приводится постановка задачи, алгоритм построения индикатрисы рассеяния д(у) по вектору д значений в узловых точках и используемые условные обозначения: ^(<7), ^ (¿г), -

соответственно вектора полной яркости, яркости однократного рассеяния (при А = 0) и яркости от "альбедного источника" ("подсветка"), наблюдаемых в альмукантарате Солнца при индикатрисе рассеяния д(ц), представляющей собой средневзвешенное от известной рэлеевской индикатрисы дт((м) и аэрозольной индикатрисы да(м)> построенной по вектору значений да. ^(<?) = Р(д) — ^(д) — - яркость многократного рассеяния при отсутствии отражения от Земли; Р* = Р{д*) - вектор измеренных яркостей.

Разделе; 2.3 содержит методы восстановления индикатрисы по наземным наблюдениям яркости небе в альмукантарате Солнца. Для оценки индикатрисы было построено несколько итерационных алгоритмов, в которых с помощью математического моделирования последовательно уточняются значения индикатрисы на основании информации об угловом распределении поля яркости на подстилающей поверхности и предположения, что доля однократно рассеянного излучения достаточно велика:

Аддитивный метод восстановления [Аитюфеев В. С., Михайлов Г. А., Лнфшпц Г. Ш., Иванов А. И., 1980], определяемый формулой

и уточненный (сравнительно с [Михайлов Г. А., 1978; Антюфеев В. С., Назаралиев М. А., 1988]) мультипликативный метод, реализуемый по формуле

Эти методы были обобщены на задачу переноса излучения с учетом поляризации и был предложен новый комбинированный метод:

9

(к)

Кратко описан процесс восстановления индикатрисы с помощью приведенных методов.

Раздел 2.4 содержит описание локальной векторной оценки метода Монте-Карло, используемой для решения "прямой" задачи, т.е. для вычисления вектора Стокса общего излучения и излучения "подсветки" при заданной матрице рассеяния.

В разделе 2.5 выведены формулы для расчета матриц Якоби для операторов перехода рассматриваемых итерационных методов.

Раздел 2.6 содержит утверждение о том, что сходимость последовательных приближений к решению следует из того, что спектральный радиус матрицы Якоби оператора перехода меньше 1. Приводится обоснование сходимости комбинированного метода при малых оптической толще атмосферы и коэффициенте отражения от Земли.

Третья глава посвящена описанию применяемых вычислительных алгоритмов.

Раздел 3.1 описывает общий весовой алгоритм моделирования траектории векторного фотона и вычисления функционалов от решения уравнения переноса поляризованного излучения. Подраздел 3.1.1 содержит алгоритм моделирования отражения от подстилающей поверхности, в подразделе 3.1.2 дано описание применяемых весовых модификаций моделирования.

В разделе 3.2 описан алгоритм вычисления интенсивности излучения и ее производной в альмукантарате с помощью локальной векторной оценки метода Монте-Карло.

В разделе 3.3 приводится алгоритм восстановления индикатрисы итерационными методами, описанными во второй главе, с учетом поляризации излучения.

Раздел 3.4 содержит алгоритм расчета матриц Якоби для итерационных методов но формулам, полученным в разделе 2.4.

В четвертой главе собраны численные результаты, полученные автором с использованием разработанных алгоритмов, п приведен их анализ.

В разделе 4.1 приведены результаты вычисления спектрального радиуса и первого собственного элемента оператора 5Р для молекулярного и рассматриваемой модели аэрозольного рассеяния.

В разделе 4.2 анализируются результаты оценки первого собственного числа оператора Кр с помощью алгоритма основанного на итерациях резольвенты. Одновременно с ним вычислялось первое собственное число оператора скалярного оператора Ьр, соответствующего переносу излучения без поляризации. Расчеты, проведенные для сред различной оптической толщины с молекулярным и аэрозольным рассеянием, показали что полученные значения р(Кр)/р(Ьр) статистически незначимо отличаются от значений р(5р), полученных аналитически и достаточно точно оцениваются с использованием уже первой итерации резольвенты..

В разделе 4.3 приводятся вычисленные ограничения на коэффициент поглощения в среде, при котором дисперсия стандартной векторной оценки конечна, показывающие, что использование весовых модификаций моделирования позволяет ослабить условия на коэффициент поглощения. Исследуется поведение дисперсии оценки векторной плотности столкновений в бесконечной однородной среде для различных значений коэффициента выживания. Результаты расчетов показали, что при коэффициенте поглощения больше критического оценка функционала сходится и ее дисперсия остается постоянной при увеличении количества траекторий, а при близком к критическому и меньше наблюдается рост дисперсии с увеличением числа траекторий, хотя оценка функционала остается достаточно точной.

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

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

В разделе 4.5 приведен анализ и сравнение результатов восстановления индикатрисы рассеяния различными методами. Расчеты проводились с помощью комбинированного, мультипликативного и аддитивного методов, при этом для каждого метода был реализован вариант моделирования переноса излучения с поляризацией и без нее.

В разделе 4.6 приводится анализ результатов расчета матриц Якоби для рассматриваемых итерационных методов. Расчеты объясняют численные результаты восстановления индикатрисы, полученные в предыдущем разделе. Для контроля матрицы Якоби были рассчитаны разностным методом и методом зависимых испытаний, результаты оказались близкими к полученным с помощью используемого алгоритма.

Заключение содержит перечень основных результатов диссертационной работы.

ЗАКЛЮЧЕНИЕ

Сформулируем основные результаты диссертационной работы.

1. Проведено исследование спектрального радиуса р матрично-интегрального оператора Кр, определяющего ковариационную матрицу, т.е. и дисперсию, стандартной векторной оценки метода Монте-Карло в задаче о переносе поляризованного излучения. Построен алгоритм оценки спектрального радиуса оператора Кр методом Монте-Карло на основе итераций соответствующей резольвенты. Показано, что величина р(Кр) для ограниченной среды приближенно равна произведению спектральных радиусов оператора, соответствующего переносу излучения без поляризации п оператора, соответствующего переносу излучения в бесконечной однородной среде, который вычислен аналитически на основе теории положительных операторов.

2. Для решения задачи восстановления индикатрисы рассеяния атмосферы по наземным наблюдениям яркости поляризованного излучения в альмукантарате Солнца предложен новый итерационный метод, представляющий собой комбинацию двух ранее известных в скалярном варианте аддитивного и мультипликативного методов и учитывающий отражение от подстилающей поверхности. Дано теоретическое обоснование сходимости предложенного метода, подтвержденное численными расчетами. В связи с этим разработан и обоснован алгоритм вычисления матриц Якоби для рассматриваемых итерационных методов.

3. На основе рассмотренных методов разработан комплекс программ, с помощью которого проведены тестовые расчеты для однородного плоского слоя. Сравнение результатов, полученных разными методами показало преимущества предложенного метода и целесообразность учета поляризации излучения. Были также проведены расчеты матриц Якоби и дан их анализ, позволивший обосновать сходимость методов для различных параметров среды.

Пользуясь случаем, автор выражает искреннюю благодарность своему научному руководителю к.ф.-м.н., с.н.с. Сергею Анатольевичу Ухинову за постоянное внимание и плодотворное руководство работой, а также члену-корреспонденту РАН Геннадию Алексеевичу Михайлову за ценные консультации.

СПИСОК РАБОТ ПО ТЕМЕ ДИССЕРТАЦИИ

Публикации в журналах списка ВАК

[1] Chimaeva A.S., Mikhailov G.A. Study of polarization estimates variance by the Monte Carlo method // Russian Journal of Numerical analysis and Mathematical Modelling. - 2006. - Vol.20, №3 - P. 305-317.

[2] Михайлов Г.А., Ухинов С.А.,Чимаева A.C. Дисперсия стандартной векторной оценки метода Монте-Карло в теории переноса поляризованного излучения //Журн. вычисл. математики и мат. физики. - 2006. Т. 46, № 11. С. 2199-2212.

[3] Трачева Н.В., Ухинов С.А., Чимаева А.С. Расчет параметров временной асимптотики выходящего из полубесконечного слоя поляризованного излучения // Вычисл. технологии. 2008. Т. 13, Специальный выпуск 4. С. 120-130.

[4] Михайлов Г.А., Ухинов С.А., Чимаева А.С. Алгоритмы метода Монте-Карло для восстановления индикатрисы рассеяния с учетом поляризации.// Доклады Академии Наук, 2008, том 423, №2, с. 161-164.

Прочие публикации

[5| Чимаева А.С. Исследование дисперсии оценок поляризации методом Монте-Карло. // Труды конференции молодых ученых ИВМиМГ СО РАН. - Новосибирск, 2004. - С. 242-248.

[6] Чимаева А. С. Численное исследование дисперсии стандартной весовой оценки метода Монте-Карло в теории переноса поляризованного излучения// Труды конференции молодых ученых ИВМиМГ СО РАН. - Новосибирск, 2006. - С. 256-262.

[7] Chimaeva A.S., Ukhiniov S.A. Monte Carlo algorithms for reconst ructing a scattering phase function with polarization taken into account // Proceedings of the 6th St. Petersburg Workshop on Simulation (2009), pp. 122-126

Чимаева Анна Сергеевна

ВЕСОВЫЕ АЛГОРИТМЫ СТАТИСТИЧЕСКОГО МОДЕЛИРОВАНИЯ ПЕРЕНОСА ПОЛЯРИЗОВАННОГО ИЗЛУЧЕНИЯ И РЕШЕНИЕ ЗАДАЧИ ВОССТАНОВЛЕНИЯ ИНДИКАТРИСЫ РАССЕЯНИЯ

05.13.18 - математическое моделирование, численные методы и комплексы программ

Автореферат

диссертации на соискание ученой степени кандидата физико-математических наук

Лицензия ИД N0 02202 от 30 июня 2000 г.

Подписано в печать 05.10.2009 г. Формат бумаги 60 X 84х/16 ОЕ Тираж 100 экз.

Объем 1,0 п.л. 0,9 уч.-изд. л.

Заказ N0 157

ООО "Омега Принт", Новосибирс:к-90, пр. Лаврентьева, 6

Введение

1. Математическая модель переноса излучения с поляризацией и соответствующие оценки метода Монте-Карло

1.1. Система уравнений переноса излучения с поляризацией

1.2. Векторная оценка метода Монте-Карло.

1.3. Уравнение для ковариационной матрицы векторной оценки сопряженного решения (с оператором = ЬР3Р).

1.4. Спектральный радиус и первый собственный элемент оператора

1.5. Оценка спектрального радиуса оператора через спектральный радиус оператора

1.6. Оценка первого собственного числа оператора Кр методом Монте-Карло.

2. Восстановление индикатрисы рассеяния с учетом поляризации

2.1. Вводная информация

2.2. Постановка задачи.

2.3. Итерационные методы для решения обратной задачи

2.4. Решение "прямой" задачи методом Монте-Карло.

2.5. Матрицы Якоби для итерационных методов.

2.6. Исследование сходимости методов.

3. Вычислительные алгоритмы 56 3.1. Моделирование переноса поляризованного излучения и стандартная векторная оценка метода Монте-Карло.

3.1.1. Моделирование отражения от подстилающей поверхности

3.1.2. Весовые модификации моделирования

3.2. Алгоритм вычисления интенсивности излучения и ее производных в альмукантарате Солнца

3.3. Алгоритм восстановления индикатрисы итерационными методами

3.4. Алгоритм расчета матриц Якоби для итерационных методов

4. Решение модельных и прикладных задач, численные результаты

4.1. Вычисление спектрального радиуса и первого собственного элемента оператора

4.2. Оценка первого собственного числа оператора Кр в средах с различными оптическими толщинами.

4.3. Исследование дисперсии стандартной векторной оценки метода Монте-Карло при значениях коэффициента выживания, близких к критическому.

4.4. Оценка и исследование интенсивностей, наблюдаемых в альмукантарате Солнца

4.5. Восстановление индикатрисы рассеяния различными методами; сравнение результатов

4.6. Вычисление матрицы Якоби, се нормы и спектрального радиуса

Имеется целый ряд физических проблем, требующих достаточно точного решения задач теории переноса излучения в атмосфере с учетом поляризации. Это, прежде всего, задачи интерпретации оптических наблюдений.

Математическое описание распространения поляризованного излучения, изложенное в работах [5, 29], предоставляет нам удобный инструмент для исследования этого процесса. Базовыми при этом являются вектор-функции Стокса, характеризующие свойства излучения в каждой конкретной точке фазового пространства, и интегрально-дифференциальное уравнение переноса описывающее процесс переноса.

Интенсивность и состояние поляризации излучения определяются чстырехкомпонентной вектор-функцией Стокса, компоненты которой определяют в совокупности интенсивность, степень поляризации, плоскость поляризации и степень эллиптичности излучения.

Процесс переноса излучения с поляризацией можно описать интегральным уравнением второго рода (см., например, [6, 14, 16, 22, 25]), оператор которого, в силу физических особенностей задачи, оставляет инвариантным множество вектор-функций Стокса. Исследованию свойств подобных операторов посвящены, например, работы [4, 11].

Изучению уравнения переноса излучения посвящена обширная литература (см., например, [8, 13, 24, 32]). В ней содержатся математические постановки задач теории переноса и вывод интегро-дифференциального уравнения переноса. Обоснование условий существования собственных значений дается в работах [7, 32]. Численные методы решения соответствующих задач рассмотрены в известных монографиях [3, 12, 13]. Среди них существенное место занимает метод статистического моделирования (метод Монте-Карло) [6, 12, 18], т.к. уравнение переноса трудно разрешимо классическими методами вычислительной математики, если учитываются реальные индикатрисы и неоднородность среды. Во многих случаях практически это можно осуществить методом Монте-Карло.

Процесс распространения света можно рассматривать как случайную марковскую цепь столкновений фотонов с веществом, которые приводят либо к рассеянию (с пересчетом вектора Стокса), либо к поглощению фотонов. Метод Мопте-Карло заключается в моделировании траекторий этой цепи на ЭВМ и вычислении статистических оценок для искомых функционалов. Построение случайных траекторий для физической модели процесса принято называть прямым моделированием. При этом веса не используются и дисперсии оценок метода Монте-Карло в скалярном случае заведомо конечны [14]. Указанный выше пересчет вектора Стокса уже предполагает использование матричного веса. В связи с этим в [14] были построены и предварительно исследованы общие матрично-весовые алгоритмы для решения систем интегральных уравнений теории переноса излучения с учетом поляризации.

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

Рассматриваемая математическая модель позволяет ставить достаточно большое множество практически интересных задач, для решения которых может быть эффективно применен метод Монте-Карло. Традиционный способ его использования заключается в следующем. Рассматривается некоторый линейный функционал ,7 от решения уравнения переноса, для него строится стандартная весовая оценка статистического моделирования математическое ожидание которой и дает нам искомое значение функционала.

Конкретный вид функционала ,7, разумеется, зависит от поставленной задачи. Так, например, для определения характеристик поляризованного излучения "в точке" нами используются "локальные оценки"".

Заметим, что для применения метода Монте-Карло к вычислению линейного функционала 3 необходимо на интегральный оператор, описывающий перенос излучения, наложить некоторые ограничения, обеспечивающие существование математического ожидания оценки ее несмещенность и конечность дисперсии. Однако, даже если дисперсия оценки £ конечна, она может оказаться довольно большой и полученный алгоритм окажется практически неприменимым. В этом случае нами применяются специальные весовые модификации моделирования переноса излучения с поляризацией, приводящие к уменьшению дисперсии оценок метода Монте-Карло.

Основными целями диссертационной работы являются:

• Исследование спектрального радиуса оператора, определяющего ковариационную матрицу векторной оценки метода Монте-Карло, построение и реализация алгоритма для его численной оценки.

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

• Разработка алгоритма вычисления матриц Якоби для рассматриваемых итерационных методов и его численная реализация.

• Проведение модельных тестовых расчетов.

Первая часть диссертационной работы посвящена изучению вопроса о конечности дисперсии соответствующих стандартных векторных оценок метода Монте-Карло, что необходимо для корректного построения требуемых доверительных интервалов. С этой целью рассматривается система интегральных уравнений, определяющая ковариационную матрицу весовой векторной оценки. С помощью полуэвристических аналитических выкладок и численных оценок на основе итерирования резольвенты показано,'что спектральный радиус соответствующего матрично-интегрального оператора близок к произведению аналогичного спектрального радиуса для бесконечной среды, который может быть вычислен аналитически, на просто оцениваемый спектральный радиус скалярного интегрального оператора, связанного с "единичной" матрицей рассеяния. Этот вывод получен как для простейшего релеевского, т.е. молекулярного, рассеяния, так и для реального аэрозольного рассеяния. Полученные результаты с большой степенью надежности показывают, что для реальных сред стандартные векторные оценки метода Монте-Карло имеют конечную дисперсию.

Вариации яркости неба главным образом определяются изменчивостью аэрозольной компоненты атмосферы, которая обладает сильной рассеивающей способностью. Теоретическое исследование проблем переноса излучения в такой среде невозможно без знания ее оптических параметров. Одним из них является аэрозольная индикатриса рассеяния. В диссертационной работе рассматривается задача восстановления индикатрисы аэрозольного рассеяния атмосферы д(/л) = гц(д) по наземным наблюдениям "яркости неба" в "альмукантарате" Солнца, то есть в различных направлениях, составляющих с зенитом тот же угол в8, что и направление на Солнце [1],[2],[17]. Поле излучения в атмосфере формируется одно- и многократно рассеяным солнечным светом, а также отраженным от подстилающей поверхности излучением. В "приближении однократного рассеяния" наблюдаемые значения яркости пропорциональны соответствующим значениям индикатрисы [1]. Для оценки индикатрисы было построено несколько итерационных алгоритмов [1],[2],[17], в которых с помощью математического моделирования последовательно уточняются значения индикатрисы на основании информации об угловом распределении поля яркости на подстилающей поверхности и предположения, что доля однократно рассеянного излучения достаточно велика. В диссертационной работе два таких алгоритма обобщены на случай поляризованного излучения, а также построен новый алгоритм такого типа, приведено его обоснование и реализация методом Монте-Карло.

Далее следует краткое содержание диссертации по главам. Первая глава посвящена вопросам теории и алгоритмов метода Монте-Карло и исследованию стандартной векторной оценки метода Монте-Карло.

В разделе 1.1 приводится вводная информация. Вектор-параметр т

Стокса интенсивности света имеет вид I = [I, Q,U,V) в четырехмерном функциональном пространстве. При этом для параметров Стокса справедливы следующие соотношения: I > О, I2 > Q2 + U2 + V2. Для стационарных задач компоненты стоксовских вектор-функций зависят от пяти переменных: I = I(x), х е X, где X = RxQ, R С R3-пространство координат, Q = {и G R3||oj| = 1} С R3 - пространство единичных векторов направлений.

Традиционной математической моделью процесса переноса излучения является интегро-дифференциальное уравнение переноса шУФ+аФ= / а3Р(ш', и) Ф(г, и') du' + f0(r, си), Jo. или в операторном виде

L Ф + стФ = БФ + f0, где

Ф(х) = (Ф^х), Ф2(х), Фз(х), Ф^х))т = 1(х)

- вектор-функция плотности потока частиц ("векторных фотонов"), иначе -интенсивности излучения; Р(ш', ш) - матричная функция рассеяния, a = сг(г) - полное сечение, a = crs 4- <тс, <тс - сечение поглощения, as - сечение рассеяния; fo = (fo'\fo'\fo3\fo^)T ~ вектор-функция плотности распределения источника частиц, /л = (о;, а/) - косинус угла рассеяния. Рассматриваемый вектор-параметр Стокса I после рассеяния преобразуется согласно формуле

1(г,ы) = Р(о/,и;,г)-1(г,а/), где P(oj', со, г) фазовая матрица рассеяния, которая определяется соотношением

P(w\w) = L(tt - i2)R(fJb)L(-i1)/2тг,

•) - специальная матрица поворота, R(-) - матрица рассеяния, которая, в общем случае, представляется как средневзвешенное от матриц молекулярного и аэрозольного рассеяния:

R( J.) = Rar)a°(r) +

J cra(r) + сгт(г) s - ось локальной сферической системы координат; - угол между плоскостью о/,s и плоскостью рассеяния а;,и/;г2 - угол между плоскостью рассеяния и, и' и плоскостью u,s.

Вектор-функция плотности столкновений Ф(г, и) с вектор-функцией плотности потока связана соотношением: Ф = (аФ1,о-Ф2,сгФз,сгФ4)т — (Фь Ф2, Фз, Ф4)Т- Известно, что она удовлетворяет векторно-интегральному уравнение переноса с учетом поляризации:

Ф(ж) = J К(х', хЩх')йх' + F(x), или Ф = КФ + F, X со следующим матричным ядром: д(г/)е~т^г''г^(г(г)Р(/и, г) „ / г - г'

К(х',х) =-:-ytj-— х 5[ш — г — г'|2 V |г — г'

Алгоритмы метода Мопте-Карло основаны на представлении решения приведенного выше векторно-интегрального уравнения рядом Неймана. Такое представление имеет место, если норма оператора К (или какой-нибудь его степени Кп°) меньше единицы. Приведено (несколько более точное, чем в [14]) обоснование сходимости ряда Неймана на основе мажорантного свойства первой компоненты вектора Стокса.

В разделе 1.2 описано построение векторных оценок метода Монте-Карло от решения рассматриваемого матричного интегрального уравнения второго рода. Общий алгоритм метода Монте-Карло для оценки таких функционалов приводится в работах [14, 16, 18, 28].

В разделе 1.3 приводится уравнение для ковариационной матрицы векторной оценки сопряженного решения. Для оценки решения методом Монте-Карло строится векторная случайная величина такая, что Е^ = Ф*(ж). Если почленное осреднение ряда для допустимо, то ковариационная матрица = ф(х) удовлетворяет матричноинтегральному уравнению (см. [19]) ф(х) = А(х) + I --<*у, где А — НФ*Т + Ф*НТ — ННТ, а у) -- переходная плотность моделируемой цепи Маркова.

Обозначим матрично-интегральный оператор из этого уравнения через Кр. Если спектральный радиус р{Кр) < 1, то указанное выше осреднение матрицы при Нт = (1,0, 0, 0) допустимо. Через обозначен оператор, получаемый из К^ при переходе к бесконечной однородной среде.

В разделе 1.4 дан детализированный вывод системы уравнений для определения величины р(Зр), которая приведена в [15] фактически без обоснования. Собственная матрица ф^ (ш) оператора представляется в виде диагональной матрицы с неотрицательными элементами (1, ах, ах, <22) на диагонали. Получена система уравнений на До, <21, <22:

СЦ + С21Й1 = А0

С12 + (С22 + с33)<21 + с43а2 = 2А0а1 Сма1 + С44«2 = АО«2 где сц = / , .а/2, и доказано, что если существует решение системы -1 с неотрицательными компонентами Ад, <21,02 , то р{8р) = Ад. Приводится общий алгоритм решения полученной системы. Здесь ^(а4) - моделируемая плотность распределения угла рассеяния, -элементы матрицы рассеяния

В разделе 1.5 сформулирована и доказана теорема об оценке спектрального радиуса оператора Кр через спектральный радиус оператора

Если матрица Я не зависит от г, то р(Кр) < р(Зр)дооо д2 где = вир где д, дх - и моделируемая о вероятности выживания частицы, г, ш),р{'(1; г, о;) - физическая и моделируемая плотности распределения длины свободного пробега. Для прямого моделирования до соответствует верхней границе выживания в среде. Выведены условия на коэффициент поглощения, обеспечивающие конечность дисперсии оценок функционалов, для различных способов моделирования процесса переноса.

Проведенные расчеты (см. раздел 4.2) показывают, что для реальных типов рассеяния р(Кр)/р(Зр) ~ р{Ьр)^ где Ьр - скалярный оператор, соответствующий переносу излучения без поляризации. Полученное приближенное равенство обосновано приближенно аналитически на основе возмущения специального функционала при переходе от бесконечной среды к конечному слою.

В разделе 1.6. приводится алгоритм вычисления спектрального радиуса оператора Кр методом Монте-Карло на основе итераций резольвенты.

Вторая глава посвящена рассмотрению задачи восстановления индикатрисы по наземным наблюдениям яркости неба в альмукантарате Солнца и методов ее решения.

Раздел 2.1 содержит вводную информацию о задаче по восстановления индикатрисы аэрозольного рассеяния атмосферы по наземным наблюдениям "яркости неба" в "альмукантарате" Солнца, то есть в различных направлениях, составляющих с зенитом тот же угол в8, что и направление на Солнце. Поле излучения в атмосфере формируется одно-и многократно рассеянным солнечным светом, а также отраженным от подстилающей поверхности излучением. В "приближении однократного рассеяния" наблюдаемые значения яркости пропорциональны соответствующим значениям индикатрисы.

В разделе 2.2 приводится постановка задачи, алгоритм построения индикатрисы рассеяния д(ц) по вектору д значений в узловых точках и используемые условные обозначения: Р(д), Р\{д), -Рк(д) ~ соответственно вектора полной яркости, яркости однократного рассеяния (при А = 0) и яркости от "альбедного источника" ("подсветка"), наблюдаемых в альмукантарате Солнца при индикатрисе рассеяния д(р), представляющей собой средневзвешенное от известной рэлеевской индикатрисы дт(ц) и аэрозольной индикатрисы да(д), построенной по вектору значений да.

9) — Р{д) ~~ Ра{э) ~ ^г(д) - яркость многократного рассеяния при отсутствии отражения от Земли; Р* = Р{д*) - вектор измеренных яркостей.

Разделе 2.3 содержит методы восстановления индикатрисы по наземным наблюдениям яркости небе в альмукантарате Солнца. Для оценки индикатрисы было построено несколько итерационных алгоритмов, в которых с помощью математического моделирования последовательно уточняются значения индикатрисы па основании информации об угловом распределении поля яркости на подстилающей поверхности и предположения, что доля однократно рассеянного излучения достаточно велика: Аддитивный метод восстановления [1], определяемый формулой д(к) = С-1[Е* р2(д^) - = Оа(1(д^), и уточненный (сравнительно с [17], /2]) мультипликативный метод, реализуемый по формуле д(к) = ) . = СтШ{д{к 1}).

Эти методы были обобщены на задачу переноса излучения с учетом поляризации и был предложен новый комбинированный метод:

Кратко описан процесс восстановления индикатрисы с помощью приведенных методов.

Раздел 2.4 содержит описание локальной векторной оценки метода Монте-Карло, используемой для решения "прямой" задачи, т.е. для вычисления вектора Стокса общего излучения и излучения "подсветки" при заданной матрице рассеяния.

В разделе 2.5 выведены формулы для расчета матриц Якоби для операторов перехода рассматриваемых итерационных методов.

Раздел 2.6 содержит утверждение о том, что сходимость последовательных приближений к решению следует из того, что спектральный радиус матрицы Якоби оператора перехода меньше 1. Приводится обоснование сходимости комбинированного метода при малых оптической толще атмосферы и коэффициенте отражения от Земли.

Третья глава посвящена описанию применяемых вычислительных алгоритмов.

Раздел 3.1 описывает общий весовой алгоритм моделирования траектории векторного фотона и вычисления функционалов от решения уравнения переноса поляризованного излучения. Подраздел 3.1.1 содержит алгоритм моделирования отражения от подстилающей поверхности, в подразделе 3.1.2 дано описание применяемых весовых модификаций моделирования.

В разделе 3.2 описан алгоритм вычисления интенсивности излучения и ее производной в альмукантарате с помощью локальной векторной оценки метода Монте-Карло.

В разделе 3.3 приводится алгоритм восстановления индикатрисы итерационными методами, описанными во второй главе, с учетом поляризации излучения.

Раздел 3.4 содержит алгоритм расчета матриц Якоби для итерационных методов по формулам, полученным в разделе 2.4.

В четвертой главе собраны численные результаты, полученные автором с использованием разработанных алгоритмов, и приведен их анализ.

В разделе 4.1 приведены результаты вычисления спектрального радиуса и первого собственного элемента оператора Зр для молекулярного и рассматриваемой модели аэрозольного рассеяния.

В разделе 4.2 анализируются результаты оценки первого собственного числа оператора Кр с помощью алгоритма основанного на итерациях резольвенты. Одновременно с ним вычислялось первое собственное число оператора скалярного оператора Ьр, соответствующего переносу излучения без поляризации. Расчеты, проведенные для сред различной оптической толщины с молекулярным и аэрозольным рассеянием, показали что полученные значения р(К.р)/р(Ьр) статистически незначимо отличаются от значений р(5,р), полученных аналитически и достаточно точно оцениваются с использованием уже первой итерации резольвенты.

В разделе 4.3 приводятся вычисленные ограничения на коэффициент поглощения в среде, при котором дисперсия стандартной векторной оценки конечна, показывающие, что использование весовых модификаций моделирования позволяет ослабить условия на коэффициент поглощения. Исследуется поведение дисперсии оценки векторной плотности столкновений в бесконечной однородной среде для различных значений коэффициента выживания. Результаты расчетов показали, что при коэффициенте поглощения больше критического оценка функционала сходится и ее дисперсия остается постоянной при увеличении количества траекторий, а при близком к критическому и меньше наблюдается рост дисперсии с увеличением числа траекторий, хотя оценка остается достаточно точной.

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

Расчеты показали, что в случае чисто аэрозольного рассеяния поляризация практически не влияет на интенсивность излучения, в то время как дня молекулярного рассеяния влияние поляризации более существенно, причем разница интенсивности поляризованного и скалярного излучений увеличивается с увеличением оптической толщины слоя и уменьшении альбедо.

Проанализированы вклады в общую яркость. яркостей однократно рассеянного, многократно рассеянного излучения и излучения "подсветки". Интенсивность рэлеевского рассеяния симметрична в силу симметрии индикатрисы и вклады различных видов рассеяния незначительно меняются в зависимости от угла наблюдения. Индикатриса же аэрозольного рассеяния сильно вытянута вперед, поэтому интенсивность излучения при небольших углах наблюдения значительно больше, чем при углах наблюдения больших 90°. Однократно рассеянное излучение дает вклад в основном в направления наблюдения "вперед", а альбедное - в направления "назад". Вклад однократно рассеянного излучения обратно пропорционален оптической толщине слоя, то есть при ее увеличении доля однократно рассеянного излучения уменьшается, а доли многократно рассеянного и альбедного излучения увеличиваются.

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

Сравнение методов показало, что при наличии отражения от подстилающей поверхности комбинированный метод снижает погрешность приближения индикатрисы, по сравнению с мультипликативным. Также сравнение расчетов с учетом поляризации и без показывает, что учет поляризации значительно повышает точность вычислений при достаточно большой доле молекулярного рассеяния, если же аэрозольного рассеяния значительно больше, чем молекулярного, то учет поляризации практпчески не влияет на точность вычислений, т.к. поляризацию аэрозолем мы не учитываем, а поляризация молекулами воздуха в этом случае мала. Аддитивный метод обеспечивает лучшую сходимость итераций при небольшой оптической толщине, однако при ее увеличении метод начинает расходиться, в отличии от комбинированного и мультипликативного методов, которые сходятся и при относительно больших оптических толщинах.

Аддитивный метод обеспечивает лучшую сходимость итераций при небольшой оптической толщине, однако при ее увеличении метод начинает расходиться, в отличии от комбинированного и мультипликативного методов, которые сходятся и при достаточно больших оптических толщинах.

В разделе 4.6 приводится анализ результатов расчета матриц Якоби для рассматриваемых итерационных методов. Расчеты проводились для задачи, описанной в предыдущем разделе. Для сходимости последовательных приближений = к д* достаточно, чтобы спектральный радиус матрицы Якоби С(д*) был меньше 1. Результаты расчетов объясняют численные результаты восстановления индикатрисы. Спектральный радиус матрицы Якоби комбинированного метода при всех рассматриваемых оптических толщинах меньше, чем мультипликативного, причем, чем меньше оптическая толщина слоя, тем больше различие спектральных радиусов.

Заметтм, что с увеличением толщины аэрозольного слоя норма и спектральный радиус матриц Якоби увеличивается для всех методов, однако для комбинированного и мультипликативного методов спектральный радиус остается меньше 1, а для аддитивного метода спектральный радиус становится больше 1, из чего следует возможность расходимости метода. Для малых оптических толщин, когда велика доля молекулярного рассеяния, спектральный радиус матрицы Якоби для аддитивного метода значительно меньше, чем для комбинированного и мультипликативного, то есть в этом случае аддитивный метод может показывать более быструю сходимость.

Расчеты для малых оптических толщин показали что нормы и спектральные радиусы матриц Якоби также малы.

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

Заключение содержит перечень основных результатов диссертационной работы.

Основные результаты диссертации опубликованы в работах [30, 33, 31, 26, 21, 34, 20].

На защиту выносятся следующие основные результаты:

1. Проведено исследование спектрального радиуса р матрично-интегрального оператора Кр, определяющего ковариационную матрицу, т.е. дисперсию, стандартной векторной оценки метода Монте-Карло в задаче о переносе поляризованного излучения. Построен алгоритм оценки спектрального радиуса оператора Кр методом Монте-Карло на основе итераций соответствующей резольвенты. Показано, что величина р(Кр) для ограниченной среды приближенно равна произведению спектральных радиусов оператора, соответствующего переносу излучения без поляризации и оператора, соответствующего переносу излучения в бесконечной однородной среде, который вычислен аналитически на основе теории положительных операторов.

2. Для решения задачи восстановления индикатрисы рассеяния атмосферы по наземным наблюдениям яркости поляризованного излучения в альмукантарате Солнца предложен новый итерационный метод, представляю!ций собой комбинацию двух ранее известных в скалярном варианте аддитивного и мультипликативного методов и учитывающий отражение от подстилающей поверхности. Дано теоретическое обоснование сходимости предложенного метода, подтвержденное численными расчетами. В связи с этим разработай и обоснован алгоритм вычисления матриц Якоби для рассматриваемых итерационных методов.

3. На основе рассмотренных методов разработан комплекс программ, с помощью которого проведены тестовые расчеты для однородного плоского слоя. Сравнение результатов, полученных разными методами показало преимущества предложенного метода и целесообразность учета поляризации излучения. Были также проведены расчеты матриц Якоби и их анализ, позволивший обосновать сходимость методов для различных параметров среды.

Результаты, изложенные в диссертационной работе, докладывались и обсуждались на семинаре Отдела статистического моделирования в физике ИВМиМГ СО РАН, на семинаре Отдела математических задач геофизики ИВМиМГ СО РАН (2009 г.), па конференциях молодых ученых ИВМиМГ СО РАН (2004, 2006, 2008 гг.), на Всероссийской конференции по вычислительной математике КВМ-2007 и КВМ-2009 в г. Новосибирске (2007,2009 гг.), на Международной Конференции по Математическим Методам в Геофизике ММГ-2008 в г. Новосибирске (2008 г.), на Шестом Санкт-Петербургском международном семинаре по стохастическому моделированию Simulation (2009 г.).

Пользуясь случаем, автор выражает искреннюю благодарность своему научному руководителю к.ф.-м.н. Сергею Анатольевичу Ухинову за постоянное внимание и плодотворное руководство работой, а также члену-корреспонденту РАН Геннадию Алексеевичу Михайлову за ценные консультации.

Заключение

Сформулируем основные результаты диссертационной работы.

1. Проведено исследование спектрального радиуса р матрично-интегрального оператора Кр, определяющего ковариационную матрицу, т.е. дисперсию, стандартной векторной оценки метода Монте-Карло в задаче о переносе поляризованного излучения. Построен алгоритм оценки спектрального радиуса оператора Кр методом Монте-Карло на основе итераций соответствующей резольвенты. Показано, что величина р(Кр) для ограниченной среды приближенно равна произведению спектральных радиусов оператора, соответствующего переносу излучения без поляризации и оператора, соответствующего переносу излучения в бесконечной однородной среде, который вычислен аналитически на основе теории положительных операторов.

2. Для решения задачи восстановления индикатрисы рассеяния атмосферы по наземным наблюдениям яркости поляризованного излучения в альмукантарате Солнца предложен новый итерационный метод, представляющий собой комбинацию двух ранее известных в скалярном варианте аддитивного и мультипликативного методов и учитывающий отражение от подстилающей поверхности. Дано теоретическое обоснование сходимости предложенного метода, подтвержденное численными расчетами. В связи с этим разработан и обоснован алгоритм вычисления матриц Якоби для рассматриваемых итерационных методов.

3. На основе рассмотренных методов разработан комплекс программ, с помощью которого проведены тестовые расчеты для однородного плоского стоя. Сравнение результатов, полученных разными методами показало преим\ щества предложенного ¿метода и целесообразность учета поляризации излучения. Были также проведены расчеты матриц Якоби и их анализ, позволивший обосновать сходимость методов для различных параметров среды.

1. Антюфеев В. С., Михаилов Г. А., Лифшиц Г. Ш., Иванов А. И. Определение аэрозольных индикатрис рассеяния безоблачной атмосферы в спектральных областях 0.55 -г- 2.4 мкм. //Изв. АН СССР. Сер. ФАО. 1980. Т. 16. № 2, С. 146-155. СО АН СССР, 1988.

2. Антюфеев В. С., Назаралиев М. А. Обратные задачи атмосферной оптики. Новосибирск: Изд-во ВЦ СО АН СССР, 1988.

3. Владимиров В. С. Математические задачи односкоростной теории переноса частиц. // Тр.МИАН СССР. 1961. - Т. 61.

4. Гермогенова Т. А., Коновалов Н. В. Спектр характеристического уравнения с учетом поляризации. ИПМ АН СССР, Препринт № 62, М., 1978.

5. Дейрменджан Д. Рассеяние электромагнитного излучения сферическими поли-диеперсными частицами. М., Мир; 1971.

6. Ермаков С. М., Михайлов Г. А. Статистическое моделирование. М., Наука; 1982.

7. Ершов Ю.И., Шихов С.Б. Математические основы теории переноса Т.1. -М.:Энергоатомиздат, 1985.

8. Дэвисон Б. Теория переноса нейтронов. -М.гАтомиздат, 1960.

9. Красносельский M. А., Вайникко Г. M., Забрейко П. П., Рутин,кий Я. Б., Сте-ценко В. Я. Приближенное решение операторных уравнений. М.: Наука, 1969.

10. Кузьмина М. Г. Общие функциональные свойства уравнения переноса поляризованного излучения // Докл. АН СССР, 1978, т. 238, № 2, стр. 314-317.

11. Марчук Г.И. Методы расчета ядерных реакторов. -М.:Госатомиздат, 1961.

12. Марчук Г.И., Лебедев В.И. Численные методы в теории переноса нейтронов. -М.:Атомиздат, 1981.

13. Марчук Г. И., Михайлов Г. А., Назаралиев М. А. и др. Метод Монте-Карло в атмосферной оптике. Новосибирск: Наука, 1976.Engl.transi.: Springer - Verlag, 1980]

14. Михайлов Г. А. Весовые алгоритмы статистического моделирования. Новосибирск: Издательство СО РАН, 2003.

15. Михайлов Г. А. Весовые методы Монте-Карло. Новосибирск: Издательство СО РАН, 2000.

16. Михайлов Г. А. Некоторые задачи теории и приложений методов Монте-Карло. Новосибирск: Препринт ВЦ СО АН СССР, 1978.

17. Михайлов Г. А., Войтишек А. В. Численное статистическое моделирование. М.: Учебно-издательский центр "Академия", 2006.

18. Михайлов Г. А. Оптимизация весовых методов Монте-Карло. М.: Наука, 1987.Engl.transi.: Springer - Verlag, 1992]

19. Михайлов Г. А., Ухинов С. А. Чимаева А. С. Дисперсия стандартной векторной оценки метода Монте-Карло в теории переноса поляризованного излучения //Журн. вычисл. математики и мат. физики. 2006. Т. 46, № 11. С. 2199-2212.

20. Михайлов Г. А., Ухинов С. А., Чимаева А. С. Алгоритмы метода Монте-Карло для восстановления индикатрисы рассеяния с учетом поляризации.// Доклады Академии Наук, 2008, том 423, №2, с. 161-164.

21. Назаралиев М. А. Статистическое моделирование радиационных процессов в атмосфере. Новосибирск, Наука; 1990.

22. Петров В. В. Предельные теоремы для сумм независимых случайных величин. М: Наука, 1987.

23. Смелов В.В. Лекции по теории переноса нейтронов. -М.:Атомиздат, 1978.

24. Сушкевич Т. А. Математические модели переноса излучения. -М.:БИНОМ, 2005.

25. Трачева Н.В., Ухинов С.А., Чимаева A.C. Расчет параметров временной асимптотики выходящего из полубесконечного слоя поляризованного излучения. // Вычислительные технологии. Спец.выпуск; 4: Избранные труды молодых ученых КВМ-2007.

26. Уилкс С. Математическая статистика. М.: Наука, 1967.

27. Ухинов С.А., Юрков Д.И. Оценки методов Монте-Карло для параметрических производных поляризованного излучения //Сиб. журн. вычисл. математики. -2002. Т.5, №1. - С.40-56.

28. Чандрасекар С. Перенос лучистой энергии. М., Наука; 1953.

29. Чимаева А.С. Исследование дисперсии оценок поляризации методом Монте-Карло. // Труды конференции молодых ученых ИВМиМГ СО РАН. Новосибирск, 2004. - С. 242-248.

30. Чимаева А.С. Численное исследование дпсперсии стандартной весовой оценки метода Монте-Карло в теории переноса поляризованного излучения// Труды конференции молодых ученых ИВМиМГ СО РАН. Новосибирск, 2006.

31. Шихов С.Б. Вопросы математической теории реакторов (линейный анализ). -М.:Атомиздат, 1973.

32. Chimaeva A.S., Mikhailov G.A. Study of polarization estimates variance by the Monte Carlo method// Russian Journal of Numerical analysis and Mathematical Modelling. 2006. - Vol.20, №3 - P. 305-317.

33. Chimaeva A.S., Mikhailov G.A., Ukhiniov S. A. Monte Carlo algorithms for reconstructing a scattering phase function with polarization taken into account // Russian Journal of Numerical analysis and Mathematical Modelling. 2009, №5.

34. Kattawar G.W., Plass G.N. Radiation and polarization of multiple scattered light from haze and clouds// Applied Optics, 1968, Vol. 7, No. 8, pp. 1519-1527.

35. Ukhiniov S. A., Yurkov D. I. Monte Carlo method of calculating the derivatives of polarized radiation // Russian Journal of Numerical analysis and Mathematical Modelling. 1998. - V. 13, № 5, pp. 425-444.

36. Ukhinov S.A., Yurkov D.I. Computation of the parametric derivatives of polarized radiation and the solution of inverse atmosphere optic problems// Russ. j. Numer. Anal. Math. Modelling. 2002. Vol. 17. No. 3. P. 283-303.