Алгоритм выделения тепловых аномалий подстилающей поверхности Земли по данным спутника NOAA Текст научной статьи по специальности «Физика»

УДК 519.226.3

A.A. Моисеев, К.Т. Протасов

Алгоритм выделения тепловых аномалий подстилающей поверхности Земли по данным спутника NOAA

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

Введение

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

В этой связи возникает необходимость привлечения всех доступных средств оперативного обнаружения пожаров на ранней стадии их развития, что объясняет возрастающую роль в этом спутниковых систем дистанционного зондирования Земли. Космический мониторинг имеет ряд преимуществ по сравнению с авиаразведкой: высокую оперативность, большую площадь охвата земной поверхности и меньшие операционные расходы. На охраняемой территории спутниковые данные служат существенным дополнением к традиционным методам обнаружения, а на неохраняемой — единственным средством мониторинга и оценки последствий лесных пожаров.

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

В качестве наиболее простых и часто применяемых на практике можно привести два алгоритма обнаружения очагов пожаров — Fixed Threshold Techniques и Spatial Analysis Techniques. Здесь для получения решающего правила используются радиационная температура Т3 в спектральном канале 3,75 мкм, Т4 — в канале 11 мкм, а также их разница Т34 = Т3- т4.

В Fixed Threshold Techniques для обнаружения очагов применяются фиксированные пороговые величины Т£шп и Т^1, для которых должны выполняться два условия вида Т3 > Т™п и ТЫ>Т.

В Spatial Analysis Techniques обнаружение события «пожар» осуществляется при выполнении двух условий вида Т3 > ц3 + Сст3 и T3i > jx34 + Сст34, где ц — средние значения Т3 или Т34, а — среднеквадратичное отклонение этих величин, а С — числовой коэффициент в диапазоне значений 1,5-4. Значения ци а вычисляются в окне размером N х N пикселей (N = 3 - 21) вокруг «потенциального» очага.

Постановка задачи

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

Описание метода

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

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

Как известно, собственное излучение объекта является функцией температуры, физических свойств и структурных характеристик излучающей поверхности. Содержательной характеристикой собственного некогерентного излучения следует считать спектральную энергетическую яркость нагретого тела, которая может быть определена по формуле Планка. Энергетическая яркость зависит от температуры, длины волны и коэффициента излучения е^ (Т), называемого спектральной степенью черноты излучающей поверхности при заданной температуре и определенном направлении визирования. Так абсолютно черное тело имеет степень черноты = 1 во всем диапазоне длин волн. Для серого тела в широком диапазоне длин волн 0 < Ех < 1. Для некоторых типов подстилающей поверхности параметр (Т) известен, но реально прибор регистрирует радиояркости их комбинаций. Кроме того, свое искажающее влияние оказывает состояние атмосферы. Это создает проблему в оценивании температуры поверхности Земли по радиояркостям, наблюдаемым в условиях влияния атмосферы.

Будем полагать, что проведена предварительная обработка изображений: устранены геометрические искажения данных; осуществлена географическая (координатная) привязка; «вырезан» фрагмент видеоданных, содержащий территорию Томской области (ТО) и ее окрестности; произведена радиационная коррекция и калибровка данных прибора АУШШ переходом к альбедо для 1-го, 2-го каналов и к термодинамической температуре в 3, 4, 5-м каналах (по обращенной формуле Планка).

Таким образом, имеем поле 1024 х 1024 5-мерных векторов, и для каждого элемента — пикселя этого изображения (поля видеоданных) заданы географические координаты. Кроме того, имеются данные метеослужб о температуре (например, температура почвы или приземного слоя воздуха) ППЗ на момент пролета спутника в некоторых, достаточно равномерно расположенных, пунктах с координатной привязкой. В идеальном случае необходимо знать температуру в контрольных точках для каждого класса (типа) поверхности (воды, поля, пашни, леса и т.п.). Эти данные служат для адаптивной настройки параметров алгоритма на конкретные условия региона с учетом оптической погоды.

Термодинамическую температуру ППЗ прогнозируемого поля будем описывать случайной величиной У е Л1, а радиояркости полей, являющихся источниками прогнозирующей информации, будем описывать случайным вектором Хей*, где Д* — А-мерное евклидово пространство, X = (X1 ,...,Хк)т ; X1— радиояркости 1-го канала прибора АУШШ, ¿ = 1 Т — знак транспонирования. Взаимосвязь прогнозируемой переменной У и вектора X будем описывать функционалом регрессии вида

Если существуют нижеследующие плотности вероятностей случайных величин X и У, то с учетом (1) имеем [1]

величины У; /(х) — плотность вероятности случайного вектора X; /"(у) — плотность вероятности случайной величины У; Ж{у) — интегральная функция распределения У.

m(x) = E(Y/X = х), где Е( ) — оператор математического ожидания, причем £(|У|) < °° .

(1)

(2)

где х е Ек, у е Л1; /(х,г/) — совместная плотность вероятностей случайных вектора X и

Пусть в нашем распоряжении имеется выборка попарно независимых одинаково распределенных случайных величин {(Х8,У4)}" , где п — количество контрольных отсчетов температур по данным метеослужб на момент фиксирования спутниковых наблюдений. В этом случае для вычисления выражения (2) естественно воспользоваться непараметрическими оценками неизвестных распределений по выборочным данным. Заменим неизвестные распределения их непараметрическими оценками ядерного типа, а Р(у) — эмпирической функцией Рп(у), тогда оценка уравнения регрессии (2) примет вид

7ЙЛ(х) = ^ д ' к--

¿П о*-гг)

/=1¡=1

(3)

1=1

1=1

где Л — ширина окна (параметр сглаживания), описываемого функцией Кп (и) = Н'1К(и/к),

в качестве К(.) может быть взято ядро Епанечникова или гауссово ядро.

Накопленный практический опыт использования подобных оценок показывает, что точностные характеристики уравнения регрессии ягЛ(х) в большей степени определяются не формой ядра, а параметром масштаба Л. Ввиду важности параметра Л естественно перейти к векторному параметру Л = (к1, к2)7 ив выражении (3) использовать модифицированное ядро

\-1 / , / 4\

(4)

К'к{и) = (к1)' К\и1 ¡Ну I = 1,

При этом возникает задача оценивания И с учетом наблюдений |(Хг, . Воспользуем-

ся для оценивания Ь методом скользящего контроля, заключающегося в построении модифицированной оценки регрессии тй ДХ; ), в которой на каждом шаге подсчета суммы исключается /-е наблюдение. Критерий качества оценивания Ь зависит от способности предсказать

набор значений {У;}. 1 по наборам подвыборок {(Х;,У;)}_:

ЛЪ) = п1 X [у;. - <У(Х;)]2 Ш(Х;.), ;=1

(5)

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

В другом случае для повышения быстродействия вычислений будем использовать ядро Епанечникова [3] с «внутренней» системой координат, обеспечивающей разворот эллипса рассеяния, согласованный с рассеянием выборочных данных:

N.

V ;=1г=1

(6)

где введена следующая вспомогательная система координат:

Г 1 % _ (О / гр\ 81

и = бх; м[ииг] = см XX вт = - А ; Л = К, , в = ет \оп У

где б — матрица декоррелирующего ортогонального преобразования; М[.] — оператор математического ожидания; X — центрированные наблюдения; А — диагональная матрица собственных значении; а = 3/4^5 ; Ь = а/5.

Особенность модифицированного ядра Епанечникова заключается в том, что единый параметр сглаживания для всех размерностей пространства наблюдений масштабируется по координатам собственными значениями А.г, * =1 ,...,п, оценочной ковариационной матрицы

Ду по выборочным данным, что существенно упрощает задачу оптимизации критерия (5). В случае векторного параметра сглаживания задачу оптимизации (5) удобно осуществить поисковым методом Цыпкина [2] с помощью двухэтапной процедуры оценивания глобального экстремума функционала. На первом этапе в области поиска, представляющей собой много-

к.адрат - - V- - —Я и верхняя о,_ _

1 = 1

параметра сглаживания соответственно, случайно с равномерным распределением бросается точка, из которой осуществляется градиентный спуск, при этом используются поисковые методы локальной оптимизации [2]. С этой целью осуществляется вариация функционала качества (5) по параметрам сглаживания следующим образом. Вычисляются значения приращений функционала (5)

= (</(А + ае^.-Лф + аек))Т; </_(Л,а) = (</(Л - аег ),..^(к - аек))т, (7)

где к — количество параметров Ь, соответствующих количеству компонентов вектора, собранных в вектор параметров Ь = ^.....А*| ; а — скалярный параметр, определяющий величину

( V

поискового шага; е; = 0,...,1,...,0 , 1 = 1,...,к, —базисные векторы поисковых направлений. Оценочное значение градиента вычислим следующим образом:

J+(h,a)-J„(h,a)

2 а

= УЛ(Ъ. ,а), (8)

где У,± — обозначение градиента. Поисковый алгоритм адаптации в рекуррентной форме н

имеет вид

Ь[;] = Ь[у-1]-у[у]7й±е7(Ь[у-1],а[;]). (9)

Выбор поискового а[ ] и рабочего у[-] шага рассмотрен в [2], причем для сходимости процедуры необходимо выполнить у [•] < а [•] .

После того как параметр Ь в выражении (3) для тЛ(х) конкретизирован, уравнение регрессии можно использовать для восстановления значений У по наблюдаемым X и для всего поля видеоданных.

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

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

Цх) = РЫх) + ЯГ1{х),

где fl(x) — распределение экстремальных значений температур (пожары); /о(х) — распределение не экстремальных значений температур; Р, <? — априорные вероятности ситуации А0иА1 соответственно, причем Р + (? = 1. Возникает задача оптимизации квадратичного критерия качества следующего вида:

т = -%{f(Xj)-PfoiXj)-QfiiXj)}2,

т

(10)

где f(Xj) —гистограмма распределения температур изображения, а 9 —вектор неизвестных

параметров, состоящий из компоненты Р и параметров функций плотности f0(x) и f1(х), принадлежащих параметрическим семействам функций. Следует заметить, что задача восстановления компонент смеси имеет решение лишь в случае её идентифицируемости. Это трудно формализуемое и проверяемое условие с геометрической точки зрения означает, что fo(x) и /х(х) должны иметь ярко выраженные моды. Поэтому доля или мера участков с экстремальными значениями температур должна быть достаточно высокой, чтобы функция плотности fl(x) могла проявить свою форму. Задача декомпозиции смешивающего распределения ввиду высокой неопределенности не всегда имеет решение. Необходимо привлечь априорные данные о форме составляющих распределений. Известно, что плотность распределения максимумов п независимых случайных величин в асимптотике растущего числа наблюдений п —» °° I типа (распределение максимальных значений Гумбеля [4]) имеет следующий вид:

где ц— параметр (мода) центра распределения; а— масштаб распределения, причем оценочные математические ожидания р и дисперсия д связаны с ц и а следующим образом: Д = ц + 0,577о , ö = 1,283а. В качестве /0(х) нами было выбрано распределение Джонсона SB с параметрами е (нижняя граница х), X (размах выборки), л, у (параметры формы). Учитывая, что часть параметров можно оценить по выборочным данным [1], фактически вектор неизвестных параметров имел лишь три компоненты и 9 = (Р, Г), у)г, где Т — знак транспонирования. Задача оптимизации критерия (10) так же решалась с привлечением адаптивных методов поиска экстремума (9) [2]. После того как смесь идентифицирована, построим Байесово решающее правило проверки двух гипотез: Н1 (тепловая аномалия) и Н0 (норма). Это решающее правило выявляет на изображении все участки видеоданных, связанных с наличием повышенных температур:

и = arg max {Pf0 (х), Qfx (*)}, (11)

{од}

где и — принимаемое решение или номер принимаемой гипотезы, и е {0,1}.

Приведенный алгоритм реализован на языке программирования С++ Builder 6.0.

Результат работы программы

Для тестирования программы взят фрагмент снимков спутника NOAA Иркутской области. Фрагмент вырезан и сохранен в формате bmp при помощи программы Scan Viewer ИТЦ Сканэкс (рис. 1).

Первым действием программы будет оценивание температуры подстилающей поверхности Земли с использованием значений реальных данных с наземных метеостанций. Так как найти реальные данные с метеостанций очень сложно, то для апробации алгоритма были взяты значения радиояркостей 3-го канала в определенных точках вырезанного фрагмента и при помощи формулы Планка пересчитаны на термодинамические температуры [5]:

о

о <т

-<» < х < > а > 0,

с

1

(12)

где Cj = 3,7415 • 108 Вт • мкм4/м2, с2 = 14388 мкм К, к = 3,74 мкм.

.ГТ7Т

.............. .............р

¿Шй'<ч

I

п

. \ Лд ' р •

'V- - ■ -1 . =■ * ? -„ ¿У-

.» -V V

-аднр* * • "V-

Ш

м 1 -о »4

шшшяш

а б

Рис. 1 — Исходное изображение (альбедо)

Полная энергия во всем интервале длин волн описывается формулой

со

\ В(\,Т)(1Х = аТЛ,

(13)

где а = 5,67 • Ю-8 Вт • м 2 • К 4. Таким образом, зная значение радиояркости в точке, значение термодинамической температуры в ней можно получить по формуле

-1п

Х5В

+ 1 -

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

/, - Л-;/

ГГАЛШ Лт у-.;

Коорл>чата Х:347 Координата __^Температура 330

Рис. 2 — Поле восстановленных температур ППЗ

После восстановления поля температур подстилающей поверхности Земли можно выделить тепловые аномалии, т.е. точки земной поверхности, в которых температура превышает допустимое пороговое значение. Аномалии выделены черным цветом (рис. 3).

Коордипата Х;253_Коордипатз У:227 _Тсутсратура 326_

Рис. 3 — Выделенные аномалии

Заключение

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

Литература

1. Хардле В. Прикладная непараметрическая регрессия / В. Хардле ; пер. с англ. - М. : Мир, 1993. - 349 с.

2. Цыпкин Я.З. Адаптация и обучение в автоматических системах / Я.З. Цыпкин. - М. : Наука, 1968. - 400 с.

3. Епанечников В.А. Непараметрическая оценка многомерной плотности вероятности. Теория вероятностей и ее применение / В.А. Епанечников. - 1969. - Т. 14, вып. 1. -С. 156-161.

4. Хан С. Статистические модели в инженерных задачах / С. Хан, С. Шапиро ; пер с англ. - М. : Мир. - 1969. - 396 с.

5. Кашкин В.Б. Дистанционное зондирование Земли из космоса. Цифровая обработка изображений : учеб. пособие / В.Б. Кашкин, А.И. Сухинин. - М. : Логос. - 2001. - 264 с.

Моисеев Артем Андреевич

Аспирант каф. компьютерных систем в управлении и пректировании ТУСУРа

Телефон: 89138060636

Эл. почта: artem0606@mail2000.ru

Протасов Константин Тихонович

Д-р техн. наук, ст. научн. сотр. Института оптики и атмосферы СО РАН Телефон: (3822) 49 11 11 доб. 1264 Эл. почта: prot@iao.ru

А.А. Moiseev, К.Т. Protasov

An algorithm of background thermal anomalies discrimination based on NOAA satellite data

In the paper, a non-parametric regression method of detection of forest fires based on NOAA satellite data is suggested. A temperature field of the Earth's background is estimated by using data of ground weather stations on surface temperature in separate points.