авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 || 3 | 4 |   ...   | 15 |

«Посвящается пионерам освоения космоса Предисловие Фактически в настоящее время закладываются основы решения фундамен- тальных проблем, связанных с ...»

-- [ Страница 2 ] --

– Алгоритм 1 удобен для анализа, чтобы показать наличие -анизотропии в направлении s0 во всех кратностях рассеяния. Но для организации расчетов он неудобен, так как на каждой итерации нужно рассчитывать n, n, n и громоздко вычисление правых частей. Однако следует ожидать, что, начиная с некоторой итерации, можно пренебречь вкладом компонент n в направлении s0 и считать n+1 гладкой функцией, а при вычислении функции источника учитывать только -вытянутость индикатрисы, т. е.

B n+1 = B n+1 + (z) c(z) n+1.

Выделив функцию рассеяния (n+2)-кратности Nn+2 = n+2 +Nn+3, получим задачи dn+ = 0, + n+2 = B n+1, n+2 n+2 = R n+1 ;

dz 0 H d Nn+3 + N = B Nn+3 + B n+2, n+ dz N = 0, Nn+3 = R Nn+3 + R n+ n+ 0 H с гладкой функцией n+2 n+2 и гладким интегралом B n+2. В результате суммарная функция рассеяния порядка кратности выше n будет гладкой 1.1. Скалярная плоская задача функцией Nn+1 = n.

n =n+ Для однородного слоя аналитически вычисляется вклад n z cz n = S exp = 0 n! n=1 n= z(1 c) z = S exp exp.

0 Эта величина совпадает со значением разности между прямыми потоками нерассеянного излучения от источника, ослабленными на оптическом пути (1 c )z и z соответственно, т. е. выделение -анизотропии в индикатрисе приводит фактически к изменению оптической толщины. С помощью выра жения (1.42) для n можно проводить оценки малости вклада и вовремя отключать учет -анизотропии в направлении s0.

АЛГОРИТМ 2. Другой способ «обрезания» ряда для функций n состоит во введении коэффициентов tot (z) = 1 c(z) (z) посредством следующей процедуры. Пусть функция рассеяния Nn кратности порядка n и выше является решением задачи (1.41) c функцией источника B Nn = B Nn + c(z)(z)Nn.

Тогда dNn + tot (z) N = B N + B n1, n n dz Nn = 0, Nn = R Nn + R n1.

0 H Выделим функцию рассеяния n-кратности как решение задачи dn = 0, + tot (z) n = B n1, n n = R n1.

dz 0 H Кратное рассеяние Nn+1 = n+1 + Nn+2, начиная с (n + 1)-го порядка и выше, определяется задачей dNn+1 + tot (z) N = B Nn+1 + Bn, n+ dz Nn+1 = 0, Nn+1 = R Nn+1 + Rn.

0 H Вычисление интеграла B n1 от двух -анизотропных функций n и необходимо регуляризовать. Приближение n = n + n ищется через 50 Глава 1. Одномерные плоские задачи решение задач d n = 0, + tot (z) n = Bn1, n n = R n1 ;

dz 0 H d n = 0, = 0;

+ tot (z) n = Bn1,, n n dz 0 H z z 1 (t) c(t) n1, (t) exp tot (z ) dz dt.

n (z, 0, 0 ) = 0 t Вклад рассеяния (n + 1)-кратности находим из задачи dn+1 = 0, + tot (z) n+1 = Bn, n+1 n+1 = R n dz 0 H с источниками – гладкими функциями – B n B n + (z) (z)n (z, cos 0 ), R n = Rn + (q 0 /)n (H).

Поэтому функция n+1 не содержит -анизотропии в угловом распределении и все следующие приближения, а значит, и Nn+2 будут гладкими функциями.

По существу, при введении tot (z) произошло замыкание последовательных приближений в составляющей решения, отвечающей -анизотропии в направ лении s0, которое свелось к тому, что в выражения для функций источника B n+1 и B N теперь не входит слагаемое с -анизотропным членом типа (z) (z) n (z, cos 0 ), т. е. произошло уменьшение источника и функции N, а в левой части уравнения появился коэффициент tot (z), дающий поправку на уменьшение оптической толщины, приводящее к уменьшению ослабления излучения и тем самым к увеличению N.

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

1. Если замена индикатрисы (z, cos ) ее аппроксимацией (1.38) вводится в исходной задаче для функции (z,, ), то решение 1 оказывается искаженным: во-первых, оно занижено в области индикатрисного максимума в направлении внешнего потока;

во-вторых, оно завышено по всем направлениям из-за того, что. Суммарная интенсивность окажется неверной, особенно для оптически тонких слоев. Такой прием используется в методе подобия, который часто применяется для решения уравнения переноса «быстрыми»

методами (двухпотоковое приближение, метод -Эдингтона и т. п.).

2. Если замена индикатрисы (z, cos ) делается в задаче для функции многократного рассеяния d = 0, + tot = B + B 0, = R + R 0, dz 0 H 1.1. Скалярная плоская задача где B вычисляется с, а B 0 – с исходной индикатрисой, то значение – 1, определяющее профиль решения при малых, вычисляется неверно с.

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

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

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

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

1.1.6.1. Учет анизотропии аэрозольно-молекулярного рассеяния. Рас смотрим задачу переноса монохроматического излучения в вертикально неоднородном плоском слое, освещаемом внешним параллельным потоком.

Интенсивность излучения (z,, ) на уровне z [0, H] в направлении 52 Глава 1. Одномерные плоские задачи s = {, } находится как решение краевой задачи d = S (s s0 ), (1.44) + t (z) = Sz, = R dz 0 H с интегралом столкновений Sz s (z) (1.45) (z,, )(z, cos ) ds.

Суммарное ослабление t (z) = s (z) + a (z). Суммарные коэффициенты рассеяния и поглощения (1.46) s (z) = a,s (z) + R (z), a (z) = a,a (z) + R,a (z) содержат аэрозольные (a,s (z) и a,a (z)) и молекулярные (R (z) и R,a (z)) компоненты. Суммарная индикатриса рассеяния определяется по формуле a,s (z) l (z) a () + R R () (1.47) (z, ) = s (z) s (z) l через индикатрисы аэрозольного a () и рэлеевского рассеяния R () с нормировкой ds = 2 (z, ) d = 1. (1.48) Вводя оптическую толщину z (z) = t (z ) dz, d = t (z) dz, в уравнении (1.44) можно перейти к пространственной переменной :

d = S (s s0 ), (1.49) + = S, = R, d 0 H S ( ) (,, )(, cos ) ds, (1.50) альбедо акта рассеяния ( ) = s ( (z))/t ( (z)).

С учетом представления (1.47) Bz () Sz = l (1.51) a,s (z) a () + R (z) R () ds, B () S = l (1.52) a,s ( ) a () + R ( ) R () ds, 1.1. Скалярная плоская задача где альбедо актов аэрозольного и молекулярного рассеяния a,s ( (z)) R ( (z)) a,s ( (z)) =, R ( (z)) =, t ( (z)) t ( (z)) так что (1.53) Sz = Saz + SRz, S = Sa + SR, Saz a,s (z) a () ds, SRz R (z) R () ds, l Sa a,s ( ) a () ds, SR R ( ) R () ds.

l В аэрозольных индикатрисах рассеяния, заданных по высотным l-зонам, l = 1 L, выделяем -анизотропию:

a () = cl (1 cos ) + l a ().

l l (1.54) l Полная исходная аэрозольная индикатриса рассеяния с нормой a,bx () Al l = a,bx () d cos a,bx перенормируется l a,bx () l a () =, 2 Al a,bx чтобы удовлетворить условию (1.48). Как-то сглаженная (обычно графиче ски или с помощью аналитической аппроксимации) исходная аэрозольная l индикатриса рассеяния a,bx () c нормой Al l = a,bx () d cos a,bx также перенормируется:

l a,bx () l a () =.

2 Al a,bx Чтобы определить весовые коэффициенты в формуле (1.54), полную исходную аэрозольную индикатрису рассеяния запишем в виде, где выделена исходная «сглаженная» аэрозольная индикатриса рассеяния:

a,bx () = f l (1 cos ) + a,bx ().

l l (1.55) Значение весового параметра f l ищем, исходя из условий нормировки индикатрис рассеяния:

2Al a () = f l (1 cos ) + 2 Al a () ;

l l a,bx a,bx 54 Глава 1. Одномерные плоские задачи Al fl (1 cos ) + a,bx a () ;

l l (1.56) a () = 2Al Al a,bx a,bx l fl 1 Aa,bx = +, 2Al 2 Al 2 a,bx a,bx откуда следует f l = Al Al 0;

(1.57) a,bx a,bx причем f l = 0, если «сглаженная» аэрозольная индикатриса рассеяния совпадает полностью с исходной аэрозольной индикатрисой рассеяния.

Сопоставляя выражения (1.54) и (1.56), получаем весовые коэффициенты:

Al Al fl 1 a,bx a,bx 1 (1 l ).

l = cl = (1.58), = = Al 2Al Al 2 a,bx a,bx a,bx Суммарная индикатриса рассеяния (1.47) с выделением -анизотропии представляется в виде, аналогичном (1.55):

(z, ) = v(z)(1 cos ) + (z, ) (1.59) с коэффициентом a,s (z) cl (1.60) v(z) = a,s (z) + R (z) и сглаженной частью a,s (z) R (z) l a () + l (1.61) (z, ) = R ().

a,s (z) + R (z) a,s (z) + R (z) Замечание: в выкладках используется функционал (1 cos ) d cos = 1. (1.62) Легко показать, что для представления (1.59) выполняется условие нормировки. Действительно, a,s (z) a,s (z) R (z) (1 l ) + l + 1= (1.63).

a,s (z) + R (z) a,s (z) + R (z) a,s (z) + R (z) Если в суммарной индикатрисе рассеяния и суммарном коэффициенте рассеяния не выделять рэлеевскую компоненту, т. е. принять l R (z) = 0, (1.64) s (z) = a,s (z), (z, ) = a (), но оставить выделение -анизотропии, то в (1.59)–(1.61) (z, ) = cl (1 cos ) + l a (), l v(z) = cl, (z, ) = l a ().

l (1.65) 1.1. Скалярная плоская задача Если -анизотропию не выделять, то l l l = 1, cl = 0, (1.66) a () = a (), l l (1.67) (z, ) = a () = a ().

Если выделить рэлеевскую индикатрису рассеяния, но -анизотропию не учитывать, то l = 1, cl = 0, v(z) = 0, (1.68) (z, ) = (z, ).

Обратим внимание на представление суммарной аэрозольно-молекулярной индикатрисы рассеяния, когда пространственной координатой является оп тическая толщина. В формулах (1.60) и (1.61) перейдем к послойной аппроксимации по высоте и умножим числители и знаменатели на величину шага zk :

a,s,k zk cl, v( (z)) = a,s,k zk + R,k zk a,s,k zk R,k zk l a () + l ( (z), ) = R ().

a,s,k zk + R,k zk a,s,k zk + R,k zk Эти выражения можно переписать в терминах оптических толщин a,s,k cl, v(k ) = a,s,k + R,k a,s,k R,k l a () + l (k, ) = R (), a,s,k + R,k a,s,k + R,k при этом s,k = a,s,k + R,k.

Но можно взвешивать аэрозольную и молекулярную индикатрисы рассе яния с помощью альбедо акта рассеяния, разделив в формулах (1.60) и (1.61) числитель и знаменатель на суммарную экстинкцию t (z). Тогда a,s ( (z)) l v( (z)) = c, s ( (z)) ( (z)) l l ( (z)) ( (z), ) = a,s a () + R R (), s ( (z)) s ( (z)) где s ( (z)) = a,s ( (z)) + R ( (z)).

1.1.6.2. Последовательность вычислений с учетом -анизотропии. Как уже отмечалось выше, представление (1.54) с меньшими погрешностями целесообразно применять, начиная с третьей кратности рассеяния. Опишем последовательность вычислений для решения задач (1.44) и (1.49). Выделяем прямое излучение 0 ( = + 0 ) – решение задачи Коши:

– d = S (s s0 ), + t (z) 0 = 0, =0 (1.69) 0 0 H dz 56 Глава 1. Одномерные плоские задачи в явном виде (z) 0 (z,, ) = S exp (s s0 ), (1.70) одинаковом для задач (1.44) и (1.49). С учетом выражения (1.54) из (1.51) получаем Bz () Sz = Bz + Bz () = (1.71) a,s (z)cl (1 cos ) + a,s (z)l a () + R (z)R () ds, l = где по определению Bz a,s (z)cl (z,, ), (1.72) Bz () a,s (z)l a () + R (z)R () ds.

l (1.73) В задаче для диффузного излучения d = 0, + t (z) = Bz () + Bz (0 ), = R + R dz 0 H (1.74) источник (z) Bz (0 ) = S exp l a,s (z)a (0 ) + R (z)R (0 ) = Bz 0 + Bz (0 ) (1.75) содержит сингулярную компоненту с -функцией B0, B0, (1 cos 0 ) Bz 0, (1.76) (z) B0, (z, 0 ) S exp a,s (z)cl, (1.77) и гладкую составляющую (z) B0 Bz (0 ) = S exp a,s (z)l a (0 ) + R (z)R (0 ).

l (1.78) Выделяем полное однократное рассеяние 1 :

(1.79) (z,, ) = 1 (z,, ) + (z,, ), которое находится из задачи (Этап 1):

d = 0, (1.80) + t (z) 1 = Bz (0 ), 1 1 = R 0 H dz 1.1. Скалярная плоская задача l с полной индикатрисой аэрозольного рассеяния a (). В силу (1.75) можно разложить 1 на две компоненты:

1, = 1, (1 cos 0 ), (1.81) 1 = 1, + 1, являющиеся решениями следующих задач (Этап 2 и Этап 3):

d 1 = 0, (1.82) + t (z) 1 = Bz (0 ), 1 1 = R0, dz 0 H d1, = 0, = 0. (1.83) 0 + t (z) 1, = B0, (z, 0 ), 1, 1, dz 0 H Разложение (1.81) используется для расчета «сглаженного» интеграла Bz (1 ) = Bz (1 ) + Bz (1, ) в задаче (1.87). В задаче d = 0, = R + R1 (1.84) + t (z) = Bz () + Bz (1 ), 0 H dz проведем преобразование коэффициента экстинкции, переходя к представле нию (1.71) для Bz ():

d + t (z) a,s (z) cl = Bz () + Bz (1 ), dz (1.85) + R1.

= 0, = R 0 H Выделяем второе рассеяние:

(1.86) (z,, ) = 2 (z,, ) + 3 (z,, ) как решение задачи (Этап 4):

d 2 + t (z) a,s (z) cl 2 = Bz (1 ), dz (1.87) 2 = 0, 2 = R1, 0 H которое рассчитывается при интеграле столкновений Bz (1 ) = Bz 1 + Bz (1 ) с полной функцией 1 – решением задачи (1.80).

– Учитывая -анизотропию в аэрозольной индикатрисе рассеяния (1.54) и в однократном приближении (1.81), распишем подробнее источник:

Bz (1 ) = Bz 1 + Bz (1 ) = (1.88) = Bz 1, + Bz 1 + Bz (1, ) + Bz (1 ) = B1, + B2, 58 Глава 1. Одномерные плоские задачи где с -анизотропией B1, Bz 1, = a,s (z) cl 1, = B1, (1 cos 0 ), (1.89) B1, (z, 0 ) = a,s (z)cl 1, (z, 0 ), гладкая составляющая B2 Bz 1 + Bz (1, ) + Bz (1 ), (1.90) Bz 1 = a,s (z)cl 1 (z,, ), (1.91) Bz (1, ) = 1, a,s (z)l a (0 ) + R (z)R (0 ).

l (1.92) Исходя из представления (1.88), второе приближение можно записать в виде суперпозиции 2, = 2, (1 cos 0 ), (1.93) 2 (z,, ) = 2 (z,, ) + 2,, члены которой находятся из следующих задач (Этап 5 и Этап 6):

d 2 + t (z) a,s (z) cl 2 = B2, dz (1.94) = 0, 2 = R1, 0 H 0 d2, + t (z) a,s (z) cl 2, = B1, (z, 0 ), dz (1.95) = 0, = 0.

2, 2, 0 H Разложение (1.93) используется для расчета источника – сглаженного – интеграла столкновений (1.96) Bz (2 ) = Bz (2 ) + Bz (2, ), Bz (2, ) = 2, a,s (z)l a (0 ) + R (z)R (0 ), l (1.97) в задаче для диффузного излучения, начиная с третьей кратности рассеяния (Этап 7):

d 3 + t (z) a,s (z) cl 3 = Bz (3 ) + Bz (2 ), dz (1.98) 3 = 0, 3 = R3 + R2, 0 H которая решается итерационно (Этап 8) по схеме (N 3 – номер итерации):

– N d3 + t (z) a,s (z) cl N = Bz (N 1 ) + Bz (2 ), 3 dz (1.99) N = RN 1 + R2.

N = 0, 3 3 0 H 1.1. Скалярная плоская задача При N = 3 полагаем Bz (N 1 ) = 0, RN 1 = 0. (1.100) 3 Аналогичная последовательность вычислений вводится при решении за дачи (1.49), когда оптическая толщина является пространственной координа той по высоте слоя. Выпишем основные формулировки задач и представления без особых комментариев.

Задача – аналог (1.69), имеющая решение в том же виде (1.70):

– d = S (s s0 ), + 0 = 0, = 0. (1.101) 0 d 0 H С учетом (1.54) из (1.52) находим:

B () S = B + B () = a,s ( )cl (1 cos ) + a,s ( )l a () + R ( )R () ds, (1.102) l = B a,s ( )cl (,, ), (1.103) B () a,s ( )l a () + R ( )R () ds.

l (1.104) Задача – аналог (1.74):

– d = 0, (1.105) + = B () + B (0 ), = R + R d 0 H содержит источник (1.106) B (0 ) = B 0 + B (0 ) = B0, + B0, B0, B0, (1 cos 0 ) B 0, (1.107) (z) B0, (, 0 ) = S exp a,s ( )cl, (1.108) (z) B0 B (0 ) = S exp a,s ( )l a (0 ) + R ( )R (0 ). (1.109) l Задача (Этап 1) — аналог (1.80):

d = 0, (1.110) + 1 = B (0 ), 1 1 = R d 0 H разбивается на две задачи (Этап 2 и Этап 3) — аналоги (1.82) и (1.83):

d 1 = 0, (1.111) + 1 = B (0 ), 1 1 = R0, 0 H d d1, = 0, = 0. (1.112) 0 + 1, = B0, (, 0 ), 1, 1, 0 H d 60 Глава 1. Одномерные плоские задачи Задача – аналог (1.84) – d =0, (1.113) + = B () + B (1 ), = R + R1, d 0 H преобразуется к аналогу (1.85) d + 1 a,s ( ) cl = B () + B (1 ), d (1.114) = 0, = R + R1.

0 H Задача для второго рассеяния (Этап 4) – аналог (1.87) – d 2 + 1 a,s ( ) cl 2 = B (1 ), d (1.115) 2 = 0, 2 = R1, 0 H разбивается на две задачи (Этап 5 и Этап 6) – аналоги (1.94) и (1.95) – d 2 + 1 a,s ( ) cl 2 = B2, d (1.116) = 0, 2 = R1, 0 H 0 d2, + 1 a,s ( ) cl 2, = B1, (, 0 ), d (1.117) = 0, = 2, 2, 0 H с источниками – аналогами (1.89) и (1.90).

– Аналогом задач (1.98) и (1.99) (Этап 7 и 8) является задача N d3 + 1 a,s ( ) cl N = B (N 1 ) + B (2 ), 3 d (1.118) N = RN 1 + R2.

N = 0, 3 3 0 H Для реализации численных алгоритмов описанные выше краевые задачи можно формализовать в виде двух задач:

а) для задач (1.44), (1.49), (1.69), (1.74), (1.80), (1.82), (1.83), (1.84), (1.101), (1.105), (1.110), (1.111), (1.112), (1.113) d = 0, (1.119) + tot (x) = F (x,, ), = f (, ), 0 H dx если x = z, t (z), (1.120) tot (x) = 1, если x = ;

1.1. Скалярная плоская задача б) для задач (1.85), (1.87), (1.94), (1.95), (1.98), (1.99), (1.114), (1.115), (1.116), (1.117), (1.118) d = 0, (1.121) + tot (x) = F (x,, ), = f (, ), dx 0 H t (z) a,s (z) cl, если x = z, (1.122) tot (x) = 1 a,s ( ) cl, если x =.

Расчет интегралов столкновений также можно формализовать, вводя следующие записи:

в) для интегралов (1.51), (1.52) B() l (1.123) s,a (x) a () + s,R (x) R () ds, г) для интегралов (1.73), (1.78), (1.92), (1.97), (1.104), (1.109) B() l (1.124) s,a(x)a () + s,R (x)R () ds, д) для интегралов (1.72), (1.76), (1.77), (1.91), (1.103), (1.107), (1.108) (1.125) B() = s, (x)(x,, ), где если x = z, если x = z, a,s (z), R (z), s,a (x) = s,R (x) = если x = ;

если x = ;

a,s ( ), R ( ), a,s (z)l, a,s (z)cl, если x = z, если x = z, s,a (x) = s, (x) = ( )l ( )cl, если x = ;

, если x =.

a,s a,s (1.126) Коэффициенты (1.120), (1.122) и (1.126) послойно аппроксимируются кусочно-постоянными значениями:

s,,k для xk tot,k, tot,k, x xk+1;

+ + + для xk x xk+1 при 0;

s,a,k, s,R,k, s,a,k для xk x xk+1 при 0.

s,a,k, s,R,k, s,a,k 1.1.6.3. Формулы расчета 1, -анизотропии. Выпишем аналитические ре шения задачи (1.83) (через z) с источником Bz 0 = B0, (z, 0 )(1 cos 0 ), (z) B0, (z, 0 ) = a,s (z)cl S exp.

62 Глава 1. Одномерные плоские задачи Для = 0 0 интегрированием по характеристике получаем z z 1 B0, (t) exp t (t ) dt dt = 1, (z, 0 ) = 0 t z S (z) exp a,s (t)cl (t) dt. (1.127) = 0 В явном виде для k = 1 K k S (zk+1 ) a,s,n c+ zn.

+ exp (1.128) 1, (zk+1 ) = n 0 n= Получим представление через рекуррентную формулу:

1, (z1 ) = 0, zk+1 zk+ 1 B0, (t) exp t (t ) dt dt = 1, (zk+1 ) = 0 t k S (zk+1 ) a,s,k c+ zk + = 1, (zk ) exp exp. (1.129) + k 0 0 Решение задачи (1.112) (через ) с источником B 0 = B0, (, 0 )(1 cos 0 ), B0, (, 0 ) = a,s ( )cl S exp, для = 0 0 находим интегрированием по характеристике:

t 1 S dt = exp B0, (t) exp a,s (t) cl (t) dt.

1, ( ) = 0 0 0 0 0 (1.130) В явном виде для k = 1 K k+ S 1, (k+1 ) = exp k+1 a,s (t)cl (t) dt = 0 k S a,s,n c+ n.

+ exp k+1 (1.131) = n 0 n= 1.1. Скалярная плоская задача Через рекуррентную формулу 1, (1 ) = 0, k+1 k+ 1 B0, (t) exp d dt = 1, (k+1 ) = 0 t k S a,s,k c+ k.

+ exp k+ = 1, (k ) exp (1.132) + k 0 0 Так как + a,s,k + + a,s,k c+ k = + c zk = a,s,k c+ zk, + + k t,k k t,k k то выражения 1, (zk ) через z и 1, (k ) через совпадают.

Введем формализованную запись, принимая, если через оптическую толщину, x= если через геометрические координаты.

z, Общая запись z (z) = S (z) 1, exp a,s (t)cl (t) dt, если x = z, 0 1, (x) = S exp a,s (t)cl (t) dt, 1, ( ) = если x =, 0 т. е. x S (x) 1, (x) = exp a,s (t)cl (t) dt. (1.133) 0 В дискретном виде для + = 0 0 имеют место явные представления:

1, (x1 ) = 0, k S s,a,k c+ xk, + exp k+ 1, (xk+1 ) = 1, (xk ) exp + k 0 0 + a,s,k, если x = z, zk, если x = z, + s,a,k = xk = + a,s,k, если x = ;

k, если x = ;

k+1 k, если x =, k = (zk+1 ) (zk ), если x = z ;

если x =, k+1, k+1 = если x = z.

(zk+1 ), 64 Глава 1. Одномерные плоские задачи 1.1.6.4. Формулы расчета 2,. Решение задачи (1.95) (через z) с источ ником B1, (z, 0 ) = a,s (z)cl (z)1, (z, 0 ) определяем интегрированием по характеристике:

z z 1 B1, (t) exp t (z ) a,s (z )cl dz dt.

(1.134) 2, (z) = 0 t Найдем выражение через рекуррентную формулу:

zk+1 zk+ 1 B1, (z) exp t (z ) cl (z )a,s (z ) dz dz = 2, (zk+1 ) = 0 z c+ + k a,s,k zk k = 2, (zk ) exp (1.135) + 2, (zk+1 ), zk+ a,s (z)cl (z)1, (z) 2, (zk+1 ) = zk zk+ exp t (z ) cl (z )a,s (z ) dz dz = I1 + I2.

z Найдем явное выражение для интеграла zk+ a,s,k c+ (zk+1 z) + a,s,k c+ zk + k k exp dz = exp a,s,k c+ + 0 k zk и получим k a,s,k c+ zk + S (zk+1 ) a,s,k c+ zk.

+ I1 = exp k exp k 0 0 n= Используем известный интеграл ax 1 ax xeax dx = e a при вычислении второго слагаемого a,s,k c+ zk + a,s,k c+ zk + (zk+1 ) k k I2 = S exp 1+ exp.

0 0 1.1. Скалярная плоская задача Так что расчеты можно проводить по явному выражению в рекуррентной форме:

+ c+ zk k + a,s,k k 2, (zk+1 ) = 2, (zk ) exp + 0 a,s,k c+ zk + S (zk+1 ) k exp 1 + exp 0 0 k (zk+1 ) a,s,k c+ zk + S exp + k n= a,s,k c+ zk + a,s,k c+ zk + k k exp 1+ (1.136).

0 Решение задачи (1.44) (через ) с источником B1, ( ) = a,s ( )cl ( )1, ( ) определяем через интеграл по характеристике 1 B1, (t) exp 1 wa,s ( )cl ( ) d dt 2, ( ) = 0 t аналогично случаю через z.

В формализованном виде явная формула расчета 2, (x) как решения задачи d2, + (x) (x) cl (x) = B (x), 0 1, tot 2, dx 2, = 0, = 0.

2, 0 H с источником B1, (x) = (x)cl (x)1, (x) находится интегрированием по характеристике:

x x 1 B1, (t) exp tot (x ) cl dx dt = 2, (x) = 0 t x x 1 (t)cl (t)1, (t) exp tot cl dt dt, = 0 t 66 Глава 1. Одномерные плоские задачи где x S (x) 1, (x) = exp (x )cl (x ) dx, 0 если x = z, t (z), tot (x) = 1, если x = ;

если x = z, a,s (z), (x) = s,a = если x =.

a,s ( ), Так что t x S (x) (t)cl (t) (x )cl (x ) dx exp 2, (x) = 2 0 x exp (t )cl (t ) dt dt. (1.137) t В рекуррентном виде для = 0 0 имеем: 2, (x1 ) = 0;

+ c+ x k + kk k 2, (xk+1 ) = 2, (xk ) exp (1.138) + 2, (xk+1 ) 0 и аналогично предыдущему получаем S (xk+1 ) 2, (xk+1 ) = 2 exp x k + c+ xk (x )cl (x ) dx exp kk 1 + + c+ xk + c+ xk (xk+1 ) kk kk + S exp 1+ exp.

0 0 1.1.6.5. Расчет однократного рассеяния. Решение задачи (1.80) находим интегрированием по характеристике для 0, = 0 :

z z 1 Bz (0 ) exp t (u) du dt + 1 (z, 0,, ) = + + t 0 z + f0 exp t (u) du. (1.139) + 1.1. Скалярная плоская задача Конечно-разностный аналог k 1 (zk+1 ) = 1 (zk ) exp (1.140) + k, + где S 0 a,s,k l ( ) + R,k R (0 ) k = t,k a 0 + t,k (zk+1 ) (zk ) + +k exp exp.

0 0 Аналогично для = k S 0 (zk ) 1 (zk ) = 1 (zk+1 ) exp exp + | | 0 + | | 1 1 a,s,k l a (0 ) + R,k R (0 ).

1 exp k | | 0 t,k t,k Для = 0 k 1 (zk+1 ) = 1 (zk ) exp + S (zk+1 ) exp l + zk a,s,k (t)a (0 ) + R,k (t)R (0 ).

0 При решении задачи (1.82) делаются следующие замены:

a,s (z) на a,s (z)l ;

l l a (0 ) на a (0 ).

При решении задачи (1.110) заменяются a,s (z) на a,s ( ), R (z) на R ( ), t (z) = 1.

При решении задачи (1.111) заменяются a,s (z) на a,s ( )l, l l R (z) на R ( ), a (0 ) на a (0 ), t (z) = 1.

1.1.6.6. О представлении рэлеевской индикатрисы рассеяния и индика трисы Хеньи-Гринстейна. В соответствии с условием нормировки R (cos ) d cos = 1 или 2 R () d = рэлеевская индикатриса рассеяния записывается в виде (1 + cos2 ), R () = где – угол рассеяния. В разложении по полиномам Лежандра – R () = n Pn (cos ) n= 68 Глава 1. Одномерные плоские задачи коэффициенты 2n + n = R (cos )Pn (cos ) d cos, 4 где P2 = (32 1)/2, P0 = 1, P1 =, [(2n 1)Pn1 (n 1)Pn2 ], Pn = n Pn ()Pm () d = nm.

2n + При этом 13 (1 + 2 ) d = = ;

4 2 16 (1 + 2 ) d = 0 ;

= 4 2 53 1 (32 1) (1 + 2 ) d = =, 4 2 16 2 2 т. е. 0 = 1, 1 = 0, 2 = 0.5, 1 (1 + 2 ).

R () = [ + 2 P2 ()] = 4 0 G. W. Kattawar предложил алгоритм аналитического представления инди катрис рассеяния с тремя параметрами HG (, g1, g2, ) = HG (, g1 ) + (1 )HG (, g2 ), в которое входят с весами и (1 ) две индикатрисы Хеньи-Гринстейна 1 g HG (, g) =, 4 (1 2g + g 2 )3/ имеющие нормировку 2 1 HG (, g) d d = 2 HG (, g) d = 1, 0 1 – косинус угла рассеяния, с параметром асимметрии – 2 g = P1 () = = HG (, g) d d, 0 1.1. Скалярная плоская задача равным g1 и g2. Дифференцируя по разложение 2 1/ (1 2g + g ) |g| 1, gn Pn (), = n= получаем представление индикатрисы Хеньи-Гринстейна через полиномы Лежандра gn (2n + 1)Pn (), HG (, g) = n= 2 Pn ()HG (, g) d d = gn.

Pn () = 0 Трехпараметрическое разложение индикатрисы рассеяния [g1 + (1 ) g2 ] (2n + 1)Pn () n n HG (, g1, g2, ) = n= содержит три неизвестных параметра, g1, g2, которые можно определить через значения g, h, t трех моментов Pn () = g1 + (1 ) g2, n n n = 1, 2, 3.

Разрешая систему уравнений g1 + (1 ) g2 = g, g2 + (1 ) g2 = h, g1 + (1 ) g2 = t, находим явные выражения параметров 1/ t h g (h g t)2 4(h g 2 ) (t g h2 ) g2 =, 2(h g 2 ) g g2 h g g (1.141) g1 =, =.

g2 g g1 g Если использовать только два момента, то из системы уравнений g1 + (1 ) g2 = g, g1 + (1 ) g2 = h 2 находим значения 1/ (h g 2 ) g (1 )g g2 = g, g1 = 1 со свободным параметром. Можно ввести новый параметр = /(1 ), с помощью которого можно варьировать погрешности в области глории и т. п., где возникают наибольшие отклонения от исходной индикатрисы рассеяния.

70 Глава 1. Одномерные плоские задачи Трехпараметрическое представление позволяет описать рэлеевскую инди катрису рассеяния с помощью двух индикатрис рассеяния Хеньи-Гринстейна.

Однопараметрической индикатрисой аппроксимировать R () нельзя, так как параметр асимметрии для индикатрисы рассеяния Хеньи-Гринстейна (1 g 2 ) d = = (1 2g + g 2 )3/ 2 n g (2n + 1) = Pn ()P1 () d = g, n=0 а для рэлеевской индикатрисы 1 2 (1 + 2 ) d = 0.

= [0 + 2 P2 ()] d = 4 1 Вычислим три момента Pn () для рэлеевской индикатрисы:

g = = P1 () = 0 ;

[0 + 2 P2 ()] P3 () d = 0 ;

t = P3 () = h = P2 () = 2 P2 ()R () d = 0.1.

Полагая g = t = 0 в формулах для параметров (1.141), находим 0.1 0. g2 = 0.1 0.3162, g1 = = 0.1 0.3162, = = 0.5.

g2 2 0. Так что с g1 = g2 = 0.1 получаем 1 g1 1 g 2 R () = + = (1 2 g1 + g1 )3/ (1 + 2g1 + g1 )3/ 4 11 0.9 0. (1.142) = + (1.1 2 0.1 ) 4 2 3/ (1.1 + 2 0.1)3/ или 1 2n + 1 n g1 [1 + (1)n ] Pn () = R () = 4 n= 1 {1 + P2 () + 0.09 P4 () +...}, (4m + 1) (0.1)m P2m () = = 4 m=0 (1.143) 1.1. Скалярная плоская задача т. е. в разложении присутствуют только четные полиномы Лежандра. При N = 5 значения индикатрис рассеяния, рассчитанные по формулам (1.142) и (1.143), практически совпадают, но обе аппроксимации, использующие равенство трех моментов, заметно отличаются от исходной рэлеевской индикатрисы рассеяния.

1.1.6.7. Двухчленная индикатриса рассеяния. В модели без разделения аэрозольного и рэлеевского рассеяния можно включить приближение Эд дингтона с индикатрисой рассеяния (cos ) E (cos ) = (1 + 3g cos ) и приближение -Эддингтона с индикатрисой рассеяния (cos ) E (cos ) = [2f (1 cos ) + (1 f )(1 + 3g cos )].

«Средний косинус» определяется через первый момент () d g = = P1 () =, () d а весовой коэффициент – через второй момент – ()P2 () d f = P2 () =, () d (32 1)/2 – нормированный где P2 () = – полином Лежандра второго порядка.

В разложении индикатрисы рассеяния (cos ) = n Pn (cos ) n= при нормировке 2 1 1 () d = (cos ) d d = 4 0 1 коэффициенты при полиномах Лежандра 1 ()Pn () d 2n+1 n = ()Pn () d = (2n+1) = (2n+1) Pn ().

2 () d 72 Глава 1. Одномерные плоские задачи Так что 0 = 1, 1 = 3 P1 () = 3g, т. е. в приближении Эддингтона в разложении индикатрисы рассеяния берется только два первых числа.

Для приближения -Эддингтона «средний косинус»

[2f (1 ) + (1 f )(1 + 3g )] d = f + g (1 f ).

E = [2f (1 ) + (1 f )(1 + 3g )] d При условии равенства первых моментов g = E «средний косинус»

полагается равным gf g=.

1f Весовой коэффициент находится из условия равенства вторых моментов P2 () = P2,E (), P2 () [2f (1 ) + (1 f )(1 + 3g )] d P2,E () = =f.

[2f (1 ) + (1 f )(1 + 3g )] d Для приближения Эддингтона суммарная аэрозольно-рэлеевская индика триса рассеяния ( ) ( ) (, ) = a,s a (, ) + R R () s ( ) s ( ) аппроксимируется двухчленным разложением по полиномам Лежандра с коэффициентами 0 = 1, 1 ( ) = 3g( ).

В случае -Эддингтон приближения с учетом того, что B() = s ( ) E (, cos ) ds = s ( )f ( )(,, ) + BE (), (1 f ( )) s ( ) (1 + 3g ( ) cos ) ds, BE () = проводится преобразование коэффициентов уравнения переноса:

d + t ( ) = s ( ) (, ) ds, d t ( ) = 1 s ( ) f ( ), s ( ) = (1 f ( )) s ( ), 1 ( ) = 3g ( ), 0 = 1.

(, ) = ( + 1 ( ) cos ), 4 1.2. Поляризационные задачи В задаче с однородным слоем можно перейти к новой оптической толщине = (1 s f ) и новому коэффициенту альбедо акта рассеяния (1 f ) s s = 1 s f и решать уравнение переноса в наиболее распространенном виде d + (,, ) = s () ds.

d Такое преобразование называют «методом подобия».

§ 1.2. Поляризационные задачи Решение задач теории переноса поляризованного излучения, как правило, предполагает применение вектор-параметра Стокса и матричного формализма.

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

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

Потребность анализа процесса распространения света в реальной атмо сфере и увеличение производительности ЭВМ обусловили развитие численных методов для решения векторного уравнения переноса. В основном развивались два метода. Один из них – метод Монте-Карло для плоской и сферической – геометрии, которым рассчитываются интенсивность и степень поляризации излучения. Больших успехов развитие этого метода достигло в ВЦ АН СССР.

Второй – это матричный метод, основанный на методе удвоения в сочетании – с методом сложения слоев. Метод удвоения требует разложения решения в ряд по азимутальным гармоникам, что порождает проблему корректного суммирования ряда, особенно острую из-за знакопеременности решений;

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

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

Успешная реализация этого метода стала возможной благодаря полному учету 74 Глава 1. Одномерные плоские задачи свойств решаемой задачи и тщательной проработке и оптимизации алгоритма.

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

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

Однако различные авторы в понятие параметров Стокса вкладывают несколько различный смысл. Под общим названием «параметры Стокса»

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

Последовательное и полное изложение способов описания поляризованного света применительно к задачам эллипсометрии дали Р. Аззам и Н. Башара.

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

Коротко введем необходимые определения. Будем рассматривать только монохроматические и квазимонохроматические электромагнитные волны. Луч света понимается в смысле геометрической оптики. Для полного описания эллиптической поляризации однородной плоской поперечно-электрической бегущей волны в плоскости, перпендикулярной направлению ее распростра нения s, требуются четыре величины. Основой для получения таких величин является разложение волны на сумму двух простых гармонических колебаний в декартовом базисе e1, e2 (тройка e1, e2, s правая, см. рис. 1.3). При этом получаются две амплитуды E1 и E2 и две фазы 1 и 2 простых гармонических колебаний:

2s 2s E(s, t) = E1 cos t e1 + E2 cos t e2. (1.144) + 1 + Здесь t – время, – круговая частота, s – координата вдоль направления – – – распространения, – длина волны.

– Эллипс поляризации лежит в плоскости волнового фронта. Вместо четырех величин E1, E2, 1, 2 можно использовать параметры, характеризующие эллипс в его плоскости:

1.2. Поляризационные задачи e a Et= e A b e s e Рис. 1.3. Правая система коорди- Рис. 1.4. Четыре параметра (,, A, ), определя нат e1, e2, s, связанная с направ- ющие эллипс поляризации в его плоскости лением распространения волны s – азимут большой оси a, отсчитываемый от фиксированного направления – оси X, ;

2 b – эллиптичность = ±, – 1;

a – угол эллиптичности : = tg, –.

4 Поляризация называется правой, если эллипс обходится по часовой стрелке при наблюдении против направления распространения пучка света;

полная амплитуда A = (a2 + b2 )1/2, где a и b – большая и малая полуоси эллипса;

– абсолютная фаза = 2 1, 2. Определение параметров поляриза ции показано на рис. 1.4. Изображена правая поляризация. Абсолютная фаза определяется начальным положением вектора электрического поля Et=0 и вспомогательной (пунктирной) окружностью.

Для описания взаимодействия между полностью поляризованной волной и оптическими элементами, изменяющими поляризацию этой волны, а также для описания недеполяризующих процессов применяются двухкомпонентные комплекснозначные декартовые векторы Джонса (комплексные амплитуды) и матрицы Джонса, преобразующие эти комплексные амплитуды. Вектором Джонса (столбцом комплексных амплитуд или столбцом Максвелла) моно хроматической однородной поперечно-электрической волны называют вектор E1 exp(i1 ) Ex E = {E1, E2 }. (1.145) Exy = = = Ey E E2 exp(i2 ) Запись вектор-столбца горизонтально в фигурных скобках также использу ется. Временная и пространственная информация о волне восстанавливается 76 Глава 1. Одномерные плоские задачи по формуле 2s E(s, t) = Re Exy exp i t (1.146).

Двумерные комплекснозначные векторы образуют комплексное векторное пространство. В этом пространстве базисные векторы можно выбирать произвольным образом, например, e1 = {1, 0}, e2 = {0, 1} – декартовы базисные векторы, – 1 e1 = {1, i}, er = {i, 1} – круговые базисные векторы.

– 2 Элементы вектора Джонса (1.145) – коэффициенты разложения волны в – декартовом базисе. В общем случае в качестве базисных векторов можно выбрать две различные эллиптические, не обязательно ортогональные, по ляризации eu и ev. Декартов вектор Джонса Exy можно представить как линейную суперпозицию базисных комплекснозначных векторов eu и ev :

(1.147) Exy = Eu eu + Ev ev.

Коэффициенты разложения образуют обобщенный вектор Джонса Euv = {Eu, Ev }.

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

Квазимонохроматическая плоская поперечно-электрическая волна с поло сой частот и центральной частотой 0 записывается в виде 2s 2s E(s, t) = E1 (t) cos 0 t e1 + E2 (t) cos 0 t e2.

+ 1 (t) + 2 (t) (1.148) В случае полностью поляризованной квазимонохроматической волны выполняются соотношения E2 (t) = 2 (t) 1 (t) = const.

= const, E1 (t) При этом амплитуда эллипса поляризации может зависеть от времени.

Различные возможные состояния плоской поперечно-электрической све товой волны представляются вектором, составленным из четырех действи тельных величин I, Q, U, V, называемых параметрами Стокса:

(s) = {I, Q, U, V }.

Каждая из этих величин имеет размерность интенсивности.

1.2. Поляризационные задачи Параметры Стокса следующим образом определяются через амплитуды и фазы простых гармонических колебаний разложения квазимонохроматической волны (1.148):

Q = E1 (t) E2 (t), 2 2 2 (1.149) I = E1 (t) + E2 (t), U = 2 E1 (t)E2 (t) cos[2 (t) 1 (t)], V = 2 E1 (t)E2 (t) sin[2 (t)1 (t)].

Угловые скобки означают усреднение по времени:

T v v(t) dt, T где T – интервал времени, достаточно большой, чтобы интеграл по времени – не зависел от выбора T. Параметры Стокса для квазимонохроматической волны имеют статистический смысл.

Выражение параметров Стокса для частично поляризованной квазимо нохроматической волны через параметры эллипса поляризации полностью поляризованной компоненты (азимут, угол эллиптичности и степень поляризации P ) имеет следующий вид:

Q = I0 P cos 2 cos 2, I = I0, (1.150) U = I0 P sin 2 cos 2, V = I0 P sin 2.

Здесь I0 A2 (t) – полная интенсивность волны.

– Степень поляризации P определяется как отношение интенсивности полностью поляризованной компоненты к общей интенсивности волны:

(Q2 + U 2 + V 2 )1/ 0 1. (1.151) P=, P I Степень поляризации изменяется от нуля для неполяризованного света до единицы для полностью поляризованного света. Параметры Стокса полностью поляризованной волны удовлетворяют соотношению I 2 = Q2 + U 2 + V 2.

Параметры Стокса для частично поляризованного и неполяризованного света удовлетворяют неравенству (соотношению «конуса») I 2 Q2 + U 2 + V 2.

Частично поляризованный свет можно разложить на неполяризованную n и поляризованную p компоненты:

= n + p, (1.152) n = {[I (Q2 + U 2 + V 2 )1/2 ], 0, 0, 0}, p = {(Q2 + U 2 + V 2 )1/2, Q, U, V }.

78 Глава 1. Одномерные плоские задачи В подпространстве Стокса (Q, U, V ) состояние поляризации квазимоно хроматической волны представляется точкой с полярными координатами (P, / 2 2, 2). Единичная сфера (P = 1) в подпространстве Стокса является сферой Пуанкаре.

Переход от параметров Стокса к параметрам эллипса поляризации – – азимуту и углу эллиптичности – происходит по формулам – 1 U V = arctg, = arcsin.

(Q + U 2 + V 2 )1/ 2 Q Знак параметра Q совпадает со знаком. Для определения угла по главному значению arctan 0 используется формула (1 sign Q) sign U.

= 0 + При переходе к параметрам Стокса теряется информация об абсолютной фазе волны.

Параметры Стокса взаимно однозначно выражаются через элементы матрицы когерентности E1 E1 E1 E2 J11 J (1.153) J= =.

E2 E1 E2 E2 J21 J Здесь E1 и E2 – элементы вектора Джонса Exy. Из элементов матрицы – когерентности формируется декартов вектор когерентности J = {J11, J12, J21, J22 }. (1.154) Определим параметры Стокса через элементы матрицы когерентности:

Q = J11 J22, V = i(J12 J21 ), (1.155) I = J11 + J12, U = J12 + J21, или в матричном виде I J Q J 1 0 0 1 = = = T J. (1.156) U 0 1 1 0 J21 0 i i V J При введении параметров Стокса в литературе иногда приводится только одно какое-либо из определений типа (1.149), (1.150), (1.156). Требуется безусловное согласование таких определений.

Определение матрицы когерентности обобщается путем подстановки в него вместо E1, E2 элементов Eu, Ev обобщенного вектора Джонса. В результате получается обобщенная матрица когерентности и, как следствие формул (1.156), обобщенный вектор Стокса uv. Г. В. Розенберг называет параметры Стокса, полученные в ортонормированном базисе eu, ev и характеризующиеся азимутом поляризации и углом эллиптичности, параметрами Стокса в 1.2. Поляризационные задачи представлении (, ). Таким образом, определения (1.149), (1.150), (1.156) есть определения параметров Стокса в декартовом представлении.

Помимо указанных формул, для точного определения параметров необ ходимо зафиксировать:

а) систему координат, связанную с волной;

б) вид разложения волны на простые гармонические колебания;

в) определение комплексных амплитуд;

г) определение параметров эллипса поляризации.

Естественно пункты а–г включать составной частью в понятие пред ставления параметров Стокса. Параметры Стокса в этом «стандартном»

представлении обозначаем (s). Для обозначения представлений вводятся верхние индексы в круглых скобках.

Нами проводилось сопоставление разных определений параметров Стокса в основном с помощью определений типа (1.149).

Рассмотрим преобразование параметров Стокса = (s) при повороте системы координат. Поворот координатных осей на угол 0 в направлении против часовой стрелки (рис. 1.5), если смотреть навстречу направлению распространения волны (направление положительного отсчета угла ), опи сывается матрицей L:

1 0 0 0 cos 2 ± sin 2 = L(), L(±) =, (1.157) 0 ± sin 2 cos 2 0 0 0 где соответствует вектору Стокса до поворота, – после поворота осей – на угол. Параметры I и V инвариантны относительно поворота. Матрица поворота обладает следующими свойствами:

3) L( ) = L(), 1) L(1 ) L(2 ) = L(1 + 2 ), L1 () 4) L( ± 2) = L().

2) = L(), Допустим, имеются два представления вектора Стокса: (1) и (2) и переход от первого представления ко второму осуществляется с помощью матрицы M : (2) = M (1). Если параметры (1) изменяются посредством некоторого процесса (поворота или замены системы координат, действия оптического прибора, взаимодействия с веществом), описываемого матрицей 1 : (1) = 1 (1), то параметры Стокса в представлении (2) будут изменяться с помощью матрицы 2 :

(2) = 2 (2), 2 M 1 M 1.

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

80 Глава 1. Одномерные плоские задачи e e e s 0 e Рис. 1.5. Схема поворота координатных осей на угол 0. Вектор распространения волны s направлен к наблюдателю (перпендикулярно к плоскости листа) Во многих работах встречаются системы параметров, для которых матрица T перехода от вектора когерентности к параметрам Стокса иная, чем матрица T в определении (1.156). Будем считать, что в таком случае параметры Стокса определены в представлении, заданном матрицей T при соответствующем базисном разложении.

Широко используется представление параметров Стокса, введенное С. Чандрасекаром. В этом случае имеет место разложение волны в декартовом базисе на простые гармонические колебания (1.144) и определения а–г совпа дают со «стандартными». Параметры Стокса в представлении С. Чандрасекара обозначаем (ch) = {I, I, U, V }, I Q I +Q I I1 = E1 = I I2 = E2 = 2,, 2 матрица поворота sin cos2 () sin () cos2 () L(ch) () = sin () sin 2 sin sin 2 cos 0 0 0 Переход от вектора когерентности к параметрам Стокса в представлении (ch) определяется матрицей T (ch) :

1 0 0 0 0 (ch) = J = T (ch) J.

0 1 0 i i 1.2. Поляризационные задачи Переход от представления параметров Стокса (s) к представлению (ch) задается матрицей Ms=ch:

1 1 0 0 1 1 1 1 0 1 1 1 1 Ms=ch = ;

M 1 = Mch=s.

s=ch = 2 2 0 0 2 0 0 0 002 00 0 Другое часто используемое представление:

I V I +V (c) = {Ir, Il, Q, U }, Ir =, Il =, 2 матрица поворота 1 0 0 0 1 0 L(c) () =.

0 0 cos 2 sin 0 0 sin 2 cos Переход от декартова вектора когерентности к параметрам Стокса в представлении (c) определяется матрицей T (c) :

1 i i 1 ii (c) = (c) 1 0 0 1 J = T J.

0 1 1 параметров Стокса (s) к представлению (c) Переход от представления задается матрицей Ms=c :

1 0 0 1 1 0 1 0 0 1 0 0 1 Ms=c = ;

= = Mc=s.

Ms=c 2 0 2 0 0 0 0 0 1 0 002 Представление вектора Стокса с комплекснозначными первым и четвертым параметрами конструируется на основе кругового вектора Джонса. Матрица поворота на угол имеет диагональный вид, так как базисные векторы являются собственными векторами оператора вращения вокруг направления распространения луча:

Q iU 1 i i 1 I V 1 i i (cpn) = = J = T (cpn)J, 2 I + V 1 i i i Q + iU i матрица поворота L(cpn)() = diag{exp (2i), 1, 1, exp(2i)}.

82 Глава 1. Одномерные плоские задачи Переход от представления параметров Стокса (cpn) к представлению (s) задается матрицей Ms=cpn:

1 i i 1 0 1 1 1 1 i i 1 1 0 0 Ms=cpn = ;

Ms=cpn = i 0 0 i = Mcpn=s.

2 1 i i i 1 0 1 1 i J. Hovenier использует систему параметров, отличающуюся от предыдущей только перестановкой параметров:

I + V 1 i i 1 I V 1 i i = J = T (cph)J, (cph) = 2 Q + iU 1 i i Q iU 1 i i матрица поворота L(cph) () = diag{1, 1, exp (2i), exp(2i)}, матрицы перехода 0 1 i 0 0 11 0 0 1 1 Ms=cph = Ms=cph = = Mcph=s.

1 0 0 1 ;

i i 0 1 i 0 1 0 1.2.2. Матрица рассеяния. Фазовая матрица. Предположим, что луч света рассеивается отдельной частицей или малым элементом объема рассеивающего вещества в точке пространства C(r) и рассеяние происходит независимо на каждой частице и без изменения частоты. Реальный световой пучок (луч) всегда ограничен в пространстве и во времени. Определим плоскость рассеяния R, содержащую падающий и рассеянный лучи (рис. 1.6). Относительно этой плоскости R определим систему координат 1, 2, s, связанную с падающим ee 1, e2, s, связанную с рассеянным лучом. Векторы лучом, и систему координат e e1 и 1 лежат в плоскости рассеяния R. Обе тройки правые (плоскости e рассеяния и референции совпадают).

Параметры эллипса поляризации в плоскости e1, e2 представлены на рис. 1.7. Векторы Стокса, описывающие падающий и рассеянный лучи, в соответствующих базисах связаны матрицей рассеяния : =. Кон кретный вид матрицы рассеяния для одного и того же акта взаимодействия определяется в зависимости от представления векторов Стокса и и выбора базисов (систем координат) для падающего и рассеянного лучей.

1.2. Поляризационные задачи / e R Падающий свет s s R Ра s s сс ея нн 1 ый e e e 2 св e e ет e e Рис. 1.6. Схема определения плоскости Рис. 1.7. Параметры эллипса поляриза рассеяния R ции и плоскость рассеяния R В общем случае 4 4 – матрица содержит 16 элементов, которые – выражаются через элементы матрицы преобразования комплексных амплитуд.

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

Если среда макроскопически оптически изотропна, т. е. отсутствуют явления линейного и циркулярного двойного лучепреломления и дихроизма, то матрица рассеяния (z, s s) является функцией угла рассеяния :

(z, s s) + (z, s s) или (z,,, ) = (z, cos ), где s = {, } – направление излучения, падающего на элементарный объем, – s = {, } – направление излучения, рассеянного элементарным объемом. В – общем случае (z, cos ) как матрица рассеяния четвертого порядка содержит 16 действительных функций mn (z, cos ).

Условие нормировки элементов матрицы (z, cos ) определяется нормиров кой первого элемента 11 (z, cos ), который является индикатрисой рассеяния в обычном понимании и удовлетворяет условию:

2 11 (z,,, ) d = 1.

d Если среда обладает свойствами инвариантности рассеяния в акте отно сительно плоскости рассеяния, то на основе качественного рассмотрения можно получить, что в фиксированной точке слоя матрица рассеяния определяется всего 6 действительными функциями и является матрицей 84 Глава 1. Одномерные плоские задачи блочно-диагонального вида:

a ( ) b( ) 0 1s s b(s ) a2 (s ) 0 (s ) =, 0 0 a3 (s ) c(s ) c(s ) a4 (s ) 0 где 0 – угол рассеяния. Такие матрицы рассеяния наиболее – s характерны для класса решаемых атмосферно-оптических задач.

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


0 a1 (cos ) b(cos ) b(cos ) a (cos ) 0 Ми (cos ) =.

a3 (cos ) c(cos ) 0 c(cos ) a4 (cos ) 0 Для этих матриц выполняются условия физической и математической корректности:

|a2 | |a3 | |b| |a4 a3 |, |c| a3 |.

a1 0, a1, a1, a1, Частным случаем матрицы рассеяния Ми является рэлеевская матрица, в которой элемент c = 0, что означает отсутствие эллиптической поляризации в акте рассеяния:

cos2 + 1 cos2 1 0 3 cos 1 cos + 2 0.

Re (cos ) = 16 0 0 2 cos 0 0 0 2 cos Если изолированные рассеиватели среды не обладают плоскостью симмет рии, но распределение рассеивателей по ориентациям полностью хаотическое, то матрица рассеяния элемента объема такой среды определяется совокуп ностью из 10 действительных функций:

11 12 13 23 12 x =.

13 23 33 24 34 При использовании матриц рассеяния в численном моделировании необхо димо проверять главное свойство физической корректности: матрицы должны сохранять свойство конуса I 2 Q2 + U 2 + V 2 и переводить параметры Стокса в параметры Стокса. Класс физически корректных матриц рассеяния 1.2. Поляризационные задачи 0 = = z= Верх s = x 90 C(r) = {x, y, z} = y X 0 Y s Z = z=H = z Дно Рис. 1.8. Система координат для плоского слоя вкладывается в класс неотрицательных матриц. Пространство параметров Стокса (s) образует выпуклый конус в действительном четырехмерном пространстве R(4).

Для описания распространения поляризованного света в плоском слое вводим декартовую правую систему координат x, y, z (рис. 1.8), связанную со слоем. Направление распространения излучения s в точке слоя C(r) = C(x, y, z), z [0, H], описываем сферическими координатами: азимутом [0, 2] и зенитным углом [0, ]. Азимут отсчитывается от оси x к оси y в направлении часовой стрелки, если смотреть в положительном направлении оси z, т. е. сверху вниз. Азимут внешнего параллельного потока на верхнюю границу z = 0 считаем равным нулю, т. е. ось x и направление внешнего потока s0 лежат в одной плоскости. Зенитный угол отсчитывается от положительного направления оси z. Направления s+ = {+, } +, s = {, }. Орты сферической системы координат:

e = { sin, cos, 0}, e = { sin cos, sin sin, cos }, r e = { cos cos, cos sin, sin }.

Состояние поляризации луча s описываем вектором Стокса (s) (r, s) (верхний индекс (s) далее опускаем), определенным в системе координат с базисными векторами e1 = e, e2 = e и s, а вектор (r, s ) – в – базисе e1, e2, s. Вектор e1 лежит в меридиональной плоскости единичной сферы, описанной в точке слоя C(r) (меридиональная плоскость содержит ось z и направление распространения луча s;

см. рис. 1.9). Вектор e перпендикулярен этой плоскости. Для определения матрицы рассеяния (см.

рис. 1.6) используем систему координат, связанную с плоскостью рассеяния, – – 86 Глава 1. Одномерные плоские задачи это две тройки векторов: e1, e2, s – для рассеянного света и 1, 2, s – для – – ee падающего.

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

Параметры Стокса падающего излучения, описанные в Этап 1.

базисе e1, e2, s, переводятся в параметры Стокса падающего излучения, описанного в базисе 1, e2, s, связанном с плоскостью рассеяния, путем e поворота на угол базиса e1, e2, s к базису 1, e2, s вокруг вектора s до e 1 :

совмещения e1 с e = L( ).

Переход от параметров Стокса падающего излучения, Этап 2.

описанных в базисе e1, 2, s, к параметрам Стокса рассеянного излучения, e заданным в базисе 1, 2, s, производится с помощью матрицы рассеяния :

ee =.

Этап 3. Обратный переход от параметров Стокса рассеянного излучения, заданных в базисе 1, 2, s, связанном с плоскостью рассеяния, к параметрам ee, заданным в базисе e1, e2, s, получается путем поворота первого базиса вокруг оси s на угол до совмещения вектора 1 с вектором e1 :

e = L().

Полученная матрица преобразования (1.158) P (z, s, s ) = P (z,,,, ) = L()(s )L( ) называется фазовой матрицей. Конкретный вид матрицы поворота L() и значения угла поворота зависят от взаимного расположения векторов s и s.

Рассмотрим случай, для которого 0 (рис. 1.9, a). Случаю 0 отвечает рис. 1.9, b. Плоскость рассеяния образует с меридиональными плоскостями и углы и – это углы между – векторами e1 и 1, e1 и e1. В данном случае, 2. Положительный e поворот на угол осуществляется против часовой стрелки, если смотреть против направления распространения луча (против вектора s).

Для того чтобы описать процесс рассеяния из направления s в направление s, необходимо повернуть вектор 1 на угол или на угол. Эта e двойственность связана с тем, что одновременный поворот базисных векторов e1 и 2 на угол вокруг вектора s не меняет значения параметров Стокса.

e Далее с помощью матрицы рассеяния осуществляется процесс рассеяния, причем рассеянный луч описывается в базисе 1, 2, s. Затем производится ee поворот на угол или на угол, т. е. переход к базису e1, e2, s. Таким 1.2. Поляризационные задачи z b e e e e D e2 A e s s e s y C(r) e x a z b e e e e A e D e s s s e e y C(r) x b Схемы определения матрицы поворота системы координат при Рис. 1.9.

0 (a) и 0 (b) 88 Глава 1. Одномерные плоские задачи B = = e 1 e e 2 1 1 e D A e e e e Рис. 1.10. Схема цепочки поворотов базисных векторов e1 = e1 = 1 = e e образом, в данном частном случае процесс рассеяния может быть описан с использованием фазовой матрицы P любым из трех способов:

P (,,, ) = L() (s ) L( ), P (,,, ) = L( + ) (s ) L( ), P (,,, ) = L() (s ) L( ).

Первый случай соответствует цепочке поворотов e1 = 1 = e1 = e e (рис. 1.10), второй – 1 = 1 = e1 = e1, третий – e1 = 1 = 1 = e1.

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

Фазовая матрица удовлетворяет следующим основным соотношениям симметрии:

P (,, ) = D1 P (,, ) D2, P (,, ) = P (,, ), P (,, ) = D2 P (,, ) D2, D1 = diag(1, 1, 1, 1), D2 = diag(1, 1, 1, 1), которые оказываются полезными при переходе к другим системам координат, а также для разработки более экономичных вычислительных алгоритмов.

Например, умножение справа и слева фазовой матрицы на D2 соответствует следующим заменам:

1) изменение направления, в котором отсчитывается азимут;

2) смена знака направляющих косинусов (замена системы координат);

3) использование параметров поляризации {I, Q, U, V } вместо {I, Q, U, V }.

1.2. Поляризационные задачи Знаки в матрицах поворота L зависят от конкретного расположения направлений падающего s и рассеянного s лучей. Всего возможны четыре варианта взаимного расположения s и s. Согласно этим вариантам знак поворота g зависит от разности азимутальных углов следующим образом:

180 180 180 1 g +1 g + Значение g можно определить, например, так: g = sign[ sin( )].

Перепишем формулу (1.158) в общем виде:

P (,,, ) = L(g) (s ) L(g ).

Согласно формулам сферической тригонометрии из решения сферического треугольника ABD (см. рис. 1.9) получаем cos s = cos cos + sin sin cos( ), cos + cos cos s + cos = =, 1 2 1 sin sin s cos + cos cos s + cos = =, sin sin s 1 2 1 cos 2 = cos2 1, 1 cos2, sin 2g = 2g cos cos 2 = cos2 1, 1 cos2.

sin 2g = 2g cos Существует несколько предельных случаев взаимного расположения s и s, когда значения cos и cos нельзя получить по этим формулам. Предельные случаи представлены в таблице cos cos 0 0 0 0 1 + 0 1 0 + = 1 ±[ ( )] cos( ) 0 ±( ) cos( ) =0 1 1 ±[ ( )] cos( ) 0 = 1 cos( ) =0 1.2.3. Векторная задача для плоского слоя с однородной отражающей границей. В предположении стационарного состояния среды и постоянства внешнего потока поле квазимонохроматического поляризованного излучения 90 Глава 1. Одномерные плоские задачи полностью описывается четырехкомпонентным вектором t (r, s), компонен тами которого являются параметры Стокса, имеющие размерность интен сивности излучения. Если среда макроскопически оптически изотропна и плоскостратифицирована, то полный вектор Стокса t находится как решение векторной краевой задачи теории переноса поляризованного излучения.

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

Основными методами решения задач с поляризацией являются: метод Монте-Карло, матричный метод (метод удвоения и метод сложения слоев), метод сферических гармоник;

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

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

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

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


Остановимся подробнее на решении векторной задачи Dz = S + F, = 0 (s), = H (s), (1.159) 0 H где вектор-функция источника S s (z) P (z,,, ) (z,, ) d d. (1.160) Ядром интегрального оператора (1.160) является фазовая матрица P.

Граничные условия записываются в общем виде:

= F0 ( 0 ) ( 0 ) + f0, 1.2. Поляризационные задачи где F0 – внешний мононаправленный поток излучения на верхней границе – слоя z = 0 в направлении 0 и f0 – диффузный источник;

– = RH + fH, H где оператор RH описывает отражение излучения от подстилающей поверх ности, fH – вектор-источник диффузного излучения на дне.

– При ламбертовой поверхности 2 q RH = R = 0 (H,, ) d tH, d 0 где q – альбедо поверхности;

0 – матрица изотропного рассеяния.

– – При квазизеркальном отражении RH (H,, ) = AH () (H,, ).

Условие AH () = AH = const соответствует зеркальному отражению. Могут использоваться законы отражения, описываемые произвольными матрицами отражения.

В силу линейности краевой задачи (1.159) ее решение можно представить в виде суперпозиции трех вектор-функций:

= 0 + 1 +, 0 01 + 02, где 01 отвечает нерассеянному излучению от мононаправленного источника, 02 – излучению от диффузных источников, 1 – однократно рассеянному, – – а – многократно рассеянному излучению. Эти вектор-функции находятся – из следующих краевых задач:

Dz 01 = 0, 01 = F0 (s s0 ), 01 = 0;

(1.161) 0 H Dz 02 = F, 02 02 (1.162) = f0 = fH ;

0 H Dz 1 = S01, 1 1 = RH 01 ;

= 0, (1.163) 0 H Dz = S + S1 + S02, (1.164) = RH + RH 1 + RH 02.

= 0, 0 H Задача (1.161) решается аналитически:

(z) s +, F0 ( 0 ) exp, 01 = s.

0, 92 Глава 1. Одномерные плоские задачи Аналитические выражения получаются для функции источника за дачи (1.163) (z) S01 = s (z) exp P (z,, 0, 0 ) F0 ;

(H) RH 01 = q0 exp 0 F при ламбертовой поверхности или 1 (H) RH 01 = AH () (0 + ) ( 0 ) exp F0, если присутствует квазизеркальное отражение.

Решения задач Коши (1.162) и (1.163) находятся интегрированием по харак теристике дифференциального оператора Dz, определяемой соотношениями:

= const, = const. В плоской задаче такие характеристики совпадают с направлениями распространения излучения s. Так что имеем:

для s + z (z) (t) 1 (z) + F(t, +, ) exp (z,, ) = + dt + f0 exp ;

+ + (1.165) для s H (t) (z) F(t,, ) exp (z,, ) = dt + | | | | z (H) (z) +fH exp (1.166) ;

| | для s + z (t) (z) (t) + P (t,+,0,0 )F0 dt;

(z,,) = + s (t) exp + (1.167) для s H (t) (t) (z) P (t,,0,0 )F0 dt+ (z,,)= s(t)exp | | | | z (H) (z) +RH 01 exp (1.168).

| | Задача (1.164) решается путем введения естественного итерационного процесса последовательных приближений по кратности рассеяния:

Dz N = SN 1 + S1 + S02, (1.169) N = RH N 1 + RH 1 + RH 02, N = 0, 0 H 1.2. Поляризационные задачи где N = 2, 3,... – номер итерации. Условно принимаем N = 0 для за – дачи (1.161), N = 1 для задач (1.162) и (1.163). В качестве критерия сходимости итерационного процесса выбираем сходимость в метрике пространства C (4) с заданной относительной точностью итераций :

N N (1.170) max.

N {z,, } Менее жестким условием сходимости итераций является условие для более гладкой функции степени поляризации [(QN )2 + (U N )2 + (V N )2 ]1/ PN = (1.171).

IN В краевой задаче (1.159) можно сделать переход от пространственной геометрической переменной z к оптической толщине (z). Этот случай является частным в нашей методике.

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

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

1.2.4. О симметрии решения плоской векторной задачи. Векторная кра евая задача для уравнения переноса поляризованного излучения чрезвычайно трудоемкая в вычислительном аспекте задача. Эффективная реализация на ЭВМ вычислительных алгоритмов потребовала детальной проработки и анализа свойств этой задачи. Особого рассмотрения заслуживает вопрос о свойствах симметрии по азимуту решения векторной краевой задачи с внешним мононаправленным потоком излучения. В скалярной задаче при условии, что поверхностные и объемные источники излучения являются четными функциями азимута, решение оказывается четной функцией. В 94 Глава 1. Одномерные плоские задачи задачах переноса поляризованного излучения даже при блочно-диагональных матрицах рассеяния такой симметрии нет.

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

Действительно, элементы фазовой матрицы состоят из некоторой ком бинации синусов и косинусов углов 2 и 2 и элементов ij матрицы рассеяния (cos s ). Элементы матрицы рассеяния как функции косинуса угла рассеяния s суть четные по азимуту функции. Комбинации синусов и косинусов, входящие в матрицы поворота, могут быть как четными, так и нечетными функциями азимутов.

Решение на n-й итерации n представим как сумму n = + n + n четной + n и нечетной n компонент. Удается получить следующее правило знаков для каждой компоненты i, i = 1 4, при переходе от (n 1)-й к n-й итерации:

+ + + + n1 = + n, n1 = n.

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

Прослеживая по итерациям изменение четности параметров Стокса, можно сделать следующие заключения о решении.

Случай 1. Источник излучения неполяризован: F = {I, 0, 0, 0}, матрица рассеяния рэлеевская:

1 = + 1, 1 = + 1, 1 = 1, 1 0, 1 1 2 2 3 n = + n, n = + n, n = n, n 0.

1 1 2 2 3 Неполяризованный источник: F0 = {I, 0, 0, 0}, матрица Случай 2.

рассеяния блочно-диагональная (n 2):

1 = + 1, 1 = + 1, 1 = 1, 1 0, 1 1 2 2 3 n = + n, n = + n, n = n, n = n.

1 1 4 2 2 3 Случай 3. Нормальное падение внешнего потока неполяризованного излу чения, слой рэлеевский. Компоненты 3 и 4 равны нулю, 1 и 2 полностью 1.2. Поляризационные задачи симметричны по азимуту. Этот случай рассмотрен в работе М. Г. Кузьминой, в которой для исследования свойств симметрии использовалось разложение в ряд по обобщенным сферическим функциям.

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

Из анализа азимутальной симметрии вектор-функции многократного рас сеяния вытекает, что в ряде случаев (среди них оказываются наиболее распространенные задачи оптики атмосферы планет) достаточно находить отдельно четные и нечетные по азимуту составляющие в интервале [0, ], а полное решение получать так:

= + + = + для [0, ], для [, 2].

1.2.5. Азимутально симметричная задача. Решение задачи (1.159) при однородной ламбертовой границе с помощью формулы учета вклада подсветки можно осуществлять в два этапа. Сначала решить задачу (1.159) при абсо лютно черном дне, а затем решить задачу для векторной функции пропускания W0 (z, s) (см. гл. 3). При большом разнообразии альбедо естественного ландшафта Земли этот путь существенно сокращает объем вычислительных работ, поскольку векторная задача для W0 обладает свойствами симметрии по азимуту. Расширим постановку этой задачи, включив в нее диффузные объемные и поверхностные источники излучения. Такие источники возникают, например, в задачах распространения длинноволнового и СВЧ-излучения.

Итак, рассмотрим следующую задачу:

(1.172) Dz W = SW + FW (z, ), W = f0 (), W = fH ().

0 H Наиболее распространенной является задача с неполяризованными источни ками излучения:

tH = {1, 0, 0, 0}.

F(z) = F (z) tH, f0 = f0 () tH, fH = fH () tH, Другой важный частный случай – задача для W0 (z, s). При любых – матрицах рассеяния вектор W(z, s) является симметричным по азимуту, а следовательно, задачу (1.172) достаточно решать для усредненной по азимуту вектор-функции W(z, ) = W(z,, ) d.

Из этого свойства вытекает вывод о том, что вклад подсветки в задаче с учетом поляризации, как и в скалярной задаче, следует учитывать только при 96 Глава 1. Одномерные плоские задачи расчете нулевой гармоники по азимуту, если использовать как первый этап предварительное разложение решения в ряд по азимутальным гармоникам.

В случае матриц рассеяния типа Ми устанавливается также, что векторная задача (1.172) для W(z, ) расщепляется на две системы уравнений или усеченное векторное уравнение переноса только для компонент W 1 и W (W 3 = W 4 = 0), если обозначить W = {W m }, m = 1 4:

Dz W 1 = s (z) 0 W 1 d + 0 W 2 d, 11 (1.173) W = 0, = 0;

W 0 H Dz W 2 = s (z) 0 W 1 d + 0 W 2 d, 21 (1.174) W W = 0, = 0;

0 H 12 = 12 cos 2, 21 = 12 cos 2, 11 = 11, 22 = 22 cos 2 cos 2 + 33 sin 2 sin 2;

ij = ij d – усредненные по азимуту элементы фазовой матрицы. Система уравне – ний (1.173)–(1.174) описывает пропускание параллельной и перпендикулярной составляющих интенсивности поляризованного излучения:

I = I Q, I = I + Q.

Нулевые азимутальные гармоники элементов фазовой матрицы определяются через интегралы ij (z,,, ) y cos( ), ij (z,, )= dy, (1 y 2 )1/ и вычисляются с помощью квадратуры Чебышева–Эрмита.

Для решения усредненных по азимуту уравнений также используется ме тод последовательных приближений по кратности рассеяния и интегрирование по характеристике = const с гауссовыми квадратурами на [1, 1].

1.2.6. Дискретный аналог векторного уравнения переноса поляризо ванного излучения. Описанное выше разбиение исходной векторной за дачи (1.159) на четыре краевые задачи (1.161)–(1.164) предполагает использо 1.2. Поляризационные задачи вание для решения уравнения переноса метода характеристик с итерациями по кратности рассеяния и квадратурами на единичной сфере. Выделение задач (1.161)–(1.164) дает возможность проводить итерации единообразно, независимо от типа и формы задания источников излучения. При этом функции источников в задаче (1.164) для вклада кратного рассеяния уже не содержат особенностей, которыми обладают некоторые исходные источники излучения. Принятые алгоритмические решения позволили в несколько раз сократить время расчетов и автоматизировать управление расчетом в зависимости от начальных данных конкретной задачи с произвольными матрицами рассеяния.

Краевые задачи (1.161)–(1.163) могут быть решены аналитически путем интегрирования по характеристике. Численное решение задачи (1.164) выпол нимо только для сеточных функций, вводимых на разностной сети задачи.

Численное решение задачи (1.164) осуществляется посредством интегриро вания по характеристикам дифференциального оператора Dz с переходом к конечно-разностному аналогу и вычислением функции источника (1.160) и операторов отражения с помощью квадратур на единичной сфере на дискретном подмножестве (разностная сеть) множества, на котором определена вектор-функция:

= [0, H] [1, 1] [0, 2], = [0, H] [1, 1] [0, 2], где [0, H] = {zk }, k = 1 K, – пространственная сеть по z;

[1, 1] = {j }, – = + [1, 0), j = 1 J +, = [1, 0), j = 1 J, – разностная сеть по – j j = { }, i = 1 I, – разностная зенитному углу [0, ], = cos ;

[0, 2] – i сеть по азимуту.

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

Рассмотрим общий случай, когда расчет проводится для [0, 2]. На множестве определим сеточные вектор-функции по образцу:

+, s +, kji (zk, ±, i ) kji (1.175) = j, s.

kji 98 Глава 1. Одномерные плоские задачи Каждая компонента сеточной вектор-функции kji имеет размерность N = K (J + + J ) I. При этом функции +, отвечающие направлениям kji нисходящего излучения с = +, имеют размерность компонент N = + j K J + I, а направлениям восходящего излучения с = отвечают j функции с размерностью компонент N = K J I, т. е. общая размерность + kji равна N = N + N. Каждый вектор типа (z,, ) будет иметь дискретный аналог размерностью N N, где N – число ненулевых компонент – вектора Стокса. К таким вектор-функциям относятся 0, 1, B1, B2, B,, для которых аналогично (1.175) вводятся сеточные вектор-функции с сохранением метки функции.

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

По существу в вычислительных алгоритмах выделяются следующие основные модули:

1) расчет функции 1, равной 1 или 02, или 1 +02, интегрированием по характеристике;

2) расчет интенсивности N it в приближении N it-й кратности рассеяния на каждой итерации интегрированием по характеристике для сеточных функций;

3) расчет интегралов на единичной сфере (функции источника B) по квадратурным формулам;

4) расчет граничных условий.

Для неазимутальной задачи добавляется операция усреднения по азимуту элементов матриц рассеяния и матриц поворота. На ЭВМ с большой оперативной памятью предварительно рассчитываются сеточные значения элементов матрицы рассеяния и элементов матриц поворота.

1.2. Поляризационные задачи Наибольшие затраты приходятся на вычисление интегралов столкновений.

Действительно, полная фазовая матрица с 16-элементной матрицей рассеяния P (,, ) = L() (cos ) L( ) = (1.176) 12 2+13 11 12 1+13 2 21 1+31 2 22 1 1+23 2 1 + 22 2 1+23 1 1 24 1+34 +32 1 2 +332 2 32 2 2+33 1 = 21 2+31 1 22 223 2 + 22 2 23 2 24 2+34 1, 1 2 +32 1 1+33 2 1 32 2 1+33 1 42 2+43 41 42 1+43 2 где = cos 2, 1 = cos 2, 2 = sin 2, 2 = sin 2, порождает 36 скалярных интегралов. Рэлеевская матрица порождает 14, а блочно-диагональная матрица Ми – 18 скалярных интегралов. А это означает, что решение одной векторной – задачи на каждой итерации эквивалентно решению 36, 14 или 18 скалярных задач соответственно.

В случае сильной анизотропии рассеяния предлагается векторный аналог выделения -анизотропии. Представим матрицу рассеяния в виде () = A (1 ) (1) + 1 (), cos s, с коэффициентом A [0, 1], подлежащим определению. С физической точки зрения это представление отвечает тому, что для закона рассеяния в акте принимается следующее приближение: доля A рассеянного излучения рассеивается точно вперед в соответствии с законом рассеяния (1), а остальное излучение рассеивается в соответствии с законом рассеяния 1 (), который отличается от () лишь в достаточно узкой области малых углов рассеяния s. Векторное уравнение переноса при этом несколько изменится:

d + t (z) G = s (z) P1 ds + F, dz G E (1), P1 = L() 1 L( ), s A s s, 0 1, a1 () d.

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

100 Глава 1. Одномерные плоские задачи 1) генерация таблиц для последовательности расчета линейных операций типа A = P при вычислении интегралов столкновений B = S, – вектор – решения в разных приближениях;

2) генерация таблиц для последовательности линейных операций типа An = (P F)n, n = 1 4, при вычислении первого рассеяния, F – вектор – Стокса источников излучения;

3) формирование массива признаков для определения нулевых и ненулевых компонент решения в разных приближениях;

4) генерация таблиц для последовательности расчета линейных операций типа 1) или 2) при вычислении граничных условий;

5) формирование массива признаков для определения нулевых и ненулевых компонент при логической разметке магнитной ленты или магнитного диска;

6) моделирование работы разных этапов задачи без выполнения расчетов (например, при продолжении счета после прерывания).

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

Для кодировки каждого слагаемого используются следующие шесть величин:

a1 = sign(ik lj ) = {+1 или 1} – знак, с которым входит слагаемое, – содержащее произведение элементов матриц поворота;

a2 – номер комбинации N S(ik lj ), принимающий значения от 1 до 9 в – соответствии с таблицей, a3 = r = (k1) 4 +l – номер элемента матрицы рассеяния kl, k, l = 14, – r = 1 16;

a4 = i – номер компоненты вектора A;

– a5 = j – номер компоненты вектора ;

– a6 – специальный признак.



Pages:     | 1 || 3 | 4 |   ...   | 15 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.