Восстановление характеристик подстилающей поверхности Сибирского региона по данным спектрорадиометра MODIS Текст научной статьи по специальности «Физика»

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

Том 11, часть 1, Специальный выпуск, 2006

ВОССТАНОВЛЕНИЕ ХАРАКТЕРИСТИК

ПОДСТИЛАЮЩЕЙ ПОВЕРХНОСТИ СИБИРСКОГО РЕГИОНА ПО ДАННЫМ СПЕКТРОРАДИОМЕТРА MODIS

A.A. Лагутин, Ю.А. Никулин, И. А. Шмаков, А. П. Жуков, Ал. А. Лагутин, А.Н. Резников, В. В. Синицин Алтайский государственный университет, Барнаул, Россия e-mail: lagutin@theory.asu.ru

In this paper we summarize the approaches and codes used at Altai State University for the atmospheric correction of visible to middle-infrared MODIS data and for retrieval of land parameters such as the bidirectional reflectance distribution function, albedo and vegetation index.

Введение

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

В силу важности данных об альбедо ПП для решения широкого спектра задач — от глобального прогнозирования погоды и моделирования обменных циклов энергией, водой, углеродом в различных типах ПП до обработки данных спутниковых и наземных детекторов при измерениях характеристик тропосферного аэрозоля — научным сообществом сформулированы требования к частоте обновления и точности этих данных [1]. Сегодня считается, что, например, для моделей общей циркуляции атмосферы необходимо ежемесячное обновление данных, причем абсолютная точность определения коротковолнового альбедо должна быть ±0.02 [1]. В то же время при восстановлении аэрозольной оптической толщины (AOT) атмосферы по спутниковым данным используются сферическое альбедо и двунаправленный коэффициент спектральной яркости (ДКСЯ) подстилающей

поверхности для состоявшейся геометрии съемки в момент ее проведения. Требование к

±

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

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

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

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

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

Необходимость проведения атмосферной коррекции данных спутникового прибора для каждого свободного от облачности пикселя приводит к тому, что спектрорадиометр, используемый для мониторинга характеристик ПП, должен также обладать возможностями восстановления по полученным в процессе съемки сцены интенсивностям характеристик основных быстро меняющихся компонентов атмосферы (прежде всего — аэрозольной оптической толщины атмосферы и общего количества водяного пара). Так, например, установленный на метеорологических спутниках системы NOAA спектрорадиометр AVHRR (Advanced Very High Resolution Radiometer) не позволяет проводить восстановление характеристик подстилающей поверхности. Одна из причин этого — невозможность получения оценки прежде всего аэрозольной оптической толщины для анализируемого пикселя сцены.

Только с выводом на полярную солнечно-синхронную орбиту спутника Terra/NASA в декабре 1999 г. со спектрорадиометром MODIS (MODerate resolution Imaging Spectroradi-ometer) па борту в рамках программы EOS (Earth Observing System)/NASA [3] с каналами, согласованными с предварительно разработанными алгоритмами восстановления ключевых элементов системы атмосфера — подстилающая поверхность Земли, появилась возможность операционного восстановления характеристик ПП. 36 каналов MODIS покрывают диапазон спектра электромагнитного излучения от 0.42 до 14.24мкм [3]. По данным 20 каналов строится маска облачности [4]. Пять каналов MODIS (2, 5, 17-19) из ближнего инфракрасного диапазона используются для определения общего количества водяного пара в атмосфере в зонах, свободных от облачности [5]. Каналы 1, 3, 7 являются ключевыми при построении оценки АОТ [2,6].

Эти уникальные возможности MODIS позволили впервые осуществлять атмосферную коррекцию в операционном режиме [7], проводить восстановление двунаправленного коэффициента спектральной яркости и альбедо подстилающей поверхности с разрешением 1 км за 16-дневный период, решать другие задачи, принципиальные для понимания функционирования системы Земля. Восстановление ДКСЯ и альбедо проводится в предназначенных для этого первых семи каналах MODIS с центрами на длинах волн 647 нм (канал 1), 858нм (канал 2), 470нм (канал 3), 555нм (канал 4), 1240нм (канал 5), 1640нм (канал 6), 2130 нм (канал 7), а также в видимом (Л = 0.3... 0.7 мкм), ближнем инфракрасном (Л = 0.7 ... 5.0 мкм) и во всем коротковолновом (Л = 0.3 ... 5.0 мкм) диапазонах длин волн [8].

В настоящей статье излагаются подходы и вычислительные технологии, используемые в Центре космического мониторинга Алтайского государственного университета при проведении атмосферной коррекции данных .\10DIS Теггн и восстановлении характеристик ПП, представлены некоторые результаты измерений двунаправленного коэффициента спектральной яркости и альбедо подстилающей поверхности Сибирского региона. Сырой поток МСШК, являющийся информационной основой представленных в работе результатов, принимался станцией ЕОСкан в Барнауле (53°21/ с.ш., 83°47/ в.д.) в режиме реального времени.

1. Коэффициент спектральной яркости подстилающей поверхности

Интенсивность солнечного излучения L(À, в, ê, p) [Вт/(м2-ср-мкм)], регистрируемая сиект-рорадиометром спутника на верхней границе атмосферы в случае ламбертовой подстилающей поверхности, может быть определена так (см, например, [9,10]):

L А, в, ê, p) = L0 А, в, ê, p) +-----. 1

п [1 - ps(À, в, ê, p)s(À)

Здесь L0(À, в, ê, ф) — вклад в показание спектрорадпометра, обусловленный отражением солнечного излучения от полубесконечной атмосферы; TÀ), T^ (pv, À) — функции про-

À

ли и поверхность Земли — спектрорадиометр спутника соответственно; F0(À)ps/n — интенсивность солнечного излучения, падающая па верхнюю границу атмосферы; ps(À, в, ê, p) — коэффициент спектральной яркости подстилающей поверхности, когда атмосфера над ней отсутствует; s(À) — сферическое альбедо атмосферы в случае, когда изотропное излуче-

в

зенитный угол Солнца (ps = cos в), ê — зенитный угол спутника (pv = cos ê), p — азимутальный угол спутника относительно Солнца, Поскольку MODIS проводит измерение интенсивности солнечного излучения, падающего на верхнюю границу атмосферы, в (1) удобно перейти к безразмерным величинам

р{А, в, ê, p) = —----, ро(А, в, ê, tp) = —----. (2)

Fo^s/п Fo^s/п

Тогда для КСЯ системы атмосфера — подстилающая поверхность Земли р получаем следующее уравнение:

л п л л п л . А)ТТ(^, А)ря(Л, в, ê, <р)

р А, в, ê, <р) = ро А, в, ê, <р) +-----. 3

1 - ps(À^,ê, p)s(À)

Хотя это уравнение справедливо лишь для однородной ламбертовой поверхности, на данном этапе ежедневного глобального мониторинга характеристик поверхности суши с использованием данных MODIS в рамках программы EOS оно положено в основу операционного алгоритма восстановления КСЯ подстилающей поверхности [7]. Из (3) видно,

что по измеренным спектрорадиометром МСШЕЗ для свободного от облачности пикселя (01, $1, ) в канале Л = ДА КСЯ системы атмосфера — подстилающая поверхность Земли р(Л, 01, $1, рI) можно получить ириближенную оценку рДЛ, 01, $I, рч), если говеетны ро, Т К Т^ и в(Л).

Ослабление солнечного излучения на трассе Солнце — поверхность Земли — спутник в оптическом и ближнем инфракрасном диапазонах обусловлено главным образом тремя процессами: молекулярным рассеянием, поглощением водяным паром, озоном, кислородом и другими газовыми компонентами атмосферы, аэрозольным рассеянием и поглощением. При заданных для рассматриваемого пикселя сцены ПП углах освещения и наблюдения ослабление излучения вследствие молекулярного рассеяния и поглощения газовыми компонентами (02, С02, Оз, Н20) в первых семи каналах МОБЕЗ, в которые не попадают основные линии поглощения этих газов, может быть определено с использованием предварительно созданных справочных таблиц и метеорологических данных МСЕР,

Ключевым элементом процедуры коррекции данных \10DIS в рассматриваемом приближении является определение аэрозольной оптической толщины. Это обусловлено не только большим вкладом аэрозоля в показания каналов \10DIS в видимом и ближнем инфракрасном диапазонах, но и отсутствием других оперативных данных по АОТ, которые могли бы быть использованы при коррекции. В реализованном в настоящее время операционном коде коррекции [7] аэрозольная оптическая толщина восстанавливается для каждого пикселя с разрешением 1 км в надире по полученным в процессе съемки интен-сивностям .\101)1Я. На этом этапе коррекции также используются справочные таблицы, построенные для различных типов аэрозолей и геометрий наблюдения с использованием кода 6Э [11].

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

2. Двунаправленный коэффициент спектральной яркости и альбедо подстилающей поверхности по данным MODIS

Период повторения орбит спутника Terra составляет 16 дней. В каждый 16-дневный цикл некоторая зона поверхности Земли наблюдается MODIS при различных углах В силу этого откорректированные на атмосферные эффекты данные MODIS в каждом цикле наблюдений могут быть использованы для восстановления усредненных ДКСЯ и последующего определения альбедо ПП по алгоритму [12].

Процедура восстановления базируется на утверждении, что линейная комбинация индикатрисы оптико-геометрической модели Кгео и расчетные данные в рамках теории переноса Ктп, дополняющие друг друга в описании прохождения квантов солнечного излучения через однородные и неоднородные растительные покровы, позволяют удовлетворительно моделировать ДКСЯ подстилающей поверхности (см., например, [8,12-15] и приведенные в них ссылки).

В операционном алгоритме для МОБК [12] двунаправленный коэффициент спектральной яркости подстилающей поверхности представляется в виде

n=3

Rs^, в, $,*>) = £ /к(Л)Кв, $, ф) = /1(Л)^0 + /2(Л)КП(в, $, ф) + /з(Л)Кгео(в, $, ф), (4)

k=i

где Ки30 = 1 — индикатриса изотропного рассеивания;

К

(п/2 — £) cos £ + sin £ п cos в + cos $ 4'

1

Ктео = 0(6*, </?) — sec в — sec # + - (1 + cos^) sec в sec

О = max ( О, — {t — sin t cos t) (sec в + sec •&) п

(5)

(6)

cos t = max

. o4AD2 + (tgfltg#sin фу 0, mm I 2-—---, 1

sec в + sec $

D = max (o, \Ag2 в + tg2 tf - 2 tg в tg & cos ф

cos £ = cos в cos $ + sin в sin $ cos ф.

Коэффициенты разложения /k(Л) в (4) являются параметрами, которые необходимо определить по результатам спутниковых наблюдений. Если рДЛ, в1, $i, ф), l = 1,..., m, — коэффициенты спектральной яркости подстилающей поверхности, полученные по резуль-m

квадратичной ошибки

2 _ 1 ^ (ps(K0u$u 4>l) - Rs( Л, eh <&h w))2

1=1

WA,l

(7)

нетрудно определить /(Л) В выражении (7) й = т — 3 — число степеней свободы; — вес, присвоенный измерению р8(Л,в1,$1 , ф).

Получаемая после минимизации система линейных алгебраических уравнений относительно коэффициентов /

n / m

tí ) tí Kk,i = Kk(в1, $i, ф1), k =1, 2, 3, приводит к решению

(8)

/к = Е

-i

^WA,l J JWKj

l=1

где I У2 ~'Ki,i

l=1

w\,i

— обратная матрица системы (8),

После определения Д и восстановления двунаправленного коэффициента спектральной яркости подстилающей поверхности по (4)-(6) можно вычислить плоское й™ (0) и сферическое аЛ альбедо:

2п

п/2

n/2

апф (9) = I dp I d$ps(A, 9, p) sin & cos af = 2 / dвапАл (9) sin в cos 9

Реальное спектральное альбедо йл (0) подстилающей поверхности представляет собой суперпозицию [9,12|

йл(0) = [1 - ^(0, та)] аЛ(0) + ^(0, тл)а?,

где 5(0, та) — доля диффузионного излучения в падающем на подстилающую поверхность потоке, которая определяется в основном оптической толщиной атмосферы та.

3. Технологии вычислений

Основой программных комплексов, использованных в данной работе при восстановлении характеристик ПП по данным MODIS, являются базовые алгоритмы [16|. Они получены Центром космического мониторинга Алтайского государственного университета из лаборатории DRL (Direct Readout Lab, GSFC/XASA) под лицензией, допускающей их изучение, использование, изменение, а также передачу измененных программ другим пользователям, Базовые алгоритмы сгруппированы в PGE (Product Generation Executive), каждый из которых содержит исходные тексты программ, справочные таблицы, обзор кода, описание требуемых входных и получаемых выходных файлов, а также краткое руководство по сборке и запуску, PGE, являющиеся копией операционных комплексов для MODIS, предназначены для работы в составе обрабатывающего дерева, вследствие чего запуск последующих PGE возможен лишь при наличии результатов работы предыдущих. Сборка и

Рис. 1. Последовательность обработки данных MODIS при проведении атмосферной коррекции.

Рис. 2. Последовательность обработки откорректированных на атмосферные эффекты данных MODIS при восстановлении ДКСЯ и альбедо.

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

Предоставляемые PGE комплектуются системой сборки па основе программы Make. Однако ввиду наличия у Make определенных недостатков для части из используемых базовых алгоритмов была разработана новая система сборки на основе GXU Autotools. В настоящее время используются как PGE, собранные штатной системой сборки, так и PGE, собранные с использованием повой системы.

Последовательность обработки данных MODIS при проведении атмосферной коррекции, восстановлении ДКСЯ и альбедо представлена на рис. 1, 2. Указанные на этих рисунках продукты вычислялись последовательно с использованием вычислительных комплексов IMAPP [17| и PGE 03, 11-13, 22, 23 четвертой версии [16|,

4. Результаты

Последовательное применение семи вычислительных комплексов позволяет восстанавливать коэффициенты спектральной яркости подстилающей поверхности для каждого пикселя сцеп за 16-дневный период и затем усредненные за этот период ДКСЯ, плоское и сферическое альбедо ПП с разрешением 1 км в первых семи каналах МОБК, а также в видимом (Л = 0.3... 0.7 мкм), ближнем инфракрасном (Л = 0.7... 5.0 мкм) и во всем коротковолновом (Л = 0.3... 5.0 мкм) диапазонах длин волн. На рис. 3 и 4 и в таблице в качестве примера показаны для 16-дневного периода (11 июпя-26 июня 2006 I'.) КСЯ зоны 54-57° с.ш,, 76-84° в.д., приведенные к надирной геометрии наблюдения при

Рис. 3. Коэффициенты спектральной яркости, приведенные к надир ной геометрии наблюдения. Каналы 1, 4, 3 МОБТБ.

Рис. 4. Плоское альбедо для 16-дневного периода (11 июня 26 июня 2006 г.) для коротковолнового диапазона 0.3.. . о.Омкм. Характерная деталь справа Обское водохранилище.

а б в

Рис. о. Доля пикселей сцены с различными значениями ХБУ1 для 13.10.04 (а), 15.06.06 (б) и 18.06.06 (в): 1 ХБУ1 без коррекции на атмосферные эффекты, 2 ХБУ1 с учетом атмосферной коррекции, 3 доля пикселей сцены, имеющая данное изменение ХБУ1 из-за учета атмосферной коррекции.

Данные по плоскому коротковолновому альбедо, полученные в настоящей работе для периода 11 июня-26 июня 2006 г. и операционные данные NASA [18] для периода 10 июня-25 июня 2006 г. для ПП в районе Новосибирска. Коды и типы ПП приведены согласно данным MOD12Q1 [18]

Код ПП Альбедо Количество пикселей

Тип подстилающей поверхности Настоящая работа Данные NASA

0 Водные объекты 0.041 0.043 1381

1 Вечнозеленые хвойные леса 0.094 0.091 194

3 Листопадные хвойные леса 0.125 0.125 143

5 Смешанные леса 0.132 0.133 5308

7 Кустарник 0.110 0.109 266

8 Лесные саванны 0.137 0.137 3004

9 Саванны 0.124 0.123 29

10 Луга 0.119 0.118 88

11 Постоянно влажные земли 0.090 0.087 5

12 Посевы 0.136 0.136 25115

13 Городские и застроенные земли 0.125 0.120 413

14 Посевы и природная растительность 0.146 0.145 910

16 Пустоши 0.106 0.105 8

среднем для каждого пикселя угле освещения, а также плоское альбедо для коротковолнового диапазона 0.3... 5.0 мкм. На рис. 5 проиллюстрировано влияние атмосферной коррекции на вегетационный индекс NDVI,

Заключение

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

Удовлетворение требований по ежемесячному режиму обновления данных по коротковолновому альбедо и его абсолютной точности [1] стимулировало развитие не только приборной базы исследования подстилающей поверхности из космоса, но и алгоритмов восстановления ее характеристик [12]. С выводом на орбиту на платформах Terra и Aqua 36-канального спектрорадиометра MODIS у исследователей появилась возможность ежедневного мониторинга подстилающей поверхности. Дополнение режима прямого вещания (Direct Broadcast) данных MODIS режимом свободного получения данных (Direct Readout), включая их получение с использованием базовых алгоритмов, позволило Центру космического мониторинга Алтайского государственного университета в 2002 г. начать мониторинговые наблюдения ряда критических параметров системы атмосфера — подстилающая поверхность Земли Сибирского региона.

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

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

Выполненные нами исследования подтверждают вывод [7] о том, что аэрозольный эффект в NDVI не устраняется полностью даже при использовании самой чистой сцены в анализируемый период (рис. 5). Более того, в Сибирском регионе в лесопожарный период этот эффект может существенно исказить картину состояния растительного покрова. Причиной являются как дымовой аэрозоль региона, так и аэрозоль, поступающий из центральной и восточной частей России. В силу сложности разделения городского и дымового аэрозолей при восстановлении АОТ атмосферы над сушей [2] для осуществления более точной атмосферной коррекции требуется подключение к вычислительным комплексам обработки спутниковых данных блока переноса дымового аэрозоля.

Из совместного анализа данных по типу подстилающей поверхности, представленного в продукте MOD12Q1 [18], и альбедо следует, что наблюдаемые сегодня альбедо в ряде зон региона не совсем точно соответствуют типу подстилающей поверхности.

Авторы выражают благодарность сотрудникам лаборатории DEL (Direct Readout Lab, GSFC/NASA) за передачу базовых алгоритмов PGE. Продукт MOD12Q1 и данные по альбедо (MOD43, ячейка (22,3)) для периода 10 июня-25 июня 2006 г., использованные авторами при сопоставлениях, представлены Центром активных распределенных архивов (DAAC) [18].

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

[1] Sellers P.J., Meeson B.W., Hall F.G. et al. Remote sensing of the land surface for studies of global change: models — algorithms — experiments // Remote Sens. Environ. 1995. Vol. 51. P. 3-26.

[2] Kaufman Y.J., Tanre D., Remer L.A. et al. Operational remote sensing of tropospheric aerosol over-land from EOS moderate resolution imaging spectroradiometer //J. Geophys. Res. 1997. Vol. 102, N D14. P. 17051-17 067.

[3] Salomonson V.V., Barnes W.L., Maymon P.W. et al. MODIS: Advanced facility instrument for studies of the Earth as a system // IEEE Trans. Geosci. Remote Sens. 1989. Vol. 27. P. 145-153.

[4] Ackerman S.A., Strabala K.I., Menzel W.P. et al. Discriminating clear sky from cloud with MODIS // J. Geophys. Res. 1998. Vol. 103, N D24. P. 32 141-32 157.

[5] Gao B.C., Kaufman Y.J. Remote sensing of water vapor in the near IR from EOS/MODIS // IEEE Trans. Geosci. Remote Sens. 1992. Vol. 30. P. 871-884.

[6] Сни D.A., Kaufman Y.J., Zibordi G. et al. Global monitoring of air pollution over land from Earth Observing System — Terra Moderate Resolution Imaging Spectroradiometer (MODIS) // J. Geophys. Res. 2003. Vol. 108, N D21. doi:10.1029/2002JD003179.

[7] Vermote E.F., Saleous N.S.El., Justice C.O. Atmospheric correction of MODIS data in the visible to middle infrared: first results // Remote Sens. Environ. 2002. Vol. 83. P. 97-111.

[8] Schaaf C.B., Gao F., Strahler A.H. et al. First operational В RDF, albedo nadir reflectance products from MODIS // Remote Sens. Environ. 2002. Vol. 83. P. 135-148.

[9] Минин И.Н. Теория переноса излучения в атмосферах планет. М.: Наука, 1988. 264 с.

[10] Vermote Е.F., Vermenlen A. Atmospheric correction algorithm: spectral reflectances (MOD09). http: /modarch. gsf с. nasa. gov/M0DIS/ATBD/atbd_mod08. pdf.

[11] Vermote E.F., Tanre D., Deuze J.L. et al. Second simulation of the satellite signal in the solar spectrum: an overview // IEEE Trans. Geosci. Remote Sens. 1997. Vol. 35. P. 675-686.

[12] Lucht W., Schaaf C.B., Strahler A.H. An algorithm for the retrieval of albedo from space using semiempirical BRDF models // IEEE Trans. Geosci. Remote Sens. 2000. Vol. 38, N 2. P. 977-998.

[13] Roujean J.L., Leron M., Deschamps P.Y. A bidirectional reflectance model of the Earth's surface for the correction of remote sensing data //J. Geophys. Res. 1992. Vol. 97, N D18. P. 20 455-20 468.

[14] Wanner W., Li X., Strahler A.H. On the derivation of kernels for kernel-driven models of bidirectional reflectance // J. Geophys. Res. 1995. Vol. 100, N D10. P. 21077-21089.

[15] Математическое моделирование переноса радиации в растительных средах / Ю. Росс, Ю. Князихин, А. Кууск и др. СПб.: Гидрометеоиздат, 1992. 200 с.

[16] http: //directreadout. gsf с. nasa. gov/dowiiload_techiiology/iiist_algorythms. cf m.

[17] International MODIS/AIRS Processing Package (IMAPP). http://cimss.ssec.wise.edu/~gumley/IMAPP/IMAPP.html.

[18] EOS Data Gateway, http://edcimswww.cr.usgs.gov/pub/imswelcome/.

Поступила в редакцию 19 октября 2006 г.