НАУКА и ОБРАЗОВАНИЕ

УДК 528.88, 528.714.1, 528.852.4

Россия, МГТУ им. Н.Э. Баумана

victoria.mayorova@gmail.com

alexm.fighter@gmail.com

dim.gr@mail.ru

igorzha@mail.ru

lv-05@mail.ru

toporkov@student.bmstu.ru

a.harlan@hotmail.com

Введение

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

–        объективностью, выражающейся в представлении данных независимо от воли заинтересованных лиц;

–        масштабностью, определяемой техническими возможностями камеры, существенно бóльшей по сравнению с потенциалом аэрофотосъёмки;

–        оперативностью, зависящей от состава орбитальной группировки КА ДЗЗ и составляющей срок от одного дня до недели.

Среди всего многообразия хозяйственной деятельности, осуществляемой в Российской Федерации, сельское хозяйство (c/х) является приоритетным направлением [1]. Огромная площадь и требования к повышению качества обработки земель определяют прямую необходимость внедрения космического сегмента в процесс контроля состояния посевов.

В настоящей работе представлены результаты проведённых исследований, направленных на установление возможности применения космической мультиспектральной и гиперспектральной съёмки для прогноза состояния с/х полей России и дана оценка точности полученных результатов моделирования.

Объекты исследования

            Для отработки моделей прогнозирования были выбраны четыре экспериментальных района (рис. 1): юг Московской области (р-н с. Мочилы), центр Московской области (р-н п. Михнево), Новосибирская область (р-н Бердска, с. Улыбино), Ростовская область (р-н г. Морозовска).

Рисунок 1. Районы прогнозирования

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

Экологическая обстановка в районе с. Мочилы оценивается, как благоприятная, а степень загрязнённости почвы и атмосферы – низкая по сравнению с территориями, расположенными ближе к Москве [2]. Экологическая ситуация в районе п. Михнево характеризуется бóльшей напряженностью по сравнению с экологической ситуацией в районе с. Мочилы, степень загрязненности и изменения природной среды здесь выше. К югу от Новосибирска экологическая ситуация неплохая. Уровень вредных выбросов здесь низкий, а уровень модернизации промышленности позволяет говорить о низкой степени ее влияния на экологию региона [3]. Выбранная территория в Ростовской области подвержена среднему антропогенному воздействию, экологическое состояние оценивается как напряжённое, что ниже среднего по области [4].

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

Выбранный район Новосибирской области характеризуется ярко выраженным лесостепным ландшафтом, представляющим собой наличие с/х полей на распаханных территориях лесных массивов. Таким образом, с/х посевы расположены на локально укрытой от действия ветров местности, что обуславливает переменную чувствительность посевов к солнечной активности, устойчивость концентрации хлорофилла в листе хорошая. Кроме того, укрытость территории и соответствующий почвенный состав [6] (подзолистые, серые лесные и чернозёмные) приводит к достаточному поглощению осадков почвой.

Московская область (ближнее Подмосковье) характеризуется высокой индустриальной насыщенностью. Таким образом, с/х посевы расположены на локально-укрытой местности рядом с промышленными и инфраструктурными объектами, что обуславливает переменную чувствительность посевов к солнечной активности, устойчивость количества хлорофилла в листе хорошая.

Методика исследования

При выполнении исследования была предложена модель прогнозирования вегетационной активности, основанная на использовании относительного разностного индекса вегетации (NDVI). Общий подход к моделированию динамики NDVI состоит в установлении функциональной зависимости между метеоданными и индексом. После установления однозначной связи эти зависимости предполагается использовать для предсказания значений NDVI на некотором интервале времени. При этом в качестве аргументов выступают начальное состояние системы, прогнозные значения метеовеличин и характеристик, определяющих специально вводимые параметры влияния. С целью апробации разработанного методического и алгоритмического обеспечения был проведён спутниковый мониторинг опытного участка – с. Мочилы, который выбирался из соображений удалённости от крупных промышленных объектов, с тем, чтобы свести влияние промышленных выбросов на результаты моделирования к минимуму. Полученные предварительные результаты описаны в [7].

Рисунок 2. Различные подходы к формированию начальных данных для моделирования

На данном этапе исследований были сформулированы три подхода к формированию начальных данных для моделирования (рис. 2). Изначально использовались данные мультиспектральных сенсоров системы LandSat, а климатические данные вводились на основании информации с наземных метеостанций (схема М-Н). Наземные метеостанции имеют некоторую географическую дискретность расположения и, таким образом, не совсем объективно отображают реальное климатическое состояние интересующего района. Кроме того, эти станции предоставляют информацию о состоянии приземного воздуха, которое характеризуется быстрой изменчивостью и не отражает в полной мере суммарный эффект воздействия климата на процесс вегетации, хотя нельзя не отметить, что значения параметров воздуха являются критичными для c/х культур в случае погодных аномалий. В качестве альтернативы было предложено использовать данные с гиперспектрального сенсора в виде температуры подстилающей поверхности (ошибка определения около 250 м) и индекса влагосодержания, устойчивого к атмосферному влиянию. Полный переход на гиперспектральный источник информации (индекс NDVI также рассчитывался по гиперспектральным снимкам) не принёс положительных результатов (схема Г-Г), однако комбинация мультиспектральных данных по NDVI и гиперспектральных данных по климатическим параметрам (схема М-Г) позволила заметно улучшить качество прогнозирования.

В исследовании использовались гиперспектральные снимки с КА Aqua и Terra. Снимки скачивались с сайта NASA[8] за даты, совпадающие с датами снимков с КА LandSat, по которым ранее рассчитывались значения индекса вегетации NDVI. Район съёмки, в котором лежат с/х поля: между 38° и 40° в.д. и между 54° и 55° с.ш.

Стандартные средства программы ScanExImageProcessor® [9] дают возможность автоматически получать маску LST –  значения температуры подстилающей поверхности. Включаемая географическая сетка позволяет безошибочно определить район наблюдений и увеличить изображение до требуемого функционального уровня. Ориентируясь на цветовую раскраску, можно выделить области примерно равных температур, а исходя из относительной малости района наблюдения и его хозяйственной однородности, можно считать колебания температуры поверхности незначительными. Для приведения их к среднему значению предлагается определять квазиминимальную и квазимаксимальную температуру (у каждого участка поверхности, соответствующего пикселю, своё значение температуры и просмотреть каждый не представляется возможным и не имеет смысла), а затем, пользуясь небольшими различиями между ними, вычислять значение температуры подстилающей поверхности для района как среднее арифметическое квазиминимального и квазимаксимального значений. Вышеописанный алгоритм приведёт к получению температуры подстилающей поверхности именно в интересующей точке на поверхности Земли. Физический смысл получаемого параметра говорит о том, что при отсутствии всякого рода термальных аномалий он должен находиться в некоторой взаимосвязи (рис. 3б) с температурой приземного воздуха, данные по которой поступают с наземных метеостанций.

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

Рисунок 3. Предварительные данные по с. Мочилы: а) изменение во времени индекса NDVI; б) изменение во времени температуры приземного воздуха и температуры подстилающей поверхности по данным с КА Terra и Aqua; в) изменение во времени влажности приземного воздуха и индекса NDWI, рассчитанного по данным с КА Terra и Aqua 

Одним из таких доступных для расчёта по данным MODISиндексов, отражающих влагосодержание в растениях и в почве, является Нормализованный разностный водный индекс (Normalized Difference Water Index, NDWI), введенный впервые в 1996 г. Гао (Gao) и определяемый по аналогии с NDVI как [10]:

                                                              (1)

где

NIR– участок ближнего инфракрасного диапазона с длинами волн в интервале  0.841 - 0.876 нм;

SWIR – участок диапазона с длинами волн в интервале 1.628-1.652 нм.

Функциональность формулы объясняется следующими соображениями: вместо использования красного диапазона, интенсивность отражения в котором определяется наличием хлорофилла, применяется коротковолновый ближний инфракрасный диапазон (SWIR), в котором происходит высокое поглощение света водой. В общем случае возможно применение более широкого интервала 1500-1750 нм. Использование же ближнего инфракрасного диапазона (NIR) как и в случае с NDVIсвязано с тем, что вода не поглощает этот участок электромагнитного спектра, тем самым индекс является устойчивым к атмосферному воздействию, что выгодно отличает его от NDVI [11]. Отметим, что применительно к лесным массивам известно, что для индекса NDWI характерно более устойчивое снижение значений после достижения некоторых критических значений уровня антропогенной нагрузки, что может служить его характеристикой, как более чувствительного, по сравнению с NDVI, индикатора экологического состояния лесов. Из сравнения динамики обоих спектральных индексов можно также предположить, что реакция лесов на возрастание антропогенной нагрузки, проявляющаяся в снижении концентрации хлорофилла происходит с некоторой задержкой после начала процессов обезвоживания древесной растительности [12].

Из рис. 3в видно, что наиболее точное совпадение с динамикой данных, полученных с наземных метеостанций, в случае с. Мочилы показывает индекс NDWI, рассчитанный по снимкам КА Terra. График индекса повторяет все узловые точки графика относительной влажности воздуха. Индекс, полученный с КА Aqua совпадает с наземными данными лишь частично. Это может быть связано с тем, что в летние месяцы влага, содержащаяся в растениях и в почве, в течение дня испаряется в условиях достаточной температуры, но к моменту пролёта КА Terra процесс ещё не успевает набрать полную силу [13].

Математическая модель прогнозирования

Учитывая то обстоятельство, что эмпирически найденная кривая изменения индекса NDVI во времени носит нелинейный характер, для осуществления прогнозирования необходимо использовать аппарат дифференциального исчисления. В случае прогнозирования на короткий срок участок кривой принято в самой простейшей схеме заменять на отрезок касательной, считая функцию локально-линейной на рассматриваемом интервале прогнозирования, для чего изначально и подбирались аналитические зависимости с использованием программного приложения Formulize.exe[7]. Однако подобный подход является неприемлемым для данного случая, так как интервалы наблюдений объекта с использованием мультиспектральных сенсоров системы LandSat составляют от недели до двух, что делает шаг по времени сравнимым со всем промежутком наблюдений. В этом случае выходом является либо использование имеющих более высокий порядок слагаемых ряда разложения, либо добавление в модель возмущающих факторов с весовыми коэффициентами. Вместе с тем, влияние возмущающих факторов (в данном исследовании – влажность и температура), может быть заложено как постоянной величиной (весовой коэффициент), так и специальной зависимостью (весовая функция). В данном исследовании предполагается, что изменения температуры и влажности пропорциональны влиянию этих факторов на NDVI с точностью до некоторых констант. 

Примем справедливой зависимость следующего вида:

                         (2)

где

t^* – момент времени прогноза,

t₀ – начальный момент времени,

N(t₀) – априорно известное значение NDVI в момент времени t₀,

N(t^*) – искомое значение NDVI к моменту времени прогноза,

 – текущее время,

T(t) – прогнозное значение среднесуточной температуры,

H(t) – прогнозное значение среднесуточной влажности,

, – некоторые функции для определения зависимости NDVI от среднесуточной температуры и среднесуточной влажности:

                                                     (3)

                                                    (4)

Методика и алгоритм определения функций  и –  следующие. На выбранном интервале времени   получают данные о NDVI и метеоданных. Предпочтительной является информация на небольшом интервале времени  (не более двух недель), с большим количеством (n) измерений параметра по назначению (не менее 4 измерений). Измерению  соответствует значение N_i в момент времени t_i.

В силу того, что исходных данных по N(t) ограниченное количество, для определения зависимостей , где j=T, Н, в первом приближении делается допущение о линейном характере  в диапазоне метеоданных (текущие значения определяются на момент 13:00 местного времени):

,                                      (5)

, ,                                  (6)

Зависимость функций  и только от одного соответствующего метеопараметра объясняется использованием допущения о суперпозиции. Если в выражении (2) для N(t^*) стандартного ряда разложения изначально прогнозируемый параметр по назначению зависит формально только от времени прогноза, то с учётом (5) и (6) N есть уже функция от T, H и t, где T, H– неявные функции времени.

Для четырёх измерений N (i=4), то есть трёх интервалов времени, составляется система линейных уравнений относительно ,  и которая решается одним из известных методов:

Описанные выше алгоритмы были реализованы в программе «ПУСК» (Прогнозирование урожайности сельскохозяйственных культур), предназначенной для вычисления прогнозируемых значений индекса NDVI и соответствующих ему в зависимости от типа исследуемой сельскохозяйственной культуры предполагаемых значений урожайности.

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

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

На начальном этапе предполагалось, что функциональная зависимость индекса NDVI от времени будет получена по имеющимся измерениям с космических снимков с помощью обработки в программных приложениях, позволяющих осуществлять подбор аналитической зависимости по экспериментальным данным. Были проведены предварительные исследования и показана теоретическая возможность подбора таких зависимостей с приемлемым качеством аппроксимации, причём с использованием стандартных полиномов, со степенью не выше шестой [7]. Построенные кривые частично подтвердили предположение о возможности классификации c/х культур по геометрическим характеристикам кривой NDVI (по углам наклона восходящей и нисходящей ветвей графика и по величине максимума) [14]. Вместе с тем, основной идеей выполнения исследований являлось прогнозирование состояния локальной экосистемы, и в этом случае использование аппроксимирующих кривых затруднительно, так как их поведение в межузловых точках  и при экстраполяции зачастую вовсе не соответствует участку с опорными данными [15], следовательно, необходим переход на численные соотношения между достоверными точками и использование этих соотношений в прогнозировании.

В качестве первого варианта модели предлагалась схема Эйлера, в которой производная по времени представлялась в виде частных производных NDVI по температуре и влажности в предположении линейной зависимости каждой из частных производных от соответствующего параметра, зависимость частных производных от времени принималась неявной. Результаты моделирования показали, что при прогнозировании на недельный срок наблюдается сильное расхождение между истинными и рассчитанными данными
(рис. 4). Приближенный к реальности результат показало прогнозирование на срок в 3 дня (рис. 5), вместе с тем, было отмечено, что имеют место вылеты значений индекса NDVI за допустимый диапазон значений, а также локальные колебания относительно истинного значения с большой амплитудой (рис. 6). По итогам спутникового мониторинга наилучший результат показало моделирование, проведённое по схеме М-Г для Ростовской области с интервалом прогноза в 3 дня (рис. 7).

Рисунок 4. Истинная и прогнозная кривые значений NDVI при прогнозировании на недельный срок с использованием изначальной математической модели по схеме М-Н

Рисунок 5. Истинная и прогнозная кривые значений NDVI при прогнозировании на срок в 3 дня с использованием изначальной математической модели по схеме М-Н (Московская область)

Рисунок 6. Истинная и прогнозная кривые значений NDVI при прогнозировании на срок в 3 дня с использованием изначальной математической модели по схеме М-Н (Ростовская область)

Рисунок 7. Истинная и прогнозная кривые значений NDVI при прогнозировании на срок в 3 дня с использованием изначальной математической модели по схеме М-Г (Ростовская область)

Атмосферная коррекция мультиспектральных изображений с КА LandSat

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

Данные дистанционного зондирования с КА серии LandSat, предоставляемые конечным пользователям, проходят определенную обработку, включающую радиометрическую коррекцию и масштабирование [16] полученных значений яркости пикселов на шкалу возможных значений яркости элемента изображения. Эти данные зависят от радиометрического разрешения матрицы (количества уровней яркости) и представляют собой величины, пропорциональные количеству приходящей радиации (т.н. DN - digital numbers).

Таким образом, для того, чтобы провести АК, необходимо предварительно подготовить снимок. Самым важным этапом подготовки снимка к АК является процедура конвертации исходных значений яркости DN снимка (которые изменяются от 0 до 255 для LandSat 5 TM) в реальные значения приходящего излучения. Значения DN безразмерны и пропорциональны количеству попадающего на сенсор излучения, которое измеряется в
Вт/ (м²·стерадиан·мкм).

Перед запуском КА определяется соотношение между DN и измеряемым потоком энергии. Этот процесс носит название калибровки сенсора. Перевод данных из DN в реальные значения приходящего излучения осуществляется с помощью специальных формул [17], главная из которых (8) представлена ниже:  

,                                      (8)

где

L_λ– количество приходящего излучения;

Q_calmax– максимальное калибровочное значение DN (=255);

Q_calmin – минимальное калибровочное значение DN (=0 или 1);

L_maxλ– количество приходящего излучения, которое после масштабирования становится Q_calmax;

L_minλ– количество приходящего излучения, которое после масштабирования становится Q_calmin;

Q_cal  – калибровочное значение DN.

Значения соответствующих величин, которые нужно подставить в (8) и другие важные параметры съёмки обычно распространяются с самими данными LandSat и находятся в так называемом файле метаданных. Пример части файла метаданных LandSat 5 TM представлен ниже

..

  GROUP = MIN_MAX_RADIANCE

    LMAX_BAND1 = 193.000

    LMIN_BAND1 = -1.520

    LMAX_BAND2 = 365.000

    LMIN_BAND2 = -2.840

    ………….

  END_GROUP = MIN_MAX_RADIANCE

  GROUP = MIN_MAX_PIXEL_VALUE

    QCALMAX_BAND1 = 255.0

    QCALMIN_BAND1 = 1.0

    QCALMAX_BAND2 = 255.0

    QCALMIN_BAND2 = 1.0

    ……

  END_GROUP = MIN_MAX_PIXEL_VALUE

  ……………..

Зная параметры и используя формулу (8) можно пересчитать полученные данные с КА LandSat. Важно то, что расчёт проводится для каждого канала отдельно. При пересчете следует учитывать, что для данных TM (LandSat 5)   0 Также нужно обратить внимание, что в зависимости от даты съёмки, коэффициенты разные (см. Таблицу 1) [18].

Таблица 1

Значения калибровочных коэффициентов для проведения АК данных КА LandSat

Для АК снимков с КА LandSat 5/TM, LandSat 7/ETM+ применялся программный модуль FLAASH программного комплекса ENVI [19], при работе с которым использовались значения представленных выше калибровочных коэффициентов и констант. В ходе обработки снимков выяснилось, что в целом динамика кривой NDVIне меняется, зато абсолютные значения индекса увеличиваются на квазипостоянную величину, повышая устойчивость индекса к влиянию климатических изменений (рис. 8). Это происходит вследствие устранения замутнённости снимка и повышения его контрастности. Эффект влияния процессов рассеяния и поглощения излучения атмосферой состоит в том, что на изображении появляется некий аналог атмосферной дымки, приводящий к уменьшению контраста изображений. При прохождении атмосферы некоторое количество излучения рассеивается и преломляется, то есть уменьшается количество излучения, попадающего на датчик оптико-электронной аппаратуры, вследствие чего сцена на снимке является слегка размытой.

На рис. 9 представлены результаты моделирования динамики NDVIс применением АК для того же с/х поля, результаты моделирования для которого были показаны на рис. 7.

Рисунок 8. График зависимости значений NDVI от времени до АК и после АК

Рисунок 9. Реальная и прогнозная кривые значений NDVI при прогнозировании на срок в 3 дня с использованием изначальной математической модели по схеме М-Г (Ростовская область) с применением АК

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

Радиометрическая коррекция

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

Процесс искажения изображения может быть математически представлен в следующем виде [20]:

                                             (9)

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

                                               (10)

В случае работы с растровыми изображениями функции  и  являются  матрицами размерностью M*N, представляющими соответственно искажённое и идеальное изображение в цифровом виде.  будет являться искажающей функцией оптической системы, представленной в виде квадратной матрицы, размерностью D*D, где  должны выполняться условия K<M и K<N. Размерность  определяет размер окрестности пикселей идеального изображения, участвующих в формировании значения одного пикселя искажённого изображения. При этом степень вклада каждого пикселя идеального изображения, принадлежащего этой окрестности, в формирование значения пикселя искажённого изображения целиком определяется видом искажающей функции. В качестве примера на рис.10 приведены графики искажающих функций размытия (а) и смаза (б) размерностью 30x30.

Рисунок 10. Искажающие функции размытия (а) и смаза (б) размерностью 30х30

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

                                                                     ,                                                   (11)

где  – радиус Земли,  – пространственное разрешение датчика ДЗЗ,  – период орбиты спутника ДЗЗ,  – временной интервал дискретизации (время выдержки) датчика ДЗЗ. Правая часть равенства в формуле (11) округляется до ближайшего целого натурального значения, поскольку .

 Фильтрация Винера

Данный метод восстановления изображений соединяет в себе учёт свойств искажающей функции и статистических свойств шума. Фильтрация Винера основана на рассмотрении изображения и шума как случайных процессов. При этом задача алгоритма состоит в том, чтобы найти такую оценку  для идеального (неискажённого) изображения f, чтобы среднеквадратическое отклонение этих величин было минимальным. Среднеквадратическое отклонение задаётся следующей формулой [20]:

                                                              (12)

где M{x}– функция математического ожидания аргумента x. Предполагается, что выполнены следующие условия:

-      шум и неискажённое изображение не коррелированны между собой;

-      либо шум, либо неискажённое изображение имеют нулевое среднее значение;

-      оценка линейно зависит от искажённого изображения.

При выполнении перечисленных условий минимум среднеквадратического отклонения достигается для функции, которая задаётся в частотной области следующим выражением [21]:

где – энергетический спектр шума, – энергетический спектр неискажённого изображения.

 Регуляризация Тихонова

Решение задачи минимизации сглаживающего функционала со связью (регуляризация Тихонова) в частотной области представляется следующим выражением [20, 21]:

где

 – параметр регуляризации, обычно принимающий значения, близкие к 10^-5;

P(u,v) – преобразование Фурье для матрицы регуляризирующих параметров

                                                   (15)

Тестирование алгоритмов

Для проведения численного эксперимента по восстановлению размытого и искажённого изображения с помощью параметрической фильтрации Винера и регуляризации Тихонова был взят фрагмент снимка ДЗЗ, смоделировано его размытие и зашумление с использованием искажающей функции с гауссовым шумом с дисперсией 10^-5, а затем выполнена процедура восстановления изображения с помощью алгоритмов, реализованных в специальной программе [22]  (рис. 11).

Рисунок 11. Пример размытия изображения и восстановления его двумя методами: а) исходное изображение; б) зашумлённое изображение; в) Результат восстановления изображения методом фильтрации Винера; г) восстановление изображения методом регуляризации Тихонова

Для осуществления радиометрической коррекции космических снимков необходимо обладать информацией об искажающей функции датчика ДЗЗ, для построения которой требуется знание параметров конкретного датчика (пространственное разрешение и временной интервал дискретизации), а также параметров орбиты КА. Наиболее приближённый к оригиналу результат восстановления размытого зашумлённого изображения (дисперсия шума 10^-5) показали алгоритмы восстановления изображений методом минимизации сглаживающего функционала (регуляризация Тихонова). Качество работы алгоритмов восстановления зависит от значения дисперсии и математического ожидания самого изображения и шума, а также от соотношения сигнал/шум.

Изменения структуры модели

Несмотря на внесённые изменения в процесс подготовки начальных данных, не была в полной мере достигнута цель – получение адекватных данных прогноза на интервале 3-7 дней. Для участков Ростовской области, где результаты моделирования показали наилучшие результаты, было выдвинуто предположение, что это во многом объясняется степным ландшафтом местности [23], при котором высоко влияние климата на вегетационные процессы. Однако результаты по остальным участкам из Московской и Новосибирской областей свидетельствовали о необходимости доработки модели. В связи с этим были проведены следующие изменения:

1)     введено ограничение на допустимые значения прогнозных величин индекса NDVI:

            |NDVI| ≤ 1;

2)     введён контроль скорости изменения значений NDVI с целью получения адекватной реакции системы на изменение климатических параметров в формате ед./день:

             ед./ день;

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

3)     добавлена компонента, описывающая изменения аэрозольной оптической толщины (АОТ) атмосферы:

Рисунок 12. Представление величины АОТ на сайте NASA в виде палитры значений от 0 до 1

АОТ является безразмерной величиной и определяет количество света, которое рассеивается или поглощается атмосферными частицами во время прохождения света через атмосферу. GoddardSpaceFlightCenterNASA [24] предоставляет ежесуточные данные MOD04_L2 об АОТ с пространственным разрешением 1000 м (рис. 12). Этот тип данных относится как ко второму уровню обработки гиперспектральных изображений, так и, будучи дополненным, к третьему уровню [24]. АОТ представляет собой параметр, чувствительный к оптической прозрачности атмосферы и является индикатором содержания в воздухе различных веществ природного и техногенного происхождения. В атмосфере постоянно присутствуют различные природные частицы, в том числе водяной пар, концентрация которого (АОТ, соответственно) может изменяться в силу климатических вариаций. Вместе с тем, постоянная составляющая аэрозоля в атмосфере и колебания её величины с учётом климата относительно малы по сравнению с антропогенным влиянием, выражающимся, в частности, в выбросах продуктов горения и в распылении химических веществ при обработке полей. Отметим, что и некоторые природные факторы (например, недавнее крупное извержение вулкана в Исландии) могут оказывать аналогичное воздействие на АОТ, но эти случаи оговариваются отдельно;

4)     предложено повышение порядка системы за счёт времени или за счёт климатических параметров:

Вариант 1. Линейная модель с добавлением АОТ. Основание прогноза – 5 временных точек, система линейных алгебраических уравнений (СЛАУ) 4-го порядка:

  Вариант 2. Нелинейная по климатическим составляющим модель. Основание прогноза – 7 временных точек, СЛАУ 6 порядка:         

Вариант 3. Нелинейная по времени модель. Основание прогноза – 8 временных точек, СЛАУ 7 порядка:

(16в)

Основные результаты

В ходе моделирования было установлено, что повышение порядка по времени в варианте №3 (16в) приводит к быстрому расхождению прогнозных данных с реальными. По факту повышение порядка этим методом приводит к фиксированному усилению (так как постоянен шаг по времени) изменения температуры и влажности.

Важным моментом является существенное повышение точности шестидневного прогноза (рис. 13, 14), полученное после всех изменений как для линейной (вар. 1), так и для нелинейной (вар. 2) моделей. В отличие от графика, представленного на рис. 9, прогнозные зависимости на рис. 13 характеризуются бóльшей достоверностью. Нельзя не отметить, что моделирование и оценка расхождений с реальной кривой проводились для так называемой кривой тенденции, отражающей отфильтрованную от каждодневных колебаний динамику изменения значений индекса. При определении опорных точек графика NDVI по имевшимся в наличии мультиспектральным снимкам не было получено неадекватных данных (в отличие от каждодневных данных в случае использования гиперспектрального сенсора КА Terra), что позволило исключить ту часть ошибки прогнозирования, которая определяется некорректностью задания начальных условий. Хорошие результаты показало повышение порядка климатических компонент, вместе с тем заметно повышение колебательности прогнозной кривой (рис. 14). Дать оценку достоверности этой колебательности не представляется возможным, так как используемые мультиспектральные снимки имели период обновления от 7-ми до 21 дня.

Рисунок 13. Результаты моделирования по варианту №1 на 3 дня  и на 6 дней: а) Михнево б) Мочилы в) Новосибирск  г) Ростов

Рисунок 14. Результаты моделирования по варианту №2 на 3 дня  и на 6 дней: а) Михнево б) Мочилы в) Новосибирск  г) Ростов

Оценка точности прогнозирования

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

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

                                                                                (17)

где  – интервал прогноза, количество дней между опорной и начальной датой,

N – количество исходных дат, равное 5 для варианта 1 и равное 7 для варианта 2.

Численные эксперименты были проведены с использованием серий №1 и №2. Далее в зависимости от интервала прогнозирования значения погрешностей для каждой из серий были сведены к одному графику. Как видно из рис. 15, все значения погрешностей расположены ниже некоторой кривой максимальной погрешности:           

Рисунок 15. Кривая максимальной погрешности, полученная по сериям №1 (А_i) и №2 (Б_i)

После осреднения значений погрешности по всем экспериментам для каждого значения ширины интервала прогнозирования и аппроксимации полученной зависимости экспонентой для каждой серии были получены кривые ожидаемой погрешности (рис. 16):

Рисунок 16. Кривые ожидаемой погрешности для линейной (вариант №1) и нелинейной моделей (вариант №2) прогнозирования

Эти законы были записаны в виде эмпирических зависимостей:

                                           (18)

где

– максимально возможная нормированная погрешность прогнозирования;

 – параметр модели, равный 1 для линейной модели, и 0,87 для нелинейной модели (вариант 2);

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

 – величина интервала прогнозирования в днях;

 – слагаемое, учитывающее разреженность основания прогноза.

Закон изменения средней (ожидаемой) погрешности выглядит аналогично (16):

                                                                    (19)

Выводы

1.   Разработаны схемы оптимального сочетания мультиспектральной и гиперспектральной космической съёмки с различным пространственным и временным разрешением в составе методики мониторинга с/х угодий;

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

3.   Проведена серия численных экспериментов по радиометрической коррекции изображений и определён наиболее подходящий для данных задач метод восстановления изображений;

4.   Проанализированы и отобраны основные факторы, влияющие на процесс вегетации;

5.   Получены достоверные данные прогнозирования состояния с/х угодий с  6-дневным интервалом прогноза;

Работа выполнена при финансовой поддержке Министерства образования и науки Российской Федерации в рамках мероприятия 1.9 Федеральной целевой программы «Исследования и разработки по приоритетным направлениям развития научно-технического комплекса России на 2007-2013 годы», государственный контракт от 18 августа 2011 года № 11.519.11.5002, шифр «2011-1.9-519-004».

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

1.               Приоритетные национальные проекты. Совет при Президенте России по реализации приоритетных национальных проектов: Официальный сайт. Режим доступа: http://www.rost.ru/news/2007/10/091614_11170.shtml (дата обращения 15.04.2013).

2.               Экологические карты Московской области. Официальный сайт посёлка Загорянка. Режим доступа:  http://zagoryansky.com/environmental-maps (дата обращения 15.06.2012).

3.               Экологическая карта Новосибирской области. Режим доступа:  http://www.mapnso.ru/  (дата обращения 15.03.2013).

4.               Экологический мониторинг Ростовской области. Донской экологический центр. Режим доступа:  http://ektor.ru/pages/mon1.asp?Idr=2&id=56  (дата обращения 12.09.2012).

5.               Нормативы и методика применения побочной продукции сельскохозяйственных культур для обеспечения бездефицитного баланса органического вещества в почвах на землях сельскохозяйственного назначения. Донской НИИСХ. Режим доступа:  www.don-agro.ru/files/nauka/gosznr/Soloma.doc (дата обращения 23.11.2012).

6.               Природа Новосибирской области. Режим доступа:  http://bober.ru/research/novosibirsk.htm  (дата обращения 15.01.2013).

7.               Майорова В.И., Гришко Д.А., Муравьёв В.В., Топорков А.Г. Использование космических средств наблюдения для мониторинга локальных экосистем // Вестник МГТУ им. Н.Э. Баумана. Сер. Машиностроение. 2012. Спец. вып. «Крупногабаритные трансформируемые космические конструкции и материалы для перспективных ракетно-космических систем». С. 182-192.

8.               Центр оперативной обработки данных ДЗЗ с прибора MODIS.NASA. Режим доступа: http://rapidfire.sci.gsfc.nasa.gov/realtime/  (дата обращения 10.02.2013).

9.               Программное обеспечение.Информационно-технологический центр СканЭкс. Режим доступа: http://www.scanex.ru/ru/software/index.html (дата обращения 20.03.2012).

10.            Remote Sensing and GIS in Agriculture. Vegetation Indices. Science Education through Earth Observation for High Schools (SEOS). Режим доступа: http://www.seos-project.eu/modules/agriculture/agriculture-c01-s03.html (дата обращения 25.06.2012).

11.            Gao B.C. NDWI - A normalized difference water index for remote sensing of vegetation liquid water from space // Remote sensing of environment. 1996. Vol. 58, no. 3. P. 257-266.

12.            Разработка методики региональной экологической оценки состояния лесов по данным спутниковых наблюдений. ГИС-анализ экологического состояния лесов Московской области.Режим доступа:  http://www.wonder-lands.net/places-284-1.html (дата обращения 15.07.2012).

13.            Прогноз погоды. Режим доступа:  http://www.rp5.ru  (дата обращения 06.05.2012).

14.            Пугачева И.Ю., Шевырногов А.П. Изучение динамики NDVI посевов сельскохозяйственных культур на территории Красноярского края и республики Хакасия. Режим доступа: http://d33.infospace.ru/d33_conf/2008_pdf/2/46.pdf   (дата обращения 20.03.13).

15.            Аппроксимация функций. Факультет физики РГПУ им. А.И. Герцена. Режим доступа:  http://physics.herzen.spb.ru/library/01/01/nm_labs/approximation.htm (дата обращения 11.12.2012).

16.            Конвертация данных TM, ETM+ в показатели излучения на сенсоре. GIS-Lab («ГИС Лаборатория»): Географические информационные системы и дистанционное зондирование. Режим доступа:   http://gis-lab.info/qa/dn2radiance.html  (дата обращения 26.09.2012).

17.            Vermote E.F., Vermeulen A. Atmospheric correction algorithm: spectral reflectances (MOD09). Version 4.0. NASA, April 1999. Режим доступа: http://modis.gsfc.nasa.gov/data/atbd/atbd_mod08.pdf  (дата обращения 21.04.2012).

18.            Chander G., Markham B. Revised Landsat-5 TM Radiometric Calibration Procedures and Postcalibration Dynamic Ranges // IEEE Transactions on geoscience and remote sensing. 2003. Vol. 41, no. 11. P. 2674-2677. DOI: 10.1109/TGRS.2003.818464

19.            Руководство по эксплуатации модуля ENVI-FLAASH. Лаборатория ГИТ и ДЗ ИГМ СО РАН. Режим доступа:  http://www.nrcgit.ru/metod/nrcgit_flaash.pdf (дата обращения 26.07.2012).

20.            Гонсалес Р., Вудс Р. Цифровая обработка изображений. М.: Техносфера, 2005. 1072 с.

21.            Шовенгердт Р.А. Дистанционное зондирование. Модели и методы обработки изображений : пер. с англ.  М.: Техносфера, 2010. 560 с.

22.            Майорова В.И., Банников А.М., Зайцев К.И. Математическое моделирование процесса радиометрической коррекции снимков ДЗЗ и сравнительная характеристика алгоритмов восстановления изображений // Вестник МГТУ им. Н.Э. Баумана. Сер. Машиностроение. 2013. Спец. вып. С. 110-121.

23.            Смагина Т.А., Кутилин В.С. Ландшафты Ростовской области. Режим доступа: http://stepnoy-sledopyt.narod.ru/geologia/landshafts.htm (дата обращения: 24.01.2013).

24.            Ежесуточные данные MOD04_L2 Центра оперативной обработки данных ДЗЗ с прибора MODIS. NASA. Режим доступа: http://ladsweb.nascom.nasa.gov/data/search.html  (дата обращения 10.01.2013).