Unidad de estudio (UE): Ubicación y características ambientales

La UE comprendió dos parcelas (P1 y P2) de siete hectáreas, cada una, ubicadas en el Rancho “18 de Julio”, de la Universidad Autónoma Chapingo (25.9037 N y 103.5821 O); este Rancho fue elegido debido a que está dedicado a la cría y explotación de bovino lechero con una estrategia agro-productiva típica de la CL basada en el cultivo de forrajes para la alimentación del ganado (Figura 1).

Figura 1: Ubicación de la UE. a) Comarca Lagunera y colindancias interestatales con referencia a México; b) Integración municipal sobre imagen aérea o mapa base; y c) parcelas monitoreadas. El círculo en rojo sobre el mapa base (b) indica la ubicación de la parcela 1 (P1) y la parcela 2 (P2). 

Figure 1: Location of the study unit. (a) Comarca Lagunera and interstate borders with reference to Mexico; (b) Municipal integration on aerial image or base map; and (c) monitored plots. The red circle on the base map (b) indicates the location of plot 1 (P1) and plot 2 (P2). 

En ambas parcelas (P1 y P2) de siete hectáreas cada una, en las que se sembró MF, Híbrido-NB777 (^{Novasem, 2022}) con una densidad de 87 944 plantas ha^-1; las fechas de siembra fueron 12 y 22 de abril de 2019 para P1 y P2, respectivamente. El clima de la UE está clasificado como muy árido, semicálido (BWhw) (^{CONABIO, 2019}). Para el periodo 1951-2010, los promedios anuales de temperatura y precipitación pluvial fueron de 21 °C y 248 mm, respectivamente; la evaporación fue 10.4 veces mayor que el promedio anual de la lluvia (2591mm/248mm) (^{CONAGUA, 2021}).

En ambas parcelas se realizó una caracterización edáfica conforme a la NOM-021-RECNAT-2000 (^{SEMARNAT, 2002}) a partir de muestreos realizados a 30 cm de profundidad, cuyos resultados se presentan en la Cuadro 1.

Obtención y procesamiento de insumos

La obtención y procesamiento de insumos requeridos para llevar a cabo la presente investigación comprendió de cinco componentes clave: 1) Localización, descarga y procesamiento de imágenes de satélite; 2) Asociación de los valores de reflectancia y los datos de campo conforme a la resolución espacial de las imágenes de satélite empleadas; 3) Estimación de IV con base en ocho ecuaciones; 4) Gestión de información de campo conforme al tamaño de muestra de la unidad de estudio; y 5) Análisis estadístico (Figura 2).

Figura 2: Diagrama de flujo que ilustra la estrategia de análisis para la modelación espacial de biomasa e índices de vegetación. (R = rojo, IRC = infrarrojo cercano). 

Figure 2: Flow chart illustrating the analysis strategy for the spatial modelling of biomass and vegetation indices. (R = red, IRC = near infrared). 

Fuente, descarga y caracterización de imágenes de satélite

Las imágenes de satélite se obtuvieron de las plataformas del Servicio Geológico de los Estados Unidos de América (^{USGS, 2022}) y de la Agencia Espacial Europea (^{ESA, 2022}); las imágenes LANDSAT-8 y SENTINEL-2 fueron elegidos como insumos para la presente investigación realizada en 2019. Ello se hizo con base en dos razones que coinciden con lo descrito por ^{Zhang et al. (2018)}: 1) por su semejanza instrumental con OLI en el primero y MSI para el segundo; y 2) por su equivalencia en la segmentación electromagnética en el espectro visible y bandas infrarrojas, conformados por 11 y 12 bandas, respectivamente. Además, su amplia disponibilidad durante el tiempo que duró la investigación fue considerada.

Las imágenes de satélite se descargaron para cinco fechas en el caso de LANDSAT-8, y 12 fechas en el caso de SENTINEL-2; la alta resolución temporal (5 días) de SENTINEL-2 permitió estudiar ambas parcelas de manera simultánea mientras que LANDSAT-8 debido su resolución temporal de 16 días permitió evaluar ambas parcelas de modo independiente; las características en cuanto a resolución radiométrica (RR), espacial (RE) y fecha de adquisición (FA) se muestran en la Cuadro 2.

Asociación entre los valores de reflectancia y los datos de campo

La fenología y rendimiento del cultivo fueron asociados con valores de reflectancia extraída de las imágenes de satélite, cuyo procedimiento se desarrolló en dos momentos (Figura 2):

1) Preproceso. El preproceso consistió en lograr un ajuste exacto entre el tamaño de los pixeles o celdas de la imagen con la superficie de las dos parcelas en estudio, para esto se realizó un remuestreo a 15 metros de la imagen LANDSAT y a 10 metros de la imagen SENTINEL-2; la convergencia lograda implicó aproximadamente 40 celdas de 15×15 metros de la imagen LANDSAT-8 y 90 celdas de 10×10 metros de la imagen SENTINEL-2.

2) Postproceso. El postproceso se realizó como segunda fase en el módulo de Clasificación Semiautomática (CS) en Qgis, v. 3.10, cuya técnica es descrita por ^{Congedo (2020)}; este módulo permite la conversión de los valores de reflectancia atmosférica (TOA, por sus siglas en inglés) con presencia de perturbaciones por ruidos o partículas a valores de reflectancia de superficie mediante la aplicación de la corrección atmosférica simple con el método Dark Object Subtraction (^{Chavez, 1996}).

Como resultado del proceso de CS se obtuvieron imágenes calibradas y remuestreadas con resoluciones espaciales mayores a las originales, siendo de 15 m para LANDSAT-8 y de 10 m para SENTINEL.

Estimación de IV

Los IV fueron calculados mediante ocho ecuaciones; seis se identificaron en el grupo de los basados en la pendiente y dos en el grupo de los basados en la distancia (Cuadro 3) (^{Chen, 1996}; ^{Eastman, 2012}). Los datos derivados de su aplicación en cada ecuación prevista permitieron evaluar la calidad del desarrollo y estado de salud de las parcelas de MF conforme a distintas referencias consultadas.

Gestión de información de campo. tamaño de muestra y mediciones realizadas

Convergencia UE vs imágenes de satélite mediante marcas física en tierra. La convergencia entre las imágenes de satélite con los límites y colindancias de la UE fue fundamental. A partir de referencias geográficas obtenidas de las imágenes de satélite se procedió a su identificación en ambas parcelas mediante el uso de un receptor satelital (GPS, Garmin), con lo cual se crearon polígonos delimitados con marcas físicas en cada vértice (Figura 3).

Figura 3: Delimitación de las unidades de muestreo mediante marcas físicas en tierra (a, b y c). Levantamiento de variables morfológicas de las parcelas (d y e). Etapa de floración (f). 

Figure 3: Delimitation of the sampling units by physical marks on the surface (a), (b), and (c). Survey of morphological variables of the plots (d) and (e). Flowering stage (f). 

Obtención de datos de biomasa. Con base en las características de resolución espacial tanto de LANDSAT-8 (15 m × 15 m), como de SENTINEL-2 (10 m × 10 m) y considerando la superficie parcelaria en estudio (70 000 m² c/u) se seleccionaron de modo aleatorio tres plantas de maíz en 10 cuadrantes, 5 cuadrantes por parcela, correspondiendo a un tamaño de muestra de 6.42% para cada parcela.

Análisis Destructivo por Etapa Fenológica. De cada cuadrante se extrajeron tres plantas que fueron trasladadas al laboratorio, donde se determinó peso fresco; los resultados obtenidos fueron extrapolados conforme a la densidad de plantas por hectárea de acuerdo con la siguiente ecuación:

Donde Pbio es el potencial de producción de biomasa verde en kg. Ppp es la biomasa promedio por planta en kg; y NPH es el número de plantas establecidas por hectárea o densidad de plantas por hectárea.

Análisis estadístico

El análisis estadístico se realizó mediante el programa R-Studio, v. 4.0.2 (^{R Core Team, 2022}), para ello se integró por dos apartados: Pruebas de Normalidad (PN) y Análisis de Regresión Lineal Múltiple (ARLM), cuyas características generales son las siguientes:

1) La PN consistió en realizar una prueba de bondad de ajuste de los datos bajo la hipótesis de normalidad, para lo cual se empleó la técnica de Kolmogorov-Smirnov (^{Kim y Jennrich, 1973}).

2) El ARLM permitió obtener un modelo empírico para predecir la producción final de biomasa verde en función de las diferentes bandas de LANDSAT-8, SENTINEL-2, así como los índices de vegetación calculados para las diferentes etapas de desarrollo del cultivo en ambas parcelas.

Los criterios para juzgar la bondad de ajuste del modelo empírico se basaron en el análisis de los errores y se utilizaron para ello los estadísticos: coeficientes de determinación ajustado (R² ajustado), factor de inflación de la varianza como indicador de multicolinealidad (FIV) (^{Lopes, Ferreira, da Silva y Cespedes, 2009}) y el estadístico de Durbin-Watson (D-W) para descartar efectos de autocorrelación. Los estadísticos para la elección del mejor modelo, análisis de varianza y prueba de falta de ajuste se presentan en el Apéndice 1, cuya interpretación se realizó con base en ^{Savin y White (1977)}.

Con el índice de vegetación que mostro mejor desempeño en la validación cruzada y la biomasa estimada, una línea limite que representa la relación o el efecto límite del índice de vegetación con la variación espacial del rendimiento fue estimada (^{López-García, Mata-González, Blanco-Macías, Méndez-Gallegos y Valdez-Cepeda, 2016}).

Relación entre las etapas de crecimiento del cultivo y su respuesta espectral

Las firmas espectrales de SENTINEL-2 para las etapas tempranas del cultivo a los 11 y 27 dds presentaron mayor reflectancia de la energía incidente en el espectro visible y mayor absorción en las bandas 8 y 8a que corresponden al infrarrojo cercano (1 y 2). Esta respuesta espectral del dosel del cultivo se redujo a los 57 dds; es decir, la mayor absorción se presentó en la banda 4 que corresponde al rojo del espectro visible y la mayor reflectancia se asoció las bandas infrarrojas (Figuras 4c, 4d).

LANDSAT-8 presentó un comportamiento similar al observado en SENTINEL-2; en la etapa inicial de crecimiento del cultivo, la absorción de energía en el espectro visible disminuyó, mientras que aumentó en el infrarrojo cercano, pero la reflectancia en el infrarrojo aumentó a medida que el cultivo creció, mientras que en el espectro visible disminuyó (Figuras 4a, 4b).

Figura 4: Respuesta espectral en diferentes días después de la siembra del maíz forrajero registradas por los sensores de LANDSAT-8 (a=P1 y b=P2) y SENTINEL-2 (c=P1 y d=P2). 

Figure 4: Spectral response on different days after forage corn sowing recorded by LANDSAT-8 (a=P1 and b=P2) and SENTINEL-2 (c=P1 and d=P2) sensors. 

Modelación LANDSAT-8 vs SENTINEL-2

El mejor modelo para estimar rendimiento de biomasa verde del MF en ambas parcelas se obtuvo del conjunto de bandas SENTINEL-2; para su definición se emplearon las firmas espectrales derivadas de las bandas 3 y 8 correspondientes a los 72 dds como variables explicativas, teniendo estimadores βi muy significativos (P < 0.05); el modelo elegido fue el que tuvo el coeficiente de determinación mayor (R² = 0.66), cuya expresión se muestra enseguida:

Donde α₃ y α₈ son las bandas 3 y 8 de SENTINEL-2 y e representa al error.

Con base en el modelo anterior (ecuación 10), los valores espectrales al sustituirse con las bandas 3 y 8 (de SENTINEL-2) se produjo un mapa de distribución espacial de biomasa para las 14 ha que suman las dos parcelas estudiadas como se muestra en la Figura 5. La producción total cosechada en P1 fue de 217.39 toneladas, mientras que para P2 fue de 120.54 toneladas. La biomasa cosechada fue de 17.22 y 31.05 Mg ha^-1 para P2 y P1, respectivamente, lo cual tuvo un fuerte contraste con la media regional de 50 toneladas por hectárea.

Figura 5: Biomasa estimada en las parcelas 1 y 2 a los 72 días después de la siembra con los valores de reflectancia de las bandas 3 y 8 de SENTINEL-2. 

Figure 5: Estimated biomass in plots 1 and 2 at 72 days after sowing with reflectance values of SENTINEL-2 bands 3 and 8. 

IV evaluados con mejor desempeño y su representación geoespacial

La relación espacial entre la producción de biomasa verde y cada uno de los IV fue analizada a los 72 dds al igual que el modelo de regresión para la estimación de la cosecha. La elección del IV con mejor desempeño se llevó a cabo mediante una evaluación cruzada (EC) entre los IV y la biomasa estimada utilizando la ^{ecuación 10} para ambas parcelas. Los IV que representaron de manera adecuada la distribución espacial de la biomasa, para ambas parcelas, fueron los del RDVI, los cuales explican 75.41% de la variabilidad en la distribución espacial de la biomasa, seguido por los de SAVI con 13.7% y MSAVI2 con 9.23% (Figura 6). Los otros índices mostraron niveles pequeños de relación con la variación de biomasa; los menores corresponden a RVI con 0.05%, RATIO con 1.66%, NDVI con 1.35% y TVI y TTVI con 1.13 por ciento.

Figura 6: Comportamiento espacial de la biomasa e índices de vegetación analizados en el estudio: (a) Biomasa, (b) RDVI, (c) SAVI, y (d) MSAVI2. 

Figure 6: Spatial behavior of biomass and vegetation indices analyzed in the study: (a) Biomass, (b) RDVI, (c) SAVI, and (d) MSAVI2. 

La Figura 6 indica que a medida que incrementa el valor de los IV, la producción de biomasa aumenta. La Figura 7 proporciona un análisis del índice RDVI que fue el que mostró mayor relación espacial con la biomasa a considerar la técnica de curva límite o envolvente.

Figura 7: Comportamiento espectral de la biomasa y RDVI a los 72 días después de la siembra del maíz forrajero.  

Figure 7: Spectral behavior of biomass and RDVI at 72 days after sowing forage corn. 

En la función cuadrática de la curva límite, lo apreciable es que la producción mínima de biomasa verde fue de 5.893 Mg ha^-1; ello indica que cuando el valor de RDVI asume valores cercanos a 0.210 la producción de biomasa es pequeña. Cuando el valor de RDVI asume el valor 0.469, la producción de biomasa se maximiza hasta 63.017 Mg ha^-1 (Figura 7).

Propiedades radiativas del MF y tiempos de predicción de rendimiento

Los resultados indicaron (Figura 4) que en edades tempranas la absortancia o capacidad de absorber la luz por el dosel del cultivo fue menor en el espectro visible (0.45-0.68 µm), de modo que a medida que creció el cultivo hubo una absorbancia estimada de 8 y 16% en las bandas azul-verde y roja, respectivamente. Este comportamiento se invirtió en el infrarrojo cercano (0.78-90 µm) en la medida en que creció el cultivo; es decir, la reflectancia aumentó y la absortancia disminuyó a una tasa aproximada de 14%, lo cual es consistente con lo descrito por ^{Hamlyn (2014)}.

Otros estudios (p. ej.: ^{Abd-El, Smith, Rowland, y Abd, 2019}; ^{Elders et al., 2022}; ^{Hamlyn, 2014}; ^{Wang et al., 2021}) han consignado que la variación de estas dos propiedades radiativas (absortancia y reflectancia) está fuertemente influenciada, además de la edad del cultivo, también por tamaño, grosor y orientación de la hoja, contenido de humedad y morfología; todos esos atributos varían a medida que crece el cultivo.

A los 72 dds, en la etapa de floración, las bandas 3 y 8 de SENTINEL-2 explican con mayor precisión la variabilidad de la producción final de biomasa con RLM en ambas parcelas (R²=0.66); esto concuerda con ^{Wang et al., (2021)}, quienes señalan que, al alcanzar el máximo crecimiento del dosel, en la etapa de floración, la respuesta espectral de sensores remotos tiene mayor correlación con el rendimiento. También, ^{Dong et al. (2020)} precisan que la radiación en el espectro visible está muy correlacionada con el rendimiento, justo en el momento que el cultivo alcanza su crecimiento máximo.

Índices de vegetación y modelación de rendimiento

Las diferencias entre el coeficiente de determinación obtenido por el modelo de regresión multivariada realizado en este estudio (R²=0.66) y los consignados por ^{Prasad, Chai, Singh y Kafatos (2006)} (R²=0.78) y ^{Xin et al., (2013)} (R²=0.77) podrían deberse al menos a tres cuestiones:

Primero. El tipo de sensor empleado por ^{Prasad et al. (2006}) y ^{Xin et al. (2013)} fue el radiómetro avanzado de muy alta resolución (AVHRR, por sus siglas en inglés) y el espectrorradiómetro de imágenes de resolución moderada (MODIS, por sus siglas en inglés), respectivamente.

Segundo. En ambos estudios se usaron datos estadísticos de rendimiento en cultivos de maíz y soya.

Tercero. Las áreas referidas en estos dos estudios citados (^{Prasad et al., 2006}; ^{Xin et al., 2013}) fueron desarrollados en regiones grandes con avances tecnológicos mayores.

Lo anterior concuerda con lo expresado por ^{Ali et al. (2022)}, quienes señalan que las características edáficas y las labores agrícolas como son: barbecho, nivelación, riego, fertilización o disponibilidad de nutrientes y malezas, entre otras, tienen impacto directo sobre el rendimiento.

Por otro lado, los mejores IV para la modelación de la producción de biomasa a nivel parcelario del MF fueron RDVI, SAVI y MSAVI2. El RDVI es un índice basado en la pendiente, mientras que SAVI y MSAVI2 son índices basados en la distancia; estos últimos proporcionan información libre de ruidos e imprecisiones ocasionados por la presencia de suelo desnudo (^{Eastman, 2012}). Al aplicar la técnica de validación cruzada, el índice RDVI mostró una mejor modelación espacial de la producción de biomasa al presentar un R² = 0.75 en la relación lineal entre el índice y la biomasa estimada (Figura 6).

^{Ahmad, Singh, Fahad y Waqas (2020)} encontraron, a nivel distrito en Pakistán, un R² = 0.78 con datos LANDSAT-8; ^{Madugundu, Al-Gaadi, Tola, Kayad y Jha (2017)} estimaron, en Arabia Saudita, un R² = 0.96 y R² = 0.93 con IV EVI y NDVI, respectivamente, derivados de LANDSAT-8. En contraste, los datos de reflectancia e IV derivados del sensor LANDSAT-8 no modelaron de manera satisfactoria la distribución espacial de biomasa en las etapas de desarrollo del cultivo del presente trabajo; las discrepancias, encontradas principalmente en el NDVI, pueden deberse a factores como las fechas de siembra, mientras que en nuestro caso fue para el ciclo de primavera en el citado estudio la siembra fue en el mes de agosto; por otro lado, el sistema de riego de pivote central contra sistema de inundación empleado en nuestro estudio y sobre todo, a que la biomasa o productividad primaria neta fue calculada por modelos empíricos existentes a diferencia de nuestro estudio donde se realizó un muestreo destructivo al momento de cosecha, entre otros factores de contraste como escalas y sistemas de cultivo.

La función cuadrática como curva límite muestra la estrecha relación entre el RDVI y los rendimientos en P1 y P2. Ello sugiere la existencia de una heterogeneidad grande en el rendimiento de ambas parcelas; ya que a un valor promedio de RDVI= 0.46 puede corresponder un rendimiento mínimo de 20 Mg ha^-1 o bien un rendimiento máximo de 63 Mg ha^-1 (Figura 7). Esta diferencia tan grande entre los rendimientos mínimos y máximos para un mismo valor de RDVI indica que los insumos agrícolas y el riego no se distribuyeron de manera homogénea en ambas parcelas.

Factores que influyeron en el bajo rendimiento en la UE

Los datos de laboratorio (Cuadro 1) permitieron identificar al suelo de la UE como Solonetz cálcico hiposálico de textura franca (SNccszw/2), cuyo predominio y distribución espacial fue confirmada con apoyo de cartografía digital (^{INEGI, 2007}); el origen de este suelo se deriva primordialmente de procesos aluviales, coluviales y eólicos (^{FAO, 2015}).

Los rendimientos pequeños observados en P1 y P2 con respecto a un rendimiento promedio regional de 50 Mg ha^-1 se debieron al desnivel en ambas parcelas; la nivelación ineficiente inhibió la distribución homogénea del riego parcelario, no obstante que se aplicó una lámina de 60 cm (el óptimo 500-900 mm), que incluyó el riego de pre-siembra y tres de auxilio; ello pudo traer consigo, como lo señalan ^{Espinosa et al. (2018)}, dificultades de disponibilidad de nutrientes para las plantas.

El estado del tiempo atmosférico a lo largo del ciclo es un factor que influye directamente sobre el crecimiento y el rendimiento del cultivo; la temperatura óptima para la germinación de maíz es de 16-32 °C; por ejemplo, ^{Granados-Niño et al. (2022)} obtuvieron 55.5 Mg ha^-1 de MF cultivado en la CL en el ciclo PV, cuyas temperaturas máxima, media y mínima fueron de 34.42, 27.05 y 19.61 °C respectivamente, y precipitación acumulada de 24.29 milimetros.

En contraste con el citado estudio, las temperaturas promedio del ciclo (abril-junio) en P1 y P2 fueron en el orden de 36.75, 26 y 15.25 °C, y precipitación acumulada de 10.90 mm; estas condiciones probablemente ocasionaron que los rendimientos de MF fueran 24 Mg ha^-1 menos en P1 y 38 Mg ha^-1 menos en P2, con respecto a lo consignado en el estudio citado (i.e. ^{Granados-Niño et al., 2022}).