CC BY
65
УДК 004.94
Л.А. Хвороба, В.М. Брыксин, Н.В. Гаврилобская
ИСПОЛЬЗОВАНИЕ ИНФОРМАЦИОННЫХ ТЕХНОЛОГИЙ ПРИ ПРОГНОЗИРОВАНИИ УРОЖАЙНОСТИ ЗЕРНОВЫХ КУЛЬТУР
Работоспособные технологии прогноза урожайности зерновых культур, интегрированные в сельскохозяйственное производство и основанные на математических моделях и космической информации, в условиях России отсутствуют. В связи с этим решение проблемы осуществления упреждающего прогноза по адекватной математической модели с использованием данных космосъемки весьма актуально.
Особое место в технологии прогноза урожайности зерновых культур должны занимать прогнозы, позволяющие давать оценку ожидаемой продуктивности посевов поэтапно. Это значительно увеличивает надежность прогнозов урожайности и эффективность использования прогностической информации. Методика поэтапного прогнозирования создавалась авторами в рамках разработки технологии определения года-аналога [4, 5].
Применение адаптированной модели бионродуктивности EPIC и космоснимков MOD IS
для прогнозирования урожайности зерновых культур на территории Западной Сибири
По мнению многих ученых, основу методического и информационного обеспечения работ по оценке урожайности зерновых культур должны составлять динамические имитационные модели биопродуктивности сельскохозяйственных кулыур, комплексы математических моделей прогнозирования урожайности с использованием данных дистанционного зондирования Земли.
Авторами разработана технология прогноза урожайности зерновых культуре использованием космических снимков и компьютерного моделирования, интеллектуальным ядром которого является модифицированная динамическая модель EPIC [6], разработанная Сельскохозяйственной научно-исследовательской службой департамента сельского хозяйства Соединенных Штатов (Soil and Water Research Laboratory US DA) [I, 3j.
Технология прогноза урожайности зерновых культур реализована в виде комплекса программных средств MIDC (модель имитации развития культур) и апробирована в условиях зерновых провинций Западной Сибири — Новосибирской и Тюменской областей. Поданной модели ежедневно рассчитываются основные параметры культуры (биомасса, листовой индекс, высота растения и др.). Это позволяет сравнивать их с данными наземных измерений и дистанционного зондирования и в случае необходимости производить корректировку. Используя фактические данные по сложившимся агрометеорологическим условиям, программный комплекс позволяет осуществлять прогноз урожайности за месяц до уборки.
Необходимо отметить, что модель EPIC разрабатывалась для оценки последствий эрозии почв на продуктивность посевов, что, естественно, наложило отпечаток на ее структурную и алгоритмическую реализацию. Поэтому на первый план вышли проблемы ее идентификации по имеющейся агрометеорологической информации и реализации цели исследования.
Этапы структурно-параметрической идентификации модели EPIC
Процедура идентификации модели EPIC и разработка на ее основе имитационно-модели-рующего комплекса MIDC включала шесть этапов:
1) тестирование модели EPIC (выявлялась реакция по блокам и модели в целом на изменение основных агрометеорологических факторов);
2) модификацию блоков и структуры модели EPIC по результатам тестирования и на основе имеющейся агрометеорологической и спутниковой информации;
3) конструирование на базе технологии EPIC имитационно-молелируюшего комплекса MIDC;
4) определение параметров, входящих в блоки: коротковолновой (КВР), длинноволновой (ДВР) радиапий, теплового режима почвы и функции суточной динамики температуры воздуха;
5) комплексную идентификацию модели по экспериментальным данным регионов;
6) апробацию методики прогноза урожайности в реальных условиях Западной Сибири.
Рассмотрим результаты отдельных этапов идентификации.
Тестирование модели EPIC проводилось с использованием данных, предоставленных Западно-Сибирским управлением гидрометслуж-бы по Новосибирской области и Тюменским центром гидрометслужбы по Тюменской области. Выяснилось, что предоставленные наборы arpo- и метеоданных не полные, отсутствуют необходимые для модели величина приходящей КВР и температура поверхности почвы. В связи с этим возникла задача восполнения неполноты исходной arpo- и метеоинформации за счет разработки блоков радиации (КВР и ДВР).
При отсутствии экспериментальных данных по приходящей коротковолновой радиации в EPIC генерируются данные, сильно отличающиеся от реальных. В результате возникают большие погрешности в расчетах [3]. Как показал опыт сбора необходимой агрометеоимформации, многие метеостанции коротковолновую радиацию не измеряют. Поэтому возникает проблема модификации и настройки модели КВР для регионов Западной Сибири.
Важный входной параметр модели - температура поверхности почвы (минимальная и максимальная). Как правило, такая информация также отсутствует на многих метеостанциях. В связи с этим возникла задача расчета теплового режима почвы (определение минимальной и максимальной температуры на поверхности почвы и расчет ее температурного профиля).
Выявленные недостатки при тестировании модели EPIC потребовали разработки дополнительного программного модуля, состоящего из трех взаимосвязанных блоков: блока ДВР, или эффективного излучения (для определения температуры на поверхности почвы в ночное время), блока КВР (для определения температуры на поверхности почвы в дневные часы) и блока теплового режима почвы. Блочная структура модели
ЕР1Сделаетее гибкой и дает возможность включать новые блоки (рис. 1).
Определение температуры на поверхности почвы. Как известно, днём поверхность почвы нагревается прямой солнечной радиацией, а ночью — охлаждается в результате эффективного излучения, и от неё охлаждаются прилегающие слои воздуха [2]. Причем минимальная температура воздуха наблюдается перед восходом Солнца, а максимальная температура отмечается через два-три часа после полудня.
По экспериментальным данным установлена зависимость между максимальной температурой почвы и суммарной радиацией:
т =13.7* +0 5-0+05
'ртах 1'J * vmax U>J «i->»
где Tp max — максимальная температура на поверхности почвы; Tvmax — максимальная температура воздуха; Q — суммарная коротковолновая радиация.
В ночное время наблюдается взаимосвязь между минимальной температурой почвы и эффективным излучением. В качестве аппроксимирующей математической функции было выбрано уравнение линейной двухфакторной регрессии
7;n,in = 0,8-nmjn+1,5-£+4,6,
где Tpmin - минимальная температура на поверхности почвы; 7^т(п — минимальная температура воздуха; Е— эффективное излучение.
Параметры уравнения регрессии оценивались по выборочному, но достаточно большому временному ряду экспериментальных значений температуры на поверхности почвы и расчетным значениям Е из модели ДВР [2] по методу наименьших квадратов. С помощью корреляционного анализа осуществлена проверка адекватности полученной модели.
Температура профиля почвы рассчитывается с использованием уравнения теплопере-носа.
Суточная динамика температуры воздуха.
При расчетах длинноволновой радиации и температурного профиля почвы необходимо использовать данные потемпературе воздуха и верхнего слоя почвы с шагом времени в один час и менее. Однако в предоставляемыхданиых вболь-шинстве случаев присутствует лишь информация с шагом в одни сутки. Для решения данной проблемы предложено воспользоваться кусочно-
Рис. 1. Блок-схема модернизированной модели EPIC (серым цветом выделены измененные блоки)
непрерывной функцией F(t,J), где t — время от полуночи в часах,у — номер дня от начала года, построенной на основе минимальной и максимальной температур воздуха исходя из следующих соображений:
минимальная температура наблюдается на рассвете, а максимальная — через два часа после полудня;
к закату температура падает до значения ТтЩ/+1+ *( Tmaxy- Wi)'гае коэффициент к е [0; 1 ] рассч иты вается по экспери ментал ьн ы м дан и ы м ;
в интервалах [tdJ; tm + 2], [t„ + 2; tsJ и [tS J; tdj+1] функция F(t,J) задается квадратичной зависимостью ar2 + bt+ с, где td, tm, ts - время рассвета, полудня и заката соответственно. При этом F{t,j) имеет следующие точки экстремума для соответствующих интервалов: tm + 2 — максимум для интервала [tdj, tm + 2); tm + 2 ~ максимум для [tm+2, ts j\\ td-минимум для \ ts j, tdJ+l]. Исходя из вышеперечисленных соображений составляем три системы линейных уравнений для поиска коэффициентов а, Ь, с. При этом учитываем, что в точках стыковки функция должна быть непрерывна:
1)
2)
F(td,j) = TmixJ
F(tm + 2,7) = 7max j при / 6 [tdJ- tm + 2]; F'(tm+2,j) = 0
F(tm + 2,j) = TmaxJ
ntsjJ) ~ ^min /+1 + ^ (^rnax j ~ ^min j+1 )
при / e [tm +2; tStJ\\ F'(tm+2J) = 0
3)
F0sjJ) = Tnmj+\ + k(Tt F(td,j+\>J) = Tmin 7+1 при t e [ts j\ îd J + ,];
ntdJ+lJ) = o
max J ^miri j+1 )
При нахождении а, Ь, с учитывалось, что интервал 11 с учетом смены дат разбивается на два интервала: [0; г^] и [г5 у; 24]. Таким образом получаем четыре набора коэффициентов а, Ь, с для соответствующих интервалов.
На рис. 2 приведен пример суточной динамики расчетной и фактической температуры воздуха, °С, для с. Ключи Алтайского края (03.07.2000). Погрешность вычислений составляет 4,7 %.
Корректировка листового индекса и биомассы по спутниковым данным. Значения вегетационного индекса (NDVI), рассчитываемого по космоснимкам сканера MODIS с космических аппаратов TERRA и AQUA, хорошо коррелируют с такой характеристикой растительного покрова, как листовой индекс (LAI). Авторами алгоритма MOD 15 (NASA) для обработки изображений MODIS были получены кусочно-линейные зависимости между LAI и NDVI лля различных ландшафтов. Расчет индекса Ж)И/проводился по контурам тестовых полей.
На основе модельных и фактических данных были получены зависимости, приведённые ниже, биомассы и массы корней от листового индекса, рассчитанные на основе индекса NDVI, что позволяет производить корректировку моделируемых значений:
Я/ОЛ/= 1,09-LAlNpyi- 0,57 при росте LAl\ ВЮМ= В10МШ+Q,26(LAIDM— LAIndv,) при падении LAI\
/?ж=0,31-я/0л/ + 0,03, где ВЮМи RW— моделируемые значения биомассы и массы корня; LAINDVI — листовой индекс, полученный из индекса NDVI\ LAIÜM и BIOMDM— моделируемый листовой индекс и биомасса вдень, когда LAIмаксимален.
Данная модель позволяет производить корректировку в процессе расчетов. Например, корректировку листового индекса LAI можно вы-
полнять по данным дистанционного зондирования через индекс Л7)К/(нормализованный вегетационный индекс), который линейно связан с LAI.
Создание технологии прогноза урожайности зерновых культур по отдельным областям Западной Сибири с использованием космических снимков и компьютерного моделирования включало следующие этапы:
1. Адаптацию модели биопродуктивности EPIC и программных средств для расчета урожайности [1]. Данный этап включает в себя разработку не только базы наземных данных и обработанных данных дистанционного зондирования, но и блок-схем алгоритмов расчетов по прогнозу урожайности.
2. Разработку методики сбора и обработки исходных наземных данных и данных дистанционного зондирования. Этап предполагает сбор необходимой агрометеоинформации для модели в гидрометеослужбах регионов и в результате наземного обследования тестовых полей.
3. Разработку алгоритмов и программных средств поддержки модели, включая проверку адекватности. Этап предполагает усовершенствование модели с включением вновь разработанных блоков: приходящей КВР и ДВР, теплового режима почвы.
4. Апробацию методики прогноза урожайности в реальных условиях Западной Сибири.
В течение вегетационного периода для контуров тестовых полей по космоснимкам вычисляется индекс NDVI, по которому рассчитывается листовой индекс LAIw передается в модель для корректировки расчетов.
1 1 1 • _ |
1 1 1 1 1 "С-Ч.
1 1 . / 1 /
^¡f Ч, г
0:00 4:00 8:00 12:00 16:00 20:00 I
Рис. 2. Суточная динамика расчетной (-) и фактической (--) температуры воздуха
По заказу Департ амента агропромышленного комплекса администрации Тюменской области в 2004 году проводились работы по мониторингу состояния сельскохозяйственных полей и прогнозированию урожайности яровых зерновых культур. Для каждого административного района Тюменской области, имеющего метеостанцию, по модифицированной модели EPIC было проведено моделирование роста и развития зерновых культур (яровая пшеница, ячмень, овес, горох). На этом этапе осуществлен подбор параметров модели для каждого района (определялись усредненные сроки сева и почвенные характеристики модели) и общие для всех районов варьируемые характеристики (индекс урожайности и другие параметры). Определение даты начала уборки проводилось по двум критериям: прекращению наращивания биомассы и падению листового индекса до нуля. Погрешность расчетов урожайности по модели д ля Тюменской области в целом не превышает 1 ц/га.
Для регионов Западной Сибири также проводилась коррекция листового индекса с использованием зависимости LAIот NDVIподанным сканера MODIS. Так, например, для Коченёвского района Новосибирской области в 2005 году урожайность по модели без коррекции листового индекса составила 13,29 ц/га, с коррекцией —15,35 ц/га. Фактическая урожайность 15,7 ц/га.
Динамика прогнозируемой урожайности в зависимости от фактических метеоусловий пред-
ставлена на рис. 3. Видно, что колебания в расчетных значениях биологической урожайности яровой пшеницы минимальны в конце июля. Это подтверждает возможность прогноза урожайности яровой пшеницы на юге Западной Сибири за месяц до начала уборки.
Технология определения года-аналога
Модели любого уровня, осуществляющие прогноз, должны опираться на предполагаемые погодные сценарии прогнозного года. Отсутствие фундаментальной теории, объясняющей наличие связи между метеорологическими процессами в отдельные сезоны и урожайностью, обусловливает применение в методиках в основном аппарата математической статистики, в частности кластерного и дискриминантного анализа. В результате использования данных процедур каждый год на основе различных характеристик может быть отнесен к одному из нескольких кластеров.
Процедура определения года-аналога осуществляется поэтапно [5]. До сева яровых (в мае) поуже известным погодным данным определяется год-аналог для данного года, по погодному сценарию которого осуществляется прогноз урожайности по разработанной эмпирической модели урожайности и динамической модели МЮС [ 1 ]. Уточнение оценки урожайности производится через две декады после сева и в другие периоды вегетационного периода с корректировкой года-аналога.
Г
ч/ \
1 N Уч.. \
-
х- 1 \
Л" \
1
■ 1 .1. . . Г — 1
iuuik ааямлм iuum ШШ1 шмш niuai 1IIUMI жмз»» «ля«
Рис. 3. Динамика прогнозируемой урожайности, листового индекса и осадков для метеостанции Здвинск Новосибирской области
Идентификация технологии оценки урожайности по эмпирической модели с использованием данных года-аналога была произведена на основе следующих агроклиматических показателей по Алтайскому краю [4]: суммы эффективных температур за вегетационный период и за первые две декады вегетационного периода; суммы осадков за вегетационный период и за первые две декады вегетационного периода; количества дней с осадками за вегетационный период и за первые две декады вегетационного периода; суммы осадков за зимний, весенний и осенний периоды, урожайность пшеницы.
Рассмотрим технологию определения года-аналога и оценку урожайности яровой пшеницы для конкретного объекта — 1997 года (считаем, что погодный сценарий данного года известен на момент определения года-аналога).
По мнению специалистов, огромное влияние на формирование будущего урожая оказывают первые две декады вегетационного периода. Поэтому в работе особое внимание уделено исследованию именно этого периода.
Классификацию осуществляли по фактическим данным первых двух декад вегетационного периода 1997 года. Процедуру кластеризации проводили несколько раз с помощью программы статистической обработки данных 51а115Цса 6.0 при различных значениях числа кластеров. Затем выбиралась лучшая группировка в смысле критерия минимума отношений средних внугрикластерных (/,) и межкластерных (/2) расстояний
1.13
где
т т 1 т
(1)
т(т~0 ;=]>=/+!
т
/=|
Здесь — среднее расстояние между точками внутри /'-го класса; — среднее расстояние между парами точек /-го иу'-го классов.
Для графической интерпретации результатов используется критерий "каменистой осыпи", согласно которому необходимо найти такую абсциссу на графике, в которой линия графика начинает визуально сглаживаться в направлении правой, пологой его части. Согласно графику лучший вариант - разбиение на шесть кластеров (рис. 4).
0,49
<г 0,30 <я
т 0,16 0,04
23456789 Количество кластеров
10 11
Рис. 4. Результаты классификации по первым двум декадам вегетационного периода
Анализ проведенной классификации позволяет предположить, чтодля 1997-го годами-ана-логами будут 1973,1976,1981,1988. Урожайность 1997 года по фактическим данным этих лет будет находиться в пределах 13,4—17,5 ц/га. По погодному сценарию наиболее близкие годы — 1976, 1981, 1988. Следовательно, интервал для ожидаемой урожайности можно сузить: 13,4— 15,9 ц/га.
На следующем этапе осуществлялся прогноз урожайности по эмпирической модели [4]:
=
^т,п+(>0-Кп11п)ПЯ;(д:),еслиК/>У;
/=1
Утп + (^тах - Гу )Г1 ВД,еСЛИ Г; < У, /=1
(2)
где У. + , — урожайность текущего года (ожидаемая); Ур Ктахи Кт|п - соответственно урожайность предыдущего года, максимальная и минимальная; У — средняя урожайность по всему временному ряду; Щх) — нормированные функции отклика; к — количество рассматриваемых факторов, влияющих на урожайность.
Общий вид нормированных функций отклика
Щх) =/(*,
где Х| и хА — минимальное и максимальное значения 1-го фактора, за пределами которых формирование урожая не происходит; х2, х3 — биологически оптимальные границы /-го фактора.
Идентификация структуры и параметров модели (2) проводилась путем спланированного
численного эксперимента, который включал подбор и корректировку структуры уравнений модели с целью минимизации отклонений расчетной урожайности от наблюдаемой. Значения параметров в уравнениях подбирались с помощью оптимизационных процедур.
В качестве параметров модели за первые две декады вегетационного периода, для которых осуществлялась процедура идентификации, выступали: пороговые значения суммы температур при переходе через 5 °С и пороговые значения суммы осадков.
Процент отклонения расчетных данных от фактических по всему временному ряду в среднем составил 20 %. Существенность различий в средней фактической урожайности и расчетной оценивали по критерию Стьюдента: /рас = 2,24; гкр=2,403. При уровне значимости а = 0,02 < /|ф. Следовательно, статистическая проверка показала адекватность модели (2) и что фактическая и прогнозируемая урожайность не отличаются статистически значимо.
Результаты численных расчетов урожайности по модели (2) выбранных годов-аналогов представлены в таблице.
Фактическая и расчетная величина урожайности годов-аналогов, ц/га
Фактическая Расчетная
Год урожайность года- урожайность
аналога по модели (3)
1976 13,7 11,4
1981 13,4 12,4
1988 15,9 12,9
1997 12,8
Тогда урожайность 1997 года по расчетным данным аналогов будет находиться в пределах 11,4—12,9 ц/га, средняя урожайность для 1997 года составит 12,2 ц/га; по фактическим данным аналогов — в пределах 13,4-15,9 ц/га, средняя урожайность для 1997 года 14,6 ц/га (фактическая урожайность — 12,8 ц/га).
Сравнение приведенных результатов показывает, что оценка средней урожайности
по расчетным данным аналогов дает ошибку 0,6 ц/га, по фактическим данным аналогов —
1,8 ц/га.
Аналогичная процедура проводилась по метеорологическим данным годов-аналогов за весь период вегетации. Сравнительный анализ различных разбиений и расчетов по модели (2) за первые две декады вегетационного периода показал, что годом-аналогом для 1997-го является 1981-й. Величина урожая по погодному сценарию 1981 года (фактическая) составит 13,4 ц/га, по модельным расчетам — 12,4 ц/га. Ошибка прогноза по модельным расчетам — 0,4 ц/га, по фактической урожайности 1981 года — 0,6 ц/га.
Итак, разработана технология прогноза урожайности зерновых культур, реализованная в имитационно-моделирующем комплексе программных средств и апробированная в условиях зерновых провинций Западной Сибири. Показано, что разработанный метод прогнозирования урожайности зерновых культур позволяет определять в течение вегетационного периода основные параметры биопродуктивности зерновых культур, давать заблаговременную (за месяц до уборки) оценку урожайности с точностью до 1 ц/га и определять сроки начала уборки зерновых культур, используя модель биопродуктивности, ежедневные данные метеостанций и данные космического мониторинга.
Анализ результатов применения технологии определения года-аналога показал хорошее совпадение фактической урожайности с расчетной по данным двух декад вегетационного периода. Это говорит, во-первых, о том, что действительно первые две декады определяют в значительной степени будущий урожай; во-вторых, что разработанная авторами эмпирическая модель может успешно применяться для текущего планирования и перспективного прогнозирования урожайности яровой пшеницы на территории Алтайского края.
СПИСОК ЛИТЕРАТУРЫ
1 Брыксин В.М. Применение адаптированной модели биопродуктивности ЕРЮ и космоснимков М(Ю15 для прогнозирования урожайности зерно-
вых культур на территории Западной Сибири // Вестник НГУ. Сер. Информационные технологии. Новосибирск: И ад-во НГУ, 2007. Т. 5. Вып. 1. С. 20-26.
2. Матвеев J1.T. Курс общей метеорологии. Физика атмосферы. JI.: Гидрометеоиздат, 1984.
3. Хворова Л .А., Нрыксин В.М, Масленников Д.В.
Оценка возможности использования программного комплекса EPIC для прогноза урожайности в Алтайском крае // Региональные проблемы развития Сибири и Дальнего Востока: IV Сибирский конгресс по прикладной математике ИНПРИМ-2000 / Ин-т математики. Новосибирск, 2000. С. 112.
4. Хворова J1.A., Гавриловская Н.В. Прогнозирование урожайности зерновых культур: методы и расчеты // Изв. АлтГУ. 2008. № 1 (57). С. 65-68.
5. Хворова Л.А., Гавриловская Н.В. Разработка алгоритма определения года-аналога для оценки урожайности зерновых культур в условиях Алтайского края // Изв. АлтГУ. 2007. № 1 (53). С. 66-67.
6. Williams J.R. The EPIC Model / U.S. Department of Agriculture, Agricultural Research Service, Grassland, Soil and Water Research Laboratory. Temple, 1984.
