В.Н.Лагуткин, А.П.Лукьянов, Е.Н.Подкорытов, В.Г.Репин

 

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

 

 

Аннотация

 

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

 

 

1. Введение. Основные явления, учитываемые в модели.

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

·        спектральная зависимость собственного теплового излучения природных объектов (подстилающей поверхности, атмосферы, облаков, аэрозолей);

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

·        спектральные зависимости ослабления излучения на трассах Солнце-природный объект -сенсор;

·        высокая вероятность наблюдения облачных систем, содержащих  облака различных ярусов; 

·        эффект взаимного затенения и экранирования одних объектов другими;

·        спектральные и пространственные характеристики оптической толщины облаков, главным образом облаков верхнего яруса;

·        динамические изменения изображений, вызванные движением облаков, наблюдателя и колебаниями оптической оси;

·        искажения изображений природных объектов при прохождении излучения в атмосфере;

·        преобразование изображения в оптической системе и фотоприемном устройстве наблюдателя;

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

 

2. Общее описание модели.

Для моделирования последовательности изображений полей яркости фона используются следующие входные данные:

·        координатно-временные данные, содержащие информацию о положении и движении наблюдателя, точке наведения, дате и времени проведения наблюдений;

·        синоптические данные, включающие информацию о структуре облачности (циклоны, фронты, волны), классах облаков (перистые, слоистые, кучевые) и их основных параметрах (бальность, высота,  спектр размеров, скорость ветра) в зоне наблюдения;

·        данные о рельефе земной поверхности в районе наблюдения и ее отражательных характеристиках;

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

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

Расчеты включают следующие основные этапы:

  • Расчет альбедо  и индикатрис однократного рассеяния  облачных частиц (капель воды и кристаллов льда) в заданном спектральном поддиапазоне. Для расчетов используются асимптотические методы решения задачи дифракции [1-4]. В качестве базы данных используются спектральные зависимости комплексных показателей преломления воды и льда и известные распределения облачных частиц по размерам и формам [5,6];

·        Расчет индикатрис излучения в облаках  с учетом явления многократного рассеяния излучения. Расчет основан на численном решении уравнения переноса излучения

,

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

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

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

 

3. Проверка адекватности модели по имеющимся экспериментальным данным.

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

Примером условий, в которых теоретическое решение известно, является случай полубесконечной среды, имеющей плоскую границу и состоящей из частиц, плотность которых постоянна, а индикатриса однократного рассеяния имеет вид  ( - угол рассеяния). Количественное представление о точности разработанной статистической модели дают графики, представленные на рис. 1, на которых изображены функции КСЯ в трех азимутальных полуплоскостях, ориентированных под углами 0°, 90° и 180° относительно проекции направления падающего излучения луча на плоскость (соответственно кривые 1, 2, 3). Зенитный угол падающего излучения составляет 72°. При таком угле падения азимутальное распределение КСЯ слоя становится существенно неравномерным и более информативным для сопоставления результатов. Результаты теоретических расчетов имеют вид плавных кривых, а результаты моделирования имеют вид ступеней. Верхняя кривая КСЯ соответствует азимутальному направлению с минимальным углом рассеяния (вперед), а нижняя - направлению с максимальным углом рассеяния (назад). Из приведенных графиков наряду с констатацией факта несомненного соответствия расчетов и теории, можно отметить ту особенность, что относительная погрешность метода Монте-Карло увеличивается при малых величинах индикатрисы рассеяния. Эту особенность наряду со сделанными выше замечаниями о размерах дискретов представления выходных данных нужно учитывать при выборе числа статистических испытаний метода Монте-Карло, если нас интересует поведение индикатрисы рассеяния вне зоны ее основного максимума.

Рис. 1

Перейдем теперь к сопоставлению СПЯ, рассчитываемых по модели, с СПЯ, экспериментально полученными при наблюдениях из космоса (рис. 2). Экспериментальные данные, охватывающие диапазон длин волн от 2.1 до 3.3 мкм, относятся к наблюдениям водяных и ледяных облаков с помощью фурье-спектрометра с достаточно высоким разрешением (0.02 мкм) [9]. Наблюдения проводились в надир при различных зенитных углах Солнца. Расчеты проводились для тех же условий. Для корректности сопоставления дополнительно учитывались поглощение излучения в атмосфере и аэрозольное рассеяние над облаками.

На рис. 2а проводится пример сопоставления расчетных (кривая 1) и экспериментальных (кривая 2) СПЯ (в дБ к Вт/м2/ср/мкм ) ледяных облаков для зенитного угла Солнца ~60°. На рис. 2б аналогичное сопоставление проводится для водяных облаков при зенитных углах Солнца ~40°. В расчетах мы полагали высоту облаков равной 14 км для ледяных облаков и 3 км - для водяных. Графики свидетельствуют о достаточно хорошем согласии расчетов и экспериментов для рассмотренных условий.

                             а)                                                               б)

Рис. 2

При больших зенитных углах Солнца, превышающих 80°, расчет дает несколько более низкие величины СПЯ, чем эксперимент, причем расчетные величины, в отличие от экспериментальных, спадают быстрее, чем по ламбертовскому закону. Это может быть объяснено особенностью эксперимента, связанной с тем, что верхняя граница облаков не строго горизонтальна и реальные углы падения солнечного излучения не равны в точности зенитному углу Солнца. При этом эффект негоризонтальности возрастает по мере приближения направления падающих лучей к скользящему и результат эксперимента становится более чувствительным к выбору участка облачной поверхности, по которому проводятся измерения. Поскольку вопрос нацеливания оптической системы при экспериментальных измерениях в [9] не обсуждался, такое объяснение, в принципе возможно. Косвенное подтверждение такого объяснения дают данные других экспериментов [10], из которых следует, что КСЯ рассеянного излучения облаков имеет минимум при углах рассеяния около 90°, соответствующих обсуждаемым условиям, и увеличивается в 2-3 раза при увеличении углов рассеяния до 180°, т.е. рассеяние явно не является ламбертовским и экспериментальная зависимость КСЯ от угла рассеяния подобна расчетной.

 

 4. Принципы использования многоспектральных спутниковых изображений.

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

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

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

2)      параметры многоспектральной аппаратуры наблюдения:

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

·      пороговая яркость обнаружения на входном зрачке для каждого рабочего спектрального поддиапазона,

·      угловые или линейные (на поверхности Земли) размеры элемента разрешения (пикселя);

3)      параметры, характеризующие условия наблюдения:

·      время и дата наблюдения,

·      географические координаты спутника,

·      географические координаты условного центра выбранного фрагмента,

·      размеры фрагмента (в пикселях),

для синтеза изображений дополнительно задаются:

4)      географические координаты наблюдающего сенсора,

5)      параметры аппаратуры наблюдения сенсора:

·        параметры функции пропускания фотоприемного устройства или условные границы спектрального поддиапазона,

·        пороговая яркость обнаружения на входном зрачке,

·        угловые или линейные размеры элемента разрешения.

В состав выходной информации обработки входят следующие данные.

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

2)      карта и гистограмма оценок коэффициентов рассеяния (альбедо) облаков и подстилающей поверхности,

3)      карта и гистограмма оценок яркостной температуры верхней границы облаков и подстилающей поверхности,

4)      карта и гистограмма оценок коэффициентов излучения облаков и подстилающей поверхности,

5)      оценка среднего по кадру высотного профиля температуры атмосферы,

6)      синтезированное изображение выбранного фрагмента из заданной точки наблюдения в заданном спектральном поддиапазоне.

Обработка входной информации проходит в три этапа.

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

·        Расчет спектральных зависимостей индикатрис яркости излучения в облаках с учетом  многократного рассеяния излучения (методом Монте-Карло);

·        Расчет спектральных зависимостей коэффициентов рассеяния солнечного излучения земной поверхностью различных типов;

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

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

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

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

Для демонстрации работы комплекса алгоритмов обработки спутниковых изображений использовались многоспектральные изображения, полученные метеорологическим спутником NOAA-17 [11]. Основные результаты работы комплекса алгоритмов отображены на рис. 3-7.

В верхнем ряду рисунков показаны фрагменты изображений, полученные в различных спектральных диапазонах радиометра AVHRR/3 (Advanced Very High Resolution Radiometer, Version 3) с разрешением 1.1 км на поверхности Земли (в подспутниковой точке). Условные нижняя и верхняя границы поддиапазонов указаны сверху фрагментов. Снизу фрагментов  указаны минимальные (min) и максимальные (max) значения и пороговая яркость обнаружения на входном зрачке (N) для каждого рабочего спектрального поддиапазона (в Вт/м2 /ср). Размер фрагментов 128´128 пикселей. Размер фрагментов является входным параметром алгоритмов и может изменяться пользователем. Географические координаты подспутниковой точки (F-широта, L-долгота) и центра фрагментов, мировое время (UT) и дата (DATA) указаны с правой стороны.

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

·        карта оценок высот верхней границы облаков и природных типов подстилающей поверхности,

·        карта оценок коэффициентов рассеяния (альбедо) облаков и подстилающей поверхности,

·        карта оценок яркостной температуры верхней границы облаков и подстилающей поверхности,

·        карта оценок коэффициентов излучения облаков и подстилающей поверхности.

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

В среднем ряду показано также синтезированное изображение выбранного фрагмента из заданной точки наблюдения в заданном спектральном поддиапазоне. Условные границы спектрального поддиапазона и географические координаты наблюдающего сенсора (широта F, долгота L, высота H) указаны справа фрагмента. Снизу этой карты указаны статистические параметры: минимальное и максимальное значения, средняя величина (Mat) и СКО (Sig).

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

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

·      гистограмма оценок высот верхней границы облаков и природных типов подстилающей поверхности, причем значения H>0 соответствуют высоте верхней границы облаков (в км), а H£0 соответствуют определенным  природным типам подстилающей поверхности: H=0 - вода, H= -1- снег,      H= -2 - растительность, H= -3 - горные породы, H= -4 - почва, H= -5 – пустыня;

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

·      реальное (эталонное) изображение выбранного фрагмента в интересующем спектральном поддиапазоне (если таковое имеется).  

·      карта невязки между  синтезированным изображением и реальным (эталонным) изображением, если таковое имеется.

На рис.3, 4 представлены результаты тестирования комплекса алгоритмов. Тестирование проводилось при условии, что положение гипотетического сенсора совпадает с положением спутника NOAA-17. Для тестирования, результаты которого представлены на рис.3, в качестве эталонного использовался спектральный поддиапазон 0.73-1 мкм ближнего ИК диапазона и алгоритм синтезировал изображение для этого поддиапазона. Как видно из карты невязки между  синтезированным изображением и реальным эталонным изображением (крайний правый фрагмент в нижнем ряду), тестирование дает хорошие результаты: средняя величина невязки составляет »8%, а величина СКО »3% от максимального значения яркости на эталонном изображении. Идея тестирования, результаты которого представлены на рис.4, заключается в следующем. Из полного состава входных фрагментов изображений (5 фрагментов изображений в случае радиометра AVHRR/3 NOAA-17) исключается фрагмент изображения в одном из поддиапазонов и этот фрагмент используется в качестве эталонного, а алгоритм производит синтез изображения для выбранного поддиапазона по неполному составу входных данных. В качестве эталонного изображения для удобства сопоставления также использовался спектральный поддиапазон 0.73-1 мкм. И в этом случае тестирование дает хорошие результаты: средняя величина невязки составляет »12%, а величина СКО »4% от максимального значения яркости на эталонном изображении.

На рис. 5 представлены результаты комплекса алгоритмов в варианте синтеза изображения в полосе атмосферного поглощения 2.5-3 мкм при условии, что положение гипотетического сенсора совпадает с положением спутника NOAA-17. Эта демонстрация интересна тем, что, как и ожидалось, изображение в полосе поглощения атмосферы 2.5-3 мкм существенно отличается и по виду и по уровню величин яркости от изображения на рис. 2, что позволяет качественно судить об адекватности алгоритмов современным экспериментальным данным и теоретическим моделям.   

На рис. 6 представлены результаты комплекса работы алгоритмов в случае синтеза изображения в поддиапазоне 0.73-1 мкм при условии, что положение гипотетического сенсора отличается от положения спутника NOAA-17. Искажение фрагмента обусловлено изменением геометрии наблюдения. Реальное изображение для этого поддиапазона в качестве эталонного помещено в нижнем ряду. Этот рисунок показывает насколько значительно изменяется поле яркости при изменении положения сенсора и демонстрирует возможности предлагаемого программно-алгоритмического обеспечения.

Следующий рис.7 демонстрирует возможности применения предлагаемого программно-алгоритмического обеспечения в варианте, когда не только положение, но и спектральный поддиапазон гипотетического сенсора отличаются от положения и спектральных поддиапазонов спутника NOAA-17.

 

6. Выводы.

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

Модель позволяет получать:

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

·        последовательности многоспектральных изображений выбранного района из заданной точки наблюдения;  

·        последовательности стереоизображений фона выбранного района из различных точек наблюдения;  

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

 

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

1. Зуев В.Е. Распространение видимых и инфракрасных волн в атмосфере. М.: Сов. радио, 1970.

2. Борен К., Хафмен Д. Поглощение и рассеяние света малыми частицами. М.: Мир, 1986.

3. Хенл Х., Мауэ А., Вестпфаль К. Теория дифракции. М.: Мир, 1964.

4. Hansen J.E. and Pollack J.B. // J. Atmos. Sci., 1970, V. 27, N 2, P. 265.

5. Облака и облачная атмосфера. / Справочник под ред. Мазина И.П. и Хргиана А.Х., Л.: Гидрометеоиздат, 1989.

6. Волковицкий О.А., Павлова Л.Н., Петрушин А.Г. Оптические свойства кристаллических облаков. Л.: Гидрометеоиздат, 1984.

7. Соболев В.В. Рассеяние света в атмосферах планет. М.: Наука, 1972.

8. Ермаков С.М., Михайлов Г.А. Курс статистического моделирования. М.: Наука, 1976.

9. Веселов Д.П., Лобанова Г.И., Попов О.И. и др. // Оптический журнал, 1997, т. 64, №10, С. 48.

10. Веселов Д.П., Лобанова Г.И., Попов О.И. и др. // Исследования Земли из космоса, 1998, №1, С. 38.

11. NOAA KLM USER'S GUIDE (http://www2.ncdc.noaa.gov/docs/klm/index.htm)