Эффективные дискретно-стохастические модификации локальных оценок метода Монте-Карло для задач лазерного зондирования рассеивающих сред Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 17, № 3, 2012

Эффективные дискретно-стохастические модификации локальных оценок метода Монте-Карло

>к

для задач лазерного зондирования рассеивающих сред*

Е.Г. КАБЛУКОВА1, Б. А. КАРГИН1'2 Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия Новосибирский государственный университет, Россия e-mail: [email protected], [email protected]

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

Ключевые слова: метод Монте-Карло, локальные оценки, лазерное зондирование.

Введение

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

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

* Работа выполнена при поддержке РФФИ (гранты № 10-01-00040 и 12-01-00034), Интеграционного проекта СО РАН № 52 и Проекта президиума РАН № 14.2.

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

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

Рассмотрим область С Е Я3, заполненную рассеивающим и поглощающим излучение веществом с коэффициентами рассеяния а3(т) и ослабления а (т), индикатрисой рассеяния д(т, такой, что

1

У д(т,^Ф = 1.

-1

Здесь ^ = (ш',ш) — скалярное произведение векторов ш', ш Е П = (ш = (а,Ь,с) : а2 + Ь2 + с2 = 1} — пространство направлений. Обозначим через q(т) величину а3(т)/а(т) — вероятность выживания кванта излучения (фотона) при столкновении с элементом вещества, с — скорость распространения излучения в среде. В точке т0 находится источник, испускающий в момент времени £ = 0 импульс излучения единичной мощности в круговом конусе направлений П0 с раствором 90 относительно оси конуса, направленной вдоль единичного вектора ш0 (рис. 1).

Требуется определить временное распределение /п*(т*,£) приходящего в точку т* излучения в направлениях ш таких, что —ш Е П*, где П* — круговой конус с раствором 9* относительно оси ш*. Таким образом, вычислению подлежит функционал

Рис. 1. Геометрическая схема задачи

в котором функция I(г, ш,£) — интенсивность излучения в точке г в момент времени £ в направлении ш. Известно (см., например, [3]), что интенсивность излучения связана с плотностью столкновений <^(х) следующим соотношением:

I (г,ш,£)а(г) = <^(г, ш,£).

Полагая а (г) = 0 при г Е Я3\С (что соответствует случаю отсутствия рассеяния и поглощения во внешней к области С части пространства Я3), плотность столкновения <^(х) в приближении лучевой оптики можно описать интегральным уравнением [3]

Г ?(г')а(г) ехр(—т(г',г))д(г', (ш',ш)) / г — г'\

^(Х) = -,-т--О Ш — 7-г X

^ ; у 2п|г' — г|2 V |г — г' |У

х^ — (V + ^^^ ^(х')^х' + /(х), (1)

х = (г,ш,£) Е X = Я3 х П х Т, Т = (0, то).

Здесь плотность первых столкновений нерассеянных частиц непосредственно от источника

/ (х) = а(г)ехр(—т (гс,г))АПо (ш)А4 ( £ — ^^

|г—г'|

/г_г'

а (г' + Зг-г Ыз,

|г — г'|

0

Ап0(ш) = ( 1, Ш Е п0, а4(£)= ( 1, £> 0,

°01 ; \ 0, иначе, \ 0, £ < 0.

2. Локальные оценки

Пусть искомый функционал 1о»(г*,£) вычисляется в виде гистограммы, т.е. оцениваются величины

и

4^*)=^ 1п* (г*,£)^, г = 1,...,п*,

^г — 1

где — узлы гистограммы, £0 = 0. Для вычисления 1(г»(г*) применяется локальная оценка [3]

N

1«(г*) = е£ д„^(1)(ж„,г*),

п=1

М1)(х„.Г*} = "(г")ехр('т'г",-г*))д2(г- А«.(,}АМ,

2п|г* — гп|2

где Е — символ математического ожидания, хп = (гп,шп,£п) — цепь Маркова с плотностью первых столкновений р1 (х) и переходной плотностью р(х',х), N — случайный номер обрыва марковской цепи, £ = £п + |гп — г* |/с, з = (г* — гп)/|г* — гп|,

* Г 1, если £ Е (¿г_ 1,£г], » / \ Г 1, если — ш Е П*,

АДг) = < „ Ао» (ш) = < „

у ' 10, иначе. о у ' 0, иначе.

Весовые множители Qn в соответствии с теорией весовых методов Монте-Карло [3, 4] определяются выражениями

ГЛ ^(х) ^ ГЛ к(Хп-1,Хп) АТ

Ql =—гт, Qn = Qn-\—t-г, П = 2,...,М,

Р1(х) р(Хп-1,Хп)

здесь к(х',х) — ядро уравнения (1).

Легко видеть, что (1) эквивалентно уравнению

ф) = к 2ф) + к! (х) + ! (х), (3)

где К — интегральный оператор с ядром к(х',х). Локальная оценка ) для расчёта /П* (г*), основанная на представлении (3), называется двойной локальной оценкой и определяется формулой

/п* (г*¿) = ! к2(х',х*)(р(х')<1х',

X

k2(x,,x*) = J к(х', р)к(р,х*)<р. (4)

Хб

Последний интеграл берется по области С = {р(в) € С : р(в) = г* + шв, ш Е П*, 5 > 0}. В случае применения двойной локальной оценки вместо (2) имеем соотношение

Ф(г*) = Е£ = Е (^^(х^г*) + £Qnh(¿)(xn,г*,pn?j , (5)

(2) * _ д(г)д(р)а(р)в-т(г'р)-т(р'г*)д(г, (ш, Шр))д (р, (шр, ^)) ^ (х,г*,р')= (2п)21р - г12р(р) Х

хАп*( ^^ А (ь +1Р - г1 + 1г* - Р1

| г* - р|

Здесь шр = --г, р(р) — произвольная плотность распределения промежуточного слу-

р - г |Р - г|

чайного узла р, который выбирается так, что-(г* —р)/|г* — р| € П*. Свобода выборар(р) позволяет оптимизировать оценку (5) с целью уменьшения трудоёмкости алгоритма.

Оценка (5), (6) имеет бесконечную дисперсию, поэтому на практике вместо нее вычисляется смещенная оценка интенсивности с конечной дисперсией. В этом случае для точки столкновения (гп,шп,Ьп) дополнительный случайный узел р разыгрывается в области С\Бе, Бе = {г Е Я3 : |гп — Г| < е}, где е — некоторое наперед заданное положительное число. Отметим, что относительную погрешность вычисления величины /п* (г*), обусловленную вырезанием шара Бе, можно приближенно оценить (см. [1]) величиной

д(г)(1 - ро),

1 - д(г)№

-(1 - ехр(-а(г)е)),

1

где ^о = / — средний косинус угла рассеяния в точке г. Ниже в табл. 4

-1

(см. раздел 5) представлена зависимость значений /п* (г*) от радиуса е вырезаемого шара Бе.

3. Модифицированная оценка

Рассмотрим некоторую область С£ Е С такую, что С Е С£. Для точек столкновения = (гга,^га,£га) таких, что гп Е С£, вычисление интеграла (4) будем проводить по некоторому наперед заданному числу К случайных узлов интегрирования. В этом случае функция Л,(2) в (6) заменяется на следующую:

1

к

(2),

К

Для жга таких, что гп Е С\С£, вклад в двойную локальную оценку вычисляется по одному случайному узлу интегрирования р (К = 1). В данной работе для тестовых расчётов линейный размер области С£ выбирался порядка длины свободного пробега частицы в веществе, заполняющем область С.

Для определения оптимального числа узлов интегрирования в области С£ рассмотрим задачу вычисления интегральной по времени интенсивности излучения, приходящего в точку приемника г* в заданном конусе направлений. Воспользуемся формулой полной дисперсии и методикой определения оптимальных параметров метода расщепления [4].

Считаем, что , г*) = 0. Представим модифицированную двойную локальную

оценку интенсивности излучения в виде

N

N

£(к) = 5] ^й(2)(ж„,г*|г„ Е С\С£) + ^ д„^(3)(х„,г*|г„ Е С£).

п=2

п=2

Обозначим через £ = (х1,...,^) последовательность фазовых координат точек столкновений фотона с элементами вещества, п = (р1,..., р1,... , рN,... , PN) — случайные узлы интегрирования р^ Е С для последовательности точек столкновения £ (з = 1, если гп Е С\С£, з = К, если гп Е С£). Формула полной дисперсии случайной величины £(к) имеет вид [4]

Б£(к) = Бс Еп (С(К)|С) + Ес (£(к)|С) =

" N N

^ ^(2)(Х„,Г*|Х Е с,г„ Е С\а) + ^ д„^(3)(х„,г*|х Е с,г„ Е а)

п=2

БЛ Е

С Еп

N

п=2

+

+ЕС Бп

^д„^(2)(ж„,г*|ж Е С,г„ Е С\С£)

п=2

+ Ес Бп

N

^д„^(3)(ж„,г*|ж е с,г„ Е а)

п=2

Последнее равенство справедливо, так как Бп(£(к)|() = Бп(£(к)|С,гга Е С\С£) + Бп(С(к)|С,гп Е С£) для независимых р^. Используя независимость и одинаковую рас-пределённость компонентов вектора п, а также формулу (7), получим

Бпй(3)(ж„,г*|ж„ Е С,Гп Е (?£) = 1 Бпй(2)(ж„,г*|ж„ Е С,Гп Е (?£).

К

Введем обозначения:

£>1 = Бс Еп (£(к)|С) + Ес Бп

N

]>]д^(2)(Хп,г*|х Е С,г„ Е С\С£)

п=1

N

Б2 = есб^ Qn^2)(хп,г*1х Е (,гп Е Се).

п=1

Тогда

Б£(к) = Б1 + Б2/К.

Пусть Ь1 — среднее время моделирования точек столкновения хп, п = 1,... , N, Ь2 — среднее время моделирования одного дополнительного узла интегрирования рп для каждого хп, п = 1,...^, и вычисления функционала &2\хп,г*), 11 и 12 — среднее отношение числа точек столкновения хп из областей С\Се и Се к общему их числу (¡1 + ¡2 = 1). Тогда среднее время вычисления одного выборочного значения случайной величины

Ь(К) = ^(к + ¡2) + Ь2(к + К12) = Ь1 + Ь2к + КЬ2\2 = Ь 1 + Кг2,

Ь-1 = Ь1 + Ь211, Г 2 = Оптимальное значение К минимизирует величину трудоёмкости

Б(к) = Ь(к)Б£(к) = (11 + КГ2 )(Б1 + В2/К).

Вычислив производную и учитывая положительность величин П1, 02, Ь1, Ь2, получим, что оптимальное число К приближенно равно целому положительному числу, наиболее близкому к выражению [4]

- М- <8)

Величины В1,В2,Ь 1,Ь2 можно оценить по результатам специальных предварительных расчётов, вычисляя оценку интенсивности дважды: для параметра К = 1 и некоторого заданного К. В этом случае, зная среднее время вычисления одного выборочного

значения ^ 1 ь(1)

и )

ь(к) случайных величин £(1) и £(к), имеем

ь(к) - ь(1) (1)

Ь2 = ТК-1%, Ь1 = Ь - ^

следовательно

г _ Ь(1)К - ь(к) г ь(к) - ь(1) ь 1 = К - 1 , Ь2 = К - 1 . Величину П1 можно оценить по формуле

Dx = Щ(к) - D2/K.

4. Дискретно-стохастическая версия "выборки по важности"

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

р(р) = -^г^-;гт, (9)

2п(1 - cos в*)' KJ

не отражающая характер поведения подынтегральной функции из (4) для некоторых точек x Е X. Имея простую процедуру, определяющую, для каких x Е X функция h( ;(x,r*,p) из (4) пропорциональна функции (9), а для каких x Е X не пропорциональна, можно для моделирования дополнительных узлов р комбинировать плотность (9) с плотностью, в некотором смысле более близкой к функции h(2)(x, r*, р).

Далее рассмотрим некоторую область G£1 Е G такую, что G Е G£1 (возможно совпадающую с областью Ge, см. раздел 3), и будем моделировать дополнительые узлы интегрирования по плотности (9) для точек столкновения {xn : rn Е G\G£1}. Для {xn : rn Е Gei} в качестве плотности распределения дополнительного узла интегрирования р выберем кусочно-линейное приближение функции h(2)(x, r*, р) по переменной р Е G и вычислим интеграл (4) с помощью метода дискретно-стохастической версии "выборки по важности" [4].

Вычисление интеграла проводится в системе полярных координат с центром в r* и осью w*. В этом случае р = (f , Д,() и область G содержится в некотором параллелепипеде со сторонами (f , Д,() Е [fmin,fmax] х [0, cos 0*] х [0, 2п]. Плотность pj(р) для фиксированной точки столкновения {xj Е X : rj Е G£1} строится как кусочно-линейная аппроксимация функции h(2)(xj, r*, р) по переменной f для заданных реализаций Д,-, (

случайных величин f, (, равномерно распределённых в области [0, cos 0*] х [0, 2п]:

р(р) =

P(f |xj,f j ,(j)

2n(1 - cos 0*)'

Пусть на отрезке [fmin, fmax] задана неравномерная сетка |rmj, m = 0,..., M, с шагом hm = rm+i — fm. Тогда для фиксированных значений Д^-, случайных величин Д, ^ плотность p(f |xj , Д,-, ) представляется как

M X „ _ „ ч ,

p(f |xj, Д, ) = Pm/m(f ), /m(f ) = Xm ( ~ j Km (10)

m=0 \ hm J /

где

Г 1 + y при — 1 < y < 0, х(у) = в(1)(y)= ^ 1 — У при 0 < y < 1,

0 иначе,

ß — сплайн первой степени,

hm

hm-1 при f - fm < 0, hm при f - fm > 0,

Кт = J Хт(г)^г = Т^т-1 + Рт = СК^^Ж,- , Г*,Гт,Д?, ^) —

К

значения аппроксимируемой функции в узлах {гт} при фиксированных Д^-, ^. Нормирующая константа

, M

C = 1 / Ymh(2)(xj, r*, fm, Дj , (j).

' m=0

Моделирование координаты г дополнительного узла р проводится методом суперпозиции: по вероятностям Рд выбирается номер моделируемой плотности Р(т = к) = Рд и строится случайная величина по плотности /т(г ).

Для неравномерной сетки и плотности /т(т) с вероятностью д = hm-l/(hm-l + hm) моделируется точка из промежутка [г т - hm-1, гт], с вероятностью (1 - д) — из промежутка [гт, Гт + hm] по формулам г ^ = гт + hm-l(^aj - 1) и г ^ = гт + hm(^aj - 1) соответственно. Здесь aj — равномерно распределённые случайные величины на (0,1).

Для равномерной сетки {гт} с шагом h существует более простая формула моделирования случайной величины по плотности fт(г): г j = + а2 - 1) + гт, здесь а^,в = 1, 2 — равномерно распределённые случайные величины на (0,1).

В численных экспериментах число узлов аппроксимации {г т} выбиралось малым, при этом один из узлов находился вблизи проекции точки Гj на луч Шj = (г, фj), а узел г0 — на пересечении границы области С, ближайшей к детектору, и луча Шj.

Для быстроты моделирования случайных узлов р число узлов сетки М должно быть мало. Вместе с тем такая аппроксимация должна отражать поведение подынтегральной функции, т. е. узлы аппроксимации {гт} должны быть близки к экстремумам приближаемой функции. Это условие осложняет использование данной плотности, но даже приближённое положение узлов аппроксимации в некоторых случаях дает выигрыш в трудоёмкости вычисления интенсивности излучения, интегрированной по времени. Это можно наблюдать в нижеприведённых тестовых расчётах.

5. Результаты численных экспериментов

Для иллюстрации эффективности предложенного алгоритма представлены результаты вычисления оценки интенсивности излучения Iq* (r*) и временного распределения интенсивности излучения lQ*(r*,t) для плоского слоя G = {(x,y,z) Е R3 : h < z < H}, h = 0.5 км, H = 0.9 км, заполненного веществом с коэффициентом ослабления a(r) = Е = 50 км-1, вероятностью выживания q = 0.95 и индикатрисой рассеяния Хеньи — Гринстейна g(ß) = (1 — а2 )/2(1 + а2 — 2aß)3/2, ß Е (—1, +1), с параметрами а = 0.6, 0.7, 0.8, 0.9.

Предполагается, что детектор излучения находится в начале координат r* = (0, 0, 0), ось детектора совпадает с ш* = (0,0,1). Источник излучения, находящийся в точке r0 = (—0.7 км, 0,0), излучает импульс единичной мощности в круговом конусе с раствором в0 = 20" в направлении ш0 = (V2/2, 0, у/2/22).

Для вещества с индикатрисой рассеяния Хеньи — Гринстейна с параметром а = 0. 8 и детектора с углом раствора В* = 10 проведено сравнение двойной локальной оценки интенсивности излучения Iq* (r*) (случай K = 1) c модифицированной двойной локальной оценкой. Во втором случае выбиралась область

G£ ={r Е G : > cos30^

для точек столкновения {xn Е X : rn Е G\Ge}, K = 1, а для точек {xn Е X : rn Е G£} число K дополнительных узлов интегрирования pk варьировалось. Все дополнительные узлы pk моделировались по плотности (9). В табл. 1 показана трудоёмкость S = a2T рассмотренных оценок, построенных по n = 109 траекториям. Здесь T — среднее время моделирования одной случайной траектории, а = \JDIq* (r*)/n — среднеквадратичное отклонение оценки, t — время вычисления оценки по n = 109 траекториям. Из приведённых данных следует, что трудоёмкости S оценок примерно одинаковы для K = 200^250 и значительно меньше трудоёмкости двойной локальной оценки (K =1). Проводились

Таблица 1. Трудоёмкость 5, среднеквадратичное отклонение а модифицированной двойной локальной оценки интенсивности излучения /л* (г*) для различного числа К

К 1 25 50 100 150 180 200 230 250

5 ■ 107 33 3.4 4.2 3.9 4.4 2.7 1.9 2.0 1.8

г ■ 10-4,с 5.34 5.97 6.9 8.6 10.5 11.4 12.4 13.4 14.16

а ■ 106 7.84 2.39 2.47 2.12 2.06 1.54 1.23 1.21 1.09

также расчёты согласно формуле (8), где величины ^ и Д2 оценивались по числу траекторий п = 108. Согласно (8) получено оптимальное КорЬ « 205, которое согласуется с вычислительным экспериментом.

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

/ г_г*

V(г) = ( 1-г ,ш*

О1 = (г е О : V(г) < 71) , О2 = (г е О : 71 < V(г) < 72) ,

Оз = (г е О : 72 < V(г) < 7э) , О4 = (г е О : V(г) > 73).

В представленных расчётах 71 = сов(0.087), 72 = сов(0.052) 73 = сов(0.009). Для точек |жга е X : гп е О1} функционал (7) оценивался по одному дополнительному случайному узлу р (К = 1). Для |жга е X : гп е О 2,О 3, О4} число дополнительных узлов К2, К3, К4 выбиралось исходя из нескольких оценок интенсивности излучения /Л*(г*, £), построенных по числу траекторий п = 107-108, с различными значениями К2,К3,К4. Модифицированная двойная локальная оценка с плотностями дополнительных узлов (9) и (10) обозначена как МДЛО1, с плотностью дополнительного узла (9) — как МД-ЛО2. Эти оценки сравнивались с двойной локальной оценкой (ДЛО) (6) и локальной оценкой (ЛО) (2).

На рис. 2 представлены оценки временного распределения интенсивности излучения /(г* (г*, £) для вещества с индикатрисой рассеяния Хеньи — Гринстейна с параметром а = 0.8 и детектора с углом раствора 0* = 30" на временном интервале [4.6 мкс, 7.2 мкс] с шагом гистограммы 0.1 мкс. Число дополнительных узлов интегрирования рд в оценках МДЛО1 и МДЛО2 равно К2 = 20, К3 = 30, К4 = 50. Жирной линией обозначены среднеквадратичные отклонения а^, г = 1,...,п^, двойной локальной оценки, тонкой линией — а, г = 1,..., п^, модифицированных двойных локальных оценок МДЛО1 и МДЛО2.

В табл. 2 приведены значения локальных оценок интенсивности излучения /л*(г*), определяемой детектором с углом раствора 0* = 1' и параметром индикатрисы рассеяния а = 0.7, 0.8, среднеквадратичных отклонений а, времени вычисления оценок £ по числу траекторий п = 109 (исключение — двойные модифицированные оценки МДЛО1 и МДЛО2 для задачи с параметром а = 0.7 — п = 108) и трудоёмкости алгоритмов Б. Видно, что для малых углов раствора детектора (порядка 0* = 1' и меньше) локальная оценка, построенная по числу траекторий п = 109, не дает удовлетворительной точности вычисления временного распределения интенсивности /(г* (г*, £), для нее велики также а и Б. Трудоемкость Б модифицированных двойных локальных оценок МДЛО1 и МДЛО2 в несколько раз меньше трудоёмкости двойной локальной оценки, а для а = 0.7 эта разница составляет более одного порядка.

1е-07-

*

-Л-

а

1е-08-

~я—д-

*

- * ж2 х ДЛО

: ± о мдло!

1

ж.

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1111111

4.9 5.4 5.9 6.4 6.9

мкс

5

1е-07-

1е-08-

* =1=:

X ДЛО о МДЛ02

-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-

4.9 5.4 5.9 6.4 6.9

мкс

Рис. 2. Оценка временного распределения интенсивности излучения /^Ц(г*,г) при а = 0.8, в* = 30": а — ДЛО, МДЛО1, б -ДЛО, МДЛО2

а

Таблица 2. Результаты вычисления интенсивности излучения /п* (г*) различными методами для в* = 1' (а — параметр индикатрисы рассеяния Хеньи — Гринстейна)

Метод /п*(г*) ■ 107 а ■ 107 г, с 5 ■ 1013

а= 0.8

ЛО 1.48 0.20 15025 60

ДЛО 1.47 0.03 30880 2.8

МДЛО1 1.44 0.009 42970 0.8

МДЛО2 1.45 0.02 42430 1.7

а= 0.7

ЛО 1.14 0.118 20200 28

ДЛО 1.19 0.043 42500 7.7

МДЛО1 1.12 0.012 7300 0.11

МДЛО2 1.14 0.010 7150 0.08

В табл. 3 показана трудоёмкость S локальных оценок для различных углов раствора детектора 0* и параметров индикатрисы рассеяния a = 0.7, 0.8. Приведенные данные подтверждают, что трудоёмкость оценок МДЛО1 и МДЛО2 меньше трудоёмкости S оценки ДЛО для всех представленных параметров 0*, a.

В ряде случаев применение аппроксимации подынтегральной функции в качестве плотности моделирования дополнительного случайного узла интегрирования р позволяет незначительно сократить трудоёмкость модифицированной оценки МДЛО1 по сравнению с МДЛО2, где в качестве плотности дополнительного узла интегрирования

Таблица 3. Трудоемкость S рассмотренных методов для 0* = {10, 30', 1', 30"} с параметром индикатрисы рассеяния Хеньи —Гринстейна a = 0.7, К2 = 40, К3 = 60, К4 = 100 и a = 0.8, К = 20, K3 = 30, К = 50

Метод 10 30' 1' 30''

a= 0.7, К = 40, Кз = 60, К4 = 100

ЛО 2.3 10-8 7.0 ■ 10-9 — —

ДЛО 54 ■ 10-8 49 ■10-9 7.7 ■ 10" 13 3.4 10" 14

МДЛО1 2.2 10-8 3.0 ■ 10-9 0.11 ■ 10" -13 0.13 10" 14

МДЛО2 1.8 10-8 2.7 ■ 10-9 0.08 ■ 10" -13 0.07 10" 14

a= = 0.8, К = 20, Кз = 30, К4 = 50

ЛО 6.9 ■ 10-8 1.6 ■ 10-8 — —

ДЛО 1.4 10-6 4.7 ■ 10-7 2.8 ■ 10" 13 1.7 10" 14

МДЛО1 2.2 10-7 2.7 ■ 10-8 8.2 ■ 10" 14 2.5 10" 15

МДЛО2 1.7 10-7 1.6 ■ 10-8 1.7 ■ 10" 13 4.1 10" 15

Таблица 4. Оценки интенсивности излучения /п* (г*) и среднеквадратичного отклонения а для различных радиусов е шара В£; 0* = 30", а = 0.6, 0.7, 0.8, 0.9

ДЛО МДЛО1 МДЛО2

е (In* (r*) ± a) ■ 108 (In* (r*) ± a) ■ 108 (In* (r*) ± a) ■ 108

a = 0.6

4 10"5 2.064 ± 0.053 2.052 ± 0.017 2.046 ± 0.025

1 10"4 2.064 ± 0.053 2.049 ± 0.016 2.045 ± 0.025

4 10"4 2.041 ± 0.051 2.022 ± 0.014 2.009 ± 0.014

a = 0.7

1 10"4 2.590 ± 0.044 2.767 ± 0.034 2.709 ± 0.025

2 10"4 2.590 ± 0.044 2.756 ± 0.033 2.703 ± 0.025

4 10"4 2.589 ± 0.044 2.728 ± 0.027 2.686 ± 0.024

8 10"4 2.561 ± 0.041 2.679 ± 0.017 2.518 ± 0.017

a = 0.8

5 10"5 3.343 ± 0.058 3.527 ± 0.132 3.395 ± 0.025

1 10"4 3.333 ± 0.058 3.395 ± 0.023 3.395 ± 0.025

5 10"4 3.317 ± 0.057 3.349 ± 0.022 3.339 ± 0.021

1 10"3 3.232 ± 0.051 3.278 ± 0.020 3.328 ± 0.019

a = 0.9

4 10"5 3.390 ± 0.250 3.028 ± 0.053 3.263 ± 0.103

1 10"4 3.390 ± 0.250 3.028 ± 0.053 3.244 ± 0.103

4 10"4 3.349 ± 0.248 3.002 ± 0.052 3.216 ± 0.103

1 10"3 3.266 ± 0.245 2.955 ± 0.051 3.129 ± 0.099

1 1 | 1 1 1 1 | 1 1 1 1 | 1 1 I I | I I I 1 | 1 1 1

5 5.5 6 6.5 7 7.5

МКС

Рис. 3. Ыодифицированная двойная локальная оценка временного распределения интенсивности излучения для различных радиусов в шара Б£; в* = 30", а = 0.9

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

Для обеспечения конечности дисперсии двойных локальных оценок вычислялись смещенные оценки интенсивности излучения, в которых для точек столкновения {хп Е X : Гп Е О} дополнительные узлы интегрирования моделировались в области О\В£, В£ = {Г Е Я3 : |Гп — г| < г}. В табл. 4 (см. с. 80) показано поведение значений двойных локальных оценок интенсивности излучения /п* (г*) в зависимости от линейного размера г шара В£. Рассмотрена задача с углом в* = 30// и параметрами индикатрисы рассеяния Хеньи — Гринстейна а = 0.6, 0.7, 0.8, 0.9. Для данной задачи г ~ 4 • 10-4 приводит к незначительному смещению (порядка погрешности а) оценки /п*(г*) и уменьшает среднеквадратичное отклонение а(г) в оценке временного распределения интенсивности /п* (г*, Это отражено на рис. 3, где представлена модифицированная двойная локальная оценка МДЛО2 временного распределения интенсивности излучения /(г*(г*,^) для задачи с параметрами в* = 30//, а = 0.9, К = 10, Кз = 20, К4 = 30.

Из сравнения двойной локальной оценки с модификациями МДЛО1 и МДЛО2 следует, что применение метода "расщепления" позволяет значительно снизить трудоёмкость алгоритма. Это наиболее заметно для задач с "гладкими" индикатрисами рассеяния (а = 0.6 -г 0.8) и достаточно большими углами раствора приемника в*. Но и для вытянутых индикатрис, а также малых углов в* предложенные модификации в сравнении с двойной локальной оценкой сохраняют свое преимущество.

Список литературы

[1] Джетыбаев Е.О. К решению нестационарного уравнения переноса излучения методом Монте-Карло. Новосибирск, 1982 (Препр. ВЦ СО АН СССР; № 346).

[2] Креков Г.А. Метод Монте-Карло в проблемах атмосферной оптики // Оптика атмосферы и океана. 2007. Т. 20, № 9. С. 826-834.

[3] Марчук Г.И., Михайлов Г.А., Назаралиев М.А. и др. Метод Монте-Карло в атмосферной оптике. М.: Наука, 1976.

[4] Михайлов Г.А., Войтишек А.В. Численное статистическое моделирование. Методы Монте-Карло. М.: Академия, 2006.

Поступила в 'редакцию 2 апреля 2012 г.