ДОКЛАДЫ АКАДЕМИИ НАУК, 2008, том 423, № 2, с. 161-164
МАТЕМАТИКА
УДК 519.676
АЛГОРИТМЫ МЕТОДА МОНТЕ-КАРЛО ДЛЯ ВОССТАНОВЛЕНИЯ ИНДИКАТРИСЫ РАССЕЯНИЯ С УЧЕТОМ ПОЛЯРИЗАЦИИ
© 2008 г. Член-корреспондент РАН Г. А. Михайлов, С. А. Ухинов, А. С. Чимаева
Поступило 01.07.2008 г.
1. Для описания поляризационных свойств света используется вектор Стокса I = (I, Q, U, У)7, компоненты которого определяют интенсивность, степень и плоскость поляризации, а также степень эллиптичности излучения. Математическая модель переноса поляризованного излучения строится на основе феноменологического предположения о том, что в результате рассеяния ассоциируемый с "фотоном" вектор Стокса преобразуется заданной матрицей рассеяния (см., например, [1]).
Рассмотрим стационарное интегродифферен-циальное уравнение переноса излучения с поляризацией
юУФ + аФ = |а5Р(ю', ю)Ф(г, ю')йю' + ^(г, ю),(1)
или
L Ф + аФ = SF + f0
где Ф = (Ф1, Ф2, Ф3, Ф4)т - вектор-функция плотности потока частиц (векторных фотонов), иначе -вектор-функция интенсивности излучения; О -пространство единичных векторов направления, ю е О, г е Б с К3; Р(ю', ю) - матричная функция рассеяния, а = а(г) - полное сечение, а = а. + ас, ас -сечение поглощения, а. - сечение рассеяния; Г0 =
= (/о1), /о2), /о3), /о4) )Т - вектор-функция плотности распределения источника частиц. Матрица Р(ю', ю) определяется соотношением Р(ю', ю) = Цтс -- 12)К(^)Ь(-11), где Ь(-) - специальная матрица поворота, Я(-) - матрица рассеяния; - угол между плоскостью ю', . и плоскостью рассеяния ю, ю'; ;2 - угол между плоскостью рассеяния ю, ю' и плоскостью ю, ц = (ю', ю), . = (0, 0, 1) (см., например, [1, 2]).
Институт вычислительной математики и математической геофизики Сибирского отделения Российской Академии наук, Новосибирск
Новосибирский государственный университет
В атмосферной оптике рассматривается задача восстановления индикатрисы рассеяния атмосферы = Яц(ц) по наземным наблюдениям яркости неба в альмукантарате Солнца, т.е. в различных направлениях ю;, составляющих с зенитом тот же угол 6., что и направление на Солнце [3-5]. В приближении однократного рассеяния наблюдаемые значения яркости пропорциональны соответствующим значениям индикатрисы [3]. Следовательно, для оценки индикатрисы можно использовать итерационный алгоритм, в котором с помощью математического моделирования последовательно уточняется вклад однократного рассеяния в наблюдаемые яркости. Настоящая работа посвящена построению и обоснованию алгоритма такого типа, а также его реализации методом Монте-Карло, с учетом поляризации рассеянного атмосферой излучения.
2. Введем следующие обозначения: аа, ат - сечения (коэффициенты) соответственно аэрозольного и молекулярного рассеяния, а. = аа + ат [1]; (ф0, ф1, ..., ф^) = (0, ф1, ..., 180) - азимутальные углы наблюдения; (ю* , ю* , ..., ю*) - векторы направлений наблюдения; g = {&} - вектор узловых зна-
чений индикатрисы рассеяния, g
а + <аж
где ga = ga(Ц), gm = gm(Ц) - соответственно индикатрисы аэрозольного и молекулярного рассеяния [1]; g* - вектор значений искомой индикатрисы; А - значение альбедо подстилающей поверхности (см. [1, 3]); Е^), РАС?) -соответственно векторы полной яркости, яркости однократного рассеяния (при А = 0) и яркости от альбедного источника (подсветка), наблюдаемых в альмукантарате Солнца при заданной инди-
катрисе; Е^) = Cg, С = Цо т.ехр I--, где Ц0 =
V М-о/
= 0080. - параметр альмукантарата, т - полная оп-
н
тическая толща атмосферы, т.е. т = |а (г)йг, т. -
о
оптическая толща атмосферы по рассеянию, Н -
п
162
МИХАИЛОВ и др.
высота атмосферы (см. [3-5]); F2(g) = F(g) - FA(g) -- Fi(g) - яркость многократного рассеяния при A = 0; F(g*) - вектор измеренных яркостей, F(g*) = = (Ф1(- ю*), ..., Ф1(- ю*)), где Ф1 - интенсивность излучения.
При вычислении яркостей F, F2, FA вектор ga = = — g--gm достраивается (назад), линейно вос-
ва ва
полняется и нормируется по ц; таким образом получается индикатриса аэрозольного рассеяния. Матрица рассеяния строится соответственно заданным начальным матрицам аэрозольного и молекулярного рассеяний.
Аддитивный метод восстановления [3] определяется формулой
g(k) = C l[F(g*) -F2(g(k-1)) - Fa(g(k-1))] = = Gad(g(k-1)) -
F2(g) (т) F (я ) - F А ( я )
-— = о (т), --- = 1 + о (т),
Ся Ся
У я е X;
2) А мало и вследствие этого
F (я *) - FA( я) = F (я *)(1- О (А)), Уя е X.
Из предположений 1), 2) следует, что О(Х) с X, т.е. О воспроизводит X.
Для существования "неподвижной точки" оператора О и обоснования сходимости к ней вектора я(к) теперь достаточно установить, что норма якобиана ] дО !■ меньше единицы (норма любая, так
Р я-' 1
как здесь все нормы эквивалентны).
Покажем, что при т, А ^ 0 все производные дО
(2) dg
Мультипликативный метод реализуется по формуле (см. [4, 5])
(к) _ (к-1) F( g* ) _ G ( <k-i), g _ g су (k - 1) ч _ Gmult(g ) -
F(g )
становятся малыми.
i
При A Ф 0 имеем G (g) _
Cg F(g*) - Fa(g)
(3)
_ Co
F(g) - Fa(g) C 1 _Cg_FA(g)
В соотношениях (2), (3) равенство покомпонентное.
1 +
F2(g) F(g) - Fa(g) C Cg
_ R - B.
Отметим, что мультипликативный метод схо- Здесь в силу предположений 1), 2) имеем С0 = О(1),
дится за один шаг, если функция
F (g)
не зависит i +
от g. Такое свойство в какой-то степени имеет место для A = 0 и достаточно малого т. Аддитивный метод сходится сразу, если F(g) - F1(g) = const. Из геометрических соображений величина FA(g) слабо зависит от g. Поэтому целесообразен новый комбинированный метод:
F2(g) ))))C))))g))))))
= 0(1). Поэтому для исследования произ-
Э Rk
достаточно оценить соответствую-
F2 ( g )
водных тт d
щие производные от функции S(g) = 1 +
(k) (k-1) F(g*) - FA(g(k-1)) g _ g
CY (k-1K 77 / (k-1K
F(g ) - Fa(g )
_ G(g(k-1)). (4)
Здесь k, i - номера компонент векторов R и g. Дифференцированием получаем
Cg
Он сходится за один шаг, если FA(g) и
я
_ , - 0 __. - - не зависят от я. F (я) - FA(я)
3. Далее дано обоснование сходимости комбинированного метода для достаточно малых т и А.
д Sk _ d i) _ ^— _ Sk _
dgi
k Ф i,
F2°k(g) O-'gg-i '
F2kk)k gk - F2, k ( g) 2 ' Cgk
k _ i.
В этом выражении С = О(т), поэтому достаточно
д F
°перат°р О рассматриваем на множест,е X показать малость производных сравнитель-
векторов я с ограниченными компонентами (я;}:
X = {я: 0 < е < < М < .
Ясно, что если а(т) > 0, то я*, F(g*) е X. Предположим, что: 1) т мало и вследствие этого
но с т, т.е. что
3*2
dgi
_ o(т).
(5)
Это можно сделать прямым дифференцированием векторного уравнения переноса (1) с источ-
АЛГОРИТМЫ МЕТОДА МОНТЕ-КАРЛО ДЛЯ ВОССТАНОВЛЕНИЯ
163
ником однократного рассеяния. В результате такого дифференцирования получаем аналогичное дФ
уравнение для __-— (компонента вектора F полуд я-
чается фиксированием ю). Изменению подвергается лишь источник с учетом перехода от вектора я к индикатрисе линейной интерполяцией и нормировкой. Ясно, что такое изменение является равномерно ограниченным, так как дифференцирование "линейной интерполянты" по дает замену на единицу. Таким образом, с учетом предположения 1) соотношение (5) выполняется и,
д Sk д Як ч Д
следовательно, —, — = О(т). Аналогичные
д В
0 _ expcosTP(ю.' -ю*)е- _
n _ 1 N
V es -в Zn т , (i)
_ I в) L(П - (i2)n)
cos 6
X
_1
X Ra (ЦП'))ва + Rm (мП°)вш f ( (-)(i))n X -I D) L (-( i 1 ) n )en'
где (гп, юп) - цепь столкновений "фотонов" с веществом, гп - это г-координата точки столкновения, г0 = Н; юп - направление "фотона" перед рассеянием, юх = ю.; = (юп, -ю*); Яа, Ят - аэро-
разом, норма якобиана \ дО I действительно
меньше единицы для достаточно малых т и А.
Заметим, что полученный результат, по-видимому, имеет место и для некоторых вариантов задачи с не слишком малым значением альбедо ввиду его слабого влияния на яркость в альмукантарате. Нетрудно заметить, что сформулированный комбинированный метод эквивалентен итерационному алгоритму решения алгебраической системы уравнений яа = Оа(яО), где
рассуждения дают оценку -р = О(А), таким об- зольна* и молекулярная матрщы рассеяния;
д ()П-) - угол между плоскостью юп, . и плоско-
стью юп, ю* ; (г2 )П-) - угол между плоскостью юп,
ю* и плоскостью ю*, . = (0, 0, 1), N - случайный номер последнего столкновения. Векторный вес Qn определяется выражениями
01 = ql 1о,
л ЧпР(юп -1,юп)0п-1
0п =-]Г77.—)-' п - 2'
^11 (Мп-1 )
где 10 = (1, 0, 0, 0)г, |!п - 1 = (юп - х, юп), qn - скалярный вес, значение которого определяется используемой модификацией метода Монте-Карло (без поглощения, без вылета и т.д.).
Математическое ожидание Е^(,) = 1(>) = Ф(-ю*),
Е = F(,). Приведенные в [2] оценки показывают, что для реальной оптической толщи атмосферы дисперсии Б ^, определяющие среднеквадрати-ческие погрешности расчетов, ограничены.
В силу того, что во всех направлениях наблюдения в альмукантарате Солнца оптическая толща атмосферы одинакова и равна толще атмосферы в направлении падения солнечного света, вектор Стокса 1х однократно рассеянного излучения можно вычислить аналитически по формуле
ь
11 = ехр (-а I) ехр (-а( Ь -1)) Р (ю., -ю*) 10 М,
Gа(ga) _ I ga +
^ F (g *) - Fa ( g ) "J F ( g ) - Fa ( g) '
gm
Все приведенные рассуждения можно отнести и к мультипликативному методу, но лишь в случае А = 0.
4. В качестве тестовой модели атмосферы рассматривается однородный плоский слой 0 < г < Н с заданными коэффициентами поглощения, аэрозольного и молекулярного рассеяния и альбедо подстилающей поверхности. Матрица молекулярного рассеяния известна [1] и задана. Матрица аэрозольного рассеяния рассчитана по теории Ми. Будем полагать, что отражение от подстилающей поверхности происходит по закону Ламберта и излучение при этом деполяризуется. Источник солнечного излучения равномерно распределен по верхней границе слоя и является неполяризован-ным, ю. - направление падения солнечных лучей. Наблюдение яркости осуществляется на нижней отражающей границе слоя г = 0.
Для решения прямой задачи, т.е. для вычисления яркостей F(,), I = 1, 2, ... К, при заданной матрице рассеяния Яа применяется векторная локальная оценка = (^, ^, , ^) по столкновениям [1]:
г
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.