Решение обратных задач переноса излучения методом Монте Карло тема автореферата и диссертации по математике, 01.01.07 ВАК РФ

^ л

А ^

^ российская академия наук

Сибирское отделение Вычислительный пентр

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

Антюфеев Виктор Степанович

решение обратных чадач переноса 1плучения .методом монте карло

01.01.0i вычислительная математика

А В т и ]) е ф е р а т дне ( ерташш на соискание ученой степени доктора фтико-матема гическпх наук

Новосибирск 1997

PauOJH ВЬШО.ШГЛа В В LI 4 J f С. Ш J i". !bllO\J i[cH 1 p<4 ( HOJipCKOJ О (НДГ-КЧШЯ PA'í

Официальные oiiiiuiKMi i !,i:

чл'^-к^ррогиондгнт РАН,

11 j ), )ij>cf'v'i 'j ) Романов Н.Г..

jiuk i(i]> фи шко-ма i(\\ia Wi'uk >. ix наук.

профессор 1 рШОрЬОВ 10. Ii.,

док юр фпшко-ма юмлшчсскпх наук.

проф(ч сор Сушксвич i .Л.

Ислупыи (»pi a muai in я: Ип< iiiivi i copc i п ческой n i qui клал« ¡oii ме-

ханики ( 'jlGupChOl о оикмишя Рос< ийг.м-п Хкадсмпи на\'к.

if in ¡и I а < »я ion и 'я " ок I я оря ! !)')? iniri и . _ v ort па tac r/i-,t . . . псьпа.in

5111 )• »Hri h 11 ' » i о сове I a Jl OD'J. 1 O.I) I m i ми in И' no < cp ! <i i,»:f, на ' от kali иг \ чопоп с i (чнчш ;\чк íopa Urin* ti ¡"i i ím'íiíc. iii i с ii,iui\i цеп j pe í"() рдн it.->()(i')(l. IIokix iionp( k-!H). np Ak. .'biHpcH ! Ы'ва. í¡)

( MI' < » ] > I <•! 1111 « ' П МО/Kilo <i illa K< >\1 И IM Я В 4111 ¡I.II.IKAI ULI С оиб I ПО I Ok (I Uli, CO l':\ii i'i.)(Hi'Hj. líohoí HOllpí k 'И1. |||>. \|x. -!{«(H¡H4I 1 !,CH<t. (t)

\ В J itj/f •/}>' p.» J p,i Их . I,:f

^. ti-ll I 1,

lililí/])« I1)1), i lijin

• |<I1!.:1! < •(Ц-1

•i"I"'

¡ii í; ¡ p' »ни m luí i > t ¡jhc i.

lo.h.kv шиши

В диссертации рассматриваются обратные задачи переноса излучения II связанные с нтш вопросы расчета характеристик поля рассеянного излучения методом Монте Карло.

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

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

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

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

1ILJX ЫОффИШНЧТГОВ.

15 настоящей работы в основном рассматриваются итерационные методы, позволяющие в предел«1 получить точные значения ИСКОМЫХ величин. Л ИШЬ В Задаче компьютерной томографии значения требуемых величин удается получить за счс1 <-цсштльного расположения приемников и источников излучения.

Решение обратных задач итерационными .методами приводи! к необходимости решать "прямые" задачи, тоепь задачи расчета тюля рассеянного излучения. На сегодняшний день лучший метод для таких расчетов это .метод Моше Карло. Он позволяет включать и модель среды практически псе реальные параметры, вычш лят ь производные по параметрам среды, корреляционные характеристики ноля яркостп.

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

Н paooie мало места уделено доказательству сходимости рассма I рпнаемых итерационных .методов, так как для практически важных случаев условия сходимости невозможно проверить. Теоретические же границы сходимости очень далеки О) Ирак Iнчееких.

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

для пло( кого слоя разработан метод час тичного "математическою ожидания", позволяющий сократить трудоемкость расчета примерно в З-ö раз. Для плоских слоев большой оптической толщины разработана модификация чтото метода, позволяющая справиться с ен> '"неустойчивостью";

на основе метода частичной» "математического ожидания" получена искусственная (нефизическая) плотность для длины пробега, которую можно использовать для отнмпза-ции расчетов в задачах со сферической геометрией:

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

ногти поля рассеянного солнечного излучения:

простое преобразование позволило частично свести расчет переноса излучения в сферической атмосфере к менее трудоемкому расчету для плоского слоя:

предложено преобразование нестандартного ("обобщенного"') интегральной) транспортного уравнения, которое позволило использовать замену сильно вытянутой индикатрисы рассеяния ее " псевдо-транспортным" приближением. Для задачи расчета рассеянного светового потока в разорванной облачности получены выражения для соответствующей плотности длины пробега фотона. Использование этого метода в практических расчетах позволило уменьшить их трудоемкость примерно в 30-40 раз:

для стохастической среды на основе приближенной линеаризации модели получены формулы для расчета" корреляционных характеристик поля яркости. Получены формулы для расчета методом Монте Карло прнзводных по коэффициентам спектрального разложения поля среды. Доказана конечность дисперсии статистической оценки производных;

разработан метод определения индикатрисы рассеяния по наблюдениям яркости рассеянного солнечного излучения в точке на поверхности Земли в направлениях альмукантарата Солнца. Теоретически- исследована сходимость этого метода. Проведены практические расчеты с использованием экспериментальных данных;

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

рассматривается применение1 метода Монте Карло к задачам переноса излучения в растительном покрове. Проведены расчеты поля рассеянного излучения, давшие результаты. близкие к экспериментальным:

предложен метод восстанавления основных параметров растительного покрова. Предложена новая схема измерения потока рассеянного излучения и метод восстановления па-

раметров плотности распределения нормали к поверхности листьев:

рассматривается задача томографци с рассеяние;,!. Предложена схема расположения источников приемников и метод восстановления распределения коэффициентов ослабления и рассеяния внутри исследуемого тела. Разработанный метод проверен расчетами:

предложен метод решения интегральных уравнений второго рода ( <>'стохастическим ядром.

Все основные результаты диссертации являются новыми.

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

' Предложенные модификации метода Мол тс Карло значительно повышают эффективность расчета характеристик поля рассеянного излучения.

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

Были проведены расчеты для Астрофизического института (Алма-Ата). Института прикладной геофизики (Москва) н Института астрофизики и физики атмосферы (Тарту).

Апробация. Основные результаты диссертации докладывались и обсуждались на международных и всесоюзных конференциях и семинарах в Алма-Ате (АФИ). Баку (ИКИ), Ереване (ИПМ). Ленинграде (ГОИ). Москве (ИНГ. ИНМ. И ФА). Новосибирске (ВН. ИМ). Тарту (ИАФА). Томске (ИОА). Коршевеле (Франция, семинар по диез, зондированию). Хельсинки (Финляндия. ЮАИБЯ'ОЗ). Индианаполисе (С'1ЛА. рабочее совещание по томографии).

Публикации. Основные результаты опубликованы в рапс пах [1-29] (в их числе 1 монография).

Структура и объем. Работа состоит из введения, семи глав, приложения, заключения и списка литературы. Объем работы 23-"> страницы.

а

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

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

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

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

В 1.2 предложена экономичная модификация метода Монте Карло для расчета переноса излучения в плоском слое.

Рассматривается интегральное уравнение переноса для плоского горизонтально-однородного с лоя

Его решение представляется в виде ряда Неймана

ос

Р = £ Гв.

»=0

Запишем подробнее п-п член этого ряда: 1пв(г. и) - /... / <1шо ... («„)7?. • / <1

<7(1-:

$2 52

ы

ехр

М

1 - с0.им) .. ./(~„_I - :„_2. )/(; - :„_1 .ш).

где

£ехр 0< = <1-.

0. : $ [0.1]

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

*£.4;'рхр(-а/ = ),

а но вычислять его метолом Монте Ка]>ло. Алгоритм расчета изменяется: теперь достаточно моделировать лишь марковскую цепь -направлений пробега частицы и при очередном "столкновении" заносить в статистическую оценку вклад, который теперь зависит не только от очередной точки столкновения (как обычно), а и от всех предыдущих (немарковская статистическая оценка).

Этот метод позволил повысить эффективность расчета для небольших значений оптической толщины слоя в несколько раз.

Однако, расчеты показали, что эффективность метода падает <• увеличением оптической толщины т, а при т > 10 программа перестает работать из-за ошибок переполнения.

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

охр .г — ехр у ------:_ ^ ОХр г

~ !)

Если .г ~ у/, то эю соотношение позволяет избавиться от малой величины в знаненател< дроби и, соответственно, от ошибок округления компьютера. Однако, при этом в линейной комбинации экспонент появляются дополнительные "базисные" функции вида ; ехр(—с-г/г,-) с соответствующими коэффициентами, которые приходится пересчитывать по новым формулам.

В 1.4 точное выражение для плотности распределения частиц по высоте плоского слоя при фиксированной истории направлений, полученное в 1.2. и» пользуется в качестве приближенной функции ценности для модификации моделирования 1).г1пил пробе! а » расчетах со сферической моделью атмосферы.

В 1.5 рассматривается другое приближение функции ценности, которое приводит к модификации моделирования ма-правлетшщю&ет'л частицы при моделировании траектории. Это приближение основано на приближенной замене экспоненты в уравнении переноса на некоторую величину ц < 1. В результате такой замены уравнение переноса превращается в простое интегральное уравнение для функции на сфере

Ь(ш) = ц [<!ш' + 4,.(а>).

о 2?г

Затем доказывается, что вследствие симметрии поворота эта функция на самом деле зависит лишь от одного переменного, н пз уравнения можно найти ее коэффициенты в разложении в ряд по полиномам Лежандра. Найденную таким образом функцию можно использовать в расчетах по сопряженной схеме (см. 1.1) как приближенную функцию ценности.

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

и 1 - с/2

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

В главе 2 рассматривается уравнение переноса радиации в стохастической среде. Этот процесс может быть описан детерминированным интегральным уравнением второго рода. Расчет характеристик светового поля методом Монте Карло приводит к большой статистической погрешности, если индикатриса рассеяния у(//) в этом транспортном уравнении имеет сильно пикообразную форму.

В 2.3 рассматривается представление пикообразной индикатрисы в виде линейной комбинации ее "пикообразной" и "регулярной" составляющей

я - Pfh + чЬ-

затем делается приближенная замена д —» д1г исходной индикатрисы на на сумму регулярной и 8 - функции:

ц,,. = ]>8j -f qfj

Функция назовается д/г"пс('вдо-трннспортным приближением" для индикатрисы рассеяния д.

В '2.4 показано, как в случае стандартного транспортного уравнения можно избавиться от больших ошибок расчета, выполнив простое алгебраическое преобразование для транспортного уравнения в интегро-дифференциальной форме. В результате замены г/ —> ;/tr исходное уравнение

ш ■ V$(/-.a>) + а(г)Ф[г,ш) =

= а[г ) j Ф(г.и>)д(ш' ■ w )</«' + Фо(г. и)

и

п]ш.\гет вид

о; • \ГФ(г.и»') + 7(т(,г)Ф(г,ы) =

= ya(r) j Ф(г. ш)д(ш' • u>)<lu>' -f Фо(г.ш).

h

и мы получим уравнение переноса для среды с новым сечением рассеяния о — f/Д" и новой индикатрисой рассеяния у. Так как функция Ф(г,о/) в уравнении при этом осталась прежней. можно утверждать, что величина потока не изменит-( я (точнее говоря, изменится незначительно), если моделировать траектории частиц в < роде с параметрами <7, <у вместо исходных параметров т.д. При этом новая среда будет оптически менее плотной, чем исходная (при прежних геометрических размерах). Поэтому среднее количество моделируемых столкновении на одну траекторию частицы сократится. Кроме того, уровень вариаций значений новой индикатрисы рнс< еяння д значительно ниже, чем у исходной индикатрисы У-

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

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

В 2.3.2 детально рассматривается моделирование процесса переноса излучения (в лучевом приближении). Точки столкновения, где частица меняет направление пробега, можно условно разделить на две группы.

В первую группу отнесем все точки столкновения, где частица изменила направление полета (то есть, направление полета ш моделировалось согласно индикатрисе рассеяния </). Назовем эти точки точками ''истинного" рассеяния.

Во вторую группу отнесем все остальные пространственные точки рассеяния, то есть те точки, где частица не изменила своего направления полета. Такие рассеяния принято называть дельта-рассеяниями. Эти точки мы будем называть также точками промежуточного рассеяния.

Теперь, моделируя траектории, можно не обращать внимания на точки промежуточного рассеяния. Надо лишь научиться моделировать случайную длину пробега от одной истинной точки столкновения до другой.

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

Теперь становится понятным, что простое алгебраическое преобразование обобщенного транспортного уравнения не даст тот же результат, что и в стандартном случае. Необходимо подвергнуть интегральное уравнение преобразованию Лапласа, совершить алгебраическое преобразование в области образа и вернуться обратно. Это сделано в 2.5.3.

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

В главе 3 изложен способ вычисления моментов корреляции

B]2 = {(h-(mi2-m).

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

В 3.3 приведена постановка задачи. Считаем, что поло а(/ ) случайного коэффициента ослабления представим« в виде суммы:

Ф) = <Г(г)+<т{г)

где

if (г) = <а(г)>.

lain а(г) = (тп > 0.

н

а случайное поле а (г) допускает конечное ])азложенне:

к

"(>') = Е

¿-и

где (,- - случайные незавпешшые величины с математическим ожиданием. равным 0 и дисперсией ], коэффициенты />, > 0. y-i(r) неслучайные базисные функции.

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

Ып)« Ik(v) +- Е Df4i = Ш + Е D^p.t,,

j=II 1=11

Отсюда для корреляционного момента Ип получаем

Вп « (Е А'М^'Д.,) = hrfDlnl

1=(/

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

В 3.'2 получены статистические оценки для расчета производных Пспя яркости по коэффициентам разложения случайном! иоля о (г) методом Монте.Карло.

В .')..'! доказана конечность дисперсии опенки этих производных.

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

В 4.1 изложена физическая постановка задачи. Здесь особо е внимание уделено необходимости'учета -эффекта многократных рассеяний при использовании данных измерений. Рассматривается кусочно-линейная аппроксимация индикатрисы. Таким образом, задача о восстановлении функции .'/(/') сводится к задаче о восстановлении ее значений (¡и в узлах аппроксимации.

В 4.2 изложены два варианта метода последовательных приближений для решения этой обратной задачи. Оба метода основаны на выделении "главной части" - однократно рассеянного потока - из полной измеренной интенсивности. Для этого пользуемся следующим наблюдением. Для выбранных направлений наблюдения значения интенсивности однократно рассеянного .Д. излучения линейно зависят от соответствующих значений индикатрисы рассеяния:

•4 = С • <7ь

причем для выбранных направлений наблюдения коэффициент пропорции С известен и не зависит от номера направления к:

ехр(-г) С = ———с о.ч0с,,

¿7Г

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

п - к-« -К - Л = сил - <)Й.

где Ц -- измеренные значения интенсивности, соответствующие искомой индикатрисе д*, а /<• - расчитйнные методом Монте Карло значения интенсивности, соответствующие некоторой индикатрисе д. которая предполагается близкой к

искомой. Отсюда получаем итерационный процесс: на г-ом шаге итерации решается тривиальная (диагональная) система линейных уравнении

СЛц. =/; - У/*', т = 1.....п,

где

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

В 4.4 рассматриваются вопросы погрешности при использования плоской модели атмосферы и расчета интенсивности метолом Монте Карло.

В 4.-3 приведены результаты расчетов н их обсуждение.

Раздел 4.0 повящен исследованию сходимости "аддитивного" итерационного процесса. Доказана следующая

Теорема Пусть = {||<?|| < 7?}, I! > 0, и пусть все координаты векторов <1* + Л, й € ограничены снизу некоторым положительным числом </. Тогда при достаточно малых значениях г равномерно на множестве {1ц выполняется оценка

||6'г(;/) -6'т(//)|| < #/||, ч< 1

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

И главе 5 расматрив'аютея вопросы регуляризации решения двух обратных задач атмосферной оптики:

1) Восстановление высотного хода коэффициента рассеяния по наблюдениям яркости дневного горизонта Земли.

2) Восстановление индикатрисы углового рассеяния света по наблюдениям яркости в альмукантарате солнца.

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

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

В 5.2 известное значение оптической толщины атмосферы

п

Т =

¡=1

используется для регуляризации при восстановлении вектора (Ti,... ,<т„. Исходный итерационный алгоритм заключается в последовательном решении систем уравнений

DU — I* - I

с матрицей D. Регуляризация сводит на каждом шаге задачу к отысканию4 условного минимума функционала Лагранжа

min |DS - (Т* - 1)\'2

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

А —

!>, Ь->

л„

К О

В .">.3 такой же прием ис пользован для регуляризации при восстановлении индикатрисы методом, описанным в главе 4.

В 5.4, 5.5 рассматривается регуляризация тгоп же задачи на основе статистической информации, точнее при известной корреляционной матрице неизвестного случайного вектора (Н...../Ун-

В 5.4 искомый вектор раскладывается в сумму /У = V + !)■

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

(у,/<;) = (), /=)/-/,' + !.....».

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

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

Конструктивно алгоритм решения выглядит так. Вместо исходной-системы линейных уравнений

Б!) = Ъ

с матрицей В следует решать линейную систему вида

(БЧГО + С-1)«* = ОЧТ'Ь

с регулярпзовашюй матрицей И*\УО + С-1, где И' - диагональная матрица с диагональными элементами

1

"'» = 72 '

^ а

а ¿.",7 - среднестатистическая ошибка компоненты вектора Ь правой части.

В главе 6 рассматривается применение метода Монте Карло к решению прямых и обратных задач переноса излучения в растительном покрове1.

В С>.1 дано краткое" описание существующих моделей растительного покрова.

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

01

-/I— + С(и>)Г(т,ш) = / (1ш Г(о»' —> ш)Цт,ш') от ■>

ш

и интегральной

п

.7(т, ш) = Л к((т', ш') (г, «)).7(т', ш')(1т>(1ы' + ^(т, и») н о

форме.

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

В 6.4 описана специальная поправка к коэффициенту осла' бдения, позволяющая учитывать эффект обратного блеска" при переносе излучения в растительном покрове. Она обусловлена тем, что реальный растительный покров представляет собой дискретную, а не сплошную среду. Луч, упавший снаружи на поверхность листа внутри слоя, может после рассеяния на поверхности листа на малый угол выйти обратно через ю же свободное пространство между листьями. Этого не учитывает исходная мутная модель, поэтому приходится вводить поправочный множитель громоздкого вида (который сильно затрудняет решение обратных задач).

В 6.5 рассматривается постановка обратной задачи об определении параметров растительного покрова, итерационный метод Ньютона для ее решения. Для этого метода необходимо знать производные интенсивности по параметрам среды.

В 6.0 оппсанО-вычпсленне этих производных методом Монте Карло.

Описанный метод не дал удовлетворительных результатов при восстановлении параметров плотности распределения нормалей листьев. В связи с этим был разработан новый метод, который изложен в 6.7.

Этот метод, так же как метод, расмотренный в главе 4. основан на выделении величины 1\ однократно рассеянного излучения. Однако, для произвольно выбранных вариантов наблюдения очень трудно получить затем искомые параметры из-за сложного выражения для величины 1\. Здесь рассмотрены специальные варианты наблюдения, когда выражение для интенсивности однократно расс еянного излучения сильно упрощается за счет эффекта оПр;'пкого блеска. Приемник располагается так, что направление наблюдения совпадает с направлением прямого солнечного излучения. При этом в его апертуру попадает преимущес твенно однократно рассеянное излучение. Для произвольно выбранного направления наблюдения выражение для 1\ имеет сложный вид:

i/< • i

х I(1техр г) (>х1'

Особенно громоздким будет выражение для 0к{t.u>l.l,L0*). Однако, для рассматриваемого варианта наблюдения эта

формула сильно упростится:

___________

[ <1ш,, !,(ш1.)Ци»,..ш*)\ И И

В 6.7 показано, как двумерные интегралы в первом сомножителе4 сводятся к одномерным, а затем вычисляются простым разложением в ряд. Из полученного выражения для 1\ искомые параметры нетрудно найти, решая т ерапиями получающееся нелинейное уравнение.

В 6.8 приведены результаты решения обратных задач для растительне^го пежрова.

В главе 7 рассматриваются задачи томе>графии в среде с сильным рассеянием.

В 7.1 описано отличие предлагаемого подхода е>т "стандартной томографии", где" эффект рассейния считается пре-

»■г. , „

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

В 7.2 приводится краткое математическое обоснование метода.

В 7.3 дается детальное изложение метода. Метод состоит в определении коэффициентов уравнения переноса в следующем порядке:

совокупный коэффициент ослабления <т(.г); коэффициент рассеяния аs (х): индикатриса рассеяния <}(■!•, oj ■ коэффициент поглощения 7a(jc).

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

В 7.4 представлены некоторые численные результаты, по-казыгакицие вычислительные применения метода и искажения восстановления образа, вызванные влиянием рассеяния при применении стандартной процедуры томографии. Объект для исследования - сфера радиуса 1, разделенная на две области: концентрическая сфера радиуса 0.5 и остальная ее часть. Первая область заполнена чисто рассеиваю-

щнм (непоглощающим) веществом, а вторая чисто поглощающим. Коэффициенты рассеяния и поглощения постоянны внутри соответствующих областей и одинаковы по величине: (т, = <та = 10. Таким образом величина общего коэффициента поглощения постоянна во всей большой сфере: а = 10. Благодаря относительно большой величине коэффициента рассеяния вклад рассеянного потока излучения составил до 40%. Восстанавливали образ центрального сечения сферы. Данные для "наблюдения" были получены расчетами методом Монте Карло.

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

В приложении рассматривается задача о решении интегральных уравнений второго рода / = А/ + у? со стохастическим ядром, для которых нельзя напрямую применять стандартную процедуру метода Монте Карло. Помогает следующее замечание. Л" оператора К есть тривиальная собственная функция А'о(.г) = 1. Теперь пространство раскладывается в прямую сумму двух подпространств, инвариантных относительно интегрального оператора уравнения: одномерное, натянутое на А'о и ортогональное дополнение к соответствующей собственной функции сопряженного оператора К*. В одном из этих подпространств уравнение имеет тривиальное решение, а для второго подространства соответствующий ряд Неймана сходится. Получена статистическая оценка, позволяющая применять для решения исходного уравнения метод Монте Карло.

Основные результаты диссертации:

1) Получены экономичные модификации метода Монте Карло для решения прямых задач переноса излучения, вычисления корреляционных характеристик поля яркости и производных по коэффициентам спектрального разложения стохастического поля среды.

2) Разработан метод определения индикатрисы рассеяния по

наблюдениям яркости рассеянного солнечного излучения в точке на поверхности Земли. Теоретически исследована сходимость этого метода. Проведены практические расчеты с использованием экспериментальных данных.

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

4) Изучается применение метода Монте Карло к решению прямых и обратных задач переноса излучения в растительном покрове.

5) Предложен метод решения задачи томографии сильно рассеивающего тела.

Публикации. Основные результаты диссертации опубликованы в работах:

1. Антюфеев B.C., Михайлов Г.А. Восстановление высотного хода коэффициента рассеяния и вычисление корреляционных характеристик яркости методом Монте Карло. // Изв. АН СССР, ФАО, 197G, т.ХП, №5

2. Антюфеев B.C. Решение некоторых стохастических задач атмосферной оптики методом Монте Карло // Изв. АН СССР, ФАО, T.XIII, N¡5, 1977 .

3. Антюфеев B.C., Авасте O.A., Вайникко Г.М., Вилл-манн У.И.^Гречко Г.М., Губарев A.A.. Климук П.И., Князи-хин Ю.В. Восстановление высотного профиля коэффициента ослабления аэрозоля по оптическим измерениям из космоса // Сб." Исследования атмоеферно-оптических явлений с борта орбитальной станции Салют-4", Тарту, 1979

4. Антюфеев B.C. Метод математических ожиданий для решения задач переноса излучения с плоской геометрией. // ЖВМ и МФ, 1979, т.19, №6

5. Антюфеев B.C., Иванов А.И., Лившиц Г.С., Михайлов Г.А. Определение аэрозольных индикатрис рассеяния безоблачной атмосферы в спектральной области 0,55-2,4 мкм. // Изв. АН СССР. Сер. ФАО, 1980, т.16, №2

G. Антюфеев B.C. О сходимости одного итерационного Процесса. // Сб. Численные методы и статистическое мо-

делированпе в теории переноса, Новосибирск, 1980

7. Антюфесв B.C.. Авасте С).А., Вайнпкко Г.М., Всйеманн У.К., Виллманн У.И., Гречко Г.М., Губарев А.А., Князихин К).В. Восстановление высотного хода коэффициента аэрозольного рассеяния и концентрации водяного пара по измерениям из космоса. // Сб.'" Атмосферно-оптические явления по набюдениям с орбитальных научных станции Салют-'. Тарту, 1981.

8. Антюфесв B.C. Сходимость итерационного процесса для восстановления индикатрисы. // Сб."Методы и алгоритмы статистического моделирования'', Новосибирск. 1983

9. Антюфесв B.C. Вычисление корреляционных характеристик поля яркости в стохастической среде методом Монте Карло. // ЖВМ и МФ. т.24, .N'-10. 1984

10.Антюфесв B.C. Решение интегральных уравнении разрыва потенциала методом Монте Карло. // Сб." Теория и алгоритмы статистического моделирования'', Новосибирск, 1984

П.Антюфеев B.C. Решение интегральных уравнении II-го рода со стохастическим ядром методом Монте Карло. // ЖВМ п МФ, №4. т.20, 198G

12.Антюфесв B.C., Назаралиев М.А. Гешение некоторых обратных задач атмосферной оптики с применением метода Монте Карло. // Сб." Численные методы и математическое моделирование'", Москва, 1987

13.Антюфсев B.C.. Назаралиев М.А. Обратные задачи атмосферной оптики. // Монография, изд-во ВИ СО АН СССР. Новосибирск, 1988

14-Антюфеев B.C. Преобразование формул в методе1 математических ожиданий. // Сб.'"Теория и приложения статистического моделирования". Новосибирск, 1988

lô.Antynfcev Y.S.. Na/araliev М.А. Solving of inverse problems iti atmosphere optics hv Monte Carlo method. // S.J-NAMM, Vol.3, M2, 1988

lO.Antyufeev V.S., Marsliak A.L. Using of Monte Carlo method for calculating of radiation transfer in plant canopy. // Proc. of Int. Scin. on Rad. in Canopy, Tartu, 1988

П

17 Air ! ?j {wb В С.. Kp|.: • ,n Y.T B<«< i.iiior •,< inic ;uni>fie;io

u, .. ..¡.,.>,,.^..-,11 П()Ь< pXH'Jl I'. ■' IJA. 1'js1'. \

1-Ч \H t..-},.-. В ЬХ\ M- ;>lii.Uv Л..I. Метод I' • 1 le-Kap 1o i. у p. ¡ii;ti:< перец«« а л ;л pai-thte-jbhom покр-.тча Метод

!'..tpio и ураВНСШГ переноса ЛИ р.к ТПТ> 1ЫЮГО по-

'.¡■ei/; i ' ■ ». 19s:», vi !

!'< SlUb"l'eef) Ii.C Mo--...¡шиацня .1 If Iipin \Ы моделирования nepi н.н-л излучения :•> •■ феричесьон атмосфере методом MoHic кар-м //" СП. ''Mi !<• u.i I iarii( шческою моделирования" . il"ii. U'*'i:pcb. l'.)9'l.

'.'4 Лн ,ь .¡леи B.C.. Марш,¡к A.JI. Восстановление параметров p.« nni' ibHoio нокро»«. // OA, 1990. \.S

21 Uilynleev V.S., Marsh,)k A.L. Monte Carlo method and transport t <jiiiition lit plaul < 'm.ipies. // Remote Set,-.. Environ., 1990. 31

22 AntyuUe\ Y .S... MurЧ,'!, A.L. In-.ise ршЫеш for plant paial/.M-tcrs olnn.iliiiM jj i;i,ite Sen-,. I liv.. 1990. V.33

2-> \i.iyntei v V S., Mai--' il. A.L. Milite Carlo mode] for tiaiispnit equation iiia plain '.mupy and its inversion. // Proe. of 'itli I hi Coll. Phys. Мены 11 Sign. in Rem.Sen.. Courchevel. Frame, I'101

2 1.Aiiiyufei V Y.S. Direct and inverse problems of the radiation fi .nsier in a plant canopy. // R.JNAMM. Vol.T. ,V3, 1992

2"j.AHiH!i|iecB B.C., Лактионов Л.Г.. Лунева H.А., Смерка-лов Pi. A.. Су» тин B.C., Тулпнов Г.Ф. Спектроэнергетнческие характеристики системы океан-атмосфера. //ЮК. 1992. Л'.З

2(i.Autyufeev V.S. Ou estimating the leaf-normal distribution for a plant canopy. // Remote Sens. Env., 1993, V.41

27.Antyufeev Y.S. Optimization in solving inverse problems of atmospheric optics by the Monte Carlo method. // RJ-NAMM. Vol.9. .У.6. 1994

2b Autynfrev Y.S. Solution of the generalized transport equa-!i•■!. I'll highly peaked phase function by the Monte Carlo и ¡et h- •( i. // R.LXAMM. Vol.11. .Y,2. 199G

2't Nu'vuteev Y.S., Bondarenko A.N. X-ray tomography in scattering media. // SIAM Journal of Applied Math., USA, Yo].5(i, Л12, 1996