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

Решение задач, связанных с оценкой состояния сельскохозяйственных растений в период вегетации, позволяет эффективно использовать удобрения, получать выгодную урожайность, улучшать качественные характеристики растений, а также экологическое состояние поля. Все более перспективным направлением в точном земледелии становится применение различных методов математической статистики. Нашей целью была оценка состояния сельскохозяйственных растений с использованием подхода, основанного на совместном применении методов кригинга и бинарной регрессии, а также определение обеспеченности посевов азотом по индексу NDVI (Normalized Difference Vegetation Index). Исследования проводили на участке опытного поля, расположенного на территории филиала Агрофизического института (д. Меньково, Ленинградская обл.) в 2015 году. С помощью аэрофотоснимков, сделанных с автоматизированного беспилотного летательного комплекса Геоскан-401 («Геоскан», Россия), был получен набор значений вегетационного индекса NDVI в произвольных точках участка. Также провели ряд наземных измерений. Предложенный подход к оценке состояния сельскохозяйственных растений заключался в совместном использовании двух методов математической статистики: ординарного кригинга и логистической регрессии. Осуществляли предварительный вариограммный анализ, выполняли построение модели вариограммы. После этого методом кригинга вычисляли ряд прогнозируемых значений исследуемого параметра. На следующем этапе устанавливали пороговое значение параметра для исследуемого участка, а также вводили фиктивную переменную, принимающую значение 1, если величина параметра превышала порог, и 0 в противном случае. Затем строили логит-модель, в которой одним из факторов был ряд оценок интересующего параметра, полученных с использованием метода ординарного кригинга. Входные данные для построения логит-моделей были следующими: N(x_i) — значение NDVI в местоположении x_i, i = 1,78; переменная T = 1, если N(x_i) ≥ 0,46, иначе T = 0; переменные X и Y — координаты наблюдений, рассматриваются в качестве объясняющих переменных; N_pred(x_i) — величины параметра, спрогнозированные с использованием метода кригинга в наблюдаемых точках. Все вычисления выполнялись с помощью языка программирования R. В результате проведенного эксперимента были построены три логит-модели с зависимой переменной T: в первой модели объясняющие переменные X и Y, во второй модели — X, Y и N_pred, в третьей модели — N_pred. Тестирование показало, что при добавлении переменной N_pred логит-модель работает лучше (в 2 раза меньше ошибочного определения исследуемого параметра).Полученные результаты дают основание полагать, что добавление в факторы бинарной регрессии набора спрогнозированных методом кригинга величин позволяет существенно повысить точность расчетов.

Ключевые слова: оценка состояния растений, Normalized Difference Vegetation Index, NDVI, кригинг, бинарная регрессия, язык R.

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

Один из новых подходов в агрофизике строится на методах, использующих бинарную регрессию. Например, в работе В.М. Буре (7) описано ее применение для прогноза урожайности. В исследованиях, проведенных в Норвегии и Нидерландах, предложены методики прогноза пространственного распределения типов почв на основе мультиномиальной логистической регрессии с использованием цифрового анализа местности (8, 9). К более изученным направлениям в точном земледелии относится геостатистика, которая позволяет картировать содержание в почве питательных веществ (азот, фосфор, калий и др.) (10, 11), а также электропроводность, pH, плотность, влажность (12-14) с целью оптимального управления хозяйством. При геостатистическом подходе почва рассматривается как набор переменных, непрерывных в пространстве, и их изменения описываются с точки зрения пространственной зависимости (15, 16). Сочетание методов геостатистики и бинарной регрессии предложено только в экономике (17), для задач точного земледелия подобные методологии не описаны.

В представленной работе мы впервые спрогнозировали пространственное распределение значений вегетационного индекса NDVI (Nor-malized Difference Vegetation Index) на опытном поле, применив регрессию на основе гауссовских процессов (кригинг) совместно с логистической регрессией как вариантом бинарной регрессии. При этом результаты тестирования показали, что предложенный подход достаточно точно оценивает интересующий параметр на исследуемом участке поля.

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

Методика. Исследования проводили на участке опытного поля на территории филиала Агрофизического института (д. Меньково, Ленинградская обл.) в 2015 году. С помощью аэрофотоснимков, сделанных с автоматизированного беспилотного летательного комплекса Геоскан-401 («Геоскан», Россия), был получен набор значений вегетационного индекса NDVI (Normalized Difference Vegetation Index) в произвольных точках участка (всего 78 значений).

Для анализа состояния сельскохозяйственных растений мы применяли совместно два метода математической статистики — ординарный кригинг и логистическую регрессию. При этом использовали логит-модель (logit-model) как подход, обеспечивающий простоту оценки параметров.

Пространственное распределение интересующего параметра прогнозировали с помощью ординарного кригинга для набора измерений (18, 19):

где n — число измерений, Z(x_i) — значение наблюдаемого параметра в местоположении x_i, λ_i — неизвестный вес для параметра, — прогноз значения параметра в местоположении x₀.

Неизвестный вес находили с помощью вариограммного анализа с построением теоретической модели вариограммы γ(h) на основании полученной экспериментальной кривой .

Для выполнения логистической регрессии фиксировали значение d (порог) для исследуемого участка поля и вводили фиктивную переменную:

При этом значения переменной y(x_i) в наблюдаемых точках были известны, поскольку в этих точках известны величины интересующего параметра и его уровень относительно порогового значения. Кроме того, на предыдущем этапе мы получали набор спрогнозированных величин параметра с помощью кригинга. В результате был cформирован набор исходных данных для логистической регрессии, которая отражала зависимость между вероятностью превышения значения порога и объясняющими переменными (20, 21):

где φ_i— факторы, объясняющие фиктивную переменную y(x_i).

В качестве одного из факторов логит-модели использовали набор значений, спрогнозированных методом кригинга. Адекватность построенной логит-модели [2] проверяли с помощью классических приемов статистического анализа — теста Вальда (Wald test, W) и теста отношения правдоподобия (likelihood ratio test, LR) (22).
Для каждой точки исследуемого участка вычисляли вероятность P(y(x) = 1), которая при приближении к 1 указывала на превышение значения интересующего параметра порогового уровня d, к 0 — на нахождение значения параметра ниже порога.

Расчеты осуществляли с помощью языка программирования R (https://www.r-project.org), активно используемого при решении задач точного земледелия (23).

Результаты. На рисунке 1 представлен аэрофотоснимок обследуемого участка поля, а также карта распределения собранных данных (их исходным значениям пропорциональны диаметры кругов), построенная в программе R-statistics. В результате был получен набор значений вегетационного индекса NDVI, величина которого, как известно, коррелирует с содержанием азота в растении в рассматриваемой точке (24, 25).

Суть предложенного нами подхода заключается в применении в качестве одного из факторов модели набора спрогнозированных значений, который получали методом кригинга. В соответствии с этим на первом этапе (при прогнозе пространственного распределения интересующего параметра) проверяли выполнение условий геостатистики — стационарности и мультинормальности (26). Обнаруженные выбросы данных отсекали по 2,5 % двусторонним квантилям. Кроме того, оценивали линейную корреляцию интересующего параметра с координатами. В результате пространственный тренд не был выявлен. Проверка с помощью критерия Кол-могорова-Смирнова не дала оснований отвергнуть гипотезу о нормальном распределении (достигаемый уровень значимости равен 89,75 %).

На следующем этапе был выполнен вариограммный анализ и построена модель вариограммы с использованием функции vgm. На рисунке 2 представлена полученная экспериментальная вариограмма по четырем направлениям (0, 90, 135 и 270°) с установленной моделью вариограммы. На ее основе проводили ординарный кригинг [1]: из набора исходных наблюдений поочередно убирали значение в одной точке, после этого осуществляли его прогноз методом кригинга с использованием функции krige.

В результате с помощью функции glm были получены три логит-модели. Величину d = 0,46 выбрали пороговой. Значимость построенных моделей оценивали критерием отношения правдоподобия. Входные данные для построения логит-моделей были следующими: N(x_i) — значение NDVI в местоположении x_i, i = 1,78; переменная T = 1, если N(x_i) ≥ 0,46, иначе T = 0; переменные X и Y — координаты наблюдений, рассматриваются в качестве объясняющих переменных; N_pred(x_i) — величины параметра, спрогнозированные методом кригинга в наблюдаемых точках.

Оценки коэффициентов второй логит-модели была в целом незначимы. Величину a при этом принимали равной 0,05 (табл. 1). Статистический анализ подтвердил, что уравнение этой логит-модели в целом незначимо, в то время как уравнения первой и третьей моделей значимы, при этом их тестирование показало, что третья модель работает лучше.

На следующем этапе для анализа адекватности трех полученных моделей из набора исходных наблюдений поочередно исключали одну точку, строили логит-модели заново и оценивали, как каждая из моделей спрогнозировала вероятность превышения порогового значения в отброшенной точке. Тестирование показало, что первая модель ошибочно определила величину параметра в 26 точках (33,3 %), а третья — в 12 точках (15,38 %) (табл. 2).

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

Похожие результаты были получены G.B. Fernandes с соавт. (17), которые изучали логит-модель кредитного скоринга с пространственной переменной в качестве объясняющей. При этом сравнивались две модели: с добавлением пространственной переменной и без нее. Результаты показали, что предложенная авторами методика эффективнее классических. В нашей работе исследовался подход для прогноза пространственного распределения интересующего параметра, основанный на совместном применении кригинга и бинарной регрессии, при этом полная модель (с включением в логит-модель набора величин параметра, спрогнозированных методом кригинга) оказалась лучше альтернативных. Следует отметить, что в России стали применять метод бинарной регрессии в задачах точного земледелия лишь с 2000-х годов. В некоторых сообщениях (7, 20) подробно рассматриваются возможности использования логит- и пробит-моделей в растениеводстве, однако пространственная переменная при этом не учитывается.

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

ЛИТЕРАТУРА

1. Bure V.M., Mitrofanova O.A. Analysis of aerial photographs to predict the spatial distribution of ecological data. Contemporary Engineering Sciences, 2017, 10(4): 157-163 (doi: 10.12988/ces.2017.611175).
2. Kim Y., Reid J.F., Han S. On-the-go nitrogen sensing and fertilizer control for site-specific crop management. International Journal of Agricultural and Biosystems Engineering, 2006, 7(1): 18-26.
3. Blackmer T.M., Schepers J.S. Aerial photography to detect nitrogen stress in corn. J. Plant Physiol., 1996, 148(3-4): 440-444 (doi: 10.1016/S0176-1617(96)80277-X).
4. Graeff S., Pfenning J., Claupein W., Liebig H.P. Evaluation of image analysis to determine the N-fertilizer demand of broccoli plants. Advances in Optical Technologies, 2008, Article ID 359760 (doi: 10.1155/2008/359760).
5. Thenkabail P.S., Lyon J.G., Huete A. Hyperspectral remote sensing of vegetation. CRC Press, MA, USA, 2011.
6. Franzen D.W., Reitmeier L., Giles J.F., Cattanach A.C. Aerial photography and satellite imagery to detect deep soil nitrogen levels in potato and sugarbeet. Proc. of the 4th International Conference «Precision Agriculture». St. Paul, MN, 1999: 281-290.
7. Буре В.М. Методология применения бинарной регрессии в точном земледелии. Математические модели в теоретической экологии и земледелии. Мат. Межд. семинара, посвященного памяти профессора Ратмира Александровича Полуэктова (Полуэктовские чтения). СПб, 2014: 118-121.
8. Debella-Gilo M., Etzelmüller B. Spatial prediction of soil classes using digital terrain analysis and multinomial logistic regression modeling integrated in GIS: Examples from Vestfold County, Norway. Catena, 2009, 77(1): 8-18 (doi: 10.1016/j.catena.2008.12.001).
9. Kempen B., Brus D.J., Heuvelink G.B.M., Stoorvogel J.J. Updating the 1:50,000 Dutch soil map using legacy soil data: a multinomial logistic regression approach. Geoderma, 2009, 151(3-4): 311-326 (doi: 10.1016/j.geoderma.2009.04.023).
10. Якушев В.П., Жуковский Е.Е., Кабанец А.Л., Петрушин А.Ф., Якушев В.В. Вариограммный анализ пространственной неоднородности сельскохозяйственных полей для целей точного земледелия. СПб, 2010.
11. Claret M.M., Urrutia R.P., Ortega R.B., Best S.S., Valderrama N.V. Quantifying nitrate leaching in irrigated wheat with different nitrogen fertilization strategies in an Alfisol. Chilean Journal of Agricultural Research, 2011, 71(1): 148-156 (doi: 10.4067/S0718-58392011000100018).
12. Krasilnikov P., Sidorova V. Geostatistical analysis of the spatial structure of acidity and organic carbon in zonal soils of the Russian plain. In: Soil geography and geostatistics: concepts and applications /P. Krasilnikov, F. Carré, L. Montanarella (eds). European Communities, Luxembourg, 2008: 55-67.
13. Dulaney W.P., Lengnick L.L., Hart G.F. Use of geostatistical techniques in the design of an agricultural field experimental. Proc. of the Survey research methods section. Alexandria, VA, 1994: 183-187.
14. Adamsen F.J., Pinter P.J., Barnes E.M. Measuring wheat senescence with a digital camera. Crop Science, 1999, 39(3): 719-724 (doi: 10.2135/cropsci1999.0011183X003900030019x).
15. Panayi E., Peters G.W., Kyriakides G. Statistical modelling for precision agriculture: a case study in optimal environmental schedules for Agaricus Bisporus production via variable domain functional regression. PLoS ONE, 2017, 12(9): e0181921 (doi: 10.1371/journal.pone.0181921).
16. Isaaks E., Srivastava M. An introduction to applied geostatistics. Oxford University Press, NY, USA, 1989.
17. Fernandes G.B., Artes R. Spatial dependence in credit risk and its improvement in credit scoring. Eur. J. Oper. Res., 2016, 249(2): 517-524 (doi: 10.1016/j.ejor.2015.07.013).
18. Демьянов В., Савельева Е. Геостатистика. Теорияипрактика. М., 2010. 
19. Webster R., Oliver M.A. Geostatistics for environmental scientists, Second edition. John Wiley and Sons Ltd, Chichester, UK, 2007 (doi: 10.1002/9780470517277).
20. Якушев В.П., Буре В.М., Парилина Е.М. Бинарная регрессия и ее применение в агрофизике. СПб, 2015.
21. Hosmer D., Lemeshow S. Applied logistic regression, Second edition. Wiley, NY, 2000 (doi: 10.1002/0471722146).
22. Буре В.М., Парилина Е.М. Теория вероятностей и математическая статистика: учебное пособие. СПб, 2013.
23. Plant R.E. Spatial data analysis in ecology and agriculture using R. CRC Press, Boca Raton, 2012.
24. Якушев В.П., Канаш Е.В., Конев А.А., Ковтюх С.Н., Лекомцев П.В., Матвеенко Д.А., Петрушин А.Ф., Якушев В.В., Буре В.М., Русаков Д.В., Осипов Ю.А. Теоретические и методические основы выделения однородных технологических зон для дифференцированного применения средств химизации по оптическим характеристикам посева: практическое пособие. СПб, 2010.
25. Kanash E.V., Osipov Ju.A. Optical signals of oxidative stress in crops physiological state diagnostics. Proc. 7th European conference on precision agriculture. Wageningen, Netherlands, 2009: 81-89.
26. Goovaerts P. Geostatistics for natural resources evaluation. Oxford University Press, NY, 1997.

Solving problems related to the assessment of the status of agricultural plants during the growing season, allows us to effectively use fertilizers, obtain favorable yields, improve the quality characteristics of plants, as well as the ecological condition of the field. To solve such problems of precision farming, the use of various methods of mathematical statistics is becoming an increasingly promising direction. The aim of our work was to assess the state of agricultural plants using an approach based on the combined use of kriging and binary regression methods, as well as the determination of nitrogen planting using the NDVI (Normalized Difference Vegetation Index) index. The studies were carried out at the site of an experimental agricultural field located on the territory of the branch of the Agrophysical Institute (Menkovo, Leningrad region) in 2015. With the help of aerial photographs taken from the automatized unmanned aerial vehicle complex Geoscan-401 (Geoscan Group of Companies, Russia), a set of NDVI (Normalized Difference Vegetation Index) vegetation index values was obtained at arbitrary points of the plot. A number of ground-based measurements were also conducted on the studied area of the field. The proposed approach to assessing the state of agricultural plants consisted in the joint use of two methods of mathematical statistics: ordinary kriging and logistic regression. A preliminary variogram analysis was carried out, and a variogram model was constructed. After this, the kriging method was used to calculate a series of predicted values of the parameter being studied. At the next stage, the threshold value of the parameter for the study area was established, and also a dummy variable was entered, taking the value 1 if the parameter value exceeded the threshold, and 0 otherwise. Then a logit model was built, in which one of the factors was a series of estimates of the parameter of interest, obtained using the ordinary kriging method. The input data for building logit models were as follows: N(x_i) is the NDVI value at the location x_i, i = 1.78; variable T = 1, if N(x_i) ≥ 0.46, otherwise T = 0; the variables X and Y are the coordinates of the observations, are considered as explanatory variables; N_pred(x_i) is parameter values, predicted using the kriging method at the observed points. All calculations were performed using the R programming language. As a result of the experiment, three logit models were built with the dependent variable T: in the first model, the explanatory variables X and Y; in the second model — X, Y and N_pred; in the third model N_pred. Testing showed that when adding the N_pred variable, the logit model works better (2 times less than the erroneous determination of the level of the parameter under study). The results obtained suggest that adding in the binary regression factors a set of values predicted by the kriging method can significantly improve the accuracy of calculations.

Keywords: plant status, Normalized Difference Vegetation Index, NDVI, kriging, binary regression, language R.