Математическое моделирование роста растений на основе экспериментальных данных Текст научной статьи по специальности «Сельское хозяйство, лесное хозяйство, рыбное хозяйство»

СЕЛЬСКОХОЗЯЙСТВЕННАЯ БИОЛОГИЯ, 2007, № 1

Методика

УДК 633:581.14:58.087

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ РОСТА РАСТЕНИЙ НА ОСНОВЕ ЭКСПЕРИМЕНТАЛЬНЫХ ДАННЫХ

И.М. МИХАЙЛЕНКО

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

Ускоренное развитие информационных технологий, вычислительной техники и успехи математических дисциплин предоставляют новые возможности ученым-физиологам и аграриям для разработки теории управления продуктивными процессами культурных растений. Эта теория давно востребована при создании современных систем управления технологическими процессами в растениеводстве. То же самое в полной мере относится и к проблеме точного земледелия, областью применения которого в большей степени является растениеводство в открытом грунте. Еще более актуально развитие теории управления продуктивностью для растениеводства в защищенном грунте и установках искусственного климата.

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

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

103

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

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

М о д е л ь р о с т а р а с т е н и й о з и м о й п ш е н и ц ы в ч а с о в о м м а с ш т а б е в р е м е н и. Для обоснования структуры математической модели при исследовании суточных (циркадных) закономерностей роста растений в часовом масштабе времени мы использовали экспериментальные данные Шевелухи (1, 2). В качестве основного состояния растений мы приняли скорость роста в единицу времени (V, мм/ч), а в качестве влияющих факторов — температуру (T, °C), относительную влажность воздуха (w, %) и продолжительность солнечного сияния (% от часа R). На основании анализа экспериментальных данных базовую форму (структуру) модели суточного роста растений c использованием алгоритмов самоорганизации можно представить в форме дифференциального уравнения следующего вида (3, 4): dx

X = d- ax + ф(Ь, T) + с^(т) + c2R(t) + f (p, x) + §(x),

[1]

X Є (0,24), x(0) = 0,

где x — время суток, ч; x — выбранный показатель линейного роста совокупности (популяции) растений, мм; a — динамический параметр, учитывающий инерционность процесса, 1/ч; ci, С2 — параметры, отображающие зависимость роста от влияния различных факторов, 1/ч; p (b, T) — нелинейная часть модели, в которой учитывается неоднозначное влияние температуры в областях ниже и выше оптимальной; f (p, т) — функция времени, задающая внутренние колебания роста; b, p — параметры функций; |(г) — случайная ошибка моделирования, имеющая нулевое среднее и неизвестную дисперсию d.

Нелинейную часть модели p (b, T) мы аппроксимировали кусочно*

линейной функцией следующего вида: cp(b, T) = bi T, если Т < Т1; cp(b, T) =

= b2(T*2 — Т)*, если Т > Т2; p(b,T) = b3, если T є (Ті, Т2) в которой T*,

*

T2 — границы области оптимальных для линейного роста значений температуры, bi, b2 и Ьз — параметры выбранной аппроксимации. В качестве функции f(p, т), задающей внутренние колебания роста, была выбрана гармоническая функция с неизвестными параметрами pi, p2 и рз — f(p, т) = Pi Sin(T/p2 + Рз).

Как обычно принято в процедурах идентификации, модель состояния объекта [i] обязательно должна быть дополнена моделью измерителя состояния (в нашем случае показателя роста растения) — у(т) = hx(T) + + |(г), где у(г) — показания ростомера (ауксанограф); h — параметр преобразования ростомера; |(г) — ошибка измерения показателя роста, имеющая нулевую среднюю и известную дисперсию ст2.

Предложенная модель может описывать процесс изменения показателя для оценки неизвестных параметров в суточном интервале времени i04

в среднем для совокупности растений на выбранном участке поля или экспериментальной установке.

М о д е л и р о с т а м н о г о л е т н и х т р а в в с у т о ч-н о м м а с ш т а б е в р е м е н и. Модели циркадных ритмов роста культурных растений представляют большой интерес с исследовательской точки зрения, но для решения задач управления продукционным процессом в суточном масштабе времени они неприменимы. Поэтому рассмотрим пример идентификации модели роста многолетних трав, которые мы использовали для решения оперативных задач управления процессами заготовки кормов в хозяйствах Ленинградской области.

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

Среди множества факторов, оказывающих влияние на состояние травостоя, мы выделяли только климатические, которые рассматривали в модели как внешние возмущения: / — среднесуточная температура воздушной среды, оС; / — среднесуточный уровень радиации, Вт/(м2-ч); / — среднесуточное количество осадков, мм. Кроме рассмотренных внешних возмущений в модель вводили еще одно дополнительное возмущение, которым является календарное время — / = т, сут, чем косвенно учитывали влияние состояния почвы, которое непосредственно не определяли. В реальных условиях вышеуказанные состояния не могут измеряться непосредственно, а вместо них мы использовали производные величины: уім — общая надземная биомасса растений, кг/м2; у2м — содержание сухого вещества в общей биомассе, %.

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

х

ХІм

2м

a a х + c с с с

a11 a12 х1м с11 с12 с13 с14

21 21J _ 2м _ _ 21 22 23 24

/11 (т)

Ш

[2]

гє(0,7); хім(0) = хм; Х2м = Х20;

У1 (т) 1м 1

1 Н 1 1

1

Щ1,«)

1 У1 (т) 1м

1 1 Н 1

[3]

где выражение [3] задает параметры уравнения оценки состояния травостоя.

Для идентификации модели роста [1] в часовом времени мы вынуждены были использовать статический (регрессионный) метод идентификации. Это связано с тем, что в используемых нами данных (1, 2) представлены только показатели скорости роста (V) и влияющих факторов, являющихся соответственно выходом и входом модели. При этом измеряемый выход — статический показатель, так как скорость роста не может служить решением динамической модели роста [1]. Решением для модели является абсолютное значение показателя роста — х(г), которое в экспериментах не измеряли. Для реализации регрессионного метода исходную

105

непрерывную модель [1] мы сначала преобразовывали к модели с дискретным временем:

x(k + l) = PTW(k), [4]

где в качестве выхода служит скорость роста — x(k) = V(k), а вектором входов являются все факторы, влияющие на рост — WT(k) =

= [T(k), w(k), R(k)], где k — дискретные отрезки времени, ч.

Для оценки качества идентификации модели [1] мы использовали квадратическую функцию штрафа:

1 24 ,

Ip = — 1 M[(y(z) - hx (т| W)2dr]. [5]

24 0

Для векторного выхода модели [5] эта функция имеет следующий

вид:

1 Т

Ib = -1 М[(¥(т) - HX(t))t(Y(t) - HX(r))]-dr, [6]

Т 0

где М[-] — операция математического ожидания во времени; Y(t) — вектор измеряемых величин; Х(т) — вектор состояния модели [2]; Н — матрица преобразования измерителя [3]; T — транспонирование вектора или матрицы.

Для оценки параметров преобразованной модели, обеспечивающих минимум критерию [5], мы использовали следующий регрессионный алгоритм идентификации (6):

G(k) = R(k - 1) W(k), z(k) = WT(k)G{k),

e(k) = y(k + 1) - PT(k - 1) W(k), [7]

е (k)

P(k) = P(k - 1) + —4-g (k),

1 + z (k)

R(k) = R(k -1) - 1 G(k)GT (k),

1 + z(k)

где G(k), R(k) и z(k) — вспомогательные соответственно вектор, матрица и скаляр, являющиеся параметрами алгоритма.

Идентификацию модели роста биомассы [5] проводили по наблюдениям за состоянием травостоя, которые осуществляли в системе мониторинга кормовых угодий Ленинградской области в течение 2003-2005 годов. Интервал наблюдения от начала отрастания травостоя до первого укоса обычно составлял 21-28 сут, начиная ежегодно с 20-25 мая. Текущее измерение состояния проводили на реперных полях хозяйств 2 раза в неделю на заранее выделенных делянках посредством взятия проб с последующей их доставкой в централизованную лабораторию, где анализировали количественные и качественные показатели травостоя (5).

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

106

тегральное выражение функции штрафа и правую часть идентифицируемой модели:

Н = [(7(т) - НХ(т | П))г (7(т) - HX(x| П))] + ут(т)Ф(Х| Д,П), [8]

где Ф(Х|Д,П) обозначает всю правую часть модели [1], содержащую вектор неизвестных параметров П; у(т) — вектор сопряженных переменных, определяемый по модели следующего вида:

у = -

дН (т)

дХ ’

[9]

т є (T, 0), т(T) = 0.

На основании функции Гамильтона [8] и модели для сопряженной переменной [9] мы реализовали алгоритм поиска значений вектора параметров модели [1], обеспечивающий минимум функции штрафа [6] (7):

— Шаг 1. Задают циклическую переменную m = 1; принимают начальное значение вектора параметров модели [1] — Пт; вводят массив внешних факторов, объединенных в вектор Дт); решают систему [2] в ин-

тервале наблюдения т є (0,T); массив решений Х(т)т фиксируют в памяти.

— Шаг 2. Массив решений для основной модели Х(т)т подставляют в гамильтониан [8] и в обратном времени решают уравнение для сопряженной переменной [9]; массив решений ¥(-т)т фиксируют в памяти.

— Шаг 3. Массивы решений Х(т)т и У(-т)т подставляют в гамильтониан [8] и рассчитывают вектор частных производных по всем неизвестным параметрам модели [1]:

Gm (т) =

дД(т)

[10]

m ' ягг ді1 m

Из всего массива решений [10] запоминают значения в начальный момент времени Gm(0).

— Шаг 4. Определяют шаговую процедуру уточнения вектора параметров модели

Пт+1 = Пт — ^mGm(0), [11]

для которой определяют оптимальную длину шага д m из условия минимума критерия [6] по величине шага.

— Шаг 5. С оптимальной длиной шага д * уточняют вектор пара-

метров модели [2] — П = П - Д* G (0) .

* L 1 m+l m m my '

— Шаг 6. Решают системы [6] и [9] и находят новый вектор градиентов Gm+1(0).

— Шаг 7. Определяют направление сопряженных градиентов:

m +1

G(0) , +

v 'm + 1

||g (0)

m +1 I

G (0)m

2

2

[12]

а затем осуществляют переход к шагу 4, вплоть до достижения сходимости.

Этот алгоритм является базовым для всех вариантов критериев качества идентификации модели и динамических моделей любой размерности формы. Различия касаются вида функции Гамильтона, в которую вхо-

107

дит подынтегральное выражение для соответствующего критерия и правая часть динамической модели объекта или системы.

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

Рис. 1. Исходные данные для идентификации суточной модели роста растений озимой пшеницы:

1, 2 и 3 — соответственно относительная влажность воздуха, %, почасовая продолжительность солнцестояния, %/ч и температура наружного воздуха, оС. Цит. по ист. 1: А и Б — соответственно стр. 64 и 71.

На рисунке 2 приведены результаты идентификации суточной модели роста [1], для которой получены следующие оценки параметров:

— Для данных на рисунке 1А: с учетом инерционности роста — a = —0,0982; b\ = 0,072; c\ = 0,0068; c2 = 0,0184; без учета инерционности роста — a = 0; b\ = 0,0627; c\ = 0,0064; С2 = 0,0171; величина критерия [5] составляет 0,0287 мм/ч2, (среднеквадратическая ошибка прогнозирования +5 %).

— Для данных на рисунке 2А в соответствии с линейной аппроксимацией (кусочно-линейная функция) при Т* = 16 оС и Т2* = 20 оС:

с учетом инерционности роста — a = —0,0783; b\ = 0,0315; b2 = —0,0763; Ьз = 0,0546; С1 = 0,0065; С2 = 0,0076; без учета инерционности роста — a = 0; bj = 0,0389; b2 = —0,0844; b3 = 0,069; c\ = 0,0077; c2 = 0,0092; величина критерия [5] составляет 0,0069 мм/ч2 (среднеквадратическая ошибка прогнозирования +1,5 %).

На рисунке 3 представлены результаты идентификации модели накопления биомассы многолетних трав в суточном масштабе времени. Величина критерия качества идентификации [6] составляет 10,507 ц/га2 (точность прогнозирования на недельный интервал времени +2,5 %).

108

Рис. 2. Результаты идентификации суточной модели роста растений озимой пшеницы: А и Б —

для моделей, представленных соответственно на рис. 1А и рис. 1Б; 1 и 2 — скорость роста соответственно моделируемая и измеренная, мм.

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

Рис. 3. Результаты идентификации модели накопления сухой массы многолетних трав (тимофеевка) в суточном масштабе времени (Волховский р-н, 2005 год): ■ и Д — соответственно данные измерений и прогноз.

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

109

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

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

1T .

Ibi = -J [№) - HXi(r\щ, nMY(r) - hx(t\ щ,щ + a№l щ, П - 1)

- Xi(r\ Щ,П i))TXi(r\ Щ,П i - 1) - Xi(r\ Щ,П i))]dz,

0

[13]

где Xi(z\ Щі,Пі - 1) — прогноз вектора состояния модели [2] в текущий период наблюдения при оценках параметров, полученных в предыдущем интервале; gi — весовой множитель, определяющий степень компромисса между точностью текущей настройки вектора параметров и регулярными свойствами модели и являющийся дополнительным оптимизируемым параметром алгоритма идентификации.

Введение дополнительной штрафной компоненты изменило только выражение для функции Гамильтона, не меняя структуру алгоритма. В то же время для определения оптимального значения весового множителя gi была введена дополнительная процедура, выполняемая между отдельными интервалами идентификации:

— Шаг 8. Вычисляют критерий качества прогнозирования:

1 N

KPi = — I[Xi(r\ Щ,Пі - 0 - Yi(r)]T[Xi (т|Щ-, Пі - 0 - Y(r)]dr.

Nt = 1

— Шаг 9. Уточняют параметр gi:

[14]

gi = gi - 1, если KPi - KPi - 1 < 0, gi = Pgi - 1, если KPi - KPi - 1 > 0,

[15]

где p = 0,1 - 0,6 — параметр алгоритма.

— Шаг 9. Если \KPt - KPt - 11| < 5, то стоп, иначе — переход к очередному интервалу измерения и идентификации.

Модификация алгоритма идентификации позволила обеспечить эффективную автоматическую адаптацию параметров модели накопления биомассы и выход на требуемые показатели точности прогнозирования (см. рис. 3).

110

По этой причине мы не приводим конкретных оценок параметров

модели. Этот прием может быть применен и к любым другим моделям

роста, когда их используют для межинтервального прогнозирования.

Л И Т Е Р А Т У Р А

1. Ш е в е л у х а B.C. Рост растений и его регуляция в онтогенезе. М., 1992.

2. Ш е в е л у х а B.C. Периодичность роста сельскохозяйственных растений и пути ее регулирования. М., 1980.

3. И в а х н е н к о А.Г., Ю р а ч к о в с к и й Ю.П. Моделирование сложных систем по экспериментальным данным. М., 1987.

4. И в а х н е н к о А.Г., М ю л л е р И.А. Самоорганизация прогнозирующих моделей. Киев, 1985.

5. М и х а й л е н к о И.М., К у р а ш в и л и А.Е. Управление качеством кормов из многолетних трав. Мат. Междунар. науч.-практ. конф. «Информационные технологии, информационные измерительные системы и приборы в исследовании сельскохозяйственных процессов». Новосибирск, 2003.

6. Современные методы идентификации систем /Под ред. П. Эйкхоффа. М., 1983.

7. С е й д ж Э.П., М э л с а Дж. Идентификация систем управления. М., 1974.

Агрофизический НИИ, Поступила в редакцию

195220, С.-Петербург, Гражданский просп., 14; 30 августа 2005 года

e-mail: offis.@agrophis.ru.

MATHEMATICAL MODELLING OF PLANT GROWTH ON THE BASIS OF EXPERIMENTAL DATA

I.M. Mikhailenko S u m m a r y

The author considers the problem of identification of mathematical growth models used for forecasting and control of plant productivity. The great attention the author devoted to effective methods and algorithms of estimation of parameters of growth models on experimental data. The method and software on identification of mathematic models of plant growth at the hour and daily of time scale were founded. The use of proposed models for solution of different research and science-practical problems was discussed.

Новые книги

П е т р и ч е н к о В.М., С у х и н и -н а Т.В. Очанки Западного Урала (фармакогностические и биологические аспекты).

Пермь: ГОУ ВПО «ПГФА Росздрава», 2006, 145 с.

В монографии обобщены и систематизированы данные литературы и результаты собственных исследований авторов растений рода Очанки (Euphrasia L.) сем. Scrophulariaceae. Описано морфолого-анатомическое строение, химический состав и фармакологические свойства наиболее распространенных на Западном Урале видов. Приведены сведения о биологических свойствах семян. Дана фитоценотипическая и ресурсоведческая характеристика очанки коротковолосистой. Рассматриваются технологические аспеты получения и стандартизации фитопрепаратов из растений очанки разных видов. Уделено внимание оценке содержания биологически активных веществ в растениях рода Очанка.

Б а й к о в а Е.В. Род шалфей: морфология, эволюция, перспективы интродукции.

Новосибирск: Сибирская издательская

фирма «Наука», 2006, 248 с.

В монографии обобщены данные отечественной и зарубежной литературы по морфологии, географии, эволюции растений рода шалфей (Salvia L.). Приведены результаты исследования 43 модельных видов при интродукции в условиях Западной Сибири и около 600 видов по материалам крупнейших гербариев России. Представлены новые данные по онтогенезу, структуре побеговых систем и соцветий растений шалфея, морфологии андроцея и трихом, ультраструктуре поверхности эремов и особенностям ослиз-нения перикарпия. Сформулирована гипотеза эволюции рода, разработана оригинальная система его жизненных форм. Рекомендовано 27 видов для использования в озеленении лесостепной зоны Западной Сибири.

111