Описание математической модели и процедуры ее идентификации

Математическое и компьютерное моделирование продуктивности растений в зависимости от динамики влажности почвы

В. М. Казиев, С. К. Кирьязева, Д. А. Кирьязев

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

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

Описание математической модели и процедуры ее идентификации

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

W (t) = q(t)P(t) + P₁(t) — E(t) — H(t), (1)

где P(t) - величина осадков; q(t) - коэффициент использования осадков (определяется, например, экспертно или по формуле Харченко С.И. [6], через W_min - наименьшую влагоёмкость почвы и W_z - влажность завядания); P₁(t) - подпитывание (приток) из грунтовых вод; E(t) - суммарное испарение с корнеобитаемого слоя; H(t) - уровень (сток) грунтовых вод, W(t) - средняя по слою влажность почвы (с учётом поливов или на межполивной период).

Оценим и учтём влияние накопившейся к некоторому моменту времени биомассы растений на экологически обоснованную величину суммарного испарения в каждый момент времени.

Величину суммарного испарения с корнеобитаемой зоны растений представим в виде суммы интенсивности транспирации растениями E₁(t) и интенсивности испарения с поверхности почвы E₀(t):

E(t) = E₀(t) + E₁(t) (2)

Известно [1], что динамика прироста биомассы в предположении, что прирост биомассы хорошо коррелирует с интенсивностью транспирации растительного покрова, описывается уравнением:

x (t) = a(t)E₁(t) — b(t)x(t), (3)

где x(t) - биомасса культуры; a(t) - эффективность транспирации; b(t) - коэффициент расхода на дыхание.

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

В рассматриваемой нами процедуре моделирования будем использовать

достаточно простую для идентификации модель Ферхюльста - Вольтерра [2]:

x (t) = [ (t) —  (t)x(t)]x(t), (4)

где  - коэффициент роста (автоприроста),  - коэффициент сопротивления среды (нехватки воды).

Известно также, что динамика прироста биомассы хорошо описывается уравнением Давидсона - Филиппа [3]:

x (t)= e₀(t)(F(t) — R(t)), (5)

где e₀ - коэффициент перехода от массы усвоенной СО₂ к сухой фитомассе; F - суммарный фотосинтез растений; R - суммарное дыхание растений.

Интенсивность дыхания за сутки зависит от величины накопившейся биомассы. Экспериментально получено [4], что

R(t) = b(t)x(t) + e₁F(t), (6)

где e₁ — коэффициент затрат на рост биомассы растений.

Коэффициенты  ₀ ,  ₁ - экспериментально определяемые; для большинства культур можно полагать  ₀ = 0,68,  ₁ = 0,27.

С учетом уравнений (1) — (4) имеем следующую модель расчета влажности почвы с учетом динамики накапливаемой биомассы:

W (t) = q(t)P(t) + P₁(t) — E(t) — H(t),

(7)

E_1 (t) =[ (t) —  (t)x(t) + b(t)]x(t)/a(t),

Из (4), (5) и (6) имеем:

b(t) = (1-e₁)F(t)/x(t) — ( (t) —  (t)x(t))/e₀, (8)

Для нахождения влажности почвы нам необходимо идентифицировать переменные  и  . Эта задача достаточно сложна из-за сложности и дороговизны проведения экспериментальных исследований (мониторинга). Мы продемонстрируем имитационную процедуру её решения для случая постоянных параметров модели (4); случай кусочно-постоянных параметров - аналогичен и влияет только на размерность задачи, а случай произвольных функции сводим к проблеме аппроксимации их некоторой системой базисных функций.

Решение уравнения (4), как легко проверить, имеет вид:

x(t)= ( + Ce ^{—  t}) . (9)

Теперь для того, чтобы найти  и  нужно, согласно метода наименьших квадратов, решить задачу минимизации квадратичного функционала вида:

n

f( ,  , c) =  (x_i0 — x_i)²  min, (1)

i=1

где i - номер фазы вегетации растения (i=1,2,...,n); n - число фаз вегетации; x_i0 - экспериментальные величины урожайности культуры за репрезентативный период времени; x_i - теоретические величины урожайности сельхозкультур, определяемые по формуле (9).

Для нахождения решения задачи (10) необходимо решить нелинейную систему уравнений:

df / d =0, df / d = 0, df / dc = 0. (11)

Решаем эту систему численно (например, методом Зейделя), с требуемой точностью  и критерием адекватности вида:

( _i+1 —  _i)² + ( _i+1 — _i)² + (с_i+1 — с_i)² < ² .

Величина фотосинтеза определяется по формуле вида:

F(t)=F_max e ^{—  [s(t) — z]}[ (t)x(t)/ (t)]^2/3,

где s(t) — текущая сумма биологически активных температур, z — сумма биологически активных температур для максимального развития листовой поверхности,  — эмпирический коэффициент, зависящий от .

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

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

 (12)

где x(W) - прогнозная урожайность; x_max - максимальная урожайность сельхозкультур; W - влагообеспеченность корнеобитаемого слоя почвы, определяемая как описано выше; W_min, W_max - соответственно нижняя и верхняя границы влагообеспеченности почвы, при которой урожай равен нулю; W_opt - влагообеспеченность, соответствующая x_max;  - параметр, характеризующий темпы роста урожая с увеличением влагообеспеченности (параметр саморегуляции системы).