Mapping soil salinity in the Yangtze delta REML and universal kriging E-BLUP revised

Mapping soil salinity in the Yangtze delta:REML and universal kriging (E-BLUP)revisitedH.Y.Li a ,b ,c ,R.Webster b ,Z.Shi c ,⁎a School of Tourism and Urban Management,Jiangxi University of Finance and Economics,Nanchang 330013,Chinab Rothamsted Research,Harpenden AL52JQ,Great Britain,UKcCollege of Environmental and Resource Sciences,Zhejiang University,Hangzhou 310058,Chinaa b s t r a c ta r t i c l e i n f o Article history:Received 17February 2014Received in revised form 19May 2014Accepted 21August 2014Available online xxxx Keywords:Saline soilElectrical conductivity Trend surface analysis REMLUniversal kriging E-BLUPRice farmers in China need accurate maps of soil salinity to make rational decisions for management.Modern sensors such as the Geonics EM38conductivity meter,which records the apparent electrical conductivity,EC a ,allied to geostatistics to convert sparse punctual measurements into digital maps can provide them.We have explored the combination in reclaimed land in the Hangzhou Gulf of the Yangtze delta in Zhejiang Prov-ince in south-east China.The EC a was measured at 525points in a 2.2-ha field that was reclaimed in 1996.The data,transformed to logarithms,were treated as the realization of a mixture of strong quadratic trend and corre-lated random residuals.We estimated the coef ficients of the trend and the parameters of the covariance of the residuals by residual maximum likelihood (REML).We then kriged the log 10EC a on to a fine grid by universal kriging (UK),and transformed the predictions back to EC a for mapping.For comparison we also include regres-sion kriging using the estimated variogram of the ordinary least squares residuals from the trend.We compared the results by cross-validation and calculated the mean errors (MEs),mean squared errors (MSEs)and mean squared deviation ratios (MSDRs).All combinations of technique gave small MEs,as expected —kriging is unbiased.The MSEs varied somewhat.The MSDRs,which ideally should equal 1,varied more.The combination with an MSDR closest to 1was UK with the spherical variogram estimated by REML;its MSDR was 0.993.We matched the predictions to the US Department of Agriculture ’s classes of soil salinity and found that approx-imately half of the field fell into its slightly saline and moderately saline classes,where rice yields would yield a pro fit to the farmer,and half were in the very saline and extremely saline classes where rice yields were so poor that the farmer would lose money by attempting to grow rice.©2014Elsevier B.V.All rights reserved.1.IntroductionDuring the past 40years,more than 400000ha of the tidelands in the Yangtze delta of China has been reclaimed for agriculture (Huang et al.,2008).When first enclosed the soil is saline,and most of the salt must be leached away before the land can be used to grow rice.Barley,wheat and cotton are more tolerant of salt,but they are less pro fitable.Farmers have relied on their experience to decide when to plant rice,but often their crops have failed because the soil was still too nd managers therefore want accurate estimates of the soil ’s salinity in map form so that they can judge when and where they can success-fully grow rice (Li et al.,2013).Fortunately,modern proximal soil sensors based on electromagnetic induction (EM)such as the EM31and EM38(McNeill,1980)enable man-agers or their advisors to measure the soil ’s apparent electrical conductiv-ity (EC a )fairly quickly from above the ground surface.The technique is now widely used in surveys of soil salinity (Akramkhanov et al.,2014;Guo et al.,2013;Trianta filis et al.,2013).Even so,the data from the instru-ments are effectively at points in the field with more or less large dis-tances between them.Interpolation is usually necessary for mapping,and nowadays kriging is used for the purpose.For example,for our earlier paper (Li et al.,2013)we had recorded the EC a at only 56points from which to map a field of 2.22ha.In a more recent survey of the field we made 525observations,and these revealed a strong trend;we describe the detail below.De Clercq et al.(2009)and Akramkhanov et al.(2014)also encountered strong trends when mapping soil salinity.Ordinary kriging,the familiar ‘workhorse ’of geostatistics and now readily accessible in computer packages,is based on the assumption that the variable of interest is intrinsically stationary.If there is trend,however,the assumption is untenable,and a more elaborate model of the variation is needed to take into account the trend.The trend might also be of interest in its own right and not simply be a nuisance.Matheron (1969)introduced what he called ‘universal kriging ’to deal with the situation.It is based on the model:Z x ðÞ¼u x ðÞþεx ðÞ:ð1ÞGeoderma 237–238(2015)71–77⁎Corresponding author.E-mail address:shizhou@ (Z.Shi)./10.1016/j.geoderma.2014.08.0080016-7061/©2014Elsevier B.V.All rightsreserved.Contents lists available at ScienceDirectGeodermaj o ur n a l h o m e p a g e :w ww.e l s e v i e r.c o m/l o c a t e /g e o de r m aIn this model u (x )is a non-stationary mean at place x representing the trend,and ε(x )is a stationary random residual with zero mean and covariance function or variogram (h ).The trend is typically treated as a low-order polynomial function of x ,f (x ),linear,quadratic or cubic.Matheron ’s technique requires a model for that variogram,but Matheron did not say how to estimate it.The problem is that until one knows the trend one cannot derive a valid model for the variogram,and without that model one cannot estimate the trend correctly.Various attempts have been made to get round this impasse.Olea (1975)devised an algo-rithm for data on regular grids.Another approach,recommended by Goovaerts (1997)and applied by Meul and Van Meirvenne (2003),is to compute and model the variogram perpendicular to a dominant trend.Neither technique is general,however.Regression kriging in which the interdependence of the trend and variation in the residuals is disregarded enjoys some popularity.Its pre-dictions are unbiased,but the prediction variances underestimate the true prediction errors,often seriously (Lark and Webster,2006);see Lark et al.(2006)and Webster and Oliver (2007)for an explanation.Minasny and McBratney (2007)make the point that its bias decreases as the number and density of data increase,and it might be the only feasible technique for handling the very large numbers of data from proximal sensors ‘on the run ’.Stein (1999)drew geostatisticians ’attention to likelihood methods,though the first statistician to recommend them seems to have been Kitinidis (1983,1993)with what he called ‘generalized covariance func-tions ’,from which any trend has been removed.The latter ’s approach is effectively that of residual (or restricted)maximum likelihood (REML),originally devised by Patterson and Thompson (1971)for estimating components of variance.The method estimates both the trend coef fi-cients and the covariance of the residuals simultaneously.The two are then combined with the data to provide empirical best linear unbiased predictions (E-BLUP)(Lark et al.,2006;Minasny and McBratney,2007).Once a covariance function has been estimated in this way E-BLUP pro-duces the same predictions as universal kriging does with that function (Stein,1999).This method is quite feasible for the few hundreds of data that accrue typically from laboratory measurements and static sensors.It is now regarded as best practice,and it is the one we use here to explore the variation in salinity in the Yangtze delta.2.Site and Methods 2.1.Study Area and SamplingThe land in the coastal zone of Zhejiang Province south of HangzhouGulf of the Yangtze delta is formed of recent marine and fluvial deposits.The soil consists predominantly of uniform pro files of light loam orsandy loam textures,with a sand content of about 60%.It is also saline,with large concentrations of Na and Mg salts (in many places N 1%).The climate is subtropical with an average temperature of 16.5°C and a mean annual rainfall of 1300mm.During the past 30years much of this zone has been enclosed and reclaimed for agriculture.The fields we describe below were reclaimed in 1996and were first used to produce irrigated cotton.Since 2006these fields have been farmed for paddy rice.However,whilst reclama-tion has been fairly successful,salinity is now increasing,and countering it is becoming problematic.Means are required to measure and monitor the dynamics of salinization so that the land can be managed for rice.To test the ability of the EM38to provide the necessary information we chose to study a field approximately 2.2ha.The fields lie to the north of Shangyu City at 30°9′N,120°48′E (Fig.1).We measured the EC a with a Geonics EM38conductivity meter with the coils con figured vertically at 525points,fairly evenly spread in the field.We did so after the rice had been harvested.Each position was georeferenced by a Trimble Global Positioning System.Fig.2is a ‘bubble plot ’of the data with ‘bub-bles ’,i.e.circles or discs,of diameter proportional to the measured values.Table 1summarizes the statistics of the data.2.2.Trend Surface Analysis and Interpolation by Ordinary Kriging The bubble plot (Fig.2)shows a zone of large values towards the south-east of the region.The pattern seems to have the general form of a quadratic trend surface,and we therefore fitted such a surface ini-tially by ordinary least squares (OLS)regression.The model expressed in Eq.(1)thus becomes:Z x ðÞ¼u x ðÞþεx ðÞ¼X K k ¼0b k f k x ðÞþεx ðÞ;ð2Þin which for K =5the set f k (x )is the expansion of a constant,b 0,plus linear and quadratic terms in x .The quadratic surface is shown in Fig.2;it accounts for 71.6%of the variance.The distribution of the measurements is fairly strongly skewed,evi-dent in Fig.3(a),with skewness coef ficient is 0.62.We therefore trans-formed the measurements of EC a to common logarithms for further analysis.The histogram of the logarithms is shown in Fig.3(b)and is symmetric with skewness coef ficient −0.05(Table 1).The model contains an error term,ε,and in regression this is as-sumed to be independently and identically distributed with mean 0and variance σ2.However,because our survey is systematic with no in-dependence by design,this assumption cannot be made.Fig.4in which the variogram of log 10EC a increases approximately linearly without bound suggests that there is a trend with correlatedresiduals.Fig.1.Location of the region studied.72H.Y.Li et al./Geoderma 237–238(2015)71–772.3.Estimation of the Variogram by REML and Interpolation by UK The mathematical details of the separation of the fixed and random ef-fects by REML are now described in the pedological literature (Lark et al.,2006;Minasny and McBratney,2007;Webster and Oliver,2007),and we do not repeat them here.It will suf fice to say that the parameter values of the chosen variogram model,θ,are obtained by maximization of the log-likelihood of the residuals from the trend and the coef ficients of the trend are found simultaneously.For a quadratic trend the coef ficients b 0,b 1,…,b K for the fixed effects are found by generalized least squares.In our study we have found that the popular spherical model with three parameters fit well,so that θ≡{c 0,c 1and r }with c 0the nugget variance,c 1the sill of the structured variance and r the range.We used GenStat (Payne,2013)for the purpose.Appendix A lists the instruction code.The second step in the procedure is the calculation of the predictions and their variances either by UK,or equivalently E-BLUP,the formulae for which are in the literature cited above.Finally in our study the predicted values,Yˆ,were back-transformed to the original scale of EC a ,Z ˆ,byZ ˆx i ðÞ¼exp Y ˆx i ðÞÂln 10þ0:5σ2Y x i ðÞÂln 10n o :ð3ÞIn which,σY 2is the prediction variance.We kriged over blocks 2m ×2m in Matlab,version R2012b,and mapped the results in ArcGIS 10.3.Results and Discussion 3.1.Summary StatisticsTable 1summarizes the statistics of 525EC a (in mS m −1)measured in the topsoil.The standard measure of a soil ’s salinity is that of thesaturation extract,EC e in dS m −1.We can convert our observations to this standard (Li et al.,2013)by EC e ¼0:000671þ0:04897ÂEC a :In Table 1we have added the mean,median,minimum and maxi-mum on this scale.The US Salinity Laboratory (USDA,1954)distinguishes five classes of salinity of increasing severity:Soil with an EC e b 2is regarded as non-saline,2≤EC e b 4is slightly saline,4≤EC e b 8is moderately saline,8≤EC e b 16is very saline,and 16≤EC e is extremely saline.In our sur-vey the numbers of observations in the above five salinity classes are 0,47,243,233and 2,respectively.According to FAO (1976),an EC e of 6dS m −1is likely to cause approximately 25%reduction in yield of rice compared with non-saline soil.We recorded 66.1%of the values larger than this threshold.3.2.Variogram of log 10EC aThe experimental variogram of log 10EC a was first computed by the method of moments:^γh ðÞ¼12m h ðÞX m h ðÞj ¼1z x j −z x j þhn o 2;ð4Þwhere z (x j )and z (x j +h )are observed values transformed to loga-rithms at places x j and x j +h separated by the lag vector h ,and m is the number of paired comparisons at that lag.We treated the variation as isotropic,so that the lag is a scalar in distance only:h =|h |.The points plotted in Fig.4show the result.They follow a steady upward curve in-creasing without bound to which we could fit by weighted least squares a power function,γh ðÞ¼0:000438Âh0:98;ð5Þin which h is in metres.Although the exponent,0.98,is well within the bounds for an intrinsic stationary process.3.3.Variogram of the OLS ResidualsFig.5(a)shows by the plotted points the experimental variogram of the OLS residuals from the quadratic trend,again computed by the method of moments,Eq.(4).These clearly follow a bounded curve,and the spherical model,the solid line,is the model well fitted to them.Its formula isγh ðÞ¼c 0þc3h 2r −12h r3for 0b h b r ¼c 0þc for h ≥r ¼0for h ¼0:ð6ÞTable 2lists the coef ficients of the quadratic trend,and Table 3lists the parameters of the fitted spherical model.3.4.REML VariogramAs above,we used the REML directive in GenStat to estimate simulta-neously the fixed effects,f (x ),of the quadratic trend and theparametersFig.2.Bubble plot showing the positions where measurements were made with ‘bubbles ’of diameter proportional to the observed EC a in mS m −1and the fitted quadratic trend surface representing the fixed effect.Table 1Summary statistics of the soil's apparent electrical conductivity (EC a )and the equivalent values of EC e ,and EC a transformed to common logarithms.MinimumMaximum Mean Median Std dev.Variance Skewness EC a /mS m −158.4357.4166.1154.571.15049.30.62EC e /dS m −1 2.8617.508.137.57log 10EC a1.772.552.182.190.1890.0356−0.0573H.Y.Li et al./Geoderma 237–238(2015)71–77of the random residuals,θ.Given the evidence from the OLS analysis we initially did so for a spherical model.The maximized value of log-likelihood is 1169.7,with 516degrees of freedom.The result appears as Fig.5(b).The plotted points in the graph represent the experimental variograms of the residuals computed by the method of moments.The coef ficients of the trend are listed in Table 2,and in Table 3are the esti-mated parameters of the spherical model of the residuals.We note that the spherical variogram from REML has a smaller nugget than that from OLS,whilst the sill is larger.Also,the range is longer.The variogram of the residuals obtained from REML and that of the OLS residuals are similar at short lag distances but diverge increasingly as the lag distance increases.As Cressie (1993)explains,pages 165–168,this bias in the latter is to be expected.3.5.The Predictions3.5.1.Cross-ValidationTable 4summarizes the comparisons among the various combina-tions of variogram and kriging by cross-validation,on punctual supports of course.The table lists mean errors (ME),mean squared errors (MSE),mean kriging variances (σK 2),σ2K ,the medians of the squared errors,the mean and median squared deviation ratio (SDR).We include in thetable the kurtosis coef ficient of the cross-validation errors.The ME,MSE and MSDR are de fined as follows.ME ¼1N X N i ¼1^Z x i ðÞ−z x ðÞn o ;MSE ¼1N X N i ¼1^Z x i ðÞ−z x i ðÞn o 2;MSDR ¼1X Ni ¼1^Z x i ðÞ−z x i ðÞn o 2σK x i ðÞ:ð7ÞIn these equations z (x i )is the observed value at x i and ^Zx i ðÞis the value predicted there.Kriging is unbiased,and so the ME should ideally be zero.The MEs are all very small compared with the mean of the logarithms,2.18.In using kriging we attempt to minimize the squared errors,so one might judge the combination of model and technique that gives the smallest MSE.As the table shows,the differences among the MSEs are small.The model that most accurately describes the variation should produce squared errors equal to the kriging variances,so the MSDR is ideally 1.Table 4reveals the best in this to be for global prediction by UK,with a MSDR of 0.993.The above results derive from the inclusion of all the data except that at the target point in each prediction.With such small nugget variances,one might think that only the nearest 20points would in fluence the prediction,and jobs would run much faster.With that in mind we ran each prediction with only that number,and we report the results for comparison.See the Appendix A for the Genstat code.A model for which the MSDR =1is generally regarded as rk (2009),however,pointed out that a better diagnostic is the median of the SDR.If the cross-validation errors are normally distributed then the expected value of the median SDR is 0.455.In our study all of the median SDRs were less,even though the mean SDR in the best combina-tion was close to 1.0.The likely explanation is that the distributions,though symmetric,are leptokurtic.Fig.6(b)shows the distribution for global cross-validation with the REML variogram model,for which the kurtosis coef ficient is 1.45(Table 4).At the same time,we also calculat-ed the kurtosis of the cross-validation errors for the raw EC a with a spherical model from REML with a quadratic trend surface.These errors are much more strongly peaked,with a coef ficient of 3.13—see Fig.6(a).This marked difference is another reason for our transforming the data tologarithms.Fig.3.Histograms of (a)observed EC a and (b)log 10EC a.Fig.4.Experimental variogram of the log 10EC a computed by the method of moments,shown as points,and the power model fitted to it:γ(h )=0.000438×h 0.98with h in metres.74H.Y.Li et al./Geoderma 237–238(2015)71–773.5.2.MappingWe used the best combination for mapping with the results in Fig.7(a)and (b).The predictions are for blocks of size 2m ×2m at 2-m intervals on a grid from the logged data and back-transformed the predictions as described above,Eq.(3).Fig.7(a)shows the back-transformed predictions.We calculated and mapped also the prediction variances,which we show without transformation in Fig.7b.Note how the pattern matches the con figuration of sampling points in Fig.2.The kriging variances are smallest close to the sampling points and small generally throughout the region away from the boundaries.The vari-ances become large only near the field boundary,beyond which there are no data.The largest variances within the field are near the northern edge,where there is a cottage.3.6.Implications for FarmingAs a final part of this study we compare yields of rice and the soil ’s EC e using data recorded earlier and summarized in Li et al.(2013).Mea-surements of 1000-grain weights at 56positions in the field were matched as closely as possible to the USDA salinity classes as interpolat-ed in the current survey in ArcGIS.Table 5lists the areas and average 1000-grain weights of the four classes present.It is clear that as the salinity increases the yield of rice decreases,because the salinity in flu-ences not only the 1000-grain weights,but also the length of spike,the number of grain on each spike,and seed setting rate.So,the 1000-grain weights on the most saline soil is about half the average on the least sa-line;the total yield is a much smaller proportion,and in some patches the crop failed completely.We can supplement the yield data summarized in Table 5with farm costs and likely returns in a rough cost –bene fit analysis for each class ac-cording to the predicted yield and the minimum purchase price of rice.The cost of growing the crop,1766.8US$ha −1,is that for the environs of nearby Shanghai (Ye et al.,2006).The results show that a farmer would gain ≈1060US$ha −1on the least saline soil and ≈450US$ha −1on the moderately saline soil.He would lose money by attempting togrow rice on the more saline soil;that would not be pro fitable.Another option for the farmer is to grow cotton,which was done in earlier years (Li et al.,2007).Cotton tolerates salinity much better than rice,but it de-mands more labour,and so the financial return is less.In either case,whatever the farmer decides to do,maps of the kind we have made en-able the farmer to make rational decisions about where and when to plant rice.4.ConclusionFarmers in the Yangtze delta need to know how saline the soil is be-fore deciding to plant rice on reclaimed bining measurements of apparent electrical conductivity,EC a ,from instruments such as the EM38conductivity meter (McNeill,1980)with kriging produces maps on which they can rely.In mapping the salinity in the field from such measurements we en-countered a strong trend,which we had to take into account in addition to correlated random variation.We separated the two sets of effects and estimated their parameters by residual maximum likelihood (REML),now recognized as best practice (Lark et al.,2006;Minasny and McBratney,2007).A spherical function seemed best to describetheFig.5.Variograms of log 10EC a ;(a)of the OLS residuals from a quadratic trend surface and (b)the REML residuals from the quadratic trend surface,both with spherical models shown by the curves.Table 2Fixed effects of quadratic OLS trend surface and REML.TermCoef ficients OLS sphericalREML spherical Constant 1.98 2.17x 1 2.80×10−3 6.45×10−3x 2 5.27×10−3 6.00×10−3x 1x 2−14.2×10−6−17.7×10−6x 210.11×10−6−19.1×10−6x 22−43.3×10−6−39.8×10−6Variance explained (%)80.376.2Table 3Parameters of spherical models of logarithmic transforms of the soil ’s apparent electrical conductivity (EC a ),(a)fitted to the OLS residuals and (b)estimated by REML.The symbols are c 0for the nugget variance,c 1for the sill of the autocorrelated variance,and r for the range.Parameters OLSREML c 00.994×10−30.757×10−3c 16.42×10−312.2×10−3c 0+c 17.41×10−313.0×10−3r /m26.553.7Table 4Cross-validation of spherical models (a)for OLS residuals and (b)universal kriging (UK)predictions a .ModelSearch method MEMSEσ2KSDR KurtosisMeanMedian OLS Local −0.000290.0024690.0029640.8550.248 1.30Global −0.000210.0024700.0029040.8710.284 1.36UKLocal 0.000270.0024280.0026200.9550.283 1.64Global−0.000220.0024180.0025430.9930.3241.45aLocal predictions are with the nearest 20data;global are with all data apart from those at the target points.ME is mean error,MSE is mean squared error,σ2K is the mean kriging variance,SDR is the squared deviation ratio and kurtosis is the kurtosis coef ficient of the cross-validation errors.75H.Y.Li et al./Geoderma 237–238(2015)71–77random effect,and its parameter values were inserted into the equa-tions for prediction by universal kriging (or equivalently E-BLUP).The more easily programmed,and therefore popular,regression kriging performed fairly well in our comparisons,perhaps because we had a large sample (525data)so that the errors in the ordinary least squares regression were much the same as those from the combination of REML and UK.The final map of salinity on the EC a scale spanned four classes of salinity,from slightly saline to extremely saline,as recognized by the U.S.Salinity Laboratory.The fact that approximately 48.2%of the land falls in the very saline and extremely saline classes (EC e ≥8mS m −1≡EC a ≥163mS m −1)is serious for the farmer whowishes to grow rice;he or she can expect poor yields on such soil.The farmer could grow cotton,which tolerates salinity much better,but the greater cost of labour for harvesting and the lack of subsidy from the government mean that it is scarcely pro fitable.AcknowledgementsThis research was supported by a grant from the National High Technology Research and Development Program of China (No 2013AA102301),National Natural Science Foundation of China (No 41101197),Ministry of Education,Humanities and Social Science project (No 10YJC910002),and Natural Science Foundation ofJiangxiFig.6.Histograms of the cross-validation prediction errors by UK with REML spherical variogram and global search:(a)on observed EC a ,mS m −1,and (b)on log 10EC a .On both graphs the curves are of the normal distribution for the calculated means and variances of theerrors.Fig.7.(a)Map of apparent electrical conductivity (EC a ,mS m −1)predicted by universal kriging on 2m ×2m blocks;(b)corresponding prediction variances on the transformed logarith-mic scale.Table 5Mean 1000-grain weights of rice at harvest for four salinity classes.Salinity class Area/ha EC e /dS m −1Mean 1000-grain weight/g Yield a /kg ha −1Gross pro fit b /US$ha −1Slightly saline 0.232–427.956650.141064.39Moderately saline 0.904–823.365209.86451.22Very saline1.028–1618.213448.25−298.76Extremely saline0.03N 1614.12247.17−1661.56a Yield =Mean spikes per m 2×Mean grain per spike ×mean 1000-grain weight ×0.01.bThe cost is about 1766.8US$ha −1in the suburbs of nearby Shanghai (Ye et al.,2006),the average rate between CHY and US$is 6.201in 2013,and the minimum grain purchase price is 0.43US$kg −1.76H.Y.Li et al./Geoderma 237–238(2015)71–77Province(No20114BAB213017).We thank the bodies mentioned above and Rothamsted Research for its hospitality.Appendix AGenStat REML codeThe GenStat instruction code is for a variate z,measured at n posi-tions with coordinates x and y,and modelled as a mixture of afixed qua-dratic trend surface and spatially correlated residuals with variogram γ(h).We choose the Euclidean metric.The spherical variogram is usually expressed as in Eq.(6)in the main text.We set a plausible value for r as45.ReferencesAkramkhanov,A.,Brus,D.J.,Walvoort,D.J.J.,2014.Geostatistical monitoring of soil salinity in Uzbekistan by repeated EMI surveys.Geoderma213,600–607.Cressie,N.A.C.,1993.Statistics for Spatial Data,Revised edition.John Wiley&Sons,New York.De Clercq,W.P.,Van Meirvenne,M.,Fey,M.V.,2009.Prediction of the soil-depth salinity-trend in a vineyard after sustained irrigation with saline water.Agric.Water Manag.96,395–404.FAO,1976.Prognosis of salinity and alkalinity.Soils Bulletin No31FAO,Rome. Goovaerts,P.,1997.Geostatistics for Natural Resources Evaluation.Oxford University Press,New York.Guo,Y.,Shi,Z.,Li,H.Y.,Triantafilis,J.,2013.Application of digital soil mapping methods to identify salinity management classes in coastal lands of central China.Soil Use Manag.29,445–456.Huang,M.X.,Shi,Z.,Gong,J.H.,2008.Potential of multitemporal ERS-2SAR imagery for land use mapping in coastal zone of Shangyu City,China.J.Coast.Res.24,170–176. Kitinidis,P.K.,1983.Statistical estimation of polynomial generalized covariance functions and hydrologic applications.Water Resour.Res.19,909–921.Kitinidis,P.K.,1993.Generalized covariance functions in estimation.Math.Geol.25, 525–540.Lark,R.M.,2009.Kriging a soil variable with a simple non-stationary variance model.J.Agric.Biol.Environ.Stat.14,301–321.Lark,R.M.,Webster,R.,2006.Geostatistical mapping of geomorphic variables in the pres-ence of trend.Earth ndf.31,862–874.Lark,R.M.,Cullis,B.R.,Welham,S.J.,2006.On optimal prediction of soil properties in the presence of spatial trend:the empirical best linear unbiased predictor(E-BLUP) with REML.Eur.J.Soil Sci.57,787–799.Li,Y.,Shi,Z.,Li,F.,Li,H.Y.,2007.Delineation of site-specific management zones using fuzzy clustering analysis in a coastal saline put.Electron.Agric.56, 174–186.Li,H.Y.,Shi,Z.,Webster,R.,Triantafilis,J.,2013.Mapping the three-dimensional variation of soil salinity in a rice-paddy soil.Geoderma195–196,31–41.Matheron,G.,1969.Le krigeage universel.Cahiers du Centre de Morphologie Mathématique, no1Ecole des Mines de Paris,Fontainebleau.McNeill,J.D.,1980.Electromagnetic terrain conductivity measurement at low induction numbers.Technical Note TN-6Geonics Limited,Mississauga,Ontario.Meul,M.,Van Meirvenne,M.,2003.Kriging soil texture under different types of nonstationarity.Geoderma112,217–233.Minasny,B.,McBratney,A.B.,2007.Spatial prediction of soil properties using EBLUP with the Matérn covariance function.Geoderma140,324–336.Olea,R.A.,1975.Optimum mapping techniques using regionalized variable theory.Series on Spatial Analysis,No2Kansas Geological Survey,Lawrence,Kansas.Patterson,H.D.,Thompson,R.,1971.Recovery of inter-block information when block sizes are unequal.Biometrika58,545–554.Payne,R.W.,2013.The Guide to GenStat Release16—Part2:Statistics.VSN International, Hemel Hempstead.Stein,M.L.,1999.Interpolation of Spatial Data:Some Theory for Kriging.Springer-Verlag, New York.Triantafilis,J.,Terhune,C.H.,Monteiro Santos,F.A.,2013.An inversion approach to gener-ate electromagnetic conductivity images from signal data.Environ.Model Softw.43, 88–95.USDA,1954.Diagnosis and improvement of saline and alkali soils.Agriculture Handbook No60United States Department of Agriculture,Washington,DC.Webster,R.,Oliver,M.A.,2007.Geostatistics for Environmental Scientists,Second edition.John Wiley&Sons,Chichester.Ye,L.A.,Wu,Y.X.,Mao,G.F.,2006.Analysis of input–output of paddy production in Shanghai suburb.J.Anhui Agric.Sci.19,5088–5090.77H.Y.Li et al./Geoderma237–238(2015)71–77。