Содержание к диссертации
Введение
1. Аналитический обзор 12
1.1. Системы оптического дистанционного зондирования 12
1.2. Методы решения уравнения переноса излучения 19
1.3. Трехмерные эффекты переноса излучения 30
Выводы по первой главе 38
2. Дискретное уравнение переноса излучения 40
2.1. Дискретизация уравнения переноса 40
2.2. Многослойные среды 48
2.3. Влияние аппаратно-программных средств на эффективность алгоритмов 58
Выводы по второй главе 65
3. Метод синтетических итераций 66
3.1. Квазидвухпотоковое приближение 66
3.2. Квазидиффузионное приближение 87
3.3. Разорванная облачность 96
Выводы по третьей главе 104
Заключение 106
Источники 108
- Методы решения уравнения переноса излучения
- Трехмерные эффекты переноса излучения
- Многослойные среды
- Квазидиффузионное приближение
Введение к работе
Актуальность работы обусловлена следующими причинами. Среди глобальных задач, отдельно можно выделить такие задачи, решение которых принципиально невозможно без использования методов и средств оптического дистанционного зондирования (ОДЗ). Наиболее актуальной задачей здесь представляется мониторинг малых газовых компонентов атмосферы, вносящих вклад в парниковый эффект. Несмотря на недоказанность глобального потепления (Сорохтин О.Г.) наличие в нашей атмосфере парникового эффекта как такового сомнений не вызывает, а потому изучение его динамики представляется актуальной задачей. Согласно исследованиям (Rayner P.J.), для точной оценки опасности глобального потепления необходим мониторинг стоков и истоков углекислого газа в огромном количестве точек по всей планете, причем точность измерений должна быть очень высокой - не хуже одного процента (Yokota T.). Такие исследования принципиально возможны только при использовании спутниковых систем ОДЗ.
Более того системы, способные проводить измерения с подобной точностью появились относительно недавно (Hamazaki T.) и являются самыми перспективными для ОДЗ, что подтверждает актуальность работы в этом направлении. Все спутниковые измерения являются косвенными, а потому для получения искомой информации необходимо решать обратные задачи ОДЗ, решение которых невозможно без построения прямой модели сигнала спутниковой ОЭС, основанной на процессах переноса излучения в атмосфере. Открывшиеся технические возможности современных измерительных ОЭС кардинально изменили требования к прямым моделям переноса. К повысившейся точности добавилась гиперспектральность – многие системы проводят измерения на десятках тысяч спектральных линий. Этот факт обусловил жесткие требования к ал-
4 горитму и по скорости вычисления – не более одной секунды для одной длины волны. Для достижения необходимой точности в модель требуется включать все известные факторы, существенно влияющие на сигнал. К таким факторам можно отнести анизотропию рассеяния, истинное поглощение газовыми компонентами, отражение подложкой, вертикальную и горизонтальные неоднородности.
Наиболее яркий пример таких неоднородностей – разорванная облачность. В каждый момент времени как минимум половина планеты покрыта облаками (Ardanuy P.E., Wood R., Hamann U.), поэтому измерения почти всегда приходится проводить в разрывы или вблизи облаков. При этом облака оказывают существенное влияние на сигнал и должны быть обязательно учтены при моделировании. Учет в математической модели всех факторов, необходимых для достижения требуемой точности, разумеется, увеличивают время счета. На сегодняшний день остается открытым очень непростой вопрос баланса между скоростью и точностью вычислений. Универсального алгоритма расчета сигнала спутниковой ОЭС ОДЗ, отвечающего требованиям обратных задач в настоящее время нет.
Цель и основные задачи работы
Цель настоящей работы – создание быстрой и точной математической модели сигнала спутниковой оптико-электронной системы с учетом трехмерных эффектов переноса излучения в атмосфере (разорванной облачности). Для достижения поставленной цели в работе были решены следующие задачи:
-
Анализ возможностей ускорения существующих алгоритмов расчета сигналов ОЭС при сохранении высокой точности вычислений;
-
Применение в модели сигнала спутниковой ОЭС принципиально новой идеи метода синтетических итераций, развиваемой в теории переноса нейтронов, в противовес, по-видимому, исчерпавшим себя классическим методам решения уравнения переноса излучения (УПИ);
-
Разработка решения для регулярной части в плоской геометрии с помощью квазидвухпотокового приближения, как эталонного решения для верификации трехмерных методов;
-
Разработка решения регулярной части в квазидиффузионном приближении для плоской геометрии;
-
Формулировка краевой задачи УПИ с граничными условиями в случае зондирования ОЭС через цилиндрическое отверстие в плоскопараллельном бесконечном облаке;
-
Обобщение квазидиффузионного приближения на произвольную геометрию среды.
Основные положения, выносимые на защиту и научная новизна
В настоящей диссертации впервые получены следующие результаты:
-
Предложено правило выбора шага дискретизации УПИ при заданном разрешении на основе теоремы Котельникова;
-
Получена связь решения дискретного УПИ матрично-операторным методом с инвариантным погружением В.А. Амбарцумяна;
-
К описанию переноса в атмосфере применен метод синтетических итераций, который является мощным источником ускорения сходимости;
-
Дано решение плоской задачи в квазидвухпотоковом приближении, которое имеет преимущество по скорости по сравнению с известными алгоритмами более чем в 120 раз;
-
Предложено решение для регулярной части в квазидиффузионном приближении, не уступающее по скорости квазидвухпотоковому приближению и которое легко позволяет произвести обобщение на произвольную геометрию среды;
-
Разработан алгоритм расчета сигнала спутниковой ОЭС при наличии в поле зрения цилиндрического отверстия в бесконечном плоскопараллельном облаке на основе квазидиффузионного приближения.
Практическая значимость диссертационной работы
-
Разработанные модели сигналов спутниковых ОЭС для плоской задачи позволяют проводить верификацию новых методов решения УПИ, как в плоской, так и в некоторых случаях для трехмерной геометрий;
-
Полученное решение в квазидиффузионном приближении дает возможность решать задачи при произвольной геометрии среды, что позволит учесть не только эффекты разорванной облачности, но также и неоднородности подстилающей поверхности;
-
Созданный алгоритм получения сигналов спутниковой ОЭС позволяет проводить обработку и интерпретацию данных пассивных систем ОДЗ;
-
Аналитическая форма полученного решения уравнения переноса излучения допускает построение обратного оператора в аналитической форме;
-
Аналогия между задачами переноса излучения, рассеяния частиц и задач теплопроводности, позволяет использовать все полученные результаты для моделирования процесса переноса частиц и тепла в веществе и интерпретации экспериментов в этих областях;
Достоверность результатов диссертационной работы
Подтверждается математической строгостью всех преобразований, сравнением результатов, полученных по предлагаемому методу, с численными результатами, полученными другими методами (модифицированный метод дискретных ординат, метод Монте-Карло, метод сферических гармоник, метод конечных разностей), а также аналитическими решениями, полученными другими исследователями, сравнением с общепризнанными в ОДЗ алгоритмами DISORT (NASA), SCIATRAN (ESA), Pstar (JAXA).
Соответствие диссертации паспорту специальности
Настоящая работа соответствует паспорту специальности 05.11.07 "Оптические и оптико-электронные приборы и комплексы" по следующим причинам. Из паспорта специальности:
7 «Области исследований:
1. Исследование и разработка новых методов и процессов, которые могут быть положены в основу создания оптических и оптико-электронных приборов, систем и комплексов различного назначения».
Подавляющее большинство оптических измерений являются косвенными, а потому получение искомой информации основано на решении обратных задач. В свою очередь решение обратных задач невозможно без построения прямой модели. В работе предложена математическая модель сигнала спутниковой оптико-электронной системы оптического дистанционного зондирования Земли, основанная на новом методе решения уравнения переноса излучения. Из паспорта специальности:
«2. Разработка, совершенствование и исследование характеристик приборов, систем и комплексов с использованием электромагнитного излучения оптического диапазона волн, предназначенных для решения задач:
- исследования и контроля параметров различных сред и объектов, в том числе при решении технологических, экологических и биологических задач».
Обратные задачи оптического дистанционного зондирования выдвигают к прямым моделям переноса очень высокие требования по скорости и точности вычислений. Для этого требуется включать в математическую модель все факторы, существенно влияющие на сигнал. Впервые полученное в работе решение на основе метода синтетических итераций учитывает все известные факторы, в том числе эффекты разорванной облачности, и при этом отвечает требованиям по скорости вычислений.
Апробация работы
Всего по результатам диссертации опубликовано работ: 6 статей и 16 тезисов докладов. Результаты работы докладывались и обсуждались на 20 конференциях, симпозиумах и школах-семинарах, среди которых: Международная научно-техническая конференция студентов и аспирантов «Радиоэлектроника, электротехника и энергетика» (XIX, февраль 2013; XX февраль 2014; XXI февраль 2015; XXII февраль 2016, Москва, Россия); Научно-техническая конфе-
8 ренция "Молодые светотехники России" (ноябрь 2012; ноябрь 2013; ноябрь 2014, Москва, Россия); Международная научно-техническая конференция «Оптические методы исследования потоков» (XII, июнь 2013; XIII, июнь 2015, Москва, Россия); Международный симпозиум "Атмосферная радиация и динамика" (июнь 2015, Санкт-Петербург, Россия); Оптика атмосферы и океана. Физика атмосферы (июль 2013; июль 2015, Томск, Россия); The 36th International Symposium on Remote Sensing of Environment (май 2015, Берлин, Германия); Современные проблемы физики и технологий. III Межд. мол. науч. школа-конф. (апрель 2014, Москва, Россия); Научн. сессия НИЯУ МИФИ-2014 (январь 2014, Москва, Россия); Междун. конф. "Аэрозоль и оптика атмосферы" (к столетию Г.В. Розенберга) (октябрь 2014, Москва, Россия); The Seventh International Symposium on Radiative Transfer, RAD-13 (июнь, 2013, Кушадасы, Турция); International Radiation Symposium (август 2012, Берлин, Германия); «Современные проблемы оптики естественных вод» (VII, сентябрь 2013; VIII, сентябрь 2015, Санкт-Петербург, Россия).
Структура и объем работы
Работа изложена на 106 страницах, содержит 14 рисунков и 1 таблицу. Список литературы включает в себя 127 наименований на 12 страницах. Общий объем работы – 120 страниц.
Методы решения уравнения переноса излучения
Оптические методы зондирования можно разделить на два больших типа: пассивные (регистрирующие отраженное солнечное излучение) и активные (использующие искусственную подсветку).
Импульсное отраженное атмосферой лазерное излучение несет информацию о распределении концентрации аэрозоля по всей трассе зондирования. При этом пространственное разрешение, определяемое длительностью зондирующего импульса, составляет единицы метров [Козинцев и др., 2002].
Пассивные оптические методы находят широкое применение в спутниковых наблюдениях за газовым составом атмосферы. Хотя в большинстве случаев эти методы дают лишь суммарное содержание газовой компоненты вдоль трассы, они являются единственно возможным инструментом для наблюдения глобальных антропогенных изменений состава атмосферы [Козинцев и др., 2002].
На сегодняшний день существует большое количество спутниковых систем, имеющих на борту аппаратуру для съемки Земли производства США, Китая, Германии, Индии, России и др. [Белов и др., 2014]. Мы же подробнее остановимся на нескольких, представляющихся нам в контексте настоящей работы наиболее интересными.
Среди спутников, основной целью которых является мониторинг малых газовых компонентов атмосферы, необходимо отметить первую в мире систему, предназначенную для измерения концентраций углекислого газа и метана, двух главных парниковых газов из космоса – Greenhouse Gases Observing Satellite (GOSAT), расположенной на спутнике «IBUKI». Программа GOSAT реализуется совместными усилиями Ministry of the Environment (MOE), the National Institute for Environmental Studies (NIES), и Japan Aerospace Exploration Agency (JAXA). Спутник запущен 23 января 2009 года. Анализ измерений GOSAT открывает возможность определения глобальной картины распределения CO2 и CH4, а также, каким образом их стоки и истоки изменяются с течением времени. Эти новые данные должны повысить научное понимание причин глобального потепления [gosat.nies.go.jp].
GOSAT измеряет отраженное от земной поверхности и атмосферы излучение в инфракрасном диапазоне. Концентрация газов в столбе вычисляется на основе этих измерений [gosat.nies.go.jp].
Измерительная система GOSAT носит название Thermal And Near-infrared Sensor for carbon Observation (TANSO) и состоит из двух компонентов: Фурье-спектрометра и картографа аэрозоля и облаков. Три из четырех каналов спектрометра чувствительны к поляризации. Информация с картографа облаков и аэрозоля используется для определения наличия облаков, попавших в поле зрения Фурье-спектрометра. Когда облака попадают в поле зрения спектрометра, производится измерение характеристик облаков и количества аэрозоля. Полученная информация используется для корректировки в измерениях Фурье-спектрометра эффектов, вносимых облаками и аэрозолем [gosat.nies.go.jp].
Важно отметить, что точность измерительной аппаратуры GOSAT очень высока – приблизительно ±1 ppm [Cogan et al., 2012]. В статье [Belikov et al., 2014] было промоделировано и проведено сравнение с данными GOSAT распределение субарктического углекислого газа. Отмечается, что результаты моделирования являются вполне приемлемыми (особенно, если учесть размеры территории), однако, при этом они остаются довольно грубыми. Еще двумя системами, чувствительными к поляризации являются системы французского космического агентства CNES POLDER-1 и POLDER-2 (POLarization and Directionality of the Earths s Reflectances), расположенные на спутниках ADEOS-1 и ADEOS-2 (ADvanced Earth Observation Satellite) соответственно. Оба спутника прекратили работу меньше, чем через год после запуска в результате повреждений солнечных батарей [smsc.cnes.fr/POLDER/]. Основное назначение POLDER – определение оптических и физических свойств аэрозоля с целью его классификации и изучения его изменений; улучшение климатологического описания определенных физических, оптических и радиационных свойств облаков; точное определение влияния аэрозоля и облаков в радиационный баланс Земли smsc.cnes.fr/POLDER/]. POLDER представляет собой широкопольный изображающий радиометр, с помощью которого впервые стало возможно систематически получать глобальные данные о спектральных и поляризационных характеристиках отраженного системой Земля-атмосфера солнечного излучения. Его оригинальные на тот момент возможности наблюдения открыли перспективы к распознаванию излучения, рассеянного в атмосфере, и излучения реально отраженного от поверхности Земли. Данная система производит измерения на 15 спектральных линиях, расположенных от 443 до 910 нм. Две из них являются линиями поглощения газов: 763 нм (О2) и 910 нм (Н2О) [smsc.cnes.fr/POLDER/].
Принципиально инструмент измерения состоит из камеры с двумерным ПЗС-приемником и вращающимся барабаном, на котором расположены спектральные и поляризационные фильтры. Меняя фильтры можно получать информацию об излучении во всем используемом спектральном диапазоне. Относительная точность измерений составляет 2-3% [smsc.cnes.fr/POLDER/].
В рамках второй французской спутниковой программы PARASOL (Polarization and Directionality of the Earth s Reflectances) в 2004 году на орбиту была выведена еще одна система POLDER. Отличием является лишь меньшая приблизительно на 100 км орбита и улучшенный канал связи (16.8 Мбит/с против 883 Кбит/с). Цель миссии PARASOL – определение радиационных и микрофизических свойств облаков и аэрозолей путем измерения поляризации отраженного солнечного света
Трехмерные эффекты переноса излучения
В виду требований, предъявляемых измерительным оборудованием ОДЗ к прямым моделям переноса, а также с учетом современного уровня развития вычислительной техники, в подавляющем большинстве случаев нас интересует численная реализация решения УПИ. Для этого требуется привести его к дискретному виду. Сходу такая процедура не представляется возможной, поскольку краевая задача УПИ всегда содержит разрыв в ГУ [Гермогенова, 1986]. Поэтому для перехода к дискретному УПИ требуется проделать следующее: 1. Аналитически выделить все особенности (разрывы в ГУ); 2. Представить интеграл рассеяния в виде суммы с конечным числом членов.
Для этого обратимся к краевой задаче УПИ для плоского слоя мутной среды. Отметим следующее: несмотря на то, что целью настоящей диссертации является создание модели, учитывающей трехмерные эффекты, плоская задача представляет немалый интерес. Во-первых, для нее существует множество известных решений, с которыми удобно провести сравнение нового. Во-вторых, некоторые случаи трехмерной геометрии при предельных значениях параметров вырождаются в плоскую задачу, что позволяет проводить сравнение и с ними. Перепишем еще раз краевую задачу для плоского слоя мутной среды толщиной т0 при освещении его плоским мононаправленным (ПМ) источником в направлении 10 = U1 - //02,0,//0 : ju dL(T, ) = -L(T,\) +—i JC(1, І Щт, X)dN, dr 4л-J (2.1) = 0. т=т0, /J 0 І L (г,1) =S(\-\0), L(т,\) г=0, ц 0 Для того чтобы выделить особенность, представим решение в виде суммы L(r,l) = L(r,l) + L(r,l), (2.2) где La{r,l) - анизотропная часть (содержит особенность, и вычисляется аналитически), Lr (г,l) - регулярная (гладкая) часть, которая является гладкой функцией, и может быть представлена конечной суммой. Подставив (2.2) в УПИ (2.1), получим новую краевую задачу для гладкой части решения есть функция источников от невязки. Изменения ГУ в (2.3) связано с тем, что приближенное решение не может удовлетворять точным ГУ. Выражение (2.4) нельзя использовать для прямых расчетов, поскольку La(r,l) содержит функцию, которая связана с прямым нерассеянным излучением. Выделим ее:
Важным и не очевидным является вопрос, какую часть тела яркости считать анизотропной? В угловом пространстве это сделать затруднительно, поскольку неясен критерий оценки величины производной, которую можно было бы считать соответствующей анизотропной части решения. Наиболее удачной здесь оказалась идея перехода к угловому спектру, обладающему известным свойством: чем больше анизотропная угловая зависимость функции, тем более пологим является ее спектр. Следовательно, нужно понимать под La(r,l) часть решения соответствующего наиболее медленной изменяющейся части углового спектра. Метод на основе данной идеи получил название малоугловой модификации метода сферических гармоник (МСГ), был разработан и подробнейшим образом описан в работах [Будак и Савенков, 1980; Будак и Савенков, 1982; Будак и др., 1983; Будак и Сармин, 1990], поэтому мы не станем останавливаться подробно на выкладках, а приведем только окончательное выражение: где Qj(M) = J-——Р/С") - перенормированные полиномы Лежандра, Р,и(//) присоединенные полиномы Лежандра, Zk(r) = e_u?ir/M)-е т , dk=\- Ахк.
Однако важно отметить, что данный метод является наиболее общим среди всех малоугловых методов [Будак и Федосов, 1985], а самое главное - он позволяет наиболее точно выделить анизотропную часть решения, тем самым делая регулярную практически изотропной по углу функцией, что важно при численной реализации [Kokhanovsky et al, 2010; Sokoletsky et al, 2014]. Исходя из этих соображений, анизотропную часть решения в настоящей диссертации мы будем представлять именно в МСГ.
Вследствие выделения анизотропной части решения, функция Lr(r,l) будет являться гладкой по углу, а потому ее представление возможно по конечному базису. Например, если для дискретизации УПИ используется МДО, то такое представление будет иметь вид м
Это позволяет заменить интеграл рассеяния конечной суммой, а краевая задача (2.3) преобразуется в краевую задачу для матричного неоднородного линейного дифференциального уравнения первого порядка с постоянными коэффициентами [Budak and Korkin, 2008; Budak et al., 2011; Budak et al, 2014]: есть решение соответствующего однородного уравнения. Она описывает связь распределений яркости в двух точках t и г среды и носит название «пропагатор» [Budak et al., 2014]. Проблемы реализации решения (2.10) связаны с наличием в пропагаторе (2.11) как отрицательных, так и положительных экспонент, что физически соответствует распространению потоков сверху-вниз и снизу-вверх. Это приводит к быстрому ухудшению обусловленности матрицы и делает невозможным расчеты для световых полей с оптическими толщами г 1. Для устранения этого эффекта используется масштабное преобразование [Karp et al, 1980]:
Многослойные среды
Одной из задач прямого моделирования, которая требует обязательного обращения к трехмерной теории переноса излучения, является учет в алгоритме эффектов разорванной облачности. Квазидвухпотоковое приближение дает хорошие результаты в случае плоской задачи. Однако в силу специфики двухпотокового метода, зависящего от горизонтальной симметрии среды, в случае произвольной геометрии удобнее обратиться к спектральному представлению яркости не только в решении для анизотропной части, но и для регулярной. С точки зрения идеи и простоты аналогом двухпотокового приближения в частотной области является диффузионное приближение. Этот метод представляет собой метод сферических гармоник при ограничении минимальным количеством членов – одним (P1-приближение). Выбор данного метода (как и двухпотокового) обусловлен теми же причинами: - один из простейших, а значит и быстрейших с точки зрения реализации; - после выделения анизотропной части решения по МСГ, регулярная часть близка к изотропной; - неточности данного этапа будут исправлены последующей итерацией.
Впервые идея решения УПИ методом синтетических итераций с определением решения на первом шаге в диффузионном приближении была сформулирована в работе [Гольдин, 1964] и получила название квазидиффузионного приближения. Этот метод был развит в [Аристова и Гольдин, 1998] на случай сильной анизотропии тела яркости, для чего в решении выделялось прямое нерассеянное излучение. Анализ решения УПИ для сред с сильно анизотропным рассеянием, проведенный в работах [Будак и др., 2012; Budak et al., 2014], показывает, что существенно более эффективно выделять анизотропию на основе МСГ. Поэтому для анизотропной части мы снова используем МСГ, а регулярную решаем на основе синтетических итераций с диффузионным приближением в качестве первого шага. Снова для начала обратимся к плоской задаче. Уравнение для гладкой части в этом случае будет иметь вид:
Далее, умножим (3.93) на 1 и снова проинтегрируем по телесному углу. Учтем соотношения: 1. фї(А-ї ії = тА, 2. фї(А-ї)(В-ї й = 0, что приведет к уравнению ІУЕ 0(r) + (l-A/7)(r) = s(r). (3.95) Возьмем дивергенцию от выражения (3.95) -АЕ0(г) + є(1 - A/7)V(r) = Vs(r), (3.96) где для простоты допущено, что є Ф є(г) - среда однородная. Выразим V(r) из (3.94) V(r) = s0(r)-s(l-A)E0(r), (3.97) и подставим его в (3.96) -Щ(г) - s\\ - Л/7)(1 - Л)Я0(г) = -є{\ - Ap)s0(r) + Vs(r). (3.98) Разделим обе части уравнения (3.98) на а2 = є2(\ - Ар)(\ - А), что приведет к уравнению диффузии: AEJr)-EJr) = ! sn(r) + ! Vs(r). (3.99) а2 Ъе{\ -Л) Зє2 (\-АД)(\-А) Рассмотрим функцию источников от невязки МСГ в предположении однородного бесконечного слоя. В этом случае функция источников имеет вид [Budak et. al, 2015]: где Zk(z) = e 4rM)-е Г/Ч , dk=\-Axk. Для простоты анализа алгоритма ограничимся аналитически простым случаем нормального облучения среды где частное неоднородное решение Ep (), исходя из вида функции источников (3.107) должно иметь вид
Анизотропная часть в квазидиффузионном и квазидвухпотоковом методах мы вычисляем с помощью одного и того же инструмента - МСГ, которая практически полностью определяет излучение в нижнюю полусферу. Поэтому с точки зрения сравнения в этом случае представляет интерес только обратная полусфера, что соответствует сигналу спутниковой оптико-электронной системы. На рис. 10-11 показаны сравнения реализаций сигнала от угла визирования квазидвухпотоковым и квазидиффузионным методом. Рисунок 10 - Сравнение квазидиффузионного и квазидвухпотокового методов при 0=0о, =0.9, g=0.9, 0=2.
В качестве случая разорванной облачности рассмотрим цилиндрическое отверстие в плоскопараллельном облаке (разрыв), облучаемом плоским мононаправленным источником под углом 0 (геометрия представлена на рис. 12). Очевидно, что на определенном расстоянии от облака (в случае, когда мы имеем дело с его краем) и при определенном размере отверстия (в случае разорванной облачности), влияние этих эффектов становится несущественным. Поэтому важным вопросом представляется оценка значений этих величин. Поэтому задача цилиндрического отверстия представляется актуальной.
Поскольку анизотропия является локальным свойством уравнения, поэтому для ее определения мы воспользуемся предположением о независимом слое мутной среды с оптической толщей, равной оптической толщи для прямых лучей. Исходя из этого, выражение для анизотропной части решения будет иметь такой же вид, как и в случае плоской задачи, однако, оптическая толща (r) будет зависеть от пространственных координат. Для определения выражения оптической толщи в случае с цилиндрическим отверстием необходимо найти пересечение луча с цилиндром. Уравнение луча вдоль направления падения мононаправленного источника имеет вид
Квазидиффузионное приближение
Верификацию алгоритма для отверстия в облаке в квазидвухпотоковом приближении мы проводили, реализовав ту же задачу методом Монте-Карло. Метод Монте-Карло [Михайлов, 1974; Ермаков, 1975; Марчук, 1976; Ермаков и Михайлов, 1982] заключается в математическом моделировании исследуемого процесса, проведении достаточно большого числа численных экспериментов и статистической обработки их результатов. При исследовании процессов переноса излучения моделируются траектории фотонов. Процесс переноса фотонов в среде можно описать [Золотухин и Ермаков, 1972] интегральным уравнением переноса /(JC)= \f(x )k(x\x)dx + \y(x), (1.131) (X) или в операторной форме где(х) - плотность столкновений в фазовом пространстве X: (г,1), х, х є X; ф) -плотность столкновений в источнике; к(х ,х) - плотность вероятности того, что частица, испытавшая рассеяние в точке х , следующее столкновение испытает в где x0 - начальное положение частицы; {xn} - однородная марковская цепь столкновений частицы; N - номер последнего столкновения перед поглощением; M - означает усреднение по ансамблю траекторий.
Поскольку для оценки точности решения требуется произвести расчет пространственной облученности в точке или малом фазовом объеме, то используется метод локальной оценки. Метод локальной оценки основывается на выражении (1.131), которое можно также рассматривать как искомый функционал, аналогично (1.133). Тогда в соответствии с (1.134) математическое ожидание величины равно облученности в точке r от потока рассеянных средой фотонов. Следовательно, практическая реализация алгоритма сводится к моделированию последовательности точек столкновения частиц в среде с последующим усреднением по всем возможным траекториям. Большим достоинством метода локальной оценки является возможность получения информации о сигнале для нескольких приемников одновременно с помощью одной и той же цепи Маркова.
Применение различных алгоритмов метода Монте-Карло изложено в [Михайлов, 1974; Ермаков, 1975; Марчук, 1976; Ермаков и Михайлов, 1982]. Основываясь на методе локальных оценок и рекомендациях указанных статей, были разработаны алгоритмы, реализованные в системе MATALB.
Цилиндрическое отверстие учитывалось при вычислении пересечения: если луч пересекал поверхность цилиндрическую границу отверстия, то он пролетал далее до второй его границы, где уже проходил остаток разыгранного при рассеянии пробега. Координата точек пересечения определялась из системы уравнений для луча и цилиндра в проекции на границу слоя облака, что уменьшает размерность задачи.
Однако в случае очень острых индикатрис метод локальных оценок имеет медленную сходимость или возможны ошибки, если выборка недостаточно велика. С целью ускорения сходимости использовался тот факт, что система ПМ-источника для плоского облака с цилиндрическим отверстием имеет цилиндрическую симметрию, т.е. все акты рассеяния не зависят от поворота источника вокруг своей оси (азимутальный угол ф). Поэтому в одном испытании разыгрывалось К траекторий, соответствующих К значениям выбранных ф. Все К траекторий разыгрывались по одной и той же последовательности случайных чисел. Обрыв траектории происходил при вылете фотона через сферу с центром в точке источника, а радиус ее выбирался экспериментально. Расчет останавливался при достижении дисперсии менее 5%.
Индикатриса рассеяния для простоты бралась в виде Heneye-Greenstein, но, поскольку было проведено несколько серий испытании для различных g, то можно говорить о точности решения для любой реальной индикатрисы. Формула для разыгрывания косинуса угла рассеяния в случае индикатрисы Henyey-Greenstein имеет вид 2 Ll = 2g где а - равномерно распределенная в [0,1] случайная величина. Результаты сравнения для различных параметров представлены на рис. 13-14. Рисунок 13 – сравнение алгоритмов задачи цилиндрического отверстия в облаке при =0.999, g=0.99, 0=5 (– МСГ, о Монте-Карло, - квазидиффузионное приближение).
Поскольку задача в данном случае обладает осевой симметрией, распределения построены от центра отверстия. Занижение Монте-Карло в правой части связаны с тем, что для численных расчетов облако принципиально закладывается конечных размеров. В большинстве остальных точек расхождение с квазидиффузионным приближением составляет менее процента, при этом время на расчет Монте-Карло - около часа, квазидифузионное – менее секунды.