Ряды экспонент в расчетах переноса излучения методом Монте-Карло в пространственно неоднородных аэрозольно-газовых средах

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

Том 7, № 5, 2002

РЯДЫ ЭКСПОНЕНТ В РАСЧЕТАХ ПЕРЕНОСА ИЗЛУЧЕНИЯ МЕТОДОМ МОНТЕ-КАРЛО В ПРОСТРАНСТВЕННО НЕОДНОРОДНЫХ АЭРОЗОЛЬНО-ГАЗОВЫХ СРЕДАХ *

К. М. Фирсов, Т. Ю. Чеснокова, В. В. Белов,

А. Б. Серебренников, Ю. Н. Пономарев Институт оптики атмосферы СО РАН, Томск, Россия

e-mail: fkm@iao.ru

A method for consideration of the molecular absorption when solving the radiation transfer equation in spatially inhomogeneous aerosol and gaseous media by the Monte-Carlo method is considered. The method is based on the expansion of the broadband transmission function in a series of exponents. The applicability of series of exponents to the spatially inhomogeneous atmosphere is analyzed both theoretically and numerically. The results of numerical modeling have demonstrated a high efficiency of the given method: the computational time decreased more than 102 times compared to the line-by-line method, and the discrepancy of the results obtained was less than 0.2 %.

В последние десятилетия достаточно широкий круг известных специалистов по атмосферной оптике обращается к решению проблемы распространения немонохроматического излучения в многокомпонентной среде со сложными и изменяющимися в пространстве спектрами поглощения отдельных ее компонентов [1-4]. Одна из причин этого связана с необходимостью более точного учета энергетического взаимодействия солнечного оптического излучения с атмосферой Земли в моделях прогноза погоды и климата. Это можно осуществить, уточняя прежде всего информацию о тонкой структуре спектров поглощения газовой составляющей атмосферы (так как спектральные характеристики аэрозольной существенно менее изменчивы и не требуют дополнительных усилий для их детализации). Если эти уточнения осуществлены, то цель достигается решением стационарного уравнения переноса для каждой длины волны излучения с последующим интегрированием по заданному спектральному диапазону. Этот метод получения решения назван line by line. К сожалению, этот простой прием, хотя и дает асимптотически точное решение, оказывается чрезвычайно трудоемким даже для современных вычислительных средств. Последнее обстоятельство стимулировало разработку нового метода учета огромных объемов спектроскопической информации, который основывается на разложении широкополосной функции пропускания в ряд экспонент. Это обеспечивает высокую эффективность расчета интегральных по спектру характеристик излучения, прошедшего пространственно-однородную

* Работа частично поддержана Российским фондом фундаментальных исследований, гранты №02-0564492 и №01-05-61152.

© К. М. Фирсов, Т. Ю. Чеснокова, В. В. Белов, А. Б. Серебренников, Ю.Н. Пономарев, 2002.

многокомпонентную аэрозольно-газовую среду. Однако для пространственно неоднородной среды последовательный анализ приведен лишь для молекулярной нерассеивающей атмосферы [1-3]. Целью данной работы являлось исследование вопроса применимости рядов экспонент для многокомпонентной пространственно-неоднородной атмосферы.

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

г

X

Здесь г(х,р) — искомая функция, которую определяют как плотность столкновений, связанную с интенсивностью излучения соотношением

I

РехН' ’ ^ )

где х = (т, П), х — фазовое пространство координат и направлений; V — волновое число;

— секанс зенитного угла солнца; 5© — солнечная постоянная; /Зех* — коэффициент ослабления. Обобщенное ядро описывается выражением [5]

к

где "(т', V) — альбедо однократного рассеяния; /(т', П, П',*) — индикатриса рассеяния; т3(г, т',*) — оптическая толща.

Плотность распределения источников Ф(х, V) для стандартной задачи равна

Ф

Решение интегрального уравнения (1) представимо в виде ряда Неймана

г =

г=1

ЕЛ'4-

Вопросы, связанные с существованием решения, обсуждаются в работе [5]. Для более компактной записи решения в (3) и (4) выделим сомножители Ф0(ж,^) и к0(х, х', V), в которые не входят характеристики молекулярного поглощения

к(х', х, V) = вехь(г, V)-—1—- ехр{-т(г, г', V)}ко(х', х, V),

РехЦг', V)

Ф(х, V) = вехі(г, V) ехр{ т(г, гте, V)}Фо(х, V),

7 / ' \ П / ' \/(т', П, П', V) ехр{ —Га<т, т', V)} г (т — т' \ , ,

где ко(х, x'v) = яс^^)------------|т—т'|2-------------* — утзту); =

ехр{—та<т, тто, v)}*<П — По); — коэффициент аэрозольного рассеяния; т, та — оптичес-

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

Рассмотрим два первых члена ряда (5):

КФ = у к(х1, х, v)Ф(x1, v)dx1 =

X

вех^х^у ко(хі, х^)Фо(хі^) х ехр{-т(гі, гто^) - т(г, Гі, V)}^хі,

X

К2Ф = ^ у>к(х1, х^)к(х2, х^ v)Ф(x2, v)dx1dx2 =

X X

= вехі(х^)^ J ко(хі, х, V)ко(х2, хі, V)Фо(х2, V) х

X X

х ехр{-т(г2, гто, V) — т(г, г1, V) — т(г1, г2, v)dx1dx2.

По аналогии можно получить выражения и для других членов ряда. Для расчета интегрального по спектру значения интенсивности излучения полученное решение усредним по интервалу — Дv/2' ^ + Дv/2]. Величину интервала Дv будем полагать такой, чтобы внутри него содержалось большое количество спектральных линий, однако солнечную постоянную, коэффициент рассеяния и индикатрису рассеяния можно было бы считать неизменными. На рис. 1 приведен типичный спектр поглощения водяного пара в ближней инфракрасной области. Ширина спектрального интервала составляет величину 100 см-1, которая удовлетворяет сформулированным выше требованиям. Меняя порядок интегрирования и вынося из-под знака интеграла по V все сомножители, которые не содержат характеристики молекулярного поглощения, можно получить выражение

Рис. 1. Спектр поглощения водяного пара.

vo+Av

IAv(r, П) = J I(r, Q,v)dv =

vo—A v

vo+Av

= Av|^o|So(vo^ ko(xi, ж,^о)Фо(ж1,^о) x J exp{-r(ri, r^,v) - r(r, rbv)}dvdxi +

X v0+Av

+Av |^o|So (v) J j ko(Xl, X,Vo)ko(X2, ЖЬ ^о)Фо(Ж2, Vo)x

X X

vo+Av

x Av J exp{-T

vo +Av

В формуле (7) в явном виде появляется функция пропускания, обусловленная молекулярным поглощением:

т=Av / exp{-Tj

Av

где Tj(v) — оптическая толща j-го члена ряда (7). В частности, T2(v) имеет вид

T2

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

ного поглощения соотношением

Tj (v) = eabs(v)LJ ,

где Lj — суммарная длина трассы. Экспоненциальная зависимость пропускания от длины трассы позволяет при помощи преобразования Лапласа осуществить переход в пространство кумулятивных волновых чисел g [6], где функция пропускания имеет вид

oi

Tj=Av / exp {-e>

o

Численная реализация этого преобразования достаточно простая: необходимо рассчитать массив {ваы (v)}Ni с постоянным шагом по v, а затем провести его сортировку по возрастанию. Коэффициент Д”ы(g) можно интерпретировать как коэффициент поглощения в пространстве кумулятивных волновых чисел g. Действительно, область значений функции em>s(v) такая же, как у функции em>s(g), обе функции — непрерывные. Для каждого значения em>s(v) в пространстве волновых чисел v существует равный ему коэффициент embs(g) в пространстве кумулятивных волновых чисел (обратное утверждение также справедливо). Существенное отличие состоит в том, что em>s(g) — монотонно возрастающая функция аргумента g, а не быстроосциллирующая, как embs(v). Применяя к (10) квадратурные формулы, легко получить короткий ряд экспонент (5-10 членов ряда):

П

Tj = £ Ci exp {-ffis

i=i

(С, — коэффициенты и узлы соответствующих квадратурных формул). Наиболее про-

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

Описанный выше прием позволяет рассчитать функцию пропускания с высокой точностью при использовании небольшого набора коэффициентов поглощения на неко-

торых специально выбранных волновых числах. Переход от (8) к (10) является интегральным преобразованием, т. е. нет явной связи между монохроматическими коэффициентами АтЫ (V) и ваЬв(^), поэтому определить эти специально выбранные частоты возможно лишь после того, как рассчитан спектр ) и преобразован в Д”^#).

Для того чтобы осуществить описанное выше интегральное преобразование, в случае неоднородной трассы можно формально ввести параметр £ в функцию пропускания, описываемую формулой (8):

т = ¿V / ехР{-т)^) х

осуществить преобразование, а затем положить £ равным единице.

Однако для неоднородной трассы есть существенное отличие, связанное с тем, что соотношение (9) в пространстве кумулятивных частот в общем случае не выполняется, так как каждому значению т) (V) можно поставить в соответствие лишь т) (#). Причем положение т) (#) в пространстве кумулятивных волновых чисел будет зависеть от ], которое определяет трассу, т. е. соотношения

Т1 Ы = т(Г1, гте,&) + т(г, Г1 ,&),

Т2Ы = т

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

Исследования, проведенные в работах [1-3], показали, что наиболее сильно эти соотношения нарушаются, если имеются большие градиенты температуры. Однако в атмосфере Земли перепад температур не превышает 100 К. Поэтому можно полагать, что формулы (12) являются достаточно хорошим приближением и оптическая толща в пространстве кумулятивных частот связана с коэффициентом поглощения так же, как в пространстве волновых чисел

т

г

Данное приближение получило название к-корреляции и неоднократно проверялось для молекулярной нерассеивающей атмосферы в работах [1-3], где было показано, что погрешности расчета функции пропускания не превышали 1 %. Используя данное приближение, в формуле (7) можно заменить интеграл по V на ряд типа (11). С учетом этого для каждого члена ряда получается выражение, подобное (6). После очевидных преобразований легко получить систему уравнений

,(ж,^) = у к(ж', ж,^),(ж',^)^ж' + Ф*(ж,^),

X

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

N

Zi

ід*(r, П, Vo) = |^o|so(vo) ^2 Ci

es“cat(r,V0) + embs(r,gi)’

i=1

где

ki(x', x, Vo) =

exp{-T(r, r', gi)}ko(x', x, Vo),

^i(x) V0) = [es“cat(r,Vo) + embs(r,gi)] exp{ T(r, гте; gi)}^o(x, Vo)•

Этот прием позволит решать уравнение переноса любым методом для небольшого числа заданных волновых чисел д^, а затем провести суммирование решений с весами С^. Следует также отметить, что условие постоянства солнечной постоянной в пределах ¿V не является обязательным. Для преодоления данного ограничения функцию пропускания можно переопределить [7]

и представить в виде ряда (10). Таким образом, для широких спектральных интервалов Av, величина которых может достигать 100 см-1 и более, достаточно точное решение уравнения переноса может быть получено при использовании 5-6 членов ряда.

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

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

Ограничимся рассмотрением наиболее простой модельной схемы, требующей относительно небольших ресурсов ПЭВМ на базе процессора Intel Pentium Celeron 333, чтобы выяснить принципиальную возможность включения изложенного выше метода “ k-распределения” в алгоритмы прямого моделирования распространения нерассеянных и диффузных световых потоков методом Монте-Карло. Кроме того, важно оценить степень эффективности новых решений по сравнению с традиционными (line by line) именно в упрощенной ситуации. Это важно для того, чтобы сделать вывод о целесообразности исследования границ их применимости в более широкой области как оптико-геометрических аспектов постановок задач, так и различных модификаций метода Монте-Карло, используемых при решении более сложных атмосферно-оптических задач.

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

[9]):

0(0) = ер1 040) + 7^ 0т(0),

7sct 7sct

где 7sct, 7m — значения коэффициентов аэрозольного и молекулярного рассеяний в j-м слое, 7sct = 7sct + emt; g(0) — “взвешенная”, g°(0) — аэрозольная, gm(0) — молекулярная индикатрисы рассеяния.

При проведении численных расчетов атмосфера разбивалась на оптически однородные слои. От 0 до 25 км толщина слоев полагалась равной 1 км, в диапазоне высот 25 ... 50 км — 5 км, а верхняя атмосфера разбивалась на два слоя с границами 50... 70 и 70... 100 км. Полагалось, что основной (поглощающей) молекулярной компонентой является водяной пар. Требовалось оценить освещенность земной поверхности в спектральном интервале 0.943 ...0.952 мкм (или 10 500... 10 600 см-1).

При выборе модели оптических свойств атмосферы мы стремились воспроизвести си-

Рис. 2. Геометрическая схема эксперимента.

О 20 40 60 80 100 h, КМ

Рис. 3. Вертикальный профиль коэффициента аэрозольного ослабления.

туацию, не обязательно близкую к реальной, а такую, чтобы рассеянное излучение составляло заметную долю в полной освещенности и, кроме того, позволяло осуществить контрольные расчеты line by line на ПЭВМ класса Intel Pentium Celeron 333 за практически реализуемый промежуток времени, свободный в том числе и от сбоев компьютера. Поэтому метеорологическая дальность видимости в приземном слое аэрозольной атмосферы была принята равной 8 км. На рис. 3 в качестве примера приведен высотный профиль коэффициентов аэрозольного ослабления, а на рис. 1 показана спектральная зависимость коэффициентов поглощения водяным паром в приземном слое атмосферы. Если характеризовать выбранную модель интегральными оптическими характеристиками, то они соответствовали следующим значениям: оптическая толщина, обусловленная аэрозольным рассеянием, та = 1.42, молекулярным рассеянием — тт = 0.0107, тогда как оптическая толща, обусловленная молекулярным поглощением в данном спектральном интервале, менялась от 0.22 до 250.

Численная имитация методом Монте-Карло процесса распространения излучения в молекулярно-аэрозольной атмосфере осуществлена прямым моделированием. Этот выбор обусловлен простотой метода и его невысокой трудоемкостью при расчетах освещенности земной поверхности на заданной длине волны Л [5].

Тестовые расчеты по методике line by line для рассмотренной выше модели атмосферы показали, что среднее время ПЭВМ, затрачиваемое на моделирование освещенности Ej (где j — индекс длины волны из интервала 0.943... 0.952 мкм), составило 13.32 с при среднеквадратичной ошибке вычисления е ~ 0.1 %. Для корректного учета спектральной зависимости коэффициентов поглощения водяного пара в области 1 мкм требуется проведение расчетов с шагом по Л = 10-2 см-1 (сопоставимым с полушириной спектральной линии), т. е. выполнение 10 000 вариантов расчетов методом Монте-Карло. При этом интегральная освещенность E рассчитывается по формуле

где N = 10 000; = 0.01 см-1; Б„ — значение солнечной постоянной, а Е^ — доля потока

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

N

j=0

£,% 1.2

Рис. 4. Погрешность расчета освещенности земной поверхности в зависимости от числа гауссовских квадратур.

Отметим, что на моделирование 10 000 значений E потребовалось 37 ч времени ПЭВМ. Расчет освещенности с использованием рядов экспонент для учета поглощения молекулами Н2О излучения в том же диапазоне длин волн 0.943 ... 0.952 мкм потребовал проведения не более 10 оценок Ej методом Монте-Карло для специально подобранного набора длин волн и суммирования по формуле

П

E = AVJ2Cj Ej Sx,,

j = 1

где Cj — коэффициенты гауссовских квадратур. Машинное время на получение этого решения при n =10 (вместе с моделированием методом Монте-Карло) не превысило 5 мин. При этом относительная погрешность оценки (12) по сравнению с эталонными расчетами методом line by line по формуле (11) не превысила значения 0.2 %. Были также проведены расчеты с различным числом членов ряда экспонент (рис. 4). Результаты моделирования показали, что для достижения высокой точности в рассматриваемом спектральном интервале достаточно 6-7 членов ряда при использовании гауссовских квадратур.

Сопоставление метода k-распределения и прямого расчета методом line by line для нерассеивающей атмосферы показало, что погрешности расчета были близки к тем, которые были получены для рассеивающей атмосферы. Насыщение, наблюдаемое на рис. 4 для n > 6, обусловлено погрешностью интегрирования по v при вычислении функции пропускания методом line by line.

В таблице приведено сопоставление расчетов потоков нисходящего излучения с использованием наших программ (метод k-распределения) и расчетами Б. А. Фомина. Из таблицы следуют два вывода. С одной стороны, наблюдается хорошее согласие с не-

Потоки нисходящего излучения в интервале 10 000... 10 500 см-1, рассчитанные методами k-распределения и line-by-line [10] при поглощении водяным паром для лета средних широт.

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

Высота, км Нисходящие потоки, Вт/м2 Относит. разность, %

Метод k-распределения line-by-line [10]

0 27.11 26.99 0.44

1 29.30 29.20 0.34

5 34.48 34.45 0.09

10 35.72 35.69 0.08

70 35.76 35.76 0.00

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

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

Таким образом, в данной постановке задачи переход к новой методике учета поглощения молекулярной составляющей атмосферы для расчета освещенности земной поверхности излучением в полосе 0.943-0.952 мкм привел к потере точности оценок в 0.2% при возрастании скорости их получения более чем в 350 раз. Этот результат дает основание считать, что рассмотренная методика параметризации спектроскопической информации о молекулярной компоненте атмосферы может перевести в класс практически решаемых проблему учета влияния рассеяния (включая многократное) в атмосфере на характеристики световых потоков в широких спектральных интервалах.

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

[1] Lacis A. A., Oinas V. A description of the K-distribution methods for modelling nongray gaseous absorption, thermal emission, and multiple scattering in vertically inhomogeneous atmospheres // J. Geoph. Research. 1991. Vol. 96, No. D5. P. 9027-9063.

[2] GoüDY R., West R., Chen L., Crisp D. The correlated-k method for radiation calculations in nonhomogeneous atmospheres // J. Quant. Spectrosc. and Radiat. Transfer. 1989. Vol. 42, No. 6. P. 539-550.

[3] Riviere Ph., Soufani A., Taine J. Correlated-k and fictious gas methods for H2O near 2.7 ^m // J. Quant. Spectrosc. and Radiat. Transfer. 1992. Vol. 48, No. 2. P. 187-203.

[4] ПЕРЕНОС радиации в рассеивающих и поглощающих атмосферах. Стандартные методы расчета / Под ред. Ж. Ленобль. Л.: Гидрометеоиздат, 1990. 262 с.

[5] МЕТОД Монте-Карло в атмосферной оптике / Под ред. Г. И. Марчука. Новосибирск: Наука, 1976. 285 c.

[6] ТВОРОГОВ С. Д. Некоторые аспекты задачи о представлении функции пропускания в ряд экспонент // Оптика атмосферы и океана. 1994. Т. 7, №3. C. 315-326.

[7] Firsov K.M., MITSEL A. A., Ponomarev Yu.N., Ptashnik I. V. Parametrization of transmittanse for application in atmospheric Optics // J. Quant. Spectr. and Radiat. Trasf. 1998. Vol. 59, No. 3-5. P. 203-213.

[8] Зуев В. Е., Белов В. В., Веретенников В. В. Теория систем в оптике дисперсных сред. Томск: Изд-во СО РАН, 1997. 402 с.

[9] Kneizys F. X. ET. AL. User Guide to LOWTRAN 7, AFGL-TR-86-01777. ERP N 1010 / Nanscom AFB, MA01731.

[10] Fomin B.A., Romanov S.V., Rublev A. N., Trotsenko A.N. Line by iline benchmark calculations of solar radiation transfer parameters in a scattering atmosphere. М., 1992 (Препр. ИАЭ. IAT 5525.1).

Поступила в редакцию 28 марта 2001 г., в переработанном виде —16 апреля 2002 г.