Crop yield forecasting based on the satellite monitoring of carbon dynamics in terrestrial ecosystems

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):

$\begin{array}{c}NP{P}_i=PA{R}_i\times FAPA{R}_i\times \\ \times LUE\left({\text{CO}}_2,T_i,VP{D}_i,W_i,\mathrm{...}\right)\times \left(1-R_i\right), \end{array}$ (1)

где i — номер расчетного периода, NPPi (Net Primary Production) — чистая первичная продукция органического вещества (гр. сухой биомассы/м²/период), PARi (Photosynthetically Active Radiation) — энергия солнечного излучения, поступающего за заданный отрезок времени на подстилающую поверхность в фотосинтетически активной области спектра, FAPARi (Fraction of Absorbed Photosynthetically Active Radiation) — доля фотоситнетически активной радиации, поглощаемая пигментами зеленого листа, LUE (Light Use Efficiency) — эффективность преобразования поглощенной растительностью фотоситнетически активной радиации в органическое вещество, зависящая от ряда метеорологических параметров, CO₂ — концентрация углекислого газа в атмосфере, которая на протяжение вегетационного периода практически не изменяется, T_i — температура приземного воздуха, VPD_i (Vapour Pressure Deficit) — дефицит давления водяного пара (разница между давлением насыщенного пара и фактическим парциальным давлением пара при данной температуре), W_i — содержание влаги в деятельном слое почвы, R_i — доля затрат органического вещества на дыхание живой биомассы.

При использовании уравнения (1) для расчетов урожаев СХК делается предположение, что фитомасса урожая Y пропорциональна общей биомассе, накопленной с начала вегетационного периода, т. е.:

$Y=C\sum _{i=1}^{N}NP{P}_i,$ (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 зададим в виде произведения двух функций η₁(T) и η₂(W), определяющих зависимость интенсивности фотосинтеза от температуры воздуха T и влажности почвы W.

Функция η₁(T) имеет вид одновершинной кривой с максимальным значением η₁(T_opt) = 1 в температурном оптимуме фотосинтеза T_opt, зависящем от типа растительности. Известно несколько функциональных зависимостей η₁(T), пригодных для разных типов растительности. Однако, для того чтобы не перегружать модель биопродуктивности слишком большим количеством свободных параметров, которые подлежат установлению на основе ограниченного временного ряда эмпирических данных, будем использовать самую простую температурную зависимость фотосинтеза, учитывающую ее наиболее характерные особенности (Сиротенко, 1981):

${\eta }_1\left(T\right)=\exp \left\{-a{\left(\left(T-T_{opt}\right)/10\right)}^2\right\}$, (3)

где a — положительный безразмерный параметр.

Для функции η₂(W) также выбрано наиболее простое выражение с единственным свободным параметром (Дмитренко, 2010):

${\eta }_2\left(W\right)=1-{\left[\left(W-W_{opt}\right)/W_{opt}\right]}^2$, (4)

где W и W_opt — фактическая и оптимальная влажность метрового слоя почвы.

Затраты биомассы на дыхание рассчитываются в зависимости от температуры воздуха по следующей формуле (Сиротенко, 1981):

$R\left(T\right)=R_0\cdot 2^{\left(T-T_R\right)/10}$ (5)

где TR — базисная температура дыхания, превышение которой температурой воздуха на 10 °C сопровождается увеличением расходов органического вещества на дыхание в два раза.

В связи со сложностью и высокой стоимостью полевых измерений влажности почвы на больших территориях, данные по ней, как правило, относятся лишь к отдельным опытно-экспериментальным участкам и недоступны в оперативном режиме. Тем не менее, ее можно рассчитать на основе стандартных метеорологических параметров (таких как температура, осадки, влажность воздуха и скорость ветра) с использованием моделей переноса влаги в системе «почва–растение–атмосфера» (Сиротенко, 1981; Полуэктов и др., 2006). Однако при мониторинге состояния и развития СХК из космоса такое моделирование необходимо проводить отдельно для каждого пиксела спутникового снимка, что приводит к неоправданно большим вычислительным затратам. Выходом из ситуации может быть использование аналитических параметризаций зависимости влажности почвы от стандартных метеорологических параметров. Для получения таких параметризаций необходимы репрезентативные данные натурных наблюдений или расчетные данные моделей переноса влаги.

Нами для расчета влажности метрового слоя почвы использовались аналитические параметризации, полученные на основе данных реанализа Европейского центра среднесрочных прогнозов погоды ERA-Interim. Исходные данные представлены для трех слоев почвы (0–7, 7–28, 28–100 см) с дискретностью по времени в 1 сутки и по географическим координатам в 0.75°. Эти данные пересчитывались в среднеобъемную влажность метрового слоя почвы с учетом вклада каждого слоя и линейно интерполировались на территорию РБ. Значения влажности почвы аппроксимировались в зависимости от стандартных метеорологических параметров следующим выражением:

$W_i=W_{i-1}-c_1\cdot VP{D}_i\cdot W_{i-1}+c_2\cdot P_i+c_3\cdot VP{D}_i$ (6)

Таблица 1. Коэффициенты уравнения (6) для областей РБ

где 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), позволяющая косвенно определять средний за сутки поток солнечной радиации по разнице между максимальной и минимальной температурой приземного воздуха и степени покрытия атмосферы облаками:

$F_{BOA}=F_{TOA}\left(t_d\right)\left[A\sqrt{T_{\max }-T_{\min }}+B\sqrt{1-CF}\right]$, (7)

где F_TOA — средний поток солнечной радиации, поступающий на верхнюю границу атмосферу за день (Top of Atmosphere), t_d — время в днях с начала года, T_min и T_max — минимальная и максимальная дневная температура воздуха, CF — степень покрытия атмосферы облаками в дневное время (Cloud Fraction), A и B — эмпирические константы, зависящие от географических координат.

Поток F_TOA рассчитывается в зависимости от времени года и зенитного угла Солнца (Sen, 2008):

, (8)

где F₀ = 1367 Вт/м² — солнечная постоянная,  — средний косинус зенитного угла Солнца в данное время года.

Среднее значение  в данное время года на данной широте β находится по следующей формуле (Sen, 2008):

Рис. 3. Моделирование динамики влажности почвы в Витебской (1), Гродненской (2) и Брестской (3) областях Республики Беларусь: прямые линии — реанализ ERA-Interim, пунктир — регрессионная модель (уравнение 6).

$\overline{\cos {\theta }_s}=\sin \beta \sin \delta +\cos \beta \cos \delta \frac{\sin {\psi }^{*}}{{\psi }^{*}}$ (9)

где $\delta =-23.45\cos \left(2\pi \left(t_d+10\right)/365\right)$ — угол склонения Солнца, ${\psi }^{*}=\text{arccos}\left(-\text{tg}\beta \text{tg}\delta \right)$ — часовой угол Солнца (между плоскостями небесного меридиана и круга склонений) при восходе или закате.

Коэффициенты A и B в формуле (7) определялись по данным реанализа ERA-Interim. Средние значения этих коэффициентов для областей РБ приведены в табл. 2. Средние за день значения FBOA, рассчитанные по формулам (7)–(9) для Минской области, сопоставлены с исходными данными реанализа на рис. 4. Коэффициент корреляции между ними ~0.95.

Таблица 2. Значения коэффициентов эмпирической формулы (7) для областей РБ

Рис. 4. Результаты расчетов среднего за день потока суммарной солнечной радиации на нижней границе атмосферы в сопоставлении с данными реанализа ERA-Interim  для Минской области РБ.

Из всего спектра солнечной радиации растениями для фотосинтеза используется только та его часть, которая лежит в полосе поглощения хлорофилла (400–700 нм). Доля солнечной радиации, поглощаемая растительным покровом (FAPAR), зависит от листового индекса, углового распределения листвы, концентраций пигментов листа, оптических свойств почвы, зенитного угла Солнца и других факторов. Полный их учет при определении FAPAR возможен только с привлечением трудоемких полевых измерений. Тем не менее, сезонный ход FAPAR вполне адекватно отражается индексом NDVI, достаточно легко определяемым по измерениям из космоса. Наличие линейной корреляции между FAPAR и NDVI показано как экспериментально, так и на основе моделирования переноса излучения в растительных покровах (Myneni, Williams, 1994; Ganguly et al., 2014). В связи с этим при моделировании динамики биомассы на основе уравнения Монтейта нами используется предположение, что энергия солнечного излучения, поглощаемая за день растительностью, пропорциональна произведению потока FBOA и NDVI.

С учетом сделанных предположений, оценка урожайности СХК на i-е сутки вегетационного периода выполняется на основе следующего уравнения:

$Y=\gamma \sum _{i=1}^N{F}_{TOA,i}\times NDV{I}_i\times {\eta }_1\left(T_i\right)\times {\eta }_2\left(W_i\right)\times \left(1-R\left(T_i\right)\right),$ (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. Средние значения параметров динамико-статистической модели биомассы яровых культур РБ

Рис. 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

Author for correspondence.
Email: lysenkorfe@gmail.com
Russian Federation, Minsk

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1. The dynamics of the yield of grain and leguminous crops cultures in the regions of the Republic of Belarus from 1990 to 2017 years (dotted line with symbols) and linear trends Jail before and after 2000 (direct).

Download (243KB)

Indexing metadata

3. Fig. 2. Map of the correlation coefficient between the crop spring wheat in the regions of the Republic of Belarus Russia and 16-day vegetation indices of thorium within these areas.

Download (583KB)

Indexing metadata

4. Fig. 3. Modeling the dynamics of soil moisture in the Vitebsk (1), Grodno (2) and Brest (3) regions of the Res- Belarus public: straight lines - reanalysis of ERA-Interim, dashed line - regression model (equation 6).

Download (176KB)

Indexing metadata

5. Fig. 4. The results of calculations of the average daily flow of sum solar radiation at the lower boundary of the atmosphere phera versus ERA-Interim reanalysis data for the Minsk region of Belarus.

Download (202KB)

Indexing metadata

6. Fig. 5. The estimated dynamics of the biomass of spring wheat in the regions of the Republic of Belarus from March 22 to July 11 2017 year

Download (126KB)

Indexing metadata

7. Fig. 6. Series of spring wheat productivity in the Republic Belarus 2000–2017 corresponding to the data social statistics (1) and model forecasts, calibrated bathtub on training data for all years (2) and for all years except for the current year on the chart (3).

Download (111KB)

Indexing metadata