CC BY
67
2006 ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА. Сер. I. Вып. 3
ФИЗИКА
УДК 551.510:551.521 А. В. Васильев
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ИНТЕНСИВНОСТИ МНОГОКРАТНО РАССЕЯННОГО СОЛНЕЧНОГО ИЗЛУЧЕНИЯ И ПРОИЗВОДНЫХ ОТ НЕЕ С УЧЕТОМ СФЕРИЧЕСКОЙ ГЕОМЕТРИИ АТМОСФЕРЫ (КОМПЬЮТЕРНЫЙ КОД 8САТ1Ш)
Введение. Одной из актуальных задач при реализации различных дистанционных методов измерений параметров атмосферы и поверхности Земли является численное (компьютерное) моделирование измеряемых величин. Наличие подобных моделей и соответствующих им компьютерных кодов позволяет исследовать зависимости этих величин от параметров атмосферы и поверхности, знание которых в дальнейшем используется в различных методах интерпретации результатов измерений. Вычисление интенсивности излучения для заданных моделей атмосферы и поверхности, а также расчет линеаризованного оператора прямой задачи относятся к основным элементам в большинстве методик обращения радиационных измерений. При этом, вследствие все возрастающей точности современной аппаратуры, важным требованием к численным моделям становится достаточно адекватный учет разных процессов, формирующих измеряемую величину, в частности - многократного рассеяния излучения и сферичности атмосферы.
В данной работе рассматривается моделирование измерений интенсивности рассеянного солнечного излучения в атмосфере Земли, которые применяются в различных современных методах изучения газового и аэрозольного составов атмосферы и характеристик подстилающей поверхности [1]. В настоящее время для решения указанной задачи создало немало компьютерных кодов, например [2-5], однако по ряду причин, подробно обсуждающихся в [6], они не обладают некоторыми необходимыми для научно-исследовательских задач свойствами, в частности возможностью гибкого использования разных моделей атмосферы, поверхности и измерительного прибора, а также расчета производных от моделируемой величины по любым параметрам атмосферы и поверхности, что особенно важно для задач интерпретации измерений [7, 8].
Модель переноса излучения. Моделируемая величина является сверткой монохроматической интенсивности и аппаратной функции прибора [1, 8, 9]. Поэтому в конечном счете задача моделирования измерений сводится к расчету именно монохроматической интенсивности, что и будет рассматриваться ниже. Обзор алгоритмов учета аппаратной функции прибора приведен в [9]. Кроме того, как известно, монохроматическое приближение достаточно точное вне селективных полос молекулярного поглощения, в частности в видимой области спектра, где находится максимум интенсивности рассеянного солнечного излучения.
© А. В. Васильев, 2006
При изложении модели переноса излучения удобно приводить основные соотношения для случая плоской атмосферы, а уже после - их обобщения для сферической геометрии.
Поле рассеянного солнечного излучения в плоской атмосфере при пренебрежении рефракцией и поляризацией описывается (см., например, [1]) уравнением
2тг 1
dI{T,ri0,Ti,tp0,<p) Uo(г) f,,f,l( \ т( / .
V-^-=-I(T>Vo>V>tPo,<P) +dtp dr¡ х(т,х)1{т,щ,г) ) +
О -1
+Fo^-x(TyXo)P(0>r,rlo). (1)
Здесь x = m' + \/(l ~ ??2)(1 - (r,í')2) eos(tp - <p'); xo = ЩП + y/(l - 7?g)(l - r¡2) cos(<p -V?o), I(0,r)o,T),(po,(p) = 0 при r¡ > 0; J(t0, 770, 77, (p0, y) = 0 при 77 < 0; / - интенсивность излучения, определяемая в точке с координатами (z, при координатах солнца ($о,<ро), z - высота, р и _ зенитные углы излучения и солнца соответственно, при
Zoo
этом r¡ = cosí?, r¡o = cosido; v3 и Vo _ азимуты излучения и солнца; г = f a(z')dz' - оп-
z
тическая глубина, соответствующая высоте z, в выражении для которой а - объемный коэффициент ослабления атмосферы, гж - высота ее верхней границы. Оптическая глубина однозначно связана с высотой подстановкой dr = —a(z)dz, поэтому в дальнейшем из соображений удобства будем записывать различные выражения в координатах как оптической глубины г, так и высоты z. Величина cjq = а/а - альбедо однократного рассеяния, или вероятность выживания фотона (а - объемный коэффициент рассеяния, а — а + к, где к - объемный коэффициент поглощения); - индикатриса
i
рассеяния, нормируемая по условию | J x(x')dx' = 1; Ро - поток солнечного излуче-
-1
ния на верхней границе атмосферы; tq - оптическая глубина всей атмосферы, соответствующая высоте ее нижней границы zq. Наконец, P(ri, 72,77) - функция пропускания плоской атмосферы между оптическими глубинами ту и т2, определяемая как P(Ti,T2,v) — ехР(— lr2 — Til/7?)- Задачей моделирования является нахождение интенсивности при заданных параметрах визирования (z, $0, <Р, <Ро) и атмосферы a(z), k(z), x(x,z).
Сумма двух последних слагаемых в (1), по определению, есть функция источников B(T,r¡o,r¡,ípo,if)- Следствием (1) является выражение интенсивности через функцию источников, которое удобно записать не в классическом виде двух формул [1, 8], а следующим образом:
т
I(T>Vo,V,Vo,<p) = ~ / B(T',rio,ri,<po,<p)P(T',T,ri)dT', (2)
V J
1-0(1-6(7)))
где 0(77) = 0 при 77 < 0 и ©(77) = 1 при 77 > 0 (функция Хэвисайда). Оставляя в функции источников только слагаемое без интеграла, получаем из (2) решение задачи в приближении однократного рассеяния [1, 8]. Для сферической геометрии (см. [6]) оно приобретает вид 4
h
i(z,vo,v,<Po,v) = 2; Ja(z(l »x(z(l )>Xo)Poo(l )P(l )dl, (3)
/1
в котором I - координата трассы распространения излучения с началом отсчета в точке, ближайшей к центру Земли, Роо{1 ) и Р{1 ) - функции пропускания атмосферы на трассах «точка рассеяния - солнце» и «точка рассеяния - точка визирования» соответственно (см. подробнее [6], где также описано нахождение пределов интегрирования /1, 12). В сферическом случае необходимо учитывать зависимость всех координат от точки трассы визирования /, что приводит к формулам
z(l)=r(l)-R, v(D = -l/r(l), {l)=(K + *b)Vo,b-((R + zb)Vb + l)Xo^
r{l )
r(i) = ^i2 + (R + zby(i-vt), cosЫП-МП) =
((R + zb)( 1 - rj%) - lm)Jl - 77^ cos(<pb - <p0,b) + ((R + zb)T)b + I )i?o,bv/l ~ *1ь
=----:-,- (4)
r(i Wi-nlU)
В (4) R - радиус Земли, а все величины с индексом b относятся к точке визирования. Для общего случая многократного рассеяния формальное выражение для решения уравнения (1) может быть записано в операторной форме [1, 8]. Интегральным оператором [10], порождаемым некоторой функцией а(х,х') и ставящим в соответствие
ь
функции f(x) функцию д(х), называется соотношение д(х) = f а(х, x')f(x')dx'. Сим-
а
волически оно записывается g = Af, функция а(х, х') называется ядром интегрального оператора А. В рассматриваемом случае переменных не одна, а три - (г, 77,9?), но в принципе это ничего не меняет - однократный интеграл заменяется тройным, а ядро является функцией шести переменных.
Решение уравнения (1) записывается как ряд Неймана
I = To(bi + Tibi + T?b! + T?bi + ...), Ti = T0Si, (5)
где I - искомая интенсивность излучения I(z,r]o,ri,(po,<p)-, To - оператор переноса прямого излучения с ядром ¿0(2,77, р, z',ri',<p') = (6(77)0(2 — z') + (1 — 0(г/))(1 — в(z - z')))a(z')P(z', z, 77)6(7?' - Ti)S(<p' - <p)l\n\\ <5()- дельта-функция [10]; Si - оператор однократного рассеяния с ядром si(z,r},<p, z',r]',ip') = u)0(z)x(z,x)S(z — z')/47r;bi-функция источников однократно рассеянного солнечного излучения ¿>1 (т, 770,77, > V7) = f0u)0(t)x(t, хо)р(0, т, 77о)/4тг.
Учет отражающей поверхности. Уравнение переноса (1) записывается для чисто рассеянного излучения, без отражения от поверхности Земли, учет которого необходимо добавить в модель переноса [1, 8].
В качестве количественной характеристики отражения от поверхности в изучаемом случае удобно пользоваться функцией отражения ¥>2>$ь¥>i)> определяемой (см. [1]) как
.¿7г
где (19ь</?1) - направление падающего излучения; ($2, </>2) ~ направление отраженного излучения; I - как и выше, интенсивность излучения (но в данном случае - для отражения при отсутствии атмосферы). Важной энергетической характеристикой поверхности является доля отражаемой ею энергии падающего излучения, характеризуемая альбедо поверхности а. Из его определения (см. подробнее [1]) следует связь
2я 7Г
= ! &Р2 J у( , VI, 1?2, Уг) вт ■ (7)
о тг/2
Для описания модели отражения достаточно указать конкретный вид функции
2,<¿>2,</>1)- Отметим две стандартные модели. Для идеальной зеркальной поверхности отраженное излучение присутствует только в одном направлении (угол отражения равен углу падения), следовательно,
у (1?2 , <р2, 1?Ь ) = 27ГГ(1?! )6(#2 ~ Ж¥>2 ~ <¿>1 )•
Здесь г(я91) - коэффициент зеркального отражения, вычисляемый (см., например, [1]) как следствие из известных формул Френеля
т2 сое г?1 — у/т2 - 1 + сой2 #1 + сое — у/т2 - 1 + соб2 1?!
т2 сов + л/т2 — 1 + сов2 1?! соб^! + у/т2 — 1 + сой2 1?!
(8)
в котором т - комплексный показатель преломления поверхности, полностью описывающий рассматриваемую модель. Для изотропной отражающей поверхности направление отраженного излучения не зависит от направления падающего, и альбедо а является константой, также полностью описывающей отражение. В этом случае
1{02,<Р2) = -АР\
(9)
поток излучения, падающий на поверхность.
Рассмотрим сначала учет отражения от поверхности в приближении однократного рассеяния.
Для идеальной зеркальной поверхности, для краткости опуская в формулах координаты (г?о,^о), будем иметь
1(г, ф) = (г, дм + /(3) ч>) + ^ ч>) + 1{5) ч>), где ,$,</?) - интенсивность однократно рассеянного в атмосфере излучения, расчет которой обсуждался выше; ,!?,</?) - интенсивность излучения, отраженного от поверхности после однократного рассеяния; /(4)(.г,19,<£>) - интенсивность излучения, однократно рассеянного после отражения от поверхности; д, (р) - интенсивность излучения, отраженного от поверхности, до этого однократно рассеянного после отражения от поверхности (два отражения).
Для плоской атмосферы (2,$, <р) = г(тт—'д)1^(го,тт—'в,<р)Р(го,г,'в). Для нахождения 19,</з) достаточно рассчитать однократно рассеянное излучение от источника, освещающего атмосферу со стороны нижней границы с интенсивностью зеркально 6
отраженного солнечного излучения, что приводит просто к замене внеатмосферного потока .Ро на ГоР(г00, ^оЖ^о), а также щ на — щ. Для последней компоненты можно записать = г(тг - тг -
В сферическом случае вся сложность заключается в вычислении компоненты д, </з), для которой нельзя заранее приписать координаты источнику при «перенесении» его на поверхность из-за их изменения на трассе распространения излучения. Выход - в модификации соотношения (3) для одновременного вычисления рассеянного и отраженного излучений. Запишем сразу итоговую формулу, а потом объясним ее. Итак, вместо (3) вычисляем
= ^ I а(г(1))Р(1)(х(г(1)уХо)Роо(П + к
+ &Ы1) - г,"(1 ))х(г(1), Хо(1 )Ш(1 )Р'(1)) Л,
„,п _ У(211 + г(1) + г0)(г(1)-г0) Г1{)~ Я + г(1)
Хо(0 = ) + у/(1-Ш1))2)(1-Ч2(1))со8(<р - М1)),
/ а + (Я + *о)2(1 - КС ))2) " я) Л'
^ (я+гоНС) )
Р^(1)=ех р
Р'(/) = ехр
I а (у(/')2 + (Д + г(0)2(1-К(0)2) - я) я'
■ (Ю)
В (10) вторая и четвертая компоненты интенсивности в сферической атмосфере рассчитываются по одной трассе визирования. Отраженный от поверхности с направлением —77о(0 луч солнца должен прийти в точку рассеяния г{1) под направлением —г)'0(1), которое для обратной трассы из г{1) с направлением (г]'0(1 ),<р'0(I )) в точке пересечения с поверхностью будет иметь зенитный угол и азимут, равный солнечному (условия зеркального отражения). Эти условия дают 1р'0(1) = 7г + </Зо(0> и все сводится к нахождению г)'0(1), для которого из формул (4) следует уравнение
(Д + г{1))щ{1) - ((Д + г(1М(1) + 1'Ыо(1)))ХоШ)) + Пч 1,(1)) = 0,
ГШ)) = -у/(л + ^о)2-(Д + ^(0)2(1-(^(0)2),
хоШО) = т(1Н(1) - у/(1-Ыт(1-Ш))2),
легко решаемое численно по итерационной формуле Ньютона. Компонента /(3+5) (z, ф) вычисляется аналогично, как и в плоском случае, при переходе к источнику на нижней границе атмосферы.
Для изотропного отражения в рамках приближения однократного рассеяния в суммарной интенсивности присутствуют только две компоненты: рассеянное атмосферой и отраженное от поверхности (без рассеяния). Вторая компонента в плоском случае, как следует из (9), равна ^AF0cos'doP(z00,zo,'do)P(zo,z,'d). Для сферического случая структура этого выражения такая же, только следует учесть, согласно (4), изменения всех углов на трассах от верхней до нижней границы атмосферы и от нижней границы до точки визирования.
Рассмотрим теперь общую схему учета отражения в случае многократного рассеяния [8, 11]. Вернемся к ряду Неймана (5) и заметим, что в нем bi = Tibo, где формально bo(z,r)o,r),(po,(p) = F05(z — Zoo)5(r) - щ)8{ч> — </?o)M/a(zoo), а это позволяет записать (5) в виде I = ТЬо, в котором Т = Т0(1 + Ti + Т2 + Tf + ...) - оператор переноса излучения. Процесс отражения эквивалентен появлению на нижней границе атмосферы дополнительных источников излучения [1, 8, 11], которые охарактеризуем функцией ВГ)о- Для нее (см. подробнее [11]), согласно (б), имеем ВГ)о = Rii, где i - интенсивность излучения до отражения, Ri - оператор однократного отражения с ядром Г1(т,77,<Лт',7/,¥>') = t1 ~ ®(r)))®(r)')r)'y(v,v,'n'уР'Жт' ~ т0)6(т - т0)/2тг. Дополнительные источники излучения на поверхности вновь изменяют падающую интенсивность i, что приводит к ряду
I = (Т + To(RiT) + To(RiT)2 + To(RiT)3 + ...)b0, (11)
который теперь вместо (5) является формальной записью решения задачи при учете отражения.
Решение (11) получено для плоской атмосферы. Однако достаточно очевидно, что его формальный вид не изменится и при учете сферичности. В самом общем случае (включая неоднородную по горизонтали атмосферу) необходимо принять во внимание зависимость всех функций и ядер от географических координат (широты и долготы), что поднимет кратность интегрирования до пяти, но явный вид ядер выражения (11) тогда можно и не выписывать.
Расчет интенсивности многократно рассеянного излучения методом Монте-Карло. Достоинствами метода Монте-Карло (ММК) [12] по сравнению с другими методами расчета поля рассеянного излучения являются [8, 13]: возможность непосредственного использования индикатрис рассеяния (без каких-либо их преобразований типа разложений в ряды [1, 8, 13], неизбежно приводящих к искажениям); простота метода и его физическая наглядность; возможность без проблем учитывать поляризацию, анизотропию отражения, горизонтальную неоднородность и сферичность атмосферы (последнее и предопределило использование метода в данной работе). Главный недостаток ММК - наличие у получаемых им результатов случайной погрешности (в данном смысле это действительно полный аналог измерений). Некоторые приемы ее уменьшения рассматриваются ниже. Наглядность ММК позволяет излагать его алгоритмы, отталкиваясь от физического смысла всех его операций [1, 8]. Однако такое изложение весьма объемно. Потому здесь вынуждены привести лишь основные используемые в ММК соотношения, в большинстве случаев без выводов и обоснований, отсылая читателей за ними к работам [1, 8].
Суть ММК составляет компьютерное моделирование некоторого случайного процесса и получение искомых результатов как статистических его характеристик. Исходными же яв-
ляются физические модели явлений природы как случайных событий. Перенос излучения в атмосфере моделируется как движение частиц - фотонов. В начале каждой траектории фотон стартует с верхней границы атмосферы с координатами z' = z00 (г' = 0), т( = 770, <р' = <ро.
Свободный пробег фотона в плоской атмосфере в единицах оптической глубины распределен с плотностью вероятности р(Дг') = ^ . Это дает модель Дт' = —т\ 1п /3. Здесь и далее /3 обозначена равномерно распределенная на интервале [0,1] случайная величина, причем при каждом появлении в формулах подразумевается ее новое значение. В силу однозначной связи высоты z и оптической глубины т, найти высоту фотона после пробега несложно (см. подробнее [8]). В сферическом случае имеем просто р(Дт') = ехр(—Дг') и
Дт' = — 1п /3, но здесь Дг' - оптический путь вдоль конкретной трассы распространения излу-(
чения: Дг' = f a(z(l'))dl'. Это соотношение позволяет численно по Дг' получать I, а по ней из
о
(4) - z (задача сводится к решению методом Ньютона уравнения, конкретный вид которого определяется квадратурной формулой и способом интерполяции значений объемного коэффициента ослабления по высотной сетке). После моделирования свободного пробега фотон может либо вылететь через верхнюю границу атмосферы (z > Zoo), тогда его траектория заг кончена, либо попасть на поверхность (z' < zo) - тогда следует моделировать взаимодействие с ней, либо остаться в атмосфере (zq < z' < Zoo) и взаимодействовать с ней.
При взаимодействии с атмосферой если /3 > ljo(z'), то фотон поглощен и его траектория закончена. Иначе, для стандартно используемой в ММК таблично заданной индикатрисы рассеяния x(xí), г = 1,..., К , находится сначала косинус угла рассеяния у (см- подробнее [1]):
, x(z', Xk) - Vx4z>,Xb) + 2ДФ')(2Ь - Xk(z')) X = Xk +-д^-, Ь = 0, (12)
71-1
где Xn(z') - £ ,Xm+i) + x(z',Xm))(Xm - Xrn+i), n = 1 номер к в (12) опреде-
ли 1
ляется из условия Xk{z) < 26 < Xk+i{z')t a &k{z') = '^^-xj+i > и далее ~ новое направление движения фотона (r¡", ip")
tj" = tj'x + л/С1 - fo')2)(l -X2)cos ф,
и
— ю + sign (ir — ф) arceos I — X V V -
Vv^-fo'^Xi-fo")2),
(ф = 27Г/3 - случайный азимут рассеяния), после чего моделируется свободный пробег с новыми координатами.
При взаимодействии с поверхностью роль вероятности отражения играет альбедо (7), поэтому если /3 > А(г)', ip'), то происходит поглощение - конец траектории фотона, иначе необходимо моделировать направление отражения. В общем случае такал модель может быть получена аналогично (12) путем табуляции функции отражения (6), однако для двух стандартных рассмотренных выше моделей поверхностей все значительно проще. Для идеальной зеркальной поверхности направление отражения детерминировано: т]" = —т]', = <р'. Для изотропной поверхности все направления равновероятны, что дает г)" = — cos (f/З), ip" = 27Г/3. Далее следует моделирование свободного пробега фотона с полученными координатами (г/", <р") и высотой z = zо (г' = го).
Рассмотренный алгоритм моделирования траектории фотона сам по себе еще не позволяет получать искомые интенсивности 7(z, 770, ?7, у>о> ¥>)• Для этого вводятся специальные переменные - «счетчики» [1, 8] с нулевыми значениями в начале траектории. Перед моделированием каждого очередного свободного пробега фотона, находящегося в атмосфере, в счетчик (см. подробно [8]) прибавляется в плоском случае следующая величина:
Mwp) = ^êà1 exp(-(x(r) - r')h) (e(r})Q(z' - z) + (1 - в(т?))(1 - Q(z' - z))),
x = m + \/(i - i?2)(i - Ы)2) cos- (13)
Для фотона, после отражения с поверхности, опять же для плоского случая имеем
" - г0)/7?)(1 - в(т,)). (14)
Соответствующие выражения для сферического случая будут приведены ниже. Для изотропной поверхности из (6), (7), (9) независимо от {г]',<р') получаем ip(z, т), tp) — ^ ехр((го — t(z))/t]){ 1 — 0(77)). Для идеальной зеркальной поверхности направление отражения детерминировано, поэтому запись в счетчик необходимо производить в каждом случае рассеяния, учитывая попадание в него при т) < 0 излучения через отражающую поверхность, что дает вместо (14) прибавление к (13) величины
V, Ч>) = X(Z¿f] ехр((г0 - т')/г, )г(-т,) ехр((т0 - r(z))/r,)( 1 - 6(77)),
х" = -W + V/(1-Í72)(1-(^)2)cos(^ - (15)
где r{—T¡) - коэффициент зеркального отражения (8).
После моделирования достаточного числа траекторий искомое значение интенсивности находится как математическое ожидание (среднее арифметическое) записанных в счетчик на каждой траектории величин. Для перевода в абсолютные единицы оно умножается на поток на верхней границе атмосферы: Forjo в плоском случае и Fo - в сферическом. Помимо математического ожидания может быть оценена и случайная среднеквадратическая погрешность интенсивности (путем подсчета сумм квадратов записываемых в счетчик величин - см. подробнее [8]). Это дает возможность контроля точности расчетов, в частности автоматического выбора числа траекторий фотонов [8].
Общей идеей минимизации случайной погрешности расчетов по ММК является уменьшение разброса записываемых в счетчики значений. Для этого вводится дополнительный параметр - вес фотона w'. В начале траектории w — 1. При каждом акте взаимодействия с атмосферой вес пересчитывается как w" = u>o(z')w' и моделирование типа взаимодействия (поглощение-рассеяние) не производится, так как в w" уже записана доля рассеянного излучения, и затем сразу находится направление рассеяния и рассматривается дальнейший свободный пробег с весом w' = w". Аналогично при взаимодействии с поверхностью пересчитываем вес фотона w" = А(т}', ip')w' и проводим сразу моделирование направления отражения с w' — w". При записи в счетчики теперь величина ф(г, 77, (р) умножается на текущее значение веса к/.
При движении фотона к верхней границе атмосферы доля остающегося в ней излучения также известна, потому пересчитаем вес в плоском случае как w" = w' ехр(—T(z)/\r)'\) и будем моделировать свободный пробег фотона так, чтобы он неизбежно остался в атмосфере. Это приводит к модели Ат' — ?)'1п(1 — /3(1 — ехр(—т(г')/\т}'\))) (для сферического случая w" = w' exp (—т') и А г = 1п(1 — /?(1 — ехр(—т'))), где т' - оптический путь по текущей трассе до границы атмосферы). Эта схема известна как моделирование «без вылета». Она требует искусственного обрыва траектории фотона, когда его вес становится меньше заданной величины [8]. Данная схема применима и при движении фотона к нижней границе атмосферы, достаточно лишь в формулах плоского случая заменить r(z') на то — t(z') (а для сферического случая вид формул не меняется). Однако, в отличие от движения к верхней границе, здесь надо учесть долю излучения, достигающего поверхности. Технически это делается организацией «очереди отражения» - сохранением (в памяти компьютера) параметров попадающих на
поверхность фотонов (с весами 1 — ю"). После окончания моделирования траектории оставшейся части фотона в атмосфере извлекается фотон из очереди и начинается новая часть траектории (с отражения от поверхности). Поскольку вес фотона при каждом взаимодействии только уменьшается, исчерпание очереди неизбежно. Удобнее всего ее строить по принципу стека («последним пришел - первым ушел»), извлекая всегда последний записанный в очередь фотон.
Выше рассматривалась схема переноса фотонов в атмосфере, соответствующая физике процесса - от солнца к приемнику излучения. Недостатком ее является малая вероятность попадания фотона в счетчики в определенных случаях (например, для визирования верхней полусферы на больших высотах). Он проявляется и в сферической геометрии, где определить возможность попадания в счетчик - достаточно сложная задача, в отличие от формулы плоского случая (13). Можно, однако, считая, что фотон уже достиг приемника, проследить его обратную траекторию - к источнику излучения. Это логическое построение приводит к схеме «моделирования от приемника» (иначе - по обратным траекториям). Ее детальное рассмотрение с учетом связи алгоритмов ММК с исходным уравнением переноса (см. [8, 12]) дает простой результат - для реализации схемы «от приемника» достаточно всего лишь формально поменять координаты источника (начала траектории) и приемника излучения, ничего более не изменяя в описанном алгоритме ММК. Это является общим следствием обратимости оптических явлений [1]. Заметим, что такая простая реализация схемы «от приемника» возможна лишь для неполяризованного излучения.
Для окончания описания алгоритма ММК осталось лишь привести выражения для функций записи в счетчики -ф(г)т],'р) в сферическом случае, учитывая использование схемы «от приемника» и пересчет координат солнца щ и <р'0 по формулам (4). Вместо (13) имеем
Ф(*,Г},Ч>) = ^Роо(^о), X = Мо + л/(1 — (ч#)2)(1 — (Чо)2) С08(у># - <р'о),
где Роо(-г', т}'о) ~ функция пропускания до верхней границы атмосферы от текущей высоты г с направлением —щ; запись производится, если фотон не находится в области тени, что дает условия щ > 0 или г)'0 < 0, но (Я + — (г)'0)2 > Я + го- При отражении от поверхности
выражение (14) заменяем на
при г)'0 > 0. Для идеальной зеркальной поверхности вместо (15) получаем
X" = « + ч/(1-(^)2)(1-(^)2)со8(^ - <р'о).
Здесь Po(z',t)q) - функция пропускания до поверхности от высоты z с направлением %, необходимым для отражения косинусом зенитного угла солнца в точке рассеяния; щ - этот же угол, но на поверхности. Смысл этих величин и метод их нахождения рассматривались выше при описании приближения однократного рассеяния (10). Запись в счетчик производится, если в (10) rj'o < —г)"(1).
Реализация алгоритма в виде кода SCATRD, расчет производных. Описанный алгоритм реализован в виде компьютерного кода SCATRD (scattering radiation with derivatives) (его ранняя версия, общие принципы и особенности подробно изложены в [6]). Отметим универсальность кода: расчеты по нему могут вестись в любой геометрии визирования и переноса излучения (плоской, сферической и смешанной, когда ее выбор определяется зенитным углом конкретной трассы).
Важное достоинство кода БСАТГШ - вычисление производных по любым параметрам входных моделей атмосферы и поверхности. Для ММК оно сводится к дифференцированию ряда Неймана (11), для плоского случал - подробно с приведением практических вычислительных формул и алгоритмов, описанного в [8]. Учет сферичности не вносит в технику дифференцирования какой-либо специфики.
Тестирование кода ЗСАТШЭ. Задачей тестирования компьютерных кодов, предназначенных для научных расчетов, является проверка соответствия результатов их работы некоторым ожидаемым значениям. Можно выделить четыре этапа тестирования: микротестирование на уровне отдельных блоков кода, техническое тестирование - проверка работы кода исходя из заявленных в нем возможностей (например, корректности зависимости результатов расчетов от параметров), проверка по системе физических тестов, сравнение с результатами независимых расчетов. Рассмотрим подробнее два последних этапа.
Физические тесты - это определенным образом сформулированные условия, которым должны подчиняться результаты вычислений. В качестве примера подобной системы тестов укажем на б тестов, предложенных в [14] для расчетов оптических характеристик аэрозольных частиц. Для задач переноса излучения в атмосфере аналогичные тесты не сформулированы (в том числе и в фундаментальной работе [13]), хотя многие из них достаточно очевидны.
Первый тест связан с появляющейся симметрией вращения при расположении солнца в зените.
Тест 1. При щ = 1 (солнце в зените) интенсивность излучения (в любой геометрии визирования) не зависит от азимута визирования.
Следующие три теста основаны на очевидном анализе зависимости интенсивности от параметров атмосферы и поверхности.
Тест 2. При индивидуальном увеличении объемного коэффициента поглощения (всей атмосферы или любого ее уровня) интенсивность может только уменьшаться. В терминах производных: производная от интенсивности рассеянного солнечного излучения по объемному коэффициенту поглощения не может быть положительной.
Тест 3. Для модели изотропной отражающей поверхности при индивидуальном увеличении ее альбедо интенсивность может только возрастать. В терминах производных: производная от интенсивности по альбедо поверхности в случае изотропного отражения не может быть отрицательной.
Тест 4■ При стремлении к нулю (последовательном уменьшении) объемного коэффициента рассеяния на всех уровнях атмосферы интенсивность должна асимптотически стремиться к интенсивности однократно рассеянного излучения.
Пятый тест основан на соотношении (9).
Тест 5. Для изотропной отражающей поверхности с достаточно большим альбедо А в плоской модели атмосферы при визировании поверхности (с любых высот под любыми углами), при стремлении к нулю (последовательном уменьшении) объемного коэффициента ослабления на всех уровнях атмосферы интенсивность должна асимптотически стремиться к значению ^АРощ.
Шестой тест базируется на известных соотношениях симметрии для коэффициентов яркости атмосферы (см., например, [1]).
Тест 6. Для плоской модели атмосферы с изотропной поверхностью (с любым альбедо) при визировании с верхней границы в направлении поверхности (77 < 0) при взаимной замене косинусов зенитных углов солнца и визирования должно выполняться
соотношение
VoI{zoo,V> Vo,<p) = -vl(zoo,-r]o,-ri,v>).
Седьмой тест основан на возможности аналитического вычисления интенсивности однократно рассеянного излучения в плоской однородной атмосфере и может быть использован для тестирования расчетов в приближении однократного рассеяния, а с учетом теста 4 - и как асимптотический тест общего случая.
Тест 7. В плоской однородной атмосфере с оптической толщиной то, вероятностью выживания кванта uiq, индикатрисой ж(х), поверхностью с нулевым альбедо (А = 0) при <ро = 0 интенсивность однократно рассеянного излучения равна
' Ч'ММ) = F^ZX[X0) если Г, > 0,
если г) < 0.
Выполнение всех приведенных тестов для кода SCATRD было подтверждено многочисленными расчетами при различных параметрах визирования, атмосферы и поверхности.
Важным этапом тестирования является сравнение с результатами вычислений, независимо выполненных другими авторами. Для плоской геометрии расчеты по коду SCATRD сравнивались с данными, приведенными в [13] (раздел 7.1), при всех вариациях точки визирования и параметров атмосферы. Среднее (по модулю) отклонение расчетов по коду SCATRD от «точных» (термин и кавычки из [13]) результатов составляет 3,6% при максимальном 17,9%. Аналогичные показатели для приведенных в [13] и выполненных по ММК расчетов - 3,0 и 14,2%. Для сферической геометрии переноса излучения расчеты по коду SCATRD сравнивались с расчетами по коду SCIATRAN [4, 5], который принимал участие в кампании сравнения шести кодов [2] и в этом смысле считался эталонным. Здесь среднее отклонение результатов составило 3,4%, что укладывается в общий интервал разброса кодов, приведенный в [2]: 2-4%. Таким образом, сравнение с результатами независимых расчетов следует признать хорошим, а в качестве общей оценки отклонения результатов кода SCATRD от других аналогичных радиационных кодов может быть взято значение 3,5%.
Заключение. Компьютерный код SCATRD, в котором реализован описанный в настоящей работе вычислительный алгоритм, успешно прошел процедуру тестирования и может быть использован для расчетов интенсивности рассеянного солнечного излучения и производных от нее по параметрам атмосферы и поверхности в различных задачах моделирования оптических измерений в атмосфере. Планируется дальнейшее развитие и совершенствование кода.
Автор благодарит сотрудников Института физики окружающей среды (Institute of Environmental Physics) университета г. Бремена, во главе с проф. Дж. Борроусом (J. P. Burrows), и А. В. Розанова за предоставление для сравнения расчетов компьютерного кода SCIATRAN.
Summary
Vasilyev A. V. Numerical simulation of solar multiple scattering radiation intensity and its derivatives in spherical geometry atmosphere (computer code SCATRD).
The physical-mathematical model of scattering radiative transfer in spherical geometry atmosphere is described. Multiple scattering and reflection of radiation form surface axe considered and two reflection models: ideal mirror and isotropic axe analyzed. Some formulas for single-scattering approximation and an algorithm for calculations of multiple scattering radiation by Monte-Carlo technique axe described in details. The main characteristics of a computer code realizing this model are considered. The system of some tests for examination of radiation code work are suggested.
Литература
1. Тимофеев Ю. M., Васильев А. В. Теоретические основы атмосферной оптики. СПб.,
2003. 2. Loughman R. P., Griffioen Е., Oikarinen L. et al. // J. Geophys. Res. 2004. Vol. 109, N D6, D06303 (doi:10.1029/2003JD003854). 3. Постыляков О. В. // Изв. РАН. Сер. Физика атмосферы и океана. 2004. Т. 40, № 3. С. 276-290. 4. Rozanov A., Rozanov V., Buchwitz М. et al. 11 Adv. in Space Research. 2005. Vol. 36, N 5. P. 1015-1019. 5. Rozanov A. Computer code SCIATRAN (User Guide, download account), http://www.iup.physik.unibremen.de/sciatran.
2004. 6. Васильев А. В., Поляков А. В. // Оптика атмосферы и океана. 2005. Т. 18, № 4. С. 352-358. 7. Тимофеев Ю. М., Поляков А. В. Математические аспекты решения обратных задач атмосферной оптики. СПб., 2001. 8. Васильев А. В., Мельникова И. Н. Коротковолновое солнечное излучение в атмосфере Земли. Расчеты. Измерения. Интерпретация. СПб., 2002. 9. Васильев А. В. // Оптика атмосферы и океана. 2004. Т. 17, № 4. С. 334-338. 10. Колмогоров А. Н., Фомин С. В. Элементы теории функций и функционального анализа. М., 1989. 11. Васильев А. В. // Оптика атмосферы и океана. 2003. Т. 16, № 8. С. 725-729. 12. Метод Монте-Карло в атмосферной оптике / Под ред. Г. И. Марчука. Новосибирск, 1976. 13. Перенос радиации в рассеивающих и поглощающих атмосферах: стандартные методы расчета / Под ред. Ж. Ленобль; Пер. с англ. Ж. К. Золотовой. Л., 1990. 14. Борен К. Ф., Хаф-мен Д. Р. Поглощение и рассеяние света малыми частицами/ Пер. с англ. 3. И. Фейдулина и др. М., 1986.
Статья поступила в редакцию 6 декабря 2005 г.
