First of all there is a question of terminology

1. "Trend surface" pertains to the data but ultimately it is not the data but rather the expected value of the random function (that presumably generated the data). Unfortunately without knowing the multivariate probability distribution there is no way to actually compute the expected value hence the Trend surface is an ad-hoc solution. Note that any statistical tests in lm() are based on specific statistical assumptions which actually contradict the assumptions implicit in the use of kriging.

2. Geometric anisotropy and zonal anisotropy are characteristics of the random function and its variogram (the model, not the sample variogram)

Having said all that there is the practical matter of what to do. First I would suggest making a plot of the data locations with each location coded by the data value as a quick way of checking to see if there is a "trend". The "trend" may have a directional effect and this could lead to thinking there is an anisotropy. Next look at a histogram of the data. Third, ask your self what you know about the phenomenon that generated the data, e.g. is it reasonable to expect a non-constant expected value and/or a directiional effect.

lm( ) should be able to compute the residuals for both a first order polynomial trend and a second order trend. Generate the coded plots for both sets of residuals, how do they compare with each other and with the coded plot of the original data. Compute the histograms for both sets of residuals, again compare with each other and with histogram of original data.

Compute and plot the sample variograms for the original data, for the residuals from first order polynomial and for second order polynomial.

If you can fit reasonable variogram models to any or all of these data sets (original, first order residuals, second order residuals), use the model with cross-validation and the associated data sets.

For each of the three data sets the cross validation should produce four values at each data location (the actual data or residual, the kriged value, the "error" and the kriging std deviation. Since the kriging estimator is unbiased, the "mean error" should be zero so compare with average error in each of the three cases, the mean square normalized error (normalized by the kriging std deviation) should be one so compare with empirical value, Look at the fraction of nomalized errors exceeding 2.5, theoreticalkly this should be less than 0.5. Make a bi-variate plot for  "data value  vs kriged value" the correlation should be high, make a bi-variate plot of kriged values vs error, these should be uncorrelated.

Did you see evidence of anisotropy in the directional sample variograms for the first and/or second order residuals?

If the cross validation results look acceptable for either the first order residuals or the second order residuals then you can proceed two ways

A. Use Ordinary kriging with the residuals data set and add to those the value of the trend surface at each location. You won't get a kriging variance this way (you can't used the estimated variance from lm( ))

B. Use Universal kriging with the original data set, the order of the trend surface but not the coefficients in the trend surface. This way you will get a kriging variance at each interpolation location.

Unless you have an explanation related to the phenomenon to justify using an anisotropy I suggest not doing it, in  particular no zonal. gstat does not include a zonal anisotropy. As for the geometric there is always the problem of an insufficient number of pairs for the plotted points in the directional sample variograms (the total number of pairs if fixed by the size of the data set and hence the pairs are split up for the different directions.

Note that the cross validation results as well as any final kriging results may be sensitive to choices such as (1) whether you use a circular or elliptial seach neighborhood (2) the minimum and maximum number of data locations used in the search neighborhood as well as the variogram model type (exponential, gaussian, spherical, etc) but not as sensitive to the variogram parameter values

I am not active in the image classification

Hi Donald Myers, As usual, you have given not only an outstanding explanation but also a complete strategy for a research article. Thank you very much.