Crop yield forecasting based on the satellite monitoring of carbon dynamics in terrestrial ecosystems
- Authors: Lysenko S.A.1
-
Affiliations:
- Institute for Nature Management of National Academy of Science of Belarus
- Issue: No 4 (2019)
- Pages: 48-59
- Section: Использование космической информации о Земле
- URL: https://journals.eco-vector.com/0205-9614/article/view/15714
- DOI: https://doi.org/10.31857/S0205-96142019448-59
- ID: 15714
Cite item
Full Text
Abstract
A low-parametric model of crop biomass dynamics using the data of satellite remote sensing of terrestrial vegetation index and routine meteorological observations is developed. Modeling is performed with a data filtered by yield-correlation image masking technique. The simulation is based on the Monteith equation for carbon dynamics in terrestrial ecosystems. Meteorological parameters that have an effect on the photosynthesis intensity but rarely measured directly are calculated on the basis of analytical parameterizations obtained from reanalysis data. A validity of the model is demonstrated by the example of satellite monitoring of the spring wheat yield in the regions of the Republic of Belarus.
Full Text
МЕТОДЫ ПРОГНОЗИРОВАНИЯ УРОЖАЙНОСТИ СЕЛЬСКОХОЗЯЙСТВЕННЫХ КУЛЬТУР
Урожайность сельскохозяйственных культур (СХК) — один из важнейших показателей эффективности всего агропромышленного комплекса, учитываемый при планировании импортно-экспортных операций на рынке сельскохозяйственно продукции. Стандартные методы сбора информации об урожайности на уровне страны и ее административных единиц часто субъективны и не обеспечивают необходимой оперативности и заблаговременности. С развитием технологий дистанционного зондирования Земли из космоса появилась возможность проводить мониторинг состояния посевов на обширных территориях оперативно и с минимальными финансовыми затратами (Atzberger, 2013; Клещенко, 2015; Савин и др., 2010; Сладких и др, 2016).
К настоящему времени в мире накоплен большой опыт использования спутниковой информации в агрометеорологии. Физической основой спутникового мониторинга СХК является тот факт, что многие структурные и биохимические изменения растений в течение вегетационного периода проявляются в спектрах отражения ими солнечного излучения, регистрируемых фотоприемником на спутниковом носителе. Измеряемые характеристики отражения растительных покровов представляются либо в виде безразмерных вегетационных индексов (NDVI, EVI и др.) (Huete et al., 2014), либо в виде параметров моделей переноса излучения (листовой индекс (LAI), доля поглощаемой солнечной радиации в фотосинтетически активной области спектра (FAPAR), концентрация хлорофилла в листьях и др.) (Baret, Buis 2008; Ganguly et al., 2014). Спутниковый мониторинг их изменений в течение вегетационного периода в совокупности с данными наземных метеонаблюдений открывают большие возможности для разработки методов агрометеорологических прогнозов.
Прогноз урожайности СХК включает экстраполяцию тренда, обусловленного преимущественно технологическими факторами, и прогноз отклонения от него. Экстраполяция временного ряда урожайности выполняется хорошо проработанными и достаточно надежными методами (наименьших квадратов, гармонических весов, экспоненциального сглаживания и др.). Основную сложность представляет прогноз отклонения урожайности от тренда, которое определяется большим числом метеорологических факторов. Учет их влияния на формирование урожая и является главной задачей методов агрометеорологических прогнозов.
Простейшие методы прогноза урожайности СХК основаны на нахождении года-аналога, в который вегетационные индексы и метеорологические условия формирования урожая были наиболее близки к их значениям, наблюдаемым в текущем году. Данный метод является достаточно грубым и требует длительного ряда спутниковых наблюдений.
Другой класс методов прогноза урожайности составляют методы, основанные на различных регрессионных и нейросетевых моделях, в которых предикторами выступают интегральные агрометеорологические характеристики (суммы температур и осадков, гидротермический коэффициент и др.) и (или) вегетационные индексы за выбранный отрезок времени (Ferencz et al., 2014; Zhang et al., 2005; Becker-Reshef et al., 2010; Куссуль и др., 2012; Страшная и др., 2014; Kouadio et al., 2014). Такие модели наиболее просты в использовании, однако их применимость ограничена достаточно узким диапазоном почвенно-климатических условий, для которых осуществлена идентификация параметров модели.
Наиболее надежные прогнозы урожаев составляются на основе моделирования биопродуктивности растений с учетом основных физических и биологических процессов в системе «почва–растение–атмосфера» (фотосинтез, дыхание, рост и отмирание органов, приток углерода в почву, почвенное дыхание, эвапотранспирация и др.) (Сиротенко, 1981; Полуэктов и др., 2006; Delécolle et al., 1992; Kasampalis et al., 2018). Настройка (инициализация) моделей биопродуктивности предполагает сравнение результатов моделирования с натурными данными в границах каждого поля и подбор ряда эмпирических констант. Такие модели обладают наибольшей универсальностью и детальностью учета метеорологических условий формирования урожая, однако требуют задания большого количества метеорологических и агрохимических характеристик, контроль которых для каждого поля является трудно выполнимой задачей. При этом многие агрометеорологические параметры (солнечная радиация, влажность почвы и др.), используемые в расчетах продукционного процесса растений, не относятся к числу стандартных параметров, измеряемых на сети метеостанций.
При сочетании моделей продукционного процесса растений с моделями переноса излучения и данными дистанционного зондирования Земли (ДЗЗ) появляется возможность корректировки и уточнения модельных параметров путем сопоставления расчетных отражательных характеристик растительного покрова с результатами их дистанционных измерений (Delécolle et al., 1992; Rembold et al., 2013; Jin et al., 2016). Альтернативный способ усовершенствования динамического моделирования биомассы растений состоит в замене расчетных параметров модели на параметры, определяемые по измерениям из космоса (Rembold et al., 2013; Клещенко, Найдина, 2011; Брыскин, 2007). Наиболее часто такие замены выполняется для блока модели, отвечающего за фотосинтез, включая такие параметры как LAI и FAPAR. При недостаточных возможностях спутникового сенсора для их определения, вместо них используются спектральные вегетационные индексы, косвенно отражающие сезонный ход фотосинтеза и фенологию растений.
Для расчетов динамики углерода в наземных экосистемах по ДЗЗ широко используются методы, основанные на уравнении Монтейта (Савин и др., 2010; Liu et al., 2010; Duchemin et al., 2008; Rembold et al., 2013; Xiao et al., 2014). Данное уравнение описывает прирост биомассы за заданный отрезок времени в следующей функциональной форме (Monteith, 1972):
(1)
где i — номер расчетного периода, NPPi (Net Primary Production) — чистая первичная продукция органического вещества (гр. сухой биомассы/м2/период), PARi (Photosynthetically Active Radiation) — энергия солнечного излучения, поступающего за заданный отрезок времени на подстилающую поверхность в фотосинтетически активной области спектра, FAPARi (Fraction of Absorbed Photosynthetically Active Radiation) — доля фотоситнетически активной радиации, поглощаемая пигментами зеленого листа, LUE (Light Use Efficiency) — эффективность преобразования поглощенной растительностью фотоситнетически активной радиации в органическое вещество, зависящая от ряда метеорологических параметров, CO2 — концентрация углекислого газа в атмосфере, которая на протяжение вегетационного периода практически не изменяется, Ti — температура приземного воздуха, VPDi (Vapour Pressure Deficit) — дефицит давления водяного пара (разница между давлением насыщенного пара и фактическим парциальным давлением пара при данной температуре), Wi — содержание влаги в деятельном слое почвы, Ri — доля затрат органического вещества на дыхание живой биомассы.
При использовании уравнения (1) для расчетов урожаев СХК делается предположение, что фитомасса урожая Y пропорциональна общей биомассе, накопленной с начала вегетационного периода, т. е.:
(2)
где N — общее количество расчетных периодов, C — коэффициент пропорциональности, определяемый путем сопоставления расчетов биомассы с фактической урожайностью.
Интегрирование приростов биомассы с учетом складывающихся термических и гидрологических условий вегетационного периода, как правило, позволяет улучшить точность оценок урожая по сравнению с эмпирическими регрессиями между урожаем и одномоментными значениями вегетационных индексов (Савин и др., 2010; Rembold et al., 2013). При этом, однако, возникает та же проблема, что и для методов динамического моделирования биопродуктивности растений: необходимость получения обучающих (калибровочных) данных отдельно для каждого поля. К тому же, при комплексировании методов расчета динамики биомассы с технологиями ДЗЗ требуется построение схемы размещения сельскохозяйственных полей с идентификацией, выращиваемых на них культур (Терехов и др., 2010; Плотников и др., 2011; Сладких и др., 2016; Hansen et al., 2000; Becker-Reshef et al., 2010), что на практике не всегда возможно, поскольку доступ к архивам сельскохозяйственных организаций с необходимой информацией в большинстве случаев отсутствует.
Для целей мониторинга состояния СХК на больших территориях (на уровне отдельных государств или в глобальном масштабе) наиболее предпочтительно использование спутниковых сенсоров среднего спектрального разрешения (~0.25–1.0 км) (Rembold et al., 2013). Однако в этом случае на спутниковых изображениях можно идентифицировать только достаточно крупные поля, площадью в десятки гектаров. Многие сельскохозяйственные поля Европы не отвечают этому требованию, что затрудняет применение к ним описанных выше технологий динамико-статистического моделирования биомассы растений.
Преодолеть проблему выделения сельскохозяйственных полей на спутниковых изображениях можно с использованием корреляционного метода, предложенного в работе (Kastens et al., 2005). В основе данного метода лежит предположение, что метеорологические условия текущего вегетационного периода схожим образом сказываются на эффективности продукционного процесса растений, как в границах сельскохозяйственных полей, так и за их пределами. Это позволяет при составлении агрометеорологических прогнозов рассматривать на аэрокосмических изображениях все пикселы, отображающие растительные покровы в исследуемом регионе, которые могут быть не менее показательны относительно потенциальной урожайности, чем пикселы, относящиеся непосредственно к сельскохозяйственным полям.
Техническая реализация корреляционного метода выделения информативных пикселов на аэрокосмическом изображении состоит в следующем. На основе многолетней серии спутниковых изображений и статистических данных по урожайности конкретной культуры в конкретном регионе вычисляется коэффициент корреляции между вегетационными индексами всех пикселов спутникового изображения и урожайностью. Задается пороговое значение коэффициента корреляции. Пикселы, для которых коэффициент корреляции между вегетационным индексом и урожайностью меньше этого порога, считаются неинформативными и исключаются из рассмотрения. Пикселы с коэффициентом корреляции, превышающим пороговое значение, усредняются и используются для прогнозирования урожайности в данном регионе. На примере ведущих сельскохозяйственных штатов США показано (Kastens et al., 2005), что использование корреляционной маски для осреднения вегетационных индексов обеспечивает примерно такую же точность прогнозирования урожайности, как и использование маски непосредственно пахотных земель. При этом очевидно, что методы прогнозирования урожайности, не требующие информации о расположении пахотных земель на исследуемой территории, являются наиболее универсальными и удобными в применении, поскольку структура посевных площадей многих государств подвержена ежегодной ротации (меняются как возделываемые культуры, так и поля, отведенные под них).
В настоящей работе описывается метод прогнозирования урожайности СХК на основе ДЗЗ и наземной метеорологической информации, использующий принципы построения корреляционных масок под конкретные культуры и моделирование их биопродуктивности на основе уравнения Монтейта с регионально адаптированными параметрами. Прирост биомассы рассчитывается с суточным разрешением с учетом потока солнечной радиации на уровне подстилающей поверхности, температуры воздуха, осадков, влажности воздуха и почвы. Учитывая практические сложности получения информации о потоке солнечной радиации и влажности почвы, они рассчитываются с использованием аналитических параметризации, полученных на основе данных реанализа Европейского центра среднесрочных прогнозов погоды. Эффективность метода демонстрируется на примере прогнозирования урожайности яровой пшеницы в областях Республики Беларусь (РБ).
СТАТИСТИЧЕСКИЙ АНАЛИЗ ИНФОРМАТИВНОСТИ ПИКСЕЛОВ АЭРОКОСМИЧЕСКОГО ИЗОБРАЖЕНИЯ
Учитывая отмеченные выше сложности идентификации пахотных земель на аэрокосмических изображениях, нами используется корреляционный метод выделения информативных пикселов. В качестве спутникового прибора для мониторинга состояний посевов и прогноза урожайности выбран американский спектрорадиометр среднего пространственного разрешения MODIS, имеющий достаточно длинный ряд наблюдений (2000–2018 гг.).
Полагаем, что состояние посевов комплексно характеризуется нормализованным разностным вегетационным индексом NDVI, рассчитываемым с пространственным разрешением 500 м по коэффициентам яркости подстилающей поверхности в красном (620–670 нм) и ближнем инфракрасном (841–876 нм) каналах прибора MODIS. В отличие от таких информационных продуктов MODIS как LAI (Leaf Area Index) и FAPAR (Fraction of Absorbed Photosynthetically Active Radiation) для определения NDVI не требуется использования априорных предположений об архитектуре растительного покрова и моделей переноса излучения, что делает его менее подверженным ошибкам измерений и расчетов. За счет этого случайная компонента во временной динамике NDVI, связанная с неточностями атмосферной коррекции спутниковых снимков, значительно менее выражена, чем случайная компонента LAI и FAPAR.
Данные по NDVI, получаемые прибором MODIS, находятся в свободном доступе и предоставляются в виде 16-дневных композиционных карт, построенных по всем доступным за это время снимкам с учетом следующих требований к каждому пикселу: минимальная облачность, наименьший угол визирования, наибольшее значение NDVI. Для целей мониторинга посевов используются композиционные карты NDVI за восемь 16-дневных периодов с конца марта по конец июля.
Для Республики Беларусь данные по урожайностям основных СХК имеются отдельно для каждой области и административного района с 1960 г. Анализ этих данных показывает наличие трендов урожайности СХК: с 1990 по 2000 гг. — урожайность падала, а после 2000 г. начала расти (рис. 1). Столь резкое изменение трендовой компоненты урожайности с 2000 г. говорит о том, что ее происхождение связано исключительно с технологическими факторами (уровнем техники, сменой сортов, изменениями доз удобрений и др.). Компонента урожайности, остающаяся после вычитания линейного тренда, определяется преимущественно метеорологическими факторами. Ниже, говоря об урожайности, будем подразумевать именно эту компоненту, прогнозирование которой и составляет цель большинства агрометеорологических исследований.
Рис. 1. Динамика урожайности зерновых и зернобобовых культур по областям Республики Беларусь с 1990 по 2017 гг. (пунктир с символами) и линейные тренды урожайности до и после 2000 г. (прямые).
Для выделения на картах NDVI пикселов, содержащих информацию о потенциальной урожайности той или иной СХК, рассчитываются коэффициенты корреляции между урожайностью этой культуры и 16-тидневными значениями NDVI. Такие корреляции рассматривались отдельно для каждой из шести областей РБ и каждого из восьми внутрисезонных композитов NDVI с учетом только тех пикселов, которые попадают в границы рассматриваемой области. Из коэффициентов корреляции, рассчитанных для каждого пиксела, выбирается максимальный по модулю коэффициент, относящийся к одному из 16-дневных интервалов в течение периода вегетации.
Карта максимального коэффициента корреляции между NDVI и урожайностью яровой пшеницы в областях РБ приведена на рис. 2. Положительные значения коэффициентов корреляции соответствуют периоду интенсивного роста растений, сопровождающегося увеличением NDVI. Отрицательные — периоду колошения яровой пшеницы, когда происходит формирование и налив зерна, а значения NDVI уменьшаются.
Рис. 2. Карта коэффициента корреляции между урожайностью яровой пшеницы в областях Республики Беларусь и 16-дневными вегетационными индексами территорий в пределах этих областей.
Пикселы для прогнозирования урожайности выбираются исходя из требования, чтобы соответствующий им коэффициент корреляции между NDVI и урожайностью превышал некоторое пороговое значение, условно выбранное равным 0.5. Значения NDVI, усредненные по маске таких пикселов для каждой области, используются для расчета динамики биомассы СХК в этих областях.
МОДЕЛИРОВАНИЕ ДИНАМИКИ БИОМАССЫ СЕЛЬСКОХОЗЯЙСТВЕННЫХ КУЛЬТУР
В основу моделирования динамики биомассы СХК положено уравнение Монтейта (1). Для его практического использования конкретизируем входящие в него величины. Эффективность использования солнечной радиации для фотосинтеза LUE зададим в виде произведения двух функций η1(T) и η2(W), определяющих зависимость интенсивности фотосинтеза от температуры воздуха T и влажности почвы W.
Функция η1(T) имеет вид одновершинной кривой с максимальным значением η1(Topt) = 1 в температурном оптимуме фотосинтеза Topt, зависящем от типа растительности. Известно несколько функциональных зависимостей η1(T), пригодных для разных типов растительности. Однако, для того чтобы не перегружать модель биопродуктивности слишком большим количеством свободных параметров, которые подлежат установлению на основе ограниченного временного ряда эмпирических данных, будем использовать самую простую температурную зависимость фотосинтеза, учитывающую ее наиболее характерные особенности (Сиротенко, 1981):
, (3)
где a — положительный безразмерный параметр.
Для функции η2(W) также выбрано наиболее простое выражение с единственным свободным параметром (Дмитренко, 2010):
, (4)
где W и Wopt — фактическая и оптимальная влажность метрового слоя почвы.
Затраты биомассы на дыхание рассчитываются в зависимости от температуры воздуха по следующей формуле (Сиротенко, 1981):
(5)
где TR — базисная температура дыхания, превышение которой температурой воздуха на 10 °C сопровождается увеличением расходов органического вещества на дыхание в два раза.
В связи со сложностью и высокой стоимостью полевых измерений влажности почвы на больших территориях, данные по ней, как правило, относятся лишь к отдельным опытно-экспериментальным участкам и недоступны в оперативном режиме. Тем не менее, ее можно рассчитать на основе стандартных метеорологических параметров (таких как температура, осадки, влажность воздуха и скорость ветра) с использованием моделей переноса влаги в системе «почва–растение–атмосфера» (Сиротенко, 1981; Полуэктов и др., 2006). Однако при мониторинге состояния и развития СХК из космоса такое моделирование необходимо проводить отдельно для каждого пиксела спутникового снимка, что приводит к неоправданно большим вычислительным затратам. Выходом из ситуации может быть использование аналитических параметризаций зависимости влажности почвы от стандартных метеорологических параметров. Для получения таких параметризаций необходимы репрезентативные данные натурных наблюдений или расчетные данные моделей переноса влаги.
Нами для расчета влажности метрового слоя почвы использовались аналитические параметризации, полученные на основе данных реанализа Европейского центра среднесрочных прогнозов погоды ERA-Interim. Исходные данные представлены для трех слоев почвы (0–7, 7–28, 28–100 см) с дискретностью по времени в 1 сутки и по географическим координатам в 0.75°. Эти данные пересчитывались в среднеобъемную влажность метрового слоя почвы с учетом вклада каждого слоя и линейно интерполировались на территорию РБ. Значения влажности почвы аппроксимировались в зависимости от стандартных метеорологических параметров следующим выражением:
(6)
Таблица 1. Коэффициенты уравнения (6) для областей РБ
Область | с1×10−3 | с2×10−3 | с3×10−3 |
Брестская | 4.162 | 2.398 | 0.841 |
Витебская | 4.559 | 2.214 | 1.049 |
Гомельская | 4.301 | 2.352 | 0.877 |
Гродненская | 4.298 | 2.321 | 0.918 |
Минская | 4.455 | 2.338 | 0.931 |
Могилевская | 4.173 | 2.305 | 0.866 |
где i — порядковый номер дня в году, VPD — дефицит давления водяного пара в гПа, P — суммарное количество осадков за сутки в мм, с1, с2, с3 — положительные определенные коэффициенты, зависящие от почвенно-климатических особенностей региона. Первое слагаемое в (6) — это запасы влаги в почве на начало i-го расчетного периода (суток), второе — учитывает уменьшение этих запасов за счет испарения, а третье и четвертое — их увеличение за счет осадков и конденсации водяного пара соответственно.
Средние значения коэффициентов уравнения (6) для областей РБ приведены в табл. 1. При практическом использовании этого уравнения необходимо задание влажности почвы в начале вегетационного периода, которое можно взять из данных реанализа, получаемых с задержкой около месяца от соответствующей им даты. Расчет последующих значений W выполняется ежесуточно в зависимости от количества выпавших за эти сутки осадков и дефицита влажности воздуха. Пример моделирования динамики W в 2017 г. для трех областей РБ приведен на рис. 3. Если исходить из данных реанализа ERA-Interim за 1979–2017 гг., то абсолютная погрешность расчетов влажности метрового слоя почвы предложенным методом не превышает 0.026 (или 2.6% от объема).
Еще одним метеорологическим параметром, представляющим значительные сложности для прямых наземных измерений, является поток полной солнечной радиации (прямой и рассеянной) на нижней границе атмосферы FBOA (Bottom of Atmosphere). Для его расчета нами используется модифицированная эмпирическая формула из (Supit, Van Kappel, 1998), позволяющая косвенно определять средний за сутки поток солнечной радиации по разнице между максимальной и минимальной температурой приземного воздуха и степени покрытия атмосферы облаками:
, (7)
где FTOA — средний поток солнечной радиации, поступающий на верхнюю границу атмосферу за день (Top of Atmosphere), td — время в днях с начала года, Tmin и Tmax — минимальная и максимальная дневная температура воздуха, CF — степень покрытия атмосферы облаками в дневное время (Cloud Fraction), A и B — эмпирические константы, зависящие от географических координат.
Поток FTOA рассчитывается в зависимости от времени года и зенитного угла Солнца (Sen, 2008):
, (8)
где F0 = 1367 Вт/м2 — солнечная постоянная, — средний косинус зенитного угла Солнца в данное время года.
Среднее значение в данное время года на данной широте β находится по следующей формуле (Sen, 2008):
Рис. 3. Моделирование динамики влажности почвы в Витебской (1), Гродненской (2) и Брестской (3) областях Республики Беларусь: прямые линии — реанализ ERA-Interim, пунктир — регрессионная модель (уравнение 6).
(9)
где — угол склонения Солнца, — часовой угол Солнца (между плоскостями небесного меридиана и круга склонений) при восходе или закате.
Коэффициенты A и B в формуле (7) определялись по данным реанализа ERA-Interim. Средние значения этих коэффициентов для областей РБ приведены в табл. 2. Средние за день значения FBOA, рассчитанные по формулам (7)–(9) для Минской области, сопоставлены с исходными данными реанализа на рис. 4. Коэффициент корреляции между ними ~0.95.
Таблица 2. Значения коэффициентов эмпирической формулы (7) для областей РБ
Область | A | B |
Брестская | 0.162 | 0.522 |
Витебская | 0.168 | 0.563 |
Гомельская | 0.169 | 0.506 |
Гродненская | 0.166 | 0.532 |
Минская | 0.166 | 0.539 |
Могилевская | 0.168 | 0.531 |
Рис. 4. Результаты расчетов среднего за день потока суммарной солнечной радиации на нижней границе атмосферы в сопоставлении с данными реанализа ERA-Interim для Минской области РБ.
Из всего спектра солнечной радиации растениями для фотосинтеза используется только та его часть, которая лежит в полосе поглощения хлорофилла (400–700 нм). Доля солнечной радиации, поглощаемая растительным покровом (FAPAR), зависит от листового индекса, углового распределения листвы, концентраций пигментов листа, оптических свойств почвы, зенитного угла Солнца и других факторов. Полный их учет при определении FAPAR возможен только с привлечением трудоемких полевых измерений. Тем не менее, сезонный ход FAPAR вполне адекватно отражается индексом NDVI, достаточно легко определяемым по измерениям из космоса. Наличие линейной корреляции между FAPAR и NDVI показано как экспериментально, так и на основе моделирования переноса излучения в растительных покровах (Myneni, Williams, 1994; Ganguly et al., 2014). В связи с этим при моделировании динамики биомассы на основе уравнения Монтейта нами используется предположение, что энергия солнечного излучения, поглощаемая за день растительностью, пропорциональна произведению потока FBOA и NDVI.
С учетом сделанных предположений, оценка урожайности СХК на i-е сутки вегетационного периода выполняется на основе следующего уравнения:
(10)
где γ — эмпирический коэффициент, зависящий от региона и СХК.
В общей сложности динамико-статистическая модель урожайности, основанная на формулах (3)–(10), включает шесть неизвестных параметров: a и TОРТ в формуле (3), WОРТ в (4), R0 и TR в (5), γ в (10). Их конкретные значения устанавливаются отдельно для каждого региона и выращиваемой в нем СХК путем сравнения результатов моделирования биомассы с фактической урожайностью в этом регионе.
СПУТНИКОВЫЙ МОНИТОРИНГ ПРОЦЕССА ФОРМИРОВАНИЯ УРОЖАЯ СЕЛЬСКОХОЗЯЙСТВЕННЫХ КУЛЬТУР
Описанный выше метод моделирования динамики растительной биомассы протестирован на примере мониторинга процесса формирования урожая яровых культур, выращиваемых в РБ. Расчеты динамики биомассы СХК выполнялись с разрешением по времени в одни сутки на основе метеорологических данных ERA-Interim (температура, осадки, дефицит влажности воздуха) и данных спутникового спектрорадиометра MODIS. Композиционные индексы NDVI, определяемые по данным прибора MODIS, полагались неизменным на протяжении всего охватываемого ими 16-дневного периода.
Моделирование динамики биомассы начинается со дня получения первого композиционного изображения NDVI, приходящегося на период вегетации (для рассматриваемого региона это 7 апреля). Дата окончания расчетного периода выбрана исходя из анализа корреляции между рассчитываемой биомассой и официальными данными по урожайности культур. Наилучшие корреляции достигаются при завершении расчетного периода в начале июля, т. е. как минимум за месяц до начала уборки.
Параметры динамико-статистической модели биомассы определялись отдельно для каждой области РБ на основе обучающих данных за 18 лет (2000–2017 гг.). Их средние значения для яровых культур приведены в табл. 3. Примеры моделирования динамики биомассы яровой пшеницы представлены на рис. 5.
Таблица 3. Средние значения параметров динамико-статистической модели биомассы яровых культур РБ
Культура | a | TОРТ | WОРТ | R0 | TR | γ |
Пшеница | 0.965 | 293 | 0.261 | 0.037 | 277 | 7.522 |
Ячмень | 0.593 | 291 | 0.247 | 0.026 | 282 | 3.390 |
Овес | 0.791 | 291 | 0.226 | 0.028 | 283 | 3.377 |
Рис. 5. Расчетная динамика биомассы яровой пшеницы в областях Республики Беларусь с 22 марта по 11 июля 2017 г.
Рис. 6. Ряды урожайности яровой пшеницы в Республике Беларусь за 2000–2017 гг. соответствующие данным официальной статистики (1) и прогнозам модели, откалиброванной на обучающих данных за все годы (2) и за все годы кроме текущего года на графике (3).
Результаты ретроспективного моделирования средней урожайности яровой пшеницы в Республике Беларусь приведены на рис. 6. Коэффициент корреляции между рассчитанными значениями урожайности и данными официальной статистики составляет 0.84. В виду небольшого объема выборки для проверки качества модели использована процедура кросс-валидации Leave-One-Out, основанная на последовательном исключении из обучающей выборки данных для одного года, калибровке модели по данным для оставшихся годов и тестировании модели на исключаемых данных. Поскольку временной ряд урожайности, приведенный на рис. 6, предварительно освобожден от трендовой компоненты, то каждая точка этого ряда является полностью независимой от остальных. Это дает возможность провести валидацию созданной модели без получения новых данных по урожайности СХК.
По результатам кросс-валидации коэффициент корреляции между фактической и предсказанной урожайностью яровой пшеницы составил ~0.70, что говорит о достаточной устойчивости модели к изменчивости метеорологических условий формирования урожая. Отклонения между фактическими и расчетными значениями урожайности могут быть связаны с недостаточно точным описанием суточных метеорологических полей данными реанализа. В связи с этим, дальнейшее повышение точности модели возможно за счет использования для ее калибровки данных прямых метеорологических измерений на сети наземных станций. Однако при этом возникают проблемы пропусков метеорологических данных и неоднородности покрытия ими территории, устранение которых представляет отдельную и достаточно непростую задачу.
ЗАКЛЮЧЕНИЕ
В работе представлена малопараметрическая модель динамики биомассы сельскохозяйственных культур, учитывающая основные метеорологические факторы формирования урожая и легко адаптируемая для различных регионов и сельскохозяйственных культур. Для моделирования динамики биомассы требуются данные спутниковых измерений вегетационного индекса подстилающей поверхности и стандартных метеорологических наблюдений (температура, влажность воздуха, осадки и облачность). Труднодос- тупные для измерений метеорологические параметры, такие как поток солнечного излучения на нижней границе атмосферы и влажность почвы, рассчитываются на основе эмпирических формул с коэффициентами, получаемыми для исследуемого региона по данным реанализа.
В расчетах динамики биомассы используются только те территории региона, для которых вегетационный индекс коррелирует с урожайностью рассматриваемой сельскохозяйственной культуры, что устраняет необходимость идентификации пахотных земель на спутниковых снимках. Моделирование динамики биомассы сразу на уровне административных единиц государства также минимизирует объем информации, необходимой для калибровки модели.
На примере отдельного государства (Республика Беларусь) показана применимость разработанной модели для прогнозирования урожайности сельскохозяйственных культур с использованием легкодоступных данных официальной статистки, реанализа Era-Interim и спутниковых измерений NDVI прибором MODIS.
About the authors
S. A. Lysenko
Institute for Nature Management of National Academy of Science of Belarus
Author for correspondence.
Email: lysenkorfe@gmail.com
Russian Federation, Minsk
References
- Atzberger C. Advances in Remote Sensing of Agriculture: Context Description, Existing Operational Monitoring Systems and Major Information Needs // Remote Sensing. 2013. V. 5. № 2. P. 949–981. doi: 10.3390/rs5020949
- Baret F., Buis S. Estimating Canopy Characteristics from Remote Sensing Observations: Review of Methods and Associated Problems / Liang S. (eds). Advances in Land Remote Sensing. Springer: Dordrecht, 2008. P. 173–201. doi: 10.1007/978–1–4020–6450–0_7
- Becker-Reshef I., Vermote E., Lindeman M., Justice C. A generalized regression-based model for forecasting winter wheat yields in Kansas and Ukraine using MODIS data // Remote Sens. Environ. 2010. V. 114. № 6. P. 1312–1323. doi: 10.1016/j.rse.2010.01.010
- Bryskin V.M., Evtyushkin A.V. Ispol'zovanie modeli bio-produktivnosti EPIC i kosmosnimkov MODIS dlya progno-zirovaniya urozhainosti zernovykh kul'tur // Sovremennye problemy distantsionnogo zondirovaniya Zemli iz kosmosa: fizicheskie osnovy, metody i tekhnologii monitoringa okruzhayushchei sredy, potentsial'no opasnykh yavlenii i ob"ektov: sb. nauch. st. [Using the model of bioproductivity of EPICand MODIS satellite imagery to predict the yield of grain crops // Modern problems of remote sensing of the Earth from space: physical bases, methods and technologies for monitoring the environment, potentially dangerous phenomena and objects] M.: 2007. T. II. № 4. P. 189–196. (In Russian).
- Delécolle R., Maas S.J., Guérif M., Baret F. Remote sensing and crop production models: Present trends // ISPRS J. Photogramm. 1992. V. 47. № 2–3. P. 145–161. doi: 10.1016/0924–2716(92)90030-D
- Dmitrenko V.P. Pogoda, klimat i urozhai polevykh kul'tur [Weather, climate and crop harvest]. Kiev: Nika-tsentr, 2010. 620 p. (In Russian).
- Duchemin B., Maisongrande P., Boulet G., Benhadj I. A simple algorithm for yield estimates: Evaluation for semi-arid irrigated winter wheat monitored with green leaf area index // Environ. Modell. Softw. 2008. V. 23. № 7. P. 876–892. doi: 10.1016/j.envsoft.2007.10.003
- Ferencz Cs., Bognár P., Lichtenberger J., Hamar D., Tarcsai Gy., Timár G., Molnár G., Pásztor Sz., Steinbach P., Székely B., Ferencz O.E., Ferencz-Árkos I. Crop yield estimation by satellite remote sensing // Int. J. Remote Sens. 2004. V. 25. № 20. P. 4113–4149. doi: 10.1080/01431160410001698870
- Ganguly S., Nemani R.R., Baret F., Bi J., Weiss M., Zhang G., Milesi C., Hashimoto H., Samanta A., Verger A., Singh K., Myneni R.B. Green Leaf Area and Fraction of Photosynthetically Active Radiation Absorbed by Vegetation / J.M. Hanes (ed.). Biophysical Applications of Satellite Remote Sensing. Berlin: Springer-Verlag Berlin Heidelberg, 2014. P. 43–61. doi: 10.1007/978–3–642–25047–7_2
- Hansen M.C., Defries R.S., Townshend J.R.G., Sohlberg R. Global land cover classification at 1 km spatial resolution using a classification tree approach // Int. J. Remote Sens. 2000. V. 21. № 6 & 7. P. 1331–1364. doi: 10.1080/014311600210209
- Huete A., Miura T., Yoshioka H., Ratana P., Broich M. Indices of Vegetation Activity / J.M. Hanes (ed.). Biophysical Applications of Satellite Remote Sensing. Berlin: Springer-Verlag Berlin Heidelberg, 2014. P. 1–41. doi: 10.1007/978–3–642–25047–7_1
- Jin X., Kumar L., Li Z., Xu X., Yang G., Wang J. Estimation of Winter Wheat Biomass and Yield by Combining the AquaCrop Model and Field Hyperspectral Data // Remote Sens. 2016. V. 8. № 12. P. 972–1–972–15. doi: 10.3390/rs8120972
- Kasampalis D.A., Alexandridis T.K., Deva C., Challinor A., Moshou D., Zalidis G. Contribution of Remote Sensing on Crop Models: A Review // J. Imaging. 2018. V. 4. № 4. P. 52–1–52–19. doi: 10.3390/jimaging4040052
- Kastens J.H., Kastens T.L., Kastens D.L.A., Price K.P., Martinko E.A., Lee R.-Y. Image masking for crop yield forecasting using AVHRR NDVI time series imagery // Remote Sens. Environ. 2005. V. 99. № 3. P. 341–356. doi: 10.1016/j.rse.2005.09.010
- Kleshchenko A.D., Lebedeva V.M., Naidina T.A., Savits-kaya O.V. Ispol'zovanie sputnikovoi informatsii MODIS v operativnoi agrometeorologii [The use of MODIS satellite information in operational agrometeorology] // Sovremennye problemy DZZ iz kosmosa. 2015. V. 12. № 2. P. 143–154. (In Russian).
- Kleshchenko A.D., Naidina T.A. Ispol'zovanie dannykh distantsionnogo zondirovaniya dlya modelirovaniya fiziologicheskikh protsessov rastenii v dinamicheskikh modelyakh prognozirovaniya urozhaya [Using remote sensing data to model plant physiological processes in dynamic yield forecast models] // Sovremennye problemy distantsionnogo zondirovaniya Zemli iz kosmosa. 2011. V. 8. № 1. P. 170–177. (In Russian).
- Kouadio L., Newlands N.K., Davidson A., Zhang Y., Chipanshi A. Assessing the Performance of MODIS NDVI and EVI for Seasonal Crop Yield Forecasting at the Ecodistrict Scale // Remote Sens. 2014. V. 6. № 10. P. 10193–10214. doi: 10.3390/rs61010193
- Kussul' N.N., Kravchenko A.N., Skakun S.V., Adamenko T.I., Shelestov A. Yu., Kolotii A.V., Gripich Yu.A. Regressionnye modeli otsenki urozhainosti sel'skokhozyaistvennykh kul'tur po dannym MODIS [Regression models for crop yield estimation according to MODIS] // Sovremennye problemy distantsionnogo zondirovaniya Zemli iz kosmosa. 2012. V. 9. № 1. P. 95–107. (In Russian).
- Liu J., Pattey E., Miller J.R., McNairn H., Smith A., Hu B. Estimating crop stresses, aboveground dry biomass and yield of corn using multi-temporal optical data combined with a radiation use efficiency model // Remote Sens. Environ. 2010. V. 114. № 6. P. 1167–1177. doi: 10.1016/j.rse.2010.01.004
- Monteith J.L. Solar radiation and productivity in tropical ecosystems // J. Appl. Ecol. 1972. V. 9. № 3. P. 747–766. doi: 10.2307/2401901
- Myneni R.B., Williams D.L. On the Relationship between FAPAR and NDVI // Remote Sens. Environ. 1994. V. 49. № 3. P. 200–211. doi: 10.1016/0034–4257(94)90016–7
- Plotnikov D.E., Bartalev S.A., Zharko V.O., Mikhailov V.V., Prosyannikova O.I. Ehksperimental'naya otsenka raspoz-navaemosti agrokul'tur po dannym sezonnykh sputnikovykh izmerenii spektral'noi yarkosti [Experimental assessment of the recognition of crops from seasonal satellite mea-surements of spectral brightness] // Sovremennye problemy distantsionnogo zondirovaniya Zemli iz kosmosa. 2011. V. 8. № 1. P. 199–208. (In Russian).
- Poluektov R.A., Smolyar E.I., Terleev V.V., Topazh A.G. Mo-deli produktsionnogo protsessa sel'skokhozyaistvennykh kul'-tur[Crop production process models]. SPb.: Izd-vo Sankt-Peterburgskogo universiteta, 2006. 392 p. (In Russian).
- Rembold F., Atzberger C., Savin I., Rojas O. Using Low Resolution Satellite Imagery for Yield Prediction and Yield Anomaly Detection // Remote Sens. 2013. V. 5. № 4. P. 1704–1733. doi: 10.3390/rs5041704
- Savin I. Yu., Bartalev S.A., Lupyan E.A., Tolpin V.A., Khvostikov S.A. Prognozirovanie urozhainosti sel'skokhozyaistvennykh kul'tur na osnove sputnikovykh dannykh: vozmozhnosti i perspektivy [Predicting crop yields based on satellite data: opportunities and prospects] // Sovremennye problemy distantsionnogo zondirovaniya Zemli iz kosmosa. 2010. V. 7. № 3. P. 275–285. (In Russian).
- Sen Z. Solar Energy Fundamentals and Modeling Techniques Atmosphere, Environment, Climate Change and Renewable Energy. London: Springer-Verlag London Limited, 2008. 276 p.
- Sirotenko O.D. Matematicheskoe modelirovanie vodno-teplovogo rezhima i produktivnosti agroehkosistem [Mathematical modeling of water-thermal regime and productivity of agroecosystems]. L.: Gidrometeoizdat, 1981. 168 p. (In Russian).
- Sladkikh L.A., Zakhvatov M.G., Saprykin E.I., Sakha-rova E. Yu. Tekhnologiya monitoringa sostoyaniya posevov po dannym distantsionnogo zondirovaniya Zemli na yuge Zapadnoi Sibiri [Technology of monitoring the state of crops according to remote sensing data in the south of Western Siberia] // Geomatika. 2016. № 2. P. 39–48. (In Russian).
- Strashnaya A.I., Bartalev S.A., Maksimenkova T.A., Chub O.V., Tolpin V.A., Plotnikov D.E., Bogomolova N.A. Agro-meteorologicheskaya otsenka sostoyaniya ozimykh zernovykh kul'tur v period prekrashcheniya vegetatsii s ispol'zo-vaniem nazemnykh i sputnikovykh dannykh na primere Privolzhskogo federal'nogo okruga [Agrometeorological assessment of the state of winter crops in the period of vegeta-tion termination using ground and satellite data on the example of the Volga Federal District] // Trudy gidrometeorolo-gicheskogo nauchno-issledovatel'skogo tsentra Rossiiskoi Federatsii. 2014. № 351. P. 85–107. (In Russian).
- Supit I., Van Kappel R.R. A simple method to estimate global radiation // Solar Energy. 1998. V. 63. № 3. P. 147–160. doi: 10.1016/S0038–092X(98)00068–1
- Terekhov A.G., Vitkovskaya I.S., Batyrbaeva M. Zh., Spivak L.F. Printsipy agrolandshaftnogo raionirovaniya pakhotnykh zemel' Severnogo Kazakhstana po dannym LANDSAT i MODIS [Principles of agrolandscape zoning of arable lands of Northern Kazakhstan according to LANDSAT and MODIS data] // Sovremennye problemy distantsionnogo zondirovaniya Zemli iz kosmosa. 2010. V. 7. № 3. P. 292–304. (In Russian).
- Xiao X., Jin C., Dong J. Gross Primary Production of Terrestri- al Vegetation / J.M. Hanes (ed.). Biophysical Applications of Satellite Remote Sensing. Berlin: Springer-Verlag Berlin Heidelberg, 2014. P. 127–148. doi: 10.1007/978–3–642–25047–7_5
- Zhang P., Anderson B., Tan B., Huang D., Myneni R. Potential monitoring of crop production using a satellite-based Climate- Variability Impact Index // Agr. Forest Meteorol. 2005. V. 132. № 3–4. P. 344–358. doi: 10.1016/j.agrformet.2005.09.004
Supplementary files
![](/img/style/loading.gif)