INTRODUÇÃO

A agricultura depende muito da pesquisa experimental para que se consigam informações necessárias ao seu desenvolvimento (^{Gomes e Garcia, 2002}). Segundo os mesmos autores, experimentos de campo são, geralmente, divididos em blocos que são parcelas ou áreas relativamente pequenas, amostradas aleatoriamente. No entanto, quando há dependência espacial dos atributos do solo, bem como a semelhança de comportamento entre eles, sugere que se analise a variabilidade espacial dos atributos do solo em conjunto, antes de adotar um delineamento experimental ao acaso, visto que qualquer tratamento adotado nesta parcela que necessite da homogeneidade do solo levará a falsos resultados (^{Grego e Vieira, 2005}).

Desta forma, uma possível solução é a definição prévia de zonas de manejo (MZs), utilizando o mapeamento do solo, uma vez que este faz um agrupamento dos solos que possuem atributos semelhantes quanto aos fatores limitantes à produção (^{Irigoin et al., 2016}). Assim, os mapas de solos de determinadas áreas possibilitam uma visualização da distribuição espacial e da variação existente na população dos solos (^{Teske et al., 2015}).

Nesse sentido, foi observado que a geoestatística pode ser usada na elaboração de mapas temáticos, a partir de semivariogramas, e na obtenção de valores característicos dos atributos dos solos, utilizando a técnica de interpolação Krigagem (^{Dai et al., 2018}). A variabilidade dos atributos do solo pode ser avaliada separadamente, mas isto dificulta a definição de MZs devido a muitas variáveis que atuam simultaneamente. Assim sendo, é relevante uma análise da variabilidade do solo baseada na análise dos componentes principais (ACP).

Segundo ^{Hongyu et al. (2015}), a ACP é uma técnica estatística de análise multivariada que transforma linearmente um conjunto original de variáveis, inicialmente correlacionadas entre si, num conjunto substancialmente menor de variáveis não correlacionadas que contém a maior parte da informação do conjunto original. Portanto, os atributos do solo analisados com a ACP são agrupados de acordo com suas semelhanças, nas áreas estudadas, que por sua vez são separadas por critérios de dissimilaridade, o que leva à formação de grupos (^{Oliveira et al., 2015}).

Alguns trabalhos têm mostrado a eficiência na definição de MZs envolvendo a ACP e a geoestatística, também chamada de geoestatística multivariada (^{Freitas et al., 2015}; ^{Castilho et al., 2017}). Contudo, são escassos trabalhos com esta temática na região semiárida brasileira, que possui uma variabilidade de solos ainda maior, devido às condições edafoclimáticas da região.

Desse modo, objetivou-se com esse estudo: (1) quantificar a variabilidade espacial dos atributos químicos e físicos do solo e (2) definir MZs para a alocação de áreas experimentais agrícolas usando a geoestatística multivariada.

MATERIAL E MÉTODOS

O estudo foi realizado em uma área de experimentação agrícola de 8,86 hectares, localizada no município de Petrolina, Pernambuco, Brasil, na região do Vale do Submédio São Francisco, na latitude 9°19'16,8" Sul, longitude 40°32'38,8" Oeste e altitude média de 400 m. A área era anteriormente ocupada por vegetação nativa do bioma caatinga, sendo removida a vegetação em 2015 por meio de trator de pneu com lâmina frontal, mantendo-se a lâmina a aproximadamente 5 cm acima do nível do solo. Atualmente existem experimentos com espécies olerículas e silvicuturais instalados em parte da área. O clima da região, segundo a classificação de Köppen e Geiger, é do tipo BSh - semiárido quente, com temperatura média de 32,7°C, umidade relativa média de 56,7% e precipitação média anual de 373,3 mm distribuídos em quatro a cinco meses durante o ano (^{Alvares et al., 2013}).

Para definição da classificação taxonômica do solo da área em estudo, abriu-se uma trincheira e procedeu-se a descrição morfológica do perfil (Quadro 1), seguindo metodologia proposta por ^{Santos et al. (2015}).

Também foram coletadas amostras de solo nos horizontes delimitados para realização das análises químicas e físicas (Quadro 2) segundo metodologia proposta por ^{Teixeira et al. (2017}). O solo foi classificado segundo o Sistema Brasileiro de Classificação de Solos (EMBRAPA, 2018) como Neossolo Quartzarênico órtico típico sendo sua equivalência no American Classification Soil Taxonomy (Soil Survey Staff, 2014) como Quartzipsamments.

Definiram-se, na área, 56 pontos de amostragem conforme verificado na Figura 1. Esse número de pontos amostrais está dentro do intervalo exigido para um mapeamento ultradetalhado segundo recomendações do IBGE (2015). Em cada ponto realizou-se a coleta deformada de solo nas camadas de 0,00 - 0,20 m e 0,40 - 0,60 m de profundidade utilizando-se um trado do tipo holandês. As profundidades de coleta foram baseadas na avaliação morfológica do perfil do solo, sendo a primeira camada (0,00 - 0,20 m) correspondendo ao horizonte A e parte do horizonte AC e a segunda camada (0,40 - 0,60 m) correspondendo ao final do horizonte AC e início do horizonte C.

Realizaram-se as análises de granulometria pelo método da pipeta, pH em água em uma relação solo:água de 1:2, acidez potencial (H+Al), acidez trocável (Al³⁺) e cátions de troca (Ca²⁺, Mg²⁺, K⁺ e Na⁺). Calcularam-se a soma das bases de troca (SB), saturação por bases (V), saturação por alumínio (m) e porcentagem de saturação por sódio (PST). Todos os procedimentos laboratoriais foram realizados conforme descrito por ^{Teixeira et al. (2017}). A determinação da condutividade elétrica (CE) foi feita em meio aquoso conforme a metodologia proposta por ^{Souza et al. (2013}).

Posteriormente, os dados foram submetidos à análise de estatística descritiva (valores médios, máximos, mínimos, desvio padrão e coeficiente de variação). As distribuições dessas variáveis foram testadas quanto à normalidade usando o teste Shapiro-Wilk (p< 0,05).

A partir da matriz de covariância, verificou-se a presença de discrepância nas dispersões relativas (coeficientes de variação) das variáveis. Assim, realizou-se uma padronização dos dados originais das variáveis utilizando sua média e desvio padrão, obtendo-se uma matriz de correlação, por meio da qual se obtiveram os componentes principais (CPs), autovalores, auto vetores, variância explicada e variância acumulada utilizando-se o programa estatístico R Versão 3.5.1 (Development Core Team).

Modelos de semivariogramas foram ajustados aos dados para a descrição da dependência espacial entre as amostras (^{Oliver e Webster, 2014}) tanto das variáveis simples como dos CPs. Para a análise dos componentes principais (ACP) foram utilizados apenas os dois primeiros CPs em cada profundidade, os quais obtiveram autovalores ≥ 1 e variância acumulada de 70% (^{Hair et al., 2005}). As informações geradas pelos semivariogramas foram usadas na interpolação de dados para a produção de mapas pelo interpolador Krigagem ordinária (Oliver e Webster, 2014). Para as variáveis que apresentaram tendência, foi realizada a sua remoção pela subtração de uma superfície ajustada por mínimos quadrados, conforme descrito por ^{Vieira et al. (2010}). O índice de dependência espacial (IDE) foi avaliado conforme a classificação proposta por ^{Seidel e Oliveira (2016}).

O algoritmo de agrupamento fuzzy c-means foi utilizado para classificar a área em zonas menores e mais homogêneas, utilizando os dois primeiros CPs das duas profundidades estudadas. O método fuzzy possibilita a identificação de elementos amostrais que estão na interface dos grupos, isto é, que se assemelham a mais de um cluster (^{Mingoti, 2007}).

Para a aplicação do algoritmo fuzzy c-means aos dados, utilizou-se o Management Zone Analyst 1.0 (MZA). Os ajustes usados no MZA para o presente estudo foram, medida de similaridade euclidiana, expoente fuzziness = 1,3, número máximo de iterações = 300, critério de convergência = 0,0001, número mínimo de zonas = 2 e número máximo de zonas = 6 (^{Fridgen et al., 2004}).

Para determinar o número ótimo de clusters, utilizaram-se os índices de desempenho de imprecisão fuzzy (FPI - fuzziness performace index) e de entropia de classificação normalizada (NCE - normalized classification entropy) propostos por ^{Odeh et al. (1992}), em que o número ideal de zonas é determinado quando os índices FPI e NCE estão no mínimo e, no caso de os mesmos diferirem, opta-se pelo menor valor entre ambos.

Após a delimitação das zonas de manejo (MZs) com o MZA foi realizado o Teste t de Student (p < 0,05) para comparar as médias dos valores de CPs das MZs. As médias foram comparadas de duas a duas, com a finalidade de confirmar a diferença entre as zonas determinadas pelo algoritmo fuzzy.

Para classificar os valores médios dos atributos do solo encontrados em cada MZ foi utilizada a classificação de interpretação de fertilidade proposta por ^{Ribeiro et al. (1999}) conforme apresentada no Quadro 3.

RESULTADOS E DISCUSSÃO

Os resultados referentes à estatística descritiva para os atributos analisados nas camadas de 0,00 - 0,20 m e 0,40 - 0,60 m de profundidade estão apresentados no Quadro 4.

É possível observar que os coeficientes de variação foram baixos (CV% < 12) para areia e pH, médios (12 ≤ CV% ≥ 60) para argila, silte, H+Al, Al³⁺, Ca²⁺ _{0,00 - 0,20 m}, Mg²⁺, SB, T, V, t e m_{0,40 - 0,60 m} e altos (CV% > 60) para CE, Na⁺, K⁺, Ca²⁺ _{0,40 - 0,60 m}, m_{0,00 - 0,20 m}, e PST conforme a classificação proposta por ^{Warrick e Nielsen (1980}). O CV é um importante indicador da variabilidade presente em uma base de dados e quanto mais elevado ele for, maiores são as possibilidades de haver heterogeneidade nessas variáveis (^{Miqueloni et al., 2015}). De acordo com ^{Sana et al. (2014}) esses resultados com altas diferenças de CV entre as variáveis indicam que a tomada de decisão quanto ao manejo do solo com base nos valores médios de seus atributos pode ser ineficiente.

Pode-se observar, ainda no Quadro 4, que alguns atributos não apresentaram distribuição normal, o que segundo ^{Carneiro et al. (2016}) pode ocorrer em virtude da variação do processo de formação do solo, ressaltando assim a importância da adoção de uma maior precisão do manejo da fertilidade do solo, a fim de tratá-lo de maneira pontual, utilizando métodos de manejo específico. Além disso, a presença de normalidade dos dados não é uma condição exigida para a aplicação de análise geoestatística (^{Mion et al., 2012}) e nem dos componentes principais (CPs) (^{Mingoti, 2007}).

Apesar de ^{Li e Heap (2011}) verificarem que a normalidade dos dados é capaz de interferir no comportamento dos métodos de interpolação espacial, ^{Wu et al. (2006}) comprovaram que, entre os mapas construídos a partir de dados com distribuição normal e os elaborados com dados não normais, as diferenças existentes são muito pequenas. Diante disso, optou-se por não normalizar os dados para interpolação neste estudo.

A partir da matriz de correlação (dados não apresentados), foi possível observar que as variáveis estavam altamente correlacionadas. Dessa forma, a análise dos componentes principais (ACP) foi realizada para resumir a variabilidade dos atributos avaliados. Para a avaliação final foram selecionados os dois primeiros CPs de cada profundidade (Quadro 5).

No Quadro 5 pode-se observar que na camada de 0,00 - 0,20 m de profundidade, o CP 1 explicou 50% da variância total e foi dominado pelos atributos Mg²⁺, Ca²⁺, SB, V, t e m. Para a mesma profundidade, o CP 2 explicou um adicional de 19% do total de variância, sendo determinado pelo teor de argila com correlação negativa. Na camada de 0,40 - 0,60 m de profundidade, o CP 1 explicou 51% da variância total havendo predomínio dos atributos Mg²⁺, Ca²⁺, SB, T, t e V e o CP 2 explicou 24% da variância total prevalecendo H+Al e Al³⁺.

Portanto, para a camada de 0,00 - 0,20 m de profundidade, a ACP agregou sete atributos do solo em dois CPs, explicando 69% da variabilidade espacial e para a camada de 0,40 - 0,60 m de profundidade a ACP reuniu oito atributos do solo explicando 75% da variância total. Segundo ^{Oliveira et al. (2015}) deve-se utilizar um número de CPs que expliquem no mínimo 70% da variância total existente na área, a fim de se obter confiança na tomada de decisão quanto ao manejo do solo.

Quanto aos parâmetros dos semivariogramas ajustados (dados não apresentados), foi possível verificar que as variáveis pH, Al³⁺, Mg²⁺ e m na camada de 0,00 - 0,20 m de profundidade e silte, Na⁺, e PST na camada de 0,40 - 0,60 m de profundidade apresentaram efeito pepita puro (EPP) e, consequentemente, não tiveram dependência espacial. O EPP remete a uma variabilidade presente em pequenas distâncias, necessitando por vezes de uma malha amostral mais densa para a sua representação (^{Sana et al., 2014}).

O menor valor de alcance (55,8 m) foi encontrado para areia na camada de 0,00 - 0,20 m de profundidade e o maior valor (218,6 m) foi encontrado para H+Al na camada de 0,40 a 0,60 m de profundidade, o que aponta maior continuidade espacial desta última variável. Além disso, verificou-se que nenhum dos atributos apresentou o mesmo valor de alcance para as duas profundidades. Isto salienta a existência de variabilidade vertical e horizontal entre as profundidades (^{Dalchiavon et al., 2012}).

Os dados dos atributos, areia nas duas profundidades e K⁺, Al³⁺, T e o CP 2 na camada de 0,40 - 0,60 m de profundidade apresentaram tendência, visto que seus semivariogramas, quando ajustados aos dados originais, apresentaram crescimento contínuo, não ocorrendo estabilização do patamar. Desta forma, foi necessária a realização da remoção das tendências dos atributos supracitados, utilizando a equação de superfície plana simples Z(X,Y) = A₀₀ + A₀₁ Y + A₁₀ X, em que X corresponde a latitude e Y a longitude.

Os mapas dos atributos avaliados que apresentaram autovetor ≥ 0,75 (Quadro 5), construídos a partir dos ajustes dos semivariogramas para as profundidades de 0,00 - 0,20 m e de 0,40 - 0,60 m, são apresentados nas Figuras 2 e 3, respectivamente. Os mapas para cada CP são mostrados na Figura 4. Os atributos Mg²⁺, Ca²⁺, SB, V, t e m apresentaram-se mais correlacionados com o CP 1 na camada de 0,00 - 0,20 m de profundidade do que os demais atributos. Por conseguinte, o mapa de Krigagem do CP 1 (Figura 4A) foi visualmente semelhante aos mapas desses atributos, exceto Mg²⁺ e m que apresentaram EPP (Figuras 2A, B, C e D). A argila foi o atributo mais correlacionado com o CP 2 na mesma profundidade.

Do mesmo modo, o mapa de Krigagem do CP 2 (Figura 4B) foi muito semelhante ao inverso do mapa de argila (Figura 2E), uma vez que esse atributo apresentou correlação negativa. O mapa do CP 1 (Figura 4C) na camada de 0,40 - 0,60 m de profundidade assemelhou-se aos mapas dos atributos Mg²⁺, Ca²⁺, SB, T, t e V (Figuras 3A, B, C, D, E e F) e o mapa do CP 2 para a segunda camada avaliada (Figura 4D) foi semelhante aos mapas de H+Al e Al³⁺ (Figuras 3G e H), já que esses foram os atributos com maior correlação. Observa-se que os atributos que mais se correlacionaram com os CPs, em ambas as camadas, foram os relacionados com a acidez do solo. ^{Oliveira et al. (2018}) também verificaram predomínio dos parâmetros relacionados com a acidez do solo em CPs, sendo o CP 1 predominado por Ca²⁺ e Mg²⁺ e o CP 2 por H+Al e Al³⁺.

Para classificar os quatro CPs em Zonas de Manejo (MZs), utilizaram-se os dados interpolados dos pontos amostrais dos CPs, os quais foram submetidos ao agrupamento fuzzy c-means para agrupar os dados por similaridade, sendo os menores valores de FPI e NCE correspondentes a seis clusters, indicando que este seria o número ótimo de zonas (Figura 5).

Os valores de FPI podem variar de 0 a 1, de modo que quando o índice está mais próximo de 0, as classes são distintas com menor compartilhamento de membros e quando está mais próximo de 1, as classes não são distintas e possuem um alto grau de compartilhamento de membros (^{Fridgen et al., 2004}). Ainda segundo os mesmos autores, o número ideal de clusters é determinado quando cada índice está no mínimo, o que representa uma menor adesão de compartilhamento (FPI) ou maior organização (NCE) como resultado do processo de agrupamento.

As zonas 1, 2, 3, 4, 5 e 6 corresponderam respectivamente a 5%, 9%, 14%, 19%, 24% e 29% da área (Figura 5). A delimitação da área em seis classes distintas implica a formação de zonas muito pequenas e fragmentadas, sendo inviável do ponto de vista do manejo do solo. ^{Nascimento et al. (2014}) afirmam que a subdivisão de uma área em mais de três zonas pode inviabilizar a aplicação de práticas específicas de manejo. Complementarmente, ^{Córdoba et al. (2016}) verificaram que as zonas resultantes de uma divisão da área em poucas classes tendem a possuir limites mais coerentes em relação às zonas decorrentes de divisões superiores a cinco classes, que por sua vez, apresentam zonas pequenas e irregulares. Além disso, ^{Carvalho et al. (2016}) confirmam que as melhores zonas são as que tornam o gerenciamento da adubação do solo e de fertilizantes em locais específicos mais práticos do ponto de vista técnico, operacional e econômico.

Dessa forma, com a finalidade de otimizar a demarcação das MZs, adotou-se a delimitação da área do presente estudo em três classes distintas entre si, visto que o resultado do teste t de Student (p < 0,05) foi significativo para o agrupamento dos dois primeiros CPs dos atributos do solo em três zonas (Quadro 6). O mapa das MZs resultante é mostrado na Figura 6. As zonas 1, 2 e 3 corresponderam a 26,06%, 40,31% e 33,63% da área, respectivamente.

Comparando a distribuição espacial das MZs na área (Figura 6) com o comportamento espacial dos atributos do solo avaliados individualmente (Figuras 2 e 3), é possível verificar que na zona 1 houveram os menores valores dos atributos do solo, tendo esta, registrado na camada de 0,00 - 0,20 m de profundidade, um predomínio de níveis baixos de SB e t e teores de Ca²⁺ variando de muito baixos a baixos, enquanto a camada de 0,40 - 0,60 m de profundidade apresentou predominância de níveis médios de Mg²⁺, baixos de SB, t e H+Al e muito baixos de Ca²⁺.

Na zona 3 foram constatados os maiores valores dos atributos avaliados, ocorrendo predominantemente teores de Ca²⁺ variando de baixos a médios, e níveis médios de SB e t na camada de 0,00 - 0,20 m de profundidade, ao passo que na camada de 0,40 - 0,60 m de profundidade, observou-se níveis bons de Mg²⁺, médios de SB e t e teores de Ca²⁺ variando de baixos a médios.

Na zona 2, houve a ocorrência dos maiores valores dos atributos relacionados a acidez do solo, sendo dominada principalmente por níveis médios de H+Al e teores de Al³⁺ variando de médios a altos na camada de 0,40 - 0,60 m de profundidade. Para os demais atributos analisados, a zona 2 abrangeu valores intermédios nas duas profundidades avaliadas.

Em relação a V, observa-se que as zonas 1 e 2 são caracterizadas predominantemente como distróficas (V < 50%) e a zona 3 é basicamente eutrófica (V ≥ 50%). Assim, dentre as três zonas avaliadas, a zona 3 apresentou as melhores condições, em termos de fertilidade do solo, para o desenvolvimento vegetal.

Os mapas de zonas obtidos neste trabalho permitem uma melhor alocação dos experimentos agrícolas na área de estudo, evitando uma possível heterogeneidade dentro dos blocos experimentais estabelecidos.

CONCLUSÕES

1- A geoestatística possibilitou verificar que a maioria dos atributos do solo apresentou dependência espacial, havendo variabilidade na área, não sendo apropriada a aplicação de práticas agrícolas com base em valores médios.

2- As duas primeiras componentes principais da camada de 0,00 - 0,20 m de profundidade explicaram 69% da variação presente na área e os da camada de 0,40 - 0,60 m de profundidade explicaram 75%, sendo dominados pelos elementos relacionados à prática da calagem.

3- O número ideal de zonas específicas de manejo para a área avaliada foi de três, uma vez que o teste de médias indicou heterogeneidade entre elas e esse número de zonas do solo torna possível a aplicação de práticas de manejo localizada.

4- O delineamento de zonas específicas de manejo com base em componentes principais se mostrou uma ferramenta importante no auxílio da alocação de experimentos agrícolas, reduzindo o erro advindo da heterogeneidade que pode ocorrer dentro dos blocos experimentais.