Контроль загрязнения природной среды аэрокосмическим методом Текст научной статьи по специальности «Нанотехнологии»

■ новых образцов целевой аппаратуры КА ДЗЗ, включая бортовой радиолокационный комплекс СВЧ-диапазона с синтезированной апертурой и бортовой гиперспектральный оптико-электронный комплекс;

■ экспериментальных образцов научного оборудования, микроэлектронных элементов, конструкционных материалов и покрытий, адаптированных для использования на борту микроспутника ДЗЗ;

■ новых информационных технологий приема, обработки, хранения и доведения до широкого круга потребителей данных ДЗЗ от КА с перспективной аппаратурой наблюдения, орбитальной группировки КА на основе применения современных достижений вычислительной техники и телекоммуникационных средств;

■ опытных аппаратно-программных комплексов и технических средств для обеспечения подготовки, переподготовки и повышения квалификации специалистов в области космического мониторинга Земли.

Рис. 5.

СО Технология

о (N и аппаратно-

>s программный

го комплекс учета

О состояния лесного

LO Z фонда и ресурсной

S оценки поврежден-

и < ных насаждений

ей О в результате

X X воздействия на них

S S неблагоприятных

< * > природно-клима-

< X тических факторов

10

Следует ожидать, что реализация новой программы «Мо-ниторинг-СГ» станет логичным продолжением долгосрочного сотрудничества двух государств по интеграции их наземных и орбитальных космических средств, а также интеллектуальных и производственных ресурсов в освоении космического пространства в мирных целях, расширении технических возможностей национальных систем ДЗЗ в интересах решения актуальных проблем экономики Беларуси и России. ■

Контроль загрязнения природной среды

аэрокосмическим методом

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

Важнейшая задача спутникового мониторинга Земли - контроль загрязнения атмосферы, в частности антропогенного, вызванного достаточно мощными выбросами локальных источников. Над решением этой проблемы работают во многих исследовательских центрах России, США (NASA), Европы (ESA), Японии и других стран.

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

Необходимо обеспечить приемлемую величину отношения сигнала к шуму в измерительных космических приборах, и это заставляет идти на компромисс при выборе пространственного разрешения и ширины спектрального канала прибора. Многозональные съемочные системы (МСС) высокого пространственного разрешения (~10-15 м) имеют сравнительно широкие спектральные каналы (-70-100 нм) (ASTER на американ-

ской платформе TERRA, Белорусский космический аппарат (БКА), российский «Канопус-В» и др.). На многоспектральных приборах (МСП), формирующих изображения с большим числом узких (~10 нм) спектральных каналов (MODIS, MERIS и др.), можно проводить атмосферную коррекцию, но со сравнительно грубым пространственным разрешением (-300-1000 м). В рамках программ Союзного государства «Космос-СГ» и «Космос-НТ» описана совместная обработка данных МСП и МСС с целью восстановления оптических характеристик аэрозольной атмосферы и поверхности Земли с высоким пространственным разрешением. Такой подход предполагается в дальнейшем использовать при анализе информации БКА.

Модель атмосферы и подстилающей поверхности

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

всех этих компонент. Система представляется состоящей из подстилающей поверхности, нижнего слоя «1» атмосферы, который считается однородным, и верхнего слоя «2» атмосферы - неоднородного.

Спектральные модели поверхностей берутся такие же, как в коде БЛЕЙ [1]. В частности, альбедо (спектральный коэффициент отражения) суши г (Л) представляется в виде линейной комбинации двух базовых функций, описывающих соответственно спектры растительности г (Л) и почвы г и(Л):

г(А)= сг^1-с)^(Л). (1)

Для альбедо водной поверхности в качестве базовых функций берутся спектры чистой океанской воды гЫтг(А) и прибрежных вод г^а,(А): С€

Г5(Л) СГс1еаг(Л) (1 C)Гcoastal('Л). (2)

Здесь с - доля поверхности в пикселе, занимаемая растительностью или чистой водой (рис. 1). Таким образом, альбедо подстилающей поверхности характеризуется одним неизвестным параметром с, который должен быть определен в итерационном процессе.

Нижний слой «1» составляет тропосфера до некоторой высоты Н (порядка 2-3 км), величина которой принимается равной для всех пикселей обрабатываемого изображения. Считается, что в слое происходит аэрозольное и молекулярное рассеяние, аэрозольное и газовое поглощение. В рассматриваемых далее алгоритмах используется модель аэрозоля типа Сопйпвп1а1 [2]. Количество аэрозоля в слое может заметно варьироваться во времени и пространстве.

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

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

Расчет СКЯ,

регистрируемого многоспектральным прибором на космическом аппарате

Для решения обратной задачи (восстановления оптических характеристик атмосферы и земной поверхности по спектральному сигналу на летательном аппарате) необходимо определить СКЯ в верхних слоях атмосферы. О распространении излучения в однородном слое «1» можно узнать по численному методу удвоения [3] или приближенному полуаналитическому методу [4] теории переноса. Для расчета сигнала от слоя «2» применяется разработанный авторами ранее быстрый и достаточно точный код RAY [5, 6]. При этом можно использовать различные модели молекулярно-газовой и аэрозольной атмосферы - как детально структурированную по высоте, включающую до 52 подслоев, так и упрощенную стандартную модель Всемирной метеорологической организации [2], состоящую из четырех подслоев.

Для учета поглощения в узких полосах поглощения атмосферных газов используется модификация так называемого метода k-distribution [7]: для заданного спектрального интервала предварительно определяется база данных некоторых средних коэффициентов пропускания, которые затем используются, чтобы рассчитать перенос излучения с учетом многократного рассеяния. Результирующий сигнал на спутнике находится путем учета многократных переотражений света между указанными слоями атмосферы и подстилающей поверхности (почвой, травой, песком, растительностью, снегом, водой и др.).

Возможности применения моделей атмосферы иллюстрирует рис. 2. На нем приведены результаты расчетов с помощью кода RAY

04 0*5 О S О V, 0.6 065 0.7

О J он 0.5 0.55 ' Ь ОI: 'I ■

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

Подготовка входных многоспектральных данных и дискриминация облачных пикселей

В рассматриваемых нами алгоритмах восстановления оптических параметров атмосферы, характеризующих ее загрязнение [8, 4], используются данные спутникового инструмента MERIS в девяти спектральных каналах (на длинах волн 412,5, 442,5, 490, 510, 560, 620, 665, 865 и 885 нм). При этом осуществляется подготовка входных данных:

■ определение коэффициентов яркости Rtoa(A) на верхней границе атмосферы в указанных спектральных каналах для всех пикселей;

Рис. 1.

Базовые спектры альбедо суши (а) и воды (б):

а) почва (1), растительность (2),

б) прибрежная вода (1),

чистая океанская вода (2)

Рис. 2.

СПЭЯ на верхней границе атмосферы, рассчитанная с помощью кода RAY для детальной (52 подслоя) и упрощенной (4 подслоя) моделей атмосферы для случая, когда подстилающей поверхностью является растительность

о х

X

5

X

о

X X

5

X

Рис. 3. Результаты восстановления АОТ на длине волны 550 нм алгоритмом ART (слева) и другими известными алгоритмами [11] (справа) в сравнении с данными AERONET

Рис. 4.

Относительная спектральная чувствительность сенсора в каналах МСС БКА

■ дискриминация (отбраковка) облачных пикселей;

■ разделение пикселей по типу «суша» или «вода»;

■ выделение блока пикселей для осреднения.

Восстановление оптических характеристик атмосферы по многоспектральным данным

Оптические параметры атмосферы и подстилающей поверхности определяются из сравнения СКЯ, измеренных спутниковым оптическим инструментом на верхней границе атмосферы в N спектральных каналах, с рассчитанными. В большинстве современных алгоритмов используется метод lookup tables (LUT) [9]: для различных моделей атмосферы и систем регистрации рассчитываются СКЯ на верхней границе атмосферы. Достоинство метода - значительное сокращение объема вычислений. Недостаток его в том, что необходимо использовать огромный массив насчитанных данных и практически полностью заменять их при уточнении каких-либо свойств модели, в частности набора используемых длин волн, а также

в том, что невозможно учитывать многопараметрические модели ввиду резкого увеличения объема информации.

В рассмотренных нами алгоритмах ART [8] и FAR [4] не используется метод LUT, а для обеспечения быстродействия применяется быстрый код RAY [5]. Кроме того, в алгоритме FAR [4] для расчета переноса излучения в слое «1» используются приближенные полуаналитические решения теории переноса. Вследствие отказа от метода LUT выбор спектральных каналов в этих алгоритмах достаточно гибок - их число и положение может изменяться в зависимости от характеристик оптического инструмента.

Для каждого центрального пикселя блока с использованием полученных средних значений R(A) в семи спектральных каналах (на длинах волн 412,5, 442,5, 490, 510, 560, 620 и 665 нм) восстанавливаются параметр подстилающей поверхности с (см. формулы (1), (2)), оптическая толщина т1(412) слоя «1» атмосферы при А = 412 нм и параметр Ангстрема а для этого же слоя. При этом предполагается, что спектральная зависимость описывается соотношением А

,5

Оптическая толщина определяется функцией

т1(А)=т1(412)ЫЪГ

(3)

г1(А) =/г(А, h) dh,

(4)

где е(А, Ь) - показатель ослабления элементарного объема атмосферы на высоте Ь, зависящий от концентрации газов и частиц аэрозоля. Оптическая толщина г1(А) связана с коэффициентом пропускания нерассеянного

солнечного излучения Т(А) соотношением Т(А) = ехр(-т1). Таким образом, т1(А) определяется общим содержанием компонент воздуха в толще Н. Параметр Ангстрема а характеризует спектральную зависимость среднего по толщине атмосферы аэрозольного показателя ослабления. Для частиц, много меньших длины волны света, он равен четырем, для больших - нулю. По величине а можно судить о происхождении и типе загрязняющих атмосферу частиц. Оптическая толщина всей атмосферы т(А) = т1(А) + т2(А), где т2(А) - оптическая толщина слоя «2». Искомые параметры т1(412), а и с, которые не зависят от длины волны, определяются по методу наименьших квадратов в итерационном процессе. Точность алгоритма [8] иллюстрирует рис. 3, где для одного из районов Германии дано сопоставление восстановленной аэрозольной оптической толщины (АОТ) с данными мировой сети AERONET [10] по измерению ослабления прямого солнечного излучения. Здесь же для сравнения приведены результаты восстановления АОТ другими известными алгоритмами [11].

Восстановление альбедо подстилающей поверхности по многоспектральным данным

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

RT(U(U,Ui),V) -

=RJß,ß0,<Ph

T(ß)T(ß0)rs

1- r* r

atm s

(5)

где Ra,m(V,Ц0,ф) - коэффициент яркости на верхней границе атмосферы при альбедо поверхности, равном нулю, Т(ц) - коэффициент пропускания атмосферы, г*т -альбедо атмосферы при диффузном освещении ее снизу. Здесь ц0 и ц - косинусы углов падения солнечного излучения и наблюдения, ф - разность азимутальных углов наблюдения и падения

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

Восстановление альбедо подстилающей поверхности с высоким пространственным разрешением по данным МСП

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

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

Применительно к МСС высокого пространственного разрешения эффективный коэффициент яркости на верхней границе атмосферы можно представить в виде:

кмТ0А,р(^^) =

= Яа1т(ЛЛ,Ц,Цо,Ф) + + Кг (ДЛ,и,и) + К (А\,ц,и,ф), (6) где ка1т (АЛ,ц,цыф) - вклад излучения, рассеянного непосредственно атмосферой, без взаимодействия с подстилающей поверхностью, К (АЛ,ц,ц.ф) -

вклад регистрируемого пикселя, который зависит от эффективного (по спектральному интервалу ЛА) альбедо поверхности fj (ЛА), Rтг(ЛА,ц,ц), - вклад подстилающей поверхности, окружающей регистрируемый пиксель. Соотношение (6) позволяет определить альбедо пикселя ^^(ЛА). При этом данные об атмосфере и окружающей подстилающей поверхности находятся путем интерполяции по спектру результатов, восстановленных в семи каналах МСП.

Оценка точности определения альбедо подстилающей поверхности

Точность атмосферной коррекции путем совместной обработки данных МСП и МСС существенно зависит от верности спутниковых измерений. Формула (6) позволяет оценить влияние погрешности калибровки аппаратуры и правильности восстановления оптических характеристик атмосферы на точность определения альбедо подстилающей поверхности.

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

г 1 aer '

что если погрешность калибровки МСС < 4%, то ошибка восстановления альбедо поверхности может быть порядка 8% и менее (при альбедо > 0,1) и около 15% (при альбедо < 0,1). Заметим, что при малых альбедо поверхности основную погрешность вносит ошибка калибровки МСС, а при больших - неточность восстановления аэрозольной оптической толщины атмосферы по многоспектральным спутниковым данным. При более высоком уровне погрешности калибровки сенсора МСС (5% и более) она практически всегда определяет результирующую погрешность Лгр / rp восстановления альбедо пикселя. При этом, хотя относительная

погрешность Arp / rp в некоторых случаях (в коротковолновой части спектра) может равняться 40-50%, абсолютная ошибка Arp составляет порядка 1% и менее, что объясняется малой величиной альбедо поверхности в этой области спектра. Полученные результаты достаточно хорошо согласуются с оценкой погрешностей для спутникового инструмента ASTER [12], характеристики которого сходны со свойствами БКА.

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

На рис. 5 представлены результаты обработки данных космического инструмента MERIS на спутниковой платформе ENVISAT. Слева показана чистая, справа - сильно замутненная дымами лесных пожаров атмосфера. В верхней части рисунка даны композитные RGB-изобра-жения по данным инструмента MERIS. На них видны облака, а также подстилающая поверхность, вуалированная более (справа) или менее (слева) загрязненной атмосферой. Количественную оценку

Рис. 5. Композитные Ю-изобра-жения по данным инструмента МЕМБ на спутниковой платформе ЕМУШ 11 октября 2005 г. (вверху слева) и 1 мая 2006 г. (вверху справа) и восстановленные аэрозольные оптические толщины на длине волны 412,5 нм (внизу)

Рис. 6.

Восстановление СКЯ дерново-глеевой среднесуг-линистой почвы.

1 - истинные СКЯ почв;

2 и 5 - СКЯ на верхней границе атмосферы и на высоте h = 6 км;

3 и 4 - результат атмосферной коррекции с региональной

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

о

X X

5

X

о

X X

5

X

Рис. 7.

Альбедо подстилающей поверхности участка Европы на длинах волн 412,5 нм (слева) и 450 нм (справа), восстановленные по данным инструмента MERIS

Рис. 8.

Синтезированное RGB-изображение района Солигор-ского водохранилища до (слева) и после (справа) «срезания» атмосферной «дымки»

f . Ct '"— : .- 1 JSP ■ , уЛ. J 1 t , ' .. f V

гтад Tj-t^ J* '¿feh. ttWi^^Klrafe,

загрязненности атмосферы дают ее восстановленные оптические толщины на длине волны 412,5 нм (в нижней части рисунка). Здесь же приведена шкала оптических толщин. Белые участки означают, что здесь значения АОТ выходят за пределы шкалы, то есть ттг > 1.

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

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

Следует заметить, что описанная процедура необходима не

..... '

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

Аркадий Иванов,

главный научный сотрудник Института физики НАН Беларуси, член-корреспондент

Иосиф Кацев,

ведущий научный сотрудник Института физики НАН Беларуси, кандидат физико-математических наук

Элеонора Зеге,

ведущий научный сотрудник Института физики НАН Беларуси, кандидат физико-математических наук

Александр Прихач,

старший научный сотрудник Института физики НАН Беларуси

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

Однако для тематической обработки спутниковой информации, в частности, для контроля

Литература

1. Von Hoyningen-Huene W., Freitag M., Burrows J. B. Retrieval of aerosol optical thickness over land surfaces from top-of-atmosphere radiance. J. Geoph. Res., 2003, №108. Р. 4260-4272.

2. WCP-112. A Preliminary Cloudless Standard Atmosphere for Radiation Computation. World Meteorological Organization, Geneva, 1986.

3. Lenoble J. Radiative Transfer in Scattering and Absorbing Atmospheres: Standard Computational Procedures. A. Deepak Publishing, Hampton, Virginia, 1985.

4. Katsev I. L., Prikhach A.S., Zege E.P., Ivanov A.P., Grudo J.O., Kokha novsky A.A. Speeding up the aerosol optical thickness retrieval using analytical solutions of radiative transfer theory. Atmos. Meas. Tech., 2010, №3. Р. 1403-1422.

5. Tynes H., Kattawar G.W., Zege E.P., Katsev I.L., Prikhach A.S., Chaikovskaya L.I. Monte Carlo and multicomponent approximation methods for vector radiative transfer by use of effective Mueller matrix calculations. Appl. Opt., 2001, №40. Р. 400-412.

6. Kokhanovsky A.A., Cornet C. , Duan M., Emde C., Katsev I.L., Labonnote L.C., Min Q., Nakajima T., Ota Y., Prikhach A.S., Rozanov V.V., Yokota T., and Zege E. P. Benchmark results in vector radiative transfer. Journal of Quantitative Spectroscopy & Radiative Transfer, 2010, №111. Р. 1931-1946.

7. Кацев И.Л., Прихач А.С., Зеге Э.П. Методика быстрого расчета переноса излучения в атмосфере в линиях поглощения газов. Материалы Третьего Белорусского космического конгресса, конгресса, ОИПИ НАН Беларуси. - Мн., 2007. С. 209-213.

8. Katsev I.L., Prikhach A.S., Zege E.P., Ivanov A.P., Kokhanovsky A.A. Iterative procedure for retrieval of spectral aerosol optical thickness and surface reflectance from satellite data using fast radiative transfer code and its application to MERIS measurements. Satellite Aerosol Remote Sensing over Land, ed. by A.A. Kokhanovsky and G. de Leeuw, Berlin: Springer-Praxis, 2009. P. 101-134.

9. Kokhanovsky A.A., Deuze J. L., Diner D. J., Dubovik O., Ducos F., Emde C., Garay M.J., Grainger R.G., Heckel A., Herman M., Katsev I.L., Keller J., Levy R., North P.R.J., Prikhach A.S., Rozanov V.V., Sayer A.M., Ota Y., Tanr'e D., Thomas G.E., and Zege E.P. The inter-comparison of major satellite aerosol retrieval algorithms using simulated intensity and polarization characteristics of reflected light. Atmos. Meas. Tech., 2010, №3. Р. 909-932.

10. Holben B.N., Eck T.I., Slutsker I., Tanr'e D., Buis J.P., Setxer A., Vemote E., Reagan J.A., Kaufman Y.J., Nakajima T., Lavenu F., Jankowiak I., Smirnov A. AERONET-A federated instrument network and data archive for aerosol characterization. Remote Sens. Environ., 1998, №66. Р. 1-16.

11. Kokhanovsky A.A. et al. Aerosol remote sensing over land: A comparison of satellite retrievals using different algorithms and instruments. Atmos. Res., 2007, №85. Р.372-294.

12. Thome J. et al. Vicarious Calibration of ASTER via the Reflectance-Based Approach. IEEE Transactions on Geoscience and Remote Sensing, 2008, №46. Р. 3285-3295.